二維透聲窗遮擋損失計算方法研究

賴喜悅，陳美霞，董文凱

（華中科技大學 船舶與海洋工程學院，湖北 武漢 430074）

0 引言

裝有聲吶設備的艦艇艏部透聲窗是水下反潛作戰的探測窗口[1]。透聲窗內部聲吶設備工作時，基陣發射的聲波經過導流罩透射到外部聲場，由于透聲窗的遮擋和導流罩內部結構的透聲損失，波束會發生畸變，并伴隨著能量損失。因此，透聲窗設計在滿足結構強度、阻力和耐用性等要求的前提下，還需要著重考慮其聲學特性[2]。

因透聲窗遮擋作用帶來的聲波能量損失被定義為遮擋損失，是透聲窗聲學特性設計關注的重點之一。然而，由于透聲窗結構比較復雜且模型尺寸一般較大，加之工程中關注的頻率范圍并不局限于低頻，透聲窗遮擋損失的準確計算較為困難。因此，在許多文獻中往往通過對平面波入射下單層或多層介質模型的透聲性能進行研究計算來粗略預報透聲窗的聲學性能。LEE 等[3]利用三層和四層介質模型對復合透聲窗材料在一定頻率下的插入損失進行了理論計算，并開展實驗驗證了計算結果的準確性。童暉等[4]利用三層介質模型對透聲窗材料在高頻下的聲學特性進行了研究，并探討了材料屬性和平面波入射角度對其插入損失的影響。易燕等[5]提出了一種低頻寬帶測試技術，解決了雙層鈦合金導流罩平板模型低頻段和大角度測試難度問題，得到了模型試樣低頻寬帶范圍內插入損失的頻率譜和角度譜。李東升等[6]利用夾芯透聲窗的簡化模型研究了夾芯透聲窗參數對其聲學性能的影響，并對比了不同頻段下夾芯透聲窗與單層玻璃鋼透聲窗透聲性能的差異。此外，李源源[7]等建造了實際艦艇球鼻艏導流罩模型，實驗得到了球鼻首的全向和定向發射透聲損失。

綜上，對透聲窗聲學性能進行準確預報是對其進行聲學設計的前提，而簡化的單層或多層介質模型無法考慮實際透聲窗的復雜結構、復雜邊界和復雜聲源，只能粗略地估算透聲窗的聲學性能。實驗方法研究透聲窗模型的聲學性能又會帶來較大的實驗成本，因此，提出一種準確預報透聲窗遮擋損失的計算方法具有切實意義。

本文提出了一種解析–數值混合方法，用于快速預報透聲窗的遮擋損失。首先，利用球諧波展開法得到了單活塞聲源的聲場分布，并通過疊加原理得到陣列聲源的聲場分布，進而根據透聲窗有限元模型各網格節點的坐標得到模型表面任意一點由聲源激勵引起的聲壓和法向振速。將之插值回COMSOL 有限元模型之中，調用 COMSOL 的Kirchhoff-Helmholtz 模塊通過再輻射的方法求解模型的遮擋損失，最后建立簡單二維透聲窗模型，驗證了方法的有效性。具體的計算流程如圖1所示。

圖1 解析–數值混合方法計算流程圖Fig.1 Flowchart of the calculation process of analysis-numerical hybrid method

1 活塞陣列聲源聲場分布

要想對透聲窗模型遮擋損失進行準確計算，聲源的建模和求解是第一步應該做的內容。圓柱形基陣是目前國內外艦殼聲吶最普遍的一種基陣形式，受水面艦吃水深度的限制，通常采用收發合置的聲基陣，使得這種基陣形式對水平 360°范圍具有相同的探測能力。實際研究中可以將這種基陣簡化為圓面活塞陣列，在定向發射時，可以根據相控規則調控各活塞振元的振動相位，使得各個活塞輻射聲場波陣面重疊，實現波束聚焦都效果[8]。

1.1 球諧波展開法求解單活塞聲場

文獻[9]和[10]采用球諧波展開法求解了圓盤形輻射體在聲場中的聲輻射，對于二維模型來說，陣列的每一個振元都可以簡化為以一定振速在自由空間中簡諧振動的線段，線段在y>0 方向的法向振速為V0，在y<0 方向的法向振速幅值為–V0，活塞半徑為a，r和θ分別為測點半徑及其與y軸正方向的夾角，模型示意圖如圖2所示。

圖2 單活塞聲場分布示意圖Fig.2 Schematic diagram of sound field distribution of single piston

