球形重質氣體物理爆炸特性*

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

(1.南京理工大學瞬態物理重點實驗室，江蘇 南京 210094；2.浙江海洋學院海運與港航建筑工程學院, 浙江 舟山 316022；3.北京機電研究所，北京 100012)

高壓氣體在運輸、儲存、加工和使用過程中常發生爆炸事故，造成巨大的人員傷亡和財產損失，同時，由于爆炸過程還蘊含了許多諸如界面不穩定性以及激波與界面相互作用等復雜流體物理現象，因此對氣云物理爆炸特性的研究不僅具有重要的實際意義，還具有重要的理論價值。

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

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

1 控制方程和數值方法

三維可壓Navier-Stokes(N-S)方程經過Favre濾波后，其標準的可壓三維LES方程作為控制方程。連續性方程、動量方程和能量方程分別為：

(1)

(2)

(3)

在大渦模擬中，亞格子湍流應力項無法直接從可解尺度中計算求得,因此亞格子應力需要引入模型進行封閉。這里采用拉伸渦亞格子模型[14]對小尺度輸運項建模。拉伸渦亞格子模型由D.I.Pullin[14]提出，他假設亞網格內的流動是由亞網格內的渦導致的，并通過對稱渦結構來構建小尺度渦量。拉伸渦亞格子模型已被成功用于亞格子能量譜的估算以及湍流混合模擬。

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

2 計算模型

在爆炸初期，流場結構呈現對稱情形，而由于RM失穩，后期的流場漸漸轉為非對稱結構。為描述RM失穩以及后期二次激波的作用過程，選擇初始狀態為球形氣云的SF6氣體作為研究對象。SF6氣體和空氣的密度分別為5.97和1.18 kg/m3，兩者的比熱比和氣體常量分別為γ(SF6)=1.09,R(SF6)=56.92;γ(air)=1.40,R(air)=287。選擇標準狀態的空氣參量，p0=101.3 kPa，T0= 287 K。SF6氣體的初始壓力為50p0,SF6氣云初始半徑為R0=0.1 m，爆炸空間徑向大小為15R0。計算采用笛卡爾網格模擬氣云球體表面，以該網格擾動作為初始隨機擾動，球體表面網格總量約為800萬，整個計算域網格總量約為4 500萬。

3 結果與討論

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

圖2 流場波系結構及界面隨時間的變化Fig.2 Evolution of the shock wave structure and interface in the flowfield

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

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

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

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

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

圖5 二次激波與流場作用過程Fig.5 Processes of the interaction between reshock and the flowfield

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

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

4 結 論

采用大渦模擬方法，結合高階混合格式，對SF6重質氣云在空氣中的爆炸過程進行了數值模擬。模擬結果清晰地再現了界面演化的整個過程。爆炸產生的激波經過氣體分界面時分為透射激波和反射稀疏波。透射激波向外傳播，導致氣體分界面處RM失穩增強，加速了2種氣體的混合，而反射稀疏波后緣低速向外傳播，其前緣在初始時向內傳播，在球心匯聚后加速向外傳播，直至與稀疏波前緣匯合，兩者隨著氣體界面一直向外膨脹。之后，由于球心處的低壓導致回流，最終在球心處形成二次激波。在后期的發展中，該強激波與氣體分界面相互作用，使整個流場區域呈現完全湍流狀態。

[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] 孫曉暉,陳志華,張煥好.激波繞射碰撞加速誘導爆轟的數值研究[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] 鄒立勇,劉金宏,譚多望,等.弱激波沖擊無膜重氣柱和氣簾界面的實驗研究[J].高壓物理學報,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不穩定性[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.