三維淺海下圓柱殼聲輻射預報方法研究

錢治文， 何元安， 商德江， 肖 妍， 史 博

（1. 天津大學 海洋科學與技術學院， 天津300072； 2. 哈爾濱工程大學 水聲工程學院，哈爾濱150001； 3. 中國船舶工業系統工程研究院， 北京100036）

0 引 言

我國海洋領土大多屬于水深小于200 m 的典型淺海環境， 潛器在淺海所處的聲場環境并非理想的自由場，存在水面和水底的反射、散射作用[1]，其聲場特性與結構本身和淺海環境兩者密切相關。 研究淺海下結構輻射聲場特性，對開展淺海環境下結構振動輻射噪聲實時監測預報、聲學隱身性能評價以及淺海環境參數對輻射聲場影響分析等研究領域具有重要的理論和應用價值， 是今后我國水聲技術領域長期關注的熱點和難點問題之一。

然而，目前關于水下彈性結構耦合振動與聲輻射問題的研究，大多考慮無界或半空間流體域下的近距離輻射聲場問題[2-3]。 對于淺海環境下結構輻射聲場問題，國內外很少涉及該方面的研究，因該研究在試驗法、解析法和數值法研究過程中均遇到了難題，有效開展工作較為困難；試驗法周期長成本高，數值法受網格劃分以及計算量的嚴重制約，解析解法目前只能針對理想淺海環境下簡單結構聲場問題[4-5]。為了解決淺海下彈性結構聲輻射研究所遇到的這些難題，本文引入了一種對結構和流體環境適應性很強且計算高效的研究方法:波疊加法（Wave Superposition Method，簡稱WSM）。Koop-mann 等[6]最先提出了利用波疊加法進行結構輻射聲場的計算，Miller 和Fahnline 等隨之對穩定性、 計算精度和效率進行了分析[7-8]。波疊加法以適應性強、計算效率高的優勢逐步地成為了研究彈性結構聲輻射的重要方法之一，近年來受到了國內學者極大的關注，并廣泛應用于自由場或半空間流體環境下的近場聲全息[9-11]、結構聲輻射預報[12-13]和噪聲源識別[14-15]等研究領域。 但采用波疊加法進行淺海下結構聲輻射的研究卻很少，因為淺海下結構輻射聲場問題涉及流體、結構和淺海上下邊界等多個物理場之間的相互耦合作用，導致結構表面振動速度難以準確獲取，從而無法進行波疊加法的準確計算。 也有部分研究把淺海下結構輻射源視為點聲源，但這樣直接忽略了結構與流體、結構與環境的耦合作用和近場聲輻射特性。 導致目前尚無高效可靠的研究方法進行淺海環境下彈性結構聲輻射預報的有效研究，但它對淺海下結構聲輻射預報、源識別和噪聲控制等具有重要的意義，急需探索一種新方法來解決淺海下結構聲輻射預報問題[16]。

本文提出了一種三維淺海波導下彈性結構輻射聲場預報的理論方法。 該方法首先利用有限元法建立了淺海波導下結構聲輻射多物理場耦合有限元模型； 計算獲取在考慮流體環境和淺海波導上下邊界耦合影響下結構表面振動數據，然后采用波疊加法和淺海聲場傳輸函數進行虛擬源源強求解，最后進行彈性結構輻射聲場的等效疊加計算。 經有限元法和解析解法驗證該方法的正確性后，研究了典型彈性圓柱殼在理想淺海波導下的輻射聲場特性， 為今后開展復雜淺海環境下大型彈性結構聲輻射預報研究提供了一種新途徑。

1 理論模型

1.1 淺海波導下點源聲傳播模型

考慮聲源為水下簡諧點源，設海水為理想介質，壓力p 滿足非齊次亥姆霍茲方程

式中:▽2為拉普拉斯算符，c 為流體聲速，r→為聲源到場點的矢量距離，δ r→（）為狄拉克函數，A 為點源強度幅值，j 為復數虛部，ω 為角頻率。

