参考：

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

前置：
OpenFOAM开发编程基础06 开发库 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

通过前面总共七篇讨论，我们已经大概熟悉了 OpenFOAM 中数值应用的架构，熟悉了数值计算的关键元素（时间、网格、场），也熟悉了一些编程的语法、一些工作流程。下面我们使用所掌握的工具，尝试写一个最简单的 OpenFOAM 求解器（不涉及数值算法）。

回忆 CFD理论基础00 基本方程 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn) 。

计算流体力学的通用基本方程为

∂∂t(ρϕ)+∇⋅(ρUϕ)=∇⋅(Γ∇ϕ)+Sϕ\frac{\partial}{\partial t}(\rho \phi) + \nabla \cdot (\rho U\phi) = \nabla\cdot(\Gamma\nabla\phi) + S_{\phi}

为了避免未知量和 OpenFOAM 中的质量通量 phi 误会，我们使用 A 表示待求物理场。

下面，我们通过 OpenFOAM 编程求解一个不可压定常无源输运问题，控制方程如下

∇⋅(UA)=∇⋅(Γ∇A)\nabla\cdot(UA) = \nabla\cdot(\Gamma\nabla A)

建立本文项目文件夹

1
2
3
4
// terminal
ofsp
mkdir 07_1stSol
cd 07_1stSol

应用准备

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

脚本和说明

略。

求解器

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
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )
);

Info<< "Reading diffusivity\n" << endl;
dimensionedScalar gamma("gamma",dimViscosity,transportProperties);


Info<< "Reading field A\n" << endl;
volScalarField A
(
    IOobject
    (
        "A",
        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
);

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

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

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

    #include "createMesh.H"

    #include "createFields.H"

    ++runTime;
    // 推进一个时间步，避免计算结果覆盖初始时间步

    solve // 求解控制方程
    (
        fvm::div(phi,A)
        - fvm::laplacian(gamma,A)
    );

    A.write(); // 写出计算值 场A


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

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

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

    return 0;
}

方程的构建

上面的主源码咋一看好像也不难理解，但仔细想一想，读者可能会忍不住好奇，方程的构建是怎么样的？

需要强调的是，想要完全理解这些问题，毫无疑问还是要从有限体积法基础开始，深入到矩阵计算中去。强烈建议没有阅读过有限体积法基础的同学还是快速找补一下，避免基础不稳困惑加倍，参考 CFD理论基础02 有限体积法初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)。

我们看一下代码中的方程构建，可以按下面方式构建这个方程（结果和上述的主源码没有什么区别）

1
2
3
4
5
6
7
fvScalarMatrix AEqn
{
	 fvm::div(phi,A)
	- fvm::laplacian(gamma,A)
	==
	0
};

对流项

对流项的构造是 fvm::div(phi,A) ，这是一个函数，我们查找这个函数

1
2
3
// terminal
find $FOAM_SRC -iname fvmdiv.c
// /usr/lib/openfoam/openfoam2306/src/finiteVolume/lnInclude/fvmDiv.h

为了帮助理解，从代码声明 fvmDiv.H 中摘取几处代码简单讨论如下（不建议继续深究代码，代码挖的太深无益于主要学习的推进）