將無限水域劃分為子域Ⅰ（0≤r≤a,0≤θ≤π /2）、Ⅱ（0≤r≤a,π /2≤θ≤π）、Ⅲ（r≥a,0≤θ≤π）后，各子域聲壓滿足時間簡諧形式如式（1）所示：

式中：ω為振動圓頻率；j為虛數單位；μ=1，2，3。

忽略時間項exp(-jωt)，各子域聲壓還應滿足如下邊界條件：Z0為介質特性阻抗；k為流體中的波數。

從而，各子域的聲壓函數可以假設為式（3）中的形式：

式中：An，Bn和Cn為未知系數；J和H分別表示一類貝塞爾函數和漢克爾函數，H的下標表示其階數，上標表示其為第一類漢克爾函數。

各子域之間需要滿足的聲壓–振速連續性條件如下：

將各子域假設的聲壓函數代入到上述連續性條件中，可以得到式（5）：

式（5）的第2、3 個方程左右兩邊同時乘以cos（mθ），m和n均為自然數，并在兩方程θ的取值區間內積分，可以得到式（6）：

其中分部積分可求得：

式中，pm（nm,n)=qmn(m,n)=,m=n，且m,n為奇數。

從而可知pmn和qmn有下面的性質，

取m=2s+1，s為自然數，將式（6）中的2 個方程相加，利用上述pmn和qmn的性質可以得到式（9）：

式（9）等式兩邊同時在r=a處對r求導，可以得到式（10）：

聯立式（9）和式（10），取截斷數s=n=N，即可以得到一個由2*（N+1）個方程組成的線性方程組，求解線性方程組，得到聲壓函數中的未知系數An，Bn和Cn，回代到式（3），即可得到單活塞輻射聲場。

1.2 疊加原理求解活塞陣列聲源聲場分布

多個活塞組成陣列時，陣列半徑為Rr，活塞個數為2*Ns+1，即以第Ns+1 個活塞為中心，左右兩邊各有Ns個活塞，活塞間隔圓心角為δθ，第i個活塞與中心活塞的夾角為φi，陣列發射方向與中心活塞的夾角為θr。為了實現定向發射的功能，不同活塞振動在時域上需要存在一定的時延關系，這種關系可以通過調節頻域上的相位差來實現，即不同活塞的速度振幅可以由如下形式表出：

式中，V0i表示第i個活塞的法向振動速度幅值。

忽略活塞之間的相互作用，陣列輻射聲場可以看作是單活塞聲場疊加而成，從而總聲場可以表示為

式中：TR，Tθ表示測點坐標；ri，θi表示測點在第i個活塞局部坐標系下的坐標。如圖3所示，對于平面上任意一點（TR，Tθ），我們都可以由余弦定理將其在第i個活塞局部坐標系下的坐標（ri，θi）表示成式（13）的形式：

圖3 局部坐標系下測點坐標示意圖Fig.3 Schematic diagram of measurement point coordinates in local coordinate system

獲取了陣列輻射聲場任意一點的聲壓P(TR,Tθ)后，結合透聲窗艙室幾何信息，就可以獲取透聲窗及艙室表面任意一點處的聲壓與法向振速。

2 透聲窗遮擋損失計算

以陣列中心為全局坐標系的原點建立透聲窗艙室輪廓模型，對于透聲窗艙室表面上任意一點，其坐標向量可以設為=[x,y]T，該點處的外法向單位矢量為=[nx,ny]T，從而該點處由陣列聲源激勵起的法向振速vn（）可以表示為式（14）的形式，式中dΔ 為一較小的值，可取對應頻率下流場聲波波長的1/20。

依據惠更斯原理，文獻[11]提出了一種再輻射的方法來計算透聲窗內聲吶基陣的遠場輻射指向性，這種再輻射的方法很好的適應了實際透聲窗內的聲場分布情況，運算比較簡單。結合這一方法，可以將遮擋損失的計算分為2 步：第一步是計算陣列在自由場中的遠場聲壓級指向性，第二步是計算透聲窗遮擋后遠場聲壓級指向性，兩者之差即為遮擋損失。

依據圖1，獲取表面振速信息后，調用有限元軟件COMSOL Multiphysics 6.0 的Pressure Acoustics，Kirchhoff-Helmholtz 模塊（壓力聲學，基爾霍夫–亥姆霍茲模塊）進行再輻射計算。通過輸入振速，該模塊基于基爾霍夫–亥姆霍茲積分公式對高頻輻射問題進行計算，而無需對周圍流體進行建模，可以極大地降低有限元模型的復雜度，使得計算較為迅速。調用Kirchhoff-Helmholtz 模塊后，在COMSOL中調用MATLAB 插值函數將振速輸入到幾何和網格表面，得到再輻射源。計算無遮擋聲場源特性時，需要將整個艙室外表面選中，給每個節點插值對應的法向振速，而計算遮擋后的輻射聲場時，只需要將透聲表面選中并賦予振速計算。最后，只需要在后處理中將2 次計算結果相減，即可得到模型的遮擋損失。

