Mike11软件包由水动力、对流~扩散、水质、降雨~径流、洪水预报等模块组成,核心模块为水动力模块。Mike11水动力模块采用6点Abbott~Ionescu有限差分格式对圣维南方程组求解。

一、圣维南方程组

1、基本要素与假设条件

Mike11模型基于以下三个要素:反映有关物理定律的微分方程组;对微分方程组进行线性化的有限差分格式;求解线性方程组的算法。并基于以下几个假定:流体为不可压缩、均质流体;一维流态; 坡降小、纵向断面变化幅度小;符合静水压力假设。

2、圣维南方程组

{∂Q∂x+∂A∂t=q∂Q∂t+∂(αQ2A)δx+gAδhδx+gQ∣Q∣C2AR=0\left\{\begin{array}{l} \frac{\partial Q}{\partial x}+\frac{\partial A}{\partial t}=q \\ \frac{\partial Q}{\partial t}+\frac{\partial\left(\alpha \frac{Q^{2}}{A}\right)}{\delta x}+g A \frac{\delta h}{\delta x}+\frac{g Q|Q|}{C^{2} A R}=0 \end{array}\right. ⎩⎨⎧​∂x∂Q​+∂t∂A​=q∂t∂Q​+δx∂(αAQ2​)​+gAδxδh​+C2ARgQ∣Q∣​=0​
式中:Q为流量,m³/s;q为侧向入流,m³/s;A为过水面积,m²;h为水位,m;R为水力半径,m;C为谢才系数;a为动量修正系数。


关于圣维南方程的方程的简洁推导可以参考:一个简洁的圣维南方程组推导过程,根据《水文学原理》,描述洪水波在河道中传播规律的物理途径主要有两类:一是利用圣维南方程组描述河道洪水波运动;二是利用河段水量平衡方程式和槽蓄方程式描述河道洪水波运动。基于前一种途径的洪水演算方法称为水力学法。基于后一种途径的洪水演算方法称为水文学法。

3、公式推导


令水的密度为ρ\rhoρ ;重力加速庶为 ggg, h(x,t)h(x, t)h(x,t) 为水在时间ttt位置xxx的深 (高) 度; $v(x, t) $ 为在时间ttt,位置xxx上水体的流速,q为t时间流入水体的侧向流,默认右为正方向。

连续性方程

考虑x轴上的一个固定区间上的水体,则对水体的高度h对x进行积分,可以求得这一单位宽度水体的质量为
∫x0x1ρh(x,t)dx\int_{x_0}^{x_1}\rho h(x,t)dx ∫x0​x1​​ρh(x,t)dx水体流动过程中,单位区间[x0,x1][x_0,x_1][x0​,x1​]内的水体质量会发生变化,单位时间变化量可以通过水体质量对ttt的导数求得,即
ddt∫x0x1ρh(x,t)dx\frac {d}{dt} \int_{x_0}^{x_1}\rho h(x,t)dx dtd​∫x0​x1​​ρh(x,t)dx即,
∫x0x1∂h(x,t)∂tdx\int_{x_0}^{x_1} \frac {\partial h(x,t)}{\partial t} dx ∫x0​x1​​∂t∂h(x,t)​dx由于水体流入、流出、以及侧向流q,单位区间[x0,x1][x_0,x_1][x0​,x1​]内的水体质量才发生变化,故水体质量的变化还可以表示为
ρ(v(x0,t)h(x0,t)−v(x1,t)h(x1,t)+q(x1−x0))dt/dt\rho (v(x_0,t)h(x_0,t)-v(x_1,t)h(x_1,t)+q(x_1-x_0))dt/dtρ(v(x0​,t)h(x0​,t)−v(x1​,t)h(x1​,t)+q(x1​−x0​))dt/dt
ρ\rhoρ为常量,故,
∫x0x1∂h(x,t)∂tdx=v(x0,t)h(x0,t)−v(x1,t)h(x1,t)+q(x1−x0)\int_{x_0}^{x_1} \frac {\partial h(x,t)}{\partial t} dx=v(x_0,t)h(x_0,t)-v(x_1,t)h(x_1,t)+q(x_1-x_0)∫x0​x1​​∂t∂h(x,t)​dx=v(x0​,t)h(x0​,t)−v(x1​,t)h(x1​,t)+q(x1​−x0​)
右边使用反向运用牛顿-莱布尼兹公式,并将时间导数移至积分号内部,即
∫x0x1f(x,t)=F(x0,t)−F(x1,t)\int_{x_0}^{x_1}f(x,t) = F(x_0,t)-F(x_1,t) ∫x0​x1​​f(x,t)=F(x0​,t)−F(x1​,t)v(x0,t)h(x0,t)−v(x1,t)h(x1,t)=−∫x0x1∂∂xv(x,t)h(x,t)dx+∫x0x1qdxv(x_0,t)h(x_0,t)-v(x_1,t)h(x_1,t) = -\int_{x_0}^{x_1}\frac{\partial}{\partial x}v(x,t)h(x,t)dx+\int_{x_0}^{x_1}q dxv(x0​,t)h(x0​,t)−v(x1​,t)h(x1​,t)=−∫x0​x1​​∂x∂​v(x,t)h(x,t)dx+∫x0​x1​​qdx故,

∫x0x1∂h(x,t)∂tdx+∫x0x1(∂∂xv(x,t)h(x,t)−q)dx=0\int_{x_0}^{x_1} \frac {\partial h(x,t)}{\partial t} dx+ \int_{x_0}^{x_1}(\frac{\partial}{\partial x}v(x,t)h(x,t)-q)dx=0∫x0​x1​​∂t∂h(x,t)​dx+∫x0​x1​​(∂x∂​v(x,t)h(x,t)−q)dx=0故,
∫x0x1(∂h(x,t)∂t+∂∂xv(x,t)h(x,t)−q)dx=0\int_{x_0}^{x_1} (\frac {\partial h(x,t)}{\partial t} + \frac{\partial}{\partial x}v(x,t)h(x,t)-q)dx=0∫x0​x1​​(∂t∂h(x,t)​+∂x∂​v(x,t)h(x,t)−q)dx=0即,
∫x0x1(δQδx+δAδt−q)dx=0\int_{x_0}^{x_1}(\frac{\delta Q}{\delta x}+\frac{\delta A}{\delta t}-q)dx=0∫x0​x1​​(δxδQ​+δtδA​−q)dx=0故,
δQδx+δAδt=q\frac{\delta Q}{\delta x}+\frac{\delta A}{\delta t}=qδxδQ​+δtδA​=q

动量方程

使用物理守恒定律,即物质守恒定律,就可以得到 dmdt=∂(ρuA)∂t+∂(ρu2A)∂x\frac{dm}{dt} = \frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^{2} A)}{\partial x}dtdm​=∂t∂(ρuA)​+∂x∂(ρu2A)​。

