THEORETICAL MODELING OF MMIC'S USING WAVELETS, PARALLEL COMPUTING AND A HYBRID MOM/FEM TECHNIQUE by Jui-Ching Cheng A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1998 Doctoral Committee: Professor Linda P. B. Katehi, Chairperson Associate Professor Hal G. Marshall Associate Professor Kamal Sarabandi Professor John Volakis RL-959 = RL-959

ABSTRACT THEORETICAL MODELING OF MMIC'S USING WAVELETS, PARALLEL COMPUTING AND A HYBRID MOM/FEM TECHNIQUE by Jui-Ching Cheng Chairperson: Professor Linda P. B. Katehi This dissertation explores three subjects which contribute to increasing the efficiency, in terms of speed and storage, of numerical methods for electromagnetic problems. The hybrid MoM/FEM technique is first applied to the design of a microstriplinefed slot-coupled cavity-backed patch antenna, a coaxial line-fed cavity-backed patch antenna, a wide bandwidth microstripline vertical cavity coupler, and a high-Q cavity resonator. The latter was fabricated and tested to verify the theoretical results. The second is a wavelet based method of moments. Instead of ordinary basis functions, a class of bases called wavelets are used in the formulation of MoM. These wavelets are locally oscillatory and have zero average. Hence, when integrated with other smooth functions, the resulting value will be very small. When used in MoM implementations, the resulting impedance matrix will have many small

1 off-diagonal terms which can be set to zero without adversely effecting the solution. Using sparse storage techniques, substantial memory is saved by using wavelet bases. These wavelet bases were applied to the analysis of patch antenna arrays, leading to large memory savings. Parallel computers were applied to speed up the computation process. The idea of parallelization is to divide a single task to many independent jobs. The communication between the jobs must be minimal to avoid the overhead incurred by the slower (compared to the speed of CPU) communication link. The matrix filling process and the fast wavelet transform (FWT) were two processes suitable for parallelization. These processes were implemented in the C++ language using a portable and standard Message Passing Interface.

ACKNO'WLED GEMENTS ii

I came to Ann Arbor in the September of 1992. I never thought I would stay here so long. I must thank all the friends I made here during the passed five years. Without them, my life would be dull and gloomy. I am deeply grateful to my advisor, Professor Linda Katehi. Without her long term support and guidance, this work would never be finished. I also admire her enthusiasm for research and education, and kindness to her students. I would also like to thank my committee member, Professor Hal Marshall, Professor Kamal Sarabandi, and Professor John Volakis for their insightful suggestions on my dissertation. I also want to express my gratitude to many former and current members of Radiation Laboratory, especially Nihad Dib, Kazem Sabet, and Jong-Gwan Yook for their technical help, and Chi-Cheng Yu, Yuei-An Liou, Eric Li, Yi-Cheng Lin and Tsen-Chieh Chiu for their companionship. Finally I want to pay my respect to Professor C. T. Tai for his accomplishments in electromagnetics and his hospitality to Chinese students. iii

TABLE OF CONTENTS ACKNOWLEDGEMENTS....... LIST OF FIGURES............ LIST OF TABLES............. LIST OF APPENDICES......... CHAPTER....... ii vi xiii Xlll... xiv I. Introduction.. 1.1 1.2 1.3 Motivation and Objectives.. Overview... Background........................ 1................. 1........................... 3................. 5....................... 8 II. Hybrid MoM/FEM Techniques 2.1 2.2 2.3 Introduction................. Method of Moments Formulation...... FEM Formulation.............. 8 9 12 15 2.4 Combining the solutions of MoM and FEM.. III. Cavity-backed Patch Antennas........... 3.1 Introduction................... 3.2 Spectral Domain Green's Functions....... 3.3 FEM Mesh Generation............. 3.4 Slot-Coupled Microstripline-fed Patch Antennas 3.4.1 Formulation.............. 3.4.2 Numerical Results........... 3.4.3 Numerical Considerations....... 3.5 Coaxial Line-fed Patch Antennas........ 3.5.1 Formulation.............. 3.5.2 Numerical Results........... 3.6 Conclusion.................... 16 16 18 21 25 27 30 32 39 39 41 43 iv

IV. Microstripline-fed Cavity Couplers and Resonators.. 4.1 Introduction...................... 4.2 Single Cavity Formulation.................. 4.3 Multi-cavity Formulation.................. 4.4 Validation............................ 4.5 Wide-band Cavity Couplers/Vertical Transitions..... 4.6 High-Q Micromachined Cavities.............. 4.6.1 Design Consideration............... 4.6.2 A High-Q Resonator............... 4.7 Conclusion...................... V. Wavelet-Based Method of Moments.............. 46 46 46 48 54 60 66 68 70 73 82 5.1 Introduction......... 5.2 Multiresolution Expansion.. 5.3 Moment Method Formulation 5.4 Fast Wavelet Transform... 5.5 Microstrip Patch Antenna.. 5.6 Patch Antenna Arrays.... 5.6.1 3 x 3 Array.. 5.6.2 6 x 6 Array.. 5.7 Conclusion.......... *..... *..,... of a Patch *..... *..... *..... *..... *..... *..... *.... Antenna ~.... ~.... ~...,. *... * * * * * * * * * * * * * * * 82 83 85 87 92 96 96 97 103 VI. Application of Parallel Computing in MoM............ 111 6.1 6.2 6.3 6.4 Introduction......................... Parallelization Model of MoM Matrix Fill Process Parallelization Model of FWT............ Conclusion...................... 111 112:118:120:126 VII. Conclusion............................... 7.1 Summary of Achievements............. 7.2 Future Work......................... APPENDICES............................ BIBLIOGRAPHY.................................126....... 127...... 128 139 v

LIST OF FIGURES Figure 2.1 A general scattering problem which consists of a cavity, an obstacle and a current source.......................... 9 2.2 The equivalent problem of Fig. 2.1..................... 10 3.1 The structure of a slot-coupled cavity-backed microstrip line-fed patch antenna.............................. 17 3.2 The structure of a coaxial line-fed cavity-backed patch antenna... 18 3.3 A hexahedron with parallel and rectangular top and bottom planes can always be decomposed to distorted bricks as show in Fig. 3.4.. 22 3.4 A distorted brick. Note that only the vertical edges are distorted. The top and bottom planes are still parallel and rectangular..... 22 3.5 Edges on a PEC plane have to be created on both sides.. In this figure, Element 1 consists of el, e2,..., e6, and Element 2 consists of e7, e8, *., el2. Note that the position of e4 to e6 and e10 to el2 are the same........................ 25 3.6 The side and top views of the patch antenna shown in Fig. 3.1.... 27 3.7 The basis functions used on the microstripline............ 29 3.8 This figure shows how the basis functions used in MoM match the tetrahedral elements used in FEM. To simplify, not all of the PWS modes are shown................................. 30 3.9 Comparison of the reflection coefficient of a microstrip line-fed cavitybacked patch antenna calculated by the method presented in this dissertation and the FDTD method. The FDTD data were contributed by Nihad I. Dib[40]....................... 31 vi

3.10 The normalized input impedance of a microstrip line-fed cavitybacked patch antenna. The parameters are the same as in Fig. 3.9. The FDTD data were contributed by Nihad I. Dib[40].......... 32 3.11 Comparison of the normalized input impedance of microstrip linefed cavity-backed patch antennas with different slot lengths. Other parameters are the same as in Fig. 3.9................... 33 3.12 Total E-field intensity in the cavity of Fig. 3.1 at resonant frequency (3.995 GHz). The parameters are the same as Fig. 3.9. (A) and (B) are the cross section views along lines aa and bb shown in the inset, respectively. (C) is the view just beneath the patch. The black and white lines indicate the boundaries of the patch and slot. This figure is drawn according to the scale of the physical dimensions...... 34 3.13 The convergence check of the same case as Fig. 3.9 using the hybrid technique. Nb denotes the number of basis functions on one side of the patch.......................................... 35 3.14 The convergence check of the same case as Fig. 3.9 using the hybrid technique. Ns denotes the number of basis functions on the slot... 35 3.15 The convergence check of the same case as Fig. 3.9 using the hybrid technique. The first number in the parentheses denotes the number of basis functions around the patch and the second denotes the total number of edges in the cavity........................... 36 3.16 The side and top views of the patch antenna shown in Fig. 3.2... 38 3.17 Comparison of the reflection coefficient of a coaxial line-fed cavitybacked patch antenna calculated by the method presented in this dissertation and the FDTD method. The FDTD data were contributed by Nihad I. Dib[40]................................ 42 3.18 The normalized input impedance of a coaxial line-fed cavity-backed patch antenna. The parameters are the same as in Fig. 3.17. The FDTD data were contributed by Nihad I. Dib[40]................ 43 3.19 Comparison of the normalized input impedance of coaxial line-fed cavity-backed patch antennas with different feeding points. Other parameters are the same as in Fig. 3.17................. 44 vii

3.20 Total E-field intensity in the cavity of Fig. 3.2 at resonant frequency (4.190 GHz). The parameters are the same as in Fig. 3.17. (A) and (B) are the cross section views along lines aa and bb shown in the inset, respectively. (C) is the view just beneath the patch. The black lines and white circles indicate the boundaries of the patch and coaxial line. This figure is drawn according to the scale of the physical dimensions............................. 45 4.1 The structure of a microstrip cavity coupler................ 48 4.2 The side view of the cavity coupler show in Fig. 4.1. Note that the slots on the cavity are closed by a PEC with an equivalent magnetic current.................................. 49 4.3 The side view of a multi-cavity coupler.................. 49 4.4 A microstrip line coupler as described in [34]............. 55 4.5 Comparison of the result from [34] and the hybrid technique of this dissertation. In the legend, EXP and MoM are the experimental and computed results of [34], and MoM/FEM is the computed result using the hybrid technique.......................... 56 4.5 Continued from previous page...................... 57 4.6 Comparison of the computed result of the single cavity and multicavity formulations. The solid line is the data computed by the single-cavity model and the dashed line is the data computed from the multi-cavity model................................ 58 4.6 Continued from previous page...................... 59 4.7 (a)the structures of a microstripline vertical coupler. The cavity size is 0.0968 x 0.968 x 0.10. All unit in mm.(b) computed S-parameters. 61 4.8 Sensitivity study of the cavity in Fig. 4.7. Three cavity widths 0.968 mm, 0.138 mm and 0.188 mm are compared. The data shown in (b) is the same as (a) except that the scale of the figure is different... 63 viii

4.9 Sensitivity study of the cavity in Fig. 4.7. Three slot positions are compared. Curve 1: both of the slots are at the center of the cavity. Curve 2: both of the slots are offset by 0.01 mm to the center of the cavity. Curve 3: the same as the second except that the offset of one slot is in the opposite direction to the other. The data shown in (b) is the same as (a) except that the scale of the figure is different... 64 4.10 Optimized insertion loss v.s. cavity height............... 65 4.11 The structures of the proposed micromachined bandpass filters. (a) input and output microstrip lines on the same substrate. (b) input and output microstrip lines on different substrates............ 67 4.12 A square cavity resonator with both feeding slots at the center.... 69 4.13 The IS111 of square cavities. The cavity sizes are shown in the legend. Other parameters are the same as in Fig. 4.12............ 69 4.14 A rectangular cavity resonator with feeding slots separated by a distance................................... 70 4.15 The computed results of the structure shown in Fig. 4.14....... 71 4.16 The z component of the electric field inside the cavity at the first resonant frequency 19.38 GHz. The figure is drawn to scale. From left to right columns are the field profile on x-y plane at the top and the bottom of the cavity, on x-z plane at the center of the slots and the cavity, on y-z plane at the center of the cavity. The dimension of the structure is as in Fig. 4.14.................... 75 4.17 The z component of the electric field inside the cavity at the second resonant frequency 23.92 GHz. Refer to Fig. 4.16 for more detail.. 76 4.18 The z component of the electric field inside the cavity at the third resonant frequency 30.04 GHz. Refer to Fig. 4.16 for more detail.. 77 4.19 The comparison of IS11i of single, two and three cavity configurations. The dimension of the structure is the same as in Fig. 4.14.. 78 4.20 The comparison of the effect that cavity shape has on the coupling between the microstriplines. The dimension of the structure is as in Fig. 4.14. Curve 1: coplanar feedlines, trapezoidal cavity. Curve 2: feedlines on different sides, rectangular cavity............. 79 ix

4.21 A high-Q resonator fabricated by micromachining[41]....... 79 4.22 The z component of the electric field inside the cavity at the first resonant frequency 10.54 GHz. The figure is drawn to scale. The dimension of the structure is as in Fig. 4.21.. 80 4.23 The computed S-parameters of the resonator shown in Fig. 4.21.. 81 4.24 The comparison of measured[41] and computed S-parameters of the resonator in Fig. 4.21 near the first resonance............ 81 5.1 The concept of multiresolution expansion. Sn and Wn are the linear space expanded by n-th level scaling functions and wavelets respectively. The lower half of the figure illustrate the concept of multiresolution expansion by pulse functions and their wavelets........ 84 5.2 The rooftop basis functions used to expand the current on a patch. (a) Profile along the direction of the current. (b) Profile along the transverse direction of the current.................. 87 5.3 (a) zero degree B-spline scaling functions, (b) zero degree B-spline wavelets, also called Haar wavelets, (c) first degree B-spline scaling functions, (d) first degree B-spline wavelets................. 88 5.4 The modified first degree B-spline scaling functions and wavelets. The scaling functions consist of only triangle functions without the two half triangles basis at the ends. The wavelets are the same as the original except the edge wavelets have zero value at the end. Only one edge wavelet is shown in figure since the other one is just a mirror. 89 5.5 The two-dimensional rooftop basis functions and its wavelets[53]. (a) the scaling function formed by the multiplication of the zero and first degree B-spline scaling functions, (b) a wavelet formed by the zero degree scaling function and the first degree wavelets, (c) a wavelet formed by the zero degree wavelet and first degree scaling function, (d) a wavelet formed by the the wavelets of zero degree and first degree................................................ 90 5.6 A microstrip patch antenna...................... 93 5.7 Multiresolution grids at different resolution levels: (a) 16 x 16, (b) 8 x 8, (c) 4 x4.............................. 94 x

5.8 The dominant component of the induced current on the patch antenna. The parameters used in the computation are a = 5.95 cm, b = 5 cm, d = 0.159 cm, and e, = 2.55. The definition of these symbols are shown in Fig. 5.6.......................... 98 5.9 Effect of thresholding on the current distribution after applying different threshold levels: (a) 0.001%, (b) 0.01%, (c) 0.05%) and (d) 0.1%.....................................99 5.10 Structure of the moment matrix generated by the one level multiresolution expansion after thresholded by levels (a) 0.001%, (b) 0.01%, (c) 0.05% and (d) 0.1%......................... 100 5.11 This figure shows the absolute values of the elements of the impedance matrix in Fig. 5.10 before thresholding v.s. the distance between the two basis functions. The figure is categorized by the types of basis functions shown at the right hand side of the figure......... 101 5.11 Continued from the previous page..................... 102 5.12 Comparison of sparsity v.s. threshold levels of one- and two-level multiresolution expansions............................. 104 5.13 Effect of thresholding on the resonant curves of (a) one-level expansion, (b) two-level expansion....................... 105 5.14 A patch antenna array.......................... 106 5.15 Comparison of the sparsity of one- and two-level expansions of the 3 x 3 array................................ 106 5.16 Structure of the sparse matrix of the 3 x 3 array at 0.01% (two-level expansion)............................... 107 5.17 Effect of thresholding on the un-normalized far-field pattern of the 3 x 3 array............................... 108 5.18 Comparison of the sparsity of one- and two-level expansions of the 6 x 6 array................................ 109 5.19 Structure of the sparse matrix of the 6 x 6 array at a threshold level of 0.01% (two-level expansion)...................... 109 xi

5.20 The effect of thresholding on the normalized far-field pattern of the 6 x 6 array................................ 110 6.1 A dynamical allocation scheme for matrix fill process........ 115 6.2 The flow chart of a dynamical allocation scheme for matrix fill process. 116 6.3 The measured turn around time of the impedance matrix fill process. The theoretical result is plotted by curve-fitting the measured data to (6.2)................................. 122 6.4 A static allocation scheme for the FWT process............. 123 6.5 The flow chart of a static allocation scheme for the FWT process.. 124 6.6 The measured turn around time of the FWT process. The theoretical result is plotted by curve-fitting the measured data to (6.6)..... 125 A.1 The node and edge labels of a tetrahedral element. nl, n2, n3 and n4 are the nodes. el, e2,..., e6 are the edges............... 130 B.1 The structure of a conductor-backed dielectric slab......... 132 xii

LIST OF TABLES Table 3.1 Number of nodes of each data set of Fig. 3.15 in each direction... 39 4.1 The meaning of matrix in(4.1) and (4.2)............... 50 4.2 The comparison of the observed resonant frequencies and the frequencies of cavity resonant modes.................... 72 4.3 The comparison of the observed resonant frequencies and the frequencies of cavity resonant modes................... 73 5.1 Comparison of unique elements v.s. expansion levels........ 91 5.2 Comparison of sparsity levels and memory usage at different threshold levels................................. 98 A.1 Labeling system of a tetrahedron................... 130 xiii

LIST OF APPENDICES Appendix A. Tetrahedral Edge Elements....................... 129 B. Green's Functions of Conductor-backed Dielectric Slabs........ 131 C. The Implementation of Parallel Programming in C++......... 133 C.1 Dynamic job allocation..................... 133 C.2 Static Job Allocation...................... 136 xiv

CHAPTER I Introduction 1.1 Motivation and Background Advances in semiconductor technology have enabled the manufacturing of highly integrated MMIC's (Monolithic Microwave Integration Circuits) which can reduce the cost, size and weight or improve the performance of high frequency communication systems. However, in the process of integrating more and more circuit elements in the same package, interference between components becomes a major factor influencing circuit performance. This necessitates full-wave modeling of the structure to accurately predict the behavior of the circuit. Full-wave modeling uses integral or differential equations derived from Maxwell,s equations. These equations are solved by the finite linear space approximation. In this process, the unknown sources or fields in the Maxwell's equations are discretized by a finite set of basis functions. One way of finding the expansion coefficients of these basis functions is by performing a suitable inner product on the integral or differential equations with the same or different set of functions. This operation results in a matrix equation of which the rank is proportional to the number of unknowns. The equation can be solved by direct inverse of the matrices or interactive linear system solvers. Examples of these techniques include the Method of Moments 1

