基于耦合邊界元法的水下目標(biāo)低頻聲散射特性

魏克難，李 威，雷 明，柴應(yīng)彬

(華中科技大學(xué)，湖北 武漢430074)

0 引 言

水下復(fù)雜目標(biāo)的聲散射特性直接影響著武器裝備的水下作戰(zhàn)性能，它對(duì)于水中目標(biāo)的識(shí)別、隱蔽和反隱身起著決定性作用，一直以來倍受人們的關(guān)注，有著重要的軍事意義。因此，對(duì)復(fù)雜目標(biāo)的散射特性進(jìn)行精確而快速的預(yù)報(bào)一直是聲學(xué)研究的熱點(diǎn)之一。目前國內(nèi)外圍繞這個(gè)主題，通過研究提出了一系列的理論解法和近似解法，從而對(duì)水下目標(biāo)的聲散射特性有了更深入的了解。

目標(biāo)的聲散射問題理論上是解決一個(gè)數(shù)學(xué)物理問題，即求解三維流體空間中的目標(biāo)受聲激勵(lì)產(chǎn)生的滿足波動(dòng)方程、表面邊界條件和輻射條件的散射聲場。圍繞這個(gè)問題，對(duì)水下目標(biāo)聲散射特性研究的方法主要分為理論解法和近似解法。理論解法主要包括積分方程法和分離變量法;近似解法主要包括統(tǒng)計(jì)能量分析法、T 矩陣法、物理聲學(xué)方法及時(shí)域有限差分法［2－4］。國內(nèi)外也開發(fā)出了很多預(yù)報(bào)潛艇目標(biāo)強(qiáng)度的軟件，例如瑞典的M.Almgren和M.Nodin［5］開發(fā)的SUBTAS 軟件，美國海軍實(shí)驗(yàn)室(NRL)的Joseph Shirron［6］開發(fā)的有限元－無限元計(jì)算軟件SONAX 等。

本文是采取聲振耦合邊界元的方法，采用商業(yè)軟件MSC.Patran 建立目標(biāo)網(wǎng)格模型，采用聲學(xué)分析軟件LMS.Virtual.Lab 進(jìn)行數(shù)值分析計(jì)算［7－8］，對(duì)比探究剛性BeTSSi-Sub和彈性BeTSSi-Sub的全向散射特性，并且通過改變BeTSSi－Sub的外形，探究了在低頻段舵和上層建筑潛艇這2個(gè)特殊的部分對(duì)潛艇目標(biāo)強(qiáng)度的影響。

1 聲學(xué)邊界元理論

1.1 邊界元法簡介

邊界元法是將問題的微分方程變成邊界上的積分方程，然后引入邊界上的有限個(gè)單元將積分方程離散，得到只含有邊界上節(jié)點(diǎn)未知量的方程組，然后進(jìn)行數(shù)值求解。由于邊界元法是將區(qū)域上的控制方程轉(zhuǎn)化為沿著區(qū)域邊界的積分方程，因此它只需要定義邊界上的單元，結(jié)合邊界條件求解，這樣就使處理問題的維數(shù)降低一維，可以有效地將復(fù)雜的三維幾何模型簡化為二維圖形。與同樣復(fù)雜的完整三維有限元模型相比，這是更小、易于創(chuàng)建、易于檢驗(yàn)，并易于處理的模型，這些簡化模型可以在更短的時(shí)間內(nèi)得出結(jié)果。但邊界元法所建立方程組的系數(shù)矩陣稠密，而且一般非對(duì)稱，矩陣元素分量的計(jì)算量很大，增加了計(jì)算時(shí)間。由于邊界元法所利用的微分算子基本解能自動(dòng)滿足無限遠(yuǎn)處的條件，因而邊界元法特別便于無限域以及半無限域問題。

當(dāng)聲波入射物體時(shí)，如果結(jié)構(gòu)和流體之間的相互作用比較大，就必須考慮結(jié)構(gòu)和流體之間的相互作用，本文計(jì)算彈性體的聲散射就屬于這一類情況，用到的是基于有限元和邊界元相結(jié)合的耦合聲學(xué)邊界元理論。

1.2 耦合聲學(xué)邊界元分析理論

