0   引　言

实现宽工作范围内超燃冲压发动机的更优性能一直是研究者们追求的目标，双模态超燃冲压发动机的提出就是为了实现这一目标。研究者可通过对流道和燃烧释热进行设计使双模态发动机根据不同工作条件调整燃烧模态，从而获得更好的工作性能。对于宽范围运行的双模态发动机而言，实际飞行时会存在一个连续变化的动态工作过程，例如在飞行器加/减速时会伴随马赫数（Ma）、迎角的变化。来流条件的改变势必会带来发动机内速度、压力等流动参数以及激波、边界层分离等流动现象的明显变化。同时，为适应飞行条件与推力需求的变化，发动机也需不断调整工作模式。在基于超燃冲压发动机的HIFiRE项目飞行试验中，在86.2 kPa的动压下，飞行器的飞行马赫数5.5加速至8.5，双模态发动机完成了从亚燃模态到超燃模态的转换^[1]。可见变来流条件下发动机流动和燃烧特性以及工作模式发生变化是实际工作时必然会遇到的情形。

针对双模态发动机流动与燃烧特性，国内外学者开展了大量地面试验与数值模拟研究。Fotia^[2]通过直连式双模态燃烧室试验对燃烧模态转换开展了研究，发现系统所处的燃烧模态是一个与边界层状态、燃烧室中火焰模式、进气道以及注油条件等都相关的复杂函数。浮强等^[3]则针对来流总温对模态转换边界的影响开展了直连式试验研究，并指出总温升高会使发动机从超燃模态向亚燃模态转换时的油气比上升。Xiao等^[4]开展了一系列双模态燃烧室直连式试验，基于壁面压力建立了可以在试验中快速判断燃烧模态的定量判别准则；并通过风洞试验对比了燃烧室构型一致、来流参数匹配的自由射流与直连式发动机在燃烧模态、流动及燃烧上的差异^[5]。Goodwin等^[6]通过试验研究了双模态发动机凹槽里的乙烯空气预混燃烧情形。晏至辉等^[7]为了获得双模态燃烧室流场内的温度分布特性，开展了直连式试验并利用可调谐的相干反斯托克斯拉曼（CARS）技术进行了流场测量，发现主喷油位注入当量比为0.2的燃料对扩张段及后部的推力性能没有明显改变，但对燃烧室区域的温度和压力分布会产生较大影响。连欢等^[8]开展了发动机模态转换及推力突变的试验研究，指出模态转换及推力变化过程的本质是释热总量与内流道匹配，发动机动态飞行轨迹气动热和燃烧热积分效应可能改变边界层特性和发动机抗反压能力。车庆丰等^[9]搭建了基于超燃冲压发动机的10 kHz高频OH-PLIF试验系统，开展了燃烧室测量试验。测量技术的不断发展使得试验能提供更多有效信息，但对于发动机全流场的流动特性、结构及参数，数值模拟仍然是一种了获取更多流动细节的手段。

Yentsch等^[10-12]较为系统地针对HIFiRE发动机开展了数值研究工作，就模态转换时发动机内边界层、激波与膨胀波波系、凹槽动力学以及燃烧化学等特性开展了研究，对比了矩形与轴对称双模态发动机在模态转换时的差异。Kouchi等^[13]通过瞬态数值模拟，指出随着燃料当量比增加，燃烧室下游的燃烧诱导压力也随之增加；当压力超过某一临界值时，火焰前锋从燃烧室的后部迅速上传至喷注点附近，发动机进入强燃烧模态。Tian等^[14]针对不同注油当量比下的煤油燃料双模态发动机开展三维数值模拟，发现热力喉道的位置与分离区域的大小均受当量比的影响。林旭阳等^[15]运用FLUENT分析了燃烧室壁面换热与超声速燃烧之间的关系，模拟了壁温降低、激波串后移、发动机由亚燃逐渐向超燃转换的过程。

