OpenFOAM求解器开发

CFD C++ OpenFOAM

---

1. OpenFOAM简介

OpenFOAM（Open Source Field Operation and Manipulation 的缩写，意为 开源的场运算和处理 软件）是对连续介质力学问题进行数值计算的C++自由软件工具包，其代码遵守GNU通用公共许可证。它可进行数据预处理、后处理和自定义求解器，常用于CFD领域。[1]

LOGO：
[2]

---

2. OpenFOAM使用方式

使用OpenFOAM进行CFD模拟，大致可分为 三种 类型：

• 直接利用OpenFOAM的标准求解器，替代商业软件。优点：免费，缺点：交互性差，对使用者要求较高；

• 自定义求解器，利用OpenFOAM的基本类库来按照自己的需求来编写针对某类应用的求解器。用户并不需要特别关心离散和求解的最底层的知识，如时间项、对流项的离散等，关注的重点是物理问题本身。在编程中，通常是顶层的求解流程的开发，不需要用户深入研习C++语言。商业软件中的所谓UDF和OpenFOAM是不能相提并论的；

• 自定义离散方法，属于更高级的应用，如果用户更关注时间项ddt，扩散项laplacian，对流项div是如何离散的，能否有更高效更高精度的离散方法，这需要修改finiteVolume库和OpenFOAM库中对应的代码。尤其是对流项，尽管OpenFOAM已经提供了基于NVD和TVD的模板和40多种有名的高阶高精度格式，但可以预见，这仍然是不够的，毕竟对流项的离散仍然是目前CFD的重点研究方向。

可以肯定的是，目前有很多人关注 第二种 类型的应用，毕竟将OpenFOAM当成Fluent或Star-CCM来使用，并不见得方便。但是将OpenFOAM作为类库来构建自己的求解器，这是其它软件无法实现的。

---

3. 基本知识

几个重要的 环境变量:

环境变量名	值	说明
$FOAM_TUTORIALS	OpenFOAM算例目录	包含OpenFOAM官方教程算例
$FOAM_SRC	OpenFOAM总库源代码目录	finiteVolume、mesh等库源代码目录
$FOAM_APP	OpenFOAM应用源代码目录	包含solvers、test和utilities
$FOAM_APPBIN	OpenFOAM标准应用目标代码目录	官方求解器和工具应用目录
$FOAM_USER_APPBIN	OpenFOAM用户自定义应用目标代码目录	用户自定义求解器和工具应用目录
$FOAM_RUN	用户算例目录	算例目录

重要的 shell命令：

OpenFOAM专有命令	等价系统命令	说明
run	cd $FOAM_RUN	打开$FOAM_RUN目录
src	cd $FOAM_SRC	打开$FOAM_SRC目录
app	cd $FOAM_APP	打开$FOAM_APP目录
util	cd $FOAM_APP/utilities	打开$FOAM_APP/utilities目录
sol	cd $FOAM_APP/solvers	打开$FOAM_APP/solvers目录
tut	cd $FOAM_TUTORIALS	打开$FOAM_TUTORIALS目录

---

4. 求解器基本架构

基于OpenFOAM的求解器以一种规范的格式来架构，每个程序的源代码放置在以这个程序newApp命名的文件夹中，同时这个文件夹必须包含一个名为Make的文件夹，用于存放编译选项，整个架构的基本文件结构如下：

newApp
├──newApp.C
├──createFields.H
└──Make
├──files
└──options

其中：

文件	用途
newApp.C	newApp求解器顶层源代码
createFields.H	变量场的声明和初始化
files	按行存储所有源代码文件名，最后一行用来指定目标代码EXE的名称和存放位置
options	设定查找头文件和库的路径EXE_INC和需要链接的库EXE_LIBS

---

5. 编译

5.1 编译配置文件 files 解析

files 文件基本内容为：

newApp.C
EXE = $(FOAM_APPBIN)/newApp