運(yùn)用邊界元法處理求解積分方程，邊界元需要把邊界Ωa離散成許多單元和節(jié)點(diǎn)，每個(gè)單元內(nèi)部任意點(diǎn)的聲壓和法向速度可以用屬于這個(gè)單元節(jié)點(diǎn)上的聲壓pi和法相速度vni與全局矩陣形函數(shù)Na來表示:

式中Na為與邊界Ωa上na個(gè)節(jié)點(diǎn)上聲壓pi和法向速度vni相關(guān)的全局形函數(shù)。

在邊界元中，還有一些不知道聲壓和振動(dòng)速度的節(jié)點(diǎn)，例如節(jié)點(diǎn)b，也就是ra= rb，則有

式中系數(shù)矩陣Ab和Bb都是(1 × na)的矩陣。

對(duì)于耦合邊界元，聲場分布和結(jié)構(gòu)的振動(dòng)位移是同時(shí)計(jì)算出來的。對(duì)于邊界元網(wǎng)格，可以分為2個(gè)部分，一個(gè)是與結(jié)構(gòu)網(wǎng)格耦合的部分，包含an1個(gè)節(jié)點(diǎn)，另一部分是不參與耦合的部分，包含an2個(gè)節(jié)點(diǎn)(an1+an2= an)，這樣邊界元上的聲壓p(ra)和速度v(ra)可以寫為

式中Na1和Na2分別為與an1節(jié)點(diǎn)和an2節(jié)點(diǎn)相關(guān)的形函數(shù)。由于聲壓也作用在結(jié)構(gòu)上，它作為一個(gè)載荷，同樣也可以引起結(jié)構(gòu)的振動(dòng)，這樣結(jié)構(gòu)部分的動(dòng)力學(xué)方程為

式中Ks，Cs和Ms分別為結(jié)構(gòu)有限元模型的剛度矩陣、阻尼矩陣和質(zhì)量矩陣;{ui}為結(jié)構(gòu)移位向量;{Fs}為作用在結(jié)構(gòu)網(wǎng)格上的外載荷(不包含聲壓載荷);Lc·{pi1}為聲壓作用在結(jié)構(gòu)網(wǎng)格上的載荷;Lc為(ns× na1)耦合矩陣，即

式中:nse為耦合的結(jié)構(gòu)網(wǎng)格的單元數(shù)量;Ns為結(jié)構(gòu)網(wǎng)格的形函數(shù);{ ne}為單元的法向方向;負(fù)號(hào)是因?yàn)榻Y(jié)構(gòu)網(wǎng)格的單元法向方向與邊界元的單元法向方向相反。

對(duì)于耦合邊界元法，在結(jié)構(gòu)與聲音耦合面Ωs處，結(jié)構(gòu)網(wǎng)格的節(jié)點(diǎn)和聲音網(wǎng)格的節(jié)點(diǎn)需要重合。在耦合面Ωs處，由結(jié)構(gòu)網(wǎng)格的振動(dòng)速度與聲音網(wǎng)格的振動(dòng)速度的連續(xù)性與位移{u}、速度{v}之間的關(guān)系{v}=jω{u}，可得

式中{ne}為結(jié)構(gòu)網(wǎng)格的法線方向，負(fù)號(hào)表示結(jié)構(gòu)的法線方向與聲學(xué)網(wǎng)格的法線方向相反。將式(8)代入式(3)可得

結(jié)合結(jié)構(gòu)有限元方程式(6)和邊界元式(11)，可以得到結(jié)構(gòu)有限元與邊界元的耦合方程為

式中:

2 算例分析

2.1 彈性球殼平面波入射下的散射特性分析

為了驗(yàn)證該方法的有效性，采用此方法計(jì)算無界流體中收發(fā)合置下彈性球殼的反向散射特性。首先在建立球殼的網(wǎng)格模型，定球殼的外徑為1 m，厚度為0.03 m，建立網(wǎng)格模型如圖1所示。

圖1 球殼的網(wǎng)格模型Fig.1 Shell mesh model

定該球殼的密度ρ = 7 860 kg/m3，楊氏模量E= 200 Gpa，泊松比ν = 0.266。因?yàn)槭窃谒?，所以流體參數(shù)的設(shè)置，密度ρ水= 1 000 kg/m3，聲音傳播的速度c水= 1 500 m/s，求解平面波入射下的目標(biāo)強(qiáng)度隨頻率的變化，與解析解對(duì)比分析圖如圖2所示。

