3.2 VOF模型分析

本节主要介绍VFRC载荷、VOF模型分析的环境条件设置及实例的详细分析。进行VOF分析时，需要同时应用初始和边界的VFRC载荷，需要设置与外界条件相关的环境条件，VFRC公差及时间步长等。

3.2.1 VFRC载荷

进行VOF分析时，需要同时应用初始和边界VFRC载荷。在大多数情况下，默认值是合适的。

VOF分析主要用来处理设计自由表面的流动问题，FLOTRAN对体积组分VFRC使用对流算法来追踪自由表面的发展。每一个单元的VFRC值在0～1之间变化，这里0表示空或无效的单元，而1表示充满的流体单元。值在0～1之间表示相应的单元是部分充满的或表面单元，自由表面可由VFRC场的分布确定。

在FLOTRAN中，VOF算法首先按拉格朗日方程计算流体多边形的运动，然后利用两个多边形相交的计算几何模型确定VFRC通量。这种方法称为计算拉格朗日-欧拉对流再变换（CLEAR）方法。这种CLEAR-VOF方法明确的按时间追踪自由表面，因此，对VOF分析有必要激活瞬态求解选项。为保证准确，应当使用流线迎风Petrov-Galcrkin（SUPG）方法。当遇到收敛困难时，减小时间步长，并激活改进的惯性松弛（MIR）程序。

为了激活FLOTRAN中的VOF分析和模拟表面张力影响，使用下面的命令或GUI菜单路径：

（1）初始VFRC载荷

按通常的方式对速度和其他自由度设置初始条件。对于初始VFRC场，可以把每个单元定义成初始完整、局部完整或者是空的。初始设置局部单元时要使每个局部单元必须与至少一个完整单元和一个空单元相邻。在分析过程中，如果一个局部单元临近任何一个完整单元，它将会被重新设置成空单元。因此，在自由表面上方至少应包括一层空单元以使表面随时间发展。

为了设置初始体积分数场，需要使用以下一种：

默认的体积分数为零，故如果在问题区域流体已经初始化，就必须相应地设定VFRC场。

为了设置一个给定几何图形的初始体积分数，需要使用以下一种：

ICVFRC会自动计算所有选定单元的体积分数值，这些选定单元或者在给定几何图形内部，或者与之相交。体积分数值如下：

1）单元完全在几何图形内部时为1。

2）单元与几何图形的边界相交时分数等于位于几何图形内部的单元面积。

通常有两种有效地几何形状，圆形（CIRC）和椭圆形（ELPT）。对于一个圆或椭圆，VAL1和VAL2分别为中心的X和Y坐标。对于一个圆，VAL3为半径，而不使用VAL4。对于一个椭圆，VAL3和VAL4分别为X半轴和Y半轴。

为了显示初始的体积分数载荷，使用以下的操作：

PBF和PSF命令优先于PICE命令。

（2）边界VFRC载荷

按通常的方式对速度和其他自由度设置边界条件。尽管流体区域外部的边界条件最初没有物理意义，但是如果流体接触到这些边界，这些条件就开始起作用。相应地，在通过流体潜在接触的所有边界上施加所有适当的边界条件。FLOTRAN使用自然条件来处理不确定的边界。

VFRC载荷的边界条件由自由边界VFRC值和润湿状况组成。对于边界VFRC值，FLOTRAN仅允许0（空的）和1（完整的）值。默认情况下，边界VFRC值设为1（完整的）。当边界与表面单元相邻时，FLOTRAN使用边界VFRC值来确定VOF重建阶段的垂直表面方向。润湿情况可以确定流体能否被输送到问题区域。0值表明没有润湿边界，而1值表明有润湿边界。默认情况下，润湿情况设为0（没有润湿边界），而相应地边界VFRC值必须设为1才有意义。在一个没有润湿的边界上，并没有发生流体的输送或者空单元进入问题区域。边界VFRC值仅用于帮助确定垂直面方向。

