预测消声器声学性能的时域非结构有限体积法

1 声学数学模型

1.1 控制方程

静止介质中流体运动基本方程适用于声波，忽略移流加速度项，可得到线性声波动方程［17］:

式中:▽2为拉普拉斯算子，c0为声速，p为声压.

1.2 数值方法

采用交错的非结构FVM求解式(1)，时间上应用二阶中心差分格式进行离散，空间上应用有限体积法进行离散，算法采用时间显格式推进求解.

式(1)左端的时间项采用二阶中心差分格式进行离散:

式中:压力p的上标表示时刻，下标表示压力的位置，Δt表示时间步长，V表示积分控制单元的体积.

式(1)右端的空间项采用交错FVM进行离散，对于空间项由高斯积分得

如图1所示为交错网格示意图，梯度在四面体单元1234内是均匀的，压力存在节点上，四面体单元1234对节点1的贡献可以表示为▽p·Ai，其中Ai为空间多边形E12O4E13O3E14O2(由3个四边形组成)的面积矢量.因此有

对于其他节点2、3、4，同样可以得到

式中:n为节点1邻近的四面体单元个数，▽pi为节点1邻近的第i个单元内梯度，任意的四面体单元1234内的梯度有

其中对于节点1有

(x2，y2，z2)、(x3，y3，z3)与(x4，y4，z4)分别为节点2，3与4的坐标.

因此有

式中:V为控制体的体积，可由下式得到

式中:Vi为第i个四面体单元体积.

图1 交错网格示意Fig.1 Sketch map of staggered-grid

1.3 边界条件

采用的边界条件有2种:

1)“绝对硬”理想边界，其法向速度v=0，即其在边界上压力的法向梯度为零，可表示为

数值计算中该边界条件自动满足，可以不用特殊处理，文中假设消声器壁面结构即为此边界条件.

2)全透射边界条件，作为对无限远处的近似.采用一阶Clayton-Engquist-Majda(C-E-M)吸收边界条件［18-19］:

在边界上差分离散有

式中:pb表示边界面上的压力，n表示时间层.认为每一个网格面上的压力不变，等于面上各节点间的线性插值.

1.4 改进的脉冲法

采用文献［13］中使用的脉冲法，并进行相应的改进，使数据处理更加简便，并且不需要限制出口管长度.以如图2所示外插管消声器为例，简要说明计算消声器传递损失的主要步骤.

图2 外插管消声器结构示意Fig.2 Expansion chamber muffler with extended inlet and outlet

主要步骤为:

1)左侧进口设定高斯压力脉冲作为噪声源p=exp［－(t－3T)2/T2］，记录此脉冲信号做为入射波(入射波中即不会存在反射波)，当完整的入射波给定后，进口边界再设定为全透射边界(可以去除进口反射波);

2)消声器的右侧出口设定透射边界条件，监测下游监测点point1的透射压力，消声器上下游均能作为无穷远的近似;

3)消声器本身设定绝对硬壁面.

如图3为消声器上游的入射波及其频谱，如图4为消声器下游的透射波及其频谱.消声器的传递损失［20］可表示为

式中:pin为消声器进口处的入射声压，ptr末端为无反射条件下消声器出口处的透射声压，Ai的A0分别为消声器进口和出口的横截面积.

图3 入射波及其频谱Fig.3 Incident wave and its spectrum

图4 透射波及其频谱Fig.4 Transmitted wave and its spectrum

2 数值色散及稳定性分析

通过讨论声波动方程(1)有限差分方程的色散关系进行稳定性分析.

如图5所示为一正立方体被划分为5个四面体，点(l，m，n)周围围绕有8个四面体单元，在图5中有一个四面体单元ABCD.点A(l，m，n)的有限差分方程为

假设介质中一平面波具有角频率ω，则任意点的压力可表示为

式中:k是波矢量，r是位置矢量，p0为声压幅值.将式(18)代入式(17)，由于p0非零，可以得到正三角形的色散关系如下:

当 Δx→0与 Δt→0时，式(19)可化为 ω2=++，这即为线性波动方程的理想色散关系.由于sin2(ωΔt/2)≤1.0，得到正三角形网格的稳定性条件为

考虑一平面波波长为λ、波数为k，以一定角度入射，参考文献［21］，取 γ =(3/2)1/2Δt/Δx及 H=Δx/λ，这样就可以得到无量纲相速度q(数值相速度v与真实波速c0的比值)如下

如图6所示为γ=0.8时，不同入射角度的色散关系.理论上q(H)应不大于1，实际上它提供一种约束空间步长的规则，对于图4中所示的几个不同的角度 θ，当 H=0.09，q≈0.991，也就是说一个波长内需要11个空间步长才能精确的描述声辐射问题.但实际中很难做到所有的网格都如图5中所考虑的那样，为保证计算精度，最大的边长应该小于最小关心波长的1/11，将这个结果与式(20)联合起来，就可以合理的选择Δx和Δt.

图5 立方体中四面体网格示意Fig.5 Sketch map of the tetrahedral grids in one cube

图6 色散常数γ=0.8时的色散曲线Fig.6 Dispersion curves for non-dimensional acoustic wave phase velocity with a dispersion parameter γ=0.8

平面波的入射方向分别为:

3 抗性消声器传递损失的数值应用

以下计算结果均是在Intel Core2 CPU主频为1.86 GHz，内存为1.0 GB的计算机上完成的.

3.1 外插管消声器

计算图2所示外插管消声器的传递损失，如图7所示为其网格模型，声速c0=346.093 m/s，时间步长为5.0 ×10－66 s，总计算时间为 0.2 s，消声器进口给高斯压力脉冲(T=10×10－4).

图7 外插管消声器网格模型Fig.7 Unstructured mesh of expansion chamber muffler with extended inlet and outlet

进行网格无关性分析，采用4种网格进行分析，采用四面体网格，网格的相关信息列于表1.Lmax为网格中最大边长，fmax为最大有效频率，它是通过下式得到

表1 网格无关性的相关信息Table 1 Details of grid independency

如图8为不同网格的传递损失计算结果比较，网格3的结果与网格4结果吻合良好，除在共振频率处存在一定的差别.因此将网格3的结果与实验结果及SYSNOISE有限元计算结果(采用网格3)对比，如图9所示.从计算结果可以看出，计算结果与文献［22］中的实验结果及SYSNOISE(FEM)计算结果在整个频段内吻合良好，但FVM方法、FEM均忽略介质粘性，因此在共振频率处预测结果与实验测量结果间均存在明显偏差;而在一些频段内，FVM方法与实验测量结果更接近，表明FVM时域方法可以比较准确的预测无流情况下抗性消声器的声学性能.

图8 不同网格的传递损失计算结果比较Fig.8 Comparison of transmission losses of different meshes

图9 外插管消声器传递损失比较Fig.9 Comparison of transmission loss of expansion chamber muffler

3.2 复杂蒸汽消声器

图10为某蒸汽管路消声器结构示意图.消声器为有5个腔的轴对称结构，气流方向如图10所示，气流由进口进入消声器后，由A腔经中间轴向孔及周围小孔进入B腔，由B腔经靠近外壁上的小孔进入C腔，由C腔经中间轴向孔进入D腔，由D腔经靠近外腔壁上的孔进入E腔，最后气流由E腔经出口流出.声速 c0=507.3 m/s，密度 ρ=11.92 kg/m3，进出口管的直径D=0.097 m，则其截止频率fcut-off=3 763 Hz.计算模型可以采用1/4结构，采用四面体网格，网格数101 742(对进口管长度有要求)，SYSNOISE采用四面体网格，网格数为97 689.如图11比较FVM计算结果与SYSNOISE(FEM)计算结果，两者吻合良好.

图12将该消声器的出口管改到侧面，则该消声器结构的进出口边界由轴对称变为只有一个对称面，这就导致计算网格会成倍的增加，同时为测试FVM方法在大型计算中的能力，划分网格数为560 567，对于在本节提到的计算机上，FEM是难以完成的，而FVM方法由于采用显格式推进，并且不需要计算大型的刚度矩阵，所以可以比较快捷的处理较大型计算问题.如图13对比出口管位置不同时消声器的传递损失.可以看出，当出口改为侧面后在共振频率处出现多个峰值:第一个共振峰是由于侧面出口管与消声器端部之间形成一个声学共振腔引起的;后面的几个共振峰是由于侧面出口使消声器变为非对称结构而激发周向模态，消声器的高阶模态提前激起而形成的，出口改到侧面后消声器的整体声学性能得到改善，特别是1 500 Hz以下的声学性能得到明显提高.

图10 轴向出口蒸汽消声器结构示意Fig.10 Sketch map of steam muffler with end-outlet

图11 端部出口蒸汽消声器传递损失比较Fig.11 Transmission loss of steam mufler

图12 侧面出口消声器结构示意Fig.12 Sketch map of steam muffler with side-outlet

图13 端部出口与侧面出口蒸汽消声器传递损失比较Fig.13 Comparison of transmission loss of end-outlet and side-outlet

4 结论

基于三维交错的非结构FVM求解线性声波动方程，综合运用全反射与全透射边界条件，采用时域脉冲法预测消声器的传递损失.在消声器的进口给定脉冲后，再设为全透射边界条件，可以有效地消除进口处的反射波，使数据处理更加简便，并且不需要限制出口管长度.详细推导出三维四面体网格的稳定性判据及色散关系.

利用三维FVM方法计算外插管消声器的声学性能，计算结果与商业软件SYSNOISE计算结果、文献中实验结果吻合良好，表明该方法能够比较准确地预测无流情况下抗性消声器的声学性能.通过分析复杂结构消声器的声学性能，展示FVM时域方法处理大型模型的能力，同时，结果表明侧面出口管相对端部出口管可以改善蒸汽消声器的声学性能.

［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］刘恒，廖振鹏.波动数值模拟的一种显式方法——二维波动［J］.力学学报，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］王福军.计算流体动力学分析——CFD软件原理与应用［M］.北京:清华大学出版社，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］杜功焕，朱哲民，龚秀芬.声学基础［M］.南京:南京大学出版社，2001:176-177.

［18］李太宝.声场的方程和计算方法［M］.北京:科学出版社，2005:59-61.

［19］居鸿宾，沈孟育.气动声学数值模拟中无反射边界处理及声源模型的建立［J］.上海交通大学学报，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.