某超燃冲压发动机尾喷管的三维流场特性

王建明，齐晓航，栾思琦，王成军，陈宝岭

(1.沈阳航空航天大学航空发动机学院,辽宁省航空推进系统先进测试技术重点实验室,沈阳 110136;2.沈阳航天新乐有限责任公司,沈阳 110034)

超燃冲压发动机作为高超声速飞行器的动力装置，相比于传统的涡喷和涡扇发动机，它在高马赫数下能够拥有较大的推重比[1]。它的结构简单，由进气道、隔离段、燃烧室和尾喷管四部分构成，构型上采用飞行器机体一体化设计，将飞行器的下表面作为进气道和尾喷管的上壁面，以减少迎风面积和重量。超燃冲压发动机尾喷管是提供发动机推力的主要部件，它产生推力的原理是：燃烧室高温高压气体在喷管入口处膨胀加速，利用尾喷管上壁面作为外膨胀面以增大燃气的膨胀面积[2]，让出口处气流和进口处产生动量差，从而使发动机获得净推力。

影响超燃冲压发动机尾喷管性能的因素众多，对尾喷管性能影响规律的研究是近年来该领域的重要研究方向。Huang等[3]做了氢燃料超燃冲压发动机喷管的化学非平衡流动分析，方法是利用CHEMKIN软件对Hyshot型超燃冲压发动机不同工况下喷管内的化学非平衡现象进行了数值分析。通过对冻结流、平衡流和非平衡流三种化学状态的测试和比较，论证了化学反应对喷嘴流场的影响。结果表明：真实非平衡流模拟计算结果介于冻结流和平衡流两种极限情况之间，更接近于冻结流。在喷管入口燃气完全燃烧的情况下，自由基复合反应会释放出大量热量，这些热量主要用于提高气体温度，对推力作用很小。在喷管入口不完全燃烧情况下，喷管内既有燃烧反应，又有自由基复合，反应放热对推力的影响取决于喷管进口处燃烧的完全程度，推力相比于冻结流则大幅度提高。Mo等[4]做了非均匀进气超燃冲压发动机喷管设计与冷态流动的试验。提出了一种考虑非均匀入流影响的非对称超燃冲压发动机喷管设计计算程序。给出了给定非均匀马赫数分布和相应的质量加权平均均匀马赫数分布的典型设计实例。然后，通过计算流体力学分析和冷态流动实验测量，对非均匀设计的喷管性能改善进行了量化。计算和实验的结果表明，非均匀马赫数喷管设计法设计出的喷管较均匀马赫数方法下整体性能更优。Lü等[5]做了考虑横向展开和几何约束的三维超燃冲压发动机喷管设计的研究，提出了一种基于准二维超音速流动和最大推力理论的考虑横向展开和几何约束的三维喷管设计方法。为了生成三维喷管的外形，采用特征线法计算无粘流场以及参考温度法修正附面层厚度，并用计算流体力学方法计算了喷管的气动性能。结果表明：初始圆弧半径对轴向推力系数影响较小，而横向轮廓、下整流罩长度和初始膨胀角的变化对轴向推力系数有显著影响。并和用流线追踪法设计的三维喷管进行了对比研究，验证了新方法的优越性。黄伟等[6]也用数值模拟的方法论述了尾喷管构型对性能的影响，研究表明对于俯仰力矩，上壁面切线角越小越好，外罩长度越长越好，而对于升力则刚好相反，所以在设计尾喷管构型时要充分平衡升力和俯仰力矩的需求。Zhang等[7]做了非均匀入口马赫数的超燃冲压发动机喷管三维流场仿真实验，分析某种构型的喷管内的马赫数和压力等值线，结果表明相比于均匀马赫数入口条件，非均匀入口条件下的推力和陀螺力矩都会降低。

目前关于超燃冲压发动机尾喷管流场的数值模拟大多都是基于二维流场，讨论的内容也是以入口参数对发动机推力、升力等性能参数的影响为主，很少有对喷管流场三维结构的研究，而喷管流场的结构特性也是影响喷管性能的重要因素。现着重于研究喷管三维流场的结构，通过观察各个截面上密度、马赫数等值线、涡量及速度矢量线图来分析三维流场中的波系和涡系的分布特性。

1 计算模型与结果验证

1.1 计算模型和网格

喷管计算模型选自文献[8]，三维模型及结构组成如图1(a)所示，二维结构示意图如图1(b)所示，主要设计参数及尺寸如表1所示，x是主流方向。

