吳紹維, 向 陽, 夏雪寶

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

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

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

圖1 波疊加法示意圖

(1)

(2)

表1 α,β值及其組合

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

(3)

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

(4)

(5)

定義Umn為

(6)

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

(7)

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

(8)

通過上述分析可知：用無單元空間離散域波疊加法進(jìn)行聲場計(jì)算，只需在實(shí)際輻射體表面選取離散的點(diǎn)，而不需要在輻射體內(nèi)部構(gòu)建虛擬聲源。

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

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

圖2 點(diǎn)聲源對小單元的輻射

2.1 單極子自輻射聲壓項(xiàng)近似解析表達(dá)式

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

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

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

(9)

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

(10)

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

(11)

對式(11)積分得：

(12)

當(dāng)δa0→0時，對gm取極限得到單極子自輻射聲壓項(xiàng)近似表達(dá)式如下：

(13)

2.2 單極子自輻射速度項(xiàng)近似解析表達(dá)式

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

(14)

vm=sm

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

圖4 單極子自輻射速度項(xiàng)積分域示意圖

vm=sm

(16)

因?yàn)樵趕0-s1上mR垂直于故采用球坐標(biāo)積分則

(17)

當(dāng)r→0時，對取極限得單極子極子自輻射速度項(xiàng)的近似解析表達(dá)式，即：

(18)

2.3 偶極子自輻射聲壓項(xiàng)近似解析表達(dá)式

(19)

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

(20)

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

圖5 偶極子自輻射聲壓項(xiàng)積分域示意圖

則：

(21)

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

(22)

在s1上采用球坐標(biāo)進(jìn)行積分，則

(23)

當(dāng)r→0時，對gd取極限得偶極子自輻射聲壓項(xiàng)的近似解析表達(dá)式，即：

(24)

2.4 偶極子自輻射速度項(xiàng)近似解析表達(dá)式

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

(25)

(26)

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

(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維空間中，點(diǎn)到坐標(biāo)原點(diǎn)的距離。需要計(jì)算的積分為

(40)

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

(41)

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

(42)

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

(43)

3 自輻射近似表達(dá)的實(shí)例驗(yàn)證

圖6 無單元脈動球源模型

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

(44)

(45)

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

(46)

(47)

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

(48)

(49)

由方程(4)很容易得到

(50)

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

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

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

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

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

4 結(jié) 論

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

參 考 文 獻(xiàn)

[1]張雄, 劉巖. 無網(wǎng)格法[M]. 北京市：清華大學(xué)出版社，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ìn)，楊超,等. 基于波束形成和波疊加法的復(fù)合聲全息技術(shù)[J]. 聲學(xué)學(xué)報(bào)，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ì)時，吳崇健，曾革委,等. 基于波疊加法的圓柱殼聲輻射計(jì)算[J]. 艦船科學(xué)技術(shù)，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]程鴻翔，尚德江，李琪,等. 聲場匹配波疊加法的水下結(jié)構(gòu)聲輻射預(yù)報(bào)[J]. 聲學(xué)學(xué)報(bào)，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]杜功煥，朱哲民，龔秀芬. 聲學(xué)基礎(chǔ)[M]. 第三版.南京市：南京大學(xué)出版社，2012.