蘇楠，龐福振，姚熊亮

結(jié)構(gòu)外聲場的映射變階波包絡(luò)無限元法

蘇楠，龐福振，姚熊亮

（哈爾濱工程大學(xué)船舶工程學(xué)院，哈爾濱150001）

基于徑向形函數(shù)可任意變階的映射波包絡(luò)聲學(xué)無限元法，以無限長圓柱殼體為研究對象，用數(shù)值計算方法對其結(jié)構(gòu)外聲場進(jìn)行研究分析，并與其解析解進(jìn)行對比，結(jié)果表明：兩者之間能夠較好地吻合，從而驗證了映射變階波包絡(luò)聲學(xué)無限元法在計算結(jié)構(gòu)聲輻射中的可行性，且具有效率高、精度好等優(yōu)點。在此基礎(chǔ)上，本文還討論虛擬極點位置對聲輻射的影響，通過對比四極子極點偏心聲壓值發(fā)現(xiàn)：極點存在偏心時，對實際工程問題有不利影響，但隨聲學(xué)無限單元階數(shù)的增加，誤差會減小。

映射變階波包絡(luò)法；聲學(xué)無限元；虛擬極點；聲輻射

1 引言

計算結(jié)構(gòu)外聲場問題可歸結(jié)為無界域聲波波動方程的數(shù)值求解問題，如水下航行器的聲輻射問題、機(jī)翼的流激噪聲等。有限元法在求解此類問題時，計算量大且效率不高，并會帶來較大的聲學(xué)截斷誤差。因此，為了克服有限元法在求解無限域結(jié)構(gòu)外聲場問題時的不足，從而提出無限元法。

關(guān)于無限元法的相關(guān)問題研究國內(nèi)外已有大量文獻(xiàn)發(fā)表，其中具有代表性有：Bettess和Zienkiewicz[1-2]在研究表面水波和折射問題時首次提出了無限元的概念，將聲學(xué)問題在無窮遠(yuǎn)處滿足的解直接移植到人工邊界上；并針對無限元的不足提出了一種映射無限單元，通過單元建立了局部坐標(biāo)系與整體坐標(biāo)系的映射關(guān)系，但在求解過程中碰到了從未定義過的指數(shù)型積分，在數(shù)學(xué)上沒有嚴(yán)格的理論推導(dǎo)證明；Astley[3]依據(jù)Bettess[4]無限元思想提出了一種新的無限單元類型，這一單元類型中取權(quán)函數(shù)為形函數(shù)的復(fù)數(shù)共軛，運用Garlerkin加權(quán)殘值法求解，該方法徹底消除了單元矩陣中的不確定積分，極大地簡化了方程，而且系統(tǒng)矩陣中頻率與剛度矩陣相互獨立，從而可使用高斯積分對方程進(jìn)行求解，后來此法被Cremers等人[5-7]將此聲學(xué)無限元法發(fā)展成為徑向方向可任意變階的映射波絡(luò)無限元法，并將幾何映射函數(shù)和形函數(shù)分開處理。楊瑞梁[8-12]對聲學(xué)無限元法的發(fā)展和研究進(jìn)行了系統(tǒng)的闡述，同時對無限元法進(jìn)行了詳細(xì)地分類，從理論上分析了收斂性、條件數(shù)和角向量數(shù)對計算精度的影響。

在上述研究的基礎(chǔ)上，本文基于徑向方向任意變階的映射波包絡(luò)聲學(xué)無限元法，較為詳細(xì)地推導(dǎo)了系統(tǒng)矩陣平衡方程，并與解析法進(jìn)行對比，驗證了映射變階波包絡(luò)聲學(xué)無限元法求解結(jié)構(gòu)外聲場的可行性，在此基礎(chǔ)上，著重分析了虛擬極點偏心對聲輻射計算精度的影響。

2 映射變階波包絡(luò)聲學(xué)無限元法

本文基于徑向形函數(shù)可任意變階的聲學(xué)無限元法，形函數(shù)沿兩個方向相互獨立構(gòu)造：在徑向上，采用任意階數(shù)的拉格朗日多項式，從而實現(xiàn)了徑向形函數(shù)的任意變階；在橫向上，保持形函數(shù)不變，這樣可實現(xiàn)與低階有限元的直接耦合。

圖1 聲場幾何模型Fig.1 The geometry model of sound field

