曹樹杰，宋春旺，米 崢

(1.中海油田服務(wù)股份有限公司， 河北 三河 065201； 2.大連理工大學(xué) 船舶工程學(xué)院， 遼寧 大連 116000)

物探船作為勘探資源領(lǐng)域的關(guān)鍵裝備，在國家戰(zhàn)略和經(jīng)濟發(fā)展問題上的研究熱度經(jīng)久不衰([1]。其工作原理是通過氣槍發(fā)射聲波且分析接收到的回波實現(xiàn)資源探測，故對船體結(jié)構(gòu)設(shè)備的聲輻射特性很敏感[2]，若輻射噪聲過大，容易降低信號采集系統(tǒng)的探測功能嚴重影響物探作業(yè)、增加勘探成本，同時作為營運性船舶其受到的外力激勵一般包括機械設(shè)備激勵、螺旋槳激勵、流體激勵，其中機械設(shè)備的周期性激勵產(chǎn)生的結(jié)構(gòu)振動聲輻射對于作業(yè)環(huán)境的安全和舒適有很大影響，故分析不同工況激勵模式下結(jié)構(gòu)輻射聲場的熱點區(qū)域以實現(xiàn)聲輻射的預(yù)報和控制有極其重要的意義[3]。數(shù)值方法求解聲學(xué)問題主要包括聲學(xué)有限元和聲學(xué)邊界元兩大類[4]。有限元法通常用來計算封閉空間中聲場，對于無限大外空間的聲場計算通常采用聲學(xué)有限元結(jié)合AML(Automatic Matched Layer)方法和聲場邊界元法。其中，邊界元法采取表面網(wǎng)格劃分、求解域內(nèi)解析解疊加的思想，能夠高效、精確地計算結(jié)構(gòu)在聲媒質(zhì)中的聲輻射貢獻[5]。本研究通過自編程序?qū)⑦吔缭椒☉?yīng)用于振動結(jié)構(gòu)的輻射聲場計算，首先驗證了邊界元程序的有效性；其次聯(lián)合結(jié)構(gòu)有限元分析方法分析不同激勵模式下的球殼結(jié)構(gòu)聲場輻射特性，為復(fù)雜結(jié)構(gòu)的聲輻射熱點分析奠定基礎(chǔ)[6]，并通過改變激勵力的大小與方向最終實現(xiàn)輻射聲壓的減小以滿足工程要求。

1 結(jié)構(gòu)輻射聲場邊界元計算方法

1.1 理論分析

具有封閉表面的結(jié)構(gòu)穩(wěn)態(tài)振動外場聲輻射滿足聲學(xué)Helmholtz方程：

▽2p+k2p=0

其中：p為聲壓；k=ω/c為波數(shù)；ω為角頻率；c為介質(zhì)聲速。

該方程基本解為g=e-ikr/4πr，

考慮邊界條件，采用Helmholtz微分方程的基本解，可得Helmholtz積分方程：

其中：R為結(jié)構(gòu)表面S上任意點；P為空間內(nèi)任意點；

其中vn(R)為結(jié)構(gòu)表面振速。

根據(jù)Helmholtz積分方程，若給定結(jié)構(gòu)表面振速作為邊界條件，將式中的P點限定成為結(jié)構(gòu)表面任意點，選擇合適的C(P)，即可求解出結(jié)構(gòu)表面的聲壓值，從而求解域內(nèi)任意一點聲壓。

1.2 數(shù)值分析

1) 離散表面

本文選用四邊形線性等參單元，形函數(shù)為Nl(l=1,2,3,4)，則其單元內(nèi)任意位置坐標(biāo)和物理量均可由頂點物理量值與相應(yīng)形函數(shù)乘積之和表示[7]，設(shè)單元局部坐標(biāo)為(ξ,η),即：

2) 積分方程離散分析

將Helmholtz積分方程中P點限定成為結(jié)構(gòu)表面任意點，其中：

?sjg*(R,Pi)p(R)dS=

式中：m為結(jié)構(gòu)表面節(jié)點數(shù)量；n為結(jié)構(gòu)表面單元數(shù)量；J(ξ,η)為坐標(biāo)系變換的雅可比式；k1、k2分別為ξ、η方向的高斯點；wk1、wk2為積分權(quán)函數(shù)。設(shè)

最終得到Hsps=Gsvns，即若給定結(jié)構(gòu)表面振動速度，便可求解相應(yīng)的結(jié)構(gòu)表面聲壓，從而求解域內(nèi)任意一點聲壓，其離散分析過程與上述類似。

