杨亮，季振林

(哈尔滨工程大学 动力与能源工程学院，黑龙江 哈尔滨 150001)

穿孔消声器声学计算的快速多极混体边界元法

杨亮，季振林

(哈尔滨工程大学 动力与能源工程学院，黑龙江 哈尔滨 150001)

为提高计算效率以及扩展计算频率范围，本文将混体边界元方法与快速多极算法结合用于穿孔消声器的声学性能计算。通过与消声器传递损失测量结果和传统混体边界元方法计算结果的比较验证了快速多极混体边界元方法的正确性。与传统混体边界元方法相比，快速多极混体边界元方法在保证计算精度的同时，能够有效缩短大尺度问题的计算时间。将快速多极混体边界元方法应用于计算穿孔消声器的传递损失，结果表明，进出口管位置能够影响两通穿孔管消声器的高频消声特性，隔板穿孔能够改善三通穿孔管消声器的中频消声性能。

穿孔消声器； 传递损失； 混体边界元方法； 快速多极算法； 计算效率

穿孔消声器被广泛应用于机械设备的噪声控制。SELAMET等使用准一维方法对三通穿孔管消声器进行了研究并分析了结构参数变化对传递损失的影响[1]。SELAMET等应用解析方法预测了直通穿孔管阻性消声器的声学性能[2]。JI使用子结构多域边界元方法对直通穿孔管消声器及具有共振腔的穿孔管消声器进行了研究[3]，此后，多域边界元方法被用于阻性消声器的声学性能预测[4]。JI等使用同样的方法对三通穿孔管消声器进行了计算，研究了穿孔率、穿孔部分长度等对消声性能的影响[5]。FANG等[6]使用有限元方法对穿孔管消声器进行了模态分析。解析方法适用于规则结构消声器的声学性能预测，对实际应用中具有复杂内部结构的消声器无法应用，有限元方法和多域边界元方法在理论上适用于任意结构。但也存在一些缺点，在应用多域边界元方法时，必须按照边界条件将内部结构划分为许多子结构，再根据声学变量的连续性条件在边界上建立各子结构之间的联系，对于复杂的内部结构，前处理过程势必过于繁琐，进而限制了其应用范围。应用有限元方法进行穿孔消声器的计算已经十分成熟，相应的商业软件已经被广泛应用，但是对于穿孔面，由于有限元方法一般采用传递阻抗建立穿孔面内外两侧的关系，所以一般要求两个面的网格一一对应，对于具有多个穿孔面的消声器，前处理也会较为繁琐。WU等提出了一种混体边界元方法[7-8]，此方法适合具有复杂内部结构消声器的声学性能计算，不需要人为划分子结构，对于穿孔面，只需要建立一层面网格，因此非常适合具有较多穿孔结构消声器的传递损失计算。

然而，由于传统边界元方法自身的局限性，在处理大尺度和高频问题时对计算机内存的要求较高并会消耗较多的计算时间，近年来，快速多极算法的出现有效地克服了传统边界元方法的局限，应用快速多极边界元方法计算大尺度声学问题也取得了较大的进展[9-13]。

本文将混体边界元方法与快速多极算法结合用于穿孔消声器传递损失的计算，考察快速多极混体边界元方法的计算精度和计算效率，研究穿孔管位置对两通穿孔管消声器传递损失的影响，分析穿孔挡板对三通穿孔消声器声学性能的影响。

1 计算方法

1.1 混体边界元方法

直接混体边界元方法源于传统的子结构边界元方法，将具有复杂内部结构的消声器分成若干具有明确边界的几个部分，对每一部分使用边界积分方程描述，然后，根据各部分的交界面的声学参数连续条件，得到整体控制方程，交界面有两个未知量时，超奇异积分方程被用来提供额外的方程。对于如图1所示的直通穿孔管消声器，混体边界元方法可以表示为[7]

式中：P和Q分别为源点和场点；G=e-jk0r/4πrPQ，为格林函数；k为波数；z为空气中的特性阻抗；ζ为穿孔阻抗。边界条件R为常规边界条件，包括进口边界R1，出口边界R2以及刚性壁，T为内部薄壁边界，P为穿孔边界。这里只列出了本文涉及到的边界的直接混体边界元的数学表达式，完整的控制方程可以从文献[7]获得，此处不再赘述。进口给定振速v=1，出口给定阻抗边界条件，式(1)～(3)组成方程组：

(5)

其中

式中：E′只有源点P在穿孔边界时存在，对刚性壁边界其值为0，δ为Kronecker′s delta函数。求解式(5)，即可得到消声器进出口节点的声压值用于传递损失的计算，本文中，使用的单元类型为常单元。

图1 直通穿孔管消声器Fig.1 Straight-through perforated tube silencer

1.2 快速多极混体边界元方法

快速多极边界元中格林函数可以表示为[9]

(6)

其中，

(7)

(8)

(9)

将式(9)对P求偏导，得到

(10)

在远场计算中，源点与场点之间的影响系数存在多级传递关系，格林函数及其导数的多极子展开式可写成：

(11)

(12)

同理，将式(12)对P求偏导，得到

(13)

其中,

(14)

(15)

其中

(17)

(18)

(19)

(20)

图2 分级结构及通用交互组Fig.2 Hierarchical cell structure and common interaction cell

快速多极混体边界元方法的计算过程如下：

1)远场区域矩阵矢量积计算过程。

(21)

(22)

式中：Wn′n为插值方法(如拉格朗日法)的插值系数，Cml为组ml+1的父辈组所包含的子组集，此步计算将在向上传递的每一级中进行，直到级数为2。

(23)

式中：Fml为组ml的交互组集。需要注意的是，为了简便，在Tλmlλm′l的计算中，使用了通用交互组，其定义划分见图2，在迭代计算之前，根据格栅的位置关系，Tλmlλm′l已经计算完毕，迭代过程中只需根据距离关系调用。

(24)

可见，式(24)为递推公式，由于在1级与2级之间不存在向下传递关系，因而有

(25)

此步计算将在向下传递的每一级中进行，直到级数为L-1。

⑤在最低级L，计算每个节点i的远场影响系数φF,i，对应式(15)，计算公式为

(26)

对应式(16)，计算公式为

(27)

2)近场区域矩阵矢量积计算过程。

在最低级，若节点i所在组为mL，其近场组集为NmL，则近场矩阵矢量积即为计算节点i与所含节点间的影响系数，对奇异积分方程，计算公式为

(28)

对超奇异积分方程，计算公式为

(29)

最后，将近场影响系数与远场影响系数相加得到每个节点的影响系数，经过迭代计算即可得到消声器进出口节点的声压值。

2 算例验证及计算效率分析

考虑一个直通穿孔管消声器，对应图1中尺寸为：r1=0.024 5 m，r2=0.082 2 m，l=0.257 2 m，穿孔率Ф=0.08，孔径dh=0.002 49 m，壁厚tw=0.000 9 m。

穿孔阻抗表示为[14]

(30)

R和X的计算表达式：

(31)

(32)

式中：μ为空气动力粘性系数，α为单一孔模型的端部修正系数，基于活塞驱动模型

(33)

式中：b和h为相邻两孔在两个方向上的距离；ξ=dh/b，η=dh/h；式中第二个求和符号的上标“′”表示m和n不能同时为零，当m≠0且n≠0时，εmn=1，否则εmn=0.5；J1为第一类一阶贝塞尔函数。

使用传统混体边界元方法与快速多极混体边界元方法分别计算传递损失并与实验值比较得到结果如图3所示。可以看出，两种方法的计算结果与实验值吻合良好，验证了快速多极混体边界元方法的正确性。

图3 直通穿孔管消声器计算结果比较Fig.3 Transmission loss prediction of straight-through perforated silencer

下面从计算效率方面比较两种方法，选定频率为350Hz，在不同的网格数量的情况下比较计算时间，结果列于表1。所有计算均在Inteli7 3632QM, 2.2GHz的计算机上进行，从表中可以看出，在给定频率下，随着网格数量的增加，相比于传统混体边界元方法，快速多极混体边界元方法的优势越来越明显，所以，在实际应用中，可以根据研究问题的尺度大小选择适当的方法。

表1 直通穿孔管消声器计算时间比较

Table 1 Computation time comparison of straight through perforated tube silencer

编号网格数MBEM/s传统快速多极11448138427602220434613309342701402558744541626433921257542449361368569076—18975

3 计算结果分析

本节应用快速多极混体边界元方法研究穿孔管位置及穿孔挡板对消声器传递损失的影响。

首先，考虑如图4所示的两通穿孔管消声器，膨胀腔长度均为l=0.257 2 m，穿孔壁厚tw=0.000 9 m，穿孔率Ф=0.08，孔径dh=0.002 49 m，穿孔管半径r1=0.02 m，膨胀腔半径r2=0.1 m。对结构a，消声器进口管处于轴线处，考虑两种情况：a1出口管与进口管之间距离σ=0.06 m，a2出口管与进口管之间距离σ=0.05 m。对结构b，进出口穿孔管位于轴线两侧，考虑两种情况：b1两穿孔管距离轴线距离均为σ=0.06 m，b2两穿孔管距离轴线距离均为σ=0.04 m。比较穿孔管位置的变化对两种结构声学性能的影响，结果如图5所示。

可见，当进(出)口穿孔管处于轴线处时，传递损失曲线更接近轴对称的直通穿孔管消声器，消声效果也更好。穿孔管位置的改变对a、b两种结构的影响是不同的，对结构a，保持进口穿孔管处于轴线位置，改变处于非轴线位置的出口穿孔管的位置对传递损失有较大影响，在高频处更加明显。而对结构b，改变两个处于非轴线位置的穿孔管的位置对

传递损失结果影响很小。

图4 两通穿孔管消声器Fig.4 Two-pass perforated tube silencer

图5 两通穿孔管消声器传递损失计算结果Fig.5 TL comparison of two-pass perforated tube silencer

在实际应用中，常采用在消声器进出口添加刚性挡板构成三通穿孔管消声器的形式，如图6所示，穿孔管半径r1=0.02 m，膨胀腔半径r2=0.1 m，穿孔孔径dh=0.002 49 m，穿孔壁及挡板壁厚tw=0.000 9 m，穿孔率Ф=0.08，la=0.1 m，lp=0.3 m，lb=0.15 m，σ=0.06 m。分别使用混体边界元方法和快速多极混体边界元方法计算传递损失，结果如图7所示。同样选定频率为350 Hz，比较两种方法的计算时间列于表2。

图6 三通穿孔管消声器Fig.6 Three-pass perforated tube silencer

一般三通穿孔管消声器中的挡板是刚性薄壁，现在考虑将刚性薄壁结构替换为穿孔板，研究其对消声器声学性能的影响。在挡板分别为刚性壁面和穿孔板的情况下计算传递损失如图8所示。在中低频处，穿孔挡板三通穿孔消声器的有效消声频率更宽，性能更好，应用穿孔挡板代替刚性壁挡板是有意义的。

图7 三通穿孔管消声器计算结果比较Fig.7 Transmission loss prediction results of three-pass perforated tube silencer

Table 2 Computation time comparison of three-pass perforated tube silencers

编号网格数MBEM/s传统快速多极1551227933948227648460551399338864-18532

图8 穿孔挡板对消声器传递损失的影响Fig.8 Effect of perforated baffle on silencer TL

现在考虑穿孔挡板的物理参数变化对传递损失计算结果的影响，在穿孔率分别为0.03、0.08以及0.16的情况下计算传递损失，结果如图9所示。

挡板穿孔率的变化对消声器的传递损失有一定程度的影响，随着穿孔率的增大，传递损失曲线向高频方向移动，与此同时，在200～800 Hz的中频范围内，穿孔率的增大会改善消声器的声学性能。

下面考虑穿孔孔径大小对传递损失的影响，在孔径分别为0.002 49 m和0.004 98 m的情况下计算传递损失，结果如图10所示。增大穿孔孔径，会使传递损失曲线向低频方向移动。

图9 挡板穿孔率对三通穿孔消声器传递损失的影响Fig.9 Effect of baffle porosity on TL of three-pass perforated tube silencer

图10 挡板孔径对三通穿孔消声器传递损失的影响Fig.10 Effect of baffle hole diameter on TL of three-pass perforated tube silencer

在穿孔挡板的情况下，考虑三通穿孔管消声器左右两端膨胀腔长度对声学性能的影响，改变膨胀腔长度，计算结果如图11所示。由图可知，增大左右膨胀腔的长度会使传递损失曲线向低频方向移动。

图11 膨胀腔长度对三通穿孔消声器传递损失的影响Fig.11 Effect of inlet/outlet expansion chamber length on TL of three-pass perforated tube silencer

4 结论

1)快速多极混体边界元方法在能够保证计算精度的情况下，对大尺度高频问题能够有效地减少计算时间，在实际应用中，可以根据研究问题的情况选择适当的方法进行计算。

2)当进口位于轴线处时，消声器的声学性能更好，出口穿孔管的位置变化对计算结果的影响应该给予考虑，而当进口位于非轴线处时，进出口穿孔管距离的变化对传递损失影响很小，基本可以忽略不计。

3)使用穿孔板代替三通穿孔消声器中的刚性隔板能改善其中低频的声学性能，减小穿孔板穿孔孔径、增大穿孔率以及进出口膨胀腔的长度会使传递损失曲线向高频方向移动。

[1]SELAMET A, EASWARAN V, FALKOWSKI A G. Three-pass mufflers with uniform perforations[J]. The journal of the acoustical society of America, 1999, 105(3): 1548-1562.

[2]SELAMET A, XU M B, LEE I J, et al. Analytical approach for sound attenuation in perforated dissipative silencers with inlet/outlet extensions[J]. The journal of the acoustical society of America, 2005, 117(4): 2078-2089.

[3]JI Z L. Acoustic attenuation characteristics of straight-through perforated tube silencers and resonators[J]. Journal of computational acoustics, 2008, 16(03): 361-379.

[4]JI Z L. Boundary element acoustic analysis of hybrid expansion chamber silencers with perforated facing[J]. Engineering analysis with boundary elements, 2010, 34(7): 690-696.

[5]JI Z L, SELAMET A. Boundary element analysis of three-pass perforated duct mufflers[J]. Noise control engineering journal, 2000, 48(5): 151-156.

[6]FANG Z, JI Z L. Finite element analysis of transversal modes and acoustic attenuation characteristics of perforated tube silencers[J]. Noise control engineering journal, 2012, 60(3): 340-349.

[7]WU T W, CHENG C Y R, ZHANG P. A direct mixed-body boundary element method for packed silencers[J]. The Journal of the acoustical society of America, 2002, 111(6): 2566-2572.

[8]WU T W, WAN G C. Muffler performance studies using a direct mixed-body boundary element method and a three-point method for evaluating transmission loss[J]. Journal of vibration and acoustics, 1996, 118(3): 479-484.

[9]SAKUMA T, YASUDA Y. Fast multipole boundary element method for large-scale steady-state sound field analysis. Part I: setup and validation[J]. Acta acustica united with acustica, 2002, 88(4): 513-525.

[10]WU H, LIU Y, JIANG W. A low-frequency fast multipole boundary element method based on analytical integration of the hypersingular integral for 3D acoustic problems[J].Engineering analysis with boundary elements, 2013, 37(2): 309-318.

[11]LI S, HUANG Q. A fast multipole boundary element method based on the improved Burton-Miller formulation for three-dimensional acoustic problems[J]. Engineering analysis with boundary elements, 2011, 35(5): 719-728.

[12]王雪仁, 季振林. 快速多极子声学边界元法及其应用研究[J]. 哈尔滨工程大学学报, 2007, 28(7): 752-757. WANG Xueren, JI Zhenlin. Fast multipole acoustic BEM and its application[J]. Journal of Harbin Engineering University, 2007, 28(7): 752-757.

[13]吴海军, 蒋伟康, 鲁文波. 三维声学多层快速多极子边界元及其应用[J]. 物理学报, 2012, 61(5): 54301-054301. WU Haijun, JIANG Weikang, LU Wenbo. Multilevel fast multipole boundary element method for 3D acoustic problems and its applications[J]. Acta physica sinica, 2012, 61(5): 54301-054301.

[14]JI Z L, FANG Z. On the acoustic impedance of perforates and its application to acoustic attenuation predictions for perforated tube silencers[C]//Proceedings of INTER-NOISE. Osaka, 2011(5): 2803-2813.

本文引用格式：

杨亮，季振林. 穿孔消声器声学计算的快速多极混体边界元法[J]. 哈尔滨工程大学学报， 2017， 38(8): 1247-1253.

YANG Liang, JI Zhenlin. Acoustic computation of perforated silencers by fast multi-pole mixed-body boundary element method[J]. Journal of Harbin Engineering University， 2017， 38(8): 1247-1253.

Acoustic computation of perforated silencers by fast multi-pole mixed-body boundary element method

YANG Liang, JI Zhenlin

(School of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China)

To improve computation efficiency and extend the computation frequency range, the mixed-body boundary element method (MBEM) was combined with the fast multi-pole algorithm to evaluate the acoustic attenuation performance of perforated silencers. The correction of fast MBEM was validated by comparing the measuring result of transmission loss with the computation result of traditional MBEM. The fast multi-pole mixed-body boundary element method (FMMBEM) may efficiently save computational time for large-scale acoustic problems without affecting accuracy. The FMMBEM was then employed to predict the transmission loss of perforated silencers. Predicted results showed that the locations of inlet and outlet tubes may affect the acoustic attenuation characteristics of the two-pass perforated tube silencer at higher frequencies, while the perforation on bulkhead(s) may improve the acoustic attenuation performance of three-pass perforated tube silencer in the middle frequency range.

perforated silencer; transmission loss; mixed body boundary element method; fast multi pole algorithm; computation efficiency

2016-06-14.

日期：2017-04-26.

国家自然科学基金项目(11174065).

杨亮(1989-), 男, 博士研究生； 季振林(1965-), 男, 教授, 博生生导师.

季振林, E-mail:jizhenlin@hrbeu.edu.cn.

10.11990/jheu.201606040

TB535

A

1006-7043(2017)08-1247-07

网络出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20170426.1801.066.html