基于Kirchhoff近似與曲面三角形網格的水下目標聲散射特性分析

薛亞強,彭子龍*,俞強,張春雨,周富霖,劉進偉

(1.江蘇科技大學 能源與動力學院, 江蘇 鎮江 212000; 2.上海海事大學 海洋科學與工程學院, 上海 200135;3.上海交通大學 海洋工程重點實驗室, 上海 200240)

0 引言

水下目標聲散射特性分析是水聲領域的熱點問題,也是實現目標探測和識別的前提,具有重要的理論價值和實際意義。對于形狀比較規則的目標,如球、無限長圓柱等,其聲散射特性可用解析法[1]求得。當目標的幾何形狀比較復雜時,一般采用數值解法,如有限元[2]、邊界元[3-4]和將二者結合的耦合方法[5]等,以及T矩陣法[6]、時域有限差分法[7]、板塊元方法[8]和聲束彈跳法[9]等近似解法。

基于Kirchhoff近似的板塊元方法最初用于雷達散射截面的計算,自引入水下目標聲散射研究以來,廣泛應用于大尺度模型的高頻回波特性預報[10]。板塊元方法的核心算法是通過傅里葉變換積分或Gordon積分[11]將面積分運算轉化為代數運算;該方法的基本原理與直接數值積分方法[12]相似,將目標表面劃分為若干網格,對所有網格的散射聲場求和得到目標總散射聲場。

Hiptmair等[13]采用Trefftz間斷Galerkin方法和局部細化網格研究了軟邊界散射體的聲散射問題。Chaillat等[14]將邊界元方法和自適應網格用于求解復雜形狀目標聲散射,具有較高的收斂速率。文獻[15-16]基于Kirchhoff近似對水下目標的高頻聲散射特性進行了研究,采用平面三角形單元離散目標表面。然而,對于表面曲率較大或形狀較為復雜的目標,如含有多個子結構和不規則曲面的水下航行器,若采用平面單元進行網格離散,需劃分大量的網格才能較準確地擬合目標外形。與平面網格相比,曲面網格可以更精確地逼近目標的幾何形狀。Wu等[17]采用曲邊四邊形及曲邊三角形單元離散圓板及開口圓柱殼,建立了目標輻射與散射的邊界元分析模型。Foote等[18]基于曲面網格模型對比了Kirchhoff近似與邊界元得到的目標強度,并指出高頻時二者的結果基本一致,但前者具有更快的計算速率。Ven?s等[19]采用高階樣條曲面單元和等幾何邊界元建立了潛艇的聲散射分析模型。

本文基于Kirchhoff近似和曲面三角形網格建立了水下目標收發合置聲散射特性分析模型,采用高斯-勒讓德求積方法對面積分進行直接數值計算。以剛性球、橢圓柱、球冠柱和Benchmark縮比模型等為例,比較了Kirchhoff近似模型與其他方法的計算結果,驗證了本文方法的準確性與高效性。

1 Kirchhoff高頻近似

由物理聲學方法或Kirchhoff高頻近似可以得到收發合置時,剛性目標的散射聲場勢函數滿足的Helmholtz積分公式[8]為

(1)

式中:k為波數;A為任意振幅;積分面S0為散射體的幾何亮區表面;r為面元dS到收發合置換能器M的距離;θ為入射聲線r與面元單位法向矢量n的夾角。如圖1所示,r0為換能器M到參考點的矢量,rs為面元dS到參考點的矢量。

圖1 Kirchhoff近似積分方法的坐標示意圖

遠場條件下,式(1)中被積函數的指數上取r=r0+Δr,r0為換能器M到參考點的距離,Δr為聲程差,式(1)分母上取r≈r0,可得

(2)

將遠場中的聲線r與r0看作平行,聲程差近似為Δr≈rs·r0/r0。式(2)中的積分可表示為

I=?S0e2ikrq·r0/r0n·r0/r0dS

(3)

