PML吸収境界条件（Perfectly Matched Layer）
理论与物理
PML是什么
PML是什么的缩写？为什么普通的边界条件不行？
Perfectly Matched Layer——完全匹配层。这是一种人工层，用于在电磁波离开分析区域时，在边界处吸收它们而不产生反射。
诶，固定壁或者对称条件不行吗？结构分析中那样做是可行的。
结构分析处理的是“盒子内的变形”，所以有壁是可以的。但电磁波是向远方辐射的。例如，分析天线的辐射方向图时，如果计算区域的边缘是金属壁，电磁波会全部反射回来。就像在房间里打开扬声器会产生大量回声，导致无法分辨原始声音一样。
原来如此！那在PML出现之前是怎么做的呢？
使用的是诸如Mur边界条件和Engquist-Majda ABC之类的解析吸收边界条件（ABC）。但这些方法有一个致命弱点：对平面波垂直入射有效，但斜入射时会有残留反射。1994年，法国气象局的Bérenger提出了“介质阻抗与内部完全匹配的吸收层”——这就是PML。理论上，它对所有入射角和频率的反射为零。这篇论文在FDTD社区引发了革命，发表30年来被引用超过一万次。
气象局的人发明了FDTD的边界条件？太意外了…！
电导率剖面设计
PML内部具体发生了什么？我想知道电磁波是如何被吸收的。
PML内部设置了人工的电导率 $\sigma(x)$。可以想象这个电导率将电磁波的能量作为焦耳热吸收。关键在于，不是在PML界面处突然设置很大的 $\sigma$，而是从零开始逐渐增加。
PML内的电导率剖面，根据到边界的距离 $x$（$0 \leq x \leq d$, $d$ 是PML厚度）通过多项式梯度定义：
这里 $m$ 是梯度次数（通常 $m = 3 \sim 4$），$\sigma_{\max}$ 是最大电导率。
为什么要逐渐增加？感觉突然增大应该能更快吸收。
问得好。如果在PML界面处 $\sigma$ 从零突然变为一个很大的值，形成“阶跃”，就会产生阻抗失配，从而引发反射。虽然在连续理论上是零反射，但在离散化（网格划分）阶段会产生阶跃。如果缓慢平滑地增加，各网格间的 $\sigma$ 差异就会变小，从而抑制数值反射。就像在高速公路收费站，比起突然停车，在减速车道上慢慢降低速度会更顺畅，对吧？
减速车道的比喻，非常易懂！也就是说 $m$ 的值越大，上升得越慢。
没错。但如果 $m$ 设置得太大，在PML深处还没来得及充分吸收，波就会碰到外壁（PEC）发生反射。$m = 3$ 或 $m = 4$ 是实际应用中的最佳选择。
反射系数理论
$\sigma_{\max}$ 是怎么确定的？随便设大一点就行吗？
$\sigma_{\max}$ 是根据目标反射系数 $R$ 反算出来的。在PML中往返一次后返回的波的反射系数由下式给出：
这里 $\theta$ 是入射角，$\eta_0 = \sqrt{\mu_0/\varepsilon_0} \approx 377\,\Omega$ 是自由空间的阻抗。垂直入射（$\theta = 0$）时 $\cos\theta = 1$，吸收效果最好。
那比如说，想把反射降到 $-60\,\text{dB}$（$R = 10^{-3}$），$\sigma_{\max}$ 应该设为多少？
将上面的公式对 $\sigma_{\max}$ 求解：
例如，设 $m=3$，$d = 10\Delta x$（$\Delta x$ 是单元尺寸），$R = 10^{-3}$，则：
$\sigma_{\max} = \frac{4 \times 377 \times \ln(10^3)}{2 \times 10\Delta x} \approx \frac{4 \times 377 \times 6.91}{20\Delta x} \approx \frac{10416}{20\Delta x} = \frac{520.8}{\Delta x}$ [S/m]
实际应用中，也常用 $\sigma_{\max} \approx 0.8 \times (m+1) / (\eta_0 \Delta x)$ 这个简化公式。$\sigma_{\max}$ 过大会因PML内的急剧变化而增加数值反射，所以通常目标 $R$ 选择在 $10^{-4} \sim 10^{-8}$ 范围内。
复坐标拉伸
能再稍微从数学角度讲讲PML的理论吗？“完全匹配”到底是什么意思？
PML的本质是复坐标拉伸（Complex Coordinate Stretching）。将实坐标 $x$ 替换为复坐标 $\tilde{x}$：
这里 $s_x$ 是拉伸系数。在PML外部（$\sigma = 0$）$s_x = 1$，是普通坐标；在PML内部 $s_x$ 变为复数。只需将麦克斯韦方程组的空间微分 $\partial/\partial x$ 替换为 $\frac{1}{s_x}\partial/\partial x$，就能自动实现对所有入射角和偏振的阻抗匹配。
把坐标本身变成复数，这个想法太天马行空了…
这正是Bérenger的天才之处。物理上相当于“在PML内部，波沿着指数衰减的方向传播”。波不会反射，而是悄然消失——数学上很优美，实现也简单。
PML的变体：Split-field / UPML / CPML
我听说PML也有种类，它们有什么区别？
主要有三个世代：
- Split-field PML（Bérenger原论文, 1994）：将电场和磁场的各分量分割成两个副分量，分别进行衰减。物理直觉易懂，但缺点是变量数量翻倍。
- UPML（Uniaxial PML）（Sacks et al., 1995）：将PML公式化为各向异性介质张量。与频域FEM的亲和性高。不改变麦克斯韦方程组的形式，只修正介质参数，因此易于集成到现有代码中。
- CPML（Convolutional PML）（Roden & Gedney, 2000）：通过卷积积分求解辅助微分方程（ADE）的方式。不像Split-field那样分割变量，能处理任意介质和倏逝波。是当前FDTD实现的事实标准。
| PML类型 | 提出年份 | 主要方法 | 变量增加 | 倏逝波 | 主要应用领域 |
|---|---|---|---|---|---|
| Split-field PML | 1994 | 场分割 | 2倍 | 可能不稳定 | 教育/早期研究 |
| UPML | 1995 | 各向异性张量 | 无 | 可能不稳定 | 频域FEM |
| CPML | 2000 | 卷积+ADE | 仅辅助变量 | 稳定吸收 | FDTD（标准） |
CPML对倏逝波更强，是因为什么不同？
CPML的拉伸系数增加了一个称为 $\kappa$ 的实数项：
$\kappa_x > 1$ 能稳定化倏逝波的衰减，$\alpha_x > 0$ 能改善低频下的PML性能。通过这三个参数（$\sigma, \kappa, \alpha$）的组合，实现了宽带且稳定的吸收。
Bérenger的“尤里卡！”时刻
Jean-Pierre Bérenger曾是法国气象局的雷达研究员。在用FDTD模拟雷电电磁脉冲时，现有的ABC对斜入射的反射太大，无法使用。有一天，他想到一个主意：“如果把电场的每个分量分成两个，并分别给予独立的衰减会怎样？” 物理上这是不自然的“场分割”，但数学上却实现了完美的阻抗匹配。1994年发表在Journal of Computational Physics上的论文最初被怀疑“这种人工层不可能有效”，但一经实现，反射比传统ABC小了2~3个数量级，迅速成为FDTD的标准技术。
PML内的麦克斯韦方程组（频域）
- 拉伸后的法拉第定律：$\nabla_s \times \mathbf{E} = -j\omega\mu_0\mathbf{H}$。这里 $\nabla_s$ 是用拉伸系数 $s_x, s_y, s_z$ 修正后的纳布拉算子。在PML内，将 $\partial/\partial x \to (1/s_x)\partial/\partial x$ 替换。物理上意味着波的传播方向分量加上了复数的衰减。
- 拉伸后的安培定律：$\nabla_s \times \mathbf{H} = j\omega\varepsilon_0\mathbf{E}$。电场和磁场被对称地修正，阻抗 $\eta = \sqrt{\mu_0/\varepsilon_0}$ 得以保持。这就是“完全匹配”名称的由来——PML界面处的阻抗不连续为零。
- UPML介质张量：定义 $[\Lambda] = \mathrm{diag}(s_y s_z/s_x,\; s_z s_x/s_y,\; s_x s_y/s_z)$，则 $\varepsilon_{\mathrm{PML}} = \varepsilon_0[\Lambda]$, $\mu_{\mathrm{PML}} = \mu_0[\Lambda]$。只需改变介质参数就能实现PML，同时保持通常的麦克斯韦方程组形式。
PML设计参数推荐值
| 参数 | 符号 | 推荐范围 | 备注 |
|---|---|---|---|
| PML层数 | $N_{\mathrm{PML}}$ | 8〜16单元 | 越多越好但计算成本增加。10层较实用 |
| 梯度次数 | $m$ | 3〜4 | $m=3$ 最常用 |
| 最大电导率 | $\sigma_{\max}$ | $\frac{0.8(m+1)}{\eta_0 \Delta x}$ | Taflofe经验公式 |
| 目标反射系数 | $R$ | $10^{-4}$〜$10^{-8}$ | $10^{-6}$ 为标准 |
| $\kappa_{\max}$（CPML） | $\kappa_{\max}$ | 5〜15 | 控制倏逝波衰减 |
| $\alpha_{\max}$（CPML） | $\alpha_{\max}$ | $0.01$〜$0.05$ | 低频稳定化。与频率无关 |
数值解法与实现
FDTD中的PML实现
在FDTD中实现PML时，具体要改哪里，怎么改？
在FDTD通常的更新式——例如 $E_z$ 的更新——中加入PML的衰减项。对于CPML，在一维（$x$ 方向的PML）情况下如下所示：
这里 $\Psi_{E_z x}$ 是CPML的辅助变量，通过以下递归式更新：
原来如此，只是在通常的更新式上加上 $\Psi$ 这个辅助变量啊。比想象的要简单。
是的，这是CPML的最大优点。几乎无需改动现有的FDTD代码，只需在PML区域的单元中加入额外的 $\Psi$ 更新即可。内存也只需增加 $\Psi$ 的部分，对计算成本的影响很小。
FEM中的PML实现
FEM的情况呢？和FDTD用同样的方法吗？
FEM通常使用UPML（Uniaxial PML）。对于频域FEM，只需在PML区域的单元中设置各向异性的复介电常数和磁导率张量即可。在建立弱形式时，PML介质张量 $[\Lambda]$ 会自动被纳入。
在基于边单元（Nedelec单元）的公式化中，PML内的旋度方程为：
这里 $k_0 = \omega\sqrt{\mu_0\varepsilon_0}$ 是自由空间的波数。$[\Lambda]$ 是由坐标拉伸决定的复张量，在PML外部退化为单位矩阵。
FEM的网格自由度很高，那也能制作圆柱形或球形的PML吗？
没错。这是FEM的一大优势。HFSS（Ansys）默认会自动生成包围分析区域的球状辐射边界 + PML。长方体PML的角落部分处理起来很麻烦，但球状PML只需根据“到中心的距离”定义拉伸系数，更容易实现匹配。
CPML（Convolutional PML）的公式化
CPML的“卷积”是在卷积什么？和信号处理的滤波器一样吗？
CPML的“卷积”指的是在时域中，通过卷积积分来处理拉伸系数 $s_x$ 的倒数 $1/s_x$ 对场导数的影响。这相当于在时域应用一个递归滤波器（IIR滤波器），从而避免了直接处理复杂的频域关系。是的，其思想类似于信号处理中的递归滤波，目的是高效地实现复坐标拉伸在时域FDTD中的效果。
なった
詳しく
報告