第三章:溃坝

本案例中,采用不可压缩的VoF模型求解器,求解简化的二维溃坝问题。本案例为瞬态流动,两相流体由一个明确的界面或自由表面分隔开。在incompressibleVoF中,两相算法是以流体体积(VoF)为基础的。在这种算法中,使用相传输方程来确定每个计算单元中各相的相对体积分数,或称为相分数α。根据这个相对体积分数,通过加权平均计算得出材料的物理属性。VoF方法的特性意味着相之间的界面并非显式计算得出,而是作为相分数场的一个特性显现。由于相分数的取值范围在0到1之间,因此不能精确定义相界面,而是占据应存在尖锐界面的区域体积。

仿真区域由位于水箱左侧隔膜后的一列静止水柱构成。在时间t=0s时刻,移除隔板,水柱崩塌。在崩塌过程中,水冲击水箱底部的障碍物,并产生一个复杂的流场结构,流场中包括几个空气泡。几何模型和初始设置如图2.15所示:

网格生成

打开案例运行目录,并从$FOAM_TUTORIALS/incompressibleVoF/damBreakLaminar目录中复制damBreak案例,即

run
cp -r $FOAM_TUTORIALS/incompressibleVoF/damBreakLaminar/damBreak .

进入damBreak文件夹,然后运行blockMesh以生成网格。damBreak网格由五个block构成,blockMeshDict的相关条目如下所示。


16convertToMeters 0.146;
17
18vertices
19(
20    (0 0 0)
21    (2 0 0)
22    (2.16438 0 0)
23    (4 0 0)
24    (0 0.32876 0)
25    (2 0.32876 0)
26    (2.16438 0.32876 0)
27    (4 0.32876 0)
28    (0 4 0)
29    (2 4 0)
30    (2.16438 4 0)
31    (4 4 0)
32    (0 0 0.1)
33    (2 0 0.1)
34    (2.16438 0 0.1)
35    (4 0 0.1)
36    (0 0.32876 0.1)
37    (2 0.32876 0.1)
38    (2.16438 0.32876 0.1)
39    (4 0.32876 0.1)
40    (0 4 0.1)
41    (2 4 0.1)
42    (2.16438 4 0.1)
43    (4 4 0.1)
44);
45
46blocks
47(
48      hex (0 1 5 4 12 13 17 16) (23 8 1) simpleGrading (1 1 1)
49      hex (2 3 7 6 14 15 19 18) (19 8 1) simpleGrading (1 1 1)
50      hex (4 5 9 8 16 17 21 20) (23 42 1) simpleGrading (1 1 1)
51      hex (5 6 10 9 17 18 22 21) (4 42 1) simpleGrading (1 1 1)
52      hex (6 7 11 10 18 19 23 22) (19 42 1) simpleGrading (1 1 1)
53);
54
55defaultPatch
56{
57       type empty;
58}
59
60boundary
61(
62        leftWall
63       {
64             type wall;  
65             faces
66             (
67                   (0 12 16 4)
68                   (4 16 20 8)
69             );
70       }
71        rightWall
72       {
73             type wall;  
74             faces
75             (
76                   (7 19 15 3)
77                   (11 23 19 7)
78             );
79       }
80        lowerWall
81       {    
82             type wall;     
83             faces 
84             (
85                   (0 1 13 12)
86                   (1 5 17 13)
87                   (5 6 18 17)
88                   (2 14 18 6)
89                   (2 3 15 14)
90             );
91       }
92        atmosphere
93       {    
94             type wall;     
95             faces 
96             (
97                   (8 20 21 9)
98                   (9 21 22 10)
99                   (10 22 23 11)
100            );
101      }  
102);
103
104
105// ************************************************************************* //

网格信息保存在constant目录下的polyMesh文件夹内的多个文件中。可以输入如下命令,列出该文件夹下的子文件,并查看:

ls constant/polyMesh

polyMesh文件夹包含如下文件:boundary(对边界域集合的描述)、faces(单元格面的列表)、neighbour和owner(包含与特定面相连的单元格索引)、points(单元格顶点的列表)。

边界条件

boundary文件的内容是可读的,可以使用文件编辑器打开,或使用cat工具在终端窗口中输出其内容。cat工具命令如下:

cat constant/polyMesh/boundary

该文件包含了五个边界,分别是:leftWall、rightWall、lowerWall、atmosphere以及defaultFaces。值得注意,关注这些边界的类型是必要的。

