2.1 压力泊松方程

第二章 压力速度耦合算法

2.1 压力泊松方程

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

  • 连续性方程

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

  • 动量方程

(1.15)Ut+(UU)=νUpρ+g{{\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}

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

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

为了方便表述,我们对方程(1.15)进行一定的处理,首先,我们令P=p/ρP=p/\rho,其次,暂时忽略重力的作用。此时,仿照方程1.36,方程1.15可以离散为:

(2.1)aPUpn+NaNUNn=SPna_P\vec U_p^n+\sum_Na_N\vec U_N^n=\vec S-\nabla P^n \tag{2.1}

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

(2.2)Upn=SNaNUNnaP1aPPn\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}

定义:

(2.3)H=SNaNUNnHbyA=H/aP\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)可得:

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

Upn\vec U_p^n使用连续性方程进行约束:

(2.5)(Upn)=0(HbyAn1aPPn)=0\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}

即:

(2.6)(1aPPn)=(HbyAn)\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中的HH算符,HbyAHbyA算符以及整个推导流程有更好的理解,接下来使用矩阵形式重新推导一遍压力泊松方程。读者应该注意对比矩阵形式和上述形式之间的联系。

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

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

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

(2.8)[Adiag][U]+[Aoffdiag][U]=[S][A_{diag}][\vec U]+[A_{offdiag}][\vec U]=[\vec S] \tag{2.8}

定义HH算符和HbyAHbyA算符:

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

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

因此,对方程Ut+(UU)=νU+P{{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U+\nabla P进行离散可得:

(2.11)[U]=HbyA([U])[Adiag]1[P][\vec U]=HbyA([\vec U])-{[A_{diag}]}^{-1}\nabla [P] \tag{2.11}

速度满足连续性方程U=0\nabla \cdot \vec U=0,带入方程(2.11)可得:

(2.12)(HbyA([U]))=([Adiag]1[P])\nabla\cdot\left(HbyA([\vec U])\right)=\nabla\cdot\left({[A_{diag}]}^{-1}\nabla [P]\right) \tag{2.12}

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