求解远场声辐射问题的一种快速边界元法

高晟耀，王德石，朱拥勇

（海军工程大学 兵器工程系，武汉 430033）

声纳是利用回波声场蕴含的目标特征信息，来实现水下目标的探测、识别和定位。为克服背景噪声对其性能的影响，通常设计具有大空间指向性的换能器基阵并对工作带宽加以限制，其形状多采用曲面状，如球冠形、抛物面形及椭球面形等。目前在结构辐射或散射声场特性研究中，一方面限于工程结构往往含有不规则的几何外形，结构边界具有复杂的声学参量分布，难以利用特殊函数、分离变量等方法获得解析解；另一方面，与有限单元法相比，边界积分方程法具有空间降维、自动满足Helmholtz方程以及远场辐射条件等优点，这就使得该数值计算方法在声学领域得到了普遍的重视。

而当边界元法求解边界积分单元上的未知量时，一旦辐射频率等于特征频率，就会面临处理不定型的极限0 0问题，造成结构表面出现奇异积分，以及特征频率处积分算法呈现振荡特性，所以如何克服奇异性的固有缺陷，并提高计算效率和精度成为边界元算法中的研究焦点。1963年，Chen和Schweikert[1]将边界积分方程方法最早应用到声学领域，在计算任意形状结构声辐射的过程中发现了一种与弹性波、电磁波研究中类似的现象，即对应Dirichlet内问题的特征频率处的解并不存在。Schenck[2]在1968年指出该解是非唯一性而非不存在问题，并提出了一种改进的算法CHIEF法，此方法选取内域点的积分方程作为Helmholtz边界积分方程的辅助约束条件。但在驻波节点分布密集的高频段，CHIEF点愈加难以选择。之后，Burton和Miller[3]提出了另一类具有代表性的改进算法，将奇异积分方程式和对其边界法向求导后的超奇异积分方程式的纯虚数倍进行线性组合，可以求得任意波数下的唯一解，但繁琐的积分处理方式制约了数值计算的速度和效率。

为寻求传统边界元有效替代方法，Koopmann等人[4]运用了Treffetz法在求解边值问题的思路，将描述声场问题解表示为一系列满足齐次偏微分控制方程的基函数的叠加。通过在振体结构的内部虚拟面上配置单偶极子源，来等效模拟结构外辐射声场。由于虚拟曲面与实际振动边界不重合，有效避免了外域数值积分中的奇异性问题，而且所建立的间接边界元模型仅要求声单元特征尺寸小于波长而非1 6波长，适用频率范围更广。为了进一步解决虚拟曲面上相应的内问题的特征频率处出现非惟一解问题，Jeans、Mathews和向宇、黄玉盈分别提出了复数形式Buton-Miller型组合层势法[5]，以及在虚拟源强系统中增加虚拟阻尼的复数矢径波叠加法[6]，并且讨论了各自算法的计算效率问题。目前，该间接边界元法已经被广泛应用于近场声全息、噪声源识别、机械故障诊断等声学逆问题之中，取得了大量机理性和相关实验成果[7―9]，但涉及到远声场辐射特性的研究工作较少。

本文基于虚拟源强技术建立了结构声辐射的间接边界元模型，通过采用在内域虚拟曲面上配置一系列耦合极子源模拟结构外辐射声场的方案，研究了一种球冠形辐射器的远场声辐射特性，并与传统边界元法的计算结果进行了比较分析，从而为分析设计大空间指向性的辐射器提供了一种新思路。

1 间接边界元模型

考虑结构在理想流体介质中微小扰动形成的谐声场，去除时间相关性e-iωt，外域问题由Helmholtz方程，Neumann边界条件和Sommerfeld辐射条件描述

式中p(r)为场点r处复声压，k=ω/c为声波数，ω为声波的角频率，vn为边界面上法向振度。该偏微分控制方程的边值问题难以给出解析解，常将其转化为积分方程的形式，得到理论上能够计算任意形状振体辐射声场的Kchhoff-Helmholtz积分方程

传统边界元方法通过离散边界面上的积分方程，得到未知量{ps}和已知量{vn}的关系式

并利用离散化后的Helmholtz边界积分方程，推导出外辐射声场r处复声压为

式中N为离散单元数，矩阵A、B、C和D的元素见文献[10]。

图1 虚拟源强系统Fig.1 Virtual source intensity system

虚拟源强技术并非直接求解积分方程(4)，而是利用内域近离散单元形心处系列单偶极子源所辐射声场的叠加，构造虚拟源强系统如图1所示，从而模拟结构外辐射声场

式中sm为虚拟源强，αm，βm为已知常数，取值1，0时即简单源，表示无限大障板面上单元；取值0，ik时即偶极子源，表示仅弯曲变形薄壳单元；取值1，ik时即耦合极子源，表示封闭辐射体面单元。上述两式表示成矩阵形式为

式中矩阵[H]和[T]为联系等效源点和场点的传递函数，一旦配置虚拟源强系统后，传递函数可通过格林函数的计算得到。根据边界面上单元体积速度

结合已知边界法向振速求得虚拟源强{}s，之后代入(7)式和(8)式，可以得到辐射声场，从而实现对结构声辐射特性的研究。

2 算例与仿真分析

