2.1 压力泊松方程

第二章压力速度耦合算法

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

连续性方程

\nabla \cdot \vec U =0 \tag{1.8}

动量方程

{{\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}

显而易见，这组方程里有两个未知量：速度 $\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_Na_N\vec U_N^n=\vec S-\nabla P^n \tag{2.1}

忽略重力出于以下两个考虑。首先，在许多计算中重力的作用本身可以忽略。其次，重力本身是一个显式的已知量，其可以归入 $\vec S$ ，也可以同压力一起计算，关于这一点我们将在之后进行讨论。对方程进行移项可得：

\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}

定义：

\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}

将方程(2.3)代入方程(2.2)可得：

\vec U_p^n=\vec{HbyA^n}-\frac{1}{a_P}\nabla P^n \tag{2.4}

对 $\vec U_p^n$ 使用连续性方程进行约束：

\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}

即：

\nabla\cdot\left(\frac{1}{a_P}\nabla P^n\right)=\nabla\cdot\left(\vec{HbyA^n}\right) \tag{2.6}

式2.6即是用于求解压力的压力泊松方程。

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

对方程 ${{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U$ 进行离散得到：

[A][\vec U]=[\vec S] \tag{2.7}

将 $[A]$ 拆分成对角线元素 $[A_{diag}]$ 和非对角元素 $[A_{offdiag}]$ :

[A_{diag}][\vec U]+[A_{offdiag}][\vec U]=[\vec S] \tag{2.8}

定义 $H$ 算符和 $HbyA$ 算符：

H([\vec U])=[\vec S]-[A_{offdiag}][\vec U] \tag{2.9}

HbyA([\vec U])={[A_{diag}]}^{-1}H(\vec U) \tag{2.10}

因此，对方程 ${{\partial \vec U}\over{\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] \tag{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) \tag{2.12}

上述矩阵形式的推导过程对与理解代码非常有帮助，在之后关于OpenFOAM的代码解读章节，我们将对上述过程进行再次解释。

CloudBird's Idea Space