牟凱龍，趙蘭浩，毛 佳

(河海大學(xué) 水利水電學(xué)院，南京 210098)

包含運(yùn)動(dòng)界面的流動(dòng)問(wèn)題廣泛存在于自然界中[1].對(duì)這類問(wèn)題進(jìn)行數(shù)值模擬的難點(diǎn)往往在于運(yùn)動(dòng)界面的描述，尤其是當(dāng)運(yùn)動(dòng)界面經(jīng)歷劇烈的拓?fù)渥兓瘯r(shí)，問(wèn)題求解難度會(huì)進(jìn)一步增加.因此，如何精確地表示運(yùn)動(dòng)界面一直是研究熱點(diǎn).

各國(guó)學(xué)者已提出多種數(shù)值方法來(lái)表示運(yùn)動(dòng)界面，大致可以分為拉格朗日方法和歐拉方法兩類.拉格朗日網(wǎng)格類方法的網(wǎng)格與流體域完全重合并跟隨流體運(yùn)動(dòng)，可以直接得到精確的運(yùn)動(dòng)界面.而且拉格朗日求解格式不存在對(duì)流項(xiàng)，因此這類方法具有良好的體積守恒性.但當(dāng)界面變形嚴(yán)重時(shí)，畸形的背景網(wǎng)格會(huì)導(dǎo)致計(jì)算精度下降，而網(wǎng)格重生成又會(huì)消耗大量的計(jì)算資源.為克服這種缺陷，有學(xué)者提出了拉格朗日粒子類方法，如光滑粒子質(zhì)點(diǎn)(Smoothed Particle Hydrodynamics,SPH)方法[2].這類方法采用粒子替代網(wǎng)格追蹤運(yùn)動(dòng)界面，盡管能夠精確追蹤任意變形的界面，但SPH方法固有的張力不穩(wěn)定性限制了該方法的應(yīng)用.此外，由于界面由離散的粒子表示，所以難以得到光滑的界面.

基于固定背景網(wǎng)格的歐拉方法根據(jù)描述界面的形式又可以分為界面追蹤方法和界面捕捉方法兩類[3].標(biāo)志網(wǎng)格(Marker and Cell,MAC)方法[4]是應(yīng)用最早的界面追蹤方法，主要思想是在背景網(wǎng)格中引入無(wú)質(zhì)量的拉格朗日粒子追蹤界面，最后根據(jù)粒子位置確定運(yùn)動(dòng)界面.因此，MAC方法也無(wú)法提供光滑的運(yùn)動(dòng)界面.此外，界面追蹤方法難以處理界面復(fù)雜的拓?fù)渥兓痆4].

另一類歐拉方法即界面捕捉方法通過(guò)設(shè)置指示函數(shù)來(lái)隱式地捕捉界面，能夠自動(dòng)處理復(fù)雜的拓?fù)渥兓瘧?yīng)用.其中，最廣泛的是流體體積(Volume of Fluid,VOF)方法[5]和水平集(Level Set,LS)方法[6-8].VOF方法通過(guò)在流場(chǎng)中定義流體體積函數(shù)，根據(jù)其函數(shù)值利用界面重構(gòu)算法構(gòu)造運(yùn)動(dòng)界面.雖然VOF方法憑借其良好的體積守恒性得到了廣泛的應(yīng)用，但是一方面需要復(fù)雜且計(jì)算成本較高的界面重構(gòu)技術(shù)，另一方面由于體積函數(shù)具有不連續(xù)性，不能準(zhǔn)確計(jì)算界面的法線方向和曲率.

相較于VOF方法，LS方法設(shè)置了一個(gè)連續(xù)的符號(hào)距離函數(shù)來(lái)隱式地捕捉界面，不僅能夠自動(dòng)處理復(fù)雜的拓?fù)渥兓夷軌蛑苯颖硎竟饣缑?，界面的幾何特征也能夠直接?jì)算.但LS方法對(duì)數(shù)值耗散極度敏感，這也導(dǎo)致其體積守恒性較差[9].在LS方法中，數(shù)值耗散主要來(lái)源于兩方面：控制方程的離散和重新初始化過(guò)程.前者在數(shù)值求解中難以避免，可通過(guò)采用高階離散格式來(lái)降低其影響.重新初始化過(guò)程是為保持LS函數(shù)的符號(hào)距離特性而引入的[10]，但界面會(huì)發(fā)生偏移，引起體積損失.針對(duì)LS方法體積守恒性差的缺陷，各國(guó)學(xué)者已提出多種改進(jìn)的方法.守恒式水平集(Conservative Level Set,CLS)方法[11]將符號(hào)距離函數(shù)替換為Heaviside函數(shù)，改善了其體積守恒性，但重新初始化過(guò)程仍然會(huì)導(dǎo)致體積損失.考慮到VOF方法和LS方法具有互補(bǔ)的優(yōu)缺點(diǎn)，有學(xué)者將二者結(jié)合提出了耦合式水平集與流體體積 (Coupled Level Set and Vo-lume of Fluid,CLSVOF)方法[12].然而，VOF方法與LS方法同屬于歐拉方法，相較于拉格朗日方法，在數(shù)值求解中具有相似的局限性[13].

為更加精確地描述運(yùn)動(dòng)界面，Enright創(chuàng)造性地提出了粒子水平集(Particle Level Set,PLS)方法[14].基本思想是根據(jù)正負(fù)逃逸粒子的位置信息通過(guò)粒子校正機(jī)制修正界面，體積守恒性和計(jì)算精度得到極大提高.但在后續(xù)的研究和應(yīng)用中，也發(fā)現(xiàn)了該方法的不足.首先，過(guò)量粒子的引入導(dǎo)致計(jì)算成本大幅增加，難以應(yīng)用于大規(guī)模的數(shù)值模擬[15].其次，正負(fù)逃逸粒子本應(yīng)該用來(lái)修正其兩側(cè)界面，但PLS僅修正一側(cè)界面，粒子校正機(jī)制并不合理[9].此外，PLS方法的計(jì)算精度對(duì)所采用的粒子重分配策略比較敏感[9,14].

針對(duì)上述PLS方法的缺陷，本文提出了單層粒子水平集(One-Layer Particle Level Set,OPLS)方法，引入單層粒子來(lái)追蹤界面特征，計(jì)算效率得到較大提升.考慮到拉格朗日本質(zhì)的精確性，粒子將會(huì)被流場(chǎng)輸送到界面的精確位置，因此，界面可以直接根據(jù)單層粒子的位置進(jìn)行校正，粒子校正機(jī)制更加合理.另外，本文提出了更加有效的粒子重分配策略，包括粒子的添加和刪除，即使在處理復(fù)雜的拓?fù)渥兓瘯r(shí)，也可以表現(xiàn)出良好的穩(wěn)健性.最后，光滑的運(yùn)動(dòng)界面可由修正后的LS函數(shù)表示.在OPLS方法中，界面直接根據(jù)粒子的位置進(jìn)行校正，合理的粒子校正機(jī)制使得所需的粒子數(shù)量大大降低，在不影響計(jì)算精度的前提下，極大地提高了計(jì)算效率.

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

1.1 水平集方法

在LS方法中，通過(guò)定義符號(hào)距離函數(shù)φ隱式捕捉界面，符號(hào)距離函數(shù)的對(duì)流過(guò)程滿足Hamilton-Jacobi方程：

(1)

式中：u表示流速；t表示時(shí)間.在計(jì)算過(guò)程中，LS函數(shù)可能喪失符號(hào)距離特性，為保持符號(hào)距離特性，文獻(xiàn)[7,10]中提出了重新初始化過(guò)程：

(2)

式中：τ表示虛擬時(shí)間；n表示時(shí)間步；sgn()表示符號(hào)函數(shù).

由于控制方程的離散和重新初始化過(guò)程引起的數(shù)值耗散會(huì)導(dǎo)致界面失真，為提高界面的精度，Enright等[14]創(chuàng)造性地提出了PLS方法.

1.2 粒子水平集方法

