close

Вход

Забыли?

вход по аккаунту

?

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.
Документ
Категория
Без категории
Просмотров
0
Размер файла
3 824 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа