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