COMSOL超表面仿真：从S参数相位到GH位移计算的完整流程

做超表面仿真的人迟早会遇到这样一个场景：透过率、反射率的谱线调得很顺眼，但对方突然一句“相位响应呢”，整个思路就开始冒汗。我最近在做一个光子晶体超表面的透反射分析时，被问到了一个更具体的量——光束斜入射后反射光在界面上的横向位移，也就是GH位移。这名字听起来有点唬人，其实就是反射光束相对几何光学预测位置发生的一个微小横移。普通光滑界面上，这个偏移小到几乎没人注意；但光子晶体和超表面里的共振会把相位随角度的变化放大，GH位移能被推到几十上百个波长，甚至变成负位移。

这篇把从COMSOL建模、S参数相位计算，到最终算出GH位移的完整路径拆开来讲。全文不会绕开坑，适合已经能跑基本电磁仿真、但还没认真处理过相位数据的同行。你要是正准备做光子晶体超表面的相位分析，这篇文章能帮你省下好几个晚上的无效调试。

1. 为什么超表面的GH位移绕不开相位计算

1.1 振幅曲线再漂亮，相位也是另一个世界

超表面仿真的常规流程是扫一条透过率或反射率曲线，看共振峰的位置、半宽、消光比。这套思路验证“这个结构有没有响应”是够用的，但一旦涉及到波前调控、相位梯度、光束位移这些量，只看振幅就完全不够了。GH位移是最典型的例子——它本质上就是由反射或透射系数相位对角度的依赖关系决定的，跟振幅大小没有直接关系。

理解这一点有个很直觉的模型。一束有限宽度的斜入射光束，可以分解成一系列传播方向略有差异的平面波，每个平面波在界面上经历一个不同的反射相位。如果反射相位随入射角变化非常剧烈，那么反射后各个角谱分量的相对相位关系会被明显改写，实空间叠加出来的波包峰值位置就必须横向平移，才能重新满足相干条件。这个过程跟士兵列队经过斜坡类似，每个人脚下延迟不一样，整列队伍自然就横着漂一段。所以，透反射相位的“斜率”，直接决定了GH位移的大小。

1.2 目标拆解：这个项目到底要算什么

给定一个周期结构，项目的最终目标是输出两个归一化数值：反射GH位移Δ_r/λ，和透射GH位移Δ_t/λ。拆开来看流程很清晰，一共四步。

第一步是建模型，在COMSOL里做出单晶胞结构，并用端口边界得到复透射系数和复反射系数，对应S参数里的S21和S11。第二步是扫参，把入射角从接近正入射一直扫到大角度，得到每一角度下的S11、S21复数值。第三步是取相位，将S11的相位φ_r(θ)和S21的相位φ_t(θ)取出来，并且必须做相位解包裹。第四步是对角度做数值微分，代入经典公式得到归一化位移。

这个流程对任何光子晶体、超表面结构都通用。哪怕你对GH位移本身不感兴趣，中间产出物——φ_r(θ)和φ_t(θ)两条曲线——本身就是很好的设计参考，后面做相位梯度器件时是现成的基础数据。

1.3 核心公式与一个容易翻车的单位问题

GH位移的经典计算公式并不复杂，反射和透射分别写成：

Δ_r / λ = -(1 / 2π) · dφ_r / dθ Δ_t / λ = -(1 / 2π) · dφ_t / dθ

其中λ是真空波长，φ_r、φ_t分别是反射和透射系数的相位，单位是弧度，θ是入射角，单位也要用弧度，Δ/λ是无量纲的归一化位移。这个形式看起来简单到让人怀疑，但它抓住了问题的本质：只要你把相位对角度的依赖关系搞准，位移量就出来了。

实际处理时最容易栽在单位上。如果直接拿度数的θ数据和相位数据做数值差分，公式里的2π就完全对不上。我见过有人算出来的GH位移达到几十米，最后查下来就是角度忘了转弧度。建议所有计算统一用弧度制，COMSOL参数扫描的时候用角度做参数也问题不大，但后处理公式里一定要写清楚角度到弧度的转换。

给一个直观量级：假设共振附近入射角变化1°时，解包后的反射相位掉了2π，那么dφ/dθ大约等于-360（弧度制下就是-2π除以0.01745），归一化反射位移就是360除以2π，约57个波长。这个数量级只有共振超表面能做得出来，普通介质界面连这个数值的零头都达不到。

2. COMSOL建模：单晶胞、周期端口与斜入射的核心设置

2.1 单晶胞建模：一个圆柱怎么代表整片超表面

以最常见的硅纳米圆柱方阵为例，结构参数可以这样设定：晶格周期700nm，圆柱半径200nm，圆柱高度400nm，背景是空气，研究波段放在近红外。需要说明，这只是一个示意几何，换成其他光子晶体或者超表面单元，整个流程同样成立。

