水下大型结构振动声辐射预报技术引论

如前文所述，利用弹性结构表面的振动预报声场是声辐射预报研究中的经典问题，近年来学者对该问题的研究主要集中在吸收边界条件法、无限元法和边界元法等。

（一）人工边界法（ABC）

大量的数值实验结果表明，利用直接在虚拟边界上施加Sommerfeld辐射条件的方法计算声场会导致严重的虚反射，而声能顺利的透过虚边界这一自然现象又很难用简单的数学表达式进行阐述，寻找简单而有效的数学表达，是此类研究的焦点。1977年，Engquist和Majada[15]率先提出了基于波动方程旁轴近似的吸收边界，通过令所有的波都在该边界上被吸收，使得截断区域上与原有的外场问题精确等价。以此为代表的低阶吸收边界成了20世纪80年代研究的主流，在众多研究中较有代表性的是廖氏边界[16，17]和Cerjan等人[18]提出的利用有效吸收层衰减外行波的海绵吸收层。低阶吸收边界的共同特点是数学表达式较为简单，易于实现，然而在对其进行数值离散时都会遇到困难，采用低阶离散则难以保证计算精度，而要引入高阶插值函数又需要处理高阶偏微分，如何在实现高阶吸收边界的同时避免高阶偏微分，仍有待进一步的研究。

上世纪九十年代中期由Berenger[19]在对电磁场计算的研究中提出的完全匹配层（PML）是吸收边界方法发展中重要的里程碑，通过沿坐标轴人为地将场变量分解为多个分量，并在各方向上分别引入衰减因子以确保人工边界能够吸收任意角度入射的波，利用多个吸收层迭代的吸收波，该方法同时具备了简洁的数学形式和无反射的特性，对倏逝波也有较好的吸收效果。由于引入了对标量场变量的分裂，因此该方法也被称为分裂式PML（SPML）。1996年，Zhao和Cangellaris[20]提出，对场变量的分裂并非必要的，并推导了与Berenger形式相等价的非分裂PML（NPML），显而易见的是，NPML降低了人工边界上未知量的数目，也就降低了问题的计算规模，然而在对其进行数值离散的过程中，考虑到近场条件下倏逝波的存在，且入射角变化较大，如果采用单一的衰减因子，就需要使人工边界与辐射源保持一定的距离以保证解的精度，反而增大了整体的计算规模。1996年，Chew等人[21]引入复数伸展坐标，通过为PML域中的笛卡尔坐标添加虚部来替代介质中的衰减因子，给出了更具一般性的PML表达。在此后的十几年中，PML方法引起了学者的广泛重视，被应用于众多领域。在Chew等人工作的基础上，Liu和Tao[22]率先提出了Helmholtz方程PML的分裂形式，该方法保留了PML简洁的数学形式，然而在实际应用中必须对边界存在的棱边、角点进行单独的处理，不利于编程实现，也限制了该方法在工程应用中的推广。1997年，在对地震波的研究中，Berenger[23]提出了基于复频移（complex frequency shift）的PML改进形式（CFS-PML），引入复波数替代阻尼因子，实质上是将复数伸展坐标系中的实轴加以放大，并将虚轴零点平移到了负半轴，改善了NPML对倏逝波的吸收能力。Roden和Gedney[24]在对CFS-PML加以实现的过程中引入了递归卷积，在保持CFS-PML的非分裂特性同时，能够更有效的截断计算各向异性介质、色散介质、有耗介质问题。

在最近几年，PML方法也正受到国内学者的重视。方能胜[25]在其博士论文中对PML的不同形式进行了深入的总结，指出虽然这些表达式从数学角度看是存在等价性的，然而在进行了数值离散后，它们的数值稳定性不尽相同，必须分别加以分析。李佩笑等人[26，27]在总结常规的NPML和SPML基础上，对二者性能以及在声学计算中的适应性进行了研究，指出相比于SPML，NPML占用内存更少，计算速度更快，在对三维问题的计算中这一优势更为明显，然而NPML对斜入射波和倏逝波具有更好的吸收能力，可以使吸收层更接近结构表面，同等吸收精度条件下，NPML方法对吸收层的厚底也更低，因此在工程应用中仍需根据实际情况选择合适的方法。

