4. 边界元声学离散原理

边界元法(BEM)求解声学问题的核心思路,是把无限域中的Helmholtz方程通过格林第二定理转化为边界上的积分方程,从而将三维问题降为二维边界问题。与有限元法需要在整个声场域内划分体网格不同,边界元只需在辐射体或散射体表面划分面网格,天然满足Sommerfeld辐射条件,特别适合外声场、无限域和开放边界问题。本节围绕边界元声学离散的四个关键环节展开:直接与间接边界积分方程、单层与双层势、奇异积分处理以及系数矩阵填充。

直接边界积分方程与间接边界积分方程

直接边界积分方程以边界上的声压 p 和法向速度 vn 为未知量,物理意义明确。对声场中任意点 x,Helmholtz方程的基本解(自由空间格林函数)为:

G(x, y) = e^(-jk|x-y|) / (4π|x-y|)

其中 k 为波数,y 为边界上的源点。将格林第二定理应用于声场与格林函数,可得到外声场的直接边界积分方程:

C(x)p(x) = ∫_S [ p(y) ∂G/∂n_y - G(x,y) ∂p/∂n_y ] dS(y)

式中 C(x) 为立体角系数,当 x 位于光滑边界上时取 0.5,位于声场内部时取 1,位于外部时取 0。直接法的未知量就是边界上的声压和法向速度,物理量直接可读,适合与结构有限元耦合。

间接边界积分方程则引入单层势 σ 和双层势 μ 作为中间变量,通过势函数叠加构造声场。间接法在求解薄壁结构、开口结构或内外声场耦合时具有优势,因为单层势和双层势可以分别对应边界两侧的声压跳跃和法向速度跳跃。两种方法在数学上等价,但数值实现和适用场景不同:直接法更适合封闭辐射体,间接法更适合薄壳、隔板和多域耦合问题。

边界元声学离散:直接法与间接法对比 直接边界积分方程 边界 S x (场点) y (源点) G(x,y) n 未知量:边界声压 p 与法向速度 vₙ C(x)p(x) = ∫[p(y)∂G/∂n − G·∂p/∂n] dS C = 0.5 (光滑边界) / 1 (域内) / 0 (域外) 间接边界积分方程 边界 S σ (单层势) μ (双层势) 未知量:单层势 σ 与双层势 μ p(x) = ∫[σ(y)G − μ(y)∂G/∂n] dS σ ↔ 声压跳跃,μ ↔ 法向速度跳跃 适用场景对比 直接法:封闭辐射体、结构-声耦合 物理量直接可读,适合与FEM联立 间接法:薄壁、开口、多域耦合 势变量灵活,适合处理跳跃条件

单层势与双层势的物理含义

单层势对应边界上分布的“声压型”源,其物理图像是边界上每一点放置一个脉动球源,声场由这些球源叠加而成。单层势在穿过边界时连续,但其法向导数发生跳跃,跳跃量等于单层势密度 σ。双层势对应边界上分布的“偶极子”源,其物理图像是边界上每一点放置一个偶极子,声场由偶极子叠加而成。双层势在穿过边界时发生跳跃,跳跃量等于双层势密度 μ,而其法向导数连续。

在间接边界元中,边界两侧的声压差和法向速度差分别由 μ 和 σ 控制。对于薄壁结构,壁面两侧声压差直接对应双层势密度,这使得间接法在处理隔声、透射和薄板振动辐射时非常自然。实际实现时,通常将 σ 和 μ 用边界单元上的形函数展开,再通过配点法或伽辽金法建立代数方程组。

奇异积分的来源与处理策略

边界积分方程中,当源点 y 与场点 x 重合时,格林函数 G 出现 1/r 型奇异,其法向导数 ∂G/∂n 出现 1/r2 型强奇异。这是边界元法数值实现中最关键的难点之一。若直接数值积分,奇异点附近积分值会发散,导致系数矩阵元素错误。

