外域聲場(chǎng)計(jì)算中的無(wú)限元方法

繆碩， 石敏， 任春雨

(1.華中科技大學(xué) 船舶與海洋工程學(xué)院，武漢 430074； 2.中國(guó)人民解放軍 第91977部隊(duì)，北京 102249)

0 引 言

外域聲場(chǎng)問(wèn)題可以歸結(jié)為無(wú)界域聲場(chǎng)的數(shù)值求解問(wèn)題，這類(lèi)問(wèn)題的分析不論是在空氣中還是在水中都有很多現(xiàn)實(shí)需求，如水下航行器的聲輻射問(wèn)題、聲散射問(wèn)題等。

目前，對(duì)于聲輻射、聲散射問(wèn)題的求解方法主要包括解析法、數(shù)值方法。解析法僅適用于較為簡(jiǎn)單的結(jié)構(gòu)形式，局限性較大；數(shù)值方法包括有限元方法、有限元-邊界元方法以及有限元-無(wú)限元耦合方法等。有限元方法主要進(jìn)行特殊的邊界層處理如設(shè)置PML等，但計(jì)算頻率越高、模型越大，其邊界層網(wǎng)格數(shù)就越多，影響計(jì)算效率。有限元-邊界元方法的方程求解存在奇異性、多值解等問(wèn)題，影響其在實(shí)際工程中的應(yīng)用。有限元-無(wú)限元耦合方法將無(wú)限域分為內(nèi)部域和外部域，內(nèi)部域采用有限元離散，在邊界處采用無(wú)限元單元求解外場(chǎng)聲學(xué)問(wèn)題，帶來(lái)的額外的計(jì)算自由度較少，適用范圍廣。

對(duì)于有限元-無(wú)限元方法，國(guó)內(nèi)外已經(jīng)有相當(dāng)多的研究。ASTLEY等在BURNETT等無(wú)限元思想基礎(chǔ)上提出一種新的無(wú)限元單元類(lèi)型，其將形函數(shù)的共軛作為伽遼金加權(quán)殘值法的權(quán)函數(shù)，消除方程中的不確定積分，簡(jiǎn)化求解方程。ASTLEY等進(jìn)一步改善其數(shù)學(xué)嚴(yán)謹(jǐn)性以及適用性，發(fā)展成為映射包絡(luò)元亦即Astley元，并求解偶極子、四極子等的聲場(chǎng)特性。CREMERS等將Astley元進(jìn)一步發(fā)展成可實(shí)現(xiàn)任意變階的波包絡(luò)單元。楊瑞梁等闡述無(wú)限元的發(fā)展現(xiàn)狀，并對(duì)無(wú)限元進(jìn)行歸納分析，細(xì)致地給出各方法的原理及問(wèn)題，并展望無(wú)限元的發(fā)展前景。龐福振利用無(wú)限元方法預(yù)報(bào)船舶的結(jié)構(gòu)噪聲，發(fā)現(xiàn)相對(duì)于無(wú)反射邊界條件，無(wú)限元方法精度更高，而且對(duì)于遠(yuǎn)場(chǎng)聲輻射求解方便、計(jì)算效率高。尹旭超提出一種可任意變階的聲學(xué)無(wú)限元，并基于MATLAB計(jì)算單極子等聲輻射以及柱殼聲散射問(wèn)題。蘇楠等和田正東等都驗(yàn)證無(wú)限元方法求解聲學(xué)問(wèn)題的準(zhǔn)確性，并利用無(wú)限元方法研究層圓柱殼層厚度、內(nèi)部結(jié)構(gòu)等對(duì)于圓柱殼的均方振速及輻射聲功率的影響。JIANG等采用有限元-無(wú)限元方法分析滾筒洗衣機(jī)的低頻噪聲，并將計(jì)算結(jié)果與實(shí)驗(yàn)方法測(cè)量的輻射噪聲進(jìn)行對(duì)比，發(fā)現(xiàn)計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。MOHEIT等為避免在低頻計(jì)算時(shí)需要構(gòu)建較大的完美匹配層(PML)，通過(guò)算例驗(yàn)證無(wú)限元的有效性，并利用其開(kāi)展聲子晶體在自由聲場(chǎng)中特性的分析；同時(shí)利用有限元-無(wú)限元耦合方法探究有限元單元大小及無(wú)限元單元在無(wú)窮遠(yuǎn)方向上的插值多項(xiàng)式及插值節(jié)點(diǎn)數(shù)對(duì)外場(chǎng)聲學(xué)模態(tài)的影響。吳健等討論采用網(wǎng)格由細(xì)到粗的過(guò)渡劃分策略對(duì)計(jì)算結(jié)果的影響，驗(yàn)證合理劃分網(wǎng)格有利于計(jì)算。目前，有限元-無(wú)限元耦合方法以聲輻射特性的研究為主，研究對(duì)象包括柱殼、錐殼、船體結(jié)構(gòu)等，但關(guān)于聲散射特性以及有限元-無(wú)限元方法計(jì)算結(jié)果影響要素的研究較少。

