本文提出了一种精确有效的方法来处理低剖面周期结构。通过拉普拉斯变换和中心微分法，将周期结构的Y矩阵参数直接嵌入到电场的更新方程中。与早期的工作相比，该方法可以对整体结构直接采用粗网格建模，时间步将不再受金属层厚度的限制；可以避免复杂的频域至时域的卷积运算，相关的频域函数可以由矢量拟合法轻松获得；计算时间快，精度高。
1 INBC-FDTD步进公式推导 图 1给出了低剖面周期结构置于两FDTD网格Γ-和Γ+之间的侧视图。图中：n-和n+分别为由Γ-指向Γ+的向量和由Γ+指向Γ-的向量。假设该结构沿着z方向上的厚度远小于所关注的最高频段下的波长和结构的横向尺寸，则结构切向电场和切向磁场之间的关系可通过引入如下Y矩阵得到：
 | (1a) |
 | (1b) |
 |
| 图 1 置于两FDTD网格Γ-和Γ+之间的低剖面周期结构侧视图 Fig. 1 Side view of low-profile periodic structure placed between two FDTD grids Γ-and Γ+ |
| 图选项 |
式中：s=jω；Y11和Y22分别为Γ-和Γ+的自导纳，Y12和Y21分别为Γ-和Γ+的互导纳。
通过将式(1a)和式(1b)写为一系列表达式，可得
 | (2a) |
 | (2b) |
 | (2c) |
 | (2d) |
在式(2a)~式(2d)中，与频率相关的Y矩阵系数Y11(s)、Y12(s)、Y21(s)和Y22(s)可以通过矢量拟合法[13]近似为由实数或共轭复数组成的有理分式之和，表示为
 | (3a) |
 | (3b) |
式中：pn和rn分别为极点和留数，其值可为实数或者共轭复数；d为实数；N为极点的个数。
为了便于简单分析，仅对式(1a)作进一步推导。通过将式(3a)和式(3b)代入式(2a)和式(2b)，得到
 | (4a) |
 | (4b) |
引入状态方程Xn, Y11(s)、Xn, Y12(s)、Xn, Y21(s)和Xn, Y22(s)：
 | (4c) |
 | (4d) |
表达式(4c)和式(4d)可重写为
 | (5a) |
 | (5b) |
 | (5c) |
 | (5d) |
上述表达式(5a)~式(5d)可以通过拉普拉斯变换
 | (6a) |
 | (6b) |
 | (6c) |
 | (6d) |
式中：Δt为采样时间。
将式(5a)和式(5b)转换至离散时域，其在时刻为k+1/2的表达式为
 | (7a) |
 | (7bs) |
通过将式(7a)~式(7b)代入至FDTD的电场Ey+k+1和Ey-k+1时域步进公式[14]，可以得到内部阻抗边界条件(INBC)的离散时域表达式为
 | (8a) |
 | (8b) |
式中：
 | (9a) |
 | (9b) |
 | (9c) |
 | (9d) |
其中：Δx、Δy和Δz分别为沿x、y和z的空间步长；i、j和k分别为FDTD电场和磁场的空间坐标；ε0为自由空间的介电常数。
总体来说，INBC-FDTD方法的场更新步骤可归纳如下：
步骤1??更新状态函数式(6a)~式(6d)。
步骤2??使用式(9c)和式(9d)更新系数F1y和F2y。
步骤3??使用式(8a)和式(8b)更新电场Ey-k+1和Ey+k+1。
Ex+k+1和Ex-k+1的时域步进公式可以按照以上方法作类似推导。对于磁场的时域步进公式，仍采用常规的FDTD公式进行计算。
2 二维周期结构计算空间 对于二维周期结构，当采用FDTD方法对其进行计算空间离散建模时，由于其在二维方向的周期性，大尺寸多周期的二维结构的模型可以等效为二维周期延拓方向均设置周期边界的一个周期单元。图 2给出了二维周期结构的计算空间示意图。可以看出，沿着二维周期延拓方向设置垂直的周期边界，各向异性完美匹配层(UPML)吸收边界平行于周期结构表面。吸收边界通过设定吸收边界上的入射电压为0，即可正确地模拟无反射截断边界条件。
 |
