一种由VOF函数快速构造符号距离函数的改进CLSVOF方法与流程

本发明涉及一种由VOF函数快速构造符号距离函数的改进CLSVOF方法，其属于运动界面追踪数值模拟领域。

背景技术：
：

运动界面追踪数值模拟是当前国内外热门研究课题，它涉及到材料、物理、化学、生物等多个研究领域，对其开展研究具有重要的理论价值和现实意义。

上世纪80年代初，美国洛斯阿拉莫斯国家实验室的Hirt和Nichols首先提出了流体体积法(Volume of Fluids，VOF)。该方法在流场中引入一个流体体积(VOF)函数(其值为0～1)，通过求解VOF函数的运输方程，实现运动界面的追踪。VOF法计算量小，界面锐利程度高，同时具有良好的质量守恒性，因此被大量商用软件及国内外学者采用，用于求解复杂的运动界面问题。遗憾的是，由于VOF函数在界面两侧存在间断，无法准确计算出界面的法向量，降低了该方法的精度。

1988年，加州大学伯克利分校的Osher和Sethian基于Hamilton-Jacobi方程的思想提出了著名的水平集法(Level Set Method)，此方法定义了一个连续光滑的Level Set函数(如符号距离函数，Signed Distance Functions，SDF)，通过求解Level Set函数的运输方程，实现运动界面追踪。Level Set法可以很方便地计算复杂界面变化如破碎、合并等；并且由于符号距离函数的连续特性，结合高分辨率类方法如WENO、TVD等格式，可以很准确地计算出界面的法向量。然而，由于数值扩散和重新初始化等原因，Level Set法在计算过程中容易产生质量不守恒。

VOF法的缺点在于VOF函数在界面两侧存在间断，无法准确计算出界面的法向量；而Level Set法的优点则是Level Set函数连续，可以准确计算出界面法向量；如果在VOF法中引入符号距离函数，用其构造界面法向量，那么必然可以有效弥补VOF法的缺点。

2000年，加州大学戴维斯分校Sussamn和Puckett等人首先将VOF法和Level Set法结合，提出了一种耦合的VOF和Level Set方法(Coupled Level Set and VOF Method，CLSVOF)。该方法采用PLIC算法求解VOF方程，再由PLIC算法建立的网格内部线性界面(二维为直线，三维为平面)构造符号距离函数。由于三维线性界面的建立过程很复杂，此方法在三维情况下难以实现，并且很难扩展到其它VOF求解方法。

2010年，西安交通大学的孙东亮等人提出了一种VOSET方法。该方法采用几何方法计算各网格中心到PLIC线性界面界面的距离，再根据VOF函数的大小，将其转换为符号距离函数。与CLSVOF法类似，此方法同样不易扩展到三维。

2013年，英国卡迪夫大学的Yoki等人提出了一种新的CLSVOF方法，该方法采用Xiao等人提出的THINC算法求解VOF传输方程。由于THINC算法中并无几何界面的建立过程，所以需要通过其他方法构造符号距离函数。为此，Yoki等人采用线性插值的方式构造出VOF函数的0.5等值面，再通过Sethian等人提出的快速推进法(FastMarching Method，FMM)及Sussman等人提出的稳态解法重新初始化(Re-initialization)构造符号距离函数，简化了符号距离函数的重构流程。然而，FMM算法逻辑过于复杂，且在构造界面网格处的符号距离函数时存在一定误差，而稳态解法的迭代步数较多，计算效率过低。

技术实现要素：
：

本发明针对CLSVOF方法中存在的符号距离函数构造过程过于复杂的问题，提出一种由VOF函数快速构造符号距离函数的改进CLSVOF方法。

本发明采用如下技术方案：一种由VOF函数快速构造符号距离函数的改进CLSVOF方法，包括如下步骤：

第一步：定义初始VOF函数

根据实际的物理界面，在网格内定义VOF函数，界面网格的VOF函数值为0～1，界面内的网格的VOF函数值为1，界面外的网格的VOF函数值为0；

第二步：获取界面网格

首先通过网格的VOF函数值，初始化符号距离函数，VOF值大于0.5的网格的符号距离函数为正，VOF值小于0.5的网格的符号距离函数为负；其次，采用子网格技术判断当前网格是否为界面网格，所述子网格技术即通过判断当前网格与周围网格的符号距离函数值是否异号来确定其是否为界面网格，若异号，说明运动界面在两网格之间， 网格为界面网格，反之，则为非界面网格；

