平面波激勵(lì)矩形板輻射聲計(jì)算中奇異積分處理方法研究

李 亞，張 楠

(1. 中國船舶科學(xué)研究中心船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，江蘇無錫214082；2. 中國船舶科學(xué)研究中心水動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，江蘇無錫214082)

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

1.1 矩陣求解

式(5)中有兩個(gè)求和符號(hào)，需用下方矩陣形式進(jìn)行展開：

1.2 輻射聲求解

2 計(jì)算對(duì)象與計(jì)算方法

2.1 計(jì)算模型

計(jì)算對(duì)象為一薄鋼板，長度為 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)，如果頻率太高，則計(jì)算階數(shù)增加，太小則考慮到近場(chǎng)效應(yīng)平板面積要增大，這樣積分時(shí)區(qū)間劃分個(gè)數(shù)又增加，因此定為這一頻率。

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

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

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

2.3.1 奇點(diǎn)處理

圖3 奇點(diǎn)計(jì)算示意圖Fig.3 Singularity computation diagram

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

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

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

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

圖5 積分計(jì)算示意圖Fig.5 Integral calculation diagram

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

2.3.2 積分處理

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

由于一般數(shù)學(xué)軟件自帶公式無法計(jì)算四重積分，在分析三重積分時(shí)，對(duì)于長為1 m、寬為1 m、厚為 0.001 m的板，x1(下面公式中用 x表示)設(shè)為0.62，窄邊劃分區(qū)間數(shù)為1，其他邊為16，進(jìn)行三重積分計(jì)算。用Mathematic軟件進(jìn)行三維驗(yàn)證，命令為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。程序計(jì)算結(jié)果為0.039 610 0+0.035 285 5j，誤差不超過1%。

2.3.3 求逆處理

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

3 計(jì)算結(jié)果討論

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

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

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

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

4 結(jié) 論

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

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