做光子晶体能带计算的人十次里有七八次都会先碰壁MPB、MEEP、Lumerical这些专用工具各有各的语法和安装成本好不容易跑通一个二维六角格子的带隙真要转到三维结构的时候又得重新折腾一遍。其实COMSOL完全可以扛下这个活而且不止一条路。这篇文章我聚焦“三维光子晶体能带结构分析”这个硬骨头重点讲讲怎么用手写弱形式的方式在COMSOL里把Maxwell方程的特征值问题变成可解的能带扫描再把整个推导、建模、边界条件、特征值配置、布里渊区路径扫描和后处理串起来。适合读者正在用COMSOL做电磁场仿真、但对弱形式一直“敬而远之”的人已经会用内置波动光学模块算二维能带、但升级到三维后不知道边界条件该怎么处理的人以及想搞清楚“弱形式到底是怎么回事、为什么COMSOL能这么散”的研究生和工程师。源码级别的参数和表达式我会尽量写细你可以照着直接抄作业也能从中看到哪些地方需要根据自己的结构调整。1. 三维光子晶体能带的物理图像与数学框架1.1 光子晶体是什么能带图在表达什么光子晶体本质上是介电常数在空间上周期性调制的结构。电子在周期势场中运动产生电子能带光子在周期介电结构中传播同样会产生能带结构这就是光子带隙的来源。周期性结构里光子的色散关系不像在均匀介质中那样是一条直线而是被折叠成一系列能带分支某些频率区间可能完全没有对应的传播模式这就是完全带隙。做三维光子晶体能带分析核心就是算出一个频率对波矢的关系曲线横轴是第一布里渊区中选定的高对称方向路径縱轴是归一化频率。这个关系一出来带隙位置、群速度、慢光效应、缺陷模频率范围都能从这个图里读出来。很多人一开始会习惯性地把思路停在二维算一个二维圆孔的带隙多简单扫描路径也就五六段。但真到了三维问题就不一样了。光子晶体在x、y、z三个方向都有周期性边界条件必须在三对相对面同时建立波的传播方向与倒格矢空间的关系也更复杂。更麻烦的是电场的每个分量都耦合在一起特征值问题的自由度要比二维大一个数量级。1.2 布洛赫定理与第一布里渊区三维场景的核心坐标系周期结构中介电常数满足[ \varepsilon(\mathbf{r}) \varepsilon(\mathbf{r}\mathbf{R}) ]其中 (\mathbf{R}) 是晶格平移矢量。根据布洛赫–弗洛凯定理电场可以写成[ \mathbf{E}(\mathbf{r}) \mathbf{u}_k(\mathbf{r}) e^{i\mathbf{k}\cdot\mathbf{r}} ]其中 (\mathbf{u}_k) 是与晶格同周期的函数(\mathbf{k}) 是倒空间波矢。这意味着在一个原胞的边界面上场不是直接相等而是差一个相位因子 (e^{i\mathbf{k}\cdot\mathbf{R}})。这正是COMSOL里Floquet周期性边界条件的数学基础也是三维光子晶体建模中最容易被忽略、最容易出错的地方。第一布里渊区是倒空间中的Wigner-Seitz原胞。对简单立方晶格倒格矢为[ \mathbf{b}_1 \frac{2\pi}{a}\hat{x}, \quad \mathbf{b}_2 \frac{2\pi}{a}\hat{y}, \quad \mathbf{b}_3 \frac{2\pi}{a}\hat{z} ]第一布里渊区就是 (-\pi/a \le k_i \le \pi/a) 的六面体。能带图上常用的几个高对称点坐标如下表以倒格矢坐标为 ((k_xa/2\pi, k_ya/2\pi, k_za/2\pi)) 表示高对称点归一化坐标物理意义Γ(0, 0, 0)布里渊区中心X(0.5, 0, 0)面心方向M(0.5, 0.5, 0)棱中点方向R(0.5, 0.5, 0.5)角点方向实际扫描时一般沿 Γ→X→M→Γ→R→X 这样的闭合路径取点每一个k点求解一次特征值问题。1.3 为什么弱形式在这里是绕不开的核心材料领域很多人对“弱形式”有天然排斥觉得数学门槛高不如直接用现成物理场接口。但COMSOL的有限元求解器本质上就是在每个单元上组装弱形式矩阵。内置的“电磁波频域”接口内部也是一套写好的Maxwell方程弱形式只不过COMSOL替你封装好了。问题在于光子晶体能带分析经常会遇到超出内置接口能力的需求非线性材料、强各向异性介电张量、磁光效应、增益损耗分布、空间色散、特殊边界约束。这些场景下内置接口不一定能直接修改而弱形式PDE接口几乎可以写任意形式的方程给了最大的自由度。另一个原因是性能。三维光子晶体单元虽然小但特征值求解需要算足够多的模式网格自由度动辄几十万甚至上百万。弱形式接口能让你精确控制被积函数的结构去掉内置接口中一些没必要的后处理变量和自动模式选择逻辑实测下来内存占用会低一些也更容易针对具体问题时做一些矩阵结构上的调整。2. 弱形式方程推导从Maxwell方程组到COMSOL表达式2.1 特征值方程的两种写法以及一种更可靠的电磁形式忽略自由电荷和传导电流线性无损非磁性介质中频域的Maxwell方程组可以化简为电场旋度方程[ \nabla\times(\nabla\times\mathbf{E}) \frac{\omega^2}{c^2}\varepsilon_r(\mathbf{r})\mathbf{E} ]这是电场形式。还有对应的磁场形式[ \nabla\times\left(\frac{1}{\varepsilon_r}\nabla\times\mathbf{H}\right) \frac{\omega^2}{c^2}\mathbf{H} ]磁场形式的好处是无散条件 (\nabla \cdot \mathbf{H}0) 自动满足但代价是方程中出现了 (1/\varepsilon_r)在介电常数跳变的界面处处理起来略微麻烦。对COMSOL弱形式接口而言写电场形式更直接需要手动处理的伪模问题我后面专门讲。这里先把特征值写清楚。令[ \lambda \left(\frac{\omega}{c}\right)^2 ]那么电场形式就是一个标准的广义特征值问题[ \nabla\times(\nabla\times\mathbf{E}) \lambda \varepsilon_r \mathbf{E} ]COMSOL特征值节点输出的特征是 (\lambda)最后通过 (\omega c\sqrt{\lambda}) 换算回频率。2.2 弱形式推导完整过程为什么要把二阶偏导拆开直接对上面的方程做伽辽金加权余量选测试函数 (\mathbf{F})两端点乘再在全域积分[ \int_V \mathbf{F}\cdot\left[\nabla\times(\nabla\times\mathbf{E})\right]dV\lambda\int_V \varepsilon_r \mathbf{F}\cdot\mathbf{E},dV ]这里的问题是方程左边含二阶导数。有限元离散后如果直接对这个形式做节点值是二阶数值精度会很差。弱形式的精髓在于利用分部积分把二阶导数降为一阶换来的代价是需要处理边界项。用矢量恒等式[ \nabla\cdot\left[\mathbf{F}\times(\nabla\times\mathbf{E})\right](\nabla\times\mathbf{F})\cdot(\nabla\times\mathbf{E})\mathbf{F}\cdot\nabla\times(\nabla\times\mathbf{E}) ]移项并对体积积分利用散度定理[ \int_V \mathbf{F}\cdot\nabla\times(\nabla\times\mathbf{E}),dV\int_V (\nabla\times\mathbf{F})\cdot(\nabla\times\mathbf{E}),dV\oint_S \mathbf{n}\cdot\left[\mathbf{F}\times(\nabla\times\mathbf{E})\right]dS ]于是弱形式变成[ \int_V (\nabla\times\mathbf{F})\cdot(\nabla\times\mathbf{E}),dV\lambda\int_V \varepsilon_r \mathbf{F}\cdot\mathbf{E},dV\oint_S \mathbf{n}\cdot\left[\mathbf{F}\times(\nabla\times\mathbf{E})\right]dS ]右侧边界积分在Floquet周期条件下自动消失因为相对的两个周期面上场满足相位匹配边界积分的贡献在求和时相互抵消。这样最终的弱形式表达式只剩两个体积分。2.3 三维矢量场景下弱表达式的COMSOL写法COMSOL弱形式PDE接口中测试函数写作test(...)。假设因变量用u、v、w分别代表 (E_x)、(E_y)、(E_z)三维旋度的内置算子为curl(u,v,w)那么上面弱形式直接写成curl(test(u),test(v),test(w))*curl(u,v,w) - lambda*epsr*(test(u)*utest(v)*vtest(w)*w)注意第一项是两个矢量的点积COMSOL中写成dot(curl(...), curl(...))更规范不过简写*在标量点积场景下也能识别。这里的epsr是相对介电常数空间分布。对简单立方介质球结构可以用坐标函数直接定义r_sp sqrt((x-0.5)^2(y-0.5)^2(z-0.5)^2) epsr if(r_sp 0.3, 12.25, 1)假设球半径0.3、背景折射率1、球的折射率3.5。这就是一根表达式不用在不同域分别设置材料后续扫描参数也方便。球心放在立方体中心整体几何是边长为1的正方体。这是三维矢量问题的关键与二维问题可以拆分TE/TM极化不同三维情况下三个电场分量完全耦合必须把它们都作为因变量放进去。上式那一个弱表达式在COMSUL中会自动生成三个分量方程分别对test(u)、test(v)、test(w)取变分。3. COMSOL落地几何、Floquet周期条件与特征值配置3.1 三维模型的几何搭建与材料变量的定义打开COMSOL新建一个三维模型。在“几何”节点下直接用块Block建立边长为1的立方体中心顶点设置在(0,0,0)。然后添加一个球Sphere球心放在(0.5,0.5,0.5)半径0.3。这样球就在立方体的中心区域。注意COMSOL里当后续要用布尔运算时“形成联合体”和“形成装配”都可但弱形式接口下仿真域一般我还是选“形成联合体”避免装配后额外的接触对复杂度。几何建好后在“全局定义”下建几个常量a 1 // 晶格常数无量纲化 R_sphere 0.3 // 球半径相对晶格常数 eps_bg 1 // 背景相对介电常数 eps_sphere 12.25 // 球的相对介电常数 c0 1 // 光速无量纲化最终计算的是归一化频率然后在“变量”里定义空间介电常数分布。COMSOL变量支持if表达式。需要特别注意if里判断条件在单元边界的连续性处理在弱形式里直接用一般没问题但网格在球表面附近必须足够细否则介电常数跳变会被模糊化。3.2 三个方向的Floquet周期条件设置对这个三维光子晶体能带问题周期性边界条件必须同时施加在三对相对面上。很多人算二维没有问题一升到三维就开始反复报错或者算出明显不合理的能带多半就是只设了一个方向的周期条件。最推荐的做法是使用“电磁波频域”接口的内置周期性条件。在物理场下右键添加“周期性条件”选择边界时把x方向相对的两个面分别设为源和目标类型选“Floquet周期”波矢分量按照kx、ky、kz填入。然后重复三次分别处理y方向两对面和z方向两对面。这里有一个极易踩的坑COMSOL的Floquet周期条件中的波矢单位。通常几何尺寸用米为单位时波矢单位是rad/m。但在光子晶体能带计算中我们往往做无量纲化几何设置为边长a1波矢就应该直接填归一化倒格矢坐标乘以(2\pi/a)。如果你之前算过二维会发现很多人喜欢直接用0到1的归一化参数然后在边界条件里乘上2*pi/a这个习惯在三维也建议保留避免在路径扫描时写一长串带π的表达式。如果坚持用“弱形式PDE”接口而不是内置电磁模块注意弱形式PDE接口不一定直接暴露Floquet周期性条件节点。这时两个选择一是改用“系数形式PDE”或“一般形式PDE”接口这两个接口带有周期性条件可以选Floquet类型方程特征值形式也支持二是自己加约束方程但操作繁琐且数值稳定性差。我的建议很直接别硬扛三维能带用内置接口跑通再回头研究自定义弱形式。3.3 弱形式PDE接口的因变量设置与复数开关如果你明确就是要用纯弱形式接口来做建模时有三个关键设置。第一因变量数。在“弱形式PDE”节点下因变量栏把三个分量全部启起来分别命名为u、v、w。也可以用Ex命名但COMSOL内置的curl算子在向量模式下习惯接受u,v,w三个分量建议保留默认形式避免混淆。第二复数开关。Floquet周期条件中波矢都是复数相位因此因变量必须设为复数。COMSOL中弱形式PDE默认因变量是实数需要在“因变量”属性里勾选复数。这个细节一旦漏掉会在计算时出现奇异行为或者在求解过程中直接报“复特征值不适合实因变量”之类的错误。第三弱表达式。在前面推导的基础上在“弱形式PDE”节点里的“弱表达式”栏输入dot(curl(u,v,w), curl(test(u),test(v),test(w))) - lambda*epsr*(u*test(u)v*test(v)w*test(w))这里epsr就是前面定义的空间介电常数变量。若你用的是“一般形式PDE”则还应在“源项”等位置做对应变化。不管哪个接口原理上都是上面推导出来的同一套弱形式。3.4 特征值研究设置与网格策略模型建好后添加一个“特征值”研究。COMSOL特征值研究中默认求解的特征变量就是lambda。对光子晶体能带每个波矢点需要求解前N个模式N取决于你想看清几条带。三维结构通常建议先取8至10个。特征值求解器设置中一个影响很大的参数是“特征值搜索基准点”。默认情况下COMSOL会在0附近搜索这意味着会找到大量零或接近零的伪模特别在弱形式接口下更明显。建议把特征值搜索区间设置为一个合理的正值范围比如从0.01开始这样自动跳过一部分零模。网格是三维能带是否准确的决定因素。我的经验是球体表面和内部要加密背景区域可以相对稀疏。用自由四面体网格球半径0.3的结构下球体表面的最小单元尺寸0.03到0.05能给出不错的低频能带收敛性如果要算较高频分支必须再加密同时网格数会迅速增长。建议先跑粗网格看前几带的位置大致在哪个频率范围再决定加密尺度。求解完成后COMSOL会在全局数据中给出每个波矢点对应的特征值lambda。归一化频率按公式freq_norm sqrt(lambda)/(2*pi)计算因为我们已经把c和a都归一化为1了。4. 布里渊区路径扫描与能带曲线的后处理4.1 参数化扫描的路径定义技巧能带曲线的横轴是沿高对称路径的波矢位置通常用一个参数s从0扫描到5来表示路径的5个线段。以简单立方为例路径为线段1Γ (0,0,0) → X (0.5,0,0)线段2X (0.5,0,0) → M (0.5,0.5,0)线段3M (0.5,0.5,0) → Γ (0,0,0)线段4Γ (0,0,0) → R (0.5,0.5,0.5)线段5R (0.5,0.5,0.5) → X (0.5,0,0)每一段内用一个局部参数t从0到1过渡全局参数s 线段号 t。波矢表达式用if函数嵌套写kx if(s1, 0.5*s, if(s2, 0.5, if(s3, 0.5*(3-s), if(s4, 0.5*(s-3), 0.5)))) ky if(s1, 0, if(s2, 0.5*(s-1), if(s3, 0.5*(s-2), if(s4, 0.5*(s-3), 0.5*(s-4))))) kz if(s4, 0, 0.5*(s-4))注意这里的波矢已经是归一化坐标乘以 (2\pi/a) 后才进入Floquet边界条件。如果你在边界条件里写波矢为kx*2*pi/a、ky*2*pi/a、kz*2*pi/a那么扫描范围就是s从0到5步长0.05或0.02。步长太小单次求解次数多太大能带曲线会不够光滑工程上先0.05跑通再局部加密。在COMSOL“研究”节点下添加“参数化扫描”选择全局参数s填入范围range(0,0.05,5)。此时每扫描一个s就对应一个新的波矢点执行一次特征值求解。4.2 从特征值结果绘制能带曲线计算完成后在“结果”里新建一维绘图组。关键是X轴设为comp1.s或直接写sY轴设为归一化频率。COMSOL特征值研究的结果中每个s点下有N个特征值。你可以直接使用“全局”绘图数据系列里选择按“本征值”索引绘制这样每个s点会产生N个点N条带就显示为N条曲线。注意不同特征值分支会自动按大小排序所以曲线之间会有交叉感的连接实际能带就是这样带交叉位置是简并点真实物理里多数是允许的简并不要以为是bug。如果觉得COMSOL绘图不够灵活可以在“派生值”节点里提取所有s点的特征值导出到表格再用Origin或Matplotlib画图。导出时记得同时导出s、lambda、以及kx、ky、kz方便后续标注高对称点位置。高对称点对应的s值要自己在图上标竖线s0对应Γs1对应Xs2对应Ms3对应Γs4对应Rs5对应X。4.3 能带图上如何识别并剔除伪模这是三维光子晶体计算中最影响使用体验的一步。内置电磁波模块用了适合电磁场问题的单元伪模情况相对可控但在弱形式PDE接口下由于节点元自由度设置原因数值上会引入相当多物理上不存在的模式。伪模通常有两个明显特征第一特征值非常小或等于零。它们对应的是满足旋度为零但梯度非零的“静态解”模式物理上没有传播意义。画出能带后在频率低到接近零的地方会看到一堆几乎平躺的线这些不是真实能带。第二在后处理中电场分布存在明显的非物理散度。内置接口由于边界条件设计得好一般不会出现这种问题弱形式接口就必须自己把关。处理伪模的办法简单粗暴地将频率低于某个阈值比如归一化频率0.05的数值点过滤掉。这个方法足够应对大多数情况。如果还想更严格可以在后处理中计算每个模式在整个原胞体积内 (\nabla \cdot (\varepsilon_r \mathbf{E})) 的积分超过某个比例就打上非物理标记绘图时排除。这样得到的能带才是干净的。5. 常见问题排查与实战心得5.1 特征值求解不收敛、负特征值从哪来三维电磁特征值问题经常遇到“算出一大堆负特征值”的情况。原因在于原始电场形式算子 (\nabla\times\nabla\times) 的零空间非常大所有纯梯度场都是零特征值模态。有限元离散后由于数值误差这些零模会被挤到负数一侧特征值求解器会把它们一并算出来。解决办法分几步。一是增大特征值搜索的目标频率不要在0附近找。二是增大所需特征值数量把你关心的那几条带放在求解结果的中间段而不是最开头。三是如果伪模已经严重污染前几个真实模式考虑改用磁场形式或者在内置“电磁波频域”接口下做计算。别在弱形式接口下硬扛已经有无数人在这上面浪费过时间。5.2 能带曲线在某个k点突然断裂或跳动这条一般不是计算错误而是特征值排序问题。每个s点求解N个特征值COMSOL自动按特征值大小排列但当两条带在某个k点靠近或交叉时模式在下一s点可能互换顺序表现在绘图上就是曲线跳变。物理上如果结构本身存在带交叉这是对的不需要干预。真正需要警惕的是另一种情况某个s点算丢了一条带。这常常是因为特征值数量N取得不够导致高频率分支在该点被漏算。解决方案很简单增大N。一般先N8如果绘图时发现高频分支出现空洞就要加大到10或12。5.3 网格密度、介电常数跳变与计算量平衡三维各向同性介质球加背景的结构球表面是介电常数跳变面网格必须足够细。我做过一个测试边长为1、球半径0.3、背景介电常数1、介质球介电常数12.25的结构当球表面最大单元尺寸从0.1减到0.025时前三个本征频率相对变化约3%到5%。粗网格的结果只能看趋势不能用于带隙精确定位。另一个容易被忽略的问题是坐标原点。布洛赫波矢的相位是在Floquet边界条件里定义的几何位置应和波矢表达式保持一致。比如球心在(0.5,0.5,0.5)时边界条件相对面之间的距离才有意义如果把球放在原点附近但波矢计算时没有加上对应的位置偏移会出现能带明显不对称的奇怪结果。5.4 内置模块和自定义弱形式如何选一个对比表最终选型我根据自己的使用经验给出一个直白的对比对比项内置电磁波频域接口自定义弱形式PDE周期性条件自带Floquet周期使用方便视接口而定弱形式接口需要额外处理伪模抑制较好需要自行过滤容易混入静态模自定义材料有一定限制几乎无限制上手难度低高需要理解弱形式推导适合场景常规光子晶体能带、缺陷模特殊物理模型、科研扩展我的个人经验是第一套能带计算无论如何先用内置接口跑通整个流程确认结构、参数、扫描路径都是对的再尝试自定义弱形式。自定义弱形式的真正价值在于给你一双“解读COMSOL黑箱”的眼睛而不是取代内置模块。如果你后面打算做缺陷模、波导耦合、非线性频率转换那么弱形式接口的灵活性会派上大用场。但三维能带这一步稳健永远是第一优先级。先让能带图漂亮地出来再去追求自定义带来的控制感。这个顺序踩过几次坑的人都懂。