在第一章中，我们推导出了不可压缩牛顿流体的控制方程：

显而易见，这组方程里有两个未知量：速度 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} aP​U pn​+N∑​aN​U Nn​=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} U pn​=aP​S −∑N​aN​U Nn​​−aP​1​∇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 −∑N​aN​U Nn​HbyA ​=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} U pn​=HbyAn ​−aP​1​∇Pn(2.4) 对 U ⃗ p n \vec U_p^n U pn​使用连续性方程进行约束： ∇ ⋅ ( 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} ∇⋅(U pn​)=0→∇⋅(HbyAn ​−aP​1​∇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} ∇⋅(aP​1​∇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 ​+∇⋅(U U )=∇⋅ν∇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 ​+∇⋅(U U )=∇⋅ν∇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的代码解读章节，我们将对上述过程进行再次解释。