基于分區徑向基函數的Level Set界面捕捉研究*

金　秋　李廷秋　林　超

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

?

基于分區徑向基函數的Level Set界面捕捉研究*

金秋李廷秋林超

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

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

水平集；徑向基函數；運動界面；分區技術；質量守恒

0　引　　言

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

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

1　理論基礎

1.1符號距離函數和Level Set方程

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

(1)

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

(2)

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

(3)

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

(4)

將界面Γ的演化嵌入至函數Φ中，通過對Φ的求解，可自動獲得界面Γ的拓撲結構變化，如界面分裂或界面融合.

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

1.2徑向基函數

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

(5)

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

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

(6)

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

為保證RBF插值有惟一解，權重系數須滿足正交性質

(7)

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

(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方程

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

(9)

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

因此，LevelSet運動方程式(4)轉變為

(10)

LevelSet運動偏微分方程由此轉變為關于時間的全微分方程，通過更新RBF權重系數，實現更新符號距離函數的目標，繼而實現動界面的隱式表達.

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

(11)

式中：

采用一階差分方法對時間項進行離散，即可得到下一時刻權重系數值.

(12)

式中：Δt為時間步長.

1.4區域劃分插值技術

由1.2可知，當插值點個數為n時，每一個時間步迭代需要求解一個(n+3)×(n+3)的滿秩矩陣.根據實際計算工況對全局計算域劃分區域，可減少遍歷全局背景網格的區域，減少插值點數量，顯著提高矩陣求解速度.同時，在界面區域劃分更細的局域網格，有利于更精確的捕捉界面拓撲結構變化.局域區域網格劃分示意圖見圖1.

圖1　區域劃分示意圖

定義一個接近零的足夠小數正值δ(一般取為1.5～2倍均勻網格寬度)，根據函數值Φ(x,y)進行判斷，在計算域內判斷Φ(x,y)與δ的關系，當Φ(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)，以此確定局部插值域的四個角點坐標為：(xmin-dx，ymin-dy)，(xmax+dx，ymin-dy)，(xmax+dx，ymax+dy)和(xmin-dx，ymax+dy).由于在計算過程中時間步長根據CFL條件選擇，在均勻網格內一般選取dx,dy為2-3倍的網格寬度.

2　計算結果及分析

2.1分區RBF-Level Set與傳統Level Set數值實驗對比

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

算例1剪切流場.

速度場為

(13)

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

圖2　分區RBF-Level Set與傳統Level Set剪切流場結果對比

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

2.2分區RBF-Level Set與全局RBF-Level Set方法對比

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

圖3　分區與全局RBF-Level Set剪切流場結果對比

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

3　分區RBF-Level Set方法界面強迫運動模擬

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

速度場

(14)

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

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

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

4　結 束 語

提出了一種基于分區RBF-Level Set方程求解方法.數值實驗結果表明，運用RBF插值技術分區求解Level Set方程繼承了傳統方法精確表達運動界面幾何形狀優點的同時，由于避免了重新初始化過程，明顯改善了質量守恒性，具有更好的數值穩定性和精度.

[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- )：女，碩士生，主要研究領域為船舶流體力學

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