总的来说，尽管针对双模态发动机的研究已开展多年，逐渐走向工程化，但对发动机内复杂波系、边界层分离、燃烧释热特性等的理解尚不透彻。马赫数和注油量的改变是发动机在实际飞行中极有可能面临的情况，可能会带来发动机模态、流动和释热特性的改变。研究不同马赫数及当量比下发动机流动与燃烧的特性，分析流场内波系、分离以及燃烧等基础现象的演化过程，有助于进一步提升对发动机在工程应用中的控制能力。本文在直连式试验的基础上验证计算方法对该类矩形截面双模态发动机的模拟能力，通过一系列三维定常数值模拟，详细分析不同工况下双模态发动机的流动和燃烧特性。

1   计算方法及验证

1.1   计算方法

本文基于AHL3D软件对发动机开展三维定常数值模拟。AHL3D软件平台可以模拟二维或三维^[16]、定常或非定常、完全气体或化学非平衡流动^[17]，可以对激波边界层干扰^[18]、稳焰燃烧^[19]等复杂现象进行模拟^[20]，其对超燃冲压发动机的模拟能力得到了大量算例的验证。求解三维直角坐标系下的雷诺平均N–S方程，方程形式如下：

式中：${\boldsymbol{F}}、{\boldsymbol{G}}、{\boldsymbol{E}} $为无黏通量；$ {{\boldsymbol{F}}}_{v}、{{\boldsymbol{G}}}_{v}、{{\boldsymbol{E}}}_{v} $ 为黏性通量；${\boldsymbol{S}} $为源项；t为时间；$ {\boldsymbol{Q}}= {（\rho ，\rho u， \rho v，\rho w，\rho {E}_{t}，} \;{\rho {Y}_{i}）}^{\text{T}} $，其中u、v、w分别为x、y、z方向速度，$\; \rho 、{Y}_{{i}} $分别为气体密度和组分的质量分数，E_t为气体的总内能，$ {E_{{t}}}{\text{ = }} e + 0.5（{u^2} + {v^2} + {w^2}） $（e为热力学内能）。

湍流模型采用Kok^[21]提出的TNT k–ω模型。求解三维化学反应控制方程时采用隐式有限体积法离散，对无黏对流项进行离散时采用Van Leer提出的MUSCL方法，无黏通量为Steger−Warming分裂格式，黏性通量的计算采用Gauss定理构造方法，煤油化学动力学模型采用中国空气动力研究与发展中心（CARDC）^[22]简化的十二步十组分煤油简化模型。壁面采用无滑移边界条件，考虑为绝热壁。入口边界按来流参数设置条件，出口采用外推法。

1.2   模型及网格

本文针对图1所示的直连式发动机模型开展研究。模型总长约1130 mm，由等直隔离段、燃烧室及扩张段组成。隔离段长470 mm，入口尺寸为30 mm×150 mm。燃烧室由第一凹槽（l₁/h₁≈5）、凹槽间的过渡段（上壁面扩张2°）以及第二凹槽（l₂/h₂≈12）组成。第二凹槽后的扩张段上下壁面扩张2°。24个注油孔均布于隔离段尾部（第一凹槽前15 mm处），孔直径0.3 mm。

图  1  直连式发动机构型

Fig.  1  Direct connect scramjet configuration

下载: 全尺寸图片 幻灯片

以隔离段入口为零点（x=0），10个小孔沿展向（z方向）均布于第一凹槽前15 mm的隔离段尾部（x=0.455 m），垂直于流向（x方向）以声速喷注室温煤油。

由于所研究的发动机构型沿展向对称，为减少计算量，采用半模计算，总网格量约为755万，如图2所示。为能准确模拟边界层流动，壁面法向（y方向）第一层网格高度为0.001 mm，以确保无量纲壁面距离y⁺<1。计算过程收敛判断准则为：流场结构不再明显改变，能量最大残差下降3个量级，入口流量与喷油量之和与出口流量相差小于2%，继续计算2万步后流量差小于0.2%。

图  2  计算采用的网格

Fig.  2  Numerical simulation mesh topology

下载: 全尺寸图片 幻灯片

1.3   计算与试验的对比验证

针对上文介绍的发动机模型，课题组前期开展了直连式试验，如图3所示。试验模拟了飞行马赫数为6.0工况下隔离段的入口参数，来流为氢气和空气加热的高焓气体，确保氧气的体积比为21%。隔离段入口来流马赫数为3.0，总温为1500 K，总压为2.1 MPa，注油当量比为0.6。

