【CFD小工坊】浅水模型的边界条件
- 前言
- 处理边界条件的原则
- 边界处水力要素的计算
- 水位边界条件
- 单宽流量边界条件
- 流量边界条件
- 固壁边界条件
- 参考文献
前言
在浅水方程的离散及求解方法一篇中,我们学习了三角形网格各边通量值及源项的求解。但仍有一个问题没有解决:对于边界处的网格,模型边界对应的网格边的通量求解。
对此,我们借鉴王志力1的研究,学习各类边界条件下,网格边的通量的求解。
处理边界条件的原则
对于浅水水域,常见的边界有水位边界与流量边界。在此,我们假设网格
i
i
i的第
j
j
j条边对应了模型的边界,设边界上的水位为
h
i
j
∗
h_{ij}^*
hij∗,垂直(外法线方向)和平行网格边的流速为
u
~
i
j
∗
\tilde{u}_{ij}^*
u~ij∗和
v
~
i
j
∗
\tilde{v}_{ij}^*
v~ij∗。为简便起见,以下我们将
h
i
j
∗
h_{ij}^*
hij∗简记为
h
∗
h^*
h∗,将
u
~
i
j
∗
\tilde{u}_{ij}^*
u~ij∗和
v
~
i
j
∗
\tilde{v}_{ij}^*
v~ij∗简记为
u
~
L
∗
\tilde{u}_L^*
u~L∗和
v
~
L
∗
\tilde{v}_L^*
v~L∗
根据一维方程的特征线理论,沿着特征线方向有特征不变量,最终会得到如下关系:
u
~
L
∗
+
2
c
L
=
u
∗
+
2
c
∗
\tilde{u}_L^*+2c_L = u^* + 2c^*
u~L∗+2cL=u∗+2c∗
上式即确定边界条件时所要满足的原则。其中,
c
L
=
g
h
L
c_L=\sqrt{gh_L}
cL=ghL,
c
∗
=
g
h
∗
c^*=\sqrt{gh^*}
c∗=gh∗。
边界处水力要素的计算
在模型中,边界处的水力要素的计算步骤如下:
- 根据笛卡尔坐标系下边界处的 u L u_L uL和 v L v_L vL转化为网格边界外法线方向和切向的 u ~ L ∗ \tilde{u}_L^* u~L∗和 v ~ L ∗ \tilde{v}_L^* v~L∗。
- 给定边界处的 u ∗ u^* u∗或 h ∗ h^* h∗。此处的 u ∗ u^* u∗值可通过流量边界条件转化而来。
- 根据式 u ~ L ∗ + 2 c L = u ∗ + 2 c ∗ \tilde{u}_L^*+2c_L = u^* + 2c^* u~L∗+2cL=u∗+2c∗来确定网格边上的其它变量。例如,对于水位条件, h ∗ h^* h∗已知,我们需要通过上式确定 u ∗ u^* u∗。
- 根据浅水方程的对流项确定通量值。
水位边界条件
对于边界条件,
h
∗
h^*
h∗已知,则:
u
∗
=
u
~
L
∗
+
2
c
L
−
2
c
∗
=
u
~
L
∗
+
2
g
h
L
−
2
g
h
∗
u^* = \tilde{u}_L^*+2c_L - 2c^*=\tilde{u}_L^*+ 2\sqrt{gh_L} -2\sqrt{gh^*}
u∗=u~L∗+2cL−2c∗=u~L∗+2ghL−2gh∗
之后将局部坐标系的
u
∗
u^*
u∗和
v
∗
v^*
v∗转化为全局笛卡尔坐标系下的
u
b
u_b
ub和
v
b
v_b
vb;记
h
b
=
h
∗
h_b=h^*
hb=h∗。则边界处的通量为:
(
F
n
)
b
=
E
(
U
b
)
n
x
+
G
(
U
b
)
n
y
=
n
x
(
h
u
b
h
u
b
2
+
g
h
b
2
2
h
u
b
v
b
)
+
n
y
(
h
v
b
z
h
u
b
v
b
h
v
b
2
+
g
h
b
2
2
)
(\bold{F}_n)_b = \bold{E(U_b)} n_x+ \bold{G(U_b)} n_y = n_x \left( \begin{matrix} hu_b \\ hu_b^2+\dfrac{gh_b^2}{2} \\ hu_b v_b \end{matrix} \right) + n_y \left( \begin{matrix} hv_bz \\ hu_b v_b \\ hv_b^2+\dfrac{gh_b^2}{2} \end{matrix} \right)
(Fn)b=E(Ub)nx+G(Ub)ny=nx
hubhub2+2ghb2hubvb
+ny
hvbzhubvbhvb2+2ghb2
式中,
(
n
x
,
n
y
)
(n_x, n_y)
(nx,ny)表示边界处外法线方向。
单宽流量边界条件
给定网格边的单宽流量
q
=
h
∗
u
∗
q=h^*u^*
q=h∗u∗,则有:
u
~
L
∗
+
2
c
L
=
u
∗
+
2
c
∗
=
q
h
∗
+
2
g
h
∗
=
q
c
∗
2
/
g
+
2
g
h
∗
\tilde{u}_L^*+2c_L = u^* + 2c^* = \dfrac{q}{h^*} + 2\sqrt{gh^*} = \dfrac{q}{{c^*}^2/g} + 2\sqrt{gh^*}
u~L∗+2cL=u∗+2c∗=h∗q+2gh∗=c∗2/gq+2gh∗
化简后,上述方程为
c
∗
c^*
c∗的一元三次方程:
2
c
∗
3
−
(
u
L
+
2
g
h
L
)
c
∗
2
−
g
q
=
0
2{c^*}^3 - (u_L + 2\sqrt{gh_L}){c^*}^2 - gq = 0
2c∗3−(uL+2ghL)c∗2−gq=0
求解后可得
h
∗
=
c
∗
2
/
g
h^*={c^*}^2/g
h∗=c∗2/g。同理,我们可求得
h
b
h_b
hb、
u
b
u_b
ub和
v
b
v_b
vb,以及通量
(
F
n
)
b
(\bold{F}_n)_b
(Fn)b。
注意:在设置边界时,我们需要设定单宽流量的方向;对于入流边界,单宽流量方向与边界外法线方向相反,则
q
<
0
q<0
q<0;反之,对于出流边界,
q
>
0
q>0
q>0。
流量边界条件
若流量给定在一个网格的边上,我们可以先求该边界的单宽流量
q
q
q,之后按照上一小节等同的办法处理边界。若指定的边界条件涉及到m条连续的网格边(如下图边界蓝色边所示),组需要先求出每个对应网格边的单宽流量,之后再按单宽流量边界条件处理方法进行计算。
对于这m条边界上的总流量
Q
Q
Q,某一网格
i
i
i边上的单宽流量
q
i
q_i
qi是:
q
i
=
L
i
′
h
i
1.5
C
i
∑
k
=
1
m
L
k
′
h
k
1.5
C
k
Q
q_i = \dfrac{L'_i h_i^{1.5}C_i}{\sum^{m}_{k=1} L'_k h_k^{1.5}C_k} Q
qi=∑k=1mLk′hk1.5CkLi′hi1.5CiQ
式中,
L
′
L'
L′表示流量边界对应网格边的边长,
h
h
h表示对应网格的水深,
C
C
C表示对应网格的摩阻力项,有
C
=
h
1
/
6
n
C = \dfrac{h^{1/6}}{n}
C=nh1/6,n为曼宁系数。
之后根据求出的单款流量,依次处理每个边界网格的通量值。
固壁边界条件
在静止的固壁边界上,我们采用无滑移边界条件,即
u
b
u_b
ub和
v
b
v_b
vb均为0,故:
(
F
n
)
b
=
(
0
g
h
L
2
2
n
x
g
h
L
2
2
n
y
)
(\bold{F}_n)_b = \left( \begin{matrix} 0 \\ \dfrac{gh_L^2}{2}n_x \\ \dfrac{gh_L^2}{2}n_y \end{matrix} \right)
(Fn)b=
02ghL2nx2ghL2ny
参考文献
王志力. 基于Godunov和Semi-Lagrangian法的二、三维浅水方程的非结构化网格离散研究[D]. 辽宁:大连理工大学,2005. ↩︎