一種徑向基函數虛擬網格法數值模擬復雜邊界流動

辛建建 石伏龍 金秋

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

一種徑向基函數虛擬網格法數值模擬復雜邊界流動

辛建建?石伏龍 金秋

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

(2016年8月7日收到;2016年11月25日收到修改稿)

提出了一種徑向基函數虛擬網格法浸入邊界法以模擬復雜或多體浸入邊界黏性繞流問題.在該方法中,采用有限差分法離散固定笛卡爾交錯網格上的不可壓縮Navier-Stokes方程,以分步法結合三階Runge-Kutta格式進行時間積分,高階TVDMUSCL(total variation diminishing monotonic upstream-centered scheme for conservation laws)格式離散對流項;一個邊界連續的虛擬網格法施加物面邊界條件以考慮尖銳邊界對流場的影響;引入具多項式基的徑向基函數描述和重構任意復雜浸入界面,并識別背景網格屬性狀態.采用Fortran 90語言開發相應的求解器,模擬了繞圓柱、機翼和交錯布置圓柱群的黏性繞流問題,驗證了本文方法的正確性、可靠性和對復雜邊界流動問題的適用性.

浸入邊界法,徑向基函數,虛擬網格法,復雜邊界

1 引 言

數值模擬繞復雜形狀物體的黏性繞流問題一直是計算流體力學(computational fluid dynamics,CFD)領域的研究重點和難點.傳統研究方法是應用貼體動網格技術,界面附近網格隨物體一起運動以精確地描述運動界面.對于小幅運動或單物體時,該方法可得到較好的結果;對于復雜形狀物體,生成單域貼體的動網格非常復雜,當物體發生大幅運動或拓撲結構大變形時,網格再生過程會顯著降低計算效率,重構過程中也會引入插值誤差.

近年來,針對復雜形狀多體動邊界的數值模擬問題,浸入邊界法[1](immersed boundary method,IBM)受到極大的關注.其基于固定的背景笛卡爾網格,網格生成簡單,無需進行網格重構,適于模擬涉及復雜幾何形狀、多體或大變形的動邊界問題.根據流-固界面的表示不同,浸入邊界法可分為擴散界面法和尖銳界面法[2].對于擴散界面法,由于引入光滑Delta函數或假設有限厚度界面,導致界面模糊化,計算精度受到一定的限制,難以模擬高雷諾數流動.近年來發展的一類尖銳界面法旨在彌補這個基本缺欠,維持界面尖銳是現代CFD技術的發展趨勢,如虛擬網格法[3]、切割網格法[4]、嵌入邊界法[5].

本文重點強調虛擬網格法,其最初由Majumdar等[6]提出,基本思想是構建物體內的虛擬網格以施加邊界條件,考慮尖銳浸入邊界對流場的影響.虛擬網格定義為至少有一個相鄰流體網格的固體網格,對于每個虛擬網格,需要一個內插模板以隱式施加物面邊界條件.本文采用雙線性(三線性對于三維)插值[7]方法重構界面內插格式,該方法實施簡單、穩定性好.然而對于高雷諾數流動,線性重構可能導致計算錯誤,一個更高精度的重構方法是在切向線性和在法向進行二次插值,或者采用基于泰勒展開的高階插值方法[8].

相比需要引入人為剛度系數的反饋力法[9],虛擬網格法避免了時間步相關的穩定性問題;相比在動量方程中引入力源項的嵌入邊界法[10],虛擬網格法維持界面尖銳,適用于高雷諾數流動問題;相比網格分類離散的切割網格法[11],虛擬網格法以統一格式離散背景網格,易擴展到三維和曲線網格.發展至今,虛擬網格法已經得到廣泛的應用研究,如繞圓柱和機翼的可壓縮流動[12]、繞圓柱和平板的不可壓縮流動[13]、熱傳輸[14]、昆蟲飛行和游魚運動[15].對于國內研究方面,胡偶等[16]模擬了繞圓柱和機翼的二維無黏可壓縮流動.

本文提出一種簡單靈活的徑向基函數虛擬網格浸入邊界法,應用于模擬任意復雜浸入界面的不可壓縮流動問題.針對復雜或多體浸入界面流動的模擬,本方法有兩個突出要點：1)采用一個邊界連續的虛擬網格法施加物面邊界條件以隱式考慮浸入邊界對流場的影響,物體作為虛擬存在的邊界,在整個計算域內采用統一的離散方法求解Navier-Stokes方程;2)以具多項式基的徑向基函數(polynomial and radial basis function,PRBF)[17]描述和重構任意復雜浸入界面,并定位和追蹤背景網格相狀態.基于Fortran 90語言自主開發的二維浸入邊界法求解器,模擬了圓柱繞流、機翼繞流和圓柱群的繞流問題.計算結果與已有的數值和試驗數據符合較好,驗證了本文方法的精度和可靠性,證明了本文方法模擬復雜或多體界面繞流問題的能力.

2 數學模型和數值方法

2.1 控制方程

非定常黏性不可壓縮流體的控制方程如下[18]：

速度向量u=(u,v),u,v分別是直角網格x-,y-方向速度分量;ρ是密度,t是時間,p是壓力,μ是流體的動力黏性系數.

2.2 數值離散方法

以上流體控制方程采用有限差分法在交錯直角網格上離散,速度定義在網格面心,壓力、密度、黏性、距離函數等標量定義在網格體心.在該求解器中,采用分步法解耦速度和壓力[19].時間推進格式采用三階Runge-Kutta(RK3)格式,對流項顯式處理,黏性項以Crank-Nicolson格式半隱式處理.對于空間離散格式,對流項采用一種高階TVD MUSCL(total variation diminishing monotonic upstream-centered scheme for conservation laws)[20]格式離散,黏性項以中心差分格式離散.

圖1給出了本文流固耦合數值求解方法的流程圖,對于該方法的兩個關鍵數值技術,物面邊界條件施加和復雜浸入界面描述分別在下文2.3節和2.4節詳細介紹.

2.3 虛擬網格浸入邊界法

在虛擬網格法中,在每個時間步,虛擬網格變量值通過浸入邊界和臨近的流體網格變量值插值計算,本文采用雙線性插值方法.物面邊界條件通過鏡像方法直接滿足,以確保數值穩定性.

虛擬網格法的關鍵一步是描述和重構物面邊界,以定位識別背景網格相狀態(流體網格、虛擬網格、固體網格),具體過程見2.4節.流體網格為網格中心位于流體區域的網格,固體網格為除去虛擬網格后網格中心位于固體區域內的網格.第二步是插值重構虛擬網格變量值,該過程由如下子步驟組成.

圖1 浸入邊界法數值算法流程圖Fig.1.Flowchart of numerical algorithmfor immersed boundary method.

1)定義鏡像點：對于每個虛擬網格,穿過浸入邊界反射虛擬網格到流體區域惟一確定一個鏡像點.本文鏡像點定義方法類似于Pan和Shen[22]的方法,其實施簡單,可惟一確保鏡像點插值模板.定義鏡像點到物體表面的間距為ΔL,ΔL是人為指定距離.對于二維問題,dx,dy分別為該流體網格x,y方向的網格尺寸.該ΔL的選擇是確保鏡像點在被四個流體網格節點形成的一個矩形區域內,而不被邊界表面穿過.鏡像點坐標計算為

(xIP,yIP)為鏡像點網格坐標,(xGC,yGC)為虛擬網格點坐標,(xBI,yBI)為沿著物面法向從虛擬網格中心延伸與物面的投影點坐標,dGB=為虛擬網格點與投影點之間的距離,相關細節如圖2.

2)插值計算鏡像點流體變量：在上一步中圍繞鏡像點的四個流體網格節點即為插值模板點,分別用1,2,3,4表示.鏡像點值是1,2,3,4節點值的權重疊加,鏡像點變量值φIP以雙線性插值方法計算：

[φ]T為1,2,3,4節點變量值;[B]為雙線性內插格式的系數矩陣,表達式為

(x1,y1),(x2,y2),(x3,y3),(x4,y4)分別為節點1,2,3,4的直角網格坐標值.系數向量[A]從(8)式計算得到,通過(7)式就可計算出鏡像點的變量值.

圖2 虛擬網格法的二維插值格式Fig.2.Interpolation stencil of ghost cell method.

圖3 鄰近的網格節點表示Fig.3.Stencil for adjacent grid nodes.

3)計算虛擬網格值：鏡像點值φIP計算后,沿著物面法向采用線性內插方法就可計算出虛擬網格節點值φGC：

