消除声学间接边界元非唯一解的方法的对比研究_梁梦辉.pdf

131梁梦辉等消除声学间接边界元非唯一解的方法的对比研究第 2 期第 42 卷 第 2 期2023 年 4 月内蒙古工业大学学报（自然科学版）Journal of Inner Mongolia University of Technology（Natural Science Edition）Vol.42 No.2Apr.2023文章编号：1001-5167（2023）02-0131-06消除声学间接边界元非唯一解的方法的对比研究梁梦辉1，陈红永2，郑昌军1，张永斌1（1.合肥工业大学 噪声振动工程研究所，合肥 230000；2.中国工程物理研究院总体工程研究所，四川 绵阳 621999）Comparison of the methods for eliminating the non-unique solution in acoustic Indirect boundary element methodLIANG Menghui1，CHEN Hongyong2，ZHENG Changjun1，ZHANG Yongbin1收稿日期：2022-11-11基金项目：国家自然科学基金项目（11872168）第一作者：梁梦辉（1993），男，博士研究生，主要从事声学边界元法及其工程应用的研究。E-mail:mh_通信作者：郑昌军（1983），男，教授，主要从事声学边界元法及其工程应用的研究。E-mail:(1.Institute of Sound and Vibration Research,Hefei University of Technology,Hefei 230000,China;2.Institute of System Engineering,China Academy of Engineering Physics,Mianyang 621999,China)Abstract：The indirect boundary element method(IBEM)faces the problem of non-uniqueness in solving exterior acoustic problems,and therefore requires specific processing methods to ensure the stability of numerical results.When solving large-scale acoustic problems,the IBEM also faces the problem of limited computational scale,and the fast multipole method(FMM)is a commonly used method to overcome this problem.However,different non-unique solution processing methods have different impacts on the applicability and efficiency of the fast multipole method.To solve the aforementioned issues,the non-singular indirect boundary integral equations(IBIEs),based on constant elements discretization,are first derived in this paper,and then several non-unique solution problem processing methods are extended into the FMM,such as internal impedance,zero internal normal vibration velocity,and mixing potential.By contrasting numerical results acquired by various methods,the accuracy and computational efficiency of the fast multipole IBEMs are then determined.These conclusions can act as guides for the development of a more efficient fast multipole IBEM.Key words：exterior acoustic problems;indirect boundary element method;fast multipole method;non-uniqueness solution problem摘 要：间接边界元法在求解外域声学问题时面临着解的非唯一性问题，因此需要引入特定的处理方法以保证数值结果的稳定性。间接边界元法在求解大规模声学问题时还面临着计算规模受限问题，快速多极算法是目前克服此问题的常用方法，然而不同的非唯一解处理方法给快速多极算法的适用性和效率带来了不同影响。首先给出基于常值单元离散推导出的无奇异间接边界积分方程，然后将内部阻抗法、内部刚性面法以及混合势等几种常用的非唯一解处理方法引入快速多极间接边界元法，通过数值结果对比得出这些方法对快速多极间接边界元法求解精度和计算效率的影响规律，为发展更准确高效的快速间接边界元法提供了参考。关键词：外域声学问题；间接边界元法；快速多极算法；非唯一解问题中图分类号：O 242.2文献标志码：A 边界元法在声学分析中已得到了非常广泛的应用。按照边界积分方程中变量的定义方式，边界元法可以分为直接边界元法和间接边界元法。不同于直接边界元法采用实际声学量作为求解变量且要求模型边界封闭，间接边界元法的求解变量为声压差或振速差，因而能够直接获取模型边界两侧的声学信息；并且因间接边界元法可以将薄壁结构视为曲面划分网格，故在分析非封闭薄壁声学模型时，通常比直接边界元法具有更高的计算效率。文献 1-5验证了间接边界元法在薄壁结构声学分析中的高效DOI:10.13785/ki.nmggydxxbzrkxb.2023.02.006132内蒙古工业大学学报（自然科学版）2023 年 式中：2表示拉普拉斯算子；k=/c 是波数；是简谐声波的圆周频率；c 为声波在媒质中的传播速度；Q(x,y)是定义在 y 点的声源，x 表示场点；Q是声源源强；(x,y)为狄拉克函数；p(x)表示声压。本文略去时间依赖项 e-it（其中）。对于式(1)，运用格林第二等式并基于常值单元离散，即可得到无奇异的离散直接边界积分方程(16)。在此基础上，利用内外域边界积分方程的关系，可以推导出如下无奇异间接边界积分方程：式中：和 分别表示边界两侧声压差和振速差，即双层势和单层势的源；表示边界两侧的声压和，n 表示法向量方向；x、R()和 为常值单元的相关参数16；pin(x)表示入射波。2 解的非唯一性问题实际上在使用式(2)和(3)的间接边界积分方程进行外域声场预测时，在某些频率附近无法得到正确的数值结果，该现象常被称为解的非唯一性问题。为了克服这一问题，多种方法4,7-12,17-18相继被提出，下文将对其中代表性的方法进行介绍。2.1内部阻抗面法在直接边界元法中，Burton-Miller 法应用最广泛，该方法利用奇异边界积分方程和超奇异边界积分方程的组合式来克服求解外域声场问题时遇到的解的非唯一性问题。在间接边界元法中，文献 10通过在模型内表面施加阻抗，同样可以克服解的非唯一性问题，该方法等效于 Burton-Miller法10。该方法可以表示为如下内外域间接边界积分公式的组合9：性和准确性。虽然间接边界元法对很多复杂声学问题都能给出很好的数值结果，但该方法在求解外域声学问题时面临着数值稳定性问题。该问题具体表现为在某些频率附近无法得到正确的数值结果，这个现象常被称为解的非唯一性问题6-10。为提高间接边界元法在该类问题中的计算精度，学者们提出了一系列方法，如双表面法7-8、内部虚拟阻抗面法4,9-11。这些方法虽通过引入内部虚拟面能克服解的非唯一性问题，但内部虚拟面也会额外带来高计算量和高内存占有量的缺陷，限制了这些方法在大规模声学问题中的应用。近期，LEE J W12通过 SVD 分解间接边界元系数矩阵并结合基本解法构建了克服解的非唯一性的方法，该方法无需构建虚拟面，因而具有更好的计算效率，数值算例表明该方法能够在虚假特征频率处准确求解声学问题。除上述数值稳定性问题外，间接边界元法还存在求解效率问题。由于计算过程中所形成的系统矩阵为非对称满阵，由其引起大计算量制约了间接边界元法在大规模实际工程问题中的应用。快速多极算法(FMM)13-15是目前克服间接边界元法计算规模受限问题的常用方法。然而上述不同的非唯一解处理方法会给FMM 的加速过程带来不同的影响。例如对于设置内部阻抗面方法，由于需要同时在模型边界两侧指定不同的声学边界条件从而会造成计算规模增大，求解效率迅速下降，并且积分方程的具体形式还与所求解的声学问题密切相关；而若使用 SVD 分解方法则会因快速多极算法未形成显式的系数矩阵，造成该方法会受到求解酉矩阵的限制12。针对上述问题，本文首先给出基于常值单元离散推导出的无奇异间接边界积分方程，克服了间接边界积分方程的强奇异和超奇异边界积分问题；然后将内部阻抗、内部零法向振速以及混合势等几种常用的非唯一解处理方法引入快速多极间接边界元法中，提高了间接边界元法数值结果的稳定性和精确性；最后通过具有解析解的数值算例，对比得出这些方法对快速多极间接边界元法求解精度和效率的影响规律，为发展更准确高效的快速间接边界元法提供了一定的参考。1 声学间接边界元法对于本文所要研究的三维声学问题，可以用如下的 Helmholtz 方程来描述声波在均匀各向同性介质中的传播：(1)223()()(,)0,pk pQ+=Rxxx yx y(2)(3)133梁梦辉等消除声学间接边界元非唯一解的方法的对比研究第 2 期 式中：a,b,c 代表内外表面上的边界状态；p1是外域声压；p2是内域声压。基于上述方程，根据声学跳跃条件，可以得到：基于上述关系式，对于脉动球声辐射问题，可通过在内表面上施加阻抗吸收振动效应，即设置：a1=0,b1=0,c1=-i 和 a2=-ik,b2=1,c2=0，从而得到相应的间接边界积分方程：2.2内部刚性面法除了在内表面上指定阻抗来克服解的非唯一性问题以外，AMBARISHA V K 等9指出通过指定内表面法向振速为零同样也能够解决在求解外域声学问题时所遇到的解的非唯一问题。根据上述公式，对于脉动球声辐射问题，设置内表面 v2=0 和v10，据此得到 进而得到如下间接边界积分方程：2.3混合势法不同于上述在内外表面上指定阻抗或法向振速的方法，混合势法17-18通过间接边界元中混合势来联立间接边界积分方程，从而克服解的非唯一性问题。根据混合势法间接边界积分方程可以表示为 式中：表示混合势耦合系数，本文中=i/k，v表示混合势的源。2.4联合额外点基本解法除上述方法外，LEE J W 等12提出通过联合对间接边界积分方程与额外点的基本解方程亦可克服解的非唯一性问题。该方法类似于 CHIEF 方法，通过构建额外的超定方程来保证解的唯一性。相较于内部阻抗面和内部刚性面方法，该方法无需内表面上的边界积分方程，因而最终构建得到的方程不仅具有较高的计算效率，而且数值算法实现过程也更加简便。例如对于 Dirichlet 边界问题，联合单层势间接边界积分方程和额外点源的基本解方程，可以得到如下的组合式：式中：j表示额外点的强度，N 表示额外点源数量。为了求解上式需要定义额外的约束方程：式中：H表示与最小奇异值对应的右酉矩阵中的向量，H 表示 Hermitian 转置。对于式(6)和式(7)或式(9)，根据已知的振速边界条件，通过联立两个积分方程就可以求得未知间接边界元的声学变量，再利用式(14)就可以求出域内任意点处的声压值。(),=1,2,.,HjjN=0y(13)(12)(4)(5)(6)(7)(8)(9)(10)(11)(14)134内蒙古工业大学学报（自然科学版）2023 年对于采用混合势法得到的边界积分方程来说，根据已知边界条件，通过式(10)或式(11)就可以求得未知混合势的源，将其带入下述域内声压积分方程中即可完成对任意场点的求解。对于采用间接边界元法联合额外点基本解法得到的边界积分方程来说，根据已知边界条件，通过式(12)、式(13)就可以求得未知单层势的源和额外点的强度，将其带回式(12)中即可完成对任意场点的求解。3声学快速多极间接边界元法为保证间接边界元法在求解外域声学问题时的数值稳定性，上文介绍了几种能够有效克服解的非唯一性问题的方法，但这些方法在求解实际工程问题时仍然面临计算效率问题。为解决间接边界元法在进行大规模声学预测时所遇到的计算规模受限问题，本文将目前应用最广泛的快速多极算法引入上述间接边界积分方程中并对比得出这些方法对快速多极算法计算效率的影响规律。经过分析可知，Lee 所提方法并不适合采用快速多极算法加速，因为该方法需要对原系统方程的系数矩阵进行 SVD 分解，但在 FMM 中并未形成显式的矩阵。针对其余几种方法，本文通过引入快速多极算法构建了相应的快速多极间接边界元法。为了应用快速多极算法，应先将格林函数展开成如下的球谐波形式：其余各阶核函数可以表示为：式中：和的表达式为(18)式中：表示向量 的球坐标；jn和分别为 n阶球贝塞尔函数和第一类球汉克尔函数；为：式中：，表示伴随勒让德多项式。将上述展开式带入相应的间接边界积分方程即可得到相应的快速多极远场计算公式；在近场计算中，仍然采用式(6)(11)完成近场系数矩阵的求解过程；最后，联立近场系数矩阵和远场展开公式，采用 GMRES 算法迭代求解，即可完成声场的预测过程。4数值算例的对比与分析本文通过脉动球声辐射算例来对比分析几种方法对快速多极间接边界元法求解精度和计算效率的影响规律。以下分析中的声介质设为空气，其密度取为=1.2 kg/m3，声波在其中的传播速度为c=340 m/s。上述脉动球声辐射算例的解析解为：式中：p(r)表示任意测点的辐射声压；r 表示测点与球心的距离；vn表示脉动球球面上各质点的径向速度；ra表示脉动球半径。在本算例中，脉动球半径为 ra=1 m，测点位置为(6 m,0,0)，波数范围为 k=0.510，步长取为0.05。球面采用 6 860 个常三角形单元离散。在快速多极算法中，单个叶子节点内允许的最大单元数设置为 100，GMRES 收敛误差设置为 1E-5，最大迭代步数设置为 500，见图 1。(1)(,)()(,)(,)()(,)mmnnnmmnnnIk aj kr YOk ahkr Y =?,图 1 脉动球的示意图Fig.1 Diagram of a pulsating sphereOx1x3x2ra(15)(16)(17)(19)i(,)(cos)emmmmnnnYc P=(20)图 2 显示的是采用上述几种方法计算得到的声压，其中黑色实线表示采用式(20)求解得到的解析解(Ana)，蓝色圆圈线表示采用快速多极算法加速式(6)和(7)(简记为：Method 1)求解得到的数值结果，绿色方格线表示采用快速多极算法加速式(9)(简记为：Method 2)求解得到的数值解，蓝色135梁梦辉等消除声学间接边界元非唯一解的方法的对比研究第 2 期菱形线表示采用快速多极算法加速式(10)和式(11)(简记为：Method 3)求解得到的数值解。从图 2 可以看出，在整个波数范围内采用Method 1、Method 2 和 Method 3 求解得到的数值结果与解析解均能很好地吻合。因此可以得出，采用这些方法能够克服求解外域声学问题时可能遇到的解的非唯一性问题。另外，从图 3 所示的误差曲线还可以发现：3 种方法数值结果的相对误差保持在同一精度水平，而 Method 1 数值结果的相对误差相较于其余两种方法在 k=610 范围内出现略微波动；Method 2 数值结果的相对误差基本保持在相对较低的误差水平，这与文献 4 中的结论相吻合。除上述对单个测点的扫频分析外，本文还在虚假特征值 k=4.514 10 和 7.725 25 处求解相应的声压，并绘制出声压分布云图 4 和图 5。观察该图可以发现：采用快速多极算法加速式(3)(简记为：FM-IBEM)在脉动球球面附近无法得到正确的数值结果，而采用 Method 1、Method 2 和 Method 3 能够准确求解该声学问题，更为重要的是采用这些方法能够避免对系数矩阵进行 SVD 分解来克服解的非唯一性问题，提高了数值求解效率。虽然采用 Method 1、Method 2 和 Method 3 均能准确求得数值结果，但这些方法在求解效率上略有差异。正如表 1 所示，Method 1 由于需要联立式(4)和式(5)，导致矩阵维度从原来 NN 的变成了 2 N2 N，致使计算时间大于采用 Method 2 和 Method 3 所需时间；并且对于不同的声学问图 2 脉动球声辐射声压图Fig.2 Comparison of the sound pressure of a pulsating sphere(a)声压实部(b)声压虚部图 3 3 种方法的相对误差Fig.3 The relative error of the three methods图 4 k=4.514 10 时声压分布云图Fig.4 Sound pressure distribution at k=4.514 10(d)Method 3(a)FM-IBEM(b)Method 1(c)Method 20246810Wave number-100-50050100p(x)(real part)AnaMethod 1Method 2Method 3k/m-10246810Wave number k/m-1-100-50050100p(x)(imag part)AnaMethod 1Method 2Method 30246810Wave number 246810121416Relative error10-3Method 2Method 3k/m-1图 5 k=7.725 25 时声压分布云图Fig.5 Sound pressure distribution at k=7.725 25(a)FM-IBEM(b)Method 1(c)Method 2(d)Method 3136内蒙古工业大学学报（自然科学版）2023 年题，Method 1 的具体形式还会改变；而 Method 2 和 Method 3 仅需联合间接边界积分方程即可克服解的非唯一性问题，因而计算效率更高。但是Method 2 仅在脉动球声学问题中表现较好，在这种情况下 Method 3 是更为简便的，而且更容易扩展到 FMM，从而提高间接边界元法的求解效率。5结论针对采用间接边界元法求解外域声学问题时可能遇到的解的非唯一性问题，本文引入了几种常用的处理方法，进一步针对间接边界元法在求解大规模声学问题时计算规模受限的缺陷，本文将上述几种处理方法拓展到快速多极算法中。通过具有解析解的数值算例对比讨论了上述几种处理方法对快速多极算法求解效率的影响规律，最后得出上述几种处理方法均能有效克服求解外域声学问题时可能遇到的解的非唯一性问题，而基于混合势方法构建的快速多级间接边界元法具有较好的计算效率。参考文献1 PISCOYA R,OCHMANN M.Acoustical greens function and boundary element techniques for 3D half-space problemsJ.Journal of Computational Acoustics,2017,25(4):1730001.2 MATHEWS I C,JEANS R A.An acoustic boundary integral formulation for open shells allowing different impedance conditions,top and bottom surfacesJ.Journal of Sound and Vibration,2007,300(3/5):580-588.3 COOX L,ATAK O,VANDEPITTE D,et al.An isogeometric indirect boundary element method for solving acoustic problems in open-boundary domainsJ.Computer Methods in Applied Mechanics and Engineering,2017,316:186-208.4 WU Y H,DONG C Y,YANG H S.A 3D isogeometric FE-IBE coupling method for acoustic-structural interaction problems with complex coupling modelsJ.Ocean Engineering,2020,218:108183.5 GUMEROV N,ADELMAN R,DURAISWAMI R.Fast multipole accelerated indirect boundary elements for the Helmholtz equationJ.Proceedings of Meetings on Acoustics,2013,19(1):015097.6 ZHENG C J,CHEN H B,GAO H F,et al.Is the BurtonMiller formulation really free of fictitious eigenfrequencies?J.Engineering Analysis With Boundary Elements,2015,59:43-51.7 MOHSEN A,HESHAM M.An efficient method for solving the nonuniqueness problem in acoustic scatteringJ.Communications in Numerical Methods in Engineering,2006,22(11):1067-1076.8 MOHSEN A,PISCOYA R,OCHMANN M.The application of the dual surface method to treat the nonuniqueness in solving acoustic exterior problemsJ.Acta Acustica United With Acustica,2011,97(4):699-707.9 AMBARISHA V K,GUNDA R,VIJAYAKAR S M.A new indirect formulation to address the non-uniqueness problem in acoustic BEMC/INTER-NOISE and NOISE-CON Congress and Conference Proceedings,NoiseCon07,Reno NV:Institute of Noise Control Engineering,2007:1046-1055.10 DAMICO R,NEHER J,WENDER B,et al.On the improvement of the solution accuracy for exterior acoustic problems with BEM and FMBEMJ.Engineering Analysis With Boundary Elements,2012,36(7):1104-1115.11 邹明松,吴有生,祁立波,等.声介质中水弹性力学消除不规则频率的方法 J.船舶力学,2013,17(10):1202-1208.12 LEE J W,CHEN J T,NIEN C F.Indirect boundary element method combining extra fundamental solutions for solving exterior acoustic problems with fictitious frequenciesJ.The Journal of the Acoustical Society of America,2019,145(5):3116-3132.13 GREENGARD L,ROKHLIN V.A fast algorithm for particle simulationsJ.Journal of Computational Physics,1987,73(2):325-348.14 ZHENG C J,MATSUMOTO T,TAKAHASHI T,et al.A wideband fast multipole boundary element method for three dimensional acoustic shape sensitivity analysis based on direct differentiation methodJ.Engineering Analysis With Boundary Elements,2012,36(3):361-371.15 ZHENG C J,CHEN H B,CHEN L L.A wideband fast multipole boundary element method for half-space/plane-symmetric acoustic wave problemsJ.Acta Mechanica Sinica,2013,29(2):219-232.16 MATSUMOTO T,ZHENG C J,HARADA S,et al.Explicit evaluation of hypersingular boundary integral equation for 3-D Helmholtz equation discretized with constant triangular elementJ.Journal of Computational Science and Technology,2010,4(3):194-206.17 KIRKUP S.The boundary element method in acoustics:a surveyJ.Applied Sciences,2019,9(8):1642.18 KIRKUP S M.The boundary element method in acousticsM/OL.Hebden Bridge:Integrated Sound Software,2007.2019-05-10.https:/clok.uclan.ac.uk/7333/1/tbemia07.pdf.19 BURTON A J,MILLER G F.The application of integral equation methods to the numerical solution of some exterior boundary-value problemsJ.Proceedings of the Royal Society of London.A.Mathematical and Physical Sciences,1971,323(1553):201-210.表 1 Method 1、Method 2 和 Method 3 计算时间的比较Table 1 Comparison of CPU time of Method 1,Method 2 and Method 3单元数Method 1计算时间/s Method 2Method 3608059.247 430.844 424.658 9626266.752 738.732 527.955 7818088.713 648.164 736.192 7975692.983 556.810 840.546 9