基于分區(qū)徑向基函數(shù)的Level Set界面捕捉研究*

金　秋　李廷秋　林　超

(武漢理工大學(xué)交通學(xué)院　武漢　430063)

?

基于分區(qū)徑向基函數(shù)的Level Set界面捕捉研究*

金秋李廷秋林超

(武漢理工大學(xué)交通學(xué)院武漢430063)

針對傳統(tǒng)的Level Set方法的固有缺陷，通過引入徑向基函數(shù)(radial basis function，RBF)，提出一種分區(qū)捕捉動界面的Level Set數(shù)值方法，將Level Set偏微分控制方程轉(zhuǎn)化為全微分方程，避免了重新初始化過程，可精確表達(dá)界面形狀并維持質(zhì)量守恒.采用徑向基函數(shù)中的基函數(shù)和權(quán)重系數(shù)來重構(gòu)符號距離函數(shù)，并結(jié)合區(qū)域劃分技術(shù)提高計算效率.3個計算算例表明，基于分區(qū)徑向基函數(shù)Level-Set方法，可明顯提高界面幾何形狀的精度和捕捉效率，具有較好的質(zhì)量守恒性.

水平集；徑向基函數(shù)；運動界面；分區(qū)技術(shù)；質(zhì)量守恒

0　引　　言

動界面捕捉技術(shù)研究是模擬海工結(jié)構(gòu)物在復(fù)雜海洋環(huán)境中運動的必要前提.雖然當(dāng)前傳統(tǒng)Level Set方法已逐步運用于自由液面及動物面捕捉的數(shù)值模擬中[1-2]，但因其空間項離散數(shù)值誤差導(dǎo)致的質(zhì)量不守恒問題，極大的限制了其應(yīng)用范圍.為了解決傳統(tǒng)Level Set方法質(zhì)量不守恒問題，Cai等[3]構(gòu)造了一種五階緊湊CWENO格式，從而提高空間項離散精度，Lsakov等[4]提出基于求解偏微分方程的快速局部Level set算法，通過減少符號距離重新初始化的迭代次數(shù)，提高Level Set算法精度.Jemison等[5]提出改進(jìn)CLSVOF方法，運用Level Set求解相界面法向方向等物理量，結(jié)合VOF方法求解流輸運問題以保證質(zhì)量守恒.Enright等[6]提出了粒子Level Set方法(PLSM)，在計算域內(nèi)設(shè)定拉格朗日標(biāo)記粒子，根據(jù)標(biāo)志粒子修正符號距離函數(shù).

文中提出一種基于分區(qū)徑向基函數(shù)(radial basis function，RBF)的Level Set方法，根據(jù)界面實際位置對全局計算域進(jìn)行分區(qū)，突破有限差分求解框架，實現(xiàn)符號距離函數(shù)的空間和時間項分離表達(dá)，將偏微分方程轉(zhuǎn)化為全微分方程進(jìn)行求解，有效限制傳統(tǒng)方法中空間項離散帶來的數(shù)值誤差，提高計算精度的同時有效維持動界面的質(zhì)量守恒.

1　理論基礎(chǔ)

1.1符號距離函數(shù)和Level Set方程

在Level Set方法中，將移動的界面作為零等值面嵌入高一維符號距離函數(shù)中，隱式地表達(dá)或追蹤.構(gòu)造函數(shù)Φ(x,t)，使得在任意時刻，運動界面Γ為Φ(x,t)的零等值面，即

(1)

定義Φ(x,t)初值為t=0時刻，x到界面的符號距離，即

(2)

在已知初始時刻界面位置后，通過對時間求全導(dǎo)數(shù)，確定描述界面運動的偏微分形式的控制方程

(3)

式中：v為目標(biāo)點處的速度，由界面處的法向速度Fn和單位外法向向量Φ表示，則上式可寫為

(4)

將界面Γ的演化嵌入至函數(shù)Φ中，通過對Φ的求解，可自動獲得界面Γ的拓?fù)浣Y(jié)構(gòu)變化，如界面分裂或界面融合.

傳統(tǒng)LevelSet方法中，運用符號距離函數(shù)賦初值并用迎風(fēng)格式的FDM求解運動控制方程，該局部近似的方法產(chǎn)生的數(shù)值誤差將逐步累積最終影響求解結(jié)果，因此需要周期性的重新初始化過程來保證數(shù)值求解的穩(wěn)定性.

1.2徑向基函數(shù)

給定函數(shù)Φ：R+→R，對于數(shù)據(jù){Xi,fi}∈Rd?R,(i=1,2,…,n)，表示為如下徑向基函數(shù)的形式

(5)

式中：f(X)為待求物理量(如空間點到界面的距離)；X為計算域中任意空間點；Xi為插值點；λ為權(quán)重系數(shù)；n為插值點個數(shù)；P(X)為一階多項式.