對于速度,需滿足無滑移邊界條件,此時C1=-d1/d2和C2=φBI(d1+d2)/d2,φBI為物面節點速度值,d1為虛擬網格到投影點距離,d2為投影點到鏡像點距離.對于壓力,需要滿足Neumann邊界條件,此時C1=1和C2=N·(d1+d2),N為物面法向壓力導數.

以虛擬網格法施加邊界條件后,在全場(包括物面內)求解離散Navier-Stokes方程即可考慮浸入邊界對流場的影響.區別于傳統虛擬網格法只在物面外計算流場,本文虛擬網格法在整個計算區域求解.與固體區域重合的流體域稱為虛擬流體區域,在該區域固體速度和流體網格速度滿足無滑移邊界條件.由于固體與流場的相互作用是通過邊界條件實現,因此固體內虛擬流體的存在對外部流場沒有影響.該方法雖然增加了些許計算量,但極大地簡化了數據結構和浸入邊界的處理,降低了編程難度.

2.4 PRBF界面描述方法

浸入邊界法的網格生成過程得到極大簡化,但由于浸入邊界本身并未包含在網格中,精確表示在背景網格中的浸入邊界極具挑戰性.本文采用基于PRBF的隱式等值面方法重構任意界面,該方法以光滑基函數擬合離散的物面控制點構造任意復雜物體的等值曲面,并且根據等值面距離函數正負號識別場點屬性.PRBF[23,24]界面描述方法簡單、精度高,適于模擬復雜浸入界面.

針對任意浸入界面,基于有限數量的界面位置信息的離散控制點,引入PRBF構造如下場函數：

其中M為樣品點總數,(xi,yi)為某點i的坐標值,rij表示兩點間距離.φ(‖rij‖)在本文采用MQ核函數表示：

形狀參數C為常數,這里取C=0.徑向基函數權重系數αj滿足正交條件：

多項式PF(xi)為

浸入界面的場函數d(xi)均具有等值距離λ(λ＞0),λ為大于零的任意常數,本文取λ=1.因此,等值面方程可轉化為求解關于系數αj的(M+3)階線性方程組.

MQ函數是條件正定的,增加多項式基后,方程(15)的系數矩陣即變為正定矩陣,確保了方程(15)存在惟一確定的解.另外,由于MQ核函數是在所有物面控制點進行(全域求解),其系數矩陣是滿陣的.PRBF隱式界面描述方法以有限的物面控制點即可獲得較高的精度,采用傳統的高斯消元法即可快速地求解等值面方程(15).

進一步構造如下等值面函數ψ(x),以有效識別網格節點ψ(x)在場中屬性：

等值面函數ψ(x)具有如下屬性：對于任意坐標節點xP(見圖3節點xP與鄰近網格的符號表示),若ψ(xP)＞0,表示點(xP)在物體內;ψ(xP)＜0表示點(xP)在物體外.在求出全場的等值函數ψ(x)后,下一步是根據圖3中網格節點相對位置的距離函數識別定位虛擬網格、流體網格、固體網格 (如圖4).

1)如果ψ(xP)≥ 0和 (ψ(xE)≥ 0和ψ(xW)≥0和ψ(xN)≥0和ψ(xS)≥0),則為固體網格(S);

2)如果ψ(xP)≥ 0和 (ψ(xE)＜0或ψ(xW)＜0或ψ(xN)＜0或ψ(xS)＜0),則為虛擬網格(G);

3)如果ψ(xP)＜0和 (ψ(xE)＜0和ψ(xW)＜0和ψ(xN)＜0和ψ(xS)＜0),則為流體網格(F).

圖4 背景網格相狀態識別(F,流體網格;G,虛擬網格;S,固體網格)Fig.4.Identifying phase state of background grid(F,fl uid cell;G,ghost cell;S,solid cell).

3 算例與結果分析

浸入邊界繞流是經典的流動問題,具有豐富的物理現象,廣泛地存在于實際的工程應用中.本文將對雷諾數Re=40,100,200的單圓柱繞流,Re=500的NACA0012機翼繞流和Re=100,200的圓柱群繞流問題進行數值模擬,以驗證其數值方法的正確性.

3.1 單圓柱的二維黏性流動

