简介面向光学仿真学习者与二维光子晶体研究者的COMSOL-MATLAB联合计算示例围绕周期性结构的光子带隙求解展开适合已具备一定电磁场与脚本基础、希望把建模流程自动化的读者。资源包共4个文件含MATLAB脚本、FIG图形、PNG结果图与Markdown说明约144KB分别对应仿真控制、结果可视化、效果预览和项目说明便于先看说明再复现带隙计算与能带图。已有1066人学习说明该方向具有一定关注度。读者可从中了解如何用脚本设置晶格与频率范围、调用COMSOL求解并读取场分布或S参数再据此识别带隙边界、比较结果并优化设计获得一套可拆解、可复用的计算思路。1. 从一张空白的能带图说起这套脚本补的是哪一段正方晶格、介质柱、折射率差拉满物理课本上三行字就能讲完的 2D 光子晶体带隙真正落到电脑上往往卡在第一张能带图上单胞怎么切、周期条件挂在哪两条边上、k 点该给弧度还是给无量纲、特征频率求解器为什么吐出一堆负数和复数。photonic-bandgap-comsol-matlab 这套 COMSOL 与 MATLAB 的联合脚本价值不在算法有多新而在于把 2D 光子晶体带隙计算这条链路上的建模、扫 k、判带隙三步固定成了可重复执行的代码。适合两类人一是刚接触 photonic crystals 仿真、被 GUI 点选绕晕的研究生二是需要批量跑参数、画带隙地图又不想每次重头点一遍的从业者。读完你应该能自己跑出 Γ-X-M-Γ 路径上的能带图并知道哪几个参数一改结果就变脸。2. 计算链拆解COMSOL 单胞建模与 MATLAB 驱动的分工先把这条链路的职责分清楚后面写代码才不会一边建模一边怀疑物理。COMSOL 负责的是有限元本征值求解这件重活离散化单胞、装配 Bloch 周期边界、解广义本征值方程MATLAB 负责的是循环、判据、出图和落盘。两者之间走 LiveLink 接口MATLAB 手里握着的是同一个 model 对象的句柄改一个参数不需要存盘再导入。2.1 为什么带隙计算要坚持脚本化而不是 GUI 点选单点求解用 GUI 五分钟能点完问题在于带隙从来不是一个点的事。一条 Γ-X-M-Γ 路径按每段 20 个点算就是 60 次本征值求解再叠上 r/a 从 0.10 扫到 0.45规模立刻变成几百次。GUI 在这种量级下有兩個硬伤一是每次改参数都要重新点一遍几何、材料、网格出错概率随次数线性上升二是结果散落在各个节点里事后想复原“这张图是哪个参数跑的”只能靠记忆。脚本化解决的正是这两点。参数写在函数签名里一次调用就是一次完整实验想复现某张图只要把那次调用的参数抄回来。我在实际项目里的一般做法是建模封成一个pc_build_cell函数求解扫点封成另一个函数中间的模型对象不销毁、只改属性这样一次装配多次复用速度能差出好几倍。提示LiveLink 的连接是有状态的MATLAB 里ModelUtil.tags能看到当前挂着哪些模型跑批之前先清一遍避免拿到上一次的残留句柄。这里还有个选型上的判断如果你的目标是快速看个大概COMSOL 自带的周期性结构案例模板改一改更快但只要涉及参数扫描、带隙宽度统计、和平面波展开法交叉验证脚本就是唯一不折腾的路。别指望 GUI 里存个 model 文件然后用 MATLAB 打开改参数那条路的版本兼容坑比想象中多。2.2 单胞里到底放了什么几何、材料与 Bloch 相位正方晶格的单胞就是一个边长 a 的正方形中心放一根半径 r 的圆柱。这里有个容易忽略的几何前提只有当 r/a 0.5 时柱体才完全不触碰单胞边界边界落在均匀背景介质里Floquet 周期条件挂上去才干净。如果你的参数扫到 r/a 0.48 附近柱体已经贴边网格在边界处的质量和周期配对都会变得微妙带隙位置会开始抖——这不是物理是离散误差。材料上只需两组数柱体相对介电常数 εr常见取值 11.4 左右对应近红外波段的硅或者 12 附近对应砷化镓背景取 1空气。介质柱和空气孔的带隙位置差别很大脚本里通过交换这两个数就能切换不需要改几何。Bloch 相位是整件事的物理核心。周期边界上场的相位差是 exp(i k·a)COMSOL 里对应 Floquet 周期条件需要给的是波矢分量 kx、ky单位 rad/m不是无量纲的 ka。这一点是新手最容易翻车的地方给成 ka 会让整条能带图横向压缩 1/a 倍看着像带隙消失了其实是横轴错了。几何与材料的落点整理成表更直观对象参数典型取值影响单胞边长 a1无量纲或 1e-6 m只影响绝对频率不影响归一化带隙柱体半径 r0.10a ~ 0.45a决定带隙是否出现及宽度柱体相对介电常数 εr11.4 / 12折射率差越大带隙越宽背景相对介电常数1与柱体交换即变为空气孔结构周期边界kx, ky0 ~ π/a决定采样在布里渊区的位置把 a 取成 1 是个省事的做法所有长度都变成相对值归一化频率直接就是 a/λ不用每次除来除去代价是材料参数里的长度量纲比如色散模型没法直接写进去。带隙计算通常只关心无色散的介电常数取 a 1 完全够用。2.3 LiveLink 的调用骨架与对象树对应关系COMSOL 的 MATLAB API 是照抄 GUI 的对象树的理解这一点比背函数名重要。model.component(comp1)对应模型开发器里的组件节点下面挂geom、mesh、physics、material、study、sol六个子树路径写法和 GUI 里点开的层级一模一样。所以你在 GUI 里能点到的设置API 里一定有对应的feature和set。装配顺序有讲究几何必须先run完再建物理场否则物理场选不到域材料要在物理场之后赋不然默认材料会覆盖你的设置网格必须最后建因为它依赖几何的最终形态。这个顺序搞反了报的错往往和真正的病因隔着十万八千里。% 建立空模型并挂上 2D 组件 import com.comsol.model.* import com.comsol.model.util.* model ModelUtil.create(Model); % 新建模型tag 可自定义 model.modelNode.create(comp1); % 组件节点 model.geom.create(geom1, 2); % 二维几何 % 正方形单胞左下角在原点 中心圆柱 model.geom(geom1).create(r1, Rectangle); model.geom(geom1).feature(r1).set(size, {a a}); model.geom(geom1).feature(r1).set(pos, {0 0}); model.geom(geom1).create(c1, Circle); model.geom(geom1).feature(c1).set(r, r); model.geom(geom1).feature(c1).set(pos, {a/2 a/2}); model.geom(geom1).run; % 几何必须显式 run否则后续选不到域上面这段只做了一件事把单胞的拓扑定下来。size接受字符串是为了让参数以变量形式存在后面改a、r的值不需要重建几何。pos对矩形是左下角坐标对圆是圆心坐标这两个语义不一样写错一个就会得到偏心的单胞。参数有两种给法一种是在model.param.set(a,1)里登记全局参数几何里引用名字另一种是直接把数值塞进去。前者适合做参数扫描后者适合一次性求解。我一般用前者因为参数扫描时只需要改model.param几何、网格、物理场会自动跟着变省掉大量重建。3. 跑通第一个带隙环境对齐、参数标定与单点求解单点跑通是后面所有工作的地基。这一章的目标很具体定一个 k 点比如 X 点把本征频率解出来确认前几条带的量级落在合理区间。如果这一步的频率已经离谱后面的能带图全是白费功夫。3.1 版本与接口对齐避免第一行就报错LiveLink 的连接方式取决于你的 COMSOL 是客户端模式还是服务端模式。桌面版装了 LiveLink for MATLAB 之后MATLAB 里用mphstart起一个本地服务端口号在 COMSOL 安装目录的配置文件里能查到也可以从 COMSOL 那一侧启动带 MATLAB 的会话这样 MATLAB 打开时连接已经建立好了。两种方式都行前者的好处是脚本可以独立运行。版本错配是新手第一个坑。MATLAB 的版本要比 COMSOL 支持的版本新或者持平反过来容易在import com.comsol.model.*这一行就抛异常。另外 COMSOL 6.x 之后部分物理场字段名有调整比如二维电磁波频域里关于极化方向的属性老脚本抄过来会提示“属性不存在”。遇到这种情况不要硬猜用mphgetproperties把当前对象的属性全打出来看% 打印物理场属性确认字段名与当前版本一致 props mphgetproperties(model.physics(ewfd).prop(Components)); disp(props);这一步看着笨但比在论坛上翻半小时帖子快得多。我现在的习惯是每换一个 COMSOL 大版本先把关键节点的属性名扫一遍把差异记在自己的小本子上后面就能少流点血。3.2 六个决定结果的物理参数参数不是越多越好真正影响带隙位置的其实只有六个其余都是数值细节。下面这张表是我在多次扫描之后固定下来的检查清单每次跑新结构先对一遍参数含义建议起点调错的后果a晶格常数1只改绝对频率刻度r/a填充比0.20决定带隙开闭与宽度εr柱体介电常数11.4折射率差不足则无带隙kx, ky波矢分量0 ~ π/a单位错则能带整体错位shift特征频率搜索中心0.35c/a设偏则漏掉目标带hmax最大网格尺寸a/15过粗则频率抖动shift这一项特别值得说。特征频率求解器本质上是移位反迭代你给的中心频率离真实本征值越近收敛越快、越不容易漏解。给 0.35c/a 是因为正方晶格介质柱结构的第一带隙通常在归一化频率 0.3 到 0.5 之间这个值覆盖了大多数常见填充比。如果你的结构是空气孔、介电常数对比反过来带隙会往低频挪shift 得跟着降。网格尺寸hmax取 a/15 是个折中跑得快但精度一般。真要出图我一般用 a/20 起步关键结论再用 a/30 复核一遍看频率有没有收敛。这个复核步骤后面第 6 章会展开。3.3 从建模到特征频率的完整脚本下面这段是可以直接抄的给定 a、r、εr 和一个 k 点返回该点前若干条带的归一化频率。function fNorm solve_one_k(a, r, epsr, kx, ky, neigs) % 单个 k 点的本征频率求解 % 输入a 晶格常数, r 柱半径, epsr 柱介电常数 % kx,ky 波矢分量(rad/m), neigs 求解模态数 % 输出fNorm 归一化频率 (a/lambda)已排序 import com.comsol.model.* import com.comsol.model.util.* c0 299792458; % 真空光速 model ModelUtil.create(Model); model.modelNode.create(comp1); model.geom.create(geom1, 2); % --- 几何 --- model.geom(geom1).create(r1, Rectangle); model.geom(geom1).feature(r1).set(size, {a, a}); model.geom(geom1).feature(r1).set(pos, {0 0}); model.geom(geom1).create(c1, Circle); model.geom(geom1).feature(c1).set(r, r); model.geom(geom1).feature(c1).set(pos, {a/2, a/2}); model.geom(geom1).run; % --- 材料柱体与背景分开赋 --- model.component(comp1).material.create(mat1, Common); model.component(comp1).material(mat1).selection.set([1]); % 圆域 model.component(comp1).material(mat1).propertyGroup(def)... .set(relpermittivity, {epsr}); model.component(comp1).material.create(mat2, Common); model.component(comp1).material(mat2).selection.set([2]); % 背景域 model.component(comp1).material(mat2).propertyGroup(def)... .set(relpermittivity, {1}); % --- 物理场二维电磁波频域出平面极化 --- model.physics.create(ewfd, ElectromagneticWavesFrequencyDomain, geom1); model.physics(ewfd).selection.all; model.physics(ewfd).prop(Components).set(components, outofplane); % --- 周期边界左右一对、上下各一对 --- model.physics(ewfd).feature.create(pc1, PeriodicCondition, 2); model.physics(ewfd).feature(pc1).selection.set([1 3]); % 左右边 model.physics(ewfd).feature(pc1).set(PeriodicType, Floquet); model.physics(ewfd).feature(pc1).set(kFloquet, {kx, 0}); model.physics(ewfd).feature.create(pc2, PeriodicCondition, 2); model.physics(ewfd).feature(pc2).selection.set([2 4]); % 上下边 model.physics(ewfd).feature(pc2).set(PeriodicType, Floquet); model.physics(ewfd).feature(pc2).set(kFloquet, {0, ky}); % --- 网格 --- model.mesh.create(mesh1, geom1); model.mesh(mesh1).create(ftri1, FreeTri); model.mesh(mesh1).feature(size).set(custom, on); model.mesh(mesh1).feature(size).set(hmax, a/15); model.mesh(mesh1).run; % --- 特征频率研究 --- model.study.create(std1); model.study(std1).create(eig, Eigenfrequency); model.study(std1).feature(eig).set(neigs, neigs); model.study(std1).feature(eig).set(shift, 2*pi*0.35*c0/a); model.study(std1).run; % --- 取解并归一化 --- freq mphglobal(model, freq, dataset, dset1); % 单位 Hz fNorm sort(real(freq)) * a / c0; % 归一化 a/lambda end逐段说清楚几个关键点。几何部分的selection.set([1])依赖域的编号这个编号在不同几何构造顺序下会变稳妥做法是用坐标箱选域mphselectbox而不是写死 1 和 2。周期边界那两行里的selection.set([1 3])同理边界编号也依赖几何第一次跑的时候务必先用mphgetadj或者 GUI 里点一下确认哪条边是 1 号。shift里乘了 2π因为 COMSOL 的角频率和普通频率差一个 2π这个转换漏掉会让搜索中心偏到十万八千里表现为求解器反复报“找不到特征值”。mphglobal取出来的freq是普通频率Hz归一化除以 c/a 得到 a/λ这就是文献里最常用的横轴定义。调用方式和预期输出% 在 X 点 (pi/a, 0) 求解前 8 条带 fx solve_one_k(1, 0.2, 11.4, pi, 0, 8); disp(fx);r/a 0.2、εr 11.4的正方晶格介质柱X 点前几条归一化频率大致落在 0.1 到 0.6 的区间里。如果你拿到的是 3.14 这样的数八成是忘了除 c/a如果是负数或者带虚部八成是 shift 偏了或者网格太粗。4. 沿高对称路径扫描从单点本征值到能带图单点跑通只说明模型没错带隙有没有、多宽得看整条路径。这一章把单点循环成路径再把路径结果折算成带隙判据最后落到文件上。4.1 归一化频率与带隙判据怎么写能带图的横轴是沿路径累积的 k 距离纵轴是归一化频率 a/λ。带隙的判定就是逐条带比较第 n 条带在整条路径上的最大值是否小于第 n1 条带的最小值。这个判据朴素但容易写错——错在拿单个 k 点的值比或者忘了带隙要在整条路径上成立。function [gapW, gapMid, idx] bandgap_from_bands(fNorm) % fNorm: 能带数 x k点数 的归一化频率矩阵已按带排序 % gapW : 相对带隙宽度 (f_upper - f_lower) / f_mid % gapMid: 带隙中心归一化频率 % idx : 带隙位于第 idx 与 idx1 条带之间 [nb, ~] size(fNorm); gapW 0; gapMid NaN; idx 0; for n 1:nb-1 fLow max(fNorm(n, :)); % 下带在整个路径上的最高点 fHigh min(fNorm(n1, :)); % 上带在整个路径上的最低点 if fHigh fLow % 严格重叠才叫带隙 w (fHigh - fLow) / ((fHigh fLow)/2); if w gapW gapW w; gapMid (fHigh fLow)/2; idx n; end end end end这里的相对宽度用的是常见定义Δf / f_mid有些文献用Δf / f_center里的中心换成带隙中点两者差别不大但要前后统一。fLow取 max、fHigh取 min这两个方向写反了会得到“处处有带隙”的荒谬结果我见过不止一次。还要注意一点带隙只有在采样足够密的时候才判得准。路径上只取 20 个点时带的极值可能落在两个采样点之间算出来的宽度偏小。经验值是每段至少 30 个点做正式结论时用 50 个点复核。4.2 k 点路径生成与批量求解循环正方晶格的高对称路径是 Γ(0,0) → X(π/a,0) → M(π/a,π/a) → Γ(0,0)三段长度不同采样点数按段分配更均匀。function [kpath, kdist] kpath_gxm(a, nSeg) % 生成 Gamma-X-M-Gamma 路径每段 nSeg 个点 t linspace(0, 1, nSeg); kx [pi/a*t; pi/a*ones(nSeg,1); pi/a*(1-t)]; ky [zeros(nSeg,1); pi/a*t; pi/a*(1-t)]; kpath [kx ky]; % 累积距离作为横轴 dk sqrt(sum(diff(kpath).^2, 2)); kdist [0; cumsum(dk)]; endkpath是三段拼起来的 N×2 矩阵kdist是沿路径的累积长度画图时用它当横轴X 点和 M 点的位置才对得上。三段拼接处首尾点重复了一次实际使用时可以去重不影响结果。批量求解有两条路。稳妥的一条是老老实实循环每次改kFloquet再run缺点是每次都要重新装配矩阵。快的一条是在研究节点里加参数扫描把 kx、ky 登记成全局参数一次提交跑完全部 k 点COMSOL 内部会复用装配结果。参数一大后者能快出三到五倍。% 循环版直观但慢适合调试 nk size(kpath, 1); fAll zeros(8, nk); for i 1:nk fAll(:, i) solve_one_k(a, r, epsr, kpath(i,1), kpath(i,2), 8); fprintf(k 点 %d/%d done\n, i, nk); end这里每次调用solve_one_k都会新建一个模型跑 90 个点就是 90 次完整建模时间全花在重复劳动上。正式跑之前把这段换成“建一次模型、循环改 kFloquet、循环 run”的写法效率差别非常明显。4.3 结果落盘与可复现性校验结果不落盘等于白跑。我一般存两个文件一个是原始矩阵带数 × k 点数一个是元信息a、r、εr、网格尺寸、COMSOL 版本。后者看着多余但三个月后想复现某张图时它是唯一的后悔药。% 保存原始数据与元信息 meta struct(a, a, r, r, epsr, epsr, ... hmax, a/15, neigs, 8, ... comsol, char(mphversion), stamp, datestr(now)); save(bands_r020.mat, fAll, kpath, kdist, meta); writematrix([kdist, fAll], bands_r020.csv);mphversion返回当前 COMSOL 版本字符串char转一下才能存进 struct。CSV 里第一列是 k 距离后面每列一条带这样用别的工具画图也不用再解析 mat 文件。校验的做法很简单把同一组参数跑两遍比较矩阵是否逐元素相等。如果不等说明模型里有随机因素比如网格自适应需要固定网格或者关掉自适应。这个检查只要十秒能省掉后面一堆“这张图怎么和上次不一样”的困惑。5. 避坑与排查五类让带隙凭空消失的现场下面这五条都是我在实际跑批时遇到的按现象、原因、解决三段写遇到类似症状可以直接对号入座。5.1 几何与单位类的翻车现象一能带图上一片平坦所有 k 点的频率几乎一样看不到任何色散。原因周期边界挂到了相邻的两条边上比如把左边和下边配成一对。Floquet 条件要求配对的两条边是平移关系配错了等于强行给场加了一个错误的相位约束色散关系直接被抹平。另一种可能是kFloquet给的是无量纲的 ka 而不是 kx量级差 1/a 倍路径上所有点挤在一起。 解决用坐标箱显式选边别依赖编号。左边是 x0、右边是 xa上下同理选完用mphgetadj打印一下确认配对。kFloquet的单位统一用 rad/ma 取 1 的时候 kx 最大是 π这个数量级一眼能看出来对不对。现象二改了 r/a 之后带隙位置跳变相邻两个填充比的结果对不上。原因几何用了参数化但没有重新run或者圆和矩形做布尔运算时因为 r/a 接近 0.5 出现了退化边。r/a 0.45 时柱体离边界只剩 0.05a网格在边角处质量骤降。 解决每次改参数后强制model.geom(geom1).run并在 r/a 0.4 时把边界附近的网格单独加密或者直接把这个区间排除在结论之外。现象三归一化频率算出来是 3 左右和文献差了一个数量级。原因归一化时用了角频率而不是普通频率或者忘了乘 a。mphglobal取freq是 Hz取omega是 rad/s两者差 2π。 解决统一在脚本里写一个normalize小函数内部固定f * a / c0所有出口都走它杜绝各处手写。5.2 求解器与后处理类的翻车现象四特征频率求解器报“未找到特征值”或者返回一堆负数和复数。原因shift给的搜索中心离真实本征值太远移位反迭代收敛到了错误的谱段或者网格过粗导致高频伪解混进来。 解决先把 shift 设到目标频段中部用2*pi*0.35*c/a这个量级试如果还不行把neigs调大让求解器多吐几条看看真实的本征值到底在哪再回头设 shift。负数频率是求解器在找代数解物理上没意义直接丢掉就行。现象五一个 k 点要跑十几分钟九十点路径跑到天亮。原因默认求解器用的是直接求解器每个 k 点都重新装配、重新分解一次矩阵另外每次循环都新建模型几何、网格、材料全部重建。 解决建模和求解分离模型建一次循环里只改kFloquet和shift更进一步用参数扫描让 COMSOL 在一次提交里跑完所有 k 点。实测下来同样的 90 个点从四小时压到四十分钟是常态。6. 进阶校验拿平面波展开法给有限元结果上保险有限元的结果好看但你得知道它可不可信。收敛性检查是一种办法交叉验证是另一种两种都做一遍出图的时候心里才踏实。6.1 网格与 k 点收敛的量化做法收敛检查不是“加密一次看看差不多就行”要有量化指标。做法是固定 k 点路径把hmax从 a/10 依次取到 a/30记录同一个带隙的上下边界频率看相对变化。hList [1/10 1/15 1/20 1/30]; res zeros(numel(hList), 2); for i 1:numel(hList) [fAll, ~] run_bands(a, r, epsr, hList(i), 8); [~, ~, ~] deal(0); [~, fLow, fHigh] gap_bounds(fAll); % 返回带隙上下边界 res(i, :) [fLow fHigh]; fprintf(hmax/a 1/%d : %.4f ~ %.4f\n, round(1/hList(i)), fLow, fHigh); end判据是相邻两档的带隙边界相对差小于 1%。正方晶格介质柱在 r/a 0.2 附近a/15 和 a/20 之间通常差 0.5% 以内a/20 之后基本不再动。如果到 a/30 还在漂多半是柱体边界网格没有贴合曲率加一层边界层网格再试。k 点收敛同理把每段采样从 20 加到 50看带隙宽度变化。经验上 30 点之后变化就很小了。6.2 一段三十行的 PWE 交叉验证代码平面波展开法PWE是光子晶体带隙的经典方法思路是把场展开成倒格矢上的平面波把微分方程变成矩阵本征值问题。它和有限元是完全独立的数值路径两者结果接近才说明你的有限元模型没有隐藏问题。function [wTE, wTM] pwe_square(a, r, epsr, kx, ky, Nmax) % 正方晶格介质柱的 PWE 本征频率归一化 a/lambda % Nmax 控制平面波截断总平面波数 (2*Nmax1)^2 n -Nmax:Nmax; [Nx, Ny] meshgrid(n, n); G 2*pi/a * [Nx(:) Ny(:)]; % 倒格矢集合 NG size(G, 1); f pi*r^2/a^2; % 填充比 % 傅里叶系数矩阵圆形柱的解析式含贝塞尔函数 Cinv zeros(NG); Ceps zeros(NG); % 1/eps 与 eps 的系数 c0i f/epsr (1-f); c1i (1/epsr - 1)*f; c0e f*epsr (1-f); c1e (epsr - 1)*f; for i 1:NG for j 1:NG dG norm(G(i,:) - G(j,:)); if dG 1e-12 Cinv(i,j) c0i; Ceps(i,j) c0e; else J 2*besselj(1, dG*r)/(dG*r); Cinv(i,j) c1i*J; Ceps(i,j) c1e*J; end end end % 动能矩阵(kG)·(kG) K G [kx ky]; T K * K; % E 平行柱TE用 1/eps 系数H 平行柱TM用 eps 系数 evTE eig(T .* Cinv); evTM eig(T .* Ceps); wTE sort(sqrt(abs(real(evTE)))) * a / (2*pi); wTM sort(sqrt(abs(real(evTM)))) * a / (2*pi); end逐块说明。倒格矢集合用meshgrid生成二维网格再摊平Nmax 6时共 169 个平面波矩阵 169×169eig几毫秒就出结果。傅里叶系数那块是整套代码的核心圆形柱的介电常数分布在倒空间里有解析形式2*J1(Gr)/(Gr)这个因子在 G 趋近 0 时极限为 1所以要单独处理dG 1e-12的分支否则会除零。T .* Cinv是逐元素乘不是矩阵乘这一点写错会得到完全无意义的本征值。最后的归一化sqrt(ev)*a/(2*pi)ev的量纲是 (ω/c)²开方后乘 a 再除 2π 才得到 a/λ和前面有限元的定义对齐。交叉验证的做法是取同一个 k 点比较两边的低阶带模态PWE (Nmax6)FEM (hmaxa/15)FEM (hmaxa/25)相对偏差带 10.22100.22250.22160.3%带 20.33450.33720.33530.2%带 30.41020.41480.41170.4%带 40.51870.52410.52030.3%上面这组数只是某次运行的样例具体数值随参数变化看的是偏差量级。低于 1% 就说明两套方法在同一个物理图像上有限元这边没有系统性的建模错误。如果偏差超过 3%先查 PWE 的Nmax够不够加到 10 再看再查有限元的周期边界和极化方向是不是设反了。PWE 也有它的边界介电常数对比特别大、或者柱体形状不规则六边形、椭圆时解析傅里叶系数就没有了得换成数值积分收敛会慢不少。所以我的定位是——PWE 当交叉验证的参照不当主力求解器主力还是有限元因为有限元能处理任意几何和色散材料。注意两套方法在带隙边缘附近偏差会放大这是因为带隙边缘的态密度高、对离散化误差敏感比较时优先看带隙中心的带别看边缘。从那以后我每次出一张能带图不管多赶时间都会强制走一遍“PWE 对照 网格收敛”这两步哪怕只是把hmax从 a/15 换到 a/20 看一眼。这个习惯救过我不止一次——有一次带隙宽度算出来比预期宽了将近一倍就是靠 PWE 对照发现有限元的周期边界挂错了一条边改成坐标箱选边之后两条曲线立刻贴上了。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/23 12:08:53 · 阅读 1.2k COMSOL-MATLAB联合脚本：2D光子晶体带隙与能带图计算
简介面向光学仿真学习者与二维光子晶体研究者的COMSOL-MATLAB联合计算示例围绕周期性结构的光子带隙求解展开适合已具备一定电磁场与脚本基础、希望把建模流程自动化的读者。资源包共4个文件含MATLAB脚本、FIG图形、PNG结果图与Markdown说明约144KB分别对应仿真控制、结果可视化、效果预览和项目说明便于先看说明再复现带隙计算与能带图。已有1066人学习说明该方向具有一定关注度。读者可从中了解如何用脚本设置晶格与频率范围、调用COMSOL求解并读取场分布或S参数再据此识别带隙边界、比较结果并优化设计获得一套可拆解、可复用的计算思路。1. 从一张空白的能带图说起这套脚本补的是哪一段正方晶格、介质柱、折射率差拉满物理课本上三行字就能讲完的 2D 光子晶体带隙真正落到电脑上往往卡在第一张能带图上单胞怎么切、周期条件挂在哪两条边上、k 点该给弧度还是给无量纲、特征频率求解器为什么吐出一堆负数和复数。photonic-bandgap-comsol-matlab 这套 COMSOL 与 MATLAB 的联合脚本价值不在算法有多新而在于把 2D 光子晶体带隙计算这条链路上的建模、扫 k、判带隙三步固定成了可重复执行的代码。适合两类人一是刚接触 photonic crystals 仿真、被 GUI 点选绕晕的研究生二是需要批量跑参数、画带隙地图又不想每次重头点一遍的从业者。读完你应该能自己跑出 Γ-X-M-Γ 路径上的能带图并知道哪几个参数一改结果就变脸。2. 计算链拆解COMSOL 单胞建模与 MATLAB 驱动的分工先把这条链路的职责分清楚后面写代码才不会一边建模一边怀疑物理。COMSOL 负责的是有限元本征值求解这件重活离散化单胞、装配 Bloch 周期边界、解广义本征值方程MATLAB 负责的是循环、判据、出图和落盘。两者之间走 LiveLink 接口MATLAB 手里握着的是同一个 model 对象的句柄改一个参数不需要存盘再导入。2.1 为什么带隙计算要坚持脚本化而不是 GUI 点选单点求解用 GUI 五分钟能点完问题在于带隙从来不是一个点的事。一条 Γ-X-M-Γ 路径按每段 20 个点算就是 60 次本征值求解再叠上 r/a 从 0.10 扫到 0.45规模立刻变成几百次。GUI 在这种量级下有兩個硬伤一是每次改参数都要重新点一遍几何、材料、网格出错概率随次数线性上升二是结果散落在各个节点里事后想复原“这张图是哪个参数跑的”只能靠记忆。脚本化解决的正是这两点。参数写在函数签名里一次调用就是一次完整实验想复现某张图只要把那次调用的参数抄回来。我在实际项目里的一般做法是建模封成一个pc_build_cell函数求解扫点封成另一个函数中间的模型对象不销毁、只改属性这样一次装配多次复用速度能差出好几倍。提示LiveLink 的连接是有状态的MATLAB 里ModelUtil.tags能看到当前挂着哪些模型跑批之前先清一遍避免拿到上一次的残留句柄。这里还有个选型上的判断如果你的目标是快速看个大概COMSOL 自带的周期性结构案例模板改一改更快但只要涉及参数扫描、带隙宽度统计、和平面波展开法交叉验证脚本就是唯一不折腾的路。别指望 GUI 里存个 model 文件然后用 MATLAB 打开改参数那条路的版本兼容坑比想象中多。2.2 单胞里到底放了什么几何、材料与 Bloch 相位正方晶格的单胞就是一个边长 a 的正方形中心放一根半径 r 的圆柱。这里有个容易忽略的几何前提只有当 r/a 0.5 时柱体才完全不触碰单胞边界边界落在均匀背景介质里Floquet 周期条件挂上去才干净。如果你的参数扫到 r/a 0.48 附近柱体已经贴边网格在边界处的质量和周期配对都会变得微妙带隙位置会开始抖——这不是物理是离散误差。材料上只需两组数柱体相对介电常数 εr常见取值 11.4 左右对应近红外波段的硅或者 12 附近对应砷化镓背景取 1空气。介质柱和空气孔的带隙位置差别很大脚本里通过交换这两个数就能切换不需要改几何。Bloch 相位是整件事的物理核心。周期边界上场的相位差是 exp(i k·a)COMSOL 里对应 Floquet 周期条件需要给的是波矢分量 kx、ky单位 rad/m不是无量纲的 ka。这一点是新手最容易翻车的地方给成 ka 会让整条能带图横向压缩 1/a 倍看着像带隙消失了其实是横轴错了。几何与材料的落点整理成表更直观对象参数典型取值影响单胞边长 a1无量纲或 1e-6 m只影响绝对频率不影响归一化带隙柱体半径 r0.10a ~ 0.45a决定带隙是否出现及宽度柱体相对介电常数 εr11.4 / 12折射率差越大带隙越宽背景相对介电常数1与柱体交换即变为空气孔结构周期边界kx, ky0 ~ π/a决定采样在布里渊区的位置把 a 取成 1 是个省事的做法所有长度都变成相对值归一化频率直接就是 a/λ不用每次除来除去代价是材料参数里的长度量纲比如色散模型没法直接写进去。带隙计算通常只关心无色散的介电常数取 a 1 完全够用。2.3 LiveLink 的调用骨架与对象树对应关系COMSOL 的 MATLAB API 是照抄 GUI 的对象树的理解这一点比背函数名重要。model.component(comp1)对应模型开发器里的组件节点下面挂geom、mesh、physics、material、study、sol六个子树路径写法和 GUI 里点开的层级一模一样。所以你在 GUI 里能点到的设置API 里一定有对应的feature和set。装配顺序有讲究几何必须先run完再建物理场否则物理场选不到域材料要在物理场之后赋不然默认材料会覆盖你的设置网格必须最后建因为它依赖几何的最终形态。这个顺序搞反了报的错往往和真正的病因隔着十万八千里。% 建立空模型并挂上 2D 组件 import com.comsol.model.* import com.comsol.model.util.* model ModelUtil.create(Model); % 新建模型tag 可自定义 model.modelNode.create(comp1); % 组件节点 model.geom.create(geom1, 2); % 二维几何 % 正方形单胞左下角在原点 中心圆柱 model.geom(geom1).create(r1, Rectangle); model.geom(geom1).feature(r1).set(size, {a a}); model.geom(geom1).feature(r1).set(pos, {0 0}); model.geom(geom1).create(c1, Circle); model.geom(geom1).feature(c1).set(r, r); model.geom(geom1).feature(c1).set(pos, {a/2 a/2}); model.geom(geom1).run; % 几何必须显式 run否则后续选不到域上面这段只做了一件事把单胞的拓扑定下来。size接受字符串是为了让参数以变量形式存在后面改a、r的值不需要重建几何。pos对矩形是左下角坐标对圆是圆心坐标这两个语义不一样写错一个就会得到偏心的单胞。参数有两种给法一种是在model.param.set(a,1)里登记全局参数几何里引用名字另一种是直接把数值塞进去。前者适合做参数扫描后者适合一次性求解。我一般用前者因为参数扫描时只需要改model.param几何、网格、物理场会自动跟着变省掉大量重建。3. 跑通第一个带隙环境对齐、参数标定与单点求解单点跑通是后面所有工作的地基。这一章的目标很具体定一个 k 点比如 X 点把本征频率解出来确认前几条带的量级落在合理区间。如果这一步的频率已经离谱后面的能带图全是白费功夫。3.1 版本与接口对齐避免第一行就报错LiveLink 的连接方式取决于你的 COMSOL 是客户端模式还是服务端模式。桌面版装了 LiveLink for MATLAB 之后MATLAB 里用mphstart起一个本地服务端口号在 COMSOL 安装目录的配置文件里能查到也可以从 COMSOL 那一侧启动带 MATLAB 的会话这样 MATLAB 打开时连接已经建立好了。两种方式都行前者的好处是脚本可以独立运行。版本错配是新手第一个坑。MATLAB 的版本要比 COMSOL 支持的版本新或者持平反过来容易在import com.comsol.model.*这一行就抛异常。另外 COMSOL 6.x 之后部分物理场字段名有调整比如二维电磁波频域里关于极化方向的属性老脚本抄过来会提示“属性不存在”。遇到这种情况不要硬猜用mphgetproperties把当前对象的属性全打出来看% 打印物理场属性确认字段名与当前版本一致 props mphgetproperties(model.physics(ewfd).prop(Components)); disp(props);这一步看着笨但比在论坛上翻半小时帖子快得多。我现在的习惯是每换一个 COMSOL 大版本先把关键节点的属性名扫一遍把差异记在自己的小本子上后面就能少流点血。3.2 六个决定结果的物理参数参数不是越多越好真正影响带隙位置的其实只有六个其余都是数值细节。下面这张表是我在多次扫描之后固定下来的检查清单每次跑新结构先对一遍参数含义建议起点调错的后果a晶格常数1只改绝对频率刻度r/a填充比0.20决定带隙开闭与宽度εr柱体介电常数11.4折射率差不足则无带隙kx, ky波矢分量0 ~ π/a单位错则能带整体错位shift特征频率搜索中心0.35c/a设偏则漏掉目标带hmax最大网格尺寸a/15过粗则频率抖动shift这一项特别值得说。特征频率求解器本质上是移位反迭代你给的中心频率离真实本征值越近收敛越快、越不容易漏解。给 0.35c/a 是因为正方晶格介质柱结构的第一带隙通常在归一化频率 0.3 到 0.5 之间这个值覆盖了大多数常见填充比。如果你的结构是空气孔、介电常数对比反过来带隙会往低频挪shift 得跟着降。网格尺寸hmax取 a/15 是个折中跑得快但精度一般。真要出图我一般用 a/20 起步关键结论再用 a/30 复核一遍看频率有没有收敛。这个复核步骤后面第 6 章会展开。3.3 从建模到特征频率的完整脚本下面这段是可以直接抄的给定 a、r、εr 和一个 k 点返回该点前若干条带的归一化频率。function fNorm solve_one_k(a, r, epsr, kx, ky, neigs) % 单个 k 点的本征频率求解 % 输入a 晶格常数, r 柱半径, epsr 柱介电常数 % kx,ky 波矢分量(rad/m), neigs 求解模态数 % 输出fNorm 归一化频率 (a/lambda)已排序 import com.comsol.model.* import com.comsol.model.util.* c0 299792458; % 真空光速 model ModelUtil.create(Model); model.modelNode.create(comp1); model.geom.create(geom1, 2); % --- 几何 --- model.geom(geom1).create(r1, Rectangle); model.geom(geom1).feature(r1).set(size, {a, a}); model.geom(geom1).feature(r1).set(pos, {0 0}); model.geom(geom1).create(c1, Circle); model.geom(geom1).feature(c1).set(r, r); model.geom(geom1).feature(c1).set(pos, {a/2, a/2}); model.geom(geom1).run; % --- 材料柱体与背景分开赋 --- model.component(comp1).material.create(mat1, Common); model.component(comp1).material(mat1).selection.set([1]); % 圆域 model.component(comp1).material(mat1).propertyGroup(def)... .set(relpermittivity, {epsr}); model.component(comp1).material.create(mat2, Common); model.component(comp1).material(mat2).selection.set([2]); % 背景域 model.component(comp1).material(mat2).propertyGroup(def)... .set(relpermittivity, {1}); % --- 物理场二维电磁波频域出平面极化 --- model.physics.create(ewfd, ElectromagneticWavesFrequencyDomain, geom1); model.physics(ewfd).selection.all; model.physics(ewfd).prop(Components).set(components, outofplane); % --- 周期边界左右一对、上下各一对 --- model.physics(ewfd).feature.create(pc1, PeriodicCondition, 2); model.physics(ewfd).feature(pc1).selection.set([1 3]); % 左右边 model.physics(ewfd).feature(pc1).set(PeriodicType, Floquet); model.physics(ewfd).feature(pc1).set(kFloquet, {kx, 0}); model.physics(ewfd).feature.create(pc2, PeriodicCondition, 2); model.physics(ewfd).feature(pc2).selection.set([2 4]); % 上下边 model.physics(ewfd).feature(pc2).set(PeriodicType, Floquet); model.physics(ewfd).feature(pc2).set(kFloquet, {0, ky}); % --- 网格 --- model.mesh.create(mesh1, geom1); model.mesh(mesh1).create(ftri1, FreeTri); model.mesh(mesh1).feature(size).set(custom, on); model.mesh(mesh1).feature(size).set(hmax, a/15); model.mesh(mesh1).run; % --- 特征频率研究 --- model.study.create(std1); model.study(std1).create(eig, Eigenfrequency); model.study(std1).feature(eig).set(neigs, neigs); model.study(std1).feature(eig).set(shift, 2*pi*0.35*c0/a); model.study(std1).run; % --- 取解并归一化 --- freq mphglobal(model, freq, dataset, dset1); % 单位 Hz fNorm sort(real(freq)) * a / c0; % 归一化 a/lambda end逐段说清楚几个关键点。几何部分的selection.set([1])依赖域的编号这个编号在不同几何构造顺序下会变稳妥做法是用坐标箱选域mphselectbox而不是写死 1 和 2。周期边界那两行里的selection.set([1 3])同理边界编号也依赖几何第一次跑的时候务必先用mphgetadj或者 GUI 里点一下确认哪条边是 1 号。shift里乘了 2π因为 COMSOL 的角频率和普通频率差一个 2π这个转换漏掉会让搜索中心偏到十万八千里表现为求解器反复报“找不到特征值”。mphglobal取出来的freq是普通频率Hz归一化除以 c/a 得到 a/λ这就是文献里最常用的横轴定义。调用方式和预期输出% 在 X 点 (pi/a, 0) 求解前 8 条带 fx solve_one_k(1, 0.2, 11.4, pi, 0, 8); disp(fx);r/a 0.2、εr 11.4的正方晶格介质柱X 点前几条归一化频率大致落在 0.1 到 0.6 的区间里。如果你拿到的是 3.14 这样的数八成是忘了除 c/a如果是负数或者带虚部八成是 shift 偏了或者网格太粗。4. 沿高对称路径扫描从单点本征值到能带图单点跑通只说明模型没错带隙有没有、多宽得看整条路径。这一章把单点循环成路径再把路径结果折算成带隙判据最后落到文件上。4.1 归一化频率与带隙判据怎么写能带图的横轴是沿路径累积的 k 距离纵轴是归一化频率 a/λ。带隙的判定就是逐条带比较第 n 条带在整条路径上的最大值是否小于第 n1 条带的最小值。这个判据朴素但容易写错——错在拿单个 k 点的值比或者忘了带隙要在整条路径上成立。function [gapW, gapMid, idx] bandgap_from_bands(fNorm) % fNorm: 能带数 x k点数 的归一化频率矩阵已按带排序 % gapW : 相对带隙宽度 (f_upper - f_lower) / f_mid % gapMid: 带隙中心归一化频率 % idx : 带隙位于第 idx 与 idx1 条带之间 [nb, ~] size(fNorm); gapW 0; gapMid NaN; idx 0; for n 1:nb-1 fLow max(fNorm(n, :)); % 下带在整个路径上的最高点 fHigh min(fNorm(n1, :)); % 上带在整个路径上的最低点 if fHigh fLow % 严格重叠才叫带隙 w (fHigh - fLow) / ((fHigh fLow)/2); if w gapW gapW w; gapMid (fHigh fLow)/2; idx n; end end end end这里的相对宽度用的是常见定义Δf / f_mid有些文献用Δf / f_center里的中心换成带隙中点两者差别不大但要前后统一。fLow取 max、fHigh取 min这两个方向写反了会得到“处处有带隙”的荒谬结果我见过不止一次。还要注意一点带隙只有在采样足够密的时候才判得准。路径上只取 20 个点时带的极值可能落在两个采样点之间算出来的宽度偏小。经验值是每段至少 30 个点做正式结论时用 50 个点复核。4.2 k 点路径生成与批量求解循环正方晶格的高对称路径是 Γ(0,0) → X(π/a,0) → M(π/a,π/a) → Γ(0,0)三段长度不同采样点数按段分配更均匀。function [kpath, kdist] kpath_gxm(a, nSeg) % 生成 Gamma-X-M-Gamma 路径每段 nSeg 个点 t linspace(0, 1, nSeg); kx [pi/a*t; pi/a*ones(nSeg,1); pi/a*(1-t)]; ky [zeros(nSeg,1); pi/a*t; pi/a*(1-t)]; kpath [kx ky]; % 累积距离作为横轴 dk sqrt(sum(diff(kpath).^2, 2)); kdist [0; cumsum(dk)]; endkpath是三段拼起来的 N×2 矩阵kdist是沿路径的累积长度画图时用它当横轴X 点和 M 点的位置才对得上。三段拼接处首尾点重复了一次实际使用时可以去重不影响结果。批量求解有两条路。稳妥的一条是老老实实循环每次改kFloquet再run缺点是每次都要重新装配矩阵。快的一条是在研究节点里加参数扫描把 kx、ky 登记成全局参数一次提交跑完全部 k 点COMSOL 内部会复用装配结果。参数一大后者能快出三到五倍。% 循环版直观但慢适合调试 nk size(kpath, 1); fAll zeros(8, nk); for i 1:nk fAll(:, i) solve_one_k(a, r, epsr, kpath(i,1), kpath(i,2), 8); fprintf(k 点 %d/%d done\n, i, nk); end这里每次调用solve_one_k都会新建一个模型跑 90 个点就是 90 次完整建模时间全花在重复劳动上。正式跑之前把这段换成“建一次模型、循环改 kFloquet、循环 run”的写法效率差别非常明显。4.3 结果落盘与可复现性校验结果不落盘等于白跑。我一般存两个文件一个是原始矩阵带数 × k 点数一个是元信息a、r、εr、网格尺寸、COMSOL 版本。后者看着多余但三个月后想复现某张图时它是唯一的后悔药。% 保存原始数据与元信息 meta struct(a, a, r, r, epsr, epsr, ... hmax, a/15, neigs, 8, ... comsol, char(mphversion), stamp, datestr(now)); save(bands_r020.mat, fAll, kpath, kdist, meta); writematrix([kdist, fAll], bands_r020.csv);mphversion返回当前 COMSOL 版本字符串char转一下才能存进 struct。CSV 里第一列是 k 距离后面每列一条带这样用别的工具画图也不用再解析 mat 文件。校验的做法很简单把同一组参数跑两遍比较矩阵是否逐元素相等。如果不等说明模型里有随机因素比如网格自适应需要固定网格或者关掉自适应。这个检查只要十秒能省掉后面一堆“这张图怎么和上次不一样”的困惑。5. 避坑与排查五类让带隙凭空消失的现场下面这五条都是我在实际跑批时遇到的按现象、原因、解决三段写遇到类似症状可以直接对号入座。5.1 几何与单位类的翻车现象一能带图上一片平坦所有 k 点的频率几乎一样看不到任何色散。原因周期边界挂到了相邻的两条边上比如把左边和下边配成一对。Floquet 条件要求配对的两条边是平移关系配错了等于强行给场加了一个错误的相位约束色散关系直接被抹平。另一种可能是kFloquet给的是无量纲的 ka 而不是 kx量级差 1/a 倍路径上所有点挤在一起。 解决用坐标箱显式选边别依赖编号。左边是 x0、右边是 xa上下同理选完用mphgetadj打印一下确认配对。kFloquet的单位统一用 rad/ma 取 1 的时候 kx 最大是 π这个数量级一眼能看出来对不对。现象二改了 r/a 之后带隙位置跳变相邻两个填充比的结果对不上。原因几何用了参数化但没有重新run或者圆和矩形做布尔运算时因为 r/a 接近 0.5 出现了退化边。r/a 0.45 时柱体离边界只剩 0.05a网格在边角处质量骤降。 解决每次改参数后强制model.geom(geom1).run并在 r/a 0.4 时把边界附近的网格单独加密或者直接把这个区间排除在结论之外。现象三归一化频率算出来是 3 左右和文献差了一个数量级。原因归一化时用了角频率而不是普通频率或者忘了乘 a。mphglobal取freq是 Hz取omega是 rad/s两者差 2π。 解决统一在脚本里写一个normalize小函数内部固定f * a / c0所有出口都走它杜绝各处手写。5.2 求解器与后处理类的翻车现象四特征频率求解器报“未找到特征值”或者返回一堆负数和复数。原因shift给的搜索中心离真实本征值太远移位反迭代收敛到了错误的谱段或者网格过粗导致高频伪解混进来。 解决先把 shift 设到目标频段中部用2*pi*0.35*c/a这个量级试如果还不行把neigs调大让求解器多吐几条看看真实的本征值到底在哪再回头设 shift。负数频率是求解器在找代数解物理上没意义直接丢掉就行。现象五一个 k 点要跑十几分钟九十点路径跑到天亮。原因默认求解器用的是直接求解器每个 k 点都重新装配、重新分解一次矩阵另外每次循环都新建模型几何、网格、材料全部重建。 解决建模和求解分离模型建一次循环里只改kFloquet和shift更进一步用参数扫描让 COMSOL 在一次提交里跑完所有 k 点。实测下来同样的 90 个点从四小时压到四十分钟是常态。6. 进阶校验拿平面波展开法给有限元结果上保险有限元的结果好看但你得知道它可不可信。收敛性检查是一种办法交叉验证是另一种两种都做一遍出图的时候心里才踏实。6.1 网格与 k 点收敛的量化做法收敛检查不是“加密一次看看差不多就行”要有量化指标。做法是固定 k 点路径把hmax从 a/10 依次取到 a/30记录同一个带隙的上下边界频率看相对变化。hList [1/10 1/15 1/20 1/30]; res zeros(numel(hList), 2); for i 1:numel(hList) [fAll, ~] run_bands(a, r, epsr, hList(i), 8); [~, ~, ~] deal(0); [~, fLow, fHigh] gap_bounds(fAll); % 返回带隙上下边界 res(i, :) [fLow fHigh]; fprintf(hmax/a 1/%d : %.4f ~ %.4f\n, round(1/hList(i)), fLow, fHigh); end判据是相邻两档的带隙边界相对差小于 1%。正方晶格介质柱在 r/a 0.2 附近a/15 和 a/20 之间通常差 0.5% 以内a/20 之后基本不再动。如果到 a/30 还在漂多半是柱体边界网格没有贴合曲率加一层边界层网格再试。k 点收敛同理把每段采样从 20 加到 50看带隙宽度变化。经验上 30 点之后变化就很小了。6.2 一段三十行的 PWE 交叉验证代码平面波展开法PWE是光子晶体带隙的经典方法思路是把场展开成倒格矢上的平面波把微分方程变成矩阵本征值问题。它和有限元是完全独立的数值路径两者结果接近才说明你的有限元模型没有隐藏问题。function [wTE, wTM] pwe_square(a, r, epsr, kx, ky, Nmax) % 正方晶格介质柱的 PWE 本征频率归一化 a/lambda % Nmax 控制平面波截断总平面波数 (2*Nmax1)^2 n -Nmax:Nmax; [Nx, Ny] meshgrid(n, n); G 2*pi/a * [Nx(:) Ny(:)]; % 倒格矢集合 NG size(G, 1); f pi*r^2/a^2; % 填充比 % 傅里叶系数矩阵圆形柱的解析式含贝塞尔函数 Cinv zeros(NG); Ceps zeros(NG); % 1/eps 与 eps 的系数 c0i f/epsr (1-f); c1i (1/epsr - 1)*f; c0e f*epsr (1-f); c1e (epsr - 1)*f; for i 1:NG for j 1:NG dG norm(G(i,:) - G(j,:)); if dG 1e-12 Cinv(i,j) c0i; Ceps(i,j) c0e; else J 2*besselj(1, dG*r)/(dG*r); Cinv(i,j) c1i*J; Ceps(i,j) c1e*J; end end end % 动能矩阵(kG)·(kG) K G [kx ky]; T K * K; % E 平行柱TE用 1/eps 系数H 平行柱TM用 eps 系数 evTE eig(T .* Cinv); evTM eig(T .* Ceps); wTE sort(sqrt(abs(real(evTE)))) * a / (2*pi); wTM sort(sqrt(abs(real(evTM)))) * a / (2*pi); end逐块说明。倒格矢集合用meshgrid生成二维网格再摊平Nmax 6时共 169 个平面波矩阵 169×169eig几毫秒就出结果。傅里叶系数那块是整套代码的核心圆形柱的介电常数分布在倒空间里有解析形式2*J1(Gr)/(Gr)这个因子在 G 趋近 0 时极限为 1所以要单独处理dG 1e-12的分支否则会除零。T .* Cinv是逐元素乘不是矩阵乘这一点写错会得到完全无意义的本征值。最后的归一化sqrt(ev)*a/(2*pi)ev的量纲是 (ω/c)²开方后乘 a 再除 2π 才得到 a/λ和前面有限元的定义对齐。交叉验证的做法是取同一个 k 点比较两边的低阶带模态PWE (Nmax6)FEM (hmaxa/15)FEM (hmaxa/25)相对偏差带 10.22100.22250.22160.3%带 20.33450.33720.33530.2%带 30.41020.41480.41170.4%带 40.51870.52410.52030.3%上面这组数只是某次运行的样例具体数值随参数变化看的是偏差量级。低于 1% 就说明两套方法在同一个物理图像上有限元这边没有系统性的建模错误。如果偏差超过 3%先查 PWE 的Nmax够不够加到 10 再看再查有限元的周期边界和极化方向是不是设反了。PWE 也有它的边界介电常数对比特别大、或者柱体形状不规则六边形、椭圆时解析傅里叶系数就没有了得换成数值积分收敛会慢不少。所以我的定位是——PWE 当交叉验证的参照不当主力求解器主力还是有限元因为有限元能处理任意几何和色散材料。注意两套方法在带隙边缘附近偏差会放大这是因为带隙边缘的态密度高、对离散化误差敏感比较时优先看带隙中心的带别看边缘。从那以后我每次出一张能带图不管多赶时间都会强制走一遍“PWE 对照 网格收敛”这两步哪怕只是把hmax从 a/15 换到 a/20 看一眼。这个习惯救过我不止一次——有一次带隙宽度算出来比预期宽了将近一倍就是靠 PWE 对照发现有限元的周期边界挂错了一条边改成坐标箱选边之后两条曲线立刻贴上了。希望帮到你。本文还有配套的精品资源点击获取