2 (MoM) and the Finite Element Methods (FEM). Another way of solving Maxwell's equations is substituting the differential operator with a difference operator, which converts the continuous differential equations to a set of linear equations governing the relationship between the discrete points in the problem domain. By solving this set of equations, the value at these discrete points can be computed. One example of this technique is the Finite Difference Time Domain (FDTD) Methods which recently attracts large interest. The method of moments uses the integral form of Maxwell's equations. The unknown equivalent sources or polarization currents are discretized, resulting in a smaller linear system with better computational efficiency. However, in most cases, knowledge of the Green's functions, either in simple functional form or series expansion form, are required. This limits the use of MoM to simple structures in which Green's functions are available. The integral equation needs enormous amount of analytical effort to derive and implement. Different formulation are needed when the structure changes. This makes the MoM unsuitable for solving general electromagnetic problems. Finite element methods solve Maxwell's equations by discretizing the electromagnetic field directly. This results in a much larger linear system that requires more computation time and storage space than does MoM. However, finite element methods have great flexibility in modeling complex structures. Since the differential form of Maxwell's equations is used directly, the formulation is independent of the structure, making FEM a good general electromagnetic problem solver. However, because the fields must be discretized in the whole problem domain, it is not suitable to be used in unbound region like antenna problems. Recent advances in absorbing boundary conditions (ABCs) [18], perfect matching layers [6] (PMLs) [54],

3 isotropic absorbers and boundary element methods (BEM) [6] have addressed this limitation by enclosing the problem domain with an artificial boundary. The first three methods use approximate boundary conditions which requires careful choice of the boundary material's constituent parameters in order to get correct results. The last method enforces rigorous boundary conditions by using the Green's function on the boundary. However, the resulting matrix is partly sparse and partly full, which causes difficulties in storing and solving the linear system. FDTD, like FEM, solves Maxwell's equations directly. Hence, it shares most of the benefits and limitations of FEM. One key difference between these two methods is that the solution space of FDTD is in the time domain. To derive the frequency domain solution, Discrete Fourier Transform (DFT) must be performed on the time domain solution. Although a range of frequency response curves can be produced in one simulation, additional error may be introduced by the DFT process. From the above discussion, we know that MoM is superior to the FEM and FDTD techniques in computational efficiency. However, it lacks the flexibility of the FEM and FDTD in modeling complex structure commonly found in MMIC's. The large size of MMIC's makes FEM and FDTD computationally inefficient using today's computer technology, especially when the MMIC integrates antenna elements. Thus, a new full-wave technique which can address the limitation of MoM in flexibility but retain its computational efficiency is necessary to simulate a MMIC faithfully in a reasonable amount of time. 1.2 Objectives The goal of this dissertation is to develop a full-wave technique which is both computationally efficient and flexible in modeling complex structures such that it can

4 facilitate the design of MMIC's. To reach this goal, three approaches are developed in this dissertation. The first is a hybrid technique which is based on MoM and FEM. In this technique, the formulation of MoM is retained to achieve computational efficiency in regions where Green's functions are available. The FEM is only used in the formulation of complex structures to limit its computational cost. The separation of the two formulation makes the choice of basis functions in MoM independent of FEM. This enables us to use the knowledge of the electromagnetic field to choose suitable basis functions such that the unknowns in MoM formulation can be reduced. The sparsity and symmetry of the FEM matrix is also retained due to the separate formulation. Since the difference between this hybrid technique and the MoM is only the replacement of the Green's functions with the solution of FEM, this technique can apply established MoM procedures to many classes of problems that can be formulated but are very difficult to solve numerically only by the MoM. The FEM can also be easily replaced with other field solution techniques. The second is a wavelet-based MoM. In this method, wavelets are used as the basis functions in the MoM formulation. The resulting impedance matrix is thresholded to reduce storage size. Due to the use of wavelets, huge saving in storage size can be achieved with only small penalty in accuracy. The third is parallel computing. Two different parallelization schemes are implemented to speed up the process of matrix-fill and fast wavelet transform (FWT). The above approaches will be tested in practical examples to show the aforementioned merits.

5 1.3 Overview In this thesis, a hybrid full-wave technique that combines MoM and FEM is developed. This technique features the use of FEM in solving the electromagnetic field in geometrically complex regions and MoM to solve Maxwell's equations in regions where the Green's functions are available. Wavelet-based bases are also used in the formulation of MoM, which can reduce storage requirement. To further improve computation speed, parallel computation models are developed for matrixfill process, FWT. In Chapter 2, rigorous analysis of the hybrid technique is presented. The formulation primarily follows the derivation of MoM technique except in regions where Green's functions are substituted by the equivalent solution of FEM. The problem domain is separated to two parts. One consists of regions in which the Green's functions are available. Another part consists of complex structures in which the Green's function is not known. By using equivalent principle, equivalent sources are added on the boundary between the two parts. Galerkin's MoM is used in the formulation of the first part, and FEM is used in the second part. By matching the boundary conditions, the solution of FEM can be incorporated to the MoM formulation. Thus, the problem is solved. In Chapter 3, cavity-backed patch antennas are analyzed using the hybrid technique. Two feeding schemes are used: one is a slot-coupled microstripline feed and the other a coaxial line feed. Spectral Domain Green's functions of a conductor backed dielectric slab were used in the MoM formulation. In the FEM formulation, we used vector-based tetrahedral elements which have the benefit of producing spurious-free solutions. The input impedance and field distribution of the patch antennas were

6 computed. Good agreement with the simulated result of FDTD was achieved. In Chapter 4, the hybrid technique is applied to microstripline-fed cavities and couplers. A micromachined resonator which is compatible with planar circuits was designed and fabricated according to simulated data. The measured performance of the resonator showed the potential of achieving a Q as high as a conventional metallic resonator. A design of a very wide band-pass microstripline vertical transition is also shown. In Chapter 5, a Wavelet-based Method of Moments is introduced. The multiresolution analysis is briefly discussed. In this dissertation, we concentrate on the use of two dimensional wavelet bases on patch antenna problems. Theses two-dimensional bases constructed by combining Haar and linear B-spline wavelets are suitable for discretizing the current on patch antennas. Due to the finite support and zero average property of the wavelet bases, the resulting impedance matrix of the patch antenna can be thresholded to reduce memory size. This wavelet based MoM is applied to patch antenna arrays. Effect of thresholding on the antenna characteristics are investigated. Simulations of 3 x 3 and 6 x 6 patch antenna arrays are performed to show the reduced memory requirements attained by thresholding the impedance matrix. Chapter 6 deals with the utilization of parallel computing to reduce computation time. Two parallelization schemes are developed for matrix-fill process and FWT using portable Message Passing Interface. Due to the non-uniform computation time of the matrix-fill process, a dynamic job allocation scheme is applied to optimize the use of CPUs. In this scheme, one node is reserved as a controller to assign jobs to available free nodes and gather completed results. For FWT, the computation time of each job is fixed and smaller comparing to matrix-fill process such that the

7 communication time between the nodes is not negligible. A static job allocation scheme is applied to reduce the communication overhead. Finally in Chapter 7, summary of the achievements in this dissertation and future research direction are given.

CHAPTER II Hybrid MoM/FEM Techniques 2.1 Introduction In this article, we present a hybrid technique, proposed by Yuan et al. [10], that combines the advantages of the MoM and the FEM (MoM/FEM), and has been successful in solving various 2-D and 3-D scattering problems [10]-[16]. This technique uses FEM to formulate the electromagnetic field in geometrically complex regions and uses the integral equation technique to solve Maxwell's equations in regions where the Green's function is known and where FEM fails to accurately satisfy the pertinent radiation conditions. This technique first divides the problem domain into two regions: one in which the Green's functions are available and one in which the Green's functions are not available. These regions will be analyzed by MoM and FEM respectively. Through the boundary conditions between the two regions, the formulations of MoM and FEM are combined to give the final solution. In the following sections, the hybrid technique is applied to a cavity scattering problem. Detailed formulations are supplied to give the reader insight into this technique. 8

9 Ei,%.Hi Aperture Sb PEC cavity PEC obstacle g A general scattering problem which consists of a cavity, an obstacle and a current source. 2.2 Method of Moments Formulation Fig. 2.1 shows a cavity scattering problem which consists of a cavity with an aperture 5a on its walls, an obstacle with surface Sb, an impressed current Ji and the electric and magnetic fields Es, Hi produced by the impressed current. Assume that the cavity and the obstacle consist of perfect electric conductors (PEC) and the source is an electric current. The following development of formulation can also be applied to the cases with multiple cavities, obstacles and apertures. By using the equivalence principle [19], the original problem is changed to an equivalent one as shown in Fig 2.2. The aperture on the cavity is replaced by a PEC with equivalent magnetic current M". The obstacle is replaced by a fictitious equivalent electric surface current Jb. Ma and Jb must satisfy the conditions, Ma -n x E (2.1) b = x H, (2.2) where n is the outward normal vector on the surfaces and E and H are the unknown

10........S.,.... ",-A Jrb=nxH vnJ= n x H Figure 2.2: The equivalent problem of Fig. 2.1. total fields on S, and Sb, respectively. In this manner, the volume of interest is separated into two parts. One is the unbounded region outside the cavity including the sources J, M' and Jb; the other is the region inside the cavity with source -a flowing on the surface initially occupied by the aperture. The total fields can be represented as E1 - E (J) + Ei(Jb) + E1(M) (2.3) H1 = H(ji)+ (Jb) + H(Ma) (2.4) outside the cavity, and E2 = E2(-M) (2.5) H2 = H2(-M) (2.6) inside the cavity, where subscript 1 denotes the region outside the cavity and subscript 2 denotes the region inside the cavity. If the forms of E~ and Hi are readily available, the two terms E1(JY) and H1(JY) in (2.3) and (2.4) can also be represented by E~ and H1, resulting in slightly different expressions.

11 Since the Green's function of the exterior region is available, E1 and H1 in (2.3) and (2.4) can be derived from integral equations. Note that the Green's function includes the effect of the completely closed PEC cavity. By matching the boundary conditions on S, and Sb, two equations linking the two regions are derived: nx H = n xH2 on Sa (2.7) n x E1 = 0 on Sb. (2.8) Let Ma and Jb be expanded by basis functions Mn and J} as follows: Na Ma mMaa (2.9) n=l Nb J = Jn (2.10) n=1 where Na, Nb are the numbers of the basis functions and mn, jb are the unknown coefficients. By substituting equations (2.3)-(2.6) and (2.9)-(2.10) into (2.7) and (2.8), and applying Galerkin's moment method, the following coupled matrix equations are obtained: -[V] - [Zb][Jb] + [Ta][Ma] = 0 (2.11) [Ci] + [Cb][Jb]] -[ya][a] = [yc][Ma]. (2.12) where [jb] and [Ma] are column vectors with the n-th element equal to jb and in, respectively. Using the notation (A,B)s= J A BdS, (2.13) S the matrix and vector elements in (2.11) and (2.12) are defined by: Ymn = (-M, Hi,(Mn))s~ Na x Na matrix (2.14) Ymn (-MM,H2(Ma))sa Na x Na matrix (2.15)

