欠膨胀超声速冲击射流的大涡模拟

2022-07-06郑枫弋赖焕新

华东理工大学学报(自然科学版) 2022年3期

郑枫弋，赖焕新

（华东理工大学机械与动力工程学院，上海 200237）

超声速冲击射流往往伴随着强烈的传热传质现象，可以提供较大的冲击力，广泛应用于航空航天等诸多领域[1]。在垂直起降飞行器（STOVL）上，可以利用机身底部的喷管喷射的高速气体对地面的冲击来获得飞行器所需要的向上推力，从而实现垂直起降和悬停等动作。超声速冲击射流的流场结构复杂，激波、剪切层与边界层存在强烈的干涉，其湍流结构是相关应用中流体动力学的关键。

根据Mitchell 等[2]的可视化研究，可以简述欠膨胀超声速冲击射流的流场结构，射流喷出后在中心线附近形成一道激波，通常被称为马赫盘，马赫盘周围有斜激波。在马赫盘后的亚音速流与斜激波后的超音速流之间会形成剪切层。在马赫盘后损失的总压大于通过斜激波的流动损失的总压，这导致了在平板中心的压力低于外侧，从而形成一个再循环区。在再循环区和马赫盘后的气流之间形成一个接触面。通过马赫盘的流体进入剪切层，远离射流的轴线。斜激波与射流边界的交点处出现扇形膨胀波，它经剪切层反射，形成环形的压缩波。这是欠膨胀超声速冲击射流可视化得到的平均流场主要特征。

关于噪声方面，超声速冲击射流会产生强烈的离散频率噪声，这种噪声被称为冲击单音，其幅值通常高于冲击射流中其他噪声，而成为冲击射流噪声的主体。Marsh[3]于1961年首先报道了射流冲击平板产生的离散频率噪声。Wanger[4]发展了冲击射流的不稳定性模型，Neuwerth[5]考虑到喷嘴处产生下行压力波动，对不稳定性模型进行了改进。Tam等[6]基于中性稳定波在喷流内部主要向上游运动的特点，发展了描述射流冲击平板的理论模型。Ho等[7]猜测离散频率的冲击单音与流场中的反馈环有关，是由平板附近的大尺度拟序结构导致。影响冲击射流流场结构的主要参数有喷管总压与环境压力之比（Nozzle Pressure Ratio，NPR）、喷口到冲击平板的距离以及冲击平板的大小。Alvi等[8]在实验中通过粒子图像测速技术（Particle Image Velocimetry ，PIV）得到了流场中的速度分布，并测试了冲击平板上的压力分布，发现压力脉动的主频与冲击单音的频率相同。Powell[9]指出，冲击单音是受喷嘴与壁面之间的反馈环控制。然而，由于流场中存在震荡，并且流-声关系十分复杂，这种反馈控制噪声的机理仍然不清晰。

本文参考Henderson等[10]的实验工况与刘海等[11-12]的数值模拟，对欠膨胀冲击射流进行数值模拟分析，着重研究喷嘴到壁面的反馈机制，考察湍流场中大尺度拟序结构的产生和发展演化过程以及冲击射流中的周期性现象。

1 计算对象与模型

欠膨胀冲击射流的大涡模拟（LES）使用经过Favre滤波后的Navier-Stokes方程作为控制方程。为使方程中的亚网格应力项封闭，本文使用了WALE亚网格模型[13]，其亚网格黏度被模化为：

式中：Ls=min(κd,CwV1/3) ；κ为von Kármán常数；d为网格点到壁面的距离；Cw=0.325，为WALE常数；V为控制体的体积。

本文计算参考Henderson的试验工况，其中收缩喷嘴的出口直径Dj=25.4 mm，喷口距壁面的高度与喷口直径之比为h/Dj=2.08，来流总压与环境压强的比值ptotal/p0=4.03 ，其完全膨胀马赫数为1.55，喷管外廓锥角为30°，唇厚度为1.27 mm。该工况对应的完全膨胀马赫数为1.55。计算域的设置如图1所示。计算网格数为228万，在冲击平板上壁面量纲为一距离y+＜5。对于喷管入口给定总压边界条件，对于远场出口给定无反射边界条件，冲击壁面和喷嘴壁面使用绝热无滑移边界条件。

图1 计算域设置Fig. 1 Computational domain

使用Fluent UDF在喷口前Dj处添加周向多模态涡环扰动，以触发流动转捩[14-15]。

计算中使用AUSM（Advection Upstream Splitting Method）通量分裂，空间上使用3阶MUSCL（Monotone Upstream-Centered Schemes for Conservation Laws）格式，时间上使用2阶精度隐式格式。初场由RANS Realizablek-ε湍流模型计算得出。大涡模拟的时间步长∆t1=1×10−6s ，可保证全场库朗数CFL=u＜1，Δx为网格间距。经过时长为400Dj/Uj的发展，认为已经排除了初始流场的影响，开始进行变量的采样统计。采样间隔为Δt2=5×10−6s，采样时间持续5×10−3s。

图2 不同涡环扰动计算结果对比Fig. 2 Comparison of calculation results by different vortex ring perturbations

2 结果与分析

2.1 时均流场

为验证网格无关性，本文分别对网格数为124万（Coarse），228万（Mid），342万（Fine）的3组网格进行雷诺时均模拟。图3示出了各组网格中心线上压力分布，从计算结果中可以看出，在网格数228万与342万的算例中，计算结果在峰值处的偏差3.8%，可满足网格无关性的要求。

图3 各组网格中心线上压力分布Fig. 3 Pressure distribution on centerline of each mesh

图4 中心线上物理量分布与实验对比Fig. 4 Comparison of the value on centerline with experiment

图5 冲击射流激波的比较Fig. 5 Comparison of the shock between the calculation and experiment in impinging jet

图6示出了冲击区域时均速度场计算与Henderson实验的比较，图6(a)是本文计算的速度云图，可以看到在喷口与激波之间的膨胀区域中，流速不断增大，并且在激波与壁面间存在低速的再循环区。这种趋势与图6(b)中Henderson实验得到的时均速度场基本一致。说明本文的模拟结果可以有效地模拟超声速冲击射流的流场。

图6 冲击射流时均速度计算与实验的比较Fig. 6 Comparison of average velocity between the calculation and experiment in impinging jet

2.2 瞬态结果与反馈机制

图7示出了轴截面瞬时合速度分布，各幅图是等时间间隔的，展示了1个完整的流场震荡周期(T=2×10-4s)内流场的变化，可以看到流场中激波的位置与强弱存在周期性。图7(a)和7(b)反映了波系沿轴向向下游运动的过程。在马赫盘与斜激波向下游运动的过程中，环形激波随之生成并沿轴向向下游发展，进入壁面射流区。在激波下行过程中，喷口喷出的欠膨胀气体有了更充分的空间进行膨胀，可以达到更高的速度，这导致了穿过马赫盘与斜激波的流体速度差增大，使得剪切作用增强，再循环区的影响范围增大。当再循环区向下游运动时，剪切层和斜激波后的超声速流动也发生畸变，斜激波后的环形激波随再循环区向下游运动而消失，其环形激波生成和消失的周期与马赫盘的震荡周期相同。图7(c)和图7(d)示出了马赫盘向上游运动的过程中，斜激波随马赫盘向上游摆动，在斜激波后方的流体速度差减小，再循环区也变小。在外剪切层与壁面射流的交界处（近壁面处y=0.8Dj）的位置产生涡结构，并有波向上游传播。

图7 一个流动震荡周期内的瞬态速度分布（Δt/T=0.25）Fig. 7 Transient velocity distribution during a flow oscillation period (Δt/T=0.25)

Q判据被广泛用于涡的辨识，其定义式为

图8 瞬态流场的Q判据等值面Fig. 8 Iso-surface of Q-criterion in unsteady flow field

图9 一个流动震荡周期内的瞬态速度脉动分布Fig. 9 Transient velocity fluctuation distribution during a flow oscillation period

速度散度 ∇·U为连续性方程中的一项，可以反映流场中的膨胀与压缩现象，也可以描述声波的传播，图10示出了一个流动震荡周期内的瞬态速度散度分布。可以看出，射流一经流出喷口，就发生高强度的膨胀，当膨胀波从喷管唇部传播到轴线附近时与对称侧的膨胀波交汇会带来更加剧烈的膨胀，在轴线上膨胀程度最高处的静压值低于环境压力。为了与射流外界的环境压力平衡，会形成环形的斜激波，气体在通过斜激波时会发生压缩。但由于本工况欠膨胀程度过高，斜激波的角度超过了激波正常反射条件，会出现类似于马赫反射的情况，在轴线附近形成了正激波马赫盘，拦截激波与正激波在交点处产生反射激波。在初始时刻，靠近壁面y=0.8Dj的位置处向上游辐射出一道较强的波，波会在喷管唇部发生反射，反射波在剪切层附近向下游传播，与新生成的波发生干涉并继续传播至冲击平板，从而形成反馈环。

图10 一个流动震荡周期内的瞬态速度散度分布Fig. 10 Transient velocity divergence distribution during a flow oscillation period

为精确探究流场震荡与反馈回路的周期，使用带有Hanning窗的快速傅里叶变换对流场中的脉动压力进行处理，样本数N=1024，采样间隔与数值模拟采样间隔Δt相同。图11和图12分别为快速傅里叶变换得到的功率谱密度（PSD）在中心线和唇线上的分布，图中示出了流场中的脉动压力在fs=5212 Hz和倍频2fs=10425 Hz时存在强烈的离散频率。这与Henderson实验中测得的冲击单音频率（10253 Hz）相符合。在图11示出的中心线脉动压力频率分布和图12示出的唇线脉动压力频率分布（图中使用对数标尺）中可以看到，在fs以及其各倍频处都有明显的强烈离散频率，该频率对应的周期为T=1.92×10−4s，与前面观测得到的周期T=2×10−4s基本相同。由于中心线上存在超声速流动，在超声速流动中脉动量较小，因此中心线上的数据存在一段低谷，而唇线上各处都有较高的脉动水平。

图11 中心线脉动压力频率分布Fig. 11 Frequency distribution of pressure fluctuation on centerline

图12 唇线脉动压力频率分布Fig. 12 Frequency distribution of pressure fluctuation on lip line

图13、图14分别示出了对应频率为fs和2fs的脉动压力分布，图中示出了在正激波马赫盘与反射激波的后方有频率为fs脉动压力分布，其位于轴线上x=1.6Dj处。并且，在壁面射流区中的涡环结构也以fs为主频。在正激波马赫盘与反射激波的前后均有频率为2fs压力脉动分布，集中于轴线上x=1.26Dj与x=1.6Dj处，并且存在反馈环的区域中，也分布有2fs的压力脉动信号，这可以直接说明，反馈环的频率与冲击单音的频率是相关的。

图13 5 212 Hz脉动压力幅值Fig. 13 Amplitude of pressure fluctuation at 5 212 Hz

图14 10 405 Hz脉动压力幅值Fig. 14 Amplitude of pressure fluctuation at 10 405 Hz

2.3 本征正交分解分析

本征正交分解（Proper Orthogonal Decomposition，POD），也被称为Karhunen-Loeve展开或主成分分析方法（Principal Components Analysis），是一种基于矩阵论的数据统计分析方法，广泛应用于数据降维，流场分析等。POD的主要目的是为流场提供分析手段且保持其物理意义，同时着眼于使用动态系统的概念来预测流场，可以实现对复杂的非线性系统的线性降维处理。考虑脉动分量u′(x,ti) ，将其分解为时间系数an(t)与POD的空间模态 φni(x) ，即

本文使用Snapshot方法进行POD分解[17-18]，对仿真计算得到由M行空间点、N列采样时间构成的原始数据矩阵AM×N，定义相关矩阵RN×N：

求解相关矩阵RN×N的特征值问题

得到相关矩阵RN×N的特征值λi与特征向量 φi。若原始数据矩阵AM×N为脉动速度，则λn可以表示第n阶模态对流场中湍动能的贡献。

定义第n阶模态的能量贡献率为

图15 POD能量贡献率Fig. 15 Energy contribution rate of POD modes

图16 POD累计能量贡献率Fig. 16 Cumulated energy contribution rate of POD modes

图17示出了脉动速度的1、2、3、4、8阶速度模态（各模态均以该模态最大幅值进行归一化处理并记为 φ¯ ）。1阶模态是流场中能量占比最高的模态，可以反映喷流流场中主要的大尺度相干结构，具有较强的对称性。在1阶模态中，射流边界、马赫盘与斜激波后、再循环区以及壁面射流均具有较强的相关性。再循环区在斜激波后超声速流的剪切作用下表现为对称的环形涡结构，再循环区前有明显的分界面。超声速流与壁面射流交界处产生的涡结构沿径向发展，并将外部环境流体卷吸到壁面射流中。在1、2阶模态的壁面射流区域中捕捉到了沿径向的模态对。3阶模态主要反映了在射流边界，再循环区与壁面射流的交界存在较强相关性。其主要现象为沿径向射出的壁面射流导致在靠近壁面y=0.8Dj～1.6Dj处产生涡结构，卷吸流体回到再循环区与壁面射流的交界处。而4阶模态表现为斜激波后的流动直接汇入壁面射流，这与前文提到的在再循环区与壁面射流的交界处的涡结构具有周期性相符合。8阶模态捕捉到了在反馈环控制区域存在的相干结构，与图10示出的反馈环相对应。