为了设定边界VFRC值，可使用以下方法之一：

对于SFL命令，VAL1确定边界VFRC值（默认为1），而VAL2确定润湿状况（默认为0）。如果Lab=VFRC，一个VAL2设为1表明一个润湿边界。对于SFE命令，当KVAL=0或1时，VAL1确定边界VFRC值（默认为1）；而当KVAL=2时，VAL1确定边界润湿状况（默认为0）。如果Lab=VFRC和KVAL=2，VAL1设为1表明一个润湿边界。在润湿边界上，流体逆流可以保持相关单元的完整。

3.2.2 输入设置

当进行VOF分析时，可以设置与外界条件相关的环境条件、VFRC公差及时间推进方法。

当研究表明张力效应时，需对表面张力系数和壁面静态接触角设定合适的值。关于这些性质类型的信息可以参考FLOTRAN流体性质的设定。

（1）外界条件

对以下的变量设置流体区域外部的外界节点值：

●VX-U 速度；

●VY-V 速度；

●VZ-W 速度；

●TEMP温度；

●ENKE湍动能；

●ENDS湍动能耗散率。

使用以下方法之一：

当复合一个Dirichlet条件的边界节点位于流体区域外部时，这个节点上赋予的外界值可以忽略不计。一旦与之相关的单元变为局部完整或完整的，所指定的Dirichlet条件就会起作用。除了压力以外，外界值仅为了后处理而设。压力作为一个实际的边界条件。

不像其他的FLOTRAN分析那样，VOF分析的压力场同时包括动压和静压，并且动压平衡于重力加速度项。例如，在VOF分析中，一个考虑重力的平行均匀流动的等压线由具有与重力平衡的梯度的平行线组成，而在其他FLOTRAN分析中这些等压线保持恒定。

在VOF分析中，出口有一定的压力或者有自然条件。当出口压力确定时，其值不同于其他FLOTRAN分析中与静压非偏移（或高程水头）相关的压力值。对于自然边界条件，出口压力会自动调整，在这种情况下，如果存在一个稳定解，出口压力就会收敛于修正值。

（2）VFRC公差

理想情况下，VFRC值为1意味着一个完整的流体单元，其值为0意味着一个空的或无效单元。然而实际上由于不同类型的误差累积，0或1的精确值在VOF水平对流运算法则中很难维持。

因此，FLOTRAN引用较小的中止值来确定完整的或空的单元。如果VFRC值大于或等于VOFU，就认为单元完整单元，其中，VOFU为VOF水平对流运算法则的上公差；如果VFRC值小于或等于VOFL，就认为单元为空单元，其中，VOFL为VOF水平对流运算法则的下公差。默认情况下，VOFL和VOFU都设为10 ^-5 。

在FLOTRAN求解过程中，与一个完整单元相邻的一个空的单元在局部缩放网格中发生了急剧的变化，这使刚度矩阵变得更病态，并严重影响期望的质量。为了避免出现这种情况，对于那些VFRC值小于或等于层流的LAML和湍流的TRBL的单元，FLOTRAN引用了另一个公差组合成有限元方程，其中，LAML和TRBL分别为层流和湍流的下偏差。类似地，当组合有限元方程时，VFRC值大于或等于层流的LAMU和湍流的TRBU的单元则认为是完整单元。

这些公差中的任何一个都可以设置。默认情况下，LAML和LAMU可设为10 ^-2 ，TRBL和TRBU设为10 ^-1 。

设置任一公差时，可使用以下方法之一：

（3）VFRC时间步

由于CLEAR-VOF在时间上是显性的，对于VOF水平对流运算法则必须使用一个时间步，这相当于一个小于1的CFL数。另一方面，由于FLOTRAN瞬态求解运算法则在时间上是隐式的，并且时间步大小上没有稳定性限制，因此，基于VOF约束上选择时间步有时可能会非常受限制。相应地，FLOTRAN每一求解步采用多重VOF水平对流步。VOF水平对流时间步可设为等于通过指定一个参数划分的求解时间步。如果VOF水平对流时间步仍然很大，FLOTRAN就会把VOF时间步自动减小1/2。VOF时间步自动减少会持续到在VOF水平对流计算过程中VFRC场的局部不平衡性小于COFL公差。通过消除一些时间步的检测，一个良好的参数估计会使计算效率更高。

