黃筱云，李紹武

(天津大學建筑工程學院，天津 300072)

空間自適應的快速粒子Level Set方法

黃筱云，李紹武

(天津大學建筑工程學院，天津 300072)

Level set作為一種有效的界面追蹤方法，已經被廣泛地應用到單相流的自由表面追蹤或多相流的交界面追蹤.快速粒子level set方法(FPLS)是一種改進的level set方法，既具有較高的界面追蹤精度，又有較高的計算效率.將其與自適應網格結合起來，提出了一種基于樹型結構空間自適應的快速粒子 level set方法(SA-FPLS).通過在交界面附近設置高分辨率網格帶，且只在高分辨率網格帶內求解level set方程，既進一步提高了交界面模擬精度，又可以節(jié)省計算機內存，同時保持較高計算效率.

快速粒子level set；自適應網格；空間自適應的快速粒子level set

在用數值方法求解帶有交界面的流體力學問題時，必須通過一定方法預先給出流體交界面的位置.level set方法是近年發(fā)展起來的一種界面追蹤方法[1]，其原理是將交界面看作滿足以下控制方程的符號距離函數φ的零等值線(面).

式中：φ稱為level set函數，其值在交界面上為零，離開交界面為正或為負；ui為速度張量.

為使交界面的解在出現界面融合或較大φ值梯度時仍有較高的精度，一般采用重新初始化的方法[2]，即在每完成一步level set方程求解過程后，通過迭代求解方程

的解，來對level set函數進行初始化，其中0φ為level set方程的計算結果，xΔ為網格尺寸.這種level set方法在整個計算域上進行計算，稱之為全局 level set方法.

為了減少計算量，Peng等[3]提出了一種局部的level set方法，它只需在交界面附近區(qū)域求解level set方程和重新初始化.事實上，level set方法中采用迭代法進行重新初始化是影響計算速度的關鍵因素，Sethian[4]提出了一種fast marching技術求解level set方程，其特點是無需重新初始化，但計算精度較低.

與VOF方法相比較，level set方法的缺點是守恒性差.為了彌補這一缺陷，Sussman等[5]將 VOF與level set方法相結合，提出了混合型的界面追蹤方法——CLSVOF，該方法利用 VOF函數來校正 level set函數結果以達到減少質量損失的目的，但該方法并不能完全消除交界面附近的噪聲，追蹤到的交界面與實際值存在相位差.另外一種改進的level set方法是粒子level set(PLS)方法[6]，它將MAC法的特性引入到 level set方法中，計算時通過粒子來修正 level set函數，它捕捉到的交界面十分光滑，幾乎無質量損失.為了加快計算速度，Enright等[7]提出了一種局部快速粒子level set方法——FPLS，該方法采用一階精度的半Lagrangian格式[8]求解level set方程，用一階精度fast marching技術[9]進行重新初始化，粒子的運算則采用二階精度Runge-Kutta格式.

在實際應用中，復雜交界面需要高分辨率網格來模擬，采用均勻網格將消耗大量計算機內存，計算速度也將大大降低.Strain[10]在采用半 Lagrangian方法求解level set方程時，首先引入了樹型結構的自適應網格，這種網格結構在交界面附近采用高分辨率網格，而在遠離交界面的區(qū)域則保持較粗的分辨率.這樣一方面提高了計算結果的精度，另一方面減少了無關區(qū)域所占用的內存.但該方法為了保持在 T型節(jié)點插值的連續(xù)性，重新初始化階段在交界面附近區(qū)域采用三角網格[11]，計算量仍較大.2002年，Sochnikov等[12]將 fast marching技術應用到四叉樹型結構的自適應網格中，使得計算量進一步減少.這種自適應level set方法是通過減小交界面附近網格尺寸來保證解的守恒性，可能會導致交界面周圍的網格與最粗網格的尺寸相差很大，形成了過于復雜的樹型結構.

筆者將樹型結構的自適應性與 FPLS方法結合起來，提出一種改進的自適應 level set方法——SAFPLS，該方法在交界面附近能自動形成一個具有固定寬度的高分辨率網格帶，并在該網格帶內求解level set方程，進行重新初始化運算，且粒子始終包含在最小網格帶內.該方法保留了自適應level set方法的優(yōu)點，同時 FPLS方法的高精度性使得該方法無需在交界面周圍采用過細的網格，這樣在保證計算結果高精度的同時，計算速度更快，消耗內存更少.SAFPLS方法的具體步驟為：

(1)按照初始化 level set函數，生成在交界面附近加密的樹型結構網格；

(2) 采用FPLS方法求解level set方程；

(3) 根據level set函數值重新調整網格；

(4) 返回第(2)步.

1 FPLS方法的基本原理

