基于無單元聲波疊加的自輻射近似解析表達研究

吳紹維, 向 陽, 夏雪寶

(1.武漢理工大學 能源與動力工程學院，武漢 430063；2.船舶動力系統(tǒng)運用技術交通行業(yè)重點實驗室，武漢 430063)

隨著有限元法近十幾年來的巨大發(fā)展，有限元法已成為工程數(shù)值分析的有效工具，解決了很多有重大意義的科學和工程問題。但是，有限元法在分析高速沖撞、動態(tài)裂紋擴展和應變局部變化等問題時遇到了因網格畸變產生的許多困難[1]，自上世紀90年代結構力學領域逐步發(fā)展了無單元計算方法[2]，這種無單元方法起源于裂紋擴展問題[3]。裂紋擴展問題的計算要求很密集的網格單元，并且在分析計算中還要求網格自適應，這增加了計算負擔。在無單元計算方法中去除了網格只考慮計算模型中的節(jié)點，每個節(jié)點給予一定的用于計算的區(qū)域，這樣很大程度上減輕了計算負擔。Koopmann[4]提出了傳統(tǒng)的波疊加法，即任何振動輻射體表面輻射的聲場可通過該輻射體內部虛擬構建的簡單球型源產生的聲場疊加進行求解，利用實際輻射體表面法向速度邊界條件確定虛擬聲源的強度。這些虛擬構建的聲源滿足波動方程，并且這種方法被證明等效于Helmholtz積分方程。這種方法對結構簡單的模型很實用，與邊界元法相比具有計算速度快，精確性高，當虛擬聲源位于輻射體內時無奇異的特點。但是虛擬聲源必須位于實際輻射體內的一定范圍內才能用于聲場計算。對于脈動球源模型，虛擬球型聲源半徑與實際球源半徑之比須在0.05—0.4[4]范圍以內才能較為準確的計算出表面聲壓；對結構復雜的模型，虛擬等效源的位置和所在的幾何形狀對聲場計算的精確性具有很大的影響，目前還處于進一步研究當中[5-7]。受到無單元計算方法的啟示，Koopmann 的學生Brian Zellers將無單元計算方法用到計算結構振動聲輻射的問題中，將無單元計算方法與波疊加法相結合使得聲場的計算大為簡化，對求解任意振動結構輻射的聲場只需矩陣運算[8-10]，不需要構建位于實際輻射體內部的虛擬聲源。但在Brian Zellers的研究中未能推導出偶極子自輻射項的解析表達，他采用文獻[11]中無障板活塞輻射阻抗的表達式來代表偶極子自輻射速度項，這種表達方式無法直接計算，需通過多個方程的復雜求解過程來獲取計算所需要的參數(shù)。針對所存在的問題，本文對無單元空間離散域的聲波疊加法進行了進一步的研究，分別通過奇異點挖去法推導單極子自輻射聲壓項近似解析表達，積分區(qū)域替換法推導單極子自輻射速度項和偶極子自輻射聲壓項的近似解析表達，不變量嵌入法推導偶極子自輻射速度項近似解析表達。最后用脈動球源作為實例，驗證了單、偶極子點聲源自輻射項的精確性。

1 無單元空間離散域的聲波疊加法

圖1 波疊加法示意圖

(1)

(2)

表1 α,β值及其組合

本文使用體積速度邊界條件來確定聲源強度sn[12]。體積速度定義為振動結構的表面法向速度在振動表面面積區(qū)域上的積分。為保證計算的精確性，這里將振動輻射表面劃分成N個區(qū)域，區(qū)域上的體積速度為法向速度在域上的積分。由歐拉方程可知聲場中任意一質點速度與聲壓的關系為：

(3)

這里ρ為聲介質密度(1.21 kg/m2)，其中m是對接收點求梯度。這里速度為方程(3)可以寫成

(4)

(5)

定義Umn為

(6)

用矩陣形式表示方程(6)為

(7)

通過體積速度條件解出聲源強度[sn]，并代回到方程(1)即可確定輻射體表面聲壓。已知聲壓和速度，將平均聲強在一個遠場處包圍輻射體的S面上積分，可求得平均聲功率，其中式(8)中的上標H表示復共軛轉置，即Hermitian轉置

(8)

通過上述分析可知：用無單元空間離散域波疊加法進行聲場計算，只需在實際輻射體表面選取離散的點，而不需要在輻射體內部構建虛擬聲源。

2 格林函數(shù)奇異性及自輻射項的近似解析表達

采用無單元法進行聲場波疊加計算時先將振動結構表面離散化，離散化的區(qū)域幾何中心稱之為節(jié)點，將點聲源置于這些節(jié)點處。本文在推導自輻射項的近似解析表達式時，取以節(jié)點為圓心，a0為半徑的圓形作為積分區(qū)域，圓形的面積等于離散化時節(jié)點所在的區(qū)域面積，如圖2所示。將位于節(jié)點處的點聲源對自身所在的圓形區(qū)域所輻射的聲壓及由輻射導致的速度在圓形區(qū)域上積分并取平均，平均聲壓和平均速度可用于表達單、偶極子自輻射的聲壓項和速度項，在本文中對單極子點聲源和偶極子點聲源的自輻射項進行了研究，接下來將給出詳細的推導過程。