| 图 2 二维周期结构的计算空间示意图 Fig. 2 Schematic diagram of computational space for two-dimensional periodic structure |
| 图选项 |
因为FDTD方法是基于时间和空间的离散方法，需要选择空间观察面/点来记录入射或反射时域波形。通过给定入射平面波激励：
 | (10) |
可以观察入射、反射和透射波形，通过对时域波形进行时频转换，可以得到入射、反射及透射波形的频谱图。通过
 | (11) |
可以计算得到该结构的频域反射和透射系数。式中：fref(s)、ftran(s)和finc(s)分别为反射波、透射波和入射波的频域函数。
3 数值计算与分析 基于对INBC-FDTD方法的详细推导，本节给出2个计算实例，用于验证所提方法的正确性和高效性。采用如图 2所示的计算空间，其中吸收边界UPML与入射波边界的距离设为10个元胞的距离，另外，观察面与低剖面周期结构的距离Dz设置为30个网格。
3.1 十字贴片型FSS 考虑如图 3所示的一个经典的十字贴片的频率选择表面(FSS)作为计算模型，激励入射波为中心频率12.27 GHz的双余弦调制高斯脉冲，其激励形式表达式见式(10)。用矢量拟合法对该结构的Y矩阵参数频域函数进行有理式拟合，对应的多项式参数在表 1和表 2中给出。采用一组网格参数Δx=Δy=Δz=0.25 mm(Δt=Δx (y, z)/(2c0))用于对FSS进行建模，并将其应用于INBC-FDTD方法中对结构的电磁特性参数进行数值计算。图 4(a)、(b)分别给出了入射波和透射波在时域的波形图；图 5给出了频域散射参数的幅度值。为了便于比较，将相同的网格参数应用于亚网格FDTD方法来对结构进行数值模拟。另外，选取一组尺寸为Δx=Δy=0.25 mm，Δz=10 μm的网格应用于传统FDTD方法，来对结构进行数值计算。这3种方法的具体计算参数在表 3中列出。可以看出，本文INBC-FDTD方法相对于其他2种方法来说具有较快的计算速度。
 |
| 图 3 十字贴片型FSS Fig. 3 Cross patched FSS |
| 图选项 |
表 1 提取得到的十字贴片型FSS结构Y矩阵参数频域函数的等效有理式参数：极点pn及留数rn Table 1 Extracted values of poles (pn) and residues (rn) for Y-matrix coefficients of cross patched FSS
| n | pn | rn(Y11) | rn(Y12/Y21) | rn(Y22) |
| 1 | -7.637 7×106 | 1.582 5×109 | -1.581 5×109 | 1.581 5×109 |
| 2 | -3.572 0×1012 | -1.812 0×1011 | -2.284 1×1010 | -4.030 7×1010 |
| 3, 4 | -5.499 6×107±5.627 9×1010i | 4.988 4×107±4.016 6×104i | 1.510 7×105±8.706 3×103i | -3.096 3×102±1.185 1×103i |
| 5, 6 | -5.759 2×108±5.660 8×1010i | 1.463 5×107±5.061 7×104i | 4.470 2×104±2.012 8×102i | 1.468 2×102±1.464 3×103i |
| 7, 8 | -7.357 7×107±9.089 0×1010i | 2.297 8×108±1.053 0×106i | 6.637 8×105±3.934 0×104i | 1.191 6×103±9.24×102i |
表选项
表 2 提取得到的十字贴片型FSS结构Y矩阵参数频域函数的等效有理式参数：d Table 2 Extracted values of d for Y-matrix coefficients of cross patched FSS
| Y | d |
| Y11 | 0.050 7 |
| Y12 | 0.006 4 |
| Y21 | 0.006 4 |
| Y22 | 0.011 3 |
表选项
 |
| 图 4 十字贴片型FSS的反射波和透射波时域波形 Fig. 4 Reflected wave and transmitted wave of cross patched FSS in time domain |
| 图选项 |
 |
