参考：

感谢原作者们的无私引路和宝贵工作。

前置：
OpenFOAM开发编程基础04 网格类初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

OpenFOAM 的计算绝大部分都是基于物理场的。基于前几次讨论，本文开始讨论 OpenFOAM 中场的基本操作。

为了操作方便，我们可以在 bashrc 中定义快捷命令（参考 OpenFOAM 环境准备 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)）

1
2
// bashrc
alias ofsp='cd /home/aerosand/aerosand/ofsp'

建立本文的项目文件夹并进入

1
2
3
4
// terminal
ofsp
mkdir 05_field
cd 05_field

场的创建

有单位量

dimensionSet

前文用到了 OpenFOAM 的单位预设，其定义在 $FOAM_SRC/OpenFOAM/dimensionSet/dimensionSets.C ，文件中还定义了很多预设单位

回忆 OpenFOAM 的单位设置为 [质量, 长度, 时间, 温度, 摩尔, 电流, 光强] 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
// $FOAM_SRC/OpenFOAM/dimensionSet/dimensionSets.C
...
const dimensionSet dimless;

const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0);
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0);
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0);
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0);
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0);
const dimensionSet dimCurrent(0, 0, 0, 0, 0, 1, 0);
const dimensionSet dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);

const dimensionSet dimArea(sqr(dimLength));
const dimensionSet dimVolume(pow3(dimLength));
const dimensionSet dimVol(dimVolume);

const dimensionSet dimVelocity(dimLength/dimTime);
const dimensionSet dimAcceleration(dimVelocity/dimTime);

const dimensionSet dimDensity(dimMass/dimVolume);
const dimensionSet dimForce(dimMass*dimAcceleration);
const dimensionSet dimEnergy(dimForce*dimLength);
const dimensionSet dimPower(dimEnergy/dimTime);

const dimensionSet dimPressure(dimForce/dimArea);
const dimensionSet dimCompressibility(dimDensity/dimPressure);
const dimensionSet dimGasConstant(dimEnergy/dimMass/dimTemperature);
const dimensionSet dimSpecificHeatCapacity(dimGasConstant);
const dimensionSet dimViscosity(dimArea/dimTime);
const dimensionSet dimDynamicViscosity(dimDensity*dimViscosity);
...

我们可以设置自定义变量的物理单位，如

1
dimensionSet dimAerosand(1,2,3,4,5,6,7); // just for example

dimensionedType

在 OpenFOAM 代码中可以看到如下类似语句

1
2
3
4
5
6
dimensionedScalar airPressure
(
	"airPressure",
	dimensionSet(1, -1, -2, 0, 0, 0, 0),
	1.0
);

查找源代码如下

1
2
// terminal
find $FOAM_SRC -iname dimensionedType.H

找到并摘取其中的构造函数如下

1
2
3
4
5
6
7
8
9
...
	//- Construct from components (name, dimensions, value).
	dimensioned
	(
		const word& name,
		const dimensionSet& dims,
		const Type& val
	);
...

OpenFOAM 中还提供了其他的构造函数，我们可以在不同情况下构造自己的有单位物理量。目前阶段来说，我们不需要继续深挖源代码。点到为止，明白确实可以这样用，并且以后按照这种方式去用就可以了。

当然，dimensionedType.H 也定义了配套的成员函数，比如

1
2
3
4
5
//- Return const reference to value.
const Type& value() const noexcept { return value_; }

//- Return non-const reference to value.
Type& value() noexcept { return value_; }

由此，我们可以便捷的使用 value() 成员函数，例如对上面的 airPressure 取值，再计算最大值

1
maxP = max(airPressure.value());

注意，OpenFOAM 会检查数学计算式等号两遍的最终单位是否相同，如果不同则会强制报错并中断程序，后续代码会面临并解决这个问题。

volFields

我们会注意到 OpenFOAM 中有典型的语句如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
...
    Info<< "Reading fieldp\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
...

那么 volScalarField 是什么呢？又是如何构造场的呢？

我们讨论一下场是怎么被创建的。

1
2
// terminal 
find $FOAM_SRC -iname volFieldsFwd.H

也许有同学疑惑“我怎么知道要查看这个文件”。参考 OpenFOAM开发编程基础03 时间类初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn) 第一节， 安装 OFextension 插件后，直接在 volScalarField 关键字上右键跳转即可。也可以打开 OpenFOAM Doxygen （OpenFOAM: User Guide: OpenFOAM®: Open source CFD : Documentation），查询相关关键字。除特殊情况，以后不再赘述。

