建模背景
记录原因
本文记录笔者在学习计算水力学课程时,对一维矩形河道建模及编程模拟的过程。一方面方便以后自己回忆学习内容,另一方面分享自己的学习成果,共同进步。
模拟工况
单一河道,断面都为规则的矩形,坡度为0,长20km,糙率n=0.025,初始水深为4m,初始流量为0,断面间距500m,时间步长300s,模拟总时间24h,重力加速度g取9.81m/s2。
边界条件
上边界条件:水深在24h内由4m逐步增加到6m;
下边界条件:水深恒定为4m。
无旁侧如流。
模型结构
圣维南方程组
连续方程: ∂ A ∂ t + ∂ Q ∂ x = 0 \frac{\partial{A}}{\partial{t}}+\frac{\partial{Q}}{\partial{x}}=0 ∂t∂A+∂x∂Q=0
动量方程: ∂ Q ∂ t + ∂ Q u ∂ x + g A ∂ Z ∂ x + g A S f = 0 \frac{\partial{Q}}{\partial{t}}+\frac{\partial{Qu}}{\partial{x}}+gA\frac{\partial{Z}}{\partial{x}}+gAS_f=0 ∂t∂Q+∂x∂Qu+gA∂x∂Z+gASf=0
式中: A A A为断面面积; Q Q Q为通过断面的流量; u u u为断面流速; g g g为重力加速度;
Z Z Z为水位; S f S_f Sf为摩阻比降, S f = n 2 Q ∣ Q ∣ A 2 R 4 3 S_f=\frac{n^2Q|Q|}{A^2R^{\frac{4}{3}}} Sf=A2R34n2Q∣Q∣。
矩形河道简化为
连续方程:
∂
h
∂
t
+
∂
q
∂
x
=
0
\frac{\partial{h}}{\partial{t}}+\frac{\partial{q}}{\partial{x}}=0
∂t∂h+∂x∂q=0
动量方程: ∂ q ∂ t + ∂ ∂ x ( q 2 h ) + g h ∂ Z ∂ x + g n 2 q ∣ q ∣ h 7 3 = 0 \frac{\partial{q}}{\partial{t}}+\frac{\partial}{\partial{x}}\left(\frac{q^2}{h}\right)+gh\frac{\partial{Z}}{\partial{x}}+g\frac{n^2q|q|}{h^{\frac{7}{3}}}=0 ∂t∂q+∂x∂(hq2)+gh∂x∂Z+gh37n2q∣q∣=0
式中: h h h为水深; q q q为单位宽度的流量,简称单宽流量; n n n为河道糙率; g g g为重力加速度; Z Z Z为水位, Z = Z b + h Z=Z_b+h Z=Zb+h, Z b Z_b Zb为河底高程。
Preissmann隐式差分格式
简化四点线性隐格式
f
∣
M
=
f
j
+
1
n
+
f
j
n
2
f\vert_M=\frac{f^n_{j+1}+f^n_j}{2}
f∣M=2fj+1n+fjn
∂ f ∂ x ∣ M = θ ( f j + 1 n + 1 − f j n + 1 Δ x ) + ( 1 − θ ) ( f j + 1 n − f j n Δ x ) \frac{\partial{f}}{\partial{x}}\vert_M=\theta(\frac{f^{n+1}_{j+1}-f^{n+1}_j}{\Delta x}) + (1-\theta)(\frac{f^{n}_{j+1}-f^{n}_j}{\Delta x}) ∂x∂f∣M=θ(Δxfj+1n+1−fjn+1)+(1−θ)(Δxfj+1n−fjn)
∂ f ∂ t ∣ M = f j + 1 n + 1 + f j n + 1 − f j + 1 n − f j n 2 Δ t \frac{\partial{f}}{\partial{t}}\vert_M=\frac{f^{n+1}_{j+1} +f^{n+1}_{j} - f^{n}_{j+1} -f^{n}_j}{2\Delta t} ∂t∂f∣M=2Δtfj+1n+1+fjn+1−fj+1n−fjn
式中: θ \theta θ为加权系数, 0 ≤ θ ≤ 1 0\leq\theta\leq1 0≤θ≤1。 θ = 0 \theta=0 θ=0为显式; 0 < θ < 1 0<\theta<1 0<θ<1为半隐式; θ = 1 \theta=1 θ=1为全隐式。 0.5 ≤ θ ≤ 1 0.5\leq\theta\leq1 0.5≤θ≤1时离散稳定。
离散形式
本文采用全隐式格式,即取
θ
=
1
\theta=1
θ=1。
∂
h
∂
t
≈
Δ
h
Δ
t
=
h
i
+
1
2
t
+
1
−
h
i
+
1
2
t
Δ
t
=
h
i
t
+
1
+
h
i
+
1
t
+
1
2
−
h
i
t
+
h
i
+
1
t
2
Δ
t
\frac{\partial{h}}{\partial{t}} \approx\frac{\Delta h}{\Delta t}=\frac{h_{i+\frac{1}{2}} ^{t+1}-h_{i+\frac{1}{2}} ^{t}}{\Delta t}=\frac{\frac{h_i^{t+1}+h_{i+1}^{t+1}}{2}-\frac{h_i^{t}+h_{i+1}^{t}}{2}}{\Delta t}
∂t∂h≈ΔtΔh=Δthi+21t+1−hi+21t=Δt2hit+1+hi+1t+1−2hit+hi+1t
∂
q
∂
x
≈
Δ
q
Δ
x
=
[
θ
q
i
+
1
t
+
1
+
(
1
−
θ
)
q
i
+
1
t
]
−
[
θ
q
i
t
+
1
+
(
1
−
θ
)
q
i
t
]
Δ
x
=
q
i
+
1
t
+
1
−
q
i
t
+
1
Δ
x
\frac{\partial{q}}{\partial{x}} \approx\frac{\Delta q}{\Delta x}=\frac{\left[\theta q_{i+1}^{t+1}+\left(1-\theta \right)q_{i+1}^t\right] - \left[\theta q_{i}^{t+1}+\left(1-\theta \right)q_{i}^t\right]}{\Delta x}=\frac{q_{i+1}^{t+1}-q_{i}^{t+1}}{\Delta x}
∂x∂q≈ΔxΔq=Δx[θqi+1t+1+(1−θ)qi+1t]−[θqit+1+(1−θ)qit]=Δxqi+1t+1−qit+1
∂
∂
x
(
q
2
h
)
=
∂
∂
x
(
q
⋅
u
)
≈
u
∂
q
∂
x
=
u
i
t
+
u
i
+
1
t
2
∂
q
∂
x
\frac{\partial}{\partial x}\left(\frac{q^2}{h}\right)=\frac{\partial}{\partial x}\left(q\cdot u\right)\approx u\frac{\partial q}{\partial x}=\frac{u_i^t+u_{i+1}^t}{2}\frac{\partial q}{\partial x}
∂x∂(hq2)=∂x∂(q⋅u)≈u∂x∂q=2uit+ui+1t∂x∂q
h
=
h
i
t
+
h
i
+
1
t
2
h=\frac{h_i^t+h_{i+1}^t}{2}
h=2hit+hi+1t
u
=
u
i
t
+
u
i
+
1
t
2
u=\frac{u_i^t+u_{i+1}^t}{2}
u=2uit+ui+1t
q
=
q
i
t
+
q
i
+
1
t
2
q=\frac{q_i^t+q_{i+1}^t}{2}
q=2qit+qi+1t
方程组离散
将离散形式代入简化的圣维南方程组中,可以得到圣维南方程组的Preissmann离散形式
连续方程:
h
i
t
+
1
+
h
i
+
1
t
+
1
−
h
i
t
−
h
i
+
1
t
2
Δ
t
+
q
i
+
1
t
+
1
−
q
i
t
+
1
Δ
x
=
0
\frac{h_i^{t+1}+h_{i+1}^{t+1}-h_i^t-h_{i+1}^t}{2\Delta t }+\frac{q_{i+1}^{t+1}-q_i^{t+1}}{\Delta x}=0
2Δthit+1+hi+1t+1−hit−hi+1t+Δxqi+1t+1−qit+1=0
h
i
t
+
1
−
2
Δ
t
Δ
x
q
i
t
+
1
+
h
i
+
1
t
+
1
+
2
Δ
t
Δ
x
q
i
+
1
t
+
1
=
h
i
t
+
h
i
+
1
t
h_i^{t+1}-\frac{2\Delta t}{\Delta x}q_i^{t+1}+h_{i+1}^{t+1}+\frac{2\Delta t}{\Delta x}q_{i+1}^{t+1}=h_i^t+h_{i+1}^t
hit+1−Δx2Δtqit+1+hi+1t+1+Δx2Δtqi+1t+1=hit+hi+1t
动量方程:
q
i
t
+
1
+
q
i
+
1
t
+
1
−
q
i
t
−
q
i
+
1
t
2
Δ
t
+
u
q
i
+
1
t
+
1
−
q
i
t
+
1
Δ
x
+
g
h
h
i
+
1
t
+
1
−
h
i
t
+
1
Δ
x
−
g
h
S
0
+
q
n
2
∣
q
∣
h
7
3
q
i
t
+
1
+
q
i
+
1
t
+
1
2
=
0
\frac{q_i^{t+1}+q_{i+1}^{t+1}-q_i^t-q_{i+1}^t}{2\Delta t }+u\frac{q_{i+1}^{t+1}-q_i^{t+1}}{\Delta x} + gh\frac{h_{i+1}^{t+1}-h_i^{t+1}}{\Delta x} - ghS_0 +\frac{qn^2|q|}{h^{\frac{7}{3}}}\frac{q_i^{t+1}+q_{i+1}^{t+1}}{2}=0
2Δtqit+1+qi+1t+1−qit−qi+1t+uΔxqi+1t+1−qit+1+ghΔxhi+1t+1−hit+1−ghS0+h37qn2∣q∣2qit+1+qi+1t+1=0
−
g
h
2
Δ
t
Δ
x
h
i
t
+
1
+
(
1
−
u
2
Δ
t
Δ
x
+
g
n
2
∣
q
∣
Δ
t
h
7
3
)
q
i
t
+
1
+
g
h
2
Δ
t
Δ
x
h
i
+
1
t
+
1
+
(
1
+
u
2
Δ
t
Δ
x
+
g
n
2
∣
q
∣
Δ
t
h
7
3
)
q
i
+
1
t
+
1
=
q
i
t
+
q
i
+
1
t
+
2
g
h
Δ
t
S
0
-gh\frac{2 \Delta t}{\Delta x}h_i^{t+1}+(1-u\frac{2 \Delta t}{\Delta x}+\frac{gn^2|q|\Delta t}{h^{\frac{7}{3}}})q_i^{t+1}+gh\frac{2\Delta t}{\Delta x}h_{i+1}^{t+1}+(1+u\frac{2 \Delta t}{\Delta x}+\frac{gn^2|q|\Delta t}{h^{\frac{7}{3}}})q_{i+1}^{t+1}=q_i^t+q_{i+1}^t + 2gh \Delta t S_0
−ghΔx2Δthit+1+(1−uΔx2Δt+h37gn2∣q∣Δt)qit+1+ghΔx2Δthi+1t+1+(1+uΔx2Δt+h37gn2∣q∣Δt)qi+1t+1=qit+qi+1t+2ghΔtS0
设方程组形式如下:
A
h
i
t
+
1
+
B
q
i
t
+
1
+
C
h
i
+
1
t
+
1
+
D
q
i
+
1
t
+
1
=
C
1
Ah_i^{t+1}+Bq_i^{t+1}+Ch_{i+1}^{t+1}+Dq_{i+1}^{t+1}=C_1
Ahit+1+Bqit+1+Chi+1t+1+Dqi+1t+1=C1
E
h
i
t
+
1
+
F
q
i
t
+
1
+
G
h
i
+
1
t
+
1
+
H
q
i
+
1
t
+
1
=
C
2
Eh_i^{t+1}+Fq_i^{t+1}+Gh_{i+1}^{t+1}+Hq_{i+1}^{t+1}=C_2
Ehit+1+Fqit+1+Ghi+1t+1+Hqi+1t+1=C2
λ
=
2
Δ
t
Δ
x
\lambda=\frac{2\Delta t}{\Delta x}
λ=Δx2Δt,则各系数为:
A
=
1
,
B
=
−
λ
,
C
=
1
,
D
=
λ
,
C
1
=
h
i
t
+
h
i
+
1
t
A=1,B=-\lambda,C=1 ,D=\lambda,C_1=h_i^t+h_{i+1}^t
A=1,B=−λ,C=1,D=λ,C1=hit+hi+1t
E
=
−
g
h
λ
,
F
=
1
−
u
λ
+
g
n
2
∣
q
∣
Δ
t
h
7
3
,
G
=
g
h
λ
,
H
=
1
+
u
λ
+
g
n
2
∣
q
∣
Δ
t
h
7
3
,
C
2
=
q
i
t
+
q
i
+
1
t
+
2
g
h
Δ
t
S
0
E=-gh\lambda,F=1-u\lambda+\frac{gn^2|q|\Delta t}{h^{\frac{7}{3}}},G=gh\lambda,H=1+u\lambda+\frac{gn^2|q|\Delta t}{h^{\frac{7}{3}}},C_2=q_i^t+q_{i+1}^t + 2gh \Delta t S_0
E=−ghλ,F=1−uλ+h37gn2∣q∣Δt,G=ghλ,H=1+uλ+h37gn2∣q∣Δt,C2=qit+qi+1t+2ghΔtS0
边界条件离散
上下边界可以是水深或单宽流量,在五对角矩阵的第一行和最后一行添加即可。
五对角矩阵
A
x
=
B
Ax=B
Ax=B形式如下:
(
1
0
0
0
…
…
…
…
…
…
0
A
1
−
2
B
1
−
2
C
1
−
2
D
1
−
2
E
1
−
2
F
1
−
2
G
1
−
2
H
1
−
2
A
2
−
3
B
2
−
3
C
2
−
3
D
2
−
3
E
2
−
3
F
2
−
3
G
2
−
3
H
2
−
3
A
(
n
−
1
)
−
n
B
(
n
−
1
)
−
n
C
(
n
−
1
)
−
n
D
(
n
−
1
)
−
n
E
(
n
−
1
)
−
n
F
(
n
−
1
)
−
n
G
(
n
−
1
)
−
n
H
(
n
−
1
)
−
n
0
…
…
…
…
…
…
0
0
1
0
)
⏟
A
(
h
1
q
1
h
2
q
2
h
3
q
3
h
n
q
n
)
⏟
x
=
(
h
1
(
t
)
C
1
(
1
−
2
)
C
2
(
1
−
2
)
C
1
(
2
−
3
)
C
2
(
2
−
3
)
C
1
[
(
n
−
1
)
−
n
]
C
2
[
(
n
−
1
)
−
n
]
h
n
(
t
)
)
⏟
B
\underbrace{\begin{pmatrix} 1&0&0&0&\dots&\dots&\dots&\dots&\dots&\dots&0\\[4pt] A_{1-2}&B_{1-2}&C_{1-2}&D_{1-2}\\[4pt] E_{1-2} & F_{1-2} & G_{1-2}&H_{1-2} \\[4pt] & & A_{2-3} & B_{2-3} &C_{2-3} &D_{2-3} \\[4pt] & & E_{2-3} & F_{2-3} &G_{2-3} &H_{2-3} \\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & &&&&&& A_{(n-1)-n} & B_{(n-1)-n} &C_{(n-1)-n} &D_{(n-1)-n} \\[4pt] & & &&&&&E_{(n-1)-n} & F_{(n-1)-n} &G_{(n-1)-n} &H_{(n-1)-n} \\[4pt] 0&\dots &\dots &\dots&\dots&\dots&\dots&0 & 0 &1 &0 \\[4pt] \end{pmatrix}}_{A} \underbrace{\begin{pmatrix} h_{1} \\[4pt] q_{1} \\[4pt] h_{2} \\[4pt] q_{2} \\[4pt] h_{3} \\[4pt] q_{3} \\[4pt] \\[4pt] \\[4pt] \\[4pt] \\[4pt] h_{n} \\[4pt] q_{n} \\[4pt] \end{pmatrix}}_{x}=\underbrace{\begin{pmatrix} h_{1}(t) \\[4pt] C_{1(1-2)} \\[4pt] C_{2(1-2)} \\[4pt] C_{1(2-3)} \\[4pt] C_{2(2-3)} \\[4pt] \\[4pt] \\[4pt] \\[4pt] \\[4pt] C_{1[(n-1)-n]} \\[4pt] C_{2[(n-1)-n]} \\[4pt] h_{n}(t) \\[4pt] \end{pmatrix}}_{B}
A
1A1−2E1−200B1−2F1−2…0C1−2G1−2A2−3E2−3…0D1−2H1−2B2−3F2−3……C2−3G2−3……D2−3H2−3…………A(n−1)−nE(n−1)−n0…B(n−1)−nF(n−1)−n0…C(n−1)−nG(n−1)−n10D(n−1)−nH(n−1)−n0
x
h1q1h2q2h3q3hnqn
=B
h1(t)C1(1−2)C2(1−2)C1(2−3)C2(2−3)C1[(n−1)−n]C2[(n−1)−n]hn(t)
模型求解
利用Matlab矩阵求解器可以得到五对角矩阵的解。
Matlab程序
在Matlab中建立脚本,粘贴以下代码可直接运行。版本:MATLAB R2016a
%========================================================
%程序说明
%2019.11.10
%模拟工况:单一河道,断面都为规则的矩形,长20km,糙率n=0.025,
% 初始水深为4m,初始流量为0,断面间距500m,时间步长300s,模拟总时间24h
%边界条件:上边界条件:水深在24h内由4m逐步增加到6m;
% 下边界条件:水深恒定为4m。
% 无旁侧如流。
%=========================================================
%变量定义及初始化
clc,clear;
%L:河长,T:时间,dx:断面间距,dt:时间步长,Hd:下游水深,h0:初始水深,q0:初始流量,n0:糙率
L=20000;T=24*3600;dx=500;dt=300;Hd=4;h0=4;q0=0;g=9.81;n0=0.025;
%n:河道分段数,N:断面数*2即hi,qi变量数,lambda:网格比,a,b,c,d:连续方程离散后的系数,t:时间分段数,Hu:上游水深
n=L/dx;N=(n+1)*2;lambda=2*dt/dx;a=1;b=-lambda;c=1;d=lambda;t=T/dt;Hu=linspace(4,6,t);
%x:hi,qi解向量,B:五对角矩阵方程的常数项,I,J,K,M,O:五对角矩阵A的五列,
x=zeros(1,N);B=zeros(1,N);I=zeros(1,N);J=zeros(1,N-1);K=zeros(1,N-2);M=zeros(1,N-1);O=zeros(1,N-2);
%H,Q分别为模拟时段内的水深和流量过程,X:水深、单宽流量组合矩阵
H=zeros(t+1,n+1);Q=zeros(t+1,n+1);X=zeros(N,t+1);
%初始化hi,qi解向量x
for i=1:2:N
x(i)=h0;x(i+1)=q0;
end
%上下边界水位初始化
x(1)=Hu(1);
x(N-1)=Hd;
%将初始时刻的解向量放入组合矩阵中存储
X(:,1)=x';
%五对角矩阵上二向量赋值,此向量为常向量,在历次迭代中不变
for i=1:2:N-2
K(i)=0;K(i+1)=d;
end
%初始化B,I,J,M向量
B(N)=Hd;
I(1)=1;I(N)=0;
J(1)=0;M(N-1)=1;
%=========================================================
%程序主体部分
for k=1:1:t %按时间层循环
%注意上边界水深是随时间变化的
x(1)=Hu(k);B(1)=Hu(k);
%B,I,J向量赋值
for i=1:2:N-2
%B(i+1),B(i+2)分别为河段上下断面水深和流量之和即C1,C2
B(i+1)=x(i)+x(i+2);B(i+2)=x(i+1)+x(i+3);
%五对角矩阵顶角向量赋值,I(i+1)、I(i+2)分别为方程组系数B,G,h为河段平均水深
I(i+1)=b;h=(x(i)+x(i+2))/2;I(i+2)=g*h*lambda;
%五对角矩阵上一向量赋值,J(i+1)、J(i+2)分别为方程组系数C,H,q为河段平均单宽流量
%注意h^(7/3)采用nthroot(h,3/7)计算,是为了避免出现复数情况
J(i+1)=c;q=(x(i+1)+x(i+3))/2;u=(x(i+1)/x(i)+x(i+3)/x(i+2))/2;J(i+2)=1+u*lambda+g*n0*n0*dt*abs(q)/nthroot(h,3/7);
end
%五对角矩阵下一向量M赋值,M(i)、M(i+1)分别为方程组系数A,F,F=H-2u*lambda
for i=1:2:N-2
M(i)=a;M(i+1)=J(i+2)-2*u*lambda;
end
%五对角矩阵下二向量O赋值,O(i)为方程组系数E,E=-G,O(i+1)为0已经初始化
for i=1:2:N-2
O(i)=-I(i+2);
end
%组合I,J,K,M,O向量得到五对角矩阵A
A=diag(I)+diag(J,1)+diag(K,2)+diag(M,-1)+diag(O,-2);
%解五对角矩阵差分方程组,并以列形式存储在组合矩阵X中
X(:,k+1)=A\B';
%这一次的解作为下一次循环的初始值
x=X(:,k+1)';
end
%=========================================================
%从组合矩阵X中提取出水深H和单宽流量Q
%X的每一列表示一个时段的计算结果
%H,Q的每一行表示一个时段的计算结果
for i=1:2:N
j=(i+1)/2;
H(:,j)=X(i,:)';Q(:,j)=X(i+1,:)';
end
%绘制出水深和单宽流量三维图
figure(1);
mesh(H);
xlabel('河段n/500m');
ylabel('时间k/300s');
zlabel('水深/m');
figure(2);
mesh(Q);
xlabel('河段n/500m');
ylabel('时间k/300s');
zlabel('单宽流量(m^2/s)');
%=========================================================
作者简介
很高兴认识您!
我叫卢家波,河海大学水文学及水资源博士研究生,研究兴趣为高效洪水淹没预测、洪水灾害预警、机器学习、替代模型和降阶模型。
变化环境下,极端洪水事件多发,我希望能通过研究为水灾害防御做出贡献,为人民服务。
欢迎交流讨论和研究合作,vx Jiabo_Lu。
主页 https://lujiabo98.github.io
简历 https://lujiabo98.github.io/file/CV_JiaboLu_zh.pdf
博客 https://blog.csdn.net/weixin_43012724?type=blog
来信请说明博客标题及链接,谢谢。