FINAL REPORT 028918-5-T A hybrid finite element-boundary integral for the analysis of cavity-backed antennas of arbitrary shape Jian Gong, John L. Volakis, A. C. Woo and H. T. G. Wang National Aeronautics and Space Administration Ames Research Center Moffett Field CA 94035 Naval Weapons Center China Lake CA 93555-6001 August 1993 28918-5-T = RL-2401

Report #028918-5-T NASA Interchange No. NCA2-653 FINAL REPORT Grant Interchange Title: Characterization of spiral patch antennas Report Title: A hybrid finite element-boundary integral for the analysis of cavity-backed antennas of arbitrary shape J. Gong, J. L. Volakis, A. C. Woo and H. T. G. Wang Report Authors: Primary University Collaborator: John L. Volakis Volakis@um.cc.umich.edu Telephone: (313) 764-0500 Alex Woo woo@ra-next.arc.nasa.gov Telephone: (415) 604-6010 Primary NASA-Ames Collaborator: University Address: Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor MI 48109-2122 Date: August 1993 Funds for the support of this study have been allocated in part by the NASA-Ames Research Center, Moffett Field, California, under interchange No. NCA2-653.

Preface This is the final report on this project which was concerned with the analysis of cavity-backed antennas and more specifically spiral antennas. The project was a continuation of a previous analysis, which employed rectangular brick elements, and was, thus, restricted to planar rectangular patch antennas. A total of five reports were submitted under this project and we expect that at least four journal papers will result from the research described in these reports. The abstracts of the four previous reports follow this page. The first of the reports (028918-1-T) is over 75 pages and describes the general formulation using tetrahedral elements and the computer program. Report 028918-2-T was written after the completion of the computer program and reviews the capability of the analysis and associated software for planar circular rectangular patches and for a rectangular planar spiral. Measurements were also done at the University of Michigan and at Mission Research Corp. for the purpose of validating the software. We are pleased to acknowledge a partial support from Mission Research Corp. in carrying out the work described in this report. The third report (028918-3-T) describes the formulation and partial validation (using 2D data) for patch antennas on a circular platform. The 3D validation and development of the formulation for patch antennas on circular platforms is still in progress. The fourth report (028918-4-T) is basically an invited journal paper which will appear in the J. Electromagnetic Waves and Applications in early 1994. It describes the application of the finite element method in electromagnetics and is primarily based on our work here at U-M. This final report describes the culmination of our efforts in characterizing complex cavity-backed antennas on planar platforms. The report describes for the first time the analysis of non-planar spirals and non-rectangular slot antennas as well as traditional planar patch antennas. The comparisons between measurements and calculations are truly impressive. Another unique aspect of this work is the incorporation of the FFT as part of the BiCG solver by overlaying a structured triangular mesh over the unstructured mesh. The implementation of this BiCG-FFT solution algorithm is important in minimizing the CPU and storage requirements. This final report will be submitted for publication in a refereed journal. i

From Univ. of Michigan Radiation Lab. Report #028918- 1-T A FINITE ELEMENT BOUNDARY INTEGRAL FORMULATION FOR RADIATION AND SCATTERING BY CAVITY ANTENNAS USING TETRAHEDRAL ELEMENTS J. Gong, J.L. Volakis, A. Chatterjee and J.M. Jin Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109-2122 ABSTRACT A hybrid finite element-boundary integral formulation is developed using tetrahedral/triangular elements for discretizing the cavity/aperture of microstrip antennas/arrays. The tetrahedral elements with edge based linear expansion functions are chosen for modeling the volume region and triangular elements are used for discretizing the aperture. The edge-based expansion functions are divergenceless thus removing the requirement to introduce a penalty term and the tetrahedral elements permit greater geometrical adaptability than the rectangular bricks. The underlying theory and resulting expressions are discussed in detail together with some numerical scattering examples for comparison and demonstration. ii

From Univ. of Michigan Radiation Lab. Report #028918-2-T Electromagnetic scattering and radiation from microstrip patch antennas and spirals residing in a cavity J. L. Volakis, J. Gong, and A. Alexanian Abstract A new hybrid method is presented for the analysis of the scattering and radiation by conformal antennas and arrays comprised of circular or rectangular elements. In addition, calculations for cavitybacked spiral antennas are given. The method employs a finite element formulation within the cavity and the boundary integral (exact boundary condition) for terminating the mesh. By virtue of the finite element discretization, the method has no restrictions on the geometry and composition of the cavity or its termination. Furthermore, because of the convolutional nature of the boundary integral and the inherent sparseness of the finite element matrix, the storage requirement is kept very low at 0(n). These unique features of the method have already been exploited in other scattering applications and have permitted the analysis of large-size structures with remarkable efficiency. In this report, we describe the method's formulation and implementation for circular and rectangular patch antennas in different superstrate and substrate configurations which may also include the presence of lumped loads and resistive sheets/cards. Also, various modelling approaches are investigated and implemented for characterizing a variety of feed structures to permit the computation of the input impedance and radiation pattern. Many computational examples for rectangular and circular patch configurations are presented which demonstrate the method's versatility, modeling capability and accuracy. iii

From Univ. of Michigan Radiation Lab. Report #028918-3-T A Finite Element-Boundary Integral Method for Cavities in a Circular Cylinder Leo C. Kempel and John L. Volakis December 22, 1992 Abstract Conformal antenna arrays offer many cost and weight advantages over conventional antenna systems. However, due to a lack of rigorous mathematical models for conformal antenna arrays, antenna designers result to measurements and planar antenna concepts for designing non-planar conformal antennas. Recently, we have found the finite element-boundary integral method to be very successful in modeling large planar arrays of arbitrary composition in a metallic plane. Herewith, we shall extend this formulation to conformal arrays on large metallic cylinders. In this report, we develop the mathematical formulation. In particular we discuss the shape functions, the resulting finite elements and the boundary integral equations, and the solution of the conformal finite element-boundary integral system. Some validation results are presented and we further show how this formulation can be applied with minimal computational and memory resources. iv

From Univ. of Michigan Radiation Lab. Report #028918-4-T A class of hybrid finite element methods for electromagnetics: a review John L. Volakis, A. Chatterjee and J. Gong Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor MI 48109-2122 Abstract Integral equation methods have generally been the workhoirse for antenna and scattering computations. In the case of antennas, they continue to be the prominent computational approach, but for scattering applications the requirement for largescale computations has turned researchers' attention to near neighbor methods such as the finite element method, which has low O(N) storage requirements and is readily adaptable in modeling complex geometrical features and material inhomogeneities. In this paper, we review three hybrid finite element methods for simulating composite scatterers, conformal microstrip antennas and finite periodic arrays. Specifically, we discuss the finite element method and its application to electromagnetic problems when combined with the boundary integral, absorbing boundary conditions and artificial absorbers for terminating the mesh. Particular attention is given to large-scale simulations, methods and solvers for achieving low memory requirements and code performance on parallel computing architectures. v

FINAL REPORT

A hybrid finite element-boundary integral method for the analysis of cavity-backed antennas of arbitrary shape Jian Gong, John L. Volakis, A.C. Woo and H.T.G. Wang Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor MI 48109-2122 Abstract An edge-based hybrid finite element-boundary integral (FE-BI) formulation using tetrahedral elements is described for scattering and radiation analysis of arbitrarily shaped cavity-backed patch antennas. By virtue of the finite element method (FEM), the cavity irregularities, the dielectric super/substrate inhomogeneities and the diverse excitation schemes inside the cavity may be readily modeled when tetrahedral elements are used to discretize the cavity. On the aperture, the volume mesh reduces to a triangular grid allowing the modeling of nonrectangular patches. Without special handling of the boundary integral system, this formulation is typically applicable to cavity-backed antenna systems with moderate aperture size. To retain an O(N) memory requirement, storage of the full matrix due to the boundary integral equation is avoided by resorting to a right triangular aperture grid on taking advantage of the integral's convolutional property. If necessary, this is achieved by overlaying a structured triangular grid on the unstructured triangular grid and relating the edge field coefficients between the two grids via a narrow banded, sparse transformation matrix. The combined linear system of equations is solved via the biconjugate gradient (BiCG) method, and the FFT algorithm is incorporated in the solver to compute the matrix-vector product efficiently, with minimal storage requirements. 1

1 Introduction Microstrip antennas have been extensively investigated experimentally, analytically and numerically for decades. By and large, numerical methods have been serving the engineers and researchers in the analysis and design of these conformal antennas for many years. Among them the moment method in conjunction with various integral equation (IE) formulations played a major role 1-3. However, IE methods are associated with field representations in which the appropriate Green's function for the specific geometry must be employed and this limits their versatility. Moreover, IE techniques are usually formulated on the assumption of an infinite substrate, a model which obviously deviates from the practical configuration leading to inaccuracies for larger bandwidth antennas. Furthermore, in the cqntext of IE methods, antenna excitations are represented using simplified models that differ more or less from the actual configurations. Also, due to the singularity of the current distribution near the patch-probe junction(s), special measures must be taken '4, not to mention additional IE complexities due to possible substrate anisotropies or inhomogeneities in the antenna substructure. In contrast, the hybrid Finite Element-Boundary Integral (FE-BI) technique alleviates these difficulties and this was demonstrated recently when the method was applied to rectangular patch antennas 5. In this paper, we present an edge-based hybrid finite element-boundary integral formulation using tetrahedral elements for a characterization of arbitrarily shaped cavity-backed antennas. An example of such a configuration is shown in Figure 1, where a cavity is recessed in a metallic ground plane enclosing the FEM volume and the antenna elements on the aperture may be excited by different schemes, such as a simple probe, a magnetic frill generator, a practical coaxial cable, microstrip lines, slots or a CPW line. In the context of the FEM, the cavity is first discretized into a number of tetrahedral elements that naturally reduce to triangles on the cavity's aperture. For nonrectangular patches this triangular gridding is, in general, non-uniform and the exact boundary integral formulation based upon this mesh applies to any patch shape. As a result, the hybrid FE-BI technique is capable of modeling arbitrarily shaped cavity-backed antenna configurations, different substrate inhomogeneities, anisotropies, as well as various practical excitation schemes. As is well known, the boundary integral (BI) equation subsystem leads to a fully populated matrix whose size is determined by the number of aperture mesh edges. For large apertures, this analysis becomes impractical in terms of storage and computation time requirements, and to overcome this inefficiency, a uniform zoning of the aperture is required as shown in Figure 2. By resorting to the structured mesh, the boundary integral matrix can be cast into a discrete convolutional form, thus permitting the computation of the matrix-vector products via the discrete Fourier transform (DFT) 2

