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

蘇楠，龐福振，姚熊亮

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

蘇楠，龐福振，姚熊亮

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

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

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

1 引言

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

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

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

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

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

圖1 聲場(chǎng)幾何模型Fig.1 The geometry model of sound field

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

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

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

式中：U（）

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

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

2.2 方程離散

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

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

如果在遠(yuǎn)場(chǎng)邊界S2上，S2距離振動(dòng)體足夠遠(yuǎn)，等效近似有限聲學(xué)邊界條件為：1、2、3、4分別對(duì)應(yīng)母單元中的四個(gè)節(jié)點(diǎn)，其他兩節(jié)點(diǎn)（1，1）和（1，-1）映射到無(wú)限單元的無(wú)窮遠(yuǎn)處。節(jié)點(diǎn)1′、2′與節(jié)點(diǎn)3、4關(guān)于節(jié)點(diǎn)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邊距聲源點(diǎn)1′或2′測(cè)得的距離。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值計(jì)算

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

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

圖3 無(wú)限元聲壓矩陣計(jì)算框圖Fig.3 Acoustic infinite element pressure matrix computation block diagram

本文基于映射變階波包絡(luò)聲學(xué)無(wú)限元法通過數(shù)值計(jì)算給出三維無(wú)限長(zhǎng)剛性體運(yùn)動(dòng)的聲輻射。取浸沒在無(wú)限域流場(chǎng)的一個(gè)二維軸對(duì)稱無(wú)限長(zhǎng)圓柱殼體為研究對(duì)象，計(jì)算求解無(wú)限聲場(chǎng)域中聲壓。將數(shù)值計(jì)算結(jié)果與文獻(xiàn)[6]解析解對(duì)比來(lái)驗(yàn)證本文方法的有效性。以下所有算例中，聲學(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ù)值解

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

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

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

在無(wú)限長(zhǎng)圓柱體振動(dòng)表面的法向速度分布為：

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

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

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

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

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

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

圖6 虛擬極點(diǎn)位置對(duì)輻射聲壓的影響（ka=π）-+-+-，虛擬極點(diǎn)未偏心；-o-o-，虛擬極點(diǎn)偏心數(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所示，虛擬極點(diǎn)的位置不同對(duì)輻射聲壓的影響不同，通過對(duì)比四極子極點(diǎn)偏心聲壓值發(fā)現(xiàn)，當(dāng)虛擬極點(diǎn)存在偏心時(shí)，對(duì)聲壓計(jì)算精度有不利影響。但隨聲學(xué)無(wú)限單元的階數(shù)增加，誤差會(huì)變小，其數(shù)值解會(huì)接近于真實(shí)解。因此，在實(shí)際工程問題中應(yīng)增加聲學(xué)無(wú)限單元的階數(shù)來(lái)提高計(jì)算精度，控制虛擬極點(diǎn)的偏心距離來(lái)減小誤差。

4 結(jié)論

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

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

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

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

[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]楊瑞梁,汪鴻振.聲無(wú)限元進(jìn)展[J].機(jī)械工程學(xué)報(bào),2003,39(11):82-87.

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

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

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

[12]楊瑞梁,汪鴻振.用長(zhǎng)旋轉(zhuǎn)橢球函數(shù)重構(gòu)聲場(chǎng)[J].噪聲與振動(dòng)控制,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

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

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

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

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