2维水动力模型中入流边界优化处理方法

 引用本文: 侯精明, 汪煜, 张兆安, 等. 2维水动力模型中入流边界优化处理方法 [J]. 工程科学与技术, 2022, 54(4): 39-46.
HOU Jingming, WANG Yu, ZHANG Zhaoan, et al. Inflow Boundary Optimized Method in Two-dimensional Hydrodynamic Model [J]. Advanced Engineering Sciences, 2022, 54(4): 39-46.
 Citation: HOU Jingming, WANG Yu, ZHANG Zhaoan, et al. Inflow Boundary Optimized Method in Two-dimensional Hydrodynamic Model [J]. Advanced Engineering Sciences, 2022, 54(4): 39-46.

## 2维水动力模型中入流边界优化处理方法

• 收稿日期:  2021-06-05
• 网络出版时间:  2022-04-25 03:43:27
• ###### 作者简介: 侯精明（1982—），男，教授，博士. 研究方向：地表水及其附随过程数值模型. E-mail：jingming.hou@xaut.edu.cn
• 中图分类号: TV131.2

## Inflow Boundary Optimized Method in Two-dimensional Hydrodynamic Model

• 摘要: 在极端暴雨事件频率和强度不断加剧的趋势下，防洪安全面临着极其严峻的考验，研究河道洪水演进、漫溢及淹没过程，有效评估洪水灾害风险就显得尤为重要。为精确模拟河道洪水演进过程，本文采用基于非结构网格有限体积法建立的2维水动力模型，构建了一种基于虚拟单元法的入流边界改进方法。该方法采用基于地形网格数据的纵断面提取法以近似确定深泓线的位置，从而计算出河道平均纵比降；在入流边界各网格上以均匀流的方式入流，计算出不同流量下的对应水位；依据各网格上的水深和入流边长计算权重以合理地分配流量，再结合单宽流量边界条件确定边界外侧虚拟单元上的水力要素，最后耦合至模型通量部分进行计算，实现入流模块的优化。将该方法应用于Toce河物理模型试验的模拟，在入流边界附近计算出的水位、流速和流态都更为合理；河道中各测点水位与试验数据吻合度较高，洪水到达测点时间及各点间洪水运动传播时间与测量值较为接近，纳什效率系数在0.82～0.95，验证了该方法的可行性和模型的准确性；并且在单机上用时86 s完成计算，效率较高，表明模型能高效准确的模拟洪水演进过程。通过对入流模块的改进优化，模型可以更好地应用于洪水演进过程的模拟预警及灾情评估，为洪水灾害管理决策提供技术支撑。

Abstract: With the increasing frequency and intensity of extreme rainstorms, the flood management situation has become extremely severe. It is particularly important to study the process of river flood propagation, overflow and inundation and evaluate the risk of flood disaster effectively. To simulate the river flood propagation accurately, a two-dimensional hydrodynamic model using the finite volume method based on unstructured grids was applied, an improved inflow boundary processing method was proposed based on Ghost-Cell Method. The method of profile extraction based on terrain grid data was used to determine the position of thalweg approximately, and thus to calculate the average gradient of the river. The corresponding water level under different discharge was calculated by defining that the flow on the grid of inflow boundary is the uniform flow. The discharge per unit width was distributed reasonably by calculating the weight according to the water depth and inflow boundary length of each grid, and coupled to flux part of the model for calculation by determining the hydraulic elements on the ghost cell outside the boundary . The method was applied to the simulation of the physical model test of the Toce river, and the flow state near the inflow boundary was more reasonable. The water level at gauging points in the river is in good agreement with the measured data. The arrival time and travel time of flood at gauging points is close to the measured values. The NSE (Nash-Sutcliffe efficiency coefficient) is between 0.82 and 0.95, which verifies the feasibility of the method and the accuracy of the model. And it takes 86 s to complete the calculation on the personal computer, which shows that the model can simulate flood propagation process efficiently and accurately. Through the improvement and optimization of the inflow module, this model can be better applied to the simulation, early warning and disaster assessment of flood propagation process, and provides more powerful technical support for flood disaster management and decision.

• 在全球气候复杂多变和极端恶劣天气频发的大背景下，洪水灾害给人类造成了巨大的经济损失，威胁着人民群众的生命财产安全[1-3]。为做好防洪减灾和预警预报工作，2维水动力学模型作为洪水演进过程模拟的一种重要技术手段[4-5]，可及时准确地为应急决策部门提供水灾害信息，便于高效开展抢险救援工作。

