1 引言

斜爆轰是超声速的可燃混合气体流过一定角度的楔面，经过斜激波压缩、诱导燃烧直至迅速放热的一种特殊爆轰现象。由于斜爆轰波可以驻定，燃烧时间极短且拥有很高的燃烧效率，因此在超燃冲压推进系统和冲压加速器等方面都有良好的应用前景^[1]。在驻定斜爆轰发动机（Oblique detonation wave engine，ODWE）中，其来流速度大于CJ速度，推力连续，且相比常规超燃冲压发动机来说，燃烧室更短，飞行马赫数范围更宽，能在更高的飞行马赫数下（Ma > 10）保持较高的燃烧效率，熵增更小，另外还具有重量轻、飞行阻力小及易于重复启动等优点^{[2, 3]}。

目前对ODWE的研究主要局限于理论分析和机理研究，国内外对斜爆轰进行了少量原理性试验^{[3, 4]}，并在数值计算方面进行了大量工作，探究了来流马赫数、斜劈角等参数对斜爆轰波的影响^[5~16]。Li等^[5]对斜爆轰波进行了数值模拟研究，发现斜爆轰波结构由未反应斜激波、诱导燃烧区域、一系列爆燃波和斜爆轰波阵面组成。Silva等^[6]对斜爆轰波进行了参数化研究，探究了来流马赫数、温度、压力、斜劈角对斜激波转变斜爆轰波形态的影响。归明月等^[7]进一步揭示了斜爆轰波的精细结构，并分别探讨了来流马赫数和尖劈角度对劈面压缩性的影响。Liu等^[8]研究了来流马赫数对斜爆轰波结构的影响，发现随着来流马赫数增加，CJ斜爆轰波反射呈现出四种不同的形态。吴伟等^[9]则对不同锥角下的斜爆轰波进行了模拟研究，发现斜爆轰波只有在一定角度范围才能形成，锥角较小只存在斜激波，过大锥角将导致爆轰波脱体，但是并未探究不同来流马赫数下斜爆轰波形成的条件。王爱峰等^{[10, 11]}，滕宏辉等^[12~16]利用欧拉方程和单步化学反应模型，对斜爆轰波进行了数值模拟，重点研究了斜爆轰波面上多维复杂结构的形成机理，包括斜激波诱导爆轰、爆轰波后的未反应气团的出现、接触间断面的扰动和胞格结构的不规则性等，但是并没有探究来流马赫数和楔面角对其共同影响。

已有的研究主要局限于楔面角和来流马赫数对驻定斜爆轰波形态结构的单一影响和定性分析上，没有探究斜爆轰波在不同楔面角和来流马赫数下的起爆条件以及楔面角和来流马赫数对斜激波诱导区域长度、接触间断面、未反应气团等形态结构的共同影响。本文采用数值模拟方法分析了斜爆轰波的起爆条件，探讨来流马赫数、楔面角对驻定斜爆轰波形态的共同影响，并揭示其作用机理。

2 数学方法 2.1 控制方程

忽略粘性和热传导，驻定斜爆轰波的控制方程采用带化学反应的二维欧拉方程。对于化学反应，采用Arrhenius动力学控制的单步不可逆化学反应模型^[17]。守恒形式下，无量纲控制方程可以写为

$\frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial y}} = \mathit{\boldsymbol{S}}$

(1)

式中守恒变量矢量U、通量矢量F和G、源项S分别为

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{l}} \rho \\ {\rho {v_x}}\\ {\rho {v_y}}\\ {\begin{array}{*{20}{l}} E\\ {\rho Y} \end{array}} \end{array}} \right], \;\mathit{\boldsymbol{F}} = \left[ {\begin{array}{*{20}{l}} {\rho {v_x}}\\ {\rho v_x^2 + p}\\ {\rho {v_x}{v_y}}\\ {\left( {E + p} \right){v_x}}\\ {\rho {v_x}Y} \end{array}} \right], }\\ {\mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{l}} {\rho {v_y}}\\ {\rho {v_y}{v_x}}\\ {\rho v_y^2 + p}\\ {\left( {E + p} \right){v_y}}\\ {\rho {v_y}Y} \end{array}} \right], \;\mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{l}} 0\\ 0\\ 0\\ {\begin{array}{*{20}{l}} 0\\ \omega \end{array}} \end{array}} \right]} \end{array} $

(2)

式中v_x，v_y分别是流体速度在x和y方向的分量，ρ是密度，p是压强，E是单位体积的总能量，Y是反应物质量分数，并基于来流未燃气体的状态量无量纲化^[18]

