-
摘要:
针对周期结构电磁特性参数求解问题，介绍了一种基于网络分析法、矢量拟合法，用来快速求解低剖面周期结构电磁特性参数的内部阻抗边界条件（INBC）与时域有限差分（FDTD）结合的INBC-FDTD计算方法。该方法将金属层的二端口频域阻抗参数曲线先通过矢量拟合法进行有理分式拟合，再对其进行时域变换后嵌入FDTD公式完成对电场、磁场的更新工作。所提方法完整地考虑了在金属层传输的电磁场，其二端口网络阻抗参数全面地考虑了端口之间的互耦问题。
-
关键词:
- 频率选择表面(FSS) /
- 矢量拟合法 /
- 内部阻抗边界条件(INBC) /
- 周期结构 /
- 时域有限差分(FDTD)
Abstract:In order to solve the problem of electromagnetic characteristic parameters of periodic structures, an INBC-FDTD calculation method based on network analysis method and vector fitting method is:introduced in this paper, which can be used to quickly solve the electromagnetic characteristic parameters of low-profile periodic structures.In this method, the two-port frequency domain impedance parameter curve of the metal layer is firstly fitted with rational fractions by the vector fitting method, and then the time domain transform is carried out to embed the FDTD formula to complete the update of the electric field and the magnetic field. The proposed method fully considers electro-magnetic fields transmitted in the metal layer, and its two-port network impedance parameter fully considers the mutual coupling between the ports.
-
采用基于传统均匀网格的数值计算方法对含金属薄膜的大尺寸低剖面周期结构进行处理时往往面临着一大挑战，即应设置足够小的网格来捕捉导电薄层内的场变化。细微的离散化网格产生了大量的计算数据，也直接导致了较长的运算时间，同时，其还对计算机的内存和处理速度都提出了较高要求。目前，国内外研究机构已经提出了一些方法[1-12]来解决上述问题：①将亚网格技术用于薄层的细微网格建模[1-2]。但是，为保证方法的稳定性，离散时间步长仍受限制，这依然会导致大量的计算时间。②共形时域有限差分(FDTD)方法将导电薄层视为完美导体(PEC)，薄层上的切向场被设置为零[3-5]。在这种方法中，忽略导体损耗的这一处理可能会导致计算错误，尤其当结构处于高频电磁场时，导体将产生明显损耗。③采用表面阻抗边界条件(SIBC)直接分析金属薄层，可避免细微FDTD网格在薄层内部离散化的过程[6-9]。但在这种方法中，需要复杂的卷积处理和高精度的有理逼近运算。另外，Feliziani[10]提出了解决屏蔽问题的亚元胞FDTD模型，其优点是使用电路理论而不是使用卷积方程来求解屏蔽的场方程，虽然有效解决了屏蔽问题，但是对于复杂的结构(如多层板)，仍需要高精度的有理近似。Nayyeri[11]和Bouzianas[12]等提出了对石墨烯结构进行数值计算的方法，对石墨烯表面电导率方程进行直接分析，通过拉普拉斯变换和中心微分法将其结合到FDTD方法中进行建模运算。
本文提出了一种精确有效的方法来处理低剖面周期结构。通过拉普拉斯变换和中心微分法，将周期结构的Y矩阵参数直接嵌入到电场的更新方程中。与早期的工作相比，该方法可以对整体结构直接采用粗网格建模，时间步将不再受金属层厚度的限制；可以避免复杂的频域至时域的卷积运算，相关的频域函数可以由矢量拟合法轻松获得；计算时间快，精度高。
1. INBC-FDTD步进公式推导
图 1给出了低剖面周期结构置于两FDTD网格Γ-和Γ+之间的侧视图。图中：n-和n+分别为由Γ-指向Γ+的向量和由Γ+指向Γ-的向量。假设该结构沿着z方向上的厚度远小于所关注的最高频段下的波长和结构的横向尺寸，则结构切向电场和切向磁场之间的关系可通过引入如下Y矩阵得到：
(1a) 
(1b) 
图 1 置于两FDTD网格Γ-和Γ+之间的低剖面周期结构侧视图Figure 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 二维周期结构的计算空间示意图Figure 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种方法来说具有较快的计算速度。
表 1 提取得到的十字贴片型FSS结构Y矩阵参数频域函数的等效有理式参数：极点pn及留数rnTable 1. Extracted values of poles (pn) and residues (rn) for Y-matrix coefficients of cross patched FSSn 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矩阵参数频域函数的等效有理式参数：dTable 2. Extracted values of d for Y-matrix coefficients of cross patched FSSY d Y11 0.050 7 Y12 0.006 4 Y21 0.006 4 Y22 0.011 3 
图 4 十字贴片型FSS的反射波和透射波时域波形Figure 4. Reflected wave and transmitted wave of cross patched FSS in time domain
图 5 不同方法得出的频域散射参数幅度值Figure 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及其单元结构示意图Figure 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及留数rnTable 4. Extracted values of poles (pn) and residues (rn) for Y-matrix coefficients of complementary FSS with square ringn 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矩阵参数频域函数的等效有理式参数：dTable 5. Extracted values of d for the Y-matrix coefficients of complementary FSS with square ringY d Y11 0.062 4 Y12 0.016 9 Y21 0.016 9 Y22 0.144 6 
图 7 方环形结构互补型FSS的透射波时域波形Figure 7. Transmitted wave of complementary FSS with square ring in time domain
图 8 透射系数的频域波形比较Figure 8. Comparison of transmission coefficient waves in frequency domain表 6 两种方法计算参数的比较Table 6. Comparison of computational parameters for two methods方法 Δt 单元个数(长×宽×高) 时间步 CPU时间 亚网格
FDTD0.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 置于两FDTD网格Γ-和Γ+之间的低剖面周期结构侧视图
Figure 1. Side view of low-profile periodic structure placed between two FDTD grids Γ-and Γ+
图 2 二维周期结构的计算空间示意图
Figure 2. Schematic diagram of computational space for two-dimensional periodic structure
图 4 十字贴片型FSS的反射波和透射波时域波形
Figure 4. Reflected wave and transmitted wave of cross patched FSS in time domain
图 5 不同方法得出的频域散射参数幅度值
Figure 5. Amplitude of scattering coefficients in frequency domain obtain by different methods
图 6 方环形结构的互补型FSS及其单元结构示意图
Figure 6. Schematic diagram of complementary FSS with square ring and its unit cell
图 7 方环形结构互补型FSS的透射波时域波形
Figure 7. Transmitted wave of complementary FSS with square ring in time domain
图 8 透射系数的频域波形比较
Figure 8. Comparison of transmission coefficient waves in frequency domain
表 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 
下载: 导出CSV
表 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
下载: 导出CSV
表 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
下载: 导出CSV
表 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
下载: 导出CSV
表 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
下载: 导出CSV
表 6 两种方法计算参数的比较
Table 6. Comparison of computational parameters for two methods
方法 Δt 单元个数(长×宽×高) 时间步 CPU时间 亚网格
FDTD0.42 ps(粗网格)
16.7 fs(细网格)41×41×80 3 000 14 h INBC-FDTD 0.42 ps 41×41×80 3 000 17 min
下载: 导出CSV
-
[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. 期刊类型引用(0)
其他类型引用(2)
-
下载:
下载:
点击查看大图
计量
- 文章访问数: 473
- HTML全文浏览量: 65
- PDF下载量: 190
- 被引次数: 2
