参考：

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

前置：
OpenFOAM基础算法 SIMPLE & PISO & PIMPLE | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

前面已经讨论了定常标量输运方程的求解、非定常波动方程的求解。如果我们想进一步学习，可能在网上很多教程都推荐从 icoFoam 开始学习，但是当我们打开 icoFoam 求解器的源代码，会发现手足无措，陌生的代码和算法把我们蒙在了窗户纸后面，阻碍我们前进。

事实上，初学者在看 OpenFOAM 标准求解器之前应该有两件事要做

  1. 理解算法 OpenFOAM算法 SIMPLE & PISO & PIMPLE | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)
    • 这篇算法讨论其实已经写了不少代码的实现部分
  2. 直观的 SIMPLE 算法实现
    • 没有任何 SIMPLE 以外功能的代码

所以我们仍然不急于看 icoFoam 求解器或者其他标准求解器。下面我们不考虑任何的优化、也不考虑收敛检查等等额外的功能，单从 SIMPLE 算法出发，直观地在 OpenFOAM 中实现该算法的求解器。

控制方程如下

∇⋅U=0(1)\nabla\cdot U = 0 \tag{1}

∂U∂t+∇⋅UU−∇⋅(μ∇U)=−∇p(2)\frac{\partial U}{\partial t} + \nabla\cdot UU - \nabla\cdot(\mu\nabla U) = -\nabla p \tag{2}

建立本文项目文件夹

1
2
3
4
// terminal
ofsp
mkdir 09_simple
cd 09_simple

应用准备

1
2
3
4
5
// terminal
foamNewApp 09_01_simple
cd 09_01_simple
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
cp /home/aerosand{caserun,caseclean} .

脚本和说明

略。

求解器

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
/*
* Contents
*   # Basic fields
*   # Transport properties
*/

/*
* # Basic fields
*/

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


Info<< "Creating old layer pressure field p_old\n" << endl;
volScalarField p_old
(
    IOobject
    (
        "p_old",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    p
);
// 因为压力场要做松弛处理，所以需要复制一份 p 场

#include "createPhi.H"

/*
* # Transport properties
*/

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

dimensionedScalar nu
(
    "nu",
    dimensionSet(0, 2, -1, 0, 0, 0, 0),
    transportProperties
);

主源码

主源码中应用 SIMPLE 算法，实现如下

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

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

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

    #include "createMesh.H"

    #include "createFields.H" // 创建场、输运参数
    
    // 读取代数求解器参数
    Info<< "Reading fvSolutino\n" << endl;
	IOdictionary fvSolution
	(
	    IOobject
	    (
	        "fvSolution",
	        runTime.system(),
	        mesh,
	        IOobject::MUST_READ_IF_MODIFIED,
	        IOobject::NO_WRITE
	    )
	);

	// scalar alpha(readScalar(fvSolution.lookup("alpha")));
	dimensionedScalar alphaD("alpha",dimless,fvSolution);
	scalar alpha(alphaD.value());
	
	// scalar pRefCell(readScalar(fvSolution.lookup("pRefCell")));
	dimensionedScalar pRefCellD("pRefCell",dimless,fvSolution);
	scalar pRefCell(pRefCellD.value());
	
	// scalar pRefValue(readScalar(fvSolution.lookup("pRefValue")));
	dimensionedScalar pRefValueD("pRefValue",dimless,fvSolution);
	scalar pRefValue(pRefValueD.value());

	// 输出刚才读取的代数求解器参数
    Info<< "fvSolution parameters" << nl
        << "\trelaxation factor \"alpha\": " << alpha << nl
        << "\tindex of cell containing reference pressure \"pRefCell\": " << pRefCell << nl
        << "\treference pressure value \"pRefValue\": " << pRefValue << nl
        << endl;

	// 时间推进
    Info<< "\nStarting time loop\n" << endl;

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

        // 组建动量方程
        fvVectorMatrix UEqn
        (
            fvm::div(phi,U)
            - fvm::laplacian(nu,U)
            == 
            - fvc::grad(p)
        );

		// 求解动量方程
        UEqn.solve();


        // 组建压力方程
        volScalarField A = UEqn.A();
        volVectorField H = UEqn.H();

        volScalarField rA = 1.0 / A;

        volVectorField HbyA = H * rA;

        fvScalarMatrix pEqn
        (
            fvm::laplacian(rA,p)
            ==
            fvc::div(HbyA)
        );

        pEqn.setReference(pRefCell,pRefValue); // 为什么要设参考值？

		// 求解压力方程
        pEqn.solve();

		// 欠松弛处理
        p = alpha*p + (1.0-alpha)*p_old;

		// 更新速度场
        U = rA*H - rA*fvc::grad(p);

		// 求解质量通量，可以直接使用 flux() 函数
        phi = fvc::interpolate(U) & mesh.Sf();
		// phi = fvc::flux(U);

		// 边界条件修正，暂不深究
        U.correctBoundaryConditions();
        p.correctBoundaryConditions();

		// 更新旧时间步场量，用于下一次欠松弛处理
        p_old = p;

        runTime.write();
    }

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

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

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

    return 0;
}


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

