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

在消声器的设计阶段多采用传递损失作为衡量消声器声学性能的指标，因此，能够准确的计算出消声器的传递损失对消声器优化设计有很大的实用价值.目前，用于预测消声器传递损失的方法包括频域法与时域法.频域法主要包括传递矩阵法［1］、解析法［2］、有限元法(FEM)［3］及边界元法(BEM)［4］，特别是一些基于FEM及BEM的商业软件已经在行业内广泛应用;时域法近年来得到快速的发展，可分为一维时域法与三维时域法，一维时域法主要采用有限差分法(FDM)［5-10］;三维时域法又称为消声器声学性能的计算流体动力学(CFD)方法，目前主要的研究方法是基于FLUENT等商业软件预测消声器的声学性能［11-13］.

时域法的优势在于易于执行，一次计算可获得整个频段内的信息，并且对计算机的内存依赖性不像频域法那样高，同时，频域法由于要计算复杂的矩阵，所以需要大量的内存，因此频域法难以计算大型的、复杂的模型，而时域法则有广阔的应用前景，完善这一数值模拟技术的主要问题是在保证结果可靠的基础上提高计算效率;就算法而言，提高计算效率的重要途径是发展显式方法［14］.

FVM是近年来发展非常迅速的一种离散化方法［15］，广泛应用于计算流体力学中，该方法能像FEM一样适用于不规则网格和复杂边界情况，且处理效率与FDM相似，远远高于FEM，所以在数值模拟中有着很大的发展潜力.因此，采用交错的非结构FVM、显格式推进求解线性声波动方程，计算消声器的声学性能，网格采用交错网格.这种交错网格即为文献［16］中提出的非结构化网格，它与FEM的网格离散方法相同，但又区别于FEM，因为它不需要计算单元刚度矩阵，节点上的变量只用来计算单元中心的空间导数;它与FDM交错网格的积分方法相同.

因此针对小振幅扰动波动方程采用三维交错非结构FVM解法，并研究这种数值方法的色散关系，进行稳定性分析，计算结果与商业软件SYSNOISE有限元计算结果、文献中的实验结果进行对比.另外，该方法已经在内部代码GTEA中实现.

1 声学数学模型

1.1 控制方程

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

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

1.2 数值方法

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

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

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

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

如图1所示为交错网格示意图，梯度在四面体单元1234内是均匀的，压力存在节点上，四面体单元1234对节点1的贡献可以表示为▽p·A i，其中A i为空间多边形E12 O4 E13 O3 E14 O2(由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］可表示为

式中:p in为消声器进口处的入射声压，p tr末端为无反射条件下消声器出口处的透射声压，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)的有限差分方程为

式中:为 A 点 t时刻的压力，h 为立方体的边长.

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

式中: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－6 6 s，总计算时间为 0.2 s，消声器进口给高斯压力脉冲(T=10×10－4).

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

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

表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，则其截止频率f cut-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时域方法处理大型模型的能力，同时，结果表明侧面出口管相对端部出口管可以改善蒸汽消声器的声学性能.