[转载]OpenFOAM中的不可压缩湍流流动求解器turbFoam的说明

su_junwei

本文谈谈OpenFOAM不可压缩湍流流动求解器turbFoam的实现细节。
   (1) 求解器位置：applications/solvers/incompressible/turbFoam
   (2) 求解器文件夹结构
   |－Make
   |     |-options
   |     |-files
   |-createFields.H
   |-turbFoam.C
   (3) 求解器功能
       任意不可压缩湍流流动，湍流模拟采用雷诺时均方法
   (4) 文件说明
       1.options  //编译选项，用于指定编译用到的头文件位置及其动态库
//文件内容
#用到的头文件文件夹
EXE_INC = \
#雷诺是均湍流模型头文件  
    -I$(LIB_SRC)/turbulenceModels/RAS \
#传输模型头文件，牛顿流体或者非牛顿流体选择。  
    -I$(LIB_SRC)/transportModels \
#有限容积方法头文件
    -I$(LIB_SRC)/finiteVolume/lnInclude

#用到的动态连接库
EXE_LIBS = \
#不可压缩雷诺时均模型库
    -lincompressibleRASModels \
#不可压缩传输模型库（牛顿流体传输模型和非牛顿流体传输模型）
    -lincompressibleTransportModels \
#有限容积库
    -lfiniteVolume \
#网格相关工具库
    -lmeshTools

    2.files  //用于指定当前要编译的文件，这里不包含头文件，都是*.C文件。
//文件内容
turbFoam.C   //主文件
//编译后求解器的名字和存放位置
EXE = $(FOAM_APPBIN)/turbFoam

    3.createFields.H //创建场
    //提示，并创建压力场，各项意义，参看本站博文“OpenFOAM>>solver>>basic>>potentialFoam的说明”
    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
//提示，并创建速度场，各项意义，参看本站博文“OpenFOAM>>solver>>basic>>potentialFoam的说明”
    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
//创建表面流律场，面心存储，用于将体积分转化面积分，文件位于
//src &#187; finiteVolume &#187; cfdTools &#187; incompressible
#   include "createPhi.H"

//设定压力参考值和参考cell
    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

//建立传输模型，这里laminarTransport并不是说本solver为层流模型。
    singlePhaseTransportModel laminarTransport(U, phi);
//创建湍流模型，返回指向incompressible::RASModel的autoPtr指针。
//autoPtr指针有点复杂，后面将会有专门博文探讨该指针，目前你就将其当作普通指针就行了。
    autoPtr<incompressible::RASModel> turbulence
    (
        incompressible::RASModel::New(U, phi, laminarTransport)
    );

   4.turbFoam.C
//头文件
//有限容积库相关的所有头文件
#include "fvCFD.H"
//单相传输模型
#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
//雷诺时均模型
#include "incompressible/RASModel/RASModel.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//主程序入口
int main(int argc, char *argv[])
{
//设置case文件夹结构，程序运行时候输入的参数
#   include "setRootCase.H"
//创建时间对象runTime用于控制程序运行
#   include "createTime.H"
//创建网格对象mesh
#   include "createMesh.H"
//创建场，参看上面的说明
#   include "createFields.H"
//初始化连续误差。
#   include "initContinuityErrs.H"

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

    Info<< "\nStarting time loop\n" << endl;
//程序主循环
    for (runTime++; !runTime.end(); runTime++)
    {
//显示物理时间
        Info<< "Time = " << runTime.timeName() << nl << endl;
//读入piso控制参数
#       include "readPISOControls.H"
//计算courant参数
#       include "CourantNo.H"

        // Pressure-velocity PISO corrector
        {
            // Momentum predictor
//创建动量方程。关于下面书写紧凑性的背后问题，清参看本站博文“OpenFOAM中的神奇方程定义方式的背后”
            fvVectorMatrix UEqn
            (
                fvm::ddt(U)
              + fvm::div(phi, U)
         //该项包含两部分：1.分子扩散项 2.雷诺应力的偏分量的散度。
         //雷诺应力主分量的散度归结到了压力项中，这是大多数雷诺时均和大涡模型实现的一贯做法。
         //因此压力比层流模型中多了一项,那就是雷诺应力的主分量，通常被称为湍动压力。
              + turbulence->divDevReff(U)  
            );
//如果需要动量预测，则求解速度方程。
            if (momentumPredictor)
            {
                solve(UEqn == -fvc::grad(p));
            }

            // --- PISO loop
//PISO流程和icoFoam完全一样，请参看本站博文“OpenFOAM>>solver>>incompressible>>icoFoam的说明”
            for (int corr=0; corr<nCorr; corr++)
            {

                volScalarField rUA = 1.0/UEqn.A();

                U = rUA*UEqn.H();
                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);
                    pEqn.solve();

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

#               include "continuityErrs.H"

                U -= rUA*fvc::grad(p);
                U.correctBoundaryConditions();
            }
        }
       //重新计算湍动黏性，也就是对turbulence->divDevReff(U)需要的量进行更新。
       //比如k－e模型中下面函数包括求解k－e方程和计算有效黏性系数。
        turbulence->correct();
       //输出计算结果到文件中
        runTime.write();
       //显示当前执行时间和cpu时间。
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    //显示程序结束
    Info<< "End\n" << endl;
    //返回0
    return(0);
}
   （5）编译求解器
    打开控制台，进入求解器文件夹，输入wmake

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

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