三維全局弱式無網格方法計算膨脹腔消聲器聲學模態

劉鳳景，方春慧

(煙臺汽車工程職業學院，山東 煙臺 265500）

消聲器在汽車發動機進排氣系統中應用較多。消聲器的聲學模態對消聲器的消聲特性有著較大影響，可以根據聲學模態特性改進消聲器結構，進而優化消聲性能。聲學模態的計算通常有解析方法和數值方法。文獻[1] 采用解析方法計算了圓形截面簡單膨脹腔消聲器的二維橫向模態，分析了消聲器進出口位置對消聲性能的影響。解析方法具有簡單快速的優點，但是只能應用于規則形狀截面的模態計算。文獻[2－3] 使用有限元方法計算了膨脹腔消聲器和穿孔管消聲器的二維橫向模態，進而分析了模態對消聲器消聲性能的影響。有限元方法可以計算任意形狀截面的橫向模態，應用范圍較廣，只是計算精度與截面單元的數目和質量相關?；诖?，無網格方法被提出用來計算消聲器的橫向模態[4－5]。無網格方法(MFM)不需要劃分網格，利用一組散布在問題域中及邊界上的節點表示該問題。無網格方法根據形函數的求解和系統方程的離散可以組合出很多方法。Chen等[6]應用基于徑向基函數的邊界配點無網格方法計算了二維聲腔的聲學本征值，文章驗證了無網格方法與有限元方法相比，具有不需要構造復雜的背景網格，容易實施的特點。張宏偉等[7]基于點插值的思想，用三角函數作為基函數構造局部近似函數，用配點法離散微分方程，提出了點插值配點型無網格法，研究了Helmholtz方程的波傳播問題和邊界層問題，通過算例比較證明了隨著節點的增加，該方法具有較高的精度和良好的收斂性。姚凌云等[8]推導了分區光滑徑向點插值無網格法（CSRPIM），將其應用到管道和二維轎車的聲學問題分析中，數值計算結果表明該方法比標準有限元法具有更高的精度和準確度。黃其柏等[9]使用改進的局部徑向配點無網格法研究了氣動聲學和腔體結構聲學特性，數值計算結果與實驗結果對比表明，徑向基函數無網格方法是一種真正的無網格方法，擺脫了對網格的依賴，具有精度高、編程簡單等優點。2016年，曾向陽等[10]提出一種基于聲粒子分布積分的無網格方法，計算矩形封閉空間中低頻聲場，計算結果與采用商用計算軟件以及經典無網格法的結果的對比，證明了該方法在中低頻段與傳統無網格方法相比具有更高的精度。2017年，方智等[11]應用強式無網格配點方法求解任意橫截面膨脹腔消聲器的橫向模態，使用模態匹配技術求解傳遞損失，提出一種基于無網格的數值模態匹配方法，計算結果與實驗測量值吻合較好，且與三維有限元方法相比，計算速度具有一定的優勢。但當消聲器壁面導數邊界條件較多時，使用強式無網格配點方法穩定性會變差。于是，基于伽遼金加權殘數法的全局弱式無網格方法[5]被提出用于計算消聲器的橫向模態，與理論解吻合較好。

以上所有對消聲器模態的計算均是二維的，只是計算了消聲器的橫向模態，并未分析三維模態效應。本文應用三維全局弱式無網格方法計算簡單膨脹腔消聲器的三維聲學模態，使用徑向基函數點插值法求解三維無網格形函數，應用伽遼金加權殘數法對系統方程進行離散，采用三維全局弱式無網格方法用來求解消聲器的三維聲學模態頻率和模態振型，進而分析軸向模態對消聲器聲學性能的影響。

1 聲腔模態

圖1為簡單膨脹腔消聲器示意圖。膨脹腔內的三維聲壓控制方程為

圖1 簡單膨脹腔消聲器

膨脹腔的邊界條件為剛性壁邊界條件

在聲場無網格計算中，選擇形函數來表示聲壓，如方程式(3)所示

應用格林公式，可得

將邊界條件表達式(2)代入上式中，方程式(5)簡化為以下線性方程

其中K、M分別為腔內的廣義剛度矩陣和質量矩陣。為了求解這些系統矩陣，需要求解域內和邊界上的積分。因此在弱式無網格方法中需要構造背景網格。背景網格與場節點可以是相互獨立的，并且對網格質量的要求不太高[5]。

2 基于全局弱式無網格方法求解三維聲學模態

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

文中采用徑向基點插值法（RPIM）求解形函數。添加多項式的RPIM插值可以表示為

其中：Ri(x)為徑向基函數，n為徑向基函數的個數，qk(x)為空間坐標x中的單項式，m為多項式基函數的個數。方程式(7)中每一項的表達式如下所示

構造計算點x的支持域來確定系數矩陣中的ai和ck，其中支持域中的場節點數為n。支持域的大小由ds=αsdc確定，αs為支持域的無量綱尺寸，在本文中取αs=2.5，dc為計算點x支持域內的平均節點間距。在徑向基函數Ri(rk)中，αc為形狀參數，在文中取αc=2.5。第i個場節點和第k個場節點之間的距離rk為

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

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

集果裝置由傘形集果筒、搭扣、合頁等組成，傘形集果筒共有4個落果端口，位于集果筒底部，通過對紅棗樹枝幾何拓撲結構的統計分析[12]，確定集果筒高度為850 mm，頂面直徑1 200 mm，底面直徑240 mm。設計傘形集果筒如圖5所示。

由于矩陣G是對稱矩陣，因此可由方程式(12)求得系數矩陣

將方程式(16)代入到方程式(7)，可以得到

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

最終，聲壓p(x)可以表示為形函數與節點聲壓向量的乘積形式

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

文中背景網格采用四面體單元，采用Hammer積分求解系統矩陣。每個單元中取5個Hammer積分點，對每一個Hammer積分點構造支持域，求解每一個Hammer積分點的形函數Φ，進而計算背景網格單元的剛度矩陣和質量矩陣。系統矩陣的求解方程如下。

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

2.3 求解三維聲學模態

假設腔內的節點數為N，通過式可以求得膨脹腔內波數向量（即模態波數向量）k和特征向量矩陣（即模態振型矩陣）Xxyz，其中k的維數為N，Xxyz的階數為N×N。

3 計算精度驗證與結果分析

為了驗證本文方法的計算精度，文中計算了簡單膨脹腔消聲器的前23階三維聲學模態頻率和模態形狀。計算結果與有限元結果進行對比。圖1為簡單膨脹腔結構圖，具體尺寸為：膨脹腔直徑D=0.1532m，膨脹腔長度l=0.2823m，進出口管直徑d=0.0486m。由于與膨脹腔體相比，進出管的直徑比較小，所以在較高的頻率范圍（0～4150 Hz）內均為平面波傳播，在所考察的頻率0～3200 Hz范圍內，消聲器的高階模態僅在膨脹腔內被激發。鑒于此只計算了膨脹腔的三維聲學模態。圖2為膨脹腔背景單元網格圖，圖3為膨脹腔場節點分布圖。

圖2 膨脹腔背景單元網格圖

圖3 膨脹腔場節點分布圖

場節點的分布是通過簡單的MATLAB程序構造的，背景網格根據商業軟件ANSYS得到。由圖3可以看出，場節點是任意分布的，與背景網格是相互獨立的。

表1給出了分別使用全局弱式無網格方法和有限元方法計算的前23階模態頻率。圖4列出了前23階高階模態中所有的模態振型。在模態振型中，表征聲壓為零的線稱為節線。三維聲學模態可以用周向節線數m、徑向節線數n和軸向節線數z共同表達，即（m，n，z）。

由表1可以看出與有限元結果對比，前23階高階模態頻率誤差均在1%以內，說明采用全局弱式無網格方法計算膨脹腔三維聲學模態的適用性以及代碼的正確性。

為了比較所用的無網格程序與三維有限元方法的計算速度，分別記錄了計算圖1結構消聲器前23階模態所需的時間。使用無網格方法時節點數取為1775，背景網格數為2505，單元類型為四面體單元。計算時間為126 s。三維有限元方法中使用ANSYS劃分三維網格，使用Virtual.lab計算聲學模態，節點數取1775，網格數為8654，單元類型為1階四面體單元，計算時間為36 s。由此可以看出與三維有限元商業軟件相比，本文編寫的無網格計算程序在計算速度上并沒有優勢，但是若計算的模型較為復雜，使用三維有限元方法時需要耗費精力劃分精度較高的網格，而無網格方法的計算節點由MATLAB生成，對網格的要求也不高，其計算優勢則會顯示出來。此外，使用無網格計算程序計算模態時，計算的模態階數與計算節點的數目相關，在此算例中，計算的模態階數為1775階，而使用Virtual.lab時計算的階數為23階，所以用無網格方法計算的時間較長也是有一定原因的。由此也說明本文使用的無網格計算程序有很大優化空間。

對照表1的模態頻率和圖4列出的18個模態振型圖可以發現，由于膨脹腔具有軸對稱特性，所有周向模態即m不等于0的模態均為重根模態，即模態頻率相等，模態振型關于中心軸對稱。對于重根模態在圖3中僅列出了一個模態振型。

聲學模態是消聲器固有特性，可以根據模態特點分析消聲器的消聲特性。圖5給出了應用聲學有限元軟件計算的膨脹腔消聲器的傳遞損失曲線。為了便于分析，在圖5中也給出了應用一維平面波理論計算的消聲器的傳遞損失，在該計算中考慮了進出口管的端部修正[3]。

由傳遞損失曲線可以看出，兩條曲線在（0，1，0）階模態頻率處出現明顯偏離，三維有限元計算結果曲線不再具有拱形特性，說明膨脹腔內出現了橫向高階模態，由于一維平面波理論沒有考慮橫向高階模態效應，所以計算結果不再精確。而在低于（0，0，1）階模態頻率范圍內，2條曲線基本吻合，這是因為在平面波理論中考慮了軸向的高維模態效應，由于同時考慮了進出口的端部修正，所以在一維平面波理論中將消聲器進出口截面處由于截面突變引起的可耗散的三維波效應也考慮了，所以兩條曲線基本吻合。在拱形處出現的偏差可能是由于端部修正并沒有將可耗散的三維波效應完全考慮所致。

圖4 膨脹腔三維聲學模態振型圖

圖5 同軸膨脹腔消聲器傳遞損失曲線

由于周向模態關于中心軸對稱，消聲器進出口管又位于中心軸上，所以周向模態均沒有被激發，而最先被激發的模態為第1階徑向模態，可以設置出口管位于該徑向模態節線上，拓寬消聲器的消聲范圍，這已在有關文獻中分析過[1,3]。通過傳遞損失曲線還可以發現消聲器的通過頻率對應著膨脹腔的軸向高階模態頻率，所以為了消除通過頻率，除了采用插管以外[3]，還可以將出口管布置在軸向模態節線處，采用側支出口形式。根據圖3中給出的模態振型圖可以看出，（0，0，1）、（0，0，2）階模態的節線分別位于膨脹腔長度的1/2、1/4和3/4處，可以設想將出口管分別設置在這些位置，可以將相應的軸向模態效應消除，即可消除第1個和第2個通過頻率。

文中分別將出口管布置在膨脹腔長度的1/2和3/4處，計算傳遞損失。側支出口消聲器示意圖如圖6所示。

圖6 側支出口消聲器示意圖

分 別 取δ=L∕2=0.14115m 和0.211725m，其他尺寸與圖1相同。應用三維有限元方法計算得到的傳遞損失曲線如圖7所示。

圖7 側支出口膨脹腔消聲器傳遞損失曲線

從圖7可以看出，由于消聲器進口管位于中心軸線上，所以周向模態無法被激發，所以在（1，0，0）、（2，0，0）這2個周向模態頻率處，兩種出口側支結構消聲器的傳遞損失曲線均出現了峰值。側支出口管距進口管偏移量為L/2的消聲器結構，出口管位于（0，0，1）模態節線上，所以在傳遞損失曲線中在該模態頻率處出現了峰值，與進出口管同軸結構消聲器相比，消除了第1個通過頻率。觀察曲線還可發現，該結構消聲器同時消除了第3階軸向模態（0，0，3）模態，這是由于（0，0，3）階模態中一個模態節線也位于L/2處。不僅如此，該結構消聲器還在（1，0，2）和（2，0，2）階模態頻率處展現了較好的消聲性能，在該頻率處傳遞損失出現峰值。同理，側支出口管距進口管偏移量為3L/4的消聲器結構，出口管位于（0，0，2）模態節線上，消除了（0，0，2）階模態的影響，分別在（0，0，2）、（1，0，1）和（2，0，1）模態頻率處展現較好的消聲性能。與進出口同軸消聲器相比，側支出口的消聲器結構，若把出口管布置在某階軸向模態的節線處，可以消除該階模態的影響，大大優化消聲性能。

4 結語

采用三維全局弱式無網格方法計算消聲器的三維聲學模態頻率和模態振型，使用徑向基函數點插值方法求解三維無網格形函數，使用伽遼金加權殘數法離散系統方程，最終求得消聲器的三維聲學模態。計算了一簡單膨脹腔的前23階三維聲學模態頻率，與有限元計算結果對比，相對誤差控制在1%以內，驗證了采用三維無網格方法計算聲學模態的正確性。進而根據聲學模態的特性，分析了軸向模態對消聲器消聲性能的影響。將出口管側支布置在軸向模態的節線處，可以消除相對應的通過頻率，進而優化消聲性能。