脉冲爆震燃烧室内缓燃向爆震转捩数值模拟

赵 炜，沈盈盈，杨福军，张培红，吴晓军

（中国空气动力研究与发展中心计算空气动力研究所，四川 绵阳 621000）

0 引言

脉冲爆震发动机（Pulse Detonation Engine，PDE）是一种间歇触发爆震波产生推力的增压燃烧系统[1-2]，具有热效率高、结构简单、宽速域的特点，可以作为传统涡轮基为主的航空发动机的低成本替代动力装置，受到各国研究人员的广泛关注[3-5]。

理想化的PDE 燃烧室可以看作是一端封闭、一端开口的带有障碍物的管道，可燃混气被点燃后通过一系列的火焰加速机制最终完成缓燃向爆震的转捩（Deflagration to Detonation Transition，DDT）。由于PDE 周期性工作的特点，循环频率保持在100 Hz 以上才能使推力准稳态化[6]，这就要求在毫秒量级的时间内完成1 个工作循环，因此如何强化火焰加速，缩短DDT 过程是PDE 工程应用的重要课题。具体来讲，是如何在更短的距离与时间内完成火焰加速过程最终达到爆震触发的临界条件。

缩短DDT 过程的1 个有效手段是在管道中增加周期性的障碍物。目前，常见的障碍物构型有环形孔板[7]、楔 形钝 体[8]及Shchelkin 螺旋[9]。Lee 等[10]指 出障碍物的堵塞比为0.3～0.6时，促进DDT过程的效果较为显著；Ciccarelli 等[11-13]对障碍物管道内的火焰加速相关试验研究进行了系统总结，根据火焰传播速度以及影响火焰传播的主要因素的不同，可以把火焰加速过程分为层流火焰传播、湍流火焰传播、激波与火焰相互作用、爆震触发与传播4 个阶段。整个火焰加速过程受当地的几何条件、化学反应、湍流、激波等多种因素的影响，并耦合复杂的火焰面/旋涡/障碍物/激波相互作用。由于受试验测量技术的限制，大量研究主要对沿程的压力、火焰面历程与局部的光学图像等宏观参数进行测量，缺乏对上述整个火焰加速全程的精细观测，难以深入理解障碍物管道内的火焰加速内在机制。

利用CFD 技术，国内外学者针对火焰传播过程与DDT 过程的相关机理开展了深入的数值模拟研究。Boeck 等[14-15]通过平面激光诱导氢氧基荧光技术（Planar Laser Induced Fluorescence of the OH molecule，OH-PLIF）测量与2维数值模拟方法研究了预混火焰在单组、多组障碍物下的传播特性，发现火焰在障碍物上缘脱落的旋涡对火焰面的结构影响较大；Emami 等[16]利用2 维大涡模拟（Large Eddy Simulations，LES）研究了化学恰当比下的氢气/空气混合物在障碍物管道内的火焰加速过程，表明在慢速火焰阶段，火焰/旋涡的相互作用导致火焰面产生褶皱、折叠（本质上是增加火焰面积），而在快速火焰阶段，火焰加速的主要机制是反复激波/火焰相互作用；Gamezo等[17-18]采用单步总包化学反应的数值模拟结果显示：单位长度内较多的障碍物能够带来更多的扰动，促使增加火焰面积，提高放热强度，当障碍物间距（足够大）能够使前导激波产生马赫反射，则使DDT 更快地触发。Heidari 等[19-20]和Chen 等[21]通过数值模拟研究表明，不同构型的障碍物、化学反应机理、湍流模型、网格尺度等对火焰加速的具体过程有较大影响。总的来看，障碍物把管道沿流向变成周期性“收缩-扩张”的通道，在障碍物后缘产生展向漩涡结构，并在障碍物之间形成回流区。在火焰加速的初始阶段，展向旋涡与火焰的相互作用使火焰面折叠、拉伸增加火焰面积提高放热强度使火焰加速；在火焰加速的后期阶段，激波与障碍物的相互作用产生“热点”来触发爆震。

