Comsol弱解法在三维光子晶体能带计算中的应用

Comsol弱形式求解三维光子晶体能带。

深夜两点盯着屏幕上扭曲的能带曲线，突然意识到三维光子晶体的数值模拟就像在量子迷宫里玩俄罗斯方块——每个晶格参数都可能让整个能带结构瞬间崩塌。传统界面操作总让我感觉戴着镣铐跳舞，直到某天偶然翻到COMSOL的弱形式文档，就像突然拿到了打开量子围栏的钥匙。

在模型树的"定义"节点右键新建变量，敲下这段代码时手指有点发抖：

epsilon_r = 4.0*(sin(2*pi*x/a)^2 + sin(2*pi*y/a)^2 + sin(2*pi*z/a)^2) + 1.0; mu_r = 1.0;

这其实是在构造三维光子晶体的介电常数周期性分布。当晶格常数a=400nm时，这个表达式会在空间中生成类似钻石结构的介电分布，就像用数学公式捏造纳米尺度的光学乐高积木。

真正的魔法发生在弱形式PDE设置里。在电磁场控制方程中，我习惯性地把麦克斯韦方程改写为：

% 弱形式表达式 test(Ex)*(1/mu_r*(curlHz - curlHy)) + test(Ey)*(1/mu_r*(curlHx - curlHz)) + ... test(Ez)*(1/mu_r*(curlHy - curlHx)) - omega^2*epsilon_r*(test(Ex)*Ex + test(Ey)*Ey + test(Ez)*Ez)

这段看似凌乱的代码其实在偷偷帮我们处理电磁场的切向连续性。特别是当处理光子晶体界面处的场匹配时，弱形式会自动完成传统伽辽金法需要手动处理的积分项，就像有个隐形的场向导在帮忙拼接电磁场碎片。

Comsol弱形式求解三维光子晶体能带。

设置Floquet周期边界时的参数扫描最让人头疼。记得在特征频率研究中设置波矢k沿着布里渊区边界扫描，比如从Γ点到X点的路径：

kx = linspace(0,pi/a,20); ky = 0; kz = 0;

这时候COMSOL会自动将边界条件转化为exp(1i*k·r)的相位因子。有次误把实部虚部分开计算，结果得到了类似分形艺术的诡异能带图——原来软件早把复数运算封装在底层，我们只需要关心物理意义上的传播常数。

当第一个完整的能带图跳出屏幕时，发现TE模在0.32c/a处出现了明显的带隙。但随后发现TM模的带隙位置总比文献值偏移5%，检查三天才发现是介电常数张量的各向异性没考虑——在弱形式中需要把epsilon_r改写为矩阵形式：

epsilon = [epsilon_xx 0 0; 0 epsilon_yy 0; 0 0 epsilon_zz];

这种张量操作在传统设置界面里需要层层菜单切换，而弱形式允许直接用矩阵语法描述材料特性，就像用不同的偏振光给每个晶格方向"上色"。

最后在结果分析时用MATLAB脚本批量提取特征频率：

data = mphgetmatrix(model,'omega'); band_gap = (min(data(2,:)) - max(data(1,:)))/mean(data(1,:));

这个简单的差值计算却能揭示光子晶体最核心的禁带特性。当看到18%的相对带隙宽度时，突然明白那些深夜的调试其实是在和麦克斯韦方程组玩一场三维的捉迷藏游戏，而弱形式就是那盏照亮隐身处的探照灯。