基于LMS VirtualLab的剛性球散射計(jì)算

朱 菡

(中國船舶集團(tuán)有限公司第七一三研究所,河南 鄭州 450000)

數(shù)值仿真技術(shù)在現(xiàn)代工程中得到廣泛應(yīng)用,在聲學(xué)仿真領(lǐng)域,球體的散射聲場仿真具有嚴(yán)格的理論解析解,經(jīng)常作為標(biāo)準(zhǔn)模型用于驗(yàn)證各種方法和理論,同時(shí),由于球體本身可以作為某些特定問題的簡化物理模型,常被作為簡化模型進(jìn)行實(shí)際問題的仿真計(jì)算。因此,球體的聲散射仿真研究一直是人們關(guān)注的問題。本文首先對剛性球體散射聲場遠(yuǎn)場理論解析解的推導(dǎo)過程進(jìn)行了說明,并計(jì)算了剛球遠(yuǎn)場散射聲壓,其次采用直接邊界元方法,應(yīng)用Virtual.Lab軟件對水下剛性球體散射聲場進(jìn)行仿真計(jì)算,并將計(jì)算結(jié)果與理論解析解進(jìn)行比較,得到了令人滿意的結(jié)果。

1 理論解析解

1.1 波動方程及其解

如圖1所示,剛性不動球位于無限流體介質(zhì)中,半徑為,表面光滑,在沿軸方向傳播的平面波(-)作用下,其散射聲場聲壓在球坐標(biāo)系中滿足波動方程：

圖1 平面波在球面上的散射

(1)

式中,、和為球坐標(biāo)系中的坐標(biāo)變量,=為波數(shù),等于入射聲波圓頻率和球體周圍介質(zhì)中的聲速之比?？紤]球體的對稱性,且入射波和散射波關(guān)于軸對稱,因此,方程應(yīng)與變量無關(guān),式(1)簡化為

(2)

采用分離變量法求解上述方程可得

(3)

式中,為待定常數(shù),由邊界條件確定。

1.2 邊界條件

剛性球體表面介質(zhì)質(zhì)點(diǎn)徑向振速為零,即

(4)

式中,是球體周圍介質(zhì)中的密度,是介質(zhì)中的總聲壓,即入射聲壓和散射聲壓之和,是介質(zhì)質(zhì)點(diǎn)振速的徑向分量,即入射波引起的介質(zhì)質(zhì)點(diǎn)振速的徑向分量和散射波引起的介質(zhì)質(zhì)點(diǎn)振速的徑向分量之和。

用勒讓德函數(shù)(Legendre Function)和球貝塞爾函數(shù)(Spherical Bessel Function)表示入射波,并利用關(guān)系式:

(5)

可求得常數(shù)為

(6)

式中，()為階球貝塞爾函數(shù)。

1.3 散射聲場的遠(yuǎn)場解

將式(6)代入式(3)并加入時(shí)間因子-得到散射聲表達(dá)式為

(7)

式(7)是散射聲場的一般表達(dá)式,但人們一般更關(guān)心它的遠(yuǎn)場特性,因此,可用球漢克爾函數(shù)(Spherical Hankel Function)在大宗量條件下的漸近展開式:

(8)

代入式(7)得

(9)

若記

(10)

則式(9)轉(zhuǎn)化為

(11)

1.4 解析解計(jì)算

上文中提到的指向性函數(shù)()是用于描述剛球散射聲場的一個(gè)重要參量,因此,為計(jì)算剛球的指向性函數(shù),計(jì)算模型選取半徑為1 m的剛性不動球,入射聲波選取單位波幅平面波,取計(jì)算階數(shù)=100,無因次系數(shù)={1,2,3,4,5,6,7,8},聲速=1 500 m/s,為滿足遠(yuǎn)場條件,需要計(jì)算獲得距剛球之外的聲壓,其中,為剛球尺度,為聲波波長,因此,計(jì)算=50 m處的散射聲壓,并利用球面衰減規(guī)律歸算到剛球等效聲中心1 m處,如圖2所示。

圖2 不同ka值下的剛球指向性圖

2 邊界元方法

2.1 Helmholtz積分方程

剛球聲散射計(jì)算屬于外部問題,即所要計(jì)算的聲場變量都在剛球封閉表面外部,滿足Helmholtz方程:

(12)

式中,為聲壓,為三維自由空間格林函數(shù):

(13)