圖2 彈性球殼反向散射特性對(duì)比圖Fig.2 Backscattering of the elastic shell

從圖2 可看出，用數(shù)值解和解析解有很好的吻合，初步論證了耦合邊界元法在計(jì)算水下目標(biāo)聲散射特性可行。

2.2 BeTSSi-Sub的全向散射特性分析

按照參考文獻(xiàn)［1］中的潛艇規(guī)格圖 (見圖3)，在MSC.Patran 中建BeTSSi-Sub的網(wǎng)格模型，如圖4所示，此模型有6 656個(gè)單元，6 650個(gè)節(jié)點(diǎn)。

圖3 BeTSSi-Sub 規(guī)格圖Fig.3 Dimensions of BeTSSi-Sub

圖4 BeTSSi-Sub的網(wǎng)格模型Fig.4 BeTSSi-Sub mesh model

對(duì)于邊界元模型來說，要求最大單元的編程要小于計(jì)算頻率波長的1/6，所以要計(jì)算到很高的頻率，勢必網(wǎng)格單元數(shù)量會(huì)急劇上升，計(jì)算時(shí)間會(huì)急劇增長，對(duì)計(jì)算機(jī)的要求也變高。所以根據(jù)本文所畫的網(wǎng)格模型，最大的計(jì)算頻率約為373 Hz。

從文獻(xiàn)［9］可知，潛艇目標(biāo)強(qiáng)度值的測量隨測量距離發(fā)生變化，所以為了得到穩(wěn)定可靠地測量結(jié)果，測量應(yīng)該在遠(yuǎn)場進(jìn)行。對(duì)于長度為L的潛艇而言，當(dāng)入射波長為λ，即測量距離r 要大于L2/λ，所以本文定義波源位置為距離潛艇幾何中心1 000 m處，計(jì)算收發(fā)合置情況下剛性和彈性(鋼)潛艇的全向散射時(shí)的目標(biāo)強(qiáng)度，定義從潛艇首部入射為0°，從潛艇尾部入射為180°，每隔2°計(jì)算1 次，計(jì)算頻率分別為100 Hz，200 Hz，300 Hz，定義彈性BeTSSi-Sub的材料參數(shù)為鋼，密度ρ = 7 860 kg/m3，楊氏模量E = 200 Gpa，泊松比ν = 0.266，鋼板厚度為35 mm，得到剛性和彈性(鋼)的目標(biāo)強(qiáng)度對(duì)比圖如圖5所示。

圖5 剛性和彈性(鋼)BeTSSi-Sub 全向散射目標(biāo)強(qiáng)度對(duì)比圖Fig.5 Omni-directional scattering of the rigid BeTSSi-Sub and the elastic (steel)BeTSSi-Sub

從圖5 可明顯看出，剛性潛艇的目標(biāo)強(qiáng)度大于彈性(鋼)的目標(biāo)強(qiáng)度，這是由于所謂剛體散射是指物體入射聲波的作用下既不發(fā)生形變，聲波也透不到物體內(nèi)部，而彈性目標(biāo)的散射聲波則會(huì)進(jìn)入物體的內(nèi)部，進(jìn)而會(huì)在物體的內(nèi)部產(chǎn)生更復(fù)雜的散射和折射，同時(shí)也會(huì)引起目標(biāo)的變形和振動(dòng)。

2.3 不同形狀的彈性(鋼)BeTSSi-Sub 散射特性的分析

從潛艇的網(wǎng)格圖可以看出，潛艇不是規(guī)則的幾何體，尤其是舵和上層建筑，因此為了探究這2個(gè)形狀不規(guī)則的部分對(duì)潛艇整體目標(biāo)強(qiáng)度的影響，分別建立無上層建筑的BeTSSi-Sub 網(wǎng)格模型和無舵的BeTSSi-Sub 網(wǎng)格模型，如圖6和圖7所示。

圖6 無上層建筑的BeTSSi-Sub 網(wǎng)格模型Fig.6 BeTSSi-Sub mesh model with no superstructure

圖7 無舵的BeTSSi-Sub 網(wǎng)格模型Fig.7 BeTSSi-Sub mesh model with no rudder

