第二章 压力速度耦合算法
2.1 压力泊松方程
在第一章中,我们推导出了不可压缩牛顿流体的控制方程:
- 连续性方程
∇ ⋅ U ⃗ = 0 (1.8) \nabla \cdot \vec U =0 \tag{1.8} ∇⋅U=0(1.8)- 动量方程
∂ U ⃗ ∂ t + ∇ ⋅ ( U ⃗ U ⃗ ) = ∇ ⋅ ν ∇ U ⃗ − ∇ p ρ + g ⃗ (1.15) {{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U-\nabla {p\over\rho}+ \vec g \tag{1.15} ∂t∂U+∇⋅(UU)=∇⋅ν∇U−∇ρp+g(1.15)
显而易见,这组方程里有两个未知量:速度 U ⃗ \vec U U和压力 p p p。未知量个数与控制方程相等,理论上方程可解,但实际求解则要困难很多,最主要的问题来自控制方程组缺少直接求解压力的方程。具体来说,速度不可能通过方程1.8确定,只能使用方程1.15来求解速度,如果使用方程1.15来求解速度,就意味着必须要使用方程1.8来求解压力,但这是不可能的,因为方程1.8根本不包含压力项。
为了解决压力和速度的耦合求解问题,历史上存在两种处理思路。其中一种思路是将方程(1.8)与(1.15)处理为一个大方程后进行整体求解,称为耦合式解法。另一类方法是将两个方程进行分别求解,称为分离式解法,本文之后要介绍的SIMPLE和PISO等压力速度耦合算法算法就属于分离式解法。所有这些压力耦合速度算法的基础都是构建求解压力的方程,即构建压力泊松方程。在这一节,我们首先进行这项工作。
为了方便表述,我们对方程(1.15)进行一定的处理,首先,我们令
P
=
p
/
ρ
P=p/\rho
P=p/ρ,其次,暂时忽略重力的作用。此时,仿照方程1.36,方程1.15可以离散为:
a
P
U
⃗
p
n
+
∑
N
a
N
U
⃗
N
n
=
S
⃗
−
∇
P
n
(2.1)
a_P\vec U_p^n+\sum_Na_N\vec U_N^n=\vec S-\nabla P^n \tag{2.1}
aPUpn+N∑aNUNn=S−∇Pn(2.1)
忽略重力出于以下两个考虑。首先,在许多计算中重力的作用本身可以忽略。其次,重力本身是一个显式的已知量,其可以归入
S
⃗
\vec S
S,也可以同压力一起计算,关于这一点我们将在之后进行讨论。对方程进行移项可得:
U
⃗
p
n
=
S
⃗
−
∑
N
a
N
U
⃗
N
n
a
P
−
1
a
P
∇
P
n
(2.2)
\vec U_p^n=\frac{\vec S-\sum_Na_N\vec U_N^n}{a_P}-\frac{1}{a_P}\nabla P^n \tag{2.2}
Upn=aPS−∑NaNUNn−aP1∇Pn(2.2)
定义:
H
⃗
=
S
⃗
−
∑
N
a
N
U
⃗
N
n
H
b
y
A
⃗
=
H
⃗
/
a
P
(2.3)
\begin{array}{l} \vec H=\vec S-\sum_Na_N\vec U_N^n \\ \vec{HbyA}=\vec H/a_P \end{array} \tag{2.3}
H=S−∑NaNUNnHbyA=H/aP(2.3)
将方程(2.3)代入方程(2.2)可得:
U
⃗
p
n
=
H
b
y
A
n
⃗
−
1
a
P
∇
P
n
(2.4)
\vec U_p^n=\vec{HbyA^n}-\frac{1}{a_P}\nabla P^n \tag{2.4}
Upn=HbyAn−aP1∇Pn(2.4)
对
U
⃗
p
n
\vec U_p^n
Upn使用连续性方程进行约束:
∇
⋅
(
U
⃗
p
n
)
=
0
→
∇
⋅
(
H
b
y
A
n
⃗
−
1
a
P
∇
P
n
)
=
0
(2.5)
\nabla\cdot(\vec U_p^n)=0 \rightarrow \nabla\cdot\left(\vec{HbyA^n}-\frac{1}{a_P}\nabla P^n\right)=0 \tag{2.5}
∇⋅(Upn)=0→∇⋅(HbyAn−aP1∇Pn)=0(2.5)
即:
∇
⋅
(
1
a
P
∇
P
n
)
=
∇
⋅
(
H
b
y
A
n
⃗
)
(2.6)
\nabla\cdot\left(\frac{1}{a_P}\nabla P^n\right)=\nabla\cdot\left(\vec{HbyA^n}\right) \tag{2.6}
∇⋅(aP1∇Pn)=∇⋅(HbyAn)(2.6)
式2.6即是用于求解压力的压力泊松方程。
为了方便大家对式2.3中的 H H H算符, H b y A HbyA HbyA算符以及整个推导流程有更好的理解,接下来使用矩阵形式重新推导一遍压力泊松方程。读者应该注意对比矩阵形式和上述形式之间的联系。
对方程
∂
U
⃗
∂
t
+
∇
⋅
(
U
⃗
U
⃗
)
=
∇
⋅
ν
∇
U
⃗
{{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U
∂t∂U+∇⋅(UU)=∇⋅ν∇U进行离散得到:
[
A
]
[
U
⃗
]
=
[
S
⃗
]
(2.7)
[A][\vec U]=[\vec S] \tag{2.7}
[A][U]=[S](2.7)
将
[
A
]
[A]
[A]拆分成对角线元素
[
A
d
i
a
g
]
[A_{diag}]
[Adiag]和非对角元素
[
A
o
f
f
d
i
a
g
]
[A_{offdiag}]
[Aoffdiag]:
[
A
d
i
a
g
]
[
U
⃗
]
+
[
A
o
f
f
d
i
a
g
]
[
U
⃗
]
=
[
S
⃗
]
(2.8)
[A_{diag}][\vec U]+[A_{offdiag}][\vec U]=[\vec S] \tag{2.8}
[Adiag][U]+[Aoffdiag][U]=[S](2.8)
定义
H
H
H算符和
H
b
y
A
HbyA
HbyA算符:
H
(
[
U
⃗
]
)
=
[
S
⃗
]
−
[
A
o
f
f
d
i
a
g
]
[
U
⃗
]
(2.9)
H([\vec U])=[\vec S]-[A_{offdiag}][\vec U] \tag{2.9}
H([U])=[S]−[Aoffdiag][U](2.9)
H
b
y
A
(
[
U
⃗
]
)
=
[
A
d
i
a
g
]
−
1
H
(
U
⃗
)
(2.10)
HbyA([\vec U])={[A_{diag}]}^{-1}H(\vec U) \tag{2.10}
HbyA([U])=[Adiag]−1H(U)(2.10)
因此,对方程
∂
U
⃗
∂
t
+
∇
⋅
(
U
⃗
U
⃗
)
=
∇
⋅
ν
∇
U
⃗
+
∇
P
{{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U+\nabla P
∂t∂U+∇⋅(UU)=∇⋅ν∇U+∇P进行离散可得:
[
U
⃗
]
=
H
b
y
A
(
[
U
⃗
]
)
−
[
A
d
i
a
g
]
−
1
∇
[
P
]
(2.11)
[\vec U]=HbyA([\vec U])-{[A_{diag}]}^{-1}\nabla [P] \tag{2.11}
[U]=HbyA([U])−[Adiag]−1∇[P](2.11)
速度满足连续性方程
∇
⋅
U
⃗
=
0
\nabla \cdot \vec U=0
∇⋅U=0,带入方程(2.11)可得:
∇
⋅
(
H
b
y
A
(
[
U
⃗
]
)
)
=
∇
⋅
(
[
A
d
i
a
g
]
−
1
∇
[
P
]
)
(2.12)
\nabla\cdot\left(HbyA([\vec U])\right)=\nabla\cdot\left({[A_{diag}]}^{-1}\nabla [P]\right) \tag{2.12}
∇⋅(HbyA([U]))=∇⋅([Adiag]−1∇[P])(2.12)
上述矩阵形式的推导过程对与理解代码非常有帮助,在之后关于OpenFOAM的代码解读章节,我们将对上述过程进行再次解释。