本文以脈動(dòng)球源和剛性小球?yàn)閷?duì)象，研究單元類(lèi)型及單元大小等對(duì)聲學(xué)問(wèn)題計(jì)算結(jié)果的影響，為工程中的聲學(xué)計(jì)算提供參考。

1 計(jì)算方法

有限元-無(wú)限元耦合方法可以在結(jié)構(gòu)表面施加載荷，模擬流-固耦合作用并利用無(wú)限元單元求解外場(chǎng)的聲輻射特性；也可以設(shè)置入射聲波，模擬聲波與結(jié)構(gòu)的相互作用，并利用無(wú)限元單元求解外場(chǎng)的聲散射特性。有限元-無(wú)限元耦合方法利用人工截?cái)噙吔鐚o(wú)限大求解域分為內(nèi)部域和外部無(wú)限大域，如圖1所示，內(nèi)部域包含流體域和結(jié)構(gòu)域，結(jié)構(gòu)域內(nèi)部可以包括復(fù)雜內(nèi)部結(jié)構(gòu)，內(nèi)部域采用有限元單元離散；外部無(wú)限大域采用無(wú)限元單元離散，并在截?cái)噙吔缣幵O(shè)置無(wú)限元單元，以消除邊界影響并求解外聲場(chǎng)。

圖 1 有限元-無(wú)限元耦合方法原理示意

1.1 有限元-無(wú)限元耦合方法

對(duì)于理想聲學(xué)介質(zhì)，聲波在其中的傳播需要滿(mǎn)足Helmholtz方程

(1)

對(duì)于有限元離散流場(chǎng)，其離散域必定有邊界的存在，當(dāng)聲波在其內(nèi)傳播時(shí)，邊界阻抗的不匹配會(huì)導(dǎo)致聲波在邊界處產(chǎn)生反射。在有限元邊界處設(shè)置聲學(xué)無(wú)限元，可以消除邊界影響，模擬無(wú)限大流域，因此有限元-無(wú)限元耦合方法需要在邊界處滿(mǎn)足阻抗匹配邊界條件：

(2)

設(shè)結(jié)構(gòu)域表面的加速度大小為，則在結(jié)構(gòu)-聲學(xué)邊界處為保證其連續(xù)性，聲場(chǎng)和結(jié)構(gòu)應(yīng)當(dāng)滿(mǎn)足：

(3)

式中：為流場(chǎng)密度；為聲結(jié)構(gòu)邊界法向。

應(yīng)用Galerkin加權(quán)殘值法對(duì)聲場(chǎng)進(jìn)行求解，首先設(shè)聲壓的近似解

(4)

式中：為基函數(shù)或形函數(shù)(以下統(tǒng)一稱(chēng)為形函數(shù))；為待定系數(shù)，=1,2,…,。

