本文计算的圆口喷管出口直径D_j=50 mm^[16-17]，基于喷口速度及远场声速的声学马赫数$Ma = {{{U_{\rm{j}}}} / {{a_\infty }}} = 0.75$，基于喷口直径的雷诺数$R{e_D} = {{\rho {U_{\rm{j}}}{D_{\rm{j}}}} / \mu } = 8.7 \times {10^5}$，计算流域示于图1中，坐标原点位于喷管出口的圆心处，计算域总长度为97.6D_j。其中喷管长7.6 D_j，渐扩段为50 D_j，缓冲区长度40 D_j，计算域入口宽度为20 D_j，出口宽度为40 D_j。采用结构化网格，对剪切层区域进行了网格加密，网格节点数为${\rm{3}}{\rm{.5}} \times {\rm{1}}{{\rm{0}}^{\rm{6}}}$，本文前期已经通过外推法评估网格质量，结果表明该套网格可以分辨出唇线上80%以上的湍流脉动动能，满足LES的质量指标要求^[18]。

图  1  计算域及测点示意图

Figure 1.  Computational domain and measuring point positions

计算选用Fluent基于密度的求解器，扩散项和对流项分别采用中心差分和高阶迎风格式，时间采用隐式推进。初场由雷诺时均模拟给出，后续LES计算采用Smagorinsky-Lilly亚网格模型，无量纲时间推进步长为0.005 D_j/U_j，对应全场CFL≤1。经过时长为300D_j/U_j的发展，认为渐扩段已排除初始流场的影响，之后开始进行湍流量的采样统计，统计采样时间持续150D_j/U_j。采用FW-H方程计算远声场，Mendez等^[19]学者提出该方法要求声源面需要包含绝大部分含能结构，但又不能过大而影响截止频率。基于湍流近场的涡量轴截面分布，取声源面流向范围25D_j，入口区域半径2D_j，出口区域半径4D_j。监测点位于r=30D_j和r=50D_j的两个同心圆面，圆心为喷管出口中心，与流向监测夹角θ取30°至150°，测点和声源面位置如图1所示，内部虚线为声源面。

本文计算了三种不同来流条件下的圆口喷流，边界条件均为：入口边界给定总温320.4 K和总压147059.9 Pa，出口为静压边界，喷管壁面采用绝热无滑移条件，其他边界均给定远场静压进口条件，环境压力p_∞=101300 Pa，环境温度T_∞=288 K，对应a_∞=340 m/s。其中工况nopert的流场中无扰动施加，是本文用作对照的基础工况，更多计算细节前期工作^[18]已给出，工况pert1中加入来自Bogey^[8]的涡环扰动，工况pert2中的为改进的多模态线性不稳定扰动。

工况pert1扰动计算式：

式(1)中${u'_x}$和${u'_r}$分别表示轴向和径向速度脉动，为涡环($u_x^{\rm{ring}}$、$u_r^{\rm{ring}}$)和周向模态的组合；α取0.007；${\varepsilon _i}$和${\phi _i}$为第i个周向模态下的随机幅值和相位，分别在−1~1和0~2π之间随机变化，每个时间步更新一次；周向模态取n=4，m=15；r₀为涡环半径，等于喷口半径0.5D_j；扰动施加的轴向位置x₀≈r₀；∆₀为特征尺寸，即涡速度在${\Delta _0} = \displaystyle\sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {r - {r_0}} \right)}^2}}$达到最大值1。

工况pert2的来流扰动即是本文提出的以平行喷流扰动^[11]为基础改进多模态线性不稳定扰动方案。Hu等^[11]基于线性不稳定性分析，得出了一种平行射流入口的单一模态扰动方法。为使其可应用于圆口喷流，本文考虑原扰动参数已经优化，因此首先对该扰动的流向分量保留；横向分量则改造为圆口喷流的径向。同时考虑多模态快速触发湍流的优势，本文将原扰动方法的单一模态修正为周向多模态组合，得到一种改进的多模态线性不稳定入口扰动，轴向和径向扰动速度的计算式如下：

式(3)和(4)中a₁=0.33, a₂=0.4, b=6, $\omega = {{2{U_{\rm{j}}}} / {{D_{\rm{j}}}}}$，t为时间，扰动施加位置为喷口处，其余参数均采用涡环扰动中的优化数值^[8]。采用Fluent UDF将两种速度扰动施加在计算域中相应位置。

为方便表述，本文做如下规定：u和$u'$表示轴向速度及其脉动，v和$v'$则表示径向速度及其脉动；$\langle {\rm{g}}·\rangle$表示平均值，下标t和θ分别表示时间和周向平均。图2和图3示出了三种工况轴向速度的流向和径向分布。可以看出三种工况的轴向速度预测结果与文献[17]中实验结果吻合良好，pert1的势流核长度L_c/D_j 为5.52，最接近实验值6.5，pert2次之，nopert差距最大，各计算工况与实验中势流核长度如表1所示。此外，径向分布中，两种扰动工况比无扰动结果也有明显改善。以上结果表明本文所采用的计算方法及边界条件合理有效。

图  2  中心线速度轴向分布

Figure 2.  Axial velocity along centerline

图  3  轴向速度径向分布.图中从左至右位置x/D_j分别为1, 2.5和5. 线形含义同图2

Figure 3.  Radial profiles of axial velocity at x/D_j=1, 2.5 and 5. Legends are the same as in Fig.2

剪切层速度半值宽${r_{{\rm{0}}{\rm{.5}}}}$是衡量剪切层发展的一个重要指标，其斜率反映剪切层的增长速率。从图4半值宽沿流向分布中可以看出，nopert中存在斜率明显不同的两个阶段，这种转变是层流剪切层中涡配对引起的转捩^[20]。加入扰动后，两扰动工况的剪切层增长速率在近喷口区域升高，而在势流核下游降低，使剪切层半值宽整体上呈现无突变的缓慢增长，表明转捩在喷口附近短距离内已经完成，符合自由湍流的发展机制，形成了接近真实喷流的湍流来流。

图  4  剪切层半值宽

Figure 4.  Jet half width

为了验证该现象，图5和图6示出了350D_j/U_j时刻的瞬态涡强的分布。从图5涡量的轴截面分布中可以看出，nopert的涡配对在流向3D_j~4D_j发生，图中用虚线框标识出了这一区域，该位置对应剪切层半值宽斜率的转变，标志着层流向湍流过渡的开始。图5(b)和图5(c)中的pert1和pert2工况并没有配对现象，在$x < 3{D_{\rm{j}}}$之前，喷流剪切层已剧烈混合，处于完全发展的湍流状态，与均衡发展的剪切层一致。图6涡量在横截面${x / {{D_{\rm{j}}}}} = 2$的分布中，也可观察到nopert的喷流仍处于势流状态，涡量呈现规则的环状分布；而两扰动工况的旋涡周向不稳定性增强，这是由于加入的速度扰动破坏了流动稳定性，使流动与周围流体之间的掺混加剧。这些现象表明加入的扰动有效触发了流动快速转捩，形成了符合实际的湍流喷流。

图  5  xy平面涡量分布

Figure 5.  Instantaneous vorticity magnitude in the planes xy

图  6  x/D_j=2截面涡量分布

Figure 6.  Instantaneous vorticity magnitudein the planes x/D_j=2

喷射流动中，相干结构主导动力学以及声学特性，为了研究扰动对涡结构演化过程的影响，采用Q判据等值面进行可视化。图7为三种工况下Q判据等值面以及近喷口流域的局部放大图，采用轴向速度进行着色。从图7(a) nopert工况中可以观察到环状涡形成、涡配对到失稳破碎的完整级串过程，在流向x/D_j=1至4范围内有轴对称的涡环。而对于图7(b)和7(c)中两扰动工况，涡环的对称特征已再不明显，分别从pert2出口和pert1出口下游$x \approx {D_{\rm{j}}}$破碎成小尺度结构。施加扰动的射流均识别出了更细致的涡结构，物理现象更符合湍流来流的实际，再次证明了施加的两种扰动的有效性。

图  7  tU_j/D_j=350时刻Q判据等值面

Figure 7.  Iso-surface of Q-criterion at non-dimensional time tU_j/D_j =350

图8(a)和8(b)分别给出了轴向速度脉动和径向速度脉动沿中心线分布。三种计算工况的速度脉动水平均高于实验结果，且峰值位置向上游移动，这与中心线上快速的速度衰减而导致的粘性耗散有关。在$x < 3{D_{\rm{j}}}$流域，pert1和pert2的速度脉动水平高于nopert，这是由于pert1、pert2和nopert的剪切层在该区域处于不同的流动状态，与2.2章节的现象一致，并且pert2喷流向湍流过渡早于pert1，这是由于扰动涡分别从不同的流向位置开始破碎。流动向下游发展，pert1和pert2的预测结果有不同程度的改善，均低于nopert，其中pert1的峰值水平最接近实验值。图9中脉动速度的径向分布也可以看出，虽然在喷口附近本文计算的脉动峰值低于实验测量值，但随着流动向下游进行，两种扰动工况的预测结果均有所改善，与测量的脉动水平相当。相对于与参考文献的LES结果，本文两种扰动的加入，使喷流湍流特性更接近实验值。

图  8  中心线上脉动速度分布

Figure 8.  Fluctuating velocity along jet centerline

图  9  脉动速度径向分布. 图中从左至右位置x/D_j分别为1, 2.5和5.线形含义同图8

Figure 9.  Radial profiles of fluctuating velocity at x/D_j=1, 2.5 and 5. Legends are the same as in Fig.8

轴向速度脉动和径向速度脉动沿唇线分布如图10(a)和10(b)所示。对比三种工况，在唇线分布上同样可以观察到，其中pert1工况峰值水平最低。在靠近喷口区域$x < 2{D_{\rm{j}}}$，pert1和pert2的量纲为一速度脉动增长至实验值0.04以上，且pert2早于pert1，而nopert在该区域发展缓慢，随后在$x > 2{D_{\rm{j}}}$开始快速增长，Zhu^[20]提出这是层流剪切层产生的自然现象。势流核末端及其下游区域，扰动工况的径向速度脉动低于无扰动工况，而轴向速度脉动区别较小。两条流向线上的预测结果表明，施加的扰动能够触发出逼近实际喷嘴出口的湍流水平，并且降低下游的脉动水平，改善流动发展的计算准确性。

图  10  唇线上脉动速度分布

Figure 10.  Fluctuating velocity along jet lipline

高频声源位于中心线速度开始衰减的位置，低频声源位于势流核末端下游^[21]，因此选取了唇线上x/D_j=2.5，5，10三个流向位置，分析速度脉动谱密度PSD。图11和12分别是唇线上轴向脉动和径向脉动的PSD图。在流向位置x/D_j=2.5，nopert工况中轴向和径向速度脉动的峰值无量纲频率为0.6左右，接近线性不稳定性理论中$f{D_{\rm{j}}}/{U_{\rm{j}}} = 0.68$，该频率下的入口剪切层处于最不稳定轴对称模态，代表着不稳定波的发展^[22]，造成与理论差距的原因是可能是非线性效应。随着流动向下游进行，脉动峰值逐渐向低频移动，这与喷射流动中两种不稳定模态有关，剪切层模态主导喷口区域，由剪切层涡卷曲的高频振荡引起。下游则受喷流优先模态控制，与低频的大尺度湍流结构的发展有关^[23]。加入扰动后，在x/D_j=2.5位置pert1和pert2的PSD值在$f{D_{\rm{j}}}/{U_{\rm{j}}} \approx 0.6$变化不大，但其他在频段上均升高，因此与来流无扰动对比，扰动工况下的频谱没有明显的峰值，全频段的脉动水平都有所增强，说明nopert中的峰值可能与层流流动中的涡配对有关，而两种扰动加入以后，快速转捩得到的流动接近自然湍流，因此呈现出明显的多尺度特征，频谱表现为宽频特性。在势流核末端附近x/D_j=5，三种工况下的PSD趋于一致；而下游x/D_j=10的流向位置，pert1和pert2的PSD在整个频段上都低于nopert，与前文2.3.1中脉动沿流向分布结果一致，快速转捩使得扰动工况的速度脉动峰值以及下游的脉动水平下降，频谱分析也进一步表明扰动能促使喷流在近喷口区域快速转捩以及在下游湍流发展。

图  11  唇线上轴向脉动频谱

Figure 11.  PSD of axial velocity fluctuation along jet lipline

图  12  唇线上径向脉动频谱

Figure 12.  PSD of radial velocity fluctuation along jet lipline

图13和图14分别示出了r=30D_j和r=50D_j两圆面的总声压级OASPL沿测角θ 分布。对比计算与实验结果，可以看出nopert的OASPL高于测量结果5~10 dB，加入扰动后，与实验吻合程度明显改善。表明无扰动工况涡配对导致转捩的附加噪声引起很大的误差，而采用扰动促成快速转捩，可以有效地消除这种误差。在r=30D_j的圆面上，pert1在θ=30°~110°的预测优于pert2，而在θ=110°~150°时，pert2更接近实验值。而在更远位置r=50D_j，两种扰动工况预测相当。因此两种扰动虽有细微差别，但都能够大幅度地改善无扰动对声场计算造成的误差。

图  13  r=30D_j测点总声压级

Figure 13.  Far-field OASPL at r=30D_j

图  14  r=50D_j测点总声压级.线型含义同图13

Figure 14.  Far-field OASPL at r=50D_j. Legends are the same as in Fig.13

为了进一步研究扰动对声场预测的影响，选取了θ=60°、100°和140°三个角度进行频谱特性分析。图15和图16分别示出了r=30D_j和r=50D_j两圆面上的声压级SPL频谱图。与实验结果对比发现，nopert的SPL在低频段（$f{D_{\rm{j}}}/{U_{\rm{j}}} < 1$）高于实验结果，而高频段（$f{D_{\rm{j}}}/{U_{\rm{j}}} > 1$）低于实验结果，由于低频大尺度湍流噪声是总声场主要贡献者，其中就包括了涡配对带来的附加噪声误差，因此在前面总声压级分布中，nopert的预测高于实验值。相对于nopert，随着扰动的加入，pert1和pert2的高频噪声有所升高，但低频噪声的降低有效地使整体噪声减弱，也证明了亚音速喷流中低频噪声占优势的地位。此外，pert1和pert2的高频SPL与实验仍有较大的视觉误差，可能是计算网格对最小尺度湍流结构的分辨仍然不是很充分而引起的。

图  15  r=30D_j测点总声压级

Figure 15.  Far-field acoustic spectra at r=30D_j

图  16  r=50D_j测点总声压级.线型含义同图15

Figure 16.  Far-field acoustic spectra at r=50 D_j. Legends are the same as in Fig.15

本文针对喷射流动大涡模拟湍流来流条件给定的难点问题，提出了一种改进的多模态线性不稳定性入口扰动方法，并与当前文献中著名的涡环扰动对比，得到以下结论：

(1)触发喷口流动的快速转捩，使之接近实际的湍流来流，是喷流大涡模拟的一个关键。本文改进的多模态线性不稳定性扰动和涡环扰动都成功地达到了这一效果。

(2)来流快速转捩的实现，有效地消除了无扰动时涡配对造成的声场误差，改善了湍流喷流远声场的预测结果。

(3)本文改进的多模态线性不稳定扰动方法，克服了涡环扰动湍流起始位置须人为移动的缺点。多模态的加入，更符合湍流来流的多尺度特性，体现了湍流频谱的宽频特征。