在 FPLS方法中，正負粒子被隨機布置在交界面周圍窄帶中(初始時刻正值區(qū)的粒子為正粒子；反之為負粒子).每個粒子有各自的半徑 rp，大小被限定在 最 小 值 rmin= 0 .1max ( Δ x, Δy )和 最 大 值 rmax=0.5max(Δ x,Δy)之間，其值可根據下列公式確定.

式中sp為粒子的符號(正粒子為+1，負粒子為-1).粒子的φ值將通過雙線性或三線性插值得到.

這種確定半徑辦法的目的是使粒子盡量與交界面相切.當粒子穿過交界面達一個半徑時，即認為是逃逸粒子.在劇烈變化的流場中，通過 level set方程得到的交界面附近的 level set函數值可能存在較大的誤差.為了消除誤差，可借助逃逸粒子來對 level set函數值進行修正.為此，定義一個與逃逸粒子有關的局部level set函數pφ，即

如果函數pφ和φ值不相等，說明level set方程得出的φ函數可能存在誤差，取兩者中絕對值較小的作為校正后的φ值.這種對φ的校正在重新初始化的過程中同樣也需要.

當界面發(fā)生拉伸和斷裂時，交界面附近網格粒子變得稀疏，需要補充粒子以確保追蹤界面的準確性.而當非逃逸粒子遠離交界面時，則需要將這些粒子刪除.

在FPLS中，采用了一種低階的半Lagrangian算法求解level set方程，該格式無條件穩(wěn)定且收斂于真解.公式為

在進行重新初始化過程中，式(3)采用以下一階迎風差分格式進行離散.

通過求解式(7)，式(3)的解可以快速地構造出來，這就是fast marching技術的主要思想，即從界面附近區(qū)域出發(fā)，逐步構造出Eikonal方程的解(見圖1).

圖1 采用快速步進法進行重新初始化Fig.1 Reinitialization with fast marching technique

2 SA-FPLS的原理

為簡單起見，這里僅以二維的四叉樹網格說明SA-FPLS的原理.三維的八叉樹網格與之類似，不再贅述.在四叉樹網格結構中，每個四叉樹網格定義4個指向其子網格的指針(如果子網格不存在，指針指向空)和一個指向其父網格(如果父網格不存在，指針指向空)的指針，每個網格另有4個指向角點的指針和4個指向邊的指針.每個節(jié)點有4個指針指向其相鄰葉網格，對于特殊的節(jié)點，如T型節(jié)點，4個指針中有2個指向同一個網格.每個節(jié)點需要保存該節(jié)點的坐標及定義在該點的物理量.最粗的網格為 0層網格，略粗的為第1層，以此類推.圖2中網格Ⅰ有4個子網格Ⅱ、Ⅲ、Ⅳ和Ⅴ和 4個節(jié)點 A、B、C 和 D，與子網格Ⅳ有共同節(jié)點C.網格Ⅳ有 4個子網格Ⅵ、Ⅶ、Ⅷ和Ⅸ和1個父網格Ⅰ.網格Ⅱ與網格Ⅲ共享1號邊.節(jié)點E有4個指針分別指向網格Ⅱ、Ⅲ、Ⅳ和Ⅴ.在這種四叉樹網格結構中，網格在每個方位上都有唯一的鄰居，定義為與其相鄰且層數不大于該網格的網格中層數最大的一個.如圖3中網格Ⅱ與網格Ⅰ、Ⅳ、Ⅴ和Ⅵ相鄰，但只有網格Ⅳ才是網格Ⅱ的鄰居.

圖2 樹型網格結構Fig.2 Tree grid structure

圖3 網格及其鄰居Fig.3 Given grid and its neighboring grid

在生成樹型結構網格的過程中，首先分裂包含交界面的網格，若一個網格的角點上φ的符號相反，交界面必定穿越該網格，該網格即需要分裂，直至網格的層數達到規(guī)定的最大層數.因此，交界面附近的網格分辨率是最高的，那些與高分辨率網格相鄰的網格也需要進行分裂，以保證高分辨率網格具有一定帶寬.在交界面周圍生成高分辨率的網格之后，還需要對網格進行平衡處理[13]，即保持相鄰網格的大小相差不會超過 1倍.新產生的網格節(jié)點上的物理量通過線性插值得到.在計算過程中，為了提高交界面追蹤精度，交界面附近的網格應始終保持高分辨率，且加密網格隨著交界面的變化而變化.這樣可以保證只在網格帶內求解level set方程和進行φ值校正，同時防止因在交界面附近分裂網格時，線性插值導致的較大計算誤差.因此，隨著交界面的運動，需要在交界面附近將預先設定的高分辨率帶的網格不斷進行分裂，同時也要將脫離高分辨率網格帶的網格進行合并(見圖 4).