图  3  直连式试验模型安装图

Fig.  3  Installation of direct-connected test model

下载: 全尺寸图片 幻灯片

为验证本文所采用计算方法的可靠性，采用前述方法及网格对试验状态进行三维模拟，将计算所得上壁面中心线静压（简称壁面压力）与试验所得压力进行对比，结果如图4所示，图中横坐标（长度，x）以隔离段口高度（h）无量纲化，纵坐标（壁面静压，p_wall）以隔离段入口静压（p_in）无量纲化。由图可见，计算与试验所得的壁面压力吻合良好，说明本文采用的计算方法对该模型具有较好的模拟能力。

图  4  试验与计算所得壁面压力对比

Fig.  4  Comparison between calculation and experiment of wall pressure

下载: 全尺寸图片 幻灯片

1.4   本文研究状态

为模拟实际飞行过程中发动机可能面临的马赫数与当量比同时改变的情况，本文针对来流条件与当量比同时改变的6个不同飞行工况对直连式发动机开展三维数值模拟研究。

采用设计点为Ma=6、飞行高度为24 km的前体进气道来计算获得符合真实条件的直连式入口参数，飞行弹道动压为73.831 kPa。对前体进气道进行三维数值模拟，计算收敛后得到喉道处一维质量平均参数并将其作为所研究的直连式发动机的入口参数。

在飞行马赫数由5.5增大至6.0的过程中，随马赫数增大逐步减小当量比（从1.0减小至0.5）。6个计算工况的命名与具体来流参数如表1所示。

表  1  计算状态及来流参数

Table  1  Research Cases and inflow conditions

工况	飞行马赫数	当量比	隔离段入口参数
			马赫数	静压/kPa	静温/K
Case 1	5.5	1.0	2.57	72.92	680.14
Case 2	5.6	0.9	2.61	72.21	688.28
Case 3	5.7	0.8	2.67	71.24	694.88
Case 4	5.8	0.7	2.73	69.91	699.11
Case 5	5.9	0.6	2.77	69.23	708.58
Case 6	6.0	0.5	2.87	68.01	702.35

下载: 导出CSV 

| 显示表格

2   结果与分析

飞行马赫数与当量比同时发生改变是发动机实际运行中极有可能遇见的情况。来流条件的改变势必会带来发动机内速度、压力等流动参数以及激波、边界层分离等流动现象发生明显变化。注入燃料并点火后，当量比的变化会带来不同的热释放，引起流动现象及流场结构的显著改变。

2.1   发动机内波系及流动特性

本小节探究了6个工况下发动机壁面压力及一维质量平均马赫数（Ma_mw）的沿程分布曲线及对称面上马赫数波系云图分布。如图5所示，不同工况下的壁面压力分布呈现出的形式类似，在出现明显的燃烧压升之前，隔离段前部（x=0～0.38 m）的压力提升相对较为平缓，无明显波动。这部分平缓压升主要由流道内不断发展增厚的边界层引起。

图  5  不同工况下壁面压力沿程分布

Fig.  5  Upper wall pressure distributions for different cases

下载: 全尺寸图片 幻灯片

本文进行数值模拟时对隔离段入口给予均匀来流，均匀入口条件与壁面无滑移假设相互作用会形成弱激波并在壁面间反射，弱波系也会带来一定的压升。在隔离段后部及燃烧室内，不同状态下的压升起点及压升幅度差异明显。具体来说，Case 1时，压升起点（x=0.38 m）最靠前，在注油位（x=0.455 m）之前75 mm处，随着飞行马赫数增大、当量比降低，压升起点逐步朝燃烧室方向后移，燃烧室及扩张段内的压升幅度也整体下降。Case 5、6时压升起点已固定在注油位处。比较所有计算工况，需要注意的是，Case 3、 4之间壁面压力存在明显的突变，尤其是在压升起点变化很大。这与图6中一维质量平均马赫数的分布一致，随着来流马赫数的增加与当量比的降低，燃烧室内一维质量平均马赫数整体不断升高，但Case 3、 4之间的变化最明显。从一维质量平均马赫数沿程分布规律上看，Case 1～3分布趋势较为相似，燃烧室内一维质量平均马赫数波动相对较小。与Case 1～3相比，Case 4～6燃烧室内马赫数波动幅度明显更大。这一定程度上反应了两类不同的流场形态。

图  6  不同工况下一维质量平均马赫数沿程分布

Fig.  6  One-dimensional mass average Mach number for different cases

下载: 全尺寸图片 幻灯片

从图7中不同工况与状态混用下流场对称面马赫数波系云图可以获得更直观的感受，由图可见，Case 1～3明显具有相似的流场形态，燃烧反压前传并诱导边界层分离，激波边界层相互作用使隔离段内形成了典型的斜激波串结构，激波串前缘位于注油位上游。同时，与隔离段内形成的较强激波串相比，燃烧室内的波系则相对较弱，这也是Case 1～3中燃烧室内一维质量平均马赫数波动相对较小的原因。

图  7  不同工况的对称面马赫数云图

Fig.  7  Centerline planes of Mach number and the shock system for different cases

下载: 全尺寸图片 幻灯片

与Case 1～3的相似波系结构相比，Case 4～6的流场则呈现出另一种形态，隔离段内受燃烧影响的区域更少，燃烧引起的激波串没有明显前传，激波串前缘基本位于隔离段末端注油位附近。同时，Case 4～6凹槽段的波系相对更强，与其一维质量平均马赫数在燃烧室内波动较大相符。

可见，从Case 1至Case 6，在这个飞行马赫数不断增大、当量比不断降低的过程中，逐渐演变的流场形态从波系结构的角度可以分为两类，其中Case 1～3为一类，Case 4～6为一类。当发动机工况从Case 3变为Case 4时，流场形态的突然变化使得图5中的壁面压力和图6一维质量平均马赫数都出现了突变。

目前，双模态发动机燃烧模态的判定还未形成统一准则。在发动机流场分析中，常采用隔离段内是否形成前传至注油位以前的激波串结构作为区分亚燃与超燃的判据^[10]。由此可认为Case 1～3的发动机处于亚燃模态，而Case 4～6的发动机则处于超燃模态。Case 3、 4之间在流场形态、壁面压力与一维质量平均马赫数上的突变是由于发动机燃烧模态发生了变化。

需要注意的是，从图6中一维质量平均马赫数的分布可以看出，在本文研究的6个工况下，发动机内都存在一维质量平均马赫数小于1的区域。若认为流道内一维平均马赫数等于1的位置（点A）出现热力喉道，则6个工况的发动机都形成了热力喉道，且对于所有的研究工况，喉道位置保持不变，都位于燃烧室与扩张段交接处。

2.2   两种典型的燃烧模态

本小节对该矩形截面发动机不同工况下的燃烧特性进行分析。图8、9分别为不同工况下三维流道内静温（T）与二氧化碳的分布（用质量分数$\omega _{{\rm{CO}}_2} $表示）。二氧化碳作为燃烧产物之一，能一定程度上反映燃烧完成的程度，温度变化更是燃烧反应的重要特性。

图  8  不同工况下流道的静温分布

Fig.  8  Distributions of temperature in the flowpath for different cases

下载: 全尺寸图片 幻灯片

图  9  不同工况下流道的二氧化碳分布

Fig.  9  Distributions of CO₂ for different cases

下载: 全尺寸图片 幻灯片

由图8、9可见，对于低马赫数、高当量比的Case 1～3状态，注油位以前的隔离段内就出现了部分高温区域，温升与二氧化碳的出现表明了燃烧的位置前移，部分燃烧在注油位以前的近壁面区域就已经完成。对于高马赫数、低当量比的Case 4～6状，燃烧反应基本都发生在注油位后的燃烧室内。