首先,atmosphere是一个标准的patch,即它不具备任何特殊属性,仅仅是一个用来指定边界条件的单元。 其次,defaultFaces是由在blockMeshDict文件的boundary子字典中被省略的区所构成的,这些被省略的区块共同组成了一个patch,其相关属性在blockMeshDict文件中的defaultPatch子字典里进行了详细规定。本案例中,由于该patch的法线方向是在这个二维案例中无法求解的方向,因此默认类型设置为空。

leftWall、rightWall和lowerWall边界均设定为wall类型。与通用的patch相似,wall类型不包含任何关于网格的几何或拓扑信息。它与标准的patch的区别仅在于,wall类型将区块标识为壁面。这对于特定模型来说是必要的,例如湍流壁函数,还有一些在计算中需要考虑距最近壁面距离的湍流模型。

defaultFaces边界代表二维域的前后平面,与atmosphere一样,是一个空类型。 特别是在VOF模拟中,表面张力模型可以模拟界面与壁面接触点之间的壁面粘附效应。壁面粘附模型可通过在alpha(α)场上应用特定的边界条件来实现,例如使用constantAlphaContactAngle边界条件,此时需要设定一个静态接触角theta0。

本例忽略壁面与界面之间的表面张力效应。要实现这一点,一种方法是设置静态接触角\(\theta _{0} =90^{\circ} \),另一种方法是在壁面上对alpha应用zeroGradient条件。

顶部边界与大气相通,因此需根据内部流动情况允许同时存在流出和流入。为了实现这一点并保持流动稳定性,需要将压力和速度的边界条件相结合。这些边界条件包括:

  • prghTotalPressure 这一边界条件应用于压力场,但要减去由公式6.3给出的静水压分量 \(p_{\rho gh} \)

  • pressureInletOutletVelocity这一边界条件应用于速度U,该条件会对速度U的所有分量设置zeroGradient(零梯度),但在流入区域,则对切向分量应用fixedValue(固定值)条件;

  • inletOutlet 这一边界条件应用于其他场,在向外流动时表现为zeroGradient条件,在向内流动时表现为fixedValue条件。

在所有壁面边界上,fixedFluxPressure边界条件应用于压力场,并通过调整压力梯度,使边界通量与求解器(包括重力和表面张力等体力)的速度边界条件相匹配。

流体相的具体信息在constant目录下的phaseProperties文件中指定。phaseProperties文件内容如下所示:


16 
17phases                       (water air);
18
19sigma                        0.07;
20
21
22// ************************************************************************* //

第17行列出了液相和气相。除了列表中的最后一个相(本例中是气相)之外,系统会为列表中的其它相求解相分数方程。由于这里只有两个相,所以只会为第一个相分数方程求解,即求解液相分数 \(\alpha _{water}\),具体数值在0目录下的alpha.water文件中指定。

PhaseProperties文件包含两相之间的表面张力,由关键字sigma指定,单位为 \(\rm N·m^{-1} \)

设定初始场

与上文的案例不同,如下案例将为相分数 \(\alpha _{water}\) 指定一个非均匀初始条件,其中:

(1)\[\begin{equation} \alpha =\begin{cases} & 1= \rm for\ the\ water\ phase\\ & 0= \rm for\ the\ air\ phase \end{cases} \end{equation}\]

这可以通过运行setFields来实现。运行该工具需要一个位于system目录下的setFieldsDict字典,针对本案例的字典条目如下所示。


16
17defaultFieldValues
18(
19       volScalarFieldValue alpha.water 0
20);
21 
22regions
23(
24    boxToCell
25    {
26           box (0 0 -1) (0.1461 0.292 1);
27           fieldValues
28           (
29                volScalarFieldValue alpha.water 1
30           ); 
31    }
32);
33
34
35// ************************************************************************* //

defaultFieldValues用于设定场的默认值,即在regions列表中没有特别指定时,场采用默认值。regions列表中包含一些带有fieldValues的子字典,这些子字典会在特定区域内覆盖默认值。区域是通过定义一组点、面或单元的函数来确定的。在本算例中,boxToCell函数创建了一个由最小和最大边界定义的区域,这个区域定义了液相域的单元集。在这个区域内,液相分数\(\alpha _{water}\)设定为1。

