激波诱导火焰变形、混合和燃烧的数值研究＊

朱跃进，董 刚，2，刘怡昕，范宝春，蒋 华

（1.南京理工大学瞬态物理重点实验室，江苏 南京210094；2.北京理工大学爆炸科学与技术国家重点实验室，北京100081）

激波掠过预混火焰界面时，会诱导界面出现大尺度的Richtmyer－Meshkov（RM）不稳定和小尺度的Kelvin－Helmholtz（KH）不稳定，导致火焰界面发生变形，进而强化未燃气与可燃气的混合，并促进流场湍流化，引起燃烧放热的增强。这一系列复杂的物理化学现象不仅会发生在某些自然现象（如超新星爆炸［1］等）中，而且还常常在超燃推进［2］、工业爆炸灾害［3］等出现，因此有关激波诱导的火焰变形、混合和燃烧等特性的研究具有重要意义。

作为可压缩化学反应流动中的一个基本问题，平面激波与火焰界面相互作用的研究得到了一定程度的开展。G.H.Markstein［4］首先开展了激波管内激波与火焰相互作用的实验研究，显示了火焰在弱激波及反射激波作用下的变形过程；V.T.Ton等［5］考察了化学反应对激波－射流火焰相互作用的影响，由于采用的激波马赫数较小（Ma＝1.093），因而化学反应对火焰的影响并不大，火焰变形过程与惰性介质密度界面的变形过程类似；Y.Ju等［6］研究了不同入射激波强度对火焰变形的影响规律，结果显示，激波强度的增加可显著增加火焰界面的长度，使已燃气和未燃气接触面增加，进而提高火焰的燃烧速率。A.M.Khokhlov等［7］采用带化学反应的二维Navier－Stokes（NS）方程对入射激波与火焰的单次作用过程进行了数值研究，结果表明，RM不稳定是火焰变形的主要机制，变形后的火焰与周围未燃气的混合强化，其表面积增加，燃烧速率和放热速率均有所提高。但他们也指出，入射激波与火焰的单次作用对火焰变形和混合、燃烧的强化是有限的。接着，他们采用相同的数值方法，进一步研究了激波管中入射激波及其在管道封闭端形成的反射激波与火焰多次作用的过程，计算发现层流火焰在激波的多次扰动下高度变形并形成湍流火焰，其火焰表面积明显增加，燃烧放热率显著提高［8－9］，这个结果也得到了G.O.Thomas等的实验证实［10］。G.Dong等［11］采用数值模拟方法分析了激波与火焰作用过程，考察了不同化学反应机理和网格尺寸对结果的影响。

由上述可以看出，激波与火焰界面的相互作用涉及火焰变形、混合与燃烧放热等基本过程，它们之间存在相互影响。一般认为，激波诱导火焰引起的变形可以强化火焰两侧未燃气与已燃气的混合，从而促进燃烧。但是，目前不仅缺少混合和燃烧两者之间关系的详细描述，而且也缺乏对两者关系的进一步认识。为了深入细致阐明激波诱导的火焰变形以及由此带来的混合与燃烧之间的相互作用关系，同时避免三维计算带来过大的计算量，本文中，对平面入射激波及其反射激波诱导球形火焰变形的现象进行了二维数值研究，通过定义有效的量化参数，分析流场中变形火焰的演变过程、火焰几何量变化、燃烧放热规律、混合与燃烧的作用关系等特性，以深入理解激波与火焰作用的基本现象和规律。

1 数理模型与实验验证

1.1 数理模型

为描述激波与火焰的相互作用过程，采用的二维带化学反应的Navier－Stokes方程表示：

式中：ρ为密度，ui为i方向速度（i＝1，2），p为压力，τij为黏性应力张量，E 为体积总能量，E＝p／（γ－1）＋0.5ρ∑2i＝1u2i＋ρqY，p为压力，γ 为比热比，q为质量总化学能，Y 为反应物质量分数，qj为热通量，qj＝－k∂T／∂xj，ω为化学反应速率，热扩散率k、扩散系数D、运动黏度ν可表示如下［3］：

式中：ν0、D0、k0是常数，T 为温度，n＝0.7，本文中设路易斯数、普朗特数和施密特数均为1，则可选取常数ν0＝D0＝k0＝3.2×10－7kg／（s·m·K0.7）。

控制方程（1）～（4）采用分裂算法求解。空间导数项的黏性部分采用二阶空间中心差分计算，无黏部分采用高精度五点TVD格式［12］计算，时间推进过程采用二阶Runge－Kutta方法求解。

采用单步化学反应模型描述体系中的燃烧过程，化学反应热力学和动力学参数的选取见文献［13］，反应速率可表达为：

式中：A 为指前因子，取A＝1.2×108m3／（kg·s），Ea为活化能，取Ea／RT0＝38.2，T0代表预混气的初始温度，R为通用气体常数。预混气的化学反应质量放热量q和比热比γ分别取qM／RT0＝39.76和γ＝1.258，其中M为预混气的分子量，这组参数能同时反映火焰燃烧和爆轰的热力学特性［13］。

图1 计算区域和初始流场Fig.1 Computational domain and initial flow field

1.2 初始及边界条件

根据文献［10］中的实验，设计如图1所示的计算区域与初始流场。计算区域流向（x方向）长170mm，法向（y方向）高38mm；反应性预混气体组成为C2H4＋3O2＋4N2，密度ρ0＝161.5g／m3，初温T0＝293K；初始球形火焰半径R0＝19mm，火焰中心（35mm，0），火焰区内密度ρ1＝15.78g／m3，火焰内外压力均为p0＝13.3kPa；入射激波马赫数1.7，初始时位于x＝12mm处，激波沿x方向从左向右传播，激波后气体状态由Rankine－Hugoniot关系式给出。计算区域上边界（y＝38mm）和x方向右端面为管壁，均采用无滑移的刚性绝热壁面边界；y＝0为对称面，采用对称边界；x方向左端面处采用零梯度边界。计算使用的均匀网格尺寸为47.5μm，由于预混气体在本文条件下的火焰面厚度为约0.785mm［11］，因此火焰面厚度可由约16个网格刻画，这个网格分辨率足以反映火焰界面的变化过程。此外，计算时间步长取8.55ns，使库朗数小于0.2，满足计算中时间推进的稳定性判据。由于本文中使用了较高的网格分辨率，计算量较大，因此采用自己发展的并行化设计的计算程序完成。

1.3 实验验证

数值计算算法和化学反应模型的可靠性，已在文献［13］的计算中予以验证。考虑到网格精度的进一步提高，依据文献［10］的实验结果对数值模型和化学反应参数做进一步验证。图2给出了各时刻激波阵面距点火位置的距离，以及激波诱导的变形火焰的宽度和高度，其中，计算和测量结果均采用火焰初始半径R0进行量纲一化处理。从图2可以看出，激波阵面距点火位置的距离与实验结果很好吻合，而火焰尺寸的变化规律也与实验结果一致，说明本文计算方法和化学反应各参数是合理和可靠的。

图2 数值计算结果与实验结果的对比Fig.2 Comparisons between calculated and experimental results

2 结果与讨论

为深入认识激波与火焰的相互作用关系，主要从火焰演变特性（可视化结果）、火焰二维效应（积分量）、火焰分形维效应（界面形态）以及混合与燃烧的关系等方面予以讨论和分析。

2.1 火焰演变特性

图3～4分别给出了入射激波和反射激波与火焰作用的流场图像，其中每幅子图的上半部分对应密度分布，下半部分对应涡量分布。由图3可见，当入射激波与火焰开始作用时，由于火焰表面的压力梯度和密度梯度方向不一致，产生斜压效应，进而在火焰表面形成涡量，因此初始时火焰迎风面发生凹陷，从而表现出RM不稳定诱导的大尺度变形（见图3（a）～3（b））；随着时间，由于火焰界面内外存在速度差，因而在火焰表面逐渐形成了KH不稳定诱导的小尺度离散涡结构（见图3（c））；接下来，在入射激波作用后的一段时间内，产生的大尺度涡量逐渐发展至火焰中心区，而火焰界面处主要以小尺度涡量为主；由于界面处形成的涡量能卷吸环境中的未燃气体，促进燃烧，因此火焰的面积明显增加（见图3（d）～3（f））。

当入射激波行至右端壁面后发生反射，反射激波与已变形的火焰再次作用，这使得火焰变形加剧。由图4可见，在反射激波与火焰作用初期，火焰表面变形加剧，弯曲皱褶增多，小尺度涡结构也明显增加（见图4（a）～4（b））；当反射激波完全扫过火焰后，反射激波的诱导作用使变形火焰向壁面翻转（见图4（c）），这与实验结果［10］一致；反射激波波后滞止的高温高压环境为火焰进一步燃烧膨胀提供了基础，随着时间，火焰面积迅速增加，变形火焰逐渐靠近上下壁面并沿流向方向膨胀（见图4（d）～4（f））。

图3～4还显示了激波的运动过程。由于火焰内密度较小，所以入射激波在火焰内部传播速度较快，能在火焰表面形成激波分歧，而反射激波在入射激波沿壁面诱导产生的边界层内也发生了激波分歧。上述数值计算结果合理地描述了激波诱导火焰变形的过程，计算结果也与实验很好符合，这为下面进一步的分析提供了可信的结果。

