做光子晶体仿真这行,绕不开能带计算,也绕不开COMSOL。我用COMSOL 5.6把一本经典光子晶体教材里的案例逐个复现,目前手里攒了40多个mph文件,从一维多层膜到二维三角晶格再到三维反蛋白石结构都有覆盖。这篇就把整个复现过程的思路、操作细节、以及那些教材里不会写的坑一次性说清楚,给正准备动手算色散关系、画能带图的朋友做个参考。
很多人拿到教材第一反应是“照着例子搭模型就行”,真上手才发现,教材给的是思路和最终结果,中间的过程全靠自己补。尤其是光子晶体这类周期性结构,边界条件一个设错,算出来的能带就可能整体错位,更不用说三维结构动辄几十万自由度的求解代价。这个项目做完,我最大的收获不是“会点了软件按钮”,而是对整个电磁波在周期介质中传播的物理图景清晰了很多。
1. 项目缘起:为什么不停留在看书,而是动手复现
1.1 复现不是“照抄”,而是逆向理解理论
教材里的能带图都是经过验证的标准结果,我们复现的目的不是为了得到一张一模一样的图,而是通过搭建模型的过程,把教材里那些一笔带过的关键步骤补全。比如,一维光子晶体的布拉格反射条件、二维光子晶体带隙的出现条件、三维结构中的全方向带隙等,每一个结论背后都对应着特定的结构参数和物理假设。
以二维光子晶体为例,教材通常会给出三角形空气孔阵列在介电背景中的能带图,看起来似乎只是“求解一个特征值问题”。但真正动手时,你得处理倒格矢基矢的选择、不可约布里渊区路径的确定、TE与TM模式的区分、以及Floquet周期性边界条件的设置。这些细节如果不亲手操作一遍,理解永远是浮在表面的。
而且COMSOL 5.6这种通用有限元软件,不像专门的光子晶体工具那样把所有东西封装好了。你需要自己决定用什么物理场接口、用什么本征求解器、怎么设置波矢扫描参数,这个过程逼着你把底层原理搞清楚,这对后续做自己研究课题的能力提升非常有价值。
1.2 40多个mph文件怎么分类管理
项目进行到中期,mph文件数量上来之后,文件管理就成了一个实际的问题。一维、二维、三维结构在模型设置上的差异很大,如果全部堆放在一个文件夹里,后期连自己都分不清楚哪个文件对应哪个案例。
我按结构维度加物理场景的方式做了分类目录:
1D_bragg_mirror/:一维多层膜、布拉格反射镜、缺陷模腔等2D_square_lattice/:正方晶格介质柱、空气孔结构2D_triangular_lattice/:三角晶格空气孔、六角晶格介质柱3D_inverse_opal/:反蛋白石结构、蛋白石结构、woodpile结构test_cases/:临时调试的模型,跑通后再搬入正式目录
每个mph文件命名时带上结构参数和求解信息,比如hex_hole_r045_eps11_TE.mph这样的格式,一眼能看出是什么结构、占空比多大、介质常数多少、算的是哪个偏振。这个习惯帮我在后续对比不同参数结果时省了大量时间。
2. 一维到三维:光子晶体建模的差异与实操要点
2.1 一维结构:周期边界设置与能带折叠现象
一维光子晶体看起来简单,实际是检验Floquet周期边界理解是否到位的最佳练习。拿交替多层膜结构来说,每个周期里包含高低折射率两层介质,沿着一个方向做周期性排列。
模型搭建时,核心在于单位胞元的选择和周期边界的设置。在COMSOL 5.6的波动光学模块中,使用电磁波、频域接口时,需要在相对的两对边界上施加Floquet周期性条件。注意,一维情况下,只需要在传输方向的两个边界上设置周期性,垂直于传输方向的两个边界可以使用完美电导体或完美磁导体条件来模拟无限大平板在横向的延展。
这里有个容易出错的地方:周期矢量的方向设置。COMSOL中Floquet周期性的波矢k是依赖边界对之间的位移向量来定义的。如果位移向量设置错误,或者边界方向弄反,算出来的能带就会出现奇怪的“多带折叠”现象,甚至完全不对。经验是用periodic边界条件时,把源边界和目标边界选在同一坐标轴的两端,然后检查位移向量方向是否与你的晶格矢量一致。
能带折叠现象的复现也是一个有趣的过程:在均匀介质中,能带图沿着波矢是第一布里渊区内的一条直线;加入周期调制后,能带在布里渊区边界发生折叠并打开带隙。你可以通过扫描波矢从0到π/a(a为周期长度),观察前几条能带的变化,能直观看到带隙的形成过程。
2.2 二维结构:晶格类型、布里渊区路径与模式偏振
二维光子晶体是复现项目中数量最多的一类案例,因为结构类型多、物理现象丰富。正方晶格介质柱和三角晶格空气孔是最常见的两种构型,这两种结构虽然建模类似,但布里渊区形状完全不同,高对称点的位置也就不同。
二维建模中,单位胞元通常是矩形(正方晶格)或平行四边形/菱形(三角晶格)。实际操作中,我建议直接用矩形或平行四边形几何,再用周期边界条件把它们拼接起来。COMSOL中二维单位胞元的几何可以很简单,但周期边界的设置必须小心:三角晶格中,x方向的平移矢量和y方向的平移矢量并不垂直,需要在周期边界设置中正确给出这两个平移向量的分量。
波矢扫描路径方面,正方晶格的不可约布里渊区路径是Γ-X-M-Γ,三角晶格是Γ-M-K-Γ。扫描时用参数化扫描,把波矢k的x和y分量分别定义为扫描参数的函数。我用sin和cos组合来构造直线路径,比如Γ到X是kx从0到π/a、ky=0,X到M是kx=π/a、ky从0到π/a。每个路径段的插值点数量选择20~30个,太少曲线不够平滑,太多求解时间翻倍。
TE和TM模式的区分是二维光子晶体的一个关键点。在COMSOL中,三维结构块做二维近似时,TE模式对应电场垂直于周期性平面(实际上是沿着孔轴方向),TM模式对应磁场沿轴向。不同结构的带隙对偏振的依赖差异巨大:正方晶格介质柱结构通常在TM模式出现较宽的带隙,而三角晶格空气孔结构在TE模式下更容易出现带隙。复现教材案例时,同一个结构起码要分别算一次TE和TM,对比能带差异,这样才能理解为什么教材上对同一种结构会给出两张不同的能带图。
2.3 三维结构:计算代价与模型简化策略
三维光子晶体是复现难度最大的一类。全三维有限元模型的计算量非常可观,尤其当结构复杂(如反蛋白石、螺旋结构)时,一个完整的能带计算可能需要求解数目众多的特征值。
解决思路有两个方向。第一个方向是严格建模但控制求解规模。利用结构的对称性简化模型,比如woodpile结构可以在某些对称面上使用对称边界条件,把计算域缩小到原来的1/2甚至1/4。网格采用自适应策略,在介质交界面附近加密,在大块均匀区域用稀疏网格。我见过一个比较极端的做法:把空气孔和介质区域分开划分网格,介质区用二阶或三阶拉格朗日单元,空气区用一阶单元,在不牺牲精度的前提下大幅减少自由度。
第二个方向是降维近似。很多教材案例中的三维结构实际具有二维周期性(比如打孔板),可以先用二维模型计算面内色散关系,再通过有效折射率方法估计垂直方向的传播特性。这种近似对某些结构的定性分析足够,但无法捕捉真正的三维带隙。所以在复现过程中,我给自己定的原则是:教材明确给出三维能带图的案例才用全三维模型计算,其他场景能用二维近似就先做二维。
三维计算的另一个痛点是特征值求解器的选择。COMSOL 5.6的本征频率求解器,默认使用MUMPS或PARDISO直接求解器。当模型自由度超过几十万时,直接求解器的内存消耗非常夸张。我的经验是优先考虑使用"特征值求解器"中的IRAM(隐式重启阿诺尔迪方法)配合稀疏直接求解器,这一步设置在5.6中比较稳定。同时把“位移截断”适当调大,能有效减少迭代次数。更有经验的做法是先做一次粗网格计算,估个量级,再在细网格上只搜索目标频率附近的特征值。
3. 能带计算的核心:本征求解设置与弱形式方程处理色散材料
3.1 标准能带计算的求解框架
在COMSOL中计算光子晶体能带,最标准的做法是使用波动光学模块的电磁波、频域接口,给单位胞元设置Floquet周期边界,然后做“特征值”研究。
这里需要解释一下特征值研究到底在解什么。电磁波在周期性介质中的传播满足布洛赫定理,电场可以写成周期函数与平面波的乘积。代入波动方程后,得到一个关于本征频率ω和波矢k的本征值问题。给定一个k(即给Floquet边界条件中的k分量赋值),求解得到的特征值就是该k对应的模式频率。把k沿着布里渊区高对称路径扫描一遍,把所有频率画出来,就是能带图。
操作上,每给一个k做一次特征值求解,效率太低。更好的做法是使用参数化扫描定义k路径,每次扫描点求解一次特征值,然后把结果叠加显示。为了加速收敛,可以在研究设置里指定搜索感兴趣的频率范围,比如“附近搜索”选项,给出目标频率中心值和搜索半径。我的习惯是每次扫描前先观察前一步结果的频率分布,据此确定下一段路径的搜索范围,这样既能保证模式连续追踪,又能减少无用的特征值计算。
网格对特征值求解精度的影响需要特别说明。有限元离散后的本征值问题,误差主要来自网格对场分布的逼近。介质折射率突变界面处的网格必须足够细密,尤其在高折射率对比度的情况下(比如硅(n≈3.5)与空气),边界附近的场梯度很大,稀疏网格会造成频率偏高的系统性误差。我通常在高折射率区域设置至少5~6层网格节点,并且相邻网格单元的生长率不超过1.3。
3.2 用弱形式方程求解色散光子晶体能带
对大多数色散材料(比如金属、半导体等),介电常数是频率的函数。此时标准电磁波本征值问题中出现了一个非线性项,直接使用COMSOL内置的“特征值”研究会出现收敛困难,因为传统特征值求解器是为线性本征问题设计的。这种情况下就需要用到弱形式方程。
弱形式方程的思路是这样的:先从波动方程出发,乘以测试函数并分部积分,得到弱形式。以无源、非磁性介质中电场E为例,波动方程写为:
[ \nabla \times (\nabla \times E) = \left(\frac{\omega}{c}\right)^2 \varepsilon(\omega) E ]
在COMSOL的“PDE”接口中选择“弱形式PDE”模式,把电场E的每个分量作为因变量,然后输入对应的被积表达式。对x方向分量,弱形式表达式大致是:
-(curl_E_x*test_curl_E_x + curl_E_y*test_curl_E_x + curl_E_z*test_curl_E_x) + (omega/c)^2*epsi*(E_x*test_E_x)实际到界面操作时,COMSOL弱形式界面会用符号表达式,需要你定义辅助变量如omega和epsi。其中epsi可以用解析函数定义,比如Drude模型的无损形式:
epsi = 1 - omega_p^2 / (omega^2 + i*gamma*omega)这里的omega_p是等离子体频率,gamma是碰撞频率。由于omega同时出现在方程中作为待求本征值以及介电常数表达式里的变量,问题变成非线性的本征值问题。直接解依然困难,一个常用技巧是引入辅助变量把非线性本征问题线性化,或者使用“辅助扫描”方法,在实频域内扫描ω,在每个固定ω下求解本征值中的波矢k。
我复现等离子体光子晶体案例时用的就是这个思路:固定频率ω,把问题转换成“k本征值”问题。此时弱形式方程中,材料参数会保持为该固定频率下的常数,原问题变成关于k的线性本征问题,COMSOL的标准特征值求解器可以处理。这个方法对大范围色散谱线的计算非常有帮助。
弱形式界面的操作细节需要注意几点。第一,因变量要单独定义,每个分量都对应一个名称,在弱形式表达式中使用分量名称与测试函数名组合。第二,周期边界条件的引入方式和内置物理场稍有不同,需要在弱形式中加入周期贡献项,或者使用COMSOL的周期边界节点与弱形式耦合。第三,弱形式里单位制容易出错,varepsilon和真空波数这些量必须有明确定义,否则算出来的频率差出几个数量级却不自知。
4. 复现路上的典型问题与排查实录
4.1 常见问题速查表
长时间的复现过程,踩过的坑五花八门,这里整理一个速查表,很多问题不是看一眼模型就能发现的,每一条都是花时间排查换来的经验。
| 现象 | 可能原因 | 排查方向 |
|---|---|---|
| 能带沿波矢方向不连续,出现跳跃 | 波矢扫描路径定义错误或参数范围重叠 | 检查参数化扫描中各段的起点和终点,确保连续性 |
| 高频模式明显偏大(相对教材数据) | 网格分辨率不足,高折射率区未加密 | 加密高折射率区域网格,选用高阶单元 |
| TE模式算出了TM模式的场型 | 边界条件中面外方向设置错误 | 检查完美电导体/完美磁导体的分配是否与偏振定义一致 |
| 特征值求解速度极慢 | 搜索范围过大或请求模式数过多 | 缩小搜索频率范围,减少模式数量,分批求解 |
| 出现频率为0的伪模式 | 本征问题存在静态解,或自由度未约束 | 固定某一场分量(如设定E_z=0对TE模式),去除退化自由度 |
| 三维模型内存溢出 | 网格自由度过大,直接求解器内存消耗过高 | 使用迭代求解器如GMRES,或对模型做对称简化 |
| 三角晶格计算时K点简并度丢失 | 单位胞元几何对称性被网格破坏 | 开启网格角度细化,确保晶格高对称点在网格中保持对称 |
| 曲线出现负群速度段但教材没有 | 模式追踪错误,高、低支发生了交叉 | 单独画各条模式的场分布,观察模式演化,按物理连续性重新归类 |
4.2 实战排查:从“算不出”到“算得对”
举一个实战中记忆犹新的例子。我在复现一个二维正方晶格金属柱阵列的色散曲线时,用内置电磁波接口计算,结果在低频段出现了一大堆虚数频率特征值,完全无法归类为物理模式。
排查过程是这样的:先检查网格,发现金属柱边界附近网格虽然加密了,但介质区域与空气区域的界面交接处生成了两个独立的边界层,导致周期边界两侧的网格节点不对齐。Floquet周期条件要求网格在周期边界上逐节点对应,网格不对齐会让周期条件的数值实现产生误差,直接产生大量伪模式。重新生成一致化网格后,这个问题立刻消失。
随后又出现新情况:金属是色散介质,用Drude模型时,低频率下介电常数为负,本征值求解器在这个区间表现不稳定。改用前面提到的弱形式+辅助扫描方法后,问题迎刃而解,而且能算出色散曲线中由表面等离激元造成的平坦带区域,这是内置接口几乎无法给出的结果。
另一个常见的“让新手崩溃”的问题是:能带图整体看起来差不多,但某些高对称点(比如K点)的简并没有出现。这与单位胞元几何构建方式直接相关。三角晶格如果建模时用矩形裁剪出一个平行四边形区域,高对称边界处的网格无法完全对称,简并就被数值误差破坏。解决办法是直接建立正六边形单位胞元,网格生成器对这种高对称区域的对称性处理要好很多。
5. 一些个人经验和建议
5.1 从二维入手是最高效的学习路径
回看这40多个案例的复现过程,我的建议是:如果你刚开始接触光子晶体仿真,老老实实从二维结构入手,把正方晶格和三角晶格的能带算通,建立对周期边界、布里渊区路径、模式偏振这些核心概念的直觉,再碰三维结构。直接上三维案例很容易被求解时间和内存问题消耗掉耐心,反而不利于理解物理。
二维算通了,你会积累起几个几乎不需要动脑的“模板能力”:比如单位胞元几何的高效画法、参数化扫描波矢路径的快捷设置、不同偏振下特征值求解的配置方法。这些模板能力迁移到三维结构建模中,能帮你把精力集中在新问题本身,而不是重复踩同样的坑。
5.2 合理使用mph文件管理你的“知识资产”
40多个mph文件,本质上是我逐行调试、逐一验证过的知识库。每次新问题的建模,我几乎都能在自己已有的mph文件里找到最接近的结构模板,然后在它的基础上修改。这种积累方式的效率远超从零开始建模。
文件管理上有一个深刻教训:早期不注意保存求解过程设置,导致很多文件换台电脑打开后需要重新配置求解器参数。所以后来每个最终的mph文件我都保持“可立即重新计算”的状态:清空旧解、保留全部参数化扫描设置、在文件描述区记录结构参数和关键物理假设。这样换电脑、换版本、甚至分享给合作者,都能很快恢复工作环境。
5.3 给初学者的三条核心建议
第一,不要太依赖教材附带的“标准答案”,自己跑出来的能带图,哪怕和教材有微小的偏差,只要量级和趋势对,这个过程的收获就远超直接下载别人的模型。第二,养成记录日志的习惯,每次修改模型参数后把结果截图和参数记录在文本文件里,避免一周后完全忘记这个文件当初是怎么算出来的。第三,遇到算不出来的情况不要反复重试同一个求解器设置,应该回头思考问题的数学本质——是用错了接口,还是周期性条件表达不完整,还是材料属性定义有误?很多时候,思路转换比参数调整更有效。
这套思路同样适用于其他周期性结构仿真场景,比如声子晶体、超表面、微纳光学器件设计。工具是相通的,核心永远是清晰地理解你正在求解的物理问题。