彈性障板下水聲基陣輻射聲場計算

張俊，陳亞林，呂林夏

（1．中國船舶重工集團(tuán)公司第七○五研究所，陜西西安710075; 2．水下信息與控制重點實驗室，陜西西安 710075）

0 引言

在聲輻射計算領(lǐng)域，比較成熟的數(shù)值方法是有限元法和邊界元法。有限元法在水聲領(lǐng)域的應(yīng)用主要是水聲換能器的設(shè)計和分析，能很準(zhǔn)確地計算出換能器的諧振頻率，導(dǎo)納曲線發(fā)射指向性圖等。其在壓電、結(jié)構(gòu)、聲學(xué)耦合分析方面的準(zhǔn)確性已經(jīng)得到了驗證。但是在聲輻射計算方面，有限元法也存在不足。主要表現(xiàn)在:一是計算量比較大，對計算機(jī)的要求比較高;二是對于工程中常見的無限外空間聲輻射問題，有限元法不能有效地模擬截止邊界，從而導(dǎo)致計算誤差比較大。邊界元方法是由Helmholtz邊界積分方程并且結(jié)合邊界元的離散技術(shù)發(fā)展起來的一種數(shù)值計算方法。它把問題的求解維數(shù)降低一維，并且把解析解和數(shù)值解相結(jié)合，具有較小的計算量和較高的計算精度。同時邊界元法適合于無限域的聲學(xué)計算問題，這彌補(bǔ)了有限元法在這方面的不足。

在聲學(xué)基陣聲輻射方面，Audoly［1］和Yokoyama［2］分別利用改進(jìn)的Helmholtz邊界積分方程計算了有限大障板上平面基陣的輻射聲場。何正耀［3］用邊界元軟件Sysnoise計算了剛性障板條件下共形陣的輻射聲場，得到了陣元間的互輻射阻抗和基陣的遠(yuǎn)場指向性。但是他們分析計算的都是在剛性障板下的基陣聲輻射問題。對于有限大的彈性障板的特性，只采用邊界元法無法實現(xiàn)計算目的。本文利用有限元法在結(jié)構(gòu)模態(tài)分析上的準(zhǔn)確性和邊界元法在處理無限域聲輻射問題上的優(yōu)勢，采用有限元和邊界元相結(jié)合的方法，建立彈性障板下平面基陣的結(jié)構(gòu)邊界元模型，對其遠(yuǎn)場指向性進(jìn)行計算，并且對比分析彈性、柔性和剛性障板對整個基陣遠(yuǎn)場指向特性的影響。為工程中基陣的工程設(shè)計提供指導(dǎo)。

1 理論模型

1.1 聲學(xué)有限元理論模型

均勻的理想流體中三維波動方程為［3］:

用有限元方法對聲波方程進(jìn)行離散，同時考慮流固耦合處能量的損耗，得到離散化的有限元聲波方程:

式（3）與式（2）相結(jié)合，就得到完整的求解聲學(xué)問題的流固耦合控制方程為:

1.2 聲學(xué)邊界元理論模型

對于簡諧波聲壓為正弦時間函數(shù)，設(shè)波動頻率為ω，則聲壓p可以設(shè)為

代入三維波動方程式（1）中，得到:

在一個由封閉表面S所圍成的體積V中，對于2個在S上和V中都有一階和二階連續(xù)有界偏導(dǎo)數(shù)的函數(shù)φ和ψ。根據(jù)格林公式可以把體積分轉(zhuǎn)化為面積分:

式（7）是格林第二積分方程，它可以用來解決三維空間的聲輻射問題。如果振動體的體積為V，表面積為S，設(shè)格林公式中的φ和ψ分別為:

結(jié)合式（6）和式（9），考慮場點在體積V內(nèi)部、表面S上和外部3種情形，并且假設(shè)輻射表面光滑，可以推導(dǎo)出開放空間內(nèi)聲輻射問題的赫姆霍茨邊界積分方程:

當(dāng)場點X在輻射體外部時，C（X）=1;當(dāng)場點X在輻射體表面時，C（X）=1/2;當(dāng)場點X在輻射體內(nèi)部時，C（X）=0。

本文應(yīng)用當(dāng)前國際上流行的有限元軟件Ansys和邊界元聲學(xué)軟件Sysnoise聯(lián)合建立彈性障板下平面基陣的結(jié)構(gòu)邊界元聲學(xué)混合模型。首先用Ansys軟件對聲學(xué)結(jié)構(gòu)進(jìn)行精確的網(wǎng)格劃分和準(zhǔn)確地模態(tài)分析，然后在SysnoiseE軟件中，利用Ansys軟件建立的網(wǎng)格模型和模態(tài)參數(shù)建立結(jié)構(gòu)邊界元計算模型，計算彈性障板下平面基陣的水平軸面內(nèi)的聲壓分布和遠(yuǎn)場指向性。

