喷流作用下的单边膨胀后体气动载荷

先进无人作战飞机对全向隐身指标均有较为苛刻的要求，为提高其侧向及后向隐身性能，飞机机身和喷管-后体结构的一体化设计尤为重要。特殊设计的喷管-后体结构一方面使飞机棱线更好地满足平行原则，有效降低雷达散射面积(Radar Cross Section，RCS)指标，另一方面，该结构加强了喷流与外流的冷却掺混，并起到一定的遮挡作用，有效减弱红外辐射[1]。美国波音公司陆续发展的X-45A/B/C[2-3]，诺斯罗普·格鲁门公司发展的X-47B[4]、RQ-180[5]，以法国为主研制的“神经元”[6]以及俄罗斯的“猎人”[7]等无人机，均采用类似设计。

喷管-后体结构在有效提高隐身特性的同时，也在气动、结构、强度等领域带来了新的挑战。从气动角度而言，飞机后体面临复杂且苛刻的载荷环境，例如，高温、高压发动机喷流的强脉动效应，气动载荷与热载荷耦合作用，激波/附面层干涉效应，多股气流的复杂剪切、掺混等。复杂的流动现象给后体气动载荷的预测带来了较大困难，且影响参数众多，依靠单一的仿真或试验手段均无法准确给出所需数据，因此，对后体喷流浸润区内的流动现象和机理的研究是非常必要的。

单勇等[8]利用计算流体力学/红外辐射数值模拟方法，研究了轴对称喷管、单边膨胀喷管的喷流区温度分布及红外辐射特性，结果表明后者可大幅降低红外辐射峰值。张少丽等[9]对不同落压比、不同几何矢量角下的单边膨胀喷管进行分析，研究了喷管推力特性和红外辐射强度的变化规律。徐嘉等[10]采用数值模拟方法研究战斗机后体绕流与喷流相互耦合作用的流场特性，对欠膨胀喷管及后体尾部波系结构进行分析，并进行减阻优化设计。杨承宇[11]利用数值模拟方法研究了膨胀边开缝、倾斜角度、面积比等参数对单边膨胀喷管气动性能和红外辐射特性的影响规律。李春鹏等[12]针对无尾布局超声速航向静稳定性不足的问题，提出一种基于超声速压缩/膨胀流动的后体超声速航向气动增稳设计方法，研究了后体脊线形状、截面形状等参数的影响。

上述研究主要基于稳态CFD分析，研究目的集中在评估隐身特性或提高气动性能，从公开发表的文献来看，现阶段还鲜有专门针对喷流作用下后体气动载荷的相关研究。

国内外对超声速射流、剪切层掺混等基本流动现象开展了相关研究，其研究方法和部分结论同样适用于喷流气动载荷问题。朱文庆等[13]使用修正长度尺度的改进的延迟分离涡模拟(Improved Delayed Detached Eddy Simulation，IDDES)方法计算了飞行速度对喷流流场和远场噪声的影响，并与试验结果[14-15]进行对比验证了方法的正确性。朱志斌等[16]利用大涡模拟(LES)方法研究了后台阶喷流混合流场，捕捉到了流场激波系以及混合剪切层，并获得了流场结构特征及密度脉动特性。刘旭亮和张树海[17]利用直接数值模拟(DNS)方法模拟了二维激波与剪切层的相互作用，对复杂的流场和声场结构进行了分析并揭示噪声产生的机理。Tam等[18-19]采用计算气动声学方法求解线化雷诺平均Navier-Stokes(N-S)方程，得到了不同马赫数下双涵同轴喷流激波系结构，并与试验结果[20]进行了相关性分析。李栋等[21]在数值模拟方法研究基础上对矢量推力飞机喷流-外流干扰流场进行了分析。

本文以高隐身无人机喷管-后体结构为研究对象，结合高/低速风洞试验和CFD仿真方法，对喷流作用下飞机后体的稳态、动态气动载荷及其影响规律进行研究，也对流动现象及机理进行了相关分析。

1 计算方法及数模

1.1 计算方法

计算求解笛卡尔坐标系下的三维积分形式雷诺平均N-S方程：

(1)

式中：Q为守恒量；Fc为对流通量；Fv为扩散通量；S为源项；Ω为控制体积；∂Ω为控制体积的边界面；n为法向量，指向外侧。空间离散采用基于非结构网格的有限体积法，对流通量采用基于Roe的迎风格式进行离散，黏性通量采用中心差分格式进行离散。

本文研究包含稳态计算和非稳态计算两部分。稳态计算选择k-ω 剪切应力输运(SST)两方程湍流模型，基于全隐式时间推进法进行求解；非稳态计算采用基于k-ω SST湍流模型的IDDES方法[22]，该方法属于RANS/LES混合方法的一种，有效结合了二者的优点，计算精度较高且计算量可接受，得到了较为广泛的应用[23]，能解决喷流计算中包含边界层流动，主、次流掺混等工程问题。

1.2 数模及网格

图1为单边膨胀喷管-后体结构示意图(半模)，喷管为圆转方形式，收缩比为0.35。后体与喷管为套接形式连接，部分气流(来自发动机舱或飞机其他舱室)由主喷流引射流出。飞机后体下膨胀边长度约为上膨胀边的2.8倍。在真实飞机上，膨胀边构型可能是较为复杂的曲面，流道面积略微扩张。

图1 单边膨胀喷管-后体结构示意图
Fig.1 Schematic diagram of single expansion nozzle after-body structure

计算采用带棱柱层加密的非结构化网格，棱柱层第1层高度为1×10-5m，棱柱层膨胀比为1.3。对重点关注的喷流影响区域进行加密，稳态计算采用半模，喷流区域网格尺度为0.02 m， 网格总量在400万量级(经验证，网格无关性良好)；非稳态计算采用全模，喷流区域网格尺度为0.01 m，网格总量在3 000万量级。不关注传热特性，计算域内气流均设为冷流，与风洞试验条件一致。

2 稳态气动载荷分析

2.1 计算结果验证

喷流试验分别在航空工业空气动力研究院某2.5 m量级低速风洞和某1.5 m量级高速风洞完成，风洞试验模型如图2所示。模型采用腹部支撑，内接杆式天平，喷管压比(NPR，喷管进口总压与环境静压的比值)连续可调。后体上、下膨胀边布置有稳态和动态传感器，对壁面压力进行测量。稳态壁面压力测量的精度满足国军标相关要求[24]，误差在5‰范围内。

图2 风洞试验模型(局部)
Fig.2 Wind tunnel test model (local)

图3为数值计算得到的下膨胀边稳态气动载荷分布，马赫数Ma=0，NPR分别为2.4、3.9。定义喷流气动载荷等于壁面静压与环境静压之差。可以看出，超声速喷流形成的复杂激波-膨胀波系，在下膨胀边投影为压力、吸力交替的气动载荷，且随着NPR增加，载荷分布形式更为复杂。

图3 下膨胀边气动载荷分布(Ma=0)
Fig.3 Aerodynamic load distribution of lower expansion edge (Ma=0)

图4为Ma=0，不同NPR时数值计算得到的气动载荷分布曲线(位于Z=0 m截面，即对称面)与风洞试验值的对比，分别对应于无来流状态以及Ma=0.85状态，实线为数值计算结果，符号为风洞试验数据。横坐标X/L表示X方向相的对位置，L为下膨胀边长度，纵坐标Cp为喷流气动载荷基于环境静压无量纲化后的系数。可以看出，计算得到的载荷曲线与试验数据吻合较好，对于不同来流条件下的波系结构预测准确，验证了计算方法的有效性，为进一步分析提供依据。

 

图4 计算结果与风洞试验结果对比
Fig.4 Comparison between calculated results and wind tunnel test results