目標位置處的入射波勢函數φi(r)=-(A/r0)·eikr0,由定義可得目標強度為

(4)

2 基于曲面單元的目標強度計算方法

首先對散射體的幾何表面進行網格劃分,離散為N個網格,然后按入射波是否直接照射將網格分為亮區和影區,最后對亮區內所有網格單元的散射聲場求和,即可得到目標散射聲場的近似值,由式(4)可得目標強度為

(5)

式中:Ii為第i個單元上的面積分。本文采取的網格類型為六節點曲面三角形單元,并建立了曲面單元上積分Ii的數值計算方法。

2.1 曲面三角形單元

考慮圖2(a)所示的參數三角形單元,該單元由參數坐標(s,t)中的6個節點(Ai)定義,則圖2(b)中曲面三角形上的點Q可表示為

圖2 參數三角形與六節點曲面三角形之間的映射

(6)

式中:Bi為曲邊三角形的頂點;Ri(s,t)為插值函數[20],簡寫為Ri,定義如下

(7)

曲面三角形單元上的積分Ii為

(8)

式中:Qs、Qt分別為點Q在s、t兩個參數方向上的偏導數,

(9)

(10)

Ri,s與Ri,t為插值函數的偏導數,Be為曲面三角形單元6個頂點的坐標,分別表示為

(11)

(12)

Be=(B1,B2,B3,B4,B5,B6)T

(13)

2.2 高斯-勒讓德積分

采用高斯-勒讓德求積方法,可將式(8)中三角形區域上的積分表示為

(14)

式中:(ξj,ηj)為積分節點;n為積分節點個數;wj為權重系數。一般地,可將積分區域從三角形[0, 1]×[0, 1-s]映射到正方形[-1, 1]×[-1, 1],作參數變換:

(15)

式中:ξ與η為參數坐標,參數變換對應的雅可比矩陣為

(16)

將式(8)中的積分映射到正方形區域上,表示為

(17)

3 數值算例及分析

首先對散射體表面進行幾何建模,然后將模型劃分為六節點曲邊三角形網格,獲得單元和節點信息后,采用數值積分對目標強度進行計算。曲邊三角形網格的最大邊長hmax[18]應滿足:

(18)

式中:c0為聲速;fmax為最大計算頻率。

算例中的水下目標表面為剛性,聲速c0=1 500 m/s,入射波為平面波,無特別強調時聲源與目標幾何中心的距離為10 000 m。

3.1 剛性球的散射

為了驗證曲面三角形單元數值積分方法的準確性,考慮半徑a=1 m的剛性球散射。計算頻段為50～8 000 Hz,步長50 Hz,數據點個數m=160,由式(18)可知,網格尺寸需滿足hmax≤62.5 mm。

由于剛性球形狀規則,可借助球Bessel函數jn和球Hankel函數hn,采用分離變量法獲得目標強度的Rayleigh簡正級數解,對式(3)中的Kirchhoff近似積分I直接求解,可得近似解析解。Rayleigh簡正級數解和Kirchhoff近似解析解分別為

(19)

(20)

取均方根誤差(L2范數):

(21)

首先將剛性球表面劃分為若干個六節點曲面三角形單元,對比不同積分方式所得的數值結果。網格尺寸hmax取60 mm與30 mm時,三角形區域上7積分點得到的L2范數為0.264、0.067,正方形區域上3×3積分點得到的L2范數為0.268、0.067。兩種積分方式所得結果相近,與正方形區域相比三角形區域的積分點數量更少,因此三角形區域7點高斯積分將用于接下來的數值算例及分析。圖3給出了目標強度TS隨頻率的變化曲線,可以看出,曲面三角形(hmax=60 mm)數值積分與Kirchhoff近似解析解吻合良好。與簡正級數解相比,在低頻區域(Rayleigh區)差距較大;當f>2 000 Hz時,三者之間的數值差距在1 dB以內,但Kirchhoff近似解與簡正級數解的起伏周期不一致。

圖3 球的目標強度

為驗證六節點曲面三角形數值積分方法的計算效率,將其與基于三節點平面三角形單元的板塊元方法對比,根據聲學網格經驗,平面三角形網格的尺寸不大于c0/(6fmax),即31.25 mm。提取曲面三角形單元的單元和頂點坐標等信息,生成平面三角形單元,此時二者的網格尺寸與網格數量相同,節點數量不同。然后分別采用數值積分法和板塊元方法計算剛性球的目標強度,表1給出了不同網格數量下兩種方法的計算時長。采用的處理器為Intel(R) Core(TM) i7-9750H,內存8 GB;MATLAB軟件版本為R2019b,為提高程序的執行效率,采用向量化的數組運算代替循環單元。

表1 計算球目標強度所需的網格及時長

由表1可知,當網格數量相同時,數值積分的計算速率大于板塊元方法。對于每個網格,數值積分中曲面三角形的積分點個數為7,大于板塊元中平面三角形的頂點個數3,但數值積分方法的時長較短,原因是板塊元方法需將全局坐標系下的頂點坐標變換到局部坐標系,涉及到龐大的矩陣運算,耗費時間較長。

圖4給出了不同網格尺寸下兩種方法的L2范數,為便于對比,hmax為30～60 mm時,板塊元的計算結果也于圖中給出。由圖4可知,目標強度的L2范數隨著網格細化而減小,當網格尺寸hmax≥25 mm時,曲面三角形數值積分方法,所得結果的誤差小于板塊元方法,具有較高的精度。原因在于與平面三角形相比,在網絡數量較少時,曲面三角形可以更好地擬合球面的形狀。雖然可以通過增加平面三角形的單元數量更準確地擬合曲面形狀,但劃分的網格數量很大。當網格尺寸hmax≤20 mm時,兩種算法得到的L2范數相近。

3.2 有限長橢圓柱的散射

考慮有限長剛性橢圓柱聲散射,取橢圓柱的中心為坐標原點,入射方向與z軸方向垂直,聲線處在Oxy平面內,如圖5所示,聲波照射亮區不包括端面,僅為側面。橢圓柱的尺寸參數為長半軸a=2.5 m,短半軸b=1 m,高度h=4 m。

圖5 有限長橢圓柱聲散射

采用穩相法積分導出的亮點模型[21]描述橢圓柱面的回波,可得

(22)

取入射角θ=45°,計算頻段2 000～10 000 Hz,步長50 Hz,數據點個數161,采用曲面三角形數值積分和平面三角形板塊元對其目標強度進行計算。表2給出了不同網格數量下的計算時長,可以看出,當網格數量相同時,數值積分方法的計算速率大于板塊元方法。

表2 計算橢圓柱目標強度所需的節點數量及時長

圖6給出了不用頻率下的目標強度曲線,計算時采用的網格尺寸為25 mm。然后計算頻率f=5 000 Hz、入射角θ為0°～180°時橢圓柱的目標強度,步長2°,得到的結果如圖7所示。由圖6與圖7可知,板塊元、數值積分和亮點模型3種計算方法的計算結果吻合良好。

圖6 不同頻率下橢圓柱的目標強度

圖7 不同入射角下橢圓柱的目標強度

3.3 球冠圓柱的散射

圖8為剛性球冠圓柱的散射特性,該模型由半球形冠部和有限長圓柱構成。幾何尺寸為球和圓柱的半徑r=0.25 m,圓柱高度h=1.75 m,坐標原點位于圓柱中心。θ=0°時聲波從半球面入射,θ=90°時聲波從正橫方向入射。

圖8 球冠圓柱聲散射

最大計算頻率為10 000 Hz時,網格尺寸需滿足hmax≤50 mm,取網格尺寸為50 mm,對應的網格數量為2 930,節點數量為5 926。圖9給出了θ=90°時球冠圓柱的目標強度,本文結果與有限元-完美匹配層法[2]得到的結果在高頻區域吻合良好,當f>2 000 Hz時,二者之間差距在1 dB左右;在中低頻較不準確,原因是本文方法基于Kirchhoff高頻近似,在低頻部分存在誤差。

