楊　亮， 季振林， Wu T W

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

?

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

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

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

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

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

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

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

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

1直接混體邊界元方法

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

(1)

(2)

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

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

(3)

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

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

(4)

2數(shù)值配點(diǎn)法

圖1的子結(jié)構(gòu)B沿軸向等截面，但由于截面形狀任意，無(wú)法用解析法獲得橫向波數(shù)，故用數(shù)值方法進(jìn)行特征值分析，考慮二維有限元及三維邊界元網(wǎng)格劃分的一致性，為便于各子結(jié)構(gòu)阻抗矩陣連接，用二維有限元計(jì)算特征值及特征向量。以阻性穿孔管道為例，圖1中Ⅰ的內(nèi)部為空氣，Ⅱ?yàn)樘畛湮暡牧?，Ⅰ與Ⅱ之間由穿孔管分隔。區(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)。對(duì)穿孔管抗性消聲器，式(10)簡(jiǎn)化[18]為

(13)

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

(14)

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

(15)

(16)

(17)

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

(18)

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

(19)

3消聲器傳遞損失計(jì)算

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

(20)

式中：

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

(21)

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

(22)

傳遞損失計(jì)算公式可表示為

(23)

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

4算例驗(yàn)證

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

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

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

模態(tài)解析解數(shù)值解相對(duì)誤差/%(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截面會(huì)劃分較多網(wǎng)格，對(duì)式(18)、(19)，求解模態(tài)幅值系數(shù)會(huì)費(fèi)較多計(jì)算時(shí)間，實(shí)際應(yīng)用中，考慮過(guò)多模態(tài)無(wú)意義。本算例計(jì)算頻率最高到2 500 Hz，(3,0)模態(tài)或(0,3)模態(tài)激發(fā)頻率約3 440 Hz，實(shí)際計(jì)算中考慮最高到第二階模態(tài)即可滿足需要。

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

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

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

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

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

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

表2　消聲器傳遞損失計(jì)算時(shí)間比較

5結(jié)論

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

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

參 考 文 獻(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ì)算與分析[J]. 振動(dòng)與沖擊, 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]. 振動(dòng)與沖擊, 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.

基金項(xiàng)目：國(guó)家自然科學(xué)基金(11174065)

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

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

中圖分類號(hào):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月生