总体看来，PML本身显然不具备明确的物理意义，但它从数值角度确保了：

1）截断边界对任意角度入射的平面波的反射系数为零；

2）波在PML中传播时呈指数率衰减。

由于以上两点特性，声波穿透PML后变得非常微弱，再由PML外沿向内反射的波已经可以忽略不计。然而在实际应用中，数值实验表明PML方法的准确性受到匹配层厚度、虚部因子的分布函数等因素的影响；收敛性对频率的变化较敏感，尤其在低频段需要设置更多的匹配层，使得计算量大幅提高。

（二）无限元法

无限元法的历史可以追溯到1977年，Bettess、Zienkiewicz等人[28]在对海面上水波的衍射和折射问题的研究中首次提出了infinite element这一名词，其基本思想是沿着波的传播方向建立无限长的带状单元，并构建指数衰减型的形函数以满足Sommerfeld辐射条件。1984年，Bettess等人[29]提出了映射无限元，也被称为Bettess单元，利用整体坐标和局部坐标之间的奇异映射，将整体坐标系下的无穷远点影射到局部坐标系下的奇点，避免了原有方法对人工设置辐射源点的需求。1998年，Bettess等人[30]进一步对原有方法进行完善，允许几何影射在单元内、单元间变化，从而使得无限元对较大纵横比求解域中问题的求解精度大为改善。(www.xing528.com)

Astley和Eversman[31]提出了波包络线法，利用形函数的共轭作为Galerkin法中的加权函数，消除了Bettess映射单元中的不确定积分，使单元矩阵可直接由Gauss积分得出。1994年，Astley和Macaulay[32]给出了Bettess单元的标准几何影射，使得任何实际的无限单元均可由局部坐标系下的标准单元影射得出，进一步将该方法发展为映射共轭无限元，也被称为Astley-Leis单元。Cremers和Fyfe[33，34]将Astley-Leis单元推广到了任意阶次，并首次明确地将几何影射与形函数分别进行研究。总体而言，Astley-Leis单元是对Bettess单元的发展，在避免了不确定积分的同时，大大减少了无限单元解决实际问题中所需的人工干预，可以方便地实现无限方向上的二阶、三阶甚至更高阶插值，代价则是刚度矩阵失去了对称性，这是由于其取形函数为加权函数的共轭导致的。此外，由于波包络线法形状函数的构造理论是基于几何、物理意义上均不存在的虚拟声源概念，缺乏严格的理论推导，其形函数带有明显的试凑性质。

Burnett等人[35-37]针对无限元求解狭长结构声辐射、散射问题精度和效率上的不足，提出了一系列新的无限单元。1994年，Burnett首先采用共焦椭圆变换研究扁长形椭球结构的声辐射问题。通过引入多极展开理论，将有限的椭球面用常规的有限元法处理，无限方向则用一维无限元法模拟。之后在1998年，Burnett又将该理论推广到了扁圆形椭球和一般椭球的声辐射问题中。相比于球坐标系，Burnett元的创新在于：（1）基于级数展开原理构建径向形函数，保证了无限元的高阶近似，拓展了形函数构建理论的研究空间；（2）实现了笛卡尔坐标系向曲面坐标系的变换，放松了高阶单元中对节点位置的限制，使得结构网格与无限元网格之间所需填充的流体有限元网格数量大幅度减少，显著提升了计算效率。（3）通过将形函数分离为半径和角度两个独立的部分，且只有径向积分与波数相关，降低了系统矩阵的计算成本。

国内对于声学无限元法的理论研究起步较晚，但近年来正逐渐引起学者的重视。2004年，杨瑞梁、汪鸿振[38]对声无限元法的发展进行了系统性的回顾，并在其博士论文[39]中对Burnett元计算误差受条件数、角向量等因素的影响。2014年以来，苏楠[40-42]等人对Astley元展开了一系列研究，并将其应用于水下加筋圆柱壳的声振特性研究中。与此同时，对声学无限元在声辐射问题中的应用研究也呈逐年增多的趋势。

（三）边界元法

