A Modified AH-FDTD Unconditionally Stable Method Based on High-Order Algorithm

Abstract

The unconditionally stable method, Associated-Hermite FDTD, has attracted more and more attentions in computational electromagnetic for its time-frequency compact property. Because of the fewer orders of AH basis needed in signal reconstruction, the computational efficiency can be improved further. In order to further improve the accuracy of the traditional AH-FDTD, a high-order algorithm is introduced. Using this method, the dispersion error induced by the space grid can be reduced, which makes it possible to set coarser grid. The simulation results show that, on the condition of coarse grid, the waveforms obtained from the proposed method are matched well with the analytic result, and the accuracy of the proposed method is higher than the traditional AH-FDTD. And the efficiency of the proposed method is higher than the traditional FDTD method in analysing 2D waveguide problems with fine-structure.

[ProQuest: [...] denotes non US-ASCII text; see PDF]

Zheng Pan 1; 2 and Zhikai Fu 1 and Lihua Shi 1 and Zhengyu Huang 1 and Zheng Sun 1

Academic Editor:Elisa Francomano

1, National Key Laboratory on Electromagnetic Environmental Effects and Electro-Optical Engineering, PLA University of Science and Technology, Nanjing 210007, China

2, Luoyang Hydraulic Engineering Technical Institute, Luoyang 471023, China

Received 8 October 2016; Revised 20 December 2016; Accepted 16 January 2017; 13 February 2017

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The unconditionally stable method is an important research field in computational electromagnetics. It has tremendous superiority, because it can break through the restraint of Courant-Friedrichs-Levy (CFL) stable condition and improve the efficiency of FDTD method greatly. In the past thirty years, a series of unconditionally stable methods spring up in computational electromagnetic field, such as Crank-Nicolson FDTD (CN-FDTD) method [1], Alternating-Direction Implicit FDTD (ADI-FDTD) method [2], and Weighted Laguerre polynomials FDTD (WLP-FDTD) method [3].

The CN-FDTD method adopted the average way and took the average value of n+1th time step and the nth time step as the value of the n+1/2th time step. In general, the former iteration calculation problem can be converted into the solution of implicit equations, which breaks the limitation of the stable condition. But this method needs to solve the large matrix equations at each time step and consume huge memory, which makes the application of this method limited. ADI-FDTD can be seen as an improved algorithm based on CN-FDTD. The main idea of ADI-FDTD is to decouple the whole matrix equations into different dimension and convert the solution of the large matrix equations in the whole area into the solution of equations in each dimension, which reduces the dimensionality of matrix greatly and improves the efficiency. However, ADI-FDTD also brings one problem: the reduction of matrix dimension is obtained by abandoning the fixed term and the result is approximate. When the time step or the space step is big, the simulation results will suffer severe errors.

Chung et al. [3] put forward a new method based on weighted Laguerre polynomials orthogonal basis, which firstly expands the time-domain Maxwell equations using weighted Laguerre polynomials orthogonal basis and then eliminates every orthogonal term by Galerkin method so as to obtain the coefficients matrix equation. Finally, each point waveform in the spatial domain can be obtained by solving the coefficients matrix equation. The accuracy can be guaranteed for the signal reconstruction is based on the orthogonal basis. Later, WLP-FDTD method achieved comparatively large improvements in efficiency and boundary conditions [4-6]. In 2014, Huang et al. [7] proposed an unconditionally stable method based on AH basis according to the most compactly supported characteristics of Associated-Hermite basis, which can reconstruct signals by fewer orders so as to increase the efficiency. In 2015, Huang et al. [8] conducted the eigenvalue decomposition on the former AH matrix and decoupled each order basis. Finally, the matrix solution can be converted into the solution of order-based submatrix equations, which increases the efficiency further and decreases the CPU memory.

Aiming at the problems of huge memory consumption and severe numerical dispersion in electrically large object simulation with conventional FDTD method, Fang firstly proposed a high-order FDTD method [9], which adopts high-order discrete scheme in space and normal FDTD scheme in time. This solution can maintain better dispersion characteristics and calculation accuracy in coarse mesh grid. Afterwards, this method is deeply researched on numerical dispersion, stability, and high-order PML absorbing boundary conditions and is further applied in analysis of waveguide field, resonance problem and electromagnetic compatibility, and so on.

This article applies the high-order algorithm to the unconditionally stable AH-FDTD method and conducts high-order discrete format in space. The accuracy can be improved and the mesh grid in non-fine-structure region can be dissected to be relatively coarse. At last, two 2D waveguide simulation models were established, which proved the validity of the proposed method.

2. The Associated-Hermite Orthogonal Basis

For the most compactly supported characteristics of AH orthogonal basis, it has a wide application in signal processing, image analysis, biological engineering, and so on. The AH orthogonal basis function can be expressed as [figure omitted; refer to PDF] where _Hn (t~)=^(-1)net~2 (^dn /d^tn )(^e-t~2 ) indicates nth order Hermite polynomial and t~=(t-_Tf )/σ is the time variable. _Tf and σ are time-shift factor and time-scale factor separately. The first four orders of AH orthogonal basis are shown in Figure 1.

Figure 1: The first four orders of AH orthogonal basis.

[figure omitted; refer to PDF]

The differential expression of the AH orthogonal basis can be written as [figure omitted; refer to PDF]

For a signalf(t)=^∑n=0∞_{αn[varphi]¯n} (t) expanded with AH orthogonal basis, the first-order time derivative of f(t) can be expressed as [figure omitted; refer to PDF]

3. AH-FDTD Method Based on High-Order Algorithm

3.1. Maxwell Equations Expanded by AH Basis

Taking into account the lossless media, the time-domain Maxwell equations of the two-dimensional TEz model are [figure omitted; refer to PDF] where [straight epsilon] is capacitivity and μ is permeability.

Substituting the time derivative term in (4) with (3), each order of the AH orthogonal basis can be eliminated using Galerkin method. By simplifying, it can be expressed as [figure omitted; refer to PDF]

Discretizing (5) in space domain, the coefficient matrix equations of the discrete point (i,j) can be written as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and _(β)N×N is the unit matrix. [figure omitted; refer to PDF] _{(Ex )i,j} ,_{(Ey )i,j} , and _{(Hz )i,j} , respectively, represent the column vector of each field composed by AH basis expansion coefficients. Substituting (6)-(7) into (8) and gathering all points in space domain, the matrix equation on magnetic field can be obtained eventually. [figure omitted; refer to PDF] where [A] is the five-diagonal band nested matrix and {[J]} is the source term. By solving the matrix equation, the field value of each point can be obtained in the space domain. Figure 2 is the diagram of the five-diagonal band matrix.

Figure 2: The five-diagonal banded coefficient matrix [A].

[figure omitted; refer to PDF]

3.2. The High-Order Algorithm

In rectangular coordinate system, we set f(x,y,z,t) as a vector of electric field or magnetic field. The difference format with 2n order accuracy in space domain can be expressed as [figure omitted; refer to PDF]

Expanding f[x+(Δx/2)(2n-1)] and f[x-(Δx/2)(2n-1)] by Taylor series, (13) can be obtained using the method of undetermined coefficients. [figure omitted; refer to PDF]

When N=2, the fourth-order accuracy can be achieved in the space domain. Here, _a1 =9/8 and _a2 =-1/24, at the same time substituting the result obtained from (12) and (13) in (5). The fourth-order space accuracy AH-FDTD format can be obtained: [figure omitted; refer to PDF] where the variables and parameters are in accordance with Section 3.1. Assembling the equations in the whole space, the matrix equation of magnetic field can eventually be expressed as [figure omitted; refer to PDF] where [A-] is the updated banded matrix and ([J]) is the source term. By solving the equation, the field value of each point in the space domain can be obtained. Figure 3 shows the updated matrix diagram.

Figure 3: The banded coefficient matrix [A-].

[figure omitted; refer to PDF]

4. Simulation

4.1. Case 1

In order to verify the effectiveness of the proposed method, a two-dimensional parallel plate waveguide model is set up. In consideration of the convenience of comparing the numerical results with analytical results, the parallel waveguide is set to be filled with the homogeneous medium, and its size is 1.2 m × 0.08 m. Its permeability is taken as 4π×^10-7 H/m and the capacitance is taken as 8.854187817×^10-12 F/m, which is the same as the free space. In simulation, the whole region is divided into uniform grid of 100 × 8, and _P1 and _P2 are two observation points. A sinusoidal-modulated Gaussian pulse is chosen as y direction excitation source, which can be expressed as [figure omitted; refer to PDF] where _td =1/(2_fc ), _tc =3_td , and _fc = 0.8 GHz. The total time duration is _Ts = 10.81 ns. Figure 4 is the sketch map of 2D parallel plate waveguide.

