长安大学 理学院, 陕西 西安 710064
收稿日期：2021-01-27
基金项目：国家自然科学基金资助项目(11901051, 11971075); 陕西省自然科学基础研究计划项目(2020JQ-338); 中央高校基本科研业务费专项资金资助项目(300102120302)。
作者简介：高普阳(1991-)，男，陕西西安人，长安大学讲师，博士。
摘要：针对自由界面问题, 构建了拉格朗日粒子和流体体积(volume of fluid, VOF)耦合算法.拉格朗日粒子方法可以准确追踪运动界面, 但是一般很难保证流体的质量守恒性.VOF方法可以保证很好的质量守恒性, 但是不容易计算界面的几何信息.因此, 本文构造了一种耦合算法, 吸收两种方法的优点.耦合算法中还引入了四叉树自适应网格技术, 可以在大变形区域提高界面的分辨率, 并能减少计算量.利用耦合算法模拟经典的Zalesak旋转盘问题和单涡剪切流动, 数值结果和文献已有结果吻合较好, 验证了耦合算法的稳定性、有效性和准确性.
关键词：拉格朗日粒子方法流体体积方法自由界面问题四叉树自适应网格Zalesak旋转盘单涡剪切流
Coupled Algorithm of Lagrangian Particle Method and Volume of Fluid Method for Free Interface Problem
GAO Pu-yang
School of Science, Chang'an University, Xi'an 710064, China
Corresponding author: GAO Pu-yang, E-mail: gaopuyang@chd.edu.cn.
Abstract: For free interface problem, a coupled algorithm of Lagrangian particle method and volume of fluid(VOF) was proposed. Although the Lagrangian particle method is able to accurately capture the interface front, there is no mechanism to guarantee the mass conservation. The volume of fluid method is known for its good mass conservation property. However, it is not easy to calculate the geometrical information of the interface. In the coupled method, the advantages of these two methods were combined. The technique of quad-tree adaptive grid was also added, which is able to track interface with high resolution and also save the computational cost. The traditional examples, i.e., Zalesak's slotted disk and reversed single vortex flow, are employed to illustrate the stability accuracy and efficiency of our coupled algorithm for solving dynamic interface problems. The numerical results agree with the existing results in the literature.
Key words: Lagrangian particle methodvolume of fluid(VOF)methodfree interface problemquad-tree adaptive gridZalesak disksingle vortex
自由界面问题出现在流体力学中的诸多领域中, 如溃坝流动[1]、聚合物充填[2]、水波晃动[3]、气泡在水中上升[4]等, 一直都是****们关注的研究热点.总的来说, 根据表示界面的方式不同, 可以将现有算法分为两大类: 拉格朗日粒子界面追踪方法[5-7]和欧拉界面捕捉方法[8-10].
VOF(volume of fluid)方法[11]属于欧拉界面捕捉类方法, 计算过程中显示追踪每个单元内的体积分数, 因此具有较好的质量守恒性, 这也成为该方法最大的优点.然而, 如何依据每个单元内的体积分数对界面进行重构, 并准确计算界面的几何信息比较困难.1981年, Hirt等[11]提出了分片线性的界面重构法; 2002年, Renardy等[12]构造了更高阶的重构法; 2003年, Scardovelli等[13]提出了最小二乘重构法.但是上述方法都需要同时重构界面位置、外法线方向以及界面曲率等几何信息.2005年, Cummins等[14]提出了直接通过体积分数估计几何信息的方法, 但是实施过程仍然相当繁琐.为了解决这些问题, 另外一种方式就是在VOF方法的基础上耦合Level Set方法, 这时需要求解一个额外的Level Set方程, 然后通过Level Set函数就可以方便计算出界面的几何信息.
2014年, Hon等[15]提出了基于单元的粒子方法(cell based particle method, CBPM)追踪运动界面, 该方法属于拉格朗日界面追踪方法.其中, 运动界面由拉格朗日粒子进行追踪, 同时还利用了固定的欧拉背景网格.在每一个和界面相交的单元内只能有一个拉格朗日粒子.这个方法可以很自然地处理界面大变形, 并且能够精确计算界面的几何信息, 如外法线方向、曲率等.然而, 这个方法最大的缺点就是不能很好地保证流体的质量守恒性.
本文主要目的就是结合已有方法的优点, 提出CBPM和VOF耦合算法, 用以准确地追踪运动界面.在耦合算法中, 界面由CBPM进行追踪, 并利用VOF方法计算每个单元的流体体积.然后, 依据单元内的流体体积修正粒子的位置, 以保证流体的质量守恒性质.相比于VOF方法, 耦合算法可以容易获取更为精确的界面几何信息; 相比于CBPM, 耦合算法能够很好地保证质量守恒性.另外, 为了提高大变形区域界面的分辨率, 同时减少一定的计算量, 耦合算法中还引入了四叉树自适应网格技术[16-17].
1 数值方法介绍1.1 基于单元的粒子方法CBPM[15]的核心思想就是采用准均匀的无网格粒子表示和追踪运动界面, 同时还利用一套固定的欧拉网格.首先找到所有和界面相交的网格单元, 然后以这些单元中心点到界面的投影点作为该单元内的无网格粒子.这时, 每一个和界面相交的单元内都仅有一个这样的粒子, 同时这些单元在计算过程中会作为粒子的欧拉参考.下面简单总结CBPM的具体计算过程, 其他更详细的信息可以参考文献[15].
CBPM具体的计算步骤如下:
1)初始化.找到所有和界面相交的网格单元, 将他们记为激活单元.然后计算这些激活单元中心点到界面的L∞投影点.显然在每一个激活单元内, 可以得到界面上一个相关的粒子点, 将它叫作角点.这时界面的形状就可以自然地由这些无网格角点来表示.
2) 移动粒子.根据已知速度, 在一个时间步内移动所有的角点.
3) 重新采样.在步骤2)移动粒子之后, 依据粒子的新位置局部重构界面形态.随后重新计算每一个激活单元中心点到局部重构界面上的L∞投影点.
4) 更新激活单元.将激活单元周围所有的单元也激活, 并根据步骤1)中的方法计算他们的角点.另外, 将那些不与界面相交的激活单元删除.
5) 自适应.如果需要, 将欧拉背景网格进行局部加密.
6) 迭代循环.重复上述步骤2)~步骤5), 直到最后时刻.
1.2 流体体积方法下面以分段线性界面为例, 介绍流体体积方法.在每一个单元中, 依据流体体积分数, 采用线性线段来近似界面.采用φij来表示单元-(i, j)内的流体体积分数, φij=0表示该单元内不含有流体, φij=1表示该单元内充满了流体, φij的取值在0和1之间, 则表示该单元内包含一部分流体界面.标量函数φ满足下面的输运方程:
 | (1) |