avoiding a need to store the full BI matrix. This memory saving scheme has already been applied to IE solutions involving rectangular surface grids 5,6 and in this paper, we describe how the BiCG-FFT solution is implemented for triangular meshes. The differences between the rectangular and triangular meshes are also described and results are presented which demonstrate the method's vesatility in computing the scattering and input impedance of various non-rectangular printed antennas. 2 Formulation In this section, we present the edge-based hybrid FE-BI formulation using variational principles, where the matrix algebra notation is employed so that we can readily extend the formulae to the general anisotropic case. We begin with the vector wave equation, V x -V x E) - krE = -jkoZoJ' - V x (-Mi) (1) \^r ) \^r where J, and M, represent interior electric and magnetic current sources within the cavity V, C. and Pr denote, respectively, the relative permittivity and permeability, and ko is the free space wave number. The solution of (1) under certain boundary conditions is equivalent to extremizing the functional F(EH) = {(V x E) -(V x E)-korE Ed + jkoZo E x H dS + JJ E (jkoZoJi +V x -Mi) dv (2) where the surface S emcompases the cavity aperture excluding the portion occupied by the antenna elements. To proceed with the extremization of F, it is necessary to write H in terms of E. This can be accomplished through the surface equivalence principle. Following this procedure, the magnetic field H in the region exterior to the cavity may be expressed in terms of the equivalent magnetic surface current M = 2E x z (3) whose support is only over the extent of the surface S and the factor of 2 is due to image theory. Specifically, we have H(r) = HI + H - 2jkoYoj (I + -VV) Go(r,r') (E x z) dS' (4) Js \ Q 3

where H' is the incident field, if any, from the exterior region, H' is the reflected field due to H2 in the absence of the cavity, I is the unit dyad, and Go(r, r') is the free space Green's function with r and r' denoting the observation and integration points. Substituting (4) into (2) yields F(E) = JJJ{(V xE) (V E) kE.E d 2 JJJV + 2jkoZo J(E x H')~ dS + JJ E (koZoJi + V x -Mi) dv (5) - 2 (E {Jj(E x I) + -VV) Go(r, r') dS'} dS which is now only in terms of electric field. In deriving this, we used the relation n x H' = n x Ht which is valid when H' is measured on the ground plane. 2.1 Cavity Volume Modeling In proceeding with the discretization of (5), it is convenient to re-express it as F = Fv + Fs (6) where Fv denotes the volume integral contributions and similarly Fs accounts for the surface integral contributions. The cavity volume is subdivided into N tetrahedral elements Ve (e=1,2,...N), and within each tetrahedron the field is expanded as E = IVI{E} (7) with [V]e = [{Vx}{Vy}{Vz}e Vul {V) =., ux,y,z (8) El E2 {E} = kE e in which Vui is the u (u = x, y or z) component of the volume vector basis functions along the ith edge. The unknown vector {E}e has six entries, one 4

for each tetrahedron edge. (In this paper, we use square brackets for matrices and curly brackets for vectors). Inserting (7) into (5), and taking the first variation of Fv with respect to {E},, yields 6Fv = {[A]e{E}e + {I}e} (9) e where [A]e /= vi ~IDVi eDV' - koCEr'V' eV } dv (10) -J.Ve r.. Ji 1 M ix {K}e = //v [Vie {jkoZo (Ji + V x- Mi, dv (11) VTa a [DV]- = - iV (12) a a a{VY} - {VX} To carry out the above integrations, it remains to introduce the volume expansion or shape functions Ve. For our implementation we employed the linear edge-based shape functions for tetrahedral elements given in 7,8. 2.2 Aperture Modeling To discretize the surface integrals in (5), the aperture is subdivided into triangular elements since these correspond to the faces of the tetrahedrals. Within each triangle, the field is represented as E = [S]T{E}e (13) where [S]e = [{Sx}{Sy}]e Sul {Su} = Su2, u=x,y (14) Su3 Esl {Es}e = Es2 Es3 e 5

where Sui is the u(u = x, y) component of the surface vector basis functions along the ith edge. On substituting (13) into the surface integrals in (5) and taking the first variation of Fs with respect to {Es}, we obtain 5Fs =- {:B:e{Es}e + {L}e} (15) e where [B]e = J {-2koSe:se + 2 [{Sy - {Sx}] [{S} *Go(r, r') dS dS' (16) and {L}e = j2koZo [Se'( Hy ) dS (17) Note that in (16) the elements of the array [Se] are functions of the observation vector r, whereas the elements of Se T are with respect to the integration point r'. A suitable set of linear edge-based surface basis functions is ii. S/(r)= 2A z x (r- ri)e(r) r E Se (18) Si(r) = 2Ae (18) 0 otherwise In this expression, li denotes the length of the ith edge and ri is the position vector of the vertex opposite to the ith edge. Since each edge shares two triangles, one is defined as the plus and the other as the minus triangle. Therefore, e(r) is given by (r) { 1 r E S (19) where Se = S+ + S7. The constant Ae in (18) denotes the area of the plus or minus triangle depending on whether r E S+ r r S_. We note that Si(r) x x yields the basis functions used by Rao, etc. [9] in their moment method solution of boundary integral equations. 2.3 System Assembly To construct the final system for the solution of the electric field components we combine (9) and (15), and after assembly we obtain the system {[A]{E} + {I}} + {[B]{E,} + {L}} = 0 (20) In this, {K} and {L} are the excitation vectors due to the interior current sources and the exterior excitation, respectively. The unknown electric field vector {E} consists of all field expansion coefficients with respect to the 6

element edges except those coinciding with perfectly electrically conducting (PEC) walls, PEC antenna element(s) or PEC pins inside the cavity. Finally, the vector {E,} represents the unknown surface fields whose entries are part of those in {E} with their corresponding edges on the aperture. Explicit expressions for the matrices and vectors in (20) are AC = J (-(V x V) (V x V) - k2cV * Ve) dv (21) ~= 2k ls [Z x S (r)] * {J [S (r x ] Go(r, r') ds'} ds (22) +2 /se Vs* (S(r) x Z) {JJ VS [S. (r) x Z] Go(r, r') ds' ds K = JJJ V {jkoZoJt + V x (-M?) } dv (23) Li = 2jkoZo S * [Hxi n] ds (24) when specialized to inhomogeneous, isotropic dielectric cavity fillings. It is evident, from (21) and (22), that [A] and [B] are symmetric as a result of the assumed isotropic medium and reciprocity. In addition, 'A' exhibits high sparsity due to the FEM formulation whereas IB' is fully populated. Two approaches may be followed in carrying out the solution of the combined subsystems when an iterative solver is employed such as the biconjugate gradient (BiCG) method [10]. These two approaches differ in the manner used for the evaluation of matrix-vector products called for in the BiCG algorithm. One could sum the coefficient matrices [A' and [B by adding up the corresponding matrix entries prior to the execution of the BiCG algorithm, or instead the resulting vectors may be summed up after carrying out the matrix-vector products. We observed that the first approach is more efficient in terms of computation time after reordering the combined matrix and storing only the non-zero elements. This is because in the context of this scheme the combination of the two matrices is performed only once outside the BiCG iteration. However, the second approach is compatible with the BiCG-FFT scheme, when the FFT algorithm is employed to exploit the convolutional property of the integral operator, thus eliminating a need to explicitly store the entire BI matrix. Below, we discuss the implementation of the matrix-vector product of the boundary integral system for the BiCG-FFT solution. 2.4 Implementation of the Boundary Integral Matrix Vector Products Using FFT To describe the execution of the matrix-vector product in the context of the BiCG-FFT solution scheme we refer to Figure 2. In this figure we display the 7

uniform triangular grid which either overlays the irregular aperture grid or is simply the aperture grid which the volume tetrahedral mesh is by design reduced to during the discretization process. We recognize that the proposed triangular grid consists of equal right triangles and thus it involves three different classes of edges (class 1, 2 and 3). These include the x-directed, y-directed and the diagonal edges, all of which are uniformly spaced. For the FFT implementation each class of edges is independently numbered in accordance with their geometric location. Specifically, the ith class will carry the numbering (m, n) if the edge is the mth along the x direction and the nth along the y direction. The indices (m, n) take the values m = 0,1,2,...,Mi n = 0,1,2,...,N' with i = 1 for the y-directed edges, i = 2 for the diagonal edges and i = 3 for the x-directed edges. Consequently, we find that M-2 i = 1 N-1 i=1 M- M- 1 i=2 N N-1 i =2 (25) M- 1 i=3 N-2 i = 3 where M and N denote the numbers of elements along the x and y directions, respectively. To perform the integrations for the evaluation of the boundary integral matrix elements, it is now convenient to rewrite the basis functions (18) in terms of the new indices (m, n). We readily find that the edge-based basis functions associated with each of the aforementioned class of edges can be rewritten as (nAy - y)5x + (x - mAx)y (x, y) e S~+ S(x,y) =(y- (n + l)y)x + ((m + 2)x - x) (x, y) e Se (26) 0 otherwise (x2+( 2 (nAy- y)x + (x - (m + 1)Ax)y e S+ mn(x Y) y (y - (n + l)Ay) + (m - e)y S; (27) 0 otherwise 1 f ((n + 2)Ay - y)x + (x - (m + l)Ax)y (x, y) E S+ ^sz( ) = a (Y - nAy)& + (mAX - x)y (, y) e S (28) -.,(x y) (x, y) e S- (28) O otherwise where the superscripts refer to the edge class. Each entry of the boundary matrix-vector product can now be calculated as 3 Mj NJ {BI subsystem} = [B] {E,} = B,E E E nmn (29) j=,l m'=O n'=O in which (m, n) are the geometric location indices for the ith class edges of observation elements whereas (m', n') are the same for the jth class edges of 8

testing elements. Thus, the specification of the indices i, m and n completely defines the entry ki = nMi + M of the column resulting after the execution of the boundary matrix-vector product. It is readily found that Bnm = -2k JJ J S. S' G(r,r')dx dy dx' dy' (30) + (A y)2 S JJ i(r) ) j(r)li Go(r, r') dx dy dx' dy' with fy ili= /(A)2+ (y)2 i= 2 (31) Ax =i-=3 More importantly, it can be shown that the BI subsystem exhibits the convolutional property B,m'n, = B( -m'nn) and thus we can rewrite (29) as 3 [B {E} = EBij * Ej (32) j=1 where the * denotes convolution. It is now seen that the computation of the boundary matrix-vector product can be performed by employing the 2-D discrete Fourier transform (DFT) thus avoiding a need to store the BI matrix other than those entries which are unique. When the symmetry property of B( -ml_,_ is also invoked, implying B1' = B31 (33) (m-m',n-n') (m'-m, n'-n) it is concluded that the total non-redundant entries in the BI matrix are 3 3 Np = N (M + M - 1) (34) i=1 j=1 This should be compared to the (z=1l MiNi)2 entries whose storage would normally be required if the BI system was not cast in convolutional form. We remark that Np is nevertheless equal to twice the number of entries required for uniform rectangular grids [6; for one class of edges. To avoid aliasing, it is necessary that B ', -n') = BpJ(imi, ) be cast in a 2-D array which has the usual periodic form and zero padding may also be required to make use of the standard FFT routines. Specifically, the matrix-vector product (29) is 9

executed by using the MFT xNFT array B j 0 < im < M' (- -) 0< i < N' MFT- M' + 1 < im < MFT B'(-m,' -4), 0 < i < Ni B(,? B ii - 1 - NFT) 0 < < Mi m, NFT- NJ + 1 < n < NFT MFT- M' + 1 < fih < MFT Bi (mn - 1 - MFT, ii - 1 - NFT), MFT - M + 1 < m < MFT NFT- NJ + 1 < n < NFT 0 otherwise (35) with the corresponding field vector given by EJ'(m, ) 0<m<M, 0O<h< N EPJ- (i' ) - 0 n n (36) E,)={ O n otherwise (36) and MFT and NFT must be powers of 2 if a radix 2 FFT is used. In the BiCG-FFT algorithm the BI subsystem vector is computed as 3 {BI subsystem} = S {DFT- {DFT{BpJ} * DFT{Ep}}} (37) i=l The presence of the operator S indicates the necessary reordering of the 2-D array which results after the inverse FFT operation into a single column with the proper indexing for addition to the FEM subsystem. 3 Mesh Overlay Scheme for Nonrectangular ]Patches As described above, the BiCG-FFT solver requires uniform aperture gridding so that the BI subsystem can be put in block circulant form. This can be always achieved during mesh generation whenever the patches are rectangular in shape or in case of radiators which are placed at some distance (usually small) below the aperture. However, for circular, triangular, or other nonrectangular patches on the aperture, it is not possible to construct a uniform mesh using the mesh generator. Typically, the aperture mesh is necessary to comform to the patch shape, leading to an unstructured surface grid. In this case, to make use of the efficient, low memory BiCG-FFT algorithm, an approach is to overlay on the unstructured aperture grid another coincident structured grid as shown in Figure 3. The boundary integral subsystem is 10

then constructed by using the overlaid uniform grid whose edge fields can be related to those on the unstructured grid via a sparse transformation matrix. That is, it is necessary to append to the system (20) the relations {E, } = ITF: {E. } {Es}nu = [TB] {E.}) (38) where the subscripts u and nu refer to the field coefficients of the uniform and non-uniform aperture grids, respectively. Also, [TF] and [TB] refer to the forward and backward transformation matrices, respectively, with Nu and Nnu denoting the numbers of the uniform and non-uniform mesh edges on the cavity aperture. To derive the elements of [TF', we begin with the expansion (13) and enforce it at three points on each edge belonging to the uniform grid. We conveniently place these three points at the center and ends of the edge (see Figure 4). Given the fields at these points, we can approximate the field along the (m, n) edge of the uniform grid using the weighted average 1 e 1 Nend1 (Eu)(m,n) -= Enu 2Nd E (i) 2 Nendi k=1 1 Nmid + Nmid Enu u(rmid) Nmid k=1 1 Nend + E (rend ) (39) 2Nend2 k=1 in which eu denotes the unit vector along x,y or the diagonal, depending on the class of edge being considered. The quantities Ek represent the fields in the non-uniform grid triangles with the superscript k being a sum variable in case rend,, rend2 or rmid specify a point shared by more than one triangle. Obviously, Nend,, Nmid and Nend2 denote the number of non-uniform grid triangles sharing the node at rend,, rmid and rend2, respectively, and will typically be equal to unity. After assembling (39) into (38) we find that the elements of the forward transformation matrix are given by ^ ( Z Nend, 3 (TF)ij - -eu Sk U2Wnd y o6 'o ijtS^r^end ') 2 j 2Nendi k=l e=1 1 Nmid 3 +Nid El E {ijtS(rmid) mid k=1 e=i 1 Nend, 3 + E eij S (rend2) (40) 2Nend2 k=l e=l 11

in which - = 1 j-fi =3t 0 otherwise and the global indices i and j correspond to the ith uniform grid edge and the jth non-uniform grid edge. The subscripts jie is the global index used in numbering the non-uniform grid edges, whereas the subscript e (= 1, 2 or 3) is the local edge index used in the definition of the basis functions St (see (13)). We remark that the explicit computation of the transformation matrix elements results in a substantial programming simplification because it avoids the usual assembly process. Following the same procedure, we can obtain the elements of the backward transformation matrix. Assuming each uniform grid edge traverses on three or less non-uniform grid triangles, the non-zero entries in each row of [TF] will be 9 or less. However, they can reach a maximum of 18 if the midpoint and endpoints reside on an edge of the non-uniform grid. The maximum non-zero entries in each row of [TB' will be 15, but the typical number will be much less. 4 Numerical Considerations Based on the presented FE-BI formulation, a computer program was written for the analysis of the radiation and scattering by cavity-backed patch antennas of arbitrary shape. The antenna geometry is supplied to this program in an input file which as a minimum, it must contain lists of (a) the nodes and their (x, y, z) coordinates, (b) the nodes forming each tetrahedron, (c) the nodes on the cavity aperture, and (d) the nodes on metallic boundaries. For arbitrary antenna geometries, it is necessary to employ a sophisticated volume mesh generation package and a number of these are available commercially. Typically each of these packages generates a "Universal file" which can be readily preprocessed to extract the aforementioned input lists. A major effort was devoted in writing the program in a manner which minimizes the storage and computational requirements. Specifically, the BCs on the metallic surfaces are enforced a priori to obtain a system which involves only nonzero field components. The sparse finite element matrix was stored as a single array of length NNnz where Nv is the total number of unknowns within the cavity volume and Nnz denotes the maximum number of nonzero row entries. The BI matrix was stored in different ways depending on whether the FFT was to be employed for the evaluation of the matrixvector products. If the BiCG solution was to be carried out without the FFT, then the NS x N, BI integral matrix was added to the FE array resulting in a 1-D array about NNnz + N2 long. For slot antennas, including cavity-backed spirals, and moderately sized systems, it was found preferable not to use the FFT, thus avoiding any discretization errors. In that case the generation of a single combined FE-BI matrix before execution of the 12

BiCG algorithm reduces the computational requirements because a number of operations associated with the repeated combinations of the FE and BI subsystems within the BiCG iteration is avoided. When the FFT is to be used as part of the BiCG solver the FE and BI matrices must be kept in separate arrays throughout the execution process. In this case the FE matrix is again stored as a single array and similarly the non-redundant elements of the BI matrix are stored in a single array of length 9N(2M - 1). The factor of 9 is due to the three classes of edges and as usual M and N denote the number of elements along the x and y directions, respectively. Because of the storage and computational efficiency of BiCG-FFT algorithm, it is necessary to resort to uniform aperture grids for conformal antennas involving substantial number of aperture edges. Of course, one should always use uniform triangular grids when the patches on the aperture are rectangular or if the array is simply covered by a superstrate. In the case of circular patches it will be necessary to overlay a structured triangular grid over the unstructured grid generated by the mesh generator. This must be done in the preprocessing stage and should be taken into account when constructing the FE matrix. For scattering computations, the overlay of the structured grid is almost always the preferred approach because it does not, generally, compromise the accuracy of the computed scattering cross section. However, for antenna parameter computations, the interpolation scheme between the structured and unstructured grid edges may be of concern, depending on the specific antenna geometry. Generally, circular slots and spirals should be treated without resorting to structured grids and to our experience this does not cause large computational burden because these antennas are associated with small apertures. In the case of circular patches, the structured grid was not seen to compromise the computational accuracy. Of course, conclusions based on one type of antenna do not necessarily apply to others, and thus the suggested alternatives must be examined separately for each antenna before choosing one approach over the other. Of importance here is that the formulation is suitable for modeling any antenna shape and feed structure. 5 Results We present below some representative numerical results for the purpose of validating and demonstrating the robustness of the tetrahedral formulation for scattering and radiation by different configurations of cavity-backed antennas. In each case the computed results via the FE-BI method are compared with reference measured or calculated data. BiCG-FFT algorithm validation: To validate the BiCG-FFT implementation for the boundary integral matrix-vector products described in section 2.4, we consider the scattering by a circular patch placed 0.325 cm be 13

low the aperture of a rectangular cavity (recessed in a ground plane) of 7.8 cm x 5.2 cm x 0.975 cm. The geometry is illustrated in Figure 5 and we remark that the patch was purposely submerged so that the structured grid using a mesh generation package was available on the aperture. A bistatic radar cross section (RCS) pattern is shown in Figare 5(a) as computed with and without use of the FFT in the BiCG algorithm for matrix-vector products. As seen, the two solution procedures give identical scattering patterns for two different cavity dielectric filling. Overlay scheme validation: To validate the structured/unstructured mesh overlay scheme discussed earlier, we calculate the scattering by a circular patch (0.1 cm radius) placed on the aperture of a rectangular cavity 0.4 cm per side and 0.1 cm deep. The bistatic radar cross section (RCS) patterns are shown in Figure 5(b) as computed with and without use of the overlay scheme for normal and 30~ off ground plane incidence. The dielectric fillings were chosen as Cr = 4 and ~r = 2 for normal and oblique incidence, respectively. Of course, the FFT was only employed in connection with the overlaid structured mesh. Scattering by a circular patch: Figure 6 illustrates f circular patch residing on the surface of a 0.406 cm thick substrate having a relative dielectric constant of ~r = 2.9. The patch's diameter is 2.6 cm and the substrate is enclosed in a circular cavity 6.292 cm wide. This cavity and the patch were recessed in a low cross section body for measuring its RCS. A comparison of the measured and calculated backscatter aos RCS as a function of the frequency is shown in Figure 7. For this computation the direction of the incident plane wave was 60~ from normal, and as seen the agreement between measurements and calculations is very good throughout the 4-9 GHz band. Input impedance measurements and calculations for the same patch are displayed in Figure 7. The feed in this case was placed 0.$ cm from the patch's center and it is again seen that the measurements and calculations are in good agreement. Scattering by a square archimedean spiral: Figure 8 displays a square archimedean spiral. A square rather than a circular spiral was selected because it allowed comparisons with calculations based on a different technique. The subject complimentary square spiral consists of strip arms, each of width 0.09375 cm, and is placed (free-standing) on the aperture of a square air-filled cavity 2.812 cm per side and 0.9375 cm deep. The q = 45~ plane ass and a,, bistatic (incident plane wave was at 30~ off normal) RCS patterns for this structure are shown in Figure 8. It is seen that the ass + cao RCS compares well with the corresponding data based on the finite difference-time domain (FD-TD) method [13]. No FD-TD data were available for the a+ RCS pattern but as expected, the ao, RCS return is much lower and vanishes at grazing. Annular slot impedance: Figure 9 shows a narrow circular (0.75 cm wide) annular slot situated in a circular cavity 24.7 cm wide and 3 cm deep. Be 14

is very small for this application and as a result there is no need to invoke the FFT in the BiCG algorithm. The FE-BI method is basically quite effective in modeling small aperture configurations without a need for special computational considerations. Input impedance calculations as a function of frequency for this radiator, excited by a probe placed across the slot, are shown in Figure 9, and agree well with the values calculated via a modalboundary integral method 14'. For these calculations, the frequency was swept from 700-1000 MHz. The dielectric constant of the material filling the cavity was set to e, = 1.35 as in 14' and this is an effective value to account for the presence of a dielectric slot cover used as part of the measurement model for holding the plate. Radiation by a one-arm conical spiral: We considered the modeling of this radiator to demonstrate the geometrical versatility of the FE-BI method. A configuration of the spiral radiator and surface mesh is illustrated in Figure 10 from different views. The top and bottom edges of the strip forming the spiral follow the lines p = 0.0503Aexp 0.221(<q i 2.66), z = a~exp(0.221q), where (p, X, z) denote the standard cylindrical coordinates, a~ are equal to 0.0832A and 0.0257A, respectively, and 0 < q < 27r. This spiral arm resides on an inverted cone (9.24 cm tall) whose bottom cross section has a diameter of 1.68 cm and the top cross section has a diameter of 21.78 cm. For our calculations A = 30 cm (f = 1 GHz) and the spiral was situated in a circular cavity 10.01 cm deep. The computed E, radiation pattern taken in the = 90~-plane, using a probe feed at the cavity base, are given in Figure 11. It is seen that the Ee principal plane pattern is in good agreement with the data given in [12]. However, the Eo pattern differs from the measured data primarily because of the circular cavity included in the analytical model. The latter was not part of the measurement configuration which consisted of the spiral antenna on a large circular plate. 6 Conclusion We presented a hybrid finite element-boundary integral (FE-BI) formulation which incorporates linear tetrahedrals. The method was specifically developed for the radiation and scattering analysis of cavity-backed printed antennas, where the FEM is used for modeling the cavity region and the BI equation acts as a global boundary condition for terrinating the mesh on the cavity aperture. The FE-BI formulation is particularly suited for the analysis of complex configurations and much emphasis was given here in developing a solution technique requiring O(N) storage in spite of the resulting full BI subsystem. The latter was achieved by makirg use of the convolutional property resulting from the structured mesh, thus permitting use of the FFT in the BiCG solver for computing the matrix-vector products. For scattering calculations associated with large aperture structures, use of the 15

FFT proved essential in minimizing the computational requirements. A number of patches, slots as well as planar and non-planar spiral antennas were analyzed for the purpose of demonstrating the versatility and accuracy of the FE-BI technique. Certainly, the need to use a sophisticated mesh generation package is deterrent to the application of the technique for the analysis of simple antenna configurations. However, this is unavoidable when dealing with complex geometries and, moreover, the pervasive use of such commercial packages on desktop computers makes the technique quite attractive. References [1] E.H. Newman and P. Tulyathan, "Analysis of microstrip antennas using moment methods," IEEE Trans. Antennas and Propagat., Vol. AP-29, pp. 47-53, Jan. 1981. [2] J.R. Mosig and F.E. Gardiol, " Analytical and numerical techniques in Green's function treatment of microstrip antennas and scatterers," IEE Proc., Pt. H, Vol. 130, pp. 175-182, 1983. [3] M.C. Bailey and M.D. Deshpande, "Analysis of elliptical and circular microstrip antennas using moment methods," IEEE Trans. Antennas Propagat., Vol. AP-33, pp. 854-959, Sept. 1985. [4] D. Zhang and K.A. Michalski, " Analysis of coaxial fed antennas of arbitrary shape with thick substrates," J. Electromagn. Waves Appl., Vol. 25, No. 12, pp. 1303-1327, 1991. 5] J.M. Jin and J.L. Volakis, "A hybrid Finite Element Method for scattering and radiation by microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas Propagat., Vol. 39, No. 11, pp. 1598 -1604, Nov. 1991. 6 J.M. Jin and J.L. Volakis, "A biconjugate Grad ent FFT solution for scattering by planar plates," Electromagnetics, Vol. 12, pp. 105-119, 1992. [7] M.L. Barton and Z.J. Cendes, "New vector finite elements for three dimensional magnetic field computations," J. Appl. Phys., Vol. 61, No. 8, pp. 3919-3921, 1987. [8] A. Chatterjee, J.M. Jin and J.L. Volakis, "Computation of cavity resonances using edge-based finite elements," IEEE Trans. Microwave Theory Techn., Vol. 40, pp. 2106-2108, Nov. 1992. 16

[9] S.M. Rao, D.R. Wilton and A.W. Glisson, "Electromagnetic Scattering by surfaces of arbitrary shape," IEEE Trans. Antennas Propagat., Vol. AP-30, pp. 409-418, May 1982. [10] T.K. Sarkar, "On the application of the biconjugate gradient method," J. Electromagn. Waves Appl., Vol. 1, No. 3, pp. 223-242, 1987. [11] A.K. Jain, Fundamentals of digital image processing, Prentice-Hall: Eaglewood Cliffs, NJ, 1989. [12] D.W. Smith and P.E. Mayes, "Numerical and experimental analysis of circularly polarized radiating line antennas," Electromagnetics Laboratory Report No. 90-4, University of Illinois, 1990. [13] C.W. Penney and R.J. Luebbers, "Radiation and scattering from a square archimedean spiral antenna using FD-TD," Electromagnetics, 1994 (to appear). [14] H. Morishita, K. Hiraswa and K. Fujimoto, "Analysis of a cavity-backed annular slot antenna with one point shorted," IEEE Trans. Antennas Propagat., pp. 1472-1478, Oct. 1991. 17

patch ground plane aperture: BI area cavity: FEM vc excitation cable Figure 1. Illustration of a cavity-backed patch antenna and the different computational regions. yt 3 2 1 0 n - x m: 0 1 2 3 4.. Figure 2. Illustration of the uniform (structured) right triangular grid. 17

, circular patch Figure 3. Overlay of a structured aperture mesh over an unstructured mesh, shown here to conform to a circular patch. 18

A e u II -I -I $1 %^ I % I \ ^ overlaid unstructured mesh triangles n2 (m, n) uniform grid edge X\ / Figure 4: Illustration of the parameters and geometry used in constructing the transformation matrix elements between the structured and unstructured mesh. 19

4y 0.325 cm 10.65 cm 7.81cm -q a -I I -T 2.6 5.2 cm — I-0 x l 0.975 cm 3.9 cmr 7.8 cm 7.8cm 4.1I 7.8 cm I perspective view top view Cd) O co...4-4 a, ~,,, 30. 20. 10. 0. -10. -20. -30. -40. -50. -60. 0. 10. 20. 30. 40. 50. 60 70. 80. 90. Figure 5(a). Comparison of the ao* bistatic RCS patterns for the illustrated circular patch as computed by the FE-BI method with and without the mesh overlay scheme. The computation frequency was f = 6 GHz with normal incidence. 20

0. -10. s,~ -:s. uA AA -20. c -25. I i n oblique incidence(no OverlaY) X i 0 | -normal incidence (no overlay) I. -35.. o.nliqun incidence (nooverlay) ^ oblique incidence (overlay) -40. norma del J -45. go 90. -50. 430. r4 20. 30. 50. 5() 20.- cm (residing ve e nc ae o at = 00 (= 4) and the other 21

u.uu...... ' ' ' i " ' "T. I ' '' ' t ' '.....'.. - Measurement -J -122 uo FE-BIMeJlod 0o -32.-2 I46=2.6ccm -2u.|I 2x3.446cm - 5.C.................... 4.03 5.C 6. 15. 7. u D 8.u frequency [GHz) Figure 65. Comparison of i;he compui;ed and measured ao9 baclscal;i;er RCS as a functlion of frequency for the shown circular patch The incidence angle was 30~ off the ground plane. 22

1.0 0.5 -— 1.0 -— 2.0 OX- 5.Q 5.0 0.5 2.0 1.0 Figure 7. Comparison of the computed and measured input impedance for the circular path shown in Figure 5. The feed was placed 0.8 cm from the center of the patch and the frequency was swept from 3 to 3.8 GHz. 23

Coordinates of cross (X at lower right corner of conducting strip): (1.406 cm, 1.406 cm) 0- 0.09375 cm 2.8125 cm 0.09375 cm I I I 2.8125 cm Il 10. C14 t) N% — m1 Cl CL V) uf 0. -10. -20. -30. -40. -50. -60. -70.. E-BI) is i... -— aoe + a,, (FE-BI) a<,* (E-BI) A aee + a++ (FDTD) 11 ~1 1............ Qn E -OU. 0. 10................................................ I. I..... 20. 30. 40. 50. 60. 70. 80. 90. e Figure 8. Geometry and bistatic RCS calculations at f=4.5GHz for a square archimedean spiral situated on the aperture of a square 2.8125 cm per side and 0.9375 cm deep cavity. 24

0.f I f t I I i t I I \ probe 0~ o r r r r r.f a = 12.35 cm b = 0.75 cm;o=7.7 cm 0.7 <f< 1GHz d=3cm C E r =1.35 I _ 1.0 - - - Calculated [14] * FE-BI Method 0 1.0 Figure 9. Comparison of input impedance calculations for the illustrated cavitybacked slot. 25

Figure 10. Illustration of the configuration and mesh of the one-arm conical spiral used for the computation of Figure 11 (the mathematical definations of the spiral edges are given in section 5). 26

-5. 0A | — 15. o ' ' E pattem (FE-BI metho.). / ----- Ee pattern (FE-BI method) v: -25. v A E pattern [12] Xv v Ee pattern [12] -35. ll -90. -60. -30. 0. 30. 60. 90. Figure 11. Comparison of the calculated Ee and Ek radiation patterns taken in the 4 = 900-plane with data in reference [12] for the one-arm conical spiral shown in Figure 10. 27