为了求解Navier-Stokes方程，往往需要解耦速度场和压力场的耦合方程组。 OpenFOAM 的求解器大多使用 PISO 或者 SIMPLE 算法，或者二者的结合体 PIMPLE 算 法。这些算法一般采用预测-校正策略，通过迭代的方式将速度和压力的计算解耦。

PISO 以及 PIMPLE 算法用来处理非稳态问题，SIMPLE 算法用来处理稳态问题。

1.SIMPLE 算法

上图为SIMPLE 算法的流程图，其求解步骤一般如下:

1.检查是否达到收敛-simple.loop()
2.使用动量预测器-UEqn.H 预测速度
3.校正压力和速度–pEqn.H
4.求解湍流模型的传输方程–湍流->校正（）
5.返回步骤1

主程序：

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "simpleControl.H"
#include "fvOptions.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //int main(int argc, char *argv[])
{#include "postProcess.H"#include "setRootCase.H"#include "createTime.H"#include "createMesh.H"#include "createControl.H"#include "createFields.H"#include "createFvOptions.H"#include "initContinuityErrs.H"//进行湍流模式判定turbulence->validate();// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Info<< "\nStarting time loop\n" << endl;while (simple.loop()){Info<< "Time = " << runTime.timeName() << nl << endl;// --- Pressure-velocity SIMPLE corrector{#include "UEqn.H"#include "pEqn.H"}laminarTransport.correct();turbulence->correct();runTime.write();Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"<< "  ClockTime = " << runTime.elapsedClockTime() << " s"<< nl << endl;}Info<< "End\n" << endl;return 0;
}

1.动量预测：（该代码来自`simpleFoam 求解器）

MRF.correctBoundaryVelocity(U);fvVectorMatrix UEqn
(fvm::ddt(U) + fvm::div(phi, U)+ MRF.DDt(U)+ turbulence->divDevReff(U)==fvOptions(U)
);UEqn.relax();fvOptions.constrain(UEqn);if (piso.momentumPredictor())
{solve(UEqn == -fvc::grad(p));fvOptions.correct(U);
}

2.校正

根据求得的预测速度校正压力场，再根据连续性方程，校正获得新的速度场。

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
("phiHbyA",fvc::flux(HbyA)+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);MRF.makeRelative(phiHbyA);adjustPhi(phiHbyA, U, p);// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU, MRF);// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{// Pressure correctorfvScalarMatrix pEqn(fvm::laplacian(rAU, p) == fvc::div(phiHbyA));pEqn.setReference(pRefCell, pRefValue);pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));if (piso.finalNonOrthogonalIter()){phi = phiHbyA - pEqn.flux();}
}#include "continuityErrs.H"U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);

2.PISO算法

PISO算法也采用预测—校正策略，并对速度场进行预测。然后根据事先设定的迭代次数进行压力场和速度场的修正。最后，求解湍流模型的传递函数。

其编程代码基本等同于SIMPLE算法，所不同的点在于SIMPLE算法解决的是稳态问题，因此不需重复地求解每个时间步的速度场。其主要差别主要体现在主程序：

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //int main(int argc, char *argv[])
{#include "postProcess.H"#include "setRootCase.H"#include "createTime.H"#include "createMesh.H"#include "createControl.H"#include "createFields.H"#include "createFvOptions.H"#include "initContinuityErrs.H"turbulence->validate();// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Info<< "\nStarting time loop\n" << endl;while (runTime.loop()){Info<< "Time = " << runTime.timeName() << nl << endl;#include "CourantNo.H"// Pressure-velocity PISO corrector{#include "UEqn.H"// --- PISO loopwhile (piso.correct()){#include "pEqn.H"}}laminarTransport.correct();turbulence->correct();runTime.write();Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"<< "  ClockTime = " << runTime.elapsedClockTime() << " s"<< nl << endl;}Info<< "End\n" << endl;return 0;
}

3.PIMPLE算法

PIMPLE算法作为前两者的结合体，它在每个时间步长内用SIMPLE稳态算法求解（也就是将每个时间步内看成稳态流动），时间步长的步进用PISO算法来完成。

##在有限容积离散中，时间项的离散仍采用的差分格式，这样做可以得到某个时间点的流场信息，而非某个时间步长的内的平均值。采用传统的piso算法求解变化较快的流动的时候，需要的时间步长较小（因为相邻两个时间点的流场不能差别太大，否则会发散），这样会造成求解的某种流动需要的耗时过长。 pimple算法将每个时间步长内看成一种稳态的流动（采用亚松驰来解决相邻两个时间段变化大的情况），当按照稳态的求解器求解到一定的时候，则采用标准的piso做最后一步求解。##(本段摘录自苏老师的博客，及参考内容1)

主程序代码:

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //int main(int argc, char *argv[])
{#include "postProcess.H"#include "setRootCase.H"#include "createTime.H"#include "createMesh.H"#include "createControl.H"#include "createTimeControls.H"#include "createFields.H"#include "createFvOptions.H"#include "initContinuityErrs.H"turbulence->validate();// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Info<< "\nStarting time loop\n" << endl;while (runTime.run()){#include "readTimeControls.H"#include "CourantNo.H"#include "setDeltaT.H" //设置时间步长runTime++;Info<< "Time = " << runTime.timeName() << nl << endl;// --- Pressure-velocity PIMPLE corrector loop//外循环，SIMPLE修正速度场while (pimple.loop())      {#include "UEqn.H"// --- Pressure corrector loopwhile (pimple.correct()){#include "pEqn.H"}if (pimple.turbCorr())  //湍流修正{laminarTransport.correct();turbulence->correct();}}runTime.write();Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"<< "  ClockTime = " << runTime.elapsedClockTime() << " s"<< nl << endl;}Info<< "End\n" << endl;return 0;
}

速度方程：

MRF.correctBoundaryVelocity(U);
//tmp瞬态对象模板类，详见(https://editor.csdn.net/md/?articleId=107122685)
tmp<fvVectorMatrix> tUEqn
(fvm::ddt(U) + fvm::div(phi, U)+ MRF.DDt(U)+ turbulence->divDevReff(U)==fvOptions(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();UEqn.relax();fvOptions.constrain(UEqn);if (pimple.momentumPredictor())
{solve(UEqn == -fvc::grad(p));fvOptions.correct(U);
}

压力方程：

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
("phiHbyA",fvc::flux(HbyA)+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);MRF.makeRelative(phiHbyA);adjustPhi(phiHbyA, U, p);tmp<volScalarField> rAtU(rAU);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);
}if (pimple.nCorrPISO() <= 1)
{tUEqn.clear();   //清空场对象
}// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);// Non-orthogonal pressure corrector loop
//内循环，PISO修正
while (pimple.correctNonOrthogonal())
{// Pressure correctorfvScalarMatrix pEqn(fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA));pEqn.setReference(pRefCell, pRefValue);pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));if (pimple.finalNonOrthogonalIter()){phi = phiHbyA - pEqn.flux();}
}#include "continuityErrs.H"// Explicitly relax pressure for momentum corrector
p.relax();U = HbyA - rAtU()*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);

参考内容：
1.http://blog.sina.com.cn/s/blog_5fdfa7e60100fbb1.html
2.http://www.openfoam.org,
3.http://openfoamwiki.net/index.php/SimpleFoam