平面波激励矩形板辐射声计算中奇异积分处理方法研究

李 亚，张 楠

(1. 中国船舶科学研究中心船舶振动噪声重点实验室，江苏无锡214082；2. 中国船舶科学研究中心水动力学重点实验室，江苏无锡214082)

0 引 言

声波透过均匀介质的现象在生活中比较普遍，对于无限延伸的介质，经典声学教材往往直接假定介质中的波是平面波，在界面处速度与压力连续，从而可以得到透声量。实际情况下，介质往往是有限的，而且存在边界上的限制。对于平面波激励矩形板的辐射声计算，在文献[1]中已有较为详细的分析：平面波激励矩形板时由于边界阻抗突变，板内波在边界处反射，因而在板中形成驻波，驻波振幅具有一定形式的分布，利用模态的正交性，可以得到板在介质中的振动方程，然后使用瑞利公式进行辐射声计算。

在文献[1]中求矩阵元素时出现了带奇点的四重积分。对于这一问题，在文献[2]中采用了变量替换与积分区域变换的办法，将四重积分变为两重积分，变换过程较为复杂，未专门分析奇点积分的求解方法。同样在文献[3]中，也采用了变量替换与积分区域变换方法，亦没有专门分析奇点的处理办法。对于一重积分可以使用LOBATTO法则计算[4]，以避免被积函数的奇异性。对于二重带奇点积分一般是以奇点为圆心作一个小圆，来研究当圆半径逼近0时的结果[5-6]，这些多用于理论分析。本文采用理论推导，给出了求解奇点积分的具体办法，然后再用数值求解。对于矩阵求解，将在正文中看到，里面出现了两个相加符号，需要巧妙设计形成矩阵再进行求解。

文献[1]中还细致介绍了带周期栅的膜振动及其声场、周期性支撑固定的薄板振动和声散射、带周期性加强肋的薄板振动和声辐射，这些内容由简到繁地涵盖了工程中常见的典型情况。应该说在2001年左右，平面波激励矩形板的辐射声理论研究已经达到一个顶峰。因此近年来，对于平面波激励矩形板的辐射声理论分析的研究并不多见，主要集中在复杂板的隔声问题[7-8]、辐射声特性[10-11]、辐射效率[12-14]、板的优化设计方面[15]。

另外，采用声学有限元、声学边界元也可以研究矩形板的辐射声。Suzuki等[16]将有限元、边界元两种方法耦合，即将声腔用边界元离散，将结构用有限元离散，计算了声波经由弹性单层板向声腔透射的传声损失。Wu等[17]用边界元法(Boundary Element Method, BEM)对薄壳结构内外声场进行了分析，对结构用有限元法(Finite Element Method,FEM)分析，建立复杂结构的内、外声场和结构耦合振动方程，分析了声波通过薄壁壳体的声传输问题。Coyette[18]基于有限元和边界元耦合方法，计算弹性和多孔材料结构的表面阻抗和传声损失。Ding等[19]计算了声源和结构载荷共同激励下弹性薄壁腔体的声振耦合特性。对于复杂情况、单一工况问题有限元或边界元法能够很好地解决，但难以发现机理与影响规律，在设计中还要从模态分析入手，以达到工程要求。

平面波激励矩形板的辐射声是很基础的问题，但由于是一个声固耦合的问题，因此这个几何形式简单的问题并不是看起来那样简单，仍有很多机理有待深入研究与挖掘，在声学领域中具有独特的研究价值。这方面研究的进步，也会给其他研究带来启发，因此有很大的研究意义。

1 基本理论公式推导

矩形板由于入射波的作用，会引起板振动，于是向两边介质中辐射声波，同时板面也受到两面介质对板面的反作用力，所以这是一个声固耦合的复杂问题。

当板较大时，可以近似用声波穿过无限延伸媒质上的透声理论进行计算。本节首先回顾一下透声理论，然后重点推导平面波激励矩形板的辐射声理论公式。

声波通过中间层的示意图如图1所示，透射波(pt, vt)通过无限大中间层的声压之比[20]tp为

图1 平面波激励矩形板辐射声波的示意图Fig.1 Schematic diagram of soundwaves from a rectangular plate excited by plane waves

平面波激励矩形板的辐射声示意图如图 2所示，其中板厚h的值较小，属于薄板。xOz平面位于板的下表面。矩形板x方向长为l1，y方向长为l2。当声波的波向量在沿xOz平面内与Oz轴交角为θ方向入射，则在z＜0半空间内除了刚硬板面入射波和无限大刚板反射波之外，还应加上有限板在声波作用下进行振动时再辐射的声波。在z＜h半空间只存在再辐射声波，于是在板面上下介质中的声关系如下面第2节中的公式。

图2 平面入射波的传播示意图Fig.2 Schematic diagram of plane incident wave propagation

其中：µmn是对应于本征函数的本征值，它们决定于板的边界条件；µ为泊松比；E为杨氏模量；D′为损耗因子；ρ为介质密度。令(x1, y1)为矩形板表面上的点坐标，(x1′, y1′)是面上积分的动点坐标，η为板表面位移幅值。Bmn为有阻尼情况下的振动位移幅值展开项的系数。

1.1 矩阵求解

式(5)中有两个求和符号，需用下方矩阵形式进行展开：

1.2 辐射声求解

2 计算对象与计算方法