物质守恒定律要求,在任意给定的时间间隔内,物质的数量必须是定值。因此,对于任意的位置,ρuA\rho u AρuA 的变化量必须与该位置的流量相等。流量定义为 ρu2A\rho u^{2} Aρu2A。因此,我们有:

dmdt=∂(ρuA)∂t+∂(ρu2A)∂x\frac{dm}{dt} = \frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^{2} A)}{\partial x}dtdm​=∂t∂(ρuA)​+∂x∂(ρu2A)​

其中 ∂(ρuA)∂t\frac{\partial (\rho u A)}{\partial t}∂t∂(ρuA)​ 表示该位置的物质变化量,∂(ρu2A)∂x\frac{\partial (\rho u^{2} A)}{\partial x}∂x∂(ρu2A)​ 表示该位置的流量。

首先,我们假设流体的速度为 uuu,流体的体积为 AAA,河床的高度为 hhh。

对于流体的任意一小部分,其动量可以表示为:

m=ρuAm = \rho u Am=ρuA

其中 ρ\rhoρ 是流体的密度。

因此,当河流在某一时刻流动,我们可以表示为:

dmdt=∂(ρuA)∂t+∂(ρu2A)∂x\frac{dm}{dt} = \frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^{2} A)}{\partial x}dtdm​=∂t∂(ρuA)​+∂x∂(ρu2A)​

将这个式子代入动量守恒原理,即:

∂(ρuA)∂t+∂(ρu2A)∂x=−∂(ρgAh)∂x−∂F∂x\frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^{2} A)}{\partial x} = -\frac{\partial (\rho g A h)}{\partial x} - \frac{\partial F}{\partial x}∂t∂(ρuA)​+∂x∂(ρu2A)​=−∂x∂(ρgAh)​−∂x∂F​

其中 ggg 是重力加速度,FFF 表示阻力。

将 u2u^{2}u2 改写为 Q2A\frac{Q^{2}}{A}AQ2​,即:

∂(ρuA)∂t+∂(ρQ2AA)∂x=−∂(ρgAh)∂x−∂F∂x\frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho \frac{Q^{2}}{A} A)}{\partial x} = -\frac{\partial (\rho g A h)}{\partial x} - \frac{\partial F}{\partial x}∂t∂(ρuA)​+∂x∂(ρAQ2​A)​=−∂x∂(ρgAh)​−∂x∂F​

化简,得到:

∂Q∂t+∂(ρQ2A)∂x+ρgA∂h∂x+∂F∂x=0\frac{\partial Q}{\partial t} + \frac{\partial (\rho \frac{Q^{2}}{A})}{\partial x} + \rho g A \frac{\partial h}{\partial x} + \frac{\partial F}{\partial x} = 0∂t∂Q​+∂x∂(ρAQ2​)​+ρgA∂x∂h​+∂x∂F​=0

在一般情况下,FFF 可以被模拟为 gQ∣Q∣C2ARg Q \frac{|Q|}{C^{2} A R}gQC2AR∣Q∣​,其中 CCC 是河流的波速,RRR 是河流的半径。

