VOF与Level Set耦合的界面计算方法及其对波浪破碎过程的模拟

贾 伟，柳淑学，李金宣，范玉平

(大连理工大学 海岸和近海工程国家重点实验室，大连 116024)

近年来，在海洋工程流体流动和波浪传播过程的数值模拟中，经常采用两相流模型(Two-Phase Fluid Model)，即将空气和水视为两种流动参数不同，但满足相同控制方程的流体。该模型采用同一组控制方程同时求解两相流体，并采用界面计算方法来捕捉两相界面，因此也被称为拟单流体法。界面计算方法在两相流模拟中至关重要，VOF方法和Level Set方法是目前两相流模型中较为常用的两种界面处理方法。根据计算方法不同，VOF方法又分为代数方法和几何方法。

最为常见的VOF代数方法是界面压缩法，例如开源代码OpenFOAM中采用的MULES(Multidimensional Universal Limiter with Explicit Solution)方法[1]。该方法通过限制器来混合一阶与高阶格式，并且在对流方程中添加压缩项来保持界面的尖锐性。代数方法精确度相对较低，但是由于没有几何计算以及界面重构步骤，因此计算效率高。几何方法计算较为精确，但是计算量大，并且往往对网格类型有所限制，其中最常见的是由Youngs[2]提出的PLIC(Piecewise Linear Interface Construction)方法。

VOF方法的优点是可以保持质量守恒，但由于采用了体积分数来表示流体在体积单元的占比，因此体积分数场不连续，无法精确计算界面处的曲率和法向量，也无法精确计算界面处的表面张力。Level Set方法采用Level Set函数表示某点与界面的距离，曲率和法向量的计算更为准确。但该方法质量不守恒，当网格较稀疏时尤为明显。

因此，一些学者将VOF方法与Level Set方法进行耦合，从而利用两种方法各自的优点并规避缺点。由于侧重点和计算方法各不相同，各耦合方法也有较大差异。Sussman等[3]采用二阶算子分裂算法，已知某时刻的体积分数、Level Set函数和速度场，根据该时刻各个单元体积面的通量求解下一时刻的体积分数场和Level Set函数。相比于Level Set方法，对于界面有尖角时，耦合方法计算效果好。Pijl等[4]通过体积分数与Level Set函数之间的代数关系，提出了二者之间简化的耦合方法，从而避免了在每个时间步内同时计算体积分数和Level Set函数的对流方程。Sun等[5]提出了耦合的VOSET方法，该方法采用PLIC方法计算体积分数和重构界面，并提出了一种迭代的几何计算方法，用于计算界面附近的Level Set函数，从而计算界面的曲率。Albadawi等[6]采用Level Set函数计算表面张力和界面曲率，通过耦合方法模拟了液体中的气泡，并与实验结果进行了对比。该方法适用于计算常接触角问题，主要对表面张力进行了修正。Cao等[7]采用多维对流算法求解体积分数，采用类似Sun等的方法几何迭代求解Level Set函数并进行耦合，该方法可用于非结构化网格。

在现有的数值波浪计算模型中，通常采用单一的VOF方法或Level Set方法计算界面，即存在无法精确计算界面法向量或无法保持质量守恒的缺点。此外，现有的大部分VOF与Level Set耦合方法不适用于非结构化网格，并且计算量较大。因此，本文建立了基于VOF与Level Set方法的耦合方法，在保证计算精度的同时提高计算效率。该方法可用于结构化网格及非结构化网格，并易于扩展到三维。随后，采用该耦合方法计算了溃坝问题，验证了方法的有效性，并进一步应用于波浪传播过程中破碎问题，与实验结果进行了对比，表明该方法可以较好地应用于波浪破碎过程的模拟。

1 数值方法

1.1 两相流模型与界面处理方法

本文的数值模型基于OpenFOAM中的不可压缩两相流模型求解器interFoam进行开发，其控制方程为连续性方程以及N-S方程

·U=0

(1)

(2)

在interFoam求解器中采用VOF方法计算界面，即采用体积分数α作为控制参数，两相交界处的流体性质参数可以表示为

ρ=αρw+(1-α)ρa

(3)

μ=αμw+(1-α)μa

(4)

式中：ρw、ρa分别为液体、气体的密度，μw、μa分别为液体、气体的动力粘滞系数。体积分数α满足相分数方程

(5)

界面处的法向量可以表示为

(6)

在本文模型中，与Level Set方法进行了耦合。Level Set方法采用Level Set函数φ(x,t)表示t时刻x位置与界面的距离，在界面处函数值为0。用Γ(t)表示界面位置，则0时刻φ(x,0)的表达式为

(7)