PLS方法引入了兩類無(wú)質(zhì)量拉格朗日粒子.在初始時(shí)刻，在距離界面3倍最大網(wǎng)格尺寸的單元內(nèi)沿每個(gè)坐標(biāo)方向布置4個(gè)粒子(二維單元內(nèi)布置16個(gè)粒子，三維單元內(nèi)布置64個(gè)粒子)，其中，正粒子布置在φ>0的區(qū)域，負(fù)粒子布置在界面另一側(cè).隨后，將粒子吸引到界面上并隨流運(yùn)動(dòng).

在計(jì)算過(guò)程中，部分粒子可能因界面變形而穿過(guò)界面，這部分粒子被定義為逃逸粒子.如圖1所示，PLS方法的粒子校正機(jī)制就是基于逃逸粒子而建立的，根據(jù)逃逸正粒子修正φ>0區(qū)域，另一側(cè)界面根據(jù)逃逸負(fù)粒子位置進(jìn)行重建.

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

隨著界面的運(yùn)動(dòng)，部分區(qū)域可能缺少足夠的粒子追蹤界面信息，也可能存在擁有過(guò)多粒子的其他區(qū)域，為保證PLS方法的精度和計(jì)算效率，提出了粒子重分配策略.粒子的添加程序比較簡(jiǎn)單，識(shí)別距離界面3倍最大網(wǎng)格尺寸且缺少粒子的單元后，添加即可.粒子的刪除程序相對(duì)比較復(fù)雜.由于逃逸粒子能夠表征界面信息，所以需要全部保留.對(duì)于未逃逸粒子，如果距離界面較遠(yuǎn)可直接刪除，若距離界面在3倍網(wǎng)格尺寸之內(nèi)，則需要建立一個(gè)堆棧，比較粒子信息后刪除多余粒子，控制單元內(nèi)粒子數(shù)不超過(guò)默認(rèn)數(shù)量.

1.3 單層粒子水平集

1.3.1基本思想 盡管PLS方法很大程度上提高了界面精度和體積守恒性，但仍然存在局限性：計(jì)算成本較高，粒子校正機(jī)制不合理，粒子重分配策略太繁瑣.為克服上述缺陷，本文提出了OPLS方法.該方法利用LS方法隱式捕捉界面，之后對(duì)界面進(jìn)行校正.全新的校正機(jī)制如圖2所示，考慮到拉格朗日本質(zhì)的精確性，界面直接根據(jù)單層粒子的位置信息進(jìn)行校正，計(jì)算效率大大提升.總體而言，OPLS方法是采用LS方法在拉格朗日粒子的輔助下捕捉光滑、精確的界面.此外，本文提出了更加簡(jiǎn)潔有效的粒子重分配策略.

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

單層粒子LS方法的計(jì)算步驟如下：

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

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

步驟3求解對(duì)流方程和重新初始化方程更新界面;

步驟4根據(jù)粒子位置修正界面;

步驟5根據(jù)修正后的界面添加或刪除粒子.

1.3.2粒子初始布置 在OPLS方法中，界面穿過(guò)的單元被定義為界面單元.在初始時(shí)刻，粒子將被布置到界面單元中，如圖3所示.首先，根據(jù)網(wǎng)格節(jié)點(diǎn)的LS函數(shù)值檢測(cè)界面單元，然后沿界面單元內(nèi)的局部坐標(biāo)生成粒子.如圖3(a)、3(c)、3(e)所示，沿單元局部坐標(biāo)每個(gè)方向分別布置了1，2，3個(gè)粒子，即每個(gè)單元內(nèi)布置1，4，9個(gè)粒子.粒子數(shù)可以根據(jù)所需要的計(jì)算精度進(jìn)行選擇，但過(guò)多的粒子不會(huì)對(duì)計(jì)算精度有明顯的提升，反而會(huì)影響計(jì)算效率，建議采用每個(gè)單元內(nèi)4個(gè)粒子進(jìn)行數(shù)值模擬.

