版本对应：
OpenFOAM-11中的incompressibleFluid模块，对应OpenFOAM-10之前的simpleFoam与pimpleFoam。OpenFOAM-10中的稳态瞬态求解器单独处理，OpenFOAM-11中都整合在了incompressibleFluid模块。
CFD: 不可压 + 稳态
Patankar and Spalding开发的SIMPLE算法主要用于处理稳态的不可压缩的压力速度耦合。本文主要介绍OpenFOAM中的稳态压力速度耦合算法。其具体的求解细节不同于原生的SIMPLE算法，但大体相同。算法主要取自OpenFOAM-11中的incompressibleFluid模块。需要注意的是，incompressibleFluid模块可以同时处理瞬态、稳态，本文仅仅关注稳态算法。
OpenFOAM-11中的incompressibleFluid模块中，压力和速度的耦合通过SIMPLE/SIMPLEC/PISO算法进行计算，且没有考虑其他体积力（如重力等）。对于流场内温度变化较大，但可忽略浮力以及重力的可压缩稳态求解器，可选用isoThermalFluid模块（OpenFOAM-10版本之前的rhoSimpleFoam）。附加重力的可压缩浮力驱动流稳态求解器可选用fluid模块（OpenFOAM-10版本之前的buoyantSimpleFoam）。本文从最基本的N-S方程在笛卡尔网格下的有限体积离散开始推导，介绍OpenFOAM中调用的SIMPLE/SIMPLEC算法。
控制方程与算法
首先有稳态的N-S方程：
首先，需要对方程(2)做体积积分并隐性离散：
其中上标\(^n\)表示为当前迭代步（已知），上标\(^*\)表示预测步（待求），下标\(_f\)表示网格单元面上的值，\(V_\rP\)表示网格单元体积，\(\bfS_f\)表示网格单元的各个面的面矢量，\(F_f\)为通量，\(\nu\)为动力粘度。依据方程(3)-(5)有：
其中的通量\(F_f^n\)采用当前已知时间步的速度\(\mathbf{U}^n\)来计算，同时保留预测速度\(\mathbf{U}^*\)为未知量。
方程(6)中\((\nabla\mathbf{U}^*)_f\cdot\bfS_f\)需要进一步处理，其中\(\nabla\mathbf{U}^*\)为定义在体心的量，即网格体心的速度梯度。\((\nabla\mathbf{U}^*)_f\)为定义在面心的量，即网格面心的速度梯度。为使用紧致基架点防止数值震荡，\((\nabla\mathbf{U}^*)_f\cdot\bfS_f\)通常表示为面法向梯度与面矢量的模的乘积：
其中\((\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} \)表示面法向梯度，有
其中\(\mathbf{d}\)表示网格单元体心之间矢量差（距离矢量），下标\(_\mathrm{N}\)表示相邻网格单元的速度，下标\(_\mathrm{P}\)表示当前网格单元的速度。在这里需要注意的是，方程(8)只对矩形网格精准，网格非正交性增加，需要特定的数值处理提高其计算精准性。
方程(6)中\(\mathbf{U}_f^*\)被定义为面上的预测速度，然而在计算中求得的均为网格单元体心上的速度。因此需要从体心速度进行插值来获得面速度，在此步可以引入各种插值格式。假设使用中心线性格式：
方程两边同时除以网格体积，并整理有：
将上式简写为
其中
求解方程(13)即可获得预测速度\(\mathbf{U}^*_\mathrm{P}\)。方程(13)即为OpenFOAM中的速度方程。也即动量预测步。
注意：
OpenFOAM中的动量预测步主要是为了组建稀疏线性系统，如方程(13)。目的是获得矩阵系数\(A\)。因此并不强制进行求解，如果求解的话，可获得\(\bfU^*\)，不求解的话，即为\(\bfU^n\)。
N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度，但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来，其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下，对连续性方程\(\nabla\cdot\bfU=0\)进行离散后的形式为
如果能用压力表示方程(16)中的\(\mathbf{U}_{\mathrm{P},f}^{n+1}\)，是不是就可以得出压力泊松方程了呢？答案是肯定的。首先，方程(13)中的压力项不建议进行离散，这是因为方程右侧的梯度项都有可能引起数值的震荡（例如双流体模型中的湍流分散力也需要特殊处理）。同时，参考方程(13)， 在收敛情况下（时间步\(t=n+1\)）为：
由于存在未知的\(A_\mathrm{P}^{n+1},A_\mathrm{N}^{n+1}\)。因为后续将引入迭代算法，因此可以\(A_\mathrm{P}^{n+1},A_\mathrm{N}^{n+1}\)存在一定的滞后，即\(A_\mathrm{P}^{n+1}\approx A_\mathrm{P}^{n},A_\mathrm{N}^{n+1}\approx A_\mathrm{N}^{n}\)：
注意：
在这里未调用任何略去临点的假定。
其中\(\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U}_\mathrm{P}^{*}\)，\(\mathbf{U}_\mathrm{N}'=\mathbf{U}_\mathrm{N}^{n+1}-\mathbf{U}_\mathrm{N}^{*}\)，\(p_f'=p_f^{n+1}-p_f^{n}\)。既压力的修正量和当前网格点的速度修正量以及相邻网格点的速度修正量有关。对方程(18)移项有
其中我们定义
面上的速度可以通过下面的公式获得：
方程(24)即下述方程的离散形式：
方程(25)即为压力泊松方程。若有\(\mathbf{HbyA}^{n+1}\)，则可求得收敛的压力。
注意：
这里\(\mathbf{HbyA}^{n+1}\)是未知的，因此方程(25)一个方程存在两个未知量，不可解。因此需要引入下面的迭代算法。
SIMPLE算法
可以看出，最后求解的速度\(\mathbf{U}_\mathrm{P}^{n+1}\)和压力\(p^{n+1}\)应该符合方程(13)和(25)，分别对应动量方程和连续性方程。然而目前，我们只有通过动量预测步骤求出来的\(\mathbf{U}_\mathrm{P}^*\)。\(\mathbf{HbyA}_\mathrm{P}^{n+1}\)和\(p^{n+1}\)均为未知的。方程(25)不能求解。SIMPLE迭代算法则可以解决这个问题。依据SIMPLE算法，首先引入略去临点影响的假定，方程(19)可以写为：
其中\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^*+\mathbf{U}_\mathrm{P}'\)，\(p_f^{*}=p_f^n+p_f'\)。
注意：
如果不进行动量预测步计算，则\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^n+\mathbf{U}_\mathrm{P}'\)。
相比方程(18)，方程(27)由于忽略了临点的影响并不是精准的。对方程(27)进行移项有
方程(30)中可参考方程(23)-(25)的步骤组建压力泊松方程，即
注意：
如果不进行动量预测步计算，则\(\mathbf{HbyA}^{*}=\mathbf{HbyA}^{n}\)。
对方程(31)求解后有\(p^{*}\)。注将\(p^{*}\)与\(\mathbf{HbyA}^{*}\)回代到方程(28)的时候，有更新后的\(\mathbf{U}_{\mathrm{P}}^{**}\)。此处的\(\mathbf{U}_{\mathrm{P}}^{**}\)满足连续性方程。同时，\(\mathbf{U}_{\mathrm{P}}^{**}\)与\(p^*\)近似满足动量方程(28)。之所以近似，是因为方程(28)在推导的过程中，调用了忽略临点假定。
注意：
如果不进行动量预测步计算，则完成一次迭代后获得\(\bfU^{*}\)与\(p^{*}\)。
至此，原变量\(\bfU^n\)与\(p^n\)更新为\(\bfU^{**}\)与\(p^{*}\)。此处即完成一次迭代。我们强调，\(\bfU^{**}\)与\(p^{*}\)近似满足动量方程(28)，理论上满足连续性方程。在这里，需要强调以下迭代算法的重要性。
方程(27)略掉了临点的影响，因此只能近似满足动量方程。但最终收敛的情况下，修正量会变得非常小，因此方程(27)同样会满足。
在实际计算过程中，尤其是最开始的迭代时，压力修正量过大，会导致难以收敛。因此，在\(p^n\)变化至\(p^*\)的过程中，要减少这种非常大的突变，要对压力进行亚松弛。这就导致方程(31)求解后的压力\(p\)，为介于\(p^n\)变化至\(p^*\)之间的某个值。
同时，方程(31)的求解并不会将残差将至最低，也就是方程(31)求解后左右两边并不是严格相等，这也就意味着存在较大的连续性误差，\(\bfU^{**}\)并不严格的满足连续性方程。
因此，通过算例的残差、松弛必要设置，导致\(\bfU^{**}\)与\(p^{*}\)不能作为最终收敛的结果。因此需要将\(\bfU^{**}\)与\(p^{*}\)作为下一个迭代步的初始值，继续进行缓慢稳定的推进。因此，SIMPLE算法中的的迭代过程可以表示为下面几个步骤：
SIMPLEC算法
SIMPLEC算法认为SIMPLE的算法太过激进。这种略去邻点的行为会使得压力修正值\(p'\)过大，在某些情况下会引起发散。虽然这可以通过亚松弛来稳定计算流程，但是会减慢收敛速度。SIMPLEC算法认为任意网格点的速度修正可以近似认为是相邻网格点速度修正的平均：
即
相对比于SIMPLE算法中的方程(26)，对于固定的\(\bfU_\rP'\)，方程(34)右侧压力修正的乘积项\(\frac{1}{A_\rP+\sum A^n_\rN}\)必然大于\(\frac{1}{A^n_\rP}\)（有限体积法离散的\(A^n_\rN<0\)），进而\(p_f'\)的值不会特别大，因此SIMPLEC算法并不需要对压力进行激进的低松弛修正。将方程(34)和方程(13)相加有：
既
令\(\nabla\cdot\bfU^{**}=0\)有：
在SIMPLEC算法中，
求解后有\(p^*\)。随后\(\bfU^{**}\)可以下式获得：
同样的，\(\bfU^{**}\)与\(p^{*}\)近似满足动量方程，理论上满足连续性方程。但在这里可以强调SIMPLEC与SIMPLE差别：SIMPLE略去临点的假定特别激进，导致压力修正量过高。SIMPLEC相对温和，因此压力修正量并不会太大，因此SIMPLEC不需要对压力进行场松弛。同样的，每个迭代步求解的\(\bfU^{**}\)与\(p^{*}\)更加满足动量方程。
关键代码
稳态压力速度耦合求解的特殊性在于需要进行松弛。松弛又分为场松弛，以及方程松弛。例如对于速度方程，incompressibleFluid模块中momentumPredictor()语句要对速度进行矩阵松弛：
UEqn().relax();
力图提高速度场矩阵的对角占有特性。对方程进行松弛并不会改变求解方程后的解，但亚松弛会增加迭代次数。
同时，incompressibleFluid模块中correctPressure()语句要对压力进行场松弛：
p.relax();
场松弛会改变压力的值，亚松弛会改变压力的变化量。simpleFoam中可以通过下述代码调用SIMPLEC算法：
if (pimple.consistent())
{
rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
phiHbyA +=
fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU())*fvc::grad(p);
}
其中\(1/(A_\rP+\sum A_\rN)\)通过下列代码定义：
rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
方程(38)通过下列代码定义：
HbyA -= (rAU - rAtU())*fvc::grad(p);
同时，在上文的讨论中，稳态算法中动量方程(28)以及连续性方程(31)的求解并不需要将残差将至最低，因此在fvSolution中进行类似的设置：
p
{
solver GAMG;
relTol 0.1;//残差降低一个数量级即可
smoother GaussSeidel;
}
"(U|k|epsilon|omega|f|v2)"
{
solver smoothSolver;
smoother symGaussSeidel;
relTol 0.1;//残差降低一个数量级即可
}
如果进一步的看foamRun代码，
{
solver.moveMesh();//网格变形
solver.fvModels().correct();//源项添加
solver.prePredictor();//关键步1
solver.momentumPredictor();//关键步2
solver.thermophysicalPredictor();//关键步3
solver.pressureCorrector();//关键步4
solver.postCorrector();//后处理
}
其中的：
关键步1不执行；
关键步2初始速度矩阵，若需要动量预测，则进行动量预测；
关键步3不执行；
关键步4进行压力求解；