朱 菡

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

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

(4)

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

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

(5)

可求得常數(shù)為

(6)

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

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

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

(7)

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

(8)

代入式(7)得

(9)

若記

(10)

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

(11)

1.4 解析解計算

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

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

2 邊界元方法

2.1 Helmholtz積分方程

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

(12)

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

(13)

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

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

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

(14)

(15)

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

2.2 邊界元方法

(16)

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

(=1,2,3,4)

(17)

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

(18)

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

(19)

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

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

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

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

(20)

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

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

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

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

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

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

3 結(jié)束語

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

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

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

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