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
其中：pta为介质Ⅲ中的透射波声压幅值，pia为入射波声压幅值；
为媒质Ⅰ、Ⅱ的密度；c1、c2为媒质Ⅰ、Ⅱ中的声速；
，ω为角频率。
平面波激励矩形板的辐射声示意图如图 2所示，其中板厚h的值较小，属于薄板。xOz平面位于板的下表面。矩形板x方向长为l1，y方向长为l2。当声波的波向量在沿xOz平面内与Oz轴交角为θ方向入射，则在z＜0半空间内除了刚硬板面入射波和无限大刚板反射波之外，还应加上有限板在声波作用下进行振动时再辐射的声波。在z＜h半空间只存在再辐射声波，于是在板面上下介质中的声关系如下面第2节中的公式。
图2 平面入射波的传播示意图
Fig.2 Schematic diagram of plane incident wave propagation
入射波声压
，为简化推导，令声压幅值pa为1，经过一系列推导可以得到[1]：
其中：µm n 是对应于本征函数
的本征值，它们决定于板的边界条件；µ为泊松比；E为杨氏模量；D′为损耗因子；ρ为介质密度。令(x1, y1)为矩形板表面上的点坐标，(x1′, y1′)是面上积分的动点坐标，η为板表面位移幅值。Bmn为有阻尼情况下的振动位移幅值展开项的系数。
1.1 矩阵求解
式(5)中有两个求和符号，需用下方矩阵形式进行展开：
1.2 辐射声求解
由于整个过程是采用复数计算的，所以
为复数，取其幅值作为平板的辐射声，由于入射波的声压幅值设为1，此结果即为透声量。
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
当第Ⅴ部分边长L设为0.0001～0.01(与厚度同量级)，对于本文中计算对象为长1 m、宽1 m、厚0.001 m 的板，计算结果与采用 Mathematic的NIntegrate命令计算二维积分十分接近。在本文计算中，选择L≤0.000 5，与MatheMatic 8.0软件求解结精度有4位有效数字一致。其在L=0.000 5时，
(式中k=2π ×261.6/344)，可见这一部分是很小的。
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 计算结果讨论
计算阶数
中矩阵元素为 284(614 656)个，形成的
矩阵为784阶方阵。区间划分个数16，较窄边长为1个。其中D′= 0 .1(值的选择对结果影响不大)。在m=n=28时，这时对应于(28, 28)模态的固有频率为310.151 Hz，略大于计算频率。在距离矩形板 6 m处的平面上布置25×25的计算点，作为辐射声的计算接收点。程序采用 Fortran语言编写，采用双精度，在 Win7操作系统(CPU为 i7-7700)采用 Intel Visual Fortran2013软件计算，结果如图7所示。
图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 结 论
平面波垂直激励矩形板的辐射声，貌似十分简单，但若要精确解决却需要很繁难的推导与大计算量的数值计算。
本文细致地推导了平面波激励矩形板的辐射声公式，重点处理了矩阵求解问题、四重带奇点积分，并进行了程序编写，其中四重带奇点积分能够划分成五个区域分别求解来最终获得。针对某一具体算例，求得的结果与无限大板的透射声理论值十分接近。这说明采用模态正交分析是可行的。另外，采用这种办法会准确得到各个位置处的辐射声。这种办法再加以拓展还可以处理曲面辐射声的问题。
[1] 何祚镛. 结构振动与声辐射[M]. 哈尔滨: 哈尔滨工程大学出版社,2001.
[2] 何祚镛. 水声作用下矩形弹性-粘弹性复合板的振动和散射声近场(Ⅱ): 矩形弹性-粘弹性复合板散射声近场研究[J]. 声学学报, 1986,11(1): 1-19.HE Zuoyong. Vibration and near scattering field op immersed rectangular elastic-viscoelastic composite plate in an underwater sound field (Ⅱ)—scattering field near by composite plate[J]. Acta Acustica, 1986, 11(1): 1-19.
[3] LI W L, GIBELING H J. Determination of the mutual radiation resistances of a rectangular plate and their impact on the radiated sound power[J]. Journal of Sound and Vibration, 2000, 229(5):1213-1233.
[4] DAVIS P, RABINOWITZ P. Methods of numerical integration[M].New York: Academic Press, 1975.
[5] 沈杰罗夫. 水声学波动问题[M]. 何祚镛, 赵晋英, 译. 北京: 国防工业出版社, 1983.
[6] 戴遗山. 舰船在波浪中运动的频域与时域势流理论[M]. 北京: 国防工业出版社, 1998.
[7] LEE J H, KIM J. Analysis of sound transmission through periodically stiffened panels by space-harmonic expansion method[J]. Journal of Sound and Vibration, 2002, 251(2): 349-366.
[8] 孙曜, 汪鸿振. 双层板声透射性能计算分析[J]. 振动与冲击, 2003,22(4): 94-96.SUN Yao, WANG Hongzhen. Performance calculaton and analysis of sound transmission for double panel structure[J].Journal of Vibration and Shock, 2003, 22(4): 94-96.
[9] 辛锋先, 卢天健, 陈常青. 轻质金属三明治板的隔声性能研究[J].声学学报, 2008, 33(4): 340-347.XIN Fengxian, LU Tianjian, CHEN Changqing. Sound transmission through lightweight metallic sandwich panel with corrugated core[J]. Acta Acustica, 2008, 33(4): 340-347.
[10] 陈美霞, 罗琦, 魏建辉. 湍流激励下简支平板振动特性的解析法研究[J]. 船舶力学, 2016, 20(4): 487-496.CHEN Meixia, LUO Qi, WEI Jianhui. Analytical method for calculating the vibration characteristics of simple-supported plate excited by turbulent boundary layer[J]. Journal of Ship Mechanics,2016, 20(4): 487-496.
[11] 金叶青, 姚熊亮, 庞福振, 等. 均匀流中剪切变形加筋层合板声与振动特性研究[J]. 物理学报, 2013, 62(13): 335-347.JIN Yeqing, YAO Xiongliang, PANG Fuzhen, et al. Vibro-acoustic characteristics of shear deformable stiffened laminated panels in mean flow[J]. Acta Physica Sinica, 2013, 62(13): 335-347.
[12] ANDERSON J S, BRATOSANDERSON M. Radiation efficiency of rectangular orthotropic plates[J]. Acta Acoustica United with Acoustica, 2004, 91(1): 61-76.
[13] XIE G, THOMPSON D J, JONES C J C. The radiation efficiency of baffled plates and strips[J]. Journal of Sound and Vibration,2005, 280(1-2): 181-209.
[14] PUTRA A, THOMPSON D J. Sound radiation from rectangular baffled and unbaffled plates[J]. Applied Acoustics, 2010, 71(12):1113-1125.
[15] 杨德庆, 柳拥军, 金咸定. 薄板减振降噪的拓扑优化设计方法[J].船舶力学, 2003, 7(5): 91-96.YANG Deqing, LIU Yongjun, JIN Xianding. Structural topology optimal design to reduce vibration and noise of thin plate[J].Journal of Ship Mechanics, 2003, 7(5): 91-96.
[16] SUZUKI S, MARUYAMA S, IDO H. Boundary element analysis of cavity noise problems with complicated boundary conditions[J].Journal of Sound and Vibration, 1989, 130(1): 79-96.
[17] WU T W, DANDAPANI A. A boundary element solution for sound transmission through thin panels[J]. Journal of Sound and Vibration, 1994, 171(2): 145-157.
[18] COYETTE J P. The use of finite-element and boundary-element models for predicting the vibro-acoustic behaviour of layered structures[J]. Advances in Engineering Software, 1999, 30(2):133-139.
[19] DING W P, CHEN H L. A symmetrical finite element model for structure-acoustic coupling analysis of an elastic, thin-walled cavity[J]. Journal of Sound and Vibration, 2001, 243(3): 547-559.
[20] 杜功焕, 朱哲民, 龚秀芬. 声学基础[M]. 南京: 南京大学出版社,2001.
[21] 常庚哲, 史济怀. 数学分析教程 下册[M]. 3版. 合肥: 中国科学技术大学出版社, 2013.
[22] 李庆扬, 王能超, 易大义. 数值分析[M]. 武汉: 华中理工大学出版社, 1981.
[23] 徐士良. FORTRAN常用算法程序集[M]. 北京: 清华大学出版社,1995.
[24] JUNGER M C, FEIT D. Sound, structures, and their interaction[M].2nd Edition. Cambridge Massachusetts: The MIT Press, 1972.
[25] 何祚镛, 赵玉芳. 声学理论基础[M]. 北京: 国防工业出版, 1981.