对于不可压缩流体而言, 上述方程也可写成如下形式:
 | (2) |
 | (3) |
 | (4) |
 | (5) |
 | (6) |
图 1(Fig. 1)
 | 图 1 x方向上流体输运示意图Fig.1 The schematic diagram of fluid transport in x-direction |
 | (7) |
2 耦合算法本文提出CBPM和VOF的耦合算法, 整个过程中利用CBPM显式追踪运动界面, 并计算界面外法线方向等几何信息; 通过VOF方法保证流体的质量守恒性以及界面拓扑变化的稳定性.下面给出耦合算法的具体步骤:
1) 初始化CBPM里面的角点和VOF中的流体体积分数.首先按照1.1节CBPM步骤1)的方式计算角点及该点处的外法线方向, 然后利用过角点切线和单元边界所组成区域的面积近似该激活单元内流体体积.
例如, 图 2中的灰色区域表示初始时刻该激活单元内的流体, 为了减少计算量, 计算时采用三角形ABC的面积来近似流体体积.
图 2(Fig. 2)
 | 图 2 激活单元内初始时刻的流体Fig.2 Initial fluid of the active cell |
2) 移动角点并且更新流体体积函数.在每一个时间步, 根据已知的速度场移动所有的角点.同时, 采用1.2节的输运算法更新流体体积函数, 这样就可以得到所有单元内最新的流体体积分布情况.
3) 重新选取角点.在步骤2)中移动角点之后, 根据角点新位置对界面进行局部重构, 并计算由每一个激活单元中心到局部重构界面的L∞投影点, 使其成为新的角点.
4) 更新激活单元.对于激活单元所有的相邻单元, 判断它们是否和界面相交.如果相交, 将它们激活并计算它们的角点, 同时将之前激活单元中不再与界面相交的单元, 从激活单元序列中删除.
5) 修正角点位置.经过上述步骤4)之后, 能够得到新的角点以及它们的外法线方向, 由此依据步骤1)中的方法, 可以容易估算出这些单元内的流体体积.由CBPM计算出的流体体积和VOF的计算结果之间可能会有一些差异, 这时需要依据下面的算式对角点的位置进行修正，
 | (8) |
 | (9) |