为什么能用一个晶胞代表整个超表面？因为默认结构是无限周期阵列，只要单元之间没有缺陷，单晶胞加上周期边界条件得到的电磁响应，与无限阵列对外部照射的响应是严格等价的。这不是工程近似，而是周期结构的天然性质。COMSOL里的周期性条件就是为这个场景设计的。

建模时我的习惯是这样：建一个长方体作为空气介质区域，尺寸是周期乘周期再乘空气层高度，然后把圆柱放进长方体中心，做一个布尔联合操作，让结构成为单一连续几何域，再分别指定硅和空气两个材料域。这里有个经验点，尽量别用“组装体”方式保留内部交界面，直接布尔求和，否则电磁场在界面上的连续性处理容易出一些很难排查的异常。

模型的顶部和底部都要留出足够的空气层。建议从端口边界到结构表面至少留0.5λ，严格一点可以到0.75λ以上。因为端口边界假定了理想平面波，结构产生的倏逝波和衍射调制波需要在到达端口前衰减掉，留得太薄会让S参数里混进伪响应。

还要区分一下2D和3D的适用范围。如果结构在某个方向无限均匀，比如一维光栅线条，完全可以用2D模型，网格量和计算量都小很多。但如果是圆盘方阵这类点阵结构，斜入射时横向波矢分量会同时出现在两个方向，用3D模型更稳妥。我在下面默认用3D展开，因为覆盖面更广。

2.2 斜入射的核心机关：周期条件里的Bloch波矢

斜入射仿真的关键设置，就是周期条件中的Bloch波矢。很多新手把端口边界里的入射角写对了，周期性条件里却完全没动，结果S11相位不随角度变化，或者算出来的结果没有任何物理意义。Bloch波矢的本质是描述相邻晶胞之间的相位差，斜入射时这个相位差自然不再是零。

设入射角为θ，真空波矢大小为k0，则Bloch波矢x分量为：

kx = k0 · sinθ k0 = 2π / λ0

在COMSOL的操作里，先在全局参数中定义lambda0和theta。lambda0定义为1um，theta以度为单位定义。然后进入周期性条件设置，选择两个相对面，在波矢分量表达式里写：

k0*sin(theta*pi/180)

注意，theta如果用了度单位，这里一定要乘pi/180转成弧度。

同时，顶部面设置端口边界时，端口类型要选“周期端口”，并开启波激励，入射方向指向模型内部；底部面同样设为周期端口，但不开启波激励，作为透射输出。端口模式里的极化方向要根据研究的是TE还是TM来定，两套极化建议分别扫描。

端口与周期条件必须协同工作。端口负责计算散射参数，周期条件负责让求解器理解斜入射时单元边界两侧场的相位关系。如果端口写了斜入射、周期条件还是默认零波矢，两者互相矛盾，解出来的S参数就是四不像。一个快速的确认方法：斜入射时后处理电场图里应该能看到波阵面带有倾斜走向，如果还是平面对称的样子，多半是Bloch波矢没有真正生效。

2.3 网格与求解设置：相位精度需要单独照顾

相位仿真对网格的容忍度比振幅低很多。振幅结果偶尔粗糙一点还能看趋势，相位只要网格不够细，共振附近的导数就会出现锯齿，GH位移直接没法看。我的经验参考值是：空气区域最大网格尺寸取λ/5，硅圆柱内部取λ/8，圆柱曲面圆周方向至少保证30个网格单元以上。如果结构有尖角或者薄的颈缩区，还要在那些位置做局部加密。

求解器方面，频域计算一般首选直接求解器。迭代求解器在共振频率附近经常出现收敛困难，换直接求解器通常能解决，同时把求解容差收紧到1e-6。频率扫描时逐点求解，扫角度时同理。需要注意，参数化扫描如果同时包含频率和角度，计算量会膨胀得很厉害，建议固定频率扫角度，这样效率最高。

网格无关性验证对相位研究尤其重要。我习惯用粗、中、细三套网格分别跑正入射的S11相位谱，对比三条相位曲线，偏差控制在0.5°以内才算合格。如果共振峰附近偏差明显，以细网格结果为准，同时检查几何特征有没有被网格正确解析。这一步虽然花时间，但它能避免后面所有结果的相位误差被微分放大。

3. 从S参数到GH位移：相位提取与数值微分的完整做法

3.1 从COMSOL导出S参数相位：三个容易翻车的地方

COMSOL在周期端口设置正确的情况下，会自动计算S参数。在“派生值-全局计算”里，反射系数用类似ewfd.S11的变量名访问，透射系数对应ewfd.S21。要拿相位就用arg()函数，写成arg(ewfd.S11)和arg(ewfd.S21)。

