李家柱 葉望青 鐘秤平 陳 劍 劉 策

(1 合肥工業(yè)大學(xué)機(jī)械工程學(xué)院噪聲振動(dòng)工程研究所 合肥 230009)

(2 安徽省汽車(chē)NVH 工程技術(shù)研究中心 合肥 230009)

(3 江西省汽車(chē)噪聲與振動(dòng)重點(diǎn)實(shí)驗(yàn)室 南昌 330001)

0 引言

孔口模態(tài)輻射阻抗的計(jì)算是矩形法蘭孔聲傳遞損失(Transmission loss, TL)計(jì)算的基礎(chǔ)，是開(kāi)展大尺寸百葉窗隔聲量計(jì)算的關(guān)鍵環(huán)節(jié)。大尺寸百葉窗結(jié)構(gòu)在建筑物外立面、暖通系統(tǒng)以及高鐵減載式聲屏障等領(lǐng)域應(yīng)用廣泛。矩形法蘭孔孔口模態(tài)輻射阻抗表達(dá)式是包含半空間格林函數(shù)和兩組模態(tài)振型的四重積分，該積分存在奇異點(diǎn)，現(xiàn)有的高階積分算法和計(jì)算軟件無(wú)法直接求解[1]。對(duì)于這類(lèi)問(wèn)題，多數(shù)學(xué)者通過(guò)坐標(biāo)變換將積分降為三重或二重積分，然后采用高斯求積來(lái)求解[2?5]。這種方法實(shí)現(xiàn)了積分的精確求解，但在積分變換后，表達(dá)式需根據(jù)相應(yīng)的階次進(jìn)行討論，且高斯求積方法雖較成熟，但并非MATLAB 等軟件的內(nèi)置函數(shù)，仍需編寫(xiě)高斯求積程序。Sha 等[1]則在將四重積分降為二重積分后，將結(jié)構(gòu)劃分為網(wǎng)格，利用邊界條件、波數(shù)近似和傅里葉變換等方法來(lái)加速求解模態(tài)輻射阻抗，該方法需根據(jù)結(jié)構(gòu)的尺寸、波數(shù)等參數(shù)改變網(wǎng)格大小，計(jì)算過(guò)程較復(fù)雜且計(jì)算速度不夠理想。Li 等[6?7]、沈蘇等[8]和Leppington 等[9]等通過(guò)坐標(biāo)變換、參數(shù)討論等方法將這類(lèi)積分降為一重積分，簡(jiǎn)化了積分最終的求解表達(dá)式并提高了計(jì)算速度。Davy等[10]在將四重積分降為一重積分的基礎(chǔ)上，通過(guò)改寫(xiě)最終表達(dá)式消除了奇異積分的影響。然而，這類(lèi)將四重積分降為一重積分的方法推導(dǎo)過(guò)程復(fù)雜且需進(jìn)行參數(shù)討論，如果初始的模態(tài)輻射阻抗表達(dá)式與文獻(xiàn)中的表達(dá)式有差異，就需重新推導(dǎo)，對(duì)理論知識(shí)要求較高且工作量大。

鑒于此，本文在吸收以上學(xué)者成果的基礎(chǔ)上，通過(guò)坐標(biāo)變換將四重積分轉(zhuǎn)換為二重積分并消除奇異積分，且將該過(guò)程直接采用MATLAB 的內(nèi)置函數(shù)來(lái)實(shí)現(xiàn)。在實(shí)現(xiàn)四重積分求解的基礎(chǔ)上，通過(guò)詳細(xì)計(jì)算和對(duì)比，驗(yàn)證方法的正確性，最后分析矩形法蘭孔孔口模態(tài)輻射阻抗的性質(zhì)。

1 矩形法蘭孔孔口模態(tài)輻射阻抗計(jì)算

1.1 矩形法蘭孔孔口模態(tài)輻射阻抗表達(dá)式

圖1 矩形法蘭孔示意圖Fig.1 Schematic diagram of flanged rectangular aperture

聲源表面的聲壓與振速的比值稱為聲源的輻射阻抗[11]。圖1為矩形法蘭孔示意圖，該孔孔口處為無(wú)限大壁面形成的半空間，孔口的振動(dòng)單元可視為空氣薄膜，由于受孔口約束，空氣薄膜的四邊僅有圖1所示的z向平動(dòng)自由度，其振動(dòng)模態(tài)對(duì)應(yīng)孔口截面模態(tài)?？諝獗∧じ麟A振動(dòng)模態(tài)對(duì)應(yīng)的輻射阻抗即為矩形法蘭孔孔口模態(tài)輻射阻抗，其表達(dá)式為[5]

式(1)中，j為虛數(shù)單位，kf和Zf分別為波數(shù)和聲傳播介質(zhì)的特征阻抗，本文以空氣為聲傳播介質(zhì)，則Zf為空氣的特征阻抗，?mn(x,y)和?pq(x0,y0) 分別為空氣薄膜的(m,n)和(p,q)階模態(tài)的振型，表達(dá)式分別為

G(x,y,x0,y0)為半空間格林函數(shù)，其表達(dá)式為

將 式(2)～(4)以 及S= 4ab代 入 式(1)可 得(m,n,p,q)階模態(tài)輻射阻抗表達(dá)式為

1.2 積分變換

式(5)是以x0、y0、x、y為積分變量的四重積分，積分區(qū)間分別為(?a,a)、(?b,b)、(?a,a)、(?b,b)，在x=x0且y=y0處，積分函數(shù)趨于無(wú)窮大，無(wú)法直接計(jì)算。為求解該積分，采用換元法，令

可得

雅可比矩陣為

同理，Jy=1/2，此時(shí)式(5)變換為

式(9)中：

式(10)～(13)均可對(duì)ζ和γ積出解析表達(dá)式，式(9)將被轉(zhuǎn)換為四個(gè)以κ和τ為積分變量的二重積分之和。然而，當(dāng)κ和τ趨于零時(shí)依然會(huì)造成積分中格林函數(shù)分母趨于零，該積分存在奇異性，需做一次極坐標(biāo)變換，令：

變換后的式(9)被轉(zhuǎn)換為極坐標(biāo)下的二重積分，該積分可直接使用MATLAB內(nèi)置的數(shù)值積分函數(shù)求解，無(wú)需采用高斯求積或?qū)δB(tài)階次進(jìn)行討論。

1.3 MATLAB求解

第1.2 節(jié)所述過(guò)程可在MATLAB中完成，主要有以下幾步：

(1) 根據(jù)式(2)～(4)，在MATLAB 中定義模態(tài)振型函數(shù)?mn(x,y)、?pq(x0,y0)以及半空間格林函數(shù)G(x,y,x0,y0)。

(2) 將變量代換式(7)帶入第(1)步定義的函數(shù)，得到變量代換后的模態(tài)振型表達(dá)式和半空間格林函數(shù)表達(dá)式。

(3) 將變量代換后的模態(tài)振型和半空間格林函數(shù)表達(dá)式代入式(5)，即得到新的模態(tài)輻射阻抗表達(dá)式(9)，要對(duì)式(9)求解，需分別對(duì)式(10)～(13)求定積分。

(4) 在MATLAB 中 采 用int 函 數(shù) 對(duì) 式(10)～(13)中的變量ζ和γ求定積分，將式(10)～(13)由四重積分轉(zhuǎn)換為僅含變量κ和τ的二重積分，此步在MATLAB中可積出解析表達(dá)式。

(5) 降維后的式(10)～(13)由于包含半空間格林函數(shù)而存在奇異性，此時(shí)在MATLAB 中根據(jù)式(14)進(jìn)行極坐標(biāo)變換，消除奇異積分，得到極坐標(biāo)下的積分表達(dá)式。

(6) 最后利用MATLAB積分函數(shù)quad2d分別對(duì)降維和極坐標(biāo)變換后的式(10)～(13)求數(shù)值積分，將式(10)～(13)求解結(jié)果代入式(9)即可求得矩形法蘭孔孔口模態(tài)輻射阻抗。

2 計(jì)算結(jié)果驗(yàn)證

2.1 (0,0,0,0)階模態(tài)輻射阻抗驗(yàn)證

矩形法蘭孔孔口的(0,0,0,0)階模態(tài)輻射阻抗就是矩形活塞聲源的輻射阻抗，矩形活塞聲源輻射阻抗表達(dá)式和相關(guān)計(jì)算數(shù)據(jù)可在文獻(xiàn)[12]中查閱。如圖2所示，法蘭孔孔口尺寸為長(zhǎng)寬均為1 m 的正方形，圖2中Z(0,0,0,0)為(0,0,0,0)階模態(tài)輻射阻抗，Z0為空氣的特征阻抗，Z(0,0,0,0)/Z0為歸一化的模態(tài)輻射阻抗，k為波數(shù)，req為等效半徑由圖2可知，本方法算得的(0,0,0,0)階模態(tài)輻射阻抗與矩形活塞聲源輻射阻抗實(shí)部和虛部一致性很好，說(shuō)明本方法計(jì)算(0,0,0,0)階模態(tài)輻射阻抗是正確的。