3) 奇異積分問題

本文采用單元子分法解決[8]，如圖1所示，通過對四邊形線性元進行單元分割和坐標(biāo)變換。將單元內(nèi)P點分別與正方形的4個頂點相連組成4個三角形子單元，當(dāng)P點位于正方形頂點時，其對應(yīng)子單元僅剩兩個。后將每個子分的三角形單元映射到另一個坐標(biāo)系(ξ′,η′)成為具有兩個P點和剩余頂點構(gòu)成的新的四邊形單元。新四邊形單元相鄰兩個頂點具有相同的原坐標(biāo)系內(nèi)源點P的坐標(biāo)，使得分割單元后坐標(biāo)變換產(chǎn)生的雅可比式J(ξ′,η′)為一階無窮小量，從而消除弱奇異積分的影響。

4) 特征頻率下解不唯一問題

求解表面Helmholtz積分方程時，特征頻率下的解不能唯一確定，本文采用CHIEF(combined Helmholtz integral equation formulation)方法[9]，利用結(jié)構(gòu)內(nèi)部的點滿足的Helmholtz積分方程作為附加約束條件，得到超定方程組，求解其最小二乘解通過求解作為滿足計算精度的數(shù)值解。假設(shè)給定結(jié)構(gòu)表面振動速度，需求解表面聲壓，設(shè)表面Helmholtz積分方程離散后求解方程組為Hsps=Gsvns，取K個CHIEF點得到其對應(yīng)的Helmholtz積分方程Hcps=Gcvns，其中Hs和Gs均為方陣(行列數(shù)相同)，其維數(shù)與結(jié)構(gòu)表面節(jié)點數(shù)目相同，Hc和Gc的行數(shù)對應(yīng)K個CHIEF點，列數(shù)對應(yīng)結(jié)構(gòu)表面節(jié)點數(shù)目。

2 數(shù)值分析算例

2.1 脈動球驗證程序有效性

本節(jié)采用脈動球的數(shù)值計算解與已知的解析解進行對比，驗證邊界元程序的有效性。

半徑為r0、球面法向振速為v的脈動球殼產(chǎn)生的輻射聲場在距離球心r處場點的聲壓p(r)的解析解[10]為：

其中：ρ為聲學(xué)介質(zhì)的密度；c為聲速；k=/c為波數(shù)。

本文分別采用96個單元和600個單元劃分脈動球結(jié)構(gòu)表面，劃分網(wǎng)格如圖2所示，均采用四邊形線性單元，計算中，取球心為CHIEF點，ρ、c分別取空氣介質(zhì)的密度和聲速，即1.293 kg/m3、340 m/s，r0、v分別取1 m和0.01 m/s。

表1給出了不同波數(shù)下脈動球解析解、96個單元劃分數(shù)值解、600個單元劃分數(shù)值解的表面聲壓對比。從中可以看出本研究采用的聲學(xué)邊界元程序計算有效，特征頻率處解不唯一問題得到解決，圖3給出了上述參數(shù)脈動球以210 Hz的頻率振動時，表面聲壓數(shù)值解與解析解誤差與網(wǎng)格劃分數(shù)量的關(guān)系曲線，說明精細網(wǎng)格有利于計算精度的提高。

為解決特征頻率下解不唯一的問題，采用CHIEF方法，將結(jié)構(gòu)內(nèi)部點滿足的Helmholtz積分方程作為附加約束條件組成超定方程組從而控制所求解聲壓的精度。表2給出不同波數(shù)條件下脈動球表面聲壓的解析解、不采用CHIEF方法的數(shù)值解、采用CHIEF方法的數(shù)值解結(jié)果對比，可以看出：在特征頻率kr0=nπ(n=1,2,3,…)時，采用CHIEF方法可以有效提高計算精度，解決特征頻率下解不唯一問題；隨著波數(shù)增大，數(shù)值方法計算精度相應(yīng)會降低，但采用CHIEF方法后，由于加入內(nèi)部結(jié)構(gòu)點的Helmholtz積分方程約束條件，其精度在高波數(shù)情況下也能保證在一定的允許范圍。

表1 不同波數(shù)下脈動球解析解、數(shù)值解表面聲壓