2.1 域內(nèi)控制方程

考慮內(nèi)邊界為S1的無限大區(qū)域V，如圖1所示。p（ x，t）為邊界S1上的法向聲壓，假設(shè)該聲波為一簡諧波，p（ x，t）可表示為p（ x，t）=p（x） eiwt，它在聲場域V內(nèi)滿足以下波動方程：

式中：V為外部無界區(qū)域，k=w/c為聲學(xué)波數(shù)，▽2為拉普拉斯算子。在結(jié)構(gòu)輻射表面S1應(yīng)滿足以下邊界條件：

式中：U（）

x為法向振速。在無窮遠(yuǎn)處應(yīng)滿足Sommerfeld條件：

無限大外域聲場問題可歸結(jié)為求解邊值問題（1）-（4）。

2.2 方程離散

為求解邊值問題（1）-（4），引入較為簡單的人工邊界S2，將域V劃分為內(nèi)部區(qū)域V1和外部區(qū)域V2。在區(qū)域V1內(nèi)采用常規(guī)的有限元法離散；在區(qū)域V2內(nèi)采用無限單元離散，無限單元為向無限遠(yuǎn)延伸的四邊形，如圖2所示。此無限單元由一正方形母單元映射而來，節(jié)點1和2位于人工邊界S2上，節(jié)點

式中：r為整體坐標(biāo)系中坐標(biāo)原點到區(qū)域V內(nèi)一點的距離。

如果在遠(yuǎn)場邊界S2上，S2距離振動體足夠遠(yuǎn)，等效近似有限聲學(xué)邊界條件為：1、2、3、4分別對應(yīng)母單元中的四個節(jié)點，其他兩節(jié)點（1，1）和（1，-1）映射到無限單元的無窮遠(yuǎn)處。節(jié)點1′、2′與節(jié)點3、4關(guān)于節(jié)點1、2鏡像。從局部坐標(biāo)系ζ-ξ到整體坐標(biāo)系x-y的映射函數(shù)為：

圖2 映射波包絡(luò)單元Fig.2 Mapped wave envelope element

式中：

沿著邊1-3和邊2-4的幾何映射關(guān)系式可以表示為：

式中：r為沿1-3或2-4邊距聲源點1′或2′測得的距離。

2.3 單元系統(tǒng)矩陣

記Wi為權(quán)函數(shù)，Ni為形函數(shù)，對Helmholtz方程（1）應(yīng)用Glerkin加權(quán)殘值法，將在V上的體積分轉(zhuǎn)化為在S1和S2上的面積分，得到等效的Helmholtz方程為：

同時，聲學(xué)載荷量可以寫成：

那么，系數(shù)qi表達(dá)式可以通過方程式｛q｝=（［K］+ik［R］-k2［M ］）-1｛F｝求出，其中：｛q｝為節(jié)點聲壓列向量，［K］為剛度矩陣，［R］為阻尼矩陣，［M］為質(zhì)量矩陣，｛F｝為載荷向量列矩陣，k為聲學(xué)波數(shù)。

形函數(shù)及權(quán)函數(shù)選取方法具體參考文獻(xiàn)[6]。對于圖2中的二維單元，各點的聲壓是基于幾何節(jié)點1和2的聲壓值線性插值求得。波包絡(luò)無限單元的插值函數(shù)Ni，其中包括模擬幅值衰減部分和聲場中聲壓相位變化部分。對于n階波包絡(luò)法中僅考慮各自形函數(shù)下n個節(jié)點的聲壓，并認(rèn)為第（n+1）個節(jié)點的聲壓值在無窮遠(yuǎn)處為零。

母單元中n階多項式形函數(shù)采用拉格朗日多項式求得，其中n階拉格朗日多項式由分列在幾何節(jié)點之間的n個聲學(xué)節(jié)點確定。n階徑向形函數(shù)可表示為：

通過使用局部映射函數(shù)ζ=ζ1+τh=τ/（n-1）-1，其中h=1/（n-1）。當(dāng)ζ1=-1時，徑向形函數(shù)可寫成以下形式：

