ADI-MRTD Algorithm for Periodic Structures Analysis

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.


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 timedomain (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][11][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.

Mathematical Formulation
2.1.Formulation of the Incident Wave.Considering that the CTW travels in the x − y 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 k x and k y are constant numbers, then exp ∓jk x T x exp ∓jk y T y will be a constant complex number.By transforming equation ( 1) into the time domain, we obtain the following: Therefore, if k x and k y 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 wavenumber (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 t, k 0 is the incident wave vector, and θ, φ is the incident angle.From Figure 1, we obtain the following: By defining that the unit direction vector l is in the incident plane and is parallel with the x − y plane, then l is expressed as follows: The unit direction vector t 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 Y TE = k z /η 0 k 0 , η 0 is the impedance of the free space, , z 0 is the reference plane of the incident wave, and exp −k 2 0 /σ 2 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 F −1 represents the inverse Fourier transform.Figure 1 indicates that sin θ = k l /k 0 , which means that if k x and k y are constant numbers, then k l 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 T x and T y are, respectively, the lengths along the and y-direction.Nx L , Ny F , and Ny L are the mesh points closest to the corresponding boundaries in Figure 2. The PBC of the MRTD can be expressed as follows: Nx F ±l+1/2,j,k+1/2 exp −jk x T x , i+1/2,Nx F ±l+1/2,k exp −jk y T y , 10 where l = 0, 1, 2 … 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.

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 A is as follows: Matrix A is a typical block periodic tri-diagonal matrix, where a i , b i , c i , t, and s are all m square matrices.We use the D2 wavelet basis to expand the basic equation of PS-ADI-MRTD and m = 5.The elements of matrices a 1 , b 1 , c 2 , t, and s are as follows: We define that p 1 and q 1 are m square matrices, q n = p −1 1 t, and p n = sq −1 1 .Suppose that p = p 1 0 ⋯ 0 p n T and q = q 1 0 ⋯ 0 q n T .Matrix A can then be written as where According to Sherman-Morrison-Woodbury formula [13], one can have Finally, A −1 can be solved The element c ij is calculated as follows.For the given reversible matrices p 1 and q 1 , where q n = p −1 1 t and p n = sq −1 1 , one can have For the given g 1 = I m and x 1 = I m , 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.

Numerical Results
The structure of the thin metal rectangle array is shown in Figure 2, where the size of each metal rectangle is 12mm × 3mm, and the thickness is assumed to be infinitely thin.The spatial step is chosen as Δx = Δy = Δz = Δ = 0 5mm, while the size of each unit cell is 15mm × 15mm.The time 4 International Journal of Antennas and Propagation step that satisfies the stability condition is set to Δt CFL = Δ/3c and CFLN = Δt/Δt CFL .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 k y = 0.Then, the incident electric field and magnetic field time domain expressions are as follows: where σ = 950.Because the observed frequency range is set between 1~16 GHz, the value range of k x 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 the three methods is not quite well when the frequency is higher than 12 GHz; this situation even gets worse while θ = 60 °.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.
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.

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.

Figure 1 :
Figure 1: Schematic diagram of the CTW incident wave.