最近在整理光学仿真案例时翻到一套很有意思的Matlab源码用FDTD方法模拟双缝干涉并且把PML和PMC两种边界条件放进了同一个项目里。这套源码还附带一份光学综合实验报告拿到手就能直接跑适合拿来练手或做课程设计。双缝干涉是波动光学里最经典的现象之一教材上讲起来很简单但真要用数值仿真复现出来涉及的细节其实不少网格怎么剖分、边界怎么吸收、光源怎么注入、条纹怎么提取每一步都能踩坑。PML负责吸收外边界反射PMC用来利用对称性压缩计算域两者配合可以把一个看似“高大上”的时域仿真变得非常可控。这篇文章就把这套项目从物理原理、参数设计、Matlab实现到报告写作完整走一遍希望对正在做光学综合实验、FDTD入门或者双缝干涉仿真的同学有帮助。1. 双缝干涉与FDTD为什么要把两个概念放在一起1.1 公式能算为什么还要做数值仿真双缝干涉的教材结论非常干净两束相干光叠加强度分布是一个cos²项乘上单缝衍射的sinc²包络主极大位置由缝间距和波长决定。可一旦落到实物上问题就变了。缝不是无限薄金属屏也不是理想平面观察屏可能离缝只有几个波长边缘绕射、倏逝波、界面反射全都混在一起。解析公式只能帮你估算亮纹角度真要看电磁场怎么在狭缝附近传播、怎么在后方重新叠加还是得求解麦克斯韦方程组。FDTD时域有限差分就是干这件事的把空间切成很多小网格时间也切成很多小步反复迭代更新电场和磁场直到电磁波在计算域里稳定下来。它的优势是直观可以随时把某一时刻的场图导出来看波前、衍射、干涉过程一清二楚。对于双缝干涉这种既包含传播又包含边界效应的场景FDTD比纯解析公式更接近真实问题。这套源码的定位就是“用数值方法复现经典光学现象”而不是简单画一条理论曲线。1.2 PML和PMC在项目里的实际分工FDTD只能在有限区域里计算直接截断的话电磁波打到边界会产生虚假反射结果会完全变样。PMLPerfectly Matched Layer就是为了解决这个问题在计算区域外面加一层吸收材料让进入PML的波被逐渐衰减几乎不反射回计算域等效于把边界推到无穷远。实际使用中PML层厚度、电导率渐变曲线都需要仔细调不然反射照样存在。PMCPerfect Magnetic Conductor则是另一种思路它不吸收波而是模拟理想磁导体边界。当仿真模型存在对称性时可以在对称面上设置PMC让计算域只保留一半甚至四分之一从而节省内存和计算时间。这个项目里很有意思的一点就是把外边界PML和内部对称面的PMC组合使用PML负责吸收开放边界PMC负责压缩模型规模两者各管一件事。如果全模型四周都用PML也能跑通但加上PMC之后网格数量直接下降跑起来明显快很多。1.3 双缝“干扰”还是“干涉”先纠正一个用词。标题里写的“双缝干扰”我理解就是双缝干涉英文interference本身兼有“干涉、干扰”两层含义有些软件界面或早期译法会把interference pattern写成“干扰图样”。国内物理教材一般统一用“干涉”所以下文都用“双缝干涉”。做仿真时候没必要纠结这个词重点在物理模型和边界条件。2. 仿真模型搭建参数选型与边界配置2.1 几何模型与光源设计这套源码搭建的是二维FDTD模型计算平面是x-z平面狭缝沿y方向无限延伸。双缝屏放在z0附近光源从z负半轴一侧入射透射后在z正半轴形成干涉条纹。几何参数不是随便拍的要保证条纹清晰、计算量可控。我复盘时用的典型参数如下参数符号数值真空波长λ633 nm缝宽a800 nm缝间距d4 μm金属屏厚度t100 nm网格步长dx dz25 nm计算域Lx × Lz16 μm × 32 μmPML厚度Npml10 层光源类型—连续正弦平面波缝宽取800 nm略大于波长单缝衍射包络比较宽不会把双缝干涉条纹“压死”。缝间距取4 μm这样在观察屏上能分出至少五六条亮纹。网格步长取25 nm大约是λ/25既满足精度又不会让Matlab数组太大。如果你用的是16GB内存的电脑这套参数跑起来非常流畅。光源我建议用连续正弦波而不是高斯脉冲。脉冲源适合宽频分析但要看单一波长下的稳态干涉条纹连续正弦波更直接。入射电场设置为沿y方向偏振也就是电场方向与狭缝长度方向一致。这样在二维FDTD里只需要更新E_y、H_x、H_z三个分量计算量最小。2.2 网格尺寸、时间步长怎么定FDTD的网格尺寸直接影响数值色散。一般经验是网格步长至少取λ/10想算得准一点就取λ/20甚至λ/30。这里取λ/25已经是比较稳妥的选择。网格太粗的话狭缝边缘会被“台阶化”缝的有效宽度都不准干涉条纹的位置自然对不上。网格太细则计算时间成倍增加需要根据电脑性能权衡。时间步长要满足CFL稳定性条件。二维情况下有[ \Delta t \le \frac{1}{c \sqrt{\frac{1}{\Delta x^2} \frac{1}{\Delta z^2}}} ]我习惯取这个上限的0.99倍留一点余量避免数值稳定性边缘出现问题。用上面参数算下来一个时间步大约是0.06 fs量级跑3000步也就180 fs左右。但要注意电磁波要从光源位置传播到缝屏再穿过整个观察区域必须等到全计算域都建立稳定场之后才能采集数据。我实际跑下来发现至少需要跑4000到5000步条纹才会干净。如果只跑一千步结果里会带着明显的“建立过程”痕迹看起来就像条纹在抖动。2.3 PML与PMC的边界组合方式这套源码里计算区域四周默认都设置PML。PML的厚度我建议从8层起步10到16层更稳。PML内部电导率按位置渐变增加阶数通常取3到4这样从计算域到PML内部阻抗过渡更平滑反射率可以压到-60 dB以下。太薄的PML虽然省内存但会引入可见反射干涉条纹上会叠一层周期性的“毛刺”。PMC边界的使用要结合模型对称性。双缝屏和入射波关于中心平面x0是对称的TM偏振正入射时对称面上的切向磁场为零正好满足PMC边界条件因此可以只模拟x≥0的一半区域。这样原来两条缝的模型就变成一条缝加一个PMC对称面计算域直接砍掉一半。算完之后再把结果镜像翻转就能得到完整干涉图样。源码里同时保留了全模型PML和半模型PMLPMC两套设置方便对比验证。3. Matlab源码实现从Yee网格到条纹输出3.1 主循环结构Matlab实现FDTD的核心是Yee网格迭代。代码主体通常由三部分构成初始化场数组和材料参数、进入时间步循环、循环结束后处理数据。这套源码的主循环结构大致如下% 初始化 Ey zeros(Nx, Nz); Hx zeros(Nx, Nz); Hz zeros(Nx, Nz); % 时间步进 for n 1:Nt % 更新 Hx 和 Hz Hx(2:end-1, 2:end-1) Hx(2:end-1, 2:end-1) ... - dt / mu0 / dz * (Ey(2:end-1, 3:end) - Ey(2:end-1, 1:end-2)); Hz(2:end-1, 2:end-1) Hz(2:end-1, 2:end-1) ... dt / mu0 / dx * (Ey(3:end, 2:end-1) - Ey(1:end-2, 2:end-1)); % 更新 Ey Ey(2:end-1, 2:end-1) Ey(2:end-1, 2:end-1) ... dt / eps0 * ( ... (Hz(2:end-1, 2:end-1) - Hz(1:end-2, 2:end-1)) / dx ... - (Hx(2:end-1, 2:end-1) - Hx(2:end-1, 1:end-2)) / dz); % 源注入与边界处理 % ... end上面只是示意实际数组维度要对齐Yee网格但核心逻辑就是这样先更新磁场再更新电场。材料参数中金属屏区域可以用大的电导率近似PEC也可以直接在电场更新时把对应位置的场强强制置零。两种写法都能出结果但后者更简单跑起来也更快。3.2 PML和PMC的代码落地PML的代码不是简单地“给边界乘一个小数”那样会造成阻抗不匹配。工程上常用分裂场PML或CPML。这套源码采用的是经典的PML实现方案在PML区域内把场分量拆分并引入辅助变量记录吸收项。电导率沿PML深度方向呈多项式渐变比如sigma_x sigma_max * (x_in_pml / thickness_pml)^order;sigma_max和阶数需要调。我的经验是阶数取3sigma_max取几个S/m量级反射率就能压得很低。如果你初次写PML可以直接套用已有模板重点看两个地方一是PML区域内的场更新公式是否加了辅助变量二是PML与内部区域交界处有没有“硬接线”式的突变。PMC的代码实现比PML简单很多。在对称边界处只需要把边界上的切向磁场分量强制置零或者把磁场更新时跨过边界的项做镜像处理。举个例子如果PMC边界在z1这一层那么需要保证这一层边界上的切向H为0更新代码里区分“内部节点”和“边界节点”即可。这种处理方式比PML直观多了但前提是物理上确实满足对称条件。3.3 从场分布到干涉条纹的提取跑完时间循环后不能直接看某一时刻的瞬态场因为瞬态场里既有入射波也有反射波条纹并不干净。需要等稳态建立后在一个完整周期内记录电场幅度值。常见做法是记录场强的峰值包络或者对时间序列取幅度field_amplitude sqrt(mean(Ey_record.^2, 3));在观察屏位置也就是距离缝屏几个波长的地方沿x方向取这一行场强就得到干涉条纹。如果你只想看相对强度分布直接用场强平方就行。这里有个小细节观察屏不要紧贴缝屏否则近场倏逝波成分会混入结果也不要太靠计算域边缘否则PML的微小反射会污染条纹。一般来说离缝屏5到10个波长最合适。如果想把近场结果换算成远场角分布可以对近场记录做空间傅里叶变换。空间频率k_x和角度θ之间满足k_x k0 sinθ。Matlab里用fft函数就能实现但要注意先对观察屏数据进行适当的窗函数处理避免截断效应在角度谱上产生旁瓣。3.4 报告撰写思路这套源码附带的光学综合实验报告写作框架其实很固定。我建议按下面这个顺序组织报告章节写什么注意事项实验目的复现双缝干涉学习PML与PMC边界条件不要只写“仿真双缝干涉”要突出FDTD和边界条件原理部分双缝干涉公式、FDTD迭代格式、PML/PMC作用公式要写清楚但不需要推导太深仿真模型参数表、几何示意图、边界设置参数表必须和源码一致方便复现结果与分析场图、条纹强度曲线、PML/PMC对比图和曲线要标坐标轴分析要对应物理现象结论总结仿真结果和边界条件作用避免泛泛而谈要写“验证了”“得到了”报告里最加分的是对比图比如PML厚度8层和16层的结果对比或者全模型和半模型PMLPMC的结果对比。把这两张图放进去老师一眼就能看出你真的理解了边界条件的作用而不是只会跑代码。4. 实操中踩过的坑与排查思路4.1 PML厚度不足导致的“伪条纹”我第一次跑这套源码时为了省内存把PML层数从10层降到了6层。结果出来的干涉条纹整体上是对的但亮纹之间多了一层细密的周期波纹像是“条纹上面叠了另一组条纹”。后来检查发现就是PML太薄波打到边界后反射回计算域和透射波再次干涉形成了驻波。把PML加到12层同时把电导率渐变阶数从2改成3那层伪条纹立刻消失。所以如果你发现亮度分布上出现均匀的高频波纹第一个怀疑对象就是PML反射。判断方法很简单把PML厚度增加到16层或者把计算域扩大一倍再看波纹是否变化。如果变了基本可以确定是边界反射问题。4.2 稳态建立时间不够FDTD是迭代求解源开始注入后波前需要一段时间才能到达观察区域。如果时间步数太少你看到的条纹其实只是波前刚经过观察屏时的“瞬态快照”条纹间距可能对但边缘会有明显的拖尾或不对称。我常用的判断标准是在观察屏附近选一个固定点记录该点场强随时间的变化当曲线变成稳定的正弦包络时才代表稳态建立完成。为了保险起见我通常把总时间步数设置成“光穿过整个计算域所需时间”的三到四倍。用前面的参数算光从左侧边界到右侧PML大约要走48μm需要160 fs左右跑5000步、大约300 fs已经足够稳。4.3 对称边界选错导致结果整体偏移PMLPMC的省内存方案很诱人但对称边界不是随便设的。如果入射偏振或者结构对称模式不对PMC会算出一个对称性错误的解。比如本应该用PEC的地方用了PMC干涉条纹的强度和位置都会发生变化甚至主极大变成次极大。排查办法很直接先用全模型PML跑一遍再把PMC半模型结果镜像对齐和全模型逐点比较。如果两幅图重合边界条件就是对的如果不重合就要检查源和结构的对称性。这里多说一句源码里针对TM偏振选用PMC是因为该偏振下对称面满足切向磁场为零。如果你换到TE偏振情况就会反过来可能需要用PEC。不要背着皮尺照抄一定要理解边界条件的物理含义。4.4 数值色散与分辨率不足网格步长太粗会让电磁波的速度产生误差表现为条纹位置偏移、窄缝传播角度不对。我试过用λ/10的网格跑同样模型结果次级明纹位置差了接近5%肉眼可见。把网格加密到λ/25后误差就控制到可以忽略的程度。当然加密网格会带来内存和时间成本。如果电脑跑不动可以先粗网格调通流程再加密网格做最终计算这是比较省事的路径。还有一个容易被忽略的点金属屏厚度在二维网格里至少要有两三个网格来表示。如果屏只有一个网格厚狭缝边缘会出现非常强的场增强导致结果对参数特别敏感。把屏厚度设为100 nm网格25 nm正好四个网格边缘效应相对稳定。5. 这次项目调试的小体会5.1 调试顺序比参数本身更重要拿到源码后不要急着改缝宽、改波长我建议先按默认参数跑一遍确认得到正确的双缝干涉条纹。然后在保持光源和结构不变的情况下分别做三组对照实验PML厚度翻倍、去掉PMC改用全模型PML、把PMC替换成PEC。这三组结果能帮你在最短时间内建立起对边界条件的直观感受。之后再回头调整几何参数你就会很清楚哪些现象是物理本身造成的哪些是数值设置造成的。很多时候参数调不出来不是数值问题而是模型和源没对。比如入射光不平整、金属屏建模有缝隙、监视器位置太近这些都会让结果看起来“不对”。所以调试顺序应该是先简单、后复杂先验证、再优化。5.2 这套源码还能怎么扩展这套项目本身是二维单频仿真但它的框架可以很容易往多个方向扩展。比如把连续正弦波改成高斯光束可以研究聚焦后的干涉图样把双缝屏改成光栅能直接看各级衍射效率把PML边界保留、对称边界改成周期性边界就能往超表面方向走。最近很多人在讨论FDTD偏振转化效率的仿真其实就是在这个框架里增加不同的偏振光源和各向异性材料模型。如果你有时间我建议你再做一个小实验保持缝宽和屏厚度不变只改变缝间距观察条纹间距怎么变化。这个实验很简单但做完之后你可以更深刻地理解双缝干涉公式里的d和条纹间隔的关系。源码里所有参数都在头部集中定义改一个变量重新跑一遍就行非常适合做参数扫描。我就拿这套源码做过一次缝间距扫描从2 μm扫到6 μm得到的主极大角位移和理论值几乎重合。那一刻才真正觉得FDTD不是黑箱只要你把参数吃透它就是一个可以反复实验的虚拟光学平台。