饶飞雄， 丁 珏， 翁培奋， 李孝伟

(上海大学上海市应用数学和力学研究所，上海 200072)

脉冲爆震发动机(pulse detonation engine, PDE)具有设计简单、热效率高、应用范围广等优点，未来有替代传统发动机而成为新型动力装置的趋势，因此越来越受到国际研究者的高度重视[1]。

爆燃波通过火焰区的质量扩散和热扩散使前方反应物着火，从而实现火焰传播[2]。而爆轰波则是一道超声速压缩激波，它通过前导激波扫过混合物时的绝热压缩加热点燃混合物。因反应过程接近定容燃烧，爆轰波产生的能量与效率远高于普通燃烧[3]。通常，研究者采用一步、两步[4-6]反应模型模拟气相爆轰波结构。张薇等[7]采用单步化学反应模型研究气相爆轰波胞格尺寸和点火延迟时间关系；张德良等[8]采用简化二阶段化学反应模型模拟爆燃转爆轰和爆轰波在斜面上的马赫反射；刘云峰等[9]采用二步化学反应模型分析爆轰波的三波点和发展过程。由于一步、两步反应模型简化了真实爆轰过程链式反应的复杂性和影响，难以给出模型中相关参数，且反应产物组成不完全导致爆轰参数产生较大的偏差，进而影响波阵面结构及其相关参数。因此，越来越多的人采用较能正确反映爆轰波化学反应本质的基元反应。Oran等[10]采用基元反应模型对二维爆轰波进行模拟来研究能量的释放规律；Fukuda 等[11]采用详细反应机理模型研究了壁面条件对爆燃转爆轰的影响；董刚等[12]基于氢气与空气详细化学反应模拟爆轰波，分析了爆轰反应后的组分浓度和温度分布；王昌建等[13]采用基元反应数值模拟研究了初始扰动对爆轰波形成的影响；张宝亮等[14]基于有限体积法分析了障碍物与爆轰波的相互作用规律。

本工作基于敏感性分析，将二十步反应简化为九步反应模型，并将详细基元反应模型应用于爆轰波数值模拟，分析爆轰波结构、传播性质及胞格。同时，本工作研究了Ar稀释和初始压强的变化对爆轰波参数的影响, 得出其影响爆轰波特性的相关规律，为深入研究爆震发动机非均相爆轰特性奠定了理论基础。

1 爆轰波形成和发展的控制方程和数值计算方法

1.1 控制方程

爆轰波包括超声速激波以及紧随的化学反应区组成的高速燃烧波。对于爆震室可燃气体爆轰及爆轰波传播过程，采用二维Euler 控制方程：

式中：u是守恒变量，Fi是对流通量，Q为源项，分别为

其中δ是克罗内克函数，ρ是密度，u1、u2是速度的分量，p是压力，为能量方程的源项，是第K个组分的质量生成率，N是化学组分数目。

YK是第K个组分的质量分数，满足

能量E定义为：

式中：hs为显焓。

化学反应速率用Arrhenius 公式来计算

式中：Ea为化学反应活化能，A为指前因子，R为通用气体常数，T为流场温度。

1.2 计算区域的网格及初始边界条件

爆震室直径为2 cm，充满初始压力为一个标准大气压、初始温度为300 K 的氢气和氧气混合气体，其中气体组分的物质的量之比为m(H2):m(O2)=2:1。爆震室左端点火，示意图如图1 所示。

图1 爆轰室示意图Fig.1 Sketch map of detonation chamber

爆轰的反应速率是普通燃烧的103～108倍[3], 为了准确地模拟爆轰波结构，网格尺寸需要满足精度要求。因此，基于CJ 理论[15]得出的爆轰相关参数，本工作采用4 种不同的网格尺寸数值模拟爆轰波形成和发展过程，并与实验和爆轰的CJ 理论值进行比较，结果见表1。

表1 网格尺寸对爆轰参数的影响Table 1 Influence of mesh size on detonation parameters

表1 中数据显示：随着计算网格尺寸的减小，爆轰波传播速度较接近实验值，但计算量成倍增大；爆轰压强变化不大，因有效小数的限制，压强在1.911 MPa 附近变化。因此，在满足计算精度的基础上考虑计算量，本工作采用0.10 mm 的计算网格。

1.3 数值计算方法

1.3.1 爆轰的超声速流场

在超音速流动中，基于密度求解器能在求解Riemann 问题时更能彰显其优越性。针对网格i+交界面，本工作计算对流通量采用Roe 格式[16]：

