呂國策

長春職業(yè)技術(shù)學(xué)院，吉林長春 130000

基于邊界元法聲場輻射預(yù)測研究

呂國策

長春職業(yè)技術(shù)學(xué)院，吉林長春 130000

采用邊界元法預(yù)測結(jié)構(gòu)振動的聲場聲壓，已經(jīng)得到廣泛研究，然而邊界積分方程建立的過程中遇到奇異積分，雖然此時也采用高斯積分，但計算收斂性降低。本文針對奇異積分問題采取廣義極坐標(biāo)變換，對比了兩種結(jié)構(gòu)網(wǎng)格對邊界元法計算結(jié)果的影響，結(jié)果表明，采用廣義極坐標(biāo)變換后，邊界元法可以預(yù)測振動輻射聲場聲壓，當(dāng)頻率升高時，某些頻率結(jié)果誤差增加，為此，在提高網(wǎng)格質(zhì)量后，在低頻情況兩種結(jié)果與解析解極為接近，高頻段較好的網(wǎng)格質(zhì)量度數(shù)值計算結(jié)果的精度有所提高，但計算時間乘倍增加，不利于工程實際，對于大型問題要綜合考慮，在不影響計算精度的同時，可以適當(dāng)降低網(wǎng)格數(shù)量，以降低計算時間。

邊界元；奇異積分；廣義極坐標(biāo)變換

由于機(jī)械振動產(chǎn)生強(qiáng)烈的噪聲，工程中振動產(chǎn)生的噪聲越來越受到人們的關(guān)注，而噪聲對現(xiàn)場工作人員產(chǎn)生極大心理和勝利危害，為解決此問題，各國都相應(yīng)出臺了相應(yīng)的措施限制現(xiàn)場噪聲，機(jī)械制造商相應(yīng)的也對其產(chǎn)生的噪聲采取了相應(yīng)的控制措施。對生產(chǎn)制造商而言，在進(jìn)行制造前提前預(yù)知機(jī)械產(chǎn)生的噪聲強(qiáng)度對后續(xù)工作有指導(dǎo)意義，所以需要一種有效地計算輻射聲場聲壓方法[1-3]。

1 邊界元方程的建立

聲在介質(zhì)中傳播，假設(shè)介質(zhì)為理想流體，流體不存在粘滯性，沒有擾動，宏觀上靜止，且介質(zhì)絕熱，并且聲波是小振幅聲波，聲壓遠(yuǎn)小于介質(zhì)靜態(tài)聲壓，由以上假設(shè)可以將非線性問題轉(zhuǎn)化為線性聲輻射問題，而我們關(guān)注的多數(shù)問題都是脈動振源產(chǎn)生的穩(wěn)定的簡諧聲波，所以有聲壓輻射公式[4-5]

根據(jù)理想流體小振幅聲波輻射的波動方程，可以得到helmholtz等式。

其中k為波數(shù)， /k cω= ，其伴隨方程為

伴隨方程的基本解為

我們僅考慮輻射聲壓，所以有B=0，有伴隨方程的基本解

利用第二格林公式建立邊界積分方程

假設(shè)基本解和聲壓滿足此方程，得到場點聲壓與邊界信息之間的關(guān)系

2 邊界奇異積分的處理

對于給定的結(jié)構(gòu)，可能僅了解邊界表面振動部分信息，然而通過以上邊界積分方程可知，要得到場點聲壓，必須要有所有邊界信息，包括表面振動速度和聲壓，當(dāng)給定邊界條件時，可以將場點p配置到邊界表面上，然而基本解這時存在奇異性，對于奇異積分的處理，本文采用極坐標(biāo)變換法引進(jìn)雅克比消除奇異積分的方法。

如圖1所示，將三角主單元轉(zhuǎn)化為正方形單位單元，由此的到轉(zhuǎn)化關(guān)系[6]

由上式帶入圖坐標(biāo)點可得

變換后的得到的雅克比

圖1 廣義極坐標(biāo)變換

3 脈動求輻射聲場驗證

脈動球源是指隨時間在進(jìn)行均勻振動的球體，表面做簡諧振動，脈動求向它存在的聲場輻射聲壓，滿足本文假設(shè)條件，輻射聲場聲壓可以通過解析公式計算解析解，作為驗證邊界元法的一種有效手段，與邊界元法計算得到的數(shù)值解進(jìn)行對比，驗證邊界元法的有效性。

對脈動球半徑r0=0.5(m)，介質(zhì)密度ρ0=1.21kg/m3，脈動求表面振動速度ua=1m/s建兩種單元模型，脈動球聲場輻射聲壓解析解計算公式有

