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

第二章 压力速度耦合算法

2.1 压力泊松方程

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

• 连续性方程
  $$\nabla \cdot \vec{U}=0\text{\hspace{1em}(1.8)}$$
• 动量方程
  $$\frac{\partial \vec{U}}{\partial t}+\nabla \cdot (\vec{U}\vec{U})=\nabla \cdot \nu \nabla \vec{U}-\nabla \frac{p}{\rho }+\vec{g}\text{\hspace{1em}(1.15)}$$

显而易见，这组方程里有两个未知量：速度$\vec{U}$和压力$p$。未知量个数与控制方程相等，理论上方程可解，但实际求解则要困难很多，最主要的问题来自控制方程组缺少直接求解压力的方程。具体来说，速度不可能通过方程1.8确定，只能使用方程1.15来求解速度，如果使用方程1.15来求解速度，就意味着必须要使用方程1.8来求解压力，但这是不可能的，因为方程1.8根本不包含压力项。

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

为了方便表述，我们对方程(1.15)进行一定的处理，首先，我们令$P=p/\rho$，其次，暂时忽略重力的作用。此时，仿照方程1.36，方程1.15可以离散为：
$$a_P{\vec{U}}_p^n+\sum _N{a}_N{\vec{U}}_N^n=\vec{S}-\nabla P^n\text{\hspace{1em}(2.1)}$$
忽略重力出于以下两个考虑。首先，在许多计算中重力的作用本身可以忽略。其次，重力本身是一个显式的已知量，其可以归入$\vec{S}$，也可以同压力一起计算，关于这一点我们将在之后进行讨论。对方程进行移项可得：
$${\vec{U}}_p^n=\frac{\vec{S}-{\sum }_N{a}_N{\vec{U}}_N^n}{a_P}-\frac{1}{a_P}\nabla P^n\text{\hspace{1em}(2.2)}$$
定义：
$$\begin{array}{l}\vec{H}=\vec{S}-{\sum }_N{a}_N{\vec{U}}_N^n\\ \overrightarrow{HbyA}=\vec{H}/a_P\end{array}\text{\hspace{1em}(2.3)}$$
将方程(2.3)代入方程(2.2)可得：
$${\vec{U}}_p^n=\overrightarrow{Hby{A}^n}-\frac{1}{a_P}\nabla P^n\text{\hspace{1em}(2.4)}$$
对${\vec{U}}_p^n$使用连续性方程进行约束：
$$\nabla \cdot ({\vec{U}}_p^n)=0\to \nabla \cdot \left(\overrightarrow{Hby{A}^n}-\frac{1}{a_P}\nabla P^n\right)=0\text{\hspace{1em}(2.5)}$$
即：
$$\nabla \cdot \left(\frac{1}{a_P}\nabla P^n\right)=\nabla \cdot \left(\overrightarrow{Hby{A}^n}\right)\text{\hspace{1em}(2.6)}$$
式2.6即是用于求解压力的压力泊松方程。

为了方便大家对式2.3中的$H$算符，$HbyA$算符以及整个推导流程有更好的理解，接下来使用矩阵形式重新推导一遍压力泊松方程。读者应该注意对比矩阵形式和上述形式之间的联系。

对方程$\frac{\partial \vec{U}}{\partial t}+\nabla \cdot (\vec{U}\vec{U})=\nabla \cdot \nu \nabla \vec{U}$进行离散得到：
$$[A][\vec{U}]=[\vec{S}]\text{\hspace{1em}(2.7)}$$
将$[A]$拆分成对角线元素$[A_{diag}]$和非对角元素$[A_{offdiag}]$:
$$[A_{diag}][\vec{U}]+[A_{offdiag}][\vec{U}]=[\vec{S}]\text{\hspace{1em}(2.8)}$$
定义$H$算符和$HbyA$算符：
$$H([\vec{U}])=[\vec{S}]-[A_{offdiag}][\vec{U}]\text{\hspace{1em}(2.9)}$$
$$HbyA([\vec{U}])={[A_{diag}]}^{-1}H(\vec{U})\text{\hspace{1em}(2.10)}$$
因此，对方程$\frac{\partial \vec{U}}{\partial t}+\nabla \cdot (\vec{U}\vec{U})=\nabla \cdot \nu \nabla \vec{U}+\nabla P$进行离散可得：
$$[\vec{U}]=HbyA([\vec{U}])-{[A_{diag}]}^{-1}\nabla [P]\text{\hspace{1em}(2.11)}$$
速度满足连续性方程$\nabla \cdot \vec{U}=0$，带入方程(2.11)可得：
$$\nabla \cdot \left(HbyA([\vec{U}])\right)=\nabla \cdot \left({[A_{diag}]}^{-1}\nabla [P]\right)\text{\hspace{1em}(2.12)}$$
上述矩阵形式的推导过程对与理解代码非常有帮助，在之后关于OpenFOAM的代码解读章节，我们将对上述过程进行再次解释。