三维全局弱式无网格方法计算膨胀腔消声器声学模态

刘凤景，方春慧

(烟台汽车工程职业学院，山东 烟台 265500）

消声器在汽车发动机进排气系统中应用较多。消声器的声学模态对消声器的消声特性有着较大影响，可以根据声学模态特性改进消声器结构，进而优化消声性能。声学模态的计算通常有解析方法和数值方法。文献[1] 采用解析方法计算了圆形截面简单膨胀腔消声器的二维横向模态，分析了消声器进出口位置对消声性能的影响。解析方法具有简单快速的优点，但是只能应用于规则形状截面的模态计算。文献[2－3] 使用有限元方法计算了膨胀腔消声器和穿孔管消声器的二维横向模态，进而分析了模态对消声器消声性能的影响。有限元方法可以计算任意形状截面的横向模态，应用范围较广，只是计算精度与截面单元的数目和质量相关。基于此，无网格方法被提出用来计算消声器的横向模态[4－5]。无网格方法(MFM)不需要划分网格，利用一组散布在问题域中及边界上的节点表示该问题。无网格方法根据形函数的求解和系统方程的离散可以组合出很多方法。Chen等[6]应用基于径向基函数的边界配点无网格方法计算了二维声腔的声学本征值，文章验证了无网格方法与有限元方法相比，具有不需要构造复杂的背景网格，容易实施的特点。张宏伟等[7]基于点插值的思想，用三角函数作为基函数构造局部近似函数，用配点法离散微分方程，提出了点插值配点型无网格法，研究了Helmholtz方程的波传播问题和边界层问题，通过算例比较证明了随着节点的增加，该方法具有较高的精度和良好的收敛性。姚凌云等[8]推导了分区光滑径向点插值无网格法（CSRPIM），将其应用到管道和二维轿车的声学问题分析中，数值计算结果表明该方法比标准有限元法具有更高的精度和准确度。黄其柏等[9]使用改进的局部径向配点无网格法研究了气动声学和腔体结构声学特性，数值计算结果与实验结果对比表明，径向基函数无网格方法是一种真正的无网格方法，摆脱了对网格的依赖，具有精度高、编程简单等优点。2016年，曾向阳等[10]提出一种基于声粒子分布积分的无网格方法，计算矩形封闭空间中低频声场，计算结果与采用商用计算软件以及经典无网格法的结果的对比，证明了该方法在中低频段与传统无网格方法相比具有更高的精度。2017年，方智等[11]应用强式无网格配点方法求解任意横截面膨胀腔消声器的横向模态，使用模态匹配技术求解传递损失，提出一种基于无网格的数值模态匹配方法，计算结果与实验测量值吻合较好，且与三维有限元方法相比，计算速度具有一定的优势。但当消声器壁面导数边界条件较多时，使用强式无网格配点方法稳定性会变差。于是，基于伽辽金加权残数法的全局弱式无网格方法[5]被提出用于计算消声器的横向模态，与理论解吻合较好。

以上所有对消声器模态的计算均是二维的，只是计算了消声器的横向模态，并未分析三维模态效应。本文应用三维全局弱式无网格方法计算简单膨胀腔消声器的三维声学模态，使用径向基函数点插值法求解三维无网格形函数，应用伽辽金加权残数法对系统方程进行离散，采用三维全局弱式无网格方法用来求解消声器的三维声学模态频率和模态振型，进而分析轴向模态对消声器声学性能的影响。

1 声腔模态

图1为简单膨胀腔消声器示意图。膨胀腔内的三维声压控制方程为

图1 简单膨胀腔消声器

膨胀腔的边界条件为刚性壁边界条件

在声场无网格计算中，选择形函数来表示声压，如方程式(3)所示

应用格林公式，可得

将边界条件表达式(2)代入上式中，方程式(5)简化为以下线性方程

其中K、M分别为腔内的广义刚度矩阵和质量矩阵。为了求解这些系统矩阵，需要求解域内和边界上的积分。因此在弱式无网格方法中需要构造背景网格。背景网格与场节点可以是相互独立的，并且对网格质量的要求不太高[5]。

2 基于全局弱式无网格方法求解三维声学模态

2.1 基于无网格径向基点插值法求解形函数

文中采用径向基点插值法（RPIM）求解形函数。添加多项式的RPIM插值可以表示为

其中：Ri(x)为径向基函数，n为径向基函数的个数，qk(x)为空间坐标x中的单项式，m为多项式基函数的个数。方程式(7)中每一项的表达式如下所示

构造计算点x的支持域来确定系数矩阵中的ai和ck，其中支持域中的场节点数为n。支持域的大小由ds=αsdc确定，αs为支持域的无量纲尺寸，在本文中取αs=2.5，dc为计算点x支持域内的平均节点间距。在径向基函数Ri(rk)中，αc为形状参数，在文中取αc=2.5。第i个场节点和第k个场节点之间的距离rk为

使方程式(7)满足计算点x周围的每一个节点，将产生n个线性方程。求解以下方程即可得到系数矩阵a0。

其中，各个矩阵表达式如下

集果装置由伞形集果筒、搭扣、合页等组成，伞形集果筒共有4个落果端口，位于集果筒底部，通过对红枣树枝几何拓扑结构的统计分析[12]，确定集果筒高度为850 mm，顶面直径1 200 mm，底面直径240 mm。设计伞形集果筒如图5所示。

由于矩阵G是对称矩阵，因此可由方程式(12)求得系数矩阵

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

对应节点声压向量的径向点插值形函数Φ表达式为

最终，声压p(x)可以表示为形函数与节点声压向量的乘积形式

2.2 基于Galerkin弱式离散系统方程

文中背景网格采用四面体单元，采用Hammer积分求解系统矩阵。每个单元中取5个Hammer积分点，对每一个Hammer积分点构造支持域，求解每一个Hammer积分点的形函数Φ，进而计算背景网格单元的刚度矩阵和质量矩阵。系统矩阵的求解方程如下。

其中：nQ为单元e内的Hammer积分点个数，wQ为相应的加权系数，Je为单元的雅克比矩阵。

2.3 求解三维声学模态

假设腔内的节点数为N，通过式可以求得膨胀腔内波数向量（即模态波数向量）k和特征向量矩阵（即模态振型矩阵）Xxyz，其中k的维数为N，Xxyz的阶数为N×N。

3 计算精度验证与结果分析

为了验证本文方法的计算精度，文中计算了简单膨胀腔消声器的前23阶三维声学模态频率和模态形状。计算结果与有限元结果进行对比。图1为简单膨胀腔结构图，具体尺寸为：膨胀腔直径D=0.1532m，膨胀腔长度l=0.2823m，进出口管直径d=0.0486m。由于与膨胀腔体相比，进出管的直径比较小，所以在较高的频率范围（0～4150 Hz）内均为平面波传播，在所考察的频率0～3200 Hz范围内，消声器的高阶模态仅在膨胀腔内被激发。鉴于此只计算了膨胀腔的三维声学模态。图2为膨胀腔背景单元网格图，图3为膨胀腔场节点分布图。

图2 膨胀腔背景单元网格图

图3 膨胀腔场节点分布图

场节点的分布是通过简单的MATLAB程序构造的，背景网格根据商业软件ANSYS得到。由图3可以看出，场节点是任意分布的，与背景网格是相互独立的。

表1给出了分别使用全局弱式无网格方法和有限元方法计算的前23阶模态频率。图4列出了前23阶高阶模态中所有的模态振型。在模态振型中，表征声压为零的线称为节线。三维声学模态可以用周向节线数m、径向节线数n和轴向节线数z共同表达，即（m，n，z）。

由表1可以看出与有限元结果对比，前23阶高阶模态频率误差均在1%以内，说明采用全局弱式无网格方法计算膨胀腔三维声学模态的适用性以及代码的正确性。

为了比较所用的无网格程序与三维有限元方法的计算速度，分别记录了计算图1结构消声器前23阶模态所需的时间。使用无网格方法时节点数取为1775，背景网格数为2505，单元类型为四面体单元。计算时间为126 s。三维有限元方法中使用ANSYS划分三维网格，使用Virtual.lab计算声学模态，节点数取1775，网格数为8654，单元类型为1阶四面体单元，计算时间为36 s。由此可以看出与三维有限元商业软件相比，本文编写的无网格计算程序在计算速度上并没有优势，但是若计算的模型较为复杂，使用三维有限元方法时需要耗费精力划分精度较高的网格，而无网格方法的计算节点由MATLAB生成，对网格的要求也不高，其计算优势则会显示出来。此外，使用无网格计算程序计算模态时，计算的模态阶数与计算节点的数目相关，在此算例中，计算的模态阶数为1775阶，而使用Virtual.lab时计算的阶数为23阶，所以用无网格方法计算的时间较长也是有一定原因的。由此也说明本文使用的无网格计算程序有很大优化空间。

对照表1的模态频率和图4列出的18个模态振型图可以发现，由于膨胀腔具有轴对称特性，所有周向模态即m不等于0的模态均为重根模态，即模态频率相等，模态振型关于中心轴对称。对于重根模态在图3中仅列出了一个模态振型。

声学模态是消声器固有特性，可以根据模态特点分析消声器的消声特性。图5给出了应用声学有限元软件计算的膨胀腔消声器的传递损失曲线。为了便于分析，在图5中也给出了应用一维平面波理论计算的消声器的传递损失，在该计算中考虑了进出口管的端部修正[3]。

由传递损失曲线可以看出，两条曲线在（0，1，0）阶模态频率处出现明显偏离，三维有限元计算结果曲线不再具有拱形特性，说明膨胀腔内出现了横向高阶模态，由于一维平面波理论没有考虑横向高阶模态效应，所以计算结果不再精确。而在低于（0，0，1）阶模态频率范围内，2条曲线基本吻合，这是因为在平面波理论中考虑了轴向的高维模态效应，由于同时考虑了进出口的端部修正，所以在一维平面波理论中将消声器进出口截面处由于截面突变引起的可耗散的三维波效应也考虑了，所以两条曲线基本吻合。在拱形处出现的偏差可能是由于端部修正并没有将可耗散的三维波效应完全考虑所致。

图4 膨胀腔三维声学模态振型图

图5 同轴膨胀腔消声器传递损失曲线

由于周向模态关于中心轴对称，消声器进出口管又位于中心轴上，所以周向模态均没有被激发，而最先被激发的模态为第1阶径向模态，可以设置出口管位于该径向模态节线上，拓宽消声器的消声范围，这已在有关文献中分析过[1,3]。通过传递损失曲线还可以发现消声器的通过频率对应着膨胀腔的轴向高阶模态频率，所以为了消除通过频率，除了采用插管以外[3]，还可以将出口管布置在轴向模态节线处，采用侧支出口形式。根据图3中给出的模态振型图可以看出，（0，0，1）、（0，0，2）阶模态的节线分别位于膨胀腔长度的1/2、1/4和3/4处，可以设想将出口管分别设置在这些位置，可以将相应的轴向模态效应消除，即可消除第1个和第2个通过频率。

文中分别将出口管布置在膨胀腔长度的1/2和3/4处，计算传递损失。侧支出口消声器示意图如图6所示。

图6 侧支出口消声器示意图

分 别 取δ=L∕2=0.14115m 和0.211725m，其他尺寸与图1相同。应用三维有限元方法计算得到的传递损失曲线如图7所示。

图7 侧支出口膨胀腔消声器传递损失曲线

从图7可以看出，由于消声器进口管位于中心轴线上，所以周向模态无法被激发，所以在（1，0，0）、（2，0，0）这2个周向模态频率处，两种出口侧支结构消声器的传递损失曲线均出现了峰值。侧支出口管距进口管偏移量为L/2的消声器结构，出口管位于（0，0，1）模态节线上，所以在传递损失曲线中在该模态频率处出现了峰值，与进出口管同轴结构消声器相比，消除了第1个通过频率。观察曲线还可发现，该结构消声器同时消除了第3阶轴向模态（0，0，3）模态，这是由于（0，0，3）阶模态中一个模态节线也位于L/2处。不仅如此，该结构消声器还在（1，0，2）和（2，0，2）阶模态频率处展现了较好的消声性能，在该频率处传递损失出现峰值。同理，侧支出口管距进口管偏移量为3L/4的消声器结构，出口管位于（0，0，2）模态节线上，消除了（0，0，2）阶模态的影响，分别在（0，0，2）、（1，0，1）和（2，0，1）模态频率处展现较好的消声性能。与进出口同轴消声器相比，侧支出口的消声器结构，若把出口管布置在某阶轴向模态的节线处，可以消除该阶模态的影响，大大优化消声性能。

4 结语

采用三维全局弱式无网格方法计算消声器的三维声学模态频率和模态振型，使用径向基函数点插值方法求解三维无网格形函数，使用伽辽金加权残数法离散系统方程，最终求得消声器的三维声学模态。计算了一简单膨胀腔的前23阶三维声学模态频率，与有限元计算结果对比，相对误差控制在1%以内，验证了采用三维无网格方法计算声学模态的正确性。进而根据声学模态的特性，分析了轴向模态对消声器消声性能的影响。将出口管侧支布置在轴向模态的节线处，可以消除相对应的通过频率，进而优化消声性能。