介质超表面的非线性谐波仿真在Comsol里属于典型的“看起来不难、做起来全是坑”的活。尤其是要把三次谐波、倍频以及功率依赖都放进一个模型里同时还要求转换效率可复现科研上常见做工程化模型的人也越来越多。这篇文章我会把自己实际搭建这套模型的过程完整拆开讲包括物理图像怎么定、非线性极化项怎么加、为什么功率依赖不能简单地写成“折射率变化”以及最终怎么从后处理里把转换效率稳定地算出来。适合正在做介质超表面非线性研究、或者想用Comsol复现相关论文结果的朋友参考。1. 仿真前先理清楚的物理图像与模型边界1.1 谐波信号来自“超表面结构”还是“材料体效应”大多数人在开始建模时第一个问题就是介质超表面的非线性到底该算在哪个位置。以我自己的项目为例选用的是高折射率介质纳米柱阵列材料是类似GaAs或者TiO₂这类三阶非线性较强的介质。这里的物理机制要区分清楚谐波信号的主体来自纳米柱内部整个体积里的非线性极化而不是单纯来自结构的表面。这是关键差别。金属超表面往往非线性主要来自表面等离激元热点区域建模时通常要单独建一层表面非线性电流而介质超表面则更像是“体传输效应”——基频光在纳米柱内部形成局域场增强局部电场强度E(ω)远高于入射场然后在3ω频率处产生三阶极化P⁽³⁾ ε₀χ⁽³⁾E(ω)³。你在Comsol里设置材料参数时如果不把这个增强过程算进去最后得到的转换效率至少要差一两个数量级。建模之前还需要判断一下你研究的物理过程是真正的远场三次谐波还是混频产生的其他频点。标题里提到“倍频模型以及转换效率计算”倍频对应二阶非线性χ⁽²⁾三倍频对应三阶χ⁽³⁾。一个模型里同时实现两者时要格外小心频率之间的耦合顺序否则后处理时很难区分2ω信号里哪部分是真正倍频产生的、哪部分是其他非线性过程的串扰产物。我建模型时的做法是先建立一个“无超表面结构的平整衬底”模型把基频入射场跑一遍存下基准结果。后续所有非线性结果都减去这个基准就能初步分离体效应和结构诱导的增强效应。1.2 单晶胞仿真与周期边界成立的隐含条件Comsol里做超表面最常见的方式是仿真单个结构单元四面加周期性边界条件。这个做法成立的前提是结构周期远小于工作波长且阵列中没有长程相位梯度。一旦你的超表面引入了相位梯度比如用于波束偏折单个晶胞用周期边界就不再成立因为出射波里已经出现了多个衍射级每个衍射级对应一个独立的端口能量。我实际遇到的情况是介质纳米柱周期约800nm工作波长在1.55μm左右这种情况下单个晶胞加Floquet周期条件是成立的因为周期比波长小高阶衍射级基本截止。如果你的设计中有相位梯度就要改用广义的Floquet边界或多端口模型否则转换效率会算出一个“假值”——主要原因是你把本该衍射到其他级次的能量错误地导走了。1.3 泵浦功率在仿真里的角色要提前想清楚标题里“包含功率依赖”这几个字实际上有两种解读第一种基频功率不同非线性极化项随之变化谐波输出与基频功率不是线性关系转换效率随功率升高而上升随后由于损耗、泵浦耗尽等原因导致饱和甚至下降。第二种强泵浦下材料折射率本身发生改变也就是克尔效应n n₀ n₂I反过来影响基频场分布。我会在后面的第3节详细介绍在Comsol里如何实现这两类功率依赖。但一开始建模时你必须先确定自己要的是哪一种“功率依赖”因为它们的建模路径完全不同前者主要是在非线性源项里做自洽迭代后者则需要修改材料的折射率表达式。2. Comsol里的材料铺路线性参数、非线性张量与坐标取向2.1 线性介电特性怎么给才不踩坑非线性仿真最忌讳的是线性背景都没设对就开始堆非线性项。介质超表面材料在基频和谐波频率处的折射率必须分别指定不能只给一个“材料折射率”。尤其当你可以考虑多物理场耦合或色散介质时建议在Comsol材料节点里建立色散模型使用插值函数定义不同波长下的折射率。如果不考虑色散把1.55μm和517nm用得同样的折射率那么相位匹配条件天然就错了转换效率的绝对值也会失真。损耗项也不可忽视。介质材料在近红外通常损耗接近零但在三倍频的可见光波段可能吸收明显。在非线性光学中损耗不仅会吸收谐波还会影响基频在结构内的能量分布。建议先查材料的消光系数κ再估计吸收是否会显著影响结果。Comsol中定义复折射率n iκ或复介电常数都可以关键是Dim域别填错否则后处理时无法区分功率损耗来源。2.2 χ⁽²⁾和χ⁽³⁾的写法张量分量与坐标取向Comsol本身没有专门的“非线性极化率”输入界面所有非线性项都要通过修改方程或添加源项来实现。在这之前材料数据必须整理成张量形式。对于三阶非线性各向同性介质中独立的χ⁽³⁾分量并不多通常写成P⁽³⁾ᵢ ε₀ ∑ⱼ,ₖ,ₗ χ⁽³⁾ᵢⱼₖₗ Eⱼ(ω) Eₖ(ω) Eₗ(ω)由于频率参量很多要在三维情况下推导所有分量比较繁琐。我建议先在草稿纸上写清楚自己关心的是哪个偏振组合再用Comsol的变量表达式直接耦合场分量。最常用的一个简化是假设χ⁽³⁾只有对角分量且三个主轴分量相等于是Pₓ⁽³⁾ ε₀χ⁽³⁾Eₓ|E|²这样的近似对很多各向同性介质已经足够。而对于倍频二阶非线性极化写为P⁽²⁾ᵢ ε₀ ∑ⱼ,ₖ χ⁽²⁾ᵢⱼₖ Eⱼ(ω) Eₖ(ω)中心对称材料中χ⁽²⁾本征为零。如果你用的介质是GaAs闪锌矿结构或者LiNbO₃铌酸锂就必须给出完整的χ⁽²⁾张量否则倍频信号就是零。GaAs的χ⁽²⁾张量只有非对角分量且方向依赖于晶体主轴与超表面结构的相对旋转。Comsol中全局坐标系默认与几何坐标轴对齐如果你的纳米柱轴向旋转了45°那χ⁽²⁾张量也要跟着旋转。这里非常容易出现“一刷新就全错”的情况。2.3 引入功率依赖的两种做法有效折射率改写 vs 非线性源项自洽先说克尔效应式功率依赖。这种模型下材料折射率变为n n₀ n₂I其中I是局部光强。在Comsol里这个表达式不能直接写到材料折射率窗口中因为光强I本身是未知场E的函数会造成“材料属性依赖待求量”的隐性关系。我自己尝试时更推荐在变量节点中定义n_eff n0 n2 * (0.5cepsilon0n0abs(E)^2)然后在材料属性中使用这个n_eff。求解器需要迭代通常用辅助扫描或自洽循环完成。而对于三阶非线性源项功率依赖其实“天然存在”——因为P⁽³⁾本身正比于E³场强一变源项立即变化。你只需要保证求解过程中源项使用的是本步场解而不是上一步场解。在Comsol中要小心变量耦合的赋值逻辑如果你在边界设置里直接引用Eₓ、Eᵧ的默认解变量那么当同时求解基频和谐波时源项是自洽的但如果分步求解先算基频固定基频场再算谐波就变成了“固定泵浦”的非自洽近似。两种方法我都在项目中使用过前者叫全耦合模型后者叫分步注入模型。全耦合模型物理上更严格计算量也大得多分步注入模型更稳定适合参数扫描与效率快速估计。如果你要做“功率依赖曲线”个人建议分步注入模型就够了因为不同泵浦功率下先算基频场再注入到谐波方程中物理逻辑清晰数值上也避免了收敛困难。3. 倍频与三次谐波的耦合方程实现全耦合 vs 分步注入3.1 为什么Comsol默认接口不能直接表达非线性谐波Comsol波动光学模块中电磁波频域接口求解的是以下形式的方程∇ × (∇ × E) - k₀²εᵣE 0它默认介质是线性的。要引入非线性极化需要额外添加源项使方程变为∇ × (∇ × E) - k₀²εᵣE ω²μ₀P_NL这里的P_NL就是前面提到的二阶、三阶极化强度。这个方程形式上很简单但当你同时考虑基频、倍频、三倍频时的实际难点是P_NL的计算需要涉及多个频率的场乘积而默认的“电磁波频域”接口中只有一个频率自变量。所以我实际采用的方案有两种方案一全耦合多变量模型。不使用默认的单一频率接口而是自定义三个电磁波物理场接口分别对应频率f、2f、3f然后在每个接口中手动添加“电流源”或“外部源项”J_NL ∂P_NL/∂t。三个接口之间通过场变量互相引用。以三阶为例f频接口的源项中不包含3f场因为和频项需要两个正频一个负频而3f接口的源项中则包含(f, f, f)组合。这样三个方程在全局求解器中实现全耦合。方案二分步注入。先单独求解基频f的线性问题得到E(ω)。然后把E(ω)作为已知激励在2f和3f两个频率接口中分别加入源项P⁽²⁾(2ω) ε₀χ⁽²⁾E(ω)E(ω) P⁽³⁾(3ω) ε₀χ⁽³⁾E(ω)E(ω)E(ω)此时2f和3f方程之间没有相互逆向耦合即谐波不会反过来影响基频。这样的物理含义是“泵浦未耗尽”在小信号效率范围内近似很好。项目标题里要“功率依赖”和“转换效率计算”这种分步模型可以用最少的计算资源把功率依赖曲线扫出来。3.2 分步注入法的具体实现步骤以我的某次项目为例子详细说明分步法的Comsol操作流程设定参数组基频f₀ 193.5 THz 对应1.55μm波长倍频2f₀三倍频3f₀。材料折射率分别在三个频率点取值。定义非线性极化变量在“定义变量”中设置P2x、P2y、P2z和P3x、P3y、P3z。表达式引用基频解的电场分量这时候基频场是已知数。添加电流源对于2f接口在“域源”中输入-J_source电流源其中J ∂P/∂t的频域形式为J_NL jωP_NL。注意符号加入的是J_NL方程右端是吸收还是注入取决于你写的符号需要提前用平板模型确认方向。设置两个谐波接口的边界条件基频接口用Port输入 Floquet周期。2f和3f接口也使用Port边界但设置为“开放”或“无入射波”只输出。扫描泵浦功率反复修改基频端口功率P_in重新运行上述流程记录每个P_in下2f、3f的输出功率。这套流程跑下来效率计算就变成纯粹的“自动后处理”了逻辑清晰。3.3 全耦合模型什么时候有必要如果你的研究问题是强泵浦下的效率饱和、谐波对基频的反馈消耗那分步法就无效了。比如功率密度达到GW/cm²级别谐波转换效率达到10%以上泵浦耗尽效应不可忽略这时必须用全耦合模型。全耦合模型在Comsol中的实现我不会建议用“多物理场耦合”界面——它并没有现成的非线性光学耦合模块。更直接的思路是用“系数型偏微分方程”或“偏微分方程接口”改写三个频率的波动方程把非线性项写入系数或源项。这样做的好处是你可以完全控制方程中的每一项包括不同频率场之间的相乘关系。缺点是调试难度高耐人寻味的是一旦几何结构变化网格尺寸收敛性就有可能突然变差。个人经验在没有充足把握之前先做分步法确认物理趋势之后再逐步升级为全耦合。直接硬上一上来全耦合容易卡在“次数不足”、“不收敛”等问题上消耗几周时间还不知道哪里错。3.4 相位匹配介质超表面在这里最占便宜的地方转换效率公式里最关键的一项是和相位失配相关的因子η ∝ χ⁽³⁾² |E(ω)|⁴ L² · sinc²(ΔkL/2)其中Δk k(3ω) - 3k(ω)L是有效作用长度。对平整介质薄膜由于材料色散Δk通常不为零sinc²因子很快衰减所以效率极低。超表面的作用则是通过几何相位、导模共振等机制人为地提供额外的动量补偿让有效波矢匹配。所以在分析仿真结果时单纯看输出功率还不够最好把基频在结构内的相位分布也导出来计算局域波矢理解为什么能量转换效率这里高那里低。这一点在模型后处理中可以这样操作在基频结果中绘制E的相位剖面观察是否在一个周期内出现了明显的相位累积异常区域那些区域往往就是非线性源项贡献最大的地方。4. 转换效率的计算从后处理功率流到归一化4.1 怎样从Comsol中提取功率转换效率的定义简洁起见可采用η₂ P(2ω) / P(ω_in) η₃ P(3ω) / P(ω_in)但难在P(2ω)和P(3ω)的数值怎么取。很多人直接在“全局计算”里选择“电动率”变量却忽略了谐波频率下的功率流是复数坡印廷矢量的实部且应在指定平面上积分。正确的做法是在结构上方的空气域中画一条水平线二维模型或一个平面三维模型积分特定频率的坡印廷矢量的法向分量P(ω) ∫_S (1/2) Re(E(ω) × H*(ω)) · n dSComsol中对应变量通常类似基频emw.Poav单物理场接口默认自己定义频率时需要用real(0.5 * Ex * conj(Hy) - 0.5 * Ey * conj(Hx))这样的手写表达式关键的一点是要仔细选择积分平面保证平面离纳米柱顶面有一段距离避免近场局域增强带来的数值不稳定。我一般取距顶面0.5到1个波长处作为输出功率提取面。输入功率则取入射端口对应的总功率也可以直接在端口边界上积分P_in。4.2 数值一致性校验转换效率要经得起“拆解”在计算完效率之后我要求自己必须通过三道校验再下结论第一能量守恒粗校验对于无吸收材料入射功率约等于反射功率加透射功率加谐波功率。在Comsol后处理里把各个功率值都导出来如果总能量偏差超过5%多半是某个边界条件或积分面选错了。第二泵浦功率标度校验在小信号范围内三倍频输出的功率P(3ω)应正比于P(ω)³倍频输出正比于P(ω)²。在log-log坐标下画出扫描结果如果斜率不符合这个规律说明模型中有非线性项表达式错误或者功率依赖被放置的位置不对。第三对称性校验对圆形对称的纳米柱结构如果入射光是线偏振那么二次谐波只有特定方向的偏振分量如果仿真结果显示正交偏振分量异常大多半是张量坐标出了问题。4.3 一个可参考的典型数量级在普通介质薄膜上三倍频转换效率通常在10⁻⁶以下介质超表面通过共振增强场之后可以做到10⁻⁴左右。我们在自己的模型里纳米柱结构下计算出特定功率密度的THG效率在2×10⁻⁴量级这个数量级和文献趋势一致。如果你的仿真结果跳到了10⁻²甚至更高先别高兴大概率是单位或者归一化有问题——不是材料、网格、边界三种之一就是符号方向写反了。需要提醒的是很多论文中的效率是“归一化到结构单元面积”有的论文是“归一化到整个光斑面积”两者差着超表面占空比。阅读文献对比时务必先确认归一化方式否则你的数据跟文献对不上不是模型错了而是定义口径不同。5. 网格、边界与求解器实测调参经验5.1 网格尺寸必须向高频场妥协如果你设置的基频波长是1.55μm三次谐波波长只有约517nm。在计算网格尺寸时不能只满足“基频波长分辨率要求”而是要以三倍频场在结构内部的振荡为准。我通常要求最大网格边长不超过λ₃/6左右即约86nm在纳米柱内部和近场区域网格甚至更细。你可以这样检查网格质量先使用一个相对粗糙网格跑一遍分步模型然后把网格加密一倍再一次运行。如果转换效率变化超过20%说明网格还未收敛继续加密直到效率值趋于稳定。这一步是最枯燥但也最必要的。5.2 边界条件组合Floquet、PML与端口周期边界的设置在Comsol里需要与端口边界配合使用。对于二维模型两侧使用Floquet周期边界周期矢方向对应晶格矢量上下边界使用端口上方端口设置为“入射”下方端口设置为“仅透射”。对于三维模型四个侧面用Floquet上下同样用端口。PML不建议设在超表面结构附近因为非线性源项的数值稳定性容易受到PML层坐标系影响。我会在主要计算域的上方加一层空气间隔至少半个波长然后再叠加PML。单个PML的厚度设为2到3个边界波长即可太厚会引入非物理寄生模式。5.3 求解器选择与内存控制的实操体会分步注入模型的求解顺序建议是先求解基频接口此时2f、3f接口不使用或设置为非激活状态避免求解器尝试计算未初始化源项。冻结基频解激活2f接口单独求解倍频。再激活3f接口求解三次谐波。每一步都独立收敛、独立保存。如果一次同时求解三个频率Comsol默认的迭代求解器常常陷入“因非线性源项过大导致场发散”的问题。我的做法是使用分离式求解器并设置为“手动选择因变量”先解基频再解2f最后解3f每步更新非线性源项。这样虽然增加了点击次数但稳定性提升明显。内存方面三维介质超表面全耦合模型在百万级网格下动辄需要64GB以上的RAM。分步模型则可以明显降低峰值内存。如果内存不够优先压缩3f接口的网格密度而不是压缩基频的。6. 最容易出错的几个位置我的踩坑记录6.1 坐标系对应关系会让张量分量全部错误我踩过最深的坑之一把GaAs的χ⁽²⁾张量直接套进模型结果倍频信号消失了。检查后发现问题出在张量方向与晶体轴方向上。GaAs是闪锌矿结构它的二阶非线性张量是d₁₄ d₂₅ d₃₆ ≠ 0也就是说只有x(yz)、y(zx)、z(xy)这三个组合分量非零。如果你把纳米柱长轴沿全局x方向放置而晶体主轴沿全局y方向就会完全对不上。所有手册上的χ⁽²⁾数值都默认是“晶体主轴坐标”下的值你在Comsol中建模时必须做坐标旋转。这个旋转既可以用“几何变换”等命令实现也可以直接把张量分量重投影到全局坐标。6.2 “输入功率”取错了位置效率曲线整体平移另一个高频错误是输入功率定义用了端口处的总功率但端口处存在反射。如果端口反射率为R实际进入结构内并参与非线性转换的功率只有P_in(1-R)。在小信号效率计算中有人用P_in归一化有人用P_in(1-R)归一化两者相差一个系数。从物理角度说转换效率应该是谐波输出功率除以“实际进入结构内部”的功率这样才反映材料与结构的非线性转换能力。而从工程角度说很多时候人们更关心“外部入射多少光能产生多少谐波”这时用P_in归一化更实用。我建议在论文与报告中同时注明两种口径或至少在参数设置中把R单独提取出来方便后续调整。6.3 如何验证模型无超表面平板对照组是底牌非线性模型比线性模型更容易“算得出结果但结果错了”。我的习惯是每次建立超表面结构时同时建立“无结构平板”对照组保持基频功率、材料参数、边界条件完全一致。平板介质的谐波输出理论上可以直接通过薄介质二次谐波/三次谐波解析公式估算η_flat ∝ (χ⁽³⁾ E₀² L)²如果对照组仿真结果与这个近似公式偏差在一个数量级以上说明问题多半出在源项符号、坐标定义或材料单位换算上。这个对照组能帮你把排查范围缩小一大半。6.4 扫描泵浦功率时的数值陷阱最后关于功率依赖扫描一个特别容易忽视的细节是随机将端口功率设置提高几个数量级后场增强因子可能触发数值饱和。Comsol默认的相对容差1e-6可以应付大多数情况但当P_in从1mW扫到1W时基频场幅值变化上千倍非线性源项变化更大。建议在扫描参数时交替提高容差粗扫时用1e-5锁定趋势后加密扫描用1e-7。否则你会得到一条“看似合理”实则发散拼凑起来的效率曲线。