mmflwmés

..a. r ,t-r,
..

.flfimm...“ A

. .5 V ,. _
, .mmmmwwfi

3... hm.

gm,»
Us

.v.
u} .

“a
L.”

m“. x

. ‘ .. w

. . '9‘

U ".935 ‘ hit?
a. .4

4 a...

5.3.
HEM...» ..
at! Y z.
.11

. .

 

V!»- .r» _ .P
u ... (4.1%.. :5: .. 1:131.“

 

figaéfifimfii.

 

3:3
7

'2 LIBRARY

50‘ Michigan State
Unlversity

 

 

 

 

This is to certify that the
dissertation entitled
QUAN'l‘IFICATION OF INDUCED ELECTROMAGNETIC
FIELDS INSIDE MATERIAL SAMPLES PLACED INSIDE

AN ENERGIZED MICROWAVE CAVITY BY FINITE-

DIFFERENCE TIME-DOMAIN (FDTD) METHOD
presented by

Yao-Chiang Kan

has been accepted towards fulfillment
of the requirements for

Ph - D - degree in EFF:

zw

Major professor

Date 2‘3! 9 m

MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771

a... o' 4

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mm W.“

 

om
mm
mm

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 of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical and Computer Engineering

2000

QL'ANTIFIC
lVSlDE w
MICRWV

lie :mcsz‘ 2.”:

11:‘~i.

4|" I
\ H i I - - '
ambition 01 (it

it sax in

The”

goalofii

I LL ,. .
ie piesencc o

79*
“ALde‘ \\ h1\|n }

std ~
A) 1: based or

~ ‘4: tax};
Kim a;
1*~Et\

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 of the heating of a material sample in an energized microwave cavity
requires the understanding of the interaction of the electromagnetic fields with the
material sample in the cavity. The key factor for this understanding is to quantify the
distribution of the induced electromagnetic field inside the material sample placed inside
the cavity.

The goal of this research is to solve Maxwell’s equations in an electromagnetic cavity
in the presence of a material sample based on the finite-difference time-domain (FDTD)
method, which has been successfully applied to several areas in electromagnetics. This
study is based on Yee algorithm, a second-order FDTD scheme, and an improved fourth-
order FDTD scheme.

The numerical dispersion equations of Yee and other fourth-order FDTD schemes are
first discussed and the disadvantages of Yee scheme are discussed. An implicit staggered
fourth-order FDTD 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 of 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

.‘ -' “
{073.53.11.

Lt. m W"? 'L’“

‘Ut'i‘o HI: MIL-

 

03 32:05:11 cm.

0L 4" "
Um :fi...&. 13 d: ;

”3

”,0. ‘ r "‘h
in; ’i'f‘"_;‘\ v.

‘0' ., ALI.
A

    

:eratrd Slammer. '.

n -

“haw-4 0. - ; ‘00-

b‘.‘...r‘uu‘ué\\r-t1u l\ L l I I \
U

+fi :AfiaOc ‘D \\‘un |
I :h-lVl . \' I-wh

" l u
‘D ‘ .., , ”I
-' s);...u4.;.1. FDTL

hasn’t)“ 0. 9‘ .. ‘
Jadr‘bufcnk‘iinK‘n $2;

consistent.

In applying the FDTD method to cylindrical coordinates, the singularities at the center
of cylindrical coordinates become the major problem. The body of revolution (BOR)
FDTD method is applied to solve the filed distributions in the cylindrical cavities loaded
with samples with symmetries. The treatment in BOR FDTD method is quite straight
forward. Moreover, the BOR FDTD method is a 2.5D FDTD method which is much more
computationally efficient than a 3D FDTD method. For a sample with any arbitrary shape,
the general 3D cylindrical FDTD method is needed to do the calculation. The traditional
3D cylindrical FDTD method encounters the difficulties of the mode-dependent source
implementation and the treatment of singularities. In this study, a general fourth-order
FDTD method is proposed to overcome the problems encountered in a traditional 3D

cylindrical FDTD method.

To my parents and family

iv

 

 

 

. .
‘r, ’1‘4! 93% -.
Faun. .U.t...- ‘

ll. Chen. for it’s-x

475- n, - :-
m-hi‘}"\h1~\. 2'
5

mac been pm.“
~39 ‘ er -
“scanner and an:

1 mid Lie to

I‘l- '
'2’1 "Hut Y
”dB-A- ROLIV‘CN g

ll} srecral t'n xix
missions. Sheet.-

31'} Similatnon pm.

IWOUld .1150 l;

ACKNOWLEDGMENTS

First and foremost, I would like to express my deep gratitude to my advisor, Dr. Kun-
Mu 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 of my committee, Dr. Dennis Nyquist, Dr.
Edward Rothwell, and Dr. Byron Drachman for many helpful comments on this research.
My 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 of Division of Engineering
Computing Services, for appointing me as a graduate assistant from 1994 to 1998. My
special thanks to Dr. Guilin Cui, owner of Vertex Computer in Lansing, for providing me a
job for the past year. Without those financial supports, I would not be able to study my
Ph.D. at Michigan State University.

My very special thanks to my wife, Ya-Ling Peng, who have been taking care of our
children and me well so I am able to focus on my research and work.

Finally, I would like to thank my parents for providing me with the opportunity to
complete this study. Their love and support accompanied me through the years of this

work.

 

 

 

EST OF TABLES
MST OF HGL'RILS
GHPTER1
BTRODU
(INTER:
SOLVIXG}
COORDIN.
2.1 Free
Till}
2.
2.1..
2'3 Tic
13 Tet.
2.3.
2.3
2.3
2.3.
3.3.
2.4 Ext
2.4
7 2.4
MS Pm
3.5.
2.5.
ML.
26
As
3.6.
2.7 2'6

TABLE OF CONTENTS

LIST OF TABLES ...................................................... ix
LIST OF FIGURES ...................................................... x
CHAPTER 1
INTRODUCTION ................................................. 1
CHAPTER 2
SOLVING MAXWELL’S EQUATIONS BY FDTD IN RECTANGULAR
COORDINATE ................................................... 5
2.1 Frequency-Dependent FDTD Formulations with Second-Order Leapfrog
Time-Stepping of Maxwell’s Equations .......................... 7
2.1.1 The Scalar Equations of Maxwell’s Equations ............... 7
2.1.2 The Finite Difference Equations .......................... 9
2.2 The Yee’s Algorithm ........................................ 17
2.2.1 Stability Condition .................................... 19
2.2.2 Numerical Dispersion ................................. 19
2.3 The Ty(2,4) (FD)2TD Algorithm .............................. 26
2.3.1 The Fourth-Order Space Derivatives ...................... 26
2.3.2 Dispersion Analysis for Explicit Staggered Scheme .......... 27
2.3.3 Dispersion Analysis for Implicit Staggered Scheme .......... 29
2.3.4 Calculation of Derivative of E_. in Ty(2,4) ................. 33
2.3.5 Calculation of Derivative of H in Ty(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 Method ............................................ 44
2.5.1 Theory ............................................. 44
2.5.2 Estimation of Resonant Frequency and Quality Factor by Prony
Method ................................................... 47
2.6 A Single Empty Cavity with PEC Walls ......................... 48
2.6. 1 Configuration ........................................ 48
2.6.2 Numerical Results of Field Distributions .................. 48
2.7 A Lossless Loaded Cavity with PEC Walls ....................... 65
2.7.1 Quasi-cubic Case ..................................... 65

vi

 

4.3

4.4

 

 

OJIJJp'J"
I I l

 

' .IJ

,) ’J .0)”
lJlJ'JlJ')!

2.7.2 Thin Square Plate Case ................................ 68

2.7.3 Narrow Strip Case .................................... 72
2.8 A Lossy Dielectric Loaded Cavity with PEC Walls ................ 76
2.8.1 Configurations ....................................... 76
2.8.2 Numerical Results and Discussions ....................... 77
CHAPTER 3
SOLVING MAXWELL’S EQUATIONS BY BODY OF REVOLUTION
FDTD .......................................................... 88
3.1 The BOR Formulation of Maxwell’s Equations[28] ................ 89
3.1.1 Mode Selection in BOR Algorithm ....................... 93
3.1.2 The BOR-FDTD Formulation[28] ........................ 95
3.1.3 Singularity in BOR-FDTD Formulation at p = 0 ............ 99
3.2 Surface Impedance Boundary Condition ........................ 100
3.2.1 Planar Surface Impedance Boundary Condition ............ 100
3.2.2 FDTD Implementation of Planar Surface Impedance ........ 103
3.2.3 Time-Domain Approximation by Prony’ Method ........... 104
3.2.4 Frequency Domain Approximation ...................... 105
3.2.5 Z Transform and Digital Filters Approach ................ 106
3.2.6 Fields Calculation on the Cavity Wall .................... 108
3.3 An Empty Cylindrical Cavity ................................ 109
3.4 A Loaded Cylindrical Cavity ................................. 118
3.4.1 Small cylindrical sample for TM012 mode ................ 118
3.4.2 Thin rod case for TM012 mode ......................... 119
3.4.3 Thin disk case for TM012 mode ........................ 119
3.4.4 Small cylindrical sample for TEl 11 mode ................ 121
CHAPTER 4
SOLVING MAXWELL’S EQUATIONS BY FDTD IN CYLINDRICAL
COORDINATES ................................................ 136
4.1 Three Dimensional FDTD Representation of Maxwell’s Equations in
Cylindrical Coordinates ..................................... 137
4.2 The Second-order Cylindrical FDTD Scheme[35] ................ 143
4.2.1 The Second-order Cylindrical FDTD Equations ............ 143
4.2.2 FDTD Calculations at p = O ............................ 145
4.2.3 Source Implementation ............................... 146
4.2.4 Numerical Results and Discussions ...................... 148
4.3 Ty(2,4) Cylindrical FDTD Scheme ............................ 165
4.3.1 Derivatives of Fields in p direction ..................... 165
4.3.2 Derivatives of Fields in (p direction ...................... 172
4.3.3 Derivatives of Fields in z direction ...................... 175
4.3.4 FDTD Calculations at p = 0 ........... ‘ ................. 176
4.4 Time Domain Finite Difference Equations of Constitutive
Relations ................................................ 179

vii

 

 

GAMERS I
COXCU '

WENDIX A
DER VIA?
APPROX".

WEXDIX B
DERVIA"
STAGGE;

1P?E.\DIX C
DERVi-x

DER\TA

BBUOGRAPI

CHAPTER 5

CONCLUSIONS ................................................ 1 82
APPENDIX A

DERVIATION OF FOURTH-ORDER FINITE DIFFERENCE

APPROXIMATION ............................................. l 83
APPENDIX B

DERVIATION OF NUNTERICAL DISPERSION RELATOIN FOR IMPLICIT

STAGGERED SCHEME .......................................... 189
APPENDIX C .

DERVIATION OF EQUATION (2.107) ............................. 201
APPENDIX D

DERVIATION OF EQUATIONS (3.7) AND (3.8) ..................... 204
BIBLIOGRAPHY ..................................................... 206

viii

 

I ‘D
a. a.
.5.“ .JL.
.i .9
91. Ti

1\

T?"

.28; ‘ ‘
Ila-l: .-.‘

?
1

BO?

.fila
A.

m

Table 2.1

Table 2.2

Table 2.3

Table 2.4

Table 3.1

Table 3.2

LIST OF TABLES

EW 5/ E n 5 ratio at different times ............................... 67
Permittivities mapping at 're = lns, (r) = 21t(2.45e9) ............. 77
The loss power and stored energy ............................... 78
Properities of four different lossy cases .......................... 79
BOR reprenstation of Maxwell’s equations ....................... 92
Rational Approximation Results ............................... 107

ix

Figure2.l

Figure2.2

Figure2.3

Figure2.4

Figure2.5

Figure2.6

Figure2.7

Figure2.8
Figure2.9

Figure2. 10

Figure2. l 1

Figure2. 12

Figure2.l3

LIST OF FIGURES

FDTD space lattice .......................................... 11
Variation of the numerical phase velocity with CFL values at R = 10

and t1) = n/ 4. ............................................. 22
Variation of the numerical phase velocity with R at a = 0.9 and

(I) = n/ 4 . 23
Variation of the numerical phase velocity with d) at a = 0.9 and

R = 10. 24
Variation of the numerical phase velocity with HR at a = 0.9 and

(b = N4. ................................................... 25
Comparison of numerical phase velocity between explicit 4th-order scheme
and Yee scheme at R = 5, a = 0.24, and q) = 1t/4. ................... 31
Comparison of numerical phase velocity between explicit 4th-order scheme
and Yee scheme at a = 0.24, and o = 1r/4 .......................... 32
The lattices of By and aEy/az along the z axis at fix it and y. ....... 33
The lattices of Hz and 3 Hz /ay along the y axis at fix x and z. ....... 35

The configuration of the empty rectangular cavity with PEC boundary. The
excitation probe is located on one of the six faces. ................. 50

The x dependence of Ex at 2: 0.025974 and t= 0.26658 ms for TEOll mode.
.......................................................... 51

The x dependence of Ex at y= 0.034632 and t= 0.26658 ms for TEOll mode.
.......................................................... 52

The y dependence of Ex at x: 0.025974 and t: 0.26658 ms for TEOll mode.
.......................................................... 53

Egrtlli

"av-93‘ 1‘

I;UA\-.
b

f.:::2.lfi

Egrsll’

I

Fzgzzl. lS

F33:63.21

Fifi-=1 11

\
.‘ a...

F.

:‘I ‘
-.-re._z3
13.
59433.24

F3:22.15

I:
I? -. i
wait

F

:‘N‘v 1
31,3 h
h..’

‘13H.

'tw-E] 38

 

 

lhln
is ii

Figure2.l4

Fi gure2. l 5

Figure2. 16

Figure2. l7

Figure2. 1 8

Figure2.19

Figure2.20

Figure2.2 l

Figure2.22

Figure2.23
Figure2.24

Figure2.25

Figure2.26

Figure2.27

Figure2.28

The z dependence of Ex at X: 0.034632 and t= 0.2665 8 ms for TE011 mode

.......................................................... 54
The Ex variation along the yz plane at x: 0.034632 and t: 0.26658 ms for
TEOll mode. .............................................. 55

The B field on the xz plane at y: 0.034632 and t: 0.26658 ms for TEOll
mode. 56

The x dependence of Ex at 2: 0.03192 and t: 0.31019 ms for TMl 11 mode.
.......................................................... 57

The x dependence of Ey at 2: 0.04256 and t: 0.31019 ms for TM111 mode.
.......................................................... 58

The x dependence of E2 at y= 0.04256 and t= 0.31019 ms for TMl ll mode.
.......................................................... 59

The E field on the xy plane at 2: 0.04265 and t= 0.31019 ms for TMl ll
mode ..................................................... 60

The E field on the xz plane at y= 0.04265 and t: 0.31019 ms for TMl 11
mode. .................................................... 61

The E field on the yz plane at X: 0.04265 and t= 0.31019 ms for WI 11
mode ..................................................... 62

The time variation of Ex at x=0.04329, y=0.04329, and z=0.034632. . . 63
The frequency response of Figure 2.23 ........................... 64

Dimensions of the rectangular cavity and the loaded material sample. The
center of the material sample is consistent with the center of the cavity.
.......................................................... 66

The variation of electric fields in the y direction at t = 0.2535 us. . . . . 69

The variation of By in the x direction at t = 024289115 . The line with star
symbol is the calculated values and the line with circle symbol is the fitted
values for the empty cavity. ................................... 70

Variations of E w 5/ E n S in the x-directions. Each curve represents this ratio
as a function of x for different locations of z. The relative permittivity of the
thin square plate material sample is 8, = 2.5 . The solid line with symbols
is the ratios at t = 0.24289us and the dash line with symbols is those at

t = 0.24301us . ........................................... 71

xi

F"":‘ “
“5...“33

t‘s .. ‘
“5.16134

532'31 ‘8

tu' 59.j-

 

 

 

s' ..
4m}.

The \

“1an

Figure2.29

Figure2.30

Figure2.3 1

Figure2.32

Figure2.33

Figure2.34

Figure2.35

Figure2.36

Figure2.37

Figure2.38
Figure3. 1
Figure3.2

Figure3.3

Figure3.4

Figure3.5

The variation of electric fields in the y directions at x=0.34839m,
z=0.057093m, and t=0.252ms .................................. 74

The ratios of Em/ E n s in y directions. .......................... 75

The flow chart of calculating the field distribution at steady state for high Q
cavities. .................................................. 80

Stored energy for the cavity with material sample of e',((n) = 2.5 and
8",(tu) = 0.1 . Note that one time step is equal to 4.25442 ps and this is a
downsampling plot. ......................................... 81

Instaneous dissipated power for the cavity with material sample of
e',((n) = 2.5 and 8",(tn) = 0.1 . Note that one time step is equal to
4.25442 ps and this is a downsampling plot. ...................... 82

The plot of the fitting and calculated data for stored energy 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.45x10’lo and the time average
stored energy is 8.905618. .................................... 83

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 with cross is the

calculated power data. The average input power is 2.45x10'10. ....... 84

The plot of the difference between fitting and dissipated power data. The

one time step is equal to 4.25442 ps. ............................ 85
Field distributions of Ey along x axis at 0.14868 ms. 8',(tr)) is equal to 2.5
and 8",(tn) is 0.1. .......................................... 86
The ratio of calculated and fitting Ey along x axis at 0.14868 ms. ..... 87
The field locations for BOR FDTD in time and space. .............. 98
FDTD configuration for cavities with FEC wall .................. 102

Coordinates for the incident and reflected plane waves upon a lossy
surface .................................................. 103

The physical configurtion of a cylindrical cavity loaded with a material
sample. .................................................. 1 l 1

The variation of Ep of TM012 along the p and 2 directions in a PEC empty
cylindrical cavity ........................................... 112

xii

 

ti

Figure3.6

Figure3.7

Figure3.8

Figure3.9

Figure3.10

ngure3.l 1

Figure3.12

Figure3. l3

Figure3. l4

Figure3.15

Figure3. 16

Figure3.17

Figure3.18

Figure3.19

Fi gure3.20

The variation of E2 of TM012 mode along the p and 2 directions in a PEC

empty cylindrical cavity. .................................... 113
Plot of B field of TM012 mode on the p-z plane in an empty PEC cavity.
......................................................... 1 14
The variation of Ep of TE111 mode along the p and 2 directions in a PEC
empty cylindrical cavity. .................................... 115
The variation of E,» of TE] 11 mode along the p and 2 directions in a PEC
empty cylindrical cavity. .................................... 116
Plot of E field of TE111 mode on the p-z plane in an empty PEC cavity.
......................................................... 1 17
Plot of 2Ep of TE] 11 mode versus time in an empty FEC cavity with
O=102pS/m and 0': 102 S/m withm = 1. ................. 120
Plot of E2 of TE111 along r direction in an empty FEC cavity with
conductivity 104 (S/m). ..................................... 122

Plot of E field of TEl 1 1 mode on the p- -z plane 1n an empty FEC cavity with
conductivity 104 (S/m). ..................................... 123

Plot of ED of TM012 mode in a PEC cavity loaded with a small cylindrical
material sample. ........................................... 124

Plot of E2 of TM012 mode in a PEC cavity loaded with a small cylindrical
material sample. ........................................... 125

Plot of the ratio of of TM012 mode along the z direction in the PEC cavity
loaded with a small cylindrical material sample. .................. 126

Plot of ED of TM012 mode in a PEC cavity loaded with a thin rod material
sample. .................................................. 127

Plot of E2 of TM012 mode in a PEC cavity loaded with a thin rod material
sample. .................................................. 128

The ratio of E Z/ E 2 of TM012 mode in a PEC cavity loaded with a thin rod
material sample. ........................................... 129

Plot of Ep of TM012 mode in a PEC cavity loaded with a thin disk material
sample. .................................................. 130

xiii

5233.13

531244

 

 

Figure3.21

Figure3.22

Figure3.23

Figure3.24

Figure3.25

Figure4. 1

Figure4.2

Figure4.3

Figure4.4

Figure4.5

Figure4.6

Plot of E2 of TM012 mode in a PEC cavity loaded with a thin disk material
sample. ............... ' ................................... 131

The ratio of Ez/E: of TM012 mode in a PEC cavity loaded with a thin disk
material sample. ........................................... 132

Plot of ED of TEl 11 mode in a PEC cavity loaded with a small cylindrical
material sample. ........................................... 133

Plot of E4, of TEl 11 mode in a PEC cavity loaded with a small cylindrical
material sample. ........................................... 134

Plot of E2 of TE111 mode in a PEC cavity loaded with a small cylindrical
material sample. ........................................... 135

FDTD lattice for cylindrical coordinates. ....................... 141

The diagram of order of FDTD calculations along time axis. The meanings
of those steps are listed below.
(1) Using (4.9) to (4.11)

(2) Desired time stepping scheme for D.

(3) Constitutive relation of D and E which depends on material models.

(4) Using (4.12) to (4.14). At this point, the boundary conditions are
applied.

(5) Desired time stepping scheme for I}.

(6) Constitutive relation of I; and H which depends on material models.
......................................................... 142

TE011, TEl 1 1, TM011, and TM012 modes excitation techniques in a
cylindrical cavity[35]. ...................................... 147

The plots of total stored energy of TM012 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 BH source implementation.
......................................................... 149

The variation of Ep of TM012 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. ........................................... 150

The variation of ED of TM012 mode along the z direction in a PEC empty
cylindrical cavity. The upper figure is calculated by using the traditional

xiv

q
INII”4
F:;Ua5 - .
5

PM,“
I‘mmb .
V

Err-$49
A l...‘5 {
5

Him. 10

k
I

9'13."
“3154. 13

R,
“st-’64“

R...
‘ $13115

 

 

 

Figure4.7

Figure4.8

Figure4.9

Figure4.10

Figure4.] l

Figure4. 12

Figure4. 13

Figure4. 14

Figure4. 15

source implementation and the lower one by using the BH source
implementation. .......................... ' ................ 151

The variation of E2 of TM012 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. ........................................... 152

The variation of E2 of TM012 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 BH source
implementation. ........................................... 153

The plots of total stored energy of TEl 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 BH source
implementation. ........................................... 154

The variation of ED of TEl ll 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. ........................................... 155

The variation of ED of TE] 11 mode along the (1) 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. ........................................... 156

The variation of ED of TE] 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 BH source
implementation. ........................................... 157

The variation of E.» of 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

The variation of E4, of TEl ll mode along the a 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. ........................................... 159

The variation of E4, of TEl 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 BH source

XV

Fzga'sl. 16

F‘ 4 "
’3
-;dr\ .1
v

Egret 18

15384.19

55’36420

EMI,
“5442]

pm»
“ERG

 

C3

' I

31
l'“" I
rill

PI

l0;

 

Figure4. 16

Figure4. 17

Figure4.18

Figure4. 19

Figure4.20

Figure4.21

Fi gureC 1

implementation. ........................................... 160

The variation of Ep of TM012 along the p and 2 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 BH source
implementation. ........................................... 162

The variation of E2 of TM012 along the p and 2 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 BH source

implementation. ............................................ 163

Plot of the ratio of of TM012 mode along the z direction in the PEC cavity
loaded with a small cylindrical material sample. .................. 164

The locations of H field and its corresponding derivatives in
p direction. ............................................... 166

The locations of B field and its corresponding derivatives in
p direction. ............................................... 171

The FDTD lattice along the (1) direction. 174

Approxmiation of intergrant in (CI) ........................... 202

xvi

CHAPTER 1

INTRODUCTION

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 of 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 (FDTD) method is employed in this dissertation to
quantify the electromagnetic field inside an EM cavity loaded with material samples.

The finite-difference time—domain (FDTD) method has been used in computational
fluid dynamics (CFD) [1] for a long time and yields very accurate results for CFD
problems. In 1966, Kane Yee originated a set of finite-difference equations for the time-
dependent Maxwell’s curl equations system for the lossless materials case [2]. The FDTD
method was not popular in (Electromagnetic) EM research area until late 1980 and
becomes a very popular method in BM area between 1993 and 1997. Regarding to the
FDTD method applied to eigenvalue problems in EM research, Choi and Hoefer [4]

published the first FDTD simulation of waveguide/cavity structures in 1986. There are

citations 11‘»

' l
mzfed 5mm)“

515;: the cunt:
. LA I , -g F
meta-d 1S .1 3.1...
.zzrg FDTD s."

1 -- -
iglrnm used if

 

 

 

 

 

 

Hewitt. the d2»:

1”"
1.... be discussed:

Most papers H

33:01:

mzi'dons u

ACCUQC) and ma".

several papers, [5]-[10], that utilized the FDTD method to investigate the fields and power
distributions in a loaded EM cavity. The paper by Torres and Jecko [9] provides a very

detailed analysis of microwave heating by combining the Maxwell’s equations and the
heat transfer equations. The FDTD method used in this paper is called (FD)2TD method

since the constitutive parameters are assumed to be frequency-dependent. This (FD)2TD
method is a standard method to investigate EM interaction with a lossy material sample
using FDTD scheme. Moreover, the combined electromagnetic and thermal FDTD
algorithm used in this paper provides a basic framework for the FDTD calculation.
However, the dissipated power model in this paper is not clear and the improved model
will be discussed in this dissertation.

Most papers in the FDTD literature are based on second order spatial and temporal
approximations which are originated from Yee algorithm. Due to the requirement of
accuracy and performance, several higher order FDTD methods have been proposed [11]-
[15]. Among these higher order FDTD methods, Ty(2,4) FDTD method, which uses the
implicit 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 (FD)2TD method, the Ty(2,4) (FD)2TD method becomes the
essential FDTD technique in the investigation for the eigenvalue problems in EM research

area.

In this dissertation, the Ty(2,4) (FD)2TD method in rectangular and cylindrical
coordinates are studied. For a cylindrical cavity with azimuthal symmetry, the body of
revolution (BOR) scheme is employed to facilitate the FDTD calculation. In this case, the

BOR FDTD method can give very accuracy results with excellent performance; hence, the

Mi} mm 7
F

 