第三步：将VOF函数转换为符号距离函数

由于VOF函数和符号距离函数代表了相同的界面，那么在界面网格处的符号距离函数Heaviside值应等同于VOF函数值，采用Newton迭代法对上述等同关系进行求解，可获得界面网格的符号距离函数，即符号距离函数值；

第四步：重构符号距离函数

采用快速扫描法对符号距离函数进行重构，所述快速扫描法是一种求解Eikonal方程的快速解法，对符号距离函数重构等同于求解Eikonal方程|▽Φ|＝1；

第五步：采用WENO格式获取符号距离函数梯度

符号距离函数重构后，采用高精度的WENO格式获取梯度；

第六步：用符号距离函数的梯度定义界面法向量

定义界面法向量时，考虑迎风特性，若流场速度为负，取符号距离函数的右导数为界面法向量，若流场速度为正，则取符号距离函数的左导数为界面法向量；

第七步：求解VOF传输方程。

进一步地，所述第四步中包括如下步骤

1、初始化；根据第三步给出的结果，将界面处的网格状态设为“0”，非界面网格状态设为“1”，若非界面网格的符号距离函数值为正，则将其初始化为极大值，否则初始化为极小值；

2、高斯赛德尔迭代；根据快速扫描法给出的格式，对整个计算区域沿不同方向进行扫描迭代；

3、更新符号距离函数；根据初始化给出的网格状态，判断其符号距离函数是否需要更新，若网格状态为“1”，将扫描迭代后的符号距离函数值与原始符号距离函数值进行比较，取最小值，如果网格状态为“0”，则取原始值。

进一步地，所述第五步中包括如下步骤

1、选取插值模板，左导数模板为{I-2、I-1、I、I+1、I+2}、右导数模板为{I+2、I+1、I、I-1、I-2}；

2、构造三个三阶逼近s₁、s₂、s₃，以左导数为例，其插值模板分别为{I-2、I-1、I}、{I-1、I、I+1}、{I、I+1、I+2}；

3、定义三个非线性权ω₁、ω₂、ω₃；

4、定义三个光滑因子β₁、β₂、β₃；

5、构造三个三阶逼近的非线性凸组合。

进一步地，所述第七步中包括如下步骤

1、根据界面速度方向，定义λ和迎风网格iup；

2、计算双曲正切函数的跳跃中心

3、计算网格面流体体积流出量f_i+1/2、f_i-1/2；

4、将一维结果分裂至多维。

本发明具有如下有益效果：

(1)本发明采用的子网格技术、牛顿迭代法、快速扫描法、WENO格式、THINC算法逻辑简单，易于编程实现；

(2)子网格技术和牛顿迭代法和对二维和三维均适用，而快速扫描法也很容易向三维扩展，故本方法可以很容易地应用到三维情况；

(3)子网格技术只需一步计算即可完成，而牛顿迭代法一般和快速扫描法也仅需几步即可收敛，故本方法计算效率极高；

(4)界面网格的符号距离函数由牛顿迭代法得到，具有机器精度，界面外网格的符号距离函数由快速扫描法得到，具有一阶精度，与现有算法相当；

(5)本发明设计的符号距离函数构造方法不依赖于VOF求解方法，它不仅适用于THINC算法，亦可适用于SLIC、PLIC等传统的VOF方法，适用性强；

(6)THINC算法是一种不依赖于时间步、数值稳定性强的算法，本发明完美地继承了这些优点，算法鲁棒性强。

附图说明：

图1为初始VOF函数分布示意图。

图2为界面网格判断示意图。

图3为界面网格的符号距离函数分布示意图。

图4为快速扫描法更新网格符号距离函数示意图。

图5为符号距离函数左右导数的插值模板示意图。

图6(a)和6(b)为一维THINC算法原理示意图。

具体实施方式：

本发明由VOF函数快速构造符号距离函数的改进CLSVOF方法，包括如下步骤：

第一步：定义初始VOF函数

根据实际的物理界面，在网格内定义VOF函数。界面网格的VOF函数值为0～1，界面内的网格的VOF函数值为1，界面外的网格的VOF函数值为0。