最后,可以得到:
δQδt+δ(αQ2A)δx+gAδhδx+gQ∣Q∣C2AR=0\frac{\delta Q}{\delta t}+\frac{\delta\left(\alpha \frac{Q^{2}}{A}\right)}{\delta x}+g A \frac{\delta h}{\delta x}+\frac{g Q|Q|}{C^{2} A R}=0δtδQ​+δxδ(αAQ2​)​+gAδxδh​+C2ARgQ∣Q∣​=0
动量方程中的第一项 δQδt\frac{\delta Q}{\delta t}δtδQ​ 表示流速 QQQ 关于时间的变化; 第二项δ(αQ2A)δx\frac{\delta\left(\alpha \frac{Q^{2}}{A}\right)}{\delta x}δxδ(αAQ2​)​ 表示流速 QQQ关于空间的变化; 第三项 gAδhδxg A \frac{\delta h}{\delta x}gAδxδh​ 表示水位 hhh 关于空间的变化; 第四项gQ∣Q∣C2AR\frac{g Q|Q|}{C^{2} A R}C2ARgQ∣Q∣​ 表示河流的动力阻力关于流速的影响。

二、方程离散

圣维南方程中的连续性方程和动量方程通过有限差分法进行离散,计算网格由流量点和水位点组成,其中流量点和水位点在同一时间步长下分别进行 计算,见图。计算网格由模型自动生成,水位点是横断面所在的位置,相邻水位点之间的距离可能不同,流量点位于两个相邻的水位点之间。计算网格点的分布遵循以下规则:①河段上下游端点为计算水位点;②支流入流点为计算水位点;③实测断面资料点为计算水位点;④模型根据max△r值自动插入的点为计算水位点;⑤建筑物点为计算水位点;⑥两个水位点之间只存在一个计算流量点。

1. 连续方程

在连续方程中,引入蓄存宽度bsb_{s}bs​:
∂A∂t=bs∂h∂t(1)\frac{\partial A}{\partial t}=b_{s}\frac{\partial h}{\partial t}\tag{1} ∂t∂A​=bs​∂t∂h​(1)
从而连续方程转变为:
∂Q∂x+bs∂h∂t=q(2)\frac{\partial Q}{\partial x}+b_{s}\frac{\partial h}{\partial t}=q\tag{2} ∂x∂Q​+bs​∂t∂h​=q(2)
由公式可以看出仅流量Q与x有关,方程很容易得到以h点为中心的6点隐式格式见图。

应用该离散格式,则连续方程变为:
∂Q∂x≈Qj+1n+1+QJ+1n2−Qj−1n+1+Qj−1n2Δ2xj(3)\frac{\partial Q}{\partial x} \approx \frac{\frac{Q_{j+1}^{n+1}+Q_{J+1}^{n}}{2}-\frac{Q_{j-1}^{n+1}+Q_{j-1}^{n}}{2}}{\Delta 2 x_{j}}\tag{3} ∂x∂Q​≈Δ2xj​2Qj+1n+1​+QJ+1n​​−2Qj−1n+1​+Qj−1n​​​(3)∂h∂t≈hjn+1−hjnΔt(4)\frac{\partial h}{\partial t} \approx \frac{h_{j}^{n+1}-h_{j}^{n}}{\Delta t}\tag{4} ∂t∂h​≈Δthjn+1​−hjn​​(4)bs≈Ae,j+Ae,j+1Δ2xj(5)b_{s} \approx\frac{A_{e, j}+A_{e, j+1}}{\Delta 2 x_{j}}\tag{5} bs​≈Δ2xj​Ae,j​+Ae,j+1​​(5)式中: Ae,jA_{e, j}Ae,j​为网格点 j-1 与 j 之间的水面面积;Ao,j+1A_{o, j+1}Ao,j+1​为网格点 j 与 j+1 之间的水面面积;Δ2xj\Delta 2 x_{j}Δ2xj​ 为网格点 j-1 与 j+1 之间的距离。
将式(3)、式(4)代入方程(2) 变为:
Qj+1n+1+Qj+1n2−Qj−1n+1+Qj−1n2△2xj+bs(hjn+1−hjn)△t=qj\frac {\frac{Q_ {j+1}^ {n+1}+Q_ {j+1}^ {n}}{2}-\frac {Q_ {j-1}^ {n+1}+Q_ {j-1}^ {n}}{2}}{\triangle 2x_{j}}+b_ {s} \frac {(h_ {j}^ {n+1}-h_{j}^ {n})}{\triangle t} =q_{j} △2xj​2Qj+1n+1​+Qj+1n​​−2Qj−1n+1​+Qj−1n​​​+bs​△t(hjn+1​−hjn​)​=qj​
−12△2xjQj−1n+1+bs△thjn+1+12△2xjQj+1n+1+(−12△2xjQj−1n−bs△thjn+12△2xjQj+1n)=qj-\frac {1}{2\triangle 2x_{j}}Q_ {j-1}^ {n+1}+ \frac {b_ {s}}{\triangle t}h_ {j}^ {n+1}+\frac {1}{2\triangle 2x_{j}}Q_ {j+1}^ {n+1} +(-\frac {1}{2\triangle 2x_{j}}Q_ {j-1}^ {n}-\frac {b_ {s}}{\triangle t}h_ {j}^ {n}+\frac {1}{2\triangle 2x_{j}}Q_ {j+1}^ {n} )=q_{j} −2△2xj​1​Qj−1n+1​+△tbs​​hjn+1​+2△2xj​1​Qj+1n+1​+(−2△2xj​1​Qj−1n​−△tbs​​hjn​+2△2xj​1​Qj+1n​)=qj​
−12△2xjQj−1n+1+bs△thjn+1+12△2xjQj+1n+1=qj+12△2xjQj−1n+bs△thjn−12△2xjQj+1n-\frac {1}{2\triangle 2x_{j}}Q_ {j-1}^ {n+1}+ \frac {b_ {s}}{\triangle t}h_ {j}^ {n+1}+\frac {1}{2\triangle 2x_{j}}Q_ {j+1}^ {n+1} =q_{j}+\frac {1}{2\triangle 2x_{j}}Q_ {j-1}^ {n}+\frac {b_ {s}}{\triangle t}h_ {j}^ {n}-\frac {1}{2\triangle 2x_{j}}Q_ {j+1}^ {n} −2△2xj​1​Qj−1n+1​+△tbs​​hjn+1​+2△2xj​1​Qj+1n+1​=qj​+2△2xj​1​Qj−1n​+△tbs​​hjn​−2△2xj​1​Qj+1n​
该方程可以简化为:
αjQj−1n+1+βjhjn+1+γjQj+1n+1=qj−αjQj−1n+βjhjn−γjQj+1n\alpha_{j}Q_{j-1}^{n+1}+\beta_{j}h_{j}^{n+1}+\gamma_{j}Q_{j+1}^{n+1}=q_{j}-\alpha_{j}Q_{j-1}^{n}+\beta_{j}h_{j}^{n}-\gamma_{j}Q_{j+1}^{n} αj​Qj−1n+1​+βj​hjn+1​+γj​Qj+1n+1​=qj​−αj​Qj−1n​+βj​hjn​−γj​Qj+1n​
αjQj−1n+1+βjhjn+1+γjQj+1n+1=δj(6)\alpha_{j}Q_{j-1}^{n+1}+\beta_{j}h_{j}^{n+1}+\gamma_{j}Q_{j+1}^{n+1}=\delta_{j}\tag{6} αj​Qj−1n+1​+βj​hjn+1​+γj​Qj+1n+1​=δj​(6)式中:a、β、γ为b和δ的函数,其值决定于h点在时间n处及Q点在时间n+1/2处的值。

