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

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

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

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

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

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

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

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

1.1 控制方程

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

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

1.2 數(shù)值方法

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

式(1)左端的時(shí)間項(xiàng)采用二階中心差分格式進(jìn)行離散:

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

式(1)右端的空間項(xiàng)采用交錯(cuò)FVM進(jìn)行離散，對(duì)于空間項(xiàng)由高斯積分得

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

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

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

其中對(duì)于節(jié)點(diǎn)1有

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

因此有

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

式中:Vi為第i個(gè)四面體單元體積.

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

1.3 邊界條件

采用的邊界條件有2種:

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

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

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

在邊界上差分離散有

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

1.4 改進(jìn)的脈沖法

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

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

主要步驟為:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 外插管消聲器

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

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

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

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

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

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

利用三維FVM方法計(jì)算外插管消聲器的聲學(xué)性能，計(jì)算結(jié)果與商業(yè)軟件SYSNOISE計(jì)算結(jié)果、文獻(xiàn)中實(shí)驗(yàn)結(jié)果吻合良好，表明該方法能夠比較準(zhǔn)確地預(yù)測無流情況下抗性消聲器的聲學(xué)性能.通過分析復(fù)雜結(jié)構(gòu)消聲器的聲學(xué)性能，展示FVM時(shí)域方法處理大型模型的能力，同時(shí)，結(jié)果表明側(cè)面出口管相對(duì)端部出口管可以改善蒸汽消聲器的聲學(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］徐航手，季振林，康鐘緒.抗性消聲器傳遞損失的三維時(shí)域計(jì)算方法［J］.振動(dòng)與沖擊，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］劉恒，廖振鵬.波動(dòng)數(shù)值模擬的一種顯式方法——二維波動(dòng)［J］.力學(xué)學(xué)報(bào)，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］王福軍.計(jì)算流體動(dòng)力學(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］李太寶.聲場的方程和計(jì)算方法［M］.北京:科學(xué)出版社，2005:59-61.

［19］居鴻賓，沈孟育.氣動(dòng)聲學(xué)數(shù)值模擬中無反射邊界處理及聲源模型的建立［J］.上海交通大學(xué)學(xué)報(bào)，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.