12 Zmn = (-JmEl (J,))Sb Nb x Nb matrix (2.16) tn = (Jm, El(Mn))sb Nb x Na matrix (2.17) C -(M, H(J))S Na x Nb matrix (2.18) m= (MMH (Ji))sa Na x 1 column vector (2.19) i b - -b v = (-JmEl(J ))Sb Nb x 1 column vector. (2.20) Note that from reciprocity, the elements of [Ta] and [Cb] are related by cb = _tam Cron nm' Solving equations (2.11) and (2.12), we have [Ma] = ([C'][Z] [Ta] - [Ya]- [Yc])1 ([Cb][Zb]-l[i]- [Ci]) (2.21) [jb] = [zb]-1 ([Ta][Ma] - [Vi]). (2.22) Equation (2.21) and (2.22) indicate that the unknown column vectors [M"] and [Jb] can be calculated if all the matrix elements in (2.14)-(2.20) are known. Since in the most general case the Green's function inside the cavity is not available, matrix [YC] is still not known at this point. In the next subsection, we will show how to calculate [YC] by using FEM. 2.3 FEM Formulation Consider the time harmonic Maxwell equations, V x E = -jwB-M (2.23) V x H = jwD +J (2.24) V-D = p, V-B = pm with the constituent relations D = E

13 E = pH where E: electric field intensity H: magnetic field intensity D: electric flux intensity B magnetic flux intensity J: impressed electric current M impressed magnetic current pe: electric charge density Pm ' magnetic charge density e: permittivity p: permeability w: angular frequency Note that ejwt time dependency is assumed. After substituting (2.24) and the two constituent relations to (2.23), the magnetic field H2 in the cavity of previous section satisfies V x (-j. V x H2) = jLwcH2 + M, (2.25) Jwec where pc is the permeability of the medium inside the cavity and sc the permittivity. Let q denote an arbitrary vector test function. The inner product of q with both sides of equation (2.25) in the whole volume of the cavity gives //[V x (- -V x H2)] * dv - (j[jcH2 j) dv = a M * ds (2.26)

14 From vector identity V - (A x B)B -V x A = A- V x B, the integrand in the first term of (2.26) becomes 1 1 — V x H) V x -V x V- V (- -V x H2))< x. (2.27) JUWc J^Wc Substitute (2.27) to (2.26) and apply Stokes Theorem to the second term of (2.27), we have ( — V x H2) V x Odv - (jwcH2 = M d (2.28) Note that the second term of (2.27) vanishes after the volume integration because 1 V x H = E = 0 on the boundary. jwe Expand the magnetic source Ma as in (2.9) and expand H2 into No basis functions: N, H22 = ji (2.29) i=l where qij are the unknown coefficients. By applying Galerkin's method to (2.28) for each Ma, we get the following matrix equation: [A][D]= [B], (2.30) where aij- = (-. V x i V x )j) dv - JJJll( ci. j)dv N, x N matrix, (2.31) bij= - MJ 4i ds N, x Na matrix, (2.32) bij -(fij No x Na matrix. (2.33) Due to the large size of [A] commonly seen in 3-D problems, the direct inversion of [A] is not always possible. Instead, iteration methods like the bi-conjugate gradient method [6], which is also applied in this dissertation, have to be used. A 3D

15 edge-based tetrahedral vector element is used through this dissertation. A detailed definition of this element is included in appendix A. 2.4 Combining the solutions of MoM and FEM After solving [4(] from (2.30), we may find [YC] by substituting (2.29) into (2.15) as follows: Ymn = (-Mm, H(Mn))s, No = in(-Mm i)sa (2.34) i=l If we define the matrix elements of [Yk] as Ymn = (- M,,~,n)Sa, (2.35) equations (2.34) - (2.35) result into the following matrix equation: [Y] = [YO][q]. (2.36) By comparing (2.32) and (2.35), we find that [Ye] is related to [B] by [YO] = -[B]. Now that all the matrices in (2.14)-(2.20) have been derived, the unknown [Ma] and [Jb] can be calculated from (2.21) and (2.22). In the following two chapters, physical structures will be analyzed by th ie formulation introduced in this chapter.

CHAPTER III Cavity-backed Patch Antennas 3.1 Introduction One serious limitation of patch antennas is its narrow bandwidth. One way to circumvent this limitation is to increase the thickness of the substrate. However, thick substrates cause the propagation of surface waves which reduces the radiation efficiency and increases the coupling between the array elements. To avoid these undesired effects, a metallic cavity is used to enclose the patch antenna in order to suppress the surface waves. A cavity-enclosed aperture-coupled circular-patch antenna operating at Ka band was built in [1]. It showed a 2:1 input VSWR bandwidth of 12%. Another example is a cavity-backed aperture-coupled rectangular-patch antenna built in [2] with a 19% bandwidth from 8.5 to 10.3 GHz. In addition to substrate mode elimination, the metallic cavity can also serve as a heat sink to improve heat dissipation. Examples of cavity-backed patch antennas with different feeding structures are given in Fig. 3.1 and Fig. 3.2 Two techniques are found in the theoretical analysis of these structures. One is the MoM; the other is the finite element/boundary integral method (FE/BIM). The MoM has been used in the analysis of cavity-backed circular-patch antennas by Aberle and Zavosh [3, 4], and in the analysis of cavity-backed rectangular patch 16

17 Slot on ground plane antennas by Lee et al. [5]. Although the MoM is computationally efficient in solving spatially unbounded problems, it can only be applied to a few cavity configurations in which the Green's function is known. On the other hand, FE/BIM- [6]-[9] can be applied to various cavity configurations, but due to the coupling between the FEM formulation and integral equations, the resulting matrix is always partially sparse and non-symmetrical. This increases the computational cost in terms of storage and computer time. In the following sections, we will use the hybrid technique described in Chapter II to analyze the cavity-backed patch antennas shown in Fig. 3.1 and Fig. 3.2. The

18 Figure 3.2: The structure of a coaxial line-fed cavity-backed patch antenna. reflection coefficient and the equivalent input impedance are shown and compared to the results obtained using the FDTD technique. The convergence of the numerical results is discussed. Computation time and memory requirements of a given structure are also included. 3.2 Spectral Domain Green's Functions The evaluation of the impedance matrix of patch antenna problems involves the integration of the Green's functions of a conductor-backed dielectric slab with the testing and weighting functions. In Cartesian coordinates, assume that the dielectric slab is parallel to x-y plane and occupy the space between z = 0 to z = d. Also assume that the conductor is located on x-y plane. The interaction of two y-directed

19 surface electric currents, say Ji and Jj for example, located at z = d plane takes the following form: zij = if if Ji(x, y) Jj(x', y') GE (x, y, d; x', y', d) dx dy dx' dy', (3.1) where G"J is the Green's function of a conductor-backed dielectric slab with thickness d. The superscript of the Green's function denotes the orientations of the field and the source. The subscript denotes the types of the field and source. In spectral domain, the Green's function is represented as 00 GJY(x, y, d; x', y', d) = f f Gj * ejk(- x')eiky(Y-Y') dkx dku. (3.2) -00 The exact form of the spectral Green's functions used in this dissertation are listed in Appendix B. Substituting (3.2) to (3.1), we have 00 zij = |f if Ji(x, y) J,(x', y') GEJ ex(xx)eky(y') dkx dk dx dy dx' dy. (3.3) -00 The above equation involves a sixfold integration which is difficult to evaluate numerically. Instead, if the testing and weighting currents are transformed to spectral domain, the spatial integration in (3.3) can be removed and reduced to a twofold integration in spectral domain. That is, 00 zij = J ky) J (k, ky) Gj'J(kx, ky) dkx dky, (3.4) -00 where Jij k ky) = f Jij (x, y)e(kxx+kY) dx dy, (3.5) and the "*" is the conjugate operator. To further facilitate numerical integration, the double infinite integration in (3.4) is converted to polar coordinate and becomes one finite and one infinite integrations as follows 27r 00 zij = J* Jj M - d/3 dO. (3.6) o o

20 Due to the symmetry of the structure in x and y directions, the spectral domain Green's functions are either even or odd to kx or ky. If Ji and Jj are real and separable in x and y, Ji and Jj can be represented as Ji(kx, Ik) =ky) J=(kx)J (ky) Jj(kx), ky) Jj (kx)Jj (ky), and have the following properties: s(J"ux*(kx,y) j)'Y (kx,y)) is even to ky,.(Jx'Y*(kxy)JixY(kxy)) is odd to kxy, where? and S mean the real part and imaginary part of its arguments respectively. Applying these even and odd relationship to (3.6), the integration range can be further reduced to the first quadrant. That is, 7r/2 oo zij 4 JJ Px(Jx*J). pY(JiY*Jy) GJ *d3 dO, (3. 7) 0 0 where R? if GEJ is even to kcy, -ZY = (3.8) ~ if GEJ is odd to kxy. The evaluation of (3.7) still poses a challenge due to the surface wave poles in the Green's function. When the integration path of (3.7) passes these surface wave poles, the integrand becomes infinite, causing the numerical integration to diverge. A pole extraction technique is used to circumvent this problem [25]. Assume that 3o is one of the surface wave pole. Let I denote the integrand of (3.7). I can be reformulated as = (/3 p) where g(/3) = 0. (3.9) s(1P'

21 To remove the singularity at /3o, I is subtracted by a function which approximates the behavior of I near /o. Let this function be Is. We choose s = f(Po, 0) 2/o2 g'(Po) P2 - o2' To show that Is indeed approximates I near PJo, we take the limit of I at /o as follows lim I lim fi O)/ #-Po p0-po (3 - /o)g'(0) l f(/fo, 0) /o P3-+o g'(0o) -/3o - lm f(/3o,) 2/o2 p-to g'(3 o) (/3- o)(o + Po) lim Is. Ago Reformulating (3.7), we have Zij = 4 1 (I-IS)d/+ I'sd/} dO (3.10) 0 0 0 ir/2 oo r/2 f = 4 / (I- IS) d d O3 + 4 -ji/o' 'o) dO. (3.11) 0 0 0 Note that the second term of (3.10) is integrated analytically with respect to /. (3.11) is numerically integrated by 16-point Adaptive Gaussian Quadrature. The infinity range is truncated at 150ko as in [24, 27] where ko is the free space propagation constant. 3.3 FEM Mesh Generation There are a few software packages, for instance, I-DEAS, which can discretize an arbitrary 3-dimensional structure by tetrahedrons. However, the process involves interactively inputing the geometry of the structure, setting a cue for mesh generation, and outputting data which needs to be translated for further processing in

22 Figure 3.3:,. x g A hexahedron with parallel and rectangular top and bottom planes can always be decomposed to distorted bricks as show in Fig. 3.4. Figure 3.4: g A distorted brick. Note that only the vertical edges are distorted. The top and bottom planes are still parallel and rectangular. FEM codes. These processes are always time-consuming, cumbersome and requiring huge amount of computational resources. Every time the parameters of the structure are changed, the interactive process has to go over again, making it difficult for a researcher to perform a comprehensive study on a given structure. Writing a general tetrahedral mesh generator is not an easy task. However, if only

23 hexahedral structures, which are common in planar circuits, are considered, mesh generation becomes very simple. In this dissertation, a tetrahedral mesh generator for a hexahedron is implemented using distorted non-uniform bricks. As shown in Fig. 3.3, a hexahedron can be subdivided to bricks. The only limitation here is that the top and bottom planes of the hexahedron have to be parallel and rectangular. Each brick can be further divided to 5 tetrahedrons as shown in Fig. 3.4. There are only two spatially unique ways of decomposing a brick to 5 tetrahedrons. Each of them is only a 90 degree rotation to the other. Note that, these two ways have to be used alternately such that the common edges of the tetrahedrons between two adjacent bricks remain the same. The input file format is as follows m XQ X1 * Xm m+l terms n YO Yl... Yn n+l terms k Z0 Z1... Zk k+l terms m 0 X 1 ~.. Xm -, m'+l terms n' yo Y... Yn n'+i terms xdl xd2... xdm m terms ydl yd2... ydn n terms zdl zd2... zdk k terms xd[ xd2... xd', m' terms yd[ y... yd *** ni terms where m: number of divisions in x direction on the bottom plane, Xi: node position in x direction on the bottom plane, xdi: number of bricks between node xi-1 and xi on the bottom plane n: number of divisions in y direction on the bottom plane, yi: node position in y direction on the bottom plane,

24 ydi: number of bricks between node yi-1 and yi on the bottom plane k: number of divisions in z direction, zi: node position in z direction, zdi: number of bricks between node zi-1 and zi. Similar definition applies to m', n', x', xd', y\ and yd' except that they are on the top plane of the hexahedron. Let the total number of nodes along x, y and z-direction be NX, Ny and Nz respectively. We have m ml Nx = Exdi+ 1 = xd+l, i=l i=l n n' Ny = Eydi+l = Eydi+l, i=1 t=1 k Nz = Ezd + 1 i=1 Let nijk = (aik, bk, Ck) be the node coordinate of node (i,j, k), where 0 < i < N,, 0 < j < Ny, and 0< k < Nz. (aik, bjk, Ck) can be found as follows aiNz - aio aik = aio+ (ck-co) * a — a CNz - C bjk = bjo + (ck -c ) bjN - bj. CN, - Co Note that aiNz aio, bjN,, bjo, and Ck can be found directly from the input file. After knowing the coordinate of each node, the tetrahedral elements and the edge bases are easily produced from each bricks as shown in Fig. 3.4. Different from electric field formulation, the edges on the surface of PEC are not removed. If there is a planar PEC in the structure, edges have to be created on both sides. As shown in Fig. 3.5, Element 1 and Element 2 are on different sides of the PEC. Element 1 consists of edges 1 - 6, and Element 2 consists of edges 10 - 12. Although the position of edges 4 - 6 and edges 10 - 12 are the same, they are treated as different unknowns.

25 nl lement l(el, e2, e3, e4, e5, e6) el: nl-n2 n2 PC e2: nl-n3 e3: nl-n4 n3.. - n4 / e4: n2-n3 / / e5: n3-n4 / e6:n4-n2 / // n5 Element 2(e7, e8, e9, elO, el 1, e12) e7: n5-n2 e8: n5-n3 e9: n5-n4 nlO: n2-n3 nll:n3=n4 n12: n4-n2 g Edges on a PEC plane have to be created on both sides. In this figure, Element 1 consists of el, e2,..., e6, and Element 2 consists of e7, e8,.., el2. Note that the position of e4 to e6 and elO to el2 are the same. The capability of assigning different node positions on the top and bottom planes gives the flexibility in defining geometry. For example, in the cavity in Fig. 3.1, the sizes of the patch and the slot can be changed independently to each other. 3.4 Slot-Coupled Microstripline-fed Patch Antennas Patch antennas are known to be lightweight, conformable, economical to fabricate, and easily integrated with microwave circuits. At low frequencies, patch antennas are effectively fed by direct contact of the patch with a coaxial line or a microstrip line. Fig. 3.2 shows an example of coaxial line-fed patch antennas with the center conductor of the coaxial line terminating on the patch. At millimeterwave frequencies, however, the size of the patch and the size of the coaxial line or

26 microstrip line become comparable, which affects the electrical performance of the patch antenna. The connection of the coaxial line to the patch also becomes difficult due to the drilling and soldering procedure required. In 1985, Pozar [20] proposed an aperture-coupled microstrip line-fed patch antenna in which the feedline does not directly contact the patch. This structure consists of two substrates separated by a ground plane. The patch and the feedline are situated on two different substrates, and coupling is achieved through an aperture on the ground plane. This structure is suitable for millimeter wave monolithic arrays since the ground plane provides good isolation between the patch and the feeding structure, preventing interference from the active devices. It also makes it possible to choose different substrates in order to optimize the performance of the patch antenna and the feedline separately. Fig. 3.1 shows a microstrip line-fed, cavity-backed patch antenna based on this aperturecoupling technique. Now we will model the patch antenna structures described above using the hybrid technique presented in Chapter II. Results will be validated through comparison to FDTD data1. The dimensions of the cavity and the patch are kept the same for both structures. Design examples in which the patch antenna is matched to the feedline are given to show the capability of this technique as a practical design tool. Calculation of the reflection coefficient, input impedance and field inside the cavity are also shown. This technique is also capable of calculating other important design data such as antenna patterns, scattering patterns, backlobes,... etc., though these are not shown here.

27 Patch Air Microstrip line.,u lengt: s; Stub length Ls Cavity CPatch Lap Microstrip line -. ---Microstr - -- ix Ls / / II Slot y microstrip feedline. The patch is fed by a microstrip line through a slot opened on the bottom wall of the cavity. Comparing to the general problem in Chapter II, the microstrip line can be considered as the obstacle b and the exterior region includes the dielectric layers and half spaces above and beneath the cavity. The aperture 1All FDTD data presented in this dissertation were contributed by Nihad I. Dib.

28 Sa includes the opening around the patch and the slot on the bottom of the cavity. Since the feeding structure is the same as that of [21], we use the same technique to model the currents on the microstrip line as in [21]. The origin is set at the center of the slot. The Ji described in Chapter II is modeled as an incident semi-infinite traveling wave mode with unit amplitude as follow: Ji(x, y) = Jc(x, y)y - jJs(x, y)y, where Jc(x,y) = w-cos k(y-Ls), - < x <, -(y + ) < y L -, Js(x,Y) = w sinke(y-Ls), -W < x <, -yoo < y < Ls, and, Wf is the width of the microstripline, ke is the phase constant of the propagation mode on the microstripline, LS is the length of the feedline section from the center of the slot to the open end, and yoo is an infinite number. In numerical integration, the semi-infinite traveling wave mode must be truncated at a finite length to avoid singularity. We found that 3 wavelengths are sufficient. That is, Ls + yoo = 3 27r/ke. Note that the range of Jc is shifted by a quarter wavelength to avoid a nonzero current at the open end of the microstripline. The Jb in Chapter II is modeled as a reflected semi-infinite traveling wave mode with amplitude F plus several piecewise sinusoidal (PWS) modes near the end of the microstrip line to account for the effect of the disturbance caused by the open end and the slot. That is, Nb Jb(x, y) = InJ(x, y)A + F(Jc(x, y) + jj,(x, y)), n=l

29 PWS modes V \ Mic /ostrir lin / Jc t —I --- Traveling wave mode Slot Figure 3.7: The basis functions used on the microstripline. where J(x y) 1 sink(hly-) W1 snk (h wn' W- sin, — chb 2 - 2, Yn -h < Y-< Yn + h hb = (L- yo)/(Nb+ 1), y,: center of the PWS modes, Nb:number of the PWS modes. as shown in Fig. 3.7. Since the gap between the patch and the cavity is small, we assume that the electric fields only exist perpendicular to the edge of the patch. Similarly, we assume the width of the slot is small, such that only the y component of the electric field exist in the slot. Thus, the M' in Chapter II can be represented as M a MPX(x, y) + MPY(x, y) + M (x, y), MPX, Mp: magnetic currents in the gap between the patch and the cavity, MS: magnetic currents in the slot. As shown in Fig. 3.8, each of MPX, Mpy, and Ms are decomposed to PWS modes as on the feedline, except for the corner of the patch where asymmetric PWS modes

30 Roof-top basis functions Tetrahedral elements Slot Ms This figure shows how the basis functions used in MoM match the tetrahedral elements used in FEM. To simplify, not all of the PWS modes are shown. are used. The boundaries of the tetrahedral elements are set to match that of each basis function used in MoM. In the dielectric layers and free space, spectral domain Green's functions are used, and the evaluation of the pertinent matrix elements in (2.14)-(2.18) follows the same procedure as in Section 3.2. The reference plane of the the normalized input impedance is set below the center of the slot as in [21] and can be calculated by 1 - Fe-j2keL in = + re-2keL, (3.12) 3.4.2 Numerical Results An antenna that is matched to the feedline is designed with parameters as follows: patch size 27.78 mm x 27.78 mm, cavity size 32.52 mm x 32.52 mm x 3 mm, dielectric constant inside the cavity 1.0, thickness of the dielectric cover 0.508 mm and dielectric constant 2.2, slot length 13.34 mm, slot width 1 mm, microstrip line

31 0.9 0.8 0.7 a FDTD ---- 0.6 0, 0.4 0.3 0.2 0.1 3.6 3.8 4 4.2 4.4 4.6 4.8 5 Frequency(GHz) Figure 3.9: 3.9 Comparison of the reflection coefficient of a microstrip line-fed cavitybacked patch antenna calculated by the method presented in this dissertation and the FDTD method. The FDTD data were contributed by Nihad I. Dib[40]. width 1.528 mm, stub length 6.5 mm, substrate thickness 0.51 mm, and dielectric constant 2.33. The reflection coefficient of this structure is shown in Fig. 3.9 and compared to the results calculated by the FDTD method. The figure shows that the resonant frequency calculated by MoM/FEM is near 3.99 GHz with a 2:1 VSWR bandwidth about 5%, while the resonant frequency calculated by the FDTD method is 4.01 GHz. This corresponds to a difference of 0.5%. Also the minimum reflection coefficient is near zero for the method of this paper, and 0.025 for FDTD method. Overall, the two curves match very well. The input impedance calculated by both techniques is shown in Fig. 3.10. The effect of the slot length on the reflection coefficient is shown in Fig. 3.11 on a Smith chart. It can be seen that the impedance locus forms approximately a circle and indicates that the coupling factor and the resonant resistance increase as

32 1.2 I Resistance/This paper - 8 / \\ Reactance/This paper ----- 0.8/ \- Resistance/FDTD ---- Reactance/FDTD -- 0.6 - 3.6 3.8 4 4.2 4.4 4.6 4.8 5 data were contributed by Nihad I. Dib[40]. -0.2 -0.4 -0.6 3.6 3.8 4 4.2 4.4 4.6 4.8 5 Frequency(GHz) Figure 3.10: The normalized input impedance of a microstrip line-fed cavity-backed patch antenna. The parameters are the same as in Fig. 3.9. The FDTD data were contributed by Nihad I. Dib[40]. antennas [21]. slot is y-directional. 3.4.3 Numerical Considerations tenna is investigated. Fig. 3.13 and Fig. 3.14 shows the convergence in terms of the

33 5 GHz 3.5 GHz Cavity _ Patch Slot — ~ Lap = 12mm - - -I- - - Lap=13.34mm Lap --- Lap=16mm Figure 3.11: Comparison of the normalized input impedance of microstrip line-fed cavity-backed patch antennas with different slot lengths. Other parameters are the same as in Fig. 3.9. number of unknowns on one side of the patch and on the slot. The results indicates that 5 basis functions on one side of the patch and 3 basis functions on the slot are sufficient. Fig. 3.15 shows the results obtained by the hybrid technique for different discretization schemes in the cavity. In the legend of the figure, the first number in the parentheses denotes the total number of basis functions around the patch and the second denotes the total number of edges in the cavity. The curve in this figure is the same as that shown in Fig. 3.9, serving as a reference. Data Set 2 shows the result of increasing the basis functions from 20 (5 on each edge) to 28 (7 on each edge). Data

34 Total E-field (A) x (B) Y >x y I >z (C) 0 0.2 0.4 0.6 0.8 1 a a b Figure 3.12: Total E-field intensity in the cavity of Fig. 3.1 at resonant frequency (3.995 GHz). The parameters are the same as Fig. 3.9. (A) and (B) are the cross section views along lines aa and bb shown in the inset, respectively. (C) is the view just beneath the patch. The black and white lines indicate the boundaries of the patch and slot. This figure is drawn according to the scale of the physical dimensions.

35 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 Comparison of S11 with Nb 0+.- - - \ +/ e +/ ~ Nb=3 v /,Nb=5 --- ' 8 Nb=7 + A / \ +, \ / 'I / /.+ O I. t i I.. I 3.6 3.8 4 4.2 4.4 4.6 Frequency (GHz) 4.8 5 Figure 3.13: The convergence check of the same case as Fig. 3.9 using the hybrid technique. Nb denotes the number of basis functions on one side of the patch. Comparison of S11 with Ns en mCO 1 0.8 0.6 0.4 0.2 0 3.6 3.8 4 4.2 4.4 4.6 4.8 Frequency (GHz) 5 Figure 3.14: Fg The convergence check of the same case as Fig. 3.9 using the hybrid technique. Ns denotes the number of basis functions on the slot.

36 0.8 a) \ / 1 (20, 8640) 0 0.4O - V ~ 2(28,11760) a) ) \ / + 3 (20,17280) 0.2 - \o 4 (20, 34560) 3.6 3.8 4 4.2 4.4 4.6 4.8 5 Frequency(GHz) g The convergence check of the same case as Fig. 3.9 using the hybrid technique. The first number in the parentheses denotes the number of basis functions around the patch and the second denotes the total number of edges in the cavity.

37 Set 3 shows the result of doubling the number of discretizations in z direction of the cavity and Data Set 4 corresponds to doubling the number of discretization in the cavity in both x and y directions (refer to Fig. 3.12 for the definition of directions). Other basis functions used are: 3 on the slot and 14 on the microstrip line. The number of nodes for each data set is also shown in Table 3.1. From the figure, it can be seen that the convergence of the numerical result is quite good except near the resonant frequency. The total CPU time for one data point in the curve is 3.3 hours on a 50 MHz Sun Sparc 20 workstation with 32 MByte main memory of which 26% is used. It is worth mentioning that 73% of the CPU time is used for filling MoM matrices and 26% is used in solving the FEM matrix by the bi-conjugate algorithm. Note that the CPU time spent on MoM is larger than that spent on FEM in this case. However, the latter is very sensitive to the uniformity of the mesh in the cavity, and for more complex cavities such as the coaxial line-fed one, the CPU time spent on FEM will be much higher. For the same case, the CPU time used by FDTD is 7.5 hours on a CRAY YMP8/864 supercomputer. The number of FDTD cells is 188 x 166 x 288 along x, y and z directions, respectively, with Ax = 0.25 mm, Ay = 0.25 mm and Az = 0.125 mm. The total number of time steps is 8000. The factor influencing convergence most is the distance between the top absorber and the dielectric cover. We found that a distance of 30 mm is necessary to get a convergent result. It seems that FDTD uses more CPU time, but we have to point out that FDTD gets the entire frequency range in one computation, while the hybrid technique must be repeated for each frequency point.

38 Conductor i-Fictitious surface S Coaxial line Coaxial line Figure 3.16: The side and top views of the patch antenna shown in Fig. 3.2.

39 Table 3.1: Number of nodes of each data set of Fig. 3.15 in each direction. Data Set x y z 1 25 25 4 2 29 29 4 3 25 25 7 4 49 49 4 3.5 Coaxial Line-fed Patch Antennas 3.5.1 Formulation Fig. 3.16 shows the side and top views of a coaxial line-fed, cavity-backed patch antenna. The patch and the cavity are the same as in the previous section. In the coaxial line, the cavity region is extended to a fictitious surface at which all higher order coaxial line modes generated by the junction to the cavity have attenuated to practically zero. Compared to the general problem in section Chapter II, this structure does not have the obstacle b. The cavity region is extended to include the section of the coaxial line above the fictitious surface. The aperture consists of two parts; one is the opening around the patch, the other is the fictitious surface S at - = zo in the coaxial line. The exterior region includes a semi-infinite section of coaxial line, the dielectric layer and the upper half space. By choosing a suitable operational frequency range, we can assume that only the fundamental mode propagates in the coaxial line. In this manner, the field pattern of the fundamental mode is used directly, instead of the Green's function when evaluating the fields inside the semiinfinite coaxial line. Due to the position of the fictitious surface, we can use the mode pattern of the fundamental coaxial line mode as the only basis function. Comparing

40 to Chapter II, we have -Ei= op (p) e Ik(ZZO)t Hi -= ++(P) e-jk(z-zo), 77 M- MPX(x, y)x + MP"(x, y)y + M8?(p)j, k - /ai, where, and e are the permittivity and permeability of the coaxial line, and MPX and MP are the same as in the previous section. Let the reflection coefficient of the fundamental mode at the fictitious surface to be F. F is related to M* as follows 1 + F = Ms. (3.13) The reflection coefficient F' at the junction between the coaxial line and the cavity is related to F as: Fr = Fe2jk (3.14) where k is the propagation constant of the dominant coaxial line mode and f the distance from the junction to surface S. The normalized equivalent input impedance at the junction can be formulated as zin = ' (3.15) Other formulations are the same as in the previous section. Since the cavity including the coaxial line section is no longer a hexahedron, the software package I-DEAS is used to generate the mesh.

41 3.5.2 Numerical Results An antenna with the same cavity and patch as in the previous section was designed. The parameters of the coaxial line are: inner radius 1.037 mm, outer radius 3.620 mm, and dielectric constant 2.25. In order to obtain a good match, the feeding point is shifted from the center of the patch by an amount of 5.7 mm. The reflection coefficient of this structure is shown in Fig. 3.17 with comparison to the results calculated by the FDTD method. The figure shows that the resonant frequency calculated by the MoM/FEM technique is near 4.19 GHz with a 2:1 VSWR bandwidth about 5%, while the resonant frequency calculated by FDTD method is 4.22 GHz, which represents a difference of 0.72%. Also the minimum reflection coefficient is near zero for MoM/FEM, and 0.078 for FDTD method. Overall, the two curves match quite well excluding the fact that one curve seems to be shifted from the other by an amount of 0.03 GHz. This difference may be due to the approximation of the circular coaxial line by a rectangular coaxial line in the FDTD method. The input impedance calculated by the two methods is also shown in Fig. 3.18. The effect of the position of the feeding point is shown in Fig. 3.19 by plotting the input impedance loci as a function of frequency in a Smith chart. The position of the feeding point is indicated in the legend and inset in the figure. The frequency range is from 3.5 GHz to 5 GHz and each mark on the curves indicates a frequency increment of 0.1 GHz. From the figure, we can see that the impedance locus forms approximately a circle. The radius of the circle increases as the feeding point moves away from the center of the patch. This indicates an increase in the coupling factor, if we define the radius of the impedance locus circle as the coupling factor between the feedline and the cavity-patch structure as in [21]. The electric field intensity inside the cavity is shown in Fig. 3.20. Plots (A) and

42 0.9 0.8 0.7 *8 0.6 This paper -- FDTD data were contributed by --- — 0.5 0.4 0.3 0.2 0.1 3.6 3.8 4 4.2 4.4 4.6 4.8 5 Frequency(GHz) Figure 3.17: Comparison of the reflection coefficient of a coaxial line-fed cavitybacked patch antenna calculated by the method presented in this dissertation and the FDTD method. The FDTD data were contributed by Nihad I. Dib[40]. (B) show the field intensities on cross sections along lines aa and bb in the inset of Fig. 3.20, respectively. Plot (C) shows the field intensity at the top of the cavity. The black lines and white circles mark the boundaries of the patch and the coaxial line. This figure is drawn according to the scale of the physical dimensions. The field intensity is interpolated between sample points and smoothed to reduce the effect of finite discretization. The figure shows that the field spreads from the feeding point to the left and right edges of the patch and reaches maximum intensity at the two edges. Comparing Fig. 3.12 with Fig. 3.20, the total field underneath the coaxial line-fed patch is greatly distorted by the conductor post protruding from the coaxial line to the patch. This is because in both cases the dominant E-field beneath the patch is z-directional and is forced to diminish near the post in order to satisfy the boundary

43 2 I 1.8 1.6 / \ 1.4 / y \ Resistance/This paper 8 /,,-: \. Reactance/This paper -- ^ //// ",\ \\ 'Resistance/FTD.. / 1.2 // \ Reactance/FDTD -- 0.8 ' -:. E 0.6....... 0.4............... -0.2 ' ' I' ' 3.6 3.8 4 4.2 4.4 4.6 4.8 5 Frequency(GHz) Figure 3.18: The normalized input impedance of a coaxial line-fed cavity-backed patch antenna. The parameters are the same as in Fig. 3.17. The FDTD data were contributed by Nihad I. Dib[40]. condition. 3.6 Conclusion In this chapter, the hybrid method presented in Chapter II is applied to the analysis of cavity-backed patch antennas. Examples of coaxial line-fed and microstrip line-fed patch antennas are given. The calculated reflection coefficient and input impedance are comared to the results calculated by the FDTD method. Good agreement is achieved. This result shows that the MoM/FEM technique is capable of analyzing various cavity-backed patch antenna structures.

44 3.5 GHz -- Xos = 3 mm Xos = 5.7 mm --- Xos = 9 mm Figure 3.19: g Comparison of the normalized input impedance of coaxial line-fed cavity-backed patch antennas with different feeding points. Other parameters are the same as in Fig. 3.17.

45 Total E-field (A) z I -1 L- "% xz (B) Y / 1 Y T >z (C) 0 0.2 0.4 0.6 0.8 1 a ~ b I I — - - - -- - a b Total E-field intensity in the cavity of Fig. 3.2 at resonant frequency (4.190 GHz). The parameters are the same as in Fig. 3.17. (A) and (B) are the cross section views along lines aa and bb shown in the inset, respectively. (C) is the view just beneath the patch. The black lines and white circles indicate the boundaries of the patch and coaxial line. This figure is drawn according to the scale of the physical dimensions.

CHAPTER IV Microstripline-fed Cavity Couplers and Resonators 4.1 Introduction Recent advances in Si micromachining techniques enable the fabrication of miniature metallic cavities in MMIC's. These cavities can be used as resonators, vertical transitions, or couplers. In contrast to conventional metallic cavities, micromachined cavities are easily integrated into planar circuits. They are also light in weight and easy to fabricate compared to conventional metallic cavities. In the following sections, the theoretical analysis of microstripline-coupled cavities based on the hybrid technique of Chapter II is presented. A design example of a broadband vertical coupler which achieves minimum insertion loss of 0.086 dB and a 80% 1 dB insertion loss bandwidth is given. A high-Q resonator was also designed and fabricated. Measured result indicates that a high-Q of 506 is achieved. 4.2 Single Cavity Formulation As show in Fig. 4.1, two slots are open on the top and bottom of a metallic cavity that is sandwiched between two layers of dielectric substrates. Either one of the microstrip lines passing the slots can be used as the input line or output line 46

