Generalized Alternating-Direction Implicit Finite-Difference Time-Domain Method in Curvilinear Coordinate System *

In this paper, a novel approach is introduced towards an efficient Finite-Difference Time-Domain (FDTD) algorithm by incorporating the Alternating Direction Implicit (ADI) technique to the Nonorthogonal FDTD (NFDTD) method. This scheme can be regarded as an extension of the conventional ADI-FDTD scheme into a generalized curvilinear coordinate system. The improvement on accuracy and the numerical efficiency of the ADI-NFDTD over the conventional nonorthogonal and the ADI-FDTD algorithms is carried out by numerical experiments. The application in the modelling of the Electromagnetic Bandgap (EBG) structure has further demonstrated the advantage of the proposed method.


Introduction
The Finite-Difference Time-Domain (FDTD) Method [1] has been proven to be an effective algorithm in computational electromagnetics.The first FDTD algorithm was introduced by Yee [2] in 1966.Since then, extensions of this method [3][4][5][6][7][8] have been made continuously.The nonorthogonal FDTD (NFDTD) method [9][10][11] is a standard generalized FDTD algorithm based on the curvilinear coordinate system.As far as a curved geometry is concerned, the NFDTD algorithm has demonstrated improved efficiency over the conventional Yee's algorithm [12], due to the fact that considerably fewer cells are needed in the former by using conformal meshes instead of employing staircase approximation.
As is the case in the Cartesian FDTD algorithm, the time interval dt max used in the NFDTD method is constrained by the Courant-Friedrich-Levy (CFL) stability condition [13].Consider a two-dimensional NFDTD cell shown in Figure 1, ξ i and η j denote the skewed coordinates, and θ ij (i and j are cell indices) denotes the rotating angle from ξ i to η j .The Courant criterion for this particular cell (i, j) is donated as dt maxi,j .The time increment dt used in the NFDTD simulation should not exceed the Courant criterion dt max given by (1), which is the minimum value of dt maxi,j , in order to guarantee that the CFL condition is satisfied for every cell in the NFDTD meshes.The orthogonal CFL condition is a special case when θ = 90º.
For orthogonal meshes, it is straightforward to calculate dt max .However, with nonorthogonal meshes, it is not trivial to find the minimum value of dt maxi,j as diverse cells exist.More to the point, when the meshes contain globally large but locally very small or skewed cells, the using of the smallest dt maxi,j may result in low efficiency in the NFDTD simulations.
The alternating direction implicit (ADI) technique has been successfully applied to the orthogonal FDTD instead of the Yee's leapfrogging scheme, resulting in the removal of the CFL stability condition and an improvement in the efficiency of FDTD method.In [19], the ADI principle as applied in [6] and [7] is extended to the curvilinear coordinates based on Lee's NFDTD scheme and a novel NFDTD method that is free of the CFL stability condition is briefly introduced.Similar idea was proposed based on Holland's NFDTD formulation in [20].In this paper, the formulation of the generalized ADI-NFDTD as in [19] is presented in details and its numerical efficiency and accuracy is further studied.Compared to the conventional NFDTD method, the numerical efficiency can be improved by using an increased dt, and the late time instability of the NFDTD scheme has been greatly reduced.Compared with the conventional orthogonal ADI-FDTD, this generalized scheme is not unconditionally stable.However, it shows improved numerical efficiency in modelling curved structures by requiring less computer memory and computer run time, because coarser conformal grids are employed.In order to demonstrate its efficiency, the proposed ADI-NFDTD is applied in calculating the bandgap diagram of the Electromagnetic Bandgap structure (EBGs).

