基于COMSOL的光子晶体光纤SPR传感器与三芯分束器复现

先说这次复现的结论两张图一张是SPR光子晶体光纤传感器另一张是三芯光子晶体光纤分束器全部在COMSOL里从零建模跑通最终和原文的损耗谱、模场分布基本对齐。整个过程踩了不少坑尤其是“论文只给折射率不给材料模型”这种经典问题。这篇就把我从课题设计、参数还原、几何建模、边界模式分析到网格和扫参调优的完整过程记录下来顺便把图左到图右的对应关系讲清楚给后面做光子晶体光纤复现的人一个可参考的路线。1. 复现前先看懂这两个仿真到底在算什么1.1 SPR传感器用损耗峰的位置变化去感知外界折射率SPR表面等离子体共振传感器这名字听着高大上拆开看就是光在光纤纤芯里传播时如果包层里镀了一层金属膜金、银或者ITO这类材料那么特定波长的光会激发金属-介质界面上的表面等离子体激元SPP。SPP是一种沿着金属表面传播的电磁波模式它有一个特点能量被强烈束缚在金属表面而且它的传播常数对接触金属表面的介质折射率极其敏感。当外界被测液体的折射率改变时模式匹配条件就变了体现在仿真结果上就是损耗谱上吸收峰的波长位置发生移动。我在复现时选的传感器结构是这样的光子晶体光纤PCF的包层空气孔里某个或某几个孔的内壁镀上金属层再在金属层里面填充待测分析物。仿真的核心就是计算这种结构的模场损耗谱也就是基模在不同波长下的限制损耗Confinement Loss。损耗峰出现的位置和强度依赖于金属层的厚度、分析物折射率、空气孔直径和孔间距也就是节距。论文里通常直接给出损耗峰波长在折射率1.33到1.38范围内的线性响应你要复现的就是这条曲线。这里有个关键认知损耗峰不是随便设个材料就冒出来的它需要满足“相位匹配”也就是纤芯基模的有效折射率实部要接近SPP模式的有效折射率实部差值越小耦合越强损耗峰越明显。所以调整几何参数时本质就是在调两个模式的色散曲线交点这个理解后面排查问题非常有用。1.2 三芯分束器把能量从一个芯分到两个芯三芯分束器比SPR结构“温和”一些不涉及金属和SPP纯粹靠光纤芯间的倏逝场耦合。它一般由三个纤芯组成通常在PCF里通过“移除某些空气孔”形成一个或三个高折射率区域光从中间芯入射由于相邻芯之间的模式耦合光会周期性地从一个芯转移到另一个芯。这里要复现的核心结果是在不同传播距离下三个芯中的功率分配情况以及对应的“周期”。耦合长度Lc由超模之间的传播常数差决定Lc π / (β₁ - β₂)。其中β₁、β₂是三芯结构中对称模和反对称模的传播常数。所以仿真中需要同时找到前几个高阶模式的色散关系而不能只看基模。在COMSOL里处理这种周期性耦合结构最常用的做法是直接做截面上的特征值分析边界模式分析得到每个模式的模场分布和有效折射率然后通过模式叠加和传播相位来推算不同纵向位置的功率分布。更直接的做法是建一个三维结构跑电磁波频域求解但三维的耗散太大论文复现阶段一般先用二维截面分析确定Lc和模式特性再适当做三维验证。1.3 图左原文和图右复现的差异点这篇复现的图左是从论文里截出来的原图图右是我在COMSOL里重新画的。两者之间的差异有几个容易忽略的细节原图的空气孔边界不是理想圆形在镀膜局部会有过渡层而仿真模型通常把膜看作均匀同心环这就导致损耗峰位置有微小偏移。原图的颜色映射色标不同特别是模场分布图COMSOL的默认色标和论文可能不一致所以不能只看颜色要看强度的空间分布形态。三芯分束器原图可能显示了传播方向的功率变化三维演示而复现时我用了截面模式分析加后处理合成呈现的是一段传播距离内的二维分布。对照图片时不要追求像素级一致重点是验证模式形态和曲线趋势的物理一致性。2. 关键参数还原与建模思路2.1 如何从论文里反向推断材料参数论文正文不一定把所有参数都列全尤其是材料折射率的公式。PCF的基底材料通常是二氧化硅熔融石英它在1550nm附近的折射率约1.444但是不同文献用的Sellmeier系数略有差异。这里我直接用最常用的Sellmeier方程n²(λ) 1 A₁λ²/(λ² - B₁) A₂λ²/(λ² - B₂) A₃λ²/(λ² - B₃)熔融石英的常见系数λ单位μm为A₁0.6961663B₁0.0684043A₂0.4079426B₂0.1162414A₃0.8974794B₃9.896161如果你要复现的论文没写Sellmeier那大概率就是用的这个标准值。但有一种情况需要警惕论文可能在某个特定波长只给单一折射率值比如“n1.45 at 1550nm”那你就直接用这个值比Sellmeier更符合原文因为原文的计算结果很可能是用这个固定值算出来的。金属材料的折射率同样如此。SPR传感器常用金Au或银Ag它们的介电常数可以用Drude模型拟合ε(ω) ε∞ - ωp² / (ω² iγω)金的典型参数在可见光和近红外区间可以取ε∞1.0ωp1.37e16 rad/sγ1.21e14 rad/s。银则是ε∞1.0ωp1.39e16 rad/sγ3.3e13 rad/s。如果你在COMSOL里用插值表也可以直接导入Johnson和Christy的实验数据这个精度更高但注意不同版本的COMSOL材料库自带的金数据可能偏旧导致损耗峰出现几百纳米偏移。我的建议是先用Drude模型快速计算确认损耗峰的大致位置再换成实验数据精调。因为实验数据插值在损耗峰附近的斜率更准确能显著提高和原文曲线的重合度。2.2 几何建模空气孔、镀膜层和分析物流道的画法COMSOL里画PCF截面推荐用“参数化曲线旋转”的方式而不是手动一个个画圆。因为PCF的周期结构就是三角形晶格空气孔的位置完全可以由两个基矢推导出来。比如六角晶格的空气孔中心坐标可以这样生成x(i, j) i * Λ j * Λ / 2 y(i, j) j * √3 * Λ / 2其中Λ是节距i、j是整数索引。在COMSOL的全局参数里定义好Λ、空气孔直径d、镀膜厚度t后在几何节点中用循环阵列Array功能把圆一次性分布到位。这样后面如果要扫参比如扫描d/Λ从0.3到0.6只需要改两个参数几何自动更新。镀膜层有两种画法圆环法添加两个同心圆一个半径r一个半径r-t用“差集”形成圆环。这样会形成独立的域可以单独赋金属材料。指定边界法直接把金属层定义为边界条件该边界上有“过渡边界条件”或“薄层”。COMSOL里用边界代替薄层可以大幅减少网格量。但如果金属层较厚超过趋肤深度最好还是用域建模否则精度会有明显损失。SPR传感器通常需要在分析物流道内壁镀膜那么流道本身是个空气孔流道里填充的是水n≈1.33~1.38。建模时要把这个大孔和周围的PCF周期孔区分开。我的做法是先画PCF的周期空气孔再用一个大圆覆盖其中一两个孔的位置做布尔操作生成流道区域再在流道边界上生成镀膜环。注意顺序不能反否则布尔操作之后膜层边界会支离破碎。2.3 材料分配折射率、吸收和频散的正确写法COMSOL的材料节点里可以手工输入折射率的实部和虚部。对SPR传感器来说复折射率的虚部决定了损耗所以必须设置成复数值。金属用复介电常数在波动光学模块中可以用“折射率”选项直接填n和k消光系数也可以填“介电常数”的实部虚部两者本质等价但填折射率更直观。流体区域的分析物折射率通常只给实部例如1.33那就不填虚部。但如果你想更精确水里其实在近红外也有很微弱的吸收量级在1e-5以下对损耗谱峰值位置影响不大我试过几乎没变化。二氧化硅在近红外的吸收系数很小但计算限制损耗时要把材料吸收包含进去吗严格来说论文里的损耗如果是“限制损耗”是指模式能量从纤芯泄漏到包层导致的损耗和材料吸收无关。很多论文仿真时就是不勾选材料吸收只让COMSOL计算出实模态中的泄漏损耗。在边界模式分析中如果材料折射率虚部设置为0得到的有效折射率虚部就代表泄漏损耗。如果设置成非零虚部那就混入了材料吸收损耗峰会更高更宽可能和论文对不上。所以复现时先不加材料损耗除非明确说明加了。3. 实操过程与核心环节实现3.1 搭建模型从二维几何到完美匹配层我采用的完整流程如下新建模型选择“二维”空间维度“电磁波频域ewfd”物理场接口研究选择“特征值”频率选在目标波段如1550nm。在全局参数节点中定义lam0扫描中心波长lambda_p、dPCF的节距和孔直径t_m镀膜厚度n_analyte分析物折射率建几何。先建一个半径约10倍于结构尺寸的圆形域作为计算区域包一圈完美匹配层PML。PML厚度通常取工作波长的1~2倍例如1550nm时取2~4μm。PML外边界用散射边界条件这样泄漏模的辐射不会被反射回来干扰模式。用“数组”功能布置空气孔阵列。选择三角形晶格范围覆盖整个计算域然后把落在流道区域内的孔删掉。添加流道大孔和镀膜环利用布尔操作生成独立的金属域。分配材料。基底用熔融石英金属域用金的复折射率或插值材料流道内用分析物折射率PML域用内置的PML材料。物理场设置上边界模式分析的求解类型选“特征值”COMSOL会求解麦克斯韦方程组在横向截面上的模式解输出结果就是有效折射率neff和模场分布。注意要指定“模式搜索基准值”比如1.44附近否则求解器可能找半天找到一堆远离目标的模式。3.2 扫参设计一个参数一个参数地动论文复现最忌讳一次性把所有参数都改了。我的顺序是先固定结构参数扫光纤归一化频率或波长确认损耗峰位置后再扫金属厚度最后扫分析物折射率。以SPR传感器为例要在COMSOL里实现“损耗谱随波长变化”有两种方式直接做频域扫描边界模式分析在波长为变量时也可以作为频域扫描处理每次求解一个特征值问题。COMSOL的“辅助扫描”功能可以让参数lam0跨越一段波长范围每步都做一次特征值求解。后处理绘制损耗曲线损耗α的定义为α(dB/m) 20 * log10(e) * (2π / λ) * Im(neff) * 1e3。也就是把有效折射率的虚部转换成传播损耗。这个公式必须自己写COMSOL不会自动给你算“限制损耗”它只会给你复有效折射率。我实际操作时把lam0扫了100步波长从600nm到1100nm或看你需要的波段每一步都是边界模式分析求解时间取决于网格粗细。如果网格太粗损耗峰会变得很钝太细单步求解时间可能超过30秒200步就要近两小时。建议先用粗网格最大单元尺寸200nm跑一遍快速定位峰的位置再用细网格最大单元尺寸100nm把峰值附近加密复算。三芯分束器也类似但要额外关注高阶模。因为三芯结构中至少存在三个超模分别对应三根芯的对称和反对称组合。在用特征值求解时要设置“所需模式数量”为6~8个保证把相关模式全部找出来。然后查看每个模式的neff实部和虚部从中找成对模式计算传播常数差。3.3 数据导出与画图如何让输出和论文长得像模式分析完成后要输出损耗谱曲线需要用“一维绘图组”绘制α随波长的变化。把上述损耗公式写到“替换表达式”里横轴选择lambda纵轴选择表达式即可得到和图左对应的谱线。模场分布图则是二维绘图组画|E|的分布并使用镜像或对称复制来展示完整结构。三芯分束器若想展示传播效果可以用“功率”端子或后处理中的“电磁波能量密度”分层染色。实际输出时注意单位论文里经常用mV⁻¹或归一化强度自己仿真得到的是V/m需要做归一化处理。常见的一个坑是COMSOL默认显示的是全局坐标系下的电场分量但PCF模式往往是混合偏振HE、EH模式等所以最好分别绘制Ex、Ey和|E|选能反映模式对称性的那个分量。比如第一次做时基模可能看起来是x偏振的类LP01模但论文用的可能是y偏振。所以画图前先把每种偏振的模式都求一遍再决定用哪个模式跟踪损耗峰。4. 常见问题与排查技巧实录4.1 模式找不到或模式数不对调整搜索基准值特征值求解时最常遇到的问题是“找不到模式”或“找到一堆无关模式”。原因通常是搜索基准值与真实模式相差太远。COMSOL的“特征值搜索基准值”默认是0这对光学模式完全没用要把基准值设置为基底折射率比如1.44附近。如果还是找不到还有一个技巧把待求模式数量增大到20个左右然后从求解结果中手动筛选物理合理的模式。因为泄漏模的虚部通常比实部小几个数量级但比数值伪模大你可以通过模场分布图来判断。那些场分布杂乱无规律、能量不在纤芯位置的就是伪模直接忽略。4.2 损耗峰不出现或者位置偏差大这是SPR复现里最常见的问题。损耗峰找不到我总结有几种原因金属层厚度偏离最佳值。金膜太薄时SPP模式损耗太大模式非常宽太厚时表面等离子体波无法与纤芯模有效耦合。常见有效厚度在30-50nm之间最佳厚度取决于材料和波长。分析物流道离纤芯太远。SPP模式的倏逝场随金属距离指数衰减如果流道离芯太远耦合强度指数下降损耗峰就消失了。材料折射率设置成了纯实数。如果金属材料没有虚部就不会有损耗峰的概念所有模式都是无损耗的。网格太粗导致SPP模式无法分辨。SPP模式的特征尺寸很小在金属表面附近的场变化剧烈网格至少要小于表面波波长的一半。排查时先扫一个很宽的折射率范围比如分析物从1.32扫到1.40看有没有共振谷出现。如果没有就换金属厚度再看。如果还没有就把几何简化成平面结构先验证SPR是否存在再从结构上排查。4.3 三芯分束器的能量不对称网格和几何对称性三芯分束器的一个典型问题按理说三个芯完全对称但仿真结果里功率分配却不对称。排除物理原因后大概率是网格划分破环了几何对称性。COMSOL在自由三角形网格划分时如果不勾选“保持对称”有时候网格在两个镜像区域不是严格对称的导致数值上有微小差异。解决方式是在网格节点中手动设置对称平面的镜像网格操作或者在几何建模时就用半截面建模再加对称性边界条件。对三芯分束器完全对称的结构用全截面建模也没问题但网格划分时可以使用“复制”网格操作把平面一侧的网格镜像复制到另一侧彻底消除不对称。4.4 扫参时模式追踪错乱扫参扫描待测折射率或厚度时COMSOL会自动按特征值大小排序但这个排序不一定与物理模式顺序一致。比如在某个波长上基模的neff最大在另一个波长上某个高阶超模的neff超过了基模那么模式编号就变化了。这时要使用模式追踪功能。在特征值研究的求解器配置里可以通过添加“模式追踪”功能根据前一步的模式重叠积分来判断当前步该取哪个模式。或者更土的办法扫参完成后把模式的模场分布导出来看哪一个模式在每一步都保持相同的场分布特征。4.5 收敛性差和计算时间过长边界模式分析通常不需要迭代求解内迭代但网格太密时内存占用会飙高。这里有个经验值二维截面网格数量在10万到50万之间一般16GB内存的电脑都能承受。超过50万就要考虑用PML的“映射网格”了因为PML区域占了大头用映射网格可以大幅度减少网格数。如果扫参确实太多建议先把波长步长变粗如20nm一步确认峰位后再在峰值附近10-20nm范围内加密。比一开始就用高精度扫全场要节省好几倍时间。5. 关于复现效果的几点心得最后说一个我实际复现多次之后才体会到的点论文里的仿真图尤其是SPR损耗谱往往是在特定条件下算出来的但论文未必把网格数、PML厚度、模式搜索范围这些细节写清楚所以不可能做到主峰和旁瓣完全一模一样。目标应该是复现出和论文相同的物理趋势和关键特征峰的位置不差超过2%峰的形状半高宽、相对高度整体相似这就够了。实际操作里我复现SPR传感器时第一轮扫出来的峰值比论文偏红了约40nm。当时第一反应是金属折射率模型不对换成实验数据后只改善了10nm左右。后来仔细查发现是PML厚度太薄用的0.5倍波长导致泄漏模边界有微弱反射影响了模式虚部的精度。把PML加到1.5倍波长后峰值位置立刻对上了。这种细节在论文里绝对不会写但就是复现成败的关键。所以排查问题不要总盯着材料参数边界条件这类“隐形设置”往往才是罪魁祸首。如果你在复现时遇到图上对不上但又找不出原因建议按这个顺序排查几何尺寸对应原文参数表→ 材料折射率 → PML参数和边界条件 → 网格细化 → 模式追踪设置。这个顺序可以帮你少走很多弯路。