對(duì)Helmholtz方程、聲結(jié)構(gòu)表面處邊界及截?cái)嗵庍吔?span id="j5i0abt0b" class="subscript">取場(chǎng)函數(shù)及邊界條件的加權(quán)積分，有

(5)

式中：為伽遼金方法權(quán)函數(shù)，=1,2,…,。

將聲壓近似解表達(dá)式寫(xiě)成形函數(shù)與待定系數(shù)的乘積，則式(5)可表示為

(+ik-)=

(6)

式中：為待定系數(shù)矩陣，其各元素為近似解表達(dá)式中的待定系數(shù)；為剛度矩陣；為阻尼矩陣；為質(zhì)量矩陣；為聲學(xué)載荷矩陣。

剛度矩陣、阻尼矩陣、質(zhì)量矩陣、聲學(xué)載荷矩陣的元素表達(dá)式為

(7)

(8)

(9)

(10)

由式(7)～(10)可知，當(dāng)選取權(quán)函數(shù)為形函數(shù)的共軛時(shí)，可以使得求解過(guò)程得到極大的簡(jiǎn)化，以下著重討論形函數(shù)的選取。

2.2 單元的形函數(shù)

以二維典型單元為例，說(shuō)明有限元及無(wú)限元單元形函數(shù)的形式，以及2種單元在相交邊界處的關(guān)系。在內(nèi)部域采用有限元進(jìn)行離散，在局部坐標(biāo)系-下形函數(shù)的表達(dá)式為

(,)=(,)

(11)

為表述方便，選取4節(jié)點(diǎn)矩形單元，單元及各點(diǎn)坐標(biāo)見(jiàn)圖2，各節(jié)點(diǎn)形函數(shù)的表達(dá)式見(jiàn)式(12)～(15)，對(duì)應(yīng)圖形見(jiàn)圖3。

圖 2 局部坐標(biāo)系4節(jié)點(diǎn)矩形單元示意

(12)

(13)

(14)

(15)

圖 3 4節(jié)點(diǎn)矩形單元各節(jié)點(diǎn)形函數(shù)

圖4為無(wú)限元單元示意，其中，圖4(a)為局部坐標(biāo)系-下的無(wú)限元單元，圖4(b)為全局坐標(biāo)系-下的無(wú)限元單元，1、2為人工截?cái)噙吔缟系墓?jié)點(diǎn)，2個(gè)坐標(biāo)系下無(wú)限元各節(jié)點(diǎn)相對(duì)應(yīng)。2種坐標(biāo)系之間的映射函數(shù)見(jiàn)文獻(xiàn)[9]，通過(guò)坐標(biāo)變換使得全局坐標(biāo)系中無(wú)窮遠(yuǎn)處對(duì)應(yīng)局部坐標(biāo)系中=1。

圖 4 聲學(xué)無(wú)限元單元

無(wú)限元單元的形函數(shù)(,)可以表示為

(16)

式中：

(17)

(18)

(19)

()為無(wú)限元單元為沿局部坐標(biāo)系軸的形函數(shù)，,為沿局部坐標(biāo)系軸方向上的拉格朗日多項(xiàng)式，為局部坐標(biāo)系軸方向上的節(jié)點(diǎn)數(shù)。

在給出的2種形函數(shù)基礎(chǔ)上，進(jìn)一步討論單元形函數(shù)在邊界處的情況。圖5為坐標(biāo)交換后的有限元單元和無(wú)限元單元，有限元單元在邊界處節(jié)點(diǎn)1、2的形函數(shù)為式(16)，無(wú)限元單元的形函數(shù)為式(20)～(21)。

圖 5 坐標(biāo)變換后有限元單元和無(wú)限元單元

(20)

(21)

將無(wú)限元域和有限元域二者形函數(shù)繪制在一起，如圖6所示。當(dāng)<-1時(shí)代表有限元單元的形函數(shù)，-1<<0時(shí)為無(wú)限元單元的形函數(shù)，二者在邊界處=-1時(shí)能夠保持連續(xù)。

