近临界马赫数条件下楔形诱导斜爆轰数值模拟

王成，吴京昌, 韩文虎, 杨同会

(北京理工大学 爆炸科学与技术重点实验室，北京 100081)

爆轰波是由激波和波后化学反应放热相耦合而成并以超声速传播的燃烧波. 由于爆轰波在高超声速推进系统中的巨大应用潜力，其相关研究得到广泛重视. 以爆轰波为基础的推进系统存在多种类别，如旋转爆轰发动机[1]及采用驻定斜爆轰波的冲压加速器、斜爆轰发动机[2]. 以斜爆轰推进系统为基础的相关研制取得了一定进步，斜爆轰发动机适应高马赫数飞行，且结构简单，由于斜爆轰反应区极短，因此燃烧室无须很长，可大大减少发动机的质量以及壁面损失. 如驻定斜爆轰波稳定，则系统无需采用火焰助稳. 然而斜爆轰波的不稳定性，使得形成稳定斜爆轰波较为困难. 因此，从内在爆轰物理机制方面对斜爆轰波系结构及不稳定性进行研究，对斜爆轰发动机的设计具有重要参考意义.

为探索ODW的基本特性及其在推进系统中的应用，已开展大量研究. Lehr等[3]在观测氢氧介质中高速弹丸飞行时，发现弹头部分产生斜激波诱导的燃烧现象，并同时观测到不同燃料当量比、弹丸结构及飞行马赫数等条件下，斜爆轰的结构不同. Li等[4]通过数值模拟发现斜爆轰波由斜激波、诱导区、斜爆轰波波面和一组爆燃波组成. 斜激波下方气体经过压缩与自燃，形成爆燃波. Viguier等[5]进行了以氢气为燃料的斜爆轰研究，实验所得斜爆轰结构稳定且实验重复性好，与Li的数值模拟结果基本一致. Palalexandris等[6]首次发现斜爆轰波波面上可能存在三波点及横波，并发现这种结构往往出现在较大化学反应活化能的可燃气体中. Choi等[7]通过高精度模拟，发现有三波点和横波构成的小尺度多维结构确实存在，且与正爆轰波的胞格结构明显不同. Silva等[8]通过氢氧混合基元反应模型，分析了不同楔面角、不同氢氧当量比及不同马赫数条件下斜爆轰波的发展演化，并对斜爆轰波诱导时间和诱导长度进行了研究. Liu和Wang等[9]从理论和数值模拟上发现Chapman-Jouguet 斜爆轰波(C-J ODW)在ODW结构中起重要作用，并解释了C-J ODW可以上行传播的原因. Teng等[10-12]对斜激波向斜爆轰波转变区域的爆燃波与激波耦合结构进行了研究，归明月等[13]则进一步解释了斜爆轰精细结构，并分别探讨了尖劈角度及来流马赫数对劈面压缩性的影响. 综上所述，关于斜爆轰的研究大多在无粘性反应流体范畴，且集中在高马赫数条件强过驱的驻定斜爆轰行为，对于近临界马赫数条件下粘性斜爆轰传播特征的研究仍然很少. 尽管Li等[14]、Fusina等[15]研究认为黏性对斜爆轰波的结构影响仅局限于薄边界层，对主流区结构几乎无影响，而近年来的研究表明，黏性扩散对爆轰波后未反应物质的燃烧及斜爆轰上行速度等方面存在影响，不可忽略. 因此开展有、无黏性扩散的斜爆轰数值模拟，揭示斜爆轰的阵面结构及其临界马赫数条件下的传播特征.

1 控制方程

控制方程采用具有黏性、热传导及分子扩散，带化学反应的2维多组分可压缩流体Navier-Stokes方程组，具体形式为

-AρYexp(-Ea/RpT)，

p=ρRpT/M.

式中：Y为可燃气体质量分数；h=QY+CpT=QY+RpγT/(γ-1)为焓；e=QY+CVT=QY+RpT/(γ-1)为内能；Q为反应热；CV和Cp分别表示定容和定压比热容；M为可燃气体摩尔质量；Rp为理想气体常数；γ为可燃气体绝热指数；u、v分别为流体速度在x方向和y方向的分量；ρ、p、T、Y分别为密度、压力、温度、反应物的质量分数.

