简介FDTD时域有限差分是电磁场数值计算中广泛使用的方法这份压缩包提供了基于MATLAB的一维、二维与三维FDTD仿真程序适合电磁场初学者、光电或微波工程方向的学生与科研人员快速上手理解从算法到图形化结果的完整流程无论是一维波导、二维平面结构还是三维复杂器件均可对照学习。资源共10个文件以6个.m脚本为主另含2个.asv备份文件、1个.fig图形文件压缩包仅87KB下载后配合MATLAB即可运行。已有703人学习下载。代码覆盖一维到三维的差分计算并带有动态图形呈现能直观展示电磁场演化过程便于观察脉冲在介质中的传播、反射和透射.fig文件提供界面交互方便调整参数查看结果。通过研读这些脚本可掌握FDTD网格设置、边界条件、迭代更新等核心步骤为后续复杂结构电磁仿真和天线、波导等器件设计打下基础。1. 一套Matlab FDTD代码能做什么从一维脉冲到三维结构某同学做微带天线设计手里只有一台普通笔记本。商业仿真软件算一次要几分钟想改个尺寸看看趋势等待时间比建模时间还长。他后来拿到一套Matlab写的FDTD代码包一维两维三维都有才真正体会到什么叫“自己掌控仿真过程”。这里说的FDTD全称是时域有限差分法思路很直白把麦克斯韦方程组里对时间和空间的偏导换成差分然后在网格上一格一格地推进。你不需要求逆矩阵不需要做频域变换只需要不断迭代电场和磁场就会在计算域里跑起来。这套标题里能拆出几个关键词matlab_fdtd说明主语言是MATLAB二维、fdtd三维、_一维FDTD说明代码覆盖了三个维度fdtd则代表了另一个方向——当Matlab算不动超大规模模型时C版具备接棒的潜力。对电磁场新手来说这套代码是理解数值仿真的最佳入口对射频工程师来说它是验证结构趋势的快速实验台对做研究的人来说它是改边界条件、加材料模型、试各种激励源的自由沙盒。2. 一维FDTDYee网格、迭代公式与最小可运行代码2.1 为什么一维是理解FDTD的最佳起点一维FDTD虽然看上去简单但FDTD最核心的两个思想都在里面空间上电场和磁场错开半个网格时间上电场和磁场错开半个时间步。这种错位排布让“用中心差分替代偏导”成为可能而且得到的是显式迭代格式——新时刻的场值只依赖上一时刻的相邻点不需要求解任何方程组。一维情况下只保留两个分量电场Ez和磁场Hy沿z轴传播的平面电磁波电场沿x方向极化磁场沿y方向。麦克斯韦方程组简化成两个标量方程离散后得到% 一维FDTD核心迭代归一化单位c01, dx1, dt0.7 % 更新磁场Hy(i) 由相邻电场差分驱动 hy(1:end-1) hy(1:end-1) - dt * (ez(2:end) - ez(1:end-1)); % 更新电场Ez(i) 由相邻磁场差分驱动 ez(2:end) ez(2:end) - dt * (hy(2:end) - hy(1:end-1));这段代码里hy(1:end-1)对应半网格位置的磁场ez(2:end)和ez(1:end-1)的差值正好是电场在空间上的差分系数只有一个dt因为归一化单位下c01、dx1。实际工程中如果dx取2mmc0取3e8dt的物理值需要按CFL条件重新计算但代码结构一模一样。时间步长dt的取值不是随便定的。Courant稳定性条件要求dt dx/c0一维时取0.8或0.9倍的dx/c0是常见做法。取小了浪费计算量取大了仿真直接发散。你会发现一个规律所有初学者都会在某个时刻把dt调大试图加速然后看到电场值指数增长这就是FDTD第一课。2.2 二十行代码跑通高斯脉冲下面是完整的一维FDTD程序模拟高斯脉冲在自由空间传播并在两端被吸收边界吸收% 一维FDTD高斯脉冲在自由空间中的传播 c0 1; % 归一化光速 dx 1; % 空间步长 dt 0.7; % 时间步长满足CFL条件 Nx 800; % 空间网格数 Nt 1000; % 时间步数 ez zeros(1, Nx); % 电场 hy zeros(1, Nx); % 磁场 % 高斯脉冲源参数 t0 50; % 峰值出现的时间步 width 15; % 脉冲宽度时间步 source_pos 200; % 源所在网格位置 for n 1:Nt % 磁场更新所有内部点一次完成 hy(1:end-1) hy(1:end-1) - dt * (ez(2:end) - ez(1:end-1)); % 电场更新注意索引错位 ez(2:end) ez(2:end) - dt * (hy(2:end) - hy(1:end-1)); % 加载高斯源在源位置叠加一个脉冲 ez(source_pos) ez(source_pos) exp(-((n - t0) / width)^2); % 一阶吸收边界让边缘场值向内部“看齐” ez(1) ez(2); ez(end) ez(end-1); end % 可视化最终波形 plot(ez, b-, LineWidth, 1.5); xlabel(网格位置); ylabel(电场幅度); title(一维FDTD高斯脉冲传播结果);代码逻辑说明磁场先更新因为电场在时间推进上依赖磁场更新完磁场后立刻更新电场此时用到的是当前时刻的磁场值这是FDTD“蛙跳式”推进的典型写法。源加载放在电场更新之后相当于在当前时刻对源位置做一个增量修正。边界处理是纯经验的——把边缘网格的场值复制给最外面的点让波“走出去”时反射小到可以忽略。参数说明t050意味着第50步时高斯函数达到峰值width15控制脉冲的空间宽度太窄的脉冲会激发高频分量而在粗网格上高频分量色散严重波形会拖尾。source_pos200选择在计算域靠中间偏左的位置保证脉冲在Nt1000步内有足够空间传播并观察边界行为。2.3 三个必调参数与它们的作用第一个是网格分辨率dx。在归一化代码里dx1但实际物理意义是“最小仿真波长至少覆盖10个网格”。如果你要仿真2.4GHz的微波信号自由空间波长约125mmdx取12.5mm以下就够用如果结构里有介质波长缩短dx还要更小。网格越细计算域越大内存和时间成本按维度指数上涨这是FDTD最残酷的物理规律。第二个是时间步长dt。前面说了CFL条件这里补充一个经验dt通常取0.5到0.9倍的极限值。取0.5安全性高但计算慢取0.9逼近稳定边界如果介质不均匀、网格不等距很容易触发局部不稳定。我一般先取0.7跑通后再试着拉到0.85看结果是否明显变化如果没有就保留。第三个是吸收边界。一阶吸收边界就是代码里把边缘值复制给相邻点的做法只对垂直入射有效斜入射时会有可见反射。到了一维以上这句话开始成为质量瓶颈。所以二维FDTD通常直接上PML一维代码里用简单吸收边界只是为了先把逻辑跑通。当你看到脉冲遇到边界后反射幅度远小于主波说明边界有效如果反射大得一眼能看到回波就需要检查边界实现。3. 二维FDTD把PML和点源一次搞清楚3.1 TM模式与三个场分量的排布二维FDTD比一维多了一个空间维度场分量从2个变成3个。这里选择TM模式横磁波电场只有一个z方向分量Ez磁场有x和y方向分量Hx、Hy。这种模式的Yee网格排布是Ez位于网格中心Hx位于Ez的y方向偏移半个网格的位置Hy位于Ez的x方向偏移半个网格的位置。迭代公式有三个% 二维TM模式FDTD核心迭代 % Hx更新由Ez在y方向的差分别驱动 hx(1:end-1, :) hx(1:end-1, :) - dt * (ez(2:end, :) - ez(1:end-1, :)); % Hy更新由Ez在x方向的差分别驱动 hy(:, 1:end-1) hy(:, 1:end-1) dt * (ez(:, 2:end) - ez(:, 1:end-1)); % Ez更新由Hx和Hy共同驱动 ez(2:end, 2:end) ez(2:end, 2:end) - dt * (hy(2:end, 2:end) - hy(2:end, 1:end-1)); ez(2:end, 2:end) ez(2:end, 2:end) dt * (hx(2:end, 2:end) - hx(1:end-1, 2:end));注意Hy更新时差分前面是正号Hx是负号这取决于你定义的波传播方向和坐标轴方向。我见过不少人把这里符号搞反结果电场不传播而是原地震荡。这个细节最容易排查跑一个点源看波前是否对称向外扩散如果出现“只在一个方向走另一个方向不动”大概率是这里的符号或索引错位。3.2 PML吸收边界参数怎么设才算合格一维的简单吸收边界到二维就不够用了。斜入射时它会产生明显反射标准做法是加PML完美匹配层。PML的原理是在计算域外围构造一个特殊的吸收介质层让入射波在层内衰减而不反射回计算域内部。要点在于电导率σ必须从内到外渐变不能突变否则突变界面本身就是反射源。% PML参数设置 npml 10; % PML层数 sigma_max 0.8; % 最大电导率归一化单位 % 渐变系数从内到外按多项式分布 for i 1:npml % 靠近计算域的PML内部sigma接近0最外层sigma最大 sigma(i) sigma_max * ((i - 0.5) / npml)^2; endPML参数里最重要的两个判断sigma_max取0.8到1.0之间层数取8到16层。sigma_max太小吸收不充分边界反射大sigma_max太大由于离散误差PML本身会反射。层数太薄渐变不够平滑同样反射。我的一般做法是先把sigma_max固定在0.8层数从10开始逐步增加到16观察电场总能量随时间的变化曲线如果平稳下降没有抖动说明边界做得不错。新手经常会犯一个错误PML区域的场更新和计算域内部完全分开写结果PML和内部的交界面因为系数突变产生反射。正确做法是PML区域也参与正常的FDTD迭代只是在原来的迭代公式后面额外加上电导率相关的衰减项。3.3 点源激励与观察点的设置二维FDTD里最典型的激励是高斯脉冲点源直接加在Ez分量的某个网格点上% 高斯脉冲点源激励 ez(src_x, src_y) ez(src_x, src_y) exp(-((n - t0) / width)^2);这个点源产生的波是圆心对称的圆柱波适合观察传播特性。但如果你想知道某个方向的远场方向图点源就无能为力了此时需要引入TFSF总场散射场技术把计算域分成总场区域和散射场区域在分界面上加载平面波。这个属于进阶内容初学阶段点源能把PML调试通就已经迈过最大的坎。观察点的设置也有讲究。不要在PML层内放观察点——那里的场已经被衰减了结果没有参考意义。应该在计算域中心区域放若干个观察点记录场值随时间变化% 在计算域中心记录场值随时间变化 probe_x 50; probe_y 50; % 观察点位置避开PML区 ez_history(n) ez(probe_x, probe_y);观察点位置的选取直接影响你判断仿真是否正确的效率。我习惯在源的正对面放一个观察点在侧面放一个这样能同时看到两路传播时间差间接验证介质区域的等效介电常数是否正确。3.4 快速验证二维代码是否可靠跑完一个二维FDTD后不要急着写报告。先做两件事第一计算总的电场能量看曲线是否呈现先上升后下降的趋势下降段是否平滑如果有平台期说明有持续反射第二加大PML层数看观察点的波形变化是否在可接受范围。如果加大PML后波形明显改变说明原来的PML反射污染了结果需要增加层数或调整sigma_max。更严格的验证是和解析解对比比如均匀介质中点源的辐射场理论解是零阶汉开尔函数但这对初学者而言比较难操作。务实做法是和一维结果对比同一个高斯脉冲在二维里沿某一条线上的传播波形与一维结果应该几乎一样。这个对比能一次性发现源加载问题、索引错位问题和PML参数问题。4. 三维FDTD六分量迭代与内存优化4.1 三维Yee元胞与六个场分量的空间排布三维FDTD是电磁仿真的完整形态电场三个分量Ex、Ey、Ez和磁场三个分量Hx、Hy、Hz交错排列在Yee元胞的棱和面上。三维Yee网格中每个电场分量沿其方向偏移半个网格磁场分量在垂直的空间对角位置偏移半个网格。三维迭代公式很长但结构一致。以Ex为例% 三维FDTD中Ex分量的更新示意 ex(i, j, k) ex(i, j, k) - dt * (hz(i, j, k) - hz(i, j-1, k)) / dy; ex(i, j, k) ex(i, j, k) dt * (hy(i, j, k) - hy(i, j, k-1)) / dz;每个分量更新都涉及两个相邻磁场分量的差分整个循环需要遍历i、j、k三个维度。在Matlab里如果直接写三重for循环运行速度会慢得让人怀疑人生。常见优化是把内层操作向量化每次更新一个二维切片% 向量化更新每次更新整个x方向的切片 ex(:, 2:end, :) ex(:, 2:end, :) - dt/dy * (hz(:, 2:end, :) - hz(:, 1:end-1, :));这种方式在网格规模不大时效率提升明显但缺点是代码可读性下降出了错不好查。我个人的习惯是一维和二维用向量化写法提升速度三维则优先保证逻辑清晰因为你花在调试上的时间往往比多跑几个小时更值钱。4.2 内存估算为什么100×100×100网格就让人头疼三维FDTD的内存消耗按网格数的三次方增长这是它最让人头疼的地方。六个场分量加上PML辅助变量每个变量都要存储一个和原网格一样大的三维数组。以100×100×100网格为例数据项变量数双精度内存约电场分量33 × 8MB 24MB磁场分量33 × 8MB 24MBPML辅助变量通常4~6个4 × 8MB 32MB合计10约80MB看着好像不多但这只是100³网格。到了200×200×200内存直接跳到640MB300×300×300就是2.2GB。实际工程里三维模型很少是简单立方体结构细节会让网格数轻松破百万内存消耗就奔着几十GB去了。所以三维FDTD的第一个准则是能用二维模拟就不要上三维。如果一个结构在一个方向上没有变化比如无限长圆柱、均匀波导截面完全可以降维成二维模型结果基本一致计算量却少了几个数量级。很多初学者一上来就想做三维全波仿真结果笔记本内存直接报错这是完全可以避免的。第二个准则是使用单精度。Matlab默认的double类型精度对FDTD绝大多数场景来说其实过剩了场值的精度主要受限于数值色散而非浮点误差。用single类型存储场分量内存减半速度经常还有提升% 用单精度减少内存占用 ez single(zeros(nx, ny, nz)); hy single(zeros(nx, ny, nz));注意源加载和观察点输出时可能需要转回double避免精度损失。4.3 三维FDTD的降维时机与混合策略实战中真正需要完整三维FDTD的场景其实有限。微带天线、贴片天线这类薄结构厚度方向的场变化很小可以在厚度方向用粗网格平面方向用细网格。波导结构如果长度远大于横截面可以只仿真一个截面再用模式展开。如果结构具备旋转对称性还可以用BOR体旋转对称方法把三维问题降成二维计算加一个角向模式展开。我见过最直接的做法是用二维FDTD先扫参锁定结构的大致尺寸范围再对最关键的几个尺寸做一次三维验证。这样既控制了计算量又保留三维结果的完整信息。这种“二维扫趋势、三维验终值”的工作流比一上来就跑三维省下太多时间。如果你真的需要大规模三维仿真Matlab脚本语言就力不从心了此时应该转向C或CUDA实现甚至考虑GPU加速。标题里的fdtd就是这个方向——C写的高性能FDTD库。Matlab版本适合学习和中小规模验证大规模生产型仿真还是交给编译型语言更合适。5. 避坑从一维到三维最常见的五个翻车场景5.1 仿真爆炸电场值疯狂增长现象跑了没几步场值变成NaN或者正负无穷大波形完全失控。原因几乎都是CFL条件不满足。dt取值超过稳定性极限数值误差随时间步进指数放大。解决先检查dt是否满足Courant条件。三维情形下dt 1 / (c * sqrt(1/dx² 1/dy² 1/dz²))很多新手只按一维公式dt dx/c来算在二维三维就直接炸了。另外一个隐藏坑介质区域中波的传播速度是c/sqrt(εr)比真空中慢有人把介质区域的dx改小了但忘记同步改小dt照样发散。5.2 二维结果和一维对不上现象同一个高斯脉冲在一维和二维中沿着传播方向观察波形形状不同。原因二维点源产生的波前是圆柱波幅度随距离衰减不像一维平面波幅度不变。另外二维代码的场分量顺序如果从TM换成TE观察点采的场分量不同波形自然不同。解决对比时先用同一条直线上的场分布且把二维幅度乘以sqrt(r)再做归一化补偿圆柱波的几何衰减。如果幅度曲线形态还是对不上检查Hx和Hy的更新符号是否搞反。5.3 PML反射残留观察点波形出现二次回波现象主脉冲过去之后波形尾部出现一个明显的起伏像是波被弹回来了。原因PML层数太少或sigma_max取值不当导致PML没有完全吸收入射波。另一个常见原因是PML和计算域之间的介电常数或电导率突变。解决把PML层数从10增加到15或20sigma_max从0.8调到1.0观察尾部反射是否变小。如果反射没有随PML变强而降低说明PML实现有逻辑错误而不是参数问题。此时先把结构简化成均匀介质只保留PML跑点源看反射排除结构散射的干扰。5.4 Matlab内存不足或运行极慢现象三维FDTD代码一运行就报内存不足或者二维代码跑了一个小时还没完。原因三维数组默认double类型内存翻倍代码里在循环中不断对数组进行重新分配每一次都会复制整个数组。解决第一步把所有场变量声明为single类型。第二步确保所有数组在循环前一次性分配好不要在循环内用扩展赋值比如a(end1) x这种写法。第三步检查是否需要三维仿真——很多问题用二维就能回答这是最有效的优化。5.5 结果和解析解偏差大现象仿真的反射系数或透射系数和理论值差了10%以上。原因网格太粗数值色散严重或者源加载有问题比如高斯脉冲的幅度归一化不对导致能量统计失真。解决把dx减半重新跑一次看结果是否明显改变。如果改变超过几个百分点说明收敛性还不行需要加密网格。此外在频域做后处理时要确认时域信号时长足够长避免截断导致的频谱泄漏。通常取Nt使波传播到计算域边缘并完全被PML吸收再截取场值做FFT。6. 验证你的FDTD结果从网格收敛到能量守恒6.1 网格细化测试检验离散误差的基本方法FDTD结果对不对最直接的办法是网格收敛测试。把dx减半同时dt按CFL关系减半重跑一次比较关键观察点处的波形。如果两次结果几乎重叠说明网格已经细化到误差可以忽略如果结果明显变化说明网格还不够细当前结果不能被信任。% 网格收敛性验证比较不同dx下观察点波形 dx1 0.02; dx2 0.01; % 两组网格尺度 % 分别运行FDTD得到ez1和ez2 % 计算两组结果的最大差异 diff_max max(abs(ez1 - ez2(1:2:end))); disp([网格细化后的最大差异: , num2str(diff_max)]);这个做法计算量不低但值得在关键仿真上做一次。通常做一组粗网格和一组细网格的对比就够了不需要做完整的三组收敛率分析。6.2 能量守恒监控发现隐藏反射的哨兵另一个低成本验证手段是监控计算域电场总能量。在没有源持续激励的情况下总能量应该随时间单调递减——波不断被PML吸收。如果能量曲线出现平段甚至上升说明PML有反射或计算域内部有非物理的增益。% 每次迭代后计算总能量 total_energy(n) sum(sum(sum(ez.^2 ex.^2 ey.^2)));这段代码在二维里对应sum(sum(ez.^2))三维再加两个电场分量。能量曲线平滑下降就是好结果如果你看到曲线尾部翘起来优先查PML参数。6.3 我自己的验证习惯每次写完三维FDTD代码我习惯先在一个一维测试用例上验证整个流程网格收敛测试、能量监控、解析解对比三个都过了才开始正式仿真。这个习惯帮我挡掉过很多次debug问题——常常是电磁方程组本身的符号问题而不是PML或边界的问题。先跑一维因为你一眼能看出波形对不对确认底层迭代无误再升到二维和三维。希望帮到你。本文还有配套的精品资源点击获取
简介FDTD时域有限差分是电磁场数值计算中广泛使用的方法这份压缩包提供了基于MATLAB的一维、二维与三维FDTD仿真程序适合电磁场初学者、光电或微波工程方向的学生与科研人员快速上手理解从算法到图形化结果的完整流程无论是一维波导、二维平面结构还是三维复杂器件均可对照学习。资源共10个文件以6个.m脚本为主另含2个.asv备份文件、1个.fig图形文件压缩包仅87KB下载后配合MATLAB即可运行。已有703人学习下载。代码覆盖一维到三维的差分计算并带有动态图形呈现能直观展示电磁场演化过程便于观察脉冲在介质中的传播、反射和透射.fig文件提供界面交互方便调整参数查看结果。通过研读这些脚本可掌握FDTD网格设置、边界条件、迭代更新等核心步骤为后续复杂结构电磁仿真和天线、波导等器件设计打下基础。1. 一套Matlab FDTD代码能做什么从一维脉冲到三维结构某同学做微带天线设计手里只有一台普通笔记本。商业仿真软件算一次要几分钟想改个尺寸看看趋势等待时间比建模时间还长。他后来拿到一套Matlab写的FDTD代码包一维两维三维都有才真正体会到什么叫“自己掌控仿真过程”。这里说的FDTD全称是时域有限差分法思路很直白把麦克斯韦方程组里对时间和空间的偏导换成差分然后在网格上一格一格地推进。你不需要求逆矩阵不需要做频域变换只需要不断迭代电场和磁场就会在计算域里跑起来。这套标题里能拆出几个关键词matlab_fdtd说明主语言是MATLAB二维、fdtd三维、_一维FDTD说明代码覆盖了三个维度fdtd则代表了另一个方向——当Matlab算不动超大规模模型时C版具备接棒的潜力。对电磁场新手来说这套代码是理解数值仿真的最佳入口对射频工程师来说它是验证结构趋势的快速实验台对做研究的人来说它是改边界条件、加材料模型、试各种激励源的自由沙盒。2. 一维FDTDYee网格、迭代公式与最小可运行代码2.1 为什么一维是理解FDTD的最佳起点一维FDTD虽然看上去简单但FDTD最核心的两个思想都在里面空间上电场和磁场错开半个网格时间上电场和磁场错开半个时间步。这种错位排布让“用中心差分替代偏导”成为可能而且得到的是显式迭代格式——新时刻的场值只依赖上一时刻的相邻点不需要求解任何方程组。一维情况下只保留两个分量电场Ez和磁场Hy沿z轴传播的平面电磁波电场沿x方向极化磁场沿y方向。麦克斯韦方程组简化成两个标量方程离散后得到% 一维FDTD核心迭代归一化单位c01, dx1, dt0.7 % 更新磁场Hy(i) 由相邻电场差分驱动 hy(1:end-1) hy(1:end-1) - dt * (ez(2:end) - ez(1:end-1)); % 更新电场Ez(i) 由相邻磁场差分驱动 ez(2:end) ez(2:end) - dt * (hy(2:end) - hy(1:end-1));这段代码里hy(1:end-1)对应半网格位置的磁场ez(2:end)和ez(1:end-1)的差值正好是电场在空间上的差分系数只有一个dt因为归一化单位下c01、dx1。实际工程中如果dx取2mmc0取3e8dt的物理值需要按CFL条件重新计算但代码结构一模一样。时间步长dt的取值不是随便定的。Courant稳定性条件要求dt dx/c0一维时取0.8或0.9倍的dx/c0是常见做法。取小了浪费计算量取大了仿真直接发散。你会发现一个规律所有初学者都会在某个时刻把dt调大试图加速然后看到电场值指数增长这就是FDTD第一课。2.2 二十行代码跑通高斯脉冲下面是完整的一维FDTD程序模拟高斯脉冲在自由空间传播并在两端被吸收边界吸收% 一维FDTD高斯脉冲在自由空间中的传播 c0 1; % 归一化光速 dx 1; % 空间步长 dt 0.7; % 时间步长满足CFL条件 Nx 800; % 空间网格数 Nt 1000; % 时间步数 ez zeros(1, Nx); % 电场 hy zeros(1, Nx); % 磁场 % 高斯脉冲源参数 t0 50; % 峰值出现的时间步 width 15; % 脉冲宽度时间步 source_pos 200; % 源所在网格位置 for n 1:Nt % 磁场更新所有内部点一次完成 hy(1:end-1) hy(1:end-1) - dt * (ez(2:end) - ez(1:end-1)); % 电场更新注意索引错位 ez(2:end) ez(2:end) - dt * (hy(2:end) - hy(1:end-1)); % 加载高斯源在源位置叠加一个脉冲 ez(source_pos) ez(source_pos) exp(-((n - t0) / width)^2); % 一阶吸收边界让边缘场值向内部“看齐” ez(1) ez(2); ez(end) ez(end-1); end % 可视化最终波形 plot(ez, b-, LineWidth, 1.5); xlabel(网格位置); ylabel(电场幅度); title(一维FDTD高斯脉冲传播结果);代码逻辑说明磁场先更新因为电场在时间推进上依赖磁场更新完磁场后立刻更新电场此时用到的是当前时刻的磁场值这是FDTD“蛙跳式”推进的典型写法。源加载放在电场更新之后相当于在当前时刻对源位置做一个增量修正。边界处理是纯经验的——把边缘网格的场值复制给最外面的点让波“走出去”时反射小到可以忽略。参数说明t050意味着第50步时高斯函数达到峰值width15控制脉冲的空间宽度太窄的脉冲会激发高频分量而在粗网格上高频分量色散严重波形会拖尾。source_pos200选择在计算域靠中间偏左的位置保证脉冲在Nt1000步内有足够空间传播并观察边界行为。2.3 三个必调参数与它们的作用第一个是网格分辨率dx。在归一化代码里dx1但实际物理意义是“最小仿真波长至少覆盖10个网格”。如果你要仿真2.4GHz的微波信号自由空间波长约125mmdx取12.5mm以下就够用如果结构里有介质波长缩短dx还要更小。网格越细计算域越大内存和时间成本按维度指数上涨这是FDTD最残酷的物理规律。第二个是时间步长dt。前面说了CFL条件这里补充一个经验dt通常取0.5到0.9倍的极限值。取0.5安全性高但计算慢取0.9逼近稳定边界如果介质不均匀、网格不等距很容易触发局部不稳定。我一般先取0.7跑通后再试着拉到0.85看结果是否明显变化如果没有就保留。第三个是吸收边界。一阶吸收边界就是代码里把边缘值复制给相邻点的做法只对垂直入射有效斜入射时会有可见反射。到了一维以上这句话开始成为质量瓶颈。所以二维FDTD通常直接上PML一维代码里用简单吸收边界只是为了先把逻辑跑通。当你看到脉冲遇到边界后反射幅度远小于主波说明边界有效如果反射大得一眼能看到回波就需要检查边界实现。3. 二维FDTD把PML和点源一次搞清楚3.1 TM模式与三个场分量的排布二维FDTD比一维多了一个空间维度场分量从2个变成3个。这里选择TM模式横磁波电场只有一个z方向分量Ez磁场有x和y方向分量Hx、Hy。这种模式的Yee网格排布是Ez位于网格中心Hx位于Ez的y方向偏移半个网格的位置Hy位于Ez的x方向偏移半个网格的位置。迭代公式有三个% 二维TM模式FDTD核心迭代 % Hx更新由Ez在y方向的差分别驱动 hx(1:end-1, :) hx(1:end-1, :) - dt * (ez(2:end, :) - ez(1:end-1, :)); % Hy更新由Ez在x方向的差分别驱动 hy(:, 1:end-1) hy(:, 1:end-1) dt * (ez(:, 2:end) - ez(:, 1:end-1)); % Ez更新由Hx和Hy共同驱动 ez(2:end, 2:end) ez(2:end, 2:end) - dt * (hy(2:end, 2:end) - hy(2:end, 1:end-1)); ez(2:end, 2:end) ez(2:end, 2:end) dt * (hx(2:end, 2:end) - hx(1:end-1, 2:end));注意Hy更新时差分前面是正号Hx是负号这取决于你定义的波传播方向和坐标轴方向。我见过不少人把这里符号搞反结果电场不传播而是原地震荡。这个细节最容易排查跑一个点源看波前是否对称向外扩散如果出现“只在一个方向走另一个方向不动”大概率是这里的符号或索引错位。3.2 PML吸收边界参数怎么设才算合格一维的简单吸收边界到二维就不够用了。斜入射时它会产生明显反射标准做法是加PML完美匹配层。PML的原理是在计算域外围构造一个特殊的吸收介质层让入射波在层内衰减而不反射回计算域内部。要点在于电导率σ必须从内到外渐变不能突变否则突变界面本身就是反射源。% PML参数设置 npml 10; % PML层数 sigma_max 0.8; % 最大电导率归一化单位 % 渐变系数从内到外按多项式分布 for i 1:npml % 靠近计算域的PML内部sigma接近0最外层sigma最大 sigma(i) sigma_max * ((i - 0.5) / npml)^2; endPML参数里最重要的两个判断sigma_max取0.8到1.0之间层数取8到16层。sigma_max太小吸收不充分边界反射大sigma_max太大由于离散误差PML本身会反射。层数太薄渐变不够平滑同样反射。我的一般做法是先把sigma_max固定在0.8层数从10开始逐步增加到16观察电场总能量随时间的变化曲线如果平稳下降没有抖动说明边界做得不错。新手经常会犯一个错误PML区域的场更新和计算域内部完全分开写结果PML和内部的交界面因为系数突变产生反射。正确做法是PML区域也参与正常的FDTD迭代只是在原来的迭代公式后面额外加上电导率相关的衰减项。3.3 点源激励与观察点的设置二维FDTD里最典型的激励是高斯脉冲点源直接加在Ez分量的某个网格点上% 高斯脉冲点源激励 ez(src_x, src_y) ez(src_x, src_y) exp(-((n - t0) / width)^2);这个点源产生的波是圆心对称的圆柱波适合观察传播特性。但如果你想知道某个方向的远场方向图点源就无能为力了此时需要引入TFSF总场散射场技术把计算域分成总场区域和散射场区域在分界面上加载平面波。这个属于进阶内容初学阶段点源能把PML调试通就已经迈过最大的坎。观察点的设置也有讲究。不要在PML层内放观察点——那里的场已经被衰减了结果没有参考意义。应该在计算域中心区域放若干个观察点记录场值随时间变化% 在计算域中心记录场值随时间变化 probe_x 50; probe_y 50; % 观察点位置避开PML区 ez_history(n) ez(probe_x, probe_y);观察点位置的选取直接影响你判断仿真是否正确的效率。我习惯在源的正对面放一个观察点在侧面放一个这样能同时看到两路传播时间差间接验证介质区域的等效介电常数是否正确。3.4 快速验证二维代码是否可靠跑完一个二维FDTD后不要急着写报告。先做两件事第一计算总的电场能量看曲线是否呈现先上升后下降的趋势下降段是否平滑如果有平台期说明有持续反射第二加大PML层数看观察点的波形变化是否在可接受范围。如果加大PML后波形明显改变说明原来的PML反射污染了结果需要增加层数或调整sigma_max。更严格的验证是和解析解对比比如均匀介质中点源的辐射场理论解是零阶汉开尔函数但这对初学者而言比较难操作。务实做法是和一维结果对比同一个高斯脉冲在二维里沿某一条线上的传播波形与一维结果应该几乎一样。这个对比能一次性发现源加载问题、索引错位问题和PML参数问题。4. 三维FDTD六分量迭代与内存优化4.1 三维Yee元胞与六个场分量的空间排布三维FDTD是电磁仿真的完整形态电场三个分量Ex、Ey、Ez和磁场三个分量Hx、Hy、Hz交错排列在Yee元胞的棱和面上。三维Yee网格中每个电场分量沿其方向偏移半个网格磁场分量在垂直的空间对角位置偏移半个网格。三维迭代公式很长但结构一致。以Ex为例% 三维FDTD中Ex分量的更新示意 ex(i, j, k) ex(i, j, k) - dt * (hz(i, j, k) - hz(i, j-1, k)) / dy; ex(i, j, k) ex(i, j, k) dt * (hy(i, j, k) - hy(i, j, k-1)) / dz;每个分量更新都涉及两个相邻磁场分量的差分整个循环需要遍历i、j、k三个维度。在Matlab里如果直接写三重for循环运行速度会慢得让人怀疑人生。常见优化是把内层操作向量化每次更新一个二维切片% 向量化更新每次更新整个x方向的切片 ex(:, 2:end, :) ex(:, 2:end, :) - dt/dy * (hz(:, 2:end, :) - hz(:, 1:end-1, :));这种方式在网格规模不大时效率提升明显但缺点是代码可读性下降出了错不好查。我个人的习惯是一维和二维用向量化写法提升速度三维则优先保证逻辑清晰因为你花在调试上的时间往往比多跑几个小时更值钱。4.2 内存估算为什么100×100×100网格就让人头疼三维FDTD的内存消耗按网格数的三次方增长这是它最让人头疼的地方。六个场分量加上PML辅助变量每个变量都要存储一个和原网格一样大的三维数组。以100×100×100网格为例数据项变量数双精度内存约电场分量33 × 8MB 24MB磁场分量33 × 8MB 24MBPML辅助变量通常4~6个4 × 8MB 32MB合计10约80MB看着好像不多但这只是100³网格。到了200×200×200内存直接跳到640MB300×300×300就是2.2GB。实际工程里三维模型很少是简单立方体结构细节会让网格数轻松破百万内存消耗就奔着几十GB去了。所以三维FDTD的第一个准则是能用二维模拟就不要上三维。如果一个结构在一个方向上没有变化比如无限长圆柱、均匀波导截面完全可以降维成二维模型结果基本一致计算量却少了几个数量级。很多初学者一上来就想做三维全波仿真结果笔记本内存直接报错这是完全可以避免的。第二个准则是使用单精度。Matlab默认的double类型精度对FDTD绝大多数场景来说其实过剩了场值的精度主要受限于数值色散而非浮点误差。用single类型存储场分量内存减半速度经常还有提升% 用单精度减少内存占用 ez single(zeros(nx, ny, nz)); hy single(zeros(nx, ny, nz));注意源加载和观察点输出时可能需要转回double避免精度损失。4.3 三维FDTD的降维时机与混合策略实战中真正需要完整三维FDTD的场景其实有限。微带天线、贴片天线这类薄结构厚度方向的场变化很小可以在厚度方向用粗网格平面方向用细网格。波导结构如果长度远大于横截面可以只仿真一个截面再用模式展开。如果结构具备旋转对称性还可以用BOR体旋转对称方法把三维问题降成二维计算加一个角向模式展开。我见过最直接的做法是用二维FDTD先扫参锁定结构的大致尺寸范围再对最关键的几个尺寸做一次三维验证。这样既控制了计算量又保留三维结果的完整信息。这种“二维扫趋势、三维验终值”的工作流比一上来就跑三维省下太多时间。如果你真的需要大规模三维仿真Matlab脚本语言就力不从心了此时应该转向C或CUDA实现甚至考虑GPU加速。标题里的fdtd就是这个方向——C写的高性能FDTD库。Matlab版本适合学习和中小规模验证大规模生产型仿真还是交给编译型语言更合适。5. 避坑从一维到三维最常见的五个翻车场景5.1 仿真爆炸电场值疯狂增长现象跑了没几步场值变成NaN或者正负无穷大波形完全失控。原因几乎都是CFL条件不满足。dt取值超过稳定性极限数值误差随时间步进指数放大。解决先检查dt是否满足Courant条件。三维情形下dt 1 / (c * sqrt(1/dx² 1/dy² 1/dz²))很多新手只按一维公式dt dx/c来算在二维三维就直接炸了。另外一个隐藏坑介质区域中波的传播速度是c/sqrt(εr)比真空中慢有人把介质区域的dx改小了但忘记同步改小dt照样发散。5.2 二维结果和一维对不上现象同一个高斯脉冲在一维和二维中沿着传播方向观察波形形状不同。原因二维点源产生的波前是圆柱波幅度随距离衰减不像一维平面波幅度不变。另外二维代码的场分量顺序如果从TM换成TE观察点采的场分量不同波形自然不同。解决对比时先用同一条直线上的场分布且把二维幅度乘以sqrt(r)再做归一化补偿圆柱波的几何衰减。如果幅度曲线形态还是对不上检查Hx和Hy的更新符号是否搞反。5.3 PML反射残留观察点波形出现二次回波现象主脉冲过去之后波形尾部出现一个明显的起伏像是波被弹回来了。原因PML层数太少或sigma_max取值不当导致PML没有完全吸收入射波。另一个常见原因是PML和计算域之间的介电常数或电导率突变。解决把PML层数从10增加到15或20sigma_max从0.8调到1.0观察尾部反射是否变小。如果反射没有随PML变强而降低说明PML实现有逻辑错误而不是参数问题。此时先把结构简化成均匀介质只保留PML跑点源看反射排除结构散射的干扰。5.4 Matlab内存不足或运行极慢现象三维FDTD代码一运行就报内存不足或者二维代码跑了一个小时还没完。原因三维数组默认double类型内存翻倍代码里在循环中不断对数组进行重新分配每一次都会复制整个数组。解决第一步把所有场变量声明为single类型。第二步确保所有数组在循环前一次性分配好不要在循环内用扩展赋值比如a(end1) x这种写法。第三步检查是否需要三维仿真——很多问题用二维就能回答这是最有效的优化。5.5 结果和解析解偏差大现象仿真的反射系数或透射系数和理论值差了10%以上。原因网格太粗数值色散严重或者源加载有问题比如高斯脉冲的幅度归一化不对导致能量统计失真。解决把dx减半重新跑一次看结果是否明显改变。如果改变超过几个百分点说明收敛性还不行需要加密网格。此外在频域做后处理时要确认时域信号时长足够长避免截断导致的频谱泄漏。通常取Nt使波传播到计算域边缘并完全被PML吸收再截取场值做FFT。6. 验证你的FDTD结果从网格收敛到能量守恒6.1 网格细化测试检验离散误差的基本方法FDTD结果对不对最直接的办法是网格收敛测试。把dx减半同时dt按CFL关系减半重跑一次比较关键观察点处的波形。如果两次结果几乎重叠说明网格已经细化到误差可以忽略如果结果明显变化说明网格还不够细当前结果不能被信任。% 网格收敛性验证比较不同dx下观察点波形 dx1 0.02; dx2 0.01; % 两组网格尺度 % 分别运行FDTD得到ez1和ez2 % 计算两组结果的最大差异 diff_max max(abs(ez1 - ez2(1:2:end))); disp([网格细化后的最大差异: , num2str(diff_max)]);这个做法计算量不低但值得在关键仿真上做一次。通常做一组粗网格和一组细网格的对比就够了不需要做完整的三组收敛率分析。6.2 能量守恒监控发现隐藏反射的哨兵另一个低成本验证手段是监控计算域电场总能量。在没有源持续激励的情况下总能量应该随时间单调递减——波不断被PML吸收。如果能量曲线出现平段甚至上升说明PML有反射或计算域内部有非物理的增益。% 每次迭代后计算总能量 total_energy(n) sum(sum(sum(ez.^2 ex.^2 ey.^2)));这段代码在二维里对应sum(sum(ez.^2))三维再加两个电场分量。能量曲线平滑下降就是好结果如果你看到曲线尾部翘起来优先查PML参数。6.3 我自己的验证习惯每次写完三维FDTD代码我习惯先在一个一维测试用例上验证整个流程网格收敛测试、能量监控、解析解对比三个都过了才开始正式仿真。这个习惯帮我挡掉过很多次debug问题——常常是电磁方程组本身的符号问题而不是PML或边界的问题。先跑一维因为你一眼能看出波形对不对确认底层迭代无误再升到二维和三维。希望帮到你。本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