47 of this structure. Without loss of generality, we choose the lower microstripline as the input line.. By closing the two slots and applying the equivalence principle as shown in Fig. 4.2, we can derive the following matrix equations from Chapter II and Chapter III: Vinc Zf 0 Jf Tf O M1 + (4.1) 0 0 Zt it 0 Tt -M2 Cinc Cf 0 iJ Y11 + Yf Y12 M1l 0 0 Ct Jt Y21 Y22 + Yt -M2 where Jf, Jt, M1 and M2 are the expansion coefficients of the unknown currents Jf, Jt, M1 and M, are as shown in Fig. 4.2. The meaning of the matrices is tabulated in Table 4.1. The above equation is related to (2.11) and (2.12) as follows: V = nc [Ci] = [cinc 0 0 [Jb] = [ ], [Ma] = Jt -M2 0 Zt 0 Ct [Ta] Tf 0 [ya] Yf 0 [] [ ] [y]_ [yc] = [_ 1:. Y21 Y22 By substituting the above equations into (2.21) and (2.22), the unknown coefficients Jf, Jt, Mi and M2 can be solved. Note that, we use the same basis functions on the microstriplines as described in Section 3.4. This allows us to derive the reflection and transmission coefficients of the input and output lines directly.

48 Microstrip line Slot on ground plane Microstri p l i ne 4.3 Multi cavity Formulation together and coupled through slots on their common walls. One way of analyzing this structure is considering the n-cavity structure as one big cavity and using the formulation of the previous section However, it is more efficient to perform FEM computation separately on each cavity instead of on the whole cavity. computation separately on each cavity instead of on the whole cavity.

49 Jt Microstripline Substrate -M2 e~~ M2 Cavity -0 M1 Substrate.T 'r.-r Microstripline inc f Figure 4.2' g The side view of the cavity coupler show in Fig. 4.1. Note that the slots on the cavity are closed by a PEC with an equivalent magnetic current. Jt Microstripline Substrate _Mn+1 Slot n+l Substrate Cavity n -M n + Jf Cavity n-l Slot n* J. M3 Cavity 2 Slot 2 M2 Substrate Microstripline Jinc + Jf Figure 4.3: The side view of a multi-cavity coupler.

50 Outside the Cavity Inside the Cavity Matrix Interaction between Vinc Jf Jinc Cinc Ml Jinc Zf Jf Jf Zt Jt Jt Tf Jf M1 Tt Jt M2 Cf M1 Jf Ct M2 Jt Yf Mi M1 Yt M2 M2 Matrix Interaction between Yll Mi M1 Y12 M1 M2 21 M2 M1 Y22 M2 M2 Table 4.1: The meaning of matrix in(4.1) and (4.2) Based on the equations in the previous section and the extra boundary conditions on each slot between the cavities, we can derive the following equations: Vinc Zf 0 Jf Tf 0 M 1 + =, (4.3) L 0 0 L ~ zt LJt. L 0 Tt -Mn+l [Cinc] + [Cf] [Jf = [Y11 + Yf] [M1] - [Y12] [M2], (4.4) [Y21] [M] -[Y22] [M2] = [Y22] [M2] -[Y23] [M3], (4.5) [Y32] [M2] -[Y33] [M3] = [Y33] [M3]- [Y34] [M4], (4.6) [Yi,i-i] [M-i]- [Y-,i] [lMi] = [Yi'i] [Mi]- [Yi,+] [M.+i], (4.7)

51 [Yn,n-1] [Mn-]- [M 1]- [Y,][M] = [Yn] [Mn] - [Yn,n+l] [Mn+l], (4.8) [Yn+l,n] [Mn]- [Yn+l,n+l + Yt] [Mn+l] = [Ct] [Jt], (4.9) where [Yi,j] is the interaction between Mi and Mj, [Yii] is the self interaction of Mi in its lower cavity, and [Yi'] the self interaction of MP in its upper cavity. Reassembling (4.3) to (4.9), we can solve this equations in a single formula: Zf -Tf 0 0 0.. 0 0 0 0 Jf -Vn Cf -Y1 Y12 0 0.. 0 0 0 0 M -Cinc 0 Y21 -'Y2 Y23 0. 0 0 0 0 M2 o 0 Y32 -YY3 Y4 " 0 0 0 0 M3 0 O 0 0 0 0 *. Ynn-1 -Yn Yn,n+l 0 Mn 0 0 0 0 ** 0 Yn+,n — Yn+l Ct Mn+l 0 0 0 0... 0 0 Tf Zt Jt 0 (4.10) where IY = Yi,i + YjI for i = 2, 3,..., n, Yi = Y11 + Yf and Yn+ = Yn+1,n+l + Yt. If the cavities and the slots are all the same, instead of solving (4.10), all the unknowns can be solved by using interaction of the single cavity structure. Assuming the cavities and the slots are all the same, we have [Y,,-1] = [Yi,,+i] - A, A, (4.11) [Yi,i] [Yi,] = B, for i = 2,...,n. To simplify these equations, we use Mi instead of [Mi]. (4.7) becomes AMi-1 - BMi = BMi - Ai+1,

52 = Mi1 - B'Mi = B'i-Mi+, (4.12) where B' = A-1B. Assume (4.12) can be recast into the following recursive form: Mi-1 - bMi = c(Mi -- bMi+l). (4.13) We find c and b satisfy b+c=2B' and cb=l1. (4.14) Solving (4.14), we have b = B't + B'2-, (4.15) (4.15) c B'~ B'2-1. Note that it maybe difficult or impossible to explicitly derive the square root of a matrix. However, we will show that the final equations will not involve the square root terms in (4.15). Substitute the recursive relationship to (4.5)-(4.8), we can derive the following formula: M1 - bM2 c(M2- bM3) (4.16) c2(M3 -- bM4) (4.17):: (4.18) cn-1 (Mn - bMn-l ). (4.19) From (4.15), exchanging c and b will still satisfy the recursive relationship, such that we also have M1 - cM2 bn-l (M - cMn-l).. (4.20) Solving (4.19) and (4.20) simultaneously, we can derive M2 and Mn in terms of M1

53 and M,+i as follows: (bn-I-cn-1)Mj +(b-c)Mn+l = aM + OM, +l, M = (b-c)Ml+(bn-l-cn-l)Mnl 1 + M(4.21) Mn= _________ = /M,+ aMn+,. Note that it seems the above solution is not unique because b and c have two sets of solutions. However, after carefully examining the form of (4.21), we find they are symmetrical to b and c and either set of the solution of b and c gives the same result. Substituting (4.21) to (4.14) and (4.15), we have Cine 7cf 0 Jf Y1 - Y12a -Y12/3 M + =. (4.22) 0 0 Ct Jt -Yn+l,n/ Yn+l - Yn+l,na Mn+i With (4.22) and (4.3) the unknowns Jf, Jt, M1 and Mn,+ can be solved easily as in the case of single cavity. Other unknowns can be derived from M1 and Mn+1 recursively. Note that in (4.21) a and,/ are evaluated in terms of c and b. However, c and b may be difficult to derive because of the square root terms in (4.15). Instead, we can derive a and /3 in terms of b + c and bc. By taking out the common factor b - c in (4.21), a and 3 become Sn-2 (4.23) Sn-i1 /: (4.24) Sn-1 where n Sn= bicn-i. (4.25) i=o If We let S-i = 0 and So = 1, we can derive the following recursive relationship for Sn= (b+ c)S,- - bcSn-2 for n > 1. (4.26)

54 With (4.23), (4.24), (4.26) and (4.14), (a and 3 can be derived solely from the A and B in (4.11) without solving b and c directly. To check the validity of these equations, set n = 1. We find a = 0 and / = 1. After Substituting this result into (4.22), this equation becomes the same as (4.2) which is the single cavity formulation. 4.4 Validation In this section, the computed result using the single-cavity formulation is compared to the computed and experimental data of [34]. In that paper, a microstripline coupler with thick ground plane as shown in Fig. 4.4 was investigated. The cavity was treated as a rectangular waveguide and only the TE1o mode was used in its MoM formulation. Outside the cavity, quasi-static Green's functions were used in its MoM formulation. Fig. 4.5 shows the comparison between the paper's result and the computed result of the hybrid technique for three different ground plane thicknesses. The dimensions of the structure are as indicated in Fig. 4.4. Good agreement is achieved for thicknesses of 300 mil and 1000 mil. In the 50 mil case, the result of the hybrid technique is closer to the experimental data than the computed data of [34]. Since 1000 mil and 300 mil are multiples of 50 mil, we can use the multi-cavity formulation to derive the computed result of 1000 mil and 300 mil cases from 50 mil case. Fig. 4.6 shows the result of considering the 1000 mil and 300 mil cavities as twenty and six 50 mil cavities stacked together respectively. The computed data is almost the same as the single cavity model except for a small deviation of IS,, I near resonance.

55 Top view 100 196 ~ -196 Microstripline Side view 62 1000 62 Unit: mil 650 650 Figure 4.4: A microstrip line coupler as described in [34]

56 50 mil Cf2 CI~ql -14 -1: Freq(GHz) 300 mil -. C. -. Freq(GHz) Figure 4.5: Comparison of the result from [34] and the hybrid technique of this dissertation. In the legend, EXP and MoM are the experimental and computed results of [34], and MoM/FEM is the computed result using the hybrid technique.

57 1000 mil Freq(GHz) Figure 4.5: Continued from previous page.

58 Validation of Multi-Cavity Formulation 6 Validation of Multi-Cavity Formulation 4 Freq (GHz) 6 Figure 4.6: Comparison of the computed result of the single cavity and multi-cavity formulations. The solid line is the data computed by the single-cavity model and the dashed line is the data computed from the multi-cavity model.

59 OF -101 -20 /' ~-30 '-4 Validation of Multi-Cavity Formulation - 300mil --- 50milx6 i I!I -40 -50 ICA -r"I X. - 3 4 Freq (GHz) 5 6 Validation of Multi-Cavity Formulation 4 Freq (GHz) 6 Figure 4.6: Continued from previous page.

60 4.5 Wide-band Cavity Couplers/Vertical Transitions Conventionally, via holes are used as a vertical transition between different layers of a MMIC. As the operating frequency increases, the fabrication of the via holes becomes very difficult. To circumvent this problem, non-direct-contact vertical transitions are necessary. In [51], broadband vertical interconnects using slot-coupled microstrip lines are presented. In that paper, two microstrip lines are coupled through a slot on their common ground plane. A 3:1 bandwidth of 1.0 dB insertion loss is achieved. In [34], a similar design is presented except that the thickness of the common ground plane of the two microstriplines is considered. A thick ground plane can be used as a heat sink or a part of the housing of the circuits. A similar design employing metallic cavity as the transition between the microstriplines is presented in Fig. 4.7. The dimension of the structure is chosen according to the feasibility of current manufacturing processes. The computed results indicates a minimum insertion loss of 0.086 dB at 72 GHz, and a 1 dB bandwidth from 54.1 GHz to 126.3 which is amount to 80%. Due to the manufacturing process of the micromachined cavity, errors in the dimension and alignment are inevitable. Sensitivity study on the dimension of the structure is necessary to ensure the feasibility of a practical design. Fig. 4.8 and Fig. 4.9 show the effect of possible errors in the width of the cavity and in the positions of the slots which may be caused by imprecise alignment. Both computed results indicate the performance of the coupler is very stable under the above error range. This is because the change in dimension is small compared to the wavelength in the operational frequency range. Fig. 4.10 shows the minimum insertion loss achievable v.s. the height of the

61 Top view 0.0725 g 0.088 0.0725 Microstripline,4-~~:;:::~~~~~j:2~:~' Side view i 0.09 0.1;0.09! 0.279 (a) 60 80 Freq (GHz) 140 (b) Figure 4.7: (a)the structures of a microstripline vertical coupler. The cavity size is 0.0968 x 0.968 x 0.10. All unit in mm.(b) computed S-parameters.

62 cavity. Each data point is computed by adjusting the other dimensions to achieve a minimum insertion loss while keeping the height of the cavity fixed. The figure shows that the insertion loss increases as the height increases.

63 O -'I -10 -20 (-30 -40 - Cavity width 0.0968mm - -- Cavity width 0.138mm.... Cavity width 0.188mm "it~l I... 20 40 60 80 Freq (GHz) 100 120 140 (a) II.. -0.1 0 CL4 -. -Cavity width 0.0968mm - - - Cavity width 0.138mm........Cavity width 0...188m... Cavity width 0. 188mm -0.4k i 6-~( i _i i i ) 70 80 90 Freq (GHz) 100 110 120 (b) Figure 4.8: g Sensitivity study of the cavity in Fig. 4.7. Three cavity widths 0.968 mm, 0.138 mm and 0.188 mm are compared. The data shown in (b) is the same as (a) except that the scale of the figure is different.

64 60 80 100 Freq (GHz) (a) 90 100 Freq (GHz) 120 (b) Figure 4.9: Sensitivity study of the cavity in Fig. 4.7. Three slot positions are compared. Curve 1: both of the slots are at the center of the cavity. Curve 2: both of the slots are offset by 0.01 mm to the center of the cavity. Curve 3: the same as the second except that the offset of one slot is in the opposite direction to the other. The data shown in (b) is the same as (a) except that the scale of the figure is different.

65 0... -.... m -H n 4a) N -0.1 -0.2 -0.3 -0.4 ", 10, I I I I I I I I I I -0.5 ' 0.1 I..... I. I 0.15 0.2 0.25 0.3 0.35 Cavity Height 0.4 0.45 0.5 0.55 (mm) Figure 4.10: Optimized insertion loss v.s. cavity height.

66 4.6 High-Q Micromachined Cavities Conventional microwave high-Q resonators made by metallic rectangular or cylindrical waveguides are heavy in weight, large in size and costly to manufacture. It is also difficult to integrate these resonators into MMIC's. Extra circuitry is alway necessary to connect the resonators to other solid state devices, resulting in increased complexity and cost. With the maturity of micromachining techniques in fabricating microwave circuits, it is now possible to make miniature waveguides or cavities [35]-[38] as building blocks of high-Q bandpass filters. The quality factor that can be achieved with this technique is much higher than the quality factor of traditional microstrip resonators either printed on a dielectric material or suspended in air with the help of a dielectric membrane [39]. Fig. 4.11 shows possible high-Q filter geometries that consist of input and output microstrip lines and rectangular cavities on different dielectric layers. The cavities are made by excavating and metallizing the dielectric using micromachining techniques as described in [41]. Coupling between the cavities and microstrip lines are achieved via the slots etched at appropriate locations with respect to the microstrip lines. Coupling between cavities is controlled by the size, position, and orientation of the corresponding coupling slots. The microstrip lines can be laid on the same plane as in Fig. 4.11-(a), or on different planes to serve as a vertical transition as in Fig. 4.11-(b). The vertical stacking of the cavities greatly reduces the occupied area when multiple cavities are needed in the filter design. In the following sections, a detailed study of the resonator is presented. A highQ micromachined resonator was designed and built. The computed and measured results are compared.

67 Top View Side View Metalli; Top View Microstr Side View Metall Mirostrip Lines Coupl Slots Microstrip Lines Coupling Slots Microstrip Lines (a) Patent Pending iI L. P. B. Katehi J. Cheng J. Papapolymerou iIl Microstrip Lines Coupling Slots Microstrip Line (b) Figure 4.11: The structures of the proposed micromachined bandpass filters. (a) input and output microstrip lines on the same substrate. (b) input and output microstrip lines on different substrates.

68 4.6.1 Design Consideration Fig. 4.12 shows a design of a square cavity resonator coupled by two slots at the center of its top and side wall. Fig. 4.13 shows the computed results of this structure with different cavity sizes. The dimensions are as shown in Fig. 4.12. Note that, the first resonant modes of the cavities in Fig. 4.13 are at 16.2 GHz, 19.1 GHz and 23.3 GHz. However, from the data we can see no resonant behavior in the frequency range. The reason is because the two slots are too close. This causes the higher order modes excited by the slots to be strongly coupled. To avoid this phenomenon, the two slots should be separated by a sufficient distance. However, to efficiently excite the first resonant mode, the slots should be near the center of the cavity. Fig. 4.14 shows a compromised design in which the slots are located at 1/4 and 3/4 of the cavity length. Fig. 4.15 shows the computed IS11I, IS121 and the loss defined by 1 - IS 112 - IS122. Since the dielectric and conductor are assumed to be lossless, the only contribution to the loss is the radiation of the two feedline layers. Three resonances are observed in the frequency range shown in the figure. The 3 observed resonant frequencies are very close to those of TM1o, TM130 and TM150 to z modes of the cavity as listed in Table 4.2. Note that, because of the symmetry of the structure, only odd modes are excited. Fig. 4.16 - 4.18 show the field plots inside the cavity these resonant frequencies. From these field patterns, we can see clearly the corresponding cavity modes. Fig. 4.15 also shows a peak loss at 28.6 GHz. Comparing the stub length and the guided wavelength of the microstripline at this frequency, we find that the stub length is equal to half the guided wavelength of the microstripline such that the stubs become resonate dipoles and strongly radiate. As long as the microstriplines are kept the same, this peak loss at 28.6 GHz is always observed for different cavity

69 Top View H I I I.. - Slot 0.635mm x 6.35mm Side View Microstripline Dielectric 10.6 Thickness 0.635 /7 Microstripline, width 0.568mm, stub 1.778mm Figure 4.12: A square cavity resonator with both feeding slots at the center. 1 I...i iI 0.9 9x9x0.5 Unit: mm lxllxO0.5 ----... 0.8 13x13x0.5 -—. --- — 0.7 0.6 0.7.-" / -5 0.4 0.2 0. - 4 6 8 10 12 14 16 18 20 Frequency (GHz) Figure 4.13: The IS11I of square cavities. The cavity sizes are shown in the legend. Other parameters are the same as in Fig. 4.12.

70 0.635 0.635 0.404 1_ 7 0.404 1 778 1 778:!... Unit: mm 32 ' 8 ' 16, 8 0. 5 0.5 Figure 4.14' Figure 4.14: A rectangular cavity resonator with feeding slots separated by a distance. sizes. To investigate the effect of multi-cavity configurations, consider stacking the same cavity of the previous paragraph. Fig. 4.19 shows the result of single, two and threecavity configurations around the first resonant frequency at 19.38 GHz. We find that the first resonance is split into two and three resonances in two-cavity and threecavity configurations, respectively. This phenomenon is expected because when N identically resonant systems are weakly coupled, the resonant frequency will split to N slightly different ones. Also notice that the resonant curves become sharper as the number of cavities increases. 4.6.2 A High-Q Resonator In this section, a high-Q resonator is designed and fabricated. For ease of fabrication, the feedlines are placed on the same dielectric layer. Due to the fabrication process, the angle between the vertical walls and the top and bottom walls is no

71 Coplanar Coupler. Cavity: 8mm x 32mm x 0.5mm.00.0 10 0ftwo 4) 3: 5 10 15 20 25 30 35 Freq (GHz) Figure 4.15: The computed results of the structure shown in Fig. 4.14. longer 90 degrees, but 54 degrees [41]. Although the cavity is no longer rectangular, its characteristics are still very close to the rectangular one, because the thickness of the cavity is much smaller than its width and length. The coupling between the two coplanar microstriplines also has little influence on its characteristics. Fig. 4.20 is provided as evidence of the above statement. A resonator with the dimensions shown in Fig. 4.21 was builtl. The computed Sparameters are shown in Fig. 4.23 and the field plot at the first resonance is shown in Fig. 4.22. Four transmission peaks are observed at the frequencies listed in Table 4.3 with comparison to the resonant frequencies of an ideal rectangular cavity of the same size. Compared to the cavity in the previous section, the discrepancy between 1Fabrication and measurement were performed by John Papapolymerou

72 7 Cavity size: 8mm x 32mm x 0.5mm No. Observed Resonant Freq. Resonant Mode Mode 2. 1 19.38 GHz TMllo 19.33 GHz TM120 20.96 GHz 2 23.92 GHz TM130 23.44 GHz TM140 26.52 GHz 3 30.40 GHz TM150 30.01 GHz TM160 33.80 GHz Table 4.2: The comparison of the observed resonant frequencies and the frequencies of cavity resonant modes. the observed resonant frequencies and those of the ideal rectangular cavity is larger. Note that, the cavity size in the previous section is 8mm x 32mm x 0.5mm compared to 16mm x 32nzm x 0.465mm in this section. This results in a lower decaying rate of the high order modes excited by the slots such that the coupling between the slots is not only dominated by the cavity modes but also slight influenced by these higher order modes. Quantitatively, consider the cavity as a waveguide with width a and length b. Let a < b. The propagation constant f of the second waveguide mode at the resonant frequency of the cavity is n2 W2 22r2 a2 2 (4.27) 3a2 r2 a2 b2 (4.28) Since a < b, '3 is pure imaginary and decreases as a increases. Thus, the second

73 Cavity size: 16mm x 32mm x 0.5mm No. Observed Resonant Freq. Resonant Mode Mode Freq. 1 10.54 GHz TMl1o 10.48 GHz TM120 13.26 GHz 2 18.29 GHz TM130 16.90 GHz TM140 20.96 GHz 3 23.00 GHz TM15o 25.24 GHz TM160 29.64 GHz 4 33.95 GHz TM170 34.13 GHz Table 4.3: The comparison of the observed resonant frequencies and the frequencies of cavity resonant modes. waveguide mode decays faster in the cavity of the previous section than in the cavity of this section. The measured results near the first resonance is shown in Fig. 4.24. The fabrication process is as described in [41]. The measurement shows that a minimum insertion loss of 0.36 dB and 5% 3-dB bandwidth are achieved. The unloaded Q of this cavity is 506 computed as in [39, 41], which is close to the theoretical value of 526 for a metallic cavity of the same size. 4.7 Conclusion In this chapter, the hybrid MoM/FEM method is successfully applied to a microstriplinefed slot-coupled cavity. Depending on the dimension and operating frequency, the cavity can serve as a vertical transition between two substrate layers, or as a high-Q

74 resonator. A wide band vertical coupler is designed and simulated by the hybrid technique. The computed results indicate a 80v% bandwidth and a minimum insertion loss of 0.086 dB. A micromachined resonator is also designed and fabricated. The computed results are in good agreement with the measurement. A high Q of 506 is achieved.

75 Bottom Top I 0 0.2 0.4 0.6 0.8 1 Maximum= 2675554 Cavity size: 8mm x 32mm x 0.5mm z-component of E field, at 19.38 GHz Figure 4.16: The z component of the electric field inside the cavity at the first resonant frequency 19.38 GHz. The figure is drawn to scale. From left to right columns are the field profile on x-y plane at the top and the bottom of the cavity, on x-z plane at the center of the slots and the cavity, on y-z plane at the center of the cavity. The dimension of the structure is as in Fig. 4.14.

76 Bottom Top 0 0.2 0.4 0.6 0.8 1 Maximum= 1.777767e+06 Cavity size: 8mm x 32mm x 0.5mm z-component of E field, at 23.92 GHz Figure 4.17: g The z component of the electric field inside the cavity at the second resonant frequency 23.92 GHz. Refer to Fig. 4.16 for more detail.

77 Bottom Top - - - X --- 0 0.2 0.4 0.6 0.8 1 Maximum= 1.902762e+06 Cavity size: 8mm x 32mm x 0.5mm z-component of E field, at 30.40 GHz Figure 4.18: The z component of the electric field inside the cavity at the third resonant frequency 30.04 GHz. Refer to Fig. 4.16 for more detail.

78 i 0.8 0.7 0.6 /o f. two-cavity 0.2 r ree-cavity -- 0.1 - 0. 18.6 18.8 19 19.2 19.4 19.6 19.8 20 Frequency (GHz) Figure 4.19: The comparison of S11 I of single, two and three cavity configurations. The dimension of the structure is the same as in Fig. 4.14

79 -4q cr-4 cn 1 0.8 0.6 0.4 0.2 0 5 10 15 20 25 30 35 Frequency (GHz) Figure 4.20: The comparison of the effect that; cavity shape has on the coupling between the microstriplines. The dimension of the structure is as in Fig. 4.14. Curve 1: coplanar feedlines, trapezoidal cavity. Curve 2: feedlines on different sides, rectangular cavity. 0.635 0.635 ----- --- —..................-..~-..... 16.354 i. 0.410 7 | 7 0.410: 1.778 1.778:, Unit: mm 32.354 8.177: 16 8.177' o0.5 I 0.465 Figure 4.21: A high-Q resonator fabricated by micromachining[41].

80 Figure 4.22: g The z component of the electric field inside the cavity at the first resonant frequency 10.54 GHz. The figure is drawn to scale. The dimension of the structure is as in Fig. 4.21.

81 -e -20 9 9 -25 9 9 9 9 9 o..' ' '9| '.' - -30- l ' - -40-, -IS1 21 -45 - Loss/, o 50' " ' ' 5 10 15 20 25 30 35 Frequency (GHz) Figure 4.23: The computed S-parameters of the resonator shown in Fig. 4.21 -35/ -10 7 8 9 111 12 13 7 8 9 10 11 12 13 Freq (GHz) Figure 4.24: The comparison of measured[41] and computed S-parameters of the resonator in Fig. 4.21 near the first resonance.

CHAPTER V Wavelet-Based Method of Moments 5.1 Introduction In recent years, we witnessed an explosion of interest in wavelet theory in various applications. For instance, wavelet theory has been applied to compression of digital images, noise reduction in signal processing, and numerical methods in electromagnetics [44]-[48]. Wavelet theory is characterized by a special group of basis functions used to expand a linear space. When combined with the method of moments, these basis functions produce a matrix with many small off-diagonal terms that can be treated as zero without seriously reducing the accuracy of the solution. Thus, the storage space occupied by these zero terms can be removed by sparse matrix storage schemes so that the memory size required to solve the problem is reduced. This feature is important for solving large scale problems like the full-wave simulation of a finite patch antenna array or a complex microwave circuit. Although computer memory size is growing rapidly, it is still a challenge to run a full-wave simulation of electromagnetic problems in the scale of a few wavelengths. In the following sections, the theory of multiresolution expansion and FWT are introduced. The thresholding effect on the characteristics of a microstrip patch antenna is investigated. 3 x 3 and 6 x 6 arrays are also used as large scale problems 82

83 to show the benefit of wavelet bases. 5.2 Multiresolution Expansion From the point of view of linear algebra, wavelet theory[49, 50] employs a special group of basis functions to expand a linear space. These basis functions consist of two sub-groups: scaling functions and wavelets. The scaling function represents the local average. In contrast, the wavelet represents the local variation and possess zero average. These two groups of functions are orthogonal to each other and form a complete basis of the linear space. Let So and Wo denote the linear space expanded by the base scaling functions and the associated wavelets respectively. Suppose we can define a contraction operation to transform these basis functions to a higher level of resolution in the solution space we seek. To qualify as a multiresolution expansion, the space expanded by the lower level scaling functions must be a subset of the space expanded by the higher level scaling functions. Also the direct sum of the space expanded by the same level of scaling functions and wavelets are equivalent to the space expanded by scaling functions one level higher. That is So C Si C - CSmr, and (5.1) Sm+l = Wm D Sm (5.2) Wm (D Wm-1 & Sm-1 (5.3) (5.4) = Wm Wm-i ) Wo So, (5.5) where Si and Wi are the space expanded by i levels of contraction of the base scaling functions and wavelets respectively, and & is the direct sum operator as defined in linear algebra. As an example, Fig. 5.1 shows a set of scaling functions consisting of

84 s"n sn-1. sn-2...... S wn-1 wn2. *W~ Scaling functions -.T- i --- — - \-I.w *I''...... Wavelets fftf - fti g The concept of multiresolution expansion. S' and Wn are the linear space expanded by n-th level scaling functions and wavelets respectively. The lower half of the figure illustrate the concept of multiresolution expansion by pulse functions and their wavelets. 8 pulses. The space expanded by these scaling functions is equivalent to the space expanded by its lower level scaling functions and wavelets. In the method of moments, each matrix element represents the interaction between two basis functions defined by an integration with a kernel. Except self interaction terms, the matrix elements tend to be very small if at least one wavelet is involved in the integration. This is because in most of the cases the kernel function is smoother than the locally fast-changing wavelets. Since the wavelets have zero average, the resulting integration tempts to cancel itself and become very small. From the multiresolution expansion, the space expanded by a higher level scaling basis functions is equivalent to that expanded by a set of base scaling functions plus a series of wavelet bases as described in (5.5). By using the later as the basis in the moment method, the resulting matrix will have many small matrix elements which can be neglected without adversely degrading the solution. By using sparse matrix storage scheme, the memory required to store the matrix will be dramatically

85 reduced. 5.3 Moment Method Formulation of a Patch Antenna In the moment method modeling of a patch antenna, the triangle rooftop functions are frequently used as the basis functions of the induced current on the patch as shown in Fig. 5.2. Along the direction of the current, the solution space is expanded by those triangle functions in Fig. 5.2-(a), which satisfy the requirement of zero current at both ends of the patch. Along the transverse direction, the solution space is expanded by those pulse functions in Fig. 5.2-(b), which can approximate the infinite current density near the two edges of the patch. Compared to B-spline scaling functions defined in [52], the pulse basis corresponds to the zero degree B-spline scaling functions as shown in Fig. 5.3-(a). Thus, we can use their associated wavelets (Fig. 5.3-(b)) directly. However, the triangle basis is not exactly the same as the first degree B-spine scaling functions shown in Fig. 5.3 -(c) which includes half triangle bases at both ends of the problem domain. The linear space expanded by this basis no longer satisfies the condition of zero current at the ends of the patch, thus the associated wavelets can not be used directly in the multiresolution expansion of our patch problems. Instead, the edge wavelets are modified to have zero value at the two ends of the problem domain as shown in Fig. 5.4. Note that the modified edge wavelets still have zero average, but are not orthonormal to the scaling functions anymore. The latter feature is not important in our problem since the matrix elements are not the direct integration of two basis functions, but with a kernel. To summarize, zero degree B-spline scaling functions and wavelets can be used directly to model the transverse current profile on a patch. To model the longitudinal

86 current profile, we use modified first degree B-spline scaling functions and wavelets. The scaling functions and wavelets are the same as the original ones except the edge scaling functions are removed and edge wavelets are modified such that their values are zero at the edges. The two-dimensional rooftop basis and its associated wavelets[53] are formed by direct multiplication of the zero and first degree B-spline scaling functions and wavelets as shown in Fig. 5.5. These basis functions are used in the following derivation of the method of moments. Let J1, J2,... J., be basis functions. Let Jp be the induced current on the patch. Jp can be expanded by the basis as follows: n Jp = E amJm (5.6) m=1 where am is the unknown expansion coefficient. Suppose the incident field is E, then the total field Et can be represented as follows Et(r) = JG(rlr') Jp(r')ds', (5.7) SI where G(rlr') is the electric field dyadic Green's function of a conductor-backed dielectric substrate as described in Chapter III. Since the total tangential electric field on the patch must vanish, we have 0 = fJ Jm (r) G(rr') Jp(r')ds'ds, for m = 1,...,n (5.8) S SI Substitute (5.7) to (5.8) and rearrange, we reach the following matrix equation: [Z][I] = [V], (5.9) where [Z]ij = Ji(r) * G(rlr') Jj(r') ds' ds, n x n matrix, (5.10) S 5'

87 Direction of current The triangle basis functions AA Ground Plane (a).................. -......................f..... f;0.................... Transverse current profile The pulse basis functions-i-........4................X.... 0.........................................,................. -....,..... t. 00:. —. 4.........................X t......................................................................................................,.... 0......................... a. 1.............................................. l an.............. ( a )........................ c u r e n..........o f.i......l..................................................... Ground Plane (b) Figure 5.2: The rooftop basis functions used to expand the current on a patch. (a) Profile along the direction of the current. (b) Profile along the transverse direction of the current. [I]i ai [V]i =- J Ji(r). E(r) ds, s n - element vector, n - element vector. (5.11) (5.12) (5.9) can be solved by direct matrix inversion when the matrix dimension is small, or by iterative algorithm such as bi-conjugate gradient when matrix dimension is large. 5.4 Fast Wavelet Transform Due to the convolutional property of the Green's functions, The integration of (5.10) only depends on the distance between the two basis functions when they have

88 I I I I (a) (b) (c) (d) 0 1 I I I I AI 0 i v V V \1 Figure 5.3: (a) zero degree B-spline scaling functions, (b) zero degree B-spline wavelets, also called Haar wavelets, (c) first degree B-spline scaling functions, (d) first degree B-spline wavelets.

89 -1 1 i 1 Scaling Function 1/12 %,T A-l sr — s. / 5/6 1/12 Normal wavelet -1/2 -1/2 5/12 1/12 I I / \. 4. B Edge Wavelet -1 V -1/2 Figure 5.4: g The modified first degree B-spline scaling functions and wavelets. The scaling functions consist of only triangle functions without the two half triangles basis at the ends. The wavelets are the same as the original except the edge wavelets have zero value at the end. Only one edge wavelet is shown in figure since the other one is just a mirror.

90 (a) (b) (c) (d) Figure 5.5: The two-dimensional rooftop basis functions and its wavelets[53]. (a) the scaling function formed by the multiplication of the zero and first degree B-spline scaling functions, (b) a wavelet formed by the zero degree scaling function and the first degree wavelets, (c) a wavelet formed by the zero degree wavelet and first degree scaling function, (d) a wavelet formed by the the wavelets of zero degree and first degree.

91 Expansion Level 0 1 2 3 Number of Unique elements 705 3568 8876 13097 Table 5.1: Comparison of unique elements v.s. expansion levels. the same form. Thus, when only one kind of scaling function is used in the aforementioned patch problem, the elements of [Z] will only have n different values, where n is the dimension of the matrix. Instead of performing n x n times of integration, only n integration are needed, which results in a great saving in computation time. As discussed in Section 5.1, we expect the use of wavelets as part of the basis functions will have the benefit of reducing the memory requirement. However, the introduction of wavelets will ruin the Toeplitz property of the matrix, which increases the computation time. As an example, a patch antenna is discretized as a 16 x 16 grid. The number of basis functions in both x and y directions are 15 x 16 respectively. Table 5.1 shows the number of unique elements in different expansion levels. From the table, we see a 5 time increase in the number of unique elements when merely changing the basis from scaling functions only to one level expansion which compromises approximately the same number of scaling functions and wavelets. Here comes the dilemma. We want to use wavelets to reduce the memory requirement. However, doing so will increase computation time. To avoid this, instead of using the wavelet basis directly, we perform the integration of the matrix elements using higher level scaling functions. Since these functions span the same linear space expanded by the lower level scaling functions and wavelets as shown in (5.5), there must exist a linear transform between them. V = {lov +25 * * * v 4 } = {li v2,. *,, On} Let 1 be the basis only consisting of

92 higher-level scaling functions, and 4' be the basis consisting of the lower-level scaling functions and wavelets. Assume {(1, (2,..., I n} = {1, ~)2***''' n} (5.13) If matrix [P] is the transformation matrix between ( and I', we have n '= ij'k., (5.14) j=l where pij is the i-th row and j-th column element of [P]. Let [Z] in(5.10) be the matrix evaluated in basis I and [Z'] evaluated in basis V', then we have [Z'] = [P]t[Z][P] (5.15) where t is the transpose operator. Since the relationship of scaling functions and wavelets are well defined, [P] is easy to derive. 5.5 Microstrip Patch Antenna In this section, the wavelet theory described in previous sections will be applied to the analysis of a microstrip patch antenna as shown in Fig. 5.6. Without loss of generality, we assume a normal incident field Ei which can be expressed as follows E(r) = (1 + F)e(k~z)y, (5.16) where F is the reflection coefficient at the interface between the substrate and free space and ko is the propagation constant of free space. Since the patch current includes both x and y components, Jp is expanded as follows Jp(r) = Jx(x, y)x + Jy(x, y)y. (5.17) Jx and J. are expanded by the scaling and wavelet basis discussed in previous sections.

93 Figure 5.6: A microstrip patch antenna Each element of i.. mpedance matrix integrated by an adaptive Gaussian matrix has at leas. t fiive-digit accuracy.ii A uniform grid of 16 x 16 was used to discretize the patch. The grid size is equivalent to 15 triangle functions ti ir ieo the current and 16 pulse functions in thee tanvre direction as shown, in Fig. 5.7-(a). Let th is set of sicaling functions be th e highest-level resolution, and called them a zero-level multiresolution expansion. Then, the one-level multiresolution expansion contains the scaling functions and wavelets at the next lower-level resolution, i.e. 8 x 8 grid as shown in Fig. 5.7-(b). Similarly, the two-level multiresolution expansion contains scaling func

94 (a) (b) (c) Figure 5.7: Multiresolution grids at different resolution levels: (a) 16 x 16, (b) 8 x 8, (c) 4 x 4. tions and wavelets defined on 4 x 4 grid as shown in Fig. 5.7-(c) plus the wavelets in the previous resolution. Fig. 5.8 shows the dominant component of the patch current distribution near the resonant frequency of the patch in Fig. 5.6. The result of using the bases of the three multiresolution expansions is exactly the same as the theory predicted. To investigate the effect of thresholding, the moment matrix produced by the level-one basis is examined. Those elements smaller than a given percentage of the maximum element are set to zero. Then, the moment matrix is inverted to solve the unknown expansion coefficients. Fig. 5.9 shows the resulting current distribution after applying thresholding levels of 0.001%, 0.01%, 0.05% and 0.1% respectively. As expected, the larger the thresholding level, the more seriously the current distribution deteriorates. At a level of 0.05%, degradation of the current distribution is already obvious. Fig. 5.10 shows the structure of the moment matrices after thresholded at the aforementioned levels. Non-zero matrix elements are marked by a black dot in the figure. Sparsity values of 65%, 76.8%, 80.4% and 81% are obtained respectively for

95 the four cases. The sparsity value is defined by dividing the number of zero elements by the total number of the matrix elements. Note that the elements involving the interaction of two scaling functions are excluded in the thresholding process. In order to show the strength of the interaction between two basis functions, Fig. 5.11 plots the absolute values of the elements of the impedance matrix in Fig. 5.10 before thresholding v.s. the distance between the two basis functions. The figure is categorized by the types of basis functions shown at the right hand side of the figure. Only 4 types are plotted as examples. We found that the maximum value of the impedance matrix is 0.02044 which is produced by the wavelets shown in figure (b). At a, distance of 0.025 meter, or equivalently, 0.128 free space wavelength, the value drops to less than 0.001% of the maximum value. Note that in (c) the maximum value produced by the interaction of a scaling function and a wavelet is much smaller than the other three cases. The maximum value in (b) is larger than (a). This is because the wavelet spans an area 1.5 times of the area spanned by the scaling function. As expected, (b) and (d) show that the interaction decays faster when only wavelets are involved. Fig. 5.12 shows the comparison of the sparsity of the moment matrix generated by the one-level and two-level multiresolution expansions. Contrary to what we expect, the sparsity of the latter case is not always higher than the former one. One possible explanation is that the larger cell size of the lower level wavelets overlays with other functions, thus reducing the canceling effect of the wavelets. Fig. 5.13 shows the magnitude of the dominant current at the center of the patch for different threshold levels. The resonant frequency can be determined by the peak of the curve.. From the figure, very accurate results are obtained by one-level expansion even at a 0.05% threshold level. However, in the two-level expansion, the

96 curve is already seriously deteriorated at the same threshold level. To investigate the effect of thresholding on the far-field pattern, the radiation pattern of the patch antennas at different threshold levels are computed from the current distribution on the patch. Fig. 5.17 shows the un-normalized far-field pattern for the case of a one-level expansion. Compared to the current distribution, the pattern is more immune to the error caused by thresholding. The result shows very good accuracy for up to a 0.1% threshold level. 5.6 Patch Antenna Arrays In this section, based on the single patch of previous section, a 3 x 3 and a 6 x 6 arrays are studied to explore the potential of multiresolution expansion in large scale problems. Fig. 5.14 shows the configuration of the arrays. By default, the distance between adjacent elements is equal to half the free space wavelength in both x and y directions. Using the same resolution as the single patch, the number of unknowns is 4320 and 17280 respectively. 5.6.1 3 x 3 Array Fig. 5.15 shows the comparison of the sparsity of the one- and two-level expansions at different threshold levels. Unlike the single patch, the sparsity of the two-level expansion is always higher than the one-level expansion at the same threshold level. The reason is that the majority of the matrix elements are from the integration of the basis functions on different patches such that the canceling effect of wavelet bases is more obvious. Note that, at 0.01% threshold level, the two-level expansion already gives a high sparsity of 96.1 %. Fig. 5.16 shows the structure of the sparse matrix of this case. Fig. 5.17 shows the un-normalized far-field pattern of the array. The "array factor" shown in the legend of the figure means that the pattern is

97 computed by multiplying the array factor with the pattern of a single patch without thresholding. The figure shows that the pattern computed by the array factor has a significant difference in the power level with the other patterns. This discrepancy may be caused by the neglect of the coupling between patch elements when using the array factor to compute the pattern. At all the threshold levels shown in the figure, good agreement is obtained in the main beam. However, noticeable differences are seen in the side lobe levels. exhibits larger difference. 5.6.2 6 x 6 Array In this section, a 6 x 6 array is studied. Fig. 5.18 shows the comparison of the sparsity of the one- and two-level expansions at different threshold levels. Similar to the results of the previous section, the sparsity of a two-level expansion is alway higher than the case of a one-level expansion. The difference is even larger than in the previous section. At a threshold level of 0.01%, the two-level case already achieves a very high sparsity of 98.84%. The structure of the sparse matrix of this case is shown in Fig. 5.19. Fig. 5.20 shows the normalized far-field pattern of the array. At all threshold levels, the pattern changes very little. This indicates that the radiation pattern is highly insensitive to the error incurred by the thresholding process. Table 5.2 shows the memory size of the sparse matrix at different threshold levels. Note that, the full matrix storage requires 4778 MBytes, which is prohibited by most of today's computers. However, even at the lowest threshold level of 0.001%, the memory size reduces to only 118 MBytes.

98 12 T 1 ^ 0.8I 0.4 02 0 so 0 1so rc —\ 00 Y axis (mm) X axis (mm) gre 5.8: The dominant component of the induced current on the patch antenna. The parameters used in the computation are a = 5.95 cm, b = 5 cm, d = 0.159 cm, and fc = 2.55. The definition of these symbols are shown in Fig. 5.6 Threshold level Sparsity Memory size (Bytes) None Full matrix 4778 M 0.00001 98.35% 118 M 0.0001 98.84% 83.1 M 0.0005 99.10% 64.5 M 0.001 99.26% 53.0 M 0.01 99.58% 30.1 M Comparison of sparsity levels and levels. memory usage at different threshold

99 12 t 0.5 0.4 02 50 o 40 30 10 X axis (mm) 0 12 14 0.8 0.4. 02. 20 10 X axis (mm) so 40Y axis (mm) Y axis (mm) Y axis (mm) (a) (b) 12 0.5 0.4m 02 0 O0 12' -0.4 02 50 Y axis (mm) so X axis (mm) Y axis (mm) X axis (mm) o U (c) (d) Figure 5.9: g Effect of thresholding on the current distribution after applying different threshold levels: (a) 0.001%, (b) 0.01%, (c) 0.05% and (d) 0.1%.

100 0 100 200 300 400 (a) (b) '100 100 200 200\\\\ 300 30 400 400 0 100 200 300 400 0 100 200 300 400 (c) (d) Figure 5.10: Structure of the moment matrix generated by the one level multiresolution expansion after thresholded by levels (a) 0.001%, (b) 0.01%, (c) 0.05% and (d) 0.1%.

101 0 0 -002 001 -04 -004 - Y(m) -005 X (m) a (m) 45 X (m).......................... (a) Scae=2.043639e-02 -2, -4 I o-2 -J r. 1 > " *, -10.-.0 0 Y(m) -0.05 X (m) (b) Figure 5.1: This figure shows the absolute values of the elements of the impedance matrix in Fig. 5.10 before thresholding v.s. the distance between the two basis functions. The figure is categorized by the types of basis functions shown at the right hand side of the figure.

102 Scale=6.948028e-05 -2, -4 0 -8.0 00 -0.04 -0.04 -0 Y (in) -0.05 X (in) Scale=1.568894e-02 0, 12 -0.03 -0.05 X (in) V (in) (d) Figure 5.1 1: Continued from the previous page.

103 5.7 Conclusion In this chapter, the multiresolution expansion was incorporated into the formulation of the method of moments. This formulation was applied to microstrip patches and arrays. The result indicates that by using wavelet bases, the impedance matrix is highly sparse when the system is large. Thus, a significant reduction in memory size is achieved.

104 -O 10 Threshold Figure 5.12: Figure 512Comparison of sparsity v.s. tiresolution expansions. threshold levels of one- and two-level mul

105 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1.. 01. 1.45 1.49 1.53 1.57 Frequency (GHz) 1.61 1.65 (a) 0.8 0.7 0.6 0.5 0.4 3 -=. -) 0.3 0.2 0.1 01.45 1.45 1.49 1.53 1.57 Frequency (GHz) 1.61 1.65 (b) Figure 5.13: Effect of thresholding on the resonant curves of (a) one-level expansion, (b) two-level expansion.

106 z /e d| Ground Plane Figure 5.14: A patch antenna array 1.... I......... 0.98 - 2-level expansion - - 1-level expansion IX 1.0.96.9, Q. 03 0.94 ) 0.92[ r ) 4...................... -5 10 10-2 Threshold Figure 5.15: g Comparison of the sparsity of one- and two-level expansions of the 3 x 3 array.

107 Figure 5.16: Structure of the sparse matrix of the 3 x 3 array at 0.01% (two-level expansion).

108 IE|l(db) at - = 0 90 90 IEel(db) at - = 90~ 90 90 - Array Factor -- 0.00001 - - 0.0001.. 0.0005 — 0.001 -- 0.01 g Effect; of thresholding on the un-normalized far-field pattern of the 3 x 3 array.

109 1.................. -...... I.-.. 0.99 0.98 Z%0.97 Cl).r4 cts Cu CD O.96 - 2-level expansion - - 1-level expansion --. --- - o -- e --- —-- 0.95 0.94 0. 3........ -.. 10-2 Threshold Fg Comparison of the sparsity of one- and two-level expansions of the 6 x 6 array.,. '..!...... ^............ '........-.........:::;:::::': ';: _:::::::::::::::::::::::::::...:.:.::.................................................................. -................... —.................-.-...::.::::::::::::::::.....::::.:.....;::::::::: Fiffure 5 19' Structure of the sparse matrix of the 6 x6 array at a threshold level of 0.01...... (two-level expansion)...................................................,. —........................... -.............................................-.-.............................................................................. E...E.....-,..-............................................ -.................,r - ~~~~~~~~~~~~~~~rp...... I............-,-.-, -::::::::::::::::::::.: rk;.:;::::................................... -.-.... rFigure 5.19:Stutroftesas marxoth 6x6ara atatrsodlvlf 0.015V, (two-level expansion.).