In chapter L
induced. Aix
mentor. :~
diurnal d
atiL‘iidzspcn; '
11 deals. A {11".}:

FDTD calcdiath '-

faion (Q) of c.-

 

j,.. .. .. '. .
th$sߣdl \OEGJSI:

n

and PEC “.1333. a..
a: of this chapter
mam.

In chapter 3. lhi

is $31le. Mil-dc

Presemcd The {red
1*: ~ .-
ha .hdpter. the

FEC Wall is Tepid;

PETD

{01711111411011

Erma}
n apPTOXImd

a!
11'

W] 081111 R

Premfd at 1h
0
m
111
TD method for C

Ty(2,4) FDTD method is not used.

In chapter 2, a set of general finite difference equations for Maxwell equations is
introduced. Also, the loaded material is modeled by Debye equations and its FDTD
representation is presented. Then a brief introduction of Yee algorithm, stability condition

and numerical dispersion is presented. The fourth order spatial derivatives are presented

and the dispersion analysis is studied. After that, the Ty(2,4) (FD)2TD method is presented
in details. A time-domain power analysis based on Poynting’s theorem is derived. In this
FDTD calculation for cavities, the Prony’s method is employed to estimate the quality
factors (Q) of cavities. The numerical results of a single empty cavity with perfect
electrical conductive(PEC) walls, a cavity loaded cavity with a lossless material sample
and PEC walls, and a lossy dielectric loaded cavity with PEC walls are presented at the
end of this chapter. The numerical results are shown to be consistent with the theoretical
estimation.

In chapter 3, the body of revolution (BOR) FDTD formulation of Maxwell’s equations
is derived. Mode selection and source implementation in BOR FDTD algorithm are
presented. The treatment for the singularity in BOR FDTD formulation is also presented.
In this chapter, the cavities with finite electrical conductive (FEC) walls is studied and the
FEC wall is replaced by a planar surface impedance boundary condition (PSIBC). The
FDTD formulation of PSIBC is achieved by three different approachs and the frequency
domain approximation is chosen due to its simplicity. Numerical results of an empty
cylindrical cavity with PEC and FEC walls and a loaded cylindrical cavity with PEC walls
are presented at the end of this chapter. Consistent results are obtained by using this BOR

FDTD method for cavities with azimuthal symmetry.

ill-[$7385 11R“

 

generaized. 8\
team: for 1h-
Lieu maids :1

V
’1'-
H.

'1 113.11 1.007117

 

Some dt’i'l‘i'dih

In chapter 4, a general 3D FDTD method in cylindrical coordinates is considered and
studied. The disadvantages of conventional second order FDTD method in cylindrical
coordinates are presented and its improvements are also proposed. The treatment for the
singularity in conventional cylindrical FDTD is mode-dependent and difficult to be
generalized. By the introduction of Ty(2,4) FDTD in cylindrical coordinates, a general

treatment for the singularity in cylindrical FDTD can be obtained. With the Debye or

Lorent models for loaded material samples, a general Ty(2,4) (FD)% method in
cylindrical coordinates is obtained.

Some derivations that are useful in this dissertation are provided in Appendices.

SOLVIM

IN

 

 

 

 

In this Ella?"

\\
.

: :l: O‘- '

fir. ‘1‘
aides. The “M"

llIii-\

seconderder 21cc:

CHAPTER 2
SOLVING MAXWELL’S EQUATIONS BY FDTD

IN RECTANGULAR COORDINATE

In this chapter we considers the frequency-dependent implicit, fourth-order FDTD
formulation in the rectangular coordinate system and its application to rectangular
cavities. The finite difference approximation to time-stepping in those formulations is of
second-order accuracy. The second-order accuracy in time-stepping combined with
implicit staggered fourth-order accuracy in space, denoted by Ty(2,4)[17], is the focus of
this chapter. Another scheme, Ty(4,4), that uses fourth-order in time-stepping with
implicit fourth-order in space-stepping is not discussed because of 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 of Ty(2,4) and Ty(4,4) are made comparable by choosing an

appropriate time step[14][17] although the stability of the FDTD formulation of
Ty(4,4) may be improved[15].
(3) The Ty(2,4) is nondissipative, while Ty(4,4) introduces a slight dissipation.

In FDTD calculations for cavities, a many time steps are usually required if the quality

  
   
  
  

razor. Q. 01' Li"
chaser meshes
11.24.11 131101 :1 ;

Tye freque'

11:32:22 of 31.11
‘ . b
tamer. lS '

The or.

no“
L; 1 Id}
§

sector. 2.2. The s

; ’ B- I I
‘ A“ c ' U r!
11.. s rigors".

In section 2.?

T’DLAO ~
.L:.t.~‘t."d€r 5:34"

\I'

515 8.1316 disperxz

5‘01? : ‘ ‘
Wh-‘Mlle 4 SUV}

.1... ier
mesh reduction. '
k5“.

dies are d: ix.

Schemes requ; n n

1‘;

rugged in Sect

33.in- .
. UCdimg

factor, Q, of 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 for a long time integration.

The frequency-dependent FDTD formulation with second-order leapfrog time-
stepping of Maxwell’s equations is discussed in section 2.1. The one-relaxation Debye
equation is used to account for the frequency-dependent properties.

The original Yee’s algorithm is derived from this general scheme for validation in
section 2.2. The stability condition of Yee’s algorithm is also presented and the dispersion

in Yee’s algorithm is then explained.

In section 2.3, the Ty(2,4) (FD)2TD 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 FDTD scheme and

require a smaller time step for stability, the computer running time is complemented by the

mesh reduction. The numerical dispersions of explicit fourth-order and Ty(2,4) (FDYTD
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 of 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 with lossless dielectric material samples, the Q values of the
cavities are infinite so only transient-state solution are obtained to validate the program.
The source used for this case is a Blackman-Harris type that has very low sidelobes. For

loaded cavities with lossy dielectric material sample, a single frequency sinusoidal source

6

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 FDTD 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 Maxwell’s Equations

2.1.1 The Scalar Equations of Maxwell’s Equations
In differential form, Maxwell’s equations in a dielectric dispersive medium can be

written as

- at?
VXE - 4th-
_\ (2.1)
VxH = 9.9+[6]E+j3

where [o] and [it], are electric conductivity, magnetic permeability which are non-

A
dispersive in tensor form, respectively. J s 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

-1

-01x 0 0
[a] = 0 cry 0 (2.2)
O O orz

 

 

where or represents the magnetic permeability or the electric conductivity.

For frequency-dependent dielectric material, the one-relaxation Debye equation is

   

‘ v a " H§\
were 3% 5““

 

aides when t! 1‘

. A ‘
Between D and 1

H“

De“ '3? equation:

s (to) = 8' (m)—je" (to) = 8' +352 (2 3)
' ' ’ ’°° l+jmTe '

where the subscript, r, denotes the word “relative”. Moreover, 8', and 8", are the real and

imaginary parts of the relative permittivity. The 18 is a new relaxation time constant .

related to the original relaxation time constant 1- by [19]

I
T _ Tend-2

e 8—,":1-5 (2.4)

where 8'” and 8',” are the real part of the complex permittivity at zero frequency and at a
very large frequency, respectively. Note that e,(tu) is equal to e',((u) and 8',” can be any

values when 1: e is equal to zero, the frequency-independent case. The constitutive relation

between D and E is

em 0 0
D(w) = so 0 5,), 0 EU») (2.5)
O 0 am
where an, ery, and e,Z are the relative perrnittivities in different directions and satisfy

Debye equation, (2.3).

The scalar equations of (2.1) in time domain can be written as

BHx 3Ey 3Ez
u — = —-—
x a: 32 By
a”? — EJE‘ (2.6)

Ply—E— - 3x 32
BHZ _ aEx 3E),

“xv—arm-

 

- '\
The eqaaiidn t__

 

:31" the ln\ wa

11:27.0“

an_aHZ 3Hy E J
at - W—fi-Gx x— ii

any arr, an,
an _ 3H 8H

77717—87 H ‘2

 

The equation (2.5) can be written as

(1 + jwtex)Dx = OerxsEx + Eoe'rxootexwax
(1 + jw‘tey)Dy = eoe'mEy + eoe'moreyijy, (2.8)
(l + jwtez)Dz = Oerstz + 808'rz°‘="cezj(‘0Ez

then the inverse Fourier transform is applied; hence, the scalar equations are obtained as

follow

an , aEx

Dx + 1"er—ET = OErxsEx + 806 rxootex—aT
D +t a—D—y = E +8 8' 1.’ £352 (29)

y ey a: Erys y 0 ry°° ey a;

an, . at:z

l)z + T€Z_5t_ = Oerstz + 808 rzootezw

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 of I? are evaluated at integer time step

but H and the temporal derivative of D and E are evaluated at half integer time step. For

_\ n + 1/2
, and H I are considered to be the field distributions.

4" A"

the time step n, the E , D

 

 

For the spatial discretion, the second-order Yee scheme or higher-order scheme can be

used. The finite difference approximation to the first derivative of those quantities is

11220151111} 5.. :

mMHEU

 

I .. 1.
10.42533 111 111C 3,.
mam *er .;

“its;
5

11:3; so are the

 

‘l‘l
00"”?! 3
AI“ 7
111.1
5.1) M1 ~
"Hal
6.0”” w
11“] "1

denoted by 80,113 where or and B represent one of x, y, z, and t parameters and A denotes

one of H, E, D fields. The spatial locations of E and H are plotted in Figure 2.1 and D ,

locates at the same location as E does. Note that E x, aEx/at , and an/at are evaluated

at half integer spatial step along the x axis, but at integer spatial step along the other two

axes; so are the Ey, Ez, and their corresponding time derivatives. The H x, H y, H z are

evaluated at integer spatial step along the x, y, and z axes, respectively, but at half integer
A
space step along the other two axes; so is the time derivatives of H . The spatial derivatives
_\ A _\
of E are evaluated at the same locations as H and that of H are evaluated at the same

locations as E . Hence, (2.6), (2.7), and (2.9) can be rewritten as

5,114" = —(SZE I" -5 EzI" )
i,j+l/2,k+l/2 1.1 y i,j+l/2,k+l/2 y i,j+l/2,k+l/2

1
1H ln = — szln _82Ex|n (2-10)
yi+1/2,j,k+1/2 11y i+1/2,j,k+l/2 i+1/2,j,k+1/2
1
5,114" = — a ,1" 43; I"
i+l/2,j+l/2,k Y i+1/2,j+1/2,k Y i+1/2,j+l/2,k
Z
1/2 n+l/2 n+l/2 n+l/2 n+l/2
so ’” = (a H —6 H )— —

‘ x|i+l/2,j,k y Zli+1/2,j,k Z y|i+l/2,j,k ’ ‘|i+1/2,j,k "li+1/2,j,k
n+l/2 n+l/2 n+l/2 n+l/2 n+1/2
OtDyl, . - (52 xl. . — x 2'. . )_ y yl. . _ yl. . £2.11)
r,]+l/2,k r,}+l/2,k r,1+1/2,k r,1+l/2,k r,1+l/2,

n+l/2 n+l/2

n+l/2
7- zli,j,k+1/2 Zli,j,k+ 1/2

x yli,j,k+l/2 y xli,j,k+l/2)

8tDZIiigi/z = (8 H n+l/2
t,1,k+1/2
n+l/2 n+1/2
Ilia/2,131: 9‘ ‘ "li+1/2,j,k
n+1/2 1: SD n+l/2
yli,j+1/2,k 3” yli,j+l/2,k
n+1/2 n+l/2 _
zli,j,k+1/2 ez ' zli,j,k+1/2 —

= 801-:

, n+l/2

n+l/2

. n+1/2

e e
0 '35 zli,j,k+1/2

10

”‘5 x|i+ 1/2,j,k

' E
'2" yli,j+l/2,k

n+1/2
rxootexst xl
i+1/2,j,k

n+l/2
0°12 8E
9’ 0‘ yli,j+l/2,k

n+1/2
0°13 6E
’2 ‘Z ‘ zli,j,k+l/2

+ £08
+ 808' (2.12)

I
+ £08

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. .
. _
. . D
_ z _ 1i
- _
9 t a.
. _ .tuw
1 11 111... 1 _
a O 1" Ix .
/
. t. x _
. H a.
_ 1 F
I _ ..
_ _ ..
_ . . . v.
_ _ _ z e
_ s. _ _ .1 E
_ . _ _ z
0 _ A .1 .
11111 .--..--4.F.--i---. . _
I I I _
Hy // E z ’y .
z _
z
I II ..x
IIIII willlnl F'llllllllul .‘E
X I O ,,
E z
z z .
z z I
z z I
z I; I )
z
IIIIII gut—Illlllllllé k
ii.
#2 a

 

 

Figure 2.1 FDTD space lattice

ll

 

    
 

332mm: ’

e
k .

m 7315!". 35 IL'

It“:
H}.
t‘l:

H:""‘3

