A Semi-discrete Scheme for Computing Two-Dimensional Electromagnetic Field in Time Domain

A Semi-discrete Scheme for Computing Two-Dimensional Electromagnetic Field in Time Domain

Shyh-Kang Jeng* (1)

(1) _{Department of Electrical Engineering and Graduate Institute of Communication} Engineering, National Taiwan University, Taipei, Taiwan

skjeng@ew.ee.ntu.edu.tw Abstract

This paper applies an unconditionally stable semi-discrete (SD) scheme to compute the two-dimensional electromagnetic field in time domain. Numerical dispersion of this scheme is derived and compared with the alternate-direction-implicit (ADI) FDTD and the Crank-Nicolson (CN) FDTD methods. The dispersion curve of the proposed scheme is found to be the lower and the upper limits of those of the explicit and the implicit FDTD methods, respectively. As a numerical example, the adaptive Runge-Kutta method is adopted to solve the semi-discrete Maxwell equations for the fields in a 2D TM PEC cavity. Numerical results reveal that the SD scheme is much accurate than the ADI FDTD method. The computation speed, however, still has to be improved.

Introduction

The semi-discrete (SD) scheme has been proven unconditionally stable, and has much better numerical dispersion properties compared with implicit FDTD methods such as the alternate-direction-implicit (ADI) FDTD and the Crank-Nicolson (CN) FDTD methods [1]. It has also been implemented for a one-dimensional problem using matrix exponential functions [1]. The computation of the matrix exponential solution takes a lot of time. Thus in this paper, we not only extend the SD scheme to two-dimensional problems, but also implement the scheme by the faster adaptive Runge-Kutta method. The numerical dispersion for the 2D infinite frees-space propagation problem is derived and compared with those of two implicit FDTD methods, the alternate-direction-implicit (ADI) FDTD [2] and the Crank-Nicolson (CN) FDTD [3] methods. A 2D TM PEC cavity problem is used as a numerical example to show that the SD scheme is much more accurate than the ADI approach, though there are still rooms for improving the computation speed.

2D Semi-discrete Scheme

Assume that the 2D space is subdivided into Yee cells, the Maxwell equations for a TM problem can be approximated as

}

)

2

/

1

,

(

)

2

/

1

,

(

)

,

2

/

1

(

)

,

2

/

1

(

{

1

)

,

(

y

j

i

H

j

i

H

x

j

i

H

j

i

H

dt

j

i

dE

_z y y _{x x}

∆

−

−

+

−

∆

−

−

+

=

ε

(1a)

}

)

,

(

)

1

,

(

{

1

)

2

/

1

,

(

y

j

i

E

j

i

E

dt

j

i

dH

_{x z z}

∆

−

+

−

=

+

µ

(1b) 1-4244-0123-2/06/$20.00 ©2006 IEEE 3833

}

)

,

(

)

,

1

