太赫兹微结构功能器件最新章节_常胜江著

3.1 数值仿真方法

3.1.1 时域有限差分法

时域有限差分法（Finite Difference Time Domain Method，FDTD算法）是一种利用计算机对电磁场进行数值计算的方法 ^[1] ，可以用来解决太赫兹波在功能器件中的传播问题。其基本思想是将麦克斯韦方程在空间和时间上做离散处理，利用周期性边界条件，将空间网格点的场分布函数随时间的推移进行递推演化，计算电磁场的传播及其与物质的相互作用。在设计太赫兹功能器件时，FDTD算法能够对器件传输特性和场分布特性进行仿真，原则上可以处理超材料、光子晶体和亚波长孔阵列等任意结构的问题，只要空间网格划分足够精细，便可保证计算结果的精确度。

常用的空间网格划分方式为Yee网格，其单元结构Yee元胞如图3.1所示，节点的坐标为（ i ， j ， k ） = （ i Δ x ， j Δ y ， k Δ z ），其中Δ x 、Δ y 和Δ z 分别为节点沿 x 、 y 和 z 方向的空间步长，于是电磁场函数可表示为

图 3.1 FDTD算法中的Yee元胞

式中，Δ t 为时间步长； n 为整数。由式（3.1）可知，电磁场函数 ϕ （ x ， y ， z ， t ）对各个变量 x 、 y 、 z 和 t 的偏微分可表示为

对于Yee元胞，电场和磁场各取样节点在空间和时间上呈交替分布，每个磁场分量由四个电场分量环绕，同时每个电场分量也由四个磁场分量环绕，这样的取样方式符合法拉第电磁感应定律和安培环路定理，且适合用麦克斯韦方程进行差分计算，三维空间的麦克斯韦方程的差分形式为

空间步长和时间步长的取值越小，对计算机配置的要求越高，计算所需时间越长，计算结果越精细。另外，时间步长Δ t 的取值与空间步长Δ x 、Δ y 和Δ z 有关，需满足如下条件：

否则将影响电磁场数值计算结果的稳定性。

此外，FDTD算法受计算机中央处理器和内存的限制，只能在有限的区域进行模拟计算，这就要求在计算区域的边界上对电磁场分布进行限制，通过在有限的空间里施加边界条件可以实现无限空间的计算。我们在器件的仿真过程中采用的边界条件均为完美匹配层（Perfectly Matched Layer，PML），它将电磁场分量在计算区域边界处分解，并对每个被分解的分量进行损耗，这样就相当于在计算网格的边界处设置了一种特殊的吸收介质，其波阻抗与相邻介质的波阻抗完美匹配，电磁波将无反射地进入PML消耗殆尽，对于计算区域内的电磁波而言等同于在无限大的空间上传播。

3.1.2 有限元法

有限元法（Finite Element Method，FEM）是一种频域算法 ^[2] ，其基本思想是把连续的几何结构离散成有限个单元，并在每一个单元中设定有限个节点，从而将连续体看作仅在节点处相连接的一组单元的集合体，同时在每一个单元中假设一个近似插值函数以表示单元中场函数的分布规律，再建立用于求解节点未知量的有限元方程组，从而将一个连续域中的无限自由度问题转化为离散域中的有限自由度问题。这种区域分割的方法与FDTD算法不同，FDTD算法采用网格切分，只要求求出子区域网格结点上的场值，实际上仍采用点逼近。而FEM用简单的子单元逼近，每一个子单元上都用一个简单函数描述，求出的结果则是小单元的平均近似解。典型的基于FEM的数值模拟软件有COMSOL中的射频模块等。

与FDTD算法的矩形或长方体网格不同，FEM的网格多是三角形或三棱体，因此其在进行局部网格加密时更加方便，特别适合于处理具有复杂几何形状物体和边界的问题。此外，频域算法是对所求域内的频率点逐点计算，在求解稳态电磁分布、计算窄带电磁波传输、处理非均匀色散介质等方面更加有利。当所求电磁场包含强烈谐振而导致其时域信号非常长时，使用FEM比使用FDTD算法更加高效。

3.1.3 平面波展开法

求解周期性结构的光子带隙是研究其电磁性质的基本方法，即求解特定波矢的本征值问题，比较常用的方法有FDTD算法和平面波展开法（PlaneWaveExpansion Method，PWE算法） ^[3] 。PWE算法以布洛赫定理为基础，先对某一入射波矢 k 用平面波的方式展开，然后将麦克斯韦方程组化为一个本征方程组，进而求解 k 的一系列频率本征值。相比于FDTD算法，PWE算法更加精确和高效。

在求解光子带隙时，只需要给出周期结构的重复单元即可。为了避免计算边界处网格划分、材料突变引起的计算误差，一般将结构继续向外延伸，使其略大于计算区域。有些仿真软件（如Rsoft软件的Bandsolve模块）在计算光子带隙时会计算出结构在 k 空间内所有的色散关系曲线，这时只需要正确地判定结构的晶胞类型，合理地划定第一布里渊区即可。如果只关心沿某一传播方向的传输带隙，则需要在这个方向设置弗洛奎特周期性边界条件，其他边界则根据所求物理模型的不同而分别设置。从0到πrad对弗洛奎特周期性边界条件的相位进行扫描，分别计算各种相位下的特征频率，就可以得到结构沿该方向传播时的色散关系曲线。位于各条色散关系曲线之间，没有特征波矢与之对应的频段就是光子带隙。