‘r(s\ ‘

a-

Apply the second-order leapfrog scheme to time stepping and use the average value to

approximate the D, E, J at n + 1/ 2 time steps, equations (2.10), (2.11), and (2.12) can be

rewritten as follow.

n+1/2
" i,j+l/2,k+1/2

n+ 1/2
yi+1/2,j,k+l/2

+1/2
H n
Z i+ l/2,j+ l/2,k

 

 

n—l/2
xli,j+1/2,k+1/2

At
+—

x

n—1/2
y|i+ l/2,j,k+1/2

At
+—
y
n—l/2
z|i+1/2,j+ l/2,k
At
+—

Z

n+1 OxAt n+1 _ n
x i+1/2,j,k x|i+l/2,j,k xli+1/2,j,k
0 At +1/2
— ‘ xl +At(6 HZI" —
2 i+l/2,j,k y i+l/2,j,k

At 1
-——-(Jx|"+ +
2 i+1/2,j,k

oyAt

n+1
yli,j+l/2,k

Gym

2

2 yli,j+1/2,k

At +1
V i,j+ l/2,k

n )
x|i+ 1/2, j,k

n+1 n

ylt‘,j+1/2,k = yli,j+l/2,k

+ At(6sz|"+ 1/2

n )
yli,j+1/2,k

12

6 " —5 E " J

p.( Z yli,j+1/2,k+1/2 y zli,j+l/2,k+1/2
n n

( x 2'. - —82Ex|. . )

p. .+1/2,1,k+1/2 1+1/2,1,k+1/2

n n
—5 E )
p.( y x|i+l/2,j+l/2,k x y|i+1/2,j+1/2,k

+1/2
SZH l" )
Y i+1/2,j,k

n+l/2 )
i,j+l/2,k x zIt,j+1/2,k

(2.13)

(2.14)

(2.15)

(2.16)

(2.17)

n+1 02A: n+1 _ D n
zi,j,k+l/2 2 zli,j,k+l/2 zlt‘,j,k+l/2

5 A‘ 1/2 1/2
——z zln +At(5xH I“ -6 Hxl'1+ ) (2.18)
2 i,j,k+l/2 yi,j,k+l/2 y i,j,k+1/2

At n+1 n
-— 2| +le
2 i,j,k+l/2 i,j,k+1/2

2Tex n + 1

21 +1
1+ “)0 " -e (8' +5' 00—)
( At Xian/2,1,1: 0 ”5 'x At x|i+l/2,j,k

 

(2.19)

21 ,, ,,
= _ 1 _ 8.!)D + 8 (so _89 J)
( At "Lu/2,131: 0 ”‘5 "‘°° At x|i+1/2,j,k

 

n + 1 2Tey

21' n+1
1+ ”)D — 8 (8' + 8' —)E
( At yli,j+ 1/2,k 0 ”5 ’y°° At yli,j+1/2,k

 

(2.20)

21 ,, 21
= — 1— ‘3ij +8 (8' —e' J)E "
( At yli,j+l/2,k 0 'y‘ 'Y” At yli,j+1/2,k

 

2Tez

21 n+1 n+1
1 + “)D — a (5' + 8' —)
( At Z'i,j,k+1/2 0 '2‘ '2” At zli,j,k+1/2

 

(2.21)

n

21 ,,
= - 1— “)0 +8 (8' —8' J)
( At zli,j,k+1/2 0 m '3” At Zli,j,k+ 1/2

 

If we solve the equations (2.16) to (2.18) and (2.19) to (2.21) simultaneously, we will

obtain the finite difference expression for B and E as follow

1 n n
D "* = D — E
k B” x|i+l/2,j,k 62" x|i+l/2,j,k

x|i+l/2,j,
n+1/2 n+1/2
+ 5 —5H )
B34 y Zita/2,1,1: 7' y|i+1/2,j,k

63x n+1 n
-— J +J )
2 ( x|i+l/2,j,k x|i+1/2,j,k

(2.22)

13

l;

 

»
.lj:.‘n
W

n

n+1 n
- B‘yDyl —BZYEy|t,j+ l/2,k

yli,j+1/2,k i,j+l/2,k

n+1/2 n+1/2
+ 8 H —8 H )
B3y( Z xli,j+l/2,k x zli,j+l/2,k

B3): n+1 n
-— 1|. . + .1. .
2 yl,j+1/2,k - t,}+1/2,k
n+1 n n
= D — E
zli,j,k+1/2 B” zli,j,k+l/2 BZz zli,j,k+1/2

n+1/2 n+1/2
+ 5 H _5 H )
B3Z( x yli,j,k+l/2 y JrIt,j,k+1/2

BBz n+1 n
—— J +J )
2 ( zli,j,k+l/2 zli,j,k+1/2

n
i+ l/2,j,k

n+1

I!
x|i+ 1/2, j,k

= D
Y” xli+l/2,j,k

—Y2xExl

+Y ( n+1/2 6H n+1/2 )
3" y Z|i+l/2,j,k Z y|i+l/2,j,k

Y3 +1 n
-—5(Jx|" +Jx| )
2 i+l/2,j,k i+1/2,j,k

n+1 n n
= D — E
yli,j+l/2,k Y‘y Y11,j+1/2,k YZY yli,j+1/2,k

n+1/2 n+1/2
+Y3y(62Hxl k]

i,j+l/2,k x zli,j+l/2,
Y +1

—fl(J I" +J I"

2 yi,j+l/2,k yi,j+l/2,k

n+1
zt',j,k+1/2

n

= D-" — E
Y” ‘i,j,k+1/2 Yzz zli,j,k+1/2

 

+Y (8 n+1/2 6H n+1/2 )
32 x yli,j,k+l/2 y xli,j,k+l/2

Y3: n+1 n
——(le +le
2 i,j,k+1/2 i,j,k+l/2

where

14

(2.23)

(2.24)

(2.25)

(2.26)

(2.27)

Bio:

B20:

B30:

Yla

72a

73a

anda = x,y,orz.

At . ' 213
Oa(Tea - 3) + 80(8 mm + 8 ,0;wa

 

At . . 21: at
00:61:“: + '2') + 80(8 ras + 8 rawj)

20 1 £08

(1 ea “100

 

At . . 21cc:
Ga 3 +13“! +80 eras+8raoo-E-

21
0 ea
50(5 rats + 5 race-E) At

 

At ' ' 21 a
Gaff + Tea) + 80(8 rats + 8 raooTj')

 

2
At ' ' 217
60(3' + Tea) + 80(8 ras + 8 race-X?)

AI . . 21cc:
Ga '2— + Tea + E"0 8 rats—8 race—AT

 

At . . 21w
Ca —2'+Tea +80 eras+8raoo—AT

At + 21m

 

At . . 21: 0:
001(3 + Tea) + 80(8 ras + 8 ramj)

(2.28)

(2.29)

(2.30)

(2.31)

(2.32)

(2.33)

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 1w is equal to zero, the frequency-

independent case, equations (2.28) to (2.33) become

15

 

 

 

 

 

 

 

808.1113 _ Ga 2
BIG = At’ (2.34)

808),,“ + (Ia-2-
1320. = 0. (2.35)

8 8' At
0 a

B3a = r 5 At 9 (2.36)

80830” + (5Cl 2

2

808m + ca 2

YZQ = 1 . (2.38)
and
At
808',” + (fa—2-
Hence, the equations (2.24) and (2.27) become the same equation as follows:
_ GxAt At
n + l 2808;,“ n 1 808.1}: n + 1 n
x _ —_Ex — -_ xl +Jxl
i+l/2,j,k oxAt i+l/2,j,k 21 0A: i+l/2,j,k i+l/2,j,k
+ —— + —
28 8' 28 8'
0 rxs 0 rxs (2.40)
At
8 8' /
+ 0 rxs (5H,n+12 —82H In+1/2 )
oxAt y “i+1/2,j,k Yul/2,1,1:
l + —,—
280E rxs

16

.‘l‘l
E.

"1‘31 2.1

 

57"

1.11.1 3

GyAt At

 

 

 

(2.41)

2)

(2.42)

n+1 = 28O8 rysE n _l 808 rys ( |n+l +J In )
yi,j+1/2,k oyAt yi,j+1/2,k 2 oyAt yi,j+1/2,k yi,j+1/2,k
E08 rys 2“308 rys
At
EOErys n+1/2 n+1/2
+—-— 523x] —6tz|
oyAt i,j+1/2,k i,j+l/2,k
l+—,——
2808 m
Fifi At
t i
E n+1 _ 28Oerzs n _1 808m: ( n+1 +1 n
Zt,j,k+1/2 ozAt z i,j,k+l/2 2 oZAt Zli,j,k+1/2 Zli,j,k+1/
808 rzs 2808 rzs
At
5 "3' 1/2 1/2
+———° '2‘ (5,11 l“ -—6 Hxl“ )
oZAt yt',j,k+l/2 y t',j,k+1/2
1+—,—-
28Oerzs

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 of Maxwell’s equations are

 

 

n+1/2 _ H n-1/2
"Ii,j+1/2,k+1/2 "lt,j+1/2,k+1/2
E " —E " E " —E "
+_A_t yli,j+l/2,k+1/2 yli,j+l/2,k_ zi,j+l,k+1/2 Zli,j,k+1/2
I“.r ‘ AZ Ay

17

(2.43)

n+1/2

n-1/2
y|i+1/2,j,k+1/2

y|i+1/2,j,k+ 1/2

 

 

 

 

 

 

 

 

 

 

 

E n _ n n —E n (2.44)
2| zl x x
At i+1,j,k+1/2 i,j,k+1/2 i+l/2,j,k+l i+1/2,j,k
+_ .....
”y Ax AZ
n+1/2 n-1/2
ZIi+1/2,j+1/2,k z|i+l/2,j+l/2,k
Exln _ 1"" E n _ E n (2.45)
+A_t i+1/2,j+l,k i+l/2,j,k_ yi+l,j+1/2,k yi+1,j+l/2,k
uz Ay Ax
CxAt At
n + l _ 2EOE rxs n + 808 rxs .
"lt+1/2,j,k oxAt x|i+ 1/2,j,k oxAt
1+——,-— 1+5—-,—
2808 rxs 808 rxs (2.46)
H n+1/2 n+1/2 n+1/2 n+1/2
z|i+l/2,j+l/2,k Zli+1/2,j-1/2,k Zi+l/2,j,k+1/2 Zi+1/2,j,k-1/2
Ay AZ
1— OyAt At
n + l _ 2€08 rys n 808 rys .
yli,j+l/2,k 1 oyAt yli,j+1/2,k 1 oyAt
+ , +——,—
2808 m 2808 m (2.47)
n+1/2 n+1/2 n+1/2 n+1/2
xli,j+l/2,k+l/2 xli,j+1/2,k-l/2_ zi+l/2,j+l/2,k zi-1/2,j+l/2,k
AZ Ax
O'ZAI At
E n+1 _ 2808 rzs n + 808 rzs .
z|i,j,k+1/2 1 ozAt zIt,j,k+1/2 1 ozAt
+ , +———,—
2808 m 2808 m (2.48)
n+1/2 n+1/2 n+1/2 n+1/2
y|i+1/2,j,k+1/2 yli—1/2,j,k+1/2_ xi,j+l/2,k+l/2 xi,j-l/2,k+l/2
Ax Ay

18

 

13,1 StabiJ

The 31‘ (1.51..

   
  

\ U . r
max-37.4.11).

primed 1.)) Ta

that: it 13 1::
sections. res."

tomes. .\o'.:

reg-on Gen-cm

2.2.1 Stability Condition

The stability condition is required to avoid numerical instability in finite difference
approximation schemes. The stability condition of Yee algorithm is first correctly

presented by Taflove[3] and is

At S 1 (2.49)

ch + 1 +-—1—
(Ax)2 my)2 (Az)2

where At is the time step and Ax, Ay, and Az are mesh sizes along the x, y, and z

 

 

directions, respectively. The (2.49) is also called the Courant-Friedrichs-Lewy (CFL)
condition. Note that the CFL condition is derived by assuming a homogenous spatial

region. Generally, the CFL value is defined as follow:

 

+ 1 + 12 (2.50)

 

 

CFL = CAtJ 2 2
(Ax) (4y) (Az)

which is less than or equal to I.

An exact stability condition for general case is usually difficult to derived 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 of iterations with the proper choice of the
time step. In the practical modeling, (2.49) is usually used. If the numerical 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 FDTD grid can differ from the

vacuum speed of light, in fact varying with the modal wavelength, the direction of

19

 

 

uiere c is the x

131g the .t. )1 :

 

 

 

Tnc numcn;
its pimc mm"

ampzcmenzanon

1‘61'13100 Yee a
(rim (0.
LC~\I~ 3:

a” - T .
“5 1,. L . .1

3t ,
\‘1 Sin, \(1
(Pg ' '-

propagation in the grid, and the grid discretization.

The dispersion relation for a plane wave in a continuous lossless medium is simply

(”2 2 2 2
—2 = kx+ky +kZ (2.51)

where c is the speed of light, to is the radian frequency, and kx , ky , kz are wavenumbers

along the x, y, z axes in this medium.

The numerical dispersion equation for FDTD scheme can be obtained by substituting
the plane monochromatic traveling-wave trial solutions into the finite-difference
implementation of Maxwell’s equations. The dispersion equation[16] of mu three-

dimension Yee algorithm is

[isms—~32 = [me-31[84%]1848531

where kx , ky , and kZ are wavenumbers along the x, y, z axes in the computational space.

 

Assume Ax = Ay = A2 = A, choose CFL = a, and define R = A/ A which is the

number of grid cells in one wavelength, then (2.52) can be rewritten as

2 - - 2 - - - 2 - 2
32[sin(a _)] _ Sin(ksrn0cos¢) +Sin(ksrn0$mq>) +Sin(kcosfl) (2.53)
J31; R R R

is speed of wave in

 

where]? = ME/ 2 which is equivalent to ftp/c = n/l-c where 17p

~ ~

computational space and 7:2 = k;2 + ky2 + kzz. The 0 and Q) 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

7—1

4h

'4

(’7)

The)

 

Dunn

cal 1r. ‘

 

The nu:
nems of

firms,~ T}

(1)

(2)
(3)

(4)

(5)

(6)

(7)

The Yee’s FDTD scheme gives a phase error. The speed of wave in the computa-
tional space is less than that in free space. However, (2.52) and (2.53) are indenti-
cal in the limit as At , Ax, Ay, and A2 all go to zero.

The smaller the CFL values, the larger the phase error from Figure 2.2.

The larger the R, the smaller the phase error from Figure 2.3. The effect on reduc-
ing the phase error due to change of R is about 100 time that due to change of
CFL values.

There is a numerical phase velocity anisotropy in Yee algorithm or other FDTD
schemes from Figure 2.3 and Figure 2.4.

The number of 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
FDTD grid from Figure 2.5.

The numerical dispersion occurs because the higher-spatial-frequency compo-
nents of wave propagate more slowly than the lower-spatial-frequency compo-
nents. This numerical dispersion causes pulse broadening due to the spatial low-
pass characters in FDTD schemes.

The numerical dispersion can lead to spurious refraction of propagation numeri-

cal modes if the grid cell size is a function of position in the grid.

21

1.—
L
\
“x. -.
3. .R-

 

.-.
.1“!-

F3‘
C1

(1994

(1992

 

    

 

 

 

 

(199
XXXXXXXXXXXXX
0986 “-
+x
+
+ x . . .

xx : i E 3 3 E 3 3 Xx
0982 l l l l l l l L

0 20 4O 60 80 100 120 140 160 180

Figure 2.2 Variation of the numerical phase velocity with
CFL values at R = 10 and q) = 71/4.

22

 

 

-
.3’
Li‘-
.1-
4’5‘

vp/c

 

 

 

 

 

 

V ,v‘“:.:‘:“——l____T——-—:_l':—-—k‘ ,
,’ +++++ 1 : +++++ \\ ‘
+ :+ , : + .’ + t
’ + : : + ' \-T

099*" 'x' '--~ .- .. .\: _

" + + - + + . x
" ++f ++ j ‘

+ ’ . + i
.++++++.j . I
+ .+ i

0.98” . . .1
+ + f
+ ; + 3

0.97— _q
+ +~

0.96“” + ..+ ......... ..

+ +
+ :+
0.951'"'+" ....+ ._,
+ ~ +
3L+ , ++
094 1 L r r r 1 i i
0 20 4O 60 80 1 00 1 20 140 160 180

Figure 2.3 Variation of the numerical phase velocity with R
at a = 0.9 and (I) = 11/4.

23

 

 

 

 

 

 

1 I l l l I I I I
Q
\ 0.998
E.
I>
0.996
0.994
0.992
0.99
0.988
1 l l l l l l 1
0 20 40 60 80 100 120 140 160 180

Figure 2.4 Variation of the numerical phase velocity
with 4) at a=0.9 and R=10.

24

9’
Hr

.n.

Fl

“1

 

 

 

 

 

 

 

o 7 l l l l l l l l l

' 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

1/R

 

Figure 2.5 Variation of the numerical phase velocity
with I/R at a=0.9 and 0:11/4.

25

13 The 1

From the
from its 5318..
mere (recurs,
intact} pr. '

approrzmtmr.

 

 

 

 

59.13.11 dam .

AI‘GL‘A

mbimme “t

1C” Star
I.

2.3 The Ty(2,4) (FD)2TD Algorithm

From the analysis in section 2.2.2, the phase error in Yee algorithm keeps this scheme
from the calculation EM fields of 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 FDTD formulation. Hence, the fourth-order
spatial derivatives is employed to reduce the 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 FDTD schemes. The larger stencil on the FDTD mesh is very
troublesome when dealing with material discontinuities. However, the implicit 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 Derivatives
The fourth-order finite difference expression can be categorized as explicit schemes
and implicit schemes as follow:

Explicit collocated scheme

 

 

8a 8(“t+1‘“i-1)‘(“i+2‘“i-2)
— E 2.54
( x),- 12Ax ( )
Explicit staggered scheme
(flu) 5 27(“t+ 1/2 ‘ “t— 1/2) " (“H 3/2 ‘ “t—3/2) (2.55)
ax 1 24M

26

t “
[Input 1!

[what .\.’.

 

 

 

The 9073px: I;

1'33 Disper
For 31mph;

(:l‘o:0

Implicit collocated scheme

 

(8—3 43—)
6x i+l Bx i-l 2 Bu ~“t+1‘“i-1
6 7(1): 21. (2'56)

Implicit staggered scheme

 

 

(1“) 43-“)
3x 1+1 ax t-r 11(8u) ~“t+1/2—“t-1/2 (2.57)

24 + E 5} f Ax '
All 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 implicit staggered schemes

of (2.57).

2.3.2 Dispersion Analysis for Explicit Staggered Scheme
For simplicity, Maxwell’s equations in a normalized region of free space with p. = 1 ,

0 , and c = 1 are considered and can be obtained as

8:1,0

<L

ijV = a

(2.58)

Q)

t

where l7 = g + j}? . Apply (2.55) to the left hand side of (2.58) and the central

difference to the right hand side of (2.58) and let

V n _ Voejuhmx + [2.1/3y + 12.1032 - umAt)

I, J, K (2.59)

then the following equations are obtained

27

 

 

3731:”;
1\

“331mm (153

1.5101193, 3
L (A;

H7510“

L

The numenk‘dl I

*L

Fig .
m FI‘SUrQ

5518,.
R The ph.

sin((r)At/2)V _ ,27sin(IZZAz/2) — sin(3IEZAz/2)V
At "II, J, K J 24Ay VII, J, K

+ ,27sin(ltyAy/2) — sin(3/2,Ay/2)V
J 24Ay ZII,J,K

 

(2.60)

 

 

 

,27sin(l:czAz/2) — sin(312zAz/2)V + Si“(‘”A’/2)v
J 24Az Jrl1, J,K At yll,J,K
,27sin(ltxAx/2) — sin(3ltxAx/2) V
_J 24Ax 3'1, J, K

(2.61)
= O

 

.27 sinu'éyAy/z) — sin(3l;yAy/2) V +
‘1 24Ay XII,J.K
,27sin(l;xAx/2) — sin(3/2,Ax/2) sinmmu) _
J 24Ax Vy|1,J,K+ At VZII,J,K _ O

 

(2.62)

 

For non-trivial solution, the determinant of the above equations is set to zero and the

numerical dispersion relation is

 

[sin((oAt/2)]2 = [27sin(ltxAx/2) — sin(3ié,,Ax/2)]2

 

 

 

cAt 24Ax 2 63)
~ ~ ~ ~ . ( .
+ 27 sin(kyAy/2) - sin(3kyAy/2) 2 + 27 sin(k,_Az/2) - sin(3szz/2) 2
[ 24Ay ] [ 24Az ]
The numerical dispersion relation can be further reduced to
1728[ . ( an H2 I: . (l-(Sin9C08¢) . (Bl—(sinecompfl2
— srn — = 27 srn —— — srn
a2 fi R R R

- . . - . . 2

+ [27 sin(k srn 231mb) _ Sin(3km: srn¢)] (2.64)

 

 

+ [2.7 SinCEcosB) _ Sin(37cc050)]2
R R

if A = Ax = Ay = A2.
From Figure 2.6, the phase error of fourth-order scheme is much less than that of Yee

scheme. The phase error of fourth-order scheme at R=5 falls between that of Yee scheme

28

35 R213 amt:

{om-order ~

 

 

l

2.33 Dispel
For mph.
553* aquatic-

$311011 rs

 

‘16:?le ‘4 |\

at R=13 and R=40 as shown in Figure 2.7. Hence, a coarser mesh can be chosen for

fourth-order scheme if the same phase error is required.

2.3.3 Dispersion Analysis for Implicit Staggered Scheme

For implicit staggered scheme given in (2.57) with numerical boundary condition, the
scalar equations of V x V can be rewritten in a matrix form. For BVZ/By, the matrix
equation is

AX = B (2.65)

where A is a (N+1) by (N+l) matrix, N is the number of partition along the y axis,

 

 

 

 

 

 

 

 

 

 

 

 

T

X:[£:Y_.Z 8—1/2 a_.VZ .a_VZ LVZ .a—Vz :l (266)

By j=an j=lay i=2 ayj=N-2ayj=N-rayj=1v

Mi: . .
andB = 737 whereMrsa(N+l) bmeatnx and

A = V ' -V . 2.7
B [V2 -1Vz -3 2 _§ Z _2N-5VZ 2N-3 ’- _ [v-1] (6)

1-2 1-2 1-2 J-——2 —— -——2

The elements of matrices A and M depend on the coefficients in (2.57) and numerical

boundary approximation. The (2.65) can be rewritten as

X = — (2.68)

where A7 = A-IM and A-1 is the inverse of matrix A.

Apply (2.59) to (2.68) and after some calculation, the j component of X becomes

29

(*ch m. '
ill
3:21.950 W

chimed as

 

The (2.70) dcp
equation is W}
(triplex! stag!
'ximher redu;

‘ l». ‘
*mh make 1'

tr ’ ~
£4.th bascc

 

a, v, 44.2"“).
__ : — m. e (2.69)

where 172]." is the (j,n) element of matrix 117. Other derivatives of I7 can be obtained

similarly and a system of equations for V x , Vy, and Vz are setup. Setting the determinant

of above equations to be zero and A = Ax = Ay = Az, the dispersion relation is

obtained as

jlcsin0cos¢(2I-(2n + 1))]2

2+ (37)]— (2;... .

[v-1 jksin0$in¢(21-(2n+ 1)) 2
+ L2 fine R :I . (2.70)

[N_1 jicos0(2K—(2n+1))]2

+ zmxne R
n=0

 

The (2.70) depends on locations of points and matrix [171; hence, the analysis of this

equation is very complex. The phase error in implicit staggered scheme is close to that in
the explicit staggered scheme. Moreover, the phase error in implicit staggered scheme can

be further reduced if the suitable optimization method is used. The optimization methods,
which make I? close to 11 or vp/c close to one while maintaining smaller R, can be

researched based on (2.70).

30

()

 

vp/c

 

 

0.98‘ ..... ....._1

 

0.94,. ..... .lw. .. --

 

 

 

0 20 40 60 80 100 120 140 160 180

Figure 2.6 Comparison of numerical phase velocity
between explicit 4th-order scheme and Yee scheme at
R=5, a=0.24, and 43:11/4.

31

. . _
-

p . w ...
”M" ”fl." .v.
. . . .

a d
.6.
3.
09

v.
. ~

av

.D

 

0.994" '

 

0.993 8 .IV

  

0.992

0.991

 

l

   

  
 
   
 

\

R=5 (4th order) X

R: 1 5 (Yee)

,..... _ ..... R=20(~Yee)-~~ ..
1 ++++++++++z R=40 (Yee) ?

J l l

    

 

 

0.99 4
00

O
8
8%-

Figure 2.7 Comparison of numerical phase velocity

100 120 140

160 180

between explicit 4th-order scheme and Yee scheme at

a=0.24, and ¢=n/4.

32

 

2.3.4 Calculatin
momma

tendered here as

uhzch is close 1‘

an: ,E,
_ ,H ,

 

Made the ph) 5m} 5
sois the value of 615
mesh points other 1h

.4. ‘ .-
mamesofE

a,”

,2 oz;

7"

H\

 

H)

Flgum 2.

x and v.

The f

OHQ ~'
“mg h
On 1}]

2.3.4 Calculation of Derivative of E in 1y(2,4)

The derivatives of E along x, y, and z are first considered. The term dEy/dz case is

considered here as an example. From the Figure 2.8, the evaluation of aEy/azlk _ 1/2
which is close to the boundary by Az/ 2 needs values of dEy/dzlk _ 1/2,

dEy/dzlk E , and Eylk

, if (2.57) is used. However, 3E /Bz| is
=3/2 yk=0 y k=—1/2

outside the physical boundary, hence, the value of aEy/dzlk _ 1/2 has to be approximated,
so is the value of dEy/dzlk _ N “2 where N Z is the number of partition along 2 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 .

By. By

 

 

 

 

 

 

 

 

 

aEy/ dz Ey
=_1,2 ————— k=Nz+l/2
\/ l I 1 V
/\ I I I /\
k=0 k=1/2 k=l k=1/2 ————— k=Nz

Figure 2.8 The lattices of Ey and BEy/az along the z axis at fix

xand y.

The following two fourth-order approximations are used to approximate the

A
derivatives of E on those two special points:

33

“"1
a .I Q”

 

 

 

 

 

 

 

 

 

46:3 :95 "535 1235 - Ey|k=r5y|k=o
24 dz 7 12 dz 5 24 dz 3 12 dz 1 A2 '
*‘i i=5 k=§ *‘2
121.1% __§_a__Ey ,_1_<”__Ey __1_i_Ey
12 dz 1 24 dz 3 12 dz 5 24 dz 7
k=Nz_§ k=Nz-i k=Nz-i k=Nz-i

= Ey|k=N,-Eylk=N,—1

 

 

 

 

 

 

 

Az
A system of equations are derived and written in matrix form,
AX = B
where
r _
26 —5 4 —l O
l 22 1 O O
0 1 22 1 0 O
A =
O O 1 22 l O
0 0 O 1 22 l
_ O O -1 4 —5 26_
X ‘ dz 1 dz 3 dz 3 dz 1
k=§ k=§ k=Nz_§ k=Nz-i
" "i
Ey|k=1-Eylk=o
Eylk =2—Eylk =1
24
B = — ,
Az

EY|k=Nz—l -Eylk = Nz—2

 

 

L Ey|k=NZ—Ey|k=N,—1J

34

(2.71)

(2.72)

(2.73)

(2.74)

(2.75)

(2.76)

id 7 denotes the t'

  

.8
E amen points t'

13.5 Calculatim
The calculation

mitigation uhtn :

 

i=0

“gum

Iand:

If me bOUn t.

(441'

Ck”‘"‘“"¢s of

530mm

and Tdenotes the transpose of a matrix. By using the LU decomposition, the derivatives of

E at every points can be determined.

23.5 Calculation of Derivative of ii in Ty(2,4)
The calculation of de/dy is shown as an example. Figure 2.9 shows the grid

configuration along y direction.

 

 

 

 

 

 

 

 

 

dHZ/dy de/dy de/dy
'=-1/2 Hz j=Ny+1/2
\l 1 i \/
/\ l I /\
j: j=1/2 j=l --------- j=Ny

Figure 2.9 The lattices of Hz and d Hz /d y along the y axis at fix

1: and 2.

If the boundary is PEC, then Ex, E Z, and H y are zero according to Figure 2.1. Hence,

the derivatives of H z on the boundary are zero, also. The matrix equations, AX = B

becomes

P22 1 0 01
1 22 1 0
A = 0 , (2.77)
0 1 22 1
L0 0 1 2g

 

 

35

hm

1’:de

l-‘zd
If [he bOUni
hprOijdtcd‘ T
UC\“t‘i'
BHSI
a?) =
J = 0
and
d[{‘
6T- .
j ‘

 

 

 

 

 

 

 

 

 

de aHz 8H, 3Hz
y j=l y j=2 y j=N,-2 y j=N,—1
- 412'
“j=3/2 j=l/2
H. -Hz|
‘j=5/2 j=3/2
24 -
B = — , .
Ay . (279)
2|. ‘Hvl.
j=N),—3/2 ‘ J=Ny-5/2
H,| —Hz|
i j=N,—l/2 j=N,-3/2_
and
aHz de
d = ’T :0, (2.80)
y j=0 y j=Ny

 

 

If the boundary is not PEC then de/dyl. 0 and de/dyl. N need to be
1: J= ,

approximated. The following fourth-order one-way finite difference approximations are

 

 

 

 

 

 

 

 

 

 

 

 

used:
31 229 75 37 11
——H —H __ __ __
8 2 ._1+ 24 z ._3 8179.;+ SHZ ._7 12HZ ._9
3H2 1‘5 '2 1‘2 1‘5 1‘5
—d_ A (2.81)
y i=0 7
and
3Hz 31 229 75
W _ (7”: _ fiat—”z. 3'r”z._ 5
J N J- )'-i J- y-i J-Ny_§
(2.82)
37 11
*3”: ._ rfi’fi. (1)/M
I'Ny-i J=Ny-_

The components of the corresponding matrix equation are

36

 

 

 

 

-l'aH:1t ‘
- T” :0 I
and
F 31
"Eh
B = 31
A)
L
Ushers

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

24 O 0 O
1 22 1 0
A = 0 , (2.83)
O 1 22 1
_O 0 0 24
7
dHZ dHZ de de de
= _a__ T ._a_ .5— —a—— , (2.84)
y 1:0 y j=l y j=N,—2 y 1=N,1 y }=N
and
_ 31 229 75 37 11 -
7”:._1+Er”z._3‘§”z._5+§'”z._rfi”z _9
1-2 1‘2 J-§ 1-2 1-5
—H
z|j=3/2 Z|j=1/2
—H
24 Zlj=5/2 z|j=3/2
B=— (2.85)
M
-H
Zlj=zv,—3/2 z|j=N,—5/2
Hzlj=N,—l/2_Hzlj=N,-3/2
- f _
where
31 229 75
f—[_8H ._ 1+EHZ._ 3—8 z.__ 5
J‘Ny—i J-Ny—i J-Ny_§
(2.86)
+—H —flH /A
Z. 7 12 z. 9 y
J=Ny—'2' I=NY-§

37

   
  

An CKCtL’iitOfl ;'
'isit}. of the prob.-

A

]: “LUSH” .‘t;

it here

wtt‘) is the time dc

For empty PU
contain a finite cm
the. “it sideloht

{
MB and provides

idilO“

 

2.4 Excitation Source and Power Analysis

2.4.1 Excitation Source

All excitation probe is placed parallelly to one of the rectangular axes. The current

density of the probe is modeled by induced EMF method and represented as
J = (1(X)5(y - yo)5(z - zo), 5(x — xo)l(y)5(z — 20), 5(x - xo)5(y - yo)1(2))W(t) (287)
where

sink(ha — 0t)
sin kha

 

1(a) = 10 CL = x, y,z, (2.88)

w(t) is the time dependence of .7 and ha is the length of probe along the on 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 of the Blackman-Harris (BH) window[20] is approximately -92
dB and provides a smooth transition of excitation. The BH function is discretized as

follow,

1I(tn — ta)
Nhalf
2115(tn —

tc)
W) In e [tc — N ha”, tc 4' N half] (2.89)

31t(tn - 5))
Nhalf
0 otherwise

0.35875 + 0.48829 cos(

h(t ) = , +0.14128cos(
n

+ 0.0ll68cos(

 

where w(t) = h(t)coscot, In is the time step, 16 is the center of the BH window, and

N h a If is the designed half width of the BH window function.

38

For a loaded

 

frequency Sinusozt.

   
  

the finite value 0t]

The excited mt.

“here the Green's t.

 

m2 denmes the
n

le‘CS. the probe mm

7

.file integm} {Own

Fifi" =

For a loaded PEC rectangular or an aperture-coupled cascaded cavity, a single
frequency sinusoidal source is used. The loss of material sample in the cavity accounts for
the finite value of Q.

The excited modes are determined by the location of the probe according to

750) = [imam #0)de (2.90)

where the Green’s function is

_ , , , E ? E ?
Ger.) = —qu2 "‘ g” 3‘ ’ <29»
n kn_k

 

and 2 denotes the sum of all the TEnm] and TMnm] modes. For example, to excite TM

n

modes, the probe must have JZ component only.

2.4.2 Power Analysis

The integral form of the Poynting’s theorem is

)dV + fa? x I?) . fidS (2.92)
V V V S

A
where J 3 is the source current. For lossless case, the power stored inside a cavity after the

source turned off is

_ A 313‘ _. 319‘
P5 — {(E W441 T)dv. (2.93)

The FDTD approximation of the first integrand, evaluated at time step 72 , in (2.93) is

39

'- I
h-i
II

and that of the SCCU

“Lia '
““3 tne “Onm'dw

me step 71 is

P2]; = _2_1A_t{Exln(Dx|n+ l—Dxln-l) +Ey|n(Dy|n+l_DyIn—l)

(2.94)
n n+1 n-l
+EZ| (Dzl —Dz| )}
and that of the second integrand, evaluated at time step n, is
H 71+12 71—12 71+12 71-12
n __ 0 2 2 2 2
2 1 2 (2.95)

N

that“? W l ll

where the nonmagnetic material is assumed. Hence, the power stored inside the cavity at

time step n is

P: z 2(1):} + PZH)Av (2.96)

where Av = AxAyAz for uniform partitions.

For material with conductive loss, the dissipated power is the first term on the right-

hand side of (2.92). The conductive dissipated power density is 0E2 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 of (2.92) and the dielectric
dissipated power density in the time domain has to be extracted. For this case, (2.93) is the
summation of the stored power of the whole cavity and dissipated power due to the lossy
material sample. Hence, only the calculation of 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.

4O

 

Assume the 03'

there 8(0)) and E

 

sign 8. stands for t
M, and the sewn

relation of the perm

rewritten as stored a

and

Assume the one-relaxation Debye medium, the Fourier transform of the x component

inE-Q—lg-is

d

2111(Ex(m)® [jw8'(m)Ex(‘”)l + Elm) 8’ [(1)8"(U))Ex(m)]) (2.97)

where 8'((1)) and 8"(w) are the real and imaginary parts of electric permittivity and the
sign, ®, stands for the convolution. The first term in (2.97) is the stored power density due
to E x and the second term is the dissipated power per volume due to E x. Substituting the

relation of the permittivity with static and infinity frequency permittivities, (2.97) can be

rewritten as stored and dissipated powers as

-1- E (t0)®[jwe' E ((0)]+E ((1))®l:(8' —8' )jwfl] (298)
2112 x x00 x x xs xoo 1+(‘Mexl2 '

and

1 (8x; - 8x4) Ex(w)
ET(EJ‘((D) ® Ex(00) — Ex((0) ® [WD. (299)

The corresponding time domain equations for above equations are

 

dE d
e'mexU) it“) + Ex(t)l:(€'xs — Edd-£8179] (2.100)
and
(E'xs _ 8.x”)
——(E.(r)E.(t) — Hangar» (2.1m)

ex

where gx(t) is the inverse Fourier transform of Ex((1))/ [1 +(wtex)2]. Hence, the

41

and

I
there

Polnttng theorem

 

Poynting theorem for nonmagnetic material can be rewritten as

133-jaw = I[o]E-EdV+ laggE-E—[ggicgmv
V V

V.
_. a}: 4 a‘ _. 88‘ _. _. .
+1([§1]E°W+[§2]E'-£+H'WJdV+§(EXH)°ndS
V s
where
8',” O 0
[£1]: 0 8'”° O
O 0 s'zoo
(e'xs—e.’too O 0
[£2] = 0 Eds-6y“
0 E'zs-fi'zm
ox 0 0
[0] = 0 0y 0 .
O 0 oz
and
8' —8' —
XS 10° 0 O
Tex
e', -8' 00
[£3] = 0 _£_L O
'c
ey
o 0 Sis-8'2”
L Tez J

 

 

(2.102)

(2.103)

(2.104)

(2.105)

(2.106)

where the V5 is the volume of the material sample. Note that the summation of the third

term in (2.102) can be obtained by subtracting the second term in (2.102) from (2.93).

42

 

Etrluating the
l.ll+(u)t,)‘].
impractical in tli.
agp-roximuted b) .

recursise equation 3

parte‘“ and 20. t

 

, Ar
gltnht = \L
) 4: I

(.1

Hence. the stored pt

version of (2100) a

12.101.) and (2.107).

Evaluating the gx(t) involves the convolution of E x(t) and the time domain version of
1/[1+((1)1.'e)2]. To get the convolution, a history of Ex(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 g x(t) into a
recursive equation by discreting the convolution directly. By using the Fourier transform
pair, e‘al’l and 2a/(a2 + 002) , the discrete version of 3,0) is

gx(nAt) = ZALExMAt) + Zé—t—KA’MHExfln — 1)At) + e'N/‘agxan — 1)At). (2.107)

Tex ex

Hence, the stored power density due to Ex at time step 71 is calculated from the discrete

version of (2.100) and (2.107) and that of dissipated power from discrete version of
(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(w) = Ex((0)/[l+(mrex)2] or [1+(m'rex)2]gx(00) = Ex((u),

then its time domain expression is

6128.0)

dt2 = Ex(t). (2.108)

 

2
gx(t) —Tex

By using the second order central difference approximation of the second derivative,

(2.108) becomes

gx((n + 1 )At) = [2 + ($1)2]gx(nm) — gx((n — 1)At) — (gjzzzxmm). (2.109)

ex ex

The (2.107) is obvious first order accurate but (2.109) is second order accurate in time.

/1'

However, the homogenous solution of (2.108) contains a divergent term, e‘ a, and its

43

fDI'D tort

general DC'

Results
response at
agenthm.
limitation t
obsersattun
of iterations
the timed):
of netting
T : Nil. .
the true 5pm-

The Pror

FDTD formulation cannot be used. More general FDTD formulation of power for a

general Debye medium will be discussed in the next chapter.

2.5 Prony’s Method

Results in the frequency domain are usually obtained by recording the time-domain
response at the selected observation points of the FDTD mesh and applying the FFT
algorithm. There are, however, several limitations in the FFT approach. The main

limitation is that of the frequency resolution Af , which is roughly the reciprocal of

observation time, i.e., A f = 1/ (N At) where At is the time step and N the total number
of iterations used in the FDTD method. A second limitation arises from the windowing of
the time-domain data because the FDTD response is truncated in time. This has the effect

of viewing the true time-domain response through a rectangular windows of duration

T = N At. In the frequency domain this windowing is translated into the convolution of

the true spectrum with the function sin f / f .

The Prony’s method is a technique for modeling sampled data as a linear combination
of complex exponentials and is particularly suitable for calculating the resonant frequency
and Q of a resonating structure, since the impulse response of such a structure is

characterized by a superposition of decaying exponentials[21].

2.5.1 Theory

The Prony’s method fits an exponential approximation with unknown exponents of the

form

fin:

toafunction ft!)
of variables has
771: 0.1.2. ...rt—

caution can be m

0111') me main; ion

 

“here f. -

alarm“

ae

We
Can get me ft

p
f(m) = CleAim+CzeA2m+ "'+CQeA°m __, 2 Ck“? (2.110)
k=l

to a function f (t) by sampling at n equally spaced points and 11k = eA". A linear change
of variables has been introduced in advance such that the data points are

m: 0, l, 2, ...n — 1. The coefficients Ck and ilk can be complex numbers. The above

equation can be written as follow,

C1+C2+...+CP fo
C1111+C2u2+...+Cpp.p = fl
C11fi+C211§+...+Cp1lfJ = f2 (2.111)

-1 -1 -l
C1117 +C2u; +...+Cpu';, =fn_1

or in the matrix form

1 1 1 1 Cl f0
111 112 113 11,, C2 f1

2 2 2 2
11.1 “2 “3 “p C3 = f2 (2'112)

 

 

 

 

 

 

where fl. = f(i) and i = O, 1, 2, ..., n — 1. Let 111, ..., up be the roots of the following
equation

p
—1 -2 k
11P+0t111p +0t211” +...+Otp_1|.t+(1p = Eaku = O (2.113)
k=1

then we can get the following equations from (2.111) and (2.113)

45

 

or the follou tr,

equal sign

 

fp+fp-lal+fp—2a2+”.+f0ap = O

f +f0t+f_0t+---+f0t=0
p+1 p l p l 2 l p (2.114)

fn-l +fn—2al +fn_30€2+ +f,,_,,_,0tp = 0
or the following one by shifting the first terms on the left hand side to right side of

equal sign

 

 

018-2 f. 30,- ”f,“

fp fp-l fl (11 fp+l

fpfl f'p f.2 :2 = (-1) fp'+2 (2.115)
f ' f. f ' 519 f.

L 71-2 71-3 n—p—l_ _ n-l_

 

 

 

 

which is linear in 0t.

The Prony’s method solves the (2.110), a nonlinear equation in p. , by solving (1,. in
(2.114), a linear equation in 01, finding the roots, 111., in (2.113), and obtaining C, in

(2.111) where i = 1, 2, 3, ..., p.

The original Prony method chooses n = 2p + l . The original Prony method works
perfectly well 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 of root locations to the coefficients of the polynomial
in (2.113). The least squares Prony method uses n > 2p + l , then (2.112) and (2.115)
become overdetermined linear equations systems and the singular value decomposition
(SVD) is used to solve those equations.

Determining p is an important issue of applying Prony method to frequency domain

46

analysis or pred:
less than the nu
poorill]. Hem C\
modes introduce

rne'iiod With the :

2.5.2 Estimati-
Method

From the detir

the FDTD time-d.

superposition of th

n7

“1th ..
k (”k/(WC

tEOblalnedle (In

and

“here .
t 18 t
he Calm

analysis or prediction since every p,- represents a frequency component. If the value of 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 of Prony’s

method with the time sequence of the sample points in reverse order[22].

2.5.2 Estimation of Resonant Frequency and Quality Factor by Prony
Method

From the definition of quality factor,

Q _ Stored Energy

_ 2.116
0 Power loss ( )

 

the FDTD time-domain response, say Ex, of a cavity can be expressed in terms of a

superposition of the resonant modes[24]

‘0 —(a. + jznmmm p
Ex(mAt) = z Cke = 2 C1111" (2.117)
k: l k=1

with ak = (Dk/(ZQ) where m = (O, ..., )N— 1 . By Prony method, the Ck and 1.1,, can

be obtained if E x( mAt) are known. Hence, the resonant frequency and damping factor are

 

 

ima8(ln(uk))
: 2. 1
k 21tAt ( 1 8)
and
real 1n
0‘1 = (Atmm (2.119)

where 11k is the calculated mode corresponding to f k, imag(ln(uk)) is the imaginary

47

  

part of lnilh) ' ar‘
Once the dam?

factorean be 641511.-

 

Xote that (2120) ts
Prony method pro

distribution at trans

2.6 A Single E

An empty recto
it '~ '
0. his cavity is int2

he .
results presente-

2.61
Configur
The configurut
excitation probe is

aéSleed
to be
Vern

2.6.

‘) »
‘Numeric;
For
TECH m(.)(

l
r al2.45
GHZ is i

the
l

ind
£1

part of ln(1lk) , and real(ln(1tk)) is the real part of ln(11k).

Once the damping factor and the resonant frequency have been determined, the quality

factor can be easily obtained as

nfk
at

Q = . (2.120)

 

 

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 of cavities from the field

distribution at transient state.

2.6 A Single Empty Cavity with PEC Walls

An empty rectangular cavity with PEC boundary is studied in this section. The Q value
of 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 of the six faces of the rectangular cavity and its length is

assumed to be very small. The dimensions of the cavity are: b1xb2xb3.

2.6.2 Numerical Results of Field Distributions

For T1301 1 mode, the excitation probe is located at the center of the y-z plane and only

I, at 2.45GHz is provided. The BH source is used and automatically turned off at 18,875

time steps, spanning over 0.2516118, with At = 13.3299 ps. The cubic cavity with b1, b2,

and b3 all equal to 0.08658 meter is assumed. For this dimension and frequency, the

48

91

gr

)0:

m-
15‘

TElOl, T1301 1, and TMl 10 modes can be excited; however, only THO” can be excited
because only J x is provided on the probe and E xs of other modes are zero. The numbers of
partitions along the x, y, and z are all 10.

The x, y, and z dependences of 13x are shown in Figure 2.11 to Figure 2.15. For TED“

mode, Ex is constant along x at fixed 2 and the numerical results are shown in Figure 2.11

and Figure 2.12. The Ex is proportional to sin(7t01/0.08658) along the Qt axis where 0t 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 Dsin(7t0t/0.08658) where D is a
amplitude factor if 0t is from 0 to 0.08658 with increasing 0.008658 for every step. The
summations of E, E2, and Hx over all grid points are 3.7386e-06, 1.9023e-06, and 0
respectively and it shows that no Ey, E2, and HK exist.

For ml 11 mode, the probe is located at the center on the bottom x-y plane and only Jz
is provided since Ez of TElll modes is zero. The length of the cubic cavity is 0.10604
meter and the BH source is off at 025131118, 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 of Hz over all grid points is -5.7932e-31 which is almost

zero, thus confirms the excited dominant mode being a TM mode.

The time dependence of Ex for TE), 1 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 TE011

mode is observed to be close to 2.5 GHz. However, the estimated resonant frequency of

THO“ mode by Prony’s method is 2.4528 GHz which is much closer to the source

frequency. The estimated resonant frequency of ml 11 is 2.4548 GHz by Prony method.

49

 

 

 

/ A
. y \
b3
5 .

 

 

 

 

 

’ b1

 

 

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

 

 

 

 

 

 

 

1 l l I W l l l l
0 __ y=0 y=0.08658 -
-1 - -<
-2 _ a
_ _ x=Q-OQB§5§ _________________ #3327222. _ _ _ -
-3 .. a
-4 h- ..(
‘5' ._._y=9-.0_173_16_._.-._ _._._._._._._._._._._Y=_0:®9§64_._._._ ‘
-6 ~ 2
"7 ’ 1 #0025974 r . r 1 . Fit-060606. . ‘
-3 _ 4
‘L 1:0:934632 :% = x 7‘: #4 11051948: %
9 _ ...... 14994329 ................................................ y.=.o.943.2.9. ...........
_10 l l l l l l l l
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
x (m)

Figure 2.11 The x dependence of Ex at z= 0.025974 and t= 0.26658 us
for TEOII mode.

51

Ex

-6

—8

 

 

 

 

 

 

 

T l j l l l I 1
2:0 z=0.08658 y
_ ._ _. _Z=_0-0_08_55_8 ________________ £45373”. .. _. _ -
._._._Z.=_0-.0_17.3_1.5_._._._._._._._._._._._._._Z-'=_0.-0_59254_._._._ ‘
, 250025914 1 1 1 , z=0.060606, , q
x $300346? .. x x x 2:9.051948;: x ..
......... 2:0043292=004329
1 l I; 4% L l l l
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
x(m)

Figure 2.12 The x dependence of Ex at y= 0.034632 and t= 0.26658 1.18
for TEou mode.

52

Ex

-2

-4

-10

-12

 

 

 

l

l

‘ -—'
§ ‘ .—

1

l

 

z=0
z=0.008658
z=0.017316

+——+—+ z=0.025974
x—-x—x z=0.034632

z=0.04329
l

   

 

 

0 0.01

0.02

0.03

0.04

Y (m)

0.05

0.06

0.07

0.08

0.09

Figure 2.13 The y dependence of Ex at x= 0. 025974 and 1: 0.26658 us
for TEOII mode.

53

Fl

Ex

 

 

 

   

 

 

 

W T T I T I V 7
I ..
\
-2 >- .1
-4 _ a
-6 _ 4
-8 _ -
y=0
- - - - y=0.008658
-10- ------ y=0.017316 -
+—s——1- y=0.025974
x—N—x y=0.034632
......... y=0.04329
_12 4 1 l l L 1 1 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
z (m)

Figure 2.14 The 2 dependence of Ex at x: 0. 034632 and k 0.26658 113

for TEou mode.

54

0.09

5/\ 'l
‘ Kirk

Fig“

0361

 

 

 

 

 

0'08 r- <— <—< < / < <—— <— <— q
0.07 r- .— <__< < \ < <__ ._ .4
0.06 "‘ ‘— <—'< \ < < 9 ‘— .1
0 05 ‘— <—< < < < < <— .—

' ” i

 

 

 

 

0.04 - 4
‘— <— < < < < e— <— ‘—

0.03 e a
«— <—<—<—<<———é—<—- <—-— «—

0.02 - -
.— <—<—e<—<——<—— <——— ‘—

0.01 r ‘_ 6 < < <_ +9 (.— ._ ‘

0 - .— <—<——é—<——<——<— <——— «— ~

_001 1 1 l 1 1 1 1 1 1
-0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

2

Figure 2.16 The E field on the xz plane at y= 0.034632 and t= 0.26658 us
for TEOII mode.

56

 

 

 

 

  

 

6 l r I 1 ‘r l l l r F
4 _
2 P-
111 o r
-2 -
y=0
y=0.01064
y=0.02128
-4 - y=0.03192
y=0.04256
y=0.0532
-6 1. i i . 1 1 1 1 1 1 1 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
x (m)

Figure 2.17 The x dependence of E" at z= 0.03192 and t= 0.31019 us

for TMlll mode.

57

 

 

ray

 

.................. y=0.0532
+ 4 1 14% y=0.03192
._._._._._.- y=0.02128
_______ y=0.01064
2- y=0

 

 

 

 

1 l 1 l l l 1 l l l

 

 

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
x (m)

Figure 2.18 The x dependence of Ey at 7,: 0. 04256 and t= 0.31019 us
for TMlll mode.

58

 

 

 

z=0

 

 

 

 

-4 ' ——————— z=0.01064
~—1------—-- z=0.02128
1 1 1 1 1 1 z=0.03192
-6 - x—x—x—x—x—x z=0.04256
l l 1 L l l ........ 1 ........ .1 2:0-p532 l
O 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
X (m)

Figure 2.19 The x dependence of E2 at y= 0.04256 and t= 0.31019 us
for TMlu mode.

59

 

F 1'

 

0.12

,1ngTT11
0“ ._\'\I‘ ff//_.
e\\\iT////fi
”a“ <___«—\\ ‘\ i f / ///7__>
é_-s:\\‘\\\/'/’/’/...—a__>
0.06 6*ka /”fi’9_—9L_>—
9K///\\\\\>——>
0.02- @///¢x\\\\—’
”M1 1M“

0” 1i¢¢ii¢¢1

 

 

 

Figure 2.20 The E field on the xy plane at 7.: 0. 04265 and t= 0.31019 us
for TMIII mode.

60

 

 

 

 

0.12 I l I I l l
1 T T T T T T 1 t
0'1” _. / / f f T \ \ \ \ ._
0.08-
____>.———;.’—? /' / \ \ KW—
0.06"
—9_.>_——>_.—’ I x “W
W W \ I K é—<——e
0.04- é”
a%\ \ \ l / (gee
0.02- —* N \ \ \ l / / / 4” ‘—
* \- \ 1 1 1 1 / / 1’ ‘—
°' * 1 1 1 1 1 1 1 *
-0.02 1 1 1 1 1 1
-0.02 0 0.02 0.04 0.06 0.08 0.1

Figure 2.21 The E field on the xz plane at y= 0.04265 and b 0.31019 us
for TMlll mode.

61

 

0.12 fii I I r I I

. 1 1‘ T T T T 1 .
0_1_ _. / / f T l \ \ \ x ._ ‘
_,/r// f ’1 \\\<\‘__
”8’ fizz/v / \ \\<\<\<_ 1

0.06 ~ 1

0.04 ~ 4

afix \ \ / / (“é—“‘6—

 

 

 

o.021 —’ N \ \ \ I / / / e” *— -
" \ \ 1 1 1 1 / / ’ “
0" I 1 l 1L l l 1 1 I ‘
493.02 o 0.02 0.04 0.06 0.08 0.1 0.12
V

Figure 2.22 The E field on the yz plane at x= 0. 04265 and t= 0.31019 us
for TMlll mode.

62

Ex

 

15-

10-

 

.5-

 

 

-15...

 

L g 1

l

 

 

 

 

l l l l

 

0.247 0.2475 0.248

0.2485
time

0.249 0.2495 0.25 0.2505

Figure 2.23 The time variation of Ex at x=0.04329, y=0.04329, and

z=0.034632 for TE011 mode.

63

 

2500 I I I I I T j I I

 

 

 

 

 

 

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Frequency in GHz

Figure 2.24 The frequency response of Figure 2.23.

(if;

TB;

2.7

Ci?

2.7 A Lossless Loaded Cavity with PEC Walls

In this section, the numerical computations of 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 of the rectangular cavity are shown in Figure 2.25.
The excitation probe is located at the center of the xz plane and is very short comparing to

the dimension of the cavity. The excited field is TE 10 1 mode and the operation frequency is

2.4SGHz with the wavelength )1. equal to 0.12245m when there is no material sample
present. In order to compare with theoretical estimations, three material 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 of the rectangular cavity. The dimensions of the material sample are set to be

x0=0.00343m, y0=0.0034m, and 20:0.00352m. This sample is assumed to be lossless and
have the relative permittivity of 8, = 2.5. In this computation, the number of 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

of nodes starts from 0. The excitation source spans over about 67,412 time steps with
At = 3.72822 ps and the computation stops at 68000 time steps.

With this electrically very small sample, xo/A = 0.028, the static electric field

induced inside of a dielectric sphere, Em = (3/ (2 + 8,))E is used to estimate the

ns’

induced electric field in the sample where Em is the electric field inside the dielectric

65

 

 

 

 

 

 

 

 

 

 

 

 

0.072m

l

l

l

l

l

: : X0

1 x 81110 yo 2

I

._____-_z(,-:l .......... —>

,r”0 1 01163m

80.110
0.034m

Figure 2.25 Dimensions of the rectangular cavity and the loaded material
sample. The center of the material sample is consistent with the center of the
cavity.

66

sphere and E M is that in the region of sphere before the dielectric sphere is placed. For

TE ,0 1 mode, only Ey exists since Ex and Ez are zeros for this mode.

The variation of Ey along y axis at x=0.034286m, z=0.056388m and at t=0.25352us is
plotted in Figure 2.26 and the calculated ratios of E w 5/ E n 5 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 n s is approximated by the electric field at first node because Ey is constant along

y axis in TE ,0 1 mode when there is no sample present. The closeness of the numerical
results and the electrostatic estimation gives confidence to the numerical accuracy. The

ratios of E m/ En 5 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 E W s/ E n 3 ratio at different times

 

 

 

 

 

 

Time Step Ews/Ens at node 9 Ews/Em at node 10
67999 0.6550 0.6352
67979 0.6550 0.6352
67959 0.6550 0.6352
67951 0.6550 0.6352

 

 

 

 

 

Due to the induced charge on the material sample surface and the induced current in the
material sample, the other components of the electric field are induced to satisfy the

boundary conditions. The Ex, By, and EZ along the y axis are shown in Figure 2.26 and

67

shows that the excited cavity mode is not a pure T1510, mode anymore since there are Ex
and Ez inside the material sample. However, the TE 10, mode still dominates inside the
material sample judging from the amplitudes of Ey versus Ex and Ez in Figure 2.26.

The estimated resonant frequency of TElm mode for PEC empty cavity with the

dimensions shown in Figure 2.25 is 2.4532 GHz while that for PEC loaded cavity with a
quasi-square cubic sample is 2.4478 GHz. 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 of a thin square plate, having its height much smaller
than its width, is placed in the center of the cavity. The numbers of partition in this FDTD
calculation are 15x17x10 and the dimensions of 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 permittivity of e, = 2.5.

The x dependence of Ey is plotted in Figure 2.27 and then a sine function is used to fit
that curve. By this way, the E m 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 of E m = (1/8,)Ens . The ratios of Ews/ E n 5 along the x direction for

different locations of z are plotted in Figure 2.28. Those ratios inside the material sample

varies from 0.39 to 0.45 which are close to the electrostatic estimation of 0.4. The ratios of

E m/ En s varies slightly at different times for this case.

68

Electric fields

 

0.6 T I I 7 I I

x=0.034286
0.5 - z=0.056388

 

 

 

 

 

 

0.1 ~ ‘1
c 1 1 1 1 1 1 1 i 1 1 = 1 fi‘ -
_01 I 1 I 41 g l
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035
V (m)

Figure 2.26 The variation of electric fields in the y direction at
t = 0.2535 us .

69

 

0.16 r I I I I T I

0.14“

0.12-

0.1+-

I

117008

      

111—111—111—H—1u Calcualted Ey

   

e—e—e—e—e—e Fitted Ey

  

0.02

 

 

 

o 0.01 0.02 0.03 0.04 0.05 0.06 0.07
x (m)

Figure 2.27 The variation of Ey in the x direction at t = 0.24289 us .

The line with star symbol is the calculated values and the line with
circle symbol is the fitted values for the empty cavity.

70

 

0.6

0.56

0.54 '-

EWS/EHS

0.52
0.5 '

0.48 ..1. .
0.46

0.44

 

 

 

0.42.... ............ I. ................ 2:0058‘15" ............. .............. .1

 

 

 

0.02 0.025 0.03 0.035 0.04 0.045 0.05

Figure 2.28 Variations of E m/ Ens in the x-directions. Each curve
represents this ratio as a function of x for different locations of z. The
relative permittivity of the thin square plate material sample is e, = 2.5 .
The solid line with symbols is the ratios at t = 0.242891“ and the dash
line with symbols is those at t = 0.24301 us .

71

The excitation source is off at 0.239118 which covers 49,049 time steps with

At = 4.862 ps and the computation stops at 50,000 time steps. The calculated resonant
frequency is 2.4382 GHz and the frequency shift is 0.61% comparing to 2.4532 GHz

which is the resonant frequency in empty cavity. The maximum of Ex and EZ inside the

material sample is in the order of 10'4 which is very small comparing to the E). This is

because the material sample is very thin in the y direction so that there is no significant

induced field Ex and Ez.

2.7.3 Narrow Strip Case

The dimensions of this narrow strip sample are 0.002322m, 0.02125m, and 0.002115m
along the x, y, and 2 directions. This sample is also assumed to be iossless and have the
relative permittivity of 8, = 2.5. The number of partitions of this FDTD calculation is
31x16x55.

Theoretical estimation of the induced electric field in the material sample may be close
to the electric field when the cavity is empty because the initial electric field is tangential

to the major part of the material sample surface, and the continuity of the tangential

component of the electric field at the material sample surface requires this estimation. The

variation of electric field along y is plotted in Figure 2.29 and the ratios of Ews/EM are
shown in Figure 2.30. The induced EJr and Ez fields inside the material sample are almost
zero and Ey is the dominant field. The ratios of Ew s/ E n 3 inside the material sample range

from 0.94 to 0.99 which are close to the theoretical estimation of 1.

The excitation source is off at 74,835 time steps which spans over 025133113 with

72

At = 3.35841 ps , and the computation stops at 75,000 time steps. The calculated resonant
frequency for this narrow strip case is 2.4464 GHz and the frequency shift is 0.28% from

the resonant frequency of the empty cavity of 2.4532 GHz.

73

Electric fields

 

 

 

 

 

0.7 I I I I I I
0.6 ~ -
0.5 - 3
0.4 >- -
x=0.034839
0.3 L z=0.057093 -
------ Ex
0.2 - i—l—fi—I—H Ey -
~1——1—+—+—+—+ Ez
0.1 P 1
GW I -:
_O.1 l I l l l l
0 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
at x=0.34839m, z=0.057093m, and t=0.252uS.

74

0.035

 

‘25 I I T f I I

 

 

Ews/Ens

z=0.057093

1i ............ ........... ............. . ............ ............. .........

 

115?. ....... . . .....

....... d

1.1... ......... HX=OO34839 ......... ........ .....

095.. ............. ............. ............ .....
03... ....... g .......... g ....... g ............ g ............. .....

0.85” ............. ............ ...... ........... ............. 3 ............

 

 

03 i 1 i 1 i '1
0 . .

Figure 2.30 The ratios of E m/ Ens in y directions.

75

0.035

SO

12*.
(IO

15

2.8 A Lossy Dielectric Loaded Cavity with PEC Walls

The field distributions of a cavity loaded with a lossy dielectric material sample with
PEC wall 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 of material sample is chosen
large enough to cause a relatively low quality factor of the cavity. The relation of cavity
quality factor and number of the time steps which is needed to reach the steady state is
also studied in this section. In order to identify the time steps to reach the steady state, the

FDTD formulation of power analysis is used.

2.8.1 Configurations

The physical configuration of 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 15, I7, and 10 and
the dimensions of the material sample are: 0. 0336x0. 014x0. 0698 m.

The real and imaginary parts of the relative permittivity in (2.3) can be obtained as

 

e —e' ..
e',(m) = e'm+—”—-’—2- (2.121)
1+(mre)
e' -e' (or
8",(w) = ( " ’°°) 2 " (2.122)
1 + (one)

and the 8', and 8", can be expressed as follow:

 

8'” = e',(m)+wtee", (2.123)
I I r

e = e (n — . 2.124
roe r( ) (D1.- ( )

6’

For the lossy material sample, the relaxation time, re , is not zero and it depends on the

76

properties of the material. In this FDTD computation, re is assumed to be 10‘9 seconds

and the 8', is set to 2.5 with four different 8", as listed in Table 2.2.

Table 2.2 Permittivities mapping at re = lns , 0) = 2n(2.45 e9)

 

 

 

 

 

 

e',((1)) 8",(w) 8'” 8',”
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 of the Ty(2,4) (FD)2TD 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.8622ps. Note that a smaller At

needs to be used if a small 8",(w) is chosen.

2.8.2 Numerical Results and Discussions

The approximate time 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’,((n) = 2.5 and 6",(tu) = 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

the induced EMF method used is not accurate to calculate the input power at the probe
location. The stored energy is calculated from the integral of the difference of 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 of this fitting data and the dissipated
power data in the period of 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

 

 

 

 

 

 

 

 

 

,, The average loss The time-average
e .(w)
power stored energy
0.1 2.45e-10 8.905e— 1 8
0.5 1.1258e-9 8.86-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 of time steps to reach the steady state is roughly equal to the Q

78

Table 2.4 Properties of four different lossy cases

 

 

 

 

 

 

 

 

 

 

 

 

Major Q .
w 1» $223212; 1:13:21 9= 08—11151? Prony's 23.322221}?
(GHz) Method
0.1 2.4503 0.12 559.5 523.6 25,000
0.5 2.45 0.13 120.3 125 5,000
2.5 2.45 0.13 38.3 37.1 1,400
5 2.45 0.13 35.7 34.2 1,300

 

factor divided by twice the resonant frequency times the period of one time step. Hence, a

very long time integration is needed if a high Q cavity is dealt with. 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 of Ey in the steady state and the fitting sine function are plotted in

Figure 2.37. The ratios of 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 middle point of 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

 

 

Run FDTD codes for a chosen tim
steps with BH windowed sinusoidal
source.

 

 

 

Use Prony’ 8 method to determine
the major resonant frequency and
the Q factor for this cavity.

 

 

 

Estimate the time steps to reach

steady state solution.

Run FDTD codes for 1/4 to 1/3
time steps of estimated steady state
time steps with a continuous

sinusoidal source.

I

 

 

 

 

Estimate the field distributiona
steady state by Prony’ s methodo
other estimation techniques.

 

Figure 2.31 The flow chart of calculating the field distribution at
steady state for high Q cavities.

80

 

1.31 ......... ........... ....... ., _

‘6... . . . f ......... .. . ‘. . . . . . . '. . . ............ .. . . . . . . . . . ...: ........ .. ..
I - - .
. . I

 

1M4,” 111-1

121.. ....... 1.....LW. 1F

Storedenefay

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . .
06'1“" .. “1......” .1, ....

 

 

0.4b.......1.. ........ 1......p 4...1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . . .
. . . . 1
. . . . . .

02p ...... . . ................ “. . . . . . .~. . ......... -..... . . . . x . . . . . .... ........ . .1
I a I ~ v n v
. . . . . .
. . . . . .

 

o i '. ‘1 i i 1‘
0 0.5 1 1.5 2 2.5 3 3.5
Time steps x 104

Figure 2.32 Stored energy for the cavity with material sample of
8',(m) = 2.5 and 8",(w) = 0.1 . Note that one time step is equal to

4.25442 ps and this is a downsampling plot.

81

9
on

Dissipated power
0
b:

0.4

0.2

 

0
0 0.5 1 1 .5 2 2.5 3 3.5
Time steps x 10

Figure 2.33 Instantaneous dissipated power for the cavity with
material sample of s',(m) = 2.5 and 8",(w) = 0.1. Note that one

time step is equal to 4.25442 ps and this is a downsampling plot.

82

 

1.6~ W

1.4

 

1.2

d

Stored Energy
3

0.4

0.2

 

................................

l l l

.........................................................................................

 

 

 

O
3.39 3.391 3.392 3.393

3.394

3.395 3.396 3.397
Time Steps

3.399 3.4
x 10‘

Figure 2.34 The plot of the fitting and calculated data for stored
energy 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

2451:1040 and the time average stored energy is 8.905e'18.

83

-lO

x10

 

...............................................................................
.....

 

 

 

...............................................

...............................................................................................

      

    

 

 

   

 

 

o ; i 4 i 9' 1' i ;
3.39 3.391 3.392 3.393 3.394 3.395 3.396 3.397 3.398 3.399 3.4
Time Steps x 104

Figure 2.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
with cross is the calculated power data. The average input power is

2451:1040.

84

 

 

 

 

 

x10-
‘ T F .' Y '
2r- . ....
o S
' %
-2 ...
8 i
C
3 5
.35: t
_4 ............................................ ............ ..
-3 _
-8 ...... ..... ... ...
40 4 A ‘1 1‘
0.5 1 1.5 2.5 3.5
TimeSteps x10‘

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

 

-0.01

-0.02

  
 

.5
0
09
T

Ey ( Relative Unit)
.4:
o
A
l

—0.05

l

-0.06 -

 

 

 

_0.07 l l l l l l l
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

x (m)

Figure 2.37 Field distributions of Ey along x axis at 0.14868 us. e',(u)) is
equal to 2.5 and 8",(m) is 0.1.

86

 

 

 

 

 

0.3 L 1 i i i i
O 0.01 0.02 0.03 0.04 0.05 0.06 0.07
X(m)

Figure 2.38 The ratio of calculated and fitting Ey along x axis at
0.14868 us.

87

CHAPTER 3
SOLVING MAXWELL’S EQUATIONS BY BODY

OF REVOLUTION FDTD

In the study of the interaction of microwave radiation with non-ionic materials, a mate-
rial sample is placed in a cylindrical EM cavity which is excited with a fundamental cavity
mode. Most of the time, objects that are symmetric about an axis are encountered; hence
the body of revolution (BOR) FDTD is examined and used to solve the problems involv-
ing cylindrical cavity loaded with symmetric material sample.

The cylindrical cavities with perfectly electrically conducting (PEC) boundaries are
first studied in this chapter. The BOR FDTD formulation of Maxwell’s equations in cylin-
drical coordinates is first constructed in section 3.1. Two sets of equations are derived for
the FDTD formulation. The selection of which set of equations to be used for field calcu-
lation in cavity problem is determined by the characteristic of cavity modes. Then the sin-
gularity problems at p = O is discussed. The Blackman-Harris (BH) function is used to
construct the excitation source in the cylindrical PEC cavity calculation since the sidelobe
0f BH function is approximately -92 dB. The excitation probe is assumed to be located at
a P0th on the cylindrical wall. In section 3.2, the surface impedance boundary condition

(SIBC) is used to simulate the finitely electrically conducting (FEC) boundary. Three

88

methods of the FDTD formulation for SIBC is presented and the frequency domain
approximation is used in the actual calculation since it gives stable numerical solutions.
For the cavity with lossy wall, a continuous source is used and the time step to reach the

steady state is studied. TM012 and T131“ modes of the empty cylindrical PEC and FEC

cavities are calculated and the results of the former case is compared to the analytical solu-

tion in section 3.3. The field distributions of TM012 and ’I‘Elll modes with a cylindrical

PEC cavity loaded by a small cylindrical material centered in p = 0 are then calculated.
The results of the former case are compared with corresponding theoretical estimation.
The field distributions of the cavities with a thin rod or a thin disk material are also calcu-
lated and compared to theoretical estimations.

After the induced electric field inside the material sample is accurately quantified, the
dissipated power density inside the material sample can be calculated. This dissipated

power density acts as the heating source to raise the temperature of the material sample.

3.1 The BOR Formulation of Maxwell’s Equations[28]

Consider Maxwell’s curl equations in non-dispersive medium written in cylindrical

coordinates in a linear material,

Vxfi = [e]%+ {013+}, (3.1)

VXE — “Heal? + [6*]7‘; (3.2)

where [8], [0‘], [u], and [0*] are electric permittivity, electric conductivity,

0 t n n e e o A 0
permeability, and magnetic conductrvrty m tensor form, respectively. J, rs the known

89

impressed current source. These constitutive parameters can be expressed in tensor form

in cylindrical coordinates as

“on am 0‘92
“w “w “w
_azp (12¢ azzj

[0!] = (3.3)

 

 

where or represents the relative permittivity, 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 of field is expressed as a Fourier series,

E = 2 (Eucosmo + Evsinmo) (3.4)
m = 0

I? = 2 (zucosmo + zvsinmo) (3.5)
m = 0

73 = 2 (;ucosm¢ + jvsinmo) (3.6)
m = 0

Where m is the mode number and is an integer because of single valueness of azimuthal
x x A 5 ¥ ¥

dependence of the field. Fourier coefficients eu , ev, hu , hp, 1“, and jv are dependent on

r. z, t, and m where u and v stand for the even and odd azimuthal dependence,

reSpectively. The fields, E and f1 , in (3.4) and (3.5) describe the fields at any point in the
6“tire space of interest because the objects considered are symmetric about the z-axis. If

the objects considered are not symmetric about the z-axis, then different pairs of model

90

expansion are needed for different regions, followed by matching boundary conditions
between those different regions. It is not suitable to use BOR scheme for this non-

symmetrical problem. Note that if the loaded material is lossy, then Maxwell’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:
"1“ 3 " aEu v a .
:64) x hv, u + Vxhu, v = [e] a" + [o]eu’ v + 1“,», (3.7)
i156) x 211,11 + Vxéu, v = 4MB]; v + [0*]I'iu, v. (3.8)

The above two vector equations can be separated into two independent groups of six scalar
equations. These groups represent modes that are azimuthally perpendicular to each other.

The first group is as follow:

 

 

 

 

 

 

3e v v .
Be a ah u 3h .
[8]—ai;—+[c]e¢’u = a: — 33“"th (3.10)
36 v 13(9h ) .
[e] a? +[G]ez,v = "TX—v“? “—12,, (3.11)
3h u 3e
[11] a? +[o*]hp’u = 3‘2“?“ (3.12)
ah¢v aerv aezv
[p] at, +[o*]h¢,v = — 32’ + ap’ (3.13)

91

ahz u 13(9341 u) m
' * -_______;__ _.

 

[it] (3.14)
The second group can be obtained by a similar procedure. Note that the signs of terms that

contain {3 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 BOR representation of Maxwell’s equations

 

 

 

 

 

 

 

 

Fields of first group Fields of second group
ep v ep’ u
3c) u 311 v
62 v ez, u
hpu hmv
ht). v hit. a
hau hav

 

 

 

 

Assume that those constitutive parameters have the biaxial tensor form in cylindrical

coordinates given by

" 'i
0190 O
[a]: 0 09,0 . (3.15)
0 001

Z

 

 

92

Then these two groups of equations can be represented in matrix form[27] as

 

 

 

 

 

 

 

 

 

 

m )- -1 i- T
0 "az :5 ep (poupa,+o*p)hp
dz 0 —3p 6,» = - (pop at+o*¢)h (3.16)
¢ ¢
m 13 ez (l1 a + at h
-— —— _ _ o“ O )
:Fp papp 0‘ _ z t Z 24
m —- — _
O “a: T- hpi (eoepa,+op)ep+ jp
im 1.8.9 0 _hz_ _(eoeza,+oz)ez+jz_
L p pap _

 

 

where 3a denotes 3% and or = p, z, t.

3.1.1 Mode Selection in BOR Algorithm

There are two important issues that need to be determined in applying BOR algorithm.
The first one is which group of 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 of modes and the group of equations to be used. Consider the incident filed[25],
E‘ = f(p)E(t—§)mcos¢—isin¢1 = bE‘p+i>E"¢ (3.18)

where c is the speed of light and p, (b , and z are the cylindrical coordinates. The Eip is

even functions with respect to 11) so is H '1), by Maxwell’s equations. Similarly, Ei¢, H 1p,

93

and H 12 are odd functions w.r.t d). The total fields denoted by a superscript t will preserve

this angular dependence[l9] because of the symmetry about the z-axis of the object.
Therefore, these total fields components can be expanded in terms of a Fourier cosine or
sine series. The Etp, H11), and E2, are even functions and H'p, Et¢, and th are odd
functions w.r.t (1); hence the second group of 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 (4)0, zo) . The source can be expressed as

7= —pf(p)5(z— zo)—— 8(¢p—;L°)gt (t) (3.19)

where f (p) and g(t) are variation of 3 along p and t, respectively. Symbol 6 is Dirac

delta function. If we let $0 to be zero then the delta function 8(4)) can be expanded by
Fourier cosine series as:

co

8
8(4)) = Z —;)-['—"cosm¢ (3.20)
m=0
where
1 = 0
30m : m , (3.21)
2 m¢0

Form theory, the lowest TE mode is TE011 and its p component, Ep , is equal to zero.

Hence, the o dependence of E¢ is cosmq) in order to keep E4, from becoming zero for

94

TED“ mode[26]. Therefore, the first group of equations need to be used in this FDTD

calculation for cavity problems and (3.19) becomes

3 = p Z jmcosmq) (3.22)
m=0
where
E m
jm = ——:r—¥8(z—zo)g(t). (3.23)

3.1.2 The BOR-FDTD Formulation[28]

In this section, the BOR-FDTD formulation will be derived by using the leapfrog

scheme with 0* = 0. Consider the first equation in (3.17),

3h
_ . 111 m
Sp-a—t — -O'pep-]p-a—Z:FE z (3.24)

where the right hand size will be evaluated at t = (n + 9A: since h field is evaluated at

t = (n + am. The standard FDTD notation will be used in this derivation, p = iAp

and z = jAz. Using field locations in Figure 3.1, (3.24) becomes

n+1/2 . 1/2

8 1 . . . . . .
Kiley (I’D-9:011” = —opep (tn-13’

_h$+'/2(i,j)-h3H/2(i1j-1); m
91 9141/2

(1.1)
(3.25)
h: + “20'. j)

 

The components, e; + 1”U, j) and jg” 1”(1, j), have to be approximated by those values

adjacent to them since they are sampled at integer time step n. The second-order

95

approximation used here is

+1. . . . —1 . .
fg+“2(1,j) = 3f; (1,1)+6f;,'8(z.J)—f2 (1,1)

 

(3.26)

where f is e or j. The FDTD representation of equation (3.24) is then

n+1.

. . . —l . .
6,, (H) = Ap-e;(t,1)+Bp-ef, (1,1)

+1/2- -_ +1/2- -_ 27)
kg (1.1) kg (1.1 1)+ m hg+l/2(i’j)T

 

n+§
—C ' j (irj)+ —
p [0 AZ 9141/2

where

 

A = —— (3.28)

 

 

= -— (3.29)

 

 

C = —— (3.30)

 

and x = p, q), or 2. Note that jg+ ”2 is not replaced by the approximation of (3.26) in
(3.27).

The other five FDTD equations can be derived by using the similar procedure and are

listed below;

96

n+1

.. . . -1..
e¢ (‘91) = A¢-e$(r,})+B¢-e; (12])

n+-— hn+l/2(i ')_hn+l/2(i_1 )
e 2 o o ,J ,J
—C¢ - [1, (t. J) + Z A; (3.31)

 

 

_hg+ l/2(ia j) — h3+1/2(i9 j- 1)]
Az

n+l .

. .. —l ..
eZ (1,1) = Az-e:(1,})+Bz-e: (1,1)

 

 

 

 

n+-
. 2 . . m . .
—Cz°l:]z (1,1)3Fahgfl/2031) (3.32)
1
91+ 1/2h&:+ l/2(i11)‘ 91—1/2h$+1/2(i-11j)]
piAp
. . . . m . . €"(i,j+1)-e"(i,j)
hg+1/2(1,)) = hg‘1/2(r,j)IFGpaie:(z,j)+Gp[ 4’ AZ 3 ] (3.33)
. . . . e"(i,j+1)-e"(i1j) e"(i+1,J')-e"(i,j)
h$+”2(t,1)= h3‘1/2(z,1)-G¢[ " AZ " — Z AD 2 ] (3.34)

m

 

 

hg+‘/2(i.j) = 122-“20.1):Gzp

pg+1eg(i+ 1117- 0183(1) f)
(3.35)
91+ 1/2Ap

ego; j)—Gz[

1'+1/2

where Gx = At/(popx), x = p, (b, or z and pi = (i—l/2)Ap and 91/2 = p0 = 0,

and the fields associated with the coordinate (i, j) are shown in Figure 3.1. In Chen’s paper

[28], those BOR FDTD equations are obtained by using a first-order approximation,

n+1 . . . .
z, + 1,
fg+l/2(i,j) = fp ( J)2 f;( 1), (3.36)

 

97

 

I \/ l \/ |
I /\ T /\ I time
n- 112 n n+1/2 n+1 n+3/2

 

 

 

 

 

 

 

FIGURE 3.1 The field locations for BOR FDTD in time and
space.

98

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 91/2 factor, which is equal to zero, in the
denominator of the fourth term when calculating the ep(0, j). Hence, ep(0, j) is infinite
at p = O 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, hp, 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 FDTD equations

can be used for time stepping.

As observed from (3.31) and (3.32), only the components h¢(0, j) and hz(0, j) are
needed to update the adjacent ez(l, j) and e¢(1, j) fields internal to the mesh. The
60(0, j) and ep(0,j+ l) are needed to evaluate h¢(0, j). The field, h¢(0, j), is not

needed actually to calculate the ez( l, j) since

n+1

. . -l .
ez (1,))=Az.e;'(1,))+Bz-e;‘ (1,1)

n+-
. 2 . m -
—Cz.[1, (1,1)iah3*‘/2(1,1) (3.37)
l

 

93/2h$+1/2(11 J) ’ 91/2111? + V2“): 1)]
plAp

With 91/2 = 0; hence, ep(0, j) is not needed, neither. In actual calculation, h¢(0, j) and

99

ep(0, j) are set to zeros to avoid cumulation error.
Now, only h Z(O, j) is needed to update the fields point in the FDTD lattice. Note that

hZ is zero at p = 0 for m at 0 by (3.17); hence, we only need to evaluate hz(0, j) for

m = O. From the integral form of Maxwell’s equation in the time domain, we can obtain

the following time update equation for hz(0, j) when m = 0

4A!

 

hg+ 1”(0, j) = hg-1/2(0,j)— e;(1,j) . (3.38)

Once hz(0, j) is known the rest of the field components can be evaluated using (3.27) and

(3.31) to (3.34).

3.2 Surface Impedance Boundary Condition

When the boundary of a cavity has a finite electrical conductivity (FEC), there are two
FDTD approaches to calculate the field distributions in the EM cavity. For a cavity with
good conductor wall which is considered in this chapter, the skin depth and the local
wavelength are very small compared to the radius of curvature of the cavity wall. Hence,
the planar surface impedance boundary condition (PSIBC) is used to approximate the
lossy conductor wall. For more accurate simulation, an absorption boundary condition
(ABC) 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 frequency-domain concept. Consider a time-

harrnonic plane wave incident at angle 6,. on the lossy dielectric half-space as shown in

100

Figure 3.3.

The surface impedance for a planar interface at z = O can be defined through the

following relation:

Em) = 2(w)laxfi.(w)1

(3.39)

where Z ((1)) is the surface impedance. The surface impedances of interface for TE and

TM plane wave[30], respectively, are

 

2mm) = 22110392

where

Z jwuk 1/2
" _ (jwsk + 0k)

2 1/2
sin 91'

e 1— j“
””( (11808,)

« 1 , then the surface impedance can be rewritten as

 

cost)2 =11—

 

 

When lsin 202

101

= (1 - sin202)“2.

(3.40)

(3.41)

(3.42)

(3.43)

 

 

 

 

 

 

 

 

 

 

/ I
/
/ ______ r / I
/ / /

/ / /| I
" “ ‘ , 'f I I
' ’ | I I

r ——————— 1
' I | I I
I I I | |
' I ' I I
' I l ) I
/
I I /I J
l /
I /
I L _____ I /
/
l I /
I... _______ _ _ _ ._ J/
Inner Empty Cavity FEC Wall
, —— — —. ABC
/ // [VIII"”’ \ \ /\
/\ ’1’; /
I — I
| l
| l
| |
I _' __ __ ~ I
l / ’ [III], ‘ ‘ \ |
/ yyl " I \ I
1 é ‘1
\ ’/’ / /
’llllll/ ,

~_"

Figure 3.2 FDTD configuration for cavities with FEC wall

102

Free Space ‘ X Lossy Dielectric

 

 

 

 

 

 

 

Figure 3.3 Coordinates for the incident and reflected plane
waves upon a lossy surface

Z(u))

jwuz )1/2

45“”) = 2mm) = (m

(ll—21)1/2 S
710 82, Si-O'z/Sz

(3.44)

Where ‘10 = /110/80 and s = jw. Equation (3.44) is the planar surface impedance used

through this chapter.

32.2 FDTD Implementation of Planar Surface Impedance

In order to do FDTD calculations, the time domain expression of (3.39) needs to be

Obtained. The convolution in time domain is involved since there is a multiplication in

frequency domain. The convolution in time domain needs to be approximated by some

103

recursive expression in order to run on a computer. Three approximations are discussed in
this section. The time domain approximation[29] is first considered; however, this
approximation encounters some divergence problems when evaluating an integral and is
not used to the FDTD implementation. The frequency—domain[31] and Z-transform

approximation[32][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 time-domain surface impedance is

“2r 1/2
z(t) = 110(5) {5(t)+ae‘"[Il(at)+Io(at)]U(t)} (3.45)

where U(t) is the unit-step function, a E—oz/(Zez) and I0 and I 1 are the modified
Bessel functions of the first kind of zero and one order, respectively. Hence, the time

domain expression of (3.39) is

A “2r 1/2 A _‘ A _‘
Ez(t) = n0(-8;-) [an,(t)]+foae‘"[11(at)+Io(at)][an,(t—'c)]d1: (3.46)

and the discretization version of (3.46) is

_, “2, 1/2 _. ""1 _.
E1(m) = 110(5) {U2 x H.(m>1+ 23 Foam x H1(m — cm} (3.47)
r q = 0
Where
Fem) = f: +1’11 —|v - qII<aAt>e<a~>'II,(aArv) + IoraAmIdI (3.48)
and

104

q) = { . (3.49)
q
The Prony’s method is used to approximate the F 0(q) and (3.47) is then

A

E1

 

Q
n _§ ’1 ..x
= n2[H,| + 2 Gk(n)] (3.50)
k=l

.A ..x n ..x
where 01(11) = CkH,| +tthk(n—1)andk = 1,2,3, ...Q.

3.2.2.2 Frequency Domain Approximation

The surface impedance, (3.44), is rewritten as

2(5) 112,/S—-:a (3.51)

where s = jo), 112 = 710 ’2—2', and a = (32/32 . Equation (3.51) can be approximated by
2r

a rational polynomials in frequency domain and is rewritten as

Z(s') = l s. zl—i C’ (352)
1+s' . l(1)i+s' '
1:

where s' = s/ a, the L is number of terms needed in that approximation and (1),- are

 

 

known poles. The L and (1) 1' 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 will accommodate most material up to several tens of Gigahertz. The (3.39) in s-

plane is then

105

L

.3 aC- ...:
E,(s)zn2[l — 2 am._;s:l[fiXH,(s)]. (3.53)

l

 

l:

Assuming the waves are piecewise linear in time, (3.39) in discrete-time domain is

 

given by
_\ n A _; n L A n
E,| = n2[n><H1I ]— 2 Ail (3.54)
i=1
where
.9 n A _‘ n A _~. n—l _~. n-l
A,- = pil[an,' ]+pi2[nXH;| ]+p,3A,| (3.55)
and
C- (1),-At
p,l = 1120—)f[1+(e*' —1)/(Ata(1)i)] (3.56)
1
Ci —a(1),-At
1712 = nzall/(Atami)—e (1+ 1/(Atawi))] (3.57)
l
—a(1),-At
pl.3 : e , (3.58)

3.2.2.3 Z Transform and Digital Filters Approach
The term, Ci/ ((1),- + s'), in (3.52) is a analog lowpass filter with a single pole. The

digital version of this term is a digital IIR filter. Hence, this become an IIR 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

Table 3.2 Rational Approximation Results

 

Number of
Terms

c.

l

 

 

2.60266906380e-10

3.31195847511e-8
1.58725612150e—6
4.35701245008e-5
8.0865763 8450e-4
1.071972546526-2
9.48589314718e-2
0.3933263026213

2.43632801 126e-7
1 .06720343642e-5
1.59612026325e-4
1.59383673429e—3
1.20627876345e-2
7.22902130819e-2
0.3267172939002
0.86807551 10248

 

s-plane to z-plane and can be linked to the trapezoidal formula for numerical

integration[34]. The substitution of variable for bilinear transformation is

 

2(1—21-1)
s = — 3.59
T 1+z‘1 ( )

where T is the period of one time step or the sampling factor. Applying the bilinear

transformation, (3.59), to (3.53), then the time discrete version of (3.39) becomes

 

L
.5 n ..A n A n
E,| = n2|:fixH,| ]— Z f.| (3.60)
i=1
where
..x n 2"awiT Aln-l aCinzT A I? n A I? n—l 36
f' _(2+a(1)iT)f' +(2+a00iT)[nx 'l +nx '1 ] ('1)

The bilinear transformation maps the entire left half of 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 will take place when

107

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. By (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 time step.
However, the tangential magnetic field is evaluated at half-time step before the current
time step and half cell in front of the cavities wall according to Yee’s FDTD lattice in
Figure 3.1. In the practical FDTD 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.60) can be rewritten as

 

_. n _, n—l/2 L _, n
E,| = n2[an,|h ]— 2A, (3.62)
i=1
where
A n 1. _‘ n-l/2 1. .3 n—3/2 A n-l
Ail = Pil["XH'|h ]+pi2[nXH,|h :|+pi3A'l (3.63)
and
L
A n _‘ n-l/2 A n
Etl = n2[fiXH(|h ]— 2 it (3.64)
1:]
where

_. n 2—a00iT _. n-l aszT . _, n-l/2 . _, n-3/2
fi| = (m)fi| +(——)[an,| +nXHtlh ]. (3.65)
1'

108

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 FDTD

implementation of Maxwell’s equations.

3.3 An Empty Cylindrical Cavity

The physical configuration of the cylindrical cavity is shown in Figure 3.4. The radius
and height of the cavity are a and h and that of loaded material are b and l. The wall of the
cylindrical cavity is assumed to be perfectly electrically conductive (PEC) or finitely
electrically conductive (FEC) and the excitation probe is located on the side of the
cylindrical cavity.

The TM012 mode of an empty cylindrical cavity with PEC boundary is first calculated

to verify the program. The g(t) in (3.19) is a Blackman-Harris (BH) windows function

with the central frequency of 2.4571GHz and 0.1GHz bandwidth. The BH function makes

the source turn off automatically and the 0le bandwidth is chosen because only

TM012 is desired. The dimensions of the cavity are a = 0.0762m and l = 0.15458m and
the dimensions of grids are Ap = 0.0006m and A2 = 0.00118 which are lees than
A/ 20 where A is the wavelength. The numbers of partition along p and z are 127 and
131, respectively, and the one time step is 1.69363 ps. The field distributions of Ep and

E2 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 of

109

TM012 mode for the FEC and the loaded cavity. For TE 1 1 1 mode, the numbers of partition

are 127 and 67 along p and 2 directions and the one time step is 1.29676 ps. The field

distributions of E p, E q) and E on the p — z plane for TEl 11 mode are plotted in Figure 3.8

to Figure 3.10. For TEm 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 for g(t) in source equation (3.19) since

there is a loss on the boundary. As observed in Figure 3.11 with At = 1.2967608 ps , the E
field achieves a steady state about 23,135 time steps for o = 102 S / m and 120,000 time
steps for o = 104 S/ m. It will need more than 600,000 time steps to evaluate the field

distribution of TEN] mode for cavity with o = 106 S/ m which is usually encountered

in regular cavities.

The field distributions of E2 and E on the p-z plane are shown in Figure 3.12 and
Figure 3.13. There is some Ez near the boundary wall form Figure 3.12; however, the
value of E2 is much smaller than that of ED. The B 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 FDTD calculation for field distributions in cavities since most cavity

wall is made of metal with conductivity larger than 0 = 106 S / m. Moreover, this PEC

assumption requires much less time steps in FDTD simulation.

110

 

 

 

 

 

|<1
PEC
or I T
F C I
2 I
”’1,IT\-.\ _
i" 31.3%., h
I gr“ I l
I ,, : ,I
\\_L-/ —
I
1:21 '
I
' L
I
I

Figure 3.4 The physical configuration of a cylindrical cavity loaded

with a material sample.

 

 

 

 

I I I I I T I

0.1» - 9 ° 3 —
m = 0
Z =0.0531
0.08 ’- ..
0.06 ~ 4
p004 *- _
0.02 _ 1
m = 1 m = 2

o—--—— —- ———-—— — ————— ——-——— ——-——— — ————-— -— ———

_002 1 l J I l l I

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07

p (m)

 

 

 

 

 

z (m)

Figure 3.5 The variation of E9 of TM012 along the p and z directions in
a PEC empty cylindrical cavity.

112

 

0.05 -

  

z =0.0295

  
   
   
 

 

 

 

 

 

 

m =
o)..____. _._____ _ __—_ __——_— _—__ _ ——___ — _—
E
Z
2 =0.0531
-0.05 r "
m = O
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)
0.1 m = 0 p =0.0372
0.05
E2 0
-0.05
-0.1

 

 

 

0 0.05 0.1 0.15
2 (m)

Figure 3.6 The variation of EZ of TM012 mode along the p and 2
directions in a PEC empty cylindrical cavity.

113

 

_

c.

I

\

\

I

\

\I ’\ ’\ '\
\ \ \ \ \ \ K ._
\\\\\\\-~—_¢—
\\\\\k&w

'\
'\

I

7V7)

\\
\\
\

r\\\2
T\\\2
TA\x\\\

TTT‘ \

 

_ .22 22

22“ .22“ M22 -m
22,..22 22 ,
’2/1\\\\ d,//nm
’//1\\\\ ”TI/1
22,1.\\2 222 -m
222/1\\\2 22/HM
(///Iv\\\\‘ ///ft.mm
K ///vrw\w\\\:///f1 0.
//N/vlv\v\\\‘d///7th
: 2313i: 221/31
NH/¢.Iv\v\w\\\2//A/A/Ahm
///v/vlv\v\v\\\2//A/A/Al

/ J/vlvlv\v\v\w\\ 2.21/Tim
3¢IvleY1vthx.,I/TAIAI

 

 

0.15 *-

_
1. 0
0

0.05 r

\.I
m

(
2

Figure 3.7 Plot of E field of TM012 mode on the p-z plane in an

empty PEC cavity.

114

 

 

 

 

 

 

p
0.5 ’- -4
l l l l l l J
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)

 

 

 

 

 

l
O 0.01 0.02 0.03 0.04 0.05 0.06
z (m)

Figure 3.8 The variation of Ep of TEln mode along the p and 2
directions in a PEC empty cylindrical cavity.

115

 

 

 

 

1.5 T r 1 fl r fl 1
m = 1
Z =0.033024
1 .- -
Z =0.04644
E
9
0.5 " .1
6 _-__. —m;2— -._ ._.__. _._ ___ _._-__. -.- _._m.=o__.
1 1 1 1 1 r 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)
1 4 r r I f r r

 

 

 

 

 

0 0.01 0.02 0.03 0.04 0.05 0.06
2 (m)

Figure 3.9 The variation of E, of TE111 mode along the p and z
directions in a PEC empty cylindrical cavity.

116

0.06

0.05

0.04

z (m)

0.03

0.02

0.01

 

l

 

—-> —-b —p —> ——> —-> —-D —-. —o —. ... —. —o —. —.
—>-—>—->—>—>—>—a>—->—>—»—>—+ 2., —. —.
—>—>—>—>—>——>—>——>-—>—>——>—->——p —>
-—>——>——>—-—>—>—->——>—>—>—>-—>-—>—>——>
—->——:>—:>—>9——>——:>——>—>—->——>—>—>—>
——>—>9—:»—>—>——>—9——>—>——>——>—->——>—->
——->—>—->—>-—>——>——>——>———>——>——>——>—>—>—>
——>——-—>—9——->——>—>——>—>—>—>—>——>—->-—->——>4
—>—+——>—>—>—>—>—>——>-—>——>——>—>—>—>
——>——>——>——>—>—>—>—>—>—>——>—->—>-—> __,‘

——>——>-—>——>—->—->—>—>—>—->—>—>—>—-> _.

 

 

—-D -—-D -—-D —> —D —D —’ --> —D —-> —. -—b -D —. —o
l l L l l l 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)

Figure 3.10 Plot of E field of TEul mode on the p-z plane in an
empty PEC cavity.

117

3.4 A Loaded Cylindrical Cavity

For a loaded cavity with PEC boundary, three cases for TM012 mode and one case for

TEm mode are calculated. The first one is a cavity loaded with a small cylindrical
material, with I z 2b, centered at p = 0 . The second one is a cavity loaded with a thin

rod material with l » 2b and a cavity loaded with a thin disk material with l « 2b is the
last case. For TE“ 1 mode, only the case of a cavity loaded with a small cylinder material

is studied.

3.4.1 Small cylindrical sample for TM012 mode

A material sample, having the diameter equal to the height, is placed in the center of

the cylindrical cavity. The diameter is 0.006 meter and the height is 0.0059
meter; that is lz2b. The induced electric field inside the material sample can be

estimated by the electrostatic field induced inside of a dielectric sphere as
1'
Ez = (3/(2+e,))Ez.

The numbers of partition used along the p and 2 directions in this FDTD 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 of the electric is significant
near the center of the empty cavity. Due to the small dimensions of the material sample, E2

is still the dominant component near the center of the loaded cavity if other components

are compared with Ez in Figure 3.14 and Figure 3.15.
The ratio of E2 component of the induced field inside the material sample to that of the

induced field in an empty cavity is plotted in Figure 3.16. The ratio is between 0.6 to 0.75

118

and it depends on the position of p and the fitting sine function. The electrostatic
estimation of that ratio is 0.667 with 8, = 2.5. The numerical results and the theoretical

estimation are in satisfactory agreement.

3.4.2 Thin rod case for TM012 mode

A material sample with the shape of a thin rod, having its height much larger than its
diameter, is placed in the center of the cylindrical cavity. The dimensions of the material
sample are l = 0.13098 meter and b = 0.0018 meter. The ratio of the diameter to the
wavelength is about 0.03. The numbers of the partition along the p and 2 directions and

one period of time step are the same as those considered in section 3.4.1. The dominant

component near the center of the cavity is still Ez if

ED and E2 components are compared in Figure 3.17 and Figure 3.18. The field distribution

of E2 along the p axis seems to be not affected at all in Figure 3.18. In fact, E2 is not

affected by the placing of sample or the effect is too small to be noticeable. From Figure
3.19, we observe that the ratio of the z-component of 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 of that ratio is 1 based on the continuity of the tangential
component of the electric field on the sample surface. This electrostatic estimation agrees

with our numerical results.

3.4.3 Thin diSk case for TM012 mode

A material sample with the shape of 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

 

 

 

 

 
  
   

J I I I I I l I
m=1
2» o=10e2S/
I I ‘I‘ ‘ I,“'II,)1;‘1‘(:I
I IIII‘1|)I‘ ‘1 1 1 II I.)
1» “NMIIMI II I I I.)
.‘I Il ‘
E
00
II I
-1- i I , ,
‘ 1 I I I
-2—
_3 1 1 1 1 1 1 1
o 0.5 1 5 2 25 3 3.5
Time(sec) 1110*
20 I I I I I I I
15~ m=1
a=10e4S/m
10
5
E
po
-5
—10
—15

