楊迎春，周其斗

(海軍工程大學(xué) 船舶與動力學(xué)院，武漢 430033)

邊界元法是研究聲輻射問題的主要數(shù)值方法之一。決定邊界元法求解精度的主要因素是邊界面上積分式的數(shù)值計算誤差。這類積分式積分核的特點通常是含有格林函數(shù)或其在某一方向上的方向?qū)?shù)。特點決定了場點與源點在同一單元時積分式將具有奇性。如何精確計算奇異積分式是邊界元法研究的重要課題之一。以往的文獻大都只討論靜止介質(zhì)條件下的計算，所用的方法主要有坐標變換法［1－3］，以及直接解析積分法［4］;Astley［5］與 Cowper［6］等也將高斯積分法用于計算。坐標變換法利用雅克比行列式與積分核分母相互約去因子來消除積分式的奇性，復(fù)雜的坐標變換會引入較多的舍入誤差和截斷誤差。直接解析積分法只適用于特定問題，實用性不大。利用高斯積分法求解具有二階弱奇性的積分式:Becker［7］認為對三維問題，積分核分母是源點與場點距離的二階無窮小量，無法直接用高斯積分法求出正確的積分值;而Astley認為可以將含有奇點的單元再次剖分為更小的單元，在每一單元上重新布置高斯求積節(jié)點可以獲得適當?shù)姆直媛?，此方法在網(wǎng)格密度較高時程序?qū)崿F(xiàn)具有很大難度。以上文獻都沒有給出介質(zhì)運動條件下處理二階弱奇性積分式的有效方法。

解決弱奇性積分式計算問題的關(guān)鍵就是要盡可能精確地求出奇點處的積分值。坐標變換和加密單元，都是為了減少奇點處的積分誤差。本文主要研究介質(zhì)運動條件下二階弱奇性積分式的計算方法，首次推導(dǎo)出了奇點處的積分解析值。本文方法較常規(guī)方法更加易于數(shù)值實現(xiàn)，而且具有很高的精度。論文雖然只以三角形單元為例，但所述方法稍加變換也適用于其它形狀單元，并可運用于光學(xué)等相似領(lǐng)域的研究。

1 二階弱奇性積分式計算

1.1 邊界積分方程中的弱奇性積分式

直角坐標系xyz中，假設(shè)介質(zhì)以流速M沿z軸正向流動，則聲傳播方程為:

式中的p為聲壓，本文假定聲波為簡諧波，時間因數(shù)是exp(itω)，其中波數(shù) k=ω/c，ω 是角頻率，c為聲速。式(1)共軛方程對應(yīng)的格林函數(shù)方程為:

上式中 Q:(x1，y12，z1)為源點，P:(x2，y2，z2)為場點，為運動介質(zhì)中的格林函數(shù)，式(2)對應(yīng)的格林函數(shù)解為:

式中:

類似于靜止介質(zhì)中Helmholtz面積分公式的推導(dǎo)方法［8］，可以得到運動介質(zhì)中的面積分公式:

C'(P)是P點位置的函數(shù)，對外域輻射問題而言，nz是反射物面單位法向量n在z軸方向的分量，并且n指離介質(zhì)。邊界元計算需要對積分面S進行離散化，式(6)的離散式中將含有以下幾類積分項:

三維空間中，場點與源點落在同一單元面SK上時，格林函數(shù)的分母R在場點附近趨于0，由于的分子有界，因此積分核將趨于無窮大，此時積分式在場點P處具有奇性。一般將式(7)的奇性定義為一階弱性;而式(8)、式(9)的奇性定義為二階弱奇性［9］。式(9)可看作是式(8)的一部分，可以用相同的方法處理。

常規(guī)數(shù)值方法，如高斯積分法，只能直接計算三維問題的一階弱奇性積分式，若計算二階弱奇性積分式則誤差較大。為此，論文推導(dǎo)了一種按照積分區(qū)域的奇性，將面積分式分兩部分積分的新方法:奇點部分由坐標變換計算出積分式的解析值;除奇點外單元面上的積分值，由高斯積分法求出。論文首先推導(dǎo)式(8)在奇點處的解析值。

1.2 二階弱奇性積分式在奇點處的解析值