volFieldFwd.H 中有类型定义如下

1
2
3
4
...
typedef GeometricField<scalar, fvPatchField, volMesh> volScalarField;
typedef GeometricField<vector, fvPatchField, volMesh> volVectorField;
...

可以看到，volFields 包括常见的 volScalarField 和 volVectorField，是体积场量。无论是 volScalarField 还是 volVectorField 其实都是 GeometricField 的模板的不同情况。所以， volFields 只是一种分类的称呼，而不是真正的类。

surfaceFields

同样的，面场量 surfaceFields 包括常见的 surfaceScalarField 和 surfaceVectorField 。

1
2
// terminal
find $FOAM_SRC -iname surfaceFieldsFwd.H

surfaceFieldFwd.H 中有类型定义如下

1
2
3
4
...
typedef GeometricField<scalar, fvsPatchField, surfaceMesh> surfaceScalarField;
typedef GeometricField<vector, fvsPatchField, surfaceMesh> surfaceVectorField;
...

可见，surfaceFields 也是 GeometricField 的模板的不同情况。

surfaceFields 的 surface 都可以写完整，为什么 volFields 的 volume 要简写成 vol 呢？也许是历史代码的原因吧。。。

GeometricField

进一步的，我们查看 GeometricField.H 文件

1
2
// terminal
find $FOAM_SRC -iname GeometricField.H

摘取类的主要结构和几个常用构造函数如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
...
template<class Type, template<class> class PatchField, class GeoMesh>
class GeometricField
:
    public DimensionedField<Type, GeoMesh>
{
public:
	...
private:
	...
public:
	...
	// Constructors

	//- Construct given IOobject, mesh, dimensioned<Type> and patch type.
	//  This assigns both dimensions and values.
	//  The name of the dimensioned\<Type\> has no influence.
	GeometricField
	(
		const IOobject& io,
		const Mesh& mesh,
		const dimensioned<Type>& dt,
		const word& patchFieldType = PatchField<Type>::calculatedType()
	);
	...
	//- Construct and read given IOobject
	GeometricField
	(
		const IOobject& io,
		const Mesh& mesh,
		const bool readOldTime = true
	);
	...
	//- Construct as copy resetting IO parameters
	GeometricField
	(
		const IOobject& io,
		const GeometricField<Type, PatchField, GeoMesh>& gf
	);
	...
...

上述构造函数，对应实际的场的构造，分别举例如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
...
	Info<< "Creating field light pressure\n" << endl;
	volScalarField lightP
	(
		IOobject
		(
			"lightP",
			runTime.timeName(),
			mesh,
			IOobject::NO_READ,
			IOobaject::AUTO_WRITE
		),
		mesh,
		dimensionedScalar("lightP",dimPressure,0.0)
		// 以构造的 dimensioned<Type>& 类型的变量初始化
		// 基于名称，单位，初始值
		// find $FOAM_SRC -iname dimensionedType.H
	);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
...
	Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
1
2
3
4
5
6
7
8
9
10
11
12
13
14
Info<< "Calculating field mixture density\n" << endl;
volScalarField rho
(
	IOobject
	(
		"rho",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE
	),
	alpha1*rho1 + alpha2*rho2
	// 传入计算场
);

这三个例子对应了最常见的场的构造函数，暂时明白这三个即可。其他构造函数以及更深入的代码层面的解析在以后会陆续介绍，无需担心。读者也可以尝试多阅读几个。

场的时间推进

我们以压力场为例，结合场的创建，讨论一下场的时间推进计算。

应用准备

1
2
3
4
5
// terminal 
foamNewApp 05_01_field
cd 05_01_field
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
code .

脚本和说明

略。

单时间步计算

主源码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"
    // 常见三个头文件，以后不再赘述

    volScalarField p // 基于 IOobject 和 mesh 创建压力场p
    (
        IOobject // 基于以下参数创建 IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE // 根据源代码中的设置自动写入
        ),
        mesh
    );
    // 读取可选：
	//     MUST_READ，必须读取，如果文件不存在则报错
	//     READ_IF_PRESENT，如果存在才读取，如果文件格式错误则报错
	//     NO_READ，不读取
	// 写入可选：
	//     AUTO_WRITE，根据要求自动写入
	//     NO_WRITE，不写入，但是可以通过代码显式写入

    runTime++; // 推进一个时间步
    // 如果省略此行，则不推进时间步，计算结果会写入到初始时间步，覆盖掉初始场

    forAll(p, i) // 遍历压力场
    {
        if (mesh.C()[i][1] < 0.5)
        // 还记得 mesh.C() 返回什么吗？
        // 返回网格中心点坐标的y坐标值（C++从0开始编号）
        {
            p[i] = i*runTime.value(); // 随便给一个可能的计算
        }
    }
    // 完成了整个压力场的计算

    Info<< "Max p = " << max(p).value() << endl; // 直接返回压力场最大值
    Info<< "Min p = " << min(p).value() << endl; // 返回压力场最小值
    Info<< "Average p = " << average(p).value() << endl; // 返回压力场平均值
    Info<< "Sum p = " << sum(p).value() << endl; // 返回压力场总和值

    p.write(); // 显式的代码，要求只写入该时间步的压力场值

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

    Info<< nl;
    runTime.printExecutionTime(Info);

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

    return 0;
}


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

