简介这组MATLAB代码聚焦五种完美匹配层PML实现与效果对比适用于FDTD等数值仿真中需要抑制边界反射的电磁波、声波或弹性波问题也适合刚接触吸收边界条件的研究者或高年级本科生对照学习。压缩包共23个文件包含12个.m脚本、10个.mat数据文件和1份PDF文档体积仅301KB代码中分别实现了Smith-PML、Steklov-Poincaré PML、多层PML、自适应PML与改进的Bermúdez PML等不同思路并配有扫参脚本和反射计算函数可直观比较反射系数、衰减效果与计算效率。目前已有631人学习包内附有PML Comparison.pdf对比说明和对应的TMz仿真数据便于快速复现五类PML在12GHz中心频率、24GHz带宽下的表现可为实际建模中选取合适PML提供直接参考。1. 为什么要在MATLAB里对比5种PML做二维FDTD电磁波仿真的工程师多半已经学会了“在区域四周加PML”来压缩反射既然PML在理论上接近完美吸收为什么还要对比因为“接近完美”只在离散网格、垂直入射、有限带宽这三个条件下成立。换一种PML实现改一行系数可能在斜入射、低频或长时间步进时反射误差相差两个数量级而这类误差会直接进入后续的近场、远场提取结果。这里用同一个2D TE波FDTD框架把Berenger分裂场PML、UPML、CPML、NPML以及一个经常被误认成PML的指数衰减层同时实现给出可运行的MATLAB代码、参数设置和量化对比方法。这组内容不是要告诉你“哪个最好”而是给出一个能自己复现、能在论文里画对比曲线的基准。2. 五种PML的原理选型和MATLAB关键递归2.1 Berenger分裂场PML作为对比基准的经典实现1994年Berenger提出的分裂场PML核心思路是把边界区域内的场分量按坐标方向劈裂让波在不同方向拥有不同的虚衰减率从而在没有真实反射界面的情况下把外向波“吞掉”。在二维TE波Ez、Hx、Hy中常见做法是对Hx和Hy各保存两个子分量Ez保持单一场靠不同方向上的电导率σ_x、σ_y驱动衰减。% 以x方向PML区域内的Hy为例Hy拆成Hy_x与Hy_z Hy_x(i, j) b_x(i) .* Hy_x(i, j) - c_x(i) .* (Ez(i1, j) - Ez(i, j)) / dx; Hy_z(i, j) Hy_z(i, j) (Ez(i, j1) - Ez(i, j)) / dy; Hy(i, j) Hy_x(i, j) Hy_z(i, j);这里的b_x、c_x是预计算好的衰减系数数组写法取自分裂场PML的常见离散格式。逻辑上Hy_x只由Ez沿x方向的差分驱动并带上σ_x衰减Hy_z由y方向差分驱动、不受σ_x影响再把两者相加还原出Hy。参数要点σ_x在同一层内用一个标量即可但必须沿层号渐变不能直接从0跳到σ_maxc_x通常取 (1 - b_x) 的某种比例不必每一层独立推导系数只要保证b_x、c_x由同一σ_x生成。这一实现是后续四种的对照基准它的弱点是角点交界处两个方向的σ同时作用容易写错且长时间仿真中会出现通常所说的late-time reflection也就是低频尾反射。2.2 CPML卷积PML工程上默认选项CPMLConvolutional PML由Roden与Gedney在2000年提出是目前电磁仿真里默认使用的PML形式。它把坐标伸缩从频域算子改写成时域卷积又用指数递归近似卷积结果于是每个方向只需维护一两个辅助变量不分裂场分量修改起来比分裂场干净得多。% CPML系数与辅助变量psi的递归以x方向、Ez的x偏导为例 b_x exp( -(sig_x ./ kappa_x alpha_x) * dt / eps0 ); c_x sig_x ./ (sig_x .* kappa_x kappa_x.^2 .* alpha_x) .* (b_x - 1); psi_ez_x b_x .* psi_ez_x c_x .* (Ez(:, 2:nx) - Ez(:, 1:nx-1)) / dx;逻辑说明第一行是卷积指数的衰减因子里面sig_x是PML电导率kappa_x是坐标拉伸系数alpha_x是频移因子第二行把衰减因子转换成卷积幅度系数c_x第三行把当前差分依次卷积进去。注意b_x、c_x是一维数组按层预计算psi_ez_x是二维辅助数组维度与PML内网格一致。参数要点起始配置建议kappa_x1、alpha_x0先用最简形式跑通再按后续章节的调参方式优化sig_x则与分裂场共用同一多项式剖面保证对比公平性。CPML对斜入射和掠射角的容忍度明显高于分裂场是网格不细时最稳的选择。2.3 UPML用各向异性介质张量吸收边界UPML的出发点是“各向异性介质”视角把PML区看作一层拥有复相对介电常数张量和复磁导率张量的介质沿边界的切向与法向分量使用不同张量分量从而在介质层内部实现阻抗连续。MATLAB里常见的做法是预先计算一维系数数组ca_x、ca_y、cb_x、cb_y在更新电场和磁场时根据当前格点位置选取对应系数。% 左边界PML内Ez更新时先取材料系数数组 Ez(i, j) ca_x(i) .* Ez(i, j) cb_x(i) .* ( (Hy(i, j) - Hy(i, j-1)) / dx ... - (Hx(i1, j) - Hx(i, j)) / dy );这段代码与常规FDTD更新式几乎一样区别只在系数ca_x、cb_x沿PML层数渐变而不是常数1。逻辑说明ca_x对应时间项cb_x对应旋度项二者组合保证了PML介质内的波阻抗与主域一致。参数要点实现时最容易出错的是“方向对应关系”——x方向PML内切向是Ey、Hz法向是Ex而本处2D TE的坐标写得再简化也不能把ca_x、cb_x互相颠倒。UPML对垂直入射的反射性能与CPML相当但推导形式对并行计算和色散介质耦合更友好作为对比组它的行为能帮判断“吸收好坏是否来自卷积近似”。2.4 NPML近PML代码最简的一类NPML由Cummer提出思路是把PML看成对场分量做复坐标映射但在每个格点只需要一组预计算好的插值权重就可以把标准FDTD的导数换成复坐标下的组合导数。它和CPML解决的是同一个问题但实现上更接近“改差分模板”而不是“维护辅助变量”。% NPML: 用预计算的interp_fx、interp_gx替换标准x方向差分 dEz_dx interp_fx(i) .* (Ez(i1, j) - Ez(i, j)) / dx ... interp_gx(i) .* (Ez(i, j) - Ez(i-1, j)) / dx;这里interp_fx、interp_gx是初始化时根据σ_x、σ_y计算出的两层权重二者之和应保持为1否则模板会引入常数偏移。逻辑说明权重随PML层位置变化在靠近主域的一侧与标准差分几乎一致在外边界一侧则加权混合从而在数学上等价于坐标拉伸。参数要点NPML不需要psi数组内存占用比CPML少一截但对掠射角的反射通常比CPML高几个分贝在二维均匀网格中它是最容易改写的PML适合先写出来做正确性验证。权重计算在初始化时完成一次主循环里没有任何指数或除法运算因此单步耗时最低。2.5 指数衰减层不是PML但必须留一个对照组还有一类代码在MATLAB问答社区里流传很广把边界区场值每步乘一个指数衰减因子exp(-sigma*dt)再继续正常更新。它在形式上“有PML三个字母”实际只是外加阻尼遇到波阻抗变化照样会反射。这里保留它不是为了宣传而是让对比图里有一条“非PML”底线用来确认其它四种确实在吸收机制上比单纯衰减高一个量级。% DECAY层在PML区直接衰减场非PML Ez(pml_x, :) Ez(pml_x, :) .* exp(-sigma_x .* dt);逻辑说明这就是一阶衰减加入后低频分量衰减慢、高频分量衰减快频谱响应不平坦与真正的PML相比它缺少阻抗匹配因此在边界处会有一次可见反射。参数要点sigma_x可以沿用同样的剖面但即使剖面对也无法获得PML级别的宽频吸收。把它作为第5种实现放进统一框架正是为了回答“如果不写PML、只加衰减损失到底多大”这个常见问题。2.6 五种实现的结构对比速查表下表从实现结构角度给出一个快速判断具体数值会随源、网格与厚度变化后面章节再给更严谨的测量方法。实现辅助数组数量内存特征实现复杂度典型短板分裂场SF-PML4个分裂场分量额外存4个数组中高角点处理与晚时反射CPML每方向2个psi中等低系数预计算要写对UPML无明显psi主要存系数数组中系数与坐标方向易对应错NPML0个psi最低最低掠射角反射偏高指数衰减层0最低最低阻抗不匹配宽频反射选型结论可以很直白不在色散介质里做研究时CPML是最省心的默认要写教学代码或快速验证算法正确性NPML能最快跑通想在论文里展示不同吸收边界的差异分裂场和指数衰减层作为对照组价值最高。3. 用同一个FDTD框架跑通5种PML3.1 最小可运行的MATLAB主循环骨架五种PML放在一起时最容易出的问题是“变量名一样但含义不同”因此先把统一接口定下来。下面是一个最小但结构完整的MATLAB主循环模式是主域更新、边界更新分离五种PML各自实现同一个函数签名。function [probe, param] run_pml_cmp(pml_type, nx, nt) % pml_type可选 SF CPML UPML NPML DECAY npml 10; dx 1e-3; dt dx / (phys_c * sqrt(2)); % CFL1 [Hx, Hy, Ez] deal(zeros(nx, nx)); probe zeros(nt, 1); px round(nx/2); py round(nx/2); for n 1:nt Ez(px, py) Ez(px, py) sin(2*pi*1e9*n*dt)^2; % 点源 Hx(1:end-2,:) Hx(1:end-2,:) - dt/(mu0*dx) * diff(Ez, 1, 1); Hy(:,1:end-2) Hy(:,1:end-2) dt/(mu0*dx) * diff(Ez, 1, 2); Ez(2:end-1,2:end-1) Ez(2:end-1,2:end-1) ... dt/eps0/dx * (diff(Hy, 1, 2) - diff(Hx, 1, 1)); switch lower(pml_type) case cpml [Ez, Hx, Hy] cpml_step(Ez, Hx, Hy, pml_c, dt, npml); case sf [Ez, Hx, Hy] sf_step(Ez, Hx, Hy, pml_c, dt, npml); case upml [Ez, Hx, Hy] upml_step(Ez, Hx, Hy, pml_c, dt, npml); case npml [Ez, Hx, Hy] npml_step(Ez, Hx, Hy, pml_c, dt, npml); case decay [Ez, Hx, Hy] decay_step(Ez, Hx, Hy, pml_c, dt, npml); end probe(n) Ez(px, py); end代码说明主循环沿用经典2D FDTD的差分顺序即先更新H、再更新EPML步骤在每次E更新后执行。pml_c是调用方预先构造的系数结构体包含sig_x、sig_y、kappa_x、kappa_y、alpha_x、alpha_y五种实现都从同一剖面给出这样单步耗时和反射误差的差异才能全部归因到算法本身而不是参数不一致。注意点源放在中心格点且只在初期激活避免持续激励把反射波淹没在直达波里探针记录点要在物理区内而非PML内部。启动一个仿真的命令是matlab -batch run_pml_cmp(CPML, 120, 1000)matlab -batch在R2019a及之后都可用适合无GUI的批处理对比。这里要求每种PML的边界函数共用同一套pml_c结构各函数内只用索引区分“下边界、上边界、左右边界和角点”这是整个对比脚本中最容易返工的地方。3.2 五种PML必须共享同一份参数剖面也许有人会在对比时为了“让每种PML表现好一点”而单独调参这会直接破坏对比公平性。正确的做法是先确定源类型、网格步长、PML厚度和σ剖面让五种实现全部使用同一组预计算参数对实现本身预留的额外自由度如CPML的kappa与alpha、UPML的系数表保持默认一致除非单独标注“这是最优参数下的性能”。这样才能画出有意义的对比曲线。3.3 从主循环里看各实现的工作量差异从上述框架可以直观看到五种PML在单步内的额外操作分别是分裂场需要对Hx、Hy各自做两次带系数递归CPML需要更新psi并对Ez加上卷积修正UPML需要把标准更新后的场再乘一次材料系数并更新边界区NPML只需用权重模板替换差分DECAY层只做一次乘法。单从MATLAB向量化角度看NPML与DECAY的边界代码最短CPML其次分裂场和UPML最长。我一般会根据步数规模来选择单次仿真在万步以内优先用CPML并接受它的psi数组要扫几千组参数时换NPML把单位步长开销压缩下来能省出几小时等待。4. 对比实验把吸收好坏变成可写进报告的数字4.1 用大小域差值算归一化反射误差吸收边界性能的标准测量方法是“同源同探针双仿真”第一次用足够大的计算域并加PML认为探针处不会有边界反射参与作为参考第二次把域缩小到目标尺寸且同样加PML探针位置不变两次结果的差就是边界反射波。归一化表达式写成峰值误差和时域平均误差两种更合理。% 假设ref_probe来自大域small_probe来自目标域二者长度一致 err small_probe - ref_probe; norm_peak 20 * log10(max(abs(err)) / max(abs(ref_probe))); norm_avg 20 * log10(sqrt(sum(err.^2)) / sqrt(sum(ref_probe.^2)));代码逻辑很清楚峰值误差反映单次反射的最大瞬时贡献适合观察瞬时脉冲平均误差反映长时间仿真的整体污染水平适合判断晚时反射。参数说明大域的边长至少比小域多出2倍PML厚度再加20格缓冲避免大域自身PML反射在观察时段内先到达探针探针离边界距离要固定否则不同PML对近场的感应差异会混入结果。量化时建议把误差信号在时域持续累积到源完全熄灭后再延长至少1000步以便让晚时反射充分进入统计。4.2 一个可复现的典型量级参考表下面的量级取自均匀网格2D TE、点源、PML厚度10层、源频带与网格满足dx不超过lambda/10时的常见结果具体项目里的绝对值会有偏移但彼此之间的相对次序一般是稳定的。实现峰值反射误差dB晚时平均误差dB额外内存占用相对主域单步耗时相对DECAY-18 ~ -28-20 ~ -3001.00分裂场SF-PML-55 ~ -75-55 ~ -65约50%1.35UPML-55 ~ -70-55 ~ -72约20%1.20NPML-50 ~ -65-50 ~ -68约5%1.05CPML-60 ~ -85-60 ~ -85约25%1.15解释下这张表怎么读DECAY的峰值在-20dB附近说明单纯衰减会造成显著边界回波SF-PML如果长时间运行峰值和平均误差之间的差距会变大NPML的差距相对小CPML在垂直入射且σ剖面给当时能稳定压到-70dB以下。参数说明如果换成40度斜入射DECAY基本不变NPML和分裂场可能各上升6到10dBCPML通常仍低于-50dB所以表格只适用于默认垂直入射。4.3 三个可以直接落地的规律第一峰值反射误差主要由靠近主域的PML区域决定因此σ剖面在前两层里必须缓慢上升而不是把σ_max直接顶到第一个格点。第二晚时误差主要由低频分量决定CPML导入alpha_x之后能显著压低分裂场即使增加厚度也很难改善这是机制上的差异。第三PML内部有什么物理场不重要探针和结果提取区域必须留在物理域内否则会把吸收率误判成反射误差。这三点在后续章节是反复出现的检查项。5. 参数该怎么设厚度、σ剖面、κ和α5.1 一组可靠起点参数PML参数不需要每次从零开始猜。常见做法是从一套经验起点出发先跑通、再根据频谱和误差曲线微调。下面的起点适用于均匀网格、无耗背景介质。参数起点值建议范围作用与风险npml108~20太薄吸收不足太厚增加内存且对晚时误差改善有限m32~4剖面多项式阶数过大时相邻层变化过缓但起始段效果差sigma_max0.8 * (m1) / (eta0 * dx)0.5 ~ 1.2倍经验式主导吸收强度过大会在边界处形成阻抗跳变kappa_max11~20拉长坐标伸缩主要帮助掠射角alpha_max00~0.05/(dt*eps0)左右抑制低频晚时反射过大会破坏匹配这个经验式的直观含义真空中eta0约377欧姆介质中要换成介质波阻抗也就是eta0除以sqrt(epsilon_r)磁介质按对应修正这样sigma_max会随介质波阻抗自动缩放。参数要点厚度不是线性决定性能的从10层加到16层收益可观从16层加到22层收益会明显变缓看曲线时可以早点停止。5.2 σ、κ、α剖面生成代码与绘图检查参数剖面是PML实现的入口也是所有代码中最值得先画出来看的部分。下面是生成一维剖面并立即用MATLAB画图的代码% 生成PML剖面npml为层数m为多项式阶数 profile (0.5 : npml-0.5) / npml; % 每层的中心位置 sigma_x sigma_max * profile .^ m; kappa_x 1 (kappa_max - 1) * profile .^ m; alpha_x alpha_max * (1 - profile) .^ m; figure; subplot(3,1,1); plot(1:npml, sigma_x); title(sigma); subplot(3,1,2); plot(1:npml, kappa_x); title(kappa); subplot(3,1,3); plot(1:npml, alpha_x); title(alpha);逻辑说明profile用每层中心位置而不是边界位置避免在PML起始处出现半格的零厚度sigma_x从最内层向外递增kappa_x同向递增alpha_x则反向递减这样在物理域与PML交界处alpha最大、sigma最小交界阻抗连续。参数说明如果画出的sigma曲线在最外层还继续增大而不饱和说明npml不够或m偏大此时反射来自截断边界而不是PML内部如果sigma前两格增长过快则主域边界会有明显的一阶反射峰值。5.3 kappa和alpha的参数踩坑这两个参数容易调反。kappa_x增大会让坐标拉伸更强但同时也放大该层内差分误差kappa_max从10往上调对垂直入射几乎无感对45度以上掠射角才有效因此不要一上来就设成大数。alpha_x的目的是把坐标伸缩的极点从0频率挪到非零频率减少直流和低频分量的晚时反射但它会削弱极低频的吸收能力——alpha_max取到经验式上限附近时近DC分量可能反而反射回来。调参顺序推荐先固定alpha0把npml和m调到垂直入射峰值达到-70dB加入掠射角测试后需要提一档kappa最后再看长时间误差频谱决定是否注入alpha每一步都复用4.1的测量脚本而不是同时调整多个参数。5.4 网格离散度决定PML性能上限PML性能不是孤立指标。网格越粗差分对斜入射和近掠射角的数值相速误差越大PML的匹配也就越不准。经验上至少保持dx不超过最短目标波长的十分之一做掠射角对比时建议到二十分之一。时间步长按CFL条件取2D均匀网格约dt dx除以c乘以sqrt(2)不满足时PML层内部的指数系数会出现虚部误差表现为仿真后期边界区域数值发散。判断网格是否够细的笨办法把同一PML配置放到两倍细网格上如果峰值反射改善量小于3dB说明当前网格已经不是主要瓶颈。6. 验证吸收质量与三个高频坑6.1 用频谱而不是只盯时域判PML好坏时域峰值误差只能反映“某一瞬间最大反射”低频晚时反射往往在波形上只是缓慢漂移肉眼很难察觉。更好的验证方式是做一次FFT看反射误差的频谱是否在全频带上都低于目标值。% 探针误差信号err做FFT归一化到峰值并换算成dB频谱 Nf numel(err); f (0:Nf-1)/Nf/dt; spec abs(fft(err)); spec_dB 20*log10(spec / max(spec(2:end))); % 忽略直流分量 semilogx(f(2:Nf/2), spec_dB(2:Nf/2)); ylim([-100 0]);代码逻辑err来自4.1的大域与小域差值它本身就是反射波序列FFT后能看出反射能量集中在低频还是高频从而判断是晚时反射还是离散误差用log横轴更适合观察低频端。参数说明采样时长要覆盖源熄灭后的多轮反射否则FFT频率分辨率不够低频段会出现虚假起伏。验证标准可以这样定目标频带内反射谱低于-60dB的PML配置用于定量仿真才算合格-40dB以上的配置只能用于波形演示。6.2 坑一sigma剖面第一格就顶满许多初版代码会把sigma_max直接赋给PML每一层导致主域与PML交界处出现阻抗阶跃波还没进入PML就先反射一部分。判断方法是看5.2那三张曲线图sigma应从0缓慢爬升不是从最高值起始。修正方法是把sigma_x sigma_max * profile.^m里的profile从0到1渐变并保证最内层sigma相对sigma_max低于1%量级。6.3 坑二角点区域沿用一个方向公式边界段与角点段必须分开处理。边段只需要一个方向的sigma和一个方向的kappa角点段位于两个方向的PML交叠里必须同时使用两套方向系数。常见做法是用mask矩阵分别标记上、下、左、右与四个角点各自套不同更新式如果只写边段不写角点角点会成为二次辐射源表现为延迟一段时间后从角落返回一圈明显的弧状反射。检查方法是把误差场画成空间云图看是否有一圈圆弧从角点方向开始传播。6.4 坑三把探针放进PML里读数想“直接看到吸收效果”而把探针放在PML区读出来的曲线更像衰减包络不是反射误差会把PML的缺点掩盖。探针必须位于物理域内离PML内边界至少留3到5格物理域余量越小近场差异越容易被误判为吸收差。最后给一个固定调参路径用CPML先跑垂直入射固定npml10、m3、alpha_max0把sigma_max按5.1公式初始化确认峰值反射低于-70dB后加45度斜入射复核若恶化超过5dB再逐格上调kappa_max直到频谱反射在目标频带稳定低于-60dB。本文还有配套的精品资源点击获取 MATLAB FDTD仿真中五种PML实现对比与参数调优
简介这组MATLAB代码聚焦五种完美匹配层PML实现与效果对比适用于FDTD等数值仿真中需要抑制边界反射的电磁波、声波或弹性波问题也适合刚接触吸收边界条件的研究者或高年级本科生对照学习。压缩包共23个文件包含12个.m脚本、10个.mat数据文件和1份PDF文档体积仅301KB代码中分别实现了Smith-PML、Steklov-Poincaré PML、多层PML、自适应PML与改进的Bermúdez PML等不同思路并配有扫参脚本和反射计算函数可直观比较反射系数、衰减效果与计算效率。目前已有631人学习包内附有PML Comparison.pdf对比说明和对应的TMz仿真数据便于快速复现五类PML在12GHz中心频率、24GHz带宽下的表现可为实际建模中选取合适PML提供直接参考。1. 为什么要在MATLAB里对比5种PML做二维FDTD电磁波仿真的工程师多半已经学会了“在区域四周加PML”来压缩反射既然PML在理论上接近完美吸收为什么还要对比因为“接近完美”只在离散网格、垂直入射、有限带宽这三个条件下成立。换一种PML实现改一行系数可能在斜入射、低频或长时间步进时反射误差相差两个数量级而这类误差会直接进入后续的近场、远场提取结果。这里用同一个2D TE波FDTD框架把Berenger分裂场PML、UPML、CPML、NPML以及一个经常被误认成PML的指数衰减层同时实现给出可运行的MATLAB代码、参数设置和量化对比方法。这组内容不是要告诉你“哪个最好”而是给出一个能自己复现、能在论文里画对比曲线的基准。2. 五种PML的原理选型和MATLAB关键递归2.1 Berenger分裂场PML作为对比基准的经典实现1994年Berenger提出的分裂场PML核心思路是把边界区域内的场分量按坐标方向劈裂让波在不同方向拥有不同的虚衰减率从而在没有真实反射界面的情况下把外向波“吞掉”。在二维TE波Ez、Hx、Hy中常见做法是对Hx和Hy各保存两个子分量Ez保持单一场靠不同方向上的电导率σ_x、σ_y驱动衰减。% 以x方向PML区域内的Hy为例Hy拆成Hy_x与Hy_z Hy_x(i, j) b_x(i) .* Hy_x(i, j) - c_x(i) .* (Ez(i1, j) - Ez(i, j)) / dx; Hy_z(i, j) Hy_z(i, j) (Ez(i, j1) - Ez(i, j)) / dy; Hy(i, j) Hy_x(i, j) Hy_z(i, j);这里的b_x、c_x是预计算好的衰减系数数组写法取自分裂场PML的常见离散格式。逻辑上Hy_x只由Ez沿x方向的差分驱动并带上σ_x衰减Hy_z由y方向差分驱动、不受σ_x影响再把两者相加还原出Hy。参数要点σ_x在同一层内用一个标量即可但必须沿层号渐变不能直接从0跳到σ_maxc_x通常取 (1 - b_x) 的某种比例不必每一层独立推导系数只要保证b_x、c_x由同一σ_x生成。这一实现是后续四种的对照基准它的弱点是角点交界处两个方向的σ同时作用容易写错且长时间仿真中会出现通常所说的late-time reflection也就是低频尾反射。2.2 CPML卷积PML工程上默认选项CPMLConvolutional PML由Roden与Gedney在2000年提出是目前电磁仿真里默认使用的PML形式。它把坐标伸缩从频域算子改写成时域卷积又用指数递归近似卷积结果于是每个方向只需维护一两个辅助变量不分裂场分量修改起来比分裂场干净得多。% CPML系数与辅助变量psi的递归以x方向、Ez的x偏导为例 b_x exp( -(sig_x ./ kappa_x alpha_x) * dt / eps0 ); c_x sig_x ./ (sig_x .* kappa_x kappa_x.^2 .* alpha_x) .* (b_x - 1); psi_ez_x b_x .* psi_ez_x c_x .* (Ez(:, 2:nx) - Ez(:, 1:nx-1)) / dx;逻辑说明第一行是卷积指数的衰减因子里面sig_x是PML电导率kappa_x是坐标拉伸系数alpha_x是频移因子第二行把衰减因子转换成卷积幅度系数c_x第三行把当前差分依次卷积进去。注意b_x、c_x是一维数组按层预计算psi_ez_x是二维辅助数组维度与PML内网格一致。参数要点起始配置建议kappa_x1、alpha_x0先用最简形式跑通再按后续章节的调参方式优化sig_x则与分裂场共用同一多项式剖面保证对比公平性。CPML对斜入射和掠射角的容忍度明显高于分裂场是网格不细时最稳的选择。2.3 UPML用各向异性介质张量吸收边界UPML的出发点是“各向异性介质”视角把PML区看作一层拥有复相对介电常数张量和复磁导率张量的介质沿边界的切向与法向分量使用不同张量分量从而在介质层内部实现阻抗连续。MATLAB里常见的做法是预先计算一维系数数组ca_x、ca_y、cb_x、cb_y在更新电场和磁场时根据当前格点位置选取对应系数。% 左边界PML内Ez更新时先取材料系数数组 Ez(i, j) ca_x(i) .* Ez(i, j) cb_x(i) .* ( (Hy(i, j) - Hy(i, j-1)) / dx ... - (Hx(i1, j) - Hx(i, j)) / dy );这段代码与常规FDTD更新式几乎一样区别只在系数ca_x、cb_x沿PML层数渐变而不是常数1。逻辑说明ca_x对应时间项cb_x对应旋度项二者组合保证了PML介质内的波阻抗与主域一致。参数要点实现时最容易出错的是“方向对应关系”——x方向PML内切向是Ey、Hz法向是Ex而本处2D TE的坐标写得再简化也不能把ca_x、cb_x互相颠倒。UPML对垂直入射的反射性能与CPML相当但推导形式对并行计算和色散介质耦合更友好作为对比组它的行为能帮判断“吸收好坏是否来自卷积近似”。2.4 NPML近PML代码最简的一类NPML由Cummer提出思路是把PML看成对场分量做复坐标映射但在每个格点只需要一组预计算好的插值权重就可以把标准FDTD的导数换成复坐标下的组合导数。它和CPML解决的是同一个问题但实现上更接近“改差分模板”而不是“维护辅助变量”。% NPML: 用预计算的interp_fx、interp_gx替换标准x方向差分 dEz_dx interp_fx(i) .* (Ez(i1, j) - Ez(i, j)) / dx ... interp_gx(i) .* (Ez(i, j) - Ez(i-1, j)) / dx;这里interp_fx、interp_gx是初始化时根据σ_x、σ_y计算出的两层权重二者之和应保持为1否则模板会引入常数偏移。逻辑说明权重随PML层位置变化在靠近主域的一侧与标准差分几乎一致在外边界一侧则加权混合从而在数学上等价于坐标拉伸。参数要点NPML不需要psi数组内存占用比CPML少一截但对掠射角的反射通常比CPML高几个分贝在二维均匀网格中它是最容易改写的PML适合先写出来做正确性验证。权重计算在初始化时完成一次主循环里没有任何指数或除法运算因此单步耗时最低。2.5 指数衰减层不是PML但必须留一个对照组还有一类代码在MATLAB问答社区里流传很广把边界区场值每步乘一个指数衰减因子exp(-sigma*dt)再继续正常更新。它在形式上“有PML三个字母”实际只是外加阻尼遇到波阻抗变化照样会反射。这里保留它不是为了宣传而是让对比图里有一条“非PML”底线用来确认其它四种确实在吸收机制上比单纯衰减高一个量级。% DECAY层在PML区直接衰减场非PML Ez(pml_x, :) Ez(pml_x, :) .* exp(-sigma_x .* dt);逻辑说明这就是一阶衰减加入后低频分量衰减慢、高频分量衰减快频谱响应不平坦与真正的PML相比它缺少阻抗匹配因此在边界处会有一次可见反射。参数要点sigma_x可以沿用同样的剖面但即使剖面对也无法获得PML级别的宽频吸收。把它作为第5种实现放进统一框架正是为了回答“如果不写PML、只加衰减损失到底多大”这个常见问题。2.6 五种实现的结构对比速查表下表从实现结构角度给出一个快速判断具体数值会随源、网格与厚度变化后面章节再给更严谨的测量方法。实现辅助数组数量内存特征实现复杂度典型短板分裂场SF-PML4个分裂场分量额外存4个数组中高角点处理与晚时反射CPML每方向2个psi中等低系数预计算要写对UPML无明显psi主要存系数数组中系数与坐标方向易对应错NPML0个psi最低最低掠射角反射偏高指数衰减层0最低最低阻抗不匹配宽频反射选型结论可以很直白不在色散介质里做研究时CPML是最省心的默认要写教学代码或快速验证算法正确性NPML能最快跑通想在论文里展示不同吸收边界的差异分裂场和指数衰减层作为对照组价值最高。3. 用同一个FDTD框架跑通5种PML3.1 最小可运行的MATLAB主循环骨架五种PML放在一起时最容易出的问题是“变量名一样但含义不同”因此先把统一接口定下来。下面是一个最小但结构完整的MATLAB主循环模式是主域更新、边界更新分离五种PML各自实现同一个函数签名。function [probe, param] run_pml_cmp(pml_type, nx, nt) % pml_type可选 SF CPML UPML NPML DECAY npml 10; dx 1e-3; dt dx / (phys_c * sqrt(2)); % CFL1 [Hx, Hy, Ez] deal(zeros(nx, nx)); probe zeros(nt, 1); px round(nx/2); py round(nx/2); for n 1:nt Ez(px, py) Ez(px, py) sin(2*pi*1e9*n*dt)^2; % 点源 Hx(1:end-2,:) Hx(1:end-2,:) - dt/(mu0*dx) * diff(Ez, 1, 1); Hy(:,1:end-2) Hy(:,1:end-2) dt/(mu0*dx) * diff(Ez, 1, 2); Ez(2:end-1,2:end-1) Ez(2:end-1,2:end-1) ... dt/eps0/dx * (diff(Hy, 1, 2) - diff(Hx, 1, 1)); switch lower(pml_type) case cpml [Ez, Hx, Hy] cpml_step(Ez, Hx, Hy, pml_c, dt, npml); case sf [Ez, Hx, Hy] sf_step(Ez, Hx, Hy, pml_c, dt, npml); case upml [Ez, Hx, Hy] upml_step(Ez, Hx, Hy, pml_c, dt, npml); case npml [Ez, Hx, Hy] npml_step(Ez, Hx, Hy, pml_c, dt, npml); case decay [Ez, Hx, Hy] decay_step(Ez, Hx, Hy, pml_c, dt, npml); end probe(n) Ez(px, py); end代码说明主循环沿用经典2D FDTD的差分顺序即先更新H、再更新EPML步骤在每次E更新后执行。pml_c是调用方预先构造的系数结构体包含sig_x、sig_y、kappa_x、kappa_y、alpha_x、alpha_y五种实现都从同一剖面给出这样单步耗时和反射误差的差异才能全部归因到算法本身而不是参数不一致。注意点源放在中心格点且只在初期激活避免持续激励把反射波淹没在直达波里探针记录点要在物理区内而非PML内部。启动一个仿真的命令是matlab -batch run_pml_cmp(CPML, 120, 1000)matlab -batch在R2019a及之后都可用适合无GUI的批处理对比。这里要求每种PML的边界函数共用同一套pml_c结构各函数内只用索引区分“下边界、上边界、左右边界和角点”这是整个对比脚本中最容易返工的地方。3.2 五种PML必须共享同一份参数剖面也许有人会在对比时为了“让每种PML表现好一点”而单独调参这会直接破坏对比公平性。正确的做法是先确定源类型、网格步长、PML厚度和σ剖面让五种实现全部使用同一组预计算参数对实现本身预留的额外自由度如CPML的kappa与alpha、UPML的系数表保持默认一致除非单独标注“这是最优参数下的性能”。这样才能画出有意义的对比曲线。3.3 从主循环里看各实现的工作量差异从上述框架可以直观看到五种PML在单步内的额外操作分别是分裂场需要对Hx、Hy各自做两次带系数递归CPML需要更新psi并对Ez加上卷积修正UPML需要把标准更新后的场再乘一次材料系数并更新边界区NPML只需用权重模板替换差分DECAY层只做一次乘法。单从MATLAB向量化角度看NPML与DECAY的边界代码最短CPML其次分裂场和UPML最长。我一般会根据步数规模来选择单次仿真在万步以内优先用CPML并接受它的psi数组要扫几千组参数时换NPML把单位步长开销压缩下来能省出几小时等待。4. 对比实验把吸收好坏变成可写进报告的数字4.1 用大小域差值算归一化反射误差吸收边界性能的标准测量方法是“同源同探针双仿真”第一次用足够大的计算域并加PML认为探针处不会有边界反射参与作为参考第二次把域缩小到目标尺寸且同样加PML探针位置不变两次结果的差就是边界反射波。归一化表达式写成峰值误差和时域平均误差两种更合理。% 假设ref_probe来自大域small_probe来自目标域二者长度一致 err small_probe - ref_probe; norm_peak 20 * log10(max(abs(err)) / max(abs(ref_probe))); norm_avg 20 * log10(sqrt(sum(err.^2)) / sqrt(sum(ref_probe.^2)));代码逻辑很清楚峰值误差反映单次反射的最大瞬时贡献适合观察瞬时脉冲平均误差反映长时间仿真的整体污染水平适合判断晚时反射。参数说明大域的边长至少比小域多出2倍PML厚度再加20格缓冲避免大域自身PML反射在观察时段内先到达探针探针离边界距离要固定否则不同PML对近场的感应差异会混入结果。量化时建议把误差信号在时域持续累积到源完全熄灭后再延长至少1000步以便让晚时反射充分进入统计。4.2 一个可复现的典型量级参考表下面的量级取自均匀网格2D TE、点源、PML厚度10层、源频带与网格满足dx不超过lambda/10时的常见结果具体项目里的绝对值会有偏移但彼此之间的相对次序一般是稳定的。实现峰值反射误差dB晚时平均误差dB额外内存占用相对主域单步耗时相对DECAY-18 ~ -28-20 ~ -3001.00分裂场SF-PML-55 ~ -75-55 ~ -65约50%1.35UPML-55 ~ -70-55 ~ -72约20%1.20NPML-50 ~ -65-50 ~ -68约5%1.05CPML-60 ~ -85-60 ~ -85约25%1.15解释下这张表怎么读DECAY的峰值在-20dB附近说明单纯衰减会造成显著边界回波SF-PML如果长时间运行峰值和平均误差之间的差距会变大NPML的差距相对小CPML在垂直入射且σ剖面给当时能稳定压到-70dB以下。参数说明如果换成40度斜入射DECAY基本不变NPML和分裂场可能各上升6到10dBCPML通常仍低于-50dB所以表格只适用于默认垂直入射。4.3 三个可以直接落地的规律第一峰值反射误差主要由靠近主域的PML区域决定因此σ剖面在前两层里必须缓慢上升而不是把σ_max直接顶到第一个格点。第二晚时误差主要由低频分量决定CPML导入alpha_x之后能显著压低分裂场即使增加厚度也很难改善这是机制上的差异。第三PML内部有什么物理场不重要探针和结果提取区域必须留在物理域内否则会把吸收率误判成反射误差。这三点在后续章节是反复出现的检查项。5. 参数该怎么设厚度、σ剖面、κ和α5.1 一组可靠起点参数PML参数不需要每次从零开始猜。常见做法是从一套经验起点出发先跑通、再根据频谱和误差曲线微调。下面的起点适用于均匀网格、无耗背景介质。参数起点值建议范围作用与风险npml108~20太薄吸收不足太厚增加内存且对晚时误差改善有限m32~4剖面多项式阶数过大时相邻层变化过缓但起始段效果差sigma_max0.8 * (m1) / (eta0 * dx)0.5 ~ 1.2倍经验式主导吸收强度过大会在边界处形成阻抗跳变kappa_max11~20拉长坐标伸缩主要帮助掠射角alpha_max00~0.05/(dt*eps0)左右抑制低频晚时反射过大会破坏匹配这个经验式的直观含义真空中eta0约377欧姆介质中要换成介质波阻抗也就是eta0除以sqrt(epsilon_r)磁介质按对应修正这样sigma_max会随介质波阻抗自动缩放。参数要点厚度不是线性决定性能的从10层加到16层收益可观从16层加到22层收益会明显变缓看曲线时可以早点停止。5.2 σ、κ、α剖面生成代码与绘图检查参数剖面是PML实现的入口也是所有代码中最值得先画出来看的部分。下面是生成一维剖面并立即用MATLAB画图的代码% 生成PML剖面npml为层数m为多项式阶数 profile (0.5 : npml-0.5) / npml; % 每层的中心位置 sigma_x sigma_max * profile .^ m; kappa_x 1 (kappa_max - 1) * profile .^ m; alpha_x alpha_max * (1 - profile) .^ m; figure; subplot(3,1,1); plot(1:npml, sigma_x); title(sigma); subplot(3,1,2); plot(1:npml, kappa_x); title(kappa); subplot(3,1,3); plot(1:npml, alpha_x); title(alpha);逻辑说明profile用每层中心位置而不是边界位置避免在PML起始处出现半格的零厚度sigma_x从最内层向外递增kappa_x同向递增alpha_x则反向递减这样在物理域与PML交界处alpha最大、sigma最小交界阻抗连续。参数说明如果画出的sigma曲线在最外层还继续增大而不饱和说明npml不够或m偏大此时反射来自截断边界而不是PML内部如果sigma前两格增长过快则主域边界会有明显的一阶反射峰值。5.3 kappa和alpha的参数踩坑这两个参数容易调反。kappa_x增大会让坐标拉伸更强但同时也放大该层内差分误差kappa_max从10往上调对垂直入射几乎无感对45度以上掠射角才有效因此不要一上来就设成大数。alpha_x的目的是把坐标伸缩的极点从0频率挪到非零频率减少直流和低频分量的晚时反射但它会削弱极低频的吸收能力——alpha_max取到经验式上限附近时近DC分量可能反而反射回来。调参顺序推荐先固定alpha0把npml和m调到垂直入射峰值达到-70dB加入掠射角测试后需要提一档kappa最后再看长时间误差频谱决定是否注入alpha每一步都复用4.1的测量脚本而不是同时调整多个参数。5.4 网格离散度决定PML性能上限PML性能不是孤立指标。网格越粗差分对斜入射和近掠射角的数值相速误差越大PML的匹配也就越不准。经验上至少保持dx不超过最短目标波长的十分之一做掠射角对比时建议到二十分之一。时间步长按CFL条件取2D均匀网格约dt dx除以c乘以sqrt(2)不满足时PML层内部的指数系数会出现虚部误差表现为仿真后期边界区域数值发散。判断网格是否够细的笨办法把同一PML配置放到两倍细网格上如果峰值反射改善量小于3dB说明当前网格已经不是主要瓶颈。6. 验证吸收质量与三个高频坑6.1 用频谱而不是只盯时域判PML好坏时域峰值误差只能反映“某一瞬间最大反射”低频晚时反射往往在波形上只是缓慢漂移肉眼很难察觉。更好的验证方式是做一次FFT看反射误差的频谱是否在全频带上都低于目标值。% 探针误差信号err做FFT归一化到峰值并换算成dB频谱 Nf numel(err); f (0:Nf-1)/Nf/dt; spec abs(fft(err)); spec_dB 20*log10(spec / max(spec(2:end))); % 忽略直流分量 semilogx(f(2:Nf/2), spec_dB(2:Nf/2)); ylim([-100 0]);代码逻辑err来自4.1的大域与小域差值它本身就是反射波序列FFT后能看出反射能量集中在低频还是高频从而判断是晚时反射还是离散误差用log横轴更适合观察低频端。参数说明采样时长要覆盖源熄灭后的多轮反射否则FFT频率分辨率不够低频段会出现虚假起伏。验证标准可以这样定目标频带内反射谱低于-60dB的PML配置用于定量仿真才算合格-40dB以上的配置只能用于波形演示。6.2 坑一sigma剖面第一格就顶满许多初版代码会把sigma_max直接赋给PML每一层导致主域与PML交界处出现阻抗阶跃波还没进入PML就先反射一部分。判断方法是看5.2那三张曲线图sigma应从0缓慢爬升不是从最高值起始。修正方法是把sigma_x sigma_max * profile.^m里的profile从0到1渐变并保证最内层sigma相对sigma_max低于1%量级。6.3 坑二角点区域沿用一个方向公式边界段与角点段必须分开处理。边段只需要一个方向的sigma和一个方向的kappa角点段位于两个方向的PML交叠里必须同时使用两套方向系数。常见做法是用mask矩阵分别标记上、下、左、右与四个角点各自套不同更新式如果只写边段不写角点角点会成为二次辐射源表现为延迟一段时间后从角落返回一圈明显的弧状反射。检查方法是把误差场画成空间云图看是否有一圈圆弧从角点方向开始传播。6.4 坑三把探针放进PML里读数想“直接看到吸收效果”而把探针放在PML区读出来的曲线更像衰减包络不是反射误差会把PML的缺点掩盖。探针必须位于物理域内离PML内边界至少留3到5格物理域余量越小近场差异越容易被误判为吸收差。最后给一个固定调参路径用CPML先跑垂直入射固定npml10、m3、alpha_max0把sigma_max按5.1公式初始化确认峰值反射低于-70dB后加45度斜入射复核若恶化超过5dB再逐格上调kappa_max直到频谱反射在目标频带稳定低于-60dB。本文还有配套的精品资源点击获取