“非线性材料Z参数扫描”这个需求我第一次听的时候以为是件特别顺手的事光学平台那么成熟COMSOL里把高斯光束、样品、位移台照着画出来设置一个参数扫描拖动样品沿光轴走一段透过率曲线不就出来了吗。真正动手才发现这套想法在软件里几乎走不通——几何一动就要重新剖网格焦点附近网格必须加密到波长级几十个扫描点算下来能熬到天亮而且后处理稍不留神还会得到一条完全没法解释的曲线。后来我把建模思路整个推翻换了“等效扫描”的方式才把问题解决。这篇笔记就是我基于COMSOL做非线性材料Z参数扫描研究和分析方法的一段完整复盘适合正在做非线性光学仿真、Z-scan模拟或者被COMSOL参数化扫描折腾得想摔鼠标的同学参考。先说一句可能颠覆你直觉的结论在COMSOL光学模型里做Z-scan类研究瓶颈往往不在非线性材料本身而在“Z参数”的定义方式。Z参数扫得对后面的扫描、提取、拟合就一路顺畅Z参数映射错了后面每一步都在补前面的坑。1. 先搞清楚Z参数扫描在COMSOL里要模拟什么1.1 一次“把实验台搬进仿真”的翻车经历我第一次做这个题目时建模思路完全是实验复刻。样品是一块厚度很小的非线性材料放在一个可以沿Z轴移动的位移台上激光经过透镜聚焦焦点在样品后面某个位置样品后方放探测器前面加了孔径光阑。我当时把整个系统都塞进了COMSOL样品所在的域写入了非线性折射率并且为了让样品真的“移动”还专门启用了移动网格Deformed Geometry功能想着只要参数扫描改变样品域的Z坐标就能得到透过率随Z坐标变化的曲线。结果非常惨烈。样品从一个位置移到另一个位置时网格自动重构的报错几乎没停过。即使勉强算过几个点单点求解时间长到让人怀疑人生。更麻烦的是只要样品进入强聚焦区域整个非线性方程组的收敛就变得极不稳定算出来的透过率跳来跳去根本没有物理规律。那时候我才意识到Z-scan实验里的“样品移动”和COMSOL仿真里的“扫描参数”并不天然等价直接把实验台上的运动方式搬到几何上是最笨也最容易翻车的一种思路。1.2 Z-scan的物理量与建模输入之间的映射后来我让自己冷静下来先把Z-scan的物理过程拆开看。Z-scan测量的本质是样品沿光轴方向相对焦点移动样品内部的非线性效应引起附加相位或吸收进而改变远场光斑的能量分布最终导致被孔径收集到的光功率发生变化。实验中扫的是样品位置但这只是个操作手段真正关心的是“入射到样品表面的聚焦光场状态随位置的变化”。所以仿真中其实并不一定要让样品真的动。只要能让样品感受到“它相对于焦点在不同位置”物理上就等价于扫描了样品位置。而改变相对焦点位置有两种等价方式移动样品或者移动入射高斯光束的焦点。对薄样品而言只要满足样品厚度远小于瑞利长度的近似条件采用后一种方式会让建模难度降低一个数量级。这也是我在多次尝试后最终采用的做法。1.3 判读参考从曲线形态反查模型是否正确开始建模前还有一件事必须做先明确Z-scan曲线的大致预期形状否则你都不知道仿出来曲线对不对。对于不含非线性吸收的纯克尔介质闭孔Z-scan加小孔光阑归一化透过率曲线特征是非线性折射率n₂大于零时焦点前后依次出现一个谷和一个峰n₂小于零时顺序反过来先出现峰再出现谷。如果是开孔Z-scan全部收集存在双光子吸收时曲线在焦点附近出现一个对称的凹陷。这是最基础的定性判据。我自己一开始就是吃了没做这个功课的亏仿真结果出来一个奇奇怪怪的形状还花了很长时间去调参数最后才发现是模型入口的相位符号设置反了。建议所有人都把这条规律贴在屏幕旁边。2. 光学模型落地时的几何、物理接口与材料设置2.1 几何选型工作平面的作用和求解域到底该画多大我见过不少初学者一上来就建一个三层楼那么大的三维几何好像不这样就不够“真实”。实际上在COMSOL光学模型里几何应该紧扣物理过程来取舍。我做Z-scan模拟用的几何非常简单一段介质通道中间划出样品区域前后是背景介质或空气区域。创建几何建议先从工作平面开始。用工作平面的好处是可以在一个二维平面里画好截面的基本轮廓再通过拉伸、旋转等方式生成三维域COMSOL工作平面的作用不是“画个草图”那么简单它还能帮你准确控制局部坐标系、实体选择和各域的材料指派。这一点在后续筛选取数区域时特别方便。如果你只是在三维环境里凭空画块体后面定义PML层或积分面时经常要反复切割非常麻烦。求解域到底该画多大我的经验是纵向长度必须是激光瑞利长度的几倍以上至少保证从入口到焦点的距离能容纳完整的波前演化横向尺寸则是焦斑尺寸的十倍以上。假设波长为532nm聚焦光斑半径w₀取5μm瑞利长度大约为πw₀²/λ ≈ 148μm那么纵向长度建议取800到1200μm横向半径取40到50μm左右。这个尺寸既能覆盖主要光场能量又不至于让网格数量爆炸。2.2 接口选择全波、波束包络、自定义方程三选一COMSOL处理光学问题有若干条路可走差别极大。直接使用“电磁波频域”的全波形式精度高但是它会要求网格解析光波长532nm的光在几十微米网格量级直接让人绝望除非研究对象本身就在纳米尺度否则我不建议在Z-scan这种宏观光束传播场景中强上全波。更合适的选择是“电磁波波束包络”Beam Envelopes。这个接口专门处理光场包络变化比光波长缓慢的问题它通过预设背景波矢把网格需求从“解析每个波长”放宽到“解析包络变化”在毫米乃至厘米级的传播问题里非常实用。Z-scan里的焦点扫描和瑞利长度范围内的传播用这个接口非常合适。还有一条路线自定义系数型偏微分方程来模拟近轴光束传播也就是把非线性薛定谔方程或菲涅耳衍射积分直接写进PDE模块。这条路的好处是自由度极高能做到轴对称降维计算量比三维波束包络小很多缺点是需要自己推导方程对物理功底要求较高。我的建议是如果只是完成任务优先选波束包络如果想把某个特定非线性效应研究透可以考虑自定义方程。2.3 非线性材料本构COMSOL材料设置里的关键表达式非线性材料在Z-scan语境里一般指折射率或吸收系数与光强相关的材料。最常见的是克尔型非线性折射率写成n(I) n0 n2 * I其中n2是克尔系数单位是m²/WI是局部光强。有非线性吸收时还需加上吸收系数的调制alpha(I) alpha0 beta * ICOMSOL里定义非线性材料时最容易翻车的点是用错强度表达式。仿真中电磁场解出来的电场幅度单位是V/m如果直接把电场模平方当成光强数值会差好几个数量级。正确做法是先定义光强变量I 0.5 * n0 * c_const * epsilon0_const * normE^2其中normE是电场模c_const是真空光速epsilon0_const是真空介电常数。然后把这个变量引用到材料的折射率/吸收系数表达式中。这一步做对了后面的非线性效应才会真正出现在解里。2.4 PML吸收层和网格设置在数值边界处聚焦光场会有旁瓣向外传播直接给固定边界条件会产生反射干扰Z轴上的透过率曲线。我的办法是在域的两端和径向外侧设置完美匹配层PML让透射波在边界处被有效吸收。PML区域不宜太薄建议至少覆盖几个有效波长或数个网格层否则高角度分量吸收不干净。网格方面焦点附近是整条光束能量最集中的地方需要比瑞利长度和焦斑尺度更细的网格。我常用的经验值是在焦斑半径上至少划分5个网格单元如果焦斑w₀为5μm网格尺寸控制在1μm以下。其余区域可以逐渐放松到几微米形成渐变网格。不要一开始就用统一细网格否则计算量会大得毫无意义。3. 让参数扫描跑起来扫描策略与加速手段3.1 固定几何、扫描等效焦点位置我最终采用的参数扫描方案是几何里不放移动部件让整个求解域固定不动把入射高斯光束的焦点位置设成全局参数。这样样品始终在计算域中的同一位置但相对于样品等效焦点一路从样品前方扫到样品后方物理上等价于样品从焦点另一侧移动到另一侧。入射边界上的高斯光束解析式是关键。当焦点离入射面距离为s时入射面上的光斑半径为w_in w0 * sqrt(1 (s / z_R)^2)波前曲率半径为R_in s * (1 (z_R / s)^2)再把这两个量写进出射面作为高斯光束的边界表达式。扫描焦点位置本质上就是在改变入射光场的曲率和光斑半径从而模拟样品在焦点前后不同位置感受到的远场条件。这个方案的好处是几何不用动网格不会重构所有扫描点共用一套网格整个计算稳定性大幅提升。只要样品厚度远小于瑞利长度误差就可以接受如果样品太厚则要考虑样品内部的光束演化那时候这个等效方案会失效。3.2 扫描范围与步长如何设计扫描范围通常取±4倍瑞利长度就够了。Z-scan曲线的有效信息集中在焦点附近的±2倍瑞利长度内再远的地方归一化透过率基本回到1。真正会决定曲线质量的细节在步长规划焦点附近要密集扫描远离焦点的区域可以放宽步长。COMSOL参数化扫描的列表语法支持分区间定义。比如可以先在远离焦点的区域取每0.8倍瑞利长度一个点到±1倍瑞利长度范围内加密到每0.1倍瑞利长度左右焦点中心区域再加密一倍。这一点步长差异看着很小实际对曲线峰谷值的提取影响特别大。峰谷差值ΔT_pv是后面反演材料系数的重要依据如果步长太粗漏掉了真正的峰或谷数据就废了。3.3 两阶段扫描与辅助扫描的使用对于计算量大的模型我习惯分两阶段。第一阶段先跑一个粗步长扫描快速确认曲线的大致形态和粗略峰谷位置第二阶段只对峰谷附近的子范围进行细扫把焦点区域的曲线精细度补上去。这样既能看到全貌又不必把所有点都用细步长从头扫到尾。另外COMSOL的“参数化扫描”支持多参数组合如果除了Z位置还要扫材料非线性系数或孔径尺寸可以把第二参数设为辅助扫描参数。但要注意两参数组合后计算量成倍增加。我建议不要把完整扫描策略写得太贪婪先把主参数Z扫完整确认结果可信后再叠加其它参数这样可以减少很多无效计算。3.4 如果真要让样品移动何时才该用移动网格等效焦点方案并非万能。当样品厚度达到瑞利长度的几倍或者样品内部存在比较明显的折射率梯度/非线性自聚焦效应时入射光在样品内部传播不可忽略必须真实建模样品的位置和内部场演化。这时候就需要在样品域上采用移动网格或变形几何让样品域相对背景域做平移。我的建议是这种情况下不要直接把移动网格和大型参数化扫描组合。可以先给移动网格设定一个较小的变形范围在某几个关键Z位置分别创建单独的几何和网格一个一个求解再把结果拼接成完整曲线。这样虽然人工介入多一些但能避开移动网格在多个位置间反复重构时带来的收敛性和网格质量恶化问题。4. 后处理提炼数据归一化透过率曲线这样得到4.1 定义积分算子求解孔径功率Z-scan闭孔模式最关键的后处理量是孔径收集到的光功率。在COMSOL中可以用“积分算子”在探测器平面上创建一个圆域然后对通过该圆域的坡印廷矢量或光强进行面积分。为了方便切换孔径大小可以把孔径半径设为参数开孔和闭孔只需调整这一参数即可。需要特别提醒的是如果不关心绝对功率一定要跑一个“线性参考样品”作为基准。Z-scan实验的归一化透过率是用有非线性时的透过功率除以线性情况下的透过功率。如果直接用入射总功率做分母样品的线性吸收、反射损耗都会混进曲线导致你看到的不是纯非线性贡献。我在仿真里加了一个非线性开关参数等于0时材料折射率虚部为0、n2为0等于1时启用真实非线性然后在后处理中用“非线性解除以线性解”来得到归一化透过率。4.2 用withsol提取每一组扫描解参数扫描结束后所有结果都在一个数据集里。很多同学直接对默认绘图节点点开发现只有第一个解的曲线就在那儿怀疑自己扫了个寂寞。正确的做法是用withsol函数从扫描数据集中抽取值。比如要在全局绘图里画出孔径功率随Z参数变化的曲线可以对积分结果用withsol进行参数索引withsol(sol1, intop_aperture(ebfd.nPoav), setval(zpos, index))这里zpos是焦点位置参数名index是当前扫描点序号。实际写的时候语法要看COMSOL版本中的函数签名来微调但思路就是通过withsol实现对扫描解集的逐点提取。如果只是看某个具体扫描位置的场分布也可以直接在“二维绘图组”里把数据集切成指定参数值这个操作会直观很多。4.3 从峰谷差推算出材料参数仿真跑通后再回头看数据思路就清晰了。闭孔Z-scan的归一化透过率峰谷差值ΔT_pv与非线性相移幅度存在近似关系ΔT_pv ≈ 0.406 * (1 - S)^0.25 * |ΔΦ0|这里的S是孔径线性透过率与总光强的比例因子ΔΦ0是焦点处样品引起的非线性相移。若假设从仿真曲线上读出ΔT_pv等于0.12S等于0.4那么可以反推|ΔΦ0| ≈ 0.12 / (0.406 * (1 - 0.4)^0.25) ≈ 0.356 rad再看ΔΦ0的定义ΔΦ0 k * n2 * I0 * L_effk是波数I0是焦点处峰值光强L_eff是样品有效厚度。如果k、I0和L_eff都已知n2就能直接算出来。这个反推过程相当于给仿真结果一个“合理性校验”最终算出的n2是否和你在材料卡里设定的n2接近如果差太多说明模型还有问题不是曲线光滑就算对。5. 实跑中踩过的几个坑对照形态和排查路径5.1 曲线整体抬高像一直在往上飘这种情况往往出现在后处理曲线中远离焦点的位置没有回到1而是整体高于1。一开始我以为是材料非线性设置错了查了一圈发现是PML没吸收干净高角度传播分量在边界弹了回来和主光束发生干涉把远场功率“顶”上去了。排查方法很简单把求解域纵向长度和PML厚度同时增加50%重跑两三个点看曲线是否明显回落。如果回落明显基本锁定是边界反射问题。另外如果用的是非归一化功率还要检查线性参考样品是否定义为α00且折射率与背景完全一致不要把线性本底吸收算成非线性抬升。5.2 非线性一点反应都没有这是另一个高频问题。曲线扫出来平平整整好像材料是纯线性的。我的排查顺序分成三步第一步看材料卡里的表达式确认折射率或吸收系数确实引用了光强变量而不是只建了一个变量却没参与方程第二步检查光强变量里的电场变量名是否和当前物理接口一致波束包络接口的电场变量名往往不是最常见的那个第三步看求解器设置如果启用了“自动忽略小项”或材料的非线性依赖没有参与雅可比矩阵很多非线性贡献会在迭代中被按零处理。说实话依赖局部电场强度的材料属性在内置材料卡里形成雅可比并不总是顺利。如果发现材料非线性始终不响应最靠谱的办法是把非线性项以源项形式写进自定义PDE而不是硬塞到材料里。5.3 参数扫描结束了绘图却是空的COMSOL提示绘图为空也困扰过我很久。原因大多是当前绘图组绑定的数据集是常规“解1”里面只存了初值或单点解参数扫描生成的完整数据都在另外的数据集里。打开数据集下拉框把数据集切换到带“参数化扫描”字样的那个数据集再选择正确的参数值索引图形就会出来。导出数据时同理导出节点也要重新选择数据集。如果不选导出的Excel或txt里只有第一组参数下的结果看着就是一条残缺曲线特别容易让人误判为求解失败。5.4 曲线峰谷不对称或者方向反了峰谷不对称通常有三个来源扫描步长非均匀导致取点密度不够、网格在焦点附近不等价于两侧、或边界反射依然没有排干净。先加密焦点附近网格和扫描步长再检查PML基本能解决大部分不对称问题。峰谷方向反了则大概率是入口相位符号问题。高斯光束在入射边界上的二次相位项符号决定了焦点位于样品前方还是后方相位符号取反等效扫描方向就反了最终曲线峰谷次序自然颠倒。这种情况不要动物理参数把边界表达式里的相位符号反过来再试即可。6. 写在最后一个很实用的调试小习惯如果你们刚接触这个方向我建议不要急着把实验室的量程全部塞进仿真。先在一个较简单的模型上跑出开孔曲线和闭孔曲线验证峰谷方向与ΔT_pv的近似解析关系能对上再慢慢加厚度、加吸收、加移动网格这些复杂条件。我自己保留了一个习惯每一轮扫描都同时跑一组n20的参考样本把它和非线性结果一起画在同一张图里。这个线性参考有没有异常比非线性曲线本身更能说明模型有没有出问题。也就是说Z参数扫描能不能成为一套可靠的分析方法关键不看你点了多少次参数化扫描按钮而看你是不是真的理解每一次参数变化在物理上代表什么。以上就是我在COMSOL非线性材料Z参数扫描建模和分析过程中踩过坑后沉淀下来的实际经验希望能让正走在这条路上的人少熬几个夜。