# Time-domain analysis of microwave devices by the finite-difference time-domain 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 free, while others may be from any type o f computer printer. The quality of this reproduction is dependent upon the quality o f 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. Each original is also photographed in one exposure and is included in reduced form at the back o f the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6” x 9” black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. UMI A Bell & Howell Information Company 300 North Zeeb Rood, Ann Arbor MI 48106-1346 USA 313/761-4700 800/321-0600 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R e p r o d u c e d with p e r m i s s io n of t h e co p y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n . R e p r o d u c e d with p e r m i s s io n of t h e co p y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n . ABSTRACT Title of Dissertation: TIME DOMAIN ANALYSIS OF MICROWAVE DEVICES BY THE FINITE DIFFERENCE TIME DOMAIN METHOD Gary Thomas Roan, Doctor o f Philosophy, 1997 Dissertation directed by: Professor Kawthar A. Zaki Department o f Electrical Engineering In this dissertation, the Finite Difference Time Domain (FDTD) numerical method is applied to the analysis o f stripline filters and discontinuities, resonant cavities, and microwave antennas. In the FDTD method, the six components of the electric and magnetic fields are calculated on a three-dimensional mesh as functions of time. The S-parameters, resonant frequencies, quality factors, and antenna patterns are calculated from the time histories o f the computed fields. In many cases, a single time history can provide an accurate response over a very wide frequency bandwidth. The sensitivities of the derived parameters to the size and uniformity of the mesh, to the boundary conditions used to truncate the mesh, and to the particular time dependence o f the input pulse used are investigated. Calculated parameters are compared both to measured values and to previously published results. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. New equivalent circuit models are developed for the stripline discontinuities using the scattering matrix derived from the FDTD calculations. It is shown that significant discrepancies exist between the FDTD results and previously published equivalent circuits. The FDTD results are validated by a comparison to measurements. The new circuit models that are derived are used to efficiently and accurately synthesize a low pass stripline filter. The FDTD method is used to analyze the parallel plate resonator technique for the determination of the surface resistance o f superconductors at microwave frequencies. The resonant frequency, quality factor, and the dielectric and conductor losses are calculated for a microwave enclosure containing a superconductive parallel plate resonator. Results of the FDTD analysis are used to determine experimental conditions that can lead to errors in the determination of the surface resistance. The FDTD method is used to examine the radiation characteristics o f biconical antennas. The effectiveness of transverse corrugation networks to reduce antenna sidelobes are also analyzed via FDTD and compared to measurements. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n . TIME DOMAIN ANALYSIS OF MICROWAVE DEVICES BY THE FINITE DIFFERENCE TIME DOMAIN METHOD by Gary Thomas Roan Dissertation submitted to the Faculty o f the Graduate School o f the University of Maryland at College Park in partial fulfillment of the requirements for the degree of Doctor o f Philosophy 1997 Advisory Committee: Professor Kawthar A. Zaki, Chair Professor Julius Goldhar Professor Victor Granatstein Professor Leonard Taylor Professor Mohamed Aggour Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 9808661 Copyright 1997 by Roan, Gary Thomas All rights reserved. UMI Microform 9808661 Copyright 1997, by UMI Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 300 North Zeeb Road Ann Arbor, MI 48103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ©Copyright by Gary Thomas Roan 1997 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n . Dedication: To my wife, Maureen ii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ACKNOWLEDGEMENTS I would like to express my appreciation to my advisor, Professor Kawthar A. Zaki, for help, guidance and encouragement towards the completion o f this dissertation. Her sound advice was invaluable in solving the many problems. I also thank Dr. Goldhar, Dr. Granatstein, Dr. Taylor and Dr. Aggour for agreeing to serve as dissertation committee members. I wish to thank my employer, the Naval Research Laboratory, and my immediate supervisor, Mr. Eugene Starbuck, for their continued support of this dissertation. I also thank Mr. Armondo Elia, for his helpful advice and comments, and Mr. Dan Hicks, Mr. Eddie Sines and Mr. Anthony Allegrezza for some excellent technical assistance. I thank my parents for impressing on me the need to work hard, to not give up, and to be thankful for what God has given me. The completion of this dissertation has been a long time in coming. Many good and bad events occurred during this period of my life. Fortunately, I was not alone. I had the love, the support, and the devotion of my beautiful wife. She contributed greatly by providing copious quantities o f encouragement, compassion, and patience. Together, we look forward from this milestone to the rest o f our lives together. iii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TABLE OF CONTENTS LIST OF FIGURES________________________________________________________________vi LIST OF TABLES_________________________________________________________________ x CHAPTER 1 INTRODUCTION______________________________________________________ 1 CHAPTER 2 DESCRIPTION OF THE NUMERICAL PROCEDURE______________________ 10 M a x w ell ’s Equ a tio n s ......................................................................................................................................... 10 F in ite D ifference For m o f M a x w ell ’s Equations ....................................................................................11 A dvantages o f FDTD fo r so lv in g M axw ell ’s Eq u a t io n s ....................................................................17 S o u r c e s .................................................................................................................................................................... 18 B o u n d a r y conditions .......................................................................................................................................... 18 Perfect Electrical Conductors (PEC) and Perfect Magnetic Conductors (PMC)............................. 18 Open (absorbing) Boundary Condition.......................................................................................... 20 St a b il it y lim ita tio n s ..........................................................................................................................................26 N u m er ic a l dispersion ..........................................................................................................................................27 CHAPTER 3 APPLICATION OF THE FDTD METHOD TO ANALYSIS OF RESONANT CAVITIES______________________________________________________________________ 30 C alcu la tio n o f Resonant frequencies using th e Fo u r ie r T ransform _____________ 30 R e so n a n t F requencies fo r a Rectangular PEC C a v it y .......................................................................32 D eterm ination o f Q uality F a c t o r ............................................................................................................... 33 A n a ly sis o f Cav ity w ith S uperconductive th in films f o r determination o f su r fa c e resistiv ity and com pa riso n w ith m e a su r e m e n t .......................................................................................37 A na ly sis o f m icrostrip du a l patch reso n a to r ........................................................................................53 C alcu la tio n o f th e r eso n a n t frequency o f dielectric resonators and co m pa riso n w ith MEASUREMENTS....................................................................................................................................................... 60 CHAPTER 4 APPLICATION OF THE FDTD METHOD TO ANALYSIS OF STRIPLINE CIRCUITS______________________________________________________________________ 65 P ropa g a tio n constant an d im pedance for sym m etric s t r ip l in e ...................................................... 65 D eterm ination o f S- param eters fo r stripline d is c o n t in u it ie s .........................................................77 Calculation o f S-parameters and equivalent circuitfor stripline openend discontinuity................. 79 Calculation ofS-parameters and equivalent circuitfor stripline step junction discontinuity...........85 Calculation o f S-parameters and equivalent circuitfor stripline tee-junction discontinuity.......... 100 D esig n o f low pass filters in c lu d in g stripline d isc o n t in u it ie s ......................................................118 C o m pa riso n of predicted fil t e r response w ith m e a s u r e m e n t s ....................................................... 129 C o m pl e t e analysis o f low pass filters using FD T D ............................................................................. 135 CHAPTER 5 APPLICATION OF THE FDTD METHOD TO ANALYSIS OF BICONICAL ANTENNAS--------------------------------------------------------------------------------------------------------140 C a lcu la tio n o f n ea r fields fo r an ax i - sym m etric b ico n ica l a n ten n a ........................................ 140 N ea r - h e l d to fa r field tra n sfo rm a tio n for determ ination o f antenna pa ttern s .................146 C o m pa r iso n of calculated a n d m easured antenna pa ttern s fo r a biconical a n t e n n a 152 C a lcu la tio n o f radiation pa tter n fo r an axi- sy m m etric bico n ica l antenna in c lu d in g a rf choice fo r side - lobe red u ctio n ......................................................................................................................156 C o m pa r iso n o f calculated a n d m easured antenna pa ttern s fo r a biconical a n ten n a in c lu d ed a rf c h o k e .......................................................................................................................................... 157 iv Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. CHAPTER 6 CURVILINEAR, CONFORMAL AND NONORTHOGONAL GRIDS__________163 C a l c u l a t io n o f r e s o n a n t f r e q u e n c ie s f o r a c y l in d r ic a l w a v e g u id e u s in g c y l in d r ic a l a n d CONFORMAL F D T D GRIDS__________________________________________________________ 163 C a l c u l a t io n o f r e s o n a n t f r e q u e n c ie s f o r a c y l in d r ic a l w a v e g u id e F D T D GRIDS__________________________________________________________ 172 u s in g n o n o r t h o g o n a l CHAPTER 7 CONCLUSIONS AND DISCUSSION OF FUTURE RESEARCH_____________ 186 APPENDIX A: A BRIEF DESCRIPTION OF THE FDTD COMPUTER PROGRAM________ 189 REFERENCES__________________________________________________________________ 192 v Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. LIST OF FIGURES F ig u r e 2 - 1. O r ientation and position o f el e c t r ic a n d m agnetic field c o m po n e n t s w ithin Y ee CELL-----------------------------------------------------------------------------------------------------------------------II Fig u r e 2-2. N u m erical phase v elo c ity v er su s g r id spacing . ....................................................................29 F ig u r e 3-1. T e st fixture fo r th e m e a su r e m e n t o f su rfa ce resista n ce ................................................ 39 Fig u r e 3-2. P l o t o f V # ? vs. spa c er t h ic k n e ss s . .......................................................................................... 44 F ig u r e 3-3. R ec ipr o c a l o f v s . s f o r d o m in a n t m ode w ith a L a A l O j s u b st r a t e ................ 50 F ig u r e 3-4. R e c ip r o c a l o f Q em u ku r e v s - s f o r d o m in a n t m ode w i t h o u t a L a A lO j s u b s t r a t e .........50 F ig u r e 3-5. R e c ip r o c a l o f Qenclosure v s . s f o r T M ^ m o d e w rra a L a A lO , s u b s t r a t e _____________ 5 1 F ig u r e 3-6. Rec ipr o c a l o f v s . 5 f o r TM oj m o d e w ithout a L a A l O j su b st r a t e ................... 5 1 F ig u r e 3-7. La y o u t o f m icrostrip patches in a PEC cavity o f d im ensions a =0.945 in ., B =(l.276+2D ) IN., AND YP=Zp=0 3 11 IN........................................................................................................ 54 F ig u r e 3-8. R eso n a n t frequency v ersus c e l l s ize f o r a m icrostrip pa tch r e so n a t o r from t h e FD TD RESULTS.................................................................................................................................................... 58 F ig u r e 3-9. R eso n a n t frequency versus p a t c h separation fo r a 5 0 x 7 4 x 1 5 1 c e l l FD TD g rid . .59 F ig u r e 3-10. Reso n a n t frequency versus p a t c h separation fo r a 101x151x302 c e l l FDTD GRID........................................................................................................................................................................ 59 F ig u r e 3 - 11. C a lcu la ted resonant fr eq u en c y fo r dual patch reso n a to r fr o m refer en c e 26 (R epr in ted w ith perm ission o f au th o r ) ...................................................................................................60 F ig u r e 3-12. D ia g ra m o f a c y l i n d r i c a l d i e l e c t r i c r e s o n a t o r in a r e c t a n g u l a r w a v e g u id e ...6 I F ig u re 4 - 1. La y o u t o f C artesian g rid fo r sy m m e t r ic stripline ...............................................................69 F ig u r e 4-2. FDTD voltage vs . tim e betw een c e n t e r co nductor o f striplin e a n d g r o u n d plane AT K=60, 120 a n d 180........................................................................................................................................71 F ig u r e 4-3. Im a g in a r y part o f pro pagation c o n st a n t v s . frequency fo r 0.036 in . w id e str iplin e ............................................................................................................................................................... 72 F ig u r e 4-4. FDTD cu rren t versus t im e a t k = 120...........................................................................................73 F ig u r e 4-5. FDTD im pedance vs . freq u en cy fo r a 0.036 in . w ide striplin e w it h 0.060 in . ground PLANE SPACING AND £*=3.0........ 75 F ig u re 4-6. C om pa riso n o f striplin e im ped a n c e v s . w idth from t h e FDTD r e su l t s a n d from pu blish ed fo rm u la for sev era l th ic k n esses o f METAL..................................................................... 76 F ig u re 4-7. st r ipl in e D iscoNnNurriEs and t h e ir eq u ivalent circu its .....................................................78 F ig u r e 4 -8 . in c id en t and reflected vo lta g es fo r a stripline o pen disco n tin u ity (W =0.020 in ., B=0.055 IN. £*=3.0)............................................................................................................................................. 80 F ig u r e 4-9. Ph a se o f S „ versus freq u en cy f o r a striplin e o pen disco n tin u ity f o r striplin e WIDTHS OF 0.010 IN, 0.020 IN. AND 0.040 IN. (B=0.055 AND £*=3.9)....................................................... 82 F ig u r e 4-10. C o m pa riso n of th e P ha se o f S „ c a lc u la ted w rra th e n ew fo r m u l a t o t h e FDTD RESULTS................................................................................................................................................................. 85 F ig u r e 4 - 1 1. I n c id e n t a n d r e f l e c t e d w a v e f o r a s t r i p l i n e s te p ju n c tio n (W , = 0.040 in. W 2=0.020 IN., B=0.055 IN., £*=9.8)..................................................................................................................................... 87 F ig u r e 4-12. T r a n s m i tt e d w a v e f o r a s t r i p l i n e s te p ju n c tio n (W ,=0.040 in. W 2=0.020 in., b=0.055 IN., £*=9.8)............................................................................................................................................................ 87 F ig u r e 4-13. M ag nitude o f s „ fo r striplin e st e p jun ctio n (W ,=0.040 in . W ,=0.020 in ., b=0.055 IN., £*=3.9)............................................................................................................................................................ 90 F ig u r e 4 -1 4 . p h a s e o f S „ f o r s t r ip liin e s te p j u n c t i o n (W ,=0.040 in. w 2=0.020 in ., b= 0.055 in., e*=3.9)....................................................................................................................................................................90 F ig u r e 4 - 15. M a g n itu d e o f s 2, S tr ip lin e s te p j u n c t i o n (W ,=0.040 in. w 2=0.020 in., b=0.055 in., Er= 3.9)....................................................................................................................................................................91 F ig u r e 4-16. P h a s e o f S^ f o r s t r i p l i n e s te p j u n c t i o n (W ,=0.040 in. W 2=0.020 in., b=0.055 in., £*=3.9)....................................................................................................................................................................91 F ig u re 4-17. M a g n itu d e of S „ fo r striplin e st e p jun ctio n (W ,=0.040 in . W ,=0.020 in ., b =0.055 in ., e*=9.8)................................................................................................................... ' .......................................92 vi Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FIGURE 4 -1 8 . PHASE OF S „ FOR STRIPLINE STEP JUNCTION (W ,=0.040 in . W ,=0.020 in ., B=0.055 in ., 92 £*=9.8)---------------------------------------------------------------------------------------F ig u r e 4-1 9 . M agnitude o f S j , fo r striplin e step junction . (W ,=0.040 in . W ,=0.020 in ., b =0.055 in ., £*=9.8)----------------------------------------------------------------------93 F ig u r e 4 -2 0 . Phase o f S2I f o r st r ipl in e step junction (W t=0.040 in . W2= 0.020 in ., b =0.055 in ., £*=9.8)----------------------------------------------------------------------------------------------------------------------------- 93 F ig u r e 4 -2 1 . M agnitude o f S „ fo r str iplin e step junction (W ,=0.040 in . w 2=0.020 in ., b=0.055 in ., £*=24)________________________________________________ 94 F ig u r e 4 -2 2 . Phase o f S m fo r st r ipl in e step junction (W ,=0.040 in . W2=0.020 in ., b= 0.055 in ., e*=24)______________________________________________________ 94 F ig u r e 4 -2 3 . M agnitude o f S,, f o r st r ipl in e step junction (W ,= 0.040 in . W2=0.020 in ., b=0.055 IN., £*=24)_______________________________________________________________________________ 95 FIGURE 4 -2 4 . Pfiase o f S2, fo r s t r ipl in e step junction (W ,=0.040 IN. W2= 0.020 in ., B=0.055 in ., e*=24)------------------------------------------------------------------------------------------------------------------------------ 95 F ig u r e 4 -2 5 . M agnitude o f s „ fo r str iplin e step junction fr o m c u r v e f it a n d FD TD (W ,=0.040 IN. W2=0.020 IN.,AND B=0.055 IN.)..................................................................................................................98 F ig u r e 4-26. P hase o f S m fo r st r ipl in e step junction from c u r v e f it a n d FDTD (W ,=0.040 in . W2=0.020 in ., a n d B=0.055 IN.).......................................................................................................................98 F ig u r e 4 -2 7 . M agnitude o f S21 fo r str iplin e step junction fro m c u r v e f it a n d FD TD (W ,=0.040 IN. W2=0.020 in ., AND B=0.055 in .).................................................................................................................99 F ig u r e 4 -2 8 . P hase o f S2, fo r s t r ipl in e step junction from c u r v e f it a n d FD TD (W ,=0.040 in . W2=0.020 IN., b=0.055 IN.)................................................................................................................................99 F ig u r e 4-29. In cid en t a n d r eflec ted w ave on line I o f a st r ipl in e T ee - ju n c t io n (W t=W 2=0.040 in ., W j=0.020 IN, B=0.055 IN., £*=3.9)..........................................................................................................102 F ig u r e 4 -3 0 . T ransm itted w a v e fr o m line 1 to line 2 o f a st r ipl in e T ee - ju n c t io n (W ,=W ,=0.040 IN., W j=0.020 IN, B=0.055 IN., e*=3.9)........................................................................................................102 F ig u r e 4 - 3 1. T ransm itted w a v e fr o m line i to line 3 o f a st r ipl in e t e e - ju n c t io n (W ,=W ,=0.040 IN., W j=0.020 IN, B=0.055 IN., £*=3.9)............................................................................................... ."........ 103 F ig u r e 4 -3 2 . Incid en t and re fl e c t e d w ave on line 3 fo r a st r ipl in e T ee - ju n c t io n (W ,= w 2=o .040 IN., W j=0.020 IN, B=0.055 IN., £*=3.9)......................................................................................................... 103 F ig u r e 4 -3 3 . M agnitude o f S „ f o r sym m etric stripline T ee ju n c t io n (W ,=W ,= 0.040 in ., W3= 0.020 IN., B=0.055 IN.).............................................................................................................................. 106 F ig u r e 4 -3 4 . P hase o f S „ fo r s y m m e t r ic stripline T ee jun ctio n (W ,= W 2= 0.040 in ., w 3=0.020 in ., B=0.055 IN.).........................................................................................................................................................106 F ig u r e 4 -3 5 . M agnitude o f S2I fo r sy m m etric stripline t e e ju n c t io n (W ,= w ,= 0 .0 4 0 in ., W j= 0.020 IN., B=0.055 IN.).............................................................................................................................. 107 FIGURE 4 -3 6 . PHASE OF Sj, FOR SYMMETRIC STRIPLINE TEE JUNCTION (W ,= W 2= 0.040 IN., W 3=0.020 IN., B=0.055 IN.).........................................................................................................................................................107 F ig u r e 4 -3 7 . M a gnitude o f S31 f o r sy m m etr ic stripline T ee ju n c t io n (W ,= w 2=0.040 in ., W3= 0.020 IN., B=0.055 IN.).............................................................................................................................. 108 F ig u r e 4 -3 8 . P hase o f S3l fo r sy m m e t r ic stripline t e e junction (W ,= W 2= 0.040 in ., w 3=0.020 in ., B=0.055 IN.).........................................................................................................................................................108 F ig u r e 4-3 9 . M agnitude o f S „ fo r sy m m etr ic stripline T ee ju n c t io n (W ,= w 2=0.040 in ., W3=0.020 in ., B=0.055 IN.)............................................................................................*.................................109 F ig u r e 4 -4 0 . P hase o f S33 fo r sy m m e t r ic stripline T ee jun ctio n (W ,= W 2=0.040 in ., W 3=0.020 in ., B=0.055 IN.).........................................................................................................................................................109 F ig u r e 4 - 4 1. N ew equivalent a R c u r r fo r sym m etric T ee- ju n c tio n ......................................................111 F ig u r e 4-4 2 . C om parison o f c u r v e fit fo r m agnitude o f s „ to FD TD resu lts (W ,=W 2=0.008 in ., B=0.060 in ., £*=3.0)...........................................................................................................................................114 F ig u r e 4-4 3 . C om parison o f c u r v e f it f o r phase of S „ t o FDTD r e su l t s (W ,= w 2=0.008 in ., B=0.060 IN., e*=3.0)........................................................................................................................................... 114 F ig u r e 4 -4 4 . C om parison o f c u r v e f it fo r m agnitude o f S2I t o FD TD resu lts (W ,= w 2=0.008 in ., B=0.060 IN., e*=3.0)........................................................................................................................................... 115 vii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. F ig u r e 4-45. C om parison o f curve fit fo r ph a se o f S*, to FDTD results (W \=W ,=0.008 in ., B=0.060 IN., e*=3.0)____________________________________________________________________ 115 F ig u r e 4-46. C om parison o f curve fit fo r m a g n itu d e o f S,, to FDTD results (W ,= w 2=0.008 in ., B=0.060 IN., £*=3.0)____________________________________________________________________ 116 F ig u r e 4-47. C o m p a ris o n o f c u r v e f i t f o r p h a s e o f S3, t o FDTD r e s u l t s (W ,=W 2=0.008 in., B=0.060 IN., £*=3.0)____________________________________________________________________ 116 F ig u r e 4-48. C o m p a ris o n o f c u r v e f i t f o r m a g n itu d e o f S 33t o FDTD r e s u l t s (W ,= w 2=0.008 in., B=0.060 IN., e*=3.0)____________________________________________________________________ 117 F ig u r e 4-49. C om parison o f curve fit fo r ph a se o f Sjj t o FDTD results (W ,=W ,=0.008 in ., B=0.060 IN., e*=3.0)____________________________________________________________________ 117 F ig u r e 4-50. (C ontinued ) (c ) st o p band r espo n se . __________________________________________ 120 F ig u r e 4-51. La y o u t o f seven elem ent str iplin e l o w pa ss filter ____________________________ 120 F ig u re 4-52. C alculated ( a ) passband S ,, r espo n se a n d ( b ) stop band S2, r espo n se fo r a seven ELEMENT STRIPLINE FILTER WITH A 4 GHZ CUTOFF FREQUENCY............................................................. 125 F ig u re 4-53. C alculated ( a ) passband s „ r espo n se a n d ( b ) stop band S2, respo n se fo r a seven ELEMENT STRIPLINE FILTER WITH A 8 GHZ CUTOFF FREQUENCY______________________________ 126 F ig u r e 4-54. c a l c u l a t e d (a ) p a ssb a n d s „ r e s p o n s e a n d (b) s to p b a n d S2, r e s p o n s e f o r a se v e n ELEMENT STRIPLINE FILTER WITH A 12 GHZ CUTOFF FREQUENCY_________________ 127 F ig u r e 4-55. C alculated ( a ) passband s „ r espo n se a n d ( b ) stop band r espo n se f o r a n ELEVEN ELEMENT STRIPLINE FILTER WITH A 8 GHZ CUTOFF FREQUENCY............................................. 128 Fig u r e 4-56. M easured and calculated ( a ) S „ response an d ( b ) Sj , response fo r a sev en ELEMENT STRIPLINE FILTER WITH A 4 GHZ CUTOFF FREQUENCY.............................................................131 F ig u r e 4-57. M ea su red and calculated ( a ) s „ respo n se and ( b ) S,, response f o r a sev en ELEMENT STRIPLINE FILTER WITH A 8 GHZ CUTOFF FREQUENCY.............................................................132 F ig u r e 4-58. M ea su red and calculated ( a ) S „ respo n se and ( b ) Sj , response f o r a se v e n ELEMENT STRIPLINE FILTER WITH A 12 GHZ CUTOFF FREQUENCY...........................................................133 Fig u r e 4 -5 9 . M ea su r ed and calculated (a ) S ,, respo n se and ( b) S2I response f o r a n eleven ELEMENT STRIPLINE FILTER WITH A 8 GHZ CUTOFF FREQUENCY.............................................................134 F ig u re 4-60. Illustration of a n o n - uniform FDTD g r id ...........................................................................137 F ig u re 4 -6 1 . co m pa r iso n o f response ca lc u la ted using FDTD t o m easured r espo n se fo r a SEVEN ELEMENT FILTER WITH A 12 GHZ CUTOFF FREQUENCY: (A) S „ RESPONSE AND (B) S2I RESPONSE............................................................................................................................................................ 138 F ig u r e 4-6 2 . c o m p a r is o n o f re s p o n s e c a l c u l a t e d u s in g FDTD t o m e a s u re d r e s p o n s e f o r a ELEVEN ELEMENT FILTER WITH A 8 GHZ CUTOFF FREQUENCY: (A) S „ RESPONSE AND (B) S2I RESPONSE............................................................................................................................................................ 139 Fig u r e 5 - 1. ( a ) D iagram o f biconical a n ten n a a n d ( b ) FDTD representation o f c o n ic a l ANTENNA ON A TWO DIMENSIONAL GRID.......................................................................................................141 Fig u r e 5-2. H* n ea r - field around a bico n ica l a n ten n a a t tim e step n =356......................................146 Fig u r e 5-3. V ir tu a l surface S for determ ination o f far - fields ........................................................... 150 F ig u re 5-4. M ea su r ed and calculated fa r field patterns for a biconical a n ten n a a t 8 GHZ.................................................................................................................................................................. 154 F ig u re 5-5. M ea su r ed and calculated fa r field patterns for a biconical a n ten n a a t 10 GHZ................................................................................................................................................................154 F ig u r e 5-6. M ea su r ed and calculated fa r f ie l d patterns for a biconical a n ten n a a t 12 G H z................................................................................................................................................................155 F ig u r e 5-7. M ea su r ed and calculated fa r fie l d patterns fo r a biconical a n t e n n a a t 14 GHZ................................................................................................................................................................155 F ig u r e 5-8. M e a su r e d and calculated fa r fie l d patterns fo r a biconical a n ten n a a t 16G H Z ................................................................................................................................................................ 156 F ig u r e 5-9. C h o k e netw ork for reducing sid elo b es o f a biconical antenna ..................................157 F ig u re 5-10. C a lcu la ted and m easured a n ten n a pa tter n fo r a b iconical a n ten n a w ith a rf CHOKE NETWORK AT 8 GHZ............................................................................................................................ 158 F ig u re 5-11. C alcu la ted and measured an ten n a pa ttern for a biconical an ten n a w ith a rf c h o k e n etw o r k at 9 G H z ............................................................................................................................ 158 • • • VUI Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. F igure 5-12. C a lcu la ted a n d measured an ten n a pa t t e r n f o r a biconical a n ten n a w ith a r f CHOKE NETWORK AT 10 G H Z .__________________________________ _________________________159 F igure 5-13. C a lc u la ted a n d measured an ten n a pa t t e r n fo r a biconical a n ten n a with a rf CHOKE NETWORK AT 11 GHZ_________________________________________________ ___________159 F igure 5-14. C a lc u la ted a n d measured a n ten n a pa t t e r n f o r a biconical a n ten n a w rra a r f CHOKE NETWORK AT 12 G H Z .___________________________________________________________ 160 F igure 5-15. C a lc u la ted a n d measured an ten n a pa t t e r n f o r a biconical a n ten n a w ith a r f CHOKE NETWORK AT 13 GHZ_________________________________________________ ___________ 160 F igure 5-16. C a lc u la ted a n d measured a n ten n a pa t t e r n f o r a biconical a n ten n a w ith a r f CHOKE NETWORK AT 14 G H Z .______________________________________________ _____________ 161 F igure 5-17. C a lc u la ted a n d measured a n ten n a pa t t e r n fo r a biconical an ten n a w ith a r f CHOKE NETWORK AT 15 G H z____________________________________________________________ 161 F igure 5-18. C alc u la ted a nd m easured an ten n a pa tter n fo r a biconical an ten n a w ith a r f CHOKE NETWORK AT 16 GHZ......................................................................................................................... 162 F igure 5 - 19. C a lc u la ted a n d m easured an ten n a pa t t e r n f o r a biconical a n ten n a w ith a r f CHOKE NETWORK AT 17 GHZ......................................................................................................................... 162 F ig u re 6- 1. S t a i r c a s e d a p p ro x im a tio n f o r a c y l i n d r i c a l r e s o n a t o r in a c o a r s e C a r t e s i a n FDTD GRID-------------------------------------------------------------------------------------------------------------------- 165 F igure 6-2. C o n to u r p a t h integration fo r updating t h e H z- field n e x t to a cu r v ed PEC boundary ..........................................................................................................................................................168 F igure 6-3. P lo t o f t h e Ex- fie l d versus t im e d u ring t h e ( a ) fir st 2500 tim e steps a nd ( b ) BETWEEN TIME STEPS 21500 AND 24000 WHERE THE INSTABILITY IS APPARENT................................ 171 F igure 6-4. N on -o r th o g o n a l g rid used perpen d icu la r t o z - a x is fo r 1.02 in. d ia m eter CYLINDRICAL RESONATOR.............................................................................................................................. 178 F igure 6-5. N o n o r th o g o n a l g r id for a cylindrical d ie l e c t r ic d isk inside o f a cylin d rica l WAVEGUIDE........................................................................................................................................................................... 181 F ig u re 6- 6 . N o n o r t h o g o n a l g r i d u se d f o r c y l i n d r i c a l d i e l e c t r i c in sid e o f a r e c t a n g u l a r WAVEGUIDE........................................................................................................................................................................... 183 F igure 6-7. R evised n o n o r th o g o n a l grid u sed fo r c y l in d r ic a l dielectric inside o f a RECTANGULAR WAVEGUIDE............................................................................................................................ 185 F igure A-7-1. F lo w d ia g r a m fo r basic FDTD alg o rith m ........................................................................ 190 ix Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. LIST OF TABLES T a b l e 3-1 reso n a n t freq u en cies f o r a 0.7 in . x 0.8 in . x 0.9 in . Re c t a n g u l a r C a v it y .................... 33 T a b l e 3-2 q u a l it y factors fo r a 0.7 in . x 0.8 in . x o .9 in . c a v it y ............................................................. 37 T a b l e 3-3. FDTD results fo r Q ,R e so n a n t F requency and S urfa ce r e s is t a n c e o f t h e Do m in a n t M o d e in a P ara llel Pl a t e R e s o n a t o r w ith a La A l O , su b str a te . .................. 45 T a b l e 3-4. FDTD results fo r Q , R e so n a n t F requency and S urfa ce R e sist a n c e o f t h e D o m in a n t m o d e in a P a r a l l e l P l a t e R e s o n a t o r w i t h o u t a L a A l 0 3s u b s t r a t e ....................46 T a b l e 3-5. C om parison o f c a l c u l a t e d a n d m easured results f o r a 0 .5 in . x 0.25 in . resonator ______________________________________________________________________________ 52 T a b le 3-6. C alculated r eso n a n t fr eq u en c y fo r a n u m ber o f d iffe r e n t FD TD g r id s ................... 58 T a b l e 3-7. Reso n a n t F req u en cies f o r a 1.00" x 1.00” x 0 .9 2 " R ec t a n g u l a r C a v it y w ith a C ylindrical D ielectric ( er =38) placed 0.275" above th e b o tto m o f t h e w a v eg u id e w rra a C a rtesian G rid o f 2 0 x 2 0 x 1 8 c e l l s ........................................................................................................... 63 T a b le 3-8. Resonant F req u en c ies f o r a 1.00" x 1.00" x 0.92" Re c t a n g u l a r C a v it y w ith a w rra a C ylindrical d ielec tr ic ( e * =38) pla ced 0.275" a bove t h e b o tto m o f t h e w av eg u id e w rra a C a rtesian G rid o f 4 0 x40x36 c e l l s ....................................................................................................... 64 T a b l e 3- 9 . r e s o n a n t F r e q u e n c ie s f o r a 1.00" x l . 0 0 " x 0.92" R e c t a n g u l a r C a v i t y w ith a w ith a C y l i n d r i c a l D i e l e c t r i c (e* =38) p la c e d 0.275" a b o v e t h e b o t t o m o f t h e w a v e g u id e w ith a C artesian G rid o f 60x60x54 c e l l s ....................................................................................................... 64 T a ble 4-1. V alues o f co effic ien ts f o r (4.13).................................................................................................. 84 T a b le 4-2. Values o f c o efficien ts in (4-22).......................................................................................................97 T able 4-3. Se t o f C oefficients f o r Eq u a tio n (4 2 9 ).................................................................................... 113 T a b le 4-4. N orm alized r eflec tio n z e r o s a n d transm ission zeros o f a se v e n e l e m e n t filter WITH otp=0.1773 DB............................................................................................................................................121 T able 4-5. N orm alized r eflec tio n z e r o s and transm ission zero s o f a elev en e l e m e n t filter WITH AP=0.1773 DB........................................................................................................................................... 122 T a b le 4-6. C haracteristic im ped a n c es a n d lin e lengths a nd w idths o f se v e n se c t io n filter FOR 4, 8 AND 12 GHZ CUTOFF FREQUENCIES................................................................................................122 T a b l e 4-7. C h a r a c t e r i s t i c im p e d a n c e s a n d l i n e l e n g t h s a n d w id th s o f e l e v e n s e c t i o n f i l t e r FOR A 8 G H z CUTOFF FREQUENCY.................................................................................................................. 122 T a b le 6-1. Resonant frequencies f o r a 1.02 in . diameter x 0.612 in . l o n g c y l in d r ic a l ca v ity using a staircased C a r tesia n FDTD g r id w ith 20x20x12 c ells , a x = ay = a z =0.05 I in ., and 32,768 TIME STEPS............................................................................................................................................. 163 T able 6-2. reso n a n t frequencies fo r a l .02 in . D iam eter x 0.612 in . l o n g C y l in d r ic a l C a v ity usin g a cylindrical FDTD G rid w ith 10x30x12 (R,<fr,z) c ells , Ar =A z =0.051 in ., a <J>=ji/15, and n =65,536 T im e St e ps .............................................................................................................................. 164 T a b l e 6-3. r e s o n a n t f r e q u e n c ie s f o r a 1.02 in. d ia m e te r x 0.612 in. l o n g c y l i n d r i c a l c a v i t y u s in g t h e c o n t o u r p a t h m e th o d w r r a 20x 20x 12 (x ,y ,z) c e l l s , a x = a y = a z = 0 .0 5 1 in., a n d N =16,384 t im e Ste ps ....................................................................................................................................... 169 T a b le 6-4. reso n a n t frequencies f o r a n em pty I .02" dia x 0.6" C y lin d r ic a l C a v it y using a N onorthogonal FDTD G r id in t o e X Y - plane ( non 20. f ) w ith 20 x 2 0 x 12 c e l l s ...................... 180 T a b l e 6-5. Resonant F req u en cies fo r a l .02" d ia x 0.600" C y lin d r ic a l C a v it y w it h a 0.68" dia . x 0.300" C ylindrical d ie l e c t r ic ( e* =35.74) using a n o n o r t h o g o n a l FDTD G rid in th e XYp l a n e w rra 20x 20 x 12 c e l l s ......................................................................................................................... 181 T a b l e 6- 6 . R e s o n a n t F r e q u e n c ie s f o r a 1.00" x 1.00" x 1.00" R e c t a n g u l a r C a v i t y w r r a a C y l i n d r i c a l D i e l e c t r i c (e* = 38) p l a c e d 0.275" a b o v e t h e Y =o p l a n e u s in g a c y l i n d r i c a l FDTD GRID INSIDE THE DIELECTRIC AND A NON-ORTHOGONAL FDTD GRID BETWEEN THE DIELECTRIC AND THE WAVEGUIDE WALL.......................................................................................................185 X Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 1 Introduction The objective of this dissertation is the application o f time domain techniques to the analysis o f microwave devices. Until relatively recently, time domain analysis has not been applied extensively to microwave problems. Generally, frequency domain methods such as the method o f moments1,2 (MOM), mode matching method (MM)3,4, and finite element method (FEM)5,6 were the dominant numerical analysis tools for high frequency simulations. Frequency domain numerical methods typically reduce the number o f unknowns for a given problem through pre-processing, thereby requiring less computer memory and less computation time. For example, a MOM solution for a three dimensional antenna problem might solve for the currents over a two dimensional surface only. The solution for the fields throughout space is then provided by a known Green’s function. FEM and Finite Difference Frequency Domain (FDFD)7 methods determine a solution implicitly by simultaneously solving a system of equations for the unknown fields or potentials. The usual rule of thumb is that a minimum o f 10 nodes per wavelength in each dimension are required. In the case of FDFD, three dimensional problems can lead to dense matrices with millions o f unknowns. The inversion o f such large densely populated matrices may not be practical or even possible. In FE methods, the system o f equations forms a sparse matrix, which is much more amenable to numerical solution then the densely populated matrices required by the FDFD. In this dissertation, the Finite Difference Time Domain (FDTD)8,9,10,11,12 numerical method is applied to the analysis o f stripline filters and discontinuities, I Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. resonant cavities, and microwave antennas. In the FDTD method, the six components o f the electric and magnetic fields on a three-dimensional grid are calculated as functions o f time. The FDTD method calculates the unknown fields explicitly, i.e., no simultaneous solution o f equations is required. The fields are solved everywhere within the domain o f the problem via a finite difference approximation o f Maxwell’s equations. At each time step, the curl equations are solved at every space point using only the fields from the previous time step in the immediate vicinity o f the grid point (i.e., nearest neighbor). Since no matrix inversion is necessary, problems involving huge numbers of unknowns can now be solved, albeit not always quickly. Kunz13 applied FDTD methods to the analysis of an entire aircraft. Several factors work to the advantage of the FDTD method: (1) The simplicity o f the update equations permits rapid calculation o f the fields on modem vectorized or super-scalar processors. Direct solution of large three dimensional problems are now practical. (2) The field update equations can be solved in parallel. Hardware is now readily available in the form of multi-processor computers that can perform multiple numerical computations in parallel. (3) The frequency response o f a device can be calculated over a large bandwidth using the Fourier transform of a time domain solution. This partially offsets the disadvantage of having to solve Maxwell’s equations directly for a large number of unknowns since a single time domain simulation provides a solution at many frequencies. (4) The study of transient phenomenon is another possible benefit of working in the time domain solution. (5) FDTD also provides a full wave solution directly from Maxwell’s equations. FDTD solutions are generally 2 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. accurate, are easily checked for convergence, and gross errors in programming or application are usually obvious. (6) Complicated structures are easily handled since a different material property can be assigned to every grid point in the problem domain without any serious effect on the speed o f solution. For example, the FDTD method has been used to calculate the absorption o f microwave energy by different organs o f the human body.14 This dissertation applies the FDTD method to obtain S-parameters, resonant frequencies, quality factors, and radiation patterns from the time histories o f the computed fields. In many cases a single time history provides a accurate response over a very wide frequency bandwidth. The sensitivities of the derived parameters to the size and uniformity o f the mesh, to the boundary conditions used to truncate the mesh, and to the particular time dependence of the input pulse used have been investigated. Calculated parameters are compared both to measured values and to previously published results. Chapter 2 provides a discussion o f the FDTD method within the context o f the applications addressed by this dissertation. The discussion includes a description o f the field update equations for Maxwell’s equations, sources, boundary conditions, stability limitations, and numerical dispersion. The finite difference update equations are developed directly from Maxwell’s equations in differential form. A key feature of the FDTD method is the fact that the finite difference equations approximate the partial derivatives in Maxwell’s equations to 2nd order in the grid spacing15. Thus, numerical solutions in free space converge to the true solution as the 2nd power of the grid 3 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. spacing, and decreasing the grid spacing has a rapid effect on the solution convergence rate. Most boundary conditions in the FDTD method are treated without difficulty. Perfect electric conductor (PEC) and perfect magnetic conductor boundary conditions are easily applied by setting the appropriate fields to zero. The continuity o f the tangential electric fields and normal magnetic fields across the grid is ensured by the basic algorithm. Development of open or absorbing boundary conditions (ABCs) continues to be important in many applications o f the FDTD method. ABCs are required in scattering, radiation and waveguide problems whenever it is assumed that the wave propagates towards infinity in one or more directions. The original work of Mur16 led to the first practical implementation of open boundary conditions for scattering problems. His 1st and 2nd order boundary conditions are still frequently applied because of their simplicity and effectiveness. For the calculation o f waveguide S-parameters, it is particularly important that the ABC provide a very small reflection coefficient (< 0.5%) or else large errors appear in the Fourier-transformed frequency domain results.17Improved ABCs such as the dispersive boundary condition1819 and the superabsorption boundary condition20 have been developed for printed waveguide analysis. A review o f ABCs for microstrip problems was presented by Betz.21 For the stripline analysis presented in Chapter 4, the superabsorption boundary provides nearly perfect absorption at the stripline terminations, and was found to be stable and easy to implement. 4 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. In the application o f the FDTD method to three dimensional scattering and antenna problems, none o f the ABCs discussed above provide reflectionless absorption over a wide range o f angle o f incidence. Recently, the Perfectly Matched Layer (PML)22,23 appeared that improved the absorption o f outgoing waves by several orders o f magnitude over previous ABCs. The PML ABC greatly enhances the usefulness o f the FDTD method for many scattering and antenna problems. However, use o f the PML comes at some cost, and it can more than double the required computer memory and computation time.23 Since other adequate ABCs were available for the applications in Chapters 3 ,4 and 5 that required little extra memory and processing time, the PML was not applied in this dissertation. In Chapter 3, the FDTD method is applied to the analysis of a variety o f microwave resonators. In the first case, the resonant frequencies and quality factors o f rectangular waveguide cavities are calculated for comparison to the known analytical results. It is shown that accurate results can be obtained using the FDTD method with only modest computer resources. In the next case, the FDTD method is applied to the analysis o f the parallel plate resonator technique o f Taber24 for the determination o f the surface resistance o f superconductors at microwave frequencies. An FDTD analysis o f the resonant frequency, quality factor, and the dielectric and conductor losses has led to the discovery of conditions that lead to significant experimental errors in the determination of the surface resistance.25 In the third example, the resonant frequencies o f a dual microstrip patch resonator are calculated and compared to both measurements and results obtained using MOM.26 Finally, the FDTD method is applied to the 5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. calculation of resonant frequencies o f cylindrical resonators inside o f rectangular waveguide. This problem has only recently been solved accurately using the mode matching method.27 In Chapter 4, equivalent circuit parameters are extracted from the FDTD analysis o f stripline discontinuities. Equivalent circuit models for discontinuities are often required in the synthesis o f filters, matching networks and couplers.28,29 Analysis and synthesis of microwave devices using equivalent circuit models are generally much faster and more efficient than full wave analysis methods. Here, the FDTD method is used to calculate the scattering parameters for stripline openend, step and Tee-junction discontinuities. The FDTD results are compared to equivalent circuits developed by Altschuler30 and Oliner.31 It is shown that significant discrepancies exist between the FDTD results and previously published equivalent circuits. New equivalent circuit models are developed for the stripline discontinuities using the impedance matrix derived from the FDTD calculations.32 The new circuit models for the stripline discontinuities are then used in the synthesis of low pass stripline filters.33,34 The new equivalent circuit model is then validated by a comparison to the measurements. Numerical analysis o f complete microwave circuits is often not performed because o f the difficulty in mathematically representing a complicated geometry. Previous published results on low pass filters considered relatively simple two pole structures.35,55 Here, a complete analysis of three and five pole low pass filters is carried out using the FDTD method. This analysis was conducted with and without a uniform Cartesian grid. A uniform grid with a cell size o f0.002 in. was required to obtain a good 6 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. match to the low pass filter stripline circuit In such as case, good agreement with measurements was obtained over a 26 GHz wide bandwidth from a single FDTD simulation that required approximately 8 hours o f central processor unit (CPU) time. A second grid was formed using a non-uniform cell spacing,36 in which case the cell size was varied according the width o f the metal striplines. When a non-uniform cell size is applied to a FDTD grid, the truncation error in the finite difference equations is 1st order in the cell spacing. However, Monk has shown that the numerical results will still be 2nd order convergent.37 Using the non-uniform grid, an accurate prediction o f the low pass filter was still obtained and the CPU run time has been reduced by factors as small as 0 . 1. In Chapter 5, the FDTD method is used to examine the radiation characteristics of biconical antennas. Biconical antennas are very popular in applications requiring a wide frequency bandwidth (> one octave) in an omni-directional antenna. A review o f available literature revealed little information regarding the elevation radiation patterns as a function o f frequency for biconical antennas.38’39,40,41’42 Maloney43 provided calculations o f the electric field for a cylindrical monopole antenna as a function o f time in the far field, but did not consider antenna patterns. In other publications, the FDTD method has been used to calculate the radiation patterns for a monopole antenna mounted on an automobile,44 for a Vivaldi Flared horn,45 and for waveguide horns and two dimensional parabolic reflectors.46 Here, the elevation patterns for a biconical antenna are calculated as a function of frequency using the FDTD method and compared 7 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. to measurements. Excellent agreement between the calculated and measured antenna patterns was obtained over a 7-17 GHz bandwidth. Another advantage o f the FDTD method is the ability to analyze complicated structures. The standard biconical antenna provides an omnidirectional beam about its vertical axis. In the elevation plane, the biconical antenna will provide a positive gain in radiation intensity near the center o f its beam relative to an isotropic radiator. The antenna pattern will display a null at 0 and 180 degrees that will be -30 to -40 dB relative to the peak gain of the main beam. Some situations require operation o f both receivers and transmitters at the same frequency from a single location. Such systems are called repeaters, and they find applications in navigational aids in harbors and in electronic warfare applications. To avoid crosstalk between antennas, the repeater’s receive and transmit antennas are often placed above and below one another in the respective nulls o f the opposite antenna. Often, insufficient isolation between antennas occurs, and isolation networks are added to the space between the antennas to aid in isolation. One form o f isolation network is called the radio frequency (rf) choke or transverse corrugation47*48'49. The rf choke consists of a series o f transverse (to the direction of propagation) one-quarter wavelength deep rectangular conducting grooves. The grooves act as transmission lines that transform the short at the bottom of the groove to a infinite impedance at the aperture. The infinite impedance at the aperture effectively prevents the flow o f any surface current transverse to the grooves. In such manner the rf choke can reduce the currents at the edges o f the biconical antenna and 8 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. reduce the radiation into the sideiobes. In Chapter 5, the effectiveness o f rf choke networks to reduce antenna sideiobes for biconical antennas is analyzed using the FDTD method and compared to measurements. Chapter 6 examines modifications to the basic Cartesian update equations for curvilinear, conformal and nonorthogonal coordinate systems. Cylindrical resonators, both with and without a dielectric insert, were analyzed using the contour path method,50 and the nonorthogonal method.51,52’53’54It was found that both o f these methods had an adverse affect on the stability of the FDTD algorithm, and the improvement in accuracy o f the calculated resonant frequency was not dramatically improved over the Cartesian grid. A discussion of conclusions and recommendations for future research is provided in Chapter 7. All o f the numerical results in this dissertation were obtained through the use o f FORTRAN computer programs. A brief description o f the main FDTD program is provided in Appendix A. The computer programs were compiled and run primarily on either a four processor Convex C-3840 with a 12 nsec clock and a gigabyte o f RAM, or a dual processor Cray YMPel II with 512 megabytes o f RAM and a 30 nsec clock. Both o f these computers are in some sense “super-computers” in that array operations are vectorized for faster execution. However, both computers were retired as o f June 30, 1997 and replaced with more modem machines. 9 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 2 Description of the Numerical Procedure Maxwell’s Equations In this dissertation, the field distributions throughout a region o f space are calculated using the Finite Difference Time Domain (FDTD) method.8*12Using FDTD, the electric and magnetic field components are calculated at specific locations on an underlying grid as functions o f time. In general, the FDTD grids can be Cartesian or curvilinear, structured or unstructured, orthogonal or non-orthogonal. The grid consists of a number of basic Yee cells on which the individual components of the electric and magnetic field vectors are situated. On a three dimensional Cartesian grid, the three components of the electric (E-) fields are usually oriented tangential to and in the middle of the edges o f a Yee cell, while the magnetic (H-) fields are oriented normal to and in the center of the faces of the Yee cell, as illustrated in Figure 2-1 . The finite difference equations can be written directly from the Maxwell’s equations in differential form. In a Cartesian coordinate system, the source-free curl equations for the electric and magnetic fields can be written as follows: fit SHy ^ fit 5HZ ^ fit fiz 3y 5EZ 3EX fix 5z ’ 3EX fiEy 3y fix’ 10 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2. 1) z T Az Ax x Figure 2-1. Orientation and position of electric and magnetic field components within Yee cell. In (2.1), E and H are the electric and magnetic fields, e is the electric permittivity, p is the magnetic permeability, and cr is the electric conductivity. E and H are functions of position and time, and e, p, and <r are functions o f position only. All materials are assumed to be Linear. Finite Difference Form of Maxwell’s Equations The following notation will be used for the finite difference equations. The spatial discretization o f the fields will be indicated using lowercase indices so that the electric field Exat grid point ((i+l/2)Ax, jAy, kAz) will be written simply as E*(ij,k), or as Ex(i+l/2 j,k) when it is necessary to emphasize half-integral spacing o f the fields. The temporal discretization of the fields will be indicated using superscript indices so 11 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. that the electric field Ex at time nAt will be written as (E*)n. The finite difference representation of the H-field curl equations in (2.1) are as follows: e(i + - ik )— a__— H T “ 0 + j . j + j , k) - H T “ (i + j , j - j , k) (2.2a) Ay 1. H“+I/2(i + i i k + i ) - H"+I/2(i + j, k 2J 2’ Az , E " '( i , j + ~ , k ) ~ E ; ( i , j + i , k ) E(i-j + 2 - k ) h i ------------ o (i,j-F ^ ,k ) U + j ’k)) 2 (2.2b) Az V j + i k + i) - H“ l/2(i, j - i k + i)'1 Ax 12 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. s(Uj + 2 ’k) At J_. 2J (E r'C i,j,k + ^ ) + E ;(i,j,k + i ) ) H " 1,2(i + - , j, k + - ) - H"+l/2(i - - j, k + i ) , (2.2c) , HjP'ftj + i k+ i)-H;(i,j + i k + i ) r t i , j + ^ , k + i ) ------------- 2------ 2__------------ 2------ 2 . 2 2 At 1 1 E ^ C i . j ■+j ,.:k ■+1) ■- E" (i, j + - , k) (2.2d) Az f E l+V\ i , j +1, k + h - E"+1/2(i, j, k + V Ay 1. 1 M(i + 2 » i k + 2^ I r V At e i +v\ i + 1, j, k + b - e *+v\ i, j, k +h (2.2e) Ax E r I,2(i + j , i, k + 1) - EJ*l/2(i + j , j, k) Az 13 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. , . 1 . 1 1 ! H ^ ( i + i j + ^ , k ) - H ; ( i + i j + i k) i + _ , j + —, k ) ----------— 2 J 2 At 1 . n + l/2 E? U\ i + i-, j + 1, k) - E"+wz(i + ^-, j, k) (2.2f) Ay 1 E“ '/i(i + l , j + j , k ) - E " ' /2Ci; j + i , k ) ' Ax In (2.2), the material properties o f the given problem are represented using stored values for e, p, and cr. Each cell can therefore have a different material property. A complex geometry can be constructed, in rectangular building block style, from many Yee cells each having (possibly) different material properties. Planar, multilayer dielectric materials are handled without difficulty. Curved interfaces between different materials are approximated by a staircased approximation to the original curved surface. Naturally, the finer the mesh, the better the approximation. For every field point on the FDTD grid, the curl equations are solved to provide explicit difference equations for updating the E-field at time (n+1 )At in terms o f the Efield at the previous time step and the surrounding H-fields at the previous half-time step. Explicit difference equations can also be written for updating the H-fields at time (n+l/2)At in terms of the H-field at the previous time step and the surrounding E-fields at the previous half-time step. The explicit update equations for the electric and magnetic fields are found by solving (2.2) as follows: 14 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. At H "+I/2(i - f - ~ , j + ~ , k ) - H “+,/2(i + I e(i + ^ - , j , k ) *c)> Ay <r(i + j , j , k ) At H "+,/2 (i + j , j, k + j ) - H "+I/2 ( i + J , j\ k - i ) Az (2.3a) e (U + ^ ,k ) 1 E"+I(i,j + i , k ) = 2 At , E J(i,j + i , k ) s(i,j + - , k ) a (i,j + - , k ) ------- — _+ ---------- _A ----At ' H r m (i, j + 1 e(i,j + j , k ) a (i,j + i , k ) At 2 ' + 2 ,J o (i,j + ^ ,k ) 2 k + i ) - H ^ C i, j + k- j)' Az 1 .. O 2 ’J’k + 2 ' Ax (2.3b) 15 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. e(i,j + ^ ,k ) 1 cr(i,j,k + ^ ) At -2— p E ^ a i k + i ) 1 e(i,j + - , k ) cy(i,j,k + - ) ^-----+.---------L At 1 e(i,j + | , k ) * V 2 ) / — At ' Hr“ ( y + _1 Ax (r(i,j,k + ^ ) + i. 2 t + - H r” ( u - f k+^ Ay (2.3c) H r I/2( i , j + p k + i ) = H r ' /2(i,j + p k + | ) + p r K U + J .k + j) (2.3d) | E ; c i , j + j , k + i ) - E ; ( i , j + Y " k;) E :( i,j+ l,k + i ) - E ; f t i k + |-) Az Ay V I J Hr,,2(i+ }, Ak + i ) / = Hr'/2(i + L j, k + i ) + ----- p j p • (iG + ^ i k + j ) [ E ;(i+ l,j,k + i ) - E ; ( i , j . k + i ) (2-3e) E ;(i + j , j , k + l ) - E ; ( i + ^ .j,k ) Ax Az v J v • J H ^ i + ~ ,j + p k ) = H r " ( i + ~ , j + L k) + -------p - j-----t*G + p j + j . k ) i ik ) Is (2.3f) E ;(i + l,j + i , k ) - E ; ( i , j + i k) Ax £ 1 E ;(i + i . j + U k ) - E ;( i + I > 16 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Advantages o f FDTD for solving Maxwell’s Equations Equations (2.3), along with a set o f prescribed initial values and boundary conditions, provide explicit finite difference formulae for the electric and magnetic field components at all points within the grid for all subsequent times. Explicit formulae offer several advantages from a computational point o f view. Since each field update is given explicitly in term o f the field variables from the previous time step, no simultaneous solution o f a set o f equations is required. In addition, all o f the field updates during a particular time step can be performed in parallel since can none o f the variables on the right hand side of (2.3) depend on any variable on the left hand side o f (2.3). Thus, (2.3) provides a set o f equations that are very well suited to computation using parallel processing techniques.. The finite difference equations offer several other advantages for solving electromagnetic problems. To the extent that the finite difference equations accurately represent the partial derivatives in (2.1), (2.3) provide an exact solution to Maxwell’s equations without any approximations. The use o f centered finite differences for the space and time derivatives ensure that (2.3) are accurate to 2nd order in the cell spacing. Thus, the departure o f the numerical solution to a given problem from the exact solution is proportional to the square of the cell spacing, and the numerical solution converges to the true solution as the square o f the cell spacing. Since Maxwell’s equations are solved in the time domain, the FDTD method provides a means to investigate transient electromagnetic problems. In addition, the 17 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FDTD method provides an efficient means to obtain a wide frequency response o f a structure from a single simulation through the Fourier transform o f the time domain response. Sources In the FDTD method, a time dependent electric or magnetic source is introduced by setting the electric field components at one or more cells equal to a function o f time. The particular time dependence o f the source used is not critical, but the bandwidth o f the applied pulse must be sufficient to include the frequency range o f interest. Often, a Gaussian time dependence is used. The time, dependent fields are then updated in a leap frog manner where the H-fields throughout the mesh are updated during the first halftime step, then the E-fields are updated throughout the mesh during the second half-time step. The input pulse is usually applied until the source decays to some small nominal value, then the source is turned off and the fields at the source locations are updated for the remainder of the time using the source-free finite difference equations. Boundary conditions Perfect Electrical Conductors (PEC) and Perfect Magnetic Conductors (PMC) Boundaries that are coincident with perfect electrical conductors (PEC) require the vanishing of tangential electric fields and normal magnetic fields. In the FDTD, PEC boundaries are most easily enforced on the faces of the basic Yee cells. The vanishing of the tangential electric fields and normal magnetic fields are accomplished by setting the fields on the boundary grid to zero. The null fields on the boundary are then excluded from the normal update equations. 18 ! Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Boundaries that are coincident with perfect magnetic conductors (PMC) require the vanishing o f tangential magnetic fields and normal electric fields. In the FDTD, PMC conditions on a boundary is easily enforced on a plane mid-way between and parallel to the opposite faces o f a Yee cell. The vanishing o f the tangential magnetic fields and normal electric fields are accomplished by setting the fields on the boundary which passes through the center of a Yee cell to zero. The null fields on the boundary are then excluded from the normal update equations. Occasionally it is inconvenient to place the a PMC boundary on a plane that passes through the center of a Yee cell. For example, for large problems with planes of symmetry, it may be convenient to divide the size of the mesh in half. The FDTD equations are then solved twice, once with PEC conditions on the symmetry plane and once with PMC conditions on the symmetry plane. Assuming that the plane of symmetry is coincident with the face of Yee cell, the update equations for the tangential E-fields on the boundary must be included in the PMC case. Equation (2.3) can be used for updating the normal H-fields on the boundary since only fields on the boundary are involved. A new update equation is required for the tangential E-fields on the boundary since values of the magnetic fields on both sides o f the boundary are required. The new update equation assumes that the tangential magnetic fields one half-cell on either side of the PMC boundary are equal in magnitude but opposite in sign (odd). For example if a PMC boundary is applied at x=0, then the new update equations for Ey and Ez are as follows: 19 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. e (0 ,j + ^ rk ) E"+l(0,j + ^ ,k ) = 2' CT(0,j + ^ ,k ) At E " y(0,j+i,k) 1 e (0 ,j + - , k ) o(0,j + - , k ) At Hrra(0,j +|,k +i)- H r 2(0,j + i , k - i ) ' 1 e ( 0 , j + i k) At Az o (0 ,j + i k) ----- + 1 . 1 2 .H r z(0 + -,j + - ,k ) Ax (2.4a) 1 E r i( o ,j ,k + - ) = e(0,j + i k) <r(0,j,k + | ) At 2 e(0,j + ^ ,k ) o-(0,j,k + | ) At * E "(0,j,k + ^ ) 2 1 e(0,j + ^ ,k ) At <x(0,j,k + i) + -------------— O 2 • H"+I/2(0 + ~ , j, k + ^ ) 2’ Ax (2.4b) ) ' H "+l/2(0, j + 1 , k + 1) - H ^ C O , j - ir1,,kku + -~1 ) ^ 2 2 Ay Open (absorbing) Boundary Condition In some problems it is necessary to truncate the grid with an open or absorbing boundary condition (ABC). Open boundaries occur in antenna problems where radiation to infinity is assumed and in transmission line problems where one or more of 20 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the transmission lines extend to infinity. Many different ABCs have been developed, and the degree to which the ABC is reflectionless depends both on the particular ABC and the nature o f the electromagnetic fields in the vicinity of the boundary. In any case, consideration must be given both to the numerical complexity and cost o f applying the ABC and the degree to which it absorbs the outgoing wave. One o f the most often used and one o f the first to be successfully applied to FDTD problems is the Mur ABC.16 Mur based his method on the concept o f an annihilation operator for locally absorbing an outgoing wave which was described in detail by Engquist and Majda55. If W(x,y,z,t) is used to represent an outgoing wave, then Engquist and Majda’s 1st order approximation for an absorbing boundary at x=0 is (2.5) where c is the speed o f wave propagation. One o f the major advantages of (2.5) is that it can be applied locally, i.e., that it only involves fields in the immediate vicinity o f the boundary, and therefore it can be easily and efficiently applied to the FDTD method. Following Mur’s discretization procedure for (2.5) with W(x,y,z,t) replaced by Ez(ij,k+l/2) and Ey(ij+l/2,k), the ABC for the tangential electric fields at x=iAx=0 can be written as: E " '( 0 ,j + i k) = E ; ( l , j + i k) + c A t-A x cAt +A x E r '( o , j , i c + j ) = E ; a j , k + i ) + c A t-A x cAt +Ax E r ' a j . k + ^ - E J C O . i k + i ) 21 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (26b) At the boundary x=imaxAx, the ABC for the tangential electric fields are: I -1 , j + 7 ,k ) 2’ cAt - Ax E f O ™ -1 , j + cAt +Ax (2.7a) k ) - E J C w j + i , k) 2’ K +i (i max, i k + t ) = E " ( i max -1 , j, k + ~) (2.7b) cAt - Ax cAt +Ax At the boundary y=jAy=0, the ABC for the tangential electric fields are: l 1 cA t-A y Ex+l(i + -r A k ) = E"(i + -r,l,k) + cAt + Ay E r'(i,0 ,k + i ) = E ;(U ,k + i ) + ^ Er'(i + ^ , U ) - E ; ( i + i o , k ) E“+l(U . k+ r-) - E"(i,0, k + 4-) (2.8a) (2.8b) At the boundary y=ymaxAy, the ABC for the tangential electric are: E r '(i W .k ) = EJ(i t j . u - l,k) cAt - Ay E r ' ( i 4 . i „ - U ) - E " ( i + i Jlml,k) cAt + Ay (2.9a) E r 'a jm « ,k + i ) = Ez" (i,j™, - l . k + ^ ) cAt - Ay 1 1. E ^ a jm a x - U k + ^ ) - E " ( i , j max,k + ^ ) cAt + Ay (2.9b) At the boundary z=kAz=0, the ABC for the tangential electric fields can be written: cA t-A z E"+1(i + Tj>0) = E"(i + r-,j,l) + 2 cAt +Az Er'(i + ^ , j , l ) - E ; ( i + i 22 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. j,0) (2.10a) 1 tx cAt - Az E r f t j + ^ 0 ) = E ; ( i , j + ^- ,! ,I)) + ----------2 cAt + Az 1 2 (2.10b) At the boundary z=kmaxAz, the ABC for the tangential electric fields can be written: E"+l(» + J,kmax) = E"(i + j , j,k max -1 ) + c A t- Az cAt-t-Az ■ (2.1la) + \ , j. K . - 1) - E ! ( i + j , i , k m ) 1 2 (2.11b) cA t-A z cAt +Az Note that (2.5) completely annihilates a plane wave W(x+ct)moving in the negative x-direction. Waves traveling off axis will not be completely absorbed. Higher order approximations have been introduced that are more successful at absorbing offaxis waves, but these in general involve more complicated expressions that require storage of the fields in the vicinity of the boundary for several time steps. The formula in (2.6)-(2.11) have the advantage that they are simple to apply, require little additional computer memory, and they have little or no effect on the stability of the FDTD algorithm. A straightforward procedure for improving the performance of any ABC used with the FDTD technique was introduced by Mei.20 In this procedure, the same ABC is applied to the tangential H-fields one-half cell inside the boundary where the ABC was applied to the tangential E-fields. The technique is based on the observation that the error in any ABC can be reduced by applying the same ABC to both the tangential magnetic and electric fields, and then re-calculating the tangential H-fields. The error is 23 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. defined as the difference between the field that would exist on the boundary in the absence of any reflected field and that which exists due to incomplete absorption. The superabsorption technique is most effective when the angle o f incidence is near zero degrees. This renders the technique most suited to absorption o f propagating modes in transmission line problems where the mode of propagation is either TEM or pseudoTEM as discussed in Chapter 4. To apply the superabsorption technique at x=iAx=0, Mur’s 1st order ABC shown in (2.6) is applied to the tangential E-fields. Then the tangential H-fields are calculated twice: first, using the normal FDTD update equations (2.3) to obtain H*1*, and second, by applying the ABC to the H-fields to obtain H(2). At x=(0+l/2)Ax, H(2) is given by: H " ra<2)( ^ ,j,k + i ) = H ;-'/2( | , j , k + C A t-A xf n+1/2 3 cAt + Ax I\ 2 TTn-l/2,1 i |c , I \ y 2 ,J’ 2 cAt —Ax cAt + Ax + ------------/»Af i A An improved value for Hy and Hz at x=(l/2)Ax is then calculated from a weighted average of Hyn+I/2(I), Hyn+I/2(2), and HzrtH/2(l), HzIt+l/2(2), respectively, as follows: irn+vis y j,k H -i)+ PH ;+l/2(2>(i j , k + i ) ; i_ ,_\ __________2_____ 2___________ 2_____ 2_ 22 l+p K +U2W(~ J + (2.13a) k) +.p H r l/2(2)( - , j + r , k) (2.13b) H r ,/2(T>j + T ’k) 24 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where p=cAt/Ax. The new values Hy and Hz in (2.13) contain less reflection error than that given by (2.12). A further improvement can then be obtained using the new values for the Hfields to refine the value o f the tangential E-fields at x=0 by inverting the update equations for Hy and Hz in (2.3) as shown below: (2.13a) H"-"2( - , j, k + - ) - H"+I/2(—, j, k + - ) E "(0,j + |- ,k ) = E "(l,j + | , k ) (2.13b) The use o f the superabsorption reduces the boundary reflection by a factor of 4 to 10 near normal incident angles. For the stripline analysis of Chapter 0, the use o f a superabsorbing boundary to terminate the striplines results in negligible reflection. In addition, the FDTD algorithms that include superabsorption have proven to stable for at least 2 IS time steps. Stability limitations In most cases, it is desirable to choose the time step parameter At as large as possible to reduce the number of calculations. However, At cannot be chosen arbitrarily since too large a value results in a numerical algorithm that is unstable, i.e., unphysical 25 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. spurious solutions to the discretized ME lead to unbounded exponential growth o f the fields. In other words, the numerical solution blows up. The stable region can be determined from an eigenvalue analysis of the discretized curl equations where the allowed eigenvalues are chosen to ensure that the fields calculated from the update equations do not grow unbounded with time.56 According to such an analysis, the time step is governed by the following equation: 1 (2.14) At = c >/3c when Ax = Ay = Az. J In (2.14), the time increment for a cubical lattice is expressed in terms o f the grid spacing Ax, the speed of light c, and a parameter s called the stability parameter. This last equation states that the time step must be less than about 0.577 multiplied by the time it takes light to propagate one ceil spacing in a cubical lattice. It expresses one of the fundamental constraints of the FDTD method, i.e., that the value o f the stability factor s in three-dimensional cubic lattices must be less than about 0.577. When structures with high quality factors are involved, a great many time steps may be required if good frequency resolution is to be obtained. The relationship between the frequency resolution obtained from the discrete Fourier transform o f a time dependent signal, the sampling interval At and the number o f time samples N is: (2.15) 26 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Suppose, for example, that a frequency resolution of Af=10 MHz is required at a center frequency of 10 GHz. According to (2.15), At should be as large as possible in order to minimize the number of time steps for a given frequency resolution. However, At is limited by (2.14). If the maximum value o f At from (2.14) is used in (2.15), and a cell spacing o f X/20 is assumed, then the required number o f time steps N is: 34,600. (2.16) Thus, 34,600 time steps of the FDTD algorithm are required for a frequency resolution o f 10 MHz at 10 GHz, assuming a cubical cell spacing o f A/20 and the maximum allowed time step. According to (2.16), the number o f time steps required is inversely proportional to the cell spacing and the frequency resolution. Numerical dispersion The grid spacing in relation to the wavelength of the propagating waveform affects the numerical phase velocity. A propagating wave in the FDTD grid will generally have a phase velocity which is slightly less than the true phase velocity. The change in phase velocity that occurs due to the finite cell size is referred to here as numerical dispersion. Tafloves6 shows a simple analysis to determine the change in numerical phase velocity as a function o f the grid spacing. Taflove assumes a vector* field traveling wave expression of the form: 27 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.17) Px=Px, Py =P-y, Px = P-z, where the italicized j in the exponent in (2.17) represents > /-f, V* represents a linear combination o f the vector electric and magnetic fields, and px, Py, and pz, are the numerical wave numbers along each of the axes. Taflove inserted (2.17) into a space time central differencing realization of ME, similar to (2.3), to obtain the following numerical dispersion relationship: 1 . f ©At cAtSml~ 2 2 2 . fM y )' + A y sm [ 2 ’ i + J _ . f PzAz^j AzSml 2 J . (2.18) The numerical phase velocity (ca/p) determined using (2.18) is always less than the speed of light c. In a cubical lattice, the lowest vPh occurs when the direction o f propagation is along one of the coordinate axes, and is a maximum along a diagonal o f the lattice cell. When Ax, Ay, and Az are all unequal, then the slowest phase velocity (worst case) occurs along the axis having the largest grid spacing. Assuming that the wave propagation is in the x-direction, in which case Py=pz=0, (2.18) can be solved for Px: n 2 . Ax . f o A t^ p* = ^ sm I s - v t JJ- (2.19) When Ax is larger than Ay or Az, (2.19) then provides a worst case assessment for the numerical dispersion. The numerical phase velocity vPh normalized by the speed of light is given by the following expression: 28 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 2-2 plots the normalized phase velocity as a function o f the wavelength divided by the cell spacing Ax using (2.19) in (2.20). A stability parameter s=0.5 was assumed to calculate the time increment At using (14). The numerical phase velocity is seen to be about 1.3% less that c when the 10 cells per wavelength are used. The phase velocity increases to just 0.1% less than c when about 36 cells per wavelength are used. 0.99 0.98 0.96 0 .95 0 .9 4 0 10 20 30 40 Cell Dimension (A/Ax) Figure 2-2. Numerical phase velocity versus grid spacing. 29 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 3 Application of the FDTD Method to Analysis of Resonant Cavities Calculation o f Resonant frequencies using the Fourier Transform The FDTD method can be applied in a straightforward manner to find the resonant frequencies (or eigenvalues) and field distributions within a PEC enclosure or cavity. PEC boundary conditions are easily enforced on the outer boundaries by setting the tangential E-fields and normal H-fields to zero, as described in the last chapter. The fields elsewhere in the computational domain are calculated using the usual update equations (2.3) for each time interval. The time-dependent fields are updated in a leap frog manner where the H-fields throughout the mesh are updated during the first halftime step, then the E-fields are updated throughout the mesh during the second half-time step. The resonant frequencies are determined by recording the time series o f the field components at one or more locations within cavity, and then calculating the discreet Fourier transform o f time series, which is defined by: Ez(i, j,k;© ) = 2 Ez(i, j, k; nAt)e"/mnAt. (3.1) 11=1 (3.1) can be calculated using any o f the readily available discrete or fast Fourier transforms (DFT or FFT). In order to avoid aliasing, the sample interval nAt must be chosen less than one-half the period of the highest desired frequency component. In Chapter 2 it was shown that due to dispersion considerations, at least ten cells per wavelength are necessary to limit dispersion to 1.3% or less. When using a time step 30 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. stability factor of 0.5 with ten cells per wavelength at the highest frequency o f interest, the sample period will be about one-twentieth (1/20) of the period at the highest frequency o f interest. Under such conditions, adequate sampling rates are ensured for the calculation of resonant frequencies. Generally, the amplitude and phase o f the Fourier transformed fields are functions o f the location, direction, and bandwidth of the applied source. In the FDTD method, a time dependent electric source is introduced into the cavity by setting the electric field or magnetic field components at one or more locations equal to a function of time. The particular time dependence o f the source used is not critical, but the bandwidth o f the applied pulse must be sufficient to excite the resonator modes of interest. For the calculations here, a sinusoidally-modulated Gaussian time dependence was typically used. For example, the following time equation can be applied to the Ezcomponent when excitation of the transverse magnetic (TM) modes with respect to the z-direction is desired: where n ^ , specify the time step at which the pulse amplitude peaks and the width of the applied pulse, respectively, and g j 0 specifies the center frequency o f the applied pulse. If only transverse electric modes are desired, then an equation similar to (3.2) can be applied to H*. Generally, to avoid introducing undesirable high frequency spectral components, no>4n, should be chosen. (3.2) represents an applied source with a frequency bandwidth of approximately l/(2n,At), centered at <o0, with a small 31 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (negligible) dc component. The location (ij,k) of the E* field to which the source is applied is not critical. Resonant Frequencies for a Rectangular PEC Cavity The resonant frequencies below 20 GHz within an empty 0.7 in. x 0.8 in. x 0.9 in. rectangular cavity are shown in Table 3-1. The theoretical frequencies for the cavity were determined using the formula: f « =5901.4. (3.3) where fres is the resonant frequency in MHz, a,b,d are the cavity dimensions in inches, and m,n,p are the mode numbers. The resonant frequencies from the FDTD results were found using the Fourier transform of the sampled time domain fields. The time domain fields were calculated using a Cartesian FDTD Grid with 14x16x18 cells and a cell dimension(Ax=Ay=Az) of 0.05 inches. Using a time step o f 2.1 psecs, the fields at selected locations within the cavity were sampled for 131,072 time steps, which when transformed to the frequency domain provided a frequency resolution o f 3.6 MHz. The processing time was about 655 cpu secs on a CONVEX C-3840 super-computer with a 12 nsec clock. The FDTD results for the first five lowest order modes agreed with the theoretical fres in Table 3-1 to within 10 MHz (0.1%). The error increased to -0.26%, -0.24%, -0.37% and -0.31% for the TEou, TE 102, TE021 and TM 121 modes, respectively. The increased error at the higher mode numbers and higher frequencies is attributed to 32 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. numerical dispersion, as discussed in Chapter 2. Due to numerical dispersion, the waves travel more slowly through the FDTD grid at the higher frequencies, which leads to a decrease in the resonant frequency. Table 3-1 Resonant frequencies for a 0.7 in. x 0.8 in. x 0.9 in. Rectangular Cavity frcs FDTD (MHz) frcs Error Mhz(%) 9,870 10,680 11,202 12,980 12,980 15,047 15,590 16,145 16,992 17,247 17,247 18,091 18,214 18,214 18,404 19,537 19,537 19,740 9,862 10,670 11,192 12,974 12,974 15,008 15,553 16,086 16,939 17,221 17,221 18,000 18,170 18,170 18,316 19,463 19,463 19,681 -10 (-0.09) -10 (-0.09) -6 (-0.05) -6 (-0.05) -39 (-0.26) -37 (-0.24) -59 (-0.37) -53 (-0.31) -26 (-0.15) -26 (-0.15) -91 (-0.5) -44 (-0.24) -44 (-0.24) -88 (-0.48) -74 (-0.38) -74 (-0.38) -59 (-0.3) . • 0 0 00 frcs Theoretical (MHz) 1 00 Resonant Mode: TEm^p or TMm.n,D TEoii TEioi TMno TEiii TMni TEoiz f t [02 TE02l TM,20 TE,|2 TM„2 TE2oi: TE,2I TM,21 TM210 TE2H TM211 TE022 Determination of Quality Factor When losses are small, the quality factor, or Q, of a high frequency enclosure can be estimated by assuming that the losses have negligible effect on the distribution o f field energy in the cavity. The Q of an enclosure at a radian frequency coo is defined by: Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. CD0(stored energy) average power loss g>0U e (3.4) PL 1 1 The stored energy UE is given by the volume integral of ^eE2 (or ^ H 2) over the volume o f the enclosure V, where the electric field E (or magnetic field H) is the amplitude o f the field in the enclosure. For an empty enclosure, the losses will be due to conduction currents in the walls of the enclosure. When the walls are constructed from good conductors, then to a very good approximation the average losses can be estimated from the surface integral o f over the interior surface o f the enclosure S, where R* is the surface resistance o f the conductor in ohms per square and Ht is the amplitude o f the tangential magnetic field at the wall of the enclosure. Assuming a empty enclosure, (3.4) becomes: (3.5) The values of eE2 andRgHt in (3.5) must be evaluated at the frequency o)q. For the purpose o f calculating the Q in low loss enclosures, the values of E or H in the integrals (3.5) will be assumed equal to their values in a lossless enclosure, which is a good approximation if the losses are not too large (high Q). In the numerical calculations, the integrations in (3.5) are replaced by summations throughout the volume of the enclosure and by summations over the surface of the conductors. In the FDTD method, the fields throughout the grid are specified as functions o f time. In order to evaluate the Q according to (3.5), the square of the amplitude o f the 34 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. electric (or magnetic) fields throughout the volume, and the square o f the amplitude of the magnetic field across the interior conductive surfaces must be determined at a frequency too- The FDTD update equations provide the total electric or magnetic field as a function of time, which usually will include many different frequencies depending on the modal distribution o f the enclosure and frequency content and spatial distribution o f the source. The time domain fields must be converted to a function o f frequency, which can be done using a FFT as discussed in “Calculation o f Resonant frequencies using the Fourier Transform” o f this chapter. The field amplitude at a particular frequency can then be squared and used in (3.5). From a practical point of view, the evaluation o f the integrals in (3.5) in the frequency domain becomes a bit unwieldy since the time series o f the electric fields at every spatial point throughout the enclosure need to be recorded. Then, each time series at every spatial point is Fourier-transformed. The amplitude o f the Fourier component at the frequency coo is squared and used in the volume and surface integrals in (3.5) to determine Q. A second, less memory-intensive alternative for calculating the enclosure Q is to carry out the summations in the time domain prior to Fourier-transforming. In this case, only the volume summation o f £E2 or the surface summation ofR^Hf are recorded at each time step for the duration of the run. This approach requires much less computer memory than storing the individual time domain E- and H-fields throughout the volume and across the interior surfaces. Now only the values of the volume summation o f sE2 and of the surface summation o f RjFlf at each time step are Fourier-transformed to 35 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. determine the amplitudes o f the total stored energy and o f the total conductive loss as functions o f frequency. Since the field quantities to be Fourier-transformed are squared, the Fourier amplitudes at a frequency 2-cd<j are used in the calculation o f the resonator Q. The quality factors for each resonant mode listed in Table 3-1 is listed in Table 3-2. The quality factors are shown both from the FDTD results and from the analytical theory. A value o f 0.0261 ohms per square was used for R*. The FDTD results were obtained by integrating eE2 and RgH? at each time step for 131,072 time steps, then transforming the time series to the frequency domain. The values o f £E2 and at the frequency 2-©o was used in (3.5), where ©o is the resonant frequency for the particular mode. The same Cartesian FDTD Grid with 14x16x18 cells and a cell dimension (Ax=Ay=Az) o f 0.05 inches was used for Tables 3-1 and 3-2. The agreement between the theoretical Q and the Q caluclated using FDTD varied between 0.12-3.5%. Some variance was expected since the cavity was excited at a single point, and more than one mode was always present. For more precise calculations, a source that approximates the expected field distribution inside o f the cavity could be used. In that case, more energy would be placed in the mode o f interest, and there would be less noise in the calculated field quantities. 36 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 3-2 Quality factors for a 0.7 in. x 0.8 in. x 0.9 in. cavity Mode (TEni.n.D Of TMnui.o) TEou TEiot TMno TEui TMni TE012 TE 102 IE 021 TM,20 TE||2 TM i ,2 TE201 TEI2 I TM 121 TM2 I0 TM211 TE022 Q fd td Q th eo ry 10,075 10,689 11,091 9,635 10,146 15,751 16,609 16,415 17,992 13,243 13,711 17,895 13,472 14,822 18,151 15,318 20,324 9,970 10,617 11,294 9,534 10,134 15,561 16,387 15,987 17,662 13,055 13,466 17,286 13,188 14,558 18,093 14,957 19,941 Error (%) 1.05 0.68 -1.80 1.06 0.12 1.22 1.35 2.68 1.87 1.44 1.82 3.52 2.15 1.81 0.32 2.41 1.92 Analysis o f Cavity with Superconductive thin films for determination of surface resistivity and comparison with measurement The parallel plate resonator method has become a popular means to determine the surface resistance o f superconductors at microwave frequencies.24 Here, the Finite Difference Time Domain (FDTD) method is used to analyze the resonant frequency and the dielectric and conductor losses within a microwave enclosure containing a parallel plate resonator. The resonator’s losses are examined for the purpose o f identifying experimental conditions that could lead to inaccurate surface resistance measurements. The parallel plate resonator technique for the measurement o f the surface resistance o f superconductors is described by Taber.24 A diagram o f the test enclosure Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. is shown in Figure 3-1. In this method, a parallel plate resonator is formed from two samples o f a superconductive material separated by a thin dielectric spacer. The parallel plate resonator is placed in a slightly larger test enclosure. No electrical connections are made to the resonator, but excitation o f the transverse electromagnetic modes between the parallel plates are accomplished by placing small electrical probes near to the open edges o f the resonator. The electromagnetic fields within the test enclosure are dominated by the transverse electromagnetic modes between the plates o f the resonator as a result o f the close spacing between the plates. In the ideal case, the Q o f the enclosure is dominated by the superconductive losses o f the parallel plate resonator. The exact placement o f the resonator within the larger enclosure is not critical. The method has become popular principally because it does not require detailed knowledge o f the modal distributions within the resonator and because sample preparation is relatively simple. Other methods such as the cavity end wall replacement method often suffer from large experimental error. In the end wall replacement method, particular cavity modes are measured with all normal conducting materials and with one end wall replaced by superconductive material. However, the loss in the super-conductor must be found by subtracting two measurements in which differences between the normal cavity and the cavity with one superconductive end wall are relatively small. 38 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.5 in. View Superconductive thin film . L aA lO , substrate Cover Sapphire spacer Figure 3-1. Test fixture for the measurement o f surface resistance. As described by Taber,24 the three most significant losses that determine the Q o f the parallel plate resonator a re : (1) the resistive loss due to the surface resistance Rssc o f the super-conductor within the parallel plate resonator, (2) the dielectric losses inside the dielectric spacer and substrate materials, and (3) the radiation loss at the edges o f the open resonator. Taber assumes the unloaded Q o f the enclosure is given by: q -i + a s + ta n 5 , (3 6 ) s where s is the resonator spacing, tan5 is the dielectric loss tangent, and a and (3 are proportionality constants that depend on the radiation loss and the modal structure o f the resonant cavity, respectively. Taber shows that p is independent o f the mode number for the open parallel plate resonator and is equal to In Taber’s experimental procedure, the surface resistance of the super-conductor Rs is found by measuring the 39 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. unloaded Q for a small resonator spacing s (typically 0.001 in. or less) and assuming the losses are dominated by the surface resistance of the super-conductor Rs . In such a case, Rs is equal to s/Qfi. The field distributions and cavity losses within the microwave enclosure and parallel plate resonator will be examined using the Finite Difference Time Domain (FDTD) method. Using FDTD, the electric and magnetic field components are calculated throughout the enclosure in the time domain on a Cartesian mesh. Due to the very small spacing used between the superconductive surfaces o f the parallel plate resonator (about 0.001 in.), a different cell dimension is used perpendicular to the plane of the plates (z-dimension) than was used in the tangential directions (x- and ydimensions). In the FDTD method, a time dependent electric source is introduced into the cavity by setting the electric field components at one or more cells equal to a function o f time. The particular time dependence o f the source used is not critical, but the bandwidth o f the applied pulse will be chosen to excite the resonator modes o f interest. For the parallel plate resonator, a Gaussian time dependence was used. The input pulse was applied until the source decayed to some small nominal value, then the source is turned off and the fields at the source locations are updated for the remainder o f the run in the usual manner. The resonant frequency was determined by recording the time series o f the field components at one or more locations within the parallel plate resonator, and then calculating the Fourier transform o f time series. Generally, the amplitude (and phase) o f 40 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the frequency components o f the Fourier transform is a function o f the locations o f the field sample points and o f the source, and the bandwidth o f the applied source. The sample locations are usually chosen with some knowledge concerning the distribution o f the field energy. It is possible to select locations for the field sampling where the amplitude o f the particular resonance is zero, for example, at a node. However, by sampling the fields at a number o f locations, the chance o f missing a particular resonance is minimized. The losses within the normal conductors in the walls o f the enclosure, the dielectric spacer, the dielectric substrate, and the superconductive thin films were calculated. The conductor losses within the enclosure were determined for a given resonator spacing by integrating the square o f the time-averaged tangential H-field over the surface and multiplying by R*. The dielectric losses were calculated by multiplying the square o f the time-averaged E-field by soo-tand, where e is the dielectric permeability, ©o is the resonant frequency and tan5 is the dielectric loss tangent. The cavity Q was then calculated using JJJeE2dv (3.7) where in the numerical calculations the volume integrals are replaced by summations throughout the volume o f the enclosure and by summations over the surface of the conductors. The summation o f the losses over a surface or volume was carried out in the time domain since less computation time and computer memory is required than if the 41 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. time domain fields throughout the enclosure were transformed to the frequency domain and then summed. When performed in the time domain, the summation o f eE^ through the resonator volume or the summation o f Rs H2 over the conductive surfaces are recorded at each time step for the duration o f the run, then the time series o f the surface and volume summations are Fourier transformed to determine the amplitude o f the field energy and the amplitude o f the conductor loss as functions o f the frequency. The timeaveraged quantities are then given by ^ eE2(©) and ^RsHf(ca), respectively, where E2(©) and Ht(co) represent the amplitude o f the Fourier-transformed squared fields at frequency ©. The amplitudes o f the squared fields at a frequency o f 2-©o are used to calculate the resonator Q at frequency ©o. If the total field energy is assumed concentrated in the parallel plate resonator and the dielectric losses and radiative losses are small, then solving (3.6) for the surface resistance o f the super-conductor yields RSs c = s/PQ, where Q is measured for the entire enclosure. However, the LaAlC>3 substrate affects the proportion o f energy located outside o f the parallel plate resonator, and therefore affects the effective Q. Fields located outside o f the parallel plate resonator are dissipated by losses in the LaAlOj substrate and by the surface resistance o f the normal conductor RsCthat forms the surrounding enclosure. When the fields outride the parallel plate resonator are included explicitly in the loss calculations, then (3.7) can be used to re-write (3.6) in the following form: 42 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. oo0JJJeE2dv ee0JJJeE2dv 1 resonator V en clo su re JJJgE2dv (3.8) 1 V dielectric In the first term on the right hand side o f (3.8), the surface integral in the numerator extends over the surface o f the super-conductor with surface resistance given by R^,.. The surface integral in the second term on the right hand side o f (3.8) extends over the inside surfaces o f the enclosure with surface resistance given by R*.. Each volume integral in (3.8) extends over the entire volume o f the enclosure, including that o f the parallel plate resonator. The calculation of the effective Q by FDTD follows the prescription given by (3.8). Figure 3-2 plots the value o f s/fiQ calculated using FDTD versus the resonator spacing s for the fundamental mode o f a 0.5 in. x. 0.25 in. parallel plate resonator. The parallel plate resonator consists o f two 0.5 in. x. 0.25 in. superconductive thin films (Rssc=0.74xl0'3 ohms per square at 10 GHz and 77°K) separated by a sapphire dielectric spacer (e=9.2, tan8=2xl0‘7) o f thickness s. The resonator is centered in a 0.75 in by 0.5 in.by (0.040+s) in. thick, copper enclosure (RyC=0.015 ohms at 10 GHz and 77°K). The calculated data in Figure 3-2 is shown with and without a 0.020 in. thick LaAlOj substrate (e=22, tan6=l0'6) for the superconductive films. The values o f s/PQ are shown for four different sapphire spacer thicknesses o f 0.001 in., 0.002 in., 0.004 in. and 0.008 in. The curves shown are quadratic curve fits to the four data points. 43 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 3-1 shows that the dielectric substrate significantly affects the shape o f the s/pQ versus s curves. The value o f s/PQ decreases much more rapidly with the spacer thickness s with a LaAlCb substrate than without a LaAlCb substrate. Ifs/PQ is used to estimate the surface resistance Rs, then the values o f Rs derived from the curve are substantially higher than values derived from the curve without LaAlOj. The discrepancy between the two curves is due to the greater field energy stored in the LaAlOj substrate which leads to larger losses in the normally conducting walls o f the enclosure. 7.00E-04 6.00E-04 5.00E-04 I 4.00E-04 § 3.00E-04 -C o *0 2.00E-04 : R^c without LaA103 1.00E-04 O.OOE+OO 0 0.001 0.002 0.003 0.004 0.005 Spacer Thickness (mils) 0.006 0.007 0.008 Figure 3-2. Plot of s/fiQ vs. spacer thickness s. Table 3-3 and 3-4 list the resonant frequencies o f the dominant mode calculated by FDTD for the four sapphire spacer thicknesses of Figure 3-2, along with the total 44 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. effective Q and the separate Qs due to the parallel plate resonator only (super-conductor losses only), the enclosure (normal conductor losses only), and the dielectric losses only, for the resonator with and without the LaAlOj substrate, respectively. As indicated in (3.8), the reciprocal o f the effective Q is equal to the sum o f the reciprocal Qs due to each o f the individual loss mechanisms acting alone. The losses in the walls o f the enclosure in Tables 3-3 and 3-4 are further hroken down according to losses in the side walls and losses in the top and bottom walls o f the enclosure. The Q for the side walls is calculated by dividing the volume integral o f the E-field amplitude throughout the enclosure by the conductive losses in the four side walls perpendicular to the plane o f the parallel plate resonator. The Q for the top and bottom walls is calculated by dividing the volume integral o f the E-field amplitude by the conductive losses in the top and bottom walls parallel to the plane o f the parallel plate resonator. In all cases, the losses in the side walls are much smaller (i.e. higher Q) than the losses in the top and bottom walls. The Q due to losses in all six walls is equal to one divided by the sum o f the reciprocal Qs for the side walls and for the top and bottom walls. Table 3-3. FDTD results for Q,Resonant Frequency and Surface Resistance o f the Dominant Mode in a Parallel Plate Resonator with a LaAlOj substrate. Spacer thickness in inches: Resonant frequency in GHz Q due to losses in walls of enclosure: In Side Walls: In Top and Bottom Walls: Q due to super-conductor losses: Q due to dielectric losses: QefTcctive Rssc= s/ PQ with LaA1 0 3 (ohms) Rssc actual (ohms) 0.008 3.434 4379 259394 4454 46975 9592890 4004 6.88E-04 8.73E-05 0.004 3.601 5686 246305 5820 20820 7896328 4464 3.24E-04 9.60E-05 0.002 3.727 9134 341488 9385 9341 6722442 4615 1.62E-04 1.03E-04 45 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.001 3.809 12187 427970 12544 4186 6054591 3114 1.23E-04 1.07E-04 Table 3-4. FDTD results for Q, Resonant Frequency and Surface Resistance of the Dominant Mode in a Parallel Plate Resonator without a LaAlCb substrate. Spacer thickness in inches: Resonant frequency in GHz Q due to losses in walls of enclosure: In Side Walls: In Top and Bottom Walls: Q due to super-conductor losses: Q due to dielectric losses: Qefftctive Rssc= s/pQ without LaAlCb (ohms) Rssc actual (ohms) 0.008 4.194 15334 723572 15666 29707 5092450 10094 3.33E-04 1.30E-04 0.004 4.064 26538 1373469 27061 14437 5049447 9333 1.75E-04 1.22E-04 0.002 3.995 50468 2696391 51431 7136 5034216 6244 1.28E-04 1.18E-04 0.001 3.956 99691 5429076 101556 3586 5029408 3459 1.15E-04 1.16E-04 The value o f the superconductive surface resistance, calculated using Rssc= j/pQ, where Q is the effective Q, and the actual value o f the super-conductive surface resistance Rssc are also shown in Tables 3-3 and 3-4. The value for Rssc actual is based on the assumed value o f 0.74xl0~3 ohms per square for a high temperature superconductor (HTSC) at 10 GHz and 77°K that was adjusted for frequency according to the formula: R-ssc = 0-74xl0-3^j—] , (3.9) w here/is the frequency in gigahertz. The values for Rssc used to calculate the losses in the super-conductor in Tables 3-3 and 3-4 were also adjusted for frequency using (3.9). The value for Rsc used to calculate the losses in the walls o f the enclosure was based on a typical value for copper of Rsc=0.015 ohms per square at 10 GHz and 77°K that was adjusted for frequency according to: 46 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Rsc = 0.015A/ £ , / i n gigahertz. (3.10) In Tables 3-3 and 3-4, the frequency used to specify the value o f Rsc and Rssc in the calculations was the resonant frequency o f the dominant mode in each case. The resonant frequency was different for each spacer thickness. The resonant frequency o f the dominant mode in the resonator increases with decreasing spacer thickness with the LaAlCh substrate, while the dominant mode frequency decreases with decreasing spacer thickness without the LaAlC>3 substrate. The corresponding unloaded, dominant mode resonant frequencies were 3.434 GHz, 3.601 GHz, 3.727 GHz and 3.809 GHz for the 0.008 in., 0.004 in., 0.002 in. and 0.001 in. spacers, respectively, with the LaAlCb substrate. Without the LaAlC>3 substrate, the dominant mode frequencies were 4.194 GHz, 4.064 GHz, 3.995 GHz, and 3.956 GHz for the 0.008 in., 0.004 in., 0.002 in. and 0.001 in. spacers, respectively. Table 3-3 shows that the value o f s/pQeffective is greater than the actual value o f Rssc The value o f s/pQ approaches the actual value o f Rssc as the spacer thickness is decreased. At a spacer thickness o f0.008 in., s/PQ is 7.9 times larger than Rssc actual. At a spacer thickness o f 0.001 in., s/pQ is 15% larger than Rssc actual. Note that Q e ffe c tiv e is principally determined by super-conductor losses only when the resonator spacing s is equal to 0.001 in. and when a LaAlCb substrate is present since the fields are concentrated in the parallel plate resonator. When s is equal to 0.002 in. or larger, the losses in the walls o f the enclosure are comparable to or greater than the losses in the resonator, in which case the Qeffective is not simply related to the losses in the super conductor. 47 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 3-4 shows that the difference between s/fiQ and Rssc actual is smaller without the LaAIQj substrate. When compared to Table 3-3, the value o f s/pQ in Table 3-4 approaches the actual value o f R^c more rapidly as the spacer thickness is decreased. At a spacer thickness o f 0.008 in., s/pQ is 2.6 times larger than Rssc actual. The value o f s/pQ is only 43% and 8.4 % larger than RssCactual at spacer thicknesses o f 0.004 and 0.002 in., respectively, and is 1% smaller than Rssc at a spacer thickness o f 0.001 in. This indicates that without the LaAlCh substrate, the fields are more concentrated in the parallel plate resonator and the effective Q is determined principally by the losses in the super-conductor over a wider range o f spacer thickness. Only when s is equal to about 0.008 in. or larger are the losses in the walls o f the enclosure without the LaA 1 0 3 substrate greater than the losses in the resonator. According to Taber’s resonator model given by (3.6), the losses external to parallel plate resonator, i.e., the radiative loss, are proportional to the spacer thickness s. Under such an assumption, the reciprocal o f Q e n c io s u rc will be proportional to s. This assumption is examined using the FDTD simulation results in Figures 3-3 and 3-4. The reciprocal o f Q e n c io s u rc for the dominant parallel plate mode versus the spacer thickness s is plotted in Figure 3-3 including a LaAlOj substrate, and in Figure 3-4 without a LaAlCb substrate. For reference purposes, the dotted curve in Figures 3-3 and 3-4 represents a linear relationship between l / Q e n ciosurc and s, with an arbitrary constant o f proportionality a . The curve in Figure 3-3 indicates that the assumption of proportionality between 1/Q e n c io s u rc and the spacer thickness is not valid when the LaAlCb substrate is present. 48 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Without the LaAlCb substrate, Figure 3-4 shows that approaches the a s curve. The relationship between l/Q e n c io s u re I/Q e n d o s u re asymptotically and s is more complicated when the LaAlC>3 substrate is included in the analysis. It is anticipated that a different relationship between l / Q end o s u rc and s will hold for each parallel plate mode. The reciprocal of Qendosure for the TM02 parallel plate mode versus the spacer thickness s is plotted in Figure 3-5 including a LaA103 substrate, and in Figure 3-6 without a LaAlQ? substrate. In the case o f the TM02 mode, the value o f 1 / Q e n d o s u r e asymptotically approaches the a s curve only when the LaAlQj substrate is not present. Although the FDTD results indicate that the Taber model is accurate only in the absence o f a high dielectric substrate, Tables 3-3 and 3-4 indicate that effects o f the enclosure are at least minimized when the spacer thickness is very thin, say about 0.001 in. When the spacer is thin, most o f the loss will be due to the super-conductive thin film. The FDTD results are compared to measurements conducted by Northrup Grumman STC57 on a 0.5 in. x 0.25 in. parallel plate resonator with a 0.001 in. sapphire spacer in a 1.0 in x 0.75 in. x 0.967 in. high enclosure in Table 3-5. Data is shown for the TMi 1, TM21, and TM30 modes o f the parallel plate resonator. For all three modes, the calculated data indicates that the Qeffective is dominated by the Q due to the super conductor. In such a case, J/^Qeffective is expected to give a good estimate o f the superconductive surface resistance Rssc- 49 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 7.00E-04 6.00E-04 5.00E-04 S 4.00E-04 2m J 2 3.00E-04 2.00E-04 1.00E-04 O.OOE+OO 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 Spacer thickness (in.) Figure 3-3. Reciprocal o f Qenciosurc vs. s for dominant mode with a LaAlCb substrate. 8.00E-05 -j--------------------------------- :------------------ ---------- ---------- -------7.00E-05------------------------------------------ ---------- ---------- — - — ^ aa ss * 6.00E-05------------------------------------SI ' 1 4.00E-05---------- ;t--------------i-S £ “ ;----- i-------- 5.00E-05-----------i---------- ---------- ;--------------- 1 ■ i! : ------ ; Jr- j ^ 1--------------------vr j i !! i i! ' ! — |---------- | !i i----- 1---------1 : I 3.00E-05 ----------- 1---------- --------— -------------------- ;---------- !-----------i----------2.00E-05------------------1.00E-05---------- — ----------------------0.00E+00 L - 0 ; ; j ; 0.002 0.003 0.004 ■ ; | ! ;-----i---------- ;---------- — ------- :----- i---------- 1ii-----------ii----------- ---------- :-----------!---------- i---------- i---------- i1-----------iI----------- 0.001 0.005 0.006 0.007 0.008 Spacer thickness (in.) Figure 3-4. Reciprocal of Qenciosurc vs. s for dominant mode without a LaAlCb substrate. SO Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6.00E-04 5.00E-04 4.00E-04 »to 3 | 3.00E-04 g 2.00E-04 1.00E-04 O.OOE+OO 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 Spacer thickness (in.) Figure 3-5. Reciprocal o f Qenciosurc vs. s for TM02 mode with a LaAlO? substrate. 3.50E-05 ---------------- :-----:---------------------- ;---------- ----------- i-----------1 i---------- ;---------"" i----- 1 3.00E-05------------i---------- ---------------- i !1 \ \ 2.50E -05--------------------------------- ---------- 1---------- !------------------ S S. ' " = a ®t 2.00E-05-------------------- - J g i ' ' 1.50E-05--------------------------------------' ■" 1.00E-05 - j------ Jr ■ ^ ! i i ------- I j-j----------- » ! i - - - - - - - i- - - - - - - :- - - - - - - - - - - - - - - - - - - - - - - - j;- - - - - - - - 1!- - - - - - - !1- - - - - - - ! 1 ; 1 1 —--1-----------1!---------- !!---------- 1!---------i i ! ;’ !! I! ; 5.00E-06---------------------------i s- - - - - - - - - - - - i- - - - - - - - - - - - 1 - - - - - - - - - - - - i- - - - - - - - - - - - - !- - - - - - - - - - - - - j- - - - - - - - - - - - 1 ------------ 0 .0 0 E + 0 0 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 Spacer Thickness (in.) Figure 3-6. Reciprocal o f Q e n c io s u rc vs. s for TM02 mode without a LaAlC>3 substrate. 51 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 3-5. Comparison o f calculated and measured results for a 0.5 in. x 0.25 in. resonator. CALCULATED (FDTD) RESULTS1: Mode (m,n) (m is larger dimension) (3,0) (2,1) (1,1) Resonant frequency in GHz 8.930 11.829 11.282 Q due to losses in walls of enclosure: 20200 23518 41306 In Side Walls: 2084524 7566361 16578722 In Top and Bottom Walls: 20254 23786 41409 Q due super-conductor losses: 1346 1736 1368 Q due to dielectric losses: 5553468 5333417 5727609 1303 1616 1281 Q e f fe c tiv e 9.26E-04 5.54E-04 8.68E-04 R s = s / P Q e f f e c t i v c (ohms) MEASURED RESULTS2: Resonant frequency in GHz 8.737 11.450 10.950 1460 1962 2244 Q m e a s u re d 6.00E-04 R s calculated from Q m e a su re d (ohms): 5.60E-04 5.12E-04 R s theoretical (ohms) 5.90E-04 9.42E-04 1.04E-03 'Calculated results assume a 0.5"x0.25" resonator with a 0.001" spacer. 2Measured resonator is 0.508"x0.249" with a 0.00084" spacer. The measured resonant frequencies are from 2.2% to 3.2% lower than the calculated resonant frequencies for the three modes shown in Table 3-5. Some o f the discrepancy can be attributed to the slightly larger resonator (in length) and part due to the larger (in height) enclosure used for the measurements. The Q m e a s u re d for the TM u mode o f 1460 is slightly smaller than the Q c a ic u ia te d o f 1616, while Q m e a s u re d for the TM 21 and TM30 modes of 1962 and 2244, respectively, were slightly larger than the calculated values o f 1303 and 1281, respectively. The Q c a ic u ia te d leads to a value o f j / ^ Q c a i c u i a t e d which is slightly smaller than Rssc- The measured Q leads to a value o f j/^Qmeasured that is smaller than s //? Q c a ic u ia te d - The value o f Rssc estimated by jZ /J Q c a ic u ia te d decreases with frequency, which is contrast to the expected or theoretical frequency dependence which is a function o f the square of the frequency.' 52 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. In summary, calculations o f conductor losses with FDTD have shown that the LaAlC>3 substrate, with a dielectric constant o f about 22, can significantly affect the percentage of field energy outside the parallel plate resonator. This leads to significant field dissipation by the normal conductors outside o f the parallel plate resonator, in which case the assumption that losses within the parallel plate resonator dominate the effective Q is invalid. In such a case, the surface resistance o f the superconductive film will not be accurately given by s/f3Q. The minimization o f losses outside o f the resonator, perhaps by the addition o f a superconductive film to both sides o f a LaAlCb substrate, should improve the experimental accuracy. Experimental data was shown to support these conclusions. Analysis of microstrip dual patch resonator A calculation o f the resonant frequencies for a pair o f microstrip patches in a PEC cavity was carried out using the FDTD method. A analysis o f the resonant frequencies of a pair o f microstrip patches using the Method o f Moments was described by Grounds.26 A diagram o f the layout for a pair o f microstrip patches is shown in Figure 3-7. Both patches were 0.331 in. x 0.331 in. on an alumina substrate o f thickness 0.025 in. with er=9.8. The pair o f patches were separated by a distance 2d and were centered in a PEC cavity with dimensions 0.945 in. wide by (1.276+2d) in. long by 0.315 in. high. 53 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Planes o f symmetry T A 1 'r « ------------------------------ B -------------------------------------- ►! Figure 3-7. Layout of microstrip patches in a PEC cavity o f dimensions A=0.945 in., B=(l-276+2d) in., and Yp=Zp=0.311 in. The FDTD calculations were conducted using a one-quarter geometry with two planes o f symmetry. Either PEC or PMC boundary conditions were applied at y=0 or z=kmaxAz. With PEC conditions applied at y=0, the currents for the lowest order resonant mode flow predominantly in the y-direction, while for PMC conditions at y=0, the currents for the dominant resonant mode flow in the z-direction. Using the one-quarter geometry with PEC or PMC boundary conditons applied at the planes o f symmetry, a single distinct dominant mode was observed in the Fourier transform o f the time domain fields. At the beginning of the time domain simulation, a Gaussaian source was applied to the Ex-field at one of the comers o f the patch, as follows: 54 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Typically, the pulse width parameter ni was chosen such that the width o f the pulse was about equal to the length (or width) of the patch. Usually, no was set equal to 4 ni. During the simulation, the E- and H-fields were sampled around the perimeter o f the patch. Since the symmetry o f the problem determined the direction o f the dominant current flow, the sample points were chosen near the location o f peak H-fields. The resonant frequencies were calculated for patch spacings from 0.0625 in. to 0.6125 in. For a maximum patch separation o f 0.6125 in., the FDTD simulation used grids of 50x74x151 cells or 101x151x302 cells. The cell dimensions were chosen by fitting an integral number o f cells across the half-width and length o f the patch. For the smaller grid, the cell dimensions were Ax=0.00625 in., Ay= 0.00636 in. and Az=0.00624 in., while for the larger grid size the cell dimensions were Ax=0.003125 in., Ay= 0.00312 in. and Az=0.00312 in. The total number o f cells in the z-direction was varied depending on the separation between the patch and the symmetry plane at z=kmaxAz. A convergence test o f the FDTD simulation was run using the maximum patch separation, which is essentially the same case as a single patch centered in a 0.945 in. by 0.945 in. PEC cavity. The size o f the Yee cell was decreased from 0.0125 in. to 0.003125 in. and the number of cells increased to maintain the same geometrical arrangement. The results are listed in Table 3-6, along with the cpu time requied on a Convex C-3840 or a Cray YMPel supercomputer. The change in resonant versus the cell dimensions is plotted in Figure 3-8. Note that the resonant frequency increases 55 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. linearly with decreasing cell size, down to a cell dimension o f 0.00357 in. The difference in resonant frequency o f the dominant mode was 3 1 MHz (0.56%) using grids o f 50x74x151 and 101x151x302 cells. Using the smallest (25x37x74 cells) grid, the resonant frequency was 1.8% lower than that obained using the largest (101x151x302 cells) grid. The smallest grid only required 737 cpu seconds as compared to 69,695 cpu seconds for the largest grid. The resonant frequency of the dominant mode as a function o f the patch spacing 2d for both PMC and PEC boundary conditions at z=kmaxAz is plotted in Figure 3-9 for a 50x74x151 cell grid, and in Figure 3-10 for a 101x151x302 cell grid. Two resonances are shown for each boundary condition at ( w corresponding to the cases where the currents are y-directed or z-directed. The mode with currents that flow parallel to the wall at kmax (y-directed) are more influenced by the wall than the mode with z-directed currents. The resonant frequencies from the method of moments analysis13 is shown in Figure 3-11. For ease of comparison to the previously published results, the patch separation is plotted in centimeters in Figures 3-9 to 3-11. The FDTD results for the resonant frequencies in Figures 3-9 and 3-10 show a very similar dependence on the patch spacing, except that the results for the smaller grid is shifted down in frequency by about 35 MHz (about 0.6%). The difference in resonant frequency can be attributed to the larger numerical dispersion caused by a larger cell size. The FDTD results for the smaller cell size (101x151x302) grid agree agree to within a few megahertz of the method o f moments results, although the agreement was closer for the case where the currents were z-directed. A few experimental points are 56 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. shown in Figure 3-11 for the case where the currents are z-directed that correlate very well with the calculations. It should be noted that a dielectric constant o f 10.11 instead o f 9.8 was used in the Method o f Moments calculations to obtain good agreement between the calculations and measurements. A value o f 9.8 was used for all o f the FDTD calculations. 57 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 3-6. Calculated resonant frequency for a number o f different FDTD grids. No. o f Cells in Grid UmaxJ maxd^max) 25x37x74 50x74x151 63x94x188 76x114x225 88x131x265 101x151x302 Cell Dimension (in.) 0.0125 0.00625 0.005 0.00417 0.00357 0.003125 Resonant CPU Frequency Time (GHz) (secs) 5.449 737 13176 5.521 5.534 19417 5.543 35505 5.552 43577 5.552 69695 5.56 N s £ > » w s 5.52 5.5 Eh S eg S aVI v as 5.48 5.46 5.44 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Cell size (in.) Figure 3-8. Resonant frequency versus cell size for a microstrip patch resonator from the FDTD results. 58 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. - e wall frequency 1 — tn wall frequency t e wall frequency 2 : m wall frequency 2 1 5.65 5.6 S S * 5.55 S 5.5 es e 5.45 a 5.4 5.35 0.5 I 1.5 Patch Separation (cm) Figure 3-9. Resonant frequency versus patch separation for a 50x74x151 cell FDTD grid. — e wall frequency 1 — e wall frequency 2 - * —m wall frequency 1 —« - m wall frequency 2 5.65 W *5.55 S s *■ 55 2 3 tt. 1 5.45 e 5.35 0 0.5 1 1.5 2 Patch Separation (cm) Figure 3-10. Resonant frequency versus patch separation for a 101x151x302 cell FDTD grid. 59 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -o —e wall frequency 1 m wall frequency 1 - o - e wall frequency 2 -Or- m wall frequency 2 -H - measured frequency 1 -4— measured frequency 2 5.65 5.6 5.55 5.5 S’ 5.45 Uh 5.4 5.35 0 0.5 1 1.5 2 Patch Separation 2d (cm) Figure 3-11. Calculated resonant frequency for dual patch resonator from reference 26 (Reprinted with permission of author). Calculation of the resonant frequency of dielectric resonators and comparison with measurements The modeling of cylindrical dielectric structures inside o f rectan g u lar waveguide is a difficult numerical problem which was first analyzed by Liang27 using mode matching techniques. Here a identical structure will be analyzed using the FDTD method for comparison to the mode matching results. A diagram o f the resonator is shown in Figure 3-12. The dielectric portion o f the resonator consists o f a support 60 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. cylinder o f radius n and height hi and a larger dieletric cylinder o f radius iz and height h2 . The dielectric cylinders are centered at the bottom o f a section o f rectangular waveguide with width o f 2a and a height o f b. The length o f waveguide is equal to L. SIDE VIEW TOP VIEW 2^1 Figure 3-12. Diagram o f a cylindrical dielectric resonator in a rectangular waveguide. Using the mode matching method, Liang27 analyzed the case where the dielectric constants o f the lower and upper dielectrics were equal to 1.0 and 38.0, respectively, and the waveguide dimensions were 2a=b=1.0 in. and L=0.92 in. For comparison to Liang’s results, the FDTD simulation was run using nominal grids o f 20x20x18 cells, 40x40x36 cells, and 60x60x54 cells for three different cylindrical resonator sizes. The results are shown in Tables 3-7 to 3-9. In Tables 3-7 to 3-9, the dimensions o f the waveguide used in the FDTD simulation differed slightly from the actual waveguide due to discretization error. Generally, an odd integral number o f cells was chosen to fit the diameter or height o f 61 of the copyright owner. Further reproduction prohibited without permission. the dielectric exactly, which usually resulted in a slightly larger o r smaller waveguide than that used by Liang.27 Only cells within the radius or the height represented by the next smallest integer were given a dielectic constant equal to 38. This was done because the effective diameter or height o f the dielectric is slightly larger than the diameter or height represented by the Yee cells. By comparing FDTD results to known results for a cylindrical dielectric in a cylindrical waveguide, it was found that the effective diameter o f the dielectric extends one-half cell beyond the edge o f the Yee cells on a boundary between different materials. The same is true for the height of the dielectric. So if seven cells fit across the diameter o f the dielectic cylinder, only cells within a radius o f six were given a dielectric constant o f 38. Zhang19 has shown that an alternative to assuming that the dielectric extends an extra one-half cell is to set the dielectric constant of a Yee cell on a boundary between different materials to the average o f the dielectric constants. Results are shown in Table 3-7 for a FDTD grid with a nominal cell dimension o f Ax«Ay~Az»0.05 in. For comparison, the calculated results using the mode matching method and the measured results are also shown. The resonant frequency calculated using FDTD agreed to within 0.7% o f the measured resonance for the three dielectric cylinders. The mode matching results were very close to the FDTD calculations. The FDTD simulation required about 200 CPU seconds for 32,768 time steps on a Convex supercomputer. The calculated resonant frequencies using a 40x40x36 FDTD grid in shown in Table 3-8 for same resonator dimensions as Table 3-7. In this case, the nominal cell 62 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. dimensions were Ax«Ay«Azw0.025 in. Surprisingly, no improvement in accuracy over the smaller grid was observed. In fact, the results for the larger grid in Table 3-8 were slightly worse than those shown in Table 3-7. It would seem that by chance, the discretization errors for the smaller grid (i.e.’larger cell size) affected the resonant frequencies in such a way as to improve the correlation with the measured results. The FDTD results shown in Table 3-8 required about 4,900 CPU seconds for 131,072 time steps on a Convex supercomputer. In Table 3-9, the two largest dielectrics are analyzed using a 60x60x54 cell grid. In this case, a value o f 0.0167 in. was used for the three dimensions o f the Yee cells. The FDTD results in Table 3-9 agree almost exactly with the measured results. Ten hours o f CPU time were required for 262,144 time steps on a Convex C-3840 supercomputer with a 12 nsec clock speed. Table 3-7. Resonant Frequencies for a 1.00" x 1.00" x 0.92" Rectangular Cavity with a Cylindrical Dielectric (sr =38) placed 0.275" above the bottom o f the waveguide with a Cartesian Grid of 20x20x18 cells. Dielectric Diameter (in.) 0.654 0.689 0.757 Dielectric Thickness (in.) 0.218 0.230 0.253 No. of Cells FDTD (MHz) 4,414 4,161 3,792 20x18x18 2 2 x2 2 x2 0 20x20x18 Mode Matching14 (MHz) 4,388 4,161 3,772 63 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Measured14 Resonance (MHz) 4,382 4,153 3,777 Table 3-8. Resonant Frequencies fora 1.00" x 1.00" x 0.92" Rectangular Cavity with a with a Cylindrical Dielectric (er =38) placed 0.275" above the bottom o f the waveguide with a Cartesian Grid o f 40x40x36 cells. Dielectric Diameter (in.) 0.654 0.689 0.757 Dielectric Thickness (in.) 0.218 0.230 0.253 No. of Cells 38x37x36 40x39x36 40x40x38 FDTD (MHz) 4,423 4,142 3,748 Mode Matching14 (MHz) 4,388 4,161 3,772 Measured14 Resonance (MHz) 4,382 4,153 3,777 Table 3-9. Resonant Frequencies for a 1.00" x 1.00" x 0.92" Rectangular Cavity with a with a Cylindrical Dielectric (er =38) placed 0.275" above the bottom o f the waveguide with a Cartesian Grid o f 60x60x54 cells. Dielectric Diameter (in.) 0.689 0.757 Dielectric Thickness (in.) 0.230 0.253 No. of Cells 60x60x56 60x60x56 FDTD (MHz) 4,152 3,779 Mode Matching14 (MHz) 4,161 3,772 64 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Measured14 Resonance (MHz) 4,153 3,777 Chapter 4 Application of the FDTD Method to Analysis of Stripline Circuits The design of stripline circuits has recently received renewed interest because o f the emergence o f improved construction techniques. The synthesis o f filters, matching networks and couplers often require equivalent circuit models for stripline circuits and stripline discontinuities. In addition, analysis and synthesis o f stripline circuits using equivalent circuit models are generally much faster and more efficient than full wave analysis methods. In this chapter, the characteristic impedance and propagation constant o f strip transmission line, and S-parameters for open, step and tee-junction stripline discontinuities, are calculated using the Finite Difference Time Domain (FDTD) method. It is shown that the FDTD results for characteristic impedance and propagation constant of strip transmission line agree with existing models. However, significant discrepancies exist between the FDTD results for the S-parameters o f stripline discontinuities and the previously published S-parameter models. New equivalent circuit models and new values for the S-parameters are developed for the stripline discontinuities. The new circuit models are then used in the synthesis of stripline low pass filters circuits, and the measured S-parameters are compared to values calculated using the new equivalent circuit models. Propagation constant and impedance for symmetric stripline The propagation o f electromagnetic energy in symmetric stripline is essentially transverse electromagnetic mode (TEM). Therefore, symmetric stripline is nondispersive with a propagation constant equal to - - -7 = , where Xq is the free space ^oV8 65 of the copyright owner. Further reproduction prohibited without permission. wavelength and e is the dielectric constant o f the substrate. The TEM condition holds up to the onset of the first higher order propagation mode, which for wide stripline occurs at an approximate frequency fnoM in GHz given by Vendelin58: 591 1 (4.1) 1h o m ~ b 4 where W is the width and b is the ground plane spacing in inches, and sr is the dielectric constant. For stripline with W=0.036 in., b=0.060 in., and Er=3.0, the frequency of the 1st higher order mode predicted by (4.1) is 41.4 GHz. The impedance o f symmetric stripline having infinitely thin metal for the center conductor was investigated by Cohn59 using conformal mapping techniques. Cohn’s formula for the impedance Z in ohms o f thin stripline is: . 30* K'(k) , Z = —r = : , where A K(k) -i- i K(k) K’(k) n I 1-Vk (for 0 < k < 0.7) (4.2) as 1 In 1+ VTj 1-Vky (for 0.7 < k < l) where k = tan h ^ -^ -j and k ' = V l - k 2. For stripline with W=0.036 in., b=0.060 in., and er=3.0, (4.2) gives an impedance of 52.35 ohms. When the thickness of the center conductor is included, Gupta29 and Wheeler60 provide the following formula for the impedance of thick stripline as a function o f the 66 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ground plane spacing b, the substrate dielectric constant Er, and the width W and the thickness t o f the center conductor: 7 = - 30Z 7 = 1 1. 1 1+ V8r W' b -t 4 b -t 8 b -t 71 W* ,7C W’ f8b-i)2 £ +\U W’ J , where W AW -+ b -t b -t’ \m AW b —t 7t(l —x) x m = 2 1+ 3 I —x 2 0.0796x l-iln -1 , and x= (4.3) t The formula in (4.3) has a stated accuracy o f 0.5% for W7(b-t) <10. For stripline with W=0.036 in., b=0.060 in., er=3.0, and t=0.00125, (4.3) gives an impedance of 49.68 ohms. The application of the FDTD method to the analysis o f microstrip tra n sm issio n was described by Zhang. 17 The determination of characteristic impedance using FDTD can be applied to any transmission line structure. Here, the method will be applied to the analysis o f symmetric strip transmission line (stripline). The FDTD analysis begins with a three-dimensional model o f the symmetric strip transmission line on a Cartesian grid, as shown in Figure 4-1. Only half of the center conductor is represented in Figure 4-1 since the stripline is assumed symmetric about the y=0 plane. The complete grid size is (imaxjmaxJw) cells in the (x,y,z) coordinate directions, respectively. The stripline center conductor is located at x=imctAx= imaxAx/2, between y=0 and y=jmetAy, and extends from z=0 to z=kmaxAz. 67 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Rather then setting the conductivity cy in (2.3), the Ey- and Ez-fields are set to very small values on the surface of the metal by setting the values of the dielectric constants in the x- and y- directions, i.e., ex and ey, to some very large value. This is often simpler that setting the fields to zero and excluding the metal cells from the update equations. In this case, the value of the metal conductivity in (2.3) is not used. The E (ijjc) for the non-metallic cells are set to the dielectric constant of the stripline. PEC conditions are applied at x=0 and x=imaxAx to represent the ground planes, and PMC conditions are applied at y=0 to represent the plane of symmetry. Super-absorbing boundary conditions are applied at z=0 and z=kmaxAz, and at y= jmaxAy, to represent open boundaries. At the beginning of the simulation, a time-dependent electric pulse is introduced near one end o f the center conductor and allowed to propagate along the length o f the transmission line. The source was introduced by adding a extra term containing a Gaussian-time dependence to the update equations for the Ex-fields, as follows: (4.4) where n is the time step number, and no and m are the determined the zero offset and the width o f the input pulse. The constant no was usually chosen equal to 4ni to keep initial transients low. 68 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TOP VIEW _L KnvA*- Az T ------------- ► y iract^y w *y FRONT VIEW SIDE VIEW t hnet^X 0 0 i I 'm«Ax i. T Metal Ax ---------------------► jmet^y Metal y jmixAy 0 ---------------------------------- ► Z o •w ^z Figure 4-1. Layout o f Cartesian grid for symmetric stripline. The spatial distribution o f the source is not critical, and (4.4) can be applied between i=0 and i=imct, and between j=0 andj=jmet. For a symmetrical excitation, a equation similar to (4.4) but with a minus sign in front o f the exponential term is applied between i=im« and i=imax, and between j=0 ahd j=jmet. The source plane can be located a few cells away from the boundary at k=0. The source described above initially has a rectangular shape in the x-y plane and therefore contains several higher order spatial harmonics instead of a single, propagating TEM mode. However, all of the higher order modes except the fundamental TEM will not propagate very far from the source plane at frequencies well below the cut-off frequency of the first higher mode. Thus, only a single propagating (TEM) mode will be observed a few cells further down the stripline in the positive z-direction. The source equation (4.4) also initiates a wave that travels 69 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. along the transmission line in the negative z-direction, but it is absorbed by the ABC at z=0 . The voltage between the main transmission line and the ground plane are recorded at several locations along the transmission line as a function o f time as the pulse propagates along the stripline. A t the nth time step, the voltage V at z=kAz is calculated as the line integral of the Ex-field between the ground plane and the center conductor as follows: -i (4.5) The voltage given by (4.5) is independent o f the particular j-location as long as the pulse is well developed. This assumes that the samples are taken far enough away from the source such that any higher order spatial modes present initially have attenuated to zero. The voltage between the center conductor o f a stripline and the ground plane is plotted versus time in Figure 4-2 for three different locations along the z-axis. In Figure 4-2, the stripline had a center conductor frill width o f0.036 in. with a ground plane spacing o f 0.060 in. and a dielectric constant o f 3.0. A 12x80x100 cell grid was used with cell dimensions of 0.05 in., 0.03 in., and 0.03 in. in the x,y,z directions, respectively. With a time step of 0.1271 ps and a value of 110 for ni in (4.1), the bandwidth of the input pulse was approximately 50 GHz. The input pulse was applied at k=5 using (4.4), and the three voltages were calculated using (4.5) at k=60,80, and 70 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 100. Note that the pulses maintain the same shape at the three different locations along the stripline, as expected since the propagation mode is TEM. 12 V(nAtJc=60) V(nAtJc=l20) V(nAt,k=l80) 10 8 o 6 > 4 2 0 200 400 600 800 1000 1200 1400 Time Step (n) Figure 4-2. FDTD voltage vs. time between center conductor of stripline and ground plane at k= 6 0 ,120 and 180. Following the procedure described by Zhang,9 the propagation constant o f the transmission line y(co) is determined from Fourier transform o f the time history o f the voltage at two locations ki and k2 along the direction o f propagation as follows: y ( a ) = a ( a ) + /P ( g) ) = 1 ( k 2 -k ,)A z 1 (k 2 -k ,)A z In VCco.^Az) V (a ,k 2 Az) In F{V (t,ktAz)} F{V(t,k2Az)} where F{} denotes the Fourier transform. 71 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.6) In general, the propagation constant found using (4.6) is complex. However, when the transmission line is lossless, the real part of y(co) in (4.6) is very small or negligible up to the frequency o f the first higher order mode. The imaginary part ofy(co) is plotted versus frequency in Figure 4-3 up to 50 GHz using the stripline voltages from Figure 4-2 at k|=60 and k2=180 for V(t4qAz) and V(tJc2Az) in (4.6). The imaginary part o f the propagation constant is proportional to frequency, as expected for a dispersionless transmission line. The real part of the propagation constant was negligible up to 50 GHz. Above 50 GHz, the real part o f the propagation constant increased with frequency, suggesting the onset of a higher order propagation mode. 50 45 40 35 30 25 20 15 10 5 0 0 5 10 15 20 25 30 35 40 45 50 Frequency (GHz) Figure 4-3. Imaginary part o f propagation constant vs. frequency for 0.036 in. wide stripline. The impedance of the stripline can be defined in terms of V(to)/I(<o), where V(ffl)and I(co) are the Fourier transforms of the voltage V(t) and current I(t), 72 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. respectively. The voltage Vk(t) at a position z=kAz is calculated using by the line integral o f the E-field between the ground plane and the center conductor o f the stripline, as shown in (4.5). The current at a position z=kAz is calculated using the loop integral o f the H-field around the center conductor o f the stripline, as follows: l(nAt, z = kAz) = 2 A y£ j+ k+D - 1, j + 1 k + D (4.7) -2 A x H "[i,jmet + ^ , k + ^ ) . The factors o f 2 on the right hand side o f (4.7) account for the symmetry about the plane y=0. The current versus time, calculated from the FDTD H-fields using (4.7) at k=120, is shown in Figure 4-4 for a 0.036 in. wide stripline. B. E at 0.15 0.1 0.05 0 200 400 600 800 1000 1200 1400 Time Step (n) Figure 4-4. FDTD current versus time at k=120. Using the Fourier transforms of Figures 4-2 and 4-4, the impedance Z(a) is given by: 73 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Zf m) - Vt(<D) - Fi V < n A l - 2 = kAz>> It (m) f4 n F{l(nAt, z = kAz)} The impedance given by (4.8) will ordinarily have a small imaginary component. Fang61 attributes this to differences in position and time at which the current I(t) and the voltage V(t) are evaluated using fields that are offset by one-half o f a cell in the zdirection. In addition, the voltage is evaluated one-half o f a time step later than the current since the E-fields are updated one-half of a time step later than the H-fields. Following the suggestions by Fang, the spatial offset in the z-direction between the voltage and current is compensated by using the average o f the Eg-fields at positions z=kAz and z=(k+l)Az to calculate the voltage, so that the voltage Vk+i/2 spatially coincides with the current which is derived from H-fields at position z=(k+l/2)Az. The H-fields are also translated by one-half cell in the both the positive and negative zdirections to coincide with the voltage. In addition, the voltage is offset by one-half time step by multiplying V(<o) by exp(/o)At/2). The corrected impedance is given by: (Vk(ffl) + Vw (ffl))e'^M ( ' " 2 I t (c»(e « " > “ 2 + ( ’ where y(<o) is the propagation constant given by (4.6). The impedance Z(o) for a 0.036 in. wide stripline with a 0.060 in ground spacing and sr=3.0 is shown in Figure 4-5. The impedance in Figure 4-5 is calculated using the time domain voltage and current shown in Figures 4-2 and 4-4, respectively, in (4.6). The FDTD results in Figure 4-5 indicate a constant impedance at 49.95 ohms up to a frequency of 50 GHz. The imaginary part o f the impedance calculated using (4.5) 74 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. begins to increase at approximately SO GHz, which suggest the onset of higher order modes o f propagation. ^ 51 ----------------------------------------------------------------- ,---------------------- 5 0 .8 ------------------------------------------------------------------- i---------------------- 5 0 .6 ---------------------- !--------------------- to 5 0 .4 - - - - - - - - - - - - - - - - - - - - - - - - - i - - - - - - - - - - - - - - - - - - - - - - - - i- - - - - - - - - - - - - - - - - - - - - - - s $ A) 1 I - 50 2 ----------------------- _ 50 i- - .----------------------1--------------------. - , ■ ■ ' l- I j _ .. . . 4 9 .8 ---------------------- :--------------------- :----------------------1--------------------- j---------------------- ' 4 9 .6 --------------------- 1---------------------- 4 9 .4 --------------------------------------------------------- :------------------------------ i---------------------4 9 2 ---------------------- ---------------------------------------------------------------------------------------49 ----------------------------------------------------------------- :--------------------- ---------------------0 10 20 30 40 50 Frequency (GHz) Figure 4-5. FDTD Impedance vs. frequency for a 0.036 in. wide stripline with 0.060 in. ground plane spacing and er=3.0. The impedance for several stripline widths between 0.006 in. and 0.036 in. calculated using FDTD are shown in Figure 4-6 for b=0.060 in. and er=3.0. The calculated results from formulas (4.2) for stripline with zero thickness metal and (4.3) for stripline with 1.25 mil (0.00125 in.) and 2.5 mil metal are also shown in Figure 4-6. As for the previous data, the FDTD calculations used a 12x80x100 cell grid with cell dimensions o f 0.05 in., 0.03 in., and 0.03 in. in the x,y,z directions, respectively. 75 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 120 - 0 mil metal — 1.25 mil metal 110 — 2.5 mil metal f 100 • FDTD results ! w s •o Si £ 0 0.01 0.02 0.03 0.04 0.05 Stripline Width W (in.) Figure 4-6. Comparison of stripline impedance vs. width from the FDTD results and from published formula for several thicknesses o f metal. The FDTD impedance data in Figure 4-6 agrees very closely to values obtained from published formula. The FDTD results agree within about 2% when compared to formula for the impedance using a metal thickness of 1.25 mils. Generally, the FDTD results do not agree as closely with the published formula for stripline with zero thickness metal, even though the metal is represented in the FDTD grid using a zero thickness approximation. Rather, the above results suggest that the FDTD results represent a metal thickness equal to about one-quarter o f a cell width. Good agreement has been obtained between the FDTD results and published formulae for impedance. The FDTD results also suggest that the mode of propagation is TEM at least up to a frequency of 50 GHz. The stated accuracy o f the published 76 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. formula has served to validate the FDTD results presented here. The data presented provide a cross-check between the FDTD results and the formula in (4.2) and (4.3). Determination of S-parameters for stripline discontinuities In this section, S-parameters for open, step and tee-junction stripline discontinuities are calculated using the Finite Difference Time Domain (FDTD). The Sparameters calculated using FDTD are compared to values obtained from the Method o f Moments and from values obtained using previously published equivalent circuits. Generally, S-parameters obtained using FDTD and the Method o f Moments were in close agreement. It is shown that significant discrepancies exist between the FDTD results and the earlier equivalent circuits. For the S-parameter calculations, a three-dimensional Cartesian grid is used to model a section o f transmission line and the stripline discontinuity. The FDTD grids for the discontinuities are similiar to the grids used for the striplines shown in Figure 4-1, except for the center conductor. Diagrams o f the center conductors for stripline open, step, and tee-junction discontinuties are shown in Figure 4-7. The reference planes for each discontinuity are marked with a “T.” The equivalent circuit for each discontinuity is also shown in Figure 4-7. 77 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. i 'oc T T w j (a) Openend T i - a n w, az Z2 T (b) Step junction (c) Tee-junction Figure 4-7. Stripline discontinuities and their equivalent circuits. 78 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Calculation of S-parameters and equivalent circuit for stripline openend discontinuity At the beginning of the simulation, a time-dependent electric pulse is introduced at one end o f the stripline and allowed to propagate to the openend discontinuity. The input pulse is the same as that used for the propagation constant and impedance calculations o f the previous section. The incident and reflected waves are sampled some distance away from the openend so that any evanescent mode has a chance to die o u t The incident and reflected voltage from an open stripline discontinuity is shown in Figure 4-8. The stripline was 0.020 in. wide between ground planes spaced 0.055 in. apart with a dielectric constant of 3.9. The FDTD simulation used a 22x120x200 grid with a cell dimension o f 0.0025 in. in all three dimenstions. An ABC was used at x=0, x=imaxAx, y= jmaxAy and z^k^axAz, and a PMC condition was applied at y=0. The stripline was 140 cells long and 4 cells wide (half-width). The openend o f the stripline was located at k=140, which was 60 cells from the ABC at kmax. The reflected voltage versus time in Figure 4-8 was recorded at k=100, which was 40 cells from the open end o f the stripline. The simulation required about 800 seconds on a Convex C-3840 supercomputer. 79 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12 Incident Reflected 10 8 BA 6 £ 4 0 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Time Step n Figure 4-8. Incident and reflected voltages for a stripline open discontinuity (W=0.020 in., b=0.055 in. er=3.0). Referring to Figure 4-8, the reflection coefficient for the stripline open discontinuity is seen to be equal to unity. The reflection coefficient versus frequency can be made more precise by taking the Fourier transform of the time series, and calculating the reflection coefficient, or Si i parameter, as follows from Zhang62: Vref((o,z = kAz)er(‘“>L " Vinc(<a ,z = k A z )e-''<“ ><-- (4.10) In (4.10), Vref(o)) and Vjnc(<a) are the Fourier-transformed reflected and incident voltages from position z=kAz, y(©) is the propagation constant, and L is the distance between the sampled voltages and stripline open. The exponential factors in (4.10) reference the Sn parameter to the plane of the openend discontinuity. In the case of Figure 4-8, L is equal to 40Az. 80 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. To separate the incident and refected voltages for use in (4-10), two simulations are required. The first simulation uses a length of stripline that includes the discontinuity, while the second simulation uses a straight length o f stripline without the discontinuity, that is terminated in an ABC. The same source is used and the time series o f the voltage is recorded at the same location along the stripline in both cases. The first simulation yields a waveform that includes both the incident wave and the wave reflected from the discontinuity, while the second simulation contains only the incident wave. The second simulation is subtracted point by point from the first simulation to remove the incident wave. The resulting reflected wave is shown in Figure 4-8 along with the incident wave. The phase of Sn(<o), calculated using (4.10), is shown in Figure 4-9 for three stripline openend discontinuities, with b=0.055 in., er=3.9, and widths o f 0.010 in., 0.020 in., and W=0.040 in. The magnitude o f Si i (not shown) was equal to one up to a frequency o f 50 GHz. Above 50 GHz, the magnitude of Si i became less than one, presumably due either to radiation from the openend or to the onset o f higher order modes o f propagation. 81 1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Altschuler and Oliner -10 FDTD 4 -20 eo TuJ c« fi -3 0 W = 0 .0 1 0 in . -4 0 W =0.020 in. -5 0 W =0.040 in. -7 0 -80 -9 0 0 10 20 30 40 50 60 Freq(GHz) Figure 4-9. Phase of S 11 versus frequency for a stripline open discontinuity for stripline widths of 0.010 in, 0.020 in. and 0.040 in. (5=0.055 and er=3.9) The S( i calculated from the equivalent circuit published by Altschuler and Oliner30 is also shown in Figure 4-9. The Su from Altschuler and Oliner’s work is based on the following equation for the apparent increase in length A1 due to the fringing capacitance at a stripline openend discontinuity: Al = —tan 1 P „ 8 S + 2W tan(pS) , where 48 + 2W (4.11) bln( 2 ) , „ 2* ^ = --------- , and (5 = —f — ■ TZ A.n The open circuit capacitance Coc o f the stripline openend discontinuity is related to Al by the following equation: C°c (4.12) cZ 82 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Using Al from (4.11), the Sn parameter at the reference plane‘T’ indicated in Figure 4-7 is given by:29 s 11 _ Zoc ~ z o _ ~ 7Z o cot(pAl) - Z0 1- /tan(PAl) Z ^ + Zq —/Z 0 cot (P Al) +- Z0 l + /tan (p A l)’ 1 j where ZoCis the impedance of the discontinuity. The phase o f the Si i given by (4.11) and (4.13) is plotted in Figure 4-9 along with the FDTD results. In comparison to the FDTD results, Altschuler and Oliner’s equivalent circuit consistently underestimates the phase o f Si i over the frequency range o f 0 to 50 GHz. The difference is equal to about 10 degrees at 40 GHz for the three cases shown in Figure 4-9. Bases on the FDTD simulations, it appears that (4.11) underestimates the apparent increase in length at the stripline openend, or equivalently, underestimates the openend fringing capacitance Coc• Similar results occurred for other stripline openend configurations with dielectric constants ranging from 3.0 to 22 and widths ranging from 0.002 in. to 0.040 in. Measurements will be shown later in this chapter to support the FDTD results. The S-parameters calculated using the FDTD results can now be used to refine the value o f the Al, or equivalently the fringing capacitance via (4.12), in the equivalent circuit by solving (4.13) for the admittance of the open discontinuity63: Y ^ m) = G ^ ) + jo X ^ a ) = Z 0 1+ Sn (©) . (4.14) The conductance Goc in (4.13) is negligible up to the onset of the first higher order mode, and will not be considered further here. The capacitance Coc(co) calculated from 83 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.13) and the value of S n from the FDTD calculations was used to define a new value for Al, which was then used to in a least squares curve fit to determine a new empirical formula for Al. The new formula was o f the form: Al = b ;r= r-(l + E .fV i7 + O f 2), (4.15) c+DhrJ where f is the frequency in GHz and A ,B ,C ,D JE ,an d G are unknown coefficients to be determined from a least squares curve fit to the FDTD data. The coefficients for (4.13) from the curve fit procedure are listed in Table 4-1. A comparison o f the new formula in (4.13) to the phase of SI 1 calculated with FDTD is shown in Figure 4-10. The new formula matches the FDTD results to within about 2 degrees of phase up to a frequency of 30 GHz. More accurate curve fitting can be achieved when the data is fit over a more narrow range of the dielectric constant, instead of from 3. to 22. as used here. Table 4-1. Values o f coefficients for (4.13) Coefficient A B C D E G Value 0.0268 0.1695 0.2717 0.7233 1.589x10'* 5.16x10'* 84 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -10 New Formula | -20 -30 -40 -50 -60 W=0.010 in. -70 f 7W=0.020 in. -80 W=0.040 in. -90 0 10 20 30 40 50 60 Frequency (GHz) Figure 4-10. Comparison o f the Phase of Sn calculated with the new formula to the FDTD results. Calculation of S-parameters and equivalent circuit for stripline step junction discontinuity The FDTD analysis o f the stripline step discontinuity follows a set o f procedures that are similiar to those for the open discontiuity. In the case o f the step discontinuity, both the reflection Sn and the transmission S21 parameters are required to characterize the junction. Therefore, both o f the reflected and transmitted wave are required from the FDTD analyis. The reflected wave and transmitted waves are recorded a few cells away from the junction so as not to include any evanescent modes. The incident and reflected voltages for a stripline step discontinuity are shown in Figure 4-11, and the transmitted wave in shown in Figure 4-12. The input line width Wj was 0.040 in., the output line width W2 was 0.020 in., and the ground plane spacing 85 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. b was 0.055 in. with a dielectric constant er o f 9.8. The FDTD simulation used a 22x100x240 grid with a cell dimension o f 0.0025 in. in all three dimenstions, with a time step At o f 0.1058 psec. An A B C was used at x= 0, x = imaxA x, y = jmaxAy and z= knuxAz, and a PMC condition was applied at y=0. The step junction was located at k=160, and the incident and reflected voltages were recorded at ki=60, while the tranmitted voltage was recorded at k2= 2 0 0 . The port 1 reflection coefficient, i.e. the Si t parameter, for the stripline step junction can be calculated using the Fourier transform o f the incident and reflected voltages versus time from Figure 4-1 1 in (4-10) with L=l00Az. The forward transmission coefficient, i.e. the S21 parameter, can be calculated using the Fourier transform o f the incident and transmitted voltages versus time from Figures 4-11 and 412, Vinc(oi) and Vmnfco) in the following formula20: Y(o>)L, (4 .1 6 ) where value o f L/=(k-ki)Az and L2 =(k 2-k) Az in (4.16) references the S21 parameter to the plane o f the step junction, which occurs at k=160. The square root factor nomalizes the S-parameter to the square root o f the respective impedance at each port. 86 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 25 Incident 20 15 ea 10 Reflected 5 0 0 400 800 1200 1600 2000 2800 2400 Time Step n Figure 4-11. Incident and reflected wave for a stripline step junction (Wi=0.040 in. W2= 0 . 0 2 0 in., b=0.055 in., er=9.8). 30 Transmitted 25 20 10 5 0 0 400 800 1200 1600 2000 2400 2800 Time Step n Figure 4-12. Transmitted wave for a stripline step junction (Wi=0.040 in. W2= 0 . 0 2 0 in., b=0.055 in., er=9.8). 87 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The S-parameters for the stripline step junctions can also be calculated using the equivalent circuit shown in Figure 4-7b and values for the reactance X and for the transmission line lengths U and I2 given by the following formulae from Altschuler and Oliner30: U Z *; _ _ i _ b ]n ( 2 ) M— I? — * 7t (4.17) 1 where Di is the equivalent width o f the stripline o f width Wj which is given by:29 D= b f o o + * [ | ~ ln(2 t / b ) I’ (f o r W / b s °-5) /h f W + — ln(2) + - [ l - I n ( 2 t / b ) l ( f o r W / b > 0 5 ) 7T (4.18) 7C where k = tanhf -^ 7-) and :/f.~ w as given in (42). V 2 bv K(k') The S-pareameters are now given by the following formula with X and h given by (417):29 (z 2 - z , f j x ) 11 „ _— 21 Z2 + Z , + j X ,—-— (4.19) Z2 - Z , + j X ’ where Z1 and Z2 are the impedances at the two ports. The magnitude and phase o f the Sn(co) and S2 i(o>) parameters for the stripline step junction with Wj=0.040 in., W2= 0 . 0 2 0 in., and b=0.055 in. are shown in Figures 4- 88 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13 to 4-16 for a substrate dielectric constant er=3.9; in Figures 4-17 to 4-20 for a substrate dielectric constant er=9.8;; and in Figures 4-21 to 4-24 for a substrate dielectric constant er=24. Three other sets of S-parameters are also in each o f the Figures 4-17 to 4-24, which were calculated using the Mode Matching method,4 a commercial computer-aided design package’called LIBRA64, and the equivalent circuits o f Altschuler and Oliner.30 The Mode Matching results are only shown up to a frequency of 20 GHz. A value of t=0.00125 in. was assumed for the thickness o f the metal in the LIBRA and the Alshuler and Oliner calculations. The four calculations for the magnitude o f the Si i and S21 parameters for the step junction were in close agreement. The largest discrepancy occurred for the magnitude o f the Si 1 parameter (er=3.9 or 9.8) in which case both the FDTD and the Altschuler and Oliner results were about 5% lower at low frequency than the S-parameter calculated using Mode Matching or LIBRA. The four calculations for the phase 0 FS21 were in fairly close agreement, with all results showing less than 6 degrees o f phase shift. The FDTD and LIBRA results ageed to within about 1 degree o f phase, but were a few degrees higher than both the Mode Matching and the Altschuler and Oliner results. The discrepancy, however, amounted to less than 2,3, or 4 degrees for er=3.9,9.8 or 24, respectively. 89 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ______________ 0.05 -— FDTD — Mode Matching — Libra ____________________________ i_____ t — Altschuler and Olinerl 10 20 30 40 SO Freq(GHz) Figure 4-13. Magnitude o f S u for stripline step junction (Wt=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., er=3.9). -20 OH « -30 -40 « -50 '60 — FDTD — Mode Matching jq j— Libra .80 1! 0 Altschuler and Oliner 10 30 20 40 50 Freq(GHz) Figure 4-14. Phase o f S 11 for stripliine step junction (Wi=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., er=3.9). 90 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ___________ ____________ i 1 ___________ 1 - i ~ ; i! 1 - i ________ . 0.95 < ff «■* e a» s sfifi 01 S 0.9 : ' ; ■ t i i r : 0.85 ; -*-FDTD i Mode Matching Libra ----- 0.8 A ltsc h u le r an d O lin er 1 8 12 16 20 Freq(GHz) Figure 4-15. Magnitude o f S21 Stripline step junction (Wi=0.040 in. W 2=0 . 0 2 0 in., b=0.055 in., er=3.9). off -•-FD T D -12 : Mode Matching I Libra t 0 10 20 30 Altschuler and Olinerj 40 50 Freq(GHz) Figure 4-16. Phase of S21 for stripline step junction (Wi=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., er=3.9). 91 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.25 FDTD -*-MM ! Libra d Altschuler and Oliner 0 5 10 15 20 25 30 Freq(GHz) Figure 4-17. Magnitude o f Su for stripline step junction (Wi=0.040 in. W2=0.020 in., b=0.055 in., sr=9.8). ’ -20 aV2 -40 V) -60 -80 FDTD MM Libra Altschuler and Oliner -100 -120 0 5 15 10 20 25 30 Freq(GHz) Figure 4-18. Phase of Su for stripline step junction (Wi=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., er=9.8). 92 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.95 «F e 4» •O 06 « S 0.85 FDTD MM Libra Altschuler and Oliner Freq(GHz) Figure 4-19. Magnitude o f Sn for stripline step junction. (Wj=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., sr=9.8). 0 -2 ■3 -4 -5 FDTD MM Libra Altschuler and Oliner -6 -7 0 5 15 10 20 25 30 Freq(GHz) Figure 4-20. Phase o f S21 for stripline step junction (Wi=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., er=9.8). 93 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.4 FDTD Mode Matching Libra Altschuler and Oliner 0J e 4» 1 '£ uat s °-2 0.1 0 0 4 8 12 16 20 Frcq(GHz) Figure 4-21. Magnitude o f Sn for stripline step junction (Wt=0.040 in. W2=0.020 in., b=0.055 in., er=24). -20 -40 -60 FDTD Mode Matching Libra Altschuler and Oliner •80 -100 0 4 8 12 16 20 Freq(GHz) Figure 4-22. Phase o f Si i for stripline step junction (W i=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., sr=24). 94 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. aFa W O <1 O C8H 2 -•-FDTD -•-M ode Matching -•-Libra Altschuler and Oliner 8 12 16 20 Freq(GHz) Figure 4-23. Magnitude of S21 for stripline step junction (W[=0.040 in. W2=0.020 in., b=0.055 in., er=24). S£ FDTD Mode Matching -5 — Libra -x - Altschuler and Oliner 0 4 12 8 16 20 Freq(GHz) Figure 4-24. Phase of S21 for stripline step junction (Wi=0.040 in. W2=0 . 0 2 0 in., b=0.055 in., er=24). 95 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A relatively large discrepancy occurs for the phase o f the Sn parameter for the three cases shown. The FDTD results and the Mode matching results agreed to within about 5 degrees up to 20 GHz for all three cases. The LIBRA and the Altschuler and Oliner calculations gave a more rapid increase of phase shift with frequency than the FDTD or Mode Matching results. The difference amounted to 30 to 50 degrees o f phase at the highest frequency shown in Figures 4-22,4-18 and 4-14. A new set of empricial formula were derived from the S-parameters calculated using the FDTD method. The formula were based approximately on the formula of Altschuler and Oliner for X and li in (4-17). It was just shown that the existing formula gave reasonable results for the magnitude of the Si i and S21 parameter, and for the phase o f S21, but was inaccurate for the phase of Sn. Therefore, the intent was to keep the functional dependence of (4-17) in tact as much as possible, and to adjust the values of the coefficients using a least squares fit to the FDTD calculations. The first step was to extract the value of the reactance X and the line length h from the FDTD S-parameter data. This was done by taking the absolute magnitude squared o f the Sn in (4-19) and solving for X: (4.20) The value of h was extracted from the phase o f Si 1 as follows: 96 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I“ (p2 - ( p ,- Z S „ 'HI » (4.21) where cp, = tan ' 1 —---- — ' 2 ~ I and q>2 = tan ' 1( — U ,+z 2 The values o f the X and It from (4-20) and (4-21) were then used as the dependent variables to derive the coefficients A,B,C,D,and E in the following formulae using a least squares curve fitting procedure with Zt, Z2, Wj, W2, and b as the independent variables: (4.22) I, = c + d (w , - w 2) + e(w , - w 2)2. The coefficients in (4-22) that were derived from the FDTD data are shown in Table 4-2. The coefficients in Table 4-2 were derived using FDTD data for step junctions over the following range: 3.9 <er<24,0.0018 in. < Wj < 0.04 in., 0.00108 in. < W2 < 0.02 in., and b=0.055 in. The magnitude and phase o f S11 and S21 are plotted in Fig. using the curve fitted equation (4-22) for the stripline step junction with Wl=0.040, W2=0.020, b=0.060, and er=3.9, 9.8 and 24. The original FDTD data is also plotted in Figs. for comparison. Generally, the curve fitted S-parameters match the original FDTD data very closely. Table 4-2. Values of coefficients in (4-22). A B C D E 2.741 1.709 6.794x104 0.4145 -8.5928 97 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 022 Curve Fit FDTD s =24. I- e=3.9 0.195 0.19 0.185 0.18 0 10 20 30 40 50 60 Freq(GHz) Figure 4-25. Magnitude o f Si i for stripline step junction from curve fit and FDTD (W,=0.040 in. W2=0.020 in.,and b=0.055 in.). Curve Fit FDTD BO 3 -12 s= 24. -16 -20 0 10 30 20 40 50 60 Freq(GHz) Figure 4-26. Phase of Si i for stripline step junction from curve fit and FDTD (Wi=0.040 in. W2=0.020 in., and b=0.055 in.). 98 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. .Curve Fit! .FDTD 0.99 «r w ■§ 0.98 e_=24. ee cs s e = 9 .8 0.97 0.96 10 30 20 40 60 50 Freq(GHz) Figure 4-27. Magnitude of S2i for stripline step junction from curve fit and FDTD (W 1=0.040 in. W2=0.020 in., and b=0.055 in.). 0 Curve Fit FDTD el w -2 ■s. <•> © vV) n -4 -5 e.=9.8 8=3.9 e=24. -6 0 10 30 20 40 50 60 Freq(GHz) Figure 4-28. Phase of S21 for stripline step junction from curve fit and FDTD (W,=0.040 in. W2=0.020 in., b=0.055 in.). 99 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Calculation of S-parameters and equivalent circuit for stripline tee-junction discontinuity The FDTD analysis o f the stripline Tee-junction discontinuity follows a set o f procedures that are similiar to those for the open and step discontiuity. In general, a total o f six S-parameters are required to completely characterize the lossless tee-junction discontinuity. Here, for the purpose o f comparing results to existing theory, only the symmetric tee-junction will be considered in which case the the width Wj o f line 1 is equal to W2 of line 2. A diagram o f the tee-junction and the Altschuler and Oliner equivalent circuit are shown in Figure 4-7. For the lossless symmetric tee-junction, there are four unique S parameters: Su (=S22), S33, S21 (=Si2), and S3 i(=S 13= 823 = 832 ). Therefore, the reflected waves from port 1 and from port 3 (S 11,833 ), the transmitted wave from port 1 to port 2 (S 12), the transmitted wave from port 1 to port 3 (S 13), and the transmitted wave from port 3 to port 1 are required from the FDTD analyis. Two simulations o f the Tee junction are required: one simulation with the wave incident to Port 1 and one simulation with the wave incident to Port 3. As before, the reflected wave and transmitted waves are recorded a few cells away from the junction so as not to include any evanescent modes. An example of the incident and reflected waves on line 1, the transmitted wave from port 1 to port 2, and the transmitted wave from port 1 to port 3 for a stripline Teejunction discontinuity are shown in Figures 4-29,4-30, and 4-31, respectively. The incident and reflected voltages on line 3 are shown in Figure 4-32. For the example shown, the port 1 and port 2 line widths Wj and W2 were both equal to 0.040 in., the port 3 line width W3 was equal to 0.020 in., and the ground plane spacing b was 0.055 100 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. in. with a dielectric constant er o f 3.9. The FDTD simulation used a 22x120x200 grid with a cell dimension o f 0.0025 in. in all three dimenstions, with a time step At of 0.1058 psec. No symmetry was used for the Tee-junction calculations, and an ABC was applied at x=0, x=imaxAx, y=0, y= jmaxAy, z=0 and z= kmaxAz. Lines 1 and 2 o f the tee junction were centered about j=60, and line 3 was centered about k=140. The incident and reflected voltages on line 1 were recorded at ji=60 and ki=80, the transmitted voltage from line 1 to line2 was recorded at j 2 = 6 0 and k2= 2 0 0 , and the transmitted voltage from line 1 to line 3 was recorded at j 3 = 1 2 0 and k3=l40. The Tee-junction simulation required about 30 minutes on a Cray YMPel supercomputer. The line 1 reflection coefficient, i.e. the Si t parameter, for the stripline Teejunction is calculated using the Fourier transform of the incident and reflected voltages versus time from Figure 4-29 in (4.10). Likewise, the line 3 reflection coefficient, i.e. the S33 parameter, for the stripline Tee-junction is calculated using the Fourier transform of the incident and reflected voltages versus time from port 3, shown in Figure 4-32, in (4.10). The forward transmission coefficient between line 1 and line 2, 1.e. the S 12 parameter, can be calculated using the Fourier transform o f the incident and transmitted voltages versus time from Figures 4-29 and 4-30, Vjnc((o) and V ^ t a ) in (416) with Li=60Az and L2=60Az. Finally, the forward transmission coefficient between line 1 and line 3, i.e. the S 13 parameter, can be calculated using the Fourier transform o f the incident and transmitted voltages versus time from Figures 4-29 and 4-31, Vj„c(©) and V|jne3(co), in the following equation: 101 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 , Incident 12 S 4 0 Reflected -4 0 400 800 1200 Time Step n 1600 2000 Figure 4-29. Incident and reflected wave on line 1 o f a stripline Tee-junction (Wl=W2=0.040 in., W3=0.020 in, b=0.055 in., er=3.9). 16 12 Transmitted to Line 2 u M =o 8 > 4 0 0 400 800 1200 1600 2000 Time Step (n) Figure 4-30. Transmitted wave from line 1 to line 2 o f a stripline Tee-junction (Wi=W2=0.040 in., W3=0.020 in, b=0.055 in., er=3.9). 102 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 12 Transmitted to Line 3 u oc se > 8 4 0 0 400 1200 800 1600 2000 Time Step (n) Figure 4-31. Transmitted wave from line 1 to line 3 o f a stripline Tee-junction (W l=W2=0.040 in., W3=0.020 in, b=0.055 in., er=3.9). 16 Incident 12 8 4 0 -4 -8 0 400 800' 1200 1600 2000 Time Step n Figure 4-32. Incident and reflected wave on line 3 for a stripline Tee-junction (W,=W2=0.040 in., W3=0.020 in, b=0.055 in., er=3.9). 103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. S3I(oj) = V„ne3 (o>,j = 120,k = l40)eT(" ^ V;nc(co;j = 60,k = 80)e (4.23) where the value o f £/=(k-k[)Az=60Az and Lj= (j 3-j i+Wt/2)Ay=52Ay in (4.23) references the S31 parameter to the Tee-junction reference planes T and T \ indicated in Figure 4-7c. The S-parameters for the symmetric stripline Tee-junction can also be calculated using the equivalent circuit shown in Figure 4-7c and values for the reactances Xa and Xb and for the transformer turns ratio n given which are given by the following formulae from Altschuler and Oliner30: (4.24) sin(7tD3 / A.) TtD/A. where Incosei The equivalent stripline width parameters Di and D3 in (4.24) are given in (4-18). The S-parameters for the stripline Tee-junction are given by:29 104 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Z r2 + 2 X aX b + X a2 - 2 y ^ f - X a S[[ = S22 = ---- SI2 “ S2 l = 5 _ 5 ( Z I + yXa)A (ZI + yXa )A _o _c b l3 ~ s 23 ~ S 3I “ S3I 2 n ’ >/ z 1Z 2 / n 2 r d ^ --------- 7----------> ( 4 .2 5 ) A Zl - 2 ^ - + y ( X a + 2 X b) c _ n S3 3 -------------- -------------- . where A = Zy + 2 Z 3 / n2 + y(X a + 2 X b). The S-parameters for the symmetric stripline Tee-junction are plotted in Figures 4-33 to 4-40 using the FDTD simulation and using the equivalent circuit parameters of Altschuler and Oliner in (4.15). Three separate cases are plotted in the figures for substrate dielectric constants o f er=3.9, er=9:8, and er=24. The widths o f the Teejunction stripline and the ground plane spacing were the same as for Figures 4-29 to 432. For the sr=3.9, sr=9.8, and er=24 cases, the frequency o f the first higher mode given by (4.1) for a 0.040 in. wide stripline was 14.5 Ghz, 22.7 Ghz, and 35.9 Ghz, respectively. At a frequency slighlty above the frequency given by (4.1), the FDTD results became erratic. The S-parameters in the figures are shown roughly up to the frequency of the first higher mode. 105 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 03 FDTD Altschuler and Olinerj 02 I 0.15 0.05 0 10 40 30 20 Freq(GHz) Figure 4-3 3. Magnitude of S 11 for symmetric stripline Tee junction (Wi=W2=0.040 in., Wj=0.020 in., b=0.055 in.). -160 j ___ FDTD 1 .......Altschuler and Oliner i -170 e=24. er=9.8 6 r= 3J *9 7 \ ee A -a -180 Vi 1mm O aCB -190 -200 6=24. \ er=9.8 6=3.9 -210 10 20 30 40 Freq(GHz) Figure 4-34. Phase of S 11 for symmetric stripline Tee junction (Wt=W2=0.040 in., W3= 0 . 0 2 0 in., b=0.055 in.). 106 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.95 0.9 / 0.85 Zf e 4t *r=24. 0.8 0.75 DO eg S 0.7 0.65 FDTD Altschuler and Oliner 0.6 0 10 20 40 30 Freq(GHz) Figure 4-35. Magnitude of S21 for symmetric stripline Tee junction (W!=W2=0.040 in., W3=0.020 in., b=0.055 in.). 16 FDTD Altschuler and Olinerl 12 BO ■oo» ' e =9.8 e=24. 6= 3.9 o eg St a. 4 0 0 10 20 30 40 Freq(GHz) Figure 4-36. Phase of S21 for symmetric stripline Tee junction (Wi=W2=0.040 in., W 3=0.020 in.. b=0.055 in.). 107 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FDTD Altschuler and Oliner 0.6 ■o 0 .5 er=9.8 0 .4 0J 0 20 10 40 30 Freq(GHz) Figure 4-37. Magnitude of S31 for symmetric stripline Tee junction (W,=W2=0.040 in., W3=0.020 in., b=0.055 in.). B 01B 2 -10 zf IV ) -15 e=24. -20 6 = 3 .9 *25 ] -30 I 0 FDTD ' Altschuler and Oliner i 20 10 30 40 Freq (GHz) Figure 4-38. Phase o f S31 for symmetric stripline Tee junction (Wi=W2=0.040 in., W3= 0 . 0 2 0 in., b=0.055 in.). 108 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FDTD Altschuler and Oliner <■*» e « S s O «fl ■8 0.8 e=3.9 I 0.7 s 0.4 0 10 20 30 40 Freq(GHz) Figure 4-39. Magnitude of S33 for symmetric stripline Tee junction (Wi=W2=0.040 in., W 3=0.020 in., b=0.055 in.). 200 FDTD Altschuler and Oliner 180 160 140 . xt a. e =24. 120 100 0 10 20 30 40 Freq(GHz) Figure 4-40. Phase o f S33 for symmetric stripline Tee junction (Wi=W2=0.040 in., W3=0 . 0 2 0 in., b=0.055 in.). 109 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Generally, the shape of the S-parameter curves were consistent between the FDTD data and equivalent circuit model. The magnitude of the S-parameters from the FDTD results and from Altschuler and Oliner’s equivalent circuit model agreed closely at low frequency, but differed by about 10 to 2 0 % at frequencies near the first higher mode. The phase o f the S21 parameter in the two cases differed by only a few degrees. However, the phase o f the Sn, S31, and S33 parameters differed by 10 to 30 degrees above a few gigahertz. Measurements will be presented later in this chapter that support the accuracy of the FDTD results. An attempt was made to derive new values for the reactances Xa ,Xb and the transformer turns ratio n in the equivalent circuit model of Altschuler and Oliner. However, the new circuit paramters derived from the FDTD results failed to accurately replicate the S-parameters. It was found necessary to develop a new equivalent circuit model with four reactances in order to accurately represent the four S-parameters. The revised equivalent circuit for the Tee-junction discontinuity is shown in Figure 4-41. 110 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. w ( T W2 T T Figure 4-41. New equivalent circuit for symmetric Tee-junction. The new equivalent circuit in Figure 4-41 has four unknown reactances Xa, Xb, Xc, and Xj. The values of the reactances were specified in terms o f the impedance [z] matrix. The impedance is derived from the S-matrix [S] following a procedure similar to one described by Liang.65 The [S]-matrix is formed from the Fourier-transformed Sparameters determined from the FDTD simulation, where [s]« S ll S12 S I3 ^21 _S3l ^22 ^23 S3 2 S 3 3 (4.26) The three port Z-parameters are computed from [S] according to (4.27) where [I] is the identity matrix. The values o f the reactances X*, Xb, Xc, and Xd in the equivalent circuit o f Figure 4-41 are given by: jX (4.28a) c = Z 3i yXb= z2i -yXc (4.28b) 111 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. jX a = Z[ I - j X b-yXc (4.28c) yXd = z33 -jX c (4.28d) The next step was to fit the reactances in (4.28) to curves which were functions o f the stripline parameters and o f the frequency g j . The functional form for the reactances was largely developed on a trial and error basis. The reactance equations in (4.24) for Xa and Xb could be used as starting points, with adjustments made until the new equations fit the FDTD data. The functional form o f the reactance equations for Xc and Xd were based solely on trial and error until a good curve fit was obtained. This involved viewing the reactance data on a plot versus one o f the stripline parameters, such as width or impedance, and trying different families of functions, usually a polynomial of the stripline parameters, having a similar functional dependence. A particular functional family would be expressed in term a set o f unknown coefficients or parameters and tested using a least squares fit to the FDTD data to determine the “best” fit. For the stripline Tee-junction reactances, the following functions were fit to the FDTD data: X ~ L x z = ~ " [a Z.2 iu + A 2 u 2 + A3U3 + A 4u 4 + A 5U 5] , = y e ^ B , + B2 v 2 + B3 v 3 + B4 v 4 + B5vs + B 6 v 6 + B 7 v 7 + B 8v8], X X ~ A, 5r [c 1y+c 2y1+c ,y3]+c ,v +c 5v2 + c 6v3+c 7v4, b = —[ d , v + D ,v 2 + D 3v 3 + D 4v 4 + Dsw + D 6w 2 + D 7w 3 j, . W3 b W3 Z, whereu-T ,v = - , w = — . a n d y ^ . 112 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.29) The parameters Aj, Bj, Q, and Di are listed in Table 4-3. The S-parameters from the curve fit equations above are plotted along with the S-parameters from the FDTD results in Figures 4-42 through 4-49 for a stripline substrate with a dielectric constant er=3.0, W|=0.008 in. and W3 as noted in the figures. This curve fit was used in the design o f the low pass filters that are described in the next section. The agreement between the curve fit and the FDTD data was generally very good. The magnitude of the S-parameters using the curve fit was within about 1-3% of the FDTD calculations, and the phase o f the S-parameters using the curve fit was within 1-2% degrees of the FDTD calculations. The curve fit equation in (4.29) applies to Tee-junctions with Wi=0.008 in. and W3 between 0.008 in. and 0.064 in. with a dielectric constant o f 3.0. For other dielectric constants and line widths, similiar curve fits to the FDTD data provided similar accuracy. Table 4-3. Set o f Coefficients for Equation (4.29). i 1 2 3 4 5 6 7 8 9 Ai -0.50637 -0.33471 6.976533 -27.4189 52.0074 Bi • Q -6.93E-02 -3.64E-02 8.03161 0.149709 -149.264 -3.68E-02 1369.888 0.636793 -6756.13 20.29001 18571.89 -92.4422 -26753.8 138.8273 15797.35 0.922449 Di 11.37851 -86.1381 248.0138 -255.806 -1.37657 1.18E-02 -1.21E-02 113 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.6 FDTD Curve fit e=3.0 © © 3 B O «K £ 0.3 0.2 - 10 0 20 30 40 50 Freq(GHz) Figure 4-42. Comparison of curve fit for magnitude o f Si i to FDTD results (W[=W2=0.008 in., b=0.060 in., er=3.0). -120 FDTD Curve fit -140 W3=0.032 in. £ -160 ♦ -»t e v <n « -180 . ,=0.008 in. -200 0 10 20 30 40 Freq(GHz) Figure 4-43. Comparison of curve fit for phase o f Su to FDTD results (W[=W2=0.008 in., b=0.060 in., er=3.0). 114 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 0.8 0.7 W, =0.064 0.5 FDTD Curve fit 0.4 20 10 0 30 50 40 Freq(GHz) Figure 4-44. Comparison of curve fit for magnitude o f S2 t to FDTD results (Wi=W2=0.008 in., b=0.060 in., sr=3.0). 16 i FDTD | Curve fit! 12 4M 1 e 8 a 4 0 0 10 20 30 40 Freq(GHz) Figure 4-45. Comparison of curve fit for phase of S21 to FDTD results (W i=W2=0.008 in., b=0.060 in., er=3.0). 115 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 0.72 FDTD Curve fit e=3.0 W,=0.032 in. 0.68 e «BS S 0.64 W,=0.008 in. 0.6 0 10 20 30 40 50 Freq(GHz) Figure 4-46. Comparison of curve fit for magnitude o f S31 to FDTD results (W[=W2=0.008 in., b=0.060 in., er=3.0). 40 FDTD i W ,=0.064in. Curve fit ; 30 BS ■U o O Si 5inS A 0. W,=0.032 in. 20 10 W3=0.008 in. 0 0 10 20 30 40 Freq(GHz) Figure 4-47. Comparison o f curve fit for phase o f S31 to FDTD results (Wi=W2=0.008 in., b=0.060 in., er=3.0). 116 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 0.6 FDTD e= 3.0 Curve fit 0.45 W ,=0.016in. B cs s 00 0 20 10 40 30 50 Freq(GHz) Figure 4-48. Comparison of curve fit for magnitude o f S3 3 to FDTD results (W 1=W2=0.008 in., b=0.060 in., er=3.0). 40 ----- FDTD ....... Curve fit ! -a 5* «S* O 4«> : -40 ; ! ’ a. e= 3.0 i W3=0.064 in. | -80 i ; 55 .S i | -120 i I W j=0.0l6 in. ^ __ W3=0.032 in. -160 - r ---------------i -200 10 20 ' I i W3=0.008 in. J ; 30 40 Freq(GHz) Figure 4-49. Comparison of curve fit for phase of S33 to FDTD results (W,=W2=0.008 in., b=0.060 in., er=3.0). 117 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 Design of low pass filters including stripline discontinuities la this section the design o f low pass stripline filters is described that is based on a procedure by Roan.32 The essence of this method is to iteratively adjust a set o f the filter’s element values until the zeros, poles, and the scale factor o f the characteristic function o f the filter match their desired values. The method requires only simple analysis, avoids complications normally associated with network synthesis o f highdegree networks, and easily takes into account the effects o f the discontinuities and parasitics. The design o f several low pass stripline filters will be designed using the equations of the previous section for the openend, step and stripline discontinuities. The low pass stripline filters are based on a lumped element prototype that is shown in Figure 4-50. The filter is designed to have a maximum pass-band attenuation of a p up to a specified pass-band cutoff frequency 4>. The lumped element prototype filter of order n posseses a passband ripple equal to Op, with transmission zeros at frequency 4k and reflection zeros at frequency 4k* For convenience, it will be assumed that the filter orden n is odd («=2N+1) and that it possesses the m axim um number o f finite transmission zeros (i.e., M=N). Under these assumptions the cicuit will be symmetric (i.e.,the two port parameters Z n = Z^). A layout o f a stripline circuit with ir=7 which approximates the behavior of the lumped-element prototype is shown in Figure 4-51. This circuit consists o f a cascade o f short lengths o f high-impedance stripline alternating with shunt open-circuit stubs. The lengths o f the shunt stubs are chosen such that they produce the transmission zeros, 4k- The pass band and stop band 118 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. transmission loss in dB is shown in Figures 4-5 la and 4-52b, where the frequency is normalized to the pass band cutoff frequency fra. i m iT T L .. jnm. j m 50 Q V/ 50 Q (a) R eflection Zeros 0 a e < tns» § i H a ,p -0.3 0 0.6 0.4 0.8 12. Normalized Frequency (b) Figure 4-50 (a) Schematic diagram o f lumped element low pass filter, (b) Pass band response (continued). 119 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -10 aa -20 e© m -30 •o •» Transmission zeros E m B b« H -40 -50 -60 I 3 2 4 Normalized Frequency (c) Figure 4-50. (Continued) (c) Stop band response. Z 3 /3 Z 5A Figure 4-51. Layout of seven element stripline low pass filter. Three 7 element (3-pole) and one 11 element (5-pole) low pass stripline filters were designed, manufactured and measured to demonstrate the accuracy of the equivalent circuits for the stripline discontinuties developed in the previous section. A n=2N+l element stripline low pass filter has n+1 discontinuities: one openend 120 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. discontinuity and one tee-junction discontinuity for each o f the N shunt branches, and two step junction discontinuities at the input and output The effects o f the discontinuities are included in the filter synthesis procedure using the equations for the reactances o f the equivalent circuit. The three 7 element filters had a pass band cutoff frequencies of f^ equal to 4 GHz, 8 GHz and 12 GHz. The 11 element filter had a pass band cutoff frequency fCT o f 8 GHz. The pass band ripple a p, the normalized pass band reflection zeros and the stop band transmission zeros were identical for the three 7 element filters, and are shown in Table 4-4. The normalized pass band reflection zeros and the stop band transmission zeros for the 11 element filter is shown in Table 4-5. The impedance, width, and length of each of the sections for the four filters are shown in Tables 4-6 and 4-7. A 0.060 in. ground plane spacing was assumed for the stripline with a dielectric constant equal to 3.0. Table 4-4. Normalized reflection zeros and transmission zeros o f a seven element filter with Op=0.1773 dB. Normalized Reflection Zeros: 0.98604561 0.8586568 0.5355906 Normalized Transmission Zeros: 1.880000 1.230000 1.500000 121 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table 4-5. Normalized reflection zeros and transmission zeros of a eleven element filter with ap=0.l773 dB. Normalized Reflection Zeros: 0.9469231 0.9944065 0.8382131 0.6463793 0.3585913 Normalized Transmission Zeros: 2.248000 1.600000 1 .2 0 0 0 0 0 1.300000 1.885000 Table 4-6. Characteristic impedances and line lengths and widths o f seven section filter for 4, 8 and 12 GHz cutoff frequencies. 4GHz cutoff 8 GHz cutoff 12 GHz cutoff Z length width Section Z Z length width length width No. (ohms) (in.) (in.) (ohms) (in.) (in.) (in.) (ohms) (in.) 1 1 0 0 .0 0 0.1699 0.0062 1 0 0 .0 0 0.0855 0.0062 1 0 0 . 0 0 0.0576 0.0062 2 50.58 0.2336 0.0348 50.88 0.1186 0.0344 50.41 0.0800 0.0350 3 1 0 0 .0 0 0.2506 0.0062 1 0 0 .0 0 0.1264 0.0062 1 0 0 . 0 0 0.0853 0.0062 4 96.05 0.3465 0.0071 93.44 0.1719 0.0079 91.44 0.1139 0.0085 5 1 0 0 .0 0 0.2312 0.0062 1 0 0 . 0 0 0.1148 0.0062 1 0 0 .0 0 0.0760 0.0062 6 70.01 0.2853 0.0179 68.58 0.1420 0.0188 67.15 0.0945 0.0197 7 1 0 0 .0 0 0.1451 0.0062 1 0 0 .0 0 0.0709 0.0062 1 0 0 . 0 0 0.0460 0.0062 Table 4-7. Characteristic impedances and line lengths and widths of eleven section filter for a 8 GHz cutoff frequency. Section No. Z (ohms) 1 1 0 0 .0 0 2 42.76 3 4 5 1 0 0 .0 0 6 91.67 length (in.) 0.0860 0.1068 0.1564 0.1371 0.1267 0.1763 7 1 0 0 .0 0 0 .1 1 0 2 8 76.12 9 1 0 0 .0 0 10 50.18 11 1 0 0 .0 0 0.1632 0.1317 0.1188 0.0772 53.93 1 0 0 .0 0 width (in.) 0.0062 0.0462 0.0062 0.0309 0.0062 0.0084 0.0062 0.0145 0.0062 0.0353 0.0062 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The calculated filter response from the new circuit model is shown in Figures 452 to 4-54 for the seven section filters with 4, 8 and 12 GHz cutoff frequencies, and in Figure 4-55 for the eleven section filter with a 8 GHz cutoff frequency. All filters were designed for a maximum pass band attenuation o f 0.1773 dB, and for a roughly equal ripple in the stop band. The minimum attenuation in the stop band was about 40 dB for the seven element filter and 70 dB for the eleven element filter near the band edge. Due to the fact66 that the z-parameters o f a length o f transmission line vary in proportion to either the cotangent(pl) or sin(|31), the response of the stripline filter is not optimum in the sense o f the lumped element prototype. As a result, the stripline filters have a unequal ripple in the pass band and the stop band attenuation decreases above about twice the passband cutoff frequency. In the case of the eleven section filter, an attempt was made to increase the attenuation in the stop band at frequencies above double the passband cutoff frequency by placing the two transmission zeros (poles) at a frequency where the stop band attenuation would otherwise display a rapid decrease. In Figure 4-55, this occurred at frequencies o f 15.080 and 17.984 GHz. It was observed however, that small shifts in resonant frequency would result in a very rapid decrease in stop band attenuation with frequency. For comparison purposes, the filter response in Figures 4-52 to 4-55 is also shown calculated by a commercial CAD package called LIBRA. The same filter dimensions were entered into the LIBRA model as were listed in Tables 4-5 and 4-6. In 123 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the case of the 4 GHz cutoff filter, the LIBRA model agreed fairly well with the circuit model developed here up to a frequency of about 8 GHz. Above 8 GHz, the discrepancy increased. The discrepancy between the two models also increased as the cutoff frequency increased. Generally, LIBRA predicted higher frequecies for the transmission zeros. Particularly in the case of the 5 pole filter, the effect of shifting the two highest frequency transmission zeros is readily observed. 124 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Circuit M o d e l LIBRA -20 A/, CO -40 -60 J 2 0 4 3 5 Frequency (GHz) (a) -20 n C/5 Circuit M o d e l -40 LIBRA -60 4 6 10 8 12 14 Frequency (GHz) (b) Figure 4-52. Calculated (a) passband S,, response and (b) stop band S2I response for a seven element stripline filter with a 4 GHz cutoff frequency. 125 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 0 Circuit M o d e l LIBRA j -20 S' •a VI -40 -60 0 4 2 6 8 10 Frequency (GHz) (a) -20 CQ ■ o Vi -40 New Model LIBRA -60 8 16 24 32 Frequency (GHz) (b) Figure 4-53. Calculated (a) passband S„ response and (b) stop band S2l response for a seven element stripline filter with a 8 GHz cutoff frequency. 126 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0 Circuit M o d e l LIBRA -20 a -40 -60 0 4 8 16 12 Frequency (GHz) (a) -20 CO -o cnM -40 New Model -60 LIBRA J 12 20 28 36 Frequency (GHz) (b) Figure 4-54. Calculated (a) passband Stl response and (b) stop band S21 response for a seven element stripline filter with a 12 GHz cutoff frequency. 127 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0 Circuit M o d e l LIBRA WJ —L. -20 . B cn -40 -60 0 4 2 8 6 10 Frequency (GHz) (a) -20 CQ T3 - -40 'w' -60 -80 Circuit Model LIBRA -100 8 12 16 20 24 Frequency (GHz) (b) Figure 4-55. Calculated (a) passband S„ response and (b) stop band S2I response for an eleven element stripline filter with a 8 GHz cutoff frequency. 128 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Comparison of predicted filter response with measurements The low pass stripline filters discussed in the last section were fabricated on a low loss, high frequency circuit board material (Rogers Corporation high frequency circuit material, Type R03003) with a dielectric constant er of 3.0±0.04, a dielectric thickness o f 0.030+0.04 in., a dissipation factor of 0.0013 and with 1 ounce copper (about 0.00135 in. thick) on both sides. The bottom half o f the stripline filter circuits were etched on one side of either 1.0 in. x 1.0 in or 1.5 in. x 1.5 in. boards using standard printed circuit board processing techniques, while allowing the copper metal on other side to remain intact. The filter circuits were processed to a dimensional accuracy o f about 0.0005 in. A second 1.0 in. x 1.0 in or 1.5 in. x 1.5 in. board was processed for the top half of the filter by removing all o f the metal from one side only. The circuit boards were then pressed together between two 0.25 in. thick aluminum blocks (unetched metal side out) held together by four #10 screws. Petroleum jelly (nom. e =2.3) was used to fill the air gap that normally would occur between the two boards due to the 0.00135 in. thickness of the copper filter elements. Connections to the center conductor were made using SMA type panel mount (Southwest Microwave Model 101 l-02sf) connectors that were attached to the 0.25 in. thick metal blocks using four screws. The center pin of the SMA connectors had a short 0.012 in. wide x 0.005 in. thick tab that was soldered to the 50 ohm (0.035 in. wide) lines used for the input and output lines of the filters. The filter response was measured using a Hewlett Packard Model 8510C vector network analyzer. 129 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The measured and calculated magnitude o f the S,, and S2, responses for the four filter circuits are shown in Figures 4-56 to 4-59. The calculated response is based on the circuit model and is identical to that shown in Figures 4-52 to 4-55. The agreement between the measured and calculated response was very good for the 4 GHz cutoff and the 8 GHz cutoff seven element filters. The reflections zeros for measured Sn response were somewhat different that that for the circuit model, but the average S,, levels were in good agreement. The measured transmission zeros in the stop band were in excellent agreement with the circuit for the 4 GHz and the 8 GHz seven element filters. The measured Stl response was reasonably close to the circuit model for the 12 GHz seven element filter and the 8 GHz eleven element filter. The differences were largest at the frequencies in the pass band closet to the cutoff frequency. In the stop band, the transmission zeros at the lower band edge were in good agreement with the circuit model, however, some discrepancy (3-5%) was observed at the higher frequencies. 130 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. — Measured (dB) Circuit mode!(dB) -20 & -40 -60 J 0 I 3 2 4 5 Frequency (GHz) (a) -20 -40 Circuit m o d el Measured -60 4 8 12 16 Frequency (GHz) (b) Figure 4-56. Measured and calculated (a) S„ response and (b) S21 response for a seven element stripline filter with a 4 GHz cutoff frequency. 131 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Measured Circuit model -20 SB -40 -60 0 4 2 6 8 10 Frequency (GHz) (a) 0 -20 -40 Circuit m o d e l Measured -60 8 16 12 20 24 Frequency (GHz) (b) Figure 4-57. Measured and calculated (a) Slt response and (b) S2l response for a seven element stripline filter with a 8 GHz cutoff frequency. 132 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Measured Circuit model i -20 S' •o -40 -60 0 4 12 8 16 Frequency (GHz) (a) Circuit model Measured -20 ^N Q fi •O -40 . -60 -80 12 20 16 24 Frequency (GHz) (b) Figure 4-58. Measured and calculated (a) Su response and (b) S2I response for seven element stripline filter with a 12 GHz cutoff frequency. 133 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Measured Circuit model -20 a■oa C O -40 - -60 +r J 0 4 2 6 8 10 Frequency (GHz) (a) -20 -40 CQ ■a -60 -80 Circuit model — Measured -100 8 12 16 20 24 Frequency (GHz) (b) Figure 4-59. Measured and calculated (a) SMresponse and (b) S2I response for an eleven element stripline filter with a 8 GHz cutoff frequency. 134 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Complete analysis of low pass filters using FDTD A complete analysis of the 12 GHz cutoff seven element filter and the 8 GHz cutoff eleven element filter was accomplished using the FDTD method. The seven element filter was analyzed using a uniform 12x160x360 cell grid with Ax=0.005 in., Ay=0.002 in., and Az=0.002 in. A 12x80x180 cell grid was used for the eleven element filter with a uniform grid dimension in the x-direction (Ax=0.005 in.), but with a nonuniform grid dimension in the y- and z-directions. For non-uniform grid in the y- and zdirections, the individual cell dimensions were adjusted in such a way that a integral number o f cells fit exactly across the width and lengths o f the filter elements, as illustrated in Figure 4-60. In the case o f the uniform grid, the individual element dimensions were accurate only to ±0.001 in. The average cell dimension for the nonuniform grid was 0.00633 in. in both the y- and z-directions (plane of the stripline). The non-uniform grid required fewer cells and allowed the use o f a larger time step, which significantly reduced computation time without sacrificing accuracy. 131,072 time steps were used for the uniform grid with a time step o f 0.0847 picoseconds, while 65,536 time steps were used for the non-uniform grid with a time step of 0.268 picoseconds. The uniform grid required 8 hours o f cpu time, while the non-uniform grid required only about 2 hours of cpu time. A comparison o f the measured S,, and S2I filter responses to the responses calculated using a FDTD analysis is shown in Figures 4-61 and 4-62 for the seven element and the eleven element filters, respectively. For comparison, a few points o f the filter response calculated using a Finite Element Method simulator, called HFSS ,67 is 135 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. shown in Figure 4-61, for frequencies between 10 and 20 GHz. The HFSS data was collected in two program runs from 10 to 13.4 GHz in 0.1 GHz steps, and from 13.6 to 20 GHz in 0.2 GHz steps. The second HFSS data run was optimized by increasing the number of finite elements (mesh size) relative to the first run. A filter response from 0 to 24 GHz was not carried out using HFSS since approximately 30 minutes o f computer time were required for every frequency point from 10 to 13.4 GHz, and approximately 2 hours o f computer time were required for every frequency point from 13.6 to 20 GHz. Overall, the agreement between the measured response and the response calculated using FDTD is excellent. The HFSS results indicate a lower cutoff frequency than the FDTD results or the measurements, but additional finite elements may further reduce the discrepancy. The discontinuity between 13.4 and 13.6 GHz shows that the HFSS results had not converged during the first data run. The measured frequencies o f the reflection zeros in the pass band S tI response agree closely with those calculated using FDTD, although the magnitude o f the Su response is slightly higher for the measurements. The agreement between the measured S2I response and the FDTD results is excellent. 136 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 1 1 1 1 1 1 1 t 1 1 1 t 1 1 1 1 1 1 1 1 1 t 1 1 t I 1 1 1 1 1 t 1 1 1 1 1 1 1 1 1 1 1 1 1 1 t 1 1 1 1 1 1 1 1 1 1 1 1 t t 1 1 1 1 1 1 1 1 1 1 1 t 1 1 1 1 1 1 * 1 1 1 1 I •• 1 1 « 1 1 _» 1 1 1 1 1 1 1 1 1 1 t 1 t 1 l t 1 1 t 1 l t • J _ l 1 1 1 1 1 i 1 l I i i l l I 1 i I t l 1 1 l 1 l l f 1 1 1 1 1 1 1 l i i J i 1 l 1 1 t l 1 1 1 1 1 t I I 1 1 1 1 1 t 1 I 1 1 1 1 1 1 1 1 1 1 1 1 « 1 I 1 1 t 1 1 I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 I 1 1 1 1 1 I I 1 1 1 1 1 1 1 1 1 • J - .- . 1 1 I Figure 4-60. Illustration o f a non-uniform FDTD grid. The filter response calculated using FDTD was in better agreement with the measurments than the response calculated using the circuit model. The most noticable differences occur near the frequency of the two highest transmission zeros in the stop band. The circuit model provided good accuracy at the other pole frequencies. It turns out that the shunt sections corresponding to the two highest transmission zeros are in close proximity to the step junction at the input and output to the filter. The shift in resonant frequency that occurs in both the measurements and the FDTD calculations is believed due to coupling between shunt sections at either end o f the filter and the wide 50 ohm input and output lines. The other shunt sections are unaffected because o f a wider separation from the input and output lines. 137 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. M easured FDTD HFSS -20 \/ -40 -60 8 4 0 12 16 Frequency (GHz) (a) — - \ ------------ Measured ________ FDTD ____HFSS i fj i -20 \ v . .v V T v L V \ •' \ U\ T V 1 „ ' -40 ' -60 V '< F I » • 1 t l > - * 1/ i i M l \rl| * * i * ' - — % — .1 1 ' 'I 1/ f W U ■ 1 \ \ ----------- c/i \p \t , 1 ! / »* 1II1 «i -80 12 20 16 Frequency (GHz) (b) Figure 4-61. Comparison o f response calculated using FDTD to measured response for a seven element filter with a 12 GHz cutoff frequency: (a) Su response and (b) S2I response. 138 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 24 -20 a ■o w -40 — Measured — FDTD -60 J 0 4 2 8 6 10 Frequency (GHz) (a) -20 /•“s m -o -40 -60 -80 - Measured : — FDTD -100 8 16 12 20 Frequency (GHz) (b) Figure 4-62. Comparison of response calculated using FDTD to measured response for a eleven element filter with a 8 GHz cutoff frequency: (a) S,, response and (b) S2I response. 139 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 24 Chapter 5 Application of the FDTD Method to Analysis o f Biconical Antennas One o f the advantages o f the FDTD method is the capability to analyze electromagnetic devices over a wide frequency bandwidth from a single time domain simulation. The biconical antenna is capable o f producing an omni-directional radiation pattern over more than an octive bandwidth. Many antenna textbooks discuss the infinite cone biconical antenna or the extreme case of small cone angles or cone angles near 90 degrees,37,3®'39 but little information is available regarding the calculation o f the antenna gain or radiation pattern from biconical antennas with cone angles in the mid ranges (0 °« v |/« 9 0 °). In this chapter the radiation patterns from biconical antennas are calculated using the FDTD method and compared to measurements. The effect of an isolation network on the radiation pattern o f a biconical antenna will also be analyzed using FDTD, and compared to measurements. Calculation of near fields for an axi-symmetric biconical antenna The biconical problem will be solved on a FDTD grid in a cylindrical coordinate system. Since the biconical antenna is axi-symmetric, variations in the azimuthal coordinate 4» are set equal to zero, which reduces the FDTD grid to two dimensions. To reduce the number o f calculations further by one-half, an H-pIane o f symmetry (PMC boundary condition) will be assumed in a plane perperdicular to the axis of the cone and located at the apex of one of the cones. A biconical antenna with a half-cone angle of \j/ is shown in Figure 5-1 along with a illustration of the two dimensional FDTD grid used to represent one-half o f the 140 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. antenna about a plane o f symmetry. Since the FDTD grid is made of rectangular cells, the conical conductor o f diameter a and height h will be approximately represented by a staircase o f conductive Yee cells. Staircased Conductor Symmetry Plane Figure 5-1. (a) Diagram of biconical antenna and (b) FDTD representation o f conical antenna on a two dimensional grid. The FDTD update equations for a two dimensional grid in cylindrical coordinates are as follows: l H"+1/2 (i,k + i ) = H"-I/2a k + i ) + - r At E 4n( i,k - H ) - E ;( i,k )^ Az 141 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (5.1a) At H J * " (i + J , k + i)= H ;-“ (i + i k + i)+ I. 1. E"(i + l,k + - ) - E " ( i , k + —) 1 1 E"(i + - , k + l ) - E ; ( i + - , k ) Ap Az 1 n-I/2/ (5.1b) At 1 p(i + ^ ,k ) ' - (5.1c) 1 1 Ap a + i) At (i + 7 , k ) = E “(i + - , k ) + e(i + - , k ) (5.1d) H «, 0+1/2(i + y , k - j ) - H * n+,/2(i + y , k + - Az E"+l(i,k ) = E j(i,k) + At e(i,k) Az n+I/2/; i . (5-le) i.Y Ap E r ‘( i ,k + - ) = E ; (i, k + 5 ) + E (. k) r 2( i - |. k + i) Ap 142 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (5. If) In (5.1a-f), it was assumed that the conductivity o=0 and that the Yee index in the radial or p-direction runs from i= 0 to i=imix and the index in the z-direction runs from k= 0 to k=kmax. The variation in the ^-direction is equal to zero. On the surface o f the biconical antenna, PEC conditions were applied by setting the permittivity e o f the tangential Efields to some very large value. At p=o, an auxiliary equation is necessary for Ej since the H* field in (5.1f) is undefined at i=-l/2. The auxiliary is equation is based on the integral form o f Ampere’s Law with o=0: I—► -*• d r~* A f H J / = — jEndS. (5.2) Based on (5.2), the auxiliary update equation for Ez at p=0 is as follows: 1 I At 2-H "+I/2(0 + ^ k + ^ ) E"+l(0, k + —) = E"(0, k + —) + --------------------------- 2 ------ 2 2 zv 2' e(0 ,k ) Ap (5.3) -I The E- and H- fields are solved using the update equations (5.1) and (5.3) together with a set of boundary conditions and a specified source. As discussed above, a PEC boundary is applied at z=0 which is a plane o f symmetry. No boundary condition is needed at p= 0 since only E* is calculated there using the auxiliary update equation (5.3). Open boundary condtions are required at p ^ ^ A p and z=kmaxAz. The boundary at k ^ can use an equation similar to the 1st order Mur radiation boundary condition discussed in “Boundary conditions” o f Chapter 2, but adapted to the cylindrical coordinates as follows: 143 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Er‘(i + T . k« « ) = E0”(i + ^ , k lrai- l ) + cA t-A z 1 k m ax -l)-E ^ (i + - , k max) [ E r ' ( i + | . cAt +Az Er'(i. = EJ(i, k„„ - 1)+ ^ 7^ [E ;*'(i, k„„ -1) - EJ(i, kra„ )] (5.4a) (5.4b) However, the 1st order Mur radiation boundary conditions cannot be used at P=iniaxAp in a cylindrical coordinate system. Instead, the 1st order Bayliss and Turkel boundary condition can be used. The l5*order Bayliss and Turkel boundary condition is based on an annhilation operator o f the following form:56 \_ d _ d_ l vPh 5t dp _1 2p) Ez =0 (5.5) The boundary condition (5.5) is applied to both of the tangential fields Ez and E^ at p=imaxAp. So that centered differences can be used for the derivatives in (5.5) to maintain the 2nd order accuracy of the FDTD equations, (5.5) should be evaluated about point p=(iniajt-l/2)Ap and t=(n+l/2) At as follows: E ;<l( w . k * | ) - E r l( w - U Q ^ — n + l/2 ," d p Ez (^max 1 ^ )~ Ap 1 E"(lmax» k + - E"(imax -1 , k + ~) Ap 144 i Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (5.6a) 1. E z^ C im ax- U k+~)•-E"(im a x- U k+^) At ir E r ^ -i k4)=i a 2 2 I I r (5.6b) 2J At When (5.6) is substituted in (5.5) and solved for E^ Vl , . J c-H/2 ), the resulting ABC is: 1. f 4imaxAp —(4imax + l)vphAt' E"+l(imax,k + ^-) = E :( w ,k + | ) k 4 im«Ap+(4i m ax + i)vphAtJ "4(imax - l)Ap - ( 4 W - 3)vphAt" E T 'C i - . - l . k + i ) , 4imaxAp + (4 imax + l)v phAt , (5.7) ' W nm - QAp - ( 4 W - 5 ) v phAtN v. 4 ^max^P + (4i m ax J E "(imax - U k + - ) - + 0 v phAt ,n+l,. An equation to update E^ (im ^ ) is found by substituting E4 at p=(imax)Ap and z=kAz for Ez in (5.7). A sinusoidal source equation with an Guassian envelope, similar to (3.2), is applied to Ez at the apex of the bicone as follows: \2 ~ \ 1 E "[°,° + ^ ) = sin [o 0 ( n - n 0 )At]exp I 0 d d f n, J _ (5.8) As discussed previously, the time constants n, must be chosen to provide sufficient bandwidth, and Oq> 4n, to avoid undesirable transients. The near H^-fields for a one-half o f a biconical antenna with i|/=60° and a=1.0 in. is shown in Figure 5-2. The H4-fields were calculated using a 160Ap byl60Az FDTD grid with Ap=Az=0.0125 in. The time step At=0.5292 picoseconds and n<,=256 145 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. and nt=32, which equate to a pulse bandwidth o f about 29 GHz. The pulse is shown at time step n=356 in Figure 5-2, which is just after the peak reaches the outer edges o f the cone at i=80 and k=48. The fields have just begun to curl around the top edge o f the bicone. Figure 5-2. near-field around a biconical antenna at time step n=356. Near-field to far field transformation for determination of antenna patterns The goal here is the calculation o f the far field antenna patterns at a particular frequency. Calculation of the far-field radiation patterns requires a near-to-far field conversion o f the fields given by the FDTD update equations in (5.1) and (5.3). This conversion is based on the surface equivalence theorem68 whereby the far-fields can be 146 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. calculated using the virtual surface currents 7 s=n x ft and 7 m= If x n over a closed surface surrounding the antenna. It will be necessary to transform the time domain fields to phasor form in order to calculate the radiation patterns at a particular frequency. As before, the discrete Fourier transform will be used for this purpose. The far fields will be calculated from a pair of vector potentials defined as follows: where ? = rr is the observation point or position (x,y,z) where the far-field is to be calculated (pattern point), =r r' is the source point on S'(x',y',z?), lt= R R =T-T ', the angle between ? and and k is the wavenumber. A far field approximation was is used for R on the far right hand side o f (5.9) where R=r for amplitude variations and R^r-r'cos \|/ for phase variations. The electric field and magnetic field I t in the far field are given in terms of the vector potentials in (5.9) as follows: E = -yco A + p V ( V - A ) e0 VxF, (5.10) 0 H = -/co p + ^ V . F ) En V x A. 147 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The far fields will be evaluated in spherical components. Taflove56 shows that if the fields that fall off as 1/r2 or faster in (5.10) are ignored, and if the radial field components that are small in comparison to the 0 and <f>components are neglected, then the far fields from (5.10) are as follows: Er =0, Ee s - /© ( a 0 + iioF*) = ~ ^ 4~~r (L» + ^oNe), E* = - / o ( a # - n oF0) = ^ ~ ( L0 + ^oN*), Hr =0, (5.H ) where N = J j j se|kr,c0Svds', S L = J p m=lb'“ s,,ds- and s I— (5.12) Vo=J— V e 0 It must be remembered that the fields used to calculate the surface currents 7 s=n x i t and 7 m= Ij x n on the virtual surface S are phasor quantities. The phasor quantities are obtained from the Fourier transform of the time domain fields. One way to accomplish this is to record the appropriate field quantities on S from the FDTD simulations as functions o f time, then calculate the Fourier transform of each field quantity on S subsequent to the FDTD simulation. Another method, which is described 148 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. by Taflove,56 is to calculate a running transform during the time domain simulation by summing the quadrature components o f the time domain fields at each time step of the simulation. The sine and cosine quadrature components are calculated by multiplying the field component at time t=nAt by sin(©i nAt) and cos(co; nAt) where is the angular frequency at which the pattern point will be calculated. This second method only requires the storage o f the running sum of each quadrature component at each point on the surface S for each frequency fj. This second method was chosen to calculate the radiation patterns for the biconical antenna since it requires far less storage when a small number o f frequencies are involved. The far-fields for the biconical antenna can now be calculated using the integrals o f the surface currents over a virtual surface S in (5.11), as illustrated in Figure 5-3. To calculate the radiation pattern in the far field, the integrals in (5.12) are calculated using the fields from FDTD simulation. For the biconical antenna, the E-fields are vertically polarized, and therefore E*=0. Also, it is sufficient to calculate Ee only since H*= and He= E /q ^ O . According to (5.11), the Ee component requires the evaluation of L+ and Ne over the virtual surfaces p=p 0 and z =±Zq. 149 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Surface Current “► Jc or r- r -► Field Point Virtual Surface PEC BC Figure 5-3. Virtual surface S for determination o f far-fields. In cylindrical coordinates, the value o f Jm4 =-Jmpsin(<j)-<j)') + Jro4cos(<{>-<|>'), and the value o f L4 at p=p0 is calculated as follows: ^ = J J (E x ip) , e * ' “ *ds' Zo 2it = J J e z(Po, z ')c o s(({.-f)e jlc(p“sin0“ s(^ '> +z'“ se)pod<|)'dz' -zo 0 Za (5.13) 2It = jE z(p0,z,)cos(kz, cosG)dz' Jcos(<j>-«j)')ej1cPoSm0cos(*~*’)d<|>', 0 0 where the integration over zf runs from z=Zo to z=-Zo due to the image fields about the PEC boundary condition at z=0 and Ez(p0,-z)= E^poj-z). The bar over the field quantities indicates phasor quantities. The contribution to L4 due to the E-fields at z=Zo and their image at z=-Zq is as follows: 150 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. L* = JJ(E x az)+ejkPoCOSM/ds' Po 2 k = - J | E + ( p , Z o ) s i n ( ( J > - ^ V ^ ’^ ^ ^ ^ d ^ p ' d p ' 0 0 P b 2 it + J JE^(p,-Z(,) sin(<|> - <j»')ejk(p,siliecos(*-*,)-zieose)d^,p'dp' oo P0 2 it (5.14) - / /Ep(p,zi)cos(«j^-«|^')ejk(p'sin0cos(♦-♦'>fzScose)d(|^'p'dp, 0 0 P o2 It + / j E p(p,-Zo) cos(<J> - <|»')ejlt(p,si,,0cos(*_*')-2"“ se)d<j>,p,dp' 0 0 Po 2k = -2cos(kZg cosG) J e p(p',Zo) Jcos(<|> -<j>')ejlcp,smecos(*_^')d(j>'p,dp', o o where E*(p',Zo)= -E*(p',-z„) and Ep(p',Zo)= -Ep(p',-Zo) was used in (5.14). In cylindrical coordinates, J^JjpCosGcosOji-^+J^cosGsinOji-ifO-J^sinG and the value of Ne at the surface p=p0 is given by: N, = J f (f tx ap)0ejkR,“ svds' Zo 2 it = - J jH 2(po,z')cos9sin((t»-(t)')ej,c(p“sin0cos(^ ' Hz'cose)pod(j),dz' -z o 0 Zo 2 k - } J ( p 0, z') sin 0ejk(p»sinecos(*_*')+z,cos0)pod^ 'dz' -Zo 0 Zo = - 2 p 0 sinG j ( p o 2 it 0 , z ') cos(kz'cos 0 )dz' 0 where H^PojZ)^ -H^Po^z) and H4(p 0,z)= H+(p 0,-z) due to symmetry. The value of Ne at z=z„ and z=-Zo is given by: 151 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (5.15) N 0 = j j( a 2 x H ^ e ^ ^ d s ' 0o 2* = - J /H +(p,z^)cos0cos((|> 0 0 Po 2 * +JJ (P,-Zq) cos 0 cos(<(>- <j>')ejk(p'sinecos(*"*'>_zicos0)d<(»'p'dp' 0 0 (5-16) + J JHp(p, z0 )cos 0 sin($ - <j>,)ejk(psinecos(*_*,>+z«cos0 )d(|>,p'dp' 0 0 Po2it - J J h p(P,-Zq ) cos 0 sin(<|> - <|»' )ejfcCp'smGcosc+-4')-z; eos®>d«j»rp'dpr 0 0 Po 2k = -y2sin(kZo cos 0) cos 9 j*Hp(p',Zo) Jcos(<(> - <j»,)ejkp,s,n0 cos(^ ,)d<t>'p'dp', o o where H*(p ,Zq)= H+(p,-Zo) and Hp(p,Zo)= Hp(p,-z„) was used. The value of E0 in the far field at a frequency fj is calculated using the phasor quantities for L* andN„ given above in (5.11). The integrals in (5.13) to (5.16) for L* and N 0 are integrated numerically using the phasor fields derived from the FDTD simulation. Comparison of calculated and measured antenna patterns for a biconical antenna The calculated far-field patterns in decibels for a 2.0 in. diameter biconical antenna with a half-cone angle i|/=60° is shown in Figures 5-4 to 5-8 for frequencies of 8, 10, 12,14 and 16 GHz. The calculated directive gain of the antenna is shown as a function o f the spherical angle 0, for angles between 0° and 180° in 4° increments. The directive gain is defined by :40 152 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where U0 is the radiation intensity given by r^Eep/n,, and P,^ is the total power radiated from the antenna found by integrating ^lEep/Tjo over a sphere. Generally, the agreement between the FDTD results and the measurements was very good between angles o f 90° and 180° . The FDTD results betwen 0° and 90° are mirror images of the results between of 90° and 180° because of the assumed symmetry. The measured results were not quite symmetric due to the antenna mount which interfered with the measurements when the anntenna was rotated beyond 90°. 153 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FDTD Bicone (meas) CO 3* _e -10 a © -20 + -30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-4. Measured and calculated far field patterns for a biconical antenna at 8 GHz. FDTD ! ■ ! Measured !i i [ ! Bfi •S S 5 © - l0 0 30 60 90 120 150 180 Theta (degrees) Figure 5-5. Measured and calculated far field patterns for a biconical antenna at 10 GHz. 154 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FDTD Measured ca a O 5 *10 0 30 90 60 120 150 180 Theta (degrees) Figure 5-6. Measured and calculated far field patterns for a biconical antenna at 12 GHz. 10 FDTD Measured; 0 ca "9 . V/ eg <3 -20 -30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-7. Measured and calculated far field patterns for a biconical antenna at 14 GHz. 155 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 FDTD Measured!: s 5 W *10 -20 -30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-8. Measured and calculated far field patterns for a biconical antenna at 16 GHz. Calculation of radiation pattern for an axi-symmetric biconical antenna including a rf choke for side-lobe reduction A rf choke network was added to the biconical antenna for the purpose of reducing the radiation in the sidelobes. The choke network consisted o f a series o f onequarter wavelength slots positioned above and below the biconical antenna as illustrated in Figure 5-9. The chokes act as electrical opens that reduce the currents that flow around the top and bottom edges o f the antenna, with the objective o f reducing the radiation at angles near 0° and 180°. In a effort to reduce the sidelobes over a broad frequency range, four chokes were used where each one was tuned to a different frequency. The four chokes were tuned to frequencies of 8 ,10,12 and 14 GHz. 156 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. @ 14 GHz X/4@ 12 GHz m @ 10 GHz X/4 @ 8 GHz 0.25 in. m Symmetry Plane Figure 5-9. Choke network for reducing sidelobes of a biconical antenna. Comparison of calculated and measured antenna patterns for a biconical antenna included a rf choke A comparison of the calculated and measured antenna patterns for a biconical antenna with a choke network are shown in Figures 5-10 to 5-19. For comparison, the measured antenna patterns without the choke network are also shown in the figures at frequencies of 8,10,12,14 and 16 GHz. The agreement between FDTD calculations and the measured antenna patterns was very good for the antennas with the rf chokes. Again, the measured results were not quite symmetric due to the antenna mount which interferred with the measurements when the anntenna was rotated beyond 90°. The reduction in gain due to the rf chokes is readily seen up to 12 GHz. 157 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 0 10 -20 — FDTD w/choke ■—Bicone w/choke (m eas) j - - Bicone w/o chokeCmeas) \ •30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-10. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 8 GHz. TV ■o s S a -20 — FDTD w/choke —Bicone w/choke(meas) ! -30 4i 0 30 60 90 120 150 180 Theta (degrees) Figure 5-11. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 9 GHz. 158 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 e -10 "3 w v -20 FDTD w/choke — —Bicone w/choke (meas) — Bicone w/o choke (meas) -30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-12. Calculated and measured antenna pattern for a biconical antenna with a r f choke network at 10 GHz. 10 ;V ea ■o seg - l0 O —FDTD w/choke —Bicone w/choke (meas) 0 30 60 90 120 150 180 Theta (degrees) Figure 5-13. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 11 GHz. 159 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 S' •9 e i FDTD w/choke I — —Bicone w/choke (meas) ! — Bicone w/o choke(meas) I -30 T 0 i 30 90 60 120 150 180 Theta (degrees) Figure 5-14. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 12 GHz. -20 j 1 J -30 0 30 FDTD w/choke Bicone w/choke (meas) 60 90 120 150 180 Theta (degrees) Figure 5-15. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 13 GHz. 160 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 A v _ - n e '5 O -20 — FDTD w/choke - —Bicone w/choke (meas) - - Bicone w/o choke(meas) -30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-16. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 14 GHz. to .---------------------------;----------------------------------------- -------------- A £ ‘3 -20 FDTD w/choke j Bicone w/choke (meas) [ -30 0 30 60 90 120 150 180 Theta (degrees) Figure 5-17. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 15 GHz. 161 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 \i' -20 FDTD w/choke j — —Bicone w/choke (meas) | Bicone w/o choke(meas) | 0 30 90 60 120 150 180 Theta (degrees) Figure 5-18. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 16 GHz. 10 0 -10 -20 | FDTD w/choke I — —Bicone w/choke (meas) j ■30 0 30 — Bicone w/o choke (meas) 60 90 Theta (degrees) 120 150 180 Figure 5-19. Calculated and measured antenna pattern for a biconical antenna with a rf choke network at 17 GHz. 162 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 6 Curvilinear, Conformal and Nonorthogonal Grids Calculation of resonant frequencies for a cylindrical waveguide using cylindrical and conformal FDTD grids While investigating the accuracy of Cartesian grids for cylindrical resonators, it was observed that the error in resonant frequency was much larger than that for rectangular waveguide. The calculated and theoretical resonant frequencies for a few of the TE and TM modes for a 1.02 in. diameter by 0.612 in. long cylindrical cavity are shown in Table 6-1. The resonant frequencies were calculated using a staircased Cartesian FDTD grid with cell dimensions of Ax=Ay=Az=0.051 in. and a time step o f 2.159 picoseconds. The calculated resonances in Table 6-1 required 300 CPU secs on a Cray YMPel supercomputer. Note that the error in the calculated resonances in Table 6 l exceeded 1% in three cases and 2% in one case. This compares to the calculated resonant frequencies in Table 3-1 of Chapter 3 for rectangular waveguide where the error never exceeded 0.5%. Table 6-1. Resonant frequencies for a 1.02 in. diameter x 0.612 in. long cylindrical cavity using a staircased Cartesian FDTD grid with 20x20x12 cells, Ax=Ay=Az=0.051 in., and 32,768 time Steps. Resonant Mode TE U, TM ju TEj,! TE0U TM,n Theoretical MHz Calculated MHz Error Mhz (%) 11,797 13,103 14,826 17,106 11,732 13,259 14,503 17,019 -65 (-0.55) 156 (1.19) -323 (-2.18) -87 (-0.51) 17,106 17,386 280(1.64) 163 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The FDTD update equations in cylindrical coordinates are easily developed from Maxwell’s equations in cylindrical coordinates, as was done in Chapter 2 for the Cartesian case. The FDTD update equations in cylindrical coordinates for the axisymmetric case were shown in Chapter 5. For comparison to the Cartesian FDTD grid, the resonant frequecies calculated using a FDTD grid in three dimensional cylindrical coordinates is shown in Table 6-2. The cell dimension in the radial and axial directions were both equal to 0.051 in., and the cell dimension in the azimuthal direction was equal to 12 degrees. For the cylindrical case, a smaller time step is required for the update equations in cylindrical since the radial lines converge at the center. A stability parameter of s=0.1 was used for the cylindrical grid, resulting in a time step o f 0.4318 picoseconds. To obtain the same frequency resolution as the Cartesian grid case, four times as many time steps were required because to the small time step. The cylindrical grid required 660 CPU seconds on a Convex supercomputer. The accuracy o f the cylindrical case in Table 6-2 was better than 0.5% for the same resonant frequencies as shown in Table 6-1 for the Cartesian case. Table 6-2. Resonant frequencies for a 1.02 in. Diameter x 0.612 in. Long Cylindrical Cavity using a cylindrical FDTD Grid with 10x30x12 (r,<|>,z) cells, Ar=Az=0.051 in., A<ji=7t/15, and N=65,536 Time Steps. Resonant Mode TE,U TM0Il TEju TE0Il TMU, Theoretical MHz 11,797 13,103 14,826 17,106 17,106 Calculated MHz 11,767 13,075 14,771 17,033 17,033 Error Mhz(%) -30 (-0.25) -28 (-0.21) -55 (-0.37) -73 (-0.43) -73 (-0.43) 164 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The error in Table 6-1 for the resonant frequencies o f a cylindrical cavity was apparently due to the use o f staircased grid approximation for the curved boundary. The staircased approximation for the geometry of a cylindrical resonator is shown in Figure 6-1 for a very coarse Cartesian grid. Only the cells in black would be used during the FDTD simulation. The staircased grid conforms more to the outer conducting boundary as the number o f grid cells is increased. However, the Cartesian grid converges to the true solution for cylindrical resonators more slowly than was the case for rectangular waveguide resonators. —*\ Ax [+— Az Figure 6-1. Staircased approximation for a cylindrical resonator in a coarse Cartesian FDTD grid. Jurgens50 proposed a method to reduce the error due to the staircased approximation in a Cartesian geometry and increase the convergence rate. In his paper, modifications to the basic Cartesian grid FDTD update equations were described to improve the curved boundary approximation in two dimensional scattering problems. These modifications were based on the integral forms of Ampere’s and Faraday’s laws, as shown below: 165 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. f H d ! = | f f eE d S , C S f i - d T = - — /JfiHds, c s (6.1) where the path C encloses the surface S. Jurgens’ procedure only required modification to the update equations for fields next to the curved boundary. The remaining update equations for the Cartesian grid were unchanged from the usual update equations based on the differential form o f Maxwell’s equations shown in (2.1). The same time stepping procedure was also used. In this dissertation, the contour path method o f Jurgens50 was adapted for calculation o f the resonant frequencies in a three dimensional right circular cylinder. Figure 6-2 illustrates the contour path approach for update of the Hz-fields at the top (or bottom) of a three dimensional Yee cell. The update equation for H*1 at time n+1/2 follows from the integral form of Faraday’s Law (6.1) as shown below: HU.U2 (i + 1 j + I k) = H ^ -',2(i + I j + i k) + p j-----S-n(i + - , j + - , k ) 2 2 (6 .2 ) ^ . i i.\ B .l l / : . i ? . ^ i . \ . a i- A j i .-. • . Ax • E" (i + —, j +1, k) —nB • rEy’n(i +1, j + —, k) + A • E*’"(i, j + —, k) 2 where A and B are the lengths of the contour path along the sides of H^, and S is the area enclosed by the contour path. Note that the value of E, does not appear in (6.2) since it is tangential to the PEC boundary and is therefore equal to zero. 166 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The normal update equations would be used for Hz everywhere except where the cells are distorted by the curved boundary. For example, the normal Cartesian grid update equation would be used for H*2 in Figure 6-2, as shown below: [ E ; ( i + i j + i , k ) - E ; ( i + i j,k> E j o + u + i k ) - E ; ( i , j + i k) Ay V Ax J \ J Note that (6.3) is identical to the update equation found by applying the contour path to a normal rectangular cell. 167 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Ax curved PEC N boundary Figure 6-2. Contour path integration for updating the Hz-field next to a curved PEC boundary. Equations similiar to (6.2) would be used for updating H* when the curved boundary cuts through the right, left or top of the cell. Similiar equations also are used for updating H* and Hy in planes that are normal to the xy-plane shown in Figure 6-2. In each o f these cases, the value of the H-field that results from the contour path integration is used in the update equations for the E-fields in the adjacent (undistorted) cells. Most of the E-fields are updated using the normal update equations. Some special cases require separate consideration. In particular, the value of EyAin Figure 6-2 cannot be updated in the normal manner since a value ofH* does not exist to its left. Jurgens suggested that either a separate update equation should be specified for EyA based on Ampere’s law, or the “nearest neighbor” E-field value, for example Ey2 in 168 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 6-2, should be substituted for EyA. The resonant frequencies o f the cylindrical resonators were calculated using both approaches. The calculated resonant frequencies for a 20x20x12 cell grid using the contour path method is shown in Table 1-3. For most o f the grid, the same cell dimenstions of Ax=Ay=Az=0.051 in. were used as for the Cartesian grid. Only the cells along the cylindrical boundary were distorted as shown in Figure 6-2. The calculations required about 130 seconds on a Cray YMPel supercomputer for 16,384 time steps with a time step o f 0.1080 picoseconds. The time step for the contour path method was smaller than that used for the Cartesian grid due to the smaller cells used along the curved boundary. With the exception of the TE02I mode, the resonant frequencies calculated using the contour path method were closer to the theoretical frequencies than the Cartesian grid, and were just slightly better than the cylindrcal grid. However, a large error o f about 2.7% was observed for the TE,,, mode using the contour path approach. Apparently, the numerical dispersion for the TE,,, mode using the contour path method is not any bettter than for the Cartesian grid. The numerical error for the TEj,, mode using the cylindrical grid was only 0.4%. Table 6-3. Resonant frequencies for a 1.02 in. diameter x 0.612 in. long cylindrical cavity using the contour path method with 2 0 x2 0 x 1 2 (x,y,z) cells, Ax=Ay=Az= 0.051 in., and N=16,384 Time Steps. Resonant Mode T E„| TM0|i TEj,i TE0II TM IU Theoretical MHz 11,797 13,103 14,826 17,106 17,106 Calculated MHz 11,758 13,071 14,432 17,082 17,043 Error Mhz(%) -39 (-0.33) •32 (-0.24) -394 (-2.66) -24 (-0.14) -63 (-0.37) 169 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The contour path method did exhibit some instability during the calculation o f the E-and H- field time series. For case just presented, the instability became apparent after approximately 20,000 time steps. Figure 6-3a shows a plot of the Enfield component during the first 2500 time steps where no instability is apparent Figure 6-3b shows the Ex-field component between time steps 21,500 and 24,000 where the magnitude of the Ex-field has grown 50 times the earlier value. Due to the instability, the time domain fields continue to grow, usually resulting in a crash of the computer program due to overflow o f the allowed numerical values o f the computer. O f course, the growth o f the fields is non-physical since the source Gaussian source had long since decayed to zero. The instability is apparently due to a slight inaccurate representation o f Maxwells integral equations by the contour path formulation. Several different schemes were tried to reduce the instability, such as linear and quadratic extrapolation o f the field instead o f using the nearest neighbor. The calculation o f EyAusing a contour path version o f Ampere’s Law was tried. None of these attempts made any significant improvement in the stability o f the algorithm. 170 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.6 0.4 - 02 - •9 «S 0 Ed" - 0.2 - -0.4 - - 0.6 -■ 0 500 1000 1500 2000 2500 23500 24000 Time step n (a) 40 20 ■o I u 0 -20 -40 21500 22000 22500 23000 Time step n (b) Figure 6-3. Plot o f the Ex-field versus time during the (a) first 2500 time steps and (b) between time steps 21500 and 24000 where the instability is apparent. The solution of Maxwell’s equations in the time domain using the contour path method was not pursued any further because of the difficulty in implementing the 171 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. algorithm for general curved structures, and because o f the instability problem discussed above. The slight increase in accuracy that was observed for the contour path method in some cases can be more easily achieved using slightly larger Cartesian or cylindrical grids (i.e., with more cells). Here, the contour path method was found to be unstable for resonant cylindrical structures. However, the method has been used successfully for scattering problems as discussed by Jurgens. In many scattering problems, the time domain simulation only requires a few thousand time steps, during which time the instability may not be observed. Calculation o f resonant frequencies for a cylindrical waveguide using non-orthogonal FDTD grids It was believed that some o f the inacurracies encoutered using the contour path method were the results of a mesh distortion being applied only along the cylindrical boundary. It was reasoned that better results might be obtained if the mesh were distorted in a more global or uniform sense. That seemed to be the case for the cylindrical FDTD grid that provided excellent accuracy for cylindrical geometries. Here the interest was in developing a model than could handle more general geometries than a strictly cylindrical grid or Cartesian grid. For example, a grid was desired that could accurately model the cylindrical resonator inside of a rectangular waveguide, which was discussed in Chapter 3. Holland50 had originally discussed the generation o f grids for FDTD simulations based on generalized non-orthogonal geometries. This method was adapted to problems involving discontinuities in waveguide by Lee52 and resonant frequencies in cylindrical 172 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. cavites by Harms.53 Lee included a discussion o f stability that indicated no instabilities were found using the non-orthogonal approach as long as a small enough time step was used. A derivation of the largest time step that could be used and still maintain stability was included in his paper. A non-orthogonal system is based on a set of unitary vectors defined by69: dr _ dr _ dr a' " 3 u ‘ ’ 1‘2 ~ 3 u 2 ’ 33 ~ 3u3 ’ (6'4) and a set o f reciprocal unitary vectors defined by: 1 _ a =— a 2 x a_ 3, a = -1a_ 3 x a- „ _ 2 where (6.5) V = ar (a2 x a3) = a 2 -(a3 x at) = a 3 -(a, x a2), and (ul, u2, u3) is a set of non-orthogonal coordinate axes. The unitary vectors and the reciprocal unitary vectors are not necessarily unit vectors, but they do satisfy a reciprocal relationship a‘xaj=5ij, where 5fj is the Kronecker delta symbol. The unitary vectors are tangent to their respective coordinate axes u', and the reciprocal unitary vectors are a set o f vectors that are normal to the plane u^onstant. Any arbitrary vector 2 can be expanded either in terms o f the unitary vectors or the reciprocal unitary vectors as follows: E = X e ai = X eia*’ where (6.6) e 1 = E •a' and e: = E *a.-. 173 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The e' are the contravariant components o f the vector Ei and the e( are the co variant components o f the vector i f . Since the unitary vectors are not generally normalized, the contravariant and the covariant components may have strange units. If normal units are desired, then an alternative set o f normalized unitary vectors can be defined as (6.7) where now, in terms o f this set o f normalized unitary vectors, the vector if can be expanded as: E= + E 2 i2 + E3 i3 , (6.8) where the E{components are related to the e‘ components by Ej (6.9) = V a i ' a i e ‘ = g n e '* The gKparameter is given in general by ^ ,J 5x 5x 5u‘ 3u 5 8 y dy dz dz 5u‘ du.* 3u‘ 5uj (6.10) The x,y,z in (6.10) are the Cartesian coordinates. Stratton69 derives the vector curl operator in a general nonorthogonal system. Using the definition o f the curl operator, Holland51 and Lee52 give the following the finite difference update equations for the contravariant components of the magnetic and electric fields, h 1 and e1: 174 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. h U * /J (j j I k +I ) J 2 2 = h '- - ''2 ( iJ x + v J i jk+I ) +----------------------------------------------------------------------- — 2 2 1 1 T 1 HO,j + 2 ’k + 2^V8(i’J”+ 2 ’k + 2) e 2 (i»j + ^ , k + l ) - e 5 (i,j + ^ ,k ) e 5 (i,j + l,k + ^ ) - e 5 (i,j,k + ^ ) Au3 V Au2 I J J (6.11) e I n+1(i + j , j , k ) = e I n( i + y , j , k ) + j V ~ i 6(i + - , J , k ) y g ( i + - , j , k ) h"+,/2(i + j , j + ^ - , k ) - h “+1/2(i + j , j ~ ^ - . k ) _ (6.12) _ h 2 "+,/2(i + j , j , k + j ) - h 2 n+,/2( i + ^- , j, k _ In (6 .11) and (6 .12), -y/g =V which was defined in (6.5). The remaining contravariant magnetic field components h2 and h3 can be written by cyclic permutation of the indexes in (6.11). The same is true for the e2 and e3 components using (6 .12). Note that (6.11) and (6.12) are very similiar to the Cartesian update equations. The contravariant and covariant components o f magnetic and electric fields can be viewed as occupying the same relative positions on the basic non-orthogonal Yee cell as in a orthogonal system. For example, the magnetic field components occupy the middle o f the Yee cell faces, and the electric field components occupy the middle o f the Yee cell edges, similiar to Figure 1-1. The shape o f the Yee cell in a nonorthogonal system is skewed along the directions of the coordinate axes u1, u2, and u3, as illustrated by Holland. 175 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Note that the update o f the contravariant h-field component in (6.11) requires the covariant e-field components. Likewise, the contravariant e-field update equation in (6 .12) requires the covariant h-field components. Thus, subsequent to the update o f the contravariant magnetic field components, the covariant components o f the magnetic field must be determined to be used in the update equations for the contravariant electric field via (6.12). Likewise, the co variant components o f the electric field must be calculated before updating the contravariant magnetic fields. The transformation from contravariant components to covariant components is given by :69 h , « I g , hi ande; = 2 g i j e i- (6.13) where the gg were given by (6 . 1 0 ). Since the h‘ (e‘) occupy different faces (edges) of a Yee cell, an average of four nearest neighbors is taken to translate the contravariant components to a common location. The transform for the covariant magnetic field component hj is as follows: h i(i,j + ^ , k + ^ ) = g u (i,j + ^ , k + j ) h I(i,j + ^ ,k + ^-) b2(i + ^ , j , k + ^ ) + h2(i + ^ , j + l,k + ^ ) (6.14) + h 2 ( i - ^ , j , k + ^ ) + h 2 ( i - y , j + l,k + ^ 0 | 8i3(1»J + 2 ’ 2 h 3(i + ^ J + ^>k) + h3(i + j , j + ^-,k + l) + h 3 ( i - ^ , j + ^ ,k ) + h 3 ( i - ^ - ,j + ^ , k + l) Equations similiar to (6 .12) are used for h2, h3, and the for the covariant e-field components e l5 e2, and e3. Other than the extra steps involved in calculating the 176 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. covariant components, the update o f the nonorthogonal FDTD fields follows the same computational procedure as for the Cartesian case. A nonorthogonal grid that was used to calculate the resonant frequencies of a 1.02 in. diameter by 0.612 in. long cylindrincal resonator is shown in Fig. 6-4. For this problem, a nonorthogonal grid was only used perpendicular to the z-axis. The z-axis was therefore coincident with the u3-axis, and the contravariant and the covariant components along the z-axis were identical. Figure 6-4 shows a 20x20 cell nonorthogonal grid in the u‘u2 plane. Twelve equally spaced divisions were used along the z-axis normal to the plane o f Figure 6-4. The grid in Figure 6-4 was largely drawn by inspection. Quadratic generating functions were used to generate some part o f the axes, but some points had to be specified individually. The strategy was to develop a set o f curves for the u, and u, axes that made a smooth transition from a nearly square Cartesian grid at the center o f the coordinate system to the cylindrical boundaries. The grid in Figure 6-4 is similiar to grids illustrated by Lee52 and Harms.53 177 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0 .6 ' 0.4- 0. 2 - -0 .2 - -0.4- -0 .6' - 0.6 -0.4 - 0.2 0 0.4 0.6 Figure 6-4. Non-orthogonal grid used perpendicular to z-axis for 1.02 in. diameter cylindrical resonator. After generating the grid, the next step was to define the metrics gfJ used in the update equations (6 . 1 1 ) and (6 . 1 2 ), and in the contravariant- to covariant-component transform (6.14). The metrics were calculated using (6.10) and the x,y,z coordinates for the particular field component. For example, to calculate g ,2 in (6.10) for h, at the grid point u‘= iAu1, u2=(j+l/2)Au2, u3=(k+l/2)Au3, the following expression was used: 178 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. dx dx ® 12 dy dy 3u‘ 3u 2 + Su1 3u 2 Xc(u l -f Au*,u2) —Xc( u ‘ - A u ^ u 2) Xc(u l, u 1 4 - Au2) —Xc(u l, u 1 - Au2) 2Aul 2Au 2 (6 l5 ) YgCu1 + A u^u 2 ) - Y c( u ‘ - A u^u2) Y ^ u ^ u 1 + Au2) - YcQ^.u 1 - Au2) 2Au‘ 2Au2 The symbol X^u^u2) and Yc(ul,u2) refers to the Cartesian x- and y-coordinates at the position u',u2in the nonothogonal system. The z-derivatives in (6.15) are equal to zero since the z- (u3) axis is perpendicular to the ul and u2axes in this example. The metric g ,2 in (6.14), and all other metrics in this example are independent of the u3 coordinate. A separate value for gu and g I2 is stored for each grid location of the h, and e t components. The value o f -\/g for the update equations also has to be stored, which for this problem is given by: (6.16) The calculated resonant frequencies for an empty 1.02 in. diameter by 0.612 in. long cylindrical cavity are shown in Table 6-4 where the nonorthogonal grid of Figure 6-4 was used. The calculation involved 20x20x12 cells which required 270 CPU seconds on a Cray YMPel supercompter for 32,768 time steps with a time step of 0.8636 picoseconds. With the exception o f the TEj,, mode, the accuracy o f the nonorthogonal grid was better than the cylindrical grid o f the same size. The error in the TE^u was 1% for the nonorthogonal grid, 0.4% for the cylindrical grid, and 2.4% for the contour path grid. The nonorthogonal grid does show some instability, but it occurs after more than 100,000 time steps. As in the case of the contour path method, 179 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. numerical dispersion that is worse for the nonorthogonal case than the cylindrical case apparently limits the accuracy for higher order modes, such as the TEj,, mode. Table 6-4. Resonant frequencies for an empty l.02"dia x 0 .6 " Cylindrical Cavity using a Nonorthogonal FDTD Grid in the XY-plane (non20.f) with 2 0 x 2 0 x 1 2 cells. Resonant Mode TE(U: TMo„: TEjU: TE0„: TMm: Calculated MHz 11,967 13,300 15,103 17,157 17,229 Theoretical MHz 11,955 13,246 14,953 17,215 17,215 Error Mhz(%) 1 2 (0 . 1 ) 54 (0.41) 150(1) -58 (-0.34) 14 (0.08) A modification was made to the nonorthogonal grid to accomodate a cylindrical dielectric disk, which is shown in Figure 6-5. The waveguide was 1.02 in. in diameter by 0.600 in. long and the dielectric disk was 0.68 in. in diameter by 0.300 in. long with a dielectric constant o f er =35.74. This case was also discussed in the Harms paper53 with a different grid. The calculation here involved a grid with 20x20x12 cells which required 130 CPU seconds on a Cray YMPel supercompter forl6,384 time steps with a time step o f 0.8636 picoseconds. The resonant frequencies are shown in Table 6-5. The theoretical results in Table 6-5 came from a private communication that was listed in the Harms paper to an accuracy of three decimal places. In the calculations here, the resonant frequencies were within about 1% o f the theoretical frequencies, which was better than the results reported by Harms. Harms reported accuracies o f only between 1.5% to 4.6% for the four modes. More importantly, Harms made no mention o f any instability problems with his algorithm, although he does not mention the size o f the time step he used. 180 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. - 0 .6 - -0.6 -0.4 i -0.2 ■r - ~ 0 t 0.2 — ' i 0.4 0.6 Figure 6-5. Nonorthogonal grid for a cylindrical dielectric disk inside o f a cylindrical waveguide. Table 6-5. Resonant Frequencies for a 1.02"dia x 0.600" Cylindrical Cavity with a 0.68" dia. x 0.300" Cylindrical Dielectric (er =35.74) using a Nonorthogonal FDTD Grid in the XY-Plane with 20x20x12 cells. Resonant Mode TE0,: HE,,: HE„: TMo,: Theoretical53 MHz 3,440 4,270 4,370 4,600 Calculated MHz 3,460 4,253 4,325 4,614 Error Mhz(%) 20 (0.58) -5 (-0.4) -45 (-1.03) 14 (0.3) In the calculations reported here, the program became unstable around 16,384 time steps. The stability was extremely sensitive to grid placement, and small changes would result in a more unstable solution. Attempts to decrease grid spacing (i.e. increase the number of cells) to improve accuracy resulted in greater instability. 181 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. An attempt to adapt the nonorthogonal grid to the case o f a cylindrical dielectric disk inside o f a rectangular waveguide is shown in Figure 6 -6 . This problem was discussed in Chapter 3 in the context o f a strictly Cartesian grid. The geometry consisted of a 1 in. x 1. in x 1. in. long waveguide with a 0.757 in. diameter by 0.275 in. long cylindrical dielectric er place 0.275 in. above the bottom of the waveguide. The grid lines in Figure 6 - 6 are shown spaced one-half o f a cell apart. Unfortunately, the algorithm for this problem was unstable after only 500 time steps. It is believed that the instability is caused by discontinuities in the grid axes, which resulted in inaccurate metrics. A second nonorthogonal grid that was developed for the cylindrical dielectric inside of a rectangular waveguide is shown in Figure 6-7. In this case, the algorithm was stable out to approximately 32,768 time steps, which required 200 CPU seconds on a Convex supercomputer. It is believed that this grid is more stable because the discontinuities along the grid axes are less pronounced than the grid shown in Figure 6 7. The calculated resonant frequencies for different dielectric disk sizes are shown in Table 6 -6 . 182 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.5- 0 .4- 0.3- v ^ ^ fS S S ia a iS S k i^ n V V M n ttiin iiiu i 0.2 H °‘1l ■ ■ ■ ifiiff««£*2B B B B B B B B B B B B B B B B B B B B B M 5 5 a a *"U '<>"<"" N iB B B B aaaaSSaSSSSSSSSSSSSSSaaaB B B B iiiiiiim . iiii °1 iiimiiiiiiiiMaa a B a a a a a f l a a a a B a B B a B B a a a a M M M a " i i i i ii i i ii i i J IlliJ llllll'.IS S S S llS P B B B B B B B B B B B B B B B B B iB * " ;;" :!!!!!!!!!!!!! -0 . 1 #«f If If fI -0 . 2 - '0 z /; jm M M 0 a m a n f i l m 0 /' ^ # # # f f « « | | f l l l > //i i -0.3- ttZgiZgwZgw******' i+Zm**Z**+Z*wmmmm\ -0.5----- ' -0.5 i-------------0.4 -0.3 1---------------------------0.2 -0.1 0 11---- 1---------0.4 0.5 - -0.4- 0.1 0.2 0.3 Figure 6 -6 . Nonorthogonal grid used for cylindrical dielectric inside o f a rectangular waveguide. The theoretical values in Table 6 - 6 were taken from Liang27 where the actual dimensions o f the cavity were 1.00” x 1.00” x 0.92” L. A cylindrical coordinate grid was used inside o f the dielectric and a non-orthogonal grid was used between the dielectric and the rectangular wall of the cavity. The symmetrical axis of the dielectric was aligned with the z-axis of the cylindrical grid. For stability reasons, the cavity was made 1.0”x l.0 ” in the r-«J> plane. The algorithm become substantially more unstable when different grid spacing were used for a rectangular 1.0”x0.92” cavity. The radial dimension o f the cells making up the dielectric was chosen to represent the dielectric radius with a half-integral number of cells. The remaining cells were shaped to fit the 183 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. empty space between the dielectric and the cavity wall. A fixed cell height was used throughout the cavity and was chosen so that an integral number o f cells would fit the dielectric length as closely as possible. This resulted in a cavity height that was slightly more or less than 1 .0 ” since both the length of the dielectric and the height o f the cavity could not be represented exactly by a integral number o f cells. An instability limited the number o f time steps to about 32,768. The resonant frequencies listed in Table 6 - 6 were about equal in accuracy to the results obtained with a nominally 20x20x18 cell Cartesian grid that was discussed in Chapter 3. The Cartesian grid was very stable, and the size and shape of the dielectric could be approximated to whatever degree required by increasing the size of the grid. The nonorthogonal grid, on the other hand, was extremely sensitive to small changes and it was very difficult to find a particular grid that was stable to any degree. Attempts to increase the accuracy o f the nonorthogonal grid by increasing the numbers o f cells failed. It is believed the problem lies in discontinuities along the coordiante axes that lead to an inaccurate definition of the metrics via the approximation in (6-15). The nonorthogonal approach will probably be more stable in cases where the coordinate axes are lacking discontinuities. 184 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 6-7. Revised nonorthogonal grid used for cylindrical dielectric inside of a rectangular waveguide. Table 6 -6 . Resonant Frequencies fo ra 1.00" x 1.00" x 1.00" Rectangular Cavity with a Cylindrical Dielectric (er =38) placed 0.275" above the Y=0 plane using a cylindrical FDTD grid inside the dielectric and a non-orthogonal FDTD grid between the dielectric and the waveguide wall. Dielectric Diameter (in.) 0.654 0.689 0.757 Dielectric Thickness (in.) 0.218 0.230 0.253 No. of Cells 9x16x18 11x16x22 11x16x22 Measured Results MHz 4,382 4,153 3,777 Calculated MHz Error Mhz (%) 4,412 4,176 3,797 30 (0.68) 23 (0.55) 20 (0.53) 185 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 7 Conclusions and Discussion of Future Research The Finite Difference Time Domain (FDTD) method was successfully applied to the analysis o f a variety o f microwave devices including resonators, stripline discontinuities, filters, and antennas. In general, excellent results were obtained in comparison to measurements, analytical values, and previously published results. The FDTD calculations were in good agreement with results from mode matching and moment methods. A principal benefit o f working in the time domain is that a solution can be obtained at many frequencies over a wide frequency bandwidth with a single simulation. One o f the drawbacks of the FDTD method is the presence o f numerical dispersion. The numerical dispersion leads to velocity o f propagation which is not only dependent on the grid spacing and geometry, but it is also anisotropic. Unless very fine meshes are used, the numerical dispersion will always result in some propagation delay. In Chapter 2, the numerical phase delay for plane waves in free space was shown to be a predictable phenomenon. Thus, it may be possible to compensate for the phase delays in some instances. The reliable prediction of numerical dispersion for different geometries and the means to apply corrections are good areas for future research. In resonator problems, numerical dispersion leads to an resonant frequency that is slightly too low, perhaps by a few tenths o f a percent. However, an easy solution to this problem is to decrease the grid spacing and increase the number of cells. It was shown in Chapter 3 that good results can be obtained using a mesh of moderate size for initial work, then using a finer mesh for more accuracy. 186 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. In Chapter 4, the FDTD method with a Cartesian grid provided excellent results for stripiine discontinuity and filter problems. It was shown that the FDTD method can easily calculate the scattering parameters for stripiine discontinuities over a wide frequency bandwidth. New equivalent circuits were derived from the scattering parameters that were in disagreement with previously published values. The new equivalent circuits were used to design a variety of low pass filters in which case the measured filter response was in excellent agreement with the FDTD predictions. It was also shown that the non-uniform Cartesian grid can be used to substantially reduce the run time for a FDTD simulation without seriously affecting the accuracy. In Chapter 5, a method was described that gave excellent results for the prediction o f the radiation patterns o f an axi-symmetric biconical antenna, with and without a transverse corrugation for sidelobe reduction. In this case, the FDTD method was efficiently applied in cylindrical coordinates to calculate antenna radiation patterns using a program which could easily run on a modem desktop computer. It was demonstrated in Chapter 6 that the numerical dispersion in a Cartesian grid affected the accuracy of the resonant frequency for cylindrical cavities. Here the use o f use of curvilinear, conformal and non-orthogonal grids was discussed. A cylindrical grid gave the best results for problems that conformed to a cylindrical geometry. Conformal and nonorthogonal methods promise to reduce numerical dispersion relative to a Cartesian grid in general curvilinear problems, but they present some additional problems. Specifically, it was shown that the contour path method and the nonorthogonal grid were unstable in several cases involving microwave resonators. 187 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Additional research is needed to address the stability problem for general curvilinear and nonorthogonal coordinate systems. For antennas problems, the requirement for a good absorbing boundary condition (ABC) will become more critical as the size o f the antenna increases. In the case of the axi-symmetric biconical antennas discussed in Chapter 5, the 1st order ABC worked adequately partly because the boundary could be placed far enough away from the antenna. As the size and the complexity o f the antennas increase, it will advantageous to move the ABC closer to reduce die overall grid size. Future antenna work in three dimensions will benefit from a good ABC that can be placed close to the antenna and applied with the least amount of computer resources. The FDTD method is most efficient when advantage is taken o f its time domain capabilities. In contrast, the calculation of a very long time domain response may not be very efficient when only a few frequency points or a very narrow bandwidth response is required. In such cases, the Finite Element method or the Method of Moments will probably prove to be faster and more efficient. If the structure is very complicated, however, the FDTD method may be advantageous even for narrow bandwidth work. For future work, a good mesh generator to create the FDTD grid will be a great benefit to anyone interested in solving a variety of electromagnetic problems. A convenient mean o f viewing and checking the mesh is also needed. All of the meshing work in this dissertation was conducted manually, using a lot of “IF, THEN” statements in the algorithms. If there is one thing this dissertation could have used, it was a good mesh generator. 188 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Appendix A: A Brief Description of the FDTD Computer Program A flow diagram for the main FDTD computer program is shown in Figure A-1. This same flow diagram applies to all of the programs discussed in Chapters 3 through 6. Any minor variations were mentioned in the discussion o f the various chapters. The main FDTD program accounts for essentially all o f the computer run time. Run time associated with the post-processing steps, such as calculating impedance and Fourier transforms, amounts to only a few seconds of CPU time. At the beginning o f execution of the main program, all o f the array variables are initialized. This includes arrays for the six components o f the E- and H- fields and the material arrays. The only real difference between any o f the FDTD programs developed for this dissertation occurs in the definition of the material arrays. For example, in the case o f the dielectric resonator problems, all of the array elements are set to the applicable dielectric constant within the dielectric, and to one everywhere else. In the case o f the stripiine problems, the value of the dielectric constant is set to some large value (e.g. 1 0 12) on the surface of the metal center conductor, and to the substrate everywhere else. After the arrays are initialized, the main loop for time marching the E- and Hfields is entered. The source fields are calculated at the beginning o f the time loop, which may entail setting the E-field (or H-field) at one or more points within the FDTD grid. Next, all o f the H-fields are updated throughout the entire FDTD grid for the first half time step using the H-field update equations in (2.3). This step also includes the calculation of any boundary conditions involving the H-fields. After all of the H-fields 189 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. are updated, the E-fields are updated for the next half-time step. Again, any boundary conditions involving the E-fields are set at this time. Initialize E-,H-field and e , P , ct materials arrays Begin time loop Specify source field for time step n Update H-fields for time step n+ 1 /2 including BC Update E-field for time step n+ 1 including BC Calculate Field Integrals V(t), I(t) Next time step End time loop Finished Figure A-7-1. Flow diagram for basic FDTD algorithm. After all o f the E- and H-fields are updated for the current time step, any required field integrals are evaluated. This includes voltage, currents, and power relationships that will be written to output as a function of time. At the bottom o f the time loop, any required time-dependent field quantities or integrated fields are written to output. The time step is incremented and the loop repeats until the required number of time steps is achieved. 190 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Post-processing generally involves calculation of the Fourier transform o f the time domain quantities. Fast Fourier Transform (FFT) algorithms, which were used here, are usually available on most large computer systems. It makes sense to use code that has been optimized to run on a specific machine. These codes are often part o f a package, such as the IMSL subroutines, that are readily available from the commercial vendors. 191 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. References 1 M. Y ey, “Method o f Moments as Applied to Electromagnetic Problems,” IEEE Trans. Micro. Theory Tech., Vol MTT-33, No. 10, Oct 1985, pp. 972-980. R. Harrington, “Matrix Methods for Field Problems,” Proc. IEEE, Vol. 55, Feb. 1967, pp. 136-149. 2 3 A. Wexler, “Solution o f Waveguide Discontinuities by Modal Analysis,” IEEE Trans. Micro. Theory Tech., Vol MTT-15, No. 9, Sept 1967, pp. 508-517. 4 T. Chu, T. Itoh and Y. Shih, “Comparative Study of Mode-Matching Formulations for Microstrip Discontinuity Problems,” IEEE Trans. Micro. Theory Tech., Vol MTT-33, No. 10, Oct 1985, pp. 1018-1023. 5 P. Silvester, “A General High-Order Finite-Element Waveguide Analysis Program,” IEEE Trans. Micro. Theory Tech., Vol MTT-17, No. 4, April 1969, pp. 204-210. 6 J. Webb, G. Maile, R. Ferrari, “Finite-Element Solution o f Three-dimensional Electromagnetic Problems,” IEE Proc., Vol. 130, pt. H, No. 2, Mar 1983, pp. 153159. 7 M. Beaubien and A. Wexler, “An Accurate Finite-Difference Method for Higher Order Waveguide Modes,” IEEE Trans. Micro. Theory Tech., Vol MTT-16, No. 12, Dec 1968, pp. 1007-1017. 8 BC. S. Yee, “ Numerical Solution o f Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media,” IEEE Trans. Ant. Prop., Vol AP-14, No. 3, May, 1966, pp. 302-307. 9 A. Taflove and M. E. Brodwin, “Numerical Solution o f Steady-State Electromagnetic Scattering Problems Using the Time-Dependent Maxwell’s Equations,” IEEE Trans. Micro. Theory Tech., Vol MTT-23, No. 8 , Aug 1986, pp. 623-630. 10 A. Taflove and BC. R. Umashankar, “The Finite-Difference Time-Domain (FD-TD) Method for Electromagnetic Scattering and Interaction Problems,” Journal o f Electro. Waves andApps., Vol. 1, No. 3,1987, pp243-267. 11 A. Taflove and K. R. Umashankar, “The Finite Difference Time-Domain Method for Numerical Modeling o f Electromagnetic Wave Interactions,” Electromagnetics, No. 10, 1990, pp. 105-126. 12 A. Taflove and K. R. Umashankar, “The Finite-Difference Time-Domain Method for Numerical Modeling o f Electromagnetic Wave Interactions with Aribitrary Structures,” Pier II: Progress in Electromagnetics Research, Michael Morgan ed., Naval Postgraduate School, Monterey, CA, pp. 287-373. 13 BC. S. BCunz, “A Three-Dimensional Finite-Difference Solution o f the External Response of an Aircraft to a Complex Transient EM Environment: Part I—The Method and Its Implementation,” IEEE Trans. Eletro. Compat., Vol. EMC-20, No. 2, May 1978, pp. 328-333. 192 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 14 M. J. Pickett, A. Taflove, W. Lin, D. Katz, V. Sathiaseelan, and B. Mitaal, “Initial Results for Automated Computational Modeling o f Patient-specific Electromagnetic Hyperthermia,” IEEE Trans. Biomedical Engrg., Vol. 39,1992,. pp226-237. ls T. Deveze, “High Order FDTD Algorithm to Reduce Numerical Dispersion and Staircasing,” 10th Annual Review o f Progress in Applied Computational Electromagnetics, Monterey, CA, March 21-26,1994, pp. 61-68. 16 Mur, G, Absorbing Boundary Conditions for the Finite-Difference Approximation o f the Time-Domain Electromagnetic-Field Equations, IEEE Trans. EMC , Vol. 23, No. 4, Nov 1981, pp. 377-382. 17 X. Zhang, J. F. Fang, K. Mei, and Y. Liu, “Calculations of the Dispersive Characteristics o f Microstrips by the Time-Domain Finite Difference Method,” IEEE Trans. Micro. Theory Tech., Vol 36, No. 2, pp. 263-267, Feb 1988. 18 Z. Bi, K. Wu, C. Wu, and J. Litva, “ A Dispersive Boundary Condition for Microstrip Component Analysis using the FDTD Method,” IEEE Trans. Micro. Theory Tech., Vol. 40, No. 4, April 1992, pp. 774-777. 19 P. Zhao, J. Litva, and K. Wu, “A New Stable and Very Dispersive Boundary Condition for the FDTD Method,” 1994 IEEE MTT Symposium Digest, TU1B-4, pp. 35-38. 20 K. K. Mei and J. Fang, “Superabsorption—A method to improve absorbing boundary conditions,” IEEE Trans. Antennas and Propagation, Vol. 40., No. 9, Sept 1992, pp. 1001- 10. 21 V. Betz and R. Mittra, “Comparison and Evaluation of Boundary Conditions for the Absorption of Guided Waves in an FDTD Simulation,” IEEE M icrowave an d G uided Wave Letters, Vol. 2, No. 12, Dec 1992, pp. 499-501. 22 J. Berenger, “A Perfectly Matched Layer for the Absorption of Electromagnetic Waves,” Journal o f Computational. Physics , 114,1994, pp. 185-200. 23 D. Katz, E. Thiele and A. Taflove, “Validation and Extension to Three Dimensions of the Berenger PML Absorbing Boundary Condition for FDTD Meshes,” IEEE Microwave an d Guided Wave Letters, Vol. 4, No. 8 , Aug 1994, pp. 268-270. 24 R. C. Taber, “A parallel plate resonator technique for microwave loss measurements on superconductor,” Rev. Sci. Instrum. Vol 61., No. 8 , Aug 1990. 25 G. Roan and K. Zaki, “Calculation of Losses in a Superconductive Resonator using FDTD,” 1997 IEEE AP-S International Symposium and URSI North American Radio Science Meeting, APS Session 19, Montreal, Quebec, Canada, July 13-18, 1997. 26 P. W. Grounds, “Accurate Analysis and Computer Aided Design of Microstrip Dual Mode Resonators and Filter,” Ph.D. dissertation, University o f Maryland, College Park, MD, 1995. 193 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 27 X. Liang and K. Zaki, “Modeling o f Cylindrical Dielectric Resonators in Rectangular Waveguides and Cavities,” IEEE Tram. Micro. Theory Tech.,, Vol. 41, No. 12, Dec 1993,pp. 2174-2181. 28 G. Matthaei, L. Young, and E. Jones, “Microwave Filters, Impedance-Matching Networks, and Coupling Structures,” Artech House, Norwood, MA, 1980. 29 K. C. Gupta, R. Garg, and R. Chadha, “Computer Aided Design o f Microwave Circuits,” Artech House, Inc., Dedham MA, 1981. 30 H. M. Altschuler and A. A. Oliner, “Discontinuities in the Center Conductor of Symmetric Strip Transmission Line,” IRE Tram. Micro. Theory Tech., Vol. MTT-3, Mar 1955, pp. 134-143. 31 A. A. Oliner, “Equivalent Circuits for Discontinuities in Balanced Strip Transmission Line,” IRE Tram, on M icrowave Theory and Tech., Vol. MTT-3, March, 1955, pp. 134-143. 32 G. T. Roan and K. A. Zaki, “Determination of Equivalent Circuits for Stripiine Discontinuities using FDTD,” Progress in Electromagnetics Research Symposium, Camgidge, MA, July 7-11, 1997. 33 G. T. Roan and K. A. Zaki, “Computer-Aided Design o f Microstrip Filters by Iterated Analysis,” IEEE Tram. Micro. Theory Tech., Vol 36, No. 11, Nov 1988. pp. 14821487. 34 G.T. Roan and K. A. Zaki, “ An Iterative Method for the CAD o f Microstrip Elliptic Low Pass Filter,” Alta Frequenza, Vol. LVII, No. 7, Sept. 1988, pp. 439-442. 35 A. Zhao, A. Raisanen, and Srba Cvetkovic, “A Fast and Efficient FDTD Algorithm for the Analysis o f Planar Microstrip Discontinuities by using a Simple Source, Excitation Scheme,” IEEE Microw. and Guided Wave Letters, Vol. 5, No. 10. Oct 1995. 36 D. H. Choi and W. Hoefer, “A Graded mesh FDTD Algorithm for Eigenvalue Problems,” in Proc, European Microwave Conf., Rome, Sept. 1987, pp. 413-417. 37 P. Monk and E. Suli, “Error Estimates for Yee’s Method on Non-uniform Grids,” IEEE Tram, on Magnetics, Vol. 30, No. 5, Sept 1994, pp. 3200-3203. 38 S. A. Schelkunoff, “Advanced Antenna Theory,” John Wiley and Sons, Inc. New York, NY, 1952. 39 E. A. Wolff, “Antenna Analysis,” Artech House, Norwood, MA, 1988. 40 C. A. Balanis, “Antenna Theory,” Harper & Row, Inc., New York, NY, 1982. 41 C. Papas and R. King, “Input Impedance of Wide-Angle Conical Antennas Fed by a Coaxial Line,” Proceed o f the IRE, Vol. 37, No. 11, 1949, pp. 1269-1271. 42 P. Smith, “The Conical Dipole o f Wide Angle,” Journal o f A pplied Physics, Vol. 19, Jan 1948, pp. 11-23. 194 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 43 J. G. Maloney, G. Smith, and W. Scott, “Accurate Computation o f the Radiation from Simple Antennas using the Finite-Difference Time-Domain Method,” IEEE Trans, on Antennas and Propagation, Vol. 38, No. 7, July 1990, pp. 1059-1068. 44 M. Hussein and A Sebak, “Application of the Finite-Difference Time-Domain Method to the Analysis o f Mobile Antenna,” IEEE Trans, on Vehicular Tech., Vol. 45, No. 3, Aug. 1996, pp. 417-426. 45 E. Thiele and Allen Taflove, “FDTD Analysis of Vivaldi Flared Horn Antennas and Arrays,” IEEE Trans. Antennas a n d Propagation, Vol. 42, No. 5, May 1994, pp. 633641. 46 D. Katz, M. Piket-May, A. Taflove, and K. Umashankar, “FDTD Analysis of Electromagnetic Wave Radiation from Systems containing Horn Antennas,” IEEE Trans. Antennas and Propagation, Vol. 39, No. 8 , Aug 1991, pp. 1203-1212. 47 P. Kildal, “Artificially Soft and Hard Surfaces in Electromagnetics,” IEEE Trans. Antennas an d Propagation, Vol. 38, No. 10, Oct 1990, pp. 1537-1544. 48 P. Kildal, “Definition o f Artificially Soft and Hard Surfaces for Electromagnetic Waves,” Electron. Letters, Vol. 24, No. 3, Feb 4,1988, pp. 168-170. 49 E. Lier and P. Kildal, “Soft and Hard Horn Antennas,” IEEE Trans. Antennas and Propagation, Vol. 36, No. 8 , Aug 1988, pp. 1152-1157. 50 T. G. Jurgens, A. Taflove, K. Umashankar, and T. G. Moor, “ Finite-difference TimeDomain Modeling of Curved Surfaces,” IEEE Trans. Antennas an d Propagation, Vol. 40, No.4, April 1992, pp. 357-366. 51 R. Holland, “ Finite-Difference Solution of Maxwell’s Equations in Generalized Nonorthogonal Coordinate,” IEEE Trans. Nuclear Science, Vol. NS-30, No. 6 , Dec 1983, pp. 4589-4591. 52 J. F. Lee, R. Palandech, and R. Mittra, “Modeling Three-Dimensional Discontinuities in Waveguides using Nonorthogonal FDTD Algorithm,” IEEE Trans. Micro. Theory Tech, Vol. 40, No. 2, Feb 1992, pp. 346-352. 53 P. H. Harms, J. F. Lee, and R. Mittra, “A Study of the Nonorthogonal FDTD Method Versus the Conventional FDTD Technique for Computing Resonant Frequencies of Cylindrical Cavities,” IEEE Trans. Micro. Theory Tech, Vol. 40, No. 4, April 1992, pp. 741-746. 54 L. Zhao, L. Tong, and R. Carter, “The Influence of Boundary Conditions on Resonant Frequencies o f Cavities in 3-D FDTD Algorithm using Non-orthogonal Coordinates,” IEEE Trans, on Magnetics, Vol. 30, No. 5, Sep 1994, pp. 3570-3573. S5B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical simulation o f waves,” Math Comp., vol. 31, pp. 629-651, July 1977. 56A . Taflove, “Computational Electrodynamics, The Finite-Difference Time Domain Method,” Artech House, Inc., Norwood, MA, pp. 81-92,1995. 195 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. "Measurements provided by Robert Weinert, Northrup Grumman STC, 1350 Beulah Road, Pittsburgh PA 15235. S8G. D. Vendelin, “Limitations on Stripiine Q,” Microwave Journal, Vol. 13, May 1970, pp. 63-69. 59S. B. Cohn, “ Problems in Strip Transmission Lines,” IRE Trans. M icrowave Theory Tech., Vol MTT-3, Mar 1955, pp 119-126. ®°H. A. Wheeler, ‘transmission Line Properties o f a Stripiine between Parallel Planes,” IEEE Trans. Micro. Theory Tech., Vol MTT-26, No. II, Nov 1978, pp. 866-876. 61J. Fang and D. Xue, “Precautions in the Calculation of Impedance in FDTD Computations,” 1994 Antennas and Propagaion Int. Symp. Digest, Seatle, WA, Vol. 3., pp. 1814-1817. 62X. Zhang and K. K. Mei, “Time-Domain Finite Difference Approach for the Calculation o f Microstrip Open-Circuit End Effect,” 1988 IEEE MTT-S Digest, OF-115, pp. 363-366. “ X. Zhang and K. K. Mei, ”Time-Domain Finite Difference Approach to the Calculation o f the Frequency-Dependent Characteristics of Microstrip Discontinuities,” IEEE Trans. Micro. Theory Tech., Vol 36, No. 12, Dec 1988. pp. 1775-1787. 64 LIBRA is a trademark of Eesof Inc., Westlake Village, CA. 65 X. Liang, K. A. Zaki, and A. E. Atia, “A Rigorous Three Plane Mode-matching Technique for Characterizing Waveguide T-junctions, and its Application in Multiplexer Design,” IEEE Trans. Micro. Theory Tech., Vol 39, No. 12, Dec 1991. pp. 2138-2147. 66 G. Matthaei, L. Young, and E. Jone, “Microwave Filters, Impedance-Matching Networks, and Coupling Structures,” Artech House, Norwood, MA, 1980. 67 HFSS, or High Frequency Structure Simulator, is a Finite Element simulator that is marketed by Hewlet-Packard. 68 R. Harrington, “Time Harmonic Electromagnetic Fields,” McGraw Hill, New York, 1961. 69 J. A. Stratton, 196 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .

1/--страниц