3. FEM数学基础与光子器件控制方程

有限元法求解光子器件的第一步,不是画网格,而是把麦克斯韦方程组整理成"可被弱形式处理"的偏微分方程。本节围绕三条主线展开:从时谐假设导出频域 Helmholtz 方程、从时域波动方程理解宽谱问题、以及通过分部积分把强形式转化为弱形式并解释各类边界条件的物理含义。这三条主线共同决定了后续网格尺寸、求解器类型与端口设置的合理性。

一、从麦克斯韦方程到频域 Helmholtz 方程

在无源、无自由电荷、线性各向同性介质中,麦克斯韦旋度方程为 ∇×E = −jωμH、∇×H = jωεE。对第一式取旋度并代入第二式,消去 H 后得到电场满足的矢量波动方程:

∇×(μ⁻¹ ∇×E) − ω²ε E = 0

若介质在传播方向上近似均匀(如直波导、光纤截面),可假设场沿 z 方向以 exp(−jβz) 传播,令 E(x,y,z)=E_t(x,y)exp(−jβz),则 ∂/∂z → −jβ。代入后得到横向截面的本征值问题:

∇_t×(μ⁻¹ ∇_t×E_t) − (ω²ε − β²/μ) E_t = 0

这就是光子器件中最常用的频域 Helmholtz 形式,β 为本征值,对应模式的有效折射率 n_eff = β/k₀。对 TE/TM 标量近似,可进一步简化为 ∇²φ + k₀²n²(x,y)φ = β²φ,即标量 Helmholtz 方程。理解这一点很关键:FEM 求解的不是"场长什么样",而是给定 ω 下哪些 β 能让方程有非零解,这正是模式求解与频域传输求解的本质区别。

二、时域形式与频域形式的适用边界

时域波动方程 ∇×(μ⁻¹∇×E) + ε ∂²E/∂t² = 0 不做单频假设,因此能描述脉冲展宽、非线性 Kerr 效应、时变调制等宽谱或瞬态过程。代价是需要在时间上逐步推进,且网格尺寸受最短波长与 Courant 条件约束。频域形式一次求解只针对一个 ω,但可高效扫描波长、直接得到模式与 S 参数。工程上常见做法是:先用频域模式求解确定波导截面与 n_eff,再用频域传输或时域 FDTD 处理器件级响应。选择依据可归纳为下表。

判据频域 Helmholtz时域波动方程
激励单频/扫频宽带脉冲
输出模式、S 参数、场分布时域波形、频谱
非线性/时变需迭代近似自然支持
典型工具FEM 频域、本征模FDTD、FEM 时域

三、弱形式推导:为什么必须分部积分

强形式要求场二阶可导,而光子器件中常出现材料界面、尖角、金属边界,二阶导数在界面处不连续。弱形式通过乘以测试函数 v 并在求解域 Ω 上积分,再用矢量恒等式 ∇·(A×B)=B·(∇×A)−A·(∇×B) 把二阶旋度降为一阶:

∫_Ω (μ⁻¹ ∇×E)·(∇×v) dΩ − ω²∫_Ω ε E·v dΩ
  − ∮_∂Ω v·(n̂ × μ⁻¹ ∇×E) dS = 0

边界积分项 ∮ v·(n̂×H) dS 正是边界条件的入口。它把"边界上场的法向行为"转化为弱形式中的自然边界项,因此 PEC、PMC、端口、PML 都可以通过修改这一项实现,而不必强行约束节点值。这也是 FEM 相比 FDTD 在处理复杂几何与开放边界时更灵活的原因。

从强形式到弱形式:边界条件的物理入口 强形式(微分形式) ∇×(μ⁻¹∇×E) − ω²εE = 0 要求 E 二阶可导 界面/尖角处导数不连续 弱形式(积分形式) ∫(μ⁻¹∇×E)·(∇×v) − ω²∫εE·v − ∮ v·(n̂×H) dS = 0 仅需一阶可导,界面自然满足 边界积分项 ∮ v·(n̂×H) dS 承载全部边界物理 PEC / PMC / 端口 / PML ×v, 积分 分部积分 边界条件在弱形式中的实现方式 PEC 理想电导体 n̂×E = 0,切向电场为零 PMC 理想磁导体 n̂×H = 0,切向磁场为零 端口 / 散射边界 n̂×H = 已知入射场 PML 完美匹配层 复坐标拉伸,吸收出射波 关键结论:弱形式把"边界上场的法向行为"转化为积分项,因此 PEC/PMC/端口/PML 都通过修改边界积分实现, 而不必强行约束节点值——这正是 FEM 处理复杂几何与开放边界的灵活性来源。

四、边界条件的物理含义与设置要点

PEC 对应理想金属,电场切向分量为零,常用于金属电极、金属包层或对称面。PMC 对应理想磁壁,磁场切向分量为零,常用于对称面或人工磁导体。端口边界则把已知入射模式作为激励,同时吸收反射波,是计算 S 参数的核心。PML 通过复坐标拉伸使出射波在层内指数衰减,模拟无限大开放区域,其厚度通常取 0.5~1 个中心波长,反射率可控制在 −40 dB 以下。

设置边界时需注意三点:第一,PEC 与 PMC 不能同时施加在同一边界上;第二,端口模式必须与波导截面本征模一致,否则会激励高阶模;第三,PML 应放在离器件足够远的位置,避免倏逝波在 PML 内产生伪反射。对于对称结构,可利用对称面把计算域缩小一半甚至四分之一,显著降低自由度。

五、一个可复现的示例:硅条形波导模式求解

以 220 nm 厚、500 nm 宽的硅条形波导为例,包层为 SiO₂,工作波长 1550 nm。材料参数取 n_Si = 3.48、n_SiO₂ = 1.44。求解步骤为:建立二维截面几何,设置 Si 与 SiO₂ 的折射率,外边界用 PML 或散射边界,求解本征值问题得到 β 与 n_eff。合理预期是 TE 基模 n_eff 约在 2.3~2.5 之间,具体数值随网格与边界设置略有变化。若 n_eff 接近包层折射率 1.44,说明模式已接近截止;若出现多个相近本征值,需检查是否激励了高阶模或伪模。

判断模式是否收敛的实用方法:把网格从 λ/10 加密到 λ/20,观察 n_eff 变化是否小于 1e-4。若变化仍大,说明网格或边界设置需要调整,而不是继续加密。

六、常见误区与排障思路

初学者常犯的错误包括:把频域 Helmholtz 当成时域方程使用,导致无法处理脉冲;在端口处直接施加 PEC 而不用模式激励,得到全反射结果;PML 厚度不足或离器件太近,导致 S 参数出现纹波。排障时可先检查边界条件类型与位置,再检查网格是否在材料界面处加密,最后核对材料折射率与波长单位是否一致。若结果出现非物理的场增强,优先怀疑 PML 反射或伪模。