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

魏克難，李 威，雷 明，柴應彬

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

0 引 言

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

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

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

1 聲學邊界元理論

1.1 邊界元法簡介

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

當聲波入射物體時，如果結構和流體之間的相互作用比較大，就必須考慮結構和流體之間的相互作用，本文計算彈性體的聲散射就屬于這一類情況，用到的是基于有限元和邊界元相結合的耦合聲學邊界元理論。

1.2 耦合聲學邊界元分析理論

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

式中Na為與邊界Ωa上na個節點上聲壓pi和法向速度vni相關的全局形函數。

在邊界元中，還有一些不知道聲壓和振動速度的節點，例如節點b，也就是ra= rb，則有

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

對于耦合邊界元，聲場分布和結構的振動位移是同時計算出來的。對于邊界元網格，可以分為2個部分，一個是與結構網格耦合的部分，包含an1個節點，另一部分是不參與耦合的部分，包含an2個節點(an1+an2= an)，這樣邊界元上的聲壓p(ra)和速度v(ra)可以寫為

式中Na1和Na2分別為與an1節點和an2節點相關的形函數。由于聲壓也作用在結構上，它作為一個載荷，同樣也可以引起結構的振動，這樣結構部分的動力學方程為

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

式中:nse為耦合的結構網格的單元數量;Ns為結構網格的形函數;{ ne}為單元的法向方向;負號是因為結構網格的單元法向方向與邊界元的單元法向方向相反。

對于耦合邊界元法，在結構與聲音耦合面Ωs處，結構網格的節點和聲音網格的節點需要重合。在耦合面Ωs處，由結構網格的振動速度與聲音網格的振動速度的連續性與位移{u}、速度{v}之間的關系{v}=jω{u}，可得

式中{ne}為結構網格的法線方向，負號表示結構的法線方向與聲學網格的法線方向相反。將式(8)代入式(3)可得

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

式中:

2 算例分析

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

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

圖1 球殼的網格模型Fig.1 Shell mesh model

定該球殼的密度ρ = 7 860 kg/m3，楊氏模量E= 200 Gpa，泊松比ν = 0.266。因為是在水中，所以流體參數的設置，密度ρ水= 1 000 kg/m3，聲音傳播的速度c水= 1 500 m/s，求解平面波入射下的目標強度隨頻率的變化，與解析解對比分析圖如圖2所示。

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

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

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

按照參考文獻［1］中的潛艇規格圖 (見圖3)，在MSC.Patran 中建BeTSSi-Sub的網格模型，如圖4所示，此模型有6 656個單元，6 650個節點。

圖3 BeTSSi-Sub 規格圖Fig.3 Dimensions of BeTSSi-Sub

圖4 BeTSSi-Sub的網格模型Fig.4 BeTSSi-Sub mesh model

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

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

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

從圖5 可明顯看出，剛性潛艇的目標強度大于彈性(鋼)的目標強度，這是由于所謂剛體散射是指物體入射聲波的作用下既不發生形變，聲波也透不到物體內部，而彈性目標的散射聲波則會進入物體的內部，進而會在物體的內部產生更復雜的散射和折射，同時也會引起目標的變形和振動。

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

從潛艇的網格圖可以看出，潛艇不是規則的幾何體，尤其是舵和上層建筑，因此為了探究這2個形狀不規則的部分對潛艇整體目標強度的影響，分別建立無上層建筑的BeTSSi-Sub 網格模型和無舵的BeTSSi-Sub 網格模型，如圖6和圖7所示。

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

圖7 無舵的BeTSSi-Sub 網格模型Fig.7 BeTSSi-Sub mesh model with no rudder

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

圖8 不同形狀的BeTSSi-Sub 正橫反向散射目標強度對比圖Fig.8 Backscattering of the elastic BeTSSi-Sub with different shapes

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

3 結 語

從本文的數值計算算例可以看出，這種結合商業軟件的耦合邊界元數值計算方法，不僅可以解決簡單的彈性目標的聲散射特性問題，對于水下像潛艇一樣有復雜結構的水下目標，也能準確且成功預報其目標強度。

通過在低頻段對不同材料，不同形狀的BeTSSi-Sub的目標強度的計算，對比發現:剛性BeTSSi-Sub的目標強度大于彈性BeTSSi-Sub的目標強度;在正橫方向上，舵對BeTSSi-Sub 目標強度的影響作用比生層建筑更加明顯。

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

［2］李威，趙耀，張濤，等.水下剛硬體聲散射場計算及分析［J］.聲學技術，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］湯渭霖.用物理聲學方法計算界面附近目標的回波［J］.聲學學報，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］.北京:機械工業出版社，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 聲學仿真計算高級應用實例［M］.北京:國防工業出版社，2010.LI Zeng-gang，ZHAN Fu-liang，Virtual.lab.Advanced acoustic simulation application examples［M］.Beijing:National Defense Industry Press，2010.

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