其中d是与界面的距离。此外，在Level Set方法中，通常采用Heaviside函数光滑界面附近的密度、动力粘滞系数等不连续物理量，即

ρ(φ)=H(φ)ρw+(1-H(φ))ρa

(8)

μ(φ)=H(φ)μw+(1-H(φ))μa

(9)

其中H(φ)为Heaviside函数

(10)

式中：φ为Level Set函数，b为界面厚度，通常取b=1.5△x，其中△x为网格尺寸。

与体积分数相似，Level Set函数需要满足对流方程

(11)

需要注意的是，通常在采用Level Set方法计算界面时，由式(11)计算所得的函数值并不能保证与该时刻点x和界面的距离相等，因此在每个时间步前需要对Level Set函数重新初始化[3]。本文采用了几何迭代方法来计算Level Set函数，避免了求解对流方程及初始化过程。

由于Level Set函数是一个连续函数，因此可以较为精确地计算各点的法向量和曲率，即

(12)

(13)

在本文模型中，曲率采用式(13)计算，则控制方程式(2)可以表示为

(14)

其中ρ和μ等不连续物理量采用式(8)和式(9)进行计算。

1.2 耦合算法

本文采用几何方法求解体积分数从而保证计算精度，只在界面附近通过几何关系初始化Level Set函数，并迭代求解，从而避免求解Level Set函数的对流方程以及重新初始化方程，提高计算效率。该算法分为3个部分：求解体积分数的对流方程、界面重构和构造Level Set函数。

1.2.1 求解体积分数的对流方程

本文在计算体积分数的对流方程时，采用了isoAdvector方法，详细内容可以参考Roenby等人的论文[8]，其开源代码已植入到OpenFOAM-v1612+及之后的版本中。该方法的简要步骤如下：

(1)寻找界面单元，即体积分数介于0与1之间的单元。

(2)对于界面单元的顶点，用共享该顶点的所有网格单元内的体积分数进行加权平均，权重系数为该节点到单元中心距离的倒数。该步骤可以将体积分数插值到顶点，假设某界面单元有N个顶点，其坐标分别表示为x1，x2，…，xN，相应的体积分数插值结果为f1，f2，…，fN。

(3)根据各顶点的f值，可以在单元内部通过线性插值构建关于f的等值线(二维问题)或等值面(三维问题)。如果该单元的某个边(xm，xn)满足fm

(15)

找到每条边的交点并连接，即可得到等值面。

(4)找到等值面，使得该等值面按照该单元的体积分数αi将单元分成2个部分，即求f*，使α(f*)=αi。

首先根据几何关系计算f1，f2，…，fN对应的体积分数值α(f1)，α(f2)，…，α(fN)。根据该单元的体积分数αi，在α(f1)，α(f2)……α(fN)中寻找最接近αi的两个值α(fm)和α(fn)，使得α(fm)≤αi≤α(fn)。在此区间内，α与f之间的关系为三次多项式，可以在此区间内通过线性插值构造两个等值面，由几何关系求得相应的α，即可求解4×4矩阵得出多项式的系数，进而确定f*的值，并得到目标等值面。

(5)根据等值面在当前时间步的运动，即可通过几何方法求出该单元下一时间步的体积分数。

1.2.2 界面重构

图1 某网格单元内的等值线示意图

需要注意的是，上述isoAdvector方法构建的等值面仅用于计算体积分数，并非真实的两相界面。由于在构造Level Set函数时需要用到界面的几何信息，而isoAdvector并未包含界面重构的步骤，因此在本模型中建立了简便的重构算法：

(1)按照式(6)，根据体积分数场计算精度较低的法向量n′。

(2)对于顶点坐标为x1，x2，…，xN的界面单元，令d′ =n′·x，即可在单元内构建关于d′的等值面，且等值面均与法向量n′垂直，对于二维网格如图2所示。

(3)此时，可以通过1.2.1(4)中的方法计算关于d′的等值面，使得α(d′)=αi。由于该界面与法向量垂直，且将该单元划分为体积分数为αi与1-αi两个部分，因此该等值面即该单元内的界面。

图3 算法流程图

1.2.3 构造Level Set函数

参考Sun 等[5]，Cao等[7]，Level Set与VOF的耦合过程采用了几何迭代求解的方法。其简要步骤如下：

(1)标记界面单元附近3△x区域内的体积单元，其中△x为网格尺寸。

(2)计算被标记区域内的Level Set函数。通过几何关系，计算网格点与1.2.2节所得重构界面的距离d，即可按照式(7)构造Level Set函数。

(3)按照式(12)计算较为精确的法向量，从而按照1.2.2节的步骤重新构造更为精确的界面。构造Level Set函数和重构界面的过程可以迭代两次，从而得到最终的界面。

