劉肖肖

(南京航空航天大學(xué) 理學(xué)院，南京 211106)

?

多介質(zhì)球(軸)對(duì)稱問(wèn)題的高精度數(shù)值模擬方法

劉肖肖

(南京航空航天大學(xué) 理學(xué)院，南京 211106)

給出了球?qū)ΨQ和軸對(duì)稱多介質(zhì)流動(dòng)問(wèn)題的高精度數(shù)值模擬方法.在界面處利用RGFM(Real Ghost Fluid Method)方法跟蹤界面并定義界面邊界條件,將多介質(zhì)流動(dòng)問(wèn)題轉(zhuǎn)化為單介質(zhì)流動(dòng)問(wèn)題,并利用高精度RKDG(Runge-Kutta Discontinuous Galerkin)法進(jìn)行計(jì)算.對(duì)于水下爆炸等問(wèn)題中出現(xiàn)的空化現(xiàn)象,采用等熵的單流體空化模型來(lái)處理.數(shù)值試驗(yàn)結(jié)果表明,該方法能夠有效捕捉各種物理現(xiàn)象,保持界面清晰,且大大提高計(jì)算效率.

多介質(zhì)；對(duì)稱；RKDG；RGFM；空化

Cockburn和Shu等構(gòu)造了求解一維守恒律方程的高階RKDG方法[1-2],并且將該方法成功地推廣到多維標(biāo)量守恒律[3]和多維守恒律方程組[4].RKDG方法具有高精度和緊性,易于處理復(fù)雜的邊界,被廣泛應(yīng)用于流體力學(xué)問(wèn)題的數(shù)值模擬.

可壓縮多介質(zhì)流動(dòng)問(wèn)題是計(jì)算流體力學(xué)研究的重要內(nèi)容之一,該類問(wèn)題的數(shù)值模擬在流體力學(xué),航空航天,爆炸等領(lǐng)域應(yīng)用非常廣泛,尤其是涉及到激波與不同介質(zhì)交互界面的相互作用的空氣爆炸,氣泡變形與坍塌以及水下爆炸等問(wèn)題.在水下爆炸等問(wèn)題中,由于低壓區(qū)有時(shí)會(huì)出現(xiàn)空化現(xiàn)象,且空化氣泡破裂對(duì)船舶,水下艦艇會(huì)產(chǎn)生非常大的破壞作用,因此準(zhǔn)確的模擬空化問(wèn)題也具有非常重要的意義.

對(duì)于多介質(zhì)流動(dòng)問(wèn)題的數(shù)值模擬,由于在界面兩邊流體的狀態(tài)方程不同,流體的密度會(huì)出現(xiàn)大梯度變化,若采用單介質(zhì)流動(dòng)問(wèn)題的數(shù)值方法,數(shù)值不穩(wěn)定情況就會(huì)出現(xiàn),在界面處會(huì)出現(xiàn)非物理振蕩,甚至使得計(jì)算難以進(jìn)行.Fedkiw等提出的GFM(Ghost Fluid Method)方法[5]是一個(gè)簡(jiǎn)單有效的方法,Liu等在此基礎(chǔ)上構(gòu)造了更加穩(wěn)健的MGFM(Modified GFM)方法[6].本文利用RGFM[7]來(lái)處理界面,通過(guò)求解Level Set方程跟蹤界面的位置,在界面處利用Riemann問(wèn)題的解定義界面邊界條件,將多介質(zhì)問(wèn)題轉(zhuǎn)化為單介質(zhì)問(wèn)題,并結(jié)合高精度RKDG方法進(jìn)行求解.

1 方程

1.1 控制方程

球?qū)ΨQ歐拉方程組為：

(1)

軸對(duì)稱歐拉方程組為：

(2)

為了使歐拉方程組封閉,需要加入如下狀態(tài)方程：

理想氣體的狀態(tài)方程：

(3)

水的狀態(tài)方程：

p=(N-1)ρe-NB

(4)

其中γ、N、β均為正常數(shù),將在具體算例中給出.

1.2 空化區(qū)域混合流體的狀態(tài)方程

水中因壓力降低會(huì)產(chǎn)生的空化現(xiàn)象,對(duì)于空化區(qū)域的混合流體,狀態(tài)方程采用Liu等提出的等熵單流體空化模型[8].

混合流體的狀態(tài)方程：

(5)

1.3 Level Set方程

本文通過(guò)求解Level Set方程跟蹤界面的位置,二維Level Set方程為

(6)

其中φ(x,y,t)為符號(hào)距離函數(shù),u和v分別為徑向和z方向的速度.利用求解Hamilton-Jacobi方程的三階RKDG方法求解方程(6),可以得到準(zhǔn)確的界面位置.

