三維全局弱式無網(wǎng)格方法計算膨脹腔消聲器聲學模態(tài)

劉鳳景，方春慧

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

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

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

1 聲腔模態(tài)

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

圖1 簡單膨脹腔消聲器

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

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

應用格林公式，可得

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

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

2 基于全局弱式無網(wǎng)格方法求解三維聲學模態(tài)

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

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

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

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

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

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

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

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

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

對應節(jié)點聲壓向量的徑向點插值形函數(shù)Φ表達式為

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

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

文中背景網(wǎng)格采用四面體單元，采用Hammer積分求解系統(tǒng)矩陣。每個單元中取5個Hammer積分點，對每一個Hammer積分點構(gòu)造支持域，求解每一個Hammer積分點的形函數(shù)Φ，進而計算背景網(wǎng)格單元的剛度矩陣和質(zhì)量矩陣。系統(tǒng)矩陣的求解方程如下。

其中：nQ為單元e內(nèi)的Hammer積分點個數(shù)，wQ為相應的加權(quán)系數(shù)，Je為單元的雅克比矩陣。

2.3 求解三維聲學模態(tài)

假設(shè)腔內(nèi)的節(jié)點數(shù)為N，通過式可以求得膨脹腔內(nèi)波數(shù)向量（即模態(tài)波數(shù)向量）k和特征向量矩陣（即模態(tài)振型矩陣）Xxyz，其中k的維數(shù)為N，Xxyz的階數(shù)為N×N。

3 計算精度驗證與結(jié)果分析

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

圖2 膨脹腔背景單元網(wǎng)格圖

圖3 膨脹腔場節(jié)點分布圖

場節(jié)點的分布是通過簡單的MATLAB程序構(gòu)造的，背景網(wǎng)格根據(jù)商業(yè)軟件ANSYS得到。由圖3可以看出，場節(jié)點是任意分布的，與背景網(wǎng)格是相互獨立的。

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

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

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

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

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

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

圖4 膨脹腔三維聲學模態(tài)振型圖

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

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

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

圖6 側(cè)支出口消聲器示意圖

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

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

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

4 結(jié)語

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