RL - 2Z Multi-resolution methods for simulation and design of antennas by Lars S. Andersen A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1999 Doctoral Committee: Professor John L. Volakis, Chair Professor Linda P.B. Katehi Professor Noboru Kikuchi Emeritus Professor Thomas B.A. Senior RL-992 = RL-992

o Lars S. Andersen 1999 All Rights Reserved

ACKNOWLEDGEMENTS I sincerely want to thank John L. Volakis for encouraging me to come to the University of Michigan and for serving as my PhD advisor for the past four years. We most certainly have our differences but I feel that we both understand those and that we have built a relationship of tremendous mutual benefit. I thank the Danish Research Academy for the generous financial support that made my stay in Michigan possible. In this context, I also want to thank Thomas B. Thrige's Fond for providing valuable travel support and Olav Breinbjerg for being my Danish contact person. I want to thank all the past and current students and post-docs in the Rad Lab in general and in my research group in particular with whom I have interacted in the past four years. Several technical discussions have been fruitful for me professionally and numerous other conversations have benefitted me at a more personal level. I especially want to acknowledge Stephane R. Legault, Mark D. Casciato, Tayfun Ozdemir and Thomas F. Eibert, four superb gentlemen whose friendship and assistance have been no less than invaluable to me. I also want to thank the entire Rad Lab and EECS support staff for providing valuable assistance. Finally, I want to thank al the friends and family members who have provided me with significant support throughout the past four years. 11

TABLE OF CONTENTS ACKNOWLEDGEMENTS.......................... ii LIST OF TABLES................................ v LIST OF FIGURES............................... vi LIST OF APPENDICES............................ x CHAPTERS 1 INTRODUCTION............................. 1 1.1 M otivation............................ 1 1.2 Fundamental concepts...................... 5 1.3 Proposed approach.................................. 7 1.4 Organization.................................... 8 2 BACKGROUND.............................. 11 2.1 Vector wave equations................................. 11 2.2 Review of TVFEs.............................. 13 2.3 Sum m ary.............................. 20 3 TWO-DIMENSIONAL TVFES.......................... 21 3.1 Hierarchical mixed-order TVFEs for triangular elements... 21 3.2 Application to closed-domain problems................ 28 3.2.1 Formulation for closed-domain problems........... 28 3.2.2 Numerical results for closed-domain problems.. 31 3.3 Application to open-domain problems................... 33 3.3.1 Formulation for open-domain problems............... 33 3.3.2 Numerical results for open-domain problems.....38 3.4 Summary....................................... 45 4 THREE-DIMENSIONAL TVFES.......................... 46 4.1 Hierarchical mixed-order TVFEs for tetrahedral elements.. 46 4.2 Application to closed-domain problems.................... 49 iii III

