单层粒子水平集方法

牟凯龙，赵兰浩，毛 佳

(河海大学 水利水电学院，南京 210098)

包含运动界面的流动问题广泛存在于自然界中[1].对这类问题进行数值模拟的难点往往在于运动界面的描述，尤其是当运动界面经历剧烈的拓扑变化时，问题求解难度会进一步增加.因此，如何精确地表示运动界面一直是研究热点.

各国学者已提出多种数值方法来表示运动界面，大致可以分为拉格朗日方法和欧拉方法两类.拉格朗日网格类方法的网格与流体域完全重合并跟随流体运动，可以直接得到精确的运动界面.而且拉格朗日求解格式不存在对流项，因此这类方法具有良好的体积守恒性.但当界面变形严重时，畸形的背景网格会导致计算精度下降，而网格重生成又会消耗大量的计算资源.为克服这种缺陷，有学者提出了拉格朗日粒子类方法，如光滑粒子质点(Smoothed Particle Hydrodynamics,SPH)方法[2].这类方法采用粒子替代网格追踪运动界面，尽管能够精确追踪任意变形的界面，但SPH方法固有的张力不稳定性限制了该方法的应用.此外，由于界面由离散的粒子表示，所以难以得到光滑的界面.

基于固定背景网格的欧拉方法根据描述界面的形式又可以分为界面追踪方法和界面捕捉方法两类[3].标志网格(Marker and Cell,MAC)方法[4]是应用最早的界面追踪方法，主要思想是在背景网格中引入无质量的拉格朗日粒子追踪界面，最后根据粒子位置确定运动界面.因此，MAC方法也无法提供光滑的运动界面.此外，界面追踪方法难以处理界面复杂的拓扑变化[4].

另一类欧拉方法即界面捕捉方法通过设置指示函数来隐式地捕捉界面，能够自动处理复杂的拓扑变化应用.其中，最广泛的是流体体积(Volume of Fluid,VOF)方法[5]和水平集(Level Set,LS)方法[6-8].VOF方法通过在流场中定义流体体积函数，根据其函数值利用界面重构算法构造运动界面.虽然VOF方法凭借其良好的体积守恒性得到了广泛的应用，但是一方面需要复杂且计算成本较高的界面重构技术，另一方面由于体积函数具有不连续性，不能准确计算界面的法线方向和曲率.

相较于VOF方法，LS方法设置了一个连续的符号距离函数来隐式地捕捉界面，不仅能够自动处理复杂的拓扑变化，而且能够直接表示光滑界面，界面的几何特征也能够直接计算.但LS方法对数值耗散极度敏感，这也导致其体积守恒性较差[9].在LS方法中，数值耗散主要来源于两方面：控制方程的离散和重新初始化过程.前者在数值求解中难以避免，可通过采用高阶离散格式来降低其影响.重新初始化过程是为保持LS函数的符号距离特性而引入的[10]，但界面会发生偏移，引起体积损失.针对LS方法体积守恒性差的缺陷，各国学者已提出多种改进的方法.守恒式水平集(Conservative Level Set,CLS)方法[11]将符号距离函数替换为Heaviside函数，改善了其体积守恒性，但重新初始化过程仍然会导致体积损失.考虑到VOF方法和LS方法具有互补的优缺点，有学者将二者结合提出了耦合式水平集与流体体积 (Coupled Level Set and Vo-lume of Fluid,CLSVOF)方法[12].然而，VOF方法与LS方法同属于欧拉方法，相较于拉格朗日方法，在数值求解中具有相似的局限性[13].

为更加精确地描述运动界面，Enright创造性地提出了粒子水平集(Particle Level Set,PLS)方法[14].基本思想是根据正负逃逸粒子的位置信息通过粒子校正机制修正界面，体积守恒性和计算精度得到极大提高.但在后续的研究和应用中，也发现了该方法的不足.首先，过量粒子的引入导致计算成本大幅增加，难以应用于大规模的数值模拟[15].其次，正负逃逸粒子本应该用来修正其两侧界面，但PLS仅修正一侧界面，粒子校正机制并不合理[9].此外，PLS方法的计算精度对所采用的粒子重分配策略比较敏感[9,14].

