基于有限体积法对洪水汇合流流场的数值模拟

Numerical Simulation of Flood Confluence Flow Field Based on Finite Volume Method

摘要: 本研究采用有限体积法，通过数值计算，模拟分析了行洪、蓄洪工程退洪、和分流洪水汇合的复杂流场分布，分析了洪闸附近的汇合流场结构，对比洪闸工程修建前、后的流场变化，提出流场分布特征和流速参数，重点分析汇合流场对蓄洪工程建筑物和附近大堤的影响，提出退洪区、分洪区附近的水力学指标，为洪闸布置、堤坝防护提供科学依据。

Abstract: In this study, the finite volume method was used to simulate and analyze the complex flow field distribution of flood discharge and flood discharge project, and the confluence flow field near the flood gate was analyzed. After the flow field changes, the flow field distribution characteristics and flow velocity parameters are proposed. The influence of confluent flow field on the flood storage project buildings and nearby levees is analyzed. The hydraulic indexes near the flood control area and the flood diversion area are proposed, which are the flood gate layout and dam. Protection provides a scientific basis.

1. 引言

天然河道中的水流多属缓流，流量也比较均匀。在发生洪水时，行洪的流量以及流速增大且分布不均匀。尤其在河道中修建了洪闸等蓄洪工程建筑物后，加大了上下游洪水的势能落差，在通过洪闸后洪水的比能也相较于下游河道中水流大得多，这样巨大的能量，若不采取消能措施，会对下游河床以及堤坝具有更强的冲击力，这样反复的作用使得坝口及分叉河道汇流处流场十分复杂。

多股洪水汇合流问题最终简化为求解控制水流运动的非恒定流的拟线性双曲型偏微分方程组。在洪水传播过程中，垂直尺度变化较小，对Navier-Stokes方程沿垂直水深作积分平均，即可得到洪水汇合流的二维浅水波模型，也称为圣维南方程组。

对于浅水波的研究，主要有有限差分法 [1] 、有限元法 [2] [3] [4] 、有限体积法 [5] [6] [7] [8] 、混合分析法 [9] 。空气动力学高精度高分辨率格式也被引入到计算水动力学的研究中，例如Godunov型有限体积格式 [10] ，高分辨率的ENO格式 [11] 。本文借助格林公式对地形源项进行处理，建立了一种和谐的WAF型有限体积格式，模拟了洪水传播过程中带激波和干湿边界的浅水波方程，以几个工况为例进行模拟，对比实际数据，研究了洪水波演进、汇合流复杂流场分布以及不规则河床地形河道形状对洪水波计算的影响。

2. 数学模型

2.1. 控制方程

通过建立蓄洪区退洪、长江河道行洪和分洪区分洪的数学物理方程，采用二维浅水波方程组研究多股洪水汇合流问题。浅水波方程组采用以水位和单宽流量为变量的守恒形式，即：

$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=S$ (1)

式中t代表时间，x和y分别表示空间变量。U表示物理守恒变量， $F\left(U\right)$ 、 $G\left(U\right)$ 分别为x和y方向的通量，S表示源项：

$U=\left(\begin{array}{c}H\\ q_x\\ q_y\end{array}\right)，F=\left(\begin{array}{c}{q}_x\\ \frac{q_x^2}{h}+\frac{1}{2}g{H}^2-gbH\\ \frac{q_x{q}_y}{h}\end{array}\right)$

$G=\left(\begin{array}{c}{q}_y\\ \frac{q_x{q}_y}{h}\\ \frac{q_y^2}{h}+\frac{1}{2}g{H}^2-gbH\end{array}\right)$

$S=S_b+S_f=\left(\begin{array}{c}0\\ -gH\frac{\partial b}{\partial x}\\ -gH\frac{\partial b}{\partial y}\end{array}\right)+\left(\begin{array}{c}0\\ -\frac{g}{\rho }{\tau }_x\\ -\frac{g}{\rho }{\tau }_y\end{array}\right)$

其中 $q_x=hu$ ， $q_y=hv$ 表示x，y方向的单宽流量，u，v分别表示x，y两个方向上的流速， $H=h+b$ 为自由流表面高程，h为水深， $b\left(x,y\right)$ 为河道地形函数。源项S由两部分构成：则表示地形底坡； $S_f$ 表示摩阻底坡，根据经验公式给出： ${\tau }_x=\rho C_{f}u\sqrt{u^2+v^2}$ ， ${\tau }_y=\rho C_{f}y\sqrt{u^2+v^2}$ 。其中， $C_f=g{n}_m^2/h^{1/3}$ 为谢才系数， $n_m$ 为曼宁系数。为简化模型，本文不考虑与地球自转相关的科氏力。