在該圓柱繞流算例中,計算區域設置為[-8D,16D]×[-8D,8D],D為圓柱直徑,大的計算域使邊界遠離圓柱域以避免外部邊界尾流影響中心流場,幾何模型如圖5.在入口處采用Dirichlet速度邊界條件,頂部和底部均采用對稱邊界,出口處采用出流邊界條件.對于壓力邊界條件,在出口指定壓力為零,其余邊界采用均勻Neumann邊界條件.直徑D設為1 m,入流速度U∞為1 m/s.Re=40,100,200時對應的運動黏性系數ν分別為0.025,0.01,0.005 m2/s(Re=U∞D/ν). 采用兩套網格系統480×360和240×180(分別稱為M1,M2)進行網格細化研究,并取固定時間步長Δt=0.001 s.

表1比較了三種雷諾數Re=40,100,200下圓柱繞流算例的時間平均阻力系數CD、平均升力系數CL和回流長度Lc,并與其他浸入邊界法結果進行了比較.從中可看出,本文虛擬網格法計算的三個關鍵參數結果(升阻力系數和回流長度)與其他數值結果符合較好,驗證了本文方法的精度.另外,兩種網格系統的收斂性研究表明低雷諾數下粗細兩種網格系統的數值結果趨近相同,說明了本文方法的可靠性.

圖5 圓柱繞流的幾何模型Fig.5.Geometry model of flow around cylinder.

圖6 不同雷諾數下圓柱繞流的瞬時渦場圖 (a)Re=40;(b)Re=100;(c)Re=200Fig.6.Instantaneous vorticity contours of flow around cylinder in three Reynolds number：(a)Re=40;(b)Re=100;(c)Re=200.

圖6展示了采用網格系統M1計算在Re=40,100和200時圓柱尾流附近的瞬時渦場.在Re=40時,我們能看到一對穩定的對稱渦附著在圓柱后側(如圖6(a)),此時流動為定常層流,升力系數為零.隨著雷諾數增加,尾流變得不穩定,渦開始周期性交替地從圓柱后面脫落,形成Karman渦街,而且渦泄現象愈發明顯,該現象可從Re=100,Re=200時的瞬時渦場圖看出(圖6(b)和圖6(c)).由于Karman渦街的產生,導致圓柱體表面的升阻力發生周期性的變化.該現象可從圖7中Re=200時M1網格系統計算的升阻力系數時間演化過程觀察到.

表1 圓柱繞流算例中的升阻力系數比較Table 1.Comparison of lift and drag coefficient in the case of flow around cylinder.

圖7 Re=200時升阻力系數的時間歷程曲線Fig.7.Evolution history of drag and lift coefficients atRe=200.

3.2 翼型繞流

本算例研究雷諾數Re=500的NACA0012機翼繞流問題,以驗證本文求解器模擬復雜浸入界面的能力.計算區域取為[-3b,9b]×[-3b,3b],采用600×400的非均勻網格系統,固定的時間步長取為dt=0.001 s,b為翼弦長度.機翼攻角取為0?,并取機翼前緣點為坐標原點.入口邊界為均勻流,在出口采用對流邊界條件,自由流邊界采用輻射邊界.根據Imamura等[26]和李秋[27]的工作選擇流體參數,雷諾數Re=U∞b/ν=500.NACA0012機翼的幾何構型復雜,難以用解析函數表示,這里采用多項式徑向基函數擬合141個樣品點表示機翼構型.

當前方法計算的阻力系數 (CD=0.1653)比Imamura等[26]在最細網格計算的值(CD=0.1725)稍小,誤差在接受范圍內.表明本文CFD求解器和PRBF方法能較好地模擬復雜浸入邊界流動問題.從圖8的壓力等值線圖可看出,壓力沿著翼弦方向對稱布置,在機翼前緣點可觀察到一個壓力駐點.圖9展示了y-軸方向的速度剖面u,v,當前結果與Imamura等[26]的結果總體符合較好.局部差別可能由于兩者邊界處理方式不同,本文采用在物面內布置虛擬網格的方式施加邊界條件.

圖8 機翼的瞬時壓力場圖Fig.8.Instantaneous pressure contour of airfoil.

3.3 交錯布置的13個圓柱繞流