為精確地追蹤界面，需要將生成的粒子吸引到界面上，如圖3(b)、3(d)、3(f)所示.與PLS方法類似，本文將粒子沿最短路徑吸引到界面上.由于所有的粒子均位于界面附近，所以其法向可認(rèn)為是最短路徑的方向.粒子的吸引過(guò)程如下：

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

(3)

式中：xnew表示粒子被吸引后的位置；xP表示粒子的當(dāng)前位置；λ表示吸引參數(shù)，一般為1，但當(dāng)粒子被移動(dòng)距離過(guò)遠(yuǎn)穿過(guò)界面時(shí)可以減半；φgoal表示理想LS函數(shù)值，一般為0；φP表示當(dāng)前位置的LS函數(shù)值，φgoal設(shè)為0以保證粒子最終被吸引到界面上；n(xP)表示當(dāng)前位置的法向.此外，由于粒子處的法向可能不精確，式(3)并不能保證所有的粒子都被吸引，需要進(jìn)行幾次迭代，并且將始終無(wú)法被吸引的粒子直接刪除.

粒子被吸引到界面上之后隨流運(yùn)動(dòng)，可根據(jù)下式更新粒子位置：

(4)

式中：x表示粒子位置坐標(biāo)；u(xp)表示粒子位置xp處的流速.

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

1.3.3界面修正 界面修正是捕捉精確界面的關(guān)鍵.OPLS方法直接根據(jù)粒子位置對(duì)界面進(jìn)行修正，失真界面可由粒子修正機(jī)制重建.

根據(jù)式(4)可計(jì)算得到粒子位置，進(jìn)而可通過(guò)如下的插值過(guò)程得到粒子處LS函數(shù)值：

(5)

(6)

考慮到拉格朗日本質(zhì)的精確性，粒子將會(huì)隨流運(yùn)動(dòng)到界面的精確位置，即粒子處的LS函數(shù)值應(yīng)為0.但由于LS方法對(duì)數(shù)值耗散極度敏感，根據(jù)式(5)計(jì)算得到的φP可能不為0，需要對(duì)粒子附近界面進(jìn)行修正保證粒子處LS函數(shù)值為0.界面修正過(guò)程是將粒子處LS函數(shù)值的誤差迭代分配到所在單元的網(wǎng)格節(jié)點(diǎn)上：

(7)

(8)

圖4 界面修正Fig.4 Interface correction

粒子修正機(jī)制可以有效降低LS方法由于數(shù)值耗散而造成的誤差，提升界面精度.相較于PLS方法，本文方法更加合理，并且計(jì)算成本大大降低.

1.3.4粒子重分配 當(dāng)界面變化劇烈時(shí)，可能會(huì)存在某些缺少粒子的界面單元.充足的粒子是保證界面精度的前提，因此需要將粒子添加至這些單元.如圖5所示(每個(gè)單元內(nèi)所需粒子數(shù)為1)，粒子的添加程序分為兩步：首先，根據(jù)LS函數(shù)值和單元內(nèi)粒子數(shù)檢索缺少粒子的界面單元；然后，在這些單元內(nèi)生成粒子并吸引到界面上.

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

當(dāng)界面合并時(shí)，部分界面會(huì)消失，但原有的粒子會(huì)留在當(dāng)前區(qū)域，需要將其刪除.如圖6所示，粒子的刪除過(guò)程分為兩步：首先，沿粒子處法線方向距其一定距離處生成虛擬粒子I1，同時(shí)在反方向相同距離處生成虛擬粒子I2；然后，比較虛擬粒子處的LS函數(shù)值符號(hào)，若二者異號(hào)，則粒子仍位于界面上，需要保留，反之則將其刪除.相較于PLS方法，本文所采用的粒子重分配策略更加簡(jiǎn)潔有效.值得注意的是，粒子的添加和刪除并不需要每個(gè)時(shí)間步都進(jìn)行，若界面不經(jīng)歷劇烈變化，甚至在整個(gè)計(jì)算過(guò)程內(nèi)都不必執(zhí)行.

圖6 粒子的刪除Fig.6 Deletion of particles

