FDTD时域有限差分这名字听起来挺唬人但把它拆开看有限差分就是把麦克斯韦方程组里那些偏微分换成加减乘除时域就是一格一格往前推进时间。用 MATLAB 写一个基于 FDTD 的二维电磁波传播模拟本质上就是搭一个数值“试验台”观察波在自由空间中怎么冒出来、往外跑、最后被边界“吞”掉。我这次做的这个小项目正是把整个过程用最直接的方式跑通了TMz 模式、二维自由空间、点源激励、正弦脉冲边界用 Mur 一阶 ABC 吸收可视化用 MATLAB 的 pcolor 逐帧刷新。项目本身不复杂但把“原理—离散—边界—源—可视化”这条链路完整走一遍之后你对 FDTD 的理解会和只看公式完全不一样。想做电磁仿真入门、需要演示电磁波传播场景或者想在课程设计里加一个有物理内涵的 MATLAB 程序的同学这篇文章应该能帮你省不少时间。不需要什么高端硬件一台普通笔记本装好 MATLAB按下面的代码跑一遍就能看到波从中心扩散、接触边界、被“吸收”的全过程。1. 项目概述——FDTD 自由空间模拟到底在做什么1.1 技术选型背后的逻辑先看模式选择。二维问题里最简单也最常用的是 TMz 模式电场只有一个垂直分量 (E_z)磁场有两个水平分量 (H_x)、(H_y)。好处是场方程只有三个标量方程不会引入矢量场的复杂耦合出图也直观。三维 FDTD 的实现差异主要是分量多了核心思想完全一致。所以用 TMz 入门性价比非常高。再看边界条件。标题里写的 ABC 是吸收边界条件我选用的是 Mur 一阶公式。为什么不用理想导体边界因为它会把出射波反射回来模拟出来的就不是自由空间而是谐振腔。为什么不用更先进的 PML因为 Mur 一阶只有一种经典的单向波公式代码量小、物理概念清楚作为学习 FDTD 的边界处理是非常合适的起点。我没让问题一开始就变得复杂这也是很多初学者最容易犯的毛病。1.2 正弦脉冲激励源的特点源选正弦脉冲不是说随便挑一个正则的 sin 函数。它的本质是一个有限持续时间的窄带信号载波频率定为 1 GHz持续 5 个周期5 ns 后源自动关闭。相比高斯脉冲的宽频特性正弦脉冲的主要成分集中在单一频率附近非常方便观察单频波在自由空间的传播和相速度表现。做演示时你能清清楚楚看到一圈一圈的等相位面往外扩散这种视觉反馈对理解 FDTD 很有帮助。这套组合方案的关键优势是“每一块都能独立验证”笛卡尔网格下的场更新是标准教科书内容Mur 边界是标准教科书内容正弦脉冲点源也是标准教科书内容。三个标准内容拼在一起就能组装出一个完整的开放边界电磁波传播演示系统。后续想加介质、加散射体、加损耗只需要在局部格点上改参数不用重构主循环。2. FDTD 原理麦克斯韦方程是怎么变成递推公式的2.1 从旋度方程到三项更新式自由空间中麦克斯韦方程的时域形式可以写作[ \nabla \times E -\mu_0 \frac{\partial H}{\partial t} ][ \nabla \times H \varepsilon_0 \frac{\partial E}{\partial t} ]在二维 TMz 模式下电场只保留 (E_z)磁场只保留 (H_x)、(H_y)展开后得到三个标量方程[ \frac{\partial H_x}{\partial t} -\frac{1}{\mu_0} \cdot \frac{\partial E_z}{\partial y} ][ \frac{\partial H_y}{\partial t} \frac{1}{\mu_0} \cdot \frac{\partial E_z}{\partial x} ][ \frac{\partial E_z}{\partial t} \frac{1}{\varepsilon_0} \cdot \left( \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} \right) ]剩下的事情就是把偏导数用差分替代。这里需要引入 Yee 网格电场 (E_z) 放在整数网格点上磁场 (H_x) 放在 y 方向半格处磁场 (H_y) 放在 x 方向半格处时间上也错开半步。这样做的好处是空间差分立刻变成相邻两点的差时间递推则是标准的蛙跳格式。用索引写出来每一步就非常机械了。假设 (E_z(i,j)) 表示整数点(H_x(i,j)) 表示 ((i, j1/2)) 的位置(H_y(i,j)) 表示 ((i1/2, j)) 的位置那么在一个时间步 (dt) 里三个更新式分别是[ H_x^{n1/2}(i,j) H_x^{n-1/2}(i,j) - \frac{\Delta t}{\mu_0 \Delta y} \left[ E_z^n(i,j1) - E_z^n(i,j) \right] ][ H_y^{n1/2}(i,j) H_y^{n-1/2}(i,j) \frac{\Delta t}{\mu_0 \Delta x} \left[ E_z^n(i1,j) - E_z^n(i,j) \right] ][ E_z^{n1}(i,j) E_z^n(i,j) \frac{\Delta t}{\varepsilon_0} \left[ \frac{H_y^{n1/2}(i,j) - H_y^{n1/2}(i-1,j)}{\Delta x} - \frac{H_x^{n1/2}(i,j) - H_x^{n1/2}(i,j-1)}{\Delta y} \right] ]注意公式里的正负号这是最容易出错的地方。(H_x) 的更新带负号(H_y) 带正号(E_z) 内部的旋度项是 (H_y) 对 x 的偏导减去 (H_x) 对 y 的偏导。很多初学者把符号搞反结果波不但不往外传播反而原地发散看上去就像系统爆炸。2.2 稳定性条件为什么时间步和空间步必须配好FDTD 是显式格式存在数值稳定性上限。二维均匀网格的 Courant 条件给的是[ \Delta t \le \frac{1}{c_0 \sqrt{\frac{1}{\Delta x^2} \frac{1}{\Delta y^2}}} ]当 (\Delta x \Delta y d) 时简化成[ \Delta t \le \frac{d}{c_0 \sqrt{2}} ]我的代码里直接取等号即dt dx / (c0 * sqrt(2))。这是很多教材推荐的做法等于完全利用稳定边界实际工程里有人会留 10% 余量取 0.9 倍避免数值色散积累。除稳定性之外空间分辨率也必须考虑。一个波长至少要有 10 到 20 个网格点否则波的形状会严重失真。我这里 (f_0 1\text{ GHz})波长 (0.3\text{ m})而 (dx 0.01\text{ m})每个波长有 30 个网格点绰绰有余。如果频率提高十倍(dx) 就要缩小十倍时间步也跟着变小计算量会指数增加。这就是高频模拟为什么贵的根本原因。2.3 一阶 Mur 边界条件的物理意义再回来看 ABC。所谓吸收边界本质是在截断边界上构造一个“只让波出去、不让波回来”的数值条件。Mur 一阶假设边界上的外向波满足单向波动方程比如左边界 (x0) 处有[ \frac{\partial u}{\partial t} - c_0 \frac{\partial u}{\partial x} 0 ]离散化后可以解出边界上的新值。设 (s c_0 \cdot \Delta t / \Delta x)令 (r (s-1)/(s1))则左边界更新为[ E_z^{n1}(1,j) E_z^{n}(2,j) r \cdot \left[ E_z^{n1}(2,j) - E_z^{n}(1,j) \right] ]右边界、下边界、上边界形式类似。系数 (r) 在这里是负的因为 (s 0.707 1)所以 (r \approx -0.1716)。边界更新本质上不是简单的“把边界清零”而是根据内部的新值和旧值外推让波近似“走出去”。必须承认一阶 Mur 对正入射波吸收效果最好斜入射时会有一定反射。后文我会专门讲怎么识别和评估这种反射。它毕竟不是 PML不能指望在所有入射角下都完美无缺。3. MATLAB 代码实现与关键细节3.1 参数设置每一格该怎么选写代码的第一步是定义物理参数和网格。代码里我给出了一组可以直接跑的默认参数但你要清楚每个数的含义改的时候才不会乱。% % FDTD 2D TMz 模拟自由空间电磁波传播 % 吸收边界Mur 一阶 ABC % 激励源正弦脉冲点源 % clear; clc; close all; %% 物理常数与网格参数 c0 3e8; % 光速 m/s eps0 8.854e-12; % 真空介电常数 F/m mu0 pi*4e-7; % 真空磁导率 H/m dx 0.01; % 空间步长 x (m) dy dx; % 空间步长 y (m) dt dx / (c0 * sqrt(2)); % 时间步长 (Courant 极限) Nx 300; % x 方向网格数量 Ny 300; % y 方向网格数量 Nt 900; % 总时间步数 % 物理坐标矢量 x_vec (0:Nx-1) * dx; y_vec (0:Ny-1) * dy; %% 场分量初始化 Ez zeros(Nx, Ny); % 电场 z 分量 Hx zeros(Nx, Ny); % 磁场 x 分量 Hy zeros(Nx, Ny); % 磁场 y 分量 %% 激励源参数 f0 1e9; % 载波频率 1 GHz T0 1 / f0; % 周期 ncycle 5; % 脉冲持续 5 个周期 t_source ncycle * T0; % 源持续时间 % 源位置取网格几何中心 src_x round(Nx / 2); src_y round(Ny / 2); %% Mur ABC 系数 s c0 * dt / dx; % 当前 Courant 数 1/sqrt(2) r (s - 1) / (s 1); % Mur 一阶边界系数c0、eps0、mu0是真空常数dx决定空间分辨率dt由 Courant 条件推出Nx、Ny决定计算区域大小Nt是总时间步数。源参数f0、ncycle、t_source定义了激励信号。src_x、src_y是源的位置我放在网格中心这样图看起来最对称。注意一个映射细节数组索引和物理坐标的对应关系是 (x (i-1)\Delta x)(y (j-1)\Delta y)。如果你把索引直接当坐标等画图时标注的横轴就会差一个 (dx)虽然整体分布不变但坐标刻度就是错的。这一点在分析结果时很关键。3.2 场更新主循环里每一步到底发生了什么先看主循环里的 Hx、Hy、Ez 三个更新式以及源的注入位置。%% 主时间步循环 figure(Color, w, Position, [80 80 640 540]); for n 1:Nt t_advance (n - 0.5) * dt; % 半时间步时刻 % --- 更新 Hx位置 (i, j1/2) --- Hx(:, 1:end-1) Hx(:, 1:end-1) - (dt / (mu0 * dy)) * ... (Ez(:, 2:end) - Ez(:, 1:end-1)); % --- 更新 Hy位置 (i1/2, j) --- Hy(1:end-1, :) Hy(1:end-1, :) (dt / (mu0 * dx)) * ... (Ez(2:end, :) - Ez(1:end-1, :)); % --- 保存更新前的 EzMur 边界需要上一时刻边界值 --- Ez_old Ez; % --- 更新 Ez 内部区域 --- Ez(2:end-1, 2:end-1) Ez(2:end-1, 2:end-1) (dt / eps0) * ... ( (Hy(2:end-1, 2:end-1) - Hy(1:end-2, 2:end-1)) / dx - ... (Hx(2:end-1, 2:end-1) - Hx(2:end-1, 1:end-2)) / dy ); % --- 注入正弦脉冲点源加性源 --- if t_advance t_source Ez(src_x, src_y) Ez(src_x, src_y) sin(2 * pi * f0 * t_advance); end % --- Mur 一阶 ABC 四边更新 --- Ez(1, :) Ez_old(2, :) r * (Ez(2, :) - Ez_old(1, :)); Ez(end, :) Ez_old(end-1, :) r * (Ez(end-1, :) - Ez_old(end, :)); Ez(:, 1) Ez_old(:, 2) r * (Ez(:, 2) - Ez_old(:, 1)); Ez(:, end) Ez_old(:, end-1) r * (Ez(:, end-1) - Ez_old(:, end)); % --- 每隔 20 步绘制快照 --- if mod(n, 20) 0 pcolor(x_vec, y_vec, Ez); shading interp; caxis([-0.3 0.3]); colormap(jet); colorbar; xlabel(x / m); ylabel(y / m); title(sprintf(Ez 分布 t %.2f ns, t_advance * 1e9)); drawnow; end endHx、Hy、Ez三个更新式要按顺序执行先更新磁场再更新电场因为蛙跳格式中电场需要用刚刚更新过的磁场。这是我这里的具体顺序如果你以后用其他库或者换语言也要保持这种时间交错的思路。更新完内部 Ez 后立刻注入激励源。注意源是在整个场更新之后才加入的这是显式 FDTD 添加点源的常见做法之一。实际物理上源也可以视为某个格点上的场值被外部驱动放在更新后不会破坏差分格式的时间一致性。这里有个细节我没有用“硬源”而是直接把正弦函数加到Ez(src_x, src_y)上这叫加性源。硬源会直接覆盖场值入射波和散射波在源位置无法叠加加性源则允许场自由通过源点更接近物理真实。提示场更新顺序一定要是先磁场后电场。如果你先更新 Ez 再更新 Hx、Hy整个蛙跳时间结构就塌了模拟结果会偏得离谱。3.3 Mur ABC四边边界更新为什么要在电场更新之后做需要特别强调的是Mur 边界的实现顺序不能乱。我用Ez_old保存了电场更新前的值然后更新内部 Ez紧接着在边界上用Ez_old和更新后的内部值算出边界新值。为什么必须这样因为 Mur 公式同时涉及上一时刻的边界值和当前时刻的内部值。如果你先更新边界再更新内部边界上新值会被内部更新覆盖如果你不保存旧值直接拿更新后的内部值当“旧值”用公式就错了。左边界和右边界关注 x 方向下边界和上边界关注 y 方向。四个角落的计算可能被两条边重复更新我的做法是让后执行的边界决定角点值。严格来说角点应单独用一维外向波方程处理但一阶方案下误差其实很小多数演示场景可以忽略。边界更新的四条语句里索引变化很有规律左边界取出原数组第二行作为旧内部值用更新后的第二行做外推边界行是第一行右边界用倒数第二行外推边界行边界行是最后一行下边界用第二列外推边界列边界列是第一列上边界用倒数第二列外推边界列边界列是最后一列。写错索引最常见的后果是边界上出现镜像条纹或者波一碰到边界就产生剧烈振荡。如果看到这种症状先把这四条语句的索引逐行打印出来对照物理坐标图看一遍。3.4 动态可视化需要注意什么% 循环结束后的最终静态图 figure; subplot(1,2,1); pcolor(x_vec, y_vec, Ez); shading interp; colorbar; axis equal tight; title(Ez 空间分布); xlabel(x / m); ylabel(y / m); subplot(1,2,2); plot(x_vec, Ez(:, src_y), b-, LineWidth, 1.2); grid on; xlabel(x / m); ylabel(Ez); title(沿 y 中心线的切片);pcolor画的是伪色彩图注意转置因为Ez是按Nx * Ny存储的x 对应第一维y 对应第二维pcolor的输入是 x、y、Z其中 Z 是Ny * Nx所以要写Ez。我每 20 步刷新一次否则每步都drawnow会让整个循环变得很慢。颜色范围固定下来不要让它自动缩放不然你会看到每一步的色标都在跳根本没法观察波前的衰减。如果你觉得逐帧动画还不过瘾可以在循环内把Ez存成三维数组也可以在drawnow同时用getframe录制成视频。不过这会让内存占用骤增我的建议是先用快照方式观察确认现象正常后再考虑录动画。4. 模拟结果与现象解读4.1 波前扩散的完整过程如果把前 200 个时间步的数据保存下来回放看到的应该是源点首先出现一个逐渐增大的正负交替环形波前波峰和波谷像同心圆一样向外扩散。因为是单频正弦脉冲波前是一圈一圈的而在脉冲结束后源位置不再有新信号外围波前继续向前传播。计算一下时间尺度(dt \approx 23.6\text{ ps})200 步约 (4.7\text{ ns})。此时 1 GHz 波完成了大约 4.7 个周期加上前导波和边界之间还有距离波没有碰到边界。做到大约 400500 步波前开始接触计算区域边缘这时最能检验 ABC 的效果。需要提醒的是正弦脉冲虽然是窄带信号但它的瞬时开通仍然会带来一个“前导波头”。这个波头从零突然跳到正弦的某个相位空间上表现为一个较陡的过渡带它比纯正弦波更容易产生数值色散。如果你把dx调大会明显看到波头变得圆钝、拖尾变长这就是数值色散的直接证据。4.2 边界反射怎么判断判断 ABC 好不好有一个很直接的观察法波前碰到边界后边界内侧不应该出现明显的“二次波纹”再往回流。如果看到一圈清晰的反射波从边界折回中心那说明边界条件没有正常工作或者实现有误。我测试时用的是 300×300 的网格900 个时间步足够让波到达边界并大部分吸收。尾部还有一点残余反射是正常现象——一阶 Mur 对斜入射本来就有局限圆形波前到达矩形边界时除正对源点的区域外大部分都是斜入射。后文我会专门讲如何进一步削弱反射。另一个容易误判的点波的柱面扩散本身就会让振幅随距离衰减。二维自由空间中理想柱面波的电场幅度按 (1/\sqrt{R}) 衰减。所以看到靠近边界的波峰比中心位置低并不代表能量被吸收了而是几何衰减的必然结果。很多人第一次看到这个现象会误以为边界有问题先记下这个物理规律能省去不少排查时间。4.3 观察相速度和频率的关系正弦单频波在均匀网格里传播时数值相速度和真实相速度会有细微差别这叫数值色散。二维 FDTD 的色散关系比较长这里我不展开公式只讲观察方法你可以在源点附近某条固定直径方向放几个探针点记录 (E_z) 随时间变化然后求出相邻探针点之间的相位差反推数值相速度。通常情况下每波长 30 格时数值相速度和 (c_0) 的偏差在 1% 以内肉眼看图基本分辨不出来。每波长低于 10 格时偏差会明显到让波前形状都开始扭曲。5. 参数敏感性与常见问题排查实录5.1 空间步长对精度的影响如果dx等于波长的 1/5数值色散会极其明显波的相速度不再等于 (c_0)而且不同方向传播的相位误差不一样原本的圆波前会变成多边形。有人把这种现象叫“各向异性色散”。如果你看到波前变形先查空间分辨率每波长至少保证 1020 个网格点。这里有个经验值想看清干涉结构或驻波特征时每波长 20 格以上比较稳妥如果只是演示波传播趋势10 格也能看。但千万别认为 10 格就“够用”了因为 FDTD 的误差是累积的跑得越久色散导致的相位误差越明显。5.2 Courant 条件不满足的典型症状把dt取成dx / c0不带根号二就会违反二维 Courant 条件。短时间看不出问题几十步之后场值开始指数膨胀图像出现高频闪烁最后直接变成 NaN。这不是边界或者源的问题而是时间步太大。遇到这种情况立刻回头检查dt的计算公式。在排查时可以把dt再缩小 10%观察是否稳定。如果缩小后数值稳定、波形正常说明就是 Courant 条件踩线了。显式格式的一个特点是“要么稳定要么爆掉”不会给你中间过渡态。5.3 Mur 边界总是有一点反射怎么办首先确认公式和更新顺序没有错。如果确认无误还想进一步减少反射可以升级到 Mur 二阶或者换上 PML。PML 是匹配层思想通过在边界区引入人工各向异性介质吸收波效果比 Mur 好很多但实现代码要多好几行。我建议先把 Mur 跑通再去看 PML这样你能真正理解为什么要做匹配层。Mur 二阶比一阶多了一项关于切向二阶导数的修正边界附近需要额外存储一行一列的数据来构造二阶差分。这个升级不复杂但代码美观度会下降。如果你只是做自由空间演示一阶 Mur 已经足够如果要研究散射体那确实应该考虑 PML。5.4 常见问题速查表现象可能原因排查思路场值几个时间步就发散dt 超过 Courant 条件检查 dt 是否满足二维不等式波前明显变形不是圆dx 太大数值色散严重加密网格保证每波长 ≥ 20 格边界出现强反射回流Mur 实现顺序或符号错重新核对四条边的索引和 Ez_old 使用位置源点不出现波源位置数组下标越界确认 src_x、src_y 在 2 到 Nx-1 内动画太卡每步都 drawnow每隔 20 或 50 步刷新一次图像上下翻转忘记转置 Ezpcolor 里用 Ez因为存储顺序是 Nx×Ny波碰到边界后明显炸开边界更新被内部更新覆盖确认 Mur 更新在 Ez 内部更新之后执行5.5 调试用的三个小技巧把Nx、Ny先设成 50×50 的小网格跑通再放大。加一个探针点记录Ez随时间变化直接观察波形。把边界换成 PEC 跑一遍对照组再和 ABC 结果对比立刻明白吸收边界的价值。这是我实际调试时最常用的三个手段。具体来说PEC 对照组就是把 Mur 边界那四条更新语句注释掉让边界值保持为 0。跑同样的源和时间步你会看到波碰到边界后产生非常干净的镜面反射反射波和原波相互叠加最后形成复杂的驻波图案。再看 ABC 版本几乎滚不出来反射波这么一对比边界条件的贡献就一目了然了。还有个小技巧检查边界吸收率时不一定要用完整的大网格。可以做一个窄长条网格比如 60×200把源放在一端观察波沿线传播到另一端后的剩余幅度。这样能清楚地看到一次边界反射的幅度比看二维平面图更直观。6. 进一步扩展方向这个项目做完之后往几个方向扩展都很顺手把源从单点改成两个同相或反相点源就能看到干涉条纹把局部介电常数改成一组数值就成了介质柱散射模拟把二维升级成三维核心代码结构基本不变只是多出几个场分量。我个人觉得最划算的扩展是加一个简单的介质圆盘在网格区域内把某个圆内的介电常数改成 2.25水或者玻璃的典型值然后重新跑同样源。你会看到波在圆盘边界发生部分反射和透射透射波在圆盘内部波长被压短因为介质中波长 (\lambda \lambda_0 / \sqrt{\varepsilon_r})。这个现象直接用肉眼就能验证对理解介质中的电磁波非常有帮助。我实际跑下来最大的体会是FDTD 的参数选择和实现顺序同等重要。写代码之前先把时间步、空间步、波长关系算清楚比写完再调参要高效得多。你如果想自己动手试建议第一步就把Nx、Ny降到 100 以下把绘图间隔调大专注观察波前和边界反射然后再慢慢加大规模。跑通一次之后你会发现电磁场模拟并不是什么神秘的东西无非就是把麦克斯韦方程“翻译”成循环和数组操作剩下的事就是让计算机按部就班地算下去。
