033762-1-T FEMATS GRAPHICAL USERS INTERFACE AND MESH TRUNCATION SCHEMES Naval Air Warfare Center China Lake, CA 93555-6001 August 1996 33762-1-T = RL-2464

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: REPORT No.: DATE SPONSOR: Further Development of FEMATS, Including Prismatic Meshes, Graphical Interface and New Mesh Truncations Schemes FEMATS graphical users interface and performance of mesh truncation schemes 033762-1-T August 1996 Dr. Helen Wang Code 472210D, Bldg 01400 Naval Air Warfare Center Weapons Div. China Lake, CA 93555-6001 phone: (619) 939-3931 SPONSOR P.O. or CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: CONTRIBUTORS TO THIS REPORT: N68936-95-M-E248 John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave Ann Arbor, MI 48109-2122 Phone: (313) 764-0500 FAX: (313) 647-2106 several, refer to individual reports 1

FEMATS Background and Overview FEMATS is a general-purpose finite element-ABC code for scattering analysis and was developed under the sponsorship of the EMCC through the NASA NRA grant NAG 2 -866. It is the first successful finite element analysis code relying on fully sparse matrices. To our knowledge, it is one of the most validated codes and rivals the sophistication and robustness of industry codes developed under the EMCC PRDA. The main advantage of FEMATS over other frequency domain codes is its inherent ease in modeling material structures, its geometrical adaptability and its low O(N) memory requirements. It can serve as the platform for developing new classes of hybrid finite element codes(as proposed here) which employ a variety of mesh termination schemes and other improvements relating to adaptive elements and error control. It can be further interfaced with other finite element codes such as FEMA-PRISM for antenna analysis and RADOME for including effects due to the nose radome antenna. Under the above NRA, U-M successfully * Developed, implemented and demonstrated the accuracy of a new class of ABCs for terminating finite element meshes. In contrast to previous ABCs, this new class of ABCs can be enforced on a surface conformal t the scatterer and yield remarkably accurate results even when the mesh is truncated 0.3k to 0.4k away from the scatterer. This implies a substantially reduced storage as well as execution time, and is consequently the main reason that our implementation has led to a very successful code. * Calculated scattering patterns for nearly all EMCC benchmark geometries and compared them with reference data (measured or calculated). For all geometries, good agreement was observed even when the ABC/PML was placed at 0.3k away from the scatterer's surface. All of our calculations and data files have been archived for future use by the EMCC and the file names and location are listed in the submitted test-case manual. FEMATS was used to successfully run some of the largest FEM problems to date with systems involving more than 500,000 degrees of freedom. * ported FEMATS on four different parallel platforms. Namely, it was parallelized and run on the KSR1, CM5, Intel Paragon and Comvex Exemplar. In fact, most of the benchmark geometries were executed on the parallel platforms. (see Table 1 for FEMATS performance on the various parallel platforms). * Implemented, incorporated and parallelized the latest sparse solvers in FEMATS. We concluded that the symmetric and unsymmetric biconjugate gradient (BCG) and Quasiminimal residual(QMR) solvers were very effective for FEMATS. Typically N/100 to N/80 iterations are needed for convergence (N=number of degrees of freedom) which is a remarkably fast convergence rate. However, the convergence rate of BCG is not monotonic and for monotonic (but equally robust) implementation, the QMR solver is also available with FEMATS. Recently, a new sparse matrix solver (referred to as 2

CVSS) was introduced and is purported to have much faster convergence rates. We will therefore examine the performance of this new solver for EM applications * interfaced FEMATS with sophisticated commercial mesh generation packages and this is quite important for the code's portability. Indeed, using such packages, FEMATS is capable of modeling geometries of any degree of complexity and material composition. This, however, is as much of an advantage as a disadvantage. Because the code has no hardware restrictions on the input geometry definition, the user is typically required to enter every aspect of the geometry definition in a step-by-step process which is often a time consuming task. To overcome this issue, FEMATS was upgraded (as part of this project) to incorporate a public-domain prismatic mesher to automatically grow the volume mesh and to avoid use of commercial meshing packages. This mesher has already been tested in connection with some of the EMCC benchmark geometries and in conjunction with the PML for mesh truncation as described in the attached documents. * As part of this project we completed a graphical users interface (GUI) for FEMATS to facilitate the code's utility by the electromagnetics community. The FEMATS GUI allows for mesh generation, EM parameter specification, code submission and monitoring while running on remote parallel platforms, graphical display of surface fields, plotting of RCS data and so on, all in a single integrated environment. Project Goals and Progress Summary The goal of the project is 1. to enhance the University of Michigan FEMATS code with a graphical interface(GUI), 2. make FEMATS compatible with a non-commercial structured mesher 3. upgrade FEMATS with improved mesh truncation schemes 4. test and validate FEMATS using EMCC benchmarks such as the almond and VFY218 Graphical Users Interface The FEMATS GUI is now operational and we are continuing to add more features to it. The items included in the GUI are listed in the brief document appended to this report along with a display indicating a typical pull down menu. At this moment, all functions necessary to run FEMATS (geometry input, meshing, compilation, platform submission, plotting, near zone field images) are included in the menu. We are presently still working on finalizing the capabilities of the GUI with respect to data post-processing. A more detailed report on the FEMATS GUI will be provided early September 1996. It is important to note that the FEMATS GUI is completely based on non-commercial packages can therefore be ported on any desired platform supporting x windows. 3