圖9 正橫方向入射時球冠圓柱的目標強度

圖10給出了不同入射角度下球冠柱的目標強度,步長為1°。當θ=0°,只有半球面的聲散射,由式(20)可得入射波頻率為6 000 Hz、8 000 Hz時球冠圓柱的目標強度分別為-18.06 dB、-17.59 dB。當θ=180°,只有圓端面的聲散射,其目標強度為

圖10 不同入射角下球冠圓柱的目標強度

(23)

由此可得,頻率為6 000 Hz、8 000 Hz的平面波入射時,球冠圓柱的目標強度分別為-2.10 dB、0.40 dB,與當前計算結果吻合。

此外,由圖10可知,由于不同入射角度下,球冠圓柱的聲散射主導區域不同,目標強度隨入射角的變化較為復雜,大致可分為3個部分,主導區域分別為半球面、圓柱面和圓端面。正橫方向入射時球冠圓柱的目標強度最大,原因是正橫方向入射時,回波成分主要由圓柱面的鏡反射回波構成,散射的面積較大,所以目標的回波很強。

3.4 Benchmark潛艇的散射

最后研究Benchmark潛艇[22]的聲散射特性,按比例縮小的潛艇模型如圖11所示,縮小比例為 1∶15。 數值計算時,頻率取15 000 Hz與20 000 Hz,入射角θ為0°～360°,步長取1°,θ=0°時聲波從艏部入射,θ=180°時聲波從艉部入射。為滿足計算要求,曲面三角形網格尺寸取25 mm,共劃分單元 22 874個,節點45 750個;平面三角形網格尺寸取12.5 mm,共劃分單元91 412個,節點45 708個。

圖11 按比例縮小的Benchmark潛艇模型

為進一步驗證數值仿真結果的準確性,開展了湖上測試實驗,實驗模型的材料為不銹鋼,厚度 3 mm。 水聽器布置在模型與換能器之間,三者中心的布放深度均為10 m,模型距離水聽器6.55 m,換能器距離水聽器7.95 m,實驗布置圖如圖12所示。測試過程中保持換能器的位置固定,水平勻速旋轉處于懸掛狀態的縮比潛艇模型,以獲得全方位聲散射信號。

圖12 試驗設備布置圖

圖13給出了Benchmark潛艇縮比模型回聲強度ES的空間分布特性,由圖可知,數值積分結果與實驗測試結果的變化趨勢基本吻合,說明了本文方法的可靠性。當聲波從艉部及其附近入射時,實測值大于仿真值,原因是實驗模型包含仿真模型不具有的螺旋槳,螺旋槳對目標的回波具有貢獻。此外,在其他入射角度下,實測值與仿真值也存在一些偏差,一方面由于在實際測量過程中模型可能會有一定程度的傾斜,很難保證完全水平或標準正橫向入射;另一方面試件在加工及焊接過程中難以避免地會產生幾何尺寸及形狀誤差,從而影響實驗模型的目標強度。

圖13 Benchmark縮比模型的回聲強度

4 結論

本文基于Kirchhoff近似和曲面三角形單元建立了水下目標收發合置目標強度計算模型,采用六節點曲面三角形單元離散目標表面,并對亮區內曲面單元的散射聲場求和,得到目標的總散射聲場。對4種典型目標的聲散射計算結果表明:

1)與平面三角形相比,曲面三角形可以更好地擬合曲率較大目標的幾何形狀,所需的網格數量更少。

2)由于采用了Kirchhoff近似,本文方法更適用于高頻段內水下目標的散射特性分析。

3)通過湖上測試得到了Benchmark潛艇縮比模型目標強度的空間分布特性,測試結果與仿真結果吻合較好。

4)基于曲面三角形單元的高斯-勒讓德數值積分方法在高頻具有良好的求解精度,具有一定的工程實用價值。