图3 入射激波与火焰作用的密度和涡量场Fig.3 Densities and vorticities for the interaction between incident shock wave and flame

图4 反射激波与火焰作用的密度和涡量场Fig.4 Densities and vorticities for the interaction between reflected shock wave and flame

2.2 火焰二维效应

采用积分量进行分析，可以从变形火焰整体发展行为考察其演变过程。采用积分形式定义火焰有效面积A和平均化学反应放热率dh／dt，用来考察火焰变形的过程：

式中：D和F分别代表整个计算区域和火焰区（定义为反应物质量分数Y≤0.99），X代表反应物的体积分数，其中式（7）采用的火焰有效面积表征了整个计算域内去除混合影响之后的火焰面积。

图5给出了量纲一火焰有效面积A／A0（A0为初始火焰面积）和平均化学反应放热率dh／dt随时间的变化过程。可见，火焰有效面积的变化经历了4个阶段：阶段Ⅰ为入射激波压缩阶段（0～94μs），与此对应的火焰有效面积在压缩作用下逐渐变小；阶段Ⅱ为火焰膨胀阶段（94～328μs），此时激波已掠过火焰，燃烧放热使火焰开始膨胀，面积逐渐增加；阶段Ⅲ为反射激波压缩阶段（328～374μs），此时变形火焰在反射激波的压缩下面积再次变小；阶段Ⅳ为火焰再次膨胀阶段，由于反射激波波后介质处于高温高密度的滞止状态，这为化学反应速率的提高提供了有利条件，所以火焰燃烧膨胀加剧、面积迅速增加。图5中平均燃烧反应放热速率的变化显示，压缩阶段（阶段Ⅰ和Ⅲ）放热速率呈现持续增长的趋势，说明激波压缩作用提供的高温高密度环境可以促进化学反应的进行，尤其反射激波与变形火焰的再次作用可以显著提高火焰的放热速率；入射激波作用后的火焰膨胀阶段（阶段Ⅱ），放热速率则基本维持在一个恒定值，而反射激波作用后的火焰膨胀阶段（阶段Ⅳ），放热速率持续上升到一定值后维持一段时间随后又开始下降，尽管如此，反射激波作用后（阶段Ⅳ）的放热速率比入射激波作用后（阶段Ⅱ）的放热速率增加约一个数量级，这反映了反射激波对燃烧过程的巨大促进作用。

图5 火焰的有效面积和平均反应放热率Fig.5 Effective area and average reaction heat release rate

2.3 火焰分形维效应

图3～4表明，火焰界面在激波的作用下呈现皱褶和小尺度涡的形态。由于预混火焰的燃烧过程主要发生在未燃气与已燃气的界面处，为了反映火焰界面的形态和几何自相似行为对燃烧过程的影响，对数值模拟得到的二维计算结果进行图像处理，获得了火焰界面量纲一长度L／R0和火焰界面平均分形维数Dm随时间变化的结果，如图6所示。其中，分形维数Dm的处理采用了基于质量分数Y＝0.99的火焰界面的数盒子处理方法，根据不同盒子尺寸ε，可以获得覆盖火焰界面的不同盒子个数N（ε），由此可计算二维火焰界面分形维数：该量的变化反映了火焰表面小尺度涡诱导的界面皱褶程度。图6中，误差棒代表不同盒子取样位置所确定的分形维数的最小和最大值，阴影区所处的分形维范围表征充分发展湍流中的几何自相似性。

图6 火焰界面长度和分形维数Fig.6 Interface length and fractal dimension