能量扩散矢量σj和应力张量ζij为

式中：μ=pv为动黏性系数；Pr和Sc分别为Prandtl和Schmidt数；v为分子运动黏度；δij为Kronecker符号.

上述物理量及均质气体混合物的初始状态为量纲一的参数，因此初始密度、压力和温度是统一的. 根据未燃气体的状态，使方程量纲归一化，

表1 量纲一化参数[16]

2 数值计算方法及收敛性

2.1 数值计算方法

在楔形斜面上设定一个量纲一化的大小800×600矩形域，由于气流以超声速传播，会在楔面上方形成一道以顶点为起点的斜激波，斜激波诱导可燃混合气体发生燃烧反应，在下游形成复杂的斜爆轰结构. 为方便数值模拟，将坐标系旋转使之与楔面方向一致(如图1). 通过改变入流马赫数来模拟固定楔形角为30°的斜爆轰. 空间离散采用基于WENO保守有限差分格式的5阶局部特征方法，时间离散采用3阶TVD Runge-Kutta时间离散. 并用数值算例验证了相对网格分辨率和保正数值格式的正确性[17].

2.2 收敛性验证

采用20，40及80 pts/L1/2(L1/2为半反应区厚度)3种网格精度对Ma=12,楔面角30°的2维无黏斜爆轰进行收敛性验证. 由图2可知，当网格精度增加到80 pts/L1/2,发现爆轰结构变化很小，表明40 pts/L1/2的模拟结果收敛于80 pts/L1/2时的结果. 因此，40 pts/L1/2网格精度能够捕捉斜爆轰结构特征及爆轰角. Choi与Kim等[18]采用3阶MUSCL-TVD方案对对流通量进行了离散化，网格精度为40 pts/L1/2能够满足斜爆轰模拟要求. 由于目前的数值模拟采用的格式为5阶WENO计算格式，占用较少网格就能达到较高精度，且经过验证其精度可以达到5阶3维爆轰[19]. 高网格精度下接触面会产生强涡，Choi及Kim等[18]进行了详细的讨论. 本模拟对于高马赫数条件下的斜爆轰进行模拟，高过驱度会抑制接触面上的不稳定特征，因此小尺度涡不会出现[20].

3 数值结果分析和讨论

3.1 马赫数对斜爆轰结构的影响

根据爆轰极理论推导得出，图3为楔面角30°时不同马赫数下的斜爆轰角与来流马赫数间的关系(Ma=7,8,9.5,12),对比理论解给出临界马赫数为8.5. 观察图像可知，数值模拟结果预测的临界马赫

数与理论值吻合较好. 随着马赫数改变，斜爆轰呈现出驻定、非驻两种不同的斜爆轰模态. 当来流马赫数高于临界值时，斜爆轰在下游起爆，上行一段距离后驻定；来流马赫数低于临界值时，斜爆轰波不断上行，最终脱体. 当马赫数大于或接近临界值8.5时，斜爆轰从下游起爆后上行至上游，最终驻定.

图4给出马赫数为7时，t=108.8，264.1，417.9时的流场温度分布与压力分布. 可以看出t=108.8时，斜爆轰处于起爆的初始阶段，斜激波与斜爆轰波的过渡通过主三波点完成，同时爆轰下游存在C-J斜爆，推动斜爆轰波面位置向上游移动. 对比图4(a)中的温度、压力分布，爆轰压力与火焰阵面之间解耦，不能形成驻定的斜爆轰；随着时间的推进，斜激波起爆区的长度逐渐缩短，斜激波向斜爆轰过渡的三波点向上游推进；在t=264.1时，可清晰地看到主三波点连接一激波，该激波经楔面反射，并透过滑移间断诱发下游流场失稳；在t=417.9时，斜爆轰脱体，没有形成驻定爆轰波. 主三波点的纵坐标最终接近x轴，且跑出计算域外，斜爆轰波后大面积燃烧反应进行，且滑移间断面部分出现一系列明显的涡结构.

