ADAPTIVE NOVEL TIME-DOMAIN TECHNIQUES FOR THE MODELING AND DESIGN OF COMPLEX RF AND WIRELESS STRUCTURES
A Dissertation Presented to The Academic Faculty
By Nathan Bushyager
In Partial Fulfillment Of the Requirements for the Degree Doctor of Philosophy in Electrical and Computer Engineering
The Georgia Institute of Technology November 2004
Copyright 2004 Nathan Bushyager
NOVEL ADAPTIVE TIME-DOMAIN TECHNIQUES FOR THE MODELING AND DESIGN OF COMPLEX RF AND WIRELESS STRUCTURES
Approved by:
Dr. Manos M. Tentzeris, Advisor Dr. Andrew Peterson Dr. Joy Laskar Dr. Ioannis Papapolymerou Dr. Fotis Sotiropoulos Date Approved: Nov. 17, 2004
ACKNOWLEDGEMENTS
First, I would like to thank my wife, Misha, who has been the foundation that supported me through most of my graduate career. Graduate school requires many sacrifices, and while this work is mine, the effort is ours. When I needed help she was there, and when I needed a smack on the head, she gave it to me. Thanks for taking this journey with me. I would like to thank my mother, Charlene, for all of her support. Beyond making sure that I survived childhood, she fostered my interest in the sciences and taught me from early on the importance of knowledge and independence. I would like to thank my father, John, for his do-it-yourself attitude. I learned more than he knows from handing him his tools as a child. To my brother Matt, thanks for being different from me (and sometimes annoyingly similar). Sometimes. Sometimes diversity of opinion is a good thing.
I would also like to thank my grandmother, Grace; theres a strong
independent streak in my family, and I have a suspicion that she has a lot to do with it. I would like to thank my advisor, Manos Tentzeris. Without Manos, I quite
literally would not have made it through grad school. I couldnt afford it (also, Id like to thank the Georgia Electronic Design Center and the Packaging Research Center for funding our research). In addition, I learned a few things along the way. Some are in this thesis, and some are better not written down. In the academic world, I would specifically like to thank Joy Laskar, John Papapolymerou, Andrew Peterson, and Fotis Sotiropoulos. First, you let me put your names on the approval page. Individually, Dr. Laskar, has taught me leadership and has given quality movie review advice, Dr. Papapolymerou has shown me how to survive eating in foreign countries (and taught me a thing or two about MEMS), and Dr. Peterson is one of the best professors Ive ever had.
iii
I would like to thank all of the members of the ATHENA group, specifically Brian McGarvey and Edan Dalton, who have been with me practically since the beginning. They also helped get me through the prelim. I would also like to thank members of the MircTECH and MAG groups, several of which have worked on various research topics with me. Finally, Id like to thank Gerald DeJean, who helped me with some of the antenna designs, and knows a lot more about antennas and sports than I do.
iv
TABLE OF CONTENTS
ACKNOWLEDGEMENTS............................................................................................... iii LIST OF FIGURES .......................................................................................................... vii LIST OF SYMBOLS AND ABBREVIATIONS .............................................................. ix SUMMARY........................................................................................................................ x CHAPTER 1 INTRODUCTION ....................................................................................... 1 CHAPTER 2 BACKGROUND ......................................................................................... 6 2.1 2.2 MRTD BASICS ................................................................................................ 6 2.1.1 Wavelet Overview.................................................................................... 6 2.1.2 Method of Moments Overview .............................................................. 18 MRTD Update Scheme................................................................................... 20 2.2.1 Time Localization................................................................................... 20 2.2.2 Space Localization ................................................................................. 23 2.2.3 Media Discretization .............................................................................. 27 2.2.4 Numerical Stability................................................................................. 28 2.2.5 Numerical Dispersion............................................................................. 32 IMAGE THEORY FOR PEC MODELING................................................... 35 EXISTING PML TECHNIQUES IN MRTD ................................................. 37 FDTD .............................................................................................................. 38
2.3 2.4 2.5
CHAPTER 3 MRTD COMPOSITE-CELL MODELING .............................................. 41 3.1 3.2 3.3 3.3 GENERAL, EXPLICIT, SUBCELL FIELD MODIFICATION.................... 46 PROPERTIES OF HAAR-WAVELET DISCRETIZATIONS...................... 50 HAAR SUBCELL PEC APPLICATION....................................................... 58 GENERAL SUBCELL EFFECTS IN HAAR-MRTD: COMPOSITE CELLS .................................................................................... 61
CHAPTER 4 FULL WAVE HAAR-MRTD WITH COMPOSITE-CELL MODLELING, UPML, VARIABLE GRIDDING, AND TIME/SPACE ADAPTIVE GRIDDING........ 70 4.1 4.2 4.3 4.4 ARBITRARY WAVELET RESOLUTION UPML....................................... 71 NON-UNIFORM GRID IN MRTD................................................................ 75 MRTD GRID EXCITATION ......................................................................... 82 TIME/SPACE ADAPTIVE GRIDDING ....................................................... 84
CHAPTER 5 MRTD SIMULATION EXAMPLES........................................................ 89 5.1 PML ABSORPTION USING MICROSTRIP LINE...................................... 92
v
5.2 LUMPED ELEMENT VERIFICATION: RESISTOR TERMINATED MICROSTRIP LINE ...................................................................................... 94 5.3 MICROSTRIP PATCH ANTENNA .............................................................. 96 5.4 DUAL MICROSTRIP PATCH ANTENNAS.............................................. 106 CHAPTER 6 CONCLUSION........................................................................................ 111 REFERENCES ............................................................................................................... 114
vi
LIST OF FIGURES
Figure 1: Figure 2: Figure 3: Figure 4: Figure 5: Figure 6: Figure 7: Figure 8: Figure 9: Figure 10: Figure 11: Figure 12: Figure 13: Figure 14: Figure 15:
Haar scaling function (top) and mother wavelet (bottom)............................. 9 Haar wavelets for resolution 1 and 2 ............................................................ 10 Two-dimensional Haar coefficients for rmax=0 ............................................. 14 2-D Haar cell, rmax=0, equivalent grid points................................................ 17 Equivalent grid points for a 2-D MRTD cell, rmax=0.................................... 18 Time localization using pulse functions ....................................................... 21 Image theory in one dimension using biorthogonal scaling functions ......... 37 Yee-Cell in 3-D............................................................................................. 39 Time localization (pulse basis functions) for lossy case............................... 43 FDTD grid intersected by PEC ................................................................. 48 Haar scaling and wavelet functions in 2-D, rmax=1................................... 52 (a) Haar scaling and rmax=0 wavelet, (b) sample function that can be represented using (a) ................................................................................. 54 Haar example with (a) scaling through rmax=1 wavelets, (b) an example function that can be represented with these functions .............................. 56 2-D FDTD non-uniform grid example...................................................... 77 Offset between electric and magnetic fields in MRTD (a) fixed grid (b) non-uniform grid (implemented incorrectly) (c) non-uniform grid (implemented correctly)............................................................................ 80 Excitation in MRTD, rmax=0 cells demonstrated by alternate shading .... 83 CPW excitation, demonstrating subcell excitation ................................... 84 Haar MRTD cell (2-D), showing FDTD grid points, rmax=0 .................... 86 Adaptive grid example.............................................................................. 87 Reflection from PML................................................................................ 93 Characteristic impedance of microstrip line ............................................. 95
Figure 16: Figure 17: Figure 18: Figure 19: Figure 20: Figure 21:
vii
Figure 22: Figure 23: Figure 24: Figure 25: Figure 26: Figure 27: Figure 28: Figure 29: Figure 30: Figure 31: Figure 32: Figure 33: Figure 34: Figure 35: Figure 36: Figure 37: Figure 38: Figure 39: Figure 40:
Time domain reflection from resistor ....................................................... 95 Frequency domain reflection from resistor............................................... 96 Microstrip fed patch antenna .................................................................. 100 Time domain response of antenna, FDTD and MRTD rmax=1 ............... 101 S11 of patch antenna, MRTD/FDTD comparison ................................... 101 Close-up of differences between MRTD/FDTD, time domain .............. 102 FDTD/MRTD comparison, constant r ................................................... 102 Comparison of rmax=1 MRTD with surround grid fixed at r=-1 ............. 103 Comparison of rmax=1 MRTD with surround grid fixed at r=0............... 103 Comparison of fixed grid MRTD with adaptive grid tr=0.1 ta=0.01 ...... 104 Comparison of fixed grid MRTD with adaptive grid tr=0.01 ta=0.001 .. 104 Comparison of fixed grid MRTD with adaptive grid -4 -5 tr=1x10 ta=1x10 ................................................................................. 105 Number of gridpoints used in simulation vs. time, tr=1x10-4 ta=1x10-5 105 Dual microstrip patch antenna (shaded area is PML)............................. 108 S11 for dual antennas ............................................................................... 108 S22 for dual antennas ............................................................................... 109 S21 for dual antennas ............................................................................... 109 Adaptive grid for initial r=1 everywhere ................................................ 110 Adaptive grid with resolution surrounding antenna fixed at r=-1........... 110
viii
LIST OF SYMBOLS AND ABBREVIATIONS
EM FDTD MEMS MMIC MRTD LTCC PEC PMC PML RF SOP UPML
Electromagnetic Finite-Difference Time-Domain Microelectromechanical Systems Monolithic Microwave Integrated Circuit Multiresolution Time-Domain Low Temperature Cofired Ceramic Perfect Electrical Conductor Perfect Magnetic Conductor Perfectly Matched Layer Radio Frequency System on Package Uniaxial Perfectly Matched Layer
ix
SUMMARY
A method is presented that allows the use of multiresolution principles in a time domain electromagnetic modeling technique that is applicable to general structures. Specifically, methods are presented that are compatible with the multiresolution timedomain (MRTD) technique using Haar basis functions that allow the modeling of general structures without limiting the cell size to the features of the modeled structure. Existing Haar techniques require that cells be homogenous in regard to PECs and other localized effects (with the exception that r can vary throughout the cell). The techniques that are presented here allow the modeling of these structures using a subcell technique that permits the modeling of these effects at individual equivalent grid points. This is
accomplished by transforming the application of the effects at individual points in the grid into the wavelet domain. There are several other contributions that are provided in this work. First, the MRTD technique is derived for a general wavelet basis using a relatively compact vector notation that both makes the technique easier to understand and allows the differences and similarities between different MRTD schemes more apparent. Second, techniques such as the uniaxial perfectly matched layer (UPML) for arbitrary wavelet resolution and non-uniform gridding are presented for the first time. Using these techniques, any
structure that can be simulated in Yee-FDTD can be modeled with Haar-MRTD. For the first time, results for the use of a time-and-space-adaptive grid in an MRTD simulation are presented.
x
CHAPTER 1 INTRODUCTION
The field of RF and microwave design is growing at an astonishing rate due the confluence of several factors. Chief among these factors are the increased demand for RF and microwave consumer devices (cellular telephony and wireless data systems) and the increasing speed of digital devices that has led to design limitations that previously only applied to traditional RF (radar and communications) circuits. There are also
several quickly growing commercial RF applications such as automotive radar and RFID. As RF devices become more predominant in the consumer market, the design time for each generation decreases while performance demands increase. Design turnaround and performance gains traditional to the semiconductor device market are being expected of RF circuits. This in turn leads to increased demands on the RF designer and all
associated tools, such as EM simulators. Modern RF devices are built on a variety of technologies for a wide array of functionality. In the attempt to reduce size and cost, multilayer substrates are
commonplace. The constraints of each substrate are different, and the devices that can be fabricated with them often have no, or highly inaccurate, theoretical or empirical models. Design is usually performed using a top-down methodology, where the system is designed at the conceptual level, and then details of the actual system are added as the design progresses. At the bottom of this design process, the actual physical layouts that represent the design blocks must be created. The characterization of these devices
usually requires a full-wave electromagnetic simulator. Full wave simulators are often required to characterize effects that cannot be predicted, or properly accounted for, at higher levels of the design process. Examples of these effects are parasitic coupling, substrate modes, radiation, and package interference.
1
These effects are often unique to the exact layout of a device and are too complicated to be treated theoretically. Any simulator that models the complete physics of
electromagnetic interaction can be used to model these devices, however, time domain techniques are popular and particularly well suited to these devices. Time-domain techniques, as the name implies, determine the electromagnetic fields in a structure in response to an incident condition. They are contrasted to frequencydomain methods, which determine the response to a harmonic source. There are several advantages to each class of methods, and the requirements of a particular problem determine which is best. The characteristic that is most often cited in recommending time-domain methods for the simulation of microwave devices is that a broadband response can be determined from a single simulation. In addition, time domain
simulators are usually derived using a differential form of Maxwells equations, relating fields at points and allowing easy discretization of complex structures; they do not require the calculation of Greens functions. In addition nonlinear effects can be
simulated easily, as the field strength is measured as a function of time and these effects can be directly applied. The two main drawbacks of time-domain techniques are The speed of these techniques, while
representing dispersive media and speed.
comparable to other full wave techniques, is not nearly fast enough for modern microwave design. The most mature and widely used time-domain simulation technique is finitedifference time-domain (FDTD) [1]. The Yee-FDTD scheme is one of the oldest, and one of the simplest, time domain EM simulation techniques. It has remained popular due to this simplicity; not only is it relatively easy to apply, but it is also highly flexible. Due to its flexibility, the method has been applied to practically every type of electromagnetic analysis, from radar-cross-section investigation to optical characterization. A number of compatible techniques have been developed to model specific effects in the method.
2
When developing other methods, the Yee-FDTD scheme is often used as a benchmark; most other techniques are designed to be faster. While a number of techniques have been developed that are at least as fast as FDTD for the same accuracy, a technique that will permit modeling at the same level of accuracy and flexibility as FDTD while offering severe (several orders of magnitude in execution time and memory efficiency) performance improvement is not likely to be developed. One promising method of decreasing design time is to identify parts of the circuit which require full-wave simulation, and characterizing the remainder of the circuit with less computationally intensive techniques. However, as integration increases, the areas of the circuit that require full-wave simulation will increase. Thus, it is necessary to use the most efficient full-wave technique possible. One technique that has been suggested is multiresolution time-domain (MRTD) [2]. MRTD uses a wavelet based discretization to represent the EM fields. The
wavelets allow the resolution of the simulation to be changed as a function of both time and space. If the resolution is changed to respond to a propagating waveform, the number of operations that must be performed to characterize a structure can be minimized. While these advantages are suggested by the MRTD method, many of the details of their application have never been shown. While several papers have been written documenting the MRTD algorithm, most details have never been published. First, a general discretization of the perfectly-matched layer (PML) has never been presented. Existing publications show the PML for a single level of wavelet resolution. Other papers have argued that the PML cannot be efficiently represented in MRTD. A general framework for understanding and implementing
MRTD is presented here which allows the efficient modeling of the PML for any level of wavelet resolution. Another required element for any realistic simulation is non-uniform gridding. This method allows the size of the cells to vary as a function of position. To accurately
3
represent most structures in any method, non-uniform gridding is required; a single cell size that can represent all features of a structure is wasteful, and often not possible. The application of non-uniform gridding in MRTD requires a relatively straightforward but non-trivial change to currently published MRTD algorithms. algorithm is presented here. One of the largest difficulties in applying the MRTD method is representing electrically small structures in the MRTD grid. Because of its multiresolution nature, the domain of each basis function is large when compared other methods. Each discrete electric/magnetic field point is represented as a sum of several coefficients. Thus, the fields at one point cannot be modified without modifying the fields at surrounding points. Existing methods for modeling discontinuities with these basis functions require the use of image theory to decouple the fields across the boundary. These methods can be difficult and time-consuming to apply. For this research a less complex but The variable gridding
mathematically rigorous method was developed to allow the modeling of intracell structures in the Haar-MRTD method. This method allows the modeling of any structure that can be represented in the Yee-FDTD technique, while keeping the advantages of multiresolution analysis. In addition, it offers a bridge between a pointwise field
representation, such as that in FDTD, to the wavelet based field representation of MRTD. Using this technique, point based techniques that have been developed for FDTD can be applied to MRTD. This is demonstrated for discrete lumped elements. The chief contribution of this thesis is the intracell Haar-MRTD modeling method, which is termed composite-cell MRTD. To accurately demonstrate and test this The other MRTD features,
technique, a general 3D MRTD code was developed.
mentioned above, were developed as a consequence of creating this code. In addition, to fully evaluate the usefulness of this code, an adaptive resolution algorithm was developed and applied to the code.
4
The next chapter presents the background necessary to understand the techniques presented in this paper and their comparison to existing methods. As such, a full
derivation of the general MRTD technique is presented, as well as a brief overview of the FDTD technique. The remainder of this document is dedicated to the presentation of the Haar-MRTD algorithm that was developed for this investigation. This includes the arbitrary resolution PML, variable gridding, and composite cell techniques. A discussion is also presented regarding the application of these techniques to other basis functions. This paper concludes with several examples of the application of the MRTD technique. It is shown in this document that the MRTD technique can be applied to any structure and that its time-and-space-adaptive grid can be used to greatly reduce execution time. A framework is presented that allows the reader to construct an efficient and quick MRTD simulator, and several original techniques are presented.
5
CHAPTER 2 BACKGROUND
This thesis presents techniques that can be used to model complex microwave structures in MRTD. In order to properly present the MRTD techniques that were developed as a part of this investigation, a background of the MRTD method, from its genesis to current research, is presented. In addition, as the MRTD technique is
presented as an adaptive alternative to FDTD, and several MRTD examples will be contrasted to FDTD, an brief overview of the FDTD method is presented. This chapter begins with a discussion of MRTD basics, including wavelet basis functions and the method of moments, followed by a discussion of two specific areas, PML and PEC modeling, that were a particular focus of this effort. After these topics are addressed, an overview of FDTD, focusing on its similarities and differences with MRTD, is presented.
2.1
MRTD BASICS
2.1.1 WAVELET OVERVIEW
The MRTD technique was first presented by Krumpholz and Katehi in 1996 [2]. At the time, multiresolution analysis, the application of wavelet bases to numerical problems, was becoming popular in a number of fields as a way of increasing both the efficiency and accuracy of numerical methods. These techniques remain popular today. Their original MRTD paper provides a general discretization technique that can be used with any wavelet basis, however, the paper focused on the Battle-Lemarie wavelet
6
scheme. Several other wavelet schemes have since been applied in the same manner, including Haar [3], Daubechies [4], and Cohen-Daubechies-Feauveau (CDF) [5]. While the above mentioned wavelet schemes are all very different, they share a number of characteristics that suit them towards numerical modeling. All wavelet
schemes used in for MRTD are orthonormal, and are characterized by a scaling function and a mother wavelet. The mother wavelet is generated using the scaling function and is in turn used to generate all other wavelets that constitute the basis. The orthogonality that exists between the wavelets and scaling function, as well as all wavelets with other wavelets, make them natural choices for numerical discretization. The key concept regarding wavelet expansions is that of levels of wavelet resolution. These levels correspond to sets of functions that can be added to the In the following
expansion to increase the accuracy of the wavelet discretization.
expressions, the scaling functions will be represented as i ( x ) [6], where
i ( x ) =
x i , x
(1)
and (x ) represents the general scaling function. Likewise, the wavelets are offset throughout the grid, but are characterized with three indices. coefficient is represented as, A typical wavelet
ir, p ( x) = 2r / 2 0 2r / 2
x i x
p .
(2)
In this function, r represents the wavelet resolution and p can take any integer value between 0 and 2r-1. Using these coefficients, each level of resolution, r, contains 2r wavelets, offset by x 2 r . For all wavelet schemes used in multiresolution analysis, the following hold [6]:
(x ) (x ) =
i j
i, j
(3)
7
(x )
i
r j, p
( x) = 0 i, j , r , p
(4)
r i, p
( x) s,q ( x) = i , j r ,s p ,q , j
(5)
where i,j is the Kronecker delta function,
i, j =
1 i = j . 0 i j
(6)
If the closed subspace represented by all wavelets of resolution r and higher is termed Vr, the following properties are required, ...V3 V2 V1 V0 V1 V2 V3 ... (7) (8) (9)
UV
iZ
i
= L2 (R )
IV
iZ
i
= 0,
where Z is the set of all integers and L2(R) is the set of all doubly differentiable functions. The consequence of these properties is that each addition of a level of wavelet resolution increases the accuracy of the representation, and arbitrary accuracy can be achieved by using the appropriate wavelet resolution. This is in contrast to methods of improving accuracy in other basis functions, where increased accuracy is achieved by contracting the domain and increasing the number of basis functions. There are several basis functions that can be used to meet these requirements. The oldest and simplest of these is the Haar basis functions. When these wavelets are used in MRTD, a scheme equivalent to FDTD can result [7]. However, when used correctly, a time and space adaptive scheme that is more numerically efficient than FDTD results. As mentioned above, a number of other wavelet schemes have been used with the MRTD technique.
8
An overview of the Haar wavelets is presented here. They are simple enough that they can be used to explain wavelet concepts relatively easily in one and two dimensions, while adequately demonstrating wavelet properties (three dimensional wavelets will be discussed, but are difficult to present in a two dimensional format). The code developed for this work uses Haar wavelets exclusively; they are preferred for a number of reasons that will be discussed. In addition, the composite cell technique presented in the next chapter is practical only for Haar wavelets. The Haar scaling function, , and mother wavelet, , are presented in Figure 1. The Haar scaling function is defined as [0,1), the characteristic function from zero to one. The mother wavelet is based on this function. It is defined as
0 ( x ) = (2 x ) (2( x 1 2 )) ,
(10)
Figure 1:
Haar scaling function (top) and mother wavelet (bottom)
9
Figure 2:
Haar wavelets for resolution 1 and 2
Higher level wavelets can be generated using (2). These wavelets for resolutions 1 and 2 are presented in Figure 2. As stated previously, and demonstrated in Figure 2, there are 2r wavelets for each level of resolution. Thus, there are 2 wavelets at r=1 and 4 at r=2. These functions represent scaled and translated mother wavelets. The scaling factor, 2r/2, presented in (2) is required to orthonormalize the functions, which guarantees (5). The domain of all of the wavelet functions for any level of resolution is identical to that of the mother wavelet and scaling function. A scheme that uses wavelets up to and including a resolution level
r will use 2r+1 functions (including the scaling function).
There are several advantages to using Haar wavelets to model electromagnetic phenomena. The first is their finite domain nature. When Haar wavelets are used to expand an electric or magnetic field, the scaling functions and wavelets from one cell to the next do not overlap. This allows hard boundary conditions (setting discrete areas to a fixed value) to be easily applied. For example, a perfect-electrical-conductor (PEC)
10
boundary condition can be applied by zeroing tangential electric fields. If the basis functions were to overlap, this condition could not be applied. The second advantage of using Haar wavelets is the relative ease of performing derivatives and integrals. Due to their pulse nature, these operations become simple arithmetic. It will be seen that this eases the generation of the MRTD scheme. The major disadvantage of Haar wavelets is in the numerical dispersion of the resultant MRTD scheme. Compared to most other wavelet schemes, denser grids, and thus more coefficients, are required for Haar-MRTD schemes. The MRTD scheme is generated by first representing the electric and magnetic fields in terms of wavelets and then applying the method of moments to Maxwells curl equations, (11) and (12), and the constitutive equations, (15) and (16). The differential form of Maxwells equations in the time domain can be expressed as,
E(t ) = B(t ) M (t ) t
(11)
H(t ) =
D(t ) + J (t ) t
(12) (13) (14)
D(t ) = e (t ) B(t ) = m (t )
where E is the electric field (V/m), H is the magnetic field (A/m), D is the electric flux density (C/m2), B is the magnetic flux density (Wb/m2), is charge density (C/m3), M is the magnetic current density (V/m2), and J is the electric current density (A/m2). These equations relate the electric and magnetic fields at points in space. This analysis is limited to isotropic, non-dispersive media, which have the constitutive equations,
D(r ) = (r )E(r )
(15) (16)
B(r ) = (r )H (r )
11
where r is the position vector. For a general material interface, the boundary conditions are n (D 2 D1 ) = s
n (B 2 B1 ) = 0
n (E 2 E1 ) = 0
(17) (18) (19) (20)
n (H 2 H1 ) = J s
where s, and Js are the charge and electric current densities on the interface and n is the
normal unit vector to the interface. Faradays and Amperes laws, (11) and (12) respectively, are vector equations and are expressed in Cartesian coordinates as,
Bx E y E z Mx = y t z
(21)
B y t
=
E z E x My x z
(22)
Bz E x E y = Mz. y t x
Dx H z H y = Jx t y z
(23)
(24)
D y t
=
H x H z Jy z x
(25)
Dz H y H x = Jz t x y
(26)
To generate the MRTD scheme, each E and H field in Faradays and Amperes laws must be expanded in terms of scaling and wavelet functions. A one dimensional field, for
12
example Ex in a 1D scheme, is expressed in scaling and wavelet functions up to resolution rmax as,
rmax 2 1 E x ( x) = hn ( t ) n Eix , i ( x) + n Eix,r,p ir, p ( x) , , n ,i r =0 p =0
r
(27)
where, n Eix , and n Eix,r,p are coefficients representing the magnitudes of the scaling and , wavelet functions. This function is the discretization in both time and space. The time discretization is performed with simple pulse functions, equivalent to the Haar scaling functions presented in Figure 1. The pulse functions are used in time to ensure causality. If the functions representing each time step overlapped, it would be possible for past events to be modified by future ones. The spatial index, i, indicates the position in cells along the x axis. The cells are defined by the domain of each scaling function, x. It is, technically, possible to use one scaling function the size of the entire domain. Practically, however, this is difficult to implement, and it removes some of the ability to locally increase resolution. In practice, several cells per the maximum excited
wavelength are used. The criteria for choosing the cell size are discussed later, as a part of numerical dispersion. In order to expand the fields in multiple dimensions, the products of all scaling/wavelet functions of each direction must be used. In two-dimensions, this means four groups of coefficients, scaling-x/scaling-y, wavelet-x/scaling-y, scaling-x/wavelet-y, and wavelet-x/wavelet-y coefficients. As in (27), the wavelet terms require sums, and the wavelet/wavelet terms have nested x and y sums. In a two-dimensional simulation, with
rmax=0, the four functions used to represent the fields for the Haar wavelets are depicted
in Figure 3. The maximum resolution is not required to be the same in each direction. For any wavelet basis, the number of coefficients required for a given maximum resolution is
13
D+
# of coefficients = 2
rmax,i i= x , y ,z
,
(28)
where D is the dimensionality of the simulator, and the sum indicates that the maximum resolution may vary by direction. In a three dimensional simulator, there are nine groups of coefficients needed to represent the field expansion, and an increase in resolution in every direction increases the number of coefficients needed by three powers of two.
Figure 3:
Two-dimensional Haar coefficients for rmax=0
The expansion of the Ex field in three dimensions, for any wavelet basis, is
14
n Eix, ,j i ( x ) j ( y ) k ( z ) + ,k r r max 2 1 x , r n Ei , j ,k ,r , p i , p ( x ) j ( y ) k ( z ) + r =0 p = 0 rmax 2r 1 x , n Ei , j ,k ,r , p i ( x ) rj , p ( y ) k ( z ) + r =0 p = 0 rmax 2r 1 E x , ( x ) ( y ) r ( z ) + k k,p r =0 p = 0 n i , j , k , r , p i r 2r 1 r 2s 1 max E x ( x ) = hn ( t ) max x , r s n i , j , k n Ei , j , k , r , p , s , q i , p ( x ) j , q ( y ) k ( z ) + r =0 p = 0 s = 0 q = 0 rmax 2r 1 rmax 2s 1 x , r s n Ei , j ,k ,r , p ,s ,q i , p ( x) j ( y ) k ,q ( z ) + r =0 p = 0 s = 0 q = 0 rmax 2r 1 rmax 2s 1 n Eix, ,j, p ,s ,q i ( x ) rj , p ( y ) ks,q ( z ) + ,k ,r r =0 p = 0 s = 0 q = 0 r 2r 1 r 2s 1 r 2v 1 max max max x , r s v n Ei , j ,k ,r , p ,s ,q ,v ,w i , p ( x) j ,q ( y ) k ,w ( z ) r =0 p = 0 s = 0 q = 0 v = 0 w=0
(29)
The eight groups of coefficients represent all of the combinations of scaling/wavelet functions in each direction. Using (28), for rmax=0 in all directions, there are eight coefficients per cell. For rmax= 2, there are 512 coefficients per cell. The field
coefficients are identified in this equation by the superscripts, denoting the field direction and scaling/wavelet component for each direction. For example, E y , denotes a ydirected electric field coefficient, with a scaling function used in the x and y directions, and a wavelet in the z direction. The coefficients have a maximum of nine indices (for the E dir , coefficient), which represent the position in space of the scaling function as well as the resolution and position of each wavelet. The 3-D MRTD scheme is generated by first representing all fields (E, H ,D, and
B) in scaling/wavelet function expansions, as in (29), and setting the time and space
steps. However, the fields are not collocated in space or time. The original MRTD scheme presented in [2] offsets electric fields one half a cell (the offset between scaling functions) along their coordinate axis (Ex by x/2 in x, for example), while magnetic
15
fields are offset by half a cell in their two normal directions (Hx by y/2 in y and z/2 in
z).
However, two papers [8,9] correct this field arrangement, suggesting that the
magnitude of the offset should depend on the maximum wavelet resolution used in the simulation. The corrected offset, sd, is,
sd = 1 2
rd , max + 2
,
(30)
where d denotes the direction (x, y, or z). The purpose of this offset is discussed in more detail as a consequence of numerical dispersion, but a brief graphical discussion is now presented which explains some of the purpose of the choice of offset. Figure 4 shows a two-dimensional Haar cell with rmax=0. The r=0 wavelets are shown as a reference. The MRTD technique is characterized by cells that contain several basis functions. As a consequence, the field variation throughout the cell is not limited to the shape of a single function, but rather the sum of several. The position of the wavelet coefficients at rmax, which can be calculated using (2), defines a number of equivalent
grid points [9]. Specifically, two equivalent grid points are defined per rmax level wavelet
positioned x (4 2 rmax ) and 3 x (4 2 rmax ) from the position specified in (2). It should be noted that the number of equivalent grid points is equal to the number of basis functions, (28). While the half cell offset used in [2] was chosen to give the cells an offset equal to Yee-FDTD [1], the offset presented in (30) arranges the equivalent grid points in the Yee manner. It should be noted that while Figure 4 is shows the wavelets for the Haar scheme, the position of the equivalent grid points is valid for any wavelet basis. Furthermore, the cell size represents the domain of the scaling function in the Haar case, but it represents the spacing between the centers of the scaling functions in the general case. The equivalent grid points in Figure 4 are marked with Xs.
16
Figure 4:
2-D Haar cell, rmax=0, equivalent grid points
If a 2-D scheme with EX, Ey, and Hz (TEz mode) is simulated with the rmax=0 cell of Figure 4, and the fields are offset using (30) in the manner suggested above, the field arrangement shown in Figure 5 results. This figure shows the equivalent grid points for all three fields, and is valid for any wavelet basis. The arrangement of the grid points in this case is the same as the Yee-FDTD cell. If the same arrangement is used for a full 3D scheme, the equivalent grid points are arranged in the same manner as Yee-FDTD.
17
Figure 5:
Equivalent grid points for a 2-D MRTD cell, rmax=0
Like FDTD, the update equations in MRTD are fully explicit. Thus, the time step must be chosen according to a stability requirement. This is discussed later. Also like FDTD, the electric and magnetic fields are offset in time. The time basis functions used for all MRTD schemes are pulses, as in FDTD. Similarly, the fields are offset by half a time step. Once the field offset is determined, the MRTD update scheme can be
determined using the method of moments.
2.1.2 METHOD OF MOMENTS OVERVIEW
The method of moments is a general method used in a variety of numerical schemes. A brief overview of techniques that relate to MRTD is presented here. First, the method is applied by representing the unknown function as a sum of unknown coefficients multiplied by known basis functions. Next, a series of testing functions are
18
chosen. These functions are multiplied into both sides of the equation and the integral is taken over the entire domain. This process is referred to as an inner product, and the following notation will be used,
f , g = f (r )g (r )r
(31)
If the number of testing functions is equal to the number of unknown coefficients, the result is a series of linear equations that can be solved to determine the coefficients. As a simple 1-D (in space) example, the equation
f (x ) =
b( x ) x
(32)
is presented, where f(x) is unknown and b(x) is a known function. expression can be expanded into a sum of basis functions,
The unknown
f ( x, t ) = a n cn ( x )
n =0
N
(33)
where an are unknown scalar coefficients and cn(x) are the known basis functions. A set of testing function, wn(x) are chosen, and the inner product of (32) with each
wn , f = wn , b
(34)
is taken. The result is a system of equations,
a1 w1 , c1 a1 w2 , c1 M a1 wN , c1
a2 w1 , c2 a2 w2 , c2
a N w1 , c N w1 , b = w2 , b , M O a N wN , c N wN , b
L
(35)
which can be solved to yield the unknown coefficients, an. If the testing functions are chosen such that,
19
wm , cn = m,n
(36)
then the matrix on the left hand side of (35) is diagonal. In this case, the scheme is explicit; a matrix does not need to be inverted to solve the equations. Each coefficient can be determined directly,
a n = wn , b .
(37)
If an orthonormal set of basis functions is chosen,
wm , wn = m,n ,
(38) This technique is called
then the basis functions can be used as testing functions. Galerkins method.
2.2
MRTD Update Scheme
2.2.1 TIME LOCALIZATION
MRTD update equations are determined using Galerkins method with wavelet discretizations of the electric and magnetic fields, such as (29). The half time step offset of the electric and magnetic field ensure that the fields are always updated from values at previous time steps, and thus the only unknowns are the updated values of the fields. Using (21)-(23), B fields are determined from previous E field values. Likewise, (24)(26) are used to determine D fields from previous H fields. In many FDTD schemes the constitutive relationship, (15) and (16), is applied discretely at each cell, and thus H values are determined directly from E values, and vice-versa. For general MRTD basis functions, which can extend over the entire domain, this is not possible. Furthermore,
20
even for non-overlapping, finite-domain wavelets, such as Haar wavelets, it is not practical because the dielectric constant cannot vary over the relatively large MRTD cell. Thus, the method of moments must also be applied to (15) and (16), to determine E fields from D fields and H fields from B fields. While this does add complexity to the method, it also allows the material to vary continuously throughout the cell. The MRTD update equations are determined by localizing the coefficients in time by testing with the time basis functions (pulses) and then localizing in space by testing with the scaling/wavelet functions. The time derivative of the pulse function yields two Dirac Delta functions, located at the edges of the pulse, hn (t ) = (t nt ) (t (n + 1)t ) , t
(39)
The time derivatives of the pulse functions that make up the time discretization are a Delta train, represented by the top line of Figure 6. In (21)-(26), one field value is differentiated in time, while the other field values are differentiated in space. As stated previously, the E and H fields are offset in time by half a time step, which means that the spatial derivative terms are offset in time by half a time step from the time derivative terms.
Figure 6:
Time localization using pulse functions
21
The time testing function is chosen to overlap with the terms in the spatial derivative. If the B (and H) coefficients are located at the n+1/2 time steps and the D (and E) functions are located at the n time steps, Figure 6 represents the time discretization of (21) (23). The delta train on the top of the figure represents the time derivative of the B fields, while the pulse functions on the bottom represent the E fields. The shaded area represents the domain of a testing function, which is collocated in time with the E fields. If the time and space basis functions that represent the E fields are separated,
E (r, t ) = hn (t )n E (r ) ,
n =0
N
(40)
where nE(r) is the wavelet/scaling space discretization for time step n, the quantities in the brackets of (29). Similarly,
B(r, t ) = hn+ 12 (t )n+ 12 B(r ) ,
n =0
N
(41)
represents the B fields. Using these expressions, the inner products with hn are,
hn (t ), E (r, t ) = t n E (r ) ,
(42) (43)
hn (t ),
B(r, t ) = n + 12 B(r ) n 12 B(r ) , t
As an example, (42) and (43) can be applied to (21). Ignoring the current terms, and taking the inner product with hn,
hn , E Bx E = hn , y z t z y
(44)
yields,
22
n + 12
E (r ) E (r ) Bx (r ) n 12 Bx (r ) = t n y n z y z
n+ 12
(45)
Bx (r ) ,
Similar expressions can be found using (22) to (26). Solving for
n + 12
E (r ) E (r ) Bx (r )= n 12 Bx (r ) + t n y n z y z
n + 12
(46)
yields an update equation for
Bx (r ) in terms of quantities from the previous time step.
The next step in determining the MRTD update equations is to test with the spatial scaling/wavelet functions.
2.2.2 SPACE LOCALIZATION
The 3-D basis functions are expressed as ( x )( y )( z ) , where each can be either
l or r p , with l representing a directional index, i, j, or k. As stated in (28), for
l,
maximum resolution in the x, y, z direction rmax,x, rmax,y, rmax,z,, respectively, each cell (denoted by the triplet i,j,k) contains 2
3+ rmax, x + rmax, y + rmax, z
basis functions. To determine an
update formula for each coefficient, the testing functions are collocated in space with the time differentiated functions in Faradays and Amperes laws. In (46) this is equivalent to collocating the scaling/wavelet testing functions with the Bx fields. A convenient notation was first introduced in [8], which represents the field discretization as a vector. This notation has been generalized for any wavelet basis, and allows the MRTD update equations to be written in a compact and easily understandable form. The wavelet/scaling coefficients for all wavelet resolutions can be represented as a vector, for example,
23
n+ 12
B x ,i , j , k
x , n + 1 2 Bi , j ,k x , n + 1 2 Bi , j ,k , 0 , 0 M x , rmax, z n + 1 2 Bi , j ,k ,r 1 max, z , 2 x , n + 1 2 Bi , j ,k , 0 , 0 , 0 , 0 = M x , r r n + 1 2 Bi , j ,k ,r , 2 max, y 1, rmax, z , 2 max, z 1 max, y x , n + 1 2 Bi , j ,k , 0 , 0 , 0 , 0 , 0 , 0 M x , n + 1 2 B i , j , k , rmax, x , 2 rmax, x 1, rmax, y , 2 rmax, y 1, rmax, z , 2 rmax, z 1
(47)
If another vector, , is defined,
i , j , k
i ( x) j ( y ) k ( z ) 0 i , 0 ( x) j ( y ) k ( z ) M rmax, z i ( x) k ( y ) k , 2 rmax,z 1 ( z ) 0 0 i , 0 ( x) j ,0 ( y ) k ( z ) = M rmax, z ( x) rmax y ( y ) k , 2 rmax,z 1 ( z ) r i j , 2 max, y 1 i0,0 ( x) 0, 0 ( y ) k0, 0 ( z ) j M r z rmax x ( x) rmaxr y ( y ) kmax,max,z 1 ( z ) r ,2 j , 2 max, y 1 i , 2 rmax,x 1
(48)
then,
n+ 1 2
Bx (r ) = T, j ,k n+ 12 B x ,i , j ,k i
i , j ,k
(49)
The update equation for each coefficient of B can then be determined by taking the inner product of (46) with each scaling/wavelet coefficient (every row of (48)). Because the basis functions are orthonormal, a separate update equation is found for each scaling/wavelet coefficient. However, due to the offset of the B and E fields, the
24
scaling/wavelet functions used for the B fields are not orthogonal with the functions used for the E fields. If the basis functions representing the Bx field are offset in the positive y and z directions from cell i,j,k (located at x, jy, kz), while the Ey fields are offset in the y direction and Ez is offset in the z direction, then, using the vector notation of (47), the update for (46) becomes,
n + 12
B x , i , j , k = n 12 B x , i , j , k +
t Bx Bx U E y , m n E y , i , j , k + m + U E z , m n E z , i , j + m , k . xyz m m
(50)
In this expression, U represents an update matrix, and consists of the inner products of the E basis functions with the B basis functions. The B and E fields are only offset in the direction of differentiation of the E field. Thus, m represents the position of all of the neighboring cells that the E and B fields overlap (depending on the domain of the chosen basis function). The U matrices are 2
3+ rmax, x + rmax, y + rmax, z
x 2
3+ rmax, x + rmax, y + rmax, z
in size, and can be computed
before simulation begins. They take the form,
= 1 F n 2 m F1 2 , F L n 2 m M F1 L , F L n 2 m
F1
1 ,
U F12 ,m F
F 2 n 2 m F1 2 , F 2 n 2 m
F1
1 ,
L
O
F2 L m n . L , F2 L F1 m n
F1
1 ,
(51)
where F1 is the field whose update is being found, and F2 is the field on which the update depends. In this case,
F1
l is the lth row of for the F1 fields, and
F2
l
m
is the lth row of term represents
for the F2 fields, offset by m in the direction of differentiation. The
n
the space derivative, in this case n=y or z. If the basis functions used for each field were collocated, (51) would be diagonal. However, due to the offset scheme used, the
functions are only collocated in two of the three Cartesian directions. For example, the
25
Bx and Ey fields have the same location in the x and y directions, but are offset in the z direction. If the field offset in the x, y, and z directions is sx, sy, sz, respectively, then the (2,2) entry of (51) for E y ,i , j , k (same cell as the B vector, m=0) is,
U Bxy1,m, 2, 2 = i0, 0 ( x) j + s y ( y ) k + sz ( z ), E
0 i , 0 ( x) j +s y ( y ) k ( z ) = z 0 0 i,0 ( x) j+sy ( y)k +sz ( z ) z i,0 ( x) j+sy ( y)k ( z ) xyz
(
)
(52)
This integral can be separated by direction,
0 i ,0
( x) j + s y ( y ) k + sz ( z )
0 i ,0 ( x) j + s y ( y ) k ( z ) xyz = z ( z ) 0 0 i,0 ( x) i,0 ( x)x j+sy ( y) j+sy ( y)y k +sz ( z ) kz z
(
)
(53)
Next, the orthogonality between the collocated basis functions can be exploited,
0 i ,0
( x) i0,0 ( x)x j + s y ( y ) j + s y ( y )y k + sz ( z ) k ( z )z = x y k + sz ( z ) k ( z ) z z
(54)
and, finally, only a single one dimensional integral must be evaluated. Depending on the functions involved, this can be done analytically or numerically. For either case, it can be tabulated, and does not have to be performed for each simulation. It should also be noted that, due to the orthogonality of the basis functions, the majority of the entries in (51) are zero. The procedure shown above can be applied to the remainder of (22) (26), yielding update equations for all B and D fields. The remaining step in determining the MRTD update scheme is to find the updates for E and H fields from the D and B fields, respectively. As these fields are collocated in space and time, the procedure is less complicated, although potentially more computationally intensive.
26
2.2.3 MEDIA DISCRETIZATION
In a general structure, the material properties can vary as a function of position. For this work, only linear, isotropic, nondispersive media are considered. The
relationship between the fields in this case is given in (15) and (16). Time localization in this case can be performed by taking the inner product with the time yielding, E d ( x, y , z ) = H d ( x, y , z ) = 1 D ( x, y , z ) ( x, y , z ) d 1 B ( x, y , z ) ( x, y , z ) d d = x, y, or z (55)
(56)
The update for each component of E and H can be determined by taking the inner product with the scaling functions. However, because the offset of each field component is different, the update for each component is different. component takes the form
n
The update of a general E
E d ,i , j , k = U E dd ,i , j , k , a ,b , c n Dd ,i + a , j + b , k + c D
a b c
(57)
where a, b, and d represent the relative positions of the surrounding D basis functions that overlap the (i,j,k) E basis function. The number of entries in these sums is dependent on the basis functions. The coefficients of (57) take the form
E U Ddd , a ,b , c ij =
1 xyz ( x, y , z ) E d i Dd j a , b , c
(58)
Unlike the U matricies that are used in the discretization of Faradays and Amperes laws, most of the entries in these update matricies are nonzero. Furthermore, the matrices vary by position. However, for regions of uniform media, they are diagonal for
l=m=n=0, and zero elsewhere.
27
These equations can be used together in a time stepping scheme with the discretizations of Amperes and Faradays laws to create a general time marching MRTD scheme: 1. 2. 3. 4. 5. Determine B fields from E fields Determine H fields from B fields Determine D fields from H fields Determine E fields from D fields Repeat until simulation is complete
2.2.4 NUMERICAL STABILITY
The scheme that has been presented allows MRTD update schemes to be determined for any basis function. These schemes are fully explicit, the matrix
expressions demonstrated only involve multiplication, and all numerical integration is performed before the start of simulation. For the majority of cases, a library of the update, equations, U, can be constructed, so that the costly numerical integration does not have to be performed before each simulation. Most of the update matrices are very sparse, and a number of clever coding techniques and libraries can be used to reduce the processor and memory load of implementing these updates. Any basis function can be used by calculating the U matricies, which can be accomplished by solving a number of integrals, either analytically, or, more likely, numerically. To this point, however, no guidance has been given on choosing the time and space step values. To properly choose these values, the stability and dispersion of the scheme must be taken into consideration. The purpose of this investigation is to determine techniques that can be used with the Haar-MRTD scheme to reduce simulation time by increasing the efficiency and applicability of the technique. The techniques that were developed are not directly
28
related to the stability or dispersion requirements, and as such they will not be covered in depth. However, knowledge of the effects of these conditions is necessary to apply the techniques, and as such a brief overview is presented here. It is a known limitation of explicit techniques that the time step must limited as a function of the space step to maintain stability. If the time step is not below this limit, the resulting scheme quickly grows without bounds. One notable exception to this is the ADI-MRTD method [10], which is a semi-implicit scheme that remains stable (but not accurate) for any time step. This technique separates the time-step into two substeps, and, as it does not alter the spatial grid, it is very similar to ADI-FDTD [11]. The techniques presented in this study are related to the spatial discretization, and are completely compatible with the ADI technique. The discussion, however, will focus on the traditional MRTD technique. A brief stability analysis of the MRTD technique is presented in [2], while additional studies were presented in [9,12]. The stability condition for MRTD depends on the choice of the wavelet basis. The general dispersion analysis, however, is performed in the same manner as traditional FDTD techniques [13]. The time update portion (the discretization of the time
derivative) and the space update portion (the discretization of the spatial derivative) are split into two separate problems. For numerical stability, the eigenvalues of the spatial problem must contain the eigenvalues of the temporal problem. For FDTD, this results in the condition [13]:
t
1 1 1 1 1 2 + 2 + 2 x y z
.
(59)
For Battle-Lemarie S-MRTD (scaling function only) [12]
29
t
1
n a 1 i =0
a
i
1 1 1 1 2 + 2 + 2 x y z
(60)
where
sMRTD =
n a 1 i =0
a
i
,
(61)
the MRTD stability factor, is the sum of the magnitudes of the inner products of the offset scaling functions from the spatial update equations (equivalent to the entries in the 1x1 U matricies using the notation presented here). In this equation, na is the stencil size, or the domain of the basis function in cells. The Battle-Lemarie basis functions are symmetric. Alternatively, the sum could be expressed as
sMRTD =
i = nb
a
na
i
,
(62)
where na and nb represent the stencil size of the function in both the positive and negative directions, respectively. In fact, if this expression is used, the stability condition is valid for any basis function. For the Battle-Lemarie scheme [12],
sMRTD = 1
n a 1 i =0
a
= 0.6371 ,
(63)
i
and the time step for the MRTD scheme is smaller than the time step for and FDTD grid with the same cell size. However, as will be shown next, the cells used in Battle-Lemarie MRTD are generally larger than those used in FDTD, and thus the technique is more efficient overall. It is demonstrated in [9] that, for any MRTD wavelet basis, the
resolution of the scheme doubles for each addition of a level of wavelet resolution. This
30
is equivalent to saying the time step of a scheme with one level of wavelet resolution in each direction will have a time step one-half of the value of S-MRTD. The resulting stability condition, valid for any MRTD basis, is,
t 1 (64)
i = nb
a
na
i
1 1 1 1 + + 2 2 2 x y z r +1 rmax, y +1 rmax,z +1 2 max,x 2 2
The Haar MRTD technique can be analyzed using the same method. Several authors [3,7,8,9], have noted the equivalence between Haar S-MRTD and FDTD, and indeed this relationship can be exploited to directly derive the general Haar MRTD stability requirement. Using this equivalence, and applying the resolution doubling
argument to (59), the stability for Haar MRTD at any resolution level is [9],
t 1 1 1 1 1 + + 2 2 2 x y z rmax,z +1 rmax, y +1 rmax,x +1 2 2 2
.
(65)
This is a special case of (64), for the Haar scheme
sMRTD = 1 .
(66)
31
2.2.5 NUMERICAL DISPERSION
Dispersion is defined as the variation of the wavelength, , with frequency, f. More commonly, it is discussed as the variation of the wavenumber, k, with angular frequency, , where
= 2f
k= 2
(67) . (68)
By substituting the solution for a monochromatic plane wave into Maxwells equations (assuming a linear, isotropic, nondispersive medium), the following dispersion relationship results k= where,
c=
1
c
,
(69)
.
(70)
In a 3-D representation, a wavevector, k, is defined, where, in Cartesian coordinates, k = k xi + k y + k z k , j and k = kx + k y + kz . Using these parameters, a phase velocity, vp and group velocity, vg can be defined, vp = (72) (71)
k
= c
(73)
32
vg =
= c . k
(74)
These parameters demonstrate that the wavelength and frequency have linear relationship, and that the phase and group velocities are independent of frequency. These relationships are more complicated in a numerical scheme. A numerical scheme, by its nature, discretizes space into small, but finite, cells. In these cells, the waves cannot propagate in any direction, but rather in the directions defined by the grid. In addition, time is not continuous, but set as a multiple of a discrete time step. With these limitations it should not be surprising that the numerical wave propagation velocity can depend on direction, and, because of the discrete spatial and time steps, frequency. As it has been shown, Haar S-MRTD is equivalent to FDTD, and, as one of the simplest MRTD schemes, it will be used as the first example. For S-MRTD (and FDTD) [14],
k y y 1 1 k x x 1 k z z t 1 ct sin 2 = x sin 2 + y sin 2 + x sin 2 . (75)
2 2 2 2
For numerical stability, the argument was advanced that increasing the resolution by one level effectively doubled the resolution of the code. This was represented in the time stability condition by dividing the space step by 2rmax +1 . For the dispersion analysis, this condition holds as well, thus [9],
+1 r 1 t 2 max,x k x x sin sin rmax,x + 2 = ct 2 x 2
2 2
2rmax, y +1 k y y 2rmax,z +1 k z z sin rmax, y + 2 + sin rmax,z + 2 + 2 y 2 x
2
2
.
(76)
A similar dispersion analysis was performed in [12] for Battle-Lemarie S-MRTD, and, when made general for any wavelet basis and resolution level,
33
2 r +1 n 1 1 2 t 2 max,x a ai sin k x (i + 1 + 2)x sin = r ct 2 x i = nb 2 max,x
2
2rmax, y +1 na 1 k (i + 1 2 )y ai sin y r + 2 . + 2 max, y y i = nb
2
(77)
2 rmax,z +1 na 1 ai sin k x ( z r+ 1 +2 )z + 2 2 max,z z i = nb
2
This is the general dispersion relationship for MRTD. As (65) is a special case of (64), (76) is a special case of (77). It is useful to express the wavenumber in terms of angular frequency. As (77) is a sum of sinusoids, an analytical solution cannot be found for the general case It can, however, be solved numerically for specific wavelet bases. Several studies have been performed using these equations [4,9,12]. For reference, the solution for the wavenumber as a function of angular frequency for Haar MRTD in one of the grid major directions (along an axis), assuming a uniform grid size , is [14], 2 k = sin and the phase velocity is, vp =
1
1 S sin S N
,
(78)
1 S N sin sin S N
1
c,
(79)
using the definitions: N =
o
=
c
(80)
S=
ct .
(81)
34
Using this expression, the error between the numerical phase velocity and c can be determined, as well as what values of N (the number of cells per wavelength) yield stable (non-damped) results. This equation is also valid for any wavelet resolution if is the spacing between equivalent grid points. To keep phase error low, it is usual for these schemes to use more than 10 equivalent cells per smallest excited wavelength. This condition is one of the major arguments for using the MRTD technique. Although it was demonstrated for Haar MRTD that the dispersion relationship, and thus the number of cells needed per wavelength, is the same for the MRTD scheme as in FDTD, other basis functions can use significantly fewer cells per wavelength. It was reported in [12] that a discretization of three to four cells per wavelength using BattleLemarie MRTD gave results comparable to FDTD with 10 cells per wavelength.
2.3
IMAGE THEORY FOR PEC MODELING
The techniques that have been presented thus far allow the modeling of arbitrary variation of the permeability and permittivity of the media under simulation. However, most structures of interest to microwave designers also contain metals. For many
simulations in the time domain, these metals are treated as PECs, and this limitation will be used in this case as well. The treatment of a metal as a PEC is necessary, because otherwise the metal itself must be simulated (resulting in significantly increased number of grid points) or frequency dependent loss characteristics must be applied to the PEC. In the prototypical time-domain technique, FDTD, PEC structures are simulated in a simple and straightforward manner. The PECs are located along grid boundaries, and the electric field values that are tangential to the PEC structures are zeroed each time step (after their updates are calculated). This explicitly enforces (17). The normal
coefficients do not require any special processing, because they are offset from the PEC
35
location. Due to the extended, and often overlapping, nature of MRTD basis functions, this simple processing is not possible in MRTD. Instead, PEC boundary conditions are traditionally applied to MRTD codes through image theory. Several papers have been published that discuss the implementation of image theory to the modeling of PEC structures in MRTD [2,9,12,15,16,17]. Using this
method, PECs are modeled through the introduction of artificial image scaling and wavelet coefficients on the opposite side of the PEC. These image wavelets have
opposite magnitude of the wavelets whose domain include the PEC. In this manner, when the real and image wavelet and scaling coefficients are summed at the PEC interface, the result is zero, the PEC boundary condition. This technique is difficult to apply in a general 3D code, because a large amount of bookkeeping is required for each PEC in the grid. For any wavelet that contains a PEC within its stencil, a unique update scheme must be applied. In addition, multiple images, caused by the close proximity of multiple PECs, can exacerbate this problem further. While it has been stated in several publications that this technique can be applied to complex structures, no publications have been made detailing its use for a general 3D structure with multiple PECs. At maximum, external boundaries and one internal PEC have been presented. A one dimensional example of the application of image theory is presented in Figure 7. This example uses biorthogonal scaling functions [5]. These functions are used because they have a relatively compact stencil size of four cells. The shaded area represents the width of the stencil. All magnetic field points within this distance of the PEC must use the images of the electric field on the opposite of the PEC when computing an update. If this PEC intersects the grid, (the simulation is performed on both sides of the PEC) the actual field values on the opposite of the PEC are not used (the PEC effectively decouples these fields). In addition, when the fields are reconstructed, the images must be used in the reconstruction algorithm. This same method can be used to
36
simulate perfect magnetic conductors (PMC) by imaging the magnetic fields across the boundary.
Figure 7:
Image theory in one dimension using biorthogonal scaling functions
2.4
EXISTING PML TECHNIQUES IN MRTD
The MRTD scheme that has been presented allows the simulation of arbitrary media terminated with PEC boundary conditions. PEC boundary conditions are
necessary because the fields on the boundary of the simulation space must be known in order to calculate the internal field values. However, the PEC boundary conditions cause complete reflection of electromagnetic waves. When used in a time domain simulation, these reflections interfere with the fields being measured, and render the results useless. In order to reduce these reflections, absorbing boundary conditions (ABCs) must be applied.
37
The most popular absorbing boundary condition used in time-domain electromagnetic simulations is the perfectly matched layer (PML). This boundary
condition was first introduced by Berenger in 1994 [18]. The PML takes the form of a material that surrounds the simulation space that is perfectly matched to all incoming wave directions and frequencies, and is also highly lossy. Berengers original
formulation of the PML used a split field formulation of Maxwells equations. However, a similar method was later introduced that was termed the uniaxial PML (UPML) [19]. This PML takes the form of a unaxial anisotropic medium; the material properties are identical in directions transverse to the UPML interface, and vary in the normal direction. This PML formulation is often implemented because it avoids the field splitting used in Berengers original PML formulation. In addition, when implemented in FDTD, it
requires two field update steps. The first step applies Faradays and Amperes law, and the second step applies the constitutive relationships. This is a natural method to
implement in MRTD because the fields are already updated in this manner. Several authors [17,20,21] have demonstrated the implementation of the PML for S-MRTD (scaling function only) schemes. In addition, a technique that interfaces HaarMRTD and FDTD domains was introduced [22], with one application aimed at interfacing an FDTD PML with an inner MRTD grid. A derivation of the UPML for the general Haar MRTD scheme (for any wavelet resolution) has not been presented. The derivation of the technique for any wavelet resolution is given in Chapter 4.
2.5
FDTD
It has already been stated that FDTD is often used as a benchmark for other time domain schemes, such as MRTD. The chief advantage of the FDTD method is its simplicity; it was originally derived using central differences with Maxwells equations.
38
However, the key aspect of the FDTD method does not come from its relatively simple derivation, but the arrangement of it grid points. FDTD was originally developed in 1967 by K.S. Yee[1]. The feature of FDTD that makes it both simple and widely applicable is the Yee-Cell, which is presented in Figure 8. A simple Cartesian grid is used, with the fields offset from the grid points in a specific manner. In Figure 8 the electric field components are offset from the grid points half a cell in their coordinate direction (Ex in x, for example) and the magnetic fields are offset in the directions normal to their field components.
Figure 8:
Yee-Cell in 3-D
Using this simple field arrangement, and by offsetting the magnetic and electric fields by half a step in time, central differences can be applied to Faradays and Amperes law, creating a system of equations that can be used in a time marching scheme to update the fields. For example, the FDTD discretization of (21) is, H x ,i , j , k = n 1 / 2 H x ,i , j , k + t n E y ,i , j , k n E y ,i , j , k 1 n Ez ,i , j , k n Ez ,i , j 1, k , z y (82)
n +1 / 2
39
if the convention that has been used in the MRTD case is applied, where all field components are offset in the positive direction, and the offsets are not explicitly shown in the field components ( E y ,i , j , k is located at ix, (j+1/2)y, kz, for example). In this equation, the constitutive relationship (15), is directly applied, as the updates are for actual field values at points, instead of scaling/wavelet coefficients that span a much wider area. It is useful to note that the FDTD update equations are exactly the same as Haar SMRTD (scaling function only) update equations. The method-of-moments, applied with pulse basis functions, offset according to the scheme that has been presented for MRTD and FDTD, is equivalent to performing a central difference approximation. Furthermore, FDTD and Haar MRTD are equivalent for any rmax if all wavelets are used for all points in the grid. This can be demonstrated by transforming the MRTD scaling/wavelet
scheme into a pointwise (where the actual field values are used) scheme, [23], or by noting that the Haar MRTD dispersion relationship for an arbitrary resolution level is the same as FDTD, if the FDTD cells are the same size as the MRTD equivalent grid points.
40
CHAPTER 3 MRTD COMPOSITE-CELL MODELING
The MRTD scheme that was presented in the previous chapter provides a framework for applying multiresolution principles to the simulation of modern microwave structures. Every paper published on MRTD discusses two major advantages of the MRTD technique: 1. The basis functions employed in the MRTD technique allow fewer cells to be used per maximum excited wavelength (compared to FDTD), allowing a more efficient overall simulation. 2. Wavelets can be added and subtracted dynamically during simulation, resulting in an adaptive scheme that automatically tailors its computational requirements to the complexity of propagating signals. The number of cells required per wavelength is discussed in the previous chapter. The dispersion relationship, (77), can be numerically solved for any wavelet basis. It is shown in several papers [2,4,9,20] that a variety of wavelet bases can be employed that allow less dense grids to be used. There are two distinct costs of applying these
techniques. The first is determining the update coefficients, such as those used in (50) and (57), that result from Amperes and Faradays laws. However, these coefficients can be calculated before the start of simulation, and it is possible to build a library of the coefficients for several common cases. The second cost of applying the technique is the field updates themselves, which, due to the stencil of the basis functions, involve several surrounding field coefficients. Thus, while fewer cells are required per wavelength when compared to FDTD, more operations are required to calculate each field update. A comprehensive comparison of the benefits from using fewer cells in MRTD compared to the more simple updates used in FDTD has not been presented in literature.
41
Another difficulty of using general wavelet bases is that of simulating lossy media. The derivations of MRTD update equations presented in the previous chapter neglected loss. When ohmic losses are added (considering both electric and magnetic loss),
Faradays and Amperes laws become, E(t ) = H (t ) = B(t ) B(t ) t D(t ) + D(t ) . t (83)
(84)
The procedure for determining the MRTD update procedures described in the previous chapter can be applied in this case as well. For example, the update for the Bx field is determined from, Bx E y E z = Bx . t z y (85)
To determine the time localized form of this equation, the inner product of both sides of (85) must be taken with hn+1/2(t). This process is slightly different than the case without loss, because both the time derivative of the B field and the B field itself are used in this equation. In Figure 9, it is noted that the pulse used in the time localization overlaps the B field from two time steps. The inner product, hn (t ), Bx =
n + 12
Bx (r )+ n 12 Bx (r ) 2
(86)
therefore, includes terms from both the time step being updated and the previous time step. This is a statement of the semi-implicit approximation used in FDTD [14], but it is derived in this case as a direct result of the chosen basis functions.
42
Figure 9:
Time localization (pulse basis functions) for lossy case
Using (86), the time localized form of (85) is, Bx (r )= n 12 Bx (r ) t
n + 12
n + 12
Bx (r )+ n 12 Bx (r ) 2
E (r ) E (r ) + t n y n z . y z
(87)
Collecting terms,
t t E (r ) E (r ) n 12 Bx (r ) + t n y n z , n + 12 Bx (r ) =1 1 + 2 2 y z
(88)
and finally solving for
n + 12
Bx (r )
n+ 1
2 t 2t n E y (r ) n Ez (r ) Bx (r ) = 2 + t n 12 Bx (r ) + 2 + t z y , 2
(89)
When the space localization is performed to determine the updates for the individual wavelet coefficients, the following results,
43
n + 12
B x ,i , j , k = U B xx , a ,b , c n 12 B x ,i + a , j + b , k + c B
a b c
U B xy , a ,b , c n 1 E y ,i + a , j + b , k + c , E 2 t a b c + Bx xyz + U E , a ,b , c n 1 E z ,i + a , j + b , k + c z 2 a b c
(90)
using the notation presented in the last chapter. In this case, however, it is noted that the basis functions can be offset every direction (denoted by a,b, and c, in the x,y, and z directions, respectively), instead of the single direction that was denoted by m in Chapter 2. The entries of the U Bxx ,a ,b ,c matricies take the form B
U ij =
Bx i , Bx j
a ,b ,c
2 t = 2 + t Bx i Bx j a ,b ,c xyz ,
(91)
and the entries of the U Bxy ,a ,b ,c matricies take the form, E U ij =
Bx
i , E j
a ,b ,c
2t = B i E j z 2 + t
a ,b ,c
xyz ,
(92)
where j
a ,b,c
denotes the scaling/wavelet function offset by a, b, and c, (cells) in the x,
y, z directions, respectively, from the coefficient update being calculated. In the general case of ( x, y, z ) varying constantly in each direction, the B x coefficients depend on surrounding B x coefficients from the previous time step. This complicates the update of the coefficients, as when implemented on a computer, the array of updated values must be kept separate from previous values. In the lossless scheme, the fields being updated depend only on their previous value, and the updated fields can be stored in the same array as the previous fields. This doubles the amount of memory required to perform the scheme. In addition, it significantly increases the computational burden. In the lossless case, the updates of any one coefficient depend on the previous value of the coefficient and the values of two other fields (the normal E fields in the B update case, and the normal H fields in the D case). For these fields, a sum of coefficients
44
(multiplied by the entries in the U matrix) from neighboring grid points is required, and the number of grid points depends on the stencil of the basis function. For example, a stencil of twelve (six in the positive, six in the negative direction) cells is normally used for Battle-Lemarie wavelets [16]. However, because the fields are only offset in one direction, the sum must only be performed in one direction (as in (50)). In the case including loss, however, the addition of into the inner products destroys the orthogonality in the inner product integrals, and the sums must be performed in all three directions. For a stencil size n, the result is a summation over n3 elements. As an example, in the lossless case, for Battle-Lemarie wavelets using a stencil size of 12, the total number of neighboring cells required for the update is: Previous Value+Normal Field 1+Normal Field 2
1 + n + n = 2n + 1
(93) (94)
1 + 12 + 12 = 25
while for the lossy case:
n3 + n3 + n3 = 3n3
1728 + 1728 + 1728 = 5184
(95) (96)
This significantly complicates the method. However, for stencil sizes of one (such as Haar basis functions) the requirements are the same in the lossy and lossless case. The final difficulty of simulating complex structures using general wavelet bases is that of representing PEC structures within an MRTD grid. In the previous chapter, image theory was discussed; the primary example presented, as well as the examples shown in the referenced literature, was for the case of a PEC wall terminating the structure. A general method for simulating arbitrarily placed internal PEC structures has not been presented in literature. This makes simulation of most microwave structures difficult. One notable exception is [22], which utilizes an FDTD/Haar-MRTD interface to simulate
45
highly detailed structure in FDTD, and MRTD to simulate open areas. However, this technique cannot be readily expanded to other wavelet bases. In this chapter a method is presented that allows arbitrarily placed PEC structures to be simulated, however, it is shown that it is only practical for Haar-MRTD schemes. As most structures cannot be readily simulated with general wavelet bases, the focus of this thesis is on Haar MRTD techniques. In this chapter a technique is presented that allows subcell structures to be modeled within Haar MRTD cells. This technique is shown to be a bridge between pointwise field updates, such as those used in FDTD, and the wavelet/scaling updates used in MRTD. This technique allows the adaptive
resolution characteristics of MRTD to be exploited for the simulation of any structure. Thus, of the two stated goals of MRTD, it is noted that the first can only be obtained for structures consisting entirely of lossless, PEC-free structures, while this thesis demonstrates how to achieve the second for general structures using Haar basis functions.
3.1
GENERAL, EXPLICIT, SUBCELL FIELD MODIFICATION
The motivation of this research is to represent arbitrary structures in MRTD so that the time-space adaptive grid can be used effectively. To accomplish this, it is necessary to determine a method of representing PEC structures that are smaller than an MRTD cell. The method that has been developed as a part of this research is based on the method that PEC structures are represented in FDTD grids. It is shown that the method effectively provides a bridge between pointwise electromagnetic effects, and the distributed wavelet field representation used in MRTD. In FDTD, which is equivalent to Haar S-MRTD, PEC structures are explicitly represented by zeroing electric field values that are tangential to PEC field locations. The update equations for FDTD can be determined in the same manner as the MRTD
46
equations in the previous chapter. In fact, the Bx update equation is a special case of (50). Similarly, the Dx update equation is [14],
H H Dx ,i , j , k = n Dx ,i , j , k + t n +1 2 z ,i , j , k n +1 2 z ,i , j 1, k y
n +1 2
n +1
H y ,i , j , k n +1 2 H y ,i , j , k 1 , (97) z
in the lossless case (for the purpose of applying PEC boundary conditions, loss is irrelevant). In FDTD, the constitutive relationships are applied at each cell, so,
Ex,i , j , k = n Ex,i , j , k +
t n +1 2 H z ,i , j , k n +1 2 H z ,i , j 1, k y
n +1 2
n +1
H y ,i , j , k n +1 2 H y ,i , j , k 1 . (98) z
In this notation, note that the spatial offsets from the grid points are not written explicitly, similar to the MRTD case in the previous chapter (but follow the same scheme as MRTD). In this case, the electric fields are located half a cell from point ix, jy, kz in the direction of their field component, and the magnetic fields are located half a cell from ix, jy, kz in the two directions normal to the field component that they represent. A two dimensional cross-section of an FDTD grid intersected by a PEC is presented in Figure 10. In FDTD, the PEC is placed along the locations of the electric field points. The practical result is that the size of the structure being simulated is constrained by the grid. To represent field variation caused by the PECs, several cells are usually placed across each PEC. The PEC condition is applied by first updating the electric fields, and then setting the field values that overlap with the PEC locations to zero. This is possible because the electric field update equations, such as (98), only use previous field values.
47
Figure 10:
FDTD grid intersected by PEC
The time update scheme used in MRTD takes the same form as FDTD (the same basis functions are used for the time discretization). Instead of using image theory, which creates artificial scaling/wavelet coefficients to apply PEC conditions, it would be convenient to explicitly enforce the boundary condition of zero tangential electric field on the MRTD grid. For a general wavelet basis, using the notation presented in the previous chapter, the electric field is reconstructed as,
n
Edir (r ) = T, j , k n E dir ,i , j , k . i
i j k
(99)
This function is used to give the fields at every point in the grid. It is useful to note that the MRTD scheme never updates the field values themselves, as in FDTD, but only updates the scaling/wavelet coefficients that then must be reconstructed to determine field values. This is important to note when probing the fields during simulation, as the fields must be reconstructed at points of interest. When probing fields during simulation, (99) is unwieldy; reconstructing the fields for the entire grid for each time step is
48
computationally prohibitive. reconstructed using
Instead, the fields for an individual cell can be
n
Edir ,i , j , k (r ) =
a = na b = na c = na
nb
nb
nb
T i + a , j + b, k + c n
E dir ,i + a , j + b , k + c .
(100)
where na and nb represent the size of the stencil (the overlap of the wavelet/scaling function in cells). Once the fields are reconstructed, a localized PEC can be applied by multiplying the reconstructed electric field by a function P(r), where,
0 P(r ) = 1
PEC location . Elsewhere
(101)
Once the PEC boundary condition has been applied, the fields must be transformed back into the wavelet domain. This can be accomplished by applying a wavelet transform,
E dir ,i , j , k = i , j , k , P(r )n Edir ,i , j , k (r ) .
(102)
Of course, the new components must be found at all locations whose stencils include a PEC location. By combining (99) and (102), the wavelet coefficient with the PEC condition applied can be found directly from other scaling/wavelet coefficients using
E dir ,i , j , k = i , j , k , P(r ) T, j , k n E dir ,i , j , k . i
i j k
(103)
In the case of a general wavelet resolution, there are several difficulties in applying this technique. The first is that, for arbitrarily placed PECs, the reconstruction, and subsequent wavelet decomposition, of the entire grid is required for each time step. This procedure likely requires more computation time than the field updates, and is therefore not practical. In addition, the implicit assumption in the above procedure is that the fields, when deconstructed back into the wavelet domain after the application of the PEC
49
boundary condition, reconstruct to the same values at all non-PEC locations, and zero at PEC locations. While all of the wavelet bases discussed in this thesis are complete in
L2(R), this condition is only true for an infinite level of wavelet resolution. For a limited
wavelet resolution, which is kept the same before and after the application of the PEC, the application of the PEC as outlined here will result in modified fields outside the PEC region, and nonzero fields inside the PEC region. This is not consistent with Maxwells equations. For this technique to be applied successfully, the wavelet basis must satisfy two conditions: 1. The scaling/wavelet functions for one cell must not overlap with a neighboring cell. 2. The application of a PEC boundary condition (zeroing the field) over a range within an MRTD cell must not affect neighboring field values. The first condition allows the PEC boundary condition to be applied locally. This means that only the scaling/wavelet coefficients in a single cell must be modified to apply the PEC boundary condition. The second ensures that the fields will not undergo nonphysical modifications when the PEC condition is applied. Both of these conditions are satisfied by the Haar wavelet basis.
3.2
PROPERTIES OF HAAR-WAVELET DISCRETIZATIONS
To demonstrate that the Haar basis can be used to explicitly apply subcell PEC boundary conditions it is first useful to study how the Haar basis functions discretize field values. The Haar scaling and wavelet functions are defined in the previous chapter. When the Haar wavelets are reconstructed, they are constantly valued over discrete regions. For example, Figure 3 presents the 2-D Haar scaling function and wavelets for
50
rmax=0 in both directions. When these functions are summed, they yield four independent
regions of constant field value. These regions are centered at the equivalent grid points. As was stated previously, the number of equivalent grid points is the same as the number of basis functions. In the 3-D case, the number of equivalent grid points (and basis functions) is, # of equivalent grid points = 2 which is equivalent to (28) for the 3-D case. As another example, the rmax=1, in all directions, wavelets for the 2-D case are presented in Figure 11. In this case there are 2
2 + rmax, x + rmax, y
3 + rmax, x + rmax, y + rmax, z
,
(104)
= 2 2 +1+1 = 16 coefficients per
cell. In this figure, each of the wavelets has two possible, identical in magnitude, values, with opposite sign. These are represented in the figure by different shading. For the wavelet/wavelet case, these represent four distinct regions. the For highest level
r r wavelet/wavelet terms (the pmax ( x ) pmax ( y ) terms), the equivalent grid points are at the
center of the areas of constant magnitude. Each of the highest level wavelet/wavelet terms represents four equivalent grid points in the 2-D scheme.
51
Figure 11:
Haar scaling and wavelet functions in 2-D, rmax=1
This discussion is easily extended to 3-D. In the rmax=1 case an additional four scaling/wavelet coefficients are required to describe the fields in the z direction. The total number of basis functions in this case is 2
3 + rmax, x + rmax, y + rmax, z
= 23+1+1+1 = 64 . In this case,
the regions of constant value are rectangular solids, and for the wavelet/wavelet terms there are eight areas of constant value. Again, the center of these regions for the highest resolution wavelet/wavelet terms represents the equivalent grid points. When the fields are reconstructed they add to independent values, located at the equivalent grid points. The domain of each value is
L=
D
2
1+ rmax, D
,
(105)
52
in each dimension, D. In the rmax=1 case, the value is constant at each equivalent grid point over a range of x 4 y 4 z 4 . One important note is that, for any cell, the fields can be represented in an equivalent manner by scaling functions alone. In this case the scaling functions must be the size of the equivalent grid points. The advantage of using the wavelets is that the resolution can be varied on a cell by cell basis. In cells containing complex structures, high resolution cells can be used. In surrounding cells, lower resolution can be used. Additionally, the wavelet resolution can be varied during simulation. For complex field variation, high resolution wavelets can be used, for less variation, lower resolution can be used. To apply a similar method with scaling functions only (FDTD), an interpolation scheme must be used to interface the high and low resolution areas. In the MRTD scheme, variable resolution is automatically applied by zeroing high resolution wavelet coefficients when they are not needed. The practicality of using Haar wavelet addition/subtraction as a method of representing fields varying at different rates can be represented analytically and graphically. For example, a 1-D Haar wavelet system with rmax=0 is presented. In this scheme there is only one scaling and one wavelet function. These functions can be used to represent any two values. Figure 12(a) shows the Haar scaling and wavelet function and Figure 12(b) shows an example of a dual valued piecewise constant function that can be represented using these values. The function f(x) in Figure 12(b) has value c in the range (0,0.5) and the value d in the range (0.5,1). The scaling function has value s in the range (0,1) and 0 elsewhere. The wavelet function is valued w0 in the range (0,0.5) and -w0 in the range (0.5,1).
53
Figure 12:
(a) Haar scaling and rmax=0 wavelet, (b) sample function that can be represented using (a)
The magnitudes of the scaling/wavelet functions can be found using the system
s + w0 = c s w0 = d
This system can be easily solved to yield,
(106)
s=
c+d 2
w0 =
cd 2
(107)
The scaling function is the average of the two discrete values of f(x). If the two values are identical, w0=0, and the wavelet coefficient is not needed. In practical simulations, wo can be neglected if it is significantly smaller than the scaling function. This same scheme holds for higher level wavelets. For contrast, a similar example is presented for rmax = 1. In Figure 13 the addition of the level one wavelets (the tails are removed for rmax = 1 to demonstrate that the wavelets represent two independent
54
functions) allows the piecewise constant valued function to have a maximum of four independent values. The wavelets have two values, equal in magnitude and opposite in sign, and for simplicity the magnitude is indicated in the center of each function. In this case a system of four equations can be constructed to determine the scaling/wavelet function magnitudes. The solution of this system,
s= c+d +e+ f 4 cd w0,0 = 2 w0 = c+d e f 4 , e f w0,1 = 2
(108)
shows that the scaling function still represents the average of the function. The sum of the scaling term and the 0th level wavelet represents the average on either half of the domain. If the variation of the field on either half of the domain is small, the high level wavelets can be ignored.
55
Figure 13:
Haar example with (a) scaling through rmax=1 wavelets, (b) an
example function that can be represented with these functions
One other interesting property of Haar wavelet expansions related to the representation of PEC structures can be demonstrated using this example. The values of the scaling/wavelet functions presented in (108) can be determined regardless of the values of f(x). If f(x) is zeroed over one of the ranges, the remainder of the values can still be represented using the scaling/wavelet basis. The values of all of the basis
functions in the cell change, but the values obtained when the functions are summed remain the same. In the two examples that have been presented of Haar wavelet decompositions, the values of the Haar wavelets were found using a system of linear equations. In the general case of a continuous function, the method of moments can be applied. In the special cases that have been presented of piecewise constant functions with constant domains equal to the equivalent grid point dimensions, the method of moments discretization
56
reduces to a system of linear equations. If the values of the function over these constant areas are represented as a vector, F, and the scaling/wavelet function magnitudes are represented as a vector Fw, then they are related by a reconstruction matrix, R, where,
F = RFw .
(109)
This representation is also valid for the two and three dimensional cases. Because there are the same number of equivalent grid points as scaling/wavelet functions, R is square. In the two and three dimensional cases, any ordering of the wavelet coefficients and equivalent grid points can be used, it only affects the positioning of the coefficients in the
R matrix.
R can be generated quickly by examining the contribution of each
wavelet/scaling coefficient to each equivalent grid point. Once the reconstruction matrix is determined, the wavelet transform can be easily performed for arbitrary values of F. A wavelet transform matrix, W is defined,
W = R 1 ,
and
(110)
Fw = WF .
(111)
These matrices provide a quick and easy transition between scaling/wavelet coefficients and allow the PEC boundary condition to be explicitly applied at any equivalent grid point in the Haar MRTD scheme without affecting neighboring field values. The only restriction is that the metals must be the size of an equivalent grid point. By increasing the resolution to the appropriate level, arbitrary structures can be represented.
57
3.3
HAAR SUBCELL PEC APPLICATION
The discussion in the previous section demonstrates that Haar wavelets can be used to apply the explicit PEC method that has been presented. The Haar wavelets are nonoverlapping, and thus a modification of the fields in one cell does not affect neighboring cells. In addition, the reconstruction matrix can be used to quickly determine the field values within one cell, which can be subsequently modified and transformed back into the wavelet domain. The MRTD update algorithm that was presented in the previous chapter can be modified: 1. 2. 3. Determine B fields from E fields Determine H fields from B fields Determine D fields from H fields
4. 5. 6.
7. 8.
Reconstruct D fields in PEC grid locations Zero fields tangential to PECs Transform D fields back to wavelet domain
Determine E fields from D fields Repeat until simulation is complete
The added steps are indicated in boldface. The PEC condition is applied directly to the D fields after they are updated. Mathematically, these steps are relatively simple. First, the fields are reconstructed in the cell where the PEC is to be applied,
Ddir ,i , j , k = RDw , dir ,i , j , k .
(112)
In this case the subscripts indicate the fields at cell i,j,k. The vector Ddir ,i , j ,k contains the reconstructed field values in the cell. The entries in this vector that correspond to positions of PECs are zeroed. This can also be represented mathematically if a matrix IP
58
is defined, where IP is the identity matrix with the rows that correspond to the PEC locations zeroed. The application of the PEC to the field becomes,
DPEC, dir ,i , j , k = I P Ddir ,i , j , k .
(113)
To continue the MRTD field updates, the D fields must be transformed back to the wavelet domain,
Dw , dir ,i , j , k = R 1DPEC, dir ,i , j , k .
(114)
The procedure that has been presented requires that the D vector be multiplied by three matrices. These matrices are 2
3 + rmax, x + rmax, y + rmax, z
2
3 + rmax, x + rmax, y + rmax, z
in size, for the
rmax=1 case the matrix is 64x64, containing 4096 entries. While many of the entries in
these matrices are zero, and thus the multiplication can be significantly simplified, this procedure still adds non-negligible computational overhead. There are two methods that can be used to apply this process in a more efficient manner. First, the three steps (112)-(114), can be combined. multiplication, so an alternate method is,
Dw , dir ,i , j , k = Pdir ,i , j , k Dw , dir ,i , j , k = R 1I P RDw , dir ,i , j , k .
Each step is a matrix
(115)
Before simulation begins, the Pdir ,i , j ,k matrices can be calculated, and then the PEC condition can be quickly applied with a single matrix multiplication. This method is a simple way to modify an existing Haar MRTD code to add local PEC modeling. An even more efficient method is to directly modify the MRTD update equations. In Haar MRTD the D field update equations can be calculated using the same procedure as (50), yielding
n +1
D x,i , j , k = n D x,i , j , k
Dx Dx t U H y ,1 n +1 2 H y ,i , j , k +1 + U H y , 0 n +1 2 H y ,i , j , k , + xyz + U D x ,1 n +1 2 H z ,i , j +1, k + U D x , 0 n +1 2 H z ,i , j , k Hz Hz
(116)
59
n +1
D y ,i , j , k = n D y ,i , j , k
Dy Dy t U H x ,1 n +1 2 H x ,i , j , k +1 + U H x , 0 n +1 2 H x ,i , j , k , + xyz + U D y ,1 n +1 2 H z ,i +1, j , k + U D y , 0 n +1 2 H z ,i , j , k Hz Hz Dz Dz t U H x ,1 n +1 2 H x ,i , j +1, k + U H x , 0 n +1 2 H x ,i , j , k . + xyz + U D zy ,1 n +1 2 H y ,i +1, j , k + U D zy , 0 n +1 2 H y ,i , j , k H H
(117)
n +1
D z ,i , j , k = n D z ,i , j , k
(118)
In the Haar MRTD case, the sums are not required, because the overlap of the D and B basis functions only extends to the nearest neighbors. All of these equations have the same format; the only differences between the equations are the field values involved and the entries of the U matrices. It should be noted that the U matrices are the same size as the reconstruction/wavelet decomposition matrices. Neglecting unneeded subscripts, any of (116)-(118) can be represented as,
D W dir1 = D W dir1 +
t U1H W dir 2,1 + U 2 H W dir 2, 2 , xyz + U 3H W dir 3,1 + U 4 H W dir 3, 2
(119)
where the w subscript indicates that the vector is scaling/wavelet magnitudes. The updated D vector can be reconstructed to give field values at specific points by multiplication with R,
Ddir1 = RD W dir1 = RD W dir1 +
t RU1H W dir 2,1 + RU 2 H W dir 2, 2 . xyz + RU 3H W dir 3,1 + RU 4 H W dir 3, 2
(120)
The D dir1 = RD W dir1 term represents the field values at the equivalent grid points, not the scaling/wavelet coefficients. Thus, (120) is an equation that gives updated D field values from H scaling/wavelet coefficients. If the PEC condition is applied,
Ddir1, PEC = I P RDWdir1 = I P RDWdir1 +
t I P RU1H Wdir 2,1 + I P RU 2 H Wdir 2, 2 .(121) xyz + I P RU 3H Wdir 3,1 + I P RU 4 H Wdir 3, 2
Transforming back into the wavelet domain,
60
D Wdir1, PEC = R 1I P RDWdir1 =
1 1 t R I P RU1H Wdir 2,1 + R I P RU 2 H Wdir 2, 2 . R I P RD Wdir1 + xyz + R 1I P RU 3H Wdir 3,1 + R 1I P RU 4 H Wdir 3, 2 1
(122)
If a new matrix UPEC is defined,
U PEC = R 1I P RU ,
then a new update equation can be defined,
(123)
D Wdir1, PEC = D Wdir , PEC1 +
t U PEC1H Wdir 2,1 + U PEC 2 H Wdir 2, 2 . xyz + U PEC3H Wdir 3,1 + U PEC 4 H Wdir 3, 2
(124)
The R 1I P RD Wdir1 condition does not have to be explicitly enforced because D Wdir1 represents the scaling/wavelet coefficients from the previous time step, where the condition has already been applied. This technique provides a method that can be used to automatically apply the PEC boundary condition within the Haar-MRTD update equations. In this method, the U PEC matrices can be calculated before the simulation begins, and therefore do not add significant overhead to the field updates. In fact, this method can be expanded to allow other subcell effects to be automatically applied in the MRTD update scheme.
3.3
GENERAL SUBCELL EFFECTS IN HAAR-MRTD: COMPOSITE CELLS
In the previous section, the reconstruction/wavelet transform matrices were used to apply the PEC boundary condition directly in the MRTD update. reconstruction/deconstruction is applied to all of the fields in the update equation,
1 1 t RU1R RH W dir 2,1 + RU 2 R RH W dir 2, 2 Ddir1 = RD W dir1 = RD W dir1 + ,(125) xyz + RU 3R 1RH W dir 3,1 + RU 4 R 1RH W dir 3, 2
If the
61
and the substitutions,
H = RH W , U L = RUR 1 ,
(126) (127)
where the subscript L is used to refer to the update equation for local, pointwise, fields, is made, then
Ddir1 = Ddir1 +
t U L1H dir 2,1 + U L 2 H dir 2, 2 . xyz + U L 3H dir 3,1 + U L 4 H dir 3, 2
(128)
This equation is an update for the field values at the equivalent grid points from field values at equivalent grid points. It has been repeatedly stated that MRTD and FDTD are equivalent schemes, and (128) is the conversion from MRTD to FDTD. It should be noted, however, that the MRTD scheme can still be used to vary the wavelet resolution, and thus the number of equivalent grid points, on a cell by cell basis. This is an inherent property of MRTD that is not available in FDTD. The pointwise field update
representation, (128), provides the ability to manipulate the fields at individual equivalent grid points and then transform the fields back to the scaling/wavelet domain. One application of this technique, and indeed the motivation for discovering this technique, is the application of the PEC boundary condition at individual equivalent grid points that has already been presented. Another simple application of this technique is the addition of a current source at equivalent grid points. To this point in this thesis, no method of applying a source to the MRTD grid has been presented. This method is equivalent to wavelet transforming a spatial source condition. In the pointwise update scheme, a source, J, can be added at each equivalent grid point.
Ddir1 = Ddir1 +
t U L1H dir 2,1 + U L 2 H dir 2, 2 +J. xyz + U L 3H dir 3,1 + U L 4 H dir 3, 2
(129)
62
Multiplying (129) by R-1 converts the equation back to the wavelet domain, yielding,
D W dir1 = R 1Ddir1 +
U L1H dir 2,1 + U L 2 H dir 2, 2 t 1 R 1 +R J xyz + U L 3H dir 3,1 + U L 4 H dir 3, 2
t U1H W dir 2,1 + U 2 H W dir 2, 2 1 = D W dir1 + +R J xyz + U 3H W dir 3,1 + U 4 H W dir 3, 2
.
(130)
This equation can be used to apply an arbitrarily placed source into the MRTD cell. One of the advantages of the Yee-FDTD scheme is that its popularity and longevity has led to the development of a large number of techniques that can extend its use. Some examples of techniques that have been developed for Yee-FDTD are the modeling of thin wires [24], narrow slots [25], curved structures (with a locally conformal grid) [26], thin material sheets [27], dispersive surfaces (such as thin metals) [28], SPICE circuits [29], local field correction [30], and lumped elements [31]. The result of all of these
techniques is modified FDTD update equations for the areas where the effect is applied. The surrounding fields are updated normally. Using the technique that has been
presented here, it is possible directly bridge the modified FDTD update equations to the MRTD technique. An example of how this technique can be applied to bridge the pointwise FDTD modifications to the MRTD technique is presented here with lumped elements. The procedure for representing lumped elements in FDTD is presented in [31], and a brief overview is given here. A lumped element, here specified to be a resistor, capacitor, or inductor, can be represented in Amperes law as a current source, JL,
H (t ) =
D(t ) + J L (t ) . t
(131)
In this case currents due to loss and impressed currents are neglected for simplicity. For each of the lumped elements listed above, a relationship can be determined for the current
63
as a function of voltage. The Dz component update will be shown here, although this procedure can be easily modified for any direction. current, the z component of (131) is, Including the lumped element
Dz H y H x I = z , t x y xy
noting that,
(132)
J L, z =
Iz . xy
(133)
Iz can be determined for each element. For resistors, I z ,R = V 1 = Ez , R R z
(134)
capacitors,
V =C =C t Ez
z
I z ,C
t
.
(135)
and inductors,
I z ,C =
1 1 Vt = I z Ezt . I
(136)
For any of these cases, the MRTD update equations can be determined by inserting the current relationship into Amperes law and applying the method of moments. However, this method is somewhat difficult. An example is presented for the resistor. The voltage across a single equivalent grid point can be simply calculated as
Vi , j ,k = Ei , j ,k z .
(137)
64
However, for Amperes law, the update is performed at the n+1/2 time step, while the E field is only known for the n and n+1 time steps. Similar to the method that is used for Ohmic losses, the semi-implicit approximation can be applied. In this case,
Vi , j , k =
n +1
Ei , j , k + n Ei , j , k z . 2
(138)
When the method of moments is applied to (132), with the voltage-current relationship (138), the resulting update equation is,
n +1
D z , i , j , k = n D z , i , j , k + L ( n +1 E z , i , j , k + n E z , i , j , k )
Dz Dz t U H x ,1 n +1 2 H x ,i , j +1, k + U H x ,0 n +1 2 H x ,i , j , k . + xyz + U D zy ,1 n +1 2 H y ,i +1, j , k + U D zy , 0 n +1 2 H y ,i , j , k H H
(139)
The entries of L have the form,
Lij = L(r )i j r ,
r
(140)
with
z L(r ) = 2 Rxy 0
resistor location elsewhere
,
(141)
if x, y, and z are the size of the equivalent grid points.
Two difficulties are
immediately apparent with the equation. First, it includes both D and E coefficients, where before E coefficients were updated directly from the D fields. Secondly, D fields are updated using E fields at the same time step. The first of these difficulties can be easily accounted for using the constitutive relationship. In (57), the electric field is determined from the D field. In the Haar case, this relationship becomes
65
n
E d ,i , j , k = U E dd ,i , j , k , a ,b , c n Dd ,i , j , k = d ,i , j , k n Dd ,i , j , k D
(142)
Because no neighboring field values are required to determine the E coefficients from the
D coefficients, the transformation can be applied directly to the update equation. In this
case,
n +1
E z ,i , j , k = d ,i , j , k n +1 D z ,i , j , k = d ,i , j , k n D z ,i , j , k + U H x ,1 n +1 2 H x ,i , j +1, k + U H x , 0 n +1 2 H x ,i , j , k t d ,i , j , k xyz + U H y ,1 n +1 2 H y ,i +1, j , k + U H y ,0 n +1 2 H y ,i , j , k U H x ,1 n +1 2 H x ,i , j +1, k + U H x , 0 n +1 2 H x ,i , j , k t d ,i , j , k xyz + U H y ,1 n +1 2 H y ,i +1, j , k + U H y , 0 n +1 2 H y ,i , j , k
. (143)
= E z ,i , j , k +
Using this update, (139) becomes,
n +1
E z , i , j , k = n E z , i , j , k + L ( n +1 E z , i , j , k + n E z , i , j , k )
+ U H x ,1 n +1 2 H x ,i , j +1, k + U H x ,0 n +1 2 H x ,i , j , k . t xyz + U H y ,1 n +1 2 H y ,i +1, j , k + U H y , 0 n +1 2 H y ,i , j , k
(144)
Collecting terms,
( I L ) n +1 E z , i , j , k = ( I + L ) n E z , i , j , k + U H x ,1 n +1 2 H x ,i , j +1, k + U H x ,0 n +1 2 H x ,i , j , k , t xyz + U H y ,1 n +1 2 H y ,i +1, j , k + U H y ,0 n +1 2 H y ,i , j , k
(145)
and solving for
n +1
n +1
E x ,i , j , k ,
E z ,i , j , k = (I L) 1 (I + L) n E z ,i , j , k + U H x ,1 n +1 2 H x ,i , j +1, k + U H x ,0 n +1 2 H x ,i , j , k . t (I L) 1 xyz + U H y ,1 n +1 2 H y ,i +1, j , k + U H y , 0 n +1 2 H y ,i , j , k
(146)
A similar procedure can be used to determine the updates for capacitors and inductors. This procedure is not difficult, and most of the operations involved are matrix
66
multiplications that can be performed before the start of simulation. The only integration required is to determine the L matrix. However, there is another method that can be used in this case that removes the need to calculate the L matrix. The Ez update equation in FDTD for a resistor is [31],
n +1
Ez ,i , j , k
tz 1 2 Rxy = E tz n z ,i , j , k 1+ 2 Rxy . n +1 / 2 H y ,i +1, j , k n +1 / 2 H y ,i , j , k t x + tz n +1 / 2 H x ,i , j +1, k n +1 / 2 H x ,i , j , k 1 + 2 Rxy y
(147)
Using (128), two additional lumped element matrices can be created that apply the lumped element formulation in (147). The two lumped element matrices L1 and L2, are diagonal, with
tz 1 2 Rxy resistor equivalent grid points = 1 + tz , 2 Rxy 1 elsewhere 1 resistor equivalent grid points . = 1 + tz 2 Rxy 1 elsewhere
L1,i ,i
(148)
L2,i ,i
(149)
L1 and L2 can be inserted into (128), Ddir1 = L1Ddir1 + U L1H dir 2,1 + U L 2 H dir 2, 2 t L2 , xyz + U L 3H dir 3,1 + U L 4 H dir 3, 2
(150)
67
where again, for simplicity, the subscripts have been dropped. By rotating the values of
x, y, and z, L can be used for a resistor in any direction. Converting back to the
wavelet domain,
D W dir1 = R 1L1R 1D W dir1 +
U L1RH W dir 2,1 + U L 2 RH W dir 2, 2 t R 1L 2 . (151) + U L 3RH W dir 3,1 + U L 4 RH W dir 3, 2 xyz
These matrices can be combined so that each field is only multiplied by one matrix for each update,
D W dir1 = U R D W dir1 +
t U R1H W dir 2,1 + U R 2 H W dir 2, 2 . xyz + U R 3H W dir 3,1 + U R 4 H W dir 3, 2
(152)
Using this method, resistors can be inserted into the MRTD grid at any equivalent grid point, and modified update equations can be found before simulation through simple matrix multiplications. One other important note is that this formulation introduced a new update matrix for the current D field. This does add one extra matrix multiplication to the field update, but, as is shown in the next chapter, the general MRTD update including UPML includes this matrix. Similar updates can be derived for capacitors and inductors. modeled in FDTD using, Capacitors are
n +1
Ez ,i , j , k = n Ez ,i , j , k
n +1 / 2 H y ,i +1, j , k n +1 / 2 H y ,i , j , k t x , + Cz n +1 / 2 H x ,i , j +1, k n +1 / 2 H x ,i , j , k 1 + xy y
(153)
so, in this case, only one L matrix needs to be determined,
68
L2,i ,i
1 capacitor equivalent grid points = 1 + Cz . xy 1 elsewhere
(154)
In FDTD, inductors are represented using,
n +1
Ez ,i , j , k = n E z ,i , j , k
n +1 / 2 H y ,i +1, j , k n +1 / 2 H y ,i , j , k z (t )2 n t x + E , (155) n +1 / 2 H x ,i , j +1, k n +1 / 2 H x ,i , j , k Lxy m =1 m z ,i , j , k y
where it is noted that the sum can be represented as a single value that is augmented at each time step. In this case, both L1 and L2 are the identity matrix, but another matrix must be created. In this case, the sum term is similar to the previously presented current source. If (128) is modified
Ddir1 = Ddir1 +
n t U L1H dir 2,1 + U L 2 H dir 2, 2 + L D . xyz + U L 3H dir 3,1 + U L 4 H dir 3, 2 m=0
(156)
with,
z (t )2 Li ,i = Lxy 0
inductor equivalent grid points elsewhere
.
(157)
Then, using the same procedure as above, (156) can be converted into the wavelet domain. In the case of the inductor, an additional vector,
m=0
D
n
W
must be used to
calculate the updates. It is noted that it is only necessary to store this vector in cells that contain inductors.
69
CHAPTER 4 FULL WAVE HAAR-MRTD WITH COMPOSITE-CELL MODLELING, UPML, VARIABLE GRIDDING, AND TIME/SPACE ADAPTIVE GRIDDING
In the previous chapter, a method is presented that allows a variety of elements that are smaller than a single MRTD cell to be simulated using the Haar MRTD method. Unlike the method presented in [22], this technique treats the entire space as an MRTD grid, using scaling and wavelet functions exclusively to represent the field. The method was created to model subcell PEC structures, as most microwave structures are simulated using only PEC and dielectric media. In addition, it is shown that the method can be generalized to allow the modeling of several other subcell effects; a specific example is provided for lumped elements. The MRTD method that includes subcell elements and varying dielectrics within individual MRTD cells is termed the composite-cell method. This is because, instead of homogenous cells that are used in most MRTD simulations, these cells can include complex structures. The advantage of this method is that these subcell structures can be treated within the MRTD framework, allowing the adaptive MRTD grid to be used to effectively model these structures. In order to verify the method, and show that the technique can be used to apply an adaptive grid to realistic structures, a 3-D Haar MRTD code was written. This code implements the 3-D MRTD scheme presented in Chapter 2, along with the subcell technique for PECs, sources, and lumped elements presented in Chapter 3. Simulations that were performed using this code are presented in Chapter 5. In this chapter, the specific elements of the code that was generated that have not before been published are presented. These elements include a UPML that is valid for arbitrary wavelet resolution, non-uniform gridding, and adaptive gridding.
70
4.1
ARBITRARY WAVELET RESOLUTION UPML
The UPML [19] was first developed as an alternative to the Berenger PML [18] that did not require non-physical field splitting. The UPML is expressed as a material with a carefully designed permittivity and permeability tensor [14],
s y sz sx ( ) = ( ) = 0 0 0 sx sz sy 0 0 0 , sx s y sz
(158)
where,
sn = k n +
n j
n = x, y , z .
(159)
The parameter n represents the loss inside the UPML, and k n is a matching parameter that can be used to fine tune the UPML performance. The permittivity and permeability tensors inside the UPML are frequency dependent, and are therefore difficult to represent in a time domain scheme. If the constitutive relationship
Dx ( ) =
sz E x ( ) sx
D y ( ) =
sx E y ( ) sy
Dz ( ) =
sy sz
E z ( ) .
(160)
is defined, then a coupled set of time domain equations can be defined to update the electric and magnetic fields. The derivation for the electric fields is shown here, the magnetic field updates are found in the same manner. The time domain differential equation for the D field is,
71
H z H y y z k y H x H z = 0 z x H y H x 0 y x
0 kz 0
0 Dx y 0 0 Dx D + 1 0 0 Dy . 0 y z t 0 0 z Dz k z Dz
(161)
This system is very similar to Amperes law when Ohmic loss is included (the only difference is the constant that precedes the time derivative). It can be discretized using the procedure in Chapter 2. For example, the Dx component of (161) is
H z H y Dx y = ky + D , y z t x
and the time localized form is,
(162)
n +1
Dx = n Dx
t y n +1 Dx + n Dx t n +1 / 2 H z n +1 / 2 H y . + z k y 2 k y y
(163)
Combining D terms,
2k t y 2k yt n +1 / 2 H z n +1 / 2 H y + . Dx = n Dx y n +1 2k + t 2k + t y z y y y y
(164)
Next, space localization can be performed, and the update matrices can be computed. The Dx update is,
n +1
D x ,i , j , k = U D xx ,, n D x ,i , j , k
D
U D xy ,1 n +1 2 H y ,i , j , k +1 + U D x , 0 n +1 2 H y ,i , j , k . H 1 Hy + xyz + U D x ,1 n +1 2 H z ,i , j +1, k + U D x , 0 n +1 2 H z ,i , j , k Hz Hz
where,
D U Dxx,i , j =
(165)
Dx
2k t y , i , y 2k + t D x j y y
(166)
72
D U H yx ,m,i , j =
Dx
2k yt i , 2k + t z y y 2k yt i , 2k + t y y y
Hy
j
m
,
(167)
D U H zx ,m,i , j =
Dx
Hz
j
m
.
(168)
Similar equations can be found for the Dy and Dz components. The relationship between the E and D fields inside the PML is also more complicated than in the isotropic non-dispersive case that has already been presented. When the Dx relationship of (160) is converted to the time domain [19],
(k x Dx ) + x Dx = (k z Ex ) + z Ex , t t
must be solved to determine the updated E fields from D fields. The time localized form of (169) is
(169)
k x ( n +1 Dx n Dx ) +
t x
n +1
Dx + n Dx t z = k z ( n +1 Ex n Ex ) + 2 Ex
n +1
Ex + n Ex . (170) 2
Collecting terms and solving for
n +1
n +1
2k z z t (2k x + x t )n +1 Dx 1 Ex = 2k + t n Ex + (2k + t ) (2k t ) D . x x n x z z z z
(171)
Solving for the Ex coefficients,
n+1
E x ,i , j ,k = U Exx ,, n E x ,i , j ,k + U Exx + n+1 D x ,i , j ,k + U Ex , 0 n D x ,i , j ,k . D D
x
E
[
]
(172)
Where,
E U E xxi , j =
Ex
2k z z t i , 2k + t E x j , z z
(173)
73
E U D xx+ i , j =
Ex
1 2k x + x t , i , 2k z + z t E x j
(174)
E U D xxi , j =
Ex
1 2k x x t , i , 2k z + z t E x j
(175)
and, again, similar equations can be found for the other directions. A process similar to the one presented here can be used to find the updates for the magnetic fields. It has been left until this point to discuss the values of and k that should be used with this scheme. The technique for implementing the UPML that has been presented here is derived in an equivalent manner to the FDTD UPML [19] (until the last step, where the updates for the scaling/wavelet coefficients are found). The values of and k can be chosen in exactly the same manner as an FDTD implementation. In this case,
n(n), where n is x, y, or z, is only nonzero within a predetermined distance d of the n
normal outer boundaries. Usually, a width of 10 cells is sufficient for the UPML. Inside this boundary, n(n) is varied from zero to a maximum value max. A polynomial grading works well,
m n (n) = n d n, max .
( )
(176)
Similiarly, k can be varied from 1 in the non-UPML region, to kmax at the outer edge of the grid. There is one important difference between the FDTD and MRTD implementations of the UPML. In FDTD, the values are discretized and only applied at FDTD grid points. In MRTD, however, the function n(n) is used when calculating the update equations, and thus the variation in n(n) across a cell is accounted for. This method could also be applied to FDTD, by determining FDTD coefficients as in S-MRTD. It is important to note the similarity between the field updates for the PML and the field updates for fields outside the PML. If =0 and k=1, the equations are the same as
74
for the PML free case. For implementation of the PML in a computer, it is convenient to use the PML formulation everywhere, and to only use nonzero values inside the PML region. As an added convenience, by using this formulation, the PML does not need to cover an entire cell. For example, in the rmax=1 case, there are four equivalent grid points per direction. If at least 10 equivalent cells of PML coverage are required, then this scheme permits exactly 10 to be used, instead of 12 (for three cells entirely covered with PML).
4.2
NON-UNIFORM GRID IN MRTD
In the FDTD method, one common technique for conformal meshing is nonuniform gridding [32]. For a fixed grid size, it is easy to pick a structure (for example parallel transmission lines) where a uniform grid (constant x, y, z) cannot align with the structure at all points. If, instead, the grid spacing can vary with position, many more structures can be accurately represented in the FDTD grid. Of course, curved structures must still be stair-stepped, but even with these structures, a more accurate grid using fewer grid points is possible. The FDTD algorithm depends on neighboring cells having identical dimensions at intersections, and if the grid is varied so that the grid size is a function of all three coordinate directions, the grid becomes incompatible with FDTD. Therefore, FDTD non-uniform grids must vary the spacing as a function of individual directions only. This means,
x = x(i ) y = y ( j ) z = z ( k ) .
(177)
One example of a grid that uses this scheme is presented in Figure 14. The 2-D grid that is presented shows the positions of the Ex, Ey, and Hz fields. The value of x or y must be computed based on the field update. For the space derivatives
75
of fields, the value of x or y is the separation between the field components. Thus, for electric field space derivatives, the x or y value is simply the grid spacing. However, for magnetic field spatial derivatives, the value is instead the average of the values for the cells containing each magnetic field. For example, if the Hz coefficient centered at
x3,y3 in Figure 14 is updated, the Ey field coefficients required are the ones centered at
76
Figure 14:
2-D FDTD non-uniform grid example
77
y3 and located on either end of x3. For this update, y= y3. But, when the Ey field
component centered at y3 and between x3 and x4 is updated, the Hz fields required are those centered at x3 ,y3, and x4 ,y3. In this case, x=( x3 + x4)/2. It would be convenient to have a similar capability for MRTD. While the MRTD adaptive grid does allow the resolution to vary on a cell-by-cell basis, the equivalent grid size is always
nequiv =
n . 21+ rn
(178)
Using this grid size, it is difficult to align the grid with realistic structures. Furthermore, increasing the resolution until an equivalent grid point is close to a feature being modeled is an inefficient use of the adaptive grid. A more effective grid would use both nonuniform and adaptive gridding. The most difficult aspect of implementing a non-uniform grid in the MRTD method is determining the correct x, y, and z values to use for each field. In MRTD, the cell size does not appear explicitly in the field update (for the formulation presented here) but is an implicit part of determining the update coefficients. In the derivation presented in Chapter 2, a formula for the offset between the electric and magnetic fields is given,
sd =
2
d rd ,max + 2
.
(179)
In this formula, rd,max is the maximum resolution in the given direction. It was discussed that in the general 3-D scheme, the fields being updated (the time derivative terms in Faradays and Amperes laws) are offset from the other fields in the update (the spatial derivative terms) only in the direction of their spatial derivative. If this scheme is applied consistently for all electric and magnetic field coefficients, there is a specific arrangement of grid points that must be used. That arrangement is outlined here.
78
The convention chosen for this work is that all fields are offset relative to a grid that is defined by (177). The field components can be indexed to each of these grid points, but the actual domain of each function, defined by the Haar-MRTD scaling function (this discussion can be easily extended to general wavelet basis by defining the domain as the spacing between the centers of the scaling functions, centered on each scaling function), begins at the offset from the grid point. The grid is consistent with the conditions above if the electric fields are offset by s in their coordinate directions, and the magnetic fields are offset by sd in the two directions normal to their coordinate direction. An example of a grid that meets these criteria is demonstrated in Figure 15. It should be noted that while the offset can be made relative to another point (for example the magnetic fields could be offset in their coordinate directions and the electric field in the normals) the resulting arrangement of the fields relative to each other must be maintained. The offset of (177) denotes the start of the scaling function for the electric or magnetic field, but does not specify the size of the scaling function. If the size of each scaling function is chosen to simply as x(i ) y ( j ) z (k ) , then the basis functions for a single field can overlap. In Figure 15, a 2-D view of an MRTD grid for rmax=0 (in both directions) is used to demonstrate variable gridding. For rmax=0, sx=sy= . In the figure, the positions of the Ez and Hz equivalent grid points, and the domain of the scaling functions, are shown. The Ez field is offset in the z direction (normal to the page) while the Hz field is offset in the x and y directions.
79
Figure 15:
Offset between electric and magnetic fields in MRTD
(a) fixed grid (b) non-uniform grid (implemented incorrectly) (c) non-uniform grid (implemented correctly)
80
Figure 15(a) demonstrates a uniform MRTD grid.
In this case, both scaling
functions have dimensions x y . The scaling functions for each field start at the same point relative to the grid, and form a non-overlapping basis covering the entire space. Figure 15(b) demonstrates a non-uniform MRTD grid. In this case, the x1 is twice x2 . The Ez field is not offset in the x direction, and thus for each cell x = xi . The Hz field, however is offset in the x direction. When x = xi is used for this field, the Hz field for the first cell continues to the center of the neighboring cell. The Hz field for the second cell, however, begins one quarter of the cell width from x1. This causes the field
coefficients that represent the Hz field to overlap. While it may be possible to use this field arrangement for MRTD, the overlap between neighboring field components would render the scheme implicit, and the wavelet basis non-orthogonal. In order to keep the wavelet scheme non-overlapping, and still allow for a variable grid, the width of the offset scaling functions must be set so that it begins sd from the current grid point, and ends s from the next grid point (note that sd is proportional to the cell spacing). This is demonstrated in Figure 15(c). A scheme for determining the grid size for the different basis functions can then be quantified: 1. For non-offset directions (directions normal to field component for electric fields, or direction of field component for magnetic fields),
n = n(l ) ,
(180)
where l is an index denoting the position along the axis 2. For offset directions (direction of field component for electric fields, or normal to field component for magnetic fields),
n = n(l ) + n(l + 1) n(l ) . 4
(181)
81
4.3
MRTD GRID EXCITATION
In Chapter 3, a method was presented that allows the addition of an impressed field to the MRTD field update equations. A modified Dx update was presented, (129), that allows a source function to be applied at specific equivalent grid points and then be transformed back into the wavelet domain. While this equation can be used to impress a field at any (or every) point in the MRTD grid, it is likely that only a relatively small number of grid points will be modified to excite the grid. In this case, the fields can be excited using
D W dir1 = D W dir1 + R 1J ,
(182)
after the fields are updated, if J is a vector of field values that is only nonzero at desired equivalent grid points. In this case, (129) is simply split into two steps. For the cases that were simulated in this research, only a microstrip excitation was required. The microstrip excitation used for this work impressed a constant electric field in a plane normal to the microstrip metal, directly under the microstrip. For any wavelet resolution, it is necessary to use (182). This is because, when wavelets are applied, there is more than one equivalent grid point in each direction. This relatively simple excitation is demonstrated in Figure 16. In this example, rmax=0 in all directions, therefore there are two equivalent grid points per cell in each direction (eight total equivalent grid points). In this example, only the fields directly under the microstrip are excited. The field strength is represented by arrows, demonstrating that the impressed field at all points is identical in magnitude and direction. The arrows are located at the Ey equivalent grid points, and the cells are represented by the alternating shading. This is meant as simple approximation of the microstrip mode; as the field propagates through the grid it quickly matches the microstrip field pattern. In this example, five Ey field points fall under the microstrip
82
line. As there are two of the field points in the y direction per cell, one cell is split in the
x direction, where only half is excited. Similarly, all cells that contain the PEC are split
in the y direction. Although it could be claimed that this is a contrived example, as it would be a relatively simple matter to keep an even number of points across and under the microstrip line, the x-z plane demonstrates why subgridding is required.
Figure 16:
Excitation in MRTD, rmax=0 cells demonstrated by alternate shading
As in all other directions, there are two Ey field points in the z direction per cell. These field points are represented by the dots in Figure 16, and the boundaries of the top metallization are the dark lines. For the excitation, however, only a single set of fields is excited. If both are excited, the fields resemble two identical, slightly offset pulses, needlessly complicating analysis. This is similar to subcell PEC analysis, where PECs can be simulated with an even number of cells in their tangential directions, but without subcell modeling, the PECs are several equivalent cells in thickness. Similar limitations apply to coplanar waveguide (CPW) excitations; this case is presented in Figure 17. For CPW excitations, the field is guided between the center conductor and the surrounding ground planes. In this case, the field is excited as a
83
voltage between the center conductor and the outer ground planes. In the direction of propagation (the z direction in the diagram) the field is only one cell thick. The
impressed voltage along a line one equivalent cell in thickness can also used for wire antennas and probe feeds (in coaxial fed microstrip antennas or waveguides, for example).
Figure 17:
CPW excitation, demonstrating subcell excitation
4.4
TIME/SPACE ADAPTIVE GRIDDING
At the beginning of Chapter 3, the two main advantages of the MRTD method are presented. The first of these advantages is that many of the wavelet basis functions that can be used in MRTD allow the use of fewer basis functions per wavelength than FDTD, resulting in fewer coefficient updates per time step. However, it can be very difficult to apply PEC boundary conditions and the UPML boundary with these cases. In this thesis,
84
the focus is placed on Haar wavelet basis functions, and while it is true that the Haar cells are larger than FDTD cells, the total number of scaling/wavelet coefficients needed per wavelength is the same as in the FDTD case. The advantage of the Haar scheme is the time-and-space adaptive grid. The methods that have been presented in this thesis that allow for the modeling of general structures using Haar MRTD are derived using an equivalence between FDTD and Haar MRTD. If the Haar MRTD method is applied, and the maximum wavelet resolution is used in all cells, the scheme is equivalent to FDTD. However, the chief advantage of the MRTD method is that the wavelet resolution can be varied from cell to cell. Most simulations are not constrained by the dispersion requirement, but instead by the structures that are simulated. For accurate simulation, several grid points must be used across each feature. Using Haar MRTD, the resolution can be locally increased to model complex structures, and reduced elsewhere for computational efficiency. This technique can be applied statically before the start of simulation, and automatically during simulation. When applied during simulation, this feature is usually referred to as the MRTD time-and-space adaptive grid [33]. The concept of the equivalent grid point was introduced in Chapter 2 to demonstrate the offset of the electric and magnetic fields in MRTD. In Haar MRTD the equivalent grid points are functionally identical to FDTD grid points (regions of constant field value). When wavelets are added and subtracted from a Haar representation, the effect is to add and remove equivalent grid points. This allows dense discretizations to be used in the area of complex structures, and coarse discretizations to be used elsewhere. Figure 18 shows the arrangement of the fields (equivalent grid points) in a single rmax=0 Haar MRTD cell, and also shows the positions of the FDTD grid points if the Yee-cell in Chapter 2 is used. These FDTD equivalent grid points can be used to demonstrate the FDTD equivalent grid for Haar MRTD adaptive resolution.
85
Figure 18:
Haar MRTD cell (2-D), showing FDTD grid points, rmax=0
As an example of the capabilities of the Haar MRTD adaptive grid, consider the grid presented in Figure 19. In this example, three wavelet resolutions are used from
rmax=-1 (scaling function only) to rmax=1. In the rmax=-1 case, there is one equivalent
FDTD cell per MRTD cell, in the rmax=1 case there are eight equivalent FDTD cells per MRTD cell (for a 2-D example). Similar to Figure 18, the dots represent the equivalent FDTD cells, and the surrounding squares represent the MRTD grid boundaries. In Figure 19, the adaptive grid is used to place a dense discretization in the center, and a coarse discretization in the surrounding area. Using this grid, a highly detailed structure can be modeled in the center, while the surrounding, homogeneous area is simulated with low resolution. This technique is more powerful than the non-uniform mesh that is presented earlier in this chapter because the spacing of the equivalent grid points is a function of all three grid dimensions, and thus can be increased locally to represent complex features. One important note is that, while this example demonstrates the grid when the resolution is the same in all dimensions, this is not a requirement.
86
Figure 19:
Adaptive grid example
While the grid that is presented in Figure 19 can be fixed for the entire simulation, the MRTD method also provides the ability to change the resolution during simulation. In Chapter 3, it is demonstrated that at any point in the grid, the highest resolution wavelet represents the deviation from the average. If the field is not changing rapidly, the high resolution wavelets can be neglected during the simulation. In [33], it is
suggested that an absolute and relative threshold be used to determine whether a wavelet is needed. This concept forms the basis for the technique that was used in this thesis. For this work, the thresholds were chosen as in [33]. The relative threshold is a fraction of the value of the scaling function. If the value of the wavelet is less than this fraction of the scaling function, it can be neglected. The absolute threshold is required for low field values, and does not change during simulation. If the fields are near zero, the difference between the scaling and wavelet coefficients will not be large. However, the wavelets are not required because the field does not vary rapidly. In this case, the wavelets will be lower than the absolute threshold, and can be neglected.
87
This scheme is not effective if it is applied every time step. For the test to be applied, the wavelet value must be calculated, which defeats the purpose of the adaptive grid. In this work a scheme was used where the wavelets were tested at a user defined period. Wavelets that are below either threshold are neglected until the next testing period. Using this scheme, only the highest level wavelets are tested. If the wavelets are removed, at the next testing period both the currently highest level wavelets are tested (one lower than the previous level) and the next highest level is tested. In this manner, it can be determined if wavelets need to be reintroduced into the grid. In the code that was used in this simulation, both time-adaptive grids and static (but variable in space) grids are used. Using this scheme, high resolution grids can be used to specify complex structures, while the resolution around the structure can be predefined to a low value. In this manner, the adaptive grid can be applied to vary the resolution in response to field propagation, while the resolution is not higher than necessary at any point in the structure.
88
CHAPTER 5 MRTD SIMULATION EXAMPLES
The techniques that are presented in this thesis can be used to simulate a wide variety of structures in MRTD. The primary contribution of this thesis, the subcell modeling of arbitrary structures, can be used to simulate any structure that can be simulated in FDTD in MRTD. Before this method was developed, PEC structures could only be the size of entire MRTD cells (the boundary condition in that case is applied by zeroing all scaling and wavelet coefficients in the cell). This severely limits the
applicability of MRTD, as the major advantage of MRTD, the time-and-space adaptive grid, cannot be used effectively. If a structure can be represented within a cell, then the cells used in an MRTD simulation can be significantly larger. Using the method
presented in this thesis, any structure that can be simulated in FDTD can be simulated with MRTD. In addition, the MRTD time-and-space adaptive grid can be applied to allow fewer grid points in the MRTD case. In order to verify the method, several simulations were run using the MRTD code that was developed in this investigation. The code utilizes all of the features that have been discussed, including the UPML, non-uniform gridding and lumped element modeling. For the practical cases, the MRTD results were verified by an FDTD code. The FDTD code that was used to verify these results was also written by the author, and its results have been verified against measurement for a variety of structures; these structures have been presented in several journals and conferences [34-36]. To measure the effectiveness of the MRTD code, both simulation time and number of equivalent grid points required to simulate the structure were recorded. The number of equivalent grid points, which is equal to the total number of grid points in the FDTD case, is useful because it directly relates to the number of calculations required. This is
89
the best measure of how one MRTD grid compares to another. In the cases where the time-adaptive grid was used, the number of equivalent grid points provided a method to measure the effectiveness of the adaptive scheme. The microstrip structures that were evaluated in the following examples were measured by computing the voltage on the microstrip lines. The voltage is calculated using
Vab = E(r ) l ,
ab
(183)
where a is the ground plane and b is the microstrip conductor. This is the same method that is used for FDTD simulations [14]. Similarly, current is calculated using,
I = H(r ) l ,
C
(184)
where C is the contour that the current passes through. These are relatively simple calculations, and are solved in the discrete case by summing the electric or magnetic field values along the path and multiplying by the space step. Of course, in MRTD, the field must be reconstructed to perform the sum, ...