声学边界元拟奇异积分计算的自适应方法

缪宇跃, 李天匀,3, 朱 翔, 张冠军, 郭文杰

(1.华中科技大学 船舶与海洋工程学院, 武汉 430074; 2.船舶与海洋水动力湖北省重点实验室, 武汉 430074;3.船舶振动噪声重点实验室, 武汉 430064)

声学边界元拟奇异积分计算的自适应方法

缪宇跃1,2, 李天匀1,2,3, 朱 翔1,2, 张冠军1,2, 郭文杰1,2

(1.华中科技大学 船舶与海洋工程学院, 武汉 430074; 2.船舶与海洋水动力湖北省重点实验室, 武汉 430074;3.船舶振动噪声重点实验室, 武汉 430064)

提出一种自适应方法计算声学边界元中的拟奇异积分，通过单元分级细分将总积分转移到子单元上以消除拟奇异性。在此方法基础上深入研究拟奇异性，进一步提出接近度的概念，其中临界接近度可作为拟奇异积分计算的理论依据，并可用于预估拟奇异性是否存在。此方法的积分精度可调控，且不受场点位置限制，相比于已有方法更加灵活高效。数值分析表明拟奇异性强弱由场点与单元的相对位置决定,单元上远离场点的区域拟奇异性很弱，无需处理。研究结果为处理边界元法中的拟奇异性问题提供了新的选择和参考。

边界元；拟奇异性；自适应；接近度

声学边界元法经过几十年的发展, 已经成为求解结构振动声辐射问题的重要方法，得到了比有限元法更广泛的应用，尤其与快速多极算法结合后[1]，体现出较强的实用性。奇异积分问题[2,3]在边界元法中由来已久，已有多种处理方法，随后拟奇异积分问题也成为研究热点。拟奇异性，又称为近奇异性[4]或几乎奇异性[5]，是由于场点十分靠近结构表面，导致场点与源点间的距离很微小。虽然其物理机理不存在奇异性，但在数值计算时，这种微距使常规的高斯积分法失效，产生类似奇异性的收敛性差问题。

为处理边界元法的拟奇异性，国内外研究者提出了多种方法。TANAKA等[6]提供了一个较详实的关于拟奇异积分正则化方法的综述文献。ZHANG等[7]提出了一种半解析法，采用曲面单元，有效地处理了薄壁结构的拟奇异性。针对二维弹性问题，NIU等[8]提出一种解析的拟奇异积分算法，随后这种方法被应用到二维的正交各向异性势问题[9]和各向异性势问题[10]中。对于三维势问题，NIU等[11]将奇异性降阶，并结合正则化方法能够更容易地计算拟奇异积分。

自适应方法源自自动控制优化领域，并在有限元法中得到发展，随后这种自适应的思想被引入边界元法的研究中，但是主要被用来进行整体或局部的单元细分以提高计算精度[12]，KITA等[13]对自适应边界元法的发展应用做了较详实的概括。根据自适应单元细分的特性，研究者开始将这种方法应用于拟奇异积分的计算中。WEI等[14]基于CFD理论，结合有限元和边界元法分析了水下航行体噪声，其中边界元的奇异性问题采用自适应单元平均细分法处理。CROAKER等[15]对气流中结构的诱导声散射进行预报研究，为了得到近场压力和压力梯度，对拟奇异性体积分和面积分采用了单元细分法。DONG等[16]提出非稳态导热问题的边界元算法，为避免小时间步长的奇异性，采用坐标变换结合体单元细分法来提高域积分精度。DAVIES等[17]提出三维弹塑性边界元法研究复杂的地质力学问题，将域积分变换为边界积分消除强奇异性，并运用自适应细分法优化积分过程，其中单元细分过程依据经验公式进行。GAO等[18]进一步将自适应单元细分法和解析法相结合计算多种类型的二维奇异积分并取得良好效果。LI等[19-22]同样依据这种经验公式进行单元细分来处理三维问题的拟奇异性，其中子单元尺寸以及子单元与源点距离需要通过反复迭代计算和误差分析确定，而这种经验公式未充分考虑场点与单元相对位置的不同所存在的差异。ZHANG等[23-24]利用边界点法和边界面法研究三维势问题，拟奇异积分通过自适应单元细分实现，其细分准则是：比较每个子单元的对角线长度和近场点到该子单元中心的距离，若前者小于后者则认为细分达到积分要求，否则继续进行细分以达到积分要求。依据这种“对角线准则”的自适应细分法能够准确计算拟奇异积分，但若运用在声学分析中则过于苛刻，增加不必要的计算量。

本文提出分级细分自适应方法，制定合理易行的细分方案，将多种情况下的拟奇异积分计算统一起来。在深入研究场点与单元不同的相对位置所存在差异的基础上提出接近度的概念，为拟奇异积分计算提供理论基础，并预估拟奇异性程度以确定是否需要处理拟奇异性。数值分析证明了这种自适应方法能够很好的计算拟奇异积分，具有较好的灵活性和高效性，研究结果可为自适应边界元法的发展提供思路和参考。

1 声学边界元法基本理论

对于自由空间的声辐射问题，其Helmholtz边界积分方程为：

C(z)p(z)=

(1)

将基本解中的指数项按Taylor级数展开：

(2)

并设:

(3)

则有:

(4)

(5)

(6)

可见基本解G的弱奇异性是由Gr引起的，Gr是与频率无关的，而Gk是无奇异性的，可以进行高斯积分。

场点在S上时，对于平面单元，r⊥ns，故∂r/∂ns=0，则∂G/∂ns=∂G/∂r·∂r/∂ns=0，只存在弱奇异性。场点不在S上时，r⊥ns这个条件不成立，故∂G/∂ns≠0，则同时存在弱奇异性和强奇异性。实际上弱奇异积分比强奇异积分更容易收敛，所以只需保证强奇异积分收敛，便能消除积分方程的奇异性。

2 自适应细分方法

本文提出分级细分的自适应方法，并与其他自适应方法对比验证精度和收敛性。由于二维高斯积分是通过等参变换建立各种不规则单元与规则四边形单元的映射关系，之后在一个局部坐标系的正方形平面区域数值求解积分，所以文中统一采用顶点坐标为[-1, -1]、[1, -1]、[1, 1]和[-1, 1]的标准正方形单元来说明自适应细分方法，如图1所示。下文中所有省略的长度单位都是米。图1中黑点表示近场点在单元上的投影点，每行从左往右是自适应细分方案，第(1)行是平均细分方案，文献[14]在处理奇异积分时便是采取这种方案；第(2)～(5)行是投影点在不同位置时的分级细分方案。

平均细分方法是将原始单元平均细分为四个子单元，在下一次细分过程中每个子单元再平均细分为四个更小的子单元，如此递推下去将原始单元上的积分转移到各个子单元上去，直到全部子单元积分总和的误差达到收敛容许值为止。

分级细分方法是围绕近场点在单元上的投影点来细分，对原始单元进行平均细分后，只将投影点所在的子单元再次平均细分，其他子单元不变，如此下去直到收敛。由于投影点位置不同，于是产生了(2)～(5)这几种分级细分方案。

图1 自适应细分方案Fig.1. The adaptive subdivision processes

细分后，原始单元上的弱奇异积分和强奇异积分就变为子单元上的积分总和：

(7)

(8)

式中:s为原始单元，sj为子单元，N为子单元总数，Ng为高斯点数，Jj为子单元Jacobian，wu和wv为权系数。

设num为细分次数和子单元级别，则图1中五种细分方案的子单元总数为

(9)

由于平面单元的X、Y和Z坐标都是线性变化的，以X坐标为例说明如何根据自适应算法得到子单元坐标。对于平均细分方案，如图1中的第(1)行， 第一级子单元与原始单元的坐标关系如下：

(10)

其中T1～T4是坐标传递矩阵，X0是原始单元坐标：

(11)

于是通过同样的坐标传递，第二级子单元坐标为

(12)

如此继续下去，细分num次获得的所有子单元坐标为

(13)

式中:m1=1～4，m2=1～4。

对于分级细分方案，如图1中第(2)行，只对投影点所在单元进行平均细分，每次细分产生三个不含投影点的子单元，最后剩下一个投影点所在单元，故第num级子单元坐标为

(14)

图1中第(3)～(5)行的分级细分方案具有很高相似性，其子单元坐标的自适应算法都可参照第(2)行。由此便可清晰地构建各种情况下的自适应细分方法，不需要迭代和判断，只需要设定细分次数。

3 近场点的接近度

在进行收敛性分析之前，先要提出接近度的概念，并根据接近度决定单元划分的次数。由于近场点与单元的相对位置并不固定，不同的接近程度对拟奇异性的影响是有差别的。文献[25]初步指出拟奇异性是近场点与高斯积分点相互作用的结果，接近高斯点的区域会产生较强的奇异性，其他区域则很弱。为了深入发掘接近度对拟奇异性的影响，本文在标准单元周围取了若干近场点进行研究，如图2所示。由于正方形的对称性，只在八分之一直角三角形OAB附近取十个近场点，采用十点高斯积分，①～⑥号点投影在OAB的边上，⑦～⑩号点在平面OAB上，坐标如表1所示，其中t为场点到单元的距离。使每一近场点沿直线不断接近单元并考察收敛性：若频率无关的强奇异积分∫s-1/r2·∂r/∂nsds在三次细分以内收敛并且不细分的积分值与细分收敛后的积分值相对误差不超过1%，那么此时的距离称为“近远场”的临界值，超过这个临界距离就不存在拟奇异性，场点不再是近场点，单元也不需要细分就可正常积分。本文定义这种场点到单元的距离t与单元边长l的比值称为近场点的接近度E，即E=t/l。通过逼近计算得到场点到单元的临界距离t0和临界接近度E0如表2所示，其中⑦～⑩号点的∂r/∂ns=0，故这些点的强奇异积分以∫s-1/r2ds表示。

图2 近场点与标准单元的相对位置Fig.2 The relative positions of near-field points to the standard element

从表2看出①～⑥号点和⑦～⑩号点的临界距离t0和临界接近度E0各不相同，明显呈现出随着与O点距离增大而减小的趋势，表明场点离高斯积分区域越远，相互作用越弱，拟奇异性程度越低。

表1 近场点坐标

表2 临界距离和临界接近度

文献[17-22]采用的经验公式为

(15)

式中:A1为近场点到积分单元最小距离，A2为单元在积分方向上的长度，A3为该积分方向上的高斯点数，A4为奇异性阶次，A5为该积分方向上的高斯积分误差。对于标准单元取A2=2，A3=10，A4=2，A5=0.01，计算得到临界距离为A1=0.168 6，并不能满足任意位置近场点的误差要求，未体现出相对位置变化所带来的差异。

由于分级细分方法是围绕近场点在单元上的投影点来细分，实际造成场点与最小子单元的位置关系等同于⑤号点与标准单元的位置关系，所以采用十点高斯积分并且误差阈值为1%时，分级细分次数num可根据接近度(场点到单元距离与最小子单元边长之比)来决定，即t/(l/2num)≥0.04，于是得到：

(16)

式中l为原始单元平均边长。式(16)可作为拟奇异积分的单元细分准则。

4 收敛性分析

首先以标准单元(l=2)的收敛性分析验证分级自适应方法和单元细分准则的正确性。采用十点高斯积分，k=1 rad/m，取四个场点(编号1～4)，坐标分别为[0, 0, 0.01]、[1, 1, 0.01]、[0, 0, 0.001]和[1, 1, 0.001]，四个场点对应的强奇异积分分别表示为IR1、IR2、IR3和IR4，每次细分的积分值与最终收敛值的误差分别表示为er1、er2、er3和er4，收敛性曲线如图3所示。

图3 标准单元的积分收敛性曲线Fig.3 The convergence curves for the standard element

图3中只展示了IR实部，因为IR虚部几乎没有变化，可见强奇异性由积分的实部来体现，而弱奇异积分收敛非常快，同样不予展示。通过误差曲线er1和er2看出经过3次以上细分，误差就能降到1%以下，而误差曲线er3和er4显示经过7次以上细分，误差也降到1%以下。可见对于t分别为0.01和0.001的场点，收敛情况完全符合式(16)的细分准则。

若按照式(15)计算得到离场点最近子单元长度分别为0.118 6和0.011 9，需要细分5次和8次；若按照文献[23-24]中的“对角线准则”计算得到离场点最近子单元长度分别为0.008 2和0.000 82，需要细分8次和12次，都增加了额外的计算量。另外，式(15)和“对角线准则”需要考察每个子单元是否满足要求，单元的细分是一种反复迭代判断过程，而根据本文的单元细分依据式(16)可以进行预判，然后按本文的自适应方法一次性划分子单元和计算坐标，更加简捷实用。

再以宽为2，长分别4和8的两种矩形单元(分别简称为“矩四”和“矩八”)进行验证。令t=0.01，近场点5和6分别投影于“矩四”的中心和顶点；近场点7和8分别投影于“矩八”的中心和顶点。相应的奇异积分分别表示为IR5～8，每次细分的积分值与最终收敛值的误差分别表示为er5～8，收敛性曲线如图4所示。

图4 矩形单元的积分收敛性曲线Fig.4 The convergence curves for the rectangular elements

图4显示对于t为0.01的场点，分别经过4次和5次以上细分，误差也降到1%以下，收敛情况同样完全符合式(16)的细分要求。若按照式(15)的要求，需要细分6次和9次；若按照“对角线准则”，需要细分11次和13次，这样也增加了额外计算量。

值得注意的是，采用分级细分方法和平均细分方法得到的计算结果并无差别，这充分说明了分级细分方法的准确性，同时也表明单元上不同区域对拟奇异性的影响不同：越靠近场点的区域对奇异性的贡献越大，而远离节点的区域对奇异性的贡献迅速减小。常规高斯积分未考虑这种差别导致积分不准，平均细分方法平均化处理整个单元会增加很多不必要得计算量，而分级细分方法梯度化地对待远近区域，符合奇异性产生的物理本质，既保证了精度也节省了计算量。根据式(9)，两种细分方法产生的子单元总数如表3所示。

表3 两种细分方法产生的子单元总数

由以上积分收敛性分析可见，自适应分级细分方法能克服弱奇异性和强奇异性，有效计算拟奇异积分。同时依据本文所提的细分准则能够准确预判，不需要反复迭代判断，产生子单元数量更少。

另外，对于边长比例不协调、形状不规则的单元，以及场点靠近单元边角的情况，将远离场点的区域做细分是不必要的。为了很快地收敛，只需将投影点所在区域划分出来，得到一个比例较为规则的小块(阴影部分)，然后按照图1中分级细分方案进行细分，小块以外区域保持不变，如图5所示。

图5 特殊情况下的细分方法Fig.5 Some especial cases for subdivision

5 数值算例

5.1 预估拟奇异性

对于任意不规则单元，可用表2中临界接近度预估此种情况下是否存在拟奇异性。单元顶点和近场点的投影点坐标如图6所示，投影点落在不规则单元的一头。

图6 任意不规则单元和近场点的投影点Fig.6 The arbitrary irregular element and the subpoint

首先根据表2中②号点的临界接近度E0作判断：

t/l≥0.130

(17)

得到临界距离为0.333 0，其中l取为长短边的平均值，那么当近场点到单元距离t不小于0.333 0时，一定不存在拟奇异性。当t小于0.333 0时，将投影点所在部分划分出来(虚线将Y轴右半部分单元平分)，再根据表2中①号点的临界接近度E0作判断：

t/l≤0.150 5

(18)

得到临界距离为0.152 8，其中l取为投影点所在子单元的长短边平均值，那么当近场点到单元距离t不超过0.152 8时，一定存在拟奇异性。采用十点高斯积分，k=1 rad/m，t分别为0.333 0、0.3、0.2和0.152 8时，对应的拟奇异积分(分别为IR-1、IR-2、IR-3、IR-4)和误差(分别为Er-1、Er-2、Er-3、Er-4)随细分次数的收敛性曲线如图7所示。

图7 任意单元和近场点的拟奇异积分收敛性曲线Fig.7 The convergence curves of the nearly singular integrals on the arbitrary element

从图7看出，不细分时IR-1和IR-2是收敛的，IR-3和IR-4是不收敛的，表明t为0.333 0和0.3时不存在拟奇异性，t为0.2和0.152 8时存在拟奇异性。这符合前文中t不小于0.333 0时一定不存在拟奇异性和t不超过0.152 8时一定存在拟奇异性的预计，也显示了0.152 8

对于某些典型问题，利用临界接近度能很快判断是否存在拟奇异性。例如薄板声辐射计算中，厚度为多少才算是“薄板”值得探究。实际声学意义上的薄板并不是由板的长宽厚相对比例决定的，而是依据计算时单元尺寸与厚度之比。将图8中的矩形平板划分一致的矩形单元(不要求单元都为正方形)，那么上下表面的单元是重合的，节点作为近场点，其相对位置固定。采用十点高斯积分，根据计算时使用单元的类型和临界接近度得出声学薄板的临界厚度如表4所示，其收敛性已验证，l为单元长短边的平均值。厚度低于此临界值的板为声学薄板，需要处理拟奇异性。

图8 矩形薄板网格Fig.8 The rectangular thin plate grid

单元类型临界厚度/m常单元0．1505l线性单元0．04l二次单元0．1325l不连续线性单元0．11l不连续二次单元0．13l

5.2 脉动球声辐射

以水中脉动球声辐射为例验证方法的正确性，将球体划分242个四边形单元，如图9所示。半径a=1 m，c=1 500 m/s，ρ=1 000 kg/m3，vn=1 m/s，k=1 rad/m。沿半径方向距离球心R=1.001～1.1 m平均分布若干近场点，近场点与球体相对距离为D=(R-a)/a，与球心距离为R的场点声压解析解为

(19)

图9 球体网格模型Fig.9 The spherical grid

图10 两种边界元的计算对比曲线Fig.10. The results and errors with two methods

从图10看出，声压幅值曲线p2与p0十分吻合，p1与p0有一定差距，而且D越小差距越大，这一点在误差曲线Er1和Er2上可以更清晰地看到。在D的整个区间上，Er2都是小于Er1的，Er1变化较剧烈，而Er2相对平稳；Er1最大值超过18%，而Er2一直在1%左右。这些情况说明常规边界元适用于计算远场声压，计算近场声压则难以保证精度，本文的分级细分自适应边界元法能有效计算近场声压。二者计算时间基本相同，因为对靠近场点的单元进行细分的计算量相对于总体计算量十分微小。对于声学分析，采用平均细分法、式(15)或“对角线准则”所得到的场点声压与本文分级细分法计算结果并无差别，却会产生更多子单元和额外计算量。

5.3 回转壳声辐射

如图11所示，回转壳由半球壳、圆柱壳和圆锥壳组成，半球体半径为1 m，圆柱体长7 m，圆锥体高2 m，上端面半径为0.16 m，表面统一振速为vn=1×10-5m/s，其它参数同上。分别沿半球、圆柱和圆锥表面单元中心法线设置3条近场点链，距离R为0.001～2 m。三个单元平均边长l分别为0.29、0.39和0.26 m，故满足式(16)的num皆为4。分别用常规高斯积分方法和自适应积分方法计算这些近场点声压，结果如图12所示，其中eR是常规方法计算结果与自适应方法计算结果的误差绝对值。

自适应方法中num取4和6的计算结果相等，说明num为4时结果已经收敛，图12中仅展示num取4的曲线。可看出两种积分方法的计算结果曲线的近场差别比脉动球算例更加显著， 而且eR降到约1%所需的R也更大，这个现象说明常规积分方法的近场不精确程度也会受到结构尺寸的影响。

图11 回转壳网格模型Fig.11 The rotative shell grid

图12 3条近场点链声压和误差曲线Fig.12 The sound pressures and absolute errors of the 3 chains of near-field points

6 结 论

对于边界元中的拟奇异积分问题，本文鉴于自适应精确积分的思想提出一种自适应分级细分方法，不需要迭代和判断就能够准确有效地计算任意近场点的拟奇异积分，所需子单元更少，且精度可以控制，其灵活性和高效性是一些已有方法所不具备的。接近度的研究为拟奇异积分的计算和预估拟奇异性程度提供理论基础，具有一定的实用价值。

