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

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

1/--страниц