T型管内流动气体中爆轰绕射过程的数值模拟*

潘振华，范宝春，归明月

(1.江苏大学能源与动力工程学院，江苏 镇江 230027；2.南京理工大学瞬态物理国家重点实验室，江苏 南京 210094)

爆轰波传入突扩通道时，如翻越障碍物或进入分叉管等，将发生绕射。爆轰绕射现象是爆轰波传播过程中经常出现的现象，与许多工业过程有关。爆轰绕射的演化过程涉及许多基本问题，如激波阵面和化学反应阵面的解耦、熄火、激波反射与相互作用、局部爆炸和二次起爆等，因此爆轰绕射现象引起广泛关注。爆轰绕射可以发生在静止介质中，也可以发生在流动介质中，两者的演化过程不尽相同。但迄今，已经进行的爆轰绕射的研究，皆限于静止系统。

I.B.Zeldovich等[1]通过实验研究了爆轰绕射，提出爆轰管临界管径的概念，即其他参数一定的条件下，存在临界管径，当爆轰管直径小于该值时，爆轰进入突扩空腔后会熄灭。此后，V.V.Mitrofanov等[2]发现，临界管径与爆轰胞格宽度有关。R.I.Soloukin等[3]通过实验得到了爆轰绕射时激波与波后化学反应区分离以及熄灭后再生的纹影照片。S.B.Murray等[4]提出了爆轰波绕射后，熄灭再生的2种机制：(1) 反射诱导起爆机制；(2) 局部爆炸引起的自起爆机制。A.Smolinska等[5]通过实验和计算，得到爆轰绕射后的清晰胞格结构。F.Pintgen等[6]通过实验观察得到了爆轰波在2种不同混合气中的绕射过程，分别研究了爆轰波绕射后的3种传播状态。C.M.Guo等[7]和C.J.Wang等[8]对爆轰波在分叉管中的传播过程进行了实验和数值研究，其中涉及到了爆轰波衰减、马赫反射、规则反射向马赫反射的转变以及爆轰波的重新起爆过程。

与静止介质不同，流动介质中传播的爆轰可分为迎风和顺风2类。爆轰波绕射进入水平管后，激波面与高速来流作用，一侧为迎风传播，另一侧为顺风传播。M.Y.Gui等[9]基于Euler方程数值模拟了尖劈诱导的斜爆轰的胞格结构，得到了迎风传播的爆轰阵面随时间和空间位置的变化特征。T.H.Yi等[10]通过数值模拟，讨论了一维流动系统中的爆轰波，得出了迎风面的传播速度小于顺风面。潘振华等[11]通过数值模拟，讨论了二维流动系统中的爆轰波的传播特性，得到了在流动介质中，爆轰流场不对称，爆轰压力强度在上游、横向和下游方向上依次递减。K.Ishii等[12]通过实验得到了爆轰波分别在600和910 m/s的高速预混气流中传播后留下的胞格轨迹，结果表明：向上游传播的迎风爆轰波，爆速高于CJ值，胞格被压缩；向下游传播的顺风爆轰波，爆速降低，胞格被拉长。

流动介质中的爆轰传播问题，迄今尚无系统研究，特别是涉及到爆轰波的绕射、再起爆以及稳定传播的问题。本文中以H2、O2、Ar体积比为2∶1∶1的预混气体为研究对象，基于带化学反应的二维Euler方程和五阶精度的WENO(weighted essentially non-oscillatory)格式，对K.Ishii等的实验中[12]爆轰波由T型管的垂直管进入通有超高速气流的水平管中时产生的爆轰绕射及其后续演化过程进行数值研究，与相应静止介质中的爆轰绕射进行比较，分析流动系统中爆轰绕射的流场变化特征和三波点烟箔轨迹的变化特征，探索爆轰再生的动力学机理。

1 数值研究方法

1.1 控制方程

假设混合气体为理想气体，忽略扩散、粘性和热传导。在贴体坐标系中，带基元化学反应的多组分二维Euler方程为：

(1)

(2)

(3)

(4)

(5)

可燃预混气体为体积比为2∶1∶1的H2、O2、Ar混合气体，基元反应模型采用9组元和48个化学反应的详细化学反应机理[13]，反应组元分别为H、O、H2、OH、H2O、O2、HO2、H2O2、Ar。

1.2 计算方法

计算过程中，为处理爆轰化学反应带来的刚性问题，采用附加半隐的龙格-库塔法(additive semi-implicit Runge-Kutta methods，ASIRK)[14]，该方法在时间上具有二阶精度，且有很好的稳定性。空间项采用五阶精度的WENO格式[15]计算。化学反应源项采用基于Gear格式的LSODE程序[16]计算。