注:这一步简化还没理解,参考chatGPT回答,揣测了一下
其中,
αj=−12△2xj\alpha_{j} = -\frac {1}{2\triangle 2x_{j}}αj​=−2△2xj​1​βj=bs△t\beta_{j} = \frac {b_ {s}}{\triangle t}βj​=△tbs​​γj=12△2xj\gamma_{j} =\frac {1}{2\triangle 2x_{j}}γj​=2△2xj​1​δj=qj−αjQj−1n+βjhjn−γjQj+1n\delta_{j} =q_{j}-\alpha_{j}Q_{j-1}^{n}+\beta_{j}h_{j}^{n}-\gamma_{j}Q_{j+1}^{n}δj​=qj​−αj​Qj−1n​+βj​hjn​−γj​Qj+1n​

2. 动量方程

动量方程集中在流量点,其网格形式为以Q点为中心点的差分格式见图3。

依据6点中心Abbott - Ionescu差分法,动量方程可以表示如下:
∂Q∂t≈Qjn+1−Qjn△t(7)\frac{\partial Q}{\partial t}\approx \frac{Q_{j}^{n+1}-Q_{j}^{n}}{\triangle t}\tag{7} ∂t∂Q​≈△tQjn+1​−Qjn​​(7)∂(αQ2A)∂x≈[αQ2A]j+1n+1/2−[αQ2A]j−1n+1/2△2xi(8)\frac {\partial (\alpha \frac {Q^2}{A})}{\partial x}\approx \frac{[\alpha \frac{Q^2}{A}]_{j+1}^{n+1/2}-[\alpha \frac{Q^2}{A}]_{j-1}^{n+1/2}}{\triangle 2xi}\tag{8} ∂x∂(αAQ2​)​≈△2xi[αAQ2​]j+1n+1/2​−[αAQ2​]j−1n+1/2​​(8)∂h∂x=hj+1n+1+hj+1n2−hj−1n+1+hj−1n2△2xj(9)\frac {\partial h}{\partial x}=\frac{\frac{h_{j+1}^{n+1}+h_{j+1}^{n}}{2}-\frac{h_{j-1}^{n+1}+h_{j-1}^{n}}{2}}{\triangle 2x_{j}}\tag{9} ∂x∂h​=△2xj​2hj+1n+1​+hj+1n​​−2hj−1n+1​+hj−1n​​​(9)
对公式(8)中的二次项,引入以下公式:
Q2≈θQjn+1Qjn−(θ−1)QjnQjn(10)Q^2\approx \theta Q_{j}^{n+1}Q_{j}^{n}-(\theta -1)Q_{j}^{n}Q_{j}^{n}\tag{10} Q2≈θQjn+1​Qjn​−(θ−1)Qjn​Qjn​(10)式 中 : θ 角的值通过HD参数文件“默认值”中的“THETA”系数来给定,默认值为1。
由上,动量方程可以表达为:
αjhj−1n+1+βjQjn+1+γjhj+1n+1=δj(11)\alpha _{j}h_{j-1}^{n+1}+\beta _{j}Q_{j}^{n+1}+\gamma _{j}h_{j+1}^{n+1}=\delta _{j}\tag{11} αj​hj−1n+1​+βj​Qjn+1​+γj​hj+1n+1​=δj​(11)其中αj=f(A)\alpha _j=f(A)αj​=f(A)
βj=f(Qjn,Δt,Δx,C,A,R)\beta _j=f(Q_{j}^{n},\Delta t,\Delta x,C,A,R)βj​=f(Qjn​,Δt,Δx,C,A,R)γj=f(A)\gamma _j=f(A)γj​=f(A)δj=f(A,Δx,Δt,α,q,v,θ,hj−1n,Qj−1n+1/2,Qjn,hj+1n,Qj+1n+1/2)\delta _j=f(A,\Delta x,\Delta t,\alpha,q,v,\theta,h _{j-1}^{n},Q _{j-1}^{n+1/2},Q _{j}^{n},h _{j+1}^{n},Q _{j+1}^{n+1/2})δj​=f(A,Δx,Δt,α,q,v,θ,hj−1n​,Qj−1n+1/2​,Qjn​,hj+1n​,Qj+1n+1/2​)
在默认的条件下,软件在一个时间步长里用两次迭代来对这些方程进行求解。初次迭代起始于第 一个时间步长,第二次迭代采用第一次计算值的中心差值来进行计算。迭代次数可以通过NoITER系数来进行修改。

