單層粒子水平集方法

牟凱龍，趙蘭浩，毛 佳

(河海大學 水利水電學院，南京 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圓盤的轉動過程驗證了該方法的體積守恒性以及對于界面尖角的描述能力.然后，通過模擬圓在剪切流場中的流動，測試了該方法處理界面復雜變形以及解析細微結構的能力.最后，通過模擬界面合并和界面分離問題驗證了本文方法能夠處理界面復雜的拓撲變化.