Oran 等[22]指出单步总包化学反应机理能够成功模拟爆震波的触发及传播，但并不能准确的对火焰加速的整个过程进行预测。其原因是单步反应模型的点火延迟时间通常小于多步反应[23]，而点火延迟时间是正确预测火焰加速及DDT 过程的1 个重要因素。因此，有必要采用详细化学反应机理对障碍物管道内的火焰加速及DDT过程进行研究。

本文通过求解耦合详细化学反应的2 维N-S方程，采用7 组分8 反应的氢气详细化学反应机理，对2种障碍物分布（对称型和交错型）的简化管道模型进行非定常模拟。

1 数值方法

1.1 计算模型与初始边界条件

计算模型详细参数如图1 所示。采用对称型和交错型2 种障碍物布置方式，分别记作Case1 和Case2。计算域总长度L=1000 mm，半宽H=20 mm，共设置了16组等间距分布的（间距为37.5 mm）障碍物，Case2 的第1 组障碍物与计算域左侧（封闭端）的距离为17.5 mm，最后1 组障碍物与计算域出口的距离为380 mm。2 种布置的障碍物尺寸相同，宽度w=2.5 mm，高h=10 mm。本文中采用氢气与空气的化学恰当比混合物，并近似认为空气由体积分数为21%的氧气与79%的氮气组成：H2+0.5（O2+3.61N2）→H2O+1.805N2。从图中可见，在计算域封闭端中央设置1个半径为2.5 mm、温度为2500 K 的高温区，用来模拟火花塞点火后的初始层流火焰。整个计算域初始压力为101325 Pa，初始温度为300 K。壁面为无滑移边界条件，出口为压力出口条件。

图1 计算模型详细参数

1.2 数值方法与验证

湍流在火焰传播过程中具有重要作用，尤其是在障碍物附近的剪切层中，高强度的湍流造成的小尺度旋涡对火焰面精细结构有重要影响。若要准确的模拟火焰传播过程中的流动机制及火焰传播机理，需要在多尺度上对3 维火焰面进行解析求解。但其计算量是不可接受的。对于障碍物管道内DDT 过程的模拟，Goodwin 等[24]认为2 维、3 维模拟得到火焰加速机制与DDT 过程图像是一致的，其主要差别为火焰加速过程中具有不同的时间尺度，具体表现为火焰传播速度、DDT时间不同。

本文忽略湍流的影响，采用基于压力的分离求解器Ansys Fluent V15 求解2 维非定常纳维-斯托克斯（Navier-Stokes，N-S）方程。压力修正方程采用压力隐式算子分裂（Pressure Implicit with Split-Operator，PISO）方法求解，在爆震燃烧计算中有成功的运用[25]。其中，时间离散采用隐式双时间步，空间对流项与扩散项分别采用2阶迎风格式与2阶中心差分格式。化学反应模型采用层流有限速率模型，为考虑化学动力学因素的影响，采用了7 组分8 步基元反应的氢气/空气多步化学反应机理[26]，此机理被认为能够较好地满足爆震波结构及其动力学特性的模拟。

层流火焰厚度δl为[27]

式中：λu、ρu、cp、Tu分别为未燃混气的导热系数、密度、定压比热与温度；Tb为已燃混气温度；Sl为层流火焰速度，文献[28]中的数据为2.1 m/s。

根据式（1）计算得到火焰厚度δl= 0.32 mm。考虑到计算量的问题，计算网格采用Δ=0.25 mm的均匀网格，能够对层流火焰面的基本形态进行捕捉。若要刻画精细的火焰面，则推荐文献[18，24]中的取值Δ =δl/60 ～δl/80。为捕捉爆震波的传播，同时考虑计算的稳定性，设置库朗（Courant-Friedrichs-Lewy，CFL）数小于0.2，则取时间步长为2×10-8s。