编译运行

注意到主源码中的场的覆写是不区分初始时间步和后续时间步的，而且 foamCleanTutorials 命令不能恢复初始时间步文件夹内的场文件的内容。所以，为了避免初始条件被覆盖，强烈推荐备份初始条件。一旦初始场被覆盖，可以随时复原。操作如下

1
2
// terminal
cp -r debug_case/0 debug_case/0_orig

编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Create time

Create mesh for time = 0

Max p = 1.995
Min p = 0
Average p = 0.9975
Sum p = 399

ExecutionTime = 0 s  ClockTime = 0 s

End


查看 debug_case/ 文件夹下，多了 0.005/ 文件夹，其中的 0.005/p 文件中写入了计算的场。最后的输出就是最新时间步的场的值。

多时间步计算

如何计算多个时间步呢？

主源码

主源码修改如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    // runTime++;

    for (int i=0; i < 3; ++i) // 计算 3 个时间步
    {
        ++runTime; // 时间推进

        forAll(p, i)
        {
            if (mesh.C()[i][1] < 0.5)
            {
                p[i] = i*runTime.value();
            }
        }
        p.write(); // 写入该时间步
    }

    Info<< "Max p = " << max(p).value() << endl;
    Info<< "Min p = " << min(p).value() << endl;
    Info<< "Average p = " << average(p).value() << endl;
    Info<< "Sum p = " << sum(p).value() << endl;

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

    Info<< nl;
    runTime.printExecutionTime(Info);

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

    return 0;
}


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

关于写成 ++runTime 而不是 runTime++ 的讨论

编译运行

注意要先清理算例，避免错误，以后不再赘述。

1
2
3
// terminal
./caseclean
wclean

编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Create time

Create mesh for time = 0

Max p = 5.985
Min p = 0
Average p = 2.9925
Sum p = 1197

ExecutionTime = 0 s  ClockTime = 0 s

End


观察 debug_case/ 文件夹下多了三个时间步的文件夹，如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
|- debug_case/
	|- 0_orig/
		|- p
		|- U
	|- 0/
		|- p
		|- U
	|- 0.005/
		|- p
	|- 0.01/
		|- p
	|- 0.015/
		|- p
	|- constant/
	|- system/

查看各个时间步的压力场值，可以知道终端最后输出的值为最新时间步的计算值。

字典控制计算

即使是多时间步计算，我们也是直接修改主源码的计算时间步个数，每次修改都要重新编译整个应用。如果是大型应用，重新编译会相当麻烦。

在实际工作中，我们希望在应用编译完成后，也可以通过外部文件去控制计算，也就是使用 OpenFOAM 的字典文件 controlDict 控制计算。

主源码

主源码修改如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    while(runTime.loop()) // 时间推进
    {
        forAll(p, i)
        {
            if (mesh.C()[i][1] < 0.5)
            {
                p[i] = i*runTime.value();
            }
        }
		runTime.write(); // 按字典要求写入
    }

    Info<< "Max p = " << max(p).value() << endl;
    Info<< "Min p = " << min(p).value() << endl;
    Info<< "Average p = " << average(p).value() << endl;
    Info<< "Sum p = " << sum(p).value() << endl;

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

    Info<< nl;
    runTime.printExecutionTime(Info);

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

    return 0;
}


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

编译运行

清理调试算例，重新编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Create time

Create mesh for time = 0

Max p = 199.5
Min p = 0
Average p = 99.75
Sum p = 39900

ExecutionTime = 0 s  ClockTime = 0 s

End