圖2 (0,0,0,0)階模態(tài)輻射阻抗驗(yàn)證Fig.2 Radiation impedance validation of order(0,0,0,0)

2.2 高階模態(tài)輻射阻抗驗(yàn)證

由于高階模態(tài)輻射阻抗難以通過(guò)仿真或?qū)嶒?yàn)直接得到，因此本文將模態(tài)輻射阻抗計(jì)算方法代入矩形法蘭孔聲傳遞損失計(jì)算中來(lái)。通過(guò)代入本文阻抗計(jì)算方法的Superposition 法與Sgard 等[5]、聲學(xué)有限元法以及實(shí)驗(yàn)[13]的對(duì)比，間接驗(yàn)證本方法的正確性。

矩形法蘭孔的聲傳遞率與其孔口模態(tài)輻射阻抗的關(guān)系為[13]

式(15)中，Zmnmn為孔口第(m,n,m,n)階模態(tài)輻射阻抗，該式已利用本文的結(jié)論之一：互模態(tài)輻射阻抗相比于自模態(tài)輻射阻抗在多數(shù)計(jì)算中可以忽略不計(jì)的性質(zhì)進(jìn)行了簡(jiǎn)化，只計(jì)算了自模態(tài)輻射阻抗部分。關(guān)于孔口模態(tài)輻射阻抗性質(zhì)的分析將在本文后續(xù)小節(jié)中加以論述。

聲傳遞損失與聲傳遞率的關(guān)系為[13]從而可以根據(jù)式(16)求得矩形法蘭孔的聲傳遞損失，式(15)和式(16)中部分參數(shù)的具體描述請(qǐng)見(jiàn)參考文獻(xiàn)[13]。

圖3和圖4對(duì)比了代入本文阻抗計(jì)算方法的Superposition法[13]和Sgard等[5]方法的計(jì)算結(jié)果，矩形孔尺寸均為L(zhǎng)x=0.4 m，Ly=0.2 m，Lz=0.3 m。圖3為法向入射時(shí)的聲傳遞損失，圖4是入射角度為斜45?時(shí)的聲傳遞損失。圖3中Sgard 等的數(shù)據(jù)為采用其論文中的算法計(jì)算得到，經(jīng)對(duì)比二者幾乎相等。圖4中Sgard 等的數(shù)據(jù)為直接從其論文的圖中提取出來(lái)，略微存在一定誤差，這主要是提取數(shù)據(jù)的誤差造成的。

圖3 法向入射條件下矩形孔聲傳遞損失驗(yàn)證Fig.3 TL validation of a rectangular aperture with normal incident

圖4 傾斜入射條件下矩形孔聲傳遞損失驗(yàn)證Fig.4 TL validation of a rectangular aperture with oblique incident

圖5為散射聲場(chǎng)中代入本文阻抗計(jì)算方法的Superposition 法[13]、Trompette 等[14]的實(shí)驗(yàn)以及聲學(xué)有限元法的對(duì)比，矩形孔尺寸為L(zhǎng)x=0.06 m，Ly=0.13 m，Lz=0.3 m。觀察圖5可知，三條曲線一致性較好。

以上通過(guò)將本文模態(tài)輻射阻抗計(jì)算方法代入Superposition 法， 與Sgard 等的計(jì)算結(jié)果、Trompette 等的實(shí)驗(yàn)結(jié)果以及聲學(xué)有限元法的結(jié)果對(duì)比，取得了很好的一致性，從側(cè)面驗(yàn)證了本方法計(jì)算高階模態(tài)輻射阻抗的正確性。

圖5 散射聲場(chǎng)中矩形孔聲傳遞損失對(duì)比Fig.5 TL validation of a rectangular aperture in diffuse acoustic field

3 計(jì)算結(jié)果分析

3.1 自模態(tài)輻射阻抗分析

當(dāng)式(1)中m=p或n=q時(shí)的模態(tài)輻射阻抗稱為自模態(tài)輻射阻抗，圖6和圖7為孔口邊長(zhǎng)為1 m的正方形法蘭孔孔口的前7 階自模態(tài)輻射阻抗實(shí)部和虛部的對(duì)比。

觀察圖6和圖7，可得如下結(jié)論：

(1)隨著模態(tài)階次的增加，自模態(tài)輻射阻和自模態(tài)輻射抗的峰值總體呈減小趨勢(shì)，但不是單調(diào)遞減的，高階次的峰值可能高于低階次的峰值。

(2)隨著模態(tài)階次的增加，自模態(tài)輻射阻和自模態(tài)輻射抗的峰值對(duì)應(yīng)的kreq均逐漸增大，由于req為等效半徑，為常數(shù)，即對(duì)應(yīng)的波數(shù)k增大，亦即對(duì)應(yīng)的頻率逐漸增大。

圖6 自模態(tài)輻射阻對(duì)比(Lx = 1 m、Ly = 1 m)Fig.6 Resistance comparison of self-modal radiation impedance(Lx = 1 m、Ly = 1 m)

圖7 自模態(tài)輻射抗對(duì)比(Lx =1 m、Ly =1 m)Fig.7 Reactance comparison of self-modal radiation impedance(Lx =1 m、Ly =1 m)

圖8 自模態(tài)輻射阻對(duì)比(Lx = 0.5 m、Ly = 1 m)Fig.8 Resistance comparison of self-modal radiation impedance(Lx = 0.5 m、Ly = 1 m)

圖9 自模態(tài)輻射抗對(duì)比(Lx = 0.5 m、Ly = 1 m)Fig.9 Reactance comparison of self-modal radiation impedance(Lx = 0.5 m、Ly = 1 m)

(3)孔口形狀為正方形時(shí)，由于對(duì)稱性，(m,n,m,n)與(n,m,n,m) 階自模態(tài)輻射阻和自模態(tài)輻射抗相等。

進(jìn)一步，選取孔口邊長(zhǎng)分別為L(zhǎng)x= 0.5 m、Ly=1 m的矩形法蘭孔，計(jì)算其孔口模態(tài)輻射阻抗并分析其規(guī)律。

觀察圖8和圖9，并結(jié)合圖6和圖7可知，當(dāng)法蘭孔孔口形狀為長(zhǎng)方形時(shí)，其模態(tài)輻射阻抗的性質(zhì)與孔口形狀為正方形時(shí)相似。不同之處在于，其不具有(m,n,m,n)與(n,m,n,m)階模態(tài)輻射阻和模態(tài)輻射抗相等的對(duì)稱性，是因?yàn)殚L(zhǎng)方形的邊長(zhǎng)不相等。

3.2 互模態(tài)輻射阻抗分析

(1)四個(gè)階次的互模態(tài)輻射阻和輻射抗的峰值相比于自模態(tài)輻射阻和輻射抗的峰值均小很多，峰值最大的(1,2,3,4)階互模態(tài)輻射阻抗也比自模態(tài)輻射阻抗小1～2個(gè)數(shù)量級(jí)。

(2)由圖10、 圖11 和表1可知， 另三階比(1,2,3,4)階模態(tài)輻射阻抗小多個(gè)數(shù)量級(jí)。

圖10 互模態(tài)輻射阻對(duì)比(Lx = 0.5 m、Ly = 1 m)Fig.10 Resistance comparison of cross-modal radiation impedance(Lx = 0.5 m、Ly = 1 m)

圖11 互模態(tài)輻射抗對(duì)比(Lx = 0.5 m、Ly = 1 m)Fig.11 Reactance comparison of cross-modal radiation impedance(Lx = 0.5 m、Ly = 1 m)

表1 互模態(tài)輻射阻抗對(duì)比Table1 Comparison of cross-modal radiation impedance

4 結(jié)論

論文詳細(xì)闡述了矩形法蘭孔孔口模態(tài)輻射阻抗積分表達(dá)式坐標(biāo)變換的過(guò)程及其在MATLAB中的實(shí)現(xiàn)方法，對(duì)計(jì)算結(jié)果進(jìn)行了驗(yàn)證并分析了模態(tài)輻射阻抗的性質(zhì)，結(jié)論如下：

(1)通過(guò)坐標(biāo)變換將四重積分轉(zhuǎn)換為二重積分并消除了奇異積分，且將計(jì)算過(guò)程直接利用MATLAB 軟件的內(nèi)置函數(shù)實(shí)現(xiàn)，不再需要推導(dǎo)公式且不再需要根據(jù)模態(tài)階次討論模態(tài)輻射阻抗的具體表達(dá)式，顯著降低了積分求解的復(fù)雜程度。

(2)隨著模態(tài)階次的增加，自模態(tài)輻射阻抗峰值總體呈減小趨勢(shì)但并非單調(diào)遞減，峰值對(duì)應(yīng)的頻率逐漸增大。

(3)自模態(tài)輻射阻抗顯著大于互模態(tài)輻射阻抗，往往相差多個(gè)數(shù)量級(jí)，結(jié)合實(shí)際情況可只計(jì)算自模態(tài)輻射阻抗部分，以顯著降低計(jì)算工作量、提高計(jì)算速度。