激波诱导火焰失稳与爆轰的条件研究＊

朱跃进，于 蕾，张彭岗，潘振华，潘剑锋，董 刚

（1．江苏大学能源与动力工程学院，江苏镇江212013；2．南京理工大学瞬态物理重点实验室，江苏南京210094）

激波诱导火焰失稳与爆轰的条件研究＊

朱跃进1，2，于 蕾1，张彭岗1，潘振华1，潘剑锋1，董 刚2

（1．江苏大学能源与动力工程学院，江苏镇江212013；2．南京理工大学瞬态物理重点实验室，江苏南京210094）

采用九阶WENO和十阶中心差分格式数值求解激波与火焰作用过程，考察了激波强度、火焰尺寸对激波与球形火焰作用过程的影响。结果表明，增大激波强度或火焰尺寸均可在流场中引发爆轰，但激波强度的影响更大，并且其引发的爆轰可使火焰迅速膨胀，放热率提高，从而影响燃烧特性；此外，爆轰波传播过程中会迅速消耗可燃预混气，合并原有的反射激波，并在流场中形成局部高压区，极大地改变流场结构。

激波；火焰；爆轰；流场结构

激波与火焰作用时，火焰界面会发生失稳变形，流场中形成的涡量能够引起未燃气与已燃火焰的混合，从而促进燃烧放热；当初始条件改变时，甚至可在流场中引发爆轰。该过程常见于超燃推进［1］和工业爆炸灾害［2］等领域，因此，开展激波与火焰相互作用过程的研究具有重要的理论和实际意义。

早在20世纪五六十年代，G．H．Markstein［3］已经在激波管内开展了激波与火焰相互作用的实验研究，得到了火焰在入射及反射激波作用下的变形过程。21世纪初，G．O．Thomas等［4］针对激波与球形火焰的相互作用过程开展了实验研究，通过高速摄影发现：较强的入射激波及平面反射激波可以通过Richtmyer－Meshkov（RM）不稳定性使火焰形状发生严重变形；当入射激波增强时，激波作用后的变形火焰面附近会形成热点并最终发展为爆轰。另一方面，随着计算机发展，数值模拟工作开始得到广泛开展。Y．Ju等［5］研究了入射激波强度对火焰变形的影响，发现火焰总燃烧速率随激波强度的增强而增大。E．S．Oran等［2，68］详细研究了激波与火焰作用，发现当条件适当时，反射激波与火焰的相互作用可以产生热点和爆轰。H．H．Teng等［9］对激波与爆燃火焰的相互作用进行了二维数值模拟，发现入射激波扫过火焰面后会导致其发生RM不稳定，进而可能引发爆轰。谷壮志等［10］采用CE／SE方法研究了激波诱导火焰的变形过程，发现当激波增强时，流场中可能引发爆轰，但是计算精度尚显不够。董刚等［11］和朱跃进等［1215］采用三阶PPM算法对激波与火焰作用问题开展了数值研究，发现流场中可能出现爆轰，但未对爆轰波的进一步发展进行分析。

上述实验与计算结果均表明，激波与火焰的作用过程易受初始条件影响，流场内可在某些特定位置引发爆轰，但目前尚缺乏精细的波系与火焰界面数据。因此，本文中，采用高精度的九阶WENO［1617］和十阶中心差分格式对激波与火焰作用过程进行分析，通过改变激波强度和初始火焰尺寸获得爆轰引发与传播的数据，以期为揭示激波诱导爆轰现象的内在机理提供基础。

1 数理模型与实验验证

1．1 数理模型

二维带化学反应的Navier－Stokes（N－S）方程为：

数值计算过程中，对控制方程式（1）中空间导数项的黏性部分采用十阶空间中心差分计算，无黏部分采用高精度九阶WENO格式［1617］进行计算；时间推进采用三阶Runge－Kutta方法求解。另外，计算过程满足CFL数值稳定条件。

1．2 实验验证