setFields工具从文件中读取各种场信息,重新计算这些场之后,将场信息写回文件。在本案例中,alpha.water场最初以原始形式存储,并命名为alpha.water.orig。如果实际的目标文件不存在,程序会读取带有.orig扩展名的场文件。因此,setFields会读取alpha.water.orig文件,但将输出结果写入alpha.water文件(如果启用压缩功能,则写入alpha.water.gz文件)。通过这种方式,原始文件不会被覆盖,可以重复使用。

执行如下命令:

setFields

使用 paraFoam,检查初始 alpha.water 场是否对应于所需的分布,如图2.16所示。

属性

气相和液相的物理特性分别在constant目录下的physicalProperties.air和physicalProperties.water文件中设置。这2个文件描述了流体在静止情况下的物理性质。并且,通过viscosityModel关键字指定粘度模型,该关键字设置为constant,表示该值是恒定的。然后,通过nu关键字以 \(\rm m^{2}·s^{-1} \) 为单位指定粘度。流体的密度也以 \(kg·m^{-3}\) 为单位通过rho关键字指定。下面以physicalProperties.air文件为例:


16
17viscosityModel     constant;
18
19nu                 1.48e-05;
20
21rho                1; 
22
23
24// ************************************************************************* //

如果粘度随剪切速率而变化,如非牛顿或粘弹性流体,则通过动量属性文件指定粘度模型,见后文第 8.3 节。

重力

整个区域内重力加速度是均匀的,并由constant目录中的g文件夹指定。与常规场文件(如U和p)不同,g是一个uniformDimensionedVectorField,因此它只包含一组维度和一个值,对于本算例来说,它代表(0, -9.81, 0) \(\rm m·s^{-2} \)


16
17dimensions           [0 1 -2 0 0 0 0];
18value                (0 -9.81 0);
19
20
21// ************************************************************************* //

湍流模型

与cavity案例相同,可以通过momentumTransport字典中的simulationType关键字来设置湍流模型。在本案例中,假定流动过程没有涡旋,因此设置为层流:


16
17simulationType     laminar;
18
19
20// ************************************************************************* //

时间步长控制

溃坝案例是瞬态仿真,因此需要注意时间步长的设定。库朗数Co(Courant number)与时间步长相关,是针对每个单元定义的无量纲参数,定义为:

(2)\[\begin{equation} Co=\frac{\delta t\left | \bfU \right | }{\delta x} \end{equation}\]

其中, \(\delta t\)是时间步长,\(\left | \bfU \right |\) 是通过该单元的速度大小,\(\delta x\) 是速度方向上的单元大小。对于显式求解,为了保持稳定性,需要满足最大Co小于1;更严格的稳定性限制取决于流体流动模型的选择。通常,隐式求解没有固定的最大库朗数稳定性限制,但当Co超过1时,时间精度将变得更加重要。

时间步长对于界面捕捉尤为重要。incompressible Volume of Fluid求解器模块使用了由Henry Weller创建的Multidimensional Universal Limiter for Explicit Solution (MULES),以保持相分数的有界性。根据所选的MULES算法,需要对库朗数进行限制。使用最初的显式MULES算法时,为了保持稳定,通常需要将库朗数的上限为 \(Co\approx 0.25\)。并且,还有一种半隐式的MULES算法,可通过在fvSolution文件中使用MULESCorr开关来指定。对于半隐式MULES算法,Co的稳定性没有上限值,而是根据时间精度的要求来决定的。

一般来说,很难指定一个固定的时间步长来满足Co标准。因为在仿真过程中,\(\left | \bfU \right |\)会随着单元的变化而变化。对于自动调整时间步长,需要将controlDict中的adjustTimeStep切换为on,并指定相场和其他场的最大库朗数maxAlphaCo和maxCo。在本例中,maxAlphaCo和maxCo都设置为1.0。时间步长的上限maxDeltaT,可以设置为本案例中不会超过的值,例如1.0。

当使用自动时间步长,时间步长本身不会被四舍五入。因此,OpenFOAM可按需要以固定的时间步长间隔保存结果。但是,通过调整自动时间步长,可以使用controlDict字典中writeControl的可调整adjustableRunTime选项在固定时间保持结果。使用该选项,自动时间步长程序会进一步调整其时间步长,以便在writeInterval指定的确切时间“命中”,本例中设置为0.05。,controlDict文件位于system文件夹中,其字典条目如下所示。