图 3(Fig. 3)
 | 图 3 修正角点位置Fig.3 Modification of the footpoint |
6) 自适应.为了更好地捕捉锋利的尖角或者界面较薄区域, 必要时需要对欧拉背景网格进行局部加密.以界面曲率作为条件, 用来判断对欧拉网格中的哪些单元进行加密.对于每一个激活单元-(i, j), 检测如下的条件是否成立[18-19]:
 | (10) |
图 4(Fig. 4)
 | 图 4 大曲率区域的自适应网格Fig.4 Adaptation for high curvature area ●—激活单元中心; ■—角点. |
7) 迭代循环.重复之前的步骤2)~步骤5), 直至计算结束.
注1 计算过程中会存在这样一种情况, 某个单元内的流体体积分数不为零, 但是在该单元内没有任何界面上的角点.因此, 需要激活这些单元, 并且合理地在每一个单元内增加一个角点.这个过程会使得最终的结果更加精确, 并能降低最后的计算误差.这个单元内角点处的外法线方向, 可以由相邻三个激活单元角点处外法线方向的平均值进行近似计算.另外如图 5所示, 可以根据外法线方向和流体体积分数φ计算出切线AB, 使得灰色区域的面积和流体面积相等, 然后将AB的中点O选取为这个单元内的角点.
图 5(Fig. 5)
 | 图 5 增加角点的位置Fig.5 The position of added footpoint |
注2 在上面耦合算法的步骤5)中, 必须保证修正后的角点仍然落在同一个网格单元内.绝大多数情况下, 因为基本都是微小修正, 这个条件很容易满足.然而, 有一个极端情况, 不一定能满足上面的要求, 因此需要额外的处理方式.如图 6所示, O是单元中之前的角点, P是依据方程(8)和(9)修正之后的角点位置.显然这时P点已经在网格单元之外了, 不符合耦合算法步骤5)的要求, 在之后的计算中会产生额外的误差.因此, 将线段AB的中点C作为这个单元最终的角点.
图 6(Fig. 6)
 | 图 6 一个坏的情形Fig.6 A bad situation |
3 数值算例3.1 Zalesak旋转圆盘Zalesak旋转圆盘问题[20], 经常被用来检测数值算法捕捉带有奇异点自由界面的能力.在这个基准算例中, 将1个Zalesak圆盘进行旋转, 其初始界面形态和位置如图 7所示, 圆盘的半径为0.15, 圆心坐标为(0.5, 0.75), 圆盘内部的空隙宽度为0.05, 旋转速度为u=1-2y, v=2x-1.根据给定的速度, 圆盘会绕着点(0.5, 0.5)进行旋转, 并且在一个周期t=π之后回到初始的位置.整个旋转过程中对于圆盘界面的捕捉, 最大的挑战就是界面上的4个奇异点.
图 7(Fig. 7)
 | 图 7 Zalesak圆盘初始时刻界面形态Fig.7 The initial profile of Zalesak's disk |
本文计算了一个周期内圆盘的旋转情况, 并在图 8中给出了256×256均匀网格下的数值结果, 图 8a表示CBPM的计算结果, 图 8b表示本文耦合算法的计算结果, 图 8c表示耦合算法中带有自适应网格技术的计算结果, 图 8d~8f分别代表三种算法得到数值结果的局部放大图, 其中黑色角点构成的图形表示旋转一个周期之后圆盘的形状, 可以看到带有自适应网格技术耦合算法的数值结果最为准确.
图 8(Fig. 8)
 | 图 8 256×256均匀网格上圆盘旋转一个周期之后的数值结果Fig.8 The numerical results after one cycle of disc in uniform grid of 256×256 (a)—CBPM; (b)—耦合算法; (c)—带有自适应网格技术的耦合算法; (d)—图a局部放大图；(e)—图b局部放大图；(f)—图c局部放大图. |