其中，通过关键字$(FOAM_APPBIN)和newApp分别指定目标代码的存放位置和名称，标准路径$(FOAM_APPBIN)亦可替换为用户路径$(FOAM_USER_APPBIN)。

5.2 编译配置文件 options 解析

options 文件基本内容为：

EXE_INC = \
    -I<directoryPath1> \
    -I<directoryPath2> \
    ... \
    -I<directoryPathN>
EXE_LIBS = \
    -L<libraryPath1> \
    -L<libraryPath2> \
    ... \
    -L<libraryPathN> \
    -l<library1> \
    -l<library2> \
    ... \
    -l<libraryN>

其中，包含的头文件路径使用-I标识符指定，库文件路径使用-L标识符指定，库名称使用-l标识符指定，注意区分。此外，每个文件夹路径或文件名的前面都必须有标识符，在EXE_INC和EXE_LIBS之后以及每个项目后面需要使用\，最后的项目后面没有\。

5.3 wmake编译系统

由于每一段代码都需要跟某些独立的OpenFOAM库来进行对接，为简化编译，OpenFOAM提供了建立在make之上wmake程序脚本，用于执行 维护依赖文件列表 和 编译源代码 的任务。

5.3.1 头文件包含

在wmake编译系统中，编译器按照以下顺序在某个文件夹下寻找被包含的 头文件：

顺序	路径	备注
1	$WM_PROJECT_DIR/src/OpenFOAM/lnInclude文件夹
2	本地lnInclude文件夹	例如newApp/lnInclude
3	本地文件夹	例如newApp
4	$WM_PROJECT_DIR/wmake/rules/-$WM_ARCH文件中设定的依赖文件路径
5	在options文件中通过-I指定的其它文件夹

5.3.2 库文件链接

在wmake编译系统中，编译器在以下路径链接 库文件：

顺序	路径
1	$FOAM_LIBBIN目录
2	$WM_DIR/rules/$WM_ARCH文件中设定的依赖文件路径
3	在options中指定的其它目录

对于链接的 单个 库文件，必须通过标示符-l指定，并且去掉lib前缀以及.so后缀，例如：libnew.so，应该使用-lnew。默认情况下， wmake 调用下面的库文件：

顺序	路径
1	$FOAM_LIBBIN文件夹下的libOpenFOAM.so库
2	$WM_DIR/rules/$WM_ARCH文件中设定的依赖库路径
3	在options中指定的其它目录

5.3.3 运行wmake

运行wmake的命令具有如下结构：

$$\color{teal}{\verb+wmake +}\color{fuchsia}{\verb+<optionalArguments> +}\color{maroon}{\verb+<optionalDirectory>+}$$

其中$\color{maroon}{\verb+<optionalDirectory>+}$用于指定被编译程序的文件路径，我们一般在程序自己的路径下编译源代码，因此$\color{maroon}{\verb+<optionalDirectory>+}$可以省略。

一般不需要指定参数$\color{fuchsia}{\verb+<optionalArguments>+}$的值，但是在编译类型为 库 等其他类型而不是可执行目标代码的情况下，需要指定参数$\color{fuchsia}{\verb+<optionalArguments>+}$的值：

$\color{fuchsia}{\verb+<optionalArguments>+}$	编译类型
lib	静态链接库
libso	动态链接库
libo	静态链接目标文件
jar	JAVA存档文件
exe	可执行文件(默认)

5.3.4 与wmake相关的环境变量

OpenFOAM中与wmake编译系统相关的环境变量具有前缀$WM_，具体如下表(以64位机器上安装的OpenFOAM-2.4.0为例)所示：

