球形重質(zhì)氣體物理爆炸特性*

薛大文，陳志華，韓珺禮,3

(1.南京理工大學(xué)瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210094；2.浙江海洋學(xué)院海運(yùn)與港航建筑工程學(xué)院, 浙江 舟山 316022；3.北京機(jī)電研究所，北京 100012)

高壓氣體在運(yùn)輸、儲(chǔ)存、加工和使用過(guò)程中常發(fā)生爆炸事故，造成巨大的人員傷亡和財(cái)產(chǎn)損失，同時(shí)，由于爆炸過(guò)程還蘊(yùn)含了許多諸如界面不穩(wěn)定性以及激波與界面相互作用等復(fù)雜流體物理現(xiàn)象，因此對(duì)氣云物理爆炸特性的研究不僅具有重要的實(shí)際意義，還具有重要的理論價(jià)值。

自從G.I.Taylor[1]提出爆炸波沖擊理論以來(lái)，不少學(xué)者做過(guò)類(lèi)似研究[2-3]。事實(shí)上，氣相爆炸中包含了極其豐富的物理現(xiàn)象，特別是高壓、高密度的氣體云爆炸，其包含了在激波作用下自由剪切層的發(fā)展、Richtmyer-Meshkov (RM)不穩(wěn)定性以及湍流的轉(zhuǎn)捩、氣動(dòng)噪聲等基本的流體物理現(xiàn)象。其中RM不穩(wěn)定性是指不同密度的流體分界面在運(yùn)動(dòng)激波的作用下，在界面附近形成不斷擴(kuò)展的湍流混合層，這種復(fù)雜的非線性現(xiàn)象在超燃沖壓發(fā)動(dòng)機(jī)的混合燃燒[4]、爆轟[5-6]、慣性約束聚變[7]以及天體物理的超新星爆發(fā)等問(wèn)題中起重要作用，因而RM失穩(wěn)現(xiàn)象得到了廣泛的關(guān)注，研究人員做了許多激波與密度界面相互作用的研究[8-9]。C.A.Zoldi等[10]在激波管上開(kāi)展了激波加載SF6氣柱RM界面不穩(wěn)定性實(shí)驗(yàn)，并采用RAGE程序模擬了SF6氣柱在空氣沖擊波作用下界面的演化、發(fā)展過(guò)程以及后期空氣和SF6氣體的混合。G.Layes等[11]利用高速攝影方法研究了激波對(duì)不同密度氣泡作用后的變形過(guò)程。鄒立勇等[12]利用高速攝影測(cè)試技術(shù)，對(duì)弱激波沖擊下空氣中的SF6氣柱和氣簾界面的演變過(guò)程進(jìn)行了實(shí)驗(yàn)研究，對(duì)演變過(guò)程中渦的產(chǎn)生和發(fā)展進(jìn)行了初步解釋。而事實(shí)上，對(duì)于氣相爆炸，人們關(guān)注更多的是由球形激波引起的RM失穩(wěn)。J.G.Zheng等[13]采用數(shù)值模擬的方法研究了由內(nèi)聚球形激波引起的RM失穩(wěn)現(xiàn)象，他們的研究顯示了平面激波和球形激波作用的不同特征。

柱形或球形氣云爆炸所產(chǎn)生的RM不穩(wěn)定性由于氣云內(nèi)部反射激波的作用變得更復(fù)雜，而且其不穩(wěn)定性后期非線性演化可加速可壓氣體的湍流轉(zhuǎn)捩。本文中,采用大渦模擬(large eddy simulation, LES)方法以及高精度混合格式對(duì)高壓、高密度的SF6氣體云冷爆炸進(jìn)行數(shù)值模擬，以研究其爆炸過(guò)程中RM失穩(wěn)以及二次激波與其作用后期非線性演化的湍流過(guò)程。

1 控制方程和數(shù)值方法

三維可壓Navier-Stokes(N-S)方程經(jīng)過(guò)Favre濾波后，其標(biāo)準(zhǔn)的可壓三維LES方程作為控制方程。連續(xù)性方程、動(dòng)量方程和能量方程分別為：

(1)

(2)

(3)

在大渦模擬中，亞格子湍流應(yīng)力項(xiàng)無(wú)法直接從可解尺度中計(jì)算求得,因此亞格子應(yīng)力需要引入模型進(jìn)行封閉。這里采用拉伸渦亞格子模型[14]對(duì)小尺度輸運(yùn)項(xiàng)建模。拉伸渦亞格子模型由D.I.Pullin[14]提出，他假設(shè)亞網(wǎng)格內(nèi)的流動(dòng)是由亞網(wǎng)格內(nèi)的渦導(dǎo)致的，并通過(guò)對(duì)稱渦結(jié)構(gòu)來(lái)構(gòu)建小尺度渦量。拉伸渦亞格子模型已被成功用于亞格子能量譜的估算以及湍流混合模擬。