针对上述PLS方法的缺陷，本文提出了单层粒子水平集(One-Layer Particle Level Set,OPLS)方法，引入单层粒子来追踪界面特征，计算效率得到较大提升.考虑到拉格朗日本质的精确性，粒子将会被流场输送到界面的精确位置，因此，界面可以直接根据单层粒子的位置进行校正，粒子校正机制更加合理.另外，本文提出了更加有效的粒子重分配策略，包括粒子的添加和删除，即使在处理复杂的拓扑变化时，也可以表现出良好的稳健性.最后，光滑的运动界面可由修正后的LS函数表示.在OPLS方法中，界面直接根据粒子的位置进行校正，合理的粒子校正机制使得所需的粒子数量大大降低，在不影响计算精度的前提下，极大地提高了计算效率.

1 粒子水平集方法基本理论

1.1 水平集方法

在LS方法中，通过定义符号距离函数φ隐式捕捉界面，符号距离函数的对流过程满足Hamilton-Jacobi方程：

(1)

式中：u表示流速；t表示时间.在计算过程中，LS函数可能丧失符号距离特性，为保持符号距离特性，文献[7,10]中提出了重新初始化过程：

(2)

式中：τ表示虚拟时间；n表示时间步；sgn()表示符号函数.

由于控制方程的离散和重新初始化过程引起的数值耗散会导致界面失真，为提高界面的精度，Enright等[14]创造性地提出了PLS方法.

1.2 粒子水平集方法

PLS方法引入了两类无质量拉格朗日粒子.在初始时刻，在距离界面3倍最大网格尺寸的单元内沿每个坐标方向布置4个粒子(二维单元内布置16个粒子，三维单元内布置64个粒子)，其中，正粒子布置在φ>0的区域，负粒子布置在界面另一侧.随后，将粒子吸引到界面上并随流运动.

在计算过程中，部分粒子可能因界面变形而穿过界面，这部分粒子被定义为逃逸粒子.如图1所示，PLS方法的粒子校正机制就是基于逃逸粒子而建立的，根据逃逸正粒子修正φ>0区域，另一侧界面根据逃逸负粒子位置进行重建.

图1 PLS方法示意图Fig.1 Schematic diagram of PLS method

随着界面的运动，部分区域可能缺少足够的粒子追踪界面信息，也可能存在拥有过多粒子的其他区域，为保证PLS方法的精度和计算效率，提出了粒子重分配策略.粒子的添加程序比较简单，识别距离界面3倍最大网格尺寸且缺少粒子的单元后，添加即可.粒子的删除程序相对比较复杂.由于逃逸粒子能够表征界面信息，所以需要全部保留.对于未逃逸粒子，如果距离界面较远可直接删除，若距离界面在3倍网格尺寸之内，则需要建立一个堆栈，比较粒子信息后删除多余粒子，控制单元内粒子数不超过默认数量.

1.3 单层粒子水平集

1.3.1基本思想 尽管PLS方法很大程度上提高了界面精度和体积守恒性，但仍然存在局限性：计算成本较高，粒子校正机制不合理，粒子重分配策略太繁琐.为克服上述缺陷，本文提出了OPLS方法.该方法利用LS方法隐式捕捉界面，之后对界面进行校正.全新的校正机制如图2所示，考虑到拉格朗日本质的精确性，界面直接根据单层粒子的位置信息进行校正，计算效率大大提升.总体而言，OPLS方法是采用LS方法在拉格朗日粒子的辅助下捕捉光滑、精确的界面.此外，本文提出了更加简洁有效的粒子重分配策略.

图2 OPLS方法示意图Fig.2 Schematic diagram of OPLS method

单层粒子LS方法的计算步骤如下：

步骤1生成拉格朗日粒子并且将其吸引到界面上;

步骤2更新粒子位置信息;

步骤3求解对流方程和重新初始化方程更新界面;

步骤4根据粒子位置修正界面;

步骤5根据修正后的界面添加或删除粒子.