收敛性分析发现：单元上靠近场点的区域对拟奇异性的影响远远大于远离场点的区域，梯度化地对待远近区域能较好地考虑奇异性的机理，节省计算量。通过对接近度的研究发现：拟奇异性强弱由场点与单元的相对位置决定，并非单纯由近场点到单元距离决定。

数值分析表明：在预估任意不规则单元和近场点是否存在拟奇异性时，临界接近度可作为判断拟奇异性存在与否的充分条件而非必要条件。对于薄板声辐射问题，声学意义上的薄与厚并不是由板的尺寸相对比例决定的，而要以接近度为参考，不同类型的计算单元所对应的临界厚度也不相同。常规边界元积分方法近场计算精度较低，其不精确程度也会受到结构尺寸的影响。

[1] 李善德, 黄其柏, 李天匀. 新的对角形式快速多极边界元法求解声学 Helmholtz 方程[J]. 物理学报, 2012,61(6): 64301. LI Shande, HUANG Qibai, LI Tianyun. A new diagonal form fast multipole boundary element method for solving acoustic Helmholtz equation[J]. Acta Physica Sinica, 2012, 61(6): 64301.

[2] 白杨, 汪鸿振. 声学-结构设计灵敏度分析[J]. 振动与冲击,2003, 22(3):43-45. BAI Yang, WANG Hongzhen. Acoustic-structural design sensitivity analysis[J]. Journal of Vibration and Shock, 2003,22(3):43-45.

[3] 杨迎春, 周其斗. 运动介质中奇异边界元积分式的精确求解[J]. 振动与冲击, 2011, 30(3): 242-245. YANG Yingchun, ZHOU Qidou. Accurate calculation of weak singular integrals in boundary element method in moving flows[J]. Journal of Vibration and Shock, 2011, 30(3): 242-245.

[4] 张耀明, 孙翠莲, 谷岩. 边界积分方程中近奇异积分计算的一种变量替换法[J]. 力学学报, 2008, 40(2): 207-214. ZHANG Yaoming, SUN Cuilian, GU Yan. The evaluation of nearly singular integrals in the boundary integral equations with variable transformation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(2): 207-214.

[5] 牛忠荣, 张晨利, 王秀喜. 边界元法分析狭长体结构[J]. 计算力学学报, 2003, 20(4): 391-396. NIU Zhongrong, ZHANG Chenli, WANG Xiuxi. Boundary element analysis of the thin-wall structures[J]. Chinese Journal of Computational Mechanics, 2003, 20(4): 391-396.

[6] TANAKA M, SLADEK V, SLADEK J. Regularization techniques applied to boundary element methods[J]. Applied Mechanics Reviews,1994, 47(10): 457-499.

[7] ZHANG Y, GU Y, CHEN J. Stress analysis for multilayered coating systems using semi-analytical BEM with geometric non-linearities[J]. Computational Mechanics, 2011, 47(5): 493-504.

[8] NIU Z, CHENG C, ZHOU H, et al. Analytic formulations for calculating nearly singular integrals in two-dimensional BEM[J]. Engineering Analysis with Boundary Elements, 2007, 31(12): 949-964.

[9] ZHOU H, NIU Z, CHENG C, et al. Analytical integral algorithm in the BEM for orthotropic potential problems of thin bodies[J]. Engineering Analysis with Boundary Elements, 2007, 31(9):739-748.

[10] ZHOU H, NIU Z, CHENG C, et al. Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for an isotropicpotential problems[J]. Computers & Structures, 2008, 86(1516):1656-1671.

[11] NIU Z, ZHOU H. The natural boundary integral equation in potential problems and regularization of the hypersingular integral[J]. Computers & Structures, 2004, 82(2):315-323.

[12] DAVIES T G, GAO X. Three-dimensional elasto-plastic analysis via the boundary element method[J]. Computers and Geotechnics,2006, 33(3):145-154.

[13] KITA E, KAMIGA N. Error estimation and adaptive mesh refinement in boundary element method, an overview[J]. Engineering Analysis with Boundary Elements, 2001,25(7):479-495.

