光子晶体仿真看起来门槛高实际上绝大多数时间都花在修正模型的细节上。我最早跟着论文复现二维空气孔光子晶体整整一周都在跟能带图里的锯齿较劲最后才发现只是材料介电常数虚部没有清零。为了彻底摆脱反复试错的局面我决定用COMSOL 5.6把《光子晶体》教材中的典型案例完整复现一遍。目前这套复现项目积累了40多个可直接运行的mph文件涵盖一维、二维、三维三种维度的光子晶体结构包括透射谱、反射谱、能带图、本征模场分布等完整输出。这篇文章既是对这套案例库的说明也是把复现过程中踩过的坑、总结的方法和调参经验一次性讲清楚。无论你是刚接触光子晶体仿真的研究生还是已经在算能带但经常对不上文献结果的工程师这套从一维到三维的完整链路都值得参考。1. 为什么我要在COMSOL 5.6里复现光子晶体案例1.1 从“看书懂”到“动手会”的鸿沟光子晶体相关书籍通常会把能带理论讲得很细但真正动手建模的时候你会发现书上的信息根本不够用。比如书里可能只写了晶格常数a600 nm、空气孔半径r180 nm却没有写清楚Floquet边界条件的波矢量该怎么设置也没有说明用TM模式还是TE模式计算。这些信息缺失导致一百个人能跑出一百种结果。所以我决定换一个思路不再零散地在网上找示例而是选一本案例最完整、参数最清晰的光子晶体专著把其中能复现的算例逐个做出来。所谓复现不是简单画出几何而是让计算得到的禁带位置、能带宽度、透射率曲线和书中的结果一致。有了这样一条校准线后续做新结构的时候我才有底气去改参数、换材料知道哪些结果是合理的哪些是模型出错了。1.2 案例库的组织方式与文件规范40多个mph文件如果随意堆在一起半年后连自己都找不到。我在整理案例库时采用“维度—结构类型—求解目标”的三级目录结构根目录下分1D、2D、3D三大类每一类再按照具体结构细分。比如2D目录下就有正方晶格空气孔、三角晶格空气孔、六角晶格介质柱、线缺陷波导、微腔谐振器等子目录。文件命名也有一套固定的规则。我会把结构类型、关键几何参数、折射率组合、计算模式四个要素写进文件名。例如2D_Triangular_r0.3_n2.34_TM_band.mph看到名字就知道这是三角晶格空气孔结构r/a0.3背景折射率2.34算的是TM模式能带。另外每个子目录里放一个README.md记录案例来源、章节号、预期结果和模型注意事项。这样即使隔几个月再打开也能快速恢复上下文。1.3 选择COMSOL 5.6的具体原因标题里的“Comsol56”指的就是COMSOL 5.6版本。我选择这个版本主要是看中它的波光学模块稳定性。5.6对Floquet周期边界、端口边界和散射边界的底层求解器做了不少优化特征值计算的收敛性比早期版本强很多。特别是对高介电常数对比度的结构早期版本经常出现“找不到模式”的提示5.6的容错明显变好。另一个重要原因是5.6的“数学模块”里弱形式PDE接口增强了非线性求解能力。后面我要专门讲到的“基于COMSOL弱形式方程求解色散光子晶体能带”正是依赖这个接口。早期版本在这个功能上偏弱很多频散材料模型需要额外写代码才能求解。5.6把弱形式的稳定性提升了一个台阶才让我能在原生界面里完成色散能带计算。2. 光子晶体仿真的三条核心方法论2.1 布里渊区、倒格子与k空间路径不管是几维结构光子晶体仿真都绕不开三个核心概念正格子、倒格子、布里渊区。正格子是你在COMSOL里画的周期性晶胞倒格子是与之对应的动量空间周期单元布里渊区则是倒格子的原胞。能带图描绘的就是本征频率在这个布里渊区内沿特定路径的变化。初学者最容易混淆的是几何尺寸和k点路径之间的关系。几何尺寸可以用实际晶格常数建模比如三角晶格a600 nm但能带计算里的k点必须沿布里渊区边界走。三角晶格的高对称路径是Γ(0,0) → M(0.5,0) → K(0.333,0.333) → Γ(0,0)这些坐标是无量纲的是相对于倒格子基矢的。如果直接把路径坐标当作普通变量输入COMSOL得到的能带会整体变形。我的处理方法是把倒格子基矢的换算关系直接写进模型的全局参数里。以三角晶格为例倒格子基矢长度为4π/(a√3)Floquet边界需要的波矢量分量定义为kx 4*pi/(a*sqrt(3)) * s1 ky -4*pi/(3*a) * s1这样扫描参数s1遍历0到1的区间k点就自动沿高对称路径移动。我现在做二维案例时已经把这套表达式做成了公共参数组复制到任意二维模型都能直接用只需根据晶格类型修改系数。2.2 Floquet周期边界条件的参数化Floquet边界条件是光子晶体能带计算的基石。它的作用是让晶胞两侧的电磁场满足一个相位关系相当于把无限周期结构的边界效应浓缩到一个单元里。在COMSOL中设置周期边界时类型必须选“Floquet周期性”不能选普通的“周期性”否则边界两侧的电场无法传播相移。边界条件界面里需要指定两个方向的波矢量分量kFloq1和kFloq2。这两个分量的单位是rad/m不是倒格子坐标。很多人在这一步出错是因为直接把k点坐标输入进去导致相位积累错误。正确做法是换算正方晶格中kFloq1 2π·kx/akFloq2 2π·ky/a三角晶格则需要考虑基矢夹角手动算出两个方向的投影系数。在参数化扫描过程中我会把kFloq1和kFloq2定义成全局参数的表达式然后用“辅助扫描”功能让k点连续遍历高对称路径。这里有个经验扫描点数量不是越多越好。我通常设置61个扫描点既能保证能带曲线平滑又不至于让求解时间翻倍。对于多维参数扫描COMSOL的“参数扫描”会为每一组参数完整求解一次扫描点过多时建议拆成两段执行方便中途查看结果。2.3 特征值求解器的目标设置与模式筛选特征值求解器是能带计算的核心引擎。COMSOL默认会计算“所需模式数”个最低频率的特征模但光子晶体往往需要特定频率区间的模式而不是最低的那几个。我的习惯是把“特征值搜索范围”设置为目标频段的1.5倍再通过模式序号和场分布图做筛选。具体参数方面我通常设置“所需模式数”为8到12搜索范围是[0, 2×f_max]。f_max是目标频率上限。如果范围太窄高频率模式会被漏掉如果范围太宽会混入无效的数值模式。求解完成后COMSOL会在日志中给出“拒收特征值”列表这些被剔除的模式往往暗示着数值伪模或边界设置问题值得仔细查看。3. 一维案例复现透射谱与一维禁带3.1 一维多层膜模型的几何参数设定一维光子晶体最常见的形式是交叠膜堆。我复现的一个典型算例是紫外波段的多层膜SiO2层与TiO2层交替排列厚度分别为95 nm和65 nm周期数10。在COMSOL里我选择用二维模型来搭建几何虽然结构是一维周期但二维模型能直观观察场分布也为后续斜入射计算留了余地。几何构建时我会先用一个矩形代表整个膜堆再用“分割面”功能按层厚切成一系列子域。如果每一层都建独立矩形再拼接后期改厚度会非常痛苦。分割面的操作在5.6里支持参数控制把层厚定义成全局参数后改一个数值整个几何自动更新。材料方面SiO2折射率设为1.46TiO2设为2.35特别注意要把材料属性里的损耗虚部清零否则禁带位置会偏移。3.2 端口、周期边界与入射波设置一维膜堆的透射和反射谱需要在结构两侧设置端口边界条件。COMSOL的“端口”特性支持多模式设置入射端口的模式类型要选“衍射级”端口宽度必须包含一个完整周期。如果端口宽度小于一个周期透射率曲线会出现莫名其妙的震荡这个问题非常隐蔽我调试了整整一天才找到原因。上下两侧需要设置Floquet周期性边界把x方向的周期落实到模型中。注意上下边界不能使用默认的PEC或者PBC否则会引入非物理的反射。频率扫描范围设置为320 nm到440 nm波长跑完结果后能看到反射谱在390 nm附近出现明显的禁带这是两材料界面布拉格反射最强烈的波长位置。3.3 结果对照与网格精度控制我最初的版本误差很大禁带边缘频率比书中值偏移了7%左右。排查后发现是网格太粗65 nm厚的薄层里默认网格只剖了一层单元边界处电磁场分布根本没解析出来。把网格最大单元尺寸调整为25 nm后偏差缩小到了0.8%这个精度已经满足大多数工程需求。这个坑让我养成了一个习惯所有薄膜结构每一层材料至少要跨4层网格。具体做法是使用“边界层网格”在每层材料介质界面处强制加密。多层结构用边界层网格增加的自由度很少但对能带位置的影响非常明显。可以说一维光子晶体仿真精度不够大概率是网格的问题而不是求解器或物理设置的问题。4. 二维案例复现能带结构中的TM/TE模式4.1 正方晶格与三角晶格的建模差异二维案例是这套案例库中数量最多的部分。二维结构既能展现周期结构的共性又能通过不同的晶格排列得到丰富的能带性质。正方晶格和三角晶格的差别不仅仅在几何排布上更关键的是它们的倒格子形状和高对称点路径完全不同。正方晶格的倒格子仍是正方布里渊区高对称路径为Γ-X-M-Γ三角晶格的倒格子是六角对称高对称路径为Γ-M-K-Γ。在几何搭建时正方晶格只需一个正方形晶胞x和y方向设两个周期性边界。三角晶格则必须用平行四边形晶胞两个基矢长度相等但夹角120°。这里有一个常见错误很多人直接在正方形外框里放一个圆形空气孔来模拟三角晶格这等于改变了晶格对称性算出的能带并不属于真正的三角晶格。我的标准做法是建一个平行四边形晶胞使用全局参数定义顶点坐标比如a600 nm、r180 nm然后利用三角函数关系算出平行四边形的斜边顶点。Floquet边界恰好映射两个基矢方向这样才保证计算模型的对称性正确。4.2 k路径扫描的参数化实现二维案例能带计算中k路径扫描是最容易出错也最耗时的环节。我在全局参数里定义了三段扫描变量s1、s2、s3分别对应Γ-M、M-K、K-Γ三段路径。每一段路径用线性插值把扫描参数映射到k点坐标。比如Γ-M段s从0到1kx从0映射到0.5ky从0映射到0这里的坐标都是相对于倒格子基矢的无量纲坐标。为了把三段路径拼接到一次研究中我会定义一个总的扫描变量s通过分段函数判断s落在哪一段再切换对应的k坐标表达式。这个写法看起来繁琐但在COMSOL中可以用“阶梯函数”或“if条件表达式”实现设置完成后整个能带扫描是一次性跑完的。4.3 参数扫描与禁带优化二维结构最有价值的应用就是禁带优化。比如设计工作在通信波长1550 nm附近的空气孔光子晶体晶格常数a和空气孔半径r是两个最关键的自由度。书里通常给一组基准参数但实际设计时需要扫描r/a比值。我在案例库中准备了两类参数扫描模型。第一类是固定a、扫描r观察TM模禁带宽度变化。典型结果是r/a从0.2增大到0.35时TM模禁带逐渐变宽超过0.4后禁带反而开始收窄因为空气孔之间的介质墙太薄高次模开始出现。第二类是固定r/a、整体缩放a观察归一化禁带位置的变化。这类模型用来验证光子晶体的缩放定律归一化频率a/λ基本保持不变这是周期结构设计的理论基础也是能带图与实验对照的关键参照。参数扫描时内存占用不小。我在32 GB内存的机器上一个三角晶格案例单次求解约2分钟20组扫描约40分钟。如果网格超过10万自由度建议用“辅助扫描”来代替“参数扫描”能大幅减少内存压力。我实测下来两种方式的精度差异可以忽略。5. 三维案例复现从几何搭建到资源调配5.1 木堆结构的几何布尔与域设置三维光子晶体案例中木堆结构非常经典。它由多层介电柱堆叠而成每层柱子方向旋转90度四层构成一个周期。在COMSOL中搭建木堆结构最让人头疼的是几何布尔运算后的材料域标记。多个柱子做布尔并集后COMSOL有时会把交叠区域识别成内部边界导致后续网格无法跨边界传播。我的处理方式是在布尔运算之前给每根柱子做“显式选择”布尔操作时选择“保留被选中的实体”把结构分为若干子域再逐个赋予材料。这样即使后续做参数扫描修改柱宽材料分配也不会被打乱。每个三维案例我都会记录几何构建顺序因为COMSOL的布尔运算是记录在模型树里的顺序错了后续修改几乎无法进行。5.2 网格策略与内存平衡三维能带计算对硬件的需求很高。木堆结构如果直接使用默认的四面体网格单个晶胞就需要大概80万到100万个自由度内存占用逼近16 GB。我的经验是两步走先用粗网格快速试算获得能带的大致位置和模式数量然后再用细化网格在目标频段精确计算。COMSOL 5.6的“自适应网格细化”在三维模型中很有效。它能自动识别场梯度大的区域在不增加整体网格数量的前提下修正局部精度。但自适应细化会增加迭代次数和应用时间所以只适合在最终求解阶段开启试算阶段要保持关闭。内存方面三维模型建议至少16 GB内存求解时开启多核并行速度提升非常明显。5.3 三维场分布的后处理技巧三维案例除了能带曲线通常还需要输出漂亮的场分布图。光子晶体场图能直观展示光与结构的相互作用线缺陷波导模式场集中在缺陷周围木堆结构的光子禁带模式场分布在介电柱之间的空隙里。在COMSOL里输出场图关键点是选对切面和位置。我常用的方法是用一个“工作平面”切过结构中间层叠加“高度图”显示场强配合透明显示介质结构。颜色表选择需要注意彩色映射在黑白打印时会失真论文投稿建议使用灰度或双色渐变。如果结果是复数场我会分别画实部和模值实部用于观察相位拓扑模值用于观察能量分布。这两种图配合起来才能判断模式是传播态还是局域态。6. 进阶弱形式方程求解色散光子晶体能带6.1 内置求解器为什么不够用前面几节的内容都是基于“电磁波频域”接口的“特征频率”研究。这个方法对线性、无频散、各向同性的介质完全够用。但遇到两类结构内置求解器就不太好办了第一类是色散材料比如金属、等离子体材料、增益介质它们的介电常数随频率变化第二类是各向异性材料比如磁性光子晶体本构关系是张量形式。“特征频率”研究的求解流程是固定频率后解出空间场。它需要先给定一个明确的介电常数值再算对应频率。如果介电常数本身就是频率的函数这个循环就变成了需要自洽求解的方程内置求解器难以处理。而弱形式方程可以做到这也是“基于COMSOL弱形式方程求解色散光子晶体能带”这个方向的初衷。6.2 弱形式方程的具体实现步骤弱形式的核心是把微分方程两边乘以一个任意检验函数再对整个求解域积分。以二维光子晶体TM模式为例控制方程是∇ × (1/ε_r(r) · ∇ × E_z) (ω²/c²) · E_z写成等效的弱积分表达式之后被积函数中包含两个部分第一项涉及电场梯度和检验函数梯度的乘积第二项是电场与检验函数相乘再乘上频率项。在COMSOL的“弱形式PDE”接口里我们可以直接把被积函数写成“Weak Expression”-(1/ε_r) * (Ex*test(Ex) Ey*test(Ey)) (ω²/c²) * Ez*test(Ez)这里的ε_r可以是任意表达式。比如Drude色散模型可以写成ε_r ε_inf - ωp²/(ω² i*γ*ω)其中ωp是等离子体频率γ是阻尼率在COMSOL里定义为全局变量即可。此时介电常数将随频率和波矢量变化迭代求解时会自动满足色散关系这正是弱形式方法相比内置求解器的核心优势。操作上我会把E_z定义为弱形式PDE的因变量在“弱表达式”中输入上面的被积函数然后添加全局方程来约束k点和本征频率之间的关系。这样在扫描k路径时能带结果能直接反映材料的频散特性。6.3 不收敛问题的调试经验弱形式求解最大的注意点是初值。弱形式方程本质是非线性的需要一个接近真实解的初始猜测。如果从零初值出发求解器几乎必然发散。我的做法是先用无频散模型算出同一结构的能带把某个特征频率作为色散模型的初始猜测再逐步增加色散项的强度。另一个问题是模式简并。当两个模式在同一个k点频率相同弱形式求解器可能同时找到两个解并发生混淆。我通常给结构加入一个很小的人工扰动比如把某个空气孔半径从180 nm改成180.1 nm专门破缺对称性。这样算出的能带会有轻微劈裂但模式是干净的。确认完模式属性后再把扰动归零重新计算。这个技巧在三维光子晶体中同样有效尤其是处理重频模式时特别管用。7. 高频踩坑点与排查速查表7.1 能带锯齿网格问题还是物理问题能带图上出现锯齿状的尖刺多数情况下是网格问题但偶尔也有物理来源。判断方法很简单把目标频段附近的网格尺寸减半重新计算同一区域。如果尖刺消失说明是网格欠分辨如果尖刺还在那可能是模式简并劈裂需要用扰动法或模式分解来确认物理性质。在二维案例中锯齿经常出现在布里渊区边界附近这是因为边界处场分布剧烈集中局部网格密度不足会带来伪频移。我专门在高对称点附近加“角点细化”因为三角晶格中这些区域场梯度最大。这个方法在三维案例中也一样有效能减少不少返工时间。7.2 模式缺失与特征值搜索范围特征值缺失是复现案例时非常头疼的问题。发现某个高对称点附近缺少一条本该存在的能带第一反应不该是改几何而是检查特征值搜索范围。COMSOL特征值研究设置中有两个关键参数一个是“所需模式数”决定求解器返回多少个特征模式另一个是“搜索基准点”决定搜索中心的频率位置。如果频带跨度大我会把搜索基准点设置到目标频带中间而搜索范围设置为整个感兴趣频段的1.5倍。同时把所需模式数提高到8以上跑完后手动过滤不需要的模式。还有一点务必注意COMSOL默认特征频率单位是rad/s如果输入以Hz为单位的值需要先做换算再填入。这个单位坑我踩过不止一次。7.3 归一化频率换算与单位陷阱光子晶体论文普遍用归一化频率a/λ而COMSOL的特征频率输出是freq单位rad/s。从freq换算到a/λ的公式很简单a/λ a * freq / (2πc)其中c是真空光速。如果你习惯用波长作为输入还要额外注意波长与频域求解频率之间的对应关系。这个换算看似简单但做40多个案例时反复手动换算极容易出错。我的解决方案是在“结果”节点中新建一个“全局变量探针”把换算公式直接定义成变量名作为一维绘图的横坐标。所有mph文件都采用这种输出方式打开模型就能直接看到归一化频率的能带图免去了每次换算的工作量。7.4 常见问题速查表错误现象可能原因解决方向能带曲线锯齿形抖动网格过疏或界面网格不均匀加密界面网格开启局部自适应细化禁带位置偏移明显材料折射率虚部未清零检查材料属性删除损耗虚部特征值缺失频率搜索范围太窄增大搜索区间并调整所需模式数模式重复或发生交叉高对称点简并或对称性太高加入人工微扰破缺对称性结果全部为0端口或边界条件错误检查是否使用Floquet周期边界而非通用周期边界三维模型内存溢出网格自由度过高关闭自适应细化简化几何或用辅助扫描弱形式方法不收敛初始猜测偏差过大先用无频散能带的计算结果做初值能带图整体变形k点坐标未做倒基矢变换重新设置Floquet边界波矢量表达式透射率曲线震荡端口宽度不等于一个完整周期将端口宽度调整为一个周期长度最后再分享一个小技巧mph文件保存前最好执行一次“文件-压缩模型”把求解过程中遗留的临时数据和无效网格节点清理掉文件体积能缩小两到三成。这不仅方便版本管理也方便和同行交换模型时减少传输压力。另外复现案例时建议把每个模型的物理单位、频段、材料参数记录下来写在README里否则三个月后回看你可能连自己的建模思路都忘了。这套案例库目前仍在扩充下一步我准备加入更多拓扑光子晶体和谷态输运的内容有新的进展会继续整理出来跟大家交流。