close

Вход

Забыли?

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

?

Quantification of induced electromagnetic fields inside material samples placed inside an energized microwave cavity by finite-difference time-domain (FDTD) method

код для вставкиСкачать
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 of
computer printer.
The quality o f this reproduction is dependent upon the quality of the
copy submitted. Broken or indistinct print, colored or poor quality illustrations
and photographs, print bleedthrough, substandard margins, and improper
alignment can adversely affect 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.
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.
Bell & Howell Information and Learning
300 North Zeeb Road, Ann Arbor, Ml 48106-1346 USA
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.
QUANTIFICATION OF INDUCED ELECTROMAGNETIC FIELDS
INSIDE MATERIAL SAMPLES PLACED INSIDE AN ENERGIZED
MICROWAVE CAVITY BY FINITE-DIFFERENCE TIME-DOMAIN
(FDTD) METHOD
By
Yao-Chiang Kan
A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment o f the requirements
for the degree o f
DOCTOR OF PHILOSOPHY
Department of Electrical and Computer Engineering
2000
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
UMI Number: 3000556
___
®
UMI
UMI Microform 3000556
Copyright 2001 by Bell & Howell Information and Learning Company.
All rights reserved. This microform edition is protected against
unauthorized copying under Title 17, United States Code.
Bell & Howell information and Learning Company
300 North Zeeb Road
P.O. Box 1346
Ann Arbor, Ml 48106-1346
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ABSTRACT
QUANTIFICATION OF INDUCED ELECTROMAGNETIC FIELDS
INSIDE MATERIAL SAMPLES PLACED INSIDE AN ENERGIZED
MICROWAVE CAVITY BY FINITE-DIFFERENCE TIME-DOMAIN
(FDTD) METHOD
By
Yao-Chiang Kan
The investigation o f the heating o f a m aterial sample in an energized microwave cavity
requires the understanding of the interaction o f the electromagnetic fields with the
material sample in the cavity. The key factor fo r this understanding is to quantify the
distribution o f the induced electromagnetic field inside the material sample placed inside
the cavity.
The goal o f this research is to solve M a x w e ll’s equations in an electromagnetic cavity
in the presence o f a material sample based on the finite-difference time-domain (F D T D )
method, which has been successfully applied to several areas in electromagnetics. This
study is based on Yee algorithm, a second-order F D T D scheme, and an improved fourthorder F D T D scheme.
The numerical dispersion equations o f Yee and other fourth-order F D T D schemes are
first discussed and the disadvantages o f Yee scheme are discussed. A n im plicit staggered
fourth-order F D T D method is then employed to calculate the filed distribution in a
rectangular cavity with lossless or lossy samples. The quality factor and the resonant
frequency o f a cavity are obtained by a derived time-domain Poynting’s theorem and also
by Prony’s method. The quality factors calculated by these two different methods are very
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
consistent.
In applying the F D T D method to cylindrical coordinates, the singularities at the center
o f cylindrical coordinates become the major problem. The body o f revolution (BOR)
FD TD method is applied to solve the filed distributions in the cylindrical cavities loaded
with samples with symmetries. The treatment in BOR FD TD method is quite straight
forward. Moreover, the BOR F D T D method is a 2.5D FD TD method which is much more
computationally efficient than a 3D F D T D method. For a sample with any arbitrary shape,
the general 3D cylindrical FD TD method is needed to do the calculation. The traditional
3D cylindrical FD TD method encounters the difficulties o f the mode-dependent source
implementation and the treatment o f singularities. In this study, a general fourth-order
FD TD method is proposed to overcome the problems encountered in a traditional 3D
cylindrical FD TD method.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
To my parents and fam ily
iv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ACKNOWLEDGMENTS
First and foremost, I would like to express m y deep gratitude to m y advisor, Dr. KunM u Chen, for introducing me to Electromagnetics, and his invaluable help and patience.
Without his guidance and expert knowledge in Electromagnetics, this dissertation would
not have been possible. I am grateful for the opportunity to learn from his example as both
a researcher and an engineer.
I would like to thank the other members o f my committee, Dr. Dennis Nyquist, Dr.
Edward Roth w ell, and D r. Byron Drachman for many helpful comments on this research.
M y special thanks to Dr. Rothwell for giving constructive criticism in many fruitful
discussions. Special thanks go to Dr. Leo Kempel for providing me a workstation to run
my simulation program.
I would also like to thank Ms. Jackie Carlson, director o f D ivision o f Engineering
Computing Services, for appointing me as a graduate assistant from 1994 to 1998. M y
special thanks to Dr. G uilin Cui, owner o f Vertex Computer in Lansing, fo r providing me a
job for the past year. W ithout those financial supports, I would not be able to study my
Ph.D. at M ichigan State University.
M y very special thanks to my w ife, Ya-Ling Peng, who have been taking care of our
children and me w ell so I am able to focus on m y research and work.
Finally, I would like to thank my parents fo r providing me with the opportunity to
complete this study. Their love and support accompanied me through the years of this
work.
v
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TABLE OF CONTENTS
L IS T OF T A B L E S ....................................................................................................................... ix
L IS T OF F IG U R E S ...................................................................................................................... x
C H A PTER 1
IN T R O D U C T IO N ........................................................................................................... 1
C H A PTER 2
S O L V IN G M A X W E L L ’S E Q U A TIO N S B Y F D T D IN R E C T A N G U LA R
C O O R D IN A T E ............................................................................................................... 5
2.1
Frequency-Dependent FD TD Formulations with Second-Order Leapfrog
Time-Stepping of M axw ell’s Equations .........................................................7
2.1.1 The Scalar Equations o f M axw ell’s Equations ................................7
2.1.2 The Finite Difference Equations ........................................................ 9
2.2
The Yee’s A lg o rith m .........................................................................................17
2.2.1 Stability C ondition...............................................................................19
2.2.2 Numerical Dispersion ........................................................................ 19
2.3
The T y(2,4) (FD )2TD Algorithm .................................................................. 26
2.3.1 The Fourth-Order Space D erivatives............................................... 26
2.3.2 Dispersion Analysis for Explicit Staggered S chem e.................... 27
2.3.3
Dispersion Analysis for Im plicit Staggered S chem e.................... 29
2.3.4 Calculation o f Derivative o f E J n T y(2,4) .................................... 33
2.3.5
Calculation o f Derivative o f H in T y (2 ,4 ) .................................. 35
2.4
Excitation Source and Power Analysis ......................................................... 38
2.4.1
Excitation Source................................................................................ 38
2.4.2 Power Analysis .................................................................................. 39
2.5
Prony’s M e th o d .................................................................................................44
2.5.1
T h e o ry ................................................................................................. 44
2.5.2
Estimation o f Resonant Frequency and Quality Factor by Prony
M e th o d ............................................................................................................... 47
2.6
A Single Empty Cavity with PEC W a lls .......................................................48
2.6.1
C onfiguration...................................................................................... 48
2.6.2
Numerical Results o f Field Distributions ...................................... 48
2.7
A Lossless Loaded Cavity with PEC W a lls ...................................................65
2.7.1
Quasi-cubic C a s e ................................................................................65
vi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.8
2.7.2
Thin Square Plate C a s e ........................................................................68
2.7.3 Narrow Strip C a s e
................................................................. 72
A Lossy Dielectric Loaded C avity with PEC Walls ..................................... 76
2.8.1 C onfigurations.......................................................................................76
2.8.2 Numerical Results and Discussions................................................... 77
C H A PTER 3
S O L V IN G M A X W E L L ’S E Q U A TIO N S B Y B O D Y O F R E V O L U T IO N '
F D T D .................................................................................................................................88
89
3.1
The BO R Formulation o f M axw ell’s E quations[28]...........................
3.1.1 M ode Selection in B O R A lg o rith m ....................................................93
3.1.2 The B O R -FD TD Form ulation[28]........................................ ... . . . . 95
3.1.3 Singularity in B O R -FD TD Formulation atp = 0 ................ - . . . . 99
3.2
Surface Impedance Boundary C o n d itio n ............................................- . . . 100
3.2.1 Planar Surface Impedance Boundary C o n d itio n ................. _ . . . 100
3.2.2
F D T D Implementation o f Planar Surface Impedance . . . . . . . . . 103
3.2.3
Tim e-Dom ain Approximation by Prony’ M e th o d .............. . . . 104
3.2.4
Frequency Domain A pproxim ation....................................... . . . 105
3.2.5
Z Transform and D igital Filters Approach ........................... . . . 106
3.2.6
Fields Calculation on the Cavity W a ll.................................... . . . 108
3.3
An Em pty Cylindrical Cavity .............................................................. _ . . . 109
3.4
A Loaded Cylindrical C a v ity .............................................................. ... . . . 118
3.4.1
Small cylindrical sample for TM 012 m o d e ........................... . . . 118
3.4.2
Thin rod case for T M 0 1 2 m o d e ........................................... ... . . . 119
3.4.3
Thin disk case fo rT M 0 1 2 m o d e ........................................... . . . 119
3.4.4
Small cylindrical sample f o r T E lll mode ........................... . . . 121
C H A PTER 4
S O L V IN G M A X W E L L ’S E Q U A TIO N S B Y FD TD IN C Y L IN D R IC A L
C O O R D IN A T E S ........................................................................................................... 136
4.1
Three Dimensional F D T D Representation of M axw ell’s Equations in
Cylindrical Coordinates................................................................................... 137
4.2
The Second-order Cylindrical F D T D Scheme[35] ...................................... 143
4.2.1
The Second-order Cylindrical FD TD E quations...........................143
4.2.2
F D T D Calculations at p = 0 ..............................................................145
4.2.3
Source Implementation .................................................................... 146
4.2.4
Numerical Results and Discussions................................................ 148
4.3
T y(2,4) Cylindrical F D T D S ch em e................................................................ 165
4.3.1
Derivatives o f Fields in p direction ...........................................165
4.3.2
Derivatives o f Fields in <)>direction.............................................172
4.3.3
Derivatives o f Fields in z direction.............................................175
4.3.4
F D T D Calculations at p = 0 ...................
176
4.4
Tim e Domain Finite Difference Equations of Constitutive
Relations ........................................................................................................... 179
v ii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
C H A PTER 5
C O N C L U S IO N S .......................................................................................................... 182
A P P E N D IX A
D E R V LA TIO N O F FO U R TH -O R D E R F IN IT E D IFFER EN C E
A P P R O X IM A T IO N ................................................................................................... 183
A P P E N D IX B
D E R V IA T IO N O F N U M E R IC A L D IS P E R S IO N R E L A T O IN FO R IM P L IC IT
STA G G ERED S C H E M E .............................................................................................189
A P P E N D IX C
D E R V IA T IO N O F E Q U A T IO N (2.107)
.............................................................. 201
A P P E N D IX D
D E R V IA T IO N O F E Q U A TIO N S (3.7) A N D (3.8) ..............................................204
B IB L IO G R A P H Y ..................................................................................................................... 206
v iii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
LIST OF TABLES
Table 2.1
E ws/ E ns ratio at different tim e s ......................................................................67
Table 2.2
Permittivities mapping at ze = In s , CD= 2 ;t(2 .4 5 e 9 )............................... 77
Table 2.3
The loss power and stored energy..................................................................... 78
Table 2.4
Properities o f four different lossy cases.......................................................... 79
Table 3.1
BOR reprenstation o f M axw ell’s equations................................................... 92
Table 3.2
Rational Approximation Results..................................................................... 107
ix
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
LIST OF FIGURES
Figure2.1
FD TD space lattice............................................................................................. 11
Figure2.2
Variation of the numerical phase velocity with C F L values at R = 10
and <j) = J t /4 ..................................................................................................... 22
Figure2.3
Variation o f the numerical phase velocity with R at a = 0.9 and
<t> = it / 4 .
23
Variation o f the numerical phase velocity with <j> at a = 0.9 and
R = 10.
24
Figure2.4
Figure2.5
Variation of the numerical phase velocity with 1/R at a = 0.9 and
= ti/4 ................................................................................................................. 25
Figure2.6
Comparison of numerical phase velocity between explicit 4th-order scheme
and Yee scheme at R = 5, a = 0.24, and <J) = n/4............................................31
Figure2.7
Comparison of numerical phase velocity between explicit 4th-order scheme
and Yee scheme at a = 0.24, and (j) = tt/4 ........................................................32
Figure2.8
The lattices of Ey and dEy/ d z along the z axis at fix x and y................... 33
Figure2.9
The lattices of Hz and d H z/ d y along the y axis at fix x and z...................35
Figure2.10
The configuration o f the empty rectangular cavity with PEC boundary. The
excitation probe is located on one of the six faces........................................50
Figure2.11
The x dependence o fE x at z= 0.025974 and t= 0.26658 m s fo rT E O ll mode.
..............................................................................................................................51
Figure2.12
The x dependence o f E x aty=0.034632 andt=0.26658 ms fo rT E O ll mode.
..............................................................................................................................52
Figure2.13
The y dependence o f Ex at x= 0.025974 and t= 0.26658 ms fo rT E O ll mode.
..............................................................................................................................53
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure2.14
Figure2.15
The z dependence o f E x at x= 0.034632 and t= 0.26658 ms for T E 0 1 1 mode
................................................................................................................................54
The Ex variation along the yz plane at x= 0.034632 and t= 0.26658 ms for
TE011 mode......................................................................................................... 55
Figure2.16
The E field on the xz plane at y= 0.034632 and t= 0.26658 ms for T E 0 1 1
mode.
56
Figure2.17
T h exd ep en d en ceo fE xatz= 0.03192an d t= 0.31019 m s fo r T M lll mode.
................................................................................................................................57
Figure2.18
The x dependence o f E y at z= 0.04256 and t= 0.31019 ms fo r T M 1 1 1 mode.
................................................................................................................................58
Figure2.19
The x dependence o f E z at y= 0.04256 and t= 0.31019 ms for T M 1 11 mode.
................................................................................................................................59
Figure2.20
The E field on the xy plane at z= 0.04265 and t= 0.31019 ms for T M 1 11
mode..................................................................................................................... 60
Figure2.21
The E field on the xz plane at y= 0.04265 and t= 0.31019 ms for T M 1 11
mode...................................................................................................................... 61
Figure2.22
The E field on the yz plane at x= 0.04265 and t= 0.31019 ms for T M 1 11
mode..................................................................................................................... 62
Figure2.23
The time variation o f Ex at x=0.04329, y=0.04329, and z=0.034632. . . 63
Figure2.24
The frequency response o f Figure 2.23........................................................... 64
Figure2.25
Dimensions o f the rectangular cavity and the loaded m aterial sample. The
center of the material sample is consistent w ith the center o f the cavity.
66
Figure2.26
The variation of electric fields in the y direction at t = 0.2535p s
Figure2.27
The variation of E y in the x direction at r = 0.24289p s . The line with star
symbol is the calculated values and the line w ith circle symbol is the fitted
values for the empty cavity............................................................................... 70
Figure2.28
Variations o f E ws/ E ns in the x-directions. Each curve represents this ratio
as a function o f x for different locations o f z. The relative perm ittivity o f the
thin square plate m aterial sample is er = 2 .5 . The solid line with symbols
is the ratios at t = 0.24289ps and the dash line w ith symbols is those at
t = 0.24301 p s .................................................................................................. 71
xi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
69
Figure2.29
The variation o f electric fields in the y directions at x=0.34839m ,
z=0.057093m , and t=0.252m s...........................................................................74
Figure2.30
The ratios o f E ws/ E ns in y directions............................................................. 75
Figure2.31
The flow chart o f calculating the field distribution at steady state for high Q
cavities...................................................................................................................80
Figure2.32
Stored energy for the cavity with material sample o f e'r(o j) = 2.5 and
e"r((o) = 0 .1 . Note that one time step is equal to 4.25442 ps and this is a
downsampling plot.............................................................................................. 81
Figure2.33
Instaneous dissipated power for the cavity w ith m aterial sample of
e'r(co) = 2.5 and e"r(co) = 0 .1 . Note that one tim e step is equal to
4.25442 ps and this is a downsampling plot................................................... 82
Figure2.34
The plot o f the fitting and calculated data for stored energy for the first case.
The line with circle is the fitting data and that w ith cross is the calculated
power data. The average input power is 2.45x10”10 and the time average
stored energy is 8.905e”18.................................................................................. 83
Figure2.35
The plot of the fitting and calculated data for dissipated power for the first
case. The line with circle is the fitting data and that w ith cross is the
calculated power data. The average input power is 2.45xlO ”10.................. 84
Figure2.36
The plot o f the difference between fitting and dissipated power data. The
one time step is equal to 4.25442 ps.................................................................85
Figure2.37
Field distributions o f E y along x axis at 0.14868 ms. £'r(co) is equal to 2.5
and e"r(co) is 0.1.................................................................................................86
Figure2.38
The ratio of calculated and fitting Ey along x axis at 0.14868 ms
Figure3.1
The field locations for B O R F D T D in time and space................................. 98
Figure3.2
F D T D configuration fo r cavities with FEC w a ll........................................102
Figure3.3
Coordinates for the incident and reflected plane waves upon a lossy
surface .............................................................................................................. 103
Figure3.4
The physical configurtion o f a cylindrical cavity loaded with a material
sample................................................................................................................. I l l
Figure3.5
The variation o f Ep o f T M 0 1 2 along the p and z directions in a PEC empty
cylindrical cavity............................................................................................... 112
Xll
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
87
Figure3.6
The variation o f Ez o f T M 012 mode along the p and z directions in a PEC
empty cylindrical cavity.................................................................................. 113
Figure3.7
Plot o f E field o f T M 012 mode on the p-z plane in an empty PEC cavity.
............................................................................................................................. 114
Figure3.8
The variation o f Ep o f T E 1 11 mode along the p and z directions in a PEC
empty cylindrical cavity.................................................................................. 115
Figure3.9
The variation o f E<j, o f T E 1 11 mode along the p and z directions in a PEC
empty cylindrical cavity.................................................................................. 116
Figure3.10
Plot of E field o f T E 1 11 mode on the p-z plane in an empty PEC cavity.
............................................................................................................................. 117
F g m re 3 .ll
Plot of Ep o f TE111 mode versus time in an empty FE C cavity with
a = 102 S /m and <r= 102 S /m with m = 1 ........................................ 120
Figure3.12
Plot of Ez o f T E 1 11 along r direction in an empty FE C cavity with
conductivity 104 (S/m ).....................................................................................122
Figure3.13
Plot of E field o f T E 1 11 mode on the p-z plane in an empty FEC cavity with
conductivity 104 (S /m ).....................................................................................123
Figure3.14
Plot of Ep o f TM 012 mode in a PEC cavity loaded w ith a small cylindrical
material sample................................................................................................. 124
Figure3.15
Plot of Ez of TM O 12 mode in a PEC cavity loaded w ith a small cylindrical
material sample................................................................................................. 125
Figure3.16
Plot of the ratio of o f T M 012 mode along the z direction in the PEC cavity
loaded with a small cylindrical material sample.......................................... 126
Figure3.17
Plot of Ep o f TMO 12 mode in a PEC cavity loaded w ith a thin rod material
sample.................................................................................................................127
Figure3.18
Plot of Ez of TMO 12 mode in a PEC cavity loaded w ith a thin rod material
sample.................................................................................................................128
Figure3.19
The ratio o f E z/ E zl o f TMO 12 mode in a PEC cavity loaded with a thin rod
material sample................................................................................................. 129
Figure3.20
Plot of Ep o f TMO 12 mode in a PEC cavity loaded w ith a thin disk material
sample.................................................................................................................130
x iii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure3.21
Plot o f E z o f TMO 12 mode in a PEC cavity loaded with a thin disk material
sample....................................
131
Figure3.22
The ratio o f E z/ E l_ o f TM 012 mode in a PEC cavity loaded with a thin disk
material sample................................................................................................. 132
Figure3.23
Plot o f Ep o f T E 1 11 mode in a PEC cavity loaded with a small cylindrical
material sample................................................................................................. 133
Figure3.24
Plot o f E^ o f T E 1 11 mode in a PEC cavity loaded with a small cylindrical
material sample................................................................................................. 134
Figure3.25
Plot o f Ez o f T E 1 11 mode in a PEC cavity loaded with a small cylindrical
material sample................................................................................................. 135
Figure4.1
FD T D lattice for cylindrical coordinates...................................................... 141
Figure4.2
The diagram o f order o f FD TD calculations along tim e axis. The meanings
o f those steps are listed below.
(1) Using (4.9) to (4.11)
(2) Desired time stepping scheme for D .
X
A
(3) Constitutive relation o f D and E which depends on material models.
(4) Using (4.12) to (4.14). A t this point, the boundary conditions are
applied.
—
-X
(5) Desired time stepping scheme for B .
-X
_x
(6) Constitutive relation o f B and H which depends on material models.
............................................................................................................................. 142
Figure4.3
TE011, TE 111, TM 011, and TM 012 modes excitation techniques in a
cylindrical cavity[35]....................................................................................... 147
Figure4.4
The plots o f total stored energy of TMO 12 mode in a PEC empty cylindrical
cavity. The upper figure is calculated by using the traditional source
implementation and the lower one by using the B H source implementation.
............................................................................................................................. 149
Figure4.5
The variation o f E p o f TMO 12 mode along the p direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation................................................................................................. 150
Figure4.6
The variation o f E p o f TMO 12 mode along the z direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
xiv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
source implementation and the lower one by using the BH source
implementation.............................................................................................. 151
Figure4.7
The variation o f Ez o f T M 0 1 2 mode along the p direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................152
Figure4.8
The variation o f Ez o f TMO 12 mode along the z direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................153
Figure4.9
The plots o f total stored energy o f T E 1 11 mode in a PEC empty cylindrical
cavity. The upper figure is calculated by using the traditional source
implementation and the lower one by using the B H source
implementation...................................................................................................154
Figure4.10
The variation o f E p o f T E 1 11 mode along the p direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................155
Figure4.11
The variation o f E p o f T E 1 11 mode along the <j) direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................156
Figure4.12
The variation o f E p o f T E 1 11 mode along the z direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................157
Figure4.13
The variation o f E^ o f TE111 mode along the p direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................158
Figure4.14
The variation o f Ep o f T E 1 1 1 mode along the <j) direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation...................................................................................................159
Figure4.15
The variation o f E p o f T E 1 11 mode along the z direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
XV
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
implementation.............................................................................................. 160
Figure4.16
The variation o f E p o f T M 012 along the p and z directions in a PEC empty
cylindrical cavity. The above figure is calculated by using the traditional
source implementation and the below one by using the B H source
implementation..................................................................................................162
Figure4.17
The variation o f Ez o f TM O 12 along the p and z directions in a PEC empty
cylindrical cavity. The above figure is calculated by using the traditional
source implementation and the below one by using the B H source
implementation..................................................................................................... 163
Figure4.18
Plot of the ratio o f o f TM O 12 mode along the z direction in the PEC cavity
loaded with a small cylindrical material sample.......................................... 164
Figure4.19
The locations o f H fie ld and its corresponding derivatives in
p direction.......................................................................................................... 166
Figure4.20
The locations o f E field and its corresponding derivatives in
p direction.......................................................................................................... 171
Figure4.21
The FD TD lattice along the <j) direction. 174
F ig u re d
Approximation of intergrant in ( C . l ) ........................................................... 202
x vi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
C H A PT E R 1
IN TRO DUCTIO N
The research reported in this dissertation was motivated by the investigation of
microwave heating of material samples. Microwave heating techniques have been widely
utilized in industrial process. Since the microwave heating o f material samples is usually
conducted within an energized electromagnetic cavity, to understand the heating
mechanism it is essential to study the interaction of the microwave field with a material
sample in an electromagnetic cavity. The key factor in understanding this interaction is to
quantify the induced electromagnetic field inside the material sample by the cavity field.
The finite-difference time-domain (F D T D ) method is employed in this dissertation to
quantify the electromagnetic field inside an E M cavity loaded w ith material samples.
The finite-difference time-domain (FD TD ) method has been used in computational
fluid dynamics (CFD) [1] for a long time and yields very accurate results for C FD
problems. In 1966, Kane Yee originated a set o f finite-difference equations for the tim edependent M axw ell’s curl equations system for the lossless materials case [2]. The F D T D
method was not popular in (Electromagnetic) E M research area until late 1980 and
becomes a very popular method in E M area between 1993 and 1997. Regarding to the
F D T D method applied to eigenvalue problems in E M research, Choi and Hoefer [4]
published the first FD TD simulation o f waveguide/cavity structures in 1986. There are
1
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
several papers, [5 ]-[10 ], that utilized the FD TD method to investigate the fields and power
distributions in a loaded E M cavity. The paper by Torres and Jecko [9] provides a very
detailed analysis o f microwave heating by combining the M a x w e ll’s equations and the
heat transfer equations. The F D T D method used in this paper is called (F D )2T D method
since the constitutive parameters are assumed to be frequency-dependent. This (F D )2T D
method is a standard method to investigate E M interaction w ith a lossy material sample
using F D T D scheme. Moreover, the combined electromagnetic and thermal F D T D
algorithm used in this paper provides a basic fram ework for the FD TD calculation.
However, the dissipated power model in this paper is not clear and the improved model
w ill be discussed in this dissertation.
M ost papers in the F D T D literature are based on second order spatial and temporal
approximations which are originated from Yee algorithm. Due to the requirement o f
accuracy and performance, several higher order F D T D methods have been proposed [11][15]. Among these higher order F D T D methods, Ty(2,4) F D T D method, which uses the
im plicit staggered fourth order finite difference approximation in space and the explicit
second order finite difference approximation in time, provides the most promising features
[17]. Combining with the (F D )2T D method, the Ty(2,4) (F D )2T D method becomes the
essential F D T D technique in the investigation for the eigenvalue problems in E M research
area.
In this dissertation, the T y(2,4) (F D )2T D method in rectangular and cylindrical
coordinates are studied. For a cylindrical cavity with azim uthal symmetry, the body o f
revolution (B O R ) scheme is employed to facilitate the F D T D calculation. In this case, the
BOR F D T D method can give very accuracy results with excellent performance; hence, the
2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Ty(2,4) F D T D method is not used.
In chapter 2, a set o f general finite difference equations for M axw ell equations is
introduced. Also, the loaded material is modeled by Debye equations and its FD TD
representation is presented. Then a brief introduction o f Yee algorithm, stability condition
and numerical dispersion is presented. The fourth order spatial derivatives are presented
and the dispersion analysis is studied. A fter that, the T y(2,4) (F D )2T D method is presented
in details. A time-domain power analysis based on Poynting’s theorem is derived. In this
F D T D calculation fo r cavities, the Prony’s method is employed to estimate the quality
factors (Q ) of cavities. The numerical results of a single empty cavity w ith perfect
electrical conductive(PEC) walls, a cavity loaded cavity w ith a lossless material sample
and PEC walls, and a lossy dielectric loaded cavity w ith PEC walls are presented at the
end o f this chapter. The numerical results are shown to be consistent with the theoretical
estimation.
In chapter 3, the body o f revolution (B O R ) F D T D form ulation o f M axw ell’s equations
is derived. Mode selection and source implementation in B O R FD TD algorithm are
presented. The treatment fo r the singularity in BOR F D T D formulation is also presented.
In this chapter, the cavities with finite electrical conductive (FE C ) walls is studied and the
FEC w all is replaced by a planar surface impedance boundary condition (P S IB C ). The
F D T D formulation o f PSIB C is achieved by three different approachs and the frequency
domain approximation is chosen due to its simplicity. Num erical results o f an empty
cylindrical cavity w ith PEC and FEC walls and a loaded cylindrical cavity with PEC walls
are presented at the end o f this chapter. Consistent results are obtained by using this BOR
FD T D method fo r cavities with azimuthal symmetry.
3
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
In chapter 4, a general 3D F D T D method in cylindrical coordinates is considered and
studied. The disadvantages o f conventional second order FD TD method in cylindrical
coordinates are presented and its improvements are also proposed. The treatment for the
singularity in conventional cylindrical F D TD is mode-dependent and difficult to be
generalized. B y the introduction o f Ty(2,4) F D T D in cylindrical coordinates, a general
treatment for the singularity in cylindrical F D T D can be obtained. W ith the Debye or
Lorent models fo r loaded material samples, a general Ty(2,4) (F D )2T D method in
cylindrical coordinates is obtained.
Some derivations that are useful in this dissertation are provided in Appendices.
4
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
C H A PT E R 2
SO LVING M A X W ELL’S EQUATIONS BY FDTD
IN RECTANG ULAR COORDINATE
In this chapter we considers the frequency-dependent im plicit, fourth-order F D T D
form ulation in the rectangular coordinate system and its application to rectangular
cavities. The finite difference approximation to time-stepping in those formulations is o f
second-order accuracy. The second-order accuracy in time-stepping combined with
im plicit staggered fourth-order accuracy in space, denoted by T y(2,4)[17], is the focus o f
this chapter. Another scheme, Ty(4,4), that uses fourth-order in time-stepping with
im plicit fourth-order in space-stepping is not discussed because o f the following reasons:
(1)
In multistage time discretization schemes (e.g., Runge-Kutta schemes with three
or more stages), boundary conditions must be applied at intermediate levels, then
memory requirement and computer running time are increased.
(2)
The accuracies o f Ty(2,4) and Ty(4,4) are made comparable by choosing an
appropriate time step[14][17] although the stability of the F D T D formulation o f
Ty(4,4) may be improved[15].
(3)
The Ty(2,4) is nondissipative, while Ty(4,4) introduces a slight dissipation.
In F D T D calculations for cavities, a many time steps are usually required if the quality
5
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
factor, Q, o f the cavity is high. The Ty(2,4) is used to speed up the computation time since
coarser meshes are chosen than that in the original Yee scheme. A dissipation scheme like
Ty(4,4) is not a good choice fo r a long time integration.
The frequency-dependent F D T D
formulation with second-order leapfrog tim e-
stepping o f M axw ell’s equations is discussed in section 2.1. The one-relaxation Debye
equation is used to account fo r the frequency-dependent properties.
The original Yee’s algorithm is derived from this general scheme for validation in
section 2.2. The stability condition o f Yee’s algorithm is also presented and the dispersion
in Yee’s algorithm is then explained.
In section 2.3, the Ty(2,4) (F D )2T D scheme is obtained from this general scheme. The
higher-oder spatial schemes requires much fewer mesh points than Yee’s scheme does for
the same dispersion error. Although the former usually complicates the F D T D scheme and
require a smaller tim e step for stability, the computer running tim e is complemented by the
mesh reduction. The numerical dispersions of explicit fourth-order and Ty(2,4) (F D )2T D
schemes are discussed in section 2.3.2 and section 2.3.3. The compact finite difference
schemes requiring special treatments at the boundary and those numerical treatments are
discussed in section 2.3.4 and section 2.3.5. The stability o f introducting the numerical
boundary treatment is discussed. The treatment of the finitely conducting (FEC) boundary
is also discussed.
The excitation source is discussed in section 2.4. For the empty cavities with PEC
boundary and loaded cavities w ith lossless dielectric material samples, the Q values o f the
cavities are infinite so only transient-state solution are obtained to validate the program.
The source used for this case is a Blackm an-Hanis type that has very low sidelobes. For
loaded cavities w ith lossy dielectric material sample, a single frequency sinusoidal source
6
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
is used. For cavities, the Q factor and the resonant frequency are most desired quantities
that need to be calculated. The Prony method in section 2.5 provides a method to estimates
those two values without running a lenghty F D T D computation. The numerical results and
discussion are presented in section 2.6, section 2.7 and section 2.8.
2.1 Frequency-Dependent FDTD Formulations with Second-Order
Leapfrog Time-Stepping of M axwell’s Equations
2.1.1 The Scalar Equations of M axwell’s Equations
In differential form, M axw ell’s equations in a dielectric dispersive medium can be
written as
VxE = - [ H l ^
( 2 - 1)
V x H = ~ + [g ]E + J S
at
where [cr] and [p .], are electric conductivity, magnetic permeability which are nondispersive in tensor form , respectively. Js is the known impressed current source. The
above constitutive parameters are further assumed to have the biaxial tensor form in the
rectangular coordinate system given by
0
[a] =
0
0 CLy 0
0
0 az
where a represents the magnetic permeability or the electric conductivity.
For frequency-dependent dielectric material, the one-relaxation Debye equation is
7
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.2)
er(<o) = e’r(to )-7 e " r(co) = 8’^ +
E rs —E roo
&
(2.3)
1 + /C O T
where the subscript, r , denotes the word “relative”. Moreover, e'r and e"r are the real and
imaginary parts o f the relative permittivity. The xe is a new relaxation time constant
related to the original relaxation time constant j
e'
x = x
e
rJ
by [19]
+2
e’roo + 2
(2-4)
where e'rs and e'roo are the real part o f the complex perm ittivity at zero frequency and at a
very large frequency, respectively. Note that er(a)) is equal to e'r(o)) and e'roo can be any
values when xe is equal to zero, the frequency-independent case. The constitutive relation
between D and E is
e ,r 0
D (to ) = e0 0
0
where erx,
0
0 E ( co)
(2-5)
0 erz
, and erz are the relative permittivities in different directions and satisfy
Debye equation, (2.3).
The scalar equations o f (2.1) in time domain can be written as
Jd L
x dt
dz
dy
dEz
dEx
dx
dz
dH
3E x
3E y
dt
dy
dx
dHy _
y~Tt
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2 .6)
The equation (2.5) can be written as
(1 +/C 0Tex) D r = zQz'rxsE x + s 0z rxoozexj a E x
(1 + j m ey) D y = ^ ' r y s E y + £0 r y -^ e y J ^ y .
(2-8)
(1 +j(£>xez) D z = t Q€ rzsE z + z0€ rzoaxezj<QEz
then the inverse Fourier transform is applied; hence, the scalar equations are obtained as
follow
dDx
E x
^ex
— e0e rxsE
x
dE x
0^
rx ^ e x
az>v
a^v
D >-+ x ^ ~ a T = eo8'o '^ y + e0e’o-=oV a T •
(2-9)
3d,
Be
D z + Xez~df = ZoZ'rzsE z + eO€ r z - Tez-df
2.1.2 The Finite Difference Equations
For the time-stepping scheme, the leapfrog second-order scheme is applied to (2.6),
(2.7), and (2.9). D , E , and the temporal derivative o f H are evaluated at integer time step
but H and the temporal derivative o f D and E are evaluated at half integer time step. For
,/z
_ x .n + l/2
the tim e step n, the E | , D | , and H |
are considered to be the field distributions.
For the spatial discretion, the second-order Yee scheme or higher-order scheme can be
used. The finite difference approximation to the first derivative o f those quantities is
9
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
denoted by SaAg where a and (3 represent one of x, y, z, and t parameters and A denotes
one o f H , E, D fields. The spatial locations o f E and H are plotted in Figure 2.1 and D ,
locates at the same location as E does. Note that E x , dE x/ d t , and d D x/ d t are evaluated
at half integer spatial step along the x axis, but at integer spatial step along the other two
axes; so are the E y , E z , and their corresponding tim e derivatives. The H x, H y, H , are
evaluated at integer spatial step along the x, y, and z axes, respectively, but at h alf integer
space step along the other two axes; so is the time derivatives o f H . The spatialderivatives
of E
areevaluated at the same locations as H and that o f H are evaluated at the same
locations as E . Hence, (2.6), (2 .7 ), and (2.9) can be rewritten as
5 tH x\n
=
5 tH y\n
=
x\i,j+l/2,k+l/2
1
y \ i + 1/ 2, j , k + 1 / 2
1
Si r \
If + 1 /2 , j + 1 /2 , k
f si tt
1 /2
° t D x\
x \i+
Si 7-\ \n
OfT>
1 /2 ,/,*
1 /2
y U , j + 1/ 2, k
c
j--v
,n +
5 tD \
1 /2
Z\ i , j , k + W 2
y^. ,«+■ 1/2
D x\
x lf+ 1 /2 ,/,/t
0_H J
V z
— {& xE \ n
-S E J "
)
\
^
j,k
1 /2
iW + 1 /2
y \i,j,k+\/2
«
i/z +■1/2
+ t ,_ 8 .D “
1
x\i+l/2,j,k
.rt+ l/2
rr .n
o
“ 5 ^**v
^
o
y t
“ 8 v# r
y
ez ‘
n
, rt + 1 / 2
Ij, j, k + 1/2
^
j
1/ 2, k j
| / i+ 1 / 2
\ti +
G V VI
y
y\i,j+l/2,k
,
,n+l/2
./z + 1 /2
If'-t-1 / 2 , / , k
z
0
_
,
U
./i-t-1 /2
zlf,y,*+l/2
, _
t „ ,8 .£ J
.71+ 1/2
ex 1 X \i + 1/ 2 , j , k
r.
,
, / i + l / 2
z\ i , j , k + 1/2
,n+l/2
„
rr.^ fiy \
0 rz
(2-11)
y\i,j+l/2,k
T
,
^ y \ . , + 1/2_t + ^
, /t +■ 1 / 2
- J v\
" J z\
z\ i , j , k + l / 2
r
+ e0e'
, / z +■ 1 / 2
- Jxx\\ i + l / 2 , j , k
r
1 /2
\ - (yzE z\
x\i,j,k + l/2 j
rxs
r
Xl * + i / 2 , / , *
^
„ ,/i + 1/2
£J
0
1
1/ 2 , j + 1 /2 , k j
1 /2
^
(2-10)
x \i + 1/ 2, J , k + 1 / 2 )
^
,
= e0s'
, | , i+l/xt = ^
j
1/2
+ 1 /2
■*” 1 / 2
“ M U
,n + l/2
j,-. . n + 1 / 2
‘■ \ i , j + l / 2 , k + l / 2 j
z
y \ i + l / 2 , j t k)
rr
o
l f , y + 1 / 2 , fe
t t
H i + l / 2 , . / , j f c + 1 /2
y
= —
fS y ^ x f
- 5* £ ^y1f;+
\ l z \. y x \ i + l / 2 , j + l / 2 , k
8 J IJ
Dy
' Vl > + V 2 , k + W
If, y, k+
)
1 /2
f c* T7 \ti +
fS i
=
-V -l"
y\i,j+ l/2 ,k + l/2
-r
= Vl ° yv HZz\U + l / 2 ,
=
—fSFl”
JJ.XV
—
s; r
1 /2 tt(2-12)
ifi-t-l/2
cz 1 z\i, j, k+ 1/2
10
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
z
Figure 2.1 F D T D space lattice
11
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Apply the second-order leapfrog scheme to tim e stepping and use the average value to
approximate the D, E, J a.t n + 1 /2 time steps, equations (2. 10), (2.11), and (2.12) can be
rewritten as follow.
.n+l/2
_
.n — 1 /2
x \ i , j + 1 /2 , k + 1 /2
X U , j + 1 /2 , k + 1 /2
(2-13)
+ —(8,E\n
- 5 JEM "
p. V Z
j + 1/2, k + 1/2
y z\ i , j + l / 2 , k + l / 2 .
.n+l/2
. n—l / 2
y \ i + l / 2 , j , k + 1 /2
y l / + 1 /2 , j , k + 1 /2
(2.14)
+ —
pA
.n+l/2
_
)
- 5 _ £ , ln
^
z l i + 1 / 2 , j , k + 1 /2
c
*1 / + 1 /2 , j , k + 1/ 2 J
ift—l / 2
z li + 1 / 2 , / + 1 /2 , k
z lf + 1 / 2 , / + 1 /2 , k
(2.15)
+ — { byE x\n
pA
i« + l
Dy\
■*
,
e
X
x\i+l/2,j,k
2
® x ^
E_
li + 1 / 2 , j , k
n
2
x\i+l/2,j,k
_A tf
l
+A/
n+l
CTvA * A
2
A tf j
x\i+l/2,j,k
.n+l/2
z\i+l/2,j,k
,n
s
'N
y \i+ 1/2,j , k )
(2.16)
"j
= DAn
y \i,j+l/2,k
./c .
k
c» Tr . n + l / 2
- 5s ^ v
x \ i + 1/ 2 , j , k )
- * vy \ .i , j + 1 / 2 , k + M
.n + 1
)
x^ y \ i + l / 2 , j + l / 2 , k J
tn
f _
y\i,j+l/2,k
" T l >’l/,y+ l/2,
- 8 ^ /
V y
+ ^ ^ E A n+ l
2
DA
=
.n
2 V x l/+ l/2 ,/,ifc
yU,j+l/2,k
x\i+l/2,j+l/2,k
.n+1
^ ~ E x\
D vr
y
„
.n + l/2
5zH x 'i,. j -+
V Z
1 /2 , k
-
TT . n + l / 2
**
z.
_____
^
z li, j + 1 /2 , k.
+ /y l
y U , j + 1 / 2 , k.
12
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-17)
D \n + l
+ (1^ L e in + 1
z\ i , j , k + \ / l
c r . A r j-r |/Z
4
1 +
z li, y , / f c + l/ 2
' k
2 v
A f
r 1
k j, k + 1 /2
z n , j , k + 1 /2
, . / S r r , rt + 1/ 2
+ A/ 5J?V
£_
~Y
= D ,n
z\ i , j , k + 1/2
2
v.
+ /j"
- It,
) DxL+ L /2 , y , k
y lt , y , f c + l / 2
£ TT ,rt + 1/2
- 5 V# J
y
"\
x l i , y , J f c + l/ 2>/
(2.18)
fc + 1/2.
8 o (e' ^
+ e ,r.roo A f ) E x | . + 1 /2 > y > /t
(2.19)
1 /2 , j, k + 8 ° ( 8 r-rs _ e rx“ " A T ) E x \ i + 1 /2 , / , .
“ "A r)D xU
( X+
A / ) ^ ! : , j + 1 /2 , k
e ° ( e’^
+ 8 ^ ° ° a / ) ^ [ £- y + 1 /2 , *
(2.20)
= “ (■1 “ 4
)
|f- y + 1/ 2, k + 8° ( 8 ^ ~ 8 ^
A t ) D z l/, y, * + 1 /2
“ a T ) E y\i, y + 1/ 2, k
8 ° ( 8 rz* + 8'rz“ A t ) E z lf, j , k + 1 /2
(2 .21)
_ “ A r ) D z L-,y,/fc+ 1 /2 + 8° ( 8 rz* “ 8 ' z~ A ? ) E z \i, y , k + 1 /2
I f we solve the equations (2.16) to (2.18) and (2.19) to (2.21) simultaneously, we w ill
^
-X
obtain the finite difference expression for D and E as follow
D lB +l
x li + 1 / 2 , y, k
P
= Pl t D l*
C O
3xl ' y
w + 1 / 2 , y,k
tt
\n + 1 / 2
<> TT , f l + 1 / 2
z \i + 1/ 2, y, k ~
p3l ( j |" +1
-h
z
+ y I"
2 V x l / + l / 2 , y , ifc
It + 1
* 1 , + 1/ 2, y, k)
"j
x li+ \ / 2 , j , k )
13
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-22)
j + 1/2, *
K3> \( ^w X \ n+l/2
n+l/2
i,j+W 2,k - w
^
U,
j + 1/ 2 , k)j
J H l ( j \ n+ l
+ J \n
2 {. y \ i , j + 1/ 2, k
^ i I I L i /2
4 ,
)
j + 1/ 2, k)
=
1 /2
+p 3M
W 1
|/i
lx^ x\i+ l/2,j,k
jpN
£ x,s
i+ l/ 2 , j,k
/i + 1 /2
+ Y 3 ^ 5^
2
I
"t" 1
- s y- h
Ii, j , * + 1/2.
(2.24)
|—
* t/I
2 x\i + 1/2, /, k
y \ i + 1 /2 ,;,* ,
(2.25)
l"
4
+
1 /2 ,;,* ,
|W
2
,n+l/2
y—
r i/Z
ly 4 , ; + 1/ 2,* ^2y 4 , ; + 1/ 2 , *
+T 3 M
^
2
+ /
4 + 1 /2 ,;,* '
4 , ; + 1/ 2, *
4
t t
- o H \
cj TT ,n + 1 /2
l( + 1 / 2 , ; , *
Y
3 x /^ yT ,.71
+ 1
^3x1
n+
1
r*
s:
c , rk + l / 2
f / r +I
+y r
2 v z\i,j,k+ 1/2
zI/,/,/:+ 1/2.
—
“
(2.23)
z
f + I/2
k ; + i/2 ,*
i” * 1
r\
ln + 1/2
1
z li,; + i/2 ,* y
+ y i"
V y \ i , J + 1/ 2, *
|/l+ l
-5 ^
(2.26)
1
y l f , ; + 1 /2 , * j
1^
n
T
-T ,/i
|W
I
i , ; , * + 1 /2 ~~ Y l r ^ z l / , ; , * + 1 /2 ^ 2z
4 . ; .,*
, + 1/2
+ Y3z f 8 A f +1/2
- 8 i / |/I+ 1/2
")
zv
y\i ,j , k + 1/2
y
! ( , ; , * + 1/ 2;
/ t lIn + 1
+. /r I 71
2 V zli, j, * + 1 /2
211, ; , * + 1/ 2.
^ 3 z f
w h e re
14
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.27)
^ea
CTa (Tm - y ) + e o ( e'ra, + e'rra°° At ,
P l a = — :------------:---------------------------- ~ ----< M
* « x
+
f ) +E«(e''
T
J +
8 0[
ra s
+
e V
a
»
^
(2-28)
a
P2a = — ------------- 2 . ^
(2.29)
^ ( y + ^ o j + ^ e ’ra s + e' raw ^^ - a
2T,
a
v
P3a = ----- :------------
’ A * '- ----- — -----
CTa ( Y + T ea ) + e0( e’ra s
t
Via = o—j ( A—1—
•
r r
(2.30)
4 - 8° ’ ra°» ^ ^
—
'a ^ 2 + T«ccJ + e o(^e ray + e r a ~
< 2 -3 , )
^
^Tett
’ Ar
—
-----------------: ------------- 7---------------------------- 57—
Tea'
CTa ( f + Tm ) + e o(e,rcc + e'ra=o A r
^ a ( y + ^ a ) + e o(e'm s-e'
Y
20C
=
y 3a = -------
<2 -3 2 >
Ar + 2x a
_ £ 2 -------------------------------------------------- (2.33)
£ rcc°°
ra s + 0
CTa( y +Tea) +8°(e'
“ ce a ’
fa
and a = x , y , o r z .
Equations (2.13) to (2.15) and (2.22) to (2.27) are the finite difference equations that
are used to calculate the field strength. When r ea is equal to zero, the frequencyindependent case, equations (2.28) to (2.33) become
15
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
At_
eOe
P la
=
CTcc 2
2
ra s
— •
(2.34)
A r’
eOe r ^ + CTa y
P2<x = ° -
P3a
e0e'
(2.35)
Ar
= -------------------------------------------------------------------------------- (2-36)
eoe'ras + (ya y
Y la
=
-------------------------------------------------------------------------------- (2-37)
eOe'ra , + ^ a y
Y 2a =
1,
(2.38)
and
Ar
Ysce = --------------------------------------------------------------------------------(2-39)
eOe ra s +
G a~2
Hence, the equations (2.24) and (2 .2 7 ) become the same equation as follows:
Ar
1
=
x|«+ i / 2 ,y,fc
£ ] n+
^ r x s E
1+
cr Ar
,n
_
x lr-t-1/ 2 , / , *
—
1
2
e 0 £ 'rxs
1+
f j . n +
< rA r
- — —
° £ rXS
1
1/ 2 ,/,£
—
° E r“
(2.40)
Ar
+ _ f o £ j2 _ r 8
1 -i
„
,« + I/2
a rAr v y zl / + i / 2,y,/fc
—
2£0e'r
.
„
' “ rl/+ 1/ 2, /,*.
*+ 1 /2
z ■yl/+ 1/ 2 ,
16
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1
,n + l
'y U,j+ 1/2, *
(J A t
1+
A/
<v^
2sne'
o ry*£ tn
=
1
808'
t s i — f j \ n+i
y \ i , j + 1/2, k 2
1+
- 2 ~.
2ene'
0 ryA
+ j r
„ vA t ^ ^ lt, y' + 1 /2 , *
C
1
y\i, j + 1/2, k)
-
2e0s' rys
(2.41)
At
62 i 2 2 _ ( s f f l" t I / 2
- S J T J " 4- 172
C ^Af v ^
if, /+ 1/2, *
^
!/', j + 1/2, *,
1+
2eoeV
crTA f
1
n+ 1
t, y ,
At
z
1
,n
zl/, j, * + 1 /2 2
2 e 0 e ' r zsE
k+
a .At
1 /2
1+
808~ r q
1+
2e0e'rq
CT A f
f j .n f t
n
V z li, / , * + 1 /2
>
z li, / , * + 1 / 2 /
2808'rzA
(2.42)
A/
808 rzj
A
a ,A f
. s ^ r + 1/2
(.
'v l t , y , f c + i / 2
1+
.n + l/2
^
- 5 V# X
y x •f. j, * + 1/2/
2 8 0 e ’r z ,
and the original Yee’s finite difference expression can be derived from the above equations
and equations (2.13) to (2.15).
2.2 The Yee’s Algorithm
Using the second-order central difference approximation to space derivative in (2.15)
and (2.42) and no current source inside the computational space then the Yee’s
representation o f M axw ell’s equations are
.n + l/2
_
U,j + 1 / 2 , * + 1 / 2
At
fE
H
.n
1 /2
li, j + 1 /2 , * + 1 /2
y\j,j+ 1/ 2,* + 1 /2
-E
y\i, j + 1/ 2, *
Az
-E.
z\i,j+ 1, *+ 1 /2
zk /,* + 1 /2
Ay
17
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.43)
_
.n+l/2
.72—1/2
y \ i + 1/ 2, j , k + 1 / 2 ~
y \ i + 1/ 2 , j , k + 1 / 2
f+ l,/,fc + 1/2 ~
xlt + 1/2, j , k + 1
Ax
'y\
ti
z ,It, /, k + 1/2
E z
=
Az
-E Jn
x \j + 1/2, /, k
Ay
E \n
xU + l/2 ,j+ l,k
\
,n+ 1
Wrxs
-E ,
x\l+l/2,j,k
—
1+
2 ^ 0 £ 'rxs
t,72+1/2
n+l/'i
Jn + l / 2
’ t'+ 1 / 2 , / + l/2,fc
G At
•
2 e 0 B rxs
V
n+ 1
y lt,/+l/2,Jfc
z \ j + 1/ 2, j , k +
eoeV
+•
y\i,j+l/2,k'
G vA t
1+
2ene'
0 rys
t ,/+ 1/2,fc+ 1/2
G ,A t
Ar
z
2 e 0 e 'r z * £
CT Ar
1 +
eOB’rz,
a,A r
,72
z \ i , j , k + 1 /2
-■ Z ■
’
1 +
2 bobW
■*" i / 2
rr
(2.47)
.72+1/2
TT ,72 + 1/2
-H A
^
gl t + l / 2 , / + 1/2,*
z \ i - 1 /2 ,/+ 1/2,fc
Ax
Az
72 + 1
y-
2 e 0 e rys
.72+1/2
xlt,/+ l/2 ,fc -l/2
1 /2
CH x\. ,
z|i , / , * + 1/2
^
zl t + 1/2,/, &-1/2
A/
1
G vA t
1-
.7 2 + 1 /2
1/2
Az
2eoe n?s£ \n
1+
(2.46)
i« + l/2
ZU+ 1/2, / —1/2, k
Ay
i _
(2-45)
1, / + 1/2, k
At
G rA t
1+ —
y \i +
Ax
^ e OB r r J i-, ,n
x \i+l/2,j,k
—E ,n
y\i+\,j+l/2,k
a xAt
- TT
k
i+ 1 /2 ,/+ 1/2, £
fEJn
/■TT In +
(2.44)
x|t + 1/2, j ,
n—1/2
+ 1 /2
i+ 1 /2 ,/ + 1/2, fc
(H
-£
Z
(2.48)
2e0e'rz,
YJ ttl +
,/Z + 1 / 2
CtLA
—/i„
y l i + l / 2 , / , f c + l / 2 -vlt - 1/2, /, k+ 1/2
Ax
rz
1 /2
x lt ,/+ l/2 , fc+ 1/2
- e j ” + 1 /2
x lt,/ - 1/2, k + 1/2
Ay
18
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
^
2.2.1 Stability Condition
The stability condition is required to avoid numerical instability in fin ite difference
approximation schemes. The stability condition o f Yee algorithm is
first correctly
presented by Taflove[3] and is
1
Ar <
(2.49)
where Ar is the time step and A x , A y, and Az are mesh sizes along rthe x, y, and z
directions, respectively. The (2.49) is also called the Courant-Friedrichis-Lewy (C FL)
condition. Note that the CFL condition is derived by assuming a homogenous spatial
region. Generally, the CFL value is defined as follow:
(2.50)
which is less than or equal to 1.
An exact stability condition for general case is usually difficult to ckerived since it
depends on numerical boundary conditions, variable and unstructured meshing, and
material properties. However, substantial modeling experience has shown that numerical
stability can be maintained for many thousands o f iterations with the proper choice of the
time step. In the practical modeling, (2.49) is usually used. I f the num erical computation
diverges then a smaller time step is used, and so on.
2.2.2 Numerical Dispersion
The phase velocity of numerical wave modes in the F D T D grid can d iffe r from the
vacuum speed o f light, in fact varying with the modal wavelength, th « direction of
19
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
propagation in the grid, and the grid discretization.
The dispersion relation for a plane wave in a continuous lossless medium is simply
^2
= k x + ky + k l
(2-51)
where c is the speed o f light, co is the radian frequency, and k K, ky , k, are wavenumbers
along the x, y, z axes in this medium.
The numerical dispersion equation for FD TD scheme can be obtained by substituting
the plane monochromatic traveling-wave trial solutions into
the finite-difference
implementation o f M axw ell’s equations. The dispersion equation[16] of full threedimension Yee algorithm is
where kx , ky , and kz are wavenumbers along the x, y, z axes in the computational space.
Assume Ax = Ay = Az = A , choose C F L = a , and define R = X / A which is the
number o f grid cells in one wavelength, then (2.52) can be rewritten as
3 f . ( an. \ ] 2
IjsJ J
. ('fcsinGcostfA2
=
— J
. f&sinGsindA2
. fitcosG ')2
+sinl —
J
<2-53)
wherefc = X k / 2 which is equivalent to vp/ c = n / k where vp is speed of wave in
computational space and k
= kx + k y + k z . The 0 and <j) are polar and azimuthal
angles in the spherical coordinate system.
Several conclusions can be observed from (2.52) and (2.53) and are summarized as
follow.
20
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(1)
The Yee’s FD TD scheme gives a phase error. The speed o f wave in the computa­
tional space is less than that in free space. However, (2.52) and (2.53) are indentical in the lim it as A r, A x , A y , and Az all go to zero.
(2)
The smaller the C FL values, the larger the phase error from Figure 2.2.
(3)
The larger the R, the smaller the phase error from Figure 2.3. The effect on reduc­
ing the phase error due to change o f R is about 100 time that due to change o f
C FL values.
(4)
There is a numerical phase velocity anisotropy in Yee algorithm or other FD TD
schemes from Figure 2.3 and Figure 2.4.
(5)
The number o f grids in one wavelength has a lower bound that makes the numer­
ical phase velocity goes to zero and the wave can no longer propagate in the
F D T D grid from Figure 2.5.
(6)
The numerical dispersion occurs because the higher-spatial-frequency compo­
nents o f wave propagate more slowly than the lower-spatial-frequency compo­
nents. This numerical dispersion causes pulse broadening due to the spatial lowpass characters in FD TD schemes.
(7)
The numerical dispersion can lead to spurious refraction o f propagation numeri­
cal modes if the grid cell size is a function o f position in the grid.
21
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
^
0.998
0.996
.+ + + .
0.994
>+++'
0.992
+x
a=l
a=0.8
0.99
a=0.4
0.988
xxxxxxxxxxxxx
a=10e
0.986
+x
x+
0.984- ■+
0.982
20
40
60
80
100
120
140
160
180
0
Figure 2.2 Variation of the numerical phase velocity w ith
C F L values at R = 10 and (j) = n / 4 .
22
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
r= ----------- -----
_ _ J
-1
4-
/■
~4-
+
4* .
4-
0.98
"\
s
4-
4. '+•
“ n.
4-
+
...... t
4-
+
4+
4-
4-
0.97
4-
D --CA
+
■
R=20
R=10
+■
0.96
—
.......4 -
" — -------
V_
4-
4-
4-
4-
------ * > ' -
+
0.99
+ + + -N
+
..
e.
'>■
r z r r .l — -= n
'
X"
+ -H -H -++-H-4
+
444*
R -?
-h
44-
0.95
+
4-
4-
4*
.4-
0.94
+ 4-
20
40
60
80
100
120
140
160
180
0
Figure 2.3 Variation o f the numerical phase velocity w ith R
a t a = 0.9 and <{> = n / 4 .
23
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.998
0.996
0.994
0.992
0.99
0.988
0.986
0
20
40
80
60
100
120
140
160
180
0
Figure 2.4 Variation o f the numerical phase velocity
w ith $ at a=0.9 and R=10.
24
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.95
0.9
0.85
0.8
0.75
0.7
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Figure 2.5 Variation of the numerical phase velocity
w ith 1/R at a-0.9 and <J>=7t/4.
25
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.5
2.3 The Ty(2,4) (FD)2TD A lgorithm
From the analysis in section 2.2.2, the p-hase error in Yee algorithm keeps this scheme
from the calculation E M fields o f an electrical larger object or from applications that need
more accuracy, such as the phase cancellation technique. In the cavities with FEC
boundary' problem, the phase error and the error from surface impedance boundary
approximation are two main errors in the F D T D formulation. Hence, the fourth-order
spatial derivatives is employed to reduce th e cumulative phase error. For fourth-order
spatial scheme, special boundary treatment and degraded stability condition are the two
main disadvantage. The degraded stability condition is not very significant since a smaller
time step is usually used in F D T D schemes. The larger stencil on the F D T D mesh is very
troublesome when dealing with material discontinuities. However, the im plicit staggered
spatial finite difference schemes used in this, chapter simplifies this problem by using two
field points and three field derivative points.
2.3.1 The Fourth-Order Space D erivatives
The fourth-order finite difference expression can be categorized as explicit schemes
and im plicit schemes as follow:
Explicit collocated scheme
8(Mf- - n —iij- L) —(«,- + 2 —“ i —2 )
L2A x
(2.54)
Explicit staggered scheme
_ 27( “ /+ 1 /2 ~ ui - 1/ 2 ) ~ ( “ <+ 3/2 ~ “ < - 3 / 2 )
dx):~
Z4Ajc
26
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.55)
Im plicit collocated scheme
dr)
Vox) i+
1
+dr)
V o x ) Ii - I
6
C
2 fdu
3Vdx
(2.56)
Im plicit staggered scheme
dr)
Vox) I +
1
+dr)
Vo x) t
Ui + l / 2 ~ “ t - l / 2
Ax
24
(2.57)
A ll above equations can be derived by Taylor expansion and the details are in Appendix A .
The compact finite difference scheme used in this thesis is the im plicit staggered schemes
o f (2.57).
2.3.2 Dispersion Analysis for Explicit Staggered Scheme
For simplicity, M axw ell’s equations in a normalized region o f free space with jx = 1 ,
e = 1 , a = 0 , and c = 1 are considered and can be obtained as
JVxV =
3^
(2.58)
where V = H + j E . Apply (2.55) to the left hand side o f (2.58) and the central
difference to the right hand side o f (2.58) and let
j{kxl Ax + kyJAy + kzKAz —OlnAt)
then the follow ing equations are obtained
27
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-59)
sin (co A f/2) r .
,27sin(£zA z /2 ) - sin (3£zA z /2 )
X
.
.
j
Af
x \i , j , k
24Ay
.2 7 s in (£ y A y /2 ) - sin(3£yA y /2 )
+J
y 'i,J,K
(2.60)
_
24 Ay
zh,J,K
,27sin(fczA z /2 ) — sin(3fczA z /2 )
s in (c o A f/2 )T. .
r----------------------------------------------- V'
H-------^------------ V
24Az
-Mf,/, a:
Af
^l/, y, /r
27 sin
(2.61)
A x /2 ) —sin (3 A x /2 )
-------------------- TAKx-------------------- V A, . j . k = 0
,2 7 s in (£ vA y /2 ) —sin(3fcvA y /2 )
~J
24Ay
x\i, j , k +
.27 sin (fcxA x /2 ) - sin(3£xA x /2 )
J
(2.62)
s in (c o A f/2 )„ .
24Ax
y\ i , j , k
Af
_ .
2lf,/, K ~
For non-trivial solution, the determinant of the above equations is set to zero and the
numerical dispersion relation is
r s in (o )A f/2 )~r
)“|2 _ |~
p27sin(£xA
27sin(fctA xr /2
/2 ) - sm(JA:xA
sin(3£xA x /Z
/2 )T
)'
L
cAf
J
L
7
r27
27 sin(£_yAy/2) —sin (3£yA y /2 ) ii2
+L
J
24Ax
~
r2
r27 sin(fczA z /2 ) —sin(3fczA z/2)"]
J +L
24Ay
7'
(2 '63)
J
24 Az
The numerical dispersion relation can be further reduced to
1728^ . ^
a n J ' 2 _ J^yg-^^sinOcostfij
s.^3fcsin0cos<|)j'
2
+
[ 2 7 s i n ( * s i n % s in < 1 ’)
r o_ . fkcosO']
r H
-
s ft/S fs in e a in + y
(2.64)
. f3kcos6\~\
-
R -J .
if A = Ax = Ay = A z.
From Figure 2.6, the phase error o f fourth-order scheme is much less than that o f Yee
scheme. The phase error o f fourth-order scheme at R =5 falls between that of Yee scheme
28
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
at R=13 and R=40 as shown in Figure 2.7. Hence, a coarser mesh can be chosen fo r
fourth-order scheme if the same phase error is required.
2.3.3 Dispersion Analysis for Implicit Staggered Scheme
For im plicit staggered scheme given in (2.57) with numerical boundary condition, the
scalar equations o f V x V can be rewritten in a matrix form . For d V z/ d y , the matrix
equation is
AX = B
(2 .6 5 )
where A is a (N + l) by (N + l) m atrix, N is the number of partition along the y axis,
rdvz
[dy
and B =
Ay
dvz
dvz
dvz
dvz
dvz
(2 .66)
j = o dy j = l dy 1=2
j=n
- 2 dy j = N - i dy ;= v _
where M is a (N + l) by N matrix and
B =
2
J
2
J
2
2N -5
2
V
z .
2N -3
r)
J
.
J
3 ^
It
.
V
(2 .6 7 )
z .
J
2N- 1
2
J
The elements o f matrices A and M depend on the coefficients in (2.57) and numerical
boundary approximation. The (2.65) can be rewritten as
X =
MB
Ay
where M = A 1M and A 1 is the inverse of matrix A.
Apply (2.59) to (2.68) and after some calculation, the j component o f X becomes
29
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.68)
where mJn is the (j,ri) element o f m atrix M . Other derivatives o f V can be obtained
sim ilarly and a system o f equations for Vx , V , and V z are setup. Setting the determinant
o f above equations to be zero and A = Ax = Ay = A z, the dispersion relation is
obtained as
/£sin9cos<j>(2/ —(2/i-f- t ) ) -i^
R
rIV — 1
,/£sin9siti<|>(2./-(2;n- l ) ) - i 2
R
(2 .7 0 )
-n = 0
rJV-l
jkcosQ{2K —( 2 n + 1 ))- j2
R
+
-n = 0
The (2 .7 0 ) depends on locations o f points and matrix M ; hence, the analysis o f this
equation is very complex. The phase error in im plicit staggered scheme is close to that in
the exp licit staggered scheme. Moreover, the phase error in im plicit staggered scheme can
be further reduced if the suitable optimization method is used. The optimization methods,
which make k close to n or v / c close to one while maintaining smaller R, can be
researched based on (2.70).
30
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.99
0.98
0.97
0.96
4th order
Yee scheme
0.95
0.94
0.93
0.92
20
40
80
60
100
120
140
160
180
0
Figure 2.6 Comparison of numerical phase velocity
between explicit 4th-order scheme and Yee scheme at
R=5, a=0.24 , and <(>=7e/4.
31
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
\
\ +*+++++I
\
+ * + .+ + .
0.999- -+-H++ + +
a.
i>
0.998 - ...............
0.997
0.996
0.995
0.994
R=5 (4th order)
R= I3 ( Yee)
R=15 (Yee)
R=20(Yee)
R=40 (Yee)
0.993
0.992
0.991
0.99
20
40
60
80
100
120
140
160
180
0
Figure 2.7 Comparison of numerical phase velocity
between explicit 4th-order scheme and Yee scheme at
a=0.24, and <J)=te/4.
32
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.3.4 Calculation of Derivative of E in iy(2,4)
The derivatives o f E along jc, y, and z are first considered. The term d E y/ d z case is
considered here as an example. From the Figure 2.8, the evaluation o f d E y/dz\^
which
is close to the boundary
by
Az/2
needs
values
of
d E y/ d z ^
,
dEydz\
, EJ
, and E vI
if (2.57) is used. However, d E / d z \
is
y
\k = 3/2
ylfc = 0
y\ k = l
y
Ifc = —1/2
outside the physical boundary, hence, the value of dEy/ dz\ ^
so is the value o f dEy/ d z \
has to be approximated,
i/2 where N . is the number o f partition along z axis. For
mesh points other than k = 1 /2 and k = N z — 1 /2 , the (2.57) is applied to yield the
derivatives of E .
Ey
Ey
dEy/ d z
k = -l/2
V
A
Ey
i1
k=0
k = l/2
k=Nz+ l/2
V
A
i
1
k=l
lc = l/2
k= N ,
Figure 2.8 The lattices o f Ey and dEy/ d z along the z axis at fix
x and y.
The following two fourth-order approximations are used to approximate the
derivatives of E on those two special points;:
33
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ld E y
1 3Ey
24 dz
,
k= k
5 *Ey
,
12 dz
,
* =1
24 dz
*=l
13 ^ y
12 dz
,
12 dz
E y \k Az
-,
(2.71)
ld E y
ldEy
5 ^ y
24 8z
k=i
EA
13^y
12 3z
, 24 dz
k = Nz- Z
E\
y \k = N .
k =
N.
(2 .7 2 )
-E J
y \lc = N . -
1
Az
A system o f equations are derived and written in matrix form,
AX = B
(2 .7 3 )
where
0
26
-5
4
-1
1
22
1
0
•
•
0
0
1
22
1
0
•
0
A =
X =
dz
(2 .7 4 )
0
0
1
22
1
0
0
0
0
1
22
1
0
0
-1
4
-5
26
dE,
dE,
dz
dz
dz
(2.75)
k = N.--
B =
24
Az
(2.76)
£ ^l/fc
vl = Mz- 1 - £ yvl\k = N z- 2
E\
y \k = N z
-EJ
y \k = N z- \
34
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
and T denotes the transpose o f a m atrix. B y using the L U decomposition, the derivatives o f
E at every points can be determined.
2.3.5 Calculation of Derivative of f t in iy(2,4)
The calculation o f d H , / d y is shown as an example. Figure 2.9 shows the grid
configuration along y direction.
dH /d y
dH /dy
dH /dy
Hz
i
1
j=-l/2
V
A
j= 0
j=Ny+l/2
v
A
j= l/2
j= l
j= N y
F ig u re 2.9 The lattices o f H , and d H / d y along th e y axis a t fix
x and z.
I f the boundary is PEC, then E x , E z , and H y are zero according to Figure 2.1. Hence,
the derivatives o f H z on the boundary are zero, also. The m atrix equations, A X = B
becomes
22
1
0
•
0
1
22
1
•
0
•
0
A =
0
1
22
1
0
0
1
22
35
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.77)
X =
1T
~dHz
dHz
dy
dy 7 = 2
7= 1
dHz
dHz
dy 7 = ^ - 2
HJ
1/ = 3/2
&
(2.78)
j = N y- l
-H .
z\ j = l / 2
H z\Ij = 5/2 ~ K z\Ij = 3/2
24
B = —
Ay
(2.79)
HJ
-H J
z\j = Ny- 3 / 2
z\ j = N y- 5 / 2
HJ
-H
z' j = Ny - 1/2
z\ j = Ny- 3 / 2
and
dy
=
0,
ay
7= 0
=
(2.80)
0 .
7 = AT,
I f the boundary is not PEC then d H z/ d y \
^ and d H z/ d y \
^
need to be
approximated. The following fourth-order one-way finite difference approximations are
used:
31h
~8 z
dH,
+ 229h
l + 24
1
1
2
37
5 + T f f *-
75 h
3 8~ z
2
J
2
k H zz
1 ~ 7Y2
7= :
Ay
dy 7 = 0
J
2
(2.81)
and
dH,
dy
7'= Ny
V
8
z
J
..
l
= Ny- 5
24
z • „
3
J~ y 7
j
4
8
h z
z j = Ny- (2.82)
37
-j^ H
/A y
+ i Hz j = Ny- ~7 1 2 z
j = Ny- ZJ
The components o f the corresponding matrix equation are
36
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
24
1
A
X =
dy
=
0
1
0
0
dHz
'd H z
j =
0
1
0
22
•
0
•
0
•
0
22 1
0 24
dy
j =
dy
N y- 2
1 T
an.
an,
d H z
0 dy j = i
(2.83)
j = N ,-
1
ay
(2 .8 4 )
y=Af,
and
37
31
~~2
i
z
24
k H zz
7 ~ l12
5+T " 2
z
1
1
2
J
2
2
J= :
-H .
j = 1 /2
z l / =» 3 / 2
j = 5 /2
24B =
Ay
-H „
Ij
= 3 /2
(2 .8 5 )
HJ
-H J
zl/
= /yy - 5 / 2
HJ
-H J
zl;
= Wy- 3 / 2
z l/= 7 V y - 3 / 2
1/ = A/y— 1 / 2
/
where
f
/
=
31 rr
j = N y-
3 - 48
» (2.86)
37
+
^ z
7
-J 5 H Z
12 z
/ Ay
J="y-V
37
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.4
Excitation Source and Power Analysis
2.4.1 Excitation Source
An excitation probe is placed parallelly to one of the rectangular axes. The current
density o f the probe is modeled by induced E M F method and represented as
J = ( / ( x ) S ( y - y 0) 8 ( z - z 0), 8 ( x - x 0) / ( y ) 8 ( z - z 0), 8 ( x - x 0) 8 ( y - y 0) / ( z ) ) v v ( 0
(2.8 7)
where
sink(ha —a )
sinfc/z,
a = x, y, z,
(2.88)
w( t) is the tim e dependence o f J and ha is the length o f probe along the a direction.
For empty PEC rectangular cavities, the source frequency is 2.45 GHz. In order to
contain a finite energy in the cavity, the signal source needs to be turned off at some given
time. The sidelobe level o f the Blackm an-Harris (B H ) window[20] is approximately -92
dB and provides a smooth transition o f excitation. The B H function is discretized as
follow,
0 .3 5 8 7 5 + 0.48829 cos
+ 0.14128 cos
+ 0.01168 cos
0
otherw ise
where w (r) = /z(r)coscor, tn is the time step, tc is the center o f the B H window, and
N halj- is the designed half width o f the B H window function.
38
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
For a loaded PEC rectangular or an aperture-coupled cascaded cavity, a single
frequency sinusoidal source is used. The loss o f material sample in the cavity accounts for
the finite value o f Q.
The excited modes are determined by the location of the probe according to
E ( P ) = j j ( P ) G ( P , P Q)dvQ
(2.90)
where the Green’s function is
^
- x
G (r, r„) =
.
^ E n(P0) E n(P)
2_ 2
n
and ^
Kn
(2-91)
K
denotes the sum o f all the T E nml and T M nm| modes. For example, to excite T M
n
modes, the probe must have Jz component only.
2.4.2 Power Analysis
The integral form of the Poynting’s theorem is
- j E ■J sdV =
V
JE - J d V + j f f ~
V
V
+ H ~ j d V + j>(E x H ) - hdS
s
(2.92)
-X
where J s is the source current. For lossless case, the power stored inside a cavity after the
source turned o ff is
<2-9 3 )
The FD TD approximation o f the first integrand, evaluated at tim e step n , in (2.93) is
39
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
K b = 2 ^ 1 " (D* f + W
‘ ') + E ,\" {D y \" * ' - D y f - ‘)
(2.94)
+ E 2|',( D £| " 1-Z )z| ' " ' ) j
and that o f the second integrand, evaluated at time step n, is
2
HJ
V
rr
H.
2
1
c
2A r
B+
V.
P s H -~
(
H x\ “
x\
K
1n2
,n + r
r-
I
r
+
_\
o
(
HJ
y\
n- h
2
HJ
y1
K 2 _,
(2.95)
_V
where the nonmagnetic m aterial is assumed. Hence, the power stored inside the cavity at
time step n is
P> 'Z ( P "
e
+ P"h)A v
(2.96)
where Av = AxAyAz for uniform partitions.
For material with conductive loss, the dissipated power is the first term on the righthand side o f (2.92). The conductive dissipated power density is o E
and its time domain
formulation is easy to determined. However, the dissipated power for a lossy dielectric
material is contained in the second term on the right-hand side o f (2.92) and the dielectric
dissipated power density in the tim e domain has to be extracted. F o r this case, (2.93) is the
summation o f the stored power o f the whole cavity and dissipated power due to the lossy
material sample. Hence, only the calculation o f the (2.93) and the dissipated power are
needed for power analysis. However, the stored power inside the loaded material region is
also presented below to form a general Poynting theorem.
40
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Assume the one-relaxation Debye medium, the Fourier transform o f the x component
- dD .
in E ■-r— is
dt
_L
(£
271:
t (CD)
® [y(0 8 '(a ))£ :t (co)] + Ex( ca) ®
[ ( D 8 " ( co) £ : ;c( co) ] )
(2.97)
where e'(co) and e”(co) are the real and imaginary parts o f electric perm ittivity and the
sign, ® , stands for the convolution. The first term in (2.97) is the stored power density due
to Ex and the second term is the dissipated power per volume due to Ex . Substituting the
relation o f the perm ittivity w ith static and infinity frequency perm ittivities, (2.97) can be
rewritten as stored and dissipated powers as
rJ E x{co) ®
[yw e’xooi?x(( 0 )] + E x( a ) ®
( £ ’xs ~ e jro o )/05
(2.98)
1 + (C 0Tex) 2.
and
1
27t
)r
r
E J a >)
t\
(2.99)
The corresponding time domain equations for above equations are
d E {t)
r
a ggAxi t) y
- [
d
+ * * (O L ( e '„ -
(2.100)
and
( e’x ,-e 'x ~ )
t E x{t ) E x( t ) - E x{t)gx(t ))
(2 . 101)
ex
where gx( t )
is the inverse Fourier transform o f isx(C D )/[l + (toxer) 2] . Hence, the
41
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Poynting theorem fo r nonmagnetic material can be rewritten as
\ E S - J d V = J[<y]£ • E d V + J ( K 3]E • E - [5 3]E - £ ) d V
vr
v
vt
(2.102)
+
^
+ [£2] £ • | | + H - | j ^
+ f ( £ x f f ) • fidS
where
e'^co
0
G il =
0
G 2] =
0
0
Z'yco 0
0
(2.103)
e’z~
e ' x , - e 'too
0
0
0
8*ys — c
8*yoo
0
0
0
e'<.S - e z°°
Gx 0
[a ] =
0
0 ay 0
0
(2.104)
(2.105)
0 a.
and
e
—ce .r« »
° XT
K 3] =
0
0
8 ys —8 y°
0
0
(2.106)
'e y
0
0
8
c zs —8
c z°°
ez
where the Vs is the volume of the material sample. Note that the summation o f the third
term in (2.102) can be obtained by subtracting the second term in (2.102) from (2.93).
42
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Evaluating the gx(t) involves the convolution of Ex(t) and the time domain version of
1 /[1 + (Q)Te)2] . To get the convolution, a history o f E x(t )
is required and this is
impractical in the program implementation. Hence, the convolution need to be
approximated by a recursive equation. The first method is to discrete gx(t) into a
recursive equation by discreting the convolution directly. By using the Fourier transform
pair, e~aM and 2 a / { a 2 + co2) , the discrete version o f gx(t) is
gx(nAt) = £ L E x(nAt) + - ^ - e ^ - E ^ n - l) A t ) + e- * t/T'*gx« n - l)A r ).
^ Tex
ex
(2.107)
Hence, the stored power density due to Ex at time step n is calculated from the discrete
version o f (2.100) and (2.107) and that of dissipated power from discrete version o f
(2.101) and (2.107).
The other method uses a procedure similar to what we did to Debye equation in
section 2.1.1. Let gx(co) = E x(c o )/[l + (coxex) 2] or [1 + (coxer) 2]g x(co) = E x( co),
then its time domain expression is
2 d2g (t)
S x it) - Tex— ^ ~
= E* (0 •
<2-108)
B y using the second order central difference approximation o f the second derivative,
(2.108) becomes
gx((,z + l)A /) = [2 + ^
2]g x(n A r )-g x( ( n - l ) A r ) - ^ 2£:x(n A r).
(2.109)
The (2.107) is obvious first order accurate but (2.109) is second order accurate in time.
However, the homogenous solution o f (2.108) contains a divergent term, et/r 'x, and its
43
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
FD T D formulation cannot be used. M ore general F D T D form ulation o f power for a
general Debye medium w ill be discussed in the next chapter.
2.5 Prony’s M ethod
Results in the frequency domain are usually obtained by recording the time-domain
response at the selected observation points of the F D T D mesh and applying the FFT
algorithm. There are, however, several limitations in the F F T approach. The main
limitation is that o f the frequency resolution A /, which is roughly the reciprocal o f
observation time, i.e., A / = 1 /(A /A r) where At is the tim e step and N the total number
o f iterations used in the F D T D method. A second lim itation arises from the windowing o f
the time-domain data because the F D T D response is truncated in tim e. This has the effect
o f viewing the true tim e-dom ain response through a rectangular windows o f duration
T = N A t . In the frequency domain this windowing is translated into the convolution o f
the true spectrum with the function s i n / / / .
The Prony’s method is a technique for modeling sampled data as a linear combination
o f complex exponentials and is particularly suitable for calculating the resonant frequency
and Q o f a resonating structure, since the impulse response o f such a structure is
characterized by a superposition o f decaying exponentials[21].
2.5.1 Theory
The Prony’s method fits an exponential approximation w ith unknown exponents of the
form
44
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
f ( m ) = C 1eAl'n + C 2e A2/n + - " + C j2<’Aem =
C^k
^
(2.110)
k= i
to a function f ( t ) by sampling at n equally spaced points and \i k = eAk. A linear change
of
variables
has
been
introduced
in
advancesuchthat
m = 0, 1, 2, ...n — 1 . The coefficients Ck and
the
data points
are
can be complex numbers. The above
equation can be written as follow,
C l + C2 + --- + C p = f Q
C 1p 1 + C2(x2 + ... + Cpy.p
C2\ i \
C jH i +
C iP i
+ C 2p2
+
... +
c p\kp
=f x
( 2 . 111)
= f 2
+ . . . + CpPp = f n - l
or in the m atrix form
1
1
Hi
h
2
Hi
1
1
2
h 2
1
1
f
n —1
2
n- 1
Hi
h
i
2
...
:
- -
Hp
-
2
Hp
C i
h
3
h
2
3
I
1
h
3
C2
c 3
- H
p
_
where /,- = f ( i ) and i = 0, 1, 2, ..., n — 1. Let
f i
f 2
1
I
1
1
i
i
_C
f o '
=
i
1
1
i
71— 1
n —1
'
/ n- I
p_
..., \x.p be the roots o f the following
equation
p p + a 1p /," l + a 2p p _ 2 + . . . + a p_ 1p + a p =
p
j a *p * = 0
k= 1
then we can get the following equations from (2.111) and (2.113)
45
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2 .U 3 )
/ p + / p - l a l + / p - 2 0t2 + - " + / o a p = 0
/p + I + / p a l + / p - l a 2 H
+ / l a p =
/ n - I + / n - 2 « l + / „ - 3 a 2 + " * + / n - p - l « p
0
=
0
or the follow ing one by shifting the first terms on the le ft hand side to right side o f
equal sign
(2.115)
which is linear in a .
The Prony’s method solves the (2.110), a nonlinear equation in f i, by solving a , in
(2.114), a linear equation in a , finding the roots, p.,-, in (2.113), and obtaining C { in
(2.111) where i = 1, 2, 3
The original Prony method chooses n = 2p + 1. The original Prony method works
perfectly w ell when no noise (truncation error or measurement imprecision) is present in
the data. However, when a noise is introduced, then this method performs very poorly,
largely due to the extreme sensitivity o f root locations to the coefficients o f the polynomial
in (2.113). The least squares Prony method uses n > 2 p + l , then (2.112) and (2.115)
become overdetermined linear equations systems and the singular value decomposition
(S V D ) is used to solve those equations.
Determ ining p is an important issue o f applying Prony method to frequency domain
46
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
analysis or prediction since every pt- represents a frequency component. I f the value o f p is
less than the number of actual modes excited in the structure, the spectral resolution is
poor[21]. However, if p is selected to be too high, nonphysical modes appear. Nonphysical
modes introduced by Prony estimation can be suppressed by the application o f Prony’s
method w ith the time sequence o f the sample points in reverse order[22].
2.5.2 Estim ation of Resonant Frequency and Q uality Factor by Prony
M ethod
From the definition o f quality factor,
_
S to re d Energy
2 “ m° P o w er loss
the F D T D time-domain response, say E
(2116)
o f a cavity can be expressed in terms o f a
superposition o f the resonant modes [24]
E ,(m A O =
£
=
k= 1
£
^
(2 H7)
k - I
with a k = a k/ ( 2 Q ) where m = (0, ...,
) N — 1. By Prony method, the Ck and
can
be obtained if E x(mAt) are known. Hence, the resonant frequency and damping factor are
im agClnCp*))
*
2 S iJ ”
(2 ' 118)
and
re a l(ln (p ^ ))
a k = ---------- — --------
(2.119)
where \ik is the calculated mode corresponding to f k , imagClnCp.^)) is the imaginary
47
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
part o f In ^ ^ ) , and re a l(ln (p fc)) is the real P31* ° f ln (P,fc) •
Once the damping factor and the resonant frequency have been determined, the quality
factor can be easily obtained as
Note that (2.120) is obtained when the cavity mode is at a steady state[23]. However, the
Prony method provides a way to estimate quality factor o f cavities from the field
distribution at transient state.
2.6 A Single Empty Cavity with PEC Walls
An em pty rectangular cavity with PEC boundary is studied in this section. The Q value
o f this cavity is infinite so the time to achieve the steady state is almost infinite. Hence, all
the results presented here are obtained in the transient state.
2.6.1 Configuration
The configuration of the empty cavity with PEC walls is shown in Figure 2.10. The
excitation probe is located at one o f the six faces o f the rectangular cavity and its length is
assumed to be very small. The dimensions o f the cavity are: bl\b2xb3.
2.6.2 Num erical Results of Field Distributions
For T E 011 mode, the excitation probe is located at the center o f the y-z plane and only
J x at 2.45G H z is provided. The B H source is used and automatically turned o ff at 18,875
time steps, spanning over 0.2516|i.s, w ith At = 13.3299 ps. The cubic cavity with bl, b2,
and b3 a ll equal to 0.08658 meter is assumed. For this dimension and frequency, the
48
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
T E 101’ T E o ib and T M 110 modes can be excited; however, only T E q h can be excited
because only J x is provided on the probe and Exs o f other modes are zero. The numbers o f
partitions along the x, y, and z are all 10.
The x, y, and z dependences o f
are shown in Figure 2.11 to Figure 2.15. For TEqU
mode, Ex is constant along x at fixed z and the numerical results are shown in Figure 2.11
and Figure 2.12. The
is proportional to sin ( tc c c /0 .08658) along the a axis where a is
y or z and are shown in Figure 2.13, Figure 2.14 and Figure 2.15. Note that the numerical
results at those grid points are exactly equal to D s in (7 ia /0 .0 8 6 5 8 ) where D is a
am plitude factor if a is from 0 to 0.08658 with increasing 0.008658 for every step. The
summations o f Ey, Ez, and H x over all grid points are 3.7386e-06, 1.9023e-06, and 0
respectively and it shows that no Ey, Ez, and H x exist.
For T M m mode, the probe is located at the center on the bottom x-y plane and only Jz
is provided since Ez of T E L11 modes is zero. The length o f the cubic cavity is 0.10604
meter and the B H source is o ff at 0.25131ns, with At = 16.3265 ps. The electric field
distributions are plotted from Figure 2.17 to Figure 2.21, and they are consistent to those
from theory[19]. The summation o f H z over all grid points is -5.7932e-31 which is almost
zero, thus confirms the excited dominant mode being a T M mode.
The tim e dependence o f £"x for TE 0n mode at given point is plotted in Figure 2.23 and
its frequency response in Figure 2.24. From Figure 2.24, the resonant frequency of T E 011
mode is observed to be close to 2.5 G H z. However, the estimated resonant frequency o f
T E q h mode by Prony’s method is 2.4528 G H z which is much closer to the source
frequency. The estimated resonant frequency o f T M L11 is 2.4548 G H z by Prony method.
49
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Probe
b3
Figure 2.10 The configuration of the empty rectangular cavity
with PEC boundary. The excitation probe is located on one of the
six faces.
50
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1
y=0____________________________________ y=0.08658
0
-1
-2
^=0.008658 ______________________________y=0.0Z7922_
-3
-4 ■
-5 ■
. _3RQ.Pmi6______________________________ y=0.Q69264_
-6 ■
-7 -
H-------- ' ~l-------
1------------- 1-------------1— ----------1------------- 1--------- 1—r_------ —-t--------------1-
-8 J(
-9 -1 0
-
0.01
v=0).p34632 K
y=0.04329
0.02
)(
0.03
u
0.04
x(m)
ic
0.05
jt
0.06
V=p.051948Jt
y=0.04329
0.07
0.08
Figure 2.11 The x dependence of E x a t z= 0.025974 and t= 0.26658 fis
for T E 011 mode.
51
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
K
0.09
“i
i
i
i
i------------1----------- 1--------
_______ z=0________________________________
z=0.08658
-2Z=0.008658
z=0.077922
Z=0.017316
z=0.069264
x
tii
-8
H--------- ---------- -r---------1----------1---- ------1---------- 1-------- x-------- 1----------- h
-1 0
]£
£ = 0 .0 3 4 6 ^
z=0.04329
j(.
u
it
ic
Z=$.051948)(
Z=0.04329
)C
-12'----------- 1------------1----------- 1----------- 1------------1_______1_______i_______i0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
x(m)
Figure 2.12 The x dependence o f E x at_y= 0.034632 and t= 0.26658 p,s
for T E 011 mode.
52
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.09
X
UJ
+•
x
•12
0
0.01
0.02
0.03
0.04
0.05
0.06
Z=0
z=0.008658
z=0.017316
z=0.025974
z=0.034632
Z=0.04329
0.07
0.08
0.09
y(m )
Figure 2.13 The y dependence o f E x a tx = 0.025974 and t= 0.26658 (is
for T E 011 mode.
53
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-2
x
UJ
-6
-8
y=o
y=0.008658
y=0.017316
y=0.025974
y=0.034632
y=0.04329
-10
-1 2
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
z(m )
Figure 2.14 The z dependence of E x a t x = 0.034632 and t= 0.26658 jxs
for T E 011 mode.
54
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.09
LU
—
8
10
-12
0.1
0.08
0.1
0.06
0.08
0.06
0.04
0.04
0.02
z (m)
0.02
0
0
y(m)
Figure 2.15 The E x variation along th e jz plane at x= 0.034632 and t 0.26658 ps fo r T E ()11 mode.
55
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
-
0.01
-0.01
0
0.01
0.02
0.03
0.04
Figure 2.16 The E field on thexz plane a t
fo r T E q h m ode.
0.05
0.06
0.07
0.08
0.09
0.034632 and t—0.26658 p.s
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
X
Ui
y=o
y=0.01064
y=0.02128
y=0.03192
y=0.04256
y=0.0532
0
0.01
0.02
0.03
0.04
0.05
0.06
x(m )
0.07
0.08
0.09
0.1
Figure 2.17 The x dependence o f E x a t z= 0.03192 and t—0.31019 ps
fo r T M j j j mode.
57
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
y=0.0532
y=0.04256
y=0.03192
y=0.02128
y=0.01064
y=0
Hi
0
0.01
0.02
0.03
0.04
0.05
0.06
x (m)
0.07
0.08
0.09
0.1
Figure 2.18 The x dependence of E y a t z —0.04256 and t= 0.31019 fis
for T M m mode.
58
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
14
12
fsl
UJ
Z=0
z=0.01064
z=0.02l 28
z=0.03192
z=0.04256
-6
Z=o!(3532
0
0.01
0.02
0.03
0.04
0.05
0.06
x(m)
0.07
0.08
0.09
0.1
Figure 2.19 T h e x dependence o f E z a t j = 0.04256 and t= 0.31019 J4s
fo r
m ode.
59
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.12
0.1
/
^
0.08
0.06
0.04
0.02
-
0.02 L —
-
0.02
0
0.02
0.04
0.06
0.08
0.1
Figure 2.20 T h e E field on the xy plane a t z= 0.04265 and t= 0.31019 ps
for T M m mode.
60
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.12
0.12
0.1
0.08
0.06
N
0.04
0.02
-
0.02 ■—
- 0.02
0.02
0.04
0.06
0.08
0.1
x
Figure 2.21 T h e E fie ld on the xz plane a ty = 0.04265 and t= 0.31019 jis
for T M m mode.
61
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.12
0.12
0.1
0.08
0.06
N
0.04
0.02
-0.021—
-0.02
0
0.02
0.04
0.06
0.08
0.1
Figure 2.22 Th e E field on theyz plane a t x = 0.04265 and t= 0.31019 (is
fo r T M m m ode.
62
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.12
20
x
Hi
-5
-1 0
-1 5
-20
0.247
0.2475
0.248
0.2485
0.249
0.2495
0.25
time
Figure 2.23 The tim e variation o f Ex at x=0.04329, y=0.04329, and
z=0.034632 for T E 011 mode.
63
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.2505
|4S
2500
Spectrum
of Ex
2000
1500
1000
500
0.5
2.5
1.5
:
3.5
Frequency in GHz
Figure 2.24 The frequency response o f Figure 2.23.
64
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4.5
2.7
A Lossless Loaded Cavity with PEC Walls
In this section, the numerical computations o f lossless loaded cavity with PEC walls
are considered. A lossless rectangular material sample is placed in the center of the
rectangular cavity and the dimensions o f the rectangular cavity are shown in Figure 2.25.
The excitation probe is located at the center o f the xz plane and is very short comparing to
the dimension o f the cavity. The excited field is TE10I mode and the operation frequency is
2.45GHz with the wavelength k equal to 0.12245m when there is no material sample
present. In order to compare w ith theoretical estimations, three m aterial samples with
selected shape and dimensions have been studied. The excitation source is off
automatically and all the numerical results are based on transient state solution.
2.7.1 Quasi-cubic Case
A quasi-cubic material sample which has almost equal dimensions is placed in the
center o f the rectangular cavity. The dimensions o f the material sample are set to be
x0=0.00343m, yo=0.0034m, and zo—0.00352m. This sample is assumed to be lossless and
have the relative perm ittivity o f er = 2 .5 . In this computation, the number o f partitions
along the x , y, and z directions are 21, 20, and 33 respectively. The material sample is
located at node 10, node 9, and node 16 along the x, y, and z directions where the number
o f nodes starts from 0. The excitation source spans over about 67,412 time steps with
At = 3.72822ps and the computation stops at 68000 time steps.
W ith this electrically very small sample, x 0/ k = 0.028, the static electric field
induced inside o f a dielectric sphere, E ws = ( 3 / ( 2 + er) ) E ns, is used to estimate the
induced electric field in the sample where
e ws
is the electric field inside the dielectric
65
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.072m
_
A
xq |
A
/
- - z o - - 4 ----------------
0.1163m
0.034m
y
Figure 2.25 Dimensions o f the rectangular cavity and the loaded m aterial
sample. The center of the m aterial sample is consistent w ith the center of the
cavity.
66
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
sphere and E ns is that in the region of sphere before the dielectric sphere is placed. For
TE joi mode, only Ey exists since Ex and Ez are zeros for this mode.
The variation o f E y along y axis at x —0.034286m, z=0.056388m and at f=0.25352ps is
plotted in Figure 2.26 and the calculated ratios o f E ws/ E ns inside the sample are 0.6550
and 0.6352 at node 9 and node 10, respectively, and the electrostatic estimated ratio is
0.6667. E ns is approximated by the electric field at first node because Ey is constant along
y axis in TE jqj mode when there is no sample present. The closeness o f the numerical
results and the electrostatic estimation gives confidence to the numerical accuracy. The
ratios of E ws/ E ns along the y axis at several different times are shown in Table 2.1 and
those ratios are almost independent with time for the lossless case.
Table 2.1 Ew
n^ ratio at different times
rvj / E flS
E ws / E
at node 10
Tim e Step
E ws/ E ns a t n o d e 9
67999
0.6550
0.6352
67979
0.6550
0.6352
67959
0.6550
0.6352
67951
0.6550
0.6352
ns
Due to the induced charge on the material sample surface and the induced current in the
material sample, the other components o f the electric field are induced to satisfy the
boundary conditions. The E x, Ey, and Ez along the y axis are shown in Figure 2.26 and
67
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
shows that the excited cavity mode is not a pure TE iqi mode anymore since there are Ex
and Ez inside the material sample. However, the
T E jq j
mode still dominates inside the
material sample judging from the amplitudes o f Ey versus Ex and Ez in Figure 2.26.
The estimated resonant frequency o f T E 10i mode for PEC empty cavity with the
dimensions shown in Figure 2.25 is 2.4532 G H z w hile that for PEC loaded cavity with a
quasi-square cubic sample is 2.4478 G H z. The resonant frequency decreases about 0.22%
after the material sample is placed inside the cavity.
2.7.2 Thin Square Plate Case
The material sample with a shape o f a thin square plate, having its height much smaller
than its width, is placed in the center o f the cavity. The numbers o f partition in this FD TD
calculation are 15x17x10 and the dimensions o f the material sample are 0.024m, 0.002m,
and 0.02326m along the x, y, and z directions. This sample is also assumed to be lossless
with the relative perm ittivity of sr = 2.5.
The x dependence o f Ey is plotted in Figure 2.27 and then a sine function is used to fit
that curve. B y this way, the E ns is obtained since the Ey versus y plot is not a constant
anymore. The induced electric field inside the material sample is estimated by the
boundary condition o f E ws = (1 / e r) E ns. The ratios o f E ws/ E ns along the x direction for
different locations o f z are plotted in Figure 2.28. Those ratios inside the m aterial sample
varies from 0.39 to 0.45 which are close to the electrostatic estimation o f 0.4. The ratios of
E ws/ E ns varies slighdy at different times for this case.
68
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.6
x=0.034286
z=0.056388
0.5
Electric fields
0.4
0.3
Ex
0.2
Ez
0.1
-
0.1
0.005
0.01
0.015
0.02
0.025
0.03
y (m)
Figure 2.26 The variation o f electric fields in the y direction at
t = 0 .2 5 3 5 |Ij.
69
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.035
0.16
0.14
0.12
0.1
S '0.08
0.06
0.04
Calcualted Ey
Fitted Ey
0.02
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Figure 2.27 The variation o f Ey in the x direction at t = 0.24289fis .
The line w ith star symbol is the calculated values and the line with
circle symbol is the fitted values fo r the empty cavity.
70
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.6
0.58
0.56
0.54
0.52
0.48
0.46
0.44
0.42
0.41—
0.02
0.025
0.035
x(m)
0.03
0.04
0.045
0.05
Figure 2.28 Variations of Ews/ E ns in the x-directions. Each curve
represents this ratio as a function o f x for different locations of z. The
relative perm ittivity of the thin square plate m aterial sample is e r = 2.5.
The solid line with symbols is the ratios at t = 0.24289p.5 and the dash
line w ith symbols is those at t — 0.24301 (is.
71
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The excitation source is o ff at 0.239ns which covers 49,049 tim e steps with
At = 4.862ps and the computation stops at 50,000 time steps. The calculated resonant
frequency is 2.4382 G H z and the frequency shift is 0.61% comparing to 2.4532 G H z
which is the resonant frequency in empty cavity. The maximum o f Ex and Ez inside the
material sample is in the order o f 10"4 which is very small comparing to the Ey This is
because the m aterial sample is very thin in the y direction so that there is no significant
induced field Ex and E z.
2.7.3 Narrow Strip Case
The dimensions o f this narrow strip sample are 0.002322m, 0.02125m, and 0.002115m
along the x, y, and z directions. This sample is also assumed to be lossless and have the
relative perm ittivity o f er = 2 .5 . The number o f partitions o f this F D T D calculation is
31x16x55.
Theoretical estimation o f the induced electric field in the m aterial sample may be close
to the electric field when the cavity is empty because the in itial electric field is tangential
to the m ajor part o f the material sample surface, and the continuity o f the tangential
component o f the electric field at the m aterial sample surface requires this estimation. The
variation o f electric field along y is plotted in Figure 2.29 and the ratios o f E ws/ E ns are
shown in Figure 2.30. The induced Ex and Ez fields inside the m aterial sample are almost
zero and Ey is the dominant field. The ratios o f E ws/ E ns inside the m aterial sample range
from 0.94 to 0.99 w hich are close to the theoretical estimation o f 1.
The excitation source is o ff at 74,835 time steps which spans over 0.25133iis with
72
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
At = 3.35841 p s , and the computation stops at 75,000 time steps. The calculated resonant
frequency for this narrow strip case is 2.4464 G H z and the frequency shift is 0.28% from
the resonant frequency o f the empty cavity o f 2.4532 GHz.
73
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.7
0.6
0.5
Electric fields
0.4
x=0.034839
0.3
- Ex
0.2
0.1
-
0.1
0.005
0.01
0.015
0.02
0.025
0.03
y (m)
Figure 2.29 The variation of electric fields in the y directions
a t x=0.34839m, z=0.057093m, and t=0.252|xs.
74
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.035
1.25
1.2
▲
1.15
x=0.034839
z=0.057093
1.05
0.95
0.9
0.85
0.8
0.005
0.01
0.015
0 .0 2
0.025
0.03
y(m)
Figure 2.30 The ratios o f Ews/ E ns in y directions.
75
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.035
2.8 A Lossy Dielectric Loaded Cavity with PEC Walls
The field distributions o f a cavity loaded with a lossy dielectric material sample with
PEC w all is discussed in this section. The excitation source for this lossy case is a single
frequency sinusoidal source turned on all the time. The size o f material sample is chosen
large enough to cause a relatively low quality factor o f the cavity. The relation o f cavity
quality factor and number o f the time steps which is needed to reach the steady state is
also studied in this section. In order to identify the tim e steps to reach the steady state, the
FD TD formulation o f power analysis is used.
2.8.1 Configurations
The physical configuration o f this case is the same as that in Figure 2.25 but with a
larger material sample. The numbers of partition along x, y, and z are 1 5 ,1 7 , and 10 and
the dimensions o f the material sample are: 0.0336x0.014x0.0698 m.
The real and imaginary parts o f the relative perm ittivity in (2.3) can be obtained as
e'r(co) = s'r„ + e -* ~ £ r^
(2.121)
1 + (C D T e )
e " r (CO) =
(e - "
e r°°)™ Z e
(2.122)
1 + (G )T e )
and the s’r and e"r can be expressed as follow:
e'rs = e’r(Q)) + coTee"r
(2.123)
(2.124)
e
For the lossy m aterial sample, the relaxation time, xe , is not zero and it depends on the
76
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
properties o f the m aterial. In this FD TD computation, xe is assumed to be 10~9 seconds
and the e'r is set to 2 .5 with four different e"r as listed in Table 2.2.
Table 2.2 Permittivities mapping at xe = In s , to = 2 n (2 A 5 e 9 )
e'r(co)
e"r(GJ)
e 'r ,
e' r©o
2.5
0.1
4.0394
2.4935
2.5
0.5
10.1969
2.4675
2.5
2.5
40.9845
2.3376
2.5
5
79.469
2.1752
According to the analysis in section 2.3.3, the stability o f the Ty(2,4) (F D )2T D scheme
depends on the material properties. Hence, for the first two permittivities in Table 2.2 we
choose At = 4.25442ps and the last two at At = 4.8622p s. Note that a smaller At
needs to be used if a small e"r(co) is chosen.
2.8.2 Numerical Results and Discussions
The approximate tim e steps to reach the steady state needs to be identified first in the
lossy case calculation. When the average input power is equal to the average dissipated
power, the system reaches the steady state. The dissipated powers for the first case with
e'r(o)) = 2.5 and e"r(co) = 0 .1 are plotted in Figure 2.33 and the corresponding stored
energy is shown in Figure 2.32. Note that the input power is calculated from (2.93) since
77
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the induced E M F method used is not accurate to calculate the input power at the probe
location. The stored energy is calculated from the integral o f the difference o f input power
and dissipated power. In order to determine the average power, a sine function is used to fit
the calculated data between 33,900 and 34,000 time steps in those two figures as shown in
Figure 2.34 and Figure 2.35. When the difference o f this fitting data and the dissipated
power data in the period o f zero to 35,000 is plotted, we obtain Figure 2.36. From Figure
2.36, it is observed that the approximate time step to reach the steady state is about 25,000
for this case. By using the same approach, the time steps for other cases are obtained. The
average loss power and the time-average stored energy are listed in Table 2.3.
Table 2.3 The loss power and stored energy
e"r(co)
The average loss
power
The time-average
stored energy
0.1
2.45e-10
8.905e-18
0.5
1.1258e-9
8.8e-18
2.5
3.035e-9
7.55e-18
5
3.205e-9
7.45e-18
In Table 2.4, the major resonant frequency and the Q factor in the fifth column are
calculated by using Prony’s method. In order to obtain the Q factor by Prony’s method, a
windowed sinusoidal source is used. The fourth column is the Q factor calculated from the
definition. The number o f time steps to reach the steady state is roughly equal to the Q
78
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table 2.4 Properties of four different lossy cases
8 " r (CD)
M ajor
Resonant
Frequency
(G H z)
Frequency
Shift (% )
0.1
2.4503
0.12
0.5
2.45
2.5
5
Q
Prony’s
Method
Approximated
Tim e Steps
559.5
523.6
25,000
0.13
120.3
125
5,000
2.45
0.13
38.3
37.1
1,400
2.45
0.13
35.7
34.2
1,300
_ _
Stored energy
0 Power loss
factor divided by twice the resonant frequency times the period o f one time step. Hence, a
very long time integration is needed i f a high Q cavity is dealt w ith. The field distribution
calculation scheme for the high Q cavity is shown in Figure 2.31. This algorithm is easy to
corporate with other temperature related equations to calculate the temperature change in
a cavity.
The field distribution o f Ey in the steady state and the fitting sine function are plotted in
Figure 2.37. The ratios o f the calculated and fitted Ey are plotted in Figure 2.38. The ratios
are no longer close to 0.4 as that in Figure 2.28; however, the ratios at those points which
are near to the m iddle point o f the material sample is still around 0.4. In this case, the
material sample is much larger than that in section 2.7, so the field distribution is no longer
easy to be fitted by a pure sinusoidal function which is shown in Figure 2.37.
79
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Run F D T D codes for a chosen time
steps with B H windowed sinusoidal
source.
A
Use Prony’s method to determine
the major resonant frequency and
the Q factor for this cavity.
Estimate the time steps to reach >\
steady state solution.__________
'
Run FD TD codes for 1/4 to 1/3
time steps of estimated steady state
time steps with a continuous
sinusoidal source.
Estimate the field distribution at
steady state by Prony’s method or
other estimation techniques.
Figure 2.31 The flow chart of calculating the field distribution at
steady state fo r high Q cavities.
80
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 2.32 Stored energy for the cavity with m aterial sample of
e'r(co) = 2.5 and e"r((o) = 0.1 . Note that one tim e step is equal to
4.25442 p s and this is a downsampling plot.
81
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
x10‘
•9
1.4
Dissipated
power
1.2
0.8
0.6
0.4
0.2
0.5
2.5
1.5
T im e steps
Figure 2.33 Instantaneous dissipated power fo r the cavity with
m aterial sample o f e'r((o) = 2.5 and e"r(co) = 0 .1 . Note that one
tim e step is equal to 4.25442 ps and this is a downsampling plot.
82
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.5
x 10
1.8
1.6
1.4
1.2
8.905e-18
2 0.8
to
0.6
0.4
0.2
3.39
3.391
3.392
3.393
3.394
3.395
3.396
Time Steps
3.397
3.398
3.399
Figure 2.34 The plot o f the fitting and calculated data for stored
energy for the first case. The line with circle is the fitting data and that
w ith cross is the calculated power data. The average input power is
2.45xlO'10 and the tim e average stored energy is 8.905e~18.
83
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.4
4.5
3.5
a>
?
45e-10
o
Q_ 2.5
tn
<n
3
1.5
0.5
3.39
3.391
3.392
3.393
3.394
3.395
3.396
Time Steps
3.397
3.398
3.399
3.4
x1Q<
Figure 2.35 The plot o f the fitting and calculated data for dissipated
power for the first case. The line with circle is the fitting data and that
with cross is the calculated power data. The average input power is
2.45xlOT10.
84
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 2.36 The plot of the difference between fitting and dissipated
power data. The one time step is equal to 4.25442 ps.
85
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-
0.01
-
0 .02
5 -0.03
Calcualted Ey
O O O O O O Fitted Ey
2=0.04652
y=0.014
-0.05
-0.06
-0.07
0.01
0.02
0.03
0.04
x(m)
0.05
0.06
0.07
0.08
Figure 2.37 Field distributions of Ey along x axis at 0.14868 |is. e'r(co) is
equal to 2.5 and s"r(co) is 0.1.
86
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.01
0.02
0.04
0.03
0.05
0.06
x(m)
Figure 2.38 The ratio o f calculated and fitting Ey along x axis at
0.14868 pis.
87
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.07
C H A PT E R 3
SOLVING M A X W ELL’S EQUATIONS BY BODY
OF R EV O LU TIO N FD TD
In the study o f the interaction o f microwave radiation w ith non-ionic materials, a mate­
rial sample is placed in a cylindrical E M cavity which is excited with a fundamental cavity
mode. M ost of the time, objects that are symmetric about an axis are encountered; hence
the body o f revolution (B O R ) F D T D is examined and used to solve the problems involv­
ing cylindrical cavity loaded w ith symmetric material sample.
The cylindrical cavities w ith perfectly electrically conducting (PEC) boundaries are
first studied in this chapter. The B O R FD TD formulation o f M axw ell’s equations in cylin­
drical coordinates is first constructed in section 3.1. Two sets o f equations are derived for
the F D T D formulation. The selection o f which set of equations to be used for field calcu­
lation in cavity problem is determined by the characteristic o f cavity modes. Then the sin­
gularity problems at p = 0 is discussed. The Blackman-Harris (B H ) function is used to
construct the excitation source in the cylindrical PEC cavity calculation since the sidelobe
o f B H function is approximately -92 dB. The excitation probe is assumed to be located at
a point on the cylindrical w all. In section 3.2, the surface impedance boundary condition
(SIB C ) is used to simulate the finitely electrically conducting (FEC) boundary. Three
88
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
methods o f the F D T D form ulation for SIBC is presented and the frequency domain
approximation is used in the actual calculation since it giwes stable numerical solutions.
For the cavity with lossy w all, a continuous source is used and the time step to reach the
steady state is studied. T M o 12 and T E ^ modes o f the ermpty cylindrical PEC and FEC
cavities are calculated and the results o f the former case is com pared to the analytical solu­
tion in section 3.3. The field distributions of T M 012 and T E t l l modes
with a cylindrical
PEC cavity loaded by a small cylindrical material centered in p = 0 are then calculated.
The results o f the form er case are compared with corresponding theoretical estimation.
The field distributions o f the cavities with a thin rod or a th in disk material are also calcu­
lated and compared to theoretical estimations.
A fter the induced electric field inside the material samplle is accurately quantified, the
dissipated power density inside the material sample can b e calculated. This dissipated
power density acts as the heating source to raise the tem perature o f the material sample.
3.1 The BOR Formulation of Maxwell’s E quations[28]
Consider M axw ell’s curl equations in non-dispersive nnedium written in cylindrical
coordinates in a linear material,
V x H = [ e ] | | + [(7]E + / r
(3.1)
V x E = _ [ p ] M + [< * * ]#
(3.2)
Ot
where
[e ],
[c ],
[p .], and
[a *]
are electric perm ittivity, electric conductivity,
permeability, and magnetic conductivity in tensor form , respectively. J s is the known
89
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
impressed current source. These constitutive parameters can be expressed in tensor form
in cylindrical coordinates as
« pp a p<l> a pz
[a ]
=
_a-P
where
a
represents
the
relative
(3.3)
a <i>p a <D«D
a zz
perm ittivity,
relative
permeability,
the
electric
conductivity, or magnetic conductivity. It is noted that if the medium is dispersive the
approach used in chapter 2 can be applied.
The azimuthal dependence o f field is expressed as a Fourier series,
oo
E =
^
( e Hcosm<{> + e„sinm <l>)
(3.4)
2 j ( h uco sm § + h vs i n m § )
(3.5)
m = 0
OO
H =
m = 0
OO
Ss =
X
(jucosmfy + j vs in m § )
(3.6)
m= 0
where m is the mode number and is an integer because o f single valueness o f azimuthal
dependence o f the field. Fourier coefficients e u , e v , h u , h v , j u , and j v are dependent on
r,
z, t,
and m where u
and v stand for the even and oddazimuthal
dependence,
respectively. The fields, E and H , in (3.4) and (3 .5 ) describe the fields atany point in the
entire space o f interest because the objects considered are symmetric about the z-axis. I f
the objects considered are not symmetric about the z-axis, then different pairs o f model
90
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
expansion are needed for different regions, followed by matching boundary conditions
between those different regions. It is not suitable to use BO R scheme for this nonsymmetrical problem. Note that if the loaded material is lossy, then M axw ell’s equations
with B and D pair are needed to perform the Fourier series expansion. Other relations for
the lossy case are described by the Debye equation as we did in chapter 2.
Substituting (3.4), (3.5), and (3.6) into (3.1) and (3.2), we have the following pair of
equations:
(3 .7 )
(3 .8 )
The above two vector equations can be separated into two independent groups o f six scalar
equations. These groups represent modes that are azimuthally perpendicular to each other.
The first group is as follow :
dz
p z' u
JP'V
(3 -9 )
(3 -1 0 )
(3 -1 1 )
^ e (j) , u
dz
m
p &z' v
+ [o * ]ft* v = -
91
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3 -1 2 )
(3 -1 3 )
The second group can be obtained by a similar procedure. Note that the signs o f terms that
contain ^ in the second group are different from those in the first group. The summary of
fields that depend on r, z, and t for those two groups are listed in Table 3.1.
Table 3.1 BO R representation o f M axw ell’s equations
Fields o f first group
e p,
Fields o f second group
V
e p, u
e Q, u
V
e z,v
^P, u
V v
^0 . V
\ u
^Z, u
K v
Assume that those constitutive parameters have the biaxial tensor form in cylindrical
coordinates given by
aP
[a ] =
0
0
0
a0
0
0
0
az
92
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Then these two groups o f equations can be represented in m atrix form [27] as
0
a.
-d , ± z
P
o
-a .
_m i a
+
a -p
P pd p
a
denotes
da
_
o
"p
0
m i a
p pap'
where
(3-16)
e <t>
e.
(P o M r +
a *z)hz
P
0
a.
(H o lip d ,+ G *p)h p
(e0epaf + o p)e p + j p
h9
(e 0e<j,af +
■3p h«
h Z_
0
+ j^
(3.17)
(e 0ezdt + a z)e z + j z
and a = p , z , t .
3.1.1 Mode Selection in BOR Algorithm
There are two important issues that need to be determined in applying B O R algorithm.
The first one is which group o f equations should be used and the second one is how many
modes should be included to solve the problem. These two issues are determined by the
incident wave or the impressed current source.
For scattering problem, the field distribution of the incident wave determines the
number o f modes and the group o f equations to be used. Consider the incident filed[25],
E = /( p )E ^ r —^j[pcos<j> —<{>sin<l>] = p E lp + § E l$
(3.18)
where c is the speed o f light and p , $ , and z are the cylindrical coordinates. The E lp is
even functions with respect to (j) so is H \ by M axw ell’s equations. Sim ilarly, E \ , H lp ,
93
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
and H lz are odd functions w.r.t <j>. The total fields denoted by a superscript t w ill preserve
this angular dependence[19] because o f the symmetry about the z-axis o f the object.
Therefore, these total fields components can be expanded in terms o f a Fourier cosine or
sine series. The £ rp , H \ , and E *z , are even functions and H lp , E \ , and H lz are odd
functions w.r.t (j); hence the second group o f equations is selected to solve this problem.
From (3.18), only cos<{> and sin<{> are involved, so only mode m = l is needed to solve this
problem.
In the cavity problem we considered, the source is assumed to be a line current source
located at a point (<J>0, z0) . The source can be expressed as
J
8(<j)-<t>0)
= - p / ( P)8 (z -Z o )g (t)
(3.19)
r
where / ( p ) and g(t) are variation o f J along p and t, respectively. Symbol 8 is Dirac
delta function. I f we let <j)0 to be zero then the delta function 8(<j>) can be expanded by
Fourier cosine series as:
oo
5(<J>) =
e"Om
—-cosm<{>
71
V
771 =
(3.20)
0
where
[1
=
2
m = 0
m *0 '
(321)
Form theory, the lowest TE mode is TE q h and its p component, E p , is equal to zero.
Hence, the <(> dependence o f E^ is cosm(J) in order to keep E^ from becoming zero for
94
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TTTou mode [26]. Therefore, the first group o f equations need to be used in this FD TD
calculation for cavity problems and (3.19) becomes
J = P X
m
=
jm cosm§
(3.22)
0
where
im =
- * o )s (0 •
|J
JL
(3.23)
3.1.2 The BOR-FDTD Formulation[28]
In this section, the B O R -FD TD formulation w ill be derived by using the leapfrog
scheme with cr* = 0 . Consider the first equation in (3.17),
-
^ eP
ep"07
m u
where the right hand size w ill be evaluated at t =
t =
+
(324)
(5f>ep ~ J p ~ ~ ^ ^ h z
Ar since h field is evaluated at
The standard F D T D notation w ill be used in this derivation, p = iAp
and z = y’A z . Using field locations in Figure 3.1, (3.24) becomes
ep r n + \ . .
..
n ..
.N,
( I , J ) - « p(l, J ) ]
_
un+W2(j
*
n+l/2,.
=
-C Tp e p
A _ /,n + l/2 /;
{ ,J)
/ _
..
. n + l / 2 , , . ..
( I, j )
(i, J ) - yp
1\
U . 7 - 1 ) ^ _ j n _ hn+
^
'
l/
( 3 ’2 5 )
2( /
P; + 1/2 2
The components, e” + 1/2 (i, j ) and y” + 1/2(z, y ) , have to be approximated by those values
adjacent to them since they are sampled at integer time step n . The second-order
95
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
approximation used here is
jJi + 1 /2 /.. .v
Jo
(*>./) =
- I,
3 / r ‘tt J) +
S ) - f T ‘O'. J)
(3.26)
where / is e o r / . The F D T D representation o f equation (3.24) is then
n + 1/ .
eD
..
n ,.
..
„
n ~ l ,.
(«, y) = A D • eQ(i, 7 ) + 5 p • ep
..
(i, j )
(3.27)
-C „
J„
TT--------------- ±------- K » l/2(<, ])
+
l* z
P t + 1 /2
where
3 °x
1 - t — ^ -A r
4 e 0£x
A. =
(3.28)
3 CTX
1+
- Ar
8e0sx
1
B. =
Ar
8 e 0 £x
(3.29)
3 CTx
1+ r ^ A r
8s0Sx
-A t
c
£0£x
=
(3.30)
3 CTx
1 + x — —Ar
8 e o£x
and x = p, <|), o r z. Note that f l + 1 / 2 is not replaced by the approximation o f (3.26) in
(3.27).
The other five F D T D equations can be derived by using the sim ilar procedure and are
listed below:
96
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
e "^
V , j)
= A^ ■e f t , j ) + B t>- e l
h ^ l / H i , j ) - h ^ l/2( i - U j )
\ n + 2 ,.
-a
7,j>
\i, j)
(*, 7) + ---------------
(3-31)
Ap
h s + i/2V ' i ) - h r l / 2 « - j - » Az
(*,/) =
A -en
Ji,j) + B - e"
\i,j)
-C .-
P i+
(3.32)
1/2^$+
7*) —P f- l/2^ > + l/2 ( l — 1» J )
p^Ap
•*$(*. 7 +
:> ,;)T G p5 ^ ( i , y ) + G p[-
u n + 1 /2
P
]
Az
]
(3.33)
.
./-i
r« s o ',y + i)-c s (*'.y )
+ 1> y ) —e? a y ) i
» { * 1/2( i.7 ) = * r 1/20 '.J )-G » [ p
^
-------------- A o "
J <3'34)
A " l/2 ( i,7 ) = * z" - l/2 (i, » ± G . —
* "( i , j ) - G
P / + 1 /2
where Gx = A r /( p 0M-x) ’ x = p,
(3.35)
P /+ l/2 ^ P
o r z and p,- = ( i - l / 2 ) A p
and p 1/2 = p0 = 0 ,
and the fields associated with the coordinate (i, j ) are shown in Figure 3.1. In Chen’s paper
[28], those B O R F D T D equations are obtained by using a first-order approximation,
rJl +
Jo
1/ 2 / .
.v
fn
p+ \ i , j ) + f n
p(ij)
(*>./) =
(3.36)
97
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
h
e
h
e
time
■ X -
n-1/2
n
n + l/2
n+1
n+3/2
I (2 ,y )
L.
J
-2
l_
(iJ+l)Sto
l
gp
L. __
_
V + l , j + 1)
(i+ tj)
F IG U R E 3.1
space.
The field locations for B O R FD TD in tim e and
98
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-I
for fields at the half time step. Smaller time steps or spacial steps are required in Chen’s
paper than that in this chapter. This problem becomes more serious in calculating the
cavity modes with n # 0 .
3.1.3 Singularity in BOR-FDTD Formulation at p = 0
As observed from (3.27), there is a p 1/2 factor, which is equal to zero, in the
denominator o f the fourth term when calculating the ep(0 , j ) . Hence, ep(0 , j ) is infinite
at p = 0 and this makes h^(0, j ) in (3.34) infinity, also. The filed, hz(0, j ) , is infinite by
the same reason in (3.35). Other fields, ep ,
, and hz at the p = 0 also exhibit
singularities; however, the actual fields must be finite in both the time and frequency
domains. Hence, the these singularities must be removed before above F D T D equations
can be used for time stepping.
As observed from (3.31) and (3.32), only the components /^(O , j ) and hz( 0, j ) are
needed to update the adjacent ez( l , j ) and e^(l, j ) fields internal to the mesh. The
ep(0, j ) and ep(0 , j + 1) are needed to evaluate h ^ ( 0 , j ) . The field, h^(0, j ) , is not
needed actually to calculate the e (1, j ) since
en
z + \ l , j ) = Az -e n
z ( l , j ) + Bz - e n
z \ l , j)
n+ l
-c z
(3 .3 7 )
P3/2*t +1/2( 1. J) - P l / 2
p ,A p
hS *
1/2(0. J)
with p I/2 = 0 ; hence, ep(0 , j ) is not needed, neither. In actual calculation, ^ ( 0 , j ) and
99
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ep(0 , j ) are set to zeros to avoid cumulation error.
Now, only h.{ 0, j ) is needed to update the fields point in the F D T D lattice. Note that
hz is zero at p = 0 fo r m ^ O by (3.17); hence, we only need to evaluate hz(Q, j ) for
m — 0 . From the integral form o f M axw ell’s equation in the time domain, we can obtain
the following time update equation for hz(0, j ) when m = 0
* r I/ 2 ( ° .y ) =
■
(3'38)
Once hz{0, j ) is known the rest o f the field components can be evaluated using (3.27) and
(3.31) to (3.34).
3.2
Surface Impedance Boundary Condition
When the boundary o f a cavity has a finite electrical conductivity (FE C ), there are two
F D T D approaches to calculate the field distributions in the E M cavity. For a cavity with
good conductor w all which is considered in this chapter, the skin depth and the local
wavelength are very sm all compared to the radius o f curvature o f the cavity wall. Hence,
the planar surface impedance boundary condition (PSIB C ) is used to approximate the
lossy conductor w all. For more accurate simulation, an absorption boundary condition
(A B C ) has to be used outside the cavity and as shown in Figure 3.2. In this chapter, only
the PSIBC case is considered.
3.2.1 Planar Surface Impedance Boundary Condition
The surface impedance is inherently a ffequency-domain concept. Consider a timeharmonic plane wave incident at angle 0 {- on the lossy dielectric half-space as shown in
100
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 3.3.
The surface impedance for a planar interface at z = 0 can be defined through the
following relation:
Et(as) = Z (c o )[n x H / (co)]
(3.39)
where Z(co) is the surface impedance. The surface impedances o f interface for T E and
T M plane w ave[30], respectively, are
(3-40)
Z r e (“ ) - ^
Z rA /(w ) — Z ocos02
(3-41)
where
z =f
k
j(DlXk )
1/2
(3.42)
{j(QEk + a kJ
1/2
cos09 =
1-
2n
sin
0.
*i
. 2n .1/2
= (1 — sin 02)
• 2,
When sin 02 « 1 , then the surface impedance can be rewritten as
101
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3-43)
^
EEC W all
Inner Em pty Cavity
ABC
ZZZZZ
zzz
zzzz
Figure 3.2 FD TD configuration for cavities with FEC wall
102
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Free Space
Lossy Dielectric
Figure 3.3 Coordinates for the incident and reflected plane
waves upon a lossy surface
r
Z (G > ) = Z TE (CO) = Z TM( CO) =
7 cop9
\ 1/2
+ CTJ
(3.44)
where q 0 = J \ i 0/ e0 and s = j ( t i . Equation (3.44) is the planar surface impedance used
through this chapter.
3.2.2 FDTD Implementation of Planar Surface Impedance
In order to do F D T D calculations, the tim e domain expression o f (3.39) needs to be
obtained. The convolution in time domain is involved since there is a m ultiplication in
frequency domain. The convolution in time domain needs to be approximated by some
103
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
recursive expression in order to run on a computer. Three approximations are discussed in
this section. The tim e domain approximation[29] is first considered; however, this
approximation encounters some divergence problems when evaluating an integral and is
not used to the F D T D implementation. The frequency-domain[31] and Z-transform
approximation [3 2] [33] are easier and more efficient to calculate than the time domain
approximation.
3.2.2.1 Time-Domain Approximation by Prony’s Method
The tim e-dom ain surface impedance is
su
z{t )
\ 1/2
= M ^ei r JJ
(3.45)
+
where U ( t ) is the unit-step function, a = -c r2/( 2 e 2) and I 0 and
are the modified
Bessel functions o f the first kind o f zero and one order, respectively. Hence, the time
domain expression o f (3 .3 9 ) is
E t(t) = -n0
[ n x f f , ( l ) ] + j aeat[ I v{at) + / 0(<zr)][n x H t( t - x ) ] dx I
(3.46)
and the discretization version of (3.46) is
(3.47)
where
F 0(<?) = f 9 1}[ 1 - |Y - < 7 l] ( a A r ) e ( aA/)r[ / 1( a A / Y ) + / 0( a A t y ) ] ^
J<t>
and
104
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3-48)
* =
0
q = 0
q—1
q> 0
(3.49)
The Prony’s method is used to approximate the F 0(q) and (3.47) is then
Et\
= r\2
(3 .5 0 )
k= 1
where Gfc(n) = CkH t | + p.fcG * (/z - 1) and k = 1, 2, 3, . . . g .
3.2.2.2 Frequency Domain Approximation
The surface impedance, (3.44), is rewritten as
zw
where s = y'cD,r|2 =
t) 0
(3 -51)
-
hi2r
/— , and a = a 2 / e 2 . Equation (3.51) can be approximated by
H e2r
a rational polynomials in frequency domain and is rewritten as
(3 -5 2 )
i
=1
where s' = 5 / a , the L is number of terms needed in that approximation and cot- are
known poles. The L and cot- are determined using a rational Chebyshev approximation and
is listed in Table 3.2[31]. This approximation is over the real axis interval s' = [0, 3 ]
which w ill accommodate most material up to several tens o f Gigahertz. The (3.39) in splane is then
105
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
L
E t( s ) ~ r \ 2
~
, .
aCa (t+ s
i = 1
[n x
.
(3.53)
'
Assuming the waves are piecewise linear in time, (3.39) in discrete-time domain is
given by
n
E tI
n
= y\2 [ n X
L
f |"J -
n
] £ A i\
/= i
(3.54)
where
and
C :
Pn =
-flfrt.A;
[1 + (e
'CO.
- l)/(Amcoz-)]
(3.56)
P12 = Ti2^[l/(AtaM t.)-e"aa3,Ar(l + l/(Afacof))]
(3.57)
—aa>:At
p i3 = e
.
(3.58)
3.2.2.3 Z Transform and Digital Filters Approach
The term, C [/ ( ( o i + s ') , in (3.52) is a analog lowpass filter with a single pole. The
digital version o f this term is a digital H R filter. Hence, this become an HR digital filter
design from the analog filter in digital signal processing terminology.
A popular IIR filter design is the bilinear transformation which is appropriate for
lowpass, bandpass, and highpass filters. The bilinear transformation is a mapping from the
106
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table 3.2 Rational Approximation Results
Number o f
Terms
8
ci
(0-i
2.60266906380e-10
3.31195847511e-8
1.58725612150e-6
4.35701245008e-5
8.08657638450e-4
1.07197254652e-2
9.48589314718e-2
0.3933263026213
2.43632801126e-7
1.06720343642e-5
1.59612026325e-4
1.59383673429e-3
1.20627876345e-2
7.22902130819e-2
0.3267172939002
0.8680755110248
s-plane to z-plane and can be linked to the trapezoidal formula for numerical
integration[34]. The substitution o f variable for bilinear transformation is
(3 -5 9 )
ru + z -0
where T is the period o f one time step or the sampling factor. Applying the bilinear
transformation, (3.59), to (3.53), then the time discrete version o f (3.39) becomes
\n
E t|
r
,n-i
= r i 2^ x / f , |
L
_x |n
- ^ f i \
(= 1
(3 .6 0 )
where
f 2 —a(0;T>
-7\-a.|rt-l
/ ' a C [X\2T \ r
,n
/«
+ x— - ^ ) n x H A + h x H l \
_ l2 + aco.frj
‘
U + fltO/rJ.
1
1
J
(3 -6 1 )
The bilinear transformation maps the entire left h alf o f the complex plane inside the
unit circle in z-plane and the imaginary axis in the complex plane becomes the unit circle
in z-plane. As a result, frequency compression or frequency warping w ill take place when
107
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
transferring from the analog system to the digital system. The bilinear transformation is
most useful when this distortion can be tolerated or compensated.
3.2.3 Fields Calculation on the Cavity Wall
From Figure 3.1, there are two tangential electric fields and one normal magnetic field
on the physical cavity walls. B y (3.50), (3.54), and (3.60) the tangential electric fields are
evaluated from the tangential magnetic fields on the cavity walls at the same tim e step.
However, the tangential magnetic field is evaluated at half-tim e step before the current
time step and half cell in front of the cavities wall according to Yee’s F D T D lattice in
Figure 3.1. In the practical F D T D calculation, these tangential magnetic fields are used in
(3.50), (3.54), and (3.60) to obtain the tangential electric fields on the cavity walls.
Equations (3.54) and (3 .6 0 ) can be rewritten as
(3.62)
where
s .n
A t-|
i-
.71—1 /2 -]
= p n \ h x H t \h
r
_ i. ,71 —3 /2 -]
j + p i2\ h x H t\h
,n—1
J + p t-3A ,
(3.63)
and
(3.64)
where
108
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Note that the subscript, h, stands for a half space step before the walls. The normal
magnetic fields on the cavity walls are obtained by using the normal FD TD
implementation o f M axw ell’s equations.
3.3 An Empty Cylindrical Cavity
The physical configuration o f the cylindrical cavity is shown in Figure 3.4. The radius
and height o f the cavity are a and h and that o f loaded material are b and I. The wall o f the
cylindrical cavity is assumed to be perfectly electrically conductive (PEC) or finitely
electrically conductive (FE C ) and the excitation probe is located on the side o f the
cylindrical cavity.
The TMqx 2 mode o f an empty cylindrical cavity w ith PEC boundary is first calculated
to verify the program. The g(t ) in (3.19) is a Blackman-Harris (B H ) windows function
with the central frequency o f 2.4571GHz and 0.1G H z bandwidth. The B H function makes
the source turn o ff automatically and the 0.1G H z bandwidth is chosen because only
TMqx2 is desired. The dimensions o f the cavity are a = 0.0762m and I = 0.15458m and
the dimensions o f grids are Ap = 0.0006m and A z = 0.00118 which are lees than
1 /2 0
where 1 is the wavelength. The numbers o f partition along p and z are 127 and
131, respectively, and the one time step is 1.69363 ps. The field distributions o f E p and
E z are plotted in
Figure 3.5 to Figure 3.6 and the total time step is 15382. Form these
figures, we conclude that the result of m = 0 mode is the dominant mode which is much
larger than that of m = 1 mode and higher modes. The field distribution of E on the p - z
plane is plotted in Figure 3.7. Hence, only m = 0 mode is considered in the calculation o f
109
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
T M q12 mode for the FEC and the loaded cavity. F o r T E U1 mode, the numbers o f partition
are 127 and 67 along p and z directions and the one time step is 1.29676 ps. The field
distributions o f E p ,
and E on the p —z plane f o r T E ^ mode are plotted in Figure 3.8
to Figure 3.10. For T E ^ mode, the m = 1 mode is dominant mode as can be seen from
these figures.
For the empty cylindrical cavity with FEC boundary, the FEC boundary is replaced by
PSIBC. The continuous sine function is used fo r g ( t ) in source equation (3.19) since
there is a loss on the boundary. As observed in Figure 3.11 with Ar = 1.2967608p s , the E
field achieves a steady state about 23,135 tim e steps for cr = 10
2
S / m and 120,000 time
steps for a = 104 S / m . It w ill need more than 600,000 time steps to evaluate the field
distribution of T E U i mode for cavity with <7 = 106 S / m which is usually encountered
in regular cavities.
The field distributions o f E z and E on the p -z plane are shown in Figure 3.12 and
Figure 3.13. There is some Ez near the boundary w all form Figure 3.12; however, the
value o f Ez is much smaller than that o f E p. The E field on the p-z plane is nearly identical
in Figure 3.13 compared to that in Figure 3.10. Therefore, it is reasonable to use PEC wall
assumption to do the F D T D calculation for field distributions in cavities since most cavity
wall is made o f metal with conductivity larger than a = 106 S / m . Moreover, this PEC
assumption requires much less time steps in F D T D simulation.
110
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 3.4 The physical configuration of a cylindrical cavity loaded
w ith a m aterial sample.
I ll
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.1
m= 0
z =0.0531
0.08
0.06
0.04
0.02
m= 1
-
0.02
0.01
m= 2
0.02
0.04
0.03
0.05
0.06
0.07
P (m)
P=!
0.1
m= 0
p =0.0372
0.05
m= 1
m= 2
-0 .05
-0.1
0
0.1
0.05
0.15
z (m)
Figure 3.5 The variation of E p o f T M 012 along the p and z directions in
a PEC empty cylindrical cavity.
112
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.05
m= 0
Z
m= 1
=0.0295
m=2
z =0.0531
-0 .05
m=0
o
0.01
0.02
0.04
0.03
0.05
0.06
0.07
P (m )
0.1
m= 0
p =0.0372
0.05
m= 1
m=2
p =0.057
-0 .0 5
-0.1
0
0.1
0.05
0.15
z (m)
Figure 3.6 The variation of Ez of TM0i 2 mode along the p and z
directions in a PEC empty cylindrical cavity.
113
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.15
0.1
z (m )
0.05
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m )
Figure 3.7 Plot o f E field of T M q ^ mode on the p-z plane in an
empty PEC cavity.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
15
m= 1
= 0 .0 3 3 0 2 4
z = 0 .0 4 6 4 4
0.5
0
0.01
0.02
0.04
0.03
0.07
0.06
0.05
P (m)
1.4
m= 1
p = 0 .0 3 7 2
p = 0 .0 5 7
0.8
,0.6
0.4
0.2
-m -=
-0.2
0.01
0.02
0.03
0.04
0.05
0.06
z (m)
Figure 3.8 The variation of E p of T E m mode along the p and z
directions in a P EC empty cylindrical cavity.
115
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.S
m= 1
z = 0 .0 3 3 0 2 4
z = 0 .0 4 6 4 4
0.5
0
0.01
0.02
0.04
0.03
0.05
0.06
0.07
P (m)
1.4
m= 1
P=!
0.8
,0.6
0.4
0.2
m-=-0—
-0.2
0.01
0.02
0.03
0.04
0.05
0.06
z(m )
Figure 3.9 The variation of E^ of T E 1U mode along the p and z
directions in a P EC empty cylindrical cavity.
116
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.06
0.05
0.04
z(m)
•3 *-
0.03
0.02
0.01
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (rn)
Figure 3.10 Plot o f E field of T E m mode on the p-z plane in an
empty PEC cavity.
117
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.4 A Loaded Cylindrical Cavity
For a loaded cavity with PEC boundary, three cases for T M 012 mode and one case for
T E U1 mode are calculated. The first one is a cavity loaded w ith a small cylindrical
m aterial, with I ~ 2 b , centered at p = 0 . The second one is a cavity loaded with a thin
rod material w ith I » 2b and a cavity loaded with a thin disk material with I « 2b is the
last case. For T E in mode, only the case o f a cavity loaded w ith a small cylinder material
is studied.
3.4.1 Sm all cylindrical sam ple for TM012 mode
A m aterial sample, having the diameter equal to the height, is placed in the center o f
the cylindrical cavity. The diameter is 0.006 meter and the height is 0.0059
meter; that is I ~ 2b. The induced electric field inside the material sample can be
estimated
by
the
electrostatic
field
induced
inside
of
a
dielectric
sphere
as
E z = ( 3 / ( 2 + sr)) E l‘ .
The numbers o f partition used along the p and z directions in this FD TD calculation
are 127 and 131, respectively, and the period of one time step is 1.29676 ps. Observing the
electric field distribution in Figure 3.7, only Ez component o f the electric is significant
near the center o f the empty cavity. Due to the small dimensions o f the material sample, E z
is still the dominant component near the center of the loaded cavity if other components
are compared w ith Ez in Figure 3.14 and Figure 3.15.
The ratio o f E z component o f the induced field inside the m aterial sample to that o f the
induced field in an empty cavity is plotted in Figure 3.16. The ratio is between 0.6 to 0.75
118
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
and it depends on the position o f p and the fitting sine function. The electrostatic
estimation o f that ratio is 0.667 with er = 2.5. The numerical results and the theoretical
estimation are in satisfactory agreement.
3.4.2 Thin rod case for TM0i 2 mode
A material sample with the shape o f a thin rod, having its height much larger than its
diameter, is placed in the center o f the cylindrical cavity. The dimensions o f the material
sample are I = 0.13098 meter and b = 0.0018 meter. The ratio o f the diameter to the
wavelength is about 0.03. The numbers of the partition along the p and z directions and
one period o f time step are the same as those considered in section 3.4.1. The dominant
component near the center o f the cavity is still Ez if
Ep and Ez components are compared in Figure 3.17 and Figure 3.18. The field distribution
o f Ez along the p axis seems to be not affected at all in Figure 3.18. In fact,
is not
affected by the placing o f sample or the effect is too small to be noticeable. From Figure
3.19, we observe that the ratio of the z-component o f the induced electric field to that of
the electric field in an empty cavity is very close to 1 or slightly less than 1. The
electrostatic estimation o f that ratio is 1 based on the continuity of the tangential
component o f the electric field on the sample surface. This electrostatic estimation agrees
w ith our numerical results.
3.4.3 Thin disk case for TM012 mode
A material sample w ith the shape o f a thin disk, having its height much smaller than its
diameter, is placed in the center of the cylindrical cavity. The dimensions of the material
119
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
_20l---------- ■______ ■______ I______ >_______I______ ■______ I______ I
0
0.2
0 .4
0.6
0.8
1
1.2
1.4
T im e (sec)
1.6
x10' 7
Figure 3.11 Plot o f E p o f T E ^ mode versus tim e in an empty F E C
cavity with cr= 102 S / m and ct= 102 S / m w ith m = 1 .
120
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
sample are I = 0.00118 meter and b = 0.006
meter and the ratio o f diameter to
wavelength is 0.1. The z component o f the induced electric field inside the sample is still
the dominant component. The induced electric field o f this thin disk geometry can be
estimated theoretically by the boundary condition o f E z = ( 1 / s r)E z .
The numbers o f partition along the p and z axes are 131 and 262 and the one time
step is 1.05171 ps. The numerical results are shown in Figure 3.22. The ratio o f E z/ E z is
about 0.42 which is consistent with the theoretic estimation.
3.4.4 Small cylindrical sample for TE111 mode
The dimensions o f a small cylindrical material are I = 0.00516 m and 2b = 0.0048
m. The numbers o f partition along the p and z and one period time step are the same as
those used in calculating the TE111 mode in an empty cavity. The numerical results of
field distributions are shown form Figure 3.23 to Figure 3.25. There are some E z near the
center o f the cylindrical cavity due to the placing o f material sample as observed in Figure
3.25; however, the magnitude o f Ez is much smaller than E p or E^. When <j> = n / 2 , E^ is
zero since the first group o f equations in Table 3.1 is used. Hence, the E^ is the dominant
filed distribution near the center o f the cavity. From Figure 3.23, we observe that the ratio
o f induced field E p to E p in an empty cavity is about 0.7 which is close to the theoretic
estimation. Same results can be observed from Figure 3.24 for E^ at <() = 0 .
121
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
m=
0.8
0.6
0.4
0.2
z =0.031992
-0.2
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (rn)
Figure 3.12 Plot of E z of T E U1 along p direction in an empty F E C cavity
with conductivity 104 (S/m).
122
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.06
0.05
0.04
z (m )
0.03
0.02
0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m)
Figure 3.13 Plot o f E field of T E 111 mode on the p-z plane in an empty
FEC cavity with conductivity 104 (S/m ).
123
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0 .0 4
m= 0
0.02
z = 0 .0 7 4 3 4
2b=10Ap
-
l= 5 A z
0.02
-0 .0 4
-0 .0 6
-0 .0 8
z = 0 .0 7 9 0 6
-0.1
0.01
0.03
0.02
0.04
0.05
0.06
0.07
P (m )
m= 0
0.15
0.1
0.05
p = 0 .0 0 1 2
p = 0 .0 0 2 4
-0 .0 5
0
0.1
0.05
0.15
z (m)
Figure 3.14 Plot of Ep of T M 0i 2 mode in a PEC cavity loaded w ith a
small cylindrical m aterial sample.
124
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.9
m = 0
0.8
0.7
z = 0 .0 7 4 3 4
0.6
0.5
z = 0 .0 7 9 0 6
0.4
0.3
0.1
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m)
m = 0
0.8
p = 0 .0 0 1 2
0.6
0.4
p = 0 .0 0 2 4
-0 .4
-
0.6
0.05
0.1
0.15
z (m)
Figure 3.15 Plot o f E z o f T M q12 mode in a PEC cavity loaded w ith a
small cylindrical m aterial sample.
125
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.4
1.3
Ratio of E /E 1
1.2
0.9
0.8
p = 0.0024
0.7
p = 0.0012
0.6
0.072
0.074
0.076
0.078
0.08
z(m )
Figure 3.16 Plot o f the ratio o f Ez/ E z o f TMqx 2 mode along the z
direction in the PEC cavity loaded with a small cylindrical m aterial
sample.
126
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.082
0 .0 4
m= 0
0.02
z =0 .0 7 9 0 6
-
0.02
-0 .0 4
-0 .0 6
-0 .0 8
z = 0 .0 7 4 3 4
-0.1
0.01
0.02
0.04
0.03
0.05
0.06
0.07
P (m )
0.15
m= 0
0.1
0.05
p = 0 .0 0 1 2
-0 .0 5
p = 0 .0 0 2 4
-0.1
0
0.1
0.05
0.15
z (m)
Figure 3.17 Plot of E p of TM o12 mode in a PEC cavity loaded w ith
a thin rod m aterial sample.
127
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.8
m = 0
0.7
= 0 .0 7 9 0 6
0.6
z = 0 .0 7 4 3 4
0.5
0.3
0.1
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m)
m= 0
0.8
p =0.0012
0.6
0.4
-0.2
-0 .4
-0.6
-0.8
0.1
0.05
0.15
z(m)
Figure 3.18 Plot of Ez o f T M 012 mode in a PEC cavity loaded with
a thin rod m aterial sample.
128
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1 .0 2
p =0.0012
1-
p =0.0024
0.98
Ratio of E /E
N
N
0.96
0.94 -
0.92 -
0 .9 ----------------------- 1-----------------------L
0.07
0.072
0.074
0.076
0.078
0.08
z(m )
Figure 3.19 The ratio of Ez/ E \ o f TMo12 mode in a PEC cavity
loaded with a thin rod material sam ple.
129
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.082
0 .0 4 5
m= 0
0.04
0.035
2b=20Ap
0.03
0.025
0.02
0.015
z =0.0767
0.01
0.005
z =0.07729
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m )
0.03
m= 0
p =0.0048
0.02
0.01
p =0.0036
-0.01
-0.02
-0 .0 3
0.1
0.05
0.15
z (m )
Figure 3.20 Plot of Ep of T M 0i 2 mode in a PEC cavity loaded with
a thin disk m aterial sample.
130
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
m= 0
-0 .0 5
-0.1
z =0.07729
-0 .2 5
z =0.0767
o
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m )
m= 0
0.2
0.15
0.1
p =0.0036
0.05
-0 .0 5
-0.1
-0 .1 5
-0.2
-0 .2 5
0
0.1
0.05
0.15
z(m)
Figure 3.21 Plot of E z of T M 0i 2 mode in a PEC cavity loaded with
a thin disk m aterial sample.
131
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.3
=0.0036
M
Ratio of E /E
—
=0.0048
0.8
0.7
0.6
0.5
0.4'—
0.072
0.074
0.076
0.078
0.08
0.082
0.084
z(m )
Figure 3.22 The ratio o f Ez/ E zl o f T M #12 mode in a PEC cavity
loaded w ith a thin disk m aterial sample.
132
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.086
1.8
m= 1
1.6
1.4
1.2
z = 0 .0 3 5 0 8 8
0.8
0.6
0.4
0.01
0.02
0.03
0.04
0.05
0.06
P (m)
1.4
_ m = 1
0.8
p =0.0012
E
0.6
0.4
0.01
0.02
0.03
0.04
0.05
0.06
z(m )
Figure 3.23 Plot of Ep of T E m mode in a PEC cavity loaded with a
small cylindrical m aterial sample.
133
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1 .4
m= 1
1 .2
0.8
z = 0 .0 3 5 0 8 8
0.6
0.4
0.01
0.02
0.03
0.04
0.06
0.05
0.07
p (m )
1.4
m= 1
1.2
0.8
p =0.0012
0.6
0.4
0.01
0.02
0.03
0.04
0.05
0.06
z (m )
Figure 3.24 Plot o f
o f T E l n mode in a PEC cavity loaded with a
sm all cylindrical m aterial sample.
134
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0 .0 3 5
m= 1
0.03
0.025
0.02
0.015
0.01
0.005
z = 0 .0 3 5 0 8 8
0.01
0.02
0.03
0.04
0.05
0.06
0.07
P (m )
0.1
0.08
m= 1
0.06
0.04
0.02
-
p = 0 .0 0 1 2
0.02
-0 .0 4
-0 .06
-0 .08
-0.1
0.01
0.02
0.03
0.05
0.06
z (m)
Figure 3.25 Plot o f E z of T E m mode in a P E C cavity loaded
w ith a small cylindrical m aterial sample.
135
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
SOLVING M A X W ELL’S EQUATIONS BY FDTD
IN C Y LIN D R IC A L COORDINATES
In chapter 3, we considered the cylindrical geometries with rotational symmetries. The
purpose o f this chapter is to develop methods for treating problems which may not have
rotational symmetries. This chapter considers the second-order and Ty(2,4) FD TD
formulation and applications in cylindrical coordinates. The general 3-D FD TD formation
is considered and discussed in section 4.1. There are two problems in applying cylindrical
coordinate FD TD in 3-D . The first one is the singularities at p = 0 . The other
complication that arises in applying cylindrical coordinate F D T D in 3-D is that the cell
size in the
(j)
dimension decreases with decreasing p . This means that very small time
steps may be necessary in order to satisfy the Courant stability criterion, unless the region
in the vicinity o f p = 0 is filled with perfect conductor or otherwise excluded. The two
problems can be solved by using a set o f FD TD approximation equations for p = 0 after
a careful examination o f the 3-D F D T D M axw ell’s equations. This special FD TD
approximations at p = 0 requires the loaded material to be bi-axial magnetic and without
impressed or conductive current near the p = 0 . For second-order method, those FD TD
136
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
approximation equations are highly mode-dependent and not easy to be generalized as is
discussed in section 4.2. In section 4.2, source implementation for traditional second-order
F D T D approximation causes incorrect field distribution and this problem is solved by
introducting the Blackman-Harris type source.
Using the higher-order spatial finite difference scheme, the approximation order can
be controlled and general F D T D approximations at p = 0 can be obtained. Higher-order
F D T D also requires less number o f partition than that in second order scheme and hence
reduces the computational tim e. These higher-order method is studied in section 4.3. A t
the end of this chapter, the F D T D formulation o f constitutive equations for general Debye
and Lorentz material, an extension of Debye material in section 2.3, is derived in section
4.4.
4.1 Three Dimensional FDTD Representation of Maxwell’s Equations
in Cylindrical Coordinates
In this section, the 3-D F D T D expressions o f M axw ell’s equations are derived. The
differential M axw ell’s equations considered here are
V xE =
dt
^
(4 .1 )
V x H = ^ . + J C+ J S
dt
where D = [£ ]is , B = [ p ] H , Jc = [o jis , and J c is the source current. The above
constitutive parameters are further assumed to have the biaxial tensor form in the
cylindrical coordinate system given by
137
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ap °
[a ]
°
(4.2)
0 a* 0
=
0
0 a,
where a represents the electric permittivity, the magnetic permeability or the electric
conductivity. For nonmagnetic dielectric m aterial considered in this chapter, the electric
conductivity o f this material is assumed to be zero and the magnetic permeability is
assumed to be the same as that in air.
The scalar M axw ell’s equations in cylindrical coordinate are
dflp _ dE+ _ 1 dEz
—
dt
pd<f>
3B*
dt
dp
dz
z = l3 g p
pd<|>
dt
p3(jl
dz
dz
3p
=
1 d
(4 -3 )
(4 -4 )
1 3 r oE i
p3p
^
|J » ;
dD?
where E p , E^, E z , H r ,
dz
3»p
dt
Tt
dt
cp
- i
(4 -5 )
sp
'♦
l dHo
=
(4 -6 )
(4 -7 )
(4 -8 )
, and H . are electric and magnetic fields along p , <|), and z ,
respectively.
Using the F D T D notations, those finite difference equations for M axw ell’s equations
are obtained as:
138
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where SpCpff^) and SpCp-E^) need special treatment on p , and it will be discussed later.
The spatial locations o f D and H are plotted in Figure 4 .1 , E has the same location as
D does and so is B and H . The sequence o f FD TD calculations along time axis is shown
—X
-X
-X
-X
_x
in Figure 4.2. The D , E , and d B / d t are evaluated at the integer time step; but B , H , and
3 D / d t are evaluated at the half-integer time step. The order o f calculation at different
tim e steps are also described in Figure 4.2. Note that the boundary conditions are applied
139
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
after the E is calculated.
d
d
Note that the terms ^ -(p is ^ ) and
are not factorized out. The first term can
be rewritten as
fe p E J
=
J T ^ P g jj* .
(4 ,5 ,
The evaluation of E A in the right hand-side is E A n
which is not where E+ is
•P
<pl/+ l/2 ,y + l/2 ,(t
‘P
located and, thus, (4.15) is not used.
140
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A z
F IG U R E 4.1
F D T D lattice for cylindrical coordinates.
141
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
n-1
n
n -1 /2
I
1
/I —1
I
1
I
n~2
£>
n+1/2
time
-* D
B
a
(3)
1
1
a ~2
r7Z—1
(5)
H
(2 )
H
(4)
n —1 /2
(1 )
(6)
ae
a?
n — 1 /2
a?
Figure 4.2 The diagram o f order of F D T D calculations along tim e axis. The
meanings o f those steps are listed below.
(1) Using (4.9) to (4.11)
(2 ) Desired time stepping scheme for D .
A.
(3 ) Constitutive relation o f D and E which depends on material models.
(4 ) Using (4.12) to (4.14). A t this point, the boundary conditions are applied.
(5 ) Desired time stepping scheme for B .
A
(6 ) Constitutive relation o f B and H which depends on material models.
142
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4.2 The Second-order Cylindrical FDTD Scheme[35]
The second-order cylindrical F D T D scheme which applies the Yee algorithm to
cylindrical M axw ell equations is discussed in this section. The cylindrical FD TD
formulation o f M axw ell equations is presented in section 4.2.1. The singularity of
cylindrical FD TD equations at p = 0 and its traditional treatments are discussed in
section 4.2.2. The traditional source F D T D implementation is discussed in section 4.2.3
and an improved source implementation is also presented in this section. Finally, several
numeric results are presented in section 4.2.4. The problems caused by the traditional
source implementation and the improved implementation by utilizing Blackman-Harris
function are also discussed in this section.
4.2.1 The Second-order Cylindrical FDTD Equations
Applying the second-order central finite difference approximation to time and spatial
derivative o f fields, the (4.9) to (4.14) become the finite difference M axw ell’s equations in
cylindrical coordinates. The finite difference equations are listed below:
fH
.n + 1/2
zl<+ 1/2,/+ 1/2, k
E l” + *
— E \n
^
x
pl/ + i/2, j,k ~
p\ i + 1/2, j,k + e\i + l / 2 j . k X
ZL' + l/2 ,/-l/2 ,fc
Pi+ 1 /2 ^
(4-16)
YJ-
1/2
<t>l/ + 1/2, j,k+
JJ j/Z-F 1/2
1 /2
+ 1 /2 ,
j, k—1 /2
j
AZ
,n
Plt + 1/2, /, k
r TT ,n + 1 /2
+1
„ ,n
At
E v\i,j+
i>\
=
^6
+
_
i-------------l/2 ,k
p\i,j+W2,k e | . y+1/2j k X
,n
jj
.n + 1/2
zl/+ 1/2, j + 1/2, k
Az
jj
tfp f
pl/,/ + 1/2,k+ 1/2
Az
TT ,n + 1 /2
"
p 11, j + 1/2, k - 1/2
.n + 1/2
zIi- l/2 ,/+ 1/2, k _ j \n
* w, / + 1/ 2 , k
143
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4-17)
n+ 1
2A t
Z i , j , k + 1/2 = E ‘ lt.y.fc+l/2
e |/ y -i i f c + l / 2
/_
X
n+1/2
,n + 1 / 2
P'C l / 2 g ( i>1i + I / 2 , j , k + I / 2
P £ ~ 1/ 2 ^ (|) It - 1/ 2 , / , jfc+ 1/ 2
2
2
P i + 1 /2 — P t — 1 /2
d
- dP i —1/2 v
Pi+1/2
2
,/i+ 1/2
H pli,y
p\ + 1/2, k +
(4-18)
.n + 1 /2
1/2
2
pli,/-l/2,fc + 1/2
A
<
f>
P / + 1 /2 ~ P i —1 /2
r\ i J , k + 1/ 2 )
„
.n + l/2
_
At
__ . n —1 / 2
plt,y+1/ 2,*+1/2 ~
pl'i , ji + 1/ 2, k +
1/2 P|,-y+ 1/ 2, k+ 1/2
I71___________ E ^
‘Hi, j + 1/2,1:+1
(Hj,y+1/2, fc
Az
(4 .1 9 )
r
—E
\
'zli,y+ 1,fc+1/2
Hi, 7,* +1/2
P,(A<{>)
/
j-r . n + l / 2
« '<+
t>l
<pli+i/2,y,*+i/2
Tr , n —1 / 2
n*
= «
<fM,-+ i/2,y,*+i/2 +
X
I™
zlt+i,y,*+i/2
Ap
1
Ar
----------------
\k\i +l / 2JJc +
_
1/2
zI i, y, k +1/2
plt+i/ 2,y,fc+i ^ plt+ 1/ 2,y,*
Az
144
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4 .2 0 )
„ .n+l/2
<-1: . i /*> : . i
Ml
t-
i+\/2,j+l/2,k
— trrLi\I:, n - l / 2
z i + l / 2 , j + 1/2,
X
2 At
------------
->
fc
p|[ +
i / 2 7j + i / 2 , k
p f+ i
P;
-
EJn
- E ' n
Q\ i + 1/2, j + \ , k
p\i+W2,j,k
(P/ + 1- Pi) X
^
.
(4.21)
4.2.2 FDTD Calculations at p = 0
There are singularities in (4.18) and (4.19) in the above F D T D equations when p
approaches to zero. Thus, those two equations can not be used in this F D T D calculation.
Appropriate approximations to E AI
and H a I
are presented in this section.
v ip = o
ip = o
Assume that there are no impressed current and conduction current inside the region S
in Figure 4.1, we can obtain the finite difference equation for E z at p = 0 by Ampere’s
law as
n
Ez
r. n —1
1
, = E‘
0, j,k + 0,j,k + \
2
2
At
P
b
1
4 __ n 2
r .
t
1 ., 1
rl
rsl Jt k ***;
2’
J’
2
(4.22)
where r i is the distance along p for the first cell as shown in Figure 4.1.
Regarding to H A
P lp = o
However, E.\
Ip = o
, it is used to calculate E A
v lp = o
in (4.17) and E J
Ip = o
is approximated by Am per’s law in (4.22). I f E A
in (4.18).
v ip = o
is also
approximated then the calculation of H p at p = 0 is not needed. Traditionally, E^ at
p = 0 is approximated by its analytical solution. For example, E^ at p = 0 o f T E ^
mode is approximated by E^ at p = 1 times 1.0015 which is the ratio o f analytical
145
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
solution at these two points. This approximation is highly dependent on the mode being
calculated and difficult to be generalized. Other drawbacks o f this approximation are the
lack o f order control and the number o f partition along p cannot be too small even if the
Ap is much smaller than X / 10. These disadvantages can be solved theoretically by
applying the T y(2,4) F D T D to cylindrical coordinates FD TD .
4.2.3 Source Implementation
Traditionally, the way chosen to numerically incorporate an electromagnetic field
excitation source inside the cavity fo r the empty cavity simulation is based on the specific
cavity mode o f resonance. Implem enting a source is to select several lattice points as
source points and assigning magnitudes o f an electric field component at these points
based on theoretical cavity field solutions[4]. These source points are driven a few cycles
at a frequency close to the resonant frequency, then turned off. This technique o f exciting a
mode inside the cavity and then turning o ff the source gives the natural frequency
response. Examples o f excitation sources are shown in Figure 4.3.
The traditional source implementation is mode-dependent and actually gives wrong
field distributions which are shown in section 4.2.4 even though the fundamental modes
are calculated. However, by using the source implementation in (3.19) and BlackmanHarris function for g(t), correct field distributions can be obtained. In this chapter, the
second-order F D T D with traditional source implementation is denoted by T -F D T D and
that F D T D with Blackman-Harris source is called B H -F D TD .
146
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Perfectly conducting walls
tv
Excitation circle
1r
Excitation line
(Hr)
Excitation surface
E r_
(+)
A
*t
xi
-*
6.' } * Cv*'r / *1*F
**■*£.
k h r
f
*
'i
f
<1 -
< V*-
*
T M qu mode
Figure 4.3 T E 0n , T E lll9 T M 011, and T M 012 modes excitation techniques
in a cylindrical cavity[35].
147
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4.2.4 Num erical Results and Discussions
For T M q 12 mode, an empty cylindrical cavity with 0.0889 meter radius and
0.14409926 m eter height are used. The number o f partitions along p, <|), and z are 29, 36,
and 29, respectively. One time step is equal to 5 x l0 ‘ 13 second and the number o f total time
steps is 35,000. The total stored energies o f T -F D T D and B H -F D T D are plotted in Figure
4.4. The result o f B H -F D T D is closer to the correct one. The field distributions o f TM q i 2
mode o f T -F D T D and B H -F D T D are plotted in Figure 4.5 to Figure 4.8. For T-FD TD ,
only the field distribution o f Ep along z in Figure 4.6 and that o f Ez along z in Figure 4.8
are correct com paring to theoretic results. On the contrary, the B H -F D T D gives the correct
field distributions o f E field in all the figures.
Using an em pty cylindrical cavity with radius equal to 0.0889 m eter and height equal
to 0.0669089 m eter, a T E m mode can be excited. The number o f partitions along p, <(>,
and z are 59, 36, and 29, respectively. One time step is equal to lxlC T 13 second and the
number o f total tim e steps is 175,000. The total stored energies o f T -F D T D and B H -FD TD
of TE m
are plotted in Figure 4.9 and the field distributions o f E field are plotted in
Figure 4.10 to Figure 4.15. Observing from these figures, we found that all the field
distributions calculated by T-FD TD are all incorrect. However, the field distributions of Ep
along <|>, Ep along z, E^ along <J>, and E§ along z are roughly sim ilar to correct ones. Again,
the field distributions calculated by B H -F D T D are all correct comparing to theoretical
results. Hence, the traditional source implementation needs to be replaced by BlackmanHarris source when calculating the fundamental modes of empty cylindrical cavity with
PEC wall. Also, from our experience the Blackman-Harris source should be used for
148
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.5
eo
05
LLI
fZ 0.5
0.5
1.5
2.5
Time (psec)
1.4
1.2
CO
e
< 0.8
l5 0.6
S5 0.4
I—
0.2
12
14
16
18
Time (psec)
Figure 4.4 The plots of total stored energy o f T M 012 mode in a P EC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
im plem entation.
149
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.6
0.4
z=0.07205
0.2
-0.2
LU
-0 .4
-0.6
-0.8
0.01
0.02
0.04
0.0 3
0.05
0.06
0.07
0.08
p (meter)
0.12
0.1
0.08
2 0.06
z=0.07205
0.0 4
0.02
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
p (meter)
Figure 4.5 The variation of E p o f T M 012 mode along the p direction in a
PEC empty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
150
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.8
0.6
0.4
p=0.04445
c 0.2
-Q
-0.2
LU
-0 .4
-0.6
-0.8
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
z (meter)
0.8
0.6
0.4
p=0.04445
0:0.2
-0 .4
-0.6
-0.8
0.02
0.04
0.06
0.08
0.1
0.12
0.14
z (meter)
Figure 4.6 The variation o f E p of T M 012 mode along the z direction in a
PEC empty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
151
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CO
'c
ZD
z=0.07205
-2
-6
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.06
0.07
0.08
0.09
p (meter)
-0.1
-0.2
z=0.07205
-0 .3
C -0 .4
2 -0 .5
-0.6
LD
-0 .7
-0.8
-0 .9
0.01
0.02
0.03
0.04
0.05
p (meter)
Figure 4.7 The variation o f E z o f T M 012 mode along the p direction in a
PEC em pty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
152
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2 .5
1.5
3
0.5
p = 0 .0 4 4 4 5
-0 .5
Ui
-1 .5
-2
-2 .5
0.02
0.04
0.06
0.08
0.1
0.12
0.1
0.12
0.14
z (m e te r)
0.3
0.6
0.4
c: 0.2
Z)
p = 0 .0 4 4 4 5
<
LU
n-0.2
-0 .4
-0.6
-0.8
0.02
0.04
0.06
0.08
0.14
z (m e te r)
Figure 4.8 The variation o f E z o f T M 012 mode along the z direction in a
PEC em pty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
153
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0
0.5
1
1.5
2
2.5
3
3.5
4
Time (psec)
0.09
0.08
<o
•■= 0.07
CO 0.0 6
0.05
05
0 0.04
"O
0)
o 0.03
O 0.02
0.01
8
10
12
14
16
Time (psec)
Figure 4.9 The plots of total stored energy o f T E l n mode in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional
source implementation and the lower one by using the B H source
implementation.
154
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
p=0.04445
CO
z=0.033454
-2
-S -6
Ui
-8
-1 0
-1 2
-1 4
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.08
0.09
-0 .3 5
-0 .4
-0 .5
p=0.04445
z=0.033454
-
0.6
-0 .6 5
0.01
0.02
0.03
0.04
0.05
0.06
0.07
p (meter)
Figure 4.10 The variation o f E p o f T E 111 mode along the p direction in a
P E C em pty cylindrical cavity. The upper figure is calculated by using the
trad itio n al source implementation and the lower one by using the B H
source im plementation.
155
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
15
10
CO
p=0.04445
z=0.033454
LU
-5
-1 0
-1 5
50
100
150
200
250
300
350
250
300
350
<j> (degree)
0.8
0 .4
0.2
p=0.04445
z=0.033454
-0 .4
-
0.6
-
0.8
100
150
200
<j>(degree)
Figure 4.11 The variation of E p of T E m mode along the <{>direction in a
P EC em pty cylindrical cavity. The upper figure is calculated by using the
trad itio nal source implementation and the lower one by using the B H
source implementation.
156
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
p=0.04445
(O
-2
ill -6
-8
-1 0
0.01
0.02
0.0 3
0.0 4
0.05
0 .0 6
0.0 7
0.05
0 .0 6
0.0 7
z (meter)
^
-
0.1
-
0.2
-0 .3
-0 .4
LU
-0 .5
-
0.6
-0 .7
0.01
0.02
0 .0 3
0.0 4
z (meter)
Figure 4.12 The variation of E p o f T E m mode along the z direction in a
PEC empty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
157
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
<t>=160°
z=0.033454
2
-2
HI
-6
<>®
-8
0.01
0.02
0.03
0 .0 4
0.05
0.06
0.0 7
0.0 8
0.09
p (meter)
<j>=160o
-0 .0 5
z=0.033454
-
^
0.1
-0 .1 5
-
0.2
-0.25
-0 .3
-0 .3 5
0.01
0.02
0.04
0 .0 3
0.0 5
0.06
0.07
0 .0 8
0.09
p (meter)
Figure 4.13 The variation of
o f TE111 mode along the p direction in a
PEC empty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the BH
source implementation.
158
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
15
10
p = 0 .0 4 4 4 5
z = 0 .0 3 3 4 5 4
LU
-5
-10
-15
50
100
150
200
250
300
350
<j> (d e g re e )
0.5
0.4
0.3
0.2
co
0.1
p = 0 .0 4 4 4 5
z = 0 .0 3 3 4 5 4
-0.1
-0.2
-0.3
-0 .4
-0 .5
50
100
150
200
250
300
350
<> (d e g re e )
Figure 4.14 The variation o f
o f T E l n mode along the (j) direction in a
PEC em pty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
159
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
J_______________ I_______________ I______________ I_______________ I_______________
0
0.01
0 .0 2
0 .0 3
0 .0 4
0 .0 5
0 .0 6
0.07
z (meter)
-0 .0 5
p = 0 .0 4 4 4 5
-0.1
-0 .1 5
111
-0.2
-0.25
0.01
0.02
0.0 3
0.04
0.0 5
0.06
0.07
z (meter)
Figure 4.15 The variation o f
o f T E m mode along the z direction in a
P EC em pty cylindrical cavity. The upper figure is calculated by using the
traditional source implementation and the lower one by using the B H
source implementation.
160
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
lossless cavity calculation whether if the cavity is empty or not.
For a loaded lossless cylindrical cavity, the T M 012 mode is considered. The
dimensions o f this empty cavity are 0.0762 meter in radius and 0.15458 meter in height.
The numbers o f partition along p, <f>, and z are 2 9 ,3 6 ,2 9 , respectively, and one time step is
equal to 5 x l0 "13. A small cylindrical material sample with 0.0105 meter in radius, 0.0107
meter in height, and 2.5 in relative perm ittivity is placed in the center o f the cavity. The
field distributions o f E p and Ez are plotted in Figure 4.16 and Figure 4.17, respectively.
These field distributions are similar to those in Figure 3.14 and Figure 3.15. The ratio of
Ez/ E lz o f T M q 12 mode along the z direction in the PEC cavity loaded with a small
cylindrical m aterial sample is plotted in Figure 4.18. The theoretical estimation for this
case is 0.667 and the numerical results is about 0.75. This numerical result can be closer to
that in Figure 3.16 if the numbers o f partition along axes are increased and the dimensions
of the loaded small cylinder is less than tenth o f those o f the cavity.
Form these numerical results, we can conclude that the B H -F D T D with traditional
treatments at p = 0 gives correct field distributions o f a lossless cylindrical cavity.
However, the drawback o f B H -FD TD is its mode-dependent due to these treatments at
p = 0 . By using the Ty(2,4) FD TD method, a general treatment can be obtained and a
smaller number of partition along each axis is required for the same accuracy. The Ty(2,4)
FD TD scheme is then expected to be the fundamental FD TD algorithm for cylindrical
coordinates.
161
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0 .1 8
0 .1 6
0 .1 4
C/3
s
0.12
2
o.i
z=0.072137
LU
0 .0 6
0 .0 4
0.02
0.01
0.02
0 .0 3
0 .0 4
0.05
0.06
0 .0 7
0.08
p (meter)
0.2
0 .1 5
0.1
CO
HI
0.0 5
p=0.0127
0=0.05
-0.1
-0 .1 5
-0.2
0.02
0 .0 4
0 .0 6
0.08
0.1
0.12
0 .1 4
0.16
z (meter)
Figure 4.16 The variation o f E p o f TM o12 along the p and z directions in a
P E C em pty cylindrical cavity. The above figure is calculated by using the
trad itio n al source implementation and the below one by using the B H
source implementation.
162
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-0.1
-
0.2
—^ —0 .3
z = 0 .0 7 2 1 3 7
-0 .4
-0 .5
-
0.6
-0 .7
0.01
0.02
0.04
0.0 3
0.05
0.0 7
0.06
0.08
p (meter)
0.8
0.6
0 .4
-ti 0.2
p = 0 .0 1 2 7
n-0.2
-0 .4
-
0.6
-
0.8
0.02
0 .0 4
0.06
0.08
0.1
0.12
0 .1 4
0.16
z (meter)
Figure 4.17 The variation o f E z o f T M 012 along the p and z directions in a
P E C em pty cylindrical cavity. The above figure is calculated by using the
traditional source im plementation and the below one by using the B H
source implementation.
163
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.15
1.05
0.95
0.9
0.85
0.8
0.75'—
0.045
0.05
0.055
0.06
0.065
0.07
0.075
0.08
0.085
0.09
0.095
z (meter)
Figure 4.18 Plot o f the ratio o f Ez/ E lz o f T M 012 mode along the z direction in
the PEC cavity loaded w ith a small cylindrical m aterial sample.
164
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4.3 iy(2,4) Cylindrical FDTD Scheme
Applying the im plicit staggered fourth-order approximation,(2.57), to spatial
derivatives o f fields in (4.9) to (4.14), the Ty(2,4) cylindrical F D T D scheme is obtained. A
general treatment at p = 0 for Ty(2,4) cylindrical F D T D method is discussed in section
4.3.4. The matrix equations along p, <j>, and z are obtained in section 4.3.1, section 4.3.2,
and section 4.3.3, respectively.
4.3.1 Derivatives of Fields in p direction
There are four terms in (4.9) to (4.14) which involve the derivative along p direction;
8pt f z , 5p(
p
, 5 pE z, and 8p( p ^ ) .
The derivatives o f H fields are considered first. Applying (2.57) to every lattice points
in Figure 4.19, the following m atrix equation can be obtained
AX = M
(4 .2 3 )
where
1
22
1
0
1
22
•
i
0
•
0
0
0
0
A =
(4 .2 4 )
0
X =
0
•
_0
0
£
d /
dp
dp
1
22
1
0
0
0
1
22
1
0
1
22
d /
P= i
1 (Afp- l ) x ( t f p + l )
K
P= 2
dp
(4 .2 5 )
p = Nf
and
165
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
5P<PH»)|p=0
'
V P H» )U
8p(PH*)|p=2
1
P ^ l p - l 12
6P(P lM p = Np- l 8 P(PH^ p=Np
'-------P ^ l p=Np-3/2 P^<>| p=N!p- 1/2
P^P I p=3/2
Figure 4.19 The locations of H field and its corresponding derivatives
in p direction.
166
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
•^lp = 3 / 2
-^"Ip = 1 / 2
•^Ip = 5 / 2 ~ - ^ l p = 3 / 2
(4.26)
where N p is the number o f partition along the p direction and / is H z or p H ^ . In order to
use L U decomposition, a square matrix is desired. Hence, the derivatives o f H fields at
p = 0 in Figure 4.19 needs to be approximated, since (4.23) is a underdetermined linear
equation. To m aintain the fourth-order approximation for the overall F D T D scheme, this
approximation at p = 0 needs to be fourth-order at least. Both fourth-order one-way
difference and six-order one-way difference approximations are listed below:
Fourth-Order
3/
24Ap
(4.27)
Six-Order
(4.28)
where
167
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
( —5.6020833')
19.024479
-30.882813
O l5 b2, b3, bA, b5, b6, b-j)
=
30.369792
(4 -2 9 )
-18.026042
5.9640625
-0.8473983
The results by applying these two different order approximations w ill be discussed later.
Hence, the new m atrix equation is A X = M where X is the same as (4.25),
24 0
0
0
0
0
0
0
0
0
0
0
•
1 22 1
0
1 22 1 •
A =
(4 .3 0 )
0
0
•
1
0
•
1 22
0
•
0 1 22 1
_0
•
0 0
0 24 (tfp+ l)x (tfp+l)
and
a p3^/
A
>
P= o
“^lp = 3 /2
24
M = =Ap
-^ip = 1 /2
•^lp = 5 / 2 - -^lp = 3 /2
(4 .3 1 )
•^Ip =
N p
- 3 / 2 _ -^ Ip =
N p
- 5 /2
-^Ip = /Vp- 1 /2 _ -^ Ip = A /p - 3 /2
a
3 /
AP^
where ^r~
dp
P= ^p
. Another expression of (4.23) can be
is the approximation o f ^
p = o
p=0
168
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
obtained as:
-
•
0
0
1 22 1 •
•
0
0
22
1
(4.32)
A =
0
•
•
-
X =
3/
dp
p = ,i
•
1
•
1 22
•
0 1 22 OVp - l)x(JVp-
3/
3p
K P= 2
3f
3p
(4.33)
P = Afp- 1
and
f\
- f\
_ Apd/
J Ip = 3/2 J Ip = 1 /2
24 dp
24
M = =Ap
•^"Ip = 5/2
p=0
-^"lp = 3/2
(4.34)
•^Ip = /Vp- 3 / 2 ~-^lp = N p -5 /2
fl
- fI
_ AP d /
Ip = A^p—1/2 ■/|p=iVp- 3 / 2
24 3p
P = wp
The condition number o f a m atrix A measures how unstable the linear system A X = M is
under the perturbation o f the data M . In the F D T D calculation, a small condition number
which is close to 1 is desirable. The later m atrix equation is used in this chapter since the
condition number o f m atrix A in (4.32) is smaller than that in (4.30) for a same number of
partition.
For derivatives o f E fields along p direction, the 5 pE z and 8p(p £ (j)) at p = 1 /2 and
p = Wp —1 /2 has to be approximated regardless o f the physical boundary condition.
The locations o f the field and its derivative is plotted in Figure 4.20. The m atrix equation
169
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
is the same as (2.73) to (2.76) and is rewritten below:
AX = M
(4.35)
where
26
-5
4
-1
•
•
0
1
22
1
0
•
-
0
0
1
22
1
0
•
0
A =
(4-36)
_
X =
.
0
1
22
1
0
0
•
0
0
1
22
1
0
•
0
-1
4
-5
26
a/
aP
a/
(4.37)
i
a
i
ii
3 ap
p = * P- f
a.
3
p
H =2i ap P = 2
■i n 1
V
ap
0
and
-^Ip = 1 -^ip = o
•^ Ip = 2
M =
~^lp = I
24
Ap
(4.38)
■^Ip =
Np—\
f \ p = N P~ :
f \ p = N ~ f \ p = Np- l
Note that the p ^ l
in Figure 4.20 is equal to zero no matter what
vlp =o
approximations used for end points are
170
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
v
is. The
p=0
p=2
P=1
h Ez p=l/2
h Ez p=3/2
3=Np—1
h Ez p=Np—3/2
8 E
Er
(?=Np
p=Np—1/2
0
p £ <t>| p=0
p £ <t>I p=i
P^O. p=2
PE<
t> p=Np—1 p£<
t>I p=Np
1
5P(pE^ | P=i/2
V pE*)|p=;
p=3/2
1---------
5p(PE<t))| p=Np—3/2 5 p(PE<p| p=Np—1/2
Figure 4.20 The locations of E field and its corresponding derivatives in
p direction.
171
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
+i.v
123p
243p
P~ o
5
P“ O
5 df
248p
133/
123p
3
P —n
f\o = 1
f\o = o
Ap
(4-39)
p -i
and
13 8 /
120p
i a/
248p
63p
! 240p
P = Np- t
P = N p~ l
P = Afp-
P = ^ p- 4
(4.40)
f \ k = Np~ f \ k = Np- l
"
Ap
_
4.3.2 Derivatives of Fields in (j> direction
For derivatives of fields in <j> direction, there is no end points or boundary points;
hence, (2.57) can be used for all the lattice points along <{> direction. The 8^ H z and
are considered first. Applying (2.57) to the F D T D lattices in Figure 4.21, the matrix
equation for this case is
AX = M
(4.41)
where
22
1
0
•
1
1
22
1
•
0
•
•
•
0
•
.
1
22
1
0
1
22
A =
0
1
(4.42)
ir
X = V
8(J)
0= 0
3<()
.* 1
3<|)
0=1
_2 a<t> 0 = ^ - l_
and
172
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.43)
-^l<t> = ATp- l / 2
f\<l? = 3/2
= 1/2
1W
24
M = A*
(4.44)
= JV* - 3/2 ~ ^l<f> = Af, - 5 /2
•^I<1> = AT, - 1/2 ~-^l<t) = AT0 —3 /2
For &qE z and &$Ep , the matrix equation is
AX = M
(4.45)
3<j)
H
3
d) = v 2
d<l>
3 d<t>
<i>= Af0 - |
1
,N)| —
V
3(f)
-©ii
where A is the same as that in (4.42),
(4.46)
and
K*
24
M = A0
(4.47)
■^I<1) = A ^ -l f U
= Ns-2
f U = o - f U = » 0- i
173
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 4.21 The FDTD lattice along the <|> direction.
174
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
4.3.3 Derivatives of Fields in z direction
The calculation of derivatives o f fields in the z direction is sim ilar to that in 2.3.4 and
2.3.5. The matrix equation for derivatives o f magnetic fields, 8 .//^ and 8 . H p, is
AX = M
(4.48)
where
22
1
0
-
0
1
22
1
-
0
•
•
•
•
0
0
•
1
22
1
0
•
0
1
A =
X =
dz
z= 1
z
22 ( A r . - l ) x ( A T . - I )
V
M
dz t . N . - l dz
¥
dz
(4.49)
=2
(4.50 )
z = N .~ I
and
•/f \ lz = 3 / 2
- /JI \ z =
/U
M =
5/ 2
1 /2
- /U
24
=0
3 /2
24
Az
(4.51)
f \ z = Nz—3/2
f \ z = Nz- 5 /2
fJ \ \z = N , - l / 2 “ J/I\z = N .-2,/2 2 4 dz
For PEC boundary, the
dz
and^
dz
z = N.
are all equal to zero, but these quantities
z = Nr
need to be approximated by (4.27) or (4.28), respectively for FEC boundary.
For derivatives of electric fields,
and &zE p , the m atrix equation is (4.48) with
175
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
26
1
-5
4
-1
-
•
22
1
0
•
-
0
0
0
1
22
1
0
•
0
(4-52)
A =
X =
dz
0
.
0
1
22
0
-
0
0
0
•
0
-1
•
-
1 d -7
3
‘2
z= 2
0
I
1
22
4
-5
26 _ ALxU.
1
d_f
dz
(4-53)
z=
* .- §
z = N z- 1
and
/ U
r / U
/ U
2
- / U
o
1
24
M = =—
Az
(4-54)
f \ z = /V. —1 ^ l / = W.-2
■^lz = Ar. - -^lz = A/.-l
4.3.4 FDTD Calculations at p = 0
D ue to the singularity at p = 0 , special treatments need to be applied to those finite
difference equations which contain divergent terms. A t the first glance, the first two terms
on the right hand side o f (4.11) and the second term on the right hand side of (4.12)
become infinite or undefined when p approaches to zero. Hence, this two equations can
not be used for field calculations at p = 0 . Thus, the FD TD equations from (4.9) to (4.14)
at p = 0 are examined one by one in this section.
Equation (4.9) is still valid for p = 0 and is used to calculate D p which is located at
176
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
half-integer space step along the p axis. For p = 0 , this equation becomes
i
n~2
Pi
\,j,k
1
n~2
,
k j.k
1
—
2’
I
n~2
1
1
n —x
n —~
2
- / cp 1 2 - jJ s p
I*/,*
(4.55)
2 ’ J ’ k
J ’ k
Equation (4.10) is valid for p = 0 if 87H p at p = 0 is known. Here w e cannot only
at p = 0 and ignore 8 . H p at p = 0 since it is required in (4.51).
approximate
Moreover, H p at p = 0 is also required in (4.51); hence, H p at p = 0 needs to be
approximated which is discussed later in this section.
Equation (4.12) also has divergent terms so this equation can not be used for p = 0 .
Bp at p = 0 is approximated by the following equation
1 dBp
24
ap
_
1 * BP
12 a P
p=5
5 ^5p
5
P=2
24 a P
3
P= 2
13^P
12 a P
5J
P |p= l
1
P=2
-B r
PLp
Ap
=
= o
(4 .5 6 )
which is a fourth-order approximation. In order to obtain the derivatives o f Bp with
respect to p , both dBp/ d p at p = 1 /2 and p = N p - 1 /2 need to be approximated by
-9 3 B
35,
P = 1
3p
2295,
- 225B
PLP =
= 2
P
24Ap
+ 111 B r
= 3
-2 2 B
!p = 4
P Lp -= 5
(4 .5 7 )
P =
and
35,
= (-9 3 5
3p
p lp
-2 2 5 5 ,
+ 2295,
= /V p - l
P|p
= Arp - 2
p lp
= AAp — 3
P = /vp- 5
( 4 .5 8 )
1115p|
—225pI
p lp = yVp - 4
)/(2 4 A p )
M lp = Afp - 5
177
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
respectively. The matrix equation is A X = M where
A =
22
1
0
-
0
1
22
1
-
0
•
-
-
-
0
0
.
1
22
1
0
1
0
(4.59)
22 (JVp-2)x(JVp-2 )
IT
X =
dBp
dBp
dBp
dp
3 3P
p = i
dBp
3 P P = A/: - |5 3 P
5
p= i
r
p
=
3
(4.60)
z~ 5J
and
A pd£p
24 3p p = 1/2
BJ
-B r
plp =2
plp = l
24
M =
Ap
Bp|p = 3 - 5 p|p = 2
(4-61)
g |
PIp = N p —2
g
plp = JVp- 3
BA
-B I
“^ 7 3 ^
p lp = A/p- l
P|p = /Vp-2 24 3p p = A T - 1/2
.
where
dBp
a P p = 1/2
a P p = /V0- l / 2
are the approximations o f
at p = 1 /2 and
p = ATp — 1 / 2 , respectively. The values o f 5 p at p = 1 to p = iVp - 1
can be
calculated by (4.12) first, vector X is solved, and then Bp at p = 0 can be obtained. The
value o f H p at p = 0 can be obtained by applying the constitutive relation.
Equation (4.11) definitely cannot be used at p = 0 due to the divergent and undefined
terms on the right hand side. Assume that there is no impressed current but with a
178
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
conduction current inside region S in Figure 4.1, we can obtain the finite difference
equation for D . at p = 0 by Ampere’s law as
1
n~2
D z\n + D z\n 1
+ a_s_-
(4 .6 2 )
0,j,k + ±
2>j’ k + 2
where r { is the distance along p for the first cell. A more general approximation can be
obtained by sim ilar approximation which is presented in the approximation o f H p at
p = 0 . Replacing H 0 by E
E
from (4.56) to (4.61), the fourth-order approximation of
at p = 0 can be obtained.
Equation (4.13) is valid for p = 0 and it becomes
1 .,
1 = Vp * z
2’
2
J ’
-5
(4 .6 3 )
2 , J, k +
2 ’ J' k + 2
2
Equation (4.14) is valid for p = 0 , and it becomes
1
StB z
2’ J + 2’ k
pi • u
2 ’ J + 2’ k
V
p
Sp(P E J
’
2’ J + 2’ k
+
i . i ,
2’ J + 2
(4 .6 4 )
where the location o f p in 5p(p E (f)) depends on that o f E ^ ,
4.4 Time Domain Finite Difference Equations of Constitutive Relations
The electric perm ittivity used in this chapter is modeled with the general Debye or
Lorentz equations.
The Debye equation is
179
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
~ sk
e(CD) = e'(Qi) - j e " ( ( o ) = e'^+ £
-
_ £ _ J L ,
(4.65)
and the general Lorentz equation is
M
e(co) = e’(co)-ye"(co) = e’^ + (er - e’^) J — ---- 1 - ^ ---where
(4 .66)
co^ is the £-th resonant frequency, x>k is the fc-th damping frequency, e’^ = e (°°),
M
e’o = e(0)
and £
G k = 1.
k=
1
.Jk.
The scalar constitutive relations modeled with Debye relation between E and D in
time domain can be obtained as:
M
= e’a ~ £ a ( 0 + X
P a * (0
<4-67)
k = 1
where
P akW
and a is one o f
dPa* ( 0
+ T’a.ek
0^ — ( e ask ~ e
(4 .68)
p , (j), or z . B y the same way, the scalar constitutive
relationsmodeled
with Lorent equation in time domain is
M
~ e a°°p a ^ + (e as —e a~)
Pc tk ^
(4 .69)
k= I
where
d2P a t ( 0 ,
0 ?2
+
3 P ak(t)
,
0^
+ ®a.kP a.k(t}
G<xkaia.kp‘a.(t)
180
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4-70)
2
2
and only the damping case w ith x>ak —4coajt < 0 is considered.
The constitutive relations for Debye material in finite difference form are
M
»af =
+ £
P a tf
(4-71)
k= I
and
(4.72)
and that for Lorentz material are
M
(4-73)
k= 1
and
(4.74)
sf Pat
I f the second-order tim e stepping is used, then (4.72) and (4.74) can be rewritten as
2Ar _
a/fc
,n — l
,n —2
a;fc
2Ar
~
ua ek
ask
(4.75)
ccek
and
= ( - 2Y - V ® L ) ^ | ' " I + ( l + 5 y ) / > ait|', - 2 - 4 Y P a t | " - 3 + Y/>a t | " - ‘‘
,n — 1
Y^ a k ^ a k ^ a
where y = 2 / ( t ) ajfcA r ) .
181
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.76)
C H A PT E R 5
CONCLUSIO NS
In this dissertation, the finite-difference tim e-dom ain(FDTD) method is employed to
quantify the induced electromagnetic field in a material sample placed in an energized
microwave cavity. Due to the discrete nature o f FD TD methods, the stability condition and
the numerical dispersion becomes major problems when it is applied to E M problems.
Generally, smaller Courant-Friedrichs-Lewy(CFL) values lead to the more stable F D TD
scheme. However, smaller C FL values w ill cause a larger phase error which is not
desirable in the F D TD calculation. Hence, the maximum C FL value for a stable FD TD
scheme which leads to a minimum phase error and also demands a smallest number o f
grid cells in one wavelength is desired for F D T D calculation. This optimization can be
achieved by the analysis o f the numerical dispersion equation for the F D T D scheme which
is demonstrated in section 2.2.2.
The traditional second-order FD TD method by Yee in rectangular coordinates is
shown to be an efficient solver for closed boundary rectangular cavities if the object within
several wavelengths is considered. For electrically large objects or cavities, the larger
number o f grids and phase error w ill dramatically slow down the calculation and destroy
the numerical accuracy. Either higher-order scheme or other methods need to be
introduced to deal with the electrically larger objects. Higher-order FD TD schemes reduce
182
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the number o f grid i f the same phase error for the second-order is required. However, the
explicit higher-order F D T D scheme not only increases the stencil with more field points
but also complicates the treatment o f lattice points near the physical boundary. The
im plicit higher-order F D T D schemes uses the same number o f field points as that in Yee’s
algorithm. The increased calculation time by the m atrix multiplication in im plicit higherorder F D T D scheme is compensated by the decrease o f the number o f grids in one
wavelength. The Ty(2,4) method, which employed the second-order approximation in
time stepping and the im plicit fourth-order in the spatial stepping, is shown to lead to an
easy implementation in section 2.3.4 and section 2.3.5. Although the dimensions of
rectangular cavities under consideration are close to the wavelength o f microwave, the
Ty(2,4) is highly applicable to higher frequency E M problems or microwave problems
with small size devices.
The quality factor, Q, and the resonant frequency are two important measurement
factors o f cavities. By using the time-domain Poynting’s theorem derived in section 2.4.2,
Q value can be obtained from the FD TD calculation and (2.116). The resonant frequency
can be obtained by applying the fast Fourier transform (F F T ) to field data verse time. For
lossless cavities case, the F F T approach may be practical since we can control the time
step to turn o ff the source. For lossy cavities case, the F F T approach may become
impractical since the number o f time steps to reach the steady state can be very large. On
the other hand, both the quality factor and the resonant frequency o f a cavity can be
obtained by using the Prony method in section 2.5. The numerical results in section 2.8
shows that the prediction o f Q and the resonant frequency by Prony method is close to
those by time domain Poynting theorem and the FFT.
183
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The singularities introduced by the FD TD method at the center o f a cylindrical cavity
is a major problem in the application o f the F D T D method. For cylindrical cavities loaded
with a symmetric material sample, the BOR F D T D method is the most efficient method to
perform the FD TD calculation as shown in chapter 3. The treatment o f singularities in the
B O R FD TD method is also discussed in section 3.1.3. Since the B O R FD TD method is a
2.5D FD TD scheme, the calculation is much faster than that in a 3D case. Hence, the
number o f partition can be set to a much larger number and more detail around
discontinuity can be obtained. A conclusion can be drawn from section 3.2.1 is that the
PEC boundary can be assumed in most cavity calculations w ith the metal w all having a
conductivity larger than 104.
For general cylindrical cavities loaded with samples o f arbitrary shape, the 3D
cylindrical F D T D method needs to be employed. The traditional second-order cylindrical
3D FD TD method suffers the problem o f the mode-dependent source implementation and
the treatment o f sigularities. By using the Blackm an-Harris(BH) source, it solves the
source problem in the traditional second-order 3D F D T D method. In order to obtain a
general 3D FD TD method, a fourth-order treatment for sigularities is proposed in section
4.3.4 to replace the mode-dependent treatment in second-order F D T D method. Combining
w ith the im plicit staggered fourth-oder F D T D method, this proposed Ty(2,4) F D T D
method in cylindrical coordinates solves those problems in the traditional second-order
F D T D method. However, the implementation o f this Ty(2.4) F D T D method in cylindrical
coordinates suffers a numerical instability at this point.
The techniques developed in this study including the Ty(2,4) scheme in rectangular
and cylindrical coordinates and the time-domain power analysis may be useful for a wide
184
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
range o f E M problems. Some topics which are relevant fo r further study in the future are:
(1)
Optimization o f the numerical dispersion equation.
(2)
Power analysis o f other material models; e.g., Lorentz m aterial.
(3)
Extend the Ty(2,4) F D T D method for magnetic material.
(4)
Investigating the implementation o f proposed Ty(2,4) F D T D method in cylindri­
cal coordinates and the analysis of the corresponding numerical dispersion equa­
tions
(5)
Incorporate with thermal equations in the Ty(2,4) forms with a Ty(2,4) FD TD
method, and perfect matched layer (P M L) for the T y(2,4) F D T D method.
185
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX A
DERIVATION OF FO U R T H -O R D E R FIN ITE
D IFFE R E N C E APPRO XIM ATIO N
In section 2.3.1, several fourth-order finite difference approximations to spatial
derivatives are presented. This appendix details the derivation o f those finite difference
approximations based on the Taylor’s expansion. The explicit scheme have only one
unknown in the approximation; on the contrary, the im plicit scheme have more than one
unknown in the approximation. The locations o f the field point and its derivative also play
a role in the approximations. The collocated scheme has the field point and its derivative at
the same location; on the contrary, the staggered scheme has the field point and its
derivative at different locations with a half step spatial difference. Generally, three field
points are involved in the second-order approximation and five points, including field
points or their derivatives, fo r the fourth-order approximation.
For explicit collocated scheme, the following Taylor’ expansions are used:
a,_,
=
+
+
186
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A .2 )
ut ^
= u, + (2Ax) ■«(!) +
• Bp> +
• « p ) + ^ 2 ^ • «<*> + . ..
(A.3)
«, - 2 = a , - ( 2 A x ) - I<( l) + 2 | ^ - « p ) - ^ 2 | £ ) ! . „ p ) + ( 2 | £ ) f . i<(4) + . . . .
(A.4)
M ultiply (A .l) by a , (A .2 ) by 3 , (A .3 ) by y , (A .4 ) by 6 and sum all these new equations
together; then setting the sum o f the coefficients fo r
mU )
u p ) , « p ) to be zero and that of
to be one. The follow ing equation is obtained,
r
a + 3+ y+5=0
a - 3 + 2 y -2 5 = 1
(A5)
a + 3 + 4y + 4 5 = 0
(. a - 3 + 8 y -8 5 = 0
The solution for (A .5 ) is a = 8 /1 2 , 3 = —8 /1 2 , y = —1 /1 2 , and 5 = 1 /1 2 . Hence,
the explicit collocated approximation is
fdu \
8 ( « 1.+ 1 - i i |. _ 1) - ( k |. + 2 - k i- _ 2 )
U ; J , S ------------------- n K i -------------------- •
(A6)
For the explicit staggered scheme, the following Taylor’s expansions are used:
Ax m
-J ’
“ «+ 1/2 =
ui +
u i - 1/2 =
ui~ ~ 2
( A x /2 ) 2
}+
2!
Ax
3 Ax
+ 3/2 = «/ + —
3Ax
ui -
3/2 = ut - —
m . ( A x /2 ) 2
' ui
rn
21
m
( A x /2 ) 3
' U‘ + " 3 !
f2)
' U‘
f3>
‘
( A x /2 )3
~
31
( A x /2 )4
+
4!
'
m . ( A x /2 ) 4
'
+
4!
+ -
(A‘7)
m .
' “ P } + '• •
( A ' 8)
(3 A x /2 )2
m
(3 A x /2 )3
m (3 A x /2 )4
m
2f
‘ “ «C * +
3r '
• “P }+ — 4;
• «.(4) + ’ <A'9)
* “ t( I) + —
m . (3A x/2)2 m
(3A x/2)3
m . (3Ax/2)4
m , ,A
■ u\ l >+ *—
’ ui ~ — 31 '
‘ ui ’ + -— 41— * ur } + -(A-10)
Use the sim ilar procedures in the derivation o f the explicit collocated scheme, the
187
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
following set o f equations is obtained
a + P + y + 8=
1
2
1Q
- a - - B
2
3
2
0
3S
2
+ - v - - o =
,
1
( A .ll)
- a + - 0 + -y + -8 = 0
4
4
4
4
1
8
l fi
8
27
8
"
The solution fo r ( A .ll ) is a = 2 7 /2 4 ,
"
27 cj
8
n
"
3 = —2 7 /2 4 , y = —1 /2 4 , and 5 = 1 /2 4
and the explicit staggered scheme is
du\
2 1 ( u i+
1 / 2
^dx)[,-:
(S
~ ui - i / 2 ) - ( ui + 3 / 2 ~ ui - 3 / 2 )
------------------------24 Ax
(A . 12)
For the im plicit collocated scheme, the following Taylor’s expansions are used:
«,+ 1 = a . + A x .B p ) + £
^
. a (2) + ^
! . Bp ) + ^
.
a|<)+ ...
(A. 13)
(A . 14)
Mt + I =
u il\
+ Ax • Mp) +
• uf3) +
• uW + ~ y f“ •
= Mp ) _ A x - 4 2) + ^ ^ - « p ) - ^ ^ - Mp ) + ^
^
+ •••
(A . 15)
- Mp ) + . . . .
( A. 16)
The corresponding equations for a , 0, y, and 8 are
188
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
a + 0= 0
(a ~ P )A x + y + 8= 1
(A . 17)
( | a + |p )A x + (Y -8 )= 0
( ia - |p ) A x
+ ( I Y+
i8 )= °
and the solution is a = 3 /(4 A x ), P = —3 /( 4 A x ) , y = —1 /4 , and 5 = —1 /4 . Hence,
the im plicit collocated scheme is
(A . 18)
For the im plicit staggered scheme, the following Taylor’s expansions are used:
u/ +
u
1/2
-
ui
21
Ax
m , (A x /2 )2
1-1/2 = “ i - ~ ' “ i(1) + ""2 !
'
4!
3!
m
( A x /2 ) 3
m , (A x /2 )4
m
3!
U‘
+
4!
'
U‘
+
3!
= « ll) + A x- w(2 ) + (Ax)2
2 , • „k tP> +
3!
- «C4) +
4!
• «(5) +
4!
The corresponding equations for a , P, y, and 8 are
189
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A . 19)
(A .2 0 )
(A .2 1 )
(A .2 2 )
a 4-(3= 0
g a - I P)A * + y + 5= 1
(A.23)
[ | a + |p ]A x + ( y - 5 ) = 0
and the solution is a = 1 2 /( llA x ) , P = —1 2 /( ll A x ) , y = —1 /2 2 , and 5 = - 1 /2 2 .
Hence, the im plicit staggered scheme is
fir)
\ o x J i + 1 +(r)
\ d x j i_ l
n / < M _ ui+ 1/2 ~ “ f —1/2
Ax
+ o la x j,-
24
(A .2 4 )
The fourth-order approximation at the points half-integer grid away from the boundary
in (2.71) and (2.72) can be obtained by assuming
(A .2 5 )
Using the follow ing Tayor’s expansions:
Ax
“ i/2 = “ o + ~2 ‘
m
}+
(A x /2 )2
2!
'
(2\
+
( A x /2 ) 3
/3n (A x /2 )4
f41
3;
‘ “ i 3) +
4,
• 4 4) + •••
«{D = 4 i ) + A x - « p ) + ^ ^ - Mp ) + ^ ^ - M|4) + ^ ^ - « ( 5) +
3 Ax
“ 3 /2 =
k5 / 2
.
“o + —
„0
m
(3 A x /2 )2
• 4 l) + s
-2 i
m
•
+ S|£. 4 ., + (5A|2L2.
u k
}
+
+
(3 A x /2 ) 3
-—
5.
•
3!
m
4 3) +
.4
(3 A x /2 )4
a
- •
4!.
...
(A-26>
(A .2 7 )
f4,
4 4 ) + •• • (A -28 )
+ (5 ^ p .4. 4 4 , + ...(A.29)
3>
and
190
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
7A x
m
“ 7/2 = “ 0 + - 9 - • 4 l) +
(7 A x /2 )2
m
( 7 A x /2 )3 m
(7 A x /2 ) 4
m
7,
• 4 2) + V— l i
• “ &3) +
• “ ^4) + • -.(A.30)
the solution is a = 1 3 /1 2 , P = - 5 /2 4 , y = 1 / 1 2 , and 5 = - 1 / 2 4 .
The one-way fourth-order finite difference approximations in (2.81) and (2.82) can be
obtained by assuming
d ll\
(
5“
„o x J q
a “ l / 2 + P “ 3 / 2 + y 115 / 2 + ^ “ 7 / 2 + ^ “ 9 / 2
= -------------------------- A
a-------------------------------x
(A-31)
with (A .26), (A .2 8 ), (A .29), (A .30) and
9 Ax
“ 9 /2 = “ 0 +
rn
(9 A x /2 )2
m
‘ uf> } + — 2!
'
(9 A x /2 )3 m
— 3!
’
(9A x/2)4
r4>
+ — 4!
*
+ - (
}
The solution is a = - 3 1 / 8 , p = 2 2 9 /2 4 , y = - 7 5 / 8 , 8 = 3 7 / 8 , and r\ = - 1 1 /1 2 .
191
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX B
DERIVATION OF N U M ER IC A L D ISPE R SIO N
RELATIO N FOR IM PL IC IT STAGGERED
SCHEM E
In (2.70), the numerical dispersion relation o f the im plicit staggered F D T D scheme is
presented. The derivation o f this equation provides fundamental derivation o f those matrix
equations in Ty(2,4) F D T D scheme. In this appendix, the derivation of (2.70) is detailed.
Consider M axw ell’s equations in a normalized region of free space with p. = 1,
e
=
I , a
= 1 , and c — 1 and obtain
-v
jV x V =
ay
where '& = !% + j £ and
Let
^
^ j ( K j A x + KyJ A y + K zKAz + OnAt)
192
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(B .l)
where K xr , K yv , and K_z are the wave numbers along x, y, and z axes, respectively in the
computational space.
Consider the x component o f (B .l) at point (/, J, K) which is
dVz
J\ dy
d V y\
(B-4)
dz J
dt
I,J, K
I,J,K
The second-order finite difference approximation o f the time derivative o f V x is
n+ l/2 _
dV,
dt
n-l/2
Il , J , K
I/,
sin(03A £/2)
1
At
J, K
At
l,J,K
(B.5)
The time derivatives o f V y and V z can be obtained by replacing index x w ith y and z. For
8 Vz/ d y component, Vz is assumed to be located at the half integer position in space and
d V z/ d y
is located at the integer position in space. The one-way fourth-order
approximation for numerical boundary used is assumed to be
dV ,
dy j = Q
= aV.
.
J
1
2
+ bV„
+ eV,
1
2
J
2
J
(B.6)
2
for j = 0 . The numerical approximation for j = N y can be obtained by sim ilarity and is
listed below:
d V.
dy
= dV.
i = N..
i
j =
,
+bV7
z
3
J = t f. -x
+cV7
z
5
+dV7
J=N .-=
z
7
j = N -L
2
(B.7)
+ eV,
Apply the im plicit staggered four-order finite difference approximation, (2.57), to every
lattice point along y axis, a set o f equations is obtained as follow:
193
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1 dV
24 dy
j.o
z=i
12
I d V ,4
j =2
24 3y
1 avz
ll3Vt
*■
VJ
- V
zly=3/2
*•!/ = 1/2
Ay
!dVz
lld V *
_
24 dy y =1 + 1 2 3 > J = 2 + 2 4 3 y y = 3
V I
-V |
zly=s/2
zly=3/2
Ay
.(B.8)
ll^ V
1 3V.
y = ^ v- 3
12 8y
1 3VZ
l3V z
l l dV z
i
12 * y
24 3y
!
v.\
-v ,\
4y = Ny~ 1/2
Hy = Afy- 3 / 2
£
= n, - 2
_
II
j
j = Ny- l
VI
- yj
Hy=A
T
y-3/2
zly=Ny- 5 / 2
A^
n
24 dy
j = Hy- 2
II
24 3y
i a * '.
24 3y
Combining with the numeric boundary condition in (B .6), the above equation can be
rewritten in a m atrix form of
(B.9)
AX = B
where
24
0
22
24
0
1
24
•
•
•
i
i
A =
X =
dV.4.
dVz
dy y=o * y
0
1
24
0
0
•
0
0
•
0
22 1
24 24
1 (Afy+l)x(Wy+l)
0
dvz
dy
dvz
dvz
j = N y~2 dy j = N y~ i dy J = V y
194
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
bV_
aV,
7
. 3
/ =
2
J
+ CV.
2
. 5
7
+ </V,
2
. 7
7
+ eV_
2
. 9
7
2
V zl
Ij = 3/2 " V z\■/ = 1/2
Ay
VI
- V I
Z1j = 5/2
zly = 3/2
Ay
B =
. (B.12)
r I
-V |
zl/ = Ny- 3 / 2
Z‘j —Ny —5/2
Ay
+ bV_
aV,
1
sT
ii
3
- V.
zly = Ny—3/2
Ay
+ cV_
z
+ dV7
1
z
1
toi^j
z|y = wy- i / 2
II
V
+ eV_
y = vy
The m atrix 5 is rewritten as follow
a
b c d e
-1
1 0 0 0
H -3 /2
0 - 1 1 0 0
B =
7 -M V
Ay
Ay
0
0 - 1 1 0
0
0 0 - 1 1
e
d c
b a
V,
(B.13)
\j = Ny- 3 / 2
\ j = N y- l / 2
where M is a (Ny+ l ) x N y matrix and V is a vector w ith Ny dimension. Hence, the vector X
can be obtained as
X = -?-(A lM ) V =
Ay
Ay
For a point with j = I , the derivative o f V , is listed below:
195
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(B.14)
Ny- l
dV,
dy
[ , j = l,K
(B.15)
= Av
- L x mv v z
-Vp = 0
where /n/p is the element of M
at Ith row and [3th column. According I
(B.3),
VJ
can be rewritten as follow:
Hy = (2P +l)/2
V,
VJ
(B.16)
-e
Z\ [ , j = l , K
Hence, (B.15) becomes
av•
_ vk i~ i.K
ay i, j = i , k
(B.17)
--------------------2^
mi&e
Ay
p= o
Sim ilarly, we obtain
V I
dV^
^Z
v z- i
V —
= — E T 5 • 2 m<ne
y \[,J,k = q
I,J,k = q
(B-18)
y=0
Hence, (B.4) leads to
Ar
V z, I
2sin (co A f/2) + - % M .
Ay
- 1
£
N y-
- j \ Ky[J ~
L 1
H Ay]
2 j J
P = o
(B-19)
_ / Kr ^ _ 2 Y iiw i
Y
m ,„e
L1
2 J j = 0
Similarly, the y component and z component o f (B .l) at point (I, J, K) are
196
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
y U,J,K
At
2sin (co A f/2) H
v
A
,k
u
%
^
Az
l_
mKye
Y —0
V
(B.20)
-
S mr<xe
a=0
Ax
= 0
.
and
l,J, K
At
2 s in (o )A r/2 ) + v A u , k
Ax
2
mio.e
a=0
(B-21)
Ary- 1
a7“ '
S
3=o
OT-/Pe
= 0
Equations (B .19), (B .20), and (B.21) can be rewritten in a matrix form ,
AV = 0
(B.22)
where A is 3x3 m atrix as shown in (B .27) and
x \l,J,K
V =
y \i,j,K
(B.23)
I/, y, k
B y taking the determinant o f A to be zero and restoring the light speed into the time related
terms, the dispersion relation is obtained as
197
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
v_
r A L -l
X
2sin(C D A r/2)~j2 _
[
cAt
mrae
X
a=0
J
mJ$e
2_=o_____________________
Ax
Ay
(B.24)
X
mKye
U zO _____________________
Az
Assuming uniform partition along x, y, and z axes, A = Ax = Ay = A z, the
dispersion relation becomes
2
[ X - 1
+
2 , mtae
\
\
a = 0
----
w
X
“
^
J
o
X
II
GO.
12r - ( an \
sm^
I '
J3RJ-
1
-^ t_
X
(B.25)
- j [ K( K- * ± ± y ]
m KYe
y=0
where a = J Z c A t / A which is the C FL number and R = A ./A which is the number of
cells in one wavelength. L e t it = KjCx + Kyy + k Jz , then (B .25) can be further simplified
as follows,
1
N z12r
. (
an
«2Lsml ^ J -
X
.icsin9cos0(2/ —2a —L).
-J
R
mlae
La = 0
nyVy- 1
_ .icsin9sin<l)(2/ —23—l)-,2
V1 7
R
mJ$e
Lp = 0
r N, - 1
.kcos9(2AT-2y —1)-, 2
X
mKye
2
l
R
y =0
198
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(B.26)
II
8 .
on
.
M "—
31
to
®°.
5’
S
R
«>
*03
J*
>
tto
\
to
w
B
■CO
+
+
o« M
>
*•>*.
:*:
[
H
i
—
wTv
I
» M
o
SI
to
co.
3*
s'
5:
i
*
■CO
>
*•>*.
*
i
\
to
to
•co
+
8 . . ^
« M "—
o
fco. . 5:
« M ".
to
CO
(B .2 7 )
8
■CD
C*
t;
8
>
w
\
I
to
>
Si
JK
I
to
■CO
° +
to
%8
199
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
where
_
A.ic
k
= — —.
2
vp/ c
K = —
(B.28)
Hence, the ratio of wave propagation speed in the computational space to that in the
physical space is equal to the ratio o f n to ic.
200
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX C
DERIVATIO N OF EQ UATIO N (2.107)
In section 2.4.2, the time domain version o f Poynting theorem is derived and g x(t) is a
key factor in calculating stored and dissipated power in time domain. In this appendix, the
derivation o f tim e domain expression o f gx(t) is detailed.
The function, gx( t ) , is the inverse Fourier transform o f E x( o o ) / [ l + (G )xec) 2] and can
be rewritten as follow:
(C-I)
r
zM
where E x(t —x) is a causal signal; that is, E x( t —x) is equal to zero when t < x and the
system is a causal system; that is, e ^ /T,X is equal to zero when x < 0 . Discreting (C .l) by
using the approximation in Figure C .l, the discrete version o f gx(t ) is
201
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
—m At
—(m ■+• l)Ar
e
mAt
-mAt
nAt
\
\
/ mm
lllllp
e ex E x(n —m)
E x(n —m — 1)
(m + l)A t
T “
^
“
E x( n - m - 1)
/
T“
------------------- -(m+ 2)At
e
“
E v(n —m —2)
ll
(III
mAt (m + l)A t (m +2)A t
F IG U R E C .l Approxim ation of integrant in (C .l)
202
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
mAl
n —1
*(-) =
(m + 1 )A t
E J n —m — 1)
e T“ Zix(n —m) + e
Ef
m = 0
mAt
n —1
1 r i At
e X'z E J n —m ) + e
2x
^ y
exm = l
Z"
E J n —m —I )
(C.1)
At
1 Ar
E c(n ) + e ezE x(n—l )
2 t ex 2
A/
Ar
= e T“ g ( n - 1 ) + £ ~ E x(n) + £ - e Z'zE x( n - l )
^ l ex
ex
A more accurate recursive approximation can be obtained as follows by using three
field points in Figure C .l:
mAt
(m + l) A r
ti — 2
g(n) =
2xr
k
3
5 >
m= 0
Z" E J n —m) +
T“
Zsx( n - m - l )
(m + 2)Ar
+ xe
Z‘z
E x( n - m - 2)
_
n —2
(m + l) A t
mA/
1
Ar
2x
^ T3 e
“ ?
t
2^
exm = 1
T.
cz Ex(n —m) + 4 e
E J n —m — 1)
(C.2)
(m + 2)A t
+ e
Z‘z
E J n —m —2)
At
2Ar
1 At
E J n ) + Ae Z'zE x( n - l ) + e ZezE x( n - 2 )
^ex 3
At
At
= e
/
*v
*-, / v 2 Af
g ( n _ i ) + _ _ £ x(n ) + - — £
°''e x
ex
where (C .2 ) needs to store one more history of E
2 At
A/
E x(n -1 ) + — e
ex
than (C .l) does.
203
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
p / ^\
Ex( n - 2)
APPENDIX D
DER IV A TIO N OF EQUATIONS (3.7) AND (3.8)
The term V x E in (3.2) is
Vx
^
m
=
(e Mcosm<{> + £vsinm<j))
0
oo
=
^ [V x (e vsinm<j>) + V x ( l Kcosm<|>)]
m- 0
oo
=
[(sin/n<j)V x ev + Vsinm<(> x ev) + (cosm<t>V x eu + Vcosm<t) x eu)]
m= 0
OO
=
^sinm(j)^V x ev —^<j> x e^j + cosm(j)^^<f) x c v + V x
m = o
and those terms in the right hand side of (3.2) lead to
2 r r
'j
= —
oo
^
^
(hucosm§ + frvsinm<t>)
m= 0
(hucosm§ 4- /i„sinm<|>)
m = 0
oo
=
(D.2)
_ -X
X c o s rn ^ m= 0
^
oo
+ £
sinm<t>^—[M']0 j V +
204
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
■
From (D .l) and (D.2), we obtain (3.8),
. Ttl J
i.
_
x
r
~ ,dhuv
r
*-!?■
±-<j> x ev u + V x e u v = - [ f x ] - ^ - + [ o * ] h UyV.
(D.3)
For V x H in (3.1), the follow ing equation is obtained,
^
^
(/i„cosm({) + /ivsinm(j))
■ ^
Vx
>
m = 0
=
^ [V x (huCosmfy) + V x (/zvsinm(f>)]
m= 0
(D.4)
^
^
^
^
= 2 ^ [(sinm<{>V x hv + Vsin/n<() x hv) + (cosm(J)V x hu + Vcosmtj) x hu)]
m= 0
oo
=
^sinm<))^V x h v ——<j>x hu j + cosm<p^—(j>x hv + V x
m= 0
and those terms in the right hand side o f (3.1) lead to
dE
°° fOBlt
'dt„
-*•
[ e ] - ^ - + [ct]£: + / j = [e]
^
dt,
d&v
\
cosm({) + - ^ sinm(j)J
m= 0
OO
+ [a ] ^ (e Mcosw<{> + !„sinm<t>)
m= o
(D.5)
V
( 7 Mcos/n<{> + y\,sinm<t>)
m
=
=
0
de„
^ I|^cosm<|)^[e]|^“
cosm<j)l [e.
+ [a]e„ +
m= o
_3e
+ sinm <t^[e]|^y + [ cj] I v +
Comparing (D .4 ) and (D .5), we obtain (3.7),
()&
x hVt u + S7xhUi v =
+ [CT]I«. - + j u.
205
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(D.6)
BIBLIOGRAPHY
[1]
John, D ., Jr. Anderson, Computational Fluid Dynamics: The Basics With
Applications, M cG raw -H ill, New York, N Y, 1995.
[2]
K . S. Yee, “Numerical solution o f inital boundary value problems involving
M axw ell’s equations in isotropic media,” IE E E Trans. Antennas and Propagation,
vol. 14, no. 3, pp. 302-307, 1966.
[3]
A . Taflove and M . E. Brodwin, “Num erical solution o f steady-state electromagnetic
scattering problems using the time-dependent M axw ell’s equations,” IE E E Trans.
Microwave Theory and Techniques, vol. 23, no. 8, pp. 623-630, 1975.
[4]
D . H . Choi and W. J. R. Hoefer, ‘T h e finite-difference time-domain and its
application to eigenvalue problems,” IE E E Trans. Microwave Theory and
Techniques, vol. 34, no. 12, pp. 1464-1469, 1986.
[5]
A . Navarro, M . J. Nunez, and E. M artin, “Finite difference time domain FFT
method applied to axially symmetrical electromagnetic resonant devices,” IE E
Proc., vol. 137, no. 3, pp. 193-196, June, 1990.
[6]
A . Navarro, M . J. Nunez, and E. M artin, “Study o f T E 0 and T M 0 modes in
dielectric resonators by a finite difference time-domain method coupled with the
discrete Fourier transform,” IE E E Trans. Microwave Theory and Techniques, vol.
39, n o .l, pp. 14-17, 1991.
[7]
F. Liu, I. Turner, and M . Bialkowski, “A finite-difference time-domain simulation
o f power density distribution in a dielectric loaded microwave cavity,” Journal o f
Microwave Power and Electromagnetic Energy, Vol. 29, p p l38-148, 1994.
[8]
C . Su and J. Guan, “Finite-Difference Analysis o f Dielectric-Loaded Cavities
Using the Simultaneous Iteration o f the Power Method w ith the Chebyshev
Acceleration Technique,” IE E E Trans. Microwave Theory and Techniques, vol. 42,
no. 10, pp. 1998-2006, 1994.
[9]
F. Torres and B. Jecko, “Complete F D T D analysis o f microwave heating processing
in frequency-dependent and temperature-dependent media,” IE E E Trans.
Microwave Theory Tech., Vol. M T T -45, ppl08-117, 1997.
[10]
J. Guan and C. Su, “Resonant Frequencies and Field Distributions for the Shielded
206
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
U niaxially Anisotropic Dielectric Resonantor by the F D -S IC Method,” IE E E Trans.
Microwave Theory Tech., Vol. M TT-45, no. 10, pp l767-1777, 1997.
[11]
J. Fang, “Time Domain Finite Difference Computation f o r M axw ell’s Equations,”
Ph.D. thesis, University o f California at Berkeley, Berkeley, C A , 1989.
[12]
T. Deveze, L . Beaulie, and W. Tabbara, “An absorbing boundary condition for the
fourth order F D T D Scheme,” IE EE Antennas and Propagation Society
International Symposium, pp. 342-345, July, 1992.
[13]
S. K . Lele, “Compact finite difference schemes with spectral-like resolution,” J. o f
Computat. Phys., vol. 103, pp. 16-42, 1992.
[14]
J. L . Young, D . Gaitonde and J. S. Shang, ‘Towards the contruction of a fourth
order difference scheme for transient E M wave simulation: Staggered grid
approach,” IE E E Trans. Antennas and Propagation, vol. 45, no. 11, pp. 1573-1580,
1997.
[15]
M . H . Carpmter, D . Gottlieb, and S. Abarbanel, “The Stability of Numerical
Boundary Treatments for Compact High-Order Finite-Difference Schemes,” /. of
Computat. Phys., vol. 108, pp. 272-295, 1993.
[16]
A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain
Method, Artech House, Boston, M A , 1995.
[17]
A . Taflove, Advances in Computational Electrodynamics: The Finite-Difference
Time-Domain Method, Artech House, Boston, M A , 1998.
[18]
F. B. Hildebrand, Introduction to Numerical Analysis, M cG raw -H ill, New York,
N Y, 1974.
[19]
C. A . Balanis, Advanced Engineering Electromagnetics, John W iley & Sons, New
York, N Y , 1989.
[20]
F. J. Harris, “On the use o f windows for harmonic analysis with discrete Fourier
transform,” Proc. IEE E , vol. 66, pp.51-83, Jan. 1978.
[21]
J. A . Pereda, L . A . Vielva, and A . Prieto, “Computation o f resonant frequencies and
quality factors of open dielectric resonantors by a combination of the finitedifference time-domain (F D T D ) and Prony’s methods,” IE E E Microwave and
Guided Wave Letters, vol. 2, 1992, pp. 431-433.
[22]
R. Kumaresan, and D . W. Tufts, “Estimating the parameters of exponentially
damped sinusoids and pole-zero modeling in noise,” IE E E Trans. Acoustics,
Speech, and Signal Processing, vol. 30, 1982, pp. 1380-1419.
[23]
K . M . Chen, EE836 Class Notes, Michigan State University.
207
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[24]
J. D . Jackson, Classical Electrodynamics, John-W iley, N ew York, N Y, 1975.
[25]
D . B . Davidson, and R . W . Ziolkowski, “Body o f revolution finite-difference tim edomain modeling o f space time focusing by a three dimensional lens,” J. Optical
Society o f America, June 1993.
[26]
K . M . Chen, Notes f o r E E 836, Michigan State University.
[27]
D . B . Davidson, “Body-of-revolution finite-difference time-domain modeling of
space-time focusing by a three-dimensional lens,” J. Optical Society o f America,
vol. 11, no.4, pp. 1471-1490, 1994.
[28]
Y. Chen, R. M ittra, and P. Harms, “Finite-Difference Tim e-Dom ain Algorithm for
Solving M axw ell’s Equations in Rotationally Symmetric Geometies,” IE E E Trans.
Microwave Theory and Techniques, vol. 44, 1996, pp. 832-839.
[29]
J. G . Maloney and G . S. Smith, “The use o f surface impedance concepts in the
finite-difference time-domain method,” IE E E Trans. Antennas Propagat., vol. 40,
no. 1, pp. 38—48, 1992.
[30]
S. K ellali, B. Jecko, and A . Reineix, “Implementation o f a surface impedance
form alism at oblique incident in F D T D method,” IE E E Trans. Electromagnetic
Compatibility, vol. 35, no. 3, pp. 347-356, 1993.
[31]
K . S. Oh and J. E. Schutt-Aine, “An efficient implementation o f surface impedance
boundary conditions fo r the finite-difference tim e-dom ain method,” IE E E Trans.
Antennas and Propagation, vol. 43, no. 7, pp. 660-666, 1995.
[32]
D . M . Sullivan, “ Z transforms theory and F D T D method,” IE E E Trans. Antennas
and Propagation, vol. 44, n o .l, pp. 28-34, 1996.
[33]
J. H . Beggs, “A F D T D Surface Impedance Boundary Condition Using Z Transforms,” Applied Computational Electromagnetics Society Journal, vol. 13, no.
1, pp. 14-24, 1998.
[34]
J. G . Proakis and D . G . Manolakis, Digital Signal Processing, Macmillan, New
York, New York, 1992.
[35]
W. Y. Tan, “M odeling The Electromagnetic Field and The Plasma Excitation in A
M oderate Pressure M icrowave Cavity Plasma Source,”
Ph.D. Dissertation,
M ichigan State University, 1994.
[36]
D . M . Sullivan, “A Frequency-Dependent F D T D M ethod for Biological
Applications,” IE E E Trans. Antennas and Propagation, vol. 40, no. 3, pp. 532-539,
1992.
[37]
Om P. Gandhi, Ben-Qing Gao, and Jin-Yuan Chen, “A Frequency-Dependent
Finite-Difference Tim e-Dom ain Formulation for General Dispersive Media,” IE E E
208
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Trans. Microwave Theory and Techniques, vol. 41, no. 4, pp. 658-665,1993.
[38]
J. L . Young, A . Kittichartphayak, Y. M . Kwok, and D . Sullivan, “On the Dispersion
Errors Related to (F D )2T D Type Schemes”, IE E E Trans. Microwave Theory and
Techniques, vol. 43, no. 8, pp. 1902-1910, 1995.
209
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Документ
Категория
Без категории
Просмотров
0
Размер файла
5 914 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа