今天下午把 COMSOL 里的扭转双层光子晶体能带和透射谱跑通了。严格说吧计算本身不算难难的是怎么把一个“扭了一下”的周期性结构老老实实描述给 COMSOL 听。整个过程让我想起之前调魔角石墨烯模型时的狼狈角度一变超胞尺寸跟着变边界条件、网格、特征值排序全都得跟着重新捋一遍。不过搞完之后再看结果确实有意思——单层光子晶体原本干净的带隙加了一个转角之后能带开始折叠、平带冒出来透射谱也出现了窄尖峰这感觉和“两层石墨烯转个角度开始整活”确实有那么点神似。这篇就完整记录一下我的建模思路、COMSOL 操作流程、关键参数的设定以及过程中踩过的一堆坑。内容会比较长但对想复现扭转光子晶体、莫尔超胞能带、透射谱这波操作的人应该能直接拿去当操作手册用。1. 动手前先想清楚扭转光子晶体到底在算什么1.1 先理清物理图像能带、带隙、透射谱各自是什么光子晶体的能带结构本质上是介电常数呈周期性变化的材料中电磁波本征频率与波矢量之间的关系。频率和波矢组成一条条曲线叫色散曲线也就是日常说的能带。如果在某个频段内结构里找不到任何满足布洛赫条件的本征模式这个频段就对应带隙。带隙的物理结果非常直接这个频段的光无法穿过结构透射谱上会对应一段深谷。透射谱则是更贴近实验的物理量。把一束平面波打在有限厚度的周期结构上在另一侧统计透过的电磁功率就得到透射率随频率变化的曲线。带隙在透射谱上是一个凹陷区域带隙边缘附近往往会有明显的透射增强和振荡。我这次做的事可以拆成两半。第一半靠本征值求解得到能带图第二半靠频域求解得到透射谱。然后重点观察引入扭转角之后能带和透射谱发生了什么变化。很多人可能觉得这类计算在光学教材里已经写过无数遍但加上“扭转”两个字后事情就完全不一样了。1.2 为什么扭转结构值得折腾从魔角石墨烯谈起魔角石墨烯在凝聚态物理里的地位不用多说两层石墨烯相对旋转大约 1.1°费米面附近的能带变得非常平电子之间强烈的相互作用把原本普通的金属变成了超导态。核心机制不在材料本身而在扭转产生的莫尔条纹周期尺度远超原胞尺寸导致低能区出现了平带。光子晶体可以理解成电磁波版本的周期势场。如果把光子晶体的概念推广到双层结构用旋转角作为自由度会出现很多类似现象两层各自的布拉格散射条件不再一致带隙位置随转角漂移带边可能出现平坦带慢光效应、局域共振随之而来。这思路在近几年的电磁超表面、双层金属网、双层介质光栅里已经有大量应用本质上都是在用几何转角调制波的传播。我做这件事的初衷不是原创研究而是想把文献里二维莫尔光子晶体模型在 COMSOL 里从头到尾复现一遍验证能带与透射谱之间的自洽性。做完之后发现COMSOL 卡脖子的地方往往不是物理接口而是几何与周期性条件的组织方式。1.3 扭转后真正的难点公共周期性和莫尔超胞单层光子晶体能带的计算完全标准化取一个原胞加周期性边界条件沿布里渊区扫描波矢。但双层结构扭转之后两层晶格的平移对称方向错开了系统不再具有单层原胞尺度的周期性。要严格使用布洛赫定理必须找两层晶格共同的最小周期单元也就是莫尔条纹决定的超胞。超胞尺寸很不稳定。扭转角越小莫尔周期通常越大。比如转角几度时超胞可能包含数十上百个原胞计算量迅速膨胀。我第一次想到“用 COMSOL 算扭转光子晶体”时第一反应也是几何和网格怎么画得动。关键一点COMSOL 不会自动帮你寻找超胞必须自己根据两层晶格的配准关系先确定合理的扭转角再构建同时容纳两层周期性的最小矩形区域。这一步看似简单却几乎决定了后面所有计算能不能跑完。我后面会专门讲具体做法因为这一步跳过的人基本都会在边界条件上翻车。2. 建模思路与方案选型二维模型、材料参数、超胞确定2.1 二维模型还是三维模型我这次选择二维电磁波模型。理由很实际二维模型计算量小、设置直观而且只要介质结构沿厚度方向均匀面内能带物理不会因为维度的简化而失真。真实的三维双层光子晶体通常是在两个平面上分别制备薄膜中间有间隔层严格仿真需要完整建立 z 方向结构。那种模型网格量很容易翻一个数量级以上周期性边界条件也要考虑膜厚方向的模式。对复现型工作我建议先用二维模型把现象跑出来确认结果没问题再升级到三维。很多文献里的莫尔光子晶体理论模型本来就是在二维框架下推导的二维计算既能抓住主要物理又能避免把时间浪费在几何网格优化上。如果你想后续做实验对比再考虑三维也不迟。2.2 材料参数与几何参数要怎么定我用的是正方形晶格介质柱光子晶体。晶格常数 a 取 1 μm介质柱半径 r 取 0.2a材料用硅相对介电常数 11.7背景是空气。这个组合在二维电磁波仿真里有非常成熟的参考结果单层带隙大概落在归一化频率 0.45 附近方便验证。模型里最难处理的不是“哪层是硅、哪层是空气”而是扭转之后两层介质柱坐标的叠加。我强烈建议把两层的介电位移看成空间分布叠加底层按底层晶格排列顶层按旋转后的坐标排列两者在 xy 平面上叠加成最终等效介电常数分布。这样可以避开 z 方向物理上的歧义也符合二维莫尔光子晶体的常见处理方式。介质柱重叠的区域在等效模型里直接按硅处理不会造成几何布尔运算的麻烦。如果你非要每个圆柱单独建模再处理交叉重叠网格会非常难控制收敛性也会变差。下面的参数表是我这次实际使用的直接可以参考参数数值说明a1 μm晶格常数r0.2a介质柱半径ε_si11.7硅的相对介电常数ε_bg1.0背景空气介电常数θ12.5°最终选定的扭转角L7.2a莫尔超胞边长f_mod0.3 ~ 0.8归一化频率扫描范围2.3 扭转角与莫尔超胞的确定方法这一步值得单独展开。寻找公共超胞的常规做法是给底层晶格一组格点给顶层施加旋转算子然后在有限范围内搜索哪些角度下两套格点经过合适的整数平移能够重合。这个搜索用 Python 脚本几十行就能完成不需要在 COMSOL 里做。我最后选定的扭转角是 12.5° 左右对应的莫尔超胞边长大约是 7.2a。在这个超胞里大约包含上百根介质柱。选择这个角度是为了平衡物理效果和计算量角度太大莫尔周期不明显角度太小超胞大到网格无法承受。脚本输出的超胞尺寸要手动写进 COMSOL 的参数表。后续几何、网格、周期性边界条件全部围绕这个超胞尺寸展开。这里必须强调超胞的边长直接影响布里渊区大小直接影响波矢换算公式一定先记录清楚否则后面边界条件会乱套。3. 能带计算实操先从单层原胞跑通再上双层扭转3.1 单层原胞能带计算的标准流程在 COMSOL 中新建二维模型物理场接口选择“电磁波频域”研究类型选择“特征值”。几何画成一个边长为 a 的正方形中心放一个半径为 r 的圆柱。左右两边加周期性条件上下两边加周期性条件波矢分别设为 kx 和 ky。周期性条件的波矢单位是 rad/m需要自己换算。方晶格布里渊区是正方形高对称路径取 Γ(0,0) → X(π/a,0) → M(π/a,π/a) → Γ(0,0)。把每一段路径的参数化关系提前写好然后用参数化扫描求解特征值。特征值求解器里有个关键设置叫“搜索基准点”。我一般把基准点设在预计的带隙中心附近让求解器优先取到该频率附近的本征模。需要求解的特征值个数设在 8 到 15 之间。数量太少高能分支会漏掉数量太多求解时间成倍增加高阶模的收敛性也容易变差。3.2 双层扭转超胞的几何构建一旦单层流程跑通就进入核心部分。根据第 2.3 节确定的超胞边长 L在 COMSOL 中画一个 L×L 正方形。按底层晶格坐标放置一组介质柱再按旋转后的顶层坐标放置另一组介质柱。这里强烈推荐用代码批量生成坐标而不是在 GUI 一个个点。把所有圆心坐标写成一个 CSV 文件在 COMSOL 的“几何”节点里导入。Python 脚本里要特别注意顶层坐标必须预先完成旋转矩阵变换旋转中心也需要与超胞中心保持一致否则导入后位置会整体偏移。几何生成后把所有圆和矩形超胞做布尔运算保留落在超胞内部的圆。材料设置很简单超胞背景设为空气圆区域设为硅。对二维等效模型来说重叠区域按硅处理不会引入额外问题。3.3 网格划分与求解器配置超胞里有上百根半径为 0.2 μm 的介质柱如果不控制网格自由剖分会在圆形边界生成大量细碎网格整体节点数轻松突破百万。我的经验是最大单元尺寸设 a/6圆边界用“边的最大单元尺寸”约束到 0.1 μm 左右其余交给默认的自由三角形。特征值问题规模在几十万自由度时PARDISO 或 MUMPS 都能胜任。参数扫描会反复求解很多次建议把“特征值搜索基准点”固定而不是依赖上一个参数的结果。特征值问题本身存在谱排序问题参数变化后基准点固定能减少模式之间的混乱。网格粗的时候能带曲线很光滑但带隙位置有偏移加密之后曲线可能出细锯齿那是数值色散效应。我是这样确定网格尺度的连续加密两次如果带隙位置变化小于 0.5%就认为当前网格足够可信。4. 透射谱计算实操让光真实穿过扭转结构4.1 透射模型怎么搭起来能带图算的是本征模式还需要验证实际照射平面波后的透射率。因此要在频域里重新建一个模型。这次不用超胞本征问题而是把一个横向宽度为一个超胞宽的带状区域作为模拟域沿 x 方向依次是空气缓冲区、扭转光子晶体层、空气缓冲区。左右边界用“端口”定义入射端口在左侧输出端口在右侧。上下边界用周期性边界条件连接模拟横向无限周期排列。入射方向选正入射波矢为零。端口会自动计算功率和 S 参数透射率直接看 S21。这个模型的最大优点是能直接看到电磁场分布。带隙频率下场会被拦截在结构内部带隙外能看到透射波从右端口溢出。把透射谱和能带图放在一起对照物理图像会非常清晰。我建议跑完频域计算后至少截取三个频率的场图带内、带边、带外各一张。4.2 频率扫描的求解设置研究选择“频域”频率范围覆盖能带图上带隙附近的频段。我的归一化频率范围选 0.3 到 0.8对应实际频率约 9×10^13 到 2.4×10^14 Hz。扫描点先设 200 个如果结果出现很尖锐的共振峰再在峰附近局部加密。扫频用直接线性求解器比迭代求解器更稳定。PARDISO 多线程模式非常适合多核 CPU。端口边界上要确认“激励”只开放左侧端口否则右侧端口也在馈能透射谱会出现奇怪的干涉条纹。如果透射率在某些频段出现周期性的 Fabry-Perot 振荡不要急着怀疑物理模型。那通常是因为两端缓冲区不够长结构自身反射在端口之间形成了微腔。把缓冲区加长到至少两个工作波长振荡幅度会明显降低。4.3 能带图和透射谱怎么互相对照我习惯把能带图画在左侧透射谱画在右侧再用虚线标出带隙上下沿。如果计算正确带隙频率范围内透射率接近零带边附近出现明显透射振荡。振荡对应能带中靠近带边的模式在 k 空间里的耦合行为。加入扭转后最有趣的现象是带隙内部出现了一段孤立的平带对应透射谱中出现窄尖峰。这个尖峰不是普通 Fabry-Perot 共振而是莫尔超胞内的局域模式。透过率在窄频带里接近 100%偏离该频率后迅速跌落。我验证过这一点现象非常直观也是整个项目里最让人兴奋的部分。5. 那些我踩过的坑与排查方案5.1 周期性边界条件的相位设错这是新手最容易犯的错误。COMSOL 的周期性条件需要指定相位差而相位差是波矢乘以边界周期长度不是直接把波矢填进去就完事。正方形原胞边长 a波矢 kx 对应的相位差应该是 kx·a。换成超胞时这个 a 要换成超胞边长 L。我之前跑出来的能带明显不对称检查半天才发现相位写错。另外扫描路径经过 M 点时kx 和 ky 要同步更新不能只改其中一个。最稳妥方式是提前写好参数化的波矢公式把路径分成三段每一段分别定义 kx、ky 关于扫描参数的表达式。5.2 特征值排序混乱导致能带断线能带图最常见的翻车场景是曲线断断续续交叉处连不起来。原因是 COMSOL 在参数扫描时每个参数点的特征值独立求解没有自动模式追踪。两条能带接近交叉时求解器可能把模式顺序弄乱曲线自然就断了。我的解决方法是固定“特征值搜索基准点”让模式尽量沿着连续的频率分支走。如果还乱就缩小参数扫描步长再用后处理按频率连续变化重排。手动重排非常痛苦所以我写了一个脚本按上一个参数点的频率最近邻原则把当前点的各个模式重新分配。5.3 超胞网格太多算不动超胞边长从 7.2 μm 增加到十几微米时网格数量会迅速膨胀。如果还按原胞网格参数自由剖分可能直接内存溢出。建议先做网格敏感性分析找到能保证带隙位置稳定的最大单元尺寸。对光子晶体通常 a/8 到 a/10 已经足够更细的网格只拉慢速度不会让物理结论发生质变。不要在整个域统一加密。只在介质柱附近用细化区域介质柱之间的空白空间用粗网格。COMSOL 的自适应网格细化在特征值问题里也能用但强烈建议关掉因为参数扫描中每次重新剖分成本太高。5.4 透射谱的毛刺和噪声透射谱出现大量高频毛刺通常是网格太粗糙。尤其是圆形介质柱边界如果节点太少入射波与边界相互作用会被严重离散化产生虚假散射共振。处理办法在圆边界上显式设置更大数量的边单元例如每个圆圆周至少 32 个节点。如果毛刺只在很窄频带出现另一个嫌疑是端口反射。让左侧缓冲区长度大于至少两个波长并把端口激励设为平面波而不是默认高斯束能有效消除端口附近的近场耦合。纯平面波入射配合周期边界是最严格也最干净的方式。5.5 算力不够时的降级方案真的算不动时最直接的降级方案是把圆柱换成方形柱。方形柱在物理上接近圆形柱但几何简单得多网格数量能省下一大截。另一个方案是暂时选一个小角度下较小的超胞先把流程跑通确认没有 bug再放大超胞。能带求解本身不是扫描扭转角而是扫描波矢。扭转角固定、波矢扫描的耗时大概每小时几十个点。如果要把扭转角也作为参数一起扫描建议把任务拆分到多个工作进程或者直接上远程服务器做批处理。6. 一些能让你少走弯路的实操技巧6.1 一定要先做单层结果验证在折腾扭转结构之前先用单层原胞把带隙位置算出来跟文献或经验值对比。如果单层都算不对双层结果没有参考价值。这一步看起来多余却能在排查问题时省下数小时。我的工作流是先跑单层原胞确认带隙位置再跑单层超胞确认超胞能带是单层能带折叠后的结果最后加入第二层。每一步都对比上一步任何误差都能尽早暴露。6.2 用脚本辅助生成莫尔坐标COMSOL 的 GUI 里手动画多个圆既慢又容易出错。强烈建议用 Python 生成所有圆心坐标再批量导入。编写脚本时要注意两套坐标必须在同一个全局坐标系下顶层坐标预先完成旋转矩阵变换。生成超胞时还要去重两个圆心距离小于 0.01a 的直接丢掉否则布尔运算会出现退化边界网格报错。6.3 能带图和透射谱要统一归一化单位能带图纵轴通常用归一化频率 f·a/c透射谱横轴是实际频率 Hz。在同一篇文章里对比时建议把透射谱横轴也换算成归一化频率否则看图还得心算单位容易把带隙位置对应错。归一化频率里的 a 在扭转结构中到底取单层晶格常数还是超胞边长也要统一。我的习惯是统一用单层晶格常数 a这样能带图能直接与单层光子晶体结果对比莫尔超胞带来的能带折叠效果也更容易看出来。6.4 用参数化扫描观察能带随扭转角的变化如果建模时把扭转角定义为全局参数在能带计算时可以加一层参数化扫描。这种做法计算量很大但结果特别像“魔角”故事有些角度下带隙宽度突然压缩有些角度下带隙中央出现几乎水平的亮线。我在实际跑的时候会把每个角度下的能带频率导出成矩阵再用 Python 画成“频率-角度-能带”的二维伪彩图。这样能很直观地看出哪个角度出现最平坦的带也就是“光子魔角”的位置。6.5 结果导出比截图更重要COMSOL 的绘图窗口虽然漂亮但对后续分析和写报告不够用。我习惯在参数扫描结束后把所有特征频率按参数点导出成文本文件并把 S21 数据也导出到同一个工作表。接下来用 Python 或 MATLAB 统一画图风格一致修改也方便。最后分享一个我自己的小习惯计算完成后把能带图、透射谱、某个典型频率下电场分布图放在同一张画布里。这样拿到结果的人能一眼看出哪些频率能透过、哪些被挡下、空间局域模式长什么样。一张好的汇总图有时候比一整页推导更有说服力。7. 后续还能怎么扩展这个模型做完之后至少还有三个方向可以继续整活。第一个方向是扫描不同扭转角下的透射谱找到最佳慢光频率然后设计一个“可调谐光子晶体滤波器”通过控制转角改变器件的工作频率。第二个方向是研究莫尔超胞里的拓扑边界态把具有拓扑性质的模式引进来看看在转角下是否还存在边缘态。第三个方向是升级到三维模型把两层光子晶体放在不同 z 平面中间加介质间隔层和实验结构完全对齐。如果只是想快速复现二维模型已经足够。如果想发表论文或者做深入科研强烈建议在二维结果基础上再做实验对比。无论做哪个方向最重要的事情都是先把 COMSOL 的能带、透射谱计算流程打磨到熟练。这一步扎实了后面所有扩展都会顺畅很多。