# ADI-MRTD Algorithm for Periodic Structures Analysis

Authors: Xiaohua Dong; Yiwang Chen; Pin Zhang; Yichao Teng; Xiaobo Zhuansun
Journal: International Journal of Antennas and Propagation
 (2018)

---

## 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.


---

## Body

## 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 thex−y plane, when the incident wave is oblique, the expression of the PBC can be shown in the frequency domain as follows:
(1)ψ~x±Tx,y±Ty,z,ω=ψ~x,y,z,ωexp∓jkxTxexp∓jkyTy.As shown in equation (1), if kx and ky are constant numbers, then exp∓jkxTxexp∓jkyTy will be a constant complex number. By transforming equation (1) into the time domain, we obtain the following:
(2)ψx±Tx,y±Ty,z,t=ψx,y,z,texp∓jkxTxexp∓jkyTy.Therefore, ifkx and ky 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 Figure1, 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 t, k0 is the incident wave vector, and θ,φ is the incident angle. From Figure 1, we obtain the following:
(3)k0=kxx+kyy−kzz,kx=k0sinθcosφ,ky=k0sinθsinφ,kz=k0cosθ.Figure 1
Schematic diagram of the CTW incident wave.Then,(4)cosφ=kxkx2+ky2=kxkl,sinφ=kykl,kl=kx2+ky2,cosθ=kzk0,φ=tan−1kykx,θ=sin−1kx2+ky2k0.By defining that the unit direction vectorl is in the incident plane and is parallel with the x−y plane, then l is expressed as follows:
(5)l=cosφx+sinφy=kxklx+kykly.The unit direction vectort that is perpendicular to the incident plane is shown as follows:
(6)t=−sinφx+cosφy=−kyklx+kxkly.The incident wave in the frequency domain is shown as follows:(7)Ei=Eit=exp−jkxx−jkyy+jkzz−z0−k02σ2t=Et.Then, we obtain the magnetic field component of the incident wave (the transverse magnetic field):(8)Hl=Hicosθl=Eiη0cosθl=kzη0k0Eil=exp−jkxx−jkyy+jkzz−z0YTEexp−k02σ2l,where YTE=kz/η0k0, η0 is the impedance of the free space, kz=k02−kx2−ky2,k0=2πf/c, z0 is the reference plane of the incident wave, and exp−k02/σ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:
(9)EtCTW=exp−jkxx−jkyyF−1expjkzz−z0exp−k02σ2,HlCTW=1η0exp−jkxx−jkyyF−1kzk0expjkzz−z0exp−k02σ2,where F−1 represents the inverse Fourier transform. Figure 1 indicates that sinθ=kl/k0, which means that if kx and ky are constant numbers, then kl 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 Figure2, where Tx and Ty are, respectively, the lengths along the x-direction and y-direction. NxL, NyF, and NyL are the mesh points closest to the corresponding boundaries in Figure 2. The PBC of the MRTD can be expressed as follows:
(10)Exni+1/2,NyF±l,k=Exni+1/2,NyL±l,kexpjkyTy,EynNxF±l,j+1/2,k=EynNxL±l,j+1/2,kexpjkxTx,Ezni,NyF±l,k+1/2=Ezni,NyL±l,k+1/2expjkyTy,EznNxF±l,j,k+1/2=EznNxL±l,j,k+1/2expjkxTx,Hxn+1/2i,NyL±l+1/2,k+1/2=Hxn+1/2i,NyF±l+1/2,k+1/2exp−jkyTy,Hyn+1/2NxL±l+1/2,j,k+1/2=Hyn+1/2NxF±l+1/2,j,k+1/2exp−jkxTx,Hzn+1/2NxL±l+1/2,j+1/2,k=Hzn+1/2NxF±l+1/2,j+1/2,kexp−jkxTx,Hzn+1/2i+1/2,NxL±l+1/2,k=Hzn+1/2i+1/2,NxF±l+1/2,kexp−jkyTy,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.Figure 2
Geometric diagram of the rectangle patches periodic structure: (a) Top view:Tx = Ty = 15 mm, l = 12 mm, w = 3 mm; (b) side view: d = 6 mm.
(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:
(11)AX=Y,where the coefficients matrix A is as follows:
(12)A=a1b10⋯0tc2a2b2⋯000c3a3b3⋯0⋮0⋱⋱⋱⋮0⋮⋯cn−1an−1bn−1s00⋯cnan.MatrixA is a typical block periodic tri-diagonal matrix, where ai,bi,ci,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 a1,b1,c2,t, and s are as follows:
(13)a1=F1eG1eH1eI1eJ1eE2eF2eG2eH2eI2eD3eE3eF3eG3eH3eC4eD4eE4eF4eG4eB5eC5eD5eE5eF5e,b1=K1e0000J2eK2e000I3eJ3eK3e00H4eI4eJ4eK4e0G5eH5eI5eJ5eK5e,c2=A6eB6eC6eD6eE6e0A7eB7eC7eD7e00A8eB8eC8e000A9eB9e0000A10e,t=A1eB1eC1eD1eE1e0A2eB2eC2eD2e00A3eB3eC3e000A4eB4e0000A5eexp−jkyTy,s=Kn−4e0000Jn−3eKn−3e000In−2eJn−2eKn−2e00Hn−1eIn−1eJn−1eKn−1e0GneHneIneJneKneexpjkyTy.We define thatp1 and q1 are m square matrices, qn=p1−1t, and pn=sq1−1. Suppose that p=p10⋯0pnT and q=q10⋯0qnT. Matrix A can then be written as
(14)A=B+pqT,where
(15)B=A=a1−p1q1b10⋯00c2a2b2⋯000c3a3b3⋯0⋮0⋱⋱⋱⋮0⋮⋯cn−1an−1bn−1s00⋯cnan−pnqn.According to Sherman-Morrison-Woodbury formula [13], one can have
(16)A−1=B+pqT−1=B−1−B−1pIm+qTB−1p−1qTB−1.Finally,A−1 can be solved
(17)A−1=cijn×n,cij=gihj−ψiζj,i≤j,xiyj−ψiζj,i>j.The elementcij is calculated as follows. For the given reversible matrices p1 and q1, where qn=p1−1t and pn=sq1−1, one can have
(18)u1=a1−p1q1,li−1=ciui−1−1,ui=ai−li−1bi−1,i=2,3,…,n−1,un=an−pnqn−ln−1bn−1,αn=an−pnqn,βi=αi+1−1ci+1,αi=αi−biβi,i=n−1,n−2,…,1.For the giveng1=Im and x1=Im, one can have
(19)h1=αi−1,hi=−bi−1hi−1αi−1,i=2,3,…,n,gn=unhn−1,gi=−bigi+1ui−1,i=n−1,n−2,…,2,y1=α1−1,yi=−βi−1yi−1,i=2,3,…,n,xn=unhn−1,xi=−lixi+1,i=n−1,n−2,…,2,ξ=Im+q1g1h1+qnynx1p1+q1g1hn+qnynxnpn,ψ=ξ−1p1x1y+pnyng,ζ=q1g1hT+qnynxT.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.


## 2.1. Formulation of the Incident Wave

Considering that the CTW travels in thex−y plane, when the incident wave is oblique, the expression of the PBC can be shown in the frequency domain as follows:
(1)ψ~x±Tx,y±Ty,z,ω=ψ~x,y,z,ωexp∓jkxTxexp∓jkyTy.As shown in equation (1), if kx and ky are constant numbers, then exp∓jkxTxexp∓jkyTy will be a constant complex number. By transforming equation (1) into the time domain, we obtain the following:
(2)ψx±Tx,y±Ty,z,t=ψx,y,z,texp∓jkxTxexp∓jkyTy.Therefore, ifkx and ky 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 Figure1, 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 t, k0 is the incident wave vector, and θ,φ is the incident angle. From Figure 1, we obtain the following:
(3)k0=kxx+kyy−kzz,kx=k0sinθcosφ,ky=k0sinθsinφ,kz=k0cosθ.Figure 1
Schematic diagram of the CTW incident wave.Then,(4)cosφ=kxkx2+ky2=kxkl,sinφ=kykl,kl=kx2+ky2,cosθ=kzk0,φ=tan−1kykx,θ=sin−1kx2+ky2k0.By defining that the unit direction vectorl is in the incident plane and is parallel with the x−y plane, then l is expressed as follows:
(5)l=cosφx+sinφy=kxklx+kykly.The unit direction vectort that is perpendicular to the incident plane is shown as follows:
(6)t=−sinφx+cosφy=−kyklx+kxkly.The incident wave in the frequency domain is shown as follows:(7)Ei=Eit=exp−jkxx−jkyy+jkzz−z0−k02σ2t=Et.Then, we obtain the magnetic field component of the incident wave (the transverse magnetic field):(8)Hl=Hicosθl=Eiη0cosθl=kzη0k0Eil=exp−jkxx−jkyy+jkzz−z0YTEexp−k02σ2l,where YTE=kz/η0k0, η0 is the impedance of the free space, kz=k02−kx2−ky2,k0=2πf/c, z0 is the reference plane of the incident wave, and exp−k02/σ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:
(9)EtCTW=exp−jkxx−jkyyF−1expjkzz−z0exp−k02σ2,HlCTW=1η0exp−jkxx−jkyyF−1kzk0expjkzz−z0exp−k02σ2,where F−1 represents the inverse Fourier transform. Figure 1 indicates that sinθ=kl/k0, which means that if kx and ky are constant numbers, then kl 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 Figure2, where Tx and Ty are, respectively, the lengths along the x-direction and y-direction. NxL, NyF, and NyL are the mesh points closest to the corresponding boundaries in Figure 2. The PBC of the MRTD can be expressed as follows:
(10)Exni+1/2,NyF±l,k=Exni+1/2,NyL±l,kexpjkyTy,EynNxF±l,j+1/2,k=EynNxL±l,j+1/2,kexpjkxTx,Ezni,NyF±l,k+1/2=Ezni,NyL±l,k+1/2expjkyTy,EznNxF±l,j,k+1/2=EznNxL±l,j,k+1/2expjkxTx,Hxn+1/2i,NyL±l+1/2,k+1/2=Hxn+1/2i,NyF±l+1/2,k+1/2exp−jkyTy,Hyn+1/2NxL±l+1/2,j,k+1/2=Hyn+1/2NxF±l+1/2,j,k+1/2exp−jkxTx,Hzn+1/2NxL±l+1/2,j+1/2,k=Hzn+1/2NxF±l+1/2,j+1/2,kexp−jkxTx,Hzn+1/2i+1/2,NxL±l+1/2,k=Hzn+1/2i+1/2,NxF±l+1/2,kexp−jkyTy,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.Figure 2
Geometric diagram of the rectangle patches periodic structure: (a) Top view:Tx = Ty = 15 mm, l = 12 mm, w = 3 mm; (b) side view: d = 6 mm.
(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:
(11)AX=Y,where the coefficients matrix A is as follows:
(12)A=a1b10⋯0tc2a2b2⋯000c3a3b3⋯0⋮0⋱⋱⋱⋮0⋮⋯cn−1an−1bn−1s00⋯cnan.MatrixA is a typical block periodic tri-diagonal matrix, where ai,bi,ci,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 a1,b1,c2,t, and s are as follows:
(13)a1=F1eG1eH1eI1eJ1eE2eF2eG2eH2eI2eD3eE3eF3eG3eH3eC4eD4eE4eF4eG4eB5eC5eD5eE5eF5e,b1=K1e0000J2eK2e000I3eJ3eK3e00H4eI4eJ4eK4e0G5eH5eI5eJ5eK5e,c2=A6eB6eC6eD6eE6e0A7eB7eC7eD7e00A8eB8eC8e000A9eB9e0000A10e,t=A1eB1eC1eD1eE1e0A2eB2eC2eD2e00A3eB3eC3e000A4eB4e0000A5eexp−jkyTy,s=Kn−4e0000Jn−3eKn−3e000In−2eJn−2eKn−2e00Hn−1eIn−1eJn−1eKn−1e0GneHneIneJneKneexpjkyTy.We define thatp1 and q1 are m square matrices, qn=p1−1t, and pn=sq1−1. Suppose that p=p10⋯0pnT and q=q10⋯0qnT. Matrix A can then be written as
(14)A=B+pqT,where
(15)B=A=a1−p1q1b10⋯00c2a2b2⋯000c3a3b3⋯0⋮0⋱⋱⋱⋮0⋮⋯cn−1an−1bn−1s00⋯cnan−pnqn.According to Sherman-Morrison-Woodbury formula [13], one can have
(16)A−1=B+pqT−1=B−1−B−1pIm+qTB−1p−1qTB−1.Finally,A−1 can be solved
(17)A−1=cijn×n,cij=gihj−ψiζj,i≤j,xiyj−ψiζj,i>j.The elementcij is calculated as follows. For the given reversible matrices p1 and q1, where qn=p1−1t and pn=sq1−1, one can have
(18)u1=a1−p1q1,li−1=ciui−1−1,ui=ai−li−1bi−1,i=2,3,…,n−1,un=an−pnqn−ln−1bn−1,αn=an−pnqn,βi=αi+1−1ci+1,αi=αi−biβi,i=n−1,n−2,…,1.For the giveng1=Im and x1=Im, one can have
(19)h1=αi−1,hi=−bi−1hi−1αi−1,i=2,3,…,n,gn=unhn−1,gi=−bigi+1ui−1,i=n−1,n−2,…,2,y1=α1−1,yi=−βi−1yi−1,i=2,3,…,n,xn=unhn−1,xi=−lixi+1,i=n−1,n−2,…,2,ξ=Im+q1g1h1+qnynx1p1+q1g1hn+qnynxnpn,ψ=ξ−1p1x1y+pnyng,ζ=q1g1hT+qnynxT.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 Figure2, 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 step that satisfies the stability condition is set to ΔtCFL=Δ/3c and CFLN=Δt/ΔtCFL. 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 byky=0. Then, the incident electric field and magnetic field time domain expressions are as follows:
(20)EyCTW=exp−jkxxF−1expjkzz−z0−k02σ2−jt0k0c,HxCTW=1η0exp−jkxxF−1kzk0expjkzz−z0−k02σ2−jt0k0c,where σ=950. Because the observed frequency range is set between 1~16 GHz, the value range of kx is 0~336.Figure3 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 θ=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.Figure 3
Frequency characteristics at different angles of incidence.
(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.



---

*Source: 9562780-2018-10-10.xml*