#include "fvCFD.H"                            //有限体积法头文件
#include "singlePhaseTransportModel.H"        //声明singlePhaseTransportModel类，基于粘性的单相流动输运模型
#include "turbulentTransportModel.H"        //湍流模型的命名空间
#include "noRadiation.H"                    //声明noRadiation类，“不考虑辐射的辐射模型”
#include "fvOptions.H"                        //声明option类，提供一系列成员函数用于执行对源项的操作，例如T方程中即将用到的constrain()和correct()
#include "simpleControl.H"                    //声明simpleControl类，提供一系列成员函数用于执行simple循环指令,例如下面的loop()

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

int main(int argc, char *argv[])
{
    #include "postProcess.H"        //声明后处理函数，以便实时调用

    #include "setRootCaseLists.H"    //检查case文件结构
    #include "createTime.H"            //根据controlDict中的设置，确认当前时间步，并定义runTime对象
    #include "createMesh.H"            //定义前述fvMesh类的对象mesh
    #include "createControl.H"        //定义前述simpleControl类的对象simple
    #include "createFields.H"        //初始化各物理量场，见下一节
    #include "initContinuityErrs.H"    //定义并初始化连续性误差

    turbulence->validate();        //在开始求解前，确认已正确计算湍流粘度，此为部分湍流模型的功能

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

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

    while (simple.loop(runTime))    //执行simple循环
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        {
            #include "UEqn.H"    //解动量离散方程
            #include "TEqn.H"    //解温度离散方程
            #include "pEqn.H"    //解压力修正值方程和速度修正值方程
        }

        laminarTransport.correct(); //修正运动粘度nu
        turbulence->correct();    //解湍流变量离散方程

        runTime.write();        //写数据到磁盘

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

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

    return 0;
}

• creatFields.H

Info<< "Reading thermophysical properties\n" << endl;