圖 6 有限元、無(wú)限元形函數(shù)

3 計(jì)算方法影響要素

在實(shí)際的外域聲場(chǎng)計(jì)算中，由于模型規(guī)模大以及計(jì)算頻段寬等原因，計(jì)算量巨大，往往無(wú)法達(dá)到精確解的計(jì)算要求，有必要討論多種要素對(duì)于計(jì)算結(jié)果的影響。本文選取單元類(lèi)型、單元大小這2種影響要素，從聲輻射、聲散射2個(gè)方面進(jìn)行分析。

3.1 聲輻射問(wèn)題精度

聲輻射問(wèn)題的研究對(duì)象為某脈動(dòng)球源，半徑0.1 m，其外部為半徑0.2 m的球形水域，見(jiàn)圖7所示，在水域表面上設(shè)置無(wú)限元單元。本次的計(jì)算頻率為10 kHz，在計(jì)算軟件Abaqus中，在球源表面施加聲壓邊界條件=1 Pa。流體密度為1 000 kg/m，聲速為1 440 m/s，文中水域材料參數(shù)均保持一致。

圖 7 聲輻射問(wèn)題示意及模型

通過(guò)改變單元類(lèi)型、單元大小探究其對(duì)輻射聲場(chǎng)計(jì)算結(jié)果的影響。單元大小為0.03和0.02 m時(shí)，三角形和四邊形單元(均包含一次單元和二次單元)的脈動(dòng)球源網(wǎng)格模型見(jiàn)圖8，一次單元與二次單元的網(wǎng)格相同。

圖 8 不同單元類(lèi)型及單元大小的脈動(dòng)球源網(wǎng)格示意

一次和二次單元網(wǎng)格數(shù)相同，為比較二者對(duì)計(jì)算結(jié)果的影響，以節(jié)點(diǎn)數(shù)代表單元大小。本次計(jì)算結(jié)果以2種形式給出，以便對(duì)比單元類(lèi)型和大小對(duì)輻射聲場(chǎng)聲壓計(jì)算結(jié)果的影響。

第一種形式：

(22)

式中：代表有限元水域外表面上第個(gè)點(diǎn)的聲壓;為節(jié)點(diǎn)總數(shù)。代表水域外表面所有節(jié)點(diǎn)聲壓均值大小。

第二種形式：

(23)

式中：0.5代表水域表面聲壓大小理論解。表示水域外表面所有節(jié)點(diǎn)聲壓的誤差均方值。

2種形式的計(jì)算結(jié)果見(jiàn)圖9。由此可知，隨著節(jié)點(diǎn)數(shù)的增多，4種單元均可以較好地求解水域表面的輻射聲場(chǎng)。為進(jìn)一步分析，以5%的誤差為標(biāo)準(zhǔn)列出各網(wǎng)格類(lèi)型下、滿(mǎn)足精度要求對(duì)應(yīng)的節(jié)點(diǎn)數(shù)及單元大小，結(jié)果見(jiàn)表1。

圖 9 不同單元類(lèi)型及節(jié)點(diǎn)總數(shù)下輻射聲壓無(wú)限元計(jì)算結(jié)果

表 1 聲輻射問(wèn)題滿(mǎn)足精度對(duì)應(yīng)的節(jié)點(diǎn)數(shù)及單元大小

由表1可知，為滿(mǎn)足精度要求，4節(jié)點(diǎn)四邊形單元、3節(jié)點(diǎn)三角形單元的單元大小需要達(dá)到計(jì)算頻率對(duì)應(yīng)聲波波長(zhǎng)的1/6，8節(jié)點(diǎn)四邊形單元、6節(jié)點(diǎn)三角形單元的單元大小僅需為對(duì)應(yīng)聲波波長(zhǎng)的1/3及1/2.6，其原因是二次單元的插值形函數(shù)階數(shù)較高，可以較好地實(shí)現(xiàn)對(duì)聲場(chǎng)的模擬。

