李 亞，張 楠

(1. 中國船舶科學(xué)研究中心船舶振動噪聲重點實驗室，江蘇無錫214082；2. 中國船舶科學(xué)研究中心水動力學(xué)重點實驗室，江蘇無錫214082)

0 引 言

聲波透過均勻介質(zhì)的現(xiàn)象在生活中比較普遍，對于無限延伸的介質(zhì)，經(jīng)典聲學(xué)教材往往直接假定介質(zhì)中的波是平面波，在界面處速度與壓力連續(xù)，從而可以得到透聲量。實際情況下，介質(zhì)往往是有限的，而且存在邊界上的限制。對于平面波激勵矩形板的輻射聲計算，在文獻[1]中已有較為詳細的分析：平面波激勵矩形板時由于邊界阻抗突變，板內(nèi)波在邊界處反射，因而在板中形成駐波，駐波振幅具有一定形式的分布，利用模態(tài)的正交性，可以得到板在介質(zhì)中的振動方程，然后使用瑞利公式進行輻射聲計算。

在文獻[1]中求矩陣元素時出現(xiàn)了帶奇點的四重積分。對于這一問題，在文獻[2]中采用了變量替換與積分區(qū)域變換的辦法，將四重積分變?yōu)閮芍胤e分，變換過程較為復(fù)雜，未專門分析奇點積分的求解方法。同樣在文獻[3]中，也采用了變量替換與積分區(qū)域變換方法，亦沒有專門分析奇點的處理辦法。對于一重積分可以使用LOBATTO法則計算[4]，以避免被積函數(shù)的奇異性。對于二重帶奇點積分一般是以奇點為圓心作一個小圓，來研究當(dāng)圓半徑逼近0時的結(jié)果[5-6]，這些多用于理論分析。本文采用理論推導(dǎo)，給出了求解奇點積分的具體辦法，然后再用數(shù)值求解。對于矩陣求解，將在正文中看到，里面出現(xiàn)了兩個相加符號，需要巧妙設(shè)計形成矩陣再進行求解。

文獻[1]中還細致介紹了帶周期柵的膜振動及其聲場、周期性支撐固定的薄板振動和聲散射、帶周期性加強肋的薄板振動和聲輻射，這些內(nèi)容由簡到繁地涵蓋了工程中常見的典型情況。應(yīng)該說在2001年左右，平面波激勵矩形板的輻射聲理論研究已經(jīng)達到一個頂峰。因此近年來，對于平面波激勵矩形板的輻射聲理論分析的研究并不多見，主要集中在復(fù)雜板的隔聲問題[7-8]、輻射聲特性[10-11]、輻射效率[12-14]、板的優(yōu)化設(shè)計方面[15]。

另外，采用聲學(xué)有限元、聲學(xué)邊界元也可以研究矩形板的輻射聲。Suzuki等[16]將有限元、邊界元兩種方法耦合，即將聲腔用邊界元離散，將結(jié)構(gòu)用有限元離散，計算了聲波經(jīng)由彈性單層板向聲腔透射的傳聲損失。Wu等[17]用邊界元法(Boundary Element Method, BEM)對薄殼結(jié)構(gòu)內(nèi)外聲場進行了分析，對結(jié)構(gòu)用有限元法(Finite Element Method,FEM)分析，建立復(fù)雜結(jié)構(gòu)的內(nèi)、外聲場和結(jié)構(gòu)耦合振動方程，分析了聲波通過薄壁殼體的聲傳輸問題。Coyette[18]基于有限元和邊界元耦合方法，計算彈性和多孔材料結(jié)構(gòu)的表面阻抗和傳聲損失。Ding等[19]計算了聲源和結(jié)構(gòu)載荷共同激勵下彈性薄壁腔體的聲振耦合特性。對于復(fù)雜情況、單一工況問題有限元或邊界元法能夠很好地解決，但難以發(fā)現(xiàn)機理與影響規(guī)律，在設(shè)計中還要從模態(tài)分析入手，以達到工程要求。

平面波激勵矩形板的輻射聲是很基礎(chǔ)的問題，但由于是一個聲固耦合的問題，因此這個幾何形式簡單的問題并不是看起來那樣簡單，仍有很多機理有待深入研究與挖掘，在聲學(xué)領(lǐng)域中具有獨特的研究價值。這方面研究的進步，也會給其他研究帶來啟發(fā)，因此有很大的研究意義。

1 基本理論公式推導(dǎo)

矩形板由于入射波的作用，會引起板振動，于是向兩邊介質(zhì)中輻射聲波，同時板面也受到兩面介質(zhì)對板面的反作用力，所以這是一個聲固耦合的復(fù)雜問題。

當(dāng)板較大時，可以近似用聲波穿過無限延伸媒質(zhì)上的透聲理論進行計算。本節(jié)首先回顧一下透聲理論，然后重點推導(dǎo)平面波激勵矩形板的輻射聲理論公式。

聲波通過中間層的示意圖如圖1所示，透射波(pt, vt)通過無限大中間層的聲壓之比[20]tp為

圖1 平面波激勵矩形板輻射聲波的示意圖Fig.1 Schematic diagram of soundwaves from a rectangular plate excited by plane waves

平面波激勵矩形板的輻射聲示意圖如圖 2所示，其中板厚h的值較小，屬于薄板。xOz平面位于板的下表面。矩形板x方向長為l1，y方向長為l2。當(dāng)聲波的波向量在沿xOz平面內(nèi)與Oz軸交角為θ方向入射，則在z＜0半空間內(nèi)除了剛硬板面入射波和無限大剛板反射波之外，還應(yīng)加上有限板在聲波作用下進行振動時再輻射的聲波。在z＜h半空間只存在再輻射聲波，于是在板面上下介質(zhì)中的聲關(guān)系如下面第2節(jié)中的公式。

圖2 平面入射波的傳播示意圖Fig.2 Schematic diagram of plane incident wave propagation

其中：μmn是對應(yīng)于本征函數(shù)的本征值，它們決定于板的邊界條件；μ為泊松比；E為楊氏模量；D′為損耗因子；ρ為介質(zhì)密度。令(x1, y1)為矩形板表面上的點坐標(biāo)，(x1′, y1′)是面上積分的動點坐標(biāo)，η為板表面位移幅值。Bmn為有阻尼情況下的振動位移幅值展開項的系數(shù)。

1.1 矩陣求解

式(5)中有兩個求和符號，需用下方矩陣形式進行展開：

1.2 輻射聲求解

2 計算對象與計算方法

2.1 計算模型

計算對象為一薄鋼板，長度為 4 m，寬度為4 m，厚度為0.001 m，放置方向如圖2所示，其中四邊簡支，鋼板四周為無限大障板。鋼的彈性模量216 GPa，密度為7 800 kg·m-3，體縱波聲速為 6 100 m·s-1[19]，泊松比為 0.28?？諝獾膮?shù)：密度為 1.21 kg·m-3，聲速為 344 m·s-1。聲波頻率為 261.6 Hz(音名 c1)，如果頻率太高，則計算階數(shù)增加，太小則考慮到近場效應(yīng)平板面積要增大，這樣積分時區(qū)間劃分個數(shù)又增加，因此定為這一頻率。

2.2 平面波通過矩形板的理論解

將2.1節(jié)中參數(shù)代入式(1)中，可以得到透聲值為0.0648，這個值為無限大板的透聲結(jié)果。由于板較大，可以近似作為平面波通過矩形板的輻射聲。

2.3 平面波激勵矩形板輻射的數(shù)值解

2.3.1 奇點處理

圖3 奇點計算示意圖Fig.3 Singularity computation diagram

圖4 不同ε時解的變化情況Fig.4 Solution variation with different ε

另外，在編程求積分時，采用矩形以便于計算，這樣π變?yōu)?，積分最終表達為

其中：SL為邊長為L的小正方形。

在具體進行積分時，采用圖5的示意圖。首先根據(jù)(x1, y1)的位置，確定好積分區(qū)域。按從上往下的順序分別考慮且這三種情況，在每種情況下又細分為的情況。對于正方形位于矩形中間的情況，可分為Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ五個區(qū)域，其中L為邊長(見圖6)。在計算時，分別計算各區(qū)域的積分值，然后全部加在一起，其中區(qū)域Ⅴ采用方括號中的方法計算。對于處于邊線位置，仍然采用同樣辦法，只不過區(qū)域Ⅰ中積分為 0，區(qū)域Ⅴ中積分只采用有填充部分的面積。

