企业官网建设流程全解析

介质超表面的非线性谐波仿真，在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操作流程：

1. 设定参数组：基频f₀ = 193.5 THz （对应1.55μm波长），倍频2f₀，三倍频3f₀。材料折射率分别在三个频率点取值。

2. 定义非线性极化变量：在“定义>变量”中设置P2x、P2y、P2z和P3x、P3y、P3z。表达式引用基频解的电场分量（这时候基频场是已知数）。

3. 添加电流源：对于2f接口，在“域源”中输入-J_source（电流源），其中J = ∂P/∂t的频域形式为J_NL = jωP_NL。注意符号：加入的是J_NL，方程右端是吸收还是注入取决于你写的符号，需要提前用平板模型确认方向。

4. 设置两个谐波接口的边界条件：基频接口用Port输入+ Floquet周期。2f和3f接口也使用Port边界，但设置为“开放”或“无入射波”（只输出）。

5. 扫描泵浦功率：反复修改基频端口功率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 dS

Comsol中对应变量通常类似：

• 基频：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 求解器选择与内存控制的实操体会

分步注入模型的求解顺序建议是：

1. 先求解基频接口，此时2f、3f接口不使用（或设置为非激活状态），避免求解器尝试计算未初始化源项。
2. 冻结基频解，激活2f接口，单独求解倍频。
3. 再激活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。否则你会得到一条“看似合理”实则发散拼凑起来的效率曲线。