2 計算實例

2.1 圓柱形彈性障板下37元平面陣聲場計算

建立37個陣元的平面陣模型，陣元位置排列如圖1所示。陣元輻射表面為正方形薄板，邊長a=3 cm，厚度b=0.3 cm，陣元垂直和水平間距相同。障板為一圓柱殼，障板高h(yuǎn)=20 cm，底面半徑R=13 cm，殼體厚度D=0.5 cm。陣元材料為硬鋁，密度ρ=2 700 kg/m3，楊氏模量E=6.85×1010N/m2，泊松比σ=0.34，聲速c=6 260 m/s;障板材料為玻璃鋼，密度ρ=1 750 kg/m3，楊氏模量E=1.8×1010N/m2，泊松比σ=0.14，聲速c=2 510 m/s。障板阻抗ρc= 4 392 500 Pa·s/m，與水的阻抗可以相比擬，可以看成是彈性障板。

圖1 基陣邊界元模型Fig.1The BEM model of arrays

建立基陣的幾何模型，用有限元法進(jìn)行網(wǎng)格劃分，對基陣進(jìn)行結(jié)構(gòu)模態(tài)計算，頻率f范圍為10～30 kHz。共得到526個結(jié)構(gòu)模態(tài)。表1為f=20 kHz前后的10個模態(tài)。

把有限元網(wǎng)格模型和模態(tài)參數(shù)分別保存為網(wǎng)格數(shù)據(jù)文件和Ansys文件。首先把網(wǎng)格數(shù)據(jù)文件導(dǎo)入到邊界元軟件Sysnoise中建立間接邊界元模型，如圖1所示。然后把保存模態(tài)參數(shù)的Ansys文件導(dǎo)入到軟件中，軟件會自動生成間接邊界元和結(jié)構(gòu)有限元耦合模型。節(jié)點振速設(shè)為v=1 m/s，設(shè)置障板和陣元的材料參數(shù)。計算得到頻率f=20 kHz時基陣的0°方向上遠(yuǎn)場指向性分別如圖2中所示，實線為Sysnoise計算值，虛線是由平面波模型計算所得到的理論值。圖3為當(dāng)頻率f=20 kHz時，主波束在0°方向時軸平面上的聲壓分布圖，所顯示平面區(qū)域長寬均為2 m。

由結(jié)果可知，－3 dB帶寬為19.03°，旁瓣級為21.34 dB。由圖中可以看出，這與理論值是比較接近的。由于障板的彈性特性和網(wǎng)格劃分不對稱的影響，由Sysnoise所計算的基陣遠(yuǎn)場指向性左右旁瓣已經(jīng)不具有理想的對稱性，特別是接近障板邊緣處的旁瓣已經(jīng)表現(xiàn)的比較凌亂，與理論值相差比較大。從聲壓分布圖中可以看出近場處聲波的衍射現(xiàn)象比較嚴(yán)重。

2.2 不同阻抗障板下37元陣遠(yuǎn)場輻射指向性對比

在水聲場中，當(dāng)障板的聲阻抗分別遠(yuǎn)小于、遠(yuǎn)大于水的聲阻抗或者與水的聲阻抗相當(dāng)時，就可以分別看成是柔性、剛性和彈性障板。它們分別對應(yīng)著聲學(xué)中的絕對軟、絕對硬和混合邊界條件。有如下的特性:

1）柔性障板，這時障板邊界上的聲壓為0，即p（x，y，z，t）=0。其物理意義為界面上的任何點上，不論時間t取何值，聲壓p總為0。

2）剛性障板，聲波不能進(jìn)入該介質(zhì)中，此時邊界上介質(zhì)質(zhì)點的法向振速應(yīng)為0。用數(shù)學(xué)模型表達(dá)為

式中:系數(shù)a，b均為常數(shù);s為邊界。當(dāng)f（s）=0時，則稱其為阻抗邊界條件。

建立37元平面陣及障板模型，幾何參數(shù)與上一小節(jié)所建立的模型一致。障板材料阻抗屬性分別設(shè)置為柔性、剛性和彈性。根據(jù)上面3種情況，對于柔性障板和剛性障板，選用間接邊界元（DBEM）計算模型，分別設(shè)置障板節(jié)點聲壓為0和節(jié)點法向振速為0;對于彈性障板模型，同樣選用玻璃鋼，采用結(jié)構(gòu)有限元和邊界元混合計算的方法。計算頻率為歸一化頻率f=15 kHz。圖4和圖5分別為束控方向0°時，3種障板下的遠(yuǎn)場指向性圖。是表示界面沿外法線方向的偏導(dǎo)數(shù)。