2 算例驗(yàn)證

本文采用了4個(gè)基準(zhǔn)算例來(lái)驗(yàn)證OPLS方法的精確性，包括Zalesak圓盤的轉(zhuǎn)動(dòng)、圓在剪切流場(chǎng)中的流動(dòng)以及界面的合并和分離.為定量描述體積守恒性，定義體積守恒相對(duì)誤差為

(9)

H(x,y,z,t)=

(10)

式中：Ω表示積分區(qū)域；x、y、z表示積分點(diǎn)的坐標(biāo)；H表示雙曲正切函數(shù)；ε表示界面寬度的一半.

2.1 Zalesak圓盤的轉(zhuǎn)動(dòng)過(guò)程

Zalesak圓盤經(jīng)常被用來(lái)測(cè)試數(shù)值算法對(duì)于界面尖角的描述能力.整個(gè)計(jì)算域是邊長(zhǎng)為1 m的正方形區(qū)域，坐標(biāo)原點(diǎn)為左下角點(diǎn).在初始時(shí)刻，半徑為0.15 m的缺口圓盤的圓心位于(0.50,0.75)m.其中，圓盤缺口長(zhǎng)度為0.25 m，寬度為0.05 m.恒定不可壓流場(chǎng)定義為

(11)

式中：ux,y、vx,y分別表示點(diǎn)(x,y)處水平流速、垂向流速.在給定流場(chǎng)作用下，圓盤在1 s內(nèi)將旋轉(zhuǎn) 1圈.

圖7 Zalesak圓盤轉(zhuǎn)動(dòng)過(guò)程圖Fig.7 Process of the rotation of Zalesak disk

圖8 Zalesak圓盤轉(zhuǎn)動(dòng)產(chǎn)生的體積損失Fig.8 Volume loss due to rotation of Zalesak disk

本次模擬中，圓盤旋轉(zhuǎn)2圈，采用3套網(wǎng)格來(lái)檢驗(yàn)本文方法的網(wǎng)格敏感性，網(wǎng)格尺寸(邊長(zhǎng)，下同)Δx分別為0.02、0.01 及0.005 m.計(jì)算結(jié)果如圖7所示，采用OPLS方法可以精確地捕捉運(yùn)動(dòng)界面，并且能夠清晰地描述界面的幾何特性.而LS方法即使采用最細(xì)的Δx=0.005 m的網(wǎng)格也無(wú)法得到精確的界面，界面已經(jīng)完全偏離原來(lái)的形狀.為定量說(shuō)明本文方法的精確性，給出了體積守恒相對(duì)誤差的時(shí)間曲線，如圖8所示.其中，LS方法產(chǎn)生的誤差在持續(xù)增加，圓盤旋轉(zhuǎn)一圈時(shí)誤差為2.56%，在計(jì)算結(jié)束時(shí)增加至6.09%.而本文方法體現(xiàn)出良好的體積守恒性，即使是最粗的Δx=0.02 m的網(wǎng)格，兩個(gè)時(shí)刻的誤差僅為0.29%與0.39%.通過(guò)加密網(wǎng)格，有效降低了守恒性誤差，其中，網(wǎng)格尺寸為0.01 m時(shí)，圓盤旋轉(zhuǎn)1周和2周的相對(duì)誤差分別為0.09%與0.10%，而網(wǎng)格尺寸為0.005 m時(shí)對(duì)應(yīng)的兩個(gè)時(shí)刻的相對(duì)誤差分別為0.03%與0.04%.在表1中，將LS方法、PLS方法[16]的計(jì)算結(jié)果與本文方法進(jìn)行了對(duì)比.可以看出，OPLS方法在粒子數(shù)減少的情況下，依然能夠清晰地描述界面的幾何特征，極大地改善了界面精度和體積守恒性.

表1 不同方法的質(zhì)量損失相對(duì)誤差Tab.1 Relative error of global conservation of different methods

2.2 圓在剪切流場(chǎng)中的流動(dòng)