an
’su

 

I 1 1 1 2L 1 1

 

0.8 1
Time (sec) 1110'7

Figure 3.11 Plot of Ep of TElu mode versus time in an empty FEC

cavity with 6: 102 S/m and 6: 102 S/m with m = 1.

120

sample are l = 0.00118 meter and b = 0.006 meter and the ratio of diameter to

wavelength is 0.1. The z component of the induced electric field inside the sample is still

the dominant component. The induced electric field of this thin disk geometry can be
estimated theoretically by the boundary condition of E z = (1/8,)E:.
The numbers of 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 of E z/ E; is

about 0.42 which is consistent with the theoretic estimation.

3.4.4 Small cylindrical sample for TE111 mode

The dimensions of a small cylindrical material are l = 0.00516 m and 2b = 0.0048

m. The numbers of partition along the p and z and one period time step are the same as

those used in calculating the TElll mode in an empty cavity. The numerical results of

field distributions are shown form Figure 3.23 to Figure 3.25. There are some Ez near the
center of the cylindrical cavity due to the placing of material sample as observed in Figure

3.25; however, the magnitude of E2 is much smaller than Ep or E¢. When 0 = 11/2 , E4, is

zero since the first group of equations in Table 3.1 is used. Hence, the E p is the dominant

filed distribution near the center of the cavity. From Figure 3.23, we observe that the ratio

