14.3.1 VOF模型概述

14.3.1 VOF模型概述

VOF模型通过求解一组动量方程并跟踪每一相在计算区域内的体积分数，可以模拟两个或多个不混溶流体。VOF模型的典型的应用包括射流破裂的预测、液体中大气泡的运动、溃坝后液体的运动以及包含气液界面的稳态或瞬态跟踪。

14.3.2 VOF模型使用限制

Fluent中的VOF模型存在以下限制：

只能用于压力基求解器。密度基求解器无法使用VOF模型

所有控制体积必须填充单个流体相或多个相的组合。VOF不允许模型中包含不存在任何类型流体的空洞区域

只能将其中某一相定义为可压缩理想气体。使用UDF指定可压缩液体则不受此限制

使用VOF模型时，无法模拟周期性流动(指定的质量流量或指定压降)

二阶隐式时间格式无法与VOF显式格式一起使用

当DPM模型与VOF模型一起使用时，不能使用Shared Memory方法。不过可以使用Message Passing或Hybrid方法

Coupled VOF Level Set模型不能用于多面体网格

VOF模型不能与非预混、预混、部分预混燃烧模型一起使用

14.3.3 稳态与瞬态VOF计算

Fluent中的VOF模型通常用于瞬态计算，但对于只关心稳态状态的问题也可以进行稳态计算。只有当求解结果与初始条件无关，且每一相都明确的流入边界时，稳态VOF计算才是合理的。例如，由于旋转杯子内自由表面的形状取决于流体的初始液位，因此必须采用瞬态计算来计算。另一个例子，顶部有空气区域和单独的进气入口的通道中的水的流动可以用稳态来求解。

VOF模型假设计算区域内的相互不渗透。对于添加到模型中的每一相，都会引入一个新的变量：计算网格中该相的体积分数。在每个网格中，所有相的体积分数之和为1。因此，根据体积分数的值，可以得到该网格内的相的分布。换句话说，如果某个网格内第q相的体积分数为αq​，则有：

αq​=0：网格内不存在q相

αq​=1：网格内全部为q相

0<αq​<0：网格内包含 q 相与其他相的分界面

14.3.4 体积分数方程

相界面的跟踪是通过求解一个(或多个)包含相的体积分数的连续性方程来完成的。对于第q相，方程具有以下形式：

ρq​1​[∂t∂​(αq​ρq​)+∇⋅(αq​ρq​vq​)]=Sαq​+p=1∑n​(m˙pq​−m˙qp​)

式中，m˙qp​为从q相到p相的传质量；m˙pq​为从p相到q相的传质量。默认情况下，方程式右侧的源项Sαq​为零，但用户可以为每一相指定常数或自定义质量源。

对于主相（Primary Phase)，不求解体积分数方程，主相的体积分数利用下式进行约束计算：

Σq=1n​αq​=1

体积分数方程可以通过隐式或显示时间格式进行求解。

14.3.4.1 隐式格式

使用隐式公式时，体积分数方程按以下方式离散：

Δtαqn+1​ρqn+1​−αqn​ρqn​​V+f∑​(ρqn+1​Ufn+1​αq,fn+1​)=[Sαq​+​p=1∑n​(m˙pq​−m˙qp​)]V

式中，n+1为当前时间步索引；n为前一个时间步索引；αqn+1​为时间步n+1时网格内的体积分数；αqn​为时间步n的网格内的体积分数；αq,fn+1​为时间步n+1网格面上第q相的体积分数；Ufn+1​为时间步n+1是通过网格面的体积通量；V为网格体积。

由于当前时间步长处的体积分数是当前时间步长处其他物理量的函数，因此在每个时间步长处迭代求解每个次相（Secondary Phase）体积分数的标量输运方程。

使用所选的空间离散格式对面通量进行插值。User Guide中体积分数的空间离散格式讨论了Fluent中隐式格式的可用方案。

隐式格式可用于瞬态和稳态计算。

14.3.4.2 显式格式

显式格式与时间有关，体积分数按以下方式离散：

Δtαqn+1​ρqn+1​−αqn​ρqn​​V+f∑​(ρq​Ufn​αq,fn​)=[p=1∑n​(m˙pq​−m˙qp​)+Sαq​​]V

式中，n+1为当前时间步索引；n为前一个时间步索引；αq,f​为网格面上第q相的体积分数；V为网格的体积分数；Uf​为基于法向速度通过网格面的体积。

由于当前时间步长的体积分数是根据前一时间步长的已知量直接计算的，因此显式格式不需要在每个时间步长内迭代求解输运方程。

