2.3 模型验证与应用
2.3.1 一维全溃坝水流
Stoker于1957年推导出平底、无阻力、矩形断面、棱柱形河道上一维瞬时全溃坝水流的解析解。尽管这种简化条件下的溃坝洪水过程与实际地形上溃坝洪水演进有较大差异，但由于存在准确解，可用于检验数学模型的精度。因此，该算例被广泛用于浅水动力学数值模型验证。
本算例的计算条件：河渠长1000m，宽100m，坝体位于x=500m处，忽略大坝厚度；t=0时大坝瞬时全溃；河道平底、无阻力；上游及左、右岸为固壁边界，下游为自由出流边界。
分别采用以下两套网格：
（1）三角形网格。采用2000个三角网格剖分计算域，计算网格如图2.9（a）所示，共1111个节点和3110条边，单元的平均面积为50.0m2。
（2）四边形网格。采用1000个三角网格剖分计算域，计算网格如图2.9（b）所示，共1111个节点和2110条边，单元的平均面积为100.0m2。
图2.9 一维全溃坝水流算例计算网格示意图
2.3.1.1 下游河床初始水深大于零
初始条件：上游水位为5.0m，下游水位为0.2m，流速为0。
基于上述条件，模拟计算了t=10s、20s、30s、40s、50s、60s时水位、流速分布情况。水位、流速的数值解与理论解对比分别如图2.10和图2.11所示。由结果对比可知，水位、流速的数值解与理论解吻合较好；在激波附近，流速计算结果存在一定的数值振荡，但振幅较小。此外，图2.10反映了大坝溃决后向下游呈“立波”形式推进的洪水波前。
2.3.1.2 下游河床初始水深为零
初始条件：上游水位为5.0m，下游水位为0m，流速为0。
基于上述条件，模拟计算了t=5s、15s、25s、35s时水位、流速分布情况。水位、流速的数值解与理论解对比分别如图2.12和图2.13所示。由结果对比可知，水位的数值解与理论解吻合较好；在水深较大区域，流速数值解与理论解吻合较好；在干湿界面附近，流速数值解与理论解的偏差较大。
图2.10 不同时刻水位数值解与理论解对比
图2.11 不同时刻流速数值解与理论解对比
图2.12 不同时刻水位数值解与理论解对比
图2.13 不同时刻流速数值解与理论解对比
2.3.2 二维非对称型局部溃坝水流
本算例为平底河道上的二维溃坝问题。如图2.14所示，计算域为200m×200m正方形区域；大坝位于y=100～115m处，坝体宽15m；溃口位于x=95～170m处。上游初始水位为10m，下游初始水位为5m；糙率n=0.020。
分别采用两套网格：
（a）三角形网格。采用8424个三角单元剖分计算域，如图2.14（a）所示。
（b）四边形网格。采用4212个三角单元剖分计算域，如图2.14（b）所示。
图2.14 计算域及网格剖分示意图
t=7.2s时的三维水面和平面流场模拟结果如图2.15和图2.16所示。由模拟结果可知，大坝溃决后，向上游传播的洪水波所到之处水位降低，向下游传播的洪水波所到之处水位上涨，水面和流场模拟结果符合物理规律。结果表明，模型的激波捕获能力强，适合模拟溃坝洪水波的传播和区域洪水淹没过程。
2.3.3 具有干湿界面的静止水流算例
本算例为非平底地形上具有干湿界面的二维静水问题，用于检验模型的和谐性。计算域为1m×1m的正方形区域，不考虑河底摩阻项，底高程为
给定初始水位0.1m，初始流速为0。计算域四周均采用固壁边界条件。结合初始水位值和底高程数据可知，计算域内存在干湿边界，而数值模型必须一直维持水位为常数、流速为零的静水状态。
图2.15 三维水面模拟结果
图2.16 平面流场模拟结果
采用三角形、四边形混合网格剖分计算域，共1276个三角形单元和612个四边形单元，如图2.17所示。模拟计算了500s内的水流运动情况。图2.18给出了t=500s的三维水面模拟结果。图2.19给出了沿直线y=0.5m的水位、单宽流量数值解与理论解对比，结果表明计算水位未发生变化、计算流速为0，模型满足和谐性要求。
图2.17 网格剖分示意图
图2.18 三维水面模拟结果
图2.19 水位、单宽流量数值解与理论解对比
2.3.4 过驼峰的溃坝波传播问题
该算例为河底有3个驼峰的溃坝波传播问题，包括了复杂地形、干湿边界、摩阻力等复杂条件下的非恒定水流运动过程，被广泛用于检验模型计算稳定性、复杂地形和动态边界处理能力等。如图2.20所示，计算域为长75m、宽30m的矩形水槽，大坝位于x=16m处，忽略大坝厚度。大坝上游初始水位为1.875m，初始流速为0；下游为干底河床。糙率n=0.018。水槽四周采用固壁边界条件。水槽内包括两个分别位于（30m，6m）和（30m，24m）、高度为1m的驼峰和一个位于（47.5m，15m）、高度为3m的驼峰，底高程为
由于底高程和初始条件均关于直线y=15m对称，因此，任意时刻的计算结果应保持关于直线y=15m对称。为消除网格结构不对称的影响，将计算域剖分为关于直线y=15m对称的三角形网格，共9000个单元、4606个节点和13605条边，单元的平均面积为0.25m2。
图2.20 计算区域及底高程示意图
模拟计算了300s内的水流运动情况。图2.21为t=2s、6s、12s、24s、30s和300s的三维水面计算结果，直观展示了复杂地形上溃坝洪水波的传播过程。
由图2.21可知，水流特征具有较好的对称性，水流运动符合物理规律。在溃坝波的传播过程中，小驼峰曾完全被水流淹没，其涨水和退水过程明显，而高驼峰顶部未被水流淹没。同时，溃坝波冲击驼峰，并由此产生向上游传播的反射波。t=2s时，小驼峰处已有水流通过，并处于涨水状态；t=6s时，小驼峰已完全被水流淹没，而溃坝波前已到达高驼峰处；t=12s时，洪水已由高驼峰两侧流过，并向下游传播，绕流现象明显；t=24s、30s时，洪水已淹没整个平底区域，下游固壁对水流的反射现象显著；t=300s时，由于溃坝波之间、溃坝波与河床、溃坝波与固壁的相互作用，以及摩阻项引起的能量耗散，水流已趋于静止状态。
图2.22给出了t=2s、6s、12s、24s、30s和300s的水深、流速等值线和平面流场模拟结果，更清晰地展示了溃坝洪水波传播过程。计算结果表明，模型有效模拟了复杂的水流运动过程，适合模拟溃坝洪水演进及其淹没过程。
2.3.5 连续弯道水流模拟
水槽由两个弧度为90°的弯道组成，弯道半径为8.53m，两弯道的连接过渡段为4.27m的直水槽，弯道的进出口由长为2.13m的直段过渡，弯道的横断面为矩形断面，断面宽为2.34m。共布置了1～13个断面，依次对应编号为S0、S1、S2、S3、C110、π/16、π/8、3π/16、π/4、5π/16、3π/8、7π/16和π/2。试验水槽平面布置如图2.23所示。
图2.21 不同时刻的三维水面模拟结果
上游流量为0.0985m3/s，平均水深为0.115m，进口平均流速为0.366m/s，水面比降为0.00035，谢才系数取60.5m1/2/s。采用有滑移固壁边界条件，计算网格如图2.24所示，共计2250个四边形网格。
图2.22 不同时刻的水深、流速等值线和流场模拟结果
图2.23 弯道平面示意图（单位：m）
图2.24 网格剖分示意图
水面线计算结果如图2.25所示。可以看出，受弯道影响，在上弯道左岸水位高于右岸水位，在下弯道右岸水位高于左岸水位。
图2.25 水面线计算结果
流速分布结果见图2.26所示，其中，图2.26（a）和图2.26（b）分别为不考虑、考虑弯道环流的影响，直观反映了弯道环流对局部流速分布的影响。
2.3.6 坡面流模拟
选取1955年Iwagaki在光滑铝板水槽内进行的坡面流实验，用于验证水动力数学模型在暴雨山洪模拟方面的适用性。光滑铝板水槽总长24m，由3段长度均为8m、底坡不同的水槽串联组成。各段水槽内的底坡保持不变，水槽连接处不存在高程突变，从上游到下游各段水槽的底坡依次0.020、0.015和0.010。上游、中游、下游段水槽的降雨强度分别为389cm/h、230cm/h、288cm/h。按降雨历时的不同共3组实验，降雨历时分别为10s、20s、30s，不考虑下渗。
图2.26 流速分布模拟结果
采用混合网格剖分计算域，共857个单元，其中，三角形单元588个，四边形单元269个。网格剖分如图2.27所示。单元糙率取0.0098。水槽两侧及上游边界为固壁，下游为自由出流边界。图2.28为3组实验的实测值与计算值对比。由结果对比可知，本模型计算结果与实测结果较为接近，其中，模型合理再现了降雨历时为10s的实验中出口断面存在流量急剧上涨的过程，表明模型可有效模拟具有大梯度解的浅水流动问题。
图2.27 网格剖分示意图
图2.28 出口断面单宽流量结果对比
2.3.7 山区小流域暴雨洪水模拟
研究区域为某山丘区小流域，该流域地处贵州高原向广西丘陵过渡的斜坡地带，地势北高南低，地形起伏大。流域内沟壑纵横，群山高耸，山谷相间，河溪交错的地貌景观十分分明。
采用三角形网格剖分计算域，共93501个单元、45187个节点，河道局部区域进行网格加密。研究区域总面积为198.6km2，最小网格单元面积为235m2，网格边的最短长度为20m。网格地形如图2.29所示，高程为550～1650m。流域出口断面采用自由出流边界。河道布设了15个统计断面（图2.30）。
图2.29 山区小流域地形示意图
图2.30 流量统计断面位置示意图
2.3.7.1 长历时暴雨洪水模拟
为了验证模型的水量守恒性和处理陡峭地形上坡面流的能力，假设均匀恒定降雨，净雨强度为100mm/h。由经验可知，当降雨历时足够长时，流域出口断面流量将保持为恒定状态，流量值为198.6km2×100mm/h=5517m3/s。出口断面（断面13）的流量过程计算结果见图2.31。由计算结果可知，出口断面流量在t=1.4～3.9h之间急剧上涨，在t=6h后基本保持恒定，流量值为5515m3/s，与理论值基本一致。表明模型能有效处理陡峭地形的干湿边界问题，可合理计算具有复杂地形的山区小流域暴雨洪水汇流过程。模拟过程中水量误差保持在10-5m3数量级，严格保证了水量守恒性。
图2.31 出口断面（断面13）的流量过程计算结果
2.3.7.2 短历时暴雨洪水模拟
假设均匀恒定降雨历时1h，净雨强度为100mm/h。模拟了3h的洪水演进过程。各断面计算流量过程见图2.32。
由图2.32可知，模型合理计算了流域内各河道沿程断面的流量过程，得到了各断面的洪水到达时间、峰现时刻、洪峰流量等洪水要素。以断面3为例，其流量过程合理反映了断面1及断面2流量过程的叠加效应。此外，由计算结果可知，断面7至断面10的洪水到达时间相差不大，而断面11至断面13的洪水到达时间存在较为明显的延迟现象。由于断面7至断面10区间河段地势较为陡峭，洪水归槽明显，洪水流速快，因此，洪水到达时间差别不大；而断面10至断面13区间河段地势相对平坦，洪水在河道及两岸传播，洪水流速较慢，因此，洪水到达时间差别较大。此外，断面14与断面15的流量相对大小合理反映了断面相应集水面积的比例。
图2.32（一） 断面流量过程计算结果
图2.32（二） 断面流量过程计算结果
此外，在整个模拟过程中，未出现非物理大流速问题，计算稳定。本算例历时3h的暴雨洪水模拟，计算时间步长取0.2s，共计算耗时1.27h。
2.3.8 圆形水域内水平风生环境流量模拟
风生流是由风对水面的摩擦剪应力作用引起的，在黏滞力作用下使表层水体带动下层水体向前运动。因而风是引起水体混合的重要动力因素。风作用水面时，剧烈地扰动水面，其切应力将形成水面漂流或风生流，在持续作用下由于岸边界和底边界的制约，水体还将形成垂向环流和水平环流。
本算例模拟在恒定风速作用下圆形水域内的水流运动状态。圆形水域半径R0为200m，静止时水深相对于圆心呈递增趋势，水深用式（2.117）表示：
式中：h为水深；rb是该点相对于圆心的距离。
计算域及地形如图2.33所示。
图2.33 计算域及地形示意图
初始时刻水面静止，水面高程为0m。水面上出现正北风，风速为2.53m/s。水面风应力拖曳系数取常数0.0026，即保持风应力为常数0.02N/m2。
图2.34为流场、流迹线和流速等深图计算结果。由模拟结果可知，垂线平均流速方向在深水区（即水域中心区）和所作用的风向相反，而在浅水区流速方向和风向基本一致。在圆心区两侧存在两个水平环流，流场图沿风向成对称形状。计算结果表明，Hydro MPM2D_FLOW合理模拟了水平风生环流现象，可有效计算水面风应力对水流运动的影响。
图2.34 水平环流模拟结果