环境变量名	值(默认)	说明
$WM_PROJECT	OpenFOAM	工程名
$WM_PROJECT_VERSION	2.4.0	工程版本号
$WM_PROJECT_INST_DIR	$HOME/OpenFOAM	OpenFOAM安装目录总路径
$WM_PROJECT_USER_DIR	$HOME/OpenFOAM/$USER-2.4.0	OpenFOAM用户目录路径
$WM_PROJECT_DIR	$HOME/OpenFOAM/OpenFOAM-2.4.0	OpenFOAM工程目录路径
$WM_ARCH	linux64	主机架构
$WM_ARCH_OPTION	64	机器位数:32位或64位
$WM_COMPILER	Gcc	wmake使用的编译器类型
$WM_COMPILER_LIB_ARCH	64	编译器库架构
$WM_DIR	$HOME/OpenFOAM/OpenFOAM-2.4.0/wmake	wmake目录路径
$WM_MPLIB	SYSTEMOPENMPI	并行库
$WM_PRECISION_OPTION	DP	编译精度:DP双精度或SP单精度
$WM_COMPILE_OPTION	Opt	编译优化开关:Debug或Opt
$WM_OPTIONS	linux64GccDPOpt	wmake编译设置汇总

其中，$WM_OPTIONS = $WM_ARCH + $WM_COMPILER + $WM_PRECISION_OPTION + Opt。一般与wmake相关的环境变量 无需更改。

如果需要临时更改编译优化设置以获得详细的Bug信息，可在终端输入以下命令临时更改编译优化选项：

$ export WM_COMPILE_OPTION=Debug

5.3.5 依赖包的移除

在编译的过程中，wmake会创建一个扩展名为.dep的依赖包文件（newApp.dep），并在 Make/$WM_OPTIONS中产生一系列文件。如果源代码修改后，用户想要删除这些文件，可以运行wclean命令删除，其具有如下结构：

$$\color{teal}{\verb+wclean +}\color{fuchsia}{\verb+<optionalArguments> +}\color{maroon}{\verb+<optionalDirectory>+}$$

其中，与wmake相同，$\color{maroon}{\verb+<optionalDirectory>+}$用于指定被编译程序的文件路径，我们一般在程序自己的路径下运行wclean，因此$\color{maroon}{\verb+<optionalDirectory>+}$可以省略。如果用户想要删除依赖文件和Make目录中的其它文件，就不需要$\color{maroon}{\verb+<optionalDirectory>+}$。 然而，如果$\color{maroon}{\verb+<optionalDirectory>+}$中指定了lib，那么本地lnInclude文件也会被删除。

另一个工具是$\color{teal}{\verb+rmdepall+}$，它从执行点开始自上而下移除所有的.dep文件，可用于升级OpenFOAM库。

---

6. 实例

以标准laplacianFoam求解器为例，介绍OpenFOAM求解器开发的基本步骤。

6.1 基本控制方程

laplacianFoam求解器用于求解热传递温度控制方程：

$$\begin{equation} \dfrac{\partial T}{\partial t} - \nabla^2 \left( D_T T \right) = 0, \end{equation}$$

其中，$T$ [K]表示温度，$t$ [s]时间，$\nabla^2$拉普拉斯算子（laplacian），$D_T$ [m²/s]热扩散系数，由下式计算得到：

$$D_T = \dfrac{k}{\rho c_p}$$

其中，$k$ [W/(m·K) = J/(m·K·s)]表示热传导系数，$\rho$ [kg/m³]密度，$c_p$ [J/(kg·K)]比热容。

6.2 拷贝源代码并清理相关编译支持文件

以下操作均在 终端（Terminal） 中完成。

• 进入laplacianFoam求解器所在目录

$ sol
$ cd basic

• 将laplacianFoam文件夹整个拷贝至$WM_PROJECT_USER_DIR并将其重命名为myLaplacianFoam

$ cp -r laplacianFoam $WM_PROJECT_USER_DIR/myLaplacianFoam

• 进入myLaplacianFoam目录

$ cd $WM_PROJECT_USER_DIR/myLaplacianFoam

• 重命名laplacianFoam.C为myLaplacianFoam.C

$ mv laplacianFoam.C myLaplacianFoam.C

• 删除write.H文件，此文件用于输出温度场$T$的梯度，这里不需要。

$ rm write.H

• 执行清理操作

$ wclean

至此，myLaplacianFoam求解器的文件结构为

myLaplacianFoam
├──myLaplacianFoam.C
├──createFields.H
└──Make
├──files
└──options

