高晟耀，王德石，朱擁勇

（海軍工程大學 兵器工程系，武漢 430033）

聲納是利用回波聲場蘊含的目標特征信息，來實現水下目標的探測、識別和定位。為克服背景噪聲對其性能的影響，通常設計具有大空間指向性的換能器基陣并對工作帶寬加以限制，其形狀多采用曲面狀，如球冠形、拋物面形及橢球面形等。目前在結構輻射或散射聲場特性研究中，一方面限于工程結構往往含有不規(guī)則的幾何外形，結構邊界具有復雜的聲學參量分布，難以利用特殊函數、分離變量等方法獲得解析解；另一方面，與有限單元法相比，邊界積分方程法具有空間降維、自動滿足Helmholtz方程以及遠場輻射條件等優(yōu)點，這就使得該數值計算方法在聲學領域得到了普遍的重視。

而當邊界元法求解邊界積分單元上的未知量時，一旦輻射頻率等于特征頻率，就會面臨處理不定型的極限0 0問題，造成結構表面出現奇異積分，以及特征頻率處積分算法呈現振蕩特性，所以如何克服奇異性的固有缺陷，并提高計算效率和精度成為邊界元算法中的研究焦點。1963年，Chen和Schweikert[1]將邊界積分方程方法最早應用到聲學領域，在計算任意形狀結構聲輻射的過程中發(fā)現了一種與彈性波、電磁波研究中類似的現象，即對應Dirichlet內問題的特征頻率處的解并不存在。Schenck[2]在1968年指出該解是非唯一性而非不存在問題，并提出了一種改進的算法CHIEF法，此方法選取內域點的積分方程作為Helmholtz邊界積分方程的輔助約束條件。但在駐波節(jié)點分布密集的高頻段，CHIEF點愈加難以選擇。之后，Burton和Miller[3]提出了另一類具有代表性的改進算法，將奇異積分方程式和對其邊界法向求導后的超奇異積分方程式的純虛數倍進行線性組合，可以求得任意波數下的唯一解，但繁瑣的積分處理方式制約了數值計算的速度和效率。

為尋求傳統邊界元有效替代方法，Koopmann等人[4]運用了Treffetz法在求解邊值問題的思路，將描述聲場問題解表示為一系列滿足齊次偏微分控制方程的基函數的疊加。通過在振體結構的內部虛擬面上配置單偶極子源，來等效模擬結構外輻射聲場。由于虛擬曲面與實際振動邊界不重合，有效避免了外域數值積分中的奇異性問題，而且所建立的間接邊界元模型僅要求聲單元特征尺寸小于波長而非1 6波長，適用頻率范圍更廣。為了進一步解決虛擬曲面上相應的內問題的特征頻率處出現非惟一解問題，Jeans、Mathews和向宇、黃玉盈分別提出了復數形式Buton-Miller型組合層勢法[5]，以及在虛擬源強系統中增加虛擬阻尼的復數矢徑波疊加法[6]，并且討論了各自算法的計算效率問題。目前，該間接邊界元法已經被廣泛應用于近場聲全息、噪聲源識別、機械故障診斷等聲學逆問題之中，取得了大量機理性和相關實驗成果[7―9]，但涉及到遠聲場輻射特性的研究工作較少。

本文基于虛擬源強技術建立了結構聲輻射的間接邊界元模型，通過采用在內域虛擬曲面上配置一系列耦合極子源模擬結構外輻射聲場的方案，研究了一種球冠形輻射器的遠場聲輻射特性，并與傳統邊界元法的計算結果進行了比較分析，從而為分析設計大空間指向性的輻射器提供了一種新思路。

1 間接邊界元模型

考慮結構在理想流體介質中微小擾動形成的諧聲場，去除時間相關性e-iωt，外域問題由Helmholtz方程，Neumann邊界條件和Sommerfeld輻射條件描述

式中p(r)為場點r處復聲壓，k=ω/c為聲波數，ω為聲波的角頻率，vn為邊界面上法向振度。該偏微分控制方程的邊值問題難以給出解析解，常將其轉化為積分方程的形式，得到理論上能夠計算任意形狀振體輻射聲場的Kchhoff-Helmholtz積分方程

傳統邊界元方法通過離散邊界面上的積分方程，得到未知量{ps}和已知量{vn}的關系式

并利用離散化后的Helmholtz邊界積分方程，推導出外輻射聲場r處復聲壓為

式中N為離散單元數，矩陣A、B、C和D的元素見文獻[10]。

圖1 虛擬源強系統Fig.1 Virtual source intensity system

虛擬源強技術并非直接求解積分方程(4)，而是利用內域近離散單元形心處系列單偶極子源所輻射聲場的疊加，構造虛擬源強系統如圖1所示，從而模擬結構外輻射聲場

式中sm為虛擬源強，αm，βm為已知常數，取值1，0時即簡單源，表示無限大障板面上單元；取值0，ik時即偶極子源，表示僅彎曲變形薄殼單元；取值1，ik時即耦合極子源，表示封閉輻射體面單元。上述兩式表示成矩陣形式為

式中矩陣[H]和[T]為聯系等效源點和場點的傳遞函數，一旦配置虛擬源強系統后，傳遞函數可通過格林函數的計算得到。根據邊界面上單元體積速度

結合已知邊界法向振速求得虛擬源強{}s，之后代入(7)式和(8)式，可以得到輻射聲場，從而實現對結構聲輻射特性的研究。

2 算例與仿真分析

