Abstract
In this paper, a novel algorithm based on the alternating direction implicit (ADI) multiresolution time-domain (MRTD) method for periodic structure simulation is proposed. By applying the multiresolution analysis in accordance with wavelet theory, the spatial sampling rate of the conventional finite-difference time-domain (FDTD) is significantly reduced by the MRTD method. The ADI method is then used to remove the Courant-Friedrich-Levy (CFL) limit that the MRTD method experiences. The periodic boundary condition (PBC) is directly implemented in the time domain using a constant transverse wave-number (CTW) wave. Numerical results are presented to confirm the efficiency and accuracy of the proposed method.
1. Introduction
Over the past decades, notable progress on the engineering applications of periodic structures has been achieved, especially on frequency selective surfaces (FSS) [1] and metamaterials [2, 3]. For an accurate analysis of electromagnetic scatter characteristics, many numerical methods have been developed. The finite-difference time-domain (FDTD) method is a typical algorithm that plays an important role in dealing with electromagnetic problems.
Normally, in the analysis of periodic structures using the FDTD method, only a unit cell is simulated and computed (instead of the whole structure) by incorporating the appropriate periodic boundary condition (PBC). However, unlike the situation of a normal incident wave, the implementation of the PBC for an oblique incident is complicated due to the time delay in the transverse plane. The split-field technique [4] is proposed for periodic structures analysis under oblique incidence circumstance. However, one has to adopt a multipart algorithm to solve the associated extra terms due to the transformed field. The spectral FDTD (SFDTD) method is a novel technique in which the constant transverse wave-number (CTW) wave is applied [5]. It is more efficient than the split-field technique due to the angle-independent stability criterion.
However, the time step is constrained by the Courant-Friedrich-Levy condition due to the explicit time-marching technique. To overcome this problem, the weakly conditionally stable finite-difference time-domain (WCS-FDTD) [6] and the locally one-dimensional finite-difference time-domain (LOD-FDTD) [7] have been widely applied to the study of periodic structures. Furthermore, the alternating direction implicit finite-difference time-domain (ADI-FDTD) method incorporates the advantages of both the explicit and implicit formats, namely, relatively simple calculation and unconditional stability [8]. Although the weighted Laguerre polynomials based spectral finite-difference time-domain (WLP-SFDTD) scheme for periodic structure analysis is efficient, solving the correspondingly huge banded system matrix is time-consuming [9]. One has to use the sparse lower-upper (LU) factorization packages or add a perturbation term.
In recent years, the multiresolution time-domain (MRTD) technique has been successfully developed due to its highly linear dispersion performance and ability to simulate complex electromagnetic structure [10–12]. Moreover, with highly linear dispersion performance, the MRTD scheme implies that a low sampling rate in space can still provide a relatively small phase error in the numerical simulation of a wave propagation problem. Therefore, it becomes possible that larger targets can be simulated without sacrificing accuracy. However, the MRTD scheme has a major drawback that the time stability condition is more rigorous than that of the FDTD scheme, which limits the computational efficiency of the MRTD algorithm.
In this paper, the spectral-FDTD (SFDTD) technique is applied to the conventional MRTD and ADI-MRTD methods, resulting in the PS-ADI-MRTD algorithms. The MRTD method and the ADI method are, respectively, used to reduce the spatial sampling rate and remove the KCL limit. The application, the SFDTD, is mainly using a constant transverse wave-number (CTW) wave to directly implement the PBC in the time domain. To verify the efficiency and accuracy of the PS-ADI-MRTD method, the numerical example is presented later.
2. Mathematical Formulation
2.1. Formulation of the Incident Wave
Considering that the CTW travels in the plane, when the incident wave is oblique, the expression of the PBC can be shown in the frequency domain as follows:
As shown in equation (1), if and are constant numbers, then will be a constant complex number. By transforming equation (1) into the time domain, we obtain the following:
Therefore, if and are constant numbers, the PBC will have no time delay. Then, the implementation of the PBC is similar to that of the normal incidence case.
As shown in Figure 1, the constant transverse wave-number (CTW) wave is employed as the incident wave in this discussion [5], and a TE wave plane is used here. The polarization direction of the TE wave is along is the incident wave vector, and is the incident angle. From Figure 1, we obtain the following:
Then,
By defining that the unit direction vector is in the incident plane and is parallel with the plane, then is expressed as follows:
The unit direction vector that is perpendicular to the incident plane is shown as follows:
The incident wave in the frequency domain is shown as follows:
Then, we obtain the magnetic field component of the incident wave (the transverse magnetic field): where is the impedance of the free space, is the reference plane of the incident wave, and is a Gaussian pulse used to limit the bandwidth of the wave.
By applying the inverse Fourier transform on (7) and (8), we obtain the expression of the CTW wave in the following time domain: where represents the inverse Fourier transform. Figure 1 indicates that , which means that if and are constant numbers, then is also constant. Then, it can be concluded that the different frequencies correspond to different incident angles in the CTW wave.
2.2. Handling of the PBC
For a better understanding of this paper, the conventional MRTD is briefly introduced here. The topology of the square patches periodic structure is depicted in Figure 2, where and are, respectively, the lengths along the x-direction and y-direction. and are the mesh points closest to the corresponding boundaries in Figure 2. The PBC of the MRTD can be expressed as follows: where is the effective support size of the basis function of the MRTD. Its value should be no greater than the number of total meshes for the minimum distance between the boundary and the edge of the scatter.
(a)
(b)
2.3. PS-ADI-MRTD Algorithm
The equations of the proposed PS-ADI-MRTD method are similar to the ADI-MRTD method [10]. The matrix form of Maxwell’s equations can be expressed as follows: where the coefficients matrix is as follows:
Matrix is a typical block periodic tri-diagonal matrix, where and are all square matrices. We use the D2 wavelet basis to expand the basic equation of PS-ADI-MRTD and . The elements of matrices and are as follows:
We define that and are square matrices, , and . Suppose that and . Matrix can then be written as where
According to Sherman-Morrison-Woodbury formula [13], one can have
Finally, can be solved
The element is calculated as follows. For the given reversible matrices and , where and , one can have
For the given and , one can have
The TF/SF boundary and absorb boundary condition of the proposed PS-ADI-MRTD is similar to the conventional MRTD.
The detail equations for the open loop case are omitted for saving the space.
3. Numerical Results
The structure of the thin metal rectangle array is shown in Figure 2, where the size of each metal rectangle is , and the thickness is assumed to be infinitely thin. The spatial step is chosen as , while the size of each unit cell is . The time step that satisfies the stability condition is set to and . In the computing progress of the PS-ADI-MRTD algorithm, CLFN is 3. The computing area is truncated with 8 layers of the CPML along the two directions of z-axis.
The TE wave of the CTW source is set as the excitation source. The wave vector is in the x-z plane and is represented by . Then, the incident electric field and magnetic field time domain expressions are as follows: where . Because the observed frequency range is set between 1~16 GHz, the value range of is 0~336.
Figure 3 shows the reflection coefficient of the electric field component as a function of frequency. The curves from the computing results of the SFDTD algorithms and commercial solver HFSS are also presented. As indicated in the figure, the curves of the three methods are basically the same from 2 to 12 GHz, which verifies the accuracy of the PS-ADI-MRTD algorithm. However, the consistency of results from the three methods is not quite well when the frequency is higher than 12 GHz; this situation even gets worse while . For a relatively fair comparison of the proposed method and HFSS, the “maximum length of element” in HFSS is set to be 1 mm (the special step of SFDTD and PS-ADI-MRTD are both 0.5 mm). The final CPU time of the SFDTD, the HFSS, and the PS-ADI-MRTD algorithms are 3255.9, 8421.7, and 2365.6 seconds, respectively. The computing efficiency of PS-ADI-MRTD is nearly 3.5 times of HFSS. This confirmed that the computational efficiency of the PS-ADI-MRTD algorithm is significantly higher than that of the SFDTD algorithm and HFSS.
(a)
(b)
As refer to the simulation about SFDTD algorithm, the implementation of SFDTD is mainly based on [5]. The numerical results for HFSS are obtained by the version HFSS 13.0. All calculations presented in this paper were performed on an Intel (R) Xeon (R) 2.80 GHz machine.
4. Conclusion
A novel algorithm referred to as the PS-AID-MRTD is presented to simulate the periodic structure. By using the CTW wave to handle the incident wave, one can directly implement the PBC. The PS-AID-MRTD algorithm has adopted each advantage of the MRTD method and the ADI technique. As a result, the spatial sampling rate is significantly reduced, and the time stability is guaranteed. Based on the numerical example, in the frequency range 2~12 GHz, the computing accuracy of the proposed method is basically the same with the SFDTD method and HFSS, and the computing effectiveness of the proposed method is superior to them.
Data Availability
The experimental data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported by China Postdoctoral Science Foundation (2016T90995) and the National Natural Science Foundation of Jiangsu Province (BK20150715).