0 引言

流体流动产生的噪声，如射流、混合层或边界层中的噪声，一直是航空航天及许多工程技术中的严重问题，气动噪声不仅是声污染的主要来源，而且还会造成设备的声疲劳。因此，如何理解、预测进而控制气动噪声，引起了学术界和工业界的日益重视。

早在1974年，Liu^[1]就指出喷气噪声由大尺度相干结构和小尺度湍流结构产生。Tam和Burton^[2] (1984)用渐进匹配法研究了流动近场的扰动演化和远场的声辐射，发现大尺度相干结构实质上是一种不稳定波。在亚声速混合层中，近场大尺度涡结构及其合并是典型的非线性过程；在超声速混合层中，声波是以马赫波的形式辐射向远场^[3]。

混合层气动声场的计算包含了流体动力学近场和声辐射远场这两个尺度相差极大的物理场，很难通过直接数值模拟方法计算得到。行之有效的方法就是分步求解，流动近场用非定常算法，如直接数值模拟(DNS)、大涡模拟(LES)、抛物化稳定性方程(PSE)等，声辐射远场则采用的积分外推方法^[4-5]。常见的积分外推方法有Kirchhoff积分方法^[6-7]和FW-H积分方法^[8]。

PSE是近年来发展迅速的研究边界层稳定性的方法^[9]，由于采用空间推进求解，相比直接数值模拟，可以大量节约计算存储量和计算时间。Piot^[10]采用基于线性欧拉方程的PSE方法和大涡模拟对亚声速及超声速射流进行研究，结果表明PSE可以有效预测超声速射流声辐射的方向和强度。Colonius ^[11]等用线性PSE研究射流声辐射问题，其结果与实验结果一致。Cheung^[12]使用PSE计算超声速混合层近场动力学性能与远场声辐射特性，其结果与DNS吻合得很好。Sinha等^[13]结合PSE与Kirchhoff积分方法计算射流中的气动噪声问题。

上述声辐射积分外推法所依据的数学模型，由于缺乏流场的动力学信息，多依赖经验假定，对声源的处理具有一定的随意性。例如，工程中常用的Kirchhoff积分方法，其控制面的选取对远场计算结果有很大影响。Wu^[14]采用渐进匹配法与三层结构理论研究超声速不稳定模态的非线性演化，及其马赫波辐射过程，推导出远场的声压可通过积分压力幅值函数获得，不再依赖于经验假设，本文首次将这种方法应用到混合层远场压力预测上。

本文以超声速平面混合层为研究对象，近场扰动的演化采用PSE方程计算，入口加入的扰动由稳定性理论得到。而声学远场，则采用Wu^[14]文中提出的积分方法，通过对扰动波演化及马赫波辐射的分析^[15-16]，获得远场积分所需的参数，预测远场声辐射。同时还考虑了来流温度对马赫波辐射的影响。

1 超声速混合层近场计算

超声速混合层近场流场计算分两步：首先基本流通过求解定常可压缩二维边界层方程得到^[17]，然后在入口加入不稳定波，用抛物化稳定性方程计算不稳定波向下游的演化。

1.1 基本流计算

以超声速平面混合层为基本流，控制方程如式(1)，表 1给出基本流的计算参数。

$ \left\{ \begin{array}{l} \frac{{\partial (\rho u)}}{{\partial x}} + \frac{{\partial (\rho v)}}{{\partial y}} = 0 \hfill \\ \rho u\frac{{\partial u}}{{\partial x}} + \rho v\frac{{\partial v}}{{\partial y}} = 0 \hfill \\ \rho u\frac{{\partial T}}{{\partial x}} + \rho v\frac{{\partial T}}{{\partial y}} = \frac{1}{{ Re Pr }}\frac{\partial }{{\partial y}}\left( {\mu \frac{{\partial T}}{{\partial y}}} \right) \hfill \\ \;\;\;\;\;\;\;\; + \frac{{Ma_1^2(\gamma - 1)}}{{ Re }}\mu {\left( {\frac{{\partial u}}{{\partial y}}} \right)^2} \hfill \\ \rho T = 1 \hfill \\ \end{array} \right. $

(1)

式(1)和表 1中：ρ为密度，u为流向速度，v为法向速度，T为温度，Ma₁与Ma₂为混合层上下层马赫数，μ为黏性系数，p_r为普朗特数，Re为雷诺数，x、y为空间坐标。

图 1给出流向不同站位处的基本流速度分布，图 2是流向速度分布沿法向的二阶导数。可以看出，沿着流向发展，剪切层厚度不断增加，流向速度沿法向的梯度逐渐减小。

1.2 PSE计算不稳定波演化 1.2.1 控制方程

将流场的瞬时量分解为基本流和扰动量的叠加：

$ q = \bar q + q' $

(2)

将式(2)带入无量纲N-S方程并减去基本流满足的无量纲N-S方程，获得扰动控制方程：

$ \begin{array}{l} \mathit{\Gamma }\frac{{\partial q'}}{{\partial t}} + A\frac{{\partial q'}}{{\partial x}} + B\frac{{\partial q'}}{{\partial y}} + Dq' + {E_{xx}}\frac{{{\partial ^2}q'}}{{\partial {x^2}}}\\ \;\;\;\; + {E_{xy}}\frac{{{\partial ^2}q'}}{{\partial x\partial y}} + {E_{yy}}\frac{{{\partial ^2}q'}}{{\partial {y^2}}} = {F_{mn}} \end{array} $

(3)

其中，等式左边为线性项，其系数矩阵Γ、A、B、D、E_xx、E_xy、E_yy为基本流的函数，所有非线性项置于右侧F_mn中。

对于扰动幅值较小的情况，方程(3)中非线性部分F_mn可以略去，得到线性化的扰动控制方程：

$ \begin{array}{l} \mathit{\Gamma }\frac{{\partial q'}}{{\partial t}} + A\frac{{\partial q'}}{{\partial x}} + B\frac{{\partial q'}}{{\partial y}} + Dq' + {E_{xx}}\frac{{{\partial ^2}q'}}{{\partial {x^2}}}\\ \;\;\;\; + {E_{xy}}\frac{{{\partial ^2}q'}}{{\partial x\partial y}} + {E_{yy}}\frac{{{\partial ^2}q'}}{{\partial {y^2}}} = 0 \end{array} $

(4)

将扰动改写为沿流向缓慢变化的特征函数和快速变化的波函数：

$ q'\left( {x,y,t} \right) = \hat \phi \left( {x,y} \right){{\rm{e}}^{i\left[ {\int_{{x_0}}^x {\alpha \left( {\bar x} \right){\rm{d}}\bar x - \omega t} } \right]}} + {\rm{c}}.{\rm{c}}. $

(5)

其中：$\hat \phi = \left( {\hat \rho , \hat u, \hat v, \hat T} \right)$为扰动波的形状函数，ω为扰动的频率，α=α_r+iα_i为复流向波数，c.c.表示共轭项。将式(5)带入式(4)中，得到线性稳定性方程：

$ \begin{array}{l} \hat A\frac{{\partial \hat \phi }}{{\partial x}} + \hat B\frac{{\partial \hat \phi }}{{\partial y}} + \hat D\hat \phi + {E_{xx}}\frac{{{\partial ^2}\hat \phi }}{{\partial {x^2}}}\\ \;\;\; + {E_{xy}}\frac{{{\partial ^2}\hat \phi }}{{\partial x\partial y}} + {E_{yy}}\frac{{{\partial ^2}\hat \phi }}{{\partial {y^2}}} = 0 \end{array} $

(6)

扰动特征函数、扰动波数以及扰动的增长率均沿流向缓慢变化。将方程抛物化，得到抛物化稳定性方程：

$ \hat A\frac{{\partial \hat \phi }}{{\partial x}} + \hat B\frac{{\partial \hat \phi }}{{\partial y}} + \hat D\hat \phi + {E_{yy}}\frac{{{\partial ^2}\hat \phi }}{{\partial {y^2}}} = 0 $

(7)

特征长度选用混合层入口处半涡厚度δ_ω0=[(U₁-U₂)/(|∂u/∂y|_max)]_x0，其他无量纲参数选取上层自由流参数。

1.2.2 差分格式与边界条件

在法向使用4阶精度的中心差分格式，同时采用网格拉伸，使得在剪切强烈的区域点数较密，流向采用一阶精度的向后差分格式。上下边界采用渐进边界条件^[17]：

$ \frac{{{\rm{d}}\hat \phi }}{{{\rm{d}}y}} = \mp q\hat \phi ,\;\;\;\;y \to \pm \infty $

(8)

其中：q为由线性稳定性理论获得的指数衰减系数，入口条件采用稳定性理论给出的不稳定波。

1.2.3 结果分析

入口所加入扰动波由稳定性理论得出，具体计算参数如表 2所示。

其中，α_r为流向波数的实部，扰动波相对上下层自由流的相对马赫数^[18]分别为：

$ M{a_{n,1}} = M{a_1}\left| {\omega /{\alpha _{\rm{r}}} - 1} \right| $

$ M{a_{n,2}} = M{a_1}\left| {\omega /{\alpha _{\rm{r}}} - {U_2}} \right| $

进一步将扰动波分为快模态和慢模态^[19]。当扰动相对于上层快速流为超声速，相对下层慢速流为亚声速时，即Ma_{n, 1}＞1、Ma_{n, 2}＜1, 称其为慢模态扰动；当扰动相对于上层快速流为亚声速，相对下层慢速流为超声速，即Ma_{n, 1}＜1、Ma_{n, 2}＞1，称其为快模态扰动，在表 2中分别以T10F、T10S代表快慢模态的演化。

图 3、图 4分别为稳定性理论计算得到的特征函数，图中用u_r、u_i、T_r、T_i、p_r、p_i来分别表示扰动速度、温度、压力的实部和虚部。图中显示，快模态特征函数在混合层上侧以指数形式迅速衰减，在混合层下侧呈现出震荡衰减；慢模态特征函数在混合层下侧以指数形式迅速衰减，在上侧呈现震荡衰减。

为了验证程序，当前用扰动方程计算了工况T10F，图 5给出扰动方程与PSE的计算结果对比，纵坐标代表扰动压力幅值演化，二者吻合的很好，表明PSE程序的正确性。

用PSE计算扰动波向下游的演化，图 6为近场扰动压力等值线。由近场压力等值线可以看出，两种模态扰动波均沿着固定方向传播，马赫波辐射具有很强的方向性。

2 超声速混合层远声场声辐射计算及结果分析 2.1 远声场积分理论介绍

Wu^[14]应用渐进匹配法与三层结构理论，建立了一种新的远声场计算方法，该方法理论推导出近场不稳定波幅值方程，进而外推出远场压力扰动的积分形式。对于单一不稳定波，远场的扰动压力有如下形式的解：

$ p = \left[ {{{\tilde p}_0}\left( {\bar \xi ,\tilde y} \right) + {R^{ - 1}}{{\tilde p}_1}\left( {\bar \xi ,\tilde y} \right) + \cdots } \right]{{\rm{e}}^{{\rm{i}}\alpha \zeta }} + {\rm{c}}.{\rm{c}}. $

(9)

$ \begin{array}{l} {{\tilde p}_0}\left( {\bar \xi ,\tilde y} \right) = \frac{{{{\rm{e}}^{ - \pi {\rm{i}}/4}}}}{{\sqrt {\tilde y} }}{\left( {\frac{{\alpha q_ \pm ^3}}{{2{\rm{ \mathsf{ π} }}Ma_{1,2}^2{c^2}}}} \right)^{0.5}} \cdot \\ \;\;\;\;\int_{ - \infty }^\infty A \left( \zeta \right)\exp \left\{ {\frac{{{\rm{i}}\alpha q_ \pm ^3}}{{2Ma_{1,2}^2{c^2}\tilde y}}{{(\bar \xi - \zeta )}^2}} \right\}{\rm{d}}\zeta \end{array} $

(10)

$ \begin{array}{l} \zeta = x - ct + {q_ \pm }y\\ {q_ \pm } = {\left[ {Ma_{1,2}^2{{\left( {{U_{1,2}} - c} \right)}^2} - 1} \right]^{0.5}} \end{array} $

(11)

$ \bar \xi = \bar x + q_ \pm ^{ - 1}\left[ {Ma_{1,2}^2{U_{1,2}}\left( {{U_{1,2}} - c} \right) - 1} \right]R{e^{ - 0.5}}y $

(12)

其中，ζ为马赫波的第一条特征线，代表马赫波辐射方向; ξ代表马赫波的包络，沿着特征线ξ，马赫波的包络$ {\tilde p_0}$为常数；$\tilde y = {{\mathop{Re}\nolimits} ^{ - 1}}y$;c为扰动相速度; U_1，2为上下层来流速度; Ma_1，2为上下层来流马赫数; α为流向波数; A(ζ)为流场中压力幅值。该方法优点在于不需要假设，可利用近场的压力扰动幅值直接积分获得远场的声压。远场马赫波辐射由式(9)计算得出，具体计算过程如下：

(1) 计算出近场扰动压力分布，得到其压力幅值A，文中由PSE计算得到。

(2) 利用积分公式(9)计算远场马赫波的压力扰动p。

2.2 不稳定波的马赫波辐射计算

分析工况T10F和T10S，将来自近场PSE的计算结果带入公式(9)中, 得到近场扰动压力p，如图 7所示。可以看到，无论从定性上还是定量上，这种积分方法的结果在近场区域与PSE结果(图 6)吻合很好，继续用该积分方法计算超声速混合层马赫波声辐射远场。

图 8给出这种积分方法计算工况T10F、T10S远场的结果。实线代表扰动能量传播方向，虚线代表扰动相位传播方向。可见，超声速平面混合层中，马赫波辐射具有很强的方向性，由图 8的压力等值线求得慢模态幅角为45°，快模态幅角为42°；马赫波能量沿特征线ξ辐射出去，相位沿特征线ζ传播。

3 上下侧来流温度比对气动声场的影响 3.1 上下侧来流温度不同时不稳定波增长特性分析

为了研究两侧来流温度对混合层声辐射特性的影响，选取2组不同温度比的工况与表 2的工况进行对比分析，具体参数见表 3。

图 9给出入口不稳定波增长率沿频率的变化。由图中可以看出: (1)在3种不同温度比的工况中，扰动均出现多模态；(2)温度比的变化对快慢模态影响不一致。快模态，温度比升高将导致不稳定波频带变宽，增长率变大，不稳定波的增长能力提高；慢模态，温度比升高则会导致不稳定波频带变窄，增长率变小，使不稳定波的增长能力下降；(3)温度比大于1时，快模态增长率大，增长频带宽，居主要地位；而在温度比小于1时，则是慢模态增长能力占据主要地位。

3.2 上下侧来流温度比对马赫波辐射的影响

为了说明来流温度对马赫波辐射的影响，图 10分别给出不同温度比下最不稳定快慢模态的压力幅值曲线及增长区间。从图中可以看出在温度比大于1的情况下，快模态初始增长率大，增长区间长，最终快模态幅值远大于慢模态；在温度比等于1的情况中，虽然初始快慢模态增长率相近，但由于快模态增长区间长，最终幅值也会大于慢模态；但在温度比小于1的情况中，最终慢模态压力幅值超越快模态，成为占据主导地位的扰动模态。

图 11和图 12分别给出工况T12F、T07F计算得到的远场扰动压力等值线。结合图 8可以看出，工况T12F在远场扰动压力大于工况T10F，而工况T07F在远场扰动压力远小于工况T10F。

结合不同温度比下，两种模态的幅值演化以及远场扰动压力的对比，可以认为，温度比大于1时，由于快模态增长快，范围广，在自然扰动的情况下，可以认为声辐射特性主要由快模态决定；相对，温度比小于1时，由于慢模态占优势，在自然扰动的情况下，可以认为声辐射特性要由慢模态决定。

4 结论

针对超声速平面混合层中马赫波辐射声场的问题，本研究采取PSE方程计算近场流动，而将不同于声比拟模型、基于稳定性理论和渐进展开的积分模型运用到远场声辐射计算。研究表明：

(1) 该方法能够准确地反映马赫波能量和相位传播的物理特征。对超声速平面混合层中存在的两种不同扰动模态，快模态和慢模态，计算给出两种模态扰动波的马赫波辐射方向角和声压能。

(2) 温度比是影响混合层稳定性、进而影响声辐射的重要参数。计算结果表明，对快模态，温度比升高导致不稳定波频率变宽，增长率变大，扰动波增长能力提高，加强了向远场辐射马赫波的能力；对慢模态，温度比升高导致不稳定波频率变窄，增长率变小，扰动波增长能力下降，减低了向远场辐射马赫波的能力。