110 IEl(db) at > = 0 90 90 IEol(db) at = 900 0 90 90 - Array Factor -- 0.00001 - - 0.0001.. 0.0005 - 0.001 - - 0.01 Figure 5.20: The effect of thresholding on the normalized far-field pattern of the 6 x 6 array.

CHAPTER VI Application of Parallel Computing in MoM 6.1 Introduction Due to advance in semiconductor industry, the price of high performance processors falls rapidly, making it economical to build parallel computers with large number of processors. The need of high performance computing not achievable by single processor also drives the development of parallel computing technology. During the past decade we saw a lot of commercial parallel computers made by IBM, SGI and INTEL, to name a few. Parallel computers are no longer the proprieties of government or university research laboratories, but a reality. However, there is one cause hindering its widespread adoption. To tap into the power of parallel computing, the user has to modify his program according to the programming interface provided by the vender. Every time the user changes the parallel computer he used, he will need to modify his program again. There is no common programming interface among parallel computers until the Message Passing Interface (MPI) was standardized in 1994 [32]. Since then, MPI has been widely adopted by the computer industry. To utilize parallel computers, first a programmer needs to design a parallelization model for his program. The idea of parallelization is to divide a single task to many independent jobs. A simple and easy way to implement is the Single Program 111

112 Multiple Data (SPMD) model. In this model, different data sets are executed by the same program on different processors. These data sets have to be independent of each other, otherwise they cannot be executed in parallel. One simple example in electromagnetic problems is frequency response. In FEM or MoM, the computation process at each frequency point will not effect the outcome of other frequencies. Thus, we can simply submit the whole computation process of one frequency point to one processor [42, 43]. However, this may not achieve the best efficiency in terms of processor utilization. For instance, if there are 8 frequency points and 7 processors, when the last frequency point is computed on one processor, all other 6 processors will be idling. To achieve better efficiency, it is necessary to go to the inside of the computation process to seek out other independent operations which can be parallelized. Many numerical techniques for solving electromagnetic problems have utilized parallel computers. Especially, the FDTD has attracted a lot of effort in parallelization due to its computer-intensive nature. A comprehensive study of its implementation can be found in [55]. A detailed description of a parallelized implementation of a MoM code is presented in [56]. A comparison of linear equation solvers for integral-equation formulations can be found in [57]. In the following sections, two processes in MoM are parallelized using MPI. One is the impedance matrix fill process; the other is the FWT process. The detail of the implementation and the measurement of the turn around time are included. 6.2 Parallelization Model of MoM Matrix Fill Process The formulation of MoM always leads to a matrix equation as follows Ax = B, (6.1)

113 where A is a two-dimensional matrix, B is a column vector, and x is the unknown column vector we want to solve by direct inversion of A or by an iterative linear equation solver. The elements of matrix A are often the result of numerical integration of two basis functions with an integration kernel. The following equation is an example in one-dimensional space: aij = f (x) g(x, x') fj(x) dx dx' where aij is the element of the matrix at i-th row and j-th column, fi(x) and fj(x) are the basis functions, and g(x, x') is the kernel. Depending on the form and the distance between the basis functions, the time needed to perform the numerical integration of one element may vary greatly. To fully utilize the processors of a parallel computer, a dynamical job allocation scheme is necessary. As shown in Fig. 6.1, a processor is reserved as a dispatcher to send jobs and receive results from the workers. Other processors are served as workers to perform the numerical computation. 4 types of messages are passed between the dispatcher and workers. They are FREE: sent by a worker to inform the dispatcher that it is ready to accept new jobs, RESULT: sent by a worker to the dispatcher, containing the result of numerical integration, TERMINATION: sent by the dispatcher to inform the worker that all jobs are sent. JOB: send by the dispatcher, containing the information of a job. Following are the procedures performed by the dispatcher: 1. create a list of the unique elements, 2. wait for messages from the workers,

114 3. if the message is a RESULT message, store the result, otherwise the message is a FREE message and a JOB message is sent to the worker, 4. repeat step 2 and 3 until all jobs are sent and returned, 5. send TERMINATION message to inform each worker to quit. For the worker, the procedures are: 1. send FREE message to the dispatcher, 2. wait for messages from the dispatcher, 3. if it is a JOB message, perform the numerical integration and return the result, otherwise it is a TERMINATION message and the worker quits. 4. goto step 2. Fig. 6.2 shows the flow chart of the above process. The dashed lines in the figure indicate the communication paths. Two kinds of processes, dispatcher and worker, are concurrently executed. First, the process of the dispatcher sorts out the unique matrix elements and builds a job list after it starts and reads in the necessary parameters. Then, it enters a loop waiting for messages from the workers. The loop terminates only after all the jobs are sent and all the data are collected from the workers. Inside the loop, the dispatcher checks the messages from the workers. If the message contains data, the dispatcher will store it and then wait for new messages. If the message is a FREE message, the dispatcher checks if all the jobs have been sent. If not, the dispatcher will send a job to the worker which sends the message. If all jobs have been sent, the dispatcher will send a TERMINATION message to the worker. If all worker have been terminated, the dispatcher exits the loop,

115 Parallelization Model of Impedance Matrix Fill -r Command path --- Data path Worker Worker Worker i Dispatcher * Send commands to Workers. x * Receive data from Workers. * Store data in impedance matrix. * Receive commands from Dispatcher. * Perform computation. * Send result to Dispatcher. Figure 6.1: A dynamical allocation scheme for matrix fill process. outputs the results and quits. If there are still active workers, the dispatcher will continue waiting for new messages. The worker process, after starts and reads in the necessary parameters, sends a FREE message to the dispatcher to indicate that it is ready to accept new job. Then, it waits for a reply from the dispatcher. If the reply is a TERMINATION message, the worker quits. Otherwise, the message contains a job and the worker starts to execute it. After the worker has finished the job, the result is sent back to the dispatcher via a RESULT message. After sending the RESULT message, the worker sends a FREE message to the dispatcher and waits for the reply from the dispatcher again. A short program section is included in Appendix C to reveal the detail of implementation. A measurement of turn around time v.s. number of processors is shown in Fig. 6.3. Here the "turn around time" means the elapsed time from the start of the program

116 Dispatcher Worker Figure 6.2: The flow chart of a dynamical allocation scheme for matrix fill process.

117 until completion. The measurement is performed on a 48 processor IBM SP2 machine. In the figure, the theoretical result is curve-fitted to the experimental result according to the following formula time(i) = P + Comm x n + Comp/(i - 1), for i > 1 (6.2) where time(i): turn around time, P: preprocessing time, Comm communication time per message, Comp total computation time of the jobs, i number of processors, n number of transmitted messages. Note that Comm includes the communication latency and the actual transmission time of one message. Since the message size is fixed, Comm is a constant. Let N be the total number of the jobs. From the flow chart, there are N + i - 1 FREE,, N JOB, N RESULT, and i - 1 TERMINATION messages. Thus, the total number of transmitted messages is 3 x N + 2(i -1). When N is much larger than i, the total communication time is proportional to N. From (6.2), it is seen that when the number of processors increases, only the third term in the formula decreases. This means that the preprocessing time plus the total communication time governs the limit. The formula also indicates that the curve of the turn around time should approximate a hyperbola when N is much larger than i. The measured results show very good agreement with the above theoretical model. The improvement of turn around time is larger when the number of processors is small. However, when the number of processors is large, only a small improvement is achieved because the preprocessing and communication time dominate. Note that

118 the turn around time of 2 processors is higher than 1 processor. This is due to the higher preprocessing time of parallelized program and the communication cost. Also note that, due to the programming model, one processor is always reserved for coordinating other processors. When the communication cost is small, this processor will mostly be idle. It is possible to also utilize the idle time on this processor to perform numerical integration by more sophisticated programming. 6.3 Parallelization Model of FWT As discussed in Section 5.4, the FWT process transforms the matrix derived in MoM from one set of basis to another. (5.14) and (5.15) are repeated below for convenience. i Z=E pPijqj (6.3) [Z'] [P]t[Z][P] (6.4) In (6.3), 4i and q$ are the members of the basis 1( and 4X' respectively, and N is the rank of the basis. In (6.4), [Z] and [Z'] are the MoM matrices evaluated in bases I and ~' respectively. The pij in (6.3) is the i-th row, j-th column element of matrix [P] in (6.4). Let z' be the m-th row, n-th column element of [Z'], from the above two equations, we can derive N N Z..j = - PmiZmnPnj, (6.5) m=1 n=1 which is the form of the single operation in the FWT. The process of FWT fits to the SPMD model which is characterized by repeating the same operation, for example, (6.5), on different data. Comparing to the previous matrix-fill process, each operation of FWT is shorter and fixed, however the total number of operations is much larger. If the same algorithm as the previous one is

119 used, the overhead of communication will be too great. In stead, a fixed job-allocation scheme as shown in Fig. 6.4 is devised. At the beginning, each processor computes the job list and equally divides it to the total number of processors. Each processor can determines its share of jobs by his "rank"1. The only communication that takes place is when each processor sends its final computed results to one designated processor for post-processing. Following are the procedures performed by each processor: 1. create a list of the operations, 2. divide the list equally into n shares, where n is the total number of processors, 3. use the rank of the processor as an index to determine which share belongs to itself, 4. perform the FWT, 5. send the computed result to the processor with rank equal to 0. Fig. 6.2 shows the flow chart of the above process. A short program section is included in Appendix C to reveal the details of the implementation. A measurement of turn around time v.s. number of processors is shown in Fig. 6.6. The definition of "turn around time" is the same as in the previous section. The measurement is also performed on the same parallel computer. In the figure, the theoretical result is curve-fitted to the experiment result according to the following formula time(i) = P + Tran x N + L x i + Comp x N/i, (6.6) 1A terminology used by MPI to number the processors. If there is n processors, the rank of the processors will be 0, 1, 2,..., n.

120 where time(i):turn around time, P:preprocessing time, Tran:actual transmission time of one unit of data, L communication latency, Comp computation time, i: number of processors, N:number of unique matrix elements. Since each processor only sends computed results once, the total latency time equals to L x i. Comparing to other terms in (6.6) the total latency is negligible since the latency of SP2 is only 40 iS. As in the previous section, when the number of processors increases, only the last term in the formula decreases, while the first and second terms remain the same. This means that the preprocessing time and the total communication time is the limit. The formula also indicates that the curve of the turn around time should be a hyperbola. The measured results show very good agreement with the above theoretical model. The improvement of turn around time is larger when the number of processors is small. However, when the number of processors is large, only diminishing returns are obtained because the preprocessing and communication time dominates. 6.4 Conclusion In this section, the matrix-fill process and FWT in MoM are parallelized. The measured improvement is consistent to the theoretical model of the turn around time. Another bottleneck not addressed in this chapter is the process of solving (6.1). When the size of matrix A is large, direct inversion of A becomes impossible.

121 Thus, iterative solvers need to be used. An iterative solver by nature is sequential. Each iteration relies on the result of last iteration. In this dissertation, a bi-conjugate gradient [6] iterative solver is used when the matrix size is large. In this solver, the most time-consuming part is the multiplication of a matrix with a column vector. Every iteration produces an improved guess of the solution vector. Thus, the communication overhead of at least one solution vector in every iteration is not avoidable. When the dimension of the matrix is large, this communication overhead outweighs the improvement incurred by parallelizing the linear solver, since the communication speed is much slower than computation speed. This problem still needs to be addressed if further improvement on the speed of MoM is necessary.

122 Impedance matrix fill, 240 matrix elements 1; 0 Q En *4 -0 E I 5 10 15 20 25 30 Number of CPUs g The measured turn around time of the impedance matrix fill process. The theoretical result is plotted by curve-fitting the measured data to (6.2)

123 Parallelizing Model of Fast Wavelet Transform o Data path ** Perform FWT. * Evenly divided by N. * Send each share to each CPU. * Gather results from each CPU. Figure 6.4: A static allocation scheme for the FWT process.

124 Worker 0 Worker 1....... Worker n Figure 6.5: The flow chart of a static allocation scheme for the FWT process.

125 Fast wavelet transform, 3568 matrix elements 0 I/q) $=! 0 H F: 2 4 6 8 10 12 14 16 Number of CPUs Figure 6.6: The measured turn around time of the FWT process. The theoretical result is plotted by curve-fitting the measured data to (6.6)

CHAPTER VII Conclusion 7.1 Summary of Achievements The goal of this dissertation is to develop efficient numerical techniques which can be used in modeling complex and large structures like MMIC's. To achieve this goal, first we have developed a hybrid MoM/FEM technique which combines the efficiency of MoM and the flexibility of FEM. As shown in Chapter 2-4, this technique has been successfully applied to the design of a microstripline-fed slot-coupled cavity-backed patch antenna, a coaxial line-fed cavity-backed patch antenna, a wide bandwidth low-insertion-loss microstripline vertical cavity coupler, and a high-Q micromachined cavity resonator. Second, we have incorporated multiresolution expansions to the method of moments to utilize the locally oscillatory and zero averaging property of wavelet bases. These properties tend to make the off-diagonal terms of the impedance matrix very small when these terms involve the integration with a wavelet. These terms can be neglected and occupy no memory space in a sparse matrix storage scheme. Thus, the memory size required to store the matrix is reduced. Our test on patch antennas shows that when the problem size is large, significant saving in memory is possible. The 6 x 6 patch array in Chapter V reduces memory size to only 1.16% of the original 126

127 under a low threshold value of 0.01%. Third, in order to utilize parallel computers, the programs for the impedance matrix fill process and the FWT process are ported to MPI. The programs have been tested on a 48 CPU IBM SP2 parallel computers. Significant improvement in speed is evidenced by the measured turn around time. 7.2 Future Work Following is the research work that can be continued in the future. First, the high-Q resonator introduced in Chapter IV has been fabricated and tested. The feasibility of the design concept has been proved. To fulfill its potential, the next step is to design a band-pass filter by stacking several cavities vertically. Second, the bottleneck in MoM or FEM is always in the matrix fill process and the iteration solver. Since the matrix fill process can be parallelized, the iteration solver remains as the primary bottleneck. However, the parallelization of the iteration solver is difficult due to the large quantity of data passed around in each iteration. The benefit of parallel processing is always outweighed by the extra communication cost. Thus, it is necessary to develop a more efficient algorithm to parallelize the iteration solver. Finally, in the multiresolution expansion, although we get good result out of thresholding the impedance matrix, how the error relates to the threshold level is not clear. Hence, a systematic study on the relationship between the error and threshold level is necessary to ensure the accuracy of the results.

APPENDICES 128

129 APPENDIX A Tetrahedral Edge Elements Fig. A.1 shows the labeling system of the nodes and edges of a tetrahedral element. Let We be the edge basis corresponding to edge Te. Referring to Table A.1, We for i = 1,2,...,6 is defined by we f7-i + 97-i x r r E Ve ^7-i 0 r Ve -- b7-i - f7-i = ril x ri2, 6Ve bb7-ii_ - 7-i ~ bT 6Ve ri2 - ril bi bi = Iri2 - rll Ve = element volume. Note that We satisfies following relationships V W =0, V x W 27i, where and is avector on j-th edge of the tetrahedron. where 8ij = < and r3 is a vector on j'-th edge of the tetrahedron. [0 i j

130 ni n4 e5 n2 n3 Figure A. 1: The node and edge labels of a tetrahedral element. ni, n2, n3 and n4 are the nodes. el, e2,..., e6 are the edges. Node number Edge Number ii i2 1 1 2 2 1 3 3 1 4 4 2 3 5 3 4 6 4 2 Table A.I: Labeling system of a tetrahedron.

131 APPENDIX B Green's Functions of Conductor-backed Dielectric Slabs Fig. B.1 shows a conductor-backed dielectric slab with thickness d and dielectric constant cr. The notation of the Green's functions are explained below dFSd ( kx, ky; Z2, Z1), dl: direction of the source. Can be x or y, d2: direction of the field. Can be x or y, F: field type. E for electric field; H for magnetic field, S: source type. J for electric current; M for magnetic current, zl: z component of the source point coordinate, z2: z component of the field point coordinate. The following definitions are used in the Green's functions. ko = w vU/o zo - V — CO d2 = k2 + k2

132 z * R y Dielectric I:.:.....i:": -~i i... 1....-...''~i i 1i''iiil1i~i~ I -ii I~ Jll~i.. - I~ i'- i........... Ground Plane Figure B.1: The structure of a conductor-backed dielectric slab. kl = rko2 /2, Q(k1) < 0, }(k1) > 0, k2 = ko -, (k2) < 0, ) (k2) > 0. Te - k1 cos k1d +j k2 sin k1d Tm = 6rk2 cos k1d + jkl sin k1d The spectral domain Green's functions used in this dissertation are listed below: G, ( k k. Zo (erk2 - ky)k2 cos(kid) + j(k o2- ky2)k sin(kid) J ykTeTm Gy f k d, ) -= Zo jkky(k2 cos(ki d) + j k sin(kid)) EJ Y ) 47r2ko TeTM G k k k ) -= 1 -eklk2cos(kid) +j(k2(er - 1) - k 2)sin(kld) GHY3(kx k ky; O, d) = H 47r2 TeTm GHM(kx k; 0,0) = 47rkoZo k T [j (kr(e - k) + ( - x {kik2(er + 1) sin(kid) cos(kid) + j(erk sin2(kd) - k2 cos2(kld))}], GX k k O ) - 1 jkxky(er- 1) sin(kid) 4HY ) 2 TeTm Note that due to reciprocity, Gyx(k ky; dd) = G (k ky; d), GEM(kx, ky; d, 0) = -GHJx(kx ky; 0 d), GxM(kx, ky; d, 0) = -G(j(kx, ky; 0, d). Also GOxx(k, kg,; d, d) can be derived from GYYJ(kx, ky; d, d) by exchanging x and y.

133 APPENDIX C The Implementation of Parallel Programming in C++ The following two sections show the implementation of the parallel programming of dynamic and static job allocation schemes in C++ language. The definition of the MPI functions can be found in [32, 33]. C.1 Dynamic job allocation int rank, processsize; mpidata mdata; /* common buffer for MPI messages */ MPI.Comm.rank(MPICOMMWORLD, &rank); if (rank==O) { /* Dispatcher Process */ int send.count=O, receivecount=O, actual_count; int *process=new intprocess-size-1]; /* for storing the status of processors */ for (int i=O;i<process.size-1;i++) process[i]=0;

134 intpair* ip;/* pointer to the job list */ int finished=O; do MPIStatus status; /* receive message from workers */ MPI_Recv(&mdata, sizeof(mdata), MPIBYTE, MPIANY.SOURCE, MPI_ANYTAG, MPICOMMWORLD, &status); int sent=0; switch (status.tag) /* check message type */ case TAGREQUEST: /* request for job */ if (finished) /* inform the work to end */ MPI_Send(&mdata, sizeof(mdata), MPI.BYTE, status.source, TAG.END, MPI-COMMWORLD); process [status.source-1]=1; sent=l; }; while (!finished &&!sent) /* send job to the worker */ mdata.j ob=ip[send-count] mdata.i=send_count++; MPI_Send(&mdata, sizeof(mdata), MPI.BYTE, status.source, TAGDATA, MPI-COMMWORLD);

135 actual_count++; if (sendcount>=count) finished=l; }; sent=l; break; case TAGDATA: /* returned results */ ip[mdata.i].v=mdata.v; receivecount++; break; } while (receivecount<actual_count I finished==O); /* inform workers to end */ for (int i=l;i<process.size;i++) { MPI.Status status; if (process[i-1]==0) MPIRecv(&mdata, sizeof(mdata), MPI-BYTE, i, TAGREQUEST, MPI.COMMWORLD,&status); MPI_Send(&mdata, sizeof(mdata), MPI-BYTE, i, TAG.END, MPI-COMMWORLD);

136 }; /* output the results */ } else /* Work Process */ MPI_Status status; while ( MPI_Send(&mdata, sizeof(mdata), MPIBYTE, 0, TAGREQUEST, MPICOMMWORLD), MPI_Recv(&mdata, sizeof(mdata), MPI_BYTE, 0, MPIANYTAG, MPI_COMMWORLD, &status), status.tag==TAG-DATA) /* compute the job stored in mdata.job */ /* store the result in mdata.v*/ /* return the result in mdata.v */ MPI_Send(&mdata, sizeof(mdata), MPIBYTE, 0, TAGDATA, MPI_COMM_WORLD); }; }; MPI_Finalize(); C.2 Static Job Allocation int size; /* total number of jobs */