考虑到评价波叠加方法的计算精度，选取了具有解析解的球冠形辐射器作为研究对象，见图2，其远场点(r,θ,ϕ)处辐射声压理论值为

式中为第二类m阶球汉克尔函数，Pm为第m阶勒让德函数，θ0为球冠活塞面仰角，a为刚性球直径，文中θ0=45°，a=1m。

图2 球冠形辐射器Fig.2 Spherical cap radiator

以半径r=18上的点为观测坐标，为反映辐射声场的指向性，以各频率下的最大声压对声压幅值进行归一化处理，得到极坐标下的理论声压幅值随频率变化关系如图3所示。由于对于所有的辐射振膜来说，波束带宽是非常重要的远场声辐射特性，往往采用半辐射声压即声压响应衰减-6d B时点间角度的间隔横定，对该球冠形辐射器的半声压带宽进行分析，理论结果如图4所示。

根据直接边界元法的面元划分要求，即单元特征尺寸不大于λ6，离散球面为4 000个四边形单元，网格划分如图5所示，若单元特征尺寸为λ3，四边形网格的数目为1 000。

为比较两种数值算法在计算远时半声压带宽的误差，定义相对误差公式

图3 k a=4,12,20时远场归一化理论声压幅值Fig.3 Normalized magnitude of the theoretical pressure in far field atk a=4,10,20

图4 r=18时理论半声压带宽Fig.4 Variation of theoretical semi-bandwidth of sound pressure with frequency atr=18

图5 球冠形辐射器表面网格划分Fig.5 Surface mesh of the spherical cap radiator

式中An为数值计算半声压带宽，Aa为理论半声压带宽，n为球面网格划分的单元数目。

在1≤k a≤20范围内，经两种数值算法求得了不同单元尺寸的波束带宽相对误差曲线，如图7和图9所示。由上述计算结果和表1可以得知：网格密度由λ6调整至λ3时，间接边界元法的单个频率计算时间将近是直接边界元法的1 6，而且各自的计算效率分别提高了20倍和10倍；从计算精度角度看，直接边界元法计算波束宽度的相对误差水平始终维持在1%以内水平，间接边界元法在k a=14处的最大相对误差则由6.77%跳变到10.13%，但在多数频率范围内该方法的计算误差并未超过4%，从图6、8可以看出，所选择的相对误差较为敏感，间接边界元法的数值计算结果能够满足工程应用所允许的误差要求。

图6 单元尺寸取λ 6时3种方法计算的半声压带宽Fig.6 Semi-bandwidth of sound pressure with mesh densityλ 6

图7 单元尺寸取λ 6时的带宽相对误差Fig.7 Relative error in the bandwidth with mesh densityλ 6

表1 k a=4时三种方法的半声压带宽及其计算时间比较Tab.1 Acomparison between three approaches for analysis of bandwidth and computation time

3 结语

针对传统边界元法奇异性的固有缺陷，提出一种求解结构辐射声场的间接边界元法，应用该方法计算了球冠形辐射器的远场声辐射特性，将其结果与解析法以及直接边界元法的结果进行比较。分析表明，该间接边界元方法可在实波数域内获取唯一解，有效克服直接边界元法的不足。在保证工程应用所允许的误差要求条件下，由于所需表面网格划分的单元数目较少，其计算效率是直接边界元法的3～6倍。同时，此方法也为分析设计具有大空间指向性的辐射器提供参考。

图8 单元尺寸取λ 3时3种方法计算的半声压带宽Fig.8 Semi-bandwidth of sound pressure with mesh densityλ 3

图9 单元尺寸取λ 3时的带宽相对误差Fig.6 Relative error in the bandwidth with mesh densityλ 3

[1]Chen L H,Schweikert D G.Sound radiation from an arbitrary body[J].J.Acoust.Soc.Am,1963，35:1626-1632.

[2]Schenck H A.Improved integral formulation for acoustic radiation problems[J].J.Acoust.Soc.Am，1968，44(1):41-58.

[3]Burton A J,Miller G F.The application of integral equation methods to the numerical solution of some exterior boundary value problems[C].Proceedings of royal society of London,1971,A323:201-210.

[4]Koopmann G H,Song L and Fahnline J B,A method for computing acoustic fields based on the principle of wave superposition[J].J.Acoust.Soc.Am，1989，86(6)：2433-2438.

[5]Jeans R A，Matthews I C.The wave superposition method as a robust technique for computing acoustic fields[J].Journal of Acoustical Society of America,1992,92(2):1156-1166.

[6]向宇，黄玉盈.基于复数矢径的波叠加法解声辐射问题[J].固体力学学报，2004，25(1)：35-40.

[7]于飞，陈心昭，李卫兵，等.空间声场全息重建的波叠加方法研究[J].物理学报，2004，53(8)：2607-2613.

[8]李加庆，陈进，杨超，等.波叠加声场重构精度的影响因素分析[J].物理学报，2008，57(7)：4258-4264.

[9]薛玮飞，郭金泉，陈进，等.波叠加法在机械噪声故障特征提取中的应用研究[J].机械强度，2007，29(6)：900-903.

[10]Wu T W.Boundary Element Acoustic:Fundamentals and Computer Codes[M].U K:WITPress,2000.