根据文献［4］中的实验设置计算区域，如图1所示：流向（x）长为0．170m，法向（y）高为0．076m。计算域内充满C2H4／3O2／4N2组成的反应性预混气体，密度ρ0＝0．161 5kg／m3，初始温度T0＝293K；初始球形火焰的半径R0＝0．019m，火焰的中心位于x＝0．038m，y＝0．038m处，火焰区域内部密度ρ1＝0．015 78kg／m3，火焰内、外初始压力p0＝13．3kPa；初始入射激波位于x＝0．012m处，沿x方向从左向右传播，入射激波马赫数为1．7，激波后气体状态由Rankine－Hugoniot关系式给出。

设计算区域的上、下边界和x方向右端面为无滑移的刚性绝热壁面边界，而x方向左端面入口处采用零梯度边界，以保证激波强度不变。计算使用的均匀网格尺寸为0．1mm，由于预混气体C2H4＋3O2＋4N2在本文所设定条件下的火焰面厚度约为0．785mm［11］，因此，火焰面厚度由7～8个网格来刻画，网格分辨率可以满足反映火焰界面变化过程的要求。

为验证上述模型及数值计算方法的可靠性，将文献［4］中的实验结果与本文所采用数学模型和计算方法进行了对比。图2给出了不同时刻的计算密度图与实验纹影的比较，其中图2（a）～（b）为入射激波与火焰作用后的情形，图2（c）～（d）为反射激波与失稳火焰再次作用后的情形。从图2可以看出，各时刻计算得到的激波阵面形状、失稳火焰结构均与实验结果定性一致，在一定程度上说明了本文数理模型和计算方法的可靠性和合理性。需要注意的是，实验结果中，火焰存在贴壁现象［4，14］，而数值计算方法无法反映该现象，这是由于二维计算模型无法刻画真实火焰的三维结构，采用二维计算获得的火焰尺寸相比于实验结果要小。

图1 计算区域示意图Fig．1 Sketch of computational domain

图2 实验纹影［4］与计算密度的比较Fig．2 Comparison between experimental schlieren images［4］and computational density images at selected times

2 结果与讨论

激波与火焰的作用过程易受初始条件影响，当初始入射激波强度或者火焰尺寸增大时，初始球形火焰的失稳程度将大幅增加，燃烧区的放热量大大升高，流场内有可能出现爆轰现象，从而极大地改变流场结构。因此，为考察激波强度和火焰尺寸两个因素对流场变化的影响，本文中共设计6组算例，参数取值如表1所示。其中Case 1为基本算例，依据设定的实验条件，Case 2、Case 3与之相比，仅不断增大初始入射激波强度，而Case 4～Case 6中的初始球形火焰尺寸统一增大为0．024m，入射激波强度分别与Case 1～Case 3对应。

表1 不同激波马赫数和火焰尺寸的4组算例Table 1 Four cases with different shock Mach numbers and flame sizes

2．1 入射激波Ma的影响

当初始入射激波Ma＝2．1时（Case 2），入射激波波后与反射激波波后典型时刻的计算密度如图3所示。由于入射激波强度增大，初始火焰受到更强的压缩作用，与Case 1相比，此时入射波和反射波后的温度和压力更高，有利于火焰区的燃烧膨胀，尤其对于反射波后情形，因此图3（b）中的失稳火焰区域大幅增大。尽管如此，Case 2在本文计算条件下未在流场内形成爆轰。

图3 计算密度（Case 2）Fig．3 Computational density images（Case 2）

