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 workhorse 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. 1

Introduction The 1980's brought to focus a variety of computational methods for electromagnetics applications. Journal publications in this area grew dramatically and a variety of numerical codes were developed in industry and academia for analysis and design purposes. The methodologies used for developing these codes vary widely, but in general, they can be categorized in (a) moment method and (b) partial differential equation (PDE) formulations. Both were actually introduced in the sixties and seventies but the dramatic increase in computing resources over the past decade led to the development of accurate, useful and in many cases practical numerical codes for electromagnetic computations. Useful computer codes for wire (antennas and scatterers) and two-dimensional structures have been available since the seventies [1, 2, 3]. However, the more arduous task of developing useful codes applicable to complex three-dimensional structures took place, primarily, in the later part of the eighties [4]-[18] and continues at a rather fast pace. Driven by the capabilities of computing resources, a variety of basis functions, discretization schemes and solvers have been employed in connection with 3D moment method codes, all aimed at reducing the storage requirements and speeding-up the solution algorithms on either serial or parallel computers. To a great degree, though, leaps in computer speeds turned reseachers' attention to method which havey the lowest memory requirements, even at the expense of speed. PDE methods have inherently low memory requirements whether executed in time or frequency domain and have, thus, been courted by many research and engineering groups in industry, government and academia. Numerous successful applications of the Finite Difference-Time Domain (FD-TD) [11]-[13], [16], [19]-[21] and Finite Element Methods (FEM) [15, 17], [22]-[27] have been reported in the literature and by all accounts, the proliferation of the methods will continue to grow. Such methods have been routinely used for electrostatic and magnetostatic analyses, but during the past decade they were also developed for large scale computations involving hundreds of thousands of unknowns. Over the past few years, the group represented by these authors has been developing frequency domain PDE methods for scattering and antenna characterizations [17], [27]-[38]. In this paper, we review our progress with emphasis on three-dimensional applications. As is the case with any PDE method, the mesh truncation scheme plays an important role in assessing the accuracy and efficiency of the implementation. Various mesh termination schemes are mentioned here and their attributes are discussed. Generally, the mesh termination scheme is either exact or approximate and each truncation scheme leads to a different hybrid FEM. Exact mesh truncation schemes have the drawback of destroying the sparsity of the finite element matrix but in most cases they also reduce the computational domain. Approximate truncation schemes rely on the use of absorbing boundary conditions or material absorbers whose purpose is to suppress wave reflections back into the computational domain. They retain the sparsity of the finite element matrix but one is required to increase the computational domain to ensure acceptable levels of accuracy. For antenna parameter prediction, there is a need for accurate near field calculations and because the absorbing boundary conditions (ABCs) lead to greater inaccuracies in the near zone, it is best that such predictions are carried out in connection with exact truncation schemes. Below, we first summarize the finite element formulation without reference to the mesh termination scheme. This formulation accounts for the presence of materials, conductors, 2

antenna feeds, resistive and impedance cards. The discretization of the weak form of the wave equation is then presented, followed by a discussion on the mesh truncation scheme and the properties of the resulting system. Both the boundary integral (BI) and ABC termination schemes are presented with specific reference to the resulting system solvers. When a boundary integral is used to terminate the mesh, the iterative solver must be combined with the fast Fourier transform to maintain the storage requirements to O(N). In the case of the finite element-ABC (FE-ABC) formulation, we discuss the code implementation and its execution and performance on parallel computing architectures. Finally, several examples are presented which demonstrate the utility and accuracy of the discussed hybrid FEM formulations for calculations pertaining to single element antennas, finite and infinite periodic arrays, and composite scatterers. Although this paper is certainly a review of hybrid finite element formulations, most of the included results have not appeared in the literature. Formulation The variational functional Consider the cavity configuration, shown in Figure 1, and the inhomogeneous scatterer shown in Figure 2. The cavity's aperture shall be denoted by So whereas the same symbol will be used to denote the termination of the PDE's computational domain in connection with the scatterer in Figure 2. As shown, the cavity is recessed in a ground plane and is intended to house the radiating elements and their feeding structure such as striplines, coaxial cables or microstrip lines. Any material filling can be handled by the presented formulation including resistive cards [34], lumped or distributed loads [33] often used for frequency and bandwidth control. Inhomogeneities are also assumed to be present in the dielectric scatterer of Figure 2. A goal with any PDE method, such as the finite element method, is the solution of the vector wave equation (1V x E)- k^6E = -jkoZoJ - V x ( Mi) (1) in which E represents the total electric field, (Er,,r) denote the relative permittivity and permeability of the computational domain, ko is the free space wavenumber and Zo is the free space intrinsic impedance. In addition, J' and M' denote the impressed electric and magnetic sources, respectively. Instead of solving (1) directly, for numerical accuracy and stability it is best to consider the weak form of (1). This is obtained by first multiplying (1) by a suitable weighting function and subsequently integrating over the computational domain [17]. An alternative approach is to consider the minimization of the functional [34] F(E) = JJ -(V x E) (V x E) - ko E rE] dv + J J E [jkoZoJ + V x -M)] dv JJJV. [ \^~~r ) 3

+ 2jkoZo fs R( x E) (n x E)ds + jkoZo J E (H x )ds (2) in which R denotes the resistivity or impedance of the surface Sr within So, and H is the corresponding magnetic field. We note that Sr must be closed when it satisfies an opaque impedance boundary condition but can be open (i.e. a finite plate) if it satisfies the penetrable resistive sheet condition [39] n x n x E =-Rn x (H+- H-) where HA denote the fields above and below the surface Sr. As seen from (2), no special care is required for the evaluation of the integral over Sr, in spite of the discontinuity in the magnetic field. The explicit knowledge of H is, however, required over the outer surface So for the construction of a discrete system of equations for the solution of E. Before proceeding with the latter, we note that for large computational domains it may be appropriate to reexpress the functional (2) in terms of the scattered field. This reduces mesh propagation errors [40] and facilitates the implementation of the ABCs on So. Setting E = E + Ecat, H H Hi + Hscat where E' denotes the excitation field impinging from the exterior of So and Escat is the scattered field, the variational functional F takes the form F(E)sc') = 21 [ l( x Eat). (V x ESat) - k2 Escat Escat dv + JJJ EScat [koZoJ + V x (M')] dv 1 1 + 1jkZo/j l R( s x Escat). (n X Escat)ds Sr + jkoZo J | (| x Escat). (An x Ei) ds + If! [-(V x Escat) (V x E) - k2crEcat E'] dv + jkoZo iEscat ( x Hi) ds - kZo I[ Escat. ( x Hscat) ds I'kO ZOJ~ 1 E SCat..(n XH)ds + f(E') (3) In this, f(E') is a function of the incident field whose specific form is of no consequence in the subsequent analysis. Most importantly, we observe in (3) the presence of the integrals over Sd and Vd, where Sd denotes the aggregate of all surfaces in V which coincide with a dielectric interface and Vd encompasses the volume of the enclosed dielectrics. The introduction of these integrals is required when working with the scattered field because at an interface Escat is associated with a jump in permeability. 4

Discrete system To construct a system of equations from (2), the volume V is first subdivided into a number of small elements such as those shown in Figure 3. Each of these small elements occupies a volume Ve (e = 1,2, 3,..., M), where M denotes the total number of elements filling V. The field (total or scattered) in each element is then approximated by a linear or higher order expansion, viz. m E = z EjWj = {We}T{Ee} (4) 3=1 in which W^ are the (vector) expansion basis functions [41]-[43]. With this representation, the unknown expansion coefficients Ej represent the field at the jth edge of the eth element having m sides/edges. The edge-based expansion (4) has zero divergence within each element as required by Maxwell's equations and avoid an explicit specification of the field at metallic corners. Corresponding node-based field representations, which do not satisfy the zero divergence condition, may yield extraneous solutions unless augmented in some way by adding, for example, a penalty term to (2) or (3) [31], [44], or by instead working with potentials rather than the fields themselves [25]. Recently, however, node-based representations have been proposed which remove this difficulty by enforcing the divergence condition in the process of constructing the node-based basis functions [45]. The calculations to be presented here were carried out using a first order edge-based expansion in connection with tetrahedrals (also referred to as 1-form Witney elements [43]) [35, 37, 27], rectangular (right-angled) bricks [32, 33] or cylindrical-shell elements [46]. Edgebased implementations using curvilinear bricks (also referred to as parabolic elements) have also been reported [47], [48]. It should be noted, nonetheless, that both edge-based and node-based expansion functions have their advantages and disadvantages [49]. Specifically, node-based functions are often found easier to implement and the first order element provides a linear representation of the field in all directions. However, apart from their difficulty with spurious solutions (a major restriction), they require double noding or other special care at dielectric interfaces and lead to larger matrix bandwidths. The edge-based elements do not require double-noding at dielectric interfaces, and usually lead to sparser matrices. However, the unknown count is generally doubled when using edge-based elements but this is balanced by the increased sparsity of the stiffness matrix [27], [49]. Another drawback appears to be their inherent assumption that the field be constant along the edges forming the geometrical elements. The system of equations to be solved for Ej is obtained from (2) or (3) by a Rayleigh-Ritz procedure [22]. This amounts to differentiating the functional F with respect to each edge field and then setting it to zero. On substituting (4) into (2) or (3), taking the first variation of F with respect to E or Escat and assembling all elements, we obtain the system f OF ) M M N8 MP Ee ) = Z[AX]{Ee + E[Bs]{Es} + E{g8} + E{CP} = (5) e=1 s==1=1 p=l The sum of the element submatrices [Ae] results in a large sparse matrix and its elements are computed from the volume integral(s) over V in (2) or (3). The column {CP} results from the discretization of the volume integrals over Vd, the source volume V, and those 5

surface integrals over Sd and Sr in (3) which involve the incident field. The matrix [B8] results from the discretization of the surface integral over Sr into M, elements, i.e. that over the impedance or resistive sheets. Finally, the column {gS} gives the contribution of the boundary integral over So when the last is subdivided into N. surface elements. The computation of the matrix elements g' depends on the representation of H (or H""at) in terms of E (or Esat) on So, and this will be addressed later. Specific expressions of the matrix elements appearing in (5) in connection with the functional (2) are Ae = II [-(V x W).(V x W)kErW We] dv (6) CP = CP + Cip = af we [jkoZoJ+ V x(M') dv + jkZo Wj w. (H90 x n) ds. (7) and B = JkoZoJj(xW ) x ( W ds Bs = jkoZo WI (HSCt x n) ds. (8) In these, H90 represents the magnetic field on So when the antenna cavity is absent. Were we to consider the discretization of the functional in (3), the second integral in the definition of Cf in (7) must then be replaced by C = jkoZo [ WP hxH')ds+ J nx P) (nx E) ds + JJJ [~(V x W). (V x E')-keW Ei] dv. (9) in which VdP is the volume of the pth element inside the dielectric. For antenna parameter computations, it is of course understood that (E', H) are zero whereas for scattering computations (J, M') are zero. If the volume V contains lumped loads at the ith edge of the eth volume element, we must then add to the Ae matrix element the term - jkoZol eWe Ani = Z..| E[ W. dv (10) t3 ZLS v where I denotes the length of theen edge load, ZL is the value of the load in Ohms, and s represents its cross-sectional area. Furthermore, for antenna parameter computations in connection with a coaxial cable feed, we must add to the left side of (5) the term [37] [DP]{EP} + {PP} 6

where DP-jko (Wj- W i)ds (11) - ZO JJSb and pX I-ko jb,*i(ds. (12) cab P In these expressions, SPab denotes the pth surface element which belong to Scab as illustrated in Figure 1. The strength of the current in the center conductor of the cable is denoted by 1o and p denotes the radial distance measured from this conductor. Further, p is the associated unit vector along the radial direction away from the center conductor. Below we consider the specifics associated with the computation of the column {gs}. It will be seen that different hybrid finite element methods result depending on the approach used for the computation of {gS}. Finite element-boundary integral method The hybrid FE-BI method was introduced in the early seventies [50, 51] but its application to electromagnetic scattering did not occur until the eighties [52]-[55]. This method combines the adaptability of the FEM with the rigor of the BI or moment method. Because of the sparsity of the FEM matrix, this method is preferable to the volume moment method formulation [6], [7], which leads to completely full systems and is thus not attractive for iterative solution. The overall FE-BI system is partly sparse and partly full, but as will be seen later, the system can be efficiently solved with only O(N) memory demand on using uniform zoning on the mesh termination boundary which can be actually chosen to coincide with the target's surface/aperture [17], [28]-[30]. In accordance with the FE-BI method, the magnetic field in the integral over So is represented by an integral equation which may take the form [17], [31], [32] H = H~ + / G(r, r') [E(r') x ] ds' (13) where G is the dyadic Green's function of the second kind such that n x (V x G) = 0 on So. For the antenna problem in Figure 1, G becomes the free space dyadic Green's function _ 1 \ e-Ako Or-ri G =-j2koYo (i+ 2VV) (14) with r and r' representing the observation and integration points, respectively, and I xx + yy + zz is the unit dyad. The factor of 2 in (14) is due to image theory. In connection with this problem, i.e. that of a cavity recessed in a ground plane, Hg~ is equal to the sum of the incident and reflected fields for scattering computations, or zero for antenna parameter computations. If So is not planar, then as noted earlier, for scattering computation HI0 is equal to the field scattered by an otherwise metallic surface So. That is, if So is cylindrical then Hg~ is the total field due to a plane wave incident on a smooth metallic cylinder. 7

The integral in (13) can be discretized by introducing the expansion (for the s'th surface element) ns E(r')= Es W' (r') (15) i=l where W" (r') are the basis for the volume element which borders the s'th surface element belonging to the surface So. As usual, ns denotes the number of edges forming the surface elements comprising So. When (15) is introduced in (13) and then into (8), we obtain Nse ns g=- E E' G'^ (16) s'=l i=l where the coefficients Gijs are, of course, dependent on the specific form of the Green's function G(r,r'). For an aperture recessed in a ground plane, G is given by (14), and on employing the divergence theorem twice, in this case we may explicitly write [31], [32] G3 = 2k2 fJj [ x Wj(r)] {J, [W'(r') x ] Go(rr')ds'} ds (17) +2 If1 Vs [W (r) x z] {f' V'. [W (r') x z] Go(r, r) d'} ds when the normal to So is n = i. From (16), it is seen that the computation of the column {g} can be expressed in matrix form as {9g} = [1]{E) } (18) where the elements of the [9] matrix are given by G-8s once the local indices of the last are transformed to global indices. Clearly, from (17), all elements of the [9] matrix are non-zero, and in spite of the FE-BI method's rigor, this destroys the sparsity of the overall system resulting from (5). Moreover, the memory requirement for storing [9] is O(N2e), where Nse denotes the total number of edges used in the discretization of So. Unless special care is taken, a direct solution of (5) will be rather computationally intensive making the FE-BI method unsuitable for structures having electrically large So. As will be discussed later, if [9] is made block Toeplitz, its memory demand is reduced to O(Nse). Further, by using the fast Fourier transform (FFT) for computing {g3}, the number of operations in carrying out the matrix product in (18) is also reduced from O(N2s) to O(Nse log2 Nse) [57]. Finite Element-ABC Method In accordance with this method, some local boundary condition is enforced on the artificial surface So which encloses V. The goal with any ABC is to eliminate backward reflections from So, and a variety of these boundary conditions have been proposed over the years beginning with those of Bayliss and Turkel [58] and Enquist and Majda [59]. More recently, other ABCs such as those by Webb and Kanellopoulos [60] and Sesques [61] (see also Givoli [62]) have been proposed. All ABCs provide an approximate relation between the E and H fields 8

on the surface So which in most cases is derived by assuming a field expansion in inverse powers of r, where r is the radial distance from the center of So. If the ABC annihilates the first (2m + 1) inverse powers of r, it is then referred to as an mth order ABC. The zeroth order ABC is the Sommerfeld radiation condition -jkoZon x Hcat = -jko x x Escat (19) and a second order vector ABC proposed in [60] is given by -jkoZon x Ha = aE< + /V x [n (V x Ecat)] +./V, (v. Ecat) (20) In this a = jko, /3 = 1/(2jko + 2/r), and the subscripts t and n denote the tangential and normal components to So, respectively. These ABCs were derived for a spherical surface, but have been found to work well when So is made piecewise planar so that it better conforms to the body's surface [27]. In that case, the coefficient,/ reduces to /3 = 1/(2jko). A second order ABC which allows for arbitrary surface curvatures was recently proposed [61], and it is given by -jkoZo (n x Hscat) = 2 [{4m - Fc + D(jko l ) + (cl )2} Et, + i. V x {i(V x E)n} 2 2 2 2 + j (jko + 3Cm - - 2cl) {i t Vt(V E,)}]D AK - 2K1 ] (21) ko \Km 2 2 2 where AK = 1 - K2 D = 2jk0o+ 5m tKm In the above equation, t1, t2 are the principal surface directions and CI, KC2 are the corresponding curvatures, Km = (C1 + n2)/2 is the mean curvature and tog = c11c2 is the Gaussian curvature. Regardless of which ABC is chosen, we shall denote the right hand sides of (19)-(21) by the function P(Escat), allowing us to rewrite (8) as g = Jj Ws. P(Escat) (22) Then, upon introducing the expansion (15), we may express the column {gS} in matrix form as in (18). However, in this case all elements of [Q] vanish when s s', and for s = s we have that Ss sd o; f [zwa - w;(t + V (V x WS),. (V X WS), - -( w;t)(V Wjst)] ds 9

in connection with the ABC (20). As a result the [5] matrix is very sparse, and this is the main advantage of the FE-ABC method. Obviously, the surface So must be placed at some distance away from the body, whereas in the case of the FE-BI method, So can be coincident with the outer surface of the structure. However, as will be demonstrated later, we have found that for scattering computations So can be placed only 0.3 wavelengths away from the target when the FE-ABC method is implemented in connection with the ABC given by (20) [27]. Finite Element-Artificial Termination Structure Instead of enforcing an ABC on So, to eliminate backward reflections and ensure the outgoing nature of the waves, we may alternatively place in front of So some material absorber which serves the same purpose. This absorber may be of the same kind as that used in anechoic chambers, but to reduce its thickness to a minimum, it is instructive to design a thin multilayer absorber having, possibly, artificial constitutive parameters. One such metal-backed mathematical absorber is illustrated in Figure 4 [36]. The artificial termination structure (ATS) has the same advantages as the ABC in that it retains the sparsity of the overall system. Specifically, for metal-backed absorbers, the boundary integral over So is zero, and thus the column {gS} in (5) is eliminated. So far, our implementations have showed that the ATS will not necessarily reduce the volume between So and the target. Moreover, the discretization of the artificial dielectric layers may cause some additional inconvenience. Mesh Generation A first step in any numerical analysis is the generation of an appropriate discrete model of the structure. For our case, this amounts to subdividing the volume enclosed by the surface So or the cavity into bricks or tetrahedrals, a procedure referred to as mesh generation. For rectangular cavities and structures, the mesh generation using rectangular bricks is fairly simple. An approach which was employed is to subdivide the volume V into layers of thickness Az. Each layer is then readily subdivided into Ns x NI rectangular bricks. If the cavity housing the antenna structure is nonrectangular, rectangular bricks can still be used, but in that case staircasing will be unavoidable. However, for structures of arbitrary shape, it is necessary to resort to other elements such as tetrahedrals or curvilinear bricks. The generation of meshes using these elements requires sophisticated pre-processing packages and commercial software is available for this purpose. We have used a package marketed by SDRC I-DEAS [63] which employs a block meshing scheme and can be run on many different workstation platforms. A variety of other CAD packages are commercially available [64]. Regardless of which mesh generation scheme is employed, the final result is a file often called the Universal file. As a minimum, this file contains the following Tables: 1. A list of the global volume element Nos. and their material code. 2. A list of the nodes forming each volume element. 10

3. A list of the (x, y,z) coordinates for each node. 4. A list of global edge Nos. and a list of the pairs of global node Nos. forming an edge. The last table is not usually contained in the Universal file but can be readily constructed from the first three tables. To the above lists we must also add information pertaining to the antenna geometry, feeding lines, resistive cards, lumped loads and about the presence of a coax cable or a probe feed. This information is generated by a preprocessor based on information provided by the user. For example, in specifying the location of a rectangular patch or stripline, the user may enter the z location of the patch along with the (, y) coordinates for two opposite corners. For non-planar metallic surfaces, it is necessary to flag or group the appropriate surface elements bordering the metallic surface. This may be done when generating the mesh. A similar procedure is required for identifying the presence of resistive cards. The resistivity must be specified for each resistive surface element, and this is conveniently done by referring to a materials identification code. The specification of lumped loads or conducting posts is relatively simple. These are assumed to be placed at the edges of the volume elements and consequently the mesh must be generated with this in mind. Once the mesh is completed, the user specifies the location or edge number and the impedance value of each load or post. If the post is perfectly conducting, i.e. a short, it is flagged as such and during matrix assembly, the field coefficient associated with the edge occupied by the short is set to zero. The center conductor of the coaxial cable is defined in this manner. Based on the structural and loading specifications, the following additional Tables may also be generated: * List of edge Nos. on a perfectly conducting surface (the fields of these edges are set to zero a priori when employing the total field formulation). * List of volume elements along with their global edge Nos. which border the outer surface. * List of volume elements and global edge Nos. bordering a resistive surface along with the associated resistivity. * List of global edge Nos. along with subglobal edge Nos. for combinations. The preprocessor can be run interactively, with or without a graphical user interface, depending on which mesh generation package is employed. The matrix elements are computed directly from the data provided by these tables using the algorithms outlined in the previous section. Matrix Assembly and System Solution Once all input tables are read, the goal is to assemble the matrix system f4 +E [} [l [0] [0J 1 { } =_ l J {En L [] [5 {n I} - () 11

where [A] is generally a sparse square symmetric matrix, [9] was discussed earlier, the column {E'} denotes the field values at the edges not bordering So, and {EB} are the corresponding coefficients for the edges on So. The columns {bV} and {b } are computed from {CP} as defined in (7) and (11). All matrix coefficients in the above system are numbered in accordance with the unique global or their own subglobal/local numbering schemes. If the pth global edge lies on a metallic surface or coincides with a metallic wire or post, Ev or EP are set to zero or equal to the negative of the tangential incident field, depending on which formulation is employed. Consequently the computation of Apn is skipped. Note that with this notation, mn and n refer to global edge numbers, whereas the subscripts i and j used earlier are local indices for the eth element. Obviously, an edge may belong to more than one volume element and in that case A is constructed by summing the element equations for that edge. It is important to note that only the non-zero elements of the [A] matrix are stored, thus, ensuring the O(N) memory demand of the system. On the average, for tetrahedrals the minimum number of non-zero elements per row is 9 and the maximum number of nonzeros per row is 30. However, these numbers vary among geometries. In our implementations, the matrix was generated using a modified ITPACK storage scheme [65] and then stored as a long vector. This resulted in a storage requirement of approximately 15Nv to 16Nv, where Nv is the number of edges in V or unknowns. The ITPACK storage scheme is attractive for generating finite element matrices since the number of comparisons required while augmenting the matrix depends only on the locality of the corresponding edge and not on the number of unknowns. However, on using the ITPACK storage scheme for our application, almost half the space is lost in storing zeros. The modified ITPACK scheme [66] does alleviate this problem to a certain degree; however, 30% of the space is still losto in zero padding. The best trade-off between storage and speed for our application is obtained by storing the non-zero matrix elements in a long complex vector, the column indices in a long integer vector and the number of non-zeros per row in another integer vector. This data structure is referred to as the Compressed Sparse Row (CSR) format. In our implementation, a map of the number of non-zeros for each row is obtained through the modified ITPACK algorithm. The main program stores the matrix in CSR format, thus minimizing storage and sacrificing a bit of speed. The required storage is 15N complex words plus two integer vectors of length 15N and N, respectively, which serve as pointers. For the FE-ABC method, the [9] matrix is also sparse and is generally combined with the [A] matrix in solving (23). However, as noted earlier, [9] is a full complex matrix when the FE-BI method is employed. In this case, if we were to solve (23) using a standard algorithm, it is appropriate to first rewrite it in the form [AV]{E } + [AI ]{EB} = {bv} (24) [A]{En }+ [-]E} = {b } (25) where [AB II] and [A"1] are sparse matrices since they were decomposed from [A], and generally [AB] = [A4B]T. Also, the superscripts B and V imply fields/matrices associated with the boundary surface So or the volume V. The first of these can now be substituted into the second to obtain the system [52] {[g] - [AB][Av] —I[AB]} {EB } = {bB} - [ABI[Av]-{bV} (26) ~ L i LI 1 1 I?I \U 1 1 - L.MI 12

where [AV] -[AB] is actually computed by solving the system [A;{1I}1 = [AB]. The fact that [AV] is a symmetric, sparse matrix is an important advantage in solving (26). We note that a substitution of (25) into (24) would require the inversion of a full complex matrix. Any of the many existing direct and iterative solvers [67] can be used for solving (26) including LU decomposition, banded matrix algorithm, and Lanczos algorithm. The last is better suited for large systems. Nevertheless, for a general full matrix, whichever solver is used will demand excessive computational intensity. To alleviate this, it is important to make [S] block Toeplitz by employing uniform gridding on the surface So. This is quite practical for flat surfaces (see Figure 1), and although it may cause some extra effort in creating the mesh and possibly staircasing, it is well worth it. Thus, if the volume element is a rectangular brick, then all elements on So should be chosen to be rectangular and equal. In the case of tetrahedrals, we must choose all surface elements to be right triangles of equal size. If [9] is block Toeplitz, it can then be recast in a block circulant form [57], [68], [69], and in that case the FFT can be used to compute the product [9]{E } in conjunction with an iterative solver. Specifically, if we denote a Nitp x Mtp circulant block of [9] by [gtp], we may then write (nm and n denote the x and y indices defining each element/edge) [57], [70]-[72] Altp-1 Ntp-1 E E GtmmnnEn = FFT2 {FFT(GtP) FFT(E } (7) m'=0 n'=O This operation requires that only one row of the [9tp] matrix block be stored reducing the storage of the overall [9] matrix to O(Nse). WXChen the biconjugate gradient (BCG) algorithm [73] for symmetric systems is used for the solution of (23), the overall storage requirement for rectangular elements is about 4Nv + 8Nse, excluding the storage needed for [A] and [9]. For triangular surface elements the BCG algorithm's storage requirement is about 4.5 times larger. The elements of [A] are extremely simple to compute and can be actually generated on the fly from the element submatrices for storage economy at the expense of computational efficiency. We have observed rather impressive convergence rates using this algorithm (with diagonal preconditioning) for FE-BI systems having more than 100,000 unknowns. For example, for 180,000 unknowns, convergence is achieved in 57 iterations. Another system of 25,000 unknowns required 66 iterations per angle with an average CPU time of 2 sec per iteration. Generally, the sampling rate was 12 to 15 elements per linear wavelength. The BCG algorithm was also used for solving the sparse system (23) in connection with the FE-ABC formulation. This formulation was exclusively used for treating targets associated with non-planar termination surfaces Co. Rather large-scale FE-ABC systems (over 250,000 unknowns) were solved via the symmetric BCG algorithm with rather impressive convergence rates. Specifically, convergence was generally achieved with N/100 to N/50 iterations. The conjugate gradient squared (CGS) algorithm is generally faster than BCG but is more unstable since the residual polynomials are merely the squared BCG polynomials. Also, the CGS algorithm fails to exploit the symmetry of the matrix system. The performance of the BCG algorithm in solving the FE-ABC sparse system (23) was examined on the vector Cray YMP. It was observed that the vector updates and dot products reached speeds of about 190 MFLOPS. However, the matrix vector products, which involve indirect addressing, ran at about 45.5 MFLOPS on a single processor of the 8-processor Cray YMP. As a rule of thumb, the BCG algorithm on the Cray YMP consumes about 13

4.06 microseconds/iteration/unknown. A simple diagonal preconditioning can improve the convergence rate by about 35%, but a traditional ILU preconditioner was found less attractive because its convergence improvement was in most cases cancelled by the additional time required for its implementation. The parallelization of the BCG solver was also carried out on a 28-processor KSR1 parallel machine [66]. Figure 5 shows the speed-up for three versions of the sparse solver with diagonal preconditioning as the number of processors is increased from 1 to 28. The maximum speed-up is observed to be about 19 and thus the solver was about 3.25 times faster on the 28-processor KSR1 than on the single-processor Cray YMP. Results In this section we present radar cross section (RCS) computations of cavity-backed antenna configurations, arrays (in a ground plane) and free-standing inhomogeneous structures. Once the fields in the computational domain are obtained, only the tangential fields are needed over a closed surface or the aperture So for an evaluation of the RCS or radiated fields. For free-standing structures, the far zone scattered field was obtained by applying the StrattonChu integral equation over the surface bounding the target. Example 1: RCS of a single patch residing in a rectangular cavity (FE-BI method) The geometry of this configuration is shown in Figure 6 along with the computed backscatter ao0 RCS pattern. The patch is 1.448" x 1.038" and resides on a dielectric substrate having thickness t = 0.057" and relative permittivity er m 4.0. The substrate is housed in a 2.89" x 2.10" x 0.057" rectangular cavity recessed in a ground plane. The excitation was a 0 -polarized plane wave at 9.2 GHz in the xz plane, and the analysis was done using rectangular bricks. As seen the computed ato pattern is in good agreement with the measured data [33]. We note that for this example the cavity terminations play a more important role because the computation frequency is not near the resonant frequency of the patch. The measured data were collected by placing the cavity-backed antenna configuration on a 5 foot long low cross section body [74]. Example 2: RCS of a circular patch near resonance (FE-BI method) The geometry of this patch configuration is shown in Figure 7 along with the calculated ass backscatter RCS as a function of frequency. The diameter of the circular patch is 1.3 cm and resides on a substrate 0.0787 cm thick, having 6cr 2.33. The excitation is a plane wave incoming at (0i = 60~, Oi = 1800), and the computations were done using tetrahedral elements with the substrate terminated at a diameter of 1.6148 cm. However, because the calculations are near the patch resonance ( — 8 GHz), the substrate terminations do not contribute significantly to the RCS of the antenna. As seen, our computations are in good agreement with the measured data [75] which were collected with the patch placed in a cavity having a different periphery. 14

Example 3: RCS of a patch with distributed resistive loading (FE-BI method) With the motivation of reducing the patch RCS at frequencies outside its operational band without substantial compromise in gain, we consider the placement of a narrow resistive ribbon (R-skirt) around the periphery of the patch as illustrated in Figure 8. Because the ribbon's electrical width increases with increasing frequency, it is expected to have a rather noticeable effect on the R.CS of the patch at higher frequencies. Thus, greater RCS reduction can be attained at higher frequencies without compromising the antenna's gain performance. If a certain level of RCS reduction is desired at resonance, lumped loads can still be used in conjunction with the R-skirt loading. The backscatter aes RCS data given in Figure 8 with and without loading correspond to a plane wave excitation incident in the xz plane and at an angle of 0 = 70 with respect to the z axis. However, as shown in Figure 9, the RCS reduction is fairly uniform over all 0 angles. As seen, the resistive skirt loading reduces the RCS in and out of the band as opposed to the narrowband performance of the lumped loads. Not surprisingly, the placement of the R-skirt is seen to slightly change the resonant frequency of the patch from 1.9522 GHz down to 1.744 GHz. From the inset in Figure 8, the RCS reduction at resonance is about 9 dB for all angles of incidences and thus the corresponding gain reduction is approximately 4.5 dB. Example 4: Scattering by an array of 3 x 30 slots (FE-BI method) Because of the low O(N) memory demand of the FE-BI method when combined with the FFT, one can consider the case of large finite arrays. The plot in Figure 10 is the oa~, backscatter RCS pattern of a 3 x 30 array of slots in a ground plane. Each slot has a rectangular 0.094" x 0.492" aperture and is 0.2059" deep. More than 100,000 unknowns were used for the simulation of this array and as seen the computed RCS pattern is in good agreement with corresponding moment method data. A characteristic of this pattern is the larger lobe at grazing which is often undesirable. As shown, this is completely suppressed when the last three slots from each end of the array are covered with a resistive card. The values of the resistive cards placed over the last three slots were 37.7 Q, 188.5 Q and 94.25 Q. Example 5: RCS of a composite cube (FE-ABC method) The geometry to be considered is an inhomogeneous cube 0.5A per side. Half of the cube (0.5A x 0.5A x 0.25A) consists of a dielectric having r, = 2 -j2 with its surfaces bordered by resistive cards having a resistivity equal to R = Zo. The other half of the cube is metallic. In Figure 11, we compare the computed principal plane ao, RCS pattern with corresponding moment method data [76]. For the FE-ABC solution the scatterer was enclosed within a cubical outer boundary placed only 0.3A away from the scatterer. This resulted in a 30,000 unknown system which converged in about 400 iterations on the Cray YMP when using the Sommerfeld radiation condition and in 1600 iterations when the second order ABC was used. For this geometry, the second order ABC did not provide a significant improvement in accuracy over the first order condition. 15

Example 5: RCS of a 1.5A x 1A x 1A cavity (FE-ABC method) The parallelized version of the FE-ABC code was used to compute the RCS of 1.5A x 1A x 1A perfectly conducting inlet. In generating the mesh, this geometry was enclosed in a sphere of radius 1.35A. The resulting system had 224,476 unknowns and converged in about 5,000 iterations on the KSR1. The computed backscatter pattern for the vertical polarization is shown in Figure 12. As seen, the computations are in good agreement with the measured data [77]. Concluding Remarks The 1980's brought to focus a variety of computational methods for electromagnetic applications, including scattering, antennas, frequency selective surfaces, conformal arrays, microwave circuits and VLSI packaging. Not a single method is likely to predominate in all of these areas but it is nonetheless certain that future developments in computational electromagnetics will take into account the new architectures of high performance computing facilities. The established integral equation methods will likely continue to be used for small and medium scale computations and optimized on high performance computing platforms. However, growing requirements to develop codes for modeling complex and large non-metallic geometrical structures; the availability of massively parallel computing architecturese and gallops in computer speeds do not seem to favor integral equation methods. For frequency domain computations, the finite element method, because of its geometrical adaptability, low O(N) memory demand and suitability for parallelization is undoubtedly a more promising computational method. The presented examples clearly demonstrated its accuracy for large scale electromagnetic computations but, needless to mention, there is much to be accomplished before the technique can be efficiently employed in modeling structures involving several millions of grid points. For example, there are needs for efficient mesh generation algorithms, improved convergence solvers on parallel architectures, robust mesh truncation schemes, techniques for handling complex geometrical details and thin films, new adaptable element shapes and possibly new finite element formulations whose objective is to converge on a certain performance parameter (such as the RCS) rather than the fields themselves. 16

References [1] J.H. Richmond, "Radiation and scattering by thin wire structures in a homogeneous conducting medium," IEEE Trans. Antennas Propagat., pp. 124-126, Jan. 1974. [2] G.J. Burke and A.J. Pogio, "Numerical Electromagnetics Code (NEC-1), Parts I, II, III: NEC Program Description," Naval Ocean Systems Center, San Diego CA, NOSC TD 116, 1980 (original version was published in 1977). [3] E.F. Knott and T.B.A. Senior, "Non-specular cross section study," Air Force Avionics Laboratory Technical Report AFAL-TR-73-422, Jan. 1974 (also University of Michigan Radiation Laboratory Report 011764-1-T). [4] E.H. Newman, P. Alexandropoulos and E.K. Walton, "Polygonal plate modeling of realistic structures," IEEE Trans. Antennas Propagat., pp. 742-747, July 1984. [5] S.M. Rao, D.R. Wilton and A.W. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," IEEE Trans. Antennas Propagat., pp. 409-418, May 1982. [6] D.H. Schaubert, D.R. WVilton and A.W. Glisson, "A tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies," IEEE Trans. Antennas Propagat., pp. 77-85, Jan. 1984. [7] C. Tsai, H. Massoudi, C.H. Durney and M.F. Iskander, "A procedure for calculating fields inside arbitrarily shaped inhomogeneous dielectric bodies using linear basis functions with the moment method," IEEE Trans. Microwave Theory Techn., Vol. 34, pp. 1131-1139, Nov. 1986. [8] K.K. Mei, "Unimoment method of solving antenna and scattering problems," IEEE Trans. Antennas Propagat., Vol. 22, pp. 760-766, Nov. 1974. [9] K. Umashankar, A. Taflove and S.M. Rao, "Electromagnetic scattering by arbitrarily shaped three-dimensional homogeneous lossy dielectric objects," IEEE Trans. Antennas Propagat., Vol. 34, pp. 758-766, June 1986. [10] P.L. Huddleston, L.N. Medgyesi-Mitschang and J.M. Putnam, "Combined field integral equation formulation for scattering by dielectrically coated conducting bodies," IEEE Trans. Antennas Propagat., Vol. 34, pp. 510-520, April 1986. [11] A. Taflove and K.R. Umashankar, "The finite difference time-domain (FD-TD) method for electromagnetic scattering and interaction problems," J. Electromagnetic Waves Appl., Vol. 1, pp. 243-267, 1987. [12] N.K. Madsen and R.W. Ziolkowski, "Numerical solution of Maxwell's equations in the time domain using irregular nonorthogonal grids," Wave Motion, Vol. 10, pp. 583-596, 1988. [13] V. Shankar, W.F. Hall and A.H. Mohammadian, "A time domain differential solver for electromagnetic scattering problems," Proc. IEEE, Vol. 77, pp. 709-721, May 1989. [14] R.D. Graglia, P.L.E. Uslenghi and R.S. Zick, "Moment method with isoparametric elements for three-dimensional anisotropic scatterers," Proc. IEEE, Vol. 77, pp. 750-760, 1989. [15] J. D'Angelo and I.D. Mayergoyz, "Finite element methods for the solution of RF radiation and scattering problems," Electromagnetics, Vol. 10, pp. 177-199, 1990. 17

[16] K.S. Kunz and R.J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics, CRC Press, 1993. [17] J.M. Jin, J.L. Volakis and J.D. Collins, "A finite element-boundary integral method for scattering and radiation by two- and three-dimensional structures," IEEE Antennas Propagat. Magazine, Vol. 33, No. 3, pp. 22-32, June 1991. [18] A.F. Peterson, "Analysis of heterogeneous electromagnetic scatterers: research progress of the past decade," Proc. IEEE, Vol. 79, No. 10, pp. 1431-1441, Oct. 1991. [19] M.A. Fusco, M.V. Smith and L.W. Gordon, "A three-dimensional FDTD algorithm in curvilinear coordinates," IEEE Trans. Antennas Propagat., Vol. 39, pp. 1463-1471, Oct. 1991. [20] C. Wu, K.-L. Wu, Zhi-Quiang Bi and J. Litva, "Accurate characterization of planar printed antennas using finite-difference time-domain method," IEEE Trans. Antennas Propagat., Vol. 40, pp. 526-534, May 1992. [21] D.M. Sheen, S.M. Ali, M.D. Abruzahra and J.A. Kong, "Application of three-dimensional finite-difference time-domain method to the analysis of planar microstrip circuits," IEEE Trans. Microwave Theory Techn., Vol. 38, No. 7, pp. 849-858, July 1990. [22] P.P. Silvester and R.L. Ferrari, Finite Elements for Electrical Engineers, 2nd ed., Cambridge, U.K.: Cambridge Univ. Press, 1990. [23] J.W. Parker, R.D. Ferraro and P.C. Liewer, "Distributed memory tangential finite elements for inhomogeneous scatterers in free space," 1992 URSI Radio Science Meeting, Chicago IL, URSI Digest, p. 516. [24] A. Chatterjee, J.M. Jin and J.L. Volakis, "Application of edge-based finite elements and vector ABCs in 3D scattering," 1992 IEEE AP-S Symposium, Digest pp. 528-531. [25] D.D. Paulsen, W.E. Boyse and D.R. Lynch, "Continuous potential Maxwell's equations on node-based finite elements," IEEE Trans. Antennas Propagat., Vol. 40, pp. 1192-1200, Oct. 1992. [26] M.A. Morgan, ed., Finite Element and Finite Difference Methods in Electromagnetic Scattering, Elsevier: New York, 1990. [27] A. Chatterjee, J.M. Jin and J.L. Volakis, "A finite element formulation with absorbing boundary conditions for three-dimensional scattering," IEEE Trans. Antennas Propagat., Feb. 1993. [28] J.-M. Jin and J.L. Volakis, "TM scattering by an inhomogeneously filled aperture in a thick ground plane," IEE Proceedings, Part H, Vol. 137, No. 3, pp. 153-159, June 1990. [29] J.-M. Jin and J.L. Volakis, "TE scattering by an inhomogeneously filled aperture in a thick conducting plane," IEEE Trans. Antennas Propagat., Vol. AP-38, pp. 1280-1286, August 1990. [30] J.D. Collins, J.-M. Jin and J.L. Volakis, "A combined finite element-boundary element formulation for solution of two-dimensional problems via CGFFT," Electromagnetics, Vol. 10, pp. 423-437, 1990. [31] J.M. Jin and J.L. Volakis, "A finite element-boundary integral formulation for scattering by three-dimensional cavity-backed apertures," IEEE Trans. Antennas Propagat., Vol. AP-39, pp. 97-104, January 1991. 18

[32] J.M. Jin and J.L. Volakis, "Electromagnetic scattering by and transmission through a threedimensional slot in a thick conducting plane," IEEE Trans. Antennas Propagat., Vol. AP-39, pp. 543-550, April 1991. [33] J.M. Jin and J.L. Volakis, "Scattering and radiation from microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas Propagat., Vol. AP-39, pp. 1598-1604, Nov. 1991. [34] J.M. Jin, J.L. Volakis, C.L. Yu and A. Woo, "Modeling of resistive sheets in the context of finite element solutions," IEEE Trans. Antennas Propagat., Vol. 40, pp. 727-731, June 1992. [35] A. Chatterjee, J.M. Jin and J.L. Volakis, "Computation of cavity resonances using edge-based finite elements," IEEE Trans. Microwave Theory Techn., Nov. 1992. [36] J.M. Jin, J.L. Volakis and V.V. Liepa, "A fictitious absorber for truncating finite element meshes in scattering," IEE Proceedings, Part H, Vol. 139, No. 5, pp. 472-476, Oct. 1992. [37] J. Gong, J.L. Volakis, A. Chatterjee and J.M. Jin, "Characterization of cavity-backed conformal antennas and arrays using a hybrid finite element method with tetrahedral elements," 1992 IEEE Antennas Propagat. Symposium, Chicago IL, Digest pp. 1629-1631. [38] J.L. Volakis and J.M. Jin, "A technique to substantially lower the resonant frequency of the microstrip patch antenna," IEEE Trans. Microwave Theory Techn. Letters, Vol. 2, No. 7, pp. 292-293, July 1992. [39] T.B.A. Senior, "Combined resistive and conductive sheets," IEEE Trans. Antennas Propagat., Vol. AP-33, pp. 577-579, 1985. [40] O.M. Ramahi and R. Mittra, "Finite element analysis of dielectric scatterers using the absorbing boundary conditions," IEEE Trans. Magnetics, Vol. 25, No. 4, pp. 3043-3045, July 1989. [41] A. Bossavit and I. Mayergoyz, "Edge-elements for scattering problems," IEEE Trans. Magnetics, Vol. 25, No. 4, pp. 2816-2821, July 1989. [42] 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-21, 1987. [43] H. Whitney, Geometric Integration Theory, Princeton Univ. Press, 1957. [44] B.M. Azizur Rahman and J. Brian Davies, "Penalty function improvement of waveguide solution by finite elements," IEEE Trans. Microwave Theory Techn., Vol. MTT-32, No. 8, Aug. 1984. [45] D. Ross and J.L. Volakis, "Numerical divergenceless shape functions for three-dimensional finite element analysis," 1993 IEEE AP-S Symposium and URSI Radio Science Meeting, Ann Arbor MI. [46] L.C. Kempel and J.L. Volakis, "A finite element-boundary integral method for cavities in a circular cylinder," 1993 IEEE AP-S Symposium and URSI Meeting, Ann Arbor MI. [47] G.E Antilla and N.G. Alexopoulos, "Radiation and scattering from complex 3D curvilinear geometries using a hybrid finite element-integral equation method," 1992 IEEE Antennas Propagat. Symposium, Chicago IL, Digest pp. 1758-1761. 19

[48] C.W. Crowley, P.P. Silvester and H. Hurwitz, Jr., "Covariant projection elements for 3D vector field problems," IEEE Trans. Magnetics, Vol. 24, No. 1, pp. 397-400, Jan. 1988. [49] A. Bossavit, "A rationale for edge-elements in 3D fields computations," IEEE Trans. Magnetics, Vol. 24, No. 1, pp. 74-79, Jan. 1988. [50] P. Silvester and M.S. Hsieh, "Finite element solution of 2-dimensional exterior field problems," Proc. Inst. Elec. Eng., Vol. 118, pp. 1743-1747, Dec. 1971. [51] B.H. McDonald and A. Wexler, "Finite element solution of unbounded field problems," IEEE Trans. Microwave Theory Techn., Vol. MTT-20, pp. 841-847, Dec. 1972. [52] J.M. Jin and V.V. Liepa, "A note on hybrid finite element method for solving scattering problems," IEEE Trans. Antennas Propagat., Vol. 36, No. 10, pp. 1486-1490, Oct. 1988. [53] J.M. Jin and J.L. Volakis, "A new technique for characterizing diffraction by inhomogeneously filled slots of arbitrary cross section in a thick conducting plane," IEE Electronics Letters, Vol. 25, No. 17, pp. 1121-1122, 17 Aug. 1989. [54] K.L. Wu, G.Y. Delisle, D.G. Fang and M. Lecours, "Coupled finite element and boundary element methods in electromagnetic scattering," Chapter 3 in Finite Element and Finite Difference Method in Electromagnetic Scattering, ed. M.A. Morgan, Elsevier: New York, 1990. [55] X. Yuan, D.R. Lynch and J.W. Strohbehn, "Coupling of finite element and moment methods for electromagnetic scattering from inhomogeneous objects," IEEE Trans. Antennas Propagat., Vol. 38, pp. 386-393, March 1990. [56] S.K. Jeng, "Scattering by a cavity-backed slit in a ground plane; TM case," IEEE Trans. Antennas Propagat., Vol. 39, pp. 661-663, May 1991. [57] K. Barkeshli and J.L. Volakis, "On the implementation and accuracy of the conjugate gradient FFT method," IEEE Antennas Propagat. Magazine, Vol. 32, No. 2, pp. 20-26, April 1990. [58] A. Bayliss and E. Turkel, "Radiation boundary conditions for wave-like equations," Comm. Pure Appl. Math., Vol. 33, pp. 707-725, 1980. [59] B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves," Math. Comp., Vol. 31, pp. 629-651, 1977. [60] J.P. Webb and V.N. Kanellopoulos, "Absorbing boundary conditions for finite element solution of the vector wave equation," Microwave and Opt. Techn. Letters, Vol. 2, No. 10, pp. 370-372, Oct. 1989. [61] A. Chatterjee and J.L. Volakis, "Conformal absorbing boundary conditions for the vector wave equation," submitted to Microwave and Opt. Techn. Letters. [62] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier: New York, 1992. [63] SDRC I-DEAS, Structural Dynamics Research Corp., Milford 011, USA. [64] J.-C. Sabonnadiere and J.-L. Coulomb, Finite Element Methods in CAD: Electric and Magnetic Fields, Springer-Verlag: New York, 1987 (English trans. by S. Salon). 20

[65] D.R. Rincaid and T.C. Oppe, "ITPACK on supercomputers," Numerical Methods, Lecture Notes in Math., pp. 151-161, Springer-Verlag: Berlin, 1982. [66] A. Chatterjee, J.L. Volakis and D. Windheiser, "On the optimization of finite element code for 3D scattering computation," 9th Annual Review of Progress in Applied Computational Electromagnetics (ACES), Monterey CA, March 1993. [67] G.H. Golub and C.F. Van Loan, AMatrix Computations, Johns Hopkins Univ. Press: Baltimore MD, 1985. [68] J.L. Volakis and K. Barkeshli, "Applications of the conjugate gradient FFT method to radiation and scattering," in Applications of the Conjugate Gradient Method to Electromagnetics and Signal Analysis, ed. T.K. Sarkar, Elsevier: New York, 1991. [69] A.F. Peterson, S.L. Ray, C.H. Chan and R. Mittra, "Numerical implementations of the conjugate gradient method and the CG-FFT for electromagnetic scattering," in Applications of the Conjugate Gradient AMethod to Electromagnetics and Signal Analysis, ed. T.K. Sarkar, Elsevier: New York, 1991. [70] J.M. Jin and J.L. Volakis, "Biconjugate gradient FFT solution for scattering by planar plates," Electromagnetics, Vol. 12, pp. 105-119, 1992. [71] M.F. Catedra, E. Gago and L. Nufio, "A numerical scheme to obtain the RCS of threedimensional bodies of resonant size using the conjugate gradient method and the Fast Fourier Transform," IEEE Trans. Antennas Propagat., Vol. 39, pp. 528-537, May 1989. [72] P. Zwa.mborn and P.M. van den Berg, "The three-dimensional weak form of the conjugate gradient FFT method for solving scattering problems," IEEE Trans. Microwave Theory Techn., Vol. 40, pp. 1757-1766, Sept. 1992. [73] T.K. Sarkar, "On the application of the generalized biconjugate gradient method," J. Electromagnetic Waves Appl., Vol. 1, No. 3, pp. 223-242, 1987. [74] J.M. Jin and J.L. Volakis, "Radiation and scattering from microstrip patch antennas via. a hybrid finite element method," Int. Conference on Electromagnetics in Aerospace Appl., Politecnico di Torino, Italy, Sept. 1991, Conference Proceedings pp. 195-198. [75] D.G. Shively, M.C. Bailey, C.R. Cockrell and M.D. Deshpande, "Scattering from microstrip patch antennas using subdomain basis functions," Electromagnetics, 1993 (to appear). [76] Courtesy of Northrop Corp., B2 Division, Pico Rivera CA. [77] A. Woo, M. Schuh, M. Simon, T.G. Wang and M.L. Sanders, "Radar cross-section measurement data of a simple rectangular cavity," Technical Report NWC TM7132, Naval Weapons Center, China Lake CA, Dec. 1991. 21

List of Figures Figure 1. Geometry of a cavity-backed antenna/array recessed in a ground plane. Figure 2. Illustration of a scatterer enclosed by an artificial surface So. Figure 3. Geometries of various finite element shapes. Figure 4. Scatterer enclosed by an artificial absorber termination. (a) Illustration of the enclose. (b) A specific ATS. Figure 5. Speedup curves for the linear equation solver on the KSR1. Figure 6. Comparison of measured and calculated RCS patterns at 9 GHz for the shown rectangular patch. The patch is on a substrate of thickness 0.057 in having Er = 4. Figure 7. Comparison of measured and calculated backscatter RCS data for the shown circular patch as a function of frequency (Oi = 60~, Oi = 180~). Figure 8. aee backscatter RCS of the shown patch as a function of frequency with and without resistive skirt loading; incident wave direction Oi = 70~, bi = 180~. The patch is on a substrate of thickness 0.17558 cm having Er = 2.17. Figure 9. ao8 backscatter RCS of the patch in Figure 8 as a function of angle with and without resistive loading. Figure 10. Backscatter aee RCS pattern of a 3 x 30 slot array at 9 GHz. Pattern is taken in the xz plane. The slot centers are separated by 0.531" in each direction, are 0.2059" deep and have a rectangular 0.094" x 0.492" aperture. Figure 11. RCS pattern in the xz plane for the shown geometry. The solid curve is the FE-ABC pattern and the black dots are MoM data [76] for the ElnC = 0 polarization. The specific geometry is a cube (a = b = 0.5A). Its lower half section is metallic and the upper half section is composed of dielectric (Er = 2 - j2, p, = 1) bordered by a resistive card having R= Zo. 22

Figure 12. aee backscatter RCS pattern of a perfectly conducting rectangular inlet (1A x 1A x 1.5A) taken in the yz plane. The black dots indicate data obtained from the KSR1 for the FE-ABC code and the solid line represents measured data [77]. 23

Cavity (a) Plan view Patches / Region I stripline W Scab (b) Cross-sectional (xz) view Figure 1. Geometry of a cavity-backed antenna/array recessed in a ground plane. Figure 2. Illustration of a scatterer enclosed by an artificial surface So. 24

Right Angled Brick Skewed Brick I' dI ~ I Skewed SkewedTriangular Prism Triangular Prism Curvilinear Brick Tetrahedron Cylindrical shell Figure 3. Geometries of various finite element shapes. 25

15 20 10 -5 -15-... 0 30' 60 90 Observation Angle 6. deg. Figure 12. aso backscatter RCS pattern of a perfectly conducting rectangular inlet (1A x 1A x 1.5A) taken in the yz plane. The black dots indicate data obtained from the KSR1 for the FE-ABC code and the solid line represents measured data [77]. 34

A n SO Er8 ~artificial absorber X7 ~ ~ ~ ~ ~~-.Vd 0000ggg:g-X........E...= =1-j2.7 TIIl X-::gr=l1 -j2.7............ metal (a) (b) Figure 4. Scatterer enclosed by an artificial absorber termination. (a) Illustration of the enclose. (b) A specific ATS. 26

30 O linear Origina ---- 25 Postst re — B ---Prifet ch ----- 105 -..10 15 Pro es sors............. Fg10. S p curvs for te l r e n s r on te Kx 0 5 10 15 20 25 30 Processors Figure 5. Speedup curves for the linear equation solver on the KSR1. 27

y c ci -w X T = 0.057 in r = 4 f - 9.2 GHz -o — 2.89 in -10. -20. 1z) "E tz) with patch -30. -40. -50. 0. 15. 30. 45. 60. 75. 90. 0 (degrees) Figure 6. Comparison of measured and calculated RCS patterns at 9 GHz for the shown rectangular patch. The patch is on a substrate of thickness 0.057 in having Er = 4. 28

-10.0 E "C: C/) bjD a) CA 4r4 -20.0 -30.0 -40.0 -50.0 -60.0 -70.0 Er =2.33 t = 0.0787 cm measured | I EQ calculated t -80.0 7..00 7.20 7.40 7.60 7.80 8.00 8.20 8.40 8.60 8.80 9.00 frequency (in GHz) Figure 7. Comparison of measured and calculated backscatter RCS data for the shown circular patch as a function of frequency (Oi = 60~, q>i = 180~). 29

9375cm y(0, (0,0) 937Scm Resistive skirt 0.15625cm wide, R=0 Lower Left coer of patch askirt y)0.15625cm wide, R9675 Lower Left corner of patch at (xy)=(2.1875cm,2.96875cm) Patch dimensions 5.00cm x 3.4375cm Feed at (x,y) = (3.4375cm,3.75cm) -10.0 -20.0 *E V).C CO) (9 -30.0 -40.0 -50.0 -60.0 L 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00 Frequency in GHz Figure 8. aee backscatter RCS of the shown patch as a function of frequency with and without resistive skirt loading; incident wave direction 0i = 700, =i = 180~. The patch is on a substrate of thickness 0.17558 cm having e, = 2.17. 30

0.000 -10.0 -20.0 U) -30.0 - - No Rskirt at resonant f=1.9522GHz.-..-... =0 With Rskirt at resonant f=1.744GHz -40.0 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00 0 in degrees Figure 9. aeo backscatter RCS of the patch in Figure 8 as a function of angle with and without resistive loading. 31

y 0....... I....l r. l x30 Empty Slot Array -10 \1 -------- With Resistive Loading on the Last Three Elements: A Moment Method -20 -30 1 * ar~ ~ JIt- t * -II 0 10 20 30 40 50 60 70 80 90 Angle [deg] Figure 10. Backscatter aco RCS pattern of a 3 x 30 slot array at 9 GHz. Pattern is taken in the xz plane. The slot centers are separated by 0.531" in each direction, are 0.2059" deep and have a rectangular 0.094" x 0.492" aperture. -60.. 32

UNIV Y MICHIAN 3 9015 03527 4961 I Resistive Surfaces ^-.-^ —. ---...i. -- -::.:...:,[ - 1:I: 1...................... (a + 0.6)| 2 ' -I - - - - _1 MeshI Mea, 's Termination Metal i..-,;I DUllunuary 10 m 0 OE,< b -10 -20 0 30 60 90 120 150 Observation Angle 0, deg. 180 Figure 11. RCS pattern in the xz plane for the shown geometry. The solid curve is the FE-ABC pattern and the black dots are MoM data [76] for the E'nC = 0 polarization. The specific geometry is a cube (a = b = 0.5A). Its lower half section is metallic and the upper half section is composed of dielectric (c, = 2 - j2, ur = 1) bordered by a resistive card having R= Zo. 33