PISO,全称为 Pressure-Implicit with Splitting of Operators,是 Issa1 于 1986 年提出的一种求解耦合速度压力的非迭代算法,也是目前应用最广的速度压力解耦方法之一。OpenFOAM 中实现的 PISO 算法与原始版本并不完全相同,本文将对原始的 PISO 算法及 OpenFOAM 中的 PISO 算法进行对比,分析二者的不同之处。

原始 PISO

原始文献2中对 PISO 的描述如下:

The method (called PISO for Pressure-Implicit with Splitting of Operators) utilises the splitting of operations in the solution of the discretised momentum and pressure equations such that the fields obtained at each time-step are close approximations of the exact solution of the difference equations with a formal order of accuracy of the order of powers of $\delta t$ depending on the number of operation-splittings used (two such splittings are proposed herein).

从上面的描述中可以看出,PISO 算法在求解离散的动量方程和压力方程时用到了算子分裂方法。

算子分裂

算子分裂(operator splitting)是用来计算偏微分方程(PDE)的一类方法,其基本思想是对不同的物理过程(对流、扩散等)采用不同的数值方法求解。即按照 PDE 方程中的微分算子将其分解成若干个更为简单的子方程,然后分别处理。常见的算子分裂方法可以分为两类:一类是差分分裂(differential splitting),另一类是代数分裂(algebraic splitting)。关于这二者的区别这里不展开叙述,有兴趣的同学可以参考这个讲义3

PISO 用到的算子分裂是代数算子分裂,其形式如下

\[ \begin{equation} \frac{\partial \phi}{\partial t}+L(\phi) = 0, \quad L = \sum_{s=1}^S L_s \end{equation} \]

上式中,$L_s$ 代表一系列对 $\phi$ 进行的离散操作符。

为了让大家更了解整个求解过程,我们以对流扩散方程为例进行说明。考虑如下形式的对流扩散方程

\[ \begin{equation} \frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf u) - \nabla \cdot (\gamma \nabla \phi) = 0 \end{equation} \]

假设我们采用一个确定的时间离散格式,那么下一时刻值 $\phi^{n+1}$ 可以表示为当前时刻值 $\phi^n$ 的确切表达式,我们用 $\mathcal F$ 来表示这个表达式,即

\[ \begin{equation} \phi^{n+1} = \mathcal F(\phi^n, \Delta t) \end{equation} \]

将对流扩散方程按照操作符进行算子分裂,得到两个微分方程

\[ \begin{equation} \frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf u) = 0 \end{equation} \]

\[ \begin{equation} \frac{\partial \phi}{\partial t} - \nabla \cdot (\gamma \nabla \phi) = 0 \end{equation} \]

对上述两个方程的空间操作符进行离散,保持时间项不变,得到两个半离散方程

\[ \begin{equation} \frac{\partial \phi}{\partial t} + \mathcal C(\phi) = 0 \end{equation} \]

\[ \begin{equation} \frac{\partial \phi}{\partial t} + \mathcal D(\phi) = 0 \end{equation} \]

式中的 $\mathcal C(\phi)$ 和 $\mathcal D(\phi)$ 为离散后的对流项和扩散项。

对这两个半离散方程进行分步求解,先求解对流方程,再求解扩散方程,得到下一时刻值 $\phi^{n+1}$

\[ \begin{equation} \begin{aligned} \phi^{n+1/2} &= \mathcal F_C(\phi^n, \Delta t) \\ \phi^{n+1} &= \mathcal F_D(\phi^{n+1/2}, \Delta t) \end{aligned} \end{equation} \]

代数分裂的目的是对半离散的方程进行分步求解,因此一些文献中也将算子分裂称为分步法(fractional step method)。但实际上这二者并不相同,算子分裂是指将 $L=\sum {L_s}$ 分解的过程,而分步法则多指对时间项离散并进行分步求解的过程。

动量方程的算子分裂

前面介绍的对流扩散方程中的算子分裂是针对单个物理量的,而动量方程中的算子分裂涉及了两个物理量,即速度和压力。

Isaa 文章中的 PISO 算法既可用于可压缩流体也可用于不可压缩流体,这里以不可压缩流体为例来看下原始文献中动量方程的算子分裂。

带有源项的动量方程可以写成以下形式