圖4 交界面附近的高分辨率網格帶Fig.4 Grid band of highest resolution around interface

在四叉樹網格中采用半 Lagrangian方法求解level set方程與在均勻網格中有所不同，網格節(jié)點在上個時刻投影點的位置需要根據四叉樹網格的結構尋找.首先，尋找投影點所在的 0層網格；然后，從該網格逐步向下找到投影點所在的子網格；最后，根據所在子網格角點上前一時刻的φ值，線性插值得到該節(jié)點新時刻的φ值.

在重新初始化的過程中，由于采用 fast marching技術，φ值是從交界面出發(fā)，迎風地向外構造得到，因此不會影響交界面的位置，當計算的某個節(jié)點的φ值大于高分辨率帶的寬度時，向外構造過程即停止.值得注意的是，在向外構造的過程中，若遇到 T型節(jié)點時，則需要分裂該點周圍所有層數小于規(guī)定最大層數的網格，以保證高精度網格帶的寬度.

3 算例及分析

在二維算例中，采用如下一階近似算法計算交界面的長度L，即

3.1 算例1

本算例模擬星型結構的擴張過程.交界面初始形狀如圖 5(a)所示，交界面以單位速度向向外法向擴張.計算區(qū)域大小為[0 ,2]×[0 ,2]，網格層數為4層，0層網格尺寸為0.04×0.04.初始level set函數表示為

式中：Ra為平均半徑；d為偏離Ra的距離；n為觸角的個數；θ為半徑與 x正向的夾角.在本算例中，Ra= 0 .3, d = 0 .2, n= 7 .圖 5(b)是交界面分別在 0～0.05 s 6個時刻的界面位置.可以看出，界面在運動中始終保持與初始形狀特征相似的形狀，只是觸角變得越來越厚.

圖5 七角星結構擴張過程的數值模擬結果Fig.5 Numerical results of expansion of seven-pointed star

3.2 算例2

本算例目的在于檢驗 SA-FPLS方法在剪切流場中對復雜交界面追蹤的有效性.計算區(qū)域為[0 ,1]×[0 ,1]，網格層數為4層，0層網格尺寸為0.04×0.04，界面初始形狀為一個半徑為0.15、圓心在[0 .5,0.75]處的圓(見圖6(a)).速度場表示為

式中，流場周期T為8 s，時間步長為0.005 s.圓的變形過程如圖 6所示.可見，由于在交界面附近采用了高分辨率的網格，尾部絲狀結構保持得相當完整.同時，由于采用自適應網格，大部分區(qū)域保持較粗的網格，從而節(jié)約了大量內存和計算時間.

3.3 算例3

本算例模擬圓在1個有16個渦旋的流場中的變形情況.圓的初始狀態(tài)及初始網格結構與算例 2一致，不同的是，該模型采用周期性邊界條件，其周期性速度場[15]表示為

式中，流場周期T為2 s，計算時間步長為0.005 s.變形過程如圖 7所示，在變形過程中，由于采用了周期性邊界條件，交界面穿過計算區(qū)域上端后，又重新出現在計算區(qū)域下端，且上下邊界的網格與節(jié)點保持一一對應.可以看出 FPLS與自適應網格相結合，可以適應十分復雜的流場結構.

圖6 圓在簡單旋轉流場中變形過程的數值模擬結果Fig.6 Numerical results of transformation of circle in single vortex flow

圖7 圓在多渦流場中變形情況的數值模擬Fig.7 Numerical results of transformation of circle in multi-vortex flow

檢驗模型守恒性優(yōu)劣的方法是對比初始時刻界面所圍圓面積與 1個周期后界面恢復到原始形狀后面積的差異.界面所圍面積的計算式為

式中：L為交界面的期望長度；H為Heaviside函數；φe為真解；φc為計算結果.算例2和算例3計算結果的誤差分析如表1所示.

3.4 算例4

算例4和算例5將討論SA-FPLS方法在三維空間中的應用.算例4模擬一個半徑為0.15的Zalesak圓球[16]在流場作用下的變化情況，球的初始位置為[0 .5,0.75,0.5]，切口的寬度為 0.05，深度為 0.125.計算區(qū)域尺度為[0 ,1] × [ 0 ,1] × [ 0 ,1]，網格層數為 4，最粗網格大小為0.0625× 0 .0625× 0 .0625，流場表示為

表1 計算結果的誤差分析Tab.1 Error analysis of numerical results

計算時間為 2 s，時間步長為 0.005 s.界面運動過程如圖 8所示，Zalesak圓球經過一個周期的運動后又回到初始位置，而缺口仍然保持尖銳.