由图6可以发现，火焰界面长度变化并不完全遵循图5中的火焰有效面积4个阶段变化规律。反射激波作用前（阶段Ⅰ和Ⅱ），界面长度持续增长，而界面分形维数在阶段Ⅰ和阶段Ⅱ前期（t＜150μs）也快速增长（Dm由1.04增至1.32），这表明，入射激波作用引起的火焰小尺度皱褶占据主导；到阶段Ⅱ后期（150～327.9μs），火焰界面分形维数维持在1.29～1.33之间，表明火焰皱褶程度不再明显增加，而此时火焰界面长度仍在持续增加，由于此阶段火焰膨胀并不明显且反应放热速率维持在较小值（见图5），这意味着入射激波诱导的大尺度火焰变形拉伸使火焰界面长度继续增加。因此，在反射激波与火焰作用之前，物理过程对火焰的变形和皱褶起主导作用。当反射激波与变形火焰再次作用时（t＝328μs），火焰界面长度经短暂的减小（阶段Ⅲ）后，又持续增长至阶段Ⅳ前期（t＝450μs），而这一过程中，界面分形维数又经历一个持续上升的过程（Dm由1.33增至1.46）。由于此阶段化学反应放热速率急剧增加至最大值（见图5），可以断定，激波诱导的火焰表面皱褶程度的增加和燃烧导致的火焰膨胀，共同促进了界面长度的增长。然而，从阶段Ⅳ后期开始（t＞450μs），火焰界面长度不再增加，界面分形维数开始下降，这说明由于界面皱褶程度的减少导致界面长度有减小的趋势，但是此阶段火焰仍在膨胀，又促使了界面的增长，因此综合的效果是火焰界面长度维持在一定的范围。阶段Ⅳ后期的变化表明，燃烧过程对火焰的皱褶有一定的抑制作用，这主要体现在后期火焰界面皱褶程度的降低，由于皱褶的减少会减小火焰界面未燃气和已燃气接触的长度，因此反过来又会抑制反应放热过程，这个特点在反应放热速率阶段Ⅳ后期的下降（见图5）中可以得到体现。

为比较火焰在不同膨胀阶段（阶段Ⅱ和Ⅳ）的发展规律，表1给出了火焰有效面积A、火焰特征长度l和火焰界面长度L等几何量在各阶段的平均增长率，其中，特征长度定义为计算区域内半圆形火焰理想状态下（即不受变形影响）的周长。表中，ηA＝ΔA／Δt，ηl＝Δl／Δt，ηL＝ΔL／Δt，ηⅣ／ηⅡ代表了阶段Ⅳ和Ⅱ火焰各几何量平均增长率的比，用于考察反射激波作用前后的影响，而ηL／ηl则表征了火焰界面在激波作用后的变形和皱褶程度。由此可以看出，反射激波作用后的火焰各几何量增长均比作用前的几何量增长有成倍增加，而在阶段Ⅱ和Ⅳ，ηL／ηl分别达到11.57和13.21，这说明激波导致火焰变形和皱褶现象十分显著，尤其在阶段Ⅳ前期，这与图3～4反映的规律相同。

表1 阶段Ⅱ和Ⅳ火焰几何量平均增长速率的变化Table 1 Variation of the average increase rate of geometrical parameters in stagesⅡandⅣ

图6中的阴影区域代表了文献［14－15］给出的RM不稳定诱导的惰性界面湍流混合的分形维数范围（Dm＝1.3～1.4），在此范围内，变形界面表现出具有充分发展湍流特性的几何自相似特征。计算结果表明，激波诱导的火焰变形在阶段Ⅱ膨胀后期（t＞150μs）已经进入该范围，并且在反射激波作用后甚至超过该范围，尽管后期Dm有所下降，但仍未低于该范围，这说明激波诱导的火焰于入射激波作用后不久便呈现了湍流燃烧的形式。

2.4 混合与燃烧的关系

预混火焰不同尺度的变形增加了火焰界面两侧未燃气和已燃气的接触表面，可以强化两者之间的混合过程，这对于促进燃烧过程是有利的。但是上述结果也表明，在反射激波作用的后期，燃烧过程会弱化火焰表面的皱褶程度，因而在一定程度下抑制了混合过程。在有关界面混合的讨论中，J.Yang等［16］首先引入拉伸率概念定量表征混合过程，并考察了激波强度、密度比、界面几何构型等对混合的影响；S.Kumar等［17］沿用前者的思路，对激波与5种不同构型的SF6气柱界面的作用进行了实验研究，发现初期（t≤220μs）质量线的拉伸率呈指数变化，但这种方法不能描述后期的混合情况；J.H.J.Niederhaus等［18］针对激波与不同密度气泡界面作用的三维数值模拟结果，定义了气泡区域内环境气体所占的体积分数为混合度，用于描述混合。但是上述研究仅仅考虑了混合作用，而在激波与火焰界面作用的过程中，燃烧放热会与混合相互影响，但目前尚还未见到对混合与燃烧相互关系的具体描述。因此，为深入分析这两者的关系，本文中定义一个量纲一数η，用以表征火焰区域内总质量变化率与燃料消耗速率的比，表达如下：

式中：分子项代表火焰区（Y≤0.99）未燃气体混入速率，分母项代表火焰区未燃气体燃烧消耗速率。因此，当η＞1时，表示混合速率大于燃烧速率，混合过程占主导地位，当η→1时，则表示进入火焰区的未燃气速率接近燃烧速率，这代表了以火焰传播为主导的燃烧模式。