图5为数值计算得到的Ma=0.85、NPR=3.5状态下膨胀边稳态气动载荷分布，可以看出，喷管出口附近区域的压力、吸力交替载荷相对较大，而之后区域尽管仍存在正、负交替的波动现象，但与Ma=0时相比(见图3(b))，其幅值明显减小。推测是由于外流Ma增加引起的，后续将针对Ma的影响进行单独分析。红色虚线位置实际为上膨胀边斜切尾缘在下膨胀边的投影，可以看出，压力云图也存在与之相吻合的三角形区域，说明上膨胀边几何形状会对载荷分布产生直接影响。

图6为不同Z值截面位置对应的载荷曲线与试验数据的对比，截面位置见图5。可以看出，Z=0.2 m和Z=0 m位置相比，载荷极值未明显减小，这是由于该截线同样穿过了出口强激波、强膨胀波作用区域(见图5)。随着展向距离增加，载荷曲线波动幅度逐渐减小，在Z=0.4 m位置已经近似均匀。在工程实际中，可以设一个波动阈值，人为划定喷流影响区域的范围，便于数据处理。上、下膨胀边不同Z截面位置的数值计算结果与试验值吻合较好，进一步证明了计算方法的准确性。

图5 下膨胀边气动载荷分布(Ma=0.85,NPR=3.5)
Fig.5 Aerodynamic load distribution of lower expansion edge (Ma=0.85,NPR=3.5)

  

图6 计算结果与风洞试验结果对比(Ma=0.85,NPR=3.5)
Fig.6 Comparison between calculated results and wind tunnel test results(Ma=0.85,NPR=3.5)

2.2 喷管压比的影响

为进一步分析喷管NPR对壁面气动载荷的影响规律，选择相对密集的NPR序列进行计算。图7为无来流NPR在1.5～4.0范围内变化的气动载荷曲线。可以看出：

图7 喷管NPR对载荷分布的影响(Ma=0)
Fig.7 Influence of nozzle NPR on load distribution(Ma=0)

1) 喷管出口附近的内膨胀段(X/L<0.24区域)分别对应壁面的压力、吸力载荷极值点，且在NPR>3.0时，载荷极值随着NPR增加近似线性增大。

2) NPR高于2.5时，载荷曲线逐渐出现正负波动的规律，对于较小NPR，载荷可忽略。

3) 外膨胀段(X/L>0.24区域)载荷方向随着NPR变化敏感，例如NPR为3.0和3.25，载荷方向几乎全部相反，原有压力区变为吸力区，反之亦然；这也说明，后体气动载荷的相关分析要与发动机工作状态相结合。

图8为典型NPR状态下喷流区域对称面静压分布，可以较为清晰地反映出喷流形成的波系结构，经分析可得如下结论：

图8 典型状态下对称面喷流区静压分布
Fig.8 Static pressure distribution of jet region on symmetrical surface under typical conditions

1) NPR为2.5时，喷流形成的波系强度较弱，高压区、低压区交替出现；当NPR为3.5和4.5时，波系强度增加，且壁面反射作用及波系间干涉作用增强，压力分布更为复杂。

2) 尾缘附近区域受NPR影响较大，如虚线方框所示，NPR为3.5时，结尾处为高压区在前低压区在后，而NPR为4.5时，变为低压区在前高压区在后，波系结构的变化导致壁面载荷方向发生变化，这与图7的结论是一致的。

2.3 喷流与外流掺混的影响

图9为NPR=2.9时不同来流速度下膨胀边气动载荷分布曲线，Ma在0～0.8范围变化。可以看出：

1) 来流Ma从0增加到0.2对载荷分布几乎无影响，Ma>0.4时，其影响逐渐增加。

2) Ma在0～0.6范围内，其作用效果是使波系逐渐前移，但载荷分布规律基本不变，而Ma为0.8时，其分布规律也有较大变化。