16
17application              foamRun; 
18
19solver                   incompressibleVoF;
20
21startFrom                startTime;
22
23startTime                0;
24
25stopAt                   endTime;
26 
27endTime                  1;
28
29deltaT                   0.001;
30
31writeControl             adjustableRunTime;
32
33writeInterval            0.05;
34
35purgeWrite               0;
36
37writeFormat              binary;
38
39writePrecision           6;
40
41writeCompression         off;
42
43timeFormat               general;
44
45timePrecision            6;
46
47runTimeModifiable        yes;
48
49adjustTimeStep           yes;
50
51maxCo                    1;
52maxAlphaCo               1;
53
54maxDeltaT                1;
55
56 
57// ************************************************************************* //

离散化

incompressibleVoF求解器所使用的MULES算法独立于底层数值方案、网格结构等来保持相分数的有界性。因此,对流方案的选择并不局限于那些强稳定或有界的方案,如upwind。

在 fvSchemes 字典的 divSchemes 子字典中,设置对流相离散化方案。本案例中,动量方程中的对流项 \(\nabla \cdot \left ( \rho \rm UU \right ) \)(用 div(rhoPhi,U) 表示)使用Gauss linearUpwind grad(U),并获得较好的离散化效果。

\(\nabla \cdot \left ( \rm U\alpha \right ) \) 项由div(phi,alpha)表示,使用了一种定制的interfaceCompression方法。其中,vanLeer系数是一个控制界面压缩的因子,具体为:0表示不压缩;1表示保守压缩;大于1的任何值则表示增强界面的压缩。通常使用的值为1.0,如本案例所示。

其它离散项使用常用的方案,因此fvSchemes字典条目如下所示。

16
17ddtSchemes
18{
19       default       Euler;
20}
21
22gradSchemes
23{
24       default       Gauss linear;
25}
26 
27divSchemes
28{
29       div(rhoPhi,U)  Gauss linearUpwind grad(U);
30       div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
31       div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
32}
33
34laplacianSchemes
35{
36       default       Gauss linear corrected;
37}
38
39interpolationSchemes
40{
41       default       linear;
42}
43
44snGradSchemes
45{
46       default       corrected;
47}
48
49
50// ************************************************************************* //

矩阵求解器

在system/fvSolution文件中,solvers内用于alpha.water的子字典包含特定于MULES算法的元素,如下所示。

    "alpha.water.*"
    {
        nAlphaCorr      2;
        nAlphaSubCycles 1;

        MULESCorr       yes;
        nLimiterIter    5;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
    }

MULES通过计算两个限制器来将相分数保持在0和1之间。限制器的是迭代计算的,迭代次数由nLimiterIter指定。根据经验,可以将nLimiterIter设置为3 + ,以保持Co数的有界性。

通过MULESCorr开关可以激活MULES的半隐式算法。在应用MULES作为高阶校正之前,首先需要计算出一个隐式格式的迎风解。然后,需要通过solver, smoother, tolerance 和 relTol等参数,为隐式迎风解配置线性求解器。

nAlphaCorr 关键字用于控制一个解步骤中相分数方程的迭代次数。迭代用于克服对流中的非线性,在这种情况下,非线性是由interfaceCompression方案引起的。

运行代码

代码的运行已在前面的教程中进行了叙述。尝试以下使用 tee 的命令,该命令允许将输出写入标准输出和文件

cd $FOAM_RUN/damBreak
foamRun | tee log

现在采用交互式运行代码,并将输出的副本存储在log文件中。

后处理

现在可按常用的方式对结果进行后处理。用户可以及时监测相分数alpha.water的变化情况,例如下图:

并行运行

上一个示例的结果是使用相对粗糙网格生成的, 现在希望增加网格分辨率并重新计算。 使用更精细的网格进而展示 OpenFOAM 的并行处理能力。

首先复制damBreak案例,例如

run
foamCloneCase damBreak damBreakFine

进入新的case目录(路径:run/damBreakFine/system/),将blockMeshDict字典中的块描述更改为

  blocks
  (
           hex (0 1 5 4 12 13 17 16) (46 10 1) simpleGrading (1 1 1) 
           hex (2 3 7 6 14 15 19 18) (40 10 1) simpleGrading (1 1 1) 
           hex (4 5 9 8 16 17 21 20) (46 76 1) simpleGrading (1 2 1) 
           hex (5 6 10 9 17 18 22 21) (4 76 1) simpleGrading (1 2 1)
           hex (6 7 11 10 18 19 23 22) (40 76 1) simpleGrading (1 2 1)
  );