2.1 计算模型

计算对象为一薄钢板，长度为 4 m，宽度为4 m，厚度为0.001 m，放置方向如图2所示，其中四边简支，钢板四周为无限大障板。钢的弹性模量216 GPa，密度为7 800 kg·m-3，体纵波声速为 6 100 m·s-1[19]，泊松比为 0.28。空气的参数：密度为 1.21 kg·m-3，声速为 344 m·s-1。声波频率为 261.6 Hz(音名 c1)，如果频率太高，则计算阶数增加，太小则考虑到近场效应平板面积要增大，这样积分时区间划分个数又增加，因此定为这一频率。

2.2 平面波通过矩形板的理论解

将2.1节中参数代入式(1)中，可以得到透声值为0.0648，这个值为无限大板的透声结果。由于板较大，可以近似作为平面波通过矩形板的辐射声。

2.3 平面波激励矩形板辐射的数值解

2.3.1 奇点处理

图3 奇点计算示意图Fig.3 Singularity computation diagram

图4 不同ε时解的变化情况Fig.4 Solution variation with different ε

另外，在编程求积分时，采用矩形以便于计算，这样π变为4，积分最终表达为

其中：SL为边长为L的小正方形。

在具体进行积分时，采用图5的示意图。首先根据(x1, y1)的位置，确定好积分区域。按从上往下的顺序分别考虑且这三种情况，在每种情况下又细分为的情况。对于正方形位于矩形中间的情况，可分为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ五个区域，其中L为边长(见图6)。在计算时，分别计算各区域的积分值，然后全部加在一起，其中区域Ⅴ采用方括号中的方法计算。对于处于边线位置，仍然采用同样办法，只不过区域Ⅰ中积分为 0，区域Ⅴ中积分只采用有填充部分的面积。

图5 积分计算示意图Fig.5 Integral calculation diagram

图6 积分计算区域划分示意图Fig.6 Partition of integral calculation area

2.3.2 积分处理

在矩阵求解时多个元素均出现多重积分，分别采用计算重积分的累次积分法。对于积分函数有振荡的情况，一般可用复合辛普森法或复合 COTES方法[22]，由于这里要考虑计算精度与计算时间，经比较后每重积分均采用复化柯特斯公式[23]，公式为

由于一般数学软件自带公式无法计算四重积分，在分析三重积分时，对于长为1 m、宽为1 m、厚为 0.001 m的板，x1(下面公式中用 x表示)设为0.62，窄边划分区间数为1，其他边为16，进行三重积分计算。用Mathematic软件进行三维验证，命令为NIntegrate[Sin[1Pi*z]*Sin[1Pi*w]*E^(-I(2*Pi*261.6)/344*Sqrt[(z-x)^2+(w-y)^2])/Sqrt[(z-x)^2+(w-y)^2]Sin[3Pi*x]*Sin[3Pi*y],{y,0,1},{z,0,1},{w,0,1}]/x-＞0.62，结果为0.0397 207+0.035 285 5j。程序计算结果为0.039 610 0+0.035 285 5j，误差不超过1%。

2.3.3 求逆处理

由于各个元素为复数，所以矩阵为复矩阵，复矩阵求逆采用全选主元高斯-约当(Gauss-Jordan)法[23]。

3 计算结果讨论

图7 辐射声的值计算结果Fig.7 Numerical calculation results of radiation sound

由图7可以看出，在正方形中间部分的结果是十分接近理论结果的，其中心处的值为0.061 462，比理论结果0.064 8(见2.2节中)少5.2%；在边线处，透声值降低。这说明采用这种办法求奇点与求辐射声是可行的。另外，在试算时，在 15～27阶变化时，对计算结果基本影响不大。这是因为频率较高时，振动模态为角模态[24]，即大部分振动区域的振动辐射声会相互抵消，角上的振动对辐射声有较大影响。

值得指出的是，在采用瑞利公式计算活塞声源辐射声时，有远场弗朗霍夫尔(Fraunhofer)区和近场费涅尔(Fresnel)区[25]，已经比较复杂。这里有多个振动模态，情况更为复杂。在本文编程计算时发现，当声场位置与平板的距离增大时，中心处与边线处的辐射声结果比较接近，但由于其远场也类似于活塞声源，具有远场特点，其辐射声值有所降低。

本文基于振动模态的计算方法有以下近似，这也是计算的误差来源：薄板厚度忽略不计；理论模型只考虑弯曲波振动辐射声；奇点计算时，在计算公式将推导中的圆区域变为正方形区域，近似认为奇点附近相等；柯特斯法求积分时区间数不够大；计算阶次选取还可以更高些；阻尼的设置对结果有一些影响；计算采用双精度，虽已很精确，但仍存在误差。

4 结 论

平面波垂直激励矩形板的辐射声，貌似十分简单，但若要精确解决却需要很繁难的推导与大计算量的数值计算。

本文细致地推导了平面波激励矩形板的辐射声公式，重点处理了矩阵求解问题、四重带奇点积分，并进行了程序编写，其中四重带奇点积分能够划分成五个区域分别求解来最终获得。针对某一具体算例，求得的结果与无限大板的透射声理论值十分接近。这说明采用模态正交分析是可行的。另外，采用这种办法会准确得到各个位置处的辐射声。这种办法再加以拓展还可以处理曲面辐射声的问题。