Formulation
The Maxwell curl equations can be written as the following for a source-free space: where μ is the medium permeability.Both of the equations can be casted into partial differential equations in the curvilinear coordinates as shown in [15].A 2-D TE wave module is used to illustrate the scheme.The formulation for the TM module can be obtained in a similar way.Two procedures are used to solve the differential equations.The first procedure is based on (4)- (6), and the second procedure is based on (7)- (9). 1) The First Procedure:   2) The Second Procedure: , 2   where E x , E y , H z are the covariant electric and magnetic field components which represent the flow of field along the grid, while D x , D y and H z are the contravariant electric fluxes and magnetic field component which represent the flow going through the facets of the grid.μ z is the medium permeability for calculating H z and g is the metric tensor computed from the local curvilinear coordinates.
The covariant field components (E x and E y ) must be calculated from the contravariant flux (D x and D y ) by interpolating equations as introduced in the conventional NFDTD scheme [9][10].As the covariant field compo-nents should be valued in the same position with their contravariant pairs, the neighbouring averaging projection scheme is employed.For a 2-D model, g xz and g yz are zeros, and g zz equals to unity for every cell.Hence we obtain ( 10)- (14).
n n y yy xy y y y where ε i (i = x, y) is the material permittivity in the x-or ydirection.By substituting ( 10), ( 11) and ( 13) into ( 6) and the resulting expression for 5), one can obtain (15).From (15) the cell indices (i,j) are used in the equations for a neat presentation, e.g., readily obtained by ( 4), ( 15) is a tri-diagonal system of equation that can be easily solved.
 can be calculated by (6).
The field components in the second procedure can be updated in a similar way, by the use of ( 11), ( 12), ( 14) and ( 7)-( 9).
The proposed ADI-NFDTD algorithm is an extension of the ADI-FDTD in the curvilinear coordinate system.As a result, if the grid is uniform and orthogonal, this algorithm will reduce to the conventional orthogonal ADI-FDTD, as illustrated in ( 16)- (21).
We apply the following relations to the derived ADI-NFDTD equations: , where ) the corresponding electric and magnetic field components in the Cartesian coordinates respectively.The g tensors have the relationship of ( 17) under 2-D orthogonal meshes.
xx yy xy xx yy g i j g i j g i j g i j g i j dx i j g i j dy i j Substituting ( 16)-( 17) into ( 15) yields (18), which is the ADI-FDTD updating equation under arbitrary orthogonal meshes.
Given that the FDTD grid is uniform, and denoting the spatial increments along the x-and y-axes as Δx and Δy, respectively, (19) are yielded from (17).
Then substituting ( 19) into ( 18) and multiplying on both side of the equation yields (20).Equation ( 20) is the general equation for orthogonal ADI-FDTD scheme with uniform grid.Equations ( 18) and ( 20) can treat homogeneous or inhomogeneous, isotropic or anisotropic media.In a medium where μ z (i, j) = μ z (i-1, j), (20) further reduces to the formula given in [6].

Comparison of the Late Time Instability and the Numerical Accuracy of the ADI-NFDTD Algorithm with the Conventional NFDTD Algorithm
The resonant modes of a cylindrical perfect electric conductor (PEC) cavity resonator are calculated by using both the conventional NFDTD and the proposed ADI-NFDTD schemes.The radius of the resonator is 0.15m and the structure is infinitely long.The cavity is filled with vacuum of permittivity ε r = 1.The outer material is PEC, and the whole computational region is 0.56 m × 0.56 m meshed by 25 × 26 cells.A sine modulated cosine pulse ( 21) is excited inside the cavity to provide a wide band excitation to excite all the possible TE modes: where A m relates to the amplitude of the signal, f is the frequency parameter, dt is the time increment and n is the iteration index.
The temporal signatures of the magnetic fields inside the resonator are recorded.The first 8,000 time steps results from both the conventional and the ADI-NFDTD algorithms are normalized by their own peak values and are plotted in Figure 2.
More simulations were performed using different dt values to study the stability and the accuracy of the ADI-NFDTD scheme.It is shown that unlike the conventional  dt dt E i j H i j H i j dx i j dy i j dx i j dy i j  ADI-FDTD, this generalized ADI-NFDTD scheme is not an unconditionally stable scheme.Instability in the temporal results is observed when small dt (i.e., comparable with or smaller than the dt max required in the NFDTD scheme) is used.However, compared to the NFDTD scheme, the unstable result in the ADI-NFDTD occurs much later and the spurious energy grows much slower than that in the conventional NFDTD scheme.This is demonstrated in Figure 3.
The removal of the CFL stability condition on dt in the ADI-NFDTD method is also observed.When dt reaches 16ps in the simulations, the conventional NFDTD cannot provide reasonable results, indicating that the CFL stability condition is violated, while ADI-NFDTD scheme can do so, even using a greater dt value.
In order to compare the accuracy of the ADI-NFDTD with the conventional NFDTD algorithm, the resonant frequencies are calculated using the Fast Fourier Transformation (FFT) and are compared with the analytical results.Since both the NFDTD and the ADI-NFDTD methods suffer numerical instability in late time steps, the temporal data are truncated to 40,000 and 100,000 time steps respectively.Such numbers are chosen in order to obtain the best frequency resolutions, and avoid too much spurious energy in the spectra at the same time.Since dt = 1 ps, the frequency resolutions of the conventional and the ADI-NFDTD results are 0.025 GHz and 0.01 GHz respectively.The calculated resonant frequency spectra are plotted in Figure 4.
As can be seen, it is possible to distinguish two closelyspaced frequency spectra from the results obtained by using the ADI-NFDTD, while it is hard to do so in the conventional NFDTD due to the limited frequency resolution.Moreover, the noise level of the ADI-NFDTD is lower than that of the conventional NFDTD, if dt used in the two simulations are comparable.In order to compare the accuracy of the two NFDTD schemes, the averaged relative error rate of the resonant frequency as a function of the relative time interval is introduced.This error rate is calculated in the following manner.Firstly, the absolute value of the difference of the calculated resonant frequency with the theoretical one is divided by the latter, and expressed in percentage, forming the relative error rate for each resonant mode.Then, the standard deviation of the relative error rates of all the modes within frequency band of interest (0-3 GHz) are calculated as the averaged relative error rate of the numerical results.The curve of the averaged relative error rate as the function of time interval is plotted in Fig- ure 5.
In theory, the ADI-NFDTD may not necessarily provide better accuracy than the conventional NFDTD.However, to the best of the authors' knowledge, the former can always provide more stable temporal results than the latter.As the FFT being performed with a sufficient frequency resolution, the ADI-NFDTD shows smaller relative error rate at some dt values in Figure 5.The removal of the CFL condition on NFDTD algorithm can be seen when dt is greater than 16 ps.However, according to the authors' simulation experience, dt should not exceed the corresponding Courant criterion dt max in order to maintain the accuracy in ADI-NFDTD simulations.By using the same dt, the computational efficiency is decreased in the ADI-NFDTD simulation than in the NFDTD one.