3) 来流Ma增加的整体效果是使吸力载荷大幅减小，使压力载荷小幅增加，对于图中标识的a、b两区域的作用尤其明显，需要结合流场进行分析。

图10为不同Ma下的喷流区域静压分布，流场参数与图9中Ma=0.2，0.6，0.8 这3个状态相对应。可以看出，来流速度对喷流区域静压分布有一定影响，主要体现在以下几个方面：

1) 整体来看，随着Ma增加，喷流区域静压值有一定程度的增加，对膨胀波后的低压区影响更大，而对激波后的高压区影响相对较小。

2) 喷管出口的第1道膨胀波-激波区域规律性更加明显，如图中红色方框所示，随着外流Ma增加，膨胀波明显减弱，低压区范围变小，但激波后的高压区差异不大；膨胀波、激波分别与图9中a、b位置对应。

3) Ma增加使喷流与外流的掺混效果增强，波系结构也有一定差别，Ma=0.2波系间的界线比较清晰，而Ma=0.8的界线相对模糊；体现在图9压力载荷上，即前者曲线拐点尖锐，载荷梯度较大，而后者更为平缓，载荷梯度较小，这在强度预测中也是需要注意的问题。

图9 来流速度对载荷分布的影响(NPR=2.9)
Fig.9 Influence of incoming flow velocity on load distribution(NPR=2.9)

图10 不同Ma下的喷流区域静压分布
Fig.10 Static pressure distribution of jet region under different Ma

2.4 次流影响及引射特性

图11展示了套接段次流流量变化对下膨胀边气动载荷的影响，对比了无次流状态，以及次流、主流的流量比在1.0%～4.0%范围内的几种情况。可以看出：

图11 次流流量对载荷分布的影响(Ma=0,NPR=2.9)
Fig.11 Influence of secondary flow rate on load distribution(Ma=0,NPR=2.9)

1) 喷管套接段(X/L<0区域)及内膨胀段，载荷受次流影响比较大，较大次流流量对应较低的吸力载荷，同时也说明，主流对次流有较强的引射能力，即无次流时套接段为较强负压。

2) 随着次流流量继续增加，波系整体有向前移动的趋势，次流的影响相对于NPR和Ma均较弱，不会引起波系结构的变化。

3) 较大的次流流量可以有效降低套接段及喷管出口附近的吸力、压力载荷极值，而使其他区域的载荷极值小幅度增加。

综合以上分析可得知：喷管压比、外流速度、次流流量等多个参数均对喷流流场产生影响，可改变喷流波系结构和空间静压分布，进而影响壁面气动载荷分布；且多个参数之间可能存在复杂的耦合作用效果，在对后体上、下膨胀边进行强度分析，尤其是疲劳、寿命相关的分析时，要注意多参数的复合作用效果，应与飞行包线、任务剖面、发动机状态等相结合。

3 动态气动载荷分析

3.1 计算设置与数据处理方法

对于动态气动载荷的计算分析，网格尺度设置为飞机特征长度的1‰量级，认为尺度更小的涡对后体壁面动态载荷影响较小，可近似忽略。综合考虑网格尺度及流动特征速度，设置非稳态计算的物理时间步长为1×10-4s，采用双时间步方法进行计算。

非稳态计算重点关注上、下膨胀边喷流作用区域内的动态载荷情况，如图12所示，在壁面上密集布置260个监控点，记录一段时间内的压力脉动，以便于处理成功率谱密度(PSD)分布或声压级(SPL)分布。

图12 壁面动态监控点布置(局部)
Fig.12 Layout of wall dynamic monitoring points (local)

喷流作用下的壁面压力脉动为随机信号，无法直接进行傅里叶变换，也不能得到理论上精确的功率谱密度函数，只能用有限的样本数据来近似，因此结果依赖于对数据的处理方法和拟合的数学模型等。经研究，认为现代谱估计方法中的AR模型[25-26]方法适用于喷流压力脉动问题，本文用该方法进行数据处理。