Prismatic Mesher The prismatic mesher has been completed and fully tested as reported in the Radiation Lab Techn. Report 033288-1-T. The FEMATS GUI permits the use of ACAD for generating the surface mesh. In turn the surface mesh can be imported into the NASA Ames code PRISM for automatic generation of the volume mesh. The user can define material regions, metallic or absorber sections, as necessary, on the basis of the original geometry specifications. At this point, the prismatic elements are subdivided into tetrahedrons before code execution. The user is allowed some control on the mesh structure by modifying a number of scaling parameters as described in the aforementioned report. Mesh Truncation Upgrades We fully upgraded FEMATS mesh termination schemes to include the new perfectly matched absorber and we are now working on importing the latest fast algorithms for mesh truncations. Currently, FEMATS employs conformal (doubly curved) absorbing boundary conditions and artificial absorbers for mesh truncation. This conformal applications for mesh truncation was introduced by U-M 2-3 year ago and resulted in substantial CPU and memory reductions. Most of our effort under this project was devoted to investigations for assessing their accuracy of the perfectly matched absorber (referred to as PML) and its application to doubly curved surfaces. We first implemented the PML for mesh truncation of microwave structures and have developed guidelines and tables for choosing its parameter to achieve broadband absorptions below -50dB without increasing the unknown count. The incorporation of the anisotropic PML into FEMATS was rather straightforward and did not require much effort. However, the evaluation of the absorber for non-planar implementations has yet to be done and therefore of interest has been the assessment of its accuracy and appropriate choice of parameters such as thickness and absorber-to-target distance as a function of curvature. This type of assessment has been carried out (numerically and analytically) and a preliminary report on the performance of the PML for cylindrical implementations is included. Basically, we found that the performance of the PML is compromised for non-planar implementations. Nevertheless, the PML still remains more attractive than absorbing boundary conditions in terms of accuracy, ease of implementation and parallelization. One of the issues in connection with PML mesh truncation has been their slower convergence rate. To address this concern, we are currently investigating convergence acceleration schemes and special preconditioning algorithms. The recent introduction of fast integral algorithms provided for a renewed interest on boundary integral mesh truncation schemes. As is well known, boundary integrals lead to fully populated submatrices which are CPU intensive. However, implementations based on the fast multipole method (FMM) and adaptive integral methods (AIM) reduce the CPU time from O(N2) down to O(N1'5) or even lower by resorting to multilevel subgroupings. We therefore began a new investigation aimed at the hybridization of the finite element and fast integrals methods. This type of implementation combines the geometrical adaptability and material generality of the finite element method with the robustness of the integral equation for mesh truncation. Thus, provided the CPU time of the boundary integral is 4

affordable, the proposed hybridization is the most attractive among the available competing hybridizations. At this point, we have completed a study of the finite element-boundary integral formulation using three different types of fast multipole methods. This was carried out for two-dimensional scattering and it was determined that a new version of the fast multipole method based on an asymptotic evaluation of the Green's function provides for the best compromise in terms of speed and accuracy. We have now began with the implementation and examination of these schemes for three dimensional structures. We have already completed a radome scattering code using a three dimensional version of the FMM and we are now implementing hybridizations of AIM and the finite element method. FEMATS Benchmarking FEMATS has been benchmarked against all EMCC geometries except for the Almond and the VFY218. These structures require conformal mesh terminations for efficient implementation. Therefore, use of conformal PML and boundary integral truncations are most appropriate for these targets. As stated later in this report, we have already validated FEMATS for the NASA almond and intend to examine the performance of FEMATS in evaluating the performance of antennas mounted on the VF218. The latter is a new application and should be of more interest to the EMCC. 5

LIST OF DOCUMENTS IN THIS REPORT Graphical Interface to FEMATS (brief report) M. Nurnberger and J. Volakis This document provides a brief description of the current FEMATS GUI along with flow charts and example GUI displays Artificial Absorbers for Truncating Finite Element Meshes J Volakis, T. Senior, S. Legault, T Ozdemir and M. Casciato This document gives an overview of the PML absorber, its capabilities and application to various antenna and scattering problems, including the NASA almond and antennas on conical platforms Design of Planar Absorbing Layers for Domain Truncation in FEM Applications S. Legault, T. Senior and J. Volakis Design guidelines are given for the PML absorber. Curves and formula are derived which can be used to select the PML parameters to yield a desired performance. Application and Design Guidelines of the PML Absorber for Finite Element Simulations of Microwave Packages J. Gong, S. Legault, Y. Botros and J. Volakis Demonstrates the validity of the PML design curves presented in the previous document for applications to microwave circuits. The latter application can be carried under a controlled environment and is therefore better suited for validating the PML design curves Performance of an Artificial Absorber for Scattering by a Metallic Cylinder A. Brown and J. Volakis A first application of the PML for mesh truncation on a non-planar boundary. The exact solution of the circular metallic cylinder is used to gauge the accuracy of the PML. Guidelines are given for the thickness and 6

appropriate distance of the PML from the cylinder surface to achieve a desired level of accuracy. Comparison of Three Fast Multipole Techniques for Solving Hybrid FE-BI Systems S. Bindiganavale and J. Volakis Various Fast Multipole Techniques(FMM) have been introduced over the past 3 years. These techniques often compromise accuracy for speed and in this document we provide a critical comparison of the fast multipole techniques in terms of accuracy and speed-up. It is found that the new asymptotic or windowed FMM provides the best compromise in terms of accuracy and speed. Hybrid Finite Element Methodologies for Antennas and Scattering J. Volakis, T. Ozdemir and J. Gong A lengthy up-to-date review of hybrid finite element methods for scattering and radiation. This document is a concise look at Michigan's contributions, including mesh truncation schemes, feed modeling and parallelization. Many applications are included ranging from antenna radiation to radome performance evaluations and engine scattering. 7

FEMATS: Finite Element MethodAbsorbing Termination Surface FEMATS ATTRIBUTES * General purpose scattering code * Employs first and second order tetrahedrals as building blocks (surfaces are faceted) * New ABCs are used to determine the mesh conformally to target only 0.3k away from target surface * Entire matrix is very sparse (storage: 8.5N complex #) * Models resistive sheets, impedance sheets, inhomogeneous and anisotropic materials, structural details * Solver has been parallelized for KSR, Intel iPSC/860 CM-5, Intel Paragon, IBM SP-2, Convex Exemplar * Convergence achieved in -N/100 (N = # of unknowns) * Solver is based on the BiCG or QMR method with ILU preconditioning * Validated for composite structure and inlets with up to 600,000 unknowns 8

ftf%. y TABLE 1 Time taken per iteration per unknown for the solver \ on FEMATS Machine name No. of Processors Time (isecs/iteration/unknown) Cray C-90 1 0.55 KSR1 28 1.28 58 0.57 Convex SPP-1000 8 3.5 Intel Paragon 8 3.42 16 1.99 32 1.38 IBM SP-1 4 1.47 II,\ // Radiation Laboratory 14 The University of Michigan

Brief Report on the Graphical Interface to FEMATS M. W. Nurnberger and J.L. Volakis 1 Introduction Graphical User Interfaces(GUIs) are have become important utilities in all user-oriented codes which provide multifunction capabilities. FEMATS is such a code, bringing together a number of different technologies and capabilities. Moreover, because the code will be able to run on serial, vector and parallel platforms a user interface to control and monitor the submission and execution on these platforms is essential. Briefly, a user will typically go through the following steps in an RCS/Radiation simulation: * defining the geometry using a solid modeling or an IGES/facet data file * specifying the material in selected regions * generating a surface or volume mesh and viewing it to examine fidelity * selecting a mesh truncation scheme and solver (Pro-Solver, BCG, FMM, AIM etc.) * specifying data runs (over frequency ranges, angles, output selection, etc.) and verifying final geometry data sets, including tessellation * submitting the problem to one of the available computating platforms * monitoring the execution (important for parallel platforms) * data extraction and display of RCS and surface fields using line plots and 3D dataset visualization software With the intent of delivering a code which can be used by entry-level engineers the use of a GUI is important as a tool for managing the various tasks related with EM analysis. Also, future code integration into a simulation environment which combines both electromagnetic, aerospace and other disciplines necessitates such a graphical interface. 10

2 GUI functionality Here is a summary of the capabilities which best characterize the FEMATS GUI: * All operations are point and click using pull down menus and with all the necessary functionalities accessible from a single windowing environment * User will have control over meshing, and ability to verify/validate mesh at any time. * Ability to set and modify material parameters easily * Keeps track of all data management and code interface with the various computing platforms allowing the user to concentrate on the physics of the problem and engineering issues. * Easy to learn and use using the visual displays but having the availability of keyboard short cuts for the experienced user. * Program files, for automated processing on most popular computing platfroms. * Supports multiple standard geometry input formats and volume/surface mesh generators (ACAI, IDEAS, PATRAN, IGES, DXF, etc...) * Keep track of various Solver packages supported by FEMATS * Supports all custom geometry generators through a standardized file form.at, especially developed for this purpose (RLGDFF - RadLab Global Data File Format) * Its modular design allows for additions and upgrades, including new solvers, mesh generators, different elements, visualization packages and so on. * Makes maximum use of public domain and fully documented file formats, for ease in developing converter scripts, integration of extra modules, etc. * All portions of code (GUI, converter scripts, visualization) available for nearly all machines, and easy to port if necessary. Maximum use of public supported packages for visualization, mesh generation pre/post-processing is important in this regard. Also, commercial meshing packages are generally available on all popular platforms. * Allows for pre/prost-processing on local machines and keep track of code submission and execution on remote parallel/vector/workstation(PVM) platforms. * Supports 2D plotting package for radiation/scattering patterns, etc. Supports integrated 3D plotting package for mesh validation using, for example, GeomView, and will have the option to easily integrate others such as SGI OpenGL-based Open Inventor, Performer, etc. 11

* Will support a fully integrated 3D surface viewer(including powefull commercial packages such as MATLAB) of surface/volume fields and will have capability to easily integrate with other packages 3 GUI format An overview of UM's graphical display depicting some of its pull-down menu capabilities is shown in Figure 1. As seen, it integrates post-processing, pre-processing and specification of the input files and formats using pull-down menus similar to the popular word-processing packages. The UM GUI is written in the Tcl/Tk (pronounced tickle-tee-kay) language, and uses the Tix widget collection. Tcl is a interpreted script language, designed particularly for the integration and simultaneous control of several application codes, each possibly running on a different platform, making it very well suited for design and analysis tools such as FEMATS. Tk is an X11 toolkit based on Tcl, designed to work alongside Tcl and provide a graphical interface to the software controlled by the Tcl script. Tk also provides the hooks to allow program intercommunication for greater amounts of interaction between code modules. The Tix widget collection is a series of big, complex widgets, built up from the smaller widgets supplied by Tk. Thus rather complex tasks can be specified using simple commands, greatly simplifying the mundane tasks. The GUI is designed in a way that allows it to be extended through the addition of new menu items for other function controls or the addition of new modules as the application program is enhanced. Given that Tcl/Tk is a well described language, a knowlegeable user can make modification to the UM GUI at any stage of the development or in the future. This is a unique feature of the proposed GUI and is crucial to its porting on any desireable platform. Moreover, it can be easily adapted for other electromagnetic analysis packages presently used by the EMCC or to be adapted in the future. Ostensibly, the whole GUI is a very large heirarchy of short (or long, depending on what they do) scripts, each one implementing a certain function of the overall program which can be activated by clicking on the appropriate button. Indeed, these scripts can be written not only in Tcl but also using languages such as Perl, C, FORTRAN and so on. The table below provides a list of the presently available and near-term capabilities of the FEMATS GUI. The items in the list are self explanatory. Some functions not yet operational are listed to indicate that they are are important and should be considered in the near future. Ability Now Near Term General Defaults File X X Status/Project File X X Working Project Management X X Intelligent Interface - X Keyboard Shortcuts 12

Program Files Graphical Window Dumps? Mesh Post-Processing X X Material ID (setting) X X Material ID (viewing) - - Runtime Parameter Input X X Solver Parameter Input X X Local Solver Initiation X X Remote Solver Initiation? 2D Visualization (GnuPlot) X 3D/Surface Visualization (xvgr, GeomView) -some Online Help Interfacing to External Codes Geometry Generators ACAD? IDEAS - X DXF IGES others... Meshing packages internal meshing - X Prism - X others... Converters RadLab Global Data File Format (RLGDFF) Graphical Output 2D Viewing, Line Plots (GnuPlot) - X 3D Viewing (xvgr, GeomView) - X Surface Field Plots (Geomview) - X Solver Control Local Solver, Single CPU X X Local Solver, Distributed Remote Solver Solver Availability Scalar Solver X X Vector Solver (Cray) Parallel Solver(s): Intel Paragon 13

Convex SPP X X IBM SP/2 X X SGI Power Challenge - - PVM 14

15

F Geomet ry * _ Generator N Convert To/From R.L. Global Data File Format (RLGDFF) FEMATS Graphical Interface Control Diagram 16

ARTIFICIAL ABSORBERS FOR TRUNCATING FINITE ELEMENT MESHES J.L. Volakis, T.B.A. Senior, S.R. Legault, T. Ozdemir and M. Casciato Radiation Laboratory Department of Electrical Engineering and Computer Science 1301 Beal Ave. University of Michigan Ann Arbor, MI 48109-2122 Abstract A metal-backed layer of absorbing material offers a number of advantages for truncating the computational domain in a finite element simulation. In this paper we examine isotropic and anisotropic absorber layers for the purpose of truncating finite element meshes. Optimal design curves are presented for these absorbers that can be used to select the various parameters (thickness, propagation and sampling rate) so the reflectivity is minimized. Applications to radiation and scattering problems are also considered, and these illustrate the accuracy and versatility of the absorber layers for large scale electromagnetic simulations. 1 Introduction In the numerical solution of electromagnetic scattering and radiation problems it is necessary to truncate the computational domain in a manner which ensures that the waves are outgoing. This is true also in the analysis for many microwave circuits, and the need to terminate the mesh is common to finite element (FEM) and finite difference-time domain (FDTD) methods. One way to do so is to enforce an absorbing boundary condition (ABC) at a surface as close as possible to the scatterer or radiator, and a review of available ABCs has been given by Senior and( Volakis [1]. Unfortunately, ABCs are limited in their ability to conform to the surface of the scatterer. They may also require an a priori knowledge of the wave's properties and, in an FEM solution, they generally reduce the rate of convergence. Another way of terminating the mesh is to use a metal-backed layer of isotropic or anisotropic absorbing material [2,3,4,5], and such layers are often referred to as artificial absorbers since their material parameters may be physically unrealizable. The implementation of artificial absorbers for finite element mesh truncation is illustrated in Fig. 2 and as expected the layer's material composition plays a major role in the performance of the artificial absorber. However, the chosen numerical discretization of the absorber has an equally important role and cannot therefore be ignored in the design of the absorber for numerical simulations. In this paper, we examine the performance and design of both isotropic and anisotropic homogeneous absorbers for truncating finite element meshes. Their application to the finite element solution of three dimensional radiation and scattering problems is also considered and results are shown which demonstrate the utility of these mesh truncation schemes. 17

2 Analytical Study Consider the metal-backed planar layer shown in Fig. 2(a). The surface x = 0 is the interface between free space (in x < 0) and a lossy material (in x > 0) backed by a perfect electric conductor at x = t. For an incident plane wave E, or H, = e-jko(xcos +Y sin ) (1) the reflection coefficient is RE(O) or RH((), and the objective is to minimize these. If the layer is composed of a homogeneous isotropic material whose relative permittivity Cr and relative permeability /r are such that Or = Fr = b = a - j3 (say), then \/1 - b-s2 sin2 - j cos & tan(kobt/l - b-2 sin2 ) ( RE,= (2) j1- b-2 sin2 + + j cos 4 tan(kobt l - b-2 sin2 ) R() /1b-2 sin2 - j cos ) cot(kobt 1 - b-2 sin2 1) RH = - (3) /1 - b-2sin2 + jcos )cot(kobt/ - b-2sin2) These differ because the presence of the metal backing has destroyed duality, and at grazing incidence (: = 7r/2), RE = RH = -1. If sin & > 1, i.e. = 7/2 + jS with S > 0 so that sin' = cosh S and cos = -j sin S, RI differs from unity by only a small amount for both polarizations. The behavior of IRE,H(O)I as a function of sin& is illustrated in Fig. 3. At normal incidence (q = 0), (2) and (3) give RE,H(O) = Te-2jkot(-j) (4) whose magnitudes are independent of a and can be made as small as desired by choosing kot/3 sufficiently large. Since large values of 3 produce rapid field changes in the dielectric, a disadvantage is that high (and often very high) sampling rates are necessary to simulate them numerically. As kot -4 oo, RE,H(b) - O0 only for normal incidence, but Sacks et al [4] have shown that a particular uniaxial anisotropic material has this property for all 4 < 7/2. The result is an example of a perfectly matched layer (PML), and if Cr, = u b - (b - -)Xx (5) where I is the identity tensor and b = a - j, then RE,H() = Te-2jkot(-j3)cos (6) which reduces to (4) in the particular case of normal incidence. If ko0 >> 1 the reflection coefficients decay exponentially for all q < 7r/2, and since (6) is also valid for sin / > 1, the choice a > 0 ensures an exponential decay for these angles as well. The behavior of IRE,H(G)I is illustrated in Fig. 3 for the same values of kot, a and 13 used for the isotropic layer. Clearly, a major advantage of the PML is that its reflection coefficient remains low for a wide range of angles of incidence. 1R

3 Absorber Design The theoretical behavior of the reflection coefficient \RI for the metal-backed layers is relatively simple. In the case of the isotropic material, an increase in /3 and/or t decreases [R(0)1. The uniaxial material has this behavior for all real angles of incidence, and while a plays little or no role, large values of a do produce higher absorption for complex angles. It follows that for a uniaxial layer of given thickness t, a and /3 can be chosen sufficiently large to produce high absorption over a broad angular spectrum, with angles near b = r/2 providing the only exception. Unfortunately, the analytical results do not immediately translate into numerical performance. Because of the discretization inherent in an FEM implementation, the fields inside the layer are reproduced only approximately, and this is particularly true for a rapidly decaying field. To design a good absorber it is necessary to understand the impact of the sampling rate on the choice of ac, 3 and t, and our objective is to find the minimum number of sampling elements ( or discrete layers) to achieve a specified IRI. It is anticipated that the errors introduced by the discretization will have a number of consequences. In particular, for a given number N of discrete layers and given t, increasing /3 will ultimately lead to an increase in \RI because of the inability to model the increasing attenuation, and an increase in a will likely produce a similar effect. To obtain some insight into the roles played by N, a, 3 and t, a simple FEM model for computing the reflection coefficient of a metal-backed absorber layer (see Fig. 1) was considered. We examine first a homogeneous isotropic layer at normal incidence for which the theoretical reflection coefficients are given in (4). In spite of the fact that the magnitudes are the same for both polarizations, a polarization dependence shows up in the FEM implementation. This is illustrated in Fig. 4, and we note that as N increases, the FEM values of IR(0)| converge to the common theoretical value for both polarizations. The nulls associated with the Hpolarization curves are characteristic of the behavior of the modeled absorber layer, and are due to the interference of the reflected fields from the dielectric interface, the metal-backing and the individual layers used to model the absorber numerically. Given the many parameters (/3, a, t, N), it is essential to consider some optimization of the proposed metal-backed absorber as a function of these parameters. For this purpose, a study was carried out using the FEM code mentioned above. Initially, an investigation was performed to examine the effect of each parameter on the absorber's performance. As can be expected, a plays little role in the performance of the absorber for relatively small values of /3t/Ao. Also, larger values of N for the same thickness t lead to lower reflection coefficients. An optimum value exists but this depends on /3, the absorber's decay parameter. However, the most important observation from this preliminary study concerns the scalability of the product /3t. It was found that for small values of a and for thicknesses up to at least 0.5 wavelenths, the reflection coefficient is a function of the product 3t/Ao alone. As a result, one can construct reflection coefficient curves as a function of 3t which are optimum for a given N. Fig. 5 gives such design curves by plotting IR(0)I and N versus 3t/Ao on the same figure, and these can be used to determine the optimum N for a given 3t/A\. To see how to use the figure, suppose that the desired reflection coefficient at normal incidence is -50 dB. In Fig. 5 we observe that the IR(0)I curve intersects the -50 dB line at 3t/Ao v 0.58, and referring now to the N curve, the number of elements required is N = 10. The value of 3 can then be found by specifying either 10

the element size or the layer thickness. Thus, for elements 0.025Ao thick, we have t = 0.25Ao and 3 = 2.32. By increasing N we can improve the performance up to the limit provided by the theoretical value of |R(0)| which has been included in Fig. 5. A good approximation to the long-dashed curves in Fig. 5 obtained by linear regression is t -0.01061RI + 0.0433 (7) Ao N 0.147 exp [7.3533t/Ao] (8) where JR\ is measured in dB and N is the smallest integer equal to or exceeding the right hand side of (8). 4 Application to Antennas and Scattering The above study dealt with planar absorber layers but the layers can also be curved for terminating finite element meshes in a conformal manner. This is illustrated in Fig. 1 where the artificial absorber layer is used to terminate the computational domain surrounding a patch antenna. Conformal mesh terminations are quite attractive in FEM modeling because they lead to a substantial reduction of the computational volume and are also compatible with structured as well as unstructured meshes. In contrast to absorbing boundary conditions, they do not require a knowledge of the closure's principal radii of curvature and uniaxial artificial absorbers offer the possibility of reflectivity control as low as -60dB. To illustrate the applicability of the artificial absorbers in FEM modeling, two examples are considered, one dealing with antenna analysis and the other with scattering by a non-canonical structure. For simplicity both cases employed a simple version of a homogeneous artificial absorber consisting of three layers 0.05Ao thick as shown in Fig. 1. The attenuation constant for each layer was f3 = 2.6, and thus ~ = 0.405. From Fig. 5 this absorber provides a normal incidence reflection coefficient of -30dB, and from the same figure we also read off that this reduction can be realized with 3 layers or more. Thus, the design of the proposed homogeneous absorber is consistent with the curves given in Fig. 5. The absorber termination shown in Fig. 1 has been used to model a variety of patch (circular and rectangular) antennas on doubly curved platforms, including circular, spherical, conical and ogival surfaces. Of particular interest is the computation of the resonant frequency which is a rather sensitive quantity and its accurate computation via the proposed FEM model provides a good test of the absorber's performance. Fig. 6 shows the results for a rectangular patch antenna mounted on the conical surface illustrated in the figure. The patch resides on a substrate having cr = 2.32 and a thickness h = 0.114cm. Its dimensions are given in Fig. 6 and on the basis of the approximate cavity model it resonates at 3.2GHz. From the computed input impedance plot, it is seen that the resonance frequency predicted by the FEM code is 3.115GHz, which is within 3 percent of the cavity model. The FEM computations were carried out using prismatic edge elements and the surface grid is also shown in Fig. 6. A total of 2358 prisms were used for this analysis resulting in 3790 degrees of freedom. As a scattering example we consider the radar cross section of the NASA metallic almond [6] shown in Fig. 7. This body is 9.936 inches long, and precise formulae for describing its surface are given by Woo et. al.[6]. Measured data are also available at several frequencies and 20

can be used to benchmark the accuracy of the simulation. Our finite element code FEMATS [7] was used to model the almond illuminated with a plane wave at a frequency 1.19GHz. This code employs edge-based tetrahedral elements and the mesh was again terminated using the aforementioned 3-layer homogeneous artificial absorber. For ease of mesh generation, a structured prismatic mesh was first generated conformal to the almond's surface and consisted of nine 0.05A layers with the outer layers occupied by the artificial absorber. The structured prismatic mesh was then turned into a tertrahedral mesh resulting in a total of 46,878 edges. Fig. 7 displays the computed radar cross section(RCS) computations with the measured ones for both polarizations of incidence. The patterns are taken in the plane most parallel to the flat side of the almond(non-symmetric plane), with zero degrees corresponding to incidence tip-on. As seen, the calculations are in good agreement with the measured data except at near 90 degrees for HH polarization. However, other reference calculations based on a moment method code are in agreement with the FEMATS data, suggesting that the discrepancy may be due to minor alignment errors in the measurement. 5 Conclusion In this paper we investigated homogeneous isotropic and uniaxial artificial absorbers for finite element mesh truncation. By properly selecting the material properties and sampling rate, it was demonstrated that almost any desired level of absorption can be attained, and typically very few samples (less than five) are needed to achieve a reflection coefficient of -40 dB over a wide range of incidence angles. Design curves were presented which can be used to select the various parameters (loss, thickness and sampling rate) on the basis of a desired reflection coefficient. As expected, a lower material loss requires a thicker absorber to produce the same reflection coefficient but, on the other hand, a higher attenuation rate requires more samples to attain a lower reflection coefficient in a numerical implementation. Most likely, inhomogeneous (tapered) artificial absorbers can lead to lower reflection coefficients, but these have not yet been investigated to any great extent. In contrast to absorbing boundary conditions, a particular advantage of the proposed absorbers is that they can be used to terminate a finite element mesh conformal to the target or radiator surface without needing a priori information about the wave's propagation characteristics. To test the performance and applicability of the proposed absorber for truncating finite element meshes in a three dimensional setting, two examples were considered.Namnely, computation were carried for the input impedance of a patch antenna on a conical surface and the radar cross section of a non-canonical slender body. In both cases the computed values were in good agreement with reference data by using the proposed artificial absorber for mesh truncation. References [1] T.B.A. Senior and J.L. Volakis, Approximate Boundary Conditions in Electromagnetics. Stevenage, UK, IEE Press, 1995. 21

Ec / - a = 1-j2.7 0.15X / s Eb= c(l-j2.7) b,= 1-j2.7 ~b= E (1-j2.7) 4b= -j2.7 4 I PEC (Actual geometry) (FEM Simulation) Figure 1: Illustration of the finite element mesh truncation using an artificial absorber. [2] C. Rappaport and L. Bahrmasel,"An absorbing boundary condition based on anechoic absorbers for EM scattering computation," J. Electromagn. Waves Appl. 6, pp. 1621-1634, 1992. [3] T. Ozdemir and J.L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Sci. 29, pp. 1255 -1263, 1994. [4] Z.S. Sacks, I.M. Kingsland, R. Lee and J.F. Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition," to appear in IEEE Trans. Antennas Propagat. [5] J.P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comp. Phys. 114, pp. 185-200, 1994. [6] A. Woo, H.T.G. Wang, M. J. Schuh and M. J. Saunders, "Benchmark radar targets for the validation of computational electromagnetics programs," IEEE Antennas and Propagat. Magazine, 35, pp. 84-89, Feb. 1993. [7] A. Chatterjee and J.L. Volakis, "Conformal absorbing boundary conditions for 3D problems: Derivation and applications," IEEE Trans. Antennas and Propagat., 43, pp. 860-866, Sept. 1995. 22

ty y. N 1... N x x Figure 2: Geometry of the metal backed absorber layer and its finite element discretization. 0 -10 -20 -30 -40 -50 -, cr -60 - -70 -80 0.0 0.2 0.4 0.6 sin o 0.8 1.0 1.2 Figure 3: Analytical results for t 0.25Ao and b= 1- j2: (...) homogeneous isotropic and anisotropic absorbing layers with isotropic E pol, (- * -) isotropic H pol. and ( —) anisotropic. 23

0 -10 -20 -30 -40 -50 -60 -70 ( 0.0 0.2 0.4 0.6 0.8 1.0 1.2 sin o Figure 4: Numerical results for a homogeneous isotropic and anisotropic absorbing layer with t = 0.15Ao and b = -j2.5. The top four curves are for Epol., the bottom four for H pol.: ( ---) exact, (- - -) N=3, ( — -) N=6, and (- — ) N=12. -20 -30 -40 -50 -60 -70 30 25 20 15 10 5 Co a: z -80 L 0.3 0.4 0.5 0.6 0.7 3t/ Xo Figure N: ( 5: Absorber design curves. The straight lines -) exact IRI; (- -) numerical \RI give IRI in dB, and the curved ones give 24

Metallic Ground Plane Patch 2 /h = 0.114cm V-,..... / 200 Ia o u = 100 -0 CA 0 0&-j 0 0 rb= 12.9cm r = 10cm 0,= 18.34~ 3.0 3.1 3.2 (GHz) 3.3.. ---- Y (a) (b) Vertically directed electric field magnitude in the substrate Metal Termination Boundary Absorber-air interface Microstrip patch Feed location (c) Figure 6: Microstrip sectoral patch on cone: (a) configuration, (b) patch input resistance as a function of frequency for different cone angles and (c)display of the mesh vertical fields below the patch and in the substrate

(a) RCS - NASA METALLIC ALMOND 80 100 120 140 160 180 Azimuth (Degrees) (b) Figure 7: Geometry and RCS results for the NASA almond[6]. (a) surface and volume prismatic conformal mesh of the 9.936 almond and (b) RCS pattern calculations and measurements at 1.19GHz. 26

Design of Planar Absorbing Layers for Domain Truncation in FEM Applications S.R. Legault, T.B.A. Senior and J.L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan, Ann Arbor, MI 48109-2122 ABSTRACT A metal-backed layer of absorbing material offers a number of advantages for truncating the computational domain in a finite element simulation. In this paper we present design curves for the optimal selection of the parameters of the layer to achieve a specified reflection coefficient. The curves are based on one-dimensional finite element simulations of the absorbers, and the optimization is therefore a function of the sampling rate. Three types of material are considered, including the recently introduced perfectly matched uniaxial material, either homogeneous or with a quadratic material profile. Two three-dimensional applications are also presented and used to examine the validity of the design curves. 1 INTRODUCTION In the numerical solution of electromagnetic scattering and radiation problems it is necessary to truncate the computational domain in a manner which ensures that the waves are outgoing. This is true also in the analysis for many microwave circuits, and the need to terminate the mesh is common to finite element (FEM) and finite differencetime domain (FDTD) methods. One way to do so is to enforce an absorbing boundary condition (ABC) at a surface as close as possible to the scatterer or radiator, and a review of available ABCs has been given by Senior and Volakis [1]. Another way is to use a metal-backed layer 27

of isotropic absorbing material [2,3], but both schemes have limitations. For example, an ABC requires a priori knowledge of the propagation constant which, in a microwave problem, may not be the same in all section of the computational domain. Also, when used to terminate an open domain, ABCs reduce the convergence rate and may be hard to implement on a surface conformal to the scatterer or radiator. An isotropic dielectric layer alleviates some of these difficulties, but its accuracy and aspect coverage are limited. Recently a new anisotropic absorber has been proposed for termninating the domain. By introducing an additional degree of freedom, Sacks et al. [4] have shown that a uniaxial material can be designed to have zero reflection coefficient at its interface for all angles of incidence. If the material is also lossy, a thin metal-backed layer can be used to terminate an FEM mesh, and though the material is no longer realizable physically, the associated fields are still Maxwellian. This is often referred to as a perfectly matched layer (PML), and its development was motivated by the non-Maxwellian layer introduced by Berenger [5] (see also [6]) for FDTD problems. By choosing the parameters appropriately, it is possible to achieve any desired level of absorption for almost all angles of incidence using only a thin layer, but its numerical simulation is a more challenging task. Because of the rapid exponential decay of the fields within the layer, there are large variations in a small distance, and it is difficult to reproduce these in a numerical simulation. Thus, for a discretized PML, the numerical sampling as well as the material properties affect the reflection coefficient that is achieved. In this paper we consider the design and performance of three types of metal-backed planar layers for terminating FEM meshes: homogeneous isotropic, homogeneous anisotropic (uniaxial), and inhomogeneous (tapered) uniaxial materials. Using one-dimensional finite element simulations, their numerical performance is examined and coinpared with their theoretical capability. Not surprisingly, the sampling rate has a major effect on the reflection coefficient. Based on a detailed numerical study, we identify scalable parameters in the numerical model and use these to generate design curves and formulas for choosing the sampling rate and material properties to achieve a specified reflection coefficient. As expected, a tapered uniaxial material proves superior to the homogeneous one. The applicability of these results to threedimensional problems is then illustrated for a simple microwave circuit and a rectangular waveguide. 28

Ay y x x (a) (b) Figure 1: Geometry of the metal-backed absorber layer (a) and its FEM implementation (b). 2 ANALYTICAL STUDY Consider the metal-backed planar layer shown in Fig. l(a). The surface x = 0 is the interface between free space (in x < 0) and a lossy material (in x > 0) backed by a PEC at x = t. For an incident plane wave Ez or Hz = e-ko(cos5+ysin7) (1) the reflection coefficient is RE(O) or RH(4), and the objective is to minimize these. If the layer is composed of a homogeneous isotropic material whose relative permittivity Or and relative permeability /7 are such that Or - ir = b = a - j3 (say) where a and 3 are real, then R(5) - / 1 - b-2 sin2 - j cos 4 tan(kobt /1 - b-2 sin2 4>) V/1 - b-2 sin2 b + j cos q tan(kobtl - b-2 sin2 (2) ) /1 - b-2 sin2 4> - j cos X cot(kobt/1 - b-2 sin2) (3) RH -(3) /1 - b-2 sin2 4/ + j cos (j cot(kobt/l - b2 sin2 4) These differ because the presence of the PEC backing has destroyed duality, and at grazing incidence (4 = 7r/2), RE = RH = -1. If sinq > 1, i.e. q = Tr/2 + j with S > 0 so that sin4 = coshS and cos = -jsin, IRI differs from unity by only a small amount for 29

0 -10 -20 -30 -40 -50 -60 \ -70 -80 I 0.0 0.2 0.4 0.6 0.8 1.0 1.2 sin ) Figure 2: Analytical results for homogeneous isotropic and anisotropic absorbing layers with t = 0.25Ao and b = 1 -j2: ( ) isotropic E pol., (- -) isotropic H pol. and ( —) anisotropic. both polarizations. The behavior of IRE,H(q)I as a function of sin k is illustrated in Fig. 2. At normal incidence ( = 0), (2) and (3) give RE,H(O) = Fe-26jkt(-j ) (4) whose magnitudes are independent of a and can be made as small as desired by choosing kot/3 sufficiently large. As kot -+ oo, RE,H() -+ 0 only for normal incidence, but Sacks et al. [4] have shown that a particular uniaxial anisotropic material has this property for all q < 7r/2. The result is an example of a perfectly matched layer (PML), and if 6r = Jr = bl - (b - b )xx (5) where I is the identity tensor, then RE,H(4) =:e-2jkot(-j3)cos (6) which reduces to (4) in the particular case of normal incidence. If kof > 1 the reflection coefficients decay exponentially for all q < 7r/2, and since (6) is also valid for sin 4> 1, the choice a > 0 ensures an exponential decay for these angles as well. The behavior of IRE,H(O)I is illustrated in Figure 2 for the same values of kot, a and 0 used for 30

the isotropic layer. Clearly, a major advantage of the PML is that its reflection coefficient remains low for a wide range of angles of incidence. Although the outer surface of the layer is reflectionless for all >, the abrupt change in the material properties at x = 0 may produce a contribution in an FEM solution. We can eliminate the discontinuity by tapering the properties as a function of x to produce an inhomogeneous anisotropic layer. As shown by Legault and Senior [7], if b = - {xT(x)} a wave propagating into the material has the form e-jko {xY(x) cos +y sin } (7) and when 3 (8) 7(x)=l+( _-j-_1) 18) which tends to unity as x -- 0+, the reflection coefficients of the layer are identical to those given in (6). With this expression for y(x), the attenuation is less where the field is larger, i.e. close to the interface, and increases as the field is absorbed. A simplified version of (8) is employed in Section 3.4. 3 NUMERICAL STUDY For all three types of layer the theoretical behavior of IRI is relatively simple. In the case of the isotropic material, an increase in 3 and/or t decreases JR(0)I. The uniaxial material has this behavior for all real angles of incidence, and while a plays little or no role, large values of a do produce higher absorption for complex angles. It follows that for a uniaxial layer of given thickness t, a and 3 can be chosen sufficiently large to produce high absorption over a broad angular spectrum, with angles near - = 7r/2 providing the only exception. Unfortunately, the analytical results do not immediately translate into numerical performance. Because of the discretization inherent in an FEM implementation, the fields inside the layer are reproduced only approximately, and this is particularly true for a rapidly decaying field. To design a good absorber it is necessary to understand the impact of the sampling rate on the choice of a, / and t, and our objective is to find the minimum number of sampling elements (or discrete layers) to achieve a specified JRI. It is anticipated that the errors introduced by the discretization will have a number of consequences. In particular, 31

0 -10 -20 c -30 -70 -60 -70,- lI —, --- ——, --- —0 - 0.0 0.2 0.4 0.6 0.8 1.0 1.2 sin ) Figure 3: Numerical results for a homogeneous isotropic layer with t = 0.15Ao and b = -j2.5. The top four curves are for E pol., the bottom four for H pol.: ( — ) exact, ( — -) N=3, ( — -) N=6 and (-) N=12. for a given number N of discrete layers and given t, increasing /3 will ultimately lead to an increase in JR\ because of the inability to model the increasing attenuation, and an increase in oa will likely produce a similar effect. To obtain some insight into the roles played by N, a, /3 and t, we now consider a simple FEM model of the layers. 3.1 Numerical Model A one-dimensional FEM code was used to examine the numerical performance of the absorbing layers. The computational domain was limited to the discretized layer structure shown in Fig. l(b), with the appropriate boundary conditions applied at the interface x = 0 and the PEC backing x = t. We consider first a homogeneous isotropic layer at normal incidence for which the the theoretical reflection coefficients are given in (4). In spite of the fact that the magnitudes are the same for both polarizations, a polarization dependence shows up in the FEM implementation. This is illustrated in Fig. 3, and we note that as N increases, the FEM values of IR(0)I converge to the common theoretical value for both polarizations. 32

3.2 Dependence on a and 3 For a layer of constant thickness the theoretical value of IR(0)l is independent of a and polarization, but in the numerical implementation the behavior is much more complicated. Figure 4 shows ]R(0)I plotted versus a and 3 for a layer of thickness t = 0.25Ao made up of 5 (=N) elements, where the darker tones indicate lower values. For small 3 the results are in close agreement with theory. As evident from the level lines, ]R(0)I is almost independent of a and decreases exponentially with /3, leading to a linear decrease on a dB scale. For large 3, however, the behavior is quite different, and the most striking feature is the series of deep minima whose spacing in a increases with increasing a and decreasing 3. These are numerical artifacts which are common to both polarizations and may depend on the particular numerical code employed. The minima for the two polarization are interlaced, and for H polarization the first minimum occurs at a = 0, 3 = 1.6. Their locations also depend on t and N. If N is fixed, the spatial sampling is inversely proportional to t. Decreasing t results in better sampling, pushing the minima to higher values of 3 and producing agreement with the theoretical values for larger 3 than before. Increasing t has the opposite effect. On the other hand, if t is fixed, increasing N improves the accuracy, and shifts the minima to higher /3. Apart from the minima, the reflection coefficients for fixed 3 increase slightly with increasing a, and it is therefore sufficient to confine attention to the lower values of a. In Figure 5 the reflection coefficients are plotted as functions of / for the same layer with a = 0 and a = 0.75. The curves correspond to vertical cuts through the patterns in figure 4, and we also show the theoretical value obtained from (4). We observe that as / increases the reflection coefficients decrease initially at almost the same rate implied by (4), but beyond a certain point they begin to increase. The deep minimum at a = 0 and / = 1.6 in Figure 4(a) is clearly seen, but for design purposes it is logical to focus on the worst case, i.e. the polarization for which the reflection coefficient is larger. The upper curves in Figure 5 are almost identical and constitute this case. Since they correspond to two different values of a, either of them would suffice, but for reasons that will become evident later, we choose a = 0. 33

.........I......................... 1 0 - 4 10 4 -20 3-30 IRI [dB] -40 -50 0 -60 0 1 2 3 4 5 a (a) 5 -0 -10 -20 3 -30 IRI [dB] -40 -50 0' -60 0 1 2 3 4 5 a (b) Figure 4: Plot of IR(O)l in dB for (a) an H polarized and (b) an E polarized wave incident on a homogeneous isotropic layer with t = 0.25Ao and N = 5. The solid curves are level lines. 34

0 -10 -20. 0.75 -30 -40 -50 -60 0 1 2 3 4 5 Figure 5p: R(0)i with t = 0.25Ao and N = 5: ( ) Exact, ( f -) ta 0 Efpol., ( — -)ca 0Hpol., (- - -) a = 0.75 Evpol. and (s. T.) a = 0.75 H pol. 3.3 Dependence on 3co N and t We now seek a connection between the values of /, N and t for which R(0)y is minimized. To this end, we first examine |R(0)[ as a function of 5 and N for constant t, and the resulting plot is shown in Figure 6 for E polarization with t = 0.25A0 and a = 0 as before. For fixed 3 the reflection coefficient tends to its exact values as N increases. This is evident from the level curves and, as expected, the convergence is better for the smaller f. Considrer now the behavior of IR(0) for fixed N. As increases from zero, the reflection coefficient decreases to a minimum and then increases. The location of the minima is indicated by the solid line. This is consistent with the behavior shown in Figure 5 and the upper curve is, in fact, just a vertical cut through Figure 6 at N =5. The solid line in Figure 6 therefore gives the value of ft at which t R(0)o is a minimum as a function of the number of elements. If the process is repeated for other layer thicknesses, it is found that for minimum |R(0)| the curve of f/t/Ao versus N is virtually the same for all thin layers. The observation that /3t/\o is a scalable parameter is an important conclusion of our study, and by choosing a constant layer thickness we can produce a universal curve for the optimal choice of N and /3 in FEM simulations. Such a curve is shown in Figure 7 and can be interpreted as giving the value of 3t/Ao for a specified N to minimize the reflection coefficient [R(0)|. For example, if t = 0.2Ao 35

.......................I... l. 0 -20 2-30 IRI [dB] -40 16 -60 2 4 6 8 10 12 14 16 18 20 N Figure 6: JR[ as a function of /3 and N for E pol. with t = 0.25Ao and a = 0. The dashed curves are level lines and the solid curve gives the location of the minimum \R\. and N -= 3, then f = 2.13. If a smaller value of /3 is chosen, IR(0)\ will be larger (see Figure 5), and this can be attributed to the fact that the field reflected from the metal backing has not been attenuated sufficiently. If f3 is set to a value larger than 2.13, |R(O)[ will still be larger because the chosen N is too small to reproduce the rapid field decay within the layer. So far we have considered only normal incidence, but for the anisotropic layer it is a simple matter to extend the results to all real angles of incidence q < 90~. As evident from the exponent in (6), the absorption at normal incidence is reproduced at any angle if the layer thickness is inversely proportional to cos q>. This can be achieved by specifying the layer thickness t as a fraction S of the wavelength A, along the normal (or x axis) to the material interface. Since t = 8Ao/ cos 4 = -A\, t/Ax is now independent of Xb and the scalable parameter /3t/Xo (at normal incidence) becomes /3t/AX. For the anisotropic layer, all the results obtained at normal incidence are made applicable for arbitrary q by substituting Ax for A0. For example, plotting the optimum /3t/Ax as a function of N duplicates the curves shown in Figure 7. This notion can also be used to account for problems where the outer medium is not free space. In such cases, we have Ax = Ao/0/67 (at normal incidence) 36

0.7 0.6 0.5 0.4 0.3 0.2 I,,, I I 2 4 6 8 10 12 14 16 18 20 N Figure 7: 3t/Ao computed from the E pol. case with a = 0: ( —) t = 0.15Ao, (- -) t= 0.25Ao and ( --- -) t = 0.5Ao. where er is the relative permittivity of the outer medium. Although the scaling property of 3t/AX has been established only for a = 0, it holds to a reasonable degree for small a $ 0, but as a increases, the Pt/Ax versus N curves become increasingly dependent on a. The scalability also extends to the associated values of IR(q)|, and this enables us to provide a simple design prescription for an absorbing layer. 3.4 Design Curves Since the quantities /3t/A and IR(0)I are the same for layer thicknesses up to about 0.5AX at least, design curves can be obtained by plotting IR(q)I and N versus it/Ax on the same figure as shown in Figure 8. To see how to use the figure, suppose that the desired reflection coefficient at normal incidence is -50 dB. At ) = 0 we now have A, = Ao. In Figure 8 we observe that the IR(0)[ curve intersects the -50 dB line at ft/Ao _- 0.58, and referring now to the N curve, the number of elements required is N = 10. The value of 3 can then be found by specifying either the element size or the layer thickness. Thus, for elements 0.025Ao thick, we have t = 0.25Ao and 3 = 2.32. By increasing N we can improve the performance up to the limit provided by the theoretical value of IR(0)I which has been included in Figure 8. A good approximation to the short-dashed curves in Figure 8 obtained by linear 37

30 -30- 25 JRI -40 -50 m., 20 -70 - ------ 1. - / 5 -80......... 0.3 0.4 0.5 0.6 0.7 0.8 3t/x Figure 8: Absorber design curves. The straight lines give IRI in dB, and the curved ones give N: ( — ) exact |RI, (- - -) homogeneous case and (- -) inhomogeneous case. regression is t = -0.01061RI + 0.0433 (9) Ax N = 0.147exp [7.353/t/Ax] (10) where Ax = Ao/ cos b, IRI is measured in dB and N is the smallest integer equal to or exceeding the right hand side of (10). These equations hold for q, 0~ < (q < 90~, for the anisotropic layer, and 4 = 0~ for the isotropic one. The design criterion provided above applies to specific angles of incidence. In the case where a specific absorption level is required over a range of angles of incidence the layer must be designed for the largest angle occurring. Doing so ensures that the absorption is superior at all smaller angles. The performance can be improved by making the anisotropic material inhomogeneous, and to illustrate this we consider the case 7y(x) = -ji (x/t)2 for which the theoretical reflection coefficient is the same as before. The scalability is still preserved and the resulting curve is shown in Figure 8. The fact that the curve for the quadratically tapered layer lies below that of the homogeneous material confirms the improvement in performance, and we can now achieve a reflection coefficient of -50 dB by choosing ft/lAx _ 0.64 corresponding to N = 9. Approximations 38

to the long-dashed curves in Figure 8 are -0.0119l1R + 0.0451 (]1) Ax N =0.298exp [5.263/t/Ax] (12) where IR| and N are as before. Compared with the homogeneous material the decrease in the number of elements required becomes more pronounced as IR[ is reduced. As previously mentioned, these equations hold for all real angles of incidence in the case of the anisotropic layer and for normal incidence if the layer is isotropic. 4 THREE-DIMENSIONAL VERIFICATION As noted earlier, a PML is particularly attractive for terminating a finite element mesh in the simulation of microwave circuits. For these applications a PML has an advantage over a traditional ABC because it does not require an a priori knowledge of the guided wave propagation constant. To demonstrate the applicability of the design equations in three dimensions, we consider a shielded microstrip line and a rectangular waveguide. The microstrip line has width w = 0.71428 cm, substrate thickness 0.12 cm and relative permittivity e = 3.2, and is enclosed in a metallic cavity whose dimensions are shown in Figure 9. It should be noted that the height of the cavity from the microstrip line is sufficiently large to suppress any reflections from the cavity walls. As a result, the characteristic impedance of the line should be that same as if the line was in free space. The microstrip line was terminated using a twosection homogeneous uniaxial absorber having material parameters:=, u. in the upper section and 3.2=, Or in the lower section to match the substrate. The calculations were carried out at several frequencies using an FEM code [8] and we show the results for 4.0 GHz. At this frequency the element width was 0.05Ao and a five layer absorber having a total thickness of t = 0.25Ao was used. With a = 0 the computed reflection coefficient of the transmission line structure as a function of d3 is shown in Figure 9. Recognizing that most of the field is confined to the substrate, we have Ax = Ao/-r = 0.559A0. Using t = 1.789AX and N = 5 in the design formulas (9) and (10) gives / = 1.07 and 39

A-r21 2 cmc - - -- 3.095 cn Uniaxial absorber Microstrip (a) 0.75 1.00 1.25 1.50 1.75 2.00 (b) Figure 9: (a) Geometry for the microstrip line, (b) computed NR\ for the microstrip line at 4 GHz with a = 0, t = 0.25Ao and N = 5. The intersecting straight lines indicate the predicted values. IRI = -41 dB (indicated on the plot by the vertical and horizontal lines, respectively). These values agree reasonably well with the numerical results shown in Figure 9 where the maximum absorption occurs for 13 1 and \RI -50 dB. The fact that the minimum jRI is lower than predicted is, perhaps, not surprising. We recall that the design curves are based on the worst case, i.e. the polarization providing the largest minimum [Rj, and the curve in Figure 9 resembles more the H polarization curve in Figure 5 than the E polarization which constitutes the worst case. The other geometry considered is the rectangular waveguide shown in Figure 10. The elements are 0.5 cm bricks and the absorbing layer is 5 cm thick (10 elements used) with a material parameter b =1 -j3. The reflection coefficient was computed at 4.0, 4.5 and 5 GHz for various values of 3 in order to determine the maximum absorption point. The resulting reflection coefficients are plotted in Figure 10. Using equations (9) and (10) with N = 10 yields f3t/AX = 0.574 and \R\ =-50 dB. The vertical and the horizontal lines in Figure 10 indicate the location of these values. Once again, the agreement with the predicted values is good, with only a slight discrepancy in the value of i3t/AX and a deviation of about 5 dB in the anticipated reflection coefficient. There are two points to keep in mind here. First, the design criteria have been applied to a non-normal incidence case in a three-dimensional setting. Secondly, the real part a of the material parameter b was set to 1, 40

-20 mf l: -40- ' 2.215cm /^ 4 4 Uniaxial Absorber | -50 ----------- - ---- - --- - - - - - - --- ------------ --- 4.775cm — 0 1 2 3 4 (a) t/x (b) Figure 10: (a) Geometry for the rectangular waveguide, (b) computed JR\ for the waveguide at 4.0 GHz (- -), 4.5 GHz (.-) and 5.0 GHz ( -) with a = 1, t = 5.0 cm and N = 10. The intersecting straight lines indicate the predicted values. demonstrating that the criteria may still apply when a 4 0. 5 CONCLUSION A uniaxial perfectly matched layer provides a powerful means for truncating finite element meshes close to the modeled structure. By properly selecting the material properties and sampling rate, almost any desired level of absorption can be attained, and typically very few samples (less than five) are needed to achieve a reflection coefficient of -40 dB over a wide range of incidence angles. In this paper we described a detailed study of three types of layer material including homogeneous and inhomogeneous uniaxial ones, and by identifying the scalable parameters of the layers, universal design curves and formulas were developed. The curves or formulas can be used to specify the numerical, geometrical and electrical parameters of the PML to achieve a desired absorption down to -60 dB or lower. They are valid for all real angles of incidence for the anisotropic layer and restricted to normal incidence for the isotropic layer. As expected, a lower material loss requires a thicker absorber to produce the same reflection coefficient. On the other hand, a higher attenuation rate requires more samples to attain a lower reflection coefficient in a numerical implementation. An inhomogeneous 41

(tapered) PML is better than a homogeneous one since the material loss can be increased to larger values close to the metal backing where the field is smallest. To test the applicability of the design criteria in a three-dimensional setting, a shielded microstrip line and a rectangular waveguide were considered. With both structures terminated with a homogeneous PML, the results were in reasonable agreement with prediction, and the discrepancies were no more than could be expected in view of the conditions under which the criteria were established. These conditions are: (i) use of a particular one-dimensional FEM code (ii) based on the worst case polarization, i.e. the polarization for which the minimum reflection coefficient is largest (iii) assumption of a pure imaginary propagation constant, i.e. a = 0, in the layer. Condition (iii) is a requirement for scalability, and though small values of a are still admissible, the condition is clearly inappropriate if there is substantial power at complex angles of incidence for which large a is required for absorption. If the polarization can be specified, (ii) is also inappropriate, and the design criteria may underestimate the performance that can be achieved. In any given problem where there is the luxury of testing a variety of layer specifications, it is probable that a performance can be achieved which is better than that predicted by the criterion, but even then the design values are a logical place to start. In the more likely situation where prior testing is not feasible, we believe that the design criteria provide a logical basis for specifying the parameters of the PML and its sampling. Acknowledgment We are indebted to Ms. Lin Zhang and Mr. Jian Gong for their assistance with the calculations for the microstrip line and the rectangular waveguide. 42

References [1] T.B.A. Senior and J.L. Volakis, Approximate Boundary Conditions in Electromagnetics. Stevenage, UK, IEE Press, 1995. [2] C. Rappaport and L. Bahrmasel,"An absorbing boundary condition based on anechoic absorbers for EM scattering computation," J. Electromagn. Waves Appl. 6, pp. 1621-1634, 1992. [3] T. Ozdemir and J.L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Sci. 29, pp. 1255-1263, 1994. [4] Z.S. Sacks, D.M. Kingsland, R. Lee and J.-F. Lee, "A Perfectly Matched Anisotropic Absorber for Use as an Absorbing Boundary Condition," IEEE Trans. Antennas Propagat. 43, pp. 1460-1463, 1995. [5] J.P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comp. Phys. 114, pp. 185-200, 1994. [6] D.S. Katz, E.T. Thiele and A. Taflove, "Validation and extension to three dimensions of the Berenger PML absorbing boundary conditions for FD-TD meshes," IEEE Microwave and Guided Wave Lett. 4, pp. 268-270, 1994. [7] S.R. Legault and T.B.A. Senior, "Matched planar surfaces and layers," University of Michigan Radiation Laboratory Report, 1995. [8] J. Gong, J.L. Volakis, A.C. Woo and H.T.G. Wang, "A hybrid finite element-boundary integral method for the analysis of cavitybacked antennas of arbitrary shape," IEEE Trans. Antennas Propagat. 42, pp. 1233-1242, 1994. 43

Application and Design Guidelines of the PML Absorber for Finite Element Simulations of Microwave Packages J. Gong*, S. Legault*, Y. Botros* and J.L. Volakis* Abstract The recently introduced perfectly matched layer(PML) uniaxial absorber for frequency domain finite element simulations has several advantages. In this paper we present the application of PML for microwave circuit simulations along with design guidelines to obtain a desired level of absorption. Different feeding techniques are also investigated for improved accuracy. I. Introduction In the numerical simulation of 3D microwave circuits using partial differential approaches, it is necessary to terminate the domain with some type of non-reflective boundary conditions. When using frequency domain PDE formulations, such as the finite element method, the standard approach is to employ some type of absorbing boundary conditions(ABCs) [1], [2], [3]. Also, the use of infinite elements [4] or port conditions [5] have been investigated. All of these mesh truncation methods require a priori knowledge of the dominant mode fields and, to a great extent, their success depends on the purity of the assumed mode expansion at the mesh truncation surface. Larger computational domains must therefore be used and the accuracy of the technique in computing the scattering parameters could be compromised. Recently, a new anisotropic (uniaxial) absorber [6] was introduced for truncating finite element meshes. This absorber is reflectionless(i.e. perfectly matched at its interface) for all incident waves, regardless of their incidence angle and propagation constants. As a result, it can be placed very close to the circuit discontinuity and is particularly attractive for terminating the computational domain of high density microwave circuits where complex field distributions could be present. Although the proposed uniaxial PML absorber has a perfectly matched interface, in practice a finite metal-backed (say) layer must be used which is no longer reflectionless due *Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109-2122 44

to the presence of the pec (see Fig. 1). It is therefore of interest to optimize the absorptivity of the layer by proper selection of the parameters to achieve a given reflectivity with a minimum layer thickness. In this paper, we present guidelines for implementing the PML absorber to truncate finite element meshes in microwave circuit simulations. Example microwave circuit calculations are also given to demonstrate the accuracy of the PML absorber, the FEM simulator and the feed model. More examples will be presented at the conference. II. Absorber Design An extensive study was carried out using two-dimensional(see Fig. 1) and three dimensional models (see Fig. 2) in order to optimize the absorber's performance using the minimum thickness and discretization rate. As expected, the absorber's thickness, material properties and the discretization rate all play an equally important role on the performance of the PML. The typical field behavior interior to the absorber is shown in Fig. 3. As seen, for small 3 values the field decay is not sufficient to eliminate reflections from the metal backing. For large /3 values, the rapid decay can no longer be accurately modeled by the FEM simulation and consequently the associated VSWR increases to unacceptable values. However, an optimum value of / which minimizes the reflection coefficient for a given layer thickness and discretization can found. The parameters /3 and t play complimentary roles and the study shows that the PML absorber's performance can be characterized in terms of the product Atg (a scalable quantity when a = 0) and the discretization rate. A two-dimensional analysis was carried out to determine the optimum values of A3 and N (the number of samples in the Ag PM\4L layer) for maximum absorption near normal incidence. It was determined that given a desired reflection coefficient IRI for the PML absorber, the optimum $! and N values are Ag approximately given by the expressions [7] /t = -0.01061RI + 0.0433 g N = 0.147 exp 7.353 ] where IRI must be given in dB and N is equal to or exceeding the right hand value. As an example, if we desire to have a value of IR\ equal to -50dB, from the above formulae we have that t 3 0.58 and N = 10. It should be noted that though the design formulae were Ag derived with a = 0 they also hold for small non-zero values of a. III. Feed Excitation Two feed models were used in conjunction with the scattering parameter extraction method. One was the horizontal current probes (Fig. 2) linking the back PEC wall with the beginning of a microstrip feed line. About 3 to 5 horizontal probes were needed for convergence and this scheme proved more accurate that the usual single vertical probe. The other feeding scheme employed here involved the specification of the quasi-static TEM mode at the microstrip line port. In the context of the FEM, the excitation is introduced by imposing boundary conditions across the entire cavity cross section through the 45

input port. These conditions also serve to suppress backward reflections from the modeled circuit discontinuity. Consequently, they can be placed close to the discontinuity without compromising the accuracy of the scattering parameter extraction. IV. 3D Modeling Examples The PML performance as predicted by the formulae was investigated by using it to truncate the domain of 3-D microwave circuits. For example, Fig. 4 shows the optimum value of ft w 0.96 obtained from the above design equations compares well with the results of the full Ag wave FEM analysis of the microstrip line shown in Fig. 2. The 3-D FEM computations were carried out using N = 5 for modeling the PML absorber across its thickness and from the given formulae, it follows that R = -41dB and this agrees well with the optimum value shown in Fig. 4. Another example is the meander line shown in Fig. 5. For the FEM simulation, the structure was placed in a rectangular cavity of size 5.8mm x 18.0mm x 3.175mm. The cavity was tessellated using 29 x 150 x 5 edges and only 150 edges were used along the y-axis. The domain was terminated with a 10 layer PML, each layer being of thickness t = 0.12mm. The S11 results are shown in Fig 6 and are in good agreement with the measured data [8]. References [1] R.L. Higdon, "Absorbing boundary conditions for acoustic and elastic waves in stratified media," J. Comp. Phys., Vol. 101, pp. 386-418, 1992. [2] T.B.A. Senior and J.L. Volakis, Approximate Boundary Conditions in Electromagnetics, IEE Press, London, 1995. [3] J-S. Wang and R. Mittra, "Finite Element Analysis of MMIC structures and electronic packages using absorbing boundary conditions," IEEE Trans. Microwave Th. and Techn., Vol. 42, pp. 441-449, March 1994. [4] S. Tsitsos, A. Gibson and A. McCormick,"Higher Order Modes in Coupled Striplines: Prediction and Measurement," IEEE Trans. Microwave Th. and Techn., Vol. 42, pp. 2071-2077, Nov. 1994. [5] J.F. Lee, "Analysis of Passive Microwave Devices by Using three-dimensional tangential vector finite elements," Int. J. Num. Model.: Electronic Net., Dev. and Fields, Vol. 3, pp. 235-246, 1990. [6] Z.S. Sacks, D.M. Kingsland, R. Lee and J.F. Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition," to appear in IEEE Trans. Antennas Propagat. [7] S. Legault, T.B.A. Senior and J.L. Volakis, "Design of Planar Absorbing Layers for Domain Truncation in FEM Applications," submitted to Electromagnetics. [8] I. Wolf, "Finite Difference Time-domain Simulation of Electromagnetic Fields and Microwave Circuits," International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering, August 1992. 46

~o~rb'flt otrb,Erx rb -j: rb rb r.:...., _..... 11 Figure 1: Illustration of wave incidence upon without metal backing. a perfectly match interface (PMIL) with and Uniaxial /" Artificial Absorber Z,x Strip Width = A/10 =0.71428cm Az=0.21cm Ax —0. 1905cm Ay=O. 13636cm Cavity Cutoff Frequency=4.4GHz. Operating Frequency-4.2GHz. ( k-7.14286cm) 10 Absorber Layers - k/5 Figure 2: Shielded microstrip line terminated by a perfectly matched uniaxial al)sorber layer. 47

Figure 3: Illustration of the field decay pattern inside the PML layer. -5 -— 1=4.0GHz -— f=4.5GHz -10 f=5.0GHz a -15 -5 -30 -35 -45 — 3u 0 0.5 1 1.5 2 2.5 3 3.5 4 2pt in kg (a=1.0) Figure 4: Reflection coefficient vs 23t/Ag with a=l, for the shielded microstrip line terminated by the perfectly matched uniaxial layer. 48

I II 1.525 mm.61 m I Figure 5: Illustration of a meander line geometry used for comparison with measurement. 14 16 18 20 22 frequency (GHz) Figure 6: Comparison of calculated and measured results for the meander line shown in Fig.5. 49

Performance of an Artificial Absorber for Scattering off of a Metallic Cylinder Arik D. Brown and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2212 June 1, 1996 Abstract This paper presents a finite element solution to the scattered field off of a metallic cylinder in two dimensions. Curved rectangle edge elements are used which conform to the boundary of the cylindrical scatterer and the cylindrical, metal-backed, absorber used to terminate the FEM mesh. H-polarization was considered in order to compare the solution to that of the well known analytical solution. An error analysis was done comparing different absorbers - homogeneous isotropic, homogeneous anisotropic, and inhomogeneous anisotropic. The performance of the absorber shows that it can be used to produce accurate results for the numerical analysis of electromagnetic scattering problems. 1 Introduction An important consideration in finite element applications is the termination of the surface(2d) or volume(3d) mesh used. The method used should be efficient so that the truncation boundary can be placed close to the scatterer in order to reduce the size of the computational domain. Another important 50

performance criteria of the termination scheme used, is that it does not produce unwanted reflected waves which corrupt the numerical solution. The truncation method should be invisible to the scattered fields incident upon the truncation boundary. ABC's have typically been used for mesh termination. There are basically two types of ABC's. One involves using Wilcox's expansions [1]-[3], while the other involves a one-way wave approximation [4], [5]. A problem encountered with ABC's is their accuracy. Because of this, higher order ABC's have been developed [6], but this introduces unwanted complexity in an FEM formulation because higher order derivatives are involved. An alternative to ABC's is an artificial absorber. The absorber is used to truncate the FEM mesh therefore simulating a traveling or out-going scattered wave. Both homogeneous [7], and anisotropic [8] absorbers have been investigated previously. This paper specifically investigates using an anisotropic absorber implemented in a cylindrical coordinate system. The absorber is metal backed which implies that the tangential electric fields are zero at the absorber boundary. Throughout this paper a scattering formulation is assumed. 2 Formulation Consider an H-polarized plane wave at normal incidence ( 0 = 0~ ) on a metallic cylinder, as depicted in Figure 1. Using Galerkin's method the appropriate residual equation for the FEM system is: j VxN= f(.VxE')-k=N EdQ N. xN (.VxEE).d (1) In this equation, =i and = denote the permeability and permitivity tensors and Ee is the electric field within each element of the FEM mesh. The elements used for this particular formulation were curved rectangles because they conformed easily to the cylinder and absorber truncation boundaries. Figure 2 shows the geometry of the elements. An edge element formulation was used to approximate the electric field within each element. The field was then represented as a summation of weighting functions and unknown electric field coefficients: 51

Metal-Backed Anisotropic - Absorber i / Free Space x Figure 1: H incident on metallic cylinder N3 4o + A4 4o / N2 po + Ap x Figure 2: Geometry of FEM edge elements 52

4 E = NE;J (2) j=1 The weighting functions are represented as: N (P + AP)- P(3) Ne (3) N - (o + AO)- PO) (4) N (P-Po) (5) N = ( )() (6) After substituting equation 2 into equation 1, the resulting FEM system becomes: N N [K']E' = T[Be] (7) e=l e=l In this equation, [Ke], represents the local element coefficient matrix values, EC is the vector of unknown electric field coefficients, and [Be] represents the contributions from the boundaries of each element. The weighting functions are divergence-less, but their curl is not constant due to the non-linear shape of the elements used. The values for [Ke] are listed below: ao[(po + Ap)21n(l + As) - 2Ap2] 4po + Ap KII - k~p'aP 12 ) (8) 2Po + (9) JK = -[1- In(l+ P) (9) = ZZ PP Po Ac [ p O n+ poAp)In( Po -- k2E, A 2po + Ap 13 =, 2[(po+p o Ap)ln( ~A 0 P( ) 2Ap (10) zzP Po +?12 K4 = -Ke (11) K2, = 12 (12) 53

Ke 22 K23 Ke 24 32 33 Kie Ke2 Ke 43 K4 we Ap (2po + AP) k2_ In(l + 20+A P 3 7n(l ) ] 2pzzAq (Po + AP)2 k pp 3 Po(+ 1 2po AP Po = -[ [ n(l + k)- ] = -K22 K23 12 + 2In(l + k)]-k2,P( 1 ) PLZZ AP Po 12 - K23 = K14 = K24 - e34 K22 (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23) The tensors r and 6T which take into account the anisotropic nature of the absorber have the form: I r = er = 1 c-jp 0 0 O O 0 a-j0 0 0 0 a-jt ) (24) When computing the boundary term for the FEM formulation, the integral due to the metal backed absorber is eliminated because the tangential electric field is zero. This leaves only the boundary integral associated with the cylinder surface where: scat _ inc \ ='c ~Et' (25) Dt ang ~tang (25) The matrix system can then be solved to find the electric field values. In order to verify the numerical results, the H field was found on the surface of the cylinder by using Maxwell's equations to find the H field in each element from the computed electric field values. 4 H = J —V x E NE; WLo ~ j=l (26) 54

Homogeneous lsotrcc;c Atscrber rad =, Sarnbda. d=0 Slambda. t=C 25'am-bda 1.8 -1.6K \ -- Analytical Solution! i - Beta= 1 1.4- ' - -Beta = 2 1 1.2 -D \. 1 - \ ' 0.8 i0i I 0.6 - \ - 0.4 I Be eta = 3 I, 0/ I '\v'I 0.21 0 I __I. _, I I 20 40 60 80 100 120 140 160 180 phi Figure 3: Magnetic Field Magnitude vs. q 3 Results The results obtained using the absorber were compared to the exact analytical free space solution. The distance between the cylinder and absorber rwa chosen as, so that reflections from the absorber would not affect the solution for the scattered fields. This distance was sampled using elements with a radial length of 2 while elements in the absorber had a radial length of 2. The sampling rate around the cylinder (phi direction) was also chosen as 20 samples/wavelength. Three different types of absorber were used - homogeneous isotropic, homogeneous anisotropic, and inhomogeneous anisotropic. In each case the magnetic field was plotted versus, while a was kept constant at 1. Figure 3 shows the results for the homogeneous isotropic absorber. Figures 4 and 5 show the results for the homogeneous anisotropic and inhomogeneous anisotropic absorbers. From the graphs, it is apparent that the inhomogeneous anisotropic absorber has the best performance. Both of the other absorbers work adequately, but the inhomogeneous anisotropic absorber tracks the field variations much better. Another interesting observation is that the accuracy of the field relies heavily on the value chosen for P. From the graphs it is clear that by changing 3 simply from 1 to 3 causes a great increase in accuracy. 55

Homogen, H eous Anisotropic Absorber rad.=0.51ambda: d=0.5;ambda, t=0.251ambda 1.8 -- 1.6~-~ 1.4 I i 0.6 -12 1 -0.4 0 0 \1 I1 Analytical Solution!- - Beta = 1 - -Beta= 2 Beta =3 -.. 2604 20 40 60 80 100 120 140 160 80 phi Figure 4: Magnetic Field Magnitude vs. b Inhomogeneous Anisotropic Absorter rad =0 51ambda: dl=0.5iambda: t=0.251ambda 1.8 -1.4 Blu.= 2 1 ~ - -Beta=2 Beta = 3 1.2 -4\ 10 _31 0-6[ ' 4 ~0.28 06t 04i02 / 0 20 40 60 80 100 120 140 160 180 phi Figure 5: Magnetic Field Magnitude vs. 56

Homogeneous Isotropic Atscrte 'tad =0 Siamtda, C=0 5Sambda 50 t -t =0.15lambda - - t =.18lambda t = 0.25lambda 40- \ 35 ~30 w 25 ' 20 - 15 -0 L I 0 0.5 1 1.5 2 2.5 3 3.5 Beta Figure 6: Error vs. / In order to determine how accurate the results were, the error as a function of f? and absorber thickness was calculated. Because the actual field value is << 1, a different error criteria was used. This is because since the field value is so close to zero, the relative error is dominated by this point and gives an rinaccurate error estimate. In order to circumveInt this, the following error criteria was chosen: / actual - Hnum. Error = 100 x \W tul |12 (27) Hactu al _ where W = 1 - e- ( Hma-H ) (28) This is very similar to an RMS error with the main difference being the addition of the weighting function. In equation 28, A is a constant. Figures 6 - 8 show the error for all three types of absorber. From the data, it is obvious that the error decreases as the absorber is made thicker, and with an increase in /. 57

50 - 45 - 40O1 10~ 0 ---30 C: Hom:)geneous Anisotropic Atsorter rad.=0 S1ambda, d=0 51ambda -t 0.lsl1ambda 1 t t0.l1 lambda t=0.251lambda 0.5 1 - ______ --- ______________________________J 1. 2 2.5 3 3.5 Beta Figure 7: Error vs. /3 Inhomnogeneous Arisotopic Absorber rad, =0 51larrbda I-1 -0 S~ambda T 45 -. 35-i 0- 30~.-e 25SF 15 ~ 10 - C - 0 I I - -- I = 11. c I-. -, a I.-, -a I t = 0.251ambda i i - I 0.5 1 1.5 2 2.5 3 Beta 3. 5 Figure 8: Error vs. /3 58

4 Conclusion Conformal elements were used in order to analyze the performance of an artificial absorber termination. Three different absorbers were examined - homogeneous isotropic, homogeneous anisotropic, and inhomogeneous anisotropic. From the data, the inhomogeneous anisotropic absorber performed the best. It was also shown that by increasing the absorber width and increasing /, the error can be reduced significantly. By specifying an appropriate,3 and absorber thickness, the absorber can be used to produce accurate numerical results for other applications. References [1] C.H. Wilcox, "An expansion theorem for electromagnetic fields," Comm. Pure Appl. Math. Vol. 9, pp. 115-134, 1956 [2] R. Mittra, 0. Ramahi, A. Khebir, R. Gordon, and A. Kouhi,"A review of absorbing boundary conditions for two and three-dimensional electromagnetic scattering problems," IEEE Trans. on Magnetics Vol. 25, No. 4, pp. 3034-3039, July 1989 [3] F.R. Cooray and G.I. Costache, "An overview of the absorbing boundary conditions," Journal of Electromagnetic Waves and Applications Vol. 5, No. 10, pp. 1041-1054, 1991 [4] B. Engquist and A. Majda, "Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation," Math. Comp. Vol. 31, pp. 629-651, July 1977 [5] Constantine A. Balanis and Weimin Sun, "Vector one - way wave absorbing boundary conditions for FEM applications," IEEE Trans. Antennas Progagat. Vol. 42, No. 6, pp. 872-878, June 1994 [6] Bruno Stempfel, "Absorbing boundary conditions on arbitrary boundaries for the scalar and vector wave equations," IEEE Trans. Antennas Progagat. Vol. 42, No. 6, pp. 773-779, June 1994 59

[7] Tayfun Ozdemir and John L. Volakis, "A comparitive study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Science Vol. 29, pp. 1255-1263, Sept. 1994 [8] Jian Gong and John L. Volakis, "Performance of an artificial absorber for truncating FEM meshes," [9] Stephane R. Legault and Thomas B. A. Senior, "Circular Cylindrical Absorbing Terminations," 60

Comparison of three FMM techniques for solving hybrid FE-BI systems Sunil S. Bindiganavale and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Abstract Three different versions of the Fast Multipole Method (FMM) are employed to reduce the storage and computational requirement of the boundary integral in the finite element-boundary integral method. By virtue of its O(N1-5) or O(N1'33) operation count, the application of the FMM, results in substantial speed-up of the boundary integral portion of the code, independent of the shape of the BI contour. The main goal of the paper is to provide a comparison of the various FMM approaches on the basis of CPU time and accuracy. We present such comparisons and draw conclusions on the basis of computing the scattering from large grooves. 61

1 Introduction The finite element-boundary integral (FE-BI) method has been quite popular and extensively applied to many applications. The method [1] combines the geometrical adaptability and material generality of the FEM with the rigorous BI for truncating the mesh. Nevertheless, although "exact", the FE-BI leads to a partly full and partly sparse system and is thus coinputationally intensive for large boundary integrals. When the boundary is rectangular or circular, the FFT can be used to reduce the memory and CPU requirements down to N logN [1],[2]. However, in general, the boundary integral is not convolutional and in that case the CPU requirements will be of 0(Nb2), where Nb denotes the unknowns on the boundary. The application of the fast multipole method enables the computation of the boundary matrix-vector product with a 0(N1 5) or 0(N1 33) operation count [3],[4]. In this paper, we apply three different versions of the Fast Multipole Method (FMM) to reduce the storage and computational requirements of the boundary integral when the size of the contour is large. By virtue of its low operation count, the application of the FMM results in substantial speed-up of the boundary integral portion of the code independent of the boundary shape. 2 Problem Definition and Formulation As an application of the proposed technique, we consider the scattering by a cavity-backed aperture shown in Figure 1. The FE-BI formulation for this problem was already outlined in [2] and results in the solution of the system [BT] [P {nc} (1)The vector and typical form of this system is given in Figure 2. {q} represents the magnetic field at the nodes within the groove and on the boundary. Also, {I1 } represents the magnetic currents on the boundary. By virtue of the finite element method, the matrices [K] and [B1] are sparse and thus the corresponding matrix-vector products are implemented using 0(N) operations. However [Pi] is a full sub-matrix and thus 0(Nb) operations are needed to perform the product [P1]{41Q} with Nb denoting the number of nodes on the cavity aperture. Consequently, in an iterative solution, this matrix-vector product becomes the computational bottleneck. To reduce the operation count we will herewith employ the FMM procedure for implementing the products [P1]{01 }. Next we examine three versions of the FMM to accelerate the boundary integral matrixvector product computation. Specifically, the exact FMM [5],[6], a windowed FMM [7] and an approximate FMM [8] are examined. 3 FMM Techniques As stated above, the FMM is a fast method for calculating the matrix-vector product. The computation of the mnatrix-vector product is illustrated with the boundary integral for TE 62

J, ii"* A.! A y A pec x.4 p w Figure 1: Geometry of the groove recessed in a ground plane incidence (H-pol). In this case we deal with the discretization of the integral equation HHFEM (2) Hu: + ~ - HZ, (2) with Hz the scattered magnetic field, Hz the incident magnetic field and HFEM the total field which is employed to enforce tangential continuity and H- = 2 j~ M(p)H2 (k0o - p') d' (3) Mz - Ex is the magnetic current and F1 indicates the extent of the boundary, ~p and p' are the vectors describing the observation and source points. 3.1 Exact FMM As illustrated in Figure 3, the Nb boundary unknowns are subdivided into groups with each group assigned Mb unknowns. Thus a total of Lb m b groups are constructed. Next, by invoking the addition theorem for cylindrical wave functions, the Hankel function is expanded as Q/2 H0(2 (k0 I pI + P-p ) ( E Jn(ko l- )H2)(kopill,)ejn(tl -p/p) (4) n=-Q/2 where Oiul and pp are the angles pull and p'-p- make with the x-axis respectively and is valid for pilf > Ip' - P\. It should be noted that the source and observation point vectors, p' and p have their origin at the center of the source and test groups respectively. The semi-empirical formula Q/2 = koD + 5 ln(koD + 7r) (5) where D is the diameter of the circle enclosing the groups is used to truncate (4) and in general, Q/2 = Mb, assures convergence. The Fourier integral form of the Bessel function Jn(k -P ) = - ej2 (PP 3)do (6) 27r J27r 63

0 21I I 22:: *1 I:. 2 22 - 2L:.2.*22: 114: 2 M h.:::.Mih t~ h::: h.ct ~Ilht 20 40 60 80 Ith' *222N. 322 11s ' "::. ':!. *2::. ':.ntth. 2:1. \:1 BI System 2. t2. 2 * 100 *22h. *2. h2. *2. 120 L 0 I. I. I 20 40 60 80 100 120 column Figure 2: Typical form of the FE-BI system matrix arising from problem of a groove in a ground plane. in conjunction with (4) is used to write (3) as H() = oYo- k j ( )T,( )e-'dko 27r Ir where the far-field pattern of the source group is given by VI(X)= M ()ejko P'dl' and the translation operator Tul,()) is given by Q/2 '(~)- = E H(2)(kopll)e-jn(0-,,'+r/2) n=-Q/2 with the propagation vector of the plane wave ko = ko(x cos q + y sin q) the scattering/radiation (7) (8) (9) (10) The integral over X is discretized into Q plane wave directions, thus yielding an expression for the fields radiated by the source group 1' to the receiving group I HS(p) = - A 0 T Vl q=l (11) 64

where AO = 2,r/Q indicates the angular spacing between propagation vectors of plane waves emanating from a group and thus Oq = qA0 q 1... Q and kq = ko(x cos 5Sq + y sin Oq). As mentioned earlier, the number of plane wave directions is set equal to twice the number of elements in the group (Q = 2Mb), thus satisfying the Nyquist sampling theorem with respect to the integration over 4. f/-=o,I I. The pattern of the source group is computed. Mathematically, this corresponds toy b x d mi ), K - ( z x GuP min~, ' Far group Ps a < d tion Near group Figure 3: Computation of the boundary integral matrix vector product using exact FMM In the exact FMM, the matrix vector-product is computed in the following sequence 1. The pattern of the source group is computed. Mathematically, this corresponds to (q) M= K( 7)e 'dl' (12) Evaluating this pattern for a single source group and in a single direction requires Mab operations, corresponding to the number of elements in the group (the integration over the line segment is performed as a summation). Consequently for Lb groups and Q directions for each group, the operation count is QMbLb. 2. The translation operation is employed to evaluate the pattern of a source group at a test group. Mathematically, AI(4q) = V'1(/q)Tr11(4'q) (13) For Q directions and Lb source and test groups, this process involves an operation count of QL. 3. At the receiving group, the fields are redistributed. Mathematically, Hk(P) = - A Aj(q)ei- (14) q=1 Thus, computing the redistributed field at a single point requires Q operations, and for Lb groups each containing Mb unknowns the operation count is QLbMb. 65

The total operation count in the above sequence is CQMbLb+ C2QLb. By choosing, Q ~ ib for convergence, the operation count is given by ClMbNb + C2 N-. and on setting Mb ~ vNb, the final operation count is NJ'5. Improvements to the operation count can be achieved by nesting groups leading to the multi-level FMM. Beyond the math, the above breakdown of operations are based on the manager-worker model. Basically, we can view each group as managed by the center element with the workers comprising the elements of the group. Communication/interaction among the groups takes place through the managers which in turn interact with the group elements. However, this type of model is based on certain simplification/decompositions of the original boundary integral. Clearly, they reduce the direct interdependence of each group member with other elements belonging to different groups. This is the essence of the CPU speed up afforded by FMM. However there are inherent approximations which must be understood in order to assess the accuracy of each FMM algorithm. 3.2 Windowed FMM In the exact FMM, the translation operation between groups was considered isotropically. But, it is suggestive that the groups would interact very strongly along the line joining them and less strongly in other directions. It was shown that by employing a high frequency analysis [7], that the translation operator could be contemplated as being composed of a geometrical optics term, along the line joining the source and test group, and two uniform diffraction terms associated with the shadow boundaries of the GO term. The translation operator for different group separation distances along the groove of width 50A is shown in Figure 4. The number of unknowns on the boundary are 750, resulting in 27 groups. It is seen that the "lit" region of the translation operator decreases with increasing group separation distance, eventually displaying the predictable sine function behavior for large group separation distances. The high frequency analysis enables the identification of a lit region even for groups which are not very widely separated. Thus the plane wave interactions can be filtered by defining a filter function as WV1iQ q) { I (I( q -,1 I| < < s) = -cC(Jqq-0qt' [-_s)2 (~q- #I/1 > /3s) (15) with =sin- ) (16) 2kopu, and a is the taper factor. The discretized plane wave expansion is now - koY () (17) q=1 The operation count now associated with (13) is now reduced to C3L2,_ Mb with the corresponding total operation count given by ClMbNb + C4. Grouping the unknowns with b 66

N1/3 per group, results in an operation count of O(N4/3). The computation of the boundary integral matrix vector product by employing the windowed FMM is depicted pictorially in Figure 5. It is seen that the filter has the effect of eliminating plane wave interactions at directions away from the line joining the interacting groups. The spectrum of interactions around the line joining the centers of the two groups which are retained, reduces with increasing group separation distance. 3.00 I-.- -, Groups 1 and 3 250 I=s4;- -= 5 ----- Q54;P =15Xk - 2.00 ff,,, =Qs4;P,/= 45 x, T G V oups 1 and 8 Groups 1 an 27 1.00 -0I I I 3.3 Fast Far Field Algorithm The third algorithm employed for hybridization was the fast far field algorithm (FAFFA) which is an approximate version of the FMM. Unlike the exact FMM, where the kernel is approximated with tile addition theorem, in this algorithm the large argument approximation ) (2) /9j ko_-jkopl -j -. H8 (ko P- P' e-J oP m e OP Pnil (18') " V 2k'opl is used, where pTr is the distance between the center of the test group I and the center of the source group r; Pnlo is the distance between the nth source element and its group center and pim is the distance between the mth test element and its group center, as depicted in Figure 6. It should be noted that the use of the large argument expansion, as an additional 67

l I Group j, Group,'Group/Group KI,, 4' K- 1 i N | Window ' —" ' ---; ' ---" X 'l I> dmin~ | I Far group P J< dmin= Near group Figure 5: Computation of the boundary integral matrix vector product using windowed FMM approximation, necessitates that the FMM procedure can now be used for groups which are only very well separated. The decoupling of the test-source element interactions in the kernel as in (18) enables the computation of the matrix-vector product for far-field groups with a reduced operation count as detailed in the following sequence. 1. For each test group, the aggregation of source elements in a single source group involves Mb operations, corresponding to the number of elements in the source group. The aggregation operation corresponds to Mb V1 = M -i (19) j=l 2. Since the above aggregation operation needs to be done for all source groups the operation count becomes O(N — Mb) - O(Nb), where -Nb represents the total number of groups. Also this operation, being dependent only on the test group rather than the test element, needs to be repeated for Nb test groups leading to a total operation count of O(-b) for aggregation. 3. The next step is a translation operation corresponding to Alli = TllVl'l (20) where Ti l — o (21) V /kopIl,i Since this needs to be done only at the group level, it involves O(N) operations for all possible test and source group combinations and is the least computationally intensive step. 68

l I I'P, Group Group,'Group Group K K + K - i/' I '__ _ dmin, ' \ I / I V \ / I// ---, _ / I I I I n!I- P lPIm m > dmin L Far group | P/|l< d in Near group Figure 6: Computation of the boundary integral matrix vector product using the FAFFA 4. The final step in the sequence would be the process of disaggregation corresponding to the operation o yo NbI Mb H (p) --- ~Y Aiie-jk~O Plm (22) ~ / 2 2 =1 Conceptually, this process is the converse of aggregation. Since this operation involves only the source group instead of the source element it needs to be done for each source group thus implying an O(N_ ) operation to generate a single row of the matrix-vector product. To generate Mb rows corresponding to a test group the operation count would be O(Nb). With -N test groups, the operation count would be of O(- ). N2 5. The near field operation count being of O(NbMb) and the far field being O(M) gives a total operation count of N2 Op.count ClNbMb + C2 (23) Mlb Typically, we can set the elements in each group, Mb = \NNb and as a result the total operation count is O~ NbJ5. 4 Results A computer code based on the above formulation was implemented and executed on a HP 9000/750 workstation with a peak flop rate of 23.7 MFLOPS. The geometry considered was the rectangular groove shown in Figure 1. The performance of the hybrid algorithms with respect to accuracy and speed were compared. The benchmark for accuracy was the RMS error [9]. Table 1 compares the execution time and error of the standard FE-BI to the FE-Exact FMM, FE-FAFFA and the FE-Windowed FMM for grooves of widths 25A, 35A 69

and 50A. The depth of the groove was 0.35A with a material filling of Or = 4 and,r=l and was illuminated at normal incidence. The data reveals that the FE-Exact FMM offers almost a 50% savings in execution time with almost no compromise in accuracy. While the FE-FAFFA is the fastest of the three algorithms, the RMS error was substantially higher (> I dB). If the maximum tolerable RMS error is set at 1 dB [9], the FE-Windowed FMNM is the most attractive option as it meets the error criterion while being only slightly slower than the FE-FAFFA. An observation of interest was the comparison of the residual error as CPU Time for BI (Minutes,seconds) Groove Width Total Unknowns BI Unknowns FE-BI (CG) FE-FAFFA FE-Exact FMM FE-WFMT 25A 2631 375 (8,48) (3,26) (5,25) (4,13) 35A 3681 525 (16,34) (5,55) (10,31) (7,22) 50A 5256 750 (45,1) (14,31) (26,18) (16,10)........ RMS error (dB) Groove Width FE-FAFFA FE-Exact FMM FE-WFMM 25A 1.12 0.0752 0.6218 35A 1.2 0.1058 0.721 50A 1.36 0.1123 0.843 Table 1: CPU Times and RMS error of the hybrid algorithms a function of the number of iterations in the CG solver for the hybrid algorithms. This is shown in Figure 7. It is seen that the convergence curves for the FE-BI and the FE-Exact FMM overlap each other to graphical accuracy while the FE-Windowed FMM shows a very small deviation. Thus, the hybridization of the FMM does not have a deleterious effect on the conditioning of the FE-BI system. The execution time of the fastest of the three hybrid techniques, the FE-FAFFA is conmpared with the FE-BI employing the special CG-FFT solver, suited for only planar apertures in Table 2. It is seen that the CG-FFT solver is substantially faster but is applicable only to convolutional boundary integrals. The performance of the hybrid algorithms at a more stressing angle of incidence is depicted in Figure 8. The RMS error follows the same trend as for normal incidence illumination. The width of the groove illuminated was 10A and this example reveals that the technique is scalable for smaller problems. The near group radius in this example was 1A; implying that the matrix vector products for groups separated by a distance less than a wavelength was computed using the exact method of moments procedure. It has to be noted that the application of the hybrid techniques for the 10A groove illustrates that the neargroup radius can be reduced as the problem size becomes smaller down to a minimum in the vicinity of 0.3A. 70

100 10-1 10 10-3 0 100 200 300 400 500 600 Iteration number Figure 7: Convergence curves for the hybrid algorithms for the groove of width 25A 4.1 Summary The hybridization of the finite element - boundary integral and fast multipole methods was examined from a computational performance point of view. Three different versions of the fast multipole method was used for executing the matrix-vector product associated with the boundary integral. It was shown that a considerable reduction in CPU time could be achieved and further reductions are likely as the boundary integral increases in size. The FE-Windowed FMM provides the best compromise in terms of speed and accuracy. The FE-BI with the CG-FFT solver is faster than any of the FEM-FMM versions and can be used if the terminating boundary is amenable to its use. References [1] J.L. Volakis, A. Chatterjee, and J. Gong. A class of hybrid finite element methods for electromagnetics: A review. Journal of Electromagnetic Waves and Applications, 8(9/10):1095-1124, 1994. [2] J.M. Jin and J.L. Volakis. TE scattering by an inhomogeneously filled aperture in a thick conducting plane. IEEE Transactions on Antennas and Propagation, 38:1280-1286. 71

CPU Time for BI (Minutes, seconds) Groove Width Total Unknowns BI Unknowns FE-BI (CG) FE-FAFFA FE-BI (CGFFT) 25A 2631 375 (8,48) (3,26) (1,41) 35A 3681 525 (16,34) (5,55) - (3,24) 50A 5256 750 (45,1) (14,31) (5,40) Table 2: CPU Times for the FE-BI (CGFFT) compared to the FE-FAFFA August 1990. [3] S. S. Bindiganavale and J. L. Volakis. A hybrid FEM-FMM technique for electromagnetic scattering. In 12th Annual Review of Progress in Applied Computational Electromagnetics (ACES) Digest, Naval Postgraduate school, Monterey CA, pages 563-570, March 1996. [4] N. Lu and J. MI. Jin. Application of the fast multipole method to finite element-boundary integral solution of scattering problems. In 12th Annual Review of Progress in Applied Computational Electromagnetics (ACES) Digest, Naval Postgraduate school, Monterey CA, pages 1182-1189, March 1996. [5] V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions. J. Comput. Phys., 86(2):414-439, 1990. [6] R. Coifman, V. Rokhlin, and S. Wandzura. The fast multipole method for the wave equation: A pedestrian prescription. IEEE Antennas and Propagation Magazine, 35(3):7 -12, 1993. [7] R.J. Burkholder and D.H. Kwon. High-frequency asymptotic acceleration of the fast multipole method. Radio Science, To appear in 1996. [8] C.C. Lu and W.C. Chew. Fast far field approximation for calculating the RCS of large objects. Micro. Opt. Tech. Lett., 8(5):238-241, April 1995. [9] S. S. Bindiganavale and J. L. Volakis. Guidelines for using the fast multipole method to calculate the RCS of large objects. Micro. Opt. Tech. Lett., 11(4):190-194, March 1996. 72

30 20 P3. 10 C-10 ci -10. w= lO (a) -20 -30 0 30 60 90 120 150 180 Observation angle (deg.) (b) RMS error (dB) Groove Width FE-FAFFA FE-Exact FMM FE-WFMM 10A 0.7627 0.1621 0.3291 I (c) Figure 8: Scalability of the hybrid techniques to smaller problems (a) Problem geometry (b) Bistatic patterns (c) Error table 73

Hybrid Finite Element Methodologies for Antennas and Scattering J. L. Volakis, T. Ozdemir and J. Gong Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, Michigan 48109-2122 Abstract This paper is an overview of the finite element method (FEM) as applied to electromagnetic scattering and radiation problems. A brief review of the methodology is given with particular emphasis on new developments over the past five years relating to feed modeling, parallelization and mesh truncation. New applications which illustrate the method's capabilities, versatility and utility for general purpose application are discussed. 1 Introduction Over the past 10 years we have witnessed an increasing reliance on computational methods for the characterization of electromagnetic problems. Although traditional integral equation methods continue to be used for many applications, one can safely state that in recent years the greatest progress in computational electromagnetics has been in the development and application of partial differential equation(PDE) methods such as the finite difference-time domain and finite element (FEM) methods, including hybridizations of these with integral equation and high frequency techniques. The major reasons for the increasing reliance on PDE methods stem from their inherent geometrical adaptability, low O(N) memory demand and their capability to model heterogeneous (isotropic or anisotropic) geometries. These attributes are essential in developing general-purpose analysis/design codes for electromagnetic scattering, antennas, microwave circuits and biomedical applications. For example, modern aircraft geometries contain large nonmetallic or composite material sections and antenna geometries may involve many layers of materials, including complex microwave circuit feeding networks. Although, the moment method continues to remain the most accurate and efficient approach for subwavelength size bodies of simple geometries, PDE methods and hybrid versions of these 74

have shown a much greater promise for large scale simulations without placing restrictions on the geometry and material composition of the structure. This paper is a selected review of hybrid finite element methods and their applications to scattering and antenna analysis. We begin by first introducing the mathematical basics of the method without reference to a specific applications and in a manner that identifies the method's inherent advantages in handling arbitrary geometrical configurations which may incorporate impedance boundary conditions and anisotropic materials. In this section we also identify the key challenges associated with the implementation of the method such as mesh truncation and feed modeling for antenna applications. Section 3 of the paper reviews the various mesh truncation schemes to be employed on the outer mesh surface for the unique solution of the vector wave equation. Absorbing boundary conditions, integral equations and artificial absorbers are discussed, all leading to different versions of hybrid FEM methodologies, and we comment on their accuracy and ease of implementation. Section 4 presents several approaches for antenna feed modeling in the context of the FEM, including coaxial as well as aperture coupled feeds. The next section is devoted to parallelization issues specific to finite element codes. We give performances of typical FEM codes and provide storage and implementation guidelines for maximizing code performance on parallel computing platforms. The final section of the paper gives some representative applications of the method to scattering and antenna problems. 2 Theory 2.1 FEM Formulation Consider the antenna and scattering configurations shown in Figure 1. In the case of a scatterer, the entire computational domain is enclosed by a fictitious surface SO that may encompass a variety of composite/dielectric volumes as well as metallic, impedance and resistive surfaces. The antenna geometry is assumed to be recessed in some doubly curved surface. In this case, the bounding surface SO may either be located as shown or can be coincident with the antenna aperture. As usual, the recessed cavity is intended to house the radiating elements and their feeding structure such as coaxial cables, striplines, microstrip lines or aperture coupled feeds. The cavity may encompass any inhomogeneous or anisotropic material including resistive cards, lumped or distributed loads and so on. The goal with any finite element formulation is to obtain the solution of the vector wave equation V[x (v (1) V X [l * (V x E)] -. E =-jkoZoJi- V x (1 M) (1) in which E represents the total electric field, (e m.) denote the relative permittivity and permeability tensors of the computational domain, ko is the free space wavenumber and Zo is the free space intrinsic impedance. In addition, Ji and M- denote the impressed electric and magnetic sources, respectively, and represent the excitation for the antenna problem. As is well known, the standard finite element (FE) solution scheme is to 75

Conformal Patch Array So Surface So enclosing the computational domain - i - Resistive or ADiet eleof impedance S Pt s............... (enclosed by Sdhicle dd Platform line Absorber termination (c) Figure 1: ComputaResistiveonal domains for FEM analysis, (a) the various regions enclosed by So, (b) typical tetrahedral mesh, (c) computational regions for antenna analysis consider the weak form of (1), instead of solving it directly. This can be achieved by extremizing the impe dartinent functional [1] + IIL E [jkoZoJi + V x (U. M4)Idv + ikoZo JJ E* (H x n)ds So +Sf + ljkoZo I]I (n x E)* (n x E)ds + d |I E EdVj (2) volwher ume dnts t Vd re surface S S stands or the junction opening to possible guided feeding structures and H is the corresponding magnetic field on So and Sf whose outer normal is given by n. Also, Vs is the source volume and Vi is the volume occupied by the load Z, whose length and crby d and (a) line Absorber termination (c) s, respectively. It is noted 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 eresistive she pertinent functional [2] r(r) = I (V x r). [%-1. (V x r)]- ( ~ E). E} dv x E x x E = -R x (H +- -) where Hi denotes the fields on the two sides of the surface Sr. As seen from (2), in spite of the discontinuity in the magnetic field, no special care is correquired for the evaluation of of the discontinuity in the magnetic field, no special care is required for the evaluation of 76

the integral over S,. The explicit knowledge of H is, however, required over the surface So and Sf (referred to as mesh truncation surfaces) for the unique solution of E. As an alternative to (2), we can instead begin with (1) by weighting it with a testing function T. Subsequently, application of the divergence theorem yields the residual < RT > = JJ {(VxT) r (V x E)-k Er ET} dV + JJJI T {jkoZoJ, + V x (-1 * M)} dV + jkoZo ff T (H x n)dS So+Sf + j Zo {J (n x T). (n x E)dS + ZS TJ. EdV' (3) This is typically referred to as the weak form of the wave equation, whereas (2) is referred to as the variational form. It will be seen later that when set equal to zero,(3) yields the same system of equations as those deduced from the functional (2). Therefore both methods are equivalent. When dealing with large computational domains as is often the case in scattering applications, it is instructive to work with the scattered rather than the total field, and also distinguish the air and dielectric regions. This approach results in reduced errors within the computational domain and to proceed we introduce the definitions E = Einc + Esat H H=in + Hscat (4) where Eic denotes the excitation field impinging from the exterior of So and Escat is the scattered field. With these substitutions, the functional F takes the form F(E ) = 2 JIVV {(V x Escat) (V x Est) - kEscat E } dV + fJ/ Ecat' [jkoZoJi + V x (.-1 Mi)] dV + jkoZo If Escat (H x n) dS +Sf + 5jkoZo { J (n x E ) (n x ) Est) dS + zd Ik E. E scdV + jkoZo { (A x Esat). (A x Ei) dS + zL Il E~. EincdV} + 1 IVd {v x Escat (1 V X Einc) - k0r' Einc Escat} dV (5) + 1 d {V x Einc (71 * V x Escat) - Escat Enc} dV + f(Einc) in which f (Einc) is a collection of terms involving only the incident field Einc. Its specific form is of no consequence since it will be eliminated in the subsequent minimization of F. Most importantly, we observe in (6) the presence of integrals over the volume Vd 77

occupied by the dielectrics and over the surface Sr. The introduction of these integrals is required when working with the scattered field because Escat is discontinuous at boundaries separating media with different permeabilities [4]. Also note that the presence of the last two integrals is necessitated because EiTnc satisfies the wave equation only in free space, and that these two volume integrals become equal when the dielectric anisotropy reduces to isotropy. The functionals (2) and (6) or the weighted residual (3) can be cast into a discrete system of equations for the solution of E. To accomplish this, it is first necessary to subdivide the volume as a collection of small elements such as those shown in Figure 2 [3]. Within each volume element, the field can then be expanded as Right Angled Brick Skewed Brick Curvilinear Brick Tetrahedron Distorted Prism Cylindrical Shell Figure 2: Illustration of the various elements for tessellating three dimensional volumes M=#of edges E = E~w = [We ]{Ee} k=l1 (6) in which We are referred to as the edge-based shape or basis functions of the eth element in the computational domain. In contrast to the traditional node-based shape functions, the edge-based shape functions are better suited [5] for simulating three dimensional electromagnetic fields at corners and edges. Moreover, edge-based shape functions overcome difficulties associated with spurious resonances [6]. They were proposed by Whitney [7] over 35 years ago and were revived in the 1980s [8],[9]. Their representation is different for each element but have the common properties [10] of being divergence free (i.e., V e Wk = 0 within the element) and normalized in such a manner so that the expansion coefficients Ek represent the average field value across the kth edge of the eth element. 78

One disadvantage of the edge-based elements is that they increase the unknown count. However, this is balanced by the increased sparsity of the resulting stiffness matrix. A detailed mathematical presentation of the edge-based shape functions for various two and three dimensional elements (bricks, hexahedra and tetrahedra) [11]-[14] is given by Graglia et. al. [15] in this issue of the Transactions. Linear as well as higher order (curvilinear) elements are discussed [16]-[19]. 2.2 Discretization and System Assembly Once the computational domain has been tessellated using appropriate elements and shape functions, the discretization of (2), (3) or (6) proceeds by introducing the expansion E or Ecat = EINe Ee, with E as given by (6) and Ne denoting the number of elements comprising the computational domain. For the functional (2), the system of equations is constructed by setting (Rayleigh-Ritz procedure) OF(E) 0 =O e=l,2,..., e; k=,2,...M (7) whereas in the case of (3) Galerkin's approach is used by setting T = We. Regardless of the procedure, the resulting system will be of the form N, N, Ne Np [Ae]{Ee} + 1[Bs]{ES} + Ke} + {CP} = O (8) e=l s=1 e=l p=l in which the brackets denote square matrices and the curly braces refer to columns. Among the various new parameters, [Ae] is referred to as the element matrix and results from the discretization of the first volume integral in (2),(3) or (6); N, is equal to the aggregate of the surface elements (quadrilaterals for hexahedra or triangles for tetrahedra and prisms) used for the tessellation of S, Sf and Sr; the column {Ke} results from the discretization of the source integral over Vs; and [Bs] is the matrix associated with the third and fourth surface integrals in (2),(3) or (6). Finally, Np denotes the aggregate of the surface and volume elements over Sr, V and Vd and the column {CP} results from the discretization of the integrals in (6) involving the external incident field Ei. Basically, {CP} provides the external excitations (scattering problem) whereas {Ke} is the corresponding excitation column for internal sources as is the case with antenna problems. The entries of the element matrix [Ae] are specific to the choice of the shape functions and are. compactly given by [AK] = Ifv {[DWe][=]-[DWe]T - k~ [We][r][We]} dV (9) with [DW ]= ax {We]-T (10) a {We} T a {We}T ')x y a ~ (lo 79

and Ve denotes the volume of the discrete element while the subscripts (x, y, z) in (10) imply the (x, y, z) components of the shape function We. If lumped loads are present (i.e., in the presence of a volume integral over Vie), the diagonal entries of [Ae] are simply modified with the addition of the term jkoZod2/ZL. Perfectly conducting posts located at the kth edge are handled by setting the corresponding total fields equal to zero and rearranging the final matrix as discussed later. The internal source excitation column is simply given by Ji x Mi x Ke = JJJ [W] {koZo JiY + V x [(r-1r { AMy }) dV (11) and the entries of the corresponding excitation column due to the incident field with the presence of dielectric materials are Ex inc {CP} - jkoZo [P] E dS i n T + - (vx { V ) [r] [DWP] - o[Wp][r] dV z ) ( Ez in which WP implies the representation of the shape function over the required pth surface or volume element. The specification of the matrix [B5] can not be completed without first introducing some a priori relationship between the E and H fields on So. This relationship (or boundary condition on So) is referred to as the mesh termination condition and its form depends on the physics of the problem. For example, in the case of problems where the entire computational domain is enclosed by a perfect electric or magnetic conductor, the unknown column {Es} is eliminated from {E}. For scattering and radiation problems, the mesh termination condition should be such that the artificial boundary So appears transparent to all waves incident from the interior of the computational domain. Clearly, if So is placed far enough from the scatterer or radiator, the simple Sommerfeld radiation condition provides the appropriate relationship between E and H. However, this is not practical since it will yield large computational domains. To bring So as close as possible to the scattering or radiating surface, more sophisticated boundary conditions must be introduced on So. These mesh termination conditions are critical to the accuracy and efficiency of the formulation and are some of the major bottlenecks in the implementation of the FEM. Regardless of the employed mesh termination approach, after carrying out the sums in (8), a process referred to as matrix assembly, the resulting matrix system 80

will be of the general form,4, <f {EV b [o] [o] 1 J {Ev} | - J {b \ ^ [A] {EB} L[8 ] a l {E} {bv } (12) In this, {EV} denotes the field unknowns within the volume enclosed by So + Sf whereas {E } represents the corresponding unknowns on the boundaries So and/or Sf. The excitation column {b'T} results from the assembly of {Ke} and similarly {bB} is associated with {CP}. The matrix [A] is very sparse (9 to 30 non-zero entries per row) and this is a major characteristic and advantage of FEM. By using special storage schemes and solvers suitable for sparse systems, the CPU and memory requirements are maintained at 0(N) and scalability can be attained on multiprocessor platforms. The boundary matrix [9] in (12) is associated with the boundary fields {EB} and its specific form is determined by the employed mesh termination condition on SO. Over the past five years, much research on FEM has concentrated on the development of mesh truncation schemes which minimize the computational burden without compromising accuracy. As will be seen later, this compromise is difficult to attain and is subject to the computational problem being considered. For the purpose of our discussion, we will subdivide the various mesh termination schemes into three categories: (1) exact boundary conditions, (2) absorbing boundary conditions (or ABCs), and (3) artificial absorbers. Exact boundary conditions provide an integral relation between the electric and magnetic fields, and because they are exact, they permit the placement of SO very near or exactly on the surface of the scatterer or radiator. The resulting formulation is referred to as the finite element-boundary integral (FE-BI) method and combines the adaptability of the FEM and the rigor of the boundary integral methods. However, it yields a fully populated matrix [9] and as a result, the FE-BI method is associated with higher computational demands even though only a small part of the overall system is full (see Figure 3a). To alleviate the higher computational demands of the rigorous FE-BI method, ABCs or artificial absorbers (AAs) can instead be used to terminate the mesh. Both of these are approximate mesh termination methods, but lead to completely sparse and scalable systems of equations. Basically, the sub-matrix [G] is eliminated with the application of ABCs or AAs (see Figure 3b). In the case of ABCs, a local boundary condition in the form of a differential equation is applied on SO to relate E and H so that SO appears as transparent as possible to the incident fields from the interior. The resulting method will be referred to as the finite element-ABC (FE-ABC) method and has so far been the primary method for general-purpose scattering computations. Finally, the use of artificial absorbers (including perfectly matched layers or PMLs) [20] have recently gained major attention because of their potential for greater accuracy and inherent implementation simplicity. In the context of the finite element-artificial absorber (FE-AA) method, the mesh is terminated by using a material absorber (typically non-realizable in practice) to absorb the outgoing waves and suppress backward reflections. Below, we briefly discuss the specifics for each of the three mesh termination schemes. 81

0 20 40 I=60 80 I I Iz. 2. — - U2 Ifu. 2:.: 2: 22::::. System 100 ~ *222:2 1201 0 X 104 0 0.5 1.5 2 2.5 IL__ _ 20 40 60 column 80 100 120 1.5 nz = 469151 X10 Figure 3: Examples of matrix systems generated by the finite element-boundary integral (top figure) and finite element-ABC methods (lower figure) 82

3 Mesh Termination Schemes 3.1 Finite Element-Boundary Integral Method The FE-BI method was introduced in the early seventies [21],[22] as a natural extension of the FEM for modeling unbounded problems. However, because of its larger computational requirements, the method did not enjoy a widespread application to electromagnetics until the late eighties [23, 24]. In accordance with the FE-BI method, the relation between E and H on So is provided by the Stratton-Chu integral equation (or its dual) n x H(r) = n x Hi'(r) + n x JJ {[n' x H(r)] x V'Go(r, r') + jkYo[,n' x E(r)]Go(r r') (13) + o V' [nt' E(r')] VGO(r, r') dS' where r and r' denote the observation and integration points, respectively, and Go (r, r') = exp(-jk lr - r'l/(47rr - r'l) is the free-space Green's function. The above is the most general form of the boundary integral and places no restrictions on the shape of So with the exception that it must be closed as shown in Figure l(a). Provided So is placed just above the outer boundary of the scatterer or radiator, any material composition which is interior to So can be handled with ease using the FEM without relation to the boundary integral. This form of the FE-BI was used by Yuan [25] and later by Angelini et. al. [26], Antilla and Alexopoulos [27], and Soudais [28] to model three dimensional scatterers with anisotropic treatments. The method has also been successfully used for biomedical simulations [29],[30]. The discretization of (14) follows by introducing the expansion (over the sth surface element.) n=# of edges E(r) EESs(R)= [S]T {Es} (14) k=l with a similar expansion for the H field (if necessary). The surface shape functions Ss(r) can be set equal to We(r) when the latter are evaluated on the surface of the volume element. For the tetrahedron, fn x Sk are simply equal to the traditional basis functions given by Rao et. al. [31]. However, as pointed out in [32], Ss can be chosen indepen'dently of We provided care is exercised when (14) is substituted/combined with (2)-(6). Regardless, when (14) is substituted into (3) or (6) after discretization we will get the partly sparse, partly full system [33] given in (12) and illustrated in Figure 3(a). Finally, we point out that since So is a closed surface, the final system is not immune to false interior resonances due to the boundary integral. In this case, the combined field formulation [34] or the simpler complexification scheme [35] may be employed. When So coincides with an aperture in a ground plane (see Figure l(b)), the integral equation (14) can be simplified substantially [36]. Specifically the integral relation between E and H on the aperture reduces to n x H = n x H~ + J n x G(r, r') [E(r) x n] dS (15) wiSot 83

where G is the dyadic Green's function of the second kind, with n x V x G -0 on So and for planar platforms it reduces to, G -\ i-+ jVVr-r (16) K+t vV) 47rlr - r'l In this, I is the unit dyad and the factor of 2 is due to image theory. Also, Hg" is equal to the sum of the incident and reflected fields (in the absence of the aperture) for scattering or zero for antenna parameter computations. For non-planar So, Hg" is equal to the field scattered by an otherwise smooth surface again with the aperture removed. In this case, the Green's function should also be modified accordingly with respect to the non-planar platform. One of the first implementations of the FE-BI for radiation and scattering from rectangular apertures/antenna recessed in a ground plane was given by Jin and Volakis [14],[33],[36] and was later extended to antenna analysis on planar [37], [38], [39] and cylindrical [40] platforms. The FE-BI method is particularly attractive in terms of CPU and memory requirements when the [9] matrix is Toeplitz and can therefore be cast in circulant form [41],[42], [43], [44]. In this case, the entire system can be solved using an iterative solver [45] in conjunction with the FFT to reduce the CPU requirements down to O(Nselog Nse) for the boundary integral sub-matrix. The FFT is simply used to carry out the matrix-vector product [G] {EB} within the iterative solver. For example, if the symmetric biconjugate gradient (BCG) method [45] is employed to solve (12) and rectangular elements are used, the storage requirement is only 4Nee + 8Nse, where Nee and Nse denote the number of edges within the computational volume and on SO, respectively. For triangular surface elements, the storage requirement is about 4.5 times larger due to the presence of the diagonal edges across the quadrilaterals. Whether the full matrix [G] is Toeplitz or not depends on the shape of the surface So and on the uniformity of the mesh. It can be shown that [9] will be Toeplitz for planar and cylindrical surfaces provided the surface grid is uniform. The Toeplitz form of [9] has been exploited with great success and systems with more than 0.5 million unknowns spanning apertures of 15A x 15A have been solved accurately on a desktop workstation. Convergence of the iterative solver is also quite impressive for these systems and as demonstrated by Ozdemir and Volakis [46] (see Figure 4) the FE-BI method is the method of choice for these applications. For example, we observed that a 180,000 unknowns system converged in 57 iterations and for another system of 25,000 unknowns convergence was achieved after 66 iterations with an average CPU time of 2 sec/iteration on an HP9000/750 rated at 24.7 MFlops. Because of these impressive performances, it is advantageous to transform [9] to a Toeplitz matrix even when the surface grid is not uniform. To do so, a uniform grid can be superimposed onto the non-uniform mesh (see Figure 5) and the edge fields on the two grids can be related via a connectivity matrix [39]. In this manner, non-rectangular antenna elements and apertures can be modeled with the full geometrical adaptability of the finite element method and without compromise in accuracy. It is important to note that similar transformation matrices can be exploited for decomposing computational domains as done, for example, in [47]. When So is not planar, the boundary integral matrix will inevitably cause large CPU and storage requirements for large Nse to the point where FE-ABC or FE-AA methods 84

100000 X X r w * X * * X w w 0 u 0 cn 0-. 0 a.) 10000 1000 100 i1111 /,l / > FE-BI --- - FE-ABC - FE-AA.................. 0 5 10 15 20 Number of Array Elements in x direction Figure 4: Comparison FE-ABC and FE-AA of the convergence behavior/CPU requirements among the FE-BI, methods In:. 1 2.3 Figure 5: Overlay of a uniform grid on an unstructured mesh for implementation of the FFT to perform the matrix-vector products. 85

become more attractive at the expense of accuracy. Recently, though, techniques have been proposed which show promise in reducing the computational requirements of the matrix-vector product [G] {EB}. Among them, the fast multipole method [48], [49] can reduce the operation count in carrying out the matrix-vector product from 0 (Nt2~) down to 0 (NA'5) and reductions down to 0 (NA133) have also been proposed. The main idea of the FMM [49] is to subdivide the surface SO into groups, each containing approximately Mse VNse unknowns. By rewriting the free space Green's function as an expansion [50] or by introducing its far field approximation [51], it can be shown that the interactions between the unknowns within the groups can be carried out in 0 (NseMse) operations whereas the interactions between groups can be carried out in 0 (Nl2/Mse) operations. The sum of these two operation counts yields a total operation count of 0 (NAe5) for the boundary integral matrix when Mse - NAe, and this must be added to the 0 (Nee) operation count associated with the sparse matrix [A]. Details for the implementation of the FMM are given in [48], [49] and when combined with the FEM, we can refer to it as the FE-FMM method [52]. A comparison of the FE-BI (with LU and conjugate gradient-FFT solvers) and the FE-FMM methods for the calculation of the scattering by a large groove in a ground plane is tabulated in Table 1. The groove was 50A long CPU Time for BI (Minutes, seconds) Groove Width Total Unknowns BI Unknowns FE-BI (CG) FE-FMM FE-BI (CGFFT) 25A 2631 375 (8,48) (3,26) (1,41) 35A - 3681 - 525 (16,34) (5,55) (3,24) 50A 5256 750 (45,1) (14,31) (5,40) Storage of BI (KB) Groove Width FE-BI (CG) FE-FMM 25A 1072 277 35A 2102 -458 50A 4291 784 RMS error (dB) Groove Width FE-FMM 25A 1.12 35A 1.2 50A 1.36 Table 1: CPU Times, storage requirement and error and 0.35A deep (all metallic) and filled with a dielectric having Cr = 4 and y, = 1. By using a sampling rate of about 15 edges per wavelength, the system unknowns were Nee = 5256 and Ns = 750. From the data in Table 1, it is clear that the FE-FMM is more than 3 times faster and requires 3-5 times less memory than the FE-BI when LU decomposition is used for solving the FE-BI system. Because of the grouping/averaging of the surface elements, the FE-FMM exhibits certain errors which are functions of many parameters [53] and these errors must be kept in mind as the system size is increased. It is expected that the FMM will exhibit greater speed-ups as the system size is increased. Nevertheless, the special implementation of the FE-BI using the FFT is still by far the most accurate approach. We close this section by noting that another approach for treating the boundary integral matrix-vector products is via the adaptive integral method (AIM) [54]. This approach relies on the introduction of multipole expansions (similar to Taylor series expansions) to replace the sources on SO by equivalent point sources placed on rectangular 86

grids. In this manner the CPU requirements can be reduced by making use of three dimensional FFT routines. It can again be shown that AIM reduces the boundary matrix CPU requirements down to O (N,15) or even down to O (N133). 3.2 Finite Element-Absorbing Boundary Condition Method 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 [55], and Engquist and Majda [56]. More recently, other ABCs such as those by Webb and Kanellopoulos [57], and Chatterjee and Volakis [4] (see also Senior and Volakis [2]) have been proposed. All ABCs provide an approximate relation between the E and H fields 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 second order ABC derived by Chatterjee and Volakis [4] is given by jkoZo(Hscat x ) = a ' E-tt + 3 a [ (V(V x Escat),] + V,(V E at), (17) where a = atlttl + a2t2t2 = 77(t^ + t2t2) ltl + (jko + 3K, -2K1 )t t + (jko + 36Km 9 2K2)t2t2 jko Km jko Km a1 = 1[4K - Kg + D(jko - K) + + mAK] a2 7= [4K - Kg + D(jko - K2) + K2 - KmAK] 1 71~ D-2}m D = 2Km AK = K - K2 D = 2jko +5m -- Klm In the above equation, the subscripts n and t denote surface normal and tangential components of the field, tl,t2 are the principal surface directions and K1i,^2 are the corresponding curvatures, Km = (Ki + K2)/2 is the mean curvature and Kn = K1K2 is the Gaussian curvature. For K1 = K2 (spherical termination boundary), (17) reduces to the ABC derived by Webb and Kanellopoulos [57]. The latter authors have recently proposed a correction [58] for the implementation of these ABCs which should be incorporated in existing codes. Upon substituting (17) into the surface integral of (6) and making use of vector differential and integral identities [59], we have jkozZo ESCt (H x n) dS = J['(Ecat)2 + E t)2]dS + Js 7(V x Escat)2dS JJo 87

Metal Patch r r= 1j2.7 Figure 6: Illustration of an antenna terminated by an isotropic absorbing material - f(V Ecat)[V ( at )t]dS. (18) Next, expansion (6) or (14) is introduced and (18) is differentiated as in (7) to obtain the finite element equations (8). Since the ABC is a local condition, the sparsity of the finite element matrix is preserved and the resulting system is again given by (12) with [G] removed. However, in the case of arbitrarily curved boundaries, the matrix will not be symmetric except for planar and spherical boundaries. The system will be also symmetric for cylindrical boundaries only if linear edge-based expansion functions are used. These ABCs have been extensively validated for a number of antenna and scattering configurations[59], [60]. Higher order ABCs have been proposed for a more effective suppression of the outgoing waves. These are often difficult to implement but as demonstrated by Senior et. al. [61], they provide greater accuracy and effort being placed so much closer to the scatterer or radiator. 3.3 Finite Element-Isotropic Absorber Method In accordance with this method, the outgoing waves are suppressed by an absorbing dielectric layer placed at some distance from the antenna as shown in Figure 6. Typically, the layer is backed by metal and the finite element region extends all the way to the metal at which point the mesh is terminated by setting the tangential component of the total electric field to zero. This is equivalent to removing the integral over SO in (2)-(6). Obviously, the accuracy of the technique depends on how well the metal-backed absorber suppresses the incoming waves without introducing backward reflections. The plane wave reflection coefficient of the absorber shown in Figure 6 has been minimized over the entire visible angular range [46]. The fact that the dielectric has the same relative permittivity and permeability ensures that there is no reflection from the air-dielectric interface at normal incidence (perfect impedance match), but the performance degrades away from normal with the reflection coefficient reaching unity at grazing. This limitation led researchers to consider perfectly matched interfaces as discussed next. Nevertheless, even though the isotropic absorbers are not perfect, they are still useful in modeling antennas on doubly curved platforms. Figure 7(a) displays a sectoral microstrip patch printed on a conical surface and Figure 7(b) shows the antenna 88

Metallic Ground h h = 0.114cm Plane i f. 200e ' ' ' 0=350 2= 18.32 /_ _ __... --- 100/ Figure 7: 12.9Sectoral microstrip patch on cone quency drops with the cone angle for the same patch dimensions. It should be also remarked that the computed resonance frequency is within 3.2 percent of that predicted by the approximate cavity model, and this is reasonable. In generating the results given in Figure 7, the computational domain was discretized using 2, 358 prisms resulting in 3, 790 degrees of freedom. One frequency run took 5.5 minutes on an HP9000 (Model 715/64) workstation rated at 22 MFlops. 3.4 Finite Element Anisotropic Absorber Method This mesh termination method is similar to that described in the previous section except that the layer is comprised of an anisotropic lossy dielectric which has zero reflection coefficient at the air-dielectric interface over all incident angles and can therefore be considered as a perfectly matched interface (PMI) [20]. The intrinsic parameters of the PMI layer are = - or - 0 j a 0 0 1 40 e-jj3 where a = (Z/Zo)2, Z is the intrinsic impedance of the medium being terminated by the absorber, whereas oa, f and the thickness of the layer are parameters which can be optimized for maximum absorption. The advantage of the anisotropic (over the isotropic) artificial absorber in terminating finite element meshes is illustrated in Figure 8 [62] where it is shown that the anisotropic absorber retains its low reflectivity at oblique incidences except near grazing (q = 0~). The performance shown in Figure 8 is based on a purely theoretical analysis but as shown in Figure 9, this performance can be realized in numerical simulations with careful choices of /3t and the sampling rate N within the absorber. It has been found [63] that 89

A 0 m 5' -10 -20 -30 -40 -50 -60 -70 Isotropic Er= r= b -80 ' I I I I I I 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Anisotropic ~r= r= b I - (b - /b) x x sin ~ (a) (b) Figure 8: Reflectivity of a A/4 thick planar metal-backed perfectly matched absorber (o = 1 - j2) as a function of incidence angle. (a) Geometry, (b) plane wave reflection coefficient vs. angle 12 cm -f 1.05 cm..i 0.21 cm -- 3.095 cm Uniaxial absorber Microstrip -" 0.75 1.00 1.25 1.50 1.75 2.00 ( (a) (b) Figure 9: Reflection coefficient of the PMI for terminating a microstrip line as extracted from a numerical implementation 90

the reflectivity curve shown in Figure 9 is typical to most situations and the location of the minimum reflection coefficient can be predicted a priori. An extensive numerical study based on a two dimensional planar absorber model demonstrated that ft/Atx (Ax = Ao/ cos q) is a scalable parameter and that given the desired reflection coefficient IRI, the formulas [62] Mt - -0.01061RI + 0.0433 (19) Ax N = 0.147exp[7.353 y] (20) Ax can be used to choose ft/Ax and the minimum number of discrete samples to achieve the desired absorption. In these, \R\ is specified in dB and AO denotes the free space wavelength. These formulas are in agreement with the actual numerical results in Figure 9 but it has yet to be determined how well they apply for curved perfectly matched layers which are placed conformal to scattering and radiating surfaces. Improvements to their absorptivity though can be attained by considering tapered layers and formulas similar to (19)-(20) are given by Legault [62] for one such tapered anisotropic absorber. 4 Feed Modeling For scattering problems where the plane wave incidence is usually the 'source', the righthand-side excitation has been explicitly given in (12) and will not be discussed further. However, for antenna parameter computations, the explicit form of {IKe} in (11) will depend on the type of feeding scheme being employed. Below we discuss specific forms of {fKe} corresponding to different feeding choices. Simple Probe Feed: For thin substrates the coaxial cable feed may be simplified as a thin current filament of length I carrying an electric current I 1. Since this filament is located inside the cavity, the first term of the integral in (2) or (3) needs to be considered for this model. Specifically, the ith (global) entry of the excitation vector Ki becomes Ki = jkoZoI I. W(r), = j3 2...,jm (21) where r is the location of the filament, m is the number of (non-metallic) element edges and jm is the global edge numbering index. In general, m such entries are associated with m element edges, and thus the index i goes from ji up to jm. This expression can be further reduced to Ki = jkoZoI 1, provided that the ith edge is coincident with the current filament. Voltage Gap Feed: This excitation is also referred to as a gap generator and amounts to specifying a priori the electric voltage V across the opening of the coax cable or any other gap. Since VI = E d, where d is a vector whose magnitude is the gap width, and V E the electric field across the gap, we have that Ei - d, where cos0i is equal to 1 dcosOi if the ith edge is parallel to d. Numerically, this gap voltage model can be realized by first setting the diagonal term Aii of the [A] matrix equal to unity and the off-diagonal 91

terms A (i - j) to zero. For the right-hand-side vector, only the entry corresponding to the ith (global) edge across the gap is specified and set equal to the value Ei whereas all other entries associated with edges not in the gap are set to zero. Coaxial Cable Feed Model: The simple probe feed of the coaxial cable is accurate only if the substrate is very thin. For thicker substrates, an improved feed model is necessary and this can be achieved by evaluating the functional Fc= -jkoZo s (E x H). dS (22) f over the aperture Sf of the coax cable. Assuming a TEM mode across Sf, the fields within the cable may be expressed as (see Figure 10) [64] E = -r, H =, (23) r r with ho = - eo+-. (24) ZO 7r In these expressions, Crc is the relative permittivity inside the cable, E and H are cavity patch T -____ ------- ------------- 2b cavity-cable junction (a) (b) Figure 10: (a) Side view of a cavity-backed antenna with a coax cable feed; (b) Illustration of the FEM mesh at the cavity-cable junction (the field is set to zero at the center conductor surface). the electric and magnetic fields, respectively, measured at z = 0 and Io is the center conductor current. Also, (r, 4, z) are the polar coordinates of a point in the cable with the center at r = 0. We observe that (24) is the desired constraint at the cable junction in terms of the new quantities ho and eo which can be used as new unknowns in place of the fields E and H. However, before introducing Fc into the system, it is necessary to relate eo and ho to the constant edge fields associated with the elements in the cavity region which border the cable aperture. Since the actual field has a 1/r behavior in the cable, we find that A/V = E(b - a) = eoln- i = Np(p = 1, 2,... Nc) (25) a 92

where AV denotes the potential difference between the inner and outer surface of the cable and Np denotes the global number for the edge across the coax cable. When this condition is used in the functional Fc, it introduces the excitation into the finite element system without a need to extend the mesh inside the cable or to employ a OFc fictitious current probe. The derivation of and its incorporation into the system is dE, then a straightforward task [64]. As can be expected, the above feed model assumes the presence of only the dominant(TEM) mode at the cavity-cable junction, an assumption which may not be suitable for certain applications. Of course, the model can be improved by extending the mesh (say, a distance d) into the cable. The equi-potential condition will then be applied at z=-d, where all higher order modes vanish. Aperture-Coupled Microstrip Line Model: As shown in Figure 11, when the microstrip antenna is fed with a microstrip line network underneath the ground plane (cavity's base) via a coupling aperture, special treatment of the feed structure must be considered in the FEM formulation. This is because the microstrip line is usually designed to have different size and shape as compared to the cavity's geometries. Hence, the conventional simulation of treating the entire 3-D domain using a single type of element is not efficient or appropriate for this feed. Referring to Figure 11, it is appropriate to separate the computational domains because of the small element size required in modeling the guided feed structure. As an example, let us consider a rectangular aperture which has been extensively employed in practice. The cavity fields may be discretized using tetrahedral elements, whereas in the microstrip line region rectangular bricks are the best candidate since the feed structure is rectangular in shape and the substrate is of constant thickness. Although both types Antenna Elments - -(I) St (11) Truncation Plane Coupling Aperture Sa Figure 11: Cross-section of an aperture coupled patch antenna, showing the cavity region I and the microstrip line region II for two different FEM computation domains. of elements employ edge-based field expansions, the meshes across the common area (coupling aperture) are different, and consequently some connectivity matrix must be introduced to relate the mesh edges across the aperture. However, since the aperture is very narrow, a 'static' field distribution may be assumed at any given frequency. Therefore, the potential concept may be again applied to relate the fields below and above 93

the aperture. Specifically, the 'equi-potential' continuity condition is enforced, and to demonstrate its implementation, let us first classify the slot edges as follows Tetrahedral Mesh (Cavity Mesh): * E 1 j = 1,2,3,... vertical edges * Ej2 j =1,2,3,... diagonal edges Brick Mesh (Feed Mesh): * E2 j =1,2,3,... vertical edges only Then the 'equi-potential' continuity condition requires that Ebl = eE t Eb2 = (ejE + ei+iEJ+i) 2d 3 in which — 1' { +1 i = 1 ' whereas t and d are the lengths of the vertical and diagonal edges, respectively. That is, t is simply the width of the narrow rectangular aperture. The coefficient ej is equal to ~1 depending on the sign conventions associated with the meshes above and below the coupling aperture. This connectivity scheme can be extended to entirely different computational domains. Also, it is apparent that the approach makes the FEM implementation straightforward for different geometry/size domains that would be significantly inefficient if only one type of elements were used for modeling the structure. In addition, the technique ensures a good system condition since the number of distorted elements in the mesh is minimized. 5 Parallelization When considering 3D problems of practical interest, the unknown count of the computational domain can easily reach several million degrees of freedom. The sparsity of the FEM system (particularly for the FE-ABC and FE-AA methods) makes possible the storage of such large scale problems but even at O(N) computational demands, their practical solution requires efficient use of parallel and vector platforms. Modern computing platforms can now deliver sustained speeds of several GFlops and CPTJ speeds in the Tflops range are within our reach. The inherent sparse matrices of PDE methods are particularly suited for execution on multiprocessor and vector platforms but the exploitation of these processors requires special storage schemes and techniques to perform the matrix-vector product required in the iterative algorithms at the Flop rates sustained on these multiprocessors. 94

Complex Operation + Matrix-vector Products nze nze-N Vector Updates 4N 3N Dot Products 3N 3N Total # of Operations 29N 3N N = # of unknowns nze = # of nonzero matrix elements Table 2: Floating Point Operations of BCG Per Iteration To parallelize and vectorize the FEM codes, it is essential to first optimize the execution of the iterative solvers which typically take-up 90% of the CPU time. Among them, the conjugate gradient algorithms (CG, BCG, CGS and QMR) have been found very attractive and a brief comparison of the pros and cons for these is given in [45]. The Generalized Minimal Residual Method (GMRES) is another iterative solver which can exhibit faster convergence rates. However, it stores the direction vectors and as a result it requires much higher storage. For the discussion below we will primarily concentrate on the BCG and QMR algorithms and we note that the symmetric form of BCG requires minimal number of arithmetic operations (see Table 2). A disadvantage of the BCG is its erratic convergence pattern whereas the QMR has smooth and monotonic convergence. However, neither BCG nor QMR can guarantee convergence and typically they both converge or not for the same problem. When considering the parallelization of a completely sparse system such as that resulting from the FE-ABC method, the following issues must be addressed: Storage of Sparse System: The performance of the code is strongly dependent on the employed storage scheme. Since a typical FEM matrix has about 8.5 Nee or so non-zero entries, it is essential that the non-zero elements be stored in a manner that keeps the storage requirements nearly equal to the non-zero entries and minimizes inter-processor communications. The ITPACK [65] and the compressed row storage (CRS) schemes are most appropriate for parallel computing. The ITPACK storage format is most efficient for retrieving the matrix elements and certainly the choice method when the number of non- zero elements are nearly equal for every matrix row. Basically, the ITPACK format casts the FEM matrix in a smaller rectangular matrix having the same rows as the original matrix and can be unpacked by referring to a pointer integer matrix of the same size. However, this rectangular matrix can contain as much as 50% zeros which results in space wastage. By using a modified ITPACK scheme, space wastage can be reduced down to 30%. Even with less wastage, the CRS format may be the most efficient storage scheme with some compromise in CPU speed. It amounts to storing [A] as a single long row which can be uncompressed using two integer pointer arrays. For the symmetric BCG algorithm, the CRS format results in only 8.5 N complex numbers and N integers. However, it should be pointed out that the CRS format is not appropriate for vector processors such as the C-90. For vectorization, it is best to organize the storage in 95

Processors # of Processors, P Tiem (u-secs/iteration/unknown) Cray C-90 1 (275 MFlops) 0.55 KSR 1 28 1.28 58 0.57 8 3.42 Intel Paragon 16 1.99 32 1.38 IBM SP-1 4 1.47 Table 3: CPU Time Per Unknown for Solving Typical FE-ABC Systems sections of long vectors and to achieve this for our type of matrices the jagged diagonal format [66] appears to work best. Using this format the rows are reordered so that the row with the maximum number of non- zeros is placed at the top of the matrix and rows with the least non-zero entries are shuffled to the bottom. This reordering greatly enhances vectorization because it allows tiling of the shorter rows to yield very long vector lengths in the matrix-vector multiplication phase. Specifically, for some problem the jagged diagonal storage format allowed the matrix-vector multiplication routine to run at about 275 MFlops on a Cray C-90 whereas the same routine ran at 60 MFlops using the CRS format. The dot product speeds and the vector updates reached 550 MFlops and 600 MFlops for the same problem. Table 3 provides a relative comparison of CPU estimates on various computers. Interprocessor Communications: For distributed memory platforms, the method of partitioning the stiffness matrix [A] among the processors, the chosen storage scheme and the inherent unstructured sparsity of [A] are all crucial to the overall speed of the code. An approach that has worked well on massively parallel processors (such as the SP-2, Intel Paragon, Convex Exemplar) is that of assigning each processor a section of the matrix and by dividing the vectors among the P processors. Thus, each processor is responsible for carrying out the matrix-vector product for the block of the matrix it owns. However, the iterate vector is subdivided among all processors, and therefore narrow-band or structured sparse matrices have an advantage because they reduce interprocessor communication. Since typical FEM matrices are unstructured, algorithms such as the Recursive Spectral Bisection (RSB) have been found very effective in reducing inter-processor communication. However, the standard Gibbs-Pool-Stockmeyer profile reduction algorithm has been found even more effective in reducing the initial FE-ABC matrix (see Figure 3) to banded form as illustrated in Figure 12. This type of matrix reordering can deliver speed-ups as close to linear as possible. Matrix Preconditioning: Preconditioned iterative solvers are intended to improve the convergence rate of the algorithm. At times, preconditioners are necessary as may be the case with some dielectrically loaded structures. However, for relatively small systems (less than 100,000 unknowns) it has been found that diagonal preconditioning is typically most effective and should always be applied. This preconditioning amounts to normalizing each row by the largest element, but even this simple operation can lead to 96

0 0.5 1 1.5 2 2.5 3 nz=469151 4 X 1 0 Figure 12: Reduced bandwidth of the FE-ABC system after application of the GibbsPool-Stockmeyer profile reduction algorithm substantial convergence speed-ups. Block and incomplete LU preconditioners are more effective in improving the convergence of the solver but are more costly to implement and one must judge on the overall CPU requirements rather than on the improved convergence alone. For example, the incomplete LU preconditioner given in [671 reduced the iterations to 1/3 of those needed with diagonal preconditioning. However, each iteration was 3 times more expensive due to the triangular solver bottleneck. 6 Additional Applications We choose three more examples to demonstrate the capability of the hybrid finite element methods. Scattering by a Large Cone-Sphere: A cone-sphere is basically a hemisphere attached to a cone. This is a difficult geometry to mesh since a surface singularity exists at the tip of the cone. The singularity can be removed in two ways: i) by creating a small region near the tip and detaching it from the surface or ii) by chopping off a small part near the tip of the cone. The second option inevitably leads to small inaccuracies for backscatter from the conical tip; however, we chose this option since the conical angle in our tested geometry was extremely small (around 7~) and the mesh generator failed to mesh the first case on numerous occasions. In Figure 13, we plot the backscatter patterns of a 4.5A long cone-sphere having a radius of 0.5A for 00 polarization. The mesh truncation surface is a rectangular box placed 0.4A from the surface of the cone-sphere. As seen, the far-field results compare extremely well with computations from a body of revolution code [68]. Frequency Selective Surfaces (FSS): FSS structures [69] are arrays of tightly packed periodic elements which are typically sandwiched between dielectric layers. The periodic 97

10 10 --l':/ -l o b -20 --30- 4 FEMATS -40-,.... I..... I..... I.......... - -40 --90 -60 -30 0 30 60 9( Observation Angle 0o, deg. Figure 13: Backscatter pattern of a perfectly conducting conesphere for q(qX and 00 polarizations. Black dots indicate computed values using the FE-ABC code (referred to as FEMATS) and the solid line represents data from a body of revolution code [68]. Mesh termination surface is a rectangular box. 98

elements may be of printed form or slot configurations designed to resonate at specific frequencies. As such, they are penetrable around the element resonances and become completely reflecting at other frequencies. To meet bandwidth design specifications, stacked element arrays may be used in conjunction with dielectric layer loading. Here we consider the analysis of FSS structures via the FE-BI method. Because of the fine geometrical detail associated with the FSS surface, the finite element method has yet to be applied for the characterization of FSS structures, but use of prismatic elements makes this a much easier task. Of particular interest in FSS design is the determination of the transmission coefficient as a function of frequency, and since the array is periodic, it suffices to consider a single cell of the FSS. For computing the transmission coefficient T, the periodic cell is placed in a cavity as shown in Figure 14 and the structure is excited by a plane wave impinging at normal incidence. Assuming that near resonance the wave transmitted through the FSS screen will retain its TEM character, the transmission coefficient of the FSS panel can be approximated as T) = 10 log dB O where a is the reflection coefficient of the absorber placed at the bottom of the cavity and should be kept small (< 0.1) to suppress higher order interactions. By adding the next higher order interaction, a more accurate expression for the transmission coefficient is TdB r TJ() + 10 log [i - a(i - T())] The above FSS modeling approach was applied for a characterization of multi-layered slot FSS structures. The geometry of the multilayer radome is given in Fig. 14. The total thickness of the FSS was 6.3072cm and is comprised of two slot arrays (of the same geometry) sandwiched within the dielectric layers. For modeling purpose, a 1.54cm thick absorber is placed below the FSS as shown in Fig 14. It is seen that the results generated by the FE-BI method are in good agreement with the measurements [70]. Scattering by Jet Engine Inlets: This example demonstrates a hybridization of the finite element and high frequency methods for scattering by jet engine structures[71]. In this hybridization, the jet engine itself is modeled by the FEM taking advantage of the geometrical adaptability of finite methods. The fields generated by the FEM are then propagated to the mouth of the jet engine using a ray technique such as the geometrical theory of diffraction or the shooting and bouncing ray methods. To analyze the scattering by the jet engine, the blade structure is placed in a computational domain which is truncated by a perfectly matched absorbers at about 0.5 wavelengths from the engine face. By exciting the engine on a mode by mode basis, the appropriate modal scattering matrix is generated using the procedure described in [71]. However, in this case the computational domain is confined to a single blade lobe (cylindrical periodic cell) using appropriate phase boundary conditions at the boundaries of the periodic cells. This reduces the computational requirements to very low orders and permits the simulation of practical jet engine geometries. Given the modal scattering matrix of the jet engine, high 99

slot array placed 60mils below top surface of 90mils layer slot array at 30mils below top of 9Orils layer r= 4 \. / 3 -90mils 9Omils || 5.85cm 2cm j absorber, ~ r 4.- j 10.4 =l -= _ 6.45cm I,... -- i -5 -10 0. -15 E c l — i 1 / measured - calculated -20 -25 - -30 0 0.2 0.4 0.6 0.8 1 1.2 frequency (GHz) 1.4 1.6 1.8 l Figure 14: Upper figure: geometry of the multilayer frequency selective surface (FSS) used for modeling; lower figure: measured and calculated transmission coefficient through the FSS structure 100

frequency techniques can be readily used for propagating the fields to the inlet opening as described by Karty et.al.[72]. A demonstration on the accuracy of this hybrid technique is given in Figure 15 where we show the computed scattered field due to a single mode excitation. The specific engine-blade geometry was placed in a cylindrical metal-backed cavity as shown in Figure 15. It consisted of 20 straight blades mounted radiallly on a cylindrical hub 9.25" in radius with the blades extending all the way to the outer wall of the cavity which was 14.375" in radius. The blades were 2 degrees thick and 2.75"long, whereastheentireinletdepthwasl9.685" and we remark that this configuration is similar to an F4 fighter engine. The shown computation was carried out at 4.348GHz and involved the solution of an FEM system with only 87,000 edges (64,000 elements) since a single periodic lobe was descretized but the overall structure spans about 4.5 wavelengths in diameter. As seen the agreement between the FEM computations and those of the reference modal solution are in excellent agreement, nearly overlaying each other. 7 Conclusion In this paper we reviewed hybrid finite element methods as applied to electromagnetic scattering and radiation problems. Much of the emphasis dealt with the various mesh truncations schemes and we presented an up-to-date account of these schemes. The usual finite element-boundary integral method was presented and new developments for reducing the CPU equireements of this technique using the fast integral methods were discussed. Antenna feed modeling in the context of the finite element method had not been discussed before and for the first time we presented an overview of the modeling approaches for the most popular antenna feeds, including aperture coupled feeds. Parallelization will continue to play an increasingly greater role and a section was included discussing our experiences for better implementation of finite element codes on distributed and vector architectures. A number of examples illustrating the successful application of the finite element method were included throughout the paper and these were intended to demonstrate the method's geometrical adaptability and inherent capability to treat highly heterogeneous structures. As can be expected, issues relating to mesh truncation, mixing of elements [73], domain decomposition[74, 75], robustness, adaptive refinement [76], accuracy, error control, feed modeling and parallelization for large scale simulations will continue to dominate future research and developments relating to partial differential equation methods. An apparent advantage of the finite element method is its potential hybridization with all other frequency domain methods. Future applications of the finite element method are likely to make greater use of hybridization techniques aimed at increasing the method's accuracy and efficiency while retaining its inherent geometrical adaptability and ease in handling materials. References [1] J. L. Volakis, A. Chatterjee, and L. C. Kempel, 'A review of the finite element; method 101

Comparison of RCS for fan termination 35.00] --- —--------------- - RCS pi -pIUJ Hidu HutuliUig 3.00 - RCS(.ph1-ph1l FK 30.00..... CS(theta-thet) Mode atehinq............... RCS(theta-theta) FEf4.0,............................................................................................ 25.00 20.00........ 15.00..............:..................... 00................................................... -5.00.. ~,. *' ' -20.00,,,!!,,..'.: -50.00 -36J00 -12.00 12.00 36.00 0.00 Theta (dog.) L..................................................................................................... Figure 15: Scattering by a cylindrical guide terminated with a straight bladed engine configuration. 102

for three dimensional scattering', J. Opt. Soc. of America, A, pp. 1422-1433, 1994. [2] T. B. A. Senior and J. L. Volakis, Approximate Boundary Conditions in Electromagnetics, London, IEE Press, 1995. [3] 0. C. Zienkiewicz, The finite element method, McGraw Hill, New York, 3rd edition, 1979. [4] 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., vol. 41, pp. 221-226, 1993. [5] J. P. Webb, "Edge elements and what they can do for you," IEEE Trans. Mag., vol. 29, pp. 1460-5, 1993. [6] D. Sun, J.Manges, X. Yuan, and Z. Cendes, "Spurious modes in Finite-element mothods," IEEE Antennas Propagat. Magaz., Vol.37, No.5, pp.12-24, Oct. 1995 [7] H. Whitney, Geometric Integration Theory, Princeton University Press, 1957. [8] J. C. Nedelec, "Mixed finite elements in R3", Numer. Math., vol. 35, pp. 315-41, 1980. [9] A. Bossavit and J. C. Verite, "A mixed FEM-BIEM method to solve 3D eddy current problems," IEEE Trans. Mag., vol. 18, pp. 431-5, March 1982. [10] A. Bossavit, "Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism," IEE Proc., vol. 135, pt. A, no. 8, Nov. 1988. [11] M. L. Barton and Z. J. Cendes, "New vector finite elements for three-dimensional magnetic field computation," J. Appl. Phys., vol. 61, no. 8, pp. 3919-21, Apr 1987. [12] J. S. van Welij, "Calculation of eddy currents in terms of H on hexahedra," IEEE Trans. Mag., vol. 21, pp. 2239-41, Nov. 1985. [13] J. F. Lee, D. K. Sun, and Z. J. Cendes, "Full-wave analysis of dielectric waveguides using tangential vector finite elements," IEEE Trans. Microwave Theory Tech., vol. 39, no. 8, pp. 1262-71, Aug 1991. [14] J. M. Jin and J. L. Volakis, "Electromagnetic scattering by and transmission through a three-dimensional slot in a thick conducting plane," IEEE Trans. Antennas Propagat., vol. 39, pp. 543-550, Apr 1991. [15] R. D. Graglia, A. F. Peterson and D. R. Wilton, "Higher order conforming vector bases on curvilinear elements," This issue. [16] J. F. Lee, D. K. Sun, and Z. J. Cendes, "Tangential vector finite elements for electromagnetic field computation," IEEE Trans. Mag., vol. 27, no. 5, pp. 4032-5, Sep. 1991. [17] G. Mur and A. T. deHoop, "A finite element method for computing three-dimensional electromagnetic fields in inhomogeneous media," IEEE Trans. Mag., vol. 21, pp. 2188-91, Nov. 1985. [18] J. S. Wang and N. Ida, "Curvilinear and higher order edge finite elements in electromagnetic field computation," IEEE Trans. Mag., vol. 29, no. 2, pp. 1491-4, March 1993. 103

[19] J. P. Webb and B. Forghani, "Hierarchal scalar and vector tetrahedra," IEEE Trans. Mag., vol. 29, no. 2, pp. 1495-8, March 1993. [20] Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition," IEEE Trans. Antennas Propagat., vol. 43, no. 12, Dec. 1995. [21] P. Silvester and M.S. Hsieh, "Finite element solution of 2-dimensional exterior field problems", Proc. IEE, vol. 118, pp. 1743-1747, Dec. 1971. [22] B.H. McDonald and A. Wexler, "Finite element solution of unbounded field problems," IEEE Trans. Microwave Theory Tech., vol. 20, pp. 841-847, Dec. 1972. [23] J. M. Jin and J. L. Volakis, "TM scattering by an inhomogeneously filled aperture in a thick ground plane", IEE Proc., Part H, vol. 137, no. 3, pp. 153-159, June 1990. [24] X. Yuan and D.R. Lynch and J.W. Strohbehn, "Coupling of finite element and moment methods for electromagnetic scattering from inhomogeneous objects," IEEE Trans. Anetnnas Propagat., vol. 38, pp. 386-393, March 1990. [25] X. Yuan, "Three dimensional electromagnetic scattering from inhomogeneous objects by the hybrid moment and finite element method," IEEE Trans. Antennas Propagat., vol. 38, pp. 1053-1058, 1990. [26] J. Angelini and C. Soize and P. Soudais, "Hybrid numerical method for harmonic 3D Maxwell equations: Scattering by a mixed conducting and inhomogeneous anisotropic dielectric medium," IEEE Trans. Antennas Propagat., vol. 41, no. 1, pp. 66-76, January 1993. [27] G.E Antilla and N.G. Alexopoulos, "Scattering from complex three-dimensional geometries by a curvilinear hybrid finite element-integral equation approach," J. Opt. Soc. Am. A, vol. 11, no. 4, pp. 1445-1457, April 1994. [28] P. Soudais, "Computation of the electromagnetic scattering from complex 3D objects by a hybrid FEM/BEM method," J. Elect. Waves Appl., vol. 9, no. 7-8, pp. 871-886, 1995. [29] K. D. Paulsen, X. Jia and J. Sullivan, "Finite element computations of specific absorption rates in anotomically conforming full-body models for hyperthermia treatment analysis," IEEE Trans. Biomedical Engr., vol.40, no.9, pp. 933-45, Sept. 1993. [30] T. Eibert and V. Hansen, "Calculation of unbounded field problems in free-space by a 3D FEM/BEM-hybrid approach,"J. Elect. Waves Appl., vol 10, no. 1, pp.61-78, 1996. [31] S. M. Rao, D. R. Wilton and A. W. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," IEEE Trans. Antennas Propagat. vol. 30, no. 3, pp. 409-418, May 1982. [32] T. Cwik, "Coupling finite element and integral equation solutions using decoupled boundary meshes (electromagnetic scattering), IEEE Trans. Antennas Propagat., vol. 40, no. 12, pp. 1496-1504, Dec. 1992. [33] J.M. Jin and 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 and Propagat. Society magazine, vol. 33, no. 3, pp. 22-32, June 1991. 104

[34] E. Arvas, A. Rahhal-Arabi, A. Sadigh and S. M. Rao, "Scattering from multiple conducting and dielectric bodies of arbitrary shape," IEEE Trans. Antennas Propagat. Soc. Mag., vol. 33, no. 2, pp. 29-36, April 1991. [35] J.D. Collins and J.M. Jin and J.L. Volakis, "Eliminating interior resonances in FE-BI methods for scattering," IEEE Trans. Antennas Propagat., vol. 40, pp. 1583-1585", December 1992. [36] J.M. Jin and J.L. Volakis, "A hybrid finite element method for scattering and radiation by microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas Propagat., vol. 39, pp. 1598-1604, 1991. [37] J.L. Volakis and J. Gong and A. Alexanian, "A finite element boundary integral method for antenna RCS analysis," Electromagnetics, vol. 14, no. 1, pp. 63-85, 1994. [38] J.M. Jin and J.L. Volakis, "Scattering and radiation analysis of three-dimensional cavity arrays via a hybrid finite element method," IEEE Trans. Antennas Propagat., vol. 41, pp. 1580-1586, Nov 1993. [39] J. Gong, J. L. Volakis, A. Woo, and H. Wang, "A hybrid finite element boundary integral method for analysis of cavity-backed antennas of arbitrary shape", IEEE Trans. Antennas Propagat., vol. 42, pp. 1233-1242, 1994. [40] L. C. Kempel, J. L. Volakis and R. Sliva, "Radiation by cavity-backed antennas on a circular cylinder," IEE Proceedings, Part H. pp. 233-239, 1995. [41] Y. Zhuang, K-L. Wu, C. Wu and J. Litva, "Acombined full-wave CG-FFT method for rigorous analysis of large microstrip antenna arrays," IEEE Trans. Antennas Propagat., vol. 44, pp. 102-109, January 1996. [42] J.D. Collins and 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. [43] K. Barkeshli and J.L. Volakis, "On the implementation and accuracy of the conjugate gradient FFT method," IEEE Trans. Antennas Propagat., vol. 32, pp. 20-26, 1990. [44] J.M. Jin and J.L. Volakis, "Biconjugate gradient FFT solution for scattering by planar plates," Electromagnetics, vol. 12, pp. 105-119, 1992. [45] J. L. Volakis, "Iterative Solvers," IEEE Antenna Propagat. Soc. Mag., vol. 37, no. 6, pp. 94-96, December 1995. [46] T. Ozdemir and J. L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for terminating finite element meshes', Radio Sci., vol 29, no. 5, pp. 1255-1263, Sept.-Oct. 1994 [47] C. Farhat and F-X. Roux, "An unconventional domain decomposition method for an efficient parallel solution of large-scale finite element systems," SIAM J. Sci. Stat. Comput., vol. 13, pp. 379-396, January 1992. [48] V. Rokhlin, "Rapid solution of integral equations of scattering theory in two dimensions," Journal of Computational Physics, vol. 86, no. 2, pp. 414-439, 1990. 105

[49] W. C. Chew, C. C. Lu, E. Michielssen and J. M. Song, "Fast solution methods in electromagnetics," This issue. [50] R. Coifman, V. Rokhlin and S. Wandzura, "The fast multipole method for the wave equation: A pedestrian prescription," IEEE Antennas and Propagat. Magazine, June 1993. [51] C.C. Lu and W.C. Chew, "Fast far field approximation for calculating the RCS of large objects," Micro. Opt. Tech. Lett., vol.8, no.5, pp.238-241, April 1995. [52] S. Bindiganavale and J. L. Volakis, "A hybrid FEM-FMM technique for electromagnetic scattering," Proceedings of the 12th annual review of progress in applied computational electromagnetics (ACES), Naval Postgraduate School, Monterey CA, pp 563-570, March 1996. [53] S.S. Bindiganavale and J.L. Volakis, "Guidelines for using the fast multipole method to calculate the RCS of large objects," Micro. Opt. Tech. Lett., vol.11, no.4, March 1996. [54] E. Bleszynski, M. Bleszynski and T. Jaroszerwicz, "A fast integral equation solver for Electromagnetic scattering problems," IEEE Antennas Propagat. Symposium Processdings, Seattle, WA, pp. 417-420, 1994 [55] A. Bayliss and E. Turkel, 'Radiation boundary conditions for wave-like equations', Comm. Pure Appl. Math., vol. 33, pp. 707-725, 1980. [56] B. Engquist and A. Majda, 'Absorbing boundary conditions for the numerical simulation of waves', Math. Comp., vol. 31, pp. 629-651, 1977. [57] 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, October 1989. [58] V. N. Kanellopoulos and J. P. Webb, "The importance of the surface divergence term in the'finite element-vector absorbing boundary condition method," IEEE Trans. Microw. Theory Tech., vol. 43, no. 9, pp. 2168-2170, Sept. 1995. [59] A. Chatterjee and J. L. Volakis, "Conformal absorbing boundary conditions for 3D problems: Derivation and applications', IEEE Trans. Antennas Propagat., vol. 43, no. 8, pp. 860-866, August 1995. [60] L. C. Kempel and J. L. Volakis, "Scattering by cavity-backed antennas on circular cylinder," IEEE Trans. Antennas Propagat., vol. 42, pp. 1268-1279, 1994. [61] T. B. A. Senior, J. L. Volakis and S. R. Legault, "Higher order impedance boundary conditions," IEEE Trans. Antennas Propagat., to appear. [62] S. R. Legault, T. B. A. Senior and J. L. Volakis, "Design of planar absorbing layers for domain truncation in FEM applications," Electromagnetics, to appear. [63] J. Gong and J. L. Volakis, "Optimal selection of uniaxial artificial absorber layer for truncating finite element meshes," Electronics Letters, vol. 31, no. 18, pp. 1559-1561, August 1995. 106

[64] J. Gong and J. L. Volakis, "An efficient and accurate model of the coax cable feeding structure for FEM simulations," IEEE Trans. Antennas Propagat., vol. 43, no. 12, pp. 1474-1478, December 1995. [65] D. R. Kincaid and T. K. Oppe, "ITPACK on supercomputers," Int. J. on Num. Methods, Lecture Notes in Math., vol. 1005, pp. 151-161, 1982. [66] E. Anderson and Y. Saad, "Solving sparse triangular linear systems on parallel computers," Int. J. of High Speed Computing, vol. 1, pp. 73-95, 1989. [67] A. Chatterjee, J. L. Volakis and D. Windheiser, "Parallel computation of 3D electromagnetic scattering using finite elements," Int. J. Num. Modeling.: Electr. Net. Dev. and Fields, vol. 7, pp. 329-342, 1994. [68] J.M. Putnam and L.N. Medgyesi-Mitschang, "Combined field integral equation formulation for axially inhomogeneous bodies of revolution," McDonnell Douglas Research Labs, MDC QA003, December 1987 [69] E.L. Pelton and B.A. Munk, "Scattering from periodic arrays of crossed dipoles," IEEE Trans. Antennas Propagat., vol. AP-27, pp. 323-330. [70] H. Wang, Personal communication, China Lake, CA, 1995 [71] D. C. Ross, J.L. Volakis and H. T. Anastassiu, A Hybrid finite element-modal analysis of jet engine scattering, IEEE Trans. Antenna Propagat., AP-43, pp. 277-285, Feb. 1995 [72] J. Karty, J. Roedder and S. Alspach, CAVERN: A prediction code f or cavity electromagnetic analysis, IEEE Ant. Propagat. Soc. Magaz., Vol. 37, No.3, pp. 68-72, June 1995. [73] N.E. Boyse and A.A Seidl, "A hyrid finite element method for 3D scattering using nodal and edge elements," IEEE Trans. Antennas Propagat., vol. 42, pp. 1436-1442, Oct. 1994 [74] T. Cwik,"Parallel decomposition methods for the solution of electromagnetic scattering problems," Electromagnetics, Vol. 42 pp. 343-357, 1992 [75] P. Le Tallec, E. Sallel and M. Vidrascu, "Solving large scale structural problems on parallel computers using domain decomposition techniques," Chapter in Advances in Parallel and Vector Processing for Structural Mechanics, (edited by B. Topping and M. Papadrakakis), CIVIL-COMP Ltd. Edinburgh, Scotland, 1994 [76] N. Golias, A. Papagiannakis and T. Tsiboukis, "Efficient mode analysis with edge elements and 3D adaptive refinement," IEEE Trans. MTT, Vol. 42 pp. 99-107, Jan. 1994 107