宣領(lǐng)寬，張文平，明平劍，龔京風(fēng)，許萬軍

(哈爾濱工程大學(xué)動力與能源工程學(xué)院，黑龍江哈爾濱150001)

對于其他節(jié)點2、3、4，同樣可以得到

在消聲器的設(shè)計階段多采用傳遞損失作為衡量消聲器聲學(xué)性能的指標，因此，能夠準確的計算出消聲器的傳遞損失對消聲器優(yōu)化設(shè)計有很大的實用價值.目前，用于預(yù)測消聲器傳遞損失的方法包括頻域法與時域法.頻域法主要包括傳遞矩陣法［1］、解析法［2］、有限元法(FEM)［3］及邊界元法(BEM)［4］，特別是一些基于FEM及BEM的商業(yè)軟件已經(jīng)在行業(yè)內(nèi)廣泛應(yīng)用;時域法近年來得到快速的發(fā)展，可分為一維時域法與三維時域法，一維時域法主要采用有限差分法(FDM)［5-10］;三維時域法又稱為消聲器聲學(xué)性能的計算流體動力學(xué)(CFD)方法，目前主要的研究方法是基于FLUENT等商業(yè)軟件預(yù)測消聲器的聲學(xué)性能［11-13］.

時域法的優(yōu)勢在于易于執(zhí)行，一次計算可獲得整個頻段內(nèi)的信息，并且對計算機的內(nèi)存依賴性不像頻域法那樣高，同時，頻域法由于要計算復(fù)雜的矩陣，所以需要大量的內(nèi)存，因此頻域法難以計算大型的、復(fù)雜的模型，而時域法則有廣闊的應(yīng)用前景，完善這一數(shù)值模擬技術(shù)的主要問題是在保證結(jié)果可靠的基礎(chǔ)上提高計算效率;就算法而言，提高計算效率的重要途徑是發(fā)展顯式方法［14］.

FVM是近年來發(fā)展非常迅速的一種離散化方法［15］，廣泛應(yīng)用于計算流體力學(xué)中，該方法能像FEM一樣適用于不規(guī)則網(wǎng)格和復(fù)雜邊界情況，且處理效率與FDM相似，遠遠高于FEM，所以在數(shù)值模擬中有著很大的發(fā)展?jié)摿?因此，采用交錯的非結(jié)構(gòu)FVM、顯格式推進求解線性聲波動方程，計算消聲器的聲學(xué)性能，網(wǎng)格采用交錯網(wǎng)格.這種交錯網(wǎng)格即為文獻［16］中提出的非結(jié)構(gòu)化網(wǎng)格，它與FEM的網(wǎng)格離散方法相同，但又區(qū)別于FEM，因為它不需要計算單元剛度矩陣，節(jié)點上的變量只用來計算單元中心的空間導(dǎo)數(shù);它與FDM交錯網(wǎng)格的積分方法相同.

因此針對小振幅擾動波動方程采用三維交錯非結(jié)構(gòu)FVM解法，并研究這種數(shù)值方法的色散關(guān)系，進行穩(wěn)定性分析，計算結(jié)果與商業(yè)軟件SYSNOISE有限元計算結(jié)果、文獻中的實驗結(jié)果進行對比.另外，該方法已經(jīng)在內(nèi)部代碼GTEA中實現(xiàn).

1 聲學(xué)數(shù)學(xué)模型

1.1 控制方程

靜止介質(zhì)中流體運動基本方程適用于聲波，忽略移流加速度項，可得到線性聲波動方程［17］:

式中:▽2為拉普拉斯算子，c0為聲速，p為聲壓.

1.2 數(shù)值方法

采用交錯的非結(jié)構(gòu)FVM求解式(1)，時間上應(yīng)用二階中心差分格式進行離散，空間上應(yīng)用有限體積法進行離散，算法采用時間顯格式推進求解.

式(1)左端的時間項采用二階中心差分格式進行離散:

式中:壓力p的上標表示時刻，下標表示壓力的位置，Δt表示時間步長，V表示積分控制單元的體積.

式(1)右端的空間項采用交錯FVM進行離散，對于空間項由高斯積分得