對于二維問題，其表達(dá)式為P(X)=λn+1+λn+2x+λn+3y，φ(‖·‖)為徑向基基函數(shù)，文中選用MQ基函數(shù)

(6)

式中:c為形狀參數(shù)，取0～1間的常數(shù).

為保證RBF插值有惟一解，權(quán)重系數(shù)須滿足正交性質(zhì)

(7)

n個插值點組成的徑向基函數(shù)組寫成矩陣形式為

(8)

Aij=φj(xi),i,j=1,2,…,N

Pij=pj(xi),i=1,2,…,N;j=1,2,3.

1.3基于RBF求解Level Set方程

基于徑向基函數(shù)的散亂點插值功能，對LevelSet函數(shù)符號距離函數(shù)進(jìn)行重新描述.在表達(dá)形式上實現(xiàn)時間變量和空間變量的分離，當(dāng)插值點固定不變時，空間項即為已知項，函數(shù)中僅有RBF權(quán)重系數(shù)與時間相關(guān).

(9)

式中：ψ(x)=[φ1(x)…φN(x)1xy]

因此，LevelSet運動方程式(4)轉(zhuǎn)變?yōu)?/p>

(10)

LevelSet運動偏微分方程由此轉(zhuǎn)變?yōu)殛P(guān)于時間的全微分方程，通過更新RBF權(quán)重系數(shù)，實現(xiàn)更新符號距離函數(shù)的目標(biāo)，繼而實現(xiàn)動界面的隱式表達(dá).

結(jié)合式(8)、(10)可以改寫為

(11)

式中：

采用一階差分方法對時間項進(jìn)行離散，即可得到下一時刻權(quán)重系數(shù)值.

(12)

式中：Δt為時間步長.

1.4區(qū)域劃分插值技術(shù)

由1.2可知，當(dāng)插值點個數(shù)為n時，每一個時間步迭代需要求解一個(n+3)×(n+3)的滿秩矩陣.根據(jù)實際計算工況對全局計算域劃分區(qū)域，可減少遍歷全局背景網(wǎng)格的區(qū)域，減少插值點數(shù)量，顯著提高矩陣求解速度.同時，在界面區(qū)域劃分更細(xì)的局域網(wǎng)格，有利于更精確的捕捉界面拓?fù)浣Y(jié)構(gòu)變化.局域區(qū)域網(wǎng)格劃分示意圖見圖1.

圖1　區(qū)域劃分示意圖

定義一個接近零的足夠小數(shù)正值δ(一般取為1.5～2倍均勻網(wǎng)格寬度)，根據(jù)函數(shù)值Φ(x,y)進(jìn)行判斷，在計算域內(nèi)判斷Φ(x,y)與δ的關(guān)系，當(dāng)Φ(x,y)≤δ時記錄所有的點集，分別求得xmin{(x,y)|Φ(x,y)≤δ},xmax{(x,y)|Φ(x,y)≤δ},ymin{(x,y)|Φ(x,y)≤δ},ymax{(x,y)|Φ(x,y)≤δ},定位界面位置，保證所劃分的局部插值域包含整個界面.引入dx,dy兩個安全范圍，滿足dx>max(ux×Δt),dy>max(vy×Δt)，以此確定局部插值域的四個角點坐標(biāo)為：(xmin-dx，ymin-dy)，(xmax+dx，ymin-dy)，(xmax+dx，ymax+dy)和(xmin-dx，ymax+dy).由于在計算過程中時間步長根據(jù)CFL條件選擇，在均勻網(wǎng)格內(nèi)一般選取dx,dy為2-3倍的網(wǎng)格寬度.

2　計算結(jié)果及分析

2.1分區(qū)RBF-Level Set與傳統(tǒng)Level Set數(shù)值實驗對比

用上述2種方法分別對典型圓形界面在剪切流中的變形現(xiàn)象進(jìn)行模擬計算，在相同的計算條件下，將文中提出的分區(qū)RBF-Level Set方法與傳統(tǒng)Level Set方法進(jìn)行對比，以驗證新方法對質(zhì)量守恒性的改善效果.

算例1剪切流場.

速度場為

(13)

由圖2a)可知，設(shè)定計算域為(0,0)×(1,1)，初始時刻圓心位置(0.35,0.5)，半徑R=0.25 m，剪切流動旋轉(zhuǎn)中心(0.5,0.5).

圖2　分區(qū)RBF-Level Set與傳統(tǒng)Level Set剪切流場結(jié)果對比

本案例中計算了剪切流中圓形界面在不同時刻的運動變形情況.對比可知，分區(qū)RBF-Level Set方法能更精確捕捉運動界面的形狀與位置.統(tǒng)計不同時間點RBF-Level Set方法和傳統(tǒng)Level Set方法界面內(nèi)質(zhì)量的相對誤差，分別為0.15%/1.05%，0.35%/1.9%和0.75%/2.5%，明顯提高了計算結(jié)果的精確度.對比驗證結(jié)果表明，基于徑向基函數(shù)的Level Set方法可以更好地維持界面形狀同時保證質(zhì)量守恒特性.