设置每一求解时间步的VOF水平对流时间步的数目时，可使用以下方法之一：

对于一个VOF分析，密度为常数，因此，质量守恒相当于总流体体积（或二维下的面积）守恒。理想情况下，进入问题区域和流出问题区域的体积之差应等于区域内总体积的增加或减少，然而在实际有限元分析中，用Galerkin弱结构表示满足连续性。因此，一般每个单元都有一个较小的质量不平衡。

局部质量不平衡通常与离散误差相似，而相似误差为VOF分析中所引用的局部VFRC不平衡性的一个主要来源。尽管与总体积相比，这个不平衡可能会非常小，但是它会随着时间步的增加呈指数积累。因此，在每个VOF时间步上FLOTRAN对VFRC值进行了局部调整，通过按比例地增加或减少局部单元的VFRC值保证局部VFRC平衡。在Jobname.PFL文件中列出了最终的VFRC不平衡。

3.2.3 后处理

体积组分在结果文件单元数据部分作为不可求和的其他记录。如果查看体积组成，使用以下命令：

因为体积组分是必须以单元为基础的量，它不连续地穿过单元。如果希望见到光滑和连续的体积组合穿过，首先应当将结果存入单元表，使用以下命令：

然后，画出平均的结果，使用以下命令：

然而，不用之前的命令，也可以使用以下命令：

PLVFRC宏列出ETABLE和PLETAB命令。当CONT=0时，PLVFRC使用0.5和1的轮廓水平指出自由面和流体区域。当CONT=1时，使用当前的轮廓设置。

也可以对结果做动画。通过在时间步之间插值获得一段时间上的动画，使用以下的PLVFRC命令：

为了产生基于时间步的动画，而不用插值，使用以下的PLVFRC命令：

如果CONT=0，动画基于0.5和1两个轮廓水平。如果CONT=1，动画基于AUTOCNTRKY设置。

实例3-3——水坝的VOF分析

在此例中，两垂直的墙壁限制了一个矩形静水柱。当去掉右面的墙壁时，水将流出到刚性的水平面上，ICE命令可以设置区域的初始容积率。计算模型如图3-37所示，其单位为m。

【光盘文件】

——附带光盘“Ch3\实例3-3_start”

——附带光盘“Ch3\实例3-3_end”

——附带光盘“AVI\Ch3\3-3.avi”

【操作步骤】

（1）设置分析选项

选择Main Menu→Preference命令，弹出如图3-38所示的Preferences for GUI Filtering对话框，在Individual discipline(s) to show in the GUI栏中选择FLOTRAN CFD选项，单击OK按钮。

（2）定义工作文件名

选择Utility Menu→File→Change Jobname命令，弹出如图3-39所示的Change Jobname对话框，在Enter new jobname后面的文本框中输入“Example 3-3”，选择“New log and error files？”复选框为Yes，单击OK按钮。

（3）定义工作标题

选择Utility Menu→File→Change Title命令，弹出如图3-40所示的Change Title对话框，在Enter new title后面的文本框中输入“CFD analysis of the dam”，单击OK按钮。

（4）定义单元类型

选择Main Menu→Preprocessor→Element Type→Add/Edit/Delete命令，弹出Element Types对话框。单击Add按钮，弹出如图3-41所示的Library of Element Types对话框，在左面的列表栏中选择FLOTRAN CFD，在右面的列表栏中选择2D FLOTRAN 141，单击OK按钮，再单击Close按钮，完成单元类型的设置。

（5）生成分析区域的几何面