此外，和結(jié)果達(dá)到精度要求時(shí)，對(duì)應(yīng)的節(jié)點(diǎn)數(shù)及單元大小基本一致，只需選擇其中一種進(jìn)行計(jì)算即可。

3.2 聲散射問(wèn)題精度

聲散射問(wèn)題選取剛性小球?yàn)閷?duì)象，小球半徑大小、水域大小、網(wǎng)格劃分與聲輻射問(wèn)題中處理一致。計(jì)算頻率為1～10 kHz，在Abaqus中設(shè)置入射波形式為平面波，聲壓大小為1 Pa。

為便于后續(xù)對(duì)比分析，選取剛體表面上一點(diǎn)的聲壓，探究單元類(lèi)型、單元大小對(duì)計(jì)算結(jié)果的影響。不同單元類(lèi)型及大小下的模型示意圖與第3.1節(jié)中圖8類(lèi)似，此處不再贅述。5 kHz下基于有限元-無(wú)限元耦合方法的聲壓計(jì)算結(jié)果見(jiàn)圖10。

圖 10 5 kHz下不同單元類(lèi)型、大小散射聲壓無(wú)限元計(jì)算結(jié)果

由圖10可知，當(dāng)4種單元類(lèi)型達(dá)到一定數(shù)量時(shí)，均可以較好地求解散射聲場(chǎng)。為進(jìn)一步分析其收斂特性，以5%誤差為標(biāo)準(zhǔn)列出滿(mǎn)足精度要求的節(jié)點(diǎn)數(shù)以及對(duì)應(yīng)單元大小，結(jié)果見(jiàn)表2。

表 2 5 kHz下聲散射問(wèn)題滿(mǎn)足精度對(duì)應(yīng)的節(jié)點(diǎn)總數(shù)及單元大小

由表2可知，為滿(mǎn)足誤差要求，在5 kHz時(shí)4節(jié)點(diǎn)四邊形單元、3節(jié)點(diǎn)三角形單元的單元大小需要達(dá)到計(jì)算頻率對(duì)應(yīng)聲波波長(zhǎng)的1/7.5及1/6；而對(duì)于8節(jié)點(diǎn)四邊形單元、6節(jié)點(diǎn)三角形單元的單元大小僅需約為對(duì)應(yīng)聲波波長(zhǎng)的1/3，聲散射相比聲輻射對(duì)網(wǎng)格的要求稍高。

4 結(jié) 論

基于有限元軟件Abaqus，采用有限元-無(wú)限元耦合的方法開(kāi)展聲輻射、聲散射計(jì)算問(wèn)題的研究，討論該方法中單元類(lèi)型及單元尺寸等對(duì)聲場(chǎng)計(jì)算求解的影響，通過(guò)仿真結(jié)果發(fā)現(xiàn)：

(1)有限元-無(wú)限元耦合方法適用于聲輻射、聲散射問(wèn)題的計(jì)算分析，且計(jì)算結(jié)果滿(mǎn)足精度要求。

(2)二次單元的單元大小僅需為求解頻率對(duì)應(yīng)波長(zhǎng)的1/3左右，計(jì)算結(jié)果即可滿(mǎn)足一定精度要求?？傊ㄟ^(guò)構(gòu)建二次單元可以提高計(jì)算效率。

需要說(shuō)明的是，本文雖然僅討論無(wú)限元方法在外域聲場(chǎng)計(jì)算中的應(yīng)用，但該方法同樣適用于內(nèi)域問(wèn)題。另外，本文的計(jì)算方法可應(yīng)用于工程結(jié)構(gòu)(如艦船與潛艇)的外域聲場(chǎng)計(jì)算，雖然許多工作還有待深入，如無(wú)限元外部形狀對(duì)計(jì)算的影響，以及考慮結(jié)構(gòu)的材料屬性時(shí)的結(jié)構(gòu)與流體耦合問(wèn)題等，這些問(wèn)題都將在后續(xù)的工作中進(jìn)行討論。