扭转光子晶体远场偏振的COMSOL仿真与出图实战

光子晶体这个方向我做了不少年仿真各种数值工具也摸过一圈。最近一个项目里需要研究扭转光子晶体twisted photonic crystal对远场偏振的调控并且最后要直接出能放论文的图。试了一圈工具之后我在COMSOL里把几何建模、参数化扫描、远场投影计算和后处理出图全部打通从EFARx/EFARy两个正交分量到三维辐射方向图一条线做下来几乎没有再导出数据到别的软件。如果你也是做微纳光子、超表面、手性光学或者偏振调制器件仿真的人这篇文章应该能帮你把“计算远场偏振并直接出图”这条路走得顺一些。1. 从科学问题到仿真思路扭转光子晶体为什么要盯住远场偏振1.1 扭转光子晶体在做什么普通光子晶体靠周期性排列的介电结构形成能带和光子带隙。扭转光子晶体则是两层或更多层周期性结构之间相对旋转一个角度这个自由度在早期研究里经常被忽略因为大多数人做光子晶体都默认结构是严格周期对齐的。直到后来有实验发现上下两层光栅或者孔阵列相对扭转之后透射谱和偏振特性会发生明显变化这个方向才热起来。扭转角在几何上并不复杂就是上层结构绕某个轴旋转θ但它给系统带来的物理后果却很深。比如原本对称的电磁模式会因为扭转而发生耦合不同偏振分量的传播常数产生差异宏观上就表现为圆双折射、偏振旋转或者圆二色性。我用COMSOL做这类结构时最关心的不是能带有多干净而是这个扭转角能否改变远场辐射的偏振态。因为对器件来说能带只是基础真正进入自由空间的是远场探测器和接收端看的也是远场。这类结构可以用在一个很实在的场景设计一个紧凑的偏振转换器或者手性光发射器。传统的偏振调节往往靠半波片、四分之一波片加上机械旋转体积大且波长依赖强。扭转光子晶体是平面器件只需要在加工时控制两层之间的扭转角就能把出射光的偏振从线偏振调成椭圆偏振甚至圆偏振。这种思路在超表面和微纳光子学里已经很常见但真正落实到仿真和出图环节很多细节需要自己一点点踩。1.2 远场偏振的物理与仿真意义远场偏振到底是什么简单说就是在一段距离之外、通常是在几个波长到无穷远的位置上电场矢量末端在垂直于传播方向平面内画出的轨迹。这个轨迹是直线就是线偏振是正圆就是圆偏振介于两者之间就是椭圆偏振。判断偏振态需要两个信息两个正交电场分量的幅度比以及它们之间的相位差。为什么要算远场而不是直接看近场因为近场电场包含大量倏逝波和局域共振场这些场不能传播到远处对实际探测没有贡献。如果你在结构表面附近画电场分布看出一堆涡旋和热点那说明局域场很强但器件对外表现的偏振态必须用远场投影来定。远场的计算本质上是把包围结构的闭合面上的等效电流做辐射积分得到空间中各个方向上的辐射幅度和相位。在COMSOL里这个积分被封装成“远场域”节点求解之后自动给出EFARx、EFARy、EFARz这些复数分量。判断偏振调控效果我习惯盯三组量方向图也就是远场强度在空间中的角度分布正交分量幅度比 |EFARx / EFARy|这个决定偏振椭圆的长短轴比例相位差 δ arg(EFARx) - arg(EFARy)这个决定偏振椭圆的形状和旋向。这三组量全部可以从COMSOL远场结果里直接提取不需要自己写积分也不需要额外导出数据。1.3 为什么选择COMSOL直接出图一开始我也考虑过用RCWA或者FDTD。RCWA算平板光子晶体确实快但处理任意扭转几何和远场投影要自己加代码后处理也不是所见即所得。FDTD的近场到远场变换需要先跑完时域仿真再做变换最后导出数据到Python或者MATLAB里画中间链条很长改一个扭转角就要重复一遍。COMSOL的好处是有限元建模、求解、远场投影、后处理绘图全在一个环境里完成而且参数化扫描之后可以直接生成不同角度下的远场图这对找趋势、做对比非常方便。更重要的是COMSOL的远场绘图组是专门为辐射类问题设计的。二维极坐标图可以画角度分布三维远场图可以画球坐标辐射方向图表达式框里直接填 ef ar.EFarx 这样的变量就能出图。论文里常用的极坐标曲线和辐射球图在COMSOL里几分钟就能导出一张清晰的结果图。当然它也不是没有短板比如偏振椭圆图这种专业工具它没有内置需要自定义表达式或者用参数化曲线来画这部分后面我会详细讲。2. 模型搭建工作平面、几何参数与边界条件的关键细节2.1 几何建模思路我做的模型是一个双层光子晶体平板底层和顶层各有一组周期性排列的空气孔顶层相对底层绕z轴旋转一个角度。为了避免全结构建模导致计算量爆炸我采用单胞建模加周期性边界条件的方案。周期单元包括上下两层介质平板、空气孔以及上下方的空气层。空气层的厚度至少要留出半个波长以上因为远场积分面要放在离结构足够远的位置。具体参数可以按你的器件工作波长来定我这里给一组常用起始值周期 a 800 nm每层平板厚度 t 200 nm空气孔半径 r 200 nm工作波长 λ0 1550 nm背景介质为空气平板材料为硅折射率 n 3.45扭转角 θ 从 0° 扫到 90°关键技巧是把所有尺寸都定义为全局参数。这样后面做扫描的时候只需要在参数化扫描里添加一个角度变量模型就会自动重建。不要直接在几何里写死数值否则每次改几何都要手动修改多处很容易出错。2.2 工作平面和扭转参数化的作用COMSOL的工作平面在三维建模里非常实用。它相当于一个可移动的二维草图画布你可以在任意空间平面上画草图再通过拉伸、旋转等操作变成三维实体。对扭转光子晶体来说工作平面的价值在于它让两层结构的定位变得特别清晰。我的习惯是这样先在全局坐标系的原点处创建一个工作平面画出底层周期单元的背景平板然后在z坐标等于t的位置再创建一个工作平面画顶层平板。画空气孔的时候先在每个工作平面上画一个圆再用“拉伸”操作生成圆柱最后从平板中减去这些圆柱。顶层孔的旋转不要手动改坐标而是用COMSOL的“变换”节点里的“旋转”功能把旋转中心设为单元中心旋转角填全局参数theta_twist。这样做的好处是扫描扭转角时几何会自动跟着变不会出现坐标错乱的问题。提到工作平面顺便说一个易错点如果你添加了多个工作平面一定要在“工作平面”设置里确认基准面是哪一个否则容易出现草图画在错误平面上的情况。检查方法很简单打开几何视图在工作平面节点上右键“启用”看看绘图区域的高亮位置是否和你预期的一致。2.3 材料、边界条件和端口设置材料设置基本是直来直去平板区域给硅的折射率空气孔和空气层给折射率1。如果有衬底也可以加一层二氧化硅但那样模型层数会变多建议先跑一个没有衬底的简化版本验证物理规律再决定要不要加。物理接口选择“电磁波频域”ewfd。单胞建模时四个侧面需要使用周期性边界条件或者如果你算的是透射/反射谱直接用“周期性端口”更合适。上下两个端面用端口边界端口1设为入射端口端口2设为出射端口。入射波可以设为x偏振这样正好便于观察远场里x分量和y分量的比例变化。如果目标不是计算透射反射谱而是分析某个光源照射结构后的远场散射偏振那可以把上下边界替换为散射边界条件或者完美匹配层并在结构附近放一个偶极子源或者高斯束激励。无论哪种场景远场计算都支持。我最常用的做法是上下加PML内部边界的包络面选为远场积分面这样的好处是背向散射也能被完整捕捉远场方向图不会出现截断伪影。边界条件这一步容易出问题的是端口模式设置。对周期性端口你需要指定模式的极化方向。如果端口极化角设置不对入射波的偏振方向就不是你想要的x偏振后面分析偏振比例时就会一头雾水。所以每次建模后我会先用一个不含扭转角的结构跑一遍确认透射光的偏振与入射一致再做扭转扫描。2.4 网格划分与远场域的处理网格是有限元仿真里最影响准确性和资源消耗的一环。光子晶体结构的孔洞边缘曲率较大需要局部加密。我一般把最大单元尺寸设为工作波长的六分之一到八分之一在孔洞边界再加边界层网格。对于1550 nm的波长硅中波长约449 nm空气孔区域可以稍微粗一些但硅内部至少要保证足够的网格密度。远场域的处理方式很关键。远场积分面要选一个包围结构的闭合面并且这个面不能贴在结构表面最好离结构至少半个波长。COMSOL会在这个面上做等效电流积分如果面积取得太小近场的倏逝波成分会直接影响远场结果导致EFARx/EFARy出现明显误差。我通常把空气层高度设置为1000 nm到2000 nm在最外侧选一个长方体盒子的表面作为远场域边界然后在这个边界上用较均匀的网格。网格质量会影响远场相位精度特别是在顶层结构旋转之后倾斜的孔壁会产生一些扭曲单元。网格生成后我习惯打开“统计”检查最小单元质量低于0.1就要考虑局部细化或者重新剖分。相位差对网格误差尤其敏感因为它是两个分量相减得到的小的单元误差可能被放大成几度的相位误差。3. 远场计算与EFARx/EFARy的提取直接把数据变成图3.1 远场计算的原理与启用方式COMSOL的远场计算原理并不神秘。它基于斯特拉顿-朱工程积分把远场积分边界上的近场解当作等效电磁流再通过自由空间格林函数辐射积分得到空间中任意方向的远场电场复数分量。之所以选择这种间接方式是因为直接模拟无穷远处的边界会消耗巨大计算资源而远场近似之后所有角度上的远场结果都可以从一次近场求解中后处理得到。启用方式很简单在“电磁波频域”接口下右键添加“远场域”节点然后选中你想作为积分面的所有边界。有些版本把这个节点叫“远场计算”或者“远场频域”本质一样。添加完之后求解器会自动在求解过程中保存远场数据。需要注意如果模型中没有添加远场域后面的远场绘图组就是空的这是很多人一上来就踩的坑。如果要在整个球面方向上看远场还需要在“远场域”的设置里指定角度分辨率。COMSOL支持设置theta和phi的步长默认值往往比较粗比如每10度一个点画出来的三维方向图会比较粗糙。论文出图时我一般把步长设为2度或者1度但角度点数增加会显著增加后处理计算量可以先粗扫定位再细扫出图。3.2 表达式选择EFARx、EFARy与相位信息求解完成之后远场数据会出现在结果的数据集里变量名通常是 ef ar.EFarx、ef ar.EFary、ef ar.EFarz还有一个表示总场的 ef ar.EFar。这些变量都是复数实部虚部分别对应电场分量的相位信息。要画幅度信息表达式用 abs(ef ar.EFarx) 或 abs(ef ar.EFary)要画辐照度可以用 ewfd.normEfar或者直接对总场取平方。我的建议是先画出各分量的幅度分布再看相位差。相位差在COMSOL里不是一个现成的量需要用 arg(ef ar.EFarx) - arg(ef ar.EFary) 这样的表达式自己定义。注意复数幅角函数返回的是主值范围在负π到π之间如果数据点在正负π附近跳变画图会出现不连续跳跃需要做展开处理。我一般会在后处理表达式里加一个判断把相位差连续化或者干脆导出数据后做unwrap。不过如果只是看趋势不连续点通常是可接受的。3.3 直接出图的几种手段“直接出图”是这篇文章的核心所以我把COMSOL里几种最适合远场偏振展示的出图方式列出来二维极坐标图。这是最常用的方式。新建绘图组时选择“二维极坐标”数据集选远场数据集然后在“极坐标图”设置里把角度坐标设为theta表达式填 |EFARx| 和 |EFARy|就能看到两个正交偏振分量在空间不同方向上的强度分布。这个图可以直接看出偏振分量的空间分离情况是判断器件偏振调控能力的第一张图。三维远场辐射图。当你需要展示整个半球空间的辐射强度时用三维绘图组添加远场绘图节点表达式用总场辐照度颜色表选“彩虹”网格分辨率调高就能得到漂亮的辐射方向球图。配合透明度设置可以看到内部辐射瓣的形状。偏振椭圆参数化曲线。COMSOL没有内置“偏振椭圆”这个节点但可以通过1D绘图组里的“参数化曲线”来实现。思路是以时间t为参数在某个固定方向(theta0, phi0)上画出 x Re(EFARx * exp(j * omega * t))y Re(EFARy * exp(j * omega * t)) 的利萨如图。这样画出来的就是一个偏振椭圆旋向和长短轴比例一目了然。这个操作稍微有点绕但确实是纯COMSOL环境的解法。全局计算加1D曲线。如果只是看某个远场方向上相位差随扭转角的变化可以先用“派生值”里的“全局计算”提取目标角度下的EFARx和EFARy再绘制成1D曲线。这种方法很适合做参数扫描后的趋势分析。出图时还有一个和论文相关的设置在绘图窗口里可以调整坐标轴范围、字体大小和颜色表导出图片时选择合适的分辨率。我一般导出PNG分辨率设300 DPI尺寸按论文栏宽设置这样直接插入LaTeX就能用不用再到别的软件里二次裁剪。4. 参数化扫描实战扭转角如何影响远场偏振4.1 参数化扫描的配置现在到了最能体现COMSOL优势的环节。在研究中新建一个“参数化扫描”步骤把扭转角设为扫描参数比如写 range(0,15,90)意思是从0度开始每15度采样一个点一直到90度。如果你希望更精细也可以改成 range(0,5,90)但要注意求解次数和计算时间会成倍增加。这里有个关键设置扫描得到的解要保存为“全部解”而不是只保存最后一步。否则后处理时只能在最后一个角度下画图无法对比不同角度的结果。默认情况下求解器有时只保存当前参数值对应的解需要在研究设置里勾选“存储所有参数解”。我建议先做单频点扫描。因为光子晶体器件往往有较强的色散如果同时扫描波长和角度数据量会非常大而且很难快速判断趋势。先固定在一两个关键波长上扫描角度找到规律后再做波长扫描。4.2 后处理中实现多曲线对比参数化扫描完成后绘图时数据集选择“参数化解”而不是“解1”。这样在绘图设置里会多出一个“参数”选择项可以通过“保留”参数或者“生成图例”的方式把不同扭转角的曲线画在同一张图里。最佳操作方式是在二维极坐标图的数据设置里选择“所有参数值”这样每条曲线分别对应一个扭转角图例会显示角度值。线条太密时可以只选几个关键角度比如0°、15°、30°、45°、60°、90°用不同线型和颜色区分。如果要画相位差随扭转角的变化建议走“全局计算”的路子。先为某个固定远场方向定义一组求解表达式再在1D绘图组里横轴设为 θ纵轴设为相位差这样得到的就是一条干净的曲线。这张图在论文里往往是最有说服力的因为它直接证明了扭转角与偏振态之间的定量关系。4.3 偏振态判读与结果分析为了帮助理解我列一组示意性的趋势数据不是某个具体模型的真实结果真实结果需要自己跑| 扭转角 θ | |EFARx/EFARy| | 相位差 Δ | |---|---|---| | 0° | 1.00 | 0° | | 15° | 0.95 | 12° | | 30° | 0.85 | 35° | | 45° | 0.80 | 62° | | 60° | 0.72 | 78° | | 75° | 0.68 | 85° | | 90° | 0.65 | 90° |从这个示意趋势可以看出当扭转角为0时两个正交分量幅度比接近1相位差接近0远场是线偏振随着扭转角增大相位差逐渐增大幅度比开始偏离1偏振态逐步从线偏振过渡到椭圆偏振当扭转角达到90度附近相位差趋于90度结构表现出更明显的圆偏振特性。这个变化趋势说明扭转角确实可以作为一个几何自由度来调控远场偏振。当然实际仿真结果不可能像我这张示意表这么漂亮真实数据里还会出现方向图旁瓣、相位差振荡等细节。你需要盯住的是主辐射方向上的数值趋势以及相位差变化是否单调。如果出现跳变先检查相位unwrap问题再检查网格收敛性最后再判断是不是物理上真的存在模式切换。5. 常见问题与排查陷阱绘图为空、导入警告、网格内存优化5.1 远场绘图为空怎么排查“绘图为空”这个问题在COMSOL论坛里几乎每周都有人问尤其是在远场出图时。我遇到的情况基本可以归为四类第一远场域节点没有启用。如果模型里压根没有添“远场域”求解时就不会产生远场数据绘图组自然什么都没有。应对方法回到物理接口添加远场域节点选中正确的边界重新求解。第二数据集选错了。远场数据有时会在结果里生成一个独立的数据集比如叫“Study 1/Solution 1 (far field)”你在绘图组里却选了“Study 1/Solution 1”这样即使表达式正确也是空的。第一时间检查数据集名称。第三绘图组类型不对。远场角度扫描要用“二维极坐标”绘图组如果用了普通的1D绘图组需要自己指定角度数组很容易配置错误。建议直接用极坐标绘图组。第四表达式里引用变量名错误。不同版本变量的名称有差异有的是 ef ar.EFarx有的是 ewfd.EFarx。最快的办法是打开“派生值”里的“计算远场”对话框里面有自动插入表达式的按钮直接用它生成的表达式不要手敲变量名。5.2 外部几何导入警示SolidWorks转STEP的常见坑很多同行习惯在SolidWorks里画好结构再导入COMSOL结果一导入就弹出一大堆警告。我处理过的案例里最有代表性的警告是“几何操作中检测到缺陷”和“边/面参数化失败”。原因通常是SolidWorks模型里有圆角、倒角、细小台阶这些特征导出STP后这些特征被描述成复杂的曲面片在COMSOL内核重建时容易出问题。光子晶体结构本身基本都是圆柱和方块反而最适合直接在COMSOL里画。如果你一定要走导入路线我的经验是在SolidWorks里另存为STEP之前先把所有对仿真无关紧要的圆角、倒角删除导出时选择“实体”而不是“装配体”避免出现多余坐标系和实例COMSOL导入时如果警告不少可以在导入设置里适当调整修复容差比如从默认的1e-6改成1e-5但不要改太大否则细小孔洞特征可能被吞掉导入后检查一下几何对象把多余的曲面对象删除只保留实体或表面。说实话对于扭转光子晶体这种参数化结构直接在COMSOL里用工作平面和旋转操作重建一遍比手动导入再修几何快得多也更容易做参数化扫描。5.3 网格与内存优化心得三维光子晶体的网格数量很容易失控。一个单胞模型如果网格过细自由度可能到几百万甚至上千万普通工作站跑起来很吃力。我的优化顺序是这样先检查能否用单胞代替全结构。这是最有效的一步。只要有周期性四个侧面用周期性边界条件计算量可以小一个数量级。网格尺寸按波长评估而不是盲目加密。每波长6个二阶单元通常能给出可以接受的结果先把粗网格跑通再在关键区域加密而不是一开始就追求全局高精度。求解器选择上默认的UMFPACK对三维问题内存占用很高我通常换成PARDISO它支持多核并行内存管理更好。如果自由度特别大可以试迭代求解器配预条件。不过远场计算这种问题我还是建议先控制自由度而不是依赖更强的求解器。还有一个容易被忽略的点扭转结构产生的倾斜网格面会降低单元质量。建议在旋转孔洞的周围加边界层网格同时开启网格质量直方图检查。单元质量低于0.1会导致求解收敛变慢甚至不收敛别等到算到一半才回头改网格。另外远场点的数量也直接影响后处理内存。三维方向图如果phi和theta都设成1度步长会有64800个远场点每个点都有三个复数分量数据量不小。建议先用粗糙角度步长扫描确定感兴趣的方向后再单独细化这些方向。我在实际跑这个课题时最后发现最花时间的其实不是求解而是把扭转角和偏振椭圆参数之间的关系整理成清晰图表。如果一开始就把参数化扫描和后处理模板搭好后面改参数、改材料、改入射条件都会非常省事。一个建议是把自己常用的远场出图配置保存成单独的模型模板下次新建模型直接套用这样就能把精力集中在物理分析上而不是反复折腾后处理设置。