如圖1所示為交錯網(wǎng)格示意圖，梯度在四面體單元1234內(nèi)是均勻的，壓力存在節(jié)點上，四面體單元1234對節(jié)點1的貢獻可以表示為▽p·Ai，其中Ai為空間多邊形E12O4E13O3E14O2(由3個四邊形組成)的面積矢量.因此有

對于其他節(jié)點2、3、4，同樣可以得到

式中:n為節(jié)點1鄰近的四面體單元個數(shù)，▽pi為節(jié)點1鄰近的第i個單元內(nèi)梯度，任意的四面體單元1234內(nèi)的梯度有

其中對于節(jié)點1有

(x2，y2，z2)、(x3，y3，z3)與(x4，y4，z4)分別為節(jié)點2，3與4的坐標.

因此有

式中:V為控制體的體積，可由下式得到

式中:Vi為第i個四面體單元體積.

圖1 交錯網(wǎng)格示意Fig.1 Sketch map of staggered-grid

1.3 邊界條件

采用的邊界條件有2種:

1)“絕對硬”理想邊界，其法向速度v=0，即其在邊界上壓力的法向梯度為零，可表示為

數(shù)值計算中該邊界條件自動滿足，可以不用特殊處理，文中假設(shè)消聲器壁面結(jié)構(gòu)即為此邊界條件.

2)全透射邊界條件，作為對無限遠處的近似.采用一階Clayton-Engquist-Majda(C-E-M)吸收邊界條件［18-19］:

在邊界上差分離散有

式中:pb表示邊界面上的壓力，n表示時間層.認為每一個網(wǎng)格面上的壓力不變，等于面上各節(jié)點間的線性插值.

1.4 改進的脈沖法

采用文獻［13］中使用的脈沖法，并進行相應(yīng)的改進，使數(shù)據(jù)處理更加簡便，并且不需要限制出口管長度.以如圖2所示外插管消聲器為例，簡要說明計算消聲器傳遞損失的主要步驟.

圖2 外插管消聲器結(jié)構(gòu)示意Fig.2 Expansion chamber muffler with extended inlet and outlet

主要步驟為:

1)左側(cè)進口設(shè)定高斯壓力脈沖作為噪聲源p=exp［－(t－3T)2/T2］，記錄此脈沖信號做為入射波(入射波中即不會存在反射波)，當(dāng)完整的入射波給定后，進口邊界再設(shè)定為全透射邊界(可以去除進口反射波);

2)消聲器的右側(cè)出口設(shè)定透射邊界條件，監(jiān)測下游監(jiān)測點point1的透射壓力，消聲器上下游均能作為無窮遠的近似;

3)消聲器本身設(shè)定絕對硬壁面.

如圖3為消聲器上游的入射波及其頻譜，如圖4為消聲器下游的透射波及其頻譜.消聲器的傳遞損失［20］可表示為

式中:pin為消聲器進口處的入射聲壓，ptr末端為無反射條件下消聲器出口處的透射聲壓，Ai的A0分別為消聲器進口和出口的橫截面積.

圖3 入射波及其頻譜Fig.3 Incident wave and its spectrum

圖4 透射波及其頻譜Fig.4 Transmitted wave and its spectrum

2 數(shù)值色散及穩(wěn)定性分析

通過討論聲波動方程(1)有限差分方程的色散關(guān)系進行穩(wěn)定性分析.

如圖5所示為一正立方體被劃分為5個四面體，點(l，m，n)周圍圍繞有8個四面體單元，在圖5中有一個四面體單元ABCD.點A(l，m，n)的有限差分方程為

假設(shè)介質(zhì)中一平面波具有角頻率ω，則任意點的壓力可表示為

式中:k是波矢量，r是位置矢量，p0為聲壓幅值.將式(18)代入式(17)，由于p0非零，可以得到正三角形的色散關(guān)系如下:

當(dāng) Δx→0與 Δt→0時，式(19)可化為 ω2=++，這即為線性波動方程的理想色散關(guān)系.由于sin2(ωΔt/2)≤1.0，得到正三角形網(wǎng)格的穩(wěn)定性條件為

考慮一平面波波長為λ、波數(shù)為k，以一定角度入射，參考文獻［21］，取 γ =(3/2)1/2Δt/Δx及 H=Δx/λ，這樣就可以得到無量綱相速度q(數(shù)值相速度v與真實波速c0的比值)如下

如圖6所示為γ=0.8時，不同入射角度的色散關(guān)系.理論上q(H)應(yīng)不大于1，實際上它提供一種約束空間步長的規(guī)則，對于圖4中所示的幾個不同的角度 θ，當(dāng) H=0.09，q≈0.991，也就是說一個波長內(nèi)需要11個空間步長才能精確的描述聲輻射問題.但實際中很難做到所有的網(wǎng)格都如圖5中所考慮的那樣，為保證計算精度，最大的邊長應(yīng)該小于最小關(guān)心波長的1/11，將這個結(jié)果與式(20)聯(lián)合起來，就可以合理的選擇Δx和Δt.

圖5 立方體中四面體網(wǎng)格示意Fig.5 Sketch map of the tetrahedral grids in one cube

圖6 色散常數(shù)γ=0.8時的色散曲線Fig.6 Dispersion curves for non-dimensional acoustic wave phase velocity with a dispersion parameter γ=0.8

平面波的入射方向分別為:

3 抗性消聲器傳遞損失的數(shù)值應(yīng)用

以下計算結(jié)果均是在Intel Core2 CPU主頻為1.86 GHz，內(nèi)存為1.0 GB的計算機上完成的.

3.1 外插管消聲器

計算圖2所示外插管消聲器的傳遞損失，如圖7所示為其網(wǎng)格模型，聲速c0=346.093 m/s，時間步長為5.0 ×10－66 s，總計算時間為 0.2 s，消聲器進口給高斯壓力脈沖(T=10×10－4).

圖7 外插管消聲器網(wǎng)格模型Fig.7 Unstructured mesh of expansion chamber muffler with extended inlet and outlet

進行網(wǎng)格無關(guān)性分析，采用4種網(wǎng)格進行分析，采用四面體網(wǎng)格，網(wǎng)格的相關(guān)信息列于表1.Lmax為網(wǎng)格中最大邊長，fmax為最大有效頻率，它是通過下式得到

表1 網(wǎng)格無關(guān)性的相關(guān)信息Table 1 Details of grid independency

如圖8為不同網(wǎng)格的傳遞損失計算結(jié)果比較，網(wǎng)格3的結(jié)果與網(wǎng)格4結(jié)果吻合良好，除在共振頻率處存在一定的差別.因此將網(wǎng)格3的結(jié)果與實驗結(jié)果及SYSNOISE有限元計算結(jié)果(采用網(wǎng)格3)對比，如圖9所示.從計算結(jié)果可以看出，計算結(jié)果與文獻［22］中的實驗結(jié)果及SYSNOISE(FEM)計算結(jié)果在整個頻段內(nèi)吻合良好，但FVM方法、FEM均忽略介質(zhì)粘性，因此在共振頻率處預(yù)測結(jié)果與實驗測量結(jié)果間均存在明顯偏差;而在一些頻段內(nèi)，F(xiàn)VM方法與實驗測量結(jié)果更接近，表明FVM時域方法可以比較準確的預(yù)測無流情況下抗性消聲器的聲學(xué)性能.

圖8 不同網(wǎng)格的傳遞損失計算結(jié)果比較Fig.8 Comparison of transmission losses of different meshes

圖9 外插管消聲器傳遞損失比較Fig.9 Comparison of transmission loss of expansion chamber muffler

3.2 復(fù)雜蒸汽消聲器