由于聲場具有圓柱對稱性，選用柱坐標下，以通過點源且垂直向下的方向為z 軸，垂直于z 軸的方向為r 軸，方程可以寫成[1]

式中:k0為流體介質中聲波波數，定義為k0=ω/c，z0為聲源在z 軸上的位置。

采用簡正波理論求解式，可得淺海波導下點源聲場表達式p （r, z ）即淺海波導傳輸函數G （r, z ）為

海面通常為Dirichlet 邊界條件，滿足的邊界條件為

海底是一個反射和散射邊界， 典型海底邊界模型有Neumann 邊界、Rayleigh 定律、Sommerfeld 模型、Cauchy 邊界和地聲模型等。

在Neumann 邊界、Rayleigh 定律和Sommerfeld 模型的交互面上，滿足聲學邊界方程為

1.2 波疊加法模型

建立如圖1 所示的三維淺海下結構聲輻射波疊加法計算理論模型，流體密度為ρ1，聲速為c1；結構中心和點源在海面的投影為o 和o′，o-xy 和o′-x′y′與海面重合，o-xyz 為全局坐標系，o′-x′y′z′為各虛擬源的局部坐標系，oo′連線延長線為y′軸，o′y′與oy 的夾角為θ； 振動體表面記為S′，n 為結構表面的外法向矢量，zs為結構中心在豎直方向的位置；S 為與結構面共形的等效源面，Q 為虛擬源，z0i為各虛擬源在局部坐標中的位置；P 為場點位置，Q 和P 間的距離記為rp。

圖1 三維淺海波導下波疊加法原理示意圖Fig.1 Principle of wave superposition method in three-dimensional shallow water waveguide

波疊加法與Helmholtz 積分公式等效，經過離散化處理后，全局坐標系下結構表面某一點的法向速度un可由N 個簡單聲源構建[18]

把（7）式寫成矩陣形式，可表示為

其中，［D］-1為［D ］的廣義逆矩陣。

為了求解傳遞矩陣的廣義逆［D］-1，需要對其進行奇異值分解

式中:S、V 分別為M 階、N 階方陣，Σ=diag （σ1, σ2, σ3…σn）為包含n 個非奇異值σi的對角陣。

通過奇異值分解求逆矩陣獲取輻射體結構內部虛擬源源強后，便可計算任意一點的聲壓

從（9）式和（11）式可知，傳遞矩陣D 和單極矩陣T 中涉及的傳遞函數G （x, y,z ）均為全局坐標下三維聲場傳遞函數，而1.1 節求解的聲場傳遞函數為局部柱對稱坐標下二維聲場傳遞函數G （r, z ），需要進行如下空間坐標轉化。

淺海下各虛擬源產生的空間聲場滿足軸對稱性， 三維局部坐標系o′-x′y′z′下聲場可降為二維軸對稱局部坐標系o′-rz′下聲場，根據局部二維坐標下各點源到場點水平距離r 和豎直距離z′求解各點源對場點的聲場作用G （r, z′ ）即G （x′ , y′, z′ ）（其中x′=rcos?，y′=rsin?，? 為r 軸與x 軸夾角），然后通過局部坐標系與全局坐標系的轉換關系建立各點源在全局坐標系下的聲場作用G （x, y,z ）。首先將局部坐標系聲場函數G （x′ , y′, z′ ）沿x 方向、y 方向和z 方向分別平移Δx、Δy 和Δz， 得到全局坐標系下聲場G （xi, yi, zi），平移變換矩陣為

然后再將全局坐標系下聲場函數G （xi, yi, zi）繞z 軸順時針旋轉角度θ， 得到與全局坐標系o-xyz下的聲場函數G （x, y,z ），旋轉變換矩陣為

通過局部坐標o′-rz′到全局坐標o-xyz 的轉換， 可求得各局部坐標下點源聲場函數G （r, z′ ）對全局坐標系下結構表面的作用G （x, y,z ）。

由（9）式和（11）式可知，該方法計算的基本思路為:首先通過有限元法建立淺海波導下彈性結構聲輻射數值模型，計算獲取結構表面法向速度U，然后通過波疊加法以及傳輸函數G 進行虛擬源源強Q 求解，最后結合單極矩陣T 便可計算在任意場點的輻射聲場P。