1.3 计算模型

图1分别为静止系统和流动系统的计算域示意图。图1(a)为静止系统，长31 cm、宽3.2 cm的水平燃烧室中心位置与长1.8 cm、宽3 cm的垂直爆轰管相连。图1(b)为流动系统，燃烧室长31 cm，宽3.2 cm，爆轰管长1.8 cm、宽3 cm，两者相连，爆轰管中心轴线在水平管9.6 cm位置上。整个计算区域采用均匀化网格，网格尺寸为0.1 mm×0.1 mm。其中燃烧室使用的网格数为3 100×320，爆轰管使用的网格数为180×300。

图1 计算域示意图Fig.1 The computational domain

对初压为10.6 kPa、初温为300 K的预混气体进行前期数值计算，在垂直爆轰管中形成稳定自持传播的爆轰波，管道内横波数为7，对应胞格数稳定在3.5，以此初压下自持胞格爆轰作为初始化计算。在其他区域中，仍充满以等当量比的H2、O2为燃料、Ar作为稀释剂的预混气体，三者体积比为2:1:1。在静止系统中，水平管内的预混气体流速为零。在流动系统中，左侧为进气端，右侧为出口端，预混气以1 km/s的速度在燃烧室内从左向右流动。设静止系统和流动系统2个算例初始压力为10.6 kPa，温度为300 K。计算中，均采用量纲一量，量纲一参考值为pf=15.4 kPa，Tf=300 K，L0=0.1 m。边界处理中，CDJ、KEF和BA用无催化、绝热固壁边界条件；边界AC和BF用外推边界条件；JK用插值外流边界条件。

2 数值验证

为了验证计算的准确性，对E.S.Oran等[17]的算例进行数值模拟，与E.S.Oran等[17]的计算结果进行对比，并将计算得到的爆轰波速与C.A.Eckett[18]的实验结果进行比较。

图2 直管中的计算胞格Fig.2 Numerical cells structure in straight tube

为了直接与文献[17]计算结果进行对比，数值验证采用与文献[17]相同的计算域和初始条件，即在直管中，充满H2、O2、Ar的预混气体，三者体积比为2∶1∶7，初始压力和温度分别为6.67 kPa和298 K。图2为本文数值模拟得到的胞格图像，表1为本文计算结果与文献[17]的对比，其中：Δx、Δy分别为单个网格横向和纵向的尺寸，va为爆轰波的平均速度，L和W分别为胞格横向和纵向的尺寸。计算得到的胞格尺寸为约54 mm×30 mm，与E.S.Oran在相同条件下得到的胞格尺寸54 mm×31 mm[17]基本一致；胞格纵横比为0.556，与文献[18]中实验得到的纵横比0.54符合较好。本文中爆轰波的传播速度为1 589 m/s，与Gordon-McBride程序算出的理论爆轰波CJ爆速(1 617 m/s)相比误差在小于2%，与文献[18]实验得到的爆轰波速度(1 550 m/s)基本一致。因此，本文中的数值模拟能够重复文献[17]的计算，并且计算结果与文献[18]的实验结果相符。

表1 数值验证结果与文献[17]计算结果的对比Table 1 Comparison between current simulation and referance [17]

3 结果与讨论

3.1 爆轰绕射的波阵面结构

图3 t=16.2 μs时静止系统爆轰波纹影图Fig.3 Numerical schlieren of detonation in a quiescent system while t=16.2 μs

图3为t=16.2 μs时，静止系统中爆轰绕射流场的计算纹影图。其中右图为左图的局部放大显示。此时，爆轰阵面尚未抵达水平管的上壁面。右侧局部放大图中，实线为激波阵面，虚线为化学反应阵面。由图可知，在未受扰动区域，激波与反应阵面耦合很好，平面形状的爆轰波以固有的速度向前发展。但在被扰动区域，阵面弯曲，激波已与反应阵面解耦，两者之间存在一个经前导激波预压的高温高压未燃气体区域。

图4为t=16.2 μs时，流动系统中爆轰绕射流场的计算纹影图。其中左图为右图的局部放大显示。对于右侧的顺风绕射，被扰动的波阵面上，激波与反应阵面解耦，大致情形与静止系统的相仿。但对于左侧的迎风绕射，从垂直管中传出的爆轰，其波后爆轰产物对水平管内的定向来流而言，具有与尖劈类似的压缩作用，这使得迎风绕射的爆轰阵面具有斜爆轰的某些特点。稳定的斜爆轰由一系列子波构成，其精细结构如图5所示[9-11]。由图5可知，激波阵面的弯曲，使阵面分为S1和S2两部分，并在三波点P处碰撞，形成横波T1，并诱发化学反应，也进而使T1和S2在三波点附近的部分成为爆轰波D1，T1称为横向爆轰波。与正爆轰的精细结构相仿，此类子波结构(精细结构)使得斜爆轰得以稳定自持。图4左侧的局部放大图中同样出现了稳定斜爆轰所具有的子波结构(如三波点P1和P2)。

