初识OpenFOAM
写在最前：
本文的目的是通过几个基本例子来对OpenFOAM求解器有一个轮廓式的认识，文章所涉及到例子源于苏军伟新浪博客，若需要了解更多资料可以拜读该博客，比较惭愧的是因能力问题本文对原文没有多大改进，最多的仅是代码高亮的改进。
OpenFOAM是一个免费、开源的CFD软件包，由OpenCFD有限责任公司出品。它有着庞大的商业和科研用户基础，涉及工程、科学等领域。OpenFOAM求解的问题范围非常广，既能求解化学反应、湍流、热传递等复杂流动，又能求解固体动力学和电磁学等问题。OpenFOAM是一个完全由C++编写的面向对象的CFD类库，采用类似于我们日常习惯的方法在软件中描述偏微分方程的有限体积离散化，支持多面体网格（比如CD-adapco公司推出的CCM+生成的多面体网格），因而可以处理复杂的几何外形，支持大型并行计算，等。
每个软件都有其自身的独特特点，OpenFOAM也不例外，相比其他CFD软件，OpenFOAM具有很好的灵活性。举一个不切当的比喻就是：你可以把它认为是专用于计算流体问题的Matlab,用户可以通过调用各种函数来实现计算目的，这个过程中可以不用特别的专注于一些很基础的问题。
在开始之前，首先需要记在心里面：OpenFOAM是基于C++的，明白这一点对于下面的进一步阅读一定有帮助。本文一共涉及到三大例子，通过总结这些例子，大家一定会发现OpenFOAM是一个优秀的软件。
例：OpenFOAM>>solver>>basic>>laplacianFoam
//createFields.H
//－屏幕提示。Info 等价于 C++中std::cout,Info间接调用cout
Info<< "Reading field T\n" << endl;
//－声明一个标量场,网格中心存储变量。该场是通过读入文件(IOobject and MUST_READ)进行设置,并根据controlDict中的设置进行自动write,由write.H中的runTime.write()来执行write();。
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//- 提示读入参数控制文件
Info<< "Reading transportProperties\n" << endl;
//- 参数控制文件声明过文件形式读入
IOdictionary transportProperties
(
IOobject
(
"transportProperties", //文件名字
runTime.constant(), //文件位置，case文件夹中constant子文件夹
mesh,
IOobject::MUST_READ,//通过read一个文件，初始化
IOobject::NO_WRITE //并不根据时间对文件进行写
)
);
//－提出读入扩散律
Info<< "Reading diffusivity DT\n" << endl;
//－通过查询参数控制文件，初始化带有单位的标量，lookup中的“DT”为关键字
dimensionedScalar DT
(
transportProperties.lookup("DT")
);
//laplacianFoam.C
#include "fvCFD.H" //-cfd头文件，包括大多数cfd计算需要的头文件，在src » finiteVolume » cfdTools » general » include
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//主程序入口
int main(int argc, char *argv[])
{
//设置rootcase，根据输入参数argc 和 argv对
# include "setRootCase.H"
//－创建时间，下面的runTime控制
# include "createTime.H"
//创建网格，根据constant文件中polyMesh文件夹中的网格数据，创建网格对象mesh，位于src » OpenFOAM » include
# include "createMesh.H"
//创建场对象，在前面已经说明
# include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//提示计算温度分布
Info<< "\nCalculating temperature distribution\n" << endl;
//计算主控制流程
for (runTime++; !runTime.end(); runTime++)
{
//提示当前计算时间
Info<< "Time = " << runTime.timeName() << nl << endl;
//读入simple算法参数，位于
src » finiteVolume » cfdTools » general » include
# include "readSIMPLEControls.H"
//对于网格非正交循环修正。
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
//求解拉普拉斯方程,这里的solve是Foam空间的全局函数，内参数为fvMatrix，fvm表示隐式离散，返回有限容积稀疏矩阵类fvMatrix对象，具体对象中内容，以后说明
solve
(
fvm::ddt(T) - fvm::laplacian(DT, T)
);
}
//对求解变量进行写
# include "write.H"
//提示执行时间及CPU耗时
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
//提示程序结束
Info<< "End\n" << endl;
return(0);
}
//write.H
if (runTime.outputTime())
{
volVectorField gradT = fvc::grad(T); //计算温度梯度，向量场
//声明3个变量，分别以gradT的三个方向的分量进行初始化。
volScalarField gradTx
(
IOobject
(
"gradTx", //变量名字
runTime.timeName(), //位置
mesh, //mesh，主要用于对象注册，根据runTime进行写
IOobject::NO_READ, //
IOobject::AUTO_WRITE //自动写
),
gradT.component(vector::X) //用场初始话，vector::X 枚举变量，可直接写0 //gradT.component(0)
);
volScalarField gradTy
(
IOobject
(
"gradTy",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::Y)
);
volScalarField gradTz
(
IOobject
(
"gradTz",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::Z)
);
//对场进行写
runTime.write();
}
Tips：runTime是类Time的一个对象，runTime.timeName()其实就是当前运行到的物理时间（非稳态物理问题的时间），你程序运行到了5.1 s，runTime.timeName()就会将字符串“5.1”返回给你，OpenFOAM在对不同时刻的数据进行存取的时候就是靠着这个字符串。runTime.constant()返回的就是case下那个constant文件夹的名字，这个名字可以改的，默认为constant。OpenFOAM根据case文件夹里面的system/controlDict里面对输出的设置来确定当前时刻时候输出，如果当前时刻输出的话，outputTime()就为true，就输出数据了。
关于OpenFAOM中类的介绍，可以到苏军伟博客。
例：OpenFOAM>>solver>>basic>>potentialFoam
//creatField.H
//提示创建压力场
Info<< "Reading field p\n" << endl;
//压力场为标量场，网格中心存储变量（OpenFOAM用的是非结构化同位网格），下面为创建标量场压力，两个参数IOobject对象和网格对象mesh，IOobject主要从事输入输出控制
volScalarField p
(
IOobject
(
"p", //压力场初始文件的名字。
runTime.timeName(), //位置，该位置由case中的system/controlDict中的startTime控制
mesh, //网格对象，主要从事对象注册，以便由runTime.write()控制输出
IOobject::MUST_READ,//该对象由文件读出创建，因此，需要READ
IOobject::NO_WRITE //不输出压力场
),
mesh //压力场所用的网格对象，在createMesh.H创建
);
//压力场初始化为0，单位为上面文件中的单位。dimensionedScalar 为带单位的标量，初始化三个参数，名字，单位，数值。也可采用如下形式
//dimensionedScalar("zero",dimensionSet(1,-1,-2,0,0 ,0, 0),0.0);
//应当注意，OpenFOAM中的大部分case对动量的求解都是求解的速度场，压力单位初始化时候，一般为除去密度后的值.
//dimensionSet中有7个参数，他们依次为质量、长度、时间、温度、摩尔数、安培、坎德拉。
p = dimensionedScalar("zero", p.dimensions(), 0.0);
//提示读入速度场
Info<< "Reading field U\n" << endl;
//创建速度场，向量场，体心存储变量。
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE //自动写，根据runTime.write()或者U.write();
),
mesh
);
//初始化速度场。这里初始化速度为0 ，如果初始化Ux=1，Uy=2，Uz=3 （单位为国际单位）可采用
//U = dimensionedVector("0", U.dimensions(), vector(1,2,3));
U = dimensionedVector("0", U.dimensions(), vector::zero);
//表面场，phi，界面流率，存储在体之间的交接面上。表面场(surface...)不能和体积场(vol...)
//直接计算，因为他们存储在不同地方，大小不同。
//可以将体积场转化为表面场（运用fvc::interpolate())
//或者由表面场转化为体积场（运用fvc::reconstruct())进行计算。
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
//U是体积场，运用插值转为表面场和表面积场进行相乘来初始化流率
//mesh.Sf()返回网格交接面面积矢量。
fvc::interpolate(U) & mesh.Sf()
);
//压力参考cell
label pRefCell = 0;
//压力参考值
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
//只有求解区域所有的压力边界都为第二类边界条件是，上面的值才会用到。如果有第一类边界条件，
//压力参考值为这点边界值。对于不可压缩流动压力值为相对值，上面的参考值的大小对结果无影响。
//potentialFoam.C
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//主程序入口
int main(int argc, char *argv[])
{
//增加新的有效输入参数。
argList::validOptions.insert("writep", "");
//设置case目录
# include "setRootCase.H"
//创建时间，对计算流程进行控制，主要是名为runTime的对象
# include "createTime.H"
//创建网格对象mesh
# include "createMesh.H"
//创建场变量，前面已经说过
# include "createFields.H"
//读simple控制参数
# include "readSIMPLEControls.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//计算势流提示
Info<< nl << "Calculating potential flow" << endl;
//对流率进行调整，以保证方程的连续性
adjustPhi(phi, U, p);
//网格非正交性循环，如果网格是正交的，可以设定nNonOrthCorr=1
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
//创建压力方程，该方程为标量稀疏矩阵类
fvScalarMatrix pEqn
(
fvm::laplacian //拉普拉斯离散，隐式
(
dimensionedScalar
(
"1",
dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
1
),
p
)
==
fvc::div(phi)//速度离散，显示。因为是压力方程，其他变量只能显示
);
//设定压力参考值
pEqn.setReference(pRefCell, pRefValue);
//求解压力方程，调用的fvMatrix成员函数
pEqn.solve();
//流率修正，应当注意OpenFOAM中对压力本身求解，而非压力变化值。
//关于simple算法和PISO算法的OpenFOAM实现，会在simpleFoam及其icoFoam中详细说明。
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
//提示连续性方程求解误差
Info<< "continuity error = "
<< mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
<< endl;
//根据表面场重建速度场
U = fvc::reconstruct(phi);
//对速度场边界进行修正，主要针对第二类边界条件下的边界场
U.correctBoundaryConditions();
Info<< "Interpolated U error = "
<< (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
/sum(mesh.magSf())).value()
<< endl;
// Force the write
//直接对速度进行输出
U.write();
//界面流率输出。注意当前场存储在界面上，phi的大小（个数）和U的大小（个数）不相同的
phi.write();
//如果用户计算是，让输出p，即输入了writep，则输出p
if (args.options().found("writep"))
{
p.write();
}
//提示执行时间，cpu时间
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
//提示程序结束
Info<< "End\n" << endl;
//返回0
return(0);
}
例：OpenFOAM>>solver>>basic>>scalarTransportFoam
//createFields.H
//提示读入温度场
Info<< "Reading field T\n" << endl;
//温度场创建，标量场，需要初始化文件，下面各项具体含义，参看以前solver的说明
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//提示读入速度场
Info<< "Reading field U\n" << endl;
//速度场创建，各项意义，前面solver说明中已经给出
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//读入参数提示
Info<< "Reading transportProperties\n" << endl;
//根据名字为transportProperties的参数文件构建参数字典，以便查找。
IOdictionary transportProperties
(
IOobject
(
"transportProperties",//参数字典文件名字
runTime.constant(),//参数字典文件位置
mesh, //网格对象
IOobject::MUST_READ,//需要读入文件
IOobject::NO_WRITE //不对文件进行重写
)
);
//查询参数字典提示
Info<< "Reading diffusivity D\n" << endl;
//参数字典查询，初始化带单位标量DT（温度扩散率）
dimensionedScalar DT
(
transportProperties.lookup("DT")
);
//创建表面流率场，该文件位于
//src » finiteVolume » cfdTools » incompressible
# include "createPhi.H"
//scalarTransportFoam.C
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//主程序入口
int main(int argc, char *argv[])
{
//设置case目录相关，位于src » OpenFOAM » include
# include "setRootCase.H"
//创建time对象runTime，位于 src » OpenFOAM » include
# include "createTime.H"
//创建网格对象mesh，位于src » OpenFOAM » include
# include "createMesh.H"
//创建场对象，前面已经详述
# include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//提示计算标量传输方程
Info<< "\nCalculating scalar transport\n" << endl;
//显示当前courant数，位于src » finiteVolume » cfdTools » incompressible
# include "CourantNo.H"
//计算主流程
for (runTime++; !runTime.end(); runTime++)
{
//显示当前时间（物理时间，非cpu耗时）
Info<< "Time = " << runTime.timeName() << nl << endl;
//读入simple算法控制参数，位于src » finiteVolume » cfdTools » general » include
# include "readSIMPLEControls.H"
//网格非正交循环
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
//构造并求解方程
solve
(
fvm::ddt(T) //非稳态项，隐式离散
+ fvm::div(phi, T) //对流项，隐式离散
- fvm::laplacian(DT, T) //扩散项，隐式离散
);
}
runTime.write(); //求解结果输出，由于采用了注册机制，所有AUTO_WRITE声明的变量，
//都会输出
}
Info<< "End\n" << endl; //提示程序结束
return(0); //返回0
}
三个例子过后，我想聪明的你一定能在其中得到不少启发，正所谓举一反三。到最后不妨再仔细回想下其中的计算流程、以及过程的实现手段；作为一款优秀的CFD计算软件，大多数人喜欢它的原因正是这种高定制能力，如果你也因为其他流体计算软件没法灵活的实现所需功能，不妨仔细去查看下相关资料，当然了，这套软件能够解决的问题并不局限于流场计算。
