做非线性光学的仿真，二次谐波产生（SHG）是绕不开的第一道坎，也是我觉得最适合拿来上手 COMSOL 的模型之一。它物理图像清晰、方程形式简单，但真在软件里搭起来，又会冒出不少反直觉的坑。我自己从二维波导倍频做到块状晶体聚焦倍频，再往后加了周期极化结构，前后掉过不少头发。这篇就专门讲 COMSOL 仿真二次谐波产生这件事：从非线性极化的物理公式，到如何在 COMSOL 里把倍频源项写进方程，再到怎么把转换效率算得和理论对得上。
目标读者是两类人：一类是刚开始用 COMSOL 做非线性光学的同学，想找一个能跑通又符合物理的参考模型；另一类是已经跑通线性仿真、想往非线性方向扩展的工程师。文章不会只给步骤，还会解释每个设置背后的原因，读完你应该能照着搭出一个自己的倍频模型，而不是只会抄参数。
1. 二次谐波产生的物理基础与仿真映射
1.1 SHG的物理本质：二阶极化率如何成为倍频的源
二次谐波产生的本质，是介质在强光作用下出现了非线性极化。基频电场E(ω)通过二阶极化率χ(2)在介质内激起极化强度：
P(2ω) = ε0 * χ(2) * E(ω) * E(ω)
这个极化强度可以理解成一个频率为2ω的振荡偶极子源，它不断向外辐射二次谐波。这里有个容易搞混的点：基频光并不是“变成”了倍频光，而是非线性极化把能量耦合到了2ω的电磁场模式上。所以在 COMSOL 里建模，最常见也最不容易出错的做法，是把P(2ω)直接当作一个等效电流源加到倍频方程里，而不是去修改材料的介电常数或折射率。
往后看转换效率，主要看三件事：泵浦光强、有效非线性系数、以及作用长度内的相位失配。相位失配Δk = k(2ω) - 2k(ω)，当Δk * L远大于π时，倍频光一边产生一边又往回转换，最后输出的功率会周期振荡甚至趋近于零。真实器件里用双折射相位匹配或准相位匹配来补偿Δk，这也是后面仿真里必须手动控制的关键变量。
1.2 为什么用频域接口而不是瞬态求解
不少教程一上来就建议用“波动光学，电磁波，频域”接口，但很少说清楚为什么可以用频域。关键在于连续光SHG是一个单频稳态过程：泵浦是单频的，输出也是固定频率的2ω，整个系统里只有这两个频率分量，完全可以用频域求解。
瞬态仿真当然也能做倍频，但它要处理脉冲包络的时间演变、色散以及数值稳定性，模型规模和计算量都会大很多。对于入门和器件设计，先不用考虑瞬态。
频域建模的思路其实很简单：把ω和2ω这两个频率看成两个独立的物理场接口，分别求解Maxwell方程，再用一个耦合项把基频场的平方作为源引到倍频接口上。一旦建立了这个“双频域接口＋非线性源项”的框架，后面所有设置都顺理成章。
1.3 单向耦合还是双向耦合：仿真策略的取舍
严格来讲，SHG是基频和倍频互相耦合的过程：倍频增强之后会反过来消耗基频功率，这就是泵浦耗尽。但在大多数工程场景里，转换效率只有百分之几甚至更低，泵浦耗尽可以忽略。这时可以用单向耦合：先求基频场，再把基频场作为已知源去算倍频场。单向耦合的好处是数值稳定、调试简单，结果也足够准。
COMSOL里做单向耦合最简单的方式，是定义两个电磁波频域接口。一个在频率f1，负责基频；另一个在频率f2 = 2*f1，负责倍频。基频接口的源是端口边界条件指定的输入功率，倍频接口的源是手动添加的体电流密度。电流源表达式里直接引用基频接口的电场分量，比如写成ewfd1.Ez。
如果转换效率超过10%甚至更高，单向耦合会明显高估输出功率。这时才需要考虑双向耦合：在基频方程里也加入与倍频场相关的反向耦合项，两个接口一起迭代。这种做法的收敛性会差很多，我一般只在做高Q微环谐振增强时才会硬着头皮去调。
2. 动手建模前，先把这几个关键决定做对
2.1 几何选型：块状晶体、波导还是微环
很多人第一反应是“我要仿谐振腔”，但第一个模型最好别碰腔体。块状晶体倍频适合模拟入射光强低、作用距离短的实验；波导结构因为有模式束缚，有效模面积很小，单位长度转换效率比块状晶体高好几个量级；微环或谐振腔则靠循环功率增强，效率可以非常高，但设计参数多，调试麻烦。
我建议从二维波导开始。二维模型几何简单、网格量小，几分钟就能算完，还能清楚看到基频模场与倍频模场的重叠情况。等二维结果和理论公式对上了，再往三维扩展或者改微环都不迟。
实际器件里，波导截面的宽度、脊高、衬底折射率差都会影响模式分布，进而影响有效模面积和重叠积分。拿铌酸锂脊波导来说，宽度和脊高的变化会让模场从泄漏到强束缚之间切换，效率差距可以到几个数量级。所以几何选择本质上是在选“模式约束强度”。
2.2 材料参数：色散、非线性系数与单位换算
COMSOL内置材料库确实有不少光学材料，但它的折射率多数是常温下的静态数据，色散拟合精度对非线性仿真不够。更好的办法是查对应晶体的Sellmeier方程，自己把折射率随波长的变化写进模型。以铌酸锂为例，常温下ne在1064nm大约2.233，在532nm大约2.320，这0.087的折射率差会让相位失配严重到几乎无法产生倍频，所以必须靠温度调谐或周期极化来补偿。
非线性系数χ(2)需要自己换算。文献里更多给的是d系数，单位是pm/V，二者关系是：
χ(2) = 2d
比如铌酸锂z向极化的d33大约是25pm/V，那么χ(2)就是5×10^-11 m/V。这个单位换算我见过太多人栽跟头：数值差一个数量级，效率就差两个数量级，最后怎么都调不回来。
另外要注意张量方向。不同晶向允许的非线性系数通道不同，不是E^2随便乘一个标量就行。对于z切铌酸锂，如果基频和倍频电场主分量都是z方向，走的是d33通道；如果是x或y偏振，用的可能是d31或其他分量，效率会明显低一截。
2.3 散射边界和网格的取舍
模型边界如果不处理，电磁波会在边界反射，反射波和入射波混在一起，最后算出来的效率全乱。对于一进一出的波导，输入截面要用端口边界条件给定模式功率，输出端和侧壁要么用散射边界条件（SBC），要么用完美匹配层（PML）。SBC对垂直入射吸收比较好，但对斜入射会有残余反射；PML更稳但会增加网格数量。
网格尺寸的基准是每个波长方向至少5到6个网格。注意这里说的波长是介质中的波长，不是真空波长。对折射率n的材料，实际波长是λ/n，所以最大网格尺寸要取到λ/n / 5，边界处最好再额外加边界层网格去刻画倏逝场。
网格太粗会引入数值色散，等于把相位失配量给偷偷改变了。你明明材料参数都设对了，最后倍频效率还是不对，很可能就是网格搞的鬼。
3. 实操全过程：1064nm到532nm的波导倍频模型
3.1 参数设定与几何构建
下面给一个可以直接照抄的参数起点。模型是x方向传播的脊波导，z切铌酸锂，宽度1.2μm，脊高0.6μm，衬底厚度取2μm，包层是空气。芯片长度先从1mm开始，验证没问题再拉长。
| 参数 | 数值 | 说明 |
|---|---|---|
| 泵浦波长 | 1064 nm | 基频波长 |
| 倍频波长 | 532 nm | 二次谐波波长 |
| 波导宽度 | 1.2 μm | 脊宽 |
| 脊高 | 0.6 μm | 脊刻蚀深度 |
| 衬底厚度 | 2 μm | 下层铌酸锂 |
| 泵浦功率 | 100 mW | 端口输入功率 |
| 波导长度 | 1 mm 起 | 可后续拉长 |
| 极化周期 | 约6.1 μm | 准相位匹配周期 |
几何构建没什么特殊技巧，画一个长方形覆盖脊波导截面，再画一个衬底矩形，材料分别赋给铌酸锂和空气即可。注意COMSOL几何里如果把整个晶体都当作铌酸锂来画，基频模式可以跨整个基底传播，这和真实刻蚀波导不太一致。实际模型要在脊和衬底之间区分开，才能得到接近真实的模式约束。
材料部分，铌酸锂折射率用Sellmeier公式拟合，ne和no分开给。非线性系数不需要在材料节点里设置，直接写在后续的体电流源表达式中即可。
关于周期极化，我的做法不是真的在几何里画一排正负畴区，那样网格会爆炸。更高效的做法是在体电流源里加一个符号调制项：
χ_eff(x) = χ(2) * sign(cos(2π x / Λ))
这样就能在均匀材料里近似模拟周期极化，而不用网格去解析每一个畴壁。
3.2 物理场设置：端口、体电流源和变量映射
两个电磁波频域接口的设置分别是：
基频接口（假设叫ewfd1）：
- 频率设为
f1 = c_const / lambda1。 - 输入截面用“端口”边界，模式设为“TE”或“TM”，取决于你要激励的模式。
- 端口输入功率设为
100mW。 - 求解后要能输出基频电场分量，方便倍频接口引用。
倍频接口（假设叫ewfd2）：
- 频率设为
f2 = 2*f1。 - 无需端口输入，但要添加体电流密度源。
- 电流源表达式沿z方向，写成：
J_z = i * 4 * pi * f2 * eps0_const * chi_eff * ewfd1.Ez^2
这里有两个关键点，稍微不注意就会全军覆没。
第一，ewfd1.Ez是频域复振幅，直接平方后，时域上对应的正好是2ω分量。如果手滑写成abs(ewfd1.Ez)^2，那就等于只取了模的平方，物理上对应的是直流分量而不是倍频分量，输出会完全错误。这个坑我踩过不止一次，每次换模型重新写表达式时都要提醒自己。
第二，电流源里带i是频域方程的习惯写法。因为时域里P(2ω)对时间求两次导，对应频域中要乘上i 2ω再乘i 2ω？严格说，倍频极化源在无泵浦耗尽近似下，等效电流密度是J = ∂P_nl / ∂t，频域中就是J = i * 2ω * P_nl。所以J_z = i * 2 * ω2 * ε0 * χ_eff * E_z^2，而ω2 = 2π f2，2 * ω2 = 4π f2，最终表达式就是我上面写的那条。
在COMSOL变量里用ewfd1.Ez引用另一个接口的场，前提是两个接口求解顺序正确。我一般用两步研究：第一步只算ewfd1，第二步再算ewfd2，这样倍频求解时基频场已经是已知量，不会形成隐式耦合。如果强行把两个接口放在同一个频域研究里一起解，一般会得到完全不收敛或者物理上不正确的解。
另外，倍频模式也有可能泄漏到衬底或包层，网格必须能解析倍频电场。如果只按基频波长设计网格，倍频波长短一半，网格很可能不够。最大网格尺寸设置时直接按倍频介质波长来设，最保险。
3.3 研究顺序与后处理：把转换效率算出来
研究设置按顺序做两步：
- 第一步：只求解ewfd1，得到基频场。
- 第二步：求解ewfd2，利用上一步结果作为源，得到倍频场。
倍频输出功率在COMSOL后处理里，可以用输出截面上的坡印廷矢量积分来算。定义一个积分算子intop施加在输出边界上，然后在全局变量里写：
P_out = intop(ewfd2.nPoav)
这里ewfd2.nPoav是倍频接口的平均功率流法向分量。如果你用的版本变量名不完全一样，可以在结果里手动插入“表面积分”，表达式选ewfd2.nPoav即可。
转换效率定义为：
η = P_out / P_pump
在低转换效率极限下，理论上倍频功率正比于泵浦功率的平方、作用长度平方、有效非线性系数平方，反比于有效模面积。所以换算成归一化效率时，常用单位是%/W/cm²。仿真结果如果和这个标度关系对不上，比如扫泵浦功率时η不是线性增长，基本可以肯定模型里有问题。
我习惯的计算流程是这样的：先固定泵浦功率100mW，扫波导长度，看倍频输出是否随L²增长；然后再固定长度，扫泵浦功率，看η是否线性变化。两个扫描都符合预期，才认为这个模型是可信的。不要第一步就直接拿一个长波导去算最高转换效率，一旦结果不对，你很难判断是物理问题还是数值问题。
4. 常见问题与排查技巧实录
4.1 倍频效率几乎为零？先手动算相位失配
如果你照上面的步骤搭完模型，发现倍频输出低到小数点后一大堆零，第一反应别去怀疑COMSOL，先手动算一次相位失配。
用前面的铌酸锂例子，1064nm处ne1 = 2.233，532nm处ne2 = 2.320，波导长度L = 1mm。相位失配近似为：
Δk = 4π * (n2 - n1) / λ1
代入数值：
Δk ≈ 4π * 0.087 / 1.064e-6 ≈ 1.03e6 rad/m
再看Δk * L ≈ 1030，远远大于π。这说明倍频光在第一个振荡周期都没走完，输出几乎归零，纯粹是相位匹配没做好，不是软件问题。
解决办法是准相位匹配。周期取：
Λ = 2π / Δk = λ1 / (2 * (n2 - n1)) ≈ 6.1μm
在体电流源里加上sign(cos(2π x / Λ))调制后，等效每走一个周期就给相位失配补一次π，倍频输出会指数级上升。做这一步时，你会非常直观地体会到“相位匹配才是SHG的灵魂”。
4.2 不收敛、场图振荡、报错：按顺序排查
倍频仿真最常见的数值问题有三个：场图出现棋盘状振荡、求解器报奇异矩阵、或者收敛迭代慢到无法接受。
遇到这些问题，排查顺序应该是：
- 网格太粗：把最大单元尺寸减半，看倍频输出是否明显变化。如果变化很大，继续加密，直到结果基本稳定。
- 源项太尖锐：
sign(cos(2πx/Λ))这种周期符号函数在畴壁处不连续，局部电场梯度极大。可以用tanh(x/δ)近似平滑符号函数，δ取波长的1/20左右，既保留物理效果，又避免数值振荡。 - 求解器不合适：倍频接口里电场突变和强源项会让迭代求解器很难收敛，换用直接求解器（比如PARDISO）通常会好很多。
- 边界反射：如果输出端和侧壁用的是SBC，而波的入射角很斜，残留反射会影响倍频场。换成PML或者把边界拉远。
这套排查顺序基本覆盖了我遇到过的绝大多数情况。先把网格和求解器搞定，再去抠物理细节，效率会高很多。
4.3 结果和文献对不上：从这三个方向找原因
仿真结果和文献实验值对不上时，先别急着怪文献，也先别怪自己参数填错。我通常按下面几个方向排查：
第一，材料色散模型是否一致。不同Sellmeier公式拟合出来的折射率在可见光和中红外波段可能差出0.001到0.01，对相干长度的影响非常大。你用的色散和文献用的不一定是同一套，相位匹配位置就会偏移。
第二，模式重叠是否一致。二维截面近似和三维真实波导的模场分布不同，尤其近截止区。倍频效率正比于基频与倍频模式的重叠积分，如果重叠很小，效率就低。可以算一下归一化重叠因子，和文献对比。
第三，泵浦功率的定义是否一致。COMSOL端口功率是输入总功率，但非线性源用的是局部电场强度。如果波导模场面积很大，单位面积光强很低，即使总功率有100mW，倍频效率也可能低得可怜。反过来，如果用了归一化模场却没有正确缩放，输出的绝对功率就毫无意义。
下面这个速查表是我自己调试时一直用的：
| 症状 | 优先检查 |
|---|---|
| 倍频输出极低 | 相位失配、Λ是否设置正确 |
| 输出随长度振荡 | 网格不够、边界反射 |
| 输出随功率非线性增长 | 模式场缩放、归一化错误 |
| 求解不收敛 | 网格、求解器、符号函数太尖锐 |
5. 进阶方向：从静态模型走向真实器件
5.1 周期极化、温度调谐和微环增强
如果你跑通了均匀波导倍频，下一步强烈建议做周期极化扫描。把Λ设成扫描参数，从5μm扫到7μm，你会看到一个明显的倍频效率峰值。这个峰值位置和Δk计算值吻合时，说明你的色散和模式都设对了。
再进一步，可以把温度也加进变量。铌酸锂的折射率随温度变化明显，不同温度下最佳极化周期也不同。这样仿真就能指导真实实验的温度控制。
微环增强的倍频则更有挑战性。微环里循环功率比输入功率大很多，所以同样的非线性系数可以产生高得多的倍频效率。但要注意，微环的谐振频率随温度飘移严重，仿真里必须把热效应和相位匹配一起考虑，否则结果和实验很难对上。
5.2 模式重叠积分与波导截面的优化
随着模型跑通，你会越来越关注模式重叠。倍频产生的模场分布和基频模场不一定完全匹配。如果基频和倍频模场在截面上的空间重叠很小，即使相位匹配完美，效率也高不起来。
可以在后处理里定义一个重叠积分：
overlap = ∫ E_z(ω) * E_z(2ω) dS / sqrt(∫|E_z(ω)|² dS * ∫|E_z(2ω)|² dS)
然后用参数扫描分别改变波导宽度、脊高、刻蚀深度，观察重叠因子和有效模面积的变化。这个优化过程很考验参数设计能力，但也是做真实器件设计最有价值的部分。
5.3 脉冲与瞬态：什么时候才需要上瞬态
连续光倍频做熟了以后，你会碰到脉冲泵浦的问题。这时频域接口已经不够用，因为脉冲包含多个频率分量，而且时间包络会随位置演化。COMSOL里做瞬态倍频，需要用含时间项的波动方程，并考虑材料色散在时间域的等价形式。
瞬态模型的计算量远大于频域，网格和时间步长都要严格验证，通常只用来研究超短脉冲倍频的波形变换，或者腔体内光场的瞬态建立过程。如果你只是要稳态转换效率，踏踏实实用频域模型就好，别给自己找麻烦。
做完了这些工作，我越来越觉得COMSOL仿真二次谐波产生，真正有价值的部分不在软件操作，而在把物理过程拆成“源项、相位匹配、模式重叠”这三个可以分别量化的环节。每一个环节对应一个参数或一个后处理量，只要一个一个地验证，模型就能从“能跑”变成“可信”。
我个人调试模型时养成的最后一个习惯，是每次只改一个参数，扫完一个变量之后记录结果，再改下一个。这个习惯看着慢，实际上帮你省掉大量因为多因素纠缠而浪费的时间。如果你照着前面流程搭好了模型，先固定其他参数，只扫一遍极化周期，看到倍频输出从几乎为零跳到峰值的那一刻，你对SHG的理解就算真正落地了。
