OpenFOAM程序开发入门（2）

五、例子：在icoFoam中加入温度场求解

准备：

能量控制方程：
$\frac{\partial \mathbf{T}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{T})=\nabla \cdot(DT \nabla \mathbf{U})$
其中DT为热扩散率。在壁面上给定值条件。

需要解决的问题：

如何创建标量场T
如何创建热扩散率DT
如何定义温度方程，并求解
如何在算例中设定T和DT
如何设定T的离散格式
如何设定T的求解器的收敛标准等

步骤：

1、创建程序需要的新物性和新变量场
打开myicoFoam.C可以看到，程序开始运行时调用CreateFields.H，创建变量场。
打开CreateFields.H，可以看到程序首先从transportProperties文件中读入物性，

Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",              //从字典文件transportProperties读入
            runTime.constant(),                 //transportProperties文件位于目录runTime.constant()中
            mesh,                               //网格对象                  
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
);                                              //创建了Iodictionary类型对象 transportProperties

dimensionedScalar nu                            //首先读入粘性系数
    (
        transportProperties.lookup("nu")
    );                                          //创建有量纲标量nu，nu通过从字典transportProperties查找”nu”来赋值
//添加新方程需要的物性DT
dimensionedScalar DT                            //首先读入热扩散率
    (
        transportProperties.lookup("DT")
    );                                                         //创建有量纲标量DT，DT通过从字典transportProperties查找”DT”来赋值
//此外还要从createFields中读入p,U场，我们要加入的新的变量场为温度场T，最快的加入温度场的方法是拷贝p场的代码，修改为
Info<< "Reading field T\n" << endl;
    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

这样，创建了新的vol标量场T，从文件T中读入。
对于T的创建具体解释如下：

创建了标量场T
T通过读(IOobject::MUST_READ)在runTime.timeName()目录下名称为“T”的文件创建，在开始计算时，runTime.timeName()是contorlDict中设定的startTime值决定的。
T将自动写入(IOobject::AUTO_WRITE)计算结果到runTime.timeName()目录中,runTime.timeName()随迭代是变化的，写入控制由contorlDict中设定。
T是定义在mesh对象上的，这意味着T在内部cell上有值internalField，在边界上还需要边界条件，这与polyMesh/boundary中要一致。

2、在求解器中加入新的求解方程
下一步回到myicoFoam.C加入新的微分方程，由于温度场依赖于速度场，可放在PISO循环后面。

# include "continuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
//Add the temperature equation
fvScalarMatrix TEqn //温度是标量方程  
(
  fvm::ddt(T) 
+ fvm::div(phi, T)      //要用到界面流量
- fvm::laplacian(DT, T)  //扩散项
);
TEqn.solve();          //求解

3、编译

wmake

4、在算例中加入新方程的初始和边界条件

拷贝一个cavity算例到mycavity
修改transportProperties字典文件，设定DT

cd constant

修改transportProperties文件，前面已提到DT要从该字典文件读入。设定DT=0.002m2/s

DT              [0 2 -1 0 0 0 0] 0.002;

修改T文件，设定初始值和边界

cd 0 #进入0目录

拷贝一个T文件

cp p T

修改T文件为

FoamFile
{
    version        2.0;
    format         ascii;
    class          volScalarField;
    object         T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions      [0 0 0 1 0 0 0];
internalField  uniform 300;    //初始内部点为300℃
movingWall      
{
    type            fixedValue;
    value uniform  350;        //边界为350℃
}
fixedWalls      
{
    type            fixedValue;
    value uniform  300;        //边界为300℃
}

5、修改离散格式和代数求解器求解控制文件

进入system目录
由于温度方程有非稳态项，对流项，扩散项，分别要在ddt，div，laplacian中设置
打开fvSchemes文件，添加

divSchemes
{
    default        none;
    div(phi,U)      Gauss upwind;
    div(phi,T)      Gauss upwind;
}

laplacianSchemes
{
    default        none;
    laplacian(nu,U) Gauss linear corrected;
    laplacian(DT,T) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
}

在fvSolution中设置代数求解器选项

    T
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    };

注意T方程形成的矩阵是非对称的，不要用PCG和DIC
6、运行

blockMesh
myicoFoam . mycavity