圖8 Zalesak圓球的旋轉過程數值模擬結果Fig.8 Numerical results of rotation of Zalesak sphere

3.5 算例5

本算例模擬一個半徑為 0.15、初始位置為[0.35,0.35,0.35]的球體在三維流場作用下的變形過程，如圖 9所示.計算區(qū)域為[0 , 1] × [ 0 , 1] × [ 0 , 1]，網格層數為 3，最粗網格大小為0.04× 0 .04× 0 .04.流場表示為式中，流場周期 T為 2,s，計算時間步長為 0.005,s.結果表明，圓球在1.0,s時達到最大拉伸后，又恢復到原樣，而其中各圖為圓球分別在 t為0～2.0,s時的形狀.

4 結 語

將自適應網格與 FPLS結合起來，探討了一種交界面追蹤方法 SA-FPLS.從算例中可以看出，SAFPLS一方面繼承了 FPLS的優(yōu)點，交界面具有良好的光滑性，總體質量有良好的守恒性；另一方面由于它只在交界面附近窄帶內求解level set函數，包括進行重新初始化，同時基于四叉樹或八叉樹的自適應網格只在交界面附近加密網格，而在遠離交界面的區(qū)域保持較粗網格，因而這種方法既可節(jié)約內存，又可提高計算效率.

［1］Osher S，Sethian J. Fronts propagating with curvature dependent speed：Algorithms based on Hamiliton-Jacobi formulations[J].Journal of Computational Physics，1988，79(1)：12-49.

［2］Sussman M，Puckett E G，Osher S. A level set approach for computing solutions to incompressible two-phase flow[J].Journal of Computational Physics，1994，114(2)：146-159.

［3］Peng D，Merriman B，Osher S，et al. A PDE-based fast local level set method[J].Journal of Computational Physics，1999，155(2)：410-438.

［4］Sethian J. A fast marching level set method for monotonically advancing fronts[J].Applied Mathematics，1996，93(4)：1591-1595.

［5］Sussman M，Puckett E G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J].Journal of Computational Physics，2000，162(2)：301-337.

［6］Enright D，Fedkiw R. A hybrid particle level set method for improved interface capturing[J].Journal of Computational Physics，2002，183(1)：83-116.

［7］Enright D，Losassom F. A fast and accurate semi-Lagrangian particle level set method[J].Computers and Structure，2005，83(6)：479-490.

［8］Strain J. Semi-Lagrangian methods for level set equations[J].Journal of Computational Physics，1999，151(2)：498-533.

［9］Sethian J. Fast marching methods[J].SIAM Review，1999，41(2)：199-235.

［10］Strain J. Tree methods for moving interfaces [J].Journal of Computational Physics，1999，151(2)：616-648.

［11］Strain J. Fast tree-based redistancing for level set computations[J].Journal of Computational Physics，1999，152(2)：664-686.

［12］Sochnikov V，Efrim S. Level set calculations of the evolution of boundaries on a dynamically adaptive grid[J].International Journal for Numerical Methods in Engineering，2002，56(13)：1913-1929.

［13］Mark de B，Mark van K，Mark O，et al.Computational Geometry Algorithms and Applications[M]. Berlin：Springer，1999.

［14］Sussman M，Fatemi E. An efficient，interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow[J].SIAM Journal on Scientific Computing，1999，20(6)：1165-1191.

［15］Smolarkiewicz P. The multi-dimensional crowley advection scheme[J].Monthly Weather Review，1982，110(12)：1968-1983.

［16］Zalesak S. Fully multidimensional flux-corrected transport algorithms for fluids[J].Journal of Computational Physics，1979，31(3)：335-362.

Spatially Adaptive Fast Particle Level Set Method

HUANG Xiao-yun，LI Shao-wu
(School of Civil Engineering，Tianjin University，Tianjin 30072，China)

As an effective method of interface tracking，the level set method has been widely applied to surface tracking of single phase flow and interface tracking of multi-phase flow. The fast particle level set(FPLS)method，as an improved version of the level set method，not only is of relatively high accuracy in surface tracking，but also possesses high efficiency of calculation. By combination of the FPLS with tree-based adaptive gridding system，a spatially adaptive fast particle level set(SA-FPLS)approach is proposed，By solving the level-set equation with higher gridding resolution in a narrow mesh band around the interface，the solution accuracy can be further improved with less computer memory consumption and higher computational efficiency.

fast particle level set；adaptive grid；spatially adaptive fast particle level set

TV131.2

A

0493-2137(2010)11-0981-07

2009-09-02；

2010-05-14.

黃筱云（1980— ），男，博士研究生，td98298@sina.com.

李紹武，lishaowu@tju.edu.cn.