吸收边界条件（PML）的设计与实现

在数值电磁仿真中，有限计算域的边界处理是影响结果精度的关键因素。
吸收边界条件（Perfectly Matched Layer，简称 PML）通过在计算域外围构造一层人工介质，使进入该层的电磁波在数值上实现指数衰减，从而模拟出无穷空间的开放边界。下面系统地介绍 PML 的理论依据、常见实现方式以及在 FDTD（时域有限差分）中的具体实现步骤，并给出完整的代码示例和参数调优建议。

• 匹配性（Perfect Matching）
  PML 通过在 Maxwell 方程中引入复坐标拉伸，将真实空间的波数向量 k 替换为
  [
  \tilde{k}_i = \frac{k_i}{1 + j\sigma_i / \omega},\quad i = x,y,z
  ]
  其中 (\sigma_i) 为各方向的吸收系数，(\omega) 为角频率。该变换保证在 PML 与主计算域交界面上阻抗完全匹配，理论上不存在反射。

• 指数衰减
  在 PML 区域，电场与磁场的幅值随深度 (d) 按 (\exp(-\alpha d)) 衰减，(\alpha) 与 (\sigma_i) 成正比。只要 (\sigma_i) 选取得当，波在有限厚度内即可衰减到机器精度以下。

• 各向异性介质等效
  拉伸坐标等价于在 PML 区域引入各向异性电导率张量，使得 Maxwell 方程仍保持原始形式，只是介质参数被修改为复数值。

2. 常见的 PML 实现方式

本文重点演示 UPML 在二维 TM 模式 FDTD 中的实现，因为它兼具实现简洁和数值稳定的优势。

3. 在 FDTD 中实现 UPML 的步骤

1. 确定 PML 区域厚度
   通常取 8~12 个网格单元；厚度越大衰减越充分，但会增加计算量。

2. 构造吸收系数分布
   采用多项式递增方式：
   [
   \sigma_i (n) = \sigma_{\max}\left(\frac{n}{N_{\text{PML}}}\right)^m,\quad n=1,\dots,N_{\text{PML}}
   ]

   • (\sigma_{\max} = \frac{(m+1)}{150\pi\Delta})（经验公式），(\Delta) 为网格间距。
   • 指数 (m) 常取 3~4，可根据需求微调。

3. 计算更新系数
   对于 TM 模式（(E_z, H_x, H_y)），在 PML 区域需要分别引入 (\kappa) 与 (\alpha) 参数：
   [
   \kappa_i = 1 + (\kappa_{\max}-1)\left(\frac{n}{N_{\text{PML}}}\right)^m,\quad
   \alpha_i = \alpha_{\max}\left(1-\frac{n}{N_{\text{PML}}}\right)
   ]
   其中 (\kappa_{\max}) 常设为 1，(\alpha_{\max}) 取 (\pi f_{\text{max}})（最高工作频率）。

4. 修改更新方程
   以 (E_z) 为例，原始更新式为
   [
   E_z^{n+1}(i,j)=E_z^{n}(i,j)+\frac{\Delta t}{\varepsilon}\left[\frac{H_y^{n+1/2}(i,j)-H_y^{n+1/2}(i-1,j)}{\Delta x}
   -\frac{H_x^{n+1/2}(i,j)-H_x^{n+1/2}(i,j-1)}{\Delta y}\right]
   ]
   在 PML 中加入分量 (\psi)（卷积项）并使用 (\sigma,\kappa,\alpha) 调整系数，得到
   [
   \begin{aligned}
   \psi_{E_x}^{n+1} &= b_{E_x}\psi_{E_x}^{n}+c_{E_x}\left(H_y^{n+1/2}(i,j)-H_y^{n+1/2}(i-1,j)\right)\
   E_z^{n+1} &= a_{E_z}E_z^{n}+d_{E_z}\left(\psi_{E_x}^{n+1} - \psi_{E_y}^{n+1}\right)
   \end{aligned}
   ]
   其中系数 (a,b,c,d) 由 (\sigma,\kappa,\alpha) 预先计算得到。

5. 边界条件的统一处理
   在 PML 与主域交界处直接使用相同的网格更新公式，无需额外的匹配层，保证“完美匹配”。

4. 示例代码（二维 TM‑FDTD + UPML）

% -------------------------------------------------
% 2D TM FDTD with UPML (MATLAB 示例)
% -------------------------------------------------
clear; clc;

% ---------- 参数 ----------
c0   = 299792458;               % 真空光速
eps0 = 8.854187817e-12;
mu0  = 4*pi*1e-7;
f0   = 3e9;                     % 工作频率
lambda0 = c0/f0;
dx = dy = lambda0/40;          % 空间步长
dt = 0.99/(c0*sqrt(1/dx^2+1/dy^2)); % 时间步长 (CFL 条件)

Nx = 200; Ny = 200;            % 主域尺寸
npml = 12;                     % PML 厚度

