第二章:后向台阶流(校对中)

后向台阶流为一个二维的、等温的、不可压缩流动。其计算域由以下特点组成:

  • 左边存在一个进口;

  • 右边存在一个出口;

  • 上下均为边界;

流体从左侧进口以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面板来进行设置。 用户可以这样操作:

  1. 选择Solid Color确定颜色;

  2. 点击 Edit按钮(位于“着色”部分内),选择一个合适的颜色,比如黑色;

  3. 在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。

湍流模型

为了确定流动是否为湍流,需要估算雷诺数。雷诺数的定义为:

(1)\[\begin{equation} Re=\frac{\left | \bfU \right | L}{\nu} \end{equation}\]

\(\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\)可以通过以下公式计算:

(2)\[\begin{equation} k=\frac{2}{3} \left ( \left | \bfU \right | I\right ) ^{2} \end{equation}\]

其中\(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\) 可计算为

(3)\[\begin{equation} \varepsilon =C_{\mu }^{0.75} \frac{k^{1.5} }{l_{\rm m} } \end{equation}\]

湍流耗散率 \(\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。其包含所有的场的信息;

离散化和线性求解器

用户在system目录下的fvSchemes文件中指定有限体积离散化方案的选择。线性方程求解器、容差和其他算法控制的规格在系统目录下的fvSolution文件中进行。

这两个文件的详细信息分别在4.5和4.6节中描述。对于这个例子,需要注意以下几点:

  • ddtSchemes 默认为 fvSchemes 中的 stableState,以调用稳态计算;

  • 稳态解使用了基于SIMPLE的算法,其控制项在fvSolution文件中的SIMPLE子字典中进行了配置;

  • SIMPLE子字典包含一个consistent开关,它被设置为yes,应用了SIMPLE算法的“consistent”形式 (SIMPLEC);

  • SIMPLEC的收敛对fvSolution文件中的relaxationFactors非常敏感;0.9的值经过测试的,并不适合标准的SIMPLE算法;

  • SIMPLEC的敏感收敛性通常使其仅适用于几何形状简单的情况。

运行应用程序

在对单一区域进行模拟时,可采取foamRun应用。它从controlDict文件中的solver加载相关的求解器模块来执行计算。这次,我们将在终端运行foamRun,通过在案例目录中输入以下内容来执行。

foamRun

模拟进度将在终端窗口进行报告。它描述了连续的求解步骤,并给出了每个方程的初始和最终残差以及守恒误差。

fvSolution中的SIMPLE包括一个residualControl子字典,其中包含一组p、U和湍流场的容差。当所有场的初始残差低于其各自的容差时,incompressibleFluid求解器将终止计算。在此示例中,求解器在287次迭代后终止,并显示以下声明。

SIMPLE solution converged in 287 iterations

模拟结果将被写入pitzDailySteady目录内的时间目录中。用户可以使用“ls”命令列出目录内容,以查看包含结果的时间目录(100、200和287)。

ParaView中的时间选择

一旦结果写入时间目录,就可以使用ParaView查看它们。第一步是激活时间选择器,包括按钮顶部一行右侧的“Time”框,如图7.4所示。

如果在没有任何时间目录(即只有一个0目录)前打开ParaView,则必须重新激活时间选择器,通过以下方式解决时间目录:

  • 在ParaView 中选择属性窗口的顶部(如有必要,向上滚动面板);

  • 切换面板顶部的“缓存网格”(Cache Mesh)按钮(位于“刷新时间”(Refresh Times))按钮下方;

  • 单击“Apply 按钮。

随后,时间选择器会进行更新,Time 变成一个下拉菜单,其中包含案例中的时间目录(0、100、200和287)。为了查看287时的解,用户可以使用 VCR Controls or Current Time Controls 将当前时间更改为287。

对于教程其余部分:

在ParaView中,如果窗口面板(例如“Properties”)不包含预期的关键词,请确保在 Pipeline Browser 中突出显示相关模块(以蓝色显示)。例如,仅当选择顶部模块(此处为pitzDailySteady.OpenFOAM)时,“Refresh Times”才会出现。

表面着色

按照图2.4所示进行压力展示,从ParaView顶部的第二行按钮中进行选择,或者滚动到“属性”中的“显示”面板,然后进行以下选择:

    1. 从“Representation”菜单中选择“Surface”;

    1. 在“Coloring”中选择

    1. 如有必要,单击“Rescale”按钮将色标设置为数据范围。

压力场应如上所示,在出口通道收缩时压力增加到最大值。使用点图标(),压力场在每个单元之间进行插值,以呈现出连续的展示。相反,如果用户从“Coloring”菜单中选择单元格图标(),则每个单元格将被赋予一个压力值,由没有过度的单一颜色表示。

颜色图例可以通过点击ParaView顶部第二行按钮左侧的“Toggle Color Legend Visibility”按钮来禁用。这些按钮是“Active Variable Controls”工具栏的一部分,如图7.4所示。“Active Variable Controls”工具栏中左侧的第二个按钮“Edit Color Map”会打开“Color Map Editor”窗口,如图2.5所示,用户可以在其中设置颜色比例尺和颜色条的一系列属性。

特别是,ParaView默认使用蓝色-白色-红色的颜色比例尺,而不是更常见的蓝色-绿色-红色(彩虹)。因此,当第一次执行ParaView时,可能需要更改颜色比例尺。可以在“Color Scale Editor”中选择“Choose Preset”按钮(带有心形图标)来完成。CFD的常规颜色比例尺是“Blue to Red Rainbow”,仅当用户在搜索栏中输入名称或选中该栏右侧的齿轮图标时才会列出。

选择“Blue to Red Rainbow”并单击“Apply”和“Close”后,可以单击面板底部的“Save as Default”按钮(文件保存符号),以便ParaView始终采用这种类型的颜色条。

用户还可以通过单击搜索栏最右侧的“Edit Color Legend Properties”来编辑颜色图例属性,例如文本大小、字体选择和比例的编号格式,如图 2.5 所示。

切面

如果通过按住图像窗口中的鼠标左键并移动光标来旋转图像,可以看到已经通过压力对整个几何表面进行着色。为了生成二维等高线图,用户应首先创建一个切割平面,或称为“切片(slice)”。在 Pipeline Browser 中突出显示 pitzDailySteady.OpenFOAM 模块后,用户应从ParaView顶部菜单的“Filters”菜单中选择“Slice”(在某些系统上可以通过屏幕顶部访问)。可以在“Common”子菜单中找到“Slice”,或者在ParaView顶部的第三行按钮中找到“Common and Data Analysis”按钮(见图7.4)。

选择“Slice”会在Pipeline Browser中创建一个新项。在“Properties”窗口中,切割平面的原点应位于z轴上的某一点,例如(0,0,0),其法线应设置为(0,0,1)(单击Z法线按钮)。

单击“Apply”后,切片将出现在“RenderView”窗口中,而原始的pitzDailySteady.OpenFOAM模块将消失。每个模块的可见性可以通过Pipeline Browser中每个模块左侧的眼睛按钮来启用和禁用。

对于教程其余部分:

在 ParaView 中,如果 RenderView 窗口中没有显示关键词,请通过打开Pipeline Browser中的eye按钮来确保相关模块可见。

矢量图

在绘制流速矢量之前,请通过在Pipeline Browser中突出显示“Slice”模块并单击其左侧的眼睛按钮来关闭其显示。目标是在每个单元的中心生成一个速度向量图。因此,首先需要从网格几何体中过滤出单元中心,如7.1.7节所述。在Pipeline Browser中高亮显示pitzDailySteady.OpenFOAM模块后,用户应从“Filters->Alphabetical”菜单中选择“单元中心”,然后单击“应用”。

在Pipeline Browser中突出显示Centers后,用户应该从“Filters->Common”菜单(或第三行按钮)中选择“Glyph”。 属性窗口面板应如图 2.6 所示。

显示速度矢量时,glyphs需要四种主要设置:

  • glyph类型,设置为Arrow;

  • 箭头方向,由 Orientation Array 设置;

  • 箭头长度,由Scale Array and Scale Factor设置;

  • 箭头的数量,由Glyph Mode设置。

单击“Apply”后,glyphs显示为单色,例如 白色的。 用户可以根据速度大小给glyphs着色,通过在第二行按钮左侧的下拉菜单中设置 U 来控制。

Legend也被显示。 用户可以通过单击“Edit Color Map”(第二行左侧第二个按钮)来配置图例。 这将打开Color Map Editor 窗口。可以通过单击搜索框最右侧的按钮来配置图例。 此按钮打开“Edit Color Legend Properties”窗口。 应选中搜索框右侧的“Advanced Properties”齿轮按钮。

标题和标签可以完全配置。对于图2.7,图例标题设置为速度U[m/s]。通过取消“Automatic Label Format”并在“Label Format”文本框中输入%-#6.1f,可以将标签指定为1位固定小数点。同时,“Add Range Labels”也未选中。

图2.7显示了样本输出,放大了台阶下游部分再循环区域。在下壁面,glyphs指向与壁面附近流动相反的方向。上述情况是由于这些glyphs被绘制在壁面边界块的中心,并且为无滑移条件,速度大小为0而引起的。向量没有任何方向,ParaView将箭头定向为默认的x方向。快速删除这些向量的方法是:

  • 返回到Pipeline Browser顶部的pitzDailySteady.OpenFOAM模块;

  • 在 Mesh Parts 面板中,取消选中 lowerWall patch;

  • 点击 Apply.

常用滤镜

ParaView包含200多个过滤器,这些过滤器可以通过“Filters->Alphabetical”菜单列出。其中只有一小部分,例如10-15个过滤器与CFD相关,因此建议将这些添加到“Filters->Favourites”菜单中。为此,请从“Filters->Favourites”菜单中选择“Manage Favourites”。搜索以下重要的过滤器,然后单击“Add>>”将它们添加到“Favourites”菜单中:

  • Extract Block,用于选择域的组件,例如边界块和内部单元;

  • Slice,用于在几何体中插入一个平面;

  • Cell Centers and Glyph, 主要用于绘制速度矢量图;

  • Stream Tracer and Tube, 用于绘制流线;

  • Contour,用于在表面上绘制等高线和等值面;

  • Feature Edges,用于捕捉表面上的特征以获得更好的图像清晰度。

等高线

在准备绘制流速的等高线前,先在“Pipeline Browser”中选择“Glyph”模块,并点击其左侧的眼睛图标以隐藏该模块的显示。随后,用户需选中“Slice”模块,并按流速U来着色。然后在切片上以U的1、2、…、10为间隔来绘制等高线。

要绘制等高线,我们可以对切片应用“Contour”过滤器。此过滤器能够依据标量的恒定值,在二维空间中绘制线条或在三维空间中绘制曲面。由于U是一个向量,因此我们不能直接根据U绘制等高线。为此,我们需要首先生成一个表示U大小的标量场。

使用函数对象进行后处理来生成mag(U)场,具体方法在7.3节中有详细描述。用户需要运行foamPostProcess工具,并通过以下方式使用“-func”选项来调用mag函数对象:

    foamPostProcess -func "mag(U)"

该程序会逐个处理所有时间目录, 并读取每个时间目录内得U值,计算其模值mag(U)。然后将结果作为字段文件存回相应的时间目录。接下来,需要将mag(U)字段加载到ParaView中,操作步骤如下:先在“Pipeline Browser”中找到pitzDailySteady.OpenFOAM,点击“Refresh Times”;然后滚动至“Fields”面板,选中mag(U)后点击“Apply”。

操作者需要在“Pipeline Browser”中再次选中“Slice”模块,并应用“Contour”过滤器。接着,在“Properties”面板中,从“Contour By”下拉菜单里选择mag(U)。在“Isosurfaces”部分,先点击“-”按钮清除默认值,然后按照图2.8所示添加10个值。最后,可以选择以Wireframe形式展示等高线,并使用纯黑色进行着色。

流线图

在绘制流线之前,请确保关闭了“Contour”模块的显示。对于本二维案例需展示的流线,用户先在“Pipeline Browser”中选择“Slice”模块,随后应用“Stream Tracer”过滤器。

将会打开一个新的“StreamTracer”模块,并通过其“Properties”窗口进行配置。追踪器通过从 seed points 开始,沿流动方向追踪线条来创建。若选择“Integration Direction BOTH”,则线条会在种子点的上游和下游都进行追踪。用户应滚动“Properties”窗口来配置“Seeds”。默认的Seed Type为“Line”,即沿着在指定点之间绘制的线条播种点。在“Line Parameters”中,用户应将两点设置为从(0, −0.025, 0)到(0.2, 0.025, 0)。分辨率则指定了沿线条分布的种子点数量,建议减少到25个。点击“Apply”后,将如图2.9所示生成追踪器。用户可以通过调整线条点和分辨率来获得不同的流线追踪输出结果。

入口边界条件

用户需要检查流动域的入口边界处的流动状况。在\(0/U\)中,将速度条件设定为固定值(10 0 0)。这个速度值会应用于入口边界的所有面。入口与无滑移条件的墙壁边界相邻,这会导致相邻边界面之间的流速U发生突变。

用户应用鼠标右键放大到几何图形的入口区域。随后,在“Pipeline Browser”中激活“Slice”模块,并根据p值的cell values来着色(请进行选择 )。为了凸显入口附近单元格的压力变化,用户应根据图2.10所示,设定一个自定义范围,即 \(-5 < p < -1\)。要实现自定义范围,请点击“Rescale to Custom Data Range”按钮(, 该按钮位于第二行按钮中,从左数第五个)。之后会弹出一个面板,在此输入范围的最小值和最大值,再点击“Rescale”按钮。若要查看速度分布图,请返回“Pipeline Browser”中的“Glyph”模块,并将其设置为可见。根据流速调整箭头大小,可以描绘出相应的轮廓。在“Glyph”的Properties中,滚动到“Scale”面板,将“Scale Array”设置为 U,选择“Scale By Magnitude”并设定一个Scale Factor为8e-05。此外,为向量分配一种白色Solid Color。

图示为入口区域与下侧壁边界的相邻部分。在图像左侧,向量显示出均匀的分布。由于墙壁上的剪切力作用,流动开始从近壁处减速。这种减速会引起压力增加,从而产生一个驱动力梯度,使流动略微偏离墙壁,以维持质量守恒。

为了减小压力的增长,通过调整边界条件,使得入口速度不再保持均一。其中,“flowRateInletVelocity”边界条件是用于指定入口流量的通用U边界条件。可以使用“foamInfo”脚本来查看这一边界条件的详细文档。这个脚本是一个通用工具,能够提供OpenFOAM中应用程序、模型和工具的文档。对于“flowRateInletVelocity”边界条件,可以按以下方式运行该脚本(请注意,有时为了简化,该名称可以缩写,比如这里使用“flowRateInlet”而非“flowRateInletVelocity”)。

    foamInfo flowRateInlet

foamInfo脚本能够定位并输出相关代码的头文件,并从中提取“Description”和“Usage”。此外,它还能识别出相关的模型(在这种情况下指的是其他边界条件),并列出使用该模型的示例案例。虽然foamInfo脚本并非尽善尽美,但在大多数情况下(十次中有九次),它都能迅速提供有用的信息。

根据文档内容,flowRateInletVelocity条件可通过massFlowRate, volumetricFlowRate 或 meanVelocity 来指定流量。用户应在文本编辑器中打开“\(0/U\)”文件,并找到入口区域的边界字段关键词。随后,可以通过修正inlet子字典来应用meanVelocity条件,修改方式如下:

    inlet
    {
          type               flowRateInletVelocity;      // modify
          meanVelocity       10;                         // insert
          value uniform      (10 0 0);                   // leave
    }

此条件会对边界上的速度进行评估,并为边界区域上的所有面设定相应的值。尽管对于OpenFOAM而言,“值”这一关键词是多余的,但ParaView需要它来显示初始化区域字段的初始值。因此,保留了“值”关键词,以便于ParaView的使用。

在保存了“\(0/U\)”文件之后,用户可以重新运行,首先确认“flowRateInletVelocity”是否成功模拟了原始的“fixedValue”条件。在“controlDict”文件中,案例的起始时间设定为0,并且它将覆盖时间文件夹中的先前结果。用户只需重新运行“foamRun”求解器,即可测试新的条件。

    foamRun

求解器的运行过程与之前相同,最终在287次迭代后结束。用户可以返回ParaView,并在“Pipeline Browser”的“pitzDailySteady.OpenFOAM”模块中点击“Refresh Times”按钮。结果未发生变化,这表明“flowRateInletVelocity”在此阶段设置了一个统一的值“(10 0 0)”。

用户通过包含速度的profile文件来修改\(0/U\)”文件中的“flowRateInletVelocity” 条件。 文档中强调了两个自定义的分布函数:turbulentBL和laminarBL,它们分别为充分发展的湍流和层流边界层提供power-law和quadtratic。在此案例中,请按以下方式将turbulentBL分布添加到边界条件中:

  inlet
   {
         type               flowRateInletVelocity;   
         meanVelocity       10;  
         profile            turbulentBL; // add        
         value uniform      (10 0 0);    
   }

在重新运行模拟之前,建议先删除先前的解决方案时间文件夹,以避免混淆。因为求解器可能会在不同的时间点终止,如果不删除先前的结果,旧案例的最终时间文件夹中的结果将不会被新结果覆盖,这可能会导致混淆。因此,现在应该删除先前的结果。

foamListTimes实用工具为用户提供了一种快速简便的方式来删除解决方案的时间文件夹,同时保留“0”文件夹。首先,请在终端中按以下方式运行foamListTimes:

    foamListTimes

执行上述命令后,会返回解决方案时间文件夹的列表,其中包括100、200和287,但不包括“0”文件夹。要删除这些列出的文件夹,可以在foamListTimes命令中添加“-rm”删除选项,即运行以下命令:

    foamListTimes -rm

现在,删除了包含之前结果的时间文件夹,接下来像之前一样运行foamRun求解器。这次,foamRun在277次迭代后终止。然后,用户应返回ParaView,并在“Pipeline Browser”的“pitzDailySteady.OpenFOAM”模块中点击“Refresh Times”按钮。

输出时间的变化可能会让ParaView造成混乱,导致它在时间选择器中显示带有“?”符号的时间查询。为了解决这个问题,用户可以选择恢复到时间0,或者点击“Cache Mesh”和“Apply”来重建选择器关键词。然后,用户可以选择序列中的最后一个时间点(277)。此时,速度剖面和压力会更新,如图2.11所示。对于压力,自定义范围调整为 \(-7<p<-3\)。入口处的速度不再保持均匀,而是根据turbulentBL函数形成了一个特定的剖面。与下壁面相邻面的速度与之前相比显著降低,因此近壁单元内由剪切力引起的减速效应也相应减弱。结果显示,左上角与其周围单元之间的压力差没有以前那么大,现在大约是 \(+2 \rm m^{2} s^{-2} \)(从-7变为-5),而之前则是 \(+4 \rm m^{2} s^{-2} \)(从-5变为-1)。

湍流模型

逆向台阶案例的设计初衷是让用户能够快速尝试多种不同的湍流模型。在constant文件夹中的momentumTransport文件内,有一条注释详细列出了在此案例上测试过的各种湍流模型。部分模型除了解决k和ε的方程外,还解决了如ω(omega)等其他场的方程。为了尽可能减少在更换模型时的工作量,我们在0文件夹中已经预先包含了这些其他场的相关文件。同样,在system文件夹中的fvSchemes和fvSolution文件里也分别加入了这些场的计算方案和求解器的相关关键词。

为了更改湍流模型,用户应在文本编辑器中打开momentumTransport文件。接着,用户可以通过以下设置,将模型切换为可实现 \(k-\varepsilon\) 模型。

    model realizableKE;

现在,若要重新运行使用此模型的模拟,请先删除解决方案时间文件夹,然后执行以下foamRun命令:

    foamListTimes -rm
    foamRun

求解器在进行了255次迭代后终止。用户返回ParaView,并点击“Refresh Times”按钮来查看第255秒的结果。结果如图2.12所示,该图使用了之前设置好的速度场切片和流线过滤器。

可实现k-ε模型的扩散性低于标准k-ε模型。它在台阶底部捕获了一个次级涡旋,其高度约为台阶高度的一半。再循环区域也更长,在下壁开始向出口逐渐变细的位置发生再附着。

相比标准 \(k-\varepsilon\) 模型,可realizable \(k-\varepsilon\) 模型的扩散性更低。它能够捕捉到台阶底部的一个次级涡旋,其高度大约是台阶高度的一半。此外,再循环区域也更长,再附着现象发生在下壁面,并且开始逐渐收窄直至出口的位置。

用户也可以尝试其它湍流模型。在工业CFD领域,k-ω SST模型备受欢迎。用户可以通过在momentumTransport文件中进行以下设置来选择该模型。

    model kOmegaSST;

这个模型需要对特定的湍流动能耗散率ω进行初始化。这个参数的计算方法与等式2.3中的ε类似,具体计算方式如下,其中会用到\(l_{m} \)

(4)\[\begin{equation} \omega = C_{\mu }^{-0.25} \frac{k^{0.5} }{l_{m} } . \end{equation}\]

使用之前的估计 \(l_{m} =2.54\)mm,则 \(\omega =0.09^{-0.25} \times 0.375^{0.5} /0.00254=440.2\rm s^{--1} \)。 在 \(0/omega \)文件中,440.2 用于初始internalFields和入口value。

当前,用户可以删除解决方案的时间文件夹,并像之前一样执行foamRun命令,并使用此模型重新模拟。值得注意的是,这次求解器不会因为fvSolution文件的residualControl子字典中指定的残差控制容差范围内收敛而提前终止。相反,它会在达到2000次迭代时终止,这是controlDict文件中设定的 endTime。

采用k-ω SST模型的结果显示,在台阶的底部存在一个更大的次级涡旋。这个涡旋并没有达到稳定状态,而是在连续的求解步骤中出现小幅度的振荡。这些振荡可以通过观察不同的求解步骤(如1000、1100、1200等)中的速度矢量来发现。图2.13展示了次级涡旋的范围,并指出了在涡旋重新附着点周围矢量振荡的位置。

可以通过修改动量平流的数值方法,以此来稳定次级涡旋。用户应从“system”文件夹中打开“fvSchemes”文件。在“divSchemes”子字典中,通过形如“div(phi,…)”的关键词来指定平流项的离散化方式。动量平流采用的是线性迎风格式,具体如下所示。

    div(phi,U)    bounded Gauss linearUpwind grad(U);

linearUpwind 格式利用单元梯度,通过外推法将场从单元中心插值到单元面上。其中的grad(U)关键词用于确定梯度计算的方式,具体采用“gradSchemes”子字典中指定的方案。在本例中,它仅包含一个default 格式,用于计算所有方程式中的所有梯度项。

    gradSchemes
    {
          default   Gauss linear;   
          grad(U)   cellLimited Gauss linear 1; 
    }

保存fvSchemes文件之后,用户可以删除之前的时间文件夹,并像之前一样执行foamRun命令来重新运行模拟。在应用了cellLimited格式后,解会正常收敛,求解器在迭代285次后终止。此时,次级涡旋会稳定下来,如图2.14所示。