這里模擬了13個規則布置的圓柱黏性繞流問題,以驗證當前方法對多體浸入邊界的適用性.在計算區域[-10D,20D]×[-10D,10D]上布置均勻網格,在垂向和橫向相鄰的圓柱間距為D,交錯圓柱間距為0.5D,時間步長取為0.001 s.圖10和圖11分別展示了Re=100和Re=200時圓柱群附近位置的壓力和渦等值線流場圖.從中可看出,在Re=100時,渦場和壓力場比較穩定,但已不再是固定不動的對稱渦.隨著雷諾數增加,在Re=200時流場變成具有復雜渦泄的非定常流動.其與圖12中圓心坐標(0,0)處的圓柱在Re=100,Re=200時阻力系數時間歷程曲線反映的現象一致.在圖10中,Re=100時的阻力系數趨向于穩定,周期性現象不明顯,Re=200時阻力系數在經歷較長一段時間不規則變化后,產生周期性的升阻力系數時間變化.而在單圓柱繞流中Re=100時周期性渦泄和阻力系數已經發生,表明圓柱之間的相互作用抑制了渦泄的產生.該結果與Gao和Tseng[28]在18圓柱群樁繞流中得出的結論是一致的.另外,在Re=100和Re=200時的阻力系數值明顯小于單圓柱繞流計算的結果(對比表1),反映了多體之間的水動力干涉對圓柱表面壓力載荷的影響.

圖10 Re=100時的瞬時流場圖 (a)壓力場圖;(b)渦場圖Fig.10.Instantaneous fluid field of cylinder group atRe=100：(a)Pressure contours;(b)vorticity contours.

這里對多圓柱群樁繞流的水動力特性和多體間水動力干涉的影響進行了初步分析,對于其中的物理機理不在本文研究范圍內.分析結果表明,本文方法能較好地考慮多個浸入物體對流場的影響,可應用于模擬多個物體繞流問題.

圖11 Re=200時的瞬時流場圖 (a)壓力場圖;(b)渦場圖Fig.11.Instantaneous fluid field of cylinder group atRe=200：(a)Pressure contours;(b)vorticity contours.

圖12 圓心坐標(0,0)的圓柱在Re=100,200時的阻力系數比較Fig.12.Comparison of drag coefficients of cylinder with center coordinates(0,0)atRe=100,200.

4 結 論

本文提出一種徑向基函數虛擬網格浸入邊界法模擬復雜邊界流動.一個PRBF隱式界面描述方法以表面樣品點描述任意復雜幾何形狀,高效地識別網格屬性狀態.在整個區域包括物體內部計算Navier-Stokes方程,以虛擬網格隱式考慮浸入邊界對流場的影響.在圓柱繞流算例中驗證了本文方法的可靠性和精度,在機翼繞流和圓柱群繞流中驗證了本文方法模擬復雜幾何形狀或多體邊界繞流的能力.本文方法能輕易地擴展到復雜動邊界流動問題,這也是下一階段我們的研究重點.

[1]Li Q,Li W M 2016Acta Phys.Sin.65 064601(in Chinese)[李強,李五明 2016物理學報 65 064601]

[2]Sotiropoulos F,Yang X 2014Prog.Aerosp.Sci.65 1

[3]Lee J,You D 2013J.Comput.Phys.233 295

[4]Dechristé D,Mieussens L 2015J.Comput.Phys.314 465

[5]Yang J,Stern F 2012J.Comput.Phys.231 5029

[6]Majumdar S,Iaccarino G,Durbin P 2001Annual Research Briefs30 353

[7]Tseng Y H,Ferziger J H 2003J.Comput.Phys.192 593

[8]Xia J,Luo K,Fan J 2015Int.J.Heat Mass Transfer89 856

[9]Yang Q,Cao S Y,Liu S Y 2014Acta Phys.Sin.63 214702(in Chinese)[楊青,曹曙陽,劉十一2014物理學報63 214702]

[10]Monasse L,Daru V,Mariotti C,Piperno S,Tenaud C 2012J.Comput.Phys.231 2977

[11]Schneiders L,Günther C,Meinke M,Schr?der W 2016J.Comput.Phys.311 62

[12]Ghias R,Mittal R,Dong H 2007J.Comput.Phys.225 528

[13]Berthelsen P A,Faltinsen O M 2008J.Comput.Phys.227 4354

[14]Xia J,Luo K,Fan J 2014Int.J.Heat Mass Transfer75 302

[15]Mittal R,Dong H,Bozkurttas M,Najjar F M,Vargas A,von Loebbecke A 2008J.Comput.Phys.227 4825

