【转载】pimple算法简述

su_junwei

终于回国了，事情比较多，blog更新速度慢了些。  今天一起看看OpenFOAM中的pimple算法流程。 pimple算法是simple算法和piso算法的结合体。

pimple的基本思想是：将每个时间步长内用simple稳态算法求解（也就是将每个时间步内看成稳态流动），时间步长的步进用piso算法来完成。

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



#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())  //时间步进
    {
        #include "readTimeControls.H"
        #include "readPIMPLEControls.H"  //pimple控制
        #include "CourantNo.H"   //courant数
        #include "setDeltaT.H"  //设置时间步长

        runTime++;  //步进时间

        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity PIMPLE corrector loop
        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)   //外循环 simple修正
        {
            if (nOuterCorr != 1)
            {
                p.storePrevIter();  //每次流程开始存储p，以便连续性方程亚松弛
            }

            #include "UEqn.H"  //速度方程

            // --- PISO loop
            for (int corr=0; corr<nCorr; corr++)  //piso修正
            {
                #include "pEqn.H"
            }

            turbulence->correct();  //湍流修正
        }

        runTime.write();  //输出

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}



tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)
);

if (oCorr == nOuterCorr-1)
{
    UEqn().relax(1);  //最后一步不进行亚松驰，以便采用标准的piso流程
}
else
{
    UEqn().relax();
}

volScalarField rUA = 1.0/UEqn().A();

if (momentumPredictor)
{

    //下面的条件语句可以使得最后一次的速度的求解残差和以前修正不一样。
    if (oCorr == nOuterCorr-1)
    {
        solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));   
    }
    else
    {
        solve(UEqn() == -fvc::grad(p));
    }
}
else
{
    U = rUA*(UEqn().H() - fvc::grad(p));
    U.correctBoundaryConditions();
}
压力方程

U = rUA*UEqn().H();

if (nCorr <= 1)
{
    UEqn.clear();
}

phi = (fvc::interpolate(U) & mesh.Sf())
    + fvc::ddtPhiCorr(rUA, U, phi);

adjustPhi(phi, U, p);

// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
    // Pressure corrector
    fvScalarMatrix pEqn
    (
        fvm::laplacian(rUA, p) == fvc::div(phi)
    );

    pEqn.setReference(pRefCell, pRefValue);

    if
    (
        oCorr == nOuterCorr-1
     && corr == nCorr-1
     && nonOrth == nNonOrthCorr
    )
    {
        pEqn.solve(mesh.solver("pFinal"));
    }
    else
    {
        pEqn.solve();
    }

    if (nonOrth == nNonOrthCorr)
    {
        phi -= pEqn.flux();
    }
}

#include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector except for last corrector
if (oCorr != nOuterCorr-1)   //最后一次不采用亚松驰，采用标准piso流程。
{
    p.relax();
}

U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();


上面没有注释的请参看本站其他求解器的说明。

转自OpenFOAM研究：http://blog.sina.com.cn/openfoamresearch

OpenFOAM

兄弟! 欢迎回来!!