简介Meep 是一款基于有限差分时域方法的免费开源电磁仿真软件面向科研人员、工程师及电磁场专业学生可用于光子晶体、波导、谐振腔等微纳光学器件的建模与模拟支持一维、二维、三维及柱坐标仿真。资源包内共含六百七十一个文件其中以 Python 脚本、C 源文件、Scheme 控制脚本和 Markdown 文档为主另有大量结果图像与演示动画压缩包整体仅三十三点六一兆字节。已有七百六十二人学习使用。解压后可通过说明文档快速掌握编译配置与接口用法借助其中的交互式笔记本学习从基础到进阶的仿真案例。配套材料覆盖各向异性材料、吸收边界、亚像素平滑等核心特性有助于掌握仿真流程、二次开发技巧及性能优化思路适合希望深入电磁仿真与开源软件定制的开发者。1. 免费的 C FDTD 电磁仿真为什么值得自己动手FDTD时域有限差分是电磁仿真里最直白的一种方法把麦克斯韦旋度方程在网格上做时空离散按时间步一步步往前推不需要组装矩阵也没有迭代求解器。商用工具如 Lumerical FDTD 的授权费用不低安装流程也常劝退个人开发者而免费、用 C 写的 FDTD 项目其实不少openEMS、Meep、gprMax 都能从 GitHub 直接下载源码覆盖天线、光子学、探地雷达等场景。C 的价值在性能和可控性同样的网格规模下C 内核比 Python 纯循环快两个数量级而且你能看清每一个更新公式、每一层吸收边界的真实实现。下文按“原理落地、选型编译、仿真实战、验证调优”的顺序把一套可复现的免费 C FDTD 方案讲透。2. FDTD 原理与 C 最小实现从 Yee 网格到可运行代码2.1 Yee 网格与蛙跳更新为什么 FDTD 天然适合 CFDTD 的离散化思路是 1966 年 Kane Yee 提出的电场和磁场分量在空间上交错半个网格间距时间上也交替半个时间步推进形成“蛙跳”结构。网格上某个 Ez 分量只依赖周围四个 H 分量的差分某个 Hy 又只依赖两个 E 分量整个更新过程没有矩阵求逆、没有隐式方程组每个时间步就是若干层规则的内存访问和显式差分运算。这种形态恰好是 C 编译器最擅长优化的对象。内层循环可以自动向量化交错网格带来的访存模式也相当规整缓存命中率远好于稀疏矩阵类求解器。Python 和 MATLAB 更适合做几何建模、网格生成和结果后处理真正的时间步进内核几乎全部用 C/C 实现这也是所有主流开源 FDTD 项目的共同选择。理解这一点后再看任何一款免费 FDTD 软件的代码主线都非常清晰初始化场数组、循环更新 H、循环更新 E、注入源、处理边界。2.2 一维 FDTD 的最小 C 代码与参数拆解2.2.1 能直接编译运行的代码先写一个零依赖、单文件的一维 FDTD代表真空中的 TEM 波传播。把代码存成fdtd1d.cpp命令行执行g -O2 fdtd1d.cpp -o fdtd1d ./fdtd1d就能跑#include cstdio #include cmath #include vector int main() { const int N 1000; // 网格总数 const double dx 0.01; // 空间步长 1 cm const double c0 3e8; // 真空光速 const double dt dx / (2.0 * c0); // 时间步取 CFL 上限的一半 const int steps 1000; // 时间步数 std::vectordouble ez(N, 0.0), hy(N, 0.0); for (int n 0; n steps; n) { // 更新 H旋度项用相邻 Ez 的差分离散 for (int i 0; i N - 1; i) hy[i] (dt / (c0 * dx)) * (ez[i 1] - ez[i]); // 更新 E用相邻 Hy 的差分离散 for (int i 1; i N; i) ez[i] (dt / (c0 * dx)) * (hy[i] - hy[i - 1]); // 软源在域中心叠加高斯脉冲而不是直接赋值 ez[N / 2] std::exp(-std::pow((n - 200) / 50.0, 2)); if (n % 200 0) std::printf(step %4d ez(center) %g\n, n, ez[N / 2]); } return 0; }逻辑说明H 更新循环里ez[i1] - ez[i]是对 x 方向的差分离散乘上系数dt/(c0*dx)后等价于用时间差分近似时间导数E 更新循环对称地使用hy[i] - hy[i-1]。交错半个网格体现在索引偏移上Hy 在 Ez 的相邻中点取值这正是 Yee 网格在 1D 下的形态。源注入放在 E 更新之后用叠加而不是用赋值后面会解释原因。2.2.2 参数表dx、dt、steps 怎么定参数取值含义与依据N1000网格总数计算域长度为N * dx 10 mdx0.01 m空间步长按最高频率对应的波长取 1/20 左右dtdx / (2.0 * c0)时间步取 CFL 稳定性上限dx / c0的一半steps1000推进步数1000 步后波前恰好传播约 5 mdx 不是随便定的它必须远小于仿真频段内的最小波长否则空间离散会产生严重色散误差脉冲会“拉宽”。先定 dx 再定 dt顺序不能反。这里把 dt 取成上限的一半是留足安全系数三维情况还要考虑 y、z 方向的贡献。2.3 CFL 稳定性条件与源项注入方式三维均匀网格的 CFL 条件表达式是dt 1 / (c * sqrt(1/dx^2 1/dy^2 1/dz^2))一维情形退化为dt dx / c。工程上一般取上限的 0.90.99 倍宁可小不要大超过上限误差会随步数指数增长很快出现 nan 或溢出。判断是否失稳有个直观特征网格里某个点的场值在源已经停止注入后仍然持续快速增大那就是 CFL 超限了。源项注入分两种。硬源直接赋值ez[src] value相当于在这个网格点上放置了一个理想导体会产生不必要的二次散射一般只用于简单的定性演示。软源用ez[src] value叠加网格电导率不变适合需要干净入射场的场景。本文代码用的是软源。如果做散射体分析更规范的做法是总场散射场TF/SF边界把入射场限制在总场区散射场区只接收目标散射信号。2.4 边界处理从 MUR 到 CPML 的取舍一维代码运行到一定步数后波会撞上网格末端。默认情况下网格边界相当于理想导体波会被反射回来污染观测点数据所以真实仿真必须加吸收边界。最简单的是 Mur 一阶吸收边界用外向行波条件近似对垂直入射的波吸收效果尚可但斜入射时残差明显适合快速验证。目前开源软件的事实标准是 CPML卷积完全匹配层在计算域外圈加 812 层介质层电导率按多项式渐变通过记忆变量做卷积更新宽频带、大角度入射都有较好的吸收效果。选型时主要关注两个参数吸收层层数和电导率渐变的阶数。层数太少了低频反射大阶数匹配不对会在层间产生数值反射。后面实战部分用到的 openEMS其官方示例默认用 MUR追求低反射时再按版本 API 切换 PML 类边界。3. 免费 FDTD 软件怎么选、怎么下载编译openEMS、Meep 与 gprMax3.1 openEMS、Meep、gprMax 三个 C/C 项目怎么选免费的 C FDTD 项目里真正活跃且资料成体系的三个是 openEMS、Meep 和 gprMax它们的定位差异很大选错方向会浪费不少时间。项目内核与接口典型场景下载方式openEMSC 内核Octave/MATLAB 接口社区有 Python 封装天线、微波电路、电磁兼容GitHub 源码或 release 归档MeepC 内核Python 3 与 Scheme 接口光子晶体、超材料、集成光学源码编译或包管理器安装gprMaxC 核心加 Python 封装探地雷达、埋地目标回波PyPI 包或源码安装openEMS 的强项是微波和天线工程几何建模基于 CSXCAD端口、S 参数、近场远场变换都有现成封装对做射频的工程师最友好。Meep 偏向纳米光学和周期性结构非常擅长色散介质和任意阶非线性Python 接口写得干净适合做物理研究。gprMax 面向探地雷达内置了土壤色散模型和粗糙地面生成器是行业专用工具。许可证方面三者都是开源但协议条款不同商用前务必读一遍各自的 LICENSE 文件。3.2 用 CMake 从源码编译 openEMS 的最小步骤openEMS 是标准的 CMake 工程最新版本依赖 Boost、HDF5、CGAL、TinyXML 等库。在 Ubuntu 22.04 或更新发行版上的编译步骤sudo apt update sudo apt install -y build-essential cmake \ libboost-all-dev libhdf5-dev libcgal-dev libtinyxml-dev git clone --depth 1 https://github.com/thliebig/openEMS.git cd openEMS mkdir -p build cd build cmake .. -DCMAKE_BUILD_TYPERelease make -j4参数说明--depth 1只拉最新提交不带历史能明显缩短克隆时间-DCMAKE_BUILD_TYPERelease必须加FDTD 内核在 Debug 模式下慢一个数量级-j4是并行编译的核数按机器实际 CPU 调整。编译完成后openEMS 的 Octave 接口还需要把代码仓库里的 matlab 目录加入 Octave 搜索路径常见做法是在~/.octaverc里写一行addpath(/path/to/openEMS/matlab);。如果你的电脑是 Windows先别急着在原生环境折腾 CMake。常见做法是用 WSL 装 Ubuntu 再按上述步骤编译然后用 VSCode 配置好 C/C 环境去阅读和调试 openEMS 源码远程挂到 WSL 里看代码、改代码都比较顺。懒得编译的话去官方 release 页面看看有没有现成二进制包比自己从零拖依赖省事得多。3.3 下载源码时的三个坑第一个坑是直接下载 master 分支的 zip。master 通常处于开发态接口可能和文档不一致官方示例脚本跑不通。正确做法是认准 release tagexamples 目录下的脚本和 tag 是配套的。用git clone --depth 1 --branch tag名拉指定发布版稳妥。第二个坑是子模块。不少 C 科学计算项目用 git submodule 管理第三方依赖直接下载 zip 会漏掉子模块编译时缺头文件。克隆命令记得加--recurse-submodules参数。第三个坑是版本与教程错配。网上很多教程截图来自旧版 openEMSAPI 已经改名。下载后先看仓库里的 examples 目录以当前 tag 附带的脚本为准而不是以搜索引擎里的旧博客为准。4. 实战用 openEMS 跑通偶极子天线的 S 参数仿真4.1 可直接运行的 Octave 脚本含端口与网格openEMS 官方生态以 Octave 为主下面是一份精简过的半波偶极子脚本覆盖建模、网格、端口和求解全流程% dipole_openEMS.m physical_constants; % 脚本内定义 c0、eps0 等常量 unit 1e-3; % 几何单位mm f0 2.4e9; % 设计频率 2.4 GHz lambda c0 / f0 / unit; % 波长单位 mm fdtd InitFDTD(endCriteria, 1e-4); fdtd SetGaussExcite(fdtd, 0.5/f0, 0.5/f0); fdtd SetBoundaryCond(fdtd, {MUR,MUR,MUR,MUR,MUR,MUR}); csx InitCSX(); csx AddMetal(csx, dipole); % 定义 PEC 材料 gap 3; % 端口间隙 3 mm csx AddBox(csx, dipole, 10, [-lambda/4, -1, -1], [-gap/2, 1, 1]); csx AddBox(csx, dipole, 10, [ gap/2, -1, -1], [ lambda/4, 1, 1]); [csx, port] AddLumpedPort(csx, 10, 1, true, [-gap/2, 0, 0], [gap/2, 0, 0], 50); mesh.x sort(unique([-lambda/2 : lambda/20 : lambda/2, -gap : 0.5 : gap])); mesh.y -lambda/2 : lambda/20 : lambda/2; mesh.z -lambda/2 : lambda/20 : lambda/2; csx DefineRectGrid(csx, unit, mesh); Sim_Path sim; Sim_CSX dipole.xml; WriteOpenEMS(Sim_Path, Sim_CSX, csx); RunOpenEMS(Sim_Path, Sim_CSX); f linspace(2e9, 3e9, 201); [Zin, S11] calcPort(port, Sim_Path, f0, f); plot(f/1e9, 20*log10(abs(S11))); grid on;脚本关键点InitFDTD(endCriteria, 1e-4)表示当场内能量衰减到初始值的万分之一时自动停止这个值直接控制仿真时长和精度SetGaussExcite的第二个参数是脉冲宽度取 0.5/f0 能在 2.4 GHz 附近有足够带宽两段偶极子臂之间留了 3 mm 间隙AddLumpedPort在这个间隙里放置 50 欧姆集总端口。网格部分特意在-gap : 0.5 : gap区间加密保证端口处有网格线经过这是很多新用户跑不出 S 参数的直接原因。RunOpenEMS启动的是 C 求解内核仿真过程中的时间步进、CPML 更新全部在编译后的本地代码里执行Octave 只负责等待结果。跑完后calcPort读取端口处存储的时域电压电流做傅里叶变换后得到输入阻抗和 S11参考阻抗就是 AddLumpedPort 里写的 50 欧姆。用这份脚本仿真2.4 GHz 附近能看到明显的谐振凹陷。4.2 网格密度、时间步与内存估算网格密度决定精度也决定内存。起点是每波长 20 个网格也就是lambda/20对应多数天线问题的工程精度如果关心场强细节或弯曲结构几何区域加密到lambda/50并不夸张。openEMS 会自动根据网格推断时间步并把 CFL 条件控制在稳定范围内不需要手算 dt但你要有能力手工复核先看mesh三个方向的最小网格再用dt 1 / (c * sqrt(sum(1/min_dx^2)))粗算一遍确认库生成的时间步没有越过上限。内存占用可以用一个简单式子估算场数组需要的字节数约等于Nx * Ny * Nz * 6 * 8系数 6 是三维的六个场分量8 是 double 的字节数python3 - EOF nx, ny, nz 250, 250, 125 # 三个方向的网格数 cells nx * ny * nz mem cells * 6 * 8 print(f场数组约 {mem/1e6:.1f} MB) print(f含 8 层 PML 与网格结构约 {mem*1.3/1e6:.1f} MB) EOF这个估算不含 PML 层的额外记忆变量和几何网格结构本身后两者通常再带来 20%40% 的开销。如果估算结果接近物理内存上限优先降 nz 方向分辨率因为 z 方向对贴片天线这类结构影响相对小。4.3 常见报错与参数修正对照表现象原因处理方式启动即报 mesh 相关错误几何体与网格没有对齐端口处无网格线在端口坐标处强制插入网格线参考 4.1 的 mesh.x 处理S11 曲线噪声大、不收敛endCriteria 过大时域截断过早把1e-4降到1e-5观察结果变化内存不足网格总数超出物理内存先按 4.2 公式估算粗网格跑通后再加密S11 出现周期性波纹吸收边界反射确认边界设置必要时换 PML 类边界并加厚吸收层脚本接口对不上版本太新或太旧API 改名以仓库 examples 目录里同版本脚本为准遇到第 3 类问题时不要一次加密到位按 20% 幅度递增网格并观察 S11 变化比直接拉满网格效率高得多。5. 结果验证与 C FDTD 内核的性能调优技巧5.1 用解析解和收敛性检查给结果把关免费软件跑出的 S 参数不能直接采信第一步是和解析解对拍。半波偶极子的谐振频率理论值约等于c / (2 * L)其中 L 是偶极子总长。仿真结果里 S11 的谐振点和理论值的偏差应在百分之几以内偏差过大先查网格密度而不是查端口设置。第二步做收敛性检查把三个方向的网格同时减半再跑一次对比两次 S11 曲线。如果谐振频率偏移超过 1%说明当前网格仍然太粗加密后再重复一次直到两次结果差异进入可接受范围。这个方法对自写的一维代码同样适用观察两个观测点的波形时延应该严格等于距离除以光速。5.2 场更新循环里的内存布局与 OpenMP 技巧如果下载源码后想改自己的内核优先优化两个地方。第一是内存布局三维场数组建议按分量分开存储例如std::vectordouble ex(N), ey(N), ez(N), hx(N), hy(N), hz(N)而不是把六个分量打包成结构体数组后者在更新循环里会引入缓存行浪费。第二是循环顺序C 默认行优先最内层循环应沿内存连续的方向三维场景下通常把网格的 x 方向作为最内层。多核并行从 OpenMP 开始最划算。在 H 更新和 E 更新的外层循环前加#pragma omp parallel for注意两个更新循环之间保持隐式同步不要在同一个循环里同时读写同一个场分量。一个更具体的技巧是用-O3 -marchnative重编 openEMS 内核编译器会针对本机 CPU 启用 AVX2 或 AVX-512 向量化这一步对时间步进循环的加速非常明显且不需要改任何业务代码代价是换机器后需要重新编译。验证边界效果时把吸收边界层数从 8 加到 12比较同一条 S11 曲线的偏移——网格加密一倍S11 曲线偏移小于 0.1 dB 之前先别急着换更复杂的边界条件把这一步做完免费 FDTD 的仿真结果才敢写进报告。本文还有配套的精品资源点击获取