图4 t=16.2 μs时流动系统中的计算纹影图Fig.4 Numerical schlieren of detonation in a flowing system while t=16.2 μs

图5 迎风面上波系结构示意图Fig.5 Schematic diagram of wave structure on windward side

3.2 爆轰绕射的再生机理

当爆轰波抵达水平管的上壁面时，将在壁面反射。图6给出了流通系统中爆轰波在上壁面反射后不同时刻流场的计算纹影、压力分布、温度分布和H组分的质量分数分布图。

图6 流动系统中爆轰波的再生及波系结构的演变Fig.6 Reinitiation event of detonation and evolution of wave structure in a flowing system

激波在壁面的反射包括2部分：未受扰动的平面爆轰在壁面的正反射R1和随后的被扰动的弯曲爆轰在壁面的马赫反射R2。反射后，在壁面附近形成向两端传播的马赫干M和向下方传播的横向反射波R(图6(a)和图6(b)纹影图中的水平白线)。反射波波后为高温高压区域(见图6的压力和温度图)。由于马赫干的强度较高，故可诱导爆轰波。而横向传播的反射波，在传播过程中，有部分阵面扫过因激波阵面与反应区解耦而形成的预压未燃区，从而使之点火，形成横向传播的爆轰波T。因此，整个反射波阵面是爆轰波和激波组成的复合波阵面，两端为弯曲的爆轰波T，中间由激波R连接。在反射波阵面向下传播的同时，左侧类斜爆轰阵面上的三波点P3和P4也作相对运动，直至碰撞(如图6(a)～(c)所示)，这与斜爆轰精细结构的情形一样。左侧的反射波(爆轰子波)在与相对稳定的类斜爆轰阵面的相互作用过程中向左传播并向下扩展，阵面趋于平整并与水平方向垂直。而右侧反射波(强度较弱的爆轰子波)在与不断衰减的激波的相互作用过程中向右传播并向下扩展，其强度也不断衰减，从而形成倾斜的波阵面。在形成斜爆轰的过程中，激波阵面与火焰阵面未很好耦合(如图6(a)～(c)所示)。由图6～9还可看出，虽然迎风爆轰的强度大于顺风爆轰，但顺风爆轰的传波速度大于迎风爆轰。此后，顺风斜爆轰在下壁面反射，由于反射爆轰的强度大于入射爆轰，故随着反射爆轰阵面的扩展，爆轰得以加强(如图6(d)所示)。经过上下壁面的若干次反射，右侧爆轰波也可逐渐发展为与水平方向垂直的正爆轰。至此，爆轰波经历了绕射导致的熄火和二次起爆，最终成为沿水平管，迎风和顺风向两端传播的正爆轰波。静止系统中爆轰波绕射过程与流动系统顺风侧的情形类似。

3.3 爆轰绕射过程的胞格结构

从精细流场角度讲，稳定传播的爆轰阵面的强度是不均匀的。该阵面是由一系列马赫碰撞事件构成的，三波点上的压力最高。三波点(即高压点)的轨迹通常具有图2所示的网状图形，称为爆轰胞格。图7为静止系统中爆轰在T型管中绕射时胞格结构演变图。爆轰在垂直管中稳定传播时，形成图7中a1区的规则胞格。爆轰进入水平管绕射时，部分阵面因扰动而熄火，b1区域为爆轰波未受扰动区域，胞格形状不变，与a1区一致。c1区为“爆轰熄火”区，胞格消失。b1区与c1区的边界为倾斜直线，倾斜角等于扰动角。绕射爆轰抵达上壁面后，在壁面反射，反射波传播区域中，仅有强度最大的横向爆轰留下传播轨迹，即图中d1区域。反射波在下壁面再次反射，形成向两端传播，强度不大的斜爆轰波。波阵面上，马赫碰撞的三波点单向运动，形成e1区平行的三波点轨迹线。反射波在上下壁面反复反射，强度不断增加，最终会重新成长为稳定传播的平面爆轰，胞格也重新形成f1区中规则的网状结构。