c为音速，U=nxu+nyv为界面法向速度，nx和ny为界面法向向量在坐标轴上的分量。另外，为了保证计算结果收敛与稳定性，时间步采用10-8s。

1.3.2 化学反应动力学

对于激波与紧随化学反应区的爆轰波，需要求解反应区中多组分的化学反应。本工作考虑m(H2):m(O2)=1:0 的8 组元二十步基元反应H2/O2模型[17]，具体参数如表2 所示。

表2 化学反应和反应参数Table 2 Chemical reactions and reaction parameters

2 九步基元反应模型

带有反应流的数值模拟中采用详细基元反应机理，计算量增大，同时基于详细反应机理的反应流的计算常常出现刚性问题，给数值模拟带来了较大的困难。因此，本工作根据敏感性分析方法[18]，分析反应机理中参数变化对计算结果的影响程度，消除详细机理中对系统主要参数值影响较小的基元反应和组分。

2.1 敏感性分析法

敏感性分析分为局部敏感性分析和总体敏感性分析。局部敏感性主要用于研究基元反应的动力学特性。局部敏感性分析正交表达式为

总体敏感性分析表达式为

式中：Bcj和Bfj分别为第j个反应对N种组分生成速率和浓度的总体正交敏感性系数。本工作主要从整体研究部分基元反应对主要计算结果影响程度，因此采用总体敏感性分析。

2.2 基元反应的敏感性分析

爆震室预混物中m(H2):m(O2)=2:1，对二十步基元反应进行敏感性分析，结果如图2 所示。敏感性系数为正数代表此基元反应对此组分起促进作用，负数代表此基元反应对此组分进行消耗。图2 显示：基元反应3 是影响整个系统反应速率的关键基元反应；对整个反应系统影响较大的基元反应序号为1～5、7、12～13 和16。

图2 总体敏感性分析Fig.2 Overall sensitivity analysis

将通过敏感性分析所得出的九步基元反应用于爆轰波数值模拟，计算所得的爆速、温度与压强与实验值的误差分别为0.236%、2.76%和4.57%，而爆轰波传播1 μs 所需计算时间平均为50.556 min。

3 二十步基元反应模型

3.1 爆轰波参数

表3 为网格为0.1 mm 二十步模型和九步模型数值计算结果与实验、CJ 理论的比较。由表3 可知：九步简化模型和二十步基元反应模型计算所得的爆速、温度和压力与实验值、理论结果较为一致。二十步模型数值计算所得的爆速、温度与压强与实验值的误差分别为0.2%、2.78%和2.74%，表明了本工作所采用的模型和数值计算方法的合理、正确性。

表3 基元反应模型模拟爆轰波参数的对比Table 3 Comparison of simulated detonation parameters in elementary reaction models

相比于九步简化模型，二十步模型数值计算时间略长。爆轰波传播1 μs 所需计算时间平均为58.25 min。

3.2 爆轰波瞬态结构

图3 爆轰场压力、温度及OH质量分数的分布云图Fig.3 Distribution of detonation pressure, temperature and OH mass fraction

图3给出了爆轰波传播中，三波点在碰撞过程压力、温度和OH 离子质量分数的分布情况，其中D 为未反应气囊，TP 为三波点，Hotspot 为热点。数值模拟显示，爆轰波传播中具有典型的三波点结构：入射激波(I)，横波(TW)，马赫干(M)。

两个三波点碰撞前的爆轰流场(见图3(a))。非平面的诱导激波波阵面是由入射激波和马赫干交替组成，入射激波强度比马赫干弱。OH 质量分数场显示：化学反应区与马赫干耦合在一起向前传播，而与入射激波之间有一段明显的距离，表现出化学反应区与入射激波解耦的性质。马赫干、入射激波和横波交汇处形成三波点，且爆轰场形成了多个未反应的气囊。

两个三波点碰撞时的爆轰流场(见图3(b))。三波点相互碰撞，形成热点，导致该位置的温度和压力急剧升高。而此处的未反应气囊中，化学反应速率急剧加快，释放大量的能量，即为二次反应。二次反应释放能量对爆轰波自持传播起了重要作用。

三波点碰撞后的爆轰流场(见图3(c))。由于三波点的二次反应产生大量的能量，未反应气体和产物被挤压形成向前的射流。由于不稳定性，射流在其头部形成蘑菇状的结构。先前的入射激波转变为马赫干，原先的马赫干退化成入射激波。强度较大的马赫干与反应区发生强耦合作用。

