基于一種改進的虛擬單元法模擬包含靜止/運動邊界的流動問題*

魏峰 金亮? 柳軍 丁峰 鄭新萍

1)(國防科技大學(xué)空天科學(xué)學(xué)院,長沙 410073)

2)(國防科技大學(xué)計算機學(xué)院,長沙 410073)

1 引 言

在計算流體力學(xué)領(lǐng)域,笛卡爾網(wǎng)格相比貼體網(wǎng)格有著高效率、高魯棒性、高靈活性等突出優(yōu)勢,能夠很好地應(yīng)對復(fù)雜的流動問題[1,2].作為笛卡爾網(wǎng)格技術(shù)的關(guān)鍵,物面處理的方式一直是國內(nèi)外學(xué)者們研究的熱點,近年來,浸入邊界法(immersed boundary method)[3]越來越受到研究者們的重視,這種方法相比切割單元法(cut-cell method)[4]不需要對網(wǎng)格進行切割,極大地減小了計算量,在模擬包含運動邊界的流動問題時有著明顯的優(yōu)勢.常規(guī)的浸入邊界法包括連續(xù)力法和離散力法,這兩種方法在模擬不可壓縮流的領(lǐng)域都獲得了極大的成功[5,6],然而在模擬高雷諾數(shù)情況下的可壓縮流動時都存在著很大的誤差.為了彌補這一缺陷,Majumdar等[7]提出了虛擬單元法(ghost cell method),即在物體內(nèi)部構(gòu)建虛擬流場,通過對虛擬流場內(nèi)的虛擬單元進行合理的賦值,隱式地表達物面邊界條件.一般情況下,虛擬單元定義為周圍至少存在一個流場單元的物體內(nèi)部單元,但是對于一些高階的虛擬單元法,則需要擴展更多的虛擬單元[8],對于存在細(xì)長結(jié)構(gòu)的幾何模型,虛擬單元還可以定義在流場內(nèi)部[9].虛擬單元法以其靈活簡單、魯棒性好等優(yōu)點得到了眾多學(xué)者的關(guān)注.Tseng和Ferziger[10]對浸入邊界虛擬單元法在各種邊界條件下的應(yīng)用進行了系統(tǒng)的研究;Dadone和Grossman[11],Farooq等[12]以及Chi等[13]都對可壓縮流動中的虛擬單元法進行了詳細(xì)研究并發(fā)展了各自的虛擬單元法.

對于含有運動邊界的復(fù)雜流動問題,早期的物面邊界處理方式為切割單元法,容易導(dǎo)致所謂的“細(xì)小網(wǎng)格”問題,嚴(yán)重限制了時間步長.為了避免這一問題,Mittal 和Iaccarino[5]將細(xì)小的切割單元與周圍的流場單元進行融合以消除細(xì)小單元,Forrer 和Berger[14]借鑒了浸入邊界法的思想,將切割單元進行了填充并且設(shè)置了一系列的虛擬單元,兩者都成功地求解了含有運動邊界的流動問題.Lee 等分析研究了求解運動邊界問題時產(chǎn)生的壓力振蕩,通過實驗發(fā)現(xiàn)減小網(wǎng)格尺寸或者增大時間步長都可以減小壓力振蕩[15],進而提出了一種能夠有效地控制壓力振蕩的隱式浸入邊界虛擬單元法[16].Tan和Shu[8]推導(dǎo)了一種復(fù)雜的能夠處理運動邊界問題的高階虛擬單元法,通過有限差分的方法分別求解了一維算例和二維算例,其精度分別達到了四階和三階.Peter和De[17]在非慣性系下對GCIBM (ghost-cell immersed boundary method)方法進行了拓展,同時能夠?qū)崿F(xiàn)并行求解,相比在慣性系下求解包含運動邊界的流動問題,避免了壓力的數(shù)值振蕩,但是也存在著無法模擬多體相對運動的局限性.在國內(nèi),趙寧等[18,19]在自適應(yīng)笛卡爾網(wǎng)格的基礎(chǔ)上對虛擬單元法進行了大量的研究,同時也對包含運動邊界的情況進行了數(shù)值模擬.

在文獻[12]提出的SGCM (simplified ghost cell method)基礎(chǔ)上,本文提出了一種改進的虛擬單元法,同樣使用了流場點作為鏡像點,但是相比之下具有更高的準(zhǔn)確性和計算效率.考慮到物面邊界變速運動會導(dǎo)致邊界附近的壓力、密度和切向速度存在梯度,給出了改進型虛擬單元法的推廣形式,用于求解包含運動邊界的流動問題.通過求解Schardin問題和激波抬升輕質(zhì)圓柱問題驗證了該虛擬單元法及其推廣形式在處理包含靜止/運動邊界的流動問題時的準(zhǔn)確性.本文第2節(jié)簡要介紹了求解過程中使用的控制方程以及空間、時間離散格式,進而詳細(xì)介紹了改進的虛擬單元法及其推廣;第3節(jié)給出了兩種經(jīng)典的流動問題的計算結(jié)果和討論;第4節(jié)給出了結(jié)論.

2 數(shù)值方法

2.1 控制方程

二維無黏可壓縮歐拉方程的有限體積形式如下[20]：

其中,

U為守恒變量;F(U)和G(U)分別為X和Y方向上的通量;ρ為密度;u和v分別為沿X和Y方向的速度分量;p為壓強;E為總能.考慮到氣體的熱力學(xué)性質(zhì),補充方程(5)封閉上述方程系統(tǒng),

其中比熱比γ取1.4.

2.2 離散

采用AUSM + (advection upstream splitting method +)格式計算對流通量[21],為了提高計算精度,通過MUSCL (monotonic upstream-centered scheme for conservation laws)方法[22]在單元交界處采用二階精度重構(gòu)格式對變量進行重構(gòu),同時選用MINMOD限制器限制重構(gòu)時使用的梯度,防止重構(gòu)格式出現(xiàn)振蕩.

時間離散采用顯式三階TVD (total-variationdiminishing)Runge-Kutta格式[23],具體形式如下：

其中R(*)為殘值,Δt為時間步長.對于非定常問題,應(yīng)采用全場統(tǒng)一時間步長,全場統(tǒng)一時間步長是所有當(dāng)?shù)貢r間步長的最小值,即

式中Δti為當(dāng)?shù)貢r間步長.

2.3 改進的虛擬單元法

在使用常規(guī)的虛擬單元法求解可壓縮流動問題時,通常需要得到虛擬單元關(guān)于物面的鏡像點,然后通過插值得到鏡像點的變量值,進而求得虛擬單元的變量值.一般情況下,鏡像點周圍都有著足夠的插值點,因此能夠采用各種插值方法得到比較準(zhǔn)確的數(shù)值,但是對于圖1中出現(xiàn)的情況,物體底部虛擬單元的鏡像點缺乏合適的插值點,因此需要進行特殊的處理.一般地,可以將距離鏡像點最近流場點的變量值賦給鏡像點,或者是距離虛擬單元最近流場點的變量值賦給鏡像點.然而,這種近似的賦值方式會引入一定的誤差,對于誤差存在積累效應(yīng)的非定常問題,如果一開始就存在較大的誤差將嚴(yán)重影響最終結(jié)果的準(zhǔn)確性.

圖1 狹縫問題(物體底部虛擬單元的鏡像點缺乏合適的插值點)Fig.1.Slit problem,which is caused by the lack of interpolation points of the mirror points.

為了不采用這種近似的處理方法,與Farooq等[12]在2013年提出的SGCM相似,選擇了X或者Y方向的流場點代替了鏡像點,避免了插值運算,因此能夠很好地處理狹縫問題.考慮到鏡像點與物面點十分接近時,速度賦值公式存在分母接近0的情況而產(chǎn)生非物理解,在SGCM中假設(shè)了物面點位于虛擬單元和鏡像點中間,但同時也引入了一定的誤差.如圖2所示,使用了精確的物面點位置求得虛擬單元的變量值,從而能夠精確地表達物面邊界,同時為了保證方法的穩(wěn)定性,在計算過程設(shè)定當(dāng)鏡像點與物面點距離小于0.2個網(wǎng)格單元長度時,就要取鏡像點向外延伸后的第二個流場點作為鏡像點,例如,圖2中虛擬單元B的鏡像點本來應(yīng)該是流場點N,但是由于N點與物面點MB的距離過近而選擇了N點上方的流場點FB.值得注意的是,對于這種改進的虛擬單元法,如果虛擬單元距離在X方向或者Y方向上的鏡像點過遠,則將只考慮距離虛擬單元最近的鏡像點,例如圖中的虛擬單元B沿X方向的第一個鏡像點為流場點,這兩點的距離達到了3倍的網(wǎng)格長度,因此在對虛擬單元B賦值時只需考慮其沿Y方向的鏡像點即可,在實際計算時,當(dāng)虛擬單元與第一個鏡像點之間距離為1到2個網(wǎng)格長度時需要考慮兩個方向的鏡像點,如虛擬單元A和C.對于非曲面物面,虛擬單元的賦值方法如下所示：

圖2 改進的虛擬單元法示意圖Fig.2.Demonstration of the improved ghost cell method.

其中Q代表壓強p、密度ρ以及沿著物面的切向速度ut;un代表物面法向速度;a和b分別代表鏡像點到物面點的距離和虛擬單元到鏡像點的距離.對于曲面物面,在計算除un以外的其他變量時不能忽略曲率半徑的影響,因此在計算這些變量時需要對(10)式進行修改,本文采用了與文獻[12]相同的壓力賦值公式.在SGCM中,鏡像點的選擇(X方向上選擇或者Y方向上選擇)依據(jù)物面的法向向量,因此在鏡像點選擇更替的區(qū)域,流場變量會可能出現(xiàn)不連續(xù),如圖3所示,在使用SGCM方法計算超聲速圓柱繞流時,在r/R=0.7附近的壓力系數(shù)會出現(xiàn)扭曲現(xiàn)象.為了解決這個問題,Farooq等[12]在鏡像點的選擇中增加了斜上(下)方的流場點,得到了MSGCM (modified simplified ghost cell method),消除了扭曲,但是使得代碼中的判斷語句變得非常多.我們則根據(jù)距離加權(quán)的方式來光滑變量值,即同時得到X方向以及Y方向的鏡像點,通過(10)和(11)式分別求出相應(yīng)的虛擬單元變量值,最終的虛擬單元變量值根據(jù)它們距離物面的長短進行加權(quán),這一過程并沒有引入過多的判斷語句.在相同的設(shè)置參數(shù)和網(wǎng)格量下我們求解了文獻[12]中的超聲速圓柱繞流算例,計算結(jié)果如圖3所示,可以看出本文發(fā)展的改進型虛擬單元法(ISGCM)同樣能夠消除計算過程中出現(xiàn)的扭曲現(xiàn)象,并且比SGCM和MSGCM更加精確.

圖3 SGCM中出現(xiàn)的扭曲現(xiàn)象(標(biāo)記區(qū)域)以及采用ISGCM求得的圓柱表面壓力系數(shù)與采用SGCM、MSGCM、貼體網(wǎng)格得到的結(jié)果的對比Fig.3.Kink noted in the picture,which is happened when using the SGCM,and the comparison of the pressure coefficients obtained by using SGCM,MSGCM,ISGCM and the body fitted mesh.

2.4 推廣的改進型虛擬單元法

對于變速度的運動邊界問題,邊界附近存在著壓力、密度和切向速度的梯度[12],因此在計算虛擬單元的變量值時我們額外選取了一層鏡像點(圖4),增加了FA,2,以及FB,2,結(jié)合對應(yīng)的第一層鏡像點FA,1,和FB,1得到近壁面的變量梯度,進而求解虛擬單元的變量值.(12)和(13)式給出了這種推廣形式的虛擬單元賦值方法：

在對虛擬單元賦值時,除了上文中提到的特殊情況,都需要分別得到沿X和Y方向求得的虛擬單元變量值,再根據(jù)虛擬單元與物面點之間的距離進行加權(quán).當(dāng)物面附近存在激波時,例如一道激波介于流場點FA,1和FA,2之間,采用(12)式得到的壓力等變量會出現(xiàn)很大或者很小的數(shù)值,甚至出現(xiàn)負(fù)值,導(dǎo)致計算發(fā)散.因此我們增加了一個激波監(jiān)測器α,用來檢測FA,1和FA,2之間是否激波,其定義如下：

圖4 推廣的改進型虛擬單元法示意圖Fig.4.The demonstration of the extended ISGCM.

3 計算結(jié)果與討論

為了驗證本文發(fā)展的改進型虛擬單元法,首先求解了Schardin問題,并與相關(guān)文獻進行對比,驗證了方法在求解包含靜止邊界的非定常流動問題時的準(zhǔn)確性,進而求解了存在狹縫問題的激波抬升輕質(zhì)圓柱問題,通過和相關(guān)文獻進行對比,驗證了在包含運動邊界情況下方法的準(zhǔn)確性以及對狹縫問題的處理效果.

3.1 Schardin問題

Schardin問題于1957年由Schardin提出[24],此后,一些學(xué)者對該問題進行了實驗和數(shù)值模擬,本算例選用了Chang等[25]提出的研究模型,模型描述如圖5所示,初始時刻(T=0),距離左邊界50 mm處有一道激波,隨著時間的推移,激波將以Ma=1.34向后掃過一個夾角約為60°的三角楔,流動過程中將出現(xiàn)馬赫桿、三波點、滑移線、膨脹波等復(fù)雜的流場結(jié)構(gòu),初始時低壓區(qū)的靜壓為0.05 MPa,更多的細(xì)節(jié)可以查看文獻[25].左邊界賦值為波后條件,右邊界以及上下邊界為超聲速出口條件,物面邊界采用本文發(fā)展的ISGCM,計算網(wǎng)格采用了均勻的笛卡爾網(wǎng)格,分別計算了網(wǎng)格量大小為142×100,426×300和710×500時的流動情況.

圖5 Chang等[25]用來研究Schardin問題的物理模型Fig.5.Physical model of the Schardin's problem proposed by Chang et al.[25].

圖6 不同時刻下的實驗結(jié)果(左)、Chang等[25]計算所得的密度等值線(中)和采用ISGCM計算所得的密度等值線(右)(a)T=28 μs;(b)T=53 μs;(c)T=102 μs;(d)T=130 μs;(e)T=172 μsFig.6.Experimental results (left),Chang et al.[25] results and the density contour computed by us (right)：(a)T=28 μs;(b)T=53 μs;(c)T=102 μs;(d)T=130 μs;(e)T=172 μs.

圖6分別給出了在網(wǎng)格量大小為710×500時采用ISGCM所得的密度等值線圖(右)、Chang等[25]使用數(shù)值計算得到的密度等值線圖(中)以及通過雙脈沖全息干涉法得到的T=28,53,102,130,172 μs時的實驗結(jié)果(左),其中文獻[25]在數(shù)值計算過程中在自適應(yīng)的四邊形貼體網(wǎng)格上使用MUSCL方法得到了二階精度的計算結(jié)果.圖7(a)給出了不同網(wǎng)格量下計算所得的沿模型中間對稱線上的密度分布與文獻[25]的數(shù)值計算結(jié)果以及部分實驗數(shù)據(jù)的對比,圖7(b)給出了不同網(wǎng)格量下計算所得的沿模型中間對稱線上的馬赫數(shù)分布與文獻[25]的數(shù)值計算結(jié)果的對比.可以看出,計算結(jié)果與文獻[25]的實驗以及計算結(jié)果基本一致,表明本文所發(fā)展的改進型虛擬單元法能夠達到與貼體網(wǎng)格相同的精度.

圖7 (a)計算所得沿模型中間對稱面上的密度分布與Chang等[25]的計算結(jié)果以及實驗數(shù)據(jù)的對比;(b)沿模型中間對稱面上的馬赫數(shù)分布與Chang等[25]的計算結(jié)果的對比Fig.7.(a)Comparison of the density distribution along the symmetry plane between the results obtained by this paper and Chang et al.[25] as well as the experiment;(b)the comparison of the Mach number distribution along the symmetry plane between the results obtained by this paper and the results of Chang et al.[25].

3.2 激波抬升輕質(zhì)圓柱

在此算例中,一道激波從左至右傳播,一個剛性輕質(zhì)圓柱放置在管道內(nèi),當(dāng)激波接觸圓柱后,會對圓柱產(chǎn)生作用力并將其抬升,文獻[8,14,26]都求解了該問題并進行了分析.本文采用文獻[8]的計算模型,模型示意圖見圖8,計算域的大小為[0,1]×[0,0.2],剛性輕質(zhì)圓柱的半徑為0.05,密度為10.77.初始時刻,圓柱中心點的位置為(0.15,0.05),一道激波在X=0.08處,并以馬赫數(shù)Ma=3的速度朝著X的正方向傳播.左邊界設(shè)置為波后條件,右邊界為超聲速出口邊界,上下邊界均為固壁.初始流場中低壓區(qū)的靜壓為1.0,密度為1.4.