3.2 计算结果分析

图13为H=0 m、Ma=0、NPR=2.9状态，Z=0 m截面稳态计算结果与非稳态时均结果的对比，其中H为飞行高度，非稳态结果的监控点与图12中Z=0 m截面监控点对应。可以看出，整体上两种方法计算结果吻合较好，一定程度上可相互验证准确性。

图13 稳态计算与非稳态时均计算载荷分布对比(Ma=0,NPR=2.9)
Fig.13 Comparison of steady computation and time-averaged unsteady computation load distribution(Ma=0,NPR=2.9)

图14为H=0 m、Ma=0、NPR=2.9状态，非稳态计算得到的某时刻喷流作用区域静压分布，分别展示了XY截面(对称面)及ZX截面(通过喷管中心的横截面)的结果。

图14 某时刻喷流作用区域静压分布
Fig.14 Static pressure distribution of jet action area at a certain time

分析流动特征，认为喷流区域存在两重剪切掺混效果：首先，喷流主流和从套接段引射而来的次流形成“第一重掺混”，主要发生在内膨胀段，但掺混形成的脱落涡会贴近壁面继续向后传播；其次，喷流主流和引射次流掺混后的喷流，再与外流形成“第二重掺混”，第二重掺混的效果往往比第一重掺混更强烈。

图14中用虚线标识出了关键的剪切掺混区域，实际为上述的“第二重掺混”。剪切掺混形成了一系列脱落涡，具有两方面的影响：一方面，高、低压脱落涡交替扫过某区域，直接形成了压力脉动；另一方面，涡的迁移又直接或间接对喷流核心区域的压力分布产生影响，即与波系产生干涉，会形成更强的压力脉动。而气流压力脉动直接作用在后体表面，形成较大的动态气动载荷。

图15为H=0、Ma=0、NPR=2.9状态，非稳态计算和稳态计算得到的XY截面静压分布，其中图15(a)是图14(a)喷流区域的局部放大。与稳态计算结果相比，非稳态计算得到的流场反映了更多的流动细节，将图示范围内的流场分为 A～E 5个区域：A为喷流与外流上掺混区，B为喷流与外流下掺混区，C为喷流核心区，D为喷流主流与次流上掺混区，E为喷流主流与次流下掺混区。

图15 稳态、非稳态流场结构对比(XY截面)
Fig.15 Comparison of steady and unsteady flow field structures(XY section)

实际上，D、E可认为是“第一重掺混”区域，而A、B则为“第二重掺混”区域。由图15(a)和图15(b)对比分析，可以得到如下结论：

1) 对比两图的C区，与稳态结果相比，非稳态结果的高压、低压区域轮廓更明显，波系结构更清晰；而稳态计算结果则存在较强的耗散，距离喷管较远的区域已无明显的高、低压波动，但实际上该区域的瞬时气动载荷仍可能较大，稳态结果并不能反映真实情况。

2) 对比两图的A、B区，稳态结果该区域压力分布较为均匀，与外流压力相当，没有反映出剪切掺混的效果；而非稳态结果则存在多个高、低压间隔分布的脱落涡，剪切掺混效果明显；其原因是，该区域流动相当不稳定，经时间平均，高、低压脱落涡的贡献近似抵消，因此稳态计算得到了均匀的压力分布。

3) 对比两图的 E区，非稳态结果中存在贴近下壁面的较薄剪切层，是喷流主流与次流掺混形成的，也存在快速向后传播脱落涡，导致壁面附近形成压力脉动；而稳态结果中则没有反映出该剪切层。

4) 对比两图的 D区，其差别不大，尽管上壁面也存在喷流主流与次流的掺混，但由于内膨胀段中主流波系较强，高、低压区域相对稳定，掺混没有形成明显剪切层。

图16为非稳态计算和稳态计算得到的XZ截面静压分布，图16(a)为图14(b)的局部放大，图中分别标识出了两侧的掺混区和喷流核心区，该掺混区主要在外膨胀段，为“第二重掺混”区域。分析图16，可以得到与图15相同的结论：一方面，稳态计算结果中喷流核心区存在较强的耗散，压力分布很快就趋于“均匀”；另一方面，稳态计算结果无法模拟掺混区的流动现象，这也是由于高、低压脱落涡的时均贡献近似抵消导致的。

图17为某时刻喷流作用区域涡量分布，涡量定义为空间速度场的旋度，可以比较清晰地反映出流场内剪切掺混的强度。可以看出，喷流核心区内涡量值较小，而喷流主流与次流的掺混区域，以及喷流与外流的掺混区，涡量值均较大，这与图15和图16的分析是吻合的。

图16 稳态、非稳态流场结构对比(XZ截面)
Fig.16 Comparison of steady and unsteady flow field structures(XZ section)

图17 某时刻喷流作用区域涡量分布
Fig.17 Vorticity distribution of jet action area at a certain time

在高、低速风洞试验中对膨胀边壁面的动态气动载荷进行测量，图18为计算得到的声压级(SPL)频域特性与对应测试点试验结果的对比，状态为Ma=0、NPR=2.9。可以看出，计算结果整体与风洞试验结果吻合较好，在f=10～1 000 Hz频率范围内，分布规律和量级预测较为准确。高于1 000 Hz的部分属于噪声范畴，已不是气动载荷所关注的范围。

图18 计算得到的SPL与试验测量结果对比
Fig.18 Comparison of calculated SPL with experimental results

图19为非稳态计算得到的后体下膨胀边总声压级分布云图，圆圈为计算时设置的监控点(见图12)，总声压级反映接收点压力脉动的强弱特性，可由脉动压力的时域信号的均方根值(即总有效值)计算得到。从图19可以看出，动态载荷极值区对应于喷流与外流强烈掺混的位置，而喷管出口及喷流核心区等流动相对稳定的区域，其动态载荷相对较低。壁面动载荷强度的分布规律与图14～图17的分析是吻合的。

图19 总声压级分布云图(半模)
Fig.19 Contour of overall SPL distribution (half model)

4 结 论

对高隐身无人机喷流作用下单边膨胀后体的稳态、动态气动载荷的分布规律，以及相关流动机理进行了研究，主要结论如下：

1) 超声速喷流往往会形成激波-膨胀波交替的复杂波系，投影到后体上、下膨胀边壁面上会成为压力、吸力交替的载荷形式，且载荷分布随NPR变化较为敏感，NPR小幅度变化可能导致载荷方向发生变化。

2) 外流与喷流的掺混会影响壁面载荷分布，外流Ma增加的影响趋势是使吸力载荷大幅减小，而使压力载荷小幅增加。

3) 引射次流主要对套接段以及内膨胀段区域产生影响，次流流量增加的整体作用是使吸力、压力载荷极值均降低。

4) 基于IDDES的非稳态计算方法适用于喷流动态气动载荷预测，捕捉到了多股流动的复杂剪切掺混现象，而剪切掺混是形成压力脉动的关键。

5) 非稳态计算与稳态计算得到的喷流流场有明显差别，前者能够捕捉到喷流与外流的剪切掺混，而后者则不能；其原因是脱落涡的时均贡献近似抵消，稳态计算只能得到近似均匀的压力分布，不能反映真实瞬态流场。

6) 将流场划分为喷流主流与次流掺混区、喷流核心区、喷流与外流掺混区等几个区域，有助于从流动机理上解释脉动压力的分布规律及其成因，且得到了试验验证；动态载荷的极值区出现在喷流与外流强烈掺混的位置，而喷管出口及喷流核心区等流动相对稳定的区域，其动态载荷相对较低。