常用的处理策略包括:

  • 坐标变换法:将奇异单元通过极坐标变换或三角形细分,使雅可比行列式引入 r 因子,从而抵消 1/r 奇异。对弱奇异积分,极坐标变换后积分变为常规可积形式。
  • 刚体位移法:利用声场在单位声压下的刚体运动条件,间接确定强奇异积分的自由项系数。该方法不需要直接计算强奇异积分,而是通过矩阵行和为零的条件反推对角元素。
  • 正则化法:将强奇异积分分解为可积部分和奇异部分,奇异部分通过解析公式计算,可积部分用常规高斯积分处理。
  • 自适应细分:在奇异点附近自动加密积分点,适用于高阶单元和复杂几何。
对于 1/r 弱奇异积分,极坐标变换是最直接有效的方法;对于 1/r2 强奇异积分,刚体位移法在工程实现中更为稳健,尤其适合与直接法配合使用。

系数矩阵填充与数值实现

边界元离散后,系数矩阵 H 和 G 的填充是计算量最大的环节。以直接法为例,对每个边界节点 i 和每个单元 j,需要计算:

H_ij = ∫_Sj N_j(y) ∂G(x_i, y)/∂n_y dS(y)
G_ij = ∫_Sj N_j(y) G(x_i, y) dS(y)

其中 Nj 为单元形函数。对于非奇异单元,采用标准高斯积分即可;对于奇异单元,按上述策略处理。矩阵填充完成后,结合边界条件(如法向速度已知、声压已知或阻抗边界条件),将未知量移至方程左侧,形成最终线性方程组 Ax = b。

实际编程时,通常采用以下流程:先循环所有场点,再循环所有单元,判断场点是否落在当前单元上,若是则调用奇异积分处理函数,否则调用常规高斯积分。对于大规模问题,可采用快速多极子算法(FMM)或层次矩阵(H-matrix)加速矩阵向量乘积,避免显式存储稠密矩阵。

示例参数(教学假设):半径 0.5 m 的脉动球,表面法向速度 0.01 m/s,频率 500 Hz,空气密度 1.21 kg/m³,声速 343 m/s。边界离散为 320 个三角形线性单元,节点数 162。采用直接边界元法,配点法建立方程。计算得到表面声压幅值约 2.3 Pa,远场声压级约 85 dB(参考 20 μPa)。该结果用于演示矩阵填充与奇异积分处理流程,非实测数据。

在系数矩阵填充过程中,需要特别注意单元法向方向的一致性。所有单元法向必须指向声场外部(或统一指向内部,取决于方程推导约定),否则积分符号会出错。建议在网格生成后统一检查法向方向,并在代码中固定法向定义,避免因方向混乱导致结果异常。

数值验证与常见问题

完成矩阵填充后,可通过以下方式验证实现正确性:

  • 解析解对比:对脉动球、平面波等有解析解的问题,比较边界元计算结果与解析解,误差应随网格加密而收敛。
  • 刚体位移检验:在单位声压、零法向速度条件下,声场应为均匀场,检验矩阵行和是否接近零。
  • 互易性检验:交换源点和场点,结果应满足互易定理。
  • 频率扫描:检查是否存在非物理的伪模态(对应内部问题特征频率),必要时采用CHIEF法或Burton-Miller法消除。

常见问题包括:奇异积分处理不当导致矩阵元素异常大或异常小;单元法向不一致导致符号错误;网格过粗导致高频结果失真;内部特征频率处出现伪峰。针对这些问题,应优先检查奇异积分处理函数和法向定义,再考虑网格加密或采用更稳定的积分方程形式。

边界元声学离散的核心在于正确理解直接法与间接法的物理含义,掌握单层势与双层势的跳跃条件,并针对奇异积分采取合适的数值处理策略。系数矩阵填充是计算量最大的环节,需要结合网格质量、积分精度和加速算法综合优化。只有把这几环节都做扎实,边界元法才能在工程声学仿真中发挥其降维、无限域和辐射条件自动满足的优势。