第二章：CFD算例之后向台阶流
后向台阶流为一个二维的、等温的、不可压缩流动。其计算域由以下特点组成：
左边存在一个进口；
右边存在一个出口；
上下均为边界；
流体从左侧进口以10米/秒的速度流入。在OpenFOAM-11中，将调用incompressibleFluid模块。
前处理
OpenFOAM所有自带算例都存储在$FOAM_TUTORIALS文件夹。用户需要首先将相应的算例，复制到一个文件夹下。在本算例中，我们将自带算例复制到$FOAM_RUN文件夹下。用户需要分3次运行下述代码：
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady .
cd pitzDailySteady
在输入第一行代码的时候，如果输出：
dyfluid@dyfluid-virtual-machine:~$ cd $FOAM_RUN
bash: cd: /home/dyfluid/OpenFOAM/dyfluid-11/run: No such file or directory
用户需要手动创立一个文件夹：
cd /home/dyfluid/OpenFOAM
mkdir dyfluid-11
cd dyfluid-11
mkdir run
注意：
上方代码中的dyfluid为你的用户名。如果你的ubuntu系统的用户名为pangmao，则需要把相应的dyfluid改为pangmao。
然后执行
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady .
cd pitzDailySteady
即可。这样在pitzDailySteady会发现3个文件夹。所有OpenFOAM算例文件夹都主要分为3个，一个是0文件夹，一个是constant文件夹，一个是system文件夹。在0文件夹里，可以设置算例的初始条件。在constant文件夹里，可以设置算例的物理性质。在system文件夹，主要设置算例的算法控制。该算例下还存在一个Allrun文件。其主要用于自动处理脚本。
网格生成
OpenFOAM使用blockMesh来进行网格生成。blockMesh在生成网格的时候需要指定特别复杂的点位置信息。因此其只能用于生成一些基本的网格。当然了，其可以生成非常复杂的网格，但在这种情况下需要花费大量的时间去写点位置信息。OpenFOAM中的blockMesh主要是结合snappyHexMesh来做网格。snappyHexMesh为OpenFOAM内置的另外一个网格生成程序。将会在以后介绍。本算例采用blockMesh来生成一个简单的后向台阶流网格。
下面是blockMesh的字典文件。其名称为blockMeshDict，其位于system文件夹下：
1/*--------------------------------*- C++ -*----------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration | Website: https://openfoam.org
5 \\ / A nd | Version: dev
6 \\/ M anipulation |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10 format ascii;
11 class dictionary;
12 object blockMeshDict;
13}
14// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
15
16// Note: this file is a Copy of $FOAM_TUTORIALS/resources/blockMesh/pitzDaily
17
18convertToMeters 0.001;
19
20vertices
21(
22 (-20.6 0 -0.5)
23 (-20.6 25.4 -0.5)
24 (0 -25.4 -0.5)
25 (0 0 -0.5)
26 (0 25.4 -0.5)
27 (206 -25.4 -0.5)
28 (206 0 -0.5)
29 (206 25.4 -0.5)
30 (290 -16.6 -0.5)
31 (290 0 -0.5)
32 (290 16.6 -0.5)
33
34 (-20.6 0 0.5)
35 (-20.6 25.4 0.5)
36 (0 -25.4 0.5)
37 (0 0 0.5)
38 (0 25.4 0.5)
39 (206 -25.4 0.5)
40 (206 0 0.5)
41 (206 25.4 0.5)
42 (290 -16.6 0.5)
43 (290 0 0.5)
44 (290 16.6 0.5)
45);
46
47negY
48(
49 (2 4 1)
50 (1 3 0.3)
51);
52
53posY
54(
55 (1 4 2)
56 (2 3 4)
57 (2 4 0.25)
58);
59
60posYR
61(
62 (2 1 1)
63 (1 1 0.25)
64);
65
66
67blocks
68(
69 hex (0 3 4 1 11 14 15 12)
70 (18 30 1)
71 simpleGrading (0.5 $posY 1)
72
73 hex (2 5 6 3 13 16 17 14)
74 (180 27 1)
75 edgeGrading (4 4 4 4 $negY 1 1 $negY 1 1 1 1)
76
77 hex (3 6 7 4 14 17 18 15)
78 (180 30 1)
79 edgeGrading (4 4 4 4 $posY $posYR $posYR $posY 1 1 1 1)
80
81 hex (5 8 9 6 16 19 20 17)
82 (25 27 1)
83 simpleGrading (2.5 1 1)
84
85 hex (6 9 10 7 17 20 21 18)
86 (25 30 1)
87 simpleGrading (2.5 $posYR 1)
88);
89
90boundary
91(
92 inlet
93 {
94 type patch;
95 faces
96 (
97 (0 1 12 11)
98 );
99 }
100 outlet
101 {
102 type patch;
103 faces
104 (
105 (8 9 20 19)
106 (9 10 21 20)
107 );
108 }
109 upperWall
110 {
111 type wall;
112 faces
113 (
114 (1 4 15 12)
115 (4 7 18 15)
116 (7 10 21 18)
117 );
118 }
119 lowerWall
120 {
121 type wall;
122 faces
123 (
124 (0 3 14 11)
125 (3 2 13 14)
126 (2 5 16 13)
127 (5 8 19 16)
128 );
129 }
130 frontAndBack
131 {
132 type empty;
133 faces
134 (
135 (0 3 4 1)
136 (2 5 6 3)
137 (3 6 7 4)
138 (5 8 9 6)
139 (6 9 10 7)
140 (11 14 15 12)
141 (13 16 17 14)
142 (14 17 18 15)
143 (16 19 20 17)
144 (17 20 21 18)
145 );
146 }
147);
148
149// ************************************************************************* //
如上图所示，在blockMeshDict字典文件中，第22行至第44行，为上图中编号从0只21的位置点信息。convertToMeters关键词用于单位缩放。比如其指定为1时，则所有的数字的单位均为米。如果其指定为10时，则所有的数字需要乘以10（单位同样为米）。
上图中很明显的计算域被分解成了5个block，因此在blockMesh中的block关键词下存在5个block的指定。blockMeshDict文件的内容将在此进行简要回顾，但相关的详细信息，请参阅以后的章节（待更新）。
blockMeshDict以block的顶点的坐标开始。然后，该文件定义了block（在这里是5个）。每个块都是一个六面体形状，由hex关键词给出。hex关键词后面列出了八个顶点标签。每个block的各个方向上的单元数量由三个整数组成的向量来指定。 例如，第一个块指定 了(18 30 1)，表明在 x方向划分出18个单元格，在 y方向上生成30个单元，在z方向（未使用的方向）上仅给出1个单元格。
这些block包含网格非均一分布。它还包含多级的非均一分布，这在后续会进行详细描述。
网格区分为入口、出口和壁面等边界，具体包含以下几个部分：upperWall和lowerWall代表壁面；inlet和outlet代表进出口。在z法线方向（即垂直于z轴的方向）上的边界，被归入一个名为frontAndBack的单独区域中。其为一个2D算例，因此是一个empty的边界条件类型。
具体来说，通过在blockMeshDict文件上运行blockMesh来生成网格。在案例目录中，只需在终端输入以下内容即可：
blockMesh
执行命令后，blockMesh的运行状态会在终端窗口中显示。如果blockMeshDict文件中有任何错误，blockMesh会输出相应的错误消息。
查看网格
在计算之前需要先看一下网格。用户只需在案例目录中执行paraFoam命令，就可以轻松地启动ParaView进行后处理。
所有linux程序都可以有两种方式运行：
前端运行：即终端会等待该程序执行完毕后才结束；
后端运行：在运行的同时，可以执行其他程序；
在这里，用户希望运行ParaView的同时，同时可以运行一些别的程序。那可以这样执行：
paraFoam &
然后在“Properties”窗口的“Mesh Parts”面板中选择要查看的几何图形。点击“Apply”。如果网格的边界数量比较少，您可以直接选中“Mesh Parts”面板上面的复选框快速选择所有边界。完成选择后可点击“Apply”。
ParaView顶部的第二行按钮可以展示网格的不同形式。其也可以通过下面的Display面板来进行设置。 用户可以这样操作：
选择Solid Color确定颜色；
点击 Edit按钮（位于“着色”部分内），选择一个合适的颜色，比如黑色；
在Representation菜单中选择Wireframe模式，这样网格看起来比较清楚；背景颜色也可以在 Properties窗口底部的 View（Render View）面板中进行设置；
另外，属性窗口中的很多高级设置，需要通过点击窗口顶部搜索框旁边的“高级属性”齿轮按钮（⚙️）才能显示出来。
边界和初始条件
网格分完后，可以对初始场进行配置。该案例从时间t = 0s开始，因此初始场数据被存储在cavity文件夹的0文件夹中。这个0文件夹包含几个文件，其中包括压力p和速度U。首先我们来看p：
16dimensions [0 2 -2 0 0 0 0];
17
18internalField uniform 0;
19
20boundaryField
21{
22 inlet
23 {
24 type zeroGradient;
25 }
26
27 outlet
28 {
29 type zeroGradient;
30 value uniform 0;
31 }
32
33 upperWall
34 {
35 type zeroGradient;
36 }
37
38 lowerWall
39 {
40 type zeroGradient;
41 }
42
43 frontAndBack
44 {
45 type empty;
46 }
47}
该文件主要有三个关键词：
dimensions 单位：这里是运动压力，即\(\mathrm{m}^2 \mathrm{s}^{-2}\)（除掉了密度的单位）；
internalField 内部场的值：其可以是均一值，也可以是非均一的，这种情况下必须为所有网格指定具体值；
boundaryField 边界场的值：包括边界条件和所有边界块的数据；
对于pitzDailySteady案例，初始场是均一的。
Note
由于是不可压缩流，这里的压力是运动压力，且其绝对值并不重要，因此为方便起见，将其设置为uniform 0。
边界包括upperWall, lowerWall, inlet和outlet。壁面和入口的压力p的边界均为zeroGradient，这意味着“压力的法向梯度为零”。出口使用了固定值fixedValue边界条件，压力值设为uniform 0。二维模拟的两侧平面的frontAndBack边界面被指定为empty。
对于速度U。内部场初始化设为零，但是速度必须通过3个分量来表示，即(0 0 0)。速度在壁面上应用无滑移条件noSlip。入口流速U在x方向为10m/s，通过fixedValue条件进行设置，具体数值为uniform (10 0 0)。出口使用zeroGradient。frontAndBack面设为empty。
物性
物性配置都存储在constant目录下的字典文件中。具体物性在physicalProperties文件中指定：
17viscosityModel constant;
18
19nu 1e-05;
其中的viscosityModel设置为常数。nu关键词来指定运动粘度，运动粘度用希腊符号\(\nu\)表示。其值设定为1e-05。
湍流模型
为了确定流动是否为湍流，需要估算雷诺数。雷诺数的定义为：
\(\bfU\)和\(L\)分别代表特征速度和特征长度，\(\nu\)表示运动粘度。本算例入口台阶高度\(L =25.4\)毫米，当速度为10米每秒时，雷诺数\(Re=25400\)。对于管道中的流动，当雷诺数\(Re > 2000\)时，流动状态通常会由层流转为湍流，所以该算例为湍流。
momentumTransport文件中可以设置牛顿流体、非牛顿流体、湍流以及粘弹性等。本算例包含了如下湍流模型配置：
17simulationType RAS;
18
19RAS
20{
21 // Tested with kEpsilon, realizableKE, kOmega, kOmega2006, kOmegaSST, v2f,
22 // ShihQuadraticKE, LienCubicKE.
23 model kEpsilon;
24
25 turbulence on;
26
27 printCoeffs on;
28
29 viscosityModel Newtonian;
30}
simulationType需指定层流或湍流。RAS表示雷诺平均湍流模型。RAS子字典中包含一个model关键词，该关键词若为kEpsilon，则表示其为\(k-\varepsilon\)模型。turbulence关键字如果是on，表示开启湍流模型。所有湍流模型都有默认值，这些默认值可以在momentumTransport文件中通过指定关键词来进行覆盖。当printCoeffs开关打开时，系数将会出现到终端。
Note
viscosityModel关键词表示流体为牛顿流体，牛顿流体是默认的流体类型，因此其可以省略。
\(k-\varepsilon\)模型需要求解湍流动能\(k\)和湍流耗散率\(\varepsilon\)的传输方程。这些场的初始和边界条件分别在0文件夹下的k和epsilon文件中进行设置。
湍流变量的边界场要尽可能的准确，这样计算才好收敛。湍流动能\(k\)可以通过以下公式计算：
其中\(I\)表示湍流强度。在本例中估算\(I=0.05\)，因此 \(k=1.5\times \left ( 10\times 0.05 \right ) ^{2}=0.375 \rm m^{2} \rm s^{-2}\)。因此在k文件中，0.375可以同时用作初始内部场和入口值。
湍流耗散率 \(\varepsilon\) 可计算为
湍流耗散率 \(\varepsilon\) 可以通过 \(C_{\mu}=0.09\) 和普朗特混合长度 \(l_{\rm m} \) 的估值来计算。在这个例子中，\(l_{m}\)为台阶高度的0.1倍，即2.54毫米, 从而得到\(\varepsilon =0.09^{0.75} \times 0.375^{1.5} /0.00254=14.855 \rm m^{2} \rm s^{-3} \)。在epsilon文件中，14.855同时用作初始内部场和入口值。
Note
\(l_{m}\)与\(I\)是用户完全需要自行预估的变量。
在模拟壁面的时候，可使用壁面函数。壁面函数通过0文件夹下的nut、k以及epsilon文件来指定。首先看nut：
16
17dimensions [0 2 -1 0 0 0 0];
18
19internalField uniform 0;
20
21boundaryField
22{
23 inlet
24 {
25 type calculated;
26 value uniform 0;
27 }
28 outlet
29 {
30 type calculated;
31 value uniform 0;
32 }
33 upperWall
34 {
35 type nutkWallFunction;
36 value uniform 0;
37 }
38 lowerWall
39 {
40 type nutkWallFunction;
41 value uniform 0;
42 }
43 frontAndBack
44 {
45 type empty;
46 }
47}
本案例使用标准壁面函数，分别在upperWall和lowerWall通过nutkWallFunction类型进行指定。其他可选的壁面函数模型还包括粗糙壁面函数，其可以通过nutRoughWallFunction关键字进行指定。同时，epsilon相应的壁面需指定为epsilonWallFunction，k需要指定为kqRWallFunction。
Note
kqRWallFunction是一种通用的边界条件，可应用于任何湍流动能类型的场，例如\(k\)、\(q\) 或雷诺应力\(\bfR\)。
求解控制
controlDict文件可以进行求解层面的相应设置。其位于system目录下：
16
17application foamRun;
18
19solver incompressibleFluid;
20
21startFrom startTime;
22
23startTime 0;
24
25stopAt endTime;
26
27endTime 2000;
28
29deltaT 1;
30
31writeControl timeStep;
32
33writeInterval 100;
34
35purgeWrite 0;
36
37writeFormat ascii;
38
39writePrecision 6;
40
41writeCompression off;
42
43timeFormat general;
44
45timePrecision 6;
46
47runTimeModifiable true;
48
49cacheTemporaryObjects
50(
51 kEpsilon:G
52);
53
54functions
55{
56 #includeFunc streamlinesLine
57 (
58 name=streamlines,
59 start=(-0.0205 0.001 0.00001),
60 end=(-0.0205 0.0251 0.00001),
61 nPoints=10,
62 fields=(p k U)
63 )
64
65 #includeFunc writeObjects(kEpsilon:G)
66}
关键词incompressibleFluid表示模拟不可压缩等温流体的稳态或瞬态流动，同时支持动网格；
本算例的模拟从0时刻开始，因此求解器从0文件夹中读取数据。这样需要将startFrom关键字设为startTime，并将startTime关键字设定为0；
因为其是一个稳态模拟。所有方程中的时间导数 \(∂/∂t\) 均为0，其可以在fvSchemes文件中，将ddtSchemes设置为steadyState来实现；在稳态下，关键字deltaT表示的时间步长仅用作时间增量（因为它不再用于离散化 \(∂/∂t\)）。它的值不会影响结果，因此对于稳态解决方案，它被设置为1，时间步仅仅表示代表求解步数；
endTime用于设置求解停止的时间。设为2000可以提供足够数量的求解步数，使稳态解能够收敛；
随着计算的进行，通常希望在特定的时间间隔内记录结果。writeControl提供了多种选项来设定记录结果的时间点；这里的timeStep选项指定了每隔n个时间步记录一次结果，n的值由writeInterval关键字指定。如果间隔设为100，那么结果将在第100步、第200步等时间点被记录；OpenFOAM在每次写入一组数据时，都会创建一个以当前时间命名的新目录，例如100。其包含所有的场的信息；