可以使用界面跟踪或捕捉算法(例如Geo-Reconstruct、CICSAM、Compressive和 Modified HRIC来内插面通量。Fluent会自动优化体积分数方程积分的时间步长，但用户可以通过修改Courant数来影响此时间步长的计算。可以选择为每个时间步长更新一次体积分数，或为每个时间步长内的每个迭代更新一次体积分数。

重要提示：当使用显式格式时，必须计算瞬态。

14.3.4.3 界面附近插值

Fluent的控制体公式要求计算通过控制体积表面的对流和扩散通量，并与控制体内的源项平衡。

在Geometry Reconstruct和donor-acceptor方法中，Fluent对两相界面附近的网格进行了特殊的插值处理，线图显示了实际界面形状以及这两种方法在计算过程中假定的界面。

显式格式和隐式格式以与完全填充一个或另一个相的网格相同的插值方式处理这些网格格(即使用标准迎风格式(一阶迎风格式)、二阶格式(二阶迎风格式)、QUICK格式)、修改的HRIC格式(修改的HRIC格式)、COMPRESSION(压缩格式和基于界面模型的变体)或CICSAM格式(任意网格的压缩界面捕获方案(CICSAM))，而不是应用特殊处理。

14.3.4.1 几何重构算法

在几何重构方法中，使用Fluent 中的标准插值格式，在网格完全填充一相或另一个相时获得面通量。当网格靠近两相界面时，采用几何重构方法。

几何重构方法用分段线性方法表示流体之间的界面。在Fluent中此格式是最精确的方法，其适用于一般的非结构化网格。该方法假设两个流体之间的界面在每个网格内存在一个线性斜率，并使用这个线性形状来计算流体通过网格面的对流通量。

该重构方案的第一步是根据网格体积分数及其导数的信息，计算线性界面相对于每个部分填充网格中心的位置。第二步是使用计算出的线性界面表示及界面的法向和切向速度分布信息，计算通过每个表面的流体对流量。第三步是使用前一步计算的通量平衡计算每个网格中的体积分数。

重要提示：使用Geometry Reconstruction方法时，必须采用瞬态计算。此外，如果使用一致网格（即两个子域相交的边界处的网格节点一致），则必须确保计算域内没有无厚度壁面。

14.3.4.2 Donor-Acceptor算法

在Donor-Acceptor方法中，Fluent使用的标准插值方法来获得网格完全充满一相或另一相时的面通量。当网格靠近两相界面时，Donor-Acceptor方法用于确定通过表面的流体对流量。该方法将一个网格识别为来自一个相的一定量流体的供体，将另一个（相邻）网格识别为相同量流体的受体，并用于防止界面处的数值扩散。可以通过网格边界对流的一个相的流体量受到两个值中的最小值的限制：供网格中的填充体积或受网格中的自由体积。

界面的方向也用于确定面通量。界面方向可以是水平的，也可以是垂直的，这取决于网格内相的体积分数梯度的方向，以及共享面的相邻单元的体积分数梯度的方向。根据界面的方向及其运动，通量值可通过纯上风、纯下风或两者的某种组合获得。

重要提示：当使用Donor-Acceptor方法时，必须计算时间相关的解。此外，Donor-Acceptor方案只能用于四边形或六面体网格。如果使用一致网格（即两个子域相交的边界处的网格节点一致），则必须确保计算域内没有无厚度壁面。

14.3.4.3 CICSAM算法

基于Ubbink工作[629]的任意网格压缩界面捕获方法（The compressive interface capturing scheme for arbitrary meshes，CICSAM）是一种高分辨率差分格式。CICSAM方案特别适用于相之间具有高粘度比的流动。CICSAM作为显式格式在Fluent中实现，其优点是生成的界面几乎与几何重建方案一样尖锐。

14.3.4.4 Compressive方法

Compressive方法是基于斜率限制器的二阶重构方法。在空间离散化方法中使用斜率限制器，以避免由于解域中的急剧变化而在高阶空间离散化格式中可能出现的虚假振荡或摆动。下面的理论适用于分区离散化和相位局部化离散化，它们使用压缩格式的框架。

αf​=αd​+β∇αd​⋅dr

式中，αf​为面上VOF值；αd​为donor网格的VOF值；β为斜率限制值；∇αd​为donor网格的VOF梯度值；dr为网格到面的距离。

斜率限制器被限制为介于0和2（包括0和2）之间的值。对于小于1的值，空间离散由低分辨率格式表示。对于介于1和2之间的值，空间离散由高分辨率格式表示。斜率限制器的值及其离散格式如下表所示：

斜率限制器的值格式

0 first order upwind

1 以体积分数的全局最小值/最大值为界的second order reconstruction

2 Compressive

0<β<1及1<β<2 混合：其中0和1之间的值表示一阶和二阶的混合，1和2之间的值表示二阶和压缩方案的混合

压缩格式的离散化取决于界面区域类型的选择。当选择sharp界面区域建模时，压缩方案仅适用于sharp界面建模。然而，当选择sharp/dispersed界面建模时，压缩方案适用于sharp和dispersed界面建模。

14.3.4.5 有界梯度最大化（BGM）

BGM方法用于获得与VOF模型的尖锐界面，与几何重建方法获得的界面相当。目前，该方案仅适用于稳态求解器，不能用于瞬态问题。在BGM方法中，离散化以这样一种方式发生，即通过最大化面值对外推顺风值的加权程度，使梯度的局部值最大化[652]。

14.3.5 材料属性

输运方程中出现的性质由每个控制体积中存在的组分相确定。例如在两相系统中，如果相由下标1和2表示，并且如果跟踪第二个相位的体积分数，则每个网格中的密度由下式给出：

ρ=α2​ρ2​(1−α2​)ρ1​

一般而言，对于n相系统，体积分数平均密度呈现以下形式：

ρ=∑αq​ρq​

其他物性参数（如粘度）计算也采用相同的方式。

14.3.6 动量方程

在整个区域内求解单个动量方程，得到的速度场在各相之间共享。动量方程（如下所示）取决于所有相的体积分数。

∂t∂​(ρv)+∇⋅(ρvv)=−∇p+∇⋅[μ(∇v+∇vT)]+ρg​+F

共享场近似的一个限制是，在相间存在较大速度差的情况下，界面附近计算的速度精度可能会受到不利影响。

注意如果粘度比大于1000，可能会导致收敛困难。CICSAM适用于具有高相间粘度比的流动，因此解决了收敛性差的问题。

14.3.7 能量方程

能量方程在各相之间也是共享的。

∂t∂​(ρE)+∇⋅(v(ρE+p))=∇⋅(keff​∇T−q∑​j∑​hj,q​j​j,q​+(τˉeff​⋅v))+Sh​

其中， keff​为有效热导率，keff​=k+kt​，其中kt​为湍流热导率，根据所使用的湍流模型定义。Jj​为组分j的扩散通量；hj,q​为相q中的组分j的焓；Jj,q​为相q中的组分j的扩散通量。方程(9)右侧的前三项分别表示由于传导、组分扩散和粘性耗散引起的能量传递。Sh​包含已定义的体积热源，但不包括由有限速度体积或表面反应生成的热源，因为总焓计算中已经包含了组分形成焓。

VOF模型将能量E视为质量平均变量：

E=∑q=1n​αq​ρq​∑q=1n​αq​ρq​Eq​​

式中，

Eq​=hq​−ρq​p​+2v2​

其中，各相的hq​基于该相的比热及共享温度。

物性ρ、keff​（有效热导率）及μeff​通过各相的体积平均值计算得到。源项Sh​包括辐射及其他体积热源的贡献。

与速度场一样，在流体相之间存在较大温差的情况下，界面附近温度的精度受到限制。在物性变化几个数量级的情况下，也会出现此类问题。例如，如果模型包含液态金属与空气的组合，则材料介质的电导率可能相差多达四个数量级。这些物性上的巨大差异导致方程组具有各向异性系数，这反过来又会导致收敛和精度下降。

14.3.8 附加标量方程

根据问题定义，求解方案中可能会涉及其他的标量方程。在湍流求解的问题中求解一组单一的输运方程，湍流变量（例如，k、ε和或雷诺应力）由计算区域中的各相共享。

14.3.9 表面张力及附着

在Fluent中，可以在模拟中沿一对相之间的界面包含表面张力的影响。可以将表面张力系数指定为常数、通过多项式指定为温度的函数或通过UDF指定为任何变量的函数。求解器将包括由于表面张力系数变化而产生的附加切向应力项（导致所谓的Marangoni对流）。可变表面张力系数效应通常仅在零重力/接近零重力条件下才重要。

壁面附着效应可以通过指定相与壁面之间的接触角以及多孔跳跃的方式来考虑。

14.3.9.1 表面张力

表面张力是由于流体中分子之间的吸引力而产生的。例如考虑水中的气泡，在气泡内相邻分子间的净作用力为零，然而在表面上其净力是径向向内的，整个球面上力的径向分量的综合作用是使表面收缩，从而增加表面凹面上的压力。表面张力是一种仅作用于表面的力，在这种情况下需要保持平衡。它的作用是平衡径向向内的分子间吸引力和穿过表面的径向向外的压力梯度力。在两种流体分离但其中一种流体不是球形气泡的区域，表面张力通过减小界面面积来最小化自由能。

在Fluent中存在两种表面张力模型：连续表面力（Continuum Surface Force，CSF）及连续表面应力（Continuum Surface Stress，CSS）。

注：三角形和四面体网格的表面张力效应的计算不如四边形和六面体网格的精确。因此，表面张力效应最重要的区域应该用四边形或六面体进行网格划分。

14.3.9.1.1 CSF模型

Brackbill等提出的连续表面力（CSF）模型将表面张力解释为穿过界面的连续三维效应，而不是界面上的边界值条件。通过在动量方程中添加源项来模拟表面张力效应。为了理解源项的起源，考虑表面张力沿表面是恒定的这一特殊情况，并且只考虑界面的法向力。可以证明，跨越表面的压降取决于表面张力系数及由两个正交方向的表面曲率半径R1和R2。

p2​−p1​=σ(R1​1​−R2​1​)

式中，p1​、p2​为界面两侧的两种流体的压力。

表面曲率通过界面上的表面法线的局部梯度计算得到。假设表面法向n定义为q相的体积分数αq​的梯度:

n=∇αq​

曲率κ定义为单位法向n^的散度：

κ=∇⋅n^

其中：

n^=∣n∣n​

表面张力可以用穿过表面的压力jiey阶跃来表示。使用散度定理可以将表面上的力表示为体积力，该体积力为增加到动量方程中的源项。其形式如下：

Fvol​=ij,i

此表达式允许在存在两个以上相的网格附近平滑地叠加力。如果一个网格中只有两相，则有ki​=−ki​及∇αi​=−∇αi​，此时上式可以简化为：

Fvol​=σij​21​(ρi​+ρj​)ρκi​∇αi​​

式中，ρ为体积平均密度。上式表明网格表面张力源项与网格内的平均密度成正比。

14.3.9.1.2 CSS模型

与连续表面力（CSF）方法的非守恒表述不同,连续表面应力（CSS）方法是以守恒的方式对表面张力进行建模的另一种方法。CSS避免了显式的曲率计算，其可以表示为基于表面应力的毛细管力建模的各向异性变体。

在CSS方法中，由表面张力引起的表面应力张量表示为：

T=σ(I−n^⊗n^)∣n∣

n=∇α

n^=∣n∣n​

式中，I为单位张量；σ为表面张力系数；⊗为向量的张量积；α为体积分数；n为体积分数梯度。

T=σ(∣∇α∣I−∣∇α∣∇α⊗∇α​)

则表面张力可表示为：

FCSS​=∇⋅T

14.3.9.1.3 CSS模型与CSF模型比较

与CSF方法相比，CSS方法几乎没有额外的优势，尤其是对于涉及可变表面张力的情况。由于压力梯度和表面张力的不平衡，CSS和CSF方法都会在界面处引入parasitic currents 。

在CSF方法中，表面张力以非守恒方式表示，如下所示：

FCSF​=σκ∇α

式中，κ为曲率。此表达式仅在表面张力恒定的情况下有效。

对于可变表面张力，CSF模型要求根据表面张力梯度在界面的切向方向上模拟一个附加项。

在CSS模型中，表面张力以守恒的形式表示：

FCSS​=∇⋅[σ(∣∇α∣I−∣∇α∣∇α⊗∇α​)]

CSS方法不需要对曲率进行任何显式计算。因此其在低分辨率区域（例如锐角）的物理性能更高。

CSS方法由于其表达式守恒，因此不需要任何额外的项来模拟可变表面张力。

14.3.9.1.4 何时需要考虑表面张力

表面张力效应的重要性取决于两个无量纲量的值：雷诺数Re和毛细管数Ca；或者雷诺数Re和韦伯数We。

当ℜ≪1时，采用毛细管数Ca进行考虑：

Ca=σμU​

当ℜ≫1时，采用韦伯数进行考虑：

We=σρLU2​

式中，U为自由流速度；当Ca≫1 或 We≫1时，可以忽略表面张力。

14.3.9.2 壁面附着

在CSF模型框架内，可以很容易地估计在平衡状态下与刚性边界接触的流体界面处的壁面粘附力的影响。CSF模型通过假设流体与壁面之间的接触角用于调整壁面附近网格中的表面法线，而非将此边界条件施加在壁面本身。这种所谓的动态边界条件导致了近壁表面曲率的调整。

若θw​为壁面接触角，则临近壁面处的活动网格的曲面法向为：

n^=n^w​cosθw​+t^w​sinθw​

式中n^w​及t^w​分别为壁面的法向单位向量与切向单位向量。该接触角与距离壁面一个网格的通常计算的表面法向的组合确定表面的局部曲率，该曲率用于调整表面张力计算中的体力项。

下图显示了在接触角为90度的情况下使用壁面粘附时(A)和不使用壁面粘附时(B)自由表面的位置。

如果使用90度接触角(A)的壁面粘附，那么在曲率计算中，接触角的边界条件迫使界面垂直于边界。

如果不使用壁面粘附(B)，那么在计算曲率时，壁面处的体积分数梯度从邻近壁面的网格复制，且接触角没有边界条件强制。

14.3.9.3 Jump Adhesion

与墙壁附着力类似，在使用VOF模型时，也可以选择提供跳跃附着力。这里接触角处理适用于多孔跳跃边界的两侧，假设两侧的接触角相同。

因此，如果θw​为多孔阶跃的接触角，则与多孔阶跃相邻的网格的面法向为：

n^=n^w​cosθw​+t^w​sinθw​

式中，n^w​与t^w​分别为多孔阶跃的法向单位向量与切向单位向量。

ANSYS Fluent提供了两种在多孔跳跃边界处处理阶跃附着的方法：

Constrained Two-Sided Adhesion Treatment

该选项选项会对粘附处理施加约束。这里接触角处理将仅应用于多孔性跳跃与非多孔性流体区域相邻的一侧。因此，接触角处理不会应用于与多孔介质区域相邻的多孔性跳跃的一侧。如果禁用此选项，则接触角处理将应用于多孔层跳跃的所有侧。此选项为默认选项。

Forced Two-Sided Adhesion

FLUENT允许对流体区使用强制双侧接触角处理，而不会施加任何约束。

14.3.10 明渠流

Fluent可以使用VOF模型和明渠边界条件模拟明渠水流的影响（例如河流、水坝及无边界流动中的表面穿透结构）。这些流动涉及流动的流体与其上方流体（通常为大气）之间存在自由表面。在这种情况下，波传播与自由表面行为变得很重要。明渠流动通常由重力和惯性力控制。此功能主要适用于海洋应用和通过排水系统的流动分析。

明渠流的特征是无量纲弗劳德数，其定义为惯性力与静水力的比值。

Fr=gy​V​

其中，V为速度；g为重力加速度；y为长度尺度，在明渠流中，通常指从通道底部到自由表面的距离。上式中的分母为波的传播速度，静止参考系中观察的波速为：

Vw​=V±gy​

基于弗劳德数，明渠流可分为三类：

当Fr<1，此时V

当Fr=1，此时Vw​=gy​，流动为临界流动，上游传播的波保持静止，在这种情况下，流动特征会发生变化。

当Fr>1，此时Vw​>gy​，流动为超临界流动，扰动不能向上游传播。在这种情况下，下游条件不会影响上游流动。

14.3.10.1 上游边界条件

明渠流中可以使用两种类型的上游边界条件：

Pressure Inlet

Mass Flow Rate

14.3.10.1.1 Pressure Inlet

入口总压通过下式计算得到：

p0​=21​ρV2+(ρ−ρ0​)∣g​∣(g^​⋅(b−a))

式中，b与a分别为网格面质心的位置向量与自由面上任意点的位置向量，这里自由面被假定为水平的，且垂直于重力的方向。g​为重力向量；∣g​∣为重力大小；g^​为重力方向单位向量；V为速度值；ρ为网格内混合物的密度；ρ0​为参考密度。

动压q可通过下式计算：

q=2ρ​V2

静压ps​计算为：

ps​=(ρ−ρ0​)∣g​∣(g^​⋅(b−a))(14-40)

可以将其展开为：

ps​=(ρ−ρ0​)∣g​∣((g^​⋅b)+ylocal​)(14-41)

其中，从自由表面到参考位置的距离ylocal​可表示为：

ylocal​=−(a⋅g^​)(14-42)

14.3.10.2 Mass Flow Rate

与明渠流相关联的每相的质量流量定义为：

m˙phase ​=ρphase ​( Area phase ​)( Velocity )(14-43)

14.3.10.1.3 Backflow Volume Fraction Specification

在明渠流中，Fluent根据Boundary Conditions对话框中指定的输入参数在内部计算体积分数，因此该选项被禁用。

对于亚临界入口流(Fr < 1)， Fluent利用相邻网格的数值在边界上重构体积分数值。这可以通过以下步骤来完成:

利用网格值计算边界处体积分数的节点值

利用节点值插值计算每个边界面的体积分数

对于超临界进口流(Fr > 1)，边界上的体积分数值可以用距离底部固定高度的自由表面来计算。

14.3.10.2 下游边界条件

14.3.10.2.1 Pressure Outlet

静压值取决于压力指定方法：

Free Surface Level：自由面静压通过式(14-40)与(14-42)计算得到。对于亚临界出口流(Fr <1)，静压取自边界上指定的压力分布，否则取自相邻网格的压力。对于超临界流体(Fr > 1)，压力总是来自相邻的网格。

From Neighboring Cell：静压值取自相邻网格

Gauge Pressure：静压值由用户指定。

14.3.10.2.2 Outflow边界

可以在明渠流出口处使用outlet边界条件来模拟流速和压力细节在求解之前未知的流动出口。如果流出边界处的条件未知，则Fluent将从计算区域内部推断所需信息。

然而，了解这种边界类型的局限性很重要：

只能在出口处使用单个流出边界，这是通过将flow rate weighting设置为1来实现的。换句话说，在具有流出边界的明渠流中，不允许出现流出分流。

模拟中应存在初始流场，以避免因流动出口处的回流而导致收敛问题，这将导致不可靠的计算结果。

outlet边界条件只能与mass-flow inlet边界一起使用，其与Pressure-inlet和Pressure-outlet边界不兼容。例如，如果选择pressure-inlet作为入口边界，则只能在出口处使用Pressure-outlet边界条件。如果入口选择mass-flow-inlet边界，则可以在出口处使用outflow或pressure-outlet边界条件。注意，这仅适用于明渠流。

注意，outflow边界条件假设在垂直于流出边界面的方向上流动充分发展。因此，应相应放置此类表面。

14.3.10.2.3 Backflow Volume Fraction Specification

由于使用相邻网格值在内部计算出口边界上的体积分数值，因此此选项处于禁用状态。

14.3.10.3 数值沙滩处理

在某些应用中，需要抑制由通过波的出口边界引起的数值反射。为避免波浪反射，在压力出口边界附近的网格区域的动量方程中增加阻尼汇项[482]，[710]。

S=−[C1​ρV+21​C2​ρ∣V∣V]f(z)f(x)

式中，z^为沿重力竖直方向；x^为流动方向；S为z^方向的动量汇项；C1​为线性阻尼阻力(1/s)，默认值为10；C2​为二次阻尼阻力(1/m)，默认值为10；V为沿z^方向的速度；z为离自由面的距离；x为沿流动方向x^的距离；f(x)为x^方向的阻尼函数；f(z)为z^方向的阻尼函数。

x^与z^方向的缩放因子通过下式进行定义：

rX​=Xe​−XS​X−XS​​

rZ​=Zb​−Zfs​Z−Zfs​​

其中 x^与z^方向的阻尼函数定义为：

f(x)=(rx​)2

f(z)=1−rz​

式中，Xs​与Xe​为X^方向上阻尼区的起点与终点；Zfs​与Zb​为沿z^方向的自由面与底部标高。

14.3.11 明渠波边界条件

明渠波浪边界条件允许您模拟规则/不规则波浪的传播，这在海洋工业中用于分析波浪运动学和波浪对移动体和海上结构的冲击载荷。

小振幅波理论一般适用于较低的波陡和较低的相对高度，而有限振幅波理论更适用于增加波陡或增加相对高度。波陡一般定义为波高与波长之比，相对高度定义为波高与水深之比。

Fluent 为通过速度入口边界条件的入射表面重力波提供以下选项：

First order Airy波理论，适用于从浅到深的水位范围内的小振幅波，其本质上是线性的。

高阶Stokes波理论，适用于中到深的水位范围内的有限振幅波，本质上是非线性的。

高阶 Cnoidal/孤波理论，适用于浅水位范围内的有限振幅波，本质上是非线性的。

线性波叠加，用于产生干扰、betas、驻波、不规则波等各种物理现象。

长波/短波谱，用于基于波能量分布函数对中深水位范围内的非线性随机波进行建模。

每种波浪理论的短重力波表达式基于无限水深，而浅波或中波表达式基于有限水深。

波高H定义为：

H=2A=At​+Ac​

式中，A为波幅，At​为波谷处的波幅；Ac​为波峰处的波幅。对于线性波理论At​=Ac​；对于非线性波理论At​=Ac​。

波数k定义为：

k=λ2π​

式中λ为波长。

波数的向量形式为：

K=kX​X^+ky​y^​

式中，X^为参考波传播方向；Z^为重力相反的方向；y^​为X^及Z^的法线方向。X^及y^​方向的波数定义为：

kx​=kcosθ

ky​=ksinθ

其中θ为波航向角（wave heading angle），其定义为波阵面与参考波波传播方向之间的角度，在x^与y^​平面内。

有效波频率ωe​定义为：

ωe​=ω+K⋅U

其中ω为波频率，U为平均流速。当相对于移动物体的运动指定流动时，移动物体的影响也可以与流动结合。

波速c定义为：

c=kω​

通过叠加所有速度分量获得的入射波的最终速度矢量表示为：

V=U+ux^+vy​+wz^

式中，u、v、w分别为x^、t^、z^方向的速度分量。

变量α用于所有的波浪理论，其定义为：

α=kx​x+ky​y−we​t+ε

其中，x、y分别为x^、y^​方向的空间坐标；ε为相位差；t为时间。

14.3.11.1 Airy波理论

线性波的波形分布为：

ζ(X,t)=Acosα

浅波/中波的波频率ω定义为：

ω=gktanh(kh)​

短重力波的波频率定义为：

ω=gk​

式中，h为水深；k为波数；g为重力加速度。

入射波边界条件的速度分量可以用浅/中波和短重力波来描述。

1、浅/中波的速度分量

(uv​)=ωgkA​cosh(kh)cosh[k(z+h)]​(cosθsinθ​)cosα

w=ωgkA​cosh(kh)sinh[k(z+h)]​sinα

2、短重力波的速度分量

(uv​)=ωgkA​ekz(cosθsinθ​)cosα

w=ωgkA​ekzsinα

其中z是z^方向上自由表面水平面的标高，与重力方向相反。

14.3.11.2 Stokes波理论

Fluent 根据 John D. Fenton [165] 的工作实现了 Stokes 波浪理论。此波浪理论适用于在中到深的水深范围内运行的高陡度有限振幅波。

高阶斯托克斯理论（二阶到五阶）的波分布的一般表达式为：

ζ(X,t)=k1​i=1∑n​j=1∑i​bij​ξicos(jα)

相关速度势的一般表达式为：

Φ(X,t)=k1​kg​tanh(kh)​i=1∑n​ξij=1∑i​aij​cosh(jjk(z+h))sin(jα)

式中，

ξ=2kH​

其中，ξ为波陡度；n为波浪理论指数（2 到 5：分别从 2 阶斯托克斯到 5 阶斯托克斯）。

波速c通过下式确定：

c=kg​tanh(kh)​(1+c3​ξ2+c5​ξ4)

对于二阶斯托克斯波，c3​=c5​=0，因此二阶斯托克斯波使用与一阶斯托克斯博相同的分散关系；对于三阶及四阶斯托克斯波，c5​=0。

aij​、bij​、Ci​为kH的负数表达式。

波频率ω定义为：

ω=kc

表面重力波的速度分量来自速度势函数：

u=∂x∂Φ​cosθ

v=∂y∂Φ​sinθ

w=∂z∂Φ​

14.3.11.3 弦波/孤波理论

Fluent 根据 John D. Fenton (1998) [166] (p. 1004) 的工作实现了使用复雅可比和椭圆函数表达的 Cnoidal/孤波理论。 cnoidal计算结果在浅水中显示出长而平坦的波谷和狭窄的波峰。在无限波长的限制下，cnoidal 解将孤立波描述为具有单个波峰的波，没有波谷。由于连续波理论的复杂性，孤波理论更广泛地用于浅水区。

注：为简单起见，作为相对高度函数的高阶项已从下面的数值细节中省略。有关详细的数值表达式，请参阅 John D. Fenton [166] (p. 1004) 的工作。

椭圆函数参数m , 通过求解以下非线性方程来计算：

hλ​=4K(m)h3H​​[1+ Higher ⋅ order ⋅ terms ]

从底部算起的高度htr​由以下关系式得到：

hhtr​​=1−hH​K(m)E(m)​+ Higher ⋅ order ⋅ terms

波数k定义为：

k=htr​1​4htr​3H​​[1+ Higher ⋅ order ⋅ terms ]

波速c定义为：

c=ghtr​​(1+hH​[21​−K(m)E(m)​]+ Higher ⋅ order ⋅ terms )

波频率ω定义为：

ω=kc

浅波的波轮廓定义为：

ζ(X,t)=−h+htr​[1+htr​H​cn2(α,m)+ Higher ⋅ order ⋅ terms ]

式中，α定义为：

α=kx​(x−x0​)+ky​y−ωe​t

式中，x0​为参考坐标系原点在x^方向的平移距离。

五阶波理论的速度分量表示为：

u=​c−ghtr​​⎩⎨⎧​1−i=1∑5​(34(khtr​)2​)ij=0∑i−1​(htr​z+h​)2jl=0∑i​cn2l(α,m)Cijl​⎭⎬⎫​​cos(θ)

v=​c−ghtr​​⎩⎨⎧​1−i=1∑5​(34(khtr​)2​)ij=0∑i−1​(htr​z+h​)2jl=0∑i​cn2l(α,m)Cijl​⎭⎬⎫​​sin(θ)

w=​amp;ghtr​​2khtr​cn(α,m)sn(α,m)dn(α,m)amp;​i=1∑5​(34(khtr​)2​)ij=0∑i−1​(htr​z+h​)2j+1l=1∑i​cn2(l−1)(α,m)2j+11​Cijl​​​

式中Cijl​为文献[166]中提到的数值系数。

孤波（Solitary wave）理论表达式是通过假设波具有无限波长而导出的。基于这个假设可以得到：

m=1

htr​=h

14.3.11.4 波理论选择

将线性与非线性波理论应用于浅波受到 Ursell 数的限制。对于非线性波浪，Ursell 数准则基于波浪是单峰的，在波谷处不产生任何次波峰。应根据波浪的陡度和相对高度来选择波浪理论，因为波浪会随着波浪的增加试图获得非线性模式波陡度或相对高度。二阶和四阶波理论更容易产生二次波峰。为特定应用选择正确的波浪理论取决于破波和稳定性限制内的输入参数，如下所述。

1、在破波极限内强制检查完整的波浪状态

破波极限内的波高与液深之比H/h（相对高度）定义为：

[hH​]max , theoretical ​=0.78

[hH​]max , Numerical ​=0.55

破波极限内的波高与波长之比H/λ（波陡）定义为：

[λH​]max,theoretical​=0.142

[λH​]max,Numerical​=0.12

2、检查波浪破碎和稳定极限内的波浪理论

对于此检查，波浪类型以下列方式用适当的波浪理论表示：

线性波：Airy波理论

Stokes波：五阶Stokes波理论

浅水波：五阶Cnoidal/孤波理论

在线性和五阶斯托克斯波之间执行二阶到四阶之间的波理论检查。

波型检查

短重力波的水深与波长之比H/λ定义为：

[λh​]min​=0.5

对于Stokes波，该比值定义为：

[λh​]min​=0.06

对于浅水波（shallow wave），比值定义为：

[λh​]max​=0.085

波陡检查

对于线性波的波高与波长的比值H/λ定义为：

[λH​]max​=0.02tanh(2πλh​)

对于Stokes波，该比值定义为：

[λH​]max​=0.142tanh(2πλh​)

注：斯托克斯波在波陡度低于 0.1 时通常是稳定的。波在 0.1 到 0.12 的范围内可能是稳定的或不稳定的。波最有可能在波陡度值高于 0.12 时破裂。

对于浅水波，该比值定义为：

[λH​]max​=0.142(2πλh​)

相对高度检查

对于线性波，波高与水深之比H/h定义为：

[hH​]max​=0.1

对于浅水区域的Stokes波，该比值定义为：

[hH​]max​=0.5

注：斯托克斯波在相对高度低于 0.4 时通常是稳定的。波浪最有可能在相对高度高于 0.4 时破裂。

对于浅水波，比值定义为：

[hH​]max​=0.55

3、Ursell数稳定性准则

Ursell数定义为：

Ur=h3Hλ2​

线性波的 Ursell 数稳定性准则定义为：

(Ur)max​=332π2​

浅水区的Stokes波的判据定义为：

(Ur)max​=38π2​

注：斯托克斯波对于Ur<10来说通常是稳定的。 10 ~26 之间的 Ursell 数代表向浅层的过渡，因此斯托克斯波可能是不稳定的。斯托克斯波更适用于中深水深区域。

浅水波的判据定义为：

(Ur)min​=38π2​

14.3.11.5 波的叠加

波的叠加原理可以应用于同时通过同一介质的两个或多个波。合成的波型是由单个波的波剖面和速度势的总和产生的。根据叠加的性质，可能出现以下现象：

干扰：相同频率和波高的波的叠加，沿相同方向传播。相长干涉发生在波以相同相位传播的情况下，而相消干涉发生在波以相反相位传播的情况下。

Beats/Wave Group：相同波高的波以几乎相同的频率沿相同方向传播的叠加。节拍是通过随时间增加/减少波的合成频率而产生的。对于海洋应用，这种波浪叠加的概念被称为波浪群。在波群中，波包（波群）以比包内每个单独波的相速度慢的群速度移动。

驻波：沿相反方向传播的相同波的叠加。这种叠加形成了一个驻波，其中每个粒子执行具有相同频率但不同幅度的简谐运动。

不规则波：具有不同波高、频率、相位差和传播方向的波的叠加。不规则海浪在海洋应用中很常见，可以通过叠加多个波浪来建模。

注：叠加原理只对线性波和小振幅波有效。有限振幅波的叠加可能是近似的或无效的，因为它们的非线性行为和波色散对波高的依赖性。

14.3.11.6 使用波谱模拟随机波

海面上的随机波浪是由风在自由表面上的作用产生的。风速越高，风的持续时间越长，风吹（取）的距离越长，产生的波浪就越大。因此可以通过指定一个描述波频率范围内能量分布的波谱来模拟明渠波边界条件的随机波作用。 Fluent 提供了几种基于风速、方向、持续时间和获取的波谱公式。

14.3.11.6.1 定义

1、完全发展的海况

当吹过海面的风将最大的能量传递给海浪时，可以认为海洋已经完全发展。这意味着海况与风吹（取）的距离和风的持续时间无关。在这种情况下，可以假设海平面在统计上是稳定的。海况通常由以下参数表征，这些参数可由风速估计：

有效波高（Hs）

最大 1/3 波浪的平均波高。

峰值波频率 (ωp) 或周期

对应于最高波能量的频率或周期。

2、Long-Crested Sea

如果观测到的波浪的不规则性仅在主导风向，以单向波为主，幅度不同但相互平行，则称该海为长顶不规则海。

3、Short-Crested Sea

当沿多个方向的波峰出现明显的不规则时，海被称为短波峰。

14.3.11.6.2 波谱实现理论

14.3.11.6.2.1 长波峰随机波（单向）

长波峰随机波谱是单向的，只是频率的函数Ansys Fluent 提供以下波谱公式：

Pierson-Moskowitz Spectrum

JONSWAP Spectrum

TMA Spectrum

1、Pierson-Moskowitz Spectrum

Pierson-Moskowitz 谱适用于完全发展的海洋，并假设波浪与风在无限制的获取 [2] (p. 1009) 中达到平衡。

SPM​(ω)=165​ω5HS2​ωp4​​e−(4ω45ωp4​​)

式中，ω为波频率；ωp​为峰值波频率；Hs​为有效波高。

2、JONSWAP Spectrum

JONSWAP 谱是 Pierson-Moskowitz 谱的获取限制版本，其中波谱被认为由于非线性波波相互作用在很长的时间和距离内继续发展。与 Pierson-Moskowitz 谱 [2] (p. 1009) 相比，波谱中的峰值更为明显。

SJS​(ω)=SPM​(ω)Aγ​γ∧exp(−0.5(σωp​ω−ωp​​)2)

式中，SPM​(ω)为Pierson-Moskowitz谱；ω为波频率；ωp​为峰值波频率；γ为峰形参数；Aγ​=1−0.287lnγ；σ为谱宽参数，当ω≤ωp​时，σ=0.07，当ω>ωp​时，σ=0.09。

3、TMA Spectrum

TMA 谱是 JONSWAP 波谱的修正版本，适用于有限水深 [2] (p. 1009)。

STMA​(ω)=SJS​(ω)ϕ(ω,h)

式中，ϕ(ω,h)为指定的深度函数：

ϕ(ω,h)=(sinh(kh))2+gω2h​(cosh(kh))2​

式中，k为通过耗散关系计算得到的波数：

ω2​=gktanh(kh)

式中，g为重力加速度；h为水深。

14.3.11.6.2.2 短波峰随机波（多向）

短波峰是多向的，是频率和方向的函数。短波峰频谱表示为：

S(ω,θ)=S(ω)G(θ,ω)

式中，S(ω)为频谱；G(θ,ω)为方向扩散函数。

需要满足以下等式：

∫θmin​θmax​​S(ω,θ)dθ=S(ω)

对扩展函数施加了以下条件：

∫θmin​θmax​​G(θ,ω)dθ=1

式中

（未完）

14.3.11.6.3 选择波谱和输入

波谱由具有两个基本输入特征的波能量分布函数表示：有效波高 ( ) 和峰值角频率 ( )。此外，您将输入最小和最大波频率 ( 和 ) 以及用于评估频谱的频率数。您应该选择频率范围，以使大部分波能量包含在指定范围内。一般来说，推荐的频率范围是：

最小角频率ωmin​=0.5ωp​

最大角频率ωmax​=2.5ωp​

在选择频谱和输入值时，重要的是要分析和考虑生成波的各种参数。 Fluent 可以根据输入参数计算各种指标，并评估所选频谱和参数是否合适。计算的指标详述如下。

1、波长估计

波长通过下式计算：

λ=k2π​

式中k为波数，通过下面的关系计算：

对于中等水深，ω2=gktanh(kh)

对流短重力波，ω2=gk

因此有：

最小波长：λmin​=kmin​2π​

最大波长：λmax​=kmax​2π​

峰值波长：λp​=kp​2π​

2、时间周期估计

峰值周期由下式计算得到：

Tp​=ωp​2π​

平均时间周期和过零时间周期近似如下：

平均周期：Tm​=0.7303+0.04936γ−0.006556γ2+0.000361γ3

过零周期：Tm​=0.7303+0.04936γ−0.006556γ2+0.000361γ3

式中，γ为峰形参数。

3、波型检查

基于短重力波假设，应当确保：

λmax​h​>0.5

一般而言，TMA 谱适用于中等深度，而 JONSWAP 和 Pierson-Moskowitz 谱在短重力波假设下有效。

4、相对高度检查

要使中等深度/深水处理有效，相对高度应满足：

hHs​​<0.25

5、海陡度检查

海陡度的计算方法有两种：

基于峰值周期：Sp​=gTp2​2πHs​​

基于过零时间周期：SZ​=gTz2​2πHs​​

与每种方法的限制关系进行比较：

SZ​≤⎩⎨⎧​101​101​612−Tz​​+151​6Tz​−6​151​​amp;TZ​≤6samp;6s

Sp​≤⎩⎨⎧​151​151​715−Tp​​+251​7Tp​−8​251​​amp;Tp​≤8samp;8s

6、峰形参数估计

峰形参数γ通过以下方式估计：

γ=⎩⎨⎧​55exp(5.75−1.15Hs​​Tp​​)1​amp;Hs​​Tp​​≤3.6amp;3.6

注：γ=1时 ，JONSWAP 频谱简化为 Pierson-Moskowitz 频谱。

14.3.11.7 明渠波变量名

（略过）

14.3.12 耦合 Level-Set 和 VOF 模型

水平集方法是计算具有拓扑复杂界面的两相流的一种常用的界面跟踪方法。这类似于 VOF 模型的界面跟踪方法。在 Level-Set 方法 479 中，界面被 Level-Set 函数捕获和跟踪，该函数定义为与界面的有符号距离。由于水平集函数是光滑连续的，因此可以准确地计算其空间梯度。这反过来将产生界面曲率和由曲率引起的表面张力力的准确估计。然而，水平集方法在保持体积守恒方面存在缺陷 467。另一方面，VOF 方法自然是体积守恒的，因为它计算和跟踪每个网格中特定相位的体积分数，而不是界面本身。 VOF 方法的缺点在于计算其空间导数，因为 VOF 函数（特定相的体积分数）在界面上是不连续的。针对 Level-Set 法和 VOF 法的不足，在 ANSYS FLUENT 中提出了 Level-Set 法和 VOF 法的耦合方法。

注：水平集和 VOF 耦合模型是专门为不涉及传质的两相流动而设计的，仅适用于瞬态流动问题。二维情况下，网格应限制为四边形、三角形或两者的组合；三维情况下，网格应限制为六面体、四面体或两者的组合。

14.3.12.1 理论

水平集函数被定义为到界面的带有符号的距离。因此，在两相流系统中，界面为零水平集，即 φ(x,t)，可表示为 Γ=x∣φ(x,t)=0 ：

φ(x,t)=⎩⎨⎧​+∣d∣0−∣d∣​ifx∈主相ifx∈Γifx∈次相​(14.109)

式中，d 为与界面的距离。

水平集函数的演变可以用类似于 VOF 模型的方式给出：

∂t∂φ​+∇⋅(uφ)=0(14.110)

式中，u 为底层速度场。

动量方程可以写成：

∂t∂(ρu)​+∇⋅(ρuu)=−∇p+∇⋅μ[∇u+(∇u)T]−Fsf​+ρg​(14.111)

14.3.12.1.1 表面张力

在方程 (14.111)中，Fsf​ 为由表面张力效应产生的力，其计算公式为

Fsf​=σκδ(φ)n(14.112)

式中，σ 为表面张力系数；κ 为局部平均界面曲率；n 为局部界面法向向量。

且有：

δ(φ)={02α1+cos(πφ/α)​​∣φ∣≥α∣φ∣<α​(14.113)

式中 α 为界面厚度；h 为网格间距。

法向向量 n 及界面曲率 κ 可以使用下式进行估算：

n=∣∇φ∣∇φ​​φ=0​(14.114)

κ=∇⋅∣∇φ∣∇φ​​φ=0​(14.114)

在某些情况下，应用公式 14.112 中的默认表面张力会导致溶液中出现假流动。为了减轻这些影响，Fluent 提供了两个加权函数，可以在界面单元中将表面张力重新分配给较重的相。

1、密度修正

在密度校正公式中，公式 14.112 通过引入密度比进行修正：

Fsf​=0.5(ρ1​+ρ2​)ρ​σκδ(φ)n(14.116)

式中，ρ 为基于体积的密度。

2、Heaviside 函数缩放

在 Heaviside 函数缩放公式中，公式 14.112 通过引入 Heaviside 函数进行了修正：

Fsf​=2Hφ​σκδ(φ)n(14.117)

其中：

Hφ​=⎩⎨⎧​0121​[1+aφ​+π1​sin(aπφ​)]​∣φ∣>a∣φ∣>a∣φ∣≤a​ (gas phase) (liquid phase) ​(14.118)

14.3.12.1.2 通过几何方法重新初始化水平集函数

根据水平集函数方程 14.110（第 625 页）的传输方程的性质，∣∇φ∣=1 的距离约束在求解后不可能保持不变。无法保持的原因是界面的变形、不均匀的轮廓以及整个界面的厚度。这些误差会在迭代过程中不断累积，导致质量和动量求解产生较大误差。因此，每个时间步都需要重新初始化。这里使用的是几何界面前沿构造方法。几何方法涉及一个简单的概念，在生成准确的界面前沿几何数据方面非常可靠。VOF 值和水平集函数值都用于重建界面前沿。也就是说，VOF 模型提供了可能的界面通过的单元中切口的大小，而水平集函数的梯度决定了界面的方向。

在构建界面前沿时，还采用了片断线性界面构建（PLIC）的概念。界面前沿重建的程序可简述如下，如图 14.5 所示

找到交替符号 sign(φ) 或体积分数值介于 0 和 1 之间的界面前沿单元格，即部分填充单元格

根据水平集函数梯度计算每个前单元中界面段的法线

定位切穿单元，确保与相邻单元相比，该单元至少有一角被指定相占据

找到单元中心线法线与界面的交点，以满足 VOF

找到界面线与单元边界的交点；这些交点被称为前沿点。

一旦重建了界面前沿，就可以按以下步骤开始最小化给定点到界面的距离：

计算域中给定点到前单元每个切面的距离。距离计算方法如下：

如果从给定点到界面的法线在切分段内相交，则将计算出的距离作为到界面的距离。

如果交点在切分段端点之外，则将给定点到切分段端点的最短距离作为到界面切分段的距离，如图 14.6 所示

最小化给定点到所有前切面的所有可能距离，从而表示给定点到界面的距离。因此，这些距离的值将用于重新初始化水平集函数。

14.3.12.2 局限性

该模型目前仅适用于两相流动体系，即两种流体不相互渗透。

水平集模型只能在启用 VOF 模型时使用。

耦合 VOF 水平集模型不支持以下功能：

质量传递

周期边界

多面体网格

动网格

重叠网格