圣维南方程组求解 python 二维圣维南方程

目录

• 说明
• 一、圣维南方程组

• 1、基本要素与假设条件
• 2、圣维南方程组
• 3、公式推导

• 连续性方程
• 动量方程

• 二、方程离散

• 1. 连续方程
• 2. 动量方程

• 三. 离散方程求解

• 1.河道方程
• 2.节点方程
• 3.边界条件

• 参考文献

说明

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

一、圣维南方程组

1、基本要素与假设条件

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

2、圣维南方程组

$\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.$
式中：Q为流量，m³/s；q为侧向入流，m³/s；A为过水面积，m²；h为水位，m；R为水力半径，m；C为谢才系数；a为动量修正系数。

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

3、公式推导

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

连续性方程

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

$\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$故，
$\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$即，
$\int_{x_0}^{x_1}(\frac{\delta Q}{\delta x}+\frac{\delta A}{\delta t}-q)dx=0$故，
$\frac{\delta Q}{\delta x}+\frac{\delta A}{\delta t}=q$

动量方程

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

$\frac{dm}{dt} = \frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^{2} A)}{\partial x}$

其中 $\frac{\partial (\rho u A)}{\partial t}$ 表示该位置的物质变化量，$\frac{\partial (\rho u^{2} A)}{\partial x}$

首先，我们假设流体的速度为 $u$，流体的体积为 $A$，河床的高度为 $h$。

对于流体的任意一小部分，其动量可以表示为：$m = \rho u A$，其中 $\rho$

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

$\frac{dm}{dt} = \frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho u^{2} A)}{\partial 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}$

其中 $g$ 是重力加速度，$F$

将 $u^{2}$ 改写为 $\frac{Q^{2}}{A^2}$，即：

$\frac{\partial (\rho u A)}{\partial t} + \frac{\partial (\rho \frac{Q^{2}}{A^2} A)}{\partial x} = -\frac{\partial (\rho g A h)}{\partial x} - \frac{\partial F}{\partial x}$

化简，得到：

$\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$

在一般情况下，$F$ 可以被模拟为 $g Q \frac{|Q|}{C^{2} A R}$，其中 $C$ 是河流的波速，$R$

最后，可以得到：
$\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$
动量方程中的第一项 $\frac{\delta Q}{\delta t}$ 表示流速 $Q$ 关于时间的变化； 第二项$\frac{\delta\left(\alpha \frac{Q^{2}}{A}\right)}{\delta x}$ 表示流速 $Q$关于空间的变化； 第三项 $g A \frac{\delta h}{\delta x}$ 表示水位 $h$ 关于空间的变化； 第四项$\frac{g Q|Q|}{C^{2} A R}$

二、方程离散

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

1. 连续方程

在连续方程中，引入蓄存宽度$b_{s}$：

$\frac{\partial A}{\partial t}=b_{s}\frac{\partial h}{\partial t}\tag{1}$

从而连续方程转变为：

$\frac{\partial Q}{\partial x}+b_{s}\frac{\partial h}{\partial t}=q\tag{2}$

由公式可以看出仅流量Q与x有关，方程很容易得到以h点为中心的6点隐式格式见图。

应用该离散格式，则连续方程变为：

$\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}$$\frac{\partial h}{\partial t} \approx \frac{h_{j}^{n+1}-h_{j}^{n}}{\Delta t}\tag{4}$$b_{s} \approx\frac{A_{e, j}+A_{e, j+1}}{\Delta 2 x_{j}}\tag{5}$式中: $A_{e, j}$为网格点 j-1 与 j 之间的水面面积；$A_{o, j+1}$为网格点 j 与 j+1 之间的水面面积；$\Delta 2 x_{j}$ 为网格点 j-1 与 j+1 之间的距离。

将式(3)、式(4)代入方程(2) 变为：

$\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}$

$-\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}$

$-\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}$

该方程可以简化为：

$\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}$

$\alpha_{j}Q_{j-1}^{n+1}+\beta_{j}h_{j}^{n+1}+\gamma_{j}Q_{j+1}^{n+1}=\delta_{j}\tag{6}$式中：a、β、γ为b和δ的函数，其值决定于h点在时间n处及Q点在时间n+1/2处的值。

注：这一步简化还没理解，参考chatGPT回答，揣测了一下
其中，
$\alpha_{j} = -\frac {1}{2\triangle 2x_{j}}$$\beta_{j} = \frac {b_ {s}}{\triangle t}$$\gamma_{j} =\frac {1}{2\triangle 2x_{j}}$$\delta_{j} =q_{j}-\alpha_{j}Q_{j-1}^{n}+\beta_{j}h_{j}^{n}-\gamma_{j}Q_{j+1}^{n}$

2. 动量方程

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

依据6点中心Abbott - Ionescu差分法，动量方程可以表示如下：