第二步：获取界面网格

首先通过网格的VOF函数值，初始化符号距离函数。VOF值大于0.5的网格的符号距离函数为正，VOF值小于0.5的网格的符号距离函数为负。其次，采用子网格技术判断当前网格是否为界面网格。所谓子网格技术，即通过判断当前网格与周围网格的符号距离函数值是否异号来确定其是否为界面网格。若异号，说明运动界面在两网格之间，网格为界面网格；反之，则为非界面网格。

第三步：将VOF函数转换为符号距离函数

由于VOF函数和符号距离函数代表了相同的界面，那么在界面网格处的符号距离函数Heaviside值应等同于VOF函数值。采用Newton迭代法对上述等同关系进行求解，可获得界面网格的符号距离函数，即符号距离函数值。

第四步：重构符号距离函数

为了高效快速重构符号距离函数，本发明采用快速扫描法对符号距离函数进行重构，快速扫描法是一种求解Eikonal方程的快速解法，对符号距离函数重构等同于求解Eikonal方程|▽Φ|＝1。

其实现步骤如下：

1、初始化；根据第三步给出的结果，将界面处的网格状态设为“0”，非界面网格状态设为“1”。若非界面网格的符号距离函数值为正，则将其初始化为极大值；否则初始化为极小值。

2、高斯赛德尔迭代；根据快速扫描法给出的格式，对整个计算区域沿不同方向进行扫描(迭代)。一般来说，二维情况下仅需扫描4次，三维情况下需扫描8次，算法 效率极高。

3、更新符号距离函数；根据初始化给出的网格状态，判断其符号距离函数是否需要更新。若网格状态为“1”，将扫描迭代后的符号距离函数值与原始符号距离函数值进行比较，取最小值。如果网格状态为“0”，则取原始值。

第五步：采用WENO格式获取符号距离函数梯度

符号距离函数重构后，可采用高精度的WENO格式获取梯度，以五阶WENO为例说明其实现过程。

具体实现步骤如下：

1、选取插值模板，左导数模板为{I-2、I-1、I、I+1、I+2}、右导数模板为{I+2、I+1、I、I-1、I-2}；

2、构造三个三阶逼近s₁、s₂、s₃，以左导数为例，其插值模板分别为{I-2、I-1、I}、{I-1、I、I+1}、{I、I+1、I+2}；

3、定义三个非线性权ω₁、ω₂、ω₃；

4、定义三个光滑因子β₁、β₂、β₃；

5、构造三个三阶逼近的非线性凸组合；

第六步：用符号距离函数的梯度定义界面法向量

定义界面法向量时，必须要考虑迎风特性。若流场速度为负，取符号距离函数的右导数为界面法向量；若流场速度为正，则取符号距离函数的左导数为界面法向量。

第七步：求解VOF传输方程

本发明中采用THINC算法求解VOF传输方程。

具体实现步骤如下：

1、根据界面速度方向，定义λ和迎风网格iup；

2、计算双曲正切函数的跳跃中心

3、计算网格面流体体积流出量f_i+1/2、f_i-1/2；

4、将一维结果分裂至多维。

原THINC法中，界面法向量计算采用了Youngs的方法，此方法只具有一阶精度， 且不易扩展到三维。本发明采用精度更高的符号距离函数梯度计算界面法向量，可以极大地提高THINC算法的精度。

下面以二维旋转流场中的Zaleska圆盘界面运动为例，说明本发明的实施过程。

本例中，网格数为100×100，单位长度为0.01mm；Zaleska圆盘半径为0.3mm，圆心位于(0.5，0.5)，开口处高度为0.7mm，具体分布如图1所示。圆盘内VOF值为设为1，外部为0，界面处为0～1。

第二步：获取界面网格

首先通过网格的VOF函数值，初始化符号距离函数。VOF函数值大于0.5的网格的符号距离函数为正，VOF函数值小于0.5的网格的符号距离函数为负，其数学表达式如下。

其中，M为计算区域边长，F_i,j为当前网格的VOF函数值。

其次，采用子网格技术获判断当前网格是否为界面网格。如图2所示，对于任意网格Φ_i,j，假设其网格中心在运动界面的内部(Φ<0)，若相邻网格的符号距离函数值为正(Φ>0)，那么运动界面必包含在两网格之间，此网格为界面网格。基于上述判断条件，可获得整个计算区域内所有的界面网格。