[14] WEI Y, WANG Y, CHANG S, et al. Numerical prediction of propeller excited acoustic response of submarine structure based on CFD, FEM and BEM[J]. Journal of Hydrodynamics, 2012, 24(2): 207-216.

[15] CROAKER P, KESSISSOGLOU N, MARBURG S. Strongly singular and hypersingular integrals for aeroacoustic incident fields[J]. International Journal for Numerical Methods in Fluids, 2015, 77(5): 274-318.

[16] DONG Y, ZHANG J, XIE G, et al. A general algorithm for the numerical evaluation of domain integrals in 3D boundary element method for transient heat conduction[J]. Engineering Analysis with Boundary Elements, 2015, 51: 30-36.

[17] DAVIES T G, GAO X. Three-dimensional elasto-plastic analysis via the boundary element method[J]. Computers and Geotechnics, 2006, 33(3):145-154.

[18] GAO X W, YANG K, WANG J. An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals[J]. Engineering Analysis with Boundary Elements, 2008, 32(8): 692-696.

[19] LI Y, YE J H. The interaction between membrane structure and wind based on the discontinuous boundary element[J]. Science China Technological Sciences, 2010, 53(2): 486-501.

[20] 王静, 高效伟. 基于单元子分法的结构多尺度边界单元[J].计算力学学报, 2010, 27(2): 258-263. WANG Jing, GAO Xiaowei. Structural multi-scale boundary element method based on element subdivision technique[J]. Chinese Journal of Computational Mechanics, 2010, 27(2): 258-263.

[21] QU S, CHEN H, LI S. Adaptive integration method based on sub-division technique for nearly singular integrals in near-field acoustics boundary element analysis[J]. Journal of Low Frequency Noise, Vibration and Active Control, 2014,33(1): 27-46.

[22] GAO X W,DAVIES T G. Adaptive integration in elasto-plastic boundary element analysis[J]. Journal of the Chinese Institute of Engineers, 2000, 23(3):349-356.

[23] ZHANG J, TANAKA M, MATSUMOTO T. Meshless analysis of potential problems in three dimensions with the hybrid boundary node method[J]. International Journal for Numerical Methods in Engineering, 2004, 59(9): 1147-1166.

[24] ZHANG J, QIN X, HAN X, et al. A boundary face method for potential problems in three dimensions[J]. International Journal for Numerical Methods in Engineering, 2009, 80(3): 320-337.

[25] KANG S C, IH J G. On the accuracy of nearfield pressure predicted by the acoustic boundary element method[J]. Journal of Sound and Vibration, 2000, 233(2): 353-358.

Adaptive method for calculating nearly singular integrals in acoustic boundary elements

MIAO Yuyue1,2, LI Tianyun1,2,3, ZHU Xiang1,2, ZHANG Guanjun1,2, GUO Wenjie1,2

(1. School of Naval Architecture and Ocean Engineering，Huazhong University of Science & Technology, Wuhan 430074,China;2.Hubei Key Laboratory of Naval Architecture & Ocean Engineering Hydrodynamics, Wuhan 430074,China;3.National Key Laboratory on Ship Vibration & Noise, Wuhan 430064,China)

An adaptive method was presented to calculate nearly singular integrals in acoustic boundary elements. The element was subdivided into subelements hierarchically so as to remove the near singularity by transforming the integral over the initial element to Gauss integrals over subelements. Based on this method, the near singularity was thoroughly studied and the concept of proximity was described. The critical proximity can be used as the criterion for calculating nearly singular integrals and to predict the existence of near singularity. Compared with previous methods, the adaptive method is more flexible and efficient. The integral precision can be adjusted and is not restricted by the positions of near-field points. It is found that the positional relationships between near-field points and the element have important effect on nearly singular integrals. The work provides more understanding and choice to the problem of near singularity in acoustic boundary elements.

boundary element; nearly singular intergral; adaptivity; proximity

国家自然科学基金(51379083;51479079);高等学校博士学科点专项科研基金(20120142110051)

2015-09-09 修改稿收到日期：2015-12-05

缪宇跃 男，博士生，1988年生

李天匀 男，教授，博士生导师，1969年生 E-mail:ltyz801@mail.hust.edu.cn

TB532

A

10.13465/j.cnki.jvs.2017.02.004