1.3.2粒子初始布置 在OPLS方法中，界面穿过的单元被定义为界面单元.在初始时刻，粒子将被布置到界面单元中，如图3所示.首先，根据网格节点的LS函数值检测界面单元，然后沿界面单元内的局部坐标生成粒子.如图3(a)、3(c)、3(e)所示，沿单元局部坐标每个方向分别布置了1，2，3个粒子，即每个单元内布置1，4，9个粒子.粒子数可以根据所需要的计算精度进行选择，但过多的粒子不会对计算精度有明显的提升，反而会影响计算效率，建议采用每个单元内4个粒子进行数值模拟.

为精确地追踪界面，需要将生成的粒子吸引到界面上，如图3(b)、3(d)、3(f)所示.与PLS方法类似，本文将粒子沿最短路径吸引到界面上.由于所有的粒子均位于界面附近，所以其法向可认为是最短路径的方向.粒子的吸引过程如下：

xnew=xP+λ(φgoal-φP)n(xP)

(3)

式中：xnew表示粒子被吸引后的位置；xP表示粒子的当前位置；λ表示吸引参数，一般为1，但当粒子被移动距离过远穿过界面时可以减半；φgoal表示理想LS函数值，一般为0；φP表示当前位置的LS函数值，φgoal设为0以保证粒子最终被吸引到界面上；n(xP)表示当前位置的法向.此外，由于粒子处的法向可能不精确，式(3)并不能保证所有的粒子都被吸引，需要进行几次迭代，并且将始终无法被吸引的粒子直接删除.

粒子被吸引到界面上之后随流运动，可根据下式更新粒子位置：

(4)

式中：x表示粒子位置坐标；u(xp)表示粒子位置xp处的流速.

图3 粒子的初始布置Fig.3 Initial allocation of particles

1.3.3界面修正 界面修正是捕捉精确界面的关键.OPLS方法直接根据粒子位置对界面进行修正，失真界面可由粒子修正机制重建.

根据式(4)可计算得到粒子位置，进而可通过如下的插值过程得到粒子处LS函数值：

(5)

(6)

考虑到拉格朗日本质的精确性，粒子将会随流运动到界面的精确位置，即粒子处的LS函数值应为0.但由于LS方法对数值耗散极度敏感，根据式(5)计算得到的φP可能不为0，需要对粒子附近界面进行修正保证粒子处LS函数值为0.界面修正过程是将粒子处LS函数值的误差迭代分配到所在单元的网格节点上：

(7)

(8)

图4 界面修正Fig.4 Interface correction

粒子修正机制可以有效降低LS方法由于数值耗散而造成的误差，提升界面精度.相较于PLS方法，本文方法更加合理，并且计算成本大大降低.

1.3.4粒子重分配 当界面变化剧烈时，可能会存在某些缺少粒子的界面单元.充足的粒子是保证界面精度的前提，因此需要将粒子添加至这些单元.如图5所示(每个单元内所需粒子数为1)，粒子的添加程序分为两步：首先，根据LS函数值和单元内粒子数检索缺少粒子的界面单元；然后，在这些单元内生成粒子并吸引到界面上.

图5 粒子的添加Fig.5 Addition of particles

当界面合并时，部分界面会消失，但原有的粒子会留在当前区域，需要将其删除.如图6所示，粒子的删除过程分为两步：首先，沿粒子处法线方向距其一定距离处生成虚拟粒子I1，同时在反方向相同距离处生成虚拟粒子I2；然后，比较虚拟粒子处的LS函数值符号，若二者异号，则粒子仍位于界面上，需要保留，反之则将其删除.相较于PLS方法，本文所采用的粒子重分配策略更加简洁有效.值得注意的是，粒子的添加和删除并不需要每个时间步都进行，若界面不经历剧烈变化，甚至在整个计算过程内都不必执行.

图6 粒子的删除Fig.6 Deletion of particles

2 算例验证

本文采用了4个基准算例来验证OPLS方法的精确性，包括Zalesak圆盘的转动、圆在剪切流场中的流动以及界面的合并和分离.为定量描述体积守恒性，定义体积守恒相对误差为

