三维FDTD电磁仿真：Yee网格、PEC边界与Python动画实现

简介基于Yee算法的三维麦克斯韦方程组数值模拟项目面向电磁场方向的学习者与科研入门者用于直观解析电磁波传播、反射与散射过程。压缩包共5个m文件体积仅6KB全部为MATLAB脚本以主程序为核心配合全局变量配置和多个PEC边界条件输出子程序整体结构紧凑适合边看代码边理解有限差分法的完整实现流程。项目采用PEC理想导体边界条件适用于波导、天线等封闭或无限大导体结构的建模模拟。目前已有220人学习运行代码后可生成电磁场在三维空间随时间的四维动画动态展示电场与磁场交替更新、边界反射等关键现象。通过分析这些代码可掌握麦克斯韦方程组在三维网格中的数值离散方法并理解PEC边界条件的程序实现进一步为天线设计、雷达系统、通信工程等领域的电磁仿真打下实用基础。1. YEE3D 解决什么问题Yee 网格、PEC 与三维电磁场动画的最小闭环第一次跑三维 FDTD 程序最让人头疼的往往不是方程本身而是六个场分量铺开后根本看不出电磁波在往哪个方向走。YEE3D.zip 这个名字里的三个关键词串起了一条完整链路用 Yee 网格把 Maxwell 旋度方程离散成可迭代的差分格式用 PEC 边界把计算域围住再把每一步的电磁场分布输出成动画帧。这里的三维分别体现在 Ex/Ey/Ez 与 Hx/Hy/Hz 在立体网格上的交替更新PEC 则决定了边界上哪一个切向分量必须为零。这篇文章从一个可运行的 30×30×30 主循环出发讲更新公式怎么落成切片代码、PEC 在哪一步生效、动画参数怎么调以及用什么标准确认结果没有算错。适合准备写电磁场可视化、或想用 FDTD 验证 PEC 边界处理是否正确的人。2. 三维 Yee 网格与麦克斯韦方程更新公式怎么变成代码三维数值模拟里Yee 网格解决的是一个很实际的离散问题电场旋度给出磁场磁场旋度又给出电场两者互相依赖同一个空间点上无法同时定义所有分量。Yee 的做法是把电场放在棱边中点、磁场放在面心两者错开半个网格这样每个旋度分量的差分都恰好作用在两个相邻场分量的中点上空间二阶精度天然成立。时间上电场和磁场再错开半步更新构成蛙跳格式这是三维 FDTD 稳定迭代的基础。下面把六条更新公式变成不加三重循环的 Python 切片代码。2.1 三维 FDTD 主循环的最小实现六个切片更新场分量均匀立方网格 dxdydzdl 时Hx 的离散形式是Hx[i,j,k] dt/mu * ((Ey[i,j,k1] - Ey[i,j,k]) - (Ez[i,j1,k] - Ez[i,j,k])) / dlEy 沿 z 方向的差对应 ∂Ey/∂zEz 沿 y 方向的差对应 ∂Ez/∂y两者相减就是旋度的 x 分量。Hy、Hz 和三个电场分量只是交换了坐标轴写成数组切片后整个程序只需要六个更新语句import numpy as np from math import sqrt N 30 c0 3.0e8 dx 1.0e-3 # 空间步长 1 mm30 格共 3 cm dt dx / (c0 * sqrt(3)) # CFL 条件取等号 mu0 4 * np.pi * 1e-7 eps0 8.854e-12 cH dt / (mu0 * dx) cE dt / (eps0 * dx) Ex np.zeros((N, N, N)); Ey np.zeros((N, N, N)); Ez np.zeros((N, N, N)) Hx np.zeros((N, N, N)); Hy np.zeros((N, N, N)); Hz np.zeros((N, N, N)) def update_H(): Hx[:-1, :-1, :-1] cH * ( (Ey[:-1, :-1, 1:] - Ey[:-1, :-1, :-1]) - (Ez[:-1, 1:, :-1] - Ez[:-1, :-1, :-1]) ) Hy[:-1, :-1, :-1] cH * ( (Ez[1:, :-1, :-1] - Ez[:-1, :-1, :-1]) - (Ex[:-1, :-1, 1:] - Ex[:-1, :-1, :-1]) ) Hz[:-1, :-1, :-1] cH * ( (Ex[:-1, 1:, :-1] - Ex[:-1, :-1, :-1]) - (Ey[1:, :-1, :-1] - Ey[:-1, :-1, :-1]) ) def update_E(): Ex[1:-1, 1:-1, 1:-1] cE * ( (Hz[1:-1, 1:-1, 1:-1] - Hz[1:-1, :-2, 1:-1]) - (Hy[1:-1, 1:-1, 1:-1] - Hy[1:-1, 1:-1, :-2]) ) Ey[1:-1, 1:-1, 1:-1] cE * ( (Hx[1:-1, 1:-1, 1:-1] - Hx[1:-1, 1:-1, :-2]) - (Hz[1:-1, 1:-1, 1:-1] - Hz[:-2, 1:-1, 1:-1]) ) Ez[1:-1, 1:-1, 1:-1] cE * ( (Hy[1:-1, 1:-1, 1:-1] - Hy[:-2, 1:-1, 1:-1]) - (Hx[1:-1, 1:-1, 1:-1] - Hx[1:-1, :-2, 1:-1]) )H 场更新索引范围是 0 到 N-2E 场更新索引范围是 1 到 N-2两个范围在空间上错开正是 Yee 网格半格偏移在数组层面的体现。所有差分都来自相邻数组切片的差没有额外插值每一条语句都和一条旋度分量直接对应。cH、cE 来自同一个 dt只要 dt 满足 CFL 条件迭代就是稳定的。主循环顺序是经典的蛙跳先更新 H再更新 E然后强制 PEC 边界最后注入源。def apply_pec(): Ey[0, :, :] 0.0; Ez[0, :, :] 0.0 Ey[-1, :, :] 0.0; Ez[-1, :, :] 0.0 Ex[:, 0, :] 0.0; Ez[:, 0, :] 0.0 Ex[:, -1, :] 0.0; Ez[:, -1, :] 0.0 Ex[:, :, 0] 0.0; Ey[:, :, 0] 0.0 Ex[:, :, -1] 0.0; Ey[:, :, -1] 0.0 for step in range(300): update_H() update_E() apply_pec() # 2.2 节注入偶极子源apply_pec 的六个面赋值对应 PEC 边界的实际含义边界上的切向电场必须为零。以 x 方向两个面为例切向分量是 Ey 和 Ez法向 Ex 不在置零范围y 面和 z 面同理。把这段代码放进时间循环里而不是只在初始化时执行一次是 PEC 最容易踩中的第一个坑。场分量Yee 网格位置切片索引范围Ex(i1/2, j, k)[1:-1]Ey(i, j1/2, k)[1:-1]Ez(i, j, k1/2)[1:-1]Hx(i, j1/2, k1/2)[:-1]Hy(i1/2, j, k1/2)[:-1]Hz(i1/2, j1/2, k)[:-1]2.2 激励源怎么加偶极子电流与脉冲参数没有源FDTD 迭代多少步都是零解。z 方向偶极子源的常见做法是每步在中心点的 Ez 上叠加一个高斯调制的正弦脉冲t0 40.0 * dt # 脉冲峰值时刻 tau 10.0 * dt # 脉冲宽度 f0 c0 / (8.0 * dx) # 中心频率对应波长 8 个网格 def source_value(step): t step * dt envelope np.exp(-((t - t0) / tau) ** 2) return envelope * np.sin(2 * np.pi * f0 * t)f0 取“8 个网格一个波长”是两个约束折中的结果最少要保证 6 格数值色散才不会把波前打散超过 20 格时计算量按立方上涨动画演示没有必要。tau 决定脉冲频谱宽度tau 越大带宽越窄连续波激励时可以直接去掉高斯包络只保留正弦项。2.3 PEC 边界置零的位置为什么放在 E 更新之后apply_pec 必须放在 update_E 之后。如果先置零再更新边界节点又会被内部差分修正回来相当于 PEC 条件在时间上晚执行了一个半步反射系数会偏小。PEC 盒子的棱边是两个面的交界顶点是三个面的交界面置零操作在棱边和顶点上会有重复赋值这是正常现象不是 bug。真正要留意的是H 场的更新范围天然比 E 场窄一圈边界最外层的 H 分量没有参与更新这也和 PEC 切向电场为零的物理约束一致。3. PEC 边界与数值稳定性三个最容易翻车的参数三维模拟跑起来之后最常见的症状是“跑几十步就 NaN”和“角落出现尖峰”。两个问题都指向同一类东西时间步长不对、空间步长过粗、或边界执行顺序不对。把这三个参数逐一说清楚大部分稳定性问题都可以在主循环内部解决。3.1 CFL 条件与时间步长三维的 sqrt(3) 别拿错三维均匀网格的 CFL 上限是dt ≤ 1 / (c · sqrt(1/dx² 1/dy² 1/dz²))dxdydz 时就是 dt ≤ dx / (c · sqrt(3))。这里 sqrt(3) 是从二维升到三维最容易被忽略的变化二维公式是 sqrt(2)。取等号是理论允许的极限数值误差会累积到临界点工程上我一般乘 0.9courant 0.9 dt courant / (c0 * np.sqrt(1/dx**2 1/dy**2 1/dz**2))CFL 系数越大每物理秒需要的步数越少但超过极限会指数发散。反过来调太小也不行同样物理时间步数变多累积截断误差反而更大。0.8 到 0.9 是三维 FDTD 常见取值区间。3.2 空间步长与每波长网格数精度和计算量的平衡对给定的中心频率 f0空间步长按最短波长决定。常见做法是每波长不低于 10 个网格动画演示可以把标准放宽到 8 格快速验证允许 6 格。每波长网格数适用场景表现6快速验证、代码调试波前有可见畸变只能看趋势8动画演示波前形态基本正确10定量比较常用折中色散较难察觉20精确谐振频率计算计算量大适合最终校核在 YEE3D 这类教学演示里我优先保证能在一分钟内跑完连续波持续激励时由于波在 PEC 盒内不断反射叠加低频长波的驻波图样对网格数更敏感这时建议直接上 12 格以上。3.3 边界置零的覆盖范围与执行顺序切向、棱边与错位存储3.3.1 切向与法向棱边顶点为什么被重复赋值PEC 只约束切向电场法向电场可以存在。面置零代码里x 面清零 Ey/Ez、y 面清零 Ex/Ez、z 面清零 Ex/Ey正是切向与法向的区分。棱边是两个 PEC 面的交界顶点是三个面的交界同一个数组元素会被多条置零语句覆盖。因为操作都是赋零覆盖顺序没有影响读者看到重复赋值不需要担心。真正的问题是“只置零一次”PEC 条件必须在每个时间步强制执行否则内部点的差分会在下一步把边界值重新写回非零。3.3.2 执行顺序与错位存储PEC 条件必须在每个时间步生效严格按 Yee 网格存储时Ex 在 x 方向偏移半格电场外边界并不正好落在数组首层而是在数组之外的半格平面上。切片代码里看不到这层偏移但注释必须写清楚哪些层代表 PEC 面、哪些层是计算域内部。检查边界是否真正生效可以跑一个不含源的盒子初始给一列平面波看反射波幅值是否接近 1。用动画波前来判断不可靠反射系数差 0.1 在视觉上基本看不出来。4. 三维电磁场动画输出从数值结果到能看的帧主循环跑完Ex/Ey/Ez 是几十层三维数组动画只是把这些数组以正确的方式投射到二维屏幕上。三维场数据可以呈现的方式有三种取舍很不一样。4.1 三种可视化方式怎么选切面、矢量箭头与等值面正交切面是三维数值模拟最稳的选择。用 matplotlib 的 imshow 或 pcolormesh 显示某个平面上的一个场分量x 轴和 y 轴都是空间坐标颜色是场强视觉直接、开销小。矢量箭头 quiver 能把电场方向画出来但每帧几百个箭头会把画面塞满适合看局部方向而不是整体传播。等值面能呈现立体波形但三维渲染依赖项较重在演示程序里不如三个切面更可靠。我的做法是默认输出三个正交切面合成一张图分别显示 Ez 的 xy、xz、yz 截面。方式信息量依赖典型用途正交切面场强弱、波前位置matplotlib动画主图矢量箭头方向、漩涡matplotlib局部放大等值面立体波前形态mayavi/plotly封面图、演示4.2 用 FuncAnimation 把 Ez 切面存成 GIF 动画常见的做法是保存每一帧的 png再用 ffmpeg 合成 mp4如果想直接得到一个 GIF用 matplotlib 的 FuncAnimation 配合 pillow 一次完成。动画回调和 FDTD 主循环分开写回调只负责把数组放入已经创建好的 Image 对象import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation fig, ax plt.subplots(figsize(5, 5)) data0 Ez[:, :, N // 2] im ax.imshow(data0, cmapRdBu, vmin-0.5, vmax0.5, originlower, extent[0, N, 0, N]) def draw_frame(step): im.set_data(Ez[:, :, N // 2]) ax.set_title(fstep {step:04d}) return im, anim FuncAnimation(fig, draw_frame, framesrange(0, 300, 5), interval50, blitTrue) anim.save(yee3d_pec.gif, writerpillow, fps15)frames 用 range(0, 300, 5) 表示每 5 个时间步取一帧一共 60 帧interval 是播放帧间隔 50msfps15 控制保存后的播放速度。vmin/vmax 固定为 ±0.5不要用自动色标否则每帧的归一化范围不同动画会一帧亮一帧暗。创建 Image 对象后反复 set_data比每帧重新 pcolormesh 快一个数量级300 帧演示也不会卡。如果要把三个切面都放进去创建一个 1×3 的子图每个子图各建一个 im 对象回调里依次 set_data。改成 mp4 时把 writer 换成 ffmpeg 即可。4.3 动画里应该观察到的传播特征球面波、反射与驻波z 方向偶极子在 PEC 盒内的传播过程有很强的可预期性高斯脉冲前半段从中心向四周扩散波前应近似球面碰到盒壁后反射波开始向中心收拢与向外传播的波叠加在切面上形成一圈圈明暗交替的驻波条纹。xz 与 yz 切面在偶极子轴向上有对称分布xy 中心切面则应看到圆环。如果动画里出现明显的菱形波前通常是每波长网格数不足如果出现从边界向内扩散的直线条纹多半是 PEC 置零只执行了一次边界在“漏波”。5. 三维模拟结果验证的三个实用技巧程序能出动画不等于程序算对了。动画看不出的误差用三种定量检查可以在几行代码内暴露出来。5.1 用 PEC 腔体谐振频率反推网格与边界矩形 PEC 腔的理论谐振频率是fmnp (1 / (2 · sqrt(μ₀ε₀))) · sqrt((m/Lx)² (n/Ly)² (p/Lz)²)以 30 mm 立方腔为例TE101 模约 7.07 GHz。验证时把 2.2 节的正弦项去掉用纯高斯脉冲做宽带激励记录中心点的 Ez 后做 FFT峰值频率应落在理论值附近。偏差超过 5% 先查 Lx 与 N·dx 是否一致再查 apply_pec 是否在每个时间步执行。这里注意动画演示时用的 37.5 GHz 激励源不一定会显著激励 7 GHz 模式做谐振验证要换更宽的脉冲。5.2 用场能量曲线判断稳定性写主循环时先加一段能量累计而不是先做动画。PEC 盒无损耗总电磁场能量在时间上应当保持恒定数值上会有缓慢漂移。如果能量曲线是单调上升的指数形状一定是 CFL 超出或边界置零顺序写反energy np.sum(eps0 * (Ex**2 Ey**2 Ez**2) mu0 * (Hx**2 Hy**2 Hz**2))把这段代码放进每个时间步的末尾打印前 50 步的能量变化率。正常时变化率在 1e-6 量级如果几大步内能量翻倍立刻停下查 dt 和 courant 系数。5.3 对称性检查与最小冒烟测试偶极子放在盒子正中心时xy 切面场强应围绕中心呈旋转对称。每次改完边界条件或源位置先跑 20 步、保存三个切面的 png 核对对称性再跑完整动画。这一步能在生成 GIF 之前把边界贴错、数组转置这类问题定位到具体坐标轴xy 切面不对称查 z 方向边界yz 切面不对称查 x 方向边界。把能量检查和对称性检查写成一个 assert放在 FDTD 循环之后输出不合格时直接终止后续换网格尺寸、换激励频率时就不用重新人肉看动画判断对错了。本文还有配套的精品资源点击获取