空間自適應(yīng)的快速粒子Level Set方法

黃筱云，李紹武

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

空間自適應(yīng)的快速粒子Level Set方法

黃筱云，李紹武

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

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

快速粒子level set；自適應(yīng)網(wǎng)格；空間自適應(yīng)的快速粒子level set

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

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

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

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

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

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

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

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

(1)按照初始化 level set函數(shù)，生成在交界面附近加密的樹型結(jié)構(gòu)網(wǎng)格；

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

(3) 根據(jù)level set函數(shù)值重新調(diào)整網(wǎng)格；

(4) 返回第(2)步.

1 FPLS方法的基本原理

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

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

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

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

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

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

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

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

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

2 SA-FPLS的原理

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

圖2 樹型網(wǎng)格結(jié)構(gòu)Fig.2 Tree grid structure

圖3 網(wǎng)格及其鄰居Fig.3 Given grid and its neighboring grid

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

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

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

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

3 算例及分析

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

3.1 算例1

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

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

圖5 七角星結(jié)構(gòu)擴張過程的數(shù)值模擬結(jié)果Fig.5 Numerical results of expansion of seven-pointed star

3.2 算例2

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

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

3.3 算例3

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

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

圖6 圓在簡單旋轉(zhuǎn)流場中變形過程的數(shù)值模擬結(jié)果Fig.6 Numerical results of transformation of circle in single vortex flow

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

檢驗?zāi)Ｐ褪睾阈詢?yōu)劣的方法是對比初始時刻界面所圍圓面積與 1個周期后界面恢復(fù)到原始形狀后面積的差異.界面所圍面積的計算式為

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

3.4 算例4

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

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

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

圖8 Zalesak圓球的旋轉(zhuǎn)過程數(shù)值模擬結(jié)果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]，網(wǎng)格層數(shù)為 3，最粗網(wǎng)格大小為0.04× 0 .04× 0 .04.流場表示為式中，流場周期 T為 2,s，計算時間步長為 0.005,s.結(jié)果表明，圓球在1.0,s時達到最大拉伸后，又恢復(fù)到原樣，而其中各圖為圓球分別在 t為0～2.0,s時的形狀.

4 結(jié) 語

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

［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，F(xiàn)edkiw 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，F(xiàn)atemi 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.