上下邊界采用了對稱反射邊界條件,物面邊界采用了ISGCM的推廣形式.初始時刻,圓柱剛好與計算域的下邊界接觸,形成狹縫問題,設(shè)定在計算距離下邊界很近的虛擬單元變量值時只考慮X方向的鏡像點.值得注意的是,這里并沒有引入對鏡像點的近似賦值,仍然是對物面邊界的精確表達.表1列出了不同網(wǎng)格量下T=0.1641和0.30085時圓柱中心的坐標(biāo),計算結(jié)果與文獻[8]符合得較好.圖9和圖10分別給出了網(wǎng)格量為600×120時計算所得T=0.1641和0.30085時的壓力等值線圖(a)和網(wǎng)格量為640×128時文獻[8]的計算結(jié)果(b),其中文獻[8]在計算過程中使用了復(fù)雜的三階精度虛擬單元法.可以看出計算結(jié)果基本吻合,表明推廣的改進型虛擬單元法能夠用于求解包含運動邊界的流動問題.

圖8 激波抬升輕質(zhì)圓柱問題的模型描述Fig.8.Physical model of the cylinder lift-off problem.

圖9 T=0.1641時的壓力等值線圖 (a)本文所推廣的改進型虛擬單元法(網(wǎng)格量為600×120);(b)Tan等[8]使用的高階虛擬單元法(網(wǎng)格量為640×128)Fig.9.Pressure contour atT=0.1641：(a)Extended ISGCM with the grid size of 600×120;(b)high-order ghost cell method proposed by Tan et al.[8] with the grid size of 640×128.

圖10 T=0.30085時的壓力等值線圖 (a)本文所推廣的改進型虛擬單元法(網(wǎng)格量為600×120);(b)Tan等[8]使用的高階虛擬單元法(網(wǎng)格量為640×128)Fig.10.Pressure contour atT=0.30085：(a)Extended ISGCM with the grid size of 600×120;(b)high-order ghost cell method proposed by Tan et al.[8] with the grid size of 640×128.

表1 不同網(wǎng)格量下T=0.1641和0.30085時圓柱中心的坐標(biāo)Table 1.Position of the center of the cylinder atT=0.1641 and 0.30085.

4 結(jié) 論

本文提出了一種改進的虛擬單元法求解包含靜止和運動邊界的流動問題.在該方法中,虛擬單元的鏡像點用沿X和Y方向的流場點代替,虛擬單元的最終值通過其沿X和Y方向得到的變量值按照距離長短加權(quán)得到.考慮到物面邊界變速運動時不能忽略近壁面處的變量梯度,我們對該虛擬單元法進行了推廣,增加了一層鏡像點,使其能夠處理包含變速運動邊界的流動問題,同時為了避免激波貼體時產(chǎn)生非物理解,加入了激波監(jiān)視器,提高了方法的穩(wěn)定性.最后通過求解Schardin問題和激波抬升輕質(zhì)圓柱兩個經(jīng)典算例驗證了方法的可靠性和準(zhǔn)確性.本文的結(jié)論如下：

1)與傳統(tǒng)的虛擬單元法相比,本文的方法不需要插值得到鏡像點的變量值,保持了笛卡爾網(wǎng)格的簡單性,減少了計算量,同時由于鏡像點就是流場點,因而能夠很好地處理狹縫問題;

2)與SGCM相比,本文的方法使用了精確的物面點求解虛擬單元的變量值,能夠準(zhǔn)確地表達物面邊界條件,并且不存在扭曲現(xiàn)象,同時還避免了過多判斷語句的引入,具有更高的準(zhǔn)確性和計算效率;

3)通過求解Schardin問題驗證了本文的方法可以較好地處理包含靜止邊界的流動問題,同時能夠達到與貼體網(wǎng)格相同的精度;通過求解激波抬升輕質(zhì)圓柱問題,并與復(fù)雜的高階虛擬單元法進行對比,驗證了本文相對簡單的虛擬單元法能夠滿足處理包含運動邊界問題的要求.

感謝上海航天科技創(chuàng)新基金提供的資金支持,同時衷心感謝審稿人的建設(shè)性意見和建議.