简介「photonic-bandgap-comsol-matlab」面向光学、光子晶体方向的研究生与工程人员聚焦如何借助 COMSOL Multiphysics 的光波光学模块配合 MATLAB 脚本完成二维光子晶体能带结构与光子带隙的自动化计算。压缩包共 4 个文件、约 144KB含 MATLAB 主脚本、fig 图形文件、结果对比图与说明文档主脚本承担建模、参数扫描与结果读取其余文件分别用于保存能带曲线、核对仿真输出与解释使用方式。目前已有 1066 人学习下载。读者可从中了解 COMSOL 与 MATLAB 联合调用的接口写法参考晶格常数、单元形状、PML 边界与频率范围的组织方式并借助结果对比图核对带隙位置与宽度为光通信、光学微腔等场景的器件设计积累可复用的脚本模板与调试思路。1. 带隙算出来是负的这套 COMSOL-MATLAB 脚本把哪几步固定住了二维光子晶体的带隙算出来是负的或者相邻两条能带在某个 k 点莫名其妙交叉多半不是物理错了而是周期边界、波矢扫描和模式排序这三步里有一环没固定住。这份 photonic-bandgap-comsol-matlab 脚本包干的正是这件事用 MATLAB 通过 LiveLink 驱动 COMSOL把二维光子晶体的元胞建模、Floquet 周期边界、布里渊区高对称路径扫描、特征频率求解和带隙判定串成一条可重复执行的流水线。它适合正在做光子晶体能带、微腔或波导设计又不想每次都在 GUI 里手点参数的人。改几个参数就能跑出 Γ-X-M-Γ 路径上的 band structure省下的是反复调边界条件的时间。2. 技术栈与目录结构MATLAB 当控制层COMSOL 当计算引擎2.1 LiveLink 的分工为什么不在 GUI 里点COMSOL 本身能算光子晶体能带GUI 里也能设周期边界和参数扫描。问题出在“重复”和“记录”上每换一个晶格常数、每换一个孔半径都要重新点一遍几何、材料、边界、研究设置点错一个边编号结果就偏了。LiveLink for MATLAB 把 COMSOL 变成一个可以被脚本调用的计算引擎MATLAB 负责参数、循环、后处理和画图COMSOL 只负责求解。分工清楚之后整个计算过程就是一段可版本管理的代码而不是一堆散落在模型树里的设置。这个项目的技术栈不复杂但每一层都有明确的职责层级工具职责控制层MATLAB参数管理、k 点循环、结果收集、带隙判定、画图求解层COMSOL Multiphysics二维频域电磁波、Floquet 周期边界、特征频率求解接口层LiveLink for MATLAB把 COMSOL 模型对象暴露成 MATLAB 可调用的 API数据层CSV / MAT保存每个 k 点的特征频率和带隙表常见做法是把所有可调参数收进一个params.m主脚本只负责按顺序调用建模、求解、后处理三个函数。这样二次开发时不用翻主脚本改参数文件就行。2.2 目录结构与每个文件的职责源码包的目录组织大概是这个形状不同版本可能略有差异但职责划分基本一致photonic-bandgap-comsol-matlab/ ├── main_bandgap.m # 主入口加载参数、建模、扫描、出图 ├── config/ │ └── params.m # 所有可调参数集中在这里 ├── src/ │ ├── build_unit_cell.m # 元胞几何与材料 │ ├── set_physics.m # 频域电磁波与 Floquet 周期边界 │ ├── set_mesh.m # 网格划分 │ ├── run_k_sweep.m # k 点循环与特征频率提取 │ └── find_bandgaps.m # 带隙判定 ├── models/ │ └── pc_2d_base.mph # 可选的空白模型模板 ├── results/ │ ├── band_structure.csv # 归一化频率矩阵 │ └── figures/ # 输出的能带图 └── utils/ ├── comsol_startup.m # 启动 COMSOL Server 并连接 └── high_symmetry_path.m # 生成 Γ-X-M-Γ 路径这个结构的好处是每个文件只做一件事。build_unit_cell.m只管几何set_physics.m只管物理场和边界run_k_sweep.m只管循环和取数。排查问题时可以单独把某一段拎出来跑不用每次都从头执行整个主脚本。我一般会把results/整个目录加到.gitignore只保留params.m和src/下的脚本这样换机器时不会把几 GB 的中间结果带过去。2.3 环境准备启动顺序和版本匹配COMSOL 和 MATLAB 的联动有一个固定的启动顺序顺序错了就会卡在连接那一步。常见做法是先在系统里启动 COMSOL Multiphysics Server记下端口号默认 2036再在 MATLAB 里调用mphstart连接。也可以在 MATLAB 里直接启动 server但需要把 COMSOL 的bin目录加到系统路径里。% utils/comsol_startup.m function model comsol_startup(port) % 连接已经启动的 COMSOL Server返回一个空的模型容器 if nargin 1 port 2036; % COMSOL Server 默认端口 end import com.comsol.model.util.* try mphstart(port); catch ME error(连接 COMSOL Server 失败先确认 server 已启动、端口 %d 没被占用%s, ... port, ME.message); end ModelUtil.showProgress(true); % 在 MATLAB 命令行打印求解进度 model ModelUtil.create(Model); end这段代码的逻辑很直接mphstart负责建立 MATLAB 到 COMSOL Server 的 TCP 连接ModelUtil.create创建一个空的模型对象后面的几何、物理场、研究都挂在这个对象上。ModelUtil.showProgress(true)不是必须的但求解特征频率时能看到进度条长时间扫描不容易以为程序卡死。参数说明port要和 COMSOL Server 启动时显示的端口一致一台机器上同时跑多个 COMSOL 实例时端口会顺延不要硬编码。如果mphstart报连接超时先检查 server 是否真的在运行再检查防火墙有没有拦本地回环端口。这一步翻车的人不少血泪经验是先把 server 启动日志看一眼比在 MATLAB 里反复试端口快得多。版本匹配也值得提一句COMSOL 的 LiveLink 接口和 MATLAB 版本之间有兼容矩阵太新的 MATLAB 配太旧的 COMSOL 可能在ModelUtil那一层就报方法找不到。拿到源码包之后先确认本机 COMSOL 版本支持当前 MATLAB再往下走。3. 物理建模与核心脚本从元胞到 Bloch 边界的完整链路3.1 元胞、晶格常数与不可约布里渊区二维光子晶体的能带计算第一步是把无限周期结构截成一个元胞。正方晶格的元胞取边长等于晶格常数 $a$ 的正方形中心挖一个圆孔或者放一根介质柱孔半径或柱半径记作 $r$。介电常数在元胞内部分布背景介质一个值散射体一个值。这个元胞在二维平面内按 $a$ 平移复制就得到整个晶体所以求解域只需要一个元胞边界上用周期条件把相邻元胞“接”起来。晶格常数 $a$ 决定了归一化频率的尺度。带隙位置通常用 $a/\lambda$ 或 $\omega a/2\pi c$ 表示这样换材料、换晶格常数时能带图可以直接对比。正方晶格的高对称点是 Γ(0,0)、X(0.5,0)、M(0.5,0.5)单位是 $2\pi/a$。不可约布里渊区边界就是 Γ-X-M-Γ 这条折线扫描这条线上的 k 点就能覆盖带隙判定需要的信息。% utils/high_symmetry_path.m function [kx, ky, s] high_symmetry_path(nk) % 生成正方晶格不可约布里渊区边界 Γ-X-M-Γ 上的 k 点 % 输入 nk每条线段上的点数含端点 % 输出 kx, ky归一化波矢单位 2*pi/a % 输出 s路径序号画图时当横坐标用 kG [0.0, 0.0]; kX [0.5, 0.0]; kM [0.5, 0.5]; seg1 [linspace(kG(1), kX(1), nk), linspace(kG(2), kX(2), nk)]; seg2 [linspace(kX(1), kM(1), nk), linspace(kX(2), kM(2), nk)]; seg3 [linspace(kM(1), kG(1), nk), linspace(kM(2), kG(2), nk)]; % 去掉相邻线段重复的端点 kpts [seg1; seg2(2:end, :); seg3(2:end, :)]; kx kpts(:, 1); ky kpts(:, 2); s (0:size(kpts, 1) - 1); end这段函数的参数只有nk它决定每条高对称线上的采样密度。nk 21是一个常见起点三条线段去重之后一共 61 个 k 点。点数太少会漏掉能带交叉附近的细节带隙边界会算宽点数太多单次扫描时间线性增长。实际做法是先跑nk 11看整体形状确认带隙大致位置后再用nk 31或51做精细扫描。s是给画图用的横坐标能带图的横轴不是真实的 k 值而是路径上的累计序号在 Γ、X、M 三个位置打刻度。3.2 Floquet 周期边界怎么落到脚本里周期边界是这个项目里最容易设错的一环。COMSOL 的Periodic边界条件有两种类型连续型和 Floquet 型。光子晶体能带计算必须用 Floquet 型因为它允许场在跨越元胞时带一个相位因子 $e^{i\mathbf{k}\cdot\mathbf{a}}$这个相位因子就是 Bloch 波矢。用连续型周期边界等于强制 $k0$扫出来的不是能带是 Γ 点的几个驻模。脚本里设置 Floquet 周期的关键是把波矢分量写成全局参数左右一对边界给 $k_x$上下对给 $k_y$% src/set_physics.m节选 function model set_physics(model, p) % 二维频域电磁波特征频率求解 model.physics.create(ewfd, ElectromagneticWavesFrequencyDomain, geom1); % 背景介质 model.physics(ewfd).material(mat1).set(relpermittivity, {num2str(p.eps_host)}); % 散射体介质域编号按实际几何调整 model.physics(ewfd).material(mat2).selection.set(p.hole_domains); model.physics(ewfd).material(mat2).set(relpermittivity, {num2str(p.eps_hole)}); % Floquet 周期边界左右一对 model.physics(ewfd).create(pcx, Periodic, 1); model.physics(ewfd).feature(pcx).selection.set(p.left_right_edges); model.physics(ewfd).feature(pcx).set(PeriodicType, Floquet); model.physics(ewfd).feature(pcx).set(kFloquet, {kx, 0, 0}); % Floquet 周期边界上下对 model.physics(ewfd).create(pcy, Periodic, 1); model.physics(ewfd).feature(pcy).selection.set(p.bottom_top_edges); model.physics(ewfd).feature(pcy).set(PeriodicType, Floquet); model.physics(ewfd).feature(pcy).set(kFloquet, {0, ky, 0}); end这里的kx和ky不是数值是模型里的参数名实际值在扫描循环里用model.param.set逐点写入。kFloquet的三个分量对应三维二维模型里第三个分量填 0。p.left_right_edges和p.bottom_top_edges是边的编号不建议硬编码成[1 3]这种字面量因为几何一变编号就漂。稳妥的做法是在params.m里按坐标筛选边或者建模之后用mphgetadj之类的接口按位置取边。硬编码边编号是新手最容易踩的坑之一换个 COMSOL 版本或者改一下几何顺序周期边界就套到错误的边上求解不报错但结果是错的。还有一个细节Floquet 周期边界的相位因子方向要和波矢方向一致。左右边界的法向是 $x$所以给 $k_x$上下边界的法向是 $y$所以给 $k_y$。如果两边都给同一个分量相当于在某个方向上强制 $k0$带隙会算得偏高或偏低。3.3 特征频率求解器与网格设置COMSOL 里算能带用的是特征频率研究Eigenfrequency不是频域扫描。频域扫描给定一个频率求场分布特征频率求解给定边界条件求系统的本征频率后者才对应能带上的一个个模式。研究设置里要指定求解多少个特征频率neigs以及一个频移shift。频移的作用是把求解器引到感兴趣的频率附近避免它去挖零模或者极低频的伪模。% 研究设置特征频率围绕目标频段做频移 model.study.create(std1); model.study(std1).create(eig, Eigenfrequency); model.study(std1).feature(eig).set(neigs, p.nmodes); model.study(std1).feature(eig).set(shift, p.fshift); model.study(std1).feature(eig).set(eigwhich, lm); % 最接近 shift 的模参数说明p.nmodes是每个 k 点求解的模式数至少要比你关心的能带数多两三个留出余量给伪模p.fshift是频移单位 Hz通常设在归一化频率 $a/\lambda \approx 0.3$ 对应的实际频率附近eigwhich lm表示找最接近频移的特征值比默认的lr最大实部更适合带隙计算因为带隙附近的模式不一定实部最大。网格设置直接决定频率算得准不准。光子晶体元胞里有圆孔孔边缘是曲面网格太粗会让特征频率整体偏移。常见做法是给孔边界加一层加密同时限制全局最大单元尺寸% src/set_mesh.m function model set_mesh(model, p) model.mesh.create(mesh1, geom1); model.mesh(mesh1).autoMeshSize(4); % 整体“较细” model.mesh(mesh1).feature(size).set(hmax, p.hmax); model.mesh(mesh1).feature(size).set(hmin, p.hmax / 10); % 孔边界加密 model.mesh(mesh1).create(edg1, Size); model.mesh(mesh1).feature(edg1).selection.set(p.hole_edges); model.mesh(mesh1).feature(edg1).set(hmax, p.hmax / 3); model.mesh(mesh1).run; endp.hmax一般取 $a/20$ 到 $a/30$。这个值不是拍脑袋定的后面第 4 章会做网格收敛测试把hmax从 $a/10$ 逐步降到 $a/40$看带隙中心频率偏移是否小于 1%。如果偏移还在漂说明网格还不够细。孔边界的局部加密比全局加密划算因为场在孔边缘变化最剧烈元胞内部相对平缓。3.4 参数入口与几何构建把参数集中到一个函数里是这个脚本包最值得抄的做法。光子晶体带隙对 $r/a$ 和介电常数对比很敏感二次开发时经常要扫这两个量。参数散落在各个脚本里扫到第三组就忘了第一组改过什么。集中之后主脚本只需要一行p params();所有函数共享同一个结构体。% config/params.m function p params() % 二维光子晶体带隙计算的参数入口 p.a 1e-6; % 晶格常数单位 m p.r 0.30 * p.a; % 孔半径0.3a 是常见起点 p.eps_host 12.0; % 背景介质相对介电常数 p.eps_hole 1.0; % 孔内相对介电常数 p.pol TE; % 偏振约定TE 表示电场沿面外 p.nk 21; % 每条高对称线上的点数 p.nmodes 8; % 每个 k 点求解的模式数 p.fshift 0.3 * 3e8 / p.a; % 频移对应 a/lambda 约 0.3 p.hmax p.a / 20; % 全局最大网格尺寸 p.hole_edges []; % 孔边界编号建模后按坐标填 p.left_right_edges []; p.bottom_top_edges []; p.outdir fullfile(pwd, results); end几何构建的函数只做三件事建正方形元胞、建圆孔、做布尔差集。圆孔的位置放在元胞中心这样周期平移之后孔阵列是对齐的。如果孔放在角上平移之后会变成四个四分之一孔拼在一起几何上等价但边界条件容易设错。% src/build_unit_cell.m function model build_unit_cell(model, p) model.geom.create(geom1, 2); % 正方形元胞 model.geom(geom1).create(cell, Rectangle); model.geom(geom1).feature(cell).set(size, {a, a}); model.geom(geom1).feature(cell).set(pos, [0 0]); % 中心圆孔 model.geom(geom1).create(hole, Circle); model.geom(geom1).feature(hole).set(r, r); model.geom(geom1).feature(hole).set(pos, {a/2, a/2}); % 布尔差集元胞减去孔 model.geom(geom1).create(dif1, Difference); model.geom(geom1).feature(dif1).selection(input).set({cell}); model.geom(geom1).feature(dif1).selection(input2).set({hole}); model.geom(geom1).run; % 按坐标筛选边避免硬编码编号 p.left_right_edges mphselectbox(model, geom1, [0 0; 0 p.a], boundary); p.bottom_top_edges mphselectbox(model, geom1, [0 0; p.a 0], boundary); endmphselectbox是按坐标框选几何实体的接口比记边编号可靠。[0 0; 0 p.a]框住的是左边那条边[0 0; p.a 0]框住的是下边那条边。实际脚本里可能用mphgetadj或直接遍历边坐标逻辑是一样的位置定边不靠编号。3.5 k 点扫描与后处理判据扫描循环是整个脚本里最耗时的部分。每个 k 点都要重新设置kx、ky重新求解特征频率然后把频率取出来存进矩阵。循环里最容易出错的是单位换算high_symmetry_path返回的是归一化波矢单位是 $2\pi/a$写进 COMSOL 参数时要乘上 $2\pi/a$ 换成 $1/m$。% src/run_k_sweep.m核心循环 function freq run_k_sweep(model, p, kx, ky) npts numel(kx); freq zeros(p.nmodes, npts); for i 1:npts % 归一化波矢转实际波矢单位 1/m kx_real kx(i) * 2 * pi / p.a; ky_real ky(i) * 2 * pi / p.a; model.param.set(kx, [num2str(kx_real), [1/m]]); model.param.set(ky, [num2str(ky_real), [1/m]]); model.study(std1).run; % 取特征频率实部虚部对应损耗无损耗时接近 0 f mphglobal(model, ewfd.freq, dataset, dset1); f sort(real(f(:))); % 模式数不够时补 NaN避免矩阵维度不一致 if numel(f) p.nmodes f(end1:p.nmodes) NaN; end freq(:, i) f(1:p.nmodes); end % 归一化频率 a/lambda f * a / c freq freq * p.a / 3e8; end这段循环的逻辑外层遍历 k 点内层由 COMSOL 求解器完成。mphglobal从结果数据集里取全局标量ewfd.freq返回的是一个列向量。sort保证每个 k 点的模式按频率从低到高排后处理判带隙时才能按行对应能带。归一化那一步放在循环外面一次性把 Hz 转成 $a/\lambda$省得在循环里反复乘除。后处理判带隙的判据不复杂把频率矩阵按行看成一条条能带第 $n$ 条带在所有 k 点上的最大值如果小于第 $n1$ 条带在所有 k 点上的最小值两条带之间就有带隙。带隙宽度是两者之差中心频率是两者之和的一半。% src/find_bandgaps.m function gaps find_bandgaps(freq_norm, min_width) % freq_norm归一化频率矩阵行是模式列是 k 点 % min_width带隙最小宽度低于这个值不算 nb size(freq_norm, 1); gaps []; for n 1:nb-1 upper max(freq_norm(n, :)); % 第 n 条带的上边界 lower min(freq_norm(n1, :)); % 第 n1 条带的下边界 if lower upper width lower - upper; center 0.5 * (lower upper); if width min_width gaps [gaps; n, upper, lower, width, center]; end end end end这个判据的前提是能带已经按频率排好序且没有出现“带交叉”——两条能带在某个 k 点交换了顺序。正方晶格的光子晶体在低阶带里交叉不常见但高阶带或者参数接近带隙关闭时会出现。稳妥做法是在排序之后再做一次连续性跟踪每个 k 点选择与上一个 k 点频率最接近的模式作为同一条带而不是简单按频率排序。脚本包里如果只做了sort高阶带的结果要打个问号这一点在第 5 章会展开。4. 跑通与结果验证一条命令跑出 Γ-X-M-Γ 能带4.1 最短跑通路径拿到源码包之后最短路径是四步启动 COMSOL Server、在 MATLAB 里cd到项目根目录、改config/params.m里的a和r、运行main_bandgap。主脚本的职责就是把前面几节的函数串起来并在最后把结果写盘。% main_bandgap.m clear; clc; addpath(genpath(fullfile(pwd, src))); addpath(genpath(fullfile(pwd, utils))); addpath(genpath(fullfile(pwd, config))); p params(); % 1. 读参数 model comsol_startup(); % 2. 连 COMSOL Server model build_unit_cell(model, p); % 3. 建元胞 model set_physics(model, p); % 4. 物理场与周期边界 model set_mesh(model, p); % 5. 网格 [kx, ky, s] high_symmetry_path(p.nk);% 6. 高对称路径 freq run_k_sweep(model, p, kx, ky); % 7. 扫描求解 gaps find_bandgaps(freq, 0.01); % 8. 判带隙 writematrix(freq, fullfile(p.outdir, band_structure.csv)); disp(gaps);第一次跑建议把p.nk设成 11p.nmodes设成 6先确认整条链路能走通、结果形状对再加密。main_bandgap.m里没有画图代码画图单独放一个函数或者直接在命令行里plot(s, freq., .-)这样求解和展示解耦批量跑的时候不会因为画图拖慢速度。跑通之后results/band_structure.csv里是一个nmodes × npts的矩阵每列对应一个 k 点每行对应一条能带。gaps矩阵每一行是一条带隙列依次是能带序号、上边界、下边界、宽度、中心频率。4.2 band structure 图怎么读能带图的横轴是 Γ-X-M-Γ 路径上的累计序号纵轴是归一化频率 $a/\lambda$。每个 k 点上有若干个点连起来就是一条条能带。带隙出现在两条能带之间没有点的频率区间图上表现为一条横向的空白带。看这张图的时候先确认三件事Γ 点横轴起点的频率是否合理X 点和 M 点的频率是否有明显的分裂带隙上下边界是否平直。如果 Γ 点的频率接近 0大概率是伪模混进来了如果 X 点和 M 点的频率几乎一样可能是周期边界的方向设反了如果带隙边界随 k 点剧烈抖动说明 k 点太稀或者网格不收敛。TE 和 TM 偏振的带隙位置不一样这是正常现象。二维光子晶体里空气孔阵列的 TE 带隙通常比 TM 带隙更明显介质柱阵列反过来。脚本里p.pol这个参数最终要落到物理场的偏振设置上TE 对应电场沿面外二维模型里可以用面外电场分量TM 对应磁场沿面外。如果脚本里只设了一个偏振换偏振时要改物理场设置不能只改参数名。4.3 网格与 k 点收敛验证带隙中心频率算出来是 0.42 还是 0.40取决于网格和 k 点够不够密。验证方法不复杂但必须做固定其他参数只变hmax看带隙中心频率怎么变固定hmax只变nk看带隙宽度怎么变。网格尺寸带隙中心 $a/\lambda$带隙宽度相对上一次变化a/100.4180.031—a/150.4120.0361.4%a/200.4090.0380.7%a/300.4080.0390.2%a/400.4080.039 0.1%这张表是典型收敛曲线网格从 $a/10$ 加密到 $a/20$中心频率偏移接近 2%从 $a/20$ 到 $a/30$ 偏移降到 0.2% 以下再加密变化很小。实际项目里取 $a/20$ 到 $a/30$ 之间的值兼顾精度和速度。k 点收敛类似nk从 11 加到 21带隙宽度会收窄一点因为边界 k 点采得更密能带极值找得更准从 21 加到 41 变化通常在 1% 以内。提示收敛测试不要只做一次。换晶格类型、换填充比、换介电常数之后收敛行为会变。填充比越高、介电常数对比越大需要的网格越细。4.4 结果导出与复用band_structure.csv是最值得保存的中间结果。有了它画图、判带隙、对比不同参数都不用重新跑 COMSOL。二次开发时常见做法是把params.m里的关键参数和band_structure.csv一起打包文件名带上r_over_a和eps比如band_r0.30_eps12.csv这样扫参数扫出来的几十组结果不会混在一起。gaps矩阵可以单独存成bandgaps.csv每行一条带隙列名固定。后面要做带隙随 $r/a$ 变化的相图直接读这个文件就行不用再碰 COMSOL。导出的时候注意归一化频率的单位统一别一组用 $a/\lambda$、一组用 THz混在一起画图会得到一条斜率莫名其妙的曲线。5. 避坑与排查特征频率求解最常见的五类翻车5.1 求解阶段的三个高频问题现象一每个 k 点都解出一堆接近 0 的频率。原因通常是频移设得太低求解器把静模和伪模一起挖出来了。COMSOL 的特征频率求解器在频移附近找模频移接近 0 时零空间里的伪模就会被选中。解决办法是把p.fshift设到目标带隙附近比如归一化频率 0.3 对应的实际频率同时在取结果时过滤掉 $a/\lambda 0.05$ 的点。如果过滤之后模式数不够把p.nmodes加大用数量换有效模式。现象二扫描到某个 k 点突然报“找不到特征频率”。原因多半是频移落在了两条能带之间的带隙里那个频率附近确实没有本征模。带隙计算本来就会经过带隙区间求解器在带隙里找不到模是正常的但脚本如果没做容错就会中断。解决办法是在循环里包一层try-catch求解失败时给该 k 点填NaN继续下一个点。后处理判带隙时跳过含NaN的列或者把NaN当作“该模式不存在”处理。现象三同一组参数跑两次带隙宽度差了几个百分点。原因可能是网格自适应每次生成的网格不完全一致也可能是特征频率求解器对初值敏感。解决办法是关掉自适应网格重划分固定hmax和边界加密层数求解器那边把eigwhich固定成lm不要用默认设置。如果还漂把p.nmodes加大到能覆盖带隙上下各两条带减少模式漏选的概率。5.2 后处理与结果解读的坑现象四带隙表里出现负宽度或者相邻带隙重叠。这是模式排序错乱导致的。简单sort在每个 k 点独立排序遇到能带交叉时会把两条带的顺序交换后处理按行相减就得到负值。解决办法是做能带跟踪从第一个 k 点开始每个新模式选择与上一个 k 点同序号模式频率最接近的那个而不是无条件按频率排序。跟踪之后再做带隙判定负宽度基本消失。如果嫌跟踪麻烦至少在判带隙前检查相邻两行的频率差是否连续出现跳变就标记出来人工看一眼。现象五TE 和 TM 的带隙位置几乎一样。二维光子晶体的 TE 和 TM 带隙通常不在同一个频率两个偏振算出来一样大概率是物理场设置里偏振没真正分开。检查set_physics.m里有没有设置面外分量TE 需要电场只有面外分量TM 需要磁场只有面外分量。如果两个偏振用的是同一套边界和材料设置只是改了p.pol这个字符串求解结果不会变。偏振要落到物理场的分量约束或者端口设置上不能只改参数名。5.3 环境与接口的坑现象六mphstart报连接被拒绝。原因有三种COMSOL Server 没启动、端口号不对、防火墙拦了本地回环。排查顺序是先看 server 启动日志里显示的端口再用telnet或者 MATLAB 的tcpclient试连一下确认端口通不通。server 启动时如果提示端口被占用它会自动顺延到下一个端口脚本里硬编码 2036 就会连错。稳妥做法是把端口号也放进params.m或者让comsol_startup支持从环境变量读。现象七几何建完周期边界套到了错误的边上。前面提过硬编码边编号是主要原因。COMSOL 的边编号在几何重建、布尔操作之后会变脚本里写死[1 3]这次能跑下次改个圆孔位置就套错。解决办法是按坐标筛选边或者建模之后先把所有边的坐标打印出来核对一遍再设边界。这个坑很隐蔽因为套错边之后求解器不一定报错只是结果不对排查起来费时间。现象八MATLAB 版本和 COMSOL LiveLink 不兼容。表现是ModelUtil.create或者mphglobal报方法不存在。原因通常是 MATLAB 太新COMSOL 的 LiveLink 还没跟上。解决办法是查 COMSOL 安装目录下的兼容说明换一个匹配的 MATLAB 版本或者反过来升级 COMSOL。这个问题没有脚本层面的绕法只能对齐版本。6. 二次开发改晶格、加缺陷、批量扫参数6.1 换晶格从正方到三角正方晶格的元胞是正方形高对称点是 Γ-X-M-Γ。换成三角晶格元胞变成菱形或者等效的六边形高对称点变成 Γ-M-K-Γ不可约布里渊区的形状也变了。改的时候不需要重写求解循环只需要换三样东西build_unit_cell.m里的几何菱形元胞加两个圆孔或者一个圆孔、high_symmetry_path.m里的高对称点坐标、set_physics.m里的周期边界对菱形有两组对边方向不是正交的。% 三角晶格的高对称点归一化单位 2*pi/a kG [0, 0]; kM [0.5, 0]; kK [1/3, 1/sqrt(3)]; % Γ-M-K-Γ 路径注意 K 点坐标和正方晶格不同几何换完之后周期边界的法向也跟着变kFloquet的分量不能再简单给{kx,0,0}和{0,ky,0}要把波矢分解到两组对边的法向方向上。这一步容易翻车建议先用一个已知带隙的三角晶格参数跑一遍和文献里的归一化频率对比确认边界方向没设错再往下做。6.2 加缺陷点缺陷微腔和线缺陷波导光子晶体带隙算出来之后下一步往往是引入缺陷。点缺陷的做法是去掉中心一个孔或者把某个孔的半径改小破坏局部周期性在带隙里形成一个局域模。线缺陷是去掉一整排孔形成一条波导带隙里的光被限制在缺陷线上传播。在脚本层面点缺陷的实现是在build_unit_cell.m里加一个条件判断如果p.defect true中心孔不建或者半径乘一个系数。但要注意单胞模型加缺陷之后周期性被破坏不能再用一个元胞加 Floquet 边界得用超胞supercell——把 $N \times N$ 个元胞拼成一个大元胞缺陷放在中心大元胞边界上仍然用 Floquet 条件。超胞的代价是几何变大、网格变多、求解变慢$本文还有配套的精品资源点击获取
简介「photonic-bandgap-comsol-matlab」面向光学、光子晶体方向的研究生与工程人员聚焦如何借助 COMSOL Multiphysics 的光波光学模块配合 MATLAB 脚本完成二维光子晶体能带结构与光子带隙的自动化计算。压缩包共 4 个文件、约 144KB含 MATLAB 主脚本、fig 图形文件、结果对比图与说明文档主脚本承担建模、参数扫描与结果读取其余文件分别用于保存能带曲线、核对仿真输出与解释使用方式。目前已有 1066 人学习下载。读者可从中了解 COMSOL 与 MATLAB 联合调用的接口写法参考晶格常数、单元形状、PML 边界与频率范围的组织方式并借助结果对比图核对带隙位置与宽度为光通信、光学微腔等场景的器件设计积累可复用的脚本模板与调试思路。1. 带隙算出来是负的这套 COMSOL-MATLAB 脚本把哪几步固定住了二维光子晶体的带隙算出来是负的或者相邻两条能带在某个 k 点莫名其妙交叉多半不是物理错了而是周期边界、波矢扫描和模式排序这三步里有一环没固定住。这份 photonic-bandgap-comsol-matlab 脚本包干的正是这件事用 MATLAB 通过 LiveLink 驱动 COMSOL把二维光子晶体的元胞建模、Floquet 周期边界、布里渊区高对称路径扫描、特征频率求解和带隙判定串成一条可重复执行的流水线。它适合正在做光子晶体能带、微腔或波导设计又不想每次都在 GUI 里手点参数的人。改几个参数就能跑出 Γ-X-M-Γ 路径上的 band structure省下的是反复调边界条件的时间。2. 技术栈与目录结构MATLAB 当控制层COMSOL 当计算引擎2.1 LiveLink 的分工为什么不在 GUI 里点COMSOL 本身能算光子晶体能带GUI 里也能设周期边界和参数扫描。问题出在“重复”和“记录”上每换一个晶格常数、每换一个孔半径都要重新点一遍几何、材料、边界、研究设置点错一个边编号结果就偏了。LiveLink for MATLAB 把 COMSOL 变成一个可以被脚本调用的计算引擎MATLAB 负责参数、循环、后处理和画图COMSOL 只负责求解。分工清楚之后整个计算过程就是一段可版本管理的代码而不是一堆散落在模型树里的设置。这个项目的技术栈不复杂但每一层都有明确的职责层级工具职责控制层MATLAB参数管理、k 点循环、结果收集、带隙判定、画图求解层COMSOL Multiphysics二维频域电磁波、Floquet 周期边界、特征频率求解接口层LiveLink for MATLAB把 COMSOL 模型对象暴露成 MATLAB 可调用的 API数据层CSV / MAT保存每个 k 点的特征频率和带隙表常见做法是把所有可调参数收进一个params.m主脚本只负责按顺序调用建模、求解、后处理三个函数。这样二次开发时不用翻主脚本改参数文件就行。2.2 目录结构与每个文件的职责源码包的目录组织大概是这个形状不同版本可能略有差异但职责划分基本一致photonic-bandgap-comsol-matlab/ ├── main_bandgap.m # 主入口加载参数、建模、扫描、出图 ├── config/ │ └── params.m # 所有可调参数集中在这里 ├── src/ │ ├── build_unit_cell.m # 元胞几何与材料 │ ├── set_physics.m # 频域电磁波与 Floquet 周期边界 │ ├── set_mesh.m # 网格划分 │ ├── run_k_sweep.m # k 点循环与特征频率提取 │ └── find_bandgaps.m # 带隙判定 ├── models/ │ └── pc_2d_base.mph # 可选的空白模型模板 ├── results/ │ ├── band_structure.csv # 归一化频率矩阵 │ └── figures/ # 输出的能带图 └── utils/ ├── comsol_startup.m # 启动 COMSOL Server 并连接 └── high_symmetry_path.m # 生成 Γ-X-M-Γ 路径这个结构的好处是每个文件只做一件事。build_unit_cell.m只管几何set_physics.m只管物理场和边界run_k_sweep.m只管循环和取数。排查问题时可以单独把某一段拎出来跑不用每次都从头执行整个主脚本。我一般会把results/整个目录加到.gitignore只保留params.m和src/下的脚本这样换机器时不会把几 GB 的中间结果带过去。2.3 环境准备启动顺序和版本匹配COMSOL 和 MATLAB 的联动有一个固定的启动顺序顺序错了就会卡在连接那一步。常见做法是先在系统里启动 COMSOL Multiphysics Server记下端口号默认 2036再在 MATLAB 里调用mphstart连接。也可以在 MATLAB 里直接启动 server但需要把 COMSOL 的bin目录加到系统路径里。% utils/comsol_startup.m function model comsol_startup(port) % 连接已经启动的 COMSOL Server返回一个空的模型容器 if nargin 1 port 2036; % COMSOL Server 默认端口 end import com.comsol.model.util.* try mphstart(port); catch ME error(连接 COMSOL Server 失败先确认 server 已启动、端口 %d 没被占用%s, ... port, ME.message); end ModelUtil.showProgress(true); % 在 MATLAB 命令行打印求解进度 model ModelUtil.create(Model); end这段代码的逻辑很直接mphstart负责建立 MATLAB 到 COMSOL Server 的 TCP 连接ModelUtil.create创建一个空的模型对象后面的几何、物理场、研究都挂在这个对象上。ModelUtil.showProgress(true)不是必须的但求解特征频率时能看到进度条长时间扫描不容易以为程序卡死。参数说明port要和 COMSOL Server 启动时显示的端口一致一台机器上同时跑多个 COMSOL 实例时端口会顺延不要硬编码。如果mphstart报连接超时先检查 server 是否真的在运行再检查防火墙有没有拦本地回环端口。这一步翻车的人不少血泪经验是先把 server 启动日志看一眼比在 MATLAB 里反复试端口快得多。版本匹配也值得提一句COMSOL 的 LiveLink 接口和 MATLAB 版本之间有兼容矩阵太新的 MATLAB 配太旧的 COMSOL 可能在ModelUtil那一层就报方法找不到。拿到源码包之后先确认本机 COMSOL 版本支持当前 MATLAB再往下走。3. 物理建模与核心脚本从元胞到 Bloch 边界的完整链路3.1 元胞、晶格常数与不可约布里渊区二维光子晶体的能带计算第一步是把无限周期结构截成一个元胞。正方晶格的元胞取边长等于晶格常数 $a$ 的正方形中心挖一个圆孔或者放一根介质柱孔半径或柱半径记作 $r$。介电常数在元胞内部分布背景介质一个值散射体一个值。这个元胞在二维平面内按 $a$ 平移复制就得到整个晶体所以求解域只需要一个元胞边界上用周期条件把相邻元胞“接”起来。晶格常数 $a$ 决定了归一化频率的尺度。带隙位置通常用 $a/\lambda$ 或 $\omega a/2\pi c$ 表示这样换材料、换晶格常数时能带图可以直接对比。正方晶格的高对称点是 Γ(0,0)、X(0.5,0)、M(0.5,0.5)单位是 $2\pi/a$。不可约布里渊区边界就是 Γ-X-M-Γ 这条折线扫描这条线上的 k 点就能覆盖带隙判定需要的信息。% utils/high_symmetry_path.m function [kx, ky, s] high_symmetry_path(nk) % 生成正方晶格不可约布里渊区边界 Γ-X-M-Γ 上的 k 点 % 输入 nk每条线段上的点数含端点 % 输出 kx, ky归一化波矢单位 2*pi/a % 输出 s路径序号画图时当横坐标用 kG [0.0, 0.0]; kX [0.5, 0.0]; kM [0.5, 0.5]; seg1 [linspace(kG(1), kX(1), nk), linspace(kG(2), kX(2), nk)]; seg2 [linspace(kX(1), kM(1), nk), linspace(kX(2), kM(2), nk)]; seg3 [linspace(kM(1), kG(1), nk), linspace(kM(2), kG(2), nk)]; % 去掉相邻线段重复的端点 kpts [seg1; seg2(2:end, :); seg3(2:end, :)]; kx kpts(:, 1); ky kpts(:, 2); s (0:size(kpts, 1) - 1); end这段函数的参数只有nk它决定每条高对称线上的采样密度。nk 21是一个常见起点三条线段去重之后一共 61 个 k 点。点数太少会漏掉能带交叉附近的细节带隙边界会算宽点数太多单次扫描时间线性增长。实际做法是先跑nk 11看整体形状确认带隙大致位置后再用nk 31或51做精细扫描。s是给画图用的横坐标能带图的横轴不是真实的 k 值而是路径上的累计序号在 Γ、X、M 三个位置打刻度。3.2 Floquet 周期边界怎么落到脚本里周期边界是这个项目里最容易设错的一环。COMSOL 的Periodic边界条件有两种类型连续型和 Floquet 型。光子晶体能带计算必须用 Floquet 型因为它允许场在跨越元胞时带一个相位因子 $e^{i\mathbf{k}\cdot\mathbf{a}}$这个相位因子就是 Bloch 波矢。用连续型周期边界等于强制 $k0$扫出来的不是能带是 Γ 点的几个驻模。脚本里设置 Floquet 周期的关键是把波矢分量写成全局参数左右一对边界给 $k_x$上下对给 $k_y$% src/set_physics.m节选 function model set_physics(model, p) % 二维频域电磁波特征频率求解 model.physics.create(ewfd, ElectromagneticWavesFrequencyDomain, geom1); % 背景介质 model.physics(ewfd).material(mat1).set(relpermittivity, {num2str(p.eps_host)}); % 散射体介质域编号按实际几何调整 model.physics(ewfd).material(mat2).selection.set(p.hole_domains); model.physics(ewfd).material(mat2).set(relpermittivity, {num2str(p.eps_hole)}); % Floquet 周期边界左右一对 model.physics(ewfd).create(pcx, Periodic, 1); model.physics(ewfd).feature(pcx).selection.set(p.left_right_edges); model.physics(ewfd).feature(pcx).set(PeriodicType, Floquet); model.physics(ewfd).feature(pcx).set(kFloquet, {kx, 0, 0}); % Floquet 周期边界上下对 model.physics(ewfd).create(pcy, Periodic, 1); model.physics(ewfd).feature(pcy).selection.set(p.bottom_top_edges); model.physics(ewfd).feature(pcy).set(PeriodicType, Floquet); model.physics(ewfd).feature(pcy).set(kFloquet, {0, ky, 0}); end这里的kx和ky不是数值是模型里的参数名实际值在扫描循环里用model.param.set逐点写入。kFloquet的三个分量对应三维二维模型里第三个分量填 0。p.left_right_edges和p.bottom_top_edges是边的编号不建议硬编码成[1 3]这种字面量因为几何一变编号就漂。稳妥的做法是在params.m里按坐标筛选边或者建模之后用mphgetadj之类的接口按位置取边。硬编码边编号是新手最容易踩的坑之一换个 COMSOL 版本或者改一下几何顺序周期边界就套到错误的边上求解不报错但结果是错的。还有一个细节Floquet 周期边界的相位因子方向要和波矢方向一致。左右边界的法向是 $x$所以给 $k_x$上下边界的法向是 $y$所以给 $k_y$。如果两边都给同一个分量相当于在某个方向上强制 $k0$带隙会算得偏高或偏低。3.3 特征频率求解器与网格设置COMSOL 里算能带用的是特征频率研究Eigenfrequency不是频域扫描。频域扫描给定一个频率求场分布特征频率求解给定边界条件求系统的本征频率后者才对应能带上的一个个模式。研究设置里要指定求解多少个特征频率neigs以及一个频移shift。频移的作用是把求解器引到感兴趣的频率附近避免它去挖零模或者极低频的伪模。% 研究设置特征频率围绕目标频段做频移 model.study.create(std1); model.study(std1).create(eig, Eigenfrequency); model.study(std1).feature(eig).set(neigs, p.nmodes); model.study(std1).feature(eig).set(shift, p.fshift); model.study(std1).feature(eig).set(eigwhich, lm); % 最接近 shift 的模参数说明p.nmodes是每个 k 点求解的模式数至少要比你关心的能带数多两三个留出余量给伪模p.fshift是频移单位 Hz通常设在归一化频率 $a/\lambda \approx 0.3$ 对应的实际频率附近eigwhich lm表示找最接近频移的特征值比默认的lr最大实部更适合带隙计算因为带隙附近的模式不一定实部最大。网格设置直接决定频率算得准不准。光子晶体元胞里有圆孔孔边缘是曲面网格太粗会让特征频率整体偏移。常见做法是给孔边界加一层加密同时限制全局最大单元尺寸% src/set_mesh.m function model set_mesh(model, p) model.mesh.create(mesh1, geom1); model.mesh(mesh1).autoMeshSize(4); % 整体“较细” model.mesh(mesh1).feature(size).set(hmax, p.hmax); model.mesh(mesh1).feature(size).set(hmin, p.hmax / 10); % 孔边界加密 model.mesh(mesh1).create(edg1, Size); model.mesh(mesh1).feature(edg1).selection.set(p.hole_edges); model.mesh(mesh1).feature(edg1).set(hmax, p.hmax / 3); model.mesh(mesh1).run; endp.hmax一般取 $a/20$ 到 $a/30$。这个值不是拍脑袋定的后面第 4 章会做网格收敛测试把hmax从 $a/10$ 逐步降到 $a/40$看带隙中心频率偏移是否小于 1%。如果偏移还在漂说明网格还不够细。孔边界的局部加密比全局加密划算因为场在孔边缘变化最剧烈元胞内部相对平缓。3.4 参数入口与几何构建把参数集中到一个函数里是这个脚本包最值得抄的做法。光子晶体带隙对 $r/a$ 和介电常数对比很敏感二次开发时经常要扫这两个量。参数散落在各个脚本里扫到第三组就忘了第一组改过什么。集中之后主脚本只需要一行p params();所有函数共享同一个结构体。% config/params.m function p params() % 二维光子晶体带隙计算的参数入口 p.a 1e-6; % 晶格常数单位 m p.r 0.30 * p.a; % 孔半径0.3a 是常见起点 p.eps_host 12.0; % 背景介质相对介电常数 p.eps_hole 1.0; % 孔内相对介电常数 p.pol TE; % 偏振约定TE 表示电场沿面外 p.nk 21; % 每条高对称线上的点数 p.nmodes 8; % 每个 k 点求解的模式数 p.fshift 0.3 * 3e8 / p.a; % 频移对应 a/lambda 约 0.3 p.hmax p.a / 20; % 全局最大网格尺寸 p.hole_edges []; % 孔边界编号建模后按坐标填 p.left_right_edges []; p.bottom_top_edges []; p.outdir fullfile(pwd, results); end几何构建的函数只做三件事建正方形元胞、建圆孔、做布尔差集。圆孔的位置放在元胞中心这样周期平移之后孔阵列是对齐的。如果孔放在角上平移之后会变成四个四分之一孔拼在一起几何上等价但边界条件容易设错。% src/build_unit_cell.m function model build_unit_cell(model, p) model.geom.create(geom1, 2); % 正方形元胞 model.geom(geom1).create(cell, Rectangle); model.geom(geom1).feature(cell).set(size, {a, a}); model.geom(geom1).feature(cell).set(pos, [0 0]); % 中心圆孔 model.geom(geom1).create(hole, Circle); model.geom(geom1).feature(hole).set(r, r); model.geom(geom1).feature(hole).set(pos, {a/2, a/2}); % 布尔差集元胞减去孔 model.geom(geom1).create(dif1, Difference); model.geom(geom1).feature(dif1).selection(input).set({cell}); model.geom(geom1).feature(dif1).selection(input2).set({hole}); model.geom(geom1).run; % 按坐标筛选边避免硬编码编号 p.left_right_edges mphselectbox(model, geom1, [0 0; 0 p.a], boundary); p.bottom_top_edges mphselectbox(model, geom1, [0 0; p.a 0], boundary); endmphselectbox是按坐标框选几何实体的接口比记边编号可靠。[0 0; 0 p.a]框住的是左边那条边[0 0; p.a 0]框住的是下边那条边。实际脚本里可能用mphgetadj或直接遍历边坐标逻辑是一样的位置定边不靠编号。3.5 k 点扫描与后处理判据扫描循环是整个脚本里最耗时的部分。每个 k 点都要重新设置kx、ky重新求解特征频率然后把频率取出来存进矩阵。循环里最容易出错的是单位换算high_symmetry_path返回的是归一化波矢单位是 $2\pi/a$写进 COMSOL 参数时要乘上 $2\pi/a$ 换成 $1/m$。% src/run_k_sweep.m核心循环 function freq run_k_sweep(model, p, kx, ky) npts numel(kx); freq zeros(p.nmodes, npts); for i 1:npts % 归一化波矢转实际波矢单位 1/m kx_real kx(i) * 2 * pi / p.a; ky_real ky(i) * 2 * pi / p.a; model.param.set(kx, [num2str(kx_real), [1/m]]); model.param.set(ky, [num2str(ky_real), [1/m]]); model.study(std1).run; % 取特征频率实部虚部对应损耗无损耗时接近 0 f mphglobal(model, ewfd.freq, dataset, dset1); f sort(real(f(:))); % 模式数不够时补 NaN避免矩阵维度不一致 if numel(f) p.nmodes f(end1:p.nmodes) NaN; end freq(:, i) f(1:p.nmodes); end % 归一化频率 a/lambda f * a / c freq freq * p.a / 3e8; end这段循环的逻辑外层遍历 k 点内层由 COMSOL 求解器完成。mphglobal从结果数据集里取全局标量ewfd.freq返回的是一个列向量。sort保证每个 k 点的模式按频率从低到高排后处理判带隙时才能按行对应能带。归一化那一步放在循环外面一次性把 Hz 转成 $a/\lambda$省得在循环里反复乘除。后处理判带隙的判据不复杂把频率矩阵按行看成一条条能带第 $n$ 条带在所有 k 点上的最大值如果小于第 $n1$ 条带在所有 k 点上的最小值两条带之间就有带隙。带隙宽度是两者之差中心频率是两者之和的一半。% src/find_bandgaps.m function gaps find_bandgaps(freq_norm, min_width) % freq_norm归一化频率矩阵行是模式列是 k 点 % min_width带隙最小宽度低于这个值不算 nb size(freq_norm, 1); gaps []; for n 1:nb-1 upper max(freq_norm(n, :)); % 第 n 条带的上边界 lower min(freq_norm(n1, :)); % 第 n1 条带的下边界 if lower upper width lower - upper; center 0.5 * (lower upper); if width min_width gaps [gaps; n, upper, lower, width, center]; end end end end这个判据的前提是能带已经按频率排好序且没有出现“带交叉”——两条能带在某个 k 点交换了顺序。正方晶格的光子晶体在低阶带里交叉不常见但高阶带或者参数接近带隙关闭时会出现。稳妥做法是在排序之后再做一次连续性跟踪每个 k 点选择与上一个 k 点频率最接近的模式作为同一条带而不是简单按频率排序。脚本包里如果只做了sort高阶带的结果要打个问号这一点在第 5 章会展开。4. 跑通与结果验证一条命令跑出 Γ-X-M-Γ 能带4.1 最短跑通路径拿到源码包之后最短路径是四步启动 COMSOL Server、在 MATLAB 里cd到项目根目录、改config/params.m里的a和r、运行main_bandgap。主脚本的职责就是把前面几节的函数串起来并在最后把结果写盘。% main_bandgap.m clear; clc; addpath(genpath(fullfile(pwd, src))); addpath(genpath(fullfile(pwd, utils))); addpath(genpath(fullfile(pwd, config))); p params(); % 1. 读参数 model comsol_startup(); % 2. 连 COMSOL Server model build_unit_cell(model, p); % 3. 建元胞 model set_physics(model, p); % 4. 物理场与周期边界 model set_mesh(model, p); % 5. 网格 [kx, ky, s] high_symmetry_path(p.nk);% 6. 高对称路径 freq run_k_sweep(model, p, kx, ky); % 7. 扫描求解 gaps find_bandgaps(freq, 0.01); % 8. 判带隙 writematrix(freq, fullfile(p.outdir, band_structure.csv)); disp(gaps);第一次跑建议把p.nk设成 11p.nmodes设成 6先确认整条链路能走通、结果形状对再加密。main_bandgap.m里没有画图代码画图单独放一个函数或者直接在命令行里plot(s, freq., .-)这样求解和展示解耦批量跑的时候不会因为画图拖慢速度。跑通之后results/band_structure.csv里是一个nmodes × npts的矩阵每列对应一个 k 点每行对应一条能带。gaps矩阵每一行是一条带隙列依次是能带序号、上边界、下边界、宽度、中心频率。4.2 band structure 图怎么读能带图的横轴是 Γ-X-M-Γ 路径上的累计序号纵轴是归一化频率 $a/\lambda$。每个 k 点上有若干个点连起来就是一条条能带。带隙出现在两条能带之间没有点的频率区间图上表现为一条横向的空白带。看这张图的时候先确认三件事Γ 点横轴起点的频率是否合理X 点和 M 点的频率是否有明显的分裂带隙上下边界是否平直。如果 Γ 点的频率接近 0大概率是伪模混进来了如果 X 点和 M 点的频率几乎一样可能是周期边界的方向设反了如果带隙边界随 k 点剧烈抖动说明 k 点太稀或者网格不收敛。TE 和 TM 偏振的带隙位置不一样这是正常现象。二维光子晶体里空气孔阵列的 TE 带隙通常比 TM 带隙更明显介质柱阵列反过来。脚本里p.pol这个参数最终要落到物理场的偏振设置上TE 对应电场沿面外二维模型里可以用面外电场分量TM 对应磁场沿面外。如果脚本里只设了一个偏振换偏振时要改物理场设置不能只改参数名。4.3 网格与 k 点收敛验证带隙中心频率算出来是 0.42 还是 0.40取决于网格和 k 点够不够密。验证方法不复杂但必须做固定其他参数只变hmax看带隙中心频率怎么变固定hmax只变nk看带隙宽度怎么变。网格尺寸带隙中心 $a/\lambda$带隙宽度相对上一次变化a/100.4180.031—a/150.4120.0361.4%a/200.4090.0380.7%a/300.4080.0390.2%a/400.4080.039 0.1%这张表是典型收敛曲线网格从 $a/10$ 加密到 $a/20$中心频率偏移接近 2%从 $a/20$ 到 $a/30$ 偏移降到 0.2% 以下再加密变化很小。实际项目里取 $a/20$ 到 $a/30$ 之间的值兼顾精度和速度。k 点收敛类似nk从 11 加到 21带隙宽度会收窄一点因为边界 k 点采得更密能带极值找得更准从 21 加到 41 变化通常在 1% 以内。提示收敛测试不要只做一次。换晶格类型、换填充比、换介电常数之后收敛行为会变。填充比越高、介电常数对比越大需要的网格越细。4.4 结果导出与复用band_structure.csv是最值得保存的中间结果。有了它画图、判带隙、对比不同参数都不用重新跑 COMSOL。二次开发时常见做法是把params.m里的关键参数和band_structure.csv一起打包文件名带上r_over_a和eps比如band_r0.30_eps12.csv这样扫参数扫出来的几十组结果不会混在一起。gaps矩阵可以单独存成bandgaps.csv每行一条带隙列名固定。后面要做带隙随 $r/a$ 变化的相图直接读这个文件就行不用再碰 COMSOL。导出的时候注意归一化频率的单位统一别一组用 $a/\lambda$、一组用 THz混在一起画图会得到一条斜率莫名其妙的曲线。5. 避坑与排查特征频率求解最常见的五类翻车5.1 求解阶段的三个高频问题现象一每个 k 点都解出一堆接近 0 的频率。原因通常是频移设得太低求解器把静模和伪模一起挖出来了。COMSOL 的特征频率求解器在频移附近找模频移接近 0 时零空间里的伪模就会被选中。解决办法是把p.fshift设到目标带隙附近比如归一化频率 0.3 对应的实际频率同时在取结果时过滤掉 $a/\lambda 0.05$ 的点。如果过滤之后模式数不够把p.nmodes加大用数量换有效模式。现象二扫描到某个 k 点突然报“找不到特征频率”。原因多半是频移落在了两条能带之间的带隙里那个频率附近确实没有本征模。带隙计算本来就会经过带隙区间求解器在带隙里找不到模是正常的但脚本如果没做容错就会中断。解决办法是在循环里包一层try-catch求解失败时给该 k 点填NaN继续下一个点。后处理判带隙时跳过含NaN的列或者把NaN当作“该模式不存在”处理。现象三同一组参数跑两次带隙宽度差了几个百分点。原因可能是网格自适应每次生成的网格不完全一致也可能是特征频率求解器对初值敏感。解决办法是关掉自适应网格重划分固定hmax和边界加密层数求解器那边把eigwhich固定成lm不要用默认设置。如果还漂把p.nmodes加大到能覆盖带隙上下各两条带减少模式漏选的概率。5.2 后处理与结果解读的坑现象四带隙表里出现负宽度或者相邻带隙重叠。这是模式排序错乱导致的。简单sort在每个 k 点独立排序遇到能带交叉时会把两条带的顺序交换后处理按行相减就得到负值。解决办法是做能带跟踪从第一个 k 点开始每个新模式选择与上一个 k 点同序号模式频率最接近的那个而不是无条件按频率排序。跟踪之后再做带隙判定负宽度基本消失。如果嫌跟踪麻烦至少在判带隙前检查相邻两行的频率差是否连续出现跳变就标记出来人工看一眼。现象五TE 和 TM 的带隙位置几乎一样。二维光子晶体的 TE 和 TM 带隙通常不在同一个频率两个偏振算出来一样大概率是物理场设置里偏振没真正分开。检查set_physics.m里有没有设置面外分量TE 需要电场只有面外分量TM 需要磁场只有面外分量。如果两个偏振用的是同一套边界和材料设置只是改了p.pol这个字符串求解结果不会变。偏振要落到物理场的分量约束或者端口设置上不能只改参数名。5.3 环境与接口的坑现象六mphstart报连接被拒绝。原因有三种COMSOL Server 没启动、端口号不对、防火墙拦了本地回环。排查顺序是先看 server 启动日志里显示的端口再用telnet或者 MATLAB 的tcpclient试连一下确认端口通不通。server 启动时如果提示端口被占用它会自动顺延到下一个端口脚本里硬编码 2036 就会连错。稳妥做法是把端口号也放进params.m或者让comsol_startup支持从环境变量读。现象七几何建完周期边界套到了错误的边上。前面提过硬编码边编号是主要原因。COMSOL 的边编号在几何重建、布尔操作之后会变脚本里写死[1 3]这次能跑下次改个圆孔位置就套错。解决办法是按坐标筛选边或者建模之后先把所有边的坐标打印出来核对一遍再设边界。这个坑很隐蔽因为套错边之后求解器不一定报错只是结果不对排查起来费时间。现象八MATLAB 版本和 COMSOL LiveLink 不兼容。表现是ModelUtil.create或者mphglobal报方法不存在。原因通常是 MATLAB 太新COMSOL 的 LiveLink 还没跟上。解决办法是查 COMSOL 安装目录下的兼容说明换一个匹配的 MATLAB 版本或者反过来升级 COMSOL。这个问题没有脚本层面的绕法只能对齐版本。6. 二次开发改晶格、加缺陷、批量扫参数6.1 换晶格从正方到三角正方晶格的元胞是正方形高对称点是 Γ-X-M-Γ。换成三角晶格元胞变成菱形或者等效的六边形高对称点变成 Γ-M-K-Γ不可约布里渊区的形状也变了。改的时候不需要重写求解循环只需要换三样东西build_unit_cell.m里的几何菱形元胞加两个圆孔或者一个圆孔、high_symmetry_path.m里的高对称点坐标、set_physics.m里的周期边界对菱形有两组对边方向不是正交的。% 三角晶格的高对称点归一化单位 2*pi/a kG [0, 0]; kM [0.5, 0]; kK [1/3, 1/sqrt(3)]; % Γ-M-K-Γ 路径注意 K 点坐标和正方晶格不同几何换完之后周期边界的法向也跟着变kFloquet的分量不能再简单给{kx,0,0}和{0,ky,0}要把波矢分解到两组对边的法向方向上。这一步容易翻车建议先用一个已知带隙的三角晶格参数跑一遍和文献里的归一化频率对比确认边界方向没设错再往下做。6.2 加缺陷点缺陷微腔和线缺陷波导光子晶体带隙算出来之后下一步往往是引入缺陷。点缺陷的做法是去掉中心一个孔或者把某个孔的半径改小破坏局部周期性在带隙里形成一个局域模。线缺陷是去掉一整排孔形成一条波导带隙里的光被限制在缺陷线上传播。在脚本层面点缺陷的实现是在build_unit_cell.m里加一个条件判断如果p.defect true中心孔不建或者半径乘一个系数。但要注意单胞模型加缺陷之后周期性被破坏不能再用一个元胞加 Floquet 边界得用超胞supercell——把 $N \times N$ 个元胞拼成一个大元胞缺陷放在中心大元胞边界上仍然用 Floquet 条件。超胞的代价是几何变大、网格变多、求解变慢$本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