需要注意的是，图8中，尽管Case 1～3的发动机内整体静温更高，高温区域更大，但Case 4～6的下凹槽段局部静温却更高，静温沿高度（y方向）分布的梯度更大。图9中二氧化碳分布的发展规律与静温一致，从Case 1到Case 6，随着马赫数的增高、当量比的降低，流道内生成的二氧化碳总量逐步减少，但Case 4～6扩张段内二氧化碳的局部浓度反而更高。二氧化碳浓度在燃烧室后部相对燃烧室前部有更明显的提升，一定程度上体现了与Case 1～3相比，Case 4～6可能有更多的燃烧反应在燃烧室后部完成。

相比于温度和二氧化碳，释热量分布对燃烧的反映更为直接。对于绝热燃烧系统而言，化学反应前后的系统总焓值是恒定不变的，总焓由生成焓与总显焓（静显焓与动能之和）组成。燃烧将燃料生成焓转化为显焓。图10中的释热模型直观描述了此过程。在化学反应之前，燃烧系统的很大部分能量以燃料生成焓的形式存在，化学反应使分子重组。燃烧改变了系统的组分，使反应物向生成物转化。产物的生成焓小于反应物的生成焓，所减少的生成焓用于提升了系统的总显焓。能量的改变和转移表现为燃烧释热，引起了温度和压力的增加，改变了发动机的流动速度。

图  10  释热模型^[14]

Fig.  10  Heat release model^[14]

下载: 全尺寸图片 幻灯片

为分析不同状态下发动机内燃烧释热的特性，本文以总显焓变化量来计算发动机模型沿着长度方向的释热量$H_{{\rm{r}}，x}$，如式（2）所示：

式中：$ {H}_{x} $与$ {H}_{\infty } $分别表示x位置与入口处截面的总显焓。总显焓$ H $可以通过式（3）计算：

式中：积分沿流向的每个截面进行，$ {c_i} $为第i个组分的浓度；n为组分数；$ {H_i}（T） $和$ {H_i}（298.0） $分别为温度为T和298.0 K时第i个组分的静焓；$ {\dfrac{1}{2}}{\overrightarrow {\left| V \right|} ^2} $为气流的动能；$ \overrightarrow V $为计算单元表面的流动速度；${\rm{d}} \overrightarrow {s} $为表面微元的面积向量。

将释热量$ H_{{\rm{r}}，x} $对x进行求导，可获得沿流向的释热变化率H_R，x，如式（4）所示：

计算得到不同工况下发动机释热量与释热变化率沿流向的分布分别如图11、12所示。从图11中可见，Case 1～3的释热量沿流向的变化趋势更接近，释热量增加的位置也更靠前，Case 4～6燃烧释热开始的位置更靠近隔离段出口处。

图  11  不同工况下释热量沿程分布

Fig.  11  The heat release distributions along the flow direction for different cases

下载: 全尺寸图片 幻灯片

图12中释热变化率的分布规律则将两种不同燃烧模态的差异展现得更加明显，Case 1～3的释热变化率从注油位以前的隔离段中就开始逐步增加，在燃烧室入口处（第一凹槽头部）迅速增加直至峰值，然后逐渐下降。而Case 4～6的释热变化率的增长则是从注油位处开始，在燃烧室内均处于相对较高的水平，存在多个峰值，直至燃烧室尾部才逐渐下降。可见在Case 1～3的发动机大部分燃烧释热主要在燃烧室头部快速完成，释热相对更集中；而Case 4～6的释热沿流向分布相对更分散、更均匀。与上文通过二氧化碳分布观察到的结论—相对Case 1～3，Case 4～6有更多的燃烧反应在燃烧室后部完成—相一致。

图  12  不同状态下释热变化率沿程分布

Fig.  12  The heat release rate distributions along the flow direction for different cases

下载: 全尺寸图片 幻灯片

2.3   流动分离的演化

由2.2小节可知，不同燃烧模态下发动机释热变化率分布规律存在明显差异，为了增加对流动和燃烧的理解，本小节通过数值流法对不同工况下发动机内流动分离进行分析。