1
2
3
4
5
6
7
8
9
10
11
12
...
namespace Foam // 属于命名空间 Foam
{
namespace fvm // 属于命名空间 fvm
...
    template<class Type> // 模板
    tmp<fvMatrix<Type>> div // 返回类型 + 函数名
    (
        const surfaceScalarField&, // 形参1 surfaceScalarField 类型的引用
        const GeometricField<Type, fvPatchField, volMesh>& // 形参2 GeometricField 类型的应用，回忆 GeometricField 有不同的类型定义
    );
...

我们了解了此函数的形参类型，明白返回的是 fvMatrix 类型的变量，数学上对应的是和待求量相关的矩阵。

扩散项

扩散项的构造是 fvm::laplacian(gamma,A) ，这也是一个函数，我们查找这个函数

1
2
3
// terminal
find $FOAM_SRC -iname fvmlaplacian.H
// /usr/lib/openfoam/openfoam2306/src/finiteVolume/lnInclude/fvmLaplacian.H

同样的，摘部分代码如下

1
2
3
4
5
6
7
8
9
10
11
namespace Foam // 属于命名空间 Foam
{
namespace fvm // 属于命名空间 fvm
...
    template<class Type, class GType> // 模板
    tmp<fvMatrix<Type>> laplacian // 返回类型 + 函数名
    (
        const dimensioned<GType>&, // 形参1 dimensioned 类型的引用
        // dimensioned 有不同类型的定义
        const GeometricField<Type, fvPatchField, volMesh>& // 形参2 同样是 GeometricField 类型的引用 
    );

我们知道该函数同样返回 fvMatrix 类型，对应数学上和待求量相关的矩阵。

构建

离散方程最后形式为

Ax=bAx = b

对流项和扩散项的矩阵相加，构成了待求左侧项，因为没有源项，所以右侧项置空为零，对应的代码就是

1
2
3
4
 fvm::div(phi,A)
- fvm::laplacian(gamma,A)
==
0

求解

构建完成后，我们使用 solve 函数求解线性代数系统

1
solve(AEqn);

暂时不用翻它的源代码。

其他疑惑

为什么散度计算中，使用的是 phi 而不是公式中的 U 呢？此外， fvm 和常见到的 fvc 有什么区别呢？

简单理解的话， 求解器中的 fvm 离散所有对系数矩阵有贡献的项（隐式离散），相对比的，fvc 离散和系数矩阵无关的项（显式离散）。

但是这么轻飘飘的一句话解释，其实很难让初学者理解。为了讨论上面两个问题，我们需要回到控制方程的偏微分方程数值求解上，从有限体积法开始，一步一步推导，直到矩阵求解，才可以真正明白为什么要使用通量进行计算，为什么会有隐式离散和显式离散的区分。这里先不深究，参见 CFD 基础系列讨论和 OpenFOAM 求解器系列讨论（ofss）。

好在求解语句的代码语法已经十分贴近数学公式，即使带着这些困惑也不要紧，我们仍然可以编写简单的求解器。

编译

1
2
// termianl
wmake

调试算例

需要对调试算例做一些修改，以便计算进行。

  • system/blockMeshDict 无需修改，我们仍然使用原几何模型和几何网格
  • system/controlDict 无需修改，我们仍然使用求解控制的设置（步长步数等）
  • decomposedParDict 和 PDRblockMeshDict 不用管，初学阶段暂时用不到

fvSchemes

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
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,A)      Gauss linear;
}

laplacianSchemes
{
    default         none;
    laplacian(gamma,A)  Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

  • 因为求解涉及到散度和拉普拉斯，所以需要指明两个计算项的离散格式

至于这些格式是在计算过程的哪里起作用的，怎么在起作用的，很难简单解释清楚，还是要参见 CFD 基础系列讨论和 ofss 系列讨论

fvSolution

只计算了 场 A

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    A
    {
        solver          GAMG;
        smoother        DILUGaussSeidel;
        tolerance       1e-06;
        relTol          0.05;
    }
}

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

  • 因为求解器只计算了场 A，所以只设置 A 的代数求解器

初始条件

在调试算例的 p 文件上修改文件名称得到 A 文件，内容修改如下

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
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      A;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    movingWall
    {
        type            zeroGradient;
    }

    fixedWalls
    {
        type            fixedValue;
        value           uniform 1;
    }

    frontAndBack
    {
        type            empty;
    }
}


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

U 文件保持不变。

运行

1
2
// terminal
./caserun

终端显式如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Create time

Create mesh for time = 0

Reading transportProperties

Reading diffusivity

Reading field A

Reading field U

Reading/calculating face flux field phi

GAMG:  Solving for A, Initial residual = 1, Final residual = 1.1539e-06, No Iterations 1

ExecutionTime = 0 s  ClockTime = 0 s

End


调试算例下多了新的一个时间步文件夹，其中有场 A 的计算结果。

通过 paraview 打开显式结果（以后不再赘述命令）

1
2
// termianl
paraFoam -case debug_case

可以看到一个时间步的计算结果。

没有时间步动画显示，这是当然的，因为这是定常问题，所以我们只计算一个时间步（仅仅是为了不覆盖初始条件）。

小结

我们终于写了第一个求解器，但是还有一些疑惑涉及到有限体积法的数值细节，其他部分的代码和以前讨论里的应用其实没有多少区别。

我们记录一下可能有的疑惑：

  1. fvSchemes 文件中对离散格式的指定是怎么反映到求解计算中的？
  2. fvSolution 文件中对求解方法的指定是怎么反映到求解计算中的？

这两个问题可以当成黑箱去设置调用，具体细节暂时放一放，需要参考 CFD 基础系列和 ofss 系列，不要担心，到时候会一个一个解决。

  1. 代码实现中为什么使用质量通量而不是直接使用速度？
  2. 所谓最后 solve 的线性代数系统在数学上是怎么组建的？fvm 和 fvc 在其中有什么区别？

这两个可能的疑惑请记下来并在后续学习中逐渐解决。如前所述，好在 OpenFOAM 的语法简明直观，即使带着这些困惑也可以实现简单的求解器。

我们下一篇讨论再求解一个简单的非定常问题。