2.2. 计算方法

本研究采用三角形网格处理计算区域，计算方法采用有限体积法模拟区域水位、水流、流量和流态。考虑到问题物理区域的复杂性，采用三角网格离散计算区域，三角形网格和物理量布置如图1a所示。数值计算方法采用目前比较成熟的方法，可满足计算精度和稳定性要求。计算方法构造如下：

引入算子 $\nabla =\frac{\partial }{\partial t}i+\frac{\partial }{\partial t}j$ 将浅水波方程改写成向量形式：

$\frac{\partial U}{\partial t}+\nabla \cdot E=S$ (2)

其中 $E={\left(F,G\right)}^{\text{T}}$ ，对上式在控制区域上积分得：

$\frac{\partial }{\partial t}{\displaystyle {\iint }_{V}U\text{d}V}+{\displaystyle {\iint }_V\nabla \cdot E\text{d}V}={\displaystyle {\iint }_{V}S\left(x,y,U\right)\text{d}V}$

利用格林公式得：

$\frac{\partial }{\partial t}{\displaystyle {\iint }_{V}U\left(x,t\right)\text{d}V}+{\displaystyle {\oint }_{V}E\cdot n^{\text{T}}\text{d}\zeta }={\displaystyle {\iint }_{V}S\left(x,y,U\right)\text{d}V}$

其中， $n={\left(\cos \theta ,\sin \theta \right)}^{\top }$ 为单元边界的单位外法向量。半离散有限体积法：

$\frac{\partial U_i}{\partial t}+\frac{1}{\left|V_i\right|}{\displaystyle \underset{j=1}{\overset{3}{\sum }}{\displaystyle {\int }_{L_j}\cos {\theta }_{j}F\left(U\right)+\sin {\theta }_{j}G\left(U\right)\text{d}\varsigma }}=S_i$ (3)

模拟的是守恒量的单元积分平均值：

$U_i=\frac{1}{\left|V_i\right|}{\displaystyle {\iint }_{V_j}U\left(x,y,t\right)\text{d}V}$

其中 $\left|V_i\right|$ 为计算单元 $V_i$ 的面积， $S_i=\frac{1}{\left|V_i\right|}{\displaystyle {\iint }_{V_i}S\left(x,y,U\right)\text{d}V}$ 为源项为单元平均。这样，问题等价于给出单

元法向通量和源项的近似，进而求解一个关于时间的常微分方程。根据二维浅水波方程组法向通量的旋转不变性，可以将转化成一个一维的Riemann问题近似解 [12] (图1b)：

(a) (b)

Figure 1. The structure of approximate solution of Riemann problem

图1. Riemann问题近似解的结构

${\tilde{U}}_t+F{\left(\tilde{U}\right)}_x=0$

$\tilde{U}\left(\tilde{x},0\right)=\{\begin{array}{l}{\tilde{U}}_L\tilde{x}<0\\ {\tilde{U}}_R\tilde{x}>0\end{array}$

旋转矩阵 $T\left(\theta \right)$ 定义为：

$T\left(\theta \right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \theta & \sin \theta \\ 0& -\sin \theta & \cos \theta \end{array}\right)$ (4)

其中， $\tilde{U}$ 为U旋转后的量，通量函数 $F\left(\tilde{U}\right)$ 满足 $E\left(U\right)\cdot n=T^{-1}\left(\theta \right)F\left(\tilde{U}\right)$ 。进一步，采用时间向前差分离散得：

$U_i^{n+1}=U_i^n-\frac{\Delta t}{\left|V_i\right|}{\displaystyle \underset{j=1}{\overset{3}{\sum }}{T}^{-1}\left({\theta }_{i,j}\right)F\left({\tilde{U}}_{i,j}\right)\left|L_{i,j}\right|}+\Delta t{S}_i$ (5)

上式中源项采用中心差分进行离散，基于局部Riemann问题自相似解波系结构分析，采用TVD-WAF格式给出数值通量为：

$F^{WAF}=\frac{1}{2}\left(c_4{F}_R+c_0{F}_L\right)-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{N}{\sum }}sign\left(c_k\right)\varphi \left(r^k,\left|c_k\right|\right)\Delta F^k}$ (6)

其中 $\Delta F^k$ 为第k个波处的数值通量的跳量， $c_k$ 为波 $S_k$ 的Courant数：

(7)

通量限制器 $\varphi \left(r^{\left(k\right)},c_k\right)=\{\begin{array}{ll}1\hfill & r^{\left(k\right)}\le 0\hfill \\ 1-\frac{\left(1-\left|c_k\right|\right)r^{\left(k\right)}\left(1+r^{\left(k\right)}\right)}{1+r^{\left(k\right)}\times r^{\left(k\right)}}\hfill & r^{\left(k\right)}>0\hfill \end{array}$ ，其中变量 $r^{\left(k\right)}$ 为：

$r^{\left(k\right)}=\{\begin{array}{cc}\frac{p_i^{\left(k\right)}-p_{i-1}^{\left(k\right)}}{p_{i+1}^{(k)}-p_i^{(k)}}& S_k\ge 0\\ \frac{p_{i+2}^{\left(k\right)}-p_{i+1}^{\left(k\right)}}{p_{i+1}^{\left(k\right)}-p_i^{\left(k\right)}}& S_k<0\end{array}$ (8)

上式中变量p为 $p^{\left(k\right)}=H$ ，， $p^{\left(2\right)}=q_y$ 。式(7)中波系间的数值通量可以有HLLC近似里面求解器给出：

${\left[F^{hllc}\right]}_{\left(k\right)}={\left[F^{hll}\right]}_{\left(k\right)},k=1,2，{\left[F^{hllc}\right]}_{\left(3\right)}=\{\begin{array}{cc}{\left[F^{hll}\right]}_{\left(1\right)}{v}_L& S_2>0\\ {\left[F^{hll}\right]}_{\left(1\right)}{v}_R& S_2<0\end{array}$ (9)

下标 $\left(k\right)$ 表示通量向量的第k个分量，其中HLL求解器定义为：

$F^{hll}=\{\begin{array}{ll}{F}_L\hfill & S_1\ge 0\hfill \\ \frac{S_2{F}_L-S_1{F}_R+S_1{S}_2\left(U_R-U_L\right)}{S_2-S_1}\hfill & S_1<0<S_2\hfill \\ F_R\hfill & S_2\le 0\hfill \end{array}$ (10)

其中，左右波速度 $S_1$ 和 $S_3$ 分别为 [13] ：

$S_1=\min \left({\tilde{u}}_L-\sqrt{g{h}_L},{\tilde{u}}_s-\sqrt{g{h}_s}\right)，S_3=\max \left({\tilde{u}}_R+\sqrt{g{h}_R},{\tilde{u}}_s+\sqrt{g{h}_s}\right)$ (11)

其中 ${\tilde{u}}_s=\left({\tilde{u}}_L+{\tilde{u}}_R\right)/2+\sqrt{g{h}_L}-\sqrt{g{h}_R}$ ， $\sqrt{g{h}_s}=\left(\sqrt{g{h}_L}+\sqrt{g{h}_R}\right)/2+\left({\tilde{u}}_L-{\tilde{u}}_R\right)/4$ 。接触间断波波速 $S_2$ 的估计：

$S_2=\frac{S_1{h}_R\left({\tilde{u}}_R-S_3\right)-S_3{h}_L\left({\tilde{u}}_L-S_1\right)}{h_R\left({\tilde{u}}_R-S_3\right)-h_L\left({\tilde{u}}_L-S_1\right)}$ (12)

数值计算需要处理动边界问题，及计算单位处理干湿单元的问题。当计算单元边界为干湿边界时，左干河床 $S_1={\tilde{u}}_R-2\sqrt{g{h}_R}$ ， $S_3={\tilde{u}}_R+\sqrt{g{h}_R}$ ；右干河床： $S_1={\tilde{u}}_L-\sqrt{g{h}_L}$ ， $S_2={\tilde{u}}_L+2\sqrt{g{h}_L}$ 。 $S_3$ 会自动和估计的干河床重合。

为了保证数值格式的稳定性，需要在进行计算时调整时间，采用：

$\Delta t=CFL\underset{1\le i\le M}{\min }\underset{1\le j\le 3}{\min }\frac{\min \left(\Delta {\tilde{x}}_{i,j,1},\Delta {\tilde{x}}_{i,j,2}\right)}{\max \left(S_{i,j,1},S_{i,j,3}\right)}$ (13)

其中 $CFL\left(0<CFL\le 1\right)$ 是Courant数。

为保证数值算法的结构的和谐性，首先用Green公式转化地形源项 $S_b$ ，然后利用中矩形公式进行离散。因此， $S_b$ 第二、三个分量为：

${\displaystyle {\iint }_{V_i}gH\nabla b\text{d}x\text{d}y}=-gH{\displaystyle \underset{j=1}{\overset{3}{\sum }}{\displaystyle {\oint }_{L_{i,j}}\left(b_j,0\right)\cdot n_j\text{d}\zeta }}=-gH{\displaystyle \underset{j=1}{\overset{3}{\sum }}{b}_j\cdot n_j\left|L_{i,j}\right|}$

可以证明该格式是和谐的 [14] 。

2.3. 计算条件

1) 计算区域：计算模型要求计算研究区域上下游需要由足够的计算长度，避免进出口条件设置的影响。结合现有河道地形资料，选取洪闸上游约4 km、下游约10 km，共计约14 km为计算区域。河道地形采用工程区附近2011年地形资料。

2) 网格剖分：考虑计算区域边界不规则形，边界条件复杂性，计算方法采用三角形剖分算法将物理区域划分为规则性三角形组成的计算区域，计算网格数为37,082个，为兼顾整个计算区域，同时为提高计算精度，在洪闸、分洪区、蓄洪区出口附近，采用了网格加密处理(如图2)。

(a)(b)

Figure 2. Computational area mesh Generation and Local encryption Generation Diagram

图2. 计算区域网格剖分与局部加密剖分图

2.4. 模型率定

模型利用洪闸附近水位~流量关系对计算河道的糙率进行了率定(表1)。通过分析糙率取值范围在0.023~0027。本研究采用河床糙率取值0.025。

Table 1. Model rate table

表1. 模型率定表

Figure 3. Distribution of flow field

图3. 流场分布

3. 数值计算

3.1. 长江与分洪汇合流计算

根据不同的计算需要设置了工况1~11。为了解蓄滞洪区工程建设前洪闸附近长江河道流场分布，考虑分洪任务要求，设置了计算工况8~11。通过该工况计算一方面掌握工程修建前，长江与分洪汇合流场分布。分流洪水在入江口处受长江主流顶托，流场向左岸偏离，混合长江洪水后进入下游河道，未直接影响长江上游河道(图3)。分洪流量大小对汇合流场结构影响不大，在较大洪水情况下，工况8和工况10流场相似，在较小洪水情况下，工况9和工况11流场分布相似。

3.2. 洪闸、长江和分流汇合洪水计算

在以上计算基础上，为研究蓄洪区退洪，长江行洪和分流复杂汇合流问题，加入东荆河分流洪水影响，设置了工况1~7。通过计算观测，三股洪水形成混合流如图4，首先上游长江洪水退洪洪水接触混合，偏向主槽，然后于内荆河排水混合，上游混合流在河道内进行调整。在东荆河入江口处，长江洪水与东荆河分流洪水汇合，最后向右转弯，进入下游河道。长江洪水量越大，对下游挤压作用越大，长江洪水大时，东荆河洪水偏向左侧较明显，洪水相对小，向左偏向减弱；蓄洪区泄洪流量越大，对上游长江洪水阻挡作用越强，洪闸左侧回流淘刷强度越大。各工况不同位置监测点流速如表2。

Figure 4. Distribution of flow field

图4. 流场分布

Table 2. Characteristic value of flow velocity at monitoring points of embankment

表2. 堤坝部分监测点流速特征值

4. 成果分析结论

1) 三股洪水同时汇合时，工程区域流场复杂，上游长江洪水同补元闸退洪洪水接触混合，受挤压偏向主槽，然后与内荆河排水混合，上游混合流在河道内进行短距离调整。在东荆河入口处，长江洪水再与东荆河分流洪水汇合，最后向右转弯，进入下游河道。

2) 通过流场分析，蓄洪区退洪过程一定程度上影响工程区长江洪水流场分布。通过不同工况数值计算，各监测点计算流速如图5。各工况测点流速大部分为0.05~0.3 m/s之间，在洪水相对较小时，堤防沿程各测点流速变化不大，长江洪水大于45,900 m³/s的工况时，在洪闸左侧长江下游堤段范围流速增加，最大值为0.31 m/s，结合流场分析，主要由于该区域为内荆河入江口处，同时叠加河道回流影响。

Figure 5. Flow velocity distribution of embankment toe in lower reaches of Yangtze River

图5. 长江下游堤脚流速分布

4. 结论

1) 蓄洪区退洪、长江行洪、东荆河分洪三股洪水汇合处流场结构非常复杂，在洪闸左侧长江下游河道出现两处较大范围的回流现象，但整体回流过程比较平稳，未出现主流偏离或主流折冲现象，整体布置比较安全。

2) 东荆河分流洪水同长江上游洪水汇合，主流偏左岸，然后流向下游，未直接影响长江上游区域，对补元闸工程、长江上游堤防影响不大。

3) 蓄洪区退洪与长江汇合流场复杂，工程附近流速较小，闸左、右两侧流速v < 0.4 m/s。蓄洪区退洪对工程附近长江下游流态有一定影响，增强了闸后长江下游堤防附近的回流强度。