030601-10-T ROBUST HYBRID FINITE ELEMENT METHODS FOR ANTENNAS AND MICROWAVE CIRCUITS J. Gong and J.L. Volakis National Aeronautics and Space Administration Langley Research Center Hampton, VA 23681-0001 October 1996 30601-10-T = RL-2429

Report #030601-10-T NASA Langley Grant NAG 1-1478 Grant Title: Simulation of Conformal Spiral Slot Antennas on Composite Platforms Report Title: Report Authors: Primary University Collaborator: Primary NASA-Langley Collaborator: University Address: Date: Robust Development of Hybrid Finite Element Methods For Antennas and Microwave Circuits J. Gong* and J.L. Volakis John L. Volakis Fred Beck Telephone: (804) 864-1829 Radiation Laboratory Department of Electrical Engineering and Computer Science Ann Arbor, Michigan 48109-2122 Email: volakis@umich.edu Phone: (313)764-0500 October 1996 *This report was also submitted by the first author(JG) toward the fulfillment of the requirements for the Ph.D. degree at the University of Michigan.

TABLE OF CONTENTS DEDICATION............................... ii ACKNOWLEDGEMENTS...................................... iii LIST OF FIGURES.............................. ix LIST OF TABLES................................. xv LIST OF APPENDICES............................... xvi CHAPTER I. Introduction................................. 1.1 O verview............................. 1 1.2 Fundamentals of Electromagnetic Theory............ 5 1.2.1 Maxwell 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.................. 36 2.4 Parameter Extraction......................... 38 2.4.1 Radiation and RCS Pattern.............. 38 2.4.2 Gain and Axial Ratio................ 40 vi

III. Edge-Based FE-BI Technique...........;13 3.1 Introduction............................ 13:3.2 Hybrid System Functional...........................-11 3.2.1 FE I Subsystem.................... 45 3.2.2 Boundary Integral Subsystem.................... IS 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 M otivation....................... 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 M odeling............................. 115 vii

7.1 Introductionl........................... I 7.2 Numerical De-embedding............................... 1 7 7.3 Truncatioil Usiilg DNIT................................. 120 7.4 Truncation UsiIlg PMIL................................. 121 7.4.1 Theory.......................... 1'.22 7.4.2 Results....................................... 124 VIII. AWE: Asymptotic Waveform Evaluation............... 133 8.1 Brief Overview of AW E.......................... 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 viii

