写这篇文章之前先交代一下来龙去脉。我之前接手一个项目要做三维反蛋白石结构的光子晶体能带计算手上正好有 COMSOL 的许可证但只有数学系列模块射频模块的授权还没批下来。后来折腾了几天发现 COMSOL 里真正扛事的其实是一套“弱形式”有限元内核只要你把这个内核用明白很多看起来只能靠专用物理接口才能做的事自己也能搭出来。这篇就把我踩过坑之后的完整思路写下来从麦克斯韦方程怎么一步步化成 COMSOL 能吃的弱形式再到三维周期性结构的边界条件、波矢扫描、伪模处理以及能带图后处理希望对正卡在光子晶体能带计算上的朋友有点帮助。内容偏硬核但我会尽量把每一步的“为什么”也讲清楚毕竟能带计算这种活光会点按钮是不够的。1. 为什么要用弱形式算三维光子晶体能带1.1 三维光子晶体能带问题到底在算什么光子晶体本质上是介电常数在空间周期变化的材料光在里面传播时会形成类似电子在晶体中那样的能带结构。能带图横轴是布洛赫波矢纵轴是归一化频率如果某个频率区间在所有波矢方向上都不存在传播模式就形成了光子带隙这就是所谓“光子绝缘体”的基础。典型结构包括反蛋白石、木堆结构、以及简单立方或面心立方排列的介质球。三维光子晶体和二维相比麻烦不在物理概念而在计算量。二维问题的求解域是一块平面变量分解成 TE/TM 两个标量方程就够了。三维结构里电场和磁场都是三维矢量每个点上有三个场分量而且它们之间通过旋度方程强耦合离散后自由度数会大一个量级。COMSOL 求解这类问题的方式基本上就是有限元它在求解过程中把原来的强形式方程转化为弱形式然后在有限维函数空间里找近似解。能不能把弱形式用明白直接决定你能不能自定义方程、改边界条件、加特殊材料模型。1.2 主流的三种算法为什么最终选有限元做光子晶体能带最经典的方法其实是平面波展开法把介电常数和电磁场都展开成傅里叶级数然后解一个矩阵特征值问题。PWE 的最大优势是代码写起来直观收敛性好尤其对光滑周期结构算能带非常快。但遇到复杂几何、多材料嵌套、带衬底的薄膜光子晶体PWE 的收敛速度就会明显下降而且它处理缺陷态、非线性材料和金属色散都比较吃力。FDTD 时域有限差分在器件模拟上用得也很多它可以直接通过傅里叶变换提取能带信息也能看透射反射谱。但 FDTD 要得到三维结构全布里渊区的能带往往要做很多次激励消耗的时间非常可观而且网格要满足数值稳定性条件格点多了之后让人很崩溃。COMSOL 有限元方案的优点在于直接在频域求解特征值问题网格可以贴合真实曲面处理复杂几何和局部细化完全没问题尤其是后来我需要在同一个框架里加入磁光耦合项和各向异性介电张量这时候通用弱形式几乎是绕不开的。COMSOL 内置的 RF 模块本质也走的是有限元框架它内部的方程视图里就把旋度-旋度算子和测试函数展开了完全符合弱形式逻辑。所以无论你用哪个模块理解弱形式都是理解 COMSOL 计算内核的第一步。1.3 弱形式思路把“求两次导”变成“求一次导”光子晶体本征方程里最麻烦的算子是旋度再旋度表达式写出来是二阶偏微分但有限元直接处理二阶算子很烫手因为要求基函数至少一阶连续这在任意网格上很难保证。弱形式的核心操作是拿一个测试函数去乘原方程再做分部积分把二阶导数导数转移到测试函数上。这样一来方程里只需要场的一阶导数基函数只需要容易实现的低连续性边界条件也天然被吸收到积分表达式中。说得再直白一点强形式像是在要求一根曲线“处处光滑可导”弱形式则允许在有限个点上出现导数不连续只要整体积分意义上满足方程就行。实际物理问题里电磁场的切向分量在介电界面处本来就可以不连续因此弱形式比强形式更贴近真实的电磁问题。2. 从麦克斯韦方程到弱形式特征值问题2.1 磁场的旋度-旋度方程电介质里不考虑自由电荷和电流取时谐因子 exp(-iωt)麦克斯韦方程组可以直接消去电场得到只有磁场 H 的本征方程∇ × (1/ε_r ∇ × H) (ω/c)^2 H这里 ε_r(r) 是随空间周期变化的相对介电常数c 是真空光速ω 是角频率。等式左边算子的分母里带着材料参数 ε_r所以磁场形式比电场形式在处理多材料结构时更方便材料界面上也不容易出额外数值问题。这就是本征值问题中(ω/c)^2 的作用COMSOL 的特征值求解器求的就是这个数。如果一定要用电场 E 来做方程形式是 ∇ × ∇ × E (ω/c)^2 ε_r E也可以但需要小心界面上电场法向分量的非连续性问题。我个人更喜欢用磁场形式因为在光子晶体这种介电常数频繁跳变的结构里H 的界面连续性条件更容易被有限元离散自然满足。2.2 Bloch 定理和波矢替换无限大的周期结构电磁场满足 Bloch 定理场量可以写成周期函数与平面波的乘积H(r) h(r) exp(i k·r)其中 h(r) 与晶格周期完全一致k 是第一布里渊区内的波矢。把上式代入旋度方程后对 r 的导数会同时作用于周期函数 h 和指数项出现交叉项形式上相当于做了一个替换∇ → ∇ i k实际操作中你可以把 H 直接代入方程再用 COMSOL 的 Floquet 周期条件处理 k也可以先定义新的周期变量 h对 k 做参数化扫描。后一种做法在自写弱形式时更顺因为 h 本身就是周期函数边界只需要最简单的周期相等条件所有与 k 有关的项都体现在方程表达式里。在 COMSOL 物理接口里你并不需要手动做这种替换在边界上指定 Floquet 周期条件并输入布洛赫波矢 kx、ky、kz 即可。理解替换过程的意义在于你将来如果要自己改方程、加非标准物理项就能清楚知道哪里会出现 i k 交叉项哪里会出现 -i k 交叉项。2.3 弱形式的完整推导把磁场形式的方程两边点乘一个测试函数 T(r)然后在原胞体积内积分∫ [ (1/ε_r)(∇ × H)·(∇ × T) - (ω/c)^2 H·T ] dV 0这里我跳过了一步分部积分时边界项在周期边界条件下会互相抵消所以最终体积分里没有额外边界项。如果是开边界或有金属边界这里就要补上相应的表面积分项COMSOL 的边界弱贡献就是干这个用的。推完之后就会发现方程里原来对 H 求旋度再被测试函数旋度作用实际微分阶数只有一次这就是弱形式名字的由来。有限元把所有连续问题离散成矩阵方程后就得到广义特征值问题 Ax λBx其中 λ(ω/c)^2A 是刚度矩阵B 是质量矩阵。2.4 在 COMSOL 方程视图里亲眼看弱形式很多人不知道 COMSOL 有一个非常实用的功能叫方程视图在物理场设置界面下方或者右键菜单中可以打开。打开之后你能看到每个物理场节点背后对应的弱贡献表达式比如电磁波频域接口会显示包含 μ0、εr、test(E) 等项的长表达式。第一次看到的时候确实又惊又喜原来 COMSOL 把 RF 模块的所有内部弱形式都摆在明面上。这意味着即便你没有编写自定义弱形式的能力也可以通过修改方程视图中某些幕后表达式来引入自己的模型项。比如我在一个磁光材料模型里需要在方程中加入非对称介电张量直接把旋度-旋度算子所对应的介电张量补充项写进弱贡献里比重新搭一套物理场接口要省事得多。所以这篇讲弱形式并不是只写给拿数学 PDE 接口硬干的人凡是想修改内置物理场方程的人都值得看。3. 实操用 COMSOL 算一个三维光子晶体能带3.1 几何建模和材料参数设置我用的是最简单的三维光子晶体模型边长为 a 的立方体单元中心放一个半径 r 的介质球背景为空气构成简单立方点阵。如果做反蛋白石结构可以把球换成空气孔背景换成高折射率材料。晶格常数 a 我取 1 μm介质球半径 r 取 0.3a球的相对介电常数 ε_r 取 12.25背景取 1这套参数默认工作在近红外波段附近且能在高频区看到较大的带隙。几何构建很直接在 3D 组件里建一个边长 a 的立方体块再在中心建一个半径 r 的球体用布尔差集或者并集做域划分。需要注意周期性光子晶体单元不需要把空气域留太多余量几何体严格等于一个原胞外面不加任何边界层。网格方面我自己习惯用自由四面体网格球面附近加密最小单元尺寸压到 0.05a 以下最大单元尺寸不超过 0.15a。三维矢量问题自由度很容易过百万这个网格尺寸在能带计算精度和计算量之间算比较平衡。3.2 物理场与 Floquet 周期条件设置物理场接口我建议先用内置的电磁波频域接口模块选择“电磁波频域”研究步骤选特征频率。这个接口的优点是自带棱边元离散比通用 PDE 接口手写矢量方程更稳后面我会讲为什么三维矢量场不建议直接用普通拉格朗日单元手写。周期性边界条件是重中之重。选中立方体所有外边界添加周期性条件类型选择 Floquet 周期。六个面自动配对成三组每组指定波矢偏移。COMSOL 里 Floquet 周期条件会建立源边界到目标边界的相位耦合背后的数学基础就是 Bloch 定理。你只需要把波矢分量设成全局参数 kx、ky、kz这样后面扫描时可以直接改变参数不需要重复修改边界条件。在域设置里别忘了把材料相对介电常数域选择正确介质球域填 12.25背景空气域填 1。磁相对磁导率都设为 1电导率保持 0。如果算的是金属光子晶体还需要加色散模型这里不展开。3.3 布里渊区高对称路径的参数扫描三维简单立方晶格的倒空间也是简单立方第一布里渊区高对称点分别是 Γ(0,0,0)、X(π/a,0,0)、M(π/a,π/a,0)、R(π/a,π/a,π/a)。能带图通常沿 Γ-X-M-Γ-R-M 这样的路径画。由于 COMSOL 的特征值求解器一次只能对一组建模参数求解你需要把这条路径转成单参数 s 的分段线性函数。我习惯把 s 定义成 0 到 4 的连续变量每一段对应一个高对称点之间的连线。为了省事情可以在 COMSOL 全局参数里直接写分段表达式。比如让 s 在 0 到 1 之间表示从 Γ 到 X此时 kx 0.5 * s * (2π/a)ky 0kz 0。然后 s 在 1 到 2 之间表示从 X 到 M此时 kx 0.5 * (2π/a)ky 0.5 * (s-1) * (2π/a)kz 0。后面 Γ 到 R、R 到 M 也类似。为了简化代码我常先定义无量纲参数 bkx、bky、bkz实际波矢分量再乘以 2π/a。在 COMSOL 参数表达式里可以直接写bkx if(s1, s*0.5, if(s2, 0.5, if(s3, 0.5-(s-2)*0.5, 0))) bky if(s1, 0, if(s2, (s-1)*0.5, if(s3, (s-2)*0.5, 0.5))) bkz if(s1, 0, if(s2, 0, if(s3, (s-2)*0.5, 0.5-(s-3)*0.5)))这个写法只给了参考框架实际使用时每个分支都要根据你的扫描点自己做分段。要注意的是 if 表达式在参数扫描里没问题但别在求解器内部反复切换否则容易造成不连续影响求解器收敛。COMSOL 的参数化扫描会把 s 按你设置的列表自动展开逐个计算特征值。3.4 特征值求解和后处理画出能带图研究选择特征值求解器设置里指定要求解的特征值个数我通常设置每个波矢点找最低的 8 到 12 个特征值。特征值搜索基准可以设置为 0也可以设置为某个目标值比如接近带隙中心频率平方的位置这样可以减少算到不感兴趣的高频模式。线性求解器用 MUMPS因为三维网格自由度很大MUMPS 在内存允许的范围内表现很稳定。求解完成后特征值列表给出的是 λ (ω/c)^2实际物理频率换算成f c √λ / (2π)能带图的纵轴常写成归一化频率 a/λ也就是 a √λ / (2π)。这里 λ 不是波长的符号容易混淆我一般直接用 f a / c 作为纵轴。导出数据的时候可以在全局计算里定义表达式 sqrt(real(λ))/(2*pi)乘上晶格常数 a 以后再画图。用二维截线图功能处理扫描结果很高效把 s 作为横坐标所有特征频率纵坐标画在同一张图上就得到能带结构。不要忘了把虚部接近非零的模式标记出来这些通常是伪模画图时可以筛掉。后面第 5 节我会讲怎么高效过滤伪模这是整个流程里最考经验的地方。4. 如果只用通用 PDE 接口手写弱形式怎么走4.1 先拿二维 TM 模型练手虽然三维矢量问题我建议用 RF 模块但如果你想真正吃透弱形式最好先用二维 TM 标量模型练手。二维光子晶体里TM 模式只有 Ez 分量电场垂直于周期平面波动方程退化成标量亥姆霍兹方程∇ · (1/ε_r ∇ E_z) (ω/c)^2 E_z 0这个方程乘测试函数再做分部积分弱形式表达式非常干净∫ [ -(1/ε_r) ∇E_z · ∇T (ω/c)^2 E_z T ] dV 0在 COMSOL 通用 PDE 接口里写这个表达式几乎就是一行的量。用这个例子你就能体会到弱形式表达式、测试函数、质量项系数之间的关系然后再去理解三维矢量问题里旋度和棱边元的概念就容易多了。从标量版本出发你也可以加入各种自定义现象比如增益介质、非线性折射率、磁光耦合等只需要在原式子上添加额外弱项不需要整个重写物理场。4.2 三维矢量场为什么不能简单照搬标量写法很多人第一次尝试用 COMSOL 的弱形式 PDE 接口做三维电磁波时会自然地把 H 拆成 Hx、Hy、Hz 三个标量因变量每个因变量写一个方程三个方程拼起来表示矢量方程。这个思路看起来顺理成章实际做起来却会面临一个很隐蔽的问题旋度-旋度算子的离散需要合适的功能空间。有限元里用 Lagrange 节点基函数表示三个直角分量并不能很好地逼近 H(curl) 空间中的场。节点基函数强制要求场分量连续但电磁场在材料界面处切向分量连续、法向分量可以跳变这种需求需要棱边元来满足。所谓棱边元就是把自由度定义在网格棱边上而不是节点上本质上是为旋度算子专门设计的一种基函数序列也就是常说的 Nédélec 元。COMSOL 的电磁波频域接口默认就使用了这类单元而通用 PDE 接口通常默认是 Lagrange 单元需要折腾离散化设置才能换成矢量元一般用户很难在 UI 层直接改。这不是说弱形式思路错了而是说实现时选择正确的离散空间和它同样重要。用件对应一个数学空间楼梯对应另一个数学空间你用错数值工具方程再漂亮也会算出一堆假模式。4.3 手写三维弱表达式的体验和补救方案如果你确实没有 RF 模块授权又必须手写三维矢量本征问题我建议参考通用 PDE 接口搭配矢量因变量的方式。COMSOL 6.x 的某些“数学 PDE”接口比旧版本更灵活可以在方程视图中手动指定更多矢量操作符。具体操作上把因变量组设为长度 3 的向量场然后在弱表达式里用 comp1.dx(Hy)-comp1.dy(Hx) 这类分量表达式构造旋度算子。因为完整的三分量弱表达式非常长很容易写错一个正负号我在实际项目里基本不会把这一步作为首选而是用内置 RF 模块 方程视图修改。真到必须自写时建议先在二维模型上逐项验证再扩展到三维。同理手写求解时也需要在弱形式里加入一个散度惩罚项强制 ∇·H 尽量接近零。这个惩罚项在方程里写成s ∫ (∇·H)(∇·T) dVs 取一个与刚度项同量级的小系数不能太大不然会压低物理模式。加了惩罚项拉格朗日单元可以抑制大部分伪模但并不能百分之百消除算完一定要检查场分布。我这边验证下来伪模的特征频率往往落在带隙里而真实的布洛赫模式会严格贴着能带曲线走。两者结合判断基本不会把伪模当成物理结果。5. 典型坑与排查笔记5.1 特征值全是虚频或者出现大量接近零的模式出现虚频和极低频率模式最常见原因是周期性边界条件没有完全配好导致本征问题缺少约束出现零能量刚体模式宽度和数量与自由度有关。检查步骤从三方面走第一确认三对相对边界都设置了 Floquet 周期条件且波矢参数表达式正常第二确认扫描参数不是零向量Γ 点处会有比较多的简并和零频模式这是物理现象但取带结构时容易被干扰第三如果有金属边界或默认绝缘边界混进来也要逐面检查。另一种情况是用通用 PDE 接口手写矢量场时由于没有散度惩罚旋度-旋度算子的零空间会被放大出现大量非物理模式。处理方法就是前面提到的散度惩罚项或者在求解器中检查“约束隐藏自由度”相关选项。5.2 能带曲线在布里渊区边界不闭合或者关键点处频率跳变能带计算中高对称点处有些模式理论上应该简并但如果几何网格或离散方式破坏了对称性简并会裂开。这个现象在 Γ 点和 M 点尤其明显。想要缓解就必须保证网格在高对称操作下也保持对称比如球体默认的网格生成器通常能保持旋转对称但如果你用了布尔操作后没有重新划分对称性可能被破坏。另外参数化扫描各分段边界处如果表达式不连续曲线就算断开所以分段表达式一定要在分段点严格连续。还有一种常见问题是波矢单位搞错。COMSOL Floquet 周期条件的波矢单位是 rad/m如果你直接拿简单立方布里渊区的高对称点坐标 π/a 乘进去实际波矢就是 π/a rad/m这没问题但如果你做的是无量纲参数扫描没乘上 2π/a得到的能带会在横向上被压缩或者错位。我的建议是在全局参数里单独定义 a、kx_scan、ky_scan然后再用 a*... 把它们组装成真实波矢。5.3 周期边界网格不匹配导致报错或收敛差COMSOL 的 Floquet 周期条件要求源边界和目标边界的网格节点位置一一对应否则会出现映射错误或者误差很大。自由四面体网格在三组对面之间往往不会自动对应需要额外处理。一种方法是对三组相对面使用“复制面”网格操作先在一个面上布置表面网格再复制到对应面最后生成体网格这样周期边界上节点严格一致。另一种方法是用扫描网格从底面扫到顶面自动保持对应关系。对于简单立方单元我用扫描网格的效果非常好自由度上升不大但各边界节点天然匹配。复杂几何比如倾斜孔洞扫描网格不好做就只能老老实实设置“一致网格”映射这在周期性结构仿真里非常实用推荐所有做微纳光子器件的人都掌握。5.4 三维模型内存不够求解时间太长三维光子晶体能带计算最烦人的就是内存。自由度数上百万后稀疏矩阵分解对内存的消耗非常吓人。几条实用经验第一在特征值求解器里限制搜素区间只求带隙附近模式不要竹篮打水捞200个高频模式。用目标特征值搜索配合特征值数量限制能显著减少变换矩阵的规模。第二网格收敛性先在一维路径上验证比如只算 Γ-X 一条线上的能带网格加密到结果基本稳定再在整条布里渊区路径上扫。第三如果算大结构优先用迭代求解器搭配合适的预处理方式比如 GMRES 配几何多重网格。当然 MUMPS 这种直接求解器在三维问题上更稳就看内存是否撑得住。我个人在反蛋白石模型上的经验是网格从 0.15a 降到 0.1a能带频率大约变化 1% 以内。如果只是定性看带隙位置验证结构参数0.15a 的默认网格已经够用不需要一开始就开到最高精度。5.5 后处理时如何快速过滤伪模伪模的过滤其实没有万能的公式但可以总结几个实用判据。首先看特征值的虚部真实无耗系统中的模式应该是实数频率虚部明显非零的大概率有问题。其次看场分布真实布洛赫模式场分布与 Bloch 相位规律一致伪模往往空间分布杂乱不符合周期结构的对称性。第三把同一波矢下的本征模场能量集中部位和介质分布对照伪模经常在网格缺陷处或者边界附近莫名其妙地集中。我在 COMSOL 里常用一个小技巧在后处理中用色标显示 Ez 或 H 场模然后按频率从小到大循环播放各模式。真实模式看起来像模像样的驻波或行波场伪模则像一团乱麻。虽然听起来有点经验主义但在工程实践中效率相当高。最后再分享一点个人体会我做了几年微纳光学仿真最大的感悟是COMSOL 里真正值钱的内核就是弱形式所有物理场接口的本质都是这套数学工具上的壳。刚开始算三维光子晶体能带时我也盯着物理接口的周期条件设置看了半天后来自己动手把弱形式推了一遍才发现原来 Floquet 周期条件并不是什么玄学它就是在布洛赫替换后的边界相位补偿。当你遇到 COMSOL 内置接口满足不了的模型时回到弱形式层面去修改往往反而能打开一条新的通路。以后无论你是想加非线性材料、磁光效应、增益介质还是做拓扑光子晶体里的规范场项弱形式都会成为你手里最趁手的工具。