以三角形單元為例，將Sk用ΔA1A2A3表示，并建立圖1所示的單元局部直角坐標系x'y'z'。

圖1 運動介質(zhì)中三角形單元上的半球積分面Fig.1 Integration hemispherical surface on triangle element in moving flow

局部坐標系原點o與三角形單元質(zhì)心(場點)P重合;z'軸與總體坐標系中的z軸平行，并且方向相同。易證ΔA1A2A3所在平面上總存在點d，直線od與z'軸垂直。以od作為x'軸，垂直于z'軸與x'軸的直線作為y'軸。由右手法則確定x'軸和y'軸正向。平面y'oz'和ΔA1A2A3的交線l與z'軸成 α角。

記Q為球面上的任意源點，將線段oQ與z'軸的夾角表示為φ，oQ在x'oy'上的投影長度表示為r，并記r與 x'軸正向所成角為 θ，則 α≤φ≤π +α，0≤θ≤π，球坐標系與局部坐標系的關(guān)系為:

Q處的單位面法向量為:

于是可將式(8)在球面∑上的積分寫為:

式(15)第二個等號后的負號是由源點與場點和面法線之間的相互作用關(guān)系確定的?！獹r與～Gz'分別是～G關(guān)于r和z'的偏導(dǎo)數(shù)。

當ε→0時E→1，此時式(15)方括號中只有第三項對結(jié)果有貢獻。做積分運算［10］:

上式即為式(8)在奇點P處的解析值，并且與α無關(guān)。當M→1時，式(17)取極限為2π;本方法不適用于介質(zhì)流速大于或等于1的情況。

1.3 高斯積分法

使用2.2節(jié)方法求出式(8)在奇點處的積分值后，單元Sk面上其他部分的面積分值可以用高斯積分法計算。使用高斯積分法的好處是，它可以用相對簡單的迭代算法計算各種形狀單元Sk上的面積分式。為了避免高斯積分法重復(fù)計算奇點處的積分值，推導(dǎo)了式(8)非奇性部分積分的變換式:

式中的r:

易于證明，式(18)中?(L/r)/?n在奇點P處的積分值與式(17)結(jié)果相同;而在其它部分的積分值為0。因此式(18)積分結(jié)果與式(8)非奇性部分的積分值相等。

論文以三角形單元的7點高斯法為例，說明如何布置高斯點以獲得足夠的精度，7點高斯法的詳細步驟可參閱文獻［10，11］。

將單元按質(zhì)心(場點P)和三個頂點連線剖分為三份，如圖2所示，在三個子單元上分別按照7點高斯積分法的規(guī)則布置求積節(jié)點。

由式(18)得到非奇性部分的高斯積分式為:

式中下標i表示相應(yīng)的子單元編號;下標k表示子單元上相應(yīng)高斯點的編號。Sk表示第i個子單元的面積，fik是式(18)的積分核在第i個子單元上的第k個高斯點處的計算值，wik是對應(yīng)此點的權(quán)系數(shù)。式(21)計算結(jié)果加上式(17)表示的解析值，即可得到式(8)的完整積分值。

圖2 三角形單元上高斯點布置方案Fig.2 The gaussian points location scheme on triangle element

對于式(7)，則將fik換為式(7)的積分核在高斯點上的計算值，式(21)計算結(jié)果即為式(7)的積分值。

對多個單元進行積分計算后發(fā)現(xiàn)，采用本文方法計算式(8)和式(9)要比Astley的結(jié)果具有更高的精度;而所用的子單元數(shù)目要比Astley的方法至少減小4個，求積節(jié)點密度至少減少2.25倍。

2 實驗驗證

為了檢驗方法，在流速為M的介質(zhì)中作一個球徑a=1 m的虛擬邊界球面。球心處有一點源，球面外的部分作為外域。按照“簡單源”理論［8］，點源將在球面上產(chǎn)生密度為σ的等效聲源。等效聲源在外域產(chǎn)生的輻射聲場與點源直接在外域產(chǎn)生的聲場，理想情況下應(yīng)該相等。具體方法如下所述:

點源強度［12］滿足N，對文獻［12］中的點源公式進行變換，得到運動介質(zhì)條件下點源在自由場中的聲壓解析式:

球面離散化后，由簡單源理論［13］，球面任意單元k上的聲壓導(dǎo)數(shù)滿足:

N為球面單元總數(shù)，nk是單元k上的法向量，與式(8)形式相同，σj是單元j上的源密度。場點P處的聲壓為:

圖3是球面邊界元模型，共有254個節(jié)點，504個三角形單元。

計算了不同流速條件下，距球心2.5 m處的聲場，0方向為流速方向，180方向為流速反向。流速為0，角頻率為1000時的聲壓指向性如圖4所示。

流速為0.5，角頻率為1000條件下的聲壓指向性如圖5所示。

圖3 球面網(wǎng)格模型Fig.3 Meshes of sphere

圖4 聲壓級圖(r=2.5 m，M=0，ω =1000)Fig.4 Sound pressure level(r=2.5 m，M=0，ω =1000)

圖5 聲壓級圖(r=2.5 m，M=0.5，ω = 1000)Fig.5 Sound pressure level(r=2.5 m，M=0.5，ω = 1000)

圖4、圖5中的實線是數(shù)值結(jié)果，三角形點是理論值。由圖可見，點聲源在流速影響下的聲壓級具有了指向性，0角附近的聲壓級增大，180附近的聲壓級減小。由于流動介質(zhì)中的計算要比靜止介質(zhì)中復(fù)雜，所以流動介質(zhì)中的數(shù)值計算結(jié)果與解析值比較的誤差相對較大。

3 結(jié)論

本文給出了三維聲學(xué)邊界元理論中，精確計算介質(zhì)流動條件下具有奇性積分式的方法;給出了二階弱奇性積分式在奇點處的解析值。數(shù)值實驗表明，此方法與理論值符合較好。新方法與常規(guī)方法相比，更加易于編程實現(xiàn)，而且精度較高。

本文雖然只介紹了新方法在平面三角形單元上的應(yīng)用，但本方法也可用于其它類型的單元。對于多邊形平面單元，可以將其劃分為多個三角形子單元分別進行積分運算;而對于曲面單元，則可借助型函數(shù)［14］先將曲面單元映射為平面單元，再用本文所述方法計算。

［1］Rizzo F J，Shippy D J.An advanced boundary integral equation method for three-dimensional thermoelectricity［J］.International Journal for Numerical Methods in Engineering，1977，11:1753－1768.

［2］Seybert A F，Soenarko B，Rizzo F J，et al.An advanced computational method for radiation and scattering of acoustic waves in three dimensions［J］.J.Aoust.Soc.Am，1985，77(2):362－368.

［3］Wu T W.Boundary Elementt Acoustics［M］.UK:WITpress，2000:51 －60.

［4］Chen JT， Chen K H. Dualintegralformulation for determining the acustic modes of a two-dimensional cavity with a degenerate boundary［J］.Engineering Analysis with Boundary Elements，1998，21(2):105 －116.

［5］Astley R J，Bain J G.A three-dimensional boundary element scheme for acoustic radiation in low mach number flows［J］.Journal of Sound and Vibration，1986，109(3):445 －465.

［6］Cowper G R.Gaussian quadrature formulas for triangles［J］.International Journal for Numerical Methods in Engineering，1973，7:405 －408.

［7］Becker A A.The boundary element method in engineering-A complete course［M］.UK:McGraw-Hill BOOK COMPANY，1992:112.

［8］Williams E G.Fourier acoustics:Sound radiation and nearfield acoustical holography［M］.San Diego:Academic Press，1999:Chapter 8.

［9］Mikhlin S G.Multi-dimensional singular integral equations［M］.Oxford:Pergamon Press，1965.

［10］數(shù)學(xué)手冊編寫組.數(shù)學(xué)手冊［M］.北京:高等教育出版社，1979:259，1053.

［11］章本照.流體力學(xué)數(shù)值方法［M］.北京:機械工業(yè)出版社，2003:196－198.

［12］張海瀾.理論聲學(xué)［M］.北京:高等教育出版社，2007:232.

［13］Zhou Q，Joseph P F.A numerical method for the calculation ofdynamic response and acoustic radiation from an underwater structure［J］.Journal of Sound and Vibration，2005，283:853－873.

［14］Zienkiewicz O C，Taylor R L.The finite element method，volume 1:the basis［M］.U.K:Butterworth-Heinemann，2000:Chapter 9.