油流法常用于显示物体表面的流动图谱，对于显示分离、旋涡等复杂流动而言是一种简单有效的手段。通过分析油流谱可以了解流体发生分离的特点、方式和位置，观察旋涡的形成^[23]。油流轨迹线可近似为气流的物面摩擦力线^[24]。图13为不同状态下的数值油流图与壁面静压分布。表面油流线展现了一些流动特性，流谱分析时的奇点是指物面上摩擦应力向量或表面涡量为零的点^[24]，由于摩擦应力分量${\tau _{x，{\rm{w}}}}$、${\tau _{y，{\rm{w}}}}$是非线性函数，物面上可能存在多个奇点。从流动特性上说，一般常见的奇点有鞍点（S）、结点（N）、焦点（F）。其中结点（N）根据轨迹线的方向分为轨迹线指向结点的分离结点（汇结点）和轨迹线背离结点的附着结点（源结点）。类似的，焦点也根据轨迹线方向分为分离焦点（汇焦点）和附着焦点（源焦点）。一般出现焦点说明存在旋涡流动^[23]。

图  13  不同工况下发动机表面油流及静压分布

Fig.  13  Surface oil flows overlaid with static pressure for different cases

下载: 全尺寸图片 幻灯片

从图13可以看到：Case 4～6隔离段内油流线分布均匀，注油位上游没有明显的流动分离，燃烧并没有前传，均在注油位下游完成。而Case 1～3隔离段内油流线的分布不再均匀，上壁面与侧壁面上均出现代表旋涡结构存在的焦点F，可见随着激波串的形成与前移，隔离段流场中形成了明显的分离旋涡。旋涡结构可能会将燃料卷至上游，这解释了部分燃烧在注油位之前已经完成的现象。

对比不同工况下第一凹槽与第二凹槽之间过渡段处的油流线可见，Case1～3分离的流动在第一凹槽倾斜后壁附着结点N处已附着，过渡段内流线均匀，并没有形成明显的分离漩涡结构，可知燃烧室内的分离主要发生于凹槽内部，并没有持续到过渡段。Case 4～6则不同，第一凹槽侧壁均可见焦点F，油流线分布反映了凹槽以及凹槽后过渡段内形成较大的旋涡结构，并且从Case 4到Case 6，旋涡结构流向跨度不断增大，Case 6时已形成了从第一凹槽前缘处持续到第二凹槽处的分离旋涡，可见2个凹槽间也形成了连续存在的流动分离。分离旋涡有助于燃料混合以及燃烧向下游影响传播，这可能是Case 4～6燃烧释热沿流向分布更分散、更均匀的原因。

根据上文所得2类不同燃烧模态下流动分离及燃烧释热变化率的差异可知，在2个凹槽之间过渡段内诱导流动分离，使得燃烧室内形成流向跨度大的分离旋涡结构有助于燃烧向下游传播，从而实现分布式释热，延长释热距离，避免释热过于集中导致激波串前传。

3   结　论

针对所选煤油燃料矩形截面双模态超燃冲压发动机，研究其在不同飞行工况下的流动及燃烧特性，包括流道内壁面压力与一维质量平均马赫数的变化规律，波系结构、释热变化率及流动分离等特性。结论如下：

1） 不同飞行工况下发动机明显工作于2类不同的燃烧模态（亚燃和超燃模态），不同燃烧模态时流场形态的明显变化使得壁面压力与一维质量平均马赫数都出现了突变。

2）当发动机处于预燃激波串前传至注油位以前的亚燃模态时，凹槽段的波系相对较弱，随着激波串的形成与前移，隔离段中形成了明显的分离旋涡结构并将燃料卷至上游，部分燃烧在注油位之前已经完成。燃烧室内流动分离主要发生于凹槽内部，燃烧释热集中在第一凹槽头部，释热相对更集中。

3）当发动机处于激波串未向隔离段上游前传的超燃模态时，凹槽段的波系相对更强，流动参数波动更大。燃烧室内形成的分离旋涡在流向跨度大，从第一凹槽前缘持续到第二凹槽处，过渡段存在连续的流动分离。旋涡有助于燃烧混合以及燃烧向下游传播，因此释热沿流向分布更均匀、更分散。

4） 在凹槽之间过渡段内诱导流动分离，使得燃烧室内形成流向跨度大的分离旋涡结构可能有助于燃烧向下游传播，从而实现分布式释热，延长释热距离，避免释热过于集中导致激波串前传。