计算网格如图2所示，图2(a)是整个计算域的网格，图2(b)是喷管处的网格。在喷管近壁面处的网格均进行了加密处理，网格原点取在喷管入口的左下角。网格大部分由O型网格组成，局部位置采用非结构网格，近壁面第一层网格高度设置为0.005 mm，y+＜1，网格总数约230万。

图1 喷管模型图Fig.1 Model of nozzle

表1 喷管主要设计参数Table 1 Main design parameters of nozzle

1.2 数值方法

流动可视为理想气体可压缩定常流动。在三维坐标系下微分形式的N-S方程如下：

(1)

(2)

(3)

式(1)～式(3)分别是质量、动量和能量的守恒形式方程[9-10]。式中：t为时间；i、j代表方向；ρ为流体密度；p为静压；u为流体速度；E为流体微团总能量；T为温度；gi和Fi分别为i方向上控制体所受的重力体积力和外部体积力；τij为应力张量；k为有效热传导系数；Sh为其他热源项。

图2 计算域与网格Fig.2 Computing domain and grid

现选用基于SST(shear-stress transport)k-ω湍流模型。其中湍动能k的输运方程为

(4)

特定湍流耗散率ω的输运方程为

(5)

式中：Γk和Γω分别为k和ω的有效耗散系数；Gk和Gω分别为k和ω的产生项；Yk和Yω分别为k和ω由于湍流引起的耗散项；Xω是交叉扩散项。湍动能的耗散项Yk可表示为

Yk=ρβkωFDES

(6)

式(6)中：FDES为湍流强度项。

数值模拟用FLUENT软件，求解器选择压力基，计算方法采用Coupled算法，选取Green-Gause节点格式，密度、压力等参数均选取二阶迎风格式做定常计算，共迭代14 000步。设置流体为理想气体，流体黏性采用Sutherland公式。所有壁面皆为绝热、无滑移壁面。根据参考文献设定边界条件为喷管入口静压35 574 Pa，静温298 K，马赫数2.5，自由流静压2 940 Pa，总温673 K，马赫数7.1。

1.3 网格稳定性和计算结果验证

为了保证计算精度的基础上加快计算速度，对网格的独立性进行验证，分别计算了180万、200万、230万和270万的网格，并对计算结果进行无量纲化，计算结果如图3所示，横坐标为喷管上壁面中心线位置坐标，纵坐标是喷管上壁面中心线上压力无量纲化值，根据文献[8],P是上壁面中心线的压力值，Pc是试验台模拟的燃烧室出口气体的压力值，按照文献[8]中无量纲化的方法取(P/Pc)以10为底的对数值。如图3可知当网格在230万时结果已经不再随网格量的增加而改变，因此选取230万的网格进行计算。

仿真结果与文献[8]的实验结果的对比图如图4所示，图4中的横坐标为喷管上壁面中心线坐标，纵坐标是喷管上壁面中心线处压力的无量纲化的数值，图4表明仿真的曲线结果与实验测得的九个测量点的结果基本吻合，进一步说明了仿真结果的准确性和可信性。

图3 网格独立性验证Fig.3 Grid independence verification

图4 数值结果与实验值对比Fig.4 Comparison of numerical and experimental values

2 三维流场结构的结果与分析

2.1 三维流场中膨胀波分布

为了观察三维流场结构，除了截取喷管对称面(z=0.025 m)以外，另取y=0.005 7 m、x=0.02 m、x=0.13 m三处截面，截面位置示意图如图5所示。

流场各个截面的密度分布图如图6所示。根据流场中密度分布情况可以确定膨胀波的位置。图6(a)中，发现喷管内存在沿流向密度下降明显的边界，流场中密度下降、压力必然下降，因此该处即为管内膨胀波的位置。图6(a)中可以看到这两道膨胀波的起始位置分别在下壁面喷管入口处和上壁面内喷管出口处，这些位置产生膨胀波的原因是这两处都是超声速气流发生外折的位置。根据图6(a)中所示两道膨胀波会在内喷管出口处附近交汇，管内气流在经过这两道膨胀波时膨胀加速，进而产生推力和升力。图6(b)中可以看到图6(a)中喷管内膨胀波在y方向上的位置，同时可以观察到产生于两侧壁面出口处的膨胀波，侧壁面产生的两道膨胀波均存在内折角。根据图6(c),可以得知膨胀波之所以发生内折，是因为侧壁尾端的射流边界与侧壁形成了外钝角,而且射流边界处的密度小于喷管出口的密度，即射流边界处压强小于入口处压强，因此会在两个侧壁尾端产生内折的膨胀波。