整体而言，上述方法可以视为一种改进的VOSET方法。在计算体积分数时，原始的VOSET方法采用了PLIC方法，该方法适用于结构化网格，拓展到非结构化网格时较为复杂，因此本文采用了isoAdvector算法，该算法计算速度较快，并且适用于各类网格；之后参考isoAdvector中计算等值面的方法，建立了重构界面的简便算法；采用VOSET方法中的迭代几何方法计算Level Set函数，从而精确计算界面处的法向量及曲率，避免了求解Level Set函数的对流方程及初始化方程。该算法的计算流程如图3所示。

2 数值模拟验证

图4 溃坝算例的布置图

Fig.4 Layout of dam break problem

图5 水体前端的无因次位置与无因次时间的关系

Fig.5 Relationship between dimensionless locations of water front and dimensionless time

溃坝问题是经典的两相流算例，因此通过模拟该问题，与Martin 等[9]的实验结果进行对比。图4为溃坝算例的布置图，整个计算域是边长为4L的正方形，其中L=0.146 m。在计算域左下方有一块宽为L，高为2L的静止水体，当计算开始后，水体由于重力作用开始流动。水的密度为1.0×103kg/m3，表面张力系数为0.075 5 N/m，重力加速度为9.8 m/s2。

在数值模拟时，采用均匀的正方形网格，共100×100个。采用图4所示的xOy坐标系，将运动水体的最前端与y轴的距离用x表示。将x和时间t无因次化，即令X*=x/L,t*=t·(2g/L)0.5，则计算所得X*与t*的关系如图5所示。图中同时给出了Martin 等[9]的实验结果以及Sun 等[5]和Cao等[7]的VOSET模型计算结果。从图中可以看出，本文模型的计算结果与实验结果及文献中VOSET模型的计算结果均较为接近，因此可以说明本文的界面计算方法是有效的。作为示例，图6给出了溃坝过程中，典型时刻的瞬时水体自由表面图。

6-at=0.0 s 6-bt=0.075 s 6-ct=0.15 s 6-dt=0.225 s 6-et=0.30 s 6-ft=0.375 s

图6 不同时刻的水体自由表面图

Fig.6 Snapshots of the water interfaces at different times

3 波浪传播破碎问题的模拟

波浪破碎是波浪传播过程中的重要现象，由于波浪破碎过程的复杂性和随机性，其界面往往较为复杂。本文模拟了规则波在斜坡地形及岛礁地形上的破碎问题，并与实验结果进行了对比。

3.1 波浪在斜坡上的破碎

图7 波浪在斜坡上破碎实验布置图(单位：m)

该实验在大连理工大学海岸和近海工程国家重点实验室的三维水池中进行[10]，水池长度40 m，宽度24 m，实验水深0.4 m，地形为1:15的斜坡，坡脚与造波机的距离为10 m，坡顶高度0.4 m，从而保证波浪破碎发生在斜坡上，浪高仪及地形布置如图7，波浪参数如表1所示。由于数值计算中仅采用了正向规则波，因此简化为二维问题进行计算。

在数值波浪水槽中，计算域总长度为20 m，高度0.57 m，在水槽边界采用速度入口法[11]进行造波。整个计算域内采用尺寸相同的矩形网格，网格尺寸为△x=0.01 m，△y=0.007 5 m，在坡面与网格交界处应用OpenFOAM中的snappyHexMesh方法，将网格切割从而使网格与坡面贴合。湍流采用k-ωSST模型进行模拟。表1为模拟波浪的参数。

表1 斜坡地形波浪参数

图8 算例1的波面历时曲线

针对算例1所模拟的波面历时曲线与试验测量结果的对比如图8所示，图中右上角标注数字是图7中的浪高仪编号，图中虚线为实验结果，实线为数值计算结果。1号浪高仪位于平底区域，可以看出数值模型中的入射波浪与实验入射波浪基本一致。由4号浪高仪到13号浪高仪，波浪受浅化作用影响，波浪的非线性及波面的不对称性逐渐增强，13号浪高仪处为波浪的破碎点，之后进入破碎区，波高迅速减小。由实验观察可知，该波浪为卷破波，波面翻卷之后拍打坡面，造成液滴飞溅，因此在破碎区内，波浪的紊流特性较为明显，波面随机性极强，同时破碎波浪存在的浪花及水气混合使得波浪测量结果存在一定的误差。整体而言，该数值模型能够较为准确地模拟波浪浅化过程，也能够较为准确地预测波浪破碎位置及破碎波高。

