[转载]OpenFOAM>>solver>>incompressible>>icoFoam的说明

su_junwei

//createFields.H  
//读入属性提示
Info<< "Reading transportProperties\n" << endl;
//声明属性字典类对象，该对象由constant 文件夹下的“transportProperties”初始化创建。
    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",//字典文件名字
            runTime.constant(), //位置
            mesh,               //网格对象
            IOobject::MUST_READ, //必须读，由文件创建
            IOobject::NO_WRITE   //不对该文件进行写
        )
    );
//字典查询黏性一遍初始化带有单位的标量
    dimensionedScalar nu
    (
        transportProperties.lookup("nu")
    );
//提示创建压力场
    Info<< "Reading field p\n" << endl;
//创建压力场
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

//提示并创建速度场
    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

//创建界面流律
#   include "createPhi.H"

//设定压力参考cell的index和参考值
    label pRefCell = 0;
    scalar pRefValue = 0.0;
//通过查询system/fvSolution下的PISO进行更新压力参考cell更新并设定参考值。
    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

//icoFoam.C
int main(int argc, char *argv[])
{

//下面是很多solver都包括的，前面solver介绍时候已经说明。
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"
#   include "createFields.H"
#   include "initContinuityErrs.H"

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

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

// 用runTime 对象控制时间推移
    for (runTime++; !runTime.end(); runTime++)
    {
//显示当前物理时间
        Info<< "Time = " << runTime.timeName() << nl << endl;
//读入piso控制参数，写在runTime循环中，以便运行过程中的修改生效
#   include "readPISOControls.H"
//计算courant数
#   include "CourantNo.H"
      //构造速度方程
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

//应当指出上面方程中并没有包含grad(p),动量方程中所有的项要么在压力方程中要么在动量方程中，而不同时在两个方程中。
//求解动量方程，下面方程中的p由上一时刻p计算，下面求解称为速度预测。

        solve(UEqn == -fvc::grad(p));
        // --- PISO loop---- 速度修正步
        for (int corr=0; corr<nCorr; corr++)
        {
        //根据预测的速度值求解系数矩阵对角元素的倒数。
            volScalarField rUA = 1.0/UEqn.A();

// 更新速度， See Hrv Jasak's thesis eqn. 3.137 and Henrik Rusche's thesis, eqn. 2.43

            U = rUA*UEqn.H();

// 计算表面流律
// The ddtPhiCorr用来考虑差值速度场计算的流律和实际流律的差别
            phi = (fvc::interpolate(U) & mesh.Sf())
                + fvc::ddtPhiCorr(rUA, U, phi);
//调整边值，以便满足连续性条件
            adjustPhi(phi, U, p);
//非正交循环
            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
  //设定压力方程并进行求解
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rUA, p) == fvc::div(phi)
                );
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
    //根据新求解的压力，更新表面流律
                if (nonOrth == nNonOrthCorr)
                {
                    phi -= pEqn.flux();
                }
            }  // end of non-orthogonality looping
//计算连续性方程误差
#           include "continuityErrs.H"

//根据新的压力，修正速度

            U -= rUA*fvc::grad(p);
//修正速度边界(主要针对第二类边界条件）
            U.correctBoundaryConditions();

        } // end of the PISO loop
//输出计算结果
        runTime.write();
//显示执行时间，CPU时间
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;

    } // end of the time step loop
//提示运行结束。
    Info<< "End\n" << endl;
//返回0
    return(0);
}
转自OpenFOAM研究：http://blog.sina.com.cn/openfoamresearch

[ 本帖最后由 su_junwei 于 2009-5-2 18:28 编辑 ]