参考：

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

前置：
OpenFOAM开发编程基础05 场的基本操作 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

OpenFOAM 中，当某一些类实现某个特定的功能，可以把它们开发成库。

对“类”“库”链接有疑问的，请参考 OpenFOAM开发编程基础00 基本实现和开发 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

我们将上篇讨论中的压力场、速度场的伪计算开发成独立的库，然后在应用中调用它，当然开发库也可以被任何应用调用。

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

1
2
3
4
// terminal
ofsp
mkdir 06_devLib
cd 06_devLib

应用准备

1
2
3
4
5
// terminal
foamNewApp 06_01_devLib
cd 06_01_devLib
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
cp -r debug_case/0 debug_case/0_orig

在应用文件夹中新建开发库

1
2
// terminal
foamNewApp computeVelocityPressure

脚本和说明

略。

文件结构

整个应用的文件结构如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
|- 06_01_devLib/
	|- computeVelocityPressure/
		|- Make/
			|- files
			|- options
		|- computeVelocityPressure.H
		|- computeVelocityPressure.C
	|- debug_case/
		|- 0/
		|- constant/
		|- system/
	|- Make/
		|- files
		|- options
	|- createFields.H
	|- 06_01_devLib.C
	|- caserun
	|- caseclean
	|- README.md

开发库

头文件

头文件 /userApp/computeVelocityPressure/computeVelocityPressure.H ，内容如下

1
2
3
4
5
6
7
8
9
#include "fvCFD.H"

scalar computeR(const fvMesh& mesh, volScalarField& r, dimensionedVector x0);

volScalarField computePressure(const fvMesh& mesh, scalar t, dimensionedVector x0, scalar scale,scalar f);

void computeVelocity(const fvMesh& mesh, volVectorField& U, word pName = "p");
// 第三个参数给了默认值

  • 各个函数的原型
  • 为什么开发库的头文件要包含 fvCFD.H ？读者可以尝试注释掉这行，重新编译，并阅读报错信息。这样可以加深对 fvCFD.H 文件内容的理解。

源代码

源文件 /userApp/computeVelocityPressure/computeVelocityPressure.C，我们采用更“OpenFOAM”风格的实现方式，内容如下

注意和 OpenFOAM开发编程基础05 场的基本操作 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn) 比较各个函数的实现

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
#include "computeVelocityPressure.H"

scalar computeR(const fvMesh& mesh, volScalarField& r, dimensionedVector x0)
{
    r = mag(mesh.C() - x0); // 场 r 保存各个单元网格中心点向量到参考向量的距离

    return max(r).value(); // 返回各个距离的最大值
}

volScalarField computePressure // 计算压力场，注意返回类型
(
    const fvMesh& mesh, // 传入 mesh 的引用，引用传递效率高
    scalar t, 
    dimensionedVector x0, // 不仅是 vector，而且是有单位制的 vector
    scalar scale, // 纯 scalar
    scalar f // 纯 scalar
)
{
    volScalarField r(mag(mesh.C()-x0)/scale); 
    // 上一篇讨论中的 r 只是 scalar 类型，这里是 volScalarField 类型
    // 注意计算中变量类型的统一

    volScalarField rR(1.0/(r+dimensionedScalar("tmp",dimLength,1e-12)));
    // 赋予 1e-12 单位，以通过OpenFOAM对计算的单位检查
    // 注意计算中结果单位的统一

    return Foam::sin(2.0*Foam::constant::mathematical::pi*f*t)*rR*dimensionedScalar("tmp",pow(dimLength,3)/pow(dimTime,2),1.0);
    // 单位可以按数学方式计算
    // 同样为了保证结果单位一致，额外乘以含单位临时量 tmp
    // 该函数返回的计算结果的单位，等于，主函数要接收该值的变量的单位
}

void computeVelocity(const fvMesh& mesh, volVectorField& U, word pName)
// 注意函数的实现不能再写上形参的默认值
{
    const volScalarField& pField = mesh.lookupObject<volScalarField>(pName);
    // 按 pName 字段查找场量，并赋值给 pField

    U = fvc::grad(pField) * dimensionedScalar("tmp",dimTime,1.0);
    // 计算 U，保证单位统一
}

也许会有同学对比数据类型有困惑，比如 scalar 类型和 volScalarField 类型。简单理解就是，scalar 只是一个值，而 volScalarField 是一组数据，包含了计算域内所有离散点的 scalar 值，相当于是一个矩阵，也就是表示了一个“场”。

库 Make

/userApp/computeVelocityPressure/Make/files 内容如下

1
2
3
4
computeVelocityPressure.C

LIB = $(FOAM_USER_LIBBIN)/libcomputeVelocityPressure

  • 注意 LIB 关键字的变化

/userApp/computeVelocityPressure/Make/options

1
2
3
4
5
6
7
8
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

  • 该库，强调是这个库，在编译的时候，不会使用到其他更多的库，所以 options 中有这些基础库就够了，不需要修改

编译

终端使用命令

1
2
3
4
// terminal
wmake computerVelocityPressure/

// wclean computerVelocityPressure/
  • 这个库也是独立的，可以被其他任何应用调用，只要写对调用路径

createFields.H

创建场 /userApp/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
113
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();
    );

主源码

主源码 /userApp/06_01_devLib.C，内容如下

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
#include "fvCFD.H"

#include "computeVelocityPressure.H"

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

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

    #include "createMesh.H"

    #include "createFields.H"

    dimensionedVector originVector("x0",dimLength,vector(0.05,0.05,0.005));
    // 构造 dimensionedVector 类的对象 originVector，并赋单位和初始值

    scalar f(1.0);

    volScalarField r // 只是计算的过程量，所以没有放进 createFields.H 中
    (
        IOobject
        (
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("r",dimLength,0.0)
    );

    const scalar rFarCell = computeR(mesh,r,originVector);
    // 计算得到距离的最大值，并传给 rFarCell 变量保存

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

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

        p = computePressure(mesh,runTime.time().value(),originVector,rFarCell,f);
        // 计算 p 场

        computeVelocity(mesh,U);
        // 计算 U 场

        phi = fvc::flux(U);
        // 计算 phi 场

        runTime.write(); // 写出场
    }

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

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

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

    return 0;
}


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

应用 Make

/userApp/Make/files

1
2
3
4
06_01_devLib.C

EXE = $(FOAM_USER_APPBIN)/06_01_devLib

/userApp/Make/options

1
2
3
4
5
6
7
8
9
10
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -IcomputeVelocityPressure/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -L$(FOAM_USER_LIBBIN) -lcomputeVelocityPressure

  • 注意开发库的处理，既要“包含”，也要“链接”。

编译运行

1
2
wmake
./caserun

运行结果如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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

...

ExecutionTime = 0.01 s  ClockTime = 0 s

End

我们也可以通过 paraview 打开查看计算结果。

1
2
// terminal
paraFoam -case debug_case

小结

研究中涉及到的通用方法可以抽象成开发库，但是初学者不能追求把方法都写成开发库，具体问题应当具体对待。通过本篇讨论，最主要的还是熟悉代码架构和代码语法。

此外，编程中注意单位统一，注意数据类型统一。