OpenFOAM中预编译了大量适用于各种应用场景的求解器，但有时可能希望向现有求解器中添加一些内容，或者基于新的研究开发新的求解器。本案例将演示如何将温度的标量输运方程添加到icoFoam求解器。

通过利用OpenFOAM的现有功能及其独特的体系结构，实现此求解器只需要对现有代码进行少量的修改即可。

本案例中，温度被模拟为一个守恒的被动标量，也就是说其不会影响流体的压力或速度，这种情况仅适用于温度变化相对较小的问题中。从速度场到温度的单向耦合通过温度方程中的对流项来考虑：

1 程序结构

求解器icoThermalFoam从icoFoam改编而来，因此先从solver目录中将icoFoam的源代码拷贝到工作目录。

利用下面的命令完成文件准备工作。

cd $FOAM_RUN
mkdir Demo && cd Demo
cp -r $FOAM_SOLVERS/incompressible/icoFoam icoThermalFoam
cd icoThermalFoam

文件结构如下图所示。

• 修改icoFoam.C的文件名为icoThermalFoam.C
• 修改files文件的内容

icoThermalFoam.C
EXE = $(FOAM_APPBIN)/icoThermalFoam

这里的$FOAM_APPBIN为环境变量，求解器编译到此位置后在全局均可直接调用。

2 修改头文件

求解器包含了一个createFields.H的头文件，该文件从OpenFOAM案例中读取各种字典文件来获取求解数据。本案例中需要从字典文件中获取参数DT以及T，其中DT为物性参数，T为场变量。

修改createFields.H头文件，在其中添加下面的代码：

// 从字典文件中读取DT，这里读的是物性参数，从transportProperties字典文件中读取
dimensionedScalar DT
(
    "DT",  // 字典文件中关键字的名称，后面在案例字典中指定的关键字必须与一致
    dimViscosity,   // 指定量纲，这里利用预置的量纲dimViscosity，也可以使用dimensionSet自己定义量纲
    transportProperties   // 前面定义的transportProperties字典
);
 

// 读取各边界的温度数据，参照p文件
Info << "Reading field Tn" << endl;
// 定义标量场T
volScalarField T  
(
    // 创建IOobject对象作为标量场的参数，IOobject对象指定了标量的参数
    IOobject
    (
        "T",      // 定义了文件名为T，后面准备案例文件时在0文件夹中需要有相应的T文件
        runTime.timeName(),  // 文件位置
        mesh,				// 注册对象
        IOobject::MUST_READ,  // 指定数据必须读取，若0文件夹中没有T文件，则报错
        IOobject::AUTO_WRITE  // 指定数据自动写入，在源文件中的write()函数会将计算结果数据写入T文件中
    ),
    mesh
);

这里是常规的文件读取，前面已经讲过，就不细述了。

3 修改源文件

需要在源文件icoThermalFoam.C中写入控制方程。

由于温度T的输运方程与压力和速度是解耦的，因此可以将温度求解放在压力速度修正之后。

修改后的源文件为（文件前半部分没有列出）：

            #include "continuityErrs.H"
            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }
        //将温度控制方程加到这里
        // --------------求解温度控制方程-------------------
        solve
        (
            fvm::ddt(T) + fvm::div(phi,T) == fvm::laplacian(DT,T)
        );
        //------------------------------------------------
        runTime.write();  // 将所有指定为写入的物理场数据写入到文件中
 
        runTime.printExecutionTime(Info);
    }
    Info<< "Endn" << endl;
    return 0;
}

修改完毕后可以通过wmake编译程序。

编译完毕后可以输入icoThermalFoam -help测试是否求解器是否能够调用。如下图所示。

4 准备测试文件

从前面的源代码可以看出，求解器icoThermalFoam需要读入物性参数DT及场变量T。这里采用修改OpenFOAM自带案例cavityClipped来进行测试。

利用下面的命令准备文件：

cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity .
cd cavity
cd 0 && cp p T && cd ..
blockMesh

案例文件结构如下图所示。

• 修改constant/transportProperties文件

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      transportProperties;
}
// * * * * * * * * * * * * * //
nu              0.01;
DT              0.001;  // 关键字DT与前面程序中的名称保持一致，这里不需要指定量纲（V9版本可能不同）

由于在createFileds.H头文件中读取DT时指定了量纲为dimViscosity，因此在字典文件中不再需要指定量纲。

• 修改0/T文件夹，指定温度边界值与初始值

FoamFile
{
        version         2.0;
        format          ascii;
        class           volScalarField;
        location        "0";
        object          T;
}
// 指定量纲与初始值
dimensions              [0 0 0 1 0 0 0];
internalField           uniform 300;
// 指定边界条件
boundaryField
{
        lid
        {
                type    fixedValue;
                value   uniform 400;
        }
 
        fixedWalls
        {
                type    fixedValue;
                value   uniform 300;
        }
 
        frontAndBack
        {
                type    empty;
        }
}
 

• 修改system/fvSchemes指定离散格式

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * //
 
ddtSchemes
{
    default         Euler;
}
 
gradSchemes
{
    default         Gauss linear;
}
 
divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    // ---添加温度方程离散格式---
    div(phi,T)      Gauss linear;
}
 
laplacianSchemes
{
    default         Gauss linear orthogonal;
}
 
interpolationSchemes
{
    default         linear;
}
 
snGradSchemes
{
    default         orthogonal;
}

• 修改system/fvSolution指定T方程的求解方式

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * //
 
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }
 
    pFinal
    {
        $p;
        relTol          0;
    }
 
    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
    // ---添加T方程求解方式-----
    T
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
    //---------------------------
}
 
PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

5 测试求解

测试文件准备完毕后，可以在案例根目录下运行命令进行求解计算：

icoThermalFoam

计算完毕后查看结果，如下图所示。

6 小结

总结一下，利用OpenFOAM开发一个新的求解器的基本步骤包括：

• 从预置求解器中找一个最接近的求解器
• 添加物性参数读写及物理场变量读写的代码
• 添加控制方程

使用的时候主要注意：

• 物性参数指定
• 边界条件指定
• 添加的物理场的离散算法指定
• 方程组求解方式指定

其他的就没什么了。本案例极为简单，仅为演示而已，不过复杂的求解器，其实套路也都是一样的。

OpenFOAM高层编程其实还是比较简单的，主要是熟悉各种模板类的调用，搞清楚数据结构以及数据之间的调用关系基本上就能上手了。高层编程其实还是属于应用范畴，谈不上开发，这部分工作并不难，真正难的是开发底层算法而不是编写程序，OpenFOAM的架构设计非常漂亮，将各种算法转化为代码其实都是体力活儿而已。但是算法的推导就比较难了，不过万幸的是有那帮子搞纯理论的人在为这事儿掉头发。

”

（完）