做这个课题的人我猜多半是在超表面方向被“相位没法连续调”卡住或者是在液晶显示/液晶光子学这边被“器件尺寸太大、调相速度太慢”困住最后两头一碰决定仿一个液晶加载超表面单元。这个方向确实是这几年可编程超表面、光束偏转、动态波前调控里最热的一条路也是论文里最容易出功能演示图的题材。但说实话COMSOL里的“张量矩阵设置”这一关第一次接触的人没有不懵的——液晶不是普通材料它不能只填一个折射率完事你得让介电张量跟着分子指向转而且不同位置的液晶方向还可能不一样。这篇文章我就从底层讲清楚为什么必须要用张量、张量怎么从分子指向矢组装出来、任意液晶分布怎么塞进COMSOL、以及最终怎么把一个液晶-超表面单元跑通拿到相位曲线。适合正在做可调超表面、液晶光子学、COMSOL光学仿真方向的研究生和工程师参考也适合想复现这类论文但又不想在材料设置上卡一个月的初学者。1. 为什么液晶和超表面要放进同一个模型里算1.1 超表面相位调制的两堵墙先泼一盆冷水纯无源超表面做相位调制相位分布一旦加工完成就焊死在芯片上了。你可以在设计阶段用几何相位、传播相位或者耦合模理论把0到2π的相位铺满但做完之后想动态改波前对不起做不到。你只能做固定功能的透镜、分束器、偏转器。这是第一堵墙。第二堵墙是调相的速度和灰度。很多人想到用MEMS、相变材料、热光效应去调但每个方案都有硬伤MEMS有可动结构可靠性存疑相变材料如GST有阈值能耗问题连续灰阶不容易控制热光效应响应时间慢到毫秒级。液晶在这时候反而成了综合性价比最好的选项——它的折射率各向异性大双折射Δn普遍在0.1~0.2以上加上电场调制的响应可以做到微秒到毫秒级而且连续可调工艺也非常成熟整个显示产业都在用代工和材料都不难找。所以把液晶盒直接叠在超表面上用液晶的可变双折射去连续调制超表面单元的谐振相位就成了很自然的“取长补短”方案。你设计好的超表面负责把相位铺满0~2π液晶层负责给这个相位加上一个连续的、可按电压切换的偏移量。仿真上这个结构就是一个典型的多层各向异性电磁问题必须在COMSOL里把液晶、超表面、基底一起算。1.2 一起仿真真正的难点在哪这个模型不是“两个仿真拼一起”那么简单。难点有三个第一是尺度跨度。超表面单元周期一般是几百纳米谐振结构长宽几百纳米、高度一百到三百纳米液晶盒厚度则是几微米比超表面的尺度大一个数量级。网格划分的时候液晶层不能太粗至少要两三层网格去分辨相位累积超表面本身又要在表面处加密否则谐振算不准。这两个需求叠加模型自由度就不小3D全波仿真要控制好网格数。第二是材料属性不是标量。超表面金属金或银的介电常数是复数并且随波长变化液晶的介电常数是一个张量方向由分子指向决定。你在模型中填的每一个材料参数最终都会影响谐振频率和相位。尤其是液晶只填一个折射率必然算错——光在一个方向看到的折射率和在另一个方向看到的不一样。必须用完整的相对介电常数张量去表示。第三是液晶分布可能处处不一样。出厂时你可以在基板上做配向让液晶分子沿某个方向排列加电压后分子转向而且不是所有位置同步同角度转向存在弹性形变的连续分布。如果只做“均匀液晶层”假设那简单但也意味着你放弃了那些利用液晶空间分布调制的更高级玩法比如液晶衍射光栅、螺旋分布的操控等。要做好“任意液晶分布”核心就是把你计算或假设的指向矢分布nx(x,y,z)、ny(x,y,z)、nz(x,y,z)映射到COMSOL材料的张量表达式里。实用上我建议把问题解耦先用向列相液晶弹性理论或直接拿现成工具如LCDSim、LC3D、JEDIMS等算出电压下的指向矢分布再把它作为已知分布导入COMSOL光学仿真里只算COMSOL电磁场不动耦合方程。这样做的好处是计算量可控、各模块独立也更方便做参数化研究——你想扫电压、扫角度、扫波长哪个都行。真要严格做液晶-电磁双向耦合COMSOL也能做但那基本是给自己找罪受工程和论文里绝大多数都走解耦路线。2. 张量矩阵设置把分子指向变成能算的介电性质2.1 为什么液晶不能只填两个折射率液晶是典型的单轴各向异性材料。以最常见的向列相液晶5CB为例它可以视作棒状分子沿一个平均方向取向这个方向叫指向矢。沿分子长轴方向光感受到的折射率是异常光折射率n_e垂直长轴方向光感受到的是寻常光折射率n_o。对一个特定的传播方向、特定的偏振你确实只要取一个等效折射率这就是做实验时常用的简化。但仿真不一样COMSOL在求解麦克斯韦方程组时每个网格点都需要完整的介电张量因为它要处理两个偏振间的耦合、斜入射、以及结构几何对偏振的混合你只给两个标量折射率物理场根本不知道这个材料在不同方向上的响应差异自然算不出真正的模式分布。用一个生活化的类比液晶分子像一排排朝向固定方向的小木棍。小木棍顺着排列的方向看过去有点“粗壮”横着看就是“纤细”电磁波也是类似——顺着光轴偏振的光和垂直光轴偏振的光看到的介电常数不同。如果只在材料里写“折射率1.5”相当于把所有小木棍全叠一起当根粗柱子当然不行。2.2 从指向矢到介电张量的数学翻译单轴液晶的介电张量可以写成非常简洁的形式ε ε_o · I (ε_e − ε_o) · (n ⊗ n)其中ε_o n_o²ε_e n_e²I是单位张量n是指向矢的单位向量n⊗n是它的外积直积。这个公式的意思很直观材料的介电响应等于一个各向同性的基底ε_o I叠加一个沿指向矢方向的“额外增量”(ε_e−ε_o)n⊗n表示沿分子长轴比垂直于长轴多出的那一部分响应。如果把n写成三个分量(nx, ny, nz)张量各项就是εxx ε_o (ε_e − ε_o)·nx²εyy ε_o (ε_e − ε_o)·ny²εzz ε_o (ε_e − ε_o)·nz²εxy ε_yx (ε_e − ε_o)·nx·nyεxz ε_zx (ε_e − ε_o)·nx·nzεyz ε_zy (ε_e − ε_o)·ny·nz重点提醒非对角线元素不是随便填的。介电张量在能量上必须是实对称的无损耗时或者Hermitian有损耗时。很多新手在COMSOL里直接把εxy填一个数、εyx填另一个数结果模式完全乱掉检查半天才发现是对称性破坏。所以上面这6个式子就是你在COMSOL里要用到的核心表达式后面对照填就行。实际建模中更常用的习惯是用球坐标描述指向矢。设θ为指向矢与z轴液晶盒法向的夹角φ为指向矢在xy平面上的方位角nx sinθ·cosφny sinθ·sinφnz cosθ这样只要用θ、φ两个角度函数就能控制全空间液晶方向的任意分布非常方便。比如初始平行取向且分子沿x方向那就是θ90°、φ0加电压后分子往z轴方向倒θ逐渐变小如果做了扭曲取向φ沿z方向线性变化。2.3 COMSOL里两种落地配置COMSOL的“电磁波频域”物理场下材料属性可以选“相对介电常数”或“折射率各向异性”。两种方式都能表达液晶但适合的场景不同我按实操经验说方式A用“相对介电常数”完整矩阵按变量表达式填写。在材料属性里选择“相对介电常数”取消“由折射率计算”把6个独立分量手工填成上面的公式。这种方式最灵活因为nx、ny、nz可以是任意空间变量、插值函数、参数。这也是实现“任意液晶分布”的最推荐方式。方式B用“折射率各向异性”加材料坐标系旋转。在材料里填nx、ny、nz三个折射率值然后定义材料坐标系的旋转欧拉角。这种方式适合液晶方向在全局范围内统一比如整个液晶盒都是沿某个固定方向但不能处理空间渐变分布因为材料坐标系的旋转角一般是常数没法让它随xyz变动。另外提一句金属超表面金、银不要用“无损折射率”记得用复折射率数据。COMSOL内置材料库里有部分金属数据但很多情况下更可靠的是自己查Johnson和Christy的实测数据表然后插值进材料里。如果波长在你的研究范围内有已知的光学常数直接填“折射率nik”是省事的做法。最后给个填写示例截图参考的描述。假设你定义了两个全局参数th_lcθ单位deg和phi_lcφ单位deg且液晶n_o1.50、n_e1.70那么先定义变量nx sin(th_lc*pi/180)cos(phi_lcpi/180)ny sin(th_lc*pi/180)sin(phi_lcpi/180)nz cos(th_lc*pi/180)然后在介电常数矩阵里填εxxn_o^2 (n_e^2 - n_o^2)*nx^2εyyn_o^2 (n_e^2 - n_o^2)*ny^2εzzn_o^2 (n_e^2 - n_o^2)*nz^2εxy与εyx(n_e^2 - n_o^2)nxnyεxz与εzx(n_e^2 - n_o^2)nxnzεyz与εzy(n_e^2 - n_o^2)nynz只要保证对称性后面随便怎么扫θ和φ都不会有材料设置问题。3. 任意液晶分布让每一根分子都有自己的方向3.1 解析式定义均匀取向、扭曲分布、弯曲分布做参数扫描研究时最简单的“任意分布”其实是解析表达式。比如平行排列的光学液晶盒在未加电压时分子躺在基板平面内取θ90°、φ0。加上电场后如果忽略液晶弹性形变带来的梯度等效介质近似可以认为整个液晶盒内部的θ都变成同一个值只是这个值随电压变化。那模型就可以简化成固定φ0扫描th_lc从90°到10°。每扫一个θ就相当于给液晶加了一个特定的有效电压。这样得到的相位随θ变化曲线就是你的“可调相位范围”。更复杂一点扭曲向列相TN液晶盒的φ随z呈线性变化φ(z) φ0 (π/2)·(z/d)。此时你只需在变量里定义phi_lc phi_0 pi/2*(z/d_lc)这里z来自COMSOL的变量d_lc是液晶盒厚度参数。配合前面的外积公式COMSOL在每一层z切片上都会自动换一个方位角不需要额外建模。这个技巧可以用来做偏振旋转器、液晶光栅等结构。还有π盒OCB模式的弯曲分布——分子在盒中心与两侧指向不同θ(z)一般是关于中心对称的弧线分布。工程上常近似写为θ(z)θ_m·cos(πz/d)实际偏弹性形变但作为仿真初值足够也可以直接用解析式驱动。我希望传达的核心思维是解析表达式适用于规律性较强的分布好处是轻量、参数可调、不占内存。3.2 导入真实的液晶指向矢分布数据如果你要做“任意液晶分布”光靠解析式是不够的。真实的液晶盒在电压下指向矢分布要考虑锚定能、弹性常数K11、K22、K33和介电各向异性Δε分布往往不是教科书正弦形状。这时你需要外接工具算出分布再导进COMSOL。具体做法分三步在LCD仿真工具或者自己写的小程序里把液晶盒划分成网格输出每个点的指向矢分量nx、ny、nz或者是θ、φ。保存成文本文件格式通常是四列x坐标、y坐标、z坐标、数值。注意COMSOL的“插值函数”导入的是“x,y,z,value”形式。在COMSOL“全局定义”里添加“插值函数”设定为3D插值x、y、z三个自变量返回值是指向矢的一个分量。建立三个插值函数nx_interp(x,y,z)、ny_interp(x,y,z)、nz_interp(x,y,z)。然后定义变量nx nx_interp(x,y,z)ny ny_interp(x,y,z)nz nz_interp(x,y,z)。用上面的外积公式填介电张量此时自然就是“每个体素都有自己的张量”。这里有个非常影响成败的细节插值函数默认是“线性插值”还好但超出数据范围的部分要么外推、要么钳制限制在范围内。如果指向矢数据只覆盖液晶盒区域而几何模型里这个区域之外还有一点点网格就可能导致插值函数报错或给出荒谬值。建议在插值函数属性里把“边界外”设为NaN处理或者把数据范围定义得比几何边界略大一圈确保每个网格点都有有效值。另一个习惯是导入数据前先把指向矢归一化。因为后面张量公式假设n是单位向量数值计算软件里你从外部导入的nx、ny、nz必须满足nx²ny²nz²≈1误差越小越好。如果外部数据没有归一可以在COMSOL里定义归一变式nx_norm nx/sqrt(nx^2ny^2nz^2)千万别图省事直接用原始数据填张量否则张量的本征值会偏谐振位置直接错掉。3.3 参数化扫描代替电压耦合工程简化的价值很多初学者上来就问COMSOL能不能同时算液晶弹性方程和电磁波方程能但你得有强大的理由才值得——比如你研究的就是液晶形变与超表面近场之间的动态相互作用且时间尺度一致。大多数情况从论文复现的角度做一个解耦的参数扫描就是最有效率的方案把th_lc或某个电压映射值作为扫描参数跑一遍频域仿真直接看S21相位怎么随θ变化。比如你想知道“这液晶加载超表面在1550nm处能不能调出180°相位差”那就把th_lc做成参数扫描77°、65°、56°、48°、40°……每个θ跑一个频点最后画出相位-θ曲线一目了然。如果显示范围不够就要去优化超表面的谐振强度或液晶层的厚度。如果你真的需要“电压-光学相位”的完整映射我建议分两段走第一段用弹性理论或LCD软件算出电压V对应的θ分布第二段把θ分布当参数导入COMSOL。这样做出来的V→θ→相位曲线照样能发表且每一步都可检验、可替换。这是很多成熟课题组的做法。4. 从零搭一个液晶-超表面单元完整实操记录4.1 几何参数与材料数据准备下面是给新手一个能直接开跑的基准设计参数参考的是近红外波段可调超表面单元的常见配置。注意这套参数不是最优的只是让你第一次跑通流程用的“保底参数”后面要按谐振位置和相移量做优化。单元结构自下而上层材料厚度/尺寸备注基底熔融石英500 nm折射率约1.4441550nm超表面谐振器金椭圆纳米柱长轴300 nm短轴140 nm高220 nm复折射率数据可用Johnson-Christy取向层可选如PVA5 nm初期可不建影响极小液晶盒向列相液晶2 μmno1.50ne1.70近似取值上基板熔融石英500 nm可画可不画看你要不要计算反射单元周期p800nm。这个周期在1550nm附近小于波长不会出现衍射级适合做相位调控。上面的金椭圆纳米柱是典型的传播相位单元谐振方式近似于局域表面等离激元共振LSPR。液晶取向初始设为平行排列分子指向沿x方向即θ90°、φ0。加电场时分子抬高θ逐渐减小。geometrically液晶层直接长在超表面上覆盖整个单元用“长方体”画就行xy范围与单元周期一致z从金柱顶部延伸到液晶盒顶部。金椭圆柱的绘制建议直接在x-y平面画一个椭圆长短轴沿x/y轴然后拉伸extrude220nm。如果你想让椭圆长轴方向可变可以用参数定义长轴与x轴的夹角。这个夹角本身也能作为扫描参数——超表面几何各向异性配合液晶取向是后面做偏振相关调控的基础。4.2 物理场、边界条件与端口配置物理场选“电磁波频域”。这个模块注意在COMSOL的版本里叫“Electromagnetic Waves, Frequency Domain”属于波动光学模块。如果你只有RF模块也能做但光学材料、端口设置的表现形式不如波动光学方便。边界条件要点如下周期方向在单元的两个方向x和y都加“周期性条件”类型选“Floquet周期”如果正入射倒格矢的相位偏移kx0、ky0。正入射时其实用最简单的“周期性边界”就行但Floquet更通用推荐直接养成用Floquet的习惯后面换斜入射不用重新建模。传播方向z向上下两端加“端口”边界。端口1在上方入射端口2在下方透射。每个端口要指定一个模式这里直接指定TE模或TM模。COMSOL会内部求解端口的模式分布但对于我们这种平面波入射电场基本是均匀的模式求解非常快。你需要格外注意端口距离超表面至少半波长以上避免高阶倏逝波污染S参数。入射场设置端口1的“激励”设为“开”端口2设为“关”。这样S21就是从端口1到端口2的透射系数S11是反射系数。相位信息就在arg(S21)里。一个常犯的错误端口模式偏振定义与预期相反。COMSOL端口默认的模式极化由波导截面几何决定而我们希望端口1的模式是x偏振或y偏振的平面波。建议在“端口”设置里的“模式”栏中手动输入电场分量例如设Ex1、Ey0并对应波长传播方向。设置完后先别急着跑扫描先跑一个频点查看端口的电场分布是否如预期均匀分布这一步能省后面几个小时。4.3 研究配置与求解器选择启动参数频率范围取中心波长1550nm对应的频率附近193.5 THz扫±30 THz。如果你想看清谐振曲线建议扫的范围足够包含谐振谷频点密度在谐振附近至少间隔0.1 THz。计算S参数在“研究”里建议加一个“S参数”的步骤或者在后处理里直接提取端口变量。COMSOL提供了“S-parameter”表达式可直接全局计算。求解器对于这种单元胞模型网格自由度通常在几十万到一两百万之间。默认的“物理场控制网格”往往过于稀疏要手动加密。推荐用“用户控制网格”液晶层至少剖两层六面体/扫掠网格超表面椭圆表面加三角形边界层网格。求解器选“直接求解器MUMPS或PARDISO”鲁棒性最好如果你的内存吃紧再考虑迭代求解器的GMRESMultigrid但收敛性需要调。关键技巧频率扫描或参数扫描一定开启“辅助扫描”Auxiliary sweep。在“研究”设置里把“频率”作为辅助扫描变量勾选“为每个扫描步骤计算初始值”。这样做的好处是每扫下一个频点/参数COMSOL会把上一步的解作为初值能显著加速而且避免每个点都从零开始求。我在默认配置下发现不开辅助扫描时一个3D液晶超表面参数扫描要跑两三个小时开了之后经常可以压缩到几十分钟。后处理取相位的路径全局计算→选择“S-parameter, port1, port2”表达式→输出“arg(S21)”。注意这是未解包裹的相位跨过±π边界时会跳变。若要连续曲线要么在COMSOL里启用相位解包裹功能部分版本支持要么导出数据后在MATLAB/Python里np.unwrap处理。这是老玩家都知道的坑后面专门说。5. 常见错误与排查实录5.1 症状对照速查表做这类液晶超表面模型我踩过的坑比踩过的坑还多这里直接给一份速查表症状可能原因处理方式网格从液晶层穿过去就报错各向异性层网格太粗张量随坐标变化剧烈液晶层至少剖2~3层如果指向矢变化较剧烈额外加密集网格S21相位出现明显的锯齿跳变相位未解包裹arg()在±π处跳变后处理里unwrap或用unwrap(S21相位)表达式谐振峰位置跟文献差几十纳米金属复折射率用了常温近似液晶双折射取值不准换用精确色散数据或Drude模型核对n_o/n_e取值偏振结果和理论矛盾如TE/TM调换端口偏振定义反了查看端口电场分布按预期修改模式分量加电压扫描θ后相位不变液晶双折射与入射偏振方向成90°或超表面对该偏振不敏感检查液晶光轴与入射偏振是否匹配调整φ或改用椭圆柱旋转调控内存直接爆掉3D全波网格太多求解器内存消耗大用SPM稀疏矩阵存储或减小网格开启辅助扫描减少重复求解插值函数导入后材料张量出现NaN插值数据越界设置插值范围外为约束/NaN处理或扩大插值范围反射曲线很大且不随θ变化液晶层厚度过大形成法布里-珀罗干涉叠加减小液晶层厚度或在结果中扣除干涉背景或厚度的1/2波片条件附近操作5.2 相位跳跃单独拿出来说这几乎是人人必踩的一坑。S21相位是从复数比值arg()算出来的数学上phases are defined mod 2π。你扫频率或扫θ时S21会绕着原点转圈arg()就不断地在π和-π之间往返跳。跳变本身物理上没有意义但画出来的相位曲线看起来很“毛刺”让人误以为相位不连续。解决办法第一招就是unwrap。在COMSOL里如果你用的是“全局计算”表格直接把结果导出到Excel里加一个解包裹逻辑或者在MATLAB用unwrap函数。如果论文需要直接在COMSOL里出图可以在“二维截线绘图”或“全局”绘图的表达式里写“unwrap(phase(S21))”部分新版本支持这种写法。谨慎起见我一般都在后处理里用Python把数据拉出来做unwrap顺便还能做平滑和归一化。模型的物理相位变化是连续、单调的只要你看到一条平滑的、覆盖0~2π的相位曲线就可以放心它的调相能力。另外扫描步长要足够密。如果θ扫描步长很大导致相位在相邻两个点之间跳了超过π那unwrap也救不回来——因为算法不知道是π还是-π方向。一般扫描θ从70°到10°分成30步以上每个频点的相位连续性才有保障。5.3 资源节省和模型扩展的进阶心得如果你发现3D全波仿真太慢有两个成熟技巧。技巧一先用2D模型验证。把椭圆柱换成无限长矩形柱模型从3D缩成2D x-z平面算出的谐振类型不一样但可以快速验证液晶张量设置对不对、端口有没有问题。确认2D跑通后再上3D调试效率高得多。技巧二谐振峰附近才细扫。先跑一个粗扫描5 THz步长把谐振峰大致位置找出来再在峰附近加密扫描省时间也便于看清相位在谐振附近剧烈变化的行为。更进一步如果你的目标是“宏观效果演示”比如偏转角度、聚焦效率没必要对所有单元逐一建模。单胞仿真拿到S参数后可以用超表面等效表面阻抗/偏振张量去搭宏模型这就是典型的“多尺度仿真”思路。液晶层对单元的影响已经包含在你算出来的相位曲线里了宏模型直接查表调用即可。我还想强调一点目前这类设计里液晶盒都较薄1~3μm所以光是直接穿过液晶层经历的几何相位延迟有限真正大范围调相主要靠你用液晶修调超表面谐振的方式实现。这意味着你在分析时要把“液晶单独导致的相位延迟”和“超表面谐振对液晶扰动的放大”分开看。我自己的习惯是同一个模型里跑三组对照——只有超表面、只有液晶、液晶超表面。三组曲线的差异能清晰揭示谐振耦合到底贡献了多少调相范围这也是审稿人最爱问的问题提前把这个分析做好后面写论文会舒服很多。5.4 想延伸做动态切换几个后续方向仿真模型跑通不是终点后续可以做很多扩展把θ扫描换成电压参数。这时你只需要外部定义V→θ的映射函数即可最简单的形式是θ(V)90°-arctan(V/Vth类似的经验关系更精确的用液晶弹性方程的数值解。把单胞复制成阵列。超表面人工设计的相位梯度如线性分布可以让透射波发生偏转这是可编程超表面最基本的应用。你需要在单胞仿真里获得相位曲线然后对不同位置赋予不同θ相当于不同单元的液晶控制电压不同再组装成大阵列。这个过程不建议直接在COMSOL里建几千个单元的大模型而是每个单元单独仿真最后用等效全波阵列理论合并。把频率从近红外换到可见光波段。可见光波段金/银谐振结构更小周期也更紧凑300~500nm液晶层厚度可能要降到几百纳米此时网格和求解器设置基本复用但材料色散数据要重新核对要特别注意液晶、基底在可见光的吸收损耗不可忽略。如果你对液晶的动态响应过程本身感兴趣——比如分子转向期间透射率的瞬态变化——那就需要引入时间维度COMSOL的“移动网格”配合“瞬态求解器”可以做液晶动力学与电磁场的耦合问题类似相场/动网格的仿真思路。不过这是另一个复杂的课题不适合与本文的单元建模混在一起做。先把静态全波模型跑通再逐步加复杂度是比较稳妥的路线。最后再分享一点个人体会液晶超表面仿真最难受的阶段不是求解器而是材料张量写得对不对、端口方向对不对、网格路径对不对。这些错误都很隐蔽跑出来结果也能画但就是物理上荒唐。我通常的做法是先用极端参数验证——比如把Δn设为0液晶变成各向同性看相位曲线是否跟单纯超表面一致或者把θ定为90°且φ0检查张量是否退化成diag(n_e²,n_o²,n_o²)再用电场分布图确认传播方向上确实存在预期的双折射相位延迟。只要这两步验证过了后面的参数扫描和优化才算真正可信。希望这篇实战记录能帮你绕开我当年走的弯路一次把模型跑通。