第一处容易翻车的是只导出模值。有人习惯直接导出abs(S11)，等到后面想解包裹的时候发现相位信息全丢了，只能重跑。正确做法是导出实部和虚部，自已在外部脚本里重组复数，这样信息不丢失，也不依赖COMSOL的相位输出格式。

第二处是端口编号。模型稍微复杂一点，端口编号可能不是默认顺序，S21到底代表哪个方向的透射要仔细核对。建议先把所有S参数列出来看一眼，确定每个参数对应哪个端口传输方向，再开始批量处理。

第三处是波激励没开。如果两个端口都没有激励，S参数可能是空值或零。顶部端口作为入射端口必须明确开启波激励，底部端口如果要看透射，可以设为接收但也要在端口设置里配置到位。

导出数据时，我喜欢导出一个包含五列的表格：角度、S11实部、S11虚部、S21实部、S21虚部。这样后面用Python或者MATLAB都能方便处理。

3.2 相位解包裹：一步不做，导数全废

COMSOL默认输出的相位会被压缩在[-π,π]区间，所以相位扫过2π边界时，看起来就是从+π直接蹦到-π。如果拿着这种折叠相位去做数值微分，dφ/dθ会出现巨大的虚假尖峰，GH位移结果直接失真到离谱。

COMSOL里没有直接的解包裹函数，所以常规做法是导出数据后在外部处理。Python里一行就能解决问题：

import numpy as np phi_r = np.angle(S11_complex) phi_r_unwrap = np.unwrap(phi_r) phi_t = np.angle(S21_complex) phi_t_unwrap = np.unwrap(phi_t)

np.unwrap的原理是，相邻数据点之间的相位跳变超过π时自动补上±2π。这也意味着，如果角度扫描步长太粗，共振区域相位变化极快，中间可能只有一个数据点，相邻跳变超过π甚至接近π，unwrap就会失效。所以解包裹之前得先确认角度的采样密度足以分辨相位变化。共振区间的步长加密到0.1°甚至0.01°是常有的事。

解包之后可以做一个直观的连续性检查：曲线整体应该平滑，除了真正的物理突变点之外没有突然的折角。如果解包后的曲线还有锯齿，要么是扫描步长不够，要么是底层S参数数据噪声太大，这两种情况的处理方向完全不同。

3.3 数值微分与结果解读：共振附近的GH位移能有多大

相位解包完成后，剩下就是数值微分。均匀角度网格下，中心差分是最直接的选择：

theta_rad = np.deg2rad(theta_deg) dphi_dtheta = np.gradient(phi_r_unwrap, theta_rad) GH_refl_norm = -dphi_dtheta / (2 * np.pi)

np.gradient在非均匀网格下也能用，会自动处理网格间距变化，这点比较省心。

相位导数对噪声极其敏感。如果dphi_dtheta曲线毛刺很多，不要直接去平滑GH位移结果，而是回到φ(θ)本身做平滑。比如用Savitzky-Golay滤波器，窗口长度不要超过总数据点的十分之一，阶数用1阶或2阶就够。窗口太长会把共振峰压平，导数反而失真。

给一个我实际见过的量级参考：某个介质超表面在固定频率下，反射相位在25°入射角附近产生导模共振，不到0.5°的角度范围内相位掉了大约2π，数值微分给出的归一化反射GH位移接近100。而同一个结构在远离共振的频率下，GH位移只有不到2个波长。这一对比很直接地说明了问题——超表面里做GH位移研究，本质上就是在找相位曲线最陡的那一段。

反射和透射的GH位移要分开算、分开看。很多情况下反射位移大得惊人，透射位移却很平庸，或者反过来。它们对应的物理过程不同，不要为了方便只算一种就下结论。

4. 常见问题与排查：六种现场实录

4.1 一个高效的验证路径：先用单层平板校准

我强烈建议，任何新模型跑GH位移之前，先做一个最简单的验证结构：一块无限大介质平板。比如500nm厚的硅膜，上下是空气，用同样的端口、周期条件、相位提取、数值微分流程跑一遍。

为什么要先跑平板？因为平板结构可以用传输矩阵法手算出精确的解析参考，不需要依赖COMSOL结果本身来验证自己。把COMSOL得到的φ_r(θ)、φ_t(θ)和解析值放在同一张图里对比，相位偏差小于1°说明建模、端口、提取、微分的整个链路是通的。链路通了，再换到真实的超表面结构，出现的异常才值得怀疑是物理特性，而不是数值流程搞砸了。

这一步多花半小时，但能省掉后面几天的无效调试。我接手一个仿真数据时，流程第一件事就是拿平板校准环境，确认整个链路没问题再跑正式结构。