圓形界面在剪切流場(chǎng)中會(huì)被拉伸，產(chǎn)生一些難以捕捉的細(xì)小界面，經(jīng)常被用來(lái)檢驗(yàn)數(shù)值方法精確解析細(xì)微結(jié)構(gòu)的能力.計(jì)算域是邊長(zhǎng)為1 m的正方形區(qū)域，坐標(biāo)原點(diǎn)為左下角點(diǎn)半徑為0.15 m的圓形界面在初始時(shí)刻圓心位于(0.5,0.75)m，然后在剪切流場(chǎng)的作用下開(kāi)始運(yùn)動(dòng)，流場(chǎng)定義為

(12)

式中：T表示剪切流場(chǎng)的周期.而界面在t=T/2時(shí)變形尺度最大，本次模擬取T=8 s.

圖9所示為不同時(shí)刻的粒子分布圖，采用Δx=0.02 m的網(wǎng)格模擬的界面頭部與尾部有細(xì)微的變形，尤其當(dāng)界面變形程度最大時(shí)，這主要是由于網(wǎng)格分辨率較低造成的，圖9(b)和9(c)的界面精度隨網(wǎng)格的加密而明顯提高.此外，本文方法在所有的網(wǎng)格分辨率下都體現(xiàn)出良好的體積守恒性.圖10中對(duì)比了本文方法和LS方法的結(jié)果，可以明顯看出LS方法模擬的界面變形嚴(yán)重，并且最后無(wú)法復(fù)原.圖11所示為體積守恒相對(duì)誤差的時(shí)間曲線，LS方法產(chǎn)生的誤差持續(xù)增加，在計(jì)算結(jié)束時(shí)為34.58%.表2分別給出了LS方法、PLS方法[16]與本文方法在網(wǎng)格數(shù)為256×256時(shí)的計(jì)算結(jié)果，可以看出，PLS方法與本文方法在粒子的協(xié)助下能夠保證較高的計(jì)算精度.同時(shí)，對(duì)比了不同粒子數(shù)對(duì)計(jì)算精度的影響，圖10表明ncell為4與9的結(jié)果幾乎一致.當(dāng)單元內(nèi)粒子數(shù)為1時(shí)，盡管在t=T時(shí)，誤差為0.10%，但在t=T/2時(shí)，誤差高達(dá)3.87%，這說(shuō)明缺少足夠的粒子來(lái)追蹤界面，而ncell=4與ncell=9時(shí)產(chǎn)生的誤差隨時(shí)間變化不大，在計(jì)算結(jié)束時(shí)分別為0.09%與0.03%.因此在同時(shí)考慮計(jì)算成本和計(jì)算精度的情況下，建議在數(shù)值模擬中設(shè)置ncell=4.通過(guò)圓在剪切流場(chǎng)中流動(dòng)的模擬，證明OPLS方法能夠處理復(fù)雜的界面變形，解析細(xì)小結(jié)構(gòu)，并保持質(zhì)量守恒.

圖9 不同網(wǎng)格尺寸下粒子的分布Fig.9 Particle distributions at different grids

圖10 圓在剪切流場(chǎng)中不同時(shí)刻的形狀Fig.10 Shape of circle in shear flow field at different times

圖11 圓在剪切流場(chǎng)流動(dòng)產(chǎn)生的體積損失Fig.11 Volume loss caused by flow of circle in shear flow field

表2 不同方法在特定時(shí)刻的相對(duì)質(zhì)量損失誤差Tab.2 Relative global mass conservation error of different methods at a specific time

2.3 界面合并

通過(guò)兩個(gè)圓形界面的合并來(lái)證明本文方法能夠精確處理界面的合并問(wèn)題.初始條件如圖12所示，計(jì)算域?yàn)檫呴L(zhǎng)為1 m的正方形區(qū)域，坐標(biāo)原點(diǎn)為左下角點(diǎn)，兩個(gè)半徑為0.30 m的圓形界面的圓心分別位于(0.35,0.50)m與(0.65,0.50)m，然后在流場(chǎng)作用下逐漸合并.流場(chǎng)定義為

(13)