[16]Hu O,Zhao N,Lin J M,Wang D H 2011Acta Aerodyn.Sin.29 491(in Chinese)[胡偶,趙寧,劉劍明,王東紅2011空氣動力學學報 29 491]

[17]Zhang X,Liu Y,Ma S 2009Adv.Mech.39 1(in Chinese)[張雄,劉巖,馬上2009力學進展 39 1]

[18]Versteeg H K,Malalasekera W 2000An Introduction to Computational Fluiddynamics:the Finite Volume Method(Second edition)(Beijing：World Book Inc)p24

[19]Lee J,Kim J,Choi H,Yang K S 2011J.Comput.Phys.230 2677

[20]ShinBR2000EuropeanCongressonComputational Methods in Applied Sciences and Engineering,Barcelona,Spain,September 2-6,2000 p1

[21]Kershaw D S 1978J.Comput.Phys.26 1

[22]Pan D,Shen T T 2009Int.J.Numer.Method Fluid60 1378

[23]Rendall T C S,Allen C B 2007Int.J.Numer.Method Eng.74 10

[24]Ye J J,Li T Q 2013Proceedings of National Congress on HydrodynamicsZhoushan,China,September 20-25,2013 p343(in Chinese)[葉俊杰,李廷秋2013全國水動力學學術會議舟山,中國,2013年9月20-25日第343頁]

[25]Frisani A,Hassan Y A 2015Comput.Fluids121 51

[26]Imamura T,Suzuki K,Nakamura T,Yoshida Masahiro 2005J.Comput.Phys.202 645

[27]Li Q 2010M.S.Thesis(Nanjing：Nanjing University of Aeronautics and Astronautics)(in Chinese)[李秋2010碩士學位論文(南京：南京航空航天大學)]

[28]Gao T,Tseng Y H,Lu X Y 2007Int.J.Numer.Method Fluid55 1189

PACS:47.85.Dh,47.11.Bc,45.50.Jf DOI:10.7498/aps.66.044704

Numerical simulation of complex immersed boundaryflow by a radial basis function ghost cell method

Xin Jian-Jian?Shi Fu-Long Jin Qiu

(School of Transportation,Wuhan University of Technology,Wuhan 430063,China)

7 August 2016;revised manuscript

25 November 2016)

A radial basis function ghost cell immersed boundary method of simulating flows around arbitrary complex or multiple immersed boundaries is proposed in this paper.In this method,incompressible Navier-Stokes equations are discretized on fixed Cartesian staggered gridby the finite difference method.A fractional step method is used for time integration,together with third order Runge-Kutta scheme.A high-order TVD MUSCL(total variation diminishing monotonic upstream-centered scheme for conservation law)scheme is used to discretize convective terms.Two salient features are emphasized in the present study.First,boundary conditions at the immersed interface are enforced by a continuous ghost cell method to consider the influence of immersed boundary on the flow field.The immersed bodies are treated as virtual boundaries immersed in the flow.And Navier-Stokes equations are solved in the entire computation domain,including solid domain.Therefore,programming complexity is greatly reduced and the treatment of immersed boundaries is simplified.Second,a polynomial and radial basis function is introduced to implicitly represent and reconstruct arbitrary complex immersed boundaries.Iso-surface distance functions about interface geometries are fitted with some sampling points of body surfaces.It is flexible and robust.Moreover,the information about interface positions on the background grid can be easily identified by the signed distance functions.Based on our in-house developed immersed boundary method solver,typical test cases are simulated to validate the proposed method.Theflows around a cylinder at Reynolds numbers of 40,100 and 200 are first simulated and a grid resolution study is carried out.Good agreement is achieved by comparing with previous numerical results,which shows that this method is accurate and reliable.In the second case of flow around airfoil,the good agreement with previous study shows that the present method has the ability to simulate complex immersed boundary flow.In the last case of flow around array of thirteen cylinders,the ability of present method for multiple immersed boundaries is well proved.And hydrodynamic interaction among multiple bodies is briefly analyzed.

immersed boundary method,radial basis function,ghost cell method,complex boundary

:47.85.Dh,47.11.Bc,45.50.Jf

10.7498/aps.66.044704

?通信作者.E-mail：xinjinjin1990@sina.com

?Corresponding author.E-mail：xinjinjin1990@sina.com