圖5 積分計算示意圖Fig.5 Integral calculation diagram

圖6 積分計算區(qū)域劃分示意圖Fig.6 Partition of integral calculation area

2.3.2 積分處理

在矩陣求解時多個元素均出現(xiàn)多重積分，分別采用計算重積分的累次積分法。對于積分函數(shù)有振蕩的情況，一般可用復(fù)合辛普森法或復(fù)合 COTES方法[22]，由于這里要考慮計算精度與計算時間，經(jīng)比較后每重積分均采用復(fù)化柯特斯公式[23]，公式為

由于一般數(shù)學(xué)軟件自帶公式無法計算四重積分，在分析三重積分時，對于長為1 m、寬為1 m、厚為 0.001 m的板，x1(下面公式中用 x表示)設(shè)為0.62，窄邊劃分區(qū)間數(shù)為1，其他邊為16，進行三重積分計算。用Mathematic軟件進行三維驗證，命令為NIntegrate[Sin[1Pi*z]*Sin[1Pi*w]*E^(-I(2*Pi*261.6)/344*Sqrt[(z-x)^2+(w-y)^2])/Sqrt[(z-x)^2+(w-y)^2]Sin[3Pi*x]*Sin[3Pi*y],{y,0,1},{z,0,1},{w,0,1}]/x-＞0.62，結(jié)果為0.0397 207+0.035 285 5j。程序計算結(jié)果為0.039 610 0+0.035 285 5j，誤差不超過1%。

2.3.3 求逆處理

由于各個元素為復(fù)數(shù)，所以矩陣為復(fù)矩陣，復(fù)矩陣求逆采用全選主元高斯-約當(dāng)(Gauss-Jordan)法[23]。

3 計算結(jié)果討論

圖7 輻射聲的值計算結(jié)果Fig.7 Numerical calculation results of radiation sound

由圖7可以看出，在正方形中間部分的結(jié)果是十分接近理論結(jié)果的，其中心處的值為0.061 462，比理論結(jié)果0.064 8(見2.2節(jié)中)少5.2%；在邊線處，透聲值降低。這說明采用這種辦法求奇點與求輻射聲是可行的。另外，在試算時，在 15～27階變化時，對計算結(jié)果基本影響不大。這是因為頻率較高時，振動模態(tài)為角模態(tài)[24]，即大部分振動區(qū)域的振動輻射聲會相互抵消，角上的振動對輻射聲有較大影響。

值得指出的是，在采用瑞利公式計算活塞聲源輻射聲時，有遠場弗朗霍夫爾(Fraunhofer)區(qū)和近場費涅爾(Fresnel)區(qū)[25]，已經(jīng)比較復(fù)雜。這里有多個振動模態(tài)，情況更為復(fù)雜。在本文編程計算時發(fā)現(xiàn)，當(dāng)聲場位置與平板的距離增大時，中心處與邊線處的輻射聲結(jié)果比較接近，但由于其遠場也類似于活塞聲源，具有遠場特點，其輻射聲值有所降低。

本文基于振動模態(tài)的計算方法有以下近似，這也是計算的誤差來源：薄板厚度忽略不計；理論模型只考慮彎曲波振動輻射聲；奇點計算時，在計算公式將推導(dǎo)中的圓區(qū)域變?yōu)檎叫螀^(qū)域，近似認(rèn)為奇點附近相等；柯特斯法求積分時區(qū)間數(shù)不夠大；計算階次選取還可以更高些；阻尼的設(shè)置對結(jié)果有一些影響；計算采用雙精度，雖已很精確，但仍存在誤差。

4 結(jié) 論

平面波垂直激勵矩形板的輻射聲，貌似十分簡單，但若要精確解決卻需要很繁難的推導(dǎo)與大計算量的數(shù)值計算。

本文細致地推導(dǎo)了平面波激勵矩形板的輻射聲公式，重點處理了矩陣求解問題、四重帶奇點積分，并進行了程序編寫，其中四重帶奇點積分能夠劃分成五個區(qū)域分別求解來最終獲得。針對某一具體算例，求得的結(jié)果與無限大板的透射聲理論值十分接近。這說明采用模態(tài)正交分析是可行的。另外，采用這種辦法會準(zhǔn)確得到各個位置處的輻射聲。這種辦法再加以拓展還可以處理曲面輻射聲的問題。