2 數(shù)值計(jì)算方法

為簡(jiǎn)單起見,本文只給出一維含源項(xiàng)歐拉方程組(1)的離散方法.

首先將計(jì)算區(qū)域[a,b]劃分成N個(gè)網(wǎng)格,a=r1/2

(7)

(8)

(9)

(10)

(11)

(12)

Uh(x,t)在網(wǎng)格Ij左右邊界處的值可表示為：

(13)

(14)

其中TVD限制器函數(shù)為：

(15)

對(duì)于方程組(1),(2),Uh為向量函數(shù),可將Uh映射到特征空間,對(duì)其每一個(gè)分量重復(fù)式(13)～(14)即可.

三維歐拉方程組經(jīng)坐標(biāo)變換后得到方程組(1)和(2),可見r=0處為奇點(diǎn),直接用數(shù)值方法求解方程組,會(huì)導(dǎo)致數(shù)值不穩(wěn)定．為了保證所采用的數(shù)值格式與方程組(1),(2)的相容性,在奇點(diǎn)附近采用如(16)的守恒型方程組

(16a)

(16b)

在奇點(diǎn)處u=0

(16c)

α=2時(shí)為軸對(duì)稱問(wèn)題在任意z處奇點(diǎn)附近的流動(dòng)控制方程,α=3為球?qū)ΨQ問(wèn)題在奇點(diǎn)附近的流動(dòng)控制方程.

對(duì)于方程組(16)采用如下差分格式離散：

(17a)

(17b)

(17c)

(17d)

其中λr=Δt/Δr,下標(biāo)“1”和“2”分別表示變量在r1=Δr和r2=2Δr處的值，“0”表示變量在r=0處的值.

3 數(shù)值實(shí)驗(yàn)

將二階(k=1)精度RKDG方法中的方程(12)結(jié)合一維和二維RGFM界面處理方法應(yīng)用于多介質(zhì)球(軸)對(duì)稱水下爆炸問(wèn)題的數(shù)值模擬.對(duì)于一維方程組(1),r方向采用均勻步長(zhǎng),根據(jù)RKDG方法的穩(wěn)定性條件,時(shí)間步長(zhǎng)取為

對(duì)于二維方程組(2)，r方向和z方向均采用均勻步長(zhǎng),時(shí)間步長(zhǎng)取為

其中CFL=0.3,c為音速.

例1 深水高壓氣泡爆炸問(wèn)題

本算例是一個(gè)深水高壓氣泡爆炸問(wèn)題,氣泡內(nèi)是高壓氣體,氣泡外是水,氣泡的半徑為1,氣體的無(wú)量綱初始狀態(tài)值為uH=0.0,pH=8 290,γ=1.4,ρH=1.27,水的無(wú)量綱初始狀態(tài)值為ρl=1.0,ul=0.0,pl=1.0,N=7.15,A=1.0,B=3 310,計(jì)算時(shí)刻t=0.004 78時(shí)的計(jì)算區(qū)域?yàn)閇0,4],網(wǎng)格數(shù)為400,t=0.015和t=0.021 2時(shí)的計(jì)算區(qū)域時(shí)[0,6],網(wǎng)格數(shù)為600,t=0.04時(shí)的計(jì)算區(qū)域?yàn)閇0,10],網(wǎng)格數(shù)是1 000.圖1分別給出了密度在不同時(shí)刻的分布圖,并給出了半徑隨時(shí)間的變化曲線(圖2).

圖1 密度在半徑方向的分布圖

圖2 半徑隨時(shí)間的分布圖

例2 圓柱形容器內(nèi)水下爆炸問(wèn)題

本算例是容器內(nèi)氣泡爆炸的軸對(duì)稱問(wèn)題,圓柱形容器的半徑為0.044 45 m,高度為0.228 6 m,在容器中心有一個(gè)高壓氣泡,氣泡半徑為0.015 m,氣泡外是水,氣泡內(nèi)氣體的初始狀態(tài)值為ρH=1 770 kg/m3,uH=0.0 m/s,vH=0.0 m/s,pH=2.0×109Pa,γ=2.0,氣泡外水的初始狀態(tài)值為ρl=1 000 kg/m3,N=7.15,ul=0.0 m/s,vl=0.0 m/s,pl=105Pa,A=105Pa,B=3 310×105Pa,計(jì)算區(qū)域?yàn)閞×z∈[0,0.044 45]×[-0.114 3,0.114 3],網(wǎng)格數(shù)為70×360,在容器壁處采用反射邊界條件,圖3分別給出了不同時(shí)刻壓力的等值線圖.

(a) t=10 μs (b) t=30 μs (c) t=50 μs圖3 不同時(shí)刻壓力的等值線圖

