企业官网建设流程全解析

光子晶体仿真看起来门槛高，实际上绝大多数时间都花在修正模型的细节上。我最早跟着论文复现二维空气孔光子晶体，整整一周都在跟能带图里的锯齿较劲，最后才发现只是材料介电常数虚部没有清零。为了彻底摆脱反复试错的局面，我决定用COMSOL 5.6把《光子晶体》教材中的典型案例完整复现一遍。目前这套复现项目积累了40多个可直接运行的mph文件，涵盖一维、二维、三维三种维度的光子晶体结构，包括透射谱、反射谱、能带图、本征模场分布等完整输出。这篇文章既是对这套案例库的说明，也是把复现过程中踩过的坑、总结的方法和调参经验一次性讲清楚。无论你是刚接触光子晶体仿真的研究生，还是已经在算能带但经常对不上文献结果的工程师，这套从一维到三维的完整链路都值得参考。

1. 为什么我要在COMSOL 5.6里复现光子晶体案例

1.1 从“看书懂”到“动手会”的鸿沟

光子晶体相关书籍通常会把能带理论讲得很细，但真正动手建模的时候，你会发现书上的信息根本不够用。比如，书里可能只写了晶格常数a=600 nm、空气孔半径r=180 nm，却没有写清楚Floquet边界条件的波矢量该怎么设置，也没有说明用TM模式还是TE模式计算。这些信息缺失，导致一百个人能跑出一百种结果。

所以，我决定换一个思路：不再零散地在网上找示例，而是选一本案例最完整、参数最清晰的光子晶体专著，把其中能复现的算例逐个做出来。所谓复现，不是简单画出几何，而是让计算得到的禁带位置、能带宽度、透射率曲线和书中的结果一致。有了这样一条校准线，后续做新结构的时候，我才有底气去改参数、换材料，知道哪些结果是合理的，哪些是模型出错了。

1.2 案例库的组织方式与文件规范

40多个mph文件如果随意堆在一起，半年后连自己都找不到。我在整理案例库时，采用“维度—结构类型—求解目标”的三级目录结构：根目录下分1D、2D、3D三大类，每一类再按照具体结构细分。比如2D目录下，就有正方晶格空气孔、三角晶格空气孔、六角晶格介质柱、线缺陷波导、微腔谐振器等子目录。

文件命名也有一套固定的规则。我会把结构类型、关键几何参数、折射率组合、计算模式四个要素写进文件名。例如2D_Triangular_r0.3_n2.34_TM_band.mph，看到名字就知道这是三角晶格空气孔结构，r/a=0.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点路径之间的关系。几何尺寸可以用实际晶格常数建模，比如三角晶格a=600 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/a，kFloq2 = 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.46，TiO2设为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°。这里有一个常见错误：很多人直接在正方形外框里放一个圆形空气孔来模拟三角晶格，这等于改变了晶格对称性，算出的能带并不属于真正的三角晶格。

我的标准做法是：建一个平行四边形晶胞，使用全局参数定义顶点坐标，比如a=600 nm、r=180 nm，然后利用三角函数关系算出平行四边形的斜边顶点。Floquet边界恰好映射两个基矢方向，这样才保证计算模型的对称性正确。

4.2 k路径扫描的参数化实现

二维案例能带计算中，k路径扫描是最容易出错也最耗时的环节。我在全局参数里定义了三段扫描变量s1、s2、s3，分别对应Γ-M、M-K、K-Γ三段路径。每一段路径用线性插值把扫描参数映射到k点坐标。

比如Γ-M段，s从0到1，kx从0映射到0.5，ky从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 常见问题速查表

最后再分享一个小技巧：mph文件保存前，最好执行一次“文件-压缩模型”，把求解过程中遗留的临时数据和无效网格节点清理掉，文件体积能缩小两到三成。这不仅方便版本管理，也方便和同行交换模型时减少传输压力。另外，复现案例时，建议把每个模型的物理单位、频段、材料参数记录下来，写在README里，否则三个月后回看，你可能连自己的建模思路都忘了。这套案例库目前仍在扩充，下一步我准备加入更多拓扑光子晶体和谷态输运的内容，有新的进展会继续整理出来跟大家交流。