基于大涡模拟的热喷流地面效应气动特性研究

陈鹏飞，史亚林，杨 华，王 衡，刘志友

(1.中国航发四川燃气涡轮研究院，四川 绵阳 621000；2.西北工业大学 动力与能源学院，西安 710129)

1 引言

近年来民航领域蓬勃发展，但由此带来的噪声严重影响了机场周围居民的生活。当前高性能航空发动机工作时燃烧室的温度可达2 000 K 左右，排气温度升高导致的喷流辐射声强相关问题十分突出[1-4]。上世纪70 年代Hoch 发现，温度升高使低速喷流的辐射声强升高和使高速喷流的辐射声强降低的现象。实际上飞机往往是在起飞和降落阶段产生的噪声最为明显，在这些阶段，喷管工作状态大，飞机近地飞行，喷管排气受到地面阻挡，为有限空间排气，其流场更为复杂。目前，关于高温下喷流近场的流动特性及近场与远场的相关性的研究还不够广泛[5-8]，研究热喷流地面效应气动特性是理解温度造成喷流辐射声场特性变化的前提。

国外学者中，Hussian 等[9]研究发现，喷射流动中相干结构的发展由卷吸、配对及破碎等一系列过程组成。Thurow 等[10]观察到，大尺度涡结构的尺寸相当或大于喷嘴直径，大尺度涡结构是主要的流体动力学含能结构。国内学者中，李经警等[11-13]为揭示二元圆转矩喷管尾喷流强化掺混的内在机制，应用大涡模拟(LES)方程对两种相同进、出口直径的喷管模型(轴对称、二元圆转矩)进行了数值模拟计算，结果表明：与轴对称喷管相比，圆转矩喷管射流掺混效应增强，速度衰减快，核心区长度和高温区域面积减小，同时尾喷流拟序结构(包含涡环、涡辫、发卡涡、螺旋涡)发生变化。因而，从机理上弄清有限空间喷流的发展过程及掺混过程的影响因素极为必要，不仅可以为研究喷流地面效应机理提供依据，也能为建立近场辐射模型提供理论基础。

本文主要采用数值计算的方法，建立发动机尾喷管排气模型，对超声速轴对称喷管在地面影响下有限空间(Condition A，喷管中心离地面高度一般小于5 倍喷管直径)和自由空间(condition B，喷管中心离地面高度一般大于5 倍喷管直径)内喷流流动特性进行数值模拟研究。研究采用可压缩流大涡模拟控制方程动态亚格子模型进行，对喷管尾喷流核心区长度与卷吸掺混面积变化、涡量和雷诺剪切应力分布，以及射流场拟序结构发展变化等方面进行了比较分析。

2 模拟算法及控制方程

数值模拟采用LES 方程的动态Smagorinsky模型，非稳态条件。利用滤波将瞬时变量φ(x,t)划分为大尺度量(x,t) 和小尺度量φ′ (x,t)，(x,t) 通过以下加权积分得到：

式中：G(x-x′,Δ) 为滤波函数；Ω 为计算区域；Δ为滤波的宽度，与网格分辨率有关。

式(1)表示的滤波函数处理瞬时状态下不可压缩流N-S 方程时，有：

连续性方程

动量方程

能量方程

式中：µSGS为涡黏系数。

3 计算模型及网格

3.1 物理模型

为简化计算工作量，研究模型为轴对称拉法尔喷管，如图1 所示。喷管出口直径为Dj，喷管中心离地高度为H，坐标原点位于喷管进口圆心处。

图1 喷管物理模型Fig.1 Physical model of nozzle

3.2 求解方法

计算域如图2 所示。长度方向为 20Dj，宽度和高度方向均为 10Dj。流动方向为X轴正向，宽度方向为Y轴方向，高度方向为Z轴方向。其中Condition A 状态H=2Dj，Condition B 状态H=5Dj。

图2 不同离地高度状态下喷管计算域Fig.2 Calculation domain of nozzle at different ground clearance

具体边界条件为：喷管进口为压力进口，总温为870 K，总压为32 300 Pa。喷管壁面为绝热无滑移条件。Far1、Far2、Far3、Far4 为压力远场。模型采用半模计算，对称面边界条件为Symmetry。Condition A 状态时，Ground 为无滑移条件。Condition B 状态时，Far5 为压力远场。声学马赫数基于喷口速度及远场声速，喷管出口马赫数为1.526。

采用可压缩流动基于密度的方法求解流动控制方程。初场由雷诺时均方法结合SSTk-ω湍流模型得出。经过时长为400 个时间步长的发展，认为流场已排除初始流场的影响，之后开始进行湍流量的采样统计。计算总时间为0.02 s，取最终0.02 s 的计算结果进行流场分析。

3.3 网格划分

为更好地捕捉射流与尾喷流相互作用区域的流动细节以及地面效应，在喷管壁面处和地面处对网格加密处理，根据Y+<1，得出第一层附面层厚度为0.01 mm。远场网格随着X值增大变大。经网格无关性检验和考虑模拟精确程度，最终网格数为600 万。不同状态下的网格划分如图3 所示。

图3 不同状态下网格划分Fig.3 Grid generation in different states

4 计算结果及分析

4.1 温度分布

轴向截面温度分布如图4 所示。中间温度高的区域为喷流核心区，由于喷管出口为超声速，导致其核心区不会与亚声速喷流核心区一样光滑。Condition A 的核心区长度大约为 11.8Dj，核心区在X/Dj=10 位置有个间断面，在此位置地面开始对喷流产生较大影响，导致地面温度升高。由于地面对喷流的影响，在喷管轴线以上的区域，喷流与环境掺混区域增大，掺混作用增强。Condition B的核心区长度大约为 10Dj，核心区在X/Dj=6 位置时有个间断面。

图4 轴向截面温度分布Fig.4 Temperature distribution of axial section

4.2 速度分布

模型喷流中心线上无量纲速度(U/Uj，Uj为喷管出口截面平均速度)分布如图5 所示。可见，在喷管内部，2 个状态速度逐渐增加；在喷管出口初始段，存在1 段速度变化不大的区域(核心区)；在核心区内，随着激波和膨胀波在自由边界的交替发展，速度交替增大和减小，但总体上在一定范围内波动；随着射流与环境之间能量、动量交换增强，射流动能降低，速度逐渐减小。在X/Dj=9 之前，Condition A 和Condition B 喷流中心线上无量纲速度几乎相同；在X/Dj=9 之后，由于地面对Condition A 喷流的影响，使得其喷流与空气的掺混相较于Condition B加剧，速度整体来说降低得更快，如图6 所示。

图5 喷流中心线无量纲速度分布Fig.5 Non-dimensional velocity distribution of jet centerline

图6 Condition A 对称面马赫数分布Fig.6 Mach number distribution of symmetry plane of condition A

4.3 掺混特性

采用大涡模型对射流掺混进行模拟，能够描绘出其射流拟序结构随时间和空间的变化，喷流剪切层的发展与旋涡运动密切相关。Q准则在涡识别中计算效率较高，效果明显，是一种常用的涡量识别方法，本文采用Q准则作为旋涡可视化的工具。Q准则定义式为：

根据Q准则显示出的尾喷流瞬时拟序结构如表1 所示，图中颜色由速度场着色得到。可见，流场中涡旋结构主要由涡环、涡辫、发卡涡、螺旋涡组成，Condition A 和Condition B 射流场涡核发展过程大致相同。

表1 不同时刻涡结构Q 等值面对比Table 1 Comparison of Q equivalent surfaces of vortex structures at different times

以Condition B为例，当射流从喷管出口流出时，与环境气流掺混作用较弱，但由于射流剪切层比较稳定，诱导出的涡结构较小，存在一段光滑区。在射流向前发展的过程中，气流向内螺旋汇聚，破坏速度剪切层，形成涡环结构；其在剪切作用下会逐渐拉伸，而后脱落(脱落频率受剪切层脉动影响)，随后向内卷吸的气流在射流剪切作用下，又形成新的涡环。同时，受涡环外侧反向速度拉伸，使连接2 个涡环之间的涡管结构生长——这部分结构被称为涡辫。该结构是在流动向下游发展的过程中，由于涡环之间的非线性不稳定作用增强，射流脉动较强的条件下逐渐形成的。随着流动进一步发展，射流柱在剪切作用下脉动特征加强，涡环之间距离缩小，涡环与涡辫之间相互作用增强，最终形成尺度较大的发卡涡，使射流柱表面呈鱼鳞状；随着射流继续向下游发展，射流脉动减弱，大尺度发卡涡被耗散成更多小尺度的螺旋涡结构。

从Condition A 可以看出，当X/Dj=9 左右时，地面开始对喷流产生影响。对比t=0.010 2 s 和0.020 0 s 时刻，Condition A 比Condition B 在轴向更早出现发卡涡。Condition A 时，在地面有限空间内的湍流脉动更加强烈，喷流在地面的作用下诱导出二次涡流，涡核分布更加集中，且其中以大尺度发卡涡为主。从Condition B 在t=0.020 0 s 时刻的涡核结构可以完整地看到涡环产生-脱落-产生，涡环反向运动将涡管拉伸为涡辫，随着涡环距离减小，发卡涡形成，但最终随着湍流脉动的耗散，发卡涡破碎为螺旋涡。

4.4 雷诺应力分布