$ \begin{array}{l} \bar \rho = \frac{\rho }{{{\rho _0}}}, \;\;\bar p = \frac{p}{{{p_0}}}, \;\;\bar T = \frac{T}{{{T_0}}}, \;\;Ma = \bar V = \frac{V}{{{c_0}}}, \\ \bar L = \frac{L}{{{L_0}}}, \;\;\bar E = \frac{E}{{{p_0}}}, \;\;\bar K = \frac{{K{L_0}}}{{{c_0}}}, \;\;\bar q = \frac{q}{{c_0^2}}, \;\;\overline {{T_{\rm{i}}}} = \frac{{{T_{\rm{i}}}}}{{{T_0}}} \end{array} $

式中p₀=101.325kPa，T₀=300K，$ {c_0} = \sqrt {\gamma R{T_0}} $，特征长度L₀为反应区长度。总能量E定义为

$ E = \frac{p}{{\gamma - 1}} + \frac{1}{2}\rho \left( {v_x^2 + v_y^2} \right) + \rho qY $

(3)

式中q是热释率，γ是比热比。源项ω采用阿累尼乌斯的形式计算

$ \omega = - K\rho Y{{\rm{e}}^{ - \left( {{T_{\rm{i}}}/T} \right)}} $

(4)

式中T是温度，T_i是活化温度，K是恒定指前因子。对于理想气体，状态方程为

$ p = \rho T $

(5)

2.2 数值方法

本文控制方程采用有限元的方法进行离散，采用均匀的空间网格，时间被离散到初始时间t₀到终止时间t_f之间的N个CFL自适应时间步；激波捕捉采用二阶的Godunov格式^[19]。Greenough等^[20]总结出在相同的计算花费上，二阶的Godunov格式比五阶WENO格式具有更高的精度。尽管很明显，更高阶的格式对于一些应用来说有优势，但是本文的算法精度还是在空间和时间上都限制到二阶精度。

初始时刻，绝热条件下，γ=1.2，一定来流马赫数Ma的超声速可燃混合气流向楔形物体运动，楔面与水平面夹角为θ角，如图 1所示。为了模拟方便，把坐标系逆时针旋转θ角，使之与楔面方向一致，x轴和y轴分别和计算区域边界平行。

初始来流中，ρ₀=1，p₀=1，T₀=1，u₀=Macosθ，v₀=-Masinθ，反应物质量分数Y₀=1。为使半反应区中有足够多的网格点，无量纲计算区域设为0.0≤x≤50.0，0.0≤y≤8.0，x₀=1.0，采用均匀的1200×300网格计算。热释率q=50，活化温度T_i=50，恒定的指前因子K取2566.42，文中所采用的变量均为无量纲。

边界条件如图 1所示，入口和上边界（深色部分）采用入口边界条件，入流条件给定为初始值ρ₀, u₀, v₀, T₀, p₀和Y₀；出口（浅色部分）采用出口边界条件；下边界，x≤x₀的区域（浅色）采用出口边界条件，x≥x₀的平板区域采用滑移反射边界条件。

3 数值验证 3.1 斜爆轰波角验证

假定平面斜爆轰：（1）来流为预混可燃气体，混合均匀，流动稳定；（2）化学反应层等效于放热，放热量为Q；（3）流动无粘绝热。根据上面的假设条件，利用流体守恒方程，可以得到下面关于放热量、楔角、马赫数和斜爆轰波角之间的关系^{[10, 13]}

$ \begin{array}{*{20}{l}} {\frac{{{\rm{tan}}\left( {\beta - \theta } \right)}}{{{\rm{tan}}\beta }} = }\\ {\frac{{1 + \gamma Ma_n^2 \pm \sqrt {{{(Ma_n^2 - 1)}^2} - 2\left( {\gamma + 1} \right)Ma_n^2q} }}{{\left( {\gamma + 1} \right)Ma_n^2}}, }\\ {Ma_n^2 = M{a^2}{\rm{si}}{{\rm{n}}^2}\beta } \end{array} $

(6)

式中Ma_n为波前法向马赫数，β为斜爆轰波角。

不同来流马赫数和楔面角情况下，斜爆轰波角的理论值和数值结果及相对误差见表 1。由表可见，数值解与理论解非常接近，计算的最大相对误差为1.25%，说明数值解的结果是可信的。

3.2 斜爆轰波流场验证

为了进一步验证驻定斜爆轰波结构中复杂的结构，将数值结果与实验结果进行对比，图 2是实验纹影图^[21]与数值温度云图的对比结果。图中斜激波、斜爆轰波与爆燃波交汇成处可以看到典型的三波点结构以及三波点后向下游发展的接触间断面。三波点后产生的激波传播到楔面发生反射，透过接触间断面，与爆轰波阵面产生的横波相互作用扰动下游流场，使得胞格呈现出不规则的结构。同时爆轰波阵面产生的横波传播到楔面后也发生反射，导致接触间断面下游出现一系列的涡结构^[10]。因此，数值模拟的结果和实验吻合较好，所采用的数值计算能够对斜爆轰波结构进行详细的解析。