如圖3所示,剛球表面為,在體積為表面為Σ半徑無限大的區(qū)域內(nèi),除點(diǎn)為格林函數(shù)的奇異點(diǎn)外,其余各點(diǎn)的聲壓在這個(gè)空間內(nèi)是平穩(wěn)和非奇異的,因此,可環(huán)繞點(diǎn)做半徑→0的表面為的小球,在表面積為++Σ的空間-內(nèi)就不再有奇異性。

圖3 計(jì)算區(qū)域定義

將式(12)代入格林公式,則有

(14)

(15)

式中,表示球體外部空間,表示球體表面。

2.2 邊界元方法

(16)

式中,(,)為單元形函數(shù),當(dāng)單元為四節(jié)點(diǎn)四邊形時(shí),有

(=1,2,3,4)

(17)

式中,和為單元第個(gè)頂點(diǎn)的局部坐標(biāo),其值為

(18)

將式(16)、(17)代入離散后的式(15),經(jīng)整理得

(19)

式中,[]和[]為系數(shù)矩陣,若表示球面聲壓,則式(19)給出了聲壓及其法向?qū)?shù)的關(guān)系,通過求解該微分方程組即可得到剛球表面聲壓分布,進(jìn)而得到空間內(nèi)任意一點(diǎn)的散射聲壓值。

2.3 Virtual.Lab Acoustics聲學(xué)仿真計(jì)算

平面波作用下,剛球散射聲場可以采用Virtual.Lab Acoustics聲學(xué)模塊進(jìn)行仿真計(jì)算,本文建立剛球四邊形邊界元網(wǎng)格,如圖4所示,網(wǎng)格劃分尺寸通常取波長的1/6,即

圖4 邊界元網(wǎng)格

(20)

式中,為網(wǎng)格最大尺寸,為聲波波長,為頻率,為水中聲速。

分別計(jì)算={1，2，3，4，5，6，7，8}時(shí)的散射聲場,圖5給出了剛球周圍2 m范圍內(nèi)近場的聲壓干涉圖。

圖5 剛球聲散射近場聲壓干涉圖

可以看出,邊界元算法能夠?qū)傂郧虮砻媛晧汉颓蝮w周圍聲場分布給出較好的計(jì)算結(jié)果,不論球體的聲學(xué)“照亮區(qū)”還是“影區(qū)”,算法均能給出準(zhǔn)確的計(jì)算結(jié)果,同時(shí),隨著ka值即計(jì)算頻率的增大,剛性球近場聲壓干涉條紋隨之加密,能夠較好地模擬水下目標(biāo)近場聲壓分布。

為了與上文中的理論解析解進(jìn)行比較,本文需要計(jì)算剛球指向性函數(shù),為此利用軟件自帶的Field Point Meshes功能下的Directivity Field Point Mesh功能,在距球心50 m處每隔1°設(shè)置一個(gè)計(jì)算場點(diǎn),用于提取該處的聲壓。將提取到的聲壓根據(jù)球面衰減規(guī)律歸算到剛球等效聲中心1 m處,并與解析解進(jìn)行比較,如圖6所示。可以看出,數(shù)值算法的計(jì)算結(jié)果與理論解析解吻合較好,能夠準(zhǔn)確地給出剛性球不同ka值下的指向性,證明該方法能夠準(zhǔn)確模擬任意分置角下水下目標(biāo)收發(fā)分置聲散射強(qiáng)度。

圖6 遠(yuǎn)場聲壓指向性數(shù)值解與解析解比較

3 結(jié)束語

本文利用LMS Virtual.Lab軟件,采用直接邊界元法計(jì)算水下剛性球體在平面聲波作用下的指向性函數(shù),并與理論解析解進(jìn)行比較,得到如下結(jié)論:

1)無論解析解還是邊界元方法數(shù)值解,計(jì)算時(shí)需滿足遠(yuǎn)場條件,否則,無法得到穩(wěn)定的、吻合度較高的計(jì)算結(jié)果;

2)采用Virtual.Lab Acoustics聲學(xué)模塊計(jì)算遠(yuǎn)場聲壓時(shí),當(dāng)計(jì)算點(diǎn)輸出聲壓數(shù)值較小時(shí),需要提高軟件參數(shù)設(shè)置中默認(rèn)聲壓輸出精度,否則,導(dǎo)出的計(jì)算結(jié)果數(shù)據(jù)精度將無法滿足需要;

3)邊界元法能夠準(zhǔn)確計(jì)算水下目標(biāo)的表面聲壓分布、近場聲場分布以及收發(fā)分置聲散射強(qiáng)度,可用于水下目標(biāo)聲散射特性預(yù)報(bào)。