1、2023年9期研究视界科技创新与应用Technology Innovation and Application二维蜂窝状光子晶体能带的数学模型与仿真研究杨之（北京邮电大学，北京 100876）光子晶体是由宏观光学介质代替微观粒子按一定规律周期性排列的人工微观结构，其周期阵列限制了电磁波在光子晶体中的传播，对特定波长的光子具有阻挡作用，从而形成了光子带隙1。根据应用需要设计光子晶体的结构，由于光子带隙的存在可以禁止光在一定频率范围内的传播，使得光子晶体有着广泛的应用，其应用包括窄带滤波器、分复用器等2-5。获得宽的光子晶体带隙是设计高性能光子晶体器件的基础6，所以如何设计光子晶体结构来获得更宽的
2、带隙一直是光子晶体领域的研究重点。已有研究结果表明，二维蜂窝状光子晶体具有较大的带隙7，且带隙大小和介质柱的半径大小与介电常数密切相关8-9。改变光子晶体的晶格结构、降低结构对称性和引入新材料等措施都可以增大光子晶体的带隙，但是由于直接进行实验来设计新型光子晶体结构的成本巨大，现在一般先采用数值模拟方法进行结构设计，以期获得结构设计的指导。目前，数值模拟二维光子晶体带隙的常用方法有平面波展开法10-11、时域有限差分法2，12和有限元方法13，其中有限元方法因其计算效率较高，常用于三角晶格光子晶体和正方晶格光子晶体的数值模拟，但目前在蜂窝状光子晶体领域研究较少。本文采用有限元法研究二维蜂窝状光
3、子晶体的能带结构，以往在对蜂窝晶格光子晶体进行建模时通常选择正六边形结构作为单元晶格，本文在数值模拟过程中直接采用了易于操作的三角晶格胞元，使建模仿真和边界处理更加容易。并利用 COMSOL 软件计算了二维蜂窝状光子晶体的能带结构，通过对比文献14的结果，验证了方法的正确性和有效性。进一步，本文采基金项目：国家自然科学基金（11671052）作者简介：杨之（1997-），女，硕士研究生。研究方向为计算数学。摘要：该文基于有限元方法研究新型的二维蜂窝状光子晶体的能带结构，并运用数值模拟方法进行不同材料和结构参数的蜂窝状光子晶体的设计。对于二维蜂窝状光子晶体，在三角晶格胞元上采用线性有限元离散结合
4、 COMSOL 软件可以实现晶格的能带结构求解，数值结果表明数学模型和仿真算法的正确及有效性。在带隙最大化的优化问题研究中，基于数值模拟方法和线搜索，探究二维蜂窝状光子晶体的参数与带隙的关系，并应用新型材料构造宽带隙的新型二维蜂窝状光子晶体结构，为新材料的应用提供数值模拟方法和设计指导。关键词：蜂窝状光子晶体；有限元方法；能带结构；数值模拟；优化设计中图分类号：O242.21文献标志码：A文章编号：2095-2945（2023）09-0098-06Abstract:In this paper,the energy band structure of a new two-dimensional
5、cellular photonic crystal is studied based on finiteelement method,and the cellular photonic crystal with different materials and structural parameters is designed by numericalsimulation method.For two-dimensional cellular photonic crystals,the solution of lattice band structure can be realized by u
6、singlinear finite element discrete and COMSOL software on triangular lattice cells.The numerical results show that the mathematicalmodel and simulation algorithm are correct and effective.In the study of the optimization problem of band gap maximization,basedon numerical simulation method and line s
7、earch,we explored the relationship between the parameters of two-dimensional honeycombphotonic crystal and the band gap,and applied new materials to construct a new two-dimensional honeycomb photonic crystalstructure with wide band gap,providing numerical simulation methods and design guidance for t
8、he application of new materials.Keywords:cellular photonic crystal;finite element method;band structure;numerical simulation;optimal designDOI：10.19981/j.CN23-1581/G3.2023.09.02498-研究视界科技创新与应用Technology Innovation and Application2023年9期用数值模拟的方法探究了二维蜂窝状光子晶体的参数与带隙的关系，并通过在蜂窝状光子晶体中应用新型材料构造了宽带隙的新型二维蜂窝状光子
9、晶体结构。1能带结构计算的数学模型在本节中，将介绍二维蜂窝状光子晶体能带求解的数学模型。考虑由周期介电圆柱组成的二维蜂窝状光子晶体阵列，如图 1（a）所示，a0为晶格常数，即相邻2 个正六边形中心的距离，a1和 a2是长度 a0的向量。放大的单个六角蜂窝晶格如图 1（b）所示，其中 R 为六边形的边长，d 为圆柱体的直径，a和 b分别是周围背景和介电圆柱的介电常数。本文采用三角晶格胞元进行数值计算，光子晶体等效至 xoy 平面后选取的三角形晶格胞元，如图 1（c）所示，其边长为 a0。注：（a）蜂窝状光子晶体；（b）蜂窝晶格局部放大图；（c）计算采用的三角单元晶格。图1蜂窝状光子晶体示意图由图
10、 1（a）可知，三角形晶格的晶格基矢为 a1=a0（1，0），a2=a0（1/2，3/2），其对应的倒格基矢为 b1=k0（3/2，-1/2），b2=k0（0，1），其中 k0=4/3a0。靠近原点有 6 个倒格点，其垂直平分线所围成的区域就是第一布里渊区，由图 2 可知 KM 所围成的三角形即为第一布里渊区，3 个顶点的坐标分别为 K=k0（3/3，0），=（0，0），M=k0（3/4，1/4）。图2三角晶格第一布里渊区下面根据图 1（c）所示的三角单元晶格，推导光子晶体特征值问题。在经典电磁理论中，电磁现象可由 Maxwell 方程组（1）描述。在光子晶体中，研究的对象是电磁波，其在光子晶
11、体中传播时满足 Maxwell 方程组，（1）式中：E 是电场强度；H 是磁场强度。对于二维光子晶体，直接对 Maxwell 方程进行求解比较困难，由于电磁波可分解为简谐波，即有 E（r，t）=E（r）e-it，H（r，t）=H（r）e-it，这里 是频率。代入式（1）可得到与时间变量无关的 Helmholtz 方程，即电场方程（2）和磁场方程（3）考虑二维光子晶体带隙问题时，可简化为 TM 和TE 两种偏振情形。在 TM 模式下，有 E（r）=（0，0，ER（r?）；在 TE 模式下，有 H（r）=（0，0，HR（r?）。这里r?是二维坐标（x，y）。代入式（2）和（3），分别可以得到 TM
12、 模的电场波动方程（4）和 TE 模的磁场波动方程（5）式（4）和（5）中的=（x，y），不同于前面的二维算子。由于二维光子晶体的周期性，满足 Bloch 定理，可将上述波动方程化至单位原胞上再进行数值实验。在二维情况下，记=c()2，k=（a，b）表示的是二维平面空间内的波矢矢量，其取值范围为不可约布里渊区 B。首先考虑 TM 模。（6）是求解的二维集合区域，在单位原胞边界=上，满足边界条件，（7）式中：1和 3是单位原胞的一组对边；2和 4是单位原胞的另一组对边。求解区域如图 1（c）三角单元晶格所示。0DHJtBEtBD|=+|=-|=|=|?()()()01()()0()E rr E
13、rcH rH rrc|-=|-=|，（2）。（3）?()()()01()()0()?Err ErcHrHrrc|+=|+=|?，（4）。（5）()()()()()0RRiiErr Er+=?kk?|?EEEE=;MKa0a1a24baR321d99-2023年9期研究视界科技创新与应用Technology Innovation and Application求解 TE 模的过程类似，为方便讨论变分过程，可将 2 种模式统一记成如下特征值形式，给定矢量 k，求 和 uR使得，（8）式中：=1，TM1（r?），T|E；（r?），TM1，TE；是待求特征值；uR表示待求场量即特征向量。只要矢量 k
14、遍历第一不可约布里渊区域的边界，就能够绘制出能带结构图。下面针对式（8），给出其有限元变分过程和离散过程。首先，定义 Sobolev 空间要求解的特征值问题为：对kB，求满足（k，u）CHk1（）对于固定的波矢 k，当求解区域是线性三角元时，u、v在单元 e 内的近似 ue、ve，取线性三角元的基函数 Nei、Nej和 Nek的线性组合，这里所说的基函数是关于（x，y）的势函数，不妨取试探函数 ve为单元基函数，则有这里 x 取值i，j，k，而 ui，uj，uk是待定常数。由此得到单元刚度矩阵，然后根据线性三角元局部编号和整体编号关系，生成总刚度矩阵 A 和 B，得到变分问题对应的离散特征值问
15、题：AU=BU。（13）在求解广义特征值问题时，还需要考虑周期边界条件（7），从而得到最终的待求解广义特征值问题AU?=BU?.（14）对于形如（14）的广义特征值问题，在二维蜂窝状光子晶体带隙模型中，矩阵 A 和 B 常具有正定性，则有 Cholesky 三角分解B=LLT，（15）则对式（14）两侧的左侧乘 L-1，并在右侧整理提出因子LT得（L-1AL-T-I）LTU?=0。（16）式（16）可以看成以 为特征值，LTU?为特征向量的特征值问题，并且式（16）和（14）具有相同的特征值，系数矩阵仍满足正定性，这样就将广义特征值问题（14）转化为特征值问题（16），可利用直接法或迭代法进行
16、求解特征值问题（16）。2数值仿真结果2.1算例1：二维蜂窝状光子晶体仿真本节考虑图 1 所示的二维蜂窝状光子晶体，周围背景的介电常数是 a=1，圆柱体的介电常数是 b=11.7，六边形的长度为 R，圆柱体截面圆形的直径为 d且 d=2/3R。选择三角单元晶格模拟得到的能带图如图3 所示，在其他参数不变的情况下，仅改变晶格参数的大小，此时三角晶格中圆形的大小和位置会发生改变。对不同参数取值下的能带结构进行仿真，结果如图 3所示，当 a0/R=3.125 时，带隙归一化频率在 0.480.51之间；当 a0/R=3 时，带隙完全消失；当 a0/R=2.8 时，带隙归一化频率在 0.440.48 之间。3 种情况下的计算结果均和文献14中的仿真结果一致，并且仿真所用时间比采用常见三角晶格进行仿真时间更短、效率更高。保持其他参数不变，仅改变晶格参数 a0/R，保证蜂窝晶格存在的条件下可使得 a0/R 取值从 2.79.7，得到TM 模式下蜂窝状光子晶体的带隙变化图，图 4 所示可以清晰地看到带隙随晶格参数 a0/R 的变化情况。当填充比为3 时，完全禁带的宽度为 0；当晶格参数 a0/R=5