\[ \begin{equation} \frac{\partial \mathbf u}{\partial t} + \nabla \cdot (\mathbf u \otimes \mathbf u) - \nabla \cdot (\gamma \nabla \mathbf u) = - \frac{1}{\rho} \nabla p + \mathbf S \end{equation} \]

对动量方程进行离散,为了方便说明,时间项用一阶精度的欧拉格式离散。引入 $\mathcal H$ 操作符,令 $-\mathcal H$ 表示对流项和扩散项离散后的成分,即

\[ \frac{1}{V}\int_V \left[\nabla \cdot (\mathbf u \otimes \mathbf u) - \nabla \cdot (\gamma \nabla \mathbf u) \right]dV= \mathbf A \mathbf u = -\mathcal H(\mathbf u) \]

其中 $\mathbf A$ 为对流项和扩散项离散后的系数矩阵。

于是动量方程可以写作

\[ \begin{equation} \frac{\mathbf u^{n+1} - \mathbf u^n}{\Delta t} - \mathcal H(\mathbf u^{n+1}) = - \frac{1}{\rho} \nabla p^{n+1} + \mathbf S \end{equation} \label{eq:mom_eqn} \]

对上式进行算子分裂,引入上标 $^{\ast}$ 、 $^{\ast\ast}$ 和 $^{\ast\ast\ast}$ 表示算子分裂中间步骤产生的速度。

  1. 动量预估

将上一时刻的压力场 $p^n$ 代入 $\eqref{eq:mom_eqn}$,用 $\mathbf u^{\ast}$ 代替 $\mathbf u^{n+1}$,得到

\[ \begin{equation} \frac{\mathbf u^{\ast} - \mathbf u^n}{\Delta t} - \mathcal H(\mathbf u^{\ast}) = - \frac{1}{\rho} \nabla p^{n} + \mathbf S \end{equation} \]

求解上式可以得到一个中间速度 $\mathbf u^{\ast}$,此时的 $\mathbf u^{\ast}$ 不满足连续性方程,需要对其进行修正。

  1. 第一次修正

由于第一步中的压力用的是上一时刻的值,因此得到的速度不满足连续性方程。我们假设一个新的压力场 $p^{\ast}$ 可以得到满足连续性方程的速度场 $\mathbf u^{\ast\ast}$

\[ \begin{equation} \nabla \mathbf u^{\ast\ast} = 0 \end{equation} \label{eq:cont_eqn} \]

则动量方程 $\eqref{eq:mom_eqn}$ 变为

\[ \begin{equation} \frac{\mathbf u^{\ast\ast} - \mathbf u^n}{\Delta t} - \mathcal H(\mathbf u^{\ast}) = - \frac{1}{\rho} \nabla p^{\ast} + \mathbf S \end{equation} \label{eq:1st_mom_eqn} \]

这里 H 操作符用到了动量预估步中得到的速度 $\mathbf u^{\ast}$。

将动量方程 $\eqref{eq:1st_mom_eqn}$ 代入连续性方程 $\eqref{eq:cont_eqn}$,可以得到

\[ \begin{equation} \frac{1}{\rho}\nabla^2 p^{\ast} = \frac{1}{\Delta t}\nabla \mathbf u^n + \nabla \mathcal H (\mathbf u^\ast) + \nabla\mathbf S \end{equation} \label{eq:ppe} \]

上式被称为压力泊松方程,求解该方程可以得到新的压力 $p^\ast$ ,将 $p^\ast$ 代入方程 $\eqref{eq:1st_mom_eqn}$ 得到新的速度 $\mathbf u^{\ast\ast}$,该速度满足连续性方程。

  1. 第二次修正

由于第一次修正步中的 H 操作符采用了显式操作,因此需要进行第二次修正,以减少误差。第二次修正的基本步骤同第一次修正相同。假设一个新的压力场 $p^{\ast\ast}$ 以及对应的满足连续性方程的速度场 $\mathbf u^{\ast\ast\ast}$

\[ \begin{equation} \nabla \mathbf u^{\ast\ast\ast} = 0 \end{equation} \]

则动量方程用以下形式表示

\[ \begin{equation} \frac{\mathbf u^{\ast\ast\ast} - \mathbf u^n}{\Delta t} - \mathcal H(\mathbf u^{\ast\ast}) = - \frac{1}{\rho} \nabla p^{\ast\ast} + \mathbf S \end{equation} \label{eq:2nd_mom_eqn} \]