边界元法又称为边界积分方程法（Boundary Integral Equation Method），它以定义在边界上的边界积分方程为控制方程，通过对边界上离散单元的插值计算，将边界积分方程化为线性代数方程组进行求解。目前应用较为广泛的边界元法大致有直接边界元法（Direct BEM）、间接边界元法（Indirect BEM）、双重边界元法（Dual BEM）、子结构边界元法（Substructure BEM） 以及快速多极子边界元法（Fast Multipole BEM ）等。

与声无限元法相比，直接边界元法对自由场中结构振动-声辐射预报问题具有更好的适应性，主要表现在它不需要对结构外部有限距离内的流体进行离散，而是直接以结构表面为求解域，建立振动位移、振速与表面声压之间的联系，在求解得到结构表面声压、振速后又可以容易地利用Helmholtz积分将解外推至声场空间中任意一点。早在1963年，Schweikert和Chen[43]成功运用简单源法，即在假设封闭结构表面分布了一层单极子源，设它们的强度为未知量，解决了一系列声辐射问题。在之后的数年中，多位学者[44-46]从虚假本征频率、奇异积分等方面对边界元法进行了完善，形成了经典的声学边界元方程。其中具有代表性的是Burton和Miller[47]提出的，考虑Helmholtz积分方程在表面上各个参考点法向求导后的方程，与原积分方程进行线性叠加，从而虚假本征频率出解不唯一的问题，得到全频率适用的边界元方程，然而对边界积分方程的求导引入了超奇异积分，为数值计算带来了困难[48，49]。

虽然边界元法可以将问题的维度降低一维，然而由于其系统矩阵为密集矩阵，随着问题规模的增加，其计算耗时和内存占用都迅速增长，因此如何降低边界元法的计算复杂度成了过去二十年中学者研究的焦点，成果主要可分为两大类：以提高矩阵-向量乘积运算速度为目标的算法[50]和矩阵压缩算法[51]。前者以快速多极子边界元法为代表，后者更接近于纯代数方法，以层级矩阵技术和自适应交叉近似算法为代表。

快速多极子边界元法的思想来自于Rokhlin[52]在对二维Laplace方程快速求解算法研究中提出的快速多极子算法，通过将核函数分成两部分，一部分仅与源点位置和配置点位置相关，另一部分仅与配置点位置及场点位置相关，由于这两部分可以分别重复利用，放松了物理场中各部分的耦合，从而使得方程的迭代求解器单次矩阵与向量相乘操作的运算量从n2级别降到了n级别。随后Greengard[53]进一步完善了Rokhlin的算法，采用树结构对空间域进行逐级剖分，引入多极展开（Multipole expansion）和局部展开（Local expansion）的概念，使多极子算法正式成型。1990年Rokhlin[54]首次将多极子算法推广到Helmholtz方程的求解中。1993年Coifman等人[55]首次完整阐述了多极子算法在波动方程中的实现过程，对算法中涉及的各个参数的选择进行了早期的探讨。同年，Rokhlin[56]提出了三维声学中的传递因子对角化理论，为Helmholtz方程的核函数展开远场传递和计算提供了更为快捷、有效的方法，从而改善了FMM算法对中、小规模声场问题的求解性能。随后在 1995年，Epton等人[57]对应用 FMM 计算Laplace方程和Helmholtz 方程时的对角化转移因子理论作了更进一步的完整的描述。1997年Fukui等人[58，59]首次将多极子算法与二维声辐射问题的边界元法结合，形成了应用于声学的快速多极子边界元法（FMBEM）。之后Sakuma和Yasuda[60]将FMBEM应用于声场预报问题，从计算量、内存消耗等方面比较分析了FMBEM与传统边界元法的性能差异。由于树状多极子算法需要在初始化过程中将网格分割成不同的枝叶，经常有单元同时处于两个不同的分枝中，而常数单元的插值点在单元的中心，为树的构建带来了便利因此学者多主张采用常数单元对声场边界进行离散。CJ Zheng等人[61]导出了常单元下边界积分方程中各项奇异积分的解析形式。在此基础上，上海交通大学的蒋伟康和吴海军[62]对基于对角传递理论的FMBEM进行了研究，导出了底层多极展开系数的解析形式，从而进一步提高了算法的效率。然而Marburg[63]在对波长、网格密度与解的收敛性的研究中指出，要获得同样的收敛性，常单元需要的网格密度要比线性单元高60%。2004年，Fischer[64]构建了基于线性单元的声学FMBEM，数值实验显示该方法对网格密度的要求明显低于常单元FMBEM，但其采用的双重面积分方法十分繁琐，使运算量大幅增加，此外由于缺少有效的预条件算法，其迭代求解过程的收敛速度也较慢。

