fvm::laplacian(rUA, p) == fvc::div(phi)
这一行程序是 icoFoam 中的一行，写成张量形式为
\frac{1}{A_p}\nabla^2 p =\frac{1}{A_p}\nabla\cdot(\nabla p)= \nabla \cdot U
\tag{1}
称为压力泊松方程。rUA 就是方程中的 \dfrac{1}{A_p} 。
A_p 是动量方程离散后，去掉压力梯度项，当前 Cell 的矩阵系数：
A_\mathrm{P} \mathbf U_\mathrm{P}^{n+1} + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^{n+1} - E_\mathrm{P}^{n+1} = 0
\tag{2}
不过，我们先来看一下 U 方程：
下面 icoFoam.C 片段：
fvVectorMatrix UEqn // 速度方程
(
fvm::ddt(U) // 非稳态项，隐式离散
+ fvm::div(phi, U) // 对流项，隐式离散
- fvm::laplacian(nu, U) // 扩散项，隐式离散
);
此为无压力梯度项的速度方程，张量形式为
\text{UEqn} = \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{U})- \nabla \cdot(\nu \nabla \mathbf{U})
\tag{3}
该方程在后台就形成离散形式 (2)。
速度方程的张量形式为
\text{UEqn} =-\nabla \frac{p}{\rho}\tag{4}
求解速度方程，得预测速度 U^*
solve(UEqn == -fvc::grad(p)); // 不可压缩时, p=压力/密度
此时求得的预测速度 U^* 并不能满足连续性方程。因此，需要一点特殊的方法来处理。
假设修正后的速度为 \mathbf U^{n+1}
A_\mathrm{P} \mathbf U_\mathrm{P}^{n+1} + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^{n+1} - E_\mathrm{P}^{n+1} = - \nabla p^{n+1}
\tag{5}
修正前的速度 (去掉压力梯度项) 为
\mathbf{U}^* = \mathbf U_\mathrm{P}^r = \frac{1}{A_\mathrm{P}}\left( - \sum A_\mathrm{N} \mathbf U_\mathrm{N}^r + E_\mathrm{P}^r \right)
\tag{6}
方程 (5)-A_\mathrm{P}(6) 为
\mathbf U_\mathrm{P}^{n+1} -\mathbf{U}^{*,n+1}= - \frac{1}{A_\mathrm{P}}\nabla p^{n+1}
\tag{7}
稍作变换
\mathbf U_\mathrm{P}^{n+1} = \mathbf{U}^{*,n+1} - \frac{1}{A_\mathrm{P}}\nabla p^{n+1}
\tag{7a}
面插值版本，
\mathbf U_f^{n+1} = \mathbf{U}_f^{*,n+1} - \frac{1}{A_\mathrm{P}}\nabla_{\! f} \,p^{n+1}
\tag{7b}
而速度 \mathbf U_\mathrm{P}^{n+1} 符合连续性方程
\nabla \cdot \mathbf U_\mathrm{P}^{n+1}=0
\tag{8}
将 (8) 代入 (7) 的散度，有
\nabla\cdot\mathbf{U}^{*,n+1}= \nabla\cdot\left(\frac{1}{A_\mathrm{P}}\nabla p^{n+1}\right)
\tag{9}
这就是压力泊松方程。将 (9) 离散，得
\sum \left[ \left( \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}}\nabla_f p^{n+1} \right) \cdot \mathbf S_f \right] = 0
\tag{10}
求解该方程将用到修正前的速度 U^* ，从而可求解 p ，然后将 p 用于修正速度场 (方程 (7a))。
以下为压力泊松方程的循环。
... // 前面已经完成预测速度 Ｕ* 的计算
for (int corr = 0; corr < nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A(); // 取得 U 方程的 1/Ap
U = rUA*UEqn.H(); // 无压力梯度项之速度值 U'
phi = (fvc::interpolate(U) & mesh.Sf()) // U' 使用线性插值，需要在算例中设定interpolate(U)的格式
+ fvc::ddtPhiCorr(rUA, U, phi); // 在 ext-4.0 版本中，这一项已经去掉
adjustPhi(phi, U, p); // 修正压力边界（梯度为0）上的 phi，使其满足连续性方程
//网格非正交循环
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi) // 压力泊松方程，此处间接使用 Rhie-Chow 插值原理
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(); // 求解压力泊松方程
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux(); // 方程 (7b)，使用求解的压力场修正通量场，即在最后一次修正的时候使通量守恒，
}// pEqn.flux()返回方程 (7) RHS，为 fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf())
} // 某些可压缩求解器中, 修正项为 phi += pEqn.flux(), 这是因为 pEqn 中的 laplacian 项为 '-' 号
# include "continuityErrs.H"
U -= rUA*fvc::grad(p); // 公式 (7a)
U.correctBoundaryConditions();
}
icoFoam 使用 PISO 算法解 N-S 方程。 其基本程序为
- 根据初始条件（压力场、速度场等）求解预测速度场 U^*
- 根据预测速度求无压力梯度项的速度场 U'
- 根据无压力梯度项的速度场 U' ，求解压力（泊松）方程，得到压力场 p
- 根据压力场 p ，修正预测速度 U^* 以满足连续性方程
- 接第 2 步，直到满足收敛要求