Hermite径向基函数点插值配点法求解消声器横向模态

方春慧，刘凤景，周 予，方 智

（1.烟台汽车工程职业学院，山东 烟台 265500；2.盛宝油压机械（烟台）有限公司，山东 烟台 265500；3.华中科技大学 船舶与海洋工程学院，武汉 430074）

求解消声器横向模态常用的方法有解析方法和数值方法。解析方法[1–2]求解简单、快速、精确，但是仅局限于有解析方程的横截面，比如圆形截面，矩形截面以及椭圆形截面。对于任意形状的横截面横向模态的求解，只能使用数值方法[3–4]，现在较为成熟的常用的数值方法为二维有限元方法（FEM）和边界元方法(BEM)。该类方法均需要依赖于网格进行求解，是基于网格的求解方法。无网格方法(MFM)[5]是一种相对较新的数值方法，求解时节点是相对于网格自由存在的，甚至不需要网格，本文所用到的径向基函数点插值配点方法是一种真正的无网格方法，场节点自由任意分布在问题域和边界上。无网格方法最先在力学领域应用较多，近十几年开始使用在声学领域。1998年，Bouillard应用基于移动最小二乘的无网格伽辽金法计算了二维内部声学问题，数值算例比较结果表明该方法比有限元法的计算精度高，计算速度快[6]。2002年，Chen等应用边界配点的无网格方法计算了二维声腔的声学特征值问题，并指出该方法较有限元法简单，容易实施[7]。2010年，大连理工大学的张宏伟等应用基于点插值的配点型无网格法求解了Helmholtz方程，通过数值算例比较证明了该方法具有较高的精度和良好的收敛性[8]。2011年，湖南大学的姚凌云等将分区光滑径向点插值无网格法应用于二维声学分析中，对管道和二维轿车声学问题的数值分析结果表明，该方法比有限元法具有更高的精度[9]。2012年，海军工程大学的胡海等将边界无网格法应用到结构声辐射计算中，通过与解析结果对比表明该方法具有更高的插值和计算精度[10]。华中科技大学的黄其柏课题组将径向配点无网格法应用到气动声学和腔体结构声学的研究中，取得很好的成果[11]。Christina研究了使用径向点插值方法（RPIM）求解二维Helmholtz方程时的色散效应[12]。但是当导数边界条件存在时，使用RPIM求解会出现不稳定。鉴于此，一种新型的处理方法：Hermite RPIM方法被提出处理导数边界条件。

本文使用多项式Hermite径向基函数点插值法求解问题域内计算点的形函数，使用配点方法离散二维横向模态方程，进而求解消声器横向本征波数以及模态形状，研究形状参数对计算精度和速度的影响。

1 Hermite型RPIM求解形函数

对于膨胀腔消声器，二维声压控制方程可以表示为

其中pxy为横向声压，为二维笛卡尔坐标系的拉普拉斯算子，kxy为横向方向的本征波数。

对于刚性壁，横向边界条件可以表示为

使用Hermite径向基函数点插值法构造形函数。场节点随机分布在问题域Ω和边界Γ上。对于任意一个计算点，可以形成一个影响域。影响域可以是任意形状，本文选取圆形。问题域和影响域的示意图如图1所示。假设在影响域内有n个场节点（包括导数边界上的节点），导数边界上的节点数目为nDB。计算点的声压ph可以写成如下形式

图1 问题域和影响域

其中Ri(x)是径向基函数，dc为位于计算点附近的节点间距。αc和q是无量纲的形状参数。xi是影响域内每个节点的坐标，ai是相应的系数。是位于导数边界上的点的坐标，bj是相应的系数。qk是多项式，m为多项式的项数，ck是多项式的系数。n为导数边界上第j个计算点的单位外法线向量，lxj和lyj为方向余弦。

将方程(3)写成如下矩阵形式

方程式(3)中的系数可采用通过影响域中的所有n个节点的函数值的插值以及等于导数边界节点的函数导数值而确定。

对局部支持域中的所有n个节点（包括导数边界上的节点）可得到

式中l=1,2…,n。

对所有导数边界上的节点可得到

式中l=1,2…,nDB。

为获得唯一解，可以施加以下约束条件

联立方程式(10)至式(12)可得到如下矩阵方程

其中，广义力矩矩阵G为对称矩阵，其中每个矩阵可以表示为

求解方程式(13)可以得到声压表达式为

其中，形函数向量φ可以写成

进而，可以得到声压近似函数及其导数表达式

2 配点法求解横向模态