從圖中可看出，正橫方向(90°)是潛艇目標(biāo)強(qiáng)度最大的方位，一般我們最關(guān)心的也是正橫附近的目標(biāo)強(qiáng)度。對(duì)這3個(gè)不同形狀的潛艇網(wǎng)格模型，本文計(jì)算了正橫方向的反向散射，這3 種模型的材料均為鋼，ρ = 7 860 kg/m3，楊氏模量E = 200 Gpa，泊松比ν = 0.266，鋼板厚度為35 mm，計(jì)算結(jié)果如圖8所示。

圖8 不同形狀的BeTSSi-Sub 正橫反向散射目標(biāo)強(qiáng)度對(duì)比圖Fig.8 Backscattering of the elastic BeTSSi-Sub with different shapes

由圖8 可知，在正橫方向上，無舵的BeTSSi-Sub 在0～150 Hz的目標(biāo)強(qiáng)度比BeTSSi-Sub 大，而在150～300 Hz 比BeTSSi-Sub 小。上層建筑的存在對(duì)BeTSSi-Sub的目標(biāo)強(qiáng)度基本沒影響，在低頻段上，舵對(duì)潛艇正橫方向目標(biāo)強(qiáng)度的影響比上層建筑大。

3 結(jié) 語

從本文的數(shù)值計(jì)算算例可以看出，這種結(jié)合商業(yè)軟件的耦合邊界元數(shù)值計(jì)算方法，不僅可以解決簡單的彈性目標(biāo)的聲散射特性問題，對(duì)于水下像潛艇一樣有復(fù)雜結(jié)構(gòu)的水下目標(biāo)，也能準(zhǔn)確且成功預(yù)報(bào)其目標(biāo)強(qiáng)度。

通過在低頻段對(duì)不同材料，不同形狀的BeTSSi-Sub的目標(biāo)強(qiáng)度的計(jì)算，對(duì)比發(fā)現(xiàn):剛性BeTSSi-Sub的目標(biāo)強(qiáng)度大于彈性BeTSSi-Sub的目標(biāo)強(qiáng)度;在正橫方向上，舵對(duì)BeTSSi-Sub 目標(biāo)強(qiáng)度的影響作用比生層建筑更加明顯。

［1］NELL C W，GILROY L E.An improved BASIS model for the BeTSSi submarine ［R］.Canada;DRDC－AT LANTIC，2003.

［2］李威，趙耀，張濤，等.水下剛硬體聲散射場計(jì)算及分析［J］.聲學(xué)技術(shù)，2007，26(5):844－849.LI Wei，ZHAO Yao，ZHANG Tao，et al.Computation and analysis of the acoustic scattering field from submerged rigid objects［J］.Technical Acoustics，2007，26(5):844－849.

［3］湯渭霖.用物理聲學(xué)方法計(jì)算界面附近目標(biāo)的回波［J］.聲學(xué)學(xué)報(bào)，1999，24(1).TANG Wei-lin.Calculation of acoustic scattering from object near an interface using physical acoustic method［J］.Acta Acustica，1999，24(1).

［4］HASTING F D，SCHNEIDER J B，BROSEHAT S，Afinitedifferenee time-domain solution to scattering from a rough pressure-release surface［J］.Acoust.Soc Am，1997，102(6):3394－3400.

［5］ALMEGRON M，NODIN M.Acoustics target strength of submarines.PRAO'S，1991.

［6］SHIRRON J.Computational simulation of elastic scattering［D］.Naval Research Laboratory.Washington D.C

［7］龍凱，賈長治，李寶峰.Pratran2010與Nastran2010 有限元分析從入門到精髓［M］.北京:機(jī)械工業(yè)出版社，2011.LONG Kai，JIA Chang-zhi，LI Bao-feng.Finite element analysis from entry to the essence of pratran2010 and Nastran2010 machinery industry press［M］.Beijing:China Machine Press，2011.

［8］李增剛，詹福良.Virtual.lab.Acoustics 聲學(xué)仿真計(jì)算高級(jí)應(yīng)用實(shí)例［M］.北京:國防工業(yè)出版社，2010.LI Zeng-gang，ZHAN Fu-liang，Virtual.lab.Advanced acoustic simulation application examples［M］.Beijing:National Defense Industry Press，2010.

［9］劉伯勝.水聲學(xué)原理［M］.哈爾濱:哈爾濱大學(xué)出版社，2011.LIU Bo-sheng.Water acoustics ［M］.Harbin:Harbin University Press，2011.