% ---------- 介质参数 ----------
eps = eps0*ones(Nx,Ny);
mu  = mu0 *ones(Nx,Ny);

% ---------- PML 参数 ----------
m = 4;                         % 多项式阶数
sigma_max = (m+1)/(150*pi*dx);
sigma_x = zeros(Nx,1);
sigma_y = zeros(Ny,1);
for i=1:npml
    sigma_x(i) = sigma_max*((npml-i+1)/npml)^m;
    sigma_x(Nx-i+1) = sigma_x(i);
    sigma_y(i) = sigma_max*((npml-i+1)/npml)^m;
    sigma_y(Ny-i+1) = sigma_y(i);
end
% 生成二维 sigma 矩阵
Sigma_x = repmat(sigma_x,1,Ny);
Sigma_y = repmat(sigma_y.',Nx,1);

% ---------- 预计算系数 ----------
% 对 Ez、Hx、Hy 分别计算 a,b,c,d 系数
% 这里仅给出 Ez 的系数示例
aEz = (1 - dt*Sigma_x./(2*eps))./(1 + dt*Sigma_x./(2*eps));
bEz = dt./(eps.*dx)./(1 + dt*Sigma_x./(2*eps));
cEz = dt./(eps.*dy)./(1 + dt*Sigma_y./(2*eps));

% ---------- 场变量 ----------
Ez = zeros(Nx,Ny);
Hx = zeros(Nx,Ny-1);
Hy = zeros(Nx-1,Ny);

% ---------- 主循环 ----------
nSteps = 1000;
for n = 1:nSteps
    % 更新 Hx
    Hx(:,1:end) = Hx(:,1:end) - (dt/mu(:,1:end)).*(Ez(:,2:end)-Ez(:,1:end-1))/dy;
    % 更新 Hy
    Hy(1:end,:) = Hy(1:end,:) + (dt/mu(1:end,:)).*(Ez(2:end,:)-Ez(1:end-1,:))/dx;
    
    % 更新 Ez（含 PML 系数）
    curlH = (Hy(2:end,:)-Hy(1:end-1,:))/dx - (Hx(:,2:end)-Hx(:,1:end-1))/dy;
    Ez(2:end-1,2:end-1) = aEz(2:end-1,2:end-1).*Ez(2:end-1,2:end-1) ...
                         + bEz(2:end-1,2:end-1).*curlH;
    
    % 软源激励（高斯脉冲）
    t0 = 30; spread = 10;
    Ez(Nx/2,Ny/2) = Ez(Nx/2,Ny/2) + exp(-((n-t0)/spread)^2);
    
    % 可视化（每 20 步绘图一次）
    if mod(n,20)==0
        imagesc(Ez.');
        colorbar; axis equal tight;
        title(['Ez at step ',num2str(n)]);
        drawnow;
    end
end

说明

• 代码中仅展示了 Ez 的 PML 系数计算，Hx、Hy 的系数可类比得到。
• Sigma_x、Sigma_y 在主域内部为零，仅在 PML 区域取非零值，实现“完美匹配”。
• 通过调节 npml、m 与 sigma_max，可以在保持计算效率的前提下将边界反射降低到 -60 dB 以下。

5. 参数选取与调优技巧

1. PML 厚度

   • 对于低频（波长较大）问题，建议使用 10~15 层；
   • 高频（波长接近网格尺度）时 8~10 层已足够。

2. 吸收系数的多项式阶数 (m)

   • (m=3) 时衰减曲线平缓，适合宽频带；
   • (m=4) 或 (5) 能在较短厚度内实现更快衰减，但对网格不均匀性更敏感。

3. 最高工作频率

   • 若仿真包含多频成分，(\sigma_{\max}) 需基于最高频率计算，以防高频分量出现显著反射。

4. 数值稳定性

   • 在极端吸收系数下，时间步长 dt 可能需要适当降低（如 0.95 CFL），防止因复数系数导致的数值不稳定。

5. 验证方法

   • 采用平面波入射或点源辐射的基准案例，比较 PML 前后场强的衰减曲线；
   • 计算反射系数 (R = |E_{\text{ref}}/E_{\text{inc}}|)，目标值一般要求低于 (10^{-3})（约 -60 dB）。

6. 小结

• PML 通过复坐标拉伸实现了理论上的零反射，是现代电磁数值仿真中最常用的开放边界处理技术。
• UPML 在时域 FDTD 中的实现相对直接，只需在每个方向上引入电导率 (\sigma) 与伸缩系数 (\kappa)，并在更新方程中加入相应的加权系数。
• 参数的合理选取（厚度、吸收系数分布、最高工作频率）决定了 PML 的衰减效率与数值稳定性。
• 通过本文提供的 MATLAB 示例，读者可以快速搭建起带 PML 的二维 FDTD 仿真平台，并在实际工程问题中进一步调优以满足精度需求。

掌握 PML 的设计与实现，不仅能显著提升仿真结果的可靠性，还为后续的多物理场耦合、非线性材料以及宽频带天线分析奠定坚实的数值基础。