第二步：将VOF函数转换为符号距离函数

由于VOF界面与Level Set界面重合，那么在界面网格处的符号距离函数的Heaviside值应等同于VOF函数值，其数学表达式如下：

H(φ)＝F

对其展开：

其中，F为VOF函数、Φ为符号距离函数、ε为界面宽度，本实施例中取1倍的网格长度。

采用牛顿迭代法可对上述等式进行求解，其迭代格式为：

其中，H(Φ)、δ(Φ)分别为符号距离函数的Heaviside和Dirac函数值。上述求解过程中，控制迭代停止的极小值为eps＝1e^-8。

第四步：重构符号距离函数

首先，根据第二步的结果对网格进行标记。若为界面网格，其网格状态为“0”；反之，网格状态则为“1”。

其次，对于“1”网格，若Φ>0，赋极大值；Φ<0，赋极大值；“0”网格保留原始值。

然后，沿4个方向进行高斯赛德尔迭代，其迭代格式如下：

其中，h为网格尺寸，a、b为相邻网格的极值，其表达式如下：

图3为不同方向扫描更新的网格区域，虚线表示物理界面。

最后，更新符号距离函数。若网格状态为“0”，不更新；否则，按下式更新：

图4所示为更新后的界面网格符号距离函数值。

第五步：获取符号距离函数梯度

首先，对于任意网格Φ_i,j，如图5所示，定义符号距离函数左右导数的插值模板。 以左导数为例，其表达式如下：

v₁＝D^-φ_i-2

v₂＝D^-φ_i-1

v₃＝D^-φ_i

v₄＝D^-φ_i+1

v₅＝D^-φ_i+2

其中，ν₁、ν₂、ν₃、ν₄、ν₅为五个插值模板，D^-为左导数。

其次，定义三个三阶逼近，其表达式如下：

再次，定义三个非线性权，其表达式如下：

其中，ε为防止分母为零的参数，本例中取1e^-10，S1、S2、S3为光滑因子，表达式如下：

最后，构造三个三阶逼近的非线性凸组合。

第六步：用符号距离函数梯度定义界面法向量

根据迎风特性，若流场速度为负，取符号距离函数的右导数为界面法向量；若流场速度为正，取符号距离函数的左导数为界面法向量。

第七步：求解VOF传输方程

首先，离散VOF方程：

其中，为当前时刻下的i网格内流体体积函数，为下一时刻的流体体积函数，f_i+1/2和f_i-1/2分别为网格右界面和左界面流体体积流出量，u_i+1/2和u_i-1/2分别为网格右界面和网格左界面速度，Δt和Δx分别为时间和空间步长。

由上式可知，VOF方法的关键在于求解网格面的流体体积流出量，而THINC方法正是通过解析方法获得网格面的流体体积流出量。如图6(a)和6(b)所示，图6(a)坐标中粗实线表示真实的物理界面，而图6(b)坐标中粗实线表示THINC构造的界面，即双曲正切函数，若已知流场速度和时间步长，对双曲正切函数求积分即可获得网格面的流体体积流出量。

双曲正切函数表达式如下：

其中，β是与法向量有关的量，γ是与速度方向有关的量，x_i-1/2为网格左界面坐标， 为双曲正切函数的跳跃中心。

β表达式为：

β_x＝2.3|n_x|+0.01

β_y＝2.3|n_y|+0.01

γ表达式为：

表达式为：

其中，Φ_iup为迎风网格，以右界面为例，若u_i+1/2<0，iup＝i+1，否则iup＝i；

值得注意的是，在原THINC算法中，nx、ny是由Youngs方法得到的界面法向量，本发明中采用第六步构造的符号距离函数梯法向量。

其次，对双曲正切函数求时间积分，可获得一定时间内流出网格面的流体体积。

以网格右界面为例，其流出流体体积为：

其中，λ为是与速度方向有关的量，其表达式如下：

左界面的流出体积为左网格右界面的流出体积，故可将上述结果代入离散的VOF方程中，获得更新后的VOF函数。

最后，采用分裂算法沿x、y、z三个方向分裂，对多维VOF函数进行更新。

分裂算法表达式如下：

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的保护范围。