计算爆震波行进问题，用于验证化学反应模型的可靠性及网格、时间步长设置的合理性。通过CHEMKIN软件计算得到氢气/空气在101325 Pa、300 K 的条件下，理想的查普曼-朱格特（Chapman-Jouguet，CJ）爆震波的速度Vcj=1978 m/s，CJ 压力Pcj=1.56 MPa。不同时刻1维爆震波压力历程如图2所示。计算得到的爆震波速度为1977 m/s，爆震波压力为1.53 MPa，与理论值的相对误差为1.9%与0.5%，可以说明本文选取的化学反应模型及计算设置能够较为准确地描述爆震波传播的整体特性。

图2 1维爆震波压力历程

2 计算结果与讨论

2.1 障碍物管道内火焰加速过程

采用对称型障碍物（Case1）时，初始火焰传播过程如图3所示。黑色等值线表示O2的质量分数，代表火焰面位置，背景为涡量场云图。从图中可见，在第7.2 ms时（图3（a）），初始的层流火焰向下游传播越过第1 个障碍物，燃烧膨胀诱导了下游未燃混气的流动，在下游每组障碍物前缘的顶部脱落形成了层流涡，这是由障碍物顶部的剪切流动（Kelvin-Helmholtz不稳定性）造成的。随着时间的发展（t=8.4～9.2 ms），燃烧速率的增大推动下游未燃混气的流动加速，增加了障碍物顶端剪切层的不稳定性，层流涡逐渐脱落成旋涡串，充满障碍物之间的区域，形成较大回流区。

图3 初始火焰传播过程（对称型）

第4、5 组障碍物之间旋涡的发展过程以及火焰与旋涡的相互作用如图4 所示，其中绿色等值线表示火焰面。在t=5 ms 与t=7.2 ms 时（图4（a）、4（b）），从障碍物前缘脱落的旋涡进入障碍物之间的区域，并随着时间的推移回流区逐渐变大。随着时间的推移，旋涡的大小、强度增加并卷入回流区中心，使回流区的变形加剧。当火焰面前锋通过第4 组障碍物（图4（e））时，障碍物之间的回流区范围继续增加，旋涡串开始变形、合并。火焰面传播接近至第5 组障碍物时，相对光滑的火焰面边缘被卷入到回流区内部，发生变形、褶皱。增加的火焰面积提高了单位面积的放热强度，使燃烧更加剧烈，从而推动了火焰的加速传播。随着时间的发展，至t=9.365 ms 时（图4（g）），局部燃烧破坏了相对规则回流区，形成几组小的旋涡，当燃烧产物几乎充满整个障碍物之间区域时，回流区也随之破坏消失。

火焰传播中期在障碍物附近产生的局部爆炸如图5 所示，左列为温度，右列为压力梯度，时间范围t=9.796～9.818 ms。此时火焰传播速度接近声速，处于雍塞状态。从图5（a）中可见，在t=9.796 ms 时，在火焰面前形成了前导激波，与火焰面的距离约为20 mm。弓形的前导激波在上下壁面反射形成横向运动的反射激波，扫过火焰面前的未燃混气提高其热力学参数，能够提高化学反应的放热强度。同时，卷入上部障碍物之间回流区中的火焰面前缘存在很强的激波，这道激波在障碍物壁面反射增强，压缩脱落在已燃混气中的未燃混气团产生“热点”（图5（b）），热点发生爆炸在当地形成局部的高温高压，可以认为触发了爆震波。但爆震波向四周传播时由于缺乏足够的新鲜混气，迅速衰减成声学激波，向下游运动的激波与障碍物作用会形成新的爆炸中心，并推动火焰面加速，使其与前导激波更加接近。从图5（d）中可见，火焰面与前导激波非常接近，前导激波在障碍物的拐角处（壁面与障碍物面）分别产生了反射，压缩当地的新鲜混气，为燃烧向爆震波的转捩创造了条件。

图5 火焰传播中期在障碍物附近产生的局部爆炸

