楊　亮， 季振林， Wu T W

(1.哈爾濱工程大學(xué) 動力與能源工程學(xué)院，哈爾濱　150001；2.肯塔基大學(xué) 機(jī)械工程學(xué)院，列克星敦　40506)

?

消聲器傳遞損失預(yù)測的邊界元數(shù)值配點混合方法

楊亮1， 季振林1， Wu T W2

(1.哈爾濱工程大學(xué) 動力與能源工程學(xué)院，哈爾濱150001；2.肯塔基大學(xué) 機(jī)械工程學(xué)院，列克星敦40506)

摘要：將邊界元法與數(shù)值配點法結(jié)合形成混合方法用于計算任意截面形狀消聲器的傳遞損失。消聲器劃分為若干子結(jié)構(gòu)，用邊界元法計算具有非規(guī)則形狀的子結(jié)構(gòu)阻抗矩陣，用二維有限元法提取等截面子結(jié)構(gòu)特征值及特征向量，用配點法獲得阻抗矩陣；將每個子結(jié)構(gòu)阻抗矩陣連接用于傳遞損失計算。為減少計算時間提出簡化方法計算消聲器傳遞損失。結(jié)果表明，混合法在保證計算精度前提下可節(jié)省計算時間。

關(guān)鍵詞：消聲器；傳遞損失；邊界元方法；數(shù)值配點方法

解析方法、數(shù)值方法廣泛用于消聲器傳遞損失計算。前者雖計算速度快、精度高，但只能用于規(guī)則形狀消聲器計算，Selamet等[1-3]用解析模態(tài)匹配法計算膨脹腔、內(nèi)插管膨脹腔及回流腔等消聲器傳遞損失。Selamet等[4]用二維解析模態(tài)匹配法計算填充吸聲材料的阻性消聲器聲學(xué)性能。數(shù)值方法適用性強(qiáng)，理論上可用于任意復(fù)雜結(jié)構(gòu)消聲器計算，Peat等[5]用有限元、Ji[6]用邊界元、Wu等[7-8]用直接混體邊界元法對不同類型消聲器進(jìn)行廣泛研究。其中，直接混體邊界元法更具靈活性，尤其適用消聲器的聲學(xué)性能計算。對大尺度問題，數(shù)值方法計算時間長、內(nèi)存需求大。因此，考慮將解析法與數(shù)值法結(jié)合用于消聲器聲學(xué)性能研究具有實際意義。Kirby等[9-12]用二維有限元法獲得消聲器或空調(diào)管道截面的特征值及特征函數(shù)，用配點法獲得相關(guān)聲學(xué)性能評價參數(shù)，并用數(shù)值模態(tài)匹配法對阻性消聲器進(jìn)行研究。Fang等[13-14]用數(shù)值模態(tài)匹配法對抗性消聲器及有流穿孔阻性消聲器進(jìn)行計算。

無論數(shù)值配點法或數(shù)值模態(tài)匹配法，對沿軸向方向消聲器或消聲管道截面形狀變化情況均不適用，因此，不同計算方法的耦合及優(yōu)化需進(jìn)一步研究。

本文擬用邊界元-數(shù)值配點混合法進(jìn)行消聲器傳遞損失計算。將消聲器劃分為若干子結(jié)構(gòu)，對沿軸向管道截面變化或具有復(fù)雜內(nèi)部結(jié)構(gòu)的子結(jié)構(gòu)用直接混體邊界元法計算其阻抗矩陣，考慮三維邊界元法與二維有限元法網(wǎng)格劃分的一致性，對截面均勻一致的子結(jié)構(gòu)考慮用二維有限元方法進(jìn)行特征值分析，獲得橫向波數(shù)，進(jìn)而用配點法計算獲得阻抗矩陣，將不同子結(jié)構(gòu)阻抗矩陣連接獲得消聲器的傳遞損失。為減少計算時間，用簡化方法壓縮子結(jié)構(gòu)阻抗矩陣，并通過算例驗證混合方法的有效性。

1直接混體邊界元方法

該方法源于傳統(tǒng)子結(jié)構(gòu)邊界元方法，即將具有復(fù)雜內(nèi)部結(jié)構(gòu)的消聲器分成若干具有明確邊界的幾部分，對每部分用邊界積分方程描述，據(jù)各部分交界面的聲學(xué)參數(shù)連續(xù)條件獲得整體控制方程。交界面有兩個未知量時，用超奇異積分方程提供額外方程。本文涉及的邊界積分方程(完整積分見文獻(xiàn)[8])為