2.2分區(qū)RBF-Level Set與全局RBF-Level Set方法對比

算例2剪切流場，流場信息同算例2.

圖3　分區(qū)與全局RBF-Level Set剪切流場結(jié)果對比

由圖3可知，2種計算方法均能較好的模擬剪切流中界面變形，計算結(jié)果基本完全吻合，但在界面變形的細(xì)節(jié)處，如界面尾端尖點處，分區(qū)RBF-LevelSet方法的捕捉更為精確細(xì)致.全局方法中每一步的計算遍歷整個計算域，計入全局速度場對界面的影響，但消耗了大量計算時間，且單個時間步內(nèi)遠(yuǎn)離界面處的速度對界面變形影響較小.文中提出的分區(qū)思想隨界面移動更新局部區(qū)域，劃分局部細(xì)網(wǎng)格，將插值點集中于界面周圍，節(jié)省時間的同時更好的模擬界面細(xì)節(jié)處的變化.

3　分區(qū)RBF-Level Set方法界面強(qiáng)迫運動模擬

算例3船球鼻首剖面的橫搖運動.

速度場

(14)

設(shè)定計算域為計算域為(0,0)×(2,2)，初始時刻球鼻艏處于平衡位置，圖3分別表示球鼻艏運動1/4,1/2,3/4,1個周期位置.

圖4　RBF-Level球鼻艏橫搖運動

利用文中所提出的分區(qū)RBF-Level Set方法對船體曲面曲率變化較為復(fù)雜的球艏剖面的強(qiáng)迫橫搖運動進(jìn)行捕捉，由圖中可以看出，RBF-Level Set界面捕捉結(jié)果與模型初始時刻形狀十分吻合，尖銳處仍能保持較高分辨率.統(tǒng)計橫搖各時刻界面內(nèi)質(zhì)量誤差，分別為0.08%，0.06%，0.1%，計算結(jié)果表明，文中所提出的RBF-Level Set方法也適用于曲率變化較大的界面運動捕捉.

4　結(jié) 束 語

提出了一種基于分區(qū)RBF-Level Set方程求解方法.數(shù)值實驗結(jié)果表明，運用RBF插值技術(shù)分區(qū)求解Level Set方程繼承了傳統(tǒng)方法精確表達(dá)運動界面幾何形狀優(yōu)點的同時，由于避免了重新初始化過程，明顯改善了質(zhì)量守恒性，具有更好的數(shù)值穩(wěn)定性和精度.

[1]CHENY Y, BOTELLA O. The LS-STAG method: a new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[J]. Computational Physics,2010,229(4):1043-1076.

[2]KRAUSE R, MOHR C. Level set based multi-scale methods for large deformation contact problems[J]. Applied Numerical Mathematics，2011，61:428-442.

[3]CAI L, FENG J H, XIE W X. A CWENO-type central-upwind scheme for ideal MHD equations [J] .Applied Mathmatics and Computation,2005(1):600-612.

[4]ISAKOV V, LEUNG S, QIAN J. A fast local level set method for inverse gravimetry[J]. Communications in Computational Physics,2011,10(4):1044-1070.

[5]JEMISON M, LOCH E, SUSSMAN M,et al. A coupled level set-moment of fluid method for incompressible two-phase flows [J]. Journal of Scientific Computing, 2013,54(3):454-491.

[6]ENRIGH T D, FEDKIW R, FERZIGER J, et al. A hybrid particle level set method for improved interface capturing [J]. Journal of Computational Physics, 2014,183(2):83-116.

Regional RBF based Level Set Method for Interface Capturing

JIN QiuLI TingqiuLIN Chao

(SchoolofTransportation,WuhanUniversityofTechnology,Wuhan430063,China)

The renewal of signed distance function usually relies on reinitialization in conventional Level Set method. An alternative method is proposed in this paper to regionally capture the moving interfaces via the introduction of the Radial Basis Functions (RBF). This approach provides an accuracy and conserved representation of the implicit function, location and shape of interfaces. The signed distance function is rebuilt by based functions and weight coefficients. Its deformation is then considered as an updating of the RBF interpolants, turning the partial differential equation (PDE) problem into an ordinary differential equation (ODE) problem. Reinitialization is found to be no longer necessary and hence it avoids large numerical errors occurring in the conventional method. Three numerical examples show that the proposed regional RBF-Based Level Set method can significantly improve both the accuracy and conservation in capturing moving interfaces.

level set; RBF; moving interfaces; regional technology; mass conservation

2016-07-05

U661.1

10.3963/j.issn.2095-3844.2016.04.013

金秋(1992- )：女，碩士生，主要研究領(lǐng)域為船舶流體力學(xué)

*國家自然科學(xué)基金項目資助(51579196、51139005)