OpenFOAM学习笔记_01_icoFoam理解_openfoam icofoam

参考资料：

1. 东岳博文：http://dyfluid.com/icoFoam.html
2. 牛奕博文：http://wap.sciencenet.cn/blog-3410930-1175782.html?mobile=1
3. cloudbird07博文：https://blog.csdn.net/cloudbird07/article/details/105441128
4. OpenFOAMwiki： https://openfoamwiki.net/index.php/IcoFoam
5. OpenFOAM C++ Source Code Guide：https://cpp.openfoam.org/v8/icoFoam_8C_source.html

OpenFOAM所提供的icoFoam求解器，所求解的是非定常不可压缩层流流动问题（连续方程+不可压缩Navier-Stokes方程），其使用PISO算法来解耦速度与压力。

1 控制方程

$$\begin{gathered} \nabla \cdot \mathbf{u}=0 \\ \frac{\partial \mathbf{u}}{\partial t}+\nabla \cdot (\mathbf{u}\mathbf{u})=-\nabla \left(\frac{p}{\rho }\right)+\nabla \cdot (\nu \nabla \mathbf{u}) \end{gathered}$$

由于不可压缩问题中密度$\rho$为常数，可以直接将其从方程中去除，即认为$\rho =1$，将$p/\rho$计作$p$。运动粘性$\nu =\mu /\rho$，同样不包含密度$\rho$。对流项$\nabla \cdot (\mathbf{u}\mathbf{u})$为非线性项，通常将头一个$\mathbf{u}$固化为上个时刻的${\mathbf{u}}^n$。

2 算法流程

PISO算法流程如下：

a. 离散和求解动量方程来获取新速度场${\mathbf{u}}^{\star }$，注意压力梯度项已经从源项中剥离出来；

$$a_C^{\mathbf{u}}{\mathbf{u}}_C+\sum _{F\sim NB(C)}{a}_F^{\mathbf{u}}{\mathbf{u}}_F={\mathbf{b}}_C^{\mathbf{u}}-\nabla p_C^{(n)}$$

b. 基于新速度场${\mathbf{u}}_C^{\ast }$，构建${\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A}}_C^{\star }$；

$${\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A}}_C^{\star }=\frac{1}{a_C^{\mathbf{u}}}\left({\mathbf{b}}_C^{\mathbf{u}}-\sum _{F\sim NB(C)}{a}_F^{\mathbf{u}}{\mathbf{u}}_F^{\star }\right)$$

c. 求解Possion方程得到新的压力场

$$\nabla \cdot ({\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A}}^{\star })=\nabla \cdot \left(\frac{1}{a_C^{\mathbf{u}}}\nabla p^{\star }\right)$$

d. 根据新的压力场$p^{\star }$，更新界面通量${\mathbf{u}}_f^{\star }$和速度场${\mathbf{u}}_C^{\star }$；

$$\begin{gathered} {\mathbf{u}}_f^{\star }={\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A}}_f^{\star }-\frac{1}{a_f^{\mathbf{u}}}\nabla p_f^{\star } \\ {\mathbf{u}}_C^{\star }={\mathbf{H}\mathbf{b}\mathbf{y}\mathbf{A}}_C^{\star }-\frac{1}{a_C^{\mathbf{u}}}\nabla p_C^{\star } \end{gathered}$$

e. 重复b-d直至收敛，然后进行下一时间步的推进直至终止时间。

3 代码分析

