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

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

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

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

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

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

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

以膨脹腔C為例，介紹全局弱式無(wú)網(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è)有吸聲材料時(shí)穿孔結(jié)構(gòu)的特性聲阻抗率

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 求解本征方程

假設(shè)橫截面上的節(jié)點(diǎn)數(shù)為n，通過(guò)式和可以分別求得膨脹腔和進(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 計(jì)算精度與效率驗(yàn)證

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

2.1 圓形截面

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

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

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

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

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

2.2 非規(guī)則截面

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

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

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

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

2.3 圓形穿孔截面

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

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

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

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

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

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

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

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

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

圖7 支持域尺寸對(duì)計(jì)算精度的影響

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

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

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

4 結(jié)語(yǔ)

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

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

[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)國(guó).基于點(diǎn)插值的配點(diǎn)型無(wú)網(wǎng)格方法解 Helmholtz問(wèn)題[J].計(jì)算力學(xué)學(xué)報(bào)，2010，27（3）：533-536.

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

[8]胡海，郭文勇，馬龍.結(jié)構(gòu)聲輻射計(jì)算的邊界無(wú)網(wǎng)格法[J].噪聲與振動(dòng)控制，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]曾向陽(yáng)，王海濤，杜博凱.基于聲粒子分布積分的無(wú)網(wǎng)格聲場(chǎng)計(jì)算方法[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ì)算與分析[J].振動(dòng)與沖擊，2014，33(7)：138-146.