3.3 爆轰场中单个胞格的传播性质

爆轰波传播中呈现典型的三波点结构：马赫干、入射激波、横波。横波与三波点周期性运动形成胞格结构。图4 给出了41.7 us 时刻爆轰场单个胞格的结构图。图中显示，爆轰波自持传播并产生了周期性的爆轰胞格，在模型宽度方向上爆轰胞格不完全对称。

图4 单个胞格图(单位: m)Fig.4 Single detonation cell pattern (unit: m)

图5 为沿着单个胞格中心线的压力分布，其中展示了3 个完整胞格结构。三波点相互碰撞导致化学反应速率增加，该位置的压力和温度升高，随后中心线上压力逐渐衰减，马赫干衰退为入射激波，直至下一次三波点碰撞，压力又达到峰值。计算初始点火位置形成的膨胀波，经过壁面反射形成了如图4 的反射波，波的相互作用导致流场中x轴上的0.053 m 至0.062 6 m之间区域压力变化平缓，形成压力平台。随着反射激波的发展，压力逐渐衰减。

图5 沿着单个胞格中心线的压力图Fig.5 Pressure profile along cell centerline versus distance

4 惰性气体Ar 的稀释作用

为了研究惰性气体Ar 的稀释作用，在相同氢氧混合气体中添加质量分数为70%的Ar。计算得到的稳定爆轰波爆速，温度和压强分别为1 639.676 m/s, 2 933.018 K，1.958 MPa。将未添加与添加了惰性气体的稳定爆轰参数进行对比，结果如图6 所示。可以看到：稀释后的爆轰波参数，如爆轰压力和温度峰值分别降低0.047 MPa 和749.738 k。分析其原因，主要是添加了惰性气体Ar 后，延缓了反应的释放，爆轰场温度变化显著。而压强变化不大，主要是惰性气体并未参与化学反应所致。

图6 Ar 稀释对爆轰波参数的影响Fig.6 Influence of Ar dilution on detonation parameters

5 爆震室初始气相压强的影响

在真实环境中，初始气相压强不恒定，因此研究不同初始压强对爆轰波参数的影响十分必要。考虑气体初始压强分别为101 325、80 000 和60 000 Pa 的情况，图7 为物质的量之比为1:1 时的氢氧混合物爆轰参数的分布，数值计算结果如表4 所示。可以看到：爆速、温度和压强随着初始压强的降低而降低，即初始气相压强每下降1 Pa，爆轰波峰值压力平均下降约18.51 Pa；而初始压强对流场温度变化影响不太明显。

图7 初始压强对爆轰参数的影响Fig.7 Influence of initial pressure on detonation parameters

表4 初始压强对爆轰波参数的影响Table 4 Influence of initial pressure on detonation parameters

6 结束语

脉冲爆震发动机爆震室气相爆轰形成与发展是气体流动和化学反应同时发生且强烈耦合的过程。本工作建立九步简化模型和二十步基元反应模型，对气相爆轰过程开展数值模拟，并将数值结果与实验数据、理论计算结果进行对比和分析。

(1) 二十步基元反应较为准确地反映流场中爆轰波精细结构：马赫干、入射激波、横波.在爆轰波传播中，入射激波与反应区是分离的，而马赫干与化学反应区紧紧耦合在一起；三波点的碰撞导致化学反应速率增加释放的能量，为爆轰波持续传播提供重要的支持。横波和三波点周期性运动形成了爆轰波的胞格结构。

(2) 化学反应动力学影响着爆轰波的结构和爆轰参数。通过总体敏感性分析，对二十步基元反应简化得出九步基元反应。并将两种基元反应模型应用于爆轰波模拟中，得出：基于二十步基元反应模拟计算的爆速、温度与压强与实验值的误差较小，分别为0.2%、2.78%和2.74%。九步简化模型的优势是计算相比详细基元反应，爆轰波传播1 μs 所需的计算时间可节省13.21%。

(3) 爆轰波传播过程中，惰性气体Ar 延缓了化学反应能量的释放，使得爆轰波传播速度和流场温度降低。此外，气相初始压强的减少使得单位体积气相燃料释放的热量降低，因此爆速、温度和压强也随之降低。计算结果显示，初始压力每下降1 Pa，则爆压峰值下降约18.51 Pa。Ar 稀释和初始压强的变化对爆轰波参数的影响和所得相关规律，为深入研究爆震发动机非均相爆轰特性奠定了理论基础。