爆震触发及爆震波传播的过程如图6 所示。在t=9.875 ms 时，位于0.45 m 下部障碍物附近的未燃混气发生了局部爆炸，产生了爆震波核；爆震波核迅速发展绕射通过障碍物，往下游传播并与前导激波耦合。随后爆震波前锋迅速扩展到整个空间，同时与下壁面碰撞反射，在已燃混气中形成横波（图6（c）），完成了缓燃向爆震的转捩。从图6（d）～6（f）中可见爆震波触发后稳定传播的过程。在障碍物区域，爆震波呈周期性运动：先绕射通过障碍物爆震波前缘呈弓形，在障碍物之间的壁面形成马赫反射，横波呈横向周期运动。在此过程中，爆震波的平均传播速度基本保持不变。当进入光滑管道区域，爆震波锋面趋于平面。

图6 爆震触发及爆震波传播的过程

爆震管内火焰锋面沿程加速的历程如图7 所示。火焰逐渐加速至燃烧产物的声速，接近1100 m/s。这时的火焰状态称为快速火焰，是爆震转换前的最终状态。如果当地的几何条件并不能满足爆震触发的条件（如没有障碍物），则火焰会以恒定的速度传播。从图中可见，在x=500 mm 附近火焰传播速度有一个极高的突跃，火焰传播速度瞬间提高到3000 m/s，产生了“过驱爆震”，压力可达12 MPa 左右。但“过驱爆震”是不稳定的，会逐渐衰减成稳定的CJ 爆震，火焰传播速度稳定在1970 m/s左右。

图7 爆震管内火焰锋面沿程加速的历程

2.2 交错障碍物中的火焰传播过程

在交错障碍物下（Case2）的初始火焰传播及涡量场如图8 所示。从图中可见，障碍物后旋涡脱落的大小随着时间的推移而增长，其变化与Case1 的旋涡脱落类似。由于障碍物交错布置，相当于使火焰沿程传播的距离增加，火焰前锋呈“S”型运动。这使得初始的火焰传播速度慢于采用对称型障碍物时的初始火焰传播速度。

图8 交错障碍物下（Case2）的初始火焰传播及涡量场

爆震波的触发（交错型）如图9 所示，在x=0.45 m的障碍物上游产生了热点。值得注意的是，这里热点的产生并没有发生在障碍物附近，而是在壁面。局部爆炸产生的激波在障碍物处发生反射，同时绕射通过障碍物。

图9 爆震波的触发（交错型）

2.3 火焰传播速度的对比

Case1、Case2中心线附近初始火焰面位置随时间的变化如图10 所示。Case1 的火焰传播速度要略高于Case2的。当通过第5组障碍物时，Case1的火焰传播速度约为200 m/s（t=9.20～9.24 ms），Case2的火焰传播速度约为180 m/s（t=9.94～10.4 ms）。

图10 Case1、Case2中心线附近初始火焰面位置随时间的变化

DDT 时间定义为可燃混气被点燃到爆震波触发这段时间，用于表征着火-层流火焰-湍流火焰-爆震触发整个阶段的快慢程度。根据定义计算得到Case1与Case2的DDT时间分别为10.67与11.49 ms，二者触发的位置都在450±30 mm之间，说明采用对称型障碍物的火焰传播速度相对交错型障碍物的稍快。分析认为，对称型障碍物的堵塞比是交错型障碍物的2倍，高堵塞比造成的强扰动对燃烧的影响更大，从而使火焰更快传播。但纯2 维模拟忽略了火焰加速过程中的真实湍流效应，与考虑湍流的3 维模拟相比，2维模拟得到的DDT过程相对较慢[29]。并且，真实爆震波具有3维、不稳定的阵面结构，其传播模式也与2维情况下有较大区别[30]。

3 结论

（1）在火焰发展的初始阶段，障碍物、旋涡与火焰的相互作用是火焰加速的主要因素。

（2）在火焰发展的中期阶段，激波与障碍物的碰撞、反射在可燃混气产生爆炸中心，完成缓燃向爆震的转捩。

（3）由于堵塞比的差异，对称障碍物得到的初始火焰传播速度比采用交错型障碍物时的更快。