图9显示了算例1在不同时刻的水面变化情况，图中灰色区域为空气，黑色区域为水体。根据实验观察，13号浪高仪为波浪的破碎点(即图中的垂线所在位置，距离造波机14.25 m处)，即卷破波传播到该点时，波前近乎为垂直面。由图示模拟结果可见，卷破波的波前在t=29.9 s时几乎垂直，而在t=29.92 s时已微微前卷，说明数值模型所预测的破碎点与实验结果一致。此外，图9展示了卷破波的波前翻卷，与水体接触之后飞溅(splash-up)的过程，说明该界面计算方法能够较为准确地计算波浪破碎过程。

9-at=29.82 s 9-bt=29.9 s 9-ct=29.92 s

9-dt=30 s 9-et=30.08 s 9-ft=30.16 s

注：图中垂线的横坐标为实验观察所得的破碎点横坐标。

图9 模拟算例1的波浪破碎过程

Fig.9 Snapshots of breaking process in case 1

算例2的计算结果如图10所示。其入射波高与第1组基本相同，但相较于第1组波浪周期增大，波长变大且水深不变，因此波浪非线性明显增强。在坡脚附近，由于波浪的非线性的影响，波浪已经存在较为明显的变形。之后波浪在斜坡地形上发生浅化变形，直到13号浪高仪处发生破碎。与第1组波浪类似，本模型较好地计算了波浪的破碎过程，包括破碎位置、破碎波高等与实验结果基本一致。

3.2 波浪在岛礁地形上的破碎

图10 算例2的波面历时曲线

实验在大连理工大学海岸和近海工程国家重点实验室的波流水槽内进行[12]，实验的水槽长度为69 m，宽0.8 m，深1.8 m，地形及浪高仪布置如图11所示。岛礁地形简化为坡度为1:5的斜坡，坡顶高度为0.5 m，水深0.715 m。模拟波浪的周期为1.6 s，其他参数如表2所示。与3.1部分的算例不同，该地形为简化的岛礁地形，波浪在经过斜坡地形之后进入顶部的礁坪段，即波浪在浅化变形之后再在平面上传播，地形对波浪的影响更为复杂。选取的2组算例均在波浪传播到礁坪之后发生破碎，分别为崩破波和卷破波。

表2 岛礁地形波浪实验参数

图11 波浪在简化岛礁地形上破碎实验布置图(单位：m)

图12 算例3的波面历时曲线

在数值波浪水槽中，计算域总长度为45 m，高度0.9 m，在水槽边界采用速度入口法进行造波。整个计算域内的初始网格均为尺寸相同的正方形，△x=△y=0.012 m，在坡面与网格交界处采用SnappyHexMesh方法，将网格切割从而使网格与坡面贴合。

算例3的数值波面历时曲线与实验结果的对比如图12所示。可以看出，数值模拟结果与实验结果吻合较好。在图中1号、3号浪高仪处，波浪在平底上传播，波形较为稳定；在6号、8号、10号浪高仪处，波浪在岛礁斜坡上传播，由于水深快速变浅，波浪浅化作用明显，波高逐渐增大，波浪的非线性、波面的不对称性逐渐增强；12号浪高仪处于坡顶的礁坪段，波高继续增大，并最终在13号浪高仪处发生波浪破碎；14号、15号、16号浪高仪位于波浪破碎点之后，该区域内湍流极强，波面的随机性也较强，在该区域内数值结果与实验结果存在一些误差，但整体上基本一致。

图13 算例4的波面历时曲线

算例4的计算结果如图13所示。可以看出，与之前组次类似，数值结果与实验结果整体上吻合较好，可以较为准确地计算波浪的浅化作用(9号)、破碎点(12号)的位置以及破碎波高。总之，根据以上4个算例，采用本文的界面计算方法可以较为准确地模拟在复杂斜坡上的波浪传播破碎问题。

4 结论

本文建立了一种VOF方法与Level Set方法耦合的界面计算方法，其中体积分数采用几何方法——isoAdvector方法进行计算，可保证计算精度；建立了简化的界面重构算法，从而便于采用几何方法建立Level Set函数；通过迭代方法计算更为准确的Level Set函数及法向量。该算法是VOSET方法的一种改进，保持VOF方法及Level Set方法各自的优点，并尽可能地保证算法的简洁以及非结构化网格的适用性。此外，虽然本文仅考虑了二维模型，但该方法的各个步骤均易于推广到三维。

基于溃坝问题的求解，证明了该计算方法的有效性。并进一步应用于波浪在斜坡地形、岛礁地形上的传播破碎问题，数值模拟结果与实验结果吻合较好，说明该方法可以很好地模拟复杂地形上的波浪传播破碎问题。