# Enhanced spatial network method with absorbing boundary conditions for the study of microwave integrated circuit discontinuities

код для вставкиСкачатьINFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter free, while others may be from any type o f computer printer. T he quality o f this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely 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 of the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6” x 9” black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. UMI A Bell & Howell Information Company 300 North Zeeb 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. ENHANCED SPATIAL NETWORK METHOD WITH ABSORBING BOUNDARY CONDITIONS FOR THE STUDY OF MICROWAVE INTEGRATED CIRCUIT DISCONTINUITIES by Dragos M. Bica Diplomat Engineer University “Politehnica” Bucharest, 1985 Master of Science University o f South Carolina, 1994 Submitted in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Department of Electrical and Computer Engineering College of Engineering University of South Carolina 1997 Major Professor Committee Member Committee Member Committee Member [airman of Examining Committee Dean of the Graduate School Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 9738212 UMI Microform 9738212 Copyright 1997, by UMI Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 300 North Zeeb Road Ann Arbor, MI 48103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. To my wife, Adina ii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Acknowledgments I would like to thank my advisor Dr. Benjamin Beker for his continuous guidance, understanding, support, and unlimited patience during my Ph.D. studies at the University of South Carolina. I would also like to thank professors George Cokkinides, Douglas Meade, and Roger Dougal for iheir support and advice during my research. My former and present laboratory colleagues at USC have also helped me finding answers to some of my questions through the discussions we always had. Finally, I would like to acknowledge the support for this research provided by the US Army Research Office under Grant DAAL03-92-G-0275. iii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Abstract A new finite difference time domain method is proposed in this dissertation. This method is formulated in terms o f the equivalent circuits similar to the model proposed in the SNM (Spatial Network Method). The new algorithm is called the Enhanced Spatial Network Method (ESN) and represents a generalization and significant improvement of the SNM. The numerical properties of ESN are examined and issues of consistency, stability and convergence are addressed. In addition to the development o f ESN, a new approach to computing the “loss” coefficients for the PML (perfectly matched layer) absorbing boundary conditions is also proposed. This technique is very useful for open boundary truncation and the method can be applied to any differencing algorithm. The numerical results show that very good attenuation o f the fields (-50dB) within the PML wall can be obtained with as little as 3 cells. To extend the versatility o f ESN, lumped active and passive circuit element models are also developed. For the first time a general method for implementing 2-port devices into finite differencing algorithms is proposed, and the methodology is demonstrated by modeling a JFET (junction field-effect transistor). The model for the JFET is the first iv Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. successful attempt to simulate such transistors in an electromagnetic field solver, considering that, to date, existing models were for bipolar junction transistors (BJTs) only. Finally, to demonstrate other uses o f ESN, analysis of planar and volumetric microstrip discontinuities was carried out. Structures of engineering interest, such as band-pass filters, were modeled. The results were compared to measured or previously published data, always showing a good level o f agreement. v Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Table of contents 1. Introduction..........................................................................................................................1 1.1 Literature R eview ......................................................................................................... I 1.2 Motivation for Research............................................................................................ 11 1.3 Proposed Research...................................................................................................... 12 2. 3-D Finite Differencing Time Domain Algorithms.......................................................15 2.1 The Spatial Network Method (SNM)....................................................................... 16 2.2 The Enhanced Spatial Network Method (ESN).......................................................21 2.3 Numerical Properties of ESN ................................................................................... 26 2.3.1 Consistency..........................................................................................................29 2.3.2 Stability................................................................................................................ 32 2.3.3 Convergence....................................................................................................... 34 2.4 Numerical Treatment o f Metals in ESN.................................................................. 36 2.4.1 Metallic structures interior to Vin......................................................................37 2.4.2 Edges, comers and walls.................................................................................... 39 3. Absorbing Boundary Conditions (ABC’s).....................................................................41 3.1 Wave Equation ABC ’s ..............................................................................................42 3.1.1 Mur ABC’s...........................................................................................................44 3.1.2 Higdon ABC’s .................................................................................................... 46 3.2 Perfect Matching Layer (PML) ABC’s ................................................................... 49 vi Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.3 Proposed New PMLs for Inhomogenous M edia.......................................................55 3.3.1 ESN-Specific Implementation............................................................................ 60 3.3.2 FD-TD Implementation.......................................................................................63 3.4 Numerical Results o f Different ABC’s for ESN ...................................................... 64 4. 5. 6. Subcell Models for Active and Passive Circuit Elements..............................................73 4.1 Introduction to subcell m odels...................................................................................74 4.2 Passive Circuit Element Models................................................................................ 77 4.2.1 Resistors.................................................................................................................77 4.2.2 Capacitors..............................................................................................................78 4.2.3 Inductors................................................................................................................80 4.2.4 Voltage sources.....................................................................................................81 4.3 One Port Active Circuit Element Models (Diodes).................................................. 82 4.4 Two Port Active Circuit Element Models.................................................................83 4.5 Numerical Results for Lumped Element Models..................................................... 88 Results of Microwave Discontinuities Simulations...................................................... 99 5.1 Microstrip step............................................................................................................. 99 5.2 Microstrip resonant stub........................................................................................... 100 5.3 3-D Via through a Ground Plane............................................................................. 101 5.4 Band pass filter.......................................................................................................... 104 Conclusions and Recommendations............................................................................. 107 6.1 Conclusions................................................................................................................107 6.2 Recommendations for Future Research................................................................... 108 Bibliography........................................................................................................................... 110 vii Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. List of figures Figure 1.1 Examples o f planar discontinuities.........................................................................3 Figure 1.2 Examples o f non-planar discontinuities................................................................ 4 Figure 2.1 Basic computation block in SNM......................................................................... 17 Figure 2.2 Equivalent circuit for the Vex(ij,k) node............................................................. 18 Figure 2.3 Computational volume and edge/side computation...........................................37 Figure 2.4 Implementation of yz-plane metallization...........................................................38 Figure 2.5 Implementation of metallic walls and edges...................................................... 40 Figure 3.1 Reflections o f one-wave ABC’s for oblique incidence from [25].................... 49 Figure 3.2 Regular/PML computational space..................................................................... 52 Figure 3.3 Implementation of PMLs for inhomogenous media...........................................53 Figure 3.4 Generalized PML problem setup..........................................................................58 Figure 3.5 Reference (a) and PML (b) problems.................................................................. 59 Figure 3.6 Mur and Higdon ABC’s reflections for open microstrip line........................... 65 Figure 3.7 Reflections from PMLs and second order Higdon ABC’s ................................67 Figure 3.8 "Loss" coefficients vs. relative permittivity....................................................... 70 Figure 3.9 Frequency domain reflection errors for the numerical PM L s...........................70 Figure 3.10 Dependence of reflection error on the w/h ratio............................................... 71 Figure 4.1 Inclusion of lumped loads along x-axis............................................................... 75 Figure 4.2 Implementation of distributed resistors............................................................... 79 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 4.3 General reperesentation o f a two port device......................................................83 Figure 4.4 JFET equivalent transistor model........................................................................ 84 Figure 4.5 Relative impedance o f the microstrip in frequency domain...............................89 Figure 4.6 Distribution of lumped resistor under microstrip............................................... 90 Figure 4.7 Reflections of a lumped resistor terminated microstrip......................................91 Figure 4.8 Lumped capacitor in series with a microstrip line.............................................. 92 Figure 4.9 Lumped capacitor between a microstrip line and ground plane........................ 93 Figure 4.10 Voltage on a lumped diode in ESN ................................................................... 94 Figure 4.11 Spice simulation of the diode circuit................................................................. 94 Figure 4.12 Schematics o f FET microwave amplifier stage................................................ 95 Figure 4.13 ESN geometry of the JFET problem................................................................. 96 Figure 4.14 Spice simulation of the JFET circuit.................................................................97 Figure 4.15 ESN simulation of the JFET circuit................................................................... 97 Figure 5.1 S-parameters o f step in width discontinuity......................................................100 Figure 5.2 S21-parameter o f a stub.........................................................................................101 Figure 5.3 Geometry of via through ground plane.............................................................. 102 Figure 5.4 S-parameters o f 0.7 mm via................................................................................ 103 Figure 5.5 S-parameters o f 1.5 mm via................................................................................ 104 Figure 5.6 Geometry of band-pass filter.............................................................................. 105 Figure 5.7 S-parameters o f microstrip filter.........................................................................106 ix Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1. Introduction l . l Literature Review During the past decade, the interest in microstrip transmission lines has substantially increased, as can be seen from the growing number of reported research activity on the subject [1-12]. The driving factors behind this trend are increasing frequencies of operation and the continuing need for more accurate design methods for microwave integrated circuits (MICs). One of the problems associated with microstrip transmission lines are their discontinuities, which arise due to the need o f interconnecting various microwave integrated circuits in order to achieve specific functionality. The irregularities in the metallization cause reflections along the propagation paths of signals, modifying the behavior o f the circuit as a whole. As a result, discontinuities have been treated as independent circuit elements and their electrical characteristics are generally described in terms of the S-parameters. There are several types o f discontinuities, which can be broadly classified into planar and non-planar types. The planar discontinuities are two-dimensional (2-D) problems, in which the thickness of the substrate and the value of the electrical permittivity are constant, such as shown in Fig. 1.1a. Although the thickness o f the dielectric substrate 1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. remains invariant, the guiding metallization layer can be varied within a plane. The most common 2-D discontinuities are: open end (Fig. I.lb), gap (Fig. 1.1c), step in width (Fig. l.ld ), T-junction (Fig. l.le), cross-junction (Fig. l.lf), band pass filter (Fig. l.lg ) and tapers (Fig. l.lh ). Unlike the planar structures, non-planar MICs represent full three-dimensional electromagnetic problems. Some examples o f general 3-D MICs discontinuities, which often occur in practice, include vias through the dielectric (Fig. 1.2a), or ground plane (Fig. 1.2b), as well as air bridges (Fig. 1.2 c ). The electrical characterization of MIC discontinuities is a fairly mature research area which has seen the development and implementation of several distinct methods of analysis. The initial efforts employed a quasi-static approach, which yielded accurate results for frequencies up to a few GHz [1,2]. The first full-wave analysis of step discontinuities has been performed using the so-called waveguide model [3]. In this model, the microstrip is approximated on both sides o f the discontinuity with a waveguide, whose effective dielectric and propagation constants are identical to those of the microstrip. The fields in the two uniform regions, away from the discontinuity, are expanded into the corresponding waveguide modes and are subsequently matched at the plane of the discontinuity. 2 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 1.1 Examples of planar discontinuities Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Dielectric 2 V ia M icrostrip 2 I \ Microstrip 1 G round p lan e Dielectric 2 M icrostrip 1 D ielectric 1 V ia M icrostrip 2 D ielectric 1 G round plane b Air b rid ge I■ Microstrip 1 M icrostrip 2 " G round p la n e il D ielectric Figure 1.2 Examples of non-planar discontinuities Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The unknown modal expansion coefficients are determined from the boundary conditions at the junction and are used to calculate the corresponding S-parameters. The waveguide model was employed in the analysis of 2-D problems (step, open end discontinuities) giving accurate results at low frequencies. Its main drawbacks, however, include the inability to take into account the radiation losses and the surface wave generation. More recent efforts into the study of discontinuities resulted in the formulation of the Spectral Domain Method [4, 5], which was employed in the analysis of planar discontinuities such as steps, in both shielded and open transmission lines. In SDM the fields and currents are Fourier transformed with respect to space variables in the plane which contains the direction of propagation. In the spectral-domain, appropriate boundary conditions are enforced to relate the currents on the metallization to the electric and magnetic fields in the plane containing the discontinuity. This yields a set of integral equations which can be solved in the spectral domain to obtain the S-parameters of the structure. SDM is an accurate numerical method, but depends heavily on the assumed expansion functions for the current distribution. The combination of SDM with the Method of Moments has been used to compute the Sparameters for 3-D problems such as vias through the ground plane and air bridges [6,22]. The initial steps in this approach require finding the dyadic Green’s function for a Hertzian dipole embedded within a multilayered medium. The next step o f the procedure is to represent the currents of the metallic structure by a weighted sum of orthonormal expansion functions with unknown coefficients. Horizontal and planar portions of the 5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. metallization are assumed to carry surface current densities. These surface currents can be represented by piecewise linear basis functions, commonly referred as “rooftops” [6]. Vertically directed currents flowing on the metal post are typically approximated by volumetric current densities. The Method o f Moments, which converts the resulting integral equations into matrix form, is employed to calculate the unknown current coefficients, and subsequently the Sparameters o f the structure. The approach described in [6] applies the integral equation method for the analysis of full 3-D problems. This method is highly accurate and incorporates physical phenomena like losses due to surface waves, radiation and non ideal dielectric materials. The Time Domain Method o f Lines is yet another approach to the analysis of MIC discontinuities [7]. However, unlike the MoM, MoL is based on the solution to wave equations cast into differential form. All derivatives, except for the one normal to the plane o f propagation, in the wave equations are replaced with central difference approximations. As a result, the EM field at any point in the discretized plane of propagation, only varies along the line perpendicular to this plane. To apply the method, the discontinuity must be placed into a resonant cavity, whose walls are far from the discontinuity. The boundary conditions, due to the metallic enclosure, are enforced by ensuring that the field lines coincide with the metallization at the edges of the discretized geometry. 6 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The boundary conditions (vanishing tangential E-field on the metal) yield an eigevalue problem, with eigenvectors denoting the modal fields. Once the eigenvectors are determined, the charge and current distributions can be obtained, from which the modal fields can be derived. The modal fields can be used to find the time response of the discontinuity to the impulse excitation. The S-parameters can finally be extracted from the Fourier transform o f the time domain field wave forms. This method was shown to offer accurate numerical results for step, open end and gap discontinuities. All the aforementioned methods are based on the frequency domain formulation of the discontinuity problem. They provide the frequency response o f the discontinuity, taking into account the boundary conditions, which are built into the Green’s functions. However, for complex, non-planar geometries, and for inhomogenous substrates, the integral equation based methods are difficult to implement and volumetric methods are often used instead. Examples o f volumetric methods are the Finite Difference Time Domain (FD-TD) technique, Transmission Line Matrix method (TLM) and Spatial Network Method (SNM). All of them have been used in the study o f microstrip discontinuities. These methods are able to handle any type of planar multilayer problem, since they are based on the differential form of Maxwell’s equations. The FD-TD method is quite mature, being first proposed by Yee [8] in 1966. It has been widely used in propagation and scattering problems. Recently, it has also been employed in the study o f a wide variety of 2-D microstrip discontinuities such as open-end, gap, T-junction and cross-junction [9]. 7 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A thorough study o f 3-D through vias, such as the one in Fig. 2b, was performed in [10]. The effects o f the rod diameter, as well as the microstrip connecting angle, were assessed. In the course o f these studies it was also shown that at high frequencies the radiated power had a significant level and it could no longer be ignored in the estimation of Sparameters. In order to extract the S-parameters from the time domain results, it was necessary to transform them to the frequency domain, using the Fourier transform. The conversion from the time domain to the frequency domain was performed under the assumption that both the physical structures as well as the numerical algorithm had small dispersion for frequencies below 20 GHz. The TLM method has been derived using a transmission line type model. The electric and magnetic fields are determined at the nodes where the transmission lines intersect. The use o f reflection coefficient matrix at each node allows for generation of the time domain algorithm. Compared to FD-TD, TLM is more computationally intensive and is susceptible to yielding artificial numerical modes. The frequency domain form has also been derived for TLM and it has been used to calculate the S-parameters of MIC interconnects such as vias and air bridges [11]. The Spatial Network Method (SNM) is derived from the differential form of Maxwell’s equations, using transmission line formalism and Bergeron’s technique[23]. The method has been used in the study of waveguides, filters, bend discontinuities and interconnects [19,20,21]. However, numerical aspects o f SNM, such as convergence, stability, dispersion, as well as open boundary truncation have not been formally addressed to date. 8 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Time domain differential equation methods can easily accommodate closed boundary conditions associated with shielded structures, but they do not have the inherent ability to simulate open boundary conditions. To overcome such limitations, several Absorbing Boundary Conditions (ABC’s) have been proposed [12,13], mostly for the FD-TD method, and more recently for TLM [14]. The best results, thus far, have been obtained from Higdon boundary conditions, which are based on the finite differencing of the one way wave equations. The ABC’s have low numerical reflectivity (<25dB), if used within appropriate limits. The errors introduced by the ABC’s can be minimized using different approaches. The Geometry Rearrangement Technique [15] places the termination planes at equal distances away from the discontinuity. Another way of minimizing the errors is to simultaneously use complimentary operators and average their results, thus reducing the error by at least one order of magnitude when compared to the regular ABC’s [16]. More recently, a new type o f boundary conditions was proposed for use in conjunction with FD-TD. They have been named Perfect Matching Layer conditions (PMLs). The idea is to introduce a set o f non-physical absorbing coefficients, which yield very low reflecting coefficients (<30dB) for problems with no more that two dielectric layers [17,18]. In order to apply the PML truncation algorithm, the computational volume, where FD-TD is employed, is surrounded with thick walls in which the waves propagating normal to their surfaces are artificially attenuated. In order to avoid reflections, the ratio between the electric losses over the dielectric constant should equal the ratio between the magnetic 9 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. losses and the magnetic permeability ( — = — ), condition which must be enforced e n locally (at every point) for homogenous media. Theoretically, for continuos fields and their derivatives, the reflections due to the introduction of PMLs should be zero. In the approximate discrete case, the method generates artificial reflections, the error depending on the maximum value of the PML absorbing coefficients. Most microwave discontinuities occur at the interface between guiding structures and circuit elements. Circuit simulators such as SPICE cannot take into account the interconnects o f MIC’s or the effects of package leads of the discrete components. These can be analyzed using full 3-D time-domain field solvers, which are capable o f including lumped circuit element models for active and passive components. Thus far, lumped circuit element models for integrated or discrete devices have been developed for FD-TD only. Initial efforts of this work concentrated on incorporating passive elements and diodes into the 2-D FD-TD algorithm [32], which were later extended to 3-D [33]. Subsequently, small signal subcell models for the BJT’s have been reported for FD-TD as well [33]. Unfortunately, no models have been proposed to date for other common microwave circuit active elements, such as FET’s, or for other types o f two-port devices. 10 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1.2 Motivation fo r Research The essential difference between the integral equation and time domain differential equation methods is that the integral equation techniques offer individual solutions for each specific problem, whereas the differential equation solvers are more general, the same algorithm being used to analyze a wide range of problems. O f all methods employed in the study of arbitrary 3-D microstrip structures and their discontinuities, time domain methods appear to be most suitable for computational engines in the development o f CAD tools for microwave circuits. Due to the fast pace o f microprocessor evolution, it is now feasible to think of industrial strength CAD tools, which can be effectively used in the distributed desktop computer environment. Until now, most time domain codes required main-frame or high performance super computing resources, limiting their availability to very few end-users. Another problem frequently associated with general-purpose EM solvers is the abstract nature of the governing equations on which they are based. In general, circuit analysis is perceived by most electrical and computer engineers as rather easy to grasp, compared to Maxwell’s equations and their solutions. Therefore, a way o f recasting Maxwell’s equations in terms o f circuit elements as well as transmission lines is highly desirable to provide a more intuitive interpretation to Maxwell’s equations. Such circuit models should naturally provide subsequent interconnection of wave guiding structures, such as microstrips, to the discrete active or passive circuit elements. Given the above considerations, it has been shown that time domain algorithms can be successfully 11 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. employed in modeling waveguiding structures and lumped elements simultaneously [32,33]. The FD-TD and TLM methods have been studied thoroughly in the past, and their numerical properties are well known. FD-TD is directly derived from Maxwell’s equations using second order approximations to time and spatial derivatives, which are responsible for sub-optimal modeling of step excitations. TLM uses the reflection coefficient arrays at each node and is rather difficult to use with inhomogenous materials. In addition, it is susceptible to yielding spurious modes, needing extra care in developing and using the algorithm. On the other hand, SNM can be derived from the transmission line formalism combined with lumped circuit elements, providing an intuitive representation to Maxwell’s equations. To date, its numerical attributes have not been discussed. Moreover, the combination of SNM and ABC’s for open boundary modeling has not been addressed at all. These arguments provide the motivation and starting point for the development of a new SNM algorithm, which will significantly improve performance in terms of memory and run-time requirements. 1.3 Proposed Research From the physical point of view, any algorithm describing the propagation of electromagnetic waves should be based on Maxwell’s equations. The intent of the research in this dissertation is to show that proposed new SNM (or ESN-enhanced spatial 12 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. network method) scheme is both consistent with Maxwell’s equations and stable, and therefore convergent. Such level o f detail is required because ESN can be formulated using transmission lines and lumped elements, without resorting to Maxwell’s equations directly. In its present form, the original SNM algorithm assigns a set of six variables to each component of the electric and magnetic fields, at two consecutive time steps. Overall, the algorithm needs 72 field related variables per unit cell, which represents a very heavy demand on computer memory. In addition, existing implementations o f SNM involve a very large number of operations per unit cell at every time step (144 additions, 18 multiplications and 6 divisions). Reducing memory requirements and increasing the computational speed can significantly enhance SNM’s performance. This can be accomplished by recasting the original equations in a simpler form, with a smaller number of variables and less operations per unit cell. A method, which makes such simplifications to the algorithm possible, involves expressing the currents in terms o f voltages at every node, and retaining their time history for several steps. Until now, the research into improving the ABC’s for use as open boundary conditions has been focused on the FD-TD method only, without serious attempts of generalizing them for larger classes of finite-difference schemes. The behavior of the classical boundary conditions (Mur, Higdon) in conjunction with the ESN method are analyzed 13 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. and a new general method o f deriving the PML’s absorbing coefficients for any finitedifference algorithm is proposed. Since ESN is a general-purpose algorithm, its computational power, and range of applicability, can be increased by incorporating the ability to model lumped circuit elements. To that end, passive and active subcell circuit element models are included into ESN. This improved version of the algorithm is employed to simulate the EM response of microwave IC’s, containing both guiding elements and lumped discrete components. In order to validate the newly developed ESN algorithm, a set of test problems are analyzed in detail. Comparison with previously published data for the S-parameters of microstrip planar discontinuities such as open end, step in width and stubs are used for validation. In addition, other arbitrary 3-D structures, such as vias through the ground planes and air bridges are examined as well, to demonstrate the capabilities of ESN. 14 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2. 3-D Finite Differencing Time Domain Algorithms The differential equation methods in the time domain have emerged as numerical simulation tools, and are now becoming synonymous with the usage o f computers in solving engineering problems. The Finite-Difference Time-Domain method was first proposed by Yee [8] in 1966. Since then, it has continuously been improved by other researchers, through the incorporation of dielectric losses, open boundary conditions and lumped circuit elements. The method has been widely used in all types o f electromagnetic simulations, truly reaching its maturity. Later other time-domain methods, such as Transmission Line Matrix Method and Spatial Network Method have been proposed, especially for the analysis o f microwave discontinuities. All finite-difference methods in electromagnetics are based on Maxwell’s curl equations, which are stated below: (2 . 1) where: E is the electric field vector, H - the magnetic field vector, s - the electric 15 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. permittivity, p - the magnetic permeability, and ae, crm represent the electric and magnetic losses, respectively. In the following paragraphs the derivation o f the proposed Enhanced Spatial Network Method (ESN) is presented and its numerical properties analyzed. However, prior to considering the details of ESN, a review of the SNM is given for completeness. 2.1 The Spatial Network Method (SNM) Yoshida, Fukai and Fukukoa [19,20] proposed the original SNM in 1979 as an alternative finite-differencing algorithm to FD-TD. The algorithm is generated using circuit elements involving transmission lines, resistors and capacitors, connected in a 3-D mesh, which is build in accordance with Maxwell’s equations. A computational block of SNM, formed of 8 cells, is shown in Fig. 2.1. The solid line denotes the basic lattice cell identified by the indices (ij,k), whereas the dotted lines are used to depict the surrounding cells. All cell nodes are identified as well. They can be classified into two types: electric nodes, corresponding to the electric field components, and magnetic nodes corresponding to the magnetic field components. The computation of the fields is interleaved, since the tangential components of the electric and/or magnetic fields along the Cartesian directions are always in the same discretization plane. 16 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. X A Z ST V ^ ( i + l , j + l ,k + l ) V^(i+l«j/k+l) W ' V „ (i,j,k + 1 ) V ml(i,j.k + 1 ). X - - - < -V „(i.i+ l,k + l) V „ r(i.j.k )----- V m4(i,j+ l,k ) V„(M,k) Figure 2.1 Basic computation block in SNM Table 2.1 shows the correspondence between the equivalent circuit variables and field quantities at every electric and magnetic node. Although there are only 6 field components for each cell, the equivalent SNM circuit introduces a set of additional 6 currents for each computational node, increasing the number of cell variables to 36 [19]. Electric Nodes V™ vm y a < it 1 m N v„ Magnetic Nodes Equivalencies V« = EX le y - H , Icz = Hv V„ = Ey Icz = -Hx I« = HZ vm z Iez = -Hx Equivalencies Vmx = -Hx Imy = -Ez Imz = Ev Vmy = Hv Imz = Ex Imx = -Ez Vmz = Hz IAmx = E IAmz = -E Table 2.1 Correspondence between field and circuit variables 17 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The actual circuit elements used to interconnect the nodes o f the mesh are shown in Fig. 2.2. The gyrators are introduced in order to be able to use the same circuit nodal equations for both the electric and magnetic variables. The wave propagation characteristics of a one-dimensional transmission line are introduced into SNM by using Bergeron’s method [23]. In the following equations the dot symbol is used to denote the product between scalars. v(k,t) + z • i(Jk, t ) = v(k —1,t —At) + z ■i(k - 1 ,t - At) (2.3a) v{k - l , f ) - z - i(k - 1 , 0 = v(k,t - At) - z •i{k,t - At) (2.3b) • .z \ ^ Vmz(i,|-l/2,k)_______ • ^ .'f' V I r '_ R . ^ ^ ,Ic ~ C / ^ v ^ ^ ^ Vmi(i,i+l/2,k) Posihve gyrator negative gyrator ---------------- Figure 2.2 Equivalent circuit for the VtI(i,j,k) node 18 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Application of (2.3) to the transmission lines surrounding the VCT node (Fig. 2.2) yields the following equations: =V TF'<nxl;2’ v« (/, j , k ) ± Z Q■C 1;2(/, j , k) = i+i (/, j , k) = - r CrU( i , j T l / 2 , k ) ± ) = VF' v«‘(/, j , k ) ± Z Q■i™ T o « 3 ;4 ’ ± Z „ - v ^ ( i ,y T l / 2 ,t ) (2.4a) (2.4b) where Z0 = ^/p0 / s0 is the characteristic impedance of the transmission lines. In addition to the currents entering the va (/,_/, k) node via the transmission lines, there are also currents through the resistor and capacitor connected to the node. The final nodal equation for va (i,j,k) is obtained by enforcing Kirchoffs current law, which yields: n , ex _ W JM -Q V L + vC 4)(/,y ,£ ) + z 0 - ^ ( h j , k ) ^ n \ /ii /• J•, k/ w ’ Z 0 + R a ( i , j , k ) - ( 4 + Z7 0. - 4^ - g a (it )) where: 4 • A/ • cn Ax ’ ft (2.5b) (2.5c) 2 - Ax with c0 - the speed of light in free space, A/ - the time step and Ax - the lattice step. The T,'„k( i , j , k ) terms, given by the right hand side of (2.4) account for the waves propagating towards the va (i,j,k) node on the four transmission lines surrounding it. % is associated with the flow of current through the capacitor and is given by: 0\ j , k) = v'a (/, j , k) + Ra (/, j , k) •i'cex(/, j , k). 19 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ( 2 .6 ) The resulting explicit finite-differencing algorithm for the va ( i, j , k ) variable is derived from (2.5) and (2.6), and the corresponding algorithm can be written as (2.7). The notations for the voltages at the nodes are self-explanatory. In the case o f the currents the notations show the type o f current (e-electric vs. m-magnetic), and whether it enters or leaves the node (see Fig. 2.2). x (2.7a) x {4+g c J ‘J'k) +ga (ij.kj)~l (2.7b) (2.7c) (2.7d) (2.7e) SNM has very large memory requirements and is very computationally intensive. For every of the six nodes comprising a cell, one voltage and 6 currents, at two consecutive time steps, are needed to update the field at each iteration. As can be seen from (2.7) the algorithm computes not only the nodal voltages, but the nodal currents as well, which imposes a large burden on the computer resources. 20 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2.2 The Enhanced Spatial Network Method (ESN) One goal of the work in this dissertation was to find a finite-difference algorithm, which is based on circuit elements, like SNM, but is closer to the computational efficiency of FD-TD. A natural approach to deriving an improved algorithm was to start with the SNM’s basic circuit cell and to reduce SNM’s nodal equations to a more compact form. Reduction of the SNM algorithm begins with modifications to equations (2.4) by recasting them in terms o f normalized voltages and currents ( Vmq = v K„ = v", ■ l m = L„ ■y f e • , where 4 = x , y , z ) . To make an and explicit distinction between the electric and magnetic variables, the (/, j , k ) terms are replaced with voltages and currents using equation (2.4) and the normalizations stated above. This leads to: i k) - ^mTUTx + ^mTOTx & « 0 J> k y v n ) j , k ) + v ea{ i j , k ) a 4 + gca( i , j \ k ) + ga ( i , j , k ) where the following notation has been used: K w n = C - 0 J + 1 / 2 ,* ) - v : z0 , j - 1 / 2 , k ) ~ (2.9) - V ^ i j , k + \ I 2 ) + V'yi{ i j , k - \ I 2 ) ImTOTx = ^mrt 0 ’7’£ + 1/ 2) + I ^ 0 ,7, k —1/ 2) (2.10) - C . 0 , 7 + 1 /2 , k ) - r my2U, 7 - 1 / 2 , k) 21 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. S c * = z o / Rc(.i’J ’k ) 311(1 = Z0 ! R(i,j,k), (2 -il) and where Icex is the current flowing through the capacitor (Fig. 2.2). To simplify the algorithm for the same electric node that was used to review the SNM, Kx(iJ,k), the currents denoted by I'mJVTx in equation (2.8) must be expressed in terms of electric and magnetic voltages only. Each of the currents comprising VmTOTx is derived from nodal equations written for the magnetic nodes surrounding the electric node Vex(i,j,k): 4 I:2( i ± 1/ 2) = + C ( / , M ± 1/ 2) ± / £ , ( / , / , * ) + V ^ \ i , j , k ) , C , ;2(/,y ± 1/ 2,k) = ± 1/ 2,fc) + (/,/, k) - V ^ \ i , j , k ) , (2.12a) (2.12b) allowing for I'mTOTx to be written as a function o f VexOJ>k) and the surrounding magnetic voltages: K y O J , k + 1 / 2) + 1 / 2) + | /' mTOTx U K i i J + 1/2, k ) ~ V L (U -1 /2 ,* ) (2.13) + K y i (i’j ’k) - I'~2 K J , k ) Application o f KirchofFs current law in (2.13) for the currents at time t-1, the present time step being t, and the usage of the result in (2.8), leads to the finite difference expression for V^fiJ.k) shown in equation (2.14). 22 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The abovementioned expression still contains a restrictive artifact from SNM, namely, the normalized algorithm velocity u = Ax / (ca •A/) = 4.0. To generalize this algorithm for any normalized velocity, the “4.0”s are replaced with u and gcex(i,j,k) is redefined as gca( i , j , £) = 2 •u •(er0',y,£) - 1 ) . As a result, the restriction on v can be removed and the stability o f the algorithm can be analyzed for a whole range of new values of the normalized velocity. (4 - a „ ( / ,/ , k)) ■V'~x(/,/, k) + gca (i , j , k ) • (i j , k) + - V ^ ( i J - 1 / 2 , k) + V ^ { i J +1 / 2 ,k) + 2 + - K y U J , k - 1/ 2) - V ^ { i, j ,k + 1/ 2 \ (2.14) x (4 + « « 0, j ,k) + g cer (/, j,k))~x. Further reduction in the number of variables per unit cell can be achieved by replacing the remaining capacitor currents in equation (2.14) with expressions containing electric voltages only. Applying a central finite difference approximation to the v-t characteristic of the capacitor i = C - d V / d t , allows for Icex to be recursively defined as shown below: C U , j \ k ) = gcex( i , j , k ) - { v ^ i , j , k ) - V ^ ( i , j , k ) ) - C ( i , j , k ) . 23 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.15) The resulting algorithm, which is obtained by repeatedly using (2.15) in the generalized form o f (2.14), is shown below: ((u + gca O’,y ,k))(\ - a „ (/,y, k))) • x + + 2- (2.16) {v:y( i , j , k - \ / 2 ) - V ; y ( i J , k + \ / 2 ) ) x [(u + ^ c„ o ‘,y ^ ))0 + a « (* .y .* ))]" l» where V ^ m(i,j, k) = (/,y, A:) - V £ m( i , j , k ) . In equation (2.16) the number of variables is the same as that in equation (2.14), but the number o f operations per unit cell decreases, as the electric and magnetic currents are no longer computed. For most practical purposes, it can be considered that all materials have the magnetic permeability equal to that o f air (p0), and in this case the number of necessary variables is only 9. In order to prepare for the subsequent development of the Perfectly Matched Layer ABC’s in ESN, the following approximations will be made: (2.17a) (2.17b) = p£X(/,y,fc). 24 Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission. The above approximations are valid for small values o f a a ( i , j , k ) , which correspond to many common materials with small conductivity values. When the terms (2.17) are placed in (2.16) the “exponential” form of ESN (in which the losses are implemented as exponentials), which will only be used for the PML layers, is obtained: P e , 0 \ / , * ) - C ( /,/,* ) + 2 X ?«(U.*)-(C(i.M>+ 2-C i.) + C \ U , k ) = $a { i J , k y - V ^ ( i , j - l l 2 , k ) + V ^ ( i , j +1 / 2,k) +') (2.18) - M 2 ) - v : , ( i , j , k + 1/ 2) v+ -I With the generalizations mentioned above, equations (2.17) and (2.18) now represent a new finite difference scheme, for updating Vex (or Ex) at regular ENS lattice nodes and within the PML layer. This new ESN algorithm is equivalent to the SNM scheme, which uses central differences for the space and time derivatives. The proposed scheme is no longer restricted to a predefined normalized velocity, which was removed by redefining terms related to the dielectric (i.e., gcex). It should be added that analogous expressions can also be derived for the remaining components o f the EM fields and used together with (2.18) to solve the field problem. 25 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Comparison of the standard SNM algorithm and the proposed ESN reveals that SNM needs one voltage and five associated currents for each o f the six nodes in a unit cell, adding up to a total of 36 variables per unit cell. ESN reduces the number of variables per unit cell from 36 to 12, utilizing only one voltage and one recursive summation for each node. In addition to saving computer memory, this approach also reduces the number of computer operations per unit cell from 6 assignments, 12 additions, 4 multiplications and one division to 2 assignments, 7 additions, 4 multiplications and one division. 2.3 Numerical Properties o f ESN When a numerical algorithm is used to simulate physical phenomena, the numerical solution should approximate reality as close as possible. In order to characterize various algorithms for the solution of systems of hyperbolic equations, the notions of consistency, stability and convergence can be used. Since such numerical attributes are generally discussed in terms of finite difference operators, ESN must be expressed in a matrix form, which contains the information about the fields, boundary conditions and the excitation. In order to obtain the aforementioned form of the governing equations of ESN, the Ve and Vm nodal variables will be merged into a single vector {£}, which will allow for equation (2.18) to be written as: {£}' = [S ,]{£('-' +[£:]{£r: + [B ]2 (-1 )* * '{£}-*♦' + (S } ', *=3 26 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.19) where the arrays [5,], [5 J and [5] can readily be determined from (2.18). If the lattice boundaries are closed, that is they are composed o f Perfect Electric Conductors (PEC) or Perfect Magnetic Conductors (PMC), then the tangential electric field components or tangential magnetic fields, respectively, are set to zero. In equation (2.19) these conditions are incorporated by assigning the rows of the [5] arrays, corresponding the boundary field components, to zero. In the case o f open lattice boundary, the rows o f the [B] arrays contributing to the computation of the field components on the boundary will be assigned according to rules o f the ABC algorithm used. The vector {S}1represents the excitation which is applied to the algorithm at the Vh time step. For simplicity, equation (2.19) may be written in a more compact form: { £ > ' = £ [ c ; k s }‘ . (2.20) *= 1 The entries o f the [C] arrays of equation (2.20) can be computed using the recursive relations in (2.19) for the time steps t=I,2,..t-I, and then, by applying the induction method for t: [C;] = [51][C r 1] + [ ^ ] [ C r 2] + [B ]X (-l)* +1C1,-*+1, (2.21a) k =3 [ c ; ] = [ c r * +i], [ c 2]= [5 ,] , [ c ;] = [ /] . (2 .2 2 b) To properly evaluate the numerical attributes of ESN, some general concepts of numerical analysis will be reviewed first, based on the concepts presented in [43]. To that 27 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. end, consider a system o f partial differential equations, defined on the domain D, which, symbolically, can be represented as: £ (“) = where P e D . (2.22a) In the above equations, all terms involving the dependent variable u are included in the left-hand side, with all terms containing independent variables {t , x , y , z grouped on the right hand side, representing the forcing function f The set of points to which the boundary and/or initial conditions are prescribed is denoted by C. Mathematically, the conditions that are to be satisfied by u on C can be represented as: B(u) = g(P), where P e C . (2.22b) Note that B is not a differential operator. It represents the operator notation for the conditions imposed at various parts o f the boundary C. In order to analyze the properties of a discrete solution, versus its continuous counterpart, a physical problem represented by (2.22) must have a unique and smooth solution u for any data in some class of smooth functions {fgj. Smooth, in this case, means ‘sufficiently differentiable’, i.e. a function having N non-zero derivatives, where N is the order of the system of differential equations. The complementary discrete problem is defined on a lattice for the independent variables which appear in (2.22) with the spacing {At, Ax, A y The points of the lattice interior to D will be denoted by the set D& and similarly, boundary lattice points will be denoted by 28 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Cj. At the points o f DAu CA a difference approximation U is defined as the solution of some set of difference equations. These may be written symbolically as: (2.23a) B±(U) = g(Pb), w hereP e D A andPb e C A. (2.23b) From the numerical standpoint, it is desired that the difference solution U of (2.23) be a close approximation to (2.22) at corresponding points of DAu CA ,when it is sufficiently smooth. Furthermore, the difference solution U should be unique solution of (2.23) and its numerical evaluation should be possible despite the loss of accuracy due to the round off errors. With the above definitions of the continuous and discrete problems, the notions of consistency, stability and convergence can be now discussed. 2.3.1 Consistency Let (j)(t,x,y„.) be any function with sufficiently many continuous partial derivatives in D lX2. For each such function and every point P e D A, let the difference between the continuous and discrete problems be: (2.24a) t M P)) = L M P ))-L M P )), and for each point P e C A, let the difference between the corresponding boundary conditions be: 29 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.24b) Then the difference problem (2.23) is consistent with the continuous problem (2.22) if (2.24c) when At—^O, Ax—>0, Ay—>0, in some manner. Above, | || represent convenient norms in the sets D j and C& with r{<f>} and f3{<p) being the local truncation errors. If a given relationship between the At, Ax, Ay, that is imposed by the algorithm, has to be maintained in order to satisfy (2.24c), then the difference formulation is conditionally consistent with the continuous problem. In order to show the above defined consistency, in the context of ESN with Maxwell’s equations, (2.16) will be used. If the time differences are separated from the space differences, then (2.16) can be written as: (2.25) / 2 a a (i, j , k ) - ( u + gca (/, j , k)) ■V'a (/, j , k ) + f - V L Q d - 1 >2 , k) + + 2 +1 / 2,kj) - {v;iy( i , j , k - l / 2 ) - V : y ( i , j , k + U 2 ) Next, the substitution o f gcex0J,k) with 2u(sr - 1) in equation (2.25) yields: 30 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (1 + 2(e„ (i , j , k ) -1)) - ( C ‘ « ,/,* ) - C V .M )) / 2 + t-2 + {£rx(i ’j ’k ) - l ) / C0 X k=2i=0 I At1 = 2A / v + C * _20’, M ) (2.26) = (<*«O '.M ) / Ax) •(w + £cer( / J , &)) • ( / , / , A:) + C O J - 1 / 2,*) + C ( i , y +1 / 2 , ^ / Ax + V ' Q J * k - l / 2 ) - V ' ( i J , k + l/2) In the limit At-*0, the sum in the second term o f the left side in equation (2.26), becomes equal to (/, j , k) / a . By applying the limits as A/ —>0 and Ax -» 0 to every term in equation (2.26), leads to: s„{ij,k) ^ |. j k ( _ {u + g ca( i , j , k ) ) - a a ( i , j , k ) y j. . k ( a Ax (2.27) When the normalization, described in (2.2, pp. 22), is taken into account, along with the fact that ( v + g cex(/, j, k) •<*„(/', j, k) / Ax = yjju^ /e0crx (/, j , k ), equation (2.27) becomes identical to the x-component o f Maxwell’s magnetic curl equation, thus showing the consistency o f ESN with Maxwell’s equations. It is now clear that the values for gcex(i,j,k) and ctex(ij, k) must be chosen such that the finite difference algorithm is consistent with Maxwell’s equations. 31 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2.3.2 Stability A difference scheme that is defined by linear difference operators LA and BA is stable, if there exists a finite positive quantity k, which is independent of the lattice spacing, such that: (2.28) for all functions U defined on DAu C A. If (2.28) holds for any spacing between the lattice nodes, then the linear difference scheme { B J is unconditionally stable. If (2.28) holds for some restricted range o f spacings between the lattice nodes, in which At, Ax, Ay, may all be made arbitrarily small, then { LK BA} is conditionally stable. For the ESN, the stability of the algorithm will be demonstrated using equation (2.19) as a starting point. The explicit forms of the finite difference operators { Lx Z?A} are: L ,(E ) = {£ } '-[B,]{E}'-1~[B2]{E}'~2 - [S ] ]T (- l)* +,{£}'-*+1 = {5}', (2.29a) Ba (E) = 0. (2.29b) The use of the above expression in (2.28), shows that for ESN to be stable the following must be true: (2.30) 32 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. This result can be rewritten in terms of ESN specific matrix operators with the help of equation (2.20): (2.31) £ r c r ‘ *,]{s>‘ *=l To somewhat simplify the discussion, {5} will be decomposed using a set of orthonormal vectors {«,}={ 1,0,0,..}T,{w2}={0,1,0,0..}T. In this case, the excitation can be written as: (2.32) The orthogonal decomposition of the excitation allows for finding individual stability criteria for each column of the [C] arrays. The norms for each component of excitation and the C vectors (columns of the [C] array) will be chosen as: i IK o ( m* } |= Ik o l H k l = (2.33) 0 |C|'"*+11 = max|C,'-i+l (k, /)], (2.34) Finally, the stability criteria is met if ||C,'~ [| < 1 and £=1, as is implied by the following result: (2.35) k =I k=\ 33 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The condition ||c ,' *+1|| < 1 was chosen, as the [C] arrays are generated from repetitive multiplications of the [B] arrays, as can be seen from equation (2.21). Due to the complexities involved in computing the [C] arrays, no attempts were made to compute them symbolically. However, for the u > 4 case, which is commonly used in ESN and SNM, it was verified that all the terms of the [C] arrays are smaller than 1, thus meeting the above stability criterion. Finally, the convergence, which is a direct consequence o f consistency and stability, will be derived next. 2.3.3 Convergence Let u be the solution to equation (2.22) and let U be the solution to (2.23). The approximate difference solution converges to the exact solution means: ||n(/>)-£ /(/> )!-> 0 (2.36) for all P e D A[jCAand t—>0, Ax—>0, Ay—>0, in some manner, with || || representing a norm in the set Da u C j. Although equation (2.36) is a formal, general definition o f convergence, a slightly different approach was used to assess the convergence o f ESN. In particular, the 34 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. convergence is examined for difference algorithms meeting the consistency and stability requirements. To that end the following approach is proposed: THEOREM (see [43], pp. 521) Let {LA, BA} be linear difference operators, defined over the discrete lattice (At, Ax, Ay,..), which are stable and consistent with {L, B}. The approximate difference solution U o f (2.23) will converge to solution u o f (2.22), when (At, Ax, Ay,..) are made arbitrarily small. PROOF For each point P e D A in the lattice, the difference between (2.22a) and (2.23a) is given by: 0 = La(U(P)) - L(u(P)) = La (U (P )) - La( u(P)) + La ( u(P)) - L (u(P)). (2.37) Since the difference operators are assumed to be linear, the definition o f the local truncation error (2.24), allows for casting the problem as: La ( U - u) = x( u( P ) ) , P z Da , (2.38a) Ba(U - u) = p(w(R)), P e C A. (2.38b) Furthermore, since the linear operators in (2.38) are stable, the difference between U and u can be written as: 35 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.39) Finally, since the operators satisfy the consistency condition, letting t —>0, Ax—>0, Ay->0,... in such manner that ||t| -» 0,||p|| -» 0 in the limit, leads to ||C/ - wj| -+ 0, which demonstrates the convergence. Since it has already been shown that ESN is consistent and stable, the above theorem indicates that it converges to the continuous solution as well. The proof that ESN converges toward the continuous solution is very important, because it assures that the equivalent circuit implementation of numerical algorithms such as SNM and ESN to solve Maxwell’s equations is valid. 2.4 Numerical Treatment o f Metals in ESN The regular ESN algorithm, which for Vex is summarized in equation (2.18), is not capable o f taking into account the metal that may be present interior to the ESN computational volume, as well as on its surfaces. The following discussion outlines the special algorithms for the modeling of metallic structures at interior lattice points (such as inside Vin shown in Fig. 2.3). 36 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I V ^ V to \ z E,(i,l,2) ■ ■ E„(i,2,2) ■------■------------ - E,(i,l,l) Ex(i,2,l) Figure 2.3 Computational volume and edge/side computation It also focuses on modeling o f metals on the surface which bounds the computational domain. In the following, it is assumed that the metals are perfect electric conductors, or (PECs). The closed surface £ v, defining the boundaries of Vm, is comprised of walls, edges and comers. ESN must be supplemented with appropriate special algorithms for the computation of the fields on such geometric discontinuities. 2.4.1 Metallic structures interior to \Z„ Metallic structures can be composed of linear (wires), planar (sheets) and volumetric shapes. If the metal is perfectly conducting, then the tangential components o f the electric 37 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. field must be zero on its surface. This means that for x, y or z-directional wires, E* Ey and E~ must be set to zero. In the case of ESN, the values for VeXl Vgy and V^, respectively, are set to zero at the beginning of the computations, and their values are no longer updated. The figure below clearly illustrates some limitations of all differential equation based methods. Specifically, since field components are defined at lattice nodes and material parameters are specified over rectangular lattice cells, ESN is best suited for modeling rectangular geometries. If an object with curved boundaries must be analyzed, then a stair-case approximation is used to model its contour. Az Ay metallization Figure 2.4 Implementation of yz-plane metallization Solid metallic objects are modeled by ESN in a similar manner as described above. Their surfaces are treated identically to the planar objects, while the fields inside them are set to 38 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. zero. The treatment o f edges and comers must be done differently, and will be described in the next section. 2.4.2 Edges, corners and walls Unlike the edges of isolated planar metallic conductors, edges defined by a junction of two planar metallic surfaces, such as two plates at right angles (see Fig. 2.5), must be treated differently. Depending on the orientation of the edge, only one component of the electric field is assigned to it, and it must be orientated along the edge, such as £ .o f zx plane in Fig. 2.5. The corresponding electric field component on the adjacent wall (yz plane) does not need to be defined, since the edge was already accounted for by setting E. to zero in the zx plane. The treatment of comer cells is not different from that o f the edge cells. As can be seen in Fig. 2.5 (next page), there is no actual electric node placed in a comer, and interpolation formulas, such as those employed in the computation of the potentials at the comers in electrostatics [41], are not needed. The treatment of comers and edges discussed above can also be employed for lattice truncation, if the EM problem under consideration is enclosed within a shielded box. For 39 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. this case, the approach to enforcing boundary conditions at edges and comers is same as outlined above, and the field tangential to the walls can be set to zero directly. Legend 0 <J) E, 0 E, <f) Ez - B ^ yz plane xy plane — B * ------- B *--1— B ^ - : 0 (|) 0 r § 0 § 0 (|) 0 (!) ^ ^ - B > ----- B*—|— B*-j § 0 ^ ------- <$>------ <$>---------<f>-------- <f) B *~ 0 § 0 ,\ (f> 0 § / 0 / <t> 0 / -- § 0 0 0 □ > - L- □ > 0 -- -B ^ - - -B*~ - - —B* 0 0 0 B*- 0 — B>— -Q *-A - 0 0 ■B*- x 0 0 □> --—B*—:—B*—-—B*—' 0 zx plane Figure 2.5 Implementation of metallic walls and edges 40 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -> 3. Absorbing Boundary Conditions (ABC’s) A fundamental problem in finite-difference time-domain algorithms, which are used to solve electromagnetic field problems, is to simulate “open” boundary conditions, where the computational domain should be of infinite extent. The storage capabilities, of all computers, are not unlimited. It is therefore desirable to limit the computational domain to a volume just large enough to contain the structure of interest. In order to simulate “open” boundary conditions, a suitable method for truncating fields incident on the lattice boundary must be added to the main differencing algorithm. These boundary algorithms are called Absorbing or Radiating Boundary Conditions (ABC’s or RBC’s, respectively). In order to be effective ABC’s must suppress spurious reflections, due to the truncation of the infinite volume, to acceptable levels, permitting the time domain solution to be valid at all time-steps. Several methods have been employed for the development o f ABC’s. Initial efforts were based on the factoring of the one-way wave equation, which were later followed by the implementation o f the Perfectly Matched Layer (PML) ABC’s. 41 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.1 Wave Equation ABC’s A partial differential equation, which mathematically describes wave propagation only along a certain direction, is called a “one-way wave equation”. The derivation o f ABC’s based on the one-way wave equation can be explained in terms of operator factoring, which will be discussed following the approach proposed in [44]. For example, the 2-D wave equation in Cartesian coordinates can be written as: (3.1) where U is a scalar function (or a single component of a field vector). The subscripts xx, yy and tt denote second partial derivatives with respect to x, y and t, respectively, and c is the wave velocity. The partial differential operator, associated with equation (3.1), is given by: (3.2a) which uses the following notation: (3.2b) The wave equation can then be written in terms of operators as: LU = 0. (3.3) 42 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The wave operator itself, can then be factored in the following manner: L U = L*L~U = [ d x + y V l - S 2) ( £ * - y V l - S 2) f /, (3.4) with S = cDy / D,. Engquist and Majda [47] showed that at a planar boundary, x=0 for example, the application of L~ to the wave function (U) will absorb a plane wave incident on the boundary at an arbitrary angle 9. Thus: L'U = 0, (3.5) applied at x=0 functions as an exact analytical ABC, which absorbs the incident waves, and prevents them from being reflected back into the interior of the volume V. The presence of the radical in L' is undesirable, as it does not allow for a direct numerical implementation of (3.5). Therefore, the radical is approximated using a two term Taylor series expansion: y l\-S 2 =l - S 2/ 2 . (3.6) Substituting (3.6) in (3.5), an intermediary result is obtained: ' D. A- + c cDl U = Q. y 2D,) (3.7) 43 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The multiplication of (3.7) by D/ and the replacement o f the differential operators in terms o f partial derivatives produces the following approximation to the one-way wave equation: UJt- ( l / c ) U „ + ( c / 2 ) U yy=0. (3.8) Equation (3.8) is a fairly good approximation to (3.5), if S' is relatively small. The derivation o f the ABC’s in three dimensions follows the methodology outlined above, with S given by: (3.9) 3.1.1 Mur ABC’s A simple finite-difference implementation of the ABC’s given in equation (3.8) was introduced by Mur [12]. For clarity, Mur’s scheme will be illustrated for the 2-D case, at jc=0 point o f the grid. The Mur scheme determines the electric field value E"(0,j) at a lattice node on the boundary, approximating the partial derivatives with central differences, expanded around the intermediary field E'y( l / 2 , j ) component, located one half-space cell from the lattice boundary. 44 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The first step in the derivation involves rewriting the mixed partial x and t derivative using central differences: £ ; - ' ( i , ; ) - £ ; l(o,y)'| ^ 2AxAt ^ 2AxAt (3.10) ^ Next, the partial tt derivative is written as an average of time derivatives at the adjacent points (OJ) and (l,j). The averaging is needed as E'y( l / 2 , j ) i s not a grid point and its value can be estimated only by computing the mean value of the field at the adjacent points: f E';' (0, / ) - 2E\ (0, J) + E';' (0, j ) (£;)(/a/2j) = i/2 (A/)2 (3.11) Similarly, the partial derivative with respect to y of £ '( 1 /2 ,/) is written as an average of y derivatives at the adjacent points (0,j) and ( //) , respectively. Substituting the above finite-difference approximations for the derivatives in equation (3.8), and solving for £ '+1(0 ,/) , leads to the Mur time-stepping algorithm for the ABC’s, shown in equation (3.12). Analogous finite-difference expressions for the Mur ABC’s on the other lattice boundaries can be derived as well. 45 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (3.12) ^ (cA/y C y {V,J + l ) - Z £ y {V ,J) + £.y( 0 , j - l ) ' 2S(cAr + S ) [ + E,y (l,j + l ) - 2 £ ; ( l J ) + E ; ( l , j - l ) J Extension o f the Mur ABC’s to 3-D follows the above discussion closely. Note that, in essence, the resulting 2-D and 3-D algorithms take into account not only the variation of the field variables along the direction of propagation, but also their variation in the transverse plane. 3.1.2 Higdon ABC’s Higdon proposed a simplified algorithm, based on the one-way wave equation. In this case, the derivatives corresponding to all directions perpendicular to the direction of propagation are neglected. The initial form of the boundary condition was given in [25]: (3.13) By numerical experiment, it was discovered that the absorption of reactive modes was improved if attenuation factors were added in equation (3.13), which lead to the expression for the lossy Higdon’s ABC’s: (3.14) 46 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. where N is the order o f the ABC, 9j are the plane wave incidence angles to the boundary, and the Q>0 are attenuation terms, which may be used to increase the absorption of the reactive fields. Equation (3.14) can be discretized in several different ways. With Oj-O and N = l, the Higdon algorithm is equivalent to the first order Mur ABC. Adding the dependence on the incidence angle and the attenuation factor, the first order Higdon ABC’s take the following form: E'+x(0, j , k ) = (l —£ ) E ‘(1, j , k ) + acAt - Ax acAt + Ax where C, is a small positive number corresponding to an energy loss, which is being chosen to achieve desired losses. The constant a is related to the assumed incidence angle on the outer lattice boundary. The extension of Higdon ABC’s to higher orders can be performed efficiently by using the coefficient matrix notation described by Luebbers in [25]. For the second order Higdon ABC, the coefficient matrix method yields: If it assumed that (1 - ) = (1 - ) = a < 1 and that b{ =b2 = a - A t —Ax then equation a • At + Ax ’ (3.16) becomes: 47 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -1 R = -2b -b2 2b -b2 2(cr +b2) - 2 a - b 2a-b (3.17) -a2 If the ABC operator R is applied to the z=n boundary, where n is the maximum lattice size along the z-axis, the second order Higdon boundary condition for the electric field components Ex and Ey, tangential to the boundary, can be written as: £ '+l = 2ct • £ ' (i,y,/i - 1) - a 2E'~x( i,j ,n - 2) + 2b (3.18) - b 2( E ,+l (i , j , n - 2) - 2 E ‘ (i , j , n -1 ) + E " 1 Similar expressions can be derived for the other boundaries of the lattice, and other components of the fields that are tangential to them. Compared to the Mur ABC’s, the Higdon ABC’s do not use values o f the electric field that are located at the neighboring nodes on the tangential planes to the direction of propagation. Since the introduction of Mur and Higdon ABC’s, several other methods, based on the factoring of the one-way wave equation were proposed, and they include the Superabsorpting ABC’s [48] and Liao’s ABC’s [49]. The reflections due to lattice boundary truncation by all the aforementioned ABC’s are in the order of -20dB to -40dB. In addition, they all depend on the angle o f incidence and on the fact that the medium has to be homogeneous at the boundary plane. The lowest reflections, according to a study 48 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. done by Luebbers [25], as shown in Fig. 3.1, were for incidence angles of about 25° for 2ndorder Mur, and for 10° for the other ABC’s. Magnitude of Reflection Coefficient vs. Incidence Angle -10 •20 O -30 © 2nd Order Mw Super Slablized 2nd Order Liao UnsUbtizcd 2nd Older Liao Sublizcd 2nd Order Higdon Unstablized 2nd Order Higdon 50 60 70 Figure 3.1 Reflections of one-wave ABC’s for oblique incidence from [25] 3.2 Perfect Matching Layer (PML) ABC’s In 1994, Berenger [23] proposed a new type o f boundary conditions for use in conjunction with FD-TD, initially implemented for the 2-D case, and immediately extended to the 3-D case by Katz and Taflove [17]. The idea is to introduce a set of non physical absorbing coefficients, which yield very low reflections for uniform dielectric 49 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. layers, which are added in front o f the metallic boundary planes. In order to apply the PML truncation algorithm, the computational volume, where FD-TD is employed, is surrounded with thick walls, in which the waves propagating normal to their surfaces are artificially attenuated. In order to avoid reflections, the equality: (3.19) must be preserved at every point within the homogenous media. The PML algorithm increases the number of computational variables, by splitting the electric and magnetic field components into two subcomponents. As a result, the two Maxwell’s curl equations must be replaced by a set o f 4 equations in the 2-D case and 12 equations in the full 3-D case. For example, the equations corresponding to the 2-D TE (transverse electric) polarization case, where the magnetic field is split as H . = + H , are: (3.20) dE y dx (3.21) dH dE X Designating ¥ to be any component of a wave propagating in the PML medium, Berenger showed that: 50 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. jr-cos++y-siii+ q r cos+ q y sin+ cG ZqcG e 0cG (3.22) where c is the speed of light, <t>is the angle between the electric field and the direction of propagation, with the remaining quantities given by: G = ijwx cos2 <j»+ wy sin2 <|>, 1 ~ j G a / COS0 1 - y ^ / ^ o ’ _ ' (3.23a) l-y ^ /C O S p (3.23b) l-y'CT^/CDUo If each pair of variables (^.s.cy^) and ,o mv,j satisfies the constraint imposed by (3.19), then G = wx= wy = 1 at any frequency. Consequently the field components and the wave impedance become: x cos^-t- v sin ^ co sf ffyStn^ (3.24) Equation (3.24) shows that the wave in the PML medium propagates with the speed of light, and decays exponentially along x and ^-directions. Notice that the wave impedance in the PML layers is the same as that o f free space. The implementation o f the PML algorithm in 2-D, for absorbing boundaries normal to the jc-axis, and no losses along the y-axis, requires a y=0, and for <rx to increase along the x-axis, with the increasing penetration depth within the wall (see Fig.3.2): 51 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Regular computational space -> x > Figure 3.2 Regular/PML computational space (3.25) ° * ( * h ■all) ~ ® x m a x { X "a ll ^ ?m ill) where t va„ is the PML thickness. This theoretically yields a reflection factor of: (3.26) R(Q) = exp(- 2 a maxflra// cosG / (n + l)s0c ) . The extension of this procedure to the 3-D case is straightforward, as was shown in [17]. Numerical results, for uniform dielectric at the boundary plane, show reflections o f -27 dB for a 4 cell PML wall, -55dB for an 8 cell wall and -75 dB for a 16 cell wall. If the medium wherein wave propagation is taking place is composed of layers with different dielectric constants (see Fig. 3.3), then the artificial magnetic losses are kept constant (am = constant) and the artificial losses in different dielectrics are forced to satisfy the following condition [18]: 52 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. FDTD volume P M L w a ll cl.nO e2,u0 E l,|i0,crel ,cjh 1 wave direction --------------- > e3,h 0 E2,H0,ac2,OHl c3,hO,oe3,ohI Figure 3.3 Implementation of PMLs for inhomogenous media Application of the split-variable PML algorithm has several drawbacks. First, the extra variables require additional memory and additional numerical operations at each lattice node. The implementation becomes rather complicated, as the extra variables have to be placed in, at least, 12 supplementary arrays. For some problems in FD-TD, and for ESN in particular, the PML algorithm is not stable. Due to these factors, a simpler PML implementation was sought for boundary truncation in this dissertation. In order to place the new PML algorithm in proper context, it is appropriate to review the present state o f the art. Several different variations o f the unsplit PML algorithm were 53 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. recently proposed in the literature [30,31]. The approach proposed by Zhao and Cangellaris [31] will be briefly described below. First, the authors have shown that the PML algorithm is equivalent to the Chew-Weedon perfectly matched medium, which is obtained from Maxwell’s equations, by stretching the coordinates. Secondly, the simplified PML algorithm, with unsplit variables, was derived. This version o f the PML algorithm introduces changes to only one component of Faraday’s law, when wave absorption along the z-axis is assumed: = -j<oHx - ^ - H z , e p i dy dz 1 (8E, — (i v dz dE.\ = -j(oH - - H dx J e 1+ - ^ \ jasj dEv \ dx (3.28a) , (3.28b) dE, dy In order to develop the PML algorithm, equations (3.28) must be transformed to time domain. In the above expressions, (3.28a) and (3.28b) have the standard frequency domain form o f the magnetic field curl equation for a lossy magnetic medium, with the “magnetic conductivity” Equation (3.28c), corresponding to the component of the H field, which is normal to the interface between the computational domain and the PML layers, when transformed in time domain, has an additional integral term, as shown in (3.29): 1 r dEv dx dE' dy) dH. dt a J7js(t)rfc. ep ' (3.29) 54 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Finally, equation (3.29) is discretized using central differences and is incorporated in the regular FD-TD algorithm. The numerical performance o f the simplified PMLs is similar to that o f the original scheme, yielding reflections of -43 dB for a 4 cell PML wall, and 87 dB for a 8 cell PML wall, respectively. It is important to point out that these reported results are for free space and homogenous media, where PMLs are expected to perform the best. However, in realistic problems, where the propagation medium is composed of several dielectrics, none of the existing PMLs can approach even close to the results that are typically quoted for the ideal problem. 3.3 Proposed New PMLs fo r Inhomogenous Media To overcome the difficulties associated with the existing PMLs, a new approach to open lattice boundary truncation is proposed. This method allows for computing the “loss” coefficients for an arbitrary geometry (1-D, 2-D, or 3-D) and excitation. A detailed implementation of the procedure is presented for ESN in one dimension. It is shown that the results remain valid for 3-D problems involving quasi-TEM modes o f propagation. The resulting reflections from lattice truncation boundaries range from -90dB to -52dB for 1-D and 3-D problems, respectively, with as few as 3 PML layers. As a starting point, consider a general form of any finite difference time domain algorithm, that can be used to solve a coupled system of first order hyperbolic differential equations (such as Maxwell’s curl equations), which can be written as: 55 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. f £ ) '= E K ] f £ > " ‘ + ('S>'- (3-30) *=l In equation (3.30), {£}' represents a vector that contains the values of all field variables at time step while [/?*] (k = 1, 2, t) contains a set of matrix operators that are specific to a particular finite difference algorithm used to approximate Maxwell’s equations. The excitation, applied at time may be located anywhere within the computational lattice (including its boundary), and its values are stored in {S}'. In time domain algorithms, equation (3.30) is used to advance the field in time, based on the field values at the previous times steps. However, the field values at any time “r” can also be expressed in terms of the excitation, {5}' ( k = 1, 2, ..., t), directly, as shown below: W = t [ C kp } k > *=i (3-31) where the [ct ] arrays can be determined by utilizing recursive relationships given in equation (3.30). To determine the absorption coefficients, two one-dimensional problems, shown in Fig. 3.5, will be used. One of them is a fairly long microstrip, which is for simplicity, terminated with a metallic wall (Fig. 3.5a). This corresponds to the uniform section o f the (output) microstrip o f Fig. 3.4, which is sufficiently far from the discontinuity (~A.g/4) and 56 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. carries a quasi-TEM wave. A reference location is selected at distance L sh o rt from the metallic wall such that the time history of the entire field { e ^ p }' can be recorded prior to the arrival o f any reflections from the wall. The second problem contains a shorter section of the same transmission line, but with three PML layers added between the field recording point and the wall, as illustrated in Fig. 3.5b. The time domain simulation is repeated and field values { f^ recorded at the same location. The initial values of the absorption coefficients within the PML layers for this computation are chosen according to the “rules” outlined in [29]. The time domain simulation of both problems is carried out over the same number of time steps and a measure o f the error due to the PMLs is computed from: (3.32) *=1 after each simulation q. Note that unlike }/ > which has to be calculated Q times until the error {e} is minimized, {£/^£F}, is computed only once. Following every simulation q, a new set of absorption coefficients is determined from equation (3.32) with the help of multivariable minimization algorithm, such as the Nedlear-Mead Downhill Simplex method [39]. 57 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. lattice boundaries incident wave transmitted wave ; i n f o n n rricnostrip/ se c tio n s ground plane w Figure 3.4 Generalized PML problem setup The resulting absorption coefficients are used in PMLs to recompute in order to repeat the minimization. This process is continued until the user-specified error (e.g., {e} =0.0001) criterion is met. It is important to point out that since these are 1-D problems, time domain computations are carried out very quickly, placing minimum burden on computer resources. The above procedure can be used in one of two ways. Families of absorption coefficients can be precomputed as functions of geometry and/or material parameters and later used in PMLs. 58 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. "Y S source y / dielectric y ✓ y ✓ ^ s ; y s / s / / s / .’ / / ' l short t Lref 4-------------------------- ------------------------------------------------ (a) 3 PM L {S }1 , AT , layers air ^ source | {EL0sS> '' '~ -l ( ' / dielectric ' y (b) Figure 3.5 Reference (a) and PML (b) problems Alternatively, this exercise can be implemented as a pre-processor to the general time domain solver, so that appropriate “loss” coefficients could be calculated for a specific problem geometry. In addition, it was found that the absorption coefficients extracted from a 1-D problem yield good results in 3-D as well. This makes the overhead of 59 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. optimizing the absorption coefficients a worth while effort, considering that very low reflections (-50 dB and lower) can be obtained with as few as three PML layers (as opposed to -80dB for 16 PML layers reported previously in literature [50]). 3.3.1 ESN-Specific Implementation To illustrate the details of the proposed technique, its implementation in ESN will be outlined. For clarity of presentation, only the 1-D formulation will be considered. Since the proposed optimization scheme simulates a TEM wave along the transmission line in 1-D, it should be understood that the results in three dimensions are expected to be valid for quasi-TEM fields. In addition, the truncation boundary in the actual problem is assumed to be normal to the microstrip, with the material parameters and excitation being the same as in the 1-D loss and reference problems (Figs. 3.5a and 3.5b). Now consider the time-stepping ESN equations in one dimension, derived from equation (2.17), with the losses implemented as in (2.18): + 2gce(/)£ '-'(/) + (3.33a) f-3 + 4g«a)2(-i)‘+lr~*+,(02 ( / / '- ’(i +1 / 2) - H'~' (/ - 1 / 2)), 60 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. H'(i) = (0/ / '- 2(/) - 2 / v (£ '-‘(/ +1 / 2) - E ' - \ i - 1 / 2)))e'°-(0, (3.33b) where gce(i) = 2v(er - l ) . When the values o f the electric and magnetic field components at every point in space are combined into a single vector {£}', the recursive time-domain ESN algorithm can be rewritten in the following matrix form: {£}' = [£,]{£}'-' +[S1]{£}"! + [ B £ ( - 1 ) * * '{£}'-**' +{£}'. (3.34) *=3 where the [5], [/?,], and [i?2] arrays can be readily obtained from equation (3.33). In addition, if the above expression is rewritten in terms of excitation, as in equation (3.31), a recursive relationship for the computation of the [C] arrays is obtained as well: [C,' ] = [£,] [Q -1] + £Bj ][C,'-2] + c—I)**1[CV'*»' ] (3.35a) k =3 [C'] = [C;-*+1] (3.35b) [C,2] = [51]and[C,l] = [/] . (3.35c) 61 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For optimization, within the ESN, the reference problem consists o f a lossless transmission line with substrate permittivity zr, whose length is L ref - nAx (see Fig. 3.5a). The corresponding loss problem is composed of a similar transmission line, whose length is L lo ss - wAx (where m < n). For both geometries, the fields due to the same excitation are calculated at the reference point one cell in front o f the first PML layer and are recorded at every time step At. It was found that collection o f field data at a single test point “ptest” (rather than on the entire ESN discretization line) was sufficient to obtain the desired “loss” coefficients that lead to low reflections. Finally, with the above notation, the minimization function can be written as: (3.36) where S(k) = e is the kfh element of {S}. The absorption coefficients a e(i) and a m(/) (see equation (3.33)) can now be calculated from equation (3.36) using the Simplex method applied to the error function s. As pointed out earlier, the initial values for a e(/) and ccm(i), which lead to quick convergence, were chosen according to “rules” outlined in [23]. As can be seen from equation (3.36), the “loss” coefficients will depend on the material properties and geometry of the problem (all o f which are embedded in the [C] arrays), as well as on the applied excitation. 62 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.3.2 FD-TD Implementation In order to show the applicability of the above proposed method for various finitedifference schemes, an outline o f a possible path towards implementing it in FD-TD is provided. For example, consider the time-stepping equation for Hx: " J _ a A/n 2M f E " ( i, j + l / 2 , k +1) - -i / Az + A/ j + l / 2 , k + l\ 2/J J ,(3.37) +■ /A y with losses depending on the material parameters cr, and a m. The [5] and [C] arrays can be obtained from the above expression in a similar manner as the procedure outlined earlier for ESN. Once the FD-TD algorithm is recast in terms o f the excitation, the corresponding error function (3.36) can be formulated. Since the attenuation in FD-TD depends on the material conductivities of the electric and magnetic media, the Simplex method can be employed to minimize the error function using material parameters and a m as variables (unlike the a ’s that were used in ESN). 63 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.4 Numerical Results o f Different ABC’s fo r ESN The first attempts to incorporate ABC’s in ESN were based on the one-way equation method. Both Mur and Higdon absorbing boundary conditions were tested for modeling open boundary conditions within the ESN algorithm. The emphasis of the study was placed on their behavior when applied to transmission line type structures. The Higdon absorbing boundary conditions performed better both for structures having uniform and non-uniform dielectric, respectively. Sample representative results, displaying both the excitation and the reflections are shown in Fig. 3.6. It can be seen that the Higdon absorbing boundary conditions have reflectivity o f about -23dB and -16 dB for uniform and non-uniform dielectric, respectively. Although these ABC’s were not very good at absorbing waves normal to the absorbing surface, they proved to be very useful at simulating open boundaries at grazing incidence, when the wave is propagating parallel to the ABC plane. In this case, it was observed that the wave was losing very little energy along its path. As a result, the second order Higdon ABC’s were used in all simulations for boundaries parallel to the direction of propagation. 64 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.8 0.6 0.4 0.2 125 - 250 375 500 625 750 875 0.2 Time steps Figure 3.6 Mur and Higdon ABC’s reflections for open microstrip line As mentioned earlier, attempts of using the spilt-variable PML algorithm in ESN failed, as it proved to be unstable. On the other hand, the PML algorithm with unsplit variables yielded better results, which are presented below. The integral term necessary for the computation of the H~ component of the magnetic field (see equation 3.29) was assumed to be zero, due to the quasi-TEM nature of the problem. The reflections at the interface were assessed by using three different uniform microstrip transmission lines. One consists of a strip suspended in free space, the second o f a strip placed on a grounded substrate with er=2.2, and the third is identical to the second, but is covered with an upper layer of er=10.2. 65 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The geometrical details o f the problem are shown in the inset of Fig. 3.7. The discretization along the x, y and z axes is uniform with A=0.25mm and the computational space is subdivided into a 32x32x60 lattice. The ESN absorption factors are chosen to be aa - a ey = = a my = 0.1 and the PML wall thickness is 32 cells. The reflections, shown in Fig. 3.7, are smallest for uniform media (R=En.flectc<1/Ein<;,denl =0.9%), and increase with the increasing number o f dielectric layers. In the case of a single substrate R is 1.9 % and increases to 3.1% with the addition of the cover layer. The performance of the unsplit-variable PML ABC’s is acceptable, but does not yield sufficiently accurate results. Since it is desirable to compare simulation data to high precision measurements, the artificial reflections must be below 1% for layered dielectric problems, like the open microstrip line. The implementation o f the new method, proposed in section 3.3, for determining the “loss” coefficients, reduced reflections significantly, bringing them below 0.25% in the range of 0.5GHzto 10GHz. 66 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1.2 ♦-FreeSt) .0 - a - ESN Che -ft-ESN Two -x-HgTwo Q) ■a c O) CO E a4 x LJJ X X-X x -x x-x-X-X-X-X-X-X-X-X X -X 'X X X - 0.2 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 0 100 200 300 400 500 600 700 000 900 1000 Mirfcer of time steps Figure 3.7 Reflections from PMLs and second order Higdon ABC’s To assess the numerical performance o f the new absorption coefficients as a function of the substrate permittivity and the form of excitation, a one-dimensional transmission line whose length is 60Ax and Ac = 0.2mm was taken to be the reference problem (see Fig. 3.5a). The time domain excitation is a Gaussian pulse s(i) = e ,(x = 3a) , which is applied (in the from of the electric field) under the strip, is allowed to last 300 time steps (A/ = 0.167 psec). The complimentary loss problem (see Fig. 3.5b) was composed o f a similar transmission line that is only 8Ax long with three absorption layers, 3Ax in thickness. In both problems, the “test” or field recording point is located one spatial cell in front of the first PML layer. 67 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The magnitude o f the frequency domain error between the fields in the reference and the loss problems is computed from: z(kAf) = ( E ^ i k A f ) - ELOSS(kAfj)f EREF(kAf) (in %) (3.38) for frequencies ranging from 0 to 20GHz. The corresponding absorption coefficients for the microstrip are listed in Table I for typical combinations of the substrate permittivity and pulse width. To give a better idea o f how the “loss” coefficients vary with the substrate permittivity, a family o f curves were generated for er ranging from 1 to 12, with all other parameters fixed (Ax = 0.2mm . <r = 30A/). The corresponding data for this case study are shown in Fig. 3.8. Note that the values o f the absorption coefficients tend to increase with increasing values o f the dielectric constant of the substrate. This, however, is not unexpected, as the higher values o f er cause greater departure of the field form the TEM distribution. To evaluate the effectiveness of the loss coefficients obtained for the 1-D case in three dimensions, an open microstrip line shown in the inset o f Fig. 3.9 was excited with a Gaussian pulse (ct = 30A/ and x = 3 a ). Three PML layers, each one A wide, with the loss coefficients selected from Table I were employed for lattice truncation. 68 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. T a b l e I. R e f l e c t i o n s a n d Lo ss Ex c it a t io n C o e f f ic ie n t s f o r Param S CT 1.02 30 elected G e o m e t r ic a l a n d eters Reflection (%) a e[ = -2.01 O'4 a m l = 5.7-10'3 1.7-10'3 a e2 = 3.3 -1O'2 a m2 = 9.68-10‘2 a ej = 2.12-10 '1 a m3 =4.752-10'' 2 30 4 30 10 30 12 30 a ej = 3.5-10'3 a m j = 1.31-10*2 a e2 = 2.43-1 O’2 a m2 = 5.66-1 O'2 a ej = 1.075-10' 1 a mj = 3.081-10'1 a ei = 2.3-10'3 a m / = 1-2-1 O'2 a e2 = 2.79-1 O'2 a m2 = 5.86-10'2 a ej = 1.075-10' 1 a mj = 2.288-10-' a ei = 1.9-10'3 °-ml = 7.9-10'3 a e2 = 1.6-1 O'2 a m 2 = 3.11-10'2 a e3 = 5.49-1 O'2 a m i = 1-376-10'1 a ei = 1.4-10'3 a ml = 7.0-10'3 a e2 = 1.68-10*2 a m2 = 3.44-10'2 a ej = 6.23-1 O'2 a m3 = 1.326-10'1 -1.8-10'3 -2 .5 -10~‘ -1.2e-10’2 -2.5-10'2 The time domain simulation was performed for 1000 time steps (A/ = 0.167 psec) for the reference and loss problems defined in Fig. 3.5, but now in three dimensions. The size of the lattice used in the reference and loss problems was (12A x 50A x 140A) and (12A x 50A x 70A), respectively, with A being 0.2 mm. Signal waveforms were recorded at the testing point “ptest” for both reference and loss problems, and the error in the frequency domain was computed. The errors are plotted in Fig. 3.9 from which it can be seen that they are higher than for the 1-D case (-51 vs. 98dB). 69 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.50 I V 0.45 -A -C l c2 X -A—c3 0.40 X \ 0.35 -o ~ c 4 -o -c 5 \ —o —c6 '« X 0.30 0.25 0.20 ' 0.15 a— 0.10 _ , ° ----- A 0.05 <--------------- ' -------------- -------a ---------------------------- -------------- A----------— 0.00 T ^ T ^ -T —------------r 4 — 1 2 ------- a ---------------■ ■ 4 ----- ,1 10 12 R elative p e rm ittiv ity Figure 3.8 "Loss" coefficients vs. relative permittivity 025 02 0.15 005 UJ •0.05 •0.1 R e la tiv e p e r m i t t i v i t y « 2 , 4 - 0.2 06 24 42 6 78 F r e q u e n c y (G H z) Figure 3.9 Frequency domain reflection errors for the numerical PMLs 70 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.6 Such behavior was expected, since the wave propagating in the uniform microstrip is quasi-TEM, unlike the pure TEM which was used in one dimension. Irrespective of the higher reflections, the errors are still well below 1%, which is impressive, since only three PMLs were used. Next, the influence o f the microstrip geometry on the artificial reflections from the lattice boundary was assessed. 0.4 03 ^ 02 0C3 o t3 0) o.i £ o ® 3 TJ LU 0 D" ,• «r \ •0.1 -O- w/h=10/3 -*-w/h=6/3 -O—w/h=2/3 -02 -0.3 0.6 3.0 54 7.8 F req u en cy in GHz Figure 3.10 Dependence of reflection error on the w/h ratio The dielectric constant o f the substrate and the absorption coefficients (er = 10, a e\ = 0.0019, a mj = 0.0079, a e2 = 0.016, a m2 = 0.0311, a ej = 0.0594, a mj = 0.1376) were held constant, while time domain simulations for several w/h ratios were carried out. Fig. 71 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3.10 shows the errors due to artificial reflections for the open microstrip o f Fig. 3.9, whose substrate thickness is h = 3A, but whose width w is 10A, 6A, and 2A, respectively. The simulation data show that at low frequencies (< 5.5 GHz) the errors increase with decreasing w/h values. At higher frequencies, for all w/h ratios, the errors are similar in value, up to 0.28%. The increase o f the error is due to departure from the TEM wave model, behavior already reported in literature [50]. 72 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4. Subcell Models for Active and Passive Circuit Elements The utility o f time-domain methods in simulating wave propagation in waveguides and microstrip transmission lines has been demonstrated in the past. In microwave amplifier circuits, such guiding elements provide interconnections between discrete passive and active components, such as resistors, capacitors, diodes and transistors. In order to simulate microwave circuits containing lumped and distributed elements, time-domain methods must be modified accordingly. To date, lumped circuit element models have been developed for FD-TD only. Initial efforts concentrated on incorporating passive elements and diodes into the 2-D FD-TD algorithm [32] and were later extended to 3-D [33]. Subsequently, small-signal subcell models for the BJT’s have been developed for FD-TD as well [33]. As it was already pointed out, ESN is a finite differencing algorithm, derived from a transmission line network to which lumped circuit elements are added. The addition of lumped circuit elements to ESN, modifies only the specific cell to which they are added. In this chapter the focus is on incorporating active and passive lumped circuit elements such as diodes, resistors and capacitors into ESN, using a similar approach to that 73 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. proposed in [33]. A new way of incorporating two-port devices (FET’s and BJT’s), that is based on their transient V-I port relations, into the finite-difference algorithm is also proposed. 4.1 Introduction to subcell models Analytically, lumped elements can be introduced into ESN using the integral form of Ampere’s Law: (4.1) This equation is applied to the appropriate components of the electric field, within the computational cell of the FD-TD or ESN algorithm, where the lumped load is to be placed. The effect of the load is included as an additional current flowing through the surface I c, shown in Fig. 4.1. Its value is given by the functional relationship I=ftV), which describes the electrical behavior of the lumped circuit element. As a specific example, consider a lumped load, placed along the x-axis, which is piercing surface Zc (see Fig. 4.1). To properly introduce lumped loads into the ESN, its finite-differencing algorithm must be related to Maxwell’s curl equations. 74 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. x load Figure 4.1 Inclusion of lumped loads along x-axis If equation (4.1) is specialized to the x-directed load, shown in Fig. 4.1, the x-component o f Ampere’s Law becomes: dHv dH.. c ~ oz r " ~ ~ oy ^ r = aE*+ , dEx + 8 ° ■'~otr r + 8° "(e - _ 8EX at (4.2) In the case of uniform air dielectric, and without any lumped loads, the ESN algorithm can be written as: f~ - 1 / 2,*) + V ^ U +1 / 2,k) \ V m v k - \ I T ) - V ' my( i , j , k + \!2) ) where u = Ax/(c0Af) and a a ( i, j ,k ) = ^Jn0/ e 0(aAx)/ (u(2er(/,y,& )- 1)). 75 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.3) After normalizing equation (4,2), the one-to-one correspondence between it and equation (4.3) can be outlined. If the dielectric is not air, then the third term on the right hand side o f equation (4.2) is approximated using central finite differences in FD-TD. In ESN, the effect of this term it can be added as a node current through a capacitors (see equations 2.5 and 2.6 and Fig. 2.2). For convenience, equation (2.16), depicting the ESN algorithm, is rewritten below: ((u + g ca (/, J\ £))(! ~ a « (/, j , k))) ■ K ~ ' O'J. * ) + 2 • gccx 0\ / » k ) • K O'. 7. k ) + C '0 \ M ) = - O ' J - 1 / 2,*) + O J +1 / 2,k) (4.4a) +2 U 4 g ca( i , j , k ) - V exsum e/-2 1 ((« + g m O 'J. ^))(! + a « 0 , J , *0)), where: (4.4b) =C U M ) - O U M ), & » 0 \M ) = 2 -u -(e r - l ) , (4.4c) u = Ax/(c0A/) , (4.4c) = V^o/®0 (aAx) 1(U(2er0‘»y» *) - 0) (4.4d) The load current has to be expressed in terms of electric field intensities and current densities. Its inclusion into the regular ESN algorithm is accomplished by using finite 76 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. difference approximations to the V-I relationships obtained from Maxwell’s equations. The specific procedure for implementing lumped active and passive circuit elements will be shown below. 4.2 Passive Circuit Element Models 4.2.1 Resistors For a resistor placed along the x-axis, the V-I relationship is given by Ohm’s law V=RI and JxL, at any given time /, can be expressed in terms o f finite differences as: (4.5) The expression for the load current density (JxL) is similar to the current density due to dielectric losses at the given node (aEx), as can be seen from equation (4.2). If the current through the lumped resistor is along the x-axis, within one cell, then the term OLex(iJ,k) must be modified, in order to include the additional losses due to R: (4.6) It should be noted that equation (4.6) is expressed in terms o f the normalized ESN variables, as described in [36]. If the lumped resistor must be distributed over several ESN cells, then an appropriate (Ohm’s Law type) V~f(I) relation should be used to 77 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. account for the current flowing in each cell. The value of R, which is to be used in a a of each cell occupied by the resistor, must be modified accordingly, taking into account the rules of series and/or parallel resistor combinations. For example, a lumped resistor terminating an open microstrip line (see Fig. 4.2a) can be implemented by distributing its value across the width of the microstrip and along the thickness of the dielectric, as in Fig. 4.2b. Since the width of the strip and the thickness of the dielectric are both equal to 3A, the value of the resistor for each cell will be equal to the original value o f the undistributed resistor. 4.2.2 Capacitors For lumped capacitors the V-I relationship is defined as I = C(dV/dt). The corresponding expression for the additional current density Jxi within a single ESN cell, where the capacitor is located, is given by: J'l L ={C/Ax)-dEx / a . (4.7) It is assumed that the capacitor is placed such that the current through it flows in the xdirection. To incorporate equation (4.7) into ESN, the x-component of Ampere’s law is used: 78 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. which suggests that the lumped capacitor should be implemented in ESN in the same way as non-air dielectrics. The effect o f the capacitor is taken into account by modifying the dielectric material parameters in the lattice cell occupied by the capacitor. microstrip (metallization) lumped resistor zoomed in area dielectrtic a) A b) Figure 4.2 Implementation of distributed resistors 79 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. As a result, in the cell where the lumped capacitor is placed, the gj(i,j,k) term takes the following form: gca{ i j , k) = 2 •u •((sr - 1 ) + C / (Ax •£„)). (4.9) As for the resistor, if the capacitor has to be distributed across several cells, then series/parallel equivalent capacitor combinations should be used. The values o f the individual distributed capacitors should be modified accordingly, such that their series/parallel combinations add up to the total (lumped) capacitor value. 4.2.3 Inductors In order to incorporate lumped inductors into ESN, the following V-I relation should be used: (4.10) 0 The current density through an ESN, in terms of normalized variables, with an inductor placed along the x-axis is given by: / y [r(/,y,it) = (p0 / z ) - X ^ O ‘J ^ ) *=i (4.11) Equation (4.11) can be used as is in equation (4.1), replacing Ekx byV£, but it has the drawback o f requiring the entire time history of the electric field V^(J ,j, k) . As the 80 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. storage o f such large amounts of data is not feasible, the sum in (4.11) can be computed at each time step by using an intermediate variable Vim( i, j,k ) = V ^ ( i , j , k ) + V‘x (i,j,k) for each o f the cells containing the inductor. 4.2.4 Voltage sources Assume that a voltage source is characterized by a time domain waveform Vs and having an internal resistance Rs. If it is placed along the x-axis, the current density that goes through is: V E J * = iRr r Ax t T + TR rlAx r- (4-12) When equation (4.12) is normalized (Va = y j ^ E x ), its addition to the time-stepping algorithm for V ^ (i,j, k ) leads to: [(« + Scex O'* . Jgs, **, )XJ + a « O'*. J'gs •k v ))] • C + 2 ■gca ( v , j „ , ^ • K O'p . / „ , * „ ) + 4gccc 0 * f - C O W * " 112’U + (If, > /*,.**) = C l . + O W * + 1r 2’d 2- \KyOg, J g ^ k g , - 1 / 2 ) - (ip J ^ k g , + 1 / 2 ) where: 81 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.13a) v = Ax/(C0A/) , (4.13b) a a (/, J, k) = V /V ^(o A * + 1 1 R ' ) / (v{2sr(i, j , k) - 1)), (4.13c) g caU ^ j ^ k gs) = 2-u-(er - 1) , (4.13d) K L ( . iv J v ^ ) = V ^ \ i ^ J p , k p ) - V £ m(iS!, j gs,kgs) . (4.13e) The lumped voltage source, especially its DC form, will be useful in the analysis of transistors, as will be shown later. 4.3 One Port Active Circuit Element Models (Diodes) In the case of diodes, their nonlinear V-I characteristic is described by the equation: ij ~ ~ l) • The diode has a nonuniform voltage distribution inside its packaging. For the model developed in ESN, it is assumed that the diode has punctual physical dimensions and is contained within a lattice cell. The current density of a diode assumed to be placed along the x-axis is given by: J ‘«. = (/0/(Axy)-(s’“ ;' ‘r - l ) , (4.14) where the voltage across it is Vd « Ax • Ex. Although ESN uses field values at three consecutive time steps, t+1, t and t-1, respectively, due to stability problems the average value o f E'x = ( £ '+l + E'~l) / 2 was used to add equation (4.14) to the expression (4.4). 82 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The transcendental equation was solved with Newton’s method. If the diode has to be spanned across several cells,, series/parallel connection schemes can be used, modifying the thermal voltage accordingly. 4.4 Two Port Active Circuit Element Models A two port device, shown in Fig. 4.3, can be described in terms of two equations relating voltages and currents at the input and output ports: 11 I2 V2 V1 Figure 4.3 General reperesentation of a two port device There are several sets of equally valid network parameters that can be used to describe the device. In order to include the 2-port model into ESN, it is convenient to work with its admittance (or Y-matrix) representation: (4.15) 83 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Since the goal is to include the transistor model into the time domain ESN algorithm, equations (4.14) have to be written in the time domain as well. Most semiconductor devices are characterized in terms o f their equivalent circuit models [45,46], which can be analyzed using the Laplace transform technique. In order to be useful in ESN, the appropriate equations for the currents flowing through the device will be found in the Sdomain. Subsequently, the required time domain I=f(V) relationship will be obtained by applying the inverse Laplace transform to the S-domain results. To demonstrate this approach, the implementation of the n-channel JFET (Junction FieldEffect Transistor) into the ESN algorithm will be illustrated. The starting point o f the analysis involves choosing the circuit model for the FET, which describes the behavior of the device for both large and small signals. The equivalent circuit model was chosen following [45,46], and is shown in Fig. 4.4 below: G VG5 D E6 c& Figure 4.4 JFET equivalent transistor model For the triode region ( vDS < vcs - Vp), the drain current ip is given by [45]: 84 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. However, when the transistor is biased into the pinch-off region ( vDS > vGS - Vp ), then the drain current is: ~ ^DSS 1 - 'G S (4.16b) P ' The time domain representation o f the Y-matrix for the network shown in Fig.(4.2) is obtained with the help o f the Laplace transform: dv dt DS '•o (0 = ( c „ + c „ ) % - c , dt drcs dt dt R (4.17a) (4.17b) The next step is to include equations (4.17) into the appropriate components o f Ampere’s Law. Assuming x-directed currents, the combination o f equation (4.17) with (4.2) yields: dHy dz dHz dy dErfs dt ~ ®ExGS + e 0 Cgj dElDS Ax dt (4.18a) (8r - l ) + dE xG S + Cgd £0 -Ax . dt 85 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The quantities Exes and Exds above, were used to denote the lattice cells where the gateto-source and drain-to-source ports were added to the field equations. The system (4.18) was discretized using the same finite-differencing scheme that is employed for ESN. The additional time derivatives were discretized using central differences with respect to time. Once fully discretized, equation (4.18a), for the gate-to-source cell becomes: [(u + g ca Op ,Jgs, k gs) J l + cta Op >Jg, >kg, ))] • C Op J p . )- - u C ^ / (e0Ax) •V £ x( ^ , 7 ^ , k j, ) = = [(u + r O p, jg, , *p ))(l - a „ (/„, j \ s,^p))]* C ' Op, j p , ) (4.19a) - u C d / ( e 0A x ) - V ^ l(ith, j lls,kds) + f-2 2 ‘ Scex Op ’Op»^p ) * Op ’Op >^p ) ^Scex Op ’Jgs »^p ) ‘ exsum - C O p .y 'p - 1/ 2 ,* „ ) + ^ O p .y ’p + 1/ 2 , ^ ) N + 2v K ,y Op»y p ^ p - i / 2 ) - c O p . j g , ^ p + 1/ 2) J where: (4.19b) y = Ax/(c0Ar) , Sc« O p,7p >* p ) = 2 • y •((*, - 1 ) + (C„ + 0. ) / (ff0A x)), y.*) = V/A)/«b(oAx) / (w(2ffr(i, y, *) - 1)) 86 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.19c) (4.19d) (4.19e) Similarly, the drain-to-source cell algorithm is modified as follows: ((” + £ « ( w * . * * ) X 1+ a « ( W * . * < J ) ) - C l(W * .* .fr) + + o C ^ / (e0Ax) • C ‘O'* ./* . ■ ) = = ((u + g ce*0* >7*. ** )X! “ a « 0 * , / * . ^ ))) ’ C 0* ./* » * * ) + + »CgJ/ ( e 0A x ) . V ^ ( i gsJ g, , k gs) + /-2 + 2 •g c„ (/* , j ds,kdt)- V'a (zA, y4 ,£ * ) + 4gfex(/*, yA, ^ ) • Feexsum (4.20a) ~ C O W * - 172,**) + C O W * + 172> 0 ' + 2- - 1 / 2 ) ~ V my{ids, j j s,kds + 1/2) , -2 -V m7 -7 " > Ax where: (4.20b) u = Ax/(c0A /), W W W * ) = 2 • v ( ( e r -1 ) + W > y. *) = / (*0Ax)), ylMo M o t e + 1 7 fl) / (o (2 ^ (r, y, &) - 1)), and V!Zum(ijs’jdS’kjs) = C 0 W W * ) ~ C W W W * ) > (4.20c) (4.20d) (4.20e) Finally, it should be added that, within the ESN lattice, the locations of the gate-to-source and drain-to-source ports are denoted by the coordinates (/ , j ^ , k ^ ) and (id, , j d,, kds). 87 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For microwave JFET’s the capacitances are given by: C ^ l p F and C gd« C^. Since Cgd is small, its effects are typically neglected in the above circuit model. The parasitics of the packaging that encase the die are not taken into account by the above model, but they can be included in ESN independently as microstrips, vias and air-bridges. 4.5 Numerical Results fo r Lumped Element Models A single microstrip geometry (w=1.0 mm, /z=0.6mm, er=4.0), shown in the inset of Fig. 4.5, has been used for most simulations of passive lumped elements. According to the empirical formulas in [42], the parameters of the transmission line are: impedance Zo=56.6£2 and propagation velocity V p = I . 71x10s m/sec. The relative normalized impedance of the line ( Z / Z 0, Z0=120tc Q), in frequency domain, has also been computed using ESN and is presented in Fig. 4.6. The computed quantity is: z ;v( / ) = r ( / ) / / ( / ) , where V(f) and 1(f) are the Fourier transforms of v(t) and i(t): [ istnp \ 2 ^ E x( i , m l 2 , n l 2 ) - dec i=i and fir/rrp+l/2 /( /) = f(/(0) = F jStripRight £ 2) \i= is tr ip - \/2 j =jStripLeft 88 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. respectively. The corresponding characteristic impedance is extracted from Z " ( / ) by multiplying it by Z„ (377 ohm). The numerical value of Zc is 59.2Q, which is in good agreement with the empirical formula [42]. The first problem o f combining lumped elements with ESN, consisted in terminating a microstrip line with its characteristic impedance. The expected result is to have the incoming waves absorbed, with low or very low reflections. 020 0.19 0.18 0.17 0.16 0.15 0.14 0.13 o NJ 0) 0.12 0 11 0.10 (0 a> rr 0.09 v 0.08 0.07 006 005 0.04 h * 0 .6 m m 003 Rel at i ve per m i t i vi t y*4. 0 0.02 0.01 0.00 5.006+08 5.506+09 1.05E+10 1.556+10 F r e q u e n c y (H z) Figure 4.5 Relative impedance of the microstrip in frequency domain 89 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The end resistor has been implemented in three different ways: a) a column o f cells placed at the middle under the microstrip, b) three resistor columns placed at and around the middle under the strip, and c) a distributed resistance placed under the entire width of the strip, as can be seen in Fig. 4.6. ... M ic ro s trip D is trib u te d r e s i s t o r s G r o u n d p la n e a) b) c) Figure 4.6 Distribution of lumped resistor under microstrip 90 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. For the above simulations, the corresponding values of g c x (i,j,k ) have been calculated using the characteristic impedance of Zc = 59.2Q in equation (4.5). The reflection coefficients Su for the three cases of differently distributed terminations are plotted in Fig. 4.7. It can be seen that the least amount of reflections is generated by the uniformly distributed resistor. The levels of reflections are rather low, being comparable to reflections from the general purpose Absorbing Boundary Condition algorithms. 10 0 CQ -10 (/) C -20 o o 0) q= -30 <D a.: ♦ ♦ ♦ A A-*'*' . ■40 1 C d im n Resistor . 3 C d i m s R esistor i 'A T ' -50 5.00E+08 . 5 C d u m s Resistor 3.00E+09 5.50E+09 a00E+O9 1.05E+10 1.30E+10 1.55E+10 1.80E+10 Frequency (Hz) Figure 4.7 Reflections of a lumped resistor terminated microstrip The S21-parameters of two microstrip lines separated by a small gap, that are electrically connected using a series capacitor was also analyzed. The series capacitor should have a 91 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. high insertion loss at low frequencies, which should decrease at higher frequencies. The ESN seems to predict this quite well, as can be seen from the data in Fig. 4.8, where the S2I-parameter is shown. The lumped series capacitor was implemented in ESN by modifying the gce: term according to equation (4.9). •to C *D o Csj C/D C«12.5pF g*0.4m m h*0.6m m R e la tiv e p e r m itivity«4.0 -60 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------5 00E+07 1.05E*09 2 0 5 E -0 9 3Q 5E *09 4.05E+09 5.05E+09 6 0 5 E -0 9 7.05E+09 B.05E+09 9.05E+09 F re q u e n c y (H z) Figure 4.8 Lumped capacitor in series with a microstrip line Another problem o f interest that involves capacitors, is a lumped capacitor connecting a microstrip line to ground, as can be seen in the inset o f Fig. 4.9. This is particularly useful in experimental characterization of microwave capacitors. Measurements were taken for the geometry and dielectric constant values shown in the inset, and are plotted along with the simulation data obtained from ESN. This time, the lumped capacitor was modeled by modifying gcex term in accordance with equation (4.9). 92 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. .10 o— - a - -40 M easured -A -E S N F re q u e n c y (Hz) Figure 4.9 Lumped capacitor between a microstrip line and ground plane Numerical simulations were performed for active circuit elements as well, such as for the diode and high mobility FET transistor. In the case of the diode, the electrical parameters of the components under investigation are shown in the inset o f Fig. 4.10. The diode is placed half-way between the excitation plane and the terminating resistor. The time domain results are plotted in Fig. 4.10. For comparison, the same problem was also simulated using SPICE [51], for which the results are shown in Fig. 4.11. The differences between the two plots can be attributed to the different models used to implement the diode and the transmission lines in ESN and SPICE, respectively. The diode model is more complicated in SPICE, compared to the equation (4.13) used in ESN. On the other hand, the transmission line model in SPICE is much simpler than the full wave representation o f the microstrip in ESN. 93 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R2 59 v U _^w v ^_ 02 OS -10 1 67E-13 837E -11 1 67E-10 2.51E-10 3.34E-10 T im e (se c o n d s ) Figure 4.10 Voltage on a lumped diode in ESN 2.0 0 0- - 2.0 - -6 . 0 - - 10 . 0 - 0 On 0 In 0 2n T IM E [sac] 0 3n 0 .4 n 0 4n Figure 4.11 Spice simulation of the diode circuit As an example o f 2-port device modeling, a microwave amplifier stage, which contains an FET transistor was simulated. The corresponding circuit schematic is shown in Fig. 4.12. 94 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R2 1k —•-AM/---R5 50 = 5 0 , td = 6 0 | IZ=50. td = 60l R3 T 1.5k ~ R4 50 V2 20V Figure 4.12 Schematics of FET microwave amplifier stage The parameters of the JFET are 7 ^ = 37.5 mA, Vp=-2V and Q s=lpF. The AC source has an offset of V0ffsel=-Q.5'V , Vg=0.5V and the frequency is 5GHz. In ESN, the microstrip is specified by w=lmm, h=0A mm and £>=2.73, which yields the characteristic impedance of ZC=5I.9Q and the propagation delay of 4.98 psec/mm. Each segment of the microstrip has a length o f 12mm. The overall structure (two microstrip sections and active device) was discretized into a 15x26x130 lattice, with Ar=0.2mm. The FET was implemented as a two port device directly under the microstrip, as depicted in Fig. 4.13. The waveforms from the SPICE simulation of the above amplifier are shown in Fig. 4.14. The simulation is at steady state, and does take into account the time delays introduced by the transmission lines. 95 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. E e 2 mm Transistor c a s in g £ relative p er m it tiv ity = 2.73 ground plane Figure 4.13 ESN geometry of the JFET problem The time domain waveforms of the ESN simulation of the same transistor circuit are plotted in Fig. 4.15. Compared to the SPICE model, several differences between the results have to be mentioned. ESN is a time domain algorithm, which cannot reproduce exactly step excitations, without generating numerical oscillations due to the nature of the algorithm. Thus the biasing voltage V2, that is applied to the drain must vary smoothly, and for this reason was implemented as a rising one-sided Gaussian pulse. Another cause for discrepancies is due to the steady state excitation of the microstrip, which cannot be absorbed very well by the ABC’s. As a result the terminating resistor of the output 96 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. microstrip will have a value which will differ from the characteristic impedance o f the line. s .e u i eu .6 n e <)ns U (R M :2) 8 .8 n s 1 . 8n t • U ( J 3 :d ) Figure 4.14 Spice simulation of the JFET circuit SV) 5 c (0 tfl -O 2 .0 0 E -0 0 1 OOE+QQ 1.6717E-10 T im e Figure 4.15 ESN simulation of the JFET circuit 97 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Irrespective o f the pointed out sources of discrepancies, the results of the ESN simulation are still close to SPICE generated data. More importantly, despite various limitations o f ESN, the results are very encouraging, as for the first time the JFET transistor sub-cell model was included into a full wave finite difference time domain type algorithm. 98 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5. Results of Microwave Discontinuities Simulations Finally, the new ESN algorithm in this dissertation was applied to analyze several microstrip discontinuities. Some of them are planar, such as the step, the resonant stub and the band-pass filter. Other discontinuities are 3-D structures, such as the via through ground plane. The simulations used as excitation a Gaussian pulse, narrow enough to contain frequencies up to 20 GHz. The incident, transmitted and reflected time domain signals were subsequently transformed into the frequency domain and the S-parameters were computed. All results have been compared against data obtained by other means, such as frequency domain methods, FD-TD simulations, or direct measurements. 5.1 Microstrip step The first microstrip discontinuity, which was simulated with ESN and was used for validation, is the step in width. The corresponding geometry is shown in the inset o f Fig. 5.1. The structure is divided into a l x m x n = 30x30x80 lattice, with A = Ax = Ay = Az = 0.423/ww. The time step is taken to be At = 0.75 x 10"12 sec, for the total run time of 1333 time steps, to obtain the frequency resolution of 1 GHz. The 99 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. numerical results, shown in Fig. 5.1, are compared to data available in the literature [18,19] with a very good agreement seen between them. Q9 2 | I? §5 Q8 < 2W, .(*>&. [W 19) .SM I Q. i -FW[I9| CO 05 O (0 *©3o Q4 I>Q (0 3 S Q2i .SM I 2W, S,, 2a Q1 5 6 Frequency (GHz) Figure 5.1 S-parameters of step in width discontinuity 5.2 Microstrip resonant stub In the second validation example, the S-parameters are computed for a symmetric stubtype structure, which is shown in the inset of Fig. 5.2. The S,,-parameter of the double stub structure is shown in Fig. 5.2, from which it can be seen that the data obtained with ESN is similar to data published in [13]. Both sets of results, as well as the measured data, predict resonant frequencies at about 8 and 19 GHz. The structure in Fig. 5.2 is discretized into a / x m x « = 1 8 x 6 1 x 8 0 lattice, with A = 0.4mm, for the total run time of /=9000 time steps, producing a frequency increment of A / = 250M H z. The S21- 100 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. parameters show a good agreement near the first resonance. For the second resonance S2, is -23dB for ESN, and -lOdB for FD-TD, compared to -32 dB measured data. oo -10 efi -15 -20 20.32mm 5.65mm 2.09mm 0 6' -35 •40 - o —Measured —A—Ref. [131 -O -SN M £ A -45 7 8 0.794 m m 9 10 11 12 13 14 15 16 17 18 19 F r e q u e n c y (G H z) Figure 5.2 S2I-parameter of a stub 5 .5 3-D Via through a Ground Plane A via through a ground plane is one of the most common discontinuities encountered in the manufacturing o f printed circuit boards (PCBs). In the industry, for frequencies below 100MHz, the effects o f a via on the signal integrity is considered to be similar to the effect of a 5pF capacitor. ESN was used to ascertain if this approximation holds for frequencies in the GHz range as well. Two via diameters of 0.7mm and 1.5mm, respectively, have been simulated, and the geometry of the via is shown in Fig. 5.3. It should be noted that since the discretization of ESN is based on the rectangular 101 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. coordinates, the circular shape of the via, pads, and the hole through the ground plane were modeled using a stair-case approximation (see Fig. 5.3). 3.2 mm 0.7 or 1.5 mm, 4 mm A 3.4 CD I T = 3.4 metal plane (a) (b) Figure 5.3 Geometry of via through ground plane Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. * 1 -o— -15 -0-S11-ESN -a-S11-ftteaiBd S11-FDTD -X-S21-ESN -x-S21-ftteared; S21-HJIL) CO-25 -30 -36 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Frequency (G-fe) Figure S.4 S-parameters of 0.7 mm via The computed S-parameters of the two vias, shown in Figs. 5.4 and 5.5, respectively, were compared against measured data and FD-TD generated data available in [10]. Note that the agreement between the measured data, FD-TD, and ESN is very good. The minor discrepancies in measured and ESN data can be attributed to the differences in the actual geometry and the stair-case approximations used in ESN. On the other hand, the differences between the ESN and FD-TD data can be related to different discretization methods, such as to the use of variable and uniform grids in [10] and ESN, respectively. To generate the data with ESN, the geometry was uniformly discretized into 32x40x184 cell lattice, with A=0.2 mm. 103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Comparing the plots in Fig 5.5 for the via, and Fig. 4.8 for the series capacitor, respectively, it can be seen that the shapes of the S2I-parameters curves are similar. It can thus be suggested that a large diameter via (« 1.5mm) can be approximated with a series capacitor, which electrically connects the two microstrip segments. For smaller diameter vias the (<0.7mm) the S2I-parameter curve in Fig. 5.4 suggests that a more complex circuit approximation is needed. -0 -S 1 1 -E S N -D -S 1 1 -M e a s u re d -20 S11-FDTD -25 - x — S 21-M easured - x — S21-E SN - 0 -S 2 1 -F D T D -30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Frequency(GHz) Figure 5.5 S-parameters of 1.5 mm via 5.4 Band pass filter Highly resonant structures, such as side coupled microstrip band pass filter, was also analyzed with ESN. The geometry o f the filter is shown in Fig. 5.6. 104 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12.7mm 6.35 mm 1.272 mm 1.272 mm 6.35 mm 1.272 mm Relative permittivity= 10.0 Figure 5.6 Geometry of band-pass filter The structure was discretized using a 20x80x100 lattice, with Ax=0.212mm. The filter was excited with a Gaussian pulse s(k) = e"0-5^ l,a) , with t = 90At and cr = 30A/ . The simulation was allowed to run for 9000 time steps. As in the case of the stub, the runtime o f the algorithm must be long, as the energy in the resonant structure is accumulated under the microstrip, and is slowly released, at specific resonant frequencies. The S2Iparameter o f the filter is shown in Fig. 5.7, where the simulation results o f ESN are compared against the previously published simulation using SNM [20], as well as measured data. 105 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. •15 • *o •25 •35 -45 •55 S21-ESN - - S21-SNM —S2l-measurad -65 0.1 1.1 2.1 3.1 4.1 5.1 6.1 7.1 8.1 9.1 1 0.1 11.1 12.1 13.1 F r e q u e n c y (G H z) Figure 5.7 S-parameters of microstrip filter The resonant frequencies are in very good agreement. The ESN simulation yields a lower maximum value for the S2l-parameter because there is still energy accumulated under the microstrip, which has not yet been released. The results would have been closer to the measured values if the simulation was allowed to run for a much larger number of timesteps (>20,000), which is a very long simulation (30 hours), compared to the usual 2000 time steps needed for non-resonant structures. For practical purposes, a difference of about -5dB can be considered as acceptable in engineering applications and design. 106 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6. Conclusions and Recommendations 6.1 Conclusions As it was stated in the introduction, full-wave simulations using time domain methods are increasingly being used to analyze microwave discontinuities and more complex circuits in which wave-guiding and lumped elements are found together. One o f the goals o f this work was to develop a circuit/transmission line based algorithm, which makes the understanding o f time domain simulations easier for non-specialists in electromagnetics. It was essential to verify the correctness of the model with respect to Maxwell’s equations. In addition, to make the new time domain solver more versatile, it was necessary to incorporate subcell models of lumped circuit elements into it. Last, but not the least important, it was required to develop an appropriate method for implementing the ABC’s for open lattice truncation. The new field solver, which was called ESN, was developed starting with the circuit model of SNM. However, the new ESN algorithm was generalized and made much more computationally efficient than SNM. All variables in the algorithm, related to material parameters, were defined according to the consistency requirements with Maxwell’s equations, thereby explaining the heuristic assumptions which were present in SNM. A numerical optimization method was developed for implementing PML absorbing 107 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. boundary conditions, aimed especially at microstrip line structures. It should be noted that the developed PML method is algorithm independent, which makes it useful not only for ESN, but for FD-TD as well. Appropriate subcell models for active and passive circuit elements were obtained for use in ESN and a new method of implementing two-port devices has been proposed. The numerical analysis of different discontinuities has shown that ESN yields results in the same precision range as FD-TD, and very close to measured data as well. Very encouraging are the results obtained for the FET model, which offer the possibility of analyzing microwave integrated circuit building blocks in the near future. 6.2 Recommendations for Future Research Although the general goals set at the beginning of the research were reached, a few topics suggested below are worthy of further investigation. One of these is the inclusion of anisotropic materials that are characterized with arbitrary forms of permittivity and permeability tensors. The need for this arises since most material substrates used in microwave and millimeter-wave applications display such properties. It should be added that although it was never explicitly mentioned in the dissertation, the present ESN algorithm incorporates biaxial (diagonal permittivity tensor) electric anisotropy. 108 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The computational speed of any time domain algorithm is always an issue for which improvements must be sought at any time. In its present form, ESN requires ten times less run-time, and about half the memory of SNM. However, its CPU performance can be further improved, if the algorithm adapted to take advantage o f multi-processor machines. Currently, even desk-top operating systems such as Windows NT® are able to utilize up to 6 parallel microprocessors. This implies that writing a 6 synchronous-thread program would increase the speed o f ESN substantially. Another area in which advances were made and described in this dissertation is the open boundary truncation. The new method for computing the PML “loss” coefficients was shown to improve the performance o f ESN for microstrip type problems significantly. In order to simulate radiation problems, it is desirable to study the behavior o f the PML ABC’s at different angles o f incidence as well. Advances in this area would expand the applicability of ESN to microstrip antenna problems as well. Finally, further work is warranted in the development and implementation of lumped circuit element models. The JFET transistor model is very sensitive to the discontinuities in the excitation applied to it, hence needing further refinement. At the same time, it is desirable to build a library of circuit models, similar to that of SPICE, which allows non specialists to utilize the algorithm for the design and simulation of integrated microwave circuits. 109 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Bibliography 1. T. Itoh, R. Mittra and D. Ward, “A method for computing edge capacitance of finite and semi-finite microstrip lines”, IEEE Trans. Microwave Theory Tech., vol. MTT20, pp. 847-849, Dec. 1972. 2. C. Gupta and A. Gopinath, “Equivalent circuit capacitance of microstrip step change in width”, IEEE Trans. Microwave Theory Tech., vol. MTT-25, pp. 819822, Oct. 1977. 3. G. Compra and R. Mehran, “Planar waveguide model for calculating microstrip components”, Electron. Lett., vol. 11, pp. 459-460, Sept. 1975. 4. R. H. Jansen, “The spectral domain approach for microwave integrated circuits”, IEEE Trans. Microwave Theory Tech., vol. MTT-33, pp. 1043-1056, Oct. 1985. 5. T. Itoh and R. Mittra, “A technique for computing dispersion characteristics o f shielded microstrip lines”, IEEE Trans. Microwave Theory Tech., vol. 22, no. 10, pp. 896-898, Oct. 1974. 6. T. Becks and I. Wolff, “Analysis of 3-D metallization structures by a full wave spectral domain technique”, IEEE Trans. Microwave Theory Tech., vol. 40, no. 12, pp. 2219-2227, Dec. 1992. 7. S. Nam, H. Ling and T. Itoh, “Characterization of microstrip line and its discontinuities using the time-domain method of lines”, IEEE Trans. Microwave Theory Tech., vol. 37, no. 12, pp. 2051-2057, Dec. 1989. 110 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 8. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media”, IEEE Trans. Antennas Propag., vol. AP14, no. 3, pp. 302-307, May 1966. 9. X. Zhang and K. K. Mei, “Time-domain finite difference approach to the calculation of the frequency-dependent characteristics of microstrip discontinuities”, IEEE Trans. Microwave Theory Tech., vol. 36, no. 12, pp. 17751787, Dec. 1988. 10. S. Maeda and T. Kashiwa, “Full wave analysis o f propagation characteristics of a through hole, using the finite difference time domain method”, IEEE Trans. Microwave Theory Tech., vol. 39, no. 12, pp. 2154-2159, Dec. 1991. 11. H. Jin and R. Vahldieck, “The frequency-domain transmission line matrix method a new concept”, IEEE Trans. Microwave Theory Tech., vol. 40, no. 12, pp. 22072218, Dec. 1992. 12. G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations”, IEEE Trans. Electromagnetic Compatibility, vol. EMC-23, no. 4, pp. 377-382, Nov. 1981. 13. R.L. Higdon, “Numerical absorbing boundary conditions for the wave equation”, Math. Computation, vol. 49, no. 179, pp. 61-91, July 1987. 14. C. Eswarappa and W. J. R. Hoefer, “One-way equation absorbing boundary conditions for 3-D TLM analysis of planar and quasiplanar structures”, IEEE Trans. Microwave Theory Tech., vol. 42, no. 9, pp. 1669-1677, Sept. 1994. Ill Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 15. X. P. Lin and K. Naishdaham, “A simple technique for minimization o f ABC induced error in the FD-TD analysis of microstrip discontinuities”, IEEE Microwave Guided Wave Lett., vol. 4, no. 12, pp. 402-404, Dec. 1994. 16. O. Ramahi, “Complementary operators: a method to annihilate artificial reflections from the truncation of the computational domain in the solution o f partial differential equations”, IEEE Trans. Antennas Propag., vol. 43, no. 7, pp. 697-704, July 1995. 17. D. Katz, E. Thiele and A. Taflove, “ Validation and extension to three dimensions o f the Berenger PML absorbing boundary conditions for FD-TD meshes”, IEEE Microwave Guided Lett., vol. 4, no. 8, pp. 268-271, Aug. 1995. 18. A. Bahr, A. Laurel and I. Wolff, “Application o f the PML absorbing boundary conditions to the FD-TD analysis of microwave circuits”, 1995 IEEE MTT-S Symp. Digest, vol. 1, pp. 27-30. 19. N. Yoshida and I. Fukai, “Transient analysis of a stripline having a comer in threedimensional space”, IEEE Trans. Microwave Theory Tech., vol. MTT-32, no. 5, pp. 491-498, May 1984. 20. T. Shibata, T. Hayashi and T. Kimura, “Analysis o f microstrip circuits using three dimensional full wave electromagnetic field analysis in time domain”, IEEE Trans. Microwave Theory Tech., vol. 36, no. 6, pp. 1029-1035, June 1988. 21. S. Sasaki, R. Konno and H. Tomimuro, “3-D electromagnetic field analysis of interconnections in copper-polymide multichip modules”, IEEE Trans. Comp. Hybrids Manuf. Technol., vol. 14, pp. 755-760, Dec. 1991. 112 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 22. M. J. Tsai, T. S. Horag and N. G. Alexopoulos, “Via hole, bond wire and shorting pin modeling for multi-layered circuits”, 1994 IEEE MTT-Symp Digest, vol. 3, pp. 1777-1780. 23. J. Bergeron, ‘'‘Water hammer in hydraulics and wave surges in electricity”, (translation), New York: Wiley & Sons, 1961. 24. K. Wu, M. Yu, R. Vahldieck, “Rigorous analysis of 3-D planar circuit discontinuities using the space-spectral domain approach (SSDA)”, IEEE Trans. Microwave Theory Tech., vol. 40, no. 7, pp. 1475-1483, July 1992. 25. R. Luebbers, “FD-TD for antennas and scattering”, Short Course Notes from the 10th Annual ACES Conference, Monterey, CA, 1994. 26. N. H. L. Koster and R. H. Jansen, “ The microstrip step discontinuity : A revised description”, IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 213-223, Feb. 1986. 27. Y. Chen, and B. Beker, “Study of microstrip step discontinuities on bianisotropic substrates using the method o f lines and transverse resonance technique”, IEEE Trans. Microwave Theory Tech., vol. 42, no. 10, pp. 1945-1950, Oct. 1994. 28. B. Katehi and N. G. Alexopoulos, “Frequency dependent characteristics of microstrip discontinuities in millimeter wave integrated circuits”, IEEE Trans. Microwave Theory Tech., vol. MTT-33, pp. 1029-1035, Oct. 1985. 29. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves”, J. Comp. Phys., vol. 114, no. 2, pp. 185-200, Oct. 1994. 113 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 30. J. C. Veihl and R. Mittra, “An efficient implementation of Berenger’s perfectly matched layer (PML) for finite-difference time-domain mesh truncation”, IEEE Microwave Guided Wave Letters, vol. 6, no. 2, pp. 94-96, Feb. 1996. 31. L. Zhao and A. Cangellaris, “ A general approach for the development o f unsplitfield time-domain implementations of perfectly matched layers for FD-TD grid truncation”, IEEE Microwave and Guided Wave Letters, vol. 6, no. 5, pp. 209-211, 1996. 32. W. 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 Tech., vol. 40, pp. 724-730, Apr. 1992. 33. M. Picket-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 Tech., vol. 42.,no. 8, ppl514-1523, Aug. 1994. 34. Y. Tsuei, A. C. Cangellaris and J. L. Prince, “Rigorous electromagnetic modeling of chip-to-package (first-level) interconnections”, IEEE Trans. Components, Hybrids and Manufacturing Tech., vol. 16, pp. 876-883, Dec. 1993. 35. D. Bica and B. Beker, “Analysis o f microstrip discontinuities using the spatial network method with absorbing boundary conditions”, IEEE Trans. Microwave Theory Tech., vol. 44, no.7, pp. 1157-1161, July 1996. 36. D. Bica and B. Beker, “ Enhanced spatial network method for the analysis o f open microstrip discontinuities”, to be published in IEEE Trans. Microwave Theory Tech., vol. 45, July 1997. 114 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 37. T. Shibata, T. Hayashi and T. Kimura, “ Analysis of microstrip circuits using threedimensional full-wave analysis in the time-domain”, IEEE Trans. Microwave Theory Tech., vol.36, no. 6, pp.1064-1070, June 1988. 38. J. Fang and Z. Wu, “Generalized perfectly matched layer for the absorption of propagating and evanescent waves in lossless and lossy media”, IEEE Trans. Microwave Theory Tech., vol. 44, no. 12, pp. 2216-2222, Dec. 1996. 39. A. Nedler and R. Mead, “A simplex method for function minimization”, Computer Journal, vol. 7, pp. 308-313,1965. 40. J. Cain, “Numerical analysis of microstructure discontinuities of millimeter wave integrated circuits printed on anisotropic substrates”, Ph.D. Dissertation, University of South Carolina, Columbia, SC, Dec. 1995. 41. L. Lapidus and G. H. Pinder, Numerical solutions o f partial differential equations in science and engineering, New York: John Wiley&Sons, 1982. 42. A. Balanis, Advanced engineering electromagnetics, New York: Wiley, 1989. 43. E. Issacson and H. Bishop Keller, Analysis o f numerical methods, Dover Publications, 1994. 44. A. Taflove, The finite-difference time-domain method fo r numerical modeling o f electromagnetic wave interactions with arbitrary structures, PIERS-2 - Progress in Electromagnetic Research, pp. 287-369, Elesevier 1990. 45. A. Sedra and K. Smith, Microelectronic circuits, Holt, Rinehart and Wilson, 1987. 46. H. Newdeck, Electronic circuit analysis and design, Houghton Mifflin Company, Boston, 1976. 115 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 47. B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical simulation of waves”, Math. Comp., vol. 31, pp. 629-651,1977. 48. K. K. Mei and J. Fang, “Superabsorption - A method to improve absorbing boundary conditions”, IEEE Tran. Antennas Propag., vol. AP-40, pp. 1001-1010, Sept. 1992. 49. Z. P. Liao, H. L. Wong, G. P. Yong and Y. F. Yuan, “A transmitting boundary for transient wave analysis”, Scientia Sinica, vol. 28, no. 19, pp. 1063-1076, Oct. 1984. 50. Z. Wu and J. Fang, “ Numerical implementation and performance of perfectly matched layer boundary condition for waveguide structures”, IEEE Tran. Microwave Theory Tech., vol. 43, no. 12, pp. 2672-2683, Dec. 1995. 51. Microsim PSpice EDA Software, 20 Fairbanks, Irvine, California, 92718 116 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

1/--страниц