以上LES方程可用有限體積法進(jìn)行數(shù)值求解。由于爆炸過(guò)程包括激波與氣體相互作用加速誘導(dǎo)湍流轉(zhuǎn)捩，因此對(duì)方程中的對(duì)流項(xiàng)離散要求很高，需數(shù)值格式同時(shí)具備對(duì)激波與湍流的高分辨率。一般采取高階的中心差分與迎風(fēng)型格式的混合形式，本文中采用TCD(tuned centered-difference)/WENO(weighted essentially nonoscillatory )混合數(shù)值方法[15]。時(shí)間推進(jìn)采用三階精度的Runge-Kutta法。

2 計(jì)算模型

在爆炸初期，流場(chǎng)結(jié)構(gòu)呈現(xiàn)對(duì)稱情形，而由于RM失穩(wěn)，后期的流場(chǎng)漸漸轉(zhuǎn)為非對(duì)稱結(jié)構(gòu)。為描述RM失穩(wěn)以及后期二次激波的作用過(guò)程，選擇初始狀態(tài)為球形氣云的SF6氣體作為研究對(duì)象。SF6氣體和空氣的密度分別為5.97和1.18 kg/m3，兩者的比熱比和氣體常量分別為γ(SF6)=1.09,R(SF6)=56.92;γ(air)=1.40,R(air)=287。選擇標(biāo)準(zhǔn)狀態(tài)的空氣參量，p0=101.3 kPa，T0= 287 K。SF6氣體的初始?jí)毫?0p0,SF6氣云初始半徑為R0=0.1 m，爆炸空間徑向大小為15R0。計(jì)算采用笛卡爾網(wǎng)格模擬氣云球體表面，以該網(wǎng)格擾動(dòng)作為初始隨機(jī)擾動(dòng)，球體表面網(wǎng)格總量約為800萬(wàn)，整個(gè)計(jì)算域網(wǎng)格總量約為4 500萬(wàn)。

3 結(jié)果與討論

圖1 不同時(shí)刻密度等值面Fig.1 Density isosurfaces shown at different times

圖2 流場(chǎng)波系結(jié)構(gòu)及界面隨時(shí)間的變化Fig.2 Evolution of the shock wave structure and interface in the flowfield

當(dāng)激波從重質(zhì)氣體傳向輕質(zhì)氣體時(shí)，會(huì)產(chǎn)生向內(nèi)傳播的稀疏波和向外傳播的透射激波，從而引起失穩(wěn)，加速界面混合。激波作用不同密度氣體界面后，界面演化一般經(jīng)歷3個(gè)階段[16-17]：(1)振幅線性增長(zhǎng)階段；(2)非線性增長(zhǎng)階段，輕質(zhì)流體發(fā)展成“氣泡”，重質(zhì)流體發(fā)展成“尖釘”；(3)湍流混合階段，“尖釘”開(kāi)始破碎，不同尺度的渦相互吞并，不同流體達(dá)到湍流混合狀態(tài)。從圖1中可以清晰地看到氣體界面的3種狀態(tài)：在膨脹初期，隨機(jī)擾動(dòng)在氣體界面發(fā)生，擾動(dòng)振幅在界面線性增長(zhǎng)，見(jiàn)圖1(a);隨著膨脹的加劇，振幅經(jīng)歷非線性增長(zhǎng)，“氣泡”與“尖釘”結(jié)構(gòu)出現(xiàn),見(jiàn)圖1(b);隨著膨脹的結(jié)束以及球心處低壓區(qū)的形成，圓球界面開(kāi)始向球心回流，“尖釘”結(jié)構(gòu)破碎，不同流體開(kāi)始形成湍流混合狀態(tài),見(jiàn)圖1(c)。

