买专利卖专利找龙图腾，真高效！ 查专利查商标用IPTOP,全免费！专利年费监控用IP管家,真方便！

摘要：针对现有CLSVOF方法中符号距离函数构造过程过于复杂、不易扩展到三维，构造效率低等问题，本发明提出一种由VOF函数快速构造符号距离函数的改进CLSVOF方法，该方法首先利用子网格技术找出界面网格，再通过牛顿迭代法将界面网格内的VOF函数转换为符号距离函数，在此基础上，采用快速扫描法构造全计算区域的符号距离函数，然后通过WENO格式和迎风格式计算获取准确的界面法向量值，最后采用修正的THINC算法对求解VOF方程。

主权项：一种由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传输方程。

全文数据：一种由VOF函数快速构造符号距离函数的改进CLSVOF方法技术领域：[0001]本发明涉及一种由V0F函数快速构造符号距离函数的改进CLSVOF方法，其属于运动界面追踪数值模拟领域。背景技术：[0002]运动界面追踪数值模拟是当前国内外热门研究课题，它涉及到材料、物理、化学、生物等多个研究领域，对其开展研究具有重要的理论价值和现实意义。[0003]上世纪80年代初，美国洛斯阿拉莫斯国家实验室的Hirt和Nichols首先提出了流体体积法VolumeofFluids，V0F。该方法在流场中引入一个流体体积V0F函数其值为0~1，通过求解V0F函数的运输方程，实现运动界面的追踪。V0F法计算量小，界面锐利程度高，同时具有良好的质量守恒性，因此被大量商用软件及国内外学者采用，用于求解复杂的运动界面问题。遗憾的是，由于V0F函数在界面两侧存在间断，无法准确计算出界面的法向量，降低了该方法的精度。[0004]1988年，加州大学伯克利分校的Osher和Sethian基于Hamilton-Jacobi方程的思想提出了著名的水平集法LevelSetMethod，此方法定义了一个连续光滑的LevelSet函数（如符号距离函数，SignedDistanceFunctions，SDF，通过求解LevelSet函数的运输方程，实现运动界面追踪。Leve1Set法可以很方便地计算复杂界面变化如破碎、合并等；并且由于符号距离函数的连续特性，结合高分辨率类方法如WEN0、TVD等格式，可以很准确地计算出界面的法向量。然而，由于数值扩散和重新初始化等原因，LevelSet法在计算过程中容易产生质量不守恒。[0005]V0F法的缺点在于V0F函数在界面两侧存在间断，无法准确计算出界面的法向量；而LevelSet法的优点则是LevelSet函数连续，可以准确计算出界面法向量;如果在V0F法中引入符号距离函数，用其构造界面法向量，那么必然可以有效弥补V0F法的缺点。[0006]2000年，加州大学戴维斯分校Sussamn和Puckett等人首先将V0F法和LevelSet法结合，提出了一种親合的V0F和LevelSet方法（CoupledLevelSetandVOFMethod，CLSVOF。该方法采用PLIC算法求解VOF方程，再由PLIC算法建立的网格内部线性界面二维为直线，三维为平面）构造符号距离函数。由于三维线性界面的建立过程很复杂，此方法在三维情况下难以实现，并且很难扩展到其它V0F求解方法。[0007]2010年，西安交通大学的孙东亮等人提出了一种V0SET方法。该方法采用几何方法计算各网格中心到PLIC线性界面界面的距离，再根据V0F函数的大小，将其转换为符号距离函数。与CLSVOF法类似，此方法同样不易扩展到三维。[0008]2013年，英国卡迪夫大学的Yoki等人提出了一种新的CLSVOF方法，该方法采用Xiao等人提出的THINC算法求解V0F传输方程。由于THINC算法中并无几何界面的建立过程，所以需要通过其他方法构造符号距离函数。为此，Yoki等人采用线性插值的方式构造出V0F函数的0.5等值面，再通过Sethian等人提出的快速推进法FastMarchingMethocUFMM及Sussman等人提出的稳态解法重新初始化Re-initialization构造符号距离函数，简化了符号距离函数的重构流程。然而，FMM算法逻辑过于复杂，且在构造界面网格处的符号距离函数时存在一定误差，而稳态解法的迭代步数较多，计算效率过低。发明内容：[0009]本发明针对CLSV0F方法中存在的符号距离函数构造过程过于复杂的问题，提出一种由V0F函数快速构造符号距离函数的改进CLSV0F方法。[0010]本发明采用如下技术方案：一种由V0F函数快速构造符号距离函数的改进CLSV0F方法，包括如下步骤：[0011]第一步:定义初始V0F函数[0012]根据实际的物理界面，在网格内定义V0F函数，界面网格的V0F函数值为0~1，界面内的网格的V0F函数值为1，界面外的网格的V0F函数值为0;[0013]第二步:获取界面网格[0014]首先通过网格的V0F函数值，初始化符号距离函数，V0F值大于0.5的网格的符号距离函数为正，V0F值小于0.5的网格的符号距离函数为负;其次，采用子网格技术判断当前网格是否为界面网格，所述子网格技术即通过判断当前网格与周围网格的符号距离函数值是否异号来确定其是否为界面网格，若异号，说明运动界面在两网格之间，网格为界面网格，反之，则为非界面网格；[0015]第三步:将V0F函数转换为符号距离函数[0016]由于V0F函数和符号距离函数代表了相同的界面，那么在界面网格处的符号距离函数Heaviside值应等同于V0F函数值，采用Newton迭代法对上述等同关系进行求解，可获得界面网格的符号距离函数，即符号距离函数值；[0017]第四步:重构符号距离函数[0018]采用快速扫描法对符号距离函数进行重构，所述快速扫描法是一种求解Eikonal方程的快速解法，对符号距离函数重构等同于求解Eikonal方程|νΦ|=1;[0019]第五步:采用WEN0格式获取符号距离函数梯度[0020]符号距离函数重构后，采用高精度的WEN0格式获取梯度；[0021]第六步:用符号距离函数的梯度定义界面法向量[0022]定义界面法向量时，考虑迎风特性，若流场速度为负，取符号距离函数的右导数为界面法向量，若流场速度为正，则取符号距离函数的左导数为界面法向量；[0023]第七步:求解V0F传输方程。[0024]进一步地，所述第四步中包括如下步骤[0025]1、初始化;根据第三步给出的结果，将界面处的网格状态设为"0"，非界面网格状态设为"1"，若非界面网格的符号距离函数值为正，则将其初始化为极大值，否则初始化为极小值；[0026]2、高斯赛德尔迭代;根据快速扫描法给出的格式，对整个计算区域沿不同方向进行扫描迭代；[0027]3、更新符号距离函数;根据初始化给出的网格状态，判断其符号距离函数是否需要更新，若网格状态为"1"，将扫描迭代后的符号距离函数值与原始符号距离函数值进行比较，取最小值，如果网格状态为"〇"，则取原始值。[0028]进一步地，所述第五步中包括如下步骤[0029]1、选取插值模板，左导数模板为{Ι-2、Ι-1、Ι、Ι+1、Ι+2}、右导数模板为{1+2、1+1、1、1-1、1-2};[0030]2、构造三个三阶逼近81、別、83，以左导数为例，其插值模板分别为{1-2、1-1、1}、{I-1、I、I+1}、{1、1+1、1+2};[0031]3、定义三个非线性权[0032]4、定义三个光滑因子&、&、&;[0033]5、构造三个三阶逼近的非线性凸组合。[0034]进一步地，所述第七步中包括如下步骤[0035]1、根据界面速度方向，定义λ和迎风网格iup;[0036]2、计算双曲正切函数的跳跃中心F;[0037]3、计算网格面流体体积流出量fi+12、fi-12;[0038]4、将一维结果分裂至多维。[0039]本发明具有如下有益效果：[0040]1本发明采用的子网格技术、牛顿迭代法、快速扫描法、WEN0格式、THINC算法逻辑简单，易于编程实现；[0041]2子网格技术和牛顿迭代法和对二维和三维均适用，而快速扫描法也很容易向三维扩展，故本方法可以很容易地应用到三维情况；[0042]3子网格技术只需一步计算即可完成，而牛顿迭代法一般和快速扫描法也仅需几步即可收敛，故本方法计算效率极高；[0043]⑷界面网格的符号距离函数由牛顿迭代法得到，具有机器精度，界面外网格的符号距离函数由快速扫描法得到，具有一阶精度，与现有算法相当；[0044]5本发明设计的符号距离函数构造方法不依赖于V0F求解方法，它不仅适用于THINC算法，亦可适用于SLIC、PLIC等传统的V0F方法，适用性强；[0045]6THINC算法是一种不依赖于时间步、数值稳定性强的算法，本发明完美地继承了这些优点，算法鲁棒性强。附图说明：[0046]图1为初始V0F函数分布示意图。[0047]图2为界面网格判断不意图。[0048]图3为界面网格的符号距离函数分布示意图。[0049]图4为快速扫描法更新网格符号距离函数示意图。[0050]图5为符号距离函数左右导数的插值模板示意图。[0051]图6a和6b为一维THINC算法原理示意图。具体实施方式：[0052]本发明由V0F函数快速构造符号距离函数的改进CLSV0F方法，包括如下步骤：[0053]第一步:定义初始V0F函数[0054]根据实际的物理界面，在网格内定义V0F函数。界面网格的V0F函数值为0~1，界面内的网格的VOF函数值为1，界面外的网格的VOF函数值为0。[0055]第二步:获取界面网格[0056]首先通过网格的V0F函数值，初始化符号距离函数。V0F值大于0.5的网格的符号距离函数为正，V0F值小于0.5的网格的符号距离函数为负。其次，采用子网格技术判断当前网格是否为界面网格。所谓子网格技术，即通过判断当前网格与周围网格的符号距离函数值是否异号来确定其是否为界面网格。若异号，说明运动界面在两网格之间，网格为界面网格;反之，则为非界面网格。[0057]第三步:将V0F函数转换为符号距离函数[0058]由于V0F函数和符号距离函数代表了相同的界面，那么在界面网格处的符号距离函数Heaviside值应等同于V0F函数值。采用Newton迭代法对上述等同关系进行求解，可获得界面网格的符号距离函数，即符号距离函数值。[0059]第四步:重构符号距离函数[0060]为了高效快速重构符号距离函数，本发明采用快速扫描法对符号距离函数进行重构，快速扫描法是一种求解Eikonal方程的快速解法，对符号距离函数重构等同于求解Eikonal方程|▽Φ|=1〇[0061]其实现步骤如下：[0062]1、初始化;根据第三步给出的结果，将界面处的网格状态设为"0"，非界面网格状态设为"1"。若非界面网格的符号距离函数值为正，则将其初始化为极大值;否则初始化为极小值。[0063]2、高斯赛德尔迭代;根据快速扫描法给出的格式，对整个计算区域沿不同方向进行扫描迭代）。一般来说，二维情况下仅需扫描4次，三维情况下需扫描8次，算法效率极尚。[0064]3、更新符号距离函数;根据初始化给出的网格状态，判断其符号距离函数是否需要更新。若网格状态为"1"，将扫描迭代后的符号距离函数值与原始符号距离函数值进行比较，取最小值。如果网格状态为"〇"，则取原始值。[0065]第五步:采用WEN0格式获取符号距离函数梯度[0066]符号距离函数重构后，可采用高精度的WEN0格式获取梯度，以五阶WEN0为例说明其实现过程。[0067]具体实现步骤如下：[0068]1、选取插值模板，左导数模板为{Ι-2、Ι-1、Ι、Ι+1、Ι+2}、右导数模板为{1+2、1+1、1、1-1、1-2};[0069]2、构造三个三阶逼近81、別、83，以左导数为例，其插值模板分别为{1-2、1-1、1}、{I-1、I、I+1}、{1、1+1、1+2};[0070]3、定义三个非线性权ω!、〇〗、ω3;[0071]4、定义三个光滑因子&、&、&;[0072]5、构造三个三阶逼近的非线性凸组合；[0073]第六步:用符号距离函数的梯度定义界面法向量[0074]定义界面法向量时，必须要考虑迎风特性。若流场速度为负，取符号距离函数的右导数为界面法向量;若流场速度为正，则取符号距离函数的左导数为界面法向量。[0075]第七步:求解V0F传输方程[0076]本发明中采用THINC算法求解V0F传输方程。[0077]具体实现步骤如下：[0078]1、根据界面速度方向，定义λ和迎风网格iup;[0079]2、计算双曲正切函数的跳跃中心X;[0080]3、计算网格面流体体积流出量fi+12、fi-12;[0081]4、将一维结果分裂至多维。[0082]原THINC法中，界面法向量计算采用了Youngs的方法，此方法只具有一阶精度，且不易扩展到三维。本发明采用精度更高的符号距离函数梯度计算界面法向量，可以极大地提高THINC算法的精度。[0083]下面以二维旋转流场中的Zaleska圆盘界面运动为例，说明本发明的实施过程。[0084]本例中，网格数为100X100，单位长度为0·01mm;Zaleska圆盘半径为0·3mm，圆心位于0.5,0.5，开口处高度为0.7_，具体分布如图1所示。圆盘内V0F值为设为1，外部为0，界面处为0~1。[0085]第二步:获取界面网格[0086]首先通过网格的V0F函数值，初始化符号距离函数。V0F函数值大于0.5的网格的符号距离函数为正，V0F函数倌小于0.5的网格的符号距离函数为负，其数学表达式如下。[0088]其中，Μ为计算区域边长，Fi,j为当前网格的V0F函数值。[0089]其次，采用子网格技术获判断当前网格是否为界面网格。如图2所示，对于任意网格Φ^，假设其网格中心在运动界面的内部（Φ〈〇，若相邻网格的符号距离函数值为正（Φ0，那么运动界面必包含在两网格之间，此网格为界面网格。基于上述判断条件，可获得整个计算区域内所有的界面网格。[0090]第二步:将V0F函数转换为符号距离函数[0091]由于V0F界面与LevelSet界面重合，那么在界面网格处的符号距离函数的Heaviside值应等同于V0F函数值，其数学表达式如下：[0092]ΗΦ=F[0093]对其展开：[0095]其中，F为V0F函数、Φ为符号距离函数、ε为界面宽度，本实施例中取1倍的网格长度。[0096]采用牛顿迭代法可对上述等式进行求解，其迭代格式为：L〇〇98J其中，ΗΦ、δΦ分别为符号距离函数的Heaviside和Dirac函数值。上述求解过程中，控制迭代停止的极小值为eps=1e4。[0099]第四步:重构符号距离函数[0100]首先，根据第二步的结果对网格进行标记。若为界面网格，其网格状态为"0"；反之，网格状态则为"1"。[0101]其次，对于"Γ网格，若φ0,赋极大值；φ〈0,赋极大值；"0"网格保留原始值。[0102]然后，沿4个方向进行高斯赛德尔迭代，其迭代格式如下：[0105]其中，h为网格尺寸，a、b为相邻网格的极值，其表达式如下：[0107]图3为不同方向扫描更新的网格区域，虚线表示物理界面。[0108]最后，更新符号距离函数。若网格状态为"〇"，不更新;否则，按下式更新：[0110]图4所示为更新后的界面网格符号距离函数值。第五步:获取符号距离函数梯度[0112]首先，对于任意网格Φ^，如图5所示，定义符号距离函数左右导数的插值模板。以左导数为例，其表达式如下：[0113]V1=D-Φi-2[0114]V2=D"Φi-i[0115]V3=D"Φi[0116]V4=D"Φi+i[0117]V5=D"Φi+2[0118]其中，乂1^2^3^4^5为五个插值模板，『为左导数。[0120]其次，定义三个三阶逼近，其表达式如下：[0124]再次，定义三个非线性权，其表达式如下：[0131]其中，ε为防止分母为零的参数，本例中取1^1，31、32、33为光滑因子，表达式如下：[0135]最后，构造三个三阶逼近的非线性凸组合。[0136]處[0137]第六步:用符号距离函数梯度定义界面法向量[0138]根据迎风特性，若流场速度为负，取符号距离函数的右导数为界面法向量;若流场速度为正，取符号距离函数的左导数为界面法向量。[0139]第七步:求解V0F传输方程[0140]首先，离散V0F方程：[0143]其中，f1为当前时刻下的i网格内流体体积函数，月8+1为下一时刻的流体体积函数，fi+l2和fi-l2分别为网格右界面和左界面流体体积流出量，Ui+l2和Ui-l2分别为网格右界面和网格左界面速度，At和△X分别为时间和空间步长。[0144]由上式可知，V0F方法的关键在于求解网格面的流体体积流出量，而THINC方法正是通过解析方法获得网格面的流体体积流出量。如图6a和6b所示，图6a坐标中粗实线表示真实的物理界面，而图6b坐标中粗实线表示THINC构造的界面，即双曲正切函数，若已知流场速度和时间步长，对双曲正切函数求积分即可获得网格面的流体体积流出量。[0145]双曲正切函数表达式如下：[0147]其中，β是与法向量有关的量，γ是与速度方向有关的量，Xi-12为网格左界面坐标，f为双曲正切函数的跳跃中心。[0148]β表达式为：[0151]γ表达式为：[0153]无表达式为：[0155]其中，4_为迎风网格，以右界面为例，若出+12〈0，：[即=丨+1，否则;[即=;[;[0156]值得注意的是，在原THINC算法中，nx、ny是由Youngs方法得到的界面法向量，本发明中采用第六步构造的符号距离函数梯法向量。[0157]其次，对双曲正切函数求时间积分，可获得一定时间内流出网格面的流体体积。[0158]以网格右界面为例，其流出流体体积为：[0160]其中，λ为是与速度方向有关的量，其表达式如下：[0162]左界面的流出体积为左网格右界面的流出体积，故可将上述结果代入离散的V0F方程中，获得更新后的V0F函数。[0163]最后，采用分裂算法沿x、y、z三个方向分裂，对多维V0F函数进行更新。[0164]分裂算法表达式如下：[0168]以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的保护范围。