圖2 點聲源對小單元的輻射

2.1 單極子自輻射聲壓項近似解析表達式

為了解決單極子模型中格林函數(shù)的弱奇異性，去掉以節(jié)點為圓心，半徑為δa0的圓形部分，如圖3所示。將位于節(jié)點處的單極子聲源對環(huán)形區(qū)域所輻射的聲壓及由輻射導致的速度在環(huán)形區(qū)域上積分并取平均，這些值將隨著δa0變化，當δa0趨近于0，平均聲壓和平均速度可用于表達單極子自輻射的聲壓項和速度項。

圖3 單極子對環(huán)形小單元的輻射

根據聲學理論和波疊加法原理，單極子輻射的聲壓可表達為第一類自由空間格林函數(shù)與聲源強度的乘積,即

(9)

這里sm為單極子聲源強度，r為單極子點聲源與環(huán)形區(qū)域上點之間的的距離。單極子對環(huán)形區(qū)域輻射的平均聲壓定義為

(10)

這里a0為圖3所示積分區(qū)域半徑，δa0為去掉的那一部分的半徑。令

(11)

對式(11)積分得：

(12)

當δa0→0時，對gm取極限得到單極子自輻射聲壓項近似表達式如下：

(13)

2.2 單極子自輻射速度項近似解析表達式

根據方程(3)，圓形區(qū)域上接收點處的速度v與聲壓p之間的關系為

(14)

vm=sm

為克服被積函數(shù)的奇異性，將積分區(qū)域分為兩部分，一部分替換為半徑為r的半球面域s1，另一部分為s0-s1的環(huán)形平面域，如圖4所示，其中r為單極子聲源與s1上的點之間的距離，R為單極子聲源與s0-s1上的點之間的距離。則:

圖4 單極子自輻射速度項積分域示意圖

vm=sm

(16)

因為在s0-s1上mR垂直于故采用球坐標積分則

(17)

當r→0時，對取極限得單極子極子自輻射速度項的近似解析表達式，即：

(18)

2.3 偶極子自輻射聲壓項近似解析表達式

(19)

偶極子對圓形區(qū)域輻射的平均聲壓定義為:

(20)

其中a0為圓形區(qū)域的半徑，s0表示圓形積分區(qū)域。為克服被積函數(shù)的奇異性，將積分區(qū)域分為兩部分，一部分替換為半徑為r的半球面域s1，另一部分為s0-s1的環(huán)形平面域，如圖5所示。

圖5 偶極子自輻射聲壓項積分域示意圖

則：

(21)

其中r為偶極子聲源與s1上的點的距離，R為偶極子聲源與s0-s1上的點的距離。在環(huán)形平面域s0-s1上，dR垂直于故則

(22)

在s1上采用球坐標進行積分，則

(23)

當r→0時，對gd取極限得偶極子自輻射聲壓項的近似解析表達式，即：

(24)

2.4 偶極子自輻射速度項近似解析表達式

根據方程(3)，圓形區(qū)域上接收點處的速度v與聲壓p的關系為

(25)

(26)

為克服被積函數(shù)中格林函數(shù)的超奇異性，對(26)中的被積函數(shù)進行恒等變形

(27)

(28)

令

(29)

(30)

則G=G1+G2。將G1展開為泰勒級數(shù)

(31)

將式(30),式(31)代入式(28)得到:

(32)

(33)

令

(34)

(35)

(36)

很容易求得

I2=πa0k2

(37)

(38)

對于I1這樣的超奇異積分，采用Invariant Imbedding Method[13]求取有限積分值。定義積分函數(shù)

(39)

其中R為n維空間中，點到坐標原點的距離。需要計算的積分為

(40)

并定義β=lj1+lj2+…+ljn-h為被積函數(shù)f的度，α=β+n為該奇異積分I的度。則對于α≠0的奇異積分的積分有限值為

(41)

對于I1,β=-3，α=-1則

(42)

?S0是圓形區(qū)域的邊界。用fpI1替代I1，將式(37)，式(38)，式(42)代入式(33)，從而偶極子自輻射速度項的近似解析表達式為：

(43)

3 自輻射近似表達的實例驗證

圖6 無單元脈動球源模型

脈動球源所輻射的聲壓解析解為[14]

(44)

(45)

當rs>a時沒有奇異性問題發(fā)生，考慮到脈動球源在整個表面均勻的輻射，每個點聲源強度均相等，用ss表示，可求得ss如下：

(46)

(47)

這里g11代表自輻射聲壓項，其中ps為式(44)中當rs=a時所確定的球面上的聲壓。解方程(47)得：

(48)

(49)

由方程(4)很容易得到

(50)

圖7 20個節(jié)點模型單極子自輻射聲壓項近似解析表達式與數(shù)值解對比

圖10 60個節(jié)點模型單極子自輻射速度項近似解析表達式與數(shù)值解對比

圖13 60個節(jié)點模型偶極子自輻射聲壓項近似解析表達式與數(shù)值解對比

圖14 60個節(jié)點模型偶極子自輻射速度項近似解析表達式與數(shù)值解對比

由以上結果對比中可以看出，20個節(jié)點單、偶極子模型自輻射聲壓項與速度項近似解析表達式與數(shù)值解在ka=0～4內具有較好的一致性，聲壓項近似解析表達式在接近ka=4時有輕微的偏離，這是因為當頻率增加時，聲波波長變小，用于計算的圓形區(qū)域尺寸已經不滿足所規(guī)定的條件(要求圓形區(qū)域的尺寸須小于λ/6，這類比于有限元計算中求解頻率越高，網格必須劃分的越細，以滿足小單元尺寸與聲波波長的關系)。對于60個節(jié)點，單、偶極子模型近似解析表達式與數(shù)值解也具有較好的一致性，當ka≥7時，自輻射聲壓項近似解析表達式與數(shù)值解的偏離相比ka取較小值時增大，這也是因為隨著k值增加，波長減小，用于積分的圓形區(qū)域尺寸不滿足所規(guī)定的條件所致。

4 結 論

本文針對無單元聲波疊加自輻射項存在的弱奇異、奇異和超奇異問題，利用部分積分區(qū)域替換等方法推導了單、偶極子的自輻射聲壓項和速度項，通過實例驗證表明所求得的自輻射項近似解析表達式在中低頻時能夠較為準確的表達自輻射項，可用于聲場的計算，即克服了奇異性問題。但在中高頻時，自輻射項與數(shù)值解存在一定程度的偏離。在較高頻率處，是因為波長變小用于積分的圓形區(qū)域尺寸不滿足條件所致，直接的方法是取更多的聲源點，但是增加聲源點數(shù)會導致相鄰聲源距離減小，其格林函數(shù)出現(xiàn)近奇異性的問題。因此解決這種矛盾還需要進一步研究，此問題解決后即可解決高頻偏離問題。

參 考 文 獻

[1]張雄, 劉巖. 無網格法[M]. 北京市：清華大學出版社，2004.

[2]Nayroles B, Touzot G,Villon P. Generalizing the finite element method: diffuse approximation and diffuse element[J]. Computational Mechanics, 1992, 10(5): 307-318.

[3]Belytschko T, Lu Y Y, Gu L. Element-free galerkin methods[J]. International Journal for Numerical Methods in Engineering, 1994,37(2):229-256.

[4]Koopmann G H, Songl,Fahnline J B, A method for computing acoustic fields based on the principle of wave superposition[J].Journal of Acoustic Society America，1989，86(6)：2433-2438,

[5]李加慶，陳進，楊超,等. 基于波束形成和波疊加法的復合聲全息技術[J]. 聲學學報，2008，33(2):152-158.

LI Jia-qing, CHEN Jin, YANG Chao. A hybrid acoustic holography technique based on beamforming and ave superposition algorithm[J]. Acta Acustica, 2008, 33(2):152-158.

[6]熊濟時，吳崇健，曾革委,等. 基于波疊加法的圓柱殼聲輻射計算[J]. 艦船科學技術，2011,33(1): 54-58.

XIONG Ji-shi, WU Chong-jian, ZENG Ge-wei. Sound radiation numeration of cylinder based on the wave superpositionmethod[J]. Ship Science and Technology, 2011,33(1): 54-58.

[7]程鴻翔，尚德江，李琪,等. 聲場匹配波疊加法的水下結構聲輻射預報[J]. 聲學學報，2013,38(2)：137-146.

CHENG Hong-xiang, SHANG De-jiang, LI Qi. Sound radiation prediction for underwater structure by field-matching wave superposition method[J]. Acta Acustica, 2013,38(2)：137-146.

[8]Zellers B. Element free structural acoustic for efficient shape optimization[C]. IMECE2005-82958.Orlando, Florida USA, 2005,125-135.

[9]Zellers B. An acoustic superposition method for computing structural radiation in spatially digitized domains[C].IMECE2007-41313.Seattle,Washington,USA,ASME,2007.

[10]Zellers B. An acoustic superposition method for computing structural radiation in spatially digitized domains[D]. Pennsylvania: Department of Mechanical and Nuclear Engineering, The Pennsylvania State University,2006.

[11]Mellow T, Karkkainen L. On the sound field of an oscillating disk in a finite open and closed circular baffle[J]. Journal of the Acoustical Society of America, 2005, 118(3), 1311-1325.

[12]Koopmann G H. Designing quiet structures[M]. 525 B Street, Suite 1900,San Diego, California 92101-4495,USA: Academic Press ,Inc 1997.

[13]Vijayakumar S. An invariant imbedding methoad for singular integral evaluation on finite domains[J]. SIAM Journal on Applied Mathematics, 1988,48(6)： 1331-1349.

[14]杜功煥，朱哲民，龔秀芬. 聲學基礎[M]. 第三版.南京市：南京大學出版社，2012.