三. 离散方程求解

方程求解采用追赶法(双扫描法)。

1.河道方程

根据水流连续方程和运动方程整理后的最终形式[式(6)和式(11)],河道内某一 水力变量Z,(水位h,或流量Q,)与相邻差分格点的水力参数的关系可以表示为一线性方程:
αjZj−1n+1+βjZjn+1+γjZj+1n+1=δj(12)\alpha _{j}Z _{j-1}^{n+1}+\beta _{j}Z _{j}^{n+1}+\gamma _{j}Z _{j+1}^{n+1}=\delta _{j}\tag{12}αj​Zj−1n+1​+βj​Zjn+1​+γj​Zj+1n+1​=δj​(12)
假设某河道由n个差分格点离散化,由于河道首末计算断面均为水位点,所以n为奇数。对于河道的所有差分格点写出式(12),可以得到n个线性方程:
α1Husn+1+β1h1n+1+γ1Q2n+1=δ1\alpha _{1}H _{us}^{n+1}+\beta _{1}h_{1}^{n+1}+\gamma _{1}Q _{2}^{n+1}=\delta _{1} α1​Husn+1​+β1​h1n+1​+γ1​Q2n+1​=δ1​α2h1n+1+β2Q2n+1+γ2h3n+1=δ2\alpha _{2}h _{1}^{n+1}+\beta _{2}Q_{2}^{n+1}+\gamma _{2}h _{3}^{n+1}=\delta _{2} α2​h1n+1​+β2​Q2n+1​+γ2​h3n+1​=δ2​α3Q2n+1+β3h3n+1+γ3Q4n+1=δ3\alpha _{3}Q _{2}^{n+1}+\beta _{3}h_{3}^{n+1}+\gamma _{3}Q _{4}^{n+1}=\delta _{3} α3​Q2n+1​+β3​h3n+1​+γ3​Q4n+1​=δ3​α4h3n+1+β4Q4n+1+γ4h5n+1=δ4\alpha _{4}h _{3}^{n+1}+\beta _{4}Q_{4}^{n+1}+\gamma _{4}h _{5}^{n+1}=\delta _{4} α4​h3n+1​+β4​Q4n+1​+γ4​h5n+1​=δ4​α5Q4n+1+β5h5n+1+γ5Q6n+1=δ5(13)\alpha _{5}Q _{4}^{n+1}+\beta _{5}h_{5}^{n+1}+\gamma _{5}Q _{6}^{n+1}=\delta _{5}\tag{13} α5​Q4n+1​+β5​h5n+1​+γ5​Q6n+1​=δ5​(13)...... ...αn−2Qn−3n+1+βn−2hn−2n+1+γn−2Qn−1n+1=δn−2\alpha _{n-2}Q _{n-3}^{n+1}+\beta _{n-2}h_{n-2}^{n+1}+\gamma _{n-2}Q _{n-1}^{n+1}=\delta _{n-2} αn−2​Qn−3n+1​+βn−2​hn−2n+1​+γn−2​Qn−1n+1​=δn−2​αn−1hn−2n+1+βn−1Qn−1n+1+γn−1hnn+1=δn−1\alpha _{n-1}h _{n-2}^{n+1}+\beta _{n-1}Q_{n-1}^{n+1}+\gamma _{n-1}h _{n}^{n+1}=\delta _{n-1} αn−1​hn−2n+1​+βn−1​Qn−1n+1​+γn−1​hnn+1​=δn−1​αnQn−1n+1+βnhnn+1+γnHdsn+1=δn\alpha _{n}Q _{n-1}^{n+1}+\beta _{n}h_{n}^{n+1}+\gamma _{n}H _{ds}^{n+1}=\delta _{n} αn​Qn−1n+1​+βn​hnn+1​+γn​Hdsn+1​=δn​
这样方程就多出2个未知量。第一个方程中的HusH_{us}Hus​和最后一个方程中的HdsH_{ds}Hds​分别代表上下游节点的水位。某一河道第一个网格点的水位等于与之相连河段上游节点的水位:h1=Hush_{1}=H_{us}h1​=Hus​,即:α1=−1,β1=1,γ1=0,δ1=0\alpha _{1}=-1, \beta _{1}=1,\gamma _{1}=0,\delta _{1}=0α1​=−1,β1​=1,γ1​=0,δ1​=0。同样,hn=Hdsh_{n}=H_{ds}hn​=Hds​, 即:αn=0,βn=1,γn=−1,δ1=0\alpha _{n}=0, \beta _{n}=1,\gamma _{n}=-1,\delta _{1}=0αn​=0,βn​=1,γn​=−1,δ1​=0。
对于单一河道,只要给出上下游水位边界,即HusH_{us}Hus​与HdsH_{ds}Hds​为已知,就可用消元法求解方程组 (13)了。
对于河网,方程组(13)用矩阵的形式表示为:
[α1β1γ1δ1α2β2γ2δ2α3β3γ3δ3α4β4γ4δ4α5β5γ5δ5......αn−2βn−2γn−2δn−2αn−1βn−1γn−1δn−1αnβnγnδn]\begin{bmatrix} \alpha _{1}&\beta _{1}&\gamma _{1}& & & & & & & & & &\delta _{1} \\ & \alpha _{2}&\beta _{2}&\gamma _{2}& & & & & & & & &\delta _{2} \\ & & \alpha _{3}&\beta _{3}&\gamma _{3}& & & & & & & &\delta _{3} \\ & & & \alpha _{4}&\beta _{4}&\gamma _{4}& & & & & & &\delta _{4} \\ & & & & \alpha _{5}&\beta _{5}&\gamma _{5}& & & & & &\delta _{5}\\ & & & & & .& .& .& & && \\ & & & & & & .& . & . && & \\ & & & & & & & \alpha _{n-2}&\beta _{n-2}&\gamma _{n-2}& & &\delta _{n-2}\\ & & & & & & & & \alpha _{n-1}&\beta _{n-1}&\gamma _{n-1}& &\delta _{n-1}\\ & & & & & & & & & \alpha _{n}&\beta _{n}&\gamma _{n}&\delta _{n} \end{bmatrix} ​α1​​β1​α2​​γ1​β2​α3​​γ2​β3​α4​​γ3​β4​α5​​γ4​β5​.​γ5​..​..αn−2​​.βn−2​αn−1​​γn−2​βn−1​αn​​γn−1​βn​​γn​​δ1​δ2​δ3​δ4​δ5​δn−2​δn−1​δn​​​
用消元法变成了
[a11b1c1a21b2c2a31b3c3a41b4c4a51b5c5........an−21bn−2cn−2an−11bn−1cn−1an1bncn]\begin{bmatrix} a _{1}&1& & & & & & & & & &b _{1}&c _{1}\\ a _{2}& &1& & & & & & & & &b _{2}&c _{2}\\ a _{3}& & &1& & & & & & & &b _{3}&c _{3}\\ a _{4}& & & &1& & & & & & &b _{4}&c _{4}\\ a _{5}& & & & &1& & & & & &b _{5}&c _{5}\\ .& & & & & &.& & & & &.&.\\ .& & & & & & &.& & & &.&.\\ a _{n-2}& & & & & & & &1& & &b _{n-2}&c _{n-2}\\ a _{n-1}& & & & & & & & &1& &b _{n-1}&c _{n-1}\\ a _{n}& & & & & & & & & &1&b _{n}&c _{n} \end{bmatrix} ​a1​a2​a3​a4​a5​..an−2​an−1​an​​1​1​1​1​1​.​.​1​1​1​b1​b2​b3​b4​b5​..bn−2​bn−1​bn​​c1​c2​c3​c4​c5​..cn−2​cn−1​cn​​​
由此便能将河道内任意点的水力变量ZjZ_{j}Zj​(水位或流量)表示为上下游节点的水位函数:
Zjn+1=cj−ajHusn+1−bjHdsn+1(14)Z _{j}^{n+1} = c_{j}-a_{j}H_{us}^{n+1}-b_{j}H_{ds}^{n+1}\tag{14} Zjn+1​=cj​−aj​Husn+1​−bj​Hdsn+1​(14)在求解过程中首先求出河道中各节点的水位值,之后便可应用式(14)求解任一河段,任一差分格点的水力参数。

