边界元与完美匹配层

壳振动到声压计算——BEM与FEM

当电容器的结构计算已经得到外壳振动后，下一步是求外壳如何推动空气，以及声音传播到指定测点后有多大。本篇只讨论这一步声辐射计算。外壳振动可以来自电磁力、Maxwell 应力、等效本征应变或实测数据；外部声学计算并不要求重新求解这些激励的来源。

常见的两条路线是：

外壳法向振动---边界元法 BEM---外部测点复声压
外壳法向振动---有限元空气域＋PML---空气域内测点复声压

两条路线都可以得到同一测点的频域声压。它们的主要区别是：BEM 通过外壳边界求开放空间声场；空气域＋PML 则在有限大小的空气体积内直接求解，并在外缘吸收向外传播的声波。

1. 两种方法共同需要的输入

结构计算给出外壳表面在频率 下的复位移 。声学计算需要的是外壳如何推动空气，因此先取位移沿外壳法向的分量，再换算成法向复速度。按 的频域约定：

这里 指向外部空气， 表示复数单位。得到 后，就把它作为外壳与空气交界处的声学边界条件：紧贴外壳的空气必须具有相同的法向速度。BEM 和空气域＋PML 从这一边界条件出发，分别计算声波传播到测点后的复声压。COMSOL 声固边界耦合说明

施加边界条件时应保留 的复数值，因为它包含外壳各位置的振动相位；只传入速度幅值会改变声波在测点的叠加结果。

如果给定的是外壳复加速度 ，也可先取其法向分量，再换算为法向复速度：。

2. 路线一：边界元法 BEM

2.1 基本原理

BEM 是 Boundary Element Method，即边界元法。它把外壳表面分成许多小块，计算各小块发出的声波如何到达测点，再把这些贡献叠加。描述“一处的声波怎样传到另一处”的工具就是格林函数。在均匀、无限的三维空气中，它可写成

这里的两个自变量都是空间位置： 是发出声波的源点，通常取外壳上的一个位置； 是接收声波的观察点，例如某个麦克风测点。令两点距离 ，则 描述单位点声源从 到 的传播。分母 表示声波向三维空间扩散时，幅值随距离减小；指数项 表示传播距离造成的相位变化，其中 是波数， 是声速。指数的负号与前面选用的 时间约定一致。对固定的源点，改变观察点就得到不同位置的传播结果；换一个源点，距离和相位也随之改变。

外壳不是单个点声源，而是一整片振动表面。对表面上一小块 ，格林函数给出它传播到观察点的距离和相位权重。为了确定这一小块实际贡献多大，还需要该处的边界声压 和声压法向梯度 。格林恒等式把所有小块的贡献相加，得到

这不是“知道边界上一点的声压，就直接推出测点声压”：必须综合整个外壳边界的信息。已知外壳法向速度后，声学边界条件确定声压法向梯度；BEM 用边界积分方程求出尚未知的边界声压。此时两项边界数据都已具备，将某个测点坐标代入上式并对整个外壳积分，就得到该点的复声压。复数相加会自动保留各处声波到达测点时的相位差。对17个测点重复评估即可，不需在测点之间绘制空气网格。COMSOL BEM 建模介绍、COMSOL BEM 后处理说明

上式的 是外部空气计算域的外法向；在外壳这个内边界上，它与前文定义的 方向相反。实际在 COMSOL 中应按所选“法向速度”边界条件的正方向输入，避免符号颠倒。

2.2 反射地面怎样处理

若地面近似为无限大声硬平面，则地面上法向声速为零。BEM 可以利用对称条件，或者等效地引入镜像声源。以 为地面时，声硬平面的示意性格林函数为

原声源和镜像声源的声压在地面同相叠加，从而满足声硬条件。实际有限尺寸地面、吸声地面或试验室墙面不能仅靠这一理想镜像条件代表。

2.3 BEM 的特点

BEM 无需绘制并网格化几米甚至几十米的外部空气域，开放空间辐射条件由积分公式体现。它特别适合均匀空气中的开放空间辐射、测点距设备较远、反射面可用理想边界表示的情况。其代价是边界积分通常产生较密的线性系统；边界网格很多时，内存和计算时间也可能迅速增加。BEM 对外部空气的空间变参数、复杂局部流动或大量实体障碍物，也不如空气体积有限元直接。COMSOL BEM 与 FEM 适用场景

2.4 为什么叫边界元

“边界元”指的是：数值计算主要把边界划分成小单元，并在这些单元上求解未知的声学量。对于本文的外部声场，已知外壳各处怎样振动，再结合空气的声速、密度和描述声波传播的格林函数，就能建立边界积分方程，求出外壳边界上的声压。然后将整个边界各单元的传播贡献积分相加，便可计算外部任一测点的声压。

它的核心思路是：**先求清楚边界如何发声，再利用已知的传播规律把边界上的声场延伸到空间中。**因此不必给外部空气的每一小块体积都划网格、逐点求解；但仍必须明确空气介质、传播条件及所有相关边界。对这台电容器而言，外壳是发声边界，17个测点只是积分结果的评价位置。

3. 路线二：有限元空气域＋PML

3.1 空气域负责求解什么

FEM 是 Finite Element Method，即有限元法。先在电容器外建立一个包含17个测点的真实空气域，并把空气体积分成许多相连的小单元。未知量是空气中的复声压 ：它在空间各处的分布，必须满足给定频率下的声波方程；在外壳边界，还必须符合已知的外壳法向速度。

有限元在每个单元内用简单的形函数近似未知声压，例如

这里的 是当前这个网格单元的节点数， 逐一指向这些节点。它不是整个空气模型的节点数，也不是17个麦克风测点数。 是该单元第 个节点的待求复声压； 是该节点对单元内部位置 的插值权重。

先看一个一维例子，便于理解插值。设一段长度为 的空气被划为一个线单元，左右端点分别位于 和 ，端点复声压为 、。这个单元有两个节点，即 。若采用一次形函数，单元内的声压近似为

在左端 ，结果就是 ；在右端 ，结果就是 ；在四分之一处 ，结果是 。这说明：知道单元节点的复声压后，就可以通过形函数估计单元内部任意位置的复声压。此处的复数必须直接参与加权，不能先把两个声压取绝对值。

再看一个真正的二维、二次插值例子。取三角形单元，三个顶点为节点1 、节点2 、节点3 ；三条边的中点依次为节点4 、节点5 、节点6 。因此这个单元有 个节点。这里的坐标是为演示插值而选的局部坐标。

先定义 ，它们叫三角形的面积坐标，可以把一个内部点看成由三个顶点按权重组合而成。节点1、2、3的坐标分别为 、、，所以对内部点 有

比较横、纵坐标便得到 、；再由三者之和为1，得到 。在节点1，三个 的值是 ；在节点2是 ；在节点3是 。例如边1—2的中点（节点4）对应 。这些 本身只是表示点在三角形中的位置，还不是二次形函数。

现在构造六个 。规则是： 在自己的节点等于1，在另外五个节点等于0，即让这六个点都位于插值函数上。先以顶点1的 为例。节点1处 ，两个相邻边中点处 ，其余三个节点处 。要用 构造一个二次式，使它在 和 时为0、在 时为1，便得到 。对顶点2和3重复同样的构造，得到 、。

再看边1—2中点的 ：它在该中点要等于1，在其他节点等于0。 在边1—2的两端及另外两条边的节点上为0；在该边中点则为 ，所以乘以4得到 。另外两条边按同样方法得到 、。六个二次形函数汇总为

于是单元内部任意点 的声压由六个节点共同决定：。例如在三角形重心 ，有 ，三个顶点的权重各为 ，三个边中点的权重各为 ：

为检验这个插值，假设六个节点的数值恰好来自 ：它们依次为 。代入上式得到重心声压 ，与直接计算 一致。这里的数值只用于说明二次插值；真实声学计算中，六个 是由方程组求出的复数。某些二次形函数在单元内部可以取负值，重心处的 属于正常的插值权重。

实际电容器外部空气是三维的。例如一个采用一次形函数的四面体单元有4个顶点节点，即 ；如果采用二次形函数，节点还会包括边上的点。17个测点分别落在某些空气单元内，求解后按各自所在单元的形函数取值。

相邻单元共用节点，因此它们的声压场连接起来。把单元内的近似声压代入声波方程的积分形式，各单元会给出一小组方程；再按共用节点将所有单元的方程组装起来，得到整个空气域的矩阵方程：

矩阵 描述空气中声压如何相互影响，右端 包含外壳法向振动施加的边界激励。求解一次方程组，就得到该频率下各节点的复声压。如果测点不恰好落在节点上，再用它所在单元的形函数插值，即可得到测点复声压。对每个需要的频率分别计算，便得到测点的声压频谱。COMSOL 6.1 声学模块手册