of 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 ,9 at 111 = 0.

121

 

0.8 -

0.6 -

z =0.031992

 

 

 

 

 

_02 1 1 1 1 1 1 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07

p (m)

Figure 3.12 Plot of E2 of TEnl along p direction in an empty FEC cavity
with conductivity 104 (S/m).

122

 

—. —. —. -. -. -. -. —.
—o —o —o —o —o —. —o —o

0.06 '- —. -—o —. —. —o —o —. —o
—o —o —. —o —-o —. —> —-o
-—. —> —-p —-D —. —p —D —.
—5 —> —> -—> —> —> —> —-p

0.05

0.04
z (m)

—>-——>——>—>—>-—>—>——>
——>—-—>-—>——>——->—>—>—>
-—>-—>——>—>—>-—>—>——>
—->—>—>—->-——>

0.03

M
iii

HI
HIHHHHHHH

""”HHIHHHIHHIHW*"'
-""‘HHIHHHiHHHIIIHII"°_

0.02

IH
111

0.01

 

‘

Jilllliiiilllliltltsa

I

"‘tHIIHHHHHHHHH“""

_""”HHHHHHHHHH“""'_

.itilll]]]ll£l

 

 

0.06 0.07

Figure 3.13 Plot of E field of TEul mode on the p-z plane in an empty