(1)

(2)

圖1　消聲器子結(jié)構(gòu)劃分Fig.1 Substructure division of silencer

源點、場點距離較近時，對公式中涉及的奇異積分及超奇異積分需妥善處理[15-16]。不同于文獻(xiàn)[7-8]求解矩陣的直接方法，本文用帶預(yù)處理的BI-CG求解器。為描述問題，將消聲器分三部分(圖1)，對子結(jié)構(gòu)A、C，考慮用邊界元方法獲得阻抗矩陣；對子結(jié)構(gòu)A，阻抗矩陣可表示為

(3)

式中：PAI，vAI分別為進(jìn)口聲壓及質(zhì)點振速向量；PAO，vAO分別為出口聲壓及質(zhì)點振速向量。

阻抗矩陣作為結(jié)構(gòu)聲學(xué)特性描述，是唯一確定的。對式(3)需在不同邊界條件下用BEM計算。如，為獲得阻抗矩陣第一列數(shù)值，對進(jìn)口第一單元設(shè)v=1，進(jìn)出口其它單元設(shè)v=0，進(jìn)出口聲壓值即為第一列阻抗矩陣值。盡管需多次運(yùn)行BEM，但每次求解的矩陣相同。對結(jié)構(gòu)C，得阻抗矩陣為

(4)

2數(shù)值配點法

圖1的子結(jié)構(gòu)B沿軸向等截面，但由于截面形狀任意，無法用解析法獲得橫向波數(shù)，故用數(shù)值方法進(jìn)行特征值分析，考慮二維有限元及三維邊界元網(wǎng)格劃分的一致性，為便于各子結(jié)構(gòu)阻抗矩陣連接，用二維有限元計算特征值及特征向量。以阻性穿孔管道為例，圖1中Ⅰ的內(nèi)部為空氣，Ⅱ為填充吸聲材料，Ⅰ與Ⅱ之間由穿孔管分隔。區(qū)域Ⅰ、Ⅱ內(nèi)橫向聲壓控制方程可表示為

(5)

(6)

式中：pxy1，pxy2分別為區(qū)域Ⅰ、Ⅱ的橫向聲壓分量；kxy1，kxy2分別為空氣、吸聲材料中橫向波數(shù)，兩區(qū)域軸向波數(shù)kz相同，與橫向波數(shù)滿足關(guān)系為

(7)

(8)

穿孔管內(nèi)外滿足邊界條件為

(9)

用伽遼金加權(quán)余量法并代入邊界條件獲得阻性穿孔管道橫截面的本征方程[17]為

(10)

(11)

(12)

求解式(10)可獲得軸向波數(shù)及本征向量Φn(x,y)。對穿孔管抗性消聲器，式(10)簡化[18]為

(13)

對剛性壁直管道，橫截面本征方程表達(dá)式為

(14)

管道內(nèi)聲壓、質(zhì)點振速表示為

(15)

(16)

(17)

在B的進(jìn)出口各選N個配點，則有N個模態(tài)，據(jù)式(16)，獲得B進(jìn)出口配點振速為

(18)

同理，據(jù)式(15)得進(jìn)出口聲壓表達(dá)式為

(19)

3消聲器傳遞損失計算

若使子結(jié)構(gòu)B進(jìn)出口配點分別與子結(jié)構(gòu)A的出口截面、子結(jié)構(gòu)C的進(jìn)口截面節(jié)點一致，可連接二者阻抗矩陣獲得整體阻抗矩陣用于傳遞損失計算。對圖1消聲器，按結(jié)構(gòu)連接順序，建立前兩子結(jié)構(gòu)A、B進(jìn)出口阻抗矩陣關(guān)系，即

(20)

式中：

用相同方法，將式(20)中阻抗矩陣與子結(jié)構(gòu)C的阻抗矩陣連接可得整個結(jié)構(gòu)進(jìn)出口的阻抗矩陣關(guān)系，即

(21)

對含更多子結(jié)構(gòu)的消聲器，阻抗矩陣計算方法相同。在消聲器進(jìn)出口，聲波可視為平面波，此時式(21)被壓縮為2×2的矩陣。消聲器四極參數(shù)表示為

(22)

傳遞損失計算公式可表示為

(23)

式中：SI，SO分別為消聲器進(jìn)出口面積。

4算例驗證

以一方形膨脹腔消聲器為例驗證本文方法的正確性，見圖2。將消聲器劃分為3個子結(jié)構(gòu)，lA=lC=0.05 m，lB=0.125 m，進(jìn)出口管道長寬均b=0.05 m，膨脹腔長寬均bc=0.15 m。子結(jié)構(gòu)A、C的阻抗矩陣由BEM獲得，子結(jié)構(gòu)B的阻抗矩陣由配點法計算，各子結(jié)構(gòu)交界面均生成144個單元，即可考慮144個模態(tài)，計算結(jié)果見圖3。由圖3看出，混合方法計算結(jié)果與三維邊界元吻合較好，該方法的正確性獲得驗證。用數(shù)值配點方法計算阻抗矩陣，只需在子結(jié)構(gòu)進(jìn)出口劃分網(wǎng)格，可節(jié)省計算時間，該時間與子結(jié)構(gòu)長度無關(guān)，適用于具有較長等截面結(jié)構(gòu)的大尺度聲學(xué)問題。需注意的是，由于橫向波數(shù)及本征向量由二維有限元方法獲得，較解析方法相比存在一定誤差，且隨模態(tài)階數(shù)提高誤差有增大趨勢。據(jù)文獻(xiàn)[9],本算例認(rèn)為最高到第29個模態(tài)結(jié)果較準(zhǔn)確，其后高階模態(tài)對應(yīng)的特征值、特征向量與準(zhǔn)確值誤差較大。因高階模態(tài)對結(jié)果貢獻(xiàn)量小，故傳遞損失計算結(jié)果相對準(zhǔn)確，解析法與數(shù)值法對前5個高階模態(tài)計算結(jié)果比較見表1。

圖2　方形膨脹腔消聲器Fig.2 Square expansion chamber silencer

圖3　方形膨脹腔消聲器傳遞損失Fig.3 TL of square expansion chamber silencer

模態(tài)解析解數(shù)值解相對誤差/%(1,0)20.94421.0040.29(1,1)29.61929.7040.29(2,0)41.88842.3681.15(1,2)46.83247.2880.97(2,2)59.23859.9171.15

若模型尺寸較大，為保證與子結(jié)構(gòu)A、C截面網(wǎng)格一致，B截面會劃分較多網(wǎng)格，對式(18)、(19)，求解模態(tài)幅值系數(shù)會費較多計算時間，實際應(yīng)用中，考慮過多模態(tài)無意義。本算例計算頻率最高到2 500 Hz，(3,0)模態(tài)或(0,3)模態(tài)激發(fā)頻率約3 440 Hz，實際計算中考慮最高到第二階模態(tài)即可滿足需要。

因此，如何保證計算精度的同時控制模態(tài)數(shù)量尚待進(jìn)一步研究。本文考慮壓縮子結(jié)構(gòu)A、C的阻抗矩陣，在A的出口截面及C的進(jìn)口截面劃分若干子區(qū)域，用區(qū)域內(nèi)部節(jié)點平均聲學(xué)變量近似表示整個區(qū)域聲學(xué)變量數(shù)值，對子結(jié)構(gòu)B的進(jìn)出口截面，選擇與子結(jié)構(gòu)A、C劃分子區(qū)域中心對應(yīng)的配點，進(jìn)而用數(shù)值配點法計算阻抗矩陣。子結(jié)構(gòu)A進(jìn)口及C出口認(rèn)為平面波傳播，只劃分1個子區(qū)域。

以方形穿孔管抗性消聲器為例，驗證簡化方法的正確性。圖4的結(jié)構(gòu)尺寸與圖2相同，區(qū)別為膨脹腔內(nèi)部添加穿孔管，穿孔率φ=0.08，穿孔孔徑dh=0.002 49 m，穿孔壁厚tw=0.000 9 m。子結(jié)構(gòu)交界面區(qū)域簡化劃分見圖5，各子結(jié)構(gòu)交界面劃分為9個子區(qū)域，bs= 1/3bc=0.05 m。簡化后式(3)、(4)中矩陣規(guī)模減小為10×10。對子結(jié)構(gòu)B，選擇與子區(qū)域1～9中心一致的配點，阻抗矩陣計算同前，式(17)中矩陣規(guī)模減小為18×18，大大減少計算時間。傳遞損失計算結(jié)果見圖6，可見與三維邊界元吻合良好，驗證所提方法實際應(yīng)用的有效性。

