子结构快速多极子边界元法声学迭代计算收敛特性分析

崔晓兵

（中国舰船研究院 北京100192）

引 言

将快速多极子算法应用到传统边界元法中，研究形成快速多极子边界元法（FMBEM），凭借其计算效率高、存储空间要求小和传统边界元法具有同等计算精度等特点，被广泛应用于大尺度声场的三维计算中［1］，但对于复杂结构声场计算将存在不足，例如对于复杂结构的排气消声器的内部声场计算，有以下问题需要解决：

一是对于抗性消声器，其往往具有进出口内插管、穿孔管、隔板等薄壁结构，在离散求解时出现奇异性边界，FMBEM无法直接应用；二是对于阻抗复合型消声器，由于结构中存在吸声阻尼材料，由于声场波数的不同，声学域需要按照不同介质进行划分，各区域间的传播与耦合问题将使一般的FMBEM应用受到限制。

为了解决上述问题，将子结构技术［2］应用到快速多极子边界元法中，发展形成子结构快速多极子边界元法，应用于含有薄壁结构和带有阻性材料的声场快速计算中，由于FMBEM直接计算各行的矩阵向量积，并不形成完整的系数矩阵，常用Krylov子空间迭代法来求解线性方程组。所以，子结构快速多极子边界元法的迭代收敛速度无疑是其在复杂结构声场中成功应用的关键。本文将以消声器内部声场计算为例，分析影响子结构FMBEM的迭代次数的主要因素，对矩阵方程的建立进行优化，找到使其快速收敛的方法，通过与传统边界元法比较来验证其计算精度。为子结构FMBEM在大尺度、高网格数声学问题中的应用提供指导意见。

1 子结构快速多极子边界元法

图1 声学域的划分

图1（a）和（b）分别为矩形支管与膨胀腔消声器示意图。按照子结构技术的解决思路，声学域可划分为如图子结构Ω1与子结构Ω2。p1、v1为子结构Ω1真实边界上的声压与质点振速，p2、v2分别为子结构Ω2真实边界上的声压与质点振速。将两个子结构的共用面称为交界面，子结构Ω1的交界面处声压与质点振速为p12与v12，对于子结构Ω2则为p21与v21。按照这样划分，所构成的两个子结构均具有完整的边界，进口、出口以及壁面的边界条件分别为振速边界、阻抗边界和刚性边界，交界面的边界条件为连续性边界，从而应用FMBEM 到每个子结构中［3］。

不同于传统边界元法，FMBEM迭代计算时直接求解矩阵行向量与列向量的积，不构成完整系数矩阵，所以，子结构FMBEM不能通过对子结构矩阵的求逆运算来构建整体矩阵。本文通过叠加各子结构矩阵各行向量积以计算结构整体矩阵的向量积，如式（1）所示：

式中：Hmn与Gmn为离散的子结构单元m中的边界n所形成的系数矩阵。交界面满足连续性条件如下：

联立以上三式，即可求解所有边界上离散节点处物理量。

2 迭代收敛的影响因素

2.1 矩阵构建方式对迭代收敛的影响

鉴于FMBEM迭代求解时无法构建完整的未知量系数矩阵，所以无法应用高斯直接消元法。在迭代算法中，工程中应用最为广泛的是Krylov子空间方法，然而，其未知系数矩阵的条件数是影响其收敛速度的重要因素，因此引入恰当的矩阵预条件处理技术对大型矩阵问题求解就尤为需要，否则迭代次数将使FMBEM的计算速度优势大打折扣。对于子结构FMBEM，其整体系数矩阵为带状，丧失了矩阵对角占优的特点，同样影响着迭代计算的收敛特性，而且，随着模型的改变以及求解方程数目的变化，求解过程的迭代次数也极不稳定，甚至会出现不能收敛的情况。

本文应用BICGSTAB（l）子空间迭代法和ILUT预条件处理技术［4-5］计算子结构FMBEM构建的整体矩阵方程，下面就以图1（a）所示矩形侧支管道内部声场计算为例，分析其迭代收敛特性，并找到影响其收敛的几个因素。

离散节点编号顺序为：（1）子结构Ω1的实边界面；（2）子结构Ω1的交界面；（3）子结构Ω2的实边界面；（4）子结构Ω2的交界面。取列向量构建次序为：P1、P12、P2、V12。构建矩阵方程如式（4）；取列向量构建次序为：P1、P2、P12、V12。构建矩阵方程如式（5）。

采用构建列向量次序为：P1、P12、V12、P2。离散节点编号顺序为：（1）Ω1真实边界面；（2）Ω1交界面；（3）Ω2交界面；（4）Ω2真实边界面。构建的矩阵方程如式（6）所示。另取离散节点编号顺序为：（1）Ω1真实边界面；（2）Ω1交界面；（3）Ω2真实边界面；（4）Ω2交界面。构建形成的矩阵方程如式（7）所示。

