全局弱式無網(wǎng)格方法求解消聲器橫向模態(tài)

趙 輝，左繼強(qiáng)，劉成洋，方 智

（1.華中科技大學(xué) 船舶與海洋工程學(xué)院，武漢 430074；2.中國艦船研究設(shè)計中心，武漢 430064）

求解消聲器橫向模態(tài)常用的方法有解析方法和數(shù)值方法。解析方法[1]求解簡單、快速、精確，但是僅局限于有解析方程的橫截面，比如圓形截面，矩形截面以及橢圓形截面。對于任意形狀的橫截面橫向模態(tài)的求解，只能使用數(shù)值方法[2]，現(xiàn)在較為成熟的常用的數(shù)值方法為二維有限元方法（FEM）和邊界元方法(BEM)。該類方法均需要依賴于網(wǎng)格進(jìn)行求解，是基于網(wǎng)格的求解方法。無網(wǎng)格方法(MFM)[3]是一種相對較新的數(shù)值方法，該方法不需要劃分網(wǎng)格，采用基于點(diǎn)的近似，利用一組散布在問題域中以及域邊界上的節(jié)點(diǎn)表示該問題域及其邊界。無網(wǎng)格方法最先在力學(xué)領(lǐng)域應(yīng)用較多，近十幾年開始在聲學(xué)領(lǐng)域被應(yīng)用。無網(wǎng)格法對聲學(xué)問題的求解，主要集中在對波動方程（Helmholtz方程）和線性歐拉方程的求解。1998年，Bouillar[4]應(yīng)用基于移動最小二乘的無網(wǎng)格伽遼金法計算了二維內(nèi)部聲學(xué)問題，數(shù)值算例比較結(jié)果表明該方法比有限元法的計算精度高，計算速度快。2002年，Chen等應(yīng)用邊界配點(diǎn)的無網(wǎng)格方法計算了二維聲腔的聲學(xué)特征值問題，并指出該方法較有限元法簡單，容易實施[5]。2010年，大連理工大學(xué)的張宏偉等應(yīng)用基于點(diǎn)插值的配點(diǎn)型無網(wǎng)格法求解了Helmholtz方程，通過數(shù)值算例比較證明了該方法具有較高的精度和良好的收斂性[6]。2011年，湖南大學(xué)的姚凌云等將分區(qū)光滑徑向點(diǎn)插值無網(wǎng)格法應(yīng)用于二維聲學(xué)分析中，對管道和二維轎車聲學(xué)問題的數(shù)值分析結(jié)果表明，該方法比有限元法具有更高的精度[7]。2012年，海軍工程大學(xué)的胡海等將邊界無網(wǎng)格法應(yīng)用到結(jié)構(gòu)聲輻射計算中，通過與解析結(jié)果對比表明該方法具有更高的插值和計算精度[8]。華中科技大學(xué)的黃其柏課題組將徑向配點(diǎn)無網(wǎng)格法應(yīng)用到氣動聲學(xué)和腔體結(jié)構(gòu)聲學(xué)的研究中，取得很好的成果[9]。2016年，曾向陽等提出一種基于聲粒子分布積分的無網(wǎng)格數(shù)值算法計算小尺度封閉空間中低頻聲場的數(shù)值方法，指出該方法特別適用于飛機(jī)座艙等關(guān)注中頻段聲場分布的情況[10]。2017，方智等提出應(yīng)用強(qiáng)式無網(wǎng)格配點(diǎn)方法耦合模態(tài)匹配技術(shù)計算任意橫截面膨脹腔消聲器的傳遞損失，計算結(jié)果與實驗測量值吻合較好，且計算速度較三維有限元方法具有一定的優(yōu)勢[11]。但當(dāng)導(dǎo)數(shù)邊界條件較多時，使用強(qiáng)式無網(wǎng)格配點(diǎn)方法求解消聲器橫向模態(tài)穩(wěn)定性會變差。

本文基于無網(wǎng)格方法本身的優(yōu)勢，使用徑向基函數(shù)點(diǎn)插值法離散本征方程，應(yīng)用伽遼金加權(quán)殘數(shù)法進(jìn)行數(shù)值積分，構(gòu)成全局弱式無網(wǎng)格方法用來求解消聲器橫向本征波數(shù)以及模態(tài)形狀，更方便添加導(dǎo)數(shù)邊界條件。進(jìn)而分析研究形狀參數(shù)對計算精度和速度的影響。