(

{

1

)

,

2

/

1

(

x

j

i

E

j

i

E

dt

j

i

dH

_{y z z}

∆

−

+

=

+

µ

(1c)

Similar equations for TE problems can be given through duality.

To analyze the stability, we assume that the field components in (1) are in the form ofE_x(i, j)=E_x0(t)e−j(iβx∆x+jβy∆y)_{. The equation (1) then becomes a system of ordinary}

differential equations

[ ] [ ][ ]

0 _{A0 X0}

dt X d

= with

[ ]

X

₀

₌

[

E

_z0

H

_x0

H

_y0

]

T . We can prove that all eigenvalues of

[ ]

A

₀ are imaginary, thus the solution will be always stable. The numerical dispersion can be obtained by assuming

[ ] [ ]

X ₌ ejnω∆t

0

0 X , and transforming (1) to a matrix eigenvalue equation. With the help of Mathematica, the eigenvalues are solved out to find the numerical dispersion relationship as







 ∆

+











 ∆

=













2

sin

2

sin

2 2 2 2 2

_y

x

c

y y x x

β

β

β

β

ω

₍₂₎

which is independent of the time step size ∆t. In other words, we are free to choose the time step size. Note that c=1/

εµ

in (2).

Discussions on Numerical Dispersions

Since SD and the implicit FDTD methods are all unconditionally stable, the numerical dispersion does not depend solely on the Courant number any more. The dispersion varies with the sampling rates in space and time. It is natural to express the numerical dispersion in terms of spatial and temporal sampling rates. A common measure for spatial sampling rate is the mesh density

N

_λ, number of spatial divisions per wavelength. The temporal sampling rate thus can be measured by

N

_T, number of sampling intervals in a periodT . We can show that the Courant number equals

N /

_λ

N

_T. Under the assumption that

_∆

x

₌

_∆

y

₌

_∆

, and

β

x

=

β

cos

φ

,

β

y =

β

sin

φ

,

β

=

ω

/

v

,

∆

= /

λ

λ

N

,

N

_T

=

2

π

/(

ω

∆

t

)

, equation (2) can be simplified as













+













=













π

π

_φ

π

_φ

λ λ λ

sin

sin

cos

sin

2 2 2

N

v

c

N

v

c

N

(3)

where v represents the phase velocity in the SD scheme.

Figure 1 compares the relative phase velocity of the SD scheme with those of the ADI and the CN methods by varying the angle

φ

of propagation direction for

N

λ

=

N

T

=

50

.

This diagram shows that the numerical anisotropy of the SD scheme is almost the same with that of the CN method, and is much better than that of the ADI method. The ADI and CN results in Fig. 1 agree with those in [3].

.

Figure 2 displays the relative phase velocity with respect to the mesh density for the SD scheme and the ADI method with different temporal sampling rates. We can observe that the SD curve behaves as an upper limit of the ADI curves as

N

_T increases. This is

reasonable, because when the time step size approaches zero, the ADI equations approach the semi-discrete Maxwell equations (1). This argument is supposed to hold for all implicit FDTD methods. The same reasoning can also be used to infer that the explicit FDTD dispersion curves approach the SD dispersion curve, too, except that the SD curve becomes a lower limit this time. Note also that for a given

N

_T, increasing the mesh density does not help too much in reducing errors. On the contrary, the SD error keeps decreasing as

N

_λincreases.

Numerical Example

The semi-discrete Maxwell equations (1) can be solved by many numerical methods. A numerical solution by putting the solution into a matrix exponential form has been tried in [1], and found unsatisfactory. This study adopts the adaptive Runge-Kutta (ARK) method to solve (1) plus boundary conditions. Some numerical experiments on one-dimensional problems show that the ARK implementation is much faster and more stable than the implementation using the matrix exponential solution.

A 2D square PEC cavity TM problem is chosen as a numerical example. Assume that the cavity is excited by the initial condition

0

)

0

,

,

(

)

0

,

,

(

,

1

2

1

1

2

1

)

0

,

,

(

=



_

−

−

_



_



−

−



_

H

x

y

=

H

x

y

=

a

y

a

x

y

x

E

_{z x y} (4)

where a is the dimension of the cavity and divided into N cells. The mesh density and temporal sampling density are defined with respect to the fundamental

TM

₁₁ mode. ThusNλ = 2N . Figure 3 shows the relative error defined by

2 2 exact exact

X

X

X

err

=

−

(5)

where X and

X

_exactrepresent column vectors whose components are all of the computed and the exact sampled E_z,

_η

H_x,

_η

H_y when t =3.5T, and X ₂is the L2 norm of X. Here

η

is the intrinsic impedance of free space. From Fig. 3 we can observe that the SD-ARK scheme is much accurate than the ADI FDTD method by which

N

_T has to be very large to achieve similar accuracy level. However, the AD-ARK scheme is still slow. Some faster implementations are being investigated.

References

[1] S. K. Jeng, “Time domain electromagnetic field computation with a semi-discrete scheme as an alternative to implicit FDTD methods, “ 2005 IEEE AP-S International Symposium and USNC/URSI National Radio Science Meeting, Washington DC, July 2005.

[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, October 1999.

[3] G. Sun and C. W. Trueman, “Unconditionally stable Crank-Nicolson scheme for solving two-dimensional Maxwell equations,” Electronics Letters, vol. 39, no. 7, pp. 595-597, April 2003.

Fig. 1. Relative phase velocity vs. propagation angle

Fig. 2. Relative phase velocity vs. mesh density

Fig. 3. Relative error norm of the cavity problem vs. mesh density 0 10 20 30 40 50 60 70 80 90 100 0.95 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 mesh density re la ti ve p h a s e v e lo c ity v /c SD ADI CN N_T=10 NT=20 NT=40 0 10 20 3 0 4 0 5 0 6 0 7 0 8 0 90 0 . 98 6 0 . 98 8 0. 9 9 0 . 99 2 0 . 99 4 0 . 99 6 0 . 99 8 1 p ro p a g a tio n a ng le p hi in d e g re e s re la ti ve p h a s e ve lo c it y v/ c S D A D I C N s=1 s=2 s=3 20 40 60 80 100 120 140 160 180 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 mesh density re la ti v e e rro r n o rm SD ADI NT=40 NT=60 NT=80 3836