3 方法有效性驗證

3.1 單活塞聲場有限元驗證

在COMSOL 中建立二維單活塞模型，活塞半徑為35.5 mm，計算頻率為10 kHz，活塞法向速度幅值為1 m/s。有限元模型如圖4所示，其中內圓的水平直徑即為活塞，內圓和中間圓為水域，水的密度為1 000 kg/m3，水中聲速為1 500 m/s，外圓為PML 層。理論與有限元計算得到的該頻率下距離圓心r=1 000 m 處一周的聲壓如圖5所示。可以看到兩者的計算結果吻合較好，說明球諧波展開法能較準確地求解單活塞在水中的輻射聲場。

圖4 單活塞有限元模型Fig.4 Finite element model of single piston

圖5 單活塞r=1 000 m 指向性圖Fig.5 Directivity pattern of single piston at r=1 000 m

3.2 活塞陣列聲源聲場有限元驗證

在COMSOL 中建立活塞陣列模型，NS=10，各活塞的直徑均為35.5 mm，活塞法向速度幅值為1 m/s，陣列半徑為1.2 m，活塞間隔圓心角為3.5°，θr=0°，計算頻率為10 kHz，有限元模型示意圖如圖6所示，理論與解析計算得到的該頻率下距離圓心r=1 000 m 處一周的聲壓如圖7所示。也可以看到理論結果能與有限元結果吻合較好，陣列定向發射時，能量主要集中在發射方向。

圖6 活塞陣列有限元模型示意圖Fig.6 Schematic diagram of finite element model for piston array

圖7 活塞陣列r=1 000 m 指向性圖Fig.7 Directivity pattern of piston array at r=1 000 m

3.3 透聲窗遮擋損失計算有限元驗證

為驗證本文提出方法的準確性，下面將建立一個簡單的二維模型，分別采用本文方法與有限元法對模型的遮擋損失進行計算，并對結果進行對比驗證。

在計算時，取陣列發射角度為0°，通過改變陣列發射方向與透聲窗中軸的夾角θx，可以在掃描扇面內進行周向掃描，取遮擋前后陣列平面內沿著陣列發射方向距陣列中心r=1 000 m 處的場點聲壓級之差，為模型在θx角度時的遮擋損失，記為SL（θx），即

式中：LP(θx)為遮擋前的場點聲壓級，LP'(θx)為遮擋后的場點聲壓級。以θx為周向坐標，SL（θx）為徑向坐標，繪制極坐標曲線，就可以表征透聲窗全向的遮擋損失特性。

以一橢圓模型為驗證模型，橢圓長半徑為2.3 m，短半徑為1.8 m，橢圓中心在原點，截取長邊頂點兩邊各5/16 橢圓弧長作為透聲窗，其余弧為吸聲側壁驗證模型的陣列參數與3.2 節一致，全局坐標系原點與陣列中心重合。模型示意圖如圖8所示。

圖8 驗證模型示意圖Fig.8 Schematic diagram of verification model

分別采用純有限元和解析–數值混合方法得到的模型全向遮擋損失計算結果如圖9所示，可以看到，本文提出的方法準確性較高，且相較于有限元，該方法無需建立出完美匹配邊界和復雜的陣列活塞聲源，模型簡單，計算效率較高。

圖9 模型全向遮擋損失結果對比Fig.9 Comparison of results of model omnidirectional obstruction loss

4 結束語

本文提出了一種用于透聲窗遮擋損失計算的解析–數值混合方法，并通過二維簡單模型進行了方法驗證，得到了以下結論。

1）球諧波展開法可以有效求解二維活塞聲場，結果與有限元吻合較好，本文對各子域的聲壓函數假設合理。

2）不考慮活塞間的相互作用，運用疊加原理，可以較準確的計算出活塞陣列的聲場分布，計算結果與有限元吻合較好，陣列定向發射時，聲能量主要集中在發射方向上。

3）本文提出的解析–數值混合方法可以較準確的預報二維透聲窗模型的全向遮擋損失，思路清晰，計算簡單，相較于純有限元方法，不涉及無限域的模擬，也不需要建出活塞陣列（純有限元模型一般在活塞陣列附近網格質量要求較高），對硬件要求較低且耗時較短。

后續可以依據這一思路，繼續將該方法拓展至三維模型，方法快速性的優點將得到更充分的體現。