图5给出了t为60.4，111.1，164.4时，来流马赫数为8时的流场温度分布与压力分布数值模拟结果. 观察模拟结果可知t=60.4时，起爆的斜激波区域长度较图4(a)短，斜激波向斜爆轰波的过渡较为平缓. 而主三波点附近的结构较为复杂，斜爆轰波阵面下游结构及阵面后横波传播均呈水平直线. 从流场整体温度均值及压强均值来看，t=60.4时的爆轰温度及爆轰压力还较低；t=111.1时，主三波点附近的结构与图4时相近时刻的爆轰场相比较而言，已基本稳定，观察图5(b)中的温度云图可以发现，爆轰波阵面后方的反应十分剧烈，并在爆轰右上部波阵面后方出现一些未反应物质. 与马赫数为7的模拟结果比较，斜爆轰角较小，且滑移线上涡结构更为复杂；t=164.4时，斜爆轰仍没有驻定，主三波点继续沿x轴上行，滑移线下游的涡结构发展充分，涡呈水平排列且尺度均匀，此时的未反应物主要集中在波阵面后右上端. 爆轰最终脱体，没有实现驻定，上行至计算域外.

3.2 与有粘状态下的流场结构比较

图6给出了Ma=9.5的无黏斜爆轰在t为81.9，178.5，274.7及324.9时的温度场模拟结果. 可以看出t=81.9时，斜爆轰起爆的上游斜激波区域长度较图4(a)和5(a)更短，斜爆轰角较大；在x=100～150之间区域，明显观察到三波点后连接一激波，在发生反射，并透过滑移间断面与爆轰波阵面产生的横波相互作用，使下游流场受到扰动；爆轰波产生的横波遇楔面后也发生较明显反射，使得接触间断面下游出现一系列涡结构.t=178.5时，涡线下方的涡结构进一步发展，斜爆轰角进一步增大，主三波点上行速度渐缓；t=274.7时，斜爆轰的波面后方出现较明显的胞格结构，三波点及滑移线的纵坐标水平变得很低，斜爆轰上行过程中出现全局波动现象；主三波点上行至x=35处，而后在t=324.9时退至x=50处；在t=324.9时，主三波点诱发的激波，及滑移线上的涡结构变得更加明显，爆轰波阵面后右上端紧贴波面处保留着部分未反应物质，此时斜爆轰已实现驻定，不再继续上行传播.

图7给出了黏性斜爆轰在t=32.9，78.6，134.9，222.2，274.7及324.9时的温度场云图. 与图4(a)及图5(a)比较发现，在t=32.9时，斜爆轰角更小，诱导区的长度更短，在t=78.6时，x=100～150之间涡线下方区域的激波及反射结构与图6(a)相似,但其爆轰波阵面后反应区内的小尺度结构更加规则和均匀. 比较图6(d)、(e)、(f)中的主三波点横坐标位置可知，黏性斜爆轰形成过程中也出现上行波的波动. 通过对比有、无粘性模拟结果，可以看出黏性对过驱情况的影响更为显著，而对于接近C-J状态的情况影响较小. 对于接近C-J状态的斜爆轰，黏性C-J爆轰上行传播速度要比非黏性情况快. 整体上，黏性不影响斜爆轰全局特征，但黏性抑制了滑移线上不稳定诱发的小尺度涡结构的产生.

4 结 论

基于一步总包化学反应模型，采用NS及Euler控制方程组对2维楔形斜爆轰进行了数值模拟，得到主要结论如下.

① 随着马赫数值增大，斜爆轰呈现出驻定、非驻斜爆轰模态. 临界马赫数约为8.5，来流马赫数高于临界值时，斜爆轰可实现驻定；来流马赫数低于临界值时，斜爆轰波不断上行，成为脱体爆轰.

② 通过对比有无黏性两种情况下的ODW模拟，可以看出黏性对过驱情况的影响更为显著. 整体上，黏性不影响整体特征，但黏度效应抑制了斜爆轰结构迹线上因不稳定引起的小规模旋涡. 对于接近C-J状态的ODW，爆轰出现上行波动现象，爆轰在粘性情况下的上行传播速度要快于非粘性情况.