為了能更清晰地展現(xiàn)爆炸過(guò)程中流場(chǎng)內(nèi)部激波與球體界面相互作用的過(guò)程，圖2給出了高壓球體爆炸過(guò)程中對(duì)稱面不同時(shí)刻的密度紋影。由圖2可看出，開(kāi)始時(shí)透射激波向外傳播、稀疏波向內(nèi)傳播以及界面失穩(wěn)過(guò)程中的詳細(xì)變化。開(kāi)始時(shí)，內(nèi)部的SF6氣體由于高壓往外膨脹，其形成的入射激波和氣體分界面作用分成透射激波和反射的稀疏波，透射激波一直往外傳播，而稀疏波則向內(nèi)傳播，見(jiàn)圖2(a)～(b)。此時(shí)為振幅線性增長(zhǎng)階段，氣體分界面在透射激波的作用下加速往外膨脹，且氣體分界面上產(chǎn)生的隨機(jī)增長(zhǎng)呈線性變化。隨后為非線性階段，重質(zhì)的SF6氣體進(jìn)入空氣中形成尖釘結(jié)構(gòu)，而輕質(zhì)的空氣演化為氣泡結(jié)構(gòu),見(jiàn)圖2(c)，且發(fā)展過(guò)程中界面RM不穩(wěn)定性呈現(xiàn)非線性增長(zhǎng)，振幅越來(lái)越大。隨著稀疏波的向內(nèi)傳播以及氣體分界面向外的膨脹，氣體分界面內(nèi)部的氣體壓力迅速降低，導(dǎo)致氣體分界面開(kāi)始向內(nèi)急劇收縮,見(jiàn)圖2(d)～(e)，從而加速界面部分流場(chǎng)的層流向湍流轉(zhuǎn)捩。至t=10 ms時(shí)，反射的稀疏波在圓心處匯聚形成反射的二次激波,見(jiàn)圖2(f)。

圖3給出了不同時(shí)刻對(duì)稱面上的徑向壓力p曲線，其中p1為透射激波前緣壓力，p2為反射稀疏波壓力。初始時(shí)刻，SF6氣體在內(nèi)部高壓下急劇膨脹，經(jīng)過(guò)氣體界面產(chǎn)生透射激波和反射稀疏波，透射激波在p1的作用下向外傳播。從t=1.4 ms到t=3.0 ms，隨著透射激波向外傳播，p1衰減較快，但是透射激波波速很高。而在t=3.0 ms之后，圖中已無(wú)p1軌跡，由此說(shuō)明透射激波此時(shí)已傳出計(jì)算域。反觀稀疏波，其壓力遠(yuǎn)低于透射激波壓力：在形成初期，稀疏波隨氣體界面向外傳播，p2逐漸減小(t=0～3.0 ms)；而隨著SF6氣體急劇膨脹，在球心形成極低壓區(qū)；隨后，稀疏波反向向內(nèi)回傳，最終在t=10 ms時(shí)在原球心位置會(huì)聚形成二次激波，此時(shí)p2達(dá)至其極大值。

圖3 徑向壓力變化Fig.3 Pressure distributions along the radius at different times

圖4 氣體界面，透射激波以及稀疏波軌跡Fig.4 Trajectories of gas interface, transmitted shock and reflected wave

圖5 二次激波與流場(chǎng)作用過(guò)程Fig.5 Processes of the interaction between reshock and the flowfield

圖4顯示了氣體分界面、透射激波，稀疏波前緣和后緣以及后期形成的二次激波的變化軌跡。從圖中可看出，透射激波以約600 m/s的速度在空氣中往外傳播，在該激波的作用下氣體分界面加速膨脹，在t=40 ms之前，氣體界面一直處于線性膨脹階段，之后，由于球心處的低壓導(dǎo)致其速度逐漸降低，最終在t=65 ms時(shí)開(kāi)始向內(nèi)收縮。同時(shí)從圖中可以看到稀疏波的傳播：稀疏波前緣在初始時(shí)向內(nèi)傳播，在球心匯聚后加速向外傳播，而其后緣則一直低速向外傳播；稀疏波的前緣約在t=2.8 ms時(shí)追上其后緣，從而匯合成同一道波，與氣體界面類(lèi)似，兩者匯聚后一直向外傳播，直到t=5 ms時(shí)開(kāi)始向內(nèi)傳播，最終在t=10 ms時(shí)在球心再次會(huì)聚碰撞，形成二次激波；該激波線性向外膨脹，直到t=124 ms時(shí)與向內(nèi)傳播的氣體界面相互作用，導(dǎo)致后期流場(chǎng)湍流發(fā)展。

圖5則刻畫(huà)了在圓心處會(huì)聚反射的二次激波往外傳播及其與流場(chǎng)的過(guò)程。此階段反射激波的強(qiáng)度很高，見(jiàn)圖3(f)，且氣體界面已不甚明顯，輕質(zhì)“氣泡”和重質(zhì)“尖釘”處的流場(chǎng)基本形成湍流。此時(shí)反射激波往外傳播，而湍流區(qū)界面則向內(nèi)傳播，因此兩者相互作用強(qiáng)烈,見(jiàn)圖5(c)，最終導(dǎo)致整個(gè)流場(chǎng)呈現(xiàn)完全湍流狀態(tài),見(jiàn)圖5(d)。

4 結(jié) 論