在2维浅水流动模拟中，边界条件处理是极其重要的部分，若边界条件处理的不恰当，就会影响计算结果，有时甚至会阻碍计算的正常进行[6-7]。所以对边界条件进行合理有效的处理是数值模拟过程中必须面对的问题。除了因极端暴雨形成的流域及城市洪水[8]和堤坝瞬间溃决形成的江河洪水[9-10]外，对于大部分洪水演进模拟而言，若在上游入口处没有输入确定的水位流量关系，通常就只是给定一个入流流量过程，通过水动力模型计算来获得水流沿河道向下游传播时的运动状态，但是断面流量在入流边界处如何合理地分配水力要素就成了研究重点。目前，入流边界处理方法已有很多，最简单的方法就是用断面流量除以入口宽度去计算单宽流量，但这样处理在入流边界附近会出现高处过水、水位间断等非物理现象，而且对入流边界的宽度有一定要求。为此，于普兵[11]提出按过水断面面积分配各网格上的流量；宋利祥[12]根据曼宁和谢才公式，以边界上的水深和网格边长计算权重来分配流量。这些方法都与入流边界上的水深有很大关系，有学者认为边界上的水深就等于相邻网格上的平均水深值[13]，但这样处理对于初始为干地形的河道洪水模拟而言无法实现流量分配；还有的学者设置初始条件时在入流边界网格上赋予一定的水深，以弥补这一缺陷。但这样给定的水深是否合理还尚未可知，而且对于不同入口宽度也会带来一定的问题。

针对上述入流边界处出现的问题，本文采用基于非结构网格有限体积法建立的2维水动力模型[14]，对入流模块做了部分改进，依据给定断面流量计算出的水位来分配入流边界各网格上的水力要素，使得入口附近水流流态更为合理，模型能更好地适用于洪水演进过程的模拟计算及应用。

模型以2维浅水方程为控制方程，其矢量形式可表示为：

 $$\frac{\partial {\boldsymbol{q}}}{\partial t}+\frac{\partial {\boldsymbol{f}}}{\partial x}+\frac{\partial {\boldsymbol{g}}}{\partial y}={\boldsymbol{S}}$$ (1)
 \begin{aligned}[b]& \boldsymbol{\text{其中，}}{q}=\left[ \begin{array}{c} h \\ q_{x} \\ q_{y} \end{array} \right], \boldsymbol{f}=\left[ \begin{array}{c} u h \\ u q_{x}+g h^{2} / 2 \\ u q_{y} \end{array} \right], \boldsymbol{g}=\left[ \begin{array}{c} v h \\ v q_{x} \\ v q_{y}+g h^{2} / 2 \end{array} \right]\text{，} \\& {{{\boldsymbol{S}}}} = {\kern 1pt} {\kern 1pt} \left[ \begin{gathered} 0 \\ - gh\partial {{\textit{z}}_{\rm{b}}}/\partial x - {C_{{\rm{f}}} }u\sqrt {{u^2} + {v^2}} \\ - gh\partial {{\textit{z}}_{{\rm{b}}} }/\partial y - {C_{{\rm{f}}} }v\sqrt {{u^2} + {v^2}} \end{gathered} \right] \\[-33pt] \end{aligned} (2)

式（1）～（2）中：h为水深；qxqy分别为xy方向的单宽流量；g为重力加速度；uv分别为xy方向的流速；fg分别为xy方向的通量矢量；S为源项矢量；zb为河床底高程；Cf为床面摩擦系数，Cf = gn2 / h1/3，其中n为曼宁系数。

该水动力模型使用三角形非结构网格剖分计算域，采用基于Godunov格式的有限体积法对2维浅水方程进行数值求解。在网格单元i内，其积分形式可写为：

 $$\int\limits_{\varOmega} \frac{\partial \boldsymbol{q}}{\partial t} \mathrm{d} \varOmega+\int\limits_{\varOmega}\left(\frac{\partial \boldsymbol{f}}{\partial x}+\frac{\partial \boldsymbol{g}}{\partial y}\right) \mathrm{d} \varOmega=\int\limits_{\varOmega} \boldsymbol{S} \mathrm{d} \varOmega$$ (3)

式中，Ω为控制体i的面积。

应用高斯定理，将式（3）中通量项的面积分转化为边界上的线积分：

 $$\int\limits_{\varOmega} \frac{\partial \boldsymbol{q}}{\partial t} \mathrm{d} \varOmega+\oint\limits_{\varGamma} \boldsymbol{F}(\boldsymbol{q}) \cdot \boldsymbol{n} \mathrm{d} \varGamma=\int\limits_{\varOmega} \boldsymbol{S} \mathrm{d} \varOmega$$ (4)

式中，Г为控制体i的边界，n为边界Г所对应外法线方向的单位向量。相应界面通量F(q)·n可以表示为：

 $$\boldsymbol{F}(\boldsymbol{q}) \cdot \boldsymbol{n}=\left[ \begin{array}{c} q_{x} n_{x}+q_{y} n_{y} \\ \left(u q_{x}+g h^{2} / 2\right) n_{x}+v q_{x} n_{y} \\ u q_{y} n_{x}+\left(v q_{y}+g h^{2} / 2\right) n_{y} \end{array} \right]$$ (5)

式中，nxny为边界外法向单位向量nxy方向的分量。

在三角形非结构网格i中，各通量如图1所示，j1j2j3分别为i的相邻网格单元，则单元i内通量项线积分可写为：

图  1  三角形网格上各边通量
Fig.  1  Flux of each side on triangular grids
 $$\oint\limits_{\varGamma} \boldsymbol{F}(\boldsymbol{q}) \cdot \boldsymbol{n} \mathrm{d} \varGamma=\sum_{k=1}^{3} \boldsymbol{F}_{k}(\boldsymbol{q}) \cdot \boldsymbol{n}_{k} l_{k}$$ (6)

式中，k为网格边的编号，lk为第i个网格单元上第k条边的长度。通量项采用HLLC近似黎曼解算器求解，底坡源项和摩阻源项的具体计算方法见文献[15-16]，本文不再赘述。此外，引入GPU加速技术以大幅提高模型计算效率[17]

水动力模型中常用的边界条件有：入流边界、自由出流边界、固壁边界、水位边界等[18-21]。本文采用虚拟单元边界处理方法[22]，如图2所示，对模型中的边界条件均采用HLLC近似黎曼求解器求解界面通量[23]，边界左侧即内部网格上信息为已知，需构造边界右侧即虚拟网格上的信息包括地形高程、水位、流速等。模型中设置所有边界左右两侧网格上地形高程相等，不同边界右侧网格上的水力要素则根据边界性质分别构造。

图  2  边界处理方法示意图
Fig.  2  Diagram of boundary processing method

在固壁边界上，采用无滑移边界条件，法向和切向流速都为0，认为左右两侧水深相等，则右侧网格上信息构造为：

 $$h_{\mathrm{R}}=h_{\mathrm{L}}, \quad u_{\mathrm{R}}^{\perp}=-u_{\mathrm{L}}^{\perp}, \quad u_{\mathrm{R}}^{\|}=-u_{\mathrm{L}}^{\|}$$ (7)

式中，下标R表示的是右侧虚拟网格，下标L表示的是左侧内部网格， $u^{\perp}=un_{x}+vn{ }_{y}$ $u^{\|}=-u n_{y}+v n_{x}$ 分别表示边界上的法向流速和切向流速。

该边界上的扰动不会对计算域内的水流流态产生影响，故而边界右侧网格的水力要素值均采用左侧网格的值：

 $$h_{\rm{R}}=h_{\mathrm{L}}, \quad u_{\rm{R}}^{\perp}=u_{\mathrm{L}}^{\perp}, \quad u_{\rm{R}}^{\|}=u_{\mathrm{L}}^{\|}$$ (8)

通常在水位边界上给定一个水位变化过程z=z(t)，可算出边界右侧每个网格上的水深hR；而边界左侧网格上水深和流速为已知，根据1维特征线理论：

 $$u_{\mathrm{R}}^{\perp}+2 c_{\mathrm{R}}=u_{\mathrm{L}}^{\perp}+2 c_{\mathrm{L}}$$ (9)

式中，波速 $c=\sqrt{g h}$

由式（9）计算得： $u_{\mathrm{R}}^{\perp}=u_{\mathrm{L}}^{\perp}+2 c_{\mathrm{L}}-2 c_{\mathrm{R}}$ 。认为两侧网格上切向流速相等，即 $u_{\mathrm{R}}^{\|}=u_{\mathrm{L}}^{\mathrm{L}}$

当右侧网格上单宽流量 $q_{\mathrm{R}}^{\perp}$ 已知，由 $q_{\mathrm{R}}^{\perp}=h_{\mathrm{R}} u_{\mathrm{R}}^{\perp}$ 和式（9）联立可得方程：

 $$2 c_{\mathrm{R}}^{3}-c_{\mathrm{R}}^{2}\left(u_{\mathrm{L}}^{\perp}+2 c_{\mathrm{L}}\right)+g q_{\mathrm{R}}^{\perp}=0$$ (10)

采用牛顿–拉夫逊方法迭代求解 $c_{{\rm{R}}}$ ，最终解得 $h_{\mathrm{R}}=c_{\mathrm{R}}^{2} / \mathrm{g}$ ，然后由已知的单宽流量计算得出法向流速：

 $$u_{\rm{R}}^{\perp}=q_{\rm{R}}^{\perp} / h_{\rm{R}}$$ (11)

对于一般的河道洪水模拟来说，通常都是在上游入口处给一个断面流量过程，但是如何合理的分配流量就成了一个问题。模型采用三角形非结构网格剖分计算域，其优势在于能够灵活处理复杂地形边界，所以模型中设置为垂直于入流边界入流。最简单的处理办法就是用流量除以入口宽度来计算单宽流量：

 $$q_{i}^{\perp}=Q^{\perp} / B$$ (12)

式中， $Q^{\perp}$ 为入口处断面流量，B为入口宽度。

由于在入流边界附近，该处理方式对水流流态有一定影响，故而本文对入流模块做了部分改进，采用一种基于地形网格数据的纵断面提取方法，确定河道纵比降，在入流边界网格上以均匀流的方式入流，确定边界上的水位，进而分配边界各网格上的水力要素，最后耦合至模型通量部分进行计算。

对于河流中的某一段来说，河底高程沿纵断面变化趋势大概一致；再考虑到实际地形获取以及网格剖分时产生的误差，导致局部河道纵比降存在较大的不确定性，故选用河道平均纵比降计算。

沿深泓线的剖面称为河道纵断面，能反映河床的沿程变化。在实际工作中，通常以河槽底部转折点高程为纵坐标，以河长为横坐标，绘制出河道纵断面图。本文引入一种基于地形网格数据的河道纵断面提取方法，从入流边界开始，首先，确定入流边界处高程最低的网格，以此为起点和搜索中心，自动搜索所在网格附近的相邻网格，找出起点网格下游处高程最低的网格编号；然后，以该网格为搜索中心，继续重复上述搜索步骤，并依次向下游搜索，直到出流边界为止，搜索完成，所有搜索到的最低高程网格中心点依次连线近似为深泓线的位置。以Toce河地形为例，如图3所示，在这一段河道中，深泓线一共经过534个网格，全长59.85 m；Toce河物理模型试验的真实地形如图4所示。通过对比图3的深泓线与图4试验地形上深泓线可以看出，该深泓线的位置较为准确。

图  3  Toce河测点分布及深泓线示意图
Fig.  3  View of the gauging point and thalweg of Toce river
图  4  Toce河物理模型试验地形[24]
Fig.  4  Physical model topography of Toce river[24]

按照河道平均比降计算方法，确定河道纵比降，考虑到从上游入流，上游附近区域的局部地形对入口处流态有一定影响，故选取上游断面代替原方法的下游断面河底高程为基点做一斜线，使得斜线以下的面积与原河底线以下面积相等，该斜线坡度即为河道的平均纵比降，计算式为：

 $$J = \frac{2 {\textit{z}}_{\mathrm{b} 0} L-\left(\left({\textit{z}}_{\mathrm{b} 0}+{\textit{z}}_{\mathrm{b} 1}\right) l_{1}+\left({\textit{z}}_{\mathrm{b} 1}+{\textit{z}}_{\mathrm{b} 2}\right) l_{2}+{\cdots}+\left({\textit{z}}_{\mathrm{b} m-1}+{\textit{z}}_{\mathrm{b} m}\right) l_{{m}}\right)}{L^{2}}$$ (13)

式中，zb0zb1 $\cdots、$ zbm 分别为自上游到下游沿程各网格中心点高程，l1l2 $\cdots、$ lm 分别为相邻网格中心点间的距离，L为河段全长。

图5为河底高程沿程变化及平均纵比降示意图，其中，J1为原方法计算出的比降，J2为本文调整后计算出的比降。从图5可以看出：对于Toce河中的某一段来说，河底高程变化趋势确实基本一致，但在下游河道出口附近变化较为剧大。与整段河道变化趋势相比，以上游断面河底高程为基点的平均纵比降计算方法更为合理，由式（13）计算出Toce河该段平均比降约为1.32%。

图  5  Toce河道纵断面及平均纵比降示意图
Fig.  5  Vertical profile and average gradient of Toce river

在入流边界各网格上以均匀流的方式入流，通过给定断面流量计算出所对应的水位。认为每个网格上水力半径Ri=hiJ统一取平均比降，则每个网格上的流量为：

 $$Q_{i}^{\perp}=A_{i} C \sqrt{R_{i} J}=\frac{1}{n} l_{i} h_{i}^{5 / 3} J^{1 / 2}$$ (14)

式中，Ai为入流边各网格上过水面积，C为谢才系数。假设入流边上的水位为z时，则：

 $$h_{i}={\textit{z}}-{\textit{z}}_ {\mathrm{b} i}$$ (15)

总断面流量可表示为：

 $$Q^{\perp}=\sum_{k=1}^{m} Q_{k}^{\perp}=\frac{1}{n} J^{1 / 2} \sum_{k=1}^{m} l_{k}\left({\textit{z}}-{\textit{z}}_{\mathrm{b} k}\right)^{5 / 3}$$ (16)

用二分法逼近求解式（16）的近似解，可计算出入流边上的水位，从而得到入流边界上的水深。误差估计为 $|Q^{\perp}-Q_{\text{*}}^{\perp}| < 1 \times 10^{-6}$ 。此时，并不直接按式（14）计算每个网格上的流量，而是通过流量分配方法以确保流量守恒。

根据谢才公式和曼宁公式，流速近似与水深的2/3次方成正比，则在入流边界第1至m个网格上[11]

 $$\frac{u_{1}^{\perp}}{h_{1}^{2 / 3}}=\frac{u_{2}^{\perp}}{h_{2}^{2 / 3}}=\frac{u_{3}^{\perp}}{h_{3}^{2 / 3}}={\cdots}=\frac{u_{m}^{\perp}}{h_{m}^{2 / 3}}$$ (17)

当给定上游入流断面流量 $Q^{\perp}$ 时，

 $$Q^{\perp}=\sum_{k=1}^{m}(h_{k} u_{k}^{\perp} l_{k})$$ (18)

式中，lk为第k条入流边的长度，可得：

 $$Q^{\perp}=\sum_{k=1}^{m}\left(h_{k} u_{k}^{\perp} l_{k}\right)=\left(\frac{u_{i}^{\perp}}{h_{i}^{2 / 3}}\right) \sum_{k=1}^{m}(h_{k}^{5 / 3} l_{k})$$ (19)

则第i条边的流量为：

 $$Q_{i}^{\perp}=h_{i} u_{i}^{\perp} l_{i}=\frac{Q^{\perp}(h_{i}^{5 / 3} l_{i})}{\displaystyle\sum_{k=1}^{m}(h_{k}^{5 / 3} l_{k})}$$ (20)

则第i条边的单宽流量为：

 $$q_{i}^{\perp}=\frac{Q_{i}^{\perp}}{l_{i}}=h_{i} u_{i}^{\perp}=\frac{Q^{\perp} h_{i}^{5 / 3}}{\displaystyle\sum_{k=1}^{m}(h_{k}^{5 / 3} l_{k})}$$ (21)

然后，结合单宽流量边界条件迭代求解边界右侧网格上的水深和流速，进而代入HLLC近似黎曼求解器计算界面通量。

以Toce河物理模型试验模拟为例，对入流边界采用了4种不同的处理方法，见表1。比较各方法所对应的入口流态，确定一种合理的入流边界处理方法。入口处局部网格划分如图6所示，流量过程如图7所示。

表  1  入流边界不同处理方法
Table  1  Different processing method of inflow boundary
 方法 入口宽度 /m 流量分配 方法 右侧网格 信息构造 M1 3.407 式（12） 式（10）+（11） M2 1.703 式（12） 式（10）+（11） M3 3.407 式（16）+（21） 式（10）+（11） M4 3.407 式（16）+（21） 式（15）+（11）
图  6  入口附近局部网格划分示意图
Fig.  6  View of local mesh generation near the inlet
图  7  Toce河入流流量过程
Fig.  7  Flow discharges of Toce river

按照流量除以入口宽度计算单宽流量，方法M1：入口宽度B=3.407 m，在入口处水位不平顺，出现了高处过水的现象，经分析是由于入口宽度太宽的原因导致的；方法M2：缩短入流边界宽度B=1.703 m，发现只在入流边界处入流，水位在两侧非入流边界处存在间断，结果如图89中入流边界内部网格上的水位和流速分布所示。而对于河流横断面来说，水位是连续的，这种高处过水和水位存在间断的现象极不合理。

图  8  入流边界上水位分布
Fig.  8  Water level distribution on inflow boundary
图  9  入流边界上流速分布
Fig.  9  Velocity distribution on inflow boundary

对于不同的入流流量来说，其过水断面的水面宽度也是不同的，这就需要根据不同流量合理地确定水面宽度，所以本文采用了一种新的水面宽度自适应确定方法。方法M3：在每个网格上以均匀流的方式入流，用二分法逼近求解对应流量下的限定水位，用该水位计算出的水深去分配入流边界上每个网格上的流量（高出限定水位的网格上不分流量），以确保流量守恒。再结合单宽流量边界条件迭代求解边界右侧网格上的水力要素，进而计算界面通量。方法M4：从方法M3已求得的入流边界右侧网格上的水位和单宽流量，不用迭代，就可直接通过式（11）计算流速，进而计算通量。

