参考：

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

前置：
OpenFOAM开发编程基础07 第一个求解器 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

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

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

求解的波动方程如下

∂2A∂t2=∇⋅(c2∇A)\frac{\partial^2 A}{\partial t^2} = \nabla \cdot (c^2 \nabla A)

我们依然用 A 表示待求的物理场。

新建本文的项目文件夹

1
2
3
4
// terminal
ofsp
mkdir 08_waveSol
cd 08_waveSol

应用准备

1
2
3
4
5
6
7
// terminal
foamNewApp 08_01_waveSol
cd 08_01_waveSol
cp /home/aerosand{caserun,caseclean} .
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity .
cp -r debug_case/0 debug_case/0_orig
code .

脚本和说明

略。

求解器

createFields.H

场文件 /userApp/createFields.H 接入 方程所需要的场 A 和参数 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
/*
* Contents
*   # Basic fields
*   # Transport properties
*/


/*
* # Basic fields
*/

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


/*
* # Transport properties
*/

Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )
);

Info<< "Reading c from transportProperties" << endl;
dimensionedScalar c("c",dimViscosity,transportProperties);

主源码

文件 /userApp/08_01_waveSol.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
#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    argList::addNote(
        "This application solves the wave equation over given mesh."
    );

    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    #include "createFields.H"

    Info<< nl << "Starting time loop..." << endl;

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

        Info<< "Solving h field" << endl;

        fvScalarMatrix hEqn
        (
            fvm::d2dt2(h)
            == fvm::laplacian(sqr(C), h)
        );

        hEqn.solve();

        runTime.write();

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

    }

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


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

    return 0;
}

编译

1
2
3
// termianl
wclean
wmake

调试算例

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

transportProperties

文件 debug_case/constant/transportProperties 内容修改如下

1
C              0.05;

初始条件

基于初始条件 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
36
FoamFile
{
    version     2.0;
    format      ascii;
    arch        "LSB;label=32;scalar=64";
    class       volScalarField;
    location    "0";
    object      A;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform 0;
    }
    fixedWalls
    {
        type            fixedValue;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
    }
}


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

setFieldsDict

使用命令创建一个模板到调试算例下

1
2
// terminal
foamGetDict -case debug_case setFieldsDict

修改文件 /debug_case/system/setFieldsDict ，内容如下

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

defaultFieldValues
(
    volScalarFieldValue A 0
);

regions
(
    // Set cell values
    // (does zerogradient on boundaries)
    boxToCell
    {
        box (0.04 0.04 0) (0.06 0.06 0.01);

        fieldValues
        (
            volScalarFieldValue A 1
        );
    }
);

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

  • 注意调试算例的 blockMeshDict 中设置了 scalar 0.1; ，几何尺寸有了缩放，所以在设置 box 的时候要注意坐标点
  • 设置 boxToCell 中 box 的两个角点坐标
  • 注意，setFields 只重设了初始场，区别于场里的 源 。

fvSchemes

增加时间项二阶偏导的数值计算离散格式，其他不变

1
2
3
4
5
...
d2dt2Schemes
{
    default         Euler;
}

fvSolution

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

solvers
{
    A
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    AFinal
    {
        $p;
        relTol          0;
    }
}

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

也可以写成 A|AFinal ，然后一起设置

controlDict

随意设置计算时间长短等。

脚本

因为额外使用了 setFields 命令，且该命令会更改初始条件，所以需要修改脚本。

caserun 内容如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#!/bin/sh
cd "${0%/*}" || exit 1                              # Run from this directory
#------------------------------------------------------------------------------

appPath=$(cd `dirname $0`; pwd)                     # Get application path
appName="${appPath##*/}"                            # Get application name

# cd debug_case

# echo $appName

blockMesh -case debug_case | tee debug_case/log.mesh
echo "Meshing done."

cp -r debug_case/0_orig/ debug_case/0/
setFields -case debug_case
echo "Setting fileds done."

$appName -case debug_case | tee debug_case/log.run

caseclean 内容如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#!/bin/sh
cd "${0%/*}" || exit 1                              # Run from this directory
#------------------------------------------------------------------------------

appPath=$(cd `dirname $0`; pwd)
appName="${casePath##*/}"

cd debug_case

rm -rf log.*
foamCleanTutorials
rm -rf 0/
echo "Cleaning done."

计算

1
2
3
4
5
// terminal
./caseclean
./caserun

paraFoam -case debug_case

使用 paraview 可以看到中心点场随着时间的波动变化。

小结

本文讨论求解了一个非定常波动问题。我们通过注释和代码拆出，尝试把求解器结构写的更加清晰易读。首次使用了 OpenFOAM 自带的工具 setFields （本质上也是一种 应用）。

需要注意的是，因为待求的波动方程简单，所以本文求解器对偏微分方程的求解并没有使用任何的算法，只是简单的在每个时间步求解。实际开发中，更多的问题是多个物理场耦合，所以我们需要一些数值算法来求解偏微分方程组，而不是单个偏微分方程。

虽然本系列是开发编程系列，但是在正式进入 OpenFOAM 标准求解器开发前，我们在下一次讨论中不得不讲一下 OpenFOAM 的基础求解算法，主要是 SIMPLE, PISO, PIMPLE 。

虽然系列讨论追求给初学者带来“连续”“逐深”“有方向”“不间断”的学习体验，但这是避不开的一环。