一般问题域内场节点是任意分布的，为了方便起见，配点取与场节点相同的点。假设问题域内部节点数为Nd，边界节点数为Nb，总的节点数为N=Nb+Nd。若所研究的影响域内包含导数边界上的节点，对于影响域内部节点xI，将式(20)代入到式(1)中可以得到以下矩阵方程式

KI和MI分别是影响域内部节点xI的声刚度和声质量矩阵。值得注意的是，在所有内部节点都需建立方程式(25)。

导数边界上的节点需要满足两个方程，一个为方程式(25)，另一个为如下形式

值得注意的是，在所有导数边界节点上均需要建立方程式(28)。

若所研究的影响域内不包含边界节点，则方程式(25)和式(28)均可简化。

组装所有节点声刚度和声质量矩阵，可得到如下系统本征方程。

与传统有限元方法不同的是，在无网格配点方法中，矩阵组装的具体方法是将其节点矩阵按行摆放在一起组成总体矩阵。

求解方程式(31)可以得到相应的本征波数和本征向量矩阵。

3 结果与讨论

由理论分析可知，本征波数计算误差主要来自三部分：径向基函数形参的选择，导数边界条件的处理以及配点布置。鉴于圆形截面的本征波数有解析解，本文采用Hermite型RPIM配点法求解图2所示的圆形截面的本征波数，并进行其误差分析。圆形横截面的半径取r=0.024 3 m。本征波数数值结果与解析结果的相对误差定义为

如图2所示，圆形截面内的节点分布是任意的，由MATLAB编写简单程序得到。

图2 圆形横截面示意图

3.1 场节点数目的影响

场节点数目对前6阶横向本征波数的计算误差的影响如图3和图4所示。

图3 问题域内部节点数目的影响

图4 问题域边界节点数目的影响

从图中可以看出，本征波数的计算精度并没有随着节点数目的增加而增大。对于本文研究的圆形横截面，当内部节点数目超过100，边界节点数目超过60，本征波数的计算误差反而会增大，这是因为随着节点增多，广义力矩阵G的条件数也会随之增大，会造成计算结果的误差增大。所以在计横截面本征波数时，需要选择合适的节点数目。对于图2所示的圆形横截面，内部节点和边界节点分别取100和40最为合适。

3.2 影响域尺寸的影响

图5给出了影响域内包含不同个数节点时的本征波数的计算误差值。与问题域内场节点数目的影响趋势大致相同，本征波数的计算精度并不是随着影响域增大而增大。

图5 影响域尺寸的影响

由图5可以看出，当影响域内节点数目超过35，计算误差反而增大，这也是由于随着影响域尺寸的增大，广义力矩阵条件数增大的原因。所以在选择计算点的影响域尺寸时，需要平衡计算精度和计算速度的关系。

3.3 径向基函数形参的影响

图6和图7给出了径向基函数形参对本征波数计算误差的影响。形参的选择对计算误差的影响比较大，本文通过多次多参数尝试，选择了几个常用的有相互关联的形状参数q和ac，并且研究了它们的影响。

由图6可以看出，当q=1时，相对误差基本上为100%，数值计算结果不正确。当q=±0.5时，复合二次径向基函数为标准径向基函数，但是数值计算结果并没达到最大的精度，反而当q=0.98和q=1.03时本征函数的计算误差最小。

由图7可以看出，计算误差随着ac的增大而减小，但是存在一个计算误差最小的形参值ac=3.5。可以得出结论，对于本文研究的二维声学问题，最优的径向基函数形参可以取为q=0.8或者q=1.03和ac=0.35。但是需要注意的是，对于不同的问题，甚至是相同的问题不同的横截面形状和尺寸，最优的形状参数的取值会发生变化。

图6 形参q的影响

图7 形参ac的影响

3.4 计算速度的比较

为了验证本文方法和编写的计算程序的正确性和适用性，本文分别使用二维有限元方法和本文的无网格方法计算了如图8所示的不规则结构跑道圆横截面的本征频率和模态形状。跑道圆的尺寸为长短轴分别为a=0.25和b=0.15 m。横截面的本征波数和模态形状分别如图9和图10所示。

图8 跑道圆横截面示意图

由图可以看出，有限元方法和无网格方法计算的本征频率的相对误差在5%以内，再次验证了本文方法的正确性。

图9 跑道圆截面本征频率

图10 跑道圆截面模态形状

