拿到“COMSOL弱形式求解三维光子晶体能带结构分析”这个题目可能很多人第一反应是“又是个科研算例”。但做过光子晶体的人都知道能带结构算得准不准、边界条件加得对不对、三个方向的Bloch因子怎么作用远比“画出几条能带曲线”复杂得多。三维光子晶体的能带计算最磨人的不是麦克斯韦方程本身而是如何把H场本征问题以稳定的弱形式写出来再在Floquet周期条件下把布里渊区高对称路径跑完。这篇文章我就从头拆一遍为什么要用弱形式、弱形式方程怎么展开、COMSOL里具体怎么搭、后处理如何判读最后再把我踩过的几个坑盘一遍。适合正打算用COMSOL算三维光子带隙、却卡在“二阶方程边界不连续”或“特征值出来一堆虚数”的同学。1. 项目概况三维光子晶体能带到底怎么算1.1 要复现的物理图像光子晶体是一种介电常数呈周期分布的结构周期尺度在光波长量级。电子在周期性势场中会形成能带光子在周期性介电环境中同样会形成光子能带与带隙。三维光子晶体的典型算例是简立方晶格在每个格点放置一个高折射率介质球背景是空气。晶格常数取a球的半径设为r球材料折射率n3.5左右对应相对介电常数约12.25背景为空气n1。这样一套结构如果球的填充比选在合适区间会在某个归一化频率范围内出现三维完全带隙从任何方向入射的光子在这个频段都无法传播。能带结构分析要做的就是在一个单位晶胞上求解麦克斯韦方程得到频率与波矢的色散关系ω(k)。这个色散关系沿着第一布里渊区的高对称路径展开就是大家常说的能带图。三维问题比二维麻烦很多。二维光子晶体只需在xy平面做周期展开z方向是平板或均匀延伸边界条件相对简单三维则是三个方向都有周期边界而且每个方向上的Bloch相位都是复数。介电常数在球面处发生突变E场分量在界面处不连续如果直接用电场方程做强形式求解界面项处理起来会很痛苦。所以更稳妥的做法是用H场写方程H在界面处的法向连续性天然更好避免很多数值振荡。1.2 为什么选弱形式而不是现成物理场COMSOL的“电磁波、频域”接口确实提供了求解本征问题的能力也能搭配周期性边界条件。但内置接口按默认真实电场或磁场构造方程很多与材料界面相关的项被封装在内部用户要对方程做修改时插不上手。更关键的是三维光子晶体需要的是带相位因子的Floquet周期条件内置接口虽然也支持但对求解域、网格、特征值搜索范围的灵活度有限。弱形式PDE接口则完全不同方程按积分弱式显式输入。你写下的每一项就是有限元离散后每个方程的贡献。优点有三个方程完全可见H场的旋度项、介电常数倒数都能精确控制界面条件自动包含在分部积分之后只要被积函数写对介质突变处的弱解天然满足连续性边界约束可以自由定义Floquet相位关系可以直接写成约束表达式。代价是启动门槛高一点。但一旦把弱形式表达式理清后续修改材料模型、加入各向异性介质、甚至扩展到非线性光子晶体都只是加一项减一项的事。这也是我最终坚持用弱形式算三维光子晶体的原因。1.3 整条仿真链路完整流程可以分为四个阶段理论准备、模型搭建、参数扫描、后处理验证。理论准备要先把本征方程改写成H场形式然后用Bloch定理把无限周期问题限制到单个晶胞再把强形式改写成弱形式。模型搭建阶段在COMSOL里建立简立方晶胞几何定义材料参数与空间函数创建弱形式PDE接口输入三个磁场分量的弱方程再施加带Bloch相位的周期边界约束。第三个阶段是扫描布里渊区高对称路径上的k向量每给一组k就求解一次特征值得到对应频率。最后一个阶段把多组特征频率汇总绘制ωa/2πc随k路径变化的曲线并识别带隙。这条链路里最容易出问题的是弱形式方程的分量拆分和Floquet边界的相位设置。下面按环节讲细节。2. 原理层从麦克斯韦方程到弱形式表达式2.1 本征值方程与Bloch边界无源、非磁性、无损耗介质中麦克斯韦方程可以约化成磁场H的本征方程∇ × (1/ε ∇ × H) (ω/c)² H其中εε(r)是空间位置相关的相对介电常数ω是角频率。由于介质周期排列ε(rR)ε(r)R是晶格矢量。Bloch定理告诉我们周期介质中的场具有调幅平面波形式H(r) u(r) e^{i k·r}u(r)是周期函数k是第一布里渊区内的波矢。于是有限元计算只需在一个原胞内进行边界条件变成H(rR) H(r) e^{i k·R}这是带相位的周期边界也就是Floquet条件。在纪立方晶格中晶格矢量为a(1,0,0)、a(0,1,0)、a(0,0,1)因此三个方向的相位因子分别是e^{i k_x a}、e^{i k_y a}、e^{i k_z a}。2.2 弱形式的降阶与连续性含义把上面本征方程乘上测试函数V并在单位晶胞上积分∫ V · [∇ × (1/ε ∇ × H)] dV (ω/c)² ∫ V · H dV直接用这个式子做有限元需要二阶导数存在而且界面处1/ε发生跳变二阶导数意义不明。对左端做一次分部积分把旋度转移到测试函数上∫ (1/ε) (∇ × V) · (∇ × H) dV (ω/c)² ∫ V · H dV这就是弱形式。好处很明显场只需一阶可导介电常数跳变被吸收到体积分内部在介质界面分部积分自然给出法向连续性条件不需要额外处理强形式的界面方程。对H场来说界面处法向分量连续性天然满足进一步减少了界面伪振荡。如果令λ(ω/c)²弱形式写成∫ (1/ε) (∇ × V) · (∇ × H) dV - λ ∫ V · H dV 0在COMSOL特征值研究里λ就对应“特征值”变量。最终能带图画的是f ω/2π而ω c√λ注意单位换算即可。2.3 三维旋度弱式的分量拆分COMSOL弱形式PDE接口支持自定义因变量。三维情况下磁场有三个分量Hx、Hy、Hz对应三个测试函数Vx、Vy、Vz。关键是旋度场的交叉耦合∇×H的分量不是某个单独分量的导数而是不同分量的导数之差。如果不把弱式展开直接写成矢量形式COMSOL的标量弱形式编辑器并不认。需要把被积函数按测试函数拆成三个方程。设rotH_x ∂Hz/∂y - ∂Hy/∂z rotH_y ∂Hx/∂z - ∂Hz/∂x rotH_z ∂Hy/∂x - ∂Hx/∂y同理写出rotV的三个分量。完整弱贡献是(1/ε)( rotH_x·rotV_x rotH_y·rotV_y rotH_z·rotV_z ) - λ( Hx·Vx Hy·Vy Hz·Vz )把它按Vx、Vy、Vz分组得到三个独立的被积表达式方程对应分量弱形式被积表达式不含测试函数导数测试函数导数需自行展开Hx方程(1/ε)( rotH_y·∂Vx/∂z - rotH_z·∂Vx/∂y ) - λ·Hx·VxHy方程(1/ε)( -rotH_x·∂Vy/∂z rotH_z·∂Vy/∂x ) - λ·Hy·VyHz方程(1/ε)( rotH_x·∂Vz/∂y - rotH_y·∂Vz/∂x ) - λ·Hz·Vz这里把质量项λH·V分摊到各分量。实际写入COMSOL时测试函数Vx、Vy、Vz分别用test(Hx)、test(Hy)、test(Hz)表达导数是d(test(Hx),z)之类。之所以这样拆分是为了让每个因变量的方程都包含正确的旋度耦合如果写成强形式再乘以测试函数有限元刚度矩阵会出现非对称或刚度缺失的问题。3. COMSOL实操建模、边界与求解配置3.1 几何与材料参数定义建模前先在全局参数里定义一组变量参数数值含义a1晶格常数归一化r0.3球半径eps_b1背景介电常数eps_s12.25介质球介电常数几何用单位立方体表示晶胞长宽高都是a。球心放在晶胞原点位置。由于简立方晶格球心实际上在格点但球在相邻晶胞之间是周期排列的球体可以放在立方体中心也即一个晶胞的原点并让球体穿透边界。为了保证周期几何完整球体必须复制到六个相邻边界面上x-a/2与xa/2处各一半球体y、z方向同理。如果建模时只放一个中心球不处理穿出部分后面的周期约束会在边界上产生几何不连续网格剖分和边界映射都会出问题。为了避免穿出球体与边界面相交处出现悬空面建议用“阵列布尔并集”构造几何在中心位置放球然后把球按照±a三个方向平移后做并集得到一个含跨界球体的晶胞几何。这样每个边界上的球截面形状与相邻晶胞严格匹配周期边界才有映射基础。材料参数不直接给材料节点而是在变量里定义一个空间函数epseps eps_b (eps_s - eps_b) * (r^2 (x^2 y^2 z^2))这个表达式用球内外的判断构造介电常数分布避免额外设置材料域弱表达式里直接引用epsInv1/eps即可。注意这里的x、y、z是空间坐标变量COMSOL中默认即为变量名。3.2 弱形式PDE节点的表达式新建物理场时选择“数学 弱形式偏微分方程(weak form PDE)”。因变量设置为三个标量Hx、Hy、Hz。COMSOL允许在一个弱形式PDE接口里添加多个因变量但为了方便检查每个分量的方程我习惯添加三个单独的弱形式PDE节点分别命名为Hx方程、Hy方程、Hz方程每个节点只控制一个因变量。先定义中间变量rotHx d(Hz,y) - d(Hy,z)rotHy d(Hx,z) - d(Hz,x)rotHz d(Hy,x) - d(Hx,y)这些可以在“变量”节点预定义也可以在弱表达式里直接写。然后三个H方程节点的弱贡献分别写成Hx方程epsInv * ( rotHy * d(test(Hx),z) - rotHz * d(test(Hx),y) ) - lambda * Hx * test(Hx)Hy方程epsInv * ( -rotHx * d(test(Hy),z) rotHz * d(test(Hy),x) ) - lambda * Hy * test(Hy)Hz方程epsInv * ( rotHx * d(test(Hz),y) - rotHy * d(test(Hz),x) ) - lambda * Hz * test(Hz)其中epsInv 1/eps是在变量节点里定义的空间函数。lambda是COMSOL特征值求解器默认的本征值变量其物理量纲为(ω/c)²。这里要注意COMSOL的test函数写法在较新版本中测试函数写作test(u)测试函数的导数写作d(test(u),x)。如果版本比较旧可能用u_test、ux_test之类的记号需要根据实际版本调整。写完弱表达式之后还要在“因变量”设置里把Hx、Hy、Hz的“复数”属性打开。Floquet边界涉及e^{ikx}复数相位所有磁场自由度都必须是复数否则无法描述传播波解。3.3 Floquet周期边界条件处理这是整个建模里最容易被卡住的环节。基本思想在x0与xa的对应面、y0与ya的对应面、z0与za的对应面之间建立映射约束H(a方向目标面) H(0方向源面) · e^{i k·R0}其中R0是该方向上的晶格平移矢量。在COMSOL的弱形式PDE接口中没有现成的“Floquet边界”按钮需要手工施加。我的做法是为每个方向添加“边界对”映射然后定义“约束”节点。具体来说以 x 方向为例选中x0面作为源边界xa面作为目标边界在模型开发树的“定义 映射”中建立一个边界映射把目标和源按“相反法向、对应坐标”对应起来在物理场部分的“约束”节点中写入Hx_tgt - Hx_src * exp(ikxa) 0Hy_tgt - Hy_src * exp(ikxa) 0Hz_tgt - Hz_src * exp(ikxa) 0kx是预先定义好的参数取值来自布里渊区扫描。同理定义y、z方向的约束并把相位因子换成exp(ikya)和exp(ikza)。需要注意约束中的Hx_tgt与Hx_src分别指目标边界和源边界上的因变量值。COMSOL源边界与目标边界的标记方式取决于映射的方向如果写反了相位因子变成共轭能带会整体偏移甚至出现错误简并。识别方法是检查约束是否满足“互换源目标后方程只差一个共轭”如果映射方向反了通常特征频率会出现负值或虚部异常。3.4 特征值研究、参数扫描与网格研究类型选择“特征值”。COMSOL特征值变量lambda默认是复数但由于我们的方程在无损耗介质中特征频率是实数求解时可以只需要实数部分。在求解器设置里将特征值搜索方式设为“靠近目标”目标值设为某个预测区间例如归一化频率0.3附近如果扫描频率范围比较大就设成“全部特征值”让求解器返回任意数目的最小特征值再人工筛选。计算能带时需要对布里渊区高对称路径逐点扫描。简立方晶格通常取Γ(0,0,0) → X(0.5,0,0) → M(0.5,0.5,0) → Γ(0,0,0) → R(0.5,0.5,0.5) → X(0.5,0,0) → M(0.5,0.5,0) → R(0.5,0.5,0.5)坐标是用2π/a归一化的k坐标。为减少重复计算一般只取Γ-X-M-Γ-R-X-M-R这段路径上的部分点。建议在全局参数里分别定义kx、ky、kz作为扫描参数然后在“扫描”设置里建立一个包含各段路径的参数表让三个参数同时变化。COMSOL的“辅助扫描”可以按参数表逐点计算每个k点求解一次特征值也可以一次性把所有k点加入参数扫描但三维模型下内存占用会很大建议分批求解每批扫一段路径。网格方面边界或曲面处用自由四面体球面附近局部细化。由于介质折射率对比高12.25比1界面附近的场变化剧烈球体表面至少需要覆盖三层网格。晶体内部区域网格可以稀疏一点但球边界处最小单元尺寸建议小于0.05a。网格单元数通常会在几万到几十万之间取决于精度要求。粗网格能跑出趋势但带边缘位置误差可能达到10%以上正式出图前最好做一次网格收敛性验证。求解器方面推荐MUMPS内存消耗可控且稳定。特征值求解时如果报告“求解失败”或“未找全特征值”优先检查目标值范围是否覆盖了所需频率区间再把“特征值数”适当调大。4. 后处理能带图绘制与判读4.1 数据导出和绘图每个k点求解完成后COMSOL派生值里可以查看特征频率列表。把扫描得到的结果导出成文本格式大概是每一行对应一个k点、一组特征频率。我习惯用Python脚本整理数据横轴是沿高对称路径的累积k距离纵轴是归一化频率ωa/2πc。因为不同k段所对应的实际波矢长度不同X-M段的长度是M-Γ段的一半不能直接把点号当作等距坐标要按K点之间的欧氏距离累加。归一化频率的计算如果COMSOL返回的特征值是lambda那么ω c·sqrt(lambda)。在a1且c设为单位的情况下纵轴直接取sqrt(lambda)/2π即可。我在实际中会把a、c全部显式写上避免忘记单位导致曲线整体平移。4.2 带隙、态密度和结构参数的影响能带图出来后首先要找带隙在某条能带的最低点与下一条能带的最高点之间是否存在一个频率区间使所有路径上都没有能带穿过。如果这个区间在所有k方向都为空就是完全带隙如果只在部分方向存在则是方向带隙不能完全阻止传输。三维简立方介质球结构的典型规律是折射率对比越高带隙越容易出现填充比f4πr³/(3a³)存在一个最佳区间。r/a在0.3到0.35附近时通常能得到相对宽的三维带隙。填充比太小散射太弱带隙闭合填充比太大介质球几乎联通等效均匀介质带隙也会消失。算完后可以改变r/a重算几条能带看带隙随填充比的变化趋势这是论文里常见的参数图。结构参数优化的核心是球半径。同样折射率下不同r会改变能带的带边位置。把所有参数扫描结果叠加到同一张图能直观看到带隙从哪里开启、在哪里闭合。这里建议不要只算对称点M到R之间往往还藏有决定带隙上限的能带极值点忽略中间点会把带隙宽度算大。4.3 验证和经典结果比三维简立方介质球能带是最经典的光子晶体算例之一很多文献都有现成数据。完成一次完整扫描后先别急着优化结构拿r0.3、ε12.25的结果和文献曲线对比差距在2%以内说明模型基本正确差距很大则需要检查是否用错了周期性边界的相位因子kx、ky、kz是否混用是否漏了某个方向的Bloch约束网格是否足够细。我自己的经验是如果Γ点附近第一条能带的频率和文献差超过5%基本都是边界约束写错了只有高频带第五条以上误差大才是网格问题。5. 常见问题与排查技巧实录5.1 伪模态与虚特征值特征值求解结果中混入虚数或实部很小的模态是很常见的情况。伪模态的来源通常是约束不完整某个面的Floquet相位忘了加导致该边界上出现不合理的反射模式网格太粗尤其是球面附近曲面离散误差产生局部振荡弱形式质量项符号写错λH·V项的符号如果和刚度项同号特征值会变成负数或虚数。排查方法很简单把虚部较大的模态对应的场分布画出来看能量是否聚集在局部区域。如果能场局域在某个边界上多半是边界约束缺失如果能场在球面上出现细密条纹多半是网格问题加密后条纹消失的就是伪模态加密后仍在的才是真实表面模态。5.2 边界映射设错的症状Floquet条件里源面与目标面的映射方向是整个计算最容易出错的地方。症状包括能带在Γ点之外出现劈裂但理论上该路径上存在二重简并特征频率在kx很小的时候随k变化异常剧烈三个独立方向的约束都能满足但结果不具备晶格对称性。一旦怀疑映射设错先把kxkykz0的极限情况跑一遍。这时Bloch相位退化为1边界应等效为普通周期边界。如果计算结果与直接使用标准周期约束的结果一致说明映射方向对。反之可通过交换源目标面并重新计算逐方向排查。5.3 低对称点简并与网格退化三维光子晶体能带在某些高对称点上存在简并比如Γ点的某些模态会因晶格对称性出现多重简并。但有限元网格如果剖分不对称会导致本应简并的模式分裂出微小频率差破坏能带图的对称性。遇到简并分裂先看网格剖分是否保持几何对称。如果网格不对称需要手动控制网格分布对球周围区域采用对称细化或者在模型里保留对称面并做镜像网格。此外特征值求解到简并模态时求解器返回的本征向量可能混叠表现为两个本征值非常接近但场分布交替重叠。这属于正常现象绘制能带图时合并相近点即可。5.4 特征值搜索设置建议我先给一组我自己常用且较稳的特征值设置求解器特征值目标值设为预期带隙附近的λ所需特征值个数每次8到12个最低能带通常只需5到6个但多算几个可以防止漏带特征值搜索方法靠近目标值区间不宜太大太大容易混入伪模态复数处理保留所有返回特征值先看虚部占比虚部大于1e-6的剔除再从实数特征值里排查伪模态。如果扫描整个布里渊区时发现某个k点附近特征值不连续大概率是这个k点处所需特征值个数不够导致一条能带分支没被算出来。此时单独在该点增加特征值数量确认是否有遗漏。5.5 三维大网格下的内存优化技巧三维计算中内存消耗是瓶颈。我有几个降内存的经验先做二维截面对比在xy平面做二维光子晶体能带验证弱形式表达式和边界约束逻辑正确再扩展到三维扫描k点时不要把所有k值放进同一个参数扫描任务每段路径单独跑利用对称性只算部分路径。例如简立方结构在kxkykz方向段和另外几个方向段有对称等价关系但前提是结构本身各向同性对称采用“求解后继续”模式每个k点求解结束后清理旧的解数据避免累积占用。网格阶数方面二阶单元通常比一阶单元精度高很多但在三维大网格里二阶也会显著增加自由度。如果内存紧张先默认一阶网格跑通整个流程再在高精度需求点局部升阶。结尾一点体会这个弱形式项目做下来我最大的感受是最难的不是写弱方程本身而是让周期边界施加得“绝对干净”。三维Bloch边界稍有偏差带隙闭合得毫无规律让人满头雾水。后来我把验证分成三步先用二维截面对照解析结果再把三维结果在Γ点极限对照普通周期边界最后才信任完整的三维能带图。现在如果重新做一遍我会把kx、ky、kz三个目标方向分成三个独立的“定义”变量并在约束表达式中分别命名方便在出问题时单独注释掉某一方向做定位。另外提醒一下不同COMSOL版本对复数因变量和特征值相位因子的支持细节有差异遇到版本特有的报错优先查弱形式接口帮助文档比满网搜答案更快。最后分享一个细节球体表面网格细化时不要只加边界层最好同时给球面内部一个低生长率的自由四面体分布否则高频能带边缘会莫名其妙出现锯齿。祝顺利跑出漂亮的能带图。