波數(shù)解析解聲壓/Pa96個單元劃分數(shù)值解聲壓/Pa無chief方法96個單元劃分數(shù)值解聲壓/Pa有chief方法12.198 1+2.198 1i2.187 1+2.183 1i2.187 0+2.183 7i23.517 0+1.758 5i3.611 2+1.728 8i3.607 4+1.732 3i33.956 6+1.318 9i5.156 9+1.307 9i4.189 7+1.345 4iπ3.991 8+1.270 6i-0.108 6-28.471 0i4.023 9+1.294 2i44.137 6+1.034 4i3.978 5+0.856 7i4.022 5+0.902 9i54.227 1+0.845 4i4.305 6+0.585 1i4.297 6+0.647 5i64.277 4+0.712 9i5.315 7-0.073 2i4.499 5+0.593 3i2π4.287 6+0.682 4i-1.480 1-13.774i4.327 1+0.684 4i74.308 3+0.615 5i3.720 7+0.506 4i4.0619+0.5772i84.328 6+0.541 1i4.131 8+0.190 3i4.251 0+0.344 4i94.342 6+0.482 5i4.506 5-0.968 7i4.424 3+0.216 2i3π4.347 3+0.461 3i-5.358 0-6.026 0i4.382 6+0.440 8i104.352 7+0.435 3i3.363 6+1.240 1i4.195 1+0.628 9i

圖4給出的是上述參數(shù)脈動球以210 Hz的頻率振動時，不同場點處聲壓幅值的600個單元劃分數(shù)值解與解析解誤差與其距球心距離的關(guān)系曲線，可以看出，距球心大于1.25 m的場點的計算結(jié)果誤差已足夠小，說明本程序可有效計算結(jié)構(gòu)振動輻射聲場；對于距離脈動球表面很近的場點，其結(jié)果存在奇異性，需采取其他方法消除。

圖5給出的是上述參數(shù)脈動球的遠場(r=10 m)輻射聲壓幅值(解析解、600個單元劃分數(shù)值解)隨波數(shù)變化曲線，可以看出，遠場聲壓數(shù)值計算結(jié)果與解析解大致相同，證明此程序能夠有效計算輻射聲場。

2.2 不同激勵模式的球殼輻射聲場分析

對于物探船甲板上一球型設(shè)備產(chǎn)生的振動聲輻射，本文通過結(jié)構(gòu)有限元/聲學(xué)邊界元方法實現(xiàn)不同激勵狀況下的設(shè)備簡諧振動聲輻射特性分析。設(shè)備球心處有一激勵源，設(shè)備底部與甲板連接，本文簡化為剛固條件展開分析。

2.2.1對球殼進行3種激勵方式的輻射聲場分析

結(jié)構(gòu)在不同激勵模式下振動特性不同，其向外輻射聲場差異明顯。本文采用結(jié)構(gòu)有限元方法計算球殼結(jié)構(gòu)在不同激勵模式下的響應(yīng)振動信息，輸入聲場邊界元程序求解其輻射聲場并分析其特性。

所選結(jié)構(gòu)為厚度為0.001 m的鋼制球殼，材料為45號鋼，鋼密度7 850 kg/m3、彈性模量210 GPa、泊松比0.3，空氣聲速340 m/s、密度1.293 kg/m3。球心位于(0,0,0)，約束簡化為點(0,0,-1)全約束，即球殼最下方與甲板接觸點剛固。如圖6(a)～圖6(c)所示，設(shè)中心激勵設(shè)備對球殼激勵模式有三類：一點激勵、兩點激勵、三點激勵，激勵位置均分布于y=0平面上，激勵力均為210 Hz簡諧力。激勵力的頻率避開了球殼結(jié)構(gòu)模態(tài)分析中得到的固有頻率，故阻尼對計算影響較小，可忽略。

通過有限元/邊界元方法分別計算球殼結(jié)構(gòu)受到3種不同激勵模式的輻射聲場中z=0平面距球心10 m處的各點聲壓幅值大小，得到其輻射聲場指向性圖，如圖7(a)～圖7(c)。

圖6 三類激勵模式

根據(jù)上述聲壓圖可以看出，不同激勵模式下的結(jié)構(gòu)振動輻射聲場存在顯著差異。位于球心的設(shè)備對球殼產(chǎn)生激勵時，若球殼結(jié)構(gòu)只受單點簡諧激勵，其遠場輻射聲場分布呈現(xiàn)激勵方向較大，在90°～120°與240°～270°范圍內(nèi)其輻射聲壓較??；受到相對兩點簡諧激勵時，其遠場輻射聲場大小分布呈現(xiàn)激勵方向較為均勻、垂直激勵方向相對較??；而受到如圖6(c)所示的三點激勵時，其遠場聲輻射大小分布呈現(xiàn)軸向較大且x軸方向幅值大于y軸方向，但整體幅值較低。

如上述分析，不同激勵模式的輻射聲場的特性迥異，實際作業(yè)中，需在設(shè)計階段根據(jù)周圍環(huán)境和儀器設(shè)備的安裝條件選擇合適的激勵模式使其負面影響降到最低。當(dāng)球心設(shè)備激勵球殼結(jié)構(gòu)時，單點激勵模式會導(dǎo)致較大的單向聲壓指向性，故該模式適用于在特定方向隔聲效果較好的作業(yè)環(huán)境，其他方向受到影響較??；兩點激勵模式相較來說有較為模糊的指向性，且整體來說聲壓強度明顯減小。若作業(yè)環(huán)境中要對設(shè)備周圍所受聲輻射影響敏感度相當(dāng)，則選擇此激勵方式效果較佳。三點激勵模式由于有一部分能量上下激勵，故其環(huán)向輻射的聲壓數(shù)值最小，其最大聲壓級低于86 dB，選擇此方法可有效降低環(huán)向聲輻射的負面影響。

2.2.2設(shè)置減振設(shè)備的三點激勵結(jié)構(gòu)輻射聲場分析

對于三點激勵模式，若加設(shè)減振裝置于某一激勵位置，改變該方向激勵的大小，其環(huán)向聲輻射特性會有顯著的變化。三點激勵方向仍如圖6(c)所示，若F2和F3幅值降為之前的四分之一，即83.325 N，F(xiàn)1仍為333.3 N時，其環(huán)向聲輻射指向性如圖8(a)；若F3幅值降為之前的四分之一，即83.325 N，F(xiàn)1、F2仍為333.3 N時，其環(huán)向聲輻射指向性如圖8(b)；若F1幅值降為之前的四分之一，即83.325 N，F(xiàn)2、F3仍為333.3 N時，其環(huán)向聲輻射指向性如圖8(c)。

與對稱三點激勵的聲壓分布圖進行對比，可以得出如下結(jié)論：對三點其中任意一點采取耗能性的減振措施，其結(jié)構(gòu)的振動聲輻射指向性將發(fā)生變化，同時也會影響輻射聲壓的幅值，幅值最大值由之前的86 dB升高至88 dB；如果針對其中對稱兩點做減振處理，其聲輻射指向性更加明顯，幅值最大值大致不變，但聲輻射最小值明顯改觀，由之前的80 dB降為74 dB。由上述結(jié)果分析得到，對于算例中的球殼結(jié)構(gòu)，盲目地采用減振措施并不一定能夠降低振動聲輻射的大小；如果對于軸向所夾角中線方向有特定的較低聲輻射要求，則可采用于F2和F3位置同時采用減振措施，否則減振的投入相較于得到的效果不理想。

2.2.3 改變?nèi)c激勵方向結(jié)構(gòu)輻射聲場計算

對于三點激勵模式，如將圖6(c)所示的激勵方向變?yōu)閳D9所示，F(xiàn)1和F2平行于x軸，其環(huán)向聲輻射分布如圖10所示，可以看出其在環(huán)向10 m處聲壓整體降低，最大聲壓級可以控制在85 dB以內(nèi)，y軸軸向方向聲壓級小于80 dB，滿足設(shè)計要求的作業(yè)過程低噪聲影響。

3 結(jié)論

本研究通過自編程序?qū)崿F(xiàn)了基于邊界元方法求解結(jié)構(gòu)振動聲輻射，并輔以脈動球表面聲壓和遠場輻射聲壓數(shù)值算例對比解析解驗證了聲學(xué)邊界元程序的有效性；計算了不同激勵模式下的鋼制球殼結(jié)構(gòu)在空氣中的輻射聲場、分析其分布特性，并通過改變激勵力的大小與方向從而實現(xiàn)輻射聲場的減小以滿足工程要求。復(fù)雜結(jié)構(gòu)在激勵源作用下產(chǎn)生的振動聲輻射熱區(qū)分布對于設(shè)備的安裝布置等工作至關(guān)重要，本文采用的有限元/邊界元耦合計算方法可有效實現(xiàn)復(fù)雜結(jié)構(gòu)受激聲輻射特性分析。