close

Вход

Забыли?

вход по аккаунту

?

Development of an accurate FDTD technique and its applications in electromagnetics and microwave circuits

код для вставкиСкачать
INFORMATION TO USERS
This manuscript has been reproduced from the microfilm master. UMI
films the text directly from the original or copy submitted. Thus, some
thesis and dissertation copies are in typewriter face, while others may be
from any type o f computer printer.
The quality of this reproduction is dependent upon the quality of the
copy subm itted.
Broken or indistinct print, colored or poor quality
illustrations and photographs, print bleedthrough, substandard margins,
and improper alignment can adversely afreet reproduction.
In the unlikely event that the author did not send UMI a complete
manuscript and there are missing pages, these will be noted.
Also, if
unauthorized copyright material had to be removed, a note will indicate
the deletion.
Oversize materials (e.g., maps, drawings, charts) are reproduced by
sectioning the original, beginning at the upper left-hand comer and
continuing from left to right in equal sections with small overlaps. Each
original is also photographed in one exposure and is included in reduced
form at the back o f the book.
Photographs included in the original manuscript have been reproduced
xerographically in this copy.
Higher quality 6” x 9” black and white
photographic prints are available for any photographs or illustrations
appearing in this copy for an additional charge. Contact UMI directly to
order.
UMI
A Beil & Howell Information Company
300 North Zed) Road, Ann Arbor MI 48106-1346 USA
313/761-4700 800/521-0600
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
DEVELOPMENT OF AN ACCURATE FDTD TECH N IQ U E AND ITS
APPLICATIONS IN ELECTROM AGNETICS AND M ICRO W A VE
CIRCUITS
A Dissertation
Presented for the
Doctor of Philosophy-Engineering Science in E lectrical Engineering
Degree
The University of Mississippi
Adel Mohamed Abdin
December 1997
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
UMI Number: 9835633
UMI Microform 9835633
Copyright 1998, by UMI Company. All rights reserved.
This microform edition is protected against unauthorized
copying under Title 17, United States Code.
UMI
300 North Zeeb Road
Ann Arbor, MI 48103
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Copyright ®Adel Mohamed Abdin. 1998
All rights reserved
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
To the Graduate Council:
I am submitting herewith a dissertation written by ADEL ABDIN entitled
“DEVELOPMENT OF AN ACCURATE FDTD TECHNIQUE AND ITS APPLICATIONS
IN ELECTROMAGNETICS AND MICROWAVE CIRCUITS ”. I have examined the final
copy of this dissertation for form and content and recommend that it be accepted in partial
fulfilment of the requirements for the degree of Doctor ofPhilosophy in Engineering Science,
Professor-of Electrical Engineering
(Director of Dissertation)
Chair and Professor of Electrical
Engineering
(Co-Director o f Dissertation)
We have read this dissertation and
recommend its acceptance:
Professor of Electrical Engineering
Professor of Mathematics
Accepted for the Council
Dean of the Graduate School
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ACKNOWLEDGMENTS
The author wishes to express his sincere gratitude and thanks to Allah the
creator. The author wishes to express his deep gratitude to his dissertation directors,
Drs. Atef Z. Elsherbeni and Charles E. Smith, for their fruitful suggestions, guidance,
patienceassistance, and encouragement throughout the course of this work.
Sincere thanks to Dr. Kishk for helpful discussion and cooperation.
Special thanks to Drs. Atef Z. Elsherbeni, Charles E. Smith, Ahmed A. Kishk,
and Eldon L. Miller for serving as members of the dissertation committee.
The author would also like to acknowledge the Egyptian Government for
providing the financial support during the course of this study and to Dr. A. Elsohly
for helping the author to come to the University of Mississippi.
Finally, the author is wholly grateful to his parents, his wife, and his kids, for
their support, encouragement, and patience throughout his graduate study.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
iv
ABSTRACT
DEVELOPMENT OF AN ACCURATE FDTD TECHNIQUE AND ITS
APPLICATIONS IN ELECTROMAGNETICS AND MICROWAVE CIRCUITS
ABDIN, ADEL. B.S., Military Technical College, Egypt, 1980. Diploma and M.S., Ain
Shams University, Egypt, 1987 and 1990 respectively. Ph. D. University of Mississippi,
1998.
Dissertation directed by Drs. Atef Z. Elsherbeni and Charless E. Smith.
The present research investigates the properties of a fourth order accuracy, finite
difference time domain technique with special attention to the electromagnetic and microwave
applications of this technique.
The total field formulations of the fourth order scheme (4x4) for electromagnetic and
microwave applications has been derived in one, two and three dimensional space in lossy and
lossless media. The fundamental performances o f this algorithm, which include accuracy,
memory requirement and dispersion are presented and compared with the Yee’s second order
accurate (2x2) scheme through several numerical experiments for different dimensional
spaces. In a three dimensional space, the fourth order FDTD technique is used to compute
the resonant frequencies of a rectangular cavity resonator with different discretization in
space, and the results are compared with that o f the second order algorithm and the exact
solution. The propagation of a plane wave with a Gaussian pulse waveform is also
investigated in a free space media and with the existence of a lossy or lossless dielectric slab
or magnetic slab. The numerical results in this case for the (2x2) and the (4x4) order FDTD
methods are analyzed and compared with the second order in time and fourth order in space
(2x4) technique. The modeling of planar microstrip circuits such as microstrip transmission
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
V
lines and low-pass filters is investigated with these FDTD schemes. The computed reflection
and transmission parameters in a filter example are compared with the available experimental
data, and the results show good agreement. The soft voltage source (voltage source with
internal resistor) and the matched load are modeled for the fourth order FDTD technique to
excite and absorb the traveling wave that reaches the load or the source terminals.
Comparisons are made between the (4x4) and the (2x2) order FDTD methods in a
two dimensional free space. Two different excitation sources are applied to test the accuracy
and level o f dispersion of the fourth order method with respect to the well known second
order Yee algorithm at different frequency bands. Comparison with the exact solution is also
provided.
In a one dimensional space, fields are computed using the (2x2) and (4x4) order
FDTD schemes in free space and in a dielectric slab due to an incident plane wave. The
numerical analysis shows that the fourth order scheme has higher accuracy and reduces the
required computer resources.
The problem due to the large memory size required in the implementation of the FDTD
method is also partially solved using the proposed higher order technique presented in this
research project.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TABLE OF CONTENTS
Page
LIST OF F IG U R E S .....................................................................................................
viiii
Chapter
1.
INTRODUCTION..........................................................................................................1
2. TOTAL FIELD ANALYSIS IN THREE DIMENSIONAL SPACE USING
HIGHER
ORDERS FDTD TECHNIQUES.............................................................. 6
2.1 Total Field Difference Equations of Fourth Order FDTD Scheme in Time and
Space (4x4) ........................................................................................................... 6
2.2 Total Field Difference Equations of Second Order in Time and Fourth order in
Space FDTD Scheme (2 x 4 )............................................................................... 25
2.3 Total Field Difference Equations of Second Order FDTD Scheme in Time and
Space (2x2) ......................................................................................................... 29
2.4 Node Distribution of Different Orders FDTD Schemes ...................................31
2.5 Applications..................................................................................................... 31
2.5.1 Cavity Resonator ....................................................................................35
2.5.2 Free Space .............................................................................................. 40
2.5.3 Lossy and Lossless Dielectric Structures ............................................... 49
2.5.4 Lossy and Lossless Magnetic Structures ............................................... 55
2.5.5 Planar Microstrip Circuits ....................................................................... 60
2.5.5.1 Modeling of a Resistive Voltage Source and a Resistive Load 60
2.5.5.1 Microstrip Transmission Lines ................................................... 67
2.5.5.2 Microstrip Low-Pass F ilte r.......................................................... 74
3.
TOTAL FIELD ANALYSIS IN TWO DIMENSIONAL SPACE USING HIGHER
ORDER FDTD TECHNIQUES ............................................................................... 78
3.1 Total Field Formulation of Fourth Order FDTD Scheme ...................................78
3.2 Total Field Difference Equations of Second Order FDTD Scheme ................. 82
3.3 Node Distribution of the Second and Fourth Order FDTD S chem es............... 82
3.4 Dispersion and Stability Analysis..................................................................... 83
3.4.1 Numerical Dispersion ..................................................................................88
3.4.2 Stability Analysis ........................................................................................ 94
3.5 Numerical Verification.......................................................................................... 95
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
vii
Chapter
4.
TOTAL FIELD ANALYSIS IN ONE DIMENSIONAL SPACE USING HIGHER
ORDER FDTD TECHNIQUES............................................................................. 113
4.1
4.2
4.3
4.4
5.
Page
Total Field Formulation of Fourth Order FDTD Schem e...................................113
Total Field Difference Equations of Second Order FDTD Scheme ..................115
Node Distribution for the Second and Fourth Order FDTD Schemes ........... 115
Applications..................................................................................................... 117
4.4.1 Free Space .............................................................................................. 119
4.4.2 Lossy and Lossless Dielectric Structures............................................... 119
CONCLUSIONS ..................................................................................................... 128
REFERENCES............................................................................................................... 130
APPENDIX...................................................................................................................... 135
A. The Finite Difference Approximations for Second and Third Order
Derivatives..................................................................................................... 135
V ita ............................................................................................................................................... 138
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
viii
LIST OF FIGURES
Fig.
Page
2.1
Spatial depositions of the electromagnetic nodes involved in the
computation o f electric field components, (a) Ex with (2x2), (b) Ex
with (4x4), (c) E y with (2x2), (d) E y with (4x4), (e) E z with (2x2)
(f) E z with (4x4)............................................................................................................ 32
2.2
Spatial depositions o f the electromagnetic nodes involved in the
computation o f magnetic field components, (a) Hx with (2x2), (b)
Hx with (4x4), (c) H y with (2x2), (d) H ywith (4x4), (e) H . with
(2x2) (£) H z with (4x4)...............................................................*................................33
2.3
E z versus time calculated in a domain o f size 100x100x100 layers
with side equal 4 mm using second and fourth order FDTD
schemes when a current source with a Gaussian pulse waveform
is applied at the point (50,50,50) in the z-direction. Ez is computed
at the points: (a) (60,50,50), (b) (65,50,50), (c) (70,50,50), (d)
(75,50,50), (e) (80,50,50), and (f) (85,50,50)............................................................ 42
2.4 The shapes o f a Gaussian pulse propagated in free space using
the same parameters mentioned in Fig. 2.3. Ez versus delta-*
position after various numbers o f time steps (n) when a current
source of a Gaussian pulse wave form is applied at(50,50,50). (a)
n=20, (b) n=35, (c) n=50, (d) n=65, (e) n=80, and (f) n=95....................................... 45
2.5 E, versus time calculated in a domain of size 50x50x50 layers with
side equal 4mm using (2x2), (2x4), and (4x4) order FDTD schemes
when a current source with a Gaussian pulse wave form is applied
at the point (25,25,25) in the z-direction. E. is computed at the
points: (a) (35,25,25), (b) (40,25,25)............ ‘ ............................................................48
2.6 E . versus time calculated in a domain of size 50x50x50 layers with
side equal 4mm using (2x2), (2x4), and (4x4) order FDTD schemes
when a Gaussian pulse as a current source is applied at the point
(25,25,25) in the z-direction when At = 5.929 ps. E. is computed
at the points:(a) (35,25,25), (b) (40,25,25)............. ‘ ............................................... 50
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ix
Fig.
Page
2.7 E 2 versus time calculated in a domain o f size 50x50x50 layers with
side equal 4mm using (2x2), (2x4), and (4x4) order FDTD schemes
when a Gaussian pulse as a current source is applied at the point
(25,25,25) in the z-direction when zlt= 9.317 ps. Ez is computed
at the points: (a) (35,25,25), (b) (40,25,25)............................................................... 51
2.8 E, at different values o f er using (2x2) FDTD order scheme: (a) E,
versus delta-x position, (b) ^ v ersu s time.................................................................... 52
2.9 E, at different values o f er using (4x4) FDTD order scheme: (a) E.
versus delta-x position, (b) Esversus time.................................................................... 53
2.10 E . as function of delta-x position after 30 time steps with different
O' using: (a) (2x2) order scheme, (b) (4x4) order scheme.............................................54
2.11 A comparison between (2x2) and (4x4) is performed when a current
source with a Gaussian pulse waveform is applied at (40,40,40)
in a domain of size 80x80x80 when a lossless dielectric slab of
er=2.0 is located between delta-x positions 50 and 60 and both
delta-y positions and delta-z positions 20and 60. (a) Ez versus time at
(45.40.40). (b) E, versus time at (55,40,40). (c) E, versus time at
(65.40.40).................................................................................................................... 56
2.12 A comparison between (2x2) and (4x4)is performed when a current
source of a Gaussian pulse wave form is applied at (25,25,25) in
a domain of size 50x50x50 when a lossy dielectric slab of o= 0.1 S/m
and €r =2.0 is located between delta-x positions 30 and 40 and both
delta-y delta-z positions 10 and 40. (a) E. versus time at (30,25,25).
(b) E. versus time at (35,25,25). (C) E . versus time at (40,25,25).......................... 57
2.13 A comparison between (2x2), (2x4), and (4x4) using a current source
with a Gaussian pulse waveform in the z-direction at the point
(50,50,50) in a domain of size 100x100x100 cells with a lossless
dielectric slab (£,=2.0) located between delta-x positions 70
and 80 and both delta-j^ and delta-z positions at 20 and 80. E , is plotted
versus time at different test points: (a) at (60,50,50), (b) at (65,50,50)
(c) at (70,50,50). (d) at (75,50,50). (e) at (80,50,50), (f) at (85,50,50)............... 58
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
X
Fig.
Page
2.14 Hz versus time at different values of /fusing (2x2) order scheme, (a)
at
(25.25.25). (b) at (35,25,25). (c) at (45,25,25)........................................................ 61
2.15 H. versus time at different values of p r using (4x4) order scheme.(a)
at
(2*5,25,25). (b) at (35,25,25). (c) at (45,25,25)........................................................ 62
2.16 H. versus time when p r =2 using (2x2) and (4x4) order schemes, (a) at
(25.25.25). (b) at (35,25,25). (c) at (45,25,25)........................................................ 63
2.17 H. versus time at (45,25,25). (a) using (2x2) order scheme, (b) using
(4x4) order scheme, (c) both (2x2) and (4x4) order schemes when
am=0. (d) both (2x2) and (4x4) order schemes when om= 5xl03....................... 64
2.18 Geometry o f a microstrip line (all units are in mm and the cross-sections
are not drawn to scale)................................................................................................. 69
2.19 Voltage versus time at different test points on the transmission line
computed using (2x2) and (4x4) order schemes where the source is
a voltage source of a Gaussian pulse waveform. The test points are: (a)
(7.29.15). (b) (7,29,45). (c) (7,29,75)....................................................................... 70
2.20 Voltage versus time at different test points on the transmission line
computed using (2x2) and (4x4) order schemes when the source is
a current source of a Gaussian pulse waveform. The test points are: (a)
(7.29.15). (b) (7,29,45). (c) (7,29,75)....................................................................... 72
2.21 Voltage versus time at different test points on the transmission line
computed using (2x2) and (4x4) order schemes when the source is
a voltage source with the derivative of a Gaussian pulse waveform. The
test points are: (a) (7,29,15). (b) (7,29,45). (c) (7,29,75)...................................... 73
2.22 Microstrip low-pass filter configuration (cross-sections are not drawn to
scale)............................................................................................................................. 75
2.23 The experimental and numerical values of Su and
in dB for the lowpass filter, (a) reflection coefficient, (b) transmission coefficient.............................. 77
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
xi
Fig.
Page
3.1a Spatial depositions of the electromagnetic nodes involved in the
computation of E /iJ ) component in two-dimensional space for
the (2x2) scheme...........................................................................................................84
3. lb Spatial depositions of the electromagnetic nodes involved in the
computation of Ez (ij) component in two-dimensional space for
the (4x4) scheme...........................................................................................................84
3.2a FDTD domain and boundary in the TM2 polarization case for the
computation of E / i j ) ................................................................................................... 85
3.2b FDTD domain and boundary in the TM2 polarization case for the
computation of H /iJ ) ...................................................................................................86
3.2c FDTD domain and boundary in the TM2 polarization case for the
computation of H / i j ) ...................................................................................................87
3.3 Normalized dispersion curves for both (2x2) and (4x4) order
FDTD schemes for a plane wave traveling at 0° with respect to
the grid lines versus PPW discretization......................................................................92
3.4 Phase error curves for the YA and (4x4) methods fo ra plane
wave traveling at 0°with respect to the grid lines versus PPW
discretization................................................................................................................. 92
3.5 Variation of the normalized phase velocity with wave propagation
angle in a two dimensional FD-TD grid using YA and (4x4)
resolutions..................................................................................................................... 93
3.6 Variation of the normalized phase velocity error with wave
propagation angle in a two dimensional FD-TD grid using YA and
(4x4) resolutions........................................................................................................... 93
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
xii
Fig.
Page
3.7 Ex versus time calculated in a domain of size 1000x1000 cells
with a side equal 4 cm using second and fourth order FDTD
schemes when a current source with a Gaussian pulse waveform is
applied at the point (500,500) in the z-direction. The shapes o f
a Gaussian pulse propagated in free space using the same parameters
mentioned in Fig. 3.7 ,E X versus delta-x position after various
numbers o f time steps, (a) n=50 (4.7 ns). (b) n=200 (18.8 ns). (c)
n=400 (37.6 ns). (d) n=600 (56.4 ns).(e) n=800 (75.2 ns), (f)n=l 100
(103.4 ns)....................................................................................................................... 99
3 .8
The shapes o f a Gaussian pulse propagated in free space using the
same parameters mentioned in Fig. 3.7. Ex versus delta-x position
after various numbers of time steps: (a) n=50 (4.7 ns). (b) n=200
(18.8 ns). (c) n=400 (37.6 ns). (d) n=600 (56.4 ns). (e) n=800
(75.2 ns).(f) n=l 100 (103.4 ns)...................................................................................100
3.9
The current source with a Gaussian a pulse waveform and with a
Gaussian pulse modulated with a sinusoidal waveform in the
frequency domain.........................................................................................................104
3.10 E. versus time calculated in a domain of size 211x211 cells with
a side equal 2.5 cm using second and fourth order FDTD
schemes when the source is applied at the point (106,106) in the
z-direction. E. is computed at the point (86,46) using (a) a current
source with a Gaussian pulse waveform, (b) a current source with
a Gaussian pulse modulatedwith a sinusoidal waveform.........................................104
3.11 The current source with a Gaussian a pulse waveform and with a
Gaussian pulse modulated with a sinusoidal waveform in the frequency
domain.......................................................................................................................... 107
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
X lll
Fig.
Page
3.12 E. versus time calculated in a domain of size of 1000x1000 cells with
a side equal 4cm using second and fourth order FDTD schemes when
the excitation is a current source with a Gaussian pulse waveform
applied at the point (500,500) in the z-direction. Ez is computed at the
following points: (a) (600,500). (b) (700,500). (c) (800,500). (d)
(900,500).................................................................................................................. 108
3.13 Ez versus time calculated in a domain size of 1000x1000 cells with a
side equal 4 cm using second and fourth order FDTD schemes when
the excitation is a current source with a Gaussian pulse modulated
with a sinusoidal waveform applied at the point (500,500) in the
z-direction. Ez is computed at the following points: (a) (600,500).
(b) (700,500). (c) (800,500). (d) (900,500)........................................................... 110
4.1
Spatial depositions of the electromagnetic nodes involved in the
computation of E /i) component in one dimensional space with the
propagation in the x-direction for: (a) (2x2) order scheme, (b) (4x4)
order scheme............................................................................................................... 116
4.2 The shapes of a Gaussian pulse propagated in free space at different
time steps (a) 100 time steps ( time is 0.034 ns), (b) 500 time steps
(time is 0.169 ns), (c) 1500 time steps (time is 0.508 ns), (d) 2000
time steps ( time is 0.678 ns)...................................................................................... 120
4.3
Interaction of a Gaussian pulse plane wave with a lossy dielectric slab,
with relative permitivity 4.0, 30 cm thick, located between delta-x
positions 250 and 309, electric conductivity 0.05, and magnetic
conductivity 1.5 calculated using FDTD. Plots are electric field versus
delta-x position through different cells after different time steps. Each
time step is 8.33 ps. (a) 100 time steps (time is 0.833 ns), (b) 200 time
steps ( time is 1.666 ns), (c) 400 time steps (time is 3.332 ns), (d) 625
time steps (time is 5.21 ns)......................................................................................... 121
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
x iv
Page
Fig4.4
Interaction of a Gaussian pulse plane wave with a lossy dielectric slab,
with relative permitivity 4.0, 30 cm thick, located between delta-x
positions 250 and 309, electric conductivity 0.05, and magnetic
conductivity 1.5 calculated using FDTD. Plots are electric field versus
delta-x position through different cells after 400 time steps (time is
3.332 ns)
124
4.5
Interaction o f a Gaussian pulse plane wave with a lossless dielectric
slab, with relative permitivity 4.0, 9 cm thick, located between delta
-x positions 250 and 309 calculated using FDTD. Plots are electric
field versus delta-x position through different cells after different time
steps. Each time step is 2.5 ps, cell size is 1.5 cm, and xp is 170. (a)
50 time steps (time is 0.125 ns), (b) 300 time steps (time is 0.75 ns),
(c) 625 time steps ( time is 1.56 ns)......................................................................... 125
4.6
Interaction o f a Gaussian pulse plane wave with a lossless dielectric
slab, with relative permitivity 4.0, 9 cm thick, located between
delta-x positions 125 and 154 calculated using FDTD. Plots are
electric field versus delta-x position through different cells after
different time steps. Each time step is 5 ps, cell size is 3 cm, and xp
is 85. (a) 50 time steps (time is 0.25 ns), (b) 175time step (time is
0.875 ns), (c) 300 time steps (time is 1.5 ns). (p = 0.5)........................................ 126
4.7
Interaction o f a Gaussian pulse plane wave with a lossless dielectric
slab, with relative permitivity 4.0, 9 cm thick, located between
delta-x positions 62 and 77 calculated using FDTD. Plots are electric
field versus delta-x position through different cells after different
time steps. Each time step is 10 ps, cell size is 6 cm, and xp is 42.5.
(a) 12 time steps (time is 0.12 ns), (b) 75 time step (time is 0.75 ns),
(c) 150 time steps (time is 1.5 ns)...........................................................................127
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1
Chapter 1
INTRODUCTION
The first paper on the Finite-Difierence Time-Domain technique (FDTD), was by Yee
in 1966 [1], This technique is a simple and elegant way to numerically solve Maxwell’s
equations in the time domain. Yee used an electric-field (E) grid, which was offset both
spatially and temporally from a magnetic-field (H) grid, to obtain update equations that yield
the present fields throughout the computational domain, in terms of the past fields. The
update equations are used in a leap-frog scheme, to incremently march the £ and H fields
forward in time. The computational domain must be large enough to have no reflection on
the region of interest. The reflected field values are set to zero at the outer boundaries. This
technique has been widely used for many electromagnetic, coupling, and propagation
problems. Its recent implementations include a wide range of applications for the analysis of
microstrip circuits, microstrip antennas, waveguides, electromagnetic biological effects, etc.
This is because the FDTD technique has many advantages. In this formulation, Maxwell’s
equations can be discretized without any analytic processes. The discretized equations are
solved in a sequential manner and are very appropriate for computer operations. This method
can be applied to problems with complex structures which can be very difficult to solve with
either analytical or other numerical methods. The FDTD approach gives a direct solution for
transient problems, and it is easy to transform the solution from the time domain into the
frequency domain. Thus, the characteristics of devices can be obtained over a wide frequency
range with one computation. The FDTD technique can be applied to lossy, inhomogeneous,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2
nonlinear, anisotropic, time varying, and dispersive media. On the other hand, most other
methods do not have this generality. However, this powerful tool has some drawbacks that
include: (1) the spurious reflections o f the outgoing wave on the outer boundary o f the
computational lattice, (2) the numerical dispersion due to the error over the frequencies in the
bandwidth of the waves propagating in the computational domain, and (3) the numerical
dispersion due to the anisotropy of the mesh which induces local deviations of the wave
vector and deformation of the wave front. Numerical dispersion and grid-an isotropy errors
can be kept small by having a sufficient number of grid spaces per wavelength. Even so, the
main disadvantage of the method is its requirement of large computer memory space and
long computational time for large scale problems. This can be overcome by the development
of more powerful the computers and the improvement of the method. According to Nyquist
sampling, the discretization in the spatial domain depends on the highest frequency related to
the analysis and the highest permitivity in the computational space. In general, the cell size
should be kept in the order o f A/20 (A is the minimum wavelength in the media) and should
fit the geometry in any direction. The time step is determined by the cell size and the velocity
of propagation in media according to Courant condition [2],
The computational domain is restricted to a limited region due to the limitation o f
computer memory. A good absorbing boundary condition (ABC) used to prevent reflections
can not only improve the overall accuracy of numerical results but also decrease the memory
requirements for the same structure. The application of ABC on the boundary is a special
formulation for updating the tangential field components on the outer boundary. Due to the
importance of ABC for the development of FDTD techniques, Mur presented the first and
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3
second order ABC for uniformly meshed homogeneous media based on the wave equation
in 1981 [3], The first order Mur’s ABC accounts for the normal incident wave and the second
order Mur’s ABC accounts for the oblique incident wave. Mur’s first order ABC has been
proved to be adequate in many applications related to microstrip circuits and antennas.
Recently, many kinds o f absorbing boundary conditions have been reported for uniformly
meshed homogeneous dielectrics. A soft voltage source (voltage source with internal resistor)
and a lumped element are used to excite and absorb the traveling wave that reaches the load
terminal and source [4,5], There were many publications applying the FDTD method on
electromagnetic scattering and propagation problems during the late seventies until the early
nineteen nineties [6-28].
The time domain finite difference technique is used in the computation of Maxwell’s
equations with fourth order accuracy in three dimensional space in 1989 by Fang [29]. This
work is done for the total field formulation in lossless medium. A fourth-order FDTD scheme
in time and space in one dimensional space is presented for radio-wave propagation in a
lossless cold plasma by Young in 1996 [30], The purpose of the research presented here is
to derive and study the fourth order scheme in time and space for lossy media in one, two and
three-dimensional space. This higher order FDTD scheme is compared with the second order
scheme in different dimensional spaces. The study will critically investigate the accuracy,
memory size requirements, and numerical dispersion o f the fourth and second order schemes.
The fourth order scheme will be used for the analysis of several applications such as the
reflection and transmission from planar interfaces, radiation problems, and planar microstrip
circuits.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4
In Chapter n , the total field FDTD formulations of the fourth order scheme are
derived in detail in a three dimensional space of lossy and lossless media. The corresponding
formulations of the second order scheme (2x2) in space and in time and the second order in
time and fourth order in space (2x4) are deduced from (4x4) formulations. The number of the
electromagnetic nodes involved in the computation for (2x2) and (4x4) formulations are
discussed for both the electric field and magnetic field components. A general FDTD
computer code is developed based on the time-dependent Maxwell’s equations to model and
solve three-dimensional problems. This code is used in many applications such as to compute
the resonant frequencies for a rectangular cavity resonator, to determine the field in the near
zone due to a current source of a Gaussian pulse waveform, to calculate the field inside and
outside a lossy or lossless dielectric or magnetic structures, to test the voltage at different
points on a microstrip transmission line and to measure numerically the transmitting and
reflection parameters of a low pass filter. Stability and accuracy constraints are also discussed.
Comparisons between (2x2) and (4x4) formulations are done in all the cases and between
(2x2), (2x4), and (4x4) in some cases. The expressions for the soft sources which are used
to excite the microstrip circuits are modified to the fourth order accuracy.
In Chapter III, the total field FDTD formulations of the fourth and second order
schemes are deduced
from the two dimensional case. The required number of the
electromagnetic nodes for the higher and lower order algorithms for (2x2) and (4x4) are
studied. A FDTD computer code for the two-dimensional problems is presented. This code
is used to obtain the responses at several test points when using two different current sources
with various waveforms, such as a Gaussian pulse and a Gaussian pulse modulated by a
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
5
sinusoidal wave. This is performed for both low frequency and high frequency sources.
Comparisons between the higher order and the Yee algorithm (YA) are presented. The
stability and accuracy constraints are discussed. The numerical dispersion is also investigated
and discussed for the two order methods.
In Chapter IV, the total field formulations are derived in one dimensional space for
the fourth and second order FDTD techniques in different media. The characteristics o f the
two algorithms are tested in free space and in a free space medium of containing a lossy or
lossless dielectric by the developed FDTD computer code.
Finally, conclusions and suggestions for future research topics are included in Chapter
V.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
6
Chapter 2
TOTAL FIELD ANALYSIS IN T H REE DIMENSIONAL SPACE
USING HIGHER ORDER FDTD TECHNIQUES
In this chapter, three finite difference algorithms o f different accuracies are discussed.
The formulation o f fourth order accurate finite difference algorithm in space and time (4x4)
is derived in lossy and lossless media. The corresponding formulation of second order in time
and fourth order in space (2x4) and second order in time and space (2x2) is derived from the
fourth order formulation. The characteristics of these different formulations for the various
order schemes are analyzed and compared through different applications such as computing
the resonant frequencies of a rectangular cavity resonator, testing the field due to a near zone
current source with a Gaussian pulse waveform in free space, a lossy or lossless dielectric
slab or magnetic slab, and modeling and measuring the parameters of a microstrip
transmission line and a low-pass filter.
2.1 Total Field Difference Equations of Fourth O rder FDTD Scheme in Time and
Space(4x4)
In a homogeneous isotropic media, the total field satisfies the following Maxwell's
equations:
BE
dt
(2-1)
£
£
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
7
= - — V * .£ - — H
ii
u
dt
(2-2)
If the electric field component is expanded in a Taylor’s series about the time t* to the
time tj + At, keeping the space point fixed at (i, j, k), [29]; then:
t*
Er (tj
k) + AA . dE
^ |v, ( t } k }
*< )m|( i j k ) - fM
+ At
E \t"(
d3-^ |
d t j
3 !
K u . j . k )
+
4
d f 4
,
A^ t2^ d2E |r>
, ( Wt k )
+ 0 ( A t5)
I
U ,A > .j.k)
(2-3)
The last term is the remainder, or error term.
Similarly, a Taylor’s series expansion at the instant (^ - A t), again keeping the space
point fixed at (i, j, k) is:
f
A vi
E
A . dE |
ft,
( » . ; . * )
A t3
3!
a3£ .
dt3 ' ( ''y’ * } +
A t2
+
a /4
4!
~
2
\ ~
d2E ,
$ t 2
'
a4£ ,
a / 4 ', ( w ' fc) +
Now, subtracting (2-4) from ( 2-3 ), we obtain:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
{ , ' J' k )
^2_4^
If one considers the interval (t* - At /2, t* + A t /2 ) instead o f the interval (t* - A t, t; + A t ),
equation (2-5 ) can be written as follows:
>U*.
BE ,
a,
1 , At .3
* 3 ( 2 )
B?E ,
r v s .*5\
a , 3 U u , . * ) + °<A' >
(2-6)
Equation ( 2-6 ) can be written more explicitly as:
E
= E * 2 + At
dt
+ - (-^ -)3
3
2
+ 0 ( A t5)
(2' 7)
Similarly, the magnetic field can be expanded as:
a
^ 1“ -s ♦ a ,
* i
a/
3
(_ 4 L )
2
. 0(A ;!)
3
a/ 3
Substituting (2-1) into (2-7), one obtains:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2. 8)
9
1
1
n*—
n -—
2
IT*
Ex OjJc) + ”
A c (tfJfc)
1
dz
1
a3# z (V*)
a3# y
(v * )
i N
/!* —
a#.
but
dz
- a.
2
0 (Ar5)
(2-9)
3 rs
1
n+—
2
H
2
a2£ .
8z 3fx
3.y 3 r
n .i
a
E
2
) - « (v^)
1
/i*—
2
1
n*—
'
(V * )
n**-—
2
A /3
24 8
dH,y
dH.z 0jX)
3y
A/
rrn
2
2
-
d %
Az
(v#
3 Az
(4 ^ -)
1
n»—
2
0 (A z5)
dz:
(2 - 10)
1
77♦—
a/f.
a^
2
H
(V-^0
3 A^
Ay
(v*)
0(A>/5)
<3yJ
(2- 11)
The second derivative in time in equation (2-9) is converted to a second derivative in space
to reduce the memory requirement for computation. This can be implemented by taking the
curl o f (2-2) and substituting (2-1) into the resulting equation. One then can get
d2E''n = J _
ye
df
d2E
ax-
d2E
dy2
2r
d2
E
8:2
t,n
g e °m
8t
£ ■ i.n
y e
(2- 12)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
10
and for the magnetic field, similarly:
&
H
1
, # H ’’n
# H
t'n
P H 1' 1 ,
a e
°m ,
= —
( -------- + -------+ -------) - (, —
+ —
)
dx2
dy2
d z2
e
\i
1*
dt2
M
' ' n
a e °m
-
dt
f j
cl
r.n
p e
(2-13)
now, substituting (2-10), (2-11), and (2-13) into (2-9), one obtains:
n
r» n * l
^
j-»n
~H
XOjJc) - Z'x QjJc)
2 4 e r2 p r
)
PAH
1
1
n -—
1
2
(c Ar)3
2,
$H ,z (V*)
a3# z
(v JO
ay-
8y 8x:
2
ay 022
y
2 i )
n -—
n-—
d2K y
H
2 i
1
1
n +—
2
l
n-—
n~—
z ( V - j.* )
n»—
11
1
1
2
O j J c)
02 d x 2
2
a3#.y (V^)
a^ ay
02 3
1
rt-r —
24 e.
2
H p, A r3 a3# > (V^)
24 e.
02'
ay-
o Ar
*
/r
1
n-—
2
a3# >- (V*)
"4
2
(v#
1
n-—
2
A r3 02£ ,'* (V^)
24 e
dt1
(2-14)
,
1
ere
11 ' 7
1
c Ar
/
p3=- A 7 ’
c Ar
"aT
a-1 5 )
I f (2-12) is substituted into (2-14) and the finite difference expressions, as discussed in
appendix I [31], are used, the final electric field updating equations are derived as:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
11
1
2
1
2
rtf—
r tf —
- n*l
E n
Sx (yjfc)
PA H
xM)
- Ps ( ^
V
1 1
2
2
2
-
2
2
z (/
1
2 (»
I t
1,
2 2
2
>
2
L
*1 Py Pr
24 n iM ?
2 #
<« ; w , ,
2
2
2
rtf—
2 tf"
2
2
2
2
1
2
rtf—
H
2 4 “H1£ r
j
!
2
2
, 1
3 . 1 ,
2
-
3 H
2
2
, 1
n Py Pz
24
1
2
2
2
2
1
2
1,
1.
+ 3 H
>
\
2
-
3 H
- H
,
, 1
1 . 1 ,
r ( i f —j f — J c f — )
-j
i
J 23
n* —
jr
2
it
1
1 . 1 ,
z ( i f —j - —Jcf—)
2 2
23
2
i
n ~—
rtf—
- 2
H
2
rtf 1
2
2
2
■ rtfl
z ( t f l j - Lj cf l )
2
2
2
z ( I * * ,- ! * - ! )
2
2
,
2>
i
,
3
2
2 2 2
1
2
( / f , * 1. 3*
,
n»—
2
2
,
2
n»—
3 . 1 ,
2
1
1
2
rtf —
rtf—
1
1
2
24 Til M r
2
)
z ( i - i 1/--I ><:»I)
1
2
rtf—
n py
~ H
rtf—
2
rtf —
2
1
1
2
rtf—
*P y
"*{
I I I
-•**-)
2
2
2
2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2
,
1 . 3 .
I,
2 2
2
z ( t f — j - —Jc f—)
12
1
nf—
2
< * 3
1 . 1.
>- ('V V * 2>
P z Px
24 W / r
- 2 H
«*-1
2
1
/i*—
2
*1 Pz P,
1
2
13
n-f—
( H
24 Tli^A
2
1
2
..
1 . 3 ,
*1 Pz
24
*1 Pz
24 ! W r
2
n. I
n»—
iff—
(■ H
1
2
2 tf
1.
"j j t - t f
y
2 2 2
y (j + L j + ± J c - L )
1
nf —
1
2
2
H
/ I
1 t 3.
yO*-jf-Jcfj)
2
1
2
2
1
nf —
3 H
2
ri p
a
a
1 O' ( _ i + _ hl )
24 T)iEr
e
li
2
2
3 H
2
1
n* —
2
2
2
1
nf —
2,
.
.
y .fl j f^
+ 3
2
H
/
J' I
1
n*—
1t I*
2
2
2
1
2
nf —
( H,
2
2
2
- H
1 1
y o*-j
2
n*i
nf —
H
2t t j )
2 2 2
Q t — J ~ — Jc ~ — 'X
2
2
2
2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
//
13
-(H
n—1
2
1
n- —
- H
z (/■
42> '42 ’* 42 )
2
z (/
42 4 *2 4 >2
1
—
2
ti p,
o
a
•—
( _1 + _?L)
24 rj1er e
ji
n --
2
2
2
'
2
2
2
n -i
- Cflr" 2!
!
! - H
2
2
2
2!
!
j)
2
2
2’
r
1
n-*-—
( c At)2 T]OeOmpy
24ri1e*nr
^
2i i i
«< 444)
1
( c AQ2 r|ggqwp2
2
2 2 2
1
n -r—
H
24ri1e ^ r
1
n-f—
ae At
n +—
\
1 1 -H
2
> -< 4 4 4 >
> < 24 42 42 )
24 e rij
-*n-l
^ (v*)
(2-16)
O A/
where
’’li = 1 +
r)2 = I +
e
(J A/
+ ——
24e
12 e
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-In\
( " 2)
(2-18)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
15
1
n-*-—
..i
3 77
77
\
1
n
-
77
1
2
n p,
Py
2 4 T i^ e J
(#
1
/ ■L
I
1
3,
2
2
2
1, - 2 7 7
-
3 77
2
t
J
-
n
2
x
2 77* 2 !
2
1
fl*—
2
77
/I
2
2
1 «.
)
In '
n4
*o-y>,4
2 2 jfc
24 )
-
, , , >
4
1
—
- 2 H
\
I
«♦—
x
x
- H
1
1
n-r —
-
2
2
77
/ I 1 i.
zC’l ' V ’l)
n p,
0.
om
1 - ( _ i + _2L )
2 4 rijE,.
e
1
tff —
2
2 2 2
- <
1
2
2
- <■•? V i >
1
n»—
/■1 ■1 t. 3.
z 0"T^’T-^r)
2
2
2
- H
1
«■*•—
3 77
"4
( ,-V 1 * 4 ,
n '- —
2
x (i
2
2 4 T i,n rE*
- 77
2
1
2
1
n» —
^ Px P.-
2
2
n *—
1
rt+—
2
< » ,. 1 ! ,! ,
2"'t 2 I
3jt
1
2
3 H
n »—
24
2
i ( ,- 1 ,4 ^ 4 , ■
n*—
2
H
+ 3 H
3 77
2 4 r i 1e r
*1 Pz
.
i
/!♦—
p,
.
* o<4;-4*4)
1 “- f ' - H '
2 4 Tiltxre r
,.i
4 2 *
2 2*
}
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
i
(2-19)
E.
n-l
cij*)
1
Jh. F n
Til
+
pAH
2
n Px
24 t i,e r
1
2
2
1
/i-—
2
H
,( ■ 4 4 * 1 )
1
n*—
3 tf
2
y
.
.
-)
+ 3 H
- H
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2
17
3
H
y
24
" \
y 4
h
< 24 42 4 2 }
1
«f—
2
( * f. 1 3> U - 2
*1 Px Py
.
( ^
y
..
1
3 ,
1.
\
y o--
2
H
1
2
( *
2
1
2
2
.
K
J
2
2
-
24
- H
y t 'v W
rtfi
2
)
/■I 1 l | J
I 'W
2
/■I U 1,
2 H
2
2
2 h ] \
3,
1
rtf—
2
1 ,-3
O'-J
2
2
2
2
2
Tj_
2
*1 P , Px
,«V? >
4
2
* 4 24 42 >2
t _1
rt f —
,
- H
1
nf—
I . 1 . 3 ,
24 T]l(ircJ
( #
y
rtf—
n Px Pr
+ 3 H
4 4 >
1
n->-—
2 H
.
1
n*—
24 TljMr
2
1
rr*—
1
n*—
^_1
4
*1 Px
y
y - H ^ y
2
2
)
2
2
2
1
rtf—
- 2
H
1
nf—
r
2
3 l . l.
X° '2 J~2 2
2 H
\
2
y
y
2
2
~
H
2
2
2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2
2
2
2
18
1
2
/if —
*1 Py Pr
x(i
24 TiiM?
1
2
1
1 . 3 .
2
2
2
-2 H
2
- H
*(/ 1 >1 U
2
2
2
2
2
2
1
2
/if —
71** —
- H
1
2
/|f —
- 2 H
2
2
2
2
2
T1 P ,
n-r —
- 3 H
H
2
1
/if—
1
n-*—
l
^Py
24 ^
2
2
2
2
*
, 1 3 . 1 ,
2
2
2
2
1
2
H
24 ’l l ^
2
- H
2
1*1,
^ Py
( ae
2
2
2
2
1
2
2
2
2
1
/|f —
n Px (
24 rll8r
e
1
2
' (/f r l
t j **—
~ H
r 2i
2
2
O
2
7 I- —
1
2
n-'—
24
e
1
2
/if —
/if—
— )
li
-
>
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
/ I
3 . 1 .
2
2
2
19
n --
-{H
n --
2x x 1 - H 2, 2 ! )
X(,-A^A^A)
XovA^vAjfc.A)
2 2 2
2 2 2
h-A
( c AQ2 Tiaegwpx
2 4 lW
„ .a
21
I 1
2
2 2
r
\
1
n-»—
( c A/)2 Tiaga mp>
H
24rli«vHr
2
2
2
1
n»—
2
/ I
1 I
2
l v 1,
2
2
- H
2
x (I
e
A*
24 e rjj
r-n i
OjJc)
(2 -20 )
Similarly, the H-field updating equations can be derived as:
1
2
n+ —
H
=
Pr ( E y (ijjc'l)
Zl
rr
H
1
n- —
2
E y {ijjc) )
PvPx
24T)?lMrer
y
( E , (jj+2jc)
3 Ez
(/>rU) + 3
E . (ljJc) - E , ( l j . ljc) )
24n5,nr
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
20
(E «
pI
3 E , f„ , i ** + 3 E z (UJc)
E z (u-uk) )
24T i^ e rfi r
Py
Pz
1
«<■—
< it-1
2
2 £ z ((V>1Jfc)
2
( E: w+ijc+1)
1
2
n+—
£z
)
24tl^nrer
-
(2 ^ 0
- 2 E ; WJt) - ^ V i ) )
1
n 1
« *-
2
Pz Px
2
24itfi|irer
( -^y (ij+ljc)
2 E y (,Vjjfc)
PzPy
E y (V. U ) )
( Ey (ij+i'X-ri)
2 E y ^ ^ i)
E
24ti ^ n rer
P_24tl^nr
p!
Ey
(ljjC+2)
^ Ey
E y o jJ c * 2)
+ 3 Ey
Ey
3 E y Q j j d ) + 3 E y (ljJc)
JJc-
E y0jJc. d
24r|^i njer
p>
( f i ♦ £= )
24ri51|ir
e
n
(
- £ ; (u. U) >
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
21
Pz
24ri l ^ r
e
n
(cAQ2a gqmP>
[ £ - (V-lJfc)
24ti$!^er
( cA/)2a eg>npz
241^^
(ij-ljc )
( £ ,> (vJt*D - £ y'
where
£, = 1 +
o At
OjJc)
o A/
1
H
+
) -
H
2
a Ar
(2-21)
^
(2-22)
2 4 [i
o „ A/
S2 = 1 + - = —
(2 -2 3 )
12p
1
/I - - —
H
and
Px (
2
e
—
H
’1
O'ljJc)
Pz P ,
2x i !
2 2 2
Ez{, j Jc) )
(y-U-1)
( * ,V
o r"
_ f '1
\
^ X{ijjc-rl)
(Ij-ljc-l) )
24TiE1Hr£r
(-^x(v-U)
2 Ez{ljJc)
Ex(lJ. lJc)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
22
P r Px
24r| 5 tn rer
( E X U -ljJ c)
+ --------
24ti^nr
2 E x (ijjc) ~ E x
( Ex
n -ijjc )
)
3 E x (ijjc,I) + 3 -Ex (ijjc) - E x Qjjc-X) )
Pr
( E x* (ijjc^z)
2
24 ti^ n r er
3 E x Qjjc-i) + 3 E z (ijjc)
_
E x (ijjc-1) )
Px
24ti^nr
( E g (, -rZjJe)
3 -^Z O-ljJc) + 3 E z VjJc)
* — T T - ( SV w i - 3
24Tl51Mrer
PxPy
24r]51nJcr
( E y o j.ijc )
( -^r
+3 ^
2 -^r
2 E y ( ,jJc)
E z (,-ljJc) )
)
-^z O-lj-ljfc) )
E y ( i J . ut) )
Px Pz2
( -^z (1 * 1 ^ 1 )
2 £ . (l<fUW
(l »!,/,£)
E
z o(/-I
- ijj-I
- ijJf c -l)
-n )
L-Z
24 t i ^ n r er
E z (,jjc+1)
2 E z {ijJc)
E Z (,jjc-1) ) ]
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
23
( — - — )[(
)-( 5 . V « - £<«*> )
W+l
24ri
e
’z
P
( c 4 0 Jo .o j)I r
~ 2
I 1 W-*’1)
( c
r*Ht l
O -ljJ c)
v -^Z
^ pfl
O -ljJc)
\
Z (<M)
'
.
1 (Vr*
A /)2 a ea mPx r c „
-
v /• r*W
Z (vJfc) /
~
™
z (.ijjc)
L z C'*1/*)
1
24
(2 -2 4 )
1
H
\
1
1
i;
Py ( Ex (V. U)
24t| ^ P r
4
Px ( -^y 0 -1 ^ )
,
E y O jJ c )
-£r (v^) )
[ E y ( ,* 2 jJ c )
Px
[
.1,
0 - 2 jJ c )
3
(/♦Wjfc)
3
+
Ey(^ i jjc)
3
^ y O j J c)
+ 3
Ey0Jtk)
E ■y
y ( ' - 1 jJc)
Ey 0 . ljjc)
24ti ^ n rer
Px Py
24 ti5 ^ rer
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
)
Px Pz
24
( Ey
' 2 Ey (ijjc)
pl
Ey
'x
-
)
2 Fxn( v > U )
Z
-
£ x" (/-
24Tl51Hr8r
-
7 F n
Z -^r (v*)
Pz2
2 Ex0jrlJc)
Ex
)
24T1
Ex OjJc*l)
2 Ex 0jJc)
Ex0jJc. 1) )
+ 3 E,’* OjJc) - E,x (,j-ljc)
- 3 E
24n
”
0 j* 2 J c )
- 'I p n
j
-^x ( v - U )
Exn(V jfc) - F"x (v-U)
24T1 ^lPrEr
.n»l
E nO
j
J c)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
One can substitute in the previous general formulations o f the total field by suitable
values of the constitutive parameters (e, ac , pi, and om) to get the total field difference
equations o f fourth order scheme (4x4) in different media (lossy and lossless dielectric, lossy
and lossless magnetic materials, and free space).
2.2 Total Field Difference Equations of Second Order in Time and Fourth Order in
Space FDTD Scheme (2x4)
By substituting (2-10) and (2-11) in the second term on right hand side of equation
(2-9) and neglecting the third term, one can get Ex in (2x4) accuracy:
1
pr n~l
OjJc)
_L
113
F n
(,jjc)
+•
^
2 2 2
13 r
H"\\
1
1
)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
( 2- 26 )
where
n3 =
aAt
- - 2-
1
(2-27)
1
1
r- n-1
and
1
" y O j J c)
~
113
T1 Pz
24 ri3£r
n Px
24 r)3er
2
y O j J c)
2
Tj
+
_
_
' liv
1
1
H
^
PA H
2 2 2
n ->—
2
3 H
n*—
2
1
n* —
x C ^ , - ^ j)
1
n*—
2
3
1 .
1.
i
n** —
x
1
n» —
,
2 2 2
"'•W
1 1
z O '-J '-J c * -)
2
n
3 //
c -f
1
n- —
2
1
1
O '- J '- 7
J c --)
2
*
1
n» —
, 1 >z. 1, - H
-- O-jJ-jS-j)
(2-28)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-29)
Similarly, the H-field equations are derived as:
1
i
n~~
H
1J-—
1JC-—
1)x= TL H * ( ' •2—
1j ' —
IIJc'—
1,)
x (i —
2
2
2
2 2 2
1
n» —
2
Pz (
^
O j J c-
24t)$3*ir
P.24tl$3nr
i)
E y O j J c)
Py ( K r (ij +IJc) - E --nO j Jc)
)
0j-2Jc) - J2 ^F’=nCm-
y 0jJc-2) - 2 F"
y Oj Jc- i)
n
+ 3 £.- (v^) -
E n
O j J c)
* (v-U)
- E •y
, O jJ c - i)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-30)
28
o Ar
fn
(2 -3 1 )
1
n»—
2
H
2 2
[ Pr ( K
2
2
2
Fn
+
J3 ^rOVJfc-l)
24ti53jx2er
Px
- K u * )
2
z (i+ljjc)
2
J
it"
x ( ijjc )
♦ 3 ^
-
r"
{tjjc-l)
- £ •;
(2 -3 2 )
24Tl53MxEr
1
H
2
—
5
z (4 ^ i^ i)
Py(^x( v- U)
24t| £3pr
24ti $3nr
i
i
n-r —
H
x ( ,211 , 4 i^ 1 i)
Tl^Pr
Pr ( ^*y (i*ljJc)
Ey (ijjc)
EX(ijJc) )
E"
*v
-'y 0~2jJc) - 33 +
■y (i*
0-1^)
(ij+2Jc) - 3 E"
- ^r " 0-1^)
+ 2 EX(ijjc) - E"r ( V - U )
(2 -3 3 )
In this case, one can substitute, in the equations of this section, suitable values of the
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
29
constitutive parameters (e, <r„ ft, and am) to get the total field difference equations of the
second order in time and fourth order in space (2x4) in different mediums (lossy and lossless
dielectric, lossy and lossless magnetic materials, and free space).
2.3 Total Field Difference Equations of Second Order FDTD Scheme in Time and
Space (2x2)
By neglecting the third term of equation (2-9) and using the central difference
approximation for the derivatives, one can get:
(e+aA/)
-?AH
(2-34)
(e+oAr)
- PA H
)
(e+oAt
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-35)
Similarly, the H„ Hy and Hz components will take the following expressions:
1
n~ —
„-i
H
1
2
2
c
H
1 1 . 1,
(P+OmA0 X(' 2 2 2
r (;
c( p +° mA0
pv ( K’z
- E.’z O jJ c )
Pr ( E y OjJc* 1) " E y (ijjc) )
(2-37)
1
n- —
n+ —
H
y
2
2
2
2
- Pr ( Z'za(i+lj
. Jc)
(p+oA/)
E z OjJc)
-
2
H
2
2
1
Py ( Ex (,-ijjc)
9 A E '}x ( i - l jJ c ) - Ez\ ^ ) )
1
2
2
C(P +OmA0
(2-38)
n -—
2
2
)
1
n*—
2
H
1
2
(p +amA0
H I I
1u ^
C(P +(JmA0
Pr ( E '}y 0 * l jJc)
Ex (ljJc) )
- E,'y
OjJc)
(2-39)
In this case, one can also substitute in the above equations the values of the constitutive
parameters (e, a0 jx, and am) to get the total field difference equations of the second order
(2x2) in different mediums (lossy and lossless dielectric, lossy and lossless magnetic
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
31
materials, and free space).
2.4 Node Distribution of Different Orders FDTD Schemes
The number of the electromagnetic nodes involved in the computation of each
component o f the field in the fourth order (4x4) and (2x4) schemes is double the number
required in the second order scheme (2x2). Figures 2.1 and 2.2 show the required number
of nodes for (2x2) and (4x4) orders FDTD schemes in three dimensional space for the electric
and magnetic field components [28].
2.5 Applications
Using FDTD, Maxwell’s equations are solved directly in time domain via finite
differences and time stepping. This is done for (2x2), (2x4), and (4x4) order schemes in free
space, dielectric, and magnetic media. A computer program based on this FDTD formulation
has been developed . This program is modified and used in several applications such as the
computation o f the resonant frequencies of a rectangular cavity resonator, investigation of the
field components in free space or in lossy or lossless dielectric slab or in lossy or lossless
magnetic slab located in a free space, when a current source with a Gaussian pulse wave form
is located at certain point in the domain. A comparison is made between (2x2) and (4x4)
formulations in all the previous applications and includes the (2x4) formulation in some cases.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
32
Z
H,
Y
H,
H.'t7
‘l
H,
(a)
(b)
|— ► Ht
H,
H,
H,
I
J L t £r
H- 1
r
t H«
H,
H,
— ► H,
(c)
(d)
*-z
I k.
▲
H,
z
H,
(e)
(0
Fig. 2.1 Spatial depositions of the electromagnetic nodes involved in the computation of
electric field components, (a) Ex with (2x2), (b) Ex with (4x4), (c) E y with (2x2),
(d) E y with (4x4), (e) E z with (2x2) (f) E 2with (4x4).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
33
Y
E, E.
E,:
(a)
(b)
E.
E,
E,
E,
_1_±^
Hu
E*
E,
E,
5.
(C)
( d)
H.
*”T
E,
Z
E,
(e)
2.2
(0
Spatial depositions of the electromagnetic nodes involved in the computation of
magnetic field components, (a) H* with (2x2), (b) H, with (4x4), (c) H v with (2x2),
( d ) H y with (4x4), (e) H* with (2x2) (f) H z with (4x4).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
34
The stability and accuracy constraints on Ax, Ay, Az, and At must be taken into
consideration. For good accuracy, Ax, Ay, and Az must be much less than the minimum
required wavelength (sA/10). For stability, the time step must be small enough so that the
field values can affect only the nearest neighbors during one time step interval. The necessary
stability condition for the three dimensional case when the discretization in space is equal in
all directions is:
rc Ar = Ax / f t
(2-38)
The excitation of the difference equations can be done by using different shapes of
pulses that depend on the frequency content. For accurate computations, the excitation source
must be absorbed at the limits of the problem space. To insure this accuracy, different
absorbing boundary conditions (ABC) are employed in this investigation, such as Mur’s and
Liao’s ABC. In other cases, a perfectly electric conductor (PEC) wall is used. In one attempt
with the fourth order scheme; the first order outer boundary condition at the comer, the
second order at the outer boundaries and fourth order at the interior region is used.
However, this procedure causes divergence to the (4x4) scheme in some applications which
have composite geometries. This problem is solved by using a small cell size in order to
represent any geometry complexity by a number of cells not less than four and applying the
forward difference approximation (FDA) at the left boundary and the backward difference
approximation (BDA) at the right boundary as an example for a one dimensional space.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
35
2.5.1 Cavity Resonator
A computer program based on FDTD has been developed to compute the resonant
frequencies (Fr) in a rectangular cavity in the time domain using (2x2) and (4x4) FDTD
algorithms. The program provides very quick information on resonant frequencies for
rectangular cavities and this allows the user to judge the feasibility of undertaking the design
of such structures. The results from the FDTD schemes are transformed from the time domain
to the frequency domain using the Fourier transform procedure and then illustrated and
compared to the results obtained from the exact solution. The results are analyzed and
conclusions are presented.
A cavity resonator of dimension 20x20x20 mm3 was chosen with PEC walls. An Ez
component with a pulse function with amplitude of 10 V/m and pulse width o f 3.85 ps is
used as a source at the center of the cavity to excite the T M ^ mode. There are three cases
discussed in this application. In the first case, the cell size to be considered is 1mm, the
number of time steps is set to 5120 and the size of each time step is chosen to be 1.9258 ps.
The domain size of the FDTD is 29x29x29 cells. The PEC walls of the cavity resonator are
chosen to be four cells in all sides to have no switching between (4x4) and (2x2) in the higher
order computation part. The frequency range and the mode numbers corresponding to the
cavity was determined from the following equation [32]:
m
m np
27t/jle ^
(— r + ( ^ )
a
b
(
c
m=1,2,3,...,5, «=1,2,3,..5, /?=0,1,2,...,5
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-41)
36
The total observed frequency range was from 10.5993 GHz to 64.9076 GHz. The results
from the FDTD higher order and the YA are transformed to the frequency domain (FD) and
the resonant frequencies are defined. These resonant frequencies are compared with the
corresponding values from equation (2-39) to determine the number o f the corresponding
modes. The results are shown in Table 2-1. The values o f the resonant frequencies computed
from FDTD algorithms ((2x2) and (4x4) )and the exact values obtained from equation (2-39)
are used into the following equation to calculate the percentage error PE,
IF r obtained using FDTD - Exact\
PE = — ------------------------------Exact
x1 0 0 %
(2-42)
The results are shown in Table 2-2. The maximum error obtained in this case when using the
YA is 1.95241, while it is 0.38304 when using the fourth order FDTD technique (4x4).
The second case that is conducted is to look at the same rectangular cavity, except
the cell size is set at 2mm. All other conditions are kept the same as in the first case except
the size of the time step was increased to 3.8516ps. This will result in a reduction in the
domain size to 19x 19x19 cells. The exact values are again compared with that obtained from
equation (2-39). The values of the resonant frequencies computed using FDTD technique
(second and fourth orders in time and space) and their PE are illustrated in Table
2 -1
and
Table 2-2, respectively. The maximum error obtained when using the YA is 1.73079 while
it is 0.85477 when using the fourth order FDTD technique.
The same program is run again with the same rectangular cavity, except that the cell
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
37
Table 2-1 Comparison o f resonant frequencies obtained from the FDTD code (second and
fourth orders) when the cell size is 1, 2, and 3.33 mm and the exact values.
Mode#
F r(GHz)
110
10.5993
112
FM3
Fr22
Fr<2
Fr21
F*i
10.5599
10.6700
10.5900
10.6899
10.5999
10.6399
18.3586
18.1399
18.3199
18.2799
18.4499
18.3400
18.4099
310
23.7001
22.2999
23.4899
23.2899
23.8099
23.5799
23.7899
311
24.8577
23.6899
24.6099
24.4499
24.9200
24.7600
24.9200
411
31.7981
31.6899
31.7399
31.6099
31.7199
31.6999
31.7699
332
35.1541
35.9199
35.2099
35.0400
35.1699
35.1299
35.2500
340
37.4744
37.8400
37.4199
37.4899
37.4799
37.4599
37.4399
334
43.7023
42.0099
43.5900
43.5200
43.6800
43.6699
43.6800
235
46.2016
45.0699
45.6299
46.1699
46.2199
46.2500
46.2899
541
48.5724
49.4399
48.7500
48.3899
48.4399
48.2099
48.6699
444
51.9261
50.6699
51.5299
51.5799
51.9599
51.8799
51.9699
543
52.9963
52.1699
52.5599
52.7999
53.0099
52.8499
52.9300
155
53.5242
53.1299
53.1499
53.2500
53.5499
53.3699
53.6199
255
55.0759
54.2099
55.0399
55.5699
55.0399
54.5399
55.1199
535
57.5693
57.9499
57.5599
57.8599
57.5399
58.3499
57.5099
554
60.8887
59.7500
60.0099
59.8899
61.0000
59.6999
60.8599
555
64.9076
64.9099
64.9899
64.9099
64.9000
64.9099
64.9199
where Fm„ is the resonant frequency for the m order FDTD algorithm and cell size of n mm
and Fr is the exact value o f the resonant frequency.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
38
Table 2-2 Comparison o f percentage error in resonant frequencies obtained from FDTD
code (second and fourth orders) when choosing different cell sizes (1, 2, and 3.33 mm).
Mode#
F IG
Hz)
PE*
PE*
PEv
p e
42
PE*
PE*
110
10.5993
0.37172
0.47833
0.08774
0.85477
0.00566
0.38304
112
18.3586
1.19127
0.21080
0.42869
0.49731
0.10132
0.27943
310
23.7001
5.90799
0.88691
1.73079
0.46329
0.50717
0.37890
311
24.8577
4.69794
0.99688
1.64054
0.25063
0.39304
0.25063
411
31.7981
0.34027
0.18303
0.59186
0.24592
0.30882
0.08869
332
35.1541
2.17841
0.15873
0.32456
0.04495
0.06883
0.27280
340
37.4744
0.97560
0.14543
0.04136
0.01468
0.03869
0.09206
334
43.7023
3.87257
0.41714
0.41714
0.05103
0.07414
0.05103
235
46.2016
2.44948
1.23740
0.06861
0.03961
0.10476
0.19112
541
48.5724
1.78599
0.36564
0.37573
0.27279
0.74631
0.20073
444
51.9261
2.41921
0.76301
0.66672
0.06509
0.08897
0.08434
543
52.9963
1.55936
0.82346
0.37059
0.02566
0.27625
0.12510
155
53.5242
0.73667
0.69931
0.51229
0.04801
0.28828
0.17880
255
55.0759
1.57238
0.06536
0.89694
0.06536
0.97320
0.07989
535
57.5693
0.66112
0.01633
0.50478
0.05107
1.35593
0.10318
554
60.8887
1.87013
1.44329
1.64037
0.18279
1.95241
0.04729
555
64.9076
0.00354
0.12679
0.00354
0.01171
0.00354
0.01895
where P E is the percentage error in resonant frequency for the m order and cell size of n
mm.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
39
size used is 3.33mm. The domain size is 15x15x15 cells, time step is set at 8192 steps, and
the size of the time step is increased to 6.4187ps. The highest frequency deviation when using
the YA is 5.90799, while it is 1.44329 when using the fourth order FDTD technique.
In Table 2-3, a comparison between FDTD techniques (second and fourth orders)
taking into consideration the cell size, the number of cells in free space used in the previous
three cases, the time taken by the CPU in the computation, the saving in space and the highest
percentage deviation in computing the resonant frequencies.
Table 2-3 A comparison between FDTD techniques taking into consideration the cell size,
the number of cells in free space, the CPU tim e, the saving in space and the maximum PE.
2 nd
Order of accuracy
2 nd
4*
2 nd
4*
4th
Cell size (mm)
3.33
3.33
2
2
Number of cells in
216
216
1000
1000
8000
8000
1
1
free-space
Saving in space
1
1
4.6
4.6
37
37
CPU time (sec.)
1.39
11.09
2.54
31.49
8.32
176.60
Max. PE
5.90799
1.44329
1.73079
0.85477
1.95241
0.38304
One can see from the data shown in Tables 2-1 through 2-3 that the maximum
percentage error of the fourth order FDTD scheme for any cell size is smaller than that of the
YA. Moreover, one can use the (4x4) scheme to get results with better accuracy than that
from YA and with reduction in the computational domain by a factor of 37
with
approximately the same CPU time. This is clearly evident from the case when the cell size is
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
1mm for YA (2x2) and the case when the discretization in space is 3.33mm for the fourth
order technique (4x4).
2.5.2 Free Space
The total field formulations concerning the propagation in free space case are
discussed in sections 2.1,2.2 and 2.3 for different order schemes (2x2), (2x4), and (4x4). In
a free space region, the corresponding equations as discussed in sections 2 .1 are programmed
for different order schemes (2x2), (2x4), and (4x4), respectively. A computer program based
on FDTD has been developed to test the fields at different points in the domain of
computation.
In this application, an electric current source with a Gaussian pulse waveform is
chosen, and it is given by
(2-43)
m Ax
2 c
(2-44)
to - n T
where Ja is the peak amplitude of the current.
JQ , m, and n are chosen to be 5
amperes/meter2, 5, and 2, respectively, m and n are chosen to control the pulse width and the
effective start time of the Gaussian pulse. Equation (2-43) can then be written as:
(2-45)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
41
This pulse is added to the expression o f the electric field in the z-direction.
Figure 2.3 shows a plot ofE z as function of time at different observation points when
the source is located at (50,50,50). The field is computed using second and fourth order
FDTD schemes. The three dimensional space is of volume of 100x100x100 layers that have
sides equal to 4 mm, and the time increment is 7.7 ps. The medium is free space. In this
figure, the (4x4) scheme gives very low distortion levels as compared with the (2x2) scheme
without using more memory size during the computation.
The shapes of a Gaussian pulse propagated in free space at several numbers o f time
steps are shown in Fig. 2.4. In this figure, E2 is plotted versus delta-x position at time steps
20, 35, 50, 65, 80, and 95 using the same parameters mentioned for Fig. 2.3.
In the following case, the current source mentioned before is used here with an
amplitude o f 10 amperes/meter2, and it is added to the expression of Ez. Figure 2.5 shows
a plot of E. as function of time computed using (2x2), (2x4), and (4x4) order FDTD
schemes. The field is measured at different test points, (35,25,25) and (40,25,25), when the
excitation point positioned at (25,25,25). The size of the chosen domain is 50x50x50 cells and
the cell length in x, y and z- directions are 4 mm, and the size o f the time step is 8 .47 ps. In
Fig. 2.5, the level of distortion decreases as the order of accuracy increases without using
more memory in the higher order cases. In this case there is no interface between higher order
(4x4) or (2x4) and the lower one (2x2) as discussed before by applying the forward difference
approximation (FDA) and the backward difference approximation (BDA) closer to the
appropriate boundaries. The (2x4) FDTD order scheme starts to diverge after reaching the
steady state. When the increment in time (At=5.929 ps) is decreased by 30% , the fields of
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
42
0.03
0.02
—
2 ’4
—
4 u order
order
0.01
0
-0.01
-0.02
-0.03
0
0.2
0.4
0 .6
0.8
1
0.8
1
time(ns)
(a)
0.03
0.02
—
2 *‘ order
—
4 “ order
0.01
0
-0.01
-0.02
-0.03
0
0.2
0.4
0 .6
time(ns)
(b)
Fig. 2.3 Ezversus time calculated in a domain o f size of 100x100x100 layers with sides equal
to 4mm using second and fourth order FDTD schemes when a current source with
a Gaussian pulse waveform is applied at the point (50,50,50) in the z-direction.
Ez is computed at the points: (a) (60,50,50), (b) (65,50,50), (c) (70,50,50), (d)
(75,50,50), (e) (80,50,50), and (f) (85,50,50).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
43
0.03
0 .0 2
—
2 ’ " o rd e r
—
4 * order
0 .0 1
0
-0 .0 1
-0 .0 2
-0.03
0
0 .2
0 .6
0.4
time(ns)
0 .8
(c)
0.03
2 *" order
4 “ order
0 .0 2
0 .0 1
0
-0 .0 1
-0 .0 2
-0.03
0
0 .2
0 .6
0.4
time(ns)
0 .8
(d)
Fig. 2.3 Ezversus time calculated in a domain of size o f 100x1OOx100 layers with sides equal
to 4mm using second and fourth order FDTD schemes when a current source with
a Gaussian pulse waveform is applied at the point (50,50,50) in the z-direction.
Ez is computed at the points: (a) (60,50,50), (b) (65,50,50), (c) (70,50,50), (d)
(75,50,50), (e) (80,50,50), and (f) (85,50,50).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
44
0.03
—
—
0 .0 2
E
>
2 14 order
4 * order
0 .0 1
N
111
-0 .0 1
-0 .0 2
-0.03
0
0 .2
0 .6
0.4
time(ns)
0 .8
1
(e)
0.03
—
—
0 .0 2
e
>
o.oi
m
-0 .0 1
2 *4 order
4 Ik order
-0 .0 2
-0.03
0
0 .2
0.4
0 .6
time(ns)
0 .8
(0
Fig. 2.3 Ez versus time calculated in a domain o f size of 100x100x100 layers with sides equal
to 4mm using second and fourth order FDTD schemes when a current source with
a Gaussian pulse waveform is applied at the point (50,50,50) in the z-direction.
Ez is computed at the points: (a) (60,50,50), (b) (65,50,50), (c) (70,50,50), (d)
(75,50,50), (e) (80,50,50), and (f) (85,50,50).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
45
E
>
0.06
0.05
0.04
0.03
—
—
2 nd order
4 ^ order
0 .0 2
0 .0 1
N
LU
0
-0 .0 1
-0 .0 2
-0.03
50
60
70
80
position
90
100
90
100
(a)
E
>
N
LU
0.06
0.05
0.04
0.03
—
—
2 no order
4 th order
0 .0 2
0 .0 1
0
-0 .0 1
-0 .0 2
-0.03
80
70
position
(b)
Fig. 2.4 The shapes of a Gaussian pulse propagated in free space using the same parameters
mentioned for Fig. 2.3. Ez versus delta-x position after various numbers of time
steps (n) when a current source o f a Gaussian pulse wave form is applied at
(50,50,50). (a) n=20, (b) n=35, (c) n=50, (d) n=65, (e) n=80, and (f) n=95.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
46
-
e
*>
4 order
N
LU
-0 . 0 2
position
E
>
M
UJ
0.06
0.05
0.04
0.03
0 .0 2
0 .0 1
—
--
2 no order
4 ^ order
U
I
0
-0 .0 1
-0 . 0 2
-0.03
|
50
60
70
80
position
90
100
(d)
Fig. 2.4 The shapes o f a Gaussian pulse propagated in free space using the same parameters
mentioned for Fig. 2.3. Ez versus delta-x position after various numbers of time
steps (n) when a current source of a Gaussian pulse wave form is applied at
(50,50,50). (a) n=20, (b) n=35, (c) n=50, (d) n=65, (e) n=80, and (f) n=95.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
47
0.06
0.05
0.04
.—
. 0.03
E
;> 0.02
3
0.01
uT
0
-0.01
-0.02
-0.03
— 2 nd order
- - 4 th order
50
60
70
80
position
90
100
90
100
(e)
e
>
0.06
0.05
0.04
0.03
— 2 nd order
- - 4 th order
0 .0 2
0 .0 1
N
LU
0
-0 .0 1
-0 .0 2
-0.03
50
60
70
80
position
(f)
Fig. 2.4 The shapes of a Gaussian pulse propagated in free space using the same parameters
mentioned for Fig. 2.3. Ez versus delta-x position after various numbers of time
steps (n) when a current source of a Gaussian pulse wave form is applied at
(50,50,50). (a) n=20, (b) n=35, (c) n=50, (d) n=65, (e) n=80, and (f) n=95.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
48
e
0.06
(2x2) order
0.04
•— ■ (2x4) order
(4x4) order
0 .0 2
0
LU
-0 . 0 2
-0.04
-0.06
0
0 .2
0.4
0 .6
time (ns)
0 .8
1
0 .8
1
(a)
0.04
—
(2x2) order
(2x4) order
0 .0 2
(4x4) order
0
-0 . 0 2
-0.04
0
0 .2
0.4
0 .6
time (ns)
(b)
Fig. 2.5 Ez versus time calculated in a domain o f size 50x50x50 layers with side equal 4mm
using (2x2), (2x4), and (4x4) order FDTD schemes when a current source with a
Gaussian pulse wave form is applied at the point (25,25,25) in the z-direction.' E z
is computed at the points: (a) (35,25,25), (b) (40,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
49
the three different orders are very near to each other as shown in Fig. 2.6. On the other hand,
when the increment in time is increased by 10% of that of the first case shown in Fig. 2.5, the
higher orders still have lower distortion than that of the (2x2) order as shown in Fig. 2.7. This
means that, the results are more stable as the order becomes higher, but the (2,4) order still
starts to diverge after the steady state. This divergence disappears when the increment in time
is reduced.
2.5.3 Lossy and Lossless Dielectric Structures
The finite difference equations of the total field in a dielectric region is discussed in
sections 2.1. In this section, these equations are programmed for the fourth order scheme
(4x4), the second order scheme (2x2), and the second order in time and fourth order in space
(2x4), respectively.
At first, the program is tested for different of erand fu sin g a current source with a
Gaussian pulse added to E. expression with amplitude 5 amperes/meter2 as mentioned
previously. Figure
2 .8
shows the variation of Ez either with delta-x position or with time
using different values o f er when the second order accuracy (2x2) is applied, while Fig. 2.9
illustrates E. versus delta-x position and time, respectively, with different er using the fourth
order scheme (4x4) with the same excitation source. The variation of £ ;with delta-x position
after 30 time steps with different at using (2x2) and (4x4) order schemes are shown in Fig.
2 .1 0 .
A comparison between (2x2) and (4x4) is performed when a current source with a
Gaussian pulse waveform is applied in the z-direction and added to the expression of E. at
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
50
0.04
0.03
—
(2x2) order
(2x4) order
0 .0 2
(4x4) order
0 .0 1
0
-0 .0 1
-0 . 0 2
0
0 .2
0.4
time (ns)
0 .6
(a)
(2x2) order
(2x4) order
(4x4) order
0 .0 2
0 .0 1
0
-0 .0 1
-0 . 0 2
0
0 .2
0.4
time (ns)
0 .6
(b)
Fig. 2.6 Ez versus time calculated in a domain size of 50x50x50 layers with sides equal to
4mm using (2x2), (2x4), and (4x4) order FDTD schemes when a Gaussian pulse
as a current source is applied at the point (25,25,25) in the z-direction when At
= 5.929ps. Ez is computed at the points: (a) (35,25,25), (b) (40,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
51
0.08
0.06
0.04
E*
~
0 .0 2
0
m~ -0 . 0 2
-0.04
-0.06
0
0.2
0 .4
0.6
time (ns)
0.8
1
0.8
1
(a)
0.08
0.06
0.04
E*
0.02
o
-0 . 0 2
-0.04
-0.06
0
0.2
0 .4
0.6
time (ns)
(b)
Fig. 2.7 Ez versus time calculated in a domain size of 50x50x50 layers with sides equal to
4mm using (2x2), (2x4), and (4x4) order FDTD schemes when a Gaussian pulse
as a current source is applied at the point (25,25,25) in the z-direction when
At = 9.317ps. E 2 is computed at the points: (a) (35,25,25), (b) (40,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
52
0.04
0 .0 2
e
>
LU
-0 . 0 2
-0.04
30
35
40
position
45
50
0.3
0.4
(a)
0.08
0.06
e
0.04
~
0 .0 2
LU
-0 . 0 2
0
0 .1
0 .2
time (ns)
(b)
Fig.
2 .8
E, at different values of e, using (2x2) FDTD order scheme: (a) E, versus delta-x
position, (b) E z versus time.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
53
0.04
Ez (v/m)
0 .0 2
-0 . 0 2
-0.04
45
40
35
30
position
(a)
0.08
Ez (V/m)
0.06
0.04
0 .0 2
-0 . 0 2
0
0 .1
0 .2
0.3
0.4
position
(b)
Fig. 2.9 Ez at different values of er using (4x4) FDTD order scheme: (a) Ezversus delta-x
position, (b) Ezversus time.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
54
0.02
L
—
\
\
e
0
>
/
\
UJ
-0.02
30
----- a ° =0'01
=0.10 (S/m)
f
\
N
<Je =0.00 (S/m)
/
35
40
45
position
(a)
0 .0 2
—
a e =0.00 (S/m)
—
a e =0.05 (S/m)
■— ■ a e =0.10 (S/m)
E
>
LU
-0 . 0 2
30
35
40
45
position
(b)
Fig. 2.10 Ez as function of delta-x position after 30 time steps with different oe using:(a)
(2x2) order scheme, (b) (4x4) order scheme.
Reproduced with permission of the copyright owner. Further reproduction prohibited w ithout permission.
55
(40,40,40) in a domain of size 80x80x80 when a lossless dielectric slab o f e=2.Q is located
between delta-x positions 50 and 60 and both delta-y positions and delta-r positions between
20 and 60. The response resulting from the (2x2) and (4x4) formulation is shown in Fig.
2.11. The same thing is done for a lossy dielectric as illustrated in Fig. 2.12. The increment
in time in this case is 7.698ps and the cell size is 4mm. The fourth order scheme (4x4) always
has low ripples compared with the second order scheme.
Finally, a comparison between (2x2), (2x4), and (4x4) formulation is done using a
current source with a Gaussian pulse waveform at the point (50,50,50) in the z-direction in
a domain size of
1 00 x 100x 100
cells with a lossless dielectric slab (^ = 2 .0 ) located between
delta-x positions 70 and 80 and both delta->> positions and delta-z positions are 20 and 80. Ez
is plotted versus time at different test points as shown in Fig. 2.13. The distortion of (2x4)
scheme is smaller than (4x4) scheme in this case but the (2x4) scheme starts to diverge after
reaching the steady state at the points preceding the dielectric interfaces.
2.5.4 Lossy and Lossless Magnetic Structures
The equations corresponding to the case of lossy and lossless magnetic structures are
discussed in sections 2.1. A computer program based on FDTD schemes for (4x4), (2x2), and
(2x4) formulations has been developed to test the field at different points in the domain.
At the beginning, the program is tested for different values of pr and omusing a
magnetic current source with a Gaussian pulse waveform which is applied at (15,25,25) and
added to the ^expression with amplitude 10*volts/meter2. This magnetic current is given by:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
56
(2x2) order
(4x4) order
-
0.02
0.6
0.8
time (ns)
(a)
0.02
E
>
—
(2x2) order
—
(4x4) order
0
fM
LU
-
0.02
0
0.2
0.4
0.6
0.8
1
time (ns)
1.2
1.4
1.6
1.4
1.6
(b)
0.01
E
>
—
(2x2) order
— • (4x4) order
0
LU
-
0.01
0
0.2
0.4
0.6
0.8
1
time (ns)
1.2
(c)
Fig. 2.11 A comparison between (2x2) and (4x4)is performed when a current source with
a Gaussian pulse waveform is applied at (40,40,40) in a domain size of 80x80x80
when a lossless dielectric slab o f er=2.0 is located between delta-x positions 50
and 60 and both delta-y positions and delta-z positions 20and 60. (a) Ezversus
time at (45,40,40). (b) Ez versus time at (55,40,40). (c) Ez versus time at
(65,40,40).
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
57
_
0.06
0.04
-i
0.02
>
— (2x2) order
(4x4) order
N
LU
-
0.02
-0.04
0
0.1
0 .2
0.3
0.4
0.5
time (ns)
0.6
0.7
0 .8
(a)
0.02
_ 0.01
E
>
0
-
0.01
-
0.02
0
0.1
0.2
0.3
0.4
0.5
time (ns)
(b)
—
(2x2) order
—
(4x4) order
0.6
0.7
0 .8
0.01
—
_ 0.005
(2x2) order
*■•• (4x4) order
m-0.005
-
0.01
0
0.1
0.2
0.3
0.4
0.5
time (ns)
0 .6
0.7
0.8
(c)
Fig. 2.12 A comparison between (2x2) and (4x4) is performed when a current source of a
Gaussian pulse waveform is applied at (25,25,25) in a domain size of 50x50x50
when a lossy dielectric slab o f ae=0 . 1 S/m and e = 2 .0 is located between delta-x
positions 30 and 40 and both delta-y and delta-z positions lOand 40. (a) Ez
versus time at (30,25,25). (b) Ez versus time at (35,25,25). (C) Ez versus time
at (40,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
58
E, (Vim)
0.04
— (2x2) order
-— ■ (2x4) order
—— (4x4) order
A
0.02
/
1
<
>
0
-
I f 'j l
1
0.02
------ ,-------- , --------,---------,--------- 1---------,
0 .2
0
0.4
. ,
0.6
,
.
i
0.8
1
I
>
! i:
1 .2
1.4
1 .2
1.4
1 .2
1.4
time (ns)
(a)
0.02
E, (Vim)
0.01
-
0.01
-
0.02
0 .2
0
0.4
0 .6
0 .8
1
E* (V/m)
time (ns)
(b)
— (2x2) order
■— ■ (2x4) order
— (4x4) order
01
0
01
0
0 .2
0.4
0 .6
0 .8
time (ns)
(c)
Fig. 2.13
A comparison between (2x2), (2x4), and (4x4) using a current source with
a Gaussian pulse waveform in the z-direction at the point (50,50,50) in a
domain size of
located
10 0 x 1 00 x 100
cells with a lossless dielectric slab (er=2 .0 )
between delta-x positions 70 and 80 and both delta-y and delta-z
positions at
20
and 80. Ezis plotted versus time at different test points: (a) at
(60.50.50). (b) at (65,50,50). (c) at (70,50,50). (d) at (75,50,50). (e) at
(80.50.50), (f) at (85,50,50).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
59
0.01
—
(2x2) order
—— (2x4) a d e r
e 0.005
•— * (4x4) order
uT -0.005 f
-
0.01
0
0.2
0.4
0.8
0.6
1
1.2
1.4
time (ns)
(d)
0.005
E
>
—
—
(2x2) order
(2x4) order
—
(4x4) order
N
UJ
-0.005
0
0.2
0.4
0.8
0.6
1
1.2
1.4
time (ns)
(e)
0.01
—
(2x2 ) order
— > (2x4) order
■— ■ (4x4) order
0
-
0.01
0
0.2
0.4
0.8
0.6
1
1.2
1.4
time (ns)
(0
Fig. 2.13 A comparison between (2x2), (2x4), and (4x4) using a current source with a
Gaussian pulse waveform in the z-direction at the point (50,50,50) in a domain
o f size 100x100x100 cells with a lossless dielectric slab (e r=2.0) located
between delta-x positions 70 and 80 and both delta-y and delta-z positions at
20 and 80. E zis plotted versus time at different test points: (a) at (60,50,50). (b)
at (65,50,50). (c) at (70,50,50). (d) at (75,50,50). (e) at (80,50,50), (f) at
(85,50,50).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where
T = —— , to=nT, m= 5, n -2
(2-46)
The domain size in this case is 50x50x50 cells with cell length 4mm in any direction and the
increment in time is 7.698ps. In addition, another case is studied where a lossless magnetic
slab is located between delta-x positions 30 and 40 and both delta-y and delta-z positions
at 10 and 40. Figures 2.14 and 2.15 show the variation of Hz with time at different testing
points with various values where fixwhen (2x2) and (4x4) FDTD order scheme formulations
are applied.
A comparison between (2x2) and (4x4) is performed with a source has the same
parameters of that used in the last two figures with the same domain. A lossless magnetic slab
where \i=2 is located between delta-x positions 30 and 40 and both delta-y positions and
delta-z positions 10 and 40. The response due to (2x2) and (4x4) is shown in Fig. 2.16. The
same thing is done for a lossy magnetic slab of different a mas illustrated in Figs. 2.17. It is
interesting to note, that the fourth order scheme (4x4) always has lower ripples compared
with the second order scheme.
2.5.5 Planar M icrostrip Circuits
In this section, the formulations of the second order (2x2) and fourth order (4x4)
accurate finite difference algorithms are applied in computer programs based on the
FDTD technique. These programs are used to analyze planar microstrip circuits such as
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
61
Hz (A/m)
4
2
0
-2
-4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1
1.2
1.4
tim e (ns)
(a)
Hz (A/m)
2
1
0
1
2
0
0.2
0 .4
0.6
0.8
tim e (ns)
Hz (A/m)
(b)
1.2
0.6
—
(1,-2
■—' H,-*
0
-
0.6
-
1.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
tim e (ns)
(c)
Fig. 2.14 Hj versus time with different values o f p r using (2x2) order scheme at different
testing points, (a) at (25,25,25). (b) at (35,25,25). (c) at (45,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
62
If
4
2
0
x" -2
-4
0.2
0.4
0.6
0.8
time (ns)
1
1.2
1.4
(a)
2
E
<
0
nT
-1
1
-2
0.2
0.4
0.6
1
1.2
1.4
0 .2
0.4
0.6
1
1.2
1.4
0.8
time (ns)
(b)
1.5
1
E 0.5
^
0
~ 0 .5
1
-1
-1.5
0.8
time (ns)
Fig. 2.15 Hj versus time with different values of pr using (4x4) order scheme at different
testing points, (a) at (25,25,25). (b) at (35,25,25). (c) at (45,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
63
Hz (A/m)
4
—
(2x2) order
—
(4x4) order
2
0
-2
-4
0
0 .2
0.4
0 .6
0 .8
1
1 .2
1.4
Hz (A/m)
time (ns)
(a)
2
—
(2x2) order
1
—
(4x4) order
0
-1
2
0
0 .2
0.4
0 .8
0 .6
1
1 .2
1.4
time (ns)
(b)
1.5
—
Hs (A/m)
1
(2x2) order
■— ■ (4x4) order
0.5
0
0.5
-1
1.5
0
0 .2
0.4
0 .6
0 .8
1
1 .2
1.4
time (ns)
(c)
Fig. 2.16 H, versus time when \ir =2 using (2x2) and (4x4) order schemes at different
testing points, (a) at (25,25,25). (b) at (35,25,25). (c) at (45,25,25).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
64
0.01
— am- 0.0 Q/m
0.015
0 .0 0 5
a m - 5 E 3 Cl/m
0.01
0. -0 0 O/m
—
a . -5 E 3 Q /m
0.0 0 5
E
E
5
x*
3
-0 .0 0 5
0
•0.005
-0.01
- 0.01
-0.015
-0 .0 1 5
0
0 .5
1
0
1.5
1
0.5
1.5
tim e (ns)
lim e (ns)
(b)
0 .0 1 5
0.01
—
(2 x2) order
(4 x 4 ) order
0 .0 0 5
0.0 0 5
E
E
3
3
(2x2) order
(4x4) order
0
-0 .0 0 5
-0 .0 0 5
•0 .0 1 5
-
0.01
0
1
0.5
tim e (ns)
1.5
tim e (ns)
(d)
Fig. 2.17 H* versus time at the test point (45,25,25) using: (a) (2x2) order scheme, (b) (4x4)
order scheme, (c) both (2x2) and (4x4) order schemes when om =0. (d) both
(2x2) and (4x4) order schemes when om=5x103 .
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
65
transmission lines and low-pass filters. The characteristics of the propagated waves on these
types of microstrip circuits are computed and analyzed.
2.5.5.1 Modeling of a Resistive Voltage Source
In this part, the fourth order accurate FDTD Maxwell’s equations are derived to
model a lumped-circuit load (resistors) and sources (resistive voltage sources) in 3-D. This
fourth order algorithm is deduced form the second order formulation discussed in [33],
When using a resistive voltage source, consider Maxwell’s curl H equation which is
suitable for the time stepping of the electric field components
^ s
7
8D
Vs n
V*H = J + — + -------dt
R AA
(2'47)
where n is the unit vector in the direction of the electric current density J„ AA is the element
area normal to Ja J = o t E, D = eE is the electric displacement and Rs is the resistance of the
voltage source. Using the central difference approximation for both space and time
coordinates, equation (2-47) reduces for the Ex field component as [33]
/
a A t)
1 — i ---2 €
o At
1 +—---2 e ;
p
( IjJc ) +
Al
1
n->—
e
2
V x H,( ijjc )
a At
1 +—
2 e )
At
R s e Ay Az
a At
1 + - *---2 e
(2-48)
since R = p L / A = L/(o, A) = A x / (a. Ay A z ) and this gives a, =Ax/ ( R Ay Az) (where R,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
66
p, and ot are the resistance, resistivity, and electric conductivity). This expression of ot can
be used in (2-48) considering R = Rs for matching. Thus, equation (2-48) becomes
1
-X( Ijjc )
At Ax
2 R s e Ay Az
At Ax
1 +2 R e Ay Az
*
J
)
E Xn( i j j c
)
At_
€
At Ax
1 +2 R e Ay Az
1
J
/
1
n+ —
2
( V* )
At
R s e Ax Ay
(2-49)
At Az
1 +2 R se Ax Ay
• 'i
In this case, the curl operator depends on the order of the finite difference accuracy. For the
second order accuracy;
1
n* —
1
n-—
H
1
n» —
h,
V x
2
- H
2
( V* )
1
n*—
2
(
)
H
y
(
1
2
n~—
2
- H
>• (
)
Ay
)
Az
(2-50)
while for fourth order accuracy,
1
z(
2 x
)
- H
1
2 x
H
z ( IJ--JC )
n
X
H
1
y(
2
i
x
)
Ay
- H
2
Az
l
H
z(
3
)
-3
h "
2 x
+ 3 H " 2 x
Z(
)
~ H* 2 i
Z(
x
y ( v ^ -- )
)
24 Ay
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Using (2-50) into (2-48), one can obtain the corresponding time-stepping relation for E J ij.k )
in second order accuracy and using (2-51) into (2-48), one can obtain the corresponding timestepping relation for Ex(iJ,k) in fourth order accuracy.
In case of a resistive load the same expressions are used with V, =0. If the resistive
voltage source or the load resistance is located in a free space region, then e =
at = 0 ,
J= 0.
2.5.S.2 M icrostrip Transmission Lines
Microstrip lines are used to connect between elements in microstrip circuits. This
transmission line is made of a uniform straight copper conductor placed on top of the dielectric
substrate backed by a ground plane [34]. Microstrip line is terminated at one end with a
resistive voltage source or a current source and on the other end with a resistive load. The
representation o f sources and matching resistive loads are discussed in section 2.5.5. 1 . The
value of the load is normally close to the characteristic impedance of the line at the center
frequency of the excitation waveform frequency spectrum. In this investigation, Mur ABC,
Liao ABC, or four cells of a perfect electric conductor (PEC) are used at the outer boundary.
Different excitation waveforms are used as a source to the microstrip line such as a
Gaussian pulse as a voltage source or current source and a sinusoidal source. Each of these
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
68
cases is discussed and analyzed. The voltage is computed at different points on the microstrip
line with the different FDTD orders of accuracy. Normally, the fourth order FDTD algorithm
is used in the computation for positions more than the first two cells away from the boundary,
but for the ground plane, it is recommended to use at least four cells, if it is connected directly
to the absorbing boundary (AB) to insure no divergence. This will enable the use of the fourth
order accuracy directly without using second order in the first four PEC cells closer to the
boundary of the domain, because the field inside these cells is zero. Generally, it is better to
have no more than one medium in the first four cells connected to the AB; otherwise there will
be a divergence in the fourth order scheme results.
In the following example, a resistive voltage source with a Gaussian pulse waveform
is used and defined as:
V (t)
lV J
V=
10 volts,
o
=
V
0
e ' [ ( r ' ,o) l,mtx]2
tm ax = 0.025 ns,’
(2-52)
to = 31max
As an example, a microstrip transmission line is modeled in a computer program based
on FDTD algorithms with Ax =0.508 mm, Ay = 0.626533 mm, Az = 0.5 mm, At = 1.033208
ps, and the total time step is 1024. The problem space is 36 x 57 x 169 cells. The transmission
line has length of 70 mm and has a dielectric substrate o f thickness o f 1.524 mm of relative
permitivity of 3.05. The strip is made of a thin PEC and connected at one end with a resistive
voltage source with a Gaussian pulse waveform of amplitude 10 volts and a source resistance
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
10
1 .5 2 4
e, - 3 .0 5
0 .5 0 8
3 .7 6
1 3 .1 2
70
port 1
port
R
V„
T
10
Fig. 2.18 Geometry of a microstrip transmission line (all units are in mm and the crosssections are not drawn to scale)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
70
volt (V)
6
—
4
(2x2) order
(4x4) order
2
0
0
0.1
0.2
0 .4
0.3
0.5
0.6
0.5
0.6
0.5
0.6
tim e (ns)
(a)
volt (V)
6
(2x2) order
4
(4x4) order
2
0
0
0.1
0 .2
0.4
0.3
time (ns)
(b)
volt (V)
6
—
4
(2x2) order
■— ■ (4x4) order
2
0
0
0.1
0.2
0.3
0.4
time (ns)
(c)
Fig. 2.19 Voltage versus time at different test points on the transmission line computed using
(2x2) and (4x4) order schemes where the source is a voltage source of a Gaussian
pulse waveform. The test points are: (a) (7,29,15). (b) (7,29,45). (c) (7,29,75).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
71
of 50 Q. The other end of the line is connected to a resistive load equal to 50 Q as shown in
Fig. 2.18. Figure 2.19 illustrates the results from the program using (2x2) and (4x4) order
schemes. These results show that the fourth order almost has the same response as the second
order scheme. There is a DC shift in the steady state part o f the result computed using (4x4)
order scheme at some test points as shown in Fig. 2.19. This DC shift in the steady state part
of the result appears only when using composite structures, where multiple reflections occur
and standing wave patterns may develop between media interfaces and due to long terms in
(4x4) scheme where the error can be accumulated rapidly.
On the other hand, when the excitation source is a current source with a Gaussian
pulse waveform is applied at the source point (the near end o f the PEC thin strip). The current
source in this case is given by
/ Xv
( /' ) = / 0
e
' l(' ‘ ° " “ J'
(2-51)
where
m ax
The DC shift in the steady state part of the result appears in this case at the test point near
the source edge as shown in Fig. 2.20. The last test case is performed by using the derivative
of the voltage source with a Gaussian pulse waveform as a source to the microstrip line and
it is given by
m ax
where
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-52)
72
2e-4
—
(2x2) order
(4x4) order
r
o>
1 e *4
0
0.1
0 .2
0 .3
0 .4
0.5
0 .6
tim e (n s)
(a)
2 e -4
— 1 e -4
0
0.1
0 .2
0 .3
—
(2x2) o rd er
■—
(4x4) order
0 .4
0.5
0 .6
tim e (n s)
(b)
2 e -4
—
(2x2) ord er
(4x4) ord er
— 1 e -4
0
0
0.1
0 .2
0 .3
0 .4
0.5
0 .6
tim e (n s)
(c)
Fig. 2.20 Voltage versus time at different test points on the transmission line computed using
(2 x2 ) and (4x4) order schemes when the source is a current source of a Gaussian
pulse waveform. The test points are: (a) (7,29,15). (b) (7,29,45). (c) (7,29,75).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
73
20-11
1e-11
0
0.2
0.1
0.3
—
(2x2) order
—
(4x4) order
0.4
0.5
0.6
tim e (ns)
(a)
2 e -1 1
1e - 1 1
>
—
(2x2) order
—
(4x4) order
o>
0
0.1
0.2
0.3
0 .4
0.5
0.6
tim e (ns)
(b)
2 e -1 1
—
1e - 1 1
(2x2) order
— ■ (4x4) order
0
0.1
0.2
0.3
0 .4
0.5
0.6
tim e (ns)
(c)
Fig. 2.21 Voltage versus time at different test points on the transmission line computed using
(2x2) and (4x4) order schemes when the source is a voltage source with the
derivative of a Gaussian pulse waveform. The test points are: (a) (7,29,15). (b)
(7,29,45). (c) (7,29,75).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
74
The response of the fourth order in this case is coincident with that of the second order
without any shift in the steady state part as illustrated in Fig. 2.21. This means that the fourth
order has less accuracy than the second order for such type o f structures and the same
conclusion is observed by C. Manry et al. [35],
2.S.5.3 Microstrip Low-Pass Filter
The second and fourth order FDTD schemes are also used for the analysis of a lowpass filter (LPF). These schemes are applied to compute the reflection and transmission
coefficients of the microstrip circuit shown in Fig. 2.22. This LPF is printed on a dielectric
substrate backed by a PEC ground plane. The number of cells for free space above the
dielectric substrate is sufficiently large to have no reflection from the top AB. In this case, four
cells of PEC are used at the outer boundary. This circuit is modeled using FDTD algorithm of
second and fourth order accuracy with Ax = 0.203 ram, Ay = 0.4233 mm, Az = 0.203 mm, and
At = 0.441 ps [36-37], The problem space size is 33 x 97 x 249 cells, and the total time steps
taken is 4096. The grid spacing based on the dimension of the horizontal strip (resonant strip),
since the resonant frequency o f this filter is sensitive to dimension error of this strip. The timedependent line voltages (power) Vt (P,) and V2 (Pj) at the input and output ports (both of
them are ten cells from the horizontal strip) are computed by second and fourth order FDTD
algorithm. These values are transformed to the frequency domain using Fourier transformation
to compute the S-parameters. The reflection coefficient (Sn) seen at portl when port 2 is
terminated in a matched load (Z 0 = 50 Q) is given by [38]:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
75
2^18 njm.
i
2.54 mm
T
i
0.794 mm
~
r
Fig. 2.22 Microstrip low-pass filter configuration (cross-sections are not drawn to scale).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where V,* (P,~) and Vf (P f) are the incident and reflected wave voltage (power),
respectively. The transmission coefficient (S21) from port
S2l ( dB ) =
20
log 10
V~
—7 1 =o
V{ 2
=
10
1
Iogio
to port 2 is given by
Pi
— Ip ;. o
Pi
Su and S21 are plotted as function of frequency in the range between 0 and 20 GHz. The results
from the second and fourth orders are basically the same. This means that the fourth order
scheme does not improve the accuracy in such guided wave structures. These results are in
close match with those computed and experimentally measured data by Sheen et al. [37]. The
results are plotted in Fig. 2.23.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
77
-10
-20
CO
-30
—
voltage rato
—
■
power ratio
experimental results
-40
F (G Hz)
(a)
04-20
CO
—
—
-30
■
voltage ratio
power ratio
experimental results
-40
0
2
4
6
8
10
12
14
16
18
F (G H z)
(b)
Fig. 2.23 The experimental [37] and numerical values of Sn and S21 in dB for the low-pass
filter, (a) reflection coefficient, (b) transmission coefficient.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
78
Chapter 3
TOTAL FIELD ANALYSIS IN TWO DIMENSIONAL SPACE
USING H IG H ER ORDER FDTD TECHNIQUES
In this chapter, the total field formulations of fourth order accurate finite difference
algorithm in space and time (4x4) are derived. The corresponding formulations of second
order in time and space (2x2) are derived from the fourth order formulations. The
characteristics of these different formulations for the two order schemes are analyzed and
compared through different applications. The field is tested in a domain when the excitation
is a current source with a Gaussian pulse waveform that is modulated with and without a
sinusoidal wave. This is done also, when the source has a low or a high frequency to check
the difference between the two FDTD algorithms and to define the regions in which the (4x4)
order scheme gives better results than the (2x2) algorithm. Numerical dispersion and stability
conditions are also discussed.
3.1 Total Field Formulation of Fourth Order FDTD Scheme
In order to reduce the complexity of programming and displaying FDTD
computations, the total field is formulated in two dimensions for TIvfmode. This means that,
H„ Hy and E. components are considered only. The total field expressions derived in the
three dimensional space can be used to develop the corresponding formulations in this case
by canceling Hv Ex, and Ey and neglecting the derivative concerning the z-direction. Using
equations (2-18), one can get:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
t~
^2
z ifl
E --vj>
Til
1
n»—
T]
_
r-n
Z(<,°
PAH /
T1ier
2
1
/if—
1
3x ~ H
2
2
\
y
)
I
/if—
r (/f —J-—)
2
2
1
1
2 4 tije r
y
2
2
H
2 4 TijHre r
3
2
2
2
/ 3
-
(jy r *
py
2
3x - 2 H / 1
lx
2J~2
- H
2 4 T i,e r
!
u)
2 ,/" 2
i. )
y (.‘*i
--j--)
2
2
i
n-*-—
- 2
H
2
/
1
2
1*
2
- H
x
2 1^ f —i ). )
2 2
1
1
^ Py
2i
y
n f—
2
« f—
« f—
x«/ V3 -T)
#
1
1
n f—
2 4 T ij^ e J
2
23
i
r tf—
Px
~ H
l
2
1
3 i
y ("W
y e n )
y
n
1
2 H
y
rtfi
,
\
1
nf—
2
1
n»—
3 H
K
y (,- ^ a >
,
2
r tfl
2
3
1
/if—
T1 Py Px
2 4 T i^ e J
.
1,
H
.1
1.
y O -^-)
1
/if—
2
1
/if—
2
2
1
/if—
n Px
1
/if—
/if—
n*<
2
*1 Px
2
2
lx -
1J - — )
x (, / f —
2
2
1
rtf—
2
1
n- —
1
rtf—
2
-
3 H
2
X (,fi,-i)|
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
80
1
2
H
/
24 i W ?
ri p
o
x- . ( _ l
2 4 t i 1e r
e
1
• 3i
2
2
3H
\
r
i
2
+ 3
2
+ -H i )
p
1
x /(ITzrJ*-z)
2
2
1
1
2
2
- H
1
(c A Q V eawpx
^
i -H
n -i.
\H
- 3
2
, 1
x (i—
^)
2
2
i
n- —
1.
. 1
<
)
2f
1
2
2
,-i
\
x
- H
\
a eA/
x
1
H
248T1J
;
(V )
1
*2
X(' V * I )
=- H
\
x
5, X(' V * I }
24ri^pr
^
[ py ( £ r (v-1)
3 E z { l J , X) + 3 £ . (y)
(v*2)
. 2)
- 3
^ (v .)+
- * W> ) ]
E , {IJ. X)
3 e:
{ij, -
X)
E z ( t _x
e:
^ X)
24ti5,8rpr
PyPx
(
(i-ij*D
i.)
rt-l
e
where rj, and rj2 are defined by equations (2-17) and (2-18), and
2
f 2i
\
i
y o-fr-j)
24T|1ptre^
n» —
- H
1
2i
(c A/)2'naeCTmp>,
1
n» —
2
n»—
a
24tl1e ^ r
1
1
n»—
*1 P y
2 E z(
X) )
24t|5iPrer
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3-1)
81
(cA/)2a gamp>
--I
gmA/
(V • « ] - 24p5,
-r (v > l)
2 4 Tl5,Prer
(3- 2)
2
2
where ^ and £2are defined by equations (2-22) and (2-23),
1
n-r—
H /
J'
2
1
1,
2
2
2 4 tl^ p r
H
—
^ z ( i - 2j )
Px
2 1
z ('*!/> - E.
+
3 ^ z (/»!/) + 3 ^ z (v)
rn
^z (i+2j)
2 4 q^P rer
i
3
+
£ zO * W )
3
E z
(V )
PxP"
0 - lv ) ]
'
^zo-w-l)
24qC,tirer
/■
24q £,pr
( ~"
e
4.
( c A/ ) 2 o<qmpx
24q^^er
p
)\
E'n
(
\
ITn *1
z 0*1,/) _
(i'l j)
E'”
z (v)
_
r n *1 \ _ /
Ez
(v) ) (
om
wAr
245,p
r ,rt
2 (i-lj)
_
r -'1
E.
(,^j \)
"‘i*
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3-3)
82
One can substitute e = eQ, // = fta, and at = cr m=0 in the equations of this section to get
the total field difference equations of fourth order scheme (4x4) in free space.
3.2 Total Field Difference Equations of Second Order FDTD Scheme
By taking into consideration the second order accurate equations derived in section
2 .1
one can obtain:
^ n+1
(V )
•.J.
Vo
(e+a At)
pn
£
(e+a At)
-J.
y - * y70-l- j ' -)l J
y O'-j*-)
(3-4)
2
2
Similarly, ^Tr and i f components will take on the following forms:
1
n- —
H
H
2
c(p+omAr)
2
1
n -—
H
2
z (v^l) - E" (V)
1
n- —
2
y o-fri)
(p+amA/)
H
2
Pz
c(P + a ™A 0
r*"
_ rr'1
X(;-!/) t-'x (ij)
(3-5)
(3-6)
The total field difference equations of second order scheme (2x2) in free space can be
obtained from the equations of this section by substituting by e = eQ, n = n Q, and ot = am=0 .
3.3 Node Distribution of the Second and Fourth Order Schemes
For the fourth order scheme (4x4) where the error term in the approximation of the
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
83
derivatives is of the form 0 (A t5) in time and 0(Ax5) in space, the number of the
electromagnetic nodes involved in the computation o f each field component is larger than that
in the case of a second order scheme (2x2). Figure 3.1 shows the difference between a (2x2)
scheme, and a (4x4) scheme in the two dimensional space. The fields on the boundary must
be computed with special boundary equations. Figures 3.2 shows the different regions of
computation. The fields at the comers are computed by using the first order Mur’s AB C while
using the second order at the boundaries. The inner nodes are computed either by the (2x2)
FDTD scheme or with the (4x4) FDTD scheme according to the required accuracy. In the
(4x4) order case the (2x2) order is used only in the computation at the first nodes next to the
boundary, but this problem is solved by using FDA at the left and top boundaries and BDA
at the right and bottom boundaries as discussed in Chapter 2.
3.4 Dispersion and Stability Analysis
In spite of the fact that FDTD algorithm is a simple and efficient means o f solving
Maxwell’s equations for a wide variety of problems, numerical dispersion and anisotropy are
inherent to FDTD methods. The study of the total field formulations indicates that numerical
noise in the computation results is often caused by numerical dispersion which can be reduced
by suitable treatments. Phase errors from dispersion and grid anisotropy are significant (more
than rt/8 ) in large computational domains unless a small spatial discretization is used [35],
For large scattering problems (forty wavelengths or larger), a higher-order method is used
to decrease the expense of FDTD simulations. The (4x4) order method has less grid
anisotropy and dispersion as compared with the Yee algorithm (YA). This means the higher
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
84
Z
x
'2)
-►X
E
?
H r (I,
(U)
f
l+1/2) , *
y,
Y
Hy (1. H-1/2 )
-UI.F-1'2 )
Fig. 3.1a Spatial depositions of the electromagnetic nodes involved in the computation of
Ez(iJ) component in two-dimensional space for the (2x2) scheme.
H ,(i. j-3/2)
----►
H ,(l, J-1/2)
E 0.D
N/
/
H,(l, !-3 /2)
H,(t, i+1/2)
h , c , K 1/2)
Hr 0, J+3/2)
H, (i, K 1/2)
H ,(l.J+ 3/2)
Fig. 3.1b Spatial depositions of the electromagnetic nodes involved in the computation of
Ez(ij) component in two-dimensional space for the (4x4) scheme.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
85
1 “ order equation is used to get this element
Z
¥
► x
---------------------------------V
2 "d order region
Y
interface between 2 ' d and
4 in order regions
H, points
Ej points
X
4 “ order region
points
Fig. 3.2a FDTD domain and boundary in the TM* polarization case for the computation of
Ez(ij).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
86
xV
2 “d order region
Y
interface b etw een 2 "d and
X
4 *" order region s
H , points
Ej points
X
4 111 order region
Fig. 3.2b FDTD domain and boundary in the TM 2 polarization case for the computation of
HWU).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
87
X
T
4 l" order region
2°a
order region
X
X
E j points
X
X
Hr points
interface between 2 “<l and
4 w order regions
Fig. 3.2c FDTD domain and boundary in the TM* polarization case for the computation of
w
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
88
order FDTD scheme provides more accurate simulation than the YA as illustrated in the test
cases examined in this dissertation. The higher order FDTD allows larger spatial cell size
which in turn enables faster computation time and decreases memory size. In this chapter,
(2x2) and (4x4) order FDTD dispersion relations are derived, examined and compared.
3.4.1 Numerical Dispersion
Dispersion is defined as the variation of a propagating wave’s wavelength X with
frequency/. For convenience, dispersion is also frequently represented as the variation of the
propagating wave’s wavenumber k=27t/A with angular frequency ay=2itf.
Consider a lossless medium and assume the following monochromatic traveling-wave
trial solution [29,36], then
r
"
J Oj-k)
ff
_
p
J\ (t» n At - V A
-
- k^^z)
(3-7)
jo
2
- jj
k (,TI ,rT*"T)
eh (" " ^ ■^ A - V A - VA)
(3-8)
where kx,k y, and k. are the x, y, and z-components o f the numerical wavevector, and i„ iy,
and i. are space indexes, and j , =V-T. The subscript j or k of the field component to indicate
it may be x or y or z-component. By substituting (3-7) and (3-8) into (2-16) to (2-23) for
a lossless medium, one can obtain the numerical dispersion relation o f the fourth order
scheme [29]:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where:
kA x
X = px sin ( ^ P )
2 . r——
2kA y
(px-l)sin: k£ * + pusu
2
2
1 -1
6
/
1-
Y = p, sin ( ^ )
i
6
kAz
Z = p. sin (------)
1-
(3-10)
(Pv_l)sin:
2 . 2k£ *
p .s in -----
(3-11)
2 • 2*«AX
2 . 2kf i y
pi sm2— - + pvsur-^—
(p .-l)sin 2——
(3-12)
2 • 2^
PxSm ~ I —
cLt
p, - — , p, -
,
where
prsin
cAt
— . p, -
5
cAt
—
(3. 13)
For simplicity, consider the two dimensional space case for the TM2 mode. Let k,=0, Ax =
Ay = h,
p
= p
1
"
= 0.5
= p =
k
27tc
.
2
, 2 7 t .
1
oj= —— , k = — cos a , k = — sm a , PPW= — , c=*
x K
y K
h
fiX
(3.14)
tz
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3- 15)
90
where PPW 'xs the number of points per wavelength discretization, nris the angle between the
direction of the propagated wave and the positive x-axis and c is the speed of light. By
substituting (3-10) to (3-15) into (3-9), one can derive the expression:
cp _ X" _ ^
C
2 . 2/ 7tCOSa
p sin (
PPW
7 tp
X
2 • 2/tcsina.
+ p sin (------- )
PPW
(4
i-i
6
/ 2 n • 2/^tCOSffx
2 • 2/
(p
-l)su r(
) +p'£
sur(-------- )
H
PPW
PPW
2 . 2/7tcosa.
2
• 2/^tsma.
p2sin2(
) +(p -l)sin (--------)
K
v PPW
PPW
(3-16)
where cp is the numerical phase velocity. Equation (3-16) gives the ratio of the velocities or
the wavelengths as function of PPW and a for fourth order algorithm. In the second order
FDTD scheme case, one can follow the same steps but equations (2-34), (2-35), and (2-36)
are used instead o f equations (2-16) to (2-23), One obtains:
cp
-
Xn
=±
ppw
Ttp
sin
. 2,rccosa.
. 2/TtsmaN i t
su r(
) + sur(
)
ppw
ppw
(3-17)
Figure 3.3 shows the variation of the numerical phase velocity with points per
wavelengths discretization (PPW) in a two-dimensional FDTD grid for both (2x2) and (4x4)
order schemes. Another way to test dispersion is to compute the phase error due to
dispersion. The phase error for the plane wave propagation for h distance is given by [29]:
6err = 360°
h_
X
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3-18)
normalized phase velocity Cp/C0
91
0.99
—
(2x2) order
—
(4x4) order
0.98
0.97
0.96
0.95
0.94
10
20
30
40
Lo
ppw
Variation of the normalized phase velocity with PPW discretization for a plane
wave traveling at 0° using (2x2) and (4x4) order FDTD schemes.
Phase velocity error
4
3
—
(2x2) order
—
(4x4) order
2
1
0
10
20
30
40
ppw
Fig. 3.4 Variation of the phase velocity error with PPW discretization for a plane wave
traveling at 0° using (2x2) and (4x4) methods.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
92
The phase error is plotted in Fig. 3.4. In Figs. 3.3 and 3.4, the plane wave traveling at angle
a=0° with respect to the positive x-axis. Both (2x2) and (4x4) order scheme curves are
always less than the ideal case which is unity as illustrated in Fig. 3.3 . The curves for both
orders increase rapidly for PPW less than 10. This means that higher frequencies causing
broadening and distortion to the traveling wave than that at lower frequencies. Furthermore,
the amount of dispersion depends on the direction of propagation. Figure 3.4 shows the
variation o f grid anisotropy for the (2x2) and (4x4) order schemes. In this figure, the phase
error is smaller in forth order than the second order specially at low frequencies. The increase
of the PPW will decrease the amount of dispersion and grid anisotropy, but this will increase
the size of the required memory and the cost of the overall computation. The fourth order
FDTD scheme is suitable to solve large problems more efficiently than the second order
algorithm since it has lower dispersion and grid anisotropy at moderate values o f PPW (10
to 20 PPW).
Figure 3.5 graphs the numerical phase velocity versus wave propagation angle using
the YA and (4x4) algorithm with different values of PPW (5, 10, and 20). Figure 3.6 also
indicates the amount of grid anisotropy for the YA and (4x4) order scheme versus wave
propagation angle with different values of PPW (5 and 20). It is clear that the phase error is
smaller in forth order than the second order especially at low frequencies. The numerical
dispersion is decreased as the electric conductivity (losses) of the medium increases for the
second order method [39], and one can conclude that the same thing will happen for the
fourth order algorithm.
As a result, FDTD schemes are inherently dispersive and anisotropic. The resource
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
93
o
o
0.99
oo
o>
ato>
a
0.98
0.97
(2x2) order w hen PPW =5
0.96
(4x4) o rd er w hen PPW =5
(2x2) o rd er w hen PPW =10
0.95
CO
(4x4) w hen PPW =10
E
w
o
e
0.94
20
10
0
30
40
50
60
w ave a n g le .a (degree)
70
80
90
Fig. 3.5 Variation o f the normalized phase velocity with wave propagation angle (a)
using (2x2) and (4x4) order FDTD schemes.
o
3.5
(2x2) o rd er w hen PPW =5
(4x4) order w hen PPW =5
o
o
<d
ac>
CoO
2.5
•C
a.
1.5
0
5
10
15
20
25
30
w ave a n g le .a (degree)
35
40
45
Fig. 3.6 Variation o f the phase velocity error with a wave propagation angle (a) using
(2x2) and (4x4) methods.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
94
requirements are reduced by using the fourth order scheme approximations for the spatial
derivatives in Maxwell’s curl equations. The higher order scheme reduces both numerical
dispersion and grid anisotropy and uses a smaller number of points per wavelength and thus
reduce the size o f the required memory and decrease the number of operations needed.
3.4.2 Stability Analysis
The stability condition for an isotropic, homogeneous and lossless medium for second
order FDTD is investigated in one, two, three dimensional space (Courant stability condition)
[29], In three dimensional space the stability condition is given by:
c Ar s(3-19)
(Ax) 2
(Ay) 2
(Az)2
when Ax = Ay = Az = h, equation (3-19) can be simplified to:
while for the two dimension case, it reduces to:
and for the one dimension case, it is:
L
c A/ ^ —
v/3
(3 -20)
’
c A/ <;— ,
v/2
(3-21)
c A/ sh
(3-22)
On the other hand, the stability criterion for a lossy dielectric is studied in the two dimensional
space using the second order accuracy [39], It is given by:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
95
-i
(3-23)
where yr is given by:
y = —
r
c
1
_
2
1 + ( — )2
\\
cue
2
-l
(3-24)
)
but the value of cosh function is always not less than 1 , this means that the stability in lossy
material is less than the lossless one. The stability condition for the fourth order FDTD in the
lossless medium is the same as the second order FDTD scheme [29].
3.5 Numerical Verification
A computer program based on FDTD has been developed for this part of the
investigation. The program directly solves Maxwell’s equations via finite differences. This is
done for the second and fourth order schemes in free space. The solution is for two
dimensions that corresponds to a line current source in the z-direction at certain point in the
domain (source point). The field is computed at certain point in the medium (observation
point) as a function of time. The program is also used to determine the field variation along
the x-axis at specified time steps.
The stability and accuracy constraints on Ax, Ay, and At must be taken into
consideration. For good accuracy, Ax and Ay must be much less than the minimum required
wavelength (usually sA/10). For stability, a time step must be small enough so that the field
values can affect only the nearest neighbors during one time step interval. The necessary
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
96
stability condition for the two-dimensional case when the cell size in x and y-directions are
equal to h is discussed in equation (3-21).
The excitation of the difference equations can be done by using different shapes of
pulses that depends on frequency content. A current source with a Gaussian pulse waveform
is chosen, and it is given by
( i -
= Jo e
-
m Ax
T T
)2
7
(3-25)
.
’
'» = n r
where Ja is the peak amplitude. Ja , m, and n are chosen to be 5 ampere,
P-26)
10,
and
2,
respectively. These parameters (m and n) are discussed before in Chapter two. Equation
(3-25) can be written in a difference form:
( n At - r )
(3-27)
Z
O
This pulse is added to the expression of the electric field in the z-direction.
The excitation pulse must be absorbed at the limits of the problem space. Mur’s
absorbing boundary condition (ABC) is used after reducing it to two dimensional
computational domain [3], In the case the boundary condition is of first order at the comers,
while it is second order at edges between the comers. For the (4x4) algorithm, E. is
computed by its second order accuracy at the outer boundaries and fourth order accuracy in
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
97
the interior region.
Figure 3.7 shows Ez as a function of time at the following observation points
(500,500), (600,500), (700,500), and (800,500) when the source is located at (500,500). The
field is computed using second and fourth order FDTD schemes after 28.2 ns (This time
corresponds to 300 time steps where each step is 94.34 ps). The two dimensional space is
square in area and contains 106 cells of free space medium. The spatial discretization in x and
y directions (Ax and Ay) are the same and of length 4 cm. It is clear from Fig. 3.7 that the
response o f the fourth order FDTD scheme has no distortion (dispersion) but the result using
second order scheme has very high level of ripples caused by dispersion. This means that the
fourth order (4x4) is of higher accuracy than YA (2x2). One can see also from the plotted
response that as the distance between the observation point and the source increases, the
amplitude of the field decreases but the level of the ripples increases relating to the maximum
value of the field.
The shapes of a Gaussian pulse propagated in free space at several numbers of time
steps are shown in Fig. 3.8. In this figure, Ez is plotted versus Ax position at time steps 50,
200, 400, 600, 800, and 1100 using the same parameters mentioned for Fig. 3 .7. These time
steps correspond to 4.7, 18.8, 37.6, 56.4, 75.2, 103.4 ns, respectively. It is clear from this
figure that the field computed by the fourth order scheme has no distortion but the second
order scheme gives high levels o f distortion. One can see also that, as the distance traveled
by the pulse increases, its amplitude decreases and also the ripple level increases compared
with the maximum amplitude o f the pulse. It is obvious from Fig. 3.8f that the pulse is
reflected at the boundary with its level becoming very small and the ripples appear in the right
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
98
15
_
10
5
0
lu
5
10
1
0
3
2
4
5
time(ns)
(a)
1. 2
2 ,d order
4 * order
0.8
E 0.4
0
lu
-0.4
-
0.8
16
18
20
24
22
time(ns)
26
28
(b)
Fig. 3.7
Ez versus time calculated in a domain size o f 1000x1000 ceils with a side equal to
4cm using second and fourth order FDTD schemes when a a current source
with a Gaussian pulse waveform is applied at the point (500,500) in the
z-direction. Ez is computed at the points: (a) (500,500), (b) (600,500), (c)
(700,500), and (d) (800,500).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
99
2 “ order
0.8
4 “ order
0.4
e
-0.4
-
0.8
30
35
40
time(ns)
45
50
(c)
0.8
0.6
"e 0.4
—
2
—
4 “ order
order
uj* -0 .2
-0.4
0.6
- 0.8
-
50
52
54
56
58 60 62
time(ns)
64
66
68
70
(d)
Fig. 3.7
Ez versus time calculated in a domain size of 1000x1000 cells with a side equal to
4cm using second and fourth order FDTD schemes when a a current source
with a Gaussian pulse waveform is applied at the point (500,500) in the
z-direction. Ez is computed at the points: (a) (500,500), (b) (600,500), (c)
(700,500), and (d) (800,500).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
100
2
—
2 '" o rd er
—
4 “ order
1
0
-1
-2
500
510
530
540
600
610
520
position
(a)
2
0.8
order
4 “ order
E 0.4
>
N
LU
-0.4
-
0.8
570
580
590
position
(b)
Fig. 3.8 The shapes o f a Gaussian pulse propagated in free space using the same parameters
mentioned in Fig. 3.7. E z versus delta-x position after various numbers o f time
steps: (a) n=50 (4.7 ns), (b) n=200 (18.8 ns), (c) n=400 (37.6 ns), (d) n=600 (56.4
ns), (e) n=800 (75.2 ns), and (f) n= l 100 (103.4 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
10 1
0.8
2
order
4
order
E* 0.4
LU
-0.4
-
0.8
660
670
680
690
position
700
710
(c)
1.2
0.8
4
order
0.4
0
-0.4
-
0.8
740
760
780
position
800
820
(d)
Fig. 3 .8 The shapes of a Gaussian pulse propagated in free space using the same parameters
mentioned in Fig. 3.7. Ez versus delta-x position after various numbers of time
steps: (a) n=50 (4.7 ns), (b) n=200 (18.8 ns), (c) n=400 (37.6 ns), (d) n=600 (56.4
ns), (e) n=800 (75.2 ns), and (f) n= 1100 (103.4 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
10 2
2 “ order
4 * order
0.8
0.4
E
>
M
LU
-0.4
-
0.8
840
860
880
900
position
920
940
(e)
0.08
—
—
0.06
_
0.04
|
0.02
-
2 “ order
4 Ik order
0.02
-0.04
-0.06
920
940
960
position
980
1000
(0
Fig. 3 .8 The shapes o f a Gaussian pulse propagated in free space using the same parameters
mentioned in Fig. 3.7. Ez versus delta-x position after various numbers o f time
steps: (a) n=50 (4.7 ns), (b) n=200 (18.8 ns), (c) n=400 (37.6 ns), (d) n=600 (56.4
ns), (e) n=800 (75.2 ns), and (f) n=1100 (103.4 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
103
side of the pulse instead of the left side.
In this part o f this section, there are two cases that will be discussed. The first case
is to compare the results of both second and fourth order methods when using an excitation
current source with a Gaussian pulse waveform or with a Gaussian pulse modulated by a
sinusoidal wave form. These waveforms are given respectively by:
_ ( f - 3f„ ) 2
Jz = J 0 e
Jz =
where
'•
- ( ' ' 3/° }2
J Q sin[27t/e(r -3 g ] e
f°
(3-28)
(3 29)
J = l6 0 0 A /m 2, t - ----- ^------, / =300 MHz, f .=450 MHz
°
(fh- f c) Jc
The signature of these sources is transformed from the time domain to the frequency domain
and sketched in Fig. 3.9. The maximum frequency of the current source with a Gaussian
pulse waveform in this case is about 0.3 GHz. Figure 3.10a shows Ez as function of time at
the observation point (86,46) when using the sources given by equation (3-28) at (106,106).
The field is computed using the second and the fourth order FDTD algorithms after 20.85 ns
(this time corresponds to 500 time steps where each step is 41.7 ps). The two dimensional
space is a square area of 211x211 cells. The medium is free space. The spatial discretization
in x and y directions (Ax and Ay) are the same and o f length 2.5 cm. It is clear from Fig.
3.10a that the response of both the YA and (4x4) order schemes are coincident. This
numerical experiment is repeated with the same parameters but using the source given by
equation (3.29) and the results are shown in Fig. 3.10b which illustrate no difference between
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
104
—
e-«'-3W ’
J Z= J 0 sin [2 n fc(t-3t0)] e '((,'3to ,/,o)‘
0.8
0 .4
1.2
1.6
F(G H z)
Fig. 3.9 The current source with a Gaussian
pulse waveform and with a Gaussian pulse
modulated with a sinusoidal waveform in the frequency domain.
60
e
—
fo u rth
4°
>20
^
second
—
exact
0
LU
-20
-4 0
8
10
12
14
16
18
20
tim efn s)
(a)
—
second
fo u rth
exact
-1 5 0
6
8
10
12
14
16
18
20
tim e(n s)
Fig. 3.10 Ez versus time calculated in a domain size of 211x211 cells with a side equal to
2.5 cm using second and fourth order FDTD schemes when the source is
applied at the point (106,106) in the z-direction. Ez is computed at the point
(86,46) using (a) a current source with a Gaussian pulse waveform, (b)
a current source with a Gaussian pulse modulated with a sinusoidal waveform.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
105
the two algorithms. The responses due to the two sources are compared with the exact
solution given in [40], The exact solution in this case for Ez in free space at the near zone due
to a line source (JJ is given by [40]:
',(P,
271
Jft-s)
^ r —
^
y]s+2x
['
Jo
(3-30)
_ Ip-pl
s is the time step, x=
where
p and p are the vectors from the origin to the field point and to the source point respectively.
The integral in equation (3-30) is solved by using a three point weighted Newton-Cotes
formula which is given by:
/ b p r ds = wl J[a) * w2 f
where
' a+b'
(3-31)
J ( t s)
nk t
fiu ) = -■
w=— «= 1, 2, 3,
\fs+2x
2
2
w .=—
1 15
.1400
(3-32)
2(5n+36)2 ~(a2+l5b 2)fa
( b - a )2
2
16 (5b-a)a 2 -(5 a-b)b
w2=15
(6 - a ) 2
and
w3=.
(15a2+Z>2)y^-2(3a-i-53)a
15
2
(6 - a ) 2
(3-33)
The next case is to compare the results of both second and fourth order methods
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
106
when using an excitation current source with a Gaussian pulse waveform or with a Gaussian
pulse modulated by a sinusoidal wave form which are given respectively by:
(3-34)
T 1
(3-35)
where
J = 125 A /m 2,
T=
10 Ax,
2c
t =27
Both sources are plotted in the frequency domain as shown in Fig. 3.11. The maximum
frequency of the current source with a Gaussian pulse waveform in this case is about 1.1
GHz. Figure 3.11a shows Ez as function o f time at different test points when using the
source given by equation (3-34) at (500,500). The field is computed using the second and the
fourth order FDTD algorithms after 93.38 ns (this time corresponds to 1400 time steps
where each step is 66.7 ps). The two dimensional space is a square o f area 1000x1000 cells.
The medium is free space. The spatial discretization in x and y directions (Ax and Ay) are the
same and are of length 4 cm. It is clear from Fig. 3.12 that there is a difference between the
response of both the YA and (4x4) order schemes. The fourth order method has no distortion
and coincident with the exact solution. On the other hand the second order (2x2) has high
level of distortion and there is a big shift in amplitude and phase compared to the exact
solution This numerical experiment is repeated with the same parameters but with the source
given by equation (3.35) and the results are shown in Fig. 3.13. In this case the higher order
method also gives a coincident response with the exact solution while the second order has
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
107
20
JZ=J0 e'({Ho,/Tmax)l
JZ=J0 sin[27i fc(t-t0) ] e ‘((t ,o,/Tm
0
0.4
0.8
1.2
1.6
2
F(GHz)
Fig. 3.11 The current source with a Gaussian a pulse waveform and with a Gaussian pulse
modulated with a sinusoidal waveform in the frequency domain.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
108
—
(2x2) order
E
(4x4) order
>
— • exact
uT
-20
-40
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
40
—
—
(2x2) order
(4x4) order
20
— •
exact
25
time (ns)
(a)
80
60
0
-20
-40
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
time (ns)
(b)
Fig. 3.12 Ezversus time calculated in a domain size of 1000x1000 cells with a side equal to
4cm using second and fourth order FDTD schemes when the excitation is a
current source with a Gaussian pulse waveform applied at the point (500,500)
in the z-direction. Ez is computed at the following points: (a) (600,500). (b)
(700,500). (c) (800,500). (d) (900,500).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
109
E,(l) (V/m)
—
(2x2) order
—
(4x4) order
—
exact
-20
-40
38
39
40
41
42
43
44
45
46
47
48
time (ns)
(c)
80
E.O) (Wm)
60
—
40
second
-*• fourth
20
-•
exact
0
-20
-40
50
52
54
56
58
60
62
time(ns)
W)
Fig. 3.12 Ez versus time calculated in a domain size o f 1000x1000 cells with a side equal to
4cm using second and fourth order FDTD schemes when the excitation is a
current source with a Gaussian pulse waveform applied at the point (500,500)
in the z-direction. Ez is computed at the following points: (a) (600,500). (b)
(700,500). (c) (800,500). (d) (900,500).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
11 0
(2x2) order
E
(4x4) order
>
exact
ilT
-20
-40
12
13
14
15
16
17
18
time (ns)
(a)
80
60
40
20
—
(2x2) order
—
(4x4) order
—
exact
0
-20
-40
24
25
26
27
28
29
30
31
32
33
34
time (ns)
(b)
Fig. 3.13 Ez versus time calculated in a domain size of 1000x1000 cells with a side equal to
4cm using second and fourth order FDTD schemes when the excitation is a
current source with a Gaussian pulse modulated with a sinusoidal waveform
applied at the point (500,500) in the z-direction. Ez is computed at the following
points: (a) (600,500). (b) (700,500). (c) (800,500). (d) (900,500).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Ill
80
60
(2x2) order
f
40
(4x4) order
>
B
uT
ex act
20
0
-20
38
39
40
41
42
43
44
45
46
47
time (ns)
(c)
80
60
— (2x2) order
I"
>
40
3
20
-• (4x4) order
- ex a c t
us
0
-20
50
51
52
53
54
55
56
57
58
59
60
61
62
time (ns)
(d )
Fig. 3.13 Ej versus time calculated in a domain size of 1000x1000 cells with a side equal to
4cm using second and fourth order FDTD schemes when the excitation is a
current source with a Gaussian pulse modulated with a sinusoidal waveform
applied at the point (500,500) in the z-direction. Ez is computed at the following
points: (a) (600,500). (b) (700,500). (c) (800,500). (d) (900,500).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
112
high level of ripples and far from the exact solution.
One can conclude that, when the frequency o f the excitation source is low (0.3 GHz)
there is no significant difference in accuracy between the YA and (4x4) algorithms, and
there is a noticeable difference when the increase o f frequency.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
113
Chapter 4
TOTAL FIELD ANALYSIS IN ONE DIMENSIONAL SPACE
USING HIGHER ORDER FDTD TECHNIQUES
The formulation of the total fields for fourth and second order accuracy is derived in
this chapter. The accuracy of the second and fourth order FDTD schemes are analyzed and
compared through numerical examples of plane wave propagating in free space and plane
wave reflection and transmission from a dielectric slab.
4.1 Total Field Formulation of Fourth Order FDTD Scheme
In this chapter the TEzmode is considered to reduce the complexity of the derivation
and programming of the total field formulation. This means that H. and Ey components are
considered only in the one dimensional problems. The total field expressions derived in the
three-dimensional space can be used to get the corresponding formulations in one dimensional
space by canceling H„ Hy, Ex, and E 2 and neglecting the terms related to the derivatives with
respect to z and^. Using equations (2-17) and (2-23), one can get:
n 1
V (i)
Jk E n(0
q,
1
n»—
2
24 ri 1er
_
-3 H
H
* («’|)
2
1
2
1
r.-r—
2
+ 3 H
- H
r(.-{)
2 (- 2 }
2}
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(c A f)V eowpx
1
H
2
H
2
24r) l ^ r
where
1
n*—
n*—
('-j)
- H
o At
2
e
r C"T)
r-lt- 1
24etij
L y{<)
and t ]2 are defined by equations (2-17) and (2-18).
\
= h .n
j~n
Px
Px
\
_ i
jr n
>(/-!) +
3£'
3 F"
E y" 0 ’ 2)
24tl£1erpr
1
J
p A/
o ffl \
/ o£
( — +
)
24q ^,pr
e
( c A/)2aeqmpx
r n
y0~i)
24
a_A/
24p^j
y
/
(
3
Ey (0
Ey 0 - 1)
"
0*1) + J3 E >" (») - ^E O
- l)
r-n*!i
(,*1)
0
" (/-l) - F"
>' (0
-( £ •yvV
/)
- TTn
y
(0
11?iPrer
i
n- —
2
(»-j)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
115
where ^ and
e = €„ fi =
$2
are defined by equations (2-22) and (2-23). In this case, one can substitute
and a, = am into equations (4-1) and (4-2) to get the total
field difference equations of fourth order scheme (4x4) in free space.
4.2 Total Field Difference Equations of Second O rder FDTD Scheme
By taking into consideration the equations derived in section 4.1, based on the second
order accurate formulation, one can obtain:
1
. /I ♦ I
£ Xfl “
E"
(e+o A t) ry ( 0
1
n -—
H
\
n+—
-H
\
(4-3)
1
1
/! +—
H
nepPx
(e+o A/)
rt- —
2
(n +o Ar)
H
2
r*n
c (p + o A 0
- y ( i-l) “
r* n
(/)
(4-4)
One can substitute e = eol [i = ji0, and at = a m=0 in equations of this section to
derive the total field difference equations of second order scheme (2 x2 ) in free space.
4.3 Node D istribution for the Second and Fourth O rder Schemes
As discussed before in Chapter Two and Three, the number of the electromagnetic
nodes involved in the computation of each component of the field in the fourth order scheme
is larger than in the second order scheme. Figure 4 .1 shows the difference between (2x2) and
(4x4) schemes in one dimensional space for the electric and magnetic field components.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
116
t
Hz(i-1/2)
A
Hz(i+1/2)
_
A
Ey(0
(a)
t
Hz(i-3/2)
A
Hz(i-1/2)
Hz(i+1/2)
A
A
Hz(i+3/2)
A
► x
x
Ey(i)
(b)
Fig. 4.1 Spatial depositions of the electromagnetic nodes involved in the computation of
Ey(i) component in one dimensional space with the propagation in the x-direction
for: (a) (2x2) order scheme, (b) (4x4) order scheme.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
117
4.4 Applications
Using FDTD, Maxwell’s equations are solved directly in time domain via finite
differences and time stepping. This is done for second and fourth order schemes in free space
and lossless and lossy dielectric medium. Here, two applications are discussed and analyzed
through using the two FDTD algorithms; the solutions correspond to a normal incidence of
a plane wave propagating through free space and through a free space medium close to a
lossy or lossless dielectric slab.
The stability and accuracy constraints on Ax and At must be taken into consideration.
For good accuracy, Ax must be much less than the minimum required wavelength (*A/1 0 ).
For stability, time steps must be small enough so that the field values can affect only nearest
neighbors during one time step interval. The necessary stability condition for the one­
dimensional case is:
c A t <Ax
(4-5)
where c is the speed of light. A conservative choice is [41]:
c A / = 0.5 A x
(4‘6)
The excitation can be done by using different shapes of pulses that depend on the
frequency content. A Gaussian pulse electric field is chosen, and it is given by [41]:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
x - x_
( » - * . -
118
,
)
1
(4-7)
Eyl = E e
T =
with
,
1,
= . ?
(4-8)
2 C
where Ea is the peak amplitude, xp is original location o f the peak, T is the pulse width, and
m is a parameter used to change the pulse width. Equation (4-7) can be written in the
following difference form:
( „ Ar - fa -
) 1
------------7 7 - ^ -------
(4‘9)
E . = Eo*
The corresponding magnetic field of the incident wave is given by its difference form:
( / » 0.5 ) Ar - x
,
( " Ar - r„ - L
L ) *
ZI
= HO
oe
Tl
(4"10)
where Ha is the peak amplitude o f the magnetic field.
The excitation pulse must be absorbed at the limits of the problem space. Mur’s
absorbing boundary condition (ABC) is used after reducing it to a one dimensional normal
incidence problem. In this case the boundary condition is:
XT' r +i
E y i =1
T7 hc A r
-
E y i- 2
A x ✓ t— /*
c
^
+
Ay
-X/~2 "
" ^7^ V
^
E y i= 2
n»—
n a
\
E y ,= l
)
(4 -1 1 )
ni
" z /.l J
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.12)
119
At
( 2* "
Ax
- E"
)
c At - Ax
(\ Ey "f=.
Ey "*v = Ey "r=//x-l
v i +
c At + Ax
1
n-—
(4-13)
(4-14)
__1_
2
H,r »=JV -1
= Hz ] 2 . i -
A/
li Ax
- E v"/=.
(4-15)
4.4.1 Free Space
In this case, the used normal incident plane wave Gaussian pulse is of amplitude
lV/m, the original location of the peak (xp) at 170, and m and n parameters are 10 and 2,
respectively. The medium is free space with 3000 cells in the x-direction. The cell size (Ax)
is 5 mm. The equations corresponding to the free space region for the (4x4) and (2x2)
schemes as discussed in sections 4.1 and 4.2 are programmed. Figure 4.2 shows a sequence
of plots of Ey as function of position x at different time steps. The exact solution, second and
fourth order FDTD schemes are illustrated in this figure. In this figure the fourth order (4x4)
FDTD scheme response is coincident with the exact solution, and, on the other hand, the
second order (2x2) FDTD scheme solution has distortion which increases as the time step
increases. This means that the fourth order algorithm gives a response with a very good
accuracy compared to that of the YA.
4.4.2 Lossy and Lossless Dielectric Structures
In this case, the normal incident plane wave Gaussian pulse used in the free space
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
120
i
- exact
0.8
- 2 “ order
- 4 * order
0.6
uT
2
order
4 * o rd e r
'
0.4
0.2
0
210
200
220
230
380
240
390
400
410
44C
1
1
—
in
430
(b)
(a)
E
420
Position
position
exact
— exact
— - 2 “ order
0.5
4 “ order
0.5
in
0
—
4 * order
0
-0.5
860
880
900
Position
(c)
920
940
1100
1120
1140
1160
1180
1200
Position
(d)
Fig. 4.2 The shapes of a Gaussian pulse propagated in free space at different time steps (a)
100 time steps ( time is 0.034 ns), (b) 500 time steps ( time is 0.169 ns), (c) 1500
time steps ( time is 0.508 ns), (d) 2000 time steps ( time is 0.678 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
121
1
—
0.8
2“ order
4“ order
0.6
0.4
0.2
0
200
240
230
220
210
position
(a)
0.6
2“ order
4th order
0.4
0.2
0
-
0.2
-0.4
210
220
230
240
250
260
270
position
(b)
Fig. 4.3 Interaction of a Gaussian pulse plane wave with a lossy dielectric slab, with relative
permitivity 4.0, 30 cm thick, located between delta-x positions 250 and 309,
electric conductivity 0.05, and magnetic conductivity 1.5 calculated using FDTD.
Plots are electric field versus delta-x position through different cells after
different time steps. Each time step is 8.33 ps. (a) 100 time steps ( time is 0.833
ns), (b) 200 time steps ( time is 1.666 ns), (c) 400 time steps ( time is 3.332 ns),
(d) 625 time steps (time is 5.21 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
122
0.2
0.1
— 2 " * order
0
- border
-0.1
-0.2
-0.3
-0.4
100
150
200
250
300
350
position
(c)
0.2
— 2“ order
--• 4thorder
0.1
0
-
0.1
0.2
-0.3
-0.4
0
50
100
150
200
250
position
300
350
400
450
(d)
Fig. 4.3 Interaction of a Gaussian pulse plane wave with a lossy dielectric slab, with relative
permitivity 4.0, 30 cm thick, located between delta-x positions 250 and 309,
electric conductivity 0.05, and magnetic conductivity 1.5 calculated using FDTD.
Plots are electric field versus delta-x position through different cells after
different time steps. Each time step is 8.33 ps. (a) 100 time steps ( time is 0.833
ns), (b) 200 time steps ( time is 1.666 ns), (c) 400 time steps ( time is 3.332 ns),
(d) 625 time steps (time is 5.21 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
123
application is applied with the same amplitude (lV/m). Figure 4.3 shows a sequence of plots
of Ey as function of position x using (2x2) and (4x4) schemes for different time steps. The
geometry has been made to be a one dimensional stratified medium with 600 layers. The
medium is free space except layers 250-309 that consists o f a lossy dielectric slab with
relative permitivity of 4, electric conductivity of 0.05 S/m, and magnetic conductivity of 1.5
S/m. The width of each layer (Ax) is 5 mm. In this case the thickness of the dielectric slab is
30 cm. The distortion in the response of the second order FDTD scheme shown in Fig. 4.3
increases as the time step increases, but the fourth order FDTD scheme solution has no
distortion. In this case, the parameters Xp, m and n of the Gaussian pulse are 170, 10, and 2,
respectively. In order to see the distortion of the second order accuracy clearly, a single case
of Fig. 4.3, is plotted with its two parts ( Incident and reflected pulses). This case is taken
at the 400 time step and Ey computed by second and fourth order schemes as shown in Fig.
4.4. Ey is plotted in Figs. 4.5, 4.6, and 4.7 for the same geometry considered in Fig.4.3, but
with different cell size, 1.5, 3, and
6
cm, respectively, and a lossless dielectric slab. From
these figures, one can conclude that, as the spatial discretization increases the size of both the
domain of computation and the required memory decreases and the distortion from the
second order FDTD scheme increases. The fourth order FDTD scheme has no distortion at
all in Figs. 4.5 and 4.6, and very low distortion begins to appear in Fig. 4.7 when the cell size
is 6 cm. The parameters used for the domain, cell size, and Gaussian pulse are different in
Figs. 4.5, 4.6, and 4.7 as illustrated in each figure.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
124
E
I
—
(2x2) o rder
E
—
(4x4) o rder
o
- 0.1
^
>»
-
Ul
0.2
-0.3
-0.4
100
200
150
300
250
350
position
(a)
0
(2x2) o rder
-
0.1
(4x4) order
0.2
-0.3
-0.4
100
120
110
130
150
140
160
position
(b)
0.2
—
(2x2) o rder
—
(4x4) o rd er
E
>
-
0.1
290
300
310
320
33 0
position
(c)
Fig. 4.4 Interaction of a Gaussian pulse plane wave with a lossy dielectric slab, with relative
permitivity 4.0, 30 cm thick, located between delta-x positions 250 and 309,
electric conductivity 0.05, and magnetic conductivity 1.5 calculated using FDTD.
Plots are electric field versus delta-x position through different cells after
400 time steps ( time is 3.332 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
125
1
0.8
^
0.6
^
0.4
Ul
0.2
0
- 0.2
2
order
4 * order
100
200
300
position
400
500
600
400
500
600
400
500
600
(a)
0.8
0.6
0.4
2 order
4 * order
0.2
0
- 0.2
ui •
-0.4
100
200
300
position
(b)
1
„ 0.8
E 0.6
>
0.4
2 ,d order
4 * order
0.2
yv
■0
-0.2
-0.4
100
200
300
position
(c)
Fig. 4.5 Interaction of a Gaussian pulse plane wave with a lossless dielectric slab, with
relative permitivity 4.0, 9 cm thick, located between delta-x positions 250
and 309 calculated using FDTD. Plots are electric field versus delta-x position
through different cells after different time steps. Each time step is 2.5 ps, cell size
is 1.5 cm, and xp is 170. (a) 50 time steps ( time is 0.125 ns), (b) 300 time steps
(time is 0.75 ns), (c) 625 time steps ( time is 1.56 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
126
0.8
0.6
0.4
0.2
—
—
2 •"o rd er
4 * order
-0.2
50
100
150
position
200
250
300
200
250
300
200
250
300
(a)
0.8
s ' 0.6
> 0.4
0.2
>»
0
LU
0.2
0.4
2 *d order
4 order
I
"
V
50
100
150
position
(b)
0.8
0.6
0.4
—
—
2
4
order
order
0.2
LU
-
0.2
-0.4
50
100
150
position
(c)
Fig. 4.6 Interaction o f a Gaussian pulse plane wave with a lossless dielectricslab, with
relative permitivity 4.0, 9 cm thick,
located between delta-x positions 125
and 154 calculated using FDTD. Plots are electric field versus delta-x position
through different cells after different time steps. Each time step is 5 ps, cell size is
3cm, and xp is 85. (a) 50 time steps ( time is 0.25 ns), (b) 175 time step ( time is
0.875 ns), (c) 300 time steps ( time is 1.5 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
127
-
0.2
50
100
150
100
150
100
150
position
(a)
0.8
e
°-6
0.4
0.2
5
r:
LL)
—
—
b
0.2
-0.4
-
2 * order
4 * order
_A
:v
50
position
(b)
0.8
0.6
0.4
—
- -
2 ,d order
4 * order
0.2
50
position
(c)
Fig. 4.7 Interaction of a Gaussian pulse plane wave with a lossless dielectric slab, with
relative permitivity 4.0, 9 cm thick, located between delta-x positions 62 and 77
calculated using FDTD. Plots are electric field versus delta-x position through
different cells after different time steps. Each time step is 10 ps, cell size is 6 cm,
and xp is 42.5. (a) 12 time steps ( time is 0.12 ns), (b) 75 time step (time is 0.75
ns), (c) 150 time steps ( time is 1.5 ns).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
128
C hapter 5
CONCLUSIONS
In this research work, a three dimensional space formulation for the fourth order
FDTD algorithm in time and space (4x4) in lossy media has been derived, and the
formulations for (2x4) and (2x2) FDTD techniques have been deduced for different
dimensional spaces and various media. Transient and continuous wave simulations are
applied to compare the accuracy of the YA and the higher order methods employing the (2x4)
and (4x4) algorithms. The higher order technique (4x4) always gives better accuracy with
the freedom of reducing the size o f the required memory compared with (2 x2 ) order
technique. However, during this investigation, the fourth order FDTD scheme is found to
have the same accuracy as the YA when applied with structures such as microstrip
transmission lines and low pass filters. The expressions for the soft sources used to excite
microstrip circuits are modified to the fourth order accuracy; thus, improved results compared
to that obtained using a second order accuracy expression for the source are obtained.
In a three dimensional space, The resonant frequencies of a rectangular cavity
resonator have been computed with different cell discretization using (2x2) and (4x4) orders
FDTD techniques. The results show that (4x4) gives more accurate results with nearly the
same execution time and a reduced memory size by a factor of 37 compared to the (2x2). A
comparison between the three different orders o f accuracy was performed in free space,
dielectric media and magnetic media . The (4x4) order FDTD algorithm always gives more
accurate results for the propagating wave with no distortion in free space and a very low
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
129
level of distortion in dielectric and magnetic structures compared to that o f the other (2x4)
and (2x2) techniques. The (2x4) algorithm gives low level of ripples compared to that o f the
(4x4) algorithm when the geometry has interfaces such as a dielectric slab or a magnetic slab
in a free space medium. However the (4x4) scheme has better stability than the (2x4) scheme.
The presented fourth order (4x4) FDTD scheme has been tested in the two
dimensional space. The results of these tests showed that the fourth order FDTD scheme
exhibited much less numerical dispersion than that o f the YA (2x2) for high frequency
applications. The numerical dispersion was also investigated for different order methods. The
higher order algorithm greatly reduces the effects o f this type dispersion. Therefore the
resource requirements (memory size and computational cost) are reduced by applying the
higher order method.
In one dimensional space, the derived fourth order FDTD technique has been
examined in free space and dielectric media. Its performance is coincident with the exact
solution in free space while the Yee algorithm gives very high level of distortion. The (4x4)
and (2 x2 ) orders of accuracy are compared for the propagation of a plane wave through a
dielectric slab using various cell discretizations. The results show that the higher order
scheme (4x4) yields more accurate results for large space discretization, and hence reduces
the size of the required memory. Finally, it is concluded that the fourth order technique can
be used for applications where high accuracy is demanded.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
130
LIST OF REFERENCES
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
131
References
[1]
K . S . Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat., vol.
AP-14, pp. 302-307, May 1966.
[2]
K. Kunz and R. Luebbers, The Finite Difference Time Domain M ethod fo r
Electromagnetics. Boca Raton, Florida: CRC Press, 1993.
[3]
G. Mur, “Absorbing boundary conditions for the finite-difference approximations
of the time-domain electromagnetic field equations,” IEEE Trans. Electromagnetic
compatibility, vol. EMC-23, pp. 1073-1077, Nov. 1981.
[4]
W. Q. Sui, D. A. Christensen, and C. H. Dumey,"Extending the two-dimensional
FDTD method to hybrid electromagnetic systems with active and passive lumped
elements,” IEEE Trans. Microwave Theory and Tech., vol. MTT-40, April 1992.
[5]
R. J. Luebbers and H. S. Langdon, "Reducing the number of time steps needed for
FDTD antenna and microstrip calculation," Proceeding o f the 11* Annual Review
of Progress in Applied Computational Electromagnetics (Monterey, CA), pp. 465471, March 1995.
[6 ]
A. Cangellaris and R. Lee, “On the accuracy of numerical wave simulations based on
finite methods,” Journal o f Electromagnetic Waves and Applications, vol. 6 , pp.
1635-1653, Aug. 1992.
[7]
A. Bretones, R. Martin and I. Garcia, “Time-Domain analysis o f magnetic-coated wire
antennas,” IEEE Trans. Antennas Propagat., vol. 43, pp. 591-596, June 1995.
[8 ]
B. Popovic and A. Nesic, “Generalization of the concept o f equivalent radius of thin
cylindrical antennas,” IEEE Trans. Antennas Propagat., vol. 43, pp. 591-596, June
1995.
[9]
N. Alexopoulos, P. Katehi, and D. Rutledge, “Substrate optimization for integrated
circuit antennas,” IEEE Trans. M T T , vol. MTT-31, pp. 550-557, July 1983.
[10] I. Rana and N. Alexopoulos, “Current distribution and input impedance of printed
dipoles,” TF.FF Trans. Antennas Propagat., vol. AP-29, pp. 99-105, Jan. 1981.
[11] S. Amari, M. Gimersky, and J. Bomemanne, “Imaginary part of antenna’s admittance
from its real part using Bode’s integral,” IEEE Trans. Antennas Propagat., vol. 43,
pp. 220-223, Feb. 1995.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
132
[12] L. Felsen, Transient Electromagnetic Fields. Springer-Verlag Berlin Heidelberg New
York, 1976.
[13] B. Toland, J. Lin, B. Houshmand, and T. Itoh, ‘TDTD analysis of an active antenna,”
TF.F.F. Microwave and Guided Wave Letters, vol. 3, Nov. 1993.
[14]
W. Sui, D. Christensen, and C. Dumey, “Extending the two-dimensional FDTD
method to hybrid electromagnetic systems with active and passive lumped elements,”
IEEE Trans. MT T , vol. 40, pp. 724-730, April 1992.
[15]
A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic
scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans.
Antennas Propagat., vol. MTT-23, pp. 623-630, Aug. 1975.
[16]
R. Holland, L. Simpson and K. S. Kunz,“Finite-difference analysis o f EMP coupling
to lossy dielectric structures,” IEEE Trans. Electromag. Compat., vol. EMC-22,
pp. 203-209, Aug. 1980.
[17]
R. Holland, L. Simpson and K. S. Kunz, “Finite-difference analysis o f EMP coupling
to thin struts and wires,” IEEE Trans. Electromag. Compat., vol. EMC-23, pp.
88-97, May 1981.
[18]
K. S. Kunz and L. Simpson,“A Technique for increasing the resolution o f finitedifference solutions of the Maxwell equation,” IEEE Trans. Electromag. Compat.,
vol. EMC-23, pp. 419-422, Nov. 1981.
[19]
B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical
simulation of waves,” Math. Comp, vol. 31, pp. 629-651, 1977.
[20]
B. Engquist and A. Majda, “Radiation boundary conditions for acoustic and elastic
wave calculations,” Comm. Pure Appl. Math., vol. 32, pp. 313-357, 1979.
[21]
D. H. Choi and W. J. R. Hoefer, “The finite-difference time-domain method and its
application to eigenvalue problems,” IEEE Trans. M T T , vol. MTT-34, pp. 14641470, Dec 1986.
[22]
Wojciech K. Gwarek, “Analysis o f arbitrary shaped two- dimensional microwave
circuits byFinite-DifferenceTime-Domain method,”IEEE Trans. M T T , vol. 36,
pp. 738-744, April 1988.
[23]
X. Zhang, J. Fang, K. Mei, and Y. Liu, "Calculation of the dispersive characteristics
of microstrips by the time-domain finite difference method," IEEE Trans.
Microwave Theory and Tech., vol. 36, pp. 263-267, Feb. 1988.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
133
[24]
J. Fang, X. Zhang and K. Mei, "Dispersion characteristics of microstrip lines in the
vicinity of a coplanar ground," Electronic letters, vol. 23, pp. 1142-1143, Oct. 1987.
[25]
X. Zhang and K. Mei, "Time-domain finite difference approach to the calculation
o f the frequency-dependent characteristics of microstrip discontinuities," IEEE
Trans. Microwave Theory and Tech., vol. 36, pp. 1775-1787, 1988.
[26]
X. Zhang and K. Mei, "Time-domain finite difference approach for the calculation
o f microstrip open-circuit end effect," IEEE MTT-S International Microwave
Symposium Digest, New York, NY, USA, pp. 363-366, May 1988.
[27]
X. Zhang and K. Mei, "Time-domain finite difference approach for the modeling of
microstrip components," IEEE AP-SInternational Symposium, San Jose, CA, USA,
pp. 1120-1123,June 1989.
[28]
Y. Liu and K. Mei, "Time-domain finite difference analysis o f electrically thin and
thick substrate micro strip antennas," Radio Science Meeting, San Jose, CA, USA,
Program and Abstract, pp. 265, June 1989.
[29]
J. Fang, “Time domain finite difference computations for Maxwell’s equations,”
Ph.D. dissertation, EECS Dept., Univ. California, Berkeley, 1989.
[30]
Jeffrey L. Young, "A Higher order FDTD for EM propagation in a collisionless
cold p l a s m a / 1 IEEE Trans. Antennas Propagat., vol. 44, pp. 1283-1289, Sept.
1996.
[31]
Matthew, N. O. Sadiku, Numerical Techniques in Electromagnetics. Boca Raton,
FL: CRC Press, 1992.
[32]
C. A Balanis, Advanced Engineering Electromagnetics. New Y ork, NY: John Wiley
and Sons, Inc., 1989.
[33]
M. Piket-May, A. Taflove, and J. Baron, "FD-TD modeling of Digital signal
propagation in 3-D circuits with passive and active loads," IEEE Trans. Microwave
Theory and Tech., vol. 42, pp. 1514-1523, August 1994.
[34]
J. Zurcher and F. Gardiol, Broadband Patch Antennas. Norwood, MA: Artech
House, Inc., 1995.
[35]
Charless W. Manry Jr., Shira L. Broschat, and John B. Schneider, "Higher-order
FDTD methods for large problems," AC ES Journal, vol. 10, pp. 17-29, July 1995.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
134
[36]
A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain
Method. Norwood, MA: Artech House, Inc., 1995.
[37]
D. Sheen, S. Ali, and J. Kong, "Application of the three-dimensional finite-difference
time-domain method to the analysis of planar microstrip circuits," IEEE Trans.
Microwave Theory and Tech., vol. 38, pp. 849-856, 1990.
[38]
D. Pozar, Microwave Engineering. Reading, MA: Addison-Wesley Publishing
Company , Inc., pp. 220-230, 1990.
[3 9] Mohammed F. Hadi, "A Modified FDTD (2,4) scheme for modeling electrically large
structures with high phase accuracy," Ph.D. Dissertation, Electrical Engineering and
Computer Engineering Dept., Univ. of Colorado, Boulder, 1992.
[40]
M. Kragalott, M. Kluskens, and W. Pala, "Time-domain fields exterior to a
two-dimensional FDTD space," IEEE Trans. Antennas Propagat., vol. 45, pp.
1655-1663, Nov. 1997.
[41]
R. Luebbers, K. Kunz, and K. Chamberlin, "An interactive demonstration of
electromagnetic wave propagation using time-domain finite differences," IEEE
Trans, on Education., vol. 33, pp. 60-68, Feb. 1990.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
135
APPENDIX
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
136
Appendix I
The Finite Difference Approximations for Second and Third Order Derivatives
To obtain the finite difference approximations for second order derivatives, one can
go through the following steps:
First, the central difference approximation (CD) for first order derivatives for an Hy field
component is given by:
dH
r-y
-k
ij *
i - Hy v*-... i )
2
+ -
2
(L I)
where pz = l/(Az). The second order derivative is given by:
&P
dH
y
dz
dH
y
(L2)
dz
d2H
(1.4)
Similarly, the third order derivative is given by:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
137
_- jdl r ^ W
y tjJc
dz
dz3
i
a 5>
dz2
Substituting equation (1.4) into (1.5), one can obtain:
y 'J*
-3
-
> </*♦-
dz3
a .6)
> V * --
The expression for the third order derivative of Hy with respect to z and then its second
derivative with respect to time is given by inspection as:
1
n»—
iL
2
a3# .
H
2 3 - 3H
2 {
1
1
2
2
1
/»♦—
- 3 H
P? Pr
1
+3 ^
y >.y.
- 3 H
yy.
3
n- —
2
y ».y.
2
3
2
- 3H
y i . j . k
3
2
yy.
--
3H
y '-y-
2 x
1
n+ —
2
v*--
i
n-»-—
H y , j\, k--3 ]
j<
n-—
3
2
n- —
y
1
1
n-—
n-—
+3 H
2
y
2
1
y ijJ c * -
n+—
(H
1
y > J **-
d t2
- 2 ( H
1
n-L
2
- H
y >.y. fc- -
3
rt~ —
^ y '-y.
2 * 3 }
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(17)
138
Vita
Adel Abdin was bom in Monofiya, Egypt, in 1957. He received the B. S. degree in
Electrical Engineering from Military Technical College, Egypt, in 1980. He received the
Diploma and the M. S. in Electrical Engineering from Ain Shams University, Egypt, in 1987
and 1990, respectively.
In August 1994, he enrolled in the graduate school at the University of Mississippi.
Adel Abdin is a student member of the IEEE.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
s ?
IMAGE EVALUATION
TEST TARGET (Q A -3 )
A
o
1.0
HI
Im 1 2 .2
■3.6 11=35
*o II 2.0
l.l
1.8
1.25
1.4
1.6
150 m m
IIW IG E . Inc
1653 East Main Street
Rochester, NY 14609 USA
Phone: 716/482-0300
Fax: 716/288-5989
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Документ
Категория
Без категории
Просмотров
0
Размер файла
4 208 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа