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

2022-07-27侯精明张兆安李丙尧王俊珲张大伟

工程科学与技术 2022年4期

侯精明，汪煜，张兆安，李丙尧，王俊珲，张大伟

(1.西安理工大学西北旱区生态水利国家重点实验室，陕西西安 710048；2.河海大学水文水资源学院，江苏南京 210024；3.中国水利水电科学研究院，北京 100038)

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

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

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

1 2维水动力模型

1.1 控制方程

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

式（1）～（2）中：

为水深；

和

分别为

、

方向的单宽流量；

为重力加速度；

、

分别为

、

方向的流速；

和

分别为

、

方向的通量矢量；

为源项矢量；

为河床底高程；

为床面摩擦系数，

，其中

为曼宁系数。

1.2 数值计算方法

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

内，其积分形式可写为：

式中，

为控制体

的面积。

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

式中，

为控制体

的边界，

为边界

所对应外法线方向的单位向量。相应界面通量

(

)

·n

可以表示为：

式中，

和

为边界外法向单位向量

在

和

方向的分量。在三角形非结构网格

中，各通量如图1所示，

、

分别为

的相邻网格单元，则单元

内通量项线积分可写为：

图1 三角形网格上各边通量Fig. 1 Flux of each side on triangular grids

式中，

为网格边的编号，

为第

个网格单元上第

条边的长度。通量项采用HLLC近似黎曼解算器求解，底坡源项和摩阻源项的具体计算方法见文献[15–16]，本文不再赘述。此外，引入GPU加速技术以大幅提高模型计算效率。

2 边界条件

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

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

2.1 固壁边界条件

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

式中，下标R表示的是右侧虚拟网格，下标L表示的是左侧内部网格，

和

分别表示边界上的法向流速和切向流速。

2.2 自由出流边界条件

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

2.3 水位边界条件

通常在水位边界上给定一个水位变化过程

(

)，可算出边界右侧每个网格上的水深

；而边界左侧网格上水深和流速为已知，根据1维特征线理论：

2.4 单宽流量边界条件

3 入流边界–断面流量边界条件

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

式中，

为入口处断面流量，

为入口宽度。

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

3.1 河道纵断面及纵比降确定

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

沿深泓线的剖面称为河道纵断面，能反映河床的沿程变化。在实际工作中，通常以河槽底部转折点高程为纵坐标，以河长为横坐标，绘制出河道纵断面图。本文引入一种基于地形网格数据的河道纵断面提取方法，从入流边界开始，首先，确定入流边界处高程最低的网格，以此为起点和搜索中心，自动搜索所在网格附近的相邻网格，找出起点网格下游处高程最低的网格编号；然后，以该网格为搜索中心，继续重复上述搜索步骤，并依次向下游搜索，直到出流边界为止，搜索完成，所有搜索到的最低高程网格中心点依次连线近似为深泓线的位置。以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]

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

图5为河底高程沿程变化及平均纵比降示意图，其中，

为原方法计算出的比降，

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

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

3.2 均匀流入流

在入流边界各网格上以均匀流的方式入流，通过给定断面流量计算出所对应的水位。认为每个网格上水力半径

，

统一取平均比降，则每个网格上的流量为：

式中，

为入流边各网格上过水面积，

为谢才系数。假设入流边上的水位为

时，则：

总断面流量可表示为：

3.3 入口流量分配

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

个网格上：

式中，

为第

条入流边的长度，可得：

则第

条边的流量为：

则第

条边的单宽流量为：

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

4 入流边界不同处理方法讨论

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

图7 Toce河入流流量过程Fig. 7 Flow discharges of Toce river

表1 入流边界不同处理方法
Tab. 1 Different processing method of inflow boundary

方法入口宽度/m流量分配方法右侧网格信息构造M13.407式（12）式（10）+（11）M21.703式（12）式（10）+（11）M33.407式（16）+（21）式（10）+（11）M43.407式（16）+（21）式（15）+（11）

图6 入口附近局部网格划分示意图Fig. 6 View of local mesh generation near the inlet

按照流量除以入口宽度计算单宽流量，方法M1：入口宽度

=3.407 m，在入口处水位不平顺，出现了高处过水的现象，经分析是由于入口宽度太宽的原因导致的；方法M2：缩短入流边界宽度

=1.703 m，发现只在入流边界处入流，水位在两侧非入流边界处存在间断，结果如图8和9中入流边界内部网格上的水位和流速分布所示。而对于河流横断面来说，水位是连续的，这种高处过水和水位存在间断的现象极不合理。

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

图9 入流边界上流速分布Fig. 9 Velocity distribution on inflow boundary

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

在入流峰值时刻

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

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

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

5 模型应用

该物理模型真实体现了Toce河上游约5 km的复杂河道地形，由ENEL在意大利米兰市按照1∶100的比尺建立，作为欧盟支持的CADAM项目的一个标准模型，常被用来检验水动力模型的精度和稳定性。该物理模型尺寸约为50 m×11 m，地形如图4所示。采用三角形网格剖分计算域，共生成了具有42 014个节点、82 895个单元的非结构网格，平均单元面积0.004 27 m。曼宁系数

=0.016 2 s/m，模型初始为干地形，水深为0，左侧为入流边界，右侧为自由出流边界，其余均为固壁边界。在河道内布置了一系列的测点，本文选择其中部分测点来检验模型精度，其测点位置如图3所示。模拟中库朗数CFL=0.8。图11、12和13分别为入流时间

=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

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