137 int rank, processsize; mpidata mdata; /* common buffer for MPI messages */ MPICommrank(MPI-COMM-WORLD, &rank); MPICommsize(MPI-COMMWORLD, &process.size); /* compute the size of jobs for each processors */ int count=size/processsize; int remain=size % process.size; /* remaining jobs */ int *recvcount, *displs, offset; if (rank==O) /* processor 0 is also in charge of collecting the results, so it needs to prepare the place to store them */ recvcount=new int[process_size]; displs=new int[processsize]; for (int i=O; i<processsize; i++) if (i<remain) displs[i]=i*(count+l)*sizeof(intpair); recvcount[i]=(count+l)*sizeof(intpair); else {

138 displs[il=(i*count+remain)*sizeof(intpair); }; }; }; /* adjustment for the remaining jobs */ if (rank<remain) count++; offset=count*rank; else else { offset=count*rank+remain; /* compute the jobs */ /* return the results */ MPI_Gatherv(&(z.ia[offset]),count*sizeof(intpair),MPIBYTE,&(z.ia[O]), recvcount, displs, MPIBYTE, 0, MPI.COMMWORLD); if (rank==O) { /* processor 0 output the total results */

BIBLIOGRAPHY 139

140 BIBLIOGRAPHY [1] J.A. Navarro, K. Chang, J. Tolleson, S. Sanzgiri and R. Q. Lee, "A 29.3 GHz cavity-enclosed aperture-coupled circular patch antenna for microwave circuit integration," IEEE Microwave and Guided Wave Letters, vol. 1, no. 7, pp. 170 -171, July 1991. [2] B.A. Brynjarsson and T. Syversen, "Cavity-backed, aperture coupled microstrip patch antenna," IEE 1993 International Conference on Antennas and Propagation Symposium Proceedings, pp. 715-718. [3] F. Zavosh and J.T. Aberle, "Infinite phased arrays of cavity-backed patches," IEEE Trans. on Antennas Propagat., vol. AP-42, no. 3, pp. 390-398, Mar. 1994. [4] J.T. Aberle and F. Zavosh, "Analysis of probe-fed circular microstrip patches backed by circular cavities," Electromagnetics, vol. 14, no. 2, pp. 239-258, Apr. 1994. [5] J. Lee, T. Horng and N.G. Alexopoulos, "Analysis of cavity-backed aperture antennas with a dielectric overlay," IEEE Trans. on Antennas Propagat., vol. AP-42, no. 11, pp. 1556-1562, Nov. 1994. [6] J. Jin, The Finite Element Method in Electromagnetics. New York: John Wiley and Son, 1993. [7] J. Jin and J.L. Volakis, "A hybrid finite element method for scattering and radiation by microstrip patch antennas and arrays residing in a cavity," IEEE Trans. on Antennas Propagat., vol. AP-39. no. 11, pp. 1598-1604, Nov. 1991. [8] J. Jin and J.L. Volakis, "Scattering and radiation analysis of three-dimensional cavity arrays via a hybrid finite-element method," IEEE Trans. on Antennas Propagat., vol. AP-41, no. 11, pp. 1580-1586, Nov. 1993. [9] J. Gong, J.L. Volakis, A.C. Woo, and H.T.G. Wang, "A hybrid finite elementboundary integral method for the analysis of cavity-backed antennas of arbitrary shape," IEEE Trans. on Antennas Propagat., vol, AP-42, no. 9, pp. 1233-1242, Sept. 1994. [10] X. Yuan, D.R. Lynch and J.W. Strohbehn, "Coupling of finite element and moment methods for electromagnetic scattering from inhomogeneous objects," IEEE Trans. Antennas Propagat., vol. AP-38, no. 3, pp. 386-393, Mar. 1990.

141 [11] W.E. Boyse and A.A. Seidl, "A hybrid finite element method for near bodies of revolution," IEEE Trans. Magn., vol. 27, no. 5, pp. 3833-3836, Sept. 1991. [12] S.D. Gedney and R. Mittra, "Analysis of the electromagnetic scattering by thick gratings using a combined FEM/MM solution," IEEE Trans. Antennas Propagat., vol. AP-39, no. 11, pp. 1605-1614, Nov. 1991. [13] T. Cwik, "Coupling finite element and integral equation solutions using decoupled boundary meshes," IEEE Trans. Antennas Propagat., vol. AP-40, no. 12, pp. 1496-1504, Dec. 1992. [14] D.J. Hoppe, L.W. Epp, and J. Lee, "A hybrid symmetric FEM/MOM formulation applied to scattering by inhomogeneous bodies of revolution," IEEE Trans. Antennas Propagat., vol. AP-42, no. 6, pp. 798-905, June 1994. [15] Dj. Jankovic, M. LaBelle, D.C. Chang, J.M. Dunn, and R. C. Booton, "A hybrid method for the solution of scattering from inhomogeneous dielectric cylinders of arbitrary shape," IEEE Trans. Antennas Propagat., vol. AP-42, no. 9, pp. 1215-1222, Sept. 1994. [16] W.E. Boyse and A.A. Seidl, "A hybrid finite element method for 3-D scattering using nodal and edge elements," IEEE Trans. Antennas Propagat., vol. AP-42, no. 10, pp. 1436-1442, Oct. 1994. [17] M.L. Barton and Z.M. Cendes, "New vector finite elements for three-dimensional magnetic field computation,' J. Appl. Phys., vol. 61, no. 8, pp. 3919-3921, April 1987. [18] A. Chatterjee, J. Jin and J.L. Volakis, "Edge-based finite elements and vector ABC's applied to 3-D scattering," IEEE Trans. Antennas Propagat., vol. AP-41, no. 2, pp. 221-226, Feb. 1993. [19] R.F. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGrawHill, 1961. [20] D.M. Pozar, "Microstrip antenna aperture-coupled to a microstripline," Electron. Lett., vol. 21, no. 2, pp. 49-50, Jan. 1985. [21] P.L. Sullivan and D.H. Schaubert, "Analysis of an aperture coupled microstrip antenna," IEEE Trans. Antennas Propagat., vol. AP-34, no. 8, pp. 977-984, Aug. 1986. [22] R.N. Simons and R.Q. Lee, "Coplanar waveguide aperture coupled patch antennas with ground plane/substrate of finite extent," Electron. Lett., vol. 28, no. 1, pp. 75-76, Jan. 1992. [23] L.P.B. Katehi, N.I. Dib and J.Cheng, "Slot coupled patch antennas," Technical Report 390442-2T, Univ. Michigan, Ann Arbor, MI., Dec. 1993.

142 [24] D.M. Pozar, "Input impedance and mutual coupling of rectangular microstrip antennas,"' IEEE Trans. Antennas Propagat., vol. AP-30, pp. 1191-1196, Nov. 1982. [25] N.K. Das and D.M. Pozar, "Multiport scattering analysis of general multilayered printed antennas fed by multiple feed ports: part I -theory," IEEE Trans. Antennas Propagat., vol. AP-40, pp. 469-481, May 1992. [26] N.K. Das and D.M. Pozar, "Multiport scattering analysis of general multilayered printed antennas fed by multiple feed ports: part II - applications," IEEE Trans. Antennas Propagat., vol. AP-40, pp. 482-490, May 1992. [27] R.W. Jackson and D.M. Pozar, "Full-wave analysis of microstrip open-end and gap discontinuities," IEEE Trans. Microwave Theory Tech., vol. MTT-33, pp. 1036-1042, Oct. 1985. [28] K. Kunz and R. Luebbers, The Finite Difference Time Domain Method for Electromagnetics, Florida: CRC press, 1993. [29] X. Zhang and K. Mei, "Time-domain finite difference approach to the calculation of the frequency-dependent characteristics of microstrip discontinuities, " IEEE Trans. Microwave Theory and Techniques, pp. 1775-1787, Dec. 1988. [30] D. Sheen, S. Ali, M. Abouzahra and J. Kong, "Application of the ThreeDimensional Finite-Difference Time-Domain Method to the Analysis of Planar Microstrip Circuits, " IEEE Trans. Microwave Theory and Techniques, pp. 849 -857, July 1990. [31] S. Visan, O. Picon and V. Hanna, "3D Characterization of Air Bridges and Via Holes in Conductor-Backed Coplanar Waveguides for MMIC Applications, " 1993 IEEE MTT-S Intl. Microwave Symp. Dig., pp. 709-712. [32] Message Passing Interface Forum, MPI: A Message-Passing Interface Standard, International Journal of Supercomputer Applications, vol. 8, nos. 3/4, 1994. [33] Peter S. Pacheco, A User's Guide to MPI, Department of Mathematics, University of San Francisco, CA, 1995. [34] M. Davidovitz, R.A. Sainati, S.J. Fraasch, "A Non-Contact Interconnection Through an Electrically Thick Ground Plate Common to Two Microstrip Lines, " IEEE Trans. Microwave Theory and Techniques, pp. 753-759, April 1995. [35] R.F. Drayton and L.P.B. Katehi, "Microwave Characterization of Microshield Lines," Digest of the 40th ARFTG Conference, Orlando, FL, Dec. 1992. [36] R.F. Drayton and L.P.B. Katehi, "Experimental study of micromachined circuits," Digest of the 1993 International Symposium on Space Terahertz Technology, Los Angeles, CA, March 1993.

143 [37] R.F. Drayton and L.P.B. Katehi, "Micromachined circuits for MM-wave applications," Digest of the 1993 European Microwave Conference, Madrid, Spain, Sept. 1993. [38] R.F. Drayton, T.M. Weller, and L.P.B. Katehi, "Development of miniaturized circuits for high-frequency applications using micromachining techniques," International Journal of Microcircuits and Electronic Packaging, third quarter, 1995. [39] C.Y. Chi, "Planar microwave and millimeter-wave components using micromachining technologies," Ph. D. dissertation, Univ. of Michigan Ann Arbor, 1995. [40] J. Cheng, N.I. Dib and L.P.B. Katehi, "Theoretical modeling of cavity-backed patch antennas using a hybrid technique," IEEE Trans. on Antennas Propagat., vol. AP-43, no. 9, pp. 1003-1013, Sept. 1995. [41] J. Papapolymerou, J. Ceng, J. East, and L.P.B. Katehi, "A micromachined high-Q X-Band resonator," IEEE Microwave and Guided Wave Letters, vol. 7, no. 6, pp. 168-170, June 1997. [42] J.-G. Yook and L.P.B. Katehi, "Characterization of MIMIC packages using parallelized 3D FEM code," Proc. of Advance in Computational Electromagnetics Symposium, pp. 954-961, 1996. [43] Jong-Gwan Yook, "Electromagnetic modeling of high-speed high-frequency interconnects," Ph. D. dissertation, Univ. of Michigan, Ann Arbor, 1996. [44] B.Z. Steinberg and Y. Leviatan,"On the use of wavelet expansions in the method of moments,"IEEE Trans. Antennas Propagat., vol. 41, pp. 610-619, May 1993. [45] K.F. Sabet and L.P.B. Katehi,"Analysis of integrated millimeter-wave and submillimeter-wave waveguides using orthonormal wavelets expansions," IEEE Trans. Microwave Theory Tech., vol. MTT-42, pp. 2412-2422, Dec. 1994. [46] G. Wang and G-W Pan,"Full wave analysis of microstrip floating line structures by wavelet expansion method," IEEE Trans. Microwave Theory Tech., vol. MTT-43, pp. 131-142, Jan. 1995. [47] G. Wang, G.-W. Pan and B.K. Gilbert, "A hybrid wavelet expansion and boundary element analysis for multiconductor transmission lines in multilayered dielectric media," IEEE Trans. Microwave Theory Tech., vol. MTT-43, pp. 664-675, March 1995. [48] K.F. Sabet and P.B. Katehi,"A study of dielectric resonators based on twodimensional fast wavelet algorithm," IEEE Microwave Guided Wave Lett., vol. 6, pp. 19-21, Jan. 1996. [49] I. Daubechies, Ten Lectures on Wavelets, Philadelphia, PA: Society for Industrial and Applied Mathematics, 1992.

144 [50] C.K. Chui, ed., Wavelets: A Tutorial in Theory and Applications. New York: Academic Press, 1992. [51] N.L. VandenBerg and L.P.B. Katehi, "Broadband Vertical Interconnects Using Slot-Coupled Shielded Microstrip Lines," IEEE Trans. Microwave Theory Tech., vol. MTT-40, no. 1, pp. 81-88, Jan. 1992. [52] E.J. Stollnitz, T.D. DeRose, and D.H. Salesin, "Wavelets for Computer Graphics: A Primer, Part 1 and 2," IEEE Trans. Computer Graphics and Applications, pp. 74-82, July 1995. [53] K.F. Sabet, L.P.B. Katehi, K. Sarabandi, "Wavelet based CAD modeling of microstrip discontinuities using least square Prony's method," 1997 IEEE MTTS Intl. Microwave Symp. Dig., pp. 1799-1801. [54] J. Jin, J.L. Volakis, and V.V. Liepa, "Fictitious absorber for truncating finite element meshes in scattering," Proc. IEE, Pt. H, vol. 139, pp. 472-476, Oct. 1992. [55] Z.M. Liu, A.S. Mohan, T.A. Aubrey, and W.R. Belcher, "Techniques for implementation of the FDTD method on a CM-5 parallel computer," IEEE Antenna and Propagation Magazine, vol. 37, No. 5, pp. 64-71, Oct. 1995. [56] J.M. Putnam, D.D. Car and J.D. Kotulskii, "Parallelization of the CARLOS-3D Method of Moments Code," Proceedings of the 11th annual review of progress in applied computational electromagnetics, Monterey, CA, pp. 848-863, March 1995. [57] K. Forsman, W. Gropp, L. Kettunen, D. Levine, and J. Salonen, "Solution of dense systems of linear equations arising from integral-equation formulations," IEEE Antenna and Propagation Magazine, vol. 37, No. 6, pp. 96-99, Dec. 1995.