在入流峰值时刻t= 20 s时，采用不同处理方法所对应入流边界内部网格上的水位和流速分布如图89所示，入流边界附近局部区域的水深分布如图10 所示。

图  10  入流边界附近局部区域水深分布
Fig.  10  Water depth distribution in local area near the inlet

图810可以看出，当在模型中设置一个固定宽度的入流边界时，对于不同的入流流量来说，若简单的按照所有入流边界网格上用断面输入流量除以入口宽度去计算单宽流量的方法处理，在入流边界内部网格及附近局部区域会出现像高处过水、水位间断等不合理的非物理现象。经过本文改进过的流量分配方法模拟结果则较为合理，而且方法M4与方法M3处理效果相似，对整体结果影响不大，避免了迭代求解过程，节省了计算时间。

该物理模型真实体现了Toce河上游约5 km的复杂河道地形，由ENEL在意大利米兰市按照1∶100的比尺建立，作为欧盟支持的CADAM项目的一个标准模型，常被用来检验水动力模型的精度和稳定性[24-29]。该物理模型尺寸约为50 m×11 m，地形如图4所示。采用三角形网格剖分计算域，共生成了具有42 014个节点、82 895个单元的非结构网格，平均单元面积0.004 27 m2。曼宁系数n=0.016 2 s/m1/3，模型初始为干地形，水深为0，左侧为入流边界，右侧为自由出流边界，其余均为固壁边界。在河道内布置了一系列的测点，本文选择其中部分测点来检验模型精度，其测点位置如图3所示。模拟中库朗数CFL=0.8。

图111213分别为入流时间t=30、45、60 s时，河道中洪水演进过程的水深、流速及流态分布。

图  11  Toce河洪水演进过程水深分布
Fig.  11  Water depth distribution of Toce river
图  12  Toce河洪水演进过程流速分布
Fig.  12  Velocity distribution of Toce river
图  13  Toce河洪水演进过程流态分布
Fig.  13  Flow state distribution of Toce river

图1113是采用优化方法处理后水深、流速和流态的计算结果，从图1113中可以看出：在Toce河物理模型溃坝试验洪水演进过程模拟中，急缓流交替，流态复杂，极其考验着模型的准确性和稳定性。图14为河道中不同位置测点的水位随时间变化过程。通过与观测数据的比较可以看出，水位的计算值与试验测量值总体上吻合较好，变化趋势基本一致；从测点P1到P26的洪水传播时间为42 s，与试验观测值40 s极为接近，采用纳什效率系数（NSE）来评价模型质量，经计算，河道中5个测点的纳什效率系数范围为0.82～0.95，说明该模型计算精度较高，验证了该入流边界处理方法的可行性和模型的准确性；而且在单机上仅用时86 s就完成了该算例的模拟计算，表明该模型可以高效准确的模拟洪水演进过程，对于洪水灾害管理决策意义重大。

图  14  Toce河道中各测点水位变化过程
Fig.  14  Water level evolution process of gauging points in Toce river

为了高效准确的模拟洪水演进过程，本文使用一种较为合理的入流边界流量分配方法，采用基于非结构网格有限体积法建立2维水动力模型，改进了模型中入流模块，并引入GPU加速技术对模型算法进行处理，以提高计算效率，将其应用于Toce河物理模型试验的模拟，结果表明：

本文采用一种基于地形网格数据的河道纵断面提取方法以确定河道比降，在入流边网格上通过均匀流的方式入流，确定水位，依据入流边界各网格上的水深和边长去计算权重以分配流量。从Toce河试验模拟结果可以看出，在上游入口处，水流流态更为合理。

从河道中各测点水位的模拟值与观测值的比较可以看出，水位变化趋势基本一致，洪水到达及传播时间较为接近，河道中各测点的纳什效率系数在0.82～0.95之间，验证了该方法的可行性和模型的准确性。

本文利用GPU加速技术对模型算法进行处理，在单机上用时86 s完成了具有82 895个计算单元的复杂地形上180 s入流洪水演进过程的模拟试验，计算效率较高。

通过本文对入流模块的改进与实现，该水动力模型可以成功的应用于流域洪水演进过程的高效准确模拟，为洪水灾害管理决策提供一个更直接有效的工具，更好地指导防洪减灾及预警预报工作。

• 图  1   三角形网格上各边通量

Fig.  1   Flux of each side on triangular grids

图  2   边界处理方法示意图

Fig.  2   Diagram of boundary processing method

图  3   Toce河测点分布及深泓线示意图

Fig.  3   View of the gauging point and thalweg of Toce river

图  4   Toce河物理模型试验地形[24]

Fig.  4   Physical model topography of Toce river[24]

图  5   Toce河道纵断面及平均纵比降示意图

Fig.  5   Vertical profile and average gradient of Toce river

图  6   入口附近局部网格划分示意图

Fig.  6   View of local mesh generation near the inlet

图  7   Toce河入流流量过程

Fig.  7   Flow discharges of Toce river

图  8   入流边界上水位分布

Fig.  8   Water level distribution on inflow boundary

图  9   入流边界上流速分布

Fig.  9   Velocity distribution on inflow boundary

图  10   入流边界附近局部区域水深分布

Fig.  10   Water depth distribution in local area near the inlet

图  11   Toce河洪水演进过程水深分布

Fig.  11   Water depth distribution of Toce river

图  12   Toce河洪水演进过程流速分布

Fig.  12   Velocity distribution of Toce river

图  13   Toce河洪水演进过程流态分布

Fig.  13   Flow state distribution of Toce river

图  14   Toce河道中各测点水位变化过程

Fig.  14   Water level evolution process of gauging points in Toce river

表  1   入流边界不同处理方法

Table  1   Different processing method of inflow boundary

 方法 入口宽度 /m 流量分配 方法 右侧网格 信息构造 M1 3.407 式（12） 式（10）+（11） M2 1.703 式（12） 式（10）+（11） M3 3.407 式（16）+（21） 式（10）+（11） M4 3.407 式（16）+（21） 式（15）+（11）
•  [1] 夏军,石卫.变化环境下中国水安全问题研究与展望[J].水利学报,2016,47(3):292–301. Xia Jun,Shi Wei.Perspective on water security issue of changing environment in China[J].Journal of Hydraulic Engineering,2016,47(3):292–301 [2] 夏军,石卫,张利平,等.气候变化对防洪安全影响研究面临的机遇与挑战[J].四川大学学报(工程科学版),2016,48(2):7–13. Xia Jun,Shi Wei,Zhang Liping,et al.Opportunity and challenge of the climate change impact on flood protection[J].Journal of Sichuan University(Engineering Science),2016,48(2):7–13 [3] Bui A T,Dungey M,Nguyen C V,et al.The impact of natural disasters on household income,expenditure,poverty and inequality:Evidence from Vietnam[J].Applied Economics,2014,46(15):1751–1766. [4] 何书会,杨慧英,杨艳玲.水文及水力学数学模型[J].河北水利水电技术,2003,27(1):24–26. He Shuhui,Yang Huiying,Yang Yanling.Mathematical model of hydrology and hydraulics[J].Hebei Water Resources and Hydropower Engineering,2003,27(1):24–26 [5] 周月华,彭涛,史瑞琴.我国暴雨洪涝灾害风险评估研究进展[J].暴雨灾害,2019,38(5):494–501. Zhou Yuehua,Peng Tao,Shi Ruiqin.Research progress on risk assessment of heavy rainfall and flood disasters in China[J].Torrential Rain and Disasters,2019,38(5):494–501 [6] 徐祖信,尹海龙.黄浦江水环境模拟计算边界条件影响分析[J].同济大学学报(自然科学版),2006,34(1):74–79. Xu Zuxin,Yin Hailong.Analysis of boundary conditions on hydrodynamic and water quality simulation for the Huangpu river[J].Journal of Tongji University(Natural Science),2006,34(1):74–79 [7] 吴红侠.基于有限体积法的二维浅水水流数值模拟技术研究[D].合肥:安徽大学,2012. Wu Hongxia.Numerical simulation on two-dimensional shallow-water flows based on finite volume method[D].Hefei:Anhui University,2012 [8] 章卫军,廖青桃,杨森,等.从郑州“2021.7.20”水灾模型推演看城市洪涝风险管理[J].中国防汛抗旱,2021,31(9):1–4. Zhang Weijun,Liao Qingtao,Yang Sen,et al.Thoughts and inspirations urban flood risk management inferred from Zhengzhou flood model[J].China Flood & Drought Management,2021,31(9):1–4 [9] Gallegos H A,Schubert J E,Sanders B F.Two-dimensional,high-resolution modeling of urban dam-break flooding:A case study of Baldwin Hills,California[J].Advances in Water Resources,2009,32(8):1323–1335. [10] Singh J,Altinakar M S,Ding Y.Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme[J].Advances in Water Resources,2011,34(10):1366–1375. [11] 于普兵.二维浅水水流数值模拟技术研究[D].南京:南京水利科学研究院,2006. Yu Pubing.Numerical simulation of two-dimensional shallow water flows[D].Nanjing:Nanjing Hydraulic Research Institute,2006. [12] 宋利祥.溃坝洪水数学模型及水动力学特性研究[D].武汉:华中科技大学,2012. Song Lixiang.Research on mathematical model and hydrodynamic characteristics of dam-break floods[D].Wuhan:Huazhong University of Science & Technology,2012. [13] 张大伟.堤坝溃决水流数学模型及其应用研究[D].北京:清华大学,2008. Zhang Dawei.Research on numerical model of levee breach and dam break water flow and its application[D].Beijing:Tsinghua University,2008. [14] 侯精明,张兆安,马利平,等.基于GPU加速技术的非结构流域雨洪数值模型[J].水科学进展,2021,32(4):567–576. Hou Jingming,Zhang Zhaoan,Ma Liping,et al.Unstructured numerical model for rainfall-runoff process in watershed based on GPU acceleration technology[J].Advances in Water Science,2021,32(4):567–576 [15] Hou Jingming,Liang Qiuhua,Simons Franz,et al.A 2D well-balanced shallow flow model for unstructured grids with novel slope source term treatment[J].Advances in Water Resources,2013,52(1):107–131. [16] Hou Jingming,Wang Tian,Li Peng,et al.An implicit friction source term treatment for overland flow simulation using shallow water flow model[J].Journal of Hydrology,2018,564:357–366. [17] 龚佳辉,侯精明,薛阳,等.城市雨洪过程模拟GPU加速计算效率研究[J].环境工程,2020,38(4):164–169. Gong Jiahui,Hou Jingming,Xue Yang,et al.Computational efficiency of GPU acceleration in numerical simulation of urban rain-flood process[J].Environmental Engineering,2020,38(4):164–169 [18] Kuiry S N,Pramanik K,Sen D.Finite volume model for shallow water equations with improved treatment of source terms[J].Journal of Hydraulic Engineering,2008,134(2):231–242. [19] Anastasiou K,Chan C T.Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes[J].International Journal for Numerical Methods in Fluids,1997,24(11):1225–1245. [20] Liang Q,Borthwick A G L.Adaptive quadtree simulation of shallow flows with wet–dry fronts over complex topography[J].Computers & Fluids,2009,38(2):221–234. [21] Sleigh P A,Gaskell P H,Berzins M,et al.An unstructured finite-volume algorithm for predicting flow in rivers and estuaries[J].Computers & Fluids,1998,27(4):479–508. [22] Dadone A,Grossman B.Ghost-cell method for inviscid two-dimensional flows on cartesian grids[J].AIAA Journal,2004,42(12):2499–2507. doi: 10.2514/1.697 [23] 马利平.二维地表水动力数值模拟内边界处理方法研究[D].西安:西安理工大学,2020. Ma Liping.Study on solving approach of inner boundary in two-dimension surface hydrodynamic numerical simulation[D].Xi’an:Xi’an University of Technology,2020. [24] Costanza A,Pasquale F,Marco S,et al.The FLO diffusive 1D-2D model for simulation of river flooding[J].Water,2016,8(5):1–29. [25] Prestininzi P.Suitability of the diffusive model for dam break simulation:Application to a CADAM experiment[J].Journal of Hydrology,2008,361(1):172–185. [26] Lai W,Khan A A.A discontinuous Galerkin method for two-dimensional shallow water flows[J].International Journal for Numerical Methods in Fluids,2012,70(8):939–960. doi: 10.1002/fld.2721 [27] Abderrezzak K,Paquier A,Mignot E.Modelling flash flood propagation in urban areas using a two-dimensional numerical model[J].Natural Hazards,2009,50(3):433–460. [28] 张大伟,李丹勋,王兴奎.基于非结构网格的溃坝水流干湿变化过程数值模拟[J].水力发电学报,2008,27(5):98–102. Zhang Dawei,Li Danxun,Wang Xingkui.Numerical modeling of dam-break water flow with wetting and drying change based on unstructured grids[J].Journal of Hydroelectric Engineering,2008,27(5):98–102 [29] 侯精明,王俊珲,同玉,等.基于非均匀网格的高效高精度洪涝过程模拟[J].工程科学与技术,2021,53(4):53–62. Hou Jingming,Wang Junhui,Tong Yu,et al.High-efficient and high-precision flood process simulation based on the non-uniform grid[J].Advanced Engineering Sciences,2021,53(4):53–62

osid

/

• 分享
• 用微信扫码二维码

分享至好友和朋友圈