1 全局弱式無網(wǎng)格方法求解橫向模態(tài)

圖1為任意形狀等截面直通穿孔管阻性消聲器示意圖，將其分為三個區(qū)域：進(jìn)口管A、膨脹腔C和出口管E，對應(yīng)的橫截面分別為S1、S3=S1+S2和S1。膨脹腔被穿孔管分為兩部分：C1和C2，對應(yīng)的橫截面分別為S1和S2，C2內(nèi)填充吸聲材料。

圖1 任意形狀等截面直通穿孔管阻性消聲器

以膨脹腔C為例，介紹全局弱式無網(wǎng)格方法求解橫向模態(tài)的基本思想。膨脹腔內(nèi)橫向聲壓控制方程滿足以下形式

其中：為二維直角坐標(biāo)系下的拉普拉斯算子，pxy1和pxy2分別為區(qū)域C1和C2內(nèi)的橫向聲壓分量，kxy和分別表示空氣和吸聲材料中的橫向波數(shù)，并且滿足以下關(guān)系

其中：kz為軸向波數(shù)k=ω/c0為空氣中的波數(shù)，ω為圓頻率，c0為空氣中的聲速，=ω為吸聲材料中的復(fù)波數(shù)，為吸聲材料中的復(fù)聲速。

穿孔膨脹腔的邊界條件可以表示為以下形式，外部環(huán)形腔剛性壁邊界條件

穿孔面邊界條件

其中：ρ0為空氣密度為多孔吸聲材料密度為一側(cè)有吸聲材料時穿孔結(jié)構(gòu)的特性聲阻抗率