调试算例 debug_case/ 文件多了几个时间步文件夹，分别是 0.1, 0.2, 0.3, 0.4, 0.5，写入的时间步正好由 system/controlDict 控制。

controlDict 字典的部分参数如下（每 0.005*20 = 0.1 写入一次）

1
2
3
4
5
6
7
...
deltaT          0.005;

writeControl    timeStep;

writeInterval   20;
...

这些参数都可以随时修改调试，无需重新编译应用。

场的基本操作

下面我们伪计算压力场、温度场、速度场。

  1. 求出网格单元和参考点的最大距离
  2. 根据时间、网格单元向量、参考向量、最大距离这 4 个参数计算压力场
  3. 根据压力梯度通过量纲处理，计算速度场

应用准备

1
2
3
4
// terminal
foamNewApp 05_02_field
cd 05_02_field
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case

脚本和说明

略。

增加物理场

拷贝过来的调试算例并没有温度场，我们给算例添加初始温度场。

温度和压力一样是标量，所以先拷贝压力场文件再修改。

1
2
3
// terminal
cp -r debug_case/0/p debug_case/0/T
code debug_case/0/T

温度场初始条件 debug_case/0/T，内容如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// dimensions      [0 2 -2 0 0 0 0];
dimensions              [0 0 0 1 0 0 0];

// internalField   uniform 0;
internalField           uniform 2023;

boundaryField
{
    movingWall
    {
        // type            zeroGradient;
        type            fixedValue;
        value           uniform 2049;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}


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

主源码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale);
// 自定义函数原型声明


int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    Info<< "Reading transportProperties\n" << endl;
    IOdictionary transportProperties // 构造字典对象，从字典文件中读入参数
    (
        IOobject
        (
            "transportProperties", // 字典文件名称
            runTime.constant(), // 字典文件路径
            mesh, // 和 mesh 建立关系
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );
    dimensionedScalar nu("nu",dimViscosity,transportProperties); // 读取字典文件中的参数

    Info<< "Reading field p\n" << endl;
    volScalarField p // 构造 p 场
    (
        IOobject // 基于 IOobject
        (
            "p", // 场文件的名称
            runTime.timeName(), // 场文件的路径
            mesh, // IOobject 和 mesh 建立关系
            IOobject::MUST_READ, // 必须要给初始条件
            IOobject::AUTO_WRITE
        ),
        mesh // 和 mesh 建立关系
    );

    Info<< "Reading field T\n" << endl;
    volScalarField T // 构造 T 场
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U // 构造 U 场
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

	// 举例创建空值场
	// 参考上文讨论的 GeometircField 构造函数
    volScalarField zeroScalarField
    (
        IOobject
        (
            "zeroScalarField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ, // 无需给初始条件的文件
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("zeroScalarField",dimless,Zero) // 初始化为零
    );

    volVectorField zeroVectorField
    (
        IOobject
        (
            "zeroVectorField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedVector("zeroVectorField",dimless,Zero) // 初始化为零
    );

    volVectorField gradT // 构造 gradT 场
    (
        IOobject
        (
            "gradT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::grad(T) // 因为是从 T 场计算而来，所以不需要再次和 mesh 建立关系
    );

    Info<< "Reading/calculating face flux field phi\n" << endl;
    surfaceScalarField phi // 构造 phi 场
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        fvc::flux(U)
        // fvc::interpolate(U)*mesh.Sf(); // 等同于此行计算
    );
    

    const vector originVector(0.05, 0.05, 0.005); // 构造带初值vector类的对象originVector

    const scalar rFarCell  = // 构造scalar类的对象rFarCell
    max // 求最大值
    (
        mag // 求数值大小
        (
            dimensionedVector("x0", dimLength, originVector) 
            // 以originVector创建dimensionedVector
            - mesh.C() // 减去网格单元中心点坐标，实际上也是向量
        )
    ).value(); // 返回值的大小

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

    while (runTime.loop()) // 时间推进
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        for (label cellI=0; cellI<mesh.C().size(); ++cellI) // 遍历所有网格单元
        {
            p[cellI] = calculatePressure // 通过自定义函数求 p 场
            (
                runTime.time().value(),
                mesh.C()[cellI],
                originVector,
                rFarCell
            );
        }

        U = fvc::grad(p) * dimensionedScalar("tmp",dimTime,1.0); // 求 U 场
        // 额外乘以一个含单位临时量 tmp ，保证两边结果单位相同

        runTime.write(); // 写出计算的场
    }

    Info<< "Calculation done." << endl;

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

    Info<< nl;
    runTime.printExecutionTime(Info);

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

    return 0;
}


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

// 自定义函数的实现
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale)
{
    scalar r(mag(x - x0) / scale); // 构造scalar类的对象r，并初始化

    scalar rR(1.0 / (r + 1e-12)); // 同上

    scalar f(1.0); // 同上

    return Foam::sin(2.0 * Foam::constant::mathematical::pi * f * t) * rR;
    // 返回一个计算值，并没有什么物理意义
}

