声学边界元方法及其快速算法

第1章 绪论 1.1 引言 对于车辆、飞行器、船舶及家用电器等机电产品，声学性能不仅是强制性要求，往往也是与对手产品竞争中取胜的关键。例如，民用客机噪声不仅是决定乘客舒适性的重要技术指标，而且可能导致机身结构产生疲劳损伤，影响飞行安全。家用电器的噪声水平则更是影响顾客选购的重要因素。在这些产品的开发过程中，噪声定量预报、噪声控制与优化是产品声学设计的重要任务。 实际工程对象比较复杂，如图1.1所示的一些大型复杂设备声学分析问题，包括噪声的定量分析计算、预测和优化，无法使用解析的方法，而实验分析只能在样品或成品阶段才能进行，且成本高、周期长。因此，数值计算是工业产品声学分析的重要手段，声学数值方法主要包括有限元方法（finite element method， FEM）、边界元方法（boundary element method，BEM）、有限差分法、声线追踪法、统计能量分析法等。声学问题中大量声学物理量可以直接通过边界元方法进行求解，或间接通过边界元方法进行分析。因此，边界元方法作为*有效的数值计算方法之一，在声学数值分析方法中占有重要的地位。 图1.1 大规模声学问题 以辐射声场及其特性为主要研究目标的声学计算，可以分为自由空间、多联通域和半空间问题等。潜艇螺旋桨噪声、飞机发动机辐射噪声等都可视为自由空间的声学问题。带声学覆盖层水下航行器的声辐射和声散射、多孔吸声材料的设计及生物组织中声学特性模拟等，属于多联通域声学问题。水下航行器靠近海平面或海底面时的声辐射和声散射分析，各类车辆的声辐射、飞机起飞或着陆过程中的噪声分析等，都要考虑边界的影响，需要建立一般阻抗平面的半空间声学模型。 边界元方法是一种基于积分方程的数值方法，可以将模型降维（三维变为二维，二维变为一维），且仅需要离散边界。声学Helmholtz边界积分方程能自动满足远场辐射条件，比较容易处理无限域的边界条件。同时作为一种半解析的数值方法，边界元方法较之其他数值方法精度要高，是一种比较适用于声学分析计算的数值方法。 虽然边界元方法降低了模型分析的维数，但是其系数矩阵是非对称，甚至是奇异的满阵，其数据存储量与模型离散自由度（N）的平方成正比，即O（N2）。在双精度数值运算下，边界元方法内存使用量和离散自由度的关系如图 1.2 所示，可以看出，当离散自由度为2×104 时，边界元方法需要使用近6GB 的内存，当离散自由度为3×104 时，其内存使用量增长到惊人的13.41GB 。根据目前个人计算机的一般性能来看，传统边界元方法很难完成超过20000 个自由度模型的声学计算。对于更大的模型，边界元方法的内存使用量将超出计算机的存储范围而无法进行计算。另外，当边界元方法中使用迭代求解器计算时，其求解复杂度为O（N2）。如果使用直接求解器（如Gauss消去法），其求解复杂度为O（N3）。因此，边界元方法的内存使用量和计算耗时严重地制约着其在声学分析问题中的应用。 图1.2 不同离散自由度下边界元方法的内存使用量 随着人们对产品声学性能的不断追求，声学分析和研究的问题逐渐向大规模、多物理场方向发展，迫切需要边界元方法完成大规模声学问题的快速计算分析。快速多极子边界元方法的出现和发展，极大地降低了边界元方法的内存使用量，提高了边界元方法的计算效率，为大规模声学问题的计算和分析提供了一种解决方法。 1.2 边界元方法 1.2.1 概述 边界元方法可分为直接法和间接法两种类型。直接法以物理意义明确的结构表面声压和法向振速为未知量，求解Helmholtz边界积分方程，适用于表面封闭结构的声辐射和声散射问题的计算。间接法以物理意义不明确的变量，如结构表面的声压差和速度差，求解Helmholtz边界积分方程，适用于表面不封闭结构的声辐射和声散射问题的计算。本节简要介绍直接边界元方法（简称边界元方法）的发展历史和特点。 边界元方法是边界积分方程（boundary integral equation， BIE）的一种数值解法，它的基本思想是基于Green函数，采用Gauss定理，把一个封闭区域上的积分转化为该区域边界上的积分。边界元方法的发展可以追溯到20 世纪60年代，Jaswon等*初用它来求解二维位势问题[1-3]。随后，Rizzo在其博士论文中将边界元方法用于二维弹性力学分析，并于1967年将成果发表在期刊上[4]。20 世纪60～70年代，众多学者在应用力学领域展开了大量关于边界积分方程及其求解方法的研究[5-13]。仿照有限元方法的命名，边界元方法术语在1977年由Brebbia 、Dominguez 、Banerjee和Butterfield在英国南安普顿大学共同商定，并在Brebbia编写的The Boundary Element Method for Engineers[14]一书中正式使用。 美国、欧洲和中国的一些学者，如Rizzo[15]、Cruse[16]、Watson[17]、Shippy[18]、Mukherjee[19]、Telles[20]、姚振汉和杜庆华[21，22]等，就各自领域边界积分方程及边界元方法的发展历史在电子期刊Electronic Journal of Boundary Elements的特刊中做了综述性的描述。2005年，A.H.D. Cheng和D.T. Cheng发表了一篇综述文献，详细描述了边界积分方程和边界元方法的早期历史和发展、18 世纪至20 世纪上半叶边界积分方程发展的重要数学贡献，以及20 世纪70年代以前边界元方法发展的历史[23]。Tanaka在1983年综述了当时边界元方法的*新发展[24]，随后，他又总结了20 世纪80～90年代奇异及超奇异边界积分方程的发展。经过几十年的发展，多本专著陆续被出版[25-27]，这可以作为不同领域边界元方法理论系统学习的资料。 在我国，边界元方法的研究工作起步较晚。1978年，清华大学杜庆华教授在国内首先开展边界元方法的研究，推动了边界元方法在力学、电磁学等学科上的迅速发展。边界元方法是继有限元方法之后的一种别具特色的新型数值方法。与有限元方法相比，边界元方法有如下三个优点： （1）降维。边界元方法可以将所分析模型的维数降低一维，而且不需像有限元方法一样离散整个求解域，仅需要离散空间模型的边界，即对于二维和三维问题，可分别采用线单元和面单元离散边界，因此减少了离散单元数，并大大降低了模型网格划分的难度。 （2）求解精度高。边界元方法是一种半解析数值方法，具有解析与离散相结合的特点，求解精度比一般的数值方法高。其误差主要来源于边界单元的离散，累积误差小，便于控制。 （3）适用于无限域问题。Helmholtz边界积分方程能自动满足无穷远处的辐射条件，无须特别处理外部问题在无限远处的边界条件，便于无限域声场分析。 但是，理论分析和工程实践应用表明，边界元方法也存在自身固有的缺陷，主要表现在以下三个方面： （1）内存使用量大和计算效率低。边界元方法所生成的线性方程组的系数矩阵是稠密、非对称满阵，其内存使用量为O（N2），N表示离散模型的自由度。当使用直接法求解时，如Gauss消去法，其计算量为O（N3），即使采用迭代法求解，其计算量也为O（N2）。虽然边界元方法降低了求解问题的维数，但仅适用于小规模问题的分析，这是阻碍边界元方法发展的主要原因，也是边界元方法的发展落后于有限元方法等其他数值方法的主要原因。 （2）非唯一性。对于外部问题，当分析频率与内部结构的共振频率重合时，解不能保证唯一性。这种不唯一现象不是现实物理问题中真实存在的，而是由采用边界积分方程的数学方法造成的。因此，需要添加辅助方程来消除非唯一性的影响，从而增加了计算量。 （3）奇异积分。当场点和积分点重合时，边界积分方程中会出现的奇异性，其计算精度对*终计算结果的影响较大。当引入边界积分的导数方程来消除外部问题的非唯一性时，方程中会出现的超奇异性，进一步增加计算的困难。 1.2.2 声学边界元方法的发展 虽然边界元方法存在上述不足，但是由于其在无限域外部问题处理上的优越性，在声学分析领域仍受青睐。经过几十年的研究，声学边界元方法也在一定程度上得到了发展，并成功应用于一些实际声学问题的分析中。 1963年，Banaugh和Goldsmith使用边界积分方程，分析了任意形状二维物体的稳态声散射问题[28]。受当时计算条件的限制，研究者仅求解了36 维的离散方程，这是边界积分方程计算声学问题的一次有意义的尝试。同年，Chen和Schweikert使用Fredholm积分方程预测了任意形状三维物体的声辐射问题[29]，显示了边界积分方程在计算任意形状模型声辐射问题的优越性。Chertock使用边界积分方程研究了振动结构的辐射声场，并指出基于边界积分方程的通用数值方法随着离散自由度的提高，计算效率变得低下[30]。Copley研究了一种简单的方法用于辐射问题的边界积分方程的计算，并提出了对称模型的非奇异边界积分方程[31]。随后，Copley发现边界积分方程在共振频率处求解会产生非唯一性问题，他是公开文献里**个发现此现象的学者[32]。 1968年，Schenck提出了一种改进的边界积分方程法，即CHIEF（combined Helmholtz integral equation formulation）法，用于解决边界积分方程的非唯一性问题[33]。CHIEF 法通过在结构内部增加Helmholtz积分点，形成超定方程系统，并通过*小二乘方法求解系统的未知量。对于简单的声学模型，CHIEF 法非常适用，只要选择合适的配置点就能有效解决非唯一性问题。然而，对于结构复杂的模型和高频声学问题，CHIEF 法配置点的数目和位置很难确定。因此，CHIEF 法在工程应用中具有一定的局限性。1971年，Burton和Miller利用边界积分方程及其方向导数方程各有一组互不重叠共振频率的特性，将两个方程进行适当的组合，消除了求解结果的非唯一性问题[34]，因为这在边界元方法领域产生重要影响，后人将此方法以两人姓名共同命名为Burton-Miller法。 虽然Burton-Miller法解决了共振频率处解的非唯一性问题，但是在方程中引入了超奇异积分项，奇异积分是边界元方法的研究热点和难点[35-37]。对于二维声学问题，存在严格的解析奇异积分方法[38]。对于三维声学问题，大量研究直接从Gauss数值积分出发来解决奇异积分的问题[39，40]。对于奇异积分，采用较多的积分点可以保证计算的精确性。但对于超奇异积分问题，即使采用更多的积分点，也不能保证计算的精度。总的来说，奇异积分技术大致可以分为正则变换法、坐标变换法和奇异值消去法三类[41-45]。Liu等在Burton-Miller方程中考虑了立体角的导数，使用对角元素重生法处理了改进的Burton-Miller方程中的奇异积分[46]。Chen等基于奇异值消去法重新推导了Burton-Miller方程的表达式[47]，使其对任意高阶单元均具有弱奇异的特性，随后Li和Huang将其用于快速多极子边界元算法中[48]。Polimeridis等将直接积分法和奇异值消去法相结合，处理了Galerkin边界元方法中的奇异积分[35]。但对于常数单元，上述计算方法仍然比较烦琐。Matsumoto等对常数单元离散的Burton-Miller方程提出了显式的无奇异积分表达式[49]，并将其应用到声学灵敏度的分析中[50]。Tadeu和António也推导了常数单元的奇异及超奇异积分的表达式[37]。Wu等将配点置于平面单元的内部，建立了线性连续单元的非奇异Burton-Miller方程表达式[51]。 Cheng等用边界元