★ FEATURED ARTICLE 简介将COMSOL Multiphysics与MATLAB结合计算二维光子晶体的光子带隙是光子晶体、集成光学与电磁仿真方向研究者和高校师生常用的联合仿真任务。光子带隙意味着特定频率范围内的电磁波无法在周期性晶体中传播直接影响光通信、光存储和光学微腔器件设计压缩包把完整流程打包为可直接运行的脚本方案。包内共4个文件、大小144KB类型涵盖m、fig、png、mdMATLAB脚本负责调用COMSOL建模、参数扫描与后处理fig保存可交互的能带结构图形png提供结果对比参考md梳理从几何设置、PML边界条件到带隙提取的思路。已有1065人学习浏览适合初次接触该接口的读者对照练习。下载后可自由调整晶格常数、单元形状与边界条件自动完成电磁场求解、能带曲线绘制和带隙边界识别结合附带对比图可验证仿真准确性为后续光子晶体滤波器、微腔或光波导器件设计提供可扩展的仿真基础和可修改的MATLAB代码。1. 用 COMSOL-MATLAB 脚本计算 2D 光子晶体的带隙把重复点击变成一条命令用 COMSOL-MATLAB 脚本计算 2D 光子晶体的带隙这件事难的不是单个 k 点上的本征值求解而是把几十个 k 点串起来的过程每改一次波矢就要重设 Floquet 周期条件、重新求解、重新取数。刚开始我手动在 GUI 里点一个 Γ-X-M-Γ 路径下来要四个小时中间还会因为漏改一个周期边界条件得到一条断成几截的假带隙。后来把整个流程写成脚本几何参数、波矢路径、特征值搜索全部塞进 MATLAB 循环同样一条路径十分钟跑完改半径也只是改一个变量。这篇笔记按这个思路展开读者是准备做光子晶体设计、或者想把带隙计算从 GUI 里解放出来的工程师。你会从这里看到怎么搭模型、怎么扫波矢、怎么取频率以及那些让结果悄悄变错的细节。2. 2D光子晶体带隙计算的物理边界为什么波矢扫描绕不开布里渊区带隙不是某个 k 点上没有模式而是整个第一布里渊区里都没有模式。这句话听起来简单但几乎所有翻车现场都源于对它理解不够。2D 光子晶体是介电常数在平面内周期性变化的平板或薄膜结构周期结构里的本征模式受 Bloch 定理约束本征频率是约化波矢 k 的函数。你算出来的所有能带本质上是波矢 k → 本征频率 ω的映射关系。只要有一个 k 点有模式落进某个频率区间那个区间就不能叫带隙。这就是为什么必须扫波矢而且不是扫几个点就算完必须扫密一点、扫对路径。2.1 2D光子晶体带隙在解什么本征频率对波矢的依赖周期介质中电场或磁场可以写成 Bloch 形式E(r) u_k(r)e^{ik·r}其中 u_k(r) 与晶格周期相同k 是约化波矢。把这个形式代入 Maxwell 方程组会得到一个以频率为特征值的算子方程。物理上等价于把无限大的周期结构约化到单个晶胞在晶胞边界上施加 Floquet 周期条件把 k 作为边界条件的相位因子然后对每个 k 求解本征频率。COMSOL 的电磁波频域接口里不需要手写 Bloch 定理它的周期性条件节点选择 Floquet 类型后会在相对边界的自由度之间建立一个复相位关系这个相位关系就是 k 的体现。2D 光子晶体还有一个便利平面内传播的波可以按偏振解耦。电场沿面外方向时称为 TM 模磁场沿面外方向时称为 TE 模。TM 模在数学上退化为一个关于 E_z 的标量方程TE 模则是关于 H_z 的标量方程。这意味着你可以用同一个几何、同一个网格分别求解两套本征问题得到两套带结构。很多论文里只给一条带隙通常默认指 TM 带隙但你复现别人结果时一定要确认对方算的是哪一套。COMSOL 的二维电磁波接口默认求解的是 Ez 分量对应 TM如果你需要 TE 带隙需要用面外磁场接口或者把材料属性做对应的对偶替换。数值上的关键是本征值问题。对每个 k有限元离散后得到一个广义本征值方程解出一组本征频率。带隙计算要求你在一个频率范围内把所有本征频率都找出来而不是只找前几个。对于常见的正方晶格空气孔结构晶格常数 a 和孔半径 R 决定了能带位置高频段模式密度大模式间距小漏掉一个模式就可能在图上造成一条假带隙。我一般每个 k 点算 8 到 12 个本征频率然后对着实部从小到大排序保证画出来的带结构是一条条连续的曲线。2.2 为什么锁定 COMSOL-MATLAB三种计算方案的取舍做 2D 光子晶体带隙业界常见有三条路线。第一条是平面波展开工具通常被叫做 MPB它把介电常数展开成 Fourier 级数在倒空间求解矩阵本征值问题速度非常快几条能带几十秒就能扫完适合快速扫描参数空间。但它的几何必须规则圆孔、方柱、三角排列没问题一旦出现渐变折射率、多层包层、非对称缺陷建模就比较痛苦。第二条是商用 FDTD 软件时域有限差分方法几何自由度很高但算带隙需要多次傅里叶变换和模式分析流程上比直接求本征值绕远而且细网格下时间步长受 Courant 条件限制跑一个完整路径不轻松。第三条就是 COMSOL 加 MATLAB 的 LiveLink 路线有限元方法几何任意各向异性、复介电常数、缺陷态都可以直接在物理场接口里设置。方案数值方法几何自由度带隙计算便利度脚本控制程度平面波展开MPB倒空间本征值规则柱/孔高命令行直接扫高商用 FDTD时域差分任意低需要额外后处理中COMSOL MATLAB有限元任意支持各向异性/损耗中但可完全脚本化高三条路线里COMSOL-MATLAB 不是最快的但它是最不容易被几何卡死的。我的习惯是如果结构是标准的圆孔或圆柱阵列先用平面波展开工具粗扫参数确定大概的带隙位置确定好半径和介电常数之后再用 COMSOL 做验证。但如果你要做的是正方形孔、椭圆柱、双层包层或者带一个缺陷腔那就直接上 COMSOL因为平面波展开在这些场景下要么建模很别扭要么根本做不了。脚本化的另一个理由是重复性设计过程中至少要扫描孔径、周期、材料三个变量每个变量十几个取值手动点击的量是几百次每次都容易出低级错误。一旦把建模和扫描写进 MATLAB 脚本改一个变量就等于重跑一遍全流程出错的概率大幅下降。3. COMSOL-MATLAB脚本搭建几何、Floquet周期条件与特征值研究这一章直接给可运行的建模代码以最常见的正方晶格空气孔为例背景是高折射率介质周期为 a 的正方形单胞中心挖一个空气圆孔。整个建模流程和 GUI 操作的顺序一一对应但所有尺寸都被提前写成变量方便后面扫参数。注意 COMSOL with MATLAB 的 API 在不同版本里属性名有细微差别如果某一行报错先右键对应节点导出 .m 脚本对照一下属性名再改这是排查脚本问题最快的办法。3.1 用 ModelUtil 建单胞几何先把参数集中在大顶几何是正方形单胞减去中心圆孔。由于后面要扫 R/a 和 a 本身所有长度都写成变量。下面这段代码创建 2D 几何并执行布尔减操作。a 0.6; % 晶格常数单位微米。这个值决定带隙落在哪个波长段 R 0.3 * a; % 空气孔半径先取 0.3a带隙宽度对 R/a 非常敏感 model ModelUtil.create(Model); model.hist.disable; % 关闭模型历史记录批量循环时能明显省内存 model.component.create(comp1, true); model.comp1.geom.create(geom1, 2); % 注意把几何长度单位设为微米后面所有波矢单位都以它为基准 model.comp1.geom(geom1).set(unit, um); model.comp1.geom(geom1).create(r1, Rectangle); model.comp1.geom(geom1).feature(r1).set(size, [a a]); model.comp1.geom(geom1).feature(r1).set(pos, [0 0]); model.comp1.geom(geom1).create(c1, Circle); model.comp1.geom(geom1).feature(c1).set(r, R); model.comp1.geom(geom1).feature(c1).set(pos, [a/2 a/2]); % 正方形域减去圆孔得到背景介质域 model.comp1.geom(geom1).create(dif1, Difference); model.comp1.geom(geom1).feature(dif1).selection(input).set(r1); model.comp1.geom(geom1).feature(dif1).selection(input2).set(c1); model.comp1.geom(geom1).run;model.hist.disable这一行很多人会忽略但它直接影响后面扫描循环的性能。模型历史记录会把每一步建模操作都存下来单次建模无所谓循环里每改一次波矢就多一层记录跑二十个 k 点之后 MATLAB 会明显变慢占用几个 GB 内存。set(unit,um)的作用是统一长度单位COMSOL 默认以米为基准如果你不在几何节点里改单位那么后面周期条件里的波矢单位就要用 rad/m配合微米尺寸的几何时数值会差 10^6 倍这是带结构断线的常见原因。孔心放在(a/2, a/2)使单胞在 x 和 y 方向对称边界条件在这个几何上最容易收敛。3.2 材料与电磁波物理场把 Floquet 周期条件变成可循环变量材料定义分两个域背景介质用相对介电常数 3.5 的平方也就是折射率 3.5孔内是空气相对介电常数为 1。物理场用电磁波频域接口然后添加 Floquet 类型周期性条件。% 背景介质高折射率材料介电常数 n^2 model.comp1.material.create(mat1, Common); model.comp1.material(mat1).propertyGroup(def).set(relpermittivity, 3.5^2); model.comp1.material(mat1).selection.set(1); % 域编号以实际几何为准 % 空气孔 model.comp1.material.create(mat2, Common); model.comp1.material(mat2).propertyGroup(def).set(relpermittivity, 1); model.comp1.material(mat2).selection.set(2); % 域编号以实际几何为准 % 电磁波频域物理场 model.comp1.physics.create(ewfd, ElectromagneticWaves, Frequency); model.comp1.physics(ewfd).create(pbc1, PeriodicCondition, 1); model.comp1.physics(ewfd).feature(pbc1).selection.set([1 2 3 4]); % 单胞四条外边界 model.comp1.physics(ewfd).feature(pbc1).set(PeriodicType, Floquet); % 先用占位波矢后面循环里由脚本更新 model.comp1.physics(ewfd).feature(pbc1).set(k0, kx, 0); model.comp1.physics(ewfd).feature(pbc1).set(k0, ky, 0);材料域编号这里要提醒一句几何执行布尔减后背景介质通常被编号为 1空气孔为 2但如果你把几何节点顺序换过编号可能反过来。第一次跑脚本前可以在 COMSOL 桌面里打开模型看一眼选择列表或者在 MATLAB 里执行model.comp1.geom(geom1).getAdjacency确认域编号。物理场的 Floquet 条件里k0是周期边界的波矢分量属性名。不同版本的 LiveLink 里它可能叫k_floquet或者显示为k_x所以如果 set 不进去不要硬猜右键模型树的周期性条件节点导出脚本看那一行真实的名字。这个步骤花两分钟能省掉后面一整天的排查时间。波矢单位要格外小心。COMSOL 的 Floquet 条件期望的波矢是物理波矢单位由几何单位决定。几何用微米时第一布里渊区边界 X 点是 (π/a, 0) rad/μm而不是无量纲的 (π, 0)。很多带结构断线问题都是因为直接填了无量纲波矢导致相对边界相位差过大有限元网格分辨不了快速振荡。3.3 网格控制与研究步骤特征值数量不是越多越好网格密度直接影响带隙位置尤其是空气孔边界的曲率。圆孔用自由三角形网格时最大单元尺寸必须能分辨孔边缘否则等效介电常数分布被磨平带隙会偏窄甚至消失。研究步骤用特征值求解器设置要算的模式数量和搜索基准。% 自由三角形网格孔边界加密 model.comp1.mesh.create(mesh1); model.comp1.mesh(mesh1).create(ftri1, FreeTri); model.comp1.mesh(mesh1).feature(ftri1).set(size, custom); model.comp1.mesh(mesh1).feature(ftri1).set(hmax, R/3); % 最大尺寸取孔半径的 1/3 model.comp1.mesh(mesh1).run; % 特征值研究期望 8 个模式搜索基准放在归一化频率 0.45 附近 model.study.create(std1); model.study(std1).create(eig, Eigenvalue); model.study(std1).feature(eig).set(neigs, 8); model.study(std1).feature(eig).set(eigenvaluesearch, around, shift, 0.45); model.sol.create(sol1); model.sol(sol1).study(std1); model.sol(sol1).create(st1, StudyStep); model.sol(sol1).runAll;hmax取 R/3 是我常用的起点。对于 0.3a 的孔半径这个尺寸大约把孔圆周分成 20 段左右能保证前三阶能带的频率误差控制在 1% 量级。如果后续发现带隙边界对网格敏感再加密到 R/6 做对比。特征值研究的两个参数neigs是希望求解器返回的特征值个数取 8 是因为低频段的带隙通常出现在第 2 到第 5 条带之间算 8 个能覆盖到带隙上边缘shift是搜索基准告诉求解器在哪个频率附近找特征值。这个值不要取太小否则会漏掉低频模式也不要离目标带隙太远否则求解器可能花很多迭代去找一堆无关的高频模式。对 a0.6 μm、折射率 3.5 的空气孔结构归一化频率 a/λ 在 0.2 到 1.0 之间有完整的带结构搜索基准 0.45 大致落在带隙区域。4. 波矢扫描主循环Γ-X-M-Γ路径参数化与特征频率批量提取模型建好后剩下的事全部在 MATLAB 里循环。核心思路是对每个 k 点更新 Floquet 条件运行求解器提取特征频率存进矩阵最后画图。这一章先讲清楚路径和横轴映射再给完整代码。4.1 把布里渊区路径写成坐标序列横轴与波矢的单位约定正方晶格的第一布里渊区是一个正方形高对称点是 Γ(0,0)、X(π/a,0)、M(π/a,π/a)。能带极值通常出现在布里渊区边界和高对称点上所以带隙计算的标准路径是 Γ → X → M → Γ。横轴不能直接用波矢长度因为 Γ 到 X 的长度和 X 到 M 的长度不一样画图时要把三段折线映射成一个连续的分段线性坐标。我习惯把横轴坐标定义为Γ 为 0X 为 1M 为 2回到 Γ 为 3。每一段内部均匀采样。这样做的好处是带结构图的横轴没有物理单位纯粹是路径顺序所有文献里的带结构图也都是这个惯例。需要说明的是X 到 M 段上 kx 保持 π/a 不变ky 从 0 变到 π/aM 回到 Γ 段则是 kx 和 ky 同时从 π/a 线性减到 0。采样密度上每段 15 个点对大多数结构足够带隙边界附近如果出现区域很窄的带隙可以加密到 25 点再确认。4.2 主循环改波矢、求解、取数、绘图下面这段是完整主循环。提取特征频率的逻辑封装成一个独立函数方便复用。% 每一段采样点数15 个点对常见结构足够分辨带边界 n1 15; n2 15; n3 15; % 路径坐标单位 rad/um因为几何单位是 um kx [linspace(0, pi/a, n1), linspace(pi/a, pi/a, n2), linspace(pi/a, 0, n3)]; ky [linspace(0, 0, n1), linspace(0, pi/a, n2), linspace(pi/a, 0, n3)]; xs [linspace(0, 1, n1), linspace(1, 2, n2), linspace(2, 3, n3)]; band zeros(length(kx), 8); % 每一行存一个 k 点的 8 个特征频率 for i 1:length(kx) % 更新 Floquet 周期条件的波矢分量 model.comp1.physics(ewfd).feature(pbc1).set(k0, kx, kx(i)); model.comp1.physics(ewfd).feature(pbc1).set(k0, ky, ky(i)); % 清掉上一 k 点保存的解避免内存累积 model.sol(sol1).clearSolution; model.sol(sol1).runAll; % 提取特征频率并按大小排序 f extractEigFreq(model); band(i, :) sort(real(f)); fprintf(k[%.3f %.3f] 前三条频率: %s\n, kx(i), ky(i), mat2str(band(i,1:3), 3)); end % 画归一化带结构图纵轴 a/lambda a*f/c normFreq a * 1e-6 * band / 3e8; % a 转米f 是 Hzc 是 m/s plot(xs, normFreq, .-, MarkerSize, 6); xlabel(波矢路径 (Γ-X-M-Γ)); ylabel(归一化频率 a/λ); grid on;function f extractEigFreq(model) % 用全局计算节点导出本征频率 try model.result.numerical(ev1).clear; catch model.result.numerical.create(ev1, EvalGlobal); end model.result.numerical(ev1).set(expr, ewfd.freq); model.result.numerical(ev1).set(solnum, all); data model.result.numerical(ev1).getData(); if iscell(data) data data{1}; end f real(data(:, 1)); end这段代码里两个地方需要重点解释。第一个是clearSolution如果不加这一行COMSOL 会把之前 k 点的解保留在求解器里跑十几轮后内存占用会持续增加求解时间也会变长。第二个是排序。特征值求解器返回的模式顺序不是按频率严格递增的尤其是存在简并模式时顺序可能跳变。如果你不 sort带结构图上会出现同一 k 点两条能带交叉后又跳回的假象sort 之后曲线才连续。ewfd.freq是 COMSOL 内置变量表示特征频率单位 Hz取实部是因为有损材料或数值误差可能带来微小的虚部带隙计算只看实部。4.3 循环内存与失败恢复跑长循环时还有一个容易忽略的问题某个 k 点可能因为网格质量或搜索基准不合适导致求解失败。这时主循环会直接抛出异常之前所有 k 点的结果全部丢失。我的习惯是把主循环包在 try-catch 里失败时保留已算出的band矩阵单独记录失败的索引然后继续跑下一个点。failed []; for i 1:length(kx) try model.comp1.physics(ewfd).feature(pbc1).set(k0, kx, kx(i)); model.comp1.physics(ewfd).feature(pbc1).set(k0, ky, ky(i)); model.sol(sol1).clearSolution; model.sol(sol1).runAll; f extractEigFreq(model); band(i, :) sort(real(f)); catch failed(end1) i; %#okSAGROW fprintf(第 %d 个 k 点求解失败: k[%.3f %.3f]\n, i, kx(i), ky(i)); end endclearSolution和runAll的搭配是这套脚本的稳定关键。clearSolution并不是重建求解器只是清掉上一轮的数值结果下一次runAll会重新开始迭代。这样每次迭代的内存占用基本恒定不会像滚雪球一样越跑越慢。失败点处理上我通常先把循环跑完再看failed里有没有点如果有单独对那些 k 点调大neigs或者改shift再算一次而不是从头全部重跑。5. 带隙计算避坑指南单位、漏根、网格与模式辨识带隙计算这个领域结果出错往往不是某一个环节的大错误而是好几个小问题叠加。下面五条是实际跑脚本最容易踩的坑每条按现象、原因、解决三个层面说清楚。5.1 Floquet波矢的单位与周期边界为什么带结构会断成几截现象带结构图上的能带不是连续曲线而是每隔几个 k 点就断一下或者出现一条带突然跳到很高的频率再跳回来。原因周期条件里的波矢单位与几何单位不匹配。几何默认用米时波矢单位是 rad/m如果几何尺寸用微米但波矢直接填了 π/a 这个量级实际相位差会扩大 10^6 倍边界条件变成完全随机的相位求解出的模式自然乱七八糟。解决建模开始时就在几何节点里显式设定单位model.comp1.geom(geom1).set(unit,um)波矢按 rad/μm 填写。还有一个隐蔽版本几何单位是 nm波矢按 μm 填同样会断线。统一原则是几何用什么长度单位波矢里 a 就用什么单位不要混用。5.2 漏根与假带隙neigs和搜索基准怎么配合现象某个 k 点算出的 8 条带里第 3 条和第 4 条之间出现一个莫名其妙的巨大间隙而且这个间隙只在部分 k 点出现整条带结构不连续。原因特征值求解器漏掉了某个本征模式。漏根最典型的触发条件是搜索基准离真实特征值太远或者neigs设置得太小求解器提前收手。解决把neigs从 8 提高到 12 到 16重新跑一次对比低频的 6 条带是否完全重合。如果提高数量后多出了之前没见过的模式说明之前的带隙是假带隙。另一种验证方式是取一个高对称点如 M 点单独用更密的网格和更大的neigs算一次确认模式数不随neigs变化。5.3 网格收敛带隙边界随网格漂移别当玄学现象同一个结构网格粗的时候带隙宽度是 0.08加密一档变成 0.11再加密变成 0.115。这个差距远超预期。原因圆孔边界的曲率是光子晶体带隙计算里最敏感的部分。太粗的三角形网格把圆弧切成了多边形等效孔径偏小带隙自然偏窄。解决把孔边界最大单元尺寸从 R/3 加密到 R/6 和 R/12分别计算观察带隙边界的变化量。变化量小于 1% 时认为网格收敛。我的经验是 R/3 到 R/6 之间带隙通常还会移动 3% 到 5%R/6 到 R/12 就稳定了。如果你的结构对带隙精度要求很高最终结果用 R/6 或更密网格保存。5.4 TE/TM偏振搞反和文献对不上的常见原因现象带结构图形态和某篇文献里的图很像但带隙位置偏移很大或者该有带隙的频段没有带隙。原因很可能偏振没对齐。COMSOL 的二维电磁波接口默认求解面外电场分量 Ez这对应 TM 偏振。如果你想复现的是 TE 带隙直接用默认物理场是不对的。解决先确认目标文献的带隙是 TE 还是 TM。如果是 TE在 COMSOL 中使用面外磁场变量 Hz也就是创建电磁波接口后把磁场的面外分量设为主变量或在物理场里选择面外磁场矢量模式。更稳妥的做法是对同一模型分别建两个物理场接口一个算 TE 一个算 TM输出两套带结构标注清楚后再和文献比对。5.5 结果提取与排序拿到手的频率为什么是乱序的现象主循环里不排序直接画图能带曲线在同一条带里出现锯齿比如第 2 条带在某几个 k 点突然跳到第 4 条带的位置再跳回来。原因特征值求解器返回的模式顺序由内部算法决定不是按频率从小到大。尤其是存在简并模式的 k 点求解器可能把两个相同频率的模式以任意顺序返回。解决提取频率后必须sort(real(f))而且画图时每列数据代表频率从低到高第几条带不要试图追踪某个模式的身份。如果你需要追踪模式比如分析带边模式场分布要用 COMSOL 的求解器选择某一个特征值而不是数组位置来定位。6. 结果自检与批量筛选把脚本封装成带隙扫描工具拿到一张带结构图先别急着写进报告。我的自检流程有三步第一看 M 点附近是否出现预期中的简并模式。正方晶格空气孔结构在 M 点通常有两重简并如果你画的图上 M 点处两条带正好重合说明网格和边界条件设置基本正确如果有一条带劈开了 1e-4 量级的小口那是网格不对称的数值误差把网格加密即可。第二把带隙上下边缘对应的模式选出来画场分布图确认场确实局域在孔或介质的特定区域且边界处的相位关系和 Floquet 条件一致而不是边界伪影。第三和平面波展开工具或文献的已知结果做一次交叉验证只对比带隙边缘两个频率误差在 2% 以内就能接受。批量筛选参数时把整个脚本封装成一个函数是最值得做的事。函数签名类似[gapWidth, gapMid] calcBandgap(a, R, epsBg, nSeg)内部完成建几何、设材料、扫路径、提取频率、计算带隙宽度。然后写两层循环扫 R/a 和背景折射率几分钟就能拿到一张带隙宽度随占空比变化的曲线。这个过程中要留意每次调用函数都创建一个全新的ModelUtil.create内存不会自动释放建议每个参数组合跑完后model.clear释放模型对象再进入下一轮。最后说一个我自己踩过的教训。有段时间我扫不同孔径的带隙宽度画出来的趋势总跟文献差一截怎么调都不对。后来逐条核对才发现是横轴归一化的问题我把波矢横轴按 2π/a 归一化纵轴归一化频率却按 a/λ两者闭锁条件完全对不上导致带隙宽度被放大了一倍。从那以后我每次跑新结构都会先在 Γ 点和 M 点各算一次把归一化频率和波矢边界手动验算一遍再信任结果。这套 COMSOL-MATLAB 脚本本身不难真正难的是让每一步都符合物理。希望你能用这套思路把自己从重复点击里解放出来把时间花在调结构、想物理上。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