权利要求：1.一种由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方程IνΦI=1;第五步:采用WENO格式获取符号距离函数梯度符号距离函数重构后，采用高精度的WENO格式获取梯度；第六步:用符号距离函数的梯度定义界面法向量定义界面法向量时，考虑迎风特性，若流场速度为负，取符号距离函数的右导数为界面法向量，若流场速度为正，则取符号距离函数的左导数为界面法向量；第七步:求解VOF传输方程。2.如权利要求1所述的由VOF函数快速构造符号距离函数的改进CLSVOF方法，其特征在于:所述第四步中包括如下步骤1初始化;根据第三步给出的结果，将界面处的网格状态设为"〇"，非界面网格状态设为"1"，若非界面网格的符号距离函数值为正，则将其初始化为极大值，否则初始化为极小值；2高斯赛德尔迭代;根据快速扫描法给出的格式，对整个计算区域沿不同方向进行扫描迭代；3更新符号距离函数;根据初始化给出的网格状态，判断其符号距离函数是否需要更新，若网格状态为"1"，将扫描迭代后的符号距离函数值与原始符号距离函数值进行比较，取最小值，如果网格状态为"〇"，则取原始值。3.如权利要求2所述的由VOF函数快速构造符号距离函数的改进CLSVOF方法，其特征在于:所述第五步中包括如下步骤1选取插值模板，左导数模板为{1-2、I-I、I、1+1、1+2}、右导数模板为{1+2、1+1、I、I-U1-2}；2构造三个三阶逼近S1、S2、S3，以左导数为例，其插值模板分别为{1-2、I-I、1}、{1-1、I、I+1}、{Ι、Ι+1、Ι+2};3定义三个非线性权ω:、ω2、ω3;4定义三个光滑因子、&、β3;5构造三个三阶逼近的非线性凸组合。4.如权利要求3所述的由VOF函数快速构造符号距离函数的改进CLSVOF方法，其特征在于:所述第七步中包括如下步骤1根据界面速度方向，定义λ和迎风网格iup;2计算双曲正切函数的跳跃中心Ϊ;⑶计算网格面流体体积流出量fi+l2、fi-V2;4将一维结果分裂至多维。

百度查询： 南京航空航天大学 一种由VOF函数快速构造符号距离函数的改进CLSVOF方法