其中：?為穿孔率，z=ρ0c0和分別空氣和吸聲材料的特性阻抗。為單個孔的特性聲阻，μ為動力黏度系數(shù)，tw為穿孔管厚度，dh為穿孔孔徑，α為單孔厚度修正系數(shù)，可以表示為[12

其中：各個參數(shù)的含義可參考文獻(xiàn)[12]。

1.1 無網(wǎng)格徑向基點(diǎn)插值法求解形函數(shù)

為了避免采用多項式基PIM所引起的奇異性問題，徑向基函數(shù)（RBF）被采用以形成徑向基點(diǎn)插值法（RPIM）形函數(shù)用于弱式無網(wǎng)格法。添加多項式的RPIM插值可以表示為

其中：Ri(x)為徑向基函數(shù)，n為RBFs的個數(shù)，qk(x)為空間坐標(biāo)xT=[x,y]中的單項式，m為多項式基函數(shù)的個數(shù)。

為了確定系數(shù)矩陣中的ai和ck，需要形成計算點(diǎn)x的支持域，其中包括n個場節(jié)點(diǎn)。支持域的大小由ds=αsdc確定，αs為支持域的無量綱尺寸，應(yīng)該由分析者事先確定，dc為計算點(diǎn)x支持域內(nèi)的平均節(jié)點(diǎn)間距。在徑向基函數(shù)Ri(rk)中，αc為形狀參數(shù)。第i個場節(jié)點(diǎn)和第k個場節(jié)點(diǎn)之間的距離rk為

使方程滿足計算點(diǎn)x周圍的n個節(jié)點(diǎn)值，一個節(jié)點(diǎn)對應(yīng)一個方程，將產(chǎn)生n個線性方程。求解以下方程即可得到系數(shù)矩陣a0。

其中，各個矩陣表達(dá)式如下

很明顯，矩陣G是對稱矩陣，所以其可逆，因此系數(shù)矩陣a0可以由方程求得

將系數(shù)矩陣代入到方程，可以得到

對應(yīng)節(jié)點(diǎn)聲壓向量的RPIM形函數(shù)Φ表達(dá)式為

最終，橫向聲壓ph()x可以表示為

1.2 基于Galerkin弱式離散系統(tǒng)方程

應(yīng)用加權(quán)殘數(shù)法構(gòu)造橫向本征方程和邊界條件的弱式無網(wǎng)格形式。采用伽遼金加權(quán)殘數(shù)方法，權(quán)函數(shù)和形函數(shù)選取一致，如同式(21)所示。本征方程和邊界條件的殘數(shù)乘以權(quán)函數(shù)并且在各自域內(nèi)積分，可以得到加權(quán)形式

應(yīng)用格林公式，方程簡化為以下線性方程

其中：K1,M1,Z1，K2,M2,Z2分別為橫截面S1和S2內(nèi)的廣義剛度矩陣、質(zhì)量矩陣和穿孔阻抗矩陣。為了求解這些系統(tǒng)矩陣，需要求解域內(nèi)和邊界上的積分。因此在弱式無網(wǎng)格方法中需要構(gòu)造背景網(wǎng)格。背景網(wǎng)格與場節(jié)點(diǎn)是相互獨(dú)立的，并且對網(wǎng)格質(zhì)量沒有太高要求。

在每一個背景網(wǎng)格單元中配置高斯點(diǎn)，對于每一個高斯節(jié)點(diǎn)形成支持域，求解每一個高斯點(diǎn)的形函數(shù)Φ，進(jìn)而計算背景網(wǎng)格單元的剛度矩陣，質(zhì)量矩陣以及阻抗矩陣。

其中：nQ為單元e內(nèi)的高斯點(diǎn)個數(shù)，wQ為相應(yīng)的高斯加權(quán)系數(shù)，Je為單元的雅克比矩陣。

在弱式無網(wǎng)格方法中，背景網(wǎng)格單元一般選取規(guī)則形狀，四邊形或者三角形。在本文中，采用四邊形單元，每個單元中高斯點(diǎn)取4個。

對于不含有穿孔的橫截面，穿孔阻抗等于零，化簡方程可以求得橫向本征方程為

求解式(24)和式(28)可以求得穿孔橫截面和非穿孔橫截面的橫向模態(tài)。

1.3 求解本征方程

假設(shè)橫截面上的節(jié)點(diǎn)數(shù)為n，通過式和可以分別求得膨脹腔和進(jìn)出口管道軸向波數(shù)向量kz和特征向量矩陣Xxy，其中kz的維數(shù)為n，Xxy的階數(shù)為n×n。橫向節(jié)點(diǎn)聲壓分量組成的向量可以表示為，其中φ代表模態(tài)幅值系數(shù)組成向量，其維數(shù)為n，i為模態(tài)階數(shù)。

橫向聲壓分量可以表示為

最終，可以得到本征函數(shù)的表達(dá)式為

2 計算精度與效率驗證

為了驗證本文方法的計算精度和效率，本文計算了圓形截面，不規(guī)則截面和穿孔截面的本征橫向波數(shù)。計算結(jié)果分別與有限元結(jié)果和解析結(jié)果進(jìn)行對比，本征值的相對誤差被定義為如下形式

2.1 圓形截面

因為規(guī)則橫截面有橫向本征波數(shù)的解析解，所以如圖2所示的圓形截面最先被分析，圓形截面半徑為r=0.0245 m。截面內(nèi)的場節(jié)點(diǎn)的分布是通過簡單的MATLAB程序構(gòu)造的，背景網(wǎng)格由商業(yè)軟件ANSYS得到。由圖2可以看出，場節(jié)點(diǎn)是任意分布的，與背景網(wǎng)格是相互獨(dú)立的。

圖2圓形橫截面內(nèi)任意分布的場節(jié)點(diǎn)和背景網(wǎng)格

圖3 給出了分別使用全局弱式無網(wǎng)格方法和有限元方法計算的前5階本征波數(shù)的相對誤差值。

圖3 圓形截面本征值計算結(jié)果比較

由圖可以看出，隨著橫截面內(nèi)節(jié)點(diǎn)數(shù)的增加，仿真結(jié)果與解析結(jié)果的誤差逐漸減小。在無網(wǎng)格方法中，節(jié)點(diǎn)數(shù)為50時，相對誤差可以控制在0.5%以內(nèi)，需要的計算時間為15 s，而在有限元方法中，節(jié)點(diǎn)數(shù)為57時，第5階模態(tài)的本征值相對誤差為1.8%。為了將前5階模態(tài)本征值的相對誤差控制在0.5%以內(nèi)，在有限元方法中需要160個節(jié)點(diǎn)，計算時間為60 s。

2.2 非規(guī)則截面

為了證明本文方法的適用性和準(zhǔn)確性，本小節(jié)計算了如圖4所示的非規(guī)則截面的本征值。四邊形的具體尺寸如圖上所標(biāo)，長度單位為m。

圖4 非規(guī)則四邊形橫截面

圖5 不規(guī)則四邊形截面本征值計算結(jié)果比較

由于該截面為非規(guī)則截面，其本征頻率沒有解析解。本文采用二維有限元方法和無網(wǎng)格方法計算其前30階本征頻率。在有限元方法中，隨著劃分節(jié)點(diǎn)數(shù)的增多，計算精度增高，經(jīng)計算得知，當(dāng)節(jié)點(diǎn)數(shù)取522時，再增加節(jié)點(diǎn)數(shù)，精度增大不大，圖線中給出了節(jié)點(diǎn)數(shù)取1000和522時，前30階本征頻率仿真值的相對誤差均在0.5%以內(nèi)，說明有限元方法中節(jié)點(diǎn)數(shù)取522，可以得到較高的計算精度，需要耗時178 s。為得到相同精度的本征值，無網(wǎng)格方法中需要節(jié)點(diǎn)數(shù)為90。若背景網(wǎng)格單元采用有限元方法形成，即每一個節(jié)點(diǎn)都被包含在單元中，背景網(wǎng)格單元數(shù)為74，計算本征值所需時間為73 s。由于無網(wǎng)格方法對單元的依賴性較小，背景網(wǎng)格單元可以任意生成，本文中針對本算例采用30個單元作為背景網(wǎng)格，可以得到相同精度的本征值計算結(jié)果，所需時間為36 s，時間減少一半。為了對比，圖中給出了有限元方法中節(jié)點(diǎn)數(shù)取110時的本征值的相對誤差，相對誤差隨著節(jié)點(diǎn)數(shù)的增加基本呈增大趨勢，在第10階本征值以后，相對誤差超過5%，計算結(jié)果不再精確。

2.3 圓形穿孔截面

為了進(jìn)一步驗證本文全局弱式無網(wǎng)格方法的適用性和準(zhǔn)確性，本小節(jié)研究了如圖6所示的圓形穿孔截面。橫截面的具體尺寸為：穿孔管直徑d=0.0249 m，膨脹腔直徑D=0.1644 m。

圖6 圓形穿孔阻性消聲器橫截面

在本文中所用吸聲材料為玻璃絲綿，密度為100 g/L。其特性聲阻抗和復(fù)波數(shù)的表達(dá)式為[14]

由公式可知，每一個激發(fā)頻率都有相應(yīng)的本征值，即軸向本征波數(shù)。本文計算了激發(fā)頻率為3000 Hz時前5階徑向模態(tài)的軸向本征波數(shù)，表1中給出了分別使用解析方法，二維有限元方法和本文無網(wǎng)格方法計算的本征值，可以看出達(dá)到一定的精度，有限元方法需要1600節(jié)點(diǎn)構(gòu)成八節(jié)點(diǎn)四邊形單元，需要時間為445 s，而無網(wǎng)格方法只需140個節(jié)點(diǎn)，需要耗時75 s即可達(dá)到相同的精度。

表1 圓形穿孔阻性消聲器橫截面的軸向本征波數(shù)

通過以上3個算例可以看出，相對于有限元方法，無網(wǎng)格方法求解橫截面本征模態(tài)，不依賴于網(wǎng)格，需要的節(jié)點(diǎn)更少，耗時更少。橫截面面積越大，邊界條件越復(fù)雜，無網(wǎng)格方法的優(yōu)勢就越突出。

3 計算結(jié)果與分析

由前面分析可以知道，每個計算節(jié)點(diǎn)支持域的尺寸由形狀參數(shù)αs決定，圖7分析了支持域尺寸對本征值計算精度的影響。

由圖7可以看出，在低頻段內(nèi)計算本征值，支持域大小對計算結(jié)果影響不大，隨著計算頻率的增高，支持域的影響越來越明顯。隨著支持域逐漸增大，相對誤差整體趨勢逐漸變小，當(dāng)支持域的大小達(dá)到一定尺寸，相對誤差變化不大。針對本文中所有算例，支持域的大小取2.5倍dc時，計算前30階本征值誤差控制在0.5%。但是計算時間會隨著支持域尺寸的增大而增多，所以支持域尺寸的確定需要考慮計算精度和計算速度之間的平衡。

圖7 支持域尺寸對計算精度的影響

徑向基函數(shù)中另外一個形狀參數(shù)為αc。圖8給出了αc取7種不同數(shù)值時，非規(guī)則四邊形截面前30階本征頻率的相對誤差。

圖8 形狀參數(shù)αc對計算精度的影響

從圖中可以看出，在低頻段內(nèi)形參αc對本征值計算精度影響不大，隨著計算頻率增大，影響逐漸明顯。在整體頻率范圍內(nèi)，αc值取太小或太大都會使得相對誤差增大，這是因為矩陣條件數(shù)變大使得計算結(jié)果不穩(wěn)定的原因。針對本文所有算例，αc取值2～3.5之間計算本征值可以取得較高的精度。

4 結(jié)語

采用全局弱式無網(wǎng)格方法提取消聲器橫向本征值（橫向波數(shù)）和本征向量（聲壓模態(tài)），使用徑向基函數(shù)點(diǎn)插值方法用來離散本征方程，使用伽遼金加權(quán)殘數(shù)法用來數(shù)值積分，最終求得消聲器橫向本征函數(shù)。通過計算圓形截面，非規(guī)則四邊形截面和穿孔阻性消聲器橫截面的本征波數(shù)，分別與解析方法和有限元方法計算結(jié)果進(jìn)行對比，驗證了全局弱式無網(wǎng)格方法求解消聲器橫向本征模態(tài)具有較高的精度，并且需要更少的節(jié)點(diǎn)和計算時間。當(dāng)支持域尺寸大于2.5倍dc時，形狀參數(shù)αc在2～3.5之間取值時，計算得到的本征值精度較高。

與有限元方法相比，無網(wǎng)格方法中背景網(wǎng)格不依賴于節(jié)點(diǎn)，在計算過程中可以隨意添加會刪除某些節(jié)點(diǎn)，這也是其優(yōu)勢所在。

[1]SELAMET A,XU M.B,LEE I-J.Analytical approach for sound attenuation in perforated dissipative silencers with inlet/outlet extension[J].J.Acoust.Soc.Am.,2005,117(4):2078-2089.

[2]ZHAO S.On the spurious solutions in the high-order finite element methods for eigenvalue problems[J].Comput MthodAppl.M.,2007,196(49-52):5031-5036.

[3]LIU G R,GU Y T.An introduction to meshfree methods and their programming[M].New York:Springer,2005.

[4]BOUILLARD P,SULEAUB S.Element-free Galerkin solutions for helmholtz problems:formulation and numerical assessment of the pollution effect[J].Comput MthodAppl.M.,1998,162:317-335.

[5]CHEN J T,CHANG M H,CHEN K H,at al.Boundary collocation method for acoustic eigen analysis of threedimensionalcavities using radialbasis function[J].Comput Mech.,2002,29(4):392-408.

[6]李美香，張宏偉，李衛(wèi)國.基于點(diǎn)插值的配點(diǎn)型無網(wǎng)格方法解 Helmholtz問題[J].計算力學(xué)學(xué)報，2010，27（3）：533-536.

[7]姚凌云，于德介，臧獻(xiàn)國.聲學(xué)數(shù)值計算的分區(qū)光滑徑向點(diǎn)插值無網(wǎng)格法[J].振動與沖擊，2011，30(10)：188-192.

[8]胡海，郭文勇，馬龍.結(jié)構(gòu)聲輻射計算的邊界無網(wǎng)格法[J].噪聲與振動控制，2012，32(5)：37-41.

[9]LI K,HUANG Q B,WANG J L,LIN L G.An improved localized radialbasisfunction meshlessmethod for computational aeroacoustics[J].Eng.Anal.Bound Elem.,2011(35):47-55.

[10]曾向陽，王海濤，杜博凱.基于聲粒子分布積分的無網(wǎng)格聲場計算方法[J].應(yīng)用聲學(xué)，2016，5（1）：84-89.

[11]FANG Z,LIU C Y.Combined mesh free method and mode matching approach for transmission loss predictions of expansion chamber silencers[J].Eng.Anal.Bound Elem.,2017(84):168-177.

[12]I Z L,FANG Z.On the acoustic impedance of perforates and its application to acoustic attenuation predictions for perforated tube silencers[C].Inter-noise 2011,Osaka,Japan,September,2011.

[13]方智，季振林.穿孔管阻性消聲器橫截面模態(tài)及消聲特性計算與分析[J].振動與沖擊，2014，33(7)：138-146.