简介这是一份面向电磁场数值模拟学习者与初研人员的MATLAB源码包围绕频域有限差分FDFD与时域有限差分FDTD方法聚焦波导与介质结构的色散特性计算。压缩包内仅含1个m脚本文件体量约1KB属于轻量级演示程序适合快速阅读代码逻辑并直接运行复现色散曲线。脚本覆盖网格定义、电磁场初始化、材料参数设置、激励源配置、迭代求解以及后处理绘图等完整环节运行后可得到折射率或介电常数随频率变化的色散曲线帮助理解不同波导模式的传播特性、截止频率与材料响应。已有543人学习下载对于光子晶体、光纤通信、微波器件设计等入门场景是一个便于拆解FDFD/FDTD算法实现细节、可作为教学与自学的实用参考样例。1. 一套包同时装下FDFD和FDTD为什么波导色散分析要双引擎并联做波导色散分析的人通常有一个共识光靠单一仿真工具很难在半天之内拿到一张可靠的色散曲线图。手里这份 FDTD.zip 把 FDFD 本征模式求解和 FDTD 时域宽带激发放到同一条流程里解决的正是从模式计算到色散曲线图绘制的完整链路。它直接对准硅光波导、光子晶体波导和 AR 光波导模组的初速设计——先用 FDFD 快速拿到模式场与有效折射率再用 FDTD 验证群速度色散和损耗趋势避免模式问题拖到后段仿真才暴露。适合需要自己跑通全套波导色散分析、又不想在商业软件里反复点界面的从业者新手可以照着示例脚本跑通第一个结构熟手则可以把 FDFD 当快速预筛、FDTD 当复核引擎边界条件、网格和 GPU 参数都在后面章节展开。2. FDFD模式求解器从亥姆霍兹方程压成矩阵特征值问题的落地路径FDFD 的思路说起来很直接假定波导沿 z 方向不变把电磁场写成横向分布乘以 e^{iβz} 之后三维麦克斯韦方程组就塌缩成一个只含横向坐标的二维算子方程。这个算子的特征值就是 β²特征向量就是模式场。难点不在概念而在怎么把连续方程变成一块能跑的稀疏矩阵。2.1 横向网格、差分格式与材料分布的三步设置第一步是确定横向计算窗口和网格密度。窗口太小会把模式场截断有效折射率偏大窗口太大则矩阵规模失控。常见做法是按工作波长取窗口边长在 6~10 μm 之间网格步长 20~40 nm。对于 λ 1.55 μm 的硅光波导这个组合能保证每个波长方向上有 50~80 个采样点色散曲线已经足够平滑。% 1_define_cross_section.m % 定义仿真区域与结构单位统一用 um Lx 6; Ly 4; nx 150; ny 100; % dx Lx/nx 0.04um, dy Ly/ny 0.04um x linspace(-Lx/2, Lx/2, nx1); y linspace(-Ly/2, Ly/2, ny1); [X, Y] meshgrid(x(1:end-1), y(1:end-1)); % 硅脊波导脊宽 0.5um厚度 0.22um下层是 SiO2上层是空气 index_map 1.444 * ones(ny, nx); % SiO2 包层 index_map(Y -0.11 Y 0.11 abs(X) 0.25) 3.47; % Si 脊 index_map(Y 0.11) 1.0; % 上包层空气 imagesc(x(1:end-1), y(1:end-1), index_map);这段代码里最容易被忽略的是单位全部用 μm 而不是 nm后面算 k0 时就不会出现 10 的负几次方级别的数值灾难。网格尺寸 0.04 μm 在大多数 1.55 μm 的硅光结构上都够用但如果算的是高折射率对比的光子晶体建议缩到 0.02 μm代价是矩阵规模变成四倍。材料分布用一个二维矩阵存后续不管做 FDFD 还是 FDTD 都共享这份折射率图改结构只需要改index_map这一段不用动主流程。2.2 组装算子矩阵与eigs特征值提取的参数细节确定完结构后第二步是组装离散化的亥姆霍兹算子。二维横向拉普拉斯算子可以用二阶中心差分构造再通过 Kronecker 积组合成一个大稀疏矩阵。折射率分布则放在对角线上与差分矩阵相加得到最终的算子 A。然后调用特征值求解器取能量最低的几个本征模。% 2_assemble_fdfd.m % 组装 FDFD 特征值问题A * E beta^2 * E lam0 1.55; % 工作波长 um k0 2*pi / lam0; dx Lx/nx; dy Ly/ny; Tx (1/dx^2) * spdiags([-1 2 -1], -1:1, nx, nx); Ty (1/dy^2) * spdiags([-1 2 -1], -1:1, ny, ny); A_op kron(speye(ny), Tx) kron(Ty, speye(nx)); A_op A_op k0^2 * spdiags(index_map(:).^2, 0, nx*ny, nx*ny); % 求解能量最低的 6 个模式 opts.tol 1e-10; opts.maxit 2000; [V, D] eigs(A_op, 6, sm, opts); beta sqrt(diag(D)); % 传播常数 neff beta / k0; % 有效折射率这里有几个参数要专门说明。eigs的sm表示寻找模最小的特征值也就是物理上最接近导波截止的模式这正好对应我们需要的最低阶模式如果用了lm你会拿回一堆高频伪模式。opts.tol 1e-10是收敛阈值调太松1e-6时特征矢量方向会有微小摇摆虽然 β 看起来没变但模式场的对称性会被破坏后面对接 FDTD 激励源时会引入偏振串扰。opts.maxit给到 2000 一般足够窗口再大一倍时建议提高到 5000。这个矩阵规模是 15000×15000在普通 PC 上 eigs 大约几十秒内能出结果。如果改用 0.02 μm 网格矩阵变成 60000×60000时间会涨到三五分钟这属于可以接受的范围。需要注意index_map(:)的展开顺序必须和 kron 算子的排列顺序一致否则折射率图会被转置——我自己踩过一次表现在算出来的 TE 模式场横竖颠倒色散曲线倒是没变化因为 β 对空间方向不敏感这种错位很难发现。2.3 边界条件如何左右色散曲线PML厚度与σ参数FDFD 的边界处理几乎是整个计算准不准的分水岭。默认包里的示例脚本常用零边界条件相当于计算窗口四周是理想导体壁模式场碰到壁会反射泄漏模会误判成导模。更常见的做法是在窗口外圈加卷积 PML完美匹配层用复坐标拉伸把外向行波吸掉。% 3_add_pml.m % 在六个网格外边带上施加 PML 复坐标拉伸 pml_thickness 20; % PML 占网格数 sigma_max 0.8; pml_scale sin(linspace(0, pi/2, pml_thickness)).^2; sigma_x [pml_scale*sigma_max, zeros(1, nx), fliplr(pml_scale*sigma_max)]; % 然后把 sigma_x 并进微分算子的分子项形成复坐标微分 % 核心做法对 Dx 的对角项做替换 Dx_pml Dx / (1 1i*sigma_x/omega);这段代码里的sigma_max是关键参数取大了低频反射变大取小了 PML 形同虚设。经验值是sigma_max ≈ 0.8 * omega窗口网格数 150 左右时 PML 厚度 20 个网格能稳定把反射压在 -60 dB 以下。PML 参数是玄学味道最重的地方不同网格密度下最优值浮动很大我一般用扫描的方式找固定结构把sigma_max从 0.4 到 1.2 按 0.1 步进比较基模 neff 随 PML 厚度的漂移量漂移小于 1e-5 就说明边界收敛。3. FDTD跑色散宽带脉冲激发、频谱相位提取与曲线拟合FDFD 给出的色散曲线是在频域单点求解的算多少个波长就要跑多少遍。FDTD 的思路完全不同往波导里打一个宽带脉冲时域推进后在两个不同 z 位置放探针记录场随时间演化再做 FFT 提取每个频率分量的相位差一次仿真拿到全波段 β(ω)。这正好和 FDFD 互补也是最容易出细节差异的环节。3.1 时域源的时宽、中心频率与监视器布点时域源通常用高斯调制脉冲数学形式是载波乘高斯包络。高斯包络的时宽 τ 直接决定频域覆盖宽度τ 越小频谱越宽能覆盖的波长范围越大但注入能量分散靠近波段边缘的信噪比会变差。% 4_fdtd_source_and_probes.m % 高斯调制脉冲中心波长 1.55um覆盖 CL 波段 lam0 1.55; f0 0.1938; % 中心频率单位为 THz 的简化表示 tau 0.6; % 脉冲时宽 fs 量级越小频带越宽 dt 0.02; % 时间步长由 CFL 条件决定 Tmax 4000*dt; t (0:dt:Tmax); src exp(-((t - 500*dt)/tau).^2) .* cos(2*pi*f0*(t - 500*dt)); % 监视器布置在 z 5um 和 z 5.8um z1 5.0; z2 5.8;监视器布置有个反直觉的点两点间距不能太大。许多人以为间距越大相位差越准实际上一旦间距超过几个波长导模的损耗和非理想色散会让两探针看到幅度严重不匹配的场分布相位提取误差反而变大。我常用的间距是 0.6~1.2 μm折中一下幅度衰减和相位分辨率。另外两个探针要避开光源平面和结构端面至少要隔开 2 μm否则入射脉冲和反射脉冲混在一起谱峰展开后相位会乱跳。3.2 两探针相位差反推传播常数β(ω)拿到两个探针的时域信号后处理过程分两步先对 e1(t) 和 e2(t) 分别做 FFT再在频域逐点计算E2(f)/E1(f)的相位角。这个相角除以探针间距就是该频率下的 β。% 5_extract_beta_from_phase.m E1 fft(probe1); E2 fft(probe2); f_axis (0:length(t)-1) / (length(t)*dt); % 频率轴 phase_diff unwrap(angle(E2 ./ E1)); beta_rad_per_um phase_diff / (z2 - z1); % 单位 rad/um lambda_axis 1 ./ f_axis(positive_freq_idx); neff_fdtd beta_rad_per_um(positive_freq_idx) ./ k0(positive_freq_idx);逻辑上要注意 E2/E1 的相位差除以间距得到的是 β而不是 neffneff 要用 β 除以真空波数 k0 换算。unwrap这一步最容易被漏掉相位角本身被限制在 [-π, π] 区间里只要两点间距超过半个导波波长真实相位差就会超过 π不 unwrap 的话整条 β 曲线呈锯齿状。代码里的positive_freq_idx对应去掉频谱负半边和直流分量奈奎斯特频率附近的数据点噪声很大一般只取频谱幅度高于峰值 1% 的频率区间。3.3 FDFD与FDTD结果互校误差通常在网格分辨率以内把 FDFD 的 β(λ) 和 FDTD 的 β(λ) 画在同一张图上后两条曲线通常在 1e-4 的相对误差量级内重合。如果偏差明显偏大优先怀疑 FDFD 的 PML 吸收没调好或者 FDTD 的探针间距超出了有效范围。实际操作中我会把两个结果相减得到一条残差曲线残差在高频端开始发散说明 FDTD 的脉冲频带边缘信噪比不足残差在全波段整体平移说明 FDFD 的窗口边界截断了模式场。这种互校要比单跑任何一种方法都可靠也就是这套包里两套引擎放在一起的真正价值。4. GPU加速与参数取舍FDTD开加速FDFD抠矩阵波导仿真越做越大之后GPU 加速是绕不开的话题。搜索热词里出现“fdtd 怎么开启 gpu”不是偶然——FDTD 的时域递推天然是数据并行的改造成 GPU 计算的收益非常直接但 FDFD 的处境不太一样不是所有环节开 GPU 都划算。4.1 从CPU到GPU时域算子的加速比与配置方式FDTD 主循环的核心操作是电场、磁场交替更新每个网格点的更新只依赖相邻网格的旧值这本质上是 stencil 运算正好是 GPU 的强项。常见做法是把场数组放到显存再用 CUDA 或 OpenCL 内核驱动 H 场到 E 场的更新。加速比通常在 20~50 倍具体取决于网格规模和 GPU 型号。% 6_gpu_fdtd.m % 将场数组移动到 GPU 显存核心循环不变 gpuDevice(1); Ex gpuArray(zeros(ny, nx, nt)); % 只保留当前时间层的场 Hy gpuArray(zeros(ny, nx, nt)); for n 1:nt-1 Ex(:,:,n1) Ex(:,:,n) c0*dt/eps_r .* ... (Hy(:,:,n1) - Hy(:,:,n)) ./ dy; % 每次时间步只更新一层避免整个时间历史占用显存 end注意这段代码里我把时间维也显式写出来了但实际应用时不要真的分配nt层的三维数组GPU 显存顶不住长时间运行。正确的做法是只保留当前层和下一层每步推进后把旧层覆盖掉。显存里真正吃空间的是近场监视区里需要整段保存的场分量那个才值得为它买 24GB 显存的卡。4.2 CFL条件与网格细化稳定性步长怎么定开 GPU 不改变一个铁律时间步长必须满足 CFL 条件否则数值不稳定跑几千步之后场值直接爆成 NaN。两维和三维的 CFL 形式略有区别# 三维 CFL 条件上限单位与网格步长保持一致 dt_max 1 / (c0 * sqrt(1/dx^2 1/dy^2 1/dz^2))网格每加密一倍dt 要缩小一倍而网格总数会变大八倍总计算量膨胀得更厉害。这就引出一个实用权衡不是网格越细越好而是要保证每个波长方向至少 20 个采样点。对 1.55 μm 的硅光波导dxdydz0.03 μm 时 dt 约 0.017 fs跑 4000 步只能覆盖约 70 fs 的时间窗口这个窗口对普通波导够用但对高 Q 微环来说远远不够微环的储能时间动辄几个皮秒必须把仿真步数提到上万步。4.3 稀疏特征值问题在GPU上的尴尬什么时候不值得开加速FDFD 的场景和 FDTD 完全不同。eigs求解稀疏特征值问题时主导运算不是稠密矩阵乘法而是稀疏矩阵向量乘这类操作极度受限于显存带宽GPU 上的稀疏格式比如 CSR访存效率反而不如 CPU 的缓存友好型实现。小规模矩阵几万阶在 GPU 上往往只快 1.5 倍甚至更慢还要花大把时间做矩阵格式转换。只有当窗口网格数超过 400×400 时 GPU 加速才开始体现出明显优势。我的做法是FDFD 保留 CPU 上的稀疏求解器只用 GPU 加速 FDTD。这套包里的示例脚本也基本是这个思路——FDFD 负责快速预筛和模式离散求解FDTD 负责宽带验证和损耗评估两者对硬件的需求正好错开。5. 避坑排查FDFD/FDTD色散计算里五个绕不开的翻车点这五个问题是我在不同项目里反复撞到过的每条按现象、原因、解决三个层次写你可以直接对号入座。5.1 色散曲线出现“平白无故”的平坦分支现象FDFD 算出的色散曲线在某段波长区间出现一条几乎水平的附件分支模式场图是密集的蜂窝状纹理。原因非交错网格下电场与磁场分量没有错开半个网格高频分量形成伪模式。解决换成 Yee 交错网格场分量分别定义在半网格位置上检查模式场时对场分布做一次快速傅里叶滤波剔除空间频率接近网格截止的分量。5.2 有效折射率低于包层折射率模式“漏”到计算窗口外现象基模的 neff 居然低于 1.444曲线整体下移。原因零边界条件把泄漏模也算进了特征解集或者 PML 吸收层太薄。解决把 PML 厚度加到 30 个网格用前面 2.3 节的扫描方法检查 neff 收敛性另外在筛选特征解时加一个条件——模式能量占比落在波导芯区 80% 以上才保留。5.3 FDTD 提取的 β 曲线呈周期锯齿现象β(ω) 曲线像锯齿一样周期跳动但 FDFD 的结果平顺。原因两个探针间距超过半个导波波长反三角运算的相位差被折叠回 [-π, π]没有正确unwrap。解决压缩探针间距到半个导波波长之内或者在提取相位后用相位连续判据做unwrap如果改用更长的探针间距必须同时在频域上做相位解缠。5.4 宽谱脉冲在波段边缘的信噪比崩掉现象FDTD 色散曲线在中心波长附近和 FDFD 很吻合但越往波段两端残差越大。原因高斯包络在频域是钟形分布边缘频率的激励能量过低。解决把脉冲时宽 τ 调小以扩展频带宽度或者直接扫两次仿真——第一次覆盖 C 波段第二次覆盖 L 波段然后拼起来。每个频段的边缘各砍掉 10% 的带宽剩余区域残差通常能压到 1e-5 以内。5.5 时域推进几百步后场值爆成 NaN现象FDTD 跑了不到一千步就开始出现 NaN越往后越密。原因时间步长越过 CFL 上限或者边界处材料参数出现不连续导致差分算子发散。解决先按 4.2 节公式计算 dt_max取 0.8 倍余量作为实际时间步长再检查结构边缘折射率是否跳变过陡必要时在边界加两层介质插值。6. 进阶脚本化addmode模式重注入与偏振转化效率自动扫描走到这一步FDFD 和 FDTD 的基本链路已经通了但真正进入器件设计阶段后还需要一个能力把前一步算出的特定模式作为激励源重新注入 FDTD并自动统计不同偏振分量的透射。这部分就要靠脚本语言来驱动了热词里常说的addmode就是干这个的。6.1 把FDFD模式场映射成FDTD激励源实际做偏振转化效率计算时以 TE 模式输入为例传统做法是直接塞一个平面波进波导但平面波的场分布和波导本征模式重叠率低大部分能量会散射进包层仿真出来的转化效率系统性偏低。正确的做法是把 FDFD 算出的模式场剖面提取出来映射成 FDTD 的模式源。% 7_addmode_and_inject.m % 由 FDFD 求解特定波长 TE0 模式并注入 FDTD mode_field fdfd_solve(struct_geo, lam_target); fsp mode_field.Ex 1i * mode_field.Ey; % 模式横截面 source_plane inject_pattern(fsp, source_dimension, 1); % 在脚本接口里这就是 addmode 背后做的映射加一个指定模式分量的激励如果是在商用 FDTD 的脚本语言里addmode一般接受两个参数模式编号和位置。默认情况下不调用它时模式源会使用基模分量导致你算了半天 TE 到 TM 转化实际注入的场里混着一半高阶模式。建议每次新建仿真文件时把端口模式列表打印出来确认注入的模式编号和 FDFD 里解出来的 mode index 一致。6.2 偏振转化效率的谱扫描与参数化偏振转化效率 PCE 定义为中心波长下输出端口目标偏振光功率占总透射光功率的比例。要得到整个波段的 PCE 谱把前面 FDTD 的宽谱仿真结果做重叠积分就行。% 8_pce_sweep.m % 自动扫描波长调用 FDFD 得模式FDTD 得透射场再算 PCE lambda_list 1.40:0.02:1.60; for idx 1:length(lambda_list) lam lambda_list(idx); mode fdfd_solve(struct_geo, lam); [T_main, T_cross] fdtd_mode_propagate(lam, mode); pce(idx) abs(T_cross)^2 / (abs(T_main)^2 abs(T_cross)^2); end plot(lambda_list, pce, -o); hold on; plot(lambda_list, 1 - pce, --);逻辑上需要注意T_main和T_cross必须是同一模式下两个正交偏振分量的透射系数不要用两个不同模式编号的结果硬凑。我在实际项目里见到的最大误区是直接拿总场做积分代替模式交叠积分这会把包层辐射和金属吸收全算进“转化效率”数值大得离谱却没有任何物理意义。正确做法是用模式监视器在输出端把场分解到 TE 和 TM 本征模式上再分别积分。6.3 自动优化循环里的一个习惯这套流程的最终形态是参数化扫描外层循环改脊宽或刻蚀深度内层调用 FDFD 快速求模式再映射到 FDTD 验证 PCE。把 6.1 和 6.2 的脚本封装成函数之后一个二维扫描大约需要跑几百次 FDFD 加几十次 FDTD在 GPU 加持下可以在一两天内完成。从那以后我每次接波导项目都强制先跑一遍“FDFD 快扫 FDTD 宽谱复核 模式注入算偏振”的完整链路色散数据和偏振数据一起要不再让模式问题拖到后段仿真才暴露希望帮到你。本文还有配套的精品资源点击获取