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

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

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

0 引言

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

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

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

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

1 Kirchhoff高頻近似

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

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

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

(5)

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

2.1 曲面三角形單元

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

圖2 參數(shù)三角形與六節(jié)點曲面三角形之間的映射

(6)

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

(7)

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

(8)

式中:Qs、Qt分別為點Q在s、t兩個參數(shù)方向上的偏導(dǎo)數(shù),

(9)

(10)

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

(11)

(12)

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

(13)

2.2 高斯-勒讓德積分

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

(14)

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

(15)

式中:ξ與η為參數(shù)坐標,參數(shù)變換對應(yīng)的雅可比矩陣為

(16)

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

(17)

3 數(shù)值算例及分析

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

(18)

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

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

3.1 剛性球的散射

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

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

(19)

(20)

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

(21)

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

圖3 球的目標強度

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

表1 計算球目標強度所需的網(wǎng)格及時長

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

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

3.2 有限長橢圓柱的散射

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

圖5 有限長橢圓柱聲散射

采用穩(wěn)相法積分導(dǎo)出的亮點模型[21]描述橢圓柱面的回波,可得

(22)

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

表2 計算橢圓柱目標強度所需的節(jié)點數(shù)量及時長

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

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

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

3.3 球冠圓柱的散射

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

圖8 球冠圓柱聲散射

最大計算頻率為10 000 Hz時,網(wǎng)格尺寸需滿足hmax≤50 mm,取網(wǎng)格尺寸為50 mm,對應(yīng)的網(wǎng)格數(shù)量為2 930,節(jié)點數(shù)量為5 926。圖9給出了θ=90°時球冠圓柱的目標強度,本文結(jié)果與有限元-完美匹配層法[2]得到的結(jié)果在高頻區(qū)域吻合良好,當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,與當前計算結(jié)果吻合。

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

3.4 Benchmark潛艇的散射

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

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

為進一步驗證數(shù)值仿真結(jié)果的準確性,開展了湖上測試實驗,實驗?zāi)Ｐ偷牟牧蠟椴讳P鋼,厚度 3 mm。 水聽器布置在模型與換能器之間,三者中心的布放深度均為10 m,模型距離水聽器6.55 m,換能器距離水聽器7.95 m,實驗布置圖如圖12所示。測試過程中保持換能器的位置固定,水平勻速旋轉(zhuǎn)處于懸掛狀態(tài)的縮比潛艇模型,以獲得全方位聲散射信號。

圖12 試驗設(shè)備布置圖

圖13給出了Benchmark潛艇縮比模型回聲強度ES的空間分布特性,由圖可知,數(shù)值積分結(jié)果與實驗測試結(jié)果的變化趨勢基本吻合,說明了本文方法的可靠性。當聲波從艉部及其附近入射時,實測值大于仿真值,原因是實驗?zāi)Ｐ桶抡婺Ｐ筒痪哂械穆菪龢?螺旋槳對目標的回波具有貢獻。此外,在其他入射角度下,實測值與仿真值也存在一些偏差,一方面由于在實際測量過程中模型可能會有一定程度的傾斜,很難保證完全水平或標準正橫向入射;另一方面試件在加工及焊接過程中難以避免地會產(chǎn)生幾何尺寸及形狀誤差,從而影響實驗?zāi)Ｐ偷哪繕藦姸取?/p>

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

4 結(jié)論

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

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

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

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

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