图8为流动系统中爆轰绕射的胞格结构演变图。受水平管中定向流动的影响，未受扰动区b2向右侧偏移，胞格也发生倾斜。在右侧(顺风侧)，胞格图像与图7中大致相同，但各区的面积被拉长。此外，横向爆轰的轨迹出现两次(图8中e2区域)，而静止系统只有一次，这说明静止系统中，绕射爆轰的增长速度更快。在左侧(迎风侧)，爆轰绕射时，基本未出现熄灭阵面(熄火仅出现在拐角附近很小的c2区域内)，绕射阵面上三波点以及发源于此的横向爆轰的运行轨迹如图中d2区所示。反射波抵达下壁面时，已基本成为正爆轰，左侧g2区胞格尺寸被显著压缩，右侧g2区胞格尺寸被显著拉长。

图8 流动系统中的胞格结构演变Fig.8 Cellular structure in a flowing system

4 结 论

对比了静止系统和流动系统中爆轰波从绕射到再生的整个过程，得到如下结论。(1) 爆轰波在流动系统中发生绕射，顺风传播的爆轰波趋于熄火；在迎风侧，三波点附近局部区域内激波阵面始终与波后化学反应区耦合在一起，形成斜爆轰结构，使爆轰得以维持。(2) 由于两侧爆轰传播速度不同，迎风侧的胞格图像明显被压缩，顺风侧被拉长。(3) 在爆轰波的传播过程中，波阵面上横向爆轰波的产生对爆轰波的起爆过程起到了关键作用。

[1] Zeldovich I B, Kogarko S M, Simonov N N.An experimental investigation of spherical detonation in gases[J].Soviet Physics: Technical Physics, 1956,1(8):1689-1713.

[2] Mitrofanov V V, Soloukhin R I.The diff'raction of mutifront detonation waves[J].Soviet Physics: Doklady, 1965,9(12):1055-1058.

[3] Soloukhin R I, Ragland K W.Ignition processes in expanding detonations[J].Combustion and Flame, 1965,13(3):295-302.

[4] Murray S B, Lee J H.On the transformation of planar detonation to cylindrical detonation[J].Combustion and Flame, 1983,52:269-289.

[5] Smolinska A, Khasainov B, Virot V, ea al.Detonation diffraction from tube to space via frontal obstacle[C]∥Proceedings of the Fourth European Combustion Meeting.Vienna, Austria, 2009.

[6] Pintgen F, Shepherd J E.Detonation diffraction in gases[J].Combustion and Flame, 2009,156(3):665-677.

[7] Guo C M, Wang C J, Xu S L, ea al.Cellular pattern evolution in gaseous detonation diffraction in a 90°-branched channel[J].Combustion and Flame, 2007,148(3):89-99.

[8] Wang C J, Xu S L, Guo C M.Gaseous detonation propagation in a bifurcated tube[J].Journal of Fluid Mechanics, 2008,599:81-110.

[9] Gui M Y, Fan B C.Wavelet structure of wedge-induced oblique detonation waves[J].Combustion Science and Technology, 2012,184(10/11):1456-1470.

[10] Yi T H, Wilson D R, Lu F K.Numerical study of unsteady detonation wave propagation in a supersonic combustion chamber[R].Arlington, Texas, USA: University of Texas at Arlington, 2004.

[11] 潘振华，范宝春，归明月，等.流动系统中爆轰波传播特性的数值模拟[J].爆炸与冲击,2010，30(6)：593-597.Pan Zhen-hua, Fan Bao-chun, Gui Ming-yue, et al.Numerical study of detonation wave propagation in a flow system[J].Explosion and Shock Waves, 2010,30(6):593-597.

[12] Ishii K, Kataoka H, Kojima T.Initiation and propagation of detonation waves in combustible high speed flows[J].Proceedings of Combustion Institute, 2009,32(2):2323-2330.

[13] Oran E S, Young T R, Boris J P, et al.Weak and strong ignition I: Numerical simulation of shock tube experiments[J].Combustion and Flame, 1982,48:135-148.

[14] Zhong X L.Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows[J].Journal of Computational Physics, 1996,128(1):19-31.

[15] Jiang G S, Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1996,126(2):202-228.

[16] Radhakrishnan K, Hindmarsh A C.Description and use of LSODE, the livermore solver for ordinary differential equations, UCRL-ID-113855[R].Washington, USA: NASA, 1993.

[17] Oran E S, Weber J W, Stefaniw E I, et al.Numerical study of a two-dimensional H2-O2-Ar detonation using a detailed chemical reaction model[J].Combustion and Flame, 1998,113(1/2):147-163.

[18] Eckett C A.Numerical and analytical studies of the dynamic of gaseous detonation[D].CA, USA: California Institute of Technology, 2001.