上述流場(chǎng)為不可壓流場(chǎng)，不會(huì)產(chǎn)生額外的體積損失.

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

圖13 界面合并的過(guò)程圖Fig.13 Process of interface merging

圖14 界面合并過(guò)程中產(chǎn)生的體積損失Fig.14 Volume loss in interface merging

圖13所示為L(zhǎng)S方法和本文方法的模擬結(jié)果.在流場(chǎng)作用下，兩個(gè)圓形界面在x方向上被壓縮而在y方向上被拉伸.當(dāng)界面合并時(shí)，部分界面會(huì)消失而粒子保留在原處，這時(shí)粒子刪除程序會(huì)將其刪除.同時(shí)，當(dāng)界面拉伸至其他區(qū)域時(shí)，會(huì)缺少足夠的粒子來(lái)追蹤界面信息，粒子添加程序會(huì)隨之生成粒子.因此，OPLS方法能夠有效處理界面的合并問(wèn)題，并能保持界面的幾何特征.而LS方法盡管能夠處理界面的合并，但無(wú)法保證體積守恒，誤差較大.如圖14所示，本文方法的體積守恒相對(duì)誤差僅為0.58%，而LS方法的誤差高達(dá)5.02%.因此，本文方法能夠在保證體積守恒的前提下精確處理界面合并的問(wèn)題.

2.4 界面分離

(14)

圖16所示為L(zhǎng)S方法和本文方法的模擬結(jié)果，可以看出，LS方法無(wú)法精確處理界面分離問(wèn)題，并且無(wú)法準(zhǔn)確傳遞界面的幾何特征.根據(jù)六角形的邊長(zhǎng)和流速，可以確定六角形會(huì)在t=1.0 s時(shí)完全分離，而LS方法所捕捉的界面甚至在t=3.0 s還沒(méi)有分離，界面也已經(jīng)發(fā)生嚴(yán)重的變形.OPLS方法在粒子的協(xié)助下，精確傳遞了界面的幾何特征，在t=1.0 s時(shí)界面完全分離為上下兩部分.如圖17所示，LS方法在計(jì)算過(guò)程中持續(xù)產(chǎn)生體積損失，最終誤差高達(dá)10.38%.而本文方法依然體現(xiàn)出良好的守恒性，ncell=4與ncell=9的模擬結(jié)果幾乎一致，誤差分別為0.46%與0.41%.因此，建議采用ncell=4來(lái)進(jìn)行計(jì)算模擬.得益于精確的粒子校正機(jī)制，OPLS方法可以精確處理界面分離的問(wèn)題，體現(xiàn)出極佳的精確性和體積守恒性.

3 結(jié)語(yǔ)

通過(guò)引入拉格朗日粒子協(xié)助水平集方法捕捉界面，提出了一種單層粒子水平集方法，能夠精確、光滑地表示運(yùn)動(dòng)界面，并且體現(xiàn)出良好的體積守恒性.首先，采用水平集方法隱式捕捉界面.然后，通過(guò)粒子校正機(jī)制利用單層粒子高度精確的位置信息直接修正界面.由于所需粒子數(shù)量大大減少，故該方法極大地提升了計(jì)算效率.此外，提出了一種簡(jiǎn)潔有效的粒子重分配策略，即使在處理界面復(fù)雜的拓?fù)渥兓瘯r(shí)也體現(xiàn)出良好的穩(wěn)健性.

采用4個(gè)基準(zhǔn)算例驗(yàn)證了單層粒子水平集方法的精確性.首先，通過(guò)模擬Zalesak圓盤的轉(zhuǎn)動(dòng)過(guò)程驗(yàn)證了該方法的體積守恒性以及對(duì)于界面尖角的描述能力.然后，通過(guò)模擬圓在剪切流場(chǎng)中的流動(dòng)，測(cè)試了該方法處理界面復(fù)雜變形以及解析細(xì)微結(jié)構(gòu)的能力.最后，通過(guò)模擬界面合并和界面分離問(wèn)題驗(yàn)證了本文方法能夠處理界面復(fù)雜的拓?fù)渥兓?