图5 截面位置示意图Fig.5 Schematic diagram of section position

2.2 三维流场激波和剪切层分布

图7是流场中各个截面的马赫数等值线，由于激波处马赫数等值线会发生汇聚，因此根据马赫数等值线汇聚的位置可以确定激波的大致位置。图7(a)中可以观察到上、下壁面出口处的羽流激波，羽流激波的存在是因为存在高空羽流效应。高空羽流效应的产生是因为超声速喷流离开喷管时会在外部的超声速介质中流动，二者的气流参数不同，此时喷管处的超声速气流在模拟状态下的超声速介质中扩散，就会产生类似于发动机在高空飞行环境下所产生的高空羽流效应，高空羽流效应会导致飞行器产生附体斜激波，这类激波被称为羽流激波[11-12]。在上壁面尾部，由于上方自由流穿过羽流激波，导致上壁面尾部附近的压力明显增高，进而导致此处附近的自由流对主流产生压缩，进而形成一道由压力决定的管内斜激波。图7(b)中，可以看到喷管侧壁面同样存在羽流激波。因此可以推断，外喷管流场周围的羽流激波是呈环状分布在喷管流场周围。并且根据图中马赫数等值线的分布情况可以看出随着流动的发展，激波处马赫数等值线变得越来越稀疏，这表明激波会随着流动的发展强度逐渐减弱。

观察马赫数等值线图，发现除了激波位置外，仍存在其他的等值线发生汇聚的位置。实际上流场中除了激波处，剪切层处的马赫数等值线也会发生汇聚，因此需要判定这些位置上是否存在剪切层。剪切层最明显的特征是它的内部涡量很高，大多数研究者根据这个特征来判定流场中是否存在剪切层。涡量云图如图8所示，根据图8(a)，发现了喷管上、下壁面尾端的剪切层包含z方向的涡。图8(b)中可以观察到喷管主流两侧剪切层则包含y方向的涡。

图6 流场密度云图Fig.6 Cloud map of density

图7 流场马赫数等值线图Fig.7 Contour map of Mach number

2.3 三维流场涡流分布

流向的涡流比较难捕捉，为了确定流场中流向涡的具体位置，做出了各个截面上的速度矢量线图，通过观察各个截面的速度矢量线图来确定涡流的位置，矢量图中速度的大小值取基于该截面下的y方向和z方向的合成速度的值。根据图9(a)，发现流场中有涡流的存在，由于该截面方向垂直于主流，因此这些涡流为流向涡，根据矢量线的分布，发现涡产生于内喷管出口附近。涡流的形成原因是由于内喷管出口处产生的高空羽流效应，使得主流的超声速气流受到羽流的影响，气流会沿横向分离并发生卷折，进而形成流向涡。图9(b)中，上壁面出口处也观察到了流向涡的存在。它的形成原因是由于上壁面的边界层存在横向分离的趋势，这一趋势导致气流离开上壁面后，进而在尾流中生成流向涡[13-14]。

图8 流场涡量云图Fig.8 Vorticity cloud map of flow field

2.4 三维流场主要结构示意图

图10是综合前面结果得到了三维流场结构示意图。图10结果与文献[8]的实验结果相比，文献[8]所得到的流场结构示意图并没有讨论流场中的膨胀波，讨论的流场结构也都是基于二维截面下的。本文揭示的流动结构则是基于三维流场，比文献[8]中的更全面，并且对所有物理现象的产生原因加以论述。

3 结论

通过数值模拟的方法对超燃冲压发动机三维定常流场结构进行分析，结论如下。

图9 速度矢量线图Fig.9 Velocity vector line diagram

(1)喷管内存在膨胀波，管内的膨胀波会发生交汇形成膨胀波面。同时侧壁也产生内折的膨胀波，该处膨胀波会随着气流膨胀，强度逐渐变弱。

(2)喷管外流场周围存在羽流激波和剪切层，羽流激波会形成一个激波面分布在外流场周围，剪切层则分布在羽流激波和主流之间。在上壁面出口处存在羽流激波和因边界层分离产生的斜激波，并且在这两道激波之间也存在剪切层。

(3)内喷管出口处气流会在喷管主流两侧形成流向涡。外喷管出口处气流也会在上壁面出口两侧形成涡流。内喷管出口附近的流向涡主要是受到羽流效应的影响，尾流中的流向涡则是由于上壁面气流沿横向分离造成的。