全局弱式無網格方法求解消聲器橫向模態

趙 輝，左繼強，劉成洋，方 智

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

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

本文基于無網格方法本身的優勢，使用徑向基函數點插值法離散本征方程，應用伽遼金加權殘數法進行數值積分，構成全局弱式無網格方法用來求解消聲器橫向本征波數以及模態形狀，更方便添加導數邊界條件。進而分析研究形狀參數對計算精度和速度的影響。

1 全局弱式無網格方法求解橫向模態

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

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

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

其中：為二維直角坐標系下的拉普拉斯算子，pxy1和pxy2分別為區域C1和C2內的橫向聲壓分量，kxy和分別表示空氣和吸聲材料中的橫向波數，并且滿足以下關系

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

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

穿孔面邊界條件

其中：ρ0為空氣密度為多孔吸聲材料密度為一側有吸聲材料時穿孔結構的特性聲阻抗率

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

其中：各個參數的含義可參考文獻[12]。

1.1 無網格徑向基點插值法求解形函數

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

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

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

使方程滿足計算點x周圍的n個節點值，一個節點對應一個方程，將產生n個線性方程。求解以下方程即可得到系數矩陣a0。

其中，各個矩陣表達式如下

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

將系數矩陣代入到方程，可以得到

對應節點聲壓向量的RPIM形函數Φ表達式為

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

1.2 基于Galerkin弱式離散系統方程

應用加權殘數法構造橫向本征方程和邊界條件的弱式無網格形式。采用伽遼金加權殘數方法，權函數和形函數選取一致，如同式(21)所示。本征方程和邊界條件的殘數乘以權函數并且在各自域內積分，可以得到加權形式

應用格林公式，方程簡化為以下線性方程

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

在每一個背景網格單元中配置高斯點，對于每一個高斯節點形成支持域，求解每一個高斯點的形函數Φ，進而計算背景網格單元的剛度矩陣，質量矩陣以及阻抗矩陣。

其中：nQ為單元e內的高斯點個數，wQ為相應的高斯加權系數，Je為單元的雅克比矩陣。

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

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

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

1.3 求解本征方程

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

橫向聲壓分量可以表示為

最終，可以得到本征函數的表達式為

2 計算精度與效率驗證

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

2.1 圓形截面

因為規則橫截面有橫向本征波數的解析解，所以如圖2所示的圓形截面最先被分析，圓形截面半徑為r=0.0245 m。截面內的場節點的分布是通過簡單的MATLAB程序構造的，背景網格由商業軟件ANSYS得到。由圖2可以看出，場節點是任意分布的，與背景網格是相互獨立的。

圖2圓形橫截面內任意分布的場節點和背景網格

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

圖3 圓形截面本征值計算結果比較

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

2.2 非規則截面

為了證明本文方法的適用性和準確性，本小節計算了如圖4所示的非規則截面的本征值。四邊形的具體尺寸如圖上所標，長度單位為m。

圖4 非規則四邊形橫截面

圖5 不規則四邊形截面本征值計算結果比較

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

2.3 圓形穿孔截面

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

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

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

由公式可知，每一個激發頻率都有相應的本征值，即軸向本征波數。本文計算了激發頻率為3000 Hz時前5階徑向模態的軸向本征波數，表1中給出了分別使用解析方法，二維有限元方法和本文無網格方法計算的本征值，可以看出達到一定的精度，有限元方法需要1600節點構成八節點四邊形單元，需要時間為445 s，而無網格方法只需140個節點，需要耗時75 s即可達到相同的精度。

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

通過以上3個算例可以看出，相對于有限元方法，無網格方法求解橫截面本征模態，不依賴于網格，需要的節點更少，耗時更少。橫截面面積越大，邊界條件越復雜，無網格方法的優勢就越突出。

3 計算結果與分析

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

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

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

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

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

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

4 結語

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

與有限元方法相比，無網格方法中背景網格不依賴于節點，在計算過程中可以隨意添加會刪除某些節點，這也是其優勢所在。

[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]李美香，張宏偉，李衛國.基于點插值的配點型無網格方法解 Helmholtz問題[J].計算力學學報，2010，27（3）：533-536.

[7]姚凌云，于德介，臧獻國.聲學數值計算的分區光滑徑向點插值無網格法[J].振動與沖擊，2011，30(10)：188-192.

[8]胡海，郭文勇，馬龍.結構聲輻射計算的邊界無網格法[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]曾向陽，王海濤，杜博凱.基于聲粒子分布積分的無網格聲場計算方法[J].應用聲學，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]方智，季振林.穿孔管阻性消聲器橫截面模態及消聲特性計算與分析[J].振動與沖擊，2014，33(7)：138-146.