采用大渦模擬方法，結(jié)合高階混合格式，對(duì)SF6重質(zhì)氣云在空氣中的爆炸過(guò)程進(jìn)行了數(shù)值模擬。模擬結(jié)果清晰地再現(xiàn)了界面演化的整個(gè)過(guò)程。爆炸產(chǎn)生的激波經(jīng)過(guò)氣體分界面時(shí)分為透射激波和反射稀疏波。透射激波向外傳播，導(dǎo)致氣體分界面處RM失穩(wěn)增強(qiáng)，加速了2種氣體的混合，而反射稀疏波后緣低速向外傳播，其前緣在初始時(shí)向內(nèi)傳播，在球心匯聚后加速向外傳播，直至與稀疏波前緣匯合，兩者隨著氣體界面一直向外膨脹。之后，由于球心處的低壓導(dǎo)致回流，最終在球心處形成二次激波。在后期的發(fā)展中，該強(qiáng)激波與氣體分界面相互作用，使整個(gè)流場(chǎng)區(qū)域呈現(xiàn)完全湍流狀態(tài)。

[1] Taylor G I.The air wave surrounding an expanding sphere[J].Proceedings of the Royal Society of London: Series A: Mathematical and Physical Sciences, 1946,186:273-292.

[2] Laumbach D D.Probstein R F.A point explosion in a cold exponential atmosphere: Part 2: Radiating flow[J].Journal of Fluid Mechanics, 1970,40(4):833-858.

[3] Singh L P, Ram S D, Singh D B.Analytical solution of the blast wave problem in a non-ideal gas[J].Chinese Physics Letters, 2011,28(11):114303-114305.

[4] Marble F E, Hendrics G J, Zukoski E E.Progress toward shock enhancement of supersonic combustion processes[R].AIAA 87-1880,1987.

[5] Oran E S, Gamezo V N.Origins of the deflagration-to-detonation transition in gas-phase combustion[J].Combustion Flame, 2007,148(1/2):4-47.

[6] 孫曉暉,陳志華,張煥好.激波繞射碰撞加速誘導(dǎo)爆轟的數(shù)值研究[J].爆炸與沖擊,2011,31(4):407-412.Sun Xiao-hui, Chen Zhi-hua, Zhang Huan-hao.Numerical investigations on detonation initiation accelerated by collision of diffracted shock waves[J].Explosion and Shock Waves, 2011,31(4):407-412.

[7] Lindl J D, McCrory R L, Campbell E M.Progress toward ignition and burn propagation in inertial confinement fusion[J].Physics Today, 1992,45(9):32-40.

[8] Zabusky N.Vortex paradigm for accelerated inhomogeneous flows: Visiometrics for the Rayleigh-Taylor and Richtmyer-Meshkov environments[J].Annual Review of Fluid Mechanics, 1999,31:495-536.

[9] Brouillette M.The Richrmyer-Meshkov instability[J].Annual Review of Fluid Mechanics, 2002,34:445-468.

[10] Zoldi C A.A numerical and experimental study of a shock-accelerated heavy gas cylinder[D].New York: State University of New York, 2002: [11] Layes G, Métayer O Le.Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[J].Physics of Fluids, 2007,19:042105.

[12] 鄒立勇,劉金宏,譚多望,等.弱激波沖擊無(wú)膜重氣柱和氣簾界面的實(shí)驗(yàn)研究[J].高壓物理學(xué)報(bào),2010,24(4):241-247.Zou Li-yong, Liu Jin-hong, Tan Duo-wang, et al.Experimental study on the membraneless heavy gas cylinder and gas curtain interfaces impacted by a weak shock wave[J].Chinesee Journal of High Pressure Physics, 2010,24(4):241-247.

[13] Zheng J G, Lee T S, Winoto S H.Numerical simulation of Richtmyer-Meshkov instability driven by imploding shocks[J].Mathematics and Computers in Simulation, 2008,79(3):749-762.

[14] Pullin D I.A vortex-based model for the subgrid flux of a passive scalar[J].Physics of Fluids, 2000,12(9):2311-2319.

[15] Hill D J, Pullin D I.Hybrid tuned center-difference-WENO method for large-eddy simulation in the presence of strong shocks[J].Journal of Computational Physics, 2004,194(2):435-450.

[16] Jourdan G, Hounas L.High-amplitude single-mode perturbation evolution at the Richtmyer-Meshkov instability[J].Physical Review Letter, 2005,95(20):4502-4505.

[17] 劉金宏,鄒立勇,柏勁松,等.激波沖擊下air/SF6界面的Richtmyer-Meshkov不穩(wěn)定性[J].爆炸與沖擊,2011,31(2):135-140.Liu Jin-hong, Zou Li-yong, Bai Jing-song, et al.Richtmyer-Meshkov instability of shock-accelerated air/SF6interfaces[J].Explosion and Shock Waves, 2011,31(2):135-140.