VOF 法与Level Set 法的比较
汪迎春1，汪德爟2
（1．2．南京河海大学环境学院，南京，210098）
摘 要：VOF 法和Level Set 法同属于界面捕获类算法，有很多相似之处。本文对两种方法进行了详细的介绍，并用剪切流场进行了测试，比较分析了两种方法的优缺点。测试的结果表明，VOF 法容易破碎，界面重构复杂，但有些VOF 实现方案能够模拟非常尖锐的界面，这一点上要比Level Set 法好。Level Set 法不需要重构界面，模拟出的界面形状要比VOF 法光滑，整体效果要比VOF 法好。 关键词：Level Set；VOF；界面
1 引言
在计算流体力学中，界面分层流问题是很普遍的。既有水和气强分层的自由面问题，也有冷热水弱分层的内界面问题。这类问题的计算因为其复杂性成为当前的热点和难点。目前解决此类问题的数值方法有很多。LevelSet 法与VOF 法同属于界面捕获类方法，因为它们在欧拉网格下构造，都用一个标量函数描述界面的几何特性，具有很强的拓扑处理能力，因而是两种很流行的界面流模拟方法。当代大学生与中华优秀传统文化自信
VOF 法在20世纪70年代末由Hirt 和Nichols 等最先提出。VOF 法是80年代国外最流行的方法，目前国内对VOF 法的研究和应用也越来越多。
高射炮打坦克水平集方法（Level Set Method）由Osher和Sethian [1]于1988年提出后，引起了人们的广泛的注意。水平集方法不需要显式地追踪运动界面，避免了对拓扑结构变化的处理，计算稳定，已经应用到诸多领域，如材料学、流体力学、图象处理、计算机视觉、网格生成等。Level Set法是90年代国外很盛行的方法，目前国内的研究和应用情况还不如VOF法充分。
2 两种方法介绍
2.1 VOF 方法
VOF 方法的基本思想是，在欧拉网格系统上定义一个函数 f(流体体积份额量），根据每个网格内所含某种物质的体积量来定义了在此网格上的值，然后用体积跟踪的方法求解方程： 0=∇⋅+∂∂f u f t VOF方法采用几何与代数相结合的方法，即在上一时刻重构界面，然后用代数方法算出下一时刻的 f 值。VOF方法在代数处理阶段常采用贡献网格方法（donor-cell) ，通量修正方法（flux 一 corrected transport , FCT )，迎风方法。界面重构策略可分为分片常数体积跟踪方法和分片线性体积跟踪方法。分片常数体积跟踪方法主要有两种：一种是 Hirt [2]等的VOF法，另一种是 Noh 和Woodward 的 SLIC 方 ———————————
作者简介：汪迎春(1976-)，男，南京河海大学环境学院，博士研究生
通讯地址：南京市西康路1号 河海大学七舍三单元506 E-MAIL：wangyingch@sohu
www.paper.edu
法[3]。SLIC方法中界面重构在同一时刻要求与某一坐标轴平行。而 Hirt方法是根据界面网格内 f 的大小，结合相邻网格 f 的值，选择与之适合的坐标轴平行，相对接近于真实界面。分片常数体积跟踪方法重构界面精度低，易产生非物理界面现象。分片线性体积跟踪法，最早是由Debar提出的，典型的有Youngs方
单因子指数法法[4]和Flair方法。
在分片线性体积跟踪方法中 Youngs方法的界面重构是在一个网格内进行，界面不穿过网格边（面）。在二维流体运动，网格内界面的法向由 9 个网格的 f 的值决定。界面的位置由法向和界面所在网格的 f 的值及相邻网格的 f 的值确定，并且保证体积守恒。界面在每一个网格内为一直线段。Flair方法是用相邻的两个网格上的有斜率的直线逼近界面。以上方法都都采用方向分裂算法向多维扩展，CICSAM [5]方法通过半隐式处理，通过求网格边界通量形式向多维扩展，便于应用于非结构网格。
VOF 还有很多其他格式[6] [7]，不再赘述。
2.2 Level Set 方法
Level Set方法把随时间运动的物质界面看作某个函数φ(x,t) 的零等值面,φ(x,t) 满足一定的方程。在每个时刻t ，我们只要求出函数的值，就可以知道其零等值面的位置，也就是运动界面的位置。构造函数φ(x,t) , 使得在任意时刻运动界面Γ(t) 恰是φ(x ,t) 的零等值面，即Γ(t) ={x ∈Ω: φ(x,t) = 0},φ(x ,t) 的初值应满足在Γ(t)附近为法向单调, 在Γ(t) 上为零。一般可取φ(x ,t) 为x 点到界面Γ(0) 的符号距离。
⎪⎩
⎪⎨⎧Ω∈Γ−Γ∈Ω∈Γ=21
x (0)
x x ))0(,(0))0(,()0,(x d x d x φ 其中d(x,Γ(0))表示x到Γ(0)的距离，Ω= Ω1∪Γ∪Ω2。
在任意时刻t, 对于活动界面Γ(t) 上的任意点x,φ(x,t) = 0，从而有 0=∇⋅+∂∂=φφφV t
甲基丙烯酸dt d ，dt dx V = （1） 方程（1）为Level Set 方程。Level Set 函数的单位外法向以及运动界面的曲率分别为
0=∇∇=φφφ
n ， 0)(=∇∇⋅∇=φφ
φk 根据（1）求得的φ( x,t)并不满足符号距离函数性质，为了保持φ( x,t) 的符号距离性质，在每个计算步后需要重新初始化，通过求解初值问题的稳定解来实现。
⎩
⎨⎧=∇−=00)0,()1)((φφφφφτx sign （2） 式中τ为虚时间。为了方便求解，将符号函数sign(φ) 光滑化为
22)(εφφφε+=sign
重新初化方程的求解常用一阶、单调的Godunov 格式，见后。
Level Set 的实现方法可归纳为两大类： PDE 实现方法和推进扫描类方法。PDE 类实现方法有： ①PDE+重新初始化，②边界元法[8]。推进扫描类方法有：①快速行进法[9]，②源点扫描法[10]。
2． 3两种方法的特性
共性：两种方法都通过求解每个网格上的一个标量函数来描述界面的几何特性，标量函数的输运方程在形式上完全相同，同属于界面捕获类算法，都有很强的拓扑表达能力。
区别：两种标量函数的定义不同，VOF 法定义的f 是相物质的流体体积份额量，它的大小在0和1之间，界面附近以外的网格单元只可能是0或1，界面附近的网格单元在0和1之间。所以在界面附近的网格单元f 是不连续的，存在很强的间断，对于低阶的对间断分辨能力不强的数值格式，f 的求解效果不是很理想，但高阶格式在非结构网格下又不易构造，这两者存在一定的矛盾。VOF 法的界面重构一般比较复杂，不同的重构策略，重构出的界面差别比较大，重构策略直接影响着VOF 法的精度。Level Set 方法定义的φ是格点离开界面的最短距离的符号距离函数，φ沿界面法向单调，根据φ的符号可判断格点在界面的内侧还是外侧。φ在计算域上一般是连续的，并且
1=∇φ，φ呈很有规则的分布，对于连续性的问
www.paper.edu
地域分布题总是有很好的数值方法，这也是Level Set 相对VOF 法的优势所在。φ直接隐含着界面的几何属性，不需要重构界面，描述的界面会比VOF 法光滑。
3 数值格式
本文对Levelset
（1） 二阶迎风格式
单元边界均值e φu u e 式中，下标u 代表迎风单元；为从迎风单元中心至单元面的距离矢量。
透传dr 其中的单元中心梯度由Gauss 定理计算
i j j j u i
A V x ∑Δ=⎟⎟⎠⎞⎜⎜⎝⎛∂∂φαφ 其中：为单元j 面的i 方向面积分量，i j A j φ可直接用线性插值求得，但是对高度扭曲的畸变网格会引
（2） Quick 格式
图3-2-2 迎风单元示意图 则单元边界均值为： )(1)()()1(21U U U U U f φφγγααφφγγαβφφ−+++−++= 定义迎风单元下标为U，相对迎风单元的上下游单元下标分别记为U1、U2，同时定义U U Uf L L 2α=，U U f U L L 22β=，U U U U L L 21γ=
（3） 三阶TVD-Runge-Kutta 格式
⎪⎪⎪⎩
⎪⎪⎪⎨⎧Δ++=Δ++=Δ+=)(323231)(414143)()2()2()0()3()1()1()0()2()0()0()1(φφφφφφφφφφφtL tL tL
（4） 一阶、单调Godunov 格式
结构网格下的一阶、单调Godunov 格式为：
)()(,,,0,1,n j i j i n j i n j i G S φφτφφεΔ−=+
2
2,,0,,0,,0)(x S j i j i j i Δ+=φφφε
[][]
⎪⎪⎩
⎪⎪⎨⎧<−+>−+=+−+−−+−+其他 0 0 1))(,)max(())(,)max(( 0 1))(,)max(())(,)max(()(j i,0,2/12222j i,0,2/12222,,0φφφd c b a d c b a G j i
其中:
a (=max(z =+a ⎩⎨<⋅+=0
)(x x n dr grad e e e p p e φφφ
b 对应的有限体积法单元边界均值插值形式
⎩⎨⎧>⋅+=<⋅+=0)(0)(x x n dr grad n dr grad e e e p p e φφφφφφ c,d根据n y 选择插值形式，类似a,b的处理 4实例比较 本文分别用VOF 法和Level Set 法计算了剪切流场中界面运动的轨迹，剪切流场是检验界面流算法性能的很好的算例，能检验算法对锐利界面的捕获能力以及算法的破碎程度。
速度场的分布为： ⎩⎨⎧−−−=−−=))(cos())(sin(),())(sin())(cos(),(0000y y x x y x v y y x x y x u ππππππ 该流场是不可压缩处处守恒的。取x 0=0.5,y 0=0.5，计算区域为[0，1]×[0，1]，初始界面为圆心在
(0.5，0.3)，半径为0.2的圆周，如图1所示。空间步长为0.01，计算时间为1s。采用单纯拉哥朗日粒子法计算的精确结果如图
4.2。
图4.1 初始界面及流场 图4.2 单纯拉哥朗日粒子法
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议