做光子晶体和超表面仿真的人几乎都绕不开一个基础操作算能带。COMSOL里跑正方晶格光子晶体的能带表面上看就是几何、材料、本征频率三件事实际趟下来发现坑不少——Floquet周期边界里波矢分量的映射、特征值数量的选择、后处理时频率归一化的写法任何一个地方不对画出来的图都不可信。这篇文章把我自己常用的完整流程拆开讲从参数设计、Floquet边界设置到波矢扫描和后处理全部是可以直接复现的关键步骤。适合刚上手光子晶体仿真、想快速搭一套能带计算模板的人如果你已经算过几条带但总感觉哪里不对这里面整理的坑也值得对一遍。1. 先想明白能带仿真到底在算什么又为什么要沿路径扫1.1 从布洛赫定理到本征频率方程光子晶体是介电常数在空间周期排布的结构在正方晶格里背景和介质柱的折射率按晶格常数 a 的周期重复。光在这种周期性介质里的行为和电子在晶体周期势场中的行为很像稳态场解可以写成布洛赫形式 E(r)e^{ik·r}u_k(r)其中 u_k(r) 和周期结构同周期。代入麦克斯韦方程组后每个波矢 k 对应一系列离散的本征频率这些本征频率的集合就是能带。用 COMSOL 求解能带本质是在给定几何、材料和边界条件下求解一个本征值问题把 k 当作已知参数把频率当作特征值。对二维正方晶格、TM 偏振电场沿 z 方向的情形方程可以简化成标量的亥姆霍兹方程数值求解比较稳定。周期条件和布洛赫定理通过 Floquet 周期边界条件引入在单胞边界两侧场的相位差等于波矢 k 与晶格矢量的点积。所以只要改变 Floquet 条件里的 k 分量就可以扫描动量空间中的不同点得到一系列本征频率把这些点按顺序连起来就是能带图。这里有个容易混淆的地方k 不是求解出来的是设定好的。COMSOL 里并不需要把布洛赫形式的解析解写进方程你只需要在周期条件的波矢分量里填一个变量然后让这个变量做参数扫描。也就是说k 的作用完全通过边界条件体现而内部的求解区域只需要一个单胞。1.2 正方晶格的布里渊区与Γ-X-M-Γ路径正方晶格的第一布里渊区是一个以原点为中心的正方形倒格子基矢大小为 2π/a。因为结构的对称性全布里渊区的色散关系可以从不可约布里渊区中的高对称路径复制出来。正方晶格常用的是 Γ-X-M-Γ 这条闭合路径Γ 点k(0,0)对应长波极限也是布里渊区中心X 点k(π/a, 0)第一布里渊区边界的边中点M 点k(π/a, π/a)第一布里渊区的角点对于 COMSOL 参数化扫描最省事的写法不是直接扫 kx 和 ky 两个参数而是定义一个路径参数把 Γ-X、X-M、M-Γ 三段映射到三段连续的参数区间。常见做法是设一个 k_scan从 0 扫到 3其中 0 到 1 对应 Γ 到 X1 到 2 对应 X 到 M2 到 3 对应 M 到 Γ。这样扫描一次其实就在整个高对称边界上走了一遍。物理上为什么要沿这条路径扫因为能带中的极值和带隙边缘通常出现在高对称点附近沿路径扫能快速掌握带隙的大致位置和宽度尤其是判断“这个结构有没有光子带隙”时绝大多数情况下用这条路径就足够了。如果后续需要更严格地确认全布里渊区没有残留的态可以再加密扫描整个二维 k 网格但一般教程级分析首选路径扫描。知道原理再看界面里的选项就不会出现“为什么这里要填 kx、ky填了有什么意义”的困惑。接下来要做的都是在模型层面把这几件事落实。2. 几何、材料与周期边界最容易埋雷的三个地方2.1 单胞几何怎么画最顺手新建一个二维组件长度单位保持米。先把全局参数写在最前面这样后面所有表达式都能引用同一套参数改起来方便参数推荐值说明a1e-6 m晶格常数典型近红外光子晶体尺度r0.25*a介质柱半径先按经典带隙较宽的位置起手n_rod3.5柱体折射率对应硅在近红外波段n_bg1背景折射率空气几何部分只画两个对象一个边长 a 的正方形一个圆心在正方形中心、半径为 r 的圆。正方形就是单胞边界圆代表介质柱区域。COMSOL 默认会把重叠区域自动分割成不同域并不需要手工做布尔运算再切回去你只要确保圆完全在正方形内部柱体和边界不相交即可。材料设定上建议直接在“材料”节点里手动输入折射率不要依赖内置材料库里的完整光学模型。内置材料库通常带有色散和损耗虽然更真实但对能带本征值问题来说实折射率模型才是默认的标准做法损耗会让特征值变复数反而干扰对带结构的判断。我的习惯是新建一个空材料把折射率实部填成 3.5另一个填成 1两个域各分配一个。2.2 材料参数背后折射率对比才是核心能带结构最敏感的参数不是绝对折射率而是折射率对比度。柱体和背景的折射率差越大布拉格散射越强带隙通常越宽。n_rod3.5、n_bg1 的对比度已经相当高非常适合教学和原型验证。如果你想算得更准确可以把 3.5 换成硅在具体工作波段的实验值比如 3.45。这里给出一个常见的材料参考区间柱体/孔材料折射率近红外适用场景硅 Si3.45~3.5最常用带隙宽工艺成熟砷化镓 GaAs3.3~3.4有源光子晶体、集成光源二氧化钛 TiO22.4~2.6可见光波段空气孔1.0背景高折射率介质中的孔阵记住一个换算关系COMSOL 光学模块里常用相对介电常数 ε_r它和折射率 n 的关系是 ε_rn²。所以填 n3.5 时也可以直接在材料属性里填相对介电常数 12.25效果一致。有的教程里让你填介电常数有的让你填折射率本质是一样的。2.3 Floquet周期条件配置细节与常见错误在建好几何和材料之后添加物理场“电磁波频域”。先设置求解区域为背景域和柱体域整体。二维 TM 模式默认求的是 Ez 分量。然后添加两个周期条件左边界和右边界配对上边界和下边界配对在周期条件中把类型改成“Floquet”展开“Floquet 周期”设置在波矢 k 分量的表达式里填上变量名 kx 和 ky。这两个变量不是固定值而是在全局定义里写好的表达式后面用参数扫描驱动。注意两对边界的“源”和“目标”对应关系要求同向边匹配。如果左右边界设反等效波矢方向会反号本征频率本身受影响通常不大但后续看模式场分布会完全对不上。新手最容易犯的错是两对边界选反或者只有一对边界设了 Floquet另一对忘设导致模型变成波导型而不是周期型模式能带图会多出很多莫名其妙的模式。检查方法很简单算完一个 k 点后看 Ez 模式的相位分布是否满足周期延拓关系。另一个常见问题是把 kx 和 ky 填成两个独立的固定值然后两个方向都设同一个变量。正方晶格需要分别指定x 方向边界用 kx 表达式y 方向边界用 ky 表达式。3. 物理场与研究设置完整流程把k_scan变成能带3.1 研究类型本征频率而非频域扫描这点很关键很多人刚开始容易把能带计算当成频域扫描在边界激发一个宽频脉冲然后看透射谱那不是能带仿真。能带计算要的是“给定 k 下系统有哪些谐振频率”所以研究要选“特征值”研究在物理场变量中选择“本征频率”。在特征值研究设置里最核心的是两个参数搜索基准点和所需特征值数。搜索基准点通常设 0也就是从零频率开始向上找所需特征值数建议先设 10 到 12。设置太少的话高频带扫到一半会缺失设置太多则计算时间明显变长。我一般从 10 个起手看前 5 条带是否完整再决定要不要加。求解器方面默认特征值求解器通常够用。如果出现收敛不稳定可以在求解器配置里把线性求解器换成 MUMPS或者把特征值搜索方式从“围绕基准点”改成“最大最小实部”但多数场景不需要手动干预。真正需要主动做的是把搜索基准点设在关心频率范围附近的实数比如 0让求解器围绕低频段找全。3.2 波矢扫描用分段函数把三条路径压进一个参数在全局参数里先定义 k_scan初始值设为 0。然后再定义 kx 和 ky 两个变量用 if 嵌套实现三段路径映射kx if(k_scan1, k_scan*pi/a, if(k_scan2, pi/a, (3-k_scan)*pi/a)) ky if(k_scan1, 0, if(k_scan2, (k_scan-1)*pi/a, (3-k_scan)*pi/a))这三段逻辑很简单k_scan 在 0 到 1kx 从 0 线性长到 π/aky 保持 0对应 Γ→Xk_scan 在 1 到 2kx 固定在 π/aky 从 0 长到 π/a对应 X→Mk_scan 在 2 到 3kx 和 ky 同时从 π/a 缩到 0对应 M→Γ然后在研究中加一个参数扫描节点选择 k_scan取值范围用 range(0,0.02,3)也就是 0 到 3 之间每隔 0.02 取一个点共 151 个点。步长越小曲线越平滑刚上手时可以先取 0.05 试跑确认流程通后再加密。这个分段表达式是整个能带扫描最核心的一步也是最容易被教程当作“魔法”带过的地方。很多案例直接用两个参数 kx 和 ky 做二维扫描然后后处理里再挑选路径上的点。那也能算但会产生大量无关的 k 点白白浪费计算量。用分段函数把路径编码进一个参数是更高效的折中。3.3 特征值排序与求解器细节参数扫描里每求解一个 k 点COMSOL 都会返回一组特征值。这些特征值默认并不保证永远是按频率升序排列尤其当两条带交叉或简并时模式序号会跳变。这是能带图断裂的根源之一不是算错了是排序不稳。为了让结果更容易整理求解器配置里通常优先使用默认特征值求解器。如果出现收敛不稳定可以尝试把线性求解器从默认的 SPOOLES 换成 MUMPS。计算量上做个估算151 个 k 点每个点求解 10 个特征值二维三角形网格几千到几万个自由度单核跑大约几分钟如果开启并行参数扫描通常十几秒到一两分钟就能完成。没必要为了省时间把网格设得太粗后面对网格收敛性有具体说明。4. 后处理从本征频率到一张规范能带图4.1 在COMSOL里快速画能带曲线参数扫描跑完后打开结果节点新建一个一维绘图组。在绘图数据里把“参数选择”设为“所有参数”或当前扫描然后添加一个“全局”绘图表达式选择特征值频率 freq。x 轴数据设为 k_scan这样每条本征频率会随着 k_scan 画成一条线。但直接画出来的 freq 是绝对频率单位 Hz数字大且不直观。光子晶体能带图习惯用归一化频率 fa/c也就是把频率乘以晶格常数再除以光速。把这个表达式写在 y 轴数据里freq*a/c_constc_const 是 COMSOL 内置光速常数可以直接引用。如果想在 COMSOL 内部带标注可以在“一维绘图组”里点“x 轴数据”选“表达式”填 k_scan然后在“轴标注”里手动写 Wave vector再在图形上手动放注释。反过来我更推荐导出后用外部工具画标注方便得多。4.2 导出数据到Python或OriginCOMSOL 自带的图形适合快速确认趋势但正式出图我习惯导出成 CSV用 Python 或者 Origin 画。操作路径是“派生值”→“全局计算”表达式选 freq在“数据”里选择全部扫描参数计算后表格会把每个 k_scan 点对应的所有特征值以多列形式列出来。右键导出为 CSV之后处理。Python 侧的核心思路很简单import pandas as pd import matplotlib.pyplot as plt df pd.read_csv(bands.csv) a 1e-6 # 晶格常数单位m c 3e8 # 光速单位m/s x df.iloc[:, 0] # 第一列是 k_scan for i in range(1, df.shape[1]): y df.iloc[:, i] * a / c # 归一化频率 fa/c plt.plot(x, y, lw1.2) plt.xlabel(Wave vector) plt.ylabel(Normalized frequency $fa/c$) plt.xticks([0, 1, 2, 3], [$\\Gamma$, $X$, $M$, $\\Gamma$]) plt.show()注意很多 CSV 的列名是自动生成的 freq1、freq2 之类直接用列号访问更省事。x 轴上 0、1、2、3 分别对应 Γ、X、M、Γ 四个位置这样图片出来就是标准的能带图。Origin 用户可以在导入时把行数设成 k 点数量把每一列拖入折线图x 轴用 k_scan 列同样可以。4.3 怎么看懂能带图带隙、简并与模式分布能带图画出来之后第一件事是找带隙在某个归一化频率范围内整条路径上都没有能带经过这就叫完全带隙。完全带隙对应光子晶体里该频率的光被禁止传播。第二个要观察的是简并点尤其在 X 点或 M 点附近两条带会碰到一起这是晶格对称性导致的不是数值误差。如果参数选择不当比如网格太粗简并处可能被错误地劈开看起来像有一个假的极小缝。判断真假的方法是将网格加密一倍看简并点附近的频率差是否趋于零。第三个要对比的是 TM 和 TE 模式。介质柱型正方晶格在小半径下 TM 的带隙往往比 TE 明显这是定性规律。如果你想看 TE 模式需要把物理场因变量从 Ez 切成 Hz 重新算一遍不能和 TM 混在一张图里。5. 常见问题与避坑清单照着排查能省半天时间5.1 曲线断裂、缺失模式和排序混乱上面提到模式排序不稳定是最常见的问题。现象是能带曲线在某个 k 点附近突然跳格或者一条带消失了、另一条带却多出一截。处理办法主要是增加特征值数量因为如果搜索到的模式数少于实际存在的模式数排序就会出问题。另一个办法是导出数据后在外部按每个 k 点对频率做一次升序排列再逐带连线。这不会改变物理结果只是把可视化修正。还有一种情况某个 k 点返回的特征值数量明显比其他点少那通常是特征值求解器在搜索基准点附近漏掉了模式尤其是简并模式。把基准点稍微抬高或者把所需特征值数从 10 加到 14往往就解决了。5.2 Floquet边界配错和幽灵模式幽灵模式一般长这样能带图里出现一些非常平、频率非常低、或者和相邻带重叠得很奇怪的带。出现原因不外乎三种周期边界没配对、边界源目标选错方向、或者网格太粗导致局部谐振模混进来。排查顺序建议检查两对 Floquet 边界是否都设置了且 kx、ky 分别正确看几组典型 k 点的 Ez 分布确认周期延拓的相位关系把网格加密一倍重新算如果幽灵带消失或明显移动说明是网格问题提示判断简并真假的最快方法就是网格加密。伪简并通常一加密就被劈开真实简并则始终保持重合。5.3 网格收敛性、扫描速度与r/a优化对能带仿真来说网格密度直接决定高频带精度。建议初始网格最大单元尺寸设为 a/10。如果对高频带感兴趣再全局细化为 a/15并对比两次结果中你关心的带隙边缘频率变化是否小于 1%。如果变化明显说明网格还不够。参数扫描的步长也值得控制。先粗扫 0.05确认大概带隙位置后再在带隙边缘附近做 0.005 的细扫能大幅减少总计算量。如果还要对比不同 r/a 的带隙可以把 r/a 设成另一个外层参数做嵌套参数扫描但注意计算量会成倍增加建议先试一两个 r/a 值不要一下扫一排。最后分享一个经验规律对于硅柱正方晶格r/a 在 0.2 到 0.3 之间通常能观察到较宽的 TM 带隙我习惯从 0.25 开始试带隙不明显就往大调带结构碎了就往小调。这是经验值不同背景折射率会有偏移。我在实际跑这一套流程时最深的体会是不要把时间耗在追求一次到位上。第一次跑通全流程用粗网格、少模式、大步长哪怕画出来的图很糙只要趋势对整个模型的骨架就已经成立了剩下的是在各自的物理问题里逐步细化。参数扫描一圈十几分钟就出结果这个迭代速度对做研究的人来说非常友好。还想分享一个判断模型是否正常的小技巧先在 Γ 点只算一次本征频率把前几个模式的频率记录下来和均匀介质色散关系做对比。Γ 点处周期结构退化为平板频率应该接近背景介质的自由空间模式线如果这里就对不上后面整条能带基本不用看。先把每个环节钉牢能带图画出来自然干净。