NUMERICAL STUDY ON THE GENERATION POSITION OF SCREECH TONE IN AXISYMMETRIC SUPERSONIC JET BASED ON MODAL DECOMPOSITION
-
摘要: 不完全膨胀超声速射流的势核中会产生准周期的激波栅格结构, 其与剪切层内拟序结构的相互作用会产生激波噪声. 啸声是主要向上游方向传播的、具有离散频率的高强度激波噪声, 其产生是受一种非线性的声反馈环机制驱动. 精确定位啸声的声源位置是定量理解啸声反馈环机制和发展准确的啸声预测模型的一个关键所在. 为了分析近场啸声, 本文采用高精度数值方法直接求解轴对称可压缩Navier-Stokes方程, 数值模拟了完全膨胀射流马赫数为1.10和1.15的圆形声速喷管欠膨胀超声速冷射流, 得到了A1和A2两种轴对称模态啸声. 通过傅里叶模态分解、本征模态分解和动态模态分解, 分析了射流时序压力场和速度场, 研究了啸声关联拟序流动结构的空间演化, 精确定位了轴对称模态啸声的声源位置. 研究表明: 啸声关联拟序流动结构存在饱和态区域, 啸声声波是在其饱和态区域产生并向外传播; 在本文所涉及的射流马赫数范围内, A1和A2两种轴对称模态啸声的有效声源位置分别是在第4和第3个激波栅格结构的尾缘.Abstract: For the imperfectly expanded supersonic jet, the quasi-periodic shock-cell structures in jet core interacts with the coherent structures in shear layer to generate shock-associated noise. Screech tone is the shock-associated noise component with discrete frequency and high intensity, and it propagates primarily toward upstream direction. Its generation is driven by a nonlinear acoustic feedback loop. The exact nature of screech-generation mechanism, including source positions, has remained an open question. Accurately locating the sound source position of screech tone is a key point to quantitatively understand the screech feedback loop mechanism and to develop exact screech prediction model. In this paper, numerical simulations of underexpanded supersonic cold jet issuing from a circular sonic nozzle are carried out through solving axisymmetric compressible Navier-Stokes equations directly, using fifth order finite difference weighted essentially non-oscillatory scheme and third order total variation diminishing Runge-Kutta scheme. The fully expanded jet Mach numbers are 1.10 and 1.15. The present numerical result is compared and in good agreement with the experimental result in the literature. The axisymmetric A1 mode and A2 mode screech tones are captured. The time sequential pressure field and velocity field of the jet are analysed through the Fourier mode decomposition, the proper orthogonal decomposition and the dynamic mode decomposition. The spatial evolution of screech-associated coherent flow structures are studied and the sound source positions of axisymmetric A1 mode and A2 mode screech tones are accurately located. The results show that each screech-associated coherent flow structure has its own saturation state region, where the screech waves generate and radiate outward. It is also found that the effective source positions of axisymmetric A1 mode and A2 mode screech tones are the trailing edges of the fourth and the third shock-cells respectively for the jet Mach numbers considered.
-
引 言
超声速射流在航空航天飞行器的推进系统中广泛存在和应用. 由于喷管构型的限制, 大多数超声速射流是不完全膨胀的, 其势核中的激波栅格结构与剪切层中的湍流拟序结构的相互作用会产生强烈的激波噪声. 啸声是激波噪声频谱上的离散频率声模态, 主要向上游传播, 声压级可高达160 dB. 高强度的啸声不仅会造成严重的环境噪声污染, 还会使飞行器控制舵面、尾喷管结构部件产生声疲劳、甚至断裂, 危及飞行安全, 是发展军用/民用超声速飞行器急需解决的重要问题之一.
Powell[1]首次发现射流啸声, 并将其产生机制解释为非线性的自激声反馈环机制. 自此先驱性工作开始, 普遍认为啸声反馈环包含4种关键机制[2-4]: 剪切层中的拟序结构的增长、激波栅格结构和拟序结构的相互作用、声反馈和喷管唇口附近的感受性过程. 此外, Powell[1]还发现啸声频率在随完全膨胀射流马赫数(Mj)变化时存在模态跳跃行为. 在圆射流中, Powell确认了A、B、C、D 4种啸声模态. Merle[5-6]发现A模态可分为A1和A2两个亚模态. Davies和Oldfield[7]和Powell等[8]发现圆射流中的啸声模态与射流的内禀不稳定性直接相关, 并且确认A1和A2模态是轴对称模态, B和D模态是摆动模态, C模态是螺旋模态. Ponton和Seiner[9]发现摆动模态是由两个同量级的反向螺旋模态叠加而成. Clem等[10]在完全膨胀射流马赫数为1.0 < Mj < 1.7的圆形收缩喷管实验还捕捉到了啸声B′模态、E 模态和F 模态. 对于低马赫数超声速射流, 轴对称模态是主导啸声模态[11-12]. Ponton等[13-14]在声速喷管欠膨胀射流实验中记录了啸声A0模态, Li和Gao[15]在轴对称URANS模拟中也捕捉到了该模态. Gao和Li[16]根据他们提出的新的啸声波长(频率)预测公式, 提出Ponton等[13-14]在实验中测量到的两种A0 模态可以归类为A1 和A2 模态, 给出了啸声A0 模态的一种新解释.
理解啸声的产生机制, 需要对啸声反馈环的各个环节进行深入研究. 作为环节之一, 射流剪切层内拟序结构与势核内激波栅格结构的相互作用涉及啸声声波的产生过程和产生位置. 为了解释激波/拟序结构相互作用如何产生啸声声波, 基于模型问题, 即剪切层与单个孤立的压缩波-膨胀波结构单元的相互作用, Manning[17]提出了“激波泄漏(shock leakage)”机制. Suzuki 和Lele[18]进一步发现局部涡量对激波起到屏障作用, 激波会在剪切层内旋涡之间的鞍点也就是局部涡量较弱的区域泄漏, 也就是说, 激波泄漏和关联啸声的产生只能发生在射流剪切层大尺度湍流结构充分发展并且具有明晰的鞍点的区域. 自此之后, 激波泄漏机制受到了越来越多的关注, 并逐渐被数值[19]和实验[20-21]所证实.
激波泄漏机制虽然指出了剪切层中的旋涡鞍点是影响射流啸声产生关键结构, 但其还有不足: 未能给出每种模态啸声的精确产生位置. 为了定量的理解啸声反馈环, 精准定位每种模态啸声的有效声源位置至关重要. Powell等[8]基于纹影发现啸声A1、A2和B模态的声波是在喷管下游5倍直径处的区域辐射, 而啸声C模态的声波辐射中心位于喷管下游六倍直径的区域. Umeda和Ishii[22]利用瞬时纹影发现螺旋C模态的主控声源位于第3个激波栅格结构的尾缘. Edginton-Mitchell等[23]采用高分辨率平面PIV也研究了螺旋C模态, 认为啸声声波可能是在第2、第3和第4个激波反射点处, 且拟序涡量经历最大脉动时产生. Gao和Li[16]基于数值分析认为前5个激波栅格是啸声A1, A2, B和C模态的有效声源区, 并且主控声波是在第2和第4个激波栅格之间产生. Panda[24]则通过啸声同相压力脉动的近场映射发现啸声A2模态的声波是从第3和第4个激波尖梢之间某个位置辐射. Shariff和Manning[25]发现与A1模态相关的激波泄漏和噪声辐射是在第3和第4个激波栅格的激波尖梢处. Edgington-Mitchell 等[21]发现在第3和第4个激波栅格结构的激波尖梢处会产生与A1 模态相关的噪声辐射. 基于近场声测量和时间分辨的纹影成像, Mercier等[26-27]认为轴对称啸声A1 模态和A2 模态都是在第4个激波栅格的尾缘(激波尖梢)处产生, 而B 模态是在第3或第4个激波尖梢处产生, 具体取决于马赫数. Li等[28-29]通过数值模拟研究发现啸声螺旋C模态的有效声源位于第4个激波栅格, A1和A2模态的有效声源位于第4个激波栅格下游, B模态的有效声源位于第3个激波栅格下游.
精确定位啸声的声源位置是定量理解啸声反馈环和预测啸声频率的一个关键所在[26]. 然而, 截止目前, 关于每种啸声模态究竟是单声源还是多声源, 以及各模态啸声具体的声源位置还存在一定的争议. 而且, 同一类型啸声模态的产生位置还可能与具体的射流马赫数相关. 现有研究涉及的射流马赫数范围还比较有限, 需要针对不同的射流马赫数开展研究. 本文采用高精度数值方法直接求解轴对称可压缩Navier-Stokes方程, 数值模拟了射流马赫数为1.10和1.15的圆形声速喷管欠膨胀超声速冷射流, 得到了A1和A2两种轴对称模态啸声. 基于精细的非定常流场数据, 通过傅里叶模态分解、本征模态分解和动态模态分解, 研究了啸声关联流动结构的空间演化, 定位了两种轴对称模态啸声的有效声源位置.
1. 物理模型和数值模拟方法
1.1 控制方程
对于低射流马赫数欠膨胀超声速射流, 啸声的主导模态是轴对称模态, 在产生啸声的流动范围内, 剪切层内的大尺度不稳定波或拟序结构也呈轴对称形态, 三维效应并不明显. 因而, 控制方程采用轴对称可压缩Navier-Stokes (N-S)方程, 该方程是从柱坐标系下的三维可压缩N-S方程退化而来.
无量纲的轴对称可压缩N-S方程具体形式如下
$$ \frac{{\partial {\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}}}}{{\partial r}} = \frac{{\partial {{\boldsymbol{F}}_\nu }}}{{\partial x}} + \frac{{\partial {{\boldsymbol{G}}_\nu }}}{{\partial r}} + {\boldsymbol{S}} $$ (1) 其中
$$ {\boldsymbol{Q}} = {[\rho \;\;\;{\kern 1pt} \rho u\;\;\;{\kern 1pt} \rho v\;\;\;{\kern 1pt} E]^{\rm{T}}} $$ (2) 是守恒变量,
$ {\boldsymbol{F}} $ 和$ {\boldsymbol{G}} $ 是x和r方向的无黏通量$$ {\boldsymbol{F}} = \left[ {\begin{array}{*{20}{l}} {\rho u} \\ {\rho uu + p} \\ {\rho uv} \\ {\left( {E + p} \right)u} \end{array}} \right] , \;\; {\boldsymbol{G}} = \left[ {\begin{array}{*{20}{l}} {\rho v} \\ {\rho vu} \\ {\rho vv + p} \\ {\left( {E + p} \right)v} \end{array}} \right] $$ (3) $ {{\boldsymbol{F}}_v} $ 和$ {{\boldsymbol{G}}_v} $ 是x和r方向的黏性通量, 即$$ {{\boldsymbol{F}}_v} = \left[ {\begin{array}{*{20}{l}} 0 \\ {{\tau _{xx}}} \\ {{\tau _{xr}}} \\ {u{\tau _{xx}} + v{\tau _{xr}} - {{\dot q}_x}} \end{array}} \right] \text{, }\;\; {{\boldsymbol{G}}_v} = \left[ {\begin{array}{*{20}{l}} 0 \\ {{\tau _{rx}}} \\ {{\tau _{rr}}} \\ {u{\tau _{rx}} + v{\tau _{rr}} - {{\dot q}_r}} \end{array}} \right] $$ (4) $ {\boldsymbol{S}} $ 是源项, 即$$ {\boldsymbol{S}} = \frac{1}{r}\left( { - \left[ {\begin{array}{*{20}{l}} {\rho v} \\ {\rho vu} \\ {\rho vv} \\ {\left( {E + p} \right)v} \end{array}} \right] + \left[ {\begin{array}{*{20}{l}} 0 \\ {{\tau _{rx}}} \\ {{\tau _{rr}} - {\tau _{\theta \theta }}} \\ {u{\tau _{rx}} + v{\tau _{rr}} - {{\dot q}_r}} \end{array}} \right]} \right)$$ (5) E是单位体积总能, 即
$$ E = \frac{p}{{\gamma - 1}} + \frac{1}{2}\rho \left( {{u^2} + {v^2}} \right) $$ (6) 黏性应力项为
$$ \left.\begin{gathered} {\tau _{xx}} = \frac{2}{3}\frac{{\mu {M_a}}}{{R{e_D}}}\left( {2\frac{{\partial u}}{{\partial x}} - \frac{{\partial v}}{{\partial r}} - \frac{v}{r}} \right) \hfill \\ {\kern 1pt} {\tau _{rr}} = \frac{2}{3}\frac{{\mu {M_a}}}{{R{e_D}}}\left( {2\frac{{\partial v}}{{\partial r}} - \frac{{\partial u}}{{\partial x}} - \frac{v}{r}} \right) \hfill \\ {\tau _{xr}} = {\tau _{rx}} = \frac{{\mu {M_a}}}{{R{e_D}}}\left( {\frac{{\partial v}}{{\partial x}} + \frac{{\partial u}}{{\partial r}}} \right)\; \hfill \\ {\kern 1pt} {\tau _{\theta \theta }} = \frac{2}{3}\frac{{\mu {M_a}}}{{R{e_D}}}\left( {2\frac{v}{r} - \frac{{\partial u}}{{\partial x}} - \frac{{\partial v}}{{\partial r}}} \right) \end{gathered}\right\} $$ (7) x和r方向的热流率为
$$ \left.\begin{gathered} {{\dot q}_x} = - \frac{{\mu {M_a}}}{{R{e_D}}} \frac{1}{{\left( {\gamma - 1} \right)Pr}}\frac{{\partial T}}{{\partial x}} \hfill \\ {{\dot q}_r} = - \frac{{\mu {M_a}}}{{R{e_D}}} \frac{1}{{\left( {\gamma - 1} \right)Pr}}\frac{{\partial T}}{{\partial r}} \end{gathered} \right\}$$ (8) 完全气体状态方程和声速表达式如下
$$ \left.\begin{split} & p = \frac{1}{\gamma }\rho T \\ &a = \sqrt {\gamma \frac{p}{\rho }} \end{split}\right\}$$ (9) 黏性系数根据Sutherland公式计算
$$\left.\begin{array}{l} \mu = {T^{3/2}}\dfrac{{1 + C/T_\infty ^*}}{{T + C/T_\infty ^*}}\\ C = 110.4\;{\rm{K}} \end{array}\right\}$$ (10) 无量纲化控制方程时所涉及的参考变量如下: 长度、速度、密度、压力及总能、温度、黏性系数和黏性应力的参考量分别是
${D^*}$ (喷管出口直径),$ a_\infty ^ * $ (声速),$\rho _\infty ^*$ ,$\rho _\infty ^*a_\infty ^{*2}$ ,$T_\infty ^*$ ,$\mu _\infty ^*$ ,${{\mu _\infty ^*a_\infty ^*} \mathord{\left/ {\vphantom {{\mu _\infty ^*a_\infty ^*} {{D^*}}}} \right. } {{D^*}}}$ ; 这里的上标“*”和下标“$\infty $ ”表示的是有量纲的环境流动条件.假定射流是从设计马赫数为1的收缩喷管喷出, 喷管构型和流动参数参考自Li和Gao[15]的工作. 喷管唇口的厚度为
$0.2{D^ * }$ ,${D^*} = 25.4 \;{\rm{mm}}$ . 环境流动变量分别为: 温度${T^ * } = 288.15\; {\rm{K}}$ , 密度$\rho _\infty ^* = 1.225 $ $ {\rm{ kg/{m^3}}}$ , 声速$a_\infty ^* = 340.294 \;{\rm{m/s}}$ , 黏性系数$\mu _\infty ^* = 1.789\;38 $ $ {\rm{kg\;(m \cdot s)^{-1}}}$ . 对于空气, 普朗特数$Pr = 0.72$ . 完全膨胀射流马赫数、射流声速马赫数和基于喷管直径的参考雷诺数定义如下$$ \left.\begin{array}{l} {M_j} = {{U_j^ * } \mathord{\left/ {\vphantom {{U_j^ * } {a_j^*}}} \right. } {a_j^*}}\\ {\text{ }}{M_a} = {{U_j^ * } \mathord{\left/ {\vphantom {{U_j^ * } {a_\infty ^ * }}} \right. } {a_\infty ^ * }}\\ {\text{ }}R{e_D} = {{\rho _\infty ^ * U_j^ * {D^ * }} \mathord{\left/ {\vphantom {{\rho _\infty ^ * U_j^ * {D^ * }} {\mu _\infty ^ * }}} \right. } {\mu _\infty ^ * }} \end{array}\right\}$$ (11) 这里,
$U_j^ * $ 是有量纲的完全膨胀射流速度. 表1列出了本文数值模拟所涉及欠膨胀超声速射流的关键流动参数, 包括参考雷诺数ReD、完全膨胀射流马赫数Mj和相应的射流声速马赫数Ma.表 1 欠膨胀超声速射流的关键流动参数Table 1. The key flow parameters in underexpanded supersonic jetMj Ma ReD 1.10 0.987 5.481 × 105 1.15 1.023 6.052 × 105 1.2 数值方法
控制方程的空间对流项离散采用五阶精度的有限差分加权本质无振荡(WENO-JS-5)格式[30], 迎风处理采用Van Leer流通量分裂方法[31]; 此外, 还采用了改善WENO格式收敛能力的技术[32-33]来消除激波的波后振荡. 空间黏性项离散采用6阶精度的有限差分中心格式, 时间导数项离散采用3阶精度的3步TVD Runge-Kutta格式[30].
5阶精度的有限差分WENO-JS-5格式形式如下:
5点总模板上的数值通量可以表示为3个3点子模板上3阶数值通量的凸组合, 即
$$ \left.\begin{split} & {\hat f_{i + 1/2}} = \sum\limits_{k = 0}^2 {{\omega _k}} \hat f_{i + 1/2}^{(k)} \\ &\sum\limits_{k = 0}^2 {{\omega _k}} = 1 \end{split}\right\}$$ (12) 子模板上的数值通量为
$$ \left.\begin{gathered} \hat f_{i + 1/2}^{(0)} = \frac{1}{3}{f_{i - 2}} - \frac{7}{6}{f_{i - 1}} + \frac{{11}}{6}{f_i} \hfill \\ \hat f_{i + 1/2}^{(1)} = - \frac{1}{6}{f_{i - 1}} + \frac{5}{6}{f_i} + \frac{1}{3}{f_{i + 1}} \hfill \\ \hat f_{i + 1/2}^{(2)} = \frac{1}{3}{f_i} + \frac{5}{6}{f_{i + 1}} - \frac{1}{6}{f_{i + 2}} \end{gathered} \right\}$$ (13) 非线性权为
$$\left.\begin{split} & {\omega _k} = {{{\alpha _k}}}\Bigg/{{\sum\limits_{l = 0}^2 {{\alpha _l}}}} \\ &{\alpha }_{k}=\frac{{d}_{k}}{{(\varepsilon + I{S}_{k})}^{p}} \end{split}\right\}$$ (14) 其中,
$\varepsilon ={10}^{-6} $ 和$ p{\text{ = 2}} $ . 线性权为$$ {d_0} = \frac{1}{{10}} \text{, }\;\; {d_1} = \frac{6}{{10}} \text{, } \;\;{d_2} = \frac{3}{{10}} $$ (15) 光滑因子为
$$ \left.\begin{gathered} I{S_0} = \frac{{13}}{{12}}{({f_{i - 2}} - 2{f_{i - 1}} + {f_i})^2} + \frac{1}{4}{({f_{i - 2}} - 4{f_{i - 1}} + 3{f_i})^2} \hfill \\ I{S_1} = \frac{{13}}{{12}}{({f_{i - 1}} - 2{f_i} + {f_{i + 1}})^2} + \frac{1}{4}{({f_{i - 1}} - {f_{i + 1}})^2} \hfill \\ I{S_2} = \frac{{13}}{{12}}{({f_i} - 2{f_{i + 1}} + {f_{i + 2}})^2} + \frac{1}{4}{(3{f_i} - 4{f_{i + 1}} + {f_{i + 2}})^2} \end{gathered} \right\}$$ (16) 控制方程的对流项和黏性项经过空间离散后, 得到一组常微分方程
$$ \frac{{{\text{d}}{\boldsymbol{Q}}}}{{{\text{d}}t}} = L({\boldsymbol{Q}}) $$ (17) 其中,
$ L({\boldsymbol{Q}}) $ 表示空间离散算子. 3阶精度的3步TVD Runge-Kutta格式形式如下$$\left. \begin{gathered} {{\boldsymbol{Q}}^{(1)}}\; = {{\boldsymbol{Q}}^n} + \Delta t \cdot L({{\boldsymbol{Q}}^n}) \hfill \\ {{\boldsymbol{Q}}^{(2)}}\; = \frac{3}{4}{{\boldsymbol{Q}}^n} + \frac{1}{4}{{\boldsymbol{Q}}^{(1)}} + \frac{1}{4}\Delta t \cdot L({{\boldsymbol{Q}}^{(1)}}) \hfill \\ {{\boldsymbol{Q}}^{n + 1}} = \frac{1}{3}{{\boldsymbol{Q}}^n} + \frac{2}{3}{{\boldsymbol{Q}}^{(2)}} + \frac{2}{3}\Delta t \cdot L({{\boldsymbol{Q}}^{(2)}}) \end{gathered} \right\}$$ (18) 1.3 计算域和初边值条件
计算域和边界条件的设置如图1所示. 在物理域
$ - 2 D \leqslant x \leqslant 15 D$ 和$0 \leqslant r \leqslant 6 D$ 内, 计算网格均匀分布,$\Delta x = \Delta r = 0.02 D$ . 在物理域之外, 网格开始拉伸以消除反射声波.![]()
图 1 计算域和边界条件示意图Figure 1. The schematic of the computational domain and boundary conditions计算域的上边界和左边界采用Thompson的特征远场边界条件[34-35]; 在下游边界区域, 采用无反射出流边界条件[36]. 在射流轴线上, 由控制方程中含“
${1 \mathord{\left/ {\vphantom {1 r}} \right. } r}$ ”的项所导致的奇性轴问题采用极限法[37]转换成偏导数进行处理. 此外, 对轴线上的偏导数进行有限差分离散时, 计算域也需拓展到的$r < 0$ 的区域, 相应的变量通过与$r > 0$ 的区域的镜像关系求得. 喷管壁面采用无滑移边界条件. 在喷管出口处, 入流平面内移6个网格以便激发出不稳定波, 减小对反馈环的数值限制和影响. 喷嘴出口处的入流平面为声速平面, 无量纲形式的流动变量设置如下[15]$$ \left.\begin{array}{l} {p_e} = \dfrac{1}{\gamma }{\left[ {\dfrac{{2 + \left( {\gamma - 1} \right)M_j^2}}{{\gamma + 1}}} \right]^{\frac{\gamma }{{\gamma - 1}}}}\\{\text{ }}{\rho _e} = \dfrac{{\gamma \left( {\gamma + 1} \right)}}{{2{T_r}}}{p_e} \hfill \\ {u_e} = \sqrt {\dfrac{{2{T_r}}}{{\gamma + 1}}} \\ {v_e} = 0 \end{array}\right\} $$ (19) 其中, 比热比
$\gamma = 1.4$ ;${T_r}$ 为储气罐温度, 根据冷射流假设,${T_r} = 1$ .需要强调的是, 在整个计算过程中, 入流条件处并没有添加任何形式的人工扰动. 因此, 啸声反馈环的产生和维持都不需要任何外部激励.
初始时刻, 除了喷管出口处, 整个计算域内都设置为环境流动条件
$$ \rho = 1,{\text{ }}u = 0,{\text{ }}v = 0,{\text{ }}p = {1 \mathord{\left/ {\vphantom {1 \gamma }} \right. } \gamma },{\text{ }}T = 1 $$ (20) 上述参数表示射流由喷管射入周围静止的空气.
为了方便对数值结果进行分析, 计算选取固定的无量纲时间步长, 具体取值为
$\Delta t = 0.002.$ 为了确保时序信号傅里叶分析的精度, 计算需要推进大量的时间步, 总的无量纲积分时长为T=2500. 本文所采用的数值方法、计算程序及针对轴对称超声速射流计算的网格无关性验证见文献[38].2. 数值模拟结果与分析
2.1 激波栅格结构
啸声声波是由射流势核中的激波栅格结构与剪切层中的大尺度拟序结构相互作用产生, 因此, 准确地模拟激波栅格结构的位置对于定位啸声的有效声源位置至关重要.
图2是完全膨胀射流马赫数为1.10和1.15的欠膨胀超声速射流轴线上的时间平均压力分布. 为了验证当前轴对称数值模拟的准确性, 将其结果与André等[39-40]的实验结果进行了对比. 图中显示, 对于前3个激波栅格, 两种状态的数值解都与实验测量值吻合得很好, 无论是激波位置还是幅值. 对于第4个激波栅格, 数值解的激波幅值与实验测量值有所偏差, 但激波位置吻合. 在第5个激波栅格结构之后的区域, 数值解与实验测量值差别较大, 其原因可能是在此区域轴对称控制方程物理上不再成立. 由前文可知, 啸声一般是在第5个激波栅格的上游区域产生, 所以这些偏差并不影响后续关于轴对称啸声模态产生位置的分析.
![]()
![]()
2.2 啸声的频率与声压级
在喷管唇口壁面上设置压力监测点[0.0D, 0.642D], 记录该位置处的时序压力信号, 采用快速傅里叶变换(单边模式)进行处理, 得到近场噪声频谱信息. 图3是Mj=1.10, 1.15这两种状态欠膨胀超声速射流的近场噪声频谱分析结果. 图中的横坐标表示频率, 单位是 Hz, 频带宽度是75 Hz; 纵坐标表示声压级, 单位是dB. 对于Mj=1.10的欠膨胀超声速射流, 有两个尖锐的声压级凸起, 分别是轴对称A1模态及其谐频. 对于Mj=1.15的欠膨胀超声速射流, 可以识别出3个有意义的声压级凸起, 分别是轴对称A2模态及其第1、第2谐频.
![]()
图 3 射流近场噪声频谱信息: 喷管唇口壁面监测点[0.0D, 0.642D]Figure 3. Spectral information of jet near-field noise: monitor [0.0D, 0.642D] located on the nozzle exit lip wall表2具体给出了上述离散频率模态的频率和声压级, 为了进行对比, 同时列出了Ponton等[13-14]的实验测量值和Loh等[41]的数值模拟结果. 对于Mj=1.10的射流, 当前数值结果与Loh等[41]的数值结果都捕捉到了轴对称A1模态及其谐频, 且频率和声压级吻合得很好. 对于Mj=1.15的射流, 不论是模态类型, 还是模态频率和声压级, 当前数值结果都与实验结果[13-14]要更为吻合.
表 2 啸声模态及其谐频的频率和声压级Table 2. The frequencies and amplitudes of screech modes and their harmonicsScreech modes Numerical results Experimental results Present Ref. [41] Ref. [13−14] Mj 1.10 1.15 1.10 1.15 1.15 A1 Freq/Hz 8560 − 9040 7690 − SPL/dB 138.4 − 139 137 − A2 Freq/Hz − 9250 − 9410 8682 SPL/dB − 135.5 − 125 147 A1
harmonic 1stFreq/Hz 17070 − 18073 15337 − SPL/dB 117 − 109 122 − A2
harmonic 1stFreq/Hz − 18526 − − 17332 SPL/dB − 115 − − 116 A2
harmonic 2ndFreq/Hz − 27879 − − − SPL/dB − 111 − − − 2.3 轴对称A1模态啸声的声源位置
图4是在完全膨胀射流马赫数Mj=1.10的欠膨胀射流中与A1模态相关的瞬时流场结构和声场. 从图中可以看到3道向上游传播的啸声声波. 采用Edgington-Mitchell等[21]的方法, 在图4中叠加一个与啸声声波重合的圆(蓝色虚线), 通过圆心来定位声源位置. 图4显示, 圆心位于流向位置x=2.29D处, 该位置恰好是第4个激波栅格结构的尾缘, 也称激波反射点或激波尖梢; 根据纹影和涡量还可以确定, 在此位置处, 剪切层内同时出现了一个旋涡鞍点. 激波泄漏机制[17-18]也指出激波会在剪切层鞍点处向外泄漏, 形成啸声声波. 因此, 可以确定第4个激波栅格结构的尾缘是轴对称啸声A1模态的有效声源位置.
![]()
图 4 Mj=1.10的欠膨胀超声速射流中与啸声A1模态相关的瞬时流动结构和声场(数值纹影$\left| {\nabla \rho } \right|$: 橙红色连续云图, 取值范围[0, 2]; 胀量$\nabla \cdot {\boldsymbol{V}} $: 灰色连续云图, 取值范围[−0.02, 0.02], 白色正值表示膨胀, 黑色负值代表压缩; 涡量: 彩色等值线, 取值范围[−0.5, 20], 15级)Figure 4. The instantaneous flow structures and acoustic field associated with screech A1 mode in the underexpanded supersonic jet with Mj=1.10(Numerical schlieren$\left| {\nabla \rho } \right|$: orange red continuous contour from 0 to 2; dilatation$\nabla \cdot {\boldsymbol{V}}$: gray continuous contour from −0.02 to 0.02, positive value with white color means expansion, negative value with black color means compression; vorticity: multi-color banded contour from −0.5 to 20, 15 levels)2.4 轴对称A2模态啸声的声源位置
通过快速傅里叶变换, 对射流压力场进行傅里叶模态分解, 即将时域空间的脉动压力场时间序列变换到频域空间, 提取啸声频率成分, 逐点计算声压级和相位, 得到啸声频率上声压级和相位的空间分布. 声压级最大的位置可认为是啸声模态的产生位置. 本节的压力场傅里叶模态分解共选取了8192个时刻的数据, 起始时刻为t = 200, 无量纲时间间隔为
$\delta t = 0.1$ .图5是Mj = 1.15的欠膨胀超声速射流中啸声A2模态频率上的声压级空间分布和势核区的时间平均压力空间分布. 在图5(a)声压级云图中可以看到沿射流剪切层外缘分布的啸声A2模态关联驻波结构. 驻波结构是向下游传播的剪切层不稳定波和向上游传播的声波相互作用所产生[24]. 图5(b)中的时间平均压力分布展现的是射流势核内的激波栅格结构. 结合图5(a)和图5(b)可知, 声压级最大的位置是在流向位置x = 2.14D处, 即第3个激波栅格结构尾缘. 因此, 可以确定啸声A2模态的声源位于第3个激波栅格结构的尾缘.
图6是Mj = 1.15的欠膨胀超声速射流中啸声A2模态频率上的相位空间分布, 从图中可以清晰地看到声波传播的波振面. 为了进一步确定A2模态的声源位置, 图7提取了直线R/D = 1.5上的相位分布. 相位斜率的符号表示声波的传播方向, 而相位斜率符号变化的位置也可指示声源位置[26,28-29]. 图中显示, 相位斜率的符号在流向位置x = 2.14D附近发生了明显变化, 在其上游, 斜率为正值, 声波向上游传播, 而在其下游, 斜率为负值, 声波则向下游传播; 该位置同时也是图5(a)中声压级最大的位置, 表明啸声A2模态的声源位置的确是在第3个激波栅格结构的尾缘.
![]()
图 6 Mj=115的欠膨胀超声速射流中啸声A2模态频率上的相位空间分布Figure 6. The spatial distribution of phase at the frequency of screech A2 mode in underexpanded supersonic jet with Mj=1.15![]()
图 7 Mj=1.15的欠膨胀超声速射流中啸声A2模态频率上的相位沿直线R/D=1.5的分布Figure 7. The distribution of phase at the frequency of screech A2 mode along the line of R/D=1.5 in underexpanded supersonic jet with Mj=1.152.5 基于本征模态分解的声源位置分析
本征正交分解(proper orthogonal decomposition, POD)是提取、分析湍流中拟序流动结构的重要工具[42]. POD方法的实质是在最小二乘意义下提供能够代表已知数据的一组正交基. POD快照法[43]的目标是寻找脉动速度场数据快照
${\boldsymbol{u}}'$ 的最小二乘表达式, 即残差$ {{\boldsymbol{u}}_{res}} $ 的模最小; 相应的数学表达式[23,44]如下$$ {\boldsymbol{u}}'\left( {{\boldsymbol{x}},{t_k}} \right) = \sum\limits_{i = 1}^M {{{\boldsymbol{b}}_i}\left( {{t_k}} \right){{\boldsymbol{\psi}} _i}\left({\boldsymbol{ x}} \right)} + {{\boldsymbol{u}}_{res}}\left( {{\boldsymbol{x}},t} \right) $$ (21) 其中,
$ {{\boldsymbol{b}}_i} $ 是时间POD模态,$ {{\boldsymbol{\psi }}_i} $ 是空间POD模态. M是模态数, M≤N, N是数据快照序列的快照总数.为了得到POD模态, 首先构造脉动速度场的快照矩阵
${\boldsymbol{V}}$ , 随后计算脉动速度场的自相关(对称正定)矩阵$ {\boldsymbol{R = }}{{\boldsymbol{V}}^{\text{T}}}{\boldsymbol{V}} $ .针对自相关矩阵进行特征分解
${\boldsymbol{R}}{{\boldsymbol{a}}_i} = {\lambda _i}{{\boldsymbol{a}}_i}$ , 求得自相关矩阵的特征值${\lambda _i}$ 和特征向量${{\boldsymbol{a}}_i}$ $$ {\boldsymbol{V}} = \left[ {{{{\boldsymbol{u}}'}^1}{\text{ }}{{{\boldsymbol{u}}'}^2} \ldots {\text{ }}{{{\boldsymbol{u}}'}^N}} \right] = \left[ {\begin{array}{*{20}{c}} {{u'}_{11}^1}&{{u'}_{11}^2}& \cdots &{{u'}_{11}^N} \\ \vdots & \vdots & \vdots & \vdots \\ {{u'}_{IJ}^1}&{{u'}_{IJ}^2}& \cdots &{{u'}_{IJ}^N} \\ {{v'}_{11}^1}&{{v'}_{11}^2}& \cdots &{{v'}_{11}^N} \\ \vdots & \vdots & \vdots & \vdots \\ {{v'}_{IJ}^1}&{{v'}_{IJ}^2}& \cdots &{{v'}_{IJ}^N} \end{array}} \right] $$ (22) 特征向量
${{\boldsymbol{a}}_i}: = \left[ {{a_i}\left( {{t_1}} \right),{{a_i}\left( {{t_2}} \right), \cdots ,{a_i}\left( {{t_k}} \right)}} \right]$ 表征了正交的时间POD模态$$ {{\boldsymbol{b}}_i} = \sqrt {{\lambda _i}} {{\boldsymbol{a}}_i} $$ (23) 空间POD模态可以表示为脉动速度场快照与时间POD模态乘积的线性组合
$$ {{\boldsymbol{\psi}} _i}\left( x \right) = \frac{1}{{\sqrt {{\lambda _i}} }}\mathop \sum \limits_{k = 1}^N {a_i}\left( {{t_k}} \right){{\boldsymbol{u}}'_k}\left( {{\boldsymbol{x}},{t_k}} \right) $$ (24) 式中特征值
${\lambda _i}$ 表征了每个空间POD模态两倍的不可压脉动动能. 因此, POD方法能够提取最富能的流动结构. 特别地, 为了表征每个POD模态的能量占比, 需要对特征值${\lambda _i}$ 进行归一化处理, 并按大小进行排序. 需要注意的是, POD模态通常是成对出现的.本节的POD分析选取了共2048个时刻的速度场数据快照, 无量纲时间间隔
$ \delta t = 0.1 $ .图8是Mj=1.10的欠膨胀超声速射流速度场的时间POD模态, 包括波形、相图和特征值. POD模态都是成对出现的, 图8(a)和图8(b)是1阶和2阶POD时间模态[
${a_1}({t_k})$ 和${a_2}({t_k})$ , k是快照时刻]及其相图. 相图能够表征各阶POD模态间的时间相关性,${r_m} = \overline {\sqrt {a_1^2 + a_2^2} } $ 表示两个模态的平均幅值. 1阶和2阶模态的波长是相匹配的; 相图中的点大致呈圆形散布, 表明前两阶模态是时间强相关的. 时间相关性以及匹配的波长说明1阶模态和2阶模态是耦合的, 两者共同表征了1个振荡过程或拟序结构. 图8(c)是各阶POD模态的能量占比, 1阶和2阶模态含能最高, 两者之和接近45%, 大于其它各阶模态.![]()
图 8 脉动速度场的POD时间模态和特征值, Mj=1.10Figure 8. Temporal POD modes and their eigenvalues of fluctuating velocity field, Mj = 1.10图8(d)是1阶和3阶POD模态的相图. 图中的点呈双扭线形分布, 表明1阶、3阶模态之间没有类似1阶、2阶模态的时间相关性. 图8(e)和图8(f)是5阶和6阶POD时间模态[
${a_5}({t_k})$ 和${a_6}({t_k})$ ]及其相图,${r_m} = \overline {\sqrt {a_5^2 + a_6^2} } $ 同样表示两者的平均幅值. 与前两阶模态类似, 5阶和6阶模态的波长相匹配, 并且是时间强相关的, 两者也共同表征了1种拟序结构.图9是Mj=1.10的欠膨胀超声速射流速度场第1对和第3对POD时间模态的傅里叶分析. 一般来说, 每个POD模态包含多种频率组分. 图9显示, 这两对POD模态都有一个很强的主导频率, 其中,第1对的主导频率是4396 Hz, 第3对的主导频率是8595 Hz. 第3对POD模态的主导频率与啸声A1模态的频率(8560 Hz)十分一致, 即第3对POD模态表征的是与啸声A1模态相关的拟序流动结构. 第1对POD模态的主导频率接近第3对POD模态的1/2, 因此, 其所表征的拟序流动结构是啸声A1模态关联拟序流动结构的亚谐频.
![]()
图 9 脉动速度场的POD时间模态的傅里叶分析, Mj=1.10Figure 9. The FFT analysis of the first and the third pairs of temporal POD modes of fluctuating velocity field, Mj = 1.10图10是Mj=1.10的欠膨胀超声速射流速度场的第1对和第3对POD空间模态, 包括u分量和v分量. 图中显示, 空间POD模态的u分量存在着非常明显的交界面, 在交界面两侧符号相反, 这一交界面是剪切层与射流主流(势核)的交界面; 此外, 每1对空间POD模态之间都存在着1/4波长的相位差. 第1对POD模态所表征的拟序结构在流向位置
$x > 3.8 D$ 的区域达到饱和态; 第3对POD模态所表征的啸声相关拟序结构发展的更快, 在流向位置$1.5 D < x < 3.8 D$ 的区域就达到饱和态, 随后开始衰减. 前文2.3节所述啸声A1模态的声源位置($x = 2.29 D$ )就位于该饱和态区域内. 在第3对POD模态中, 还可以清楚地看到从饱和态区域发出的, 规律性地向上游和下游辐射的条纹, 而向上游传播的条纹即为A1模态的啸声. 根据主导频率和空间演化, 第1对POD模态所表征的拟序结构是由啸声A1模态关联拟序结构在下游发生旋涡合并所产生.![]()
图 10 脉动速度场的POD空间模态, Mj=1.10Figure 10. The first and the third pairs of spatial POD modes of fluctuating velocity field, Mj = 1.10通过对不同流向位置处的脉动速度场切片进行POD分析, 可以研究拟序结构能量占比沿流向的发展变化. 在每一流向站位, 提取特定时间间隔内的切片速度场时间序列组成局部数据快照矩阵, 然后对局部数据快照矩阵进行POD分析.
图11是Mj=1.10的欠膨胀超声速射流速度场POD模态特征值沿流向的变化曲线. 图中显示, 在流向位置
$x = 3.1 D$ 处, 第2对POD模态特征值之和的极大值点恰好对应第1对POD模态特征值之和在此处的极小值点; 第1对POD模态所表征的拟序结构在此位置前经历了较为强烈的能量释放, 而第2对POD模态所表征的拟序结构在此位置前一直在吸收能量, 表明这两种拟序结构之间存在着能量传递. 第3对POD模态特征值之和沿流向的变化表明, 啸声A1模态关联拟序结构在其饱和态区域(1.5D$ < $ $ x < 3.8 D$ )内, 先后经历了能量吸收再释放的过程, 并且在能量释放子区间内, 其能量占比, 即第3对POD模态的特征值之和, 在流向位置$x = 2.29 D$ 处产生了一个局部极小值点, 这个局部极小值点就是前文2.3节中确定的啸声A1模态的产生位置.![]()
图 11 脉动速度场轴向切片POD分析中相关模态的特征值沿流向的变化, Mj = 1.10Figure 11. The variations of eigenvalues along streamwise of the POD modes based on the POD analysis of velocity field’s slices in different axial position, Mj = 1.102.6 基于动态模态分解的声源位置分析
动态模态分解[45](dynamic mode decomposition, DMD)方法则是基于动力系统Koopman分析发展而来的分析工具, 其通过时间正交的特征模态来重构数据快照序列. POD模态强制空间正交性(结构去相关), 每个独立的POD模态都包含了多种频率成分, 而DMD模态是时间正交的(单一频率), 一般空间不正交. 如果数据序列足够长, DMD方法能够识别主频, 提取相关的流动特征. 对源自于饱和非线性过程的实验或数值数据来说, 相较于提取按能量排序的空间去相关结构, 提取带有特定频率信息的时间去相关结构更为重要.
基于DMD模态的数据快照重构表达式[46]为
$$ \begin{gathered} {{\boldsymbol{x}}_{{\rm{DMD}}}}\left( t \right) = \sum\limits_{k = 1}^K {{b_k}\left( 0 \right){{\boldsymbol{\phi}} _k}\left( {\boldsymbol{x}} \right){\rm{exp}}\left( {{\omega _k}t} \right)} =\hfill \\ {\text{ }} {\boldsymbol{\varPhi}} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\rm{diag}}\left\{ {{\rm{exp}}\left( {{\omega _1}t} \right),{\rm{exp}}\left( {{\omega _2}t} \right), \cdots ,{\rm{exp}}\left( {{\omega _K}t} \right)} \right\}{\boldsymbol{b}} \hfill \\ \end{gathered} $$ (25) 式中,
$ {{\boldsymbol{\phi}} _k} $ 是DMD模态, K是DMD模态的数目,$ {\omega _k} = {{\ln{{\mu _k}} } \mathord{\left/ {\vphantom {{\ln {{\mu _k}} } {\Delta t}}} \right. } {\Delta t}} $ 的虚部是DMD模态的圆频率, 而实部则表示指数增长或衰减, 具体取决于符号;$ {b_k}(0) $ 是每个DMD模态的初始幅值. 初始时刻的数据快照$ {{\boldsymbol{x}}_1} = {\boldsymbol{\varPhi}} {\text{ }}{\boldsymbol{b}} $ , 所以DMD模态的初始幅值可通过${\boldsymbol{b}} = {{\boldsymbol{\varPhi}} ^ + }{{\boldsymbol{x}}_1}$ 求得.本节针对压力场的DMD分析选取与2.5节POD分析中相同时刻的数据快照序列. 具体求解过程如下.
给定数据快照矩阵
$ {{\boldsymbol{V}}}_1^{N + 1} $ $$ {{\boldsymbol{V}}}_1^N = \{ {{{\boldsymbol{x}}}_1},{\text{ }}{{{\boldsymbol{x}}}_2},{\text{ }} \cdots ,{\text{ }}{{{\boldsymbol{x}}}_N},{\text{ }}{{{\boldsymbol{x}}}_{N + 1}}\} $$ (26) 令矩阵
$ {{\boldsymbol{S}}} = {\left( {{{\boldsymbol{V}}}_1^N} \right)^ + }{{\boldsymbol{V}}}_2^{N + 1} $ , 对其进行特征分解$$ {\boldsymbol{S}}{{\boldsymbol{y}}_j} = {\mu _j}{{\boldsymbol{y}}_j} $$ (27) 其中,
${\mu _j}$ 和${{{\boldsymbol{y}}}_j}$ 分别是特征值和特征向量.对数据快照矩阵
${{\boldsymbol{V}}}_1^N$ 进行奇异值分解$$ {{\boldsymbol{V}}}_1^N = {{\boldsymbol{U}}}{\boldsymbol{\varSigma }}{{{\boldsymbol{W}}}^*} $$ (28) 其中,
${{\boldsymbol{U}}}$ 和$ {{\boldsymbol{W}}} $ 是酉矩阵,${\boldsymbol{\varSigma }}$ 是奇异值对角矩阵.最终可得DMD模态
$$ \left.\begin{split} & {{\boldsymbol{\phi}} _j} = {{\boldsymbol{U}}}{{{\boldsymbol{y}}}_j}\\ &{\boldsymbol{\varPhi}} = \left( {{{\boldsymbol{\phi }}_1},{\text{ }}{{\boldsymbol{\phi }}_2}, \cdots ,{\text{ }}{{\boldsymbol{\phi }}_N}} \right) \end{split} \right\}$$ (29) 需要说明的是, 本节后续分析中的DMD模态都是包括了初始幅值的DMD模态, 即
$ {b_j}{{\boldsymbol{\phi }}_j} $ .图12是Mj = 1.10的欠膨胀超声速射流的压力场DMD模态. 各阶DMD模态的频率分布见图12(a), 其中, 幅值最高的前4阶DMD模态的频率依次是4383 Hz, 8570 Hz, 4252 Hz和4186 Hz. 对比2.5节中Mj = 1.10的欠膨胀射流速度场POD分析结果可知, 压力场的1阶、3阶和4阶DMD模态的频率都十分接近速度场第1对POD模态的主导频率(4396 Hz), 并且三者具有相似的空间形态, 这里只给出了3阶DMD模态, 见图12(b), 其所表征的拟序结构也是在流向位置
$x > 3.8 D$ 达到饱和态. 2阶DMD模态的频率接近第3对速度场POD模态的主导频率(8595 Hz), 此频率也是啸声A1模态的频率(8560 Hz), 因此, 2阶DMD模态表征的是啸声A1模态的关联拟序结构, 其饱和态区域为$1.5 D < x < 3.8 D$ , 见图12(c), 从中可以观察到从饱和态区域向外辐射的啸声A1模态声波条纹. 根据频率, 3阶DMD模态表征的拟序结构是2阶DMD模态所表征拟序结构的亚谐频; 而且在空间演化上, 两者的饱和态区域是接续的, 因此, 前者是由后者在下游发生旋涡合并后所产生.![]()
图 12 Mj=1.10的欠膨胀超声速射流的脉动压力场DMD模态Figure 12. The DMD modes of fluctuating pressure field in the underexpanded supersonic jet with Mj = 1.10图12(d)是Mj = 1.10的欠膨胀超声速射流压力场2阶DMD模态的相位空间分布, 清晰地表征了A1模态啸声声波传播的波振面. 图13提取了2阶DMD模态在直线R/D = 2.1上的相位分布. 前文2.4节中已经说明相位斜率符号改变的位置可以指示啸声的声源位置. 2阶DMD模态相位斜率的符号是在流向位置
$x = 2.29 D$ 处发生改变, 该位置也是前文2.3节和2.5节中所确定的啸声A1模态的产生位置, 这再次证实啸声A1模态的有效声源位置是在第四个激波栅格结构的尾缘.![]()
图 13 2阶DMD模态的相位在直线(R/D = 2.1)上的分布, Mj=1.10Figure 13. The distribution of phase along the line R/D = 2.1 for the second DMD mode, Mj=1.10图14是Mj = 1.15的欠膨胀超声速射流的压力场DMD模态. 图14(a)是各DMD模态的频率分布, 其中, 幅值最高的前4阶DMD模态的频率依次是4579 Hz, 4644 Hz, 4710 Hz和9289 Hz. 前3阶DMD模态的频率近似相等, 而且也具有相似的空间形态, 因此, 这里只给出了2阶模态, 见图14(b). 4阶DMD模态接近啸声A2模态的频率(9250 Hz), 因此, 它表征的是啸声A2模态关联拟序结构, 图14(c)显示其饱和态区域为
$1.65 < x < 3.25$ , 并且可以观察到自此区域向外辐射的啸声A2模态声波条纹. 2阶DMD模态是4阶DMD模态的亚谐频, 其所表征的拟序结构在流向位置$x > 3.25$ 处达到饱和态. 同样地, 根据频率和空间演化可知, 前者也是由后者在下游发生旋涡合并所产生.![]()
图 14 Mj=1.15的欠膨胀超声速射流的脉动压力场DMD模态Figure 14. The DMD modes of fluctuating pressure field in the underexpanded supersonic jet with Mj = 1.15图14(d)是Mj = 1.15的欠膨胀超声速射流压力场4阶DMD模态相位的空间分布. 与图12(d)类似, 相位空间分布清晰地表征了A2模态啸声声波传播的波振面. 图15提取了4阶DMD模态在直线R/D = 1.5上的相位分布. 沿流向, 2阶DMD模态有多个相位斜率符号发生改变的位置, 但只有一个位于饱和态区域之内, 即
$x = 2.14 D$ . 图14(c)也显示啸声声波是在饱和态区域内向外辐射, 因此, 饱和态区域外的相位斜率符号转变位置与啸声的产生没有直接关联. 此外,$x = 2.14 D$ 这一位置也与2.4节中傅里叶模态分解得到的结果一致, 进一步证实啸声A2模态的产生位置是在第3个激波栅格结构的尾缘.![]()
图 15 4阶DMD模态的相位在直线(R/D = 1.5)上的分布, Mj = 1.15Figure 15. The distribution of phase along the line (R/D = 1.5) for the fourth DMD mode, Mj = 1.153. 结 论
为了研究轴对称模态啸声的产生位置, 本文采用高阶精度有限差分方法直接求解轴对称可压缩Navier-Stokes方程, 数值模拟了完全膨胀射流马赫数分别为1.10和1.15的圆形声速喷管欠膨胀超声速冷射流, 得到了轴对称啸声A1模态和A2模态. 通过激波栅格结构位置和啸声模态频谱与文献实验结果的对比, 验证了数值模拟的准确性.
基于傅里叶模态分解、本征模态分解和动态模态分解, 分析了射流时序压力场和速度场, 研究了啸声关联拟序流动结构的空间演化, 精确定位了轴对称模态啸声的有效声源位置.
研究表明: 啸声关联拟序流动结构存在饱和态区域, 并且啸声声波是在其饱和态区域内产生并向外辐射. 在本文所涉及的射流马赫数范围内, A1和A2两种轴对称模态啸声的有效声源位置分别是在第4和第3个激波栅格结构的尾缘. 因此, 本文的数值模拟结果支持A1和A2两种啸声模态是单声源的观点.
需要注意的是, 本文只研究了低马赫数超声速射流中的轴对称模态啸声, 对于出现在高马赫数射流中的摆动模态啸声和螺旋模态啸声没有涉及, 它们具体的声源位置以及是否为单声源, 仍有待进一步研究.
-
![]()
图 1 计算域和边界条件示意图
Figure 1. The schematic of the computational domain and boundary conditions
![]()
![]()
![]()
图 3 射流近场噪声频谱信息: 喷管唇口壁面监测点[0.0D, 0.642D]
Figure 3. Spectral information of jet near-field noise: monitor [0.0D, 0.642D] located on the nozzle exit lip wall
![]()
图 4 Mj=1.10的欠膨胀超声速射流中与啸声A1模态相关的瞬时流动结构和声场
(数值纹影$\left| {\nabla \rho } \right|$: 橙红色连续云图, 取值范围[0, 2]; 胀量$\nabla \cdot {\boldsymbol{V}} $: 灰色连续云图, 取值范围[−0.02, 0.02], 白色正值表示膨胀, 黑色负值代表压缩; 涡量: 彩色等值线, 取值范围[−0.5, 20], 15级)
Figure 4. The instantaneous flow structures and acoustic field associated with screech A1 mode in the underexpanded supersonic jet with Mj=1.10
(Numerical schlieren$\left| {\nabla \rho } \right|$: orange red continuous contour from 0 to 2; dilatation$\nabla \cdot {\boldsymbol{V}}$: gray continuous contour from −0.02 to 0.02, positive value with white color means expansion, negative value with black color means compression; vorticity: multi-color banded contour from −0.5 to 20, 15 levels)
![]()
图 6 Mj=115的欠膨胀超声速射流中啸声A2模态频率上的相位空间分布
Figure 6. The spatial distribution of phase at the frequency of screech A2 mode in underexpanded supersonic jet with Mj=1.15
![]()
图 7 Mj=1.15的欠膨胀超声速射流中啸声A2模态频率上的相位沿直线R/D=1.5的分布
Figure 7. The distribution of phase at the frequency of screech A2 mode along the line of R/D=1.5 in underexpanded supersonic jet with Mj=1.15
![]()
图 8 脉动速度场的POD时间模态和特征值, Mj=1.10
Figure 8. Temporal POD modes and their eigenvalues of fluctuating velocity field, Mj = 1.10
![]()
图 9 脉动速度场的POD时间模态的傅里叶分析, Mj=1.10
Figure 9. The FFT analysis of the first and the third pairs of temporal POD modes of fluctuating velocity field, Mj = 1.10
![]()
图 10 脉动速度场的POD空间模态, Mj=1.10
Figure 10. The first and the third pairs of spatial POD modes of fluctuating velocity field, Mj = 1.10
![]()
图 11 脉动速度场轴向切片POD分析中相关模态的特征值沿流向的变化, Mj = 1.10
Figure 11. The variations of eigenvalues along streamwise of the POD modes based on the POD analysis of velocity field’s slices in different axial position, Mj = 1.10
![]()
图 12 Mj=1.10的欠膨胀超声速射流的脉动压力场DMD模态
Figure 12. The DMD modes of fluctuating pressure field in the underexpanded supersonic jet with Mj = 1.10
![]()
图 13 2阶DMD模态的相位在直线(R/D = 2.1)上的分布, Mj=1.10
Figure 13. The distribution of phase along the line R/D = 2.1 for the second DMD mode, Mj=1.10
![]()
图 14 Mj=1.15的欠膨胀超声速射流的脉动压力场DMD模态
Figure 14. The DMD modes of fluctuating pressure field in the underexpanded supersonic jet with Mj = 1.15
![]()
图 15 4阶DMD模态的相位在直线(R/D = 1.5)上的分布, Mj = 1.15
Figure 15. The distribution of phase along the line (R/D = 1.5) for the fourth DMD mode, Mj = 1.15
表 1 欠膨胀超声速射流的关键流动参数
Table 1 The key flow parameters in underexpanded supersonic jet
Mj Ma ReD 1.10 0.987 5.481 × 105 1.15 1.023 6.052 × 105 
下载: 导出CSV
表 2 啸声模态及其谐频的频率和声压级
Table 2 The frequencies and amplitudes of screech modes and their harmonics
Screech modes Numerical results Experimental results Present Ref. [41] Ref. [13−14] Mj 1.10 1.15 1.10 1.15 1.15 A1 Freq/Hz 8560 − 9040 7690 − SPL/dB 138.4 − 139 137 − A2 Freq/Hz − 9250 − 9410 8682 SPL/dB − 135.5 − 125 147 A1
harmonic 1stFreq/Hz 17070 − 18073 15337 − SPL/dB 117 − 109 122 − A2
harmonic 1stFreq/Hz − 18526 − − 17332 SPL/dB − 115 − − 116 A2
harmonic 2ndFreq/Hz − 27879 − − − SPL/dB − 111 − − −
下载: 导出CSV
-
[1] Powell A. On the mechanism of choked jet noise. Proceedings of the Physical Society of London, 1953, 66: 1039-1056 doi: 10.1088/0370-1301/66/12/306
[2] Raman G. Advances in understanding supersonic jet screech: Review and perspective. Progress in Aerospace Sciences, 1998, 34: 45-106 doi: 10.1016/S0376-0421(98)00002-5
[3] Raman G. Supersonic jet screech: Half-century from powell to the present. Journal of Sound and Vibration, 1999, 225(3): 543-571 doi: 10.1006/jsvi.1999.2181
[4] Edgington-Mitchell D. Aeroacoustic resonance and self-excitation in screeching and impinging supersonic jets—A review. International Journal of Aeroacoustics, 2019, 18(2-3): 118-188 doi: 10.1177/1475472X19834521
[5] Merle M. Sur la fréquence des ondes sonores émises par un jet d'air a grande vitesse. Comptes-Rendus de l'Académie des Sciences de Paris, 1956, 243: 490-493 (in German)
[6] Merle M. Nouvelles recherches sur les fréquences ultrasonores émises par les jets d'air. Annales des Télécommunications, 1957, 12(12): 424-426
[7] Davies MG, Oldfield DES. Tones from a choked axisymmetric jet. ii. the self-excited loop and mode of oscillation. Acustica, 1962, 12(4): 267-277
[8] Powell A, Umeda Y, Ishii R. Observations of the oscillation modes of choked circular jets. Journal of the Acoustical Society of America, 1992, 92(5): 2823-2836 doi: 10.1121/1.404398
[9] Ponton MK, Seiner JM. Acoustic study of B helical mode for choked axisymmetric nozzle. AIAA Journal, 1995, 33(3): 413-420 doi: 10.2514/3.12402
[10] Clem MM, Zaman K, Fagan AF. Variation of shock-spacing during screech stage-jumps. International Journal of Aeroacoustics, 2016, 15(3): 324-335 doi: 10.1177/1475472X16630888
[11] Tam CKW. Supersonic jet noise. Annual Review of Fluid Mechanics, 1995, 27: 17-43 doi: 10.1146/annurev.fl.27.010195.000313
[12] Shen H, Tam CKW. Three-dimensional numerical simulation of the jet screech phenomenon. AIAA Journal, 2002, 40(1): 33-41 doi: 10.2514/2.1638
[13] Ponton MK, Seiner JM. The effects of nozzle exit lip thickness on plume resonance. Journal of Sound and Vibration, 1992, 154: 531-549 doi: 10.1016/0022-460X(92)90784-U
[14] Ponton MK, Seiner JM, Brown MC. Near field pressure fluctuations in the exit plane of a choked axisymmetric nozzle. NASA Technical Memorandum 113137, 1997: 1-91
[15] Li XD, Gao JH. Numerical simulation of the generation mechanism of axisymmetric supersonic jet screech tones. Physics of Fluids, 2005, 17: 085105
[16] Gao JH, Li XD. A multi-mode screech frequency prediction formula for circular supersonic jets. Journal of the Acoustical Society of America, 2010, 127(3): 1251-1257 doi: 10.1121/1.3291001
[17] Manning TA. A Numerical investigation of sound generation in supersonic jet screech. [PhD Thesis]. Stanford: Stanford University, California, 1999
[18] Suzuki T, Lele SK. Shock leakage through an unsteady vortex-laden mixing layer: application to jet screech. Journal of Fluid Mechanics, 2003, 490: 139-167 doi: 10.1017/S0022112003005214
[19] Berland J, Bogey C, Bailly C. Numerical study of screech generation in a planar supersonic jet. Physics of Fluids, 2007, 19: 075105
[20] Edgington-Mitchell D, Weightman J, Soria J, et al. Sound production by shock leakage in supersonic jet screech//2018 AIAA/CEAS Aeroacoustics Conference, June 25-29, 2018, Atlanta, Georgia AIAA-2018-3147
[21] Edgington-Mitchell D, Weightman J, Lock S, et al. The generation of screech tones by shock leakage. Journal of Fluid Mechanics, 2021, 908: A46
[22] Umeda Y, Ishii R. On the sound sources of screech tones radiated from choked circular jets. Journal of the Acoustical Society of America, 2001, 110(4): 1845-1858 doi: 10.1121/1.1402620
[23] Edgington-Mitchell D, Oberleithner K, Honnery DR, et al. Coherent structure and sound production in the helical mode of a screeching axisymmetric jet. Journal of Fluid Mechanics. 2014, 748, 822-847
[24] Panda J. An experimental investigation of screech noise generation. Journal of Fluid Mechanics, 1999, 378: 71-96
[25] Shariff K, Manning TA. A ray tracing study of shock leakage in a model supersonic jet. Physics of Fluids, 2013, 25: 076103
[26] Mercier B, Castelain T, Bailly C. Experimental characterization of the screech feedback loop in underexpanded round jets. Journal of Fluid Mechanics, 2017, 824: 202-229 doi: 10.1017/jfm.2017.336
[27] Mercier B, Castelain T, Bailly C. Experimental study of the coherent vorticity in slightly under-expanded supersonic screeching jets. International Journal of Aeroacoustics, 2019, 18(2-3): 207-230 doi: 10.1177/1475472X19834530
[28] Li XR, He F, Zhang XW, et al. Shock motion and flow structure of an underexpanded jet in the helical mode. AIAA Journal, 2019, 57(9): 3943-3953 doi: 10.2514/1.J058024
[29] Li XR, Zhang XW, Hao PF, et al. Acoustic feedback loops for screech tones of underexpanded free round jets at different modes. Journal of Fluid Mechanics, 2020, 902: A17
[30] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. Journal of Computational Physics, 1996, 126: 202-228 doi: 10.1006/jcph.1996.0130
[31] Leer BV. Flux vector splitting for the Euler equations. ICASE Report No. 82-30, 1982
[32] Zhang SH, Shu CW. A new smoothness indicator for the WENO schemes and its effect on the convergence to steady state solutions. Journal of Scientific Computing, 2007, 31: 273-305
[33] Zhang SH, Jiang SF, Shu CW. Improvement of convergence to steady state solutions of Euler equations with the WENO schemes. Journal of Scientific Computing, 2011, 47: 216-238 doi: 10.1007/s10915-010-9435-5
[34] Thompson KW. Time dependent boundary conditions for hyperbolic systems. Journal of Computational Physics, 1987, 68: 1-24 doi: 10.1016/0021-9991(87)90041-6
[35] Thompson KW. Time dependent boundary conditions for hyperbolic systems, II. Journal of Computational Physics, 1990, 89: 439-461 doi: 10.1016/0021-9991(90)90152-Q
[36] Poinsot TJ, Lele SK. Boundary conditions for direct simulation of compressible viscous flows. Journal of Computational Physics, 1992, 101: 104-129 doi: 10.1016/0021-9991(92)90046-2
[37] Tam CKW. Advances in numerical boundary conditions for computational aeroacoustics. Journal of Computational Acoustics, 1998, 6(4): 377-402 doi: 10.1142/S0218396X98000259
[38] Li H, Luo Y, Zhang SH. Assessment of upwind/symmetric WENO schemes for direct numerical simulation of screech tone in supersonic jet. Journal of Scientific Computing, 2021, 87(3): 1-39
[39] André B, Castelain T, Bailly C. Broadband shock-associated noise in screeching and non-screeching underexpanded supersonic jets. AIAA Journal, 2013, 51(3): 665-673 doi: 10.2514/1.J052058
[40] André B, Castelain T, Bailly C. Experimental exploration of underexpanded supersonic jets. Shock Waves, 2014, 24: 21-32 doi: 10.1007/s00193-013-0457-4
[41] Loh CY, Hultgren LS, Jorgenson PCE. Near field screech noise computation for an underexpanded supersonic jet by the CE/SE method//7th AIAA/CEAS Aeroacoustics Conference, Maastricht, The Netherlands, May 28-30, 2001, AIAA-2001- 2252
[42] Lumley JL. The structure of inhomogeneous turbulence flows//Yaglom AM, Tatarsky VI, eds. Proceedings of Atmospheric Turbulence and Radio Wave Propagation. Moscow, USSR, Nauka: 166-178, 196
[43] Sirovich L. Turbulence and dynamics of coherent structures. Part I: Coherent structures. Quarterly of Applied Mathematics, 1987, 45(3): 561-571
[44] Edgington-Mitchell D, Honnery DR, Soria J. Multi-modal instability in the weakly underexpanded elliptic jet. AIAA Journal, 2015, 53(9): 2739-2749 doi: 10.2514/1.J053738
[45] Schmid PJ. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics, 2010, 656: 5-28 doi: 10.1017/S0022112010001217
[46] Kutz JN, Brunton S, Fu X. Multi-resolution analysis of dynamical systems using dynamic mode decomposition//Proceedings of World Congress on Engineering 2015, Vol I, WCE 2015, London, July 1-3, 2015
-
期刊类型引用(15)
1. Min YI，Ming XUE，Peihong CONG，Yang SONG，Haiyang ZHANG，Lingfeng WANG，Liucheng ZHOU，Yinghong LI，Wanlin GUO. Machine learning for predicting fatigue properties of additively manufactured materials. Chinese Journal of Aeronautics. 2024(04): 1-22 . 
必应学术
2. 刘翠丽，胡剑南，李建中，石俊杰，高宣雯，于凯. 醇盐体系增材制造钛合金的电化学抛光机理. 材料与冶金学报. 2024(02): 190-196 . 百度学术
3. 于飞，廉艳平，李明健，高汝鑫. 金属增材制造晶体塑性有限胞元自洽聚类分析方法. 力学学报. 2024(07): 1916-1930 . 本站查看
4. 胡雅楠，余欢，吴圣川，奥妮，阚前华，吴正凯，康国政. 基于机器学习的增材制造合金材料力学性能预测研究进展与挑战. 力学学报. 2024(07): 1892-1915 . 本站查看
5. 罗诚，袁荒. 基于张量化微结构表征的筏化镍基单晶高温合金疲劳寿命评估. 力学学报. 2024(07): 2029-2050 . 本站查看
6. Ya-qing Hou，Xiao-qun Li，Wei-dong Cai，Qing Chen，Wei-ce Gao，Du-peng He，Xue-hui Chen，Hang Su. Research progress in CALPHAD assisted metal additive manufacturing. China Foundry. 2024(04): 295-310 . 必应学术
7. 於之杰，郭玉佩，孙汉斌，张京楠，孙侠生. 先进材料及工艺的结构完整性研究进展. 航空学报. 2024(18): 33-54 . 百度学术
8. 常珂，梁晨光，易敏. 基于离散元与相场法的激光选区熔化数值模拟. 计算力学学报. 2024(05): 830-836 . 百度学术
9. 肖庆晖，张仁嘉，刘士杰，胡文轩，吕晨晞，朱思瑛，易敏. 增材铜合金拉伸力学行为的卷积神经网络预测. 计算力学学报. 2024(05): 843-850 . 百度学术
10. 易敏，胡文轩. 晶体塑性模型及其在金属疲劳寿命预测中的应用. 南京航空航天大学学报. 2023(01): 12-27 . 百度学术
11. 叶万蓉，曾国伟，闫相木，邱乙，张博文. 3D打印金属材料疲劳实验及数值模拟研究综述. 科技资讯. 2023(02): 135-138 . 百度学术
12. 詹志新，高同州，刘传奇，吴圣川. 基于数据驱动的增材制造铝合金的疲劳寿命预测. 固体力学学报. 2023(03): 381-394 . 百度学术
13. 刘海林，易敏，王建祥，易新. 激光选区熔化铺粉过程的数值模拟及粉层表征. 力学学报. 2023(09): 1921-1938 . 本站查看
14. 易敏，张璇，胡文轩，周留成，刘士杰，郑大勇. 激光冲击强化改善增材制造金属疲劳性能. 航空制造技术. 2023(20): 38-49 . 百度学术
15. 高柏森，黄玮，王生楠，张霜银，陈先民. 增材制造Ti-6Al-4V合金断裂行为与应力三轴度关系研究. 西北工业大学学报. 2022(05): 962-969 . 百度学术
其他类型引用(4)
计量
- 文章访问数: 734
- HTML全文浏览量: 313
- PDF下载量: 107
- 被引次数: 19
