引 言
边界元法(BEM)是一种以边界积分方程(BIE)为理论基础，结合数值离散技术的高效工程数值方法．由于积分与插值都只需要针对边界变量实施，BEM被认为是一类降维方法，亦即三维问题只需要对二维边界进行离散，其计算规模相较于以有限元法为代表的域离散类工程数值方法更小[1-3]．同时，由于在方程导出过程中采用了两次分步积分，对试函数的连续性要求降至C0连续，对非连续单元的适用性好，使得其网格离散十分灵活便捷[4-6]．边界积分方程采用了物理问题的基本解，被认为是一种半解析的方法，实施时使用较少的计算代价就能获得较高的计算精度[7-9]．
BEM中基本解的奇异性保证了系数矩阵良好的性态，使其特别适合于裂纹问题分析[10-11]．同时，由于基本解的奇异性，边界积分方程在无穷远处自然满足声学问题中的Summerfield条件，因此BEM也非常适合用于求解声学等无限域问题，分析过程中无需使用人工边界等[7，12-13]．然而也正是因为基本解的奇异性，BEM在实施过程中，会遇到诸如奇异积分计算难题和近奇异积分计算难题[14-21]．BEM中的配置点往往设置在边界节点上，因此，边界积分属于瑕积分，瑕点处的积分为奇异积分，为避免直接计算奇异积分，在稳态问题分析中，常采用一种正则化方法间接计算奇异积分[14-15]．而在瞬态问题中，奇异积分的计算必须辅以非线性变换之后再实施数值积分[16-17]．当前，BEM中奇异积分计算难题已经得到较好的解决，而对近奇异积分计算则还未出现为研究人员公认的较好解决方案．
在应用BEM分析薄型乃至超薄型结构时，由于边界积分方程中一般积分项为
(1)
式中y表示配置点，x表示被积点，S为边界区域，r(x,y)表示两点的距离，f为光滑函数，p为r的指数，由相应物理问题的基本解决定．在应用边界积分方程分析薄型结构时，大量存在y点和积分单元接近的情况，当积分单元面积较大时，积分核函数在积分单元内的分布相对陡峭，如果采用常规数值积分方案，则需将积分单元细分以保证在积分单元内积分核函数值变化较小．如此一来，计算量将增加．通过细分积分单元来提高近奇异积分计算的精度是一种数值稳定的方法，其实质是通过增加积分点数量以提升计算精度．这种方法虽稳定，但牺牲了计算效率，在分析薄型特别是超薄型结构时，近奇异积分计算将耗费大量计算时间．而计算近奇异积分的另一种思路则是通过调整积分点的位置使其适合基本解的分布形式，进而以少量的积分点获得较高的积分计算精度．由于基本解是关于位置的非线性函数，因此，常应用非线性变换将积分空间映射到另一空间，在新的参数空间中再进行常规数值积分[18-23]．
计算近奇异积分的典型非线性变换方法包括角度变换[22]、指数变换[18]、多项式变换[23]、sinh变换[24]以及距离变换[25]．这些变换各有特点，其中距离变换被认为是一种非常有效且通用性好的变换方法．距离变换在实施的过程中，首先需要采用极坐标变换将积分空间映射至极坐标空间，再通过定义极坐标空间中的距离函数以实施距离变换．值得注意的是，在极坐标空间中所定义的距离函数是距离空间两个独立变量的函数，而距离变换通常在单个变量(ρ)上实施．然而，研究发现，当配置点位于积分单元边界附近时，积分核函数除了在ρ方向具有近奇异性之外，在θ方向同样具有近奇异性．对极坐标空间的θ变量进一步采用角度变换，消除积分核在该方向的近奇异性能进一步提升近奇异积分计算精度．
为高效准确地计算近奇异积分，本文在对参数空间中的距离变量实施距离变换的基础上，再对角度变量实施角度变换．通过角度-距离复合变换，使积分点在实际空间中更加聚拢于源点附近，进而使插值型积分公式能够更好地拟合积分核函数值在源点附近区域的剧烈变化情况，从而获得更高的数值积分精度．两个数值算例分别关注平面三角形积分单元与曲面三角形积分单元上的近奇异积分，算例结果验证了角度-距离复合变换法计算近奇异积分时的有效性和准确性．
1 近奇异积分
以三维稳态热传导问题的边界积分方程为例，其积分项为[3,5,25]
(2)
其配置点y和积分单元S的位置关系如图1所示，图中xc表示积分单元中距离配置点y最近的点，n表示积分单元内xc处的外法线方向．当r0与S的面积比值较小时，随着x在S内移动，1/r的值变化较大，不利于数值积分．Gauss积分是基于有理多项式插值的内插型数值积分方案，被认为在插值型积分公式中具有最高的代数精度，在一维情况下k个积分点理论可以达到2k+1次代数精度，即对次数不超过2k+1次多项式函数积分可以获得精确积分值．采用Gauss积分计算边界积分公式(2)，则
(3)
式中xi表示第i个Gauss积分点, ωi表示第i个积分点对应的权系数，Ji表示正则Gauss积分空间到实际空间转换的Jacobi在第i个积分点上的值．
图1 配置点与积分单元
Fig.1 The collocation point and the integral element
由式(2)可以看出，当配置点位于积分区域外时，积分核函数是光滑连续的，Weierstrass函数逼近定理保证了任意闭区间上的连续函数都能使用有理多项式级数以任意精度逼近．然而当配置点距离积分区域很近时，由于其核函数在积分区域内(特别是最近点附近)变化率较高，采用有理多项式逼近时，同样阶的余项值更大，因此需使用更高阶多项式才能得到足够的逼近精度，对其直接进行Gauss积分时，需要使用更多积分点．
直接使用Gauss积分不能在不降低计算效率的基础上获得理想的计算精度，众多学者通过借鉴边界元奇异积分计算的思路，通过非线性变换消除奇异性，提出了许多适合消去近奇异性的非线性变换方案，其中距离变换被认为是较理想的方案之一．
2 极坐标变换
在进行近奇异性消去之前，Gauss积分需要在正则空间中进行，边界积分单元往往需要变换到正则参数空间．以图2(a)所示的三角形积分单元为例，其正则参数空间通常为参数化正则三角形(如图2(b)所示)．
图2 三角形积分单元正则化
Fig.2 Normalization of a triangular element
实际空间中距离r可在正则化参数空间进行表达，如图1所示配置点与积分单元的相对位置，距离r在正则化参数空间中的表达式为
O((s-s0)2+(t-t0)2),
(4)
其中，下标k表示实际空间中k方向的分量，(s0,t0)是最近点xc在参数空间中的坐标．在参数空间(s, t)中，进一步对积分核实施极坐标变换：
(5)
则积分可表示为
(6)
其中f, f1为光滑函数．此时，将变换(5)代入展开式(4)，得到r在极坐标空间的表达式为
(7)
(8)
其中Ak(θ)是只与角度θ相关的项，
略去式(8)中关于ρ的高阶项并将其代入式(6)可以得到
(9)
值得注意的是，积分式(9)不仅有关于ρ的奇异性，关于θ同样存在奇异性，而关于θ的奇异性在最近点xc靠近积分单元边界时表现得很明显，采用距离变换可对关于ρ的奇异性进行处理．
3 角度-距离组合变换
在进行距离变换之前，需定义参数空间中的距离函数：
(10)
对式(9)中变量ρ做如下距离变换：
(11)
通过该变换，积分(9)可进一步转化为
(12)
此时，积分核关于ζ已经不再具有奇异性．当最近点落在积分单元内部时，距离变换能够消除整体积分核的近奇异性．然而，更一般的情况是，最近点落在积分单元的边界附近或边界上，而非落在积分单元内部．在这种情况下，积分核除了ρ方向存在近奇异性外，在θ方向也存在近奇异性．因此，对θ方向的近奇异性也可以采用非线性变换进行处理．
针对θ方向的近奇异性，可采用一种角度变换进行处理．针对θ变量，采用如下角度变换：
ϑ
(13)
其中h为正则化参数空间中最近点到积分单元边界之间的最近距离，α是正则化参数空间中过最近点到其最近边界的垂线相对于横轴的倾角．通过式(13)所示角度变换，关于θ方向的近奇异性可消除．
值得注意的是，本文中出现的最近点xc并非等于配置点y在积分单元的投影点，当投影落在积分单元外时，xc指积分单元中距离配置点的最近点(往往落在积分单元的边界上)，此时y与xc连线方向与xc上法向n并不平行．在此情况下，式(4)仍成立，只是后续在构造非线性变换时需根据最近点的不同位置进行特殊构造．针对距离变换的构造方法可以参考文献[26]，该文献对此投影点在积分单元中不同位置的情况进行了详细讨论并给出了验证算例．
4 数 值 算 例
为验证角度-距离复合变换的计算精度，算例中分别考虑最近点在积分单元内的不同位置、配置点距积分单元的不同距离以及式(14)和式(15)两种近奇异积分：
(14)
(15)
算例中考虑20种不同位置的配置点分别对应4个不同的最近点，为简单起见，算例只关注最近点落在积分单元内部的情况，假设配置点的位置可表示为
y=xc+r0n．
(16)
两个算例中数值积分分别在线性三角形单元和二次三角形单元上进行．为评估计算精度，相对误差可表示为
(17)
在综合评估非线性变换效果的过程中，式(18)所示配置点到积分单元的最近距离与积分单元面积的比值是一个非常重要的参数．
λ=r0/Sarea．
(18)
在变换实施后，每个方向采用5个Gauss积分点亦即单元内一共25个Gauss积分点计算数值积分．
4.1 三角形线性单元上的近奇异积分
考虑标准三角形单元， 三个顶点空间坐标分别为(0,0,0)，(1,0,0)和(0,1,0)， 配置点对应的4个最近点xc1， xc2， xc3以及xc4在积分平面参数空间中坐标分别为(0.1,0.1)，(0.01,0.01)，(0.3,0.4)以及(0.1,0.89)．4个最近点中，xc3最靠近参数空间中单元的几何中心，xc2靠近积分单元顶点．
表1～4列出了采用角度-距离复合变换和单独采用距离变换两种情况下近奇异积分的计算误差．
由表1～4可以看出，计算I1时应用复合变换得到的结果精度相较于单独使用距离变换得到的结果精度要高．然而在最近点落在积分单元中心附近时，单独使用距离变换亦可以获得较高计算精度．而在最近点靠近积分单元边界时，复合变换获得的计算精度较单独距离变换获得的计算精度高出3个以上的数量级．结果表明，在参数空间中最近点靠近积分单元边界时，积分核在环向也具有近奇异性，通过对环向变量进行角度变换之后，能有效提升计算精度．同时，也可以发现，在计算I1时，计算结果对特征比例大小并不敏感．
表1 最近点为xc1时两种计算方法计算I1的相对误差
Table 1 Relative errors of two methods in the computation of I1 in the case of xc1
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation7.719E-71.869E-61.901E-61.886E-61.875E-6combined transformation6.055E-102.015E-83.109E-81.631E-84.931E-9
表2 最近点为xc2时两种计算方法计算I1的相对误差
Table 2 Relative errors of two methods in the computation of I1 in the case of xc2
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation6.111 9E-73.001E-59.718E-41.128E-31.144E-3combined transformation9.889E-102.638E-82.844E-81.315E-84.168E-9
表3 最近点为xc3时两种计算方法计算I1的相对误差
Table 3 Relative errors of two methods in the computation of I1 in the case of xc3
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation4.814E-81.240E-71.232E-71.411E-71.529E-7combined transformation1.857E-101.894E-83.295E-81.642E-84.722E-9
表4 最近点为xc4时两种计算方法计算I1的相对误差
Table 4 Relative errors of two methods in the computation of I1 in the case of xc4
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation2.736E-53.020E-33.460E-33.440E-33.437E-3 combined transformation2.656E-88.180E-82.475E-87.913E-93.996E-9
为考察复合变换法计算更高阶近奇异积分的效果，表5～8列出了采用两种方法在不同最近点情况下计算I2时的计算误差．
表5 最近点为xc1时两种计算方法计算I2的相对误差
Table 5 Relative errors of two methods in the computation of I2 in the case of xc1
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation4.443E-67.388E-61.784E-55.100E-51.090E-4combined transformation1.156E-81.504E-61.331E-54.711E-51.054E-4
表6 最近点为xc2时两种计算方法计算I2的相对误差
Table 6 Relative errors of two methods in the computation of I2 in the case of xc2
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation2.616E-67.041E-34.036E-25.379E-26.039E-2combined transformation3.010E-82.332E-61.145E-53.301E-57.475E-5
表7 最近点为xc3时两种计算方法计算I2的相对误差
Table 7 Relative errors of two methods in the computation of I2 in the case of xc3
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation2.382E-78.135E-71.510E-55.372E-51.169E-4combined transformation3.614E-91.575E-61.606E-55.477E-51.180E-4
表8 最近点为xc4时两种计算方法计算I2的相对误差
Table 8 Relative errors of two methods in the computation of I2 in the case of xc4
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation3.411E-43.334E-21.886E-29.524E-34.702E-3combined transformation2.718E-73.249E-61.247E-54.509E-59.500E-5
从表5～8数据来看，在计算I2时复合变换效果略好于单独距离变换，但差距不如计算I1时明显．同时，特征比例对计算结果的影响大于计算I1时的情况．对比表1～8，在相同最近点情况下， I1计算误差低于I2计算误差，与常规情况一致．
4.2 三角形二次单元上的近奇异积分
第二个算例中考虑曲面单元上的近奇异积分计算情况，4个最近点在参数空间中的位置与线性单元算例一致．本算例关注的二次单元如图3所示．二次三角形单元6个节点的空间坐标分别为：(4, 0, 0.2)，(0, 4 ,0)，(0, 0, 0.5)，(2, 2, 0.1)，(0, 2, 0.2)以及(2, 0, 0.3)．
图3 二次三角形单元
Fig.3 The quadratic triangle unit
表9～16列出了两种方法在不同最近点情况下计算I1和I2的相对误差情况．
表9 最近点为xc1时两种计算方法计算二次单元I1的相对误差
Table 9 Relative errors of two methods in the computation of I1 over the quadratic element in the case of xc1
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation6.986E-102.554E-63.952E-64.104E-64.121E-6combined transformation5.071E-101.124E-102.217E-98.623E-108.055E-10
表10 最近点为xc2时两种计算方法计算二次单元I1的相对误差
Table 10 Relative errors of two methods in the computation of I1 over the quadratic element in the case of xc2
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation9.407E-108.036E-54.397E-49.347E-49.941E-4combined transformation8.776E-106.202E-103.414E-93.754E-102.515E-10
表11 最近点为xc3时两种计算方法计算二次单元I1的相对误差
Table 11 Relative errors of two methods in the computation of I1 over the quadratic element in the case of xc3
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation1.710E-101.104E-71.575E-71.659E-71.618E-7combined transformation4.781E-111.313E-104.299E-113.158E-91.462E-9
表12 最近点为xc4时两种计算方法计算二次单元I1的相对误差
Table 12 Relative errors of two methods in the computation of I1 over the quadratic element in the case of xc4
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation3.046E-96.959E-43.401E-33.419E-33.409E-3combined transformation5.205E-93.288E-86.251E-82.189E-88.978E-10
从表9～12可以看出，在计算I1时，除λ=1外,复合变换表现更优，在最近点靠近积分单元边界时，其精确度提升大约4～6个数量级，在最近点靠近积分单元中心时，其精确度提升大约2个数量级．说明在二次积分单元内采用复合变换计算近强奇异积分，能切实提升计算精确度．表13～16列出了二次单元上的近超奇异积分计算情况．
表13 最近点为xc1时两种计算方法计算二次单元I2的相对误差
Table 13 Relative errors of two methods in the computation of I2 over the quadratic element in the case of xc1
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation2.365E-91.320E-52.054E-52.297E-52.406E-5combined transformation8.742E-108.645E-92.179E-74.523E-74.954E-7
表14 最近点为xc2时两种计算方法计算二次单元I2的相对误差
Table 14 Relative errors of two methods in the computation of I2 over the quadratic element in the case of xc2
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation2.122E-91.437E-32.116E-24.311E-25.229E-2combined transformation1.466E-93.010E-82.194E-71.560E-64.759E-6
表15 最近点为xc3时两种计算方法计算二次单元I2的相对误差
Table 15 Relative errors of two methods in the computation of I2 over the quadratic element in the case of xc3
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation9.091E-105.451E-76.450E-72.001E-61.018E-6combined transformation1.038E-108.042E-102.424E-79.738E-78.408E-8
表16 最近点为xc4时两种计算方法计算二次单元I2的相对误差
Table 16 Relative errors of two methods in the computation of I2 over the quadratic element in the case of xc4
methodλ=0.1λ=0.01λ=0.001λ=0.000 1λ=0.000 01distance transformation2.866E-88.638E-32.878E-21.400E-27.011E-3combined transformation2.935E-87.965E-73.877E-67.290E-61.205E-5
从表13～16可以看出，在计算近超奇异积分I2时，除λ=1外,复合变换的表现全面优于距离变换的表现，采用复合变换较之单独距离变换计算精度提升约2～4个数量级．此外还可以看到，计算误差随着配置点靠近积分单元不断增大．
5 结 论
在应用BEM分析薄型结构上的物理问题时，当配置点靠近积分单元时，边界积分方程中的边界积分在极坐标参数空间中径向和环向都具有近奇异性，环向的近奇异性在配置点距积分单元最近点靠近积分单元边界时表现明显．通过在极坐标空间中，分别对径向变量实施距离变换、环向变量实施角度变换，能有效消除两个方向的近奇异性，较之单独对径向使用距离变换的情况计算精度有大幅提升．
[1] 高锁文, 汪越胜, 章梓茂, 等.含孔薄板弯曲波动的双互易边界元法[J].应用数学和力学, 2005, 26(12): 1417-1424.(GAO Suowen, WANG Yuesheng, ZHANG Zimao, et al.Dual reciprocity boundary element method for flexural waves in thin plate with cutout[J].Applied Mathematics and Mechanics, 2005, 26(12): 1417-1424.(in Chinese))
[2] 张耀明, 吕和祥, 王利民.位势平面问题的新的规则化边界积分方程[J].应用数学和力学, 2006, 27(9): 1017-1022.(ZHANG Yaoming, LÜ Hexiang, WANG Limin.Novel regularized boundary integral equations for potential plane problems[J].Applied Mathematics and Mechanics, 2006, 27(9): 1017-1022.(in Chinese))
[3] 周枫林, 李光, 孙晓, 等.水坝瞬态热传导分析中的拟初始条件边界元法[J].计算力学学报, 2016, 33(6): 826-833.(ZHOU Fenglin, LI Guang, SUN Xiao, et al.Application of quasi-initial condition boundary element method in transient heat conduction problem on gravity dams[J].Chinese Journal of Computational Mechanics, 2016, 33(6): 826-833.(in Chinese))
[4] 汪攀, 张见明, 韩磊, 等.基于带约束前沿推进的四边形网格生成方法[J].湖南大学学报(自然科学版), 2017, 44(8): 29-34.(WANG Pan, ZHANG Jianming, HAN Lei, et al.An advancing front quadrilateral mesh generation method with constraint[J].Journal of Hunan University(Natural Sciences), 2017, 44(8): 29-34.(in Chinese))
[5] ZHANG J M, LIN W C, DONG Y Q, et al.A double-layer interpolation method for implementation of BEM analysis of problems in potential theory[J].Applied Mathematical Modelling, 2017, 51: 250-269.
[6] ZHANG J M, SHU X M, TREVELYAN J, et al.A solution approach for contact problems based on the dual interpolation boundary face method[J].Applied Mathematical Modelling, 2019, 70: 643-658.
[7] 王守信, 刘喜平, 彭天国, 等.求解变系数非齐次亥姆霍茨方程的边界单元法[J].应用数学和力学, 1996, 17(1): 81-85.(WANG Shouxin, LIU Xiping, PENG Tianguo, et al.The BEM for solving the nonhomogeneous Helmholtz equation with variable coefficients[J].Applied Mathematics and Mechanics, 1996, 17(1): 81-85.(in Chinese))
[8] 丁睿, 朱正佑, 程昌钧.粘弹性薄板动力响应的边界元方法(Ⅰ)[J].应用数学和力学, 1997, 18(3): 211-216.(DING Rui, ZHU Zhengyou, CHENG Changjun.Boundary element method for solving dynamical response of viscoelastic thin plate(Ⅰ)[J].Applied Mathematics and Mechanics, 1997, 18(3): 211-216.(in Chinese))
[9] 马杭, 郭钊, 秦庆华.二维多项式本征应变边界积分方程及其数值验证[J].应用数学和力学, 2011, 32(5): 522-532.(MA Hang, GUO Zhao, QIN Qinghua.Two-dimensional polynomial eigenstrain formulation of boundary integral equation with numerical verification[J].Applied Mathematics and Mechanics, 2011, 32(5): 522-532.(in Chinese))
[10] XIE G Z, ZHANG J M, HUANG C, et al.A direct traction boundary integral equation method for three-dimension crack problems in infinite and finite domains[J].Computational Mechanics, 2014, 53(4): 575-586.
[11] ZHANG J M, DONG Y Q, LIN W C, et al.A singular element based on dual interpolation BFM for V-shaped notches[J].Applied Mathematical Modelling, 2019, 71: 208-222.
[12] WANG X H, ZHANG J M, ZHOU F L, et al.An adaptive fast multipole boundary face method with higher order elements for acoustic problems in three-dimension[J].Engineering Analysis With Boundary Elements, 2013, 37(1): 114-152.
[13] 吴海军, 蒋伟康, 刘轶军.基于Burton-Miller边界积分方程的二维声学波动问题对角形式快速多极子边界元及其应用[J].应用数学和力学, 2011, 32(8): 920-933.(WU Haijun, JIANG Weikang, LIU Yijun.Diagonal form fast multipole boundary element method for 2D acoustic problems based on Burton-Miller BIE formulation and its applications[J].Applied Mathematics and Mechanics, 2011, 32(8): 920-933.(in Chinese))
[14] LIU Y J.On the simple-solution method and non-singular nature of the BIE/BEM: a review and some new results[J].Engineering Analysis With Boundary Elements, 2000, 24(10): 789-795.
[15] LIU Y J, RUDOLPHI T J.Some identities for fundamental solutions and their applications to weakly-singular boundary element formulations[J].Engineering Analysis With Boundary Elements, 1991, 8(6): 301-311.
[16] 高效伟, 冯伟哲, 杨恺.边界元中计算任意高阶奇异线积分的直接法[J].力学学报, 2014, 46(3): 428-435.(GAO Xiaowei, FENG Weizhe, YANG Kai.A direct method for evaluating line integrals with arbitrary high order of singularities[J].Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(3): 428-435.(in Chinese))
[17] 李俊, 冯伟哲, 高效伟.一种基于直接计算高阶奇异积分的断裂力学双边界积分方程分析法[J].力学学报, 2016, 48(2): 387-398.(LI Jun, FENG Weizhe, GAO Xiaowei.A dual boundary integral equation method based on direct evaluation of higher order singular integral for crack problems[J].Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2): 387-398.(in Chinese))
[18] XIE G Z, ZHANG J M, DONG Y Q, et al.An improved exponential transformation for nearly singular boundary element integrals in elasticity problems[J].International Journal of Solids and Structures, 2014, 51(6): 1322-1329.
[19] 孙锐, 胡宗军, 牛忠荣, 等.边界元法计算声辐射时几乎奇异积分的处理方法[J].中国科学: 物理学 力学 天文学, 2017, 47(9): 22-31.(SUN Rui, HU Zongjun, NIU Zhongrong, et al.A method of treating the nearly singular integral in calculation of sound radiation with BEM[J].Scientia Sinica Physica, Mechanica & Astronomica, 2017, 47(9): 22-31.(in Chinese))
[20] 胡宗军, 牛忠荣, 程长征, 等.薄体结构温度场的高阶边界元分析[J].应用数学和力学, 2015, 36(2): 149-158.(HU Zongjun, NIU Zhongrong, CHENG Changzheng, et al.High-order boundary element analysis of temperature fields in thin-walled structures[J].Applied Mathematics and Mechanics, 2015, 36(2): 149-158.(in Chinese))
[21] XIE G Z, ZHOU F L, ZHANG J M, et al.New variable transformations for evaluating nearly singular integrals in 3D boundary element method[J].Engineering Analysis With Boundary Elements, 2013, 37(9): 1169-1178.
[22] HAYAMI K.Variable transformations for nearly singular integrals in the boundary element method[J].Publications of the Research Institute for Mathematical Sciences, 2005, 41(4): 821-842.
[23] TELLES J C F.A self-adaptive coordinate transformation for efficient numerical evluation of general boundary element integrals[J].International Journal for Numerical Methods in Engineering, 1987, 24(5): 959-973.
[24] JOHNSTON B M, JOHNSTON P R, ELLIOTT D.A sinh transformation for evaluating two-dimensional nearly singular boundary element integrals[J].International Journal for Numerical Methods in Engineering, 2007, 69(4): 1460-1479.
[25] MA H, KAMIYA N.A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two- and three-dimensional elasticity[J].Computational Mechanics, 2002, 29(4): 277-288.
[26] XIE G Z, ZHOU F L, ZHANG J M, et al.New variable transformations for evaluating nearly singular integrals in 3D boundary element method[J].Engineering Analysis With Boundary Elements, 2013, 37(9): 1169-1178.