雷诺应力能表征湍流脉动强弱，因此通过比较2 个状态雷诺剪切应力特征，可以反映其射流剪切层内流体脉动状况，进一步反映地面对喷流掺混特性的区别。水平面无量纲雷诺剪切应力为，垂直面无量纲雷诺剪切应力为。其中，u为X向脉动速度，v为Y向脉动速度，ω为Z向脉动速度。

图7(a)、图7(b) 分别给出了Condition A 和Condition B 在各截面水平中心线的雷诺应力分布。可见，两者雷诺应力分布趋势类似。在不同截面，雷诺应力从中心轴线随着Y值的增大先增大后减小，在中心线处值为0，极值只有1 个，分布曲线呈现单峰状态。

图7 各工况截面中心线雷诺应力分布Fig.7 Reynolds stress distribution at the center line of section under various working conditions

在Condition B 状态下，随着X/Dj增大，雷诺应力的极大值先增大后减小，极大值对应的位置向中心线方向移动。雷诺应力分布趋势的内在机理为：以X/Dj=5 为例，中心线(Y/Dj=0)上的雷诺应力最小，随着Y值增大而逐渐增大，在剪切层边界附近达到最大。在这个区域附近，射流与外界气流发生剧烈掺混，能量耗散加剧，而后剪切应力逐渐减小。随着X/Dj增大，射流雷诺应力的极大值先增大后减小。X/Dj=3 时，该截面距离喷口很近，射流处于光滑段，雷诺应力很小，符合前文提到的光滑段拟序结构少的特征。X/Dj=5 时，雷诺应力达到最大，此时射流处于过渡段，射流脉动最强，也佐证了射流在过渡段表面旋涡尺度增大、主要成鱼鳞状的特征。X/Dj=12 以后，雷诺应力随着大尺度涡被耗散，射流脉动减小，也验证了此区域以较小尺度的螺旋涡结构为主的拟序特征。在Condition A 状态下，X/Dj=9 时的雷诺应力极大值大于X/Dj=5 时的极大值。其内在机理为：X/Dj=9 时，地面对喷流的影响较大，射流脉动在有限空间内被挤压，地面对射流产生了扰动。

图7(b)、图7(c)分别为Condition B 各截面水平中心线雷诺应力和Condition A 各截面垂直中心线雷诺应力分布图。Condition B 为自由射流，各截面水平和垂直中心线雷诺应力保持一致。Condition A 各截面垂直中心线雷诺应力分布呈现双峰状态。Z/Dj<0 时，两者雷诺应力分布趋势一致，但Condition A 各个截面垂直中心线位置处的雷诺应力比Condition B 的大，其分布曲线更加饱满。Z/Dj<0 时，Condition A 各个截面垂直中心线位置处的雷诺应力比Condition B 的大，且各截面雷诺应力在Z向范围更大。计算结果表明，由于地面有限空间对喷流的影响，使得湍流脉动更强，并且使得湍流脉动被挤压至喷管中心线以上的区域，这也佐证了X/Dj=12 时垂直面速度分布特征。

5 试验验证

通过发动机热喷流场整机试验，测得发动机离地2Dj高时喷管排气地面X/Dj=12 处地表的最大温升约为30 K，其排气流场瞬时红外热像图如图8所示。可以看出，试验和仿真结果一致，充分表明了本文计算模型的有效性。

图8 热喷流场瞬时红外热像图Fig.8 Instantaneous infrared thermal image of hot jet field

6 结论

采用大涡模拟计算研究地面热喷流气动特性，建立不同离地高度的超声速喷管模型，对总温为870 K、马赫数为1.526 的热喷流状态开展了数值模拟，分析了有限空间喷流的发展过程。主要结论如下：

(1) 自由空间喷流的核心区长度大约为10Dj，在X/Dj=6 位置时有个间断面。发动机有限空间喷流受地面阻挡，排气能量在地面一侧堆积，能量释放区域变长，核心区长度增大为11.8Dj，核心区间断面推后到X/Dj=10 位置，在此位置地面开始对喷流产生较大影响，导致地面温度升高；在X/Dj=9 之前，有限空间喷流与自由空间喷流的无量纲速度几乎相同，在X/Dj=9 之后，由于地面对喷流的影响，使得喷流与空气的掺混加剧，其速度整体来说降低得更快。

(2) 有限空间与自由空间喷流场涡核发展过程大致相同。在流场中，涡旋结构主要由涡环、涡辫、发卡涡、螺旋涡组成；在X/Dj=9 左右位置时，地面开始对喷流产生影响，在地面有限空间内的湍流脉动更加强烈；喷流在地面的作用下诱导出二次涡流，其涡核分布更加集中，且其中以大尺度发卡涡为主。

(3) 通过雷诺应力和试验验证表明了计算模型的有效性。喷管排气地面X/Dj=12 处地表最大温升约为30 K，为研究热喷流声场测试提供了参考。