相应地,压力泊松方程变为

\[ \begin{equation} \frac{1}{\rho}\nabla^2 p^{\ast\ast} = \frac{1}{\Delta t}\nabla \mathbf u^n + \nabla \mathcal H (\mathbf u^{\ast\ast}) + \mathbf S \end{equation} \]

求解压力泊松方程得到新的压力场 $p^{\ast\ast}$,代入方程 $\eqref{eq:2nd_mom_eqn}$ 得到新的速度场 $\mathbf u^{\ast\ast\ast}$。

当然还可以进行第三、第四次甚至更多的修正步计算,然而 Issa 认为对于大部分流动问题,两次修正已经可以得到足够精确的解,无需进行更多修正步计算。

OpenFOAM 中的 PISO

OpenFOAM 中的 PISO 与原始 PISO 不同,主要体现在以下两个方面。

动量预估可选

OpenFOAM 中的 PISO 算法的动量预估步骤是可选的,通过字典文件 constant/fvSolutionPISOPIMPLE 子字典的关键词 momentomPredictor 控制。参考 icoFoam 求解器中的代码:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
        // Momentum predictor

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

只有当 momentumPredictor 设置为 yes 时才会求解动量方程。而当其设为 no 时,跳过动量预估步骤,不求解动量方程,压力泊松方程 $\eqref{eq:ppe}$ 中 $\mathcal H$ 操作符中的速度直接用上一时刻或上一迭代步的值4

momentumPredictor 的默认值是 true。但是,OpenFOAM 自带算例中的绝大部分 interFoam 算例中都关闭了动量预估,对此,Henry 的解释是5

Doing a momentum predictor is not an essential step for the convergence of the PISO corrector loop although it is sometimed benefitial but not always. For example in very low-Re flows the momentum predictor step can be severely detrimental to the convergence. For interFoam I found that the momentum predictor step did not improve the convergence behaviour and is a bit complicated to include so I removed it for simplicity. I have reinstated it for the 1.2 release just so people can find out for themselves if it is helpful or not.

Jasak 认为如果动量方程中的扩散项很重要,则需要求解动量方程6

It depends on the importance of the laplacian in the momentum equation: if you have important diffusivity in the momentum, you should solve the momentum equation implicitly.

H 操作符中的系数矩阵

OpenFOAM 中 $\mathcal H$ 操作符包括时间项、对流项和扩散项的影响,而 Issa 提出的 $\mathcal H$ 操作符只包含对流项和扩散项。求解器在构建 UEqn 时都计入了 fvm::ddt 时间项的影响。

然而,在 foam-extend 中,在计算 $\mathcal H$ 操作符时并没有考虑时间项,这个实现同 Issa 的原始版本一致。同样以 icoFoam 为例:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
        // Convection-diffusion matrix
        fvVectorMatrix HUEqn
        (
            fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        // Time derivative matrix
        fvVectorMatrix ddtUEqn(fvm::ddt(U));

        if (piso.momentumPredictor())
        {
            solve(ddtUEqn + HUEqn == -fvc::grad(p));
        }

        // Prepare clean Ap without time derivative contribution
        // HJ, 26/Oct/2015
        volScalarField aU = HUEqn.A();

此外,OpenFOAM 采用的是非结构化的同位网格,非结构化网格可能出现非正交网格,需要进行非正交修正。而同位网格为了消除压力振荡,需要进行动量插值。当然非正交修正和动量插值已经不属于 PISO 算法的范畴,关于这两个话题我们以后有机会再详聊。

参考


  1. Issa, R. I. (1986). Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics, 62(1), 40–65. [return]
  2. Issa, R. I. (1986). Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics, 62(1), 40–65. [return]
  3. http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture11.pdf [return]
  4. https://www.cfd-online.com/Forums/openfoam/90500-momentum-predictor.html#post316780 [return]
  5. https://www.cfd-online.com/Forums/openfoam-solving/60558-why-sometimes-momentum-predictor-step-not-performed.html#post203065 [return]
  6. https://www.cfd-online.com/Forums/openfoam-solving/58067-when-set-momentumpredictor-yes.html#post191166 [return]