(9)

H(x,y,z,t)=

(10)

式中：Ω表示积分区域；x、y、z表示积分点的坐标；H表示双曲正切函数；ε表示界面宽度的一半.

2.1 Zalesak圆盘的转动过程

Zalesak圆盘经常被用来测试数值算法对于界面尖角的描述能力.整个计算域是边长为1 m的正方形区域，坐标原点为左下角点.在初始时刻，半径为0.15 m的缺口圆盘的圆心位于(0.50,0.75)m.其中，圆盘缺口长度为0.25 m，宽度为0.05 m.恒定不可压流场定义为

(11)

式中：ux,y、vx,y分别表示点(x,y)处水平流速、垂向流速.在给定流场作用下，圆盘在1 s内将旋转 1圈.

图7 Zalesak圆盘转动过程图Fig.7 Process of the rotation of Zalesak disk

图8 Zalesak圆盘转动产生的体积损失Fig.8 Volume loss due to rotation of Zalesak disk

本次模拟中，圆盘旋转2圈，采用3套网格来检验本文方法的网格敏感性，网格尺寸(边长，下同)Δx分别为0.02、0.01 及0.005 m.计算结果如图7所示，采用OPLS方法可以精确地捕捉运动界面，并且能够清晰地描述界面的几何特性.而LS方法即使采用最细的Δx=0.005 m的网格也无法得到精确的界面，界面已经完全偏离原来的形状.为定量说明本文方法的精确性，给出了体积守恒相对误差的时间曲线，如图8所示.其中，LS方法产生的误差在持续增加，圆盘旋转一圈时误差为2.56%，在计算结束时增加至6.09%.而本文方法体现出良好的体积守恒性，即使是最粗的Δx=0.02 m的网格，两个时刻的误差仅为0.29%与0.39%.通过加密网格，有效降低了守恒性误差，其中，网格尺寸为0.01 m时，圆盘旋转1周和2周的相对误差分别为0.09%与0.10%，而网格尺寸为0.005 m时对应的两个时刻的相对误差分别为0.03%与0.04%.在表1中，将LS方法、PLS方法[16]的计算结果与本文方法进行了对比.可以看出，OPLS方法在粒子数减少的情况下，依然能够清晰地描述界面的几何特征，极大地改善了界面精度和体积守恒性.

表1 不同方法的质量损失相对误差Tab.1 Relative error of global conservation of different methods

2.2 圆在剪切流场中的流动

圆形界面在剪切流场中会被拉伸，产生一些难以捕捉的细小界面，经常被用来检验数值方法精确解析细微结构的能力.计算域是边长为1 m的正方形区域，坐标原点为左下角点半径为0.15 m的圆形界面在初始时刻圆心位于(0.5,0.75)m，然后在剪切流场的作用下开始运动，流场定义为

(12)

式中：T表示剪切流场的周期.而界面在t=T/2时变形尺度最大，本次模拟取T=8 s.

图9所示为不同时刻的粒子分布图，采用Δx=0.02 m的网格模拟的界面头部与尾部有细微的变形，尤其当界面变形程度最大时，这主要是由于网格分辨率较低造成的，图9(b)和9(c)的界面精度随网格的加密而明显提高.此外，本文方法在所有的网格分辨率下都体现出良好的体积守恒性.图10中对比了本文方法和LS方法的结果，可以明显看出LS方法模拟的界面变形严重，并且最后无法复原.图11所示为体积守恒相对误差的时间曲线，LS方法产生的误差持续增加，在计算结束时为34.58%.表2分别给出了LS方法、PLS方法[16]与本文方法在网格数为256×256时的计算结果，可以看出，PLS方法与本文方法在粒子的协助下能够保证较高的计算精度.同时，对比了不同粒子数对计算精度的影响，图10表明ncell为4与9的结果几乎一致.当单元内粒子数为1时，尽管在t=T时，误差为0.10%，但在t=T/2时，误差高达3.87%，这说明缺少足够的粒子来追踪界面，而ncell=4与ncell=9时产生的误差随时间变化不大，在计算结束时分别为0.09%与0.03%.因此在同时考虑计算成本和计算精度的情况下，建议在数值模拟中设置ncell=4.通过圆在剪切流场中流动的模拟，证明OPLS方法能够处理复杂的界面变形，解析细小结构，并保持质量守恒.

