2.1 压力泊松方程（OpenFOAM理论笔记系列）

第二章压力速度耦合算法

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 UU和压力ppp。未知量个数与控制方程相等，理论上方程可解，但实际求解则要困难很多，最主要的问题来自控制方程组缺少直接求解压力的方程。具体来说，速度不可能通过方程1.8确定，只能使用方程1.15来求解速度，如果使用方程1.15来求解速度，就意味着必须要使用方程1.8来求解压力，但这是不可能的，因为方程1.8根本不包含压力项。

为了解决压力和速度的耦合求解问题，历史上存在两种处理思路。其中一种思路是将方程(1.8)与(1.15)处理为一个大方程后进行整体求解，称为耦合式解法。另一类方法是将两个方程进行分别求解，称为分离式解法，本文之后要介绍的SIMPLE和PISO等压力速度耦合算法算法就属于分离式解法。所有这些压力耦合速度算法的基础都是构建求解压力的方程，即构建压力泊松方程。在这一节，我们首先进行这项工作。

为了方便表述，我们对方程(1.15)进行一定的处理，首先，我们令P=p/ρP=p/\rhoP=p/ρ，其次，暂时忽略重力的作用。此时，仿照方程1.36，方程1.15可以离散为：
aPU⃗pn+∑NaNU⃗Nn=S⃗−∇Pn(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 SS，也可以同压力一起计算，关于这一点我们将在之后进行讨论。对方程进行移项可得：
U⃗pn=S⃗−∑NaNU⃗NnaP−1aP∇Pn(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⃗−∑NaNU⃗NnHbyA⃗=H⃗/aP(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⃗pn=HbyAn⃗−1aP∇Pn(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⃗pn\vec U_p^nUpn使用连续性方程进行约束：
∇⋅(U⃗pn)=0→∇⋅(HbyAn⃗−1aP∇Pn)=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)
即：
∇⋅(1aP∇Pn)=∇⋅(HbyAn⃗)(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中的HHH算符，HbyAHbyAHbyA算符以及整个推导流程有更好的理解，接下来使用矩阵形式重新推导一遍压力泊松方程。读者应该注意对比矩阵形式和上述形式之间的联系。

对方程∂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]拆分成对角线元素[Adiag][A_{diag}][Adiag]和非对角元素[Aoffdiag][A_{offdiag}][Aoffdiag]:
[Adiag][U⃗]+[Aoffdiag][U⃗]=[S⃗](2.8)[A_{diag}][\vec U]+[A_{offdiag}][\vec U]=[\vec S] \tag{2.8} [Adiag][U]+[Aoffdiag][U]=[S](2.8)
定义HHH算符和HbyAHbyAHbyA算符：
H([U⃗])=[S⃗]−[Aoffdiag][U⃗](2.9)H([\vec U])=[\vec S]-[A_{offdiag}][\vec U] \tag{2.9} H([U])=[S]−[Aoffdiag][U](2.9)
HbyA([U⃗])=[Adiag]−1H(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⃗]=HbyA([U⃗])−[Adiag]−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)可得：
∇⋅(HbyA([U⃗]))=∇⋅([Adiag]−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的代码解读章节，我们将对上述过程进行再次解释。