Robust Development of Hybrid Finite Element Methods For Antennas and Microwave Circuits by Jian Gong A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1996 Doctoral Committee: Professor John L. Volakis, Chairperson Professor Linda P.B. Katehi Professor Thomas B.A. Senior Professor Kamal Sarabandi Professor Gregory M. Hulbert RL-942 = RL-942

~ Jian Gong 1996 All Rights Reserved

T o 9NT ~3arents tmg memortj an?. mgj encouragement To ODongli canZ' Q~endcen MOj inspircation ctnzMO m u~f iUlmcnt 11

ACKNOWLEDGEMENTS Many have helped me during my academic life and along my Ph.D. journey and they definitely deserve my deep appreciations. First of all, I am grateful to Prof. John L. Volakis who has been my advisor since 1991 when he 'grabbed' me from RIT in Rochester, New York. I would not have stayed here without his persuasive effort to the admission administration, where tedious and ineffective processes were usually required. His advisory, assistance, encouragement and continuous support made the dissertation research possible. My sincere gratitude also goes to my committee members, Profs. Thomas Senior, Linda Katehi, Kamal Sarabandi and Gregory Hulbert, for their evaluations and time spent to review the manuscript to improve the quality of the dissertation. Having stayed in the Radiation Laboratory for almost five years, I am full of pleasant memories of various events and activities, academically and socially. Many colleagues of mineduring this period, who may have already left or are still hanging around, have helped me here and there in variety of means. We all enjoy discussions of diverse and mutually interesting (cultural, political and technical) issues and we've learned from each other a great deal in many aspects (at least from my point of view). This atmosphere and environment is what I always prefer and what I will miss in the future. Among many of them in the Radiation Laboratory, the Department of EECS, or even the University, are (I chose to only list the most current JLV's group in alphabetic order) Hristos Anastassiu, Lars Andersen, Brown Arik, Yunus Erdemli, iii

Sunil Bindigavanale, Youssry Botros, Mark Casciato, Stephane Legault,1 Zhifang Li, Mike Nurnberger and Tayfun Ozdemir. I'd like to mention another fellow, Prof. T.S.M. Maclean, who was with the University of Birmingham, England and is now retired. During 1988-1990, I had the chance to work with him to investigate certain interesting radiowave propagation and scattering problems. What I was impressed is his excellent scientific research attitude and his personal model. His trust and instrumental recommendations on many occasions provided me with the important opportunities. Without his effort, I suspect that I would not have achieved what I have today. I am deeply indebted to my family. My father died when I was little during the Great Culture Revolution (GCR) in China. (He was physically tortured and beaten to death within five days by the Red Guards.) The custody burden of the children was then entirely loaded to my mother's shoulders. In spite of the incredible amount of spiritual and physical pressures, she took good care of her children and paid particular attention to their education. It was her effort and high expectation of us that encouraged me to pursue excellence in my school and college study and later in my academic career. Unfortunately she passed away at her early fifties because of the long term emotional depression. My sisters and I experienced this deplorable incident occurring in our family, and this incident, however, established an unusually close relation among us. I am so proud of all my five elder sisters! Last, but not the least, my wife Dongli deserves my sincere gratitude and appreciation. Our lovely son Davin (Chenchen) was two years old when this dissertation started being prepared. He often shows his clear preference between his mother and me because of the much less time I spent with him. My in-law parents have helped 1Officially Stephane does not belong to the JLV's group. He is however always included as an active member of the group. iv

us a great deal, especially in the later stage of this draft preparation. For years and years their patience, understanding, support and love have motivated and inspired me in my research work. Without them, this dissertation would not be so successful. v