Figure 4: The 2D parallel plate waveguide.

[figure omitted; refer to PDF]

4.1.1. Result Analysis

The parallel waveguide model is calculated, respectively, with AH-FDTD method and high-order AH-FDTD method. From the simulation, the results of AH-FDTD method and its modified method are both matched well with the analytical solution at points _P1 and _P2 . Figure 5 shows the waveform of _P1 , where we can see that the numerical dispersion is smaller using the high-order AH-FDTD method, and its waveform is more close to the analytical results. Figure 6 shows the waveform of _P2 , and the numerical dispersion is smaller using the high-order AH-FDTD method.

Figure 5: The waveform of _Ey on point _P1 .

[figure omitted; refer to PDF]

Figure 6: The waveform of _Ey on point _P2 .

[figure omitted; refer to PDF]

Figure 7 shows the relative error of the results of AH-FDTD and high-order AH-FDTD to the analytic solution. Among them, Figure 7(a) gives the relative error results of _P1 , and Figure 7(b) shows the result of point _P2 . The relative error is defined by R=20_log10 [...]((^Eya (t)-^Ey (t))/max[...](^Eya (t))), where ^Ey (t) is Y direction electric field obtained by AH-FDTD method or the modified AH-FDTD method, and ^Eya (t) is Y direction electric field obtained by the analytical solution. From Figure 7(a), we can see that their error is small, and the relative error is under -40 dB. At the same time, the relative error of the proposed method in this paper is smaller.

Figure 7: The relative error between AH-FDTD and high-order AH-FDTD.

(a) _{P 1}

[figure omitted; refer to PDF]

(b) _{P 2}

[figure omitted; refer to PDF]

From Figure 7(b), we can see that the numerical dispersion is more serious than _P1 because of the greater distance between point _P2 and the excitation source, and the relative error is below -33 dB. The relative error of the proposed method in this paper is smaller.

4.2. Case 2

In this case, a two-dimensional parallel plate waveguide with thin dielectric layer is investigated. The parallel waveguide size is 1.2 m × 0.08 m and a thin dielectric layer is embedded in the middle region shown in Figure 8. In order to clearly show the feature of the proposed algorithm, the relative permittivity of the thin layer is 2 and the width is set to three types: 2.4 × 10^-2 m, 2.4 × 10^-3 m, and 2.4 × 10^-5 m. The minimum cell size is 0.012 m × 0.01 m, 0.0012 m × 0.01 m, and 0.000012 m × 0.01 m, respectively. In the simulation, the nonuniform grid is used, and the whole calculation area is divided into 100 × 8, 108 × 8 and 110 × 8 grid. _P1 and _P2 are the two observation points, and other settings are shown in Section 4.1.

Figure 8: The 2D waveguide with dielectric material embedded in the middle region.

[figure omitted; refer to PDF]

4.2.1. Result Analysis

The simulation model is calculated, respectively, by traditional FDTD method and modified AH-FDTD method. Considering the similarity of these three results, we just show the simulation result under the width d=2.4×^10-2 m and Figures 9, 10, and 11 show the comparison result between traditional FDTD method and modified AH-FDTD method. As Figures 9 and 10 show, the results of traditional FDTD method and the modified AH-FDTD method are matched well with each other. Figure 9 shows the waveform of _Ey at _P1 , where the reflected wave from the dielectric interface can be clearly seen. And Figure 10 is the waveform of _Ey at _P2 , which displays mainly the transmission wave.

Figure 9: The waveform of _Ey on point _P1 .

[figure omitted; refer to PDF]

Figure 10: The waveform of _Ey on point _P2 .

[figure omitted; refer to PDF]

Figure 11: The relative error between traditional FDTD and modified AH-FDTD.

[figure omitted; refer to PDF]

The relative error between traditional FDTD and high-order AH-FDTD is shown in Figure 11. The relative error size is defined by R=20_log10 [...]((^Eya (t)-^Ey (t))/max[...](^Eya (t))), where ^Ey (t) is y direction electric field obtained by the modified AH-FDTD method and ^Eya (t) is y direction electric field obtained by the traditional FDTD method. From Figure 11, we can see that the relative errors of _P1 and _P2 are both under -50 dB, and the relative error of _P1 is smaller than the result of _P2 .