圖4　方形穿孔管消聲器結(jié)構(gòu)尺寸Fig.4 Structure dimension of square perforated tube silencer

圖5　子結(jié)構(gòu)A、C截面劃分Fig.5 Cross-section division of substructures A and C

圖6　方形穿孔管消聲器傳遞損失Fig.6 TL of square perforated tube silencer

通過兩個算例驗證所提方法的有效性，并與傳統(tǒng)子結(jié)構(gòu)方法計算時間比較，說明本文方法的優(yōu)越性?？紤]穿孔管阻性消聲器，模型結(jié)構(gòu)尺寸同圖4，區(qū)別為膨脹腔內(nèi)填充吸聲材料，流阻率R=4 896 Rayls/m。子結(jié)構(gòu)劃分形式與交界面子區(qū)域劃分同圖5，簡化的混合方法與三維邊界元計算結(jié)果見圖7，可見吻合良好。

圖7　方形穿孔管阻性消聲器傳遞損失Fig.7 TL of square dissipative perforated tube silencer

橢圓截面消聲器亦廣泛應(yīng)用，故以橢圓膨脹腔消聲器為例說明簡化混合方法的應(yīng)用。模型尺寸見圖8，lA=lC=0.05 m，lB=0.182 3 m，橢圓長軸b1=0.2 m，短軸b2=0.14 m，進(jìn)出口插管長l1=l2=0.025 m，消聲器進(jìn)出口管道截面為圓形，半徑r=0.025 m。各子結(jié)構(gòu)交界面子區(qū)域劃分見圖9，傳遞損失計算結(jié)果見圖10，可見計算結(jié)果與三維邊界元吻合良好。

圖8　橢圓膨脹腔消聲器結(jié)構(gòu)尺寸Fig.8 Structure dimension of elliptical expansion chamber silencer

圖9　橢圓截面子區(qū)域劃分Fig.9 Cross-section division of elliptical interface

圖10　橢圓膨脹腔消聲器傳遞損失Fig.10 TL of elliptical expansion chamber silencer

簡化混合方法與傳統(tǒng)子結(jié)構(gòu)邊界元方法在頻率1 000 Hz時計算時間比較見表2。以方形穿孔管阻性消聲器與橢圓膨脹腔消聲器為例，所有計算均在Intel i7 3632QM, 2.2GHz計算機(jī)上進(jìn)行。由于用簡化混合方法計算B的阻抗矩陣，整個子結(jié)構(gòu)無需劃分網(wǎng)格，只選適當(dāng)?shù)呐潼c，每次計算求解的矩陣規(guī)模遠(yuǎn)小于傳統(tǒng)邊界元方法，子結(jié)構(gòu)B的阻抗矩陣計算時間基本可忽略不計，因此簡化混合方法能節(jié)省大量機(jī)時。實際應(yīng)用中，對尺寸較大消聲器，劃分子結(jié)構(gòu)進(jìn)行傳遞損失計算時，混合方法優(yōu)勢較明顯。

表2　消聲器傳遞損失計算時間比較

5結(jié)論

(1) 對含若干等截面、非等截面子結(jié)構(gòu)消聲器，可分別考慮用配點法、邊界元方法計算各自阻抗矩陣，并將其連接進(jìn)行傳遞損失計算。

(2) 壓縮阻抗矩陣可節(jié)省機(jī)時。簡化混合方法在保證計算精度同時能有效減少計算時間，尤其適用具有較長等截面子結(jié)構(gòu)大尺度聲學(xué)問題計算。

參 考 文 獻(xiàn)

[1] Selamet A, Ji Z L. Acoustic attenuation performance of circular expansion chambers with offset inlet/outlet: I.analytical approach[J]. Journal of Sound and Vibration, 1998, 213(4): 601-617.

[2] Selamet A, Ji Z L. Acoustic attenuation performance of circular expansion chambers with extended inlet/outlet[J]. Journal of Sound and Vibration, 1999, 223(2): 197-212.

[3] Selamet A, Ji Z L. Acoustic attenuation performance of circular flow-reversing chambers[J]. Journal of the Acoustical Society of America, 1998, 104(5): 2867-2877.

[4] Selamet A, Xu M B, Lee I J, et al. Analytical approach for sound attenuation in perforated dissipative silencers with inlet/outlet extensions[J]. Journal of the Acoustical Society of America, 2005, 117(4): 2078-2089.