FEC cavity with conductivity 104 (S/m).

123

 

0.04 r r T I I r I

0.02 ~ .
z =0.07434

 

0: 2b=10Ap ~

 

=5Az

 

 

 

 

 

 

0.01 0.02 0.03 0.04 0.05 0.06 0.07

 

0.2 I T I

0.151-

   

p =0.0012

  

-0.05 r

 

 

 

o 0.05 0.1 0.15
Z (m)

Figure 3.14 Plot of Ep of TM012 mode in a PEC cavity loaded with a
small cylindrical material sample.

124

 

0.9 T r 17 I ‘ I I

0.8

0.7

0.6

    

z =0.07906

0.3

0.1%

 

 

 

 

 

 

 

 

 

0.05 0.1 0.15
2 (m)

Figure 3.15 Plot of E2 of TM012 mode in a PEC cavity loaded with a
small cylindrical material sample.

125

Ratio Of E /Ei

1.4

1.3

Zr‘

2

.0
to

9
a:

0.7

0.6

0.5

 

 

 

 

I I I I I
2 I \ 7
’1
/
"’ /
/ \
/ v
\
\ I
- \ p =0.0024 , -
\ / " ’ I.- a ‘ ‘ ‘ I
\ / I \ \ I
— \. / \ I .-
\
\l
- p =0.0012 -
l l l l l
07 0.072 0.074 0.076 0.078 0.08
2 (m)

Figure 3.16 Plot of the ratio of E Z/ E 'Z of TM012 mode along the z

direction in the PEC cavity loaded with a small cylindrical material
sample.

126

0.082

 

004 I I f I I I I

0.02 '- _.

 

 

 

-0.08

z =0.07434

~01 1 1 1 21 1 1 1

 

 

 

0.01 0.02 0.03 0.04 0.05 0.06 0.07
IO (N!)

0.15 TI 1 r
m = 0

 

0.1 ~ .

   
 

p =0.0012

p =0.0024

 

 

 

 

0 0.05 0.1 0.15
2 (m)

Figure 3.17 Plot of ED of TM012 mode in a PEC cavity loaded with
a thin rod material sample.

127

 

 

 

 

 

0.8 I I I I r I W
m = 0
0.7 )- 2
z =0.07906
0.6 )- 1 .1" . -I
z =0.07434
0.5 *-
E...
0.3 1-
0.2 —
0.1 -
o 1 l 1 I l l l
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)
1 I f T

 

 

 

 

 

0.05 0.1 0.15
2 (m)

Figure 3.18 Plot of E2 of TM012 mode in a PEC cavity loaded with
a thin rod material sample.

128

1.02

 

 

.......................................................................................

(O
00
I
.0
ll
.0
N
.h

229

Ratio of E IEi
28

 

 

 

0.94 ~ 2
0.92 ~ 2
09 1 1 l l 1

0.07 0.072 0.074 0.076 0.078 0.08 0.082
2 (m)

Figure 3.19 The ratio of 132/E: of TM012 mode in a PEC cavity
loaded with a thin rod material sample.

129

 

0.045 r I I I r I r

 

 

   
    

 

 

 

 

    

 

 

 

0.04~ m=0 .
0035— .
0.03» 2b=20Ap 1
l 1 =2Az
0025— if -
p
0.02- -
0.015» 2
0.01~ z=0.0767 i
9005" , ‘
.. ..... 2:0.0772
oo 0 01 0.02 0 03 0.04 0 05 0.06 o 07
p (m)
0.03 I r I
p=0.0048 m=0
0.02- ‘
0.01 p=0.0036
E
p o
-0.01
-0.02~ -
-0.03 ‘ 1 1
0 0.05 0.1 0.15
2 (m)

Figure 3.20 Plot of Ep of TM012 mode in a PEC cavity loaded with
a thin disk material sample.

130

 

~0.05 *-

   

z =0.07729

    

—0.25 ’-

   

z =0.0767

 

 

 

 

1 1 L 1 1 l l
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)
l I T

0.25

0.2

0.15

0.1

0.05

-0.05

—0.1

—0.15

-0.2

 

-0.25

 

 

 

 

o 0.05 0.1 0.15
Z (m)

Figure 3.21 Plot of E2 of TM012 mode in a PEC cavity loaded with
a thin disk material sample.

131

 

1.2r —

oz 2
to
I

Ratio of E /Ei

.0
\l
I

0.6 r

0.5 -

 

 

 

0.4 1 1 1 1 1 1
0.072 0.074 0.076 0.078 0.08 0.082 0.084 0.086

2 (m)

 

Figure 3.22 The ratio of Ez/E: of TM012 mode in a PEC cavity
loaded with a thin disk material sample.

132

 

    
 

 

 

 

 

 

 

 

 

1.8 I r I r r f T
m = 1
1.5 ‘- r
1.4 h -‘
1.2 2 Z =0.035088 ‘
1 " 2
p I
0.8 '- '4
0.6 " .1
0.4 " -
0.2 — _
0 1 1 41 1 1 1 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)
1 4 I I m I I T
_m31
1.2 - ’ \ J
/
/
/
/
/
1 h / .1
/ \
/ \
/ \
/ \
o.a . , p =0.0012 \ .
E I \\
/
p / \
/ \
0.6 — \ .1
\
\
\
.. \
0.4 \ i
\
\
\
0.2 '- \ _i
\
\
1
0 1 1 1 1 1 1
0 0.01 0.02 0.03 0.04 0.05 0.06
2 (m)

Figure 3.23 Plot of E0 of TEul mode in a PEC cavity loaded with a
small cylindrical material sample.

133

 

E, z =0.035088

0.6 r- “I
0.4 1- ~

0.2 . J

 

 

0 l l l i l 1 l

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07

 

 

0.8 — p =0.0012 .

0.6 -— r

0.2 *- -4

 

 

o l I l l i l

0 0.01 0.02 0.03 0.04 0.05 0.06
z (m)

 

Figure 3.24 Plot of E,» of TElu mode in a PEC cavity loaded with a
small cylindrical material sample.

134

0.035

0.03

0.025

0.02

0.015

0.01

0.005

0.1

0.08

0.06

0.04

0.02

 

 

 

 

 

 

m = 1 ‘
z =0.035088 *
l l 1 l l l J
0.01 0.02 0.03 0.04 0.05 0.06 0.07
p (m)
f T f I I I

 

 

 

 

1 I I l

 

 

l
0.01 0.02 0.03 0.04 0.05 0.06
Z (m)

Figure 3.25 Plot of E2 of TEln mode in a PEC cavity loaded
with a small cylindrical material sample.

135

CHAPTER 4
SOLVING MAXWELL’S EQUATIONS BY FDTD

IN CYLINDRICAL COORDINATES

In chapter 3, we considered the cylindrical geometries with rotational symmetries. The
purpose of 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) FDTD
formulation and applications in cylindrical coordinates. The general 3-D FDTD formation

is considered and discussed in section 4.1. There are two problems in applying cylindrical
coordinate FDTD in 3-D. The first one is the singularities at p = 0. The other
complication that arises in applying cylindrical coordinate FDTD in 3-D is that the cell
size in the o 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 of p = 0 is filled with perfect conductor or otherwise excluded. The two
problems can be solved by using a set of FDTD approximation equations for p = 0 after
a careful examination of the 3-D FDTD Maxwell’s equations. This special FDTD

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 FDTD

136

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
FDTD 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 FDTD approximations at p = 0 can be obtained. Higher-order
FDTD also requires less number of partition than that in second order scheme and hence
reduces the computational time. These higher-order method is studied in section 4.3. At
the end of this chapter, the FDTD formulation of 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 FDTD expressions of Maxwell’s equations are derived. The

differential Maxwell’s equations considered here are

1 013
VxE -—
a:
... (4.1)
VXH = 92+jc+js
a:

where D = [8]E, B = [p]H, Jc = [o]E, and Jc 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

' I
up 0 0
[0t] = 0 01¢ O (4.2)
O 0 01a

 

 

where or represents the electric permittivity, the magnetic permeability or the electric
conductivity. For nonmagnetic dielectric material considered in this chapter, the electric
conductivity of this material is assumed to be zero and the magnetic permeability is
assumed to be the same as that in air.

The scalar Maxwell’s equations in cylindrical coordinate are

83 _ 9.191218%

679 - a. 6% “'3’
g‘” = g—gz—g—fp (4.4)
g?” = ggp—g—LIz—JC¢—Js¢ (4.7)

where Ep, E¢, Ez, Hr, H¢, and Hz are electric and magnetic fields along p, 0, and z,

respectively.
Using the FDTD notations, those finite difference equations for Maxwell’s equations

are obtained as:

138

 

 

 

 

 

 

 

 

 

 

 

     

 

 

 

 

 

 

13"— n-- n-- n—- n__
atDp- 1 P 5¢HZ. 1 . —5zH¢. l . _JC" 1 . —JSP 1 . (4'9)
1+§11k Hid-’1‘ z+§J.k 1+-2-.J.k 5,},k ' —,j,k
n—l n--l- "_l "_1 "-1
2 2 2
8,D¢| = 52H p, —5sz 1 —JC¢ —JS¢ 1 (4,10)
ij+'2',k ij+-,k l]+§,k l,J+-,k 11+5’k
"—1 l "—1 "-1
2 2
5.0.. _ 1 — ‘7— 6,(pH,)l —5¢le. .
"LIME i,j,k+§ 131/(+5 lj,k+-
(4.11)
"-1 "-1
2 2
-ch -Js:
ljrk+" lj,k+-
II n l
5130 . 1 I — 5on 1_p 1 5 . 1 l (4.12)
(91+§,k+2 I]+'k+§ lj+’,k+l ‘j+i’k+§
n n n
Sth, 1 l — 89E, 1 52Ep l 1 (4.13)
+-j.k+- I+-,j,k+- i+-,j,k+-
2 2
n 1 n l n
5th I 1 = 6——8¢Ep 1 1 —p——Sp(pE¢) l (4.14)
' —' — . l . 1 - _- _ . 1 . 1 . _. _
1+2,J+2,k l+§,j+§,k 1+2,]+2,k 1+2’J+2’k 1+2,j+2,k

where 5p(pH¢) and 69(pE¢) need special treatment on p , and it will be discussed later.

The spatial locations of B and E are plotted in Figure 4.1, E has the same location as
B does and so is E and E . The sequence of FDTD calculations along time axis is shown
in Figure 4.2. The D, E, and 813/131 are evaluated at the integer time step; but E , 1:1 , and

BD/Bt are evaluated at the half-integer time step. The order of calculation at different

time steps are also described in Figure 4.2. Note that the boundary conditions are applied

139

after the E is calculated.

Note that the terms 5%(pE¢) and 3%(pH¢) are not factorized out. The first term can

be rewritten as

3 BE,
$(pE¢) = E¢+p§6 . (4.15)

n

The evaluation of E ¢ in the right hand-side is E¢|
1+ l/2,j+1/2,k

which is not where E¢ is

located and, thus, (4.15) is not used.

140

 

 

 

FIGURE 4.1 FDTD lattice for cylindrical coordinates.

141

 

 

 

 

 

 

 

 

I I I I
I I I I time>
_‘n-l _‘n
D J‘i —>D| ——>_\"+-
B 3|
3
11) 1(6)
3"“ "-1 is" +1
_—_‘ 2 An -
HI <5> HI ’
(2) 1(4)
(1) jln-l/Z Eu
3:
_‘n-l/Z
.313
a: ‘

 

 

Figure 4.2 The diagram of order of FDTD calculations along time axis. The
meanings of those steps are listed below.
(1) Using (4.9) to (4.11)

(2) Desired time stepping scheme for B.

(3) Constitutive relation of E and E which depends on material models.
(4) Using (4.12) to (4.14). At this point, the boundary conditions are applied.

(5) Desired time stepping scheme for E.

(6) Constitutive relation of E and E1 which depends on material models.

142

4.2 The Second-order Cylindrical FDTD Scheme[35]

The second-order cylindrical FDTD scheme which applies the Yee algorithm to
cylindrical Maxwell equations is discussed in this section. The cylindrical FDTD
formulation of Maxwell equations is presented in section 4.2.1. The singularity of
cylindrical FDTD equations at p = O and its traditional treatments are discussed in
section 4.2.2. The traditional source FDTD 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 of fields, the (4.9) to (4.14) become the finite difference Maxwell’s equations in

cylindrical coordinates. The finite difference equations are listed below:

 

 

 

 

H n+1/2 H n+1/2
n+1 -E In + At x zi+l/2,j+l/2,k Z|i+l/2,j—1/2,k
pi+l/2,j,k pi+l/2,j,k 8|i+l/2,j,k le/on
(4.16)
n+1/2 n+1/2
¢|i+I/2,j,k+I/2 ¢|i+I/2,j,k—I/2_J n
Az pli+l/2,j,k
n+1/2 n+1/2
n+1 _ n + At x pli,j+l/2,k+1/2 pli,j+l/2,k-1/2
¢i,j+l/2,k ¢li,j+l/2,k 8|i,j+l/2,k Az
(4.17)
n+1/2 H n+1/2
z|i+1/2,j+l/2,k z'i-1/2,j+l/2,k_J n
Az ¢|i,j+I/2,k

143

E n+1 n 2A!

z = EZ +
i,j,k+1/2 i,j,k+1/2 Eliij/z

p H n+1/2 p H n+1/2
1+1/2 ¢|i+1/2,j,k+l/2 "”2 ¢l,~-1/2,,-,k+1/2

 

 

 

 

 

 

 

 

x 2 2
91+1/2‘Pi-1/2
(4.18)
n+1/2 H n+1/2
91+1/2‘pi—1/2 " i,j+1/2,k+l/2 pli,j—1/2,k+l/2
— x
2 2 A4,
Pi+I/2‘pi-1/2
_JZ'" )
i,j,k+1/2
|n+l/2 _ H In-l/Z At
p. . " p. .
1,j+l/2,k+l/2 I,}+1/2,k+l/2 “I1,j+1/2,k+1/2
n n
¢|i,j+l/2,k+1 ¢|1,j+1/2,k_
A2 (4.19)
2'” _Ezln
_ i,j+1,k+l/2 i,j,k+l/2
pi(A¢)
n+1/2 _ H In—l/Z At
(1). . - q). .
1+1/2,},k+l/2 1+1/2,),k+1/2 “li+1/2,j,k+1/2
zln _Ezln
i+1,j,k+1/2 i,j,k+1/2
Ap (4.20)

n

n
p|i+ 1/2, j, k

E
p|i+1/2,j,k+1
Az

 

144

H n+1/2 _ H n-1/2 2A1 1

z . . — 2|. 0 X 2 2

1+1/2,,+1/2,k 141/2,1+1/2.k MRI/214”“ p. 1_p.
’ ! 1+ 1

 

n
p|i+ 1/2,j+ 1,1:

M

n

plI'M/2,131: . (4.21)

 

>< (pi+l_pi)x

—(p,-+ 1E¢li+1,j+ 1/2,k—piE¢|i,j+ V211»
4.2.2 FDTD Calculations at p = 0

There are singularities in (4.18) and (4.19) in the above FDTD equations when p
approaches to zero. Thus, those two equations can not be used in this FDTD calculation.

Appropriate approximations to E ,9 and H p are presented in this section.

D=0 O=0

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

l

—1 At 4 "-2

Ez " = Ez " + — - —H¢ (4.22)
Ojk+1 0jk+1 E '1 1jk+1
’ 9 2 7 ’ 2 2’ ’ 2

where rl is the distance along p for the first cell as shown in Figure 4.1.

Regarding to H 9| , it is used to calculate E¢| in (4.17) and E Zl in (4.18).
p p=0 p=0

However, EZI 0 is approximated by Amper’s law in (4.22). If E¢| 0 is also
p: p:

approximated then the calculation of H p at p = O is not needed. Traditionally, E,» at
p = O is approximated by its analytical solution. For example, E,» at p = 0 of TEl 11

mode is approximated by Ed, at p = 1 times 1.0015 which is the ratio of analytical

145

solution at these two points. This approximation is highly dependent on the mode being

calculated and difficult to be generalized. Other drawbacks of this approximation are the
lack of order control and the number of partition along p cannot be too small even if the

Ap is much smaller than A/ 10. These disadvantages can be solved theoretically by

applying the Ty(2,4) FDTD to cylindrical coordinates FDTD.

4.2.3 Source Implementation

Traditionally, the way chosen to numerically incorporate an electromagnetic field
excitation source inside the cavity for the empty cavity simulation is based on the specific
cavity mode of resonance. Implementing a source is to select several lattice points as
source points and assigning magnitudes of 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 of exciting a
mode inside the cavity and then turning off the source gives the natural frequency
response. Examples of 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 Blackman-
Harris function for g(t), correct field distributions can be obtained. In this chapter, the
second-order FDTD with traditional source implementation is denoted by T-FDTD and

that FDTD with Blackman-Harris source is called BH-FDTD.

146

Perfectly conducting walls

 

 

 

 

 

 

 

 

 

 

 

7r / N
Excitation circle
/
. o“ ., ..-. ...}
‘ -...s..-.....47
Excitatioré
I Q
TEo" mode TE!“ mode
Er
I (+)
TM01] mode TMou mode

Figure 4.3 TE“ 1, TEnl, TM01], and TM012 modes excitation techniques
in a cylindrical cavity[35].

147

4.2.4 Numerical Results and Discussions

For TM012 mode, an empty cylindrical cavity with 0.0889 meter radius and
0.14409926 meter height are used. The number of partitions along p, (1), and z are 29, 36,
and 29, respectively. One time step is equal to 5x10“13 second and the number of total time

steps is 35,000. The total stored energies of T-FDTD and BH-FDTD are plotted in Figure

4.4. The result of BH-FDTD is closer to the correct one. The field distributions of TM012

mode of T—FDTD and BH-FDTD are plotted in Figure 4.5 to Figure 4.8. For T-FDTD,

only the field distribution of ED along 2 in Figure 4.6 and that of E2 along 2 in Figure 4.8

are correct comparing to theoretic results. On the contrary, the BH-FDTD gives the correct

field distributions of E field in all the figures.
Using an empty cylindrical cavity with radius equal to 0.0889 meter and height equal

to 0.0669089 meter, a TE 1“ mode can be excited. The number of partitions along p, q),

and z are 59, 36, and 29, respectively. One time step is equal to 1x10'13 second and the

number of total time steps is 175,000. The total stored energies of T-FDTD and BH-FDTD

of TE111 are plotted in Figure 4.9 and the field distributions of 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-FDTD are all incorrect. However, the field distributions of ED
along 4), Ep along 2, 15¢ along (1), and E,» along 2 are roughly similar to correct ones. Again,
the field distributions calculated by BH-FDTD are all correct comparing to theoretical
results. Hence, the traditional source implementation needs to be replaced by Blackman-
Harris 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

 

2.5 1 r r r 1

Total Stored Energy (Arbitrary Units)

 

 

 

 

 

 

 

0 0 5 1 1.5 2.5 3
Time (psec)

1-4 I I I T I l I 1—

1.2 ~ -
’8
3:
C
3 1 _ q
a
‘3
.2
e
< 0.8 ’ .I
V
E
8
LU 0.6 ~ *
E
o
(75 0.4 b A
E
O
.—

0.2 ~ ~

0 l 1 l l 1 l l
o 2 4 e a 10 12 14 16 13
Time (psec)

Figure 4.4 The plots of total stored energy of TM012 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 BH source
implementation.

149

 

    
  
   
 
 

 

 

 

 

 

 

 

0.6 r r r 1 I r r r
0.4 *- d
z=0.07205
0.2 -‘
71?
J:
C
2) oc .
Z‘
.2
'9
$-02 -‘
O.
UJ
“04 _(
-0.6 ..
_08 l l l l l l l 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
p (meter)
0.12 r I r I I l I I
0.1 *- ..
50.08 " ..
'c
D
Z‘
S 0.06 . -
g
'E
S.
UJQ
0.04 - z=0.07205 ‘
0.02 ..
0 I 4 r 1 l r 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
9 (meter)

Figure 4.5 The variation of Ep of TM012 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.

150

 

0.8 *-

   
  
 

.0 .o
# O)
I I

.9
N

E (Arbitrary Units)
,5

 

 

 

0 0.02 0.04 0 06 0.08 0.1 0.12 0.14

 

0.8

0.6

0.4

0.2

(Arbitrary Units)
52

Q02

E

 

 

 

l l l l l

o 0.02 0.04 0.06 0.08 0.1 0.12 0.14
z (meter)

 

-0.8 1

Figure 4.6 The variation of Ep of TM012 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 BH

source implementation.

151

 

    
   

=0.07205

2

E (Arbitrary Units)

l
N
I

-41-

 

 

-6 L 1 1

o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
p (meter)

 

 

0 I I I I I I I I

-o.3 _ z=0.07205

E (Arbitrary Units)
4 '3

 

 

l l 1 l l l

l l
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
p (meter)

Figure 4.7 The variation of E2 of TM012 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.

152

 

Ez (Arbitrary Unit)

 

 

 

 

0.02 0.04 0.06 0.08 0.1 0.12 0.14
2 (meter)

 

0.4

.0
to

(Arbitrary Units)

114-0.2

E

 

 

 

l l l l

0.02 0.04 0.06 0.08 0.1 0.12 0.14
2 (meter)

 

-O.8 ' ‘
0

Figure 4.8 The variation of E2 of TM012 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 BH
source implementation.

153

 

8

 

 

 

 

Total Stored Energy(Abritrary Units)
6‘. 8 3’. 8 8

_n

O
I
l

 

 

-
-

0 L l I l l

25 3 35 4

 

Time 2(psec)

 

0.09 I I I I I I I I

.o .o 9 .o .o
9 8 8 3 8
fly I I *I I
l l l l 1

Total Stored Energy (Abritrary Units)
0
'8

 

 

.0

8

T
1

F3

O

.5
I
l

 

 

 

o L L l L l 1 L
6

8 1O 12 14 16 18
Time (psec)

Figure 4.9 The plots of total stored energy of TEul 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 BH source
implementation.

154

 

0 ' p=0.04445 _
z=0.033454 ‘

I I
b N

I
0)

Ep (Arbitrary Units)

I
on

—10

-12

 

 

 

 

 

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

 

(Arbitrary Units)

E

 

 

 

I l

 

“065 l l l J 1 l
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

p (meter)

Figure 4.10 The variation of ED of TEnl 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.

155

 

 

 

 

15 l I I I T I
II
10 -
0 Q a 9
. o J L‘ ‘
A 5 ’- ~ ' d)
g I‘ 0
c
D
2:. o 0
g 0 - - p=0.04445 . -
\\ '1

2 z=0.033454
V
uJQ 0 0

-5 — c “ " e d

f \
b - a v . 0
-10 ~ -
U
_15 L 1 L L l l
0 50 100 150 200 250 300 350
o (degree)
08 I I I I T I

 

0.4

F3
N

(Arbitrary Units)

c:0.2

E

 

 

l l
50 1 00 1 50 200 250 300 350
¢ (degree)

 

 

Figure 4.11 The variation of ED of TEul mode along the 4) 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.

156

 

 

 

 

 

 

 

 

