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

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

1/--страниц