2 數值分析

2.1 方法準確性與高效性分析

如圖2 所示，建立了淺海波導下脈動球有限元模型，波導深度h=50 m，海面為Dirichlet 邊界，海底為Neumann 邊界，流體四周無限大邊界。 小球中心深度為zs=25 m，半徑為r0=3 m，表面施加均勻振速ua，數值計算獲取表面振動數據，然后采用波疊加法計算場點聲場信息。

圖2 淺海波導下脈動球聲輻射有限元模型Fig.2 The finite element model of spherical shell in shallow water waveguide

如圖3 所示， 利用本文淺海波導下輻射聲場波疊加計算式即（11） 式計算了淺海波導下脈動球輻射聲壓級隨距離的變化曲線（各場點深度為h2=30 m），并分別與解析解法即（14）式、有限元法計算結果進行了對比分析，驗證了該方法計算聲場的可靠性。

圖3 波疊加法計算結果驗證Fig.3 The verification of wave superposition method

圖4 FFP 與Kraken 計算結果對比Fig.4 The comparison between FFP and Kraken

雖然采用簡正波法存在一定的近場聲場計算精度問題， 但從遠場波疊加法聲場計算結果驗證來看，其計算精度是滿足要求的。為了說明采用簡正波法的聲場波疊加計算精度問題，本節采用最為精準的波數積分法計算程序FFP 與簡正波計算程序Kraken 進行了計算結果對比分析。 計算環境以Pekeris波導為例，海水密度ρw=1 024 kg/m3，聲速cw=1 500 m/s；半無限液態海底密度ρw=1 800 kg/m3，海底聲速cw=2 500 m/s；海水深度h=200 m，聲源深度zs=100 m，接收點深度zr=50 m，如圖4 所示，分別采用FFP 和Kraken 計算了點源在淺海波導內的傳播損失曲線對比圖。

從圖4 可知，采用簡正波法與波數積分法計算結果在近場存在一定偏差，且波數積分法計算結果更為準確，當達到一定距離后，兩者的計算結果基本吻合。 但因簡正波的推導簡單， 低頻計算效率高，所以本文采用簡正波進行波疊加法中聲傳輸函數的求解，且聲場波疊加法計算精度也是滿足要求的。 為了進一步完善淺海下結構聲輻射預報理論體系，提高波疊加法計算精度，后續將通過從波數積分法改進淺海波導聲傳輸函數、優化結構表面測點和虛擬源布放方式進行深入研究。

表1 不同距離范圍聲場計算時間Tab.1 The calculation time of different distance range

為了顯示該方法在進行淺海波導下彈性結構輻射聲場計算方面的高效性， 以下采用該方法對不同計算距離范圍進行了時間測試，如表1 所示。

從圖3、圖4 和表1 可看出，該方法計算結果不僅與解析解法、有限元法計算結果吻合很好，而且計算效率高，能夠快速準確進行淺海波導下任意結構輻射聲場預報。

2.2 典型淺海波導下圓柱殼聲輻射分析

建立如圖5 所示的理想淺海下彈性圓柱殼聲輻射多物理場耦合有限元模型，淺海深度h=50 m，海面為Dirichlet 邊界，海底為Neumann 邊界。結構為彈性圓柱殼，半徑a=3 m，長度l=30 m，厚度d=0.01 m，材料為鋼（密度ρs=7 850 kg/m3，楊氏模量Es=2.05×1011Pa，泊松比us=0.28），圓柱殼中心距水面為25 m，在圓柱殼中間施加徑向簡諧力Fr=1 000 N，結構邊界為自由邊界，h2為圓柱殼軸線方向上各個場點距離水面的距離。

圖5 淺海波導下結構聲輻射模型Fig.5 Acoustic radiation model of structure in shallow water waveguide