当入射激波强度进一步增大（Case 3），激波与火焰的作用过程更加复杂。图4为Case 3条件下反射激波与失稳火焰作用并引发爆轰的情形，其中每幅子图的上半部分为计算密度，下半部分为预混气组分（反应物质量分数，用“S”表示）。图4（a）对应于反射激波刚扫过失稳火焰的时刻，由于火焰内部声阻抗较小，激波阵面出现一定程度的弯曲，并在流场中心对称面附近碰撞形成高压和高密度区，有利于反应放热过程。由于高压区内有大量的未燃气，因此随着时间推进，该高压区内出现的局部热点可发展形成向周围传播的爆轰波（由数值计算结果得到，该时刻波阵面的传播速度约为1 990m／s，波后压力约为1．1MPa；而根据波前气体状态，由Gordon－McBrid程序［18］计算可得，理论C－J爆速DCJ约为2 041m／s，爆压pCJ约为1．2MPa，可以看出数值计算结果与理论吻合较好），高压区内及其附近大量未燃气被迅速消耗，如图4（b）～（c）所示。当t＝237μs时，爆轰波与反射激波阵面完全合并，并继续向左传播，其中向上、下壁面传播的爆轰波与壁面作用并发生反射，形成新的反射激波（reflected shock wave，RSW），如图4（d）所示。该RSW从壁面向流场中心传播，一方面与向左传播的爆轰波碰撞形成马赫杆，其波后压力比爆轰波阵面后压力高；另一方面逐渐分化为向左下方传播的RSW1和向流场中心传播的RSW2（见图4（e））。流场区域内上下对称的RSW1和RSW2很快在流场中心对称面处发生部分碰撞并再次反射，形成向左上方传播的RSW3，如图4（f）所示，同时波后形成新的高压区。不难看出，增强入射激波可以显著改变流场结构。随着时间的进一步发展，向左传播的爆轰波阵面将逐渐变成平面爆轰波，且爆轰波阵面上存在许多三波点与横波，波后可燃组分为零，体现了本文高计算精度的优势。

图4 计算密度与质量分数（Case 3）Fig．4 Images of computational density and mass fraction（Case 3）

2．2 初始火焰尺寸的影响

为考察初始火焰尺寸对流场变化的影响，本节中仅选取表1中的典型算例Case 5进行分析，其入射激波强度与Case 2一致。图5给出了Case 5条件下反射激波扫过失稳火焰后的情形，其中每幅子图的上半部分为计算密度图，下半部分为预混气体组分。与图3（b）相比，图5（a）中反射激波在相对短的时间内运动的距离更长，表明此时反射激波运动速度更高，且波后的失稳火焰面积更大，这是因为初始球形火焰尺寸增大，更有利于激波的透射传播与火焰发展。同时随着反射激波继续向左运动，图5（b）在反射波后出现了3处局部爆轰，压力均高于1MPa，温度也均高于3 000K，而且从对应的组分图中可以发现，这3处爆轰应形成在反射波后的未燃区，可能是由反射激波与失稳火焰面之间的局部热点形成的。值得注意的是，上壁面的爆轰波附近存在大量的可燃气，因此，该处的爆轰能够向周围持续传播，如图5（c）所示。当t＝524μs时，该爆轰波部分与左传的反射激波合并，部分沿流向向右传播，此外还有部分在流场区域的中心对称面处发生碰撞，形成两处局部高压区。因此，增大初始火焰尺寸也可以显著改变流场结构，整个反射波后流场只在爆轰波尚未传播到的地方存在少量未燃气。图5中爆轰波阵面的三波点与横波结构十分清楚，也体现了高精度计算格式的优势。

图5 计算密度与质量分数（Case 5）Fig．5 Images of computational density and mass fraction（Case 5）

2．3 积分性质分析

为进一步分析各初始参数对激波诱导火焰失稳与爆轰现象的影响，采用积分形式定义火焰有效面积（A）和平均化学反应放热率（λ）如下：