圖10為某蒸汽管路消聲器結(jié)構(gòu)示意圖.消聲器為有5個腔的軸對稱結(jié)構(gòu)，氣流方向如圖10所示，氣流由進口進入消聲器后，由A腔經(jīng)中間軸向孔及周圍小孔進入B腔，由B腔經(jīng)靠近外壁上的小孔進入C腔，由C腔經(jīng)中間軸向孔進入D腔，由D腔經(jīng)靠近外腔壁上的孔進入E腔，最后氣流由E腔經(jīng)出口流出.聲速 c0=507.3 m/s，密度 ρ=11.92 kg/m3，進出口管的直徑D=0.097 m，則其截止頻率fcut-off=3 763 Hz.計算模型可以采用1/4結(jié)構(gòu)，采用四面體網(wǎng)格，網(wǎng)格數(shù)101 742(對進口管長度有要求)，SYSNOISE采用四面體網(wǎng)格，網(wǎng)格數(shù)為97 689.如圖11比較FVM計算結(jié)果與SYSNOISE(FEM)計算結(jié)果，兩者吻合良好.

圖12將該消聲器的出口管改到側(cè)面，則該消聲器結(jié)構(gòu)的進出口邊界由軸對稱變?yōu)橹挥幸粋€對稱面，這就導(dǎo)致計算網(wǎng)格會成倍的增加，同時為測試FVM方法在大型計算中的能力，劃分網(wǎng)格數(shù)為560 567，對于在本節(jié)提到的計算機上，F(xiàn)EM是難以完成的，而FVM方法由于采用顯格式推進，并且不需要計算大型的剛度矩陣，所以可以比較快捷的處理較大型計算問題.如圖13對比出口管位置不同時消聲器的傳遞損失.可以看出，當(dāng)出口改為側(cè)面后在共振頻率處出現(xiàn)多個峰值:第一個共振峰是由于側(cè)面出口管與消聲器端部之間形成一個聲學(xué)共振腔引起的;后面的幾個共振峰是由于側(cè)面出口使消聲器變?yōu)榉菍ΨQ結(jié)構(gòu)而激發(fā)周向模態(tài)，消聲器的高階模態(tài)提前激起而形成的，出口改到側(cè)面后消聲器的整體聲學(xué)性能得到改善，特別是1 500 Hz以下的聲學(xué)性能得到明顯提高.

圖10 軸向出口蒸汽消聲器結(jié)構(gòu)示意Fig.10 Sketch map of steam muffler with end-outlet

圖11 端部出口蒸汽消聲器傳遞損失比較Fig.11 Transmission loss of steam mufler

圖12 側(cè)面出口消聲器結(jié)構(gòu)示意Fig.12 Sketch map of steam muffler with side-outlet

圖13 端部出口與側(cè)面出口蒸汽消聲器傳遞損失比較Fig.13 Comparison of transmission loss of end-outlet and side-outlet

4 結(jié)論

基于三維交錯的非結(jié)構(gòu)FVM求解線性聲波動方程，綜合運用全反射與全透射邊界條件，采用時域脈沖法預(yù)測消聲器的傳遞損失.在消聲器的進口給定脈沖后，再設(shè)為全透射邊界條件，可以有效地消除進口處的反射波，使數(shù)據(jù)處理更加簡便，并且不需要限制出口管長度.詳細推導(dǎo)出三維四面體網(wǎng)格的穩(wěn)定性判據(jù)及色散關(guān)系.

利用三維FVM方法計算外插管消聲器的聲學(xué)性能，計算結(jié)果與商業(yè)軟件SYSNOISE計算結(jié)果、文獻中實驗結(jié)果吻合良好，表明該方法能夠比較準確地預(yù)測無流情況下抗性消聲器的聲學(xué)性能.通過分析復(fù)雜結(jié)構(gòu)消聲器的聲學(xué)性能，展示FVM時域方法處理大型模型的能力，同時，結(jié)果表明側(cè)面出口管相對端部出口管可以改善蒸汽消聲器的聲學(xué)性能.

［1］MUNJAL M L，RAO K N，SAHASRABUDHE A D.Aeroacoustic analysis of perforated muffler components［J］.Journal of Sound and Vibration，1987，114(2):173-188.

［2］MILES J.The reflection of sound due to change in cross section of circular tube［J］.Journal of the Acoustical Society of America，1944，16(1):14-19.

［3］PEAT K S.Evaluation of four-pole parameters for ductswith flow by the finite element method［J］.Journal of Sound and Vibration，1982，84(3):389-395.