1）选择Main Menu→Preprocessor→Modeling→Create→Areas→Rectangle→By Dimensions命令，弹出如图3-42所示的Create Rectangle by Dimensions对话框。在“X1，X2 X-coordinates”文本框中依次输入“0”，“6”，在“Y1，Y2 Y-coordinates”文本框中依次输入“0”，“3”，单击OK按钮。

2）选择Utility Menu→PlotCtrls→Numbering命令，弹出如图3-43所示的Plot Numbering Controls对话框，单击LINE Line numbers后面的单选按钮为On，单击OK按钮。生成的几何模型如图3-44所示。

（6）定义单元形状

设置单元尺寸。选择Main Menu→Preprocessor→Meshing→Size Cntrls→ManualSize→Lines→Picked Lines，选择上下两条直线L1和L3，在选择菜单中单击Apply按钮。在弹出如图3-45所示的Element Sizes on Picked Lines对话框中，在SIZE Element edge length文本框中输入“0.2”，单击Apply按钮。选择左右两条直线L2和L4，并单击OK按钮，在弹出如图3-45所示的Element Sizes on Picked Lines对话框中，在SIZE Element edge length文本框中输入“0.05”，单击OK按钮。

（7）划分有限元网格

选择Main→Preprocessor→Meshing→MeshTool，在弹出的MeshTool对话框中，采取映射网格划分，划分四边形网格，单击Mesh按钮，在弹出的Mesh Areas的拾取框中，单击Pick All按钮，单击OK按钮。生成网格如图3-46所示。

（8）施加边界条件

1）单击Utility Menu→Select→Everything命令。

2）单击Utility Menu→Plot→Lines命令，重新绘制直线。

3）单击Utility Menu→Select→Entities命令，弹出如图3-47所示的Select Enti…对话框，在第一个下拉列表框中选择Nodes，在第二个下拉列表框中选择By Location，在第三个单选框中选择X coordinates，在Min，Max文本框中输入“0”，单击OK按钮。

4）单击Utility Menu→Select→Entities命令，弹出如图3-47所示的Select Enti…对话框，在第一个下拉列表框中选择Nodes，在第二个下拉列表框中选择By Location，在第三个单选框中选择X coordinates，在Min，Max文本框中输入“0”，在第四个单选框中选择Also Select，单击OK按钮。

5）选择Main Menu→Preprocessor→Loads→Define Loads→Apply→Fluid/CFD→Velocity→On Nodes命令，弹出如图3-48所示的Apply VELO load on Nodes对话框，在VX Load value文本框中输入“0”，在VY a load value文本框中输入“0”，单击OK按钮。

6）选择Utility Menu→Select→Everything命令。

7）选择Utility Menu→Select→Entities命令，弹出如图3-47所示的Select Enti…对话框，在第一个下拉列表框中选择Nodes，在第二个下拉列表框中选择By Location，在第三个单选框中选择X coordinates，在Min，Max文本框中输入“0，1.5”，单击OK按钮。

8）选择Utility Menu→Select→Entities命令，弹出如图3-47所示的Select Enti…对话框，在第一个下拉列表框中选择Nodes，在第二个下拉列表框中选择By Location，在第三个单选框中选择Y coordinates，在Min，Max文本框中输入“0，2.5”，单击OK按钮。

9）选择Utility Menu→Select→Entities命令，弹出如图3-47所示的Select Enti…对话框，在第一个下拉列表框中选择Elements，在第二个下拉列表框中选择Attached to，在第三个单选框中选择Nodes all，单击OK按钮。

10）选择Main Menu→Preprocessor→Loads→Define Loads→Apply→Fluid/CFD→Volume Fract→Init Loads→On Elements命令，在弹出的Apply VFRC on Elems对话框中，单击Pick All按钮，弹出如图3-49所示的Initialize VFRC on Elems对话框，在VFRC文本框中输入“1”，单击OK按钮。

