繆碩， 石敏， 任春雨

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

0 引 言

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

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

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

本文以脈動球源和剛性小球為對象，研究單元類型及單元大小等對聲學問題計算結果的影響，為工程中的聲學計算提供參考。

1 計算方法

有限元-無限元耦合方法可以在結構表面施加載荷，模擬流-固耦合作用并利用無限元單元求解外場的聲輻射特性；也可以設置入射聲波，模擬聲波與結構的相互作用，并利用無限元單元求解外場的聲散射特性。有限元-無限元耦合方法利用人工截斷邊界將無限大求解域分為內部域和外部無限大域，如圖1所示，內部域包含流體域和結構域，結構域內部可以包括復雜內部結構，內部域采用有限元單元離散；外部無限大域采用無限元單元離散，并在截斷邊界處設置無限元單元，以消除邊界影響并求解外聲場。

圖 1 有限元-無限元耦合方法原理示意

1.1 有限元-無限元耦合方法

對于理想聲學介質，聲波在其中的傳播需要滿足Helmholtz方程

(1)

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

(2)

設結構域表面的加速度大小為，則在結構-聲學邊界處為保證其連續(xù)性，聲場和結構應當滿足：

(3)

式中：為流場密度；為聲結構邊界法向。

應用Galerkin加權殘值法對聲場進行求解，首先設聲壓的近似解

(4)

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

對Helmholtz方程、聲結構表面處邊界及截斷處邊界取場函數(shù)及邊界條件的加權積分，有

(5)

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

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

(+ik-)=

(6)

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

剛度矩陣、阻尼矩陣、質量矩陣、聲學載荷矩陣的元素表達式為

(7)

(8)

(9)

(10)

由式(7)～(10)可知，當選取權函數(shù)為形函數(shù)的共軛時，可以使得求解過程得到極大的簡化，以下著重討論形函數(shù)的選取。

2.2 單元的形函數(shù)

以二維典型單元為例，說明有限元及無限元單元形函數(shù)的形式，以及2種單元在相交邊界處的關系。在內部域采用有限元進行離散，在局部坐標系-下形函數(shù)的表達式為

(,)=(,)

(11)

為表述方便，選取4節(jié)點矩形單元，單元及各點坐標見圖2，各節(jié)點形函數(shù)的表達式見式(12)～(15)，對應圖形見圖3。

圖 2 局部坐標系4節(jié)點矩形單元示意

(12)

(13)

(14)

(15)

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

圖4為無限元單元示意，其中，圖4(a)為局部坐標系-下的無限元單元，圖4(b)為全局坐標系-下的無限元單元，1、2為人工截斷邊界上的節(jié)點，2個坐標系下無限元各節(jié)點相對應。2種坐標系之間的映射函數(shù)見文獻[9]，通過坐標變換使得全局坐標系中無窮遠處對應局部坐標系中=1。

圖 4 聲學無限元單元

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

(16)

式中：

(17)

(18)

(19)

()為無限元單元為沿局部坐標系軸的形函數(shù)，,為沿局部坐標系軸方向上的拉格朗日多項式，為局部坐標系軸方向上的節(jié)點數(shù)。

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

圖 5 坐標變換后有限元單元和無限元單元

(20)

(21)

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

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

3 計算方法影響要素

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

3.1 聲輻射問題精度

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

圖 7 聲輻射問題示意及模型

通過改變單元類型、單元大小探究其對輻射聲場計算結果的影響。單元大小為0.03和0.02 m時，三角形和四邊形單元(均包含一次單元和二次單元)的脈動球源網(wǎng)格模型見圖8，一次單元與二次單元的網(wǎng)格相同。

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

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

第一種形式：

(22)

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

第二種形式：

(23)

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

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

圖 9 不同單元類型及節(jié)點總數(shù)下輻射聲壓無限元計算結果

表 1 聲輻射問題滿足精度對應的節(jié)點數(shù)及單元大小

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

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

3.2 聲散射問題精度

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

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

圖 10 5 kHz下不同單元類型、大小散射聲壓無限元計算結果

由圖10可知，當4種單元類型達到一定數(shù)量時，均可以較好地求解散射聲場。為進一步分析其收斂特性，以5%誤差為標準列出滿足精度要求的節(jié)點數(shù)以及對應單元大小，結果見表2。

表 2 5 kHz下聲散射問題滿足精度對應的節(jié)點總數(shù)及單元大小

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

4 結 論

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

(1)有限元-無限元耦合方法適用于聲輻射、聲散射問題的計算分析，且計算結果滿足精度要求。

(2)二次單元的單元大小僅需為求解頻率對應波長的1/3左右，計算結果即可滿足一定精度要求?？傊ㄟ^構建二次單元可以提高計算效率。

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