2.节点方程


如图所示,绕节点的控制体连续性方程为:
Hn+1−HnΔtAfl=12(QA,n−1n+QB,n−1n−QC,2n)+12(QA,n−1n+1+QB,n−1n+1−QC,2n+1)\frac{H^{n+1}-H^{n}}{\Delta t} A_{fl}=\frac{1}{2}\left(Q_{A, n-1}^{n}+Q_{B, n-1}^{n}-Q_{C, 2}^{n}\right)+\frac{1}{2}\left(Q_{A, n-1}^{n+1}+Q_{B, n-1}^{n+1}-Q_{C, 2}^{n+1}\right)ΔtHn+1−Hn​Afl​=21​(QA,n−1n​+QB,n−1n​−QC,2n​)+21​(QA,n−1n+1​+QB,n−1n+1​−QC,2n+1​)
将上述方程中右边第2式的3项分别以式(14)替代,可以得到:
Hn+1−HnΔtAf=12(QAn+QBn−QCn)+12(cA,n−1−aA,n−1HA,usn+1−bA,n−1Hn+1)+cB,n−1−aB,n−1HB,usn+1−bB,n−1Hn+1+cC,2−aC,2Hn+1−bC,2HC,dsn+1(15)\frac{H^{n+1}-H^{n}}{\Delta t} A_{f}=\frac{1}{2}\left(Q_{A}^{n}+Q_{B}^{n}-Q_{C}^{n}\right)+\frac{1}{2}\left(c_{A, n-1}-a_{A, n-1} H_{A, us}^{n+1}-b_{A, n-1} H^{n+1}\right)+\\ c_{B, n-1} -a_{B, n-1} H_{B, us}^{n+1}-b_{B, n-1} H^{n+1}+c_{C, 2}-a_{C, 2} H^{n+1}-b_{C, 2} H_{C, ds}^{n+1} \tag{15}ΔtHn+1−Hn​Af​=21​(QAn​+QBn​−QCn​)+21​(cA,n−1​−aA,n−1​HA,usn+1​−bA,n−1​Hn+1)+cB,n−1​−aB,n−1​HB,usn+1​−bB,n−1​Hn+1+cC,2​−aC,2​Hn+1−bC,2​HC,dsn+1​(15)