4 结果与讨论

驻定斜爆轰波结构受到很多因素的影响，其中来流马赫数和楔面角是两个主要的因素。探究来流马赫数与楔面角对斜爆轰波的影响有很大的研究价值，尤其是对于驻定斜爆轰波发动机的研究，不同的斜爆轰波结构也相应要求驻定斜爆轰波发动机应有不同的构造。

4.1 驻定斜爆轰波的起爆条件

在一定情况下，超声速可燃预混气流遇到楔面，形成不了驻定斜爆轰波。通过研究不同来流马赫数和不同楔面角对驻定斜爆轰波的影响，发现斜爆轰波只有在一定的来流马赫数和楔面角范围内才能形成。较小的来流马赫数和楔面角情况下只存在斜激波，如图 3所示。较大的来流马赫数和楔面角情况下则会导致脱体爆轰的生成，如图 4所示，脱体爆轰波在楔面之前就会形成，爆轰波阵面为弧形。因为较小的来流马赫数和楔面角，y方向速度较小，所产生的斜激波较弱，不易产生燃烧；而较大的来流马赫数和楔面角会导致弓形脱体激波，进而引发脱体爆轰。

进一步的研究发现，在给定的q，T_i和K情况下，Ma=9.1，30° < θ < 40°才能形成斜爆轰波，而Ma=13.7，20° < θ < 50°就能形成斜爆轰波。来流马赫数越小，形成斜爆轰波需要的楔面角范围则越窄。因为较小的来流马赫数，y方向速度较小，所形成的斜激波较弱，不易产生燃烧，因此形成爆轰需要较大的楔面角；而较小的来流马赫数，最大可能折转角θ_max也较小，楔面角超过θ_max就会形成脱体激波^[22]，进而会生成脱体爆轰。

4.2 来流马赫数与楔面角对驻定斜爆轰波后形态结构的影响

选取不同的来流马赫数和楔面角，模拟得到的驻定斜爆轰波温度云图如图 5（a）~（f）所示，无量纲温度区间为1~21，颜色越深的地方表示温度越高。

从上述几幅图中，可以发现，来流马赫数越高、楔面角越大，斜激波诱导区域长度越短，斜爆轰波阵面越平滑，波后混合气体温度越高，未反应气团越少，胞格结构越规则，接触间断面越靠近楔面，其下游出现的涡结构也越少。因为，来流马赫数越高、楔面角越大，y方向速度越大，通过楔面反射，压强也越高，通过压缩获得的温度也越高，所产生的斜激波诱导区域长度则越短，三波点后产生的激波对波后流场的影响也越小，所以胞格结构越规则，接触点断面越低，与经过楔面反射的横波作用产生的涡结构也越少；同时，由于压缩后的混合气体温度升高，化学反应速率加快，即混合气体燃烧更剧烈，所以斜爆轰波阵面更平滑，燃烧产物温度更高，未反应气团减少。

量取不同来流马赫数、楔面角下驻定斜爆轰波中斜激波诱导区域的长度，定量研究斜激波诱导区域长度与楔面角和来流马赫数的关系，如图 6和图 7所示。图 6为来流马赫数Ma=13.7时，楔面角与无量纲斜激波诱导区域长度的关系折线图；图 7为楔面角θ=30°时，来流马赫数与无量纲斜激波诱导区域长度的关系折线图。

从图 6，图 7中可以发现，随着来流马赫数和楔面角的增大，斜激波诱导区域的长度减小，且来流马赫数、楔面角越大，斜激波诱导区域的长度减小的速率越小，斜激波诱导长度和来流马赫数、楔面角近似满足幂函数关系。

5 结论

本文对驻定斜爆轰波的起爆条件和形态结构进行了数值模拟研究，得到如下主要结论：

（1）通过模拟得到的驻定斜爆轰波角与理论分析值进行对比，误差较小，且数值模拟结果与实验纹影图吻合较好，验证了数值模拟方法的有效性。

（2）来流马赫数为9.1时，楔面角处于30°~40°才能形成驻定斜爆轰波，而来流马赫数为13.7时，楔面角处于20°~50°就能形成斜爆轰波。较小的来流马赫数和楔面角情况下，只存在斜激波；较大的来流马赫数和楔面角情况下，会导致脱体爆轰的生成。来流马赫数越小，斜爆轰波的起爆条件越苛刻。

（3）随着来流马赫数和楔面角增加，波后混合气体温度升高，爆轰反应变得更加剧烈，使得斜激波诱导区域变短，斜爆轰波阵面更加平滑，未反应气团减少，胞格结构更规则，接触间断面更靠近楔面，其下游出现的涡结构也减少。斜激波诱导长度和来流马赫数、楔面角近似满足幂函数关系。