Info<< "Reading field T\n" << endl;
volScalarField T                //通过volScalarField类的构造函数，定义一个对象T
(                                //该构造函数有两个参数
    IOobject                    //一个是对IOobject构造函数的调用
    (                              //该构造函数有五个参数
        "T",                      //变量名称
        runTime.timeName(),          //由runTime对象的成员函数timeName()获取当前时间步，与时间关联
        mesh,                      //网格对象mesh，与空间关联
        IOobject::MUST_READ,      //读写属性1：必读，因而需要用户在case中提供T文件
        IOobject::AUTO_WRITE      //读写属性2：自动写
    ),
    mesh                        //另一个参数是对象mesh，通过mesh对象间接读取数据进行初始化
);

Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh            //类似的，定义标量场p_rgh，必读，自动写，需要用户提供p文件
(
    IOobject
    (
        "p_rgh",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U                //类似的，定义向量场U，必读，自动写，需要用户提供U文件
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

#include "createPhi.H"            //调用头文件createPhi.H
/* 头文件中的语句为
surfaceScalarField phi            //定义面向量场phi, 即面通量，自动写，无需要用户提供phi文件
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    fvc::flux(U)
); */

#include "readTransportProperties.H"    //调用头文件 readTransportProperties.H
        //其中使用singlePhaseTransportModel类的构造函数定义了laminarTransport对象，

Info<< "Creating turbulence model\n" << endl;
autoPtr<incompressible::turbulenceModel> turbulence  //使用智能指针实例化一个对象turbulence
( //用incompressible::turbulenceModel类的构造函数，进行初始化
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

// Kinematic density for buoyancy force
volScalarField rhok                    //定义标量场rhok，默认读写属性，无需提供相应文件
(
    IOobject
    (
        "rhok",
        runTime.timeName(),
        mesh
    ),
    1.0 - beta*(T - TRef)            //rhok值由表达式（10）进行更新     
);

// kinematic turbulent thermal thermal conductivity m2/s
Info<< "Reading field alphat\n" << endl; 
volScalarField alphat                 //定义标量场alphat，必读，自动写，需用户提供alphat文件
(
    IOobject
    (
        "alphat",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


#include "readGravitationalAcceleration.H" //从case文件中读取指定的重力加速度
#include "readhRef.H"                       //调用头文件readhRef.H
/*执行的语句为
Info<< "\nReading hRef" << endl;
    uniformDimensionedScalarField hRef     //读取一个均一的带单位的标量场hRef，同式(6)
    (
        IOobject
        (
            "hRef",
            runTime.constant(),            //不随时间变化
            mesh,
            IOobject::READ_IF_PRESENT,     //可指定，也可不给
            IOobject::NO_WRITE               //不会自动生成文件
        ),
        dimensionedScalar(dimLength, 0)    //传一个值为0，单位为长度单位的标量作初始化
    );
*/

#include "gh.H"                                //调用头文件gh.H
/*执行的语句为
dimensionedScalar ghRef(- mag(g)*hRef);     // 定义标量场 ghRef = -|g|×hRef 
volScalarField gh("gh", (g & mesh.C()) - ghRef); //gh 见式(7)，体心值，
surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef); //gh 见式（7），面心值
            // & 为重载的运算符，代表点乘，mesh.C()为体心位置矢，mesh.Cf()为面心位置矢
*/

volScalarField p                        //标量场p,不需要提供p文件
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    p_rgh + rhok*gh                        //用p_rgh, rhok和gh进行初始化
);

label pRefCell = 0;                        //默认压力参照点
scalar pRefValue = 0.0;
setRefCell
(
    p,
    p_rgh,
    simple.dict(),
    pRefCell,
    pRefValue
);

if (p_rgh.needReference())           
{
    p += dimensionedScalar           //整体加减，以使得默认参照点处的压力等于默认参考压力
    (
        "p",
        p.dimensions(),
        pRefValue - getRefCellValue(p, pRefCell)
    );
}

mesh.setFluxRequired(p_rgh.name());   

#include "createIncompressibleRadiationModel.H" //辐射模型
#include "createFvOptions.H"        //定义option类的对象fvOptions

• readTransportProperties.H

singlePhaseTransportModel laminarTransport(U, phi); //定义了前述singlePhaseTransportModel类的对象laminarTransport

// Thermal expansion coefficient [1/K]
dimensionedScalar beta            //定义一个有单位的标量beta
(
    "beta",
    dimless/dimTemperature,        //单位为 无量纲/温度单位  即1/K
    laminarTransport            //从上述laminarTransport对象中读取，laminarTransport会从transportProperties字典中读取
);

// Reference temperature [K]
dimensionedScalar TRef("TRef", dimTemperature, laminarTransport);

// Laminar Prandtl number
dimensionedScalar Pr("Pr", dimless, laminarTransport);

// Turbulent Prandtl number
dimensionedScalar Prt("Prt", dimless, laminarTransport);

• UEqn.H

    tmp<fvVectorMatrix> UEqn    
    (
        fvm::div(phi, U)          //对速度求散度(隐式加入等号左侧系数矩阵)
      + turbulence->divDevReff(U) //见式（3）第二项
                                    //对偏应力(deviatoric stress)求散度（divergence）
                                    //其中的R代表雷诺时均，eff代表总的扩散系数,
                                    //即运动粘度nu及湍流扩散系数nu_t的和nu_eff，
     ==
        fvOptions(U)              //源项
    );

    UEqn.relax();                //速度亚松弛

    fvOptions.constrain(UEqn);  //通过源项实现特定区域速度恒定(若指定)

    if (simple.momentumPredictor()) //若fvSolution文件中SIMPLE下的momentumPredictor设为YES
    {
        solve
        (
            UEqn
          ==
            fvc::reconstruct        //由以下的面心值得到体心值，显式加入等号右侧
            (
                (
                  - ghf*fvc::snGrad(rhok) //ghf代表gh的面心值，再与 “折算密度” rhok的梯度相乘
                  - fvc::snGrad(p_rgh)    //p_rgh的梯度
                )*mesh.magSf()            //乘以面积   同式（9）
            )
        );

        fvOptions.correct(U);   //通过源项实现速度限值修正(若指定)
    }

• TEqn.H

{
    alphat = turbulence->nut()/Prt;        //求湍流热扩散系数 alpha_t = nut / Prt
    alphat.correctBoundaryConditions(); //更新边界

    volScalarField alphaEff("alphaEff", turbulence->nu()/Pr + alphat); //求总热扩散系数

    fvScalarMatrix TEqn                    //温度方程同(4)
    (
        fvm::div(phi, T)
      - fvm::laplacian(alphaEff, T)
     ==
        radiation->ST(rhoCpRef, T)        //实际上没有考虑辐射，因为头文件中默认选取的辐射模型为无
      + fvOptions(T)                    //源项
    );

    TEqn.relax();                        //温度亚松弛

    fvOptions.constrain(TEqn);            //通过源项实现指定区域温度恒定(若设置)

    TEqn.solve();                        //解温度方程

    radiation->correct();                //实际上什么都不做

    fvOptions.correct(T);               //通过源项实现温度限值修正(若指定)


    rhok = 1.0 - beta*(T - TRef);        //根据新的温度场更新rhok
}

• pEqn.H

{
    volScalarField rAU("rAU", 1.0/UEqn.A());      // rAU = 1/A
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));   //rAUf = (1/A)_f
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));  // HbyA = H / A 见(12)

    UEqn.clear();

    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
                                 //phig = b/A = - gh * grad(rhok) / A  见（12）

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::flux(HbyA)
    );                          //phiHbyA = (HbyA)_f  见（16）

    adjustPhi(phiHbyA, U, p_rgh);

    phiHbyA += phig;            //  把 b/A 并入 H/A

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p_rgh, U, phiHbyA, rAUf);  

    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) //参考压力修正值方程 同（16）
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

        p_rghEqn.solve();       //求解 参考压力修正值方程

        if (simple.finalNonOrthogonalIter())
        {
            // 计算面通量
            phi = phiHbyA - p_rghEqn.flux();
            p_rgh.relax();       //压力方程亚松弛

            // 速度修正值方程，可见下式中速度由参考压力直接得到，没有求解线性方程组，所谓显式修正
            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }

    #include "continuityErrs.H"

    p = p_rgh + rhok*gh;            //更新p

    if (p_rgh.needReference())      
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rhok*gh;
    }
}