4.2 极化、空气层、瑞利反常：三个容易被漏掉的细节

极化方向对GH位移影响极大。同一个结构，TE和TM激发的共振类型完全不同，相位随角度的响应也可能天差地别。在COMSOL里，端口极化角设置决定了入射波的极化，所以做GH位移仿真时，TE和TM最好分别建模、分别扫描、分别分析。如果偷懒只算一种极化，结果很可能以偏概全。

空气层厚度在前面提到了，这里再强调一次。端口边界离结构太近，高次调制波在端口处没有衰减干净，会在S参数里表现为额外的伪共振。我自己的经验标准是至少0.5λ，如果结构支持较强的衍射耦合，留到1λ以上更稳妥。

瑞利反常角是另一个经常被误判的点。对周期结构来说，在某些入射角度下，某个衍射级次恰好从传播模式切换到倏逝模式，S参数的相位会在这个角度附近出现一个非共振的转折。它有几个特征：出现在特定的周期和波长关系下，解包后的相位曲线呈现平台式转折，和极化关系不大，而且并没有伴随明显的振幅共振峰。如果不认识它，很容易把这个转折当成超表面共振导致的大斜率，算出一个不真实的GH位移。

4.3 材料色散与参数化扫描的建议顺序

如果研究波段较宽，材料的折射率色散不可忽略。我的建议是先用常数折射率把整个流程和代码跑通，确认无误后再切换到色散模型。因为色散会让S参数相位随频率的变化更复杂，一旦出现问题很难判断是流程bug还是色散引起的新物理。常数折射率阶段可以把所有调试变量降到最少，等流程稳定了，再一步一步加入复杂度。

参数化扫描时，频率和角度不要一开始就做成大而全的全网格扫描。推荐先固定一个远离共振的频率，粗扫角度，看相位曲线的大致形态；找到可能的共振区后，再做精细角度扫描。频率维度单独扫一次共振谱，确认好感兴趣的频率点，再回到那个频率做高密度角度扫描。这样能把计算量压到最低，同时保证共振区的采样密度足够。

5. 让GH位移落地：从仿真数值到器件方向

5.1 GH位移的传感灵敏度思路

GH位移对入射角和环境折射率非常敏感，这是它作为传感探针的核心优势。在折射率传感器设计里，常规思路是看共振峰的移动，但利用GH位移时，可以在相位曲线最陡的那个角度选定工作点，环境折射率稍微变化一点，工作点附近的相位斜率变化就会被放大成很大的位移跳变。

在COMSOL里做这个延伸并不复杂，无非是在参数扫描里增加一个环境折射率维度，然后提取不同折射率下的GH位移曲线，计算灵敏度S=d(Δ_r/λ)/dn。这套数据可以直接用于评估结构的传感潜力，不需要额外的物理场接口。

5.2 负位移与相位奇点的潜在用途

当相位随角度上升时，dφ/dθ为正，归一化GH位移变成负值，反射光束会向几何预测相反的方向偏移。负位移在普通界面上几乎看不到，但在超表面里可以做得很大，而且可以通过结构设计来控制方向和大小。

我在仿真中确实见过很强的负位移，当时第一反应是怀疑解包错误。后来回到原始相位数据，用局部斜率确认了它的真实存在。这类现象通常和相位边上的奇点临近有关，对光束路由、精密定位这类应用，等于多了一个可设计的自由度。当然前提是先把负位移的来源分析清楚，别把数值假象当物理。

5.3 值得继续深挖的几个方向

GH位移这个方向如果往下走，有几个路径比较顺。一是结合耦合模理论，用解析模型先预判共振位置和相位斜率，再拿COMSOL当数值基准验证，这种双轨方式在写报告时很有说服力。二是做多层或级联超表面，研究透射GH位移和多层之间的相位协同，目前公开资料里的系统结论还不多。三是关注圆偏振光入射时的自旋依赖位移，也就是光学自旋霍尔效应，它的计算思路和GH位移相近，但后处理时还需要提取偏振相关的场信息。四是把相位梯度超表面和GH位移结合起来，设计出人工可调的光束位移器件，这类结构在COMSOL里做数值验证并不复杂，但物理机制上比单纯超表面共振更有意思。

最后说点实操层面的体会。GH位移的公式不难，难的是把φ(θ)曲线算得干净。相位这个量对网格、端口设置、扫描步长都极其敏感，任何一环糙一点，微分出来的位移就会完全失真。所以给新接手这个方向的人，我的建议永远是同一句话：先拿单层平板把整条链路跑通，画出一条平滑的dφ/dθ，再上超表面结构。这套流程真正落地之后会发现，所谓GH位移，不过是一条相位曲线的斜率而已。