在母單元中使用此徑向形函數(shù)映射到實際單元中為1/r的展開式，它適用于模擬三維或軸對稱聲波傳播中的幅值衰減問題。但在二維聲學(xué)問題中，聲波幅值近似以形式衰減。對于二維問題，它的徑向形函數(shù)由乘以系數(shù)因子相當(dāng)于在母單元中乘以相應(yīng)系數(shù)項。此時，二維徑向形函數(shù)可以表示為：

在形函數(shù)中引入形如e-ikr的周期分量描述波形相位變化。為保證聲學(xué)單元的連續(xù)性，在有限單元與無限單元的交界面ζ=-1處相位必須設(shè)定為零。因此，在局部坐標(biāo)系中，相位變化可描述為：在二維聲學(xué)無限單元中橫向形函數(shù)可采用線性單元，不同邊界對應(yīng)的形函數(shù)為（如圖2所示）：

相位函數(shù)μ（ ξ，ζ）可由極點位置a（ξ）插值得到，a（ξ）為如下表式：

從（16）式中可看出，向外傳播的行進(jìn)波是由在近場中的虛擬極點發(fā)散出來的，虛擬極點位于有限單元與無界單元的交界面a（）ξ處。將模擬聲壓幅值衰減函數(shù)和相位變化函數(shù)結(jié)合起來，給出n階聲學(xué)單元的形函數(shù)如下：

波包絡(luò)法中權(quán)函數(shù)取為形函數(shù)的復(fù)數(shù)共軛乘以一加權(quán)因子（a/r）2，能保證無限單元中的剛度陣和質(zhì)量矩陣的積分有限，同時能夠滿足無窮遠(yuǎn)處輻射條件。故權(quán)函數(shù)的表達(dá)式為：

其中幾何權(quán)函數(shù)加權(quán)因子在局部坐標(biāo)系中可表示為：

基于文獻(xiàn)[6]可得到二維n階無限波包絡(luò)單元的剛度矩陣和質(zhì)量矩陣的系數(shù)分別為：

單元矩陣中，復(fù)雜指數(shù)形式在被積函數(shù)中已被消除，故可采用數(shù)值高斯積分得到各個單元矩陣中的系數(shù)。單元矩陣中的各項系數(shù)經(jīng)過組裝得到系統(tǒng)矩陣，帶入系統(tǒng)平衡方程即可求解得到所求節(jié)點聲壓值。

3 數(shù)值計算

3.1 映射變階波包絡(luò)聲學(xué)無限元法有效性驗證

