COMSOL-MATLAB二维光子晶体带隙脚本：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 只负责求解。分工清楚之后，整个计算过程就是一段可版本管理的代码，而不是一堆散落在模型树里的设置。

这个项目的技术栈不复杂，但每一层都有明确的职责：

常见做法是把所有可调参数收进一个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 波矢。用连续型周期边界等于强制 $k=0$，扫出来的不是能带，是 Γ 点的几个驻模。

脚本里设置 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$。如果两边都给同一个分量，相当于在某个方向上强制 $k=0$，带隙会算得偏高或偏低。

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; end

p.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'); end

mphselectbox是按坐标框选几何实体的接口，比记边编号可靠。[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(end+1: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 点上的最大值，如果小于第 $n+1$ 条带在所有 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(n+1, :)); % 第 n+1 条带的下边界 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设成 11，p.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/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 条件。超胞的代价是几何变大、网格变多、求解变慢，$

本文还有配套的精品资源，点击获取