4.2.1 Formulation for closed-domain problems......... 49 4.2.2 Numerical results for closed-domain problems.... 50 4.3 Application to open-domain problems............... 54 4.3.1 Formulation for open-domain problems........ 55 4.3.2 Numerical results for open-domain problems..... 61 4.4 Summary..................................... 6(5 5 ANALYSIS OF CONDITION NUMBERS.................. 66 5.1 Background.......................... 6...... 6 5.2 Formulation for closed-domain problems................. 5.3 Homogeneous application of TVFEs................ 69 5.3.1 Two-dimensional case.................. 69 5.3.2 Three-dimensional case........................ 71 5.4 Inhomogeneous application of TVFEs................. 73 5.5 Summary..................................... 77 6 ADAPTIVE TVFE REFINEMENT.................. 79 6.1 Background........................... 79 6.2 Adaptive refinement strategies...................... 82 6.3 Numerical results........................ 84 6.3.1 Square patch antenna........................ 84 6.3.2 Printed bowtie antenna.......................93 6.3.3 Discussion......................... 99 6.4 Summary............................. 104 7 ANALYSIS OF TAPERED SLOT ANTENNAS............... 105 7.1 Background.......................................... 105 7.2 Tapered slot antenna analysis.................. 108 7.3 Summary................................... 1:6 8 SUMMARY, CONCLUSIONS AND FUTURE WORK........... 127 8.1 Summary and conclusions......................... 127 8.2 Future work.......................................... 130 APPENDICES............................. 132 BIBLIOGRAPHY..................... 140 iv

LIST OF TABLES Table 3.1 Comparison of relevant parameters for the FEM results in Fig. 3.8.. 13 3.2 Comparison of relevant parameters for the FEM results in Fig. 3.10.. 43 3.3 Comparison of relevant parameters for the FEM results in Fig. 3.12.. 13 4.1 Definition of Case 1-3................................... 53 4.2 Computational effort for Case 1-6 for the antenna in Fig. 4.10-4.11.. 64 5.1 Condition numbers for the global matrices resulting from FEM analysis of the waveguide in Fig. 5.1......................... 70 5.2 Condition numbers for the global matrices resulting from FEM analysis of the cavity in Fig. 5.6........................... 73 5.3 Condition numbers for the global matrices resulting from FEM analysis of the cavity in Fig. 5.11............................ 77 6.1 Definition of Case 1-9................................... 85 6.2 Definition of Case 10-12.......................... 96 7.1 Definition of Case 1-4............................ 117 v

LIST OF FIGURES Figure 2.1 Illustration of a triangular element........................... 14 2.2 Illustration of a tetrahedral element.......................... 14 3.1 Geometry of a triangular element and illustration of the vectors n x (r - rn), n {1,2, 3}, describing the directions of the vector basis functions at the point P............................. 22 3.2 Illustration of the proposed vector basis function (CVj - (jVCi for a triangular element................................... 29 3.3 Illustration of the proposed vector basis function ((i - (j)(iV(Cj - (iV(i) for a triangular element................................ 29 3.4 Illustration of the proposed vector basis function (k((7V(j - C('V) for a triangular element............................ 29 3.5 Error of the dominant transverse propagation constant for a rectangular waveguide as a function of the number of unknowns. TE modes (left) and TM modes (right)........................ 32 3.6 Geometry of a scatterer illuminated by a TE polarized incident electromagnetic field and giving rise to a TE polarized scattered electromagnetic field.................................. 34 3.7 Square cylinder with a groove illuminated by a TE polarized plane wave. 40 3.8 Bistatic RCS of the cylinder in Fig. 3.7................ 40 3.9 Square cylinder with a groove loaded by a dielectric slab illuminated by a TE polarized plane wave.................................. 42 3.10 Bistatic RCS of the cylinder in Fig. 3.9.................. 412 3.11 Grating structure on top of a grounded dielectric illuminated by a TE polarized plane wave........................... 44 3.12 Bistatic RCS of the cylinder in Fig. 3.11................. 44 4.1 Homogeneous and isotropic rectangular cavity.............. 51 4.2 Convergence rate for expansion of the field within the homogeneous and isotropic rectangular cavity in Fig. 4.1 using mixed-order TVFEs of order 0.5 and 1.5...........................51 4.3 Inhomogeneous and isotropic rectangular cavity............. 52 vi

4.4 Eigenvalue error (E-formulation) for the inhomogeneous and isotropic rectangular cavity in Fig. 4.3.............................. 54 4.5 Eigenvalue error (H-formulation) for the inhomogeneous and isotropic rectangular cavity in Fig. 4.3.................................... 54 4.6 Magnetic field on the back PEC wall for the TM111 mode of the inhomogeneous and isotropic rectangular cavity in Fig. 4.3 with mixed-order TVFE of order 0.5 applied (Case 1)............................... 55 4.7 Magnetic field on the back PEC wall for the TM111 mode of the inhomogeneous and isotropic rectangular cavity in Fig. 4.3 with mixed-order TVFEs of order 0.5 and 1.5 applied (Case 3)............. 55 4.8 Side view of a cavity-backed patch antenna recessed in an infinite PEC ground plane................................................. 56 4.9 Top view of a cavity-backed patch antenna recessed in an infinite PEC ground plane for the case of a triangular patch and a circular cylindrical cavity............................................. 56 4.10 Side view of a square metallic patch antenna backed by a dielectricfilled rectangular cavity recessed in an infinite metallic ground plane. 62 4.11 Top view of a square metallic patch antenna backed by a dielectricfilled rectangular cavity recessed in an infinite metallic ground plane. 62 4.12 Real part of the input impedance of the antenna in Fig. 4.10-4.11 for C ase 1-6.................................. 65 4.13 Imaginary part of the input impedance of the antenna in Fig. 4.10-4.11 for Case 1-6..................................65 5.1 Illustration of a rectangular waveguide............................69 5.2 Condition numbers for the individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with unnormalized vector basis functions............................ 71 5.3 Condition numbers for the individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with normalized vector basis functions....................................... 71 5.4 Ratios of condition numbers the for individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with unnormalized vector basis functions...................................... 72 5.5 Ratios of condition numbers for the individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with normalized vector basis functions................................... 72 5.6 Illustration of a rectangular cavity........................... 73 5.7 Condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with unnormalized vector basis functions........................................ 74 5.8 Condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with normalized vector basis functions................................................. 74 vii

5.9 Ratios of condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with unnormalized vector basis functions........................... 5.10 Ratios of condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with normalized vector basis functions............................... 5.11 Illustration of an inhomogeneous rectangular cavity........... 75 76 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17 6.18 Real and imaginary part of input impedances Real and imaginary Real and imaginary Real and imaginary Real and imaginary Real and imaginary Real and imaginary Real and imaginary Real and imaginary part part part part part part part part of input of input of input of input of input of input of input of input impedances impedances impedances impedances impedances impedances impedances impedances for Case 1 for Case 2 for Case 3 for Case 4 for Case 5 for Case 6 for Case 7 for Case 8 for Case 9 (EIl (EIl (EI1 (EI2 (El2 (EI2 (EI3 (EI3 (EI3, Tetra). 87,Prism). 87,Brick). 87, Tetra). 88,Prism). 88,Brick). 88, Tetra). 89, Prism). 89, Brick). 89..... 90..... 90..... 90.... 91.... 91.... 91.... 92..... 92..... 92 Regions of' refinement Regions of refinement Regions of' refinement Regions of' refinement Regions of' refinement Regions of refinement Regions of' refinement Regions of refinement Regions of' refinement for Case 1 for Case 2 for Case 3 for Case 4 for Case 5 for Case 6 for Case 7 for Case 8 for Case 9 (EI1, (EI1, (EII1, (E I2, (E I2, (EI2, (EI3, (E13, (EI3, Tetra) at 4.15 GHz.. Prism) at 4.30 GHz.. Brick) at 4.30 GHz.. Tetra) at 4.10 GHz.. Prism) at 4.30 GHz.. Brick) at 4.30 GHz.. Tetra) at 4.10 GHz.. Prism) at 4.30 GHz.. Brick) at 4.30 GHz.. 6.19 Side view of a metallic printed bowtie antenna backed by an air- and absorber-filled rectangular cavity recessed in an infinite metallic ground plane................................... 6.20 Top view of a metallic printed bowtie antenna backed by an air- and absorber-filled rectangular cavity recessed in an infinite metallic ground plane.................................... 6.21 Real and imaginary part of the input impedance for the metallic bowtie patch antenna in Fig. 6.19-6.20...................... 6.22 Real and imaginary part of the input impedance for the metallic bowtie patch antenna in Fig. 6.19-6.20..................... 6.23 Regions of refinement for the metallic bowtie patch antenna in Fig. 6.19 -6.20 for Case 10 (left), Case 11 (middle) and Case 12 (right)... 6.24 Final number of iterations as function of frequency for different iterative solutions............................... 6.25 Final number of iterations as a function of frequency for adaptive solutions for Case 10 (EIl, Prism) with different starting guesses for the iterative solver.............................. 94 94 98 98 99 101 102 viii

6.26 Final number of iterations as a function of frequency for adaptive solutions for Case 11 (EI2, Prism) with different starting guesses for the iterative solver.............................. 102 6.27 Final number of iterations as a function of frequency for adaptive solutions for Case 12 (EI3, Prism) with different starting guesses for the iterative solver.................................... 102 6.28 Relative residual as a function of the iteration number for adaptive solutions for Case 10 (El1, Prism) with different starting guesses for the iterative solver.............................. 103 6.29 Relative residual as a function of the iteration number for adaptive solutions for Case 11 (EI2, Prism) with different starting guesses for the iterative solver....................................... 103 6.30 Relative residual as a function of the iteration number for adaptive solutions for Case 12 (EI3, Prism) with different starting guesses for the iterative solver............................. 103 7.1 Illustration of a metallic (dark gray) ETSA (left), LTSA (middle) and CWSA (right) backed by a dielectric (light gray)................ 106 7.2 Geometrical parameters for a LTSA (not drawn to scale)........... 109 7.3 E-plane patterns for the LTSA in Fig. 7.2...................11 7.4 H-plane patterns for the LTSA in Fig. 7.2................. 111 7.5 Co- and cross-polarized E-plane patterns for the LTSA in Fig. 7.2... 114 7.6 Co- and cross-polarized H-plane patterns for the LTSA in Fig. 7.2... 114 7.7 Co- and cross-polarized D-plane patterns for the LTSA in Fig. 7.2.. 114 7.8 Geometrical parameters for a LTSA (not, drawn to scale)......... 115 7.9 Real and imaginary part of the input impedance of the LTSA in Fig. 7.8 for different refinement schemes..................... 118 7.10 Number of iterations to reach a relative tolerance of l0-3 with a QMR solver for the LTSA in Fig. 7.8 for different refinement schemes.... 119 7.11 Regions of refinement for the LTSA in Fig. 7.8 for Case 1 with 3% refinem ent............................... 121 7.12 Regions of refinement for the LTSA in Fig. 7.8 for Case 2 with 3% refinement......................................... 121 7.13 Regions of refinement for the LTSA in Fig. 7.8 for Case 3 with 3% refinement.................................. 1'22 7.14 Regions of refinement for the LTSA in Fig. 7.8 for Case 4 with 3% refinem ent............................... 122 7.15 E-plane patterns for the LTSA in Fig. 7.8..................... 124 7.16 H-plane patterns for the LTSA in Fig. 7.8................ 124 7.17 E-plane patterns for the LTSA in Fig. 7.8......................125 7.18 H-plane patterns for the LTSA in Fig. 7.8.................. 125 A.1 Illustration of the variation of ul and v1 over a triangle.......... 134 ix

LIST OF APPENDICES Appendix A Explicit expressions for Wi, W2 and W3................ 133 B Expressions for vector basis functions.................. ]136 x

CHAPTER 1 INTRODUCTION In this chapter, a brief motivation for the dissertation is offered, some fundamental concepts are presented, a high-level description of the proposed approach is given and the organization of the dissertation is outlined. 1.1 Motivation Most practical electromagnetic problems dealing with microwave circuit characterization, scattering analysis or wireless communication contain regions where the field exhibits rapid variations and others where it varies slower. Rapid spatial field variations are often encountered in dense materials, at material boundaries and in the vicinity of corners, edges and small geometrical details while the opposite typically is the case away from such regions. For microwave circuits, examples are stripline, microstrip, slotline or coplanar line structures used extensively for constructing complex microwave components such as inductors, capacitors, resonators, filters and transformers of use in various system circuitry. For scattering analysis, I

examples are different airborne or ground vehicles for which accurate radar signature predictions are necessary for target recognition and observability minimization in commercial and military frameworks. For wireless communication, examples are numerous narrowband and broadband metallic printed antennas (rectangular or circular patches, bowties, spirals, log-periodics and tapered slots) possibly backed by dielectric or magnetic substrates of interest for antenna miniaturization and multifrequency communication purposes. Brute force application of traditional sub- or entire-domain method of moments (MoM) is of limited applicability for large-scale three-dimensional electromagnetic problems involving inhomogeneous media. Whether solving for volumetric polarization currents or repeatedly invoking equivalence principles to solve for equivalent and induced surface currents at material boundaries, the MoM leads to a set of linear equations described by a fully populated system matrix. As the number of unknowns increases, the central processing unit (CPU) time and memory requirements for cornputing and storing matrix elements become unrealistic. Moreover, for electromagnetic problems characterized by regions with different degrees of field variation, uniform discretization at a rate appropriate in the regions where the field varies the most typically leads to excess discretization elsewhere contributing even further to an Unmanageable computational problem. Acceleration methods like the fast multipole method [71] or the adaptive integral method (AIM) [12] can be invoked for avoiding the storage of a full matrix and for speeding up matrix-vector products within the iterative solver but these approaches add significant complexity to the formulation and implementation. The problems associated with uniform discretization through2

out the computational domain can partially be remedied by changing the sampling rate within the computational domain. However, such non-uniform meshing does not address the problems pertaining to a full system matrix and is unattractive in practice since robust software packages for non-uniform meshing are rare and expensive. Brute force application of lowest order scalar or vector finite element (FE) approaches only partially alleviates the problems related to traditional sub- or entiredomain MoM. The finite element method (FEM) as well as hybrid finite element / boundary integral (FE/BI) methods have been demonstrated to be attractive for modeling complex materials and fine geometrical details [91] and FE approaches for solving partial differential equations inherently overcome the MoM problem of a full system matrix by leading to sets of linear equations described by sparsely populated system matrices. However, the problems related to uniform versus non-uniform discretization are as serious for FE approaches as for MoM approaches. Robust and efficient modeling of electromagnetic problems characterized by regions with different degrees of field variation calls for a numerical method that is capable of accurately and efficiently modeling complex materials and geometric details, leads to a set of linear equations described by a (partly) sparse system matrix and in addition offers the possibility of locally increasing resolution. To this end, various MoM solutions utilizing wavelets have been proposed. Wavelets are classes of functions based on the elegant mathematical theory of a multi-resolution analysis (MRA) [21, 90, 99] and therefore have several attractive properties such as containment, upward and downward completeness, scaling, translation and, possibly, orthogonality or semi-orthogonality. They have been applied extensively in mathematics, 3

electromagnetics and signal processing for the past decade [21, 73, 90]. Within the context of the MoM in electromagnetics, wavelets have been applied as basis and weighting functions to directly produce a set of linear equations described by a system matrix that can be sparsified via thresholding [84]. Alternatively, one can use traditional sub-domain MoM to generate a set of linear equations described by a full system matrix and subsequently utilize a wavelet based transformation matrix and thresholding for sparsification [102]. Similar to the fast computation of a discrete Fourier transform via the fast Fourier transform [90], wavelet based MoM matrix elements can be computed in an extremely efficient manner using Mallat's fast wavelet algorithm [56, 73]. In addition, matrix elements small enough for thresholding can under certain circumstances be predicted a priori [73], lowering the computational and storage requirements even further. Moreover, wavelet based MoM_ matrices are inherently well-conditioned [73]. Despite all these attractive features, wavelet based MoM analysis is not free of serious drawbacks. The definition and implementation of wave-lets in the vicinity of discontinuities and material boundaries is tedious and unattractive [73]. Also, wavelets on arbitrarily curved three-dimensional surfaces are notoriously difficult to implement. Furthermore, the choice between a host of different wvavelets (Haar, Daubechies, Shannon, Meyer, Battle-Lemarie and Littlewood-Paley to mention just a few) is far from obvious when considering such factors as trade-off in resolution between domains, threshold limit, infinite versus compact support and orthogonality versus semi-orthogonality. In summary, disjoint regions with different degrees of field variation are a nearly inevitable part of any practical electromagnetic application. Accurate simulation of 4

electrornagnetic problems comprised of such regions is therefore of significant practical importance but traditional numerical methods often have difficulties doing so. It is the goal of this dissertation to develop a FE based multi-resolution method for accurate and efficient analysis of electromagnetic problems of that particular nature. 1.2 Fundamental concepts In order to understand the approach proposed in this dissertation, certain fundamental concepts must be familiar. These are presented in the following. For the initial application of the FEM in mechanical and civil engineering [108], the FE expansion of an unknown quantity was based on groups of scalar basis functions whose unknown expansion coefficients are associated with nodes of a finite element mesh. These groups known as node based finite elements are suitable for modeling scalar quantities but typically not so for simulating vector electromagnetic fields. When assigning vector field values to element nodes, values may need to be specified at locations where the true field is undefined (corners, edges), spurious modes can be generated and the enforcement of the boundary conditions occurring in electromagnetics can be a challenging task. These drawbacks prompted the introduction of groups of vector basis functions whose unknown expansion coefficients are associated with edges, faces and cells of a finite element mesh. These groups are known as tangential vector finite elements (TVFEs) [15] and have been shown to be free of the shortcomings of node based finite elements [96]. A TVFE is referred to as polynomial-complete to a given order n if all possible 5

polynomial variations up to and including order n are captured within the element and on the element boundary. Nedelec pointed out [61, 62] that it is not necessarily advantageous to employ polynomial-complete TVFEs. It was proven that a polynomial-complete expansion of a vector field A can be decomposed into a part representing the range space of the curl operator (V x A 7 0, A Z VO) and a part representing the null space of the curl operator (V x A = 0, A = V7). For representation of electromagnetic fields in a source free region, it suffices to employ a TVFE that is complete only in the range space of the curl operator. Since such a TVFE captures polynomial variations of order n interior to the element and polynomial variations of order n - 1 for the tangential field along edges, it is referred to as a mixed-order TVFE. In line with the syntax adopted by Webb [97], it will be referred to as complete to order n - 0.5 in the remainder of this dissertation. F'or extensive discussions of mixed-order TVFEs versus polynomnial-complete TVFEs, see [15, 61, 62, 67, 68, 85]. A TVFE is referred to as interpolatory if the vector basis functions forming the TVFE interpolate to (tangential) field values within or on the boundary of the element. For an interpolatory TVFE, the meaning of the expansion coefficient corresponding to a given vector basis function is typically easy to interpret physically. A class of TVFEs is referred to as hierarchical if the vector basis functions forming the TVFE of order n - 0.5 are a subset of the vector basis functions forming the TVFE of order n + 0.5. This desirable property allows use of different order TVFEs in different regions of the computational domain for efficient discretization of the unknown quantity - an approach we will refer to as selective field expansion. This 6

selective choice of TVFEs over the computational domain can lead to a memory and CPU time reduction as well as improved accuracy. 1.3 Proposed approach The FEM as well as hybrid FE/BI methods have been demonstrated to be attractive for modeling complex materials and fine geometrical details [91]. They implicitly allow a local increase of resolution, either by locally increasing the mesh density (h-refinement), the polynomial order of the expansion (p-refinement) or both (liprefinement). In the context of FE or FE/BI methods, a discretization scheme employing hierarchical mixed-order TVFEs permits robust combination of expansions of different orders within a computational domain. In a straightforward manner, it allows use of lowest order expansions where the field varies slowly and use of (different) higher order expansions where the field varies rapidly for an appropriate discretization of the unknown electromagnetic field. Although not formulated within the stringent framework of a MRA, such a discretization scheme is truly a multi-resolution discretization scheme in the sense that it allows different levels of field modeling within a computational domain. The concept of hierarchality is widely known but the applications hereof in electromagnetics seem few and far between. Nevertheless, FE or FE/BI discretization exploiting hierarchality appears to be a nearly ideal candidate for simulating and designing a large class of complex electromagnetic systems. In this dissertation, hierarchical mixed-order TVFEs for triangular (two-dimensional problems) and tetrahedral (three-dimensional problems) elements are developed and 7

tested and the effectiveness of selective field expansion is investigated for solution of electromagnetic scattering problems as well as for analysis of various narrowband and broadband antennas. 1.4 Organization Chapter 2 introduces background material. Vector wave equations used throughout the dissertation are presented and TVFEs for triangular and tetrahedral elements used for discretizing partial differential equations are reviewed. Chapter 3 deals with two-dimensional TVFEs. Hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for triangular elements are proposed and tested for solution of closed- as well as open-domain problems. For solution of certain classes of electromagnetic problems, field expansion using hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively is demonstrated to be a very promising approach in terms of accuracy, memory and CPU time requirements as compared to a more traditional approach. Chapter 4 extends the work in chapter 3 to three-dimensional TVFEs. Hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for tetrahedral elements are proposed and tested for solution of closed- as well as open-domain problems. Again, selective field expansion using hierarchical mixed-order TVFEs of order 0.5 and 1.5 is demonstrated to be a very promising approach for accurate and efficient solution of certain classes of electromagnetic problems. Chapter 5 discusses matrix condition numbers. The condition numbers resulting 8

fromn FEM analysis using the hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements proposed in chapter 3 and chapter 4 are contrasted to those of existing interpolatory and hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements. The proposed hierarchical mixed-order TVFEs of order 1.5 are proven to be better conditioned than existing hierarchical mixedorder TVFEs of order 1.5 and thus the analysis fosters no concerns for potential future convergence problems due to excessive matrix condition numbers. In addition, an approach for improving the condition numbers of FEM matrices resulting from selective field expansion is suggested and tested. The improvement comes at the expense of a more complicated formulation and computer code but does not alter accuracy. Chapter 6 focuses on adaptive TVFE refinement. A review of existing error estimators and indicators is given and the effectiveness of the proposed hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra is investigated when some of the reviewed error indicators are applied in the context of a very simple adaptive refinement strategy. The results are extremely promising for both narrowband and broadband antennas provided the refinement is carried out on a sub-domain by sub-domain basis as opposed to an element by element basis. Chapter 7 integrates the work in the previous chapters and presents an advanced application of the approach proposed in this dissertation. Specifically, adaptive TVFE refinement with hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra is used to analyze tapered slot antennas (TSAs) with large impedance and pattern bandwidths. The adaptive inclusion of a very small percentage of higher order TVFEs 9

is found to have a dramatic effect on the accuracy of the computed input impedances and far field patterns, thus justifying the approach proposed in this dissertation for large and complex problems. Chapter 8 closes the dissertation. Brief summaries and the most important conclusions for the individual chapters are given and several future tasks to be completed are suggested. 10

CHAPTER 2 BACKGROUND In this chapter, vector wave equations used throughout the dissertation are presented and tangential vector finite elements (TVFEs) for triangular and tetrahedral elements used for discretizing partial differential equations are reviewed. 2.1 Vector wave equations Consider a possibly inhomogeneous and anisotropic region V characterized by the permittivity tensor E and the permeability tensor /I. The electric and magnetic field intensities E and H, respectively, then fulfil Maxwell's equations [40] 1 V'xE=-jw/7.H-M (2.1) V x H= j E+J (2.'2) where J and M, respectively, denote the impressed electric and magnetic volume current densities. Elimination of H in (2.1) or E in (2.2) then leads to the vector 1Throughout this dissertation, a time factor eijt is assumed and suppressed. 11

wave equations V x (f-1 -V x E)- 2. E E= -V x (I- M)-jJ (2.3) V x (E-1 V x H)- w2~ H = V x ('-1. J)-jwM (2.4) for E and H, respectively. By applying Galerkin's method with the vector weighting function T and invoking elementary vector identities, we derive their weak forms T(( ) (i-t V xE) - w2 T ( E)) dV =-j T n x(7-1 Vx E) cdS JV T.Vx (- t - M) dV - w T.JdV (2.5) JV((V x T). ( -1 V x H) -w2 T~ (7- H)) dV =- Jf Tn x (-1 * x H) (lS + / T V x (-1 J) dV -j' T MdV (2.6) where n, denotes the outward directed unit normal vector to the boundary Fv of V. Boundary conditions at FV are incorporated in (2.5)-(2.6) through the first term on each right-hand side. These terms vanish at a perfectly electrically conducting (PEC) or perfectly magnetically conducting (PMC) surface. The tensors E and,u are constant if V is homogeneous while they are functions of the spatial variables if V is inhomogeneous. Letting 6o and,o denote the permittivity and permeability of free space, respectively, we introduce the relative permittivity and permeability tensors ~ and dr by c = c o (2.7) 12

Yu = oir UO. (2.8) For an isotropic material, E, and Plr are given by r = er I (2.9) Ir =/- I (2.]0) Twhere I denotes the 3 x 3 identity tensor while Er and,vr are the scalar relative permittivity and permeability, respectively. Specialized weak vector wave equations can be obtained by substituting (2.7)-(2.8) or (2.7)-(2.10) into the general weak vector wave equations (2.5)-(2.6) but for brevity such weak vector wave equations are not given here. 2.2 Review of TVFEs Consider a triangular element with nodes 1, 2 and 3 and a tetrahedral element with nodes 1, 2, 3 and 4, as illustrated in Fig. 2.1-2.2. The area of the triangle is denoted by A and the volume of the tetrahedron is denoted by V. Simplex coordinates at a point P of the triangle or the tetrahedron are defined in the usual manner [91], i.e. Ai for the triangle and i {1, 2, 3, 4} (2.12) -1; 13

for the tetrahedron where Ai denotes the area of the triangle formed by P and the nodes of the edge opposite to node Z and VK denotes the volume of the tetrahedron formed by P and the nodes of the triangular face opposite to node i. 3 P (l1,~2, 3) 1 2 Figure 2.1: Illustration of a triangular element. 4 3 1 P((1, ~2, 3,)4) 2 Figure 2.2: Illustration of a tetrahedral element. In his seminal papers, Nedelec determined the number of vector basis functions required for a mixed-order TVFE of order n - 0.5, n E N [61, 62]. The starting point is the number of vector basis functions required for a polynomial-complete expansion of order n. By excluding the degrees of freedom representing the null space of the curl operator, one arrives at =n(n + 2) (2.13) and -~, =n(n + 2)(n + 3) (2.) for triangular and tetrahedral elements, respectively. These equations hold for interpolatory as well as hierarchical mixed-order TVFEs. The vector basis functions can 14

be separated into three types, namely edge-based, face-based and cell-based vector basis functions. General verbal descriptions of each of these types of vector basis functions are given in the following while pictorial illustrations of specific vector basis functions appear in section 3.1. An edge-based vector basis function and its corresponding edge-based expansion coefficient can be associated with an edge in a, triangular or a tetrahedral mesh. Edge-based vector basis functions affect tangential as well as normal field components at edges and field continuity requirements across edges force expansion coefficients in adjacent elements corresponding to an interior edge to be related. Edge-based expansion coefficients are therefore referred to as global unknowns. They provide tangential field continuity across edges but allow normal field variation across edges. A face-based vector basis function and its corresponding face-based expansion coefficient can be associated with a face in a triangular or a tetrahedral mesh. Facebased vector basis functions affect normal field components at edges and tangential as well as normal field components at faces. For triangular elements, the corresponding expansion coefficients are not related across element boundaries. Face-based expansion coefficients for triangular elements are therefore referred to as local unknowns. They allow normal field variation across edges. For tetrahedral elements, field continuity requirements across faces force expansion coefficients in adjacent elements corresponding to an interior face to be related. Face-based expansion coefficients for tetrahedral elements are therefore referred to as global unknowns. They provide tangential field continuity across faces but allow normal field variation across edges and faces. 15

A cell-based vector basis function and its corresponding cell-based expansion coefficient can be associated with a tetrahedron in a tetrahedral mesh. Cell-based vector basis functions affect only normal field components at edges and faces and the corresponding expansion coefficients are not related across element boundaries. Cell-based expansion coefficients are therefore referred to as local unknowns. Thev allow normal field variation across edges and faces. Given a number N of vector basis functions W, j 1.., N, for anelement e, the expansion of the unknown electric or magnetic field intensity EU or He within the element is N E- = EVW (2.15) j=l or A' H6 = - H7Ww (2.16) j=l where E and H6 are unknown expansion coefficients. Their physical interpretation is determined by the vector basis functions WJ. In the following, vector basis functions for different mixed-order TVFEs are reviewed. The vector basis functions presented in this dissertation are not normalized. To normalize them, each edge-based vector basis function must be multiplied by the length of the edge it is associated with. Furthermore, the indices i, j and k are imnplicitly assumed to belong to the set {1, 2, 3} for vector basis functions for triangular elements and the set {1, 2, 3, 4} for vector basis functions for tetrahedral elements. For n - 1 leading to Nti = 3 and Net 6, a, mixed-order TVFE of order 0.5 is obtained for triangular and tetrahedral elements. Such a TVFE provides a constant 16

tangential / linear normal (CT/LN) field along edges and a linear field at faces and inside an element. Whitney [98] initially presented a mixed-order TVFE of order 0.5. It is characterized by the edge-based vector basis functions vj - CV, i < j. (2.17) The vector basis function iV(j - (jV( provides a constant tangential component along the edge between node i and node j, zero tangential component along all other edges and a linearly varying normal component along all edges. A detailed explanation of the characteristics of this vector basis function is given in [13]. For n = 2 leading to Nri = 8 and Net = 20, a mixed-order TVFE of order 1.5 is obtained for triangular and tetrahedral elements. Such a TVFE provides a linear tangential / quadratic normal (LT/QN) field along edges and a quadratic field at faces and inside an element. Peterson [66] (triangular elements) and Savage and Peterson [77] (tetrahedral eLements) proposed an interpolatory mixed-order TVFE of order 1.5. It is characterized by the edge-based vector basis functions 0VCj, i < (2.1 8) C(jVi, < (2.19) and the face-based vector basis functions C(CVVCj - jVi), i < j < k (2.20) 17

((kV(1 - (CVk), i < j < k. (2.21) Graglia et al. 137] proposed an interpolatory mixed-order TVFE of order 1.5. It is characterized by the edge-based vector basis functions (3(C - 1)(CVCj - jV), i < j (2.22) (3(j - 1)((iV(j - C), i < j (2.23) and the face-based vector basis functions (2.20)-(2.21). Webb and Forghani [97] proposed a hierarchical mixed-order TVFE of order 1.5. It is characterized by the edge-based vector basis functions C<VCj - jV~, i < j (2.24) C(V(9 + (V(, i < j (2.2_5) and the face-based vector basis functions (2.20)-(2.21). For n = 3 leading to N~r = 15 and N3*t = 45, a mixed-order TVFE of order 2.5 is obtained for triangular and tetrahedral elements. Such a TVFE provides a quadratic tangential / cubic normal (QT/CuN) field along edges and a cubic field at faces and inside an element. Peterson and Wilton [68] (triangular elements) and Savage and Peterson [77] (tetrahedral elements) proposed an interpolatory mixed-order TVFE of order 2.5. A correction of the vector basis functions in [68, 77] was later given by Peterson [65]. This corrected set of vector basis functions is the one presented here. It is 18

characterized by the edge-based vector basis functions i(2C- 1)V(j, i + j (2.26) Cj(Vi- Vj), i <j, (2.27) the face-based vector basis functions (k(2 k - 1)((Cv, - jVC), i < j < k (2.28) j(2(, - 1)((kV( - (CVC(), i < j < k (2.29) V(CjC), i < < k (2.30) k((cVj - vci), <j, i j k i (2.31) and (for tetrahedral elements only) the cell-based vector basis functions C(j((CV(i - iV,), i,j,k > 1, i # j I k < i, j < k. (2.32) The vector basis functions for the above mixed-order TVFEs of order 0.5, 1.5 and 2.5 were given explicitly since we shall use them in this dissertation. However, several other TVFEs for triangular and tetrahedral elements have been presented in the literature. We mention the work of Carrie and Webb [14], Cendes [15], GarciaCastillo and Salazar-Palma [27], Lee et al. [53, 54], Mur and de Hoop [60], Wang [94], Webb [95], Wu and Lee [100] and Yioultsis and Tsiboukis [103, 104, 105]. We note that a vast number of TVFEs for several other element shapes have also been presented in the literature. A comprehensive review of these is beyond the scope 19

of this dissertation. However, we stress that triangular and tetrahedral elements offer a, geometrical modeling flexibility that does not place serious restrictions on the classes of problems to which they can be applied. For this reason, it is not unreasonable to devote our attention to triangular and tetrahedral elements only. TVFEs for curved triangular and tetrahedral elements would provide even greater modeling flexibility but such TVFEs can be constructed from those proposed in this dissertation for straight triangular and tetrahedral elements via a straightforward mapping, see for instance [37]. We also note that Mur [59] throughout the past decade has strongly criticized the use of TVFEs in general. However, his viewpoints are not accepted within the computational electromagnetics community. 2.3 Summary In this chapter, vector wave equations used throughout the dissertation were presented and TVFEs -for triangular and tetrahedral elements used for discretizing partial differential equations were reviewed. 20

CHAPTER 3 TWO-DIMENSIONAL TVFES In this chapter, hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for triangular elements are proposed and tested for solution of closed- as well as opendomain problems. The work in this chapter is published in [7, 8]. 3.1 Hierarchical mixed-order TVFEs for triangular elements The proposed class of hierarchical mixed-order TVFEs is based on an expansion introduced by Popovic and Kolundzija [50, 69] for the surface current on a PEC generalized quadrilateral. For this expansion, it is demonstrated in [50, 69] that the surface current can be expanded using approximately ten unknowns per square wavelength as opposed to approximately one hundred unknowns per square wavelength for traditional sub-domain pulse basis functions. This suggests that the expansion introduced by Popovic and Kolundzija is very efficient. Below, corresponding vec 21

tor basis functions applicable for finite element (FE) expansion are constructed and hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 are presented. As a degenerate case of the generalized quadrilateral considered in [50, 69], we consider a triangular element of area A with nodes 1, 2 and 3 described by position vectors ri, r2 and r3 with respect to the origin 0 of a rectangular coordinate system, see Fig. 3.1. The edges from node 1 to node 2, node 2 to node 3 and node 3 to node 1 are referred to as edge #1, #2 and #3, respectively. Simplex coordinates (1, (2 and '3 at a point P described uniquely by a position vector r are defined in the usual manner, see section 2.2. We let n denote a unit normal vector to the surface of the triangle. ( #2 ( r3 r2 \ nx(r-r,) x(r-r3) / w rh Figure 3.1: Geometry of a triangular element and illustration of the vectors n x (r - rn), n C {1,2, 3}, describing the directions of tile vector basis functions at the point P. Popovic and Kolundzija expands the surface current J, over the triangle as [69] 3 3 Js = AJsn-= nV, (3.1L) n=l n=l where r - rn V = (3.2) 22 22

is a vector whose direction is from node n to P and Tn is a polynomial function of position that provides the amplitude variation of the vector current component Js, = nVn,. The polynomial On contains a number of unknown expansion coefficients. Its specific form is irrelevant at this point and will be given later. As in the RaoWilton-Glisson expansion [70], Js has no normal component along the two edges sharing, node n and Js has both a normal and a tangential component along the edge opposite to node n. Thus, the quantity F, = n x Jn = Only x V, = n x rW) (3..3) 2A - (3.3) with W, - x (r- r) (3.4) 2A has no tangential component along the two edges sharing node n and has both a tangential and a normal component along the edge opposite to node n. This suggests that the vector basis functions multiplying the expansion coefficients in the expansion (3.3) of F,z can be employed as vector basis functions for the edge opposite to node n when applying the finite element method (FEM). Considering all three edges, the FEM expansion of an unknown vector quantity F becomes 3 3 3 F = x J = n x Ja, S F, =,nWn (3.5) n=l n=l n=l where expressions for tI (depending upon the order of the expansion) and Wn (independent of the order of the expansion) are to be presented. Introducing normalized coordinates over the triangle and using relations from [69], 23

it is shown in Appendix A that Wl = ~2V3 - ((3V.2 (3.6) W2 = C3V(l - (1V(3 (3.7) W3 -=1V2 - C2V1i. (3.8) The polynomial.n is a function of position that in terms of a number of unknown expansion coefficients provides the amplitude variation of the vector component W,. It can be defined using normalized coordinates u, and vn over the triangle. Specifically, we choose u, = 0 at node n, un = 1 along the edge opposite to node n and Un = ~1 along the two edges sharing node n. A detailed description of the variation of un and vn is given in Appendix A. From [69], we have = 2 b +Z a (u-2- 1) v — (3.9) j=l i=3 where n, and n, are integer constants determining the order of the approximation and bJn and ai< are the expansion coefficients. Also, l = (2 + (31 ulVl - (2- (3 U2 - 43 + 4C 1 U2V2 == 43 - 41, U3 = 4 + 42 and u3z,3 = 4 - 42, as shown in Appendix A. The expansion (3.5) for F along with the expressions (3.6)-(3.8) for Wi, W2 arnd W3 and the expression (3.9) for 'n describes the proposed vector basis functions. However, a certain simplification provides a more familiar form. By regrouping tihe terms in (3.9) for fn, the expansion (3.5) for F can be cast into 3 \ntv n u F = E E c nu vnu (3.10) k=l m=l 24

where N"'vu = n(inu - 1) denotes the total number of vector basis functions per edge for the given values of nu and nr. Also, c j~ are expansion coefficients corresponding to edge #k while Wk, V are vector basis functions associated with edge #k. Wk"m is given as a function of (1, <2 and <3 times a direction vector (C1V(2 - C2V(1 for k = I, 2V(3 -- C3VC2 for k = 2 and C3VC1 - C1V(3 for k = 3). Except for normalization constants, the vector basis functions W":u are directly used for forming hierarchical mixed-order TVFEs. From (3.10), we recover for (n,r n) = (1,2) the three vector basis functions introduced by Whitney [98], see section 2.2. For larger values of nv and nu, (3.10) includes additional vector basis functions all of which maintain the same fundamental direction vectors ((iV(2 - (2V(1 for edge #1, C2V(3 - (3V(2 for edge #2 and C3V(i - (IV 3 for edge #3). Thus, the proposed higher order vector basis functions differ from the lowest order vector basis functions only in magnitude and hence in a given point of the triangle, the field is represented as a linear combination of vector basis functions having only three fundamental directions. This is one of the major differences between the proposed and traditional hierarchical TVFEs [14, 97]. For the latter, the higher order vector basis functions differ from the lowest order vector basis functions in both magnitude and direction. The field in a given point of the triangle is again represented as a linear combination of vector basis functions but in this case the number of fundamental vector directions used for representing the field grows with the order of the TVFE. An important property of (3.10) is that the vector basis functions W,u^ for k = 1, 2,3 and n = 1. *., Nn u are a subset of the vector basis functions W(n,+l)(n+u1) 25

for k =1,2,3 and m = 1,-, N,(+l)(ni+l). This shows that TVFEs based on the above presented vector basis functions are hierarchical. Based on the vector basis functions in (3.10) for different values of n, and n, along with knowledge of Nedelec spaces [61, 62] and existing interpolatory mixedorder TVFEs, hierarchierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 will now be proposed and compared to existing interpolatory mixed-order TVFEs. The method has the potential of providing hierarchical mixed-order TVFEs of even higher orders if so desired. Explicit expressions for vector basis functions are given in Appendix B. For the case (nPnu) = (1,2), we obtain from (3.10) a set of three vector basis functions forming a mixed-order TVFE of order 0.5 identical to that introduced by Whitney, see section 2.2. This result is expected since the lowest order expansion adopted by Popovic and Kolundzija is identical to the Rao-Wilton-Glisson expansion [70] whose vector basis functions reduce to the vector basis functions introduced by Whitney when converted using the procedure applied above. For the case (n,,nu) = (2,3), we obtain from (3.10) a set of twelve vector basis functions. To make a comparison to the interpolatory mixed-order TVFEI of order 1.5 presented by Peterson [66], see section 2.2, the twelve vector basis functions are reduced to eight vector basis functions forming a hierarchical mixed-order TVFE of order 1.5, see Appendix B. Peterson's interpolatory mixed-order TVFE of order 1.5 [66] has the desirable property of being complete to second order in the range space of the curl operator. This property ensures a complete second order expansion of a field with non-zero curl and guarantees eigenvalue solutions free of spurious non-zero eigenvalues. Since 26

Peterson's interpolatory mixed-order TVFE of order 1.5 and the proposed hierarchical mixed-order TVFE of order 1.5 span the same space (the existence of a linear transformation from Peterson's eight vector basis functions to the proposed eight vector basis functions is demonstrated in Appendix B), the proposed hierarchical mixed-order TVFE of order 1.5 has the same desirable property. However, the two mixed-order TVFEs are not identical as they have different properties and may not be equally efficient numerically. For both mixed-order TVFEs, six edge-based vector basis functions provide a linearly varying tangential component along edges while the remaining two face-based vector basis functions (identical for the two different mixed-order TVFEs - added to provide a complete linear representation of the curl of the field that is expanded) provide a quadratic variation of the normal component along edges. However, the linear variation of the tangential component along edges is obtained in two different ways. For Peterson's interpolatory mixed-order TVFE, the two unknowns per edge represent the magnitude of the field at edge endpoints. F'or the proposed hierarchical mixed-order TVFE, the two unknowns per edge represent the average field value along the edge and the deviation from this average value at edge endpoints. A generalization to even higher order hierarchical mixed-order TVFEs is possible. For the special case of (n1, nu) = (3,4), (3.10) gives vector basis functions that based on knowledge of Nedelec spaces [61, 62] and the interpolatory mixed-order TVFE of order 2.5 presented by Peterson and Wilton [65, 68], see section 2.2, can be used to form a hierarchical mixed-order TVFE of order 2.5, see Appendix B. The similarity between the hierarchical mixed-order TVFEs of order 1.5 and 2.5 is apparent and one 27

gets an impression of the needed generalizations to obtain even higher order TVF Es. However, use of TVFEs beyond order 2.5 does not seem to be of practical interest. The proposed vector basis function (V(j -~jV( providing a constant tangential component along an edge and a linearly varying normal component along all edges is pictorially illustrated in Fig. 3.2. The proposed vector basis function (i - j)( (- jV( i) providing a linearly varying tangential component along an edge and a quadratically varying normal component along all edges is pictorially illustrated in Fig. 3.3. The proposed vector basis function (k((iCVj - jV(i) providing zero tangential component along all edges and a quadratically varying normal component along two edges is pictorially illustrated in Fig. 3.4. 3.2 Application to closed-domain problems In this section, various TVFE options are employed for solution of closed-domain problems. The pertinent formulation is given and numerical results are presented. 3.2.1 Formulation for closed-domain problems Consider a source free homogeneous and isotropic waveguide with PEC walls. The waveguide region is denoted by V and its cross section is denoted by S with Yv and Fs being the boundaries of V and S, respectively. The waveguide is filled with a material characterized by a constant relative permittivity E, related to the complex permittivity tensor ~ via (2.7) and (2.9) and a constant relative permeability Fi, related to the complex permeability tensor 7 via (2.8) and (2.10). 28

Figure 3.2: Illustration of the proposed vector basis function VCVj - QCVjC for a triangular element.:L:::::::::::::l:l::::::::::::::::: 9iit:!i:: ii:i i;ii: i;ii i i:t Figure 3.3: Illustration of the proposed vector basis function ((i - j)((:i7q — (3Vi) for a triangular element. 40 -4*04.4 -40 -4 - Figure 3.4: Illustration of the proposed vector basis function Ck((iV(j -- (jV() for a triangular element. The weak vector wave equations (2.5)-(2.6) hold for the electric and magnetic field intensities E and H. However, a simplification is provided by reducing these to 29

weak vector wave equations for the transverse field components Et (transverse electric (TE) modes) and Ht (transverse magnetic (TM) modes) over S only. Assuming the waveguide is in the z-direction, we can let E = Ete-jz, H = Hte-jsz, V -= Vt + zt and T = Tt in the general weak vector wave equations (2.5)-(2.6) and substitute (2.7)-(2.10) for e and ju to arrive at the TE and TM weak vector wave equations s((Vt x Tt) (Vt x Et) - 72 Tt Et) dS 0 (3.11) s ((Vt x T,) (Vt x H,) - 72 Tt Ht) dS= 0 (3.12) with 7y -= /2 being the transverse propagation constant and k == w:6rLrCo( being the wave number. We now discretize S into Ns triangular elements via Ns S - SC (3.13) e=l and enforce the TE and TM weak vector wave equations (3.11)-(3.12) in each element S6. Next, we expand the transverse electric and magnetic fields EC and He in Se via (2.15)-(2.16). Choosing the vector weighting functions Tt = We, i = 1,, N, then leads to N N (V, x W) * (V, x Wp) dSE We * WS dSE`' (3.14) J-1 = N N (Vt x> W) _ (Vt x W )dSH = 7 / w dSHW C. (3.15) These element equations can be formulated in the matrix form [A']{x } = - [B] { } (3.16) 30

with the element matrices [Ae] and [Be] having the matrix elements A j- = f (Vt x W?). (V, x W;)dS (3.17) Se jB = WX w w dS (3.18) Se for both TE and TM modes. Assembly of element equations then leads to a global matrix equation system of the form [A]{x} = y2[B]{x} (3.19) for both TE and TM modes. We note that the assembly process includes the enforcement of boundary conditions along Fs. For TE modes, this leads to a condensation of the equation system making the equation system for TE modes smaller than that for TM modes. 3.2.2 Numerical results for closed-domain problems To validate the proposed hierarchical mixed-order TVFE of order 1.5 and to make a comparison with the interpolatory mixed-order TVFE of order 1.5 presented by Peterson [66], see section 2.2, we consider a homogeneous rectangular waveguide of side lengths a and a/2. The dominant transverse propagation constants for this waveguide are 7) =< w/a and T1M -= -/a for TE and TM modes, respectively [40]. The percentage of error in determining 7TE and 7TM as a function of the number of unknowns is given in Fig. 3.5 (log-log plot) for the mixed-order TVFE of order 0.5 (denoted 'Basis 1'), Peterson's interpolatory mixed-order TVFE of order 31

10~ 10~ -e- Basis 1 - Basis 1 -o- Basis 2 --- Basis 2 Bai3 - Basis 3- Basis 3 10-2 10 10) -2 1 0) 210 1) 02)10 10 -5 10 -3 101 10 10 10 1 Number of unknowns [] Number of unknowns [] Figure 3.5: Error of the dominant transverse propagation constant for a rectangular waveguide as a function of the number of unknowns. TE modes (left) and TM modes (right). 1.5 (denoted 'Basis 2') and the proposed hierarchical mixed-order TVFE of order 1.5 (denoted 'Basis 3'). In all cases, we observe the expected behavior that as the number of unknowns increases, the error decreases. Furthermore, higher order TVFEs are seen to be superior to the lowest order TVFE for a given number of unknowns. Matrices based on the proposed hierarchical mixed-order TVFE of order 1.5 were observed to be better conditioned than matrices based on Peterson's interpolatory mixed-order TVFE of order 1.5, something we discuss in more detail in chapter 5. More importantly, Basis 2 and Basis 3 give indistinguishable results. This is not surprising since the two mixed-order TVFEs are related through a linear transformation, see Appendix 1B, and consequently span the same space within each triangle. Thus, the proposed hierarchical mixed-order TVFE gives results that are identical to those obtained by using a traditional interpolatory mixed-order TVFE and in addition has the advantages of being hierarchical and (for this particular application) leading to better conditioned matrices. Similar observations were made for higher order modes. 32

3.3 Application to open-domain problems In this section, various TVFE options are employed for solution of open-domain problems. The pertinent formulation is given and numerical results are presented. 3.3.1 Formulation for open-domain problems Consider a two-dimensional scattering problem where the scatterer and all fields are independent of a spatial coordinate, say z. The scatterer can be composed of PElC and isotropic dielectric and / or magnetic materials and is situated in free space. The relative permittivity and permeability is denoted by Er and [,r respectively. The scatterer is illuminated by a TE polarized ' (the derivation for TM polarization closely parallels the one for TE polarization) incident electromagnetic field (EHZ) (subscript 't' denotes 'transverse to z', subscript 'z' denotes 'z-directed' and superscript 'i' denotes 'incident') and the scatterer then gives rise to a TE polarized scattered electromagnetic field (Es, Hs) (superscript 's' denotes 'scattered'). The configuration is illustrated in Fig. 3.6. The total electric field (transverse to z) is then Et = El + Es (3.20) and the total magnetic field (z-directed) is Hz = Hz + Hs. (3.21) Similar to the derivation in section 3.2.1, a weak vector wave equation for the total electric field Et can be derived from the general weak vector wave equation (2.5) 'In this section, polarization is with respect to the z-axis. 33

Hi E;S EE l. /H: H Et Scatterer Figure 3.6: Geometry of a scatterer illuminated by a TE polarized incident electromagnetic field and giving rise to a TE polarized scattered electromagnetic field. and (2.7)-(2.10) for E and 7 by carrying out the substitutions E -~ Et, T -+ Tt, V -~ Vt + ^Z = Vt, V -4 S = rFv - Fs, dV -~ dS and dS -~ dL and letting J - 0 and M = 0. We arrive at ((Vt x T, (V x t) V ) - k 2 Tt (r,Et)) dS fi T= Tt, V x ( E) dL (32) with ko -= VC0/oo( being the free space wave number. The computational domain S is of infinite extent at this point but will become finite through the introduction of a mesh truncation scheme, see below. Introducing the separation (3.20) of Et in the weak vector wave equation (3.22) gives - Tt' (Vt x (-Vt X E') - k ErE) dS (3.23) where the right-hand side vanishes in free space since the incident field is Maxwellian and therefore fulfils the vector wave equation in free space. We now discretize S into Ns triangular elements via Ns S Z (3.24) e-=l 34

and enforce (3.23) in each element,S. Next, we expand the scattered electric field E6t in 5' via (2.15) and choose the vector weighting functions Tt = W, i = 1,., N, to arrive at N 1 / ((Vt x W). (-Vt x We) - k W (rWJ)) d E+ w x (n V x Es) (dL =- W - (Vt x (-Vt x E)- r2 E dS. (3.25) JSe Pr Above, the boundary integral over Fse was not considered. This terms requires special attention as it is through this term boundary conditions, including the radiation condition, are incorporated. Two different mesh truncation schemes are considered, an artificial absorber (AA) and a boundary integral (BI). With an AA mesh truncation scheme, the mesh is terminated by a layer of a fictitious homogeneous and isotropic material backed by a PEC surface and placed a certain distance from the scatterer. It is important to note that the AA influences only the scattered field, not the incident field. In this case, evaluation of the boundary integral in (3.25) along all pertinent boundaries is trivial. In conclusion, each element leads to a set of element equations that can be formulated in matrix form as [A'] {EC} = {gC} (3.26) where [Ae] is a known N x N matrix, {Ee} is an unknown N x 1 vector and {g&} is a known N x 1 excitation vector. Subsequent assembly of element equations then leads to a sparse global matrix equation system of the form [A] {E} = {g} (3.27) 35

where [A] and {g} are a known matrix and vector, respectively, and {E} is a vector containing the unknown expansion coefficients EJ for the electric field as introduced by (2.15). With a BI mesh truncation scheme, the mesh is terminated by a contour where an exact integral equation incorporating the radiation condition is enforced. We denote this contour by ro. Consistent with the discretization (3.24) of S into triangular elements, we discretize Fo into No piecewise straight segments via No Fo -=. (3.28) -v=1 For an element Se that does not bound FO, the discretization of (3.25) is again straightforward and we arrive at a set of element equations of the form [A6] {Ef} = {g6} (3.29) where [Ae] is a known N x N matrix, {Ee} is an unknown N x 1 vector and {ge} is a known N x 1 excitation vector. For an element Se that bounds Fo, we are faced with the incorporation of the integral W- n x (-Vt x E) dL = —jwto f W. *. x Hz dL (3.30) with the index 7y (describing a segment of Fo) being a known function of the index e (describing a triangular element bounding Fo). To discretize (3.30), we introduce an expansion of Hs over Fo via No H - HZ Z X HZ WZ Z (3.3 1) Y=1 where Hz are unknown expansion coefficients and Wj are scalar pulse basis functions. In this case, we arrive at a set of element equations that can be formulated in 36

matrix form as [Ao] {E} + {Bo} H0 {go} (3..32) where [A]j is a known N x N matrix, {Eoe} is an unknown N x 1 vector, {Boe} is a known N x 1 vector, Hoe is the unknown expansion coefficient Hzs, for the yth segment of Fo corresponding to the eth triangular element of S and {gO} is a known N x 1 excitation vector. The expansion (3.31) of Hs introduces No additional unknowns and hence No additional equations must be constructed. These are provided via the exact integral equation [3] H(p) =- 4 j H (ko p - p'l (E(p') + E (p)) x dn' dS' 4 ( Jro +. / V x (H (kop - ) x (H(p) + H(p')) d S' (3.33) 3 o relating Es and Hs on Fo. In this equation, Co = VU//o is the intrinsic impedance of free space, H(2) is the Hankel function of zeroth order and second kind and p is the usual polar vector in a circular cylindrical coordinate system. Discretization of Es and Hs leads via testing to No equations that can be formulated in matrix form as {Ho} = [Ao] {Eo} + [Bo] {Ho} + {go} (3.31) where {Ho} is an unknown No x 1 vector, [Ao] is a known No x No matrix, {Eo} is an unknown No x 1 vector, [Bo] is a known No x No matrix and {go} is a known No x 1 excitation vector. Upon assembly of (3.29), (3.32) and (3.34), we arrive at a partly sparse / partly full global matrix equation system of the form [A] {E H} {g} (3.35,) 37

where [A] and {g} are a known matrix and vector, respectively, and {E H}' is a, vector containing the unknown expansion coefficients E. for the electric field as introduced by (2.15) and the unknown expansion coefficients H"s3 for the magnetic field as introduced by (3.31). 3.3.2 Numerical results for open-domain problems Finite element / artificial absorber (FE/AA) and finite element / boundary integral (FE/BI) computer codes were developed to evaluate the scattering of a TE or TM polarized plane wave by an arbitrary infinite cylinder composed of PEC and isotropic dielectric and / or magnetic materials. These are based on the FEM in conjunction with an AA or BI mesh truncation scheme as detailed above. The AA termination scheme is approximate but attractive for reasons of simplicity. We use a fictitious material of relative permittivity and permeability 1 - j2.7 and thickness (.25Ao (A0 denotes the free space wavelength) backed by a PEC surface and placed a distance 0.5A0 from the scatterer. The BI termination scheme is exact until discretized and coupled with a FE system and hence it is attractive for rigorous truncation of 1F'E meshes. For our test, the integral contour is situated a slight distance away from the scatterer so that piecewise constant (lowest order) expansions can be employed for discretizing and testing the BI. The resulting matrix equation system is solved using a quasi minimal residual (QMR) solver [72] with a relative tolerance of 10. In the following, we compare the scattering by various cylinders using different TVFE options and different uniform discretizations to demonstrate the merits of

the proposed hierarchical mixed-order TVFEs of order 0.5 and 1.5 when the field is selectively expanded over the computational domain. The reference results in this chapter are obtained using the method of moments (MoM) code RAM2D developed by Northrop. They are all found with a very fine discretization and thus can be considered accurate. To test the AA termination scheme, we consider a square PEC cylinder of side length A0 situated in free space. Centered on the upper side of the cylinder is a rectangular groove of length Ao/2 and height Ao/4. The groove is filled with a material characterized by the relative permittivity Or 2 - j2 and the relative permeability /ar = 2- j2. The cylinder is illuminated by a TE polarized homogeneous plane wave whose propagation vector forms a 45~ angle with all sides of the cylinder, as illustrated in Fig. 3.7. In Fig. 3.8, we give a comparison of results for the twodimensional radar cross section (RCS) or echo width J2-1) normalized to Ao0 as a function of the observation angle 4 2. The MoM result is denoted 'MoM'. For a mesh where the generic edge length is 0.15Ao, the FE/AA result using the mixed-order TVFE of order 0.5 is denoted 'FEM - 1 TVFE - Coarse mesh' and the FE/AA result using selective field expansion (with the groove and a layer surrounding the scatterer as the region in which the proposed hierarchical mixed-order TVFE of order 1.5 is employed) is denoted 'FEM - 2 TVFEs - Coarse mesh'. For a, mesh where the generic edge length is O.1o, the FE/AA result using the mixed-order TVFE of order 0.5 is denoted TFEM - 1 TVFE - Denser mesh'. The 'FEM - 1 TVFE - Coarse mesh' result is seen to compare reasonably well with 2 =_ 450 corresponds to backscatter and < = 225~ corresponds to forward scattering, see Fig. 3.7. 39

H1,o xo xo 4 2 4 P 4 4 - ----------- --- PEC r =2-j2 tr =2-j2 to Figure 3.7: Square cylinder with a groove illuminated by a TE polarized plane wave. 20 1,,, MoM FEM - 1 TVFE - Coarse mesh.-. — FEM - 1 TVFE - Denser mesh - FEM - 2 TVFEs- Coarse mesh -4 45 65 85 105 125 145 165 185 205 225 Observation angle ~ [degrees] Figure 3.8: Bistatic RCS of the cylinder in Fig. 3.7. the MoM result. However, discrepancies can be seen and this is not surprising since the mesh is relatively coarse. For the denser mesh, the 'FEM - 1 TVFE - Denser mesh' result shows a slight improvement. However, by keeping the original mesh and employing the proposed hierarchical mixed-order TVFE of order 1.5 close to the scatterer where the field can be expected to vary rapidly and accurate modeling is therefore necessary. the 'FEM - 2 TVFEs - Coarse mesh' result shows a significant 40

improvement. It matches the MoM result exactly except in regions surrounding nulls and it was obtained using less computational resources (less unknowns, less non-zero matrix entries and less matrix solution time) than the 'FEM - 1 TVFE -- Denser mesh' result. In conclusion, we observe selective field expansion to be superior to the more traditional approach of using a denser mesh and the mixed-order TVFE of order 0.5 throughout the computational domain. We now consider a slightly different cylinder geometry by introducing a slab of length Ao and height Ao/4 on top of the cylinder. As depicted in Fig. 3.9, the groove is filled with free space and the slab has the relative permittivity r = 2 - j2 and the relative permeability,r = 1. For the same illumination as before, results similar to those in Fig. 3.8 are given in Fig. 3.10 and they reinforce the conclusions from the previous case: The 'FEM - 1 TVFE - Coarse mesh' result compares reasonably well with the MoM result and the 'FEM - 2 TVFEs - Coarse mesh' result is, though found using less computational resources than the 'FEM - 1 TVFE - Denser mesh' result, significantly more accurate than the 'FEM - 1 TVFE - Denser mesh' result. Explicit parameter values quantifying the computational savings for the results in Fig. 3.8 and Fig. 3.10 are given in Tab. 3.1-3.2, respectively. In both cases, improved accuracy is obtained for less non-zero matrix entries (i.e., less memory) and less solution time. To test the BI termination scheme, we consider a rectangular PEC cylinder of width 3Ao and height 0.25Ao covered by a dielectric material of width 3Ao and height Ao whose relative permittivity is or = 2 - j0.5. On top of the dielectric is a grating structure of height 0.25Ao consisting of three PEC strips of lengths 0.75A0, 0.5A0 41

H1 ^ i t h 4 l r r P.. 4 EC 2 ^..-_-. --- --—. -. — PEC r =2-j2 Pr=1 To Figure 3.9: Square cylinder with a groove loaded by a dielectric slab illuminated by a TE polarized plane wave. 20 l l l l MoM FEM - 1 TVFE - Coarse mesh 10 --- - FEM - 1 TVFE - Denser mesh - FEM - 2 TVFEs - Coarse mesh:a 0 -10- r -20 -,i t; -30 45 65 85 105 125 145 165 185 205 225 Observation angle ( [degrees] Figure 3.10: Bistatic RCS of the cylinder in Fig. 3.9. and 0.75A0, respectively, separated by dielectric inserts of length 0.5Ao having the relative permittivity r, = 10. The structure is illustrated in Fig. 3.11. A structure of this type (but of different size and different material composition) is of practical interest for guiding electromagnetic waves and below we demonstrate how a selective field expansion can lead to accurate modeling of the fields in and near the grating structure and hereby accurate prediction of the scattered field. The structure is 42

Non-zero Matrix RMS FE/AA approach Unknowns matrix entries solution time error 1 TVFE - Coarse mesh 811 3955 9 seconds 3.8657 d] 1 TVFE - Denser mesh 3977 19664 111 seconds 2.1059 d] 2 TVFEs - Coarse mesh 1070 7292 21 seconds 1.5614 d] Table 3.1: Comparison of relevant parameters for the FEM results in Fig. 3.8. Non-zero Matrix RMS FE/AA approach Unknowns matrix entries solution time error B B - 1 TVFE - Coarse mesh 900 4398 11 seconds 2.1952 d]B 1 TVFE - Denser mesh 4469 22118 133 seconds 1.6213 d]3 2 TVFEs - Coarse mesh [ 1280 9466 ] 34 seconds | 0.6114 dB Table 3.2: Comparison of relevant parameters for the FEM results in Fig. 3.10. Non-zero Matrix RMS FE/BI approach Unknowns matrix entries solution time error 1 TVFE - Coarse mesh 1230 19662 172 seconds 5.3980 dB3 1 TVFE - Denser mesh 2421 43961 301 seconds 0.9807 d13 2 TVFEs - Coarse mesh 1716 25884 534 seconds 0.7705 d13 Table 3.3: Comparison of relevant parameters for the FEM results in Fig. 3.12. situated in free space and illuminated as the previous two cylinders. Results similar to those in Fig. 3.8 and Tab. 3.1 are given in Fig. 3.12 and Tab. 3.3. The results again reinforce the conclusions reported above, except that the matrix solution time for the 'FEM - 2 TVFEs - Coarse mesh' result is larger than that for the 'FEM - 1 TVFE - Denser mesh' result. We note that the pre-processing, matrix element evaluation and assembly time is significantly larger for the 'FEM - 1 TVFE - Denser mesh' result than the 'FEM - 2 TVFEs - Coarse mesh' result due to the larger BI system. This must be kept in mind when interpreting Tab. 3.3. We note that higher order TVFEs (either Peterson's interpolatory mixed-order TVFE of order 1.5 or the proposed hierarchical mixed-order TVFE of order 1.5) could 43

6, ILo 4ko 4 Hi 3ko 4 2o 2o o 2 2- 2 4, ko PI — PEC plrasas PEC IB1~p -BB P C I 7 I i Ira rI \ S I S r sit ilii 5i~U~i I:::~- II":: 'L:: ei;?~i!E I.~: ~ "ii d: *: IPs P 7 - k0 4 I / Irk; \ I I / ~r=10 ~r=2-j0.5 Figure 3.11: Grating structure on top of a grounded dielectric illuminated by a TE polarized plane wave.:10 -. 0 I 0 -10 -20... 45 75 105 135 165 195 225 Observation angle o [degrees] Figure 3.12: Bistatic RCS of the cylinder in Fig. 3.11. alternatively be applied throughout the computational domain. This approach was tested and the two mixed-order TVFEs gave similar and accurate results. However, the approach could not measure up with the selective one in terms of computational resources. 44

3.4 Summary In this chapter, hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for triangular elements were proposed and tested for solution of closed- as well as opendomain problems. For solution of certain classes of electromagnetic problems, field expansion using hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively was found to be a very promising approach in terms of accuracy, memory and central processing unit (CIPU) time requirements as compared to a more traditional approach. 45

CHAPTER 4 THREE-DIMENSIONAL TVFES In this chapter, hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for tetrahedral elements are proposed and tested for solution of closed- as well as opendomain problems. The work in this chapter is published in [5, 9] and expected to be published in [25]. 4.1 Hierarchical mixed-order TVFEs for tetrahedral elements Hierarchical mixed-order TVFEs for tetrahedral elements have been proposed up to and including order 1.5 by Webb and Forghani [97], see section 2.2. These were written up by inspection. The purpose of this section is to propose a set of hierarchical mixed-order TVFEs for tetrahedral elements beyond order 1.5. Specifically, hierarchical mixed-order TVFEs are presented up to and including order 2. 5 where the mixed-order TVFE of order 1.5 differs from the one presented by Webb and Forghaii 46

[97]. We derive the hierarchical mixed-order TVFEs from the interpolatory mixedorder TVFEs for tetrahedral elements proposed by Savage and Peterson [65, 77], see section 2.2, and the hierarchical mixed-order TVFEs for triangular elements proposed in section 3.1 in a fashion that makes tht he proposed set of hierarchical mixed-order TVFEs for tetrahedral elements the direct three-dimensional equivalent of the set of hierarchical mixed-order TVFEs for triangular elements proposed in section 3.1. ]-ierarchical mixed-order TVFEs for higher orders than 2.5 can be derived by modifying the TVFEs proposed by Graglia et al. [37]. Consider a tetrahedral element with nodes 1, 2, 3 and 4 for which simplex coordinates are defined in the usual manner, see section 2.2. Below, vector basis functions are formulated in terms of these coordinates. A mixed-order TVFE of order 0.5 providing CT/LN variation along edges and linear variation at faces and inside the element is characterized by NIt = 6 linearly independent vector basis functions. The three-dimensional equivalent of the twodimensional CT/LN vector basis functions presented in section 3.1 is identical to the vector basis functions presented by Whitney [98], see section 2.2. The 6 e(lge-based vector basis functions are given by (2.17). A mixed-order TVFE of order 1.5 providing LT/QN variation along edges and quadratic variation at faces and inside the element is characterized by N2ett = 90 linearly independent vector basis functions. Savage and Peterson [77] proposed the 12 edge-based vector basis functions (2.18)-(2.19) and the 8 face-based vector basis functions (2.20)-(2.21). The 20 linearly independent vector basis functions (2.18)(2.21) do not compare to the vector basis functions (2.17) in a hierarchical fashion. 47

We propose to replace the 12 edge-based vector basis functions (2.18)-(2.19) by (2.17) and (- Cj)((iV - (jVi), i < j. (.1) The 20 linearly independent vector basis functions (2.17), (4.1) and (2.20)-(2.21) form a mixed-order TVFE of order 1.5 that compares hierarchically to the proposed mixed-order TVFE of order 0.5. A mixed-order TVFE of order 2.5 providing QT/CuN variation along edges and cubic variation at faces and inside the element is characterized by Net = 45 linearly independent vector basis functions. Savage and Peterson [65, 77] proposed the 18 edge-based vector basis functions (2.26)-(2.27), the 24 face-based vector basis functions (2.28)-(2.31) and the 3 cell-based vector basis functions (2.32). The 45 linearly independent vector basis functions (2.26)-(2.32) do not compare to the vector basis functions (2.17) in a hierarchical fashion. We propose to replace the 18 edge-based vector basis functions (2.26)-(2.27) by (2.17), (4.1) and ( - j)2((.V(o - QV), < j. (4.2) Further, we propose to replace the 8 face-based vector basis functions (2.28)-(2.29) by (2.20)-(2.21). The 45 linearly independent vector basis functions (2.17), (4.1)-(4.2), (2.20)-(2.21) and (2.30)-(2.32) form a mixed-order TVFE of order 2.5 that compares hierarchically to the proposed mixed-order TVFEs of order 0.5 and 1.5. 48

4.2 Application to closed-domain problems In this section, various TVFE options are employed for solution of closed-domain problems. The pertinent formulation is given and numerical results are presented. 4.2.1 Formulation for closed-domain problems Consider a possibly inhomogeneous and anisotropic cavity V characterized by the permittivity tensor E and the permeability tensor 7J. The cavity V is Eassumed to be free of sources and the boundary IV of V is assumed to be a PEC surface. The electric and magnetic field intensities E and H, respectively, then fiulfil the weak vector wave equations (2.5)-(2.6) with the surface integrals over Fv eliminated. Introducing now the relative permittivity and permeability tensors gr and 'r by (2.7)-(2.8), we obtain ((V x T) ( -.1 V x E) -k2 T ( * E)) dV = 0 (4.3) j((V x T) (j-1' V x H) -k2 T ( r H)) dV 0 (4.4) with k =- wv/ o/o being the wave number. We now discretize V into Nv tetrahedral elements via ANv e=l and enforce (4.3)-(4.4) in each element V'. Next, we expand the electric and magnetic fields Ee and He in Ve via (2.15)-(2.16). Choosing the vector weighting functions T = We, i = ],., N, then leads to N N (V X W) ( 1 V x WJ) dlE = k2 We *j ( d (4.6) r z Jz:: 49

N N / (V XW,).( -1 *V XHWE) Jdll =k2 Wt ( r W) dl/H;. (4.7) These element equations can be formulated in the matrix form [Ae]xz} = k2[Be]{x} (4.8) for both electric field formulation (E-formulation) and magnetic field formulation ('Hformulation). Assembly of element equations then leads to a global matrix equation system of the form [A]{x} = k2[B]{x (4.9) for both E- and H-formulation. We note that the assembly process includes the enforcement of boundary conditions along Fv. For E-formulation, this leads to a condensation of the equation system making the equation system for E-formulation smaller than that; for H-formulation. 4.2.2 Numerical results for closed-domain problems Homogeneous and isotropic rectangular cavity Consider a homogeneous and isotropic rectangular cavity of normalized dimensions 1 x 0.75 x 0.5, see Fig. 4.1. The exact eigenvalues for this geometry are well-known [40]. A FEM solution (E-formulation) for the eigenvalues of the cavity is carried out for various tetrahedral meshes of different average edge length with the proposed mixed-order TVFEs of order 0.5 and 1.5 used for field expansion. The convergence rate for the two cases is illustrated in Fig. 4.2 where the average error of the first eight eigenvalues is plotted in percent as a function of the average edge 50

0.5 0. |i,........... 0.75 Figure 4.1: Homogeneous and isotropic rectangular cavity. length in the mesh (log-log plot). The approximate distribution around a straight line suggests that the average error decreases as xn for a decreasing average edge length. For the mixed-order TVFE of order 0.5, the exponent is 2.37 which is slightly larger than the expected value 2 [77]. This is due to the very low average error 0.56% for the average edge length 0.175. Similarly, for the mixed-order TVFE of order 1.5, the exponent is 4.66 which is again larger than the expected value 4 [77] and the exponent 3.86 found in [77] for a different mixed-order TVFE of order 1.5. This demonstrates that the mixed-order TVFE of order 1.5 proposed in this dissertation has slightly better convergence properties than the one in [77] for this particular geometry and for the employed meshes. 30 v v TVFE of order 0.5 10 TVFE of order 0.5 - curve fit o o TVFE of order 1.5 - - - TVFE of order 1.5 - curve fit 0 ' 1 - 1 0 / < 0.3 o o /, 0.1 0.0 o 0.1 0.2 0.3 0.4 0.5 0.6 Average edge length Figure 4.2: Convergence rate for expansion of the field within the homogeneous and isotropic rectangular cavity in Fig. 4.1 using mixed-order TVFEs of order 0.5 and 1.5. 51

Inhomogeneous and isotropic rectangular cavity Consider a rectangular cavity of normalized dimensions 1 x 0.75 x 0.5. It is composed of a rectangular cavity of normalized dimensions 0.25 x 0.75 x 0.5 filled with a homogeneous and isotropic dielectric of relative permittivity Er = 10 and a rectangular cavity of normalized dimensions 0.75 x 0.75 x 0.5 filled with a homogeneous and isotropic dielectric of relative permittivity r, = 1, see Fig. 4.3. The true eigenvalues can be obtained via accurate numerical solution of an exact analytical transcendental equation, see for instance [40]. 0.25= 0.5 Q'............... 05 "'" -10: 0.75 Figure 4.3: Inhomogeneous and isotropic rectangular cavity. We compare the error in computing the eigenvalues for the first few modes for three different cases. In Case 1, we apply the mixed-order TVFE of order 0.5 for a fine mesh consisting of 778 tetrahedral elements. In Case 2, we apply the mixed-order TVFE of order 1.5 for a coarse mesh consisting of 130 tetrahedral elements. In Case 3, we apply the mixed-order TVFEs of order 0.5 and 1.5 selectively for the same coarse mesh consisting of 130 tetrahedral elements. The dielectric material is modeled with the mixed-order TVFE of order 1.5. In free space, the region away from the dielectric is modeled with the mixed-order TVFE of order 0.5. This makes the region close to the dielectric a transition region where incomplete mixed-order TVFEs (complete to order 0.5 but not to order 1.5) are applied. The three cases are summarized in 52

Tab. 4.1 where also the number of unknowns and number of non-zero matrix entries are given for E- and H-formulation. E-formulation H-formulation TVFE Matrix Matrix Case order(s) Elements Unknowns entries Unknowns entries I 1 0.5 778 695 9161 1151 16051 2 1.5 130 604 17876 1084 - 38080 i i i..... I 3 0.5/1.5 130 354 9120 668 20060 Table 4.1: Definition of Case 1-3. The eigenvalue error of the first six modes for the three different cases is given in percent in Fig. 4.4-4.5 for E- and H-formulation, respectively. Comparing Case 2 and Case 3, we see that the average error is approximately the same. This is the case for both E- and H-formulation. Thus, we do not compromise accuracy by modeling only part of the cavity with the mixed-order TVFE of order 1.5 and the remainder of the cavity with the mixed-order TVFE of order 0.5 as compared to using the mixed-order TVFE of order 1.5 throughout the cavity. However, the storage and CPU time requirements drop significantly making selective field expansicn an attractive option. Comparing Case 1 and Case 3, we see that Case 1 is best for E-formulation while Case 3 is best for H-formulation. This demonstrates that the choice of TVFE(s) is not necessarily unambiguous. This unambiguity can be further demonstrated by considering the eigenmode corresponding to a given eigenvalue. As noted above, Case 3 is generally better than Case 1 for H-formulation. However, for the TM111 mode, Case 1 gives a more accurate eigenvalue than Case 3. The tangential magnetic field on the back PEC wall for the TM111 mode is plotted for Case 1 and Case 3 in Fig. 4.6-4.7, 53

1 1 Case 1 /... Case 2 2 Case 3 TE110 TM111 TE101 TE111 TM121 TE120 Figure 4.4: Eigenvalue error (E-formulation) for the inhomogeneous and isotropic rectangular cavity in Fig. 4.3. 3 Case 1 L, ~~.Case 2 i 2 Case 3 0a A r1'~ L: H.. 41 Il TE110 TM111 TE101 TE111 TM121 TE120 Figure 4.5: Eigenvalue error (H-formulation) for the inhomogeneous and isotropic rectangular cavity in Fig. 4.3. respectively. Clearly, transitions are much smoother for Case 3 than for Case 1 and significantly more accurate fields at edges can be observed for Case 3 as compared to Case 1 (the tangential magnetic field for Case 1 possesses normal components at edges that do not exist for Case 3). Thus, for a mode where Case 1 gives a more accurate eigenvalue than Case 3, Case 3 gives a more accurate eigenmode than Case 1. Further discussion of TVFE ambiguity is given in [75]. 4.3 Application to open-domain problems In this section, various TVFE options are employed for solution of open-domain problems. The pertinent formulation is given and numerical results are presented. 54

Figure 4.6: Magnetic field on the back PEC wall for the TM11l mode of the inhomogeneous and isotropic rectangular cavity in Fig. 4.3 with mixed-order TVFE of order 0.5 applied (Case 1). -R.-, ~;s-:. i:- -- S-. Figure 4.7: Magnetic field on the back PEC wall for the TM,11 mode of the inhomogeneous and isotropic rectangular cavity in Fig. 4.3 with mixed-order TVFEs of order 0.5 and 1.5 applied (Case 3). 4.3.1 Formulation for open-domain problems Consider a PEC patch antenna backed by a PEC cavity and recessed in an infinite PEC ground plane. The cavity-backed patch antenna is situated in free space characterized by the permittivity eo and the permeability do as illustrated in Fig. 4.8 (side view) and Fig. 4.9 (top view) for the case of a triangular patch backed by a finite circular cylindrical cavity. The volume of the possibly inhomogeneous and anisotropic cavity is denoted by V and characterized by the permittivity tensor E and the permeability tensor 1u. The region V can include internal PEC surfaces. Internal resistive and impedance surfaces can easily be incorporated [92] but we restrict ourselves to internal PEC surfaces in this dissertation. The boundary of the cavity V defined as 55

O toin PEC A Figure 4.8: Side view of a cavity-backed patch antenna recessed in an infinite PEC ground plane. Figure 4.9: Top view of a cavity-backed patch antenna recessed in an infinite PEC ground plane for the case of a triangular patch and a circular cylindrical cavity. the top surface (the PEC patch or the non-PEC part between the PEC ground plane and the PEC patch), the PEC side surfaces, the PEC bottom surface and any internal PEC surfaces is denoted by Fv. The non-PEC part of Fv is denoted by S while the PEC part of Fv is denoted by Fv S. Assuming the antenna feed is described by an electric volume current density J within 1, the electric and magnetic field intensities E and H, respectively, fulfil Maxwell's equations (2.1)-(2.2) for M = 0 and E fulfils the weak form (2.5) of the vector wave equation (2.3) for M = 0. Combination of the weak vector wave equation (2.5) with Maxwell's equation (2.1) for M = 0 and introduction of the rel 56

ative permnittivity and permeability tensors ar and r via (2.7)-(2.8) then leads to E fulfilling ((Vx T). (-1 V xE) - T (r E)) dV = jWo I T (i x H)dS -ji ofT TJdf (4.10) Jv where ko = wvtoo1 is the free space wave number, where n denotes the unit normal vector to S directed out of V and where a surface integral over the PEC part Fv\5,S of Fv has been eliminated since it vanishes. We now proceed to discretize (4.10) using the hybrid FE/BI method. We initially neglect the surface integral in (4.10) and discretize j ((VxT r V x T) ( r-. V xE T E)- T)dV- ( j E) JdV (4.11) to obtain a FE system of linear equations. Upon discretization of the surface integral jio / T (n x H) dS (4.1) in (4.10) using a BI (an expression for H in terms of E) and subsequent correction of the FE system of equations according to (4.10), the final FE/BI system of linear equations is obtained. We consider (4.11). Let us discretize V into Nv tetrahedral elements via Nv/ V VCe (4.13) eeeee eeeee e and enforce (4.11) in each element V6". Next, we expand the electric field E" in Ve via 57

(2.15) and choose the vector weighting functions T = We, i = 1,, N, to arrive at N 3 Z / ((V x W) * ( V X Wj) j=1 - k We (. W)) dllE = - -jio ' We J dV. (4.14) e The N element equations (4.14) for the eth element can be formulated in the matrix form [A6]{E6} = {g} (4.15) where [A6] is the element matrix with entries A7 j/ ((V x we) (,-1 V x We) - ko2 W. ( r W;)) (4.16) 33v<3 {ge} is the excitation vector with entries g = -jwo i We J dV (4.17) and {Ee} is a vector containing the unknown expansion coefficients EJ. Assembly of the element equations (4.15) and elimination of unknowns along Fv\S then leads to a global matrix equation system that can be formulated Al, AIS E g I (4.18) ASI AsS As, As Es gS where [A,], [AIS], [AsI] and [Ass] are known matrices, {g1} and {gS} are known vectors and {EI} and {ES} are vectors containing the unknown expansion coefficients associated with the interior of V (superscript I) and the boundary S of V (superscript S). This is the desired FE system of linear equations resulting from discretization of (4.11). 58

We now consider the surface integral (4.12) which is in terms of the magnetic field H. As mentioned earlier, we seek an expression in terms of the electric field E. This is obtained by using the surface equivalence principle [45] to relate H to E just outside S via H = -2jwso Go(r, r'). (E' x h') dS' (4.19) where Go is the dyadic free space Green's function given by Go(r, r') = (I + -VV ) Go(r, r') (4.20) 0 with I being the identity tensor and Go being the scalar free space Green's function given by - jko r-r'I Go(r, r) 47ir- r' (4.21) Substitution of (4.19) for H into the surface integral (4.12) yields jWbo jT (h x H) dS = 2 k 2 T (n x j Go(r, r') (E' x n') dS' ) dS (4.22) which is in terms of E and can be directly discretized. Before doing so, it is advantageous to substitute Go from (4.20) and use elementary vector identities to obtain jwuo J T. (n x H) dS =2 ko j j Go(r, r') (E' x nh') (T x n) dS' dS -2 is Go(r, r') V'. (E' x n') V -(T x n) dS' dS. (4.23) The singulatities of the integrands on the right-hand side of (4.23) are of a lower order than the singulatity of the integrand on the right-hand side of (4.22.) and hence the tangential electric field in (4.23) can be expanded using an expansion of a lower 59

order than that in (4.22). Consistent with the discretization of V into Nvi elements via (4.13), we now discretize S into Ns triangular faces via Ns S S. (4.24) Next, we expand the tangential electric field Etg linearly in Se' via Egn S = E EVJ (4.25) j=l where Etj are unknown expansion coefficients and V' =- W ' - (W ~' n) n) (4.26) are known vector basis functions, respectively, for the e'th segment. With the field representation (4.25) over S, the field uniqueness requirements at S mandate that all higher order edge- and face-based vector basis functions associated with edges and faces of S are eliminated. This implies that mixed-order TVFEs for elements with a triangular face bounding S can be complete to order 0.5 or incomplete to order 1.5 but not complete to order 1.5. Further discussion is given in section 8.2. Choosing the vector weighting functions T = Ve, i = 1,.., M, we obtain the equations Ns M j W/oJ T (n l H) dS = j / / Go(r, r') (2 ko VJ V. -2 V' (VJ' x n') V. (V x n)) dS' dS Ej. (4.27) Upon assembly of these equations, we end up with 0 0 E' jWto /f T (h x H) dS' (4.28) 0s Bss Es 60

which is the desired discretization of (4.12). The FE/BI discretization of (4.10) now follows from the FE discretization (4.18) of (4.11) representing the FE part of (4.10) by combining it with the BI discretization (4.28) of (4.12) representing the BI part of (4.10). We obtain the final FE/BI system All AIS El g: A {. (4.29) AS ASS - Bss Es gS We note that the adaptive integral method (AIM) [11, 12] can be invoked for avoiding the storage of a full BI matrix and for speeding up matrix-vector products within an iterative solver. For conservative choices of AIM parameters (dense AIM grid and large AIM near-zone threshold), these computational improvements come with virtually no compromise in accuracy. The AIM was conservatively taken advantage of for some of the computations in this dissertation but since the AIM is not essential to this dissertation the formulation is not repeated here. 4.3.2 Numerical results for open-domain problems Consider a square metallic patch antenna backed by a rectangular cavity recessed in an infinite metallic ground plane, as illustrated in Fig. 4.10 (side view) and Fig. 4.11 (top view). The cavity-backed patch antenna is situated in free space characterized by the permittivity Co and the permeability,/o. The cavity is of dimensions 1.85 cm x 1.85 cm x 0.15 cm and filled with a dielectric material of permittivity 10 Eo and conductivity 0.0003 S/cm. The patch is of side length 0.925 cm and centered in the cavity aperture. It is fed by a vertical coaxial line whose outer conductor is attached 61

Co P0 0.4625 cm 0.925 cm 0.4625 cm 0.15 cm 10~% 0 0.0003 S/cm Figure 4.10: Side view of a square metallic patch antenna backed by a dielectric-filled rectangular cavity recessed in an infinite metallic ground plane. Figure 4.11: Top view of a square metallic patch antenna backed by a dielectric-filled rectangular cavity recessed in an infinite metallic ground plane. to the ground plane and whose inner conductor is attached to the patch at the mid point of an edge, as illustrated in Fig. 4.10-4.11. The coaxial feed is modeled as a vertical probe of constant current. An almost identical antenna was considered by Schuster and Luebbers [81]. In [8-1], the cavity walls and the ground plane was removed and a similar patch on a similar but finite grounded dielectric substrate was analyzed using the finite difference time domain (FDTD) method. In spite of these geometrical differences, the two antennas are expected to have the same input impedance and, consequently, the same resonant frequency since the dominant fields are confined to a volume under and in the vicinity of the patch. The resonant frequency was found in [81] to be 4.43 GHz. The resistance 62

at resonance was found to be 400 Q while the reactance was in the range of 230 Q to -170 Q close to resonance. We note that the results in [81] were found with an extremely fine discretization and hence can be considered accurate. The patch antenna, is analyzed using the FE/BI method in conjunction with an iterative QMR solver. We discretize the cavity into tetrahedral elements and consequently discretize the surface forming the boundary between the cavity and free space into triangular faces. Two different TVFE options are applied. The first TVFE option is to use the mixed-order TVFE of order 0.5 throughout the mesh. For a mesh of average edge length 0.260 cm (Case 1), the input impedance is determined as a function of frequency and the resonant frequency of the patch is predicted. The coarse discretization of Case 1 means that this resonant frequency is most likely not accurate. For meshes of average edge lengths of 0.188 cm (Case 2), 0.153 cm (Case 3) and 0.133 cm (Case 4), more accurate resonant frequencies but also higher conmputational costs can be expected. The second TVFE option is to use the hierarchical mixed-order TVFE of order 1.5 close to the radiating edges and the mixed-order TVFE of order 0.5 elsewhere. For the meshes of average edge length 0.260 cm (Case 5) and 0.188 cm (Case 6), the input impedance is again determined and the resonant frequency is again predicted. The effectiveness of this approach (Case 5-6) in terms of accuracy, CPU time and memory requirements is compared to the previous one (Case 1-4). The six cases are summarized in Tab. 4.2. Real and imaginary parts of the input impedance as a function of frequency are given in Fig. 4.12-4.13 for Case 1-6 and corresponding resonant frequencies are provided in Tab. 4.2. For Case 1-4, a larger and larger resonant frequency is observed 63

Average Time per edge Resonant frequency TVFE length frequency BI Matrix point Case order(s) [cm] [GHz] Unknowns unknowns entries [sec] 1 0.5 0.260 3.974 345 120 17119 7.52 2 0.5 0.188 4.147 817 288 89695, 44.78 3 0.5 0.153 4.258 1489 528 291359 222.92 4 0.5 0.133 4.302 2361 840 725791 771.59 5 0.5/1.5 0.260 4.323 827 120 30675 17.33 6 0.5/1.5 ~ 0.188 4.437 1467 288 107963 77.28 -6 - Table 4.2: Computational effort for Case 1-6 for the antenna in Fig. 4.10-4.11. as the mesh becomes denser and denser. However, even for Case 4, the error as compared to the result obtained by Schuster and Luebbers is quite large (2.98 X%) for resonant frequency computation. Use of selective field expansion (Case 5-6) leads to a, significant accuracy improvement. Case 5 (error 2.42 %) gives a more accurate result than Case 1-4 and Case 6 (error 0.16 %) matches the result by Schuster and Luebbers almost exactly. The computational cost (number of unknowns, number of BI unknowns, number of non-zero matrix entries (memory usage) and CPU time per frequency point) to obtain these results are also given in Tab. 4.2. It is evident that the second TVFE option corresponding to Case 5-6 is significantly more attractive than the first TVFE option corresponding to Case 1-4. Case 5 gives a nmore accurate result than Case 4 but uses only 4.22 % of the memory and 2.15 % of the CPU time that Case 4 does. The accuracy of Case 6 is vastly superior to that of Case 4 and yet Case 6 uses only 14.88 V% of the memory and 10.02 %o of the CPU time that Case 4 does. We note that the savings in Case 5-6 are reached in part because coarlse meshes with higher order TVFEs lead to significantly smaller BI portions of' the resulting matrix equation systems than fine meshes with lowest order TVFEs. 64

3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 Frequency [GHz] Figure 4.12: Real part of the input impedance of the antenna in Fig. 4.10-4.11 for Case 1-6. 300 -C ase 1-6.Case 5 the results are not included in this dissertation. 200 Case Case 2 promising approach for accurate and efficient -20solution7 38 39 4 41 42 43 44 4'5 4'.6 4.7 Frequency [GHz] Figure 4.13: Imaginary part of the input impedance of the antenna in Fig. 4.10-4.11 for Case 1-6. the results are not included in this dissertation. 4.4 Summary In this chapter, hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 -for solution of certain classes of electromagnetic problems. 65

CHAPTER 5 ANALYSIS OF CONDITION NUMBERS In this chapter, the condition numbers resulting from FEM analysis using the hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements proposed in chapter 3 and chapter 4 are contrasted to those of existing interpolatory and hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements. In addition, an approach for improving the condition numbers of FEM matrices resulting from selective field expansion is suggested and tested. The work in this chapter is published in [6]. 5.1 Background Linear equation systems associated with the FEM and hybrid versions hereof are often solved using iterative solvers which are known to perform relatively poorly (require many iterations) when the system matrices are badly conditioned. This is the case for (generalized) eigenvalue problems as well as excitation problems. This drawback makes it desirable to construct system matrices having small condition 66

numbers. It is generally postulated that higher order TVFEs lead to larger condition numbers than the lowest order TVXFE and that hierarchical higher order TVFEs lead to larger condition numbers than interpolatory higher order TVFEs [75]. These are two postulates that make selective field expansion using hierarchical higher order TVFEs significantly less attractive. The condition numbers of element matrices based on the mixed-order TVFEs of order 1.5 developed by Graglia et al. [37] (interpolatory) and Webb and Forghani [97] (hierarchical) have been studied before [75]. However, no detailed study has been presented that examines the inter-relationships between the condition numbers of element as well as global matrices based on various interpolatory and hierarchical mixed-order TVFEs. It is the aim of this chapter to carry out such a comparison. Specifically, we compare the interpolatory mixed-order TVFEs of order 1.5 developed by Peterson [66] (two-dimensional problems), Savage and Peterson [77] (three-dimensional problems), Graglia et al. [37] and Webb and Forghani [97], see section 2.2, with the hierarchical mixed-order TVFEs of order 1.5 proposed in this dissertation. The study includes results for two- as well as three-dimensional problems, TE field formulation (TE-formulation) as well as TM field formulation (TM-formulation) (two-dimensional problems), E- as well as H-formulation (three-dimensional problems), element as well as global matrices and unnormalized as well as normalized vector basis functions. Based on the study, an approach for improving the condition numbers of system matrices resulting from selective field expansion is suggested and tested. 67

5.2 Formulation for closed-domain problems Given a mixed-order TVFE, different applications lead to structurally different global matrix systems and hence it is impossible to uniquely define a global matrix whose condition nurber characterizes the mixed-order TVFE for all applications. If we consider a waveguide or a cavity with metallic walls, FEM analysis based on a given mixed-order TVFE leads to element matrix equation systems of the form (for specifics, see chapter 3 and chapter 4) [A6]{x6} = A[B6]{x} (5.1) where each entry of [Ae] is the integration of the dot product of the curl of two vector basis functions over the element and each entry of [B6] is the integration of the dot product of two vector basis functions over the element. Assembly of element equations then leads to a global matrix equation system of the form [A]{x}= X[B]{x}. (5.2) The matrices [A6] and [A] are singular while [B6] and [B] are non-singular. In this chapter, we perform FEM analysis of waveguides and cavities with metallic walls based on different mixed-order TVFEs and use the condition numbers of [Be] and [B] as indicators of the matrix condition numbers that the different mixed-order TVFEs lead to. We note that the condition number of a matrix can be defined in a variety of ways [36]. The condition number used in this dissertation is the absolute value of the ratio of the maximum to the minimum eigenvalue. 68

5.3 Homogeneous application of TVFEs In this section, condition numbers for FEM matrices are computed when mixedorder TVFEs of order 0.5 and 1.5 are used throughout computational domains. 5.3.1 Two-dimensional case Consider a rectangular waveguide of normalized dimensions 1 x 0.5, as illustrated in Fig. 5.1. For TE- and TM-formulation, the eigenvalues are determined using the FEM (64 triangular elements) with the mixed-order TVFEs by Whitney, Peterson, Graglia et a/., Webb and Forghani and that proposed in this dissertation used for field expansion (unnormalized as well as normalized vector basis functions). 0.5, Figure 5.1: Illustration of a rectangular waveguide. In Tab. 5.1, the condition number of the global matrix [B] is given for TE- and TM-formulation with unnormalized as well as normalized vector basis functions. The abbreviations 'Wh' 'Pe' 'Gr' 'An' and 'We' denote the mixed-order TVFEs developed byr Whitney, Peterson, Graglia et a/., that proposed in this dissertation and that developed by Webb and Forghani, respectively. The higher order TVFEs are seen to lead to much larger condition numbers than the lowest order TVFE. Also, unnormalized vector basis functions are seen to lead to much larger condition numbers than normalized vector basis functions for all TVFEs. For interpolatory TVFEs, 69

the TVNFE by Graglia et al. leads to better conditioned matrices than the TVFE by Peterson while for hierarchical TVFEs, the TVFE proposed in this dissertation leads to better conditioned matrices than the TNFE by Webb and Forghani. Tlhe hierarchical TVFE proposed in this dissertation leads to better conditioned matrices than the interpolatory TVFE by Peterson, especially when the vector basis functions are normalized. The interpolatory TVFE by Graglia et al. leads to slightly better conditioned matrices than the hierarchical TVFE proposed in this dissertation. TE-formulation TM-formulation TVFE Unnormalized Normalized Unnormalized Norma]lized Wh 3 3 5 6 Pe 428 90 479 149 Gr 144 37 144 61 An 402 43 451 72 We 827 71 926 101 Table 5.1: Condition numbers for the global matrices resulting from FEM analysis of the waveguide in Fig. 5.1. To possibly correlate the global matrix condition numbers to those of the individual element matrices, condition numbers for the element matrix [Be] for each of the 64 elements used to discretize the waveguide are displayed in Fig. 5.2-5.3 for unnormalized andl normalized vector basis functions, respectively. For higher order TVFEs, the results are presented in a slightly different form in Fig. 5.4-5.5 where for each element the condition numbers of [Be] obtained with the mixed-order TVFEs by Peterson, Webb and Forghani and those proposed in this dissertation are given relative to the one obtained with the mixed-order TVFE by Graglia et al. Fig. 5.2-5.5 show that all conclusions for global matrices also hold for the individual element matrices. 70

TE/TM Unnormalized TVFE of order 0.5 - Whitney — * TVFE of order 1.5 - Peterson..........>.......... TVFE of order 1.5 - Graglia - TVFE of order 1.5 - Andersen --- TVFE of order 1.5 - Webb c ==..ro h -1g IO 0 10 20 30 40 50 60 Element number Figure 5.2: Condition numbers for the individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with unnormalized vector basis functions. TE/TM Normalized TVFE of order 0.5 - Whitney * TVFE of order 1.5 - Peterson....TVFE of order 1.5 - Graglia... v....-. TVFE of order 1.5 - Andersen 4 ---- TVFE of order 1.5 - Webb c.2.a a a 10 20 30 40 50 60 Element number Figure 5.3: Condition numbers for the individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with normalized vector basis functions. 5.3.2 Three-dimensional case Consider a rectangular cavity of normalized dimensions I x 0.75 x 0.5, as illustrated in Fig. 5.6. For E- and H-formulation, the eigenvalues are determined using the FEM (130 tetrahedral elements) with the mixed-order TVFEs by Whitney, Savage and 71

TE/TM Unnormalized E 10 0.5 0 10 20 30 40 50 60 Element number Figure 5.4: Ratios of condition numbers the for individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with unnormalized vector basis functions. TE/TM Normalized 4 * C(Peterson)/C(Graglia)........... -. C(Andersen)/C(Graglia) -3 —.5 ---- C(Webb)/C(Graglia) C 3 " 2.5 - -2.5 0.5 0 10 20 30 40 50 60 Element number Figure 5.5: Ratios of condition numbers for the individual element matrices for the waveguide illustrated in Fig. 5.1; TE/TM-formulation with normalized vector basis functions. Peterson, Graglia et al., Webb and Forghani and that proposed in this dissertation used for field expansion (unnormalized as well as normalized vector basis functions). Condition numbers similar to those in Tab. 5.1 and Fig. 5.2-5.5 are given in Tab. 5.2 and Fig. 5.7-5.10 for E- and H-formulation with unnormalized as well as normalized vector basis functions. The conclusions for the two-dimensional case are 72

05.................... 0.75 0.75 Figure 5.6: Illustration of a rectangular cavity. E-formulation H-formulation TVFE Unnormalized Normalized Unnormalized Normalized Wh 8 12 23 31 Pe 2173 287 5939 958 Gr 684 99 1387 239 An 1834 195 5341 656. We 4238 472 11564 141.4 Table 5.2: Condition numbers for the global matrices resulting from FEM analysis of the cavity in Fig. 5.6. also seen to be valid in the three-dimensional case. 5.4 Inhomogeneous application of TVFEs In this section, condition numbers for FEM matrices are computed when mixedorder TVFEs of order 0.5 and 1.5 are combined selectively within a computational domain. Consider a rectangular cavity of normalized dimensions 1 x 0.75 x 0.5. It is composed of a rectangular cavity of dimensions 0.25 x 0.75 x 0.5 filled with a dielectric of relative permittivity er = 10 and a rectangular cavity of dimensions 0.75 x 0.75 x 0.5 filled with a dielectric of relative permittivity sr = 1, as illustrated in Fig. 5.11. The field and its derivatives within the cavity are expected to be largest in the dielectric and in the immediate vicinity hereof. An effective field expansion would 73

D E 1.5 = 1 0. 0 20 40 60 80 100 120 Element number Figure 5.7: Condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with unnormalized vector basis functions. E/H Normalized 3000 I - TVFE of order 0.5 - Whitney --- TVFE of order 1.5 - Peterson ----- TVFE of order 1.5 - Graglia 2500 j..........:. --- TVFE of order 1.5 - Andersen ---- TVFE of order 1.5 - Webb 2000 ~1500 - C: 500 0 20 40 60 80 100 120 Element number Figure 5.8: Condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with normalized vector basis functions. therefore use a higher order TVFE within the dielectric and in the part of the air-filled region closest to the dielectric and the lowest order TVFE elsewhere. The field continuity requirements between the lowest and higher order regions make hierarchical higher order TVFEs attractive and the detailed investigation of matrix conditioning in the previous section suggests using the hierarchical mixed-order TVFE of order ].5 74

E/H Unnormalized - -- C(Peterson)/C(Graglia)............................ C (A ndersen)/C (G raglia) -- - C(Webb)/C(Graglia) 12 j' 4 -* l l l l 0 0 20 40 60 80 100 120 Element number Figure 5.9: Ratios of condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with unnormnalized vector basis functions. E/H Nonrmalized 20c — * — C(Peterson)/C(Graglia) 18 -.................... C(Andersen)/C(Graglia)16 - C(Webb)/C(Graglia) 16 - -2 ' '" " " " '. 2 - 0 20 40 60 80 100 120 Element number Figure 5.10: Ratios of condition numbers for the individual element matrices for the cavity illustrated in Fig. 5.6; E/H-formulation with normalized vector basis functions. proposed in this dissertation rather than the one developed by Webb and Forghani. However, the field continuity requirements between the lowest and higher order regions only mandate use of a hierarchical higher order TVFE at the boundary to the lowest order region and hence away from this boundary an interpolatory higher or 75

-,r =10 0.5 1 0.75 Figure 5.11: Illustration of an inhomogeneous rectangular cavity. der TVFE could be used. When the field within an empty cavity is expanded using the mixed-order TVFEs of order 1.5 proposed in this dissertation as well as those developed by Savage and Peterson and Graglia et al. throughout the cavity, the interpolatory TVFE by Savage and Peterson leads to larger condition numbers than the hierarchical TVFE proposed in this dissertation while the interpolatory TVFE by Graglia et al. leads to smaller condition numbers than the hierarchical TVFE proposed in this dissertation. This suggests that combination of the hierarchical TVFE proposed in this dissertation with the interpolatory TVFE by Savage and Peterson within the higher order region will lead to larger condition numbers than those when the hierarchical TVFE proposed in this dissertation is used throughout the higher order region whereas combination of the hierarchical TVFE proposed in this dissertation with the interpolatory TVFE by Graglia et al. within the higher order region will lead to smaller condition numbers than those when the hierarchical TVFE proposed in this dissertation is used throughout the higher order region. Below, we examine the three different approaches for modeling the inhomogeneous cavity. They span the same space within each element and give identical eigenvalues. However, they are not equally efficient numerically as the different TVFEs provide matrices of different condition numbers.

Condition numbers are given in Tab. 5.3 for E- and H-formulation. The abbreviations 'Wh/An', 'Wrh/An/Pe' and 'Wh/An/Gr' denote the three different approaches. For both E- and 11 —formulation, the expected relative size of the different condition numbers is observed. Hence, an approach for improving the condition numbers of FEM matrices resulting from selective field expansion has been suggested and tested. The improvement comes at the expense of a more complicated formulation and computer code but does not alter accuracy. In essence, the suggested approach uses hierarchical TVFEs to transition between regions with lowest and higher order interpolatory TVFEs. Such use of transition TVFEs is directly equivalent to the use of scalar transition elements for efficient transition between coarse and fine meshes or between regions with lowest and higher order expansions [38, 87, 108]. This is a wellknown concept that has been applied successfully in mechanical and civil engineering for decades. TVFEs E-formulation H-formulation Wh/An 1188 654 Wh/An/Pe 1747 996 Wh/An/Gr 534 282 Table 5.3: Condition numbers for the global matrices resulting from FEM analysis of the cavity in Fig. 5.11. 5.5 Summary In this chapter, the condition numbers resulting from FEM analysis using the hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements 77

proposed in chapter 3 and chapter 4 were contrasted to those of existing interpolatory and hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements. The proposed hierarchical mixed-order TVFEs of order 1.5 proved better conditioned than existing hierarchical mixed-order TVFEs of order 1.5 and thus the analysis fostered no concerns for potential future convergence problems due to excessive matrix condition numbers. In addition, an approach for improving the condition numbers of FEM matrices resulting from selective field expansion was suggested and tested. The improvement comes at the expense of a more complicated formulation and computer code but does not alter accuracy. 78 P17

CHAPTER 6 ADAPTIVE TVFE REFINEMENT In this chapter, a review of existing error estimators and indicators is given and the effectiveness of the proposed hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra is investigated when some of the reviewed error indicators are applied in the context of a very simple adaptive refinement strategy. The work in this chapter is expected to be published in [4]. 6.1 Background Numerous methods for a posteriori error estimation or indication have been studied extensively in mathematics and engineering for decades and a vast amount of literature exists on the subject. For reviews of the various methods, see for instance [2, 16, 23, 46, 63, 74]. Following the approximate solution of a partial differential equation using a lowest order FE or hybrid FE/BI approach, each method seeks to identify local regions with large error for subsequent adaptive refinement of the mesh (h-refinement), basis functions (p-refinement) or mesh and basis functions simulta79

neously (hp-refinement) leading to an improved FE or FE/BI solution. This process is repeated until a desired accuracy is deemed to be reached. The identification of local regions can be performed by estimating the actual local error (error estimation) or by determining a local quantity that indicates whether the local error is small or large without estimating the actual local error (error indication). The local error estimation or indication can be carried out on an element by element basis or clusters of elements can be grouped together in sub-domains whereby the local error estimation or indication can be carried out on a sub-domain by sub-domain basis. Methods for which the error estimator or indicator can be computed directly from the initial solution are referred to as explicit methods whereas methods for which computation of the error estimator or indicator requires the solution of a local boundary value problem are referred to as implicit methods. The various methods for a posteriori error estimation or indication can be grouped in several different ways and hence a unique classification of these is not possible. The one presented here closely follows that of [74]. Implicit residual methods are based on the solution of a local Dirichlet or Neumann boundary value problem constructed from the lowest order FE solution [1, 76]. Explicit residual methods are based on local error estimation or indication by computation of a residual directly from the lowest order FE solution [10, 32, 34, 39, 47, 48, 93]. Among these, explicit complete residual methods take into account both local interior and boundary effects, explicit incomplete residual methods take into account only local interior effects and explicit interface residual methods take into account only local boundary effects. Recovery / gradient / average / smoothing methods are based on local comparison of the gradi80

ent of the original lowest order FE solution to a smoothened version of this gradient [109]. Dual / complementary / variational / mixed / hybrid methods are based on local comparison of solutions of two dual / complementary problems [16, 32, 33]. Perturbation methods are based on the estimation of local errors from differences between FE solutions of different orders [42, 57]. Interpolation and extrapolation methods are based on interpolation and extrapolation theory to compute approximations to higher order derivatives [31, 86]. For comparisons of the different methods, see for instance [16, 20, 24, 29, 33, 35, 41]. We note that the above classification is not complete. In addition to the main classes of methods presented above, other methods have been presented [74]. Hierarchical mixed-order TVFEs for tetrahedral elements have been proposed up to and including order 1.5 by Webb and Forghani [97], see section 2.2, and up to and including order 2.5 in this dissertation. The hierarchical mixed-order TVFE of order 1.5 proposed by Webb and Forghani was tested [76] for adaptive refinement using an implicit residual method. The purpose of this chapter is to investigate the merits of various explicit residual methods for a posteriori error indication. This is done via adaptive p-refinement of hybrid FE/BI solutions using the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedral elements proposed in this dissertation. We restrict ourselves to an adaptive p-refinement approach in this dissertation and refer to [22, 74, 82, 89] or additional references in [74] for adaptive h- or hp-refinement approaches. 81

6.2 Adaptive refinement strategies Consider a general three-dimensional electromagnetic problem. We perform a lowest order FE or FE/BI analysis (mixed-order TVFEs of order 0.5 applied for field expansion) with the computational domain discretized into N6 tetrahedral elements denoted by Te, = 1,., Ne, each having 4 faces denoted by F/, i = 1,-.,4. The center of Te and Fi is denoted by C(Te) and C(Fi), respectively., and the unit normal vector to Ft directed out of the element Te is denoted by fn. The lowest order FE or FE/1I solution leads to approximations of the electric fie-ld intensity E and the electric flux density D within and on the boundary of each element Te. On the face Fie we let D"in denote the value of D evaluated in Te and let Diout denote the value of D evaluated in the element bounding Te. Based on these, we present three different methods for indicating the error in a given region of the computational domain. For each method, we give the error indicator corresponding to an element. The corresponding error indicator for a sub-domain is simply the maximum of the error indicators for the elements comprising the sub-domain. The magnitude of the electric flux density generally does not correlate with the error associated with the electric flux density. Nevertheless, regions with high flux densities often give the dominant contributions to the physical response of a given electromagnetic eigenvalue, radiation or scattering problem. This justifies accurate modeling of such regions and hereby use of the simple error indicator EII(e)= [D(C(T6))[. (6.1) Although we are strictly not computing a residual, we will refer to this method as an 82

explicit incomplete residual method. The trivial three-dimensional extension of the two-dirnensional error indicator applied by Wang and Webb [93] for adaptive refinement in surface Mo)M problems is the error indicator EI2(e) = max [n (D" - Do )C(Fe)} (6.2) This method is an explicit interface residual method. A slightly different and computationally more expensive error indicator initially proposed by Golias and Tsiboukis [34, 35] is the error indicator E3(e) = max { >, ~ (D'in - Dut)2d. (6.3) This method is also an explicit interface residual method. In this dissertation, a conceptually very simple adaptive refinement strategy is adopted. Following the lowest order FE or FE/BI analysis, we determine the degree of error in each element via an error indicator and compute a refined solution where a certain pre-specified percentage of the elements having the highest degree of error are modeled with mixed-order TVFEs of order 1.5 and the remaining elements are again modeled with mixed-order TVFEs of order 0.5. A more advanced refinement strategy would estimate the optimal percentage of refinement for the improved solution, use TVFEs of more orders for refinement and incorporate a feedback loop leading to multiple error indications and refined solutions. However, given the lack of previous applications of adaptive refinement for practical electromagnetic problems, the simple adaptive refinement strategy described above was deemed sufficient in this dissertation. 83

6.3 Numerical results In this section, the merits of the error indicators EIji, ElI2 and EI3 presented in the previous section are investigated for determining input impedances of metallic patch antennas backed by material-filled cavities recessed in infinite metallic ground planes. We apply a standard hybrid FE/BI formulation with mixed-order TVFEs of different orders used for field expansion, as detailed in section 4.3.1. 6.3.1 Square patch antenna In section 4.3.2, we considered a square metallic patch antenna backed by a rectangular cavity recessed in an infinite metallic ground plane. In this section, we consider the same antenna and we therefore refer to the beginning of section 4.3.2 for a presentation of the antenna. We discretize the BI surface and the patch into 8 x 8 = 64 squares each of which are broken into two triangles. This surface mesh is extruded into the cavity to form a prism layer and each prism is broken into three tetrahedra. Four different TVFE options are applied: First, the mixed-order TVFE of order 0.5 is applied throughout the cavity. Second, the mixed-order TVFE of order 1.5 is applied throughout the cavity. Third, the mixed-order TVFE of order 0.5 is applied in conjunction with the mixed-order TVFE of order 1.5 in the vicinity of the radiating edges of the patch (38% of the TVFEs are of order 1.5). Fourth, the mixed-order TVFE of order 0.5 is applied in conjunction with the mixed-order TVFE of order 1.5 in regions found adaptively (40% of the TVFEs are of order 1.5). The refinement is carried out for 84

each of the error indicators EIl, EI2 and EI3 given by (6.1)-(6.3) on a tetrahedron by tetrahedron (384 elements), prism by prism (128 sub-domains) as well as brick by brick (64 sub-domains) basis. The nine different cases are defined in Tab. 6(.1. Error Base of adaptive Percentage of Case indicator refinement refinement 1 EI1 Tetra 40 2 EI1 Prism 40 3 EIL Brick 40 4 E12 Tetra 40 5 EI2 Prism 40 6 EI2 Brick 40 7 EI3 Tetra 40 8 EI3 Prism 40 9 EI3 Brick 40 Table 6.1: Definition of Case 1-9. For the mixed-order TVFEs of order 0.5 and 1.5, the particular mesh is too coarse to yield the correct resonant frequency of 4.43 GHz as obtained by Schuster and Luebbers [81] and confirmed in chapter 4 for finer meshes. Nevertheless, the mesh is very useful for evaluating the merits of the various error indicators. The dynamic range of the differences in resonant frequency between solutions where mixed-order TVFEs of order 0.5 and 1.5 are applied throughout the computational domain are the largest for coarse meshes and hence coarse meshes are ideal for investigating how well shifts in resonant frequency are predicted by adaptively refined solutions based on various error indicators. This approach can of course only be justified for problems where proper convergence has been ensured by very accurately predicting the correct resonant frequency by using a more accurate approach (finer meshes, higher order TVFEs). As demonstrated in section 4.3.2, such convergence was observed for this 85

particular antenna. Real and imaginary parts of the input impedance as a function of frequency are given in Fig. 6.1-6.9 for each of the four TVFE options for Case 1-9. The results for the first three TVFE options are the same for Case 1-9. The resonant frequency predicted with the mixed-order TVFE of order 1.5 applied throughout the cavity is larger than that when the mixed-order TVFE of order 0.5 is applied throughout the cavity. Application of the mixed-order TVFE of order 1.5 in the vicinity of the radiating edges only and the mixed-order TVFE of order 0.5 elsewhere is seen to predict this shift very well. Most importantly, we observe that adaptive refinement on an element by element basis (Case 1, 4, 7) is seen to predict an inaccurate upward shift while adaptive refinement on a sub-domain by sub-domain basis (Case 2-3, 5 -6, 8-9) is seen to predict the upward shift very well. For EI2 and E13, prism by prism refinement is slightly more accurate than brick by brick refinement but more complicated antennas must be studied before decisive conclusions can be reached regarding the relative merits of these two sub-domain refinement schemes. Corresponding regions of refinement at the resonant frequency are displayed in Fig. 6.10-6.18. These figures show a top view of the mesh where the cavity and patch boundaries are marked with thick black lines. A white / light gray / average gray / dark gray triangle indicates that 0 / 1 / 2 / 3 of the three elements in the prism beneath the triangle are being refined with a mixed-order TVFE of order 1.5. Adaptive refinement on an element by element basis is seen to lead to very sporad(lic regions where the mixed-order TVFE of order 1.5 is applied. Such combination of lowest and higher order TVFEs is inefficient and consequently leads to less accurate 86

400, 1 I I I I 300 8 0 200u Nbc 200 01 3~-0 -200 --- 0.5 - - - - 1.5 - - 0.5+1.5(edges) 0.5+1.5(adaptive) 95 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4. Frequency [GHz] -4 — 0.5 -* — 1.5 --- 0.5+1.5(edges) 0.5+1.5(adaptive);requency [GHz] Figure 6.1: Real and imaginary part of input impedances for Case 1 (El,, Tetra). 4001 1 1 1 1 1 I I I - 300 8 0 N100 200 '0' 100 0 - 100 -200 -e- 0.5+1.5(edges) 0.5+1.5(adaptive) 05 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4~ Frequency [GHz] -4- 0.5 1.5 ae- 0.5+1.5(edges) 0.5+1.5(adaptive).5 3.6 317 3.8 3.9 4 4.1 4.2 4.3 4.4 4. Frequency [GHz] 5 Figure 6.2: RearI and imaginary part of input impedances for Case 2 (El,, Prism). 40 0 r -1 1 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 300 0 i200 N100 -4- 0.5 - -l — 1.5 -e- 0.5+1.5(edges) -a- 0.5+1.5(adaptive) - - - - -. - - - - -.5 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4 Frequency [GHz] 208........ 3.1 0 3.6.... 373.3..1.. 3 4 4 FrqecEGz Fiue.: el n iainr pr o nptimeane frCae3 E1,Bic) 08

400, N I I I I I I 300 S E 200 N1 100 0 — 0.5 K --- 1.5 --- '0.5+1.5(edges) 0.5+1.5(adaptive) _J r - -- - - ~ * - -- 5 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4 Frequency [GHz] Figure 6.4: Real and imaginary part of input impedances for Case 4 (El2, Tetra). 400,,,, 1 I I I 30C E -C 20 rr 100 ---- 0.5: 1.5 - 0.5+1.5(edges) 0.5+1.5(adaptive).. X,, ^.;.5 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 Frequency [GHz] Figure 6.5: Real and imaginary part of input impedances for Case 5 (EI2, Prism). 3.5 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 Frequency [GHz] Figure 6.6: Real and imaginary part of input impedances for Case 6 (EI2, Brick). 88

4001-,, I I I I I I 300 0 E 0 I200 N0 100 - 0.5 - * — 1.5 -E 0.5+1.5(edges)\ - 0.5+1.5(adaptive) \ \ \\ rA. a - - -. 1 3 5 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 Frequency [GHz] 4.5 Figure 6.7: Real and imaginary part of input impedances for Case 7 (EI3, Tetra). 400,....... 300 E 200 100 1001 - - 0.5 A 1.5 -E- 0.5+1.5(edges) - 0.5+1.5(adaptive). \ IL.I 5 3.6 3.7 3.8 3.9 4 4.1 4.2 Frequency [GHz] 4.3 4.4 4 5 Figure 6.8: Real and imaginary part of input impedances for Case 8 (EI3, Prism). Frequency [GHz] Frequency [GHz] Figure 6.9: Real and imaginary part of input impedances for Case 9 (El, Brick). 89

Figure 6.10: Regions of refinement for Case 1 (El1, Tetra) at 4.15 GHz. Figure 6.11: Regions of refinement for Case 2 (EI1, Prism) at 4.30 GHz. Figure 6.12: Regions of refinement for Case 3 (Ell, Brick) at 4.30 GHz. 90

Figure 6.13: Regions of refinement for Case 4 (EI2, Tetra) at 4.10 GHz. Figure 6.14: Regions of refinement for Case 5 (EI2, Prism) at 4.30 GHz. Figure 6.15: Regions of refinement for Case 6 (EI2, Brick) at 4.30 GHz. 91

~_ ~_ ~_ ~_ ~_ ~ Figure 6.16: Regionm s of refinement for Case 7 (EI3, Tetra) at 4.10 GHz. Figure 6.17: Regions of refinement for Case 8 (EI3, Prism) at 4.30 GHz. Figure 6.18: Regions of refinement for Case 9 (EI3, Brick) at 4.30 GHz. 92

results, as observed. Adaptive refinement on a sub-domain by sub-domain basis inherently circumvents this problem and leads to very accurate results, as observed. We note that all adaptive refinement schemes were tested for less than 40% refinement. 30% refinement was found to lead to very accurate results with EI but slightly less accurate results with EI2 and EI3. Generally, 20% refinement was not enough to significantly improve the lowest order solution. For brevity, these results are not included. 6.3.2 Printed bowtie antenna Consider a metallic printed bowtie antenna backed by a rectangular cavity recessed in an infinite metallic ground plane, as illustrated in Fig. 6.19 (side view) and Fig. 6.20 (top view). The cavity-backed patch antenna is situated in free space characterized by the permittivity so and the permeability,uo. The cavity dimensions are 48 cm x 64 cm x 12 cm. The interior metallic cavity walls are covered with an AA of permittivity (1 -j2.7):o, permeability (1 -j2.7)[to and thickness 4 cm. This absorber layer would not be part of an actual antenna but is merely a well-established computational tool [64] that serves to approximately simulate a metallic printed bowtie antenna situated in free space. Such an antenna is expected to exhibit broadband behavior in contrast to the narrowband square patch antenna considered above. Consequently, the printed bowtie antenna can provide a second and independent evaluation of the error indicators presented in section 6.2. The specific printed bowtie antenna consists of two isosceles triangular patches characterized by the opening angle 67.38~ and the 93

~o ho 4cm 8cm 24cm 8cm 4cm 0o 1o 8 cm I(j2.7) (I-j2 4cm Figure 6.19: Side view of a metallic printed bowtie antenna backed by an air- and absorber-filled rectangular cavity recessed in an infinite metallic ground plane. --- -- --- -- --- -- --- -- --- 1 1 8 cm 24 cm 8 cm i i:A 18 cm I: I 18 cm absorber-filled rectangular cavity recessed in an infinite metallic ground plane. maximum width 24 cm. The bowtie patches are centered in the cavity aperture and fed by a probe of constant current connecting the two triangular patches. An antenna very similar to the above was discussed by Collin [18] and is expected to cover the UHF channels 14 to 83 spanning the frequency range [450 MHz, 900 MHz] when used with a 300 Q feed line. That is, the real and imaginary parts of the input impedance are expected to hover around 300 Q and 0 Q. respectively, in this frequency range.!!:::1 Is 32~~~~'g 71';i~~>~ol'~IlgC. ysmla ote bv vs icse b oln 1]adisepce 94

The distance from the antenna, to the absorber is 8 cm which translates into Ao/8.33 at 450 MHz and Ao/4.17 at 900 MHz. Similarly, the thickness of the absorber is 4 cm which translates into Ao/16.67 at 450 MHz and Ao/8.33 at 900' MHz. These distances should be sufficient to simulate an antenna in free space. That the absorber actually works is confirmed by the fact that the behavior of the antenna is significantly altered if the absorber is removed and the antenna is backed by a purely metallic cavity. The main purpose of the absorber is to prevent excitation of modes norral to and bounded by the metallic antenna and the bottom of the metallic cavity. If such modes are excited, the antenna ceases to operate as an antenna situated in free space but starts to function in the presence of the cavity which is not intended here. The absorbers covering the side walls are not essential for this purpose and therefore can be expected to have less influence on the behavior of the antenna. Indeed., it was confirmed that results strikingly similar to those presented in the following are obtained when the absorbers covering the side walls are removed. The extent of the feed region is 4 cm (Ao/16.67 at 450 MHz and Ao/8.33 at 900 MHz) and hence it is reasonable to assume a constant current over the probe in the frequency range in which the antenna is operated. For larger feed regions, a phase variation of the probe current will have to be accounted for. Doing so is trivial within the context of the FE/BI method but it was deemed unnecessary in this case. For analysis and evaluation of the proposed error indicators, we discretize the BI surface and the patch into a coarse (56 antenna triangles and 330 BI triangles) and a fine (160 antenna triangles and 1030 BI triangles) mesh. These surface meshes are extruded into the cavity to form three prism layers. Each prism is broken into 95

three tetrahedra and four different TVFE options are applied: First, the mixed-order TVFE of order 0.5 is applied throughout the cavity for the coarse mesh. Second, the mixed-order TVFE of order 0.5 is applied throughout the cavity for the fine mesh. Third, the mixed-order TVFE of order 1.5 is applied throughout the cavity for the coarse mesh. Fourth, the mixed-order TVFE of order 0.5 is applied in conjunction with the mixed-order TVFE of order 1.5 in regions found adaptively (20% of the TVFEs are of order 1.5) for the coarse mesh. Supported by the findings for the square metallic patch antenna, the refinement is carried out for each of the error indicators EIl, EI2 and El3 given by (6.1)-(6.3) on a sub-domain by sub-domain basis only with a sub-domain being three adjacent prisms extruded from a surface triangle (386 sub-domains for the coarse mesh). That is, we opt not to examine adaptive refinement on an element by element basis since it did not prove efficient in the previous analysis and we further limit ourselves to only one type of sub-domain by sub-domain refinement. The three different cases are defined in Tab. 6.2. Error Base of adaptive Percentage of Case indicator refinement refinement 10 E,1 Prisms (3 adjacent) 20 11 EI2 Prisms (3 adjacent) 20 12 EI3 Prisms (3 adjacent) 20 Table 6.2: Definition of Case 10-12. Real and imaginary parts of the input impedance as a function of frequency are given in Fig. 6.21 for the first three TVFE options (no adaptivity). Instead of the strongly resonant behavior characterizing cavity-backed patch antennas, we observe the expected slightly oscillatory behavior of the real and imaginary parts. Mixed 96

order TVFEs of order 0.5 for the fine mesh and mixed-order TVFEs of order 1.5 for the coarse mesh give similar results that are better (real parts closer to 300 Q and imaginary parts closer to 0 Q) than those with mixed-order T''VFEs of order 0.5 for the coarse mesh. Although mixed-order TVFEs of order 0.5 for the fine mesh and mixedorder TVFEs of order 1.5 for the coarse mesh give similar results, the latter approach is more attractive than the former in terms of memory and CPU time requirements. To demonstrate the merits of adaptive refinement, Fig. 6.22 shows the real and imaginary parts of' the input impedance as a function of frequency for the coarse mesh with mixed-order TVFEs of order 1.5 throughout the cavity and with 80% mixed-order TVFEs of order 0.5 and 20% mixed-order TVFEs of order 1.5 found via adaptive refinement using EIl, EI2 and EI3 on a sub-domain by sub-domain basis (Case 10-12). The results are almost indistinguishable expressing that we can accurately predict the behavior of the antenna using only 20% mixed-order TVFEs of order 1.5. This presents a significant memory and CPU time improvement at virtually no cost. To illustrate the regions of refinement for Case 10-12, we consider a cross-section of the mesh parallel to the antenna with the boundaries of the metallic cavity and antenna, marked with thick black lines. As mentioned previously, the tetrahedral volume mesh is grown from a triangular surface mesh in a cut in the plane of the bowtie patches by extruding it into three prism layers and breaking each prism into three tetrahedra. Since we consider sub-domain by sub-domain refinement only, 0 or 9 tetrahedra can be refined corresponding to a given triangle in the cross-section. In the following, a white / dark gray triangle indicates that 0 / 9 tetrahedra are being 97

400 300 0 N 0- 0.5 (coarse) r 100 - -- 0.5 (fine) 0.4 0.5 0.6 0.7 0.8 0.9 1 Frequency [GHz] 300 0- 0.5 (coarse) 200............. - 0.5 (fine)........ E - - - E 1.5 (coarse) 0 100 N E 0 -10). 0.4 0.5 0.6 0.7 0.8 0.9 1 Frequency [GHz] Figure 6.21: Real and imaginary part of the input impedance for the metallic bowtie patch antenna in Fig. 6.19-6.20. 300 - 100 - -- 20%- E12 E_- 20% - E13. -100 0.4 0.5 0.6 0.7 0.8 0.9 1 Frequency [GHz] 3001!i 0; - 100%.................:............ 2 0 % - E ll................ 100 -E3!-. -...... 0.4 0.5 0.6 0.7 0.8 0.9 1 Frequency [GHz] Figure 6.22: Real and imaginary part of the input impedance for the metallic bowtie patch antenna in Fig. 6.19-6.20. refined. The regions of refinement at 0.7 GHz (close to the center of the frequency band of operation) are shown in Fig. 6.23 for Case 10-12. Although they are seen to differ slightly for the three different error indicators, they all show the general trend of predicting the feed area and, to a lesser extent, the corners of the triangular patches as the regions where mixed-order TVFEs of order 1.5 are applied. We note that all approaches predict almost the same far field patterns a.s well as a 98

Figure 6.23: Regions of refinement for the metallic bowtie patch antenna in Fig. 6.19 -6.20 for Case 10 (left), Case 11 (middle) and Case 12 (right). linearly polarized far field (parallel to the feed) in the direction normal to the bowtie patches. For brevity, these expected results are not included. 6.3.3 Discussion The error indicator EIl is an explicit incomplete residual indicator taking into account only local interior effects while the error indicators EI2 and EI3 are explicit interface residual indicators taking into account only local boundary effects. Explicit complete residual indicators take into account both local interior and boundary effects and thus are expected to offer superior performance. An obvious construction of an explicit complete residual indicator is to simply combine EIl with E12 or E'3. That is, a certain fraction of the elements or sub-domains are refined based on El11 and the remaining elements or sub-domains are refined based on EI2 or EI3. This approach was tested and found accurate for analysis of square metallic patch antennas. However, more complicated antennas must be studied before decisive conclusions can be reached regarding the efficiency of this approach. 99

An apparent drawback of adaptive refinement is that it requires multiple (in this case two) FE/BI solutions for each frequency. This seems computationally expensive since it involves the consecutive solution of different equation systems using an iterative solver for which zero is the traditional starting guess. A possible remedy might be use of the previous FE/BI solution at a given frequency as the starting guess for the iterative solver all adaptive FE/BI solutions. As stated by Golias and Tsiboukis [33], this should greatly accelerate the convergence of the adaptive FE/BI solutions. We note that each adaptive FE/BI solution requires a starting guess for numerous unknown expansion coefficients that have not previously been solved for. For these unknowns, we shall continue to use zero as the starting guess. To investigate the effectiveness of this approach for the metallic printed bowtie antenna considered above, Fig. 6.24 shows the final number of iterations to reachl a 10-3 relative residual with a QMR solver as a function of frequency. Convergence curves are given for the approaches where TVFEs of order 0.5 are used for a coarse mesh, TVFEs of order 0.5 are used for a fine mesh and TVFEs of order 1.5 are used for a coarse mesh. These three cases correspond to 3704, 11825 and 19745 unknowns, respectively. The number of iterations is seen to be slightly under 500 for all frequencies when mixed-order TVFEs of order 0.5 are applied for the coarse mesh. As expected, this number increases when a finer mesh or mixed-order TVFEs of order 1.5 are applied. It increases more in the latter case but relative to the number of unknowns this case has the best convergence properties. Fig. 6.25-6.27 shows the same information for Case 10-12 with zero as well as the solution with mixed-order TVFEs of order 0.5 used as the starting guess for the 100

9rnnI N" - 0.5 (coarse)., + 0.5 (fine) 2 1.5 (coarse) | 15000 -, 500 0.4 0.5 0.6 0.7 0.8 0.9 1 Frequency [GHz] Figure 6.24: Final number of iterations as function of frequency for different iterative solutions. iterative solver. The number of unknowns is around 6400 for all frequencies and error indicators. The number of iterations is lower than when mixed-order TVFEs of order 1.5 are applied throughout the coarse mesh. Also, the number of iterations is consistently seen to drop when the solution with mixed-order TVFEs of order 0.5 is used as the starting guess for the iterative solver instead of zero. The savings come at no cost since we need to carry out the solution with mixed-order TVFEs of order 0.5 in order to determine the regions of refinement. However, even if we pre-select the higher order regions and thus have no need to find a solution with mixed-order TVFEs of order 0.5, it is still worthwhile to do so. Subsequent use of such a solution as the starting guess for the iterative solver can reduce the iteration count by several hundred iterations for systems with approximately 6400 unknowns for the price of less than 500 iterations for systems with 3704 unknowns. Fig. 6.28-6.30 shows the relative residual as a function of the iteration number at the frequency 0.7 GHz (close to the center of the frequency band of operation) for Case 10-12. The results reinforce the conclusions from above. 101

en C: 0._ o ca z E= 2500 Zero Order 0.5 solution 2000 150d 1000 500 i1 6.4 0.5 0.6 0.7 0.8 0.9 1 Frequency [GHz] Figure 6.25: Final number of iterations as a function of frequency for adaptive solutions for Case 10 (EIl, Prism) with different starting guesses for the iterative solver. 9.(nn II 2000 o 'Z- 1500 150C a) | 100C z 50C - Zero Order 0.5 solution:: 0.4 0.5 0.6 0.7 Frequency [GHz] 0.8 0.9 Figure 6.26: Final number of iterations as a function of frequency for adaptive solutions for Case 11 (EI2, Prism) with different starting guesses for the iterative solver.,n(..... 2000 I 1500< o E 1000 z - Zero - Order 0.5 solution * 0 J:F -4: ' -- '-.:: 500 i I b.4 0.5 0.6 0.7 0.8 Frequency [GHz] 0.9 Figure 6.27: Final number of iterations as a function of frequency for adaptive solutions for Case 12 (EI3, Prism) with different starting guesses for the iterative solver. 102

10~.-o E3 a) a1) nrr 500 1000 Iteration number 1500 Figure 6.28: Relative residual as a function of the iteration number for adaptive solutions for Case 10 (El1, Prism) with different starting guesses for the iterative solver. 10~ 1 0.=====................................ Zero........... Order 0.5 solution 10 -~..........:... co 10-3 10-4 l-l 0 500 1000 1500 Iteration number Figure 6.29: Relative residual as a function of the iteration number for adaptive solutions for Case 11 (EI2, Prism) with different starting guesses for the iterative solver. 1 0::::::::::.:::::::::: Zero... - - - Order 0.5 solution: 1 0 10-..................... 104 0 500 1000 1500 Iteration number Figure 6.30: Relative residual as a function of the iteration number for adaptive solutions for Case 12 (EI3, Prism) with different starting guesses for the iterative solver. 103

We note that the above comments for reducing the iteration count do not apply to resonant antennas like the square metallic patch antenna considered previously. This is due to the fact that the field values throughout the computational domain experience a strong frequency dependence. 6.4 Summary In this chapter, a review of existing error estimators and indicators was given and the effectiveness of the proposed hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra, was investigated when some of the reviewed error indicators were applied in the context of a very simple adaptive refinement strategy. The results were extremely promising for both narrowband and broadband antennas provided the refinement was carried out on a sub-domain by sub-domain basis as opposed to an element by element basis. 104

CHAPTER 7 ANALYSIS OF TAPERED SLOT ANTENNAS In this chapter, the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra proposed( in this dissertation are used in conjunction with a simple adaptive refinement strategy to analyze the impedance and pattern characteristics of tapered slot antennas (TSAs) using the hybrid FE/BI method. The work in this chapter is expected to be published in [4]. 7.1 Background Metallic TSAs such as exponentially tapered slot antennas (ETSAs, also called Vivaldi antennas), linearly tapered slot antennas (LTSAs) and constant width slot antennas (CWSAs), see Fig. 7.1, are of practical interest as reflector antenna or lens feeds, as elements of broadband antenna arrays or (for larger TSAs) directly as transmit or receive antennas [106]. Whether isolated or backed by thin dielectrics, TSAs have a number of appealing properties. They are broadband with typical impedance and pattern bandwidths of 5:1 [83] although values as large as 40:1 have 105

been reported [52], they can provide an almost symmetric main beam despite their entirely planar nature [44] with typical -10 dB beamwidths of 30-40~ [106i], they have a moderate directivity of approximately ten times the antenna, length normalized to the free space wavelength [107] and they have reasonable cross-polarization characteristics with a typical isolation of -20 dB in the principal planes [44]. Their beamwidth and directivity are somewhat controllable when the length of the antenna is varied [106] and a narrow pattern beamwidth versus broad pattern bandwidth trade-off exists [107]. In addition, they have the practical advantages that they are inexpensive and (feed region excluded) easy to fabricate [106], that their transmitting or receiving portion is well separated from the necessary circuitry providing ample space for the latter [106] and that the transverse spacing between TSA array elements can be made very small [106]. Figure 7.1: Illustration of a metallic (dark gray) ETSA (left), LTSA (middle) and CWSA (right) backed by a dielectric (light gray). TSAs are traveling wave antennas (TWAs) of the surface wave type, i.e. a traveling wave propagates along the antenna structure with a phase velocity smaller than the speed of light. A review of TWAs was given by Zucker [110]. During the past two decades, the characteristics of TSAs have been measured as well as simulated using 106

approximate methods and more rigorous numerical techniques such as the MoM, the FDTD method and the transmission line method. TSAs in air have been characterized in [43, 44, 19] while TSAs backed by thin homogeneous dielectrics to better confine the field to the slot and thereby provide narrower beamwidths have been treated in [30, 43, 44, 49, 106, 107]. TSAs backed by inhomogeneous dielectrics (homogeneous dielectrics with holes) to lower the effective dielectric constant and provide a directivity improvement have been measured [58] and simulated using the FD'I'D method [17]. The bandwidth bottleneck of the ETSA is the transition between the feed and slot region [52]. The antipodal ETSA was introduced in [28] and comprehensively studied in [26] to overcome this difficulty while still providing beamwidths and directivity comparable to the ETSA. The antipodal ETSA, however, gives rise to a high degree of cross-polarization and experiences a significant polarization tilt for high frequencies [52]. The balanced antipodal ETSA introduced in [52] removes these limitations while maintaining the attractive features of the antipodal ETSA. Parametric studies of various TSAs are given in [26, 51, 83, 88]. TSAs are often used as elements of large phased arrays [19, 78, 79, 80, 83, 88]. When used as array elements, TSAs are usually much smaller than stand-alone TSAs and in many cases so small that they cease to work as broadband TWAs [107]. Such small TSAs are useless as stand-alone antennas but can be useful within an array framework [80]. Accurate determination of the phase velocity of the traveling wave within the tapered slot is important for accurate prediction of the beamwidth in the E-plane and, especially, the H-plane of a TSA [44]. This suggests that TSAs can be more 107

accurately characterized by improving the field modeling within the slot. It is the aim of this chapter to investigate the validity of this hypothesis for TSAs situated in free space. To this end, the TSAs are placed within a metallic cavity whose bottom and side walls are covered with an AA that, as discussed in section 6.3, would not be part of an actual antenna but is merely a well-established computational tool to simulate an antenna situated in free space. After a validation against results and trends reported in the literature, the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra proposed in this dissertation are used in conjunction with a simple adaptive refinement strategy to analyze the impedance and pattern characteristics of a LTSA using the hybrid FE/BI method. 7.2 Tapered slot antenna analysis We consider a LTSA uniquely characterized by the height 15.46 cm, the width 2.68 cm + 3.98 cm + 2.68 cm = 9.34 cm, the opening angle 14~ and the narrowest slot width 0.18 cm. The antenna operates at 10 GHz and is fed by a probe of constant current situated at the narrow end of the tapered slot, as illustrated in Fig. 7.2. The probe excites a traveling wave in the tapered slot resulting in endfire radiation. Apart from some minor geometrical differences that do "not change the physics that governs the radiation mechanism of the antenna" [43], this LTSA was measured and simulated using the MoM in [43] 1. e use the F/B method to simulate the LTSA situated in free space by placing it 1The LT'SA in [43] is actually a slightly larger antenna operating at 9 GHz. It is the same size in wavelengths as the one considered here. 108

2.68 cm 3.98 cm 2.68 cm 15.46 cm 0.18 cm Figure 7.2: Geometrical parameters for a LTSA (not drawn to scale). in a rectangular cavity with an open top surface and metallic bottom and side surfaces. The interior metallic cavity walls are covered by an AA of relative permittivity and permeability 1 - j2.7 and thickness 0.15 cm. Thicker absorbers were found to give similar results. The top of the antenna (the wide end of the tapered slot) is aligned with the open boundary of the cavity. The vertical distance from the bottom of the cavity to the bottom of the antenna is 1.65 cm while the horizontal distance from the side wall of the cavity to the antenna is 1.65 cm parallel to the antenna and 2.45 cm perpendicular to the antenna. The perpendicular dimension is larger since it is the main purpose of the absorber to prevent excitation of modes normal to and bouncded by the metallic antenna and the metallic cavity walls parallel to the antenna. If such modes are excited, the antenna ceases to operate as an antenna situated in free space but starts to function in the presence of the cavity which is not intended here. Such operation can, however, be useful for arrays where suppression of scan blindness is of importance [101]. The measured and MoM E-plane (parallel to the antenna) and H-plane (perpendicular to the antenna) far field patterns (polar cuts) reported in [43] were scanned 109

and digitized and are given in Fig. 7.3-7.4 along with FE/BI1 results using the mixedorder TVFEs of order 0.5 for a coarse mesh (129492 elements, 264513 faces and 159460 edges). Since the geometrical symmetry of the LTSA ideally translates into symmetric E- and H-plane patterns, only half of the E- and H-plane patterns are actually given in [43]. The measured and MoM patterns from [43] were therefore mirrored for comparison with the full E- and H-plane patterns found via the FE/BI method as given in Fig. 7.3-7.4. However, had the full patterns been measured and computed in [43], minor discrepancies due to measurement and mesh asymmetries would have been present. Also, since the measured and MoM data in [43] are normalized to 0 dB at endfire (0 = 0), the same normalization is applied to the FE/BI data. Overall, the FE/BI data agree well with the measured and MoM data in both the E- and the H-plane. The FE/BI main beam in the E-plane is slightly narrower than the measured and MoM main beams while the FE/1BI main beam in the Hplane is almost identical to the measured main beam. This is quantified by the approximate E-plane -10 dB beamwidths 22~ (FE/BI), 23~ (MoM) and 26~ (measured) and the approximate H-plane -10 dB beamwidths 370 (FE/BI), 35~ (MoM) and 36~ (measured). In the E-plane, the measured data have a shoulder at 30 - 40~ and a null at 50~ while the MoM data have nulls at 30~ and 50~. The FE/BI data fall in between with a null at 40~. The MoM and measured data stay below -15 dB beyond 50~ whereas the FE/BI data stay below -12 dB for these angles. In the H-plane, the FE/BI data follow the MoM and measured data very well down to -20 dB (50~) at which point the measured and MoM data level off while the FE/BI data continue to decrease. 110

E-plane 0 0 [degrees] Figure 7.3: E-plane patterns for the LTSA in Fig. 7.2. H-plane -- Li 0 0 [degrees] Figure 7.4: H-plane patterns for the LTSA in Fig. 7.2. 111

Considering the fact that the above FE/BI results were found using mixed-order TVFEs of order 0.5 with a coarse discretization, better overall agreement cannot be expected. The results show acceptable agreement with independent measured and MoM data, and display trends (such as an almost symmetric main beam with the H-plane being somewhat broader than the E-plane) that agree with the literature [43, 49, 106]. However, several other trends that we can confirm with the above antenna but cannot confirm from the above plots have been reported for LTSAs. Before proceeding to demonstrate the merits of hierarchical mixed-order TVFEs for LTSA modeling, we shall therefore discuss some of these trends to provide further validation. LTSAs have very low cross-polarization in the principal E- and H-planes and larger cross-polarization in the diagonal D-plane [49]. The actual values of the crosspolarization in the D-plane will inevitably depend on how co- and cross-polarization is defined as this definition is not unique [55] 2. Following Ludwig's Definition III [55], it is demonstrated in [49] that for a LTSA backed by a dielectric, the cross-polarized field in the D-plane is significantly smaller than the co-polarized field close to the endfire direction 0 = 0 but rises to a level comparable to and in some cases higher than the co-polarized field off-axis. To investigate the polarization characteristics of the above antenna, we plot the co- and cross-polarized patterns as defined by Ludwig's Definition III [55] in the E-, H- and D-planes, see Fig. 7.5-7.7. They are seen to be in full agreement with the trends reported in [49]. The lack of symmetry for the 2Ludwig's Definitions II and III [55] suitable for antenna patterns lead to the same definitions for co- and cross-polarization in the E- and H-planes but different definitions in the D-plane. 112

cross-polarized fields in the E- and H-planes is not unexpected and is due to the fact that the cross-polarized field values in these planes are very small and thus extremely sensitive to any asymmetries introduced by the FE/BI analysis process (asymmetrical meshes, integrations and so on). We note that similar patterns were not given in [43]. LTSAs are inherently broadband antennas with typical bandwidths up to 5:1 [83]. Using only one mesh for accurate simulation over such a wide range of frequencies is inefficient. It is necessary to keep the distance to the absorber large enough in terms of wavelengths at low frequencies (large wavelengths) which in combination with the requirement of a discretization rate small enough in terms of wavelengths at high frequencies (small wavelengths) results in an unmanageable computational domain. The solution is to re-mesh the geometry so the analysis at each frequency is carried out with a mesh that provides approximately the same electrical distance to the cavity wall and approximately the same electrical discretization rate. Since it is not the purpose of this chapter to investigate LTSA bandwidth limitations but rather to demonstrate the effectiveness of modeling LTSAs using a multi-resolution FE/BI approach, this re-meshing approach was not adopted for the LTSA in Fig. 7.2. However, the antenna in Fig. 7.2 was simulated at 8 GHz, 9 GHz, 11 GHz and 12 G- Hz with the mesh used at 10 GHz to show essentially similar patterns at all frequencies. Although this does not demonstrate the extreme bandwidths reported in the literature [52], sufficient bandwidth for validation purposes has been observed. To demonstrate the merits of hierarchical mixed-order TVFEs for multi-resolution FE/BI modeling of LTSAs, extensive numerical simulations are carried out. We therefore opt to consider a slightly smaller antenna and cavity than above. Specifically, we 113

-30 -. -: -- -' -50 -60 o-pol 1,- _Cross-pol -70 -80 -60 -40 -20 0 20 40 60 80 e [degrees] Figure 7.5: Co- and cross-polarized E-plane patterns for the LTSA in Fig. 7.2. H-plane -10 -20 -30 -40 -50 -60 - -60 Co-pol l Cross-pol -70 -80 -60 -40 -20 0 20 40 60 80 e [degrees] Figure 7.6: Co- and cross-polarized H-plane patterns for the LTSA in Fig. 7.2. D-plane mo" - --. -30,; -40.. -6-0 - Co-polt Cross-pol -70 -80 -60 -40 -20 0 20 40 60 80 0 [degrees] Figure 7.7: Co- and cross-polarized D-plane patterns for the LTSA in Fig. 7.2. 114

consider a LTSA uniquely characterized by the height 9.00 cm, the width 1.50 cm + 2.55 cm + 1.50 cm = 5.55 cm, the opening angle 15.20 and the narrowest slot width 0.15 cm, as illustrated in Fig. 7.8. As above, the antenna operates at 10 GEHz and is fed by a probe of constant current situated at the narrow end of the tapered slot, see Fig. 7.8, for excitation of a traveling wave in the tapered slot. For FE/BI analysis of this antenna situated in free space, we place it in a rectangular cavity with an open top surface and metallic bottom and side surfaces. The interior metallic cavity wails are covered by an AA of relative permittivity and permeability 1-j2.7 and thickness 0.15 cm and the top of the antenna (the wide end of the tapered slot) is aligned with the open boundary of the cavity. The vertical distance from the bottom of the cavity to the bottom of the antenna as well as the parallel and perpendicular horizontal distances from the side walls of the cavity to the antenna are all 1.65 cm. 1.50cm 2.55 cm 1.50cm 9.00 cm 0. 15cm Figure 7.8: Geometrical parameters for a LTSA (not drawn to scale). For analysis of the LTSA, we employ a coarse and a dense tetrahedral volume mesh grown from a coarse and a dense triangular surface mesh in a cut in the plane of the antenna. Each surface mesh is extruded into a number of prism layers (ten for the coarse mesh and fourteen for the dense mesh) and each prism is broken into three 115

tetrahedra. The coarse mesh has 44370 elements, 91269 faces and 55621 edges while the dense mesh has 148638 elements, 303013 faces and 182109 edges. Four different TVFE options are applied for FE/BI analysis: First, the mixed-order TVFE of order 0.5 is applied throughout the coarse mesh. Second, the mixed-order TVFE of order 0.5 is applied throughout the dense mesh. Third, the mixed-order T VFE of order 0.5 is applied within the coarse mesh in conjunction with the mixed-order TVFE of order 1.5 in the four prism layers closest to the metallic antenna (40% of the elements within the mesh) where the dominant fields are expected and accurate field modeling is therefore necessary. Fourth, the mixed-order TVFE of order 0.5 is applied within the coarse mesh in conjunction with the mixed-order TVFE of order 1.5 in regions found adaptively (Case 1-4). The refinement is carried out by confining the higher order TVFEs to the two (Case 1-2) or four (Case 3-4) prism layers closest to the metallic antenna and performing the refinement with the error indicator El1 (shown in chapter 6 to be at least as good as EI2 or EI3) on an element by element (Case 1,3) or sub-domain by sub-domain (Case 2,4) basis with 0%, 1%,..., 9% of the 44370 elements in the entire mesh being refined, as defined in Tab. 7.1. With the refinement confined to the two or four prism layers closest to the metallic antenna, a sub-domain is defined as the two or four adjacent prisms extruded from a given surface triangle, respectively. Since the input impedance of an antenna is an extremely sensitive parameter, it is expected to be influenced more than any other through improved field modeling. We therefore expect the different TVFE options to result in different input impedances. With the mixed-order TVFE of order 0.5 applied throughout the coarse mesh, the real 116

Error Base of adaptive Percentage of Constraint on Case indicator refinement refinement refinement 1 E I Tetra 0-9 2 middle layers 2 EI, Prisms (2 adjacent) 0-9 2 middle layers 3 Eh Tetra 0-9 4 middle layers 4 El1 Prisms (4 adjacent) 0-9 4 middle layers Table 7.1: Definition of Case 1-4. part of the input impedance is 77 Q while the corresponding imaginary part is -78 Q. The latter value is comparable in magnitude to the real part which is physically unrealistic for a broadband antenna like the LTSA considered here. With the mixedorder TVFE of order 0.5 applied throughout the dense mesh, the real part of the input impedance increases to 162 Q whereas the corresponding imaginary part of -81 Q is almost unchanged. However, if the mixed-order TVFE of order 0.5 is applied within the coarse mesh in conjunction with the mixed-order TVFE of order 1.5 in the four prism layers closest to the metallic antenna (40% of the elements), the real part of the input impedance increases to 234 Q while the corresponding imaginary part of 36 Q is almost an order of magnitude smaller. The real and imaginary parts of the input impedance for Case 1-4 are given in Fig. 7.9 as a function of the percentage of higher order TVFEs within the computational domain. The addition of just a few percent of higher order TVFEs is seen to have a dramatic influence on the input impedance. It converges very quickly to a real part of around 220 - 260 Q and an imaginary part that is significantly smaller. This is very close to the result when 40'% of the carvity is modeled with the mixed-order TVFE of order 1.5, i.e. a very small percentage of higher order TVFEs is needed provided these are placed properly. This is true for all the Cases 1-4, i.e. the specific refinement scheme is unimportant for this antenna. 117

E 200 15 //: rr - Case 2 50 -: - -: Case 3: Case 4 C L 0 1 2 3 4 5 6 7 8 9 Percentage of higher order TVFEs [%] 50 E 0 0 1 2 3 4 5 6 7 8 9 Percentage of higher order TVFEs [%] Figure 7.a: Real and imaginary part of the input impedance of the LTSA in Fig. v.8 for different refinement schemes. When mixed-order TVFEs of order 0.5 are applied throughout the coarse mesh, the system with 17823 unknowns converges in 2143 iterations (10-3 relative tolerance with a QMR solver). Improvement of accuracy via a denser mesh or (adaptive) refinement with mixed-order TVFEs of order 1.5 leads to increased memory and C'U time requirements. When mixed-order TVFEs of order 0.5 are applied throughout the dense mesh, the system with 164222 unknowns converges in 8858 iterations (10-3 relative tolerance with a QMR solver). When mixed-order TVFEs of order 0.5 are applied within the coarse mesh in conjunction with mixed-order TVFEs of order 1].5 in the four prism layers closest to the metallic antenna (40% of the elements), the system with 130967 unknowns converges in 7801 iterations (10-3 relative tolerance with a QMR solver). When mixed-order TVFEs of order 0.5 are applied within the coarse mesh in conjunction with mixed-order TVFEs of order 1.5 in regions found adaptively (Case 1-4), an initial solution with mixed-order TVFEs of order 0.5 must be computed (47823 unknowns and 2143 iterations for a 10i3 relative toler ance with 118

a QMR solver), regions of refinement must be determined and a refined solution must be computed. The number of iterations to reach this refined solution (10-3 relative tolerance with a QMR solver) with the previous solution used as the starting guess is given in Fig. 7.10 for Case 1-4 as a function of the percentage of higher order TVFEs within the computational domain. As for the input impedance, the number of iterations required for convergence reaches a plateau (6000-7000 iterations) after a few percent of higher order TVFEs have been added. That the number of iterations goes up as higher order TVFEs are added is fully expected and in agreement with previous observations in the literature and in this dissertation. In fact, the numbers of iterations when mixed-order TVFEs of order 1.5 are applied are by no means excessive since we are dealing with problems of 49155 (Case 1, 1% refinement) to 130967 (40% refinement) unknowns. Note that all approaches applying hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively lead to less unknowns and less iterations than when mixed-order TVFEs of order 0.5 are applied throughout a denser mesh for accuracy improvement. 8000 4 0 0 0 -..................................... 5000 CO 3 0 0 0 -....................... o 4000 1 3000 -....................... 1 I000:: ': ' i. Case 2 H u:-O- Case 3 0 1 2 3 4 5 6 7 8 9 Percentage of higher order TVFEs [%] Figure 7.10: Number of iterations to reach a relative tolerance of 10t with a QMR solver for the LTSA in Fig. 7.8 for different refinement schemes. 119

To illustrate the regions of refinement for Case 1-4, we restrict ourselves to 3%o refinement. We consider a cross-section of the mesh parallel to the antenna with the boundaries of the metallic cavity and antenna marked with thick black lines. As mentioned previously, the coarse tetrahedral volume mesh is grown from a coarse triangular surface mesh in a cut in the plane of the antenna by extruding it into ten prism layers and breaking each prism into three tetrahedra. With the refinement constrained to at most four prism layers around the antenna, anywhere between 0 and 12 tetrahedra can be refined corresponding to a given triangle in the cross-section. In the following, a white / light gray / average gray / dark gray triangle indicates that 0 i/ 1-4 / 5-8 / 9-12 tetrahedra are being refined. The regions of refinement for Case 1-4 (3% refinement) are given in Fig. 7.11-7.14. For Case 1 (Fig. 7.11), 0-6 elements can be refined and hence white, light gray or average gray is used. For Case 2 (Fig. 7.12), 0 or 6 elements can be refined and hence white or average gray is used. For Case 3 (Fig. 7.13), 0-12 elements can be refined and hence white, light gray, average gray or dark gray is used. For Case 4 (Fig. 7.14), 0 or 12 elements can be refined and hence white or dark gray is used. The figures show that the elements within and immediately bounding the tapered slot as well as elements close to the bottom metallic antenna edges in the vicinity of the feed are being refined. This is fully in agreement with the facts that the antenna works as a TWA and that certain fringing effects can be expected due to the sharp metallic edges. Electromagnetic near field parameters are generally more sensitive than far field parameters. A well-known example is electromagnetic scattering problems where fairly crude approximations to induced currents can often yield very accurate scat120

Figure 7.11: Regions of refinement for the LTSA in Fig. 7.8 for Case 1 with 3% refinement. Figure 7.12: Regions of refinement for refinement. the LTSA in Fig. 7.8 for Case 2 with -3% 121

Figure 7.13: Regions of refinement for the LTSA in Fig. 7.8 for Case 3 with 3% refinement. Figure 7.14: Regions of refinement for refinement. the LTSA in Fig. 7.8 for Case 4 with 3% 122

tered fields since integration of induced currents via radiation integrals is a very forgiving process. This also holds for certain antenna problems where crude approximations to equivalent near field currents can be integrated to provide very accurate far field patterns. To investigate the far field characteristics of the above LTSA as predicted with the different TVFE options, E- and H-plane patterns with mixedorder TVFEs of order 0.5 applied throughout the coarse mesh, mixed-order TVFEs of order 0.5 applied throughout the dense mesh and mixed-order TVFEs of order 0.5 applied within the coarse mesh in conjunction with mixed-order TVFEs of order 1.5 in the four prism layers closest to the antenna are given in Fig. 7.15-7.16 with the last pattern normalized to 0 dB at endfire. The need for accurate field modeling within the tapered slot is obvious from these figures as the pattern found with mixed-order TVFEs of order 0.5 applied throughout the coarse mesh differs significantly from the two others. The reason is that the poor field modeling within the tapered slot accumulates and leads to aperture fields so inaccurate that the integration of equivalent aperture currents cannot provide the correct far field patterns. More accurate field modeling via a denser mesh or addition of higher order TVFEs provides more accurate patterns. They are similar although the levels at endfire (0 = 0) and close to grazing are different. To investigate the merits of adaptive refinement and determine whether 40% higher order TVFEs are really needed, the E- and H-plane patterns with mixed-order TVFEs of order 0.5 applied within the coarse mesh in conjunction with mixed-order TVFEs of order 1.5 in the four prism layers closest to the antenna are repeated in Fig. 7.17-7.18 where also E- and H-plane patterns for Case 1-4 with 3% refinement 123

E-plane -80 -60 -40 -20 0 20 40 60 80 0 [degrees] Figure 7.15: E-plane patterns for the LTSA in Fig. 7.8. H-plane -80 -60 -40 -20 0 20 40 60 80 0 [degrees] Figure 7.16: H-plane patterns for the LTSA in Fig. 7.8. 124

E-plane -80 -60 -40 -20 0 e [degrees] 20 40 60 80 Figure 7.17: E-plane patterns for the LTSA in Fig. 7.8. H-plane -80 -60 -40 -20 0 20 40 60 80 0 [degrees] Figure 7.18: H-plane patterns for the LTSA in Fig. 7.8. 125

are given (same normalization as above). We observe that Case 1-4 give almost identical patterns. Discrepancies around the E-plane shoulders at 40 - 60~ can be viewed but overall they are of similar value and shape, i.e. the specific refinement scheme is unimportant for this particular application. More importantly, the patterns for Case 1-4 with 3% refinement are seen to be very similar to that corresponding to 40% higher order TVFIEs, expressing again a need for very few higher order TVFEs for accurate field modeling provided these are placed properly. 7.3 Summary In this chapter, the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra proposed in this dissertation were used in conjunction with a simple adaptive refinement strategy to analyze the impedance and pattern characteristics of TSAs using the hybrid FE/BI method. The adaptive inclusion of a very small percentage of higher orderor TVFEs was found to have a dramatic effect on the accuracy of the computed input impedances and far field patterns, thus justifying the approach proposed in this dissertation for large and complex problems. 126

CHAPTER 8 SUMMARY, CONCLUSIONS AND FUTURE WORK In this chapter, brief summaries and the most important conclusions for the individual chapters are given and several future tasks to be completed are suggested. 8.1 Summary and conclusions In chapter 1, the work presented in this dissertation was introduced. After a brief motivation, some fundamental concepts were presented, a high-level description of the proposed approach was given and the organization of the dissertation was outlined. In chapter 2, background material was given. Vector wave equations used throughout the dissertation were presented and tangential vector finite elements (TVFEs) for triangular and tetrahedral elements used for discretizing partial differential equations were reviewed. In chapter 3, hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for tri 127

angular elements were proposed. An efficient set of vector basis functions for the expansion of the surface current on a perfectly electrically conducting (PEC) generalized quadrilateral was converted to vector basis functions applicable for finite element (FE) analysis and hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 were proposed for triangular elements. The proposed hierarchical mixed-order TVFEs of order 0.5 and 1.5 were tested for solution of closed- as well as open-domain problems. For solution of certain classes of electromagnetic problems, field expansion using hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively was found to be a very promising approach in terms of accuracy, memory and central processing unit (CPU) time requirements as compared to a more traditional approach. In chapter 4, hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for tetrahedral elements were proposed. They were constructed as the direct three-dimensional equivalents of the hierarchical mixed-order TVFEs of order 0.5, 1.5 and 2.5 for triangular elements proposed in chapter 3. The proposed hierarchical mixed-order TVIEs of order 0.5 and 1.5 were tested for solution of closed- as well as open-domain problems. Again, selective field expansion was found to be a very promising approach for accurate and efficient solution of certain classes of electromagnetic problems. In chapter 5, the condition numbers resulting from finite element method (FEM) analysis using the hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements proposed in chapter 3 and chapter 4 were contrasted to those of existing interpolatory and hierarchical mixed-order TVFEs of order 1.5 for triangular and tetrahedral elements. The proposed hierarchical mixed-order TVFEs of order -1.5 proved better conditioned than existing hierarchical mixed-order TVFEs of order 1.5 128

and thus the analysis fostered no concerns for potential future convergence problems due to excessive matrix condition numbers. In addition, an approach for improving the condition numbers of FEM matrices resulting from selective field expansion was suggested and tested. The improvement comes at the expense of a more complicated formulation and computer code but does not alter accuracy. In chapter 6, a review of existing error estimators and indicators was given and the effectiveness of the proposed hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra was investigated when some of the reviewed error indicators were applied in the context of a, very simple adaptive refinement strategy. The results were extremely promising for both narrowband and broadband antennas provided the refinement was carried out on a sub-domain by sub-domain basis as opposed to an element by element basis. In chapter 7, the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedra proposed in this dissertation were used in conjunction with a simple adaptive refinement strategy to analyze the impedance and pattern characteristics of tapered slot antennas (TSAs) using the hybrid finite element / boundary integral (FE/BI) method. The adaptive inclusion of a very small percentage of higher order TVFEs was found to have a dramatic effect on the accuracy of the computed input impedances and far field patterns, thus justifying the approach proposed in this dissertation for large and complex problems. 129

8.2 Future work The hybrid F'E/BI formulation for three-dimensional open-domain problems presented in this dissertation is based on a lowest order boundary integral (BI) formulation even when higher order TVFEs are (partly or fully) used within the interior of the computational domain. This implies that even when higher order TVFEs are (partly or fully) used within the cavity, higher order vector basis functions associated with edges and faces on the BI surface are eliminated to maintain field uniqueness. Since the BI surface is partly in the near field of the metallic antenna, a (partly or fully) higher order BI formulation is expected to provide superior accuracy. Such a formulation would require a (partly or fully) higher order testing scheme for the pertinent integral equation. This formulation could be carried out and the results could be compared to those based on the lowest order BI formulation. Motivated by the superiority of hierarchical mixed-order TVFEs of order 1.5 over those of order 0.5, the merits of even higher order hierarchical mixed-order TVFEs could be investigated. Those of order 2.5 for triangular and tetrahedral elements proposed in this dissertation could be implemented and their effectiveness could be assessed. If deemed necessary, even higher order TVFEs could be suggested, implemented and evaluated following the principles outlined in this dissertation. Although methods for adaptive refinement have been studied in mathematics and engineering for decades, use of such methods for solution of practical engineering problems is only beginning to emerge. This dissertation includes a study of adaptive TVFE refinement for the hierarchical mixed-order TVFEs proposed in this disserta 130

tion using three different error indicators in the context of a very simple adaptive refinement strategy. Alternative hierarchical mixed-order TVFEs, other error indicators, various error estimators and more complex adaptive refinement strategies could be studied to provide an understanding of the inter-relations between the different approaches and their applicability for solution of practical electromagnetic problems. This dissertation is limited solely to polynomial order refinement (p-refinement) techniques with a uniform mesh density. The effectiveness of mesh density refinement (Ih-refinement) techniques with a uniform polynomial order could be studied for similar problems and contrasted to that offered by p-refinement techniques with a uniform mesh density. Subsequently, hp-refinement techniques could be developed. Dual hp-refinement techniques are theoretically always superior to isolated h-refinement or p-refinement techniques but often difficult to implement practically [29]. Useful progress to this end has been reported in [22, 89]. As justified previously, this dissertation focuses on triangular and tetrahedral elements due to their geometrical modeling flexibility. However, hierarchical mixed-order TVFEs for alternative element shapes could be developed and contrasted to those proposed in this dissertation. Hierarchical mixed-order TVFEs for curved triangular and tetrahedral elements could be constructed from those proposed in this dissertation for straight triangular and tetrahedral elements via a straightforward mapping, see for instance [37], and similar hierarchical mixed-order TVFEs for other elements (rectangles, bricks, prisms, pyramids) could be constructed by direct analogy. A comparative study of all these hierarchical mixed-order TVFEs could be carried out for solution of practical electromagnetic circuit, scattering or radiation problems. 131

APPENDICES 132

Appendix A Explicit expressions for Wi, W2 and W3 To derive an expression for W1, we introduce two coordinates (u1, v1) over the triangle. These are degenerates of similar coordinates for a generalized quadrilateral [69]. ul takes its minimum value ul,min 0 at node 1 and its maximum value Ul,max = 1 along edge #2 while v1 takes its minimum value vl,i — 1 along edge #3 and its maximum value Vrmax = 1 along edge #1. u1 is constant and v, is linear along straight lines parallel to edge #2 while u1 is linear and v1 is constant along straight lines starting a node 1 and ending at edge #2, as illustrated in Fig. A.I. Using these coordinates, the position vector r defining P can be expressed as [69] r r= + ul ru1 + uivl ruil (A.1) where r, - [(r3- ri) + (r2- ri)] (A.2) rul = (r2-r3). (A.3) Further, u1 and v1 can be shown to be related to the simplex coordinates (1, (2 and 133

Figure A.1: Illustration of the variation of u1 and v1 over a triangle. 43 via U1 = (2 + ~3 (A.4) (2 - C3 V1 -- (A.5) C2 + (3 From (3.4) for n = 1, trivial algebra then leads to W1 = C2V(3 - C3VC2. (A.6) To derive expressions for W2 and W3, we can similarly introduce coordinates (u2, v2) and (U3, V3) where U2,3 = 0 at node 2, 3, U2,3 = 1 at the edge opposite to node 2, 3 and V2,3 = ~1 along the two edges shared by node 2,3. The algebra is similar and we arrive at U2 = C3 + C1 (A.7) V2 3 - A.8) 3 + (A.S) 13 + 4I W2 = C3VCi - C1VC3 (A.9) 134

U3 (=,i + (2 V3 _ _ _ _ 2 w3 -I ( (A.1O ) (A.1I1) (A.[-12) 135

Appendix B Expressions for vector basis functions Peterson's interpolatory mixed-order TVFE of order 1.5 Peterson's interpolatory mixed-order TVFE of order 1.5 is characterized by the eight vector basis functions W: = (2V(3 (B.1) WI = C3VC2 (B.2) W1 = C3V(C (B.3) W = CVC3 (B.4) W = lVC2 (B.5) W1 = C2VCl (B.6) W = C3(C1VC2 - (2VCi) (B.7) W = Cl (2VC33 - C3V2). (B.8) 136

Proposed hierarchical mixed-order TVFE of order 1.5 The proposed hierarchical mixed-order TVFE of or der 1.5 is characterized by the eight vector basis functions W1- = 2VC3 - C3V(2 w2 = C3V(l - (iV(3 (B.9) (B.10) w2- - 3 - 1V2 - (2V(l (B. l 1) W4 = ((2 - -3)(C2V(3 -,3V(2) W2 = ((3 - l1)(3V(l - CiV(3) (B.12) (B. 13) W2 - ( - 2)(iV(2 - 2V(l) (B.1L4) W7' -- (3 ((I V (2 - (2 V (1) (B.15) W' = (l((2V(3 - (3V(2) (B.] 6) Transformation between mixed-order TVFEs of order 1.5 The two mixed-order TVFEs of order 1.5 presented above are related through a linear transformation. Let [W1] be a column vector containing Peterson's eight vector basis functions WJ, j = 1,.., 8, and [W2] be a column vector containing the 137

proposed eight vector basis functions W?, i 1,., 8. In this case, [W2] is related to [WI] via jVLIj II [W2] =[Aj][Wl] (B.17) where [A i] is the sparse 8 x 8 transformation matrix [A.j] = 1 0 0 1 0 0 0 0 -1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 -1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 C C C 0 1 0 0 O ) 0 0 1 0 ) -2 0 0 0 -1 -1 2 0 1 -I (B. L8) 1 1 1 0 Proposed hierarchical mixed-order TVFE of order 2.5 The proposed hierarchical mixed-order TVFE of order 2.5 is characterized by the fifteen vector basis functions W3I = C2V(3 - C3VC2 w3 = (3VC, - Cl V3 W3 = GCV~2 -- 2V(C (B.19) (B.20) (B.2 ) w3= (2 - C3)(82V(3 - (3V(2) 138 (B.22)

W3 - ((3 _ (1)((3V(l _ 5 - O V(3)(B2) (B. 2, J) - ((I ((2)(V(2 - w7 ((2 ((2V( -3 9NV12 - (2 (( V (2 - (B.24) (3V(2) (B.25) (B. 2)(~ - (2V(l) (B. 27) (B.28) (B.29) 2V6) (B.3 0) w13 ( 1((2V7(3 - (3V(2) w3 -(2((V I- (I V(3) (B.3-[) (B.32) (B.3 33) 139

BIBLIOGRAPHY 140

BIBLIOGRAPHY [1] M. Ainsworth and J.T. Oden. 'A procedure for a posteriori error estimation for h-p finite element methods'. Computer Methods in Applied Mechanics and Engineering, vol. 101, pp. 73-96, December 1992. [2] M. Ainsworth and J.T. Oden. 'A posteriori error estimation in finite element analysis'. Computer Methods in Applied Mechanics and Engineering, vol. 142, pp. 1-88, March 1997. [3] L.S. Andersen. Scattering by non-perfectly conducting structures. MSc thesis, Technical University of Denmark, 1995. [4] L.S. Andersen and J.L. Volakis. 'Adaptive multi-resolution antenna modeling using hierarchical mixed-order tangential vector finite elements'. Submitted to IEEE Transactions on Antennas and Propagation. [5] L.S. Andersen and J.L. Volakis. 'Accurate and efficient simulation of antennas using hierarchical mixed-order tangential vector finite elements for tetrahedra'. IEEE Transactions on Antennas and Propagation, vol. AP-47, pp. 1240-1243, August 1999. [6] L.S. Andersen and J.L. Volakis. 'Condition numbers for various FEM matrices'. Journal of Electromagnetic Waves and Applications, vol. 13, pp. 1661-1677, December 1999. [7] L.S. Andersen and J.L. Volakis. 'Mixed-order tangential vector finite elements for triangular elements'. IEEE Antennas and Propagation Magazine, vol. APAM40, pp. 104-108, February 1998. [8] L.S. Andersen and J.L. Volakis. 'Development and application of a novel class of hierarchical tangential vector finite elements for electromagnetics'. IEEE Transactions on Antennas and Propagation, vol. AP-47, pp. 112-120, January 1999. [9] L.S. Andersen and J.L. Volakis. 'Hierarchical tangential vector finite elements for tetrahedra'. IEEE Microwave and Guided Wave Letters, vol. 8, pp. 127-129, March 1998. 141

[10] I. Babuska, L. Planck, and R. Rodriguez. 'Basic problems of a posteriori error estimation'. Computer Methods in Applied Mechanics and Engineering, vol. 101, pp. 97-112, December 1992. [11] S.S. Bindiganavale. Fast memory-saving hybrid algorithms for electromagnetic scattering and radiation. PhD thesis, University of Michigan, 1997. [12] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz. 'AIM: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems'. Radio Science, vol. 31, pp. 1225-1251, September-October 1996. [13] A. Bossavit. 'Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism'. Proceedings of the IEE, Part H, vol. 135, pp. 493-500, February 1988. [14] C. Carrie and J.P. Webb. 'Hierarchal triangular edge elements using orthogonal polynomials'. In Proc. of the IEEE Antennas and Propagation Society Interinational Symposium 1997, Montreal, Quebec, Canada, vol. 2, pp. 1310-1313, July 1997. [15] Z.J. Cendes. 'Vector finite elements for electromagnetic field computation'. IEEE Transactions on Magnetics, vol. MAG-27, pp. 3958-3966, September 1991. [16] Z.J. Cendes and D.N. Shenton. 'Adaptive mesh refinement in the finite element computation of magnetic fields'. IEEE Transactions on Magnetics, vol. MAG21, pp. 1811-1816, September 1985. [17] J.S. Colburn and Y. Rahmat-Samii. 'Linear taper slot antenna directivity improvement via substrate perforation: A FDTD evaluation'. In Proc. of the IEEE Antennas and Propagation Society International Symposium 1998, Atlanta, Georgia, USA, vol. 2, pp. 1176-1179, June 1998. [18] R.E. Collin. Antennas and radiowave propagation. McGraw-Hill Book Company, USA, 1985. [19] M.E. Cooley. D.H. Schaubert, N.E. Buris, and E.A. Urbanik. 'Radiation and scattering analysis of infinite arrays of endfire slot antennas with a ground plane'. IEEE Transactions on Antennas and Propagation, vol. AP-39, pp. 1615 -1624, November 1991. [20] W. Daigang and J. Kexun. 'P-version adaptive computation of FEM'. IEEE Transactions on Magnetics, vol. MAG-30, pp. 3515-3518, September 1994. [21] I. Daubechies. Ten lectures on wavelets. Society for Industrial and Applied Mathematics, USA, 1992. 142

[22] L. Demkowicz and L. Vardapetyan. 'Modeling of electromagnetic absorption/scattering problems using lip-adaptive finite elements'. Computer Methods in Applied Mllechanics and Engineering, vol. 152, pp. 103-124, January 1998. [23] R.E. Ewing. 'A posteriori error estimation'. Computer Methods in Applied Mechanics and Engineering, vol. 82, pp. 59-72, September 1990. [24] P. Fernandes, P. Girdinio, P. Molfino, and M. Repetto. 'Local error estimates for adaptive mesh refinement'. IEEE Transactions on Magnetics, vol. MAG-24, pp. 299-302, January 1988. [25] D.S. Filipovic, L.S. Andersen, and J.L. Volakis. 'A multi-resolution method for simulating infinite periodic arrays'. Submitted to IEEE Transactions on Antennas and Propagation. [26] N. Fourikis, N. Lioutas, and N.V. Shuley. 'Parametric study of the co- and crosspolarisation characteristics of tapered planar and antipodal slotline antennas'. Proceedings of the IEE, Part H, vol. 140, pp. 17-22, February 1993. [27] L.E. Garcia-Castillo and M. Salazar-Palma. 'Second order Nedelec tetrahedral element for computational electromagnetics'. In Proc. of the 4th International Workshop on Finite Elements for Microwave Engineering 1998, FuturoscopePoitiers, France, July 1998. [28] E. Gazit. 'Improved design of the Vivaldi antenna'. Proceedings of the IEE, Part H, vol. 135, pp. 89-92, April 1988. [29] D. Giannacopoulos and S. McFee. 'Towards optimal h-p adaption near singularities in finite element electromagnetics'. IEEE Transactions on Magnetics, vol. MAG-30, pp. 3523-3526, September 1994. [30] P.J. Gibson. 'The Vivaldi aerial'. Proc. 9th European Microwave Conference, Brighton, UK, pp. 101-105, 1979. [31] P. Girdinio, P. Molfino, G. Molinari, L. Puglisi, and A. Viviani. 'Finite difference and finite element grid optimization by the grid iteration method'. IEEE Transactions on Magnetics, vol. MAG-19, pp. 2543-2546, November 1983. [32] N.A. Golias, A.G. Papagiannakis, and T.D. Tsiboukis. 'Efficient mode analysis with edge elements and 3-D adaptive refinement'. IEEE Transactions on Microwave Theory and Techniques, vol. MTT-42, pp. 99-107, January 1994. [33] N.A. Golias and T.D. Tsiboukis. 'Adaptive methods in computational magnetics'. International Journal of Numerical Modeling: Electronic Networks, Devices and Fields, vol. 9, pp. 71-80, January-April 1996. [34] N.A. Golias and T.D. Tsiboukis. 'Adaptive refinement in 2-D finite element applications'. International Journal of Numerical Modeling: Electronic Networks, Devices and Fields, vol. 4, pp. 81-95, June 1991. 143

[35] N.A. Golias and T.D. Tsiboukis. 'Adaptive refinement strategies in three dimensions'. IEEE Transactions on Magnetics, vol. MAG-29, pp. 1886-1889, March 1993. [36] G.H. Golub and C.F. Van Loan. Matrix computations. The Johns Hopkins University Press, USA, 1983. [37] R. Graglia, D.R. Wilton, and A.F. Peterson. 'Higher order interpolatory vector bases for computational electromagnetics'. IEEE Transactions on Antennas and Propagation, vol. AP-45, pp. 329-342, March 1997. [38] A.K. Gupta. 'A finite element for transition from a fine to a coarse grid'. International Journal for Numerical Methods in Engineering, vol. 12, pp. 35 -45, 1978. [39] S. Hahn, C. Calmels, G. Meunier, and J.L. Coulomb. 'A posteriori error estimate for adaptive finite element mesh generation'. IEEE Transactions on AMlagnetics, vol. MAG-24, pp. 315-317, January 1998. [40] R.F. Harrington. Time harmonic electromagnetic fields. McGraw-Hill Book Company, USA, 1961. [41] S.R.H. Hoole, S. Jayakumaran, A.W. Ananadaraj, and P.R.P. Hoole. 'Relevant, purpose based error criteria for adaptive finite element mesh generation'. Journal of Electromagnetic Waves and Applications, vol. 3, pp. 167-177, February 1989. [42] S.R.H. Hoole, S. Jayakumaran, and N.R.G. Hoole. 'Flux density and energy perturbation in adaptive finite element mesh generation'. IEEE Transactions on Magnetics, vol. MAG-24, pp. 322-325, January 1988. [43] R. Janaswamy. 'An accurate moment method model for the tapered slot antenna'. IEEE Transactions on Antennas and Propagation, vol. AP-37, pp. 1523 -1528, December 1989. [44] R. Janaswamy and D.H. Schaubert. 'Analysis of the tapered slot antenna'. IEEE Transactions on Antennas and Propagation, vol. AP-35, pp. 1058-1064, September 1987. [45] J.M. Jin and J.L. Volakis. 'A finite element - boundary integral formulation for scattering by three-dimensional cavity-backed apertures'. IEEE Transactions on Antennas and Propagation, vol. AP-39, pp. 97-104, January 1991. [46] C. Johnson and P.Johnson and P. Hansbo. 'Adaptive finite element methods in computational mechanics'. Computer Methods in Applied Mechanics and Engineering, vol. 101, pp. 143-181, December 1992. 144

[47] D.W. Kelly, J.P.S.R. Gago, O.C. Zienkiewicz, and I. Babuska. 'A posteriori error analysis and adaptive processes in the finite element method: Part I - Error analysis'. International Journal for Numerical Methods in Engineering, vol. 19, pp. 1593-1619, November 1983. [48] H. Kim, S. Hong, K. Choi, H. Jung, and S. Hahn. 'A three dimensional adaptive finite element method for magnetostatic problems'. IEEE Transactions on Mlagnetics, vol. MAG-27, pp. 4081-4084, September 1991. [49] A. Koksal and J.F. Kauffman. 'Moment method analysis of linearly tapered slot antennas'. International Journal of Microwave and Millimeter- Wave ClomputerAided Engineering, vol. 4, pp. 76-87, January 1994. [50] B.M. Kolundzija and B.D. Popovic. 'Entire-domain galerkin method for analysis of metallic antennas and scatterers'. Proceedings of the IEE, Part H, vol. 140, pp. 1-10, February 1993. [51] P.S. Kooi, T.S. Yeo, and M.S. Leong. 'Parametric studies of the linearly tapered slot antenna (LTSA)'. Microwave and Optical Technology Letters, vol. 4, pp. 200-207, April 1991. [52] J.D.S. Langley, P.S. Hall, and P. Newham. 'Novel ultrawide-bandwidth Vivaldi antenna with low crosspolarisation'. Electronics Letters, vol. 29, pp. 2004-2005, 11th November 1993. [53] J.F. Lee, D.K. Sun, and Z.J. Cendes. 'Full-wave analysis of dielectric waveguides using tangential vector finite elements'. IEEE Transactions on Microwave Theory and Techniques, vol. MTT-39, pp. 1262-1271, August 1991. [54] J.F. Lee, D.K. Sun, and Z.J. Cendes. 'Tangential vector finite elements for electromagnetic field computation'. IEEE Transactions on Magnetics, vol. MAG27, pp. 4032-4035, September 1991. [55] A.C. Ludwig. 'The definition of cross polarization'. IEEE Transactions on Antennas and Propagation, vol. AP-31, pp. 116-119, January 1973. [56] S.G. Mallat. 'A theory for multiresolution signal decomposition: The wavelet representation'. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 11, pp. 674-693, July 1989. [57] F.J.C. Meyer and D.B. Davidson. 'Error estimates and adaptive procedures for the two-dimensional finite element method'. Electronics Letters. vol. 30, pp. 936-938, 9th June 1994. [58] J.B. Muldavin, T.J. Ellis, and G.M. Rebeiz. 'Tapered slot antennas on thick dielectric substrates using micromachining techniques'. In Proc. of the IEEE Antennas and Propagation Society International Symposium 1997, MIontreal, Quebec, Canada, vol. 2, pp. 1110-1113, July 1997. 145

[59] G. Mur. 'The fallacy of edge elements'. IEEE Transactions on Magnetics, vol. MAG-34, pp. 3244-3247, September 1998. [60] G. Mur and A.T. de Hoop. 'A finite-element method for computing threedimensional electromagnetic fields in inhomogeneous media'. IEEE Transactions on Magnetics, vol. MAG-21, pp. 2188-2191, November 1985. [61] J.C. Nedelec. 'Mixed finite elements in R3,. Num. Math., vol. 35, pp. 315-341, 1980. [62] J.C. Nedelec. 'A new family of mixed finite elements in R3'. Num. Math., vol. 50, pp. 57-81, 1986. [63] J.T. Oden, L. Demkowicz, W. Rachowicz, and T.A. Westermann. 'Toward a universal h-p adaptive finite element strategy, Part 2. A posteriori error estimation'. Computer Methods in Applied Mechanics and Engineering, vol. 77, pp. 113-180, December 1989. [64] T. Ozdemir. Finite element analysis of conformal antennas. PhD thesis, University of Michigan, 1998. [65] A.F. Peterson. Private communication, e-mail, August 1997. [66] A.F. Peterson. 'Vector finite element formulation for scattering from twodimensional heterogeneous bodies'. IEEE Transactions on Antennas and Propagation, vol. AP-42, pp. 357-365, March 1994. [67] A.F. Peterson and D.R. Wilton. 'Curl-conforming mixed-order edge elements for discretizing the 2D and 3D vector Helmholtz equation'. In T. Itoh, G. Pelosi, and P.P. Silvester, editors, Finite element software for microwave engineering, chapter 5, pp. 101-124. John Wiley and Sons, Inc., USA, 1996. [68] A.F. Peterson and D.R. Wilton. 'A rationale for the use of mixed-order basis functions within finite element solutions of the vector Helmholtz equation'. In Proc. of the 11th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, USA, vol. 2, pp. 1077-1084, March 1995. [69] B.D. Popovic and B.M. Kolundzija. Analysis of metallic antennas and scatterers. IEE Electromagnetic Waves Series, vol. 38, 1994. [70] S.M. Rao, ).R. Wilton, and A.W. Glisson. 'Electromagnetic scattering by surfaces of arbitrary shape'. IEEE Transactions on Antennas and Propagation, vol. AP-30, pp. 409-418, May 1982. [71] V. Rocklin. 'Rapid solution of integral equations for scattering theory in two dimensions'. Journal of Computational Physics, vol. 86, pp. 414-439, 1990. [72] Y. Saad. Iterative methods for sparse linear systems. PWS Publishing Company, USA, 1996. 146

[73] K. Sabetfakhri. Novel efficient integral-based techniques for characterization of planar microwave structures. PhD thesis, University of Michigan, 1995. [74] M. Salazar-Palma, T.K. Sarkar, L.E. Garcia-Castillo, T. Roy, and A. I)jordjevic. Iterative and self-adaptive finite-elements in electromagnetic modeling. Artech House, Inc., USA, 1998. [75] J.S. Savage. 'Comparing high order vector basis functions'. In Proc. of the 14th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, USA, vol. 2, pp. 742-749, March 1998. [76] J.S. Savage and J. Manges. 'Local error estimation for high-frequency problems using hierarchical tangential vector finite elements'. In Proc. of the 15th Annual Review of Progress in Applied Conmputational Electromagnetics, Monterey, CA, USA, pp. 524-529, March 1999. [77] J.S. Savage and A.F. Peterson. 'Higher-order vector finite elements for tetrahedral cells'. IEEE Transactions on Microwave Theory and Techniques, vol. MTT-44, pp. 874-879, June 1996. [78] D.H. Schaubert. 'A class of E-plane scan blindnesses in single-polarized arrays of tapered-slot antennas with a ground plane'. IEEE Transactions on Antennas and Propagation, vol. AP-44, pp. 954-959, July 1996. [79] D.H. Schaubert, J.A. Aas, M.E. Cooley, and N.E. Buris. 'Moment method analysis of infinite stripline-fed tapered slot antenna arrays with a ground plane'. IEEE Transactions on Antennas and Propagation, vol. AP-42, pp. 1161-1166, August 1994. [80] D.H. Schaubert, J. Shin, and G. Wunsch. 'Characteristics of single-polarized phased array of tapered slot antennas'. In Proc. of the IEEE International Symposium on Phased Array Systems and Technology 1996, Boston, MA, USA, pp. 102-106, October 1996. [81] J.W. Schuster and R. Luebbers. Private communication. [82] M.S. Shephard. 'Automatic and adaptive mesh generation'. IEEE Transactions on Malagnetics, vol. MAG-21, pp. 2484-2489, November 1985. [83] J. Shin and I).H. Schaubert. 'A parameter study of stripline-fed Vivaldi notchantenna arrays'. IEEE Transactions on Antennas and Propagation, vol. AP-47, pp. 879-886, May 1999. [84] B.Z. Steinberg and Y. Leviatan. 'On the use of wavelet expansions in the method of moments'. IEEE Transactions on Antennas and Propagation, vol. AP-41, pp. 610-619, May 1993. 147

[85] D. Sun, J. Manges, X. Yuan, and Z. Cendes. 'Spurious modes in finite-element methods'. IEEE Antennas and Propagation Magazine, vol. APM-37, pp. 12-24, October 1995. [86] B.A. Szabo. 'Estimation and control of error based on p convergence'. In I. Babuska, O.C. Zienkiewicz, J. Gago, and E.R.A. Oliveira, editors, Accuracy estimates and adaptive refinements in finite element computations. John Wiley and Sons, Inc., USA, 1986. [87] R.L. Taylor. 'On completeness of shape functions for finite element analysis'. International Journal for Numerical Methods in Engineering, vol. 4, pp. 12-22, 1972. [88] E. Thiele and A. Taflove. 'FD-TD analysis of Vivaldi flared horn antennas and arrays'. IEEE Transactions on Antennas and Propagation, vol. AP-42, PP. 633-641, May 1994. [89] L. Vardapetyan and L. Demkowicz. 'hp-adaptive finite elements in electromagnetics'. Computer Methods in Applied Mechanics and Engineering, vol. 169, pp. 331-344, February 1999. [90] M. Vetterli and J. Kovacevic. Wavelets and subband coding. Prentice-Hall, USA, 1995. [91] J.L. Volakis, A. Chatterjee, and L.C. Kempel. Finite element method for electromagnetics. IEEE Press, USA, 1998. [92] J.L. Volakis, T. Ozdemir, and J. Gong. 'Hybrid finite element methodologies for antennas and scattering'. IEEE Transactions on Antennas and Propagation, vol. AP-45, pp. 493-507, March 1997. [93] J. Wang and J.P. Webb. 'Hierarchal vector boundary elements and p-adaption for 3-D electromagnetic scattering'. IEEE Transactions on Antennas and Propagation, vol. AP-45, pp. 1869-1879, December 1997. [94] J.S. Wang. On "edge" based finite elements and methods of moments solutions of electromagnetic scattering and coupling. PhD thesis, University of Akron, 1992. [95] J.P. Webb. 'Hierarchal vector basis functions of arbitrary order for triangular and tetrahedral finite elements'. IEEE Transactions on Antennas and Propagation, vol. AP-47, pp. 1244-1253, August 1999. [96] J.P. Webb. 'Edge elements and what they can do for you'. IEEE 'Transactions on Magnetics, vol. MAG-29, pp. 1460-1465, March 1993. [97] J.P. Webb and B. Forghani. 'Hierarchal scalar and vector tetrahedra'. IEEE Transactions on Magnetics, vol. MAG-29, pp. 1495-1498, March 1993. 148

[98] H. Whitney. Geometric integration theory. Princeton University Press, USA, 1957. [99] M.V. Wickerhauser. Adapted wavelet analysis from theory to software. IEEE Press, USA, 1994. [100] J.Y. Wu and R. Lee. 'Construction of the basis functions for Nedelec's finite element spaces for triangular and tetrahedral elements'. In Proc. of the North Anmerican Radio Science AMeeting 1997, Montreal, Quebec, Canada, pp. 38, July 1997. [101] G.J. Wunsch and D.H. Schaubert. 'Effects on scan blindness of full and partial crosswalls between notch antenna array unit cells'. In Proc. of the IESEE Antennas and Propagation Society International Symposium 1995, Newport Beach, California, USA, vol. 4, pp. 1818-1821, June 1995. [102] Z. Xiang and Y. Lu. 'An effective wavelet matrix transformation approach for efficient solutions of electromagnetic integral equations'. IEEE Transactions on Antennas and Propagation, vol. AP-45, pp. 1205-1213, August 1997. [103] T.V. Yioultsis and T.D. Tsiboukis. 'Development and implementation of second and third order vector finite elements in various 3-D electromagnetic field problems'. IEEE Transactions on Magnetics, vol. MAG-33, pp. 1812-1815, March 1997. [104] T.V. Yioultsis and T.D. Tsiboukis. 'Multiparametric vector finite elements: A systematic approach to the construction of 3-D higher order tangential vector shape functions'. IEEE Transactions on Magnetics, vol. MAG-32, pp. 1389 -1392, May 1996. [105] T.V. Yioultsis and T.D. Tsiboukis. 'A generalized theory of higher order vector finite elements and its applications in three-dimensional electromagnetic scattering problems'. Electromagnetics, vol. 18, pp. 467-480, September-October 1998. [106] K.S. Yngvesson, T.L. Korzeniowski, Y.-S. Kim, E.L. Kollberg, and J.F. Johansson. 'The tapered slot antenna - a new integrated element for millimeter-wave applications'. IEEE Transactions on Microwave Theory and Techniques, vol. MTT-37, pp. 365-374, February 1989. [107] K.S. Yngvesson, D.H. Schaubert, T.L. Korzeniowski, E.L. Kollberg, T. Thungren, and J.F. Johansson. 'Endfire tapered slot antennas on dielectric substrates'. IEEE Transactions on Antennas and Propagation, vol. AP-33, pp. 1392-1400, Iecember 1985. [108] O.C. Zienkiewicz and R.L. Taylor. The finite element method - Volumne 1 - Basic formulation and linear problems. McGraw-Hill Book Company (UK) Ltd., UK, 1989. 149

[109] O.C. Zienkiewicz and J.Z. Zhu. 'A simple error estimator and adaptive procedure for practical engineering analysis'. International Journal for Numerical Methods in, Engineering, vol. 24, pp. 337-357, February 1987. [110] F.J. Zucker. 'Surface-wave antennas and surface-wave-excited arrays'. In R.C. Johnson and H. Jasik, editors, Antenna engineering handbook, chapter 12, pp. 12.1-12.36. McGraw-Hill Book Company, USA, 1984. 150