PISO 算法

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

原始 PISO

原始文献中对 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）。关于这二者的区别这里不展开叙述，有兴趣的同学可以参考这个讲义²。

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) \newline \phi^{n+1} &= \mathcal F_D(\phi^{n+1/2}, \Delta t) \end{aligned} \end{equation} $$

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

动量方程的算子分裂

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

Issa 文章中的 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(\mathbf u)$ 表示对流项和扩散项离散后的通量加和，即

$$ \begin{equation} \mathcal H(\mathbf u) = A_m \mathbf u_{m} \end{equation} $$

其中，下标 $m$ 表示有限差分的节点编号，上式表示对节点的所有通量进行加和。

对于有限体积法，$\mathcal H$ 的定义不太一样。参考 Issa 的另一篇文献³，采用有限体积法进行离散后，对流项和扩散项形成了系数矩阵 $\mathbf A$

$$ \begin{equation} \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 \end{equation} $$

用 $-\mathcal H(\mathbf u)$ 表示 $\mathbf A \mathbf u$ 的非对角成分，即

$$ \begin{equation} -\mathcal H(\mathbf u) = \mathbf A \mathbf u - \mathbf D \mathbf u \end{equation} $$

其中 $\mathbf D$ 为系数矩阵 $\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 \label{eq:mom_eqn} \end{equation} $$

对上式进行算子分裂，引入上标 $^{\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}$ 不满足连续性方程，需要对其进行修正。

2. 第一次修正

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

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

则动量方程 $\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 \label{eq:1st_mom_eqn} \end{equation} $$

这里 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 \label{eq:ppe} \end{equation} $$

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

3. 第二次修正

由于第一次修正步中的 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 \label{eq:2nd_mom_eqn} \end{equation} $$

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

$$ \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/fvSolution 中 PISO 或 PIMPLE 子字典的关键词 momentomPredictor 控制。参考 icoFoam 求解器中的代码：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

        // 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$ 操作符中的速度直接用上一时刻或上一迭代步的值⁴。

momentumPredictor 的默认值是 true。但是，OpenFOAM 随两相流求解器 interFoam 提供的绝大部分示例算例中都关闭了动量预估，对此，Henry 的解释是⁵：

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 认为如果动量方程中的扩散项很重要，则需要求解动量方程⁶：

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 操作符中的系数矩阵

在 Issa 提出的 PISO 算法中，$\mathcal H$ 操作符只包含对流项和扩散项的离散成分，而 OpenFOAM 中 $\mathcal H$ 操作符除了包含对流项和扩散项之外，还计入了时间项的影响。求解器在构建 UEqn 时，把 fvm::ddt 时间项也放到了其中。

而 foam-extend 的实现与原始版本的 PISO 一致，在计算 $\mathcal H$ 操作符时并没有考虑时间项。以 foam-extend 4.1 中的 icoFoam 为例：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

        // 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 算法的范畴，关于这两个话题我们以后有机会再详聊。

参考