［4］JIZ L，MA Q，ZHANG Z H.Application of the boundary element method to predicting acoustic performance of expansion chamber mufflers with mean flow［J］.Journal of Sound and Vibration，1994，173(1):57-71.

［5］CHANG IJ，CUMMINGSA.A time domain solution for the attenuation，at high amplitudes，of perforated tube silencers and comparison with experiment［J］.Journal of Sound and Vibration，1988，122:243-259.

［6］MOREL T，MOREL J，BLASER D A.Fluid dynamic and acoustic modeling of concentric-tube resonators/silencers［J］.SAE Transactions，1991，100:35-51.

［7］SELAMET A，DICKEY N S，NOVAK JM.A time-domain computational simulation of acoustic silencers［J］.Journal of Vibration and Acoustics，1995，223:323-331.

［8］ONORATIA，PEROTTIM，REBAY S.Modelling one dimensional unsteady flows in ducts:symmetric finite difference schemes versus Galerkin discontinuous finite element methods［J］.International Journal of Mechanical Sciences，1997，39(11):1213-1236.

［9］DICKEY N S，SELAMET A.Effects of numerical dissipation and dispersion on acoustic predictions from a time-domain finite difference technique for non-linear wave dynamics［J］.Journal of Sound and Vibration，2003，259(1):193-208.

［10］BROATCH A，SERRANO JR，MOYA D.Time-domain computation ofmuffler frequency response:comparison of different numerical schemes［J］.Journal of Sound and Vibration，2007，305:333-347.

［11］MIDDELBERG JM，BARBER T J，LEONG SS.Computational fluid dynamics analysis of the acoustic performance of various simple expansion chamber mufflers［C］//Proceedings of Acoustics 2004.Gold Coast，2004:123-127.

［12］BROATCH A，MARGOT X，GIL A.A CFD approach to the computation of the acoustic response of exhaust mufflers［J］.Journal of Computational Acoustics，2005，13(2):301-316.

［13］徐航手，季振林，康鐘緒.抗性消聲器傳遞損失的三維時域計算方法［J］.振動與沖擊，2010，29(4):107-110.XU Hangshou，JI Zhenlin，KANG Zhongxu.Three dimensional time domain computational approach for predicting transmission loss of reactive silencers［J］.Journal Vibration and Shock，2010，29(4):107-110.

［14］劉恒，廖振鵬.波動數(shù)值模擬的一種顯式方法——二維波動［J］.力學(xué)學(xué)報，2010，42(6):1104-1116.LIU Heng，LIAO Zhenpeng.An explicit method for numerical simulation of wave motion——2-D wave motion［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(6):1104-1116.

［15］王福軍.計算流體動力學(xué)分析——CFD軟件原理與應(yīng)用［M］.北京:清華大學(xué)出版社，2004:1-2，25-28.

［16］LIU T L，LIU K S，ZHANG JX.Unstructured grid method for stress wave propagation in elastic media［J］.Computer Methods in Applied Mechanics and Engineering，2004，193:2427-2453.

［17］杜功煥，朱哲民，龔秀芬.聲學(xué)基礎(chǔ)［M］.南京:南京大學(xué)出版社，2001:176-177.

［18］李太寶.聲場的方程和計算方法［M］.北京:科學(xué)出版社，2005:59-61.

［19］居鴻賓，沈孟育.氣動聲學(xué)數(shù)值模擬中無反射邊界處理及聲源模型的建立［J］.上海交通大學(xué)學(xué)報，1998，32(7):36-39.JU Hongbin，SHEN Mengyu.Non-reflective boundary conditions and sound source modeling in computational aeroacoustic［J］.Journal of Shanghai Jiao Tong University，1998，32(7):36-39.

［20］MUNJALM L.Acoustics of ducts and mufflers［M］.New York:Wiley Inter Science，1987:92-95.

［21］VIRIEUX J.P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method［J］.Geophysics，1986，51:889-901.

［22］SELAMET A，JIZ L.Acoustic attenuation performance of circular expansion chambers with extended inlet/outlet［J］.Journal of Sound and Vibration，1999，223(2):197-212.