版本对应：
OpenFOAM-11中的incompressibleFluid模块，对应OpenFOAM-10之前的simpleFoam与pimpleFoam。OpenFOAM-10中的稳态瞬态求解器单独处理，OpenFOAM-11中都整合在了incompressibleFluid模块。
CFD: 不可压 + 瞬态
Issa在1986年提出的PISO算法，最重要的特点在于其为一个非迭代的瞬态压力速度耦合算法。PISO算法适用于不可压缩流以及可压缩流。OpenFOAM中的瞬态PISO算法与原始PISO算法大同小异。本文讨论的算法主要取自OpenFOAM-11中的incompressibleFluid模块。需要注意的是，incompressibleFluid模块可以同时处理瞬态、稳态，本文仅仅关注瞬态算法。
在老版本的OpenFOAM中存在3个瞬态不可压缩单相流求解器：icoFoam、pimpleFoam、pisoFoam。icoFoam为一个demo式的求解器，主要便于初学者理解代码，其不可调用湍流模型。pisoFoam为一个纯正的瞬态不可压缩单相流求解器。pimpleFoam可以调用pimple算法（将在后续介绍）。在OpenFOAM-11中，三者不再进行区分。incompressibleFluid模块可以涵盖所有的瞬态不可压缩压力速度耦合算法。
另外需要注意的是，OpenFOAM-11中的incompressibleFluid模块可以进行DNS模拟。其可用于模拟层流流动或流场剪切力引致的湍流准直接模拟 (quasi-DNS)，也被称之为准DNS。在这里需要注意的是，得益于谱方法以及有限差分法高阶精度的易用性，DNS计算通常采用谱方法或有限差分法。有限体积法由于本身在计算通量时引入的中点插值法则，大大限制了有限体积法的精度。然而Komen等通过研究圆管内的湍流流动，认为在网格分辨率足够的情况下，使用有限体积法同样可以获得高质量的平均速度以及脉动速度流场 。Komen等进一步的使用有限体积法进行准DNS计算，并与有限体积法大涡模拟进行标定，认为大涡模拟方法会引致一定的数值耗散。Axtmann and Rist也使用OpenFOAM进行了准DNS直接模拟，并研究了编译器与并行计算的特性。Liu等通过icoFoam对开口顶盖驱动流进行了准直接模拟 ，进一步验证了使用OpenFOAM进行准直接模拟的可行性。
控制方程与算法
给定瞬态的N-S方程：
在这里我们不考虑湍流粘度的影响。首先对方程(2)中的时间项进行对速度\(\bfU\)关于时间\(t\)的欧拉全隐离散有 ：
对流项隐性离散：
拉普拉斯项隐性离散：
压力项显性离散：
其中上标\(^t\)表示为当前的时间步（已知），上标\(^*\)表示预测步（待求），下标\(_f\)表示网格单元面上的值，\(V_\rP\)表示网格单元体积，\(\bfS_f\)表示网格单元的各个面的面矢量，\(\Delta t\)表示时间步长，\(F_f^t\)为通量，\(\nu\)为动力粘度。方程(5)中\((\nabla\mathbf{U}^*)_f\cdot\bfS_f \)需要进一步处理，其中\(\nabla\mathbf{U}^*\)为定义在体心的量，即网格体心的速度梯度。\((\nabla\mathbf{U}^*)_f\)为定义在面心的量，即网格面心的速度梯度。
需要注意的是，方程(2)中的左边第二项（对流项）是非线性的。在求解的过程中，要么选择非线性求解器，要么将对流项线性化。由于非线性求解器非常复杂，因此在OpenFOAM均采用线性化处理。具体的，在对方程(7)中的左边第二项对流项离散中，其中的通量\(F_f^t\)采用当前已知时间步的速度\(\mathbf{U}^t\)来计算，同时保留\(\mathbf{U}^*_\mathrm{P}\)为未知量。这种将高阶非线性项降为一阶线性项的过程即为线性化操作。线性化的一个问题即为通量速度的信息有所滞后。
方程(7)中\(\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}\)表示当前网格单元的速度。速度的面法向梯度在三维情况下为一个矢量（变量的面法向梯度与变量的阶数相同）。在这里需要注意的是，方程(11)只对矩形网格精准，网格非正交性增加，需要特定的数值处理提高其计算精准性。
注意：
紧致基架点与拓展基架点的英文被称之为compact stencil与extended stencil，有关基架点的详细描述可参考无痛苦NS方程笔记中有关laplacian(ϕ) = div(grad(ϕ))?的讨论。
将上式简写为
其中
注意：
如果单单看扩散项的系数，有\(A_\mathrm{P}^t=\frac{1}{V_\mathrm{P}}\sum { {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} }\)，\(-\sum A_\mathrm{N}^t=\sum \frac{1}{V_\mathrm{P}} { {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} }\)。很明显\(A_\mathrm{P}^t=-\sum A_\mathrm{N}^t\)。此对应关系大量的出现在CFD教材中。对流项也有类似的特征，只不过不是相等的关系，而是相减的关系。
注意：
另外，上面的离散都可以被认为是内部面的离散，对于边界面，可能会导致不同分量离散出来的矩阵系数可能不同。单OpenFOAM里面的UEqn.A()只有一个。这方面的介绍请参考 这个链接。边界面的离散不影响理解整个算法流程。
PISO算法
可以看出在某个时间步内\(A_\mathrm{P}^t\)、\(A_\mathrm{N}^t\)、和\(S_\mathrm{P}^t\)保持不变（注意\(A_\mathrm{N}^t\)、\(\sum A_\mathrm{N}^t\bfU_\rN\)、\(-\sum A_\mathrm{N}^t\bfU_\rN+S_\rP^t\)的区别）。求解方程(14)即可获得预测速度\(\mathbf{U}^*_\mathrm{P}\)。方程(14)即为OpenFOAM中的动量预测方程。需要注意的是，动量预测步骤并不是必须的，这主要和速度压力耦合求解策略有关。如果不调用动量预测步骤，则\(\mathbf{U}^*_\mathrm{P}=\mathbf{U}^t_\mathrm{P}\)。
N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度，但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来，其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下，对连续性方程\(\nabla\cdot\bfU=0\)进行离散后的形式为
如果能用压力表示方程(18)中的\(\mathbf{U}_{\mathrm{P},f}^{t+\dt}\)，是不是就可以得出压力泊松方程了呢？答案是肯定的。首先，方程(14)可以写为：
在当前时间步收敛情况下为：
方程(20)为一个非线性方程组。使其线性化后可以写为：
方程(20)相对于方程(21)存在一定的滞后。在\(\dt\)较小的情况下这被认为是较小的。若将方程(21)减去(19)有：
其中\(\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{t+\dt}-\mathbf{U}_\mathrm{P}^{*}\)（其他同理）。在这里未调用任何略去临点的假定。对方程(21)移项有
类似的，面上的速度可以通过下面的公式获得：
方程(26)即下述方程的离散形式：
方程(27)即为压力泊松方程。若有\(\mathbf{HbyA}^{t+\dt}\)，则可求得收敛的压力。需要注意的是，\(\nabla \cdot\left(\frac{1}{A^t} \nabla p^{t+\dt}\right)\)可以通过先求梯度再求散度进行，也可以直接通过拉普拉斯积分进行。方程(26)右端即先求梯度再求散度的方法，这种方法会引起数值震荡。若通过拉普拉斯并进行高斯积分进行离散，则会调用紧致格式防止数值震荡 (Li, 2019)。
可以看出，最后求解的速度\(\mathbf{U}_\mathrm{P}^{t+\dt}\)和压力\(p^{t+\dt}\)应该符合方程(14)和(27)，分别对应动量方程和连续性方程。然而目前，我们只有通过动量预测步骤求出来的\(\mathbf{U}_\mathrm{P}^*\)。\(\mathbf{HbyA}_\mathrm{P}^{t+\dt}\)和\(p^{t+\dt}\)均为未知的。压力泊松方程不能求解。1986年Issa提出了PISO算法可以解决这个问题 (Issa, 1986)。PISO算法在提出时主要针对不可压缩非稳态计算，其为一种在时间步内非迭代的算法，在随后也被拓展用于稳态问题以及可压问题中。PISO算法通过对当前时间步内的多次修正来获得最终的\(\mathbf{U}_\mathrm{P}^{t+\dt}\)和\(p^{t+\dt}\)。依据PISO算法，在方程(22)的基础上引入略去临点影响的假定有：
其中\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^{*}+\mathbf{U}_\mathrm{P}'\)，\(p_f^{*}=p_f^{t}+p_f'\)。相比方程(21)，方程(29)由于忽略了临点的影响并不是精准的。
注意：
如果不进行动量预测步计算，则\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^{n}+\mathbf{U}_\mathrm{P}'\)。
对方程(29)进行移项有
方程(32)中可参考方程(25)-(27)的步骤组建压力泊松方程，即
对方程(27)求解后有\(p^{*}\)。将\(p^*\)回代到方程(30)的时候，有\(\mathbf{U}_{\mathrm{P}}^{**}\)。这里的\(\mathbf{U}_{\mathrm{P}}^{**}\)满足连续性方程。但由于方程(29)忽略了临点的影响，\(\mathbf{U}_{\mathrm{P}}^{**}\)与\(p^{*}\)并不严格满足动量方程。因此在实际操作中，需要以计算后的\(\mathbf{U}_{\mathrm{P}}^{**}\)继续组建\(\bfHbyA^{**}\)，然后继续求解压力方程获得新的压力场\(p^{**}\)，然后进一步获得新的速度场\(\mathbf{U}_{\mathrm{P}}^{***}\)。这即为时间步内的第二次迭代。当然也可以在时间步内继续进行迭代。
因此对于PISO算法，每个时间步内的迭代可以表示为：
依据初始条件求解预测速度\(\mathbf{U}^*\)。通过方程(32)组建\(\mathbf{HbyA}^*\)。
求解方程(33)求解压力\(p^*\)。
通过方程(32)依据压力\(p^*\)以及\(\mathbf{HbyA}^*\)更新速度有\(\mathbf{U}^{**}\)。
回到第二步迭代求解几次。
PISO算法在时间步内的迭代一般只需要3步即可。相对于稳态SIMPLE算法的若干次迭代，瞬态PISO算法的区别主要在于：
瞬态PISO算法的时间步长要保证库朗数足够小，因此每一步的流场变动相对来说并不是很激进。因此，在求解方程(33)与方程(29)的时候，可以将残差设定为非常小的值。保证动量方程与连续性方程都是收敛的。这样，误差仅仅来自于略去临点的误差。大大增加每个时间步内的收敛速度。
瞬态PISO算法不需要对压力场进行松弛。但可以对速度、湍流变量的矩阵进行少量的松弛，过量的松弛会导致时间步进不精准。因此，PISO算法时间步内的算法可以理解为没有松弛的SIMPLE算法。
理论上，瞬态PISO算法也可以与SIMPLEC相结合，形成单时间步内迭代更快的算法。但由于瞬态求解器的库朗数的严格限制。使用SIMPLE算法就可以在每时间步内保证3次以内收敛，因此实践上，并没有使用SIMPLEC算法的实际意义。
PIMPLE算法
PIMPLE算法是OpenFOAM基金会提出特有的计算方法。并没有见过在文献里面当做专门的算法而发表的文章。PIMPLE算法更像是一种算法的融合，而非重新开发。PIMPLE算法有一些非常吸引人的特性：
调用库郎数远大于1的时间步长且同时进行瞬态步进；
在时间精度不重要的时候，获得一种拟瞬态的稳态解；
PIMPLE算法将PISO算法与SIMPLE算法融合，获得大时间步长步进的计算能力。在瞬态PISO算法中（例如上文讨论的icoFoam），需要保证时间步长要保证库朗数足够小，这样才能不需要加松弛，且保证每个时间步内计算的结果是收敛的。在PIMPLE算法中，可以调用很大的时间步长。在这种情况下，每个时间步长内的流场或许会变化的很激进容易发散，因此在PIMPLE算法中，可以在每个时间步内添加松弛。可以把PIMPLE理解为，在较大时间步长下，每个时间步内去添加松弛来获得收敛解，然后向后进行推进。PIMPLE算法将低松弛去掉，并关闭PIMPLE循环，可以完全的演变为瞬态PISO算法。
若使用PIMPLE算法，建议使用比较高的时间精度格式来减少高库郎数导致的误差。然而在全局库郎数比较大的情况下，时间历史步进的不精准是难以避免的，同时松弛因子对结果的影响也会加大。因此计算结果如果需要非常高精度的时间步进历史，需要使用PISO算法。若对时间步进要求存在要求且网格大小分布非常不均一，可以使用PIMPLE算法增大库郎数限制，但要关注全局库郎数的变化。若全局库郎数的变化特别大。不可否认的会牺牲流场变化的时间精度。若对时间步进的历史过程不关心，仅仅关注最后稳态的情况，且同时算例存在收敛解。完全可以使用SIMPLE算法来进行。也可以通过PIMPLE算法来增加到足够的库郎数来进行时间推进。也可以通过局部时间步Local Time Step的方式来进行。
下面是一个典型的PIMPLE设置。其中的nOuterCorrectors用来控制PIMPLE循环次数，也被称之为外循环。表示在每个时间步长内进行100次迭代（可以理解为100次SIMPLE迭代）。nCorrectors表示内循环。表示每次外循环内的压力方程求解次数，可以理解为PISO的迭代次数。在PIMPLE算法中，一般内循环设置为1次即可，即不需要多次PISO内循环。
因此，nOuterCorrectors设置为1的时候，PIMPLE完全的变为了PISO算法。outerCorrectorResidualControl用来控制PIMPLE外循环的收敛标准（类似SIMPLE算法的收敛标准）。比如本算例设置100次外循环，但如果满足相应的收敛标准，则停止迭代。对于一个存在近期时间稳态的算例，在求解后期，PIMPLE外循环可以降低为10次以下甚至更少，且保持非常高的库郎数递进。relaxationFactors用来设置外循环的松弛因子。其设置完全可以参考SIMPLE算法。
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 100;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
outerCorrectorResidualControl
{
"(h|p)"
{
relTol 0;
tolerance 0.0001;
}
}
}
relaxationFactors
{
fields
{
p 0.3;
pFinal 1;
}
equations
{
"U|k|epsilon|h" 0.7;
"(U|k|epsilon|h)Final" 1;
}
}
关键代码
incompressibleFluid模块中的速度方程通过下述代码来进行组建：
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
其中fvm::ddt(U)对应方程(3)，这里的fvm表示隐形离散，反之，若写为fvc则表示显性离散。fvm::div(phi, U)对应方程(4)，注意其中的phi，在OpenFOAM中，phi大部分情况下表示为通量。进一步的，在不可压缩流体中，phi表示定义在面上的体积通量\(\bfU_f\cdot\bfS_f\)，在可压缩流体中，phi表示定义在面上的质量通量\(\rho_f\bfU_f\cdot\bfS_f\)。fvm::laplacian(nu, U)则对应方程(5)。
方程(15)中定义的\(A_\rP^t\)为一个volScalarField，其可以通过UEqn.A()来表示。\(A_\rP\)的倒数，在代码中被表示为rAU，其代码如下：
volScalarField rAU(1.0/UEqn.A());
在方程(23)中，定义了\(\mathbf{HbyA}\)，其可以理解为进行一次迭代后获得的一个速度，更激进的，可以把HbyA理解为速度中间量。方程(23)右侧括号内的\({ - \sum {{A_\mathrm{N}^t}\mathbf{U}_\mathrm{N}^{t+\dt}} + S_\mathrm{P}^t}\)可以用代码UEqn.H()来表示。因此进一步的，\(\mathbf{HbyA}\)可以表示为
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
其中第一行仅仅为声明场，第二行为对其进行赋值。在获得\(\mathbf{HbyA}\)之后，方程(27)可以用来求解压力，相应的代码如下：
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
如果进一步的看foamRun代码，
{
solver.moveMesh();//网格变形
solver.fvModels().correct();//源项添加
solver.prePredictor();//关键步1
solver.momentumPredictor();//关键步2
solver.thermophysicalPredictor();//关键步3
solver.pressureCorrector();//关键步4
solver.postCorrector();//后处理
}
其中的：
关键步1不执行；
关键步2初始速度矩阵，若需要动量预测，则进行动量预测；
关键步3不执行；
关键步4进行压力求解；