例3 近水面水下爆炸問(wèn)題

此問(wèn)題是一個(gè)軸對(duì)稱問(wèn)題，用固壁代替自由面，計(jì)算區(qū)域是[0,9]×[-9,3]，網(wǎng)格數(shù)是270×360,y=3處為固壁,高壓氣泡位于(0,0)點(diǎn),半徑為1,固壁處采用反射邊界,其他邊界為自由邊界.氣泡的初始狀態(tài)值為ρH=1 270 kg/m3,uH=0.0 m/s,vH=0.0 m/s,pH=8.29×108Pa,γ=2.0,周圍水的初始狀態(tài)值為ρl=1 000 kg/m3,ul=0.0 m/s,vl=0.0 m/s,pl=105Pa,N=7.15,A=105Pa,B=3 310×105Pa.圖4～圖9分別給出了不同時(shí)刻壓力,r方向和z方向的速度等值線圖.

圖4 t=2.0 ms時(shí)刻壓力的等值線圖

圖5 t=2.0 ms時(shí)刻r方向速度的等值線圖

圖6 t=2.0 ms時(shí)刻z方向速度的等值線圖

圖7 t=4.0 ms時(shí)刻壓力的等值線圖

圖8 t=4.0 ms時(shí)刻r方向速度的等值線圖

圖9 t=4.0 ms時(shí)刻z方向速度的等值線圖

4 結(jié)論

本文給出了可壓縮多介質(zhì)流動(dòng)球?qū)ΨQ和軸對(duì)稱問(wèn)題的數(shù)值模擬方法,利用RGFM方法將多介質(zhì)問(wèn)題轉(zhuǎn)化為單介質(zhì)問(wèn)題,并利用RKDG方法求解歐拉方程.對(duì)于水下爆炸問(wèn)題中出現(xiàn)的空化現(xiàn)象,采用等熵的單流體空化模型進(jìn)行處理.數(shù)值試驗(yàn)結(jié)果表明,該算法可大大減少計(jì)算量,有效捕捉各種物理現(xiàn)象.

[1] COCKBURN B,SHU C W.TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws.II.General framework[J].Mathematics of Computation,1989,52(186)：411-435.

[2] COCKBURN B,LIN S Y,SHU C W.TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III：one-dimensional systems[J].Journal of Computational Physics,1989,84(1)：90-113.

[3] COCKBURN B,HOU S,SHU C W.The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws.IV.The multidimensional case[J].Mathematics of Computation,1990,54(190)：545-581.

[4] COCKBURN B,SHU C W.The Runge-Kutta discontinuous Galerkin method for conservation laws V：multidimensional systems[J].Journal of Computational Physics,1998,141(2)：199-224.

[5] FEDKIW R P,ASLAM T,MERRIMAN B,et al.A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[J].Journal of Computational Physics,1999,152(2)：457-492.

[6] LIU T G,KHOO B C,YEO K S.Ghost fluid method for strong shock impacting on material interface[J].Journal of Computational Physics,2003,190(2)：651-681.

[7] WANG C W,LIU T G,KHOO B C.A real ghost fluid method for the simulation of multimedium compressible flow[J].SIAM Journal on Scientific Computing,2006,28(1)：278-302.

[8] LIU T G,KHOO B C,XIE W F.Isentropic one-fluid modelling of unsteady cavitating flow[J].Journal of Computational Physics,2004,201(1)：80-108.

[責(zé)任編輯 王新奇]

A High Precision Numerical Simulation Method for theSymmetric Problem of Multi Media Balls (Shaft)

LIU Xiao-xiao

(School of science, Nanjing University of Aeronautics & Astronautics, Nanjing 211106, China)

In this paper, a high precision numerical simulation method for the flow problem of spherical symmetric and axisymmetric multimedia is presented. At the interface, the RGFM (Real Ghost Fluid Method) method is used to track the interface and define the boundary conditions. The multi medium flow problem is transformed into the single medium flow problem, and the calculation is carried out by using the high precision RKDG (Runge-Kutta Discontinuous Galerkin) method. The single fluid cavitation model with equal entropy is used to deal with the problems of the phenomenon of cavitation during the underwater explosion. The numerical experiment results show that the method can effectively capture all kinds of physical phenomena, keep the interface clear, and greatly improve the computational efficiency.

multimedia; symmetry; RKDG (Runge-Kutta Discontinuous Galerkin); RGFM (Real Ghost Fluid Method); cavitation

1008-5564(2016)04-0011-06

2016-01-15

劉肖肖(1988—)，女，山東濟(jì)寧人，南京航空航天大學(xué)理學(xué)院碩士研究生，主要從事計(jì)算流體力學(xué)研究.

O351

A