考慮到評價波疊加方法的計算精度，選取了具有解析解的球冠形輻射器作為研究對象，見圖2，其遠場點(r,θ,?)處輻射聲壓理論值為

式中為第二類m階球漢克爾函數，Pm為第m階勒讓德函數，θ0為球冠活塞面仰角，a為剛性球直徑，文中θ0=45°，a=1m。

圖2 球冠形輻射器Fig.2 Spherical cap radiator

以半徑r=18上的點為觀測坐標，為反映輻射聲場的指向性，以各頻率下的最大聲壓對聲壓幅值進行歸一化處理，得到極坐標下的理論聲壓幅值隨頻率變化關系如圖3所示。由于對于所有的輻射振膜來說，波束帶寬是非常重要的遠場聲輻射特性，往往采用半輻射聲壓即聲壓響應衰減-6d B時點間角度的間隔橫定，對該球冠形輻射器的半聲壓帶寬進行分析，理論結果如圖4所示。

根據直接邊界元法的面元劃分要求，即單元特征尺寸不大于λ6，離散球面為4 000個四邊形單元，網格劃分如圖5所示，若單元特征尺寸為λ3，四邊形網格的數目為1 000。

為比較兩種數值算法在計算遠時半聲壓帶寬的誤差，定義相對誤差公式

圖3 k a=4,12,20時遠場歸一化理論聲壓幅值Fig.3 Normalized magnitude of the theoretical pressure in far field atk a=4,10,20

圖4 r=18時理論半聲壓帶寬Fig.4 Variation of theoretical semi-bandwidth of sound pressure with frequency atr=18

圖5 球冠形輻射器表面網格劃分Fig.5 Surface mesh of the spherical cap radiator

式中An為數值計算半聲壓帶寬，Aa為理論半聲壓帶寬，n為球面網格劃分的單元數目。

在1≤k a≤20范圍內，經兩種數值算法求得了不同單元尺寸的波束帶寬相對誤差曲線，如圖7和圖9所示。由上述計算結果和表1可以得知：網格密度由λ6調整至λ3時，間接邊界元法的單個頻率計算時間將近是直接邊界元法的1 6，而且各自的計算效率分別提高了20倍和10倍；從計算精度角度看，直接邊界元法計算波束寬度的相對誤差水平始終維持在1%以內水平，間接邊界元法在k a=14處的最大相對誤差則由6.77%跳變到10.13%，但在多數頻率范圍內該方法的計算誤差并未超過4%，從圖6、8可以看出，所選擇的相對誤差較為敏感，間接邊界元法的數值計算結果能夠滿足工程應用所允許的誤差要求。

圖6 單元尺寸取λ 6時3種方法計算的半聲壓帶寬Fig.6 Semi-bandwidth of sound pressure with mesh densityλ 6

圖7 單元尺寸取λ 6時的帶寬相對誤差Fig.7 Relative error in the bandwidth with mesh densityλ 6

表1 k a=4時三種方法的半聲壓帶寬及其計算時間比較Tab.1 Acomparison between three approaches for analysis of bandwidth and computation time

3 結語

針對傳統邊界元法奇異性的固有缺陷，提出一種求解結構輻射聲場的間接邊界元法，應用該方法計算了球冠形輻射器的遠場聲輻射特性，將其結果與解析法以及直接邊界元法的結果進行比較。分析表明，該間接邊界元方法可在實波數域內獲取唯一解，有效克服直接邊界元法的不足。在保證工程應用所允許的誤差要求條件下，由于所需表面網格劃分的單元數目較少，其計算效率是直接邊界元法的3～6倍。同時，此方法也為分析設計具有大空間指向性的輻射器提供參考。

圖8 單元尺寸取λ 3時3種方法計算的半聲壓帶寬Fig.8 Semi-bandwidth of sound pressure with mesh densityλ 3

圖9 單元尺寸取λ 3時的帶寬相對誤差Fig.6 Relative error in the bandwidth with mesh densityλ 3

[1]Chen L H,Schweikert D G.Sound radiation from an arbitrary body[J].J.Acoust.Soc.Am,1963，35:1626-1632.

[2]Schenck H A.Improved integral formulation for acoustic radiation problems[J].J.Acoust.Soc.Am，1968，44(1):41-58.

[3]Burton A J,Miller G F.The application of integral equation methods to the numerical solution of some exterior boundary value problems[C].Proceedings of royal society of London,1971,A323:201-210.

[4]Koopmann G H,Song L and Fahnline J B,A method for computing acoustic fields based on the principle of wave superposition[J].J.Acoust.Soc.Am，1989，86(6)：2433-2438.

[5]Jeans R A，Matthews I C.The wave superposition method as a robust technique for computing acoustic fields[J].Journal of Acoustical Society of America,1992,92(2):1156-1166.

[6]向宇，黃玉盈.基于復數矢徑的波疊加法解聲輻射問題[J].固體力學學報，2004，25(1)：35-40.

[7]于飛，陳心昭，李衛(wèi)兵，等.空間聲場全息重建的波疊加方法研究[J].物理學報，2004，53(8)：2607-2613.

[8]李加慶，陳進，楊超，等.波疊加聲場重構精度的影響因素分析[J].物理學報，2008，57(7)：4258-4264.

[9]薛瑋飛，郭金泉，陳進，等.波疊加法在機械噪聲故障特征提取中的應用研究[J].機械強度，2007，29(6)：900-903.

[10]Wu T W.Boundary Element Acoustic:Fundamentals and Computer Codes[M].U K:WITPress,2000.