基于邊界元方法的物探船設(shè)備聲輻射特性分析

曹樹杰，宋春旺，米 崢

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

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

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

1.1 理論分析

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

▽2p+k2p=0

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

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

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

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

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

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

1.2 數(shù)值分析

1) 離散表面

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

2) 積分方程離散分析

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

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

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

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

3) 奇異積分問題

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

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

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

2 數(shù)值分析算例

2.1 脈動(dòng)球驗(yàn)證程序有效性

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

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

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

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

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

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

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

波數(shù)解析解聲壓/Pa96個(gè)單元?jiǎng)澐謹(jǐn)?shù)值解聲壓/Pa無chief方法96個(gè)單元?jiǎng)澐謹(jǐn)?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ù)脈動(dòng)球以210 Hz的頻率振動(dòng)時(shí)，不同場(chǎng)點(diǎn)處聲壓幅值的600個(gè)單元?jiǎng)澐謹(jǐn)?shù)值解與解析解誤差與其距球心距離的關(guān)系曲線，可以看出，距球心大于1.25 m的場(chǎng)點(diǎn)的計(jì)算結(jié)果誤差已足夠小，說明本程序可有效計(jì)算結(jié)構(gòu)振動(dòng)輻射聲場(chǎng)；對(duì)于距離脈動(dòng)球表面很近的場(chǎng)點(diǎn)，其結(jié)果存在奇異性，需采取其他方法消除。

圖5給出的是上述參數(shù)脈動(dòng)球的遠(yuǎn)場(chǎng)(r=10 m)輻射聲壓幅值(解析解、600個(gè)單元?jiǎng)澐謹(jǐn)?shù)值解)隨波數(shù)變化曲線，可以看出，遠(yuǎn)場(chǎng)聲壓數(shù)值計(jì)算結(jié)果與解析解大致相同，證明此程序能夠有效計(jì)算輻射聲場(chǎng)。

2.2 不同激勵(lì)模式的球殼輻射聲場(chǎng)分析

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

2.2.1對(duì)球殼進(jìn)行3種激勵(lì)方式的輻射聲場(chǎng)分析

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

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

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

圖6 三類激勵(lì)模式

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

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

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

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

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

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

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

3 結(jié)論

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