LIST OF FIGURES Figure 1.1 (a) Recessed cavity in a PEC ground plane. (b) Recessed cavity il a dielectrically coated PEC ground plane................. 12 1.2 Illustrations of equivalence principle when applied to the structure shown in fig. 1.a................................... 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...... 99 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 u0Q 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 Conmparison of tlie conlmputed aIin nmeasured( inplut initpedallce for t 11ie circular patch shtown in fig. 3.6. The feed was placed O.S cm'1 fromll t tie center of the patchl and tile frequency was swept froll 3 to 3.8 GIz. 56 3.8 Illustration of the configuration and mesh of the one-arml colical spiral used for the computation of fig. 3.9.................... 58 3.9 Comparison of the calculated radiation pattern (E,). taken in the 0 = 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

-1.8 ('OIll)IarisoIs of th e nlOIlostatic radar cross sectio Il scattering 1!' an enml)tv aperture Nwith tlhe same cavity size and( dielectric filling (( ) as tie structure shown in fig. 4.6. Agail. tile results were calculated using the regular BiCG F E-BI technique described in clha)ter: 3 aind using the BiC(G-FFT proposed in this clhapter............ 78 4.9 Bistatic RCS scattering by a crcular patch antenna nodeled usillg the regular BiCG FE-BI and the BiCG-FFT algorithm wNith 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 Geomretry 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 X = 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 RadiatioIl Pattern at f=0.944G Hz (lower end of freq(tuency ralige). It can be seen that the axial ratio of the i)atternl b)ecoies larger compared to that at the center frequency. but still remains withill 3dB for a wide angle range. This indicates that the iumlnber of tlie 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 e0 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. 111 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; C=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) Im aginary part............................ 112 xii

6.7 Measured and calculated input illtpe(ldaice for a circullar patcli alntenna having the following specifications: pIatlc ra(tius Ir=2cmIl: sul11,strate thickness dl=0.21814clm: feed location froil ccnlter.rlf=O.i7CIl: r-=9 33: ta/(7= —0.001 measurelent: xxx: this Illetod: o o o: probe modlel. (a) Real part: (b) Imaginary part................ 1 13 7.1 Illustration of a shielded microstrip linle................ 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........................................... 1.23 7.5 Typical field values of TElo mode inside a rectangular waveguide terminated by a perfectly matched uniaxial layer............ 128 7.6 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........................ 128 7.7 Reflection coefficient vs 23t/Ag (c = 1) for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6... 129 7.8 Reflection coefficient vs 2,3t/Aq, with a = /, 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/3t/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 measurement............................... 132 7.13 Comparison of calculated and measured results for the meander line shown in fig.7.12............................... 132 xiii Xlll

8.1 Illustration of the shielded microstrip) still) excited w\itli a cur(rentt probe...................................-10 8.2 Impedance calculations usiIg traditional FENI frequency analysis for a shielded microstrip stub shown in figure 8.1. Solid line is tle real part and the dashed line denotes the imaginary part of tile 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, 6 or z) along with the definition of the tilded param eters in (6.18)........................... 108 xv

LIST OF APPENDICES Appendix A. Evaluation of Matrix Elements for Tetrahedrals............. 1-1 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 developmnent 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

proxinity or aperture couple(I circuit network. etc. anll iiltt'ral '(lalletio Il mt'tlt)(ls are not easily adcaptablle to ilodeling these structures. esplecially ill tlie lpreseice of finitely sized dielectric loadings. Partial differential e(luation ('DIE) technli(llqes (e.g. finite element and finite difference methods) imay also exp)erience dlifficulties 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 sinmulate conformal antennas of arbitrary shape with diverse feeding schemes. \NVith 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 thle geometrical coIlI)lexity Nwitlhout coinllrollisilng accntrac!. IThis hybrid methodology appears to be very attractive for conformlll antennIa Ilode(liit'l. However. its development and application to more pIractical aIld ellmergilng anteIlnlas presents us with many theoretical and numerical challenges. which Awill be exteisi~vel investigated in the work. Specifically, mesh termination plays an important role in FENI 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 sim

-1 I)lified probe feed nlodel fails to accurately predict the iili)lt ii(peda('ce. () t l'he otiler hand. the numerical systenl can becone ill-coii(litioiled wlieii a feed lnetwork is modeled witlhout careful considerations. In tills dissertation various feed m11odels will be investigated in consideration of accuracy and efficiency. They inclule current and voltage gap generators, stripline. microstrip line. coaxial cable, aperture coupled microstrip. etc. In regards to the development and applications of the simulation technicines. 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 proceed 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-BI 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 where specialized truncations suited for guide wave structures are presented. The perfectly matched layer (PML). an

anisotropic artificial ab)sorber used for mesh truIicattion. is investigated iII ter'InII o(f performanace and applications. XN ideband system responses promIpt us to look at morie efficieIlt analysis tools to replace the current brute force frequency domain analysis appI)oaches. (Ch'aptcr 8 discusses a preliminary development of the FEMI 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 FEAM developmnent. 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.

1.2.1 Maxwell Equations Timle-hlarmonic Maxwell equations of differenti al form ill a liIneari. aiisot roic and uniform medium are gi'vei by [10] VxE = -j^'L H-M, (1.1) V xH = jam-E+J; (1.2) V..E = pe (1.3) V..H = (1.4) where E and H are the electric and magnetic field intensity, respectively, M is the radian frequency and the factor ejwt is assumed and suppressed throughout this dissertation; M; and J; 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 Ai represent the permittivity and permeability, respectively, and may be written, in general, as 611 C12 C13 6 = EOCr C= C21 22 623 (1.5) C31 C32 633 UL1 112 813 U = bObr = 10 /121 /22 /123 (1.6) /131 J132 /133 with Eo and /uo being the free space permittivity and permeability.

The procedulre to derive the vector wave eqcuation l)betgiIs by elinliIlat ill oine of the two field quantities from (1.1) and (1.2). To do so. Awe first take a ldot p)rodtllc of (1.1) with the tensor IT- and then take the curl on bothi sides of (1.1) to obtaiI V xM -'E-x E -, x -P x - ~ Mi (1.7) Substitution of (1.2) into (1.7) yields V x (1r V x E) =:2Ooer E-j 'poJi- V x '(1 Mi) or V x (U1 * V x E) - o E = - joJ - Vx ( 1 M (1.8) where ko = wJ/7J1oo is the free space wave number. The dual of (1.8) is given by V x ( V x H)-k 2 H = - woM + V x (1 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 ni points from medium 1 to medium 2. The fields on two sides of the interface satisfy the relation n x (E2 -E) = -Ms (1.10)

S where Ms is a fictitious magnetic surface currellt and El aIl( E.e are tlie electric iie'l inside mnedium 1 and medium 2. respectively. If nmediumll 1 is a perfectlx mallgnelttic conductor (PMC''). then E1 vanishes and (1.10) becomes 11 x E2 - -M,. The surface magnetic current Ms can either be an impressed source (excitation) or may irepl)reseilt a secondary (induced) current. If medium 1 is a perfectly electric conductor (PE(C). n x E2 also vanishes and thus Ms = 0 on the PEC surface. Similarly, for the magnetic field, 7I x (H2 -H) =Js (1.11) where Js denotes an electric surface current. The PEC surface can support electric currents, given by n x H2 = Js, since HI 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 M, = 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 [(v-1. V x E) - (~1. V x E) 1] =-jwJs (1.12) where (/~1. 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 Neuimann I)ound(lar coIndition. by w\\icl a (conlstraaint on the derivative of E at the interface is defined. IThe d(ual of (1.12) is gi\ven by i x [( xH) - r(1 v xH) - (I.I; 2 1 and this condition is used when Nworking with the H field. In many applications. tlie single field formulation is often desired since the system size may b)e kept Iininilmll 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 x nxE+Rn x[H] =0 (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-Yz + FEieY (1.15) H = Hie- - FHiWeIz (1.16) and n x E = -ZH (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 transissionI liIle nlod( of tile guide wave structure. Elilllniiatilln I froIl (1. 1)) and (1.16) in view of (1.17) yields tlce relation i x E- ZH = 2f x E-^' (1.1S) 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 FEMJ 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 'Care 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 (lefine(l plhsical p)rob)leis. The (quest ion t lien arises as to how\ to relate the solutions and hlow mlani colnditions are sufficient to achieve tilt'e correct solution. Uniqueness theorems offer the answer to thlis question. Specifically iII electromagnetics. the E.M solutions are uniquely dete rmined by thl sourcus in a gi e n region plus the tangential components of the (lectric field on boundaries, or plus the tangential components of the magnetic field on boundaries. 2 Equivalence Principle From the uniqueness theorem, an ENM problem can be uniquely solved if tile 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 EM 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 e the icurrent 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].

1 2 equivalent currents cavity ground plane (a) equivalent currents — 40 - 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 neededl and most specifically. tlie coul)liing of tlie fields ill tlie nnIe(I ai(ld outer cavity regions is desired. It is therefore colnvenient to selec tt e total tange'ltial electric or magnetic field to specify thie equivalent currents as M, E x i and J, = i x H. (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 ENM 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') = -I(r - r') (1.20)

1-1 equivalent currents RO (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.2: Illustrations of equivalence principle when applied to the structure shown in fig. l.la.

1 ) equivalent currents RO - so - m (E,H)=(O,0) * ground plane Region of Interest (ROI): Exterior (a) equivalent currents (E,H)=(O,O) cavity CROI Region of Interest (ROI): Interior (b) Figure 1.3: Illustrations of equivalence principle when applied to the structure shown in fig. 1.lb.

. Me Je (a) Me Je (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 M = 0), and I is the idem factor defined as I = xxi + yy + zz. Also, note the identity IJ {P (v x v x Q)- (V x V x P) Q} dV = J [P x V Q+(V xP) xQ] dS (1.21) and upon setting P = H and Q = G, we get JJ {(H (Vx x G)-(V x V x H) G} dV = - * [H xV xG+ (V xH) x ] dS (1.2) From (1.9) and (1.20), the left hand side (LHS) of (1.22) reduces to LHS = -H(r')- JJ V x J G(rlr') dV and the right hand side can be rearranged as RHS I H [7 x V x G] + (V x H) [7 X GI dS JSs

17 Equating the LHS and RiHS yields H(r) = - V' x J G(r'lr)r11 - J {H [i x V x G - (V' x H). [;, x G]} c/S" (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 the entire upper half space. To eliminate the curl on J, we use the dvadic identity. V. (J x G) = (V x J) G - J. (V x G and the divergence theorem to get I7tvoOlume = - JJ J (V x G) dV -- f- * (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) =: - Jf(V'xG)dV' {H- (n x V' x G) + (V' x H- J). (n x G)} dS' = ff J (V x G)dV' - {H. (n x V x G) + jLeoE (n x G)} dS' (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.

IS For ilnstanIce. in the preseIlce of a PE' I)latforml. (1.25 ) is valiId for protiridiIlog coilfigurations as shown in fig. 1.4. In thelse cases. thle surface ilt cgrations are carried out over the platform plus the outer contouirs 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 Intsurface =- ( x G)T ( x H)dS' + x ( x G). (x E)} dS' V= x J/ {H -( x G)dS' il E (1. 29) 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(rr') = Go(rlr') - Go(rlr') + 2izGo(rlr;) (1.27) where Go is the free space Green's function given by -j/k0 Ir-r Go(rlr') -= 471r - r'[

1.9 anid Go(rlr')= (i +,- ) Go(rr') Inserting (1.26) and (1.27) into (1.25). we obtain H(r) Hic(r) + Hrf(r) + 2jkYo i Go(r r') ( (i x E),d' (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 (EMI) 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 -1ill be used whenlever )possible in tliet necessary derivat ions for the FE.M. WNO ith this type of formllulations. isotroI)ic systemis rnmay 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. Thle 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 FEM4. 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.

Considler a t ypI)ical radiation or scatterillg 'prob)lem shown in fig. 2. v1. hlere t lie 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 L given by ~(I = K (2.1) where (I denotes the field E or H, and C = (VX.lVX)-(kcr-) V r Vx - (k2r.) for electric field for magnetic field (2.2) (2.3) Kn is the source term associated with the impressed electric and magnetic currents and may be explicitly given by K, =-jWoJi-V x (- 1 * Mt) Ki = -jwcoMi + V x ( ' Ji) 2 d (=r 1 I for electric field for magnetic field (2.4) (2.5)

As already menlltioned. L is a linear operator aIl(l one call re'adil shlow\ tliat for svmn-etric dielectric tensors At an( T. the pertinent functional of tlie origiiial PDI)1 has the form T(q)) = - < ). ~Lb > - < D.Ki > ('2.6) where the inner product <, > is defined as < A,B>= A- Bd (2.7) JQ (with B* being the complex conjugate of B) for lossless media, or more generally as < A,B >= A BV (B.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 6.(~) = 0 (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.

21 IThe existence of the functional (2.6) requires the operator LC be slf (I djoit.nt a property usually defined to satisfy < Lb. T >=< I. > (2.10) where 4 and AI represent any two admissible functions. If (2.10) hlolds. not only does the numerical system derived from the functional (2.6) remain symmniletric. 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 re)ports wNere also seenI in later papers [IS-231] for tlie probleiIs of wa\ve propagation inside anisotropic Inmedia. No radiation and scattering aiial\sis lias beell reported iiI this context. In what follows, we generalize the FEMNI forimulation to lossv andI anlisotropic electromagnetic problems. Specifically, we show that two different methods. one withl the aid of an adjoint auxiliary system and the other with the Lagrange multiplier caI be used to construct the pertinent functional. Galerkin's method is then conipare(l 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-Hermit ian 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 F =< ~L, > - <,Ka > - < T,K > (2.11) where D is the unknown solution function of the original PDE problem and K is the right-hand-side function as in (2.1). Similarly, I is the solution function of the adjoint PDE such that La = Ka (2.12) where Ca can be derived from < ~:, >=< (,DLa > (2.13) with Ca 4. It should be remarked that the functional (2.11) reduces to (2.6) (except for a possible constant coefficient) if L is self-adjoint. Also, the original PDE and its

ad(joint countIerl)art (2.12) cain be reco\vere( througlIl tlie \-ariatioiial )I'(('cess wit I respect to thle functions 4) and 4P. respectively (tili )le sili rivatioIl is ollitted here). After discretization is carried out. tlie filial numerical sy'stcIll is sVylIelltr'ic. This can be shown as follows. Let ) =.Vi. = yjV (2.11) i j where Vi is the basis function used for both unknown functions and.r;. y; are the corresponding expansion coefficients. Inserting (2.14) into (2.11) yields "= xiyj < Vi,Vj > - < ViKa>- y < Vj,K > (2.15) I 3 I 3 i j i j Upon performing the differentiation with respect to xi and yj individually, we get the two decoupled systems of linear equations - (2.16) where the matrices Qx, QY and the column vectors Kx, KY are given by 2j = <V,, V > Qj = < V, Vj > Kx = < V,K > Kyi = < V,Ka> In general, Q?. Q-, Qj + Qi and Qi 7 Qi. However, Qi = Q, (QOy)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.

2.1.2 Pertinent Functional For Lossy/Anisotropic Media II An alternative to using an adjoiIlt systeni is to emp)loy thle Lagraiioe imultiiliier technique in constructing the pertinent functional. The Lagrange Illmltil)lier is tlsllally used to incorporate additional constraints to a system. To illustrate the techniquiiCe. consider the same PDE model as described in (2.1). and we first rewrite it as I - Ki = 0 (2.17) Next, we assume an expansion space where the solution is defined and solved. The unknown function 4) 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 FT((, A) =< A, L-K > (2.18) This functional is now used to formulate the FEM. As described above, on applying the Rayleigh-Ritz procedure to both O and A using the same set of basis functions, viz. = xiVi, A= yjV (2.19) i j we obtain F(b, A) = E E xiyj < V, LVi -Ki > (2.20) i j Carrying out differentiation with respect to xi and yj, individually, yields (o ) () ( ) (2 21) where Qx = QC" 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

similar to those in (2.16). We further observe tlhat tile Lagracnge Im Iultiplier techllni(Ite may be regarded as a special case to (2.11) where a hoIlloogeleous adjoillt PDl) is now virtually assumed (i.e. Ka = 0). Again. it should be remarked tlhat wlheI a self-adjoint problem is considered. the above formulations (both the adjoinllt systemll 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 4 expanded using the same basis fli(unctions. 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. Traditionally, 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 <,~lb >=< T,Ki > (2.22) where 1 and I are both defined in the same function space. Specifically, in Galerkin's method one seeks the solution for the unknown function (I which satisfies certain prescribed constraints and (2.22), with the aid of another arbitrarily chosen function XV.

}29 Similarly to the variational approach. the RaVlcighl -Ritz procedure l11ma also() used to project the continuous space onto a finite discrete separabtle 1lill)ert space. The mathematical problem is theni rephrased to seek a iscrete solution set whlose entries are the coefficients of the expansion. The testing functionIl Illust. 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 j = ( results in a symmetric numerical system. Otherwise, no matter howthe 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 XiYJ < Vj,, V >= yi < Vj,K, > i j j or E3y { i < Vj,CV> - < Vj, K = 0 (2.23) As assumed, 1 is an arbitrary function. Thus the term in the curly bracket should vanish, yielding Qxx = Kx (2.24) which is exactly the same as the subsystem derived from the variational approach.

It is noted t hat if one chose the Lagrailge imultip)lier (inll tlte \variatiolial approaIch) as the testing function TJ. the samne numerical systell would result. Alalti(cally. though. the entire system is thwice the size of that derived via GCalerkinls mlethod. 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 il fig. 2.2. where a perfectly conducting electric (PEC) body is coated with a, dielectric layer whose relative permeability and permittivity are 1d and Ed, respectively. (Note that, r-ei rp ` rf Qd: dielectric coated region (Ed,T d) ff: free space (f = -/f = 1) fa: absorbing layer (ea, ia) 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 Fo: 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

: 1 with this t ype of probllems hlave not beenI carefullyv addressel in tile litIerat urc. 011e of them is the equivalence between the xariational and the Galerkinls Iletliod \\wlieii the scattered field is used as the working variablle. Proof of this ecluivalence cani b)e tedious and cumbersome but it is nevertheless an important issue. I Infortunately. 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~t as Htot = Hscat + Hinc (2.25) where HSCat and Hinc are the scattered and incident field, respectively. Next, weighting the source free wave equation with the testing function V yields JV { x e V x H- kor H} dQ = 0 (2.26)

where = Qdj + 9.f + Qa encompasses the entire coputational re'gionI. T proceed with the derivation of the weak formi wave e(quation. it is necessary to introdice certain constraints on the scattered and incident fields within Q9 and on thie bol(lundary. First, since the incident field is not allowed to pass through the albsorber layer as well as the metal back wall Fo, we note that Htot(r) H t + Hi r Qd + Qf(2.27) Hscat r C otherwise with the incident field satisfying x -0 r e Qa + Qf {Vx.-Vx-k} Hinc{ (2.28) | 0 r C 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 Hscat 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 = Escat + 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 Hinc 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 H'nd indeed satisfied Mlaxwell ecluations in Qu. a contradictorv bouiidair collnditionl would immediately result. This canl be seen by taking a curl op)eration on b)othl sides of (2.25). In view of (2.29) and Mlaxwell equations in dielectric media. we woul(d iha\e d1. xH = f. xH'nc + d v x Hs r d (2.31) Imposing the condition of total field tangential continuity at the boundary Ed woulld yield i X { dI1 V x HsctIF+ - rf V x HsC 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 n x {Ed v ' X H >ar+ - ff. V x Hsct} (2 x {33)=-1 } i - = -f x.VxH (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 media. Mathematically, this implies the condition jw'oEin = ed. V x Hin r fQd (2.34)

:3-1 Consistentl one can derive the other boulndarv colnditions required for tlie H1 F:lI formulation following thie same procedure. They are classified as Diricllet aIndl Neiimann conditions as follows. * Dirichlet Conditions Boundaries Conditions Fp n x Hsa = K - h x Hin (K, unknown current) Fd n x {(Hcat + - Hscat|} 0 rf x {HscatI+ - Hscat = -n x Hinc Fo n x Hcat = Kb (Ks unknown current) * Neumann Conditions Boundaries Conditions 1 V Hscat = - 1 inc Fp xEd V x H-x } -V Xd H Ff x {1 V x Hscat } n X:f V x Hinc r AX(V - V H = o x —1 a Vx Hscat = where {} |+ denotes {} +- {} - and r and Ha"" positive and negative sides of a specific boundary. should take the values at 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 = 0 (2.36)

where Iltd = [ XV V (* -- V X(Hsca + H") - AkV 7 (Hs-st + H")1 (/ dd ( + jV- 0 x! ~ V x H) S' + V. (i x ~ x H'') dS' p d + V ( l x ~ V x Hscat) ds (2.37) d Intf = j V x V.1 V x Hscat dQ + j V (-nhl x i1. x HsCt) dS + V. (2 x1 V x Hst) d5' (2.;38) rf Inta = V x V 1 V x Hsca dQV+ v ( X. * V X Hx) ds + V. (^ x1 V x.Hst) dS (2.:39) with no, hi, n2 and h3 being the unit normals at the boundaries Fp, Fd, Ff 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 Jf j V J2 x V:f v x Hinc dS (2.40) Similarly, the sum of the surface integrals (involving HWCt) on Fd in (2.37) and (2.38) becomes xV *1) V x HtCdS (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

836 I)erneablility. TIhis observation does Iot hold for (2.40). \\wich arises fromll tliet illlerent incident field definition on the two sides of the b)oundary I f. as Ilmat leatll callvi shown in (2.27). Furthermore. the second term il (2.41) w\ill be canlcellcd out witlh the integral on Fd involving Hizf' in (2.37). After a sinlmple manipulationl. the final system reduces to I (v x V*;-1 Vx H cat- kV,l Hscat) dQ XJd (v x vcr v x Hif kV l Hnc) dQ V x V C,1 - V x H - k ' d Q d - V * n1 X Ef1 V X HiCd + J V *n2 x 71 V x Hicd (2.42) 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 J1(H)VxH (.VxHkH H d= ( V x H H d (2.43) where, as before, r and jI denote the corresponding relative permittivity and permebility, respectively. Also similar to Galerkin's method, one would now logically proceed with the field decompostion H = Hsat + Hinc 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 7ihe discrep)ancy arises fromI thle assumIption of tlie fulIctioilal (2.-43). whlich is nlot a valid expression when the corresponding PDI)E operator is no longer self-adjoilnt. lRecall in section 2.1.1 that the presence of genleral anisotropic and(l lossy Iel(dlia reI(sliirs thle application n of aauxilar adjoint systeni. together withi which the pertinnllt functional is given by a generalized form in (2.11). It is therefore necessary to begill with this generalized functional rather than (2.43). In a source-free region. (2.11) explicitly takes the form F(Ha, H)= Ha v x = x H - k0 H dV (2.44) and the variation is imposed on the adjoint variable Ha to get 6F(Ha,H) = 0 (2.45) 6H=o 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).

:3is 2.4 Parameter Extraction An accurate full wave analysis can only predict the near field (for I 1)DE t ype methods) or current (for integral based techniques) dlistriblutiolns. whichl can )e used to obtain certain practical parameters depending oni applications. For intricate svsterns, involved numerical models may be needed for output data extraction. includiong 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 in terms of the dyadic Green's function must be used to consider the equivalent currents and the free 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 6). 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

as M = -- x E' whlerle EB is tie electric field e\valuated oil tliet a)tertllrt andl is normal to thle a)ertture surface. Thle electric vector potelitial is thleIl gixvel I)v F(0. o) = M- fM dS r477 r 4 7, r W = (oJk J(Es x ^)eJkr'dS" (2.46) where the electric field is typically expanded in terms of the surface vector b)asis function S'. Introducing this expansion. (2.46) becomes FOQ=) =(4 k E JJ s xfjk. dS' O-jkr - 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 4 projections are given by oe r Ho = -jw _4 {cos 0 cos OV, + cos 0 sin /AVy - sin OVz} 47rr 47 ' (2.49) H = -jw --- {- sin 4V, + cos <,Vy} _ 0~-3j kr 47,rP in which Po = -j {'} and Po = -j {}, and {.} stands for the corresponding terms in the curly brackets. The RCS of the t (t = 0 or 0) component can be represented in terms of Po and P, to yield rcs 4r2 lpHt Pt2 P12 (2.50) lt 4 -, 4, At20

-10 w\Nhere A and =,o/ut are the free space \wav\elenogt h and thle intrinsic cimpedance. respectively. In (2.50). the incident wave jil' was set to unit\ \v wit hout loss of geinerality. In practice. it is customary to express RC'S in (lB. w\Nhere (a is usually norIlalizedl 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 arcs to avoid a repetition of post processing. Specifically, the formula rad= t (2.51) (arcs)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 100o efficiency) is given by 27r Umax G = (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.)

-11 Since the laxilmulni powNer density ( niax canl also be explressed as Unax = - (o~ + O7o) ('2.. )) 4w the antenna gain becomes - Zo(o9 + o) (2.)-1 2Prad Thus, the computation of G is rather trivially done once ua 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 = I2Ri,, where I is the known current source on the feed, and Rjn 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 = IJ U dQ = ( J + (o+ ) dQ (2.55) over the half space. It is obvious that if a certain symmetry of the pattern remains such as circular about the vertical Z axis, then U(O, () reduces to U(0) and the above

-1'2 integration canl be effectively evaluated. Othlerwise. a Inullltrical '- 1) ilntegrationI is required. Axial ratio (AR) is another impiortant antenna paranmeter. especiallvy wNheni a circular polarization (CP) of antenna's performance is of the primlary conIceirn. SiIlce AR also features the far zone effect of the antenna, it is desirable to determine A.-? with a minimum computational load. It is noted that given the above pre-calculate(l oa and o6. one is again able to determine AR uniquely by AR n ('2l(i) ' short with 1 { 2 + 2 + [4 + g + 2 2 COS 3} Viong gV/ -{-70 -- (2.57)3 r 21 {a2 + a2 - [a4 + a + 2 gu cos d] } I short 0 {c1o+Udj where f = 2((IjH - (H,), 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 = -jm{ln P (2.58) I 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

-1-1 to extend the met hod to antennas of arl)it rar a sllap)e. IIn this chapt er. an edge-blased hybrid finite elenment-boundaru iIntegral formulation is pIresenlt ed for tlie chlaracterization of arbitrarily shaped cavity-backed antennas [29]. An exalmle of sluch 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 CPWN 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 pertinent to the scattering and radiation by a cavity-backed configuration shown in fig. 3.1 may be written as

-15 z A Patch Ground plane Aperture CavitCoax cable opening. (surface C) X - Coax cable Figure 3.1: Illustration of a typical radiation and scattering problem. F(E) = {(V xE) -(VxE)- kCrE.E dv 2 JJJv I hr + 2jkoZo j(E x Hi) *^dS + E. (jkoZoJZ + V x — M) dv (3.1) - 2ko (E x ).J (E x ). I+ -VV) Go(r, r) dS' dS where Ji and M, represent interior electric and magnetic current sources within the cavity V; HI 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; Cr and jur 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)

-1 6C whlere PLu (lelotes the volumie intecrral contribuitionis anld simi1larly I.~,; accounlts for t lie surface integral contributions. Thec cavity volumec is sull)(livide(l Into N eletrahiedial elements V' (e= 1.2...). and withiln each tetrahiedron as sliown In fig. 3.2 t lie fieldl nodes/vertices Figure 3.2: A tetrahedron and its local node/edge numbering scheme is expanded using edge-based elements as E z:[V]Tj{Ele (3.3) wi'th [Ve = If{Vx}{K}{~fVz}]e vu 2 {VU} =, u xIyIz (3.4) vu 6 ~2 {Eje ~6 / e

-1 in which 'I U is the u ( u = x. y or - ) component of th le volumIle vect or basis funct ions along thle ith edge. The unknown vector {E}l hlas six entries. one for eac1 ttetralledron edge. (Here we use square brackets for matrices and curly brackets for vectors). Inserting (3.3) into (3.1). and taking the first -ariation of F1 with respect to {E}. yields 6SF, = {[A4]{E}Y + {I}e} (3.5) where [A T - 2 {k[ c[V][ } &dv (3.6) {K}e J [V]e jkoo + V x Mi, y dv (3.7) a {; - a ay az [DV]T = a (3.8) 0 0O a {Vy} -a f{V} 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.

3.2.2 Boundary Integral Subsystem To discretize the surface integrals in (3.1). tlle apert ure is stll)(divide( inlto triallgular elements since these correspond to the faces of the tetrahledrals. \\ ithin each triangle, the field is represented as E- [S]{Esc (3.9) where [S]e = [{S}{2}]e. {S} = S2, t= x,y (3.10) Su3 {Es}e = Es2 Es3 e in which Sui 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}e, we obtain 6Fs -= {[B]e{Es}e + {Le}j (3.11) e where [B]e = - {2k[Se] [Se] + 2 [{s} - {s}] [a{sY}- a- {sx}T]} *Go(r, r') dS dS' (3.12) and {L}e = j2koZo j/s [Se] ( H' dS (3.13) 77se v -^

-19 Note that in (3.12) the elements of thie array [.s] are funlctions of tlie observatiloll vector r. whereas the elements of [,']T are with resI)ect to tile inltegrationl poiilt r'. A suitable set of linear edge-based surface basis functions is - x (r- ri)(r) r E Sc6 S,(r) = 2. (3.14) 0 otherwise In this expression (referring to fig. 3.3). 1- denotes the length of the ith edge and r, is ith edge S 0 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 rE S; where Se = S + + S. The constant Ae 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 Si(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 comIponlents cwe conlbine (3.5) and (3.11). and after assembly we obtain the systell {[A]{E} + {}} )+ {[B]{Es} + {L} = 0 (3.16) In this. {K} 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 {E,} 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 [B] 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

5 1 this scliheme. tile combinationi of the two matrices is perforilied only onIce oltside t lIe iteration. However, the second approach is compatil)le wit itlhe BiCG-FFT1' scleIIe. where the FFT algorithm is employed to exploit the convolnttionial pro)perty of tlie integral operator, thus eliminating a need to explicitly store thle 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 ca,vity-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 cavity-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;

5)7 ((d) tihe nodes identifying metallic boullndaries. antennIla patlch(es). feed(s). and(l possible vertical posts. For arbitrary antenna geometries. it is necessary to emplllo a sophlisticated vollillme 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. \Ve 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

IMPLEMENTATION FLOW CHART MESHER,?T Universal File re —v - --------- - / Pre-Processors/ FOR FEM System BI System Feed Models Combinations - FEM System ----------- -------- BI System Core Program / Feed Models BiCG Solver/ I L Combinations -----------— 4 ------ 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.

exterior regions. A significant effort was (l\evoted to (develop)illg a programll ill such a manner so tlhat both the storage and computational req(uirenllelts caln le in(llized(l. Specifically, the boundary conditions on the metallic surfaces are eIlnforce(l a pri ori to condense the system which involves only nonzero field comp)onents. To this point. the sparse finite element matrix was stored as a single array of length X\,A\Z. where N\'V is the total number of unknowns within the cavity volume and N\l 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 NS x Ns BI integral matrix is added to the FE array resulting in a 1-D array about NvNn +N NA 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 BI 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 oni the surface of a 0.406- cI11 thick substrate having a relative dielectric constant of tr = 2.9. The patcll's dialieter is 2.6 cm and the substrate is enclosed in a circular cavity 6.292 cm wide. Tllis cavity and the patch were recessed in a low cross section body for mneasuring its RCS. A comparison of the measured and calculated backscatter ar00 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-13I 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, X, z) denote the standard cylindrical coordinates, a~ are equal to 0.0832A and 0.0257A, respectively, and 0 < 5 < 27r. 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 Er principal plane radiation pattern taken in the X = 90~-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

u.uuu E CQ an U c 0 1C) -o mu m, OQJ -10.0 -20.0 -30.0 -40.0 -50.0 -60.0 L 4.00 5.00 6.00 7.00 8.00 frequency (GHz) Figure 3.6: Comparison of the computed and measured 600 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.

data given in [35]. As cain be expected. thle Ei, patterIl (nlot shIownI) differted from the mneasured data near the horizon because of interference fromi thle finite circuilar cavity housing the spiral which was included in the analytical mnodel. -lie latter w\as not part of the measurement configuration which consisted of the spiral alntennlla oni 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 or = 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 demonstrate 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

Figure 3.8: Illustration of the configuration and used for the computation of fig. 3.9. mesh of the one-arm conical spiral cS E 4-A o c ~W,., -10. -20. -30. -40. -50....., I *,.,.. I,.,.. W J I TI I I 1 I V I I. I I V I I * * Et (FE-BT) E Et measured. _... 1.. 1.. I........... -90. -6. 0. - 0. 3... 90. -90. -60. -30. 0. 30. 60. 90. 8 Figure 3.9: Comparison of the calculated radiation pattern (E~), taken in the k = 90~-plane, with data in reference [35] for the one-arm conical spiral shown in fig. 3.8

0 If 0 0 f 0 0 1 4 I I % \ probe, t r r f f a= 12.35 cm b = 0.75 cm Po=7.7 cm 0.7 <f< 1 GHz d = 3 cm I 0.5.- Calculated / 0 FE-BI Melhod / re, E = 1.35 ~. I. ~. -' 1.0 Figure 3.10: Comparison of input impedance calculations for the illustrated cavitybacked slot.

6(i) rectangular p)atcdies l)ecause tiley occuIpy a small area whenli opertated at tliet same frecluency [37]. I nIfortunately. no sufficient researchi oin tllis geonlMett r ias een repIorlted in the literature due to a lack of aInal\sis 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 antennla 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. --— lracatn 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

G2 pIroducts via the discrete Fourier transform ()DFT ) anld av\oiding a Ileed to stoic tilie full BI matrix. This memory savingi schemIe lias alrealdy beei app)lied to 11 solut ions involving rectangular grids [39.34]. and a similar impllementation was also rel)orted for triangular surface grids [40] involving inherent a)proximations. In this clhater. we first show that the BiCG-FFT solver can be precisely implemented on uniforml 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 computational problems which may be impossible for direct system solvers. It is especially desirable to employ the algorithms when one seeks the solutions of large-scale 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.

6:3 As suggested in [45.46]. the algoritlill ulse(l in this work is as follows: Given pi = rl = b - A xi P1 - rl For k= 1,2.3.... r rk P 'A Pk rk+l = rk - akA Pk rk+l = rk - A'T k Pk rk+l1 rk+l rk - *rk Pk+l = rk + /kPk Pk+l = rk+l /3kPk Xk+1 = Xk + AckPk 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 = r, 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,

(;-1 the algorithm minimizes the function f(x) = x A x - b x Hence x obtained from the CG algorithm after the k iteration steps b)ecomes 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 0(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 accomplished 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.

As described in the next section. tlhe bouIIndary illtegral eq(lation can (be cast inllto a convolutional form if a uniform grid is applied for discretizatioll. Iis is Inot surprising since the Green's function is involved in the integratioI. To solve tie eqluatiol using the CG algorithm, it remains to carry out the conxolution at eachl 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 O(N2) CPU requirement for a matrix-vector product in a traditional fashion, the scheme needs O(N log2 N) 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. WAe 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 (m, 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 m = 0,1,2,... M n = 0,1,2,...,P

660 z N I 0 m: Figure 4.1: 0 1 2 3... Structured mesh consists am a M of equal right triangles with i = for the y-directed edges, i = 2 for the diagonal edges and i = 3 for the x-directed edges. Consequently, we find that Mi= M-2 i= 1 M-1 i= 2 M-1 i=3 Ni= 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

(i7 withl each of tlle aforementioiied class of ed(ges can }e rcwrit tell as (,A-y.- /)- + (x -,,A.r)l (a.. y.) C + ^n = - (y-( + 1 )Ay).I + ((, + 2)A. \-.r)J (.r. y) St- (4 I0 ot lierwise (nAy - y)- + (x. - (77I + l)Ax) C "'5 SM'(x.y) (Ax + (Ay (y - (1 + 1)Ay)+ - + (Ax- (4.3) 0 othlerwise ((n + 2)Ay - y)x' + (x - (m + 1)A.x) (x. y) C 5+ S.n(,y) Xy (y- nAy) + (mA- x)y) (-. y) C SX (4.4) 10 otherwise where the superscripts refer to the edge class. After the discretization and assenmbly 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 Mi NJ {BI subsystem} = [B] {E,} S= E E mB.nm,E (4.5) j=l m'-0 n'=O in which (in, 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 = nMi + M of the column resulting after the execution of the boundary matrix-vector product. It is readily found that B2 = 2k2J Jj S S, Go(r, r') dx dy dx' dy' (4.6) + (A8 xy)2 fj fj (r) ejr) lil Go(r, r') dx dy dx' dy' $ /xyL/L 3

wit I _A i = 1 I=i v()2+(A )2 1 2i=2 (4.7) A _x i = 3 More importantly, it can be shown that the BI subsystem (4.6) exhibits the convolutional property BJn,, = B _J _ ) and thus we can rewrite (4.5) as 3 [B]{ES)} B'I EJ (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 BJ,, nn, ) is also invoked, implying B" I B) - (4.9) (mr-m',n-n') (= m'-m, n'-n) it is concluded that the total non-redundant entries in the BI matrix are 3 3 Np= Ni(Mi + Mj- 1) (4.10) i=1 j=l This should be compared to the (E,1 Mi N) 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 B13 -rn, ') = B(Fiin, i) be cast in a 2-D array which has the usual periodic form, (m-m',n-n'} P 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 MFTxNFT

(;9 array B o< ii <.\rT<I 0 < T7 < XI B'J(-/. 0 <). B B"(,- 1-NFT). ' mFT - B (m -l-MIFT, - - NF 0 otherwise with the corresponding field vector given by 1 < n < SIFT i < 11P' - Nj + 1 < n < NFT MrIFT - AI + 1 < M < MIFT NFT- -1E + 1< < n < NFT (4.11) E-( E(im,n) 0 < m < M3, 0 < n < N3 Ep(,) = - (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} = E S {DFT-1 {DFT{B} DFT{Ep}}} i=l (4.13) 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 mn, = B1S _-, n_) 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 hell) to prove the convolutional property of the boundary integral. (y-directed) in S+ for the trial function and the i = 2 class edges (diagonal) il Sfor the testing function will be used. To this point, substitutng S1 j(x, y) in S+ anc S2,, (x y) in S- into the first term in (4.6) yields Intmn,mln' CJJ+ JJC [(n - y)(n'a - y )+ (x - m x)(m x + 2A x- x')] +y 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 [m\x, (m + l)Ax] and [nAy, (n + 1),Ay] for the unprimed coordinates and similarly for the primed ones. Therefore, (4.14) will be simplified if the following transforms are introduced, viz. x = mAx + x' = max +,' (4.15) y = nay + y' = nay + r Indeed, on substituting the transforms, one obtains Intmnm'n' = C ((A - y ()]A Go {) [(7m - ')Ax+ G[( - ')3 +[(- '))Ay + (rd - 1x')]7} ddry d'd1 (4.16) (^- i')] + y[(n - n')Ay + (ri - ij')]} d^dr d^'dri' (4.16)

71l Altlhoiilgh this integral imay not be exp)ressible inll tersll of an elielletarv fullictionI. 01o is however assure(d that the resulting expressiol Illmust be a fullctioll of [(/l - ll'). (n - 7')] no matter w*hat form the Green's functioin Go awill take. Matihemiaticall\y, it mieans Intmnm'n' = Iltm-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 gridding 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 cavity recessed patch (a) Circular Patch Recessed Underneath the Ground Plane 30..... 'l..... ' l.|..............I "0 E a-) LZ C". Q V= u u: cil M3 CA 20. 10. 0. -10. -20. -30. -40. -50. 7 — liyl- -V ---- 1 __,I Id I 1. I II:.......I permittivity = 4 -- - permittivity = 1 A permittivity = 4 (with FFT) v permittivity = 1 (with FFT) -60. ''''"'''.......... 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 6 (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.

: () ITV () 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. That is, it is necessary to append to the system (3.11) the relations {EsU = [TF]{JES}U {EnU = [TB] {Es}u (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

71 /N eu overlaid unstructured ' mesh triangles rend2 mid - (m,n) uniform grid edge 0 \ / \ / 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 (Eu)(m,n) = { 2Nn E Eu (rendi) 2/~endy n Nmid + S- Enu(rmid) Nmid k=l + 2N Nend2 E u(rend2) (4.19) 2Nend2 k-i in which eu denotes the unit vector along x, y or the diagonal, depending on the class of edge being considered. The quantities Ek represent the fields in the non-uniform grid triangles with the superscript k being a sum variable in case rend1, rend2 or rmid specify a point shared by more than one triangle. Obviously, Nendi, Nmid and Nend2 denote the number of non-uniform grid triangles sharing the node at rend, 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

transformation miiatrix are given bv -\', nd d,3 1 f Iin, -I (TF)i = - ) - (I.S (rend,) -end k=!=1 - mid 3 + - sijS (rmid) - mnid k k=l t=l e nd, 3 + - E E Send (rend (4.20) 2 A end2 =l =S 1 fend2 k=l f=1 in which f 1 j =J6ij = s j 0 otherwise and the global indices i and j correspond to the ith uniform grid edge and the jth nonuniform grid edge. The subscript je 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 Se. 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 patch 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 formulation. 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 maximum non-zero entries in each row of [TB] will be 15, but the typical number is much less.

T7( 4.4 Results Figure 4.6 shows a cavity-backed 2 x 2 patch array. where eaclh patcli is a righlt angled triangle. Since this geometry is adaptable to a uniforll nmesll with rigllt angled triangles, it is used to verify our proposed FE-BI technique incorporatedc with the BiCG-FFT system solver. The developed FE-BI code with the BiCG-FFT is first compared with the original version of the h-ybrid 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 = O 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 cr = 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.

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.

7S 10. C4 CJ -S.s E 9 cam CE 6 Pi ar m -9 0. -10. -20. -30. -40. -50. -60. -70. -80. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 0 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. -10.0 -20.0 -30.0 -40.0 -50.0. -60.0 -70.0 -800 - (X) 10.0() 20.0) 30.00 40.(X) 50.(X) 60.00 70.(X) 80.t0 t,)()H Figure 4.8: Comparisons of the monostatic radar cross section scattering by an empty aperture with the same cavity size and dielectric filling (r=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. - 7 C 4 C-, a) L. VT) u -16. -24. -32. -40. -48. -56. -64. -72. -80. v 7v7vvx — R-7-t- Ad ---- to > - With No Transform (permittivity= 10) i With No Transform (permittivity=4) A With Transform (permittivity= 10) v With Transform (permittivity=4)..I...............I.........I........... 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. e Figure 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.

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

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 is recessed in a ground plane. To solve for the E-field inside and on the aperture of the

cavit X. a staiidard aIppIroacli is to extreillize the fulct iolial (:3. 1) Nwhich. tfor radliat ion aIld scattering problems. may be generalized as F(E) = j {(Vx E) -sr (V x E) - -2E r E}dVl + JJI E (jkoZoJi + V x 1 M) dl + jkoZo fJ E (H x) dS (1) So+Sf where or and ir denote the relative tensor constitutive parameters of the cavitmedium, 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, I' is the volume occupied by the source(s) (J3 or MZ) 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 + 2jkolY JJ Go(r, r') (z 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).

For t hie antenna problem shownI in fig. 5.1 where Hlie platformi is a p)lalnar "IroU1i(l plane. G becomes the half space dyadic Green's funIctionI /- -+k1 \ jr-r'l G I + VI (5-3) k 2-o 4/> Jr - r'l with r and r' being the observation and integration points respectively. and I = xr + yy + ' is the unit dyad. In connection with our problem, i.e. that of a cavity recessed in a ground plane. H9~ 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 = EjeV = [V]T{ Ee} (5.4) j=l where [V]e = [{1}, {Vy}, {Vz}], and {Ee} = {Ele, E,..., Eg}'. The vectors {1V}j, 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 V are chosen to be edgebased functions, the unknown coefficients E; represent the average field along the jth edge of the ethl element. A corresponding representation for the aperture fields is 3 E(r) = E' S(r) = [Sf{S}, (5.5) i=1 where [S]s = [S,-, Sy] and V(r) reduces to Ss(r) when the position vector is on the slot. To generate the discrete system for E;, (5.4) and (5.5) are substituted into (5.1) and subsequently F(E) is differentiated with respect to each unknown EJ. With the

S-i Figure 5.2: Illustration of tessellation using prisms A z 0 80 0' Z=Z c 6 i=1,2,3 j=4,5,6 k=7,8,9 Figure 5.3: Right angled prism

understaan(Iing that the surface field coefficients 1s are a sublset of 1. we obtainl {[F3} =.-]E + Z[Bs]{Es + \' + 5. =;) -e=l (=1 e=l s where the sums are over the total number of volume or surface elements. IIn this. tihe matrix elements are given by = j {V x V. V x V - V V} d (5.7) '< = jkZo /fJJ vi [J+ V x 71 x M'] dt (5.8) B9 = - 2 k S(r) S (r')Go(r, r') ds ds' +2J Jj [v x S (r)]z[V' x S'(r')]zGo(r,r') dsds' (5.9) Ls = 2jkoZo Ss. (H x Z)ds (5.10) where the subscript z in B.j denotes to take the z component. It is noted that Ls is removed in case of radiation problems and that the same holds for KI 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 1i h- = e z x (r-r,), S- = hi (5.11) where ri denotes the location of the ith node, ei is the unit vector along the ith triangular edge opposite to the ith node, li denotes the length of this edge and r is any position vector terminated inside the triangle. One way to obtain an edgebased field representation for the prism is to utilize the nodal basis functions [52]

SG and thlen apply)I the procedure discussed ill [-9.5.5. ]. Towever'. all alterIlative aind more )lphysically meaningful aIpproacll calln e emplloyed for t lie coiist ructio1 of tlie edge elements. Referring to fig. 5.2. it is evident that if r is ill th x —y I)Ianc. thel 5'fc in (5.11) gives the area of another triangle 12'3' such that the lengthls of edges joining the nodes 2 - 3 and 2' - 3' are equal. \Nith this definition of r. we define a vector a vector S = x (r-ri) (5.12) 69 where i * S = C. That is. the vector component along Ce 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. {Z - Az) Vi =z- Z Si i=1,2,3 &\z -. j 4,5,6 (5.13) Vj ( -) j=4,5,6 (5.13) Az Vk =,k k = 7,8,9 where (k is the triangle simplex coordinate associated with the kth prism vertex at (xk, Yk). As illustrated in fig. 5.3, z, and h = Az represent the offset coordinate and 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].

S 5.4 Applications Thin slot antenna structures have been treated usinlg the above fornimulated E — 13I 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. Moreover. 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

ss probe /', \, a=12.35 cm \, b=0.75 cm po= 7.7 cm a 0.7 < f <1 GHz d=3cm I E.tF1.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 aCT2 Es _R 1 -ctR

I1(). I V I T I 0 c -10 V (I, CD 0 Ca a. Ca -20 -30 -40 ' -50 -60 - -70 -8I I-I --- I I J --- —-- I I -100 -80 -60 -40 -20 0 20 theta (degree) 40 60 80 100 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 4 = 0 plane and the observation cut is perpendicular to that plane. Radiation: X-pol and Co-pol radiation patterns in the 4 = 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.

90 t Es radome cover (on c=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 7' is the transmission coefficient of tile FSS. R = 1 - 7 ailld o is tile reflection coefficient associated with the cavity base. To reduce tlie Illlltip)le iInteracti Ions within the cavity. it is appropriate to terminate the cavity witli some absorber. tllhus reducing the value of a to less than 0.1. Since R is also small near resonance. a good approximation for T is T() 10 log E and upon considering the next higher order cavity interactions, we have TdB T) + 0 log [1 - ao( - r))] 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 FEl-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 PMM Ca.c...,o - ^ '......... Measuren:' ~*.- FE-BI Et - -, — -...................,.. [........ 0 / - -15......... 30. -20....................... -25.......,J I 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 + 3, 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%o, 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 Er = 7 7r -_Er./~~~.....- -I a a 0 (n Co Vn-155 E cn ca PO 1I Figure 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

9-1 spiral shape OIl its p)erformance. A benchmark test model is designated to opIerate from 118 NMH11z to 157 M111z. to replace the conventional protruding blade antenna. The size however is mIllucl compact with its conformality property. Our simulation model is scaled by 1/S to operate at 944MXIHz to 1256MIHz with the center frequency 1100M1Hz. The values of a and /3 above were determined based on this frequency band and also the numlber of turns (4.5). The cavity is filled with a dielectric slab ( 2r = 2.2) of 0.3 cm depth, corresponding to approximately 0.011 free space wavelength at the center frequen cy. 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 ( = 900 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. At 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 wide azimuthal angle (as wide as 60~ in the optimum case). In practice, we notice that absorbing materials may be needed to

9.5 regulate the magnetic currents at tlle beginnling or elnds of tllie slot spiral. 'esI)'iallv wnhen the number of turns is minimized. Frequency (GHz) 0.944 1.100 1.256 Gain (dB) 7.22 6.66 5.2' 3 Axial Ratio (dB) 2.7 1.0 3 Table 5.1: Comparisons of gain and axial ratio at different operating frequencies 5.5 Concluding Remarks A hybrid finite element-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.

Figure 5.9: Illustration of a typical 2-arm slot-spiral design frequency=1.1 GHz / -20. \ 1 ".. "; '.I" -30.-' solid line: Ephi, dashed line: E_theta Figure 5.10: Radiation Pattern at f=l.lGHz (center frequency ratio is achieved up to 60~ degree. design). A good axial

97 frequency=0.944GHz -10 -20 -30 1 I. 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 -.,- \ —.. 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 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.

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. 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 Ki = jkoZoI. V(r), = jl, j2,...,jm where r is the location of the filament, m is the number of (non-metallic) element edges and jm is the global edge numbering index. In general, m such entries are associated with m element edges, and thus the index i goes from j 0 up to jm. This 98

99 expression can be further reduced to AK' = jko0Zl 1. provided thiat tIle tIli 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 1 across the opening of the coax cable or any other gap. Since V = E cl. where d is a vector whose magnitude is the gap width, and E the electric field across the gap, we have that Ei = - where cosOi is equal to 1 if the dcosO0 ith edge is parallel to d. Numerically, this gap roltage model can be realized by first setting the diagonal term Aii equal to unity and the off-diagonal terms Aij (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 all 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

1 ()0 apertulre wlhich has been extiensivel ely eplo)el inl ractice. I' li caxvity ields imay he (liscretize(l using tetrahedral elemients. whIereas in the mnicrostrip liCne regionF rectalngutlar bricks are the best candlidate since the feed structure is rectallngular in slhape and the substrate is of constant thickness. Although )both typ)es of elemients eml)loy Antenna Elments St (iir^< 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 'e(1li-pIotenltial continuity condition is eilforced. aii(l to I)roc(e(l to d(1 so. let tus first classify the slot edges as follows Tetrahedral Mesh (Cavity Region I): * Ebl j = 1.2,3..... vertical edges * Eb2 j = 1,2.3... diagonal edges Brick Mesh (Feed Region II): * EU j = 1,2,3,... vertical edges Then the 'equi —potential' continuity condition requires that Ebl = jE' (6.1) Eb = -(j E; + j+l E ) (6.2) in which +1 ej = -1 whereas t and d are the lengths of the vertical and diagonal edges, respectively. That is, t is simply the width of the narrow rectangular aperture. The coefficient Ce 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 1microstrip or cavi\tybacked patch antennas because of its simplicity and low spurious radliation. Ini(leed. abundant literature exists on the theoretical and experimental investigation of coax cable feeds [62-64]. Most of these papers present integral equation techniques ill 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 TEM mode excitation. Model improvements are also discussed for the case when the TEM assumption at the cavity-cable junction does not hold.

10:3 6.3.2 Hybrid FE-BI System The functional pertinent to the radiation I)by a cavity-backedt alt e'inna wit 1 a coax cable feed (as shown in fig. 6.3) is given by F(E) = V { xE). (V x E)-ko2c.E dIr — 2kg JJ (E x 5). { (E x ). I(+ - Vv) Go(r. r)dS dS jkoZo Jj(E x H) 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; Cr and Pr 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 INse Nce (E H {[Ae]{E}} + E {[Bi]{E}} + E) = 0 (6.4) e=l eES eEC where the explicit expressions for Aij and Bij 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 / (E x H). dS, (6.5)

1()-1 a boundary constraint relating E to H is Ilectded. To thlis eInd. w\\ asst1Ill' a 'I MNI nlode on C and consequently tile fields withiIi tlie cavity may be expIressed(l as (see. fig. 6.4) IoZo 1 o E ( + F)-. H (1 - F)-, (6.6) 27/6 r 27Ti ' where Crc 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, tile same location. Also, (raiz) 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 o = 2/ (1 + r), h = -( 1-) (6.7) Hence, o^0 ho E = r, H =, (6.8) r r and from (6.7) it follows ho = eo+-, (6.9) Zo T which is the desired constraint at the cable junction in terms of the new quantities ho and eo. 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 = -27jkoZoeohircln( ), (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 \Ne choose the linear edlge-based tetrahedral elemients to (liscretiZe tlie (avity and the corresponding mesh on the cross section C is shown in fig. (.-l(l). Ill tis formulation. the field across the pth edge, p=l2.... Nc ( '\- is the InulIl)er of cavity mesh edges on C), is set to a constant as dictated by the linear e(dge-based expansion function inside the cavity. However, the cable TENM 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)= eon-, i- = Np(p= 1,2,, Nc) (6.11) 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 do 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 Fc =- koZo(b-a) { -E E (6. 1) 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

10(; of the equal element fields on the cable's apertiure (becauIse of tlie axisvi!llletric property. all elements fields at the cable's aperture are equal). TIle factor i1Isi(e the curly brackets of (6.12). with the superscript src. functionls as a source ill the extremization process. Hence, the extremization of (6.12) yields 9Fc 1 ____o __b \ aF = — 7L jkoZo(b-a) { I / b - a } Ei 3 7 T Zo In, ' = U;E - V. i = A(p = 1.2...c), (6.13) where (b a) (6.14) (i j-7ko _V-Cr c (6-14) 3 in a k1 = jkoZo(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 equi-potential constraint (6.11). To proceed, we assemble the FEM system together with (6.13). Specifically, each Ui 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 n = H.r rd/ 21o (616) eoln (6.16) where Zc 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

1 07 certain apI)lications. To overcome tllis linlitation. on()i a)pproach is to t e'x(' t lie mesh (say. a distance d) into the cable. I he ecqui-p)otential con(litionl wtill tllecn he appliecl at z=-dc. where all higher order modes vanish. Tiis scheme requires a llor(e 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 tlihe 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 TEMI mode. With this understanding, a possible suitable expansion for the field in the coax cable (using shell elements rather than tetrahedrals) is 4 E = E q (6.17) q i where q=r, 6 or z, i=l,2,3 or 4 and N (r, /(, z) is the shape function for each of the 12 edges (3 directions x 4 edges per direction). They are given by Na = - - (qb - qb)(qc - qc)a (6.18) with q, Qb and qc representing r, b and z in cyclic rotation and correspondingly qa, qb and qC represent the parameters, < and. Also, i denotes the edge number along each coordinate, and Aqa is the width of the edge along the ^a direction. The correspondence between the edge numbers and the node pairs for each coordinate(r, T or z) is given in Table 6.1 along with the definition of the tilded parameters in (6.18). When an axisynmmetric 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

I 0S TABLE1 j'l ((rl i ttl. -! i ll ol |i p>ara i Ket -r!|. i is I.1:3 ( i, —: [ 1 (- -3 + Ao zI i (' C'(;' __) o-C, + A o _ 1 7 (1- ) + rl +r l + 0 2 (U-7). +Ar c 1 (1 1-.) -r l + A+r o1 + Ao~ '2 (5-S) + OAr ~ 3 (6-7'I + rl ~l -1 (2-3) I 'l- z!__ c +- X 3 (3-7) + rl A1 a1 (1-8 ) r7- + A~ ____ z r2 -I 2 I 02 x Table 6.1: The correspondence between the edge numbers and the node pairs for each coordinate(r, d or -) 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 (cr=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

10\9 other patch antenna whose input iIll)edance was melasured by Aberle ail(l lPozar [(i5]. This patch had a radius of 2.0 cmn and the 0.218 cm thick subl})strate liad(1 — 2.33 antl a loss tangent tan6=0.0012. The feed was located 0.7 cm fronl thle ('entert. and( for our FE-BI calculation the patch was placed in a circular cavity of 2.44 cmi in radiuis. As shown in fig. 6.7, the equi-potential model is again in excellent agreemient witlh 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 implemient 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 Y t Cavity Cavity -- Coax cable opening (surface C) x Coax cable Figure 6.3: Illustration of a cavity-backed patch antenna with a coax cable feed. cavity A patch - r4 - -9P 2b 2 2a 1 cavity-cable junction (a) (b) Figure6.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 Companson of E-Field Along (1-Wavelength) Coax Cable c o -0.0 *T 0 ID. - ULJ segment number (a) Comparison of E-Field Distribution Along R 6000, 5500 ' 5000 4500 N ' 4000 to < 3500 "o 300oo0 D 2500 2,o, 2000 1500 I 0... 11... '06.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 180 140. 100. 60.0 20.0 -20.0. 1.00 2.00 3.00 4.00 5.00 frequency (GHz) (a) 180. c CC s CE C.c:._E 140. 100. 60.0 20.0 -20.0 -60.0 L1.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; e6=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 100 E 80 C o 60 C40 n2 5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 frequency (GHz) (a) x100... 0 X 20 65 (a) 0 - E 60 a. -20._ 0 0 -60 -8 5 8.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 frequency (GHz) (b) Figure 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; c,=2.33; tan6=0.0012. [65]. -: measurement; xxx: this method; o o o: probe model [65] (a) Real part; (b) Imaginary part.

11 1 CiV'Ctit. Nc Shall leave t le topic to t he next cIhapter III colij uiict ionl withI ot her related SuIjc~t s.

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

116f volumie as required by) any PDE solution. Amiong tliein. a class of AMB(s is based on the one-way wave equation metlho(l [66.67] and anIotlher is deriived startinig wit l tlile \Vilcox Expansion [68.31]. Also. higher order ABCs using Htigdon's [69. 70]fotrulation and problem specific numerical ABCs have been successfully used. pIarticularly for truncating meshes in guided structures [71]. There are several difficulties with traditional ABCs. Among them is accuracy control. conformalitv. ease of parallelization and implementation difficulties when dealing with higher order ABC's. Also. tle 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 ABC, 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(PML) 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 equationls andl cannot be easily imlpleienltedl in finite eleieiit (FEM) solution. A new anisotropic( uniaxial) artificial absorl)ber [77] was introduIcedl recently for truncating FEM meshes. This artificial absorber is also reflectionless at all incidence angles. Basically, by making appropriate choices for the constituttive 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 sufficient 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- + Vre2Z (7.1) where V is proportional to the magnitude of E with 1V being the incoming and 1r the reflected wave amplitude, z is measured from the reference plane Sref, 7 is the

11 Microstrip Line Shielded FEM Region 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 S1. Since in (7.1) a, Vt 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(),-) + - zo V(zo) V(z+) = ZO - 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 from FEM computations, it follows V(d) = Ve-d + VreYd (7.3) V(0) = Vi + V V(-d) = ed + Vre-d (7.4) (7.5) To solve for y, we first add (7.3) to (7.5) to get (K + Vr) (eYd + -^d) = V(d) + V(-d) (7.6)

1 19 Thenr eliminatiIg 1; + \I fronm (7.4) and (7.6). woe obtain V'(d) + '(-(/) cosh (dj, = ( 2.( ) from which '/ can be determined. The effective guided wavelength A and effective dielectric constant eff may then be calculated by A - 2r eff = V A with /3 = Im{} and oA being the free space wavelength. From(7.3) and (7.4), I; and 1V are expressed as V(0)ed - V(d) 2 sinh(7-d) Vr = V(0)-Vi (7.9) Therefore, the reflection coefficient becomes Sl = r- (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-paraimeters from a possible discont illuit N regi onl alongI a transmission line mav be extracted at a distant reference plane. \lwhere thlere 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-7'2 + Re1z) z C Si7 Ey(x)= - (7.11) TEo(x)e Y Sout where Eo(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 S,,t (port 2), and 71>, 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 4 = 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 Asinh -ya z > d iEn=l,odd An cos sinh ( ) z (x,y)- = (7.13) O__~ Y sinhi (n~-d) nn n\ 1 3 n=lodd An ( ))cos ( ) sinh (b-y) z<d sinl b - d)) a a where sin (nr 1 2a 1 An - ton2fn F

121 W o 0 ~-a/2 a/2 Figure 7.2: Illustration of the cross section of a shielded microstrip line. with 2 f,-^ sin ( sni) d r1 ) s = i - sn h(i (b d))( ) csnh ( -d) at Sif to truncate the computational domain. This truncation simultaneously in+ioned 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 a function of thickness/frequency and for different loss factors.

1 2' 7.4.1 Theory Consider the waveguide. sllieldced ic'ost ri line aIld scatterer shlown\II ill ig. 7.3. Of interest is to model the wave propagation in tliese sticructulres using tlie tinite Electric Probe / II ~r = (- ). at ui - (a). waveguide -1 ~x=%y= e=a-jp n,=ly= p, a-jp d+t=40cm t= 5 cm cross-section: 4.755x2.215 cm d Absorbing layer Microstrip line = \ \ Er- r - - ~ = 3.2 d + t = 12.0 cm L = 2.38 cm h = 0.21 cm H = 1.06 cm w = 0.548 cm (b). Microstrip Line 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. = V x x E(1) - kor E EdV Sin + Sou t in which i7r and fr denote the permeability and permittivity tensors whereas E is the total electric field in the medium. The surface integrals over Sin and Sot 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

ext ract the reflect ion coefficient compuIte(d ii1 tli presenc)IsIe of tlie absorb)ing layer 1 se(I to terminate tl1e comp)utational domain. In this st(uy \\e( considerI tl Ie perforIn()aI i(ac of a thin uniaxial laver for terminating tihe FE mesIl in a rectanllgular wav(cegui1(e and a microstrip line. Such a uniaxial layer was proposed by Sacks Ot.al. [77] wlio considered the plane wave reflection from an anisotropic interface (see fig. 7.- ). If Region I Region 2 Reflected wave " _ Transmitted wave ------------ ------------- a l 0 0 a,O O0 Incident wave 0 b, 0 O b, 0 0 C, 0 c c v z Figure 7.4: Plane wave incidence on an interface between two diagonally anisotropic half —spaces. Pr and -r are the relative constitutive parameter tensors of the form a2 0 0 r =r= 0 2 0 (7.15) O 0 c2 the TE and TMl reflection coefficients at the interface (assuming free space as the background material) become cosTi - cosOt ]TE T E cosOi + &T cosOt a2 (7.16) a cosOt - cosOi RTM - a cosoi + J2 cosOt

12-1 aIld b choosiIlng = b2 anId ec = - it follows that IR' =- I? - 0 for all iIcile('1t1 angles. imnplying a perfectly miatchied material interface. It'f we set (a) = o-ji. tlhe reflected field for a metal-backed uniaxial laver is R(o0) = I-23kt0 (.17) where t is the thickness of the layer and 0, is the plane wave incidence angle. The parameter ok is simply the wavenumber in the absorber. Basically, the proposed metal backed uniaxial layer has a reflectivity of -30 dB if 3tcos0 = 0.2T7A or - 55dB if /3tcosOi = 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 uniaxial absorbing layer as a function of the parameter j3, 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 = rb(a - j3) for the permittivity tensor and a2 = rb(i - j/3) for the permeability tensor, where Orb and /rb 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 crosssection has dimensions 4.755 cm x 2.215 cm 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 strerigth ifnside the wavegui(le \whlich hias beenI termilinate(d by a p)erfectly lmatch(led uniaxial laver. As expected. the field decay inside the absorberI is eCxpoIlelltial adll for 3 values less than unity the wave does not have sufficient decay to sipl)Iress reflections from the metal backing of this 5cm laver. Consequently. a \'S\W\ of about 1.1 is observed for 3 = 0.5. However, as 3 is increased to unity. the \VS\\ t is nearly 1.0 and the wave decay is precisely given by e-3t0cos01 = v'7C0os P wh}Cere t is the wave travel distance measured from the absorber interface. P = 23t/A and here Oi = 44.5". It is noted that when 3 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 TElo mode reflection coefficient as a function of 2/3t/Ag, 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 3 changes with frequency but the broadband properties of the absorber are still maintained since acceptable low reflections can still be achieved for unoptimized 3 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 / = 3 (corresponding to 2ftt/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

1 2(; the discretization rate Iplays a role in finding tle ol)tiullll aluet' of.t/A,1 andl tlinls the presenited curves refer to a sampling rate of arouind 18/AS for tlie w\a\velli(d example. Not surprisingly (see (6.6)). for this example. the v-alue of o does not play an important role in the performance of the absorber and this is demonstrated in fig. 7.S. As seen. setting always o = 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 o will play a role in the presence of attenuating modes and it is therefore recommended to choose ac = A3 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 2Pt/Ag, where Ag = Ao//eff7 and eff 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 d3 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/3t/A9. For the waveguide and microstrip examples, we observe that the absorption is maximized for approximately the same value of 2't/Ag

1l7 (about 0.S8). Thus these cutrves can b)e used for other applications as well. alth llouglit should be note(l tlhat the discretization rate p)lays an equally inpI)ortalit ()role aIl(d this needs further investigation. The accuracy and validity of the PMIL applications for circuit p)aramieter coimputations can also be seen from the result illustrated in fig. 7.11. It is seen that use of the optimized 4.5 cm PAIL layer, with a = 1 and 3 = 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.175mm. 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].

12 S aOr I l,4 I0 0 > 200 m c 150 0 I so *o -50 L - - - Real Part - - - - Imaginary Part Magnitude, \, / ' ' I \' \\,, d -1 Ex= y= 1e = = 1-jp -1 x= y,= 4C= 1-j d+t=40cm t=lcm _1 ~ - %, I Figure 7.5: -,0 0 10 20 30 40 50 60 70 80 Samples Along Waveguide (14samplesl/) Typical field values of TE10 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 TE1o 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( [) r,TI -5 --- f=4.0GHz --- f=4.5GHz - =5.0GHz -10F X -15.3 0) > -20 (c m -25 1* ' -30 -35 -35 -40 -45 _ 0.2n -...... Figure 7.7: 0 2 4 6 8 10 12 14 2pt in kg (a=l.0) Reflection coefficient vs 2/3t/Ag (a = 1) for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6. 0 r- I I I I -5 -10h '0 -15 -20 m-25 C: -35 -40 -45 -50n 1 - Figure 7.8: 0 2 4 6 8 10 12 14 2pt in Xg (ax=P) Reflection coefficient vs 23t/A,, with a = i, for the perfectly matched uniaxial layer used to terminate the waveguide shown in fig. 7.6.

13() 0 0.5 1 1.5 2 2.5 3 3.5 4 23t in kg (a=1.0) Figure 7.9: Reflection coefficient vs 2m3t/Ag with o=1, for the shielded microstrip line terminated by the perfectly matched uniaxial layer. A 0) L.so 0 E III c.JC -5 --- f=4.0GHz -— f=4.5GHz -10 f=5.0GHz -15 -20 - -25 - -30\ -35 - ~, -40 -, -45 - - -II II50 l lt -- I Figure 7.10: 0 0.5 1 1.5 2 2.5 3 3.5 4 213t in Xg (ca=3) Reflection coefficient vs 23t/Ag with a = 3, for the shielded microstrip line terminated by the perfectly matched uniaxial layer.

131 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 width of microstrip (cm) Figure 7.11: Input impedance calculations for the PML terminated microstrip as compared to the theoretical reference data.

1:32.305.305 _.61J_1_ 1.525 _1_.61J r -T w I V A0 x (unit: mm) in CO Coj Figure 7.12: Illustration of surement. 0.9 - a meander line geometry used for comparison with mea I I I 0.8 0.7 0.6 m -0.5 I co + / 0.4 0.3F 0.2 0.1 1I v 5 10 15 20 25 frequency (GHz) 30 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 133

1:3-1 velop ani approximate transfer function of a given linear svst(In froil a1 lillited set of spectral solutions. Typically. a Pade expanlsion of thle transfer function is i)ostulateld whose coefficients are then determined 1)by matching the Iad' representation to tlhe available spectral solutions of the complete system. In this chapter we investigate the applicability of the AWNE 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 o V x E > -k2 < E, 3-E > +kb.t. (8.1)

1:3) \Nwhlere <. > deIlotes an ilnner-product and b.t. is a I)ssiblle bollnlar tIerin whose specific form is ot rei r t is ot e l fr this ci. Also. tle d(ladics Ot aiid aire riaterial related coefficients. k = - is the wavenuinber anld; is the corres)c poiind(i A c operating frequency with c being the speed of' light. A discretized formI of (8.1) incorporating the appropriate boundary conditions is [29] (Ao + kA + k2A2) {X} = {f} (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 A, 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} = k{fi} (8.3) with {fi} being independent of frequency. This observation will be specifically used in our subsequent presentation.

13:( 8.2.2 Asymptotic Waveform Evaluation To describe the basic idea of.-\\E ill coIljulnction withl the Fi'M\. \we- be('gin 1)! first expanding the solution { } il a Taylor series al)out Ao as {X} = {X-} + (k - ko) {X1} + (k - ko)2 {2,} +... +(k- ko) {X} + {(k- ko)+ } (8.4) where {Xo0} 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} = koAo {f1} {X1} Ao [{fi} - A1 {Xo} - 2koA2 {Xo}] {X2} = -Ao [A1 {X1} + A2({Xo} + 2ko {Xl})] (8.5) {X1} = -AOl [A1 {X1-l} + A2({X1-2} + 2ko {XI-l})] 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 Ao1 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 interested 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) X'(kc). The above moments can then be reduced to scalar form and I I ~vLILIIIL3L~~I~II V~LULU CV~C~LI III ~~I

137 the expansions (8.5) becoime a scalar represent at ion of \'(A ' aboutt t lie( correspoI(1it iIl' solution at ko. To yield a mlore convergent expressioll wte cain instcat rever t tlie moments to a Pacle expansion. whicil is a conventional ratioIial fllnctioll ill formll. special case of the qth order of such an expansion is given by \ ao + al(k - ko) + a2(k - ko)2 + * * + aq(k - ko)q q ) 1 + bl(k - ko) + b2(k - o)2 +... + bq(k - k.o) () where as and bi (i = 0, 1.... q) are referred to as the Pade coefficients. For transient analysis, it is observed that the Pade expansion can be reformulat ed by partial fraction decomposition [82,84] as q XP(k) = XP + k -0 (8.8) i=1 -- 0+ k - kO - k. where Xqo 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

13S of the exponential lboundary integral kernel to facilitate t lie extractiion of tlie syst emi moments. his approach does increase tlhe complications for ipIl)leInctin g AWE.\\ It however remains far more efficient in terms of CPV requiremenlts whenI compI)ared(l to the conventional approach to continuously repeating the solution of the entire svstem. 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.35cm 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 e6 = (3.2, -3.2) and r =- (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 ANWE representation given in (8.8) to approximate the reference response. As seen in fig. 8.3, the 4th order AWE

repIresenttat ion is inI agreemieInt NWith1 the real a(nd reactive I)arts of tlie reference iInpIit imIpedance solution over about 56/( and 33:(/( b>anIl(i(lt 1. res)cctiv'elv!. l lis ('learly shows that the contributions of tihe system Ipoles inI thle coI1Ilex A plaine lead( to ani accuracy difference to the real and reactive components. Surprisingly. thle Stl ordler AWE representation recovers the reference solution over the entire 1-3GHz bandl for both impedance components. \Ne also observed that the C'P requirements for 4th and 8th order computations are nearlyN 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.

1 -1 1.06 cm Figure 8.1: Illustration of the shielded microstrip stub excited with a current probe. 1.8 2 2.2 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.

1-11 e, E - I a) a. E (a) 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

I-1:3 * ITERATIVE SY'STEM SOLVER A memorv saving algorithm ITPACIK was intertwinled witli tlie lIhybrid E N subsystem to register only the non-zero FENI entries. The 13iC( iterative solver was independently developed for partially sparse and partially full matrices 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-f) 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

14-1 slotline. microstrip line. coax cable. etc. * PRISMATIC FENI ELEMENTS INCORPORATION A major problem in any hylibrid FEMI analysis is the tediolus pr)Ie-processilig for mesh generation. Thin layer substrates in the presence of thlick 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 other 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

1-15 was incorporated into thle finiite elemenIt nlethod. It was indeedl observed\( tliat the attractive features of A\E are nlainltaiIle(d wileii ulse(l in electl(romlagnettic 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 FEIM 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.

I-1 \With a whole set of specially deNveloped techniiques and Illetllodologies. one shlould then consider to create an integration env\ironment. As shown in fig. 9.1. we I) propose this FENI modular environment for future computational electromagnetic applications. A well designed modular finite element methods w-1ill be the most capable anld robust in the future! FEM Multi-Module Environment User Interface User's Modules Kernel FEM Modules FeedLine Mesh bricks I Truncation prisms Modules Modules tetras.................... --- —--------- 1 -—...... u Post Processor Figure 9.1: Multi-modular FEM environment

APPENDICES 147

1-lS 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 i i2 D 1 2 ~ 1 3 ~ 1 4 (2 3 4 2 _ 3 4 (b) Figure A.1: (a) A tetrahedron. (b) its local node/edge numbering scheme

1 -19 are expanded as where the basis functions; 6 E=ZEWe i=l I7 are given iby x r rC I' outside element ri, ri2: position vectors of vertices i1 and i2 (see Table) W>_i(r) f7-I g7-I ei bi Ve f7-i + g7-, b ~} b7-i = 6 —, ri x ri2 bib,_i 6Ve (r2 - rl ) bi = ri,2 - ri1 I = length of the ith edge (see Table) = element's volume We note that V.w = I~V V x Wx = 2gi indicating that W. are divergenceless. Furthermore, W(r) - j = ^; i - 0 iZj where r3 has its tip on the jth edge of the tetrahedron. This last property ensures that the coefficients Ee = E ei 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 A. We have 233 f (V x VW). (V x W-) g, JJJVe / r

150 Also. j ' W.W d/ = r /j/ {(f fj) +(r D)+ (g x r) (g, x r)} (d (r(I + 12 + 13) where D=(f, xgj,)+(fJ x g) and 1 = J fi fj dv '2 = r. Ddv 13 = J(g xr) (gj x r)dv Since f is a constant vector, I1 reduces to i =fi- fjVe To evaluate I2 we first set 4 4 4 x = Lxz; y= Ly; z = Lzi i=1 i=1 i=1 where x, yI, z (i = 1,...,4) denote the (x, y, z) coordinates of the tetrahedron's vertices and Li are the simplex coordinates or shape functions for the same element. That is, Li 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= i Xi + Dy Yi+ Diz zE i= =1 '=

151 whlere Dn is the rntl comllponenI t of D. I lie evaluation of 1S: caii be sipljIfi jedl 1I- t ie use of basic vector identities. \\e hlave 13 = g- gj Ir|2dc-; (gi*r)(gj r)dv = (giygjy - gigj) /j x2 dv + (gSg + gigjz) I y2 dv + (girgj +.Vgjii) /j - (ld - (gixgy + gjxgiy) xy dv - (gigjz + gjrgiz) x dt' — (giQgj -Q- +gjg) f where g(7 represents the 7mth component of the vector g.. Each of the aboNve integirals can be easily evaluated analytically and the result can be expressed in the general form I alam dv = I- ariami ali am a 20 20 1 i=1 i=1l for 1,m = 1,....,3. The parameters al or am can represent any of the rectilinear variables x, y, z.

1') 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 J J/ (-kSi. S, + V x S. V x S,) 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 Si(r) = c(r)i Z (B.2) Ap 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 k~21 f- -lkoR Bpq_ - -)(r dS' dS A r). (r r (r-)j(r) dS' 87A ii (B3 3 2 JJTJJT R

in which f? = Jr - r'. These integrals cain l)e readily evaltuated for IiOIi- self cell terms by numerical integrations. It is also observed t}hat once 7' (coi Icid(lts wit II Tj. the integrands become singular because of tile Green's function. IIn this case. the singularity should be removed. For the secondl integral, thllis is accomlp)lislhed(l by subtracting and adding an additional term. That is J]r //] r -j~oR J f -joR -_ + J - i(r)c-(r) dSS dS (B.4) 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 r r r r -jko R J.JJ -(r-r) ' (r (r)cj((r) dS'dS (r - r',)-(r - r)(r) (r) -jkoR 1S'S (+ j r - r-). (r- rj)6i(r)6j(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].

1 -1 APPENDIX C Formulation for Right Angle Prisms For FEM implementation, the following quantities are required Pmn = X VxVmVx Vd (C.1) Je Qmn = j Vm ndV (C.2) where the curls are given by V x Vi = -2S7 [(x - X) + (y -Yi) - ( - )] -i = 1 2,3 V x V =[(x - x j) + (y - yj) + (zc + -- z)n)] j=4,5,6 (.3) 2SeAz Y V x Vk 2Se [(Xk2 - Xkl)X + (Yk2 - Ykl)] 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, j,j'=4,5,6 denote the bottom triangle edges and 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 Pi = Cii [Diz + S( (C.4) Pjj, = Cjj, [Djj,/A + s(z)] (C.5) [ 3 (

1 Y3 Pkk'- 4k -~' [D1Az S'~(A)3I(.) 12 - pi 1.. 1 6 JS A (.7 Fkk= Pk Pik Pk- = Tk(Sv - xJS') + 1. 5 P-yS'] (c.1I0) 4(S'e)2 _ A)3 '3 2 2 (C. 14) Q kk' = AZSCTkkl (C.1i5) Q =Q = (Az.0 (C16 Qlk =Qkz Q-k = Qk1 = 0 (C.17) where Tkk' 1/6 for k = k; 1/ 12 for k:~k =i '~(C. 18) 4(5AZ)2 D = SX'X - (x 2 + x -)SX + x x S + SYY - (yi + y3)SY +P y-yjSe The remaining quantities in the above list of the expressions are defined as se dxchy SX = xdxdy

15(i S. = -X dIdy y S =.r'2 dr dy SYY = J y2 dxdy SXY = J xy dxdy e These integrals can be expressed in terms of the global coordinates of the three nodes (X -I, Y), (X-, Yj). (XNm Ym ). Specifically, assuming that the three nodes ij and m of a triangle are in counterclockwise rotation, we then have, 1 xi yz Se = dxdy = 1 yj 1 xm Ym sX = x dxdy s- (Xi + X + X) Jse 3 S = x2dxdy = {(i +Xj + Nm)2 + (x2 + X2 + x)} SYYx dxdy = ( + + + 2 ) 2 s {2 + (yr2 ~ yr2 ~ y) SXY = j xy dxdy = 12{(Xi + Xj + Xm) (YV + YV + I') + (XY I + xjN; + XmYm)}

157 APPENDIX D System Derivation From A Functional Referring to fig. 2.2, we begin with the functional fT(H)- + + (v x H ~1 (V x H - H * ofr H) d (D.1) to derive the system in terms of the scattered field. On inserting the field decomposition H = Hscat + Hinc, the functional becomes F(Hscat H ic) = (Hscat HS cat)d+f + a,+(HScat Hinc)lnd+(Hinc H Cat)nd + (Hscat, Hinc)Kf +(Hinc, Hscat)lQ (D.2) where the fact that HinC does not exist in Qa 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 = (sHscat inc (D.3) Also, in free space Qf, (H HS cat) = (Hscat Hinc)I Upon invoking the divergence theory, we have 2 (Hscat inc) Hscat x f V x H ) dS -- Hscat (n2 x rf * V x Hin') dS!T (D.4) (D.5)

1 5S Evidlently. for a self-adjoint operator. one readily recovers the systei (.-12) ob)t ainie via Galerkin's mlethod. It should be noted thlat )esides tle hoinldlary /t rallsition conditions. the self-adjoint property of a systenl oIerator sipllly re(qlires — 1 (ITh T (1.)6) In the case of a non-self-adjoint operator, it is generally not possiblle to recover the system given by (2.42) in the same manner. This is because the functional (D.2) ill terms of the scattered field is of the form (H) = / ] (V x HcV - kHS Hsca d Qd+Qf+Qa V ax HMct scIat + d (v C x H Hin at,d Hinc) dQ 2 d + V x H C.d. V X Hscat - k2H id H t dQ + Hscat (nix V x Hic) dS d - Hscat (n2 x f. V x Hin) dS (D.7) 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 lhomllogIeneous wavegui(de p)robl(llls"..Alta Frtquia —a. vol. 38.. pp. 31:3-317. 1969. [13] P.P. Silvester and G. Pelosi. Finith Ele mnmts for I-ar(l Elerctortagnetics Miethods and Techniques. IEEE Press. 1994. [14] A. Konrad. "Vector variational formulation of electromIagnetic fields in anisotropic media". IEEE Trans. MIicrowav'e Theory Tchl.. \ol. 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.

1 (i2 [2-1] A. Chatterjee. "Investigatilon of finite element - abc Iethod(s fo)r electIromllagnetic field simulation. Ph.D. Disse ratfioni. RIaditil on Laboratory. I ni( r.it/ of Michigan. 1994. [25] Jianming Jin. "The Finite eletrnen method in le(ctroitagn( tics " John \\iley Yk Sons. Inc. 1993. [26] D. Zheng and K.A. Michalski. "Analysis of coaxial fed antennas of arbitrary shape with thick substrates'". J. Electromagn. 'avt es 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.

[38] T.P. Budka. "microwave circuit electric field imagingi systems. ',.l). I)i.,st 1 -tation, Radiation Laboratory. I'nitlrrsity of.lichiqan. 1995. [39].J.I. Jin and.J.L. Volakis. "Biconjugate gradient fft solution for scattcring b1) planar plates". Electromagnetics. vol. 12. pp. 105-119. 1992. [40] A.F. Peterson. S.L.R-a. C.H. Chan. and R1. Mittra. ".\Numerical lipl( mnutations of The Conjugate GCradienf lethod and the CG-FFT for El(ctromagynetic Scattering". Elsevier Science Publishing Co. Inc.. 1991. [41] Yousef Saad. "Iterative lMethods for Sparse Linear Systems". P\WB Publishing Company, 1996. [42] R. Barrett,,I. Berry, T. F. Chan. J. Demmel, V. Eijkhout, R. Pozo. "Templates for the Solution of Linear Systems: Blocks for Iterative Mlethods". SIAM1. 1994. [43] C. Lanczos. "An iteration method for the solution of the eigenvalue of linear differential and integral operators'7. 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 R3". 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.

16-1 [53] A. Bossavit and..C. \erite. A mixedl FENI-BIENI mlethlod to solve 31) eddy current problems". IEEE Trans..lagnetic,. vol. 18. Ipp. 431-5. MNar 1S'2. [54] Z.S. Sacks and J.F. Lee. "A finite elenent time domailn methlod using prisll elements for microwave cavities". IEEE Trains. EIIC. Vol.:37. No. -1. 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 polygonal an poledral domains''. IEEE Trans. Antennas Propagat.. Vol. 32. pp. 276-281. Mlarch 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 femr 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 microstrip 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. \\ebb and \.. aiiellopoulos. ".Absorbing I)boulI(arv conl(litioIs for finite element solution of the vector wave equatiol..llicro'a( and Opt. 1(1 chn. Letters. Vol. 2.. No. 10. pp. 370 —372. Oct 1989. [69] R.L. Higdon. '"Absorbing boundary conditions for acoustic alnd elastic waves ill stratified media". J. Cormp. Phys.. Vol. 101. pp.:386-418. 1992. [70] J. Fang. "'ABCs applied to model wave propagation in mnicrowave integratedcircuits'. IEEE Trans. 1ITT. Vol. 42. No. 8. pp. 1506-1513. Aug. 1994. [71] J.-S. W\ang and R. Mittra. "Finite element analysis of mmic structures and electronic packages using absorbing boundary conditions". IEEE Tra ns. Micro ia r 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 Sciecle. [73] C. Rappaport and L. Bahrmasel. "An absorbing boundary condition based on anechoic absorber for em scattering computation". J. Electromagn. 1 Wavtes 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- ave Computer-Aided Engineering, Aug. 1992. [79] S. Kumiashiro, 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. "A\\E: Asyinmtotic \\a\forIIl EstilatioI". 1(cnical Report. SRC-C:II Research Ceinter for Colrpute r-.4idld L)esignr. (0rnluli(e-.Mlellon (Unircrsity. 1988. [82] E. Chiprout and M. Nakhla. Asymptotic IWalrcformL Ecaltuation and.loment Matching for Interconnect Analysis. Norwell. KIluwer Acad. Pubs.. 199-1. [83] J. L. Volakis, A Chatterjee, and J. Gong. "A class of hyb)rid finite eleneInt methods for electromagnetics: a review". J. Electrom agn. 1- a e Is 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 guidelies of the pml absorber for finite element simulations of microwave packages". IEEE/MTT Symposium Digenst, 1996.