[ 5] Peat K S, Rathi K L. A finite element analysis of the convected acoustic wave motion in dissipative silencers[J]. Journal of Sound and Vibration, 1995, 184(3): 529-545.

[6] 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.

[7] Wu T W, Cheng C Y R, Zhang P. A direct mixed-body boundary element method for packed silencers[J]. Journal ofthe Acoustical Society of America, 2002, 111(6): 2566-2572.

[8] Lou G, Wu T W, Cheng C Y R. Boundary element analysis of packed silencers with a substructuring technique[J]. Engineering Analysis with Boundary Elements, 2003,27(7): 643-653.

[9] Kirby R, Lawrie J B. A point collocation approach to modelling large dissipative silencers[J]. Journal of Sound and Vibration, 2005, 286(1): 313-339.

[10] Kirby R, Williams P T, Hill J. A three dimensional investigation into the acoustic performance of dissipative splitter silencers[J]. The Journal of the Acoustical Society of America, 2014, 135(5): 2727-2737.

[11] Kirby R, Amott K, Williams P T, et al. On the acoustic performance of rectangular splitter silencers in the presence of mean flow[J]. Journal of Sound and Vibration, 2014, 333(24):6295-6311.

[12] Kirby R. A comparison between analytic and numerical methods for modelling automotive dissipative silencers with mean flow[J]. Journal of Sound and Vibration, 2009, 325(3):565-582.

[13] Fang Z, Ji Z L. Acoustic attenuation analysis of expansion chambers with extended inlet/outlet[J]. Noise Control Engineering Journal, 2013, 61(2):240-249.

[14] Fang Z, Ji Z L. Numerical mode matching approach for acoustic attenuation predictions of double-chamber perforated tube dissipative silencers with mean flow[J]. Journal of Computational Acoustics, 2014, 22(2):1450004-1-1450004-15.

[15] Kirkup S M. Fortran codes for computing the discrete Helmholtz integral operators[J]. Advances in Computational Mathematics, 1998, 9(3/4): 391-409.

[16] Krishnasamy G, Schmerr L W, Rudolphi T J, et al. Hypersingular boundary integral equations: some applications in acoustic and elastic wave scattering[J]. Journal of Applied Mechanics, 1990, 57(2): 404-414.

[17] 方智,季振林. 穿孔管阻性消聲器橫向模態(tài)和聲學(xué)特性計算與分析[J]. 振動與沖擊, 2014, 33 (7):138-146.

FANG Zhi, JI Zhen-lin. Transversal modes and acoustic attenuation performance of a perforated tube dissipative silencer[J]. Journal of Vibration and Shock, 2014, 33(7): 138-146.

[18] 方智,季振林. 穿孔管消聲器橫截面模態(tài)及消聲特性的有限元分析[J]. 振動與沖擊, 2012, 31(17): 190-194.

FANG Zhi, JI Zhen-lin. Finite element analysis of cross-section modes and acoustic attenuation characteristics of a perforated tube silencer[J]. Journal of Vibration and Shock, 2012, 31(17):190-194.

基金項目：國家自然科學(xué)基金(11174065)

收稿日期：2014-11-12修改稿收到日期：2015-01-20

通信作者季振林 男，教授，博士生導(dǎo)師，1965年生

中圖分類號:TB535

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.02.025

Transmission loss prediction for silencers by using the combination of boundary element method and numerical collocation approach

YANG Liang1, JI Zhen-lin1, WU T W2

(1. School of Power and Energy, Harbin Engineering University, Harbin 150001, China; 2. College of Mechanical Engineering, University of Kentucky, Lexington 40506, Lexington)

Abstract:A technique combining the boundary element method (BEM) and numerical collocation approach was proposed to calculate the transmission loss of silencers with arbitrary cross-section. The silencer was divided into several substructures. The BEM and the numerical collocation approach were employed combinedly to evaluate the impedance matrices of substructures with non-regular cross-section and axially uniform cross-section respectively, and then the impedance matrices for all substructures were achieved. The transmission loss predictions of several silencers agree well with the results by BEM. The computational efficiency was analyzed and compared with that of the traditional substructure method, the combined technique may save computational time efficiently.

Key words:silencer; transmission loss; boundary element method; numerical collocation approach

第一作者 楊亮 男，博士生，1989年9月生