TABLE OF CONTENTS DEDICATION................................... ii ACKNOWLEDGEMENTS.......................... iii LIST OF FIGURES............................... ix LIST OF TABLES................................ xv LIST OF APPENDICES............................ xvi CHAPTER I. Introduction.............................. 1 1.1 Overview....................................... 1 1.2 Fundamentals of Electromagnetic Theory........... 5 1.2.1 M axwell Equations.................. 6 1.2.2 Boundary Conditions and Boundary Value Problems 7 1.2.3 Uniqueness Theorem and Equivalence Principle... 10 1.2.4 Integral Equation and Dyadic Green's Function.. 13 II. Finite Element Analysis in Electromagnetics.......... 20 2.1 Functional Formulation.......................... 21 2.1.1 Pertinent Functional For Lossy/Anisotropic Media I.......................................... 25 2.1.2 Pertinent Functional For Lossy/Anisotropic Media -II................................... 27 2.2 Galerkin Formulation.............................. 28 2.3 Total Field and Scattered Field Formulations.............. 30 2.3.1 Scattered/Incident Fields and Boundary Conditions 31 2.3.2 Galerkin's Method.......................... 34 2.3.3 Variational Method....................... 336 2.4 Parameter Extraction............................. 38 2.4.1 Radiation and RCS Pattern..................... 38 2.4.2 Gain and Axial Ratio...................... 4(0 vi

III. Edge-Based FE-BI Technique................ 43 3.1 Introduction........................... 43 3.2 Hybrid System Functional.................... 44 3.2.1 FEM Subsystem........................... 45 3.2.2 Boundary Integral Subsystem............ 48 3.2.3 Combined FE-BI System.................... 50 3.3 Numerical Implementation...................... 51 3.4 Selected Numerical Results......................... 54 IV. Efficient Boundary Integral Subsystem - I............ 61 4.1 Introduction........................ 61 4.2 Application of Conjugate Gradient Algorithms........ 62 4.2.1 BiCG Algorithm With Preconditioning....... 62 4.2.2 BiCG-FFT Algorithm For Linear System....... 64 4.2.3 Convolutional Form of Boundary Integral......... 65 4.3 Mesh Overlay Scheme.............................. 71 4.3.1 Field Transformations................. 71 4.4 Results............................... 76 V. Efficient Finite Element Subsystem -II............... 80 5.1 Introduction............................ 80 5.2 Hybrid FE-BI Formulation....................... 81 5.3 Edge-Based Prismatic Elements....................... 85 5.4 Applications........................... 87 5.5 Concluding Remarks............................... 95 VI. Antenna Feed Modeling............................ 98 6.1 Probe Feed...................................... 98 6.1.1 Simple Probe Feed......................... 98 6.1.2 Voltage Gap Feed.......................... 99 6.2 Aperture-coupled Microstrip Model.................... 99 6.3 Coax Cable Feed................................ 102 6.3.1 Motivation....................... 102 6.3.2 Hybrid FE-BI System................ 103 6.3.3 Proposed Coax Feed Model.............. 103 6.3.4 Results and Conclusion................ 108 6.4 Conclusion............................ 109 VII. Circuit Modeling................................... 115 vii

7.1 Introduction.......................... 115 7.2 Numerical De-embedding....................... 117 7.3 Truncation Using DMT..................... 120 7.4 Truncation Using PML........................... 121 7.4.1 Theory......................... 122 7.4.2 Results............................ 124 VIII. AWE: Asymptotic Waveform Evaluation............. 133 8.1 Brief Overview of AWE..................... 133 8.2 Theory.............................. 134 8.2.1 FEM System Recast........................... 134 8.2.2 Asymptotic Waveform Evaluation.......... 136 8.3 Numerical Implementation......................... 138 IX. Conclusions................................. 142 9.1 Discussion on the Research Work................ 142 9.2 Suggestions for Future Tasks.................. 145 9.3 Modular Development.................................. 145 APPENDICES................................. 147 BIBLIOGRAPHY............................ 159 ill

LIST OF FIGURES Figure 1.1 (a) Recessed cavity in a PEC ground plane. (b) Recessed cavity in a dielectrically coated PEC ground plane........................ 12 1.2 Illustrations of equivalence principle when applied to the structure shown in fig. 1.la............................. 14 1.3 Illustrations of equivalence principle when applied to the structure shown in fig. 1. lb............................ 15 1.4 Examples of protruding configurations (a) on a planar platform (b) on a curved platform in consideration of the equivalence principle.. 16 2.1 Illustration of a typical conformal antenna configuration....... 22 2.2 Illustration of a scattering problem setup for scattered field formulation.................................. 30 3.1 Illustration of a typical radiation and scattering problem...... 45 3.2 A tetrahedron and its local node/edge numbering scheme...... 46 3.3 Pair of triangles sharing the ith edge................. 49 3.4 A typical geometry/mesh for a cavity-backed circular patch antenna. 51 3.5 A flow chart describes the major implementation procedures from the mesh generation, a few data preprocessors, the FE-BI kernel, to the BiCG solution and finally the data output............. 53 3.6 Comparison of the computed and measured oa-0 backscatter RCS as a function of frequency for the shown circular patch. The incidence angle was 30~ off the ground plane............................ 56 ix

3.7 Comparison of the computed and measured input impedance for the circular patch shown in fig. 3.6. The feed was placed 0.8 cm from the center of the patch and the frequency was swept from 3 to 3.8 GHz. 56 3.8 Illustration of the configuration and mesh of the one-arm conical spiral used for the computation of fig. 3.9................ 58 3.9 Comparison of the calculated radiation pattern (E4), taken in the - = 90~ —plane, with data in reference for the one-arm conical spiral shown in fig. 3.8............................. 58 3.10 Comparison of input impedance calculations for the illustrated cavitybacked slot................................ 59 3.11 Visualization of the near field distribution at the lower layer of a stacked circular patch antenna..................... 60 4.1 Structured mesh consists of equal right triangles........... 66 4.2 Illustration of two triangles with the corresponding indices to help to prove the convolutional property of the boundary integral..... 70 4.3 Printed circular patch antenna is modeled using the recessed scheme to incorporate the BiCG-FFT algorithm. (a) Illustration of the configuration; (b) Comparisons of the BiCG-FFT result with that of the ordinary FE-BI technique presented in chapter 3....... 72 4.4 Overlay of a structured triangular aperture mesh over an unstructured mesh, shown here to conform to a circular patch....... 73 4.5 Illustration of the parameters and geometry used in constructing the transformation matrix elements between the structured and unstructured mesh................................. 74 4.6 Illustration of the cavity-backed triangular patch array.........77 4.7 Comparisons of the monostatic radar cross section scattering by a 2 x 2 triangular patch array shown in fig. 4.6. The reults were computed using the regular BiCG FE-BI technique described in chapter 3 and using the BiCG-FFT proposed in this chapter............. 78 x

4.8 Comparisons of the monostatic radar cross section scattering by an empty aperture with the same cavity size and dielectric filling (6,=1) as the structure shown in fig. 4.6. Again, the results were calculated using the regular BiCG FE-BI technique described in chapter 3 and using the BiCG-FFT proposed in this chapter............ 78 4.9 Bistatic RCS scattering by a crcular patch antenna modeled using the regular BiCG FE-BI and the BiCG-FFT algorithm with overlaying transform............................. 79 5.1 Geometry of cavity-backed microstrip antennas........... 81 5.2 Illustration of tessellation using prisms................ 84 5.3 Right angled prism........................... 84 5.4 Geometry of the annular slot antenna backed by a cavity 23.7 cm in diameter and 3 cm deep........................ 88 5.5 Scattering: Bistatic (co-pol) RCS patterns computed using the tetrahedral FE-BI code and the prismatic FE-BI code. The normally incident plane wave is polarized along the b = 0 plane and the observation cut is perpendicular to that plane. Radiation: X-pol and Co-pol radiation patterns in the = 0 plane from the annular slot antenna shown in fig. 5.4. The solid lines are computed using the tetrahedral FE-BI code whereas the dotted lines are computed using the prismatic FE-BI code. The excitation probe is placed at the point (y==0) marked in fig. 5.4.................... 89 5.6 Illustration of the setup for computing the FSS transmission coefficient Upper figure: periodic element (top view); Lower figure: periodic element in cavity (cross-sectional view).......... 90 5.7 Calculations and comparisons of transmission through the FSS structure shown in fig. 5.6........................ 92 5.8 Upper figure: geometry of the multilayer frequency selective surface (FSS) used for modeling; lower figure: measured and calculated transmission coefficient through the FSS structure.......... 93 5.9 Illustration of a typical 2-arm slot-spiral design........... 96 5.10 Radiation Pattern at f=1.1GHz (center frequency design). A good axial ratio is achieved up to 60~ degree................ 96 xi

5.11 Radiation Pattern at f=0.944GHz (lower end of frequency range). It can be seen that the axial ratio of the pattern becomes larger compared to that at the center frequency, but still remains within 3dB for a wide angle range. This indicates that the number of the outer turns in the spiral contour design is most likely sufficient.... 97 5.12 Radiation Pattern at f=1.256GHz (higher end of frequency range). It can be seen that the axial ratio of the pattern is deteriorated compared to those at the center frequency and lower frequency. This certainly shows that the number of inner loops still needs to be increased to insure a good quality pattern.....................97 6.1 Cross-section of an aperture coupled patch antenna, showing the cavity region I and the microstrip line region II for two different FEM computation domains....................... 100 6.2 slot and its discretization (a) slot aperture; (b) typical mesh from cavity region; (c) uniform mesh from microstrip line region........ 100 6.3 Illustration of a cavity-backed patch antenna with a coax cable feed. 110 6.4 (a) Side view of a cavity-backed antenna with a coax cable feed; (b) Illustration of the FEM mesh at the cavity-cable junction (the field is set to zero at the center conductor surface)............. 110 6.5 Field distribution in a shorted coax cable as computed by the finite element method using the expansion (6.18).: analytical; xxx: numerical. (a) Field coefficient eo along the length of the cable (leftmost point is the location of the short); (b) Field along the radial coordinate calculated at a distance A/4 from the shorted termination. ll 6.6 Measured and calculated input impedance for a cavity-backed circular patch antenna having the following specifications: patch radius r=13mm: cavity radius R=21.1mm; substrate thickness t=4.1mm; 6r=2.4; and feed location XfO0.8 cm distance from center. Results based on the simple probe model are also shown for comparison. Our modeling retains the vertical wire connection to the patch and uses the incoming coaxial mode field for excitation. (a) Real part; (b) Imaginary part...............................112 xii

6.7 Measured and calculated input impedance for a circular patch antenna having the following specifications: patch radius r=2cm; substrate thickness d=0.21844cm; feed location from center xf-=0.7cm; 6e=2.33; tan6=0.0012. -: measurement; xxx: this method; o o o: probe model. (a) Real part; (b) Imaginary part............. 113 7.1 Illustration of a shielded microstrip line............... 118 7.2 Illustration of the cross section of a shielded microstrip line. 121 7.3 A rectangular waveguide (a) and a microstrip line (b) truncated using the perfectly matched uniaxial absorbing layer........... 122 7.4 Plane wave incidence on an interface between two diagonally anisotropic half-spaces............................... 123 7.5 Typical field values of TE10 mode inside a rectangular waveguide terminated by a perfectly matched uniaxial layer.............. 128 7.6 Field values of the TE10 mode inside a waveguide terminated by a perfectly matched uniaxial layer. The absorber is 10 elements thick and each element was 0.5 cm which translates to about 13 samples per wavelength at 4.5 GHz....................... 128 7.7 Reflection coefficient vs 2/3t/Ag (a = 1) for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6... 129 7.8 Reflection coefficient vs 23_t/A\, with a -= 3, for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6... 129 7.9 Reflection coefficient vs 2/3t/Ag with a=l, for the shielded microstrip line terminated by the perfectly matched uniaxial layer........... 130 7.10 Reflection coefficient vs 2/t/Ag with a = 3, for the shielded microstrip line terminated by the perfectly matched uniaxial layer... 130 7.11 Input impedance calculations for the PML terminated microstrip as compared to the theoretical reference data...................... 131 7.12 Illustration of a meander line geometry used for comparison with m easurem ent............................... 132 7.13 Comparison of calculated and measured results for the meander line shown in fig.7.12............................ 132 Xlii

8.1 Illustration of the shielded microstrip stub excited with a current probe................................... 140 8.2 Impedance calculations using traditional FEM frequency analysis for a shielded microstrip stub shown in figure 8.1. Solid line is the real part and the dashed line denotes the imaginary part of the solutions. These computations are used as reference for comparisons....... 140 8.3 4th order and 8th order AWE implementations using one point expansion at 1.78 GHz are shown to compare with the reference data. With the 4th order AWE solutions, 56% and 33% bandwidth agreement can be achieved for the real (a) and imaginary (b) parts of impedance computations, respectively. It is also shown that the 8th order solutions agree excellently with the reference data over the entire band. (a) Real Part (b) Imaginary Part computations.... 141 9.1 Multi-modular FEM environment.................. 146 A.1 (a) A tetrahedron. (b) its local node/edge numbering scheme.... 148 xiv

Table 5.1 6.1 LIST OF TABLES Comparisons of gain and axial ratio at different operating frequencies 95 The correspondence between the edge numbers and the node pairs for each coordinate(r, 4 or z) along with the definition of the tilded parameters in (6.18)........................... 108 xv

LIST OF APPENDICES Appendix A. Evaluation of Matrix Elements for Tetrahedrals.................... 148 B. Evaluation of the Boundary Integral System Matrix.......... 152 C. Formulation for Right Angle Prisms........................... 154 D. System Derivation From A Functional....................... 157 xvi

CHAPTER I Introduction One of the primary goals in this dissertation is concerned with the development of robust hybrid finite element-boundary integral (FE-BI) techniques for modeling and design of conformal antennas of arbitrary shape. Both the finite element and integral equation methods will be first overviewed in this chapter with an emphasis on recently developed hybrid FE-BI methodologies for antennas, microwave and millimeter wave applications. The structure of the dissertation is then outlined. We conclude the chapter with discussions of certain fundamental concepts and methods in electromagnetics, which are important to this study. 1.1 Overview The development of simulation techniques for conformal antennas typically mounted on vehicles is a challenging task. By and large, existing analysis and design methods are restricted to planar and mostly rectangular patch antennas. These techniques have difficulty in being extended to non-rectangular/non-planar configurations loaded with dielectrics and comprised of intricate shapes to attain larger bandwidth and gain performance [1-4]. Moreover, practical antenna designs may also require a sophisticated feeding structure, such as coaxial cable, microstrip line, stripline, 1

proximity or aperture coupled circuit network, etc. and integral equation methods are not easily adaptable to modeling these structures, especially in the presence of finitely sized dielectric loadings. Partial differential equation (PDE) techniques (e.g. finite element and finite difference methods) may also experience difficulties in modeling unbounded field problems, such as those found in antenna radiation and scattering. The motivation of this dissertation is therefore based on the need to develop general-purpose analysis techniques which can accurately simulate conformal antennas of arbitrary shape with diverse feeding schemes. With the rapid growth of personal cellular GPS and other communication systems, there is an increasing need for such techniques since even traditional and protruding low frequency antennas are being re-designed for conformality and to meet requirements for a host of new applications [5,6]. Besides, most development of computational electromagnetics in this subject can be applied to medical diagnosis and treatment which have shown a tremendous research and application potential [7,8]. The complexity of new antennas demands that analysis and design software be developed based on methodologies that are robust, -versatile, and geometrically adaptable. Recently it has been demonstrated that the finite element method when coupled with the more traditional integral equation approach becomes quite attractive for modeling a wide variety of existing and emerging antenna configurations [9]. The finite element method is indeed ideal for modeling the interior volume of the antenna structure (multi-layer substrate, finite size dielectric loading, stacked element design, feed network and cavity volume, etc.) and is one of the most celebrated analysis methods in engineering. On the other hand, the boundary integral offers the most accurate representation of the fields exterior to the antenna. Thus the combination of the finite element and the boundary integral (FE-BI) methods provides

3 for the handling of the geometrical complexity without compromising accuracy. This hybrid methodology appears to be very attractive for conformal antenna modeling. However, its development and application to more practical and emerging antennas presents us with many theoretical and numerical challenges, which will be extensively investigated in the work. Specifically, mesh termination plays an important role in FEM simulations and, in many cases, the accuracy is subject to the performance of the domain truncation scheme. For conformal antenna modeling, a boundary integral (BI) equation has been employed in this dissertation for terminating the antenna's radiating surface and this method is theoretically free of approximation. Thus, a desired accuracy can be achieved without fundamental limitations. Antenna configurations of arbitrary shape can be readily tessellated using mesh generation packages in the context of the FE —BI technique. In modeling the interior region or the feed network, a superior artificial absorbing material perfectly matched layer (PML) has been used to ensure a minimum impact due to truncation walls. An intensive study of the PML's performance has been carried out and the optimal selection of PML parameters has been designed and employed herewith in shielded structure modeling. Frequency domain methods provide the necessary information for engineering design. However, when wideband responses are needed, they can quickly become expensive compared to time domain techniques. A method, referred to as the asymptotic waveform evaluation (AWE), can be used to alleviate this issue. It has already been successfully used in VLSI and circuit analysis. In the context of the FEM, we shall investigate the suitability and validity of AWE for simulating MMIC devices. One of the important issues in antenna analysis is the feed design. Modeling a feed using the finite element method is indeed a challenging problem, and a sinm

4 plified probe feed model fails to accurately predict the input impedance. On the other hand, the numerical system can become ill-conditioned when a feed network is modeled without careful considerations. In this dissertation various feed models will be investigated in consideration of accuracy and efficiency. They include current and voltage gap generators, stripline, microstrip line, coaxial cable, aperture coupled microstrip, etc. In regards to the development and applications of the simulation techniques, the test and design benchmark models of particular interest are microstrip (rectangular and circular) patch antennas, dual-stacked patch antenna, ring slot antenna, and cone antenna, etc. It is noted that some of them are not necessarily planar or conformal. Referring to the dissertation structure, we begin with a description of electromagnetic fundamentals and then procee(d to discuss the boundary conditions, equivalence principle, Dyadic Green's functions and the related theorems. The finite element method as applied to time-harmonic electromagnetic fields and waves is subsequently described and the basic FEM equations are derived from both variational and Galerkin techniques. The derivation is given in algebraic form allowing the inclusion of general anisotropy. The emphasis of the discussion is on the generalization of the variational functional and Galerkin techniques when anisotropic and lossy materials are present. Chapter 3,4,5 and 6 discuss the development of edge-based FE-13BI techniques with significant efficiency improvement for antennas and feed network modeling. The emphasis in these chapters is on developing novel methodologies to minimize the required computing resources. Chapter 7 is devoted to circuit modeling whe pecialized ltruncations suited for guide wave structures are presented. The perfectly matched layer (PML), an

5 anisotropic artificial absorber used for mesh truncation, is investigated in terms of performance and applications. Wideband system responses prompt us to look at more efficient analysis tools to replace the current brute force frequency domain analysis approaches. Chapter 8 discusses a preliminary development of the FEM in connection with the AWE. In the last chapter, we summarize and discuss the anticipated future research work to extend the capability and applications of the robust FEM development. A list of suggested topics is included with specific recommendations. 1.2 Fundamentals of Electromagnetic Theory Since many fundamental concepts and theorems of electromagnetics will be employed, we will describe the pertinent ones in this section for reference purpose. This will also ensure consistency in nomenclature and conventions throughout the dissertation. The vector wave equation -the only partial differential equation (PDE) considered in this research will be first derived from Maxwell equations. Various boundary conditions will be studied to establish the general mathematical models of boundary value problems (BVP). The equivalence principle, uniqueness theorem and the half-space dyadic Green's functions are then briefly discussed for EM solutions in radiation and scattering problems.

6 1.2.1 Maxwell Equations Time-harmonic Maxwell equations of differential form in a linear, anisotropic and uniform medium are given by [10] VxE = -jc -H-Mi (1.1) V x H = jcw E+J (1.2) V. e E = p, (1.3) V.uH = pm (1.4) where E and H are the electric and magnetic field intensity, respectively, w is the radian frequency and the factor ej&t is assumed and suppressed throughout this dissertation; M2 and Ji are the impressed magnetic and electric current, respectively, to serve as possible sources in the medium under consideration; finally pe and p, denote the electric and magnetic charge density. Both M- and Pm are fictitious and non-physical quantities, which facilitate the formulation of physical problems when the equivalence principle is employed. The material tensors e and Fu represent the permittivity and permeability, respectively, and may be written, in general, as E11 /12 C13 = Or = oC) 621 /22 C23 (1.5) C31 632 633 /Ill1 /112 /3I13 f-= /tObUr /= 0 /121 /122 /123 (1.6) /131 /132 /133 with c0 and /o being the free space permittivity and permeability.

7 The procedure to derive the vector wave equation begins by eliminating one of the two field quantities from (1.1) and (1.2). To do so, we first take a dot product of (1.1) with the tensor ur and then take the curl on both sides of (1.1) to obtain V x tT. V x E = -jwo(V x H - V x, Mi (1.7) Substitution of (1.2) into (1.7) yields V x (1 V x E) = w2o0Co E-jwoJi - V x ('1 Mi) or V x ( 1. V x E) -kor E = -j J-V x (U 1 M) (1.8) where k( = wU//0co is the free space wave number. The dual of (1.8) is given by V x (1 V x H) -kor H = -j oM + V x (- J (1.9) and can be similarly derived starting with (1.2). Equations (1.8) and (1.9) are the vector wave equations of the desired form. 1.2.2 Boundary Conditions and Boundary Value Problems Three types of boundary conditions are typically encountered, and in the context of the finite element method, these boundary conditions must be considered and carefully treated. In what follows we shall discuss these conditions. Dirichlet Boundary Condition Consider two media separated by a surface F whose unit normal n points from medium 1 to medium 2. The fields on two sides of the interface satisfy the relation n x (E2-E ) = -Ms (1. 1o)

8 where Ms is a fictitious magnetic surface current and E1 and E2 are the electric field inside medium 1 and medium 2, respectively. If medium 1 is a perfectly magnetic conductor (PMC), then E1 vanishes and (1.10) becomes n x E2 = -Ms. The surface magnetic current Ms can either be an impressed source (excitation) or may represent a secondary (induced) current. If medium 1 is a perfectly electric conductor (PEC), n x E2 also vanishes and thus Ms = 0 on the PEC surface. Similarly, for the magnetic field, n x (H2 -H1) = Js (1.11) where Js denotes an electric surface current. The PEC surface can support electric currents, given by n x H2 = Js, since H1 is zero within the conductor. By duality, the PMC surface does not support electric currents, i.e. V x H2. The relations (1.10) and (1.11) are inhomogeneous Dirichlet boundary conditions. They become homogeneous when Ms = 0 and Js = 0, and in those cases they imply the tangential field continuity across the dielectric interfaces. Often, Js and Ms are introduced as fictitious currents when applying the equivalence principle (except in special cases where they are specified a priori). The implication of this issue will be discussed later in the development. Neumann and Mixed Boundary Condition In formulating a physical problem using hybrid finite element methods, we usually work with either E or H field. If, for instance, we choose E as the working quantity, then (1.11) must be rewritten as n x [(1 V x E)- (1 V X E)1 =-jiJs (1.12) where (- V x E) i = 1,2 are evaluated just inside the ith medium approaching the boundary (from the ith medium). If medium 1 is a PEC, then V x E1 = 0, and

9 (1.12) reduces to a standard Neumann boundary condition, by which a constraint on the derivative of E at the interface is defined. The dual of (1.12) is given by ( )2 ( -1 x H)] -IM -j, (1.13) and this condition is used when working with the H field. In many applications, the single field formulation is often desired since the system size may be kept minimum in this manner. However, it is already seen that the single field formulation implies use of the second order conditions referred to as natural conditions. Fortunately, it is rather straightforward to impose this type of conditions in regard to finite element simulations. As for mixed boundary conditions, an example is the resistive surface where the electric and magnetic fields satisfy the condition + n xnxE+Rn x [H] =O (1.14) with R being the effective resistivity of the surface and [H] H+ - H- the field difference above and below the surface. This is a typical mixed (third type) homogeneous boundary condition. Another example of a mixed condition occurs in transmission line problem (e.g. a coax cable, or other guided wave structures), where the electric and magnetic fields at a cross-section of the line are given by E = Eie-Z + FEiez (1.15) H = HieZ - - Hie/Z (1.16) and n x Ei = -ZHi (1.17) where n = -z and (E, Hi) are the incoming fields before encountering a discontinuity or load along the transmission line. Also, Z is the wave impedance associated with

10 the transmission line mode of the guide wave structure. Eliminating F from (1.15) and (1.16) in view of (1.17) yields the relation n x E- ZH = 2n x Eie-~z (1.18) which is an example of inhomogeneous mixed (third type) boundary condition. This becomes apparent when H is expressed in terms of the derivative (curl) of E,1 and therefore the left hand side contains both differentiated and undifferentiated quantities. (In this case, the right hand side is usually considered as a known function.) The mixed boundary condition (1.18) is found very useful when applying the FEM to guided wave structures for truncation and excitation simultaneously. It is basically a form of absorbing boundary condition (ABC). 1.2.3 Uniqueness Theorem and Equivalence Principle The uniqueness theorem and the equivalence principle will be explicitly or implicitly applied to this work when dealing with integral equations to terminate the FEM mesh and when evaluating the far-field pattern. Together with dyadic Green's functions, it becomes convenient to apply these concepts to construct integral equations associated with various geometries in radiation and scattering problems. It is our intent to discuss the theorem and the principle (without proof) for later applications. Uniqueness Theorem Partial differential equations (PDE) can be solved using various approaches and the corresponding results can also be represented in numerous forms given certain boundary conditions. Moreover, many (boundary, initial, natural, essential, radiation, etc.) conditions of PDE models can be extracted from the mathematical 1Care must be taken when a curl operation is performed at a boundary discontinuity. It should be appropriate to evaluate the field derivative at a distance from a discontinuity and then let the distance tend to zero.

11 specifications of well defined physical problems. The question then arises as to how to relate the solutions and how many conditions are sufficient to achieve the 'correct' solution. Uniqueness theorems offer the answer to this question. Specifically in electromagnetics, the EMA solutions are uniquely determined by the sources in a given region plus the tangential components of the electric field on boundaries, or plus the tangential components of the magnetic field on boundaries. 2 Equivalence Principle From the uniqueness theorem, an EM problem can be uniquely solved if the tangential (either the electric or magnetic) field component at the boundary is prescribed. In this work, of interest is an EM problem where a dielectric inhomogeneous region exists in the presence of a large PEC platform, probably coated with a dielectric slab. The typical geometries are shown in fig. 1.1, where we consider the upper half space to be the exterior region and the cavity the interior region. In 1EM analysis, the fields in the exterior region can be represented in integral form containing the equivalent current sources. From (1.10) or (1.11) the tangential electric or magnetic field near the aperture (or the discontinuity region) may be equivalently expressed in terms of the surface currents Mi and/or Ji. By 'equivalence', we demand the field distribution remain the same when the fictitious surface currents are used to replace the interior region (cavity volume). It can be shown through the uniqueness theorem that this substitution indeed ensures an identical EM field distribution in the exterior region. When the interior region is excluded from consideration, the current sources in (1.10) and (1.11) may be arbitrarily chosen leading to an infinite number of choices for the equivalent currents. However, in our work the field behavior in the interior 2See the proof in reference [10].

12 equivalent currents ground plane cavity (a) equivalent currents 0' 0- - 0'- - cavity coated ground plane (b) Figure 1.1: (a) Recessed cavity in a PEC ground plane. (b) Recessed cavity in a dielectrically coated PEC ground plane.

13 region is also needed, and most specifically, the coupling of the fields in the inner and outer cavity regions is desired. It is therefore convenient to select the total tangential electric or magnetic field to specify the equivalent currents as M,=Exn and J = n xH. (1.19) This choice implies the assumption of zero interior fields when the exterior region is considered, and zero exterior fields when the interior region is needed. Fig. (1.2) and (1.3) illustrate the details of applying the principle, where the fictitious currents affect the region of interest (ROI) only with zero EM fields outside of the ROI. It is observed that this choice of equivalent currents permits a convenient interior/exterior system coupling for the "total field formulation" in hybrid FEM applications. 1.2.4 Integral Equation and Dyadic Green's Function The Dyadic Green's functions are particularly convenient for constructing integral equations in the presence of certain canonical platforms. For a planar structure, the platform of particular interest is the PEC infinite ground plane in which a cavity is recessed with dielectric loading or absorption depending on applications. The choice of the dyadic Green's function varies depending on the FEM formulations. For the electric field formulation, we are seeking an appropriate integral representation to find the magnetic field in the exterior region using the information on or near the region of the aperture. To this end, let us start with the structure containing a possible protrusion as shown in fig. 1.4, where the equivalence principle has been used on the outer contours of the structures to obtain the equivalent currents. Consider the wave equation V x V x G(r, r') - w2uocoG(r, r') = -I6(r - r') (1.20)

14 equivalent currents R01 (E,H)=(O,0) ground plane Region of Interest (ROI): Exterior (a) equivalent currents (E,H)=(O,O) cavity Region of Interest (ROI): Interior (b) Figure 1.2: Illustrations of equivalence principle when applied to the structure shown in fig. l.la.

15 equivalent currents RO - - w - -a~~aa aaraasaaeaaaasarrsanas aaaaarmanaaam (E,H)=(O,O) * ground plane Region of Interest (ROI): Exterior (a) equivalent currents (E,H)=(O,O) cavity ( Region of Interest (ROI): Interior (b) Figure 1.3: Illustrations of equivalence principle when applied to the structure shown in fig. l.lb.

16 Me Je (a) Me J e (b) Figure 1.4: Examples of protruding configurations (a) on a planar platform (b) on a curved platform in consideration of the equivalence principle. where G is the dyadic Green's function G in association with (1.9) (assuming Mi0), and I is the idem factor defined as I - X^ + + *y Z. Also, note the identity IJf{P(VxVxQ)-(VxVxP) Q} dV = - f [P x V x+( P Q +d (x Q] dS) and upon setting P H and Q = G, we get /J{H - (V x V x G) - (V x V x H). G} dV Me J e - f [H x V x G + (v x )G dS (1.22) From (1.9) and (1.20), the left hand side (LHS) of (1.22) reduces to LHS - H(r')- j V x J G(rlr') dV and the right hand side can be rearranged as RHS -II H [n x V x G + (V x H) - [n x G dS JJS

17 Equating the LHS and RHS yields H(r) V-JJ x J G(r' r) d' - {H [n x V x g] + (V' x H) [x x G- } dS' (1.23) where r and r' have been interchanged without loss of generality. As can be realized, V' is the volume containing the distributed electric current source and S' is the surface enclosing t[he entire upper half space. To eliminate the curl on J, we use the dyadic identity, V. (J x G) (V x J) G-J (V x G) and the divergence theorem to get Intolue =-J J (V x G) dV -J I (J x G) dS' (1.24) where the Sommerfeld radiation condition was invoked to eliminate the integral at infinity. Therefore, S' is only over the outer surface of the body. It remains to represent the surface integrals in terms of the electric field near the cavity since this field is typically the computable quantity. This is carried out by inserting (1.24) into (1.23), yielding H(r) = - J (V x G) d1' J {H. (n x V' x G) + ( x H- J) (n x G)} dS' J= jJ-(V'x G)dl' - J/j {H. (n x V x G) + jweoE. (n x G)} dS1 (1.25) where the Maxwell equation (1.2) has been used. It should be remarked that the above field representation is general, i.e. not restrictive to planar or conformal cases.

18 For instance, in the presence of a PEC platform, (1.25) is valid for protruding configurations as shown in fig. 1.4. In these cases, the surface integrations are carried out over the platform plus the outer contours of the structures. The field representation (1.25) shall be examined and compared for conformal and protruding structures. To this end, we rewrite the surface integral as Intsztrfa/ce = (V x G) (n x H) dS' + V x (V x G) (n x E)} dS' x H ( x G) dS' - (n x E) (V'x G) dS' (1.26) JJs' I /o J where T denotes a transpose operation of the dyadic and the integral in the last step is proportional to the electric field. If G(rlr') is the electric dyadic Green's function of the first kind defined as n x G = 0, the first term in the integrand of (1.26) vanishes on the platform, provided S' is coincident with the platform. For dielectric protrusion, this term reduces to the integration only over the outer contour not conformal to the platform. An alternative is to define an electric dyadic Green's function which satisfies the condition n x (V' x G)T 0. As can be seen, this definition of the Green's function equivalently leads to the same vanishing term in (1.26). G is referred to as the dyadic Green's function of first kind. The equivalence of both definitions can be proved from the symmetry properties of the dyadic Green's functions [11]. For a planar PEC platform, G reduces to G(rlr') = Go(rr')- Go(rlr') + 2$3Go(rlr') (1.27) where Go is the free space Green's function given by o-jkorr-r'\ Go( lr') - -r 47fir - r'[

and Go(rlr') (I + -V);o(rlr') Inserting (1.26) and (1.27) into (1.25), we obtain H(r) = Hi(r) + Hr(r) + jkYo J Go(r r') (n x E) dS' (1.28) This is the desired form of the magnetic field representation used to establish the boundary integral equation for a planar platform.

CHAPTER II Finite Element Analysis in Electromagnetics The finite element method (FEM) has been applied to electromagnetics (EM) since several decades ago [12]. Especially in the late eighties and early nineties, it is observed that the publication volume associated with the FEM in electrical engineering grew in a fairly rapid pace [13]. This is primarily because electromagnetic problems in engineering designs become increasingly complex and analytical approaches or other numerical techniques no longer meet practical needs. With its numerous attractive features over other numerical techniques, the FEM has been extensively investigated and exploited for various EM applications [13]. This chapter is organized as follows. Section 1 and 2 describe the theoretical formulations to construct the FEM equations. These are usually considered the indispensable fundamentals of the technique, even though some interesting issues associated with these basics are still in development stage, especially in terms of numerical implementation. Of interest in this context is the discussion of the variational functional and Galerkin's techniques when applied to general anisotropic and lossy electromagnetic problems. This topic is one of the least studied and documented in the literature related to computational electromagnetics. Anisotropic materials have been used for domain truncations (refer to Chapter 7) and therefore the general 20

21 vector or tensor form will be used whenever possible in the necessary derivations for the FEM. With this type of formulations, isotropic systems may be regarded as special cases. The chapter is concluded with the discussion of the physical quantities for antenna analysis in association with the computation of electromagnetic fields. The formulas given in this context require minimum amount of effort for computations. 2.1 Functional Formulation The FEM was first developed with the aid of functional analysis. Traditionally, many standard boundary value problems (BVP) encountered in practice can be equivalently related to the extremization of a certain variational functional. With the Rayleigh-Ritz procedure to project a continuous function space onto a discrete finite expansion space, the variational functional method can be used to solve those physical problems and therefore becomes one of the two important approaches to formulate the FEM. A functional version of the FEM for the vector wave equations (1.8) or (1.9) is discussed in this section, which can readily incorporate boundary conditions, sources, resistive cards and other constraints into the formulation. It is regarded as a natural, convenient and sometimes physically meaningful approach. Furthermore, the functional may represent a true physical quantity (e.g. in low frequency, power transmission applications) and hence this formulation provides a feature of merit for its evaluation. Also, as can be seen in this chapter, the variational method in general non-self-adjoint cases may be rigorously treated to result in a final symmetric system, a subset of which is identical to that obtained from Galerkin's technique. Last, but not the least, the variational functional formulation can be used to validate the expressions based on Galerkin's method.

Consider a typical radiation or scattering problem shown in fig. 2.1, where the z A Radiating element Aperture ~. Ground plane Y x Figure 2.1: Illustration of a typical conformal antenna configuration. radiating elements (or array) are enclosed in a region Q. The platform surrounding the radiation/scattering geometry can be a planar ground plane, certain canonical shape (cylinder/sphere), or even a doubly curved surface in which case the Green's function is not available. In Q, electromagnetic fields satisfy the wave equation (1.8) or (1.9), which can be concisely described using a linear operator C given by ~C = Ki (2.1) where 1 denotes the field E or H, and L (vx.V1 Vx)-(kxC-.) L (Vx v 1 * Vx)-(D jr-) for electric field for magnetic field (2.2) (2.3) K is the source term associated with the impressed electr aic and magnetic currents and may be explicitly given by Ki = -jc/oJi- V x (1r * Mi) Ki = -jo6o0M + V x (i.1 * Ji) for electric field for magnetic field (2.4) (2.5)

23 As already mentioned, C is a linear operator and one can readily show that for symmetric dielectric tensors p and (, the pertinent functional of the original PDE has the form <(~) = < <,) > - <, K; > (2.6) where the inner product <, > is defined as <A,B >= A-B*dV (2.7) (with B* being the complex conjugate of B) for lossless media, or more generally as < A,B >= j ABdV (2.8) for both lossless and lossy media. The equivalent variational problem can now be stated as the extremization of the functional (2.6) in conjunction with the essential boundary conditions (e.g. Dirichlet BC's). Specifically, the boundary value problem is equivalent to the following variational model an(b) = O (2.9) Essential Boundary Conditions Because the effect of complex materials on the resultant system is of primary interest to us, we restrict most of our discussions in this chapter to homogeneous Dirichlet boundary conditions unless otherwise specified. Therefore, the variational approach of (2.9) ensures a symmetric numerical system. This is significant since many physical problems retain a certain symmetry property and the corresponding mathematical models should therefore reflect this property. Moreover, the symmetry of a numerical system is always desirable since it leads to more efficient solution and less storage requirements.

24 The existence of the functional (2.6) requires the operator C be self-adjoint, a property usually defined to satisfy < C, P >=< C,L~ > (2.10) where <D and I represent any two admissible functions. If (2.10) holds, not only does the numerical system derived from the functional (2.6) remain symmetric, the minimization/maximization becomes physically meaningful. Of most concern is the situation where the partial differential operator of a system is no longer self-adjoint. Mathematically there exists no such a natural functional in the case similar to (2.6). A typical example is the presence of a lossy and anisotropic medium, whose dielectric material tensors are not symmetric or Hermitian, and this type of problems is more often seen nowadays. The development of finite element methods for those problems is still at an early stage because it involves numerous challenges. Traditionally, these physical problems were fictitiously simplified and dealt with using available numerical approaches. Konrad [14] first tried to formulate a 3-D FEM with three vector components to represent electromagnetic fields in anisotropic but loss-free media. The tensors were therefore assumed to be Hermitian in his study. A few years later in 1980's, the number of publications in this subject increased typically with applications to waveguide structures. Unfortunately, the variational approaches reported by different authors during that period consistently led to nonstandard and non-Hermitian eigenvalue systems (even with the aid of an adjoint system [15, 16]). Even worse, the numerical systems derived in this manner were usually doubled in size. As indicated in [17], when a non-standard eigenvalue system was manipulated to reduce to the standard form, the size of the system was doubled

25 again. Similar rep)orts were also seen in later papers [18-23] for the problems of wave propagation inside anisotropic media. No radiation and scattering analysis has been reported in this context. In what follows, we generalize the FEM formulation to lossy and anisotropic electromagnetic problems. Specifically, we show that two different methods, one with the aid of an adjoint auxiliary system and the other with the Lagrange multiplier, can be used to construct the pertinent functional. Galerkin's method is then compared to these two variational techniques. 2.1.1 Pertinent Functional For Lossy/Anisotropic Media - I As is known, the natural variational functional no longer exists for non-Hermitian operators since no matter what definition is given for inner products (see (2.7) or (2.8)), one cannot obtain a self-adjoint operator necessary for natural functional design. In these cases, we consider a generalized functional SF:< L@ > - < (),Ka> - < VK > (2.11) where '( is the unknown solution function of the original PDE problem and K is the right-hand-side function as in (2.1). Similarly, 1 is the solution function of the adjoint PDE such that Cay = Ka (2.12) where La can be derived from < ~LA, >=< KCLa, > (2.13) with L 7: La. It should be functionthalt the functional (2.11) reduces to (2.) reduces to (2.6) (except for a possible constant coefficient) if L is self-adjoint. Also, the original PDE and its

26 adjoint counterpart (2.12) can be recovered through the variational process with respect to the functions 1 and 4, respectively (the simple derivation is omitted here). After discretization is carried out, the final numerical system is symmetric. This can be shown as follows. Let = xiVi, - = yjVj (2.14) i j where V, is the basis function used for both unknown functions and xi, yi are the corresponding expansion coefficients. Inserting (2.14) into (2.11) yields S -F= H xy < LVi,Vj > -Ei < Vi,,Ka > - yj < Vj, K > (2.15) i J Upon performing the differentiation with respect to xi and yj individually, we get the two decoupled systems of linear equations Q(2.16) CY 0 y KY where the matrices Q, Q{Y and the column vectors KX, KY are given by Q^ = <~2Vj,Vi> Q3 = <~rVi,Vj > K: = <V V,K> Ky = <Vi, Ka > In general, Qx Q ji Q QYi and QZ + QC. However, Qx = QY. (QY)T. These relations indicate a loss of symmetry of the original problem, but the symmetry holds for the overall system! The storage requirement is a function of N/2, where N is the dimension of (2.16). Even though there is an auxiliary system needed to complete the analysis, in practice this system does not require storage.

27 2.1.2 Pertinent Functional For Lossy/Anisotropic Media -II An alternative to using an adjoint system is to employ the Lagrange multiplier technique in constructing the pertinent functional. The Lagrange multiplier is usually used to incorporate additional constraints to a system. To illustrate the technique, consider the same PDE model as described in (2.1), and we first rewrite it as ~D-K;- =0 (2.17) Next, we assume an expansion space where the solution is defined and solved. The unknown function ) and the multiplier function A are expanded in the same space. If (2.17) is regarded as a "constraint", we try to add the constraint to a "null" system and get (, A) =< A, A -K, > (2.18) This functional is now used to formulate the FEM. As described above, on applying the Rayleigh-Ritz procedure to both () and A using the same set of basis functions, viz. = ExiV, A = - y (2.19) i 3 we obtain F(1, A) = xiyj < Vj, CVi - K, > (2.20) i 3 Carrying out differentiation with respect to xi and yj, individually, yields ) ( ) (2.21) where Qs = QY as in (2.16). We observe that (2.21) is similar to (2.16) with two decoupled subsystems of the same size. The properties of the subsystems are also

28 similar to those in (2.16). We further observe that the Lagrange multiplier technique may be regarded as a special case to (2.11) where a homogeneous adjoint PDE is now virtually assumed (i.e. Ka = 0). Again, it should be remarked that when a self-ad joint problem is considered, the above formulations (both the adjoint system approach and the Lagrange multiplier technique) will result in a symmetric system and the auxiliary subsystem becomes either redundant (for adjoint approach) or unnecessary (for the multiplier technique). In the later case, the multiplier can be considered identical to the unknown function ' expanded using the same basis functions. 'This results in an interesting coincidence with Galerkin technique described next. 2.2 Galerkin Formulation Galerkin's method is now considered to formulate the finite element method. Traditiona]ly, Galerkin's technique used in conjunction with integral equation employs the same testing and expansion functions to obtain a symmetric dense numerical system. However, in the case of the FEM, Galerkin's method does not always lead to a symmetric system. Apart from boundary conditions, the linear operator of a PDE problem determines the symmetry feature of the resulting system. Also, unlike the variational approach, Galerkin's method solves the weighted PDE by a testing process as < P,~L >=< W,Ki > (2.22) where d) and P are )both defined in the same function space. Specifically, in Galerkin's method one seeks the solution for the unknown function ID which satisfies certain prescribed constraints and (2.22), with the aid of another arbitrarily chosen function T.

29 Similarly to the variational approach, the Rayleigh-Ritz procedure may also be used to project the continuous space onto a finite discrete separable Hilbert space. The mathematical problem is then rephrased to seek a discrete solution set whose entries are the coefficients of the expansion. The testing function 1 must, of course, be defined in the same discrete space to ensure that the original PDE is solved with proper boundary conditions. It is obvious that if the linear operator C is self-adjoint, the choice of the testing function. = 1 results in a symmetric numerical system. Otherwise, no matter how the inner product is defined, the final discrete system in general does not exhibit symmetry property. We observe that Galerkin's method can usually be applied to any linear operators even when the corresponding natural functional does not exist. Also in the general cases (as considered when describing the functional approach), Galerkin's method leads to the same numerical system as the desired portion in (2.16) and (2.21). This can be demonstrated as follows. Inserting (2.14) into (2.22), we readily obtain A_ A x.yj < Vj, ~Vi >-= y. < Vj, Ki > 3i j or yji{ ZX < Vj, Vi > - < Vj K>} = (2.23) As assumed, I is an arbitrary function. Thus the term in the curly bracket should vanish, yielding Qx = Kx (2.24) which is exactly the same as the subsystem derived from the variational approach.

30 It is noted that if one chose the Lagrange multiplier (in the variational approach) as the testing function I, the same numerical system would result. Analytically, though, the entire system is twice the size of that derived via Galerkin's method. This is due to the fundamental difference of the two techniques. Regardless, as expected, the obtained numerical systems of interest are virtually identical! 2.3 Total Field and Scattered Field Formulations In this section, we focus on a general scattering problem as illustrated in fig. 2.2, where a perfectly conducting electric (PEC) body is coated with a dielectric layer whose relative permeability and permittivity are Jid and (d, respectively. (Note that d Fo Qd: dielectric coated region (td,dd) 'j: free space ( = - i- = 1) Qa: absorbing layer (~a,P7) Fp: boundary of the PEC body Fd: boundary of the dielectric coating and free space region Ff: boundary of the absorber and free space region F0: PEC boundary of the outer absorber Figure 2.2: Illustration of a scattering problem setup for scattered field formulation. for the purpose of generality, the medium is assumed anisotropic.) The situation with absorbing boundary conditions for truncating the FEM domain has been analyzed before (see e.g. [24] or [25]). However, two issues associated

31 with this type of problems have not been carefully addressed in the literature. One of them is the equivalence between the variational and the Galerkin's method when the scattered field is used as the working variable. Proof of this equivalence can be tedious and cumbersome but it is nevertheless an important issue. Unfortunately, one is used to assuming that these two formulations are equivalent without proof. Another issue relates to the recently introduced perfectly matched absorbing material. As it turns out, there are several advantages to use artificial absorbing materials, including accuracy control, conformality, ease of boundary treatment, etc. However, their inclusion introduces additional artificial conditions inside the absorber layer and care must be taken when those conditions are enforced in the FEM formulation. Moreover, although a metal-backed absorber layer simplifies the FEM implementation, the multilayered FEM region contains high inhomogeneity, which again requires a careful presentation of the formulation. To this end, we extend our theoretical discussions on the FEM to scattered field representations, where the treatment of the boundary and transition conditions will also be described. 2.3.1 Scattered/Incident Fields and Boundary Conditions Referring to fig. 2.2, we begin with the wave equation in terms of H (the E formulation may be readily handled by duality). To proceed with Galerkin's method, we first write Ht~ as Htt = Hscat + Hinc (2.25) where Hat and Hi'c are the scattered and incident field, respectively. Next, weighting the source free wave equation with the testing function V yields V V {Vx -Vx H-k r H} dQ 0 (2.26)

32 where Q = Qd + Qf + Qa encompasses the entire computational region. To proceed with the derivation of the weak form wave equation, it is necessary to introduce certain constraints on the scattered and incident fields within Q and on the boundary. First, since the incident field is not allowed to pass through the absorber layer as well as the metal back wall Fo, we note that scat + Hinc r C Qd + Qf Htot(r) - (2.27) Hscat r C otherwise with the incident field satisfying { -1 = 0 r E, + Qf {Vx r Vx-kor.H { rf (2.28) 0 r otherwise It is thus evident that the scattered field satisfies the homogeneous wave equation in regions Qa + Qf and the inhomogeneous wave equation in Qd. The boundary conditions on Hsat can be readily derived by consistently applying the field decomposition (2.25). Note that in accordance with (2.27), the electric field is likewise decomposed as Etot Eat + Einc (2.29) However, one should be cautioned that EinC and Hinc do not satisfy Maxwell equations in the dielectric region. That is, jwcoEinc - V x H2nc r C Qd (2.30) which conflicts with what one would intuitively assume. Conventionally, the incident field is assumed to exist in the dielectrically coated region Qd as if there was no dielectric there. After a quick glance, one would immediately arrive at a conclusion that (2.30) is against Maxwell theory. In reality, it can be proved that if Einc and

33 Htic indeed satisfied Maxwell equations in fd, a contradictory boundary condition would immediately result. This can be seen by taking a curl operation on both sides of (2.25). In view of (2.29) and Maxwell equations in dielectric media, we would have Cd V x HtOV x H V H + * r (2.3) Imposing the condition of total field tangential continuity at the boundary Fd would yield n X { H1 * V X H car+ - - V x H t 0 (2.32) which is a homogeneous Neumann boundary condition. However, if we start with (2.27), upon taking the curl operation and imposing the condition of total field tangential continuity at the boundary Fd, we get dx { v - V x Hca tr+ - e1 Vx H r- = -n x { - V x Hi (2.33) which is an inhomogeneous Neumann boundary condition. This inconsistency is because (2.32) was derived on the basis of the decomposition (2.29) and the assumption that the incident field within dielectric regions also satisfies Maxwell equations. However, the decomposition (2.29) is artificial and it is therefore necessary to keep in mind that only (2.27) holds true when deriving boundary conditions. As a rule of thumb, an appropriate interpretation of the phenomenon should read: the incident field inside dielectric media existed in the same fashion as in free space as if the free space was replaced with the nmedia. Mathematically, this implies the condition JieoEi = E. V x HTn r C Qd (2.34)

34 Consistently one can derive the other boundary conditions required for the FEM formulation following the same procedure. They are classified as Dirichlet and Neumann conditions as follows. * Dirichlet Conditions Boundaries Conditions p n X Hscat = - n x Hic (Ks unknown current) rd x {Hscat+ - Hscat I} Ff n x {Hsc - Hscat - -n x Hnc or n x Hsat Kb (Ks unknown current) * Neumann Conditions Boundaries Conditions =-1.X Hscat -1 Hic Fp nX 6d V x V x H = -V x H + [+ Fd n x I V x H - +h x x -V x H Ff n x {. V x HSCat} n x. V x Hc Fo n x 71 'V X HSat = where {} + denotes {} +-{} - and -r and HSCat should take the values at positive and negative sides of a specific boundary. 2.3.2 Galerkin's Method Returning to (2.26) and in view of the vector identity A V x B = (V x A). B - V(A x B) (2.35) and the divergence theorem, we obtain the corresponding weak form wave equation Intd + Intf + Inta = ( (2.36)

35 where I11td [V x V- 1 V x (Hsca + Hic) - koV * y (H + Hn] d d +I V. V x e * V xH dS + j V. (hi x. V x Hc) ) cS V.d X X Hsca dS (2.37) Intf -- j x VX-V 1 ~ VxHct dQ+ V- - x j1 ~ VXHst dS + V (n^ X f H V x H ) dS (2.38) Inta = V x V a 1 V x Hscat d + jV (-2 x 1 * V x HSt) dS + j V (n3 x,1 V X Hs) dS (2.39) with hno i1, n2 and 23 being the unit normals at the boundaries Fp, Fd, If and Fo, respectively. They all point away from the center of the PEC body (i.e. outwards). Invoking the boundary conditions as tabulated above, we observe that the surface integrals on Fp in (2.37) and on Fo in (2.39) vanish. Also, the sum of the surface integrals on Ff in (:2.38) and (2.39) reduces to — n - V. n2 x f1 'V x Hinc dS (2.40) Similarly, the sum of the surface integrals (involving Hat) on Fd in (2.37) and (2.38) becomes V * x (f1 -() V x Hi dS (2.41) Note that both integrals (2.40) and (2.41) will not contribute to the system, but to the excitation on the right hand side. Of course, the excitation term (2.41) would disappear if the permitivity is continuous across the boundary Fd. By duality, for an electric field formulation, a term similar to (2.41) will appear for discontinuous

36 permeability. This observation does not hold for (2.40), which arises from the inherent incident field definition on the two sides of the boundary Ff, as mathematically shown in (2.27). Furthermore, the second term in (2.41) will be cancelled out with the integral on Fd involving Hinc in (2.37). After a simple manipulation, the final system reduces to (v x r V x H"ct - kV * ' * Hs^c d) -J ( v - VxHic Voi Hir Hnc) dQ - I n. x V xHin dS+ V 2 x e1 V xHin dS (2.42) f where r is again the relative permittivity with respect to the specific region. This is the desired weak: form of the wave equation, and a similar equation was obtained in [24] using the functional formulation for isotropic media. 2.3.3 Variational Method It is now of interest to employ the functional formulation to obtain the equation corresponding to (2.42). To this end, it is intuitive to begin with the total field representation and use the functional (H) = (v x H V x H - kH H) dQ (2.43) where, as before, T, and P, denote the corresponding relative permittivity and permebility, respectively. Also similar to Galerkin's method, one would now logically proceed with the field decompostion H = H""t + H0c and express the functional in terms of the scattered field. It is unfortunate that this approach will result in a different form of linear system than (2.42) and the detailed mathematical proof is presented in Appendix D.

37 The discrepancy arises from the assumption of the functional (2.43), which is not a valid expression when the corresponding PDE operator is no longer self-adjoint. Recall in section 2.1.1 that the presence of general anisotropic and lossy media requires the application of an auxilary adjoint system, together with which the pertinent functional is given by a generalized form in (2.11). It is therefore necessary to begin with this generalized functional rather than (2.43). In a source-free region, (2.11) explicitly takes the form F(HaH)= Ha V x.1. V x H- ko, H} dV (2.44) and the variation is imposed on the adjoint variable Ha to get 6T(Ha,H) =0 (2.45) 6H=0 The variational functional method of this version proves valid and identical to Galerkin's technique for any linear operators in electromagnetic problems. It is also interesting to note that once compared to Galerkin's method, the adjoint field quantity Ha in (2.44) seems to take place of the testing function V in (2.11). However, Ha theoretically differs from V in that the former is defined as a solution function for the adjoint system of the original PDE problem, whereas the latter is just an arbitrary admissible testing function that does not have to be a solution to the adjoint system. Apart from the concept difference as mentioned above, the mathematical procedure required to derive the FEM system is similar to that presented in the previous subsection, and will not be discussed here to avoid repetition. One can be assured that the final system obtained from (2.44) and (2.45) is of course identical to (2.42).

38 2.4 Parameter Extraction An accurate full wave analysis can only predict the near field (for PDE type methods) or current (for integral based techniques) distributions, which can be used to obtain certain practical parameters depending on applications. For intricate systems, involved numerical models may be needed for output data extraction, including far field evaluation in the presence of non-canonical platforms and a de-embedding process for antenna feed network or circuit simulations. Antenna parameters can be readily evaluated after the near field distribution is achieved via full wave analysis. The de-embedding process is required for feed network or circuit modeling and will be discussed in chapter 7. For a non-planar platform the far field evaluation can be obtained from the general discussion in chapter 1, where the formulation i terms of the dyadic Green's function must be used to consider the equivalent currents and the free frespae space Green's function. 2.4.1 Radiation and RCS Pattern In the case of antennas, we are mostly interested in their radiation and scattering patterns and other related parameters such as gain and axial ratio. (The near field quantities such as input impedance, feedline S-parameters, etc., will be disccussed in later chapters). Both radiation and radar cross section (RCS) patterns can be readily characterized with respect to the 3-D spherical coordinates 0 and b. Consider the planar cavity-backed antenna as shown in fig. 2.1. Once the field distribution on the aperture S of the conformal antenna is obtained from the full wave analysis, we can then proceed to evaluate the far field pattern. The most straightforward approach for this computation in the presence of an infinite conducting ground plane is use of the equivalence principle. To do so, we define the magnetic current

39 as M = -z x Es, where Es is the electric field evaluated on the aperture and Z is normal to the aperture surface. The electric vector potential is then given by F(O, ) - o It MeJkr' dS' 47rr JJs C-jk r/ = IO Jj(Es x 2)eikr dS' (2.46) 47rr JJS where the electric field is typically expanded in terms of the surface vector basis function S'. Introducing this expansion, (2.46) becomes F(4) - J S x ejk.r' dS' 4 e. -V (2.47) 47rr where V denotes the sum of the surface integral over each of the discretization elements on the aperture. The far zone magnetic field H becomes HT = -jWFT (2.48) where HT represents the transverse component of the far zone magnetic field, whose 0 and b projections are given by 3kr He = -jw — {cos 0 cos OVA + cos 0 sin )Vy - sin VVz} 47rr WCO j-3kr = P0 47rr 47rr in which Po =-j {.} and P1 =-j {-}, and {.} stands for the corresponding terms in the curly brackets. The RCS of the t (t = 0 or 4) component can be represented in terms of Po and P% to yield rrcs=4722 (2.50It2 rtQ- 47r2 rH2 = _ 2 - Z IPL (2.50) (Hil2 4A 2 02I

40 where A and Z o = Co are the free space wavelength and the intrinsic impedance, respectively. In (2.50), the incident wave IH2I was set to unity without loss of generality. In practice, it iis customary to express RCS in dB, where at is usually normalized to A2 or to squared meter first. The radiation pattern may also be represented in the same manner as mentioned above. The difference in procedure here is the normalization with respect to a maximum radiation field value. The reason is that in radiation mode, the antenna is excited by an interior source rather than the incident plane wave Hi as in the scattering case. Of interest therefore is the relative field intensity in the far zone. To get this, we represent the far field intensity in terms of the above calculated quantity crcs to avoid a repetition of post processing. Specifically, the formula grad- t (2.51) (o rcS)max is used for radiation analysis. 2.4.2 Gain and Axial Ratio Gain (G) and Axial Ratio (AR) are two important parameters which indicate the antenna's performance. It is also noted that these two quantities typically characterize the far zone features of the antenna. By definition, the gain G for a lossless antenna (with 100% efficiency) is given by G = Um (2.52) Prad for a cavity-backed structure, where Umax is the maximum power density and Prad denotes the total radiated power from the antenna in the upper half space. (It is noted that this definition of G is identical to that of directivity for lossless antennas.)

41 Since the maximum power density Lmax can also be expressed as Umax Z= (0 + ) (2.5;3) the antenna gain becomes G - Zo(o + go) (2.54) 2Prad Thus, the computation of G is rather trivially done once a is found as given in the previous subsection. In reality, Prad may be evaluated on an assumption that all input at the feed is transferred to antenna element(s) and radiated. In this case, one has Prad = I2R, where I is the known current source on the feed, and Rin is the input resistance of the antenna measured at the reference plane. However, this scheme may not always work since the gain (more precisely directivity) reflects the far field behavior, while the input resistance computation relies on an accurate full wave model for near field prediction. It is well known that the far field and the near field computations offer different accuracy. This accuracy inconsistency arises if Rin is used to determine the gain G, a far zone pattern whose accuracy is directly governed by the near field computation without averaging effect. To avoid this accuracy inconsistency, one may calculate the gain or the directivity by evaluating Prad from the far field radiation pattern. Specifically, Prad can be obtained by integrating the radiation intensity Prad J= / U dQ = o jI( + 0+ ) dQ (2.55) over the half space. It is obvious the half spaa cesrtain symmetry of the pattern remains such as circular about the vertical Z axis, then U(O, I) reduces to UI(0) and the above

42 integration can be effectively evaluated. Otherwise, a numerical 2-D integration is required. Axial ratio (AR) is another important antenna parameter, especially when a circular polarization (CP) of antenna's performance is of the primary concern. Since AR also features the far zone effect of the antenna, it is desirable to determine AR with a minimum computational load. It is noted that given the above pre-calculated au and (O, one is again able to determine AR uniquely by AR = on (2.56) Vshort with Hiong = 2 + 7 + [<T4 + q + 2(72 COS 13 vzl —n — = [u +o~+2o-cosf } (2.57) Vshort + u2 - [e + u4 + 220, cos 33] } where = 2(1H0 - IH,), the twice of the phase difference between the two magnetic field components. It can be obtained from the quantities Po and PO defined in (2.50). Since these two quantities are complex numbers, the phase difference is readily represented as 3 -jm{ n P (2.58) 3 Po with Im { } being the imaginary part.

CHAPTER III Edge-Based FE-BI Technique 3.1 Introduction Numerical methods have been serving the engineers and researchers for many years in antenna analysis and design. Among them, the moment method in conjunction with various integral equation (IE) formulations played a major role [1-3]. However, IE methods are associated with field representations in which the appropriate Green's function for the specific geometry must be employed and this limits their versatility. Moreover, IE techniques are usually formulated on the assumption of an infinite layered (not inhomogeneous) substrate, a model which deviates from the practical configuration and leads to inaccuracies for larger bandwidth antennas. Furthermore, in the context of IE methods, antenna excitations are represented using simplified models that differ more or less from the actual configurations. Also, due to the singularity of the current distribution near the feed junction(s), special measures must be taken [26] for proper modeling. In contrast, the hybrid Finite Element-Boundary Integral (FE-BI) technique alleviates these difficulties and this was recently demonstrated when the method was applied to inhomogeneous objects of canonical shape scattering [27,28] and rectangular patch antennas [9]. Based on the past success of FE-BI methods for antenna analysis, it is desirable 43

44 to extend the method to antennas of arbitrary shape. In this chapter, an edge-based hybrid finite element-boundary integral formulation is presented for the characterization of arbitrarily shaped cavity-backed antennas [29]. An example of such a configuration is shown in fig. 3.1, where a cavity is recessed in a metallic ground plane enclosing the FEM volume. The antenna elements on the aperture may be excited by different schemes, such as a simple probe, a magnetic frill generator, a practical coaxial cable, microstrip line, slot or a CPW line. In the context of the FEM, the cavity is first discretized into a number of tetrahedral elements that naturally reduce to triangles on the cavity's aperture. For non-rectangular patches this triangular gridding is, in general, non-uniform and the exact boundary integral formulation based upon this mesh applies to any patch shape. As a result, the hybrid FE-BI technique is capable of modeling arbitrarily shaped cavity-backed antenna configurations, different substrate inhomogeneities, anisotropies, as well as various practical excitation schemes. 3.2 Hybrid System Functional In this section, the edge-based hybrid FE-BI method will be formulated using the variational principle, where matrix algebra notation is employed so that one can readily extend the formulae to the general anisotropic case. As presented in [9], the complete functional pertionnent to the scattering and radiation by a cavity-backed configuration shown in fig. 3.1 may be written as

45 z A Patch Ground plane Aperture _ -.... --—...... —^ Cavity Coax cable opening (surface C) X -- ~Coax cable Figure 3.1: Illustration of a typical radiation and scattering problem. F(E) = 2 tx {(VxE) -(V x E)-k(,E. E dv +2jikoZo f(E x Hi) zdS + t E jkoZoJ + V x M) dv (3.1) 2k fj(E x z) { (E x i). ( + -VV) G(r, r') dS' dS where Ji and Mi represent interior electric and magnetic current sources within the cavity V; Ht is the incident field, if any, from the exterior region; the surface S encompasses the cavity aperture excluding the portion occupied by the antenna elements; Or and pr denote, respectively, the relative permittivity and permeability; ko is the free space wave number, I the unit dyad, and Go(r, r') the free space Green's function with r and r' denoting the observation and integration points. 3.2.1 FEM Subsystem In proceeding with the discretization of (3.1), it is convenient to decompose it as F = Fv + Fs (3.2)

46 where Fv denotes the volume integral contributions and similarly Fs accounts for the surface integral contributions. The cavity volume is subdivided into N tetrahedral elements Ve (e=I1,2,...N), and within each tetrahedron as shown in fig. 3.2 the field 1 2 3 nodes/vertices Figure 3.2: A tetrahedron and its local node/edge numbering scheme is expanded using edge-based elements as E = [V]T{E}e with (3.3) [V}e {E}e = [{vx}{ v}{ VZ}] Vu1 Vu2.u U E2 6, = x,y,z (3.4) \ / e

47 in which V,i is the u (u = x, y or z) component of the volume vector basis functions along the ith edge. The unknown vector {E}l has six entries, one for each tetrahedron edge. (Here we use square brackets for matrices and curly brackets for vectors). Inserting (3.3) into (3.1), and taking the first variation of Fv with respect to {E}, yields Fv =- {([A]e{E}e + {K}e} (3.5) e where [A]e = {J { [DV]e[DV - kr[V]e[V]} dv (3.6) Ji x Mix {K}e = - [V]e jkoZo j( + vV x } dv (3.7) J - Mi ay aa1 - V z} {Vy} [DV]T = av} k | (3.8) az] {V -,a{vx} a9 8 To carry out the above integrations, it remains to introduce the volume expansion or shape functions Ve. For our implementation we employed the linear edge-based shape functions for tetrahedral elements given in [30,31]. The explicit finite element matrix entries associated with a typical tetrahedron (as shown in fig. 3.2) are given in Appendix A for reference.

48 3.2.2 Boundary Integral Subsystem To discretize the surface integrals in (3.1), the aperture is subdivided into triangular elements since these correspond to the faces of the tetrahedrals. Within each triangle, the field is represented as E =[S]T{E}e (3.9) where [S]e [{SX}{Sy}] Sul {SJ} = Su2 = xy (3.10) Su3 Es\ {Es}e = Es2 Es3 in which Su,, is the u(u = x, y) component of the surface vector basis functions along the ith edge. On substituting (3.9) into the surface integrals of (3.1) and taking the first variation of Fs with respect to {Es}I, we obtain sFs = E {[B]e{Es}e + {L}e} (3.11) e where [B]e - 2ko2[Se][Se] + 2 [{S} - {S} [ }T,{SiT] } *Go(r, r') dS dS' (3.12) and {L}e - j2koZo [] ( dS (3.13)

49 Note that in (3.12) the elements of the array [Se] are functions of the observation vector r, whereas the elements of [Se]T are with respect to the integration point r'. A suitable set of linear edge-based surface basis functions is 2 A z x (r - ri)c(r) r C Se Si(r)= 2Ae (3.1L4) 0 otherwise In this expression (referring to fig. 3.3), li denotes the length of the ith edge and r, is ith edge +x Figure 3.3: Pair of triangles sharing the ith edge the position vector of the vertex opposite to the ith edge. Since each edge shares two triangles, one is defined as the plus and the other as the minus triangle. Therefore, c(r) is given by 1 r E S+ 1(r) S (3.15) [ -1 rC 5 -where Se = S+ + 57. The constant A, in (3.14) denotes the area of the plus or minus triangle depending on whether r C S+ or r C Se-. We note that S-(r) x z yields the basis functions used by Rao, et al. [32] in their moment method solution of boundary integral equations. The explicit expressions for the boundary integral equation subsystem is given in Appendix B for reference.

50 3.2.3 Combined FE-BI System To construct the final system for the solution of the electric field components we combine (3.5) and (3.11), and after assembly we obtain the system {[A]{E} + {K}} + {[B]{Es} + {L}} = 0 (3.16) In this, {Ki} and {L} are the excitation vectors due to the interior current sources and the exterior excitation, respectively. The unknown electric field vector {E} consists of all field expansion coefficients with respect to the element edges except those coinciding with perfectly electrically conducting (PEC) walls, PEC antenna element(s) or PEC pins inside the cavity. Finally, the vector {Es} represents the unknown surface fields whose entries are part of those in {E} with their corresponding edges on the aperture. The explicit expressions for the matrices and vectors in (3.16) can be readily extracted from (3.6), (3.7) and (3.12) (see also [33]). It is evident that [A] and [B] are symmetric as a result of the assumed isotropic medium and reciprocity. In addition, [A] exhibits high sparsity due to the FEM formulation whereas [B] is fully populated. Two approaches may be followed in carrying out the solution of the combined subsystems when an iterative solver is employed such as the biconjugate gradient (BiCG) method [34]. These approaches differ in the manner used for the evaluation of matrix-vector products called for in the iteration steps. One could sum the coefficient matrices [A] and [13] by adding up the corresponding matrix entries prior to the execution of the BiCG algorithm, or instead the resulting vectors may be summed after carrying out the individual matrix-vector products. We observe that the first approach is more efficient in terms of computation time after reordering the combined matrix and storing only the non-zero elements. This is because, in the context of

51 this scheme, the combination of the two matrices is performed only once outside the iteration. However, the second approach is compatible with the BiCG-FFT scheme, where the FFT algorithm is employed to exploit the convolutional property of the integral operator, thus eliminating a need to explicitly store the entire BI matrix. The BiCG-FFT technique will be discussed in chapter 4. 3.3 Numerical Implementation Based on the above presented FE-BI formulation, the hybrid method was implemented and a computer program was developed for the analysis of radiation and scattering by cavity-backed patch antennas of arbitrary shape. The antenna geometry/mesh is first generated as shown in fig. 3.4 and supplied to this program in an Figure 3.4: A typical geometry/mesh for a cavitv-backed circular patch antenna. input file which, as a minimum, contains (a) the nodes and their (x, y, z) coordinates; (b) the tetrahedral elements and the corresponding nodes forming each element; (c) the nodes identifying the cavity aperture;

52 (d) the nodes identifying metallic boundaries, antenna patch(es), feed(s), and possible vertical posts. For arbitrary antenna geometries, it is necessary to employ a sophisticated volume mesh generation package and a number of these are available commercially. Typically each of these packages generates a "universal file" that can be readily preprocessed to extract the aforementioned input list. Given the above list of data, an interpretation routine is used to convert the information from node-elements to edge-elements. We usually refer this procedure as data preprocessing. The flow chart shown in fig. 3.5 describes the major implementation procedures from the mesh generation, a few data preprocessors, the FE-BI kernel, to the BiCG solution and finally the data output. One of the primary complications in the hybrid technique implementation is the efficient treatment of combining the two separate subsystems. It is noted that the FEM sparse matrix is large in dimension but requires less storage space, while the boundary integral subsystem is always small in size but can be dominant in terms of memory demand. This is particularly true when the non-metallic portion on the cavity aperture predominates over the antenna radiating elements. Furthermore, the boundary integral subsystem in a general purpose hybrid FE-BI implementation is entirely independent from that of the FEM and even the basis functions can be independently developed. This also accounts for the two arbitrary numbering systems and combining them is a relatively complex task. One major advantage of the method however is that these two subsystems can be developed and validated individually. Once both of the subsystems are verified, the coupling of the subsystems is accomplished by enforcing the boundary conditions implicitly on tangential H fields via the integral representation and explicitly on tangential E fields over the interior and

53 IMPLEMENTATION FLOW CHART f FOR FEM System BI System Feed Models Combinations Data File Set FEM System r --- —-- --— 4 ------- BI System Feed Models Combinations / BiCG Solver / OUTPUT Figure 3.5: A flow chart describes the major implementation procedures from the mesh generation, a few data preprocessors, the FE-BI kernel, to the BiCG solution and finally the data output.

54 exterior regions. A significant effort was devoted to developing a program in such a manner so that both the storage and computational requirements can be minimizedl. Specifically, the boundary conditions on the metallic surfaces are enforced a priori to condense the system which involves only nonzero field components. To this point1, the sparse finite element matrix was stored as a single array of length NpNnz, where Ni is the total number of unknowns within the cavity volume and Nnz denotes the maximum number of nonzero row entries. The BI matrix was stored in different ways for the evaluation of the matrix-vector products. If the BiCG solution was to be carried out without special treatments (such as incorporating the FFT), then the N's x Ns BI integral matrix is added to the FE array resulting in a 1-D array about NAnz + N2 long. For slot antennas, including cavity-backed spirals, and moderately sized systems, it was found preferable to use this scheme. In that case the generation of a single combined FE-BI matrix before execution of the BiCG algorithm reduces the computational requirements. This is because a number of operations associated with the repeated combinations of the FE and BI subsystems within the BiCG iteration is avoided. The alternative is to carry out the matrix-vector products for the FEM and 13I subsystems separately. The advantage of the scheme becomes apparent when a special treatment is performed on to the numerical system for efficiency consideration and this will be investigated at certain depth in chapter 5. 3.4 Selected Numerical Results We present below some representative numerical results for the purpose of validating and demonstrating the robustness of the tetrahedral formulation for scattering and radiation by different configurations of cavity-backed antennas. In each case the computed results via the FE-BI method are compared with reference measured or

55 calculated data. Scattering and radiation by a circular patch: Fig. 3.6 illustrates a circular patch residing on the surface of a 0.406 cm thick substrate having a relative dielectric constant of cr = 2.9. The patch's diameter is 2.6 cm and the substrate is enclosed in a, circular cavity 6.292 cm wide. This cavity and the patch were recessed in a low cross section body for measuring its RCS. A comparison of the measured and calculated backscatter aso RCS as a function of frequency is also shown in fig. 3.6. For this computation the direction of the incident plane wave was 60~ from normal, and as seen the agreement between measurements and calculations is very good throughout the 4-9 GHz band. Input impedance measurements and calculations for the same patch are displayed in fig. 3.7. The probe feed in this case was placed 0.8 cm from the patch's center and it is again seen that the measurements and calculations are in good agreement. Radiation by a one-arm conical spiral: We considered the modeling of this radiator to demonstrate the geometrical versatility of the FE-BI method. Two projections of the spiral radiator and surface mesh are illustrated in fig. 3.8. The top and bottom edges of the strip forming the spiral follow the lines p = 0.0503A exp[0.221(~+ 2.66)], z = a+ exp(0.2210), where (p,, z) denote the standard cylindrical coordinates, a~ are equal to 0.0832A and 0.0257A, respectively, and 0 < b < 27T. This spiral arm resides on an inverted cone (9.24 cm tall) whose bottom cross section has a diameter of 1.68 cm and the top cross section has a diameter of 21.78 cm. For our calculations A = 30 cm (f = 1 GHz) and the spiral was situated in a circular cavity 10.01 cm deep. The computed Ey principal plane radiation pattern taken in the ( = 900-plane, using a probe feed at the cavity base, is given in fig. 3.9. It is seen that this pattern is in good agreement with the

56 U.UU) I —, U 00 I-C CO 00 0 CO -10.0 -20.0 -30.0 -40.0 -50.0 -60.0 ~4.00 5.00 6.00 7.00 8.00 frequency (GHz) Figure 3.6: Comparison of the computed and measured u(o0 backscatter RCS as a function of frequency for the shown circular patch. The incidence angle was 30~ off the ground plane. 1.0 0 1.0 Figure 3.7: Comparison of the computed and measured input impedance for the circular patch shown in fig. 3.6. The feed was placed 0.8 cm from the center of the patch and the frequency was swept from 3 to 3.8 GHz.

57 data given in [35]. As can be expected, the Ee pattern (not shown) differed from the measured data near the horizon because of interference from the finite circular cavity housing the spiral which was included in the analytical model. The latter was not part of the measurement configuration which consisted of the spiral antenna on a large circular plate. Annular slot impedance: Fig. 3.10 shows a narrow circular (0.75 cm wide) annular slot situated in a circular cavity 24.7 cm wide and 3 cm deep. Because the annular slot is narrow, the implementation of the BI subsystem is very small for this application and as a result there is no need to invoke the FFT in the BiCG algorithm. The FE-BI method is basically quite effective in modeling small aperture configurations without a need for special computational considerations. Input impedance calculations as a function of frequency for this radiator, excited by a probe placed across the slot, are shown in fig. 3.10, and agree well with the values calculated via a modal-boundary integral method [14]. For these calculations, the frequency was swept from 700-1000 MHz. The dielectric constant of the material filling the cavity was set to cr = 1.35 as in [36] and this is an effective value to account for the presence of a dielectric slot cover used as part of the measurement model for holding the plate. Stacked circular patch antenna: To derlonstrate the capability of the developed hybrid technique, we now present a qualitative study and visualization of the near field distribution inside a cavitybacked, stacked circular patch antenna as shown in fig. 3.11. Note that the similar configuration with rectangular patch shape has been investigated and found a significant bandwidth increase. This is because of the dual frequency resonance due to the two stacked patches. The circular patches are more attractive than stacked

58 Figure 3.8: Illustration of the configuration and mesh of the one-arm conical spiral used for the computation of fig. 3.9. I I I I ~ 7. i. - -. I L.. I... II.. r. I l I I IV T I I I II I I V T I I I I I I I T, II I V T I -1 I. I I i I 0..5 E 4 a> 0.,4.5 md d -10. -20. -30. -40. -50. E — (FE-BI) a Eo measured —.......-.....I...... I..I................I -90. -60. -30. 0. 30. 60. 90. 0 Figure 3.9: Comparison of the calculated radiation pattern (E4), taken in the ( -= 90~-plane, with data in reference [35] for the one-arm conical spiral shown in fig. 3.8

59 'o 0.0 1 10 f.4 0 f I I 4 J i I b a = 12.35 cm b = 0.75 cm po= 7.7cm 0.7 < f < GHz d =3 cm i 1.0 Figure 3.10: Comparison of input impedance calculations for the illustrated cavity backed slot.

60 rectangular patches because they occupy a small area when operated at the same frequency [37]. Unfortunately, no sufficient research on this geometry has been reported in the literature due to a lack of analysis tool. It turns out that the above presented hybrid FEM technique is well suited for this study. To show this, we chose a cavity-backed, stacked circular patch antenna fed with an offset single vertical post underneath the lower patch to link the cavity base (viewed as a ground plane). No direct electric contact between the upper and the lower patches exists and thus the power transfer has to rely completely on the electromagnetic coupling from the lower to the upper patch. This can be clearly verified from the near field visualization, which is available and complete only via the PDE related techniques. (The laboratory measurement may provide the image of the field distribution above and near a microstrip [38].) Another interesting point is that though the antenna was fed with a single offset probe, the energy is concentrated at two distinct regions. One is around the probe feed, and the other is near the opposite location with respect to the center of the patch. The two regions act as out-of-phase electric pole to effectively excite the antenna. Although the patches are circular in shape, the offset excitation ensures the linear polarization in radiated fields. Figure 3.11: Visualization of the near field distribution at the lower layer of a stacked circular patch antenna.

CHAPTER IV Efficient Boundary Integral Subsystem- I 4.1 Introduction As is known, the hybrid finite element-boundary integral method is accurate and capable of handling a variety of conformal antennas. However, the drawback associated with this technique and any other global truncation approaches can make it less attractive. This is especially true if one is interested in modeling large antenna systems (arrays). Although the FE-BI method is particularly suited for the configurations with relatively small aperture size and possibly complex cavity design (including feedlines, isotropy/anisotropy, other layers of metals, etc.), it would be much more useful to accelerate the speed and reduce the CPU requirement for the hybrid approach. One possible solution is the CG-FFT technique discussed in this chapter. The boundary integral (BI) equation subsystem leads to a fully populated matrix whose size is determined by the number of aperture mesh edges. For large apertures, this analysis becomes impractical in terms of storage and computation time requirements, and to overcome this inefficiency, a uniform zoning of the aperture is required. By resorting to the structured mesh, the boundary integral matrix can be cast into a discrete convolutional form, thus permitting the computation of the matrix-vector 61

62 products via the discrete Fourier transform (DFT) and avoiding a need to store tie full BI matrix. This memory saving scheme has already been applied to IE solutions involving rectangular grids [39,34], and a similar implementation was also reported for triangular surface grids [40] involving inherent approximations. In this chapter, we first show that the BiCG-FFT solver can be precisely implemented on uniform triangular meshes. The differences between the rectangular and triangular meshes are also described. For non-rectangular antenna geometries, a special treatment referred to as the overlaying scheme is proposed and discussed in section 4.3. A few results are presented which demonstrate the method's validity. 4.2 Application of Conjugate Gradient Algorithms The Conjugate Gradient (CG) iterative solution of linear systems of equations has been extensively investigated and the representative references are collected in [41]. Although the state-of-the-art CG algorithms do not lend themselves as a robust input/output "black-box" [42], they are indeed capable of handling large-scale comnputational problems which may be impossible for direct system solvers. It is especially desirable to employ the algorithms when one seeks the solutions of large —scJle numerical systems without resorting to costly computing resources. 4.2.1 BiCG Algorithm With Preconditioning Conjugate gradient (CG) algorithms have been developed for over forty years [43,44] and one of the primary applications nowadays is to solve large scale linear systems, as aforementioned. It is noteworthy that there exist various versions of the CG algorithms, taking advantage of different properties of the matrix such as symmetric and sparsity. Also, preconditioning is often used to speed up convergence.

63 As suggested in [45,46], the algorithm used in this work is as follows: Given p= ri = b - A x Pi = r For k = 1,2,3,... rk ~ rk Ck - k — __ Pk* A Pk rk+1 = rk - CkA * pk -* dk+l = rk - =A Pk rk * rk Pk+1 = rk + 3kpk Xk+ = Xk + c(kPk where * denotes the complex conjugate and T is the transpose. This version of the iterative algorithm is quite general in terms of the system matrix to be solved and is usually referred to as the Biconjugate gradient (BiCG) method for unsymmetric systems. If the matrix A is symmetric and the initial value is chosen as ri = r7, the algorithm can then be shortened to require only one matrix-vector product per iteration, since in each step rk and Pk are complex conjugate of rk and Pk, respectively. The ordinary conjugate gradient algorithm can be considered as a special case of the BiCG when A is Hermitian (i.e. A = A*T). Again in this case, the algorithm can be shortened to have about 50% less computational effort. The CG algorithm is also amenable to a straightforward interpretation of its convergence principle. Basically,

64 the algorithm minimizes the function f(x) = -x A x- b x Hence x obtained from the CG algorithm after the k iteration steps becomes the solution of the equation Af(x)= A x - b = 0 The CG type of algorithms did not become competitive until preconditioning was introduced to improve the original system condition and significantly reduce the convergence rate. One simple preconditioner is the inverse of the original diagonal matrix. In our applications this preconditioning has been quite successful and is used in conjunction with the BiCG algorithm. 4.2.2 BiCG-FFT Algorithm For Linear System In our work, the linear system of equations is usually large in size, partially full and partially sparse. The conjugate gradient (CG) type of algorithms can be used to alleviate the memory requirement since the LU decomposition requires excessive memory and( CPU time. However, the partially dense matrix due to the boundary integral equation may still dominate the CPU demand. This is because the method of moment (MoM) always leads to a dense system by its nature. Solving the dense system in a traditional manner requires O(N2) order of operations per iteration, where N is the boundary integral system dimension. Reduction of the operation counts will of course significantly decrease the solution time and this can be accoinplished by recasting the BI system onto a few Toeplitz submatrices and making use of the fast Fourier transform (FFT) to carry out the matrix-vector products in the BiCG solver.

65 As described in the next section, the boundary integral equation can be cast into a convolutional form if a uniform grid is applied for discretization. This is not surprising since the Green's function is involved in the integration. To solve the equation using the CG algorithm, it remains to carry out the convolution at each iteration. To this end, one may calculate the convolution by taking the Fourier Transform of two spatial data sequences (arrays) in which the convolution becomes the product of the two 'frequency' sequences. An inverse transform of the product yields the result. In contrast to the order of 0(N2) CPU requirement for a matrix-vector product in a traditional fashion, the scheme needs 0(Nlog2N) operation counts if the FFT algorithm is employed. The operation reduction is indeed significant and the technique has the lowest CPU demand among integral equation solvers, including the fast multipole method (FMM) [47], thus is always preferred. 4.2.3 Convolutional Form of Boundary Integral The boundary integral equation is discretized using the structured triangular grid, and the relation between the unstructured and structured mesh is described in the next section. We recognize that the triangular grid consists of equal right triangles as shown in fig. 4.1 and thus involves three different classes of edges (class 1, 2 and 3). These include the x-directed, y-directed and the diagonal edges, all of which are uniformly spaced. For the FFT implementation each class of edges is independently numbered in accordance with their geometric location. Specifically, the ith class will carry the numbering (rn, n) if the edge is the mth along the x direction and the nth along the y direction. The indices (m, n) take the values in = 0,1,2,...,M n =0,1,2,...N

66 z N 0 - m: 0 1 2 3........ M Figure 4.1: Structured mesh consists of equal right triangles with = 1 for the y-directed edges, i = 2 for the diagonal edges and i = 3 for the x-directed edges. Consequently, we find that m I = k M-2 i=1 M-1 i=2 M-1 = 3 i N - k N-1 i 1 N-1 i=2 N-2 i =3 (4.1) where M and N denote the numbers of elements along the x and y directions, respectively. To perform the integrations for the evaluation of the boundary integral matrix elements, it is now convenient to rewrite the basis functions (3.14) in terms of the new indices (m, n). We readily find that the edge-based basis functions associated

67 with each of the aforementioned class of edges can be rewritten as (nAy - y)x + (x - mAx)y (x,y) 6 S+ Smn(X, y) (Y - (n + 1)y)x + ((m + 2)Ax - x)y (x, y) E S 0 otherwise (4.2) SMn(X, y) Ax)2 + (Ay AwAy I ((n + 2)Ay - SM(x.y) 1= ( (- nAy)x I4 0 (nAy-y)x + (x-(m + l)Ax)y c S+ (y - (n + l)Ay)x + (mAx - x)y C S[ (4.3) 0 otherwise y)x + (x - (m + l)Ax)y (x, y) S+ -(mAx - x)y (x,y) S,- (4.4) otherwise where the superscripts refer to the edge class. After the discretization and assembly processes, one obtains a discretized version of the BI system, from which each entry of the boundary matrix-vector product can now be calculated as 3 MP NJ {EBI subsystem} = [B] {E )= > > B, mm,n m, (4.5) j=-1 m'=0 n'=O in which (rn, n) are the geometric location indices for the ith class observation edges whereas (m', n') are the same for the jth class edges belonging to integration elements. Thus, the specification of the indices i, m and n completely defines the entry ki =n: niM + M of the column resulting after the execution of the bounda]^y matrix-vector product. It is readily found that Bmn mI = -2-o f S. Smm Go(r, r') dx dy dx' dy' (4.6) + (A y)2 (r) I/(r) lil Go(r, r') dx dy dx' dy' (AXAy)2

68 with Ay i=1 = /(Ax)2 + (Ay)2 i =2 (4.) Ax i 3 More importantly, it can be shown that the BI subsystem (4.6) exhibits the convolutional property B",, = B(m,, _-n) and thus we can rewrite (4.5) as 3 [B]{Es}= BJ * E3 (4.8) j=1 where the * denotes convolution. The proof will be presented in the end of this section to ensure the smooth discussion of the remaining procedure. It is now seen that the computation of the boundary matrix-vector product can be performed by employing the 2-D discrete Fourier transform (DFT), thus avoiding a need to store the BI matrix other than those entries which are unique. When the symmetry property of B(_-,, _,-n) is also invoked, implying B _ 3 - (4.9) (m-m',n-n') (m'-m, n '-n) it is concluded that the total non-redundant entries in the BI matrix are 3 3 Np = Z Z N(M' + M - 1) (4.10) i=1 j=l This should be compared to the (Z3_ M'iNi) entries whose storage would normally be required if the BI system was not cast in convolutional form and it is noteworthy that Np is much smaller in number. To avoid aliasing, it is necessary that B(-m B (mFi, ai) bbe cast in a 2-D array which has the usual periodic form, and zero padding may also be required to make use of the standard FFT routines. Specifically, the matrix-vector product (4.5) is executed by using the MFTxNF'T

69 array K0 < in < Mt 0,< i < N' - MFT-Mt + 1 < m< MFT B"J(m:n)=<,i- -NFT) ~ < mf < Mi (4) ' ' NFT- NJ + 1 < n < NFT Bi( -- 1 - MFT, n - 1- NFT), MFT- MI + I < < MFT 0 otherwise with the corresponding field vector given by E( Effimn) 0 < < M, 0 < n < NJ Ep(Qm, n) = (4.12) 0 otherwise and MFT and NFT must be powers of 2 if a radix 2 FFT algorithm is used. In the BiCG-FFT algorithm the BI subsystem vector is symbolically computed as 3 {BI subsystem} = S {DFT-1 {DFT{BP'} DFT{E }}} (4.13) i=1 The presence of the operator S indicates the necessary reordering of the 2-D array which results after the inverse FFT operation into a single column with the proper indexing for addition to the FEM subsystem. It should be remarked that in contrast to [40] the integrals (4.6) are evaluated without introducing any approximation. This is necessary to preserve the symmetry feature of the global combined system. As promised, we now show (4.8), or the relation B",,, = BJf,,_, to conclude this section. To simplify the proof, we refer to fig. 4.2 and consider only the first integral in (4.6). The same proof can be appied to the second integral in (4.6). In addition, with no loss of generality for the proof, the i = 1 class edges

70 (m,n) ' (m',n') Figure4.2: Illustration of two triangles with the corresponding indices to help to prove the convolutional property of the boundary integral. (y-directed) in S+ for the trial function and the i = 2 class edges (diagonal) in Sfor the testing function will be used. To this point, substitutng Sl,(x, y') in S+ and ~2 n,(x: y) in Se into the first term in (4.6) yields Int nmnl J C[(n - -y)(n - y) + (,x - mx)(m/X + 2 - ')] Go [(x - x')x + (y - y')y] dxdy dx'dy' (4.14) where C is a constant coefficient and its detailed form is not of our concern for this proof. Note that the integration limits should be set as [mAx, (m + l )Ax] and [nAy, (n + ),Ay] for the unprimed coordinates and similarly for the primed one,. Therefore, (4.14) will be simplified if the following transforms are introduced, viz. x = max + x' = mAx + ~' (4..,) y = nay + rI yq = nAy + ql' Indeed, on substituting the transforms, one obtains r rpn R r raxny Intmnml = `C / / // - (A') Go -mx[(m - m')Ax+ (G - ')] + [( - n') xy + ( - ')]} d-dr dd'dq' (4.1C) wher C s a onsantcoefcint nd t ealdfr snto u ocr o thisprof. Nte hat he ntegatinliissodbeeta mx( +1A]an [n~, ( + )Ay fo th unrim cordnatsadsmlryfrtepie ns

71 Although this integral may not be expressible in terms of an elementary function, one is however assured that the resulting expression must be a function of [(m - nm'), (n - n')] no matter what form the Green's function Go will take. Mathematically, it means Intm,m'n' = Intm-m',n-n' (4.17) This is the desired relation. From the proof, we can conclude that the convolution in (4.8) holds. 4.3 Mesh Overlay Scheme As described above, the BiCG-FFT solver requires uniform aperture griddi:ng so that the BI subsystem can be put in block circulant form. This can be always achieved during mesh generation whenever the patches are rectangular in shape or in case of radiators which are placed at some distance (usually small) below the aperture as shown in fig. 4.3(a). In the latter case, one might need to add an appropriate absorbing material around the edge/corner of the cavity near the aperture to avoid possible edge/corner effects, especially when the aperture size is not large enough. Fig. 4.3(b) shows the example of this implementation. 4.3.1 Field Transformations However, for circular, triangular, or other non-rectangular patches on the aperture, it is not natural to construct a uniform mesh using the mesh generator. Typically, the aperture mesh is necessary to conform to the patch shape, leading to an unstructured free surface grid. In this case, to make use of the efficient, low memory BiCG-FFT algorithm, an approach is proposed to overlay on the unstructured aperture grid another coincident structured grid as shown in fig. 4.4. The boundary integral subsystem is then constructed by using the overlaid uniform grid whose edge

72 uniform grid -.4 / recessed patch recessed patch cavity (a) Circular Patch Recessed Underneath the Ground Plane 3 0. F-,.,, II,,, I,,..,,, I- I...I,,,,,,,,,,,.,,,,,,,,....,,,.,,,,,,,,,,.,,,,,,,,1,~.. o,cl c< 1-,-... r.j (de) c..-. A4 c5 ~...C4 PQ 20. 10. 0. -10. -20. -30. -40. -50. -VI I -, permittivity = 4 -------- permittivity = 1 Q permittivity = 4 (with FF1T) v permittivity = 1 (with FFT) '71, i IJ -60.......... I......... I,,, 1,... I- 1, I....I....I.... I......... I- I...... I........ I I....... U, I LL 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 0 (b) Figure 4.3: Printed circular patch antenna is modeled using the recessed scheme to incorporate the BiCG-FFT algorithm. (a) Illustration of the configuration; (b) Comparisons of the BiCG-FFT result with that of the ordinary FE-BI technique presented in chapter 3.

73 3 2,: 0 2 3 Figure 4.4: Overlay of a structured triangular aperture mesh over an unstructured mesh, shown here to conform to a circular patch. fields can be related to those on the unstructured grid via two sparse transformation -matrices. Trhat is, it is necessary to append to the system (3.11) the relations OEs}U = [TF]{ES}nu {Es,} = [TB]EU (4.18) where the subscripts u and nu refer to the field coefficients of the uniform and non-uniform aperture grids, respectively. Also, [TF] and [TB] refer to the forward and backward transformation matrices, respectively, with Nu and Nnu denoting the numbers of the uniform and non-uniform mesh edges on the cavity aperture. To derive the elements of [TF], we begin with the expansion (3.14) and enforce it at three points on each edge belonging to the uniform grid. We conveniently place these three points at the center and ends of the edge (see fig. 4.5). Given the fields at these points, we can interpolate the field along the (m, n) edge

74 I s. I, overlaid unstructured Y r e mesh triangles rend2 rmid- (m,n) uniform grid edge x \ / rendl / \ / Figure 4.5: Illustration of the parameters and geometry used in constructing the transformation matrix elements between the structured and unstructured mesh. of the uniform grid using the weighted average f INend 1 2(E 2 \ end, E Enu(rendi) Nmid + - Emid E$gu(rmid) mid k=l 1 Nend2 + end2 E (rend2) (4.19) 2Nend2 k-i in which &. denotes the unit vector along x, y or the diagonal, depending on the class of edge being considered. The quantities E k represent the fields in the non-uniform grid triangles with the superscript k being a sum variable in case rendi, rend2 or rmid specify a point shared by more than one triangle. Obviously, NAndi, Nmid and Nen,12 denote the number of non-uniform grid triangles sharing the node at rendi, rmid and rend2, respectively, and will typically be equal to unity. After assembling (4.19) into (4.18) we find that the elements of the forward

75 transformation matrix are given by I Nendl 3 k=e 1 =1 Nmid 3 + z 1 ij(S (rmid) mk=1 =1 Nend2 3 ( + d E 6ijS(rend) (4.20) end2 k=1 =1 in which | 1 j = je Hijt= - { 0 otherwise and the global indices i and j correspond to the ith uniform grid edge and the jth nonuniform grid edge. The subscript ja is the global index used in numbering the nonuniform grid edges, whereas the subscript f (= 1, 2 or 3) is the local edge index used in the definition of the basis functions So. We remark that the explicit computation and storage of the transformation matrix elements results in a substantial increase in efficiency because it avoids the usual assembly process during each iteration step and that the proposed overlay scheme allows the analysis of large non-rectangular patc: arrays because storage of fully populated BI system matrix is avoided. The user needs to only provide an additional data file which flags the uniform grid edges lying on a PEC element and this is an important user-oriented feature of the formulatior. Following the same procedure, we can obtain the expression for the entries of the backward transformation matrix. It should be noted that assuming each uniform grid edge traverses three or less non-uniform grid triangles, the non-zero entries in each row of [TF] will be 9 or less. However, they can reach a maximum of 18 if the midpoint and endpoints reside on an edge of the non-uniform grid. The rnaxinnur non-zero entries in each row of [TB] will be 15, but the typical number is much less.

76 4.4 Results Figure 4.6 shows a cavity-backed 2 x 2 patch array, where each patch is a rightangled triangle. Since this geometry is adaptable to a uniform mesh with rightangled triangles, it is used to verify our proposed FE —BI technique incorporated with the BiCG-FFT system solver. The developed FE-BI code with the BiCG-FFT is first compared with the original version of the hybrid FE-BI technique described in chapter 3. As shown in fig. 4.7, the monostatic radar cross section (RCS) pattern over the space 0 < 0 < 90~ at the 0 = 0 plane agrees very well with that computed using the regular BiCG solver. It is also informative to compare the scattered fields by the same cavity-backed structure without the patch array to find the contribution of the array to scattering. Figure 4.8 shows the monostatic RCS patterns by the aperture with the absence of the patch array. Again, the computations were obtained using both the regular BiCG and BiCG —FFT versions of the FE-BI methods. As can be seen, the level of the scattered field at the normal incidence reaches above zero in dB/A2 with the presence of the patch array, whereas the scattering by the aperture at the same incidence is about 23 dB/A2 lower. To varify the overlaying scheme for nonrectangular geometry, we evaluated a bistatic RCS' scattering as shown in fig. 4.9 by a cavity-backed circular patch antenna. In this case, the dielectric fillings of or = 4 and cr = 10 inside the cavity were used, respectively, and the results obtained using the regular BiCG FE-BI method are compared with those computed using the BiCG-FFT with the overlaying transformations. It is observed that the agreement is quite satisfactory in scattering analysis.

77 Patch Array _ _ Figure 4.6: Illustration of the cavity-backed triangular patch array. For radiation analysis (e.g. input impedance) where an accurate prediction of the near-zone fields may be required, the accuracy of the overlaying scheme can be significantly enhanced by considering the trial-testing element's interactions in the boundary integral. Specifically, it is suggested to separate the interactions between the closed —region elements from the far zone weak couplings. The strong closeregion couplings are treated using the normal method of moments, whereas the weak coulings are computed using the fast algorithm. This approach has been reported (see e.g. [47]), and once combined with the overlaying scheme, it can be used to control the accuracy of the FE-BI technique in an adaptive manner.

78 10. 0. lo0. -20. -30. -40. -50. -60. -70. -80. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. ft Figure 4.7: Comparisons of the monostatic radar cross section scattering by a 2 x 2 triangular patch array shown in fig. 4.6. The reults were computed using the regular BiCG FE-BI technique described in chapter 3 and using the BiCG-FFT proposed in this chapter. CC1 ec ru x -10.0 -20.0 -30.0 -40.0 -50.0 -60.0 -70.0 -8)00 t0 ((( I).0() 20.00) 30.00 40.00 50.(X 60.00 70.(00 80().0X t,).0 e Figure 4.8: Comparisons of the monostatic radar cross section scattering by an empty aperture with the same cavity size and dielectric filling (e=l1) as the structure shown in fig. 4.6. Again, the results were calculated using the regular BiCG FE-BI technique described in chapter 3 and using the BiCG-FFT proposed in this chapter.

79 RCS Comparison of Transform With No Transform -8. - -16. -24. - u 4 -CZ CL4 C/ u C4.2 cl 4-1 C'n.... 11 V V V V i -- '' _- _-_ - 6), A --- A " v _ w v. -32. -40. -48. -56. -64. -72. -80. With No Transform (permittivity=l 0) -------- With No Transform (permittivity=4) A With Transform (permittivity=10) v With Transform (permittivity=4) I. I. I I..- 11 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. e Figure 4.9:: B3istatic RCS scattering by a crcular patch antenna modeled using the regular BiCG FE-BI and the BiCG-FFT algorithm with overlaying transform.

CHAPTER V Efficient Finite Element Subsystem II 5.1 Introduction As demonstrated in the previous two chapters (also reported in [13, 29, 48]), a hybrid finite element-boundary integral technique [48,13] can be employed for characterizing conformal antennas of arbitrary shape [29]. Indeed, planar/non-planar, rectangular/non-rectangular designs, ring slot or spiral slot antennas with probe, coax cable or microstrip line feeds can be simulated with this formulation. This is because of the geometrical adaptability of tetrahedral elements used for the implementation. However, in practice, certain configurations require extremely high sampling rates due to the presence of fine geometrical details. Among them are a variety of slot antennas (spirals, rings, slot spirals, cross slots, log-periodic slots, etc.), where the slot width is much smaller than the other dimensions (cavity diameter or inter-distance of slots). In these cases, the mesh is extremely dense (with over 50, 100 or even higher samples per wavelength), whereas typical discretizations involve only 10-20 elements per wavelength. This dense sampling rate is especially severe for 3-D tetrahedral meshes, where the geometrical details usually distort the tetrahedrals. The numerical system assembled from this type of mesh often leads to a large system condition due to the degraded mesh quality. Also, mesh generation is 80

81 Ground Plane Thin Slot V Sf Figure 5.1: Geometry of cavity-backed microstrip antennas tedious and the solution CPU time is unacceptably large. In this chapter, we propose a finite element-boundary integral formulation using edge-based triangular prism elements. It can be shown that this element choice is ideally suited for planar antenna configurations containing spirals, circular and triangular slots.. Among many advantages of the prismatic elements, the most important is the simplicity of mesh generation. Also, much smaller number of unknowns is required for an accurate and efficient modeling of complex geometries. Below, we begin by first outlining the finite element-boundary integral (FE-BI) formulation for slot antenna modeling. A new, physically meaningful, set of edge-based functions for prisms is then presented to generate the discrete system of equations. Finally, the applications of the proposed hybrid FE-BI method to various antenna radiation and scattering problems are given to conclude the chapter. 5.2 Hybrid FE-BI Formulation Consider the cavity —backed slot antenna shown in fig. 5.1 where the cavity 's recessed in a ground plane. To solve for the E-field inside and on the aperture of the

82 cavity, a standard approach is to extremize the functional (3.1) which, for radiation and scattering problems, may be generalized as F(E) - I I {(V x E): 1, ((V x E) -ko2E. -E}dV + fjI E (jkoZoJi + V x Mi dV + jkoZo E - (H x )dS (5.1) o+Sf where Or and i7,. denote the relative tensor constitutive parameters of the cavity medium, Zo and ko are the free space impedance and propagation constant, respectively, So represents the aperture (or slots) excluding the metallic portions and Sf denotes the junction opening to the guided feeding structures. Also, V, is the volume occupied by the source(s) (Ji or Mi) and H is the corresponding magnetic field on So and Sf whose outer normal is given by n. As before, the explicit knowledge of H in (5.1) is required over the surface So and Sf (also referred to as mesh truncation surfaces) for a unique solution of E. Specifically, the magnetic field H over So may be replaced in terms of E via a boundary integral (BI) or absorbing boundary condition (ABC), whereas H on Sf is determined on the basis of the given feeding structure. This version of the functional as compared to (3.1) allows an easy inclusion of various feed models, such as aperture coupled slot, coax cable, etc. (see chapter 6 for details). In this chapter since we concentrate on improving the FEM efficiency, therefore the boundary integral method will be employed as in chapter 3. It will be seen that this implementation indeed meets the accuracy need without extra CPU burden. In the context of the FE-BI, H is represented by the integral H = H~0 + 2jkoYo I Go(r, r') (^ x E(r')) ds' (5.2) where G is the electric dyadic Green's function of the first kind such that n x G = 0 is satisfied on the (planar, spherical or cylindrical) metallic platform (refer to chapter 1).

83 For the antenna problem shown in fig. 5.1 where the platform is a planar ground plane, G becomes the half space dyadic Green's function /- 1 e-jko Ir-r'I G= I+ - VV --- (5.3) 47flr - rj'[ 0 with r and r' being the observation and integration points, respectively, and I Kxx + Zy + $Z$ is the unit dyad. In connection with our problem, i.e. that of a cavity recessed in a ground plane, Hg0 is equal to the sum of the incident and reflected fields for scattering computations, or zero for antenna parameter evaluations. To discretize the functional (5.1), we choose to subdivide the volume region using prismatic elements as shown in fig. 5.2 and fig. 5.3, The field in each of the prisms can be approximated using the linear edge-based expansion [49-51] 9 E = EJV [V]{Ee}, (5.4) j=1 where [V] =- [{l, {}, {V}, {V}], and {Eel {E E.... E }'T. The vectors {Vu}, u x, y, z, are of dimension m = 9 and they simply represent the x, y, z components of V, associated with the jth edge of the eth element. Since V6 are chosen to be edgebased functions, the unknown coefficients ES represent the average field along the jth edge of the eth element. A corresponding representation for the aperture fields is 3 E(r) = Es Ss(r) [S]{Es}, (5.5) i=1 where [S]s =- [Sx, SJ], and V(r) reduces to SS(r) when the position vector is on the slot. To generate the discrete system for ES, (5.4) and (5.5) are substituted into (5.1) and subsequently F(E) is differentiated with respect to each unknown E6. With the

84 6 Figure 5.2: Illustration of tessellation using prisms A z 21 Z=Zc 6 i=1,2,3 j=4,5,6 k=7,8,9 Figure 5.3: Right angled prism

85 understanding that the surface field coefficients EJ are a subset of EJ we obtain N, N, N, N, { HE = [Ae]{E + [BS]{Es} + {I } + L} = (5.6) e=l e=l esl where the sums are over the total number of volume or surface elements. In this, the matrix elements are given by A = {V x. V x V - Vj-koV V} dv (5.7) I = jkoZo I Vi. [i+ V x 1. V x Mi dv (5.8) Bs = - JJj 2kS(r)- S'(r')Go(r, r') dsds +2 J f [V x Sz(r)],[V' x So(r/)]zGo(r, r') ds ds' (5.9) LZ = 2jkoZo S (Hi x z) ds (5.10) where the subscript z in Bs denotes to take the z component. It is noted that Ls is removed in case of radiation problems and that the same holds for Ke when the scattering problem is considered. 5.3 Edge-Based Prismatic Elements Consider the right angled prism shown in fig. 5.3 whose vertical (z-directed) sides are parallel (right-angled prism). We now design two geometric quantities as ^= ei - x (r-ri S=1 (i.11) where ri denotes the location of the ith node, in is the unit vector along the ithl triangular edge opposite to the ith node, i. denotes the length of this edge and r is any position vector terminated inside the triangle. One way to obtain an edge — based field representation for the prism is to utilize the nodal basis functions [52]

86 and then apply the procedure discussed in [49,53,54]. However, an alternative and more physically meaningful approach can be employed for the construction of the edge elements. Referring to fig. 5.2, it is evident that if r is in the x-y plane, then Se in (5.11) gives the area of another triangle 12'3' such that the lengths of edges joining the nodes 2 - 3 and 2' - 3' are equal. With this definition of r, we define a vector a vector Si = x (r-r,) (5.12) S2S se where ei S6 i = -. That is, the vector component along ei has a magnitude which is equal to the ratio of the areas of the triangle 12'3' to that of 123. We observe that (5.12) is simply the edge-based expansion for the triangular elements [32] and is the appropriate expansion to be used in (5.5). The corresponding volumetric basis functions can be obtained by inspection, viz. - (Az i1,2,3 (z0 + Az - z) VJ Sz j = 4,5,6 (5.13) Vk = (k k =7,8,9 where (k is the triangle simplex coordinate associated with the kth prism vertex at (xk, oyk). As illustrated in fig. 5.3, z, and h = Az represent the offset coordinate arind the prism height, respectively. When (5.13) are substituted into (5.7), the resulting integrals can be evaluated in closed form as given in the Appendix C. However, the integrals resulting from the substitution of (5.12) into (5.9) must be carried out numerically, except the self-cells which can be performed analytically as discussed by Wilton [55].

87 5.4 Applications Thin slot antenna structures have been treated using the above formulated FE-BI technique and certain modeling results will be presented in this section to demonstrate the validity and capability of the approach. Radiation and scattering by an Annular Slot: To evaluate the accuracy and efficiency of the prismatic mesh and the aforementioned implementation, we first consider the analysis of the narrow annular slot (0.75cm wide) shown in fig. 5.4. The slot is backed by a metallic circular cavity 24.7 cm in diameter and 3 cm deep. The FE-BI method is quite attractive for this geometry because the slot is very narrow and most of the computational requirements are shifted on the finite element portion of the system. The calculation shown in fig. 5.5 were carried out using the prismatic and tetrahedral elements [29]. As seen, they overlay each other. However, only 1024 prisms were needed for modeling the cavity, whereas the number of the tetrahedral elements for this homogeneously filled cavity were 2898 for acceptable element distortion. If a multi-layered structure was considered, or a similar system condition was used as a criterion for mesh generation, then much more tetrahedrals than prisms would be needed for modeling such a structure. Mloreover, the prismatic mesh is trivially generated given the slot outline. In contrast, substantial time investment is required for generating and post-processing the tetrahedral mesh. Frequency Selective Surfaces: Frequency selective surfaces (FSS) structures [56,57] are arrays of tightly packed periodic elements which are typically sandwiched between dielectric layers. The periodic elements may be of printed form or slot configurations designed to resonate

88 probe b a=12.35 cm b=0.75 cm ',, / po= 7.7 cm a,' 0.7 < f < 1 GHz d=3 cm i ~e,p1.35 Figure 5.4: Geometry of the annular slot antenna backed by a cavity 23.7 cm "in diameter and 3 cm deep at specific frequencies. As such, they are penetrable around the element resonances and become completely reflecting at other frequencies. To meet bandwidth design specifications, stacked element arrays may be used in conjunction with dielectric layer loading. Here we shall consider the analysis of FSS structures (with slot elements) via the FE-BI method. Because of the fine geometrical detail associated with the FSS surface, the finite element method has yet to be applied for the characterization of FSS structures, but use of prismatic elements makes this a much easier task. Of particular interest in FSS design is the determination of the transmission coefficient as a function of frequency, and since the array is periodic, it suffices to consider a single cell of the FSS. For computing the transmission coefficient T, the periodic cell is placed in a cavity as shown in fig. 5.6 and the structure is excited by a plane wave impinging at normal incidence. Assuming that near resonance the wave transmitted through the FSS screen will retain its TEM character, the transmission line concept can be used to find the scattered field CaT2 ES 1 -caR

89 -10 --20 - ^-3 0 Prism Elements O - Tetrahedral Elements *-3 Mn -40 -.0 m -50 - -60 -70 -80 0 10 20 30 40 50 60 70 80 90 100 theta (degree) 10 -,- I I,,,0 --10 -.m -20 - "-30. -40 co Ir -50 - -60 - -70 - -80 -l --- -100 -80 -60 -40 -20 0 20 40 60 80 100 theta (degree) Figure 5.5: Scattering: Bistatic (co-pol) RCS patterns computed using the tetrahedral FE-BI code and the prismatic FE-BI code. The normally incident plane wave is polarized along the q = 0 plane and the observation cut is perpendicular to that plane. Radiation: X —pol and Co-pol radiation patterns in the ( = 0 plane from the annular slot antenna shown in fig. 5.4. The solid lines are computed using the tetrahedral FE-BI code whereas the dotted lines are computed using the prismatic FE-BI code. The excitation probe is placed at the point (y=O) marked in fig. 5.4.

90 tf radome cover / (on e=4.5 substrate) absorber metal backwall Figure 5.6: Illustration of the setup for computing the FSS transmission coefficient Upper figure: periodic element (top view); Lower figure: periodic element in cavity (cross-sectional view)

91 where T is the transmission coefficient of the FSS, R = 1 - T and a is the reflection coefficient associated with the cavity base. To reduce the multiple interactions within the cavity, it is appropriate to terminate the cavity with some absorber, thus reducing the value of a to less than 0.1. Since R is also small near resonance, a good approximation for T is T() - lO log E dBa and upon considering the next higher order cavity interactions, we have TCB -T + 10log [1 -ca(l -T(0))]. A more direct and traditional computation of TdB would involve the placement of the FSS element in a thick slot [58]. However, this requires enforcement of the boundary integral over the entire lower surface of the slot, leading to a much more computationally intensive implementation. The above FSS modeling approach was applied for a characterization of single layer and multi-layer FSS structures. In both cases, the periodic element was a slot configuration. The geometry of the single layer periodic element is shown in fig. 5.6 and consists of a planar slot array on a dielectric layer 0.0762 cm thick and having 6r = 4.5. The FE-BI calculation using prismatic elements is given in fig. 5.7. Clearly, our calculations are in good agreement with the measurements and data based on the more traditional method of moments [59, 60]. The geometry of the multilayer radome considered in our study is given in fig. 5.8. The total thickness of the FSS was 6.3072 cm and is comprised of two slot arrays (of the same geometry) sandwiched within the dielectric layers. For modeling purpose, a 1.54cm thick absorber is placed below the FSS as shown in fig 5.8. From the

92 3 S II X -20.............1c -25 0 2 4 6 8 10 12 14 16 18 20 Frequency (GHz) Figure 5.7: Calculations and comparisons of transmission through the FSS structure shown in fig. 5.6 calculated results, it is seen that the results generated by the FE-BI method are in good agreement with the measurements. Radiation Property study of Conformal Slot Spiral Antenna: Consider a typical Archemidean slot-spiral antenna shown in fig. 5.9. This antenna is built on a double-sided PCB with its two arms following the expression: r = aO + S3, where a = 0.1333cm and 3 = 2.8595cm. One arm can be determined from the other by rotating 180~ counterclockwisely. It is noted that this structure differs from the conventional design in that the central portion of the spiral is not fabricated. The reasoning for it relies on the facts that the antenna is designed with a bandwidth less than 30%, and that the central portion usually requires a careful fabrication because of the geometric details, and still that the central space may be used for possibly complex feed network. One of our goal is to study the effect of this

93 slot array placed 60mils below top surface of 90mils layer \slot array at 30mils below top of 90mils layer ~r \ / 90mils 90mils: 5.85cm 2cm -I \ absorber, ~ r=. - j 6.45cm ~0 ~,E e — 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figure 5.8: Upper figure: geometry of the multilayer frequency (FSS) used for modeling; lower figure: measured and mission coefficient through the FSS structure selective surface calculated trans

94 spiral shape on its performance. A benchmark test model is designated to operate from 118 MHz to 157 MHz. to replace the conventional protruding blade antenna. The size however is much compact with its conformality property. Our simulation model is scaled by 1/8 to operate at 944MHz to 1256MHz with the center frequency 1100MHz. The values of a and d3 above were determined based on this frequency band and also the number of turns (4.5). The cavity is filled with a dielectric slab (e, = 2.2) of 0.3 cm depth, corresponding to approximately 0.011 free space wavelength at the center frequency. The antenna's directivity is analyzed from the radiated pattern at lower, center and higher frequencies and the results are tabulated in Table 5.1. Figure 5.10 - 5.12 show the radiation patterns for frequency 944, 1100 and 1256 MHz, respectively. The Es and En at the principle plane & = 90~ are plotted. It is understood that when the frequency varies, the active region travels along the slot spiral. Thus the principle plane may not be coincident with the E-plane. (In fact, the E-plane is not clearly defined in this case.) The optimum axial ratio for the three cases are tabulated also in Table 5.1, and it shows that the spiral shape design really plays an important role to insure a good quality radiation pattern. iAt both center and lower frequencies, less than 3 dB AR has been achieved. When the frequency increases, the active region moves inwards to the center and becomes closer to the feeds where the EM fields exhibit a comparatively strong profile. The radiated pattern therefore is most likely affected and this explains why the AR increases at the high frequency. The AR deterioration can be avoided by adding a couple of spiral turns inside. It is seen, nevertheless, that a CP mode can be achieved within the entire designated bandwidth and with a wie azimuthal angle (as wide as 60~ in the optimum case). In practice, we notice that absorbing materials may be needed to

95 regulate the magnetic currents at the beginning or ends of the slot spiral, especially when the number of turns is minimized. Frequency (GHz) 0.944 1.100 1.256 Gain (dB) 7.22 6.66 5.23 Axial Ratio (dB) 2.7 1.0 3 Table 5.1: Comparisons of gain and axial ratio at different operating freqcuencies 5.5 Concluding Remarks A hybrid finite elemnent-boundary integral (FE-BI) formulation was presented for modeling narrow slots in metal backed cavities. Prismatic elements were used in connection with the FE-BI implementation, and in contrast to the tetrahedral elements, these offer several advantages. Among them, low sampling rates are needed for generating meshes and the mesh generation process is substantially simplified. Other advantages of the prismatic elements over the tetrahedral elements include better system conditions and faster pre/post data processing. The explicit expressions for FE-BI implementation of prismatic elements were tabulated and numerical results for slot antennas and frequency selective surfaces were presented to demonstrate the validity and capability of the technique.

96 Figure 5.9: Illustration of a typical 2-arm slot-spiral design. solid line: Ephi, dashed line: E_theta Figure 5.10: Radiation Pattern at f=1. GHz (center frequency design). A good axial ratio is achieved up to 60~ degree.

97 solid line: Ephi, dashed line: Etheta Figure 5.11: Radiation Pattern at f=0.944GHz (lower end of frequency range). It can be seen that the axial ratio of the pattern becomes larger compared to that at the center frequency, but still remains within 3dB for a wide angle range. This indicates that the number of the outer turns in the spiral contour design is most likely sufficient. frequency=1.256 GHz... — - 0 - -.... I ~~. -10 -20 \ \ 30 -30: X *. / solid line: Ephi, dashed line: Etheta Figure 5.12: Radiation Pattern at f=1.256GHz (higher end of frequency range). It can be seen that the axial ratio of the pattern is deteriorated comparec to those at the center frequency and lower frequency. This certainly shows that the number of inner loops still needs to be increased to insure a good quality pattern.

CHAPTER VI Antenna Feed Modeling For scattering problems where the plane wave incidence is usually considered as the 'source', the right-hand-side excitation has been explicitly expressed in (3.7) and (5.10) and will not be discussed further. However, for antenna impedance evaluations, we have proposed and reported several feeding schemes [61] associated with various practical feed designs for microstrip antennas. Some of these are discussed below. 6.1 Probe Feed 6.1.1 Simple Probe Feed For thin substrates the coaxial cable feed may be simplified as a thin current filament of length I carrying an electric current I 1. Since this filament is located inside the cavity, the first term of the integral in (3.7) needs to be considered for this model. Specifically, the ith (global) entry of the excitation vector KI^ becomes IKZi = jkoZoI I V (r), i z= ji i,2 * jm where r is the location of the filament, m is the number of (non-metallic) element edges and jn is the global edge numbering index. In general, m such entries are associated with m element edges, and thus the index i goes from jl up to jm. This 98

99 expression can be further reduced to Ki = jkoZoI 1, provided that the ith edge is coincident with the current filament. 6.1.2 Voltage Gap Feed This excitation is also referred to as a gap generator and amounts to specifying a priori the electric voltage V across the opening of the coax cable or any other gap. Since V = E d, where d is a vector whose magnitude is the gap width, and E the V electric field across the gap, we have that Ei = -cosO, where cosOi is equal to 1 if the dcos 0 ith edge is parallel to d. Numerically, this gap voltage model can be realized by first setting the diagonal term A1 equal to unity and the off-diagonal terms A (i - j) to zero. For the right-hand-side vector, only the entry corresponding to the ith (global) edge across the gap is specified and set equal to the value Ei whereas a]l other entries associated with edges not in the gap are set to zero. 6.2 Aperture-coupled Microstrip Model As shown in fig. 6.1, when the microstrip antenna is fed with a microstrip line network underneath the ground plane (cavity's base) via a coupling aperture, special treatment of the feed structure must be considered in the FEM formulation. This is because the microstrip line is usually designed to have different size and shape as compared to the cavity's geometries. Hence, the conventional simulation of treating the entire 3-D domain using a single type of elements is not efficient or appropriate for this feed. Referring to fig. 6.1, it is appropriate to separate the computational domains because of the small element size required in modeling the guided feed structure. One difficulty encountered when this decomposition is implemented is how to model the coupling through the aperture. As an example let us consider a rectangular

100 aperture which has been extensively employed in practice. The cavity fields may be discretized( using tetrahedral elements, whereas in the microstrip line region rectangular bricks are the best candidate since the feed structure is rectangular in shape and the substrate is of constant thickness. Although both types of elements employ Antenna Elments I I I' St (II) Truncation Plane Coupling Aperture Figure 6.1: Cross-section of an aperture coupled patch antenna, showing the cavity region I and the microstrip line region II for two different FEM computation domains. (a) (b) (c) Figure 6.2: slot and its discretization (a) slot aperture; (b) typical mesh from cavity region; (c) uniform mesh from microstrip line region. edge-based field expansions, the meshes across the common area (coupling aperture) are different, and this causes difficulty in enforcing field continuity across the slot aperture. However, since the aperture is very narrow, a 'static' field distribution may be assumed at any given frequency. Therefore, the potential concept may be again applied to relate the fields below and above the aperture. Specifically, the

101 'equi-potential' continuity condition is enforced, and to proceed to do so, let us first classify the slot edges as follows Tetrahedral Mesh (Cavity Region I): * E1 j 1,2, 3,... vertical edges * E2 j =1,2,3,... diagonal edges Brick Mesh (Feed Region II): * E j - 1,2,3,... vertical edges Then the 'equi-potential' continuity condition requires that Ebl = CjE (6.1) t E2 t (jE + ej+ E ) (6.2) - 2d 3 3) in which +1 ~j = -, — 1 whereas t and d are the lengths of the vertical and diagonal edges, respectivelv. That is, t is simply the width of the narrow rectangular aperture. The coefficient 6j is equal to ~1 depending on the sign conventions associated with the meshes above and below the coupling aperture. The connectivity scheme for entirely different computational domains may be extended by generalizing this concept. It is apparent that this approach makes the FEM implementation straightforward for different geometry/size domains that would be significantly inefficient if only one type of elements were used for modeling the structure. In addition, the technique ensures a good system condition since the number of distorted elements in the mesh are minimized.

102 6.3 Coax Cable Feed 6.3.1 Motivation The coax cable is widely used as a feeding structure for microstrip or cavitybacked patch antennas because of its simplicity and low spurious radiation. Indeed, abundant literature exists on the theoretical and experimental investigation of coax cable feeds [62-64]. Most of these papers present integral equation techniques in conjunction with the pertinent Green's functions. However, the Green's function is only available for a certain class of geometries, and this limits the application of the integral techniques to those geometry designs. Also, the formulation must be modified and recoded for different antenna configurations corresponding to Green's function variations. To avoid the complexity of the Green's function, we recently proposed a hybrid finite element - boundary integral approach [29] which is described in chapter 3 and 4. For antenna radiation, it is observed that a simple probe model with a constant current along the inner conductor linking the grounded base to the antenna element is straightforward and efficient. But the probe feed is only valid for thin substrates and this is consistent with the Moment Method (MM) results. To model an electrically thick substrate, in this section a more sophisticated feed modeling scheme is proposed in the context of the finite element method (FEM) using linear edge-based tetrahedral elements. The formulation of the entire hybrid numerical system will be first described in the presence of the necessary functional term for feedline. The proposed feed model is then presented on the basis of a TEMVI mode excitation. Model improvements are also discussed for the case when the TEM assumption at the cavity-cable junction does not hold.

103 6.3.2 Hybrid FE-BI System The functional pertinent to the radiation by a cavity-backed antenna with a coax cable feed (as shown in fig. 6.3) is given by F(E) = - Jli {(VxE) -(V x E) -k orE EE dv 2J -2kfJ (E x ). { J (E x $) ( + VV) Go(rr')dSl} dS 0 k2Js \ o / -jkoZo (E xH) dS, (6.3) where V refers to the cavity volume and the surface S encompasses the cavity aperture excluding the portion occupied by the metallic antenna elements; or and tr denote, respectively., the relative permittivity and permeability; ko is the free space wave number, Zo is the free space intrinsic wave impedance, I is the unit dyad, and Go(r, r') is the free space Green's function with r and r' denoting the observation cavity-cable junction. Following the standard discretization procedure [29], we obtain a system of equations of the form Ne se NNe sE J { [ AB] E }} + 0 [E} (6. ) e=l eES eEC where the explicit expressions for Aij and B.j may be found in [29] and the functional term Fc is the surface integral on C in (6.3). 6.3.3 Proposed Coax Feed Model To proceed with the evaluation of Fc =-jkoZo t (E x H) dS, (6.)

104 a boundary constraint relating E to H is needed. To this end, we assume a TEMI mode on C and consequently the fields within the cavity may be expressed as (see fig. 6.4) Io Zo 1 Io 1, E = (1( + F) r, H -(1 (6.6) 2wV c6 r 2w r where c,, is the relative permittivity inside the coax cable; F denotes the reflection coefficient measured at z = 0 and Io is the given input current source at the same location. Also, (r, I, z) are the polar coordinates of a point in the cable with the center at r == 0. To simplify the analysis, we introduce the quantities loZo Io eo = Io ( + F), ho = (1 - F. (6. 7) 2/- — rc 27 6 Hence, E =-r, H =-, (6.8) and from (6.7) it follows ho =- eo + - (6.9) Zo 7r which is the desired constraint at the cable junction in terms of the new quantities ho and Ce. Note that eo and ho are field coefficients as new unknowns in place of the fields E and H, and it is therefore appropriate to rewrite Fc in terms of these new coefficients. To do so, we substitute (6.7) and (6.9) into (6.5) and upon making use of the axisymmetric field property we obtain Fc - -27jkoZoeohorln(-), (6.10) a where a and b are the radii of the inner and outer cable conductors. The superscript src stands to indicate that ho is treated as a source term in the extremization of the functional.

105 We choose the linear edge-based tetrahedral elements to discretize the cavity and the corresponding mesh on the cross section C is shown in fig. 6.4(b). In this formulation, the field across the pth edge, p=1,2,..., Nc (Nc is the number of cavity mesh edges on C), is set to a constant as dictated by the linear edge-based expansion function inside the cavity. However, the cable TEM modal fields (6.6) behave as 1/r and this modeling inconsistency makes it difficult to apply the tangential field continuity condition at the cable junction ( i.e. over the aperture C). To overcome these difficulties, we can relate the fields across the cable junction by recognizing that the potential difference between the inner and outer conductors must be the same as computed by the fields of the cavity or those in the cable region. Specifically, if the Nph edge of the cavity region borders and is also across the coax cable, from (6.6)-(6.8) it then follows that the appropriate equi-potential condition is b AV = E(b - a) = eoln-, i = Np(p = 1, 2,..., Nc). (6.1t1) a where AV denotes the potential difference between the inner and outer surface of the cable. This condition simply provides a relation between the constant cavity edge field and the coax cable modal field. When used into the functional Fc, it introduces the excitation into the hybrid finite element system without a need to extend the mesh inside the cable or to employ a fictitious current probe. It remains to rewrite FC in terms of E, i.e. the field value of the edges bordering the cable and to d(o so we substitute (6.9) and (6.11) into (6.10). Then upon taking into account all Nc cavity mesh edges on the cable junction, we obtain 7- Io ' -rcb - a Fc -- jkoZo(b -a) -- b E E (6.12) a i=Np In this expression, rather than representing the functional Fc in terms of a single edge field, we made use of the average field across the cable as computed by the totality

106 of the equal element fields on the cable's aperture (because of the axisymmetric property, all elements fields at the cable's aperture are equal). The factor inside the curly brackets of (6.12), with the superscript src, functions as a source in the extremization process. Hence, the extremization of (6.12) yields &Fc 1 k Z j o ccb-a A 3_ = -1rj koZo(b - a) { - - / b - EJ } -EJ 3 - 1 - Zo ln j i = hE,- V, i Np(p = 1,2,...Nc), (6.13) where Uj = j7ko (b )2 (6.14) a Vi = j koZo(b - a)Io. (6. 15) We observe that the 'constant cavity field' along each mesh edge at the cable junction is just a fictitious field representation and its meaningful physical interpretation is governed by the equ-i-potential constraint (6.11). To proceed, we assemble the FEM system together with (6.13). Specifically, each Uj is added to the Nc diagonal entries of the finite element matrix which is associated with the NC edges bordering the coax cable. Also, the excitation column of the hybrid system is nullified everywhere except for the NC entries which are set to Vi. Once the hybrid FE-BI system is solved [29], the input admittance at z=0 is calculated from Yin = H Mr. rd 2Io 1 e oln ~ ~- *(6.16) eoln b_ Zc where Z, is the characteristic impedance of the coax cable. In the above feed model we assumed the presence of only the dominant(TEM) mode at the cavity-cable junction, an assumption which may not be suitable for

107 certain applications. To overcome this limitation, one approach is to extend the mesh (say, a distance d) into the cable. The equi-potential condition will then be applied at z=-d, where all higher order modes vanish. This scheme requires a more suitable expansion for the fields in the -d < z < 0 section to avoid the complication of extending the tetrahedral mesh into the cable and, thus, retain the efficiency of the equi-potential feed model. Since in most cases the antenna is operated in a frequency range far below the cut-off of the first higher order mode of the coax cable, the field distribution near the junction C will still be dominated by the fundamental TEM mode. With this understanding, a possible suitable expansion for the field in the coax cable (using shell elements rather than tetrahedrals) is N-(r, (), z) E E E ) (6.17) q i where q=r,, or z, i=1,2,3 or 4 and N'(r, 4, z) is the shape function for each of the 12 edges (3 directionsx 4 edges per direction). They are given by Nq = (qb - q)(q - q)a (6.18) with q,,qb and qc representing r, i and z in cyclic rotation and correspondingly qa, qb and qc represent the parameters r, ( and 5. Also, i denotes the edge number along each coordinate, and Aqa is the width of the edge along the qa direction. The correspondence between the edge numbers and the node pairs for each coordinate(r, ( or z) is given in Table 6.1 along with the definition of the tilded parameters in (6.18). When an axisymmetric field property is assumed, the numerator of the expansion in (6.17) reduces to the standard brick element format for the radial and z components, independent of the ( variable. Note also that the particular property of this expansion is the introduction of the 1/r factor, simulating the coaxial cable mode. The accuracy of (6.17) is demonstrated in fig. 6.5, where we show that only 2-3 elements

108 TA.BLE1 c(oordila.t es rnode parTamleters 1 (1-2) + )1+ _ zI + A z r 2 (4-3) ~ i z1+ Az 3 (8-7) + _i z_ 4 (5-6) -< + A z0 z 1 (1-4) 4 rl+Ar z1 +A: 2 (5-8) r- + 'Ar z 3 (6-7) -- rt1 4 (2-3) - 1l z1 + Az z 1 (1-5) +- -r + A - r ( + A 0 2 (2-6) T- rI -i- + A — 3 (3- 7) {- r / 01 L (4-8) - ri +A r 0i ___ Table 6.1: The correspondence between the edge numbers and the node pairs for each coordinate(r, 0 or z) along with the definition of the tilded parameters in (6.18). are needed along the radial direction for the accurate prediction of the dominant field distribution. When compared to the conventional linear tetrahedral elements, the efficiency of this expansion is apparent (i.e., many more tetrahedrals are needed to model the same cable region). 6.3.4 Results and Conclusion To validate our proposed feed simulation, two circular patch antenna configurations were used for calculation. One patch antenna was of radius 1.3 cm printed atop of a circular cavity (radius=2.1 cm) filled with a dielectric (6,=2.9) material 0.41 cm deep. For this patch, the feed was placed 0.8 cm from the center and the input impedance was measured over the band 2 - 5 GHz. In fig. 6.6 we compare the measured input impedance with data computed on the basis of the proposed equipotential feed model. Clearly, the results from measurements and the equi-potential model are in excellent agreement whereas the probe model yields substantially inaccurate results near resonance. Figure 6.7 shows the comparison between measurements and calculations for an

109 other patch antenna whose input impedance was measured by Aberle and Pozar [65]. This patch had a radius of 2.0 cm and the 0.218 cm thick substrate had c6=2.33 ard a loss tangent tan6=0.0012. The feed was located 0.7 cm from the center, and for our FE-BI calculation the patch was placed in a circular cavity of 2.44 cm in radius. As shown in fig. 6.7, the equi-potential model is again in excellent agreement with measurements, as opposed to the results by the probe model in [65]. In conclusion, the presented equi-potential feed model has been shown to be extremely accurate in modeling coax feed structures. Moreover, its implementation in the context of a finite element formulation is very simple and as easy to implement as the probe feed. It was also demonstrated how the proposed feed model can be generalized to the case of asymmetric feed structures where evanescent modes may be present. 6.4 Conclusion In developing numerical techniques for antennas, the feed network is one of challenging problems to solve in consideration of accuracy, efficiency and simplicity. This is primarily because the antenna feed in fabrication has certain instrumental uncertainties on one side. On the other, as aforementioned, the accuracy of numerical results is usually extremely sensitive to the feed model, the feed location, sampling rate around the feed point, and so on. The proposed numerical feeds in this chapter resemble the practical systems as closely as possible, and 'with a thorough consideration of their numerical implementations, we realize that they can be used for mostly encountered antenna problems. As an addition to the group of feed models, we also developed a circuit modal feed which coincides with domain truncations. Since this model has to do with microwave

110 z A Patch Ground plane Aperture No Y Cavity Cavity - ~~ Coax cable opening (surface C) X x Coax cable Figure 6.3: Illustration of a cavity-backed patch antenna with a coax cable feed. cavity A patch I, hh r4 2b 2 2a T cavity-cable junction (a) (b) Figure 6.4: (a) Side view of a cavity-backed antenna with a coax cable feed; (b) Illustration of the FEM mesh at the cavity-cable junction (the field is set to zero at the center conductor surface).

111 Comparison of E-Field Along (1-Wavelength) Coax Cable co 0: -0 ken v] a) U-. LJ LL segment number (a) Comparison of E-Field Distribution Along R PMn( 5500 t, 5000 iI 4500 N ir 4000 c -o < 350CI a LL - 2500 o~ 2000 1500 I IK 1001'' 0.02 0.04 0.06 distance in cm from the center (b) 0.08 Figure 6.5: Field distribution in a shorted coax cable as computed by the finite element method using the expansion (6.18). -: analytical; xxx: numerical. (a) Field coefficient eo along the length of the cable (leftmost point is the location of the short); (b) Field along the radial coordinate calculated at a distance A/4 from the shorted termination.

112 0 0 u 0a q) E -E ct C. C; 180. 140. 100. 60.0 20.0 -20.0 '1.00 2.00 3.00 4.00 5.00 frequency (GHz) (a) 180. U u r0 0l C. a 0~ C S >1 ct to Cd L, 140. 100. 60.0 20.0 -20.0 1-.00 1.00 2.00 3.00 4.00 5.00 frequency (GHz) (b) Figure 6.6: Measured and calculated input impedance for a cavity-backed circular patch antenna having the following specifications: patch radius r=13mm; cavity radius R=21.1mm; substrate thickness t=4.1mm; ec=2.4; and feed location Xf=0.8 cm distance from center. Results based on the simple probe model are also shown for comparison. Our modeling retains the vertical wire connection to the patch and uses the incoming coaxial mode field for excitation. (a) Real part; (b) Imaginary part.

113 120. 0 100 E0 - 80 0 0 x | 60- x\ C 40 o 1x 0 -~ o o 620 0 - -~20.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 frequency (GHz) (a) 100. 80 - E 60 - o x ~ 40 ~ E 20 -, 00 modeL [65] (a) Real part; (b) Imaginary part. -20 I-40 O O 0 -60 -8,,,, S5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 frequency (GHz) (b) Figure6(3.7: Measured and calculated input impedance for a circular patch antenna having the following specifications: patch radius r-2cm; substrate thickness d=0.21844cm; feed location from center xf-0.7cm; cr-2.33; tan6-O.0012. [65]. ~ measurement; xxx: this method; o o o: probe model [65] (a) Real part; (b) Imaginary part.

114 circuit, we shall leave the topic to the next chapter in conjunction with other related subjects.

CHAPTER VII Circuit Modeling Many designs of microstrip antennas require certain feedlines to carry electromagnetic signals from the source. Finite element methods are also suited and applicable to these wave propagation problems. This chapter is devoted to circuit modeling. After an overview of recent mesh termination techniques, we discuss a numerical de-embedding process appropriate for the finite element analysis, and then turn to the topic of mesh truncations for circuit simulation. 7.1 Introduction One of the most important aspects of finite element implementations is the truncation of the computational volume. An ideal truncation scheme must ensure that outgoing waves are not reflected backwards at the mesh termination surface, i.e. the mesh truncation scheme must simulate a surface which actually does not exist. To date, a variety of non-reflecting or absorbing boundary conditions (ABCs) have been employed for truncating the computational volume at some distance from the radiating or scattering surface, and applications to microwave circuits and devices have also been reported. The ABCs are typically second or higher order boundary conditions and are applied at the mesh termination surface to truncate the computational 115

116 volume as required by any PDE solution. Among them, a class of ABCs is based on the one-way wave equation method [66,67] and another is derived starting with the Wilcox Expansion [68,31]. Also, higher order ABCs using Higdon's [69, 70]formulation and problem specific numerical ABCs have been successfully used, particularly for truncating meshes in guided structures [71]. There are several difficulties with traditional ABCs. Among them is accuracy control, conformality, ease of parallelization and implementation difficulties when dealing with higher order ABC's. Also, the applications of ABC's in microwave circuit modeling requires a priori knowledge of the propagation constants which are typically not available for high density packages. An alternative to traditional ABCs is to employ an artificial absorber for mesh truncation. Basically, instead of an AB(C, a thin layer of absorbing material is used to truncate the mesh, and the performance for a variety of such absorbers have been considered [72, 73]. Nevertheless, these lossy artificial absorbers (homogeneous or not) still exhibit a non-zero reflection at incidence angles away from normal. Recently, though, Berenger [74] introduced a new approach for modeling an artificial absorber that is reflectionless at its interface for all incidence angles. In two dimensions, his approach requires the splitting of the field components involving normal (to the boundary) derivatives and assigning to each component a different conductivity. In this manner an additional degree of freedom is introduced that is chosen to simulate a reflectionless medium at all incidence angles. Provided the medium is lossy, this property is maintained for a finite thickness layer. Berenger refers to the latter as a perfectly matched layer(PMIL) and generalization of his idea to three dimensions have already been considered [75, 76]. Also, implementations of the absorber for truncating finite difference-time domain(FDTD) solutions has so far been found highly successful. Nevertheless, it should be noted that Berenger's PML does

117 not satisfy Maxwell equations and cannot be easily implemented in finite element (FEM) solution. A new anisotropic(uniaxial) artificial absorber [77] was introduced recently for truncating FEM meshes. This artificial absorber is also reflectionless at all incidence angles. Basically, by making appropriate choices for the constitutive parameter tensors, the medium impedance can be made independent of frequency, polarization, and wave incidence angle. A PML layer can then be constructed by introducing sufficiernt loss in the material properties. The implementation of this artificial absorber for truncating finite element meshes is straightforward and, moreover, the absorber is Maxwellian. 7.2 Numerical De-embedding De-embedding presented here is a numerical process used to extract certain circuit quantities. Specifically, we are interested in S-parameters for a uniform transmission line terminated with any loads denoting the possible discontinuities, which may arise from line-to-line or line-to-antenna couplings. The dominant transmission line mode is assumed at and near a reference plane Sref in this discussion. Consider a transmission line of certain length as shown in fig. 7.1. With an appropriate shielding scheme, the line is included in the computational domain. The full wave analysis provides the E field distribution anywhere including the region along the line. One is therefore able to represent E field along the transmission line with respect to the locations to get V(z) Vie- + I ez (7.1) where V is proportional to the magnitude of E with Vi being the incoming and Vr the reflected wave amplitude, z is measured from the reference plane Sref, 7 is the

-118 Microstrip Line Shielded FEM Region l Substrate Probes Figure 7.1: Illustration of a shielded microstrip line. propagation constant to be determined and R = Vr/Vi is the reflection coefficient, or S11. Since in (7.1) 7y, Vi and Vr are three independent quantities that characterize the wave propagation and reflection, to determine them one needs to specify the field values for V(z) at three points z_, zo and z+ of equal inter-distances V(z-) V(zo) V(z+) z+ - Z = z - z_ (7.2) along the line. To simplify the problem, we choose the reference plane right at the center point such that zo = 0 and z+ -z- = d. Given the three field values front FEM computations, it follows V(d) = ie-d + t1ed (7.3) V(0) = V + V (7.4) V(-d) = Ve' + I e —d (7.5) To solve for y, we first add (7.3) to (7.5) to get (V1 + Vr) (ed + e-d) = V(d) + V(-d) (7.6)

119 Then eliminating ~V + Vr from (7.4) and (7.6), we obtain cosh(ad) (d) + V(-d) (7.7) 2V(0) from which 7 can be determined. The effective guided wavelength Ag and effective dielectric constant eeff may then be calculated by 27r A0 Ceff = with = Im {(} and Ao being the free space wavelength. From(7.3) and (7.4), Vi and Vr are expressed as i = V(0)ed - V(d) (7.8) 2 sinh(7d) Vr = V(0)- 1 (7.9) Therefore, the reflection coefficient becomes S11 V (7.10) This de-embedding process is suited for one port network analysis. However, the technique may be readily extended to two-port networks. For instance, on the assumption of a perfect termination or match at port 2, once Vi is determined, S21 can be obtained by V/Vi, where VO is the outgoing wave at port 2 predicted by FEM. As mentioned before, S-parameter evaluations depend on termination methods. Low quality terminations result in prediction errors and make the analysis less reliable. Therefore, high performance termination methods are always desirable and we next discuss this issue.

120 7.3 Truncation Using DMT As already indicated, S-parameters from a possible discontinuity region along a transmission line may be extracted at a distant reference plane, where there exists only a dominant mode. For shielded microstrip lines at the input port (#1) and output port (#2), (similar to that in fig. 7.1), the modes underneath the lines are given by Eo(x)(e-1z + Reoz) z C Sn Ey(x) - (7.11) TEo(x)e2z z C Sout where Eo0(x) denotes the field distribution of incident wave at the incident plane (port 1) and R is the reflection coefficient at the same plane, T represents the transmission coefficient measured at the plane Sout (port 2), and y7, 72 are the effective propagation constants at port 1 and port 2, respectively. For truncating the FEM mesh at a specific port, it is necessary to first determine the E field pattern across the shielded structure. This can be accomplished by assuming a static model shown in fig. 7.2, where the static potential satisfies V2 = 0 0 = Vo on metallic line (7.12) where E = -V<. Sove this standard PDE model, and with a tedious mathematical derivation, it is finally found that En=, An c x) sinh (-y z > d (, y) - (7.13) sinh(-d) ( _n, =,odd An r cos (-xsih -( -y)) < d sinh (b - d)) a a where sin (W) 4, -V - 2a jA - FVo2 n~fn F

121 Y! ~ d x 0 ~-a/2 a/2 Figure 7.2: Illustration of the cross section of a shielded microstrip line. with F S= ( 2) sinh -d) n2 n2fn a n fn = -r2 a )Co0sh ( (b-d) snh(b-dfvb d)) a + cr1COSh (ayd A complete FEM system may now be constructed by introducing the EM fields at Sin to truncate the computational domain. This truncation simultaneously introduces an excitation to the numerical system and the S-parameters may then be extracted by measuring the field distributions at the input and output ports as mentioned before. 7.4 Truncation Using PML Below, we begin with a brief presentation of the artificial absorber, and this is followed by an examination of the absorber's performance in terminating guided structures and volume meshes in scattering problems. Results are presented which show the absorber's performance as a function of thickness/frequency and for differ ent loss factors.

122 7.4.1 Theory Consider the waveguide, shielded microstrip line and scatterer shown in fig. 7.3. Of interest is to model the wave propagation in these structures using the finite Electric Probe / = -j Z = x — -Jl I./ I' d (a). waveguide (a). waveguide =xy= y = a -jip d+t=40cm t= 5 cm cross-section: 4.755x2.215 cm Microstrline Absorbing layer - -- --- ed=3.2 (b) Microstrip line i --....- t(b). Microstrip Line d + t = 12.0 cm L = 2.38 cm h = 0.21 cm H = 1.06 cm w = 0.548 cm Figure 7.3: A rectangular waveguide (a) and a microstrip line (b) truncated using the perfectly matched uniaxial absorbing layer. element method. For a general uniaxial medium, the functional to be minimized is jF = V x E. (-1 * V X E) - k2~ E E dV Jv J +ot Ex (r V x E) - dS, (7.14') in +Sout in which fl, and i denote the permeability and permittivity tensors whereas E is the total electric field in the medium. The surface integrals over Sin and Sout must be evaluated by introducing an independent boundary condition and the ABC serves for this purpose but alternatively an absorbing layer may be used. An approach to evaluate the performance of an absorbing layer for terminating the FE mesh is to

123 extract the reflection coefficient computed in the presence of the absorbing layer used to terminate the computational domain. In this study we consider the performance of a thin uniaxial layer for terminating the FE mesh in a rectangular waveguide and a microstrip line. Such a uniaxial layer was proposed by Sacks et.al. [77] who considered the plane wave reflection from an anisotropic interface (see fig. 7.4 ). If Region I Region 2 Reflected wave " Transmitted wave ----------- -------------- a, O O a 0 0 Incident wave b b2 O 0 C, 0 0 C2 v z Figure 7.4: Plane wave incidence on an interface between two diagonally anisotropic half-spaces. ptr and Cr are the relative constitutive parameter tensors of the form a2 0 0 /r= r b2 0 (7.15) 0 0 C2 the TE and TM reflection coefficients at the interface (assuming free space as the background material) become cos0, - ~ cosOt RTE V cosO, + -VcosO 2 2 (7J16) cosOt - cosi ( RTM a cosOi + J2 oOt

124 1 and by choosing a2 b2 and c2 -= it follows that RTE = RTM = 0 for all incidence angles, implying a perfectly matched material interface. If we set a2 = a- j3, the reflected field for a metal-backed uniaxial layer is IR(0)l = -23ktcosOt (7.17) where t is the thickness of the layer and Oi is the plane wave incidence angle. The parameter ak is simply the wavenumber in the absorber. Basically, the proposed metal backed uniaxial layer has a reflectivity of -30 dB if /3tcosO^ = 0.275A or - 55dB if 3tcosO- = 0.5A, where A is the wavelength of the background material. The reflection coefficient (7.17) can be reduced further by backing the layer with an ABC, rather than a PEC. However, the PEC backing is more attractive because it eliminates altogether the integrals over the surfaces. Clearly, although the interface is reflectionless, the finite thickness layer is not and this is also true for Barenger's PML absorber. Below we present a number of results which show the performance of the proposed uniaxia] absorbing layer as a function of the parameter /3, the layer thickness t and frequency for the guided structures shown in Figure 7.3. We remark that for the microstrip line it is necessary to let a2 = eb( - j/3) for the permittivity tensor and a2 = rb(i - ji,3) for the permeability tensor, where erb and Prb are the relative constitutive parameters of the background material (i.e. the substrate). 7.4.2 Results Rectangular Waveguide Let us first consider the rectangular waveguide shown in fig. 7.3. The guide's cross — section has dimensions 4.755 cm x 2.215 cman and is chosen to propagate only the TE10 mode. It is excited by an electric probe at the left, and fig. 7.6 shows the mode field

125 strength inside the waveguide which has been terminated by a perfectly matched uniaxial layer. As expected, the field decay inside the absorber is exponential and for /3 values less than unity the wave does not have sufficient decay to suppress reflections from the metal backing of this 5cm layer. Consequently, a VSWR of about 1.1 is observed for / = 0.5. However, as 13 is increased to unity, the VSWR is nearly 1.0 and the wave decay is precisely given by e-/ktcos0 = e /~cosSP where t is the wave travel distance measured from the absorber interface, P = 2/3t/Ag and here 0i = 44.5~. It is noted that when / is increased to larger values, the rapid decay is seen to cause unacceptable VSWR's. One is therefore prompted to look for an optimum decay factor for a given absorber thickness and fig. 7.7 provides a plot of the TE1o mode reflection coefficient as a function of 2,3t/A,, where we chose to normalize with respect to the guided wavelength Ag. Figure 7.7 is typical of the absorber performance and demonstrates its broadband nature and the existence of an optimum value of /3 for minimizing the reflection coefficient. Basically, the results suggest that 3 must be chosen for a given absorber thickness to provide the slowest decay without causing reflections from the absorber backing. That is, the lowest reflection may be achieved when the entire absorber width is used to reduce the wave amplitude before it reaches the absorber's backing. As expected, this optimum value of 0 changes with frequency but the broadband properties of the absorber are still maintained since acceptable low reflections can still be achieved for unoptimized / values. For example, in the case of f 4.5GHz (dashed line) the optimum value of /3 = 1 gives a reflection coefficient of -45dB whereas the value of d = 3 (corresponding to 23t/Ag = 2.3) gives a reflection coefficient of -37 dB which is still acceptable for many applications. It should be noted though that setting 3 = 3 allows use of an absorber which is about 2cm or 1/3 free space wavelengths. Also, as can be realized

126 the discretization rate plays a role in finding the optimum value of /3t/Ag and thus the presented curves refer to a sampling rate of around 18/Ag for the waveguide example. Not surprisingly (see (6.6)), for this example, the value of a does not play an important role in the performance of the absorber and this is demonstrated in fig. 7.8. As seen, setting always a = 3 gives the same performance as the case of a = 1 shown in fig. 7.7. Our tests also show that other choices of a give the similar absorber performance. However, it is expected that a will play a role in the presence of attenuating modes and it is therefore recommended to choose a = a to ensure that all modes are absorbed. Microstrip line 'The performance of the perfectly matched uniaxial layer in absorbing the shielded microstrip line mode is illustrated in fig. 7.9 where the reflection coefficient is plotted as a function of 23tt/Ag, where AX = AX/ /eff and Ceff is the effective dielectric constant. In this case, the microstrip line is terminated with a 1.87 cm thick, 5 -layered absorber and the line is extended up to 4 layers inside the absorber to avoid an electric contact with the metallic wall. Similarly to the waveguide, we again observe that an optimum /3 value exists and it was verified that in the absorber the wave exhibits the same attenuation behavior as shown in fig. 7.6. The reflection coefficient at the optimum 3 = 1 is now -42dB and if better performance is required, a thicker absorbing layer may be required. Again as in the case of the waveguide example, the value of a plays little role in the performance of the absorber and this is illustrated in fig. 7.10. However, of importance is the behavior of the reflection coefficient as a function of 2/t/Ag. For the waveguide and microstrip examples, we observe that the absorption is maximized for approximately the same value of 2/3t/Ag

127 (about 0.8). Thus these curves can be used for other applications as well, although it should be noted that the discretization rate plays an equally important role and this needs further investigation. The accuracy and validity of the PML applications for circuit parameter computations can also be seen from the result illustrated in fig. 7.11. It is seen that use of the optimized 4.5 cm PML layer, with a = 1 and d = 1, yields very accurate input impedance values. The shown microstrip line impedances were computed by measuring the vertical field at the probe's location without a need to extract the VSWR which is often difficult with unstructured finite element meshes. Note that the shielded microstrip line dimensions for the data are given in fig. 7.11. Meanderline Another example is the meander line shown in fig. 7.12. For the FEM simulation, the structure was placed in a rectangular cavity of size 5.8mm x 18.0mm x 3.175mrnm. The cavity was tessellated using 29 x 150 x 5 edges and only 150 edges were used along the y-axis. The domain was terminated with a 10 layer PML, each layer being of thickness t = 0.12mm. The S11 results are shown in fig 7.13 and are in good agreement with the measured data [78].

128 300 ----- I (D.0 0c C O3 0 (A.U -.0 'I.L U. 250 200 150 100 50 0 -50 - - Real Part - - Imaginary Part Magnitude '\\ / \ /\ / \ /,,/ \ _ / \x I \,-< I/ / Electric Probe -1 -Ex= y= ( 1- ijp -1,c=4y= z= 1-jp d+t=40cm t=lcm 1WIIf3 0 10 20 30 40 50 60 70 Samples Along Waveguide (14samples/4) Figure 7.5: Typical field values of TElo mode inside a rectangular waveguide terminated by a perfectly matched uniaxial layer. Figure 7.6: 60 62 64 66 68 70 72 74 76 78 80 Segments number along waveguide (alpha=1.0, f=4.5GHz) Field values of the TElo mode inside a waveguide terminated by a perfectly matched uniaxial layer. The absorber is 10 elements thick and each element was 0.5 cm which translates to about 13 samples per wavelength at 4.5 GHz.

:129 C ~ -5 --- f=4.0GHz --- f=4.5GHz - f=5.0GHz -10F.- -15 -20 m-25 *o - -30 cr -35 F -40 F -45 -50 ', 8 1 1 Figure 7.7: 0 2 4 6 8 10 12 14 2pt in Xg (o=1.0) Reflection coefficient vs 2o3t/A9 = 1) for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6.... I I I I I -5F -10 'a -15.-5 E3 > -20 m -25 -: -30 -35 -40 -45 -50 - I I I I I I 0 2 4 6 8 10 12 14 2pt in Xg (cc=P) Figure 7.8: Reflection coefficient vs 2,3t/Ag, with a = /, for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6.

-130 0 0.5 1 1.5 2 2.5 3 3.5 4 2pt in kg (o=1.0) Figure 7.9: Reflection coefficient vs 2Pt/A9 with a=1, for the shielded microstrip line terminated by the perfectly matched uniaxial layer. -5 --- f=4.0GHz --- f=4.5GHz -10 f=5.0GHz._ -15 0 I -20 E m-25 - -30 - o: 0 0.5 1 1.5 2 2.5 3 3.5 4 2pt in kg (oa=P) Figure 7.10: Reflection coefficient vs 2,3t/Ag with c = 3, for the shielded microstrip line terminated by the perfectly matched uniaxial layer.

1131 u r I I I I I I 65 60 55 - Theory -— FEM Absorbing layer Microstnp line h dt. 12.0 cm, n......~i i }-p..... ii H L. 238cm 11h1 0.21 cm I -H 1.06 cm 3.2 w 0.47 cm I~ 3.2 50 E -45 N 40 - 35 - 30 - 25,), I i I I Ii 0.2 0.4 0.6 0.8 1 width of microstrip (cm) 1.2 1.4 1.6 Figure 7.11: Input impedance calculations for the PML terminated microstrip as compared to the theoretical reference data.

132.305.305 161 _1_ 1.525 _1_1.61J......... I I I I Y (unit: mm) LO co Figure 7.12: Illustration surement. 1 0.9 - 0.8 - of a meander line geometry used for comparison with mea I I I I 0.7h / 7 + 0.6 -. -0.5 0.4 0.4 + 0.3 0.2 0.1 I I 5 10 15 20 25 30 frequency (GHz) Figure 7.13: Comparison of calculated and measured results shown in fig.7.12. for the meander line

CHAPTER VIII AWE: Asymptotic Waveform Evaluation 8.1 Brief Overview of AWE Although full wave electromagnetic systems are large and cumbersome to solve, typically only a few parameters are needed by the designer or analyst. A reduced order modeling of these parameters (input impedance, S parameters, far field pattern, etc) is therefore an important consideration in minimizing the CPU requirements needed for generating the frequency response of the parameter. The Asymptotic Waveform Evaluation(AWE) method is one approach to construct a reduced ordering model of the input impedance or other useful electromagnetic parameters. AWE relies on a Pade approximation of the given parameters to avoid the repeated solution of the system at each frequency value. It has already been applied to problems in circuit analysis and in this paper we demonstrate its application and validity when used in conjunction with the finite element method to simulate full wave electromagnetic problems. The method of Asymptotic Waveform Evaluation (AWE) is a reduced-order modeling of a linear system and has already been successfully used in VLSI and circuit analysis to approximate the transfer function associated with a given set of ports/variables in circuit networks [79-82]. The basic idea of the method is to de 1_33

134 velop an approximate transfer function of a given linear system from a limited set of spectral solutions. Typically, a Pade expansion of the transfer function is postulated whose coefficients are then determined by matching the Pade representation to the available spectral solutions of the complete system. In this chapter we investigate the applicability of the AWE method for approximating the response of a given parameter in full wave simulation of radiation or scattering problems in electromagnetics. Of particular interest is to use AWE for evaluating the input impedance of the antenna over an entire bandwidth from a knowledge of the full wave solution at a few (even a single) frequency points. Also, the method can be used to fill-in a backscattering pattern with respect to frequency from a priori knowledge of the simulation system with a few data samples of that pattern. Given that practical partial differential equation(PDE) systems involve several thousand unknowns, AWE can indeed have a dramatic reduction of CPU requirements in generating a response for a given system parameter (state variable) without a, need to resolve the system for the fields in the entire computational grid. Below we first describe the recasting of the FEM system for application of the AWE. We then proceed to describe the AWE method and demonstrate its application, accuracy and efficiency for computing the input impedance of a shielded microstrip stub. 8.2 Theory 8.2.1 FEM System Recast The application of the finite element method to full wave electromagnetic solutions amounts to generating a linear system of equations by extremizing the functional [83] F =< V x E,. V x E> -k2 < E, E > +kb.t. (8.1)

135 where <,> denotes an inner-product and b.t. is a possible boundary term whose specific form is not required for this discussion. Also, the dyadics a and A are material related coefficients, k = - is the wavenumber and w is the corresponding, c operating frequency with c being the speed of light. A discretized form of (8.1) incorporating the appropriate boundary conditions is [29] 2(8.2) (Ao + kA1 + k2A2) {X} = {} (8.2) where Ai denote the usual square (sparse) matrices and {f} is a column matrix describing the specific excitation. Clearly (8.2) can be solved using direct or iterative methods for a given value of the wavenumber. Even though Ai is sparse, the solution of the system (8.2) is computationally intensive and must be further repeated for each k to obtain a frequency response. Also, certain analyses and designs may require both temporal and frequency responses placing additional computational burdens and a repeated solution of (8.2) is not an efficient approach in generating these responses. An application of AWE to achieve an approximation to these responses is an attractive alternative and below we formulate AWE in connection with the FEM system (8.2), whose implementation is considered in connection with antenna and microwave circuit problem. For these problems it turns out that the excitation column {f} is a linear function of the wavenumber and can therefore be stated as {f} - {fi} (8.3) with {fi} being independent of frequency. This observation will be specifically used in our subsequent presentation.

136 8.2.2 Asymptotic Waveform Evaluation To describe the basic idea of AWE in conjunction with the FEM, we begin by first expanding the solution {X} in a Taylor series about ko as {X} = {Xo} + (k - ko) {X } + (k - ko)2 {X2} + +(k- ko) {X} + 0 {(k - ko)} (8.4) where {,Xo} is the solution of (8.2) corresponding to the wavenumber ko. By introducing this expansion into (8.2) and equating equal powers of k in conjunction with (8.3), after some manipulations we find that {Xo} = koA1 f1} {Xi} = Ao [{fl} - A {Xo}- 2koA2 {Xo}] {X2} = -Ao [A1 {X1} + A2({Xo} + 2ko {X1})] (8.5) {X = -Ao [A1 {X1_} + A2({X-2} + 2ko {X/_})] with Ao = Ao + koA1 + koA2 (8.6) Expressions (8.5) are referred to as the system moments whereas (8.6) is the system at the prescribed wavenumber (ko). Although an explicit inversion of Ao'1 may be needed as indicated in (8.5), this inversion is used repeatedly and can thus be stored out-of-core for the implementation of AWE. Also, given that for input impedance computations we are typically interestedinterested in the field value at one location of the computational domain, only a single entry of {Xi(k)} needs be considered, say (the pth entry) XI(k). The above moments can then be reduced to scalar form and

137 the expansions (8.5) become a scalar representation of X'f(k) about the corresponding solution at k0. To yield a more convergent expression, we can instead revert the moments to a Pade expansion, which is a conventional rational function in form. A special case of the qth order of such an expansion is given by Xyp() - ao + al(k - k) + a2(k - ko)2 +.+ + aq(k - (87) 1 + bi(k - ko) + b2(k- ko)2 + -...- + bq(k - ko)q where ai and bi ( = 0 1,..., q) are referred to as the Pade coefficients. For transient analysis, it is observed that the Pade expansion can be reformulated by partial fraction decomposition [82,84] as q Xp(k) = AX^ ~ E r (8.8) X ) q + k - (8.8) i=1 where Xgqo is the limiting value when k tends to infinity. Clearly, this is the representation suitable for time/frequency domain transformation. The residues and poles (r, and ko + ki) in (8.7) or (8.8) correspond to those of the original physical system and play important roles in determining the accuracy of the approximation. In general, higher order expansion contains more system residues and poles and usually provides a better approximation. Since the accuracy of AWE relies on the dominant residues and poles located in a complex plane closest to the point on the real axis ko from the origin, in practice the number of poles (and residues) needed to obtain a sufficiently accurate expansion can be much smaller than that of the original numerical system, which is the beauty of AWE method. (Detailed analysis and theory of Pade expansion can be found for instance in [85].) For hybrid finite element - boundary integral system, the implementation of AWE is more involved because the fully populated submatrix of the overall system may be associated with a more complex dependence on frequency. In this case it is attractive to instead generate the full submatrix by introducing a spectral expansion

138 of the exponential boundary integral kernel to facilitate the extraction of the system moments. This approach does increase the complications for implementing AWE. It however remains far more efficient in terms of CPU requirements when compared to the conventional approach to continuously repeating the solution of the entire system. 8.3 Numerical Implementation As an application of AWE to a full wave electromagnetic simulation, we consider the evaluation of the input impedance for microstrip stub shielded in a metallic rectangular cavity as shown in fig. 8.1. As expected, the stub's input impedance is a strong function of frequency from 1-3 GHz and this example is therefore a good demonstration of AWE's capability. The shielded cavity is 2.38cm x 6.00cm x 1.06cm in size and the microstrip stub resides on a 0.35cmr thick substrate having a dielectric constant of 3.2. The stub is 0.79cm wide and A/2 long at about 1.8 GHz. We note that the cavity is terminated at the perfectly electric conductor (PEC) back wall by an artificial absorber having relative constants of er = (3.2, -3.2) and,J = (1.0, -1.0). In this study the artificial absorber was used for setting up an appropriate forced problem rather than to establish a perfectly matched interface. Nevertheless the numerical FEM system was already demonstrated valid and accurate for microwave circuit analysis [86]. The frequency response of the shielded stub was first computed using a full wave finite element code from 1 to 3 GHz at 40MHz intervals (50 points) to serve as the reference solution. We then chose the single input impedance solution at 1.78GHz in conjunction with the 4th order and 8th order AWE representation given in (8.8) to approximate the reference response. As seen in fig. 8.3, the 4th order AWE

139 representation is in agreement with the real and reactive parts of the reference input impedance solution over about 56% and 33% bandwidth, respectively. This clearly shows that the contributions of the system poles in the complex k plane lead to an accuracy difference to the real and reactive components. Surprisingly, the 8th order AWE representation recovers the reference solution over the entire 1-3GHz band for both impedance components. We also observed that the CPU requirements for 4th and 8th order computations are nearly the same except for a few more times of matrix-vector products. The number of these product operations is in the order of the AWE approximation order q and therefore much smaller than the size of the original numerical system. It is also apparent that to demonstrate the AWE efficiency we only solved the system once at one frequency point. The save of CPU time can be easily estimated when compared to solve the system conventionally for each frequency over the entire band. Thus, the AWE representation is an extremely useful addition to electromagnetic simulation codes and packages when a wideband frequency response of the system is required. The development and utility of the method for more complex numerical systems and multiple parameter simulation can be readily extended and will be considered in the future.

140 1.06 cm 6cm I 2.38cm A Figure 8.1: Illustration of the shielded microstrip stub excited with a current probe. 50 - 45 -40 -/ 035 / -30/ \ E / 25 / 20' / 15 10,,-,,,,1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 frequency in GHz Figure 8.2: Impedance calculations using traditional FEM frequency analysis for a shielded microstrip stub shown in figure 8.1. Solid line is the real part and the dashed line denotes the imaginary part of the solutions. These computations are used as reference for comparisons.

141 E!30 0- 2 E / - Exact " 20 / - - AWE-8th order - /- - -AWE-4th order" 15 // /,'/ 10 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 frequency in GHz (a) 50, I I i I I I,,, -- 45 - Exact - - AWE-8th Order 40 - - - AWE-4th Order 20 // I 30, i E 25 - 20~ 15 -10 —,-,-,-,-, -,-, — 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 frequency in GHz (b) Figure 8.3: 4th order and 8th order AWE implementations using one point expansion at 1.78 GHz are shown to compare with the reference data. With the 4th order AWE solutions, 56% and 33% bandwidth agreement can be achieved for the real (a) and imaginary (b) parts of impedance computations, respectively. It is also shown that the 8th order solutions agree excellently with the reference data over the entire band. (a) Real Part (b) Imaginary Part computations

CHAPTER IX Conclusions 9.1 Discussion on the Research Work During the period of developing the hybrid finite element methods, many expected and unexpected issues were frequently encountered. Among them are the understanding of physical systems, development of mathematical models, interpretation of results, lack of measurement data for comparison, and increased computational demands, etc. We can comfortably state that significant progress has been made during the course of this work. Some of our accomplishments are summarized below. * GENERAL PURPOSE HYBRID FE-BI METHOD DEVELOPMENT Once the FE and BI subsystems and the hybrid method were mathematically formulated, a major effort was then devoted to the integration of the two subsystems. The interface between the FE-BI program and a commercial (SDRC-IDEAS) mesh generator was developed with minimum but sufficient geometry and meshing data. The latter task was important in permitting the geometrical modeling and meshing of printed antenna configurations of arbitrary shape. It is this general version of the FE-BI code that can (in theory) be used to simulate any planar conformal antenna. 142

143 * ITERATIVE SYSTEM SOLVER A memory saving algorithm ITPACK was intertwined with the hybrid FEM subsystem to register only the non-zero FEM entries. The BiCG iterative solver was independently developed for partially sparse and partially full mnatrices in conjunction with the ITPACK algorithm. * UNIFORM GRID BI SUBSYSTEM - BiCG-FFT To facilitate the efficient storage and evaluation of the BI sub-system, a uniform right triangular zoning scheme for discretization of the boundary integral equation was introduced by re-numbering the triangle edges as dictated by their geometrical locations. This approach leads to a BI sub-system which could be cast as a 2-D discrete convolution, thus allowing use of FFT for fast execution of the iterative solver. This truncation/termination the "exact" evaluation of rectangular and right-triangular patches. * NON-UNIFORM BI SUBSYSTEM - OVERLAY-BICG-FFT For non-rectangular patches, an interpolation scheme was proposed to make use of the efficient BiCG-FFT technique by overlaying a fictitious uniform grid with the original arbitrary mesh. The forward/backward transformation matrices to account for field interpolations using localized basis functions were derived and they were indeed highly sparse. * FEED MODELING Feed modeling is one of the most important and challenging tasks in the context of the general purpose FEM. To this end, a series of commonly used feed structures were modeled using the hybrid technique, especially in consideration of efficiency and accuracy. These include probes/generators, aperture coupled

144 slotline, microstrip line, coax cable. etc. * PRISMATIC FEM ELEMENTS INCORPORATION A major problem in any hybrid FEM analysis is the tedious pre-processing for mesh generation. Thin layer substrates in the presence of thick spacer(s) are often found in practical conformal antenna designs. However, this typical configuration leads to large numerical systems when tetrahedral elements are used. To alleviate these difficulties, the prismatic edge-based elements were developed and incorporated in the hybrid system. This formulation exhibits certain features/advantages that tetrahedral FEM does not. It can therefore be used to compensate the tetrahedral FEM as a subsystem module. * MESH TRUNCATIONS WITH DMT AND PML The uniaxial or otherr anisotropic medium simulation may be readily accomplished using the proposed hybrid FEM technique due to the geometrical adaptability of the tetrahedral elements. Hence the PML was first introduced into the 3-D FEM. Various performance studies were carried out to optimize the application of the PML to microwave circuit simulations. In the meantime, an analytical approach, dominant mode truncation (DMT), was proposed and implemented as an alternative mesh truncation of the FEM domain for microstrip lines and similar structures. * REDUCED ORDER APPROXIMATION-AWE AWE has been reported useful in RLC and VLSI applications. For wideband and highly varying frequency responses, this technique is particularly efficient. Given the promise of the method for broadband simulations of VLSI circuits, we consider its application to electromagnetic system. In particular, AWE

145 was incorporated into the finite element method. It was indeed observed that the attractive features of AWE are maintained when used in electromagnetic problems. 9.2 Suggestions for Future Tasks The following is a list of suggested tasks for further development of the finite element methods * HIGHER ORDER EDGE-BASED FEM DEVELOPMENT * ADAPTIVE ELEMENTS * MIXED ELEMENTS AND INTERFACE * ANISOTROPY (WITH LOSS) FEM INVESTIGATIONS/APPLICATIONS * INCORPORATION OF MORE ROBUST TRUNCATIONS * MODULAR DEVELOPMENT AND INTEGRATION (WITH USER INTERFACE) 9.3 Modular Development Hybrid finite element methods for the analysis of various electromagnetic problems encountered in practice are still on the way to reach its maturity. As well known, any general purpose technique (such as the commercially available software in electromagnetics) either loses its efficiency or becomes incapable when simulating intricate problems. It is anticipated that at the current stage of the FEM development with the limited capacity of computing resources, more and more specialized techniques will be desired, particularly when efficiency and speed become a key consideration in large scale computations and in engineering design.

146 With a whole set of specially developed techniques and methodologies, one should then consider to create an integration environment. As shown in fig. 9.1, we propose this FEM modular environment for future computational electromagnetic applications. A well designed modular finite element methods will be the most capable and robust in the future! FEM Multi-Module Environment User Interface User's Modules ---------- --------- Mesh FEM Modules FeedLine ' vbricks Truncation prisms Modules Modules tetras I -------------—... —... — Post Processor. Figure 9.1: Multi-modular FEM environment

APPENDICES 147

148 APPENDIX A Evaluation of Matrix Elements for Tetrahedrals Referring to Fig. A.1 and the associated table, the fields in the eth tetrahedron 1 2 3 nodes/vertices (a) Table Edge Vertex Numbering Numbering il ir2 (D 1 2 ( 1 3 a 1 4 (~ 2 3 4 2 _ 3 4 (b) Figure A.1: (a) A tetrahedron. (b) its local node/edge numbering scheme

149 are expanded as 6 E = y EEWT 2-1 where the basis functions W7 are given by f7-i + g7-i x r r C Ve w~_,(r) = 0 outside element b7-, f,-i = 6 rV r, r2 position vectors of vertices i1 and i2 (see Table) 6le b bi g7-i = 6ee i I (r*2 — r ) t2 bi = Iri -- r = length of the ith edge (see Table) = element's volume We note that V. W = 0 V x W7 = 2gi indicating that W7 are divergenceless. Furthermore, f 1 i=j W7 (r) ej - j; i j 6 -0 i j where r3 has its tip on the jth edge of the tetrahedron. This last property ensures that the coefficients Ei = E. e^ represent the average field value at the ith edge of the tetrahedron. Using the above basis functions, we now proceed with the derivation of the matrix elements Ai. We have IJ -(v x W )- (V x W - )= -gi~ glVe JJJv r J P~r 3

150 Also, Er W> WJ d = - 1 {(fi fj) + (r.D) + (g x r) *(g. x r)} dv = r(II + I2 + 3) where D -= (f x gj) + (fj x gi) and = fJj f.f dv 2 = J r Ddv 13 = JJ (g x r) (gj x r)dv Since f is a constant vector, I1 reduces to I1 fi' fj To evaluate I2 we first set 4 4 4 x = Lixi; y = ~ Liyi; z = izi i=1 i=l i=1 where xi, y, z (i = 1,...,4) denote the (x, y, z) coordinates of the tetrahedron's vertices and Lz are the simplex coordinates or shape functions for the same element. That is, L, is the normalized volume of the tetrahedron formed by its three corners other than the ith, and the point (x, y, z) located within the tetrahedron. Using the standard formula for volume integration within a tetrahedral element and simplifying, we have I2 Dx x + Dy yi + i ili=1=

151 where Dm is the mth component of D. The evaluation of I3 can be simplified by the use of basic vector identities. We have 13 = gi gj/ Ir2 dv - / (gi r)( r)dv = (giygjy + gzgjz) j x2 dv + (gixgr + gjzgz) 2 dv + (gigjx + giyg3y) j 2 dv -(gzxgjy +gjxgiy) xydv - (glgjz + gjgzz) j zx dv - (gsg + gjygiz) J yzdv where gim represents the mth component of the vector g.. Each of the above integrals can be easily evaluated analytically and the result can be expressed in the general form I alar dv = I aliami + ali am 20 a z z 20I =l i=l i=l for 1, = 1,...,3. The parameters al or an can represent any of the rectilinear variables x, y, z.

152 APPENDIX B Evaluation of the Boundary Integral System Matrix The explicit representation of the boundary integral subsystem matrix is given by (3.12) and can be rewritten as Bpqi = 2 I f (-kSi. Sj + V x Si. V x Sj) Go(r, r') dS' dS (B.1) where Go(r, r') is the free space Green's function, and TiP is the pth triangle of the triangle pair Sp as shown in fig. 3.3. Similar to the finite element assembly procedure, it should be recognized that the definition of (B.1) virtually involves an assembling over the triangles. To proceed, (3.14) is used to discritize the field region and thus its curl is given by V x S,(r) =6(r) A (B.2) p where c(r) is defined by (3.15). Note that when deriving (B.2), the fact that r is located inside the pth triangle in a planar surface is considered and therefore V r = 2. Given the Green's function, it is straightforward to express the matrix entries as -, j (r - r)-(r- r ) e (r) j ( r)0 d S' 8 7rA AJp J R + 2AA f-/- W -(r)Ej2(r) dS'dS (B.3) 27-A 7JT4 R

153 in which R = r -- r'l. These integrals can be readily evaluated for non-self-cell terms by numerical integrations. It is also observed that once Tf coincides with T, the integrands become singular because of the Green's function. In this case, the singularity should be removed. For the second integral, this is accomplished by subtracting and adding an additional term. That is rr rr eko0~R f n R-:sw_\ iRJ J3iJ - i (r)c (r) dS' JdS + J r - e(r)Ji(r) dS' dS (B.4) R3 The first integral in (B.4) is evaluated using numerical integrations and the second one is carried out analytically [55]. Similarly, the first integral in (B.3) is rearranged as rr r r -jko R fJr r ) (r -r i(r)-e _ (r) RdS' dS ' - 1R f(r - r) (r - rj(r)(r) j(r) dS' dS 3 R + f f(r-r (r -rc(r)c(r dS'dS (B.5) in which the first integral on the right hand side is numerically integratable with singularity removed and the second one again may be expressed in a simple analytical form [55].

154 APPENDIX C Formulation for Right Angle Prisms For FEM implementation, the following quantities are required PImn V= X Vm VxV dV (c.1) QO Vm2 Vn dV (C.2) emn where the curls are given by V x V =-S- [(x- x)x + (y- y.) -z(z-z )] i -,2,3 Vx:Vj, =2 Az [(x-x)x + (y-yj)Y + z(z + Az-z)] j=4,5,6 (C.3) V X Vk = [(Xk2 - Xkl) + (Yk2 - Ykl)Y] k = 7, 8, 9 To this end, we follow the notation defined in (4.13) and (4.14), where ii'1,2,3 represent the top triangle edges, jj'=4,5,6 denote the bottom triangle edges and k, k'7,8,9 stand for the vertical three edges. It is found that (C.2) and (C.3) can be analytically evaluated and we tabulate the results as follows P = C D,Az + Se(Az)3] (C.4) pi 4 33 [ 4 3

155 Az1 (C.6) 4kk 34S P. P=-C-. [D.Az - ~Se(Az) (C.-7) (C.8) P~~k Pki - A.[sik(SX -~Se) + l'SY-yS) (C) Itkp 1 k(SX _ X.Se) + A2 tk(SyY.-ySe)] (C.I() (C.I1) - (A )30D(C. 12) (Z3 -C( -,D-(C.-13) (C. 14) Q kk' AZSeTkk, (C.1I5) Q (,Az)O. D. (C.-16) Qik Q k Q k- Qkj - 0 (C.17) where -kk 1/6 'for k - k; 1/12 for k: V C_ (C.18) 13 4(SeAZ)2 Dj ~SXX - (x~ + x -)SX + x ix.S' + SYY -(y. + y3-)SY + y-y -Se The remaining quantities in the above list of the (expressions are defined as -jdxdy SX jxdxdy

156 SY = y dxdy SXX = x2 dxdy Se SYY = j y2dxdy SXY Se SXY = j xydxdy These integrals can be expressed in terms of the global coordinates of the three nodes ('Xi, Yi), (Xj, Yj), (Xm, Ym). Specifically, assuming that the three nodes ij and m of a triangle are in counterclockwise rotation, we then have, 1 xi yi S dxdy 1 x yj 1 Xm Ym SX = f xdxdy = (Xi + iXj + XMn) SY = y dxdy = -(Y + s( + X;) Use 1 SXX = 2dxdy = 12 {(Xi + +X2+(X2 + Xj+Xm) } SYY dy {(Yx + xj + x)2 + (y2 + Y2 + 2)} SXY = xydxdy =-{(X ( + X + Xm) (Y + Y + m) + (XYi + XjYj + Xy~m~)}

157 APPENDIX D System Derivation From A Functional Referring to fig. 2.2, we begin with the functional (H) = 1 / (V x H V x H-k H H) dfQ (D.I) Q2 d+Qf+ +Qa to derive the system in terms of the scattered field. On inserting the field decomposition H = Hsat + Hinc, the functional becomes (HSCt Hinc) (Hscat Hscat) +(H scf H") d (H zc H S) + (Hsct, Hinc)f +(HiC, Hsct)lQf (D.2) where the fact that H iC does not exist in 2Qa has been considered, and (,) represents the integral of the same form as in (D.1). Once a self-adjoint system operator is assumed, it then follows that (HinC Hscat) /= (Hscat, Hinc) d (D.3) Also, in free space Qf, (Hi~, Hat) IQ-= (H"t n Hi )| (D.4) Upon invoking the divergence theory, we have 2(Hsc, H H)t = HSia x cf V x H"c) dS x H ) dS (D / H cat (n2 x ef ~ V x Hc) dS (D.5) f

Evidently, for a self-adjoint operator, one readily recovers the system (2.42) obtained via Galerkin's method. It should be noted that besides the boundary/transition conditions, the self —adjoint property of a system operator simply requires - (-1) r =(r) (D.6) In the case of a non-self-adjoint operator, it is generally not possible to recover the system given by (2.42) in the same manner. This is because the functional (D.2) in terms of the scattered field is of the form f(H) = 1 ( V x Hsc - tc2HS at. H.sca d 2 JQd+Qf+Qa + (V x Hscat * V x Hinc- k2Hscat - Hinc) d + 1 (v x H * ( *V Hscat - 2OHic I ' Hscat) dQ + I t' H (i x I V x HifC) dS - Ha. (n2 x e V x Hic) dS (D.7).f:1 It is observed that the first integral shows the same form of the FEM system as that in (2.42). All other integrals in (D.7) contribute to the system excitation. Apparently, the two integrals over domain Qd are not identical, leading to a different FEM system than (2.42).

BIBLIOGRAPHY 159

160 BIBLIOGRAPHY [1] E.H. Newman and P. Tulyathan. "Analysis of microstrip antennas using moment methods". IEEE Trans. Antennas and Propagat., Jan. 1981. [2] J.R. Mosig and F.E. Gardiol. "Analytical and numerical techniques in green's function treatment of microstrip antennas and scatterers". IEE Proc., Pt. H, 1983. [3] M.C. Bailey and M.D. Deshpande. "Analysis of elliptical and circular microstrip antennas using moment methods". IEEE Trans. Antennas Propagat., Sept. 1985. [4] C. Tsai, H. Massoudi, C.H. Durney, and M.F. Iskander. "A procedure for calculating fields inside arbitrarily shaped inhomogeneous dielectric bodies using linear basis functions with the moment method". IEEE Trans. Microwave Theory Tech., vol. 34, pp. 1131-1139, Nov 1986. [5] T.S. Rappaport. "Wireless personal communications: Trends and challenges". IEEE-AP Magazine, Vol. 33, No. 5, October 1991. [6] Naohisa Goto. "The impact of mobile radio communications". IEEE-AP Magazine, April 1992. [7] Christopher R. Johnson. "Computational inverse problems in medicine". IEEE Computational Science and Engineering, Winter, 1995. [8] D. Colton and P. Monk. "A new approach to detecting leukemia: Using computational electromagnetics". IEEE Computational Science and Engineering, Winter, 1995. [9] J.M. 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. Antennas Propagat., vol. 39, pp. 1598-1604, Nov. 1991. [10] R.F. Harrington. Time-harmonic Electromagnetic fields. McGraw-Hill Book Co., 1961. [11] C.T. Tai. Dyadic Green Functions in Electromagnetic Theory. IEEE Press, New York, 1994.

161 [12] P. Silvester. "Finite element solution of homogeneous waveguide problems". Alta Frequenza, vol. 38, pp. 313-317, 1969. [13] P.P. Silvester and G. Pelosi. Finite Elements for Wave Electromagnetics - Methods and Techniques. IEEE Press, 1994. [14] A. Konrad. "Vector variational formulation of electromagnetic fields in anisotropic media". IEEE Trans. Microwave Theory Tech., Vol. 24, No. 9, pp. 553-559, 1976. [15] C.H. Chen and C.D. Lien. "The variational principle for nonself-adjoint electromagnetic problems". IEEE Trans. Microwave Theory Tech., Vol. 28, pp. 878 -886, 1980. [16] G. Jeng and A. Wexler. "Self-adjoint variational formulation of problems having nonself-adjoint operators". IEEE Trans. Microwave Theory Tech., Vol. 26, No. 1, pp. 91-94, 1978. [17] S.R. Cvetkovic and J.B. Davis. "Self-adjoint vector variational formulation for lossy anisotropic dielectric waveguide". IEEE Trans. Microwave Theory Tech., Vol. 34, No. 1., pp. 129-134, 1986. [18] K. Ise and M. Koshiba. "Numerical analysis of h-plane waveguide junctions by combination of finite and boundary elements". IEEE Trans. Microwave Theory Tech., Vol. 36, No. 9, pp. 1343-1351, 1988. [19] K. Hayata, K. Miura, and M. Koshiba. "Full vector finite element formulism for lossy anisotropic waveguide". IEEE Trans. Microwave Theory Tech., Vol. 37, No. 5., pp. 875-883, 1989. [20] Jan A.M. Svedin. "Propagation analysis of chirowaveguides using the finiteelement method". IEEE Trans. Microwave Theory Tech., Vol. 38, No. 10, pp. 1488-1496, 1990. [21] A.A.P. Gibson and J.Helszajn. "Finite element solution of longitudinally magnetized elliptical gyromagnetic waveguides". IEEE Trans. Microwave Theory Tech., Vol. 37, No. 6,pp. 999-1005, 1989. [22] Luis Valor and Juan Zapata. "Efficient finite element analysis of waveguides with lossy inhomogeneous anisotropic materials characterized by arbitrary permittivity and permeability tensors". IEEE Trans. Microwave Theory Tech., Vol. 43, No. 10, October, 1995. [23] Luis Valor and Juan Zapata. "An efficient finite element formulation to analyze waveguides with lossy inhomogeneous bi-anisotropic materials". IEEE Trans. Microwave Theory Tech., Vol. 44, No. 2, February, 1996.

162 [24] A. Chatterjee. "Investigation of finite element - abc methods for electromagnetic field simulation". Ph.D. Dissertation, Radiation Laboratory, University of Michigan, 1994. [25] Jianming Jin. "The Finite element method in electromagnetics". John Wiley & Sons, Inc, 1993. [26] D. Zheng and K.A. Michalski. "Analysis of coaxial fed antennas of arbitrary shape with thick substrates". J. Electromagn. Waves Appl., 1991. [27] X. Yuan. "Three dimensional electromagnetic scattering from inhomogeneous objects by the hybrid moment and finite element method". IEEE Trans. Antennas Propagat., vol. 38, pp. 1053-1058, 1990. [28] 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. 38, pp. 386-393, Mar 1990. [29] J. Gong, J.L. Volakis, A.C. Woo, and H. G. Wang. "A hybrid finite elementboundary integral method for the analysis of cavity-backed antennas of arbitrary shape". IEEE Trans Antenna and Propagat., Sept. 1994. [30] M.L. Barton and Z.J. Cendes. "New vector finite elements for three-dimensional magnetic field computation". J. Appl. Phys., vol. 61, no. 8, pp. 3919-21, Apr 1987. [31] A. Chatterjee, J.M. Jin, and J.L. Volakis. "Computation of cavity resonances using edge-based finite elements". IEEE Trans. Microwave Theory Techn., Nov. 1992. [32] S.M. Rao, D.R. Wilton, and A.W. Glisson. "Electromagnetic scattering by surfaces of arbitrary shape". IEEE Trans. Antennas Propagat., May 1982. [33] J.L. Volakis, J. Gong, and A. Alexanian. "A finite element boundary integral method for antenna RCS analysis". Electromagnetics, vol. 14, no. 1, pp. 63-85, 1994. [34] T.K. Sarkar. "On the application of the biconjugate gradient method". J. Electromagn. Waves Appl., 1987. [35] D.W. Smith and P.E. Mayes. "Numerical and experimental analysis of circularly polarized radiating line antennas". University of Illinois, 1990. [36] H. Morishita, K. Hiraswa, and K. Fujimoto. "Analysis of a cavity-backed annular slot antenna with one point shorted". IEEE Trans. Antennas Propagat., Oct. 1991. [37] F. Croq and D. M. Pozar. "Millimeter-wave design of wide-band aperturecoupled stacked microstrip antennas". IEEE Trans. Antennas Propagat., Vol. 39, pp. 1770-1776, Dec. 1991.

163 [38] T.P. Budka. "microwave circuit electric field imaging systems". Ph.D. Dissertation, Radiation Laboratory, University of Michigan, 1995. [39] J.M. Jin and J.L. Volakis. "Biconjugate gradient fft solution for scattering by planar plates". Electromagnetics, vol. 12, pp. 105-119, 1992. [40] A.F. Peterson, S.L.Ray, C.H. Chan, and R. Mittra. "Numerical Implementations of The Conjugate Gradient Method and the CG-FFT for Electromagnetic Scattering". Elsevier Science Publishing Co. Inc., 1991. [41] Yousef Saad. "Iterative Methods for Sparse Linear Systems". PWB Publishing Company, 1996. [42] R. Barrett, M. Berry, T. F. Chan, J. Demmel, V. Eijkhout, R. Pozo, "Templates for the Solution of Linear Systems: Blocks for Iterative Methods". SIAM, 1994. [43] C. Lanczos. "An iteration method for the solution of the eigenvalue of linear differential and integral operators". J. Res. Nat. Bur. Stand., Vol. 45, pp. 255 -282, 1950. [44] C. Lanczos. "Solution of systems of linear equations by minimized iterations". J. Res. Nat. Bur. Stand., Vol. 49, pp. 33-53, 1952. [45] C.F. Smith, A.F. Peterson, and R. Mittra. "The biconjugate gradient method for electromagnetic scattering". IEEE Trans. Antennas Propagat., Vol. 38, No. 6, pp. 938-1493, 1990. [46] W.H. Press, S.A. Teukolsky, and W.T. Vetterling an B.P. Flannery. "Numerical Recipes the art of sicentific computing". Cambridge University Press, 2nd., 1992. [47] V. Rokhlin, S. Wandzura, and R. Coifman. "The fast multipole method for the wave equation: a pedestrian prescription". IEEE Antennas and Propagat. Magn., Vol. 35, pp.7-12, June 1993. [48] J.M. Jin, J.L. Volakis, and J.D. Collins. "A finite element-boundary integral method for scattering and radiation by two- and three-dimensional structures". IEEE Antennas Propagat. Magazine, vol. 33, no. 3, pp. 22-32, Jun 1991. [49] J.C. Nedelec. 'Mixed finite elements in R/3. Numer. Math., vol. 35, pp. 315-41, 1980. [50] J.P. Webb. "Edge elements and what they can do for you". IEEE Trans. Magnetics, vol. 29, pp. 1460-1465, 1993. [51] A. Bossavit and I.Mayergoyz. "Edge-elements for scattering problems". IEEE Trans. Magnetics, 1989. [52] O.C. Zienkiewicz. The finite Element Method. McGraw Hill, New York, 3rd edition, 1979.

164 [53] A. Bossavit and J.C. Verite. "A mixed FEM-BIEM method to solve 3D eddy current problems'. IEEE Trans. Magnetics, vol. 18, pp. 431-5, Mar 1982. [54] Z.S. Sacks and J.F. Lee. "A finite element time domain method using prism elements for microwave cavities". IEEE Trans. EMC, Vol. 37, No. 4. pp. 519 -527, Nov. 1995. [55] D.R. Wilton, S.M. Rao, A.W. Glisson, and D.H. Schaubert. "Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains". IEEE Trans. Antennas Propagat.. Vol. 32, pp. 276-281, March 1984. [56] E.L. Pelton and B.A. Munk. "Scattering from periodic arrays of crossed dipoles". IEEE Trans. Antennas Propagat., 1979. [57] R. Mittra, C.H. Chen, and T. Cwik. "Techniques for analyzing frequency selective surfaces -a review". Proc. IEEE, 1988. [58] J.M. Jin and J.L. Volakis. "Electromagnetic scattering by and transmission through a three dimensional slot in a thick conducting plane". IEEE Trans. Antennas Propagat., vol. 39, pp. 534-550, 1991. [59] J. Berrie. Personal communication. Mission Research Corp., 1995. [60] L.W. Henderson. "The scattering of planar array of arbitrary shaped slot and/or wire elements in stratified dielectric medium". Ph.D. dissertation, ElectroScience Lab, Ohio State University, 1983. [61] Jian Gong and John L. Volakis. "An efficient and accurate model of the coax cable feeding structure for fem simulations". IEEE Trans Antenna and Propagat., Dec. 1995. [62] J.P. Damiano and A. Papiernik. "Survey of analytical and numerical models for probe-fed microstrip antennas". IEE Proc. H., Feb., 1994. [63] J. Zheng and D.C. Chang. "End-Correction network of a coaxial probe for microstrip paths antennas". IEEE Trans. Antennas Propagat., January 1991. [64] W.C. Chew, Z. Nie, Q.H. Liu, and Y.T. Lo. "Analysis of a probe-fed microstrip disk antennas". IEE Proc. H, 1991. [65] J.T. Aberle, D.M. Pozar, and C.R. Birtcher. "Evaluation of input impedance and radar cross section of probe-fed rnicrostrip patch elements using an accurate feed model". IEEE Trans. Antennas Propagat., Dec. 1991. [66] B. Engquist and A. Majda. "Absorbing boundary conditions for the numerical simulation of waves". Math. Comp., Vol. 31,pp. 629-651, 1977. [67] L.Halpern and L.N. Trefethen. "Wide-angle one-way wave equations". J. Acoust. Soc. Amer., Vol. 84, pp. 1397-1404, 1988.

165 [68] J.P. Webb and V.N. Kanellopoulos. "Absorbing boundary conditions for finite element solution of the vector wave equation. Microwave and Opt. Techn. Letters, Vol. 2, No. 10, pp. 370-372, Oct 1989. [69] R.L. Higdon. "Absorbing boundary conditions for acoustic and elastic waves in stratified media". J. Comp. Phys., Vol. 101, pp. 386-418, 1992. [70] J. Fang. "'ABCs applied to model wave propagation in microwave integratedcircuits". IEEE Trans. MTT, Vol. 42, No. 8, pp. 1506-1513, Aug. 1994. [71] J.-S. Wang and R. Mittra. "Finite element analysis of mmic structures and electronic packages using absorbing boundary conditions". IEEE Trans. Microwave Theory and Techn., Vol. 42, pp. 441-449, March, 1994. [72] T. Ozdemir and J.L. Volakis. "A comparative study of an ABC and an artificial absorber for truncating finite element meshes". submitted to Radio Science. [73] C. Rappaport and L. Bahrmasel. "An absorbing boundary condition based on anechoic absorber for em scattering computation". J. Electromagn. Waves Appl., Vol. 6, No. 12, pp. 1621-1634, Dec., 1992. [74] J.P. Berenger. "A perfectly matched layer for the absorption of electromagnetic waves". J. Comp. Physics, vol. 114, pp. 185-200, 1994. [75] D.S. Katz, E.T. Thiele, and A. Taflove. "Validation and extension to three dimensions of the berenger pml absorbing boundary condition for fd-td meshes". IEEE Microwave and Guided Wave Letters, pp. 268-270, August, 1994. [76] W.C. Chew and W.H. Weedon. "A 3-d perfectly matched medium from modified maxwell's equations with stretched coordinates". Microwave and Optical Technology Letters, pp.599-604, September, 1994. [77] Z.S. Sacks, D.M. Kingsland R. Lee, and J.F. Lee. "A perfectly matched anisotropic absorber for use as an absorbing boundary condition". IEEE Trans. Antennas Propagat., vol. 43, pp. 1460-1463, 1995. [78] I. Wolf. "finite difference time-domain simulation of electromagnetic fields and microwave circuits". International Journal of Microwave and Millimeter- Wave Computer-Aided Engineering, Aug. 1992. [79] S. Kumashiro, R. Rohrer, and A. Strojwas. "Asymptotic waveform evaluation for transient analysis of 3-d interconnect structures". IEEE Trans. ComputerAided Design of Integrated Circuits and Systems, vol. 12, No. 7, pp 988-996, Jul, 1993. [80] R. Kao and M. Horowitz. "Asymptotic waveform evaluation for circuits with redundant dc equations". Technical Report, Computer Systems Laboratory, Stanford University, 1996.

166 [81] L. Pillage and R. Rohrer. "AWE: Asymptotic Waveform Estimation". Technical Report, SRIC-CMU Research Center for Computer-Aided Design, CarnegieMellon University, 1988. [82] E. Chiprout and M. Nakhla. Asymptotic Waveform Evaluation and Moment Matching for Interconnect Analysis. Norwell, Kluwer Acad. Pubs., 1994. [83] J. L. Volakis, A Chatterjee, and J. Gong. "A class of hybrid finite element methods for electromagnetics: a review". J. Electromagn. Waves Applications, 1994. [84] Joseph Lehner. "Partial fraction decompositions and expansions of zero". Trans. Amer. Math. Soc., 1958. [85] E. Chiprout and M. Nakhla. Pade Approximants. Addison-Wesley Pub. Co., 1981. [86] J. Gong, S. Legualt, Y. Botros, J. L. Volakis, and P. Petre. "Application and design guidelines of the pml absorber for finite element simulations of microwave packages". IEEE/MTT Symposium Digenst, 1996.