/*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
     \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
 
     OpenFOAM is free software: you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by
     the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
 
     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
 
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
     icoFoam
 
 Description
     Transient solver for incompressible, laminar flow of Newtonian fluids.
 
 \*---------------------------------------------------------------------------*/
 
 // fvCFD.h头文件中包含了许多计算所需要的头文件，里面定义了大量类和函数
 // pisoControl.h头文件中则定义了pisoControl类，用于控制piso算法流程，
 // 有time-loop和piso-loop信息。
 #include "fvCFD.H"
 #include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
     // setRootCaseLists.H头文件中设置算例case的根root目录信息
     // createTime.h中定义了时间类Time的对象runTime
     // createMesh.h中定义了网格类fvMesh的对象mesh
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createMesh.H"
 
     // 定义piso算法控制类pisoControl的对象piso，用mesh构造
     pisoControl piso(mesh);
 
     // createFields.H中定义粘性标量nu、体心标量场压力p、体心矢量场速度U、面心标量场通量phi，
     // 以及指标型参考压力单元pRefCell和标量型参考压力值pRefValue等。
     #include "createFields.H"
     // initContinuityErrs.H中定义标量型累计连续误差cumulativeContErr，初始化为0。
     #include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     // 输出“开始时间循环了！”
     Info<< "\nStarting time loop\n" << endl;
 
     // 时间步逐个推进
     while (runTime.loop())
     {
         // 输出当前时刻
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // 计算Courant数
         #include "CourantNo.H"
 
         // Momentum predictor
         // a. 动量预测，u_C^star
 
         // 离散动量方程，获得矢量方程UEqn，注意压力梯度项并未包含其中。
         // 这里的fvm表明均为隐式处理，即离散到线性代数方程的左端项系数矩阵中，以及右端项中。
         fvVectorMatrix UEqn
         (
             fvm::ddt(U)
           + fvm::div(phi, U)
           - fvm::laplacian(nu, U)
         );
 
         // 如果有动量预测步骤的话，求解动量方程，用所得速度来作为预测速度，
         // 不然的话直接用上一时刻的速度作为预测速度（即省略了步骤a，最开始不求动量方程了）
         // 注意，在求解动量方程的时候把压力梯度项-fvc::grad(p))显式拿了进来
         if (piso.momentumPredictor())
         {
             solve(UEqn == -fvc::grad(p));
         }
 
         // --- PISO loop
         // PISO循环的压力速度解耦迭代过程
         while (piso.correct())
         {
             // b. 构建HbyA_C^star
             // 体心标量场rAU，存储着动量方程对角线系数的倒数，即1/a_C^U
             volScalarField rAU(1.0/UEqn.A());
             // 体心矢量场HbyA，UEqn.H()为动量方程离散后的右端项（不含压力梯度项）
             // 减去 左端非对角元系数与预测速度的乘积，rAU*UEqn.H()则再除以对角元系数
             // constrainHbyA(..., U, p)则再对HbyA的边界值进行一定的修正。
             volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
             // 面形心标量场phiHbyA，为HbyA的面通量，
             // ddtPhiCorr项通过提取插值速度与面通量的差值，引入了面速度场的散度，
             // 第二项可能是做的某些修正，不是很清楚……
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
                 fvc::flux(HbyA)
               + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
             );
 
             // 调整进口和出口通量使其符合连续性，这对压力解的存在性，很有帮助。
             adjustPhi(phiHbyA, U, p);
 
             // Update the pressure BCs to ensure flux consistency
             // 更新压力边界条件来确保通量一致性。
             constrainPressure(p, U, phiHbyA, rAU);
 
             // Non-orthogonal pressure corrector loop
             // 非正交压力修正循环
             // Laplacian的非正交部分是由最近求得的压力值来计算的，使用延迟修正方法引入。
             // 这里多说两句：
             // 对于非正交网格而言，扩散项离散的时候会有交叉扩散项的引入，其一般采用当次迭代的
             // 变量值计算，并当做源项处理，即延迟修正引入，再进行2-3次迭代求解以使得该非正交
             // 修正充分引入，同时注意不要引入过多的非正交修正，其既不会收敛又耗费计算时间。
             // 还有就是即便是正交网格，下面的循环体部分也会执行一次，即至少求解一次该压力方程。
             while (piso.correctNonOrthogonal())
             {
                 // Pressure corrector
                 // c. 压力修正
 
                 // 组装压力修正方程，获得标量方程pEqn
                 fvScalarMatrix pEqn
                 (
                     fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                 );
 
                 // 在不可压缩流动中，只有压力的相对值才是有效的，所以需要设置某个单元为
                 // 压力参考单元，其压力值为参考压力值，来确保解的唯一性。
                 pEqn.setReference(pRefCell, pRefValue);
 
                 // 求解压力修正方程
                 pEqn.solve();
 
                 // d. 更新界面通量
                 // 如果是最后一步非正交修正，那么修正界面通量phi值，即u_f^star修正。
                 if (piso.finalNonOrthogonalIter())
                 {
                     phi = phiHbyA - pEqn.flux();
                 }
             }
 
             #include "continuityErrs.H"
 
             // d. 更新速度，修正速度u_C^star
             U = HbyA - rAU*fvc::grad(p);
             // 修正速度场u的边界条件
             U.correctBoundaryConditions();
         }
 
         // 视情况输出
         runTime.write();
 
         // 显示计算时间、计算耗时
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }
 
     // 显示终止
     Info<< "End\n" << endl;
 
     return 0;
 }
 
 
 // ************************************************************************* //
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