6.3 createFields.H 文件解析

createFields.H 文件用于创建温度场$T$并初始化，读取物理参数热扩散系数$D_T$。文件内容如下：

    //提示读入温度场T
    //Info等价于C++中std::cout,C++标准库的开发晚于OpenFOAM,info语句也因此得以出现并被保留至今
    Info<< "Reading field T\n" << endl;
    //volScalarField用于创建名为T的体心(变量T存储于控制体中心)标量场
    volScalarField T
    (
        //OpenFOAM采用IOobject定义输入输出类,但凡定义的场大多需要此类定义,底层代码可先不深究
        IOobject
        (
            //场的名称(T),用于确定初始文件名称(T)以及变量名称(T)
            "T",
            //存储位置为运行时间(文件夹)
            //该位置由$case/system/controlDict中的startTime控制,位于$case/startTime
            runTime.timeName(),
            //注册于网格对象mesh中
            mesh,
            //该对象通过读取文件创建,必须进行读取
            //如果某个场通过计算得到,可以不进行读取,例如通量phi
            IOobject::MUST_READ,
            //根据控制字典文件(controlDict)中的设置自动(AUTO_WRITE)写入结果文件
            IOobject::AUTO_WRITE
        ),
        //新建体心标量场(T)所用的网格对象,在createMesh.H创建
        mesh
    );
    //提示读入参数控制文件
    Info<< "Reading transportProperties\n" << endl;
    //IOdictionary用于新建字典文件    
    //参数控制文件声明通过文件(transportProperties)读取
    IOdictionary transportProperties
    (
        IOobject
        (
            //文件名(transportProperties)
            "transportProperties",
            //文件位置,位于$case/constant
            runTime.constant(),
            //网格对象,主要从事对象注册,以便由runTime.write()控制输出
            mesh,
            //在字典文件被更改的时候进行读取
            IOobject::MUST_READ_IF_MODIFIED,
            //不输出,字典文件不需要输出
            IOobject::NO_WRITE
        )
    );
    //提示读入热扩散系数DT
    Info<< "Reading diffusivity DT\n" << endl;
    //dimensionedScalar用于新建带单位标量,新建扩散率DT,因为各向同性,DT为常数
    //通过查询参数控制文件($case/constant/transportProperties),初始化带单位标量
    dimensionedScalar DT
    (
        //查询关键字(DT)并读取
        transportProperties.lookup("DT")
    );

6.4 myLaplacianFoam.C 文件解析

myLaplacianFoam.C 文件为myLaplacianFoam求解器顶层源代码文件。文件内容（略去文件头部许可声明及简介）如下：