为了比较本文方法的计算速度，表1给出了分别使用有限元方法和无网格方法计算图2和图8所示横截面的时间比较。为了得到相同的精度，有限元方法和无网格方法计算本征模态时选取相同的节点数，本文中圆形截面节点数设置为140，跑道圆截面节点数为1 140。表1中所示的计算时间为求解本征波数的时间，不包括前处理时间，比如有限元方法中生成网格的时间，无网格方法中前处理时间即节点生成的时间可以忽略不计，因为节点可以通过简单的MATLAB语句生成。由表1可以知道，无网格方法比有限元方法节省时间，尤其是计算较大复杂的横截面时无网格方法的快速性更为突出。

但是，由前面的分析可知，影响无网格方法计算精度的因素很多，尤其是径向基函数形状参数的选择，不同的计算结构最优的形参取值不同，域内节点数目以及每个计算点影响域的尺寸的选取也将不同，这也是与有限元方法相比无网格方法的不足之处。

表1 本征波数求解速度比较/min

3.5 跑道圆形双出口膨胀腔消声器横向模态计算

为了体现本文方法的有效性，使用该方法计算了如图11所示的跑道圆形双出口插管膨胀腔消声器每个特征横截面的本征模态。膨胀腔跑道圆横截面的长短轴长度分别为a=0.218 m和b=0.109 m，进出口管直径为d1=d2=d3=0.048 6 m，两出口偏移量均为δ1=δ2=0.054 5 m 。该消声器共有SA、SB、SC、SD、SE、SF、SG7个特征横截面，由于对称性以及进出口管尺寸相同，所以特征横截面SA、SE、SG具有相同的本征模态。表2列出了本文无网格方法计算的本征频率和模态形状。

图11 跑道圆形双出口插管膨胀腔消声器

表2 消声器特征截面横向模态

4 结语

本文提出Hermite型径向基函数点插值配点方法求解消声器的横向模态。通过比较圆形和跑道圆形横向本征波数，验证了本文提出的方法的计算精度和速度。进而分析了问题域内场节点数目，影响域的尺寸以及径向基函数的形状参数对本征波数计算误差的影响。对于每一个影响因素的选取，都存在一个最优值。虽然相比有限元方法，本文提出的径向基函数插值配点法计算精度和计算速度均具有优势，但是对于每一个影响因素最优取值的选取，则需要重复性地尝试，这也是该方法不足之处。

[1]SELAMET A,XU M B,LEE I-J.Analytical approach for sound attenuation in perforated dissipative silencers with inlet/outlet extension[J].J.Acoust.Soc.Am.,2005,117(4):2078-2089.

[2]KIRBY R,DENIA F D.Analytic mode matching for a circular dissipative silencer containing mean flow and a perforated pipe[J].J.Acoust.Soc.Am.,2007,122(6),3471-3482.

[3]ZHAO S.On the spurious solutions in the high-order finite element methods for eigenvalue problems[J].Computer Methods in Applied Mechanics and Engineering,2007,196(49-52):5031-5036.

[4]PETYT M,KOOPMAN G H,PINNINGTON R J.The acoustic modes of a rectangular cavity with a rigid incomplete partition[J].J.Sound.Vib,1997,53(1):71-82.

[5]LIU G R,GU Y T.An introduction to meshfree methods and their programming[M].New York:Springer,2005.

[6]BOUILLARD P,SULEAUB S.Element-free galerkin solutions for helmholtz problems: formulation and numerical assessment of the pollution effect[J].Computer Methods in Applied Mechanics and Engineering,1998,162,317-335，

[7]CHEN J T,CHANG M H,CHEN K H,et al.Boundary collocation method for acoustic eigen analysis of threedimensionalcavities using radialbasis function[J].Computational Mechanics,2002,29(4):392-408.

[8]李美香，张宏伟，李卫国，基于点插值的配点型无网格方法解 Helmholtz问题[J].计算力学学报，2010，27（3）：533-536.

[9]姚凌云，于德介，臧献国.声学数值计算的分区光滑径向点插值无网格法[J].振动与冲击，2011，30(10)：188-192.

[10]胡海，郭文勇，马龙.结构声辐射计算的边界无网格法[J].噪声与振动控制，2012，32(5)：37-41.

[11]LI K,HUANG Q B,WANG J L,et al.An improved localized radial basis function meshless method for computational aeroacoustics[J].Engineering Analysis with Boundary Elements,2011,35(1):47-55.

[12]CHRISTINA W,OTTO V E.Dispersion analysis of the meshfree radialpointinterpolation method forthe Helmholtz equation[J].Int J Numer Meth Eng,2009,77(12):1670-89.