icoFoam主程序代码分析

在分析这一部分代码之前，需要先了解该求解器的控制方程以及其基本离散格式，具体可参考东岳老师的网站：

在理解了icoFoam求解器所要解决的问题以及相应方法后，需要对主程序代码做进一步分析：本文参考了[openfoamwiki.net]和CloudBird07的文章做了进一步整理。

// 有限体积库头文件的集合
#include "fvCFD.H"  
//定义PISO循环,使用PISO循环必备头文件
#include "pisoControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//进入主程序
int main(int argc, char *argv[])
{
    #include "setRootCase.H"  //根据输入参数argc和argv设置算例根目录rootcase,必备头文件
    #include "createTime.H"  //创建时间对象,由runTime控制,非定常求解器必备头文件
    #include "createMesh.H"  //由 $case/constant/polyMesh,创建网格对象
//从网格mesh对象构造类pisoControl(基类:solutionControl)的对象piso
    pisoControl piso(mesh);   

    #include "createFields.H"  //创建场对象
    #include "initContinuityErrs.H"  //初始化连续误差

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;  //显示进入时间循环

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;//显示当前运行时间

        #include "CourantNo.H"          //计算库朗数并输出至屏幕
// Momentum predictor，去除压力项，构造动量方程算子，隐式离散
        fvVectorMatrix UEqn    //构造体积矢量场矩阵
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
//如果进行速度预测（on），则求解完整的动量方程得到预测速度
//如果不进行速度预测(off），则预测速度直接取当前已知时间步的速度
//momentumPredictor：动量预测求解开关，对多相流以及低雷诺数一般设置为off；
        if (piso.momentumPredictor())
        {
//求解动量方程，由上一时刻的压力p进行计算，因此又称速度预测
            solve(UEqn == -fvc::grad(p));//solve是FOAM名称空间的全局函数,参数为矩阵(fvMatrix)
        }

        // --- PISO loop
        while (piso.correct())
        {
// rAU:在速度的最后一个解中，从矩阵中提取对角项并存储倒数
//由于对流的非线性，矩阵系数是U的函数
            volScalarField rAU(1.0/UEqn.A());
//UEqn.H是UEqn减去非对角项和U的乘积
//Jocobi 迭代
//可参考constrainHbyA的源码，这部分主要是代码上的考虑
//返回结果:rAU*UEqn.H()
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
//将上述体积矢量场转化为面心标量场（参考有限体积法）
//fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)的作用是保证速度通量的全局守恒，以确保压力方程有解
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)  //(fvc::interpolate(U) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );
//调整边界值，保证速度边界条件守恒
            adjustPhi(phiHbyA, U, p);
//更正压力边界条件，保证通量守恒
            constrainPressure(p, U, phiHbyA, rAU);

// 网格非正交压力修正循环 
//如果在fvSolution字典文件中设置为0,就只求解控制方程一次,如果设置为n,则求解控制方程n-1次         
// Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector
// 构造压力方程进行求解
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );
//在不可压缩流动中，只有相对压力才起作用。除非存在边界压力
//一个网格的压力可以任意设定，以产生独特的压力解决方案
                pEqn.setReference(pRefCell, pRefValue);
//通过查询system/fvSolution下的PISO更新压力参考网格并重新设定参考值
                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();
        }
//输出计算结果，写入时间文件夹
        runTime.write();
//显示求解器执行时间及CPU耗时
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
//提示运行结束。
    Info<< "End\n" << endl;
//返回 0
    return 0;
}