简介这是一份基于MATLAB实现的二维FDTD时域有限差分源代码核心针对平面波斜入射场景并在边界处引入PML完美匹配层以模拟开放空间适合电磁场数值计算学习者、科研人员和工程师快速搭建二维散射与传播仿真实验。代码将时间与空间离散化通过迭代更新电场和磁场完整展示了Maxwell方程在二维网格上的数值求解流程同时也体现了斜入射激励源的设置、PML吸收边界构建以及Mur型边界处理等关键技术。资源压缩包共1个文件文件类型为m脚本整体大小仅3KB小巧精炼便于直接运行、逐行修改和扩展。目前已有320人学习下载足见其在二维FDTD教学与入门实践中的参考价值。获取后可重点研读主程序中的网格配置、源信号加入和PML边界更新三大部分对理解FDTD算法原理及实现细节有直接帮助。1. 从文件名到模块拆解二维FDTD斜入射代码到底先看哪里单看Tot_Src_ball_Mur.zip_FDTD 二维 matlab_Oblique incidence_PML_二维FDTD这一长串文件名真正值得先读的并不是二维 FDTD 主更新循环——主循环在任何教材里都能抄到——而是Tot、Src、ball、Mur这四个缩写背后对应的模块。Tot指 Total field / Scattered field总场/散射场分离Src是入射源ball在二维里对应圆柱散射体的圆截面Mur则是一种如今很少单独用于斜入射仿真的吸收边界。斜入射平面波一旦进入三维空间的全矢量问题网格色散、TF/SF 修正项符号、PML 各向异性都会同时发作。反直觉的结论是这类包里保留一阶 Mur 边界不是作者偷懒而是把 TF/SF 斜入射写对之后拿 Mur 的低吸收率当探针——如果连 Mur 的反射曲线都异常平直多半是连接条件符号反了而不是边界本身出了问题。2. 二维FDTD建斜入射平面波TMz极化、CFL与TF/SF连接条件二维 FDTD 没有真正的“球”ball在这个场景里是介质圆柱或 PEC 圆柱的截面。做这类仿真极化选择直接决定代码长度TMz 极化下电场只有Ez一个标量分量磁场有Hx、Hy两个分量存储和更新都比 TEz 少一套斜入射时也不需要额外处理 z 方向的波矢分量。以下内容按 TMz 展开。2.1 TMz极化下Yee网格的场摆放与更新方程Yee 网格里Ez(i,j)定义在网格点(i,j)上Hx(i,j)定义在(i, j1/2)Hy(i,j)定义在(i1/2, j)。这种半格错位让旋度方程的差分天然中心对称不需要插值。无耗介质中的更新方程写成% TMz 二维 FDTD 内部更新dx 与 dy 可以不等 Hx(i,j) Hx(i,j) - (dt / mu0 / dy) * (Ez(i,j1) - Ez(i,j)); Hy(i,j) Hy(i,j) (dt / mu0 / dx) * (Ez(i1,j) - Ez(i,j)); Ez(i,j) Ez(i,j) (dt / eps0 / dx) * (Hy(i,j) - Hy(i-1,j)) ... - (dt / eps0 / dy) * (Hx(i,j) - Hx(i,j-1));Hx的差分发生在 y 方向Hy的差分发生在 x 方向Ez的更新则由两个旋度项共同组成。若dx dy deltaCFL 条件收紧为dt delta / (c * sqrt(2))实际取 0.99 倍以保证数值稳定。空间步长常见取delta lambda / 20斜入射角度越大波前在网格对角线方向上的投影越长需要把步长加密到lambda / 40或更细否则沿边界的相位误差会在远场方向图里表现为非物理的旁瓣。表格里列出 TMz 下各分量的相对位置和更新方向写代码前先对齐这张表能省掉大量下标调试时间分量位置约定更新依赖边界修正对象Ez(i,j)网格节点相邻Hx、Hy吸收边界直接改该值Hx(i,j)y 方向半格Ez(i,j1)-Ez(i,j)TF/SF 左右边界Hy(i,j)x 方向半格Ez(i1,j)-Ez(i,j)TF/SF 上下边界2.2 斜入射波的波矢量拆分与Matlab源项注入斜入射平面波与正入射的本质区别在于波矢量k不再垂直于 TF/SF 矩形边界。设入射角theta相对 y 轴正向传播方向在 x-y 平面内为 30 度则波矢量在两个坐标轴上的投影为theta 30 * pi / 180; kx 2 * pi * f0 / c0 * sin(theta); ky 2 * pi * f0 / c0 * cos(theta); % 入射方向设为 y坐标轴不同时符号要整体翻转 % 硬源注入只在源点叠加不污染整个计算域 Ez(nx0, ny0) Ez(nx0, ny0) amp * sin(2*pi*f0*n*dt - kx*x(nx0) - ky*y(ny0));f0是载频amp是幅度x、y是网格坐标数组。源点要放在 TF/SF 边界内侧不远处而不是整个计算域逐点赋值——逐点赋解析解会导致网格色散与 FDTD 本征色散叠加波形传播几百步后高频分量明显畸变。斜入射下尤其不能只注入Ez而忽略Hx、Hy的入射分量否则总场区里的等效阻抗不匹配会在源点附近激发二次辐射。2.3 TF/SF边界在斜入射角度下的修正式符号与下标怎么对齐TF/SF 的思想是计算域内层存总场外层存散射场两者之间用一层虚拟边界连接。斜入射时最容易写错的是连接条件的方向。以左边界x ix0为例Hy(ix0-1, iy0:iy1)位于边界外侧它更新时本应使用散射区的Ez但数组里那个位置可能存的是总场因此要减掉入射场% 左边界 xix0Hy 的差分跨过 TF/SF 面减掉入射 Ez Hy(ix0-1, iy0:iy1) Hy(ix0-1, iy0:iy1) ... - dt / (mu0 * dx) * Ez_inc(ix0, iy0:iy1); % 右边界 xix1符号相反加上入射 Ez Hy(ix1, iy0:iy1) Hy(ix1, iy0:iy1) ... dt / (mu0 * dx) * Ez_inc(ix1, iy0:iy1);这里Ez_inc是解析入射场在对应网格上的采样可以预生成一个与Ez同尺寸的数组也可以当场按 2.2 的公式计算。四条边各有独立的修正表达式左右边界修正Hy上下边界修正Hx四个角上两条修正叠加。若漏掉某个角的叠加项斜入射下会看到沿 TF/SF 边界向外辐射一圈扇形杂散波角度越大越明显。所有修正式统一用“差分终点在边界内则减终点在边界外则加”的口诀做符号校验。3. Matlab主循环里把Mur和PML都接上边界更新的三种写法二维 FDTD 的主循环只有六个步骤更新Hx、更新Hy、TF/SF 修正、更新Ez、注入源、处理吸收边界。吸收边界有两种接法一阶 Mur 适合快速验证逻辑CPML 适合正式出数据。代码组织上先把内部更新写成独立函数或固定循环再把边界处理放在内层循环之外否则每步迭代都要判断ifMatlab 解释执行会慢到不可接受。3.1 主循环骨架H更新、E更新、源注入、边界的执行顺序for n 1 : Nt % 1) 更新 Hx、Hy内部区域边界由 Mur/PML 之后单独覆盖 for j 2 : Ny-1 for i 2 : Nx-1 Hx(i,j) Hx(i,j) - dt/(mu0*dy) * (Ez(i,j1) - Ez(i,j)); Hy(i,j) Hy(i,j) dt/(mu0*dx) * (Ez(i1,j) - Ez(i,j)); end end % 2) TF/SF 修正见 2.3 节 apply_tfsf_correction(Hx, Hy, Ez_inc, ...); % 3) 更新 Ez for j 2 : Ny-1 for i 2 : Nx-1 Ez(i,j) Ez(i,j) dt/(eps0*dx) * (Hy(i,j) - Hy(i-1,j)) ... - dt/(eps0*dy) * (Hx(i,j) - Hx(i,j-1)); end end % 4) 注入斜入射源源点避免落在 PML 和 TF/SF 边界上 Ez(nx0, ny0) Ez(nx0, ny0) amp * sin(2*pi*f0*n*dt - kx*x(nx0) - ky*y(ny0)); % 5) 吸收边界一阶 Mur 或 CPML apply_boundary(Ez, Hx, Hy, n); % 6) 记录观测点 probe(n) Ez(px, py); end循环从2到N-1是刻意留出最外圈让吸收边界在H与E更新完毕后再覆盖边界值。apply_tfsf_correction如果写成函数调用每次循环都有函数压栈开销建议在正式版本里直接把修正代码内联到主循环或把覆盖范围做成列向量批量运算。apply_boundary同样要避免逐点if判断Mur 和 PML 都只作用于固定索引范围可以拆成边界索引向量一次更新。3.2 一阶Mur边界公式、实现与适用角度一阶 Mur 在 x 方向左边界x0处的离散形式是% 左边界一阶 Mur外法向为 -x % Ez_new(1,j) 由内部相邻列的场外推得到 Ez(1, 2:Ny-1) Ez(2, 2:Ny-1) ... (c*dt - dx) / (c*dt dx) * (Ez(2, 2:Ny-1) - Ez_old(1, 2:Ny-1));Ez_old是上一步的边界值Ez(2, 2:Ny-1)是当前步内部第一列的新值。这个式子等效于对沿-x方向传播的波做一阶泰勒外推正入射时反射系数还能压到 1% 量级斜入射超过 30 度后反射明显上升。到 60 度接近掠射时一阶 Mur 几乎不吸收波会沿着边界表面爬行并重新折回散射场区域在观测点上表现为一个滞后约两倍边界距离的大尾巴。斜入射下Mur 的唯一优点是不需要额外参数和辅助变量几行代码就能跑。如果在调 TF/SF 阶段发现散射场有异常大的边界反射先别急着怀疑 PML——把边界临时换成 Mur若 Mur 的反射曲线在预期范围内说明问题出在 TF/SF 连接项若 Mur 反射比理论值大得多则 TF/SF 边界本身就有符号错误。这是这个 zip 文件名里同时出现Mur和PML的工程意义。3.3 CPML系数计算与辅助变量更新CPML 用坐标拉伸加数字递归卷积替代传统分裂场不需要把Ez分裂成多个子分量内存开销只在每个 PML 层多存几个辅助变量。系数预计算通常写成独立函数function [bx, cx] cpml_profile_1d(Npml, dx, dt, m, sigma_max, kappa_max, alpha_max) % Npml 为 PML 层数m 为渐变指数通常取 3~4 eps0 8.8541878128e-12; bx zeros(1, Npml); cx zeros(1, Npml); for k 1 : Npml eta (k - 0.5) / Npml; % 取单元中心位置避免边界阶跃 sigma sigma_max * eta^m; % 电导率从内到外渐变 kappa 1 (kappa_max - 1) * eta^m; alpha alpha_max * (1 - eta)^m; % alpha 在内部大、外部趋近 0 bx(k) exp(-(sigma/kappa alpha) * dt / eps0); cx(k) sigma / (kappa * (sigma kappa * alpha)) * (bx(k) - 1); end endeta取(k-0.5)/Npml而不是k/Npml是为了让每个 PML 单元的采样点落在单元中间避免电导率在界面处突然跳变。kappa从 1 渐变到kappa_max作用是补偿大角度入射下的相速误差alpha从内部到外部递减用于吸收低频分量但alpha过大会削弱对直流分量的衰减。PML 内Ez的更新要在普通旋度项后面追加辅助变量贡献% 以 x 方向右侧 PML 为例j 为内部行索引 for i Nx-Npml : Nx-1 idx i - (Nx-Npml) 1; psz(idx) bx(idx) * psz(idx) ... cx(idx) * (Ez(i1, j) - Ez(i, j)); Ez(i, j) Ez(i, j) dt/eps0 * psz(idx); endpsz记录的是Ez沿 x 方向差分的历史卷积值bx、cx逐层不同。Hx、Hy在 PML 内也需要各自的辅助变量且四个角和两条边方向交叉处要叠加两个方向的psi项。完整 CPML 的实现细节较多但不要为省事把 PML 厚度设成 4 层以下8 到 12 层是兼顾内存和吸收率的常用起点。4. 斜入射下PML和Mur哪家强用反射曲线调参避免白设边界做了斜入射仿真边界吸收效果必须用数据说话肉眼看场图判断边界好坏会误判。Mur 的反射率随入射角恶化这是理论决定的PML 的反射率则依赖sigma_max、kappa_max、alpha_max三个参数和厚度Npml调参顺序错了会得到一条高频振铃或低频拖尾的能量曲线。4.1 反射误差的度量域内能量时域斜率与观测点残差% 源停止注入后统计计算域内总能量随时间的变化 E_field 0; H_field 0; for n n_start : Nt E_field E_field eps0 * sum(Ez(2:Nx-1, 2:Ny-1).^2, all); H_field H_field mu0 * sum(Hx(2:Nx-1, :).^2 Hy(:, 2:Ny-1).^2, all); energy(n) (E_field H_field) * dx * dy; end semilogy(n_start:Nt, energy(n_start:Nt));能量统计区域要避开 PML 本身只统计内部物理域。把源在n_start时刻关闭后理想情况下总能量应该指数衰减曲线在对数坐标下接近一条直线说明 PML 在工作。若曲线末端反弹上翘是边界反射波已经回到统计区若曲线呈现阶梯状通常是某个方向的psi变量没更新或者 TF/SF 修正泄漏。观测点残差法更直接在散射体后方布置探针记录加 PML 与不加 PML 两种情况下Ez的差用峰值差除以入射波峰值得到反射系数估计。参考解可以用扩大两倍网格、避开反射波到达时间的“截断真实解”替代。4.2 PML参数表与斜入射调参顺序参数常用初值调高效果调过头了Npml8~12 层反射率下降内存和每步耗时线性上涨m3~4渐变更平滑内部电导率过低薄 PML 失效sigma_max(m1)/(150*pi*dx)中频吸收更强离散误差增大反射回升kappa_max1改善掠射角吸收波阻抗偏移正入射反射上升alpha_max0压低低频尾巴直流和近直流分量吸收变差每改一个参数只用一组斜入射角度跑一遍能量曲线先让sigma_max和Npml固定再加kappa最后用alpha修低频。网格dx lambda/20时sigma_max取(m1)/(150*pi*dx)算出来的量级已经足够追求更低反射可以按-ln(R0)*(m1)/(2*eta0*d*Npml)反推R0取 0.001 到 0.0001。4.3 容易翻车的三个现场第一个翻车点是 PML 区域里混入了 TF/SF 边界或激励源。TF/SF 修正要求边界两侧一个在总场区一个在散射场区PML 区域里既没有清晰的总场/散射场之分又有大梯度场辅助变量会产生虚假累积。解决方法是把 TF/SF 矩形缩小保证四周边界距离 PML 至少 5 到 10 个网格。第二个翻车点是单位不统一。sigma_max公式里的dx要同时带入射波频率、光速和介电常数用角频率还是普通频率算kx错了斜入射角就全偏了。第三个翻车点是用 Matlab 的diff函数做边界内更新时数组错位。diff(Ez, 1, 2)返回的数组比原数组在第二维少一列直接加到Hx上会让下标偏移一格正入射时不易察觉斜入射会表现为上下两条边界的反射不对称。5. 一张能量曲线验证二维FDTD斜入射代码直通检验与高斯束进阶边界调完最后做一次直通检验计算域内不放ball散射体斜入射平面波穿过 TF/SF 区域直达 PML。在域的另一侧设一个探针把探针时域波形与解析入射波相减残差峰值应低于入射峰值的 0.1% 到 0.5%取决于网格精度和 PML 厚度。如果残差里出现一个明显滞后于主峰的次级脉冲其滞后时间约等于边界往返距离除以群速度就可以锁定是哪条边界反射回来的。% 直通检验无散射体比较探针场与解析解 err max(abs(probe(:) - probe_analytic(:))) / max(abs(probe_analytic(:)));残差超过 1% 时优先检查 2.3 节四个角上的修正是否都执行了斜入射下四个角的误差会沿着对角线方向向外传播场图上表现为两道斜向条纹。做入射角扫描时不要只测一个角度分别跑 0 度、30 度、60 度三者能量曲线的斜率应当接近一旦某个角度反射率骤增通常是kappa_max不足或 TF/SF 边界离 PML 太近。进阶做法是把平面波源替换成高斯束用不同束腰和入射角合成一组准直波束覆盖一定角度谱。二维代码里高斯束源同样走 TF/SF 连接条件只需要把 2.2 节的Ez_inc从单频平面波换成高斯束的解析表达式PML 参数不用重调。这个技巧在做介质圆柱、多层平板和周期结构的宽角度响应时非常实用能一次算出多个角度下的散射系数而不是逐个角度重跑整个仿真。本文还有配套的精品资源点击获取
二维FDTD斜入射仿真：TF/SF边界与PML实现要点解析
简介这是一份基于MATLAB实现的二维FDTD时域有限差分源代码核心针对平面波斜入射场景并在边界处引入PML完美匹配层以模拟开放空间适合电磁场数值计算学习者、科研人员和工程师快速搭建二维散射与传播仿真实验。代码将时间与空间离散化通过迭代更新电场和磁场完整展示了Maxwell方程在二维网格上的数值求解流程同时也体现了斜入射激励源的设置、PML吸收边界构建以及Mur型边界处理等关键技术。资源压缩包共1个文件文件类型为m脚本整体大小仅3KB小巧精炼便于直接运行、逐行修改和扩展。目前已有320人学习下载足见其在二维FDTD教学与入门实践中的参考价值。获取后可重点研读主程序中的网格配置、源信号加入和PML边界更新三大部分对理解FDTD算法原理及实现细节有直接帮助。1. 从文件名到模块拆解二维FDTD斜入射代码到底先看哪里单看Tot_Src_ball_Mur.zip_FDTD 二维 matlab_Oblique incidence_PML_二维FDTD这一长串文件名真正值得先读的并不是二维 FDTD 主更新循环——主循环在任何教材里都能抄到——而是Tot、Src、ball、Mur这四个缩写背后对应的模块。Tot指 Total field / Scattered field总场/散射场分离Src是入射源ball在二维里对应圆柱散射体的圆截面Mur则是一种如今很少单独用于斜入射仿真的吸收边界。斜入射平面波一旦进入三维空间的全矢量问题网格色散、TF/SF 修正项符号、PML 各向异性都会同时发作。反直觉的结论是这类包里保留一阶 Mur 边界不是作者偷懒而是把 TF/SF 斜入射写对之后拿 Mur 的低吸收率当探针——如果连 Mur 的反射曲线都异常平直多半是连接条件符号反了而不是边界本身出了问题。2. 二维FDTD建斜入射平面波TMz极化、CFL与TF/SF连接条件二维 FDTD 没有真正的“球”ball在这个场景里是介质圆柱或 PEC 圆柱的截面。做这类仿真极化选择直接决定代码长度TMz 极化下电场只有Ez一个标量分量磁场有Hx、Hy两个分量存储和更新都比 TEz 少一套斜入射时也不需要额外处理 z 方向的波矢分量。以下内容按 TMz 展开。2.1 TMz极化下Yee网格的场摆放与更新方程Yee 网格里Ez(i,j)定义在网格点(i,j)上Hx(i,j)定义在(i, j1/2)Hy(i,j)定义在(i1/2, j)。这种半格错位让旋度方程的差分天然中心对称不需要插值。无耗介质中的更新方程写成% TMz 二维 FDTD 内部更新dx 与 dy 可以不等 Hx(i,j) Hx(i,j) - (dt / mu0 / dy) * (Ez(i,j1) - Ez(i,j)); Hy(i,j) Hy(i,j) (dt / mu0 / dx) * (Ez(i1,j) - Ez(i,j)); Ez(i,j) Ez(i,j) (dt / eps0 / dx) * (Hy(i,j) - Hy(i-1,j)) ... - (dt / eps0 / dy) * (Hx(i,j) - Hx(i,j-1));Hx的差分发生在 y 方向Hy的差分发生在 x 方向Ez的更新则由两个旋度项共同组成。若dx dy deltaCFL 条件收紧为dt delta / (c * sqrt(2))实际取 0.99 倍以保证数值稳定。空间步长常见取delta lambda / 20斜入射角度越大波前在网格对角线方向上的投影越长需要把步长加密到lambda / 40或更细否则沿边界的相位误差会在远场方向图里表现为非物理的旁瓣。表格里列出 TMz 下各分量的相对位置和更新方向写代码前先对齐这张表能省掉大量下标调试时间分量位置约定更新依赖边界修正对象Ez(i,j)网格节点相邻Hx、Hy吸收边界直接改该值Hx(i,j)y 方向半格Ez(i,j1)-Ez(i,j)TF/SF 左右边界Hy(i,j)x 方向半格Ez(i1,j)-Ez(i,j)TF/SF 上下边界2.2 斜入射波的波矢量拆分与Matlab源项注入斜入射平面波与正入射的本质区别在于波矢量k不再垂直于 TF/SF 矩形边界。设入射角theta相对 y 轴正向传播方向在 x-y 平面内为 30 度则波矢量在两个坐标轴上的投影为theta 30 * pi / 180; kx 2 * pi * f0 / c0 * sin(theta); ky 2 * pi * f0 / c0 * cos(theta); % 入射方向设为 y坐标轴不同时符号要整体翻转 % 硬源注入只在源点叠加不污染整个计算域 Ez(nx0, ny0) Ez(nx0, ny0) amp * sin(2*pi*f0*n*dt - kx*x(nx0) - ky*y(ny0));f0是载频amp是幅度x、y是网格坐标数组。源点要放在 TF/SF 边界内侧不远处而不是整个计算域逐点赋值——逐点赋解析解会导致网格色散与 FDTD 本征色散叠加波形传播几百步后高频分量明显畸变。斜入射下尤其不能只注入Ez而忽略Hx、Hy的入射分量否则总场区里的等效阻抗不匹配会在源点附近激发二次辐射。2.3 TF/SF边界在斜入射角度下的修正式符号与下标怎么对齐TF/SF 的思想是计算域内层存总场外层存散射场两者之间用一层虚拟边界连接。斜入射时最容易写错的是连接条件的方向。以左边界x ix0为例Hy(ix0-1, iy0:iy1)位于边界外侧它更新时本应使用散射区的Ez但数组里那个位置可能存的是总场因此要减掉入射场% 左边界 xix0Hy 的差分跨过 TF/SF 面减掉入射 Ez Hy(ix0-1, iy0:iy1) Hy(ix0-1, iy0:iy1) ... - dt / (mu0 * dx) * Ez_inc(ix0, iy0:iy1); % 右边界 xix1符号相反加上入射 Ez Hy(ix1, iy0:iy1) Hy(ix1, iy0:iy1) ... dt / (mu0 * dx) * Ez_inc(ix1, iy0:iy1);这里Ez_inc是解析入射场在对应网格上的采样可以预生成一个与Ez同尺寸的数组也可以当场按 2.2 的公式计算。四条边各有独立的修正表达式左右边界修正Hy上下边界修正Hx四个角上两条修正叠加。若漏掉某个角的叠加项斜入射下会看到沿 TF/SF 边界向外辐射一圈扇形杂散波角度越大越明显。所有修正式统一用“差分终点在边界内则减终点在边界外则加”的口诀做符号校验。3. Matlab主循环里把Mur和PML都接上边界更新的三种写法二维 FDTD 的主循环只有六个步骤更新Hx、更新Hy、TF/SF 修正、更新Ez、注入源、处理吸收边界。吸收边界有两种接法一阶 Mur 适合快速验证逻辑CPML 适合正式出数据。代码组织上先把内部更新写成独立函数或固定循环再把边界处理放在内层循环之外否则每步迭代都要判断ifMatlab 解释执行会慢到不可接受。3.1 主循环骨架H更新、E更新、源注入、边界的执行顺序for n 1 : Nt % 1) 更新 Hx、Hy内部区域边界由 Mur/PML 之后单独覆盖 for j 2 : Ny-1 for i 2 : Nx-1 Hx(i,j) Hx(i,j) - dt/(mu0*dy) * (Ez(i,j1) - Ez(i,j)); Hy(i,j) Hy(i,j) dt/(mu0*dx) * (Ez(i1,j) - Ez(i,j)); end end % 2) TF/SF 修正见 2.3 节 apply_tfsf_correction(Hx, Hy, Ez_inc, ...); % 3) 更新 Ez for j 2 : Ny-1 for i 2 : Nx-1 Ez(i,j) Ez(i,j) dt/(eps0*dx) * (Hy(i,j) - Hy(i-1,j)) ... - dt/(eps0*dy) * (Hx(i,j) - Hx(i,j-1)); end end % 4) 注入斜入射源源点避免落在 PML 和 TF/SF 边界上 Ez(nx0, ny0) Ez(nx0, ny0) amp * sin(2*pi*f0*n*dt - kx*x(nx0) - ky*y(ny0)); % 5) 吸收边界一阶 Mur 或 CPML apply_boundary(Ez, Hx, Hy, n); % 6) 记录观测点 probe(n) Ez(px, py); end循环从2到N-1是刻意留出最外圈让吸收边界在H与E更新完毕后再覆盖边界值。apply_tfsf_correction如果写成函数调用每次循环都有函数压栈开销建议在正式版本里直接把修正代码内联到主循环或把覆盖范围做成列向量批量运算。apply_boundary同样要避免逐点if判断Mur 和 PML 都只作用于固定索引范围可以拆成边界索引向量一次更新。3.2 一阶Mur边界公式、实现与适用角度一阶 Mur 在 x 方向左边界x0处的离散形式是% 左边界一阶 Mur外法向为 -x % Ez_new(1,j) 由内部相邻列的场外推得到 Ez(1, 2:Ny-1) Ez(2, 2:Ny-1) ... (c*dt - dx) / (c*dt dx) * (Ez(2, 2:Ny-1) - Ez_old(1, 2:Ny-1));Ez_old是上一步的边界值Ez(2, 2:Ny-1)是当前步内部第一列的新值。这个式子等效于对沿-x方向传播的波做一阶泰勒外推正入射时反射系数还能压到 1% 量级斜入射超过 30 度后反射明显上升。到 60 度接近掠射时一阶 Mur 几乎不吸收波会沿着边界表面爬行并重新折回散射场区域在观测点上表现为一个滞后约两倍边界距离的大尾巴。斜入射下Mur 的唯一优点是不需要额外参数和辅助变量几行代码就能跑。如果在调 TF/SF 阶段发现散射场有异常大的边界反射先别急着怀疑 PML——把边界临时换成 Mur若 Mur 的反射曲线在预期范围内说明问题出在 TF/SF 连接项若 Mur 反射比理论值大得多则 TF/SF 边界本身就有符号错误。这是这个 zip 文件名里同时出现Mur和PML的工程意义。3.3 CPML系数计算与辅助变量更新CPML 用坐标拉伸加数字递归卷积替代传统分裂场不需要把Ez分裂成多个子分量内存开销只在每个 PML 层多存几个辅助变量。系数预计算通常写成独立函数function [bx, cx] cpml_profile_1d(Npml, dx, dt, m, sigma_max, kappa_max, alpha_max) % Npml 为 PML 层数m 为渐变指数通常取 3~4 eps0 8.8541878128e-12; bx zeros(1, Npml); cx zeros(1, Npml); for k 1 : Npml eta (k - 0.5) / Npml; % 取单元中心位置避免边界阶跃 sigma sigma_max * eta^m; % 电导率从内到外渐变 kappa 1 (kappa_max - 1) * eta^m; alpha alpha_max * (1 - eta)^m; % alpha 在内部大、外部趋近 0 bx(k) exp(-(sigma/kappa alpha) * dt / eps0); cx(k) sigma / (kappa * (sigma kappa * alpha)) * (bx(k) - 1); end endeta取(k-0.5)/Npml而不是k/Npml是为了让每个 PML 单元的采样点落在单元中间避免电导率在界面处突然跳变。kappa从 1 渐变到kappa_max作用是补偿大角度入射下的相速误差alpha从内部到外部递减用于吸收低频分量但alpha过大会削弱对直流分量的衰减。PML 内Ez的更新要在普通旋度项后面追加辅助变量贡献% 以 x 方向右侧 PML 为例j 为内部行索引 for i Nx-Npml : Nx-1 idx i - (Nx-Npml) 1; psz(idx) bx(idx) * psz(idx) ... cx(idx) * (Ez(i1, j) - Ez(i, j)); Ez(i, j) Ez(i, j) dt/eps0 * psz(idx); endpsz记录的是Ez沿 x 方向差分的历史卷积值bx、cx逐层不同。Hx、Hy在 PML 内也需要各自的辅助变量且四个角和两条边方向交叉处要叠加两个方向的psi项。完整 CPML 的实现细节较多但不要为省事把 PML 厚度设成 4 层以下8 到 12 层是兼顾内存和吸收率的常用起点。4. 斜入射下PML和Mur哪家强用反射曲线调参避免白设边界做了斜入射仿真边界吸收效果必须用数据说话肉眼看场图判断边界好坏会误判。Mur 的反射率随入射角恶化这是理论决定的PML 的反射率则依赖sigma_max、kappa_max、alpha_max三个参数和厚度Npml调参顺序错了会得到一条高频振铃或低频拖尾的能量曲线。4.1 反射误差的度量域内能量时域斜率与观测点残差% 源停止注入后统计计算域内总能量随时间的变化 E_field 0; H_field 0; for n n_start : Nt E_field E_field eps0 * sum(Ez(2:Nx-1, 2:Ny-1).^2, all); H_field H_field mu0 * sum(Hx(2:Nx-1, :).^2 Hy(:, 2:Ny-1).^2, all); energy(n) (E_field H_field) * dx * dy; end semilogy(n_start:Nt, energy(n_start:Nt));能量统计区域要避开 PML 本身只统计内部物理域。把源在n_start时刻关闭后理想情况下总能量应该指数衰减曲线在对数坐标下接近一条直线说明 PML 在工作。若曲线末端反弹上翘是边界反射波已经回到统计区若曲线呈现阶梯状通常是某个方向的psi变量没更新或者 TF/SF 修正泄漏。观测点残差法更直接在散射体后方布置探针记录加 PML 与不加 PML 两种情况下Ez的差用峰值差除以入射波峰值得到反射系数估计。参考解可以用扩大两倍网格、避开反射波到达时间的“截断真实解”替代。4.2 PML参数表与斜入射调参顺序参数常用初值调高效果调过头了Npml8~12 层反射率下降内存和每步耗时线性上涨m3~4渐变更平滑内部电导率过低薄 PML 失效sigma_max(m1)/(150*pi*dx)中频吸收更强离散误差增大反射回升kappa_max1改善掠射角吸收波阻抗偏移正入射反射上升alpha_max0压低低频尾巴直流和近直流分量吸收变差每改一个参数只用一组斜入射角度跑一遍能量曲线先让sigma_max和Npml固定再加kappa最后用alpha修低频。网格dx lambda/20时sigma_max取(m1)/(150*pi*dx)算出来的量级已经足够追求更低反射可以按-ln(R0)*(m1)/(2*eta0*d*Npml)反推R0取 0.001 到 0.0001。4.3 容易翻车的三个现场第一个翻车点是 PML 区域里混入了 TF/SF 边界或激励源。TF/SF 修正要求边界两侧一个在总场区一个在散射场区PML 区域里既没有清晰的总场/散射场之分又有大梯度场辅助变量会产生虚假累积。解决方法是把 TF/SF 矩形缩小保证四周边界距离 PML 至少 5 到 10 个网格。第二个翻车点是单位不统一。sigma_max公式里的dx要同时带入射波频率、光速和介电常数用角频率还是普通频率算kx错了斜入射角就全偏了。第三个翻车点是用 Matlab 的diff函数做边界内更新时数组错位。diff(Ez, 1, 2)返回的数组比原数组在第二维少一列直接加到Hx上会让下标偏移一格正入射时不易察觉斜入射会表现为上下两条边界的反射不对称。5. 一张能量曲线验证二维FDTD斜入射代码直通检验与高斯束进阶边界调完最后做一次直通检验计算域内不放ball散射体斜入射平面波穿过 TF/SF 区域直达 PML。在域的另一侧设一个探针把探针时域波形与解析入射波相减残差峰值应低于入射峰值的 0.1% 到 0.5%取决于网格精度和 PML 厚度。如果残差里出现一个明显滞后于主峰的次级脉冲其滞后时间约等于边界往返距离除以群速度就可以锁定是哪条边界反射回来的。% 直通检验无散射体比较探针场与解析解 err max(abs(probe(:) - probe_analytic(:))) / max(abs(probe_analytic(:)));残差超过 1% 时优先检查 2.3 节四个角上的修正是否都执行了斜入射下四个角的误差会沿着对角线方向向外传播场图上表现为两道斜向条纹。做入射角扫描时不要只测一个角度分别跑 0 度、30 度、60 度三者能量曲线的斜率应当接近一旦某个角度反射率骤增通常是kappa_max不足或 TF/SF 边界离 PML 太近。进阶做法是把平面波源替换成高斯束用不同束腰和入射角合成一组准直波束覆盖一定角度谱。二维代码里高斯束源同样走 TF/SF 连接条件只需要把 2.2 节的Ez_inc从单频平面波换成高斯束的解析表达式PML 参数不用重调。这个技巧在做介质圆柱、多层平板和周期结构的宽角度响应时非常实用能一次算出多个角度下的散射系数而不是逐个角度重跑整个仿真。本文还有配套的精品资源点击获取