| 图 5 不同方法得出的频域散射参数幅度值 Fig. 5 Amplitude of scattering coefficients in frequency domain obtain by different methods |
| 图选项 |
表 3 不同方法的计算参数 Table 3 Computational parameters for different methods
| 方法 | Δt | 单元个数(长×宽×高) | 时间步 | CPU时间 |
| 传统FDTD | 16.7 fs | 41×41×1 040 | 60 000 | 18 d |
| 亚网格FDTD | 0.42 ps(粗网格) 16.7 fs(细网格) | 41×41×80 | 3 000 | 14 h |
| INBC-FDTD | 0.42 ps | 41×41×80 | 3 000 | 17 min |
表选项
3.2 方环形结构的互补型FSS 图 6为单元上下表面采用方环形结构的互补型FSS[15]。其上表面为方环形贴片，下表面为在接地面上开的尺寸相同的方环形缝隙，2个表面之间为厚度h=0.5 mm，相对介电常数εr=2.65的介质板。其结构参数为：单元周期Dx=Dy=D=8 mm，方形结构长度l1=l2=7 mm，宽度w1=w2=0.8 mm，方环之间的间距g1=g2=0.2 mm。
 |
| 图 6 方环形结构的互补型FSS及其单元结构示意图 Fig. 6 Schematic diagram of complementary FSS with square ring and its unit cell |
| 图选项 |
对于如图 2所示的计算模型，将激励入射波设置为中心频率为8 GHz的双余弦调制高斯脉冲。用矢量拟合法对该结构的Y矩阵参数频域函数进行有理式拟合，对应的多项式参数在表 4和表 5中给出。使用Δx=Δy=Δz=0.25 mm(Δt=Δx(y, z)/(2c0))的网格用于对该FSS进行建模，并将其应用于INBC-FDTD方法中来对结构的电磁特性参数进行数值计算。图 7给出了透射波在时域的波形图；图 8给出了经INBC-FDTD方法计算得出的频域透射参数的幅度值。同样，为了便于比较，将相同的网格参数应用于亚网格FDTD方法来对结构进行数值模拟。这2种方法的具体计算参数在表 6中列出。可以看出，本文INBC-FDTD方法相对于亚网格FDTD方法来说具有比较快的计算速度。另外，在图 8中也给出了以上结构的透射系数(S21)测试数据。可以清楚地看到，各曲线吻合良好。
表 4 提取得到的方环形结构互补型FSS结构Y矩阵参数频域函数的等效有理式参数：极点pn及留数rn Table 4 Extracted values of poles (pn) and residues (rn) for Y-matrix coefficients of complementary FSS with square ring
| n | pn | rn(Y11) | rn(Y12/Y21) | rn(Y22) |
| 1 | -5.589 0×106 | 2.253 8×109 | -1.846 6×109 | 1.757 6×109 |
| 2 | -1.665 6×1012 | -1.036 4×1011 | -2.825 0×1010 | -2.410 3×1011 |
| 3, 4 | -1.991 8×108±3.937 3×1010i | 4.629 1×106±1.976 9×105i | -1.193 3×107±3.071 4×105i | 3.011 3×107±3.314 0×105i |
| 5, 6 | -2.030 5×108±3.973 3×1010i | 2.964 7×107±1.607 1×105i | -7.124 1×107±3.997 4×105i | 1.755 8×108±3.953 6×105i |
| 7, 8 | -3.195 7×108±6.586 5×1010i | 4.565 9×108±3.611 6×106i | -2.665 9×108±2.023 7×106i | 1.555 8×108±1.25×106i |
表选项
表 5 提取得到的方环形结构互补型FSS结构Y矩阵参数频域函数的等效有理式参数：d Table 5 Extracted values of d for the Y-matrix coefficients of complementary FSS with square ring
| Y | d |
| Y11 | 0.062 4 |
| Y12 | 0.016 9 |
| Y21 | 0.016 9 |
| Y22 | 0.144 6 |
表选项
 |
| 图 7 方环形结构互补型FSS的透射波时域波形 Fig. 7 Transmitted wave of complementary FSS with square ring in time domain |
| 图选项 |
 |
| 图 8 透射系数的频域波形比较 Fig. 8 Comparison of transmission coefficient waves in frequency domain |
| 图选项 |
表 6 两种方法计算参数的比较 Table 6 Comparison of computational parameters for two methods
| 方法 | Δt | 单元个数(长×宽×高) | 时间步 | CPU时间 |
| 亚网格 FDTD | 0.42 ps(粗网格) 16.7 fs(细网格) | 41×41×80 | 3 000 | 14 h |
| INBC-FDTD | 0.42 ps | 41×41×80 | 3 000 | 17 min |
表选项
4 结论 本文提出了一种高效求解含有金属薄涂层周期结构时域电磁特性参数的INBC-FDTD方法。考虑到金属涂覆层的厚度很薄，如果使用尺寸统一的传统FDTD元胞对结构进行总体建模，对元胞的尺寸要求往往很苛刻，需要足够小来捕捉薄层中场的变化。
本文方法有以下优点：
1) 对整体结构采用粗网格建模，从而时间步将不再受金属层的厚度的限制。
2) 通过矢量拟合法拟合金属层的网络频域函数，从而可以避免复杂的频域至时域的卷积运算。
3) 精度高，计算速度快，可有效帮助工程师进行结构设计和优化分析。
参考文献
| [1] | OKONIEWSKI M. Three-dimensional subgridding algorithm for FDTD[J]. IEEE Transactions on Antennas and Propagation, 1997, 45(3): 422-427. DOI:10.1109/8.558657 |
| [2] | WHITE M J, ISKONDER M F, HUANG Z. Development of multigrid FDTD code for three dimensional applications[J]. IEEE Transactions on Antennas and Propagation, 1997, 45(10): 1512-1517. DOI:10.1109/8.633859 |
| [3] | JURGENS T, TAFLOVE A, UMASHANKAR K, et al. Finite-difference time-domain modeling of curved surfaces[J]. IEEE Transactions on Antennas and Propagation, 1992, 40(4): 357-366. DOI:10.1109/8.138836 |
| [4] | DEY S, MITTRA R. A locally conformal finite difference time-domain (FDTD) algorithm for modeling three-dimensional perfectly conducting objects[J]. IEEE Microwave and Optical Technology Letters, 1997, 7(9): 273-275. |
| [5] | JUNKIN G. Conformal FDTD modeling of imperfect conductors at millimeter wave bands[J]. IEEE Transactions on Antennas and Propagation, 2011, 59(1): 199-205. DOI:10.1109/TAP.2010.2090490 |
| [6] | MALONEY J G, SMITH G S. The use of surface impedance concepts in the finite-difference time domain method[J]. IEEE Transactions on Antennas and Propagation, 1992, 40(1): 38-48. DOI:10.1109/8.123351 |
| [7] | BEGGS J H, LUEBBERS R J, YEE K S, et al. Finite difference time-domain implementation of surface impedance boundary conditions[J]. IEEE Transactions on Antennas and Propagation, 1992, 40(1): 49-56. DOI:10.1109/8.123352 |
| [8] | KARKKAINEN M K. FDTD surface impedance model for coated conductors[J]. IEEE Transactions on Electromagnetic Compatibility, 2004, 46(2): 222-233. DOI:10.1109/TEMC.2004.826891 |
| [9] | SHI L, YANG L, MA H, et al. Collocated SIBC FDTD method for coated conductors at oblique incidence[J]. Progress in Electromagnetics Research M, 2013, 30: 239-252. DOI:10.2528/PIERM13022512 |
| [10] | FELIZIANI M. Subcell FDTD modeling of field penetration through lossy shields[J]. IEEE Transactions on Electromagetic Compatibility, 2012, 54(2): 299-307. DOI:10.1109/TEMC.2011.2160982 |
| [11] | NAYYERI V, SOLEIMANI M, RAMAHI O. Modeling graphene in the finite-difference time domain method using a surface boundary condition[J]. IEEE Transactions on Antennas and Propagation, 2013, 61(8): 4176-4182. DOI:10.1109/TAP.2013.2260517 |
| [12] | BOUZIANAS G D, KANTARTZIS N V, YIOULTSIS T V, et al. Consistent study of graphene structures through the direct incorporation of surface conductivity[J]. IEEE Transactions on Magnetics, 2014, 50(2): 161-164. DOI:10.1109/TMAG.2013.2282332 |
| [13] | GUSTAVSEN B, SEMLYEN A. Rational approximation of frequency domain responses by vector fitting[J]. IEEE Transactions on Power Delivery, 1999, 14(3): 1052-1061. DOI:10.1109/61.772353 |
| [14] | KUNZ K S, LUEBBERS R J. The finite difference time domain method for electromagnetics[M]. Boca Raton: CRC Press, 1993. |
| [15] | WANG D, CHE W, CHANG Y, et al. A low-profile frequency selective surface with controllable tri-band characteristics[J]. IEEE Microwave and Wireless Components Letters, 2013, 12: 468-471. |