参考压力解释可见 关于fluent中的压力（一） - 希望先生 - 博客园 (cnblogs.com)

编译

1
2
3
// terminal 
wclean
wmake

调试算例

  • system/blockMeshDict 无需修改，我们仍然使用原几何模型和几何网格
  • system/controlDict 无需修改，我们仍然使用求解控制的设置（步长步数等）
  • decomposedParDict 和 PDRblockMeshDict 不用管，本文用不到
  • constant/transportProperties.H 无需修改，我们沿用粘度参数
  • debug_case/0/ 初始条件无需修改

fvSchemes

注意各个格式需要至少提供默认值，例如

1
2
3
4
divSchemes
{
    default         Gauss linear;
}

fvSolution

最后添加语句

1
2
3
4
5
alpha       0.01;

pRefCell    0;

pRefValue   0;

controlDict

自行尝试更改计算参数

计算

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

paraFoam -case _case

可以看到计算结果。

整理代码

我们可以把代码整理为如下文件结构

1
2
3
4
5
6
7
8
9
10
11
12
|- userApp/
	|- debug_case/
	|- Make/
		|- files/
		|- options/
	|- 09_01_simple.C
	|- createFields.H
	|- readfvSolution.H
	|- pEqn.H
	|- UEqn.H
	|- caserun
	|- caseclean

主源码

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

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

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

    #include "createMesh.H"

    #include "createFields.H"
    #include "readfvSolution.H"

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

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

        #include "UEqn.H"

        #include "pEqn.H"

        runTime.write();
    }

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

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

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

    return 0;
}


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

readfvSolution.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
Info<< "Reading fvSolutino\n" << endl;
IOdictionary fvSolution
(
    IOobject
    (
        "fvSolution",
        runTime.system(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

// scalar alpha(readScalar(fvSolution.lookup("alpha")));
dimensionedScalar alphaD("alpha",dimless,fvSolution);
scalar alpha(alphaD.value());

// scalar pRefCell(readScalar(fvSolution.lookup("pRefCell")));
dimensionedScalar pRefCellD("pRefCell",dimless,fvSolution);
scalar pRefCell(pRefCellD.value());

// scalar pRefValue(readScalar(fvSolution.lookup("pRefValue")));
dimensionedScalar pRefValueD("pRefValue",dimless,fvSolution);
scalar pRefValue(pRefValueD.value());

Info<< "fvSolution parameters" << nl
    << "\trelaxation factor \"alpha\": " << alpha << nl
    << "\tindex of cell containing reference pressure \"pRefCell\": " << pRefCell << nl
    << "\treference pressure value \"pRefValue\": " << pRefValue << nl
    << endl;

UEqn.H

1
2
3
4
5
6
7
8
9
10
fvVectorMatrix UEqn
(
    fvm::div(phi,U)
    - fvm::laplacian(nu,U)
    == 
    - fvc::grad(p)
);

UEqn.solve();

pEqn.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
volScalarField A = UEqn.A();
volVectorField H = UEqn.H();

volScalarField rA = 1.0 / A;

volVectorField HbyA = H * rA;

fvScalarMatrix pEqn
(
    fvm::laplacian(rA,p)
    ==
    fvc::div(HbyA)
);

pEqn.setReference(pRefCell,pRefValue);

pEqn.solve();

p = alpha*p + (1.0-alpha)*p_old; // 压力松弛

U = rA*H - rA*fvc::grad(p);

phi = fvc::interpolate(U) & mesh.Sf();

U.correctBoundaryConditions();
p.correctBoundaryConditions();

p_old = p;

编译计算

仍然使用上面的调试算例，结果是相同的。

小结

我们通过百分百还原数学方程，没有加入其他任何功能，最后我们得到了一个单纯基于 SIMPLE 算法的求解器。