编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field T

Reading field U

Reading/calculating face flux field phi


Starting time loop

Time = 0.005

Time = 0.01

Time = 0.015

Time = 0.02

...


ExecutionTime = 0.01 s  ClockTime = 0 s

End

使用 paraview 将计算结果可视化。注意打开需要路径

1
2
// terminal
paraFoam -case debug_case

注意到 debug_case/ 的时间步文件夹下新建了很多场文件，以 debug_case/0.3/ 为例，其文件结构如下

1
2
3
4
5
6
7
8
9
| - 0.3/
	|- uniform/
	|- gradT
	|- p
	|- phi
	|- T
	|- U
	|- zeroScalarField
	|- zeroVectorField

可以查看各个场文件中的数值。

注意：

打开 0.3/phi 可以看到并没有计算更新值，其他非 p,U 文件也没有计算更新值。这是因为场在创建时候给定的计算公式只是初始化，并不会随着时间推进而计算更新。所以，需要在主源码的时间循环中去计算需要的场。

代码整理

在 OpenFOAM 实际应用中，一般的

  • “场的接入”放入单独的文件 createFields.H
  • OpenFOAM 本身提供 createPhi.H
  • 自定义方法也放入单独文件，如 calculatePressure.H

所以该应用的文件结构调整后为

1
2
3
4
5
6
7
8
|- userApp/
	|- debug_case/
	|- Make/
		|- files
		|- options
	|- createFields.H
	|- calculatePressure.H
	|- 05_02_field.C

场接入文件 createFields.H 为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
Info<< "Reading transportProperties\n" << endl;
    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    dimensionedScalar nu("nu",dimViscosity,transportProperties);

    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field T\n" << endl;
    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField zeroScalarField
    (
        IOobject
        (
            "zeroScalarField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("zeroScalarField",dimless,Zero)
    );

    volVectorField zeroVectorField
    (
        IOobject
        (
            "zeroVectorField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedVector("zeroVectorField",dimless,Zero)
    );

    volVectorField gradT
    (
        IOobject
        (
            "gradT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::grad(T)
    );

    Info<< "Reading/calculating face flux field phi\n" << endl;
    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        fvc::flux(U)
        // fvc::interpolate(U)*mesh.Sf();
    );

自定义方法 calculatePressure.H 为

1
2
3
4
5
6
7
8
9
10
11
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale)
{
    scalar r(mag(x - x0) / scale);

    scalar rR(1.0 / (r + 1e-12));

    scalar f(1.0);

    return Foam::sin(2.0 * Foam::constant::mathematical::pi * f * t) * rR;
}

最终主源码整理为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include "fvCFD.H"

#include  "calculatePressure.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    #include "createFields.H"
    
    const vector originVector(0.05, 0.05, 0.005);

    const scalar rFarCell  = 
    max
    (
        mag
        (
            dimensionedVector("x0", dimLength, originVector)
            - mesh.C()
        )
    ).value();

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

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        for (label cellI=0; cellI<mesh.C().size(); ++cellI)
        {
            p[cellI] = calculatePressure
            (
                runTime.time().value(),
                mesh.C()[cellI],
                originVector,
                rFarCell
            );
        }

        U = fvc::grad(p) * dimensionedScalar("tmp",dimTime,1.0);

		phi = fvc::flux(U); // 增加 phi 的计算

        runTime.write();
    }

    Info<< "Calculation done." << endl;

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

    Info<< nl;
    runTime.printExecutionTime(Info);

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

    return 0;
}


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

整理后的这种形式更接近 OpenFOAM 求解器等应用的代码组织形式。

  • 主源码功能划分非常清晰
  • 每个功能都更加易于维护

编译运行结果当然还是一样的。

小结

回顾这 5 篇讨论，我们大概可以感受到数值计算的核心要素所在——时间、网格、物理场。有了这些对象之后，我们便可以组建线性代数方程，进行数值求解。不要着急，我们距离求解器编程越来越近了。