11）选择Main Menu→Preprocessor→Loads→Define Loads→Apply→Initial Condit’n→Define命令，在弹出的“Defined Initial C…”选择对话框中单击Pick All按钮，弹出如图3-50所示的Define Initial Conditions对话框，在Lab DOF to be specified选择列表中选择VX，在VALUE Initial value of DOF文本框中输入“0”，单击Apply按钮会再弹出Defined Initial C…选择对话框，在对话框中单击Pick All按钮，弹出如图3-50所示的Defined Initial Conditions对话框，在Lab DOF to specified选择列表中选择VY，在VALUE Initial value of DOF文本框中输入“0”。加载后的几何模型如图3-51所示。

（9）求解

1）选择Main Menu→Solution→FLOTRAN Set Up→Flow Environment→Gravity命令，弹出如图3-52所示的Gravity Specification对话框，在ACELY Accel in Y direction文本框中输入“9.81”，单击OK按钮。

2）选择Main Menu→Solution→FLOTRAN SetUp→Fluid Properties命令，弹出Fluid Properties对话框，采用默认设置，单击OK按钮，弹出如图3-53所示的CFD Flow Properties对话框，在Density property type CONSTANT Constant value文本框中输入“1000”，在Viscosity property Types CONSTANT Constant value文本框中输入“0.001”，其余选项采用默认设置，单击OK按钮。

3）选择Main Menu→Solution→FLOTRAN Set Up→Solution Options命令，弹出如图3-54所示的FLOTRAN Solution Options对话框，在“TRAN Steady state or transient？”选择列表中选择Transient，单击“VOF Activate VOF advection？”后的方框，使No状态变为Yes，单击OK按钮。

4）选择Main Menu→Solution→FLOTRAN Set Up→Execution Ctrl命令，弹出如图3-55所示的Time Step Controls对话框，在STEP Select time step control后的选择列表框中选择User Defined。单击OK按钮，弹出如图3-56所示的Transient Controls对话框，在STEP User-defined time step文本框中输入“0.01”，在NUMB Number of time steps文本框中输入“70”，在GLOB Global iter per time step 文本框中输入“5”，在STEP，APPE Number of time steps文本框中输入“0.1”，单击OK按钮。

5）选择Main Menu→Soltion→FLOTRAN Set Up→Advction命令，弹出如图3-57所示的Advection Option对话框，确认MOME Momentum equations后要选择SUPG，单击OK按钮。

6）选择Utility Menu→Select→Everything，选择Utility Menu→file→Save as命令，弹出Save Dtabase对话框，在Save Dtabase to文本框中输入“3-3_start.db”，保存上述操作过程，单击OK按钮关闭该对话框。

7）求解。选择Main Menu→Solution→Run FLOTRAN，开始进行求解。当弹出Solution is Done对话框时，表示求解过程已经结束，单击Close按钮，关闭该对话框。

8）选择Utility Menu→file→Save as命令，弹出Save Dtabase对话框，在Save Dtabase to文本框中输入“3-3_end.db”，保存求解结果，单击OK按钮关闭该对话框。

（10）后处理

1）调入最后迭代步的结算结果。选择Main Menu→General Postproc→Read Results→Last Set。

2）以矢量方式显示速度分布结果。选择Main Menu→General Postproc→Plot Results→Vector Plot→Predefined，弹出如图3-58所示的对话框。在Vector item to be plotted后面的第一列表中选择DOF solution，在第二列表中选择Velocity V，单击OK按钮，所得到的结果如图3-59所示。

3）绘制压力等值线图。选择Main Menu→General Postproc→Plot Results→Contour Plot→Nodal Solu，弹出如图3-60所示的对话框，在Item to be contoured后面的两个列表框中分别选择DOF solution和Pressure，单击OK按钮，所得到的结果如图3-61所示。

4）绘制速度等值线图。重复上步操作，在Item to be contoured后面的两个列表框中分别选择DOF solution和Fluid velocity，单击OK按钮，所得到的结果如图3-62所示。