2 I I I I I I
0
p=0.04445
0< . “ ¢=160° .. .. ~
0
fa ‘2 '- 0 "
3:
c C a
D u
b a
e .4 i- -4
.2 " a r
2 ~ 9 .
v 9
G. c '4
Lu _6 _ .. a g a
g
'2 3 'J
..8 t- “ -1
_10 1 1 1 1 1 1
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07
2 (meter)
m I I l I I I
—o.1 — -
p=0.04445
-o.2 ~ ¢=160° _
A
.52
C
3 ~03 L -
2‘
t!
.1:
g-.. -
Q
[Li
—o.5 - -
-O.6 r- -4
-0.7 1 1 1 1 1 1
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07
2 (meter)

Figure 4.12 The variation of E9 of TEul 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 BH
source implementation.

157

 

 

 

 

 

4 I I I I I I I I
n
0 ‘ ~ 0
2 '- _.
n
U
U
A 0 " - O 0 "
a: 41-1 60 a
.E
:3 " z=0.033454 "
Z‘ “
e ‘2 " n n a 0 0 4
-G ” 5 a 0 “ 0 . v 9
n
“J -4 _ 0 '1 “ I, 9 0 .1
0
0 U c . ¢ 0 u
o
u
-6 1- ..
‘J
-8 1 1 1 1 1 1 1 1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
p (meter)
0 I I I I I I I I
-o.05 — ¢=1eo° ~
=0.033454
-0.1 _
73
g
C
D -o.1s ~
2‘
E
:5
s -0.2 ‘
0
DJ
—0.25 _

 

 

 

1 l l

l l I 1

-o.35 1
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
p (meter)

Figure 4.13 The variation of E.» of TElu 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

 

   
 

5 ~ p=0.04445

  
   

z=0.033454

E ¢ (Arbitrary Units)
0

-5

-10

 

 

 

 

-15 1
0 50 100 150 200 250 300 350
1) (degree)
0.5 I r r I I I
(14~ .
03- s

C’
A:

OJ

E (Arbitrary Units)
,5
- c:

 

  

1 l l l

1
50 100 150 200 250 300 350
«p (degree

Figure 4.14 The variation of E» of TEul mode along the ()1 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.

 

 

159

 

 

 

 

 

 

 

 

m I I T I I I
-1 - p=0.04445 1
¢=160°
A-2 "' .1
a:
.t
C
:3
2‘
g _3 .. -
e
:5
O
m _4 _ _
-5 - d
-6 1 1 1 1 1 1
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07
2 (meter)
T I I T I _I—
-o.05 - p=0.04445 ‘
¢=160°
A
01
3:
c -o.1 ~ .
D
2‘
.2
E
5-0.15 ~ ~
9
m
—o.2 - -
_025 L 1 1 l 1 1
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07
2 (meter)

Figure 4.15 The variation of E4, of TElu 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 BH

source implementation.

160

lossless cavity calculation whether if the cavity is empty or not.

For a loaded lossless cylindrical cavity, the TM012 mode is considered. The
dimensions of this empty cavity are 0.0762 meter in radius and 0.15458 meter in height.

The numbers of partition along p, (t), and z are 29, 36, 29, respectively, and one time step is

equal to 5x10’”. A small cylindrical material sample with 0.0105 meter in radius, 0.0107
meter in height, and 2.5 in relative permittivity is placed in the center of the cavity. The

field distributions of ED and B2 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: of TM012 mode along the z direction in the PEC cavity loaded with a small

cylindrical material 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 of partition along axes are increased and the dimensions
of the loaded small cylinder is less than tenth of those of the cavity.

Form these numerical results, we can conclude that the BH-FDTD with traditional
treatments at p = 0 gives correct field distributions of a lossless cylindrical cavity.

However, the drawback of BH-FDTD is its mode-dependent due to these treatments at
p = 0. By using the Ty(2,4) FDTD 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)
FDTD scheme is then expected to be the fundamental FDTD algorithm for cylindrical

coordinates.

161

 

0.2

0.181-

0.16 r

  
    
 

O

L.

h
I

.0

_e

N
I

z=0.0721 37 -

9.0
8

E (Arbitrary Units)

9
8

 

 

L L 1 l L 1 1

o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.09
p (meter)

 

(Arbitrary Units)

«0.05 ..l

E

 

 

 

 

’ o 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
2 (meter)

Figure 4.16 The variation of Ep of TM012 along the p and 2 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 BH
source implementation.

162

 

    
 

 

 

 

 

 

 

 

 

0 I l I I I I I
-o.1 - .
—o.2 ~ ~
A
m
:9:
C
D -o.3 z=0.072137 ‘
Z‘
.33
2-... -
11.1N <
.05 .
-o.6 ~
_07 1 1 1 1 1 1 1
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.011
p (meter)
0.8 T I l I I I I
A
(a
_O: 4
C
D
2‘
g .
E
:5
~02 ‘
LIJ
-O.8 1 1 1 1 1 1 1
o 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

2 (meter)

Figure 4.17 The variation of E2 of TM012 along the p and 2 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 BH

source implementation.

163

 

1.1 - -

1.05 - -

2

Ratio of E
O
8

.0
to
I
1

0.85 - -

 

 

J 1 1 l 1 1 1 1 1

0.75
0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0.005 0.09 0.095
2 (meter)

 

Figure 4.18 Plot of the ratio of E z/ E; of TM012 mode along the z direction in
the PEC cavity loaded with a small cylindrical material sample.

164

 

 

4.3 Ty(2,4) Cylindrical FDTD Scheme

Applying the implicit staggered fourth-order approximation,(2.57), to spatial
derivatives of fields in (4.9) to (4.14), the Ty(2,4) cylindrical FDTD scheme is obtained. A
general treatment at p = O for Ty(2,4) cylindrical FDTD method is discussed in section
4.3.4. The matrix equations along p, 0, 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;

59112, 69(pH¢) , 5pE and 5p(pE¢).

z,

The derivatives of H fields are considered first. Applying (2.57) to every lattice points

in Figure 4.19, the following matrix equation can be obtained

 

 

 

 

 

 

 

AX = M (4.23)
where
'1 22 1 0 0 0
0 1 22 1- O O
A = ' ° ' , (4.24)
O - l 22 1 O O
O O 1 22 1 0
_0 0 1 221_(Np-1)x(1vp+1)
T
x: 91 91 3_f .91 if , (4.25)
app=Oapp=l app=2 8pp=Np-lapp-Np
and

16S

59112.

0:0 801’; l p=2 . . . 691.121 p=Np—l 50H;
I I I l

”I ”I
Hzl p=l/2 Hzl p=3/2 z p=Np—3/2 z p=Np—1/2

 

 

 

 

p=1 =Np

50H:
1

5p(p11¢1p=0 59(DH¢1D=1 8‘4qu

p:2 o o o
l I . . . | l
pH¢| p=1/2 pH¢| p=3/2 pH¢| p=Np—3/2 pH¢l p=Np—1/2

69(pH4’1 pTNp-l 5"“) H11 p=Np

Figure 4.19 The locations of H field and its corresponding derivatives
in p direction.

166

 

flp=3/2—f|p=l/2

f|p=5/2—f|p=3/2
24

M = — ' .
Ap . (4 26)

f|p=Np-3/2_flp=Np-5/2

 

 

bf|p=Np-1/2-f|p=Np'3/22

where N p is the number of partition along the p direction and f is H z or pH ¢. In order to

use LU decomposition, a square matrix is desired. Hence, the derivatives of H fields at
p = O in Figure 4.19 needs to be approximated since (4.23) is a underdeterrnined linear
equation. To maintain the fourth-order approximation for the overall FDTD 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

 

 

 

3f —93f|p=%+229f|p=%—225f|p=%+lllflp=%-22f|p g
— = (4.27)
app=o 24Ap
Six-Order
if. _ b b —b b
app=0 f|p=-+ ZfID—g 3f|0=§+ 4f! ‘2 (428)
+b5fl M+b6f| l1+b7f|
_ =7 9:2
where

167

{—5.60208331
19.024479
-30.882813
(b1,b2,b3,b4,b5,b6,b7)T = 30.369792 . (4.29)
-18.026042
5.9640625
{-0.8473983)

 

 

The results by applying these two different order approximations will be discussed later.

Hence, the new matrix equation is AX = M where X is the same as (4.25),

 

 

 

 

 

"24 0 0 0 0 0
1 22 1 O O O
0 l 221- O O O
A = . . . , (4.30)
O O - - - l 22 1 0
O O 1 22 1
_0 ° ' ' 0 0 0 24_(Np+1)x(Np+1)
and
Apg—f
P p=0
f|p=3/2_flp=l/2
24 f|p=5/2-flp=3/2
M = — . (4.31)
A9
f|p=Np-3/2—f|p=Np-5/2
flpsz—l/Z-f|p=Np-3/2
Ap§—£
. p=Np _
where 3% is the approximation of 3%: . Another expression of (4.23) can be
p=0 p=0

 

168

 

obtained as:

 

 

 

 

 

r22 1 0 0
l 221- O O
A = . . . , (4.32)
0 - - - 122 1
_ 0 1 22.1~.-1)><(~.—1)
T
X: Q]: 3_f 3_f (4.33)
app=l pp=2 app=Np-1
and
r _ -
Apaf
24 f|p=5/2‘flp=3/2
M = — - . (4-34)
AD
f|p=Np—3/2—f|p=Np-5/2
A937
f|p=Np_1/2’f|p=Np—3/2_fia_pp_NJ
‘ p

 

 

—

The condition number of a matrix A measures how unstable the linear system AX = M is

under the perturbation of the data M. In the FDTD calculation, a small condition number
which is close to 1 is desirable. The later matrix equation is used in this chapter since the

condition number of matrix A in (4.32) is smaller than that in (4.30) for a same number of
partition.

For derivatives of E fields along p direction, the SpEz and 89(pE¢) at p = 1/ 2 and
p = N p — 1/ 2 has to be approximated regardless of the physical boundary condition.

The locations of the field and its derivative is plotted in Figure 4.20. The matrix equation

169

 

is the same as (2.73) to (2.76) and is rewritten below:

 

 

 

 

 

 

 

 

AX=M (4.35)
where
P26 -5 4 —1 0-
1 22 l O - - O
O 1 22 1 O - 0 1
A: . . . . . . . , (4.36)
0 0 l 22 l 0
0 O 0 l 22 1 p
_0 O —1 4 —5 26 prNp '3
T
21 if if 3_f
X: 0p lap BP 339 1 , (4.37)
9‘5 ”=6 “Np-5 “”97
and
flp.1-f|p=o
f|p=2_f|p=l
24 -
=— . 4.8
M Ap . (3)
fip=Np-l—f|p=Np—2
bflpsz—flpsz—ld

 

 

Note that the pE¢ in Figure 4.20 is equal to zero no matter what E9 is. The
p

approximations used for end points are

170

 

 

 

 

 

 

E.
” p=0 EZ 10:1 Ezlp= Ezl =Np-l EZlFFNP
I 1 ‘l
6 El I I
p z p=1/2 5DEZIp:3/2 SPEZ p=Np—3/2 SPEZ p=Np—l/2
0
.J
”Eflpa) 9591164 9591' p=2 pEtlp-er-l pE‘blTNp
l l l l

5p(pE¢)l p=1/2 5p(pE¢)| p=3/2

59“) E1) I p=Np—3/2 89(pE¢)| p=NP-1/2

Figure 4.20 The locations of E field and its corresponding derivatives in

p direction.

171

 

 

 

 

 

 

 

 

 

_iif +1591 _iif £521 f|p=l fl
248p 7 128p 5 248p 3 128p 1 Ap
”=2 9:5 “=2 ”=2
and
£31 521 +131 .1521
128p 1 248p 3 68p 5 248p 7
p”’9‘2 “”07 9:”9‘5 p= p-i

 

_ f|k=Np_f|k=Np-l
_ Ap

4.3.2 Derivatives of Fields in 0 direction

(4.40) 1

 

For derivatives of fields in 11) direction, there is no end points or boundary points;

hence, (2.57) can be used for all the lattice points along ((1 direction. The 8¢Hz and 8¢Hp

are considered first. Applying (2.57) to the FDTD lattices in Figure 4.21, the matrix

equation for this case is

 

 

AX = M
where
”22 1 0 1-
1 22 1 0
A = 0
0 l 22 1
_1 O l 22_N°)(N.
T
X _ 91 91' . . .91) 31'
a¢¢=oa¢¢zl a¢¢=N°_2 a¢¢=N°-1
and

172

(4.41)

(4.42)

(4.43)

f|¢=l/2—f|¢=Np-l/2

f|¢=3/2_f|¢=1/2
_ 24
_ A9

 

For 505; and 5¢E the matrix equation is

p9
AX=M

where A is the same as that in (4.42),

 

 

 

f|¢=~,—3/2”fl¢=N,—5/2

_fl¢=N¢-l/2_fl¢=N,—3/2_J

!

 

21 3_f 21
X= 61> .69 3 61> .
¢=§ 9:5 9=~,—§
and
i f|¢=1‘f|¢=o
f|0=2_f10=1
_24 -
“-11;
f|¢=N,—l—f|¢=N,-2
_ f|¢=0"f|¢=N,—1

 

 

173

 

(4.44)

(4.45)

(4.46)

(4.47)

Hz, Hp E

,,,,,,, ( 5,139, 5,1,5, )
_ ————— 0:0 0:19,,
~~~~~~~~ 84,112, 8,11p
~~~~ ( E. E.)

¢=N¢— 1

 

 

Figure 4.21 The FDTD lattice along the 0 direction.

174

4.3.3 Derivatives of Fields in z direction

The calculation of derivatives of fields in the z direction is similar to that in 2.3.4 and

2.3.5. The matrix equation for derivatives of magnetic fields, 52H¢ and SZHp , is

 

 

 

 

 

 

 

 

 

(4.49)

T
z=Nz—l:|

 

For PEC boundary, the a—f and 3_f
32 z = 0 az

 

M
where
_22 1 0-
l 22 0
A: 0
0 1
_ 0 22.1~.—1)x(~.—1)
(if (if 0f (if
X- a— 8‘ 8‘ a—
zz=l Zz=2 zz=Nz—2 Z
and
F ...
Azaf
fz=3/2—fz=1/2"2_43_Z
z=0
fz=5/2_f|z=3/2

flz=N,—3/2'flz=~,—5/2
42%?
Lflz: ~,—1/2—f|z=~,—3/2_2—4§E

 

 

z=N

Z

need to be approximated by (4.27) or (4.28), respectively for FEC boundary.

zp’

175

 

z=N

(4.48)

 

(4.50)

(4.51)

are all equal to zero, but these quantities

For derivatives of electric fields, BZE¢ and 5 E the matrix equation is (4.48) with

 

 

 

 

 

 

F26 —5 4 —1 . . 0'
1 22 l 0 - - O
0 1 22 1 O - O
A = . . . . . . . , (4.52)
O 0 l 22 1 O
O 0 O 1 22 1
_ 0 - O -1 4 —5 26JNZxNz
a_f a_f <21: 31
X - 32 1 82 3 02 3 dz 1 (453)
2:5 2:5 z=Nz-i z=Nz-i
and
f|z=l—f|z=0
f|z=2—flz=1
24 -
M = — . 4.54
A2 ( )

f|z=~,—1_f|f=N,—2
_ f|z=N,_f|z=Nz-l _

 

 

4.3.4 FDTD Calculations at p = 0
Due to the singularity at p = 0 , special treatments need to be applied to those finite

difference equations which contain divergent terms. At the first glance, the first two terms

on the right hand side of (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 = O . Thus, the FDTD equations from (4.9) to (4. 14)
at p = O are examined one by one in this section.

Equation (4.9) is still valid for p = 0 and is used to calculate D0 which is located at

176

half-integer space step along the p axis. For p = 0 , this equation becomes

 

 

 

 

I l 1 l 1
n-- n-- n—" ’1'“ "—-
2 1 2 2 2 2
i9 j,k 5, j,k i9 11k i 11k 51 11k it Jrk

Equation (4.10) is valid for p = 0 if 25sz at p = O is known. Here we cannot only
approximate 8,04, at p = O and ignore 52Hp at p = 0 since it is required in (4.51).
Moreover, HD at p = 0 is also required in (4.51); hence, HD 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.

BD at p = O is approximated by the following equation

1308p
+ _—
12 0p

3 —B
- pp=1 pp=0

 

 

 

 

 

13
ll
NIL»

p:

NI—‘

which is a fourth-order approximation. In order to obtain the derivatives of Bp with

respect to p, both aBp/ap at p = 1/ 2 and p = Np — 1/2 need to be approximated by

 

 

 

— — 2 B 11 — 2

3&9 = 93Bp|p=l+2293pp=2 25 pp=3+ lBPp=4 2 Bpp=5 (457)
0p 1 24Ap '

p—-

and

38" 93B 2293 2253

— = _ + _

8p 1 ( p|p=Np—l p|p=Np-2 p|p=Np—3

“NP—2 , (4.58)
1113 —228 / 24A
+ p|p=Np_4 p|p=Np_5) ( P)

177

2‘.

respectively. The matrix equation is AX=M where

 

 

 

 

 

 

 

 

 

 

 

 

22 1 0 O
1 22 1 - 0
A = 0 (4.59)
0 1 22 1
- 0 ' 0 1 2%(Np-21x1Np-2)
T
X _ 3'32 LB. 61% LB.
" 80 3 89 5 89 5 31) 3 ' (4'60)
P=§ 9:5 p_Nz-§ p=Nz-é
and
r _ -
B I _B I _Agf’fl»
24 Bp|p=3-Bp|p-2
M = _ . (4.61)
Ap B
p|p=Np-2- p|p=Np—3
B _ 3%
plp=~,-1 plp=~,_2 24 ap p=~,-1/2
a—Bp 63—, . . B.
where 0— and 0— are the approxrmatrons of ~3— at p = 1/ 2 and
p p=l/2 p p=Np—l/2 p
p = Np-1/2, respectively. The values of Bp at p = l to p = Np—l can be

calculated by (4.12) first, vector X is solved, and then Bp at p = O can be obtained. The
value of 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

 

tr

conduction current inside region S in Figure 4.1, we can obtain the finite difference

equation for DZ at p = 0 by Ampere’s law as

 

1 n n-l 1
"‘5 DZ" +DZI 4 "‘5
5D +08 = —H (4.62)
’ 30 .k 1 Z Z 2 r1 1’1 . 1
’4 +2 24**5

where r1 is the distance along p for the first cell. A more general approximation can be
obtained by similar approximation which is presented in the approximation of HD at
p = 0. Replacing H p by E z from (4.56) to (4.61), the fourth-order approximation of
£2 at p = 0 can be obtained.

Equation (4.13) is valid for p = 0 and it becomes

n n n
8’3‘1’1 . 1 = 69521 . l—SzEpl . 1' (4.63)
i,],k+§ 5,],k'1'5 i,j,k+§

Equation (4.14) is valid for p = 0, and it becomes

n l n
8,3 = —————8¢Ep _ ——z‘1p(pr:q,)1 1 (4.64)

P . - _
9k ’1... k 211+2)k

Nil—-
NI—
N."—

lj+lk 111+
2’ 2’ 2

where the location of p in 59(pE¢) depends on that of E¢.

4.4 Time Domain Finite Difference Equations of Constitutive Relations

The electric permittivity used in this chapter is modeled with the general Debye or
Lorentz equations.

The Debye equation is

179

- n t 8 — 8 °°
6(6)) = 8'((1)) — )6 ((0) = 8 .. + 2 fir, (4.65)
k = 1 9k
and the general Lorentz equation is
kag

 

M
8((0) = e'(w)—je"(w) = e'm+(e's—8'°°) 2 (4.66)
k=l

where wk is the k-th resonant frequency, uk is the k-th damping frequency, 8'00 = 8(00),
M
8'0 = 8(0) and 2 Gk =1.

k=l

A _\
The scalar constitutive relations modeled with Debye relation between E and D in

time domain can be obtained as:

 

M
Dam = e'amEaU) + 2 Pak(t) (4.67)
k = 1
where
aPakU) ' '

and 0t is one of p , q), or 2. By the same way, the scalar constitutive relations modeled

with Lorent equation in time domain is

M
Dam = E'amEaU) + (s'as—e'am) 2 Pak(t) (4.69)
k = 1

where

azpaku) apako)

2 2
T ”(x/(T +wakpak(’) = Gama/6150:“) (4'70)

180

 

 

and only the damping case with 031k — 401)?“ < 0 is considered.

The constitutive relations for Debye material in finite difference form are

M
Du " = E'amEaln’t 2 Pakl" (4.71)
k = 1
and

n n v v n
Pakl + Tackstpakl = (8 ask _ E aoo)Ea| (4'72)

and that for Lorentz material are

M
Da " -_- e'amEaln + (e'as — e'am) 2 PakI" (4.73)
k = 1

and

n

n 2 n 2
oakfitPakl +wakpakl = GakmakEoc . (4.74)

71
5,210

01k

If the second-order time stepping is used, then (4.72) and (4.74) can be rewritten as

-2 2At , .
" +r-(8654-864Wal" (4.75)

aek

n 2A1 n-l
- " Pak

 

P
ak 1:

+ Pakl
aek

and

2 —l —2 —3 —4
Pakln = (— 2y—ymak)Pak|" +(1+5y)Pak|" -4YPak|" +YPakln
(4.76)
n—l
+YGakw121kEa

 

where y = 2/(uakAt).

181

 

 

CHAPTER 5

CONCLUSIONS

In this dissertation, the finite-difference time-domain(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 of FDTD methods, the stability condition and
the numerical dispersion becomes major problems when it is applied to EM problems.
Generally, smaller Courant-Friedrichs-Lewy(CFL) values lead to the more stable FDTD
scheme. However, smaller CFL values will cause a larger phase error which is not
desirable in the FDTD calculation. Hence, the maximum CFL value for a stable FDTD
scheme which leads to a minimum phase error and also demands a smallest number of
grid cells in one wavelength is desired for FDTD calculation. This optimization can be
achieved by the analysis of the numerical dispersion equation for the FDTD scheme which
is demonstrated in section 2.2.2.

The traditional second-order FDTD 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 of grids and phase error will 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 FDTD schemes reduce

182

u.
u:—

 

the number of grid if the same phase error for the second-order is required. However, the
explicit higher-order FDTD scheme not only increases the stencil with more field points
but also complicates the treatment of lattice points near the physical boundary. The
implicit higher-order FDTD schemes uses the same number of field points as that in Yee’s
algorithm. The increased calculation time by the matrix multiplication in implicit higher-
order FDTD scheme is compensated by the decrease of the number of grids in one
wavelength. The Ty(2,4) method, which employed the second-order approximation in
time stepping and the implicit 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 of microwave, the
Ty(2,4) is highly applicable to higher frequency EM problems or microwave problems
with small size devices.

The quality factor, Q, and the resonant frequency are two important measurement
factors of cavities. By using the time-domain Poynting’s theorem derived in section 2.4.2,
Q value can be obtained from the FDTD calculation and (2.116). The resonant frequency
can be obtained by applying the fast Fourier transform (FFT) to field data verse time. For
lossless cavities case, the FFT approach may be practical since we can control the time
step to turn off the source. For lossy cavities case, the FPT approach may become
impractical since the number of time steps to reach the steady state can be very large. On
the other hand, both the quality factor and the resonant frequency of 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 of Q and the resonant frequency by Prony method is close to

those by time domain Poynting theorem and the FFT.

183

 

 

The singularities introduced by the FDTD method at the center of a cylindrical cavity
is a major problem in the application of the FDTD method. For cylindrical cavities loaded
with a symmetric material sample, the BOR FDTD method is the most efficient method to
perform the FDTD calculation as shown in chapter 3. The treatment of singularities in the
BOR FDTD method is also discussed in section 3.1.3. Since the BOR FDTD method is a
2.5D FDTD scheme, the calculation is much faster than that in a 3D case. Hence, the
number of 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 with the metal wall having a

conductivity larger than 104.

For general cylindrical cavities loaded with samples of arbitrary shape, the 3D
cylindrical FDTD method needs to be employed. The traditional second-order cylindrical
3D FDTD method suffers the problem of the mode-dependent source implementation and
the treatment of sigularities. By using the Blackman-Harris(BH) source, it solves the
source problem in the traditional second-order 3D FDTD method. In order to obtain a
general 3D FDTD method, a fourth-order treatment for sigularities is proposed in section
4.3.4 to replace the mode-dependent treatment in second-order FDTD method. Combining
with the implicit staggered fourth-Oder FDTD method, this proposed Ty(2,4) FDTD
method in cylindrical coordinates solves those problems in the traditional second-order
FDTD method. However, the implementation of this Ty(2.4) FDTD 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

 

range of EM problems. Some topics which are relevant for further study in the future are:

(1) Optimization of the numerical dispersion equation.

(2) Power analysis of other material models; e. g., Lorentz material.

(3) Extend the Ty(2,4) FDTD method for magnetic material.

(4) Investigating the implementation of proposed Ty(2,4) FDTD 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) FDTD

method, and perfect matched layer (PML) for the Ty(2,4) FDTD method.

185

 

 

 

APPENDIX A
DERIVATION OF FOURTH-ORDER FINITE

DIFFERENCE APPROXIMATION

In section 2.3.1, several fourth-order finite difference approximations to spatial
derivatives are presented. This appendix details the derivation of those finite difference
approximations based on the Taylor’s expansion. The explicit scheme have only one
unknown in the approximation; on the contrary, the implicit scheme have more than one
unknown in the approximation. The locations of 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, for the fourth-order approximation.

For explicit collocated scheme, the following Taylor’ expansions are used:

3
u- = u +Ax 14(1)+(___A2"!)2 .u(2)+(___Ax)

1+1 3! “1(3)+

-(-———A4x)4- uf4) + (A.l)

3 4
(Ax)? “(2) _ 93",; . “1(3) + 913+). - uf“) + (A-Z)

u!

ui_1= u —Ax u(1)+—

186

 

l

ZAJC2 2Ax3 ZAJC4
um = u,+(2Ax)-u,<1>+(—§—!l-u§2>+L§!_)_.u<3>+(—4!L.u,<4>+... (4.3)

2Ax 2 2Ax 3 ZAx 4
ui_2 = ui—(ZAx).ul(l)+g-T)~ul(2)—(—T)—-u,(3)+£—4!—)-u,(4)+.... (A.4)

Multiply (A.1) by or, (A2) by [3, (A3) by y , (A.4) by 8 and sum all these new equations

together; then setting the sum of the coefficients for “1', “(2) , up) to be zero and that of

u f 1) to be one. The following equation is obtained,

( a+B+y+6=0
a-B+2y—25=l
a+B+4y+45=0.
\ or—B+8y-86=0

(A5)

 

The solution for (A.5) is or = 8/12, [3 = —8/12, 7 = —1/12, and 8 = 1/12. Hence,

the explicit collocated approximation is

Bu ~ 8(“i+1—“i-1)"(“i+2‘ “(u-2)
(37),: 124x ' M

For the explicit staggered scheme, the following Taylor’s expansions are used:

Ax

“1+1/2 = “1*“2—‘1‘;

I l

2 3 4
1)+%/T2_)_ . “(2)+(_Ax3_/'22_ . “(3)+(_A_XZ£'2_) . “(4)4, (14.7)

Ax (Ax/2)2 (Ax/2)3 (Ax/2)4
_2.. . up) + T . "(2)—T . u(3) + T - 15(4) + (A.8)

ui—l/Z = “i 1 1