式中：下标“D”和“F”分别表示整个计算区域和火焰区（反应物质量分数Y≤0．99），X为反应物体积分数，A为整个计算域内去除混合影响后的火焰面积。图6和图7给出了量纲一火焰有效面积（A／A0）和平均化学反应放热率随时间的变化过程，其中A0为Case 1中初始火焰面积，垂直的点划线表示各算例中反射激波与变形火焰作用前的时刻。由图6和图7可知，Case 3和Case 6中的变形火焰与反射激波作用后急剧膨胀，且火焰区平均放热率均较高，与图4中反射激波波后火焰区形成爆轰一致。此时两者的平均放热率较接近，表明该马赫数下初始火焰尺寸的影响很小；而Case 2和Case 5中反射激波波后的变形火焰面积增长趋势明显放缓，同时火焰区的平均放热率均小于Case 3和Case 6，体现了入射激波强度减弱后对流场的影响，由于Case 5中形成了爆轰（见图5），因此Case 2与Case 5的火焰区面积和平均放热率在500μs之后开始显现差异。综上所述，改变初始条件可使火焰面积和平均放热率产生较大差异，特别是入射激波强度的影响更显著，而反射激波后的放热率相比于入射激波作用后的放热率呈现数量级的增加，说明反射激波能极大地促进燃烧过程，提高放热率，从而影响火焰区燃烧特性，改变流场结构。

图6 火焰有效面积随时间的变化Fig．6 Time histories of flame effective area

图7 火焰平均反应放热率随时间的变化Fig．7 Time histories of average reaction heat release rate

3 结 论

采用高精度计算格式对平面入射激波诱导火焰失稳与爆轰的现象进行了详细的数值研究，获得了较精细的波系和火焰结构，并考察了激波强度、初始火焰尺寸对激波与火焰作用过程的影响，结果表明：（1）激波强度越高，火焰受到的压缩作用越强，反射激波后的火焰面积越小，但反射波后的高温、高压条件有利于尽早形成爆轰；（2）若固定激波强度，仅增大火焰尺寸，则反射激波后的火焰面积明显增大，但引发爆轰需要较长时间；（3）流场内引发的爆轰可使火焰迅速膨胀，放热率大大提高，从而影响燃烧特性，并且爆轰波向四周传播时，会迅速消耗周围的可燃预混气，一部分与反射激波面发生合并，另一部分在中心对称面附近发生碰撞形成局部高压区，从而极大地改变流场结构。

［1］ Marble F E，Hendricks G J，Zukoski E E．Progress toward shock enhancement of supersonic combustion processes［C］∥23rd Joint Propulsion Conference．San Diego，CA，1987．

［2］ Oran E S，Gamezo V N．Origins of the deflagration－to－detonation transition in gas－phase combustion［J］．Combustion and Flame，2007，148（1／2）：4－47．

［3］ Markstein G H．A shock－tube study of flame front－pressure wave interaction［J］．Symposium（International）on Combustion，1957，6（1）：387－398．

［4］ Thomas G O，Bambrey R，Brown C．Experimental observations of flame acceleration and transition to detonation following shock－flame interaction［J］．Combustion Theory and Modelling，2001，5（4）：573－594．

［5］ Ju Y，Shimano A，Inoue O．Vorticity generation and flame distortion induced by shock flame interaction［J］．Symposium（International）on Combustion，1998，27（1）：735－741．

［6］ Khokhlov A M，Oran E S，Thomas G O．Numerical simulation of deflagration－to－detonation transition：The role of shock－flame interactions in turbulent flames［J］．Combustion and Flame，1999，117（1／2）：323－339．

［7］ Khokhlov A M，Oran E S．Numerical simulation of detonation initiation in a flame brush：The role of hot spots［J］．Combustion and Flame，1999，119（4）：400－416．

［8］ Gamezo V N，Oran E S，Khokhlov A M．Three－dimensional reactive shock bifurcations［J］．Proceedings of the Combustion Institute，2005，30（2）：1841－1847．

［9］ Teng H H，Jiang Z L，Hu Z M．Detonation initiation developing from the Richtmyer－Meshkov instability［J］．Acta Mechanica Sinica，2007，23（4）：343－349．

［10］ 谷壮志，王超，施红辉，等．激波诱导火焰变形的数值模拟［J］．浙江理工大学学报，2011，28（4）：529－533．Gu Zhuangzhi，Wang Chao，Shi Honghui，et al．The numerical simulation of the flame deformation induced by shock wave［J］．Journal of Zhejiang Sci－Tech University，2011，28（4）：529－533．

［11］ Dong G，Fan B，Ye J．Numerical investigation of ethylene flame bubble instability induced by shock waves［J］．Shock Waves，2008，17（6）：409－419．

［12］ 朱跃进，董刚，范宝春．受限空间内激波与火焰作用的三维计算［J］．推进技术，2012，33（3）：405－411．Zhu Yuejin，Dong Gang，Fan Baochun．Three－dimensional computation of the interactions between shock waves and flame in a confined space［J］．Journal of Propulsion Technology，2012，33（3）：405－411．

［13］ 朱跃进，董刚，刘怡昕，等．激波诱导火焰变形、混合和燃烧的数值研究［J］．爆炸与冲击，2013，33（4）：430－437．Zhu Yuejin，Dong Gang，Liu Yixin，et al．A numerical study on shock induced distortion，mixing and combustion of flame［J］．Explosion and Shock Waves，2013，33（4）：430－437．

［14］ Zhu Y J，Dong G，Liu Y X．Three－dimensional numerical simulations of spherical flame evolutions in shock and reshock accelerated flows［J］．Combustion Science and Technology，2013，185（10）：1415－1440．

［15］ 朱跃进，董刚．激波冲击火焰的涡量特性研究［J］．爆炸与冲击，2015，35（6）：839－845．Zhu Yuejin，Dong Gang．A study of vorticity characteristics of shock－flame interaction［J］．Explosion and Shock Waves，2015，35（6）：839－845．

［16］ Jiang G S，Shu C W．Efficient implementation of weighted ENO schemes［J］．Journal of Computational Physics，1996，126：202－228．

［17］ 王革，关奔．激波作用下R22气泡射流现象研究［J］．力学学报，2013，45（5）：707－715．Wang Ge，Guan Ben．A study on jet phenomenon of R22gas cylinder under the impact of shock wave［J］．Chinese Journal of Theoretical and Applied Mechanics，2013，45（5）：707－715．

［18］ Gordon S，Mcbride B J．Computer program for calculation of complex chemical equilibrium compositions and applicationⅠ：analysis［Z］．NASA Reference Publication 1311，1994．

Conditions for shock wave induced flame instability and detonation

Zhu Yuejin1，2，Yu Lei1，Zhang Penggang1，Pan Zhenhua1，Pan Jianfeng1，Dong Gang2

（1．School of Energy and Power Engineering，Jiangsu University，Zhenjiang212013，Jiangsu，China；2．Key Laboratory of Transient Physics，Nanjing University of Science and Technology，Nanjing210094，Jiangsu，China）

A computational study of the interaction between shock waves and a spherical flame was carried out using the ninth－order WENO and the tenth－order central difference schemes，and the influence of shock intensity and flame size on the interaction process was investigated．It can be found from the results of our study that the increase of the shock intensity and the flame size can both induce detonation in the flow field，but the influence of the shock intensity is relatively stronger．Further，the detonation induced by shock wave can lead to quick flame expansion and increase its heat release rate，thereby affecting the combustion characteristics．Besides，the detonation wave will quickly burn out the combustible gas，merge the previously existing reflected shock waves in the propagation process，and form local high pressure zones，which can significantly alter the flow field structure．

shock wave；flame；detonation；flow field structure

O381国标学科代码：13035

A

10．11883／1001－1455（2017）04－0741－07

（责任编辑 王玉锋）

2015－12－23；

2016－05－03

国家自然科学基金项目（11402102，11372140）；江苏省自然科学基金青年项目（BK20140524）；江苏省博士后基金项目（1402013B）；江苏大学高级专业人才科研启动基金项目（14JDG031）

朱跃进（1986－ ），男，博士，讲师，zyjwind＠163．com。