图9 不同网格尺寸下粒子的分布Fig.9 Particle distributions at different grids

图10 圆在剪切流场中不同时刻的形状Fig.10 Shape of circle in shear flow field at different times

图11 圆在剪切流场流动产生的体积损失Fig.11 Volume loss caused by flow of circle in shear flow field

表2 不同方法在特定时刻的相对质量损失误差Tab.2 Relative global mass conservation error of different methods at a specific time

2.3 界面合并

通过两个圆形界面的合并来证明本文方法能够精确处理界面的合并问题.初始条件如图12所示，计算域为边长为1 m的正方形区域，坐标原点为左下角点，两个半径为0.30 m的圆形界面的圆心分别位于(0.35,0.50)m与(0.65,0.50)m，然后在流场作用下逐渐合并.流场定义为

(13)

上述流场为不可压流场，不会产生额外的体积损失.

图12 界面合并的示意图(m)Fig.12 Schematic for interface merging (m)

图13 界面合并的过程图Fig.13 Process of interface merging

图14 界面合并过程中产生的体积损失Fig.14 Volume loss in interface merging

图13所示为LS方法和本文方法的模拟结果.在流场作用下，两个圆形界面在x方向上被压缩而在y方向上被拉伸.当界面合并时，部分界面会消失而粒子保留在原处，这时粒子删除程序会将其删除.同时，当界面拉伸至其他区域时，会缺少足够的粒子来追踪界面信息，粒子添加程序会随之生成粒子.因此，OPLS方法能够有效处理界面的合并问题，并能保持界面的几何特征.而LS方法尽管能够处理界面的合并，但无法保证体积守恒，误差较大.如图14所示，本文方法的体积守恒相对误差仅为0.58%，而LS方法的误差高达5.02%.因此，本文方法能够在保证体积守恒的前提下精确处理界面合并的问题.

2.4 界面分离

(14)

图16所示为LS方法和本文方法的模拟结果，可以看出，LS方法无法精确处理界面分离问题，并且无法准确传递界面的几何特征.根据六角形的边长和流速，可以确定六角形会在t=1.0 s时完全分离，而LS方法所捕捉的界面甚至在t=3.0 s还没有分离，界面也已经发生严重的变形.OPLS方法在粒子的协助下，精确传递了界面的几何特征，在t=1.0 s时界面完全分离为上下两部分.如图17所示，LS方法在计算过程中持续产生体积损失，最终误差高达10.38%.而本文方法依然体现出良好的守恒性，ncell=4与ncell=9的模拟结果几乎一致，误差分别为0.46%与0.41%.因此，建议采用ncell=4来进行计算模拟.得益于精确的粒子校正机制，OPLS方法可以精确处理界面分离的问题，体现出极佳的精确性和体积守恒性.

3 结语

通过引入拉格朗日粒子协助水平集方法捕捉界面，提出了一种单层粒子水平集方法，能够精确、光滑地表示运动界面，并且体现出良好的体积守恒性.首先，采用水平集方法隐式捕捉界面.然后，通过粒子校正机制利用单层粒子高度精确的位置信息直接修正界面.由于所需粒子数量大大减少，故该方法极大地提升了计算效率.此外，提出了一种简洁有效的粒子重分配策略，即使在处理界面复杂的拓扑变化时也体现出良好的稳健性.

采用4个基准算例验证了单层粒子水平集方法的精确性.首先，通过模拟Zalesak圆盘的转动过程验证了该方法的体积守恒性以及对于界面尖角的描述能力.然后，通过模拟圆在剪切流场中的流动，测试了该方法处理界面复杂变形以及解析细微结构的能力.最后，通过模拟界面合并和界面分离问题验证了本文方法能够处理界面复杂的拓扑变化.