误差定义为
 | (11) |
表 1(Table 1)
 表 1 不同算法在方程(11)定义下的误差Table 1 Errors defined in Eq.(11) under different algorithms
| 表 1 不同算法在方程(11)定义下的误差 Table 1 Errors defined in Eq.(11) under different algorithms |
3.2 单涡剪切流动本节数值模拟单涡剪切流动问题[21], 用来验证耦合算法追踪大变形界面以及较薄区域的能力.考虑图 9中的初始单涡圆界面, 圆心位于(0.5, 0.75), 半径为0.15.速度场取为
图 9(Fig. 9)
 | 图 9 初始单涡圆的位置Fig.9 The position of the single vortex at the initial time |
 |
本文计算了一个周期内单涡的运动过程, 图 10a和图 10b分别表示在512×512均匀网格上, 由CBPM计算得到的t=4和t=8时刻单涡的界面形态.实线线条表示数值结果, 可以明显看到最后时刻数值结果和初始界面之间有一些差异, 流体的质量出现了损失.
图 10(Fig. 10)
 | 图 10 不同时刻CBPM计算得到的单涡界面的形状Fig.10 The interface of single vortex via CBPM at different time instants (a)— t=4; (b)— t=8. |
图 11和图 12分别为CBPM-VOF耦合算法及自适应CBPM-VOF耦合算法的数值结果.由图 12可知, 数值结果在t=4时捕捉界面尾端较薄区域的能力更强, 而且在t=8时和初始界面的吻合程度也更好.
图 11(Fig. 11)
 | 图 11 不同时刻CBPM和VOF耦合算法的单涡界面的形状Fig.11 The interface of single vortex via CBPM coupled with VOF method at different time instants (a)— t=4; (b)— t=8. |
图 12(Fig. 12)
 | 图 12 不同时刻自适应耦合算法的单涡界面的形状Fig.12 The interface of single vortex via adaptive coupled method at different time instants (a)—t=4; (b)—t=8. |
