球形重质气体物理爆炸特性*

薛大文，陈志华，韩珺礼,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.