图7 火焰区域内总质量变化率与燃烧消耗速率的比Fig.7 Ratio of the total mass change rate to the burning consumption rate within the flame

图7显示了η随时间的变化过程，可以看出，η在阶段Ⅰ初期混合的影响逐渐下降、燃烧的影响逐渐增强，但混合过程仍占主导；在阶段Ⅰ后期，由于入射激波与火焰界面相互作用产生涡量，使得混合大大增强，所以η又开始上升；在随后的各阶段中，η缓慢震荡降低，表明火焰变形同时受到混合与燃烧作用的影响，且燃烧作用慢慢增强；在阶段Ⅳ后期，该值基本在1附近震荡，说明反射激波后的火焰发展主要以火焰传播的方式进行，混合作用已经很小。上述结果表明，未燃气体与可燃气体的混合会促进燃烧，但逐渐增强的燃烧作用能迅速地消耗混入火焰内的未燃气，反过来抑制混合。

3 结 论

采用带单步化学反应的NS方程和高网格分辨率，对平面入射激波及其反射激波诱导火焰变形过程进行数值研究，利用实验结果验证了计算模型和参数的可靠性。提出并采用流场可视化、火焰积分量、火焰界面几何特性等分析手段探讨了火焰变形、燃烧放热和混合等过程的变化规律及其相互作用关系。结果表明，初始球形火焰在激波及其反射波的各次作用下均表现出先压缩再膨胀的特点，反射激波作用前，激波的诱导作用对火焰变形和表面皱褶起到了主要作用，此阶段以物理过程为主导；反射激波作用后，变形火焰放热速率增加约1个数量级，表征火焰形态的各几何量均成倍增加，主导火焰发展的上述物理过程逐渐由化学反应过程所替代，尤其是在火焰发展后期，化学反应的增强会反过来抑制火焰界面的皱褶，进而弱化了界面两侧已燃气和未燃气之间的混合，使火焰发展主要受燃烧过程控制而非混合过程控制。

［1］Maran S P，Sonneborn G，Pun C S J，et al.Physical conditions in circumstellar gas surrounding SN 1987A12years after outburst［J］.Astrophysical Journal，2000，545（1）：390－398.

［2］Marble F E，Hendrick G J，Zukoski E E.Progress toward shock enhancement of supersonic combustion process［R］.AIAA 87－1880，1987.

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

［4］Markstein G H.A shock－tube study of flame front－pressure wave interaction［C］∥6th Symposium （International）on Combustion.Pittsburgh，USA：The Combustion Institute，1957：387－398.

［5］Ton V T，Karagozian A R，Marble F F，et al.Numerical simulations of high speed chemically reactive flow［J］.Theoretical and Computational Fluid Dynamics，1994，6（2／3）：161－179.

［6］Ju Y，Shimano A，Inoue O.Vorticity generation and flame distortion induced by shock flame interaction［C］∥27th Symposium （International）on Combustion.Pittsburgh，USA：The Combustion Institute，1998：735－741.

［7］Khokhlov A M，Oran E S，Chtchelkanova A Y，et al.Interaction of a shock with a sinusoidally perturbed flame［J］.Combustion and Flame，1999，117（1／2）：99－116.

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

［9］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.

［10］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.

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

［12］Chakravarthy S R，Osher S.A new class of high accuracy TVD schemes for hyperbolic conservation laws［R］.AIAA 85－0363，1985.

［13］朱跃进，董刚，范宝春.受限空间内激波与火焰作用的三维计算［J］.推进技术，2012，33（3）：405－411.Zhu Yue－jin，Dong Gang，Fan Bao－chun.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.

［14］Vorobieff P，Rightley P M，Benjamin R F.Shock driven gas curtain：Fractal dimension evolution in transition to turbulence［J］.Physica D：Nonlinear Phenomena，1999，133（1）：469－476.

［15］Ng H D，Abderrahmane H A，Bates K R，et al.The growth of fractal dimension of an interface evolution from the interaction of a shock wave with a rectangular block of SF6［J］.Communications in Nonlinear Science and Nu－merical Simulation，2011，16（11）：4158－4162.

［16］Yang J，Kubota T，Zukoski E E.Applications of shock－induced mixing to supersonic combustion［J］.AIAA Journal，1993，31（5）：854－862.

［17］Kumar S，Orlicz G，Tomkins C，et al.Stretching of material lines in shock－accelerated gaseous flows［J］.Physics of Fluids，2005，17（8）：082107.

［18］Niederhaus J H J，Greenough J A，Oakley J G，et al.A computational parameter study for the three－dimensional shock－bubble interaction［J］.Journal of Fluid Mechanics，2008，594：85－124.