★ FEATURED ARTICLE 简介利用COMSOL Multiphysics与MATLAB联合仿真计算二维光子晶体带隙的完整脚本包特别适合光学仿真研究者、研究生以及从事光通信、光学微腔器件设计的工程人员快速上手参考。包内共4个文件、约144KB包含可运行的MATLAB脚本、可交互查看的模型图形文件、带隙与能带结构的结果对比图以及详细说明文档完整覆盖从几何建模、边界条件设置、频率扫描到后处理分析的全流程。光子带隙决定特定频段光的传播行为该脚本包能帮助读者直观理解这一物理现象。资源目前已有1065人学习下载可供同类研究借鉴参考。借助脚本包读者可复现二维光子晶体的带隙计算掌握COMSOL-MATLAB接口调用方式学会提取能带边界并分析带隙位置与宽度还能结合PML等边界条件优化仿真参数为实际器件设计提供可靠的仿真基础。1. 手动点GUI翻车之后我改用COMSOL-MATLAB脚本计算2D光子晶体带隙手动在COMSOL图形界面里算2D光子晶体带隙我一开始是拒绝脚本的建几何、设Floquet周期边界、沿不可约布里渊区扫k路径拖拖拽拽也能出色散图。直到要把介质柱半径从0.1a扫到0.4a、共8组参数做带隙优化时我彻底绷不住了——每组参数都要重复改半径、重新网格化、重建求解任务再把频率数据拷进Excel画带隙图一晚只够跑两组。换到COMSOL-MATLAB脚本后我用循环把几何更新、k路径采样、特征值求解与归一化后处理全部串起来8组参数半小时跑完顺手还能按模式交叉自动分类。这份资源解决的就是这个痛点适合正在赶仿真进度的研究生和工程人员把单次仿真任务升级成可重复、可批量、可改结构的参数化流程。2. 带隙计算原理与选型为什么是COMSOL-MATLAB而不是纯GUI或纯MATLAB2.1 光子晶体带隙的数学本质从Maxwell方程到Bloch特征值光子晶体的带隙本质上是周期性介电常数结构中电磁波色散关系的禁区。无源无损耗介质中Maxwell方程可以化为电场或磁场矢量本征方程以磁场为例∇ × (1/ε(r) ∇ × H) (ω/c)² H因为ε(r)具有晶格周期性Bloch定理把解写成H(r) e^(i k·r) u_k(r)其中k是第一布里渊区内的波矢u_k是周期函数。这样求解区域被限制在单个原胞内对应边界条件用Floquet周期边界表达。二维情况下如果介质柱沿z方向无限长且几何不随z变化TM和TE模式解耦通常用面外分量Ez或Hz标量方程来算计算量远小于三维全矢量模型。带隙就是指在某一归一化频率范围内对第一布里渊区边界上的所有k都不存在特征值。实际计算中不可能扫整个布里渊区面只会沿不可约布里渊区边界选点扫描得到ω(k)后看相邻能带之间有没有间隔。三角格子取Γ-M-K-Γ正方格子取Γ-X-M-Γ采样点足够密时带隙边缘才能看准。频率轴我习惯用a/λ表示也就是f·a/c0这样不同晶格常数的结构可以做横向对比这也是文献里最常见的做法。2.2 三种仿真路线对比GUI、纯MATLAB、COMSOL-MATLAB脚本COMSOL图形界面适合单次探索。快速看一眼某个结构有没有带隙手动搭模型并不慢但只要是批量改半径、批量改介电常数或配合优化函数做逆向设计GUI的重复劳动和维护成本都会直线上升。同一个模型有几十个参数组合时误改某一个参数而没注意整组结果就废了。纯MATLAB可以用平面波展开法或FDTD实现带隙计算。平面波展开在均匀网格上做得很漂亮但碰到高介电对比、异形孔、多层衬底这类几何就需要逐个开发算法数值稳定性和网格收敛都要自己把关。对于只懂结构、不想重写求解器的从业者这条路的学习成本明显高于直接使用成熟有限元工具的参数化接口。COMSOL-MATLAB脚本组合恰好把效率落差补上了。LiveLink for MATLAB把COMSOL模型对象暴露给MATLAB几何构建、材料赋值、网格划分、特征值求解、数据提取全部可以在脚本里完成同时还能在COMSOL桌面端打开同一个模型检查场分布。跟纯COMSOL辅助扫描相比它最大的优势是可以在MATLAB里做模式追踪和自动排序而且能直接接入全局优化循环这是图形界面批处理模式很难替代的。我见过不少课题组最后都停在“GUI建模型导出数据外部画图”的半自动状态瓶颈就在模式排序这一步手动处理几十个k点还能忍几百个点就完全没有可复现性了。2.3 特征值求解与频域研究的参数设置基准在COMSOL中做光子晶体带隙一般用“电磁波频域”接口研究步骤选择“特征值”。“所需特征值数量”是最关键的参数设置太少会漏掉高频分支设置太多会引入数值模并拖慢求解。我的经验是先按12起步某个k点不收敛再加到16或20。无损耗介质中特征频率为实数但脚本里经常拿到复数小虚部取实部即可。另一个容易被忽略的选项是“搜索基准点频率”。默认搜索在0附近对基模没问题但高频带需要把搜索基准设到目标频段否则求解器会一直找到最低频的若干个模而漏掉目标区间。特征值问题类型方面COMSOL默认处理线性特征值但含频散材料时应切换为二次特征值。容差我一般给1e-6低于1e-7在细网格上会有收敛假象高于1e-4则带隙边缘误差明显。这些设置直接影响脚本计算结果后面第五章会展开讲对应的失败现场。简单说带隙计算不是把模型跑完就结束特征值数量、搜索基准、容差三个参数在批处理时必须固化成脚本里的显式赋值否则半径扫描时换了参数结果是否可信很难判断。3. 脚本包拆解几何构建、材料参数与网格控制该怎么传3.1 脚本结构从mphopen到参数更新的执行顺序拿到这套脚本资源后第一件事是看它的目录组织。常见做法是把预建的COMSOL模板文件与四个MATLAB脚本分开模板里只有几何和物理框架所有数值参数都留成全局参数脚本负责注入参数、跑研究、取数据。这样做的好处是几何修改不必每次都在MATLAB里重建一遍只需要改模板里对应特征即可。% main_bandgap.m % 主控脚本其他三个子函数分别负责注入参数、运行、数据提取 clc; clear; close all; a 500e-9; % 晶格常数米 r 0.25 * a; % 介质柱半径米 eps_r 11.56; % 介质柱相对介电常数硅在红外波段附近常用值 Nmode 12; % 每个k点计算的基模数量 model mphopen(fullfile(template,pc2d_template.mph)); model.param.set(a, a); model.param.set(r, r); model.param.set(eps_r, eps_r); model.param.set(Nmode, Nmode);逻辑说明mphopen负责把COMSOL模板加载成模型对象param.set把MATLAB变量同步到COMSOL全局参数表。此后所有几何和物理设置都引用a、r、eps_r三个参数名脚本不需要逐个改界面。参数说明a、r的单位直接以米写入COMSOL在几何和物理里会自动保持一致eps_r是无量纲参数Nmode不是COMSOL原生参数名把它设进参数表只是方便后续study里引用。模型模板里研究步骤的“所需特征值数”应该绑定到Nmode否则这个参数不会生效。3.2 几何与材料参数晶格常数、半径、介电常数的传递方式如果模板里几何没有预建脚本也可以用COMSOL的Java API直接创建。以正方格子单原胞为例核心是画一个空气或介质背景方块再在中心放一个半径可变的圆柱model.component(comp1).geom(geom1).create(blk1,Block); model.component(comp1).geom(geom1).feature(blk1).set(size, {a,a,1e-9}); model.component(comp1).geom(geom1).create(cyl1,Cylinder); model.component(comp1).geom(geom1).feature(cyl1).set(r, r); model.component(comp1).geom(geom1).feature(cyl1).set(h, 1e-9); model.component(comp1).geom(geom1).run();逻辑说明Block的z方向尺寸设成1e-9并非物理厚度而是2D组件对面外维度的一种处理方式圆柱同样用极薄高度构造配合二维求解域使用。run()方法执行几何操作序列并更新网格框架。参数说明边长必须等于晶格常数a圆柱在中心处的相对位置由几何默认决定如果做的是空气孔型光子晶体只需把圆柱区域从材料域改成空气孔材料赋予逻辑反过来即可。材料赋值在MATLAB里同样走参数表而不是逐个界面选区域model.component(comp1).material(mat1).propertyGroup(def).set(epsilon_r, eps_r);这里的eps_r可以是标量也可以写成3x3矩阵。各向同性介质用标量最简单脚本包里默认即这种方式。如果要做液晶等各项异性材料改成对角矩阵并保留其他分量为1即可。我在实际中常被提醒的一点是不要忘了空气区域也有材料域COMSOL默认会补上一个“材料1”给所有域如果不显式把介质柱区域覆盖掉会出现所有域都是同一材料的情况带隙自然为零。3.3 网格与求解器配置控制数值误差的几个关键选项带隙数值计算对网格质量敏感尤其是高介电对比情形介质柱界面处的场变化剧烈太粗的网格会抬高甚至完全抹掉带隙。我一般先固定最大单元尺寸为a/4查看带隙上沿再加密到a/6和a/8对比如果上沿频率变化不超过0.5%就认为网格收敛。model.component(comp1).mesh(mesh1).feature(size).set(hmax, a/5); model.component(comp1).mesh(mesh1).run();逻辑说明mesh.feature(size).set直接覆写全局最大网格尺寸在参数扫描循环里第一次设置后后续几何变化较小COMSOL可以重用网格节约大量时间。参数说明hmax的取值决定网格密度a/5适合快速预览a/8适合最终计算更细的网格在2D模型里还撑得住3D全矢量模型会明显吃力。使用默认自由三角形网格时还可以在介质柱边界加边界层选项来降低界面附近的数值误差脚本模板里通常默认关闭需要时手动打开。求解器设置中最重要的是特征值数量和搜索基准这两个参数必须绑定到全局参数而不是写死在研究步骤里否则批处理扫描半径时无法自动调整。代码上可以这样控制model.study(std1).feature(eig).set(eigmethod, auto); model.study(std1).feature(eig).set(eigcount, Nmode);逻辑说明这两行直接修改研究步骤中的“特征值”节点。eigmethod设为auto允许COMSOL自动选择线性或非线性求解器eigcount与全局参数Nmode绑定保证每个k点都请求同样数量的本征值。参数说明如果是频散材料建议把eigmethod切到nonlinear否则自动模式可能走线性求解器导致结果错误。实际执行中求解器日志会打印“本征值数量不足”这类警告一旦出现优先增加eigcount而不是加密网格因为漏模问题通常不是网格精度造成的。4. 跑通第一个带隙计算从LiveLink连接到k路径扫描与数据导出4.1 环境配置启动LiveLink并确认mph文件可加载在MATLAB里运行任何COMSOL脚本前先确认本机装了对应的LiveLink for MATLAB模块。启动方式有两种第一种是直接在MATLAB命令窗口调用mphstart这是脚本包最常采用的方式mphstart(2036); import com.comsol.model.* import com.comsol.model.util.*逻辑说明mphstart会在指定端口创建COMSOL服务器进程后面import三行是加载COMSOL Java API的固定写法缺少任何一行ModelUtil.create或mphopen都会报“找不到类”。端口2036是默认值冲突时可以换成其他空闲端口。参数说明本机只装一个COMSOL时用默认端口即可如果电脑上同时跑多个服务或与远程服务器连接端口号必须与启动参数一致。Windows下还要检查COMSOL安装时是否勾选“LiveLink for MATLAB”这个选项不会默认安装缺了会直接找不到模块。加载模板模型用mphopenmodel mphopen(fullfile(pwd,template,pc2d_template.mph));逻辑说明fullfile把当前路径和相对子目录拼成绝对路径避免MATLAB工作目录变化导致找不到文件。加载成功后model就是一个COMSOL模型对象后续所有操作都通过model调用。参数说明mph文件版本与COMSOL版本不一致时open过程可能升级模型或报版本警告一般不影响参数扫描但建议保留一份原始模板作为后悔药。我自己的习惯是复制一份模板到工作目录再打开绝不直接改原始模板否则改坏了连恢复的余地都没有。4.2 定义k路径倒格矢坐标下的Γ-X-M-Γ扫描正方格子的不可约布里渊区取Γ-X-M-Γ高对称点坐标需要换算成实际波数。以晶格常数a计算Γ在(0,0)X在(π/a,0)M在(π/a,π/a)。采样时对每一段分别做等距插值Nseg 30; Gamma [0, 0]; X [pi/a, 0]; M [pi/a, pi/a]; path1 Gamma (X - Gamma) .* linspace(0,1,Nseg).; path2 X (M - X) .* linspace(0,1,Nseg).; path3 M (Gamma - M) .* linspace(0,1,Nseg).; kpath [path1; path2; path3];逻辑说明每段用30个点整个路径共90个k点跑一次特征值扫描需要90次求解。如果只做粗略预览可以把Nseg降到15但要判断带隙是否存在建议保持30以上。参数说明linspace(0,1,Nseg)生成均匀采样系数乘到向量上得到每段的线性插值点。kpath中每一行的两个数值会分别传给全局参数kx和ky用来更新Floquet周期边界条件中的Bloch波矢。提示在COMSOL的Floquet条件里k分量必须以波数单位1/m输入。如果模板中的边界条件是相位型参数脚本包需要在参数表中额外乘以a再赋值。把k点写入模型并执行求解for i 1:size(kpath,1) model.param.set(kx, kpath(i,1)); model.param.set(ky, kpath(i,2)); model.study(std1).run(); end逻辑说明循环里参数变化只影响Floquet边界条件几何、网格、材料都不变所以COMSOL会复用网格和大部分内部分解结果单点求解比重新建模快很多。参数说明这个循环没有显式调用model.sol(sol1).runAll()因为study.run本身就包含求解步骤。如果脚本包内研究设置了独立的解序列则要改成先运行解序列再更新数据集否则后处理取不到数据。三角格子的处理方式同理只是把高对称点换成Γ-K-M-Γ倒格矢基矢多一个120度夹角。脚本包里单独有一个triangular_bandgap.m用的就是那套基矢。第一次跑的时候不要混用正方格子的kx/ky直接塞进三角格子模板色散关系一定会错。4.3 特征频率提取与色散曲线绘制特征值求解完成后把每个k点的本征频率取出来按升序排列再归一化。无损耗情况直接取实部完整写法如下freq_all zeros(Nmode, size(kpath,1)); for i 1:size(kpath,1) model.param.set(kx, kpath(i,1)); model.param.set(ky, kpath(i,2)); model.study(std1).run(); f mphglobal(model,freq); freq_all(:,i) sort(real(f(:))); end % 横轴为路径累计距离 k_dist [0; cumsum(sqrt(sum(diff(kpath).^2,2)))]; % 归一化频率 a/lambda f*a/c0 plot(k_dist, freq_all * a / 3e8, o);逻辑说明mphglobal取全局变量freq返回的结构体里直接包含特征频率取实部并sort是必须的因为求解器返回的频率顺序不保证从低到高。k_dist用相邻k点笛卡尔距离累加得到这样横轴才和色散关系图谱一致。参数说明归一化频率用f*a/c0表示c0取3e8 m/s如果介质参数里用的是角频率ω还要除以2π本脚本包里存储的一律是普通频率f直接乘a除以光速即可。绘制时纵轴的单位是a/λ这是光子晶体文献中通用的无量纲表达。画图时还有一个习惯把前两条能带用不同颜色标出来带隙边缘就落在第一条能带的上沿和第二条能带的下沿之间。如果只看单色点图很难分清哪条是第一条、哪条是第二条。脚本包里会额外输出两个变量gap_lower和gap_upper就是这两个边界频率。5. 带隙计算常见问题排查从特征值漏检到模式交叉的五个坑前四章把标准流程走通了但真正上手后大概率会遇到各种“带隙图不对”的状况。以下五个坑是我跑这套脚本时实际踩过的逐条按现象-原因-解决列出来。5.1 特征值数量不够导致色散曲线断裂现象色散曲线在中段突然少了几条分支或者某些k点只有8条频率而其他点有12条带隙位置肉眼可见地对不齐。原因固定Nmode12时个别k点因为简并或高密度模式收敛变慢求解器只返回了部分特征值sort之后矩阵缺行。解决先加大Nmode到20重跑一次如果曲线补全说明是漏检而不是物理带隙同时把求解器容差从1e-4改为1e-6并检查日志中是否有“本征值未收敛”警告。另一种常见做法是给研究步骤设置“最少特征值数量”让它在该k点主动搜索更多特征值。5.2 伪模混入检查面外厚度与极化设置现象色散曲线底部出现一条几乎为0的平带或者整个频谱多出一组与结构无关的超低频率模式。原因2D组件中物理场没有锁定为面外分量求解器把面内和面外偏振混在一起也可能是几何中圆柱的h设置得太大导致面外方向产生非物理谐振。解决在“电磁波频域”接口里把求解变量改成“面外电场”或“面外磁场”对应TM/TE中的一种不要使用默认的全矢量公式。h统一设成1e-9量级远小于晶格常数a。5.3 Floquet边界k矢量方向与倒格矢对不上现象Γ点频率不是最低色散曲线整体左移或右移甚至出现不正常的镜像对称。原因Floquet周期边界条件里输入的不是波矢而是相位偏移参数。相位θ与波矢k的关系是θk·a直接把k值填到相位输入框等于乘了一个a因子有时还差2π。解决先确认模板边界条件类型。如果是“波矢”类型单位是1/mkx和ky直接传波数即可如果是“相位”类型要把kxa和kya作为输入。脚本包里两种写法都预留了参数注释里会标明当前模板属于哪种。5.4 批量扫描内存溢出或网格反复重画现象k路径循环跑到第10个点后MATLAB开始卡顿再过几个点直接Out of Memory不得不重启软件。原因每次model.study(std1).run()都会保留求解数据和解集90次扫描累积起来占掉大量内存另外如果每次循环都重新执行网格划分临时网格数据也会堆积。解决在循环开始前设置保存策略为“只保存特征值不保存场”循环内每个k点求解完成后用model.result().clear()清理临时结果。网格只划分一次几何和材料参数变化不需要重新划分确认网格节点没有出现在循环体内即可。5.5 模式交叉导致带隙边缘误判现象带隙上下沿频率在相邻k点之间跳动画出的曲线在某处突然从低频跳到高频带隙宽度看起来忽大忽小。原因不同能带在布里渊区边界附近会发生交叉单纯按频率升序sort会把两条不同模式的曲线连成一条造成伪带隙或伪闭合。解决在MATLAB里增加简单的模式追踪。先保存上一个k点的特征向量在当前k点的多个特征向量中用场分布交叠积分找与上一点重叠度最高的模式作为同一支的延续。代码层面上可以简化成先按频率排序再做相邻点最近邻匹配但对高简并度结构仍然推荐做交叠积分。还有一个取巧办法把TE和TM模式分别计算、分别绘制避免偏振混叠造成的模式交叉假象。6. 进阶把带隙计算做成参数扫描工具并做收敛性验证当单条色散曲线能稳定复现以后真正的生产力是把k路径扫描封装成函数再对几何参数做批处理。我的做法是把第四章的循环整理成一个函数function [freq_norm, k_dist] run_bandgap(model, r, a, Nseg, Nmode) model.param.set(r, r); % ... 构建kpath并循环求解 ... end然后在外层循环半径参数输出带隙宽度随r/a的变化r_ratio 0.15:0.05:0.4; gap_width zeros(size(r_ratio)); for j 1:length(r_ratio) [freq_norm, k_dist] run_bandgap(model, r_ratio(j)*a, a, 40, 16); gap_width(j) compute_gap(freq_norm, k_dist); end plot(r_ratio, gap_width/a, o-);逻辑说明这样的封装让脚本包变成一个小工具换材料、换晶格常数、换结构类型都只需要改参数。需要提醒的是每次改r/a后都要验一下最高能带是否收敛因为半径变化会改变模式密度。参数说明Nseg从30提到40是为了让带隙边缘更精确网格也应同时从a/5提升到a/6以上。compute_gap是脚本包内置的后处理函数它先剔除被模式追踪判定为伪模的频率再取第一条能带下沿与第二条能带上沿的差。网格收敛性验证我一般固定在一个带隙最宽的半径上做对同一结构分别用a/4、a/5、a/6、a/8四套网格跑一次记录带隙上沿频率。只要上沿频率在相邻网格间的相对变化小于0.5%就认为该网格密度可以用于当前参数扫描。工程上这种做法比盲目追求细网格更高效也能在投稿前给审稿人一个可重复性说明。从那以后我每次换一个新的光子晶体结构都会先花十几分钟做完网格收敛性检查再跑正式参数扫描避免把数值假带隙当物理发现输出。希望这篇笔记能把你在COMSOL-MATLAB脚本带隙计算这条路上遇到的最常见障碍扫掉一部分帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