1998年，Hackbusch[65]提出了层级矩阵（hierarchical matrix，简称H-矩阵）的思想，其研究的目标是利用较少的数值来描述密集矩阵，从而达到压缩矩阵存储空间的目的，而压缩后的矩阵形式又不同于一般的稀疏矩阵，因此被称为数据稀疏型矩阵[66]（data-sparse matrix）。层级矩阵的核心思想是对密集矩阵进行层层分割，将子块的位置、大小等信息按树状的数据结构加以保存，对于其中秩较低的子块进行奇异值分解（SVD），保留其余分块的密集格式。在之后的数年中，该方法被应用于声场、电磁场问题的计算中。在层级矩阵的基础上，Bebendorf[67]于2004年提出了自适应交叉近似（ACA）算法，该算法主要从两个方面完善并改进了层级矩阵法，首先是利用交叉近似方法对矩阵中的低秩子块进行逼近，替代了SVD算法，从而降低了对矩阵进行压缩的计算复杂度，其次是通过考察积分核的退化特性使算法具备了自适应的能力，不再需要对矩阵子块的秩进行事先估计，甚至不要求矩阵元素为已知，而是根据迭代逼近过程中的中间结果自适应的找出必要的矩阵元素位置，再由积分算法给出相应的值，使得算法的复杂度真正有O（n2）级别降低到了O（n*logn）。Grasedyck[68]在交叉逼近的过程中引入了参考行和参考列以改进ACA算法的逼近策略，避免了原有算法在少数情况下崩溃的风险，也使得ACA算法的收敛速度有了轻微的提高，同时首次将ACA算法应用于边界元的快速计算[69]。在之后的十年中，ACA算法被应用于弹性力学[70]、电磁学[71]问题的边界元快速计算中，取得了良好的效果。从积分方程的角度看，H-矩阵与FMM算法都是基于这样一个思想，即对于一个积分方程中的积分算子L，在积分区域中具有退化性，则可以用一系列的分解对积分进行近似，从而降低积分的运算量，达到快速求解积分方程的目的。由于FMM直接利用了积分核函数多极展开的方式，因此需要对不同的问题、不同的插值阶次分别进行研究，而基于H-矩阵理论的方法着眼于矩阵本身的特性，对物理问题的类型、插值阶次相对不敏感，因而具有更好的通用性。

除了上述两种方法，双重边界元法[72]和子结构边界元法[73，74]也是较有代表性的快速边界元方法。这两种方法均出自于20世纪80年代学者对舱室、消声腔体等拓扑形状复杂的内部声场问题的研究。其中双重边界元法最早作为一种将Burton-Miller方法推广到内场问题中、处理虚假本征频率的方法被提出来，通过在厚度可忽略的退化边界或虚拟边界上考虑振速连续条件，并在退化边界、虚拟边界处以两侧声压差作为系统矩阵的未知量，在缩减自由度的同时简化了分析过程。子结构法是根据模型结构或声传播介质分布特点，将声场内部区域划分为若干个独立的有完整边界的子区域，离散各个子区域并应用边界元法建立系数矩阵，然后利用交界面上声压与质点振速的连续性条件，把各子区域的矩阵方程联系起来并重新组合，构成整体系数矩阵方程组，最后计算求解即可得到包括交界面在内所有边界上各节点处未知量。1987年Terao和Sekine[75]首次将子结构边界元技术应用于分析空气管道的声学性质。1991年Cheng 和 Seybert[76]将子结构边界元法应用到了具有内插管结构的消声器声场计算中。崔晓兵[77]指出，由于子结构法降低了求解域各个子域耦合程度，使得系数矩阵成带状稀疏结构，尤其对于交界面节点较少的细长结构，子结构法对边界元的数值性能有十分显著改善。