• 边界条件

OpenFOAM安装文件中自带了两个算例，位于$FOAM_TUTORIALS/heatTransfer/buoyantBoussinesqSimpleFoam/路径。

遗憾的是这两个算例都为密闭空间自然对流，而对于有开口存在的问题，提供以下参考。

一般对于单纯的速度-压力耦合问题，OpenFOAM对边界条件组合的评价如下：

而在具体应用时，这些边界条件没有绝对的优劣之分，越是复杂的问题越是视阴晴而定。举例来说，在一些回流存在但并不占优的问题中巧妙地使用充分发展假设可以让问题简化而达到收敛，而有时充分照顾到回流求解反而不收敛，所谓给脸不要敛😑。

对于有浮升力存在的换热问题，问题显得稍有不同，因为采用了压力变换的思路，引出了p_rgh，与速度的耦合关系有相应的变动。研读OpenFOAM自带的算例可以看到，对p_rgh 采用的是fixedFluxPressure边界，所求得的压力梯度满足动量方程，有助于计算的稳定。很多的反馈表明在使用其他出口边界条件都会发散的情况下，使用fixedFluxPressure边界能使问题更好地收敛。这一边界条件同样适用于壁面。

对于温度场，需要指定T和α_t的边界条件，在流域内部，α_t由湍流模型算得，在壁面则需要指定壁面函数。

⚠值得注意的是，从热力学的角度出发，对于一个合理的稳态解，能量进出应是匹配的，这在设置边界条件时应该给予考虑。如果是有内热源的问题，对于闭口系可设置与热源强度相当的热流边界，对于开口系也要考虑热力平衡。

• 差分格式

在求解过程中，会有各种问题造成求解发散。除了结合物理实际使用合适的边界条件外，还应使用合适的差分格式和松弛因子以保证稳定性。

👉🏼对于差分格式，CFD理论告诉我们高阶格式会带来over-shoot，这在温度场的计算中有直观的体现，温度常常溢出应有的范围；有时如果初始场设置不合理，使用高阶格式也会让求解在一开始就发散。在不妥协阶数的前提下，可以把二阶迎风换为带限的差分格式如limitedLinear和vanLeer看是否有好转，进一步的妥协是采用不会造成over-shoot的一阶迎风(upwind)，这一格式具有假扩散的特性，精度低，但能够帮助排查问题究竟是不是出在差分格式上，如果还是不稳定甚至发散，则检查边界条件、松弛因子以及思考对应的物理问题是否能发展到稳态。

👉🏼对于难伺候的问题，可以先算得等温流动的初始流场，而后设置所需的温度边界条件，用一阶迎风让温度场参与计算，残差稳定或收敛后再对温度使用limitedLinear或vanLeer格式，并使用limitTemperature源项对温度范围进行限制。

• 松弛因子

对于松弛因子，SIMPLE算法的修正思路明确对速度和压力的亚松弛处理提出要求；同时，由于流场与温度场的耦合，在循环体内，ρ_k 是先使用速度场的预测值计算温度场得到的，而后再进入压力方程参与计算，所以在温度场求解中也应做亚松弛处理。对于湍流模型也是类似。在实践中，可以先从这样的松弛因子开始：

如果发现某个量不稳定或发散，则可适当调小它的松弛因子。松弛因子只是保证求解器的稳定性，但并不能让由不合理的边界条件或物理描述带来的发散得到解决，更多时候还是要多从其他设置上找原因。

本文介绍了OpenFOAM中用于求解有浮升力存在时不可压流换热问题的求解器buoyantBoussinesqFoam，列出了理论推导的关键步骤，介绍了SIMPLE算法的迭代思路，对程序的大部分语句给出了注释及理论对照，最后分享了一些求解器使用上的经验和技巧，希望对大家有所帮助。由于本人水平和时间有限，难免有错漏荒唐之处，还请不吝指正。