1. 一维光子晶体到底在算什么收费站模型和它的物理边界第一次听到“一维光子晶体像光子的高速公路收费站”这个说法是在一次组会上。我当时的反应是形象但不完全准确。后来自己拿COMSOL做光子能带计算越做越觉得这个类比其实抓到了很核心的东西——周期性排列的介电结构对不同频率的光子来说就是一套“能不能通行”的筛选规则。你调的晶格常数、层厚比例、材料折射率对比都是在给这个收费站设杆、定价、排车道。所谓一维光子晶体就是折射率只在一个方向上做周期性变化。最典型的形态是两种不同介电常数的材料交替堆叠比如空气和硅一层一间地排。光沿垂直层面向的方向入射在每个界面上都会发生部分反射一系列反射光之间会产生干涉。当每个界面的反射相位满足特定条件时所有反射波相干叠加某个频率段的光就会被“锁”在结构里传不出去只能被反射回原介质。这就是光子带隙的雏形也是分布式布拉格反射镜能反射窄带光、能当激光腔镜的根本原因。但这里有个容易混淆的地方一维光子晶体算能带算的并不是“一束光打上去的反射率”而是“无限周期结构内部允许存在哪些本征模式”。换句话说能带图里每条曲线上的一个点都代表一个能在周期介质中稳定传播的电磁本征态。而那些不存在解的频率范围就是带隙光子在这个区域里没有可传播模式结构对它就表现为“禁行”。这个思路和半导体物理里的电子能带完全平行只是把薛定谔方程换成了麦克斯韦方程。真正动手用COMSOL之前建议先把几个基本设计参数定下来否则后面参数化扫描会乱成一团晶格常数 a一个周期的总长度也就是相邻同种层的间隔。填充比 f高折射率层厚度在单个周期中的占比f d_high / a。高折射率层材料这里选硅近红外波段折射率 n ≈ 3.42介电常数 ε ≈ 11.68。低折射率层材料空气 n ≈ 1.0或者在绝缘体上硅平台上就是 SiO₂n ≈ 1.45。工作波段由 a 决定。一维光子晶体的能带结构满足标度律把 a 和波长同时放大或缩小相同倍数归一化能带完全不变。我的推荐组合是 a 600 nmf 0.5也就是硅层 300 nm、空气层 300 nm。在这个尺子下第一带隙的中心频率对应的归一化频率 ωa/2πc 大约在 0.25 附近换算成真实波长在 1.2 μm 左右。做近红外实验的组用这个参数起步很顺手。确定了这些下一步就是建模。下面我会按COMSOL实际操作顺序展开从几何到求解器一项一项说最后再把我踩过的坑单独列一节。如果你刚接触光子晶体仿真建议从头顺序做一遍不要直接跳到自己“以为会”的步骤。2. COMSOL建模仿真全流程参数、几何、边界条件与网格的逐步配置2.1 选择模型维度二维模型是一维光子晶体的最佳打开方式你可能想问一维光子晶体用一维几何模型不就行了吗确实可以但我不建议一上来就做纯 1D。原因很现实COMSOL 的电磁波模块在处理 1D 模型时虽然也能算本征频率但后处理能看的东西太少了——你很难直接看到本征模电场分布确认模式是不是横电模还是横磁模也费劲排查错误时两眼一抹黑。我习惯做一个 2D 模型几何上取一个包含完整周期的原胞x 方向是层叠方向y 方向取一段有限的长度。这样既能通过周期性边界条件模拟无限周期又能清清楚楚看到每个 k 点上的场形图判断模式性质。二维并不会显著增加计算量因为整个 y 方向是均匀的网格可以拉得很疏自由度往往比你想的要低。在 COMSOL 里新建模型时选择“模型向导”空间维度选“二维”添加物理场时选“波动光学模块 电磁波、频域”。这个物理场默认求解的是全矢量麦克斯韦方程适用于无源介质结构。研究类型先放在“特征频率”后面生成能带时再加“参数化扫描”。这里的“特征频率”求解的是电场本征方程解出来的特征值就是角频率的平方 ω²。COMSOL 的本征频率理所当然地输出为 Hz而我们做能带分析时习惯用归一化频率这个转换要留到后处理阶段做不要在求解器里手忙脚乱地改单位。2.2 几何与参数定义把设计变量写进全局定义后面会谢天谢地进入“全局定义 参数”表把下面这些变量先全部写进去。这一步看起来简单但直接用数字硬编码几何尺寸是新手最容易犯的错误一旦后面要扫 f 或者改 a你会在几何、材料、网格三个阶段来回改烦到怀疑人生。参数名表达式说明a600[nm]晶格常数f0.5高折射率层填充比d_sif*a硅层厚度d_gap(1-f)*a空气间隙厚度y_max0.5*a计算域 y 方向半高取整数个周期即可n_si3.42硅折射率eps_sin_si^2硅相对介电常数n_air1.0空气折射率然后在几何里画一个矩形宽度设为 a高度设为 2*y_max这就是原胞整体。再画一个硅层矩形宽度设为 d_si高度同样为 2*y_max让它左边界对齐到整体矩形的左边界。最后用布尔分割把整体矩形分成两个域左边是硅右边是空气。这样做的好处是网格和材料指派都各归各域边界处不会因为几何重叠产生缝隙或重复面。COMSOL 的几何精度默认会用绝对容差处理布尔操作600 nm 这种纳米尺度下理论模型和几何之间的长度单位没有问题但记得在“几何 单位”里确认长度单位设成了 nm。如果你保留默认的 m输入 600[nm] 后模型里显示的几何大小会小到看不见初学者经常在这里浪费半小时屏幕上一片空白还以为是建模失败。2.3 材料指派与物理场细节介电常数对比是能带的生死线在“材料”节点里新建两个空材料一个命名“Si”相对介电常数设为 eps_si另一个命名“Air”相对介电常数设为 1。把材料指派到对应的域上。这里容易踩的第一个坑COMSOL 电磁波频域物理场默认会检测“折射率”或“相对介电常数”但有些版本的材料库会要求填完整的色散数据否则提示空材料不可用。如果你只是想快速验证能带直接自己定义空材料填一个常数介电值就够了不必强行调用材料库中带色散的硅模型。原因是本征频率计算是基于时谐方程色散会让每个频率点上的介电常数不同复杂度上升不少而能带结构本身最核心的物理由介电常数对比决定常数近似足以把带隙位置和宽度算得很准。物理场设置里电磁波频域的默认域方程是∇ × μ_r⁻¹ ∇ × E - k₀² ε_r E 0其中 k₀ ω/c 是真空波数。COMSOL 在这个物理场里默认电场变量为 E特征频率研究下求解的是使上述方程存在非零解的频率。不需要额外添加电流源或端口激励纯本征问题。2.4 周期性边界条件和布洛赫波矢能带图的灵魂所在能带计算的关键不在于域内方程而在于边界条件。一维光子晶体在 x 方向无限延伸但计算机只能处理有限域所以 x 方向两条边界必须用 Floquet 周期性边界条件连接起来模拟无限周期的效果。COMSOL 里在“电磁波、频域”下右键添加“周期条件”类型选“Floquet 周期”。在“周期边界条件”设置里需要指定边界对源边界选左侧 x 0目标边界选右侧 x a。角色默认都是“源”和“目标”不能选反。接下来的关键项是布洛赫波矢 k第一个 Floquet 周期化方向填入kx第二个方向填0这里 kx 是一个全局参数稍后参数化扫描时会逐点赋不同的值。布洛赫定理决定了周期性边界上电场的相位关系E(xa) E(x) · exp(-i·kx·a)kx 从 0 扫到 π/a对应第一布里渊区的边界点。由于是一维结构能带只在 x 方向这个路径上有色散y 方向不设周期条件直接将上下边界设为“完美磁导体”或“完美电导体”即可两者都不会在 y 方向引入额外的带隙。交替用这两种边界条件各跑一遍可以快速区分横电模和横磁模横电模电场沿 y 方向在上边界为完美磁导体时更自然而横磁模磁场沿 y 方向在完美电导体下贡献不同解。如果你想让设置更严谨也可以让 y 方向也用周期性边界条件但把第二个波矢分量子设为 0本质效果一样。真正的区别是周期边界允许 y 方向的 evanescent 场形态更完整可带来的自由度增加会影响求解器性能。我对这个具体结构用 PEC/PMC 做验证结果差异在 0.1% 以内可以放心省算力。网格剖分上每个介电层内至少要有 3~5 层六面体或四边形单元。在 COMSOL 中可以用“映射”网格x 方向设 20 个单元y 方向设 4 个单元总共 80 个矩形单元自由度非常小但精度已经足够好。我试过把 x 方向单元数加倍到 40带边频率偏移不到 0.05%这就是网格收敛了。后面我在单独一节会讲网格不足时你会看到什么“假带隙”这里先不展开。3. 从波矢扫掠到能带图参数化扫描、数据整理与带隙定标3.1 参数化扫描到底扫的什么第一布里渊区路径能带图上的横轴是波矢 k纵轴是频率 ω。所谓扫掠就是把之前定义的 Floquet 波矢 kx 从一个值变到另一个值在每一个值处各求一次特征频率。这一系列特征频率点对应连起来就形成能带曲线。在 COMSOL 中“研究 1”的类型选特征频率。在“研究 特征频率”设置里把“所需模式数”设成 8也就是每个 k 点求最低的 8 个本征模式。想一想我们关心的是低频段的前几条带足够了。如果设太多模式高频模式往往包含不希望的奇异解反而干扰识别带边。然后右键“研究 1”添加“参数化扫描”。扫描参数选全局定义的 kx扫掠列表里新建一个等差数列。kx 取值范围是 0 到 π/a。输入的时候要注意COMSOL 的扫描列表支持直接写表达式例如range(0, pi/(6*a), pi/a)也就是从 0 开始每次递增 π/(6a)一直取到 π/a一共 7 个 k 点。想要分辨率高一点就把分母的 6 改成 12能带曲线会更平滑但求解时间线性增加。我一般先用 7 个点跑通全流程确认没有建模错误再加密到 13 或 25 个点出图。细心的人会注意到这里只扫了第一布里渊区的一半即 0 → π/a。为什么不全扫 -π/a → π/a原因是一维系统具有时间反演对称性能带关于 k 0 对称也就是 ω(k) ω(-k)。你扫一半再镜像对称展开和全扫得到的是同一张图。少扫一半计算量减半这才是真正理解布里渊区压缩的好处而不只是照着教程抄。3.2 求解器选择与求解结果速度不算关键稳最重要特征频率问题的求解器默认会用 MUMPS 或 PARDISO。COMSOL 默认选择通常够用但在参数化扫描多个 k 点时我建议手动固定用 MUMPS原因是它在处理小规模三维自由度问题时鲁棒性更好遇到部分模式收敛失败的概率略低。打开“求解器配置 特征值求解器”把线性求解器改为“MUMPS”容差保持默认。扫描跑完后你可以先直接看图。COMSOL 会自动生成一个“全局一维图组”把特征频率随 kx 变化的散点画出来。但由于 kx 是参数默认绘图横轴是参数索引而不是波矢值看起来会很别扭。更好的做法是“派生值 全局计算”把每个 kx 下的 8 个特征频率导成表格。在“全局计算”的设置里表达式写freq也就是 COMSOL 内部的特征频率变量单位是 Hz。勾选“从参数化扫描中评估所有参数”然后点击计算。得到的数据是一个 7 行 8 列的表格每一行对应一个 kx 扫描点每一列对应模式序号。把这个表格拷贝到 Excel 或者 MATLAB 里作为绘图的原始数据。我能给你的最实在建议是绘图阶段不要用 COMSOL 内置的全局绘图功能直接导出图因为它默认对每个参数点只显示最新一组值散点连线的逻辑也经常跟你要的不一致。把这个表导出来自己写个十几行脚本横轴是 kx纵轴是 freq画 8 条不同颜色或符号的散点线一条条能带就出来了彻底排除 COMSOL 后处理的黑盒干扰。3.3 归一化频率与带隙提取仿真值怎么跟论文对标光子晶体领域发论文、对比实验、对照已有结构几乎所有人都在用归一化频率而不是绝对频率。原因是绝对频率和晶格常数直接挂钩同样的带隙设计a 是 600 nm 还是 1500 nm绝对频率差很多倍但结构本质完全一样。归一化的做法是f_norm f · a / c也可以用波长关系写成f_norm a / λ₀所以当你算出某个带边频率 f只要把它乘以晶格常数 a除以真空光速 c就得到归一化频率。以 a 600 nm 为例若你在某个 k 点求出本征频率 f 250 THz则 f_norm 250e12 × 600e-9 / 3e8 0.5。这个数直接可以和文献上的能带图对比单位是“无量纲”但大家习惯标成 ωa/2πc。从能带图上找带隙方法是看同一 k 值下相邻能带之间是否存在空隙特别是布里渊区边界 k π/a 附近。因为周期性结构中带隙最容易在能带折叠处张开。一维光子晶体典型的第一带隙通常是从第一条带的上边沿到第二条带的下边沿在 k π/a 处这两条带之间有一段空白范围那个范围就是你结构的带隙。带隙宽度用相对宽度表示Δω / ω₀其中 ω₀ 是带隙中心频率。用硅/空气 f 0.5 的一维结构第一带隙相对宽度大约能到 25% 上下这已经算很可观了。想加宽带隙核心手段是提高两种材料的介电常数对比想让带隙往高频挪就整体缩小 a想让带隙位置在更宽频带里交织出多段就调整填充比使其偏离 0.5。4. 能带图会骗人模式识别、透射谱交叉验证与缺陷态拓展4.1 每条带上的“点”都是谁电场场形图告诉你答案能带图给出了一堆离散点但如果你只看图而不去看每个点的场分布你根本不知道这些模式长什么样。我强烈建议你在布里渊区高对称点k 0 和 k π/a处各打开 2~3 个本征模的电场模图。操作方法是在“结果 电场模”里画图然后把“参数选择”固定到某个 k 点再通过模式序号切换查看不同带上的解。你会看到两类明显不同的场形一种电场集中分布在硅层内部另一种电场主要待在空气层中。前者类似于经典导波模中的“高折射率层受限模”后者则是空气层占优的模式。带隙形成的直观解释是当某个频率的光想要在结构中传播时它必须同时匹配这两种空间分布的交替要求而这个匹配条件在某个频段内无法满足于是传播被禁止。这里要特别提醒有些模式是数值上解析的“伪解”通常出现在计算域边界上形状特殊比如场的极大值完全束缚在 PEC 边界上。判断标准很简单——看它的场是否能自然延伸到域内部以及在加密网格后该频率是否几乎不变。伪解对网格极其敏感加密后频率漂移几个百分点甚至消失。遇到这种模式直接忽略不用在能带图上画出来。4.2 交叉验证用透射谱互相对照比上十次反射率都管用能带计算得到的带隙要和透射谱结果彼此印证。单独一张能带图只能说明“这个频率没有周期介质中的本征解”但它没有直接回答“如果光从外部入射会被反射多少”。两者虽然物理上相关工程验证时最好两条腿走路。验证思路同样保持原胞结构把 x 方向两端各加上一段匀质硅材料作入射/出射通道再在两端设置端口边界条件用频域研究计算 S 参数。然后在带隙频率范围内扫频如果你测到透射率 S21 呈深凹陷跟能带图上带隙范围完全吻合那说明能带计算没有出错。在 COMSOL 中实现时几何上做成一个长的计算域入射通道硅、N 个周期比如 10 个、出射通道硅两端用散射边界条件或端口边界。在“电磁波频域”研究中频率扫描范围设成覆盖在带隙附近的一段——比如从 150 THz 到 300 THz。计算透射率曲线后你会发现下降沿和带边位置几乎没有偏移。如果相差很大那基本可以断定是建模时周期边界或材料指派出了错。这里必须补充一句能带计算是理想无限周期结构的结果而透射谱对应的是有限个周期的结构。周期数少比如少于 5 个透射谱的禁带凹陷会比较浅不明显周期数多20 个以上透射率在带隙内接近零下降沿也更陡。如果只做 3 个周期就发现带隙没有完全闭合这属于物理正常现象不是错误。4.3 打破周期在结构中插入缺陷层看看带隙里如何长出“允许态”如果只学了完整的能带图就收手那确实只玩了一半。一维光子晶体最有工程价值的地方在于周期被破坏后带隙内会局域一个缺陷模式。它的物理图像是周期结构本应禁止频率 f 传播但你在其中插入一个厚度和介电常数与周期层不同的缺陷层相当于收费站车道里开了一个“特批通道”f 这个频率的光就可以非常局域地被束缚在缺陷处。这在光学微腔、滤波器、激光器中都有广泛应用。操作上把之前的 2D 几何扩展为超胞比如取 5 个完整周期再把中间那个周期里的硅层厚度改成 1.5 倍。这样超胞是周期性的吗严格说不是因为缺陷打破了周期。如果你想在一个有限结构里求本地模式只需在超胞两端也加上周期性边界条件——周期变成 5a——然后重新跑特征频率。解出来的模式会观察到带隙内出现一个孤立的频率点对应场能量高度集中在缺陷层中两侧周期结构里场的衰减极快。把这个模的 Q 值算出来工程滤波器设计的第一步就迈出去了。缺陷层的厚度偏离 0.5 的倍数越多缺陷模在带隙中偏离中心越远向一个带边靠近。这是一个非常好用的设计旋钮我实测过厚度从 1.0 倍增加到 1.5 倍缺陷模从带隙中心移动到接近下带边附近频率移动非常灵敏。5. 跑了十几次才总结出来的坑网格、模式识别、参数范围5.1 网格太小“画不出带”网格太细“画出假模式”这个坑几乎每个新手都会踩一次。网格太疏时高折射率层内部的电磁场变化根本分辨不出来算出来的特征频率整体偏高带边位置完全不准。我用 a 600 nm 的硅/空气结构测过x 方向每个硅层里只放 2 个网格单元时带隙上下边比加密后的结果偏了超过几个百分点放到 8 个单元频谱接近收敛。但网格太密也不是什么好事。在 2D 模型里自由度上升到几十万以后MUMPS 对内存的占用会陡然上升不少情况下会出现部分高频模式不收敛表现为警告信息后丢解。我从实用角度建议先粗网格跑通、确认物理过程再局部加密。尤其要注意的是在硅/空气界面附近网格允许跨越材料边界生成三角形单元但我用映射网格把两种介质分得干干净净边界处频率反而更稳定。瞬态压力峰值解耦缝隙都解决了。判断网格是否收敛的一个简单标准把单元数加倍重新算最低几个模式频率相对变化小于 0.1%就认为这个网格量是可信的。挨个 k 点检测一遍最好只在高对称点抽查不需要全部扫描点都做。5.2 模式编号不对齐交叉识别省一天时间参数化扫描时你可能会发现一个尴尬问题在 k 点 1 的第 3 个模式到了 k 点 2 可能跑到第 4 个位置因为两条带在色散中发生了交叉。如果直接按模式序号对相邻 k 点连线你会画出一些莫名其妙的折线——能带似乎突然断开又重新接上其实只是模式顺序发生了交换。解决办法有两个。第一在每个 k 点访问 eigensolution 时同时导出该模式的电场分布手动判断它属于哪条带。第二用 COMSOL 的“特征频率”研究里自带的模式追踪功能某些版本支持或者在后处理时按电场能量的分布特征来重新排序。笨办法也很有效只取最低 4 条带画图时交叉概率低取 8 条带交叉基本必然出现但低 4 条带混序的风险仍可控只要你的 k 点间隔足够密。间隔太疏时两条带夹成锐角的区域会出现“鸡生蛋”式的顺序错位这也是为什么我不建议 k 点取得太少。5.3 参数范围设错有些“带隙”只是你在边界上人为造出来的最后一个坑比较隐蔽。如果 Floquet 周期性边界里 kx 扫描范围设错了或者在研究设置里的特征频率范围过窄你可能会在目标频段内漏算模式导致能带图看起来有一条很宽的“假带隙”——其实只是你没有让求解器去那里找解。在 COMSOL 特征频率研究中“所需模式数”设置为 8 意味着求解器会求最低的几个特征值但实际求解范围受绝对容差和初始猜测影响。想更稳妥的话可以把模式数增加到 12并把“搜索频率的参考点”设在低频区附近。这之后你会发现大部分模式都重复了——同一结构的高频模式收敛过程会反复初始化不用担心只取前几条有效。还有一点是关于 y 方向边界条件。如果你在 y 方向使用 PEC那么某些模式会被强制为零边界这会导致模式缺失。你可以同时跑一次 PEC 和 PMC 两种边界条件把两次的结果合并就能把横电横磁全拿到。如果不做这一步能带图中你少了一张关键的模式分支怎么分析都会觉得缺了什么。回到文章开头那个“高速公路收费站”的比喻——现在你应该有更深的体感了一维光子晶体的能带图就是收费站的通行规则表COMSOL 就是那张规则表背后的“交通调度模拟器”。你调整晶格常数、层厚比例就像是在改收费站的限高杆和车道数量。这篇文章算是我自己从零跑通这个仿真流程的完整总结按这个顺序做下来你不仅能得到一条漂亮的能带图还能搞明白图里每根线每段空隙到底在说什么。之后再去碰二维三维光子晶体、拓扑光子结构这类更复杂的方向就不会再被建模仿真本身绊住脚了。 一维光子晶体COMSOL能带仿真：从建模到带隙分析
1. 一维光子晶体到底在算什么收费站模型和它的物理边界第一次听到“一维光子晶体像光子的高速公路收费站”这个说法是在一次组会上。我当时的反应是形象但不完全准确。后来自己拿COMSOL做光子能带计算越做越觉得这个类比其实抓到了很核心的东西——周期性排列的介电结构对不同频率的光子来说就是一套“能不能通行”的筛选规则。你调的晶格常数、层厚比例、材料折射率对比都是在给这个收费站设杆、定价、排车道。所谓一维光子晶体就是折射率只在一个方向上做周期性变化。最典型的形态是两种不同介电常数的材料交替堆叠比如空气和硅一层一间地排。光沿垂直层面向的方向入射在每个界面上都会发生部分反射一系列反射光之间会产生干涉。当每个界面的反射相位满足特定条件时所有反射波相干叠加某个频率段的光就会被“锁”在结构里传不出去只能被反射回原介质。这就是光子带隙的雏形也是分布式布拉格反射镜能反射窄带光、能当激光腔镜的根本原因。但这里有个容易混淆的地方一维光子晶体算能带算的并不是“一束光打上去的反射率”而是“无限周期结构内部允许存在哪些本征模式”。换句话说能带图里每条曲线上的一个点都代表一个能在周期介质中稳定传播的电磁本征态。而那些不存在解的频率范围就是带隙光子在这个区域里没有可传播模式结构对它就表现为“禁行”。这个思路和半导体物理里的电子能带完全平行只是把薛定谔方程换成了麦克斯韦方程。真正动手用COMSOL之前建议先把几个基本设计参数定下来否则后面参数化扫描会乱成一团晶格常数 a一个周期的总长度也就是相邻同种层的间隔。填充比 f高折射率层厚度在单个周期中的占比f d_high / a。高折射率层材料这里选硅近红外波段折射率 n ≈ 3.42介电常数 ε ≈ 11.68。低折射率层材料空气 n ≈ 1.0或者在绝缘体上硅平台上就是 SiO₂n ≈ 1.45。工作波段由 a 决定。一维光子晶体的能带结构满足标度律把 a 和波长同时放大或缩小相同倍数归一化能带完全不变。我的推荐组合是 a 600 nmf 0.5也就是硅层 300 nm、空气层 300 nm。在这个尺子下第一带隙的中心频率对应的归一化频率 ωa/2πc 大约在 0.25 附近换算成真实波长在 1.2 μm 左右。做近红外实验的组用这个参数起步很顺手。确定了这些下一步就是建模。下面我会按COMSOL实际操作顺序展开从几何到求解器一项一项说最后再把我踩过的坑单独列一节。如果你刚接触光子晶体仿真建议从头顺序做一遍不要直接跳到自己“以为会”的步骤。2. COMSOL建模仿真全流程参数、几何、边界条件与网格的逐步配置2.1 选择模型维度二维模型是一维光子晶体的最佳打开方式你可能想问一维光子晶体用一维几何模型不就行了吗确实可以但我不建议一上来就做纯 1D。原因很现实COMSOL 的电磁波模块在处理 1D 模型时虽然也能算本征频率但后处理能看的东西太少了——你很难直接看到本征模电场分布确认模式是不是横电模还是横磁模也费劲排查错误时两眼一抹黑。我习惯做一个 2D 模型几何上取一个包含完整周期的原胞x 方向是层叠方向y 方向取一段有限的长度。这样既能通过周期性边界条件模拟无限周期又能清清楚楚看到每个 k 点上的场形图判断模式性质。二维并不会显著增加计算量因为整个 y 方向是均匀的网格可以拉得很疏自由度往往比你想的要低。在 COMSOL 里新建模型时选择“模型向导”空间维度选“二维”添加物理场时选“波动光学模块 电磁波、频域”。这个物理场默认求解的是全矢量麦克斯韦方程适用于无源介质结构。研究类型先放在“特征频率”后面生成能带时再加“参数化扫描”。这里的“特征频率”求解的是电场本征方程解出来的特征值就是角频率的平方 ω²。COMSOL 的本征频率理所当然地输出为 Hz而我们做能带分析时习惯用归一化频率这个转换要留到后处理阶段做不要在求解器里手忙脚乱地改单位。2.2 几何与参数定义把设计变量写进全局定义后面会谢天谢地进入“全局定义 参数”表把下面这些变量先全部写进去。这一步看起来简单但直接用数字硬编码几何尺寸是新手最容易犯的错误一旦后面要扫 f 或者改 a你会在几何、材料、网格三个阶段来回改烦到怀疑人生。参数名表达式说明a600[nm]晶格常数f0.5高折射率层填充比d_sif*a硅层厚度d_gap(1-f)*a空气间隙厚度y_max0.5*a计算域 y 方向半高取整数个周期即可n_si3.42硅折射率eps_sin_si^2硅相对介电常数n_air1.0空气折射率然后在几何里画一个矩形宽度设为 a高度设为 2*y_max这就是原胞整体。再画一个硅层矩形宽度设为 d_si高度同样为 2*y_max让它左边界对齐到整体矩形的左边界。最后用布尔分割把整体矩形分成两个域左边是硅右边是空气。这样做的好处是网格和材料指派都各归各域边界处不会因为几何重叠产生缝隙或重复面。COMSOL 的几何精度默认会用绝对容差处理布尔操作600 nm 这种纳米尺度下理论模型和几何之间的长度单位没有问题但记得在“几何 单位”里确认长度单位设成了 nm。如果你保留默认的 m输入 600[nm] 后模型里显示的几何大小会小到看不见初学者经常在这里浪费半小时屏幕上一片空白还以为是建模失败。2.3 材料指派与物理场细节介电常数对比是能带的生死线在“材料”节点里新建两个空材料一个命名“Si”相对介电常数设为 eps_si另一个命名“Air”相对介电常数设为 1。把材料指派到对应的域上。这里容易踩的第一个坑COMSOL 电磁波频域物理场默认会检测“折射率”或“相对介电常数”但有些版本的材料库会要求填完整的色散数据否则提示空材料不可用。如果你只是想快速验证能带直接自己定义空材料填一个常数介电值就够了不必强行调用材料库中带色散的硅模型。原因是本征频率计算是基于时谐方程色散会让每个频率点上的介电常数不同复杂度上升不少而能带结构本身最核心的物理由介电常数对比决定常数近似足以把带隙位置和宽度算得很准。物理场设置里电磁波频域的默认域方程是∇ × μ_r⁻¹ ∇ × E - k₀² ε_r E 0其中 k₀ ω/c 是真空波数。COMSOL 在这个物理场里默认电场变量为 E特征频率研究下求解的是使上述方程存在非零解的频率。不需要额外添加电流源或端口激励纯本征问题。2.4 周期性边界条件和布洛赫波矢能带图的灵魂所在能带计算的关键不在于域内方程而在于边界条件。一维光子晶体在 x 方向无限延伸但计算机只能处理有限域所以 x 方向两条边界必须用 Floquet 周期性边界条件连接起来模拟无限周期的效果。COMSOL 里在“电磁波、频域”下右键添加“周期条件”类型选“Floquet 周期”。在“周期边界条件”设置里需要指定边界对源边界选左侧 x 0目标边界选右侧 x a。角色默认都是“源”和“目标”不能选反。接下来的关键项是布洛赫波矢 k第一个 Floquet 周期化方向填入kx第二个方向填0这里 kx 是一个全局参数稍后参数化扫描时会逐点赋不同的值。布洛赫定理决定了周期性边界上电场的相位关系E(xa) E(x) · exp(-i·kx·a)kx 从 0 扫到 π/a对应第一布里渊区的边界点。由于是一维结构能带只在 x 方向这个路径上有色散y 方向不设周期条件直接将上下边界设为“完美磁导体”或“完美电导体”即可两者都不会在 y 方向引入额外的带隙。交替用这两种边界条件各跑一遍可以快速区分横电模和横磁模横电模电场沿 y 方向在上边界为完美磁导体时更自然而横磁模磁场沿 y 方向在完美电导体下贡献不同解。如果你想让设置更严谨也可以让 y 方向也用周期性边界条件但把第二个波矢分量子设为 0本质效果一样。真正的区别是周期边界允许 y 方向的 evanescent 场形态更完整可带来的自由度增加会影响求解器性能。我对这个具体结构用 PEC/PMC 做验证结果差异在 0.1% 以内可以放心省算力。网格剖分上每个介电层内至少要有 3~5 层六面体或四边形单元。在 COMSOL 中可以用“映射”网格x 方向设 20 个单元y 方向设 4 个单元总共 80 个矩形单元自由度非常小但精度已经足够好。我试过把 x 方向单元数加倍到 40带边频率偏移不到 0.05%这就是网格收敛了。后面我在单独一节会讲网格不足时你会看到什么“假带隙”这里先不展开。3. 从波矢扫掠到能带图参数化扫描、数据整理与带隙定标3.1 参数化扫描到底扫的什么第一布里渊区路径能带图上的横轴是波矢 k纵轴是频率 ω。所谓扫掠就是把之前定义的 Floquet 波矢 kx 从一个值变到另一个值在每一个值处各求一次特征频率。这一系列特征频率点对应连起来就形成能带曲线。在 COMSOL 中“研究 1”的类型选特征频率。在“研究 特征频率”设置里把“所需模式数”设成 8也就是每个 k 点求最低的 8 个本征模式。想一想我们关心的是低频段的前几条带足够了。如果设太多模式高频模式往往包含不希望的奇异解反而干扰识别带边。然后右键“研究 1”添加“参数化扫描”。扫描参数选全局定义的 kx扫掠列表里新建一个等差数列。kx 取值范围是 0 到 π/a。输入的时候要注意COMSOL 的扫描列表支持直接写表达式例如range(0, pi/(6*a), pi/a)也就是从 0 开始每次递增 π/(6a)一直取到 π/a一共 7 个 k 点。想要分辨率高一点就把分母的 6 改成 12能带曲线会更平滑但求解时间线性增加。我一般先用 7 个点跑通全流程确认没有建模错误再加密到 13 或 25 个点出图。细心的人会注意到这里只扫了第一布里渊区的一半即 0 → π/a。为什么不全扫 -π/a → π/a原因是一维系统具有时间反演对称性能带关于 k 0 对称也就是 ω(k) ω(-k)。你扫一半再镜像对称展开和全扫得到的是同一张图。少扫一半计算量减半这才是真正理解布里渊区压缩的好处而不只是照着教程抄。3.2 求解器选择与求解结果速度不算关键稳最重要特征频率问题的求解器默认会用 MUMPS 或 PARDISO。COMSOL 默认选择通常够用但在参数化扫描多个 k 点时我建议手动固定用 MUMPS原因是它在处理小规模三维自由度问题时鲁棒性更好遇到部分模式收敛失败的概率略低。打开“求解器配置 特征值求解器”把线性求解器改为“MUMPS”容差保持默认。扫描跑完后你可以先直接看图。COMSOL 会自动生成一个“全局一维图组”把特征频率随 kx 变化的散点画出来。但由于 kx 是参数默认绘图横轴是参数索引而不是波矢值看起来会很别扭。更好的做法是“派生值 全局计算”把每个 kx 下的 8 个特征频率导成表格。在“全局计算”的设置里表达式写freq也就是 COMSOL 内部的特征频率变量单位是 Hz。勾选“从参数化扫描中评估所有参数”然后点击计算。得到的数据是一个 7 行 8 列的表格每一行对应一个 kx 扫描点每一列对应模式序号。把这个表格拷贝到 Excel 或者 MATLAB 里作为绘图的原始数据。我能给你的最实在建议是绘图阶段不要用 COMSOL 内置的全局绘图功能直接导出图因为它默认对每个参数点只显示最新一组值散点连线的逻辑也经常跟你要的不一致。把这个表导出来自己写个十几行脚本横轴是 kx纵轴是 freq画 8 条不同颜色或符号的散点线一条条能带就出来了彻底排除 COMSOL 后处理的黑盒干扰。3.3 归一化频率与带隙提取仿真值怎么跟论文对标光子晶体领域发论文、对比实验、对照已有结构几乎所有人都在用归一化频率而不是绝对频率。原因是绝对频率和晶格常数直接挂钩同样的带隙设计a 是 600 nm 还是 1500 nm绝对频率差很多倍但结构本质完全一样。归一化的做法是f_norm f · a / c也可以用波长关系写成f_norm a / λ₀所以当你算出某个带边频率 f只要把它乘以晶格常数 a除以真空光速 c就得到归一化频率。以 a 600 nm 为例若你在某个 k 点求出本征频率 f 250 THz则 f_norm 250e12 × 600e-9 / 3e8 0.5。这个数直接可以和文献上的能带图对比单位是“无量纲”但大家习惯标成 ωa/2πc。从能带图上找带隙方法是看同一 k 值下相邻能带之间是否存在空隙特别是布里渊区边界 k π/a 附近。因为周期性结构中带隙最容易在能带折叠处张开。一维光子晶体典型的第一带隙通常是从第一条带的上边沿到第二条带的下边沿在 k π/a 处这两条带之间有一段空白范围那个范围就是你结构的带隙。带隙宽度用相对宽度表示Δω / ω₀其中 ω₀ 是带隙中心频率。用硅/空气 f 0.5 的一维结构第一带隙相对宽度大约能到 25% 上下这已经算很可观了。想加宽带隙核心手段是提高两种材料的介电常数对比想让带隙往高频挪就整体缩小 a想让带隙位置在更宽频带里交织出多段就调整填充比使其偏离 0.5。4. 能带图会骗人模式识别、透射谱交叉验证与缺陷态拓展4.1 每条带上的“点”都是谁电场场形图告诉你答案能带图给出了一堆离散点但如果你只看图而不去看每个点的场分布你根本不知道这些模式长什么样。我强烈建议你在布里渊区高对称点k 0 和 k π/a处各打开 2~3 个本征模的电场模图。操作方法是在“结果 电场模”里画图然后把“参数选择”固定到某个 k 点再通过模式序号切换查看不同带上的解。你会看到两类明显不同的场形一种电场集中分布在硅层内部另一种电场主要待在空气层中。前者类似于经典导波模中的“高折射率层受限模”后者则是空气层占优的模式。带隙形成的直观解释是当某个频率的光想要在结构中传播时它必须同时匹配这两种空间分布的交替要求而这个匹配条件在某个频段内无法满足于是传播被禁止。这里要特别提醒有些模式是数值上解析的“伪解”通常出现在计算域边界上形状特殊比如场的极大值完全束缚在 PEC 边界上。判断标准很简单——看它的场是否能自然延伸到域内部以及在加密网格后该频率是否几乎不变。伪解对网格极其敏感加密后频率漂移几个百分点甚至消失。遇到这种模式直接忽略不用在能带图上画出来。4.2 交叉验证用透射谱互相对照比上十次反射率都管用能带计算得到的带隙要和透射谱结果彼此印证。单独一张能带图只能说明“这个频率没有周期介质中的本征解”但它没有直接回答“如果光从外部入射会被反射多少”。两者虽然物理上相关工程验证时最好两条腿走路。验证思路同样保持原胞结构把 x 方向两端各加上一段匀质硅材料作入射/出射通道再在两端设置端口边界条件用频域研究计算 S 参数。然后在带隙频率范围内扫频如果你测到透射率 S21 呈深凹陷跟能带图上带隙范围完全吻合那说明能带计算没有出错。在 COMSOL 中实现时几何上做成一个长的计算域入射通道硅、N 个周期比如 10 个、出射通道硅两端用散射边界条件或端口边界。在“电磁波频域”研究中频率扫描范围设成覆盖在带隙附近的一段——比如从 150 THz 到 300 THz。计算透射率曲线后你会发现下降沿和带边位置几乎没有偏移。如果相差很大那基本可以断定是建模时周期边界或材料指派出了错。这里必须补充一句能带计算是理想无限周期结构的结果而透射谱对应的是有限个周期的结构。周期数少比如少于 5 个透射谱的禁带凹陷会比较浅不明显周期数多20 个以上透射率在带隙内接近零下降沿也更陡。如果只做 3 个周期就发现带隙没有完全闭合这属于物理正常现象不是错误。4.3 打破周期在结构中插入缺陷层看看带隙里如何长出“允许态”如果只学了完整的能带图就收手那确实只玩了一半。一维光子晶体最有工程价值的地方在于周期被破坏后带隙内会局域一个缺陷模式。它的物理图像是周期结构本应禁止频率 f 传播但你在其中插入一个厚度和介电常数与周期层不同的缺陷层相当于收费站车道里开了一个“特批通道”f 这个频率的光就可以非常局域地被束缚在缺陷处。这在光学微腔、滤波器、激光器中都有广泛应用。操作上把之前的 2D 几何扩展为超胞比如取 5 个完整周期再把中间那个周期里的硅层厚度改成 1.5 倍。这样超胞是周期性的吗严格说不是因为缺陷打破了周期。如果你想在一个有限结构里求本地模式只需在超胞两端也加上周期性边界条件——周期变成 5a——然后重新跑特征频率。解出来的模式会观察到带隙内出现一个孤立的频率点对应场能量高度集中在缺陷层中两侧周期结构里场的衰减极快。把这个模的 Q 值算出来工程滤波器设计的第一步就迈出去了。缺陷层的厚度偏离 0.5 的倍数越多缺陷模在带隙中偏离中心越远向一个带边靠近。这是一个非常好用的设计旋钮我实测过厚度从 1.0 倍增加到 1.5 倍缺陷模从带隙中心移动到接近下带边附近频率移动非常灵敏。5. 跑了十几次才总结出来的坑网格、模式识别、参数范围5.1 网格太小“画不出带”网格太细“画出假模式”这个坑几乎每个新手都会踩一次。网格太疏时高折射率层内部的电磁场变化根本分辨不出来算出来的特征频率整体偏高带边位置完全不准。我用 a 600 nm 的硅/空气结构测过x 方向每个硅层里只放 2 个网格单元时带隙上下边比加密后的结果偏了超过几个百分点放到 8 个单元频谱接近收敛。但网格太密也不是什么好事。在 2D 模型里自由度上升到几十万以后MUMPS 对内存的占用会陡然上升不少情况下会出现部分高频模式不收敛表现为警告信息后丢解。我从实用角度建议先粗网格跑通、确认物理过程再局部加密。尤其要注意的是在硅/空气界面附近网格允许跨越材料边界生成三角形单元但我用映射网格把两种介质分得干干净净边界处频率反而更稳定。瞬态压力峰值解耦缝隙都解决了。判断网格是否收敛的一个简单标准把单元数加倍重新算最低几个模式频率相对变化小于 0.1%就认为这个网格量是可信的。挨个 k 点检测一遍最好只在高对称点抽查不需要全部扫描点都做。5.2 模式编号不对齐交叉识别省一天时间参数化扫描时你可能会发现一个尴尬问题在 k 点 1 的第 3 个模式到了 k 点 2 可能跑到第 4 个位置因为两条带在色散中发生了交叉。如果直接按模式序号对相邻 k 点连线你会画出一些莫名其妙的折线——能带似乎突然断开又重新接上其实只是模式顺序发生了交换。解决办法有两个。第一在每个 k 点访问 eigensolution 时同时导出该模式的电场分布手动判断它属于哪条带。第二用 COMSOL 的“特征频率”研究里自带的模式追踪功能某些版本支持或者在后处理时按电场能量的分布特征来重新排序。笨办法也很有效只取最低 4 条带画图时交叉概率低取 8 条带交叉基本必然出现但低 4 条带混序的风险仍可控只要你的 k 点间隔足够密。间隔太疏时两条带夹成锐角的区域会出现“鸡生蛋”式的顺序错位这也是为什么我不建议 k 点取得太少。5.3 参数范围设错有些“带隙”只是你在边界上人为造出来的最后一个坑比较隐蔽。如果 Floquet 周期性边界里 kx 扫描范围设错了或者在研究设置里的特征频率范围过窄你可能会在目标频段内漏算模式导致能带图看起来有一条很宽的“假带隙”——其实只是你没有让求解器去那里找解。在 COMSOL 特征频率研究中“所需模式数”设置为 8 意味着求解器会求最低的几个特征值但实际求解范围受绝对容差和初始猜测影响。想更稳妥的话可以把模式数增加到 12并把“搜索频率的参考点”设在低频区附近。这之后你会发现大部分模式都重复了——同一结构的高频模式收敛过程会反复初始化不用担心只取前几条有效。还有一点是关于 y 方向边界条件。如果你在 y 方向使用 PEC那么某些模式会被强制为零边界这会导致模式缺失。你可以同时跑一次 PEC 和 PMC 两种边界条件把两次的结果合并就能把横电横磁全拿到。如果不做这一步能带图中你少了一张关键的模式分支怎么分析都会觉得缺了什么。回到文章开头那个“高速公路收费站”的比喻——现在你应该有更深的体感了一维光子晶体的能带图就是收费站的通行规则表COMSOL 就是那张规则表背后的“交通调度模拟器”。你调整晶格常数、层厚比例就像是在改收费站的限高杆和车道数量。这篇文章算是我自己从零跑通这个仿真流程的完整总结按这个顺序做下来你不仅能得到一条漂亮的能带图还能搞明白图里每根线每段空隙到底在说什么。之后再去碰二维三维光子晶体、拓扑光子结构这类更复杂的方向就不会再被建模仿真本身绊住脚了。