3Ax (3Ax/2)2 (3AJt/2)3 3Ax/2 4
“1+3/2 = “1+7'“il)+’T'“i2)+—3!——'uf3)+g-——z§—)—-u,(4)+.(A.9)

3Ax

__.u§

ui-3/2 = “i 2 1 r

2 3 4
1) + 9.13%31. .u12)_ $.3ng . “(3) + gL‘Z/fl. - 141(4) + (44.10)

Use the similar procedures in the derivation of the explicit collocated scheme, the

187

 

following set of equations is obtained

a+B+y+6=O

1 1 3 3
2 55:77—55“

1 1 9 9 (AH)
4a+-B+ZY+-5- 0
1 27 y_375_0

(8a—8B+8Y

 

The solution for (A.11) is or = 27/24, B = —27/24 , y = —-1/24, and 5 = 1/24

and the explicit staggered scheme is

27 u — u. —u.
(1“)F( (41/2 “1- 1/2) ( 1+3/2 1—3/2) (A12)
3x 24Ax

For the implicit collocated scheme, the following Taylor’s expansions are used:

a,“ = u +Ax u;1>+(——A") u;2>+%x!i3.u,<3>+£iff!—)4.ug4>+... (4.13)
u,_, = u,_m.u;1>+(i‘2£!—)2.u§2>_%l—3-u;3>+%x!—)3.u,<4>+... (4.14)
611>1=u<1>44x “(2)45—2—3’3") ug3>+£i‘3j—)3.u§4>+(—%"T)3.u§5>+... (4.15)
(4,19,: <I>_Ax-u;2>+(i‘zi‘!._)2.u;3>-(i‘3j—)3.u,<4>+(£“f!—)4-u;5>+.... (A.l6)

The corresponding equations for or, B, y, and 8 are

188

 

a+B=O
(or—B)Ax+y+5=1

1 l
_ ._ _ - A.l7
(01+ lB)Ax+(y 5) 0 ( )

l l l 1
K (ga— EB)AX 4' (EV ‘1" 2'6)-

and the solution is Ot = 3/(4Ax), B = —3/(4Ax), y = —l/4, and 5 = —1/4. Hence,

 

the implicit collocated scheme is

 

 

(123.1) ..(1131)
84141 ax 1—1 2 Bu ~“1+1‘“i—1
6 231371,: 2Ax ' (“8)

For the implicit staggered scheme, the following Taylor’s expansions are used:

I

Ax (Ax/2)2 Ax/2)3 (Ax/2)4
ui+ 1/2 = “1' 3- ' “(1) + T— - 141(2) + (—-3—!-—— ' “1(3) "1' T ' “1(4) + (A.l9)

2 3 4
92.x . “(1)...Ex12—L . “(2)_(_1§x_/_g)_ . “(3).6954fi . “(4)... (A20)

ui-l/2 = “1' 1 2! 1 3! 1
4
6,11), = u(')+Ax u(2)+(——Ax) 6,13>+(——A")2. 614>+(—A—"—)—.u13>+... (4.21)
3. 4! 3
_ (Ax)2 (A102 (A304
uflll — ufll—Ax-ufz)+—2!—~uf3)--§—!—-ul(4)+—Z-!—-u115)+.... (A22)

The corresponding equations for or, B, y, and 8 are

189

2 a+B=O
Gor-éfi)Ax+y+6=l

l 1 (A23)
(gt! + §B)Ax + (y—8)= O

1 1
([4—802— 4—8 )Ax+(§y+1—=05)

and the solution is a = 12/(11Ax), (3 = —12/(11Ax), y = —1/22, and 5 = —1/22.

 

Hence, the implicit staggered scheme is

 

(A.24)

 

(‘3’!) +035)
64141 ax1-1+1_1(E)u)i5“1+1/2“1—-1/2
24 12 8x Ax °

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

Bu Bu au Bu “1"“0
5(—) + (—) + (—) +61 61(— —) . (4.25)
ax 7/2 Y 8x 5/2 B 31‘ 3/2 ax 1/2 Ax

Using the following Tayor’s expansions:

 

Ill

11,)2 = 1104—11 111,1)+L_—Ax2/2)2 11211, +E‘gflfi-1133>+(A§a§2—)2.u54>+u. (A26)
“(1) = ugh-Ax u12)+(——A2x)2- ui3)+-(——A32:)3- uj4)+-(—é42:—)4-uj5)+... (A.27)

113/2 = “0+3A_x. 11-1-81)+(-§————A222/2)2- “82)+(-3———A;/2)3- 1136) +QA—ZT/y-u34)+...(A.28)
us/2 = uo+§%—x-u31)+(—5Ax—2!/—2)—2-u32)+(—5—A—232—!/2)3-u63)+(—5A:-—!/—2):-u341+...(A.29)

and

190

71'

 

. . .(A.30)

7Ax/2 2 (7Ax/2)3 7Ax/2 4
u7/2-u0+——-u “+5” (___ ) ”‘82)+—_31_—"“83)+L4—!‘)"'“84)+

the solution is 0: =13/12, B = —5/24, y =1/12, and 5 = —l/24.
The one-way fourth-order finite difference approximations in (2.81) and (2.82) can be

obtained by assuming

 

(Bu) 0m1/2 + Bus/2 + Y“5/2 + 6“7/2 + ““9/2
_ = (A.3l)
8x 0 Ax
with (A.26), (A.28), (A.29), (A30) and
2 3 4
119/2 = u0 + 9_Ax 1131) + (—————9A:/2) 1152) + (—————9A;/2) M6 + ———(9A:'/2) 4484) + ...(A.32)

The solution is 01 = —3l/8, B = 229/24, y = —75/8, 8 = 37/8,andn = -11/12.

191

I .
.' nfl

APPENDIX B

DERIVATION OF NUMERICAL DISPERSION

RELATION FOR IMPLICIT STAGGERED

SCHEME

In (2.70), the numerical dispersion relation of the implicit staggered FDTD scheme is
presented. The derivation of this equation provides fundamental derivation of those matrix

equations in Ty(2,4) FDTD scheme. In this appendix, the derivation of (2.70) is detailed.
Consider Maxwell’s equations in a normalized region of free space with u = 1,

8 =1,0' =1,andc = l andobtain

<L

a

ijV = (B.1)

Q.)

t

where V = H +jE and

vV—~avz avy ‘an avz Aavy avx 132
x ‘ (W‘s?)”(az‘$i”(x-wi ‘"

Let

V" V j(KxIAx+K).JAy+KzKAZ+0)nAI)
lI,J,K = 0 (3'3)

192

where xx, K and Kz are the wave numbers along x, y, and z axes, respectively in the

y 9
computational space.

Consider the x component of (B.1) at point (I, J, K) which is

 

 

 

 

 

8V 3V 8V
iii-711i = 715 03-4)
y z I, J, K I, J, K
The second-order finite difference approximation of the time derivative of V11 is
11 + 1/2 11 - 1/2
an = VII,J,K ' VII,J,K = _2 .sin(wAt/2)V (B 5)
31 ”K Ar 1 At Jrl1,J,K' '

 

The time derivatives of Vy and VZ can be obtained by replacing index x with y and 2. For
BVZ/ a y component, VZ is assumed to be located at the half integer position in space and
a V2/ 3 y is located at the integer position in space. The one-way fourth—order

approximation for numerical boundary used is assumed to be

av,
W

5V

z (8.6)

 

' 1+bVZ +cVz +dVz +eVZ
J="
2

 

 

 

 

 

J =

1:

MIL»
NIL]!
NI“

for j = O. The numerical approximation for j = N y can be obtained by similarity and is

 

 

 

 

listedbelow:
avz
— =c‘zV +bV +cV +dV
dy._~ Z-_Nl Z._N z._N 5 z _N_
1- 1 1- 1-2 ,_ ,_ r5 1- ,5 111.11
+eVZ
9
I=Ny_§

Apply the implicit staggered four-order finite difference approximation, (2.57), to every

lattice point along y axis, a set of equations is obtained as follow:

193

 

 

 

 

 

 

 

 

 

.(B.8)

: Vzlj=N_,—3/2_ Vzlj=N,-5/2

 

Ay

~ Vz|j=~,—1/2- Vzlj=N,—3/2

 

23:2 .1321: +1511: 2
24 By j=0 12 By i=1 24 By i=2
2315 ,1_13_V_z +1.32: 2
24 ay 1=1 12 By j=2 24 By F3
51% +21% .1112
24 By 1:11,-3 12 ay 1=~1-2 24 ay 12er
.2va +2312 .2912 -
24 By I‘M—2 12 By j=N,-1 24 By j=N,

 

 

 

Combining with the numeric boundary condition in

rewritten in a matrix form of

 

 

 

 

 

Ay

(B.6), the above equation can be

 

 

AX = B (3.9)
where
P 1 0 o 0-
_1_ 2_2_ i O
24 24 24
A : 0 (3.10)
0 i 22 i
24 24 24
_ 0 0 0 1_<~y+1>x<~y+1>
T
X_ avz avz avz avz avz (1111)
8y j=0 8y j=l ayj=1~/,—2 3y j=N,—1 3y j=N ,

194

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

c‘zVz +sz +cV.’ +de +eVZ
._1 ._§ ~._§ ._2 ._g
1’2 1’2 "2 "2 "2
- —V
‘j=3/2 Z|j=1/2
Ay
- —V
‘j=5/2 Z|j=3/2
M
B = .(B.12)
V —V
zlijy—3/2 zlijy-S/Z
M
v- —Vz|
Lj=N),-I/2 j=Ny—3/2
Ay
ZzV, +bV, +cV,l +de +eVz
&'=N-1 ~‘.—.N-.3. ”ex/-2 ‘=N_Z '=N-2
_ J .v 2 1 .v 2 J .v 2 J .v 2 J y 2‘
The matrix B is rewritten as follow
a b c d e zl1‘=1/2
_ V
1 1 0 0 0 zlj=31/2
O -l 10 O
B _ .1. — —l-MV (B 13)
‘ Ay ’ Ay ‘
0 0—1 1 O
O 0 0 —l l z|j=N,—3/2
e d c b Z1 VI
' ‘ _ Zj=N,-1/2J

 

 

 

where M is a (Ny+1)xNy matrix and l7 is a vector with N), dimension. Hence, the vector X

can be obtained as

-1.

11717.
Ay

X _ Aym M)V _

For a point with j = l, the derivative of VZ is listed below:

195

(3.14)

(B.15)

171,3 2B 1
+
= zlij— 2 9K

3VZ 1N
3y I, '=1,K =yA

 

 

yB

where 717,3 is the element of M— at 1th row and [3th column. According to (B 3)

can be rewritten as follow:

 

 

 

 

 

 

       

Z j=(ZB+ 1)/2
_. _2B+l
J[K’(I 2 )Ay] (
VZ = z ._ 'e . (B.16) T
111=2B2+1K l,j—I,K #1
Hence, (B.15) becomes “‘
N _1 _. _2B+I

8V: Z 1.1: 1K ’2 ,7, e J[Ky(l 2 )Ay] (B17)

__ _ IB . .

3y l,j=l,K Ay 5:0
Similarly, we obtain

N —1 _- .3114
BV YII J k: z _ J[Kz(q 2 )AZ]
73—1 _ ’ ’ q We . (B.18)
Z 11‘]! k = q AZ Y = 0
Hence, (B.4) leads to
Vxl Vzl ”y" ‘i‘rlJ-ZB + ‘1‘ l
1 J K . I J K _ 2
—-—’1—-231n (nAt/Z +-—’-’—- m e
At ( ) Ay B - 0 JB
- . (B.19)
_VYII,J,KN mm Y2+1) AZ] -0
AZ Y :0

Similarly, the y component and 2 component of (B.1) at point (I, J, K) are

196

 

y VXI N,—1
4&5 - 23in((1)At/2) + 45. 2 171W:
Y=0

 

 

 

 

A1 A2
2 1 (B20)
N _1 _. _ (1+
VZ|1,J,K x _ J[K‘(l 2 )Ax] _
— Ax - 2 mlae — 0
01:0
and
- V Nx-l —j 1c l-2a+l Ax
. y x 2
——-———XIJ’K~ZSin((1)At/2)+ LI’XLK- Z Filae [ ( ) ]
“=0 2 1 . (3.21)
V-‘|1,J,K Nil_ -J[ "(J- B2+ )A] 0
— - m e =
M 6:0 ’3
Equations (B. 19), (B20), and (B21) can be rewritten in a matrix form,
AV = O (8.22)
where A is 3x3 matrix as shown in (B27) and
VxlI,J,K
V = Vy|[,J,K . (3.23)
V
_ ZII,J,K_J

 

 

By taking the determinant of A to be zero and restoring the light speed into the time related

terms, the dispersion relation is obtained as

197

' ‘ l

 

 

 

 

 

 

 

 

 

_N -1_ —)[1<,(l-2a2+ IJAX] 2 -N,-l_ -J|:Ky(J 262+1)A] 2
[23in((1)At/2)]2 20 mla + B = 0771186
cAt Ax Ay
_ _ u ‘(B.24)
-N —1 ej[K(K-2Y2+1)AZ] 2

 

 

 

Assuming uniform partition along x, y, and z axes, A = Ax = Ay = Az, the

dispersion relation becomes

 

 

 

2 2
12[Sin(_ _):l2a1t ___ ”231111“, ‘j['<1("2a2+1)“i + Nilmjfie‘jixi’ ‘26; 1M
“/31? 01:0 B=0
(3.25)
Nz‘l : K_2Y+l A
+ Zmeej[K( 2 H
y=0

where a = fiery/A which is the CFL number and R = x/A which is the number of
cells in one wavelength. Let K = xxx + Kyy + KZE, then (B25) can be further simplified
as follows,

1 j_Ksin9cos¢(21-2a—l)-2

i—ilsmtn‘g—Rll= {-244 .

a 0 _
—1_stinBSin¢(ZJ-ZB—l)q

' R
+:=0[ 2 m me (326)

I:Nz._ 1 _jt'ccose(2K- 2y - 1)]2

+ 2 fizKYe R
y=0

 

N

 

 

198

 

(B.27)

 

 

 

 

 

 

 

199

 

 

 

 

 

 

21L 1. “Rina: E 2 L 1. 113:
M mafia \fixfl N U W M NCQQ \Tnfl N
532,523 muc puo
2 E Dve
1\_quflk 1 ~Q~+ JEN— Zu 1 _ 4.70.? 1 ~m~+ _
e534 maiekbv eminesan
E 9 Pa
241— 1. <s1~@+_ D 3.1. 1. e~1~Q+_ E
M flxum ‘7.“ N V L M $38 CTN N V g .
who Quo 553538
E R

 

 

 

where

= . (13.28)

 

Hence, the ratio of wave propagation speed in the computational space to that in the

physical space is equal to the ratio of TC to K.

1?“ Jay ~

1 h

l:
'1»

 

APPENDIX C

DERIVATION OF EQUATION (2.107)

In section 2.4.2, the time domain version of Poynting theorem is derived and g Jr(t) is a
key factor in calculating stored and dissipated power in time domain. In this appendix, the
derivation of time domain expression of g Jr(t) is detailed.

The function, 811(1) , is the inverse Fourier transform of E x((o)/ [l + (0)1“)2] and can

be rewritten as follow:

g(x) = S"{ E,(m) }

l + ((D‘tex)2

_ _1 -1 ___1___
_ 3 {Ex(w)}®3 [{1 +(w1ex)2}]
:13

Ex(t) ® {211 fl}

(C.l)

 

ex

 

1 1 fl

T
= e "E (t—T)d‘r
I x

where E x(t — 1:) is a causal signal; that is, E Jr(t — 1:) is equal to zero when t< T and the

. -r/n,. . .
system is a causal system; that 1s, 13 I I 1s equal to zero when 1: < 0. Discretmg (C1) by

using the approximation in Figure C.1, the discrete version of g x(t) is

201

 

 

 

mAt (m+1)At HA!

\ -(111 + 2)A1
1(1

e Ex(n—m—2)

 

 

 

mAt (m+1)At (m+2)At

FIGURE C.1 Approximation of integrant in (CI)

202

1 mAt _(m +1)A1

 

 

 

 

 

 

 

 

n— F- _ --
g(n) = 2: Age I”Ex(n—m)+e 1" Ex(n—m—l)
exm=0 - 4
1 "_1 A _ _mAt _(m+1)At _
= —t e t“'Ex(n—m)+te 1” Ex(n-m—1)
21am_ 1 2
_ ' " . (C.l)
_fl
1 At 1.1
+ ZTex-Z—[Ex(n) +e Ex(n—l)]
_A’ .91.
T11 At At T111
— e g(n— l)+4texEx(n)+4texe Ex(n—1)

A more accurate recursive approximation can be obtained as follows by using three

field points in Figure C.l:

 

 

 

"_2 J33 _(111+1)A1
_ l l T“ 4 Tu
g(n) — 21 At 3e Ex(n—m)+§e Ex(n—m—l)
exm=0
_(m+2)At
1 Ta
+§e Ex(n—m—2):|
11-2 _"W 422.112
= 2: 2: 1%[6 t”Et('I-m)+4e 1" Ex("-m-1)
“m =1 (C.2)

(m+2)At
T

+19 ” Ex(n—m—2):|

 

 

 

 

_E -23
+ 1 I-S—I[Ex(n) + 4e I"'Ex(n—l) + e T"‘Ex(n—2)]
21“ 3
_E -91 -291
I At 2A! 1 At 1
= (I - (X — (IE _2
e g(n l) +6TexEx(n)+ 3712116 Ex(n 1) +6Texe x(n )

where (C.2) needs to store one more history of E x than (C.l) does.

203

APPENDIX D

DERIVATION OF EQUATIONS (3.7) AND (3.8)

The term V x 75‘ in (3.2) is

V x 2 (Eucosmo + Evsinmo)

= 2 [V x (Evsinmqn + V x (Eucosmml
(D.l)

[(sinm¢V x Ev + Vsinm¢ x 2v) + (cosm¢V x Eu + Vcosm¢ x 2a)]
. .\ "IA .3 "I" .x A
[smm¢(V x ev - 6(1) x e“) + cosm¢(a¢ x ev + V x 4%)]

and those terms in the right hand side of (3.2) lead to

—[ul%f7l + [0*]1‘9’ = 410% 230(3ucosmtb + szinmtb)

+ [0*] 2 (zucosmq) + zvsinmtb)

m=0

A (D2)
hu A
= zocosm¢(—[u]% +[o*]hu)

+ 2 sinm¢(— [1&ng [0*]Zv) -

m=0

204

From (DI) and (D2), we obtain (3.8),

A

igffi x 2v, u + Vxéu, v = {we}; " + [0*]72u, v. (D.3)

 

For V x I? in (3.1), the following equation is obtained,

V x Z (flucosmo +71vsinm¢)

m=0

[V x (i‘zucosm¢) + V x (szinmm)
0

u
Ms

m

(D.4)

[(sinm¢V X 71v + Vsinmtb X in) + (cosmtpV x it“ + Vcosm¢ x7210]
0

M
Ma

m

[sinm¢(V x12, — gt?) x Z“) + cosm¢('-g$ x in, + V x ha]
0

M
Ma

m

and those terms in the right hand side of (3. 1) lead to

[8]%—f+[o]i+js = [e] EGG-:"cosmo-t-Eaa—jvsinmcp)

+ [o] 2 (Eucosmt’p + Epsinmo)
0° "' = 0 (13.5)
+ 2 (iucosm¢+jvsinm¢)

m=0

= i [cosm¢([e]g—f“ + [0]?“ +12) + Sinm¢([8]%v + [612v + 12)].

0

Comparing (DA) and (D5), we obtain (3.7),

a?“ r .
at +[O]eu’v+]u’v. (D.6)

 

£56) x 72v, u + szu, v = [e]

205

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

BIBLIOGRAPHY

John, D., Jr. Anderson, Computational Fluid Dynamics: The Basics With
Applications, McGraw-Hill, New York, NY, 1995.

K. S. Yee, “Numerical solution of inital boundary value problems involving
Maxwell’s equations in isotropic media,” IEEE Trans. Antennas and Propagation,
vol. 14, no. 3, pp. 302-307, 1966.

A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic
scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans.
Microwave Theory and Techniques, vol. 23, no. 8, pp. 623-630, 1975.

D. H. Choi and W. J. R. Hoefer, “The finite-difference time-domain and its
application to eigenvalue problems,” IEEE Trans. Microwave Theory and
Techniques, vol. 34, no. 12, pp. 1464-1469, 1986.

A. Navarro, M. J. Nunez, and E. Martin, “Finite difference time domain FFI‘
method applied to axially symmetrical electromagnetic resonant devices,” IEE
Proc., vol. 137, no. 3, pp. 193-196, June, 1990.

A. Navarro, M. J. Nunez, and E. Martin, “Study of TEO and TMO modes in

dielectric resonators by a finite difference time-domain method coupled with the
discrete Fourier transform,” IEEE Trans. Microwave Theory and Techniques, vol.
39, no.1, pp. 14-17, 1991.

F. Liu, I. Turner, and M. Bialkowski, “A finite-difference time-domain simulation
of power density distribution in a dielectric loaded microwave cavity,” Journal of
Microwave Power and Electromagnetic Energy, Vol. 29, pp138-148, 1994.

C. Su and J. Guan, “Finite-Difference Analysis of Dielectric-Loaded Cavities
Using the Simultaneous Iteration of the Power Method with the Chebyshev
Acceleration Technique,” IEEE Trans. Microwave Theory and Techniques, vol. 42,

no. 10. pp. 1998-2006, 1994.

F. Torres and B. Jecko, “Complete FDTD analysis of microwave heating processing
in frequency-dependent and temperature-dependent media,” IEEE Trans.
Microwave Theory Tech., Vol. MIT-45, pp108-117, 1997.

[10] J. Guan and C. Su, “Resonant Frequencies and Field Distributions for the Shielded

206

[ll]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

Uniaxially Anisotropic Dielectric Resonantor by the FD-SIC Method,” IEEE Trans.
Microwave Theory Tech, Vol. M'IT-45, no. 10, ppl767-1777, 1997.

J. Fang, “Time Domain Finite Difierence Computation for Maxwell ’s Equations,”
PhD. thesis, University of California at Berkeley, Berkeley, CA, 1989.

T. Deveze, L. Beaulie, and W. Tabbara, “An absorbing boundary condition for the
fourth order FDTD Scheme,” IEEE Antennas and Propagation Society
International Symposium, pp. 342-345, July, 1992.

S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” J. of
Computat. Phys., vol. 103, pp. 16-42, 1992.

J. L. Young, D. Gaitonde and J. S. Shang, “Towards the contruction of a fourth
order difference scheme for transient EM wave simulation: Staggered grid
approach,” IEEE Trans. Antennas and Propagation, vol. 45, no. 11, pp. 1573-1580,
1997.

M. H. Carpmter, D. Gottlieb, and S. Abarbanel, “The Stability of Numerical
Boundary Treatments for Compact High-Order Finite-Difference Schemes,” J. of
Computat. Phys., vol. 108, pp. 272-295, 1993.

A. Taflove, Computational Electrodynamics: The Finite-Dzfirerence Time-Domain
Method, Artech House, Boston, MA, 1995.

A. Taflove, Advances in Computational Electrodynamics: The Finite-DWerence
Time-Domain Method, Artech House, Boston, MA, 1998.

F. B. Hildebrand, Introduction to Numerical Analysis, McGraw-Hill, NewYork,
NY, 1974.

C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley & Sons, New
York, NY, 1989.

F. J. Harris, “On the use of windows for harmonic analysis with discrete Fourier
transform,” Proc. IEEE, vol. 66, pp.51-83, Jan. 1978.

J. A. Pereda, L. A. Vrelva, and A. Prieto, “Computation of resonant frequencies and
quality factors of open dielectric resonantors by a combination of the finite-
difference time-domain (FDTD) and Prony’s methods,” IEEE Microwave and
Guided Wave Letters, vol. 2, 1992, pp. 431-433.

R. Kumaresan, and D. W. Tufts, “Estimating the parameters of exponentially
damped sinusoids and pole-zero modeling in noise,” IEEE Trans. Acoustics,
Speech, and Signal Processing, vol. 30, 1982, pp. 1380-1419.

K. M. Chen, EE836 Class Notes, Michigan State University.

207

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

J. D. Jackson, Classical Electrodynamics, John-Wiley, New York, NY, 1975.

D. B. Davidson, and R. W. Ziolkowski, “Body of revolution finite-difference time-
domain modeling of space time focusing by a three dimensional lens,” J. Optical
Society of America, June 1993.

K. M. Chen, Notes for EE 836, Michigan State University.

D. B. Davidson, “Body-of-revolution finite-difference time-domain modeling of
space-time focusing by a three-dimensional lens,” J. Optical Society of America,
vol. 11, no.4, pp. 1471-1490, 1994.

Y. Chen, R. Mittra, and P. Harms, “Finite-Difference T ime-Domain Algorithm for
Solving Maxwell’s Equations in Rotationally Symmetric Geometies,” IEEE Trans.
Microwave Theory and Techniques, vol. 44, 1996, pp. 832-839.

J. G. Maloney and G. S. Smith, “The use of surface impedance concepts in the
finite-difference time-domain method,” IEEE Trans. Antennas Propagat., vol. 40,
no. 1, pp. 38--48, 1992.

S. Kellali, B. Jecko, and A. Reineix, “Implementation of a surface impedance
formalism at oblique incident in FDTD method,” IEEE Trans. Electromagnetic
Compatibility, vol. 35, no. 3, pp. 347-356, 1993.

K. S. Oh and J. E. Schutt-Aine, “An efficient implementation of surface impedance
boundary conditions for the finite-difference time-domain method,” IEEE Trans.
Antennas and Propagation, vol. 43, no. 7, pp. 660-666, 1995.

D. M. Sullivan, “ Z transforms theory and FDTD method,” IEEE Trans. Antennas
and Propagation, vol. 44, no.1, pp. 28-34, 1996.

J. H. Beggs, “A FDTD Surface Impedance Boundary Condition Using Z-
Transforms,” Applied Computational Electromagnetics Society Journal, vol. 13, no.
1, pp.l4-24, 1998.

J. G. Proakis and D. G. Manolakis, Digital Signal Processing, Macmillan, New
York, New York, 1992.

W. Y. Tan, “Modeling The Electromagnetic Field and The Plasma Excitation in A
Moderate Pressure Microwave Cavity Plasma Source,” Ph.D. Dissertation,
Michigan State University, 1994.

D. M. Sullivan, “A Frequency-Dependent FDTD Method for Biological
Applications,” IEEE Trans. Antennas and Propagation, vol. 40, no. 3, pp. 532-539,
1992.

Om P. Gandhi, Ben-Qing Gao, and Jin-Yuan Chen, “A Frequency-Dependent
Finite-Difference Time-Domain Formulation for General Dispersive Media,” IEEE

208

 

Trans. Microwave Theory and Techniques, vol. 41, no. 4, pp. 658-665,]993.

[38] J. L. Young, A. Kittichartphayak, Y. M. Kwok, and D. Sullivan, “On the Dispersion

Errors Related to (FD)2TD Type Schemes”, IEEE Trans. Microwave Theory and
Techniques, vol. 43, no. 8, pp. 1902-1910, 1995.

209

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

willwillitwwuymuun