式中:HHH为该节点的水位,HA,usH_{A,us}HA,us​、HB,usH_{B,us}HB,us​分别为支流A,B上游端节点水位,HC,dsH_{C,ds}HC,ds​为支流C下游端节 点水位。
在式(15)中,将某个节点的水位表示为与之直接相连的河道节点水位的线性函数。同样,对于河网所有节点(假设为N个),可以得到N个类似的方程(节点方程组)。在边界水位或流量为已知的情况下,可以利用高斯消元法直接求解节点方程组,得到各个节点的水位,进而回代式(15)求解任意河道任意网格点的水位或流量。
原则上节点可以任意编码,但对于大型复杂河网,这样得到的节点方程组的系数矩阵将是一个阶数很高的稀疏矩阵,存贮量大,运算非常耗时。大型稀疏矩阵求解计算时间主要取决于矩阵主对角线非零元素的宽度。可以通过对河网节点进行优化编码的方法来降低节点方程组系数矩阵的带宽,使之成为主对角线元素占优的矩阵,从而方便方程组的求解,并大大减少计算耗时。

3.边界条件

若在河道边界节点上给出水位的时间变化过程:h=h(t)h=h(t)h=h(t)。此时,边界上的节点方程为(假设边界所在河道编号为jjj):
hj,1n+1=Husn+1,或hj,nn+1=Hdsn+1h_{j,1}^{n+1}=H_{us}^{n+1},或h_{j,n}^{n+1}=H_{ds}^{n+1}hj,1n+1​=Husn+1​,或hj,nn+1​=Hdsn+1​
若在河道边界节点上给出流量的时间变化过程:Q=Q(t)Q=Q(t)Q=Q(t)。对控制体应用连续方程可以得到:
Hn+1−HnΔtAfl=12(Qbn−Q2n)+12(Qbn+1−Q2n+1)\frac {H^{n+1}-H^{n}}{\Delta t}{A_{fl}} = \frac {1}{2}(Q _{b}^{n}-Q_{2}^{n})+ \frac {1}{2}(Q _{b}^{n+1}-Q_{2}^{n+1}) ΔtHn+1−Hn​Afl​=21​(Qbn​−Q2n​)+21​(Qbn+1​−Q2n+1​)根据(14)将上式中的Q2n+1Q_{2}^{n+1}Q2n+1​代入,可以得到
Hn+1−HnΔtAfl=12(Qbn−Q2n)+12(Qbn+1−c2+a2Hn+1+b2Hdsn+1)(16)\frac {H^{n+1}-H^{n}}{\Delta t}{A_{fl}} = \frac {1}{2}(Q _{b}^{n}-Q_{2}^{n})+ \frac {1}{2}(Q _{b}^{n+1}-c_{2}+a_{2}H^{n+1}+b_{2}H_{ds}^{n+1})\tag{16} ΔtHn+1−Hn​Afl​=21​(Qbn​−Q2n​)+21​(Qbn+1​−c2​+a2​Hn+1+b2​Hdsn+1​)(16)


若在河道边界节点上给出流量的时间变化过程:Q=Q(h)Q=Q(h)Q=Q(h),其处理方法同流量边界,得到与式子(16)类似的方程,只是方程中的QbnQ_{b}^{n}Qbn​和Qbn+1Q_{b}^{n+1}Qbn+1​由流量水位关系得到。

参考文献

  • 《DHI MIKE FLOOD 洪水模拟技术应用与研究》,衣秀勇等编著
  • https://zhuanlan.zhihu.com/p/372837854
  • chatGPT

【MIKE水动力】MIKE11基本原理相关推荐

  1. MIKE水动力笔记13_数字化海图2之克里金插值

    本文目录 前言 Step 1 调出地统计分析工具 Step 2 克里金插值设置 Step 3 调整图幅范围及裁剪 Step 4 转为栅格文件并保存 前言 在进行MIKE水动力建模之初,需要准备好水深数 ...

  2. MIKE水动力笔记7_实测数据与模型输出结果的拟合对比

    本文目录 前言 Step 1 拟合对比前的准备工作 Step 2 从模型输出结果dfsu文件提取出站位点处的模拟潮位dfs0文件 Step 3 将两个dfs0文件插进绘图板 Step 4 对图面进行必 ...

  3. MIKE水动力笔记9_大潮小潮对应的涨急落急时刻流场图

    本文目录 前言 Step 1 确定大潮日.小潮日.涨急时刻.落急时刻 Step 2 计算四个时刻对应的时间步数并得到对应流场图 Step 3 调整图面 后记 前言 在水动力模型运行完.模拟出结果之后, ...

  4. MIKE水动力笔记2_水动力基础理论知识

    本文目录 前言 [第2章 水动力学]重点知识 2.1 水动力过程 2.1.1 水的密度 2.1.2 守恒律 2.1.3 对流和扩散 2.1.4 质量守恒方程 2.1.5 大气驱动力 2.1.6 科氏力 ...

  5. MIKE 21 教程 2.1水动力模型介绍

    前面的第一章节相关博文中,我们讲解了MIKE21入门操作与网格文件的制作,接下来进入第二章节,水动力模型. 1 水动力模型简介 水动力模型是MIKE一切模拟的基础,用来模拟水的流速流量,高程方向等物理 ...

  6. MIKE 21 教程 2.3 水动力模块教学:求解方程与参数设置(Solution Technique),水深校正设置(Depth Correction)

    上一节讲解了水动力模块的基本设置:网格导入,时间步长设置与模块选择 MIKE 21 教程 2.2 Domain, Time, Module Selection设置教学 目录 1 Solution Te ...

  7. 【MIKE HYDRO】某河道MIKE HYDRO水质水动力模拟项目-水动力模型构建

    水动力模型构建 1.导入地图 由于只提供所属区域对应坐标系和左下角.右上角坐标,所以: 选择Map configurations->Coordinate system,Map view coor ...

  8. 【HEC-RAS水动力】HEC-RAS 1D基本原理(恒定流及非恒定流)

    一.数据说明 HEC-RAS模型主要由工程文件 (.prj) 文 件 . 河道地形数据文件 ( .g01).运行文件(p01).非恒定流文件 ( .u01) 等部分组成. 1. 一般数据 在创建并保存 ...

  9. MIKE 21 教程 2.4 水动力模块教学:干湿边界(Flood and Dry),密度关系(Density)

    目录 1. 干湿边界(Flood and Dry) 1.1 含义 1.2 数值设置 2 Density(密度) Barotropic:密度不变.

最新文章

  1. 标签的for循环和if_SO面试题08:如何从一个多层嵌套循环中直接跳出?
  2. linux read while 变量运算
  3. 大数据统计分析毕业设计_基于大数据分析的电子信息类专业毕业设计成绩影响因素研究...
  4. Java类类getClassLoader()方法及示例
  5. java调用 restapi 乱码_Java HttpURLConnection模拟请求Rest接口解决中文乱码问题
  6. python编写函数showmsg(n、name)_Python语言答案
  7. php商品秒杀时间代码,Thinkphp5+Redis实现商品秒杀代码实例讲解
  8. 32线性空间06——行空间和左零空间
  9. rk399_android7.1的mipi驱动代码追踪(部分)
  10. VMware Workstation环境下的Linux网络设置
  11. 2019.03.04【ZJOI2018】【BZOJ5212】【洛谷P4338】历史(假LCT)
  12. 微信小程序:setData 数据传输长度为 1678 KB,存在有性能问题!
  13. 燕大计算机研究生毕业待遇,研究生人均“月薪上万”是真是假,过来人坦言:想想就好,别认真...
  14. YTU OJ 2458: 换啤酒
  15. Springboot 项目学习
  16. 优秀成绩标记—— 小王是班级干部,对于即将到来的三好学生评选,负责统计平均成绩超过85分的同学
  17. 苹果屏幕自动变暗_苹果iOS 14震撼发布 全新功能对标安卓
  18. python有哪些模块安全方向_Python 常用模块
  19. Jenkins 流水线自动化部署 Go 项目
  20. Android TimeoutException治理

热门文章

  1. 第三章-集合论 3.2-Russell 悖论(选读)
  2. 外部css样式不生效的原因
  3. poj 1830 开关问题
  4. esp8266微信wifi配置AIRKISS
  5. 电脑内存条的选配与安装详述
  6. ESP32超详细学习记录:wifi连接最基础方法
  7. 生成微信小程序二维码,可跳转到小程序指定页面。
  8. 通俗易懂了解50个IT专业术语
  9. 新手指南: Linux 新手应该知道的 26 个命令
  10. android 序列化存储对象,android中对象序列化存储