$\frac{\partial Q}{\partial t}\approx \frac{Q_{j}^{n+1}-Q_{j}^{n}}{\triangle t}\tag{7}$$\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}$$\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}$

对公式(8)中的二次项，引入以下公式：

$Q^2\approx \theta Q_{j}^{n+1}Q_{j}^{n}-(\theta -1)Q_{j}^{n}Q_{j}^{n}\tag{10}$式 中 ： θ 角的值通过HD参数文件“默认值”中的“THETA”系数来给定，默认值为1。

由上，动量方程可以表达为：

$\alpha _{j}h_{j-1}^{n+1}+\beta _{j}Q_{j}^{n+1}+\gamma _{j}h_{j+1}^{n+1}=\delta _{j}\tag{11}$其中$\alpha _j=f(A)$

$\beta _j=f(Q_{j}^{n},\Delta t,\Delta x,C,A,R)$$\gamma _j=f(A)$$\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})$

在默认的条件下，软件在一个时间步长里用两次迭代来对这些方程进行求解。初次迭代起始于第 一个时间步长，第二次迭代采用第一次计算值的中心差值来进行计算。迭代次数可以通过NoITER系数来进行修改。

三. 离散方程求解

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

1.河道方程

根据水流连续方程和运动方程整理后的最终形式[式(6)和式(11)],河道内某一 水力变量Z,(水位h,或流量Q,)与相邻差分格点的水力参数的关系可以表示为一线性方程：
$\alpha _{j}Z _{j-1}^{n+1}+\beta _{j}Z _{j}^{n+1}+\gamma _{j}Z _{j+1}^{n+1}=\delta _{j}\tag{12}$
假设某河道由n个差分格点离散化，由于河道首末计算断面均为水位点，所以n为奇数。对于河道的所有差分格点写出式(12),可以得到n个线性方程：
$\alpha _{1}H _{us}^{n+1}+\beta _{1}h_{1}^{n+1}+\gamma _{1}Q _{2}^{n+1}=\delta _{1}$$\alpha _{2}h _{1}^{n+1}+\beta _{2}Q_{2}^{n+1}+\gamma _{2}h _{3}^{n+1}=\delta _{2}$$\alpha _{3}Q _{2}^{n+1}+\beta _{3}h_{3}^{n+1}+\gamma _{3}Q _{4}^{n+1}=\delta _{3}$$\alpha _{4}h _{3}^{n+1}+\beta _{4}Q_{4}^{n+1}+\gamma _{4}h _{5}^{n+1}=\delta _{4}$$\alpha _{5}Q _{4}^{n+1}+\beta _{5}h_{5}^{n+1}+\gamma _{5}Q _{6}^{n+1}=\delta _{5}\tag{13}$$...$$\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}$$\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}$$\alpha _{n}Q _{n-1}^{n+1}+\beta _{n}h_{n}^{n+1}+\gamma _{n}H _{ds}^{n+1}=\delta _{n}$
这样方程就多出2个未知量。第一个方程中的$H_{us}$和最后一个方程中的$H_{ds}$分别代表上下游节点的水位。某一河道第一个网格点的水位等于与之相连河段上游节点的水位：$h_{1}=H_{us}$，即：$\alpha _{1}=-1, \beta _{1}=1,\gamma _{1}=0,\delta _{1}=0$。同样，$h_{n}=H_{ds}$, 即：$\alpha _{n}=0, \beta _{n}=1,\gamma _{n}=-1,\delta _{1}=0$。
对于单一河道，只要给出上下游水位边界，即$H_{us}$与$H_{ds}$为已知，就可用消元法求解方程组 (13)了。
对于河网，方程组（13）用矩阵的形式表示为：
$\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}$
用消元法变成了
$\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}$
由此便能将河道内任意点的水力变量$Z_{j}$(水位或流量)表示为上下游节点的水位函数：
$Z _{j}^{n+1} = c_{j}-a_{j}H_{us}^{n+1}-b_{j}H_{ds}^{n+1}\tag{14}$在求解过程中首先求出河道中各节点的水位值，之后便可应用式(14)求解任一河段，任一差分格点的水力参数。

2.节点方程

如图所示，绕节点的控制体连续性方程为：

$\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)$

将上述方程中右边第2式的3项分别以式(14)替代，可以得到：

$\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}$

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

3.边界条件

若在河道边界节点上给出水位的时间变化过程：$h=h(t)$。此时，边界上的节点方程为(假设边界所在河道编号为$j$)：
$h_{j,1}^{n+1}=H_{us}^{n+1}，或h_{j,n}^{n+1}=H_{ds}^{n+1}$
若在河道边界节点上给出流量的时间变化过程：$Q=Q(t)$。对控制体应用连续方程可以得到：
$\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})$根据（14）将上式中的$Q_{2}^{n+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}-c_{2}+a_{2}H^{n+1}+b_{2}H_{ds}^{n+1})\tag{16}$

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