由表 2中的数据可以看到, 根据方程(11)定义的误差公式, 耦合算法的误差要稍微高于CBPM一些.然而, 耦合算法的质量误差要远远小于CBPM.另外, 相比于CBPM, 自适应耦合算法不但可以保证较好的质量守恒性, 同时也能减少最终计算结果和精确解在图形上的几何误差.
表 2(Table 2)
表 2 不同算法在方程(11)定义下的误差Table 2 Errors defined in Eq.(11) under different algorithms
| 表 2 不同算法在方程(11)定义下的误差 Table 2 Errors defined in Eq.(11) under different algorithms |
4 结论本文提出了CBPM-VOF耦合算法, 采用拉格朗日粒子追踪运动界面, 并通过VOF方法保证流体的质量守恒性.耦合算法可以稳定、准确地追踪带有奇异点以及较薄区域的自由界面, 所有的数值结果和文献中已有结果吻合较好.耦合算法所产生的质量误差要小于CBPM, 另外自适应耦合算法不但能够保证较好的质量守恒性, 而且可以降低数值结果的几何误差.
参考文献
| [1] | Issakhov A, Imanberdiyeva M. Numerical simulation of the movement of water surface of dam break flow by VOF methods for various obstacles[J]. International Journal of Heat and Mass Transfer, 2019, 136: 1030-1051. DOI:10.1016/j.ijheatmasstransfer.2019.03.034 |
| [2] | Gao P, Wang X, Ouyang J. Numerical investigation of non-isothermal viscoelastic filling process by a coupled finite element and discontinuous Galerkin method[J]. International Journal of Heat and Mass Transfer, 2019, 140: 227-242. DOI:10.1016/j.ijheatmasstransfer.2019.05.115 |
| [3] | Green M D, Peiro J. Long duration SPH simulations of sloshing in tanks with a low fill ratio and high stretching[J]. Computers & Fluids, 2018, 174: 179-199. |
| [4] | Gao P, Ouyang J, Zhou W. Development of a finite element/discontinuous Galerkin/level set approach for the simulation of incompressible two phase flow[J]. Advances in Engineering Software, 2018, 118: 45-59. DOI:10.1016/j.advengsoft.2018.01.006 |
| [5] | Kahouadji L, Nowak E, Kovalchuk N M, et al. Simulation of immiscible liquid-liquid flows in complex microchannel geometries using a front-tracking scheme[J]. Microfluidics and Nanofluidics, 2018, 22(11): 1-12. |
| [6] | Zhang S, Wang S, Zhang A, et al. Numerical study on motion of the air-gun bubble based on boundary integral method[J]. Ocean Engineering, 2018, 154: 70-80. DOI:10.1016/j.oceaneng.2018.02.008 |
| [7] | Rowlatt C F, Lind S J. Bubble collapse near a fluid-fluid interface using the spectral element marker particle method with applications in bioengineering[J]. International Journal of Multiphase Flow, 2017, 90: 118-143. DOI:10.1016/j.ijmultiphaseflow.2016.11.010 |
| [8] | Gibou F, Fedkiw R, Osher S. A review of level-set methods and some recent applications[J]. Journal of Computational Physics, 2018, 353: 82-109. DOI:10.1016/j.jcp.2017.10.006 |
| [9] | Zografos K, Afonso A M, Poole R J, et al. A viscoelastic two-phase solver using a phase-field approach[J]. Journal of Non-Newtonian Fluid Mechanics, 2020, 234: 104364. |
| [10] | Shams M, Raeini A Q, Blunt M J, et al. A numerical model of two-phase flow at the micro-scale using the volume-of-fluid method[J]. Journal of Computational Physics, 2018, 357: 159-182. DOI:10.1016/j.jcp.2017.12.027 |
| [11] | Hirt C W, Nichols B D. Volume of fluid(VOF)method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201-225. DOI:10.1016/0021-9991(81)90145-5 |
| [12] | Renardy Y, Renardy M. PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method[J]. Journal of Computational Physics, 2002, 183(2): 400-421. DOI:10.1006/jcph.2002.7190 |
| [13] | Scardovelli R, Zaleski S. Interface reconstruction with least‐square fit and split Eulerian-Lagrangian advection[J]. International Journal for Numerical Methods in Fluids, 2003, 41(3): 251-274. DOI:10.1002/fld.431 |
| [14] | Cummins S J, Francois M M, Kothe D B. Estimating curvature from volume fractions[J]. Computers & Structures, 2005, 83(6/7): 425-434. |
| [15] | Hon S Y, Leung S, Zhao H. A cell based particle method for modeling dynamic interfaces[J]. Journal of Computational Physics, 2014, 272: 279-306. DOI:10.1016/j.jcp.2014.04.032 |
| [16] | Ghazizadeh M A, Mohammadian A, Kurganov A. An adaptive well-balanced positivity preserving central-upwind scheme on quadtree grids for shallow water equations[J]. Computers & Fluids, 2020, 208: 104633. |
| [17] | Teunissen J, Keppens R. A geometric multigrid library for quadtree/octree AMR grids coupled to MPI-AMRVAC[J]. Computer Physics Communications, 2019, 245: 106866. DOI:10.1016/j.cpc.2019.106866 |
| [18] | Malik M, Fan E S C, Bussmann M. Adaptive VOF with curvature‐based refinement[J]. International Journal for Numerical Methods in Fluids, 2007, 55(7): 693-712. DOI:10.1002/fld.1490 |
| [19] | Leung S, Zhao H. A grid based particle method for moving interface problems[J]. Journal of Computational Physics, 2009, 228(8): 2993-3024. DOI:10.1016/j.jcp.2009.01.005 |
| [20] | An R, Yu C. A level set redistancing algorithm for simulation of two-phase flow[J]. Numerical Heat Transfer, Part B: Fundamentals, 2020, 78(1): 1-24. |
| [21] | Zhou W, Ouyang J, Zhang L, et al. Development of new finite volume schemes on unstructured triangular grid for simulating the gas-liquid two‐phase flow[J]. International Journal for Numerical Methods in Fluids, 2016, 81(1): 45-67. DOI:10.1002/fld.4174 |