首先采用有限元法對淺海下彈性圓柱殼聲輻射進行建模計算，獲取了圓柱殼表面附近192 個測點的振動數據。根據文獻[18]虛擬源布放建議，在離散的結構內部空間均勻布放162 個虛擬源，虛擬源源面為與彈性結構共形的圓柱面，其中虛擬圓柱面長為24 m，半徑為2.4 m。 通過（9）式計算虛擬源源強后，采用（11）式計算了各個虛擬源在任意場點的疊加聲場。 等效計算彈性結構輻射聲場后，便可研究淺海下彈性圓柱殼輻射聲場特性。

2.2.1 遠距離聲輻射曲線

圖6 不同頻率的輻射曲線Fig.6 Radiation curves at different frequencies

采用圖5 所示的有限元數值模型計算獲取圓柱殼結構在典型淺海波導下的圓柱殼表面振動數據，再通過波疊加法和淺海波導傳輸函數計算淺海下任意場點的聲場信息。 如圖6 所示，分別給出了淺海波導下彈性圓柱殼不同頻率所對應的輻射聲壓級隨水平距離變化曲線， 各場點連線與圓柱殼軸線平行，各個場點的深度取為h2=30 m。在結構近場附近，輻射聲場受到結構尺寸影響較大，輻射聲壓級曲線變化較為復雜。 由于場點聲場受到淺海波導上下界面的干涉疊加影響， 圓柱殼在整個距離上的輻射聲場發生了相長和相消的疊加，聲壓級曲線波動明顯。 但當場點距源中心達到一定距離后，輻射聲場受結構尺寸影響較小，其輻射聲壓級衰減規律按柱面波即滿足-10lg（r）衰減。 因整個聲場計算是由多個虛擬源產生聲場疊加而來的，所以總聲場的波動規律相對于單個點源聲場波動細節較多。 在近場，場點與各個虛擬源之間距離的大小相差較大，導致在采用（11）式進行聲場計算時，所涉及的單極矩陣T 相差也較大，即進行了不同相位和幅度的聲場疊加計算，導致結構近場輻射聲場波動比較復雜。 在遠場，由于場點距虛擬源的距離遠大于各個虛擬源之間的距離，各個虛擬源與場點之間的距離相差較小，所以各虛擬源在場點產生聲場近似為同相疊加，疊加的總聲場波動規律較為穩定，聲壓級按-10lg（r）遞減規律進行衰減。

2.2.2 空間聲場分布

圖7 空間場點選取示意圖Fig.7 The diagram of space field selection

圖8 不同流體環境下輻射聲壓對比Fig.8 The comparison of sound pressures in different fluid environments

為顯示圓柱殼在垂直方向上的近場聲場分布特性，按圖7（a）所示的場點選取方式（場點選擇在距離結構中心h2=25 m 的圓周上，圓周所在平面垂直圓柱殼軸線，極角90°和270°處分別對應淺海海面和海底），利用波疊加法計算了50 Hz、100 Hz、200 Hz 和300 Hz 頻率下，圓柱殼輻射聲壓在不同流體環境中的空間聲場分布，并分別與相同條件下圓柱殼在自由場、半空間流體環境中計算結果進行比對分析，如圖8 所示。通過對比三種流體環境下輻射聲場空間分布可看出， 邊界的存在會使在靠近邊界處的聲場分布出現較大波動，淺海波導下結構輻射聲場由于同時受淺海波導上下界面的影響，在深度方向上聲場為駐波場，各個駐波干涉疊加導致在整個圓周上均出現明顯的旁瓣。 在不同頻率時，下界面和上界面對整個圓柱殼輻射聲場分布影響作用不同，如當頻率為200 Hz 時，在上半圓周上的聲場分布已經受到了淺海下邊界反射聲的嚴重影響，而不會出現在其他頻率如50 Hz、100 Hz 和300 Hz 時，淺海和半空間環境下輻射聲場在上半圓周的聲場分布情況趨于一致。 相對于較高頻率，根據簡正波理論[1]，50 Hz頻率下在深度方向上只有3 階簡正波且對應波長較長， 聲場干涉不明顯， 輻射聲壓空間分布波動較小，在豎直方向上的空間分布曲線類似為圓形。 隨著頻率的增加，在深度方向上的簡正波數目和聲場波動細節增加，輻射聲壓空間分布曲線的旁瓣數量增加。 而由于淺海下場點聲場由不同成分的反射聲與直達聲干涉疊加構成，結構在淺海下疊加的總聲壓幅值均大于自由場以及半空間下的聲壓幅值。

圖9 為在水平面上選取了不同深度場點下的輻射聲壓級（單位:dB）空間分布，場點選取如圖7（b）所示，場點在水平距圓柱殼中心hs=1 000 m 處的水平圓周上，各場點距離水面距離為h2，極角0°和90°處分別對應圓柱殼端面與側面。 可看出，在水平方向上，界面反射聲對圓周上各個場點的影響作用相同，所以聲場分布特性主要與結構本身的輻射聲場相關。 在低頻段時，結構輻射聲場指向性分布不明顯，且結構在水平方向的輻射聲場為擴散聲場，無駐波的干涉作用。 頻率為30 Hz 時，圓柱殼在不同場點深度下的空間聲場分布均為圓形，說明在30 Hz 處的輻射聲場可等效為單極點源聲場，而其他頻率下聲場空間分布特性在不同場點深度下會出現明顯的分布差異，結構源不可等效為點源。

2.2.3 簡正波模式分析

淺海波導下圓柱殼輻射聲場受到了淺海上下界面的限制，在深度方向存在簡正波。 為了研究駐波中簡正波模式對圓柱殼聲輻射的影響， 圖10 分別計算了60 Hz 頻率下只保留1 個、2 個、3 個和4 個簡正模式（（3）式中n=1, 2, 3, 4）的圓柱殼聲輻射曲線。 由（11）式可知，圓柱殼在任意場點的輻射聲場是由多模式疊加作用的結果，隨著簡正波模式的增加，圓柱殼聲輻射曲線的細節也在增加。

圖10 不同簡正波模式的輻射曲線Fig.10 Radiation curves at different normal modes

為了分析各模式輻射聲場特性，可以把淺海下聲場傳輸函數（（3）式）表示為簡化的聲強形式

當采用（11）式進行結構聲輻射聲場計算時，在傳輸函數G（r, z）中只保留一個模式，復聲壓包含振蕩項ejξrr，包絡曲線是平滑的；保留兩個模式時，聲強包含了cos ［（ξr1-ξr2）r ］，聲輻射曲線出現了波動，而且干涉曲線的重復尺度遠大于該頻率下的波長；當G（r, z）保留3 個以上的輻射模式，聲強表達式包含的級數項進一步增多，干涉結構的復雜性加大。

3 結 論

本文采用WSM 進行了三維理想淺海下彈性圓柱殼輻射聲場特性的研究，主要結論如下:

（1） 在近距離處，淺海波導下圓柱殼各頻率所對應的聲輻射曲線受結構尺寸影響較大，聲場分布較為復雜；在遠距離處，由于結構尺寸對聲場影響減小，輻射聲場波動規律穩定，按柱面波擴展規律波動，即輻射聲壓級隨水平距離r 的變化按-10lg（r）進行擴展。

（2） 在豎直方向空間聲場分布上，由于聲場受上下界面影響，輻射聲場為駐波場。 淺海波導下的聲場分布不同于自由場以及半空間下聲場分布，其聲場分布差異隨著頻率的增加而增大；在水平方向的空間聲場分布上，輻射聲場為擴散場。 且在低頻處（如本文分析頻段30 Hz 處），遠場圓柱殼空間聲場分布可等效為無指向性點源聲場，而其他頻率下輻射聲場出現明顯的指向性分布，不可等效為點源聲場，必須以圓柱殼作為實際輻射源進行聲場計算，以獲取其真實的聲場分布特征。

（3） 結構內部某一虛擬源在場點的聲場是由各個模式簡正波疊加計算獲取的，且不同簡正波模式下的疊加數目不同，導致最終結構輻射聲場波疊加計算的干涉曲線波動規律不同，隨著傳輸函數中簡正模式的增加，輻射曲線波動細節增多，波疊加法的計算結果更為精確。