最近不少同行来问我同一个问题COMSOL里模型跑完了近场的Ex、Ey、Ez都很好看但我想拿到某个观察方向上的远场偏振态比如到底线偏振还是椭圆偏振、方向角多少、圆偏振度多高怎么搞这一下就把人问住了。原因很实在COMSOL的远场后处理默认输出的是电场幅值和相位不会自动给你一组干净利落的偏振分量组合而偏振恰恰是远场中那种“两个正交分量的复振幅关系”才能定义出来的东西。我摸索了一套相对通用的远场偏振计算方法从理论图层层拆到COMSOL实操再到导出数据用脚本计算偏振参数前后花了不少时间踩坑。这篇把完整思路和可直接照抄的实现流程都写出来希望能帮有同样需求的人少走弯路。1. 远场偏振计算到底难在哪1.1 COMSOL的远场后处理默认给不了你想要的COMSOL电磁波频域接口里加入远场域Far-Field Domain之后计算完可以很方便地画远场辐射方向图查看远场模值emw.farfield.Efar或者归一化辐射方向。但要注意这个模值是标量它丢掉了两个正交分量的相位关系。而偏振描述的核心恰恰是相位差同样是两个等幅度的电场分量相位差为0时是线偏振相位差为$\pi/2$时就变成圆偏振。只看模值根本区分不了这两者。所以要算偏振第一件事就是改变思路不推荐直接在COMSUL图形界面上看远场图而是把远场电场分量作为数据提取出来用表达式或者外部脚本处理。远场计算本质上是一种外推近场求解完成后在远场域内部或边界上把等效电流、等效磁流通过格林函数积分外推到任意观察角度。这个过程在COMSOL里已经内置完成但它给出的是远处球面上每个观察点的复电场分量一般可访问为emw.farfield.Ex、emw.farfield.Ey、emw.farfield.Ez不同版本变量名有差异老版本常见emw.Efarx这种写法建议先在后处理表达式里验证一下再批量导出。这个复分量信息已经是计算偏振的全部原料了只是很少有人去把它展开成完整的数据流水线。1.2 偏振不是“横竖都能算”相位精度必须先过关有些朋友会想我直接看Ex的实部和Ey的实部比值不就是偏振方向吗但如果场是椭圆偏振相位差在中间状态只取实部显然会得到错误的“短轴”投影。偏振计算需要的是复振幅而且这个复振幅的相位必须足够稳定。这带来两个隐藏要求近场求解的网格精度必须吃得消相位误差。偏振对相位是否精确非常敏感尤其是结构尺度接近波长时网格稍微粗一点远场相位可以偏掉好几度偏振度马上就变了。远场外推所用积分边界的位置也会影响相位基准。COMSOL里远场域通常依赖PML完美匹配层或散射边界包裹外推边界离结构太近会截断倏逝波太远又会增加计算域。实操上我一般建议边界离源区至少半个波长同时确保该边界上的网格和背后PML过渡层连续。也就是说算偏振的成败并不是从后处理开始而是在建模阶段就已经埋下伏笔。这一点后面第5节还会细说。2. 理论层远场外推与偏振参数的数学地图2.1 近场到远场变换本质是“面天线等效源”的叠加COMSOL远场域所做的事情可以类比成天线理论里面经典的面等效源法把一个包围着散射体/辐射源的闭合面看成新的源面上每个位置都有一组等效电流和等效磁流这些等效源辐射出去在远方某个观察方向上对远场电场的贡献满足[ \mathbf{E}{far}(\theta,\phi) \propto \frac{e^{jk r}}{r} \int{\Sigma} \left[ j k (\mathbf{J} \times \hat{r}) \times \hat{r} j k \mathbf{M} \times \hat{r} \right] e^{j k \hat{r} \cdot \mathbf{r}} , dS ]公式看着吓人实际上意思很朴素把近场面上每个点的等效源做一个加权叠加权重由观察方向与源位置的相位差决定。远场就是近场信息经过一个“傅里叶变换式”的传播映射。落到COMSOL操作上你不需要自己实现这个积分只要做两步在电磁波频域接口里添加远场域节点并指定积分面通常直接选整个外部边界。求解后在结果里用远场电场分量表达式来取值。这里的“观察点”在COMSOL远场计算里是一组离散角度的采样。默认情况下它输出的是以原点为参考球的远场球坐标分量但我们做偏振分析时一定要把数据转成全局直角坐标系分量或者直接在直角坐标系下定义观察方向否则后续Jones矢量的构造很容易出错。2.2 偏振的数学工具其实就是三行公式有了远场复电场分量 $E_x$、$E_y$以及需要时 $E_z$偏振的全部信息就浓缩在以下几个量上。我习惯固定使用这套流程通用性很强。首先构造Jones矢量[ \mathbf{J} \begin{bmatrix} E_x e^{j\delta_x} \ E_y e^{j\delta_y} \end{bmatrix} ]由于COMSOL给出的分量本身已经是复数$E_x$、$E_y$的模和辐角都已经包含在复值里面直接拿复数操作即可。然后计算Stokes参数。Stokes参数的好处在于它是实数且与坐标系旋转有明确变换关系适合做后续可视化[ S_0 |E_x|^2 |E_y|^2 ] [ S_1 |E_x|^2 - |E_y|^2 ] [ S_2 2 \mathrm{Re}(E_x E_y^{}) ] [ S_3 2 \mathrm{Im}(E_x E_y^{}) ]接着可以定义两个直观的几何量。第一个是偏振椭圆的长轴方位角 $\psi$代表椭圆长轴相对坐标系x轴的偏转角度[ \psi \frac{1}{2} \mathrm{atan2}(S_2, S_1) ]$\psi$的取值范围一般在 $-90^{\circ}$ 到 $90^{\circ}$ 之间。第二个是椭圆率角 $\chi$表示椭圆“胖瘦”取正值时代表右旋圆偏振成分占优取负值代表左旋[ \chi \frac{1}{2} \mathrm{atan2}\left(S_3, \sqrt{S_1^2 S_2^2}\right) ]$\chi$的取值范围在 $-45^{\circ}$ 到 $45^{\circ}$ 之间。当 $\chi \pm 45^{\circ}$ 时就是纯圆偏振当 $\chi 0$ 时就是纯线偏振。用生活类比来说$\psi$ 是椭圆在纸上画出的那条主轴朝哪边倾斜$\chi$ 是这张椭圆变形得“扁不扁”、旋转方向如何。单一角度下如果场不是完全偏振的比如包含随机非相干背景还可以计算偏振度[ DOP \frac{\sqrt{S_1^2 S_2^2 S_3^2}}{S_0} ]对一个完全偏振的理想电磁波DOP等于1。实际仿真中如果计算域边界处理不好远场混入非物理数值噪声DOP会略小于1这个指标也可以反过来用来检查远场外推做得好不好。3. 通用计算方法从模型设计到数据后处理的完整流水线3.1 建模阶段就要为远场偏振留好口子我把这套方法定位成“通用”就是因为它和具体结构无关——不管是纳米天线、光栅、超表面还是射频天线阵列只要COMSOL频域求解能给出近场这套流程都能用。前提是建模阶段满足几个基本要求使用电磁波频域接口而且求解的是完整频点不是本征模。计算域外围加载PML或足够可靠的吸波边界。远场计算依赖外推面如果外推面外面还有较强的反射波远场方向图会被污染。在物理场设置中显式添加远场域并把外推面选在PML层之前的那一圈内部边界上。结果保存时勾选“所有场的解”或者足够密的角度采样否则导出远场数据时角度分辨率不足后续偏振扫描就没法看。通用流程我整理成了四步频域求解近场。基于近场结果做远场外推得到多个角度下的复电场分量。把复电场分量按直角坐标导出成结构化数据如文本文件或CSV。用Python或MATLAB脚本读入数据计算Stokes参数与偏振椭圆参数。这个流水线的好处是解耦。COMSOL内做的事情只管出“准”的场数据偏振计算全部放在脚本端换模型时不需要反复修改脚本逻辑只需要改数据列映射关系。3.2 远场域的关键配置项添加远场域节点后有几个参数要特别留意积分面积分的边界通常选择包围结构的最外层计算域边界也就是PML内部的边界。直接选整个计算域外表面是没问题的。参考波矢方向即默认的平面波入射方向它对远场的相位基准有影响。如果不是垂直入射要注意观察角度的定义。输出球坐标还是直角坐标COMSOL后处理里远场结果既有球坐标分量如Efarphi、Efartheta也有直角坐标分量。偏振分析我建议固定用直角坐标分量也就是用表达式把emw.farfield.Ex、emw.farfield.Ey取出来。因为球坐标分量在不同谱面的定义跟观察方向绑定后续写通用脚本时非常容易搞混。角度采样通常后处理远场就可以设置theta和phi的扫描点数。默认可能只有几十个点做偏振角度谱时会显得太粗建议至少设置到200到400个点。数据导出这一步可以直接在COMSOL“结果”里用一维绘图组以theta或者phi作为扫描变量把远场电场分量画成线图然后将绘图数据导出为文本。也可以用“全局计算”节点用一个表格列出多个表达式的值。两种方式各有利弊我后文会聊更推荐的方案。3.3 指标体系一个脚本吃透全部偏振参数写脚本时最重要的是把Stokes参数计算封装成一个可复用的函数。以下给出一个清醒的函数骨架Python风格理论上可用于任何远端数据import numpy as np def polarization_from_farfield(Ex, Ey): 输入复数远场电场分量 Ex, Ey直角坐标系 返回 S0~S3 与偏振椭圆参数 psi, chi, DOP Ex np.asarray(Ex, dtypecomplex) Ey np.asarray(Ey, dtypecomplex) S0 np.abs(Ex)**2 np.abs(Ey)**2 S1 np.abs(Ex)**2 - np.abs(Ey)**2 S2 2 * np.real(Ex * np.conj(Ey)) S3 2 * np.imag(Ex * np.conj(Ey)) psi 0.5 * np.arctan2(S2, S1) # 椭圆长轴方位角 chi 0.5 * np.arctan2(S3, np.sqrt(S1**2 S2**2)) # 椭圆率角 DOP np.sqrt(S1**2 S2**2 S3**2) / np.maximum(S0, 1e-12) return S0, S1, S2, S3, psi, chi, DOP这段代码读取的是复值数组如果你的导出数据是分实部、虚部列出进函数之前先拼成复数Ex Ex_real 1j * Ex_imag Ey Ey_real 1j * Ey_imag这样整个计算流程就只剩下“把COMSOL的数据喂给这个函数”这一步了。我实际使用时还会顺手把角度数据一并输出这样就能画出“偏振椭圆参数随方位角变化”或者“随波长变化”的曲线。4. 实操案例亚波长光栅零级透射远场偏振分析4.1 模型搭建与关键参数用一个非常常见的场景来演示可见光波段二氧化硅基底上做周期性亚波长介质光栅正入射平面波计算零级透射远场的偏振态。光栅周期设为600 nm入射波长633 nm光栅层选TiO2高度200 nm占空比0.5。这个结构在零级透射下会表现出明显的偏振转换特性很适合用来验证计算方法。建模时的核心设置我列一下物理场选电磁波频域接口频域研究。横向边界用周期性边界条件模拟无限大周期阵列。COMSOL里周期性边界条件和Floquet周期端口要配合使用。纵向光传播方向边界用PML包裹端口激励设为TE偏振或TM偏振均可。网格方面光栅层内部最大单元尺寸控制在波长/10左右材料界面处加密外推面附近网格不需要过分加密均匀即可。求解后先在二维绘图里看近场分布确认光栅内部网格没有造成明显的锯齿状伪影。这套模型核心是用Floquet端口得到透射系数配合远场域得到透射零级远场。注意这里的零级透射远场并不是从整个周期结构外推的单个源而是等效为一个具有周期相位基准的辐射结果。由于我们关心的是零级方向的偏振远场计算时可以把角度固定在正前方附近。4.2 执行远场计算与导出复分量求解完成之后在“结果”里新建二维远场绘图。远场绘图中设定角度范围时重点观察零级透射方向$\theta0^{\circ}$为中心展开到 $\pm 15^{\circ}$。如果中心方向正好是法线方向COMSOL默认表达式可以直接给出该方向的远场分量。但我要提醒一下二维远场绘图适合看方向图形态不适合做偏振数据提取。我的做法是新建一个一维绘图组把全局变量表达式写成emw.farfield.Ex emw.farfield.Ey扫描参数选择上面的theta或phi然后生成表格数据。这里需要确认表达式的分量是复数且是直角坐标分量。检查方法很简单在“结果→派生值→全局计算”里设置一个单点评估表达式输入real(emw.farfield.Ex)和imag(emw.farfield.Ex)分别算出数值确认有实虚部输出。导出时COMSOL可以将绘图数据保存为文本文件。每一行是一个角度采样点列包含theta、phi、Ex实部、Ex虚部、Ey实部、Ey虚部。这种格式对后续脚本最友好。如果你用LiveLink for Python直接驱动COMSOL还可以在脚本里直接调用model.result().numerical().create()获取数据跳过文件导入环节自动化程度更高。4.3 用Python脚本计算偏振态并绘图导出文件假设叫farfield_data.txt每列依次为theta, phi, Ex_re, Ex_im, Ey_re, Ey_im。处理脚本可以这样写import numpy as np import matplotlib.pyplot as plt data np.loadtxt(farfield_data.txt) theta data[:, 0] Ex data[:, 2] 1j*data[:, 3] Ey data[:, 4] 1j*data[:, 5] S0, S1, S2, S3, psi, chi, DOP polarization_from_farfield(Ex, Ey) plt.figure(figsize(8, 5)) plt.plot(theta, np.degrees(psi), labelorientation angle psi) plt.plot(theta, np.degrees(chi), labelellipticity angle chi) plt.xlabel(theta (deg)) plt.ylabel(angle (deg)) plt.legend() plt.grid(True) plt.show()这个图一出来偏振转换特性就一目了然如果 $\chi$ 从0变成接近 $\pm 45^{\circ}$说明该角度下透射远场已是接近圆偏振如果 $\psi$ 在某角度发生跳变往往对应S1过零附近属于正常的偏振方向翻折现象画图时要注意连续性。我实际做过一个占空比优化扫描占空比从0.3到0.7每个占空比下都计算632.9 nm透射零级远场的 $\chi$。结果在占空比约0.55时$\chi$ 达到约41.5度很接近圆偏振状态。而这只看近场分布根本看不出来必须借助这套远场偏振流程。5. 常见问题与避坑指南5.1 远场外推边界与PML的配合问题这是最容易出问题也最容易被忽视的点。如果外推面选取的是PML外部边界而不是PML内部边界远场结果极可能出现系统性偏差。原因是PML本身会吸收并残留少量数值反射把外推面放在PML外部等于把这些不干净的“伪远场”也算了进去。我建议外推面放在PML内侧边界上而且计算域形状尽量贴近结构的外接矩形或球形。如果结构是长条形外推面用扁长方体可能会让某个方向的远场精度变差这时候宁可把计算域做大一点也要保证外推面形状规则。另外PML厚度不要省。常见设置为8到12层网格每层使用映射网格。太薄的PML在中低频段也许够用但在光学波段、高入射角情况下反射系数可能显著超标远场偏振也会跟着漂移。5.2 坐标系与相位基准的坑远场偏振计算中最隐蔽的错误来自坐标系定义。COMSOL默认的球坐标远场分量Efarphi、Efartheta是基于每个观察方向局部定义的而偏振的Jones矢量一定需要一个全局统一的直角基矢。如果你直接用泰勒展开里看到的球坐标分量去做Jones矢量在theta0附近phi方向定义会退化结果完全没意义。正确的做法是强制把观察方向对应的电场投影到全局x、y基矢上计算。这也是为什么前文建议直接提取emw.farfield.Ex、emw.farfield.Ey而不是球坐标分量的原因。另一个隐蔽问题是相位基准。如果模型里存在周期性边界条件和Floquet端口COMSOL默认的相位零点在端口面上远场相位会参考该端口。如果你期望的相位基准是某个结构中心点可能需要在后处理里对数据做一次相位平移。数值上很简单将远场复分量乘以 $e^{-j k \Delta r}$其中 $\Delta r$ 是基准点沿观察方向的投影距离。理论清楚这个操作才能解释为什么有些偏振椭圆方向角跟文献对不上。5.3 数据导出与脚本的对接细节使用LiveLink for MATLAB或Python控制COMSOL时有两类方法可以拿远场数据。一种是直接在COMSOL内通过“全局计算”节点得到数值表格然后用LiveLink的model.result().export()接口导出文件另一种是绕过文件直接用model.evaluate()类接口在脚本内取数组。后者少一步文件读写速度快很多适合大批量参数扫描但对COMSOL API的熟悉程度要求更高。我在COMSOL 6.4上测试比较多6.4的Java API和旧版本相比export模块的接口有一定调整。如果你是从6.3或6.2升级过来的老脚本建议先在官方文档里确认一下接口名字是否仍有效否则很容易报“找不到标识符”的错误。还有一个很实用的小技巧不要把所有偏振计算都塞进COMSOL后处理的全局计算表达式里也不要一开始就用LiveLink跑完整脚本。先手动在一个典型模型上导出数据用自己写的Python脚本算一遍确认物理结果合理后再考虑自动化封装。这样能节省大量调试时间。5.4 一改结构就“翻车”的网格稳定性问题亚波长光栅、超表面这类结构一旦扫描某个几何参数网格就跟着重新剖分远场结果可能抖动。这个抖动未必是物理的而是网格变化导致的数值相位偏移。我做参数扫描时会刻意把光栅层网格划分为固定数量的映射网格层数比如高度方向固定分6层、周期方向固定每周期40个单元保证扫描过程中网格拓扑不发生突变。这类问题在远场偏振计算中特别明显因为偏振对相位差的要求比对单一幅度曲线高得多。我的判断标准是把网格加密一倍如果 $\psi$ 和 $\chi$ 的变化小于0.1度当前网格就算及格如果变化大于1度必须继续加密。这个方法简单粗暴但很有效。6. 适用范围与后续扩展这套方法虽然用光栅案例来演示但它完全可以直接迁移到其他场景。射频天线阵列可以分析每个观察方向辐射波的偏振纯度和轴比。等离子体纳米天线研究单个金纳米球的散射远场偏振特别是手性结构中圆二色性信号的仿真。超表面偏振控制器加角度扫描之后可以得到“偏振转换效率-入射角”图谱这是很多光学实验设计需要的数据。波导耦合器辐射只要问题能转成外部辐射问题远场偏振计算流程同样适用。后续如果想把自动化程度再提升一截可以结合前面提到的LiveLink for Python把整个建模、求解、远场外推、偏振参数提取封装成一个函数输入材料参数和几何尺寸直接输出Stokes参数曲线。这就是把这套方法从“一次性的计算”变成“可复用工具”的关键一步。我个人在实际操作中的一个体会是远场偏振计算真正难的不是数学也不是COMSOL某个按钮而是“把流程拆开、让每一步都可验证”。先保证近场相位足够准再固定远场外推边界然后把数据处理放到脚本层最后用网格无关性检查来收尾。这套顺序我基本不换每次都能在比较短的调试周期内拿到稳定可靠的偏振结果。最后再分享一个小技巧如果输出数据里Stokes参数算出来DOP明显小于1比如小于0.99先不要怀疑偏振模型回去检查计算域外推面附近是否有反射波进场多半是边界吸波没做好。把这一点理顺远场偏振计算你就能跑得又快又稳了。
从理论到实操：COMSOL远场偏振计算完整指南
最近不少同行来问我同一个问题COMSOL里模型跑完了近场的Ex、Ey、Ez都很好看但我想拿到某个观察方向上的远场偏振态比如到底线偏振还是椭圆偏振、方向角多少、圆偏振度多高怎么搞这一下就把人问住了。原因很实在COMSOL的远场后处理默认输出的是电场幅值和相位不会自动给你一组干净利落的偏振分量组合而偏振恰恰是远场中那种“两个正交分量的复振幅关系”才能定义出来的东西。我摸索了一套相对通用的远场偏振计算方法从理论图层层拆到COMSOL实操再到导出数据用脚本计算偏振参数前后花了不少时间踩坑。这篇把完整思路和可直接照抄的实现流程都写出来希望能帮有同样需求的人少走弯路。1. 远场偏振计算到底难在哪1.1 COMSOL的远场后处理默认给不了你想要的COMSOL电磁波频域接口里加入远场域Far-Field Domain之后计算完可以很方便地画远场辐射方向图查看远场模值emw.farfield.Efar或者归一化辐射方向。但要注意这个模值是标量它丢掉了两个正交分量的相位关系。而偏振描述的核心恰恰是相位差同样是两个等幅度的电场分量相位差为0时是线偏振相位差为$\pi/2$时就变成圆偏振。只看模值根本区分不了这两者。所以要算偏振第一件事就是改变思路不推荐直接在COMSUL图形界面上看远场图而是把远场电场分量作为数据提取出来用表达式或者外部脚本处理。远场计算本质上是一种外推近场求解完成后在远场域内部或边界上把等效电流、等效磁流通过格林函数积分外推到任意观察角度。这个过程在COMSOL里已经内置完成但它给出的是远处球面上每个观察点的复电场分量一般可访问为emw.farfield.Ex、emw.farfield.Ey、emw.farfield.Ez不同版本变量名有差异老版本常见emw.Efarx这种写法建议先在后处理表达式里验证一下再批量导出。这个复分量信息已经是计算偏振的全部原料了只是很少有人去把它展开成完整的数据流水线。1.2 偏振不是“横竖都能算”相位精度必须先过关有些朋友会想我直接看Ex的实部和Ey的实部比值不就是偏振方向吗但如果场是椭圆偏振相位差在中间状态只取实部显然会得到错误的“短轴”投影。偏振计算需要的是复振幅而且这个复振幅的相位必须足够稳定。这带来两个隐藏要求近场求解的网格精度必须吃得消相位误差。偏振对相位是否精确非常敏感尤其是结构尺度接近波长时网格稍微粗一点远场相位可以偏掉好几度偏振度马上就变了。远场外推所用积分边界的位置也会影响相位基准。COMSOL里远场域通常依赖PML完美匹配层或散射边界包裹外推边界离结构太近会截断倏逝波太远又会增加计算域。实操上我一般建议边界离源区至少半个波长同时确保该边界上的网格和背后PML过渡层连续。也就是说算偏振的成败并不是从后处理开始而是在建模阶段就已经埋下伏笔。这一点后面第5节还会细说。2. 理论层远场外推与偏振参数的数学地图2.1 近场到远场变换本质是“面天线等效源”的叠加COMSOL远场域所做的事情可以类比成天线理论里面经典的面等效源法把一个包围着散射体/辐射源的闭合面看成新的源面上每个位置都有一组等效电流和等效磁流这些等效源辐射出去在远方某个观察方向上对远场电场的贡献满足[ \mathbf{E}{far}(\theta,\phi) \propto \frac{e^{jk r}}{r} \int{\Sigma} \left[ j k (\mathbf{J} \times \hat{r}) \times \hat{r} j k \mathbf{M} \times \hat{r} \right] e^{j k \hat{r} \cdot \mathbf{r}} , dS ]公式看着吓人实际上意思很朴素把近场面上每个点的等效源做一个加权叠加权重由观察方向与源位置的相位差决定。远场就是近场信息经过一个“傅里叶变换式”的传播映射。落到COMSOL操作上你不需要自己实现这个积分只要做两步在电磁波频域接口里添加远场域节点并指定积分面通常直接选整个外部边界。求解后在结果里用远场电场分量表达式来取值。这里的“观察点”在COMSOL远场计算里是一组离散角度的采样。默认情况下它输出的是以原点为参考球的远场球坐标分量但我们做偏振分析时一定要把数据转成全局直角坐标系分量或者直接在直角坐标系下定义观察方向否则后续Jones矢量的构造很容易出错。2.2 偏振的数学工具其实就是三行公式有了远场复电场分量 $E_x$、$E_y$以及需要时 $E_z$偏振的全部信息就浓缩在以下几个量上。我习惯固定使用这套流程通用性很强。首先构造Jones矢量[ \mathbf{J} \begin{bmatrix} E_x e^{j\delta_x} \ E_y e^{j\delta_y} \end{bmatrix} ]由于COMSOL给出的分量本身已经是复数$E_x$、$E_y$的模和辐角都已经包含在复值里面直接拿复数操作即可。然后计算Stokes参数。Stokes参数的好处在于它是实数且与坐标系旋转有明确变换关系适合做后续可视化[ S_0 |E_x|^2 |E_y|^2 ] [ S_1 |E_x|^2 - |E_y|^2 ] [ S_2 2 \mathrm{Re}(E_x E_y^{}) ] [ S_3 2 \mathrm{Im}(E_x E_y^{}) ]接着可以定义两个直观的几何量。第一个是偏振椭圆的长轴方位角 $\psi$代表椭圆长轴相对坐标系x轴的偏转角度[ \psi \frac{1}{2} \mathrm{atan2}(S_2, S_1) ]$\psi$的取值范围一般在 $-90^{\circ}$ 到 $90^{\circ}$ 之间。第二个是椭圆率角 $\chi$表示椭圆“胖瘦”取正值时代表右旋圆偏振成分占优取负值代表左旋[ \chi \frac{1}{2} \mathrm{atan2}\left(S_3, \sqrt{S_1^2 S_2^2}\right) ]$\chi$的取值范围在 $-45^{\circ}$ 到 $45^{\circ}$ 之间。当 $\chi \pm 45^{\circ}$ 时就是纯圆偏振当 $\chi 0$ 时就是纯线偏振。用生活类比来说$\psi$ 是椭圆在纸上画出的那条主轴朝哪边倾斜$\chi$ 是这张椭圆变形得“扁不扁”、旋转方向如何。单一角度下如果场不是完全偏振的比如包含随机非相干背景还可以计算偏振度[ DOP \frac{\sqrt{S_1^2 S_2^2 S_3^2}}{S_0} ]对一个完全偏振的理想电磁波DOP等于1。实际仿真中如果计算域边界处理不好远场混入非物理数值噪声DOP会略小于1这个指标也可以反过来用来检查远场外推做得好不好。3. 通用计算方法从模型设计到数据后处理的完整流水线3.1 建模阶段就要为远场偏振留好口子我把这套方法定位成“通用”就是因为它和具体结构无关——不管是纳米天线、光栅、超表面还是射频天线阵列只要COMSOL频域求解能给出近场这套流程都能用。前提是建模阶段满足几个基本要求使用电磁波频域接口而且求解的是完整频点不是本征模。计算域外围加载PML或足够可靠的吸波边界。远场计算依赖外推面如果外推面外面还有较强的反射波远场方向图会被污染。在物理场设置中显式添加远场域并把外推面选在PML层之前的那一圈内部边界上。结果保存时勾选“所有场的解”或者足够密的角度采样否则导出远场数据时角度分辨率不足后续偏振扫描就没法看。通用流程我整理成了四步频域求解近场。基于近场结果做远场外推得到多个角度下的复电场分量。把复电场分量按直角坐标导出成结构化数据如文本文件或CSV。用Python或MATLAB脚本读入数据计算Stokes参数与偏振椭圆参数。这个流水线的好处是解耦。COMSOL内做的事情只管出“准”的场数据偏振计算全部放在脚本端换模型时不需要反复修改脚本逻辑只需要改数据列映射关系。3.2 远场域的关键配置项添加远场域节点后有几个参数要特别留意积分面积分的边界通常选择包围结构的最外层计算域边界也就是PML内部的边界。直接选整个计算域外表面是没问题的。参考波矢方向即默认的平面波入射方向它对远场的相位基准有影响。如果不是垂直入射要注意观察角度的定义。输出球坐标还是直角坐标COMSOL后处理里远场结果既有球坐标分量如Efarphi、Efartheta也有直角坐标分量。偏振分析我建议固定用直角坐标分量也就是用表达式把emw.farfield.Ex、emw.farfield.Ey取出来。因为球坐标分量在不同谱面的定义跟观察方向绑定后续写通用脚本时非常容易搞混。角度采样通常后处理远场就可以设置theta和phi的扫描点数。默认可能只有几十个点做偏振角度谱时会显得太粗建议至少设置到200到400个点。数据导出这一步可以直接在COMSOL“结果”里用一维绘图组以theta或者phi作为扫描变量把远场电场分量画成线图然后将绘图数据导出为文本。也可以用“全局计算”节点用一个表格列出多个表达式的值。两种方式各有利弊我后文会聊更推荐的方案。3.3 指标体系一个脚本吃透全部偏振参数写脚本时最重要的是把Stokes参数计算封装成一个可复用的函数。以下给出一个清醒的函数骨架Python风格理论上可用于任何远端数据import numpy as np def polarization_from_farfield(Ex, Ey): 输入复数远场电场分量 Ex, Ey直角坐标系 返回 S0~S3 与偏振椭圆参数 psi, chi, DOP Ex np.asarray(Ex, dtypecomplex) Ey np.asarray(Ey, dtypecomplex) S0 np.abs(Ex)**2 np.abs(Ey)**2 S1 np.abs(Ex)**2 - np.abs(Ey)**2 S2 2 * np.real(Ex * np.conj(Ey)) S3 2 * np.imag(Ex * np.conj(Ey)) psi 0.5 * np.arctan2(S2, S1) # 椭圆长轴方位角 chi 0.5 * np.arctan2(S3, np.sqrt(S1**2 S2**2)) # 椭圆率角 DOP np.sqrt(S1**2 S2**2 S3**2) / np.maximum(S0, 1e-12) return S0, S1, S2, S3, psi, chi, DOP这段代码读取的是复值数组如果你的导出数据是分实部、虚部列出进函数之前先拼成复数Ex Ex_real 1j * Ex_imag Ey Ey_real 1j * Ey_imag这样整个计算流程就只剩下“把COMSOL的数据喂给这个函数”这一步了。我实际使用时还会顺手把角度数据一并输出这样就能画出“偏振椭圆参数随方位角变化”或者“随波长变化”的曲线。4. 实操案例亚波长光栅零级透射远场偏振分析4.1 模型搭建与关键参数用一个非常常见的场景来演示可见光波段二氧化硅基底上做周期性亚波长介质光栅正入射平面波计算零级透射远场的偏振态。光栅周期设为600 nm入射波长633 nm光栅层选TiO2高度200 nm占空比0.5。这个结构在零级透射下会表现出明显的偏振转换特性很适合用来验证计算方法。建模时的核心设置我列一下物理场选电磁波频域接口频域研究。横向边界用周期性边界条件模拟无限大周期阵列。COMSOL里周期性边界条件和Floquet周期端口要配合使用。纵向光传播方向边界用PML包裹端口激励设为TE偏振或TM偏振均可。网格方面光栅层内部最大单元尺寸控制在波长/10左右材料界面处加密外推面附近网格不需要过分加密均匀即可。求解后先在二维绘图里看近场分布确认光栅内部网格没有造成明显的锯齿状伪影。这套模型核心是用Floquet端口得到透射系数配合远场域得到透射零级远场。注意这里的零级透射远场并不是从整个周期结构外推的单个源而是等效为一个具有周期相位基准的辐射结果。由于我们关心的是零级方向的偏振远场计算时可以把角度固定在正前方附近。4.2 执行远场计算与导出复分量求解完成之后在“结果”里新建二维远场绘图。远场绘图中设定角度范围时重点观察零级透射方向$\theta0^{\circ}$为中心展开到 $\pm 15^{\circ}$。如果中心方向正好是法线方向COMSOL默认表达式可以直接给出该方向的远场分量。但我要提醒一下二维远场绘图适合看方向图形态不适合做偏振数据提取。我的做法是新建一个一维绘图组把全局变量表达式写成emw.farfield.Ex emw.farfield.Ey扫描参数选择上面的theta或phi然后生成表格数据。这里需要确认表达式的分量是复数且是直角坐标分量。检查方法很简单在“结果→派生值→全局计算”里设置一个单点评估表达式输入real(emw.farfield.Ex)和imag(emw.farfield.Ex)分别算出数值确认有实虚部输出。导出时COMSOL可以将绘图数据保存为文本文件。每一行是一个角度采样点列包含theta、phi、Ex实部、Ex虚部、Ey实部、Ey虚部。这种格式对后续脚本最友好。如果你用LiveLink for Python直接驱动COMSOL还可以在脚本里直接调用model.result().numerical().create()获取数据跳过文件导入环节自动化程度更高。4.3 用Python脚本计算偏振态并绘图导出文件假设叫farfield_data.txt每列依次为theta, phi, Ex_re, Ex_im, Ey_re, Ey_im。处理脚本可以这样写import numpy as np import matplotlib.pyplot as plt data np.loadtxt(farfield_data.txt) theta data[:, 0] Ex data[:, 2] 1j*data[:, 3] Ey data[:, 4] 1j*data[:, 5] S0, S1, S2, S3, psi, chi, DOP polarization_from_farfield(Ex, Ey) plt.figure(figsize(8, 5)) plt.plot(theta, np.degrees(psi), labelorientation angle psi) plt.plot(theta, np.degrees(chi), labelellipticity angle chi) plt.xlabel(theta (deg)) plt.ylabel(angle (deg)) plt.legend() plt.grid(True) plt.show()这个图一出来偏振转换特性就一目了然如果 $\chi$ 从0变成接近 $\pm 45^{\circ}$说明该角度下透射远场已是接近圆偏振如果 $\psi$ 在某角度发生跳变往往对应S1过零附近属于正常的偏振方向翻折现象画图时要注意连续性。我实际做过一个占空比优化扫描占空比从0.3到0.7每个占空比下都计算632.9 nm透射零级远场的 $\chi$。结果在占空比约0.55时$\chi$ 达到约41.5度很接近圆偏振状态。而这只看近场分布根本看不出来必须借助这套远场偏振流程。5. 常见问题与避坑指南5.1 远场外推边界与PML的配合问题这是最容易出问题也最容易被忽视的点。如果外推面选取的是PML外部边界而不是PML内部边界远场结果极可能出现系统性偏差。原因是PML本身会吸收并残留少量数值反射把外推面放在PML外部等于把这些不干净的“伪远场”也算了进去。我建议外推面放在PML内侧边界上而且计算域形状尽量贴近结构的外接矩形或球形。如果结构是长条形外推面用扁长方体可能会让某个方向的远场精度变差这时候宁可把计算域做大一点也要保证外推面形状规则。另外PML厚度不要省。常见设置为8到12层网格每层使用映射网格。太薄的PML在中低频段也许够用但在光学波段、高入射角情况下反射系数可能显著超标远场偏振也会跟着漂移。5.2 坐标系与相位基准的坑远场偏振计算中最隐蔽的错误来自坐标系定义。COMSOL默认的球坐标远场分量Efarphi、Efartheta是基于每个观察方向局部定义的而偏振的Jones矢量一定需要一个全局统一的直角基矢。如果你直接用泰勒展开里看到的球坐标分量去做Jones矢量在theta0附近phi方向定义会退化结果完全没意义。正确的做法是强制把观察方向对应的电场投影到全局x、y基矢上计算。这也是为什么前文建议直接提取emw.farfield.Ex、emw.farfield.Ey而不是球坐标分量的原因。另一个隐蔽问题是相位基准。如果模型里存在周期性边界条件和Floquet端口COMSOL默认的相位零点在端口面上远场相位会参考该端口。如果你期望的相位基准是某个结构中心点可能需要在后处理里对数据做一次相位平移。数值上很简单将远场复分量乘以 $e^{-j k \Delta r}$其中 $\Delta r$ 是基准点沿观察方向的投影距离。理论清楚这个操作才能解释为什么有些偏振椭圆方向角跟文献对不上。5.3 数据导出与脚本的对接细节使用LiveLink for MATLAB或Python控制COMSOL时有两类方法可以拿远场数据。一种是直接在COMSOL内通过“全局计算”节点得到数值表格然后用LiveLink的model.result().export()接口导出文件另一种是绕过文件直接用model.evaluate()类接口在脚本内取数组。后者少一步文件读写速度快很多适合大批量参数扫描但对COMSOL API的熟悉程度要求更高。我在COMSOL 6.4上测试比较多6.4的Java API和旧版本相比export模块的接口有一定调整。如果你是从6.3或6.2升级过来的老脚本建议先在官方文档里确认一下接口名字是否仍有效否则很容易报“找不到标识符”的错误。还有一个很实用的小技巧不要把所有偏振计算都塞进COMSOL后处理的全局计算表达式里也不要一开始就用LiveLink跑完整脚本。先手动在一个典型模型上导出数据用自己写的Python脚本算一遍确认物理结果合理后再考虑自动化封装。这样能节省大量调试时间。5.4 一改结构就“翻车”的网格稳定性问题亚波长光栅、超表面这类结构一旦扫描某个几何参数网格就跟着重新剖分远场结果可能抖动。这个抖动未必是物理的而是网格变化导致的数值相位偏移。我做参数扫描时会刻意把光栅层网格划分为固定数量的映射网格层数比如高度方向固定分6层、周期方向固定每周期40个单元保证扫描过程中网格拓扑不发生突变。这类问题在远场偏振计算中特别明显因为偏振对相位差的要求比对单一幅度曲线高得多。我的判断标准是把网格加密一倍如果 $\psi$ 和 $\chi$ 的变化小于0.1度当前网格就算及格如果变化大于1度必须继续加密。这个方法简单粗暴但很有效。6. 适用范围与后续扩展这套方法虽然用光栅案例来演示但它完全可以直接迁移到其他场景。射频天线阵列可以分析每个观察方向辐射波的偏振纯度和轴比。等离子体纳米天线研究单个金纳米球的散射远场偏振特别是手性结构中圆二色性信号的仿真。超表面偏振控制器加角度扫描之后可以得到“偏振转换效率-入射角”图谱这是很多光学实验设计需要的数据。波导耦合器辐射只要问题能转成外部辐射问题远场偏振计算流程同样适用。后续如果想把自动化程度再提升一截可以结合前面提到的LiveLink for Python把整个建模、求解、远场外推、偏振参数提取封装成一个函数输入材料参数和几何尺寸直接输出Stokes参数曲线。这就是把这套方法从“一次性的计算”变成“可复用工具”的关键一步。我个人在实际操作中的一个体会是远场偏振计算真正难的不是数学也不是COMSOL某个按钮而是“把流程拆开、让每一步都可验证”。先保证近场相位足够准再固定远场外推边界然后把数据处理放到脚本层最后用网格无关性检查来收尾。这套顺序我基本不换每次都能在比较短的调试周期内拿到稳定可靠的偏振结果。最后再分享一个小技巧如果输出数据里Stokes参数算出来DOP明显小于1比如小于0.99先不要怀疑偏振模型回去检查计算域外推面附近是否有反射波进场多半是边界吸波没做好。把这一点理顺远场偏振计算你就能跑得又快又稳了。