此处显示从 blockMeshDict 文件中摘取的内容。需要更改设置以调整网格密度,例如设置 (46 10 1) 和 (1 2 1)。随后,运行 blockMesh 生成网格。

由于对网格进行了修改,因此damBreakFine案例里包含一些与新网格不一致的元素,所以需要在0时间目录中重新初始化相场alpha.water。需要注意,因为U和p_rgh场指定为均一值,与场中的单元数量无关,所以U和p_rgh场不需更改。

然后,重新运行setFields。首先,由于现在的网格大小与0目录中alpha.water文件中的元素数量不一致,必须删除alpha.water文件,以便setFields使用原始的alpha.water.orig文件。

rm 0/alpha.water
setFields

使用区域分解法进行并行计算。其中,将几何形状和相关场分解成多个部分,并分配给单独的处理器进行求解。因此,运行并行案例的第一步是使用decomposePar程序对区域进行分解。与decomposePar相关的是名为decomposeParDict的字典,其位于system目录中。此外,OpenFOAM中的etc目录中也可以找到案例字典,可以通过运行foamGet脚本将其复制到案例目录中:

foamGet decomposeParDict

输入上面命令后弹出下列内容:

Multiple files with "decomposeParDict" prefix found:
1) /home/dyfluid/OpenFOAM/OpenFOAM-11/etc/caseDicts/annotated/decomposeParDict
2) /home/dyfluid/OpenFOAM/OpenFOAM-11/etc/caseDicts/preProcessing/decomposeParDict
** Note: it is easier to use files NOT in the "annotated" directory

Enter file number (1-2) to obtain description (suggest 2) : 2
Copying /home/dyfluid/OpenFOAM/OpenFOAM-11/etc/caseDicts/preProcessing/decomposeParDict to system
*****************************

建议输入2。然后再system文件夹下面,会出现decomposeParDict的字典。其中的numberOfSubdomains指定案例分解成子域的数量,通常与可用的处理器数量相对应。比如2个核求解,numberOfSubdomains=2;8个核求解,numberOfSubdomains=8。

本案例使用simple method。它需要根据以下步骤来常规配置simpleCoeffs。域在x、y和z方向上分割成多个部分或子域,每个方向上的子域数量由向量n给出。由于此几何形状是二维的,因此第三个方向z无法分割,所以 \(n_z\)等于1。n的\(n_x\)\(n_y\)分量必须指定在x和y方向上分割域,这样\(n_x\)\(n_y\)指定的子域数量可以等于指定的numberOfSubdomains,即\(n_xn_y\) = numberOfSubdomains。

推荐将与子域相邻的单元面的数量保持在最低限度,因此对于方形几何形状,最好保持x和y方向之间的分割均匀。例如,假设希望在4个处理器上运行。会将numberOfSubdomains设置为4,n设置为(2, 2, 1)。用户应通过以下命令运行decomposePar。

decomposePar

终端输出显示分解的子域之间为均匀分布。分解将每个子域的网格和字段写入名为processor的独立子目录中,其中N是子域ID,例如0、1、2等。需要注意,应列出案例目录中的文件,以确认存在四个目录:processor0、processor1、processor2和processor3。

本案例展示了使用标准message-passing interface(MPI)的openMPI执行并行计算。若要使用多核心CPU的4个内核执行并行计算,可运行以下命令。

mpirun -np 4 foamRun -parallel

有关如何并行运行案例的更多详细信息,可以查阅第 3.4 节获取。例如,通过创建一个文件在网络上的更多节点上运行,该文件列出了要运行案例的机器的主机名,见第 3.4.3 节。

对并行计算的案例进行后处理

当案例并行运行时,结果将写入到processor 子目录内的时间目录中。通过列出processor0 目录的时间目录进行确认:

ls processor0

通过将单独的处理器目录视为一个独立的案例,进而对单个子域进行后处理。例如,要在ParaView中查看processor1子域,运行以下命令来启动paraFoam:

paraFoam -case processor1

图 2.18 显示该子域的网格,使用simple方法对域进行分解。

或者,对分解后的场和网格进行重新组合,以便进行串行后处理。reconstructPar工具从处理器子域的时间目录中获取字段文件,并为整个域构建等效的场文件。运行如下命令:

reconstructPar

重构后的场被写入案例目录中的解决方案时间目录。这些字段可以在ParaView中正常可视化。细网格的结果如图2.19所示。可以看到,与粗网格相比,气液界面的分辨率有了显著提高。