选取计算频率为500 Hz，以空气为传播介质，模型结构尺寸：l1= 0.35 m，l2= 0.08 m，c= 0.15 m，a1=a2=b=0.05 m，离散Ω1与Ω2得到总节点数Na为1 812。表1为求解各矩阵方程迭代次数的比较，可见，若要改善收敛速度，应使尽量多的子块矩阵的对角线元素处在总体矩阵的对角线上。此外发现，预处理技术中的数值丢弃阀值和填充的非零元素个数也一定程度影响着迭代收敛次数。

表1 求解各矩阵方程迭代次数

2.2 交界面离散单元数对迭代收敛的影响

以图1（b）所示膨胀腔消声器内部声场计算为例，分析频率为500 Hz，传播介质为空气，结构尺寸为d= 0.049 m、D= 0.164 4 m、l= 0.257 2 m，按照上述要求构建矩阵方程，选取ILUT（10-5，50）预处理方案，根据不同网格尺寸离散其边界，应用子结构FMBEM分析迭代次数，并计算传递损失。

下页图2为迭代次数随交界面离散单元数变化柱形图。由图2可见，交界面的离散单元个数与其迭代收敛次数息息相关，随着单元个数的增加，迭代次数逐渐增多。为了提高计算效率，在划分子结构时，应减少产生大尺寸的交界面；另外，在满足足够的计算精度的前提下，应减少交界面的离散单元数。图3为应用子结构FMBEM与传统BEM在传递损失计算上的比较，结果展现良好的吻合性，也证实子结构FMBEM的准确性。

图2 迭代次数随交界面离散单元数变化柱形图

图3 膨胀腔消声器传递损失曲线

2.3 吸声材料对迭代收敛的影响

在消声器的工程实际应用中，吸声材料的运用十分广泛，而对于含有多种传播介质的声场计算，传统的单一区域FMBEM将无法适用，因而子结构FMBEM的发展使多区域多介质的声场计算成为可能。假设图1（b）所示为膨胀腔阻性消声器，子结构Ω1中为空气，Ω1中填充长纤维玻璃丝绵，取分析频率为500 Hz，结构尺寸与3.2节相同，穿孔管只考虑其支撑吸声材料的作用。考察吸声材料对迭代计算的影响，并计算其传递损失。

取空气密度ρa=1.156 kg/m3，测量到长纤维玻璃丝绵的特性阻抗和波数表达式为：

式中：za为空气的特性阻抗，ka为空气波数。当材料填充密度为200 g/L时，实验测得此材料的流阻率为 17 378 Rayls/m［6］。

图4针对消声器声场中有无吸声材料比较其迭代次数。可见，当声场中填充吸声材料后，由于其阻尼作用，其迭代次数普遍较低，可见将子结构FMBEM应用于阻性消声器声场的计算，更能体现其计算速度上的优势。图5较好地展现了子结构FMBEM与子结构BEM两者在频域计算结果的一致性，证实子结构FMBEM良好的计算精度。

图4 迭代次数随交界面离散单元数变化柱形图

图5 膨胀腔阻性消声器传递损失曲线

3 结 论

本文通过应用子结构FMBEM对各种结构消声器内部声场的计算，在选取恰当的预处理方法下，矩阵方程未知量矩阵向量的构建次序，交界面的选取与离散以及吸声材料的填充均对迭代次数具有一定影响。其中，模型求解过程构建矩阵方程时要注意：一是安排对角占优子矩阵尽量位于总体矩阵的长对角线上；二是采取恰当的预处理方法有效地减少迭代次数；三是过大的交界面会降低迭代收敛速度；四是当求解结构中带有吸声阻尼材料时，子结构FMBEM的迭代次数将会明显降低，且随着频率增大，其迭代特性更为稳定。因而，子结构FMBEM更适用于求解阻性和阻抗复合型消声器三维声场计算。此外，通过子结构FMBEM与BEM在消声器传递损失计算的比较，验证了该方法的正确性与计算精度。

［1］ Nishimura N.Fast multipole accelerated boundary integral equation methods［J］.Appl Mech Rev，2002（4）：299-324.

［2］ Ji Zhen-lin.Acoustic attenuation performance prediction and analysis of perforated tube dissipative silencers［J］.Journal of Vibration Engineering，2005（4）：453-457.

［3］ Sakuma T，Yasuda Y.Fast multipole boundary element method for large-scale steady-state sound field analysis.part I：Setup and Validation［J］.Acta Acustica united with Acustica，2002，88 ：513-525.

［4］ Saad Y ILUT：A dual threshold incomplete LU factorization［J］.Numer.Linear Algebra Appl.，1994（4）：387-402.

［5］ Saad Y.Iterative Methods for Sparse Linear Systems ［M］.Boston：PWS Publishing Company，1996.

［6］ Standard test method for airflow resistance of acoustical materials ［R］.American Society for Testing and Materials ASTM C 522-87，1993，Philadelphia PA.