先说这次复现的结论：两张图，一张是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.6961663，B₁=0.0684043
- A₂=0.4079426，B₂=0.1162414
- A₃=0.8974794，B₃=9.896161
如果你要复现的论文没写Sellmeier，那大概率就是用的这个标准值。但有一种情况需要警惕：论文可能在某个特定波长只给单一折射率值，比如“n=1.45 at 1550nm”，那你就直接用这个值，比Sellmeier更符合原文，因为原文的计算结果很可能是用这个固定值算出来的。
金属材料的折射率同样如此。SPR传感器常用金（Au）或银（Ag），它们的介电常数可以用Drude模型拟合：
ε(ω) = ε∞ - ωp² / (ω² + iγω)
金的典型参数在可见光和近红外区间可以取ε∞=1.0，ωp=1.37e16 rad/s，γ=1.21e14 rad/s。银则是ε∞=1.0，ωp=1.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、d：PCF的节距和孔直径
- 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参数和边界条件 → 网格细化 → 模式追踪设置。这个顺序可以帮你少走很多弯路。