因此，FEM 与前述 BEM 在计算位置上的区别很直观：FEM 在空气体积内设置未知声压并求解；BEM 主要在边界上设置未知量，再由传播积分计算空气中的测点声压。

只截取有限空气域会带来一个问题：声波到达模型外边界后，普通边界条件可能使它反射回来，污染测点结果。于是需要在真实空气域外再包一层 PML。

3.2 PML 吸收声波的原理

PML 是 Perfectly Matched Layer，即完美匹配层。它不是现实中的吸声材料，而是计算中对波动方程做特殊的复坐标伸缩：理想情况下，向外传播的声波穿过真实空气与 PML 的交界面时没有反射，进入 PML 后幅值逐渐衰减。这样有限大小的模型就能近似无限开放空间。COMSOL 声学模块手册：PML 与网格

模型布局应为：

振动外壳 │ 真实空气域（17点放在这里）│ PML │ 计算域最外边界

测点不要放在 PML 内，因为 PML 中的场已被人为衰减，它不是要报告的物理空气声压。

3.3 网格与计算成本

空气体积的网格必须解析最高频率的波长：

以 为例，400 Hz 时波长约0.86 m，4000 Hz 时只有约0.086 m。若空气域包住距离电容器约1 m的所有测点，频率越高，体网格需要越密；PML 自身也需要合适厚度和网格层数。COMSOL 的声学网格说明给出按每波长若干个二阶单元控制网格的做法，具体精度仍应由网格收敛检验。COMSOL 6.1 声学模块手册

空气域＋PML 的优点是整个真实空气域内都有直接求解的声压，便于查看声压云图、局部反射和复杂空间结构。当空气域不大、网格能装入内存时，稀疏有限元系统有时比 BEM 更快。COMSOL 方法比较

4. 怎样从复声压得到测点 SPL

两种方法最终都应给出每个测点 、每个频率 的复声压幅值 。若复幅值采用本文的峰值幅值约定，单条频线的均方根声压是

以空气中的参考声压 计算声压级：

如果直接使用 COMSOL 已经定义好的声压级变量，例如模型中的 pabe.Lp，不要再额外除以 ；只有从原始复压力自行计算时才需要核对该变量的幅值约定。

不同频率的线谱在时间平均意义下按均方声压相加，而不是将 dB 数值直接相加：

这里 必须明确。例如只算100、200、300、400 Hz，得到的只是四条频线合成声压级，不是100–4000 Hz的完整总声压级。若同一频率存在多个声源，须先将它们的复声压相加，再取幅值；不能把同频声源当作互不相干的不同频率相加。

17点的能量平均声压级可按

计算。这个量是指定测量面上的平均声压级；声功率级还需要测量面积及相应评价方法的环境修正，不能把某一个测点的 SPL 直接称为声功率级。

5. 怎样选择两种方法

条件或目标	较合适的方法	原因
单台电容器向均匀、开放空气辐射；需要1 m外的17点	BEM	无需建立覆盖测点的三维空气体网格
计算机内存较少，同时需要较宽频率范围	优先评估 BEM	高波数下空气体网格可能很大；仍需检查BEM边界网格规模
空气域较小，希望直接显示整个域内声压	空气域＋PML	有限元直接给出真实空气域内的场
外部空间存在具体墙面、挡板、设备或局部材料变化	空气域＋PML，或混合方法	可以显式建立这些对象与非均匀区域
外部是大范围均匀空气，但设备附近有复杂局部空间	局部空气有限元＋外部 BEM	各方法分别处理适合的区域

对当前电容器，内部油是有限封闭域，适合用有限元；外部若按均匀空气、理想反射地面和开放空间处理，BEM 是合理的选择。BEM 与 PML 都不能修复一个为零的外壳速度：若声学边界没有有效振动，两种方法都会得到零声压。要比较两种方法的数值精度，应给它们完全相同的复法向速度、空气参数、地面条件和测点位置，再分别做网格收敛检查。COMSOL 声学 BEM 实例

6. 总结

1. 结构输出是外壳复法向速度；声学输出是测点复声压。相位对声场很重要。
2. BEM 在边界上求开放声场；空气域＋PML 在有限空气体积中求声场，并在外缘吸收出射波。
3. SPL 必须由有效的复声压计算；不同频率按均方声压合成，所报“总声压级”必须注明覆盖的频率范围。