對脈動球邊界離散采用四邊形網(wǎng)格，選用兩種網(wǎng)格尺寸，得到場點網(wǎng)格，第一種為四邊形網(wǎng)格單元，單元尺寸為50mm，節(jié)點數(shù)為1269，單元數(shù)為1267，第二種為單元尺寸30mm，節(jié)點數(shù)和網(wǎng)格數(shù)分別為3164、3162，對比網(wǎng)格質(zhì)量提高后對計算聲場聲壓精度的影響。

由于Fortran語言計算速度較快，適合大型程序開發(fā)與維護(hù)，針對大型計算問題，可以將模型分配幾種情況單獨計算，這樣可以縮減計算時間，本文邊界元積分方程程序采Fortran語言編寫。

對于脈動球場聲壓本文分別計算不同場點和不同頻率下聲壓進(jìn)行對比，表1為場點網(wǎng)格1對應(yīng)的不同頻率下場點聲壓解析解和數(shù)值解對比結(jié)果，從表可以的得出，在低頻段，解析解與數(shù)值解吻合較好，當(dāng)頻率升高后，數(shù)值解相應(yīng)誤差增加。

表2為場點網(wǎng)格2對應(yīng)解析解與數(shù)值解對別結(jié)果，與表3變化趨勢類似，隨頻率增加后聲壓幅值誤差有所增加，對比表1，在頻率700Hz以上，計算效果明顯改善，分析原因由于脈動球建模劃分網(wǎng)格質(zhì)量原因造成表4結(jié)果更近解析解。

以上結(jié)果分別對比了不同網(wǎng)格不同頻率對比結(jié)果，表3為網(wǎng)格1計算的到不同場點聲壓對比結(jié)果，頻率為50Hz，不同場點得到的聲壓在低頻情況邊界元法計算得到的數(shù)值結(jié)果與解析結(jié)果誤差較小，并且場點與脈動球源距離變化后誤差并沒有明顯變

表4為在1500Hz頻率下計算場點聲壓結(jié)果對比，當(dāng)頻率升高后，場點聲壓解析解與數(shù)值解仍相差較小，誤差在1%以下，滿足工程應(yīng)用需要。

對于兩種網(wǎng)格，計算時間也有所區(qū)別，網(wǎng)格1節(jié)點較少，單元較少，而網(wǎng)格2節(jié)點數(shù)和單元數(shù)為網(wǎng)格1的2.5倍，兩種網(wǎng)格計算時間分別為277.479 （s）、1568.574（s），網(wǎng)格2的計算時間是1的5.65倍，嚴(yán)重增加了計算時間。

4 結(jié)論

場點聲壓隨頻率變化數(shù)值解與解析解吻合較好，當(dāng)頻率增加聲壓略微產(chǎn)生波動。場點與脈動球距離變化不大時，聲壓誤差沒有明顯變化。

網(wǎng)格質(zhì)量提高后對聲場聲壓預(yù)測低頻精度沒有明顯提高，高頻誤差減小，但計算時間數(shù)倍增長。

存在誤差可能原因是在離散邊界是認(rèn)為邊界單元為平面單元，造成誤差。在計算過程中考慮到實際問題不可能將邊界離散很多單元，所以本文離散邊界單元較大，也造成誤差。但這不影響對實際問題的分析。

表1 網(wǎng)格1不同頻率對比結(jié)果

表2 網(wǎng)格2不同頻率對比結(jié)果

表3 網(wǎng)格1場點聲壓結(jié)果

表4 網(wǎng)格2場點聲壓結(jié)果

[1]高貴兵，張道兵，岳文輝，等.一種具有唯一解的邊界元法計算三維聲散射問題[J].機(jī)械科學(xué)與技術(shù)，2013，32（8）：1172-1175.

[2]郭偉，蔡明，馬俊.聲輻射邊界元數(shù)值計算[J].船海工程2011，40（1）：101-103.

[3]黎勝，趙德有.用邊界元法計算結(jié)構(gòu)振動輻射聲場[J]，大連理工大學(xué)學(xué)報，2000，40（4）：391-394.

[4]杜功煥，朱哲民，龔秀芬.聲學(xué)基礎(chǔ)[M].南京：南京大學(xué)出版社.

[5]Ursell F. On the exterior problems of acoustics[C]//Mathematical Proceedings of the Cambridge Philosophical Society. Cambridge University Press， 1973， 74(01): 117-125.

[6] Boundary element acoustics: Fundamentals and computer codes[M]. Wit Pr/Computational Mechanics，2000.

O29

A

1674-6708（2015）145-0113-02