//fvCFD头文件
//包含大多数CFD计算需要的头文件
//涉及到时间构建,矩阵组建,有限体积离散,网格组建,量纲设置等大量内容,为自定义求解器必备头文件
#include "fvCFD.H"
//simpleControl头文件
//定义SIMPLE循环,使用SIMPLE循环必备头文件
#include "simpleControl.H"
//主程序入口
int main(int argc, char *argv[])
{
    //setRootCase头文件
    //根据输入参数argc和argv设置算例根目录rootcase,必备头文件
    #include "setRootCase.H"
    //createTime头文件
    //创建时间对象,涉及到runTime控制,非定常求解器必备头文件
    #include "createTime.H"
    //createMesh头文件
    //创建网格对象,根据$case/constant/polyMesh文件夹中的网格数据创建对象mesh,必备头文件
    #include "createMesh.H"
    //创建场对象(T),位于求解器根目录
    #include "createFields.H"
    //从网格mesh对象构造类simpleControl(基类:solutionControl)的对象simple
    simpleControl simple(mesh);
    //提示计算温度分布    
    Info<< "\nCalculating temperature distribution\n" << endl;
    //开始时间循环,采用SIMPLE算法必备语句
    while (simple.loop())
    {
        //提示当前时间
        Info<< "Time = " << runTime.timeName() << nl << endl;
        //是否进行非正交修正
        //如果在fvSolution字典文件中设置为0,就只求解控制方程一次,如果设置为n,则求解控制方程n-1次        
        while (simple.correctNonOrthogonal())
        {
            //求解方程(1)
            //solve是FOAM名称空间的全局函数,参数为矩阵(fvMatrix)
            //fvm表示隐式离散,返回有限体积稀疏矩阵类fvMatrix对象
            solve
            (
                //隐式离散(fvm关键字)时间项和Laplace项
                fvm::ddt(T) - fvm::laplacian(DT, T)     
            );
        }
        //#include "write.H"    //注释或者删除本行
        runTime.write();        //添加本行用来输出T标量场
        //提示求解器执行时间及CPU耗时
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    //提示求解器执行结束
    Info<< "End\n" << endl;
    return 0;
}

再次说明一下，注释或者删除行：

#include "write.H"

同时添加：

runTime.write();

6.5 编译准备

6.5.1 修改配置文件 files

修改配置文件 files 内容如下：

myLaplacianFoam.C
EXE = $(FOAM_USER_APPBIN)/myLaplacianFoam

6.5.2 修改配置文件 options

无需修改配置文件 options，其内容如下：

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

6.6 编译

以下操作均在 终端（Terminal） 中完成。

• 进入myLaplacianFoam目录并执行编译

$ cd $WM_PROJECT_USER_DIR/myLaplacianFoam
$ wmake

• 进入$FOAM_USER_APPBIN目录查看编译好的可执行文件myLaplacianFoam

$ cd $FOAM_USER_APPBIN
$ ls

• 查看自定义myLaplacianFoam求解器的基本信息

$ myLaplacianFoam -help

显示如下信息

至此，编译完成。

---

7. 算例验证

7.1 物理问题描述

问题物理域和初始条件如下图所示：

7.2 网格

这里使用OpenFOAM二维多面体网格准备中生成的多边形棱柱体网格，网格如下图所示：

7.3 算例文件结构

一个标准的OpenFOAM算例文件结构如下：

$case 算例根目录
├──constant 网格和输运参数目录
│ ├──polyMesh 网格目录，网格格式详见CFD多面体网格数据结构—OpenFOAM
│ │ ├──points
│ │ ├──faces
│ │ ├──owner
│ │ ├──neighbour
│ │ └──boundary
│ └──transportProperties 输运参数设置
├──0 初始和边界条件目录，每个待求解变量都需要一个单独文件设定其初始和边界条件
│ └──T 变量T初始和边界条件设置
└──system 求解器设置目录
├──controlDict 计算控制参数设置，设置起始终止时间，时间步长，输出控制等
├──fvSchemes 微分算子离散格式设置
└──fvSolution 代数方程组求解和算法设置

7.4 各个设置文件具体内容

7.4.1 transportProperties 文件

transportProperties 文件内容（略去文件头部注释）如下：

//文件说明
FoamFile
{
    version     2.0;                    //版本号
    format      ascii;                  //文本格式
    class       dictionary;             //类型为字典文件
    location    "constant";             //所在目录
    object      transportProperties;    //对象名
}
//设置热扩散系数DT的单位和值
DT              DT [ 0 2 -1 0 0 0 0 ] 0.1;

其中，[ 0 2 -1 0 0 0 0 ]用于设置DT的单位，依次表示质量[kg]、米[m]、时间[s]、开尔文[K]、摩尔质量[kgmol]、电流[A]和光强[cd]的幂。热扩散系数DT的单位为[m²/s]，因此分别设置[m]和[s]的幂为2和-1，其余均为0。

7.4.2 T 文件

T 文件内容（略去文件头部注释）如下：

//文件说明
FoamFile
{
    version     2.0;                //版本号
    format      ascii;              //文本格式
    class       volScalarField;     //场类型
    object      T;                  //对象名
}
dimensions      [0 0 0 1 0 0 0];    //设置单位
internalField   uniform 0.0;        //内部场初始化为0.0
//设置边界条件
boundaryField
{
    //2D算例,back和front面边界条件均设置为empty
    back
    {
        type            empty;
    }
    //同back面
    front
    {
        type            empty;
    }
    //固定值1.0
    top
    {
        type            fixedValue;
        value           uniform 1.0;
    }
    //固定值0.0
    left
    {
        type            fixedValue;
        value           uniform 0.0;
    }
    //固定值0.0
    bottom
    {
        type            fixedValue;
        value           uniform 0.0;
    }
    //固定值0.0
    right
    {
        type            fixedValue;
        value           uniform 0.0;
    }
}

7.4.3 controlDict 文件

controlDict 文件内容（略去文件头部注释）如下：

//文件说明
FoamFile
{
    version     2.0;                //版本号
    format      ascii;              //文本格式
    class       dictionary;         //类型为字典文件
    location    "system";           //所在目录
    object      controlDict;        //对象名
}
application     myLaplacianFoam;    //求解器名称
//时间控制部分
startFrom       startTime;          //设置计算开始时间为startTime
startTime       0;                  //定义startTime
stopAt          endTime;            //设置计算结束时间为endTime
endTime         4;                  //定义endTime
deltaT          0.005;              //设置时间步长
//计算结果输出控制部分
writeControl    runTime;            //按照计算时间控制输出
writeInterval   0.05;               //输出时间间隔,单位秒[s],与writeControl配合使用
purgeWrite      0;                  //输出是否覆盖先前记录,0为不覆盖,1可用于稳态计算,默认0
writeFormat     ascii;              //计算结果输出文件格式指定
writePrecision  10;                 //控制输出有效数字位数,与writeFormat配合使用,默认6
writeCompression uncompressed;      //是否压缩
timeFormat      general;            //时间目录名称指定
timePrecision   6;                  //时间目录名称有效数字位数,与timeFormat配合使用,默认6
//字典读取控制
runTimeModifiable yes;              //每个时间步是否读取所修改的字典文件设定

各个参数和可选取值详见OpenFOAM中 controlDict 字典文件解析。

7.4.4 fvSchemes 文件

fvSchemes 文件内容（略去文件头部注释）如下：

//文件说明
FoamFile
{
    version     2.0;                //版本号
    format      ascii;              //文本格式
    class       dictionary;         //类型为字典文件
    location    "system";           //所在目录
    object      fvSchemes;          //对象名
}
//时间一阶导项(非定常项)格式
ddtSchemes
{
    default         Euler;
}
//梯度项格式
gradSchemes
{
    default         Gauss linear;
    grad(T)         Gauss linear;
}
//散度项格式
divSchemes
{
    default         none;
}
//拉普拉斯项格式
laplacianSchemes
{
    default         none;
    laplacian(DT,T) Gauss linear corrected;
}
//插值格式
interpolationSchemes
{
    default         linear;
}
//面法向梯度格式
snGradSchemes
{
    default         corrected;
}
//需要用来计算通量的场
fluxRequired
{
    default         no;
    T;
}

各个子字典和格式详见OpenFOAM中 fvSchemes 字典文件解析。

7.4.5 fvSolution 文件

fvSchemes 文件内容（略去文件头部注释）如下：

//文件说明
FoamFile
{
    version     2.0;                //版本号
    format      ascii;              //文本格式
    class       dictionary;         //类型为字典文件
    location    "system";           //所在目录
    object      fvSolution;         //对象名
}
//代数方程组求解设置
solvers
{
    T
    {
        solver          PCG;        //求解器
        preconditioner  DIC;        //预处理器
        tolerance       1e-08;      //残值
        relTol          0;          //相对残值
    }
}
//SIMPLE算法设置
SIMPLE
{
    nNonOrthogonalCorrectors 4;     //非正交修正次数
}

各个子字典和格式详见OpenFOAM中 fvSolution 字典文件解析。

7.5 计算结果

计算得到的温度场T云图如下所示：

---

感谢您的阅读，欢迎讨论和批评指正。

作者：戴得志