3. 有限元声学离散原理

声学有限元的核心任务，是把连续介质中无法直接求解的波动方程，转化为计算机可以求解的有限维代数方程组。这一转化过程由四个关键环节构成：用加权残差法把强形式方程变为弱形式积分方程；用形函数把连续声压场离散为节点自由度；在单元层面计算声学质量矩阵与刚度矩阵；最后把所有单元矩阵按节点编号组装成全局系统方程。理解这条链路，才能真正看懂仿真软件内部在做什么，也才能在结果异常时判断问题出在哪个环节。

从强形式到弱形式：Galerkin加权残差法

考虑稳态声场，声压满足亥姆霍兹方程。在区域 Ω 内写为强形式：

∇²p + k²p = 0        (Ω 内)
∂p/∂n = -jωρ₀vₙ      (边界 Γ 上)

其中 k = ω/c 为波数，ρ₀ 为介质密度，vₙ 为边界法向振速。强形式要求声压二阶可导，对复杂几何几乎无法直接构造近似解。Galerkin法的思路是：先假设近似解 p̃ = Σ Nᵢ pᵢ，代入方程后必然产生残差 R = ∇²p̃ + k²p̃。我们要求残差在某种加权意义下为零，取权函数与形函数相同（即 Galerkin 选择），得到加权残差积分：

∫Ω Wᵢ (∇²p̃ + k²p̃) dΩ = 0

对二阶导数项使用格林公式（分部积分），把二阶导转移到权函数上，同时自然引入边界法向导数：

∫Ω ∇Wᵢ·∇p̃ dΩ - k²∫Ω Wᵢ p̃ dΩ = ∫Γ Wᵢ (∂p̃/∂n) dΓ

这一步是有限元声学离散的关键转折：它把对声压的苛刻可导要求降低为一阶可导，同时把边界振速条件以自然边界条件的形式纳入方程右端。左端两项分别对应声学刚度项和声学质量项，右端则是外部激励项。

形函数构造：把连续场变成节点自由度

形函数决定了单元内声压如何由节点值插值。以一维两节点线性单元为例，单元长度 L，局部坐标 ξ ∈ [-1, 1]，两个形函数为：

N₁(ξ) = (1 - ξ)/2
N₂(ξ) = (1 + ξ)/2

它们满足 Nᵢ(ξⱼ) = δᵢⱼ，即节点处取 1、另一节点处取 0，且 N₁ + N₂ = 1，保证常值声压可被精确表示。单元内声压 p(ξ) = N₁p₁ + N₂p₂。二维三角形单元则用面积坐标 L₁、L₂、L₃ 作为形函数，三维四面体单元用体积坐标，形式完全类似。形函数的阶次决定单元精度：线性单元每个波长至少需要约 6～10 个单元才能控制色散误差，二次单元可适当放宽，但节点数增加。

单元矩阵：质量矩阵与刚度矩阵

把 p̃ = N pᵉ 代入弱形式，对每个单元得到两个基本矩阵。声学刚度矩阵来自 ∇N·∇N 项：

Kᵉ = ∫Ωᵉ (∇N)ᵀ(∇N) dΩ

声学质量矩阵来自 N·N 项，并带有 -k² 因子：

Mᵉ = ∫Ωᵉ NᵀN dΩ

对一维线性单元，解析积分结果为：

Kᵉ = (1/L) [ 1  -1 ]
          [-1   1 ]

Mᵉ = (L/6) [ 2   1 ]
           [ 1   2 ]

可以看到刚度矩阵只与几何有关，质量矩阵与单元尺寸成正比。二维三角形线性单元的质量矩阵常用集中质量近似，把总质量平均分配到三个节点，得到对角矩阵，可显著降低求解规模并改善条件数。实际编程中，这些积分通常用高斯数值积分完成，一维两点高斯即可精确积分线性单元，二维三角形用三点或七点积分。

全局系统方程的形成

单元矩阵计算完成后，需要按节点编号把各单元贡献叠加到全局矩阵中。设全局节点声压向量为 p，组装规则是：单元第 i 个局部节点对应全局节点编号 g(i)，则 Kᵉ(i,j) 累加到 K(g(i), g(j))，Mᵉ 同理。所有单元遍历完成后，得到全局方程：

(K - k²M) p = f

其中 f 为边界振速激励向量，由 ∫Γ Nᵀ(-jωρ₀vₙ) dΓ 计算。若存在声学软边界（p = 0）或阻抗边界，还需在组装后修改对应自由度或增加阻尼项。最终这是一个大型稀疏复系数线性方程组，规模可达百万自由度，通常用稀疏直接求解器或迭代求解器处理。

组装与边界条件的实际处理

组装过程在编程中通常用一个“单元-节点映射表”驱动。以两个一维单元、三个节点为例，单元1连接节点1-2，单元2连接节点2-3。单元1的 Kᵉ 贡献到全局 K 的 (1,1)、(1,2)、(2,1)、(2,2) 位置；单元2贡献到 (2,2)、(2,3)、(3,2)、(3,3)。节点2同时属于两个单元，其对角项被累加两次，这正是有限元“节点共享”的体现。组装完成后，若节点1为刚性壁面（vₙ=0），则 f 中对应项为零；若节点3为声压释放边界（p=0），则需把第3行第3列处理为单位对角，右端置零，或采用罚函数法近似施加。

一个可复现的教学算例

设一维声管长 1 m，截面积 0.01 m²，空气密度 ρ₀ = 1.21 kg/m³，声速 c = 343 m/s。左端为振速激励 vₙ = 0.01 m/s，右端为声压释放边界。用 10 个线性单元离散，单元长度 L = 0.1 m。每个单元的质量矩阵和刚度矩阵按上述公式计算，组装后得到 11×11 的全局矩阵。施加右端 p₁₁ = 0 后，求解复系数方程组，可得到各节点声压幅值。在 500 Hz 下，理论一维驻波声压分布呈余弦形，数值解与解析解的相对误差应小于 2%；若把单元数减到 4 个，误差会明显增大，说明网格密度必须与最高分析频率匹配。

掌握这条从强形式到全局方程的完整链路后，面对仿真软件中的“声学单元类型”“质量矩阵形式”“边界条件施加方式”等选项，就能判断它们分别作用于哪个环节，也能在结果异常时快速定位是形函数阶次不足、积分精度不够，还是组装或边界处理出了问题。