基于映射變階波包絡(luò)聲學(xué)無限元法理論，求解聲壓計算框圖實現(xiàn)步驟見圖3。為滿足調(diào)整單元階數(shù)的靈活性，即通過增加徑向節(jié)點數(shù)目來實現(xiàn)聲學(xué)無限映射單元中徑向形函數(shù)的任意變階。然而，母單元中任意增加節(jié)點數(shù)，必須保證映射到實際單元中時是沿著無窮遠(yuǎn)邊界發(fā)散的，絕對不能出現(xiàn)交叉現(xiàn)象，即不存在當(dāng)ζ∈[-1，1）時，通過方程（7）得到點的坐標(biāo)不能位于1-4或2-3連線上。n階聲學(xué)映射無限元單元的位置是由它的四個幾何節(jié)點和n階徑向插值形函數(shù)確定。在n階無限單元徑向上至少需要（n+1）個高斯積分點，橫向方向上至少需要3個高斯積分點才能滿足計算精度要求。

圖3 無限元聲壓矩陣計算框圖Fig.3 Acoustic infinite element pressure matrix computation block diagram

本文基于映射變階波包絡(luò)聲學(xué)無限元法通過數(shù)值計算給出三維無限長剛性體運動的聲輻射。取浸沒在無限域流場的一個二維軸對稱無限長圓柱殼體為研究對象，計算求解無限聲場域中聲壓。將數(shù)值計算結(jié)果與文獻(xiàn)[6]解析解對比來驗證本文方法的有效性。以下所有算例中，聲學(xué)介質(zhì)密度為ρ= 1 024 kg/m3，聲傳播速度為c=1 500 m/s，如不作特殊說明，這些參數(shù)保持不變。

3.1.1 解析解

二維單極子、偶極子和四極子在柱坐標(biāo)系中的聲輻射解析解可表示為：

其中：V0為振速幅值，V0=1e-6 m/s。

3.1.2 數(shù)值解

在無量綱化參數(shù)ka=π a=1（）m下，偶極子、四極子聲輻射在聲場中的聲壓分布圖如圖4、圖5所示。確定聲學(xué)波數(shù)k使圓柱體直徑2a對應(yīng)一個波長。法向速度邊界條件V0=1e-6 m/s。

由圖4可知，通過將偶極子聲輻射數(shù)值解與解析解對比發(fā)現(xiàn)二者吻合程度很好且計算效率高；隨著波包絡(luò)單元階數(shù)提高，誤差變得更小。由此可判定聲學(xué)無限元方法能夠很好地滿足計算精度要求。

其中：Hn為第二類n階漢克爾函數(shù)，n=0，1，2分別對應(yīng)單極子、偶極子和四極子聲輻射。

在無限長圓柱體振動表面的法向速度分布為：

圖4 振動圓柱體輻射聲壓（ka=π）-+-+-，解析解；-o-o-，數(shù)值解Fig.4 The radiated sound pressure of vibration cylinder(ka=π).-+-+-,analytical solution;-o-o-,numerical solution

圖5 振動圓柱體輻射聲壓（ka=π）-o-o-，解析解；+*+*，數(shù)值解Fig.5 The radiated sound pressure of vibration cylinder(ka=π).-o-o-,analytical solution;+*+*,numerical solution

由圖5可知，通過將四極子聲輻射數(shù)值解與解析解對比發(fā)現(xiàn)二者吻合程度也很好；同時發(fā)現(xiàn)徑向形函數(shù)的階數(shù)越高，二者的吻合程度相當(dāng)高。從偶極子和四極子算例中，分析發(fā)現(xiàn)將映射變階波包絡(luò)聲學(xué)無限元法應(yīng)用到求解結(jié)構(gòu)外聲場聲壓是可行的，且具有較高的精度。

3.2 虛擬極點位置對結(jié)構(gòu)聲輻射的影響

在映射波包絡(luò)無限元法中，假設(shè)了一個在物理意義上不存在的虛擬聲源（即極點），其中極點作為聲源的擾動起點具有重要作用。一般來說，虛擬極點聲源都放置在幾何中心處，恰似結(jié)構(gòu)的輻射聲波從虛擬極點發(fā)出，但在實際工程問題中，很難滿足這樣的要求。因此有必要討論虛擬極點位置對結(jié)構(gòu)聲輻射的影響。本節(jié)以簡單聲源（四極子）為例分析了虛擬極點位置對聲輻射的影響。

在無量綱常數(shù)ka=π a=1（）m下，極點偏心距離取X=-0.5 m時，四極子聲學(xué)無限元聲輻射虛擬極點偏心與未偏心數(shù)值解對比如圖6所示。

圖6 虛擬極點位置對輻射聲壓的影響（ka=π）-+-+-，虛擬極點未偏心；-o-o-，虛擬極點偏心數(shù)值解Fig.6 The virtual sound source’s position influence on radiated sound pressure(ka=π).-+-+-,virtual sound source in center;-o-o-,numerical solution of virtual sound source with decentration

如圖6所示，虛擬極點的位置不同對輻射聲壓的影響不同，通過對比四極子極點偏心聲壓值發(fā)現(xiàn)，當(dāng)虛擬極點存在偏心時，對聲壓計算精度有不利影響。但隨聲學(xué)無限單元的階數(shù)增加，誤差會變小，其數(shù)值解會接近于真實解。因此，在實際工程問題中應(yīng)增加聲學(xué)無限單元的階數(shù)來提高計算精度，控制虛擬極點的偏心距離來減小誤差。

4 結(jié)論

本文基于映射變階波包絡(luò)聲學(xué)無限元法，以無限長圓柱殼體為研究對象，將映射型聲學(xué)無限元法與有限元法耦合求解無界域遠(yuǎn)場的輻射聲壓，提出了結(jié)構(gòu)近場可用有限元模擬，遠(yuǎn)場可用聲學(xué)無限元法離散的模型，運用映射型聲學(xué)無限元法數(shù)值求解了研究對象的聲輻射，與解析解對比驗證了映射變階波包絡(luò)法的有效性，并進(jìn)一步討論虛擬極點位置對聲輻射的影響。結(jié)論如下：

（1）在無限流場域中具有法向振速的結(jié)構(gòu)外取一簡單的人工邊界，在人工邊界內(nèi)采用有限元離散求解邊界節(jié)點信息，而在人工邊界外敷設(shè)一層聲學(xué)無限單元，在無限單元中通過耦合邊界節(jié)點信息，從而采用徑向映射變階的波包絡(luò)無限元法求解遠(yuǎn)場輻射聲壓。

（2）基于映射波包絡(luò)無限元法的理論，通過數(shù)值計算對比驗證了映射變階波包絡(luò)聲學(xué)無限元法在計算結(jié)構(gòu)外聲場聲壓的有效性，分析發(fā)現(xiàn)映射變階波包絡(luò)聲學(xué)無限元法具有較高的計算精度及求解效率；

（3）虛擬極點偏心對聲輻射計算精度有較大影響，因此在實際工程問題中，應(yīng)增加聲學(xué)無限單元的階數(shù)來滿足計算精度要求，同時也要控制減小虛擬極點的偏心距離。

[1]Burnett D S.A three-dimensional acoustic Infinite Element based on a prolate spheroidal multipole expansion[J].Acoust. Sot.Am,1996.

[2]Bettess P.Infinite Elements[M].Penshaw,Sunderland,U.K,1992.

[3]Astley R J,Eversman W.Wave envelope elements for acoustical radiation in inhomogeneous media[J].Computers and Structures,1988,30(4):801-810.

[4]Burnett D S,Holford R L Prolate and oblate spheroidal acoustic infinite element[J].Comput.Methods.Appl.Mech.Engrg, 1998,158:117-141.

[5]Tsynkov S V.Numerical solution of problems on unbounded domains[J].A Review Applied Numerical Mathematics,1998, 27:465-532.

[6]Cremers L,Fyfe K R,Coyette J P.A variable order infinite acoustic wave envelope element[J].Sound Vib,1994,171 (114).

[7]Cremers L,Fyfe K R.On the use of variable order infinite wave envelope element for acoustic radiation and scattering[J]. Acoust.Soc.Am.,1995,97(4):2028-2040.

[8]楊瑞梁,汪鴻振.聲無限元進(jìn)展[J].機(jī)械工程學(xué)報,2003,39(11):82-87.

[9]楊瑞梁.聲學(xué)無限元法的研究[D].上海:上海交通大學(xué),2004.

[10]楊瑞梁,汪鴻振.聲無限元進(jìn)展[J].機(jī)械工程學(xué)報,2003,39(11):82-87.

[11]楊瑞梁,范曉偉.使用有限元和無限元耦合求解聲輻射問題[J].振動工程學(xué)報,2004,17:1007-1009.

[12]楊瑞梁,汪鴻振.用長旋轉(zhuǎn)橢球函數(shù)重構(gòu)聲場[J].噪聲與振動控制,2003,5:15-17.

Variable order mapped wave envelope element of outer structural sound field

SU Nan,PANG Fu-zhen,YAO Xiong-liang
(Harbin Engineering University,Harbin 150001,China)

This paper is on the basis of a variable order in radial direction infinite acoustic wave envelope element(WEE).Numerical method is adopted to study the sound radiation pressure of an infinitely long cylinder as the subject investigated,and obtained results agree well with analytic solutions.The research presents that the variable order infinite acoustic wave envelope element which can be successfully applied to sound radiation pressure computation field is proven.The WEE has obvious advantages such as high efficiency and high degree of accuracy.Based on the above results,this paper also describes the influence of virtual sound source’s location in sound radiation,and it can be concluded that the decentration of virtual source has bad effects on the actual project by comparing the sound pressure values of quadrapole.However,with the increased orders of infinite acoustic element,the inaccuracy of sound pressure is decreasing.

variable order mapped wave envelope elements;infinite acoustic element; virtual sound source;sound radiation

O42

A

10.3969/j.issn.1007-7294.2014.07.016

1007-7294（2014）07-0856-08

2014-03-02

國家自然科學(xué)基金項目（51209052）；黑龍江省青年科學(xué)基金資助項目（QC2011C013）；哈爾濱市

科技創(chuàng)新人才研究專項資金項目（2011RFQXG021）；國防預(yù)研重點項目（401040XXX0103）

蘇楠（1976-），女，碩士研究生，E-mail：sulanlan1988@126.com；

龐福振（1981-），男，博士，哈爾濱工程大學(xué)副教授。