做光子晶体仿真的同学基本都会在能带计算这道坎上卡一阵子。一二维还好一旦进入三维电场是矢量场三个分量互相耦合布里渊区变成高维结构常用的平面波展开法经常遇到收敛慢的问题而 Comsol 弱形式接口让我在方程层面完全打开局面——边界条件、材料张量、增益损耗都可以直接写进积分项里。这篇就完整拆解一次三维光子晶体能带计算的弱形式路线物理背景、方程推导、Comsol 实操步骤以及我从真实计算里踩出来的一堆坑。适合已经会基础有限元、但正被三维光子晶体本征值问题折腾的仿真工程师。1. 三维光子晶体能带计算的问题定位1.1 我们要算的到底是什么光子晶体本质上是一种介电常数周期性变化的结构。光在这种结构里传播时会像电子在晶格中一样形成能带结构给定一个波矢 k只存在若干离散的本征频率 ω(k)这些频率和波矢的关系就是能带图。能带图上出现的频率间隔就是光子禁带——这个频率范围内的光无法在晶体中传播。三维光子晶体比一维和二维难在一个关键点上电磁场是矢量场TE、TM 这些在二维里能分开讨论的偏振模式在三维原胞里完全耦合在一起。再加上三维结构往往需要较多的网格自由度计算规模上升很快。所以很多在二维里管用的简化方法到三维就得重新掂量。能带计算的横轴是波矢纵轴是频率通常用归一化频率 a/λ晶格常数除以波长来表示。一个完整的能带图需要沿着不可约布里渊区边界扫掠波矢在每个波矢点求出一组本征频率。这里面的工作量主要来自两点一是本征值求解本身二是波矢扫描的数量。三维情况下高对称点之间的路径更长简并度也更高对求解器提出了额外要求。1.2 为什么我不建议直接用平面波展开法平面波展开法是光子晶体能带计算的经典方法把介电常数和电场都展开成平面波的叠加将麦克斯韦方程转成矩阵本征值问题。这个方法在一维、二维结构里很直观但到了三维就有不少别扭的地方。首先是收敛问题。折射率差越大要准确表示电场的不连续变化需要的平面波数量就越多。三维情况下平面波数量增大一个量级矩阵规模迅速膨胀内存和耗时都让人头疼。其次是处理各向异性或色散介质比较麻烦介电函数一旦变成张量或者依赖频率平面波展开的公式推导就要重新来一遍。最后如果想算缺陷态、超胞等非理想周期结构平面波法虽然也能嵌套但实现复杂度会明显上升。有限元方法在几何适应性上有天然优势。复杂的三维原胞球体、孔洞、交叉结构可以用四面体网格准确离散材料属性可以定义为空间坐标的任意函数。更重要的是有限元求解的是积分形式的弱方程这给了使用者很大的自由度——我可以在方程里直接加入各向异性张量、表面电流项、非线性项而不用重新推导整套算法。1.3 从内置模块到弱形式我的选择思路Comsol 的波动光学模块自带 Floquet 周期边界做标准三维光子晶体能带计算非常方便图形界面上设置一下就能跑。那为什么还要费劲写弱形式我的真实体会是内置模块把很多细节封装好了但遇到以下情况就束手束脚。比如想在能带方程里加入磁光材料的介电张量偏置项想在原胞里加一层表面电导率边界想用另一种形式的特征值方程来研究增益介质或者想把能带计算和电路模型耦合。这些场景下内置接口要么不支持要么需要用户用“弱贡献”去叠加反而比直接建一个弱形式模型更绕。弱形式接口相当于把有限元方程的“公式源文件”暴露给我我可以自己写控制方程、自定边界约束、定义特征值变量。代价是设置过程变长且需要理解弱形式的数学原理。这篇文章要送的“鱼”是一套能跑的三维光子晶体弱形式模型“渔”则是掌握怎么把麦克斯韦方程变成 Comsol 能吃的弱形式。2. 从麦克斯韦方程到弱形式推导2.1 控制方程与归一化坐标先写出无源区域电磁波的控制方程。在频域下电场满足∇ × (1/μ_r ∇ × E) - k₀² ε_r E 0其中 k₀ ω/c 是自由空间波数μ_r 是相对磁导率ε_r 是相对介电常数。对于大多数光子晶体材料μ_r 1方程可以简化为∇ × (∇ × E) - k₀² ε_r E 0如果直接用真实坐标 r 和真实波数 k₀Comsol 里特征值 λ 的单位会是 1/m²后处理还要反复换算。我习惯把坐标归一化令 s r/aa 是晶格常数。这样所有几何尺寸都变成了无量纲量方程变成∇_s × (∇_s × E) - (a k₀)² ε_r E 0记 λ (a k₀)² (ω a/c)²本征方程就是∇_s × (∇_s × E) λ ε_r E这个 λ 是无量纲的后处理算归一化频率时直接取 sqrt(λ)/(2π)就是 a/λ_0。这个小技巧能省去很多单位换算的麻烦数据导出也干净。2.2 加权余量与分部积分弱形式的核心思想是把偏微分方程从“在每个点上精确成立”放宽成“在一组权重函数下加权平均成立”。如果场 E 满足原来的强方程那么对于任意权重函数 F下面的积分一定为零∫Ω F · [∇_s × (∇_s × E) - λ ε_r E] dΩ 0反过来如果这个式子对很多“足够多”的 F 都成立那么 E 就是方程的解。有限元方法就是用有限个局部定义的 F 来逼近这个“无限个条件”。直接对这个积分做处理会遇到一个问题方程里有两阶导数对 E 的光滑性要求很高。好在有矢量恒等式F · (∇ × A) (∇ × F) · A - ∇ · (F × A)把它代入第一项并在整个原胞体积上积分利用散度定理就得到∫Ω [(∇ × F) · (∇ × E) - λ ε_r F · E] dΩ - ∮∂Ω F · [n × (∇ × E)] dS 0第一次看到这个推导的人可能会觉得只是数学变形但实际上这一步很重要导数阶数从二阶降到了一阶边界上出现了通量项 n × (∇ × E)。有限元网格只需要保证场本身连续不必要求导数一阶连续网格适应性大大增加。至于边界积分项在周期边界条件下相对面上的贡献会成对相消或者由Bloch相位关系统一处理。2.3 Bloch定理与弱形式里的周期条件能带计算必须体现晶格周期性。Bloch 定理说周期结构中的本征场可以写成E_k(r) u_k(r) e^{i k·r}其中 u_k 是与晶格周期相同的函数。这个表达式隐含了一个重要关系在原胞的相对面上场分量的取值满足E(r a) E(r) e^{i k·a}这里 a 是晶格平移矢量。在弱形式中这个条件不会自动满足必须作为约束条件显式加入。这也是整个模型里最容易出错的部分。Comsol 里处理周期约束有两种路径。一种是用内置的“周期条件”节点把一对相对面设为源边界和目标边界再附加 Floquet 相移因子。另一种是在弱形式方程里手动写约束表达式。无论哪种都要注意三个分量 Ex、Ey、Ez 各自独立满足相位关系方向不能弄错。还有一个容易忽略的细节Bloch 相位是复数因此电场因变量必须设置为复数因变量。如果忘了勾选复数选项求解器会把相位因子拆成实部和虚部处理得到的结果往往混乱。3. Comsol 弱形式模型实操步骤3.1 建立三维原胞几何以简单立方晶格为例原胞是边长 a 的正方体。为了配合归一化坐标我直接把几何建成无量纲形式正方体边长 1中心有一个半径为 0.25 的介质球。背景为空气ε_r 1球体用类硅材料ε_r 12.25。实际操作时Comsol 的几何单位可以保持 m在全局参数里设 a 1 μm球半径 r 0.25 μm最后再在物理场设置里用坐标除以 a 做归一化也是可以的。不过我用无量纲几何更顺手因为后面写弱表达式和扫掠波矢的时候所有量都一致。网格方面自由四面体就够了。球体附近需要加密背景区域可以相对稀疏。我通常把最大单元尺寸设为 0.05曲面网格细化至 0.02这样能在保证精度的情况下兼顾内存。第一次试跑建议更粗一些等确认结果合理了再加密。3.2 因变量与弱表达式的设定在组件定义里添加物理场数学 → PDE 接口 → 弱形式 PDE。因变量类型设为矢量场 E分量 Ex、Ey、Ez。由于 Bloch 相位的原因在因变量设置中必须勾选“复数因变量”。弱形式方程是核心。对应前面推导的无量纲形式弱表达式为-(1/mu_r)*dot(curl(E), curl(test(E))) lambda*eps_r*dot(E, test(E))拆开看这个表达式curl(E)计算电场旋度test(E)是权重函数 F 的表示dot(a,b)是矢量内积lambda是 Comsol 弱形式接口自动识别的特征值变量对应上面的 (ωa/c)²eps_r在这里写成分段表达式例如if(sqrt(x^2y^2z^2)r, 12.25, 1)mu_r通常设为 1。由于方程是对矢量场写的Comsol 会把这个弱表达式自动展开成针对 Ex、Ey、Ez 的三个方程。测试函数的分量分别对应三个方向的权重函数这正是矢量有限元处理耦合的方式。写这一步时我建议每次只加一项先只保存旋度项确保无源模求解正常再加上特征值项便于排查错误。3.3 周期边界与 Floquet 相移的具体实现在弱形式PDE节点的子菜单中添加三个“周期条件”节点分别处理 x 方向、y 方向和 z 方向的相对面。每对边界中一个设为源边界另一个设为目标边界。普通周期条件是目标边界上的场等于源边界上的场。Floquet 相移则要额外乘一个相位因子对于 x 方向的一对面约束表达式是Ex_dest - Ex_src*exp(i*kx*a) 0 Ey_dest - Ey_src*exp(i*kx*a) 0 Ez_dest - Ez_src*exp(i*kx*a) 0y 和 z 方向同理相位因子分别用 ky 和 kz。这里的 kx、ky、kz 是全局参数后续通过参数扫描改变。有一个实操上的坑周期边界的源/目标面选择必须和网格一一对应。如果两个面网格不一致Comsol 会自动做插值但精度会下降还可能出现伪模式。我在实际项目中遇到过一次莫名其妙的杂散模式最后发现是相对面网格不一致导致的。解决办法是在“周期条件”设置里强制使用相同的边界网格或者手动设置源和目标面的网格分布。3.4 特征值求解与波矢扫描研究节点选择“特征值”特征值变量就是 lambda。由于 lambda 是无量纲的 (ωa/c)²搜索范围要按预期频率设置。比如想看归一化频率在 0.3 到 0.8 之间的能带lambda 就在 (2π×0.3)² ≈ 3.55 到 (2π×0.8)² ≈ 25.3 之间。搜索范围稍放宽到 1 到 30特征值数量设为 12 到 15 个这样可以避免漏掉简并模。波矢扫描沿不可约布里渊区边界进行。对于简单立方晶格典型路径是Γ → X → M → R → Γ对应的 k 点坐标分别为Γ: (0, 0, 0)X: (π, 0, 0)M: (π, π, 0)R: (π, π, π)注意这里因为几何归一化坐标乘的是 π 而不是 π/a。如果用真实坐标乘的是 π/a。把路径用一个参数 t 进行线性插值在 Comsol 的“参数化扫描”里对 kx、ky、kz 同时扫掠即可。扫描步长一般取 20 到 40 个中间点能带轮廓就能看清了。参数扫描时有一个提升效率的细节不要每个 k 点都从零开始求解而是在研究设置里开启“在上一解附近搜索特征值”或使用前一个解的初始猜测。这样可以减少迭代次数尤其在高对称点附近能明显提速。3.5 能带后处理与绘图计算完成后每个 k 点会得到一组特征值。为了绘制能带图我通常先把数据导出为文本然后用脚本画图而不是完全依赖 Comsol 自带的绘图工具。原因是特征值本身不带“能带编号”含义简并点附近排序可能跳跃脚本里可以按频率排序并作平滑处理。导出格式可以是路径参数 t, kx, ky, kz, 特征值1, 特征值2, ...脚本里将特征值转换为归一化频率f_norm sqrt(lambda) / (2*pi)然后以 t 为横轴以 f_norm 为纵轴绘制多条曲线。每条能带颜色统一断开处通常对应简并点或高对称点附近的模式切换需要人工判断。4. 常见问题与排查技巧4.1 均匀介质验证法整个模型搭建完成后第一步不是急着算光子晶体而是把原胞里的介电常数全部设为 1跑一次均匀介质下的能带。这个情形有解析解色散关系就是 ω c k用归一化单位就是 f_norm k_norm/(2π)。如果弱形式、边界条件和求解器都正确算出来的第一能带应该和这条直线高度吻合。我每次新建模型都会做这个验证。它能在 10 分钟内排除掉 80% 的模型设置问题。如果均匀介质检不过后面光子晶体的结果基本不能信。4.2 周期边界造成的非物理模式周期边界设错的表现是能带图上出现大量低频杂散模式或者模式频率分布没有规律。最常见的原因有三个源边界和目标边界选反导致相位方向反了三个分量的相位约束没有全部加上只约束了法向分量网格不一致造成插值误差。检查方法很直接先看均匀介质验证能不能过再在某个 k 点单独求解用后处理画电场模的分布确认场在边界处相位匹配。4.3 特征值漏解与简并问题三维光子晶体在高对称点处简并度很高比如 Γ 点附近可能有多条模式频率相同。如果特征值搜索范围设置得太窄或者特征值数量设置得太少某些简并分支就会漏掉能带看起来断断续续。解决办法是适当增加特征值数量并把搜索范围扩大到目标频率上限的 1.5 倍。如果计算资源允许可以把每个 k 点的特征值数量设到 30 个甚至更多。漏解还有一种表现是能带在低对称点断开但高对称点连续这通常是网格不够细、简并被数值破缺导致的加密网格可以改善。4.4 复特征值和有损材料很多三维光子晶体会用到金属或掺杂半导体介电常数是复数或者负值。此时特征值可能变为复值而弱形式默认求解器可能只搜索实特征值。需要在特征值设置中允许复数特征值并将搜索范围改为复平面上的一个区域。另外如果介电常数本身依赖频率比如 Drude 金属模型那么特征值方程变成非线性特征值问题λ 不能简单地作为系数放在 ε_r 的表达式里。这种情况下我的经验是不要硬扛弱形式而是改用 Comsol 内置的电磁波频域接口配合“特征频率”研究。弱形式虽自由但不能违背求解器的基本约束。4.5 内存不足与速度优化三维原胞的网格量很容易突破百万自由度特征值求解又特别吃内存。几个我实测有效的优化手段先用粗网格跑通流程再用细网格精算扫描路径分段提交不要一次把 40 个 k 点全堆在一起求解器选择“SPOOLES”或“MUMPS”并开启多线程特征值搜索范围分段扫描避免求解器在空频段浪费时间。5. 扩展从能带走向光子器件设计能带算出来只是第一步。算准了能带接下来可以做很多事。比如加一个缺陷在原胞中去掉一个球再用 3×3×3 或 4×4×4 的超胞去算缺陷态能级。弱形式模型的好处此时就体现出来了不需要改任何方程只要把几何改成超胞并调整周期边界即可。再比如磁光光子晶体把介电常数设成反对称张量形式在弱表达式里加上对应项能带分裂和拓扑边界态的计算就能直接跑起来。另外还有一个很实用的小技巧在弱表达式里加入一个小参数来微调某一项比如一个微弱的增益项通过扫描参数观察特征值的虚部从负变正就能定量找到阈值增益。这类参数化分析在标准内置接口里做起来很绕但在弱形式框架里就是改一行表达式的事。我觉得弱形式的真正价值不是替代内置接口而是在你需要“理解自己在解什么方程”的时候给你一条不被封装逻辑困住的路。用熟之后你会觉得光子晶体能带计算不再是个黑箱而是自己亲手搭起来的一个物理模型。
