辛建建 石伏龍 金秋

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

一種徑向基函數(shù)虛擬網(wǎng)格法數(shù)值模擬復雜邊界流動

辛建建?石伏龍 金秋

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

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

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

浸入邊界法,徑向基函數(shù),虛擬網(wǎng)格法,復雜邊界

1 引 言

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

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

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

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

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

2 數(shù)學模型和數(shù)值方法

2.1 控制方程

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

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

2.2 數(shù)值離散方法

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

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

2.3 虛擬網(wǎng)格浸入邊界法

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

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

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

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

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

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

[φ]T為1,2,3,4節(jié)點變量值;[B]為雙線性內(nèi)插格式的系數(shù)矩陣,表達式為

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

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

圖3 鄰近的網(wǎng)格節(jié)點表示Fig.3.Stencil for adjacent grid nodes.

3)計算虛擬網(wǎng)格值：鏡像點值φIP計算后,沿著物面法向采用線性內(nèi)插方法就可計算出虛擬網(wǎng)格節(jié)點值φGC：

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

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

2.4 PRBF界面描述方法

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

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

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

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

多項式PF(xi)為

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

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

進一步構(gòu)造如下等值面函數(shù)ψ(x),以有效識別網(wǎng)格節(jié)點ψ(x)在場中屬性：

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

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

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

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

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

3 算例與結(jié)果分析

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

3.1 單圓柱的二維黏性流動

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

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

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

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

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

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

3.2 翼型繞流

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

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

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

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

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

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

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

圖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時的阻力系數(shù)比較Fig.12.Comparison of drag coefficients of cylinder with center coordinates(0,0)atRe=100,200.

4 結(jié) 論

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

[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全國水動力學學術(shù)會議舟山,中國,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