Comparison of the Numerical Efficiency and Accuracy of the ADI-NFDTD Algorithm with the ADI-FDTD Algorithm
A two-dimensional cavity resonator is modelled using both the ADI-NFDTD and the conventional orthogonal ADI-FDTD algorithms, in order to demonstrate the numerical efficiency improvement of the former scheme.The radius of the cavity is 0.15 m.The cavity is filled with vacuum of relative permittivity ε r = 1.The outer material is copper, with relative permittivity ε r = 1 and conductivity σ = 5.8 × 10 7 S/m.The whole computational region is 0.6 m × 0.6 m, meshed by the ADI-NFDTD method using 46 × 46 cells, and by the orthogonal ADI-FDTD method using 60 × 60, 70 × 70 and 80 × 80 cells respectively.The time step is chosen to be 2 ps and 25,000 iterations are run in each simulation.A modulated Gaussian pulse as shown in ( 22) is excited inside the cavity.The temporal responses of probed points inside the cavity are recorded and Fourier Transformed, after which the resonant frequencies is identified as peaks in the frequency spectra.
where A m is the maximum amplitude, f is the frequency parameter, n is the iteration index, dt is the time increment,   is the delay of the pulse in time step, and is the pulse half-duration at the 1/e point.In the simulation, A m =1 and f = 0.5 GHz.To calculate the relative error rates (RERs), the frequencies of all the genuine modes are compared with the theoretical results [22], and the results are shown in Ta- ble 1.
Table 2 presents a comparison of numerical errors, computer resources used in all simulations in the same computer under the same programming environment of Matlab 6.The averaged relative error rate calculated from the genuine and the spurious modes and that from only the genuine modes are denoted as 'ARER' and 'ARER-G', respectively.It can be seen that the ADI-NFDTD method can provide results with less numerical error while using less computer memory and less simulation run time.This demonstrates the accuracy and computing efficiency improvement of the ADI-NFDTD algorithm over the orthogonal ADI-FDTD.However, compared to the orthogonal ADI-FDTD algorithm, the ADI-NFDTD method is not an unconditionally stable method.The late time instability inherited from the NFDTD scheme is a major drawback of the ADI-NFDTD method.

Application of ADI-NFDTD Algorithm in Bandgap Structure Simulations
In this section, the ADI-NFDTD scheme is applied to model an electromagnetic bandgap (EBG) structure and the simulation results are compared with those from the NFDTD modelling.The two-dimensional EBG structure comprises of metallic rods periodically loaded in square lattice in free space.The unit cell approach [23] is used in this model.This EBG structure is infinite in the x-and the y-direction with a lattice constant (period) a.In the z-direction, the rods are infinitely long.Each rod is made from copper with relative permittivity ε r = 1 and conductivity σ = 5.8 × 10 7 S/m.The ratio of the radius r to the lattice constant a is chosen to be r/a = 0.2.The unit cell, the Brillouin Zone and the nonorthogonal mesh file used in both the ADI-NFDTD and NFDTD algorithms are shown in Figure 6.
The dispersion diagram calculated by NFDTD is plotted in Figure 7.In the NFDTD simulation, while calculating the frequency spectra for vectors near vector k = ΓM, the time domain signals exhibit instability quickly.This will cause difficulty in obtaining the dispersion diagram.However, this problem can be easily solved by using the ADI-NFDTD algorithm.Consequently, the frequency spectra obtained using the ADI-NFDTD show a reduced noise level and hence lead to more accurate solutions for identifying eigenmodes in the dispersion diagram, as shown in Figure 9.For instance, the mode at a normalized frequency of 1.36 (highlighted using the datatip in Figure 7) can be hardly detected from the NFDTD results as the NFDTD spectrum is heavily buried by the noise.However, this mode is clearly shown as a distinct peak in the ADI-NFDTD spectra.

Conclusions
In this paper, the ADI-FDTD scheme is extended into the nonorthogonal coordinate system to achieve an ADI-NFDTD algorithm.The stability, accuracy and efficiency of the proposed scheme are validated by the cylindrical resonator simulations.It is shown that the CFL condition of NFDTD is removed by the use of the ADI scheme.Secondly, although the ADI-NFDTD is not unconditionally stable, the unstable result occurs at a much later stage compared with the conventional NFDTD algorithm, which is a significant improvement in terms of late time instability.
In this way, the ADI-NFDTD is able to provide more precise frequency spectra with a smaller FFT resolution and a lower noise level than those obtained in the conventional NFDTD.This is of great benefit in distinguishing resonant or eigen-modes closely located in the frequency spectra.As a result, the ADI-NFDTD is quite suitable in the modelling of resonating or periodic structures, in which the energy is hardly dissipated, and in which the resonant or the eigen-modes are key factors to be examined in the study.
Without having to employ staircasing approximation, the curved structures or oblique surfaces are more accurately and efficiently modelled using the ADI-NFDTD scheme than using the orthogonal ADI-FDTD.
Although dt is not limited by the CFL condition and the conformal meshes can be coarser than those in the orthogonal FDTD, the increase of dt value and the decrease of the spatial resolution will result in decreased accuracy in the ADI-NFDTD results.Therefore, in ADI-NFDTD, the dt value and the mesh profile should be carefully chosen for the combined consideration of the efficiency and accuracy.

Figure 2 .
Figure 2. The normalized H field results of the first 8000 time steps from the two schemes

Figure 3 .Figure 4 .
Figure 3. Hz field temporal results of (a) the conventional NFDTD scheme, and (b) the ADI-NFDTD scheme, with time interval dt = 1 ps

Figure 5 .
Figure 5.The averaged relative error rate of the resonant frequency calculated by the conventional NFDTD and the proposed ADI-NFDTD schemes as the function of relative time interval dt/dt0 (dt0 = 1ps).After [19]

Figure 6 . 8 .
Figure 6.(a) The unit cell of the EBGs; (b) the Brillouin Zone in the reciprocal lattice with k = k 1 shown by a red vector; (c) the conformal mesh of the real lattice Take the eigen-mode calculation on one wave vector k = k1 (Figure 6(b)) for example.The temporal results at the same probe position obtained by using the ADI-NFDTD and the conventional NFDTD are compared in Figure 8.As can be seen from Figure 8(b), the ADI-NFDTD result shows resemblance with the NFDTD one at the beginning of simulations.It is noted that, after the first 50 nanoseconds, the results obtained via the ADI-NFDTD still remain stable while those from the NFDTD experience an exponential growth signaling numerical instability.

Figure 7 .Figure 8 .
Figure 7.The dispersion diagram of the square lattice EBG, TE mode results, calculated using the NFDTD method

Figure 9 .
Figure 9.The comparison of the frequency spectra in calculating k = k 1 using the conventional and the ADI-NFDTD algorithms