基于PE算法的管道噪声声传播预测技术

Prediction of Propagation of Pipeline Noise Based on Parabolic Equation Algorithm

摘要: 本文针对风扇外涵道声传播数值预测问题，基于PE方法，开发了声阻抗边界条件模型。针对二维与三维的声传播问题、多种声学边界模型，应用有限差分法进行离散，采用四阶龙格库塔法进行推进，开发了管道噪声声传播预测软件。标模算例验证结果表明，此软件与传统声学有限元软件相比，可以有效缩短风扇外涵道声传播数值预测时间，其既能保证足够的预测精度，又提高了计算效率。

Abstract: In this paper, based on the PE method, the acoustic impedance boundary condition model is developed for the numerical prediction of the noise propagation in the duct outside the fan. Aiming at the two-dimensional and three-dimensional sound propagation problems and various acoustic boundary models, the finite difference method is used for discretization, and the fourth-order Runge-Kutta method is used for advancement, and the pipeline noise sound propagation prediction software is developed. The verification results of the standard model calculation example show that compared with the traditional acoustic finite element software, this software can effectively shorten the numerical prediction time of the ducted sound propagation outside the fan. It can not only ensure sufficient prediction accuracy, but also improve the calculation efficiency.

1. 引言

常用的风扇气动噪声控制手段包括控制风扇气动噪声声源的低噪声风扇设计和控制风扇噪声通过风扇机匣、进气道和外涵道向外传播的消声声衬设计，其中消声声衬的应用可以获得非常显著的降噪效果，也是当代民用客机最常用的被动式消声设计方法。民用航空发动机进气道、外涵道的声衬优化设计需要建立考虑三维复杂管道声传播、声阻抗边界、声源边界的声传播模型，并在此基础上进行快速多次优化迭代设计，因此满足工程精度要求且计算效率足够高的噪声声传播预测模型是工程应用中开展风扇气动噪声评估与声衬结构优化设计的关键。

常用的噪声声传播预测技术包括传递单元法、声线方法、有限元法和CAA方法。基于解析预测模型的传递单元法在计算效率和预测精度上可以满足工程需要，但是只能开展二维声传播分析；基于声线理论的声传播计算虽然在计算效率上可以满足工程需要，但是预测精度不适用与外涵道声传播的预测；基于数值方法的有限元法和CAA方法虽然计算精度满足工程需求，也能够考虑三维声传播问题，但是计算效率低，需要大量的计算资源。因此需要既能考虑三维声传播问题，又能满足工程精度和计算效率的外涵道声传播预测技术。

追求精确和高效的方法来预测复杂的三维发动机风扇外涵道的声传播已经有很长的历史，从最初一般高效的近似方法到高保证，高计算强度的椭圆型方程数值求解方法 [1] - [10]。

国内方面，吴鑫波等人 [11]，利用噪声源计算软件，基于抛物方程(PE)方法，针对风电场，分析声传播过程中的衰减特性，结果表明抛物方程算法在模拟风电场的噪声传播分布时具备较高的准确性和可靠性。

朱凤芹等人 [12] 研究了海洋暖涡对声传播的影响。主要是利用涡旋动力学模型求解流函数，结合PE方法计算中尺度暖涡对声传播的影响，得到声源处于不同位置处声传播效应。

上述研究要么计算量很大，要么只研究了开阔空间的声传播问题，而本文则基于PE方法研究了声衬结构的管道中声传播算法，针对内流情况下的声传播问题，在计算速度远大于数值模拟的同时，还能保证较高的计算精度。

本研究主要通过建立考虑外涵道单方向传播的噪声声传播物理机制的线化欧拉方程，经过数学变换，构建PE方法的噪声声传播预测模型的控制方程。基于控制方程，开发了管道壁面的声阻抗边界条件模型。针对二维与三维声传播问题，应用有限差分法进行离散，采用四阶龙格库塔法进行推进，开发了管道噪声声传播预测软件。通过标模算例的对比测试，验证了采用抛物线方程形式的声传播数值方法可以有效缩短风扇外涵道声传播数值预测时间。

2. 理论模型

本文从噪声声传播的线化欧拉方程出发，通过数学变换，建立了噪声声传播预测模型的控制方程。通过将线化欧拉方程的偏微分形式简化为抛物线方程的常微分形式，可以把双向声传播耦合问题近似成单向声传播问题。通过建立控制方程下声源边界条件模型、管道口声学边界条件模型与声衬阻抗边界模型，结合频域空间离散和推进格式数值算法，开发了管道噪声声传播预测软件。

与以往的研究相比，本文中讨论的提法有几个独特的方面。首先，根据流体动力学原理，系统地推导了控制方程和边界条件的公式。这些推导明确地揭示了该方法所涉及的近似和假设，为该方法今后的改进奠定了基础。这些近似和假设可以针对性放宽，以提高控制方程的准确性和适用性。其次，利用泰勒展开法和Pade逼近法，对抛物型方程进行系统渐近展开，得到精度更高的抛物型方程。

2.1. 抛物线形式的控制方程

对于外涵道声传播问题，考虑到流道方向上是声音的主要传播方向，而环境与流动变量只是缓慢变化。频域下的波动方程如下：

$\frac{{\partial }^{2}p}{\partial x_1^2}+\frac{2i{k}_{0}M}{{\beta }^2}\frac{\partial p}{\partial x_1}+\frac{k_0^{2}p+{\nabla }_{\alpha }^{2}p}{{\beta }^2}=0$ (1)

p表示为频域上的声压， $k_0$ 为声学波数。其中，横截面上的拉普拉斯算子写成

${\nabla }_{\alpha }^2=\frac{{\partial }^{2}p}{\partial x_2^2}+\frac{{\partial }^{2}p}{\partial x_3^2}$ (2)

针对全三维外涵管道，声音沿着弯曲的外涵道传播，因此控制方程可以在正交曲线坐标系中求解方程组，设 ${\xi }_1$ 平行于外涵道流动坐标， ${\xi }_2,{\xi }_3$ 为管道截面坐标

$\xi =\left\{{\xi }_1,{\xi }_2,{\xi }_3\right\}$ (3)

左边变换系数，

$h_i={\left({{\displaystyle \underset{j=1}{\overset{3}{\sum }}\left(\frac{\partial x_j}{\partial {\xi }_i}\right)}}^2\right)}^{1/2}$ (4)

在正交曲线坐标系中的拉普拉斯算子可写成

${\nabla }^2=\frac{1}{h_1{h}_2{h}_3}\frac{\partial }{\partial {\xi }_1}\left(\frac{h_2{h}_3}{h_1}\frac{\partial }{\partial {\xi }_1}\right)+{\nabla }_{\alpha }^2$ (5)

因此，对于在正交曲线坐标系下的风扇噪声的传播问题，控制方程为

$\left\{P+\frac{h_2{h}_{3}M\left(i{k}_0-\mu \right)}{{\beta }^2}+\frac{i{h}_2{h}_3\lambda }{{\beta }^2}\right\}p=0$ (6)

其中， $P=\frac{h_2{h}_3}{h_1}\frac{\partial }{\partial {\xi }_1}$, ${\beta }^2=1-M^2$。

l定义为：

$\lambda ={\left\{k_0^2+M^2\mu \left(i{k}_0-0.25\mu \right)+{\beta }^2{\nabla }_{\alpha }^2\right\}}^{1/2}$ (7)

2.2. 逐次逼近方法

抛物型方程只是一个象征性的方程，无法求解，因为横截面上的微分是由参量l定义的平方根函数。平方根函数需要转换成允许对声压进行微分的形式。为此，参量l可以写成

$\lambda ={\lambda }_0{\left(1+q\right)}^{1/2}$ (8)

根据泰勒级数展开，控制方程可以转化为：

$\left\{P+\frac{h_2{h}_{3}M\left(i{k}_0-\mu \right)}{{\beta }^2}+\frac{i{h}_2{h}_3{\lambda }_0}{{\beta }^2}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}{a}_n{q}^n}\right\}p=0$ (9)

这里可以写成另一种形式

$\frac{\partial p}{\partial {\xi }_1}=-\frac{h_1}{{\beta }^2}\left\{M\left(i{k}_0-\mu \right)+i{\lambda }_0+i{\lambda }_0{\displaystyle \underset{n=1}{\overset{N}{\sum }}{a}_n{q}^n}\right\}p$ (10)

这是一个轴向坐标下的一阶常微分方程，所有微分运算都是显式的，可以用标准的数值方法求解。

2.3. 阻抗边界

控制方程在流道内部有效，并且根据管道壁的边界条件进行求解。对于硬壁面，边界条件仅仅是壁上法向速度为零。对于带有声衬的情况在墙壁上，边界条件考虑了衬板的位移，因此动态条件和运动条件都包括在内。前者是在声压作用下的衬管，将衬管上的法向速度与管内的压力相关壁上的流体流动，后者将壁描述为连续的并且没有流动分离。通常定义动态条件，也称为阻抗边界条件并用于频域或单个频率，因为衬管的阻抗为通常指定为频率的函数。

这些可以结合成单一的方程，通过消除表面位移，得到

$\left(i\omega -U\frac{\partial }{\partial x_1}\right)p\mp i\omega Z{u}_n=0$ (11)

它只涉及声学参量，可用于控制方程。结合动量方程，声速可以用声压代替，阻抗边界条件呈如下形式：

${\left(i\omega -U\frac{\partial }{\partial x_1}\right)}^{2}p\mp \frac{i\omega Z}{{\rho }_0}\frac{\partial p}{\partial x_n}=0$ (12)

2.4. 离散与推进方法

截面离散的方法是有限差分格式，其使用泰勒展开来表示导数。

本软件用的是四阶龙格库塔法进行时间和空间的推进计算，其精度高，易于实现。它将解从起始位置匹配到管道中的任何位置。其必须指定描述后风机管道入口噪声源的初始声压。因此噪声源必须使用类似于平面波源，高斯分布，或一系列的管道声模态的近似方法。

2.5. 程序设计

基于控制方程开发了噪声预测程序，程序核心功能是计算飞机发动机后风扇噪声。传播过程从后风扇管道入口处的噪声源开始，通过风扇管道，在风扇管道中声波与管道几何形状和声衬(如果存在的话)相互作用。主要考虑发动机结构、声衬以及风扇管道流场的影响。

软件的核心任务是或者噪声传播和衰减计算数据。噪声传播包括定义噪声源，声传播的网格，平均流量，离散化，线性处理和流动方向的推进算法。

3. 程序验证

本文通过已有的公开文献算例、试验数据与标模算例，进行模型算法的验证。

3.1. 在2000 HZ下带声衬的矩形管道

通过数值计算与试验数据对比，验证算例来自NASA Langley掠流阻抗管道，其通常用于评估声学在受控的流动和声学环境中的声衬样本 [13]。试件截面尺寸分别为50.8 mm和50.8 mm，试件长度为406.4 mm。马赫数可达0.6，可用频率范围为0.3~3.0 kHz。试件前缘的声压级可达155 dB。

表1为该测试的总体参数，其中包括设备的几何参数和测试条件。假定声速为340.3 m/s。频率此处给出的示例分别为2000 Hz和3000 Hz。将测试部分离散化为均匀间隔的三个尺寸网格。x₁ = 0时的初始声源是一个平面波，其振幅设置为在第一个遍历麦克风位置测量的值。

Table 1. Parameters for GFIT test

表1. GFIT测试参数

图1给出了处理后的无平均流动的管道壁面声压的数值计算，其中分别给出了声压的实部和虚部。

Figure 1. Sound pressure on GFIT treated wall duct at 2000 Hz

图1. GFIT处理的壁管在2000 Hz处的声压

对于这些情况，将声压级与图2中的实验数据进行比较，其中声压级与实心曲线的轴向距离成函数关系代表计算结果，符号代表测试数据，两者之间的整体误差小于5%。

Figure 2. Comparison of SPL for treated wall duct at 2000 Hz

图2. 在2000 Hz下带声衬的壁管的SPL对比

这种情况下的传递损失如图3所示，其中以dB为单位绘制了传递损失作为轴向距离的函数。

Figure 3. Transmission Loss of SPL for treated wall duct at 2000 Hz

图3. 在2000 Hz下带声衬的壁管的传递损失

下图4为3000 Hz情况下，声压级与实验数据的对比，两者趋势基本相同，整体误差小于10%。

3000 Hz情况下的传递损失如图5所示，其中以dB为单位绘制了传递损失作为轴向距离的函数。

通过上述与NASA标模数据进行对比的算例，在2000 Hz与3000 Hz情况下均与实验数据误差控制在10%以内，表明软件计算结果具备很高的精度。

Figure 4. Comparison of SPL for treated wall duct at 3000 Hz

图4. 在3000 Hz下带声衬的壁管的SPL对比

Figure 5. Transmission Loss of SPL for treated wall duct at 3000 Hz

图5. 在3000 Hz下带声衬的壁管的传递损失

3.2. 具有均匀流动的环形管道

环形管道中具有均匀流动的验证案例如表2，环形管道的几何形状为横截面恒定的管道，内半径为0.25 m，外半径为0.5 m。管道的长度为1 m。情况下的频率选择为1000 Hz，此时初始声压由第一圆周模式和第四径向模式定义。流动马赫数为0.2，由于其横截面恒定，因此在管道中是均匀的。

Table 2. Parameters for mode in annular duct with static medium

表2. 带有静态介质的环形管道中的模式参数

如图6所示，分别压力的实部和虚部。

Figure 6. Source pressure for mode in annular duct

图6. 环形管道中模式的声源压力

管道外壁上数值计算的声压如图7所示。

Figure 7. Results of mode propagation on duct outer surface at M = 0.2

图7. M = 0.2时管道外表面模式传播结果

管道内壁数值计算的声压如图8所示。

Figure 8. Results of mode propagation on duct inner surface at M = 0.2

图8. M = 0.2处的管道内表面模式传播结果

在第一和最后一个横截面的声压如图9。

Figure 9. Results of mode propagation on duct cross sections

图9. 管道横截面模式传播的结果

上述结果与预期结果相符，其优点是计算效率高，在1min以内完成计算，这是传统数值计算软件无法到达的效果。

4. 结论

本文从流体动力学控制方程出发，建立噪声声传播线化欧拉方程，通过数学变换，建立噪声声传播预测模型的控制方程。通过PE算法，把双向声传播耦合问题近似成单向声传播问题，并基于此模型开发了管道声传播预测工具。

通过与已有的公开文献算例、试验数据与标模算例，进行模型算法的验证。结果表明，软件计算具备传统声传播计算软件不具备的计算效率，同时计算结果具有较好的精度，此软件可以用来开展风扇气动噪声评估和声衬结构优化设计。

基金项目

本文来自上海市科学技术委员会科研计划项目，基金号：19511104100。