0 引言

激波和旋涡是可压缩流动的基本结构。在许多流动中，比如超声速混合层、超声速射流和燃烧等，激波与旋涡共存，它们之间会发生相互作用，如激波与激波之间、激波与旋涡之间以及旋涡与旋涡之间的相互作用，这种作用，不仅会引起激波和旋涡的变形，还会产生噪声。因此，吸引了很多学者开展理论、实验和数值模拟研究。

1955年，Hollingsworth和Richards^[1]首先研究了激波与旋涡的相互作用，他们利用激波管产生平面激波，该激波扫过一个具有一定攻角的机翼后，产生一个旋涡，激波碰到激波管端部的固壁后反射回来，与旋涡发生干扰。该实验发现，激波与旋涡干扰后，伴随着激波与旋涡变形会产生声波。Dosanjh和Week^[2]对此问题进行了较细致的实验，测量了第一道声波的强度，成为后来实验研究和数值模拟研究的一个标准。Barbosa和Skews^[3]设计了一个分叉型激波管，激波管所产生的激波被分成两支，一支通过激波管内的拐角产生旋涡，另一支直接与旋涡发生干扰，研究发现，激波与旋涡干扰后，在旋涡中心附近一个很小的区域，出现压力脉冲区，其脉冲压力是周围压力的两倍以上。

为了解释声波的产生机理，早在20世纪60年代，Ribner^[4]就发展了线性分析理论。他利用傅立叶变换将旋涡分解成平面正弦剪切波，这种单正弦波与激波相互作用后产生声波，经过合成后得到的激波与旋涡干扰产生的声波在垂直于激波的方向反对称分布，其理论结果与Dosanjh和Week^[2]的测试结果的反对称部分一致。后来，Week和Dosanjh^[5]发展了Lightill^[6]的声比拟理论，将周向压力分布展开成四极、两极和单极声源之和。单极的作用使得声波偏离了单纯反对称态，他们所得到的声波与其实验一致。

到20世纪90年代，伴随着高阶精度捕捉激波格式的成熟，采用数值模拟研究激波和旋涡干扰并捕捉声波已成为研究声波产生机理的重要手段。Ellzey等^[7-9]采用四阶差分格式，通过求解非定常可压缩Euler方程，系统地模拟了激波与旋涡的干扰，得到的第一道声波强度与实验测得的声波强度基本一致。Inoue等^[10-12]利用六阶紧致格式直接模拟了不同强度激波与旋涡相互作用及声波的产生过程，给出了弱激波与旋涡相互作用声波的产生及传播机理。他们通过精细的模拟发现，当激波碰到旋涡的瞬间，旋涡边缘发生局部压缩和膨胀，形成双极的第一道声波，随着激波与旋涡干扰的继续，出现第二组压缩和膨胀区，此压缩和膨胀区与前面出现的压缩和膨胀区域交替出现，形成第一道完整的声波，第一道声波演变成四极波。激波通过旋涡后，在旋涡两侧，形成两道反射激波，反射激波后伴随第二道和第三道声波。Grasso和Pirrozli^[13-16]也系统地模拟了较大范围激波与旋涡的干扰问题，并根据激波变形的形状，将激波与旋涡的干扰分成三类：弱干扰，具有规则反射的强干扰, 具有马赫反射的强干扰。

除了激波与旋涡之间的相互作用外，在包含激波和旋涡的复杂流场中，还存在旋涡与旋涡之间的相互作用，两个或多个旋涡之间的相互作用也是噪声产生的重要机制。比如，两个同向旋转的旋涡的合并过程，是剪切层噪声的重要部分。

这种包含激波、旋涡的复杂流场的非定常运动产生的噪声是激波噪声。目前，虽然气动声学理论已有相当的发展，但是关于激波噪声的产生机制认识还很不够。伴随高阶精度数值格式和大型并行计算机的快速发展，数值模拟已成为研究激波噪声产生机制的重要手段，但是激波噪声的计算仍然是计算气动声学所面临的相当严峻的挑战，具体表现在：

(1) 流场的强非定常特性。激波噪声是激波与其它流场结构相互作用产生的，其中存在激波变形和产生过程，与激波相互作用的流场结构也随之发生变化。比如，在激波与旋涡的相互作用过程中^[17-20]，入射激波会变形并产生反射激波，在旋涡中心附近还产生小激波串，同时旋涡被压扁，而且出现扭转运动，当激波强度达到一定程度之后，会引起旋涡的破碎。在激波与剪切层相互作用过程中，存在有激波的局部震荡^[21]。因此，传统的定常流体动力学问题计算方法或是大尺度问题的非定常方法已不再适用。尽管已发展了多种高阶精度Runge-Kutta方法^[22]，但大多Runge-Kutta方法是显式格式，由于稳定性条件的限制，允许的时间步长非常小，对包含边界层的流动计算效率非常低，因此仍需要发展时间离散精度和效率更高的计算方法。

(2) 流场的多尺度特性。激波噪声是流场的一部分，但是在很多情况下，它的扰动尺度比流体动力学尺度小很多^[23]。比如在马赫数为1.5的喷流中，距喷流直径40倍处测得的噪声强度是124 dB，声场的脉动速度与喷流的速度之比是1.5×10^-4。特别地，当流场中存在强激波时，声波引起的物理量脉动比激波引起的变化尺度差距更大。这就要求不仅能精确计算宏观尺度的流体动力学量，能够捕捉激波，还能精确计算微弱尺度的声波，对计算方法和计算网格的要求远比动力学计算高。

(3) 模拟激波噪声所需的计算区域远大于模拟气动力所需的计算区域。研究激波噪声既关心近场流体动力学的变化，挖掘噪声产生的机制，又关心远场激波噪声的传播特性，计算区域边界与声源的距离非常远，这是与传统计算流体动力学的鲜明差异。

(4) 激波与声波的计算对格式需求存在矛盾的因素。捕捉激波需要格式具有一定的耗散，而声波是近乎无色散、无耗散的等熵运动，声波的模拟需要高分辨率和无耗散格式。因此，现在还不存在一种较为理想的高阶精度、高分辨率、低耗散的激波捕捉格式。目前，声波的计算普遍采用线性紧致格式，精度和分辨率非常高，耗散性小，适合于声波的计算。但是线性格式不能捕捉流场中的强激波，而捕捉激波普遍采用WENO格式^[24]。尽管WENO格式的鲁棒性非常好，可以模拟很强的激波，但同所有的捕捉激波格式一样，在激波下游会产生微弱的非物理波动。虽然这种波动不影响流场的整体特性，但其强度常常会比流场中的噪声高。另外，虽然可以设计很高精度的WENO格式^[25]，但格式的耗散性比较大，对声波的分辨率不满意。尽管已有很多种简单的混合格式^[26]、优化WENO格式^[27]、非线性加权紧致格式^[28-29]以及DRP格式^[30]，试图提高捕捉激波类格式的分辨率，但实际应用比WENO格式并没有明显的优势，因此目前还不存在一种既能捕捉激波又适用于声波计算的较理想数值方法。

针对这些问题，我们近年开展了较为系统性的研究，解决了计算激波噪声的一些方法上存在的问题，对包含有激波、旋涡和声波的一系列问题进行了直接数值模拟，揭示了激波噪声产生机制。本文对这些研究做简要回顾。

1 针对激波噪声计算的数值方法的研究 1.1 WENO格式的不收敛问题^[31-32]

WENO格式是一种高阶精度捕捉激波的格式^[24], 被广泛用于工程实际和科学计算。但是像所有其它的激波捕捉格式一样，用WENO格式计算定常激波时，在激波下游会出现微弱的非物理波动，导致残差不能收敛到满意程度。虽然这种波动不会影响整体的气动特性，但是至少会造成两方面的影响，第一，残差是结束计算的重要判据，残差无法收敛到很小的值会使计算无法判断何时停止，只能靠经验结束进程；第二，虽然激波后的非物理波动是微弱的，不影响气动力等宏观流动特性，但是对于激波噪声，这种微弱的非物理波动则非常大，有时比激波噪声的幅度高几个数量级，因此计算气动噪声问题，应该消除这种微弱的非物理波动。经过系统的数学物理分析和大量的数值实验，我们提出了两种方法消除这种微弱的非物理波动^[31-32]。

1) 针对五阶WENO格式，我们提出了一种新的光滑测试因子如下：

$ \begin{array}{l} I{S_r} = 4\int_{{x_{i - \frac{1}{2}}}}^{{x_{i + \frac{1}{2}}}} {\Delta x{{\left( {\frac{{\partial {p_r}\left( x \right)}}{{\partial x}}} \right)}^2}{\rm{d}}x} - \\ \;\;\;\;\;\;\;\;\;\;\frac{1}{3}\int_{{x_{i - \frac{1}{2}}}}^{{x_{i + \frac{1}{2}}}} {\Delta {x^3}{{\left( {\frac{{{\partial ^2}{p_r}\left( x \right)}}{{{\partial ^2}x}}} \right)}^2}{\rm{d}}x} \end{array} $

(1)

具体形式为：

$ \left\{ \begin{array}{l} I{S_0} = {\left( {{f_{i - 2}} - 4{f_{i - 1}} + 3{f_i}} \right)^2}\\ I{S_1} = {\left( {{f_{i - 1}} - {f_{i + 1}}} \right)^2}\\ I{S_0} = {\left( {3{f_i} - 4{f_{i + 1}} + {f_{i + 2}}} \right)^2} \end{array} \right. $

(2)

采用这种光滑测试因子后，消除了激波下游的微弱非物理波动，对典型问题计算可以收敛到机器零。

2) 前面的方法只适合于五阶精度的格式，对更高阶精度格式不适用，而且会引起格式在高阶极值点附近降阶。经过大量的数值实验，我们发现WENO格式在特征投影过程中，采用Roe平均方法计算半点上的物理量确定特征值和特征矩阵，是影响WENO格式不收敛的重要原因。为此，我们提出采用迎风插值方法代替原来的Roe平均方法。即：

当${u_{i + \frac{1}{2}}} \ge 0$时，取${u_{i + \frac{1}{2}}} = {U^{\left( 1 \right)}}$；

当${u_{i + \frac{1}{2}}} < 0$时，取${u_{i + \frac{1}{2}}} = {U^{\left( 2 \right)}}$。

具体的迎风插值公式如下：

$ {\rm{一阶插值:}}\;\;\;\;\;\;{U^{\left( 1 \right)}} = {U_i}\;\;\;{U^{\left( 2 \right)}} = {U_{i + 1}} $

(3)

$ {\rm{二阶插值:}}\;\;\;\;\;\;\left\{ {\begin{array}{*{20}{l}} {{U^{\left( 1 \right)}} = \left( {3{U_i} - {U_{i - 1}}} \right)/2}\\ {{U^{\left( 2 \right)}} = \left( {3{U_{i + 1}} - {U_{i + 2}}} \right)/2} \end{array}} \right. $

(4)

最优加权插值：

$ U_{i + \frac{1}{2}}^{{\rm{WENOI}}} = \sum\limits_{k = 0}^{r - 1} {{\omega _k}U_{i + \frac{1}{2}}^k} $

(5)

$ U_{i + \frac{1}{2}}^k = \sum\limits_{l = 0}^{r - 1} {{c_{kl}}{U_{i - k + l}}} $

(6)

最优插值精度与WENO重构过程精度一致，如五阶精度的WENO格式采用五阶插值，七阶精度的WENO格式采用七阶插值，九阶精度的WENO格式采用九阶插值。图 1和图 2是采用两种方法对马赫数2的一维定常激波计算得到的密度分布和收敛历程，图 1中Res_A为残差的平均值。可以看到，改进的WENO格式可以收敛到机器零。值得指出的是，采用最优加权插值方法求解半点上的值，使WENO格式具有自相似特性，并一致高阶精度，详见文献[31-32]。

1.2 一种类谱分辨率的高阶精度紧致格式^[34-35]

WENO格式是一种高阶精度捕捉激波格式，虽然格式的精度可以很高，但是对声波的分辨率并不理想，而且耗散性也比较大，采用WENO格式模拟湍流结构、激波噪声等复杂多尺度流动并不理想。为了提高格式的分辨率，降低耗散，并保持捕捉激波的鲁棒性，我们将WENO格式的设计思想应用到紧致格式中，设计了一种类谱分辨率的高阶精度紧致格式。这一格式是基于Lele的半点型线性紧致格式^[33]设计的, Lele的半点型线性紧致格式形式是：

$ \begin{array}{l} \beta {{f'}_{j - 2}} + \alpha {{f'}_{j - 1}} + {{f'}_j} + \alpha {{f'}_{j + 1}} + \beta {{f'}_{j + 2}} = \\ \;\;\;\;\;\;\;a\frac{{{f_{j + \frac{1}{2}}} - {f_{j - \frac{1}{2}}}}}{{\Delta x}} + b\frac{{{f_{j + \frac{3}{2}}} - {f_{j - \frac{3}{2}}}}}{{3\Delta x}} + c\frac{{{f_{j + \frac{5}{2}}} - {f_{j - \frac{5}{2}}}}}{{5\Delta x}} \end{array} $

(7)

该格式左端是网格点上的函数导数值，右端是半点的函数值，它最大的优点是分辨率非常高，接近谱方法分辨率。但是右端半点的函数值是未知的，Lele^[33]设计了一种线性紧致插值的方法求得半点的值。图 3是Lele的半点型线性紧致格式的设计模板，可以看出，该格式有两个缺点：(1)设计模板包括了整点和半点的函数值，但Lele仅用了部分信息，格式没有达到该模板上的最优精度；(2)半点上的值需要插值方法求得，任何插值方法都存在传递误差，导致格式对高波数的分辨率大幅度下降。

为克服Lele格式的缺点，我们设计了一种类谱分辨率的高阶精度紧致格式。在网格节点上，我们采用如下格式计算其一阶导数：

$ \begin{array}{l} \beta {{f'}_{j - 2}} + \alpha {{f'}_{j - 1}} + {{f'}_j} + \alpha {{f'}_{j + 1}} + \beta {{f'}_{j + 2}} = \\ \;\;\;\;\;\;\;a\frac{{{f_{j + \frac{1}{2}}} - {f_{j - \frac{1}{2}}}}}{{\Delta x}} + b\frac{{{f_{j + 1}} - {f_{j - 1}}}}{{2\Delta x}} + \\ \;\;\;\;\;\;\;c\frac{{{f_{j + \frac{3}{2}}} - {f_{j - \frac{3}{2}}}}}{{3\Delta x}} + d\frac{{{f_{j + 2}} - {f_{j - 2}}}}{{4\Delta x}} + e\frac{{{f_{j + \frac{5}{2}}} - {f_{j - \frac{5}{2}}}}}{{5\Delta x}} \end{array} $

(8)

采用同样的格式计算半点上的一阶导数值：

$ \begin{array}{l} \beta {{f'}_{j - \frac{5}{2}}} + \alpha {{f'}_{j - \frac{3}{2}}} + {{f'}_{j - \frac{1}{2}}} + \alpha {{f'}_{j + \frac{1}{2}}} + \beta {{f'}_{j + \frac{3}{2}}} = \\ \;\;\;\;\;\;\;a\frac{{{f_j} - {f_{j - 1}}}}{{\Delta x}} + b\frac{{{f_{j + \frac{1}{2}}} - {f_{j - \frac{3}{2}}}}}{{2\Delta x}} + \\ \;\;\;\;\;\;\;c\frac{{{f_{j + 1}} - {f_{j - 2}}}}{{3\Delta x}} + d\frac{{{f_{j + \frac{3}{2}}} - {f_{j - \frac{5}{2}}}}}{{4\Delta x}} + e\frac{{{f_{j + 2}} - {f_{j - 3}}}}{{5\Delta x}} \end{array} $

(9)

求得网格点和半点处的一阶导数之后，通过求解相应的控制方程，即可求得下一时刻的网格点和半点处函数值。

与Lele的格式相比，式(8)、式(9)给出的数值格式有两个优点：1)格式的精度有一定的提高，Lele格式最高达到十阶精度，我们设计的格式可以达到十四阶；2)网格点和半点上的函数值都是通过同一格式计算得到，消除了原格式紧致插值引起的传递误差，保证了计算过程的精度和分辨率。图 4给出了六阶和八阶格式修正波数及其与Lele格式的比较，相比之下我们改进的格式分辨率与精确值相差无几，详见文献[34]。

线性格式不能捕捉流场中的强激波，为了使格式具有捕捉强激波的功能，我们将WENO重构思想应用到插值技术，发展了高阶WENO插值方法，并将该方法应用于式(8)、式(9)右端，初步发展了一种非线性类谱分辨率的紧致格式。具体的WENO插值方法是：在格式的设计模板S^2r－1={x_i－r+1, ..., x_i+r－1}上，采用插值函数:

$ {{\hat f}^{2r - 1}}\left( x \right) = {f_i} + \sum\limits_{l = 1}^{2\left( {r - 1} \right)} {{a_l}{{\left( {x - {x_i}} \right)}^l}} $

(10)

可以得到半点上的通量值：

$ \hat f_{i + 1/2}^L = {q^{2r - 1}}\left( {{f_{i - r + 1}}, \ldots ,{f_{i + r - 1}}} \right) $

(11)

与WENO格式重构类似，将模板S^2r－1分成三个子模板，在每个子模板上，存在r阶插值：

$ \hat f_{i + 1/2}^{\left( k \right)} = q_k^r\left( {{f_{i + k - r + 1}}, \ldots ,{f_{i + k}}} \right) $

(12)

采用加权思路，可以得到半点上的加权插值：

$ \hat f_{i + 1/2}^{{\rm{WENO}}} = \sum\limits_{k = 0}^{r - 1} {\omega _k^rq_k^r} $

(13)

在式(8)的右端，采用如下混合插值方法替代原来的精确值f_i+1/2，即：

$ \hat f_{i + 1/2}^{{\rm{hybrid}}} = \left( {1 - \sigma } \right)\hat f_{i + 1/2}^{{\rm{WENO}}} + \sigma {f_{i + 1/2}} $

(14)

而在式(9)中，采用类似的混合插值方法。得到的非线性格式，能够保证格式的精度，具有较好的捕捉激波能力，同时格式的分辨率比现有的WENO明显提高，耗散性也大幅度降低，详见文献[35]。

2 包含激波、旋涡和声波的复杂多尺度流动数值模拟 2.1 两个同向旋转的Gaussian涡合并过程^[20]

Gaussian涡是一种典型的旋涡，在一定的条件下，两个同向旋转的Gaussian涡会合并^[36]，对于合并条件，已有很多研究^[37]。1995年，Mitchell等^[36]采用直接数值模拟方法，研究了两个同向旋转的Gaussian旋涡的合并过程，给出了合并过程中所产生的声。自此，这一问题就成了计算气动声学一个标准问题。采用前述我们发展的类谱分辨率的线性紧致格式，我们也直接数值模拟了两个同向旋转的Gaussian旋涡的合并过程^[20]。Gaussian涡的初始涡量场为$\omega = \frac{\mathit{\Gamma } }{{2\pi r}}\left( {1 - {e^{ - \alpha {r^2}/R_c^2}}} \right)$，其初始结构见图 5。图 6是两个同向旋转的Gaussian涡的合并历程。耦合作用引起两个旋涡相互长时间绕旋，在这一过程中，旋涡间的距离基本保持不变。在合并前，两个旋涡的耦合作用将旋涡压缩成椭圆形，并且使得它们互相靠近。随后突然间合并为单个旋涡。一开始旋涡具有两个长臂结构，最终逐渐演化成圆形旋涡。在两个同向旋涡合并的过程中会产生强噪声。Mitchell等^[36]的数值结果显示，声波具有正弦属性，并且振幅在合并的瞬间达到峰值。图 7给出了远场r=λ/2和r=2λ处声波及其与Mohring声比拟理论的对比^[38]，结果非常一致。

值得指出的是，我们的结果与Mitchell的结果并不一致，其原因可能是由于初始参数的微弱不一致性，导致两个旋涡合并时间的差异。为了验证我们结果的正确性，我们进行了网格收敛性研究，初始条件的影响以及采用不同方法对此问题进行了模拟，结果一致，详见文献[20]。

2.2 两个Taylor旋涡的相互作用^[20]

Taylor涡也是一种典型的旋涡，与Gaussian涡结构不同，Taylor涡具有双层结构。为了深入了解涡声的产生机理，我们系统研究了两个Taylor旋涡的相互作用，并根据旋涡特性，将两个Taylor旋涡的相互作用分成四类：(1)两个反向旋转的Taylor涡的相互作用；(2)两个同向旋转Taylor涡的相互作用；(3)两个同向旋转Taylor涡的合并；(4)两个尺度或强度相差较大的旋涡的相互作用。

2.2.1 两个反向旋转的旋涡的相互作用

Taylor涡具有双层结构，内层外层涡量符号相反。如果两个反向旋转的旋涡距离足够近，它们之间会发生耦合，结果导致内层外层分离，两个涡核逐渐靠拢，形成一个涡偶极子，两个Taylor涡的外层逐渐靠拢，形成一个反向运动的涡偶极子，如图 8所示, 是强度为M_v=0.5、初始距离为4倍旋涡半径的两个反向涡相互作用的涡量演化历程。详见文献[20]。

2.2.2 两个同向旋转Taylor旋涡的相互作用

图 9是两个强度为M_v=0.5、初始旋涡间距是4倍旋涡半径的两个同向旋转的Taylor旋涡相互作用过程的涡量演化历程。这种相互作用存在两种特性，1)两个同向旋转的旋涡将演化成两个不对称的旋涡区，每一个涡区存在三个涡核结构；2)与两个反向旋转的旋涡相互作用具有本质区别的是，每一个涡区的两个涡核来自同一个初始旋涡，强的涡核来自初始旋涡的内层，而初始旋涡的外层逐渐从初始旋涡中分离形成围绕涡核的两个独立的涡核。

2.2.3 两个同向旋转Taylor旋涡的合并过程

与两个同向旋转的Gaussian旋涡相似，在一定条件下两个同向旋转的Taylor旋涡也会合并。如图 10所示，是两个强度为M_v=0.5的同向旋转的Taylor旋涡的合并过程，两个旋涡的初始距离是旋涡半径的2倍。两个涡核合并形成主涡核，而在涡核的外部，形成两个旋转臂结构，两个涡的外层逐渐集中于旋转臂附近，形成一个三核涡结构。

2.2.4 两个强度或尺度相差较大的旋涡的相互作用

两个强度相差很大的旋涡的耦合过程是最令人感兴趣的。1)两个强度或尺度相差较大的旋涡的相互作用会产生多核涡结构；2)强度或尺度的较大差异，导致它们之间相互作用所产生的多核涡结构，具有非对称特性。这种相互作用包括有两种典型状态，一是强度相差较大的两个Taylor涡的相互作用，另一是尺度相差较大的Taylor涡的相互作用。

图 11是两个强度相差较大的Taylor涡的相互作用过程，其中强旋涡的强度为M_vu=－0.8，顺时针旋转，弱旋涡的强度为M_vd=0.25，逆时针旋转，它们之间的初始距离为半径的4倍。可以看到，与两个等强度的反向旋涡相互作用类似，两个强度相差较大的Taylor涡相互作用会产生两个涡偶极子。两个旋涡的涡核逐渐靠拢，形成一个强的涡偶极子, 两个旋涡的外层逐渐从原旋涡中分离，形成一个弱的涡偶极子。由于其强度相差较大，两个涡偶极子都呈现非对称特性，导致弱的涡核围绕强涡核旋转，较弱的涡偶极子围绕强涡偶极子旋转。

图 12是两个尺度相差较大的反向旋涡相互作用的过程，其中两个旋涡的强度都是0.5，一个逆时针旋转，一个顺时针旋转。两个旋涡的尺寸相差较大，上面的旋涡半径为1，而下面的旋涡半径为0.2，两个旋涡的初始距离为2.2。可以看到，小旋涡的涡核穿过大旋涡的外层，进入到大旋涡的涡核外围，并围绕大旋涡的涡核不断旋转，最终形成一个椭圆形的涡核，而旋涡的外层形成两个涡核结构，围绕大的涡核旋转，形成一个三核涡结构，并持续做蛙跳运动。

2.2.5 两个Taylor涡相互作用过程中声波产生机理分析

两个Taylor涡的相互作用包括丰富的动力学特性，这种过程会产生噪声。图 13和图 14是两个典型状态两个Taylor涡相互作用的声压Δp=(p－p₀)/p₀等值线和其径向分布。旋涡的强度均是0.5，两个旋涡初始距离都是半径的4倍。其中图 13是两个反向Taylor涡相互作用所产生的噪声，图 14是两个同向Taylor涡相互作用所产生的噪声。可以看到，两种旋涡相互作用产生的噪声具有很大区别，两个反向旋涡相互作用仅产生几个声脉冲，而两个同向旋转的Taylor涡相互作用，与两个同向旋转的Gaussian涡相互作用很类似，会产生持续的强噪声。

以两个尺度相差较大的旋涡相互作用为例，分析噪声产生机理。我们计算了Lamb量$\nabla \cdot \left( {\rho \omega \times u} \right)$，它代表着噪声产生的源。如图 15所示，是Lamb量的等值线演化过程。可以看到，旋涡相互作用具有三个明显的运动:1)初始圆形旋涡被扭曲，弱旋涡被撕扯成两个涡带，他们与强旋涡继续相互作用；2)尽管强旋涡的涡心位置几乎不变，但整个涡核被压成椭圆形；3)两个涡带与强旋涡的涡核之间相互作用产生一种蛙跳式运动。图 16是监测点(x, y)=(100, 0)和(0, 100)的声压随时间变化曲线，压力峰值P1、P2、P3、P4、P5、P6和P7位于时间轴103、121、139、163、190、223和262。考虑到声从涡心附近产生声波到达监测点所需要的时间，它们对应于长轴在y轴的方向，而声波谷M1、M2、M3、M4、M5、M6和M7位于时间轴114、130、151、176、206、241和285，它们对应于长轴位于x的方向。这意味着，噪声主要由旋涡的蛙跳运动产生。

2.3 激波与旋涡的相互作用^[17-19]

激波与旋涡相互作用是激波噪声的一个简化模型，通过研究激波与旋涡相互作用，可以深入了解激波噪声产生机制。通过对激波与旋涡相互作用系统的模拟，我们发现激波与强旋涡相互作用存在多级干扰。第一级干扰是入射激波和初始旋涡的相互作用，第二级干扰是反射激波和变形旋涡的相互作用，第三级干扰是涡心附近产生的小激波串和变形旋涡的干扰, 如图 17所示。与此同时，旋涡的变形也存在多级特征。在第一级干扰中，初始圆形旋涡被压成椭圆形；在第二级干扰中，椭圆形旋涡被压回成圆形旋涡；在第三级干扰中，圆形旋涡又被压缩成椭圆形。此外，这种多级干扰会引起旋涡在水平线附近做振荡运动。根据激波与强旋涡相互作用过程中的多级干扰过程，当时我们预测，这种多级干扰过程会产生更为复杂的声波。这一预测被Chatterjee和Vijayaraj^[39]的研究所证实，如图 18所示。此外，我们还系统研究了激波与涡列相互作用过程流场结构和声波产生机制^[17-19]。

2.4 激波与剪切层相互作用^[40]

激波与剪切层相互作用是一个典型的问题，也可以看做是研究超声速喷流激波噪声的一个简化模型。为了了解超声速喷流激波噪声的产生机制，我们设计了一个激波与剪切层相互作用的模型。如图 19所示，一道斜激波打到一个剪切层上，下边界是一个反射固壁，透射激波打到固壁上后，形成反射激波，反射激波再次与剪切层相互作用发生反射，形成类似于超声速喷流的激波串结构。通过较为系统的数值模拟，我们揭示了两种噪声产生机制，一是激波与旋涡相互作用，发生在图 20所示的区域1，另一是激波泄露机制，发生在图 20所示的区域2。

图 21是第一个区域入射激波与剪切层中旋涡相互作用过程的一个演化周期中几个典型时刻的胀量图，其中红色虚线标记的是剪切层下层流体中被追踪的两道声波。我们看到由于流体以超声速向下游传播，声波也随着流体向下游方向对流，并且声波的辐射半径在逐渐增大。当t=91.69时，产生了第二道声波，可以明确地看到声波的产生源在涡与小激波相互作用的位置，这和张树海等^[17-19]的激波与旋涡相互作用的研究结果相同。剪切层的涡列穿过激波，每个涡与激波发生干扰时，根据这种激波噪声产生和传播的机制，在噪声源处向外辐射出一道道的声波，同时我们可以看到当声波传播到下壁面时声波反射回内场，从而声波再次与激波和剪切层相互干扰。

图 22是第二个区域反射激波与剪切层相互作用过程中一个演化周期中几个典型时刻的胀量图，其中红色虚线标记的是剪切层下层流体中的一道声波。我们看到由于流体以超声速向下游传播，声波也随着流体向下游方向对流，并且声波的辐射半径在逐渐增大。当t=91.69时，噪声源处产生了第二道声波，可以明确地看到声波的产生源在旋涡之间的鞍点处，即辫子区的位置，这和Manning等^[41-43]的激波泄漏机制完全相同。在剪切层的涡列穿过激波时，涡对与激波发生干扰，在涡对之间的鞍点处激波透过剪切层，同时在鞍点处激波被辫子区的旋臂所阻碍，激波与辫子区的相互干扰使得激波强度变弱，部分激波穿过剪切层，还有一部分激波没有穿过剪切层，而是在鞍点处以声波的形式发生泄漏。这说明了强激波与剪切层相互作用中也存在激波泄漏机制，在鞍点处激波泄漏并向外辐射声波。根据这种激波噪声产生和传播的机制，当剪切层穿过激波时，在噪声源处向外辐射出一道道的声波，同时我们可以看到当声波传播到下壁面时声波反射回内场，从而声波再次与激波和剪切层相互干扰。

2.5 轴对称超声速喷流^[44]

超声速喷流噪声是一种典型的激波噪声，与很多工程密切相关。采用五阶WENO格式直接求解轴对称Navier-Stokes方程，数值模拟了射流马赫数(完全膨胀)为M_j=1.19(啸声呈现轴对称模态，三维效应不明显)的轴对称欠膨胀射流，其对应的射流声速马赫数为Ma=1.05，雷诺数为Re=6.216×10⁵，喷管直径为D=25.4 mm，喷管厚度为0.4D，环境温度为288.15 K。

图 23给出了喷管唇口壁面[0.0, 0.642D]和喷管出口平面[0.0, 2.0D]上压力监测点处噪声信号的频谱信息，包括频率和声压级。结果显示，在大于5000 Hz的频率区间内，有4个尖锐的声压级突起：前两个声压级突起的频率为6567 Hz(128 dB)和8637 Hz(123 dB)，分别是啸声A1模态和啸声A2模态；另外两个声压级突起的频率为12087 Hz(125 dB)和14310 Hz(123 dB)，分别是啸声B模态的谐频和啸声A0模态。

综合流场结构和声场信息可以得到的啸声模态的图像及其产生位置，其中流场结构用数值纹影(密度梯度的模)表示，声场由胀量场表示。根据啸声频率可以得到啸声的波长，利用波长可以在胀量场中清晰地分辨出啸声A1模态和A2模态，相应的波长分别为λ=2.04D和λ=1.51D。从图 24中可以清晰地分辨出向上游传播的啸声，并且啸声产生位置都是在第三个激波栅格和第五个激波栅格之间的区域。

2.6 声波穿过激波的畸变过程

声波产生在近流场区域，流体介质的非定常运动产生的噪声经过复杂流场结构传播达到远场。当噪声穿过复杂流场结构过程中，会与流场相互作用，发生畸变。其中，声波穿过激波是一个典型问题，我们采用直接数值模拟方法，研究了声波穿过激波的畸变特性。

初始条件为：

$ \rho = 1.0 + \varepsilon \sin \left[ {\omega \left( {\frac{x}{{\frac{1}{M} + 1}} - t} \right)} \right] $

(15)

$ u = 1.0 + \frac{\varepsilon }{M}\sin \left[ {\omega \left( {\frac{x}{{\frac{1}{M} + 1}} - t} \right)} \right] $

(16)

$ p = \frac{1}{{\gamma {M^2}}} + \frac{\varepsilon }{{{M^2}}}\sin \left[ {\omega \left( {\frac{x}{{\frac{1}{M} + 1}} - t} \right)} \right] $

(17)

计算区域为x∈[－10, 10]，马赫数M=2.0，声波幅值ε=1.0d－5，频率ω=0.6π，比热比取γ=1.4，计算结果如图 25所示，由图可以看出，声波穿过马赫数M=2.0的激波，振幅放大了约3.64倍多。

2.7 声波与旋涡相互作用

当声波穿过旋涡时，声波的幅值和频率都会发生变化。采用我们发展的线性紧致格式，通过求解Navier-Stokes方程，我们研究了平面声波穿过等熵涡的散射特性。

图 26是计算模型，平面声波与一个二维Taylor涡相互作用。计算区域的左端入射声波为：

$ \left( {\begin{array}{*{20}{c}} {\rho '}\\ {u'}\\ {v'}\\ {p'} \end{array}} \right) = \varepsilon \left[ {\begin{array}{*{20}{c}} 1\\ 1\\ 0\\ 1 \end{array}} \right]\sin \left( {\omega \left( {x - t} \right)} \right) $

(18)

二维Taylor涡置于计算区域的中心。图 27是入射声波频率ω=0.6π与强度为M_v=0.25的Taylor涡相互作用的瞬态脉动压力场云图。可以看到，旋涡对声波的影响非常明显。在旋涡下游，形成两条强噪声干扰条带，在强条带外侧，有二次条带，由于旋涡以逆时针方向旋转，干扰条带关于y=0不对称。图 28是散射声压的均方根(RMS)等值线图，图 29是离涡心r=10的圆周上散射声压的指向性分布。

3 结束语

近年来，针对激波噪声计算的需求，我们开展了数值方法研究，通过对典型问题的直接数值模拟，揭示了激波噪声产生机理。

(1) 在数值方法方面，针对高阶精度捕捉激波的WENO格式不收敛问题，提出了一种新的光滑测试因子，发展了一种迎风插值技术计算半点上物理量，用以确定WENO重构过程中特征投影所需物理量。所改进的WENO格式，消除了激波下游的微弱非物理波动，对典型问题计算可以收敛到机器零。发展了一种类谱方法的紧致格式，这一格式包括线性和非线性两部分，其中线性格式具有高阶精度、高分辨率和低耗散特性，是模拟低速气动声学问题的理想算法。非线性格式捕捉激波能力与WENO格式相当，分辨率明显提高，耗散性大幅度降低。

(2) 采用直接数值模拟方法，系统研究了两个旋涡相互作用、激波与旋涡/涡列相互作用、激波与剪切层相互作用以及超声速喷流激波噪声。发现激波与强旋涡相互作用具有多级干扰特性，提出每级干扰都有自己的声波产生。

当然，对激波噪声的认识还相当有限，有许多问题有待进一步解决。