Comsol仿真研究：蜂窝晶格光子晶体能带结构及陈数拓扑的MATLAB与MPH脚本实现

Comsol计算蜂窝晶格光子晶体能带拓扑陈数 包含mph与matlab脚本六边形蜂窝结构的光子晶体总能玩出意想不到的拓扑特性。最近在折腾ComsolMatlab联合计算陈数的流程发现几个关键操作手册上没写的坑位顺手记录点野路子经验。1. 建模六边形蜂窝结构怎么画COMSOL建模最省事的办法是用阵列复制功能。先画单个介质柱参数化半径r和晶格常数a这里容易踩的坑是旋转角度设置% 直接在Comsol的模型参数定义 a 1e-6; % 晶格常数 r 0.3*a; % 介质柱半径 theta [0, pi/3, 2*pi/3]; % 蜂窝基矢方向然后创建三个矩形阵列分别沿三个基矢方向偏移半个周期。重点是要勾选保持原始对象选项否则会生成不连续的蜂窝结构。网格划分建议用自由四面体边界层注意在介质柱表面设置至少5层边界网格。Comsol计算蜂窝晶格光子晶体能带拓扑陈数 包含mph与matlab脚本2. 能带计算扫参怎么不翻车使用Floquet边界条件时端口方向必须与倒格子矢量对齐。这里有个隐藏设定扫描波矢范围要覆盖整个不可约布里渊区。建议用参数化扫描kx和ky比如kx linspace(-pi/a, pi/a, 50); ky kx.*sqrt(3)/2; % 六边形结构的ky范围更窄计算结果导出成txt时注意检查数据列的顺序。曾经因为列顺序搞反导致后续计算全错建议导出时添加注释行% kx ky freq1 freq2 ... freqN3. 陈数计算Matlab暴力积分法读取Comsol导出的频率数据后先要构造贝里曲率场。这里有个偷懒技巧——用离散差分代替微分function C ChernNumber(kx, ky, berryCurvature) [dky, dkx] gradient(berryCurvature, ky(1,:), kx(:,1)); C sum(sum(dkx - dky)) * (kx(2)-kx(1)) * (ky(2)-ky(1)) / (2*pi); end实际计算时发现直接用本征态波函数计算联络更稳定。核心循环长这样for i 1:Nk for j 1:Nk % 计算相邻k点的波函数交叠 U psi(:,:,i,j) * psi(:,:,circ(i1,Nk),j); V psi(:,:,i,j) * psi(:,:,i,circ(j1,Nk)); berryCurvature(i,j) angle(U*V/(V*U)); end end其中circ是处理循环边界的自定义函数。注意波函数相位必须做规范固定否则计算结果会抽风。4. 翻车实录曾因COMSOL导出的频率单位是Hz而Matlab脚本默认用归一化频率导致量纲混乱六边形结构的高对称点必须包含K和K点否则陈数计算结果少0.5并行计算时发现贝里曲率在GPU加速下会出现符号错误建议先用单线程调试最后验证时拿真空和完全光子带隙的情况做基准测试。当蜂窝结构的Chern Number跳变到±1时意味着出现了拓扑边界态——这时候就可以去冲杯咖啡庆祝了。