3）彈性障板，在界面上聲壓和振速是線性組合關(guān)系，其數(shù)學(xué)表達(dá)式為

圖4 0°束控時基陣水平遠(yuǎn)場發(fā)射指向性圖Fig.4The far-field directivity of arrays when 0°beam control

由圖6可以看出，當(dāng)主波束束控在0°角上時，柔性、剛性、彈性3種障板下，基陣的主波束吻合的非常好，第一旁瓣也基本一致，第二、三旁瓣有差異。這說明不同材料障板對基陣的影響主要集中在靠近障板邊緣處。當(dāng)主波束束控在30°和45°時，剛性障板和柔性障板下基陣遠(yuǎn)場指向性趨勢一致。但是彈性障板下基陣的主瓣和旁瓣都有較大的改變，30°和45°束控時均在－20°方向出現(xiàn)了較高的旁瓣。30°束控時主瓣寬度較之剛性和柔性障板條件下變小，45°束控時第一旁瓣已經(jīng)不明顯，主瓣下端變胖。

3 結(jié)語

1）本文通過邊界元和有限元相結(jié)合的數(shù)值方法計算了彈性障板下水聲基陣的輻射遠(yuǎn)場指向性和近場聲壓分布，所得結(jié)果與解析解相吻合。但是解析解不能很好地考慮障板條件對基陣的影響，特別是復(fù)雜的障板條件下基陣的聲場計算問題，數(shù)值方法比解析方法有明顯的優(yōu)勢。

2）由結(jié)果可以看出，障板對基陣輻射聲場的影響主要體現(xiàn)在對其遠(yuǎn)場指向性旁瓣的影響和對其近場聲場的影響。越接近障板邊緣，旁瓣的變形越大。在近場處，聲波的衍射現(xiàn)象也比較嚴(yán)重。特別是彈性障板條件下，這種影響更為嚴(yán)重，特別是做波束偏移時，主瓣和旁瓣都有比較明顯的改變。這一點，在基陣的設(shè)計時值得考慮。

3）邊界元聲學(xué)數(shù)值方法中，模型網(wǎng)格的劃分質(zhì)量對結(jié)果也有明顯的影響。對于本文中建立的圓柱體彈性障板，網(wǎng)格單元的分布不對稱，對其遠(yuǎn)場指向性旁瓣也有一定的影響。優(yōu)化網(wǎng)格劃分，進(jìn)行精確快速的建模和計算是有待進(jìn)一步開展的工作。

4）本文所用到的分析方法不只適用于平面陣，對于陣元排列和障板結(jié)構(gòu)更為復(fù)雜的體積陣也同樣可以適用。這比平面波理論模型更接近工程實際，為聲吶基陣的工程設(shè)計，特別是在多種材料混合使用下基陣設(shè)計提供了一種預(yù)報基陣聲學(xué)性能的數(shù)值方法。

［1］CHRISTIAN A．Some aspedts of acoustic interactions in sonar transducer arrays［J］．J．Acoust．Soc．Am．1991，89 （3）:1428－1433．

［2］YOKOYAMA T，HENMI M．Effects of mutual interactions on a phased transducer array［J］．Jpn．J．Appl．phys．，1998，38（5）:3356－3360．

［3］SYSNOISE Rev 5.5 Manuals［M］．LMS Company，2005．

［4］杜功煥．聲學(xué)基礎(chǔ)［M］．南京:南京大學(xué)出版社，2001．

DU Gong-huan．The Foundation of Acoustics［M］．Nanjing: Nanjing University Press，2001．

［5］周福洪．水聲換能器及基陣［M］．北京:國防工業(yè)出版社，1984．

ZHOU Hong-fu．Underwater Acoustics Transducer and Arrays［M］．Beijing:National Defence Industry Press，1984．

［6］MOAVENI S．有限元分析［M］．北京:電子工業(yè)出版社，2008．

MOAVENI S．Finiteelementanalysis［M］．Beijing: Electronic Industy Press，2008．

［7］何祚鏞．結(jié)構(gòu)振動與聲輻射［M］．哈爾濱:哈爾濱工程大學(xué)出版社，2001．

HE Zuo-yong．Structure Vibration and Sound Radiation［M］．Harbin:Harbin Engineering University Press，2001．

［8］CISKOWSKI R D，BREBBIA C A．Boundary Element MethodsinAcoustics［M］．SouthamptonBoston Computational Mechanics Publications，1991．

［9］COPLEY L G．Integral equation method for radiation from vibrating bodies［J］．Journal of the Acoustical Society of America，1967，41（4）:807－816．