Table 1 shows the CPU resources cost under different width types. As Table 1 shows, the time step of the proposed method is not limited by the CLF stability condition and can be set as the same time step at these three situations. When the width of thin dielectric layer is d=2.4×^10-5 m, the time step of the proposed method is 402 times that of the traditional FDTD method, and the total CPU time for the proposed method can be reduced to about 0.63% of the traditional FDTD method, with the accuracy still being guaranteed. And this advantage will be further enhanced as the width of the dielectric layer decreased. Considering the complexity of solving a banded nested matrix, the proposed method consumes more storage space than the traditional FDTD method, which is about 291 times the traditional FDTD method.

Table 1: Comparison of the CPU resources for 2D waveguide.

	d = 2.4 × 10 - 2 m	d = 2.4 × 10-4 m	d = 2.4 × 10-5 m
Traditional FDTD	Modified AH-FDTD	Traditional FDTD	Modified AH-FDTD	Traditional FDTD	Modified AH-FDTD
Δ t	1.6 × 10-11 s	1.6 × 10-11 s	3.98 × 10-13 s	1.6 × 10-11 s	3.98 × 10-14 s	1.6 × 10-11 s
CPU time	0.0963 s	2.5423 s	5.0036 s	2.8457 s	312.7724 s	2.9816 s
Memory	0.79 M	216.72 M	0.82 M	232.69 M	0.83 M	236.81 M

5. Conclusion

The AH-FDTD unconditionally stable method combined with high-order algorithm is proposed in this paper. Using this modified AH-FDTD method, the dispersion error induced by the space grid can be reduced, which makes it possible to set coarser grid and generally reduces the total number of grids. In simulation, two 2D parallel waveguide models are established and the simulation results show that the accuracy of the proposed method is higher than the traditional AH-FDTD when compared with the analytical solution and the efficiency of the proposed method is higher than the traditional FDTD method in analysing 2D waveguide problems with fine-structure.

References

[1] G. Sun, C. W. Trueman, "Unconditionally stable Crank-Nicolson scheme for solving two-dimensional Maxwell's equations,", Electronics Letters , vol. 39, no. 7, pp. 595-597, 2003.

[2] T. Namiki, "A New FDTD algorithm based on alternating-direction implicit method,", IEEE Transactions on Microwave Theory and Techniques , vol. 47, no. 10, pp. 2003-2007, 1999.

[3] Y.-S. Chung, T. K. Sarkar, B. H. Jung, M. Salazar-Palma, "An unconditionally stable scheme for the finite-difference time-domain method,", IEEE Transactions on Microwave Theory and Techniques , vol. 51, no. 3, pp. 697-704, 2003.

[4] Y.-T. Duan, B. Chen, Y. Yi, "Efficient implementation for the unconditionally stable 2-D WLP-FDTD method,", IEEE Microwave & Wireless Components Letters , vol. 19, no. 11, pp. 677-679, 2009.

[5] H. Lin, G. Wang, F. Liang, "A novel unconditionally stable PSTD method based on weighted laguerre polynomial expansion,", Journal of Electromagnetic Waves and Applications , vol. 23, no. 8-9, pp. 1011-1020, 2009.

[6] X.-K. Wei, W. Shao, H. Ou, B.-Z. Wang, "An efficient higher-order PML in WLP-FDTD method for time reversed wave simulation,", Journal of Computational Physics , vol. 321, pp. 1206-1216, 2016.

[7] Z.-Y. Huang, L.-H. Shi, B. Chen, Y.-H. Zhou, "A new unconditionally stable scheme for FDTD method using associated Hermite orthogonal functions,", IEEE Transactions on Antennas and Propagation , vol. 62, no. 9, pp. 4804-4809, 2014.

[8] Z.-Y. Huang, L.-H. Shi, Y.-H. Zhou, B. Chen, "An improved paralleling-in-order solving scheme for AH-FDTD method using eigenvalue transformation,", IEEE Transactions on Antennas & Propagation , vol. 63, no. 5, pp. 2135-2140, 2015.

[9] J. Fang, Time Domain Finite Difference Computation for Maxwell's Equations , University of California, Berkeley, Calif, USA, 1989.