PISO，全称为 Pressure-Implicit with Splitting of Operators，是 Issa1 于 1986 年提出的一种求解耦合速度压力的非迭代算法，也是目前应用最广的速度压力解耦方法之一。OpenFOAM 中实现的 PISO 算法与原始版本并不完全相同，本文将对原始的 PISO 算法及 OpenFOAM 中的 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 用到的算子分裂是代数算子分裂，其形式如下

$$$\frac{\partial \phi}{\partial t}+L(\phi) = 0, \quad L = \sum_{s=1}^S L_s$$$

$$$\frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf u) - \nabla \cdot (\gamma \nabla \phi) = 0$$$

$$$\phi^{n+1} = \mathcal F(\phi^n, \Delta t)$$$

$$$\frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf u) = 0$$$

$$$\frac{\partial \phi}{\partial t} - \nabla \cdot (\gamma \nabla \phi) = 0$$$

$$$\frac{\partial \phi}{\partial t} + \mathcal C(\phi) = 0$$$

$$$\frac{\partial \phi}{\partial t} + \mathcal D(\phi) = 0$$$

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

### 动量方程的算子分裂

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

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

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

$$$\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 \mathbf u - \mathbf D \mathbf u$$$

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

1. 动量预估

$$$\frac{\mathbf u^{\ast} - \mathbf u^n}{\Delta t} - \mathcal H(\mathbf u^{\ast}) = - \frac{1}{\rho} \nabla p^{n} + \mathbf S$$$

1. 第一次修正

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

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

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

1. 第二次修正

$$$\nabla \mathbf u^{\ast\ast\ast} = 0$$$

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

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

## OpenFOAM 中的 PISO

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

### 动量预估可选

OpenFOAM 中的 PISO 算法的动量预估步骤是可选的，通过字典文件 constant/fvSolutionPISOPIMPLE 子字典的关键词 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 的默认值是 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 操作符中的系数矩阵

  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(); 

## 参考

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. http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture11.pdf [return]
3. Issa, R. I., Gosman, A. D., & Watkins, A. P. (1986). The computation of compressible and incompressible recirculating flows by a non-iterative implicit scheme. Journal of Computational Physics, 62(1), 66–82. [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]