第2章　淮河中游河道一维、二维耦合水动力数学模型

2.1　洪水演进数学模型

随着计算机和计算技术的不断发展，水动力数学模型发展迅速，在工程上得到了广泛应用，能较好地模拟复杂边界条件下的水流运动，并与实体模型耦合运用相互验证，相互提供边界，可以较好地解决工程实际问题。

淮河中游河势复杂，有单一河道，有分汊河道，还有行洪区和湖泊。低水期时，水流主要在河道中运动。洪水期，随着水位的上升，上游来流超过河道泄流能力时，通过分洪、滞洪等措施，使水流通过闸门、口门、漫堤等方式进入到行洪区内。为了正确地模拟上述水流的运动情况，需要根据不同的水流特征采用不同的模拟方法，以达到提高精度和效率的目的:①蓄洪洼地由于流速较小，水面比降小，主要考虑其水量调蓄的作用，采用水库型水量平衡法;②干支流河道水流运动，需要计算水流的过流能力及水位，采用一维非恒定河网、堰闸联合模型求解;③行洪区及湖泊由于流速大小及流向随空间和时间变化，采用平面二维浅水方程求解。

2.1.1　一维水流数学模型

2.1.1.1　控制方程组及离散

河网一维水动力模型的控制方程为Saint Venant方程组。

连续方程:

动量方程:

式中:t为时间坐标;x为河道沿程坐标;Q为流量;Z为水位;A为过水断面的面积;B为水面宽度;K为流量模数;g为重力加速度;q为旁侧入流流量;n为河道糙率系数;R为水力半径。

利用Abbott六点隐式差分格式离散上述控制方程组，该离散格式在每个网格点并不同时计算水位和流量，而是按顺序交替计算水位和流量，分别称为h点和Q点，如图2.1-1所示。该格式无条件稳定，可以在相当大的Courant数下保持稳定，可以取较长的时间步长以节约计算时间[1-7]。

采用如图2.1-2所示的Abbott六点隐式差分格式，连续性方程中的各项可以写为

于是连续性方程可以写为

同样动量方程中各项可以写为

于是动量方程在流量点上的差分格式为

式(2.1-3)整理后可记为

式(2.1-4)整理后可记为

以上式中:αj、βj、γj、δj分别为离散方程的系数。

2.1.1.2　定解条件

初始条件:给定初始时刻t=0时，所有计算节点Q和h值。非恒定流数值计算表明，初始条件对于计算的初期阶段会显示影响，但这种影响将随着计算时间的延伸逐步消失。

边界条件分三类:

(1)水位边界条件:边界处水位h=h(t)。

(2)流量边界条件:边界处流量Q=Q(t)。

(3)水位—流量关系边界条件:边界处Q=f(h)给定。

2.1.1.3　求解方法

如前所述，河道内任一点的水力参数Z(水位h或流量Q)与相邻的网格点的水力参数的关系可以表示为统一的线性方程:

上式中的系数可分别由式(2.1-5)或式(2.1-6)计算。

假设一河道有n个网格点，因为河道的首末网格点总是水位点，所以n是奇数。对于河网的所有网格点写出式(2.1-7)，可以得到n个线性方程:

其中第一个方程的Hus和最后一个方程中Hds分别是上、下游汊点的水位。某一河道第一个网格点的水位等于与之相连河段上游汊点的水位:α1=-1，β1=1，γ1=0，δ1=0。同样，α1=0，β1=1，γ1=-1，δ1=0。

对于单一河道，只要给出上下游水位边界，即Hus和Hds为已知，就可用消元求解方程组式(2.1-8)。对于河网问题，由方程组式(2.1-8)，通过消元法可以将河道内任意点的水力参数(水位或流量)表示为上下游汊点水位的函数:

只要先求河网各汊点的水位，就可用式(2.1-9)求解河段任意网格点的水力参数。

图2.1-3为河网汊点方程示意图，围绕汊点的控制体连续方程为

将式(2.1-10)右边第二式的三项分别以式(2.1-9)替代，可以得到

式中:H为该汊点的水位;HA，us、HB，us分别为支流A、B上游端汊点水位;HC，ds为支流C下游端汊点水位。

在式(2.1-11)中，将某个汊点水位表示为与之直接相连的河道汊点水位的线性函数。同样，对于河网所有汊点(假设为N个)，可以得到N个类似的方程(汊点方程组)。在边界水位或流量为已知的情况下，可以利用高斯消元法直接求解汊点方程组，得到各个汊点的水位，进而回代式(2.1-9)求解河道任意网格点的水位和流量。

大型稀疏矩阵求解计算时间主要取决于矩阵主对角线非零元素的宽度，可通过对河网节点进行优化编码的方法来降低汊点方程组系数矩阵的带宽，使之称为主对角元素占优的矩阵，从而方便了方程组的求解，并大大减少了计算时间。

若在河道边界节点上给出水位的时间变化过程，即h=h(t)，此时，假设边界所在河道编号为j，则边界上的汊点方程为

若在河道边界节点上给出流量的时间变化过程，即Q=Q(t)。

对如图2.1-4所示的控制体，应用连续方程可以得到:

将以式(2.1-9)代入式(2.1-13)，可以得到

若在河道边界节点上给出的是水位流量关系Q=Q(h)，其处理方法同流量边界，得到与式(2.1-14)类似的方程，只是方程中的由水位流量关系计算得到。

平原河网大多有堰、闸等水工建筑物，此时Saint Venant方程已经不再适用，必须根据堰、闸的水力学特性作特殊处理。在模型中堰、闸通常作为流量点处理，根据相邻水位点水位关系采用宽顶堰或孔口出流计算过闸流量，得到式(2.1-6)类似的方程[1-4]。

2.1.2　平面二维水流数学模型

2.1.2.1　基本方程

大部分河流湖泊都是浅水水流，满足以下假定:①具有自由表面;②以重力为主要的驱动力，以水流和固体边界之间及水流内部摩阻力为主要的耗散力;③水平流速沿水深近似成均匀分布;④垂向流速和垂向加速度可忽略，水压接近静压分布。基于以上假定，对三维水流的运动方程沿水深进行积分，可得二维浅水水流控制方程。

连续性方程:

动量方程:

式中:ξ为水位;h为水深，h=ξ-Zb，Zb为河床底高程;u、v分别为x，y方向的垂向平均流速;f为科氏力系数，f=2ωsinφ，ω为地球自转的角速度，φ为计算水域的地理纬度;τbx、τby分别为x方向和y方向的底部摩阻，=n为曼宁系数;τsx、τsy分别为风对自由表面x方向和y方向的剪切力，cd为风阻力系数，ρa为空气密度，uw，vw为x方向和y方向的风速;ρ为水密度;Ex、Ey分别为x方向和y方向的水流紊动黏性系数;q为源或汇单位面积的流量，源取正，汇取负;u*、v*分别为源或汇输入输出时在x和y方向的流速。

上述二维浅水方程可以写成如下统一的形式:

式中:

q=(h，hu，hv)

2.1.2.2　定解条件

初始条件:给定初始时刻t=0时，计算域内所有计算变量(u，v，ζ)的初始值。

u(x，y，t)t=0=u0(x，y)

v(x，y，t)t=0=v0(x，y)

ξ(x，y，t)t=0=ξ0(x，y)

边界条件:包括固壁边界和开边界。

(1)固壁边界:近壁区因为分子黏性较大，雷诺数较低，为避免在壁面处加密网格，通常采用壁面函数法和滑移边界来处理近壁处的流速。壁面函数法通过一组半经验公式直接将壁面上的物理量与湍流核心区内待求的未知量直接联系起来;滑移边界条件，即=0。本文计算采用滑移边界条件处理固壁边界。

(2)开边界:可以给定水位、流量或流速的过程。

给定水位边界:ξ=ξ(t);

给定流量边界:Q=Q(t);

给定法向流速过程:un=un(t)。

2.1.2.3　离散方法

Mike21 FM采用非结构有限体积法离散控制方程[1-4]。有限体积法中使用的非结构网格通常由三角形或四边形网格构成，为了准确的逼近水下地形，这里仅采用三角形网格。

采用网格中心式布置计算节点，见图2.1-5。将式(2.1-18)在控制体上进行积分得:

式中:Ω为控制体平面域;dA为面积分微元。采用高斯格林公式将面积分化为沿其周界的线积分，得

式中:∂Ω为控制体周界(逆时针方向);d S为线积分微元;nx、ny为周界上外法向单位向量;分别为守恒物理量、源汇项在控制体内平均，当精度小于或等于二阶时，它们分别等于变量和源、汇项在三角形形心处的值。

记跨控制体边界的法向通量Fn(q)=F(q)nx+G(q)ny，考虑到本次计算采用三角网格，故式(2.1-19)又可写为

式中:Lj为第j边的边长;Fn(q)j为第j边平均法向数值通量。

由于Fn(q)沿控制体边界一般为非线性分布，用数值求积公式计算可取得较高的精度。因此式(2.1-20)可写为

式中:为数值求积公式计算的第j边的平均法向数值通量Fn(qj，k);ωk为积分权;n为积分权总数;qj，k为守恒物理量在j边，第k积分点的值。

对于空间离散后得到的常微分方程，本文采用以下两种常微分方程解算器进行时间离散。记式(2.1-21)等号右端项为L(q)，则式(2.1-21)可写为

(1)一阶显格式:

式中:为控制体内守恒物理量的平均值，q由重构的控制体内守恒物理量，为(x，y)的函数。

(2)二阶Runge-kutta法:

至此，所有问题归结为如何求解跨界面处的数值通量。目前常用的途径是基于间断的思想，认为在控制体界面两侧的物理量存在间断，这样在每个控制界面处可以形成一个黎曼问题。可将F(qn)取为图2.1-6所示的一维黎曼问题的解。这样一维黎曼问题可以表示为

初始条件为

其中qL和qR分别为间断左右的物理量值。

Mike21 FM采用Roe格式求解以上黎曼问题，为避免数值振荡，使用了二阶TVD限制器。

2.1.3　一维、二维耦合水流数学模型

一维模型计算简便，适合模拟河道内水流运动情况;二维模型需要占用更多的计算资源，适合模拟行洪区和湖泊内水流运动情况。通过建立一维、二维耦合水流数学模型，将上述两种模型连成一体，联合求解，既可以发挥出一维模型快速方便的特点，同时又能获得局部范围的细部信息。Mikeflood包括以下三种类型的连接方式:

(1)标准连接。标准连接主要用于将一个或多个二维网格单元连接到一个一维模型的末端。这种连接方式适用于河口复杂水流的计算中，具体应用如图2.1-7所示。

(2)侧向连接。侧向连接主要用于一串二维网格侧向连接到一维河道中(一部分河道或整条河道)。这种连接方式适用于模拟水流从河道进入行洪区的过程，具体应用如图2.1-8所示。

(3)结构物连接。结构物连接主要应用于将两个分开的二维网格通过一维结构物连接成一个整体。这种连接方式适用于在二维网格中模拟结构物对水流的影响，具体应用如图2.1-9所示。

对于上述连接方式，由于连接点不属于边界点，因此不论对于一维问题、还是二维问题都缺少边界条件，无法独立求解。需要通过在连接点处补充物理量之间的关系(水位、流量相等)，从而实现一维、二维模型的耦合。

按照耦合求解方式的不同，将其分为两大类。标准连接和侧向连接采用交替计算的方法，该方法是在一个时步内把一个模型计算出来的物理量(水位或流量)作为另一个模型的边界条件，两模型交替计算，互赋边界的方法完成整个计算过程，这种方式在时间步长不大时，省去了迭代的过程，节省了计算的时间;结构物连接是采用直接求解的方法，该方法将连接点处的连接条件作为定解条件引入方程组中，通过补充迭代关系式，从而实现一维、二维的耦合计算。这种方式每一时步都需完成迭代计算，计算量大但稳定性较高，可以适用较大的时间步长。