030601-9-T ANNUAL PROGRESS REPORT National Aeronautics and Space Administration Langley Research Center Hampton, VA 23681-0001 August 1996 30601 -9-T = RL-2428

NASA Grant No. NAG-1-1478 Simulation of Spiral Slot Antennas on Composite Platforms FORWARD The project goals, plan and accomplishments up to this point are summarized in the viewgraphs following this page. Among the various accomplishments, the most important have been: 1. The development of the prismatic finite element code for doubly curved platforms and its validation with many different antenna configurations 2. The design and fabrication of a new slot spiral antennas suitable for automobile cellular, GPS and PCS communications 3. The investigation and development of various mesh truncation schemes., including the perfectly matched absorber and various fast integral equation methods. 4. The introduction of a frequency domain extrapolation technique (AWE) for predicting broadband responses using only a few samples of the response. This report contains several individual reports most of which have been submitted for publication to refereed journals. For a report on the frequency extrapolation technique, the reader is referred to the UM Radiation Laboratory report A total of 14 papers have been published or accepted for publication with the full or partial support of this grant. Several more papers are in preparation. Published Papers or Accepted for Publication since Aug. 1994 1. T. Ozdemir and J.L. Volakis, "Triangular prisms for edge-based vector finite element analysis," Submitted to IEEE AP Trans. 2. J. Gong and J.L. Volakis, "An efficient and accurate model of the coax cable feeding structure for FEM simulations," IEEE Trans. Antennas and Propagat., Vol. 43, pp. 1474-1478, Dec. 1995. 3. D.M. Kingsland, J. Gong, J.L. Volakis and J-F. Lee, "Performance of an anisotropic artificial absorber for truncating finite element meshes, " Vol. 44, No. 7, IEEE Trans. Antennas and Propagat.., pp. 975-982, July 1996 4. J. Gong and J.L. Volakis," Optimal Selection of a Uniaxial Artificial Absorber Layer for Truncating Finite Element Meshes" IEE Electronics Lett.., pp. 1559-1561, Vol. 31, No. 18, 31st August 1995. 5. J. Gong and J.L. Volakis, " FEM modeling of slot spirals using prismatic elements" IEE Electronics Lett.., pp. 1559-1561, Vol. 31, No. 18, 31st August 1995. 6. S. Legault, T.B.A. Senior and J.L. Volakis, "Design of Planar Absorbing Layers for Domain Truncation in FEM Applications," Vol. 16, No. 4 pp. 451-464, Electromagnetics, July-August 1996. 7. S. Bindiganavale and J.L. Volakis, "A hybrid FE-FMM technique for electromagnetic scattering," to appear in IEEE Trans. Antenna Propagat. 1

8. D.M. Kingsland, J. Gong, J.L. Volakis and J-F. Lee, "Performance of an anisotropic artificial absorber for truncating finite element meshes, " Vol. 44, No. 7, IEEE Trans. Antennas and Propagat.., pp. 975-982, July 1996 9. J.L. Volakis, J. Gong and T. Ozdemir, "Applications of the FEM to Conformal Antennas," to appear in book edited by Itoh, Silvester and Pelosi (Wiley, 1996). 10. J.L. Volakis, T.B.A. Senior, S.R. Legault, T. Ozdemir and M. Casciato," Artificial absorbers for truncating finite element meshes," to appear in Direct and Inverse Electromagnetic Scattering, eds. Serbest and Cloude, Pitman Research Series in Mathematics. 11. Volakis, J. Gong and T. Ozdemir, "Hybrid Finite Element Methods for Electromagnetics: Applications to Antennas and Scattering," to appear in book by ICASE 12. M.W.Numberger and J.L. Volakis, "A new planar feed for slot spiral antennas," IEEE Trans. Antennas and Propagat., Vol. 44, pp. 130-131, Jan. 1996. 13. J.Gong and J.L. Volakis, "An AWE Implementation oor Electromagnetic FEM Analysis," submitted for publication 14. J.Gong, J.L. Volakis and H. Wang "Efficient finite element simulation of slot antennas using prismatic elements," to appear in Radio Sci. LIST OF DOCUMENTS IN THIS REPORT A Planar Slot Spiral for Conformal Vehicle Applications M. Nurnberger, J. Volakis, D. T. Fralick and F.B. Beck This document provides a brief description of the design, fabrication and measurement of the full scale slot spiral antenna and infinite balun feed for this antenna. The antenna covers the Cellular and GPS bands and is intended for automobile deployment. The measurements at the NASA Langley facility demonstrate the CP performance as well as improvements obtained with the application of a resistive coating for terminating the arms. Triangular Prisms for Edge-Based Vector Finite Element Analysis of Conformal Antennas T. Ozdemir and J. Volakis This paper describes the theory of the antenna analysis code FEMA-PRISM. Several validations are given to demonstrate the code's suitability for doubly curved antennas. This is the first code for antenna analysis on doubly curved platforms and was jointly supported by the Rome Laboratory. Comparison of Three FMM Techniques for Solving Hybrid FE-BI Systems S. Bindiganavale and J. L. Volakis The paper provides a critical look at the fast multipole methods using speed and accuracy for benchmarking purposes. This work is towards our efforts to improve the accuracy of the finite element mesh truncation schemes while maintaining memory and CPU speed at tolerable levels. Artificial Absorbers for Truncating Finite Element Meshes 2

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 derived/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 Hybrid Finite Element Methods for Electromagnetics: Applications to Antennas and Scattering J. Volakis, J. Gong and T. Ozdemir A lengthy up-to-date review of the 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 scattering. An introduction to the AWE frequency extrapolation technique is also included for modeling microwave circuits using the finite element method. This document was submitted for inclusion to a book prepared by ICASE. 3

4

(.A

9

Conformal Antenna and Radome Geometries Simulated via the Finite Element Method Conformal Antenna/Array Conformal Artificial Termination Surface FEM Mesh Te rmin ated by Assumin g Dominant Mode \0.18128 \d-t =1.00 Wp, ~=1.12 cm ' 2rT7 0:::.::^. I" E 28 cm I 0.16764cm air 0 II 0.055 88 cm rCCU slot array placed 60mils below top surface of 90mils layer slot array at 30mils below top of 90mils layer ~. = 4. 0 mils Planar Ground Plane - \absorber, E r=4. -j10.4 6.45cm University of Michigan 1

Project Goals * Develop Robust Finite Element Codes for Narrow Slots and Other Conformal Antennas on Doubly Curved Surfaces (new elements, mesh truncation, solvers, feed models etc.) * Develop Miniaturation and Design Techniques for Conformal Small Scale UHF and VHF Antennas * Investigate Feeding Techniques for Slots Antennas...Integrate Antenna and Feed Networks * Develop Specialized Preprocessing Tools for FEM Antenna Analysis University of Michigan 2

Plan * Year 1 - Develop FEM code(s) for narrow slots and spirals - Examine new feed modeling techniques for slot antennas (modeling and fabrication methods) - Preprocessing tools - Mesh Truncation methods * Year 2 - Examine miniaturization approaches and measure performance - Analysis on doubly curved platforms - Examine curvature effects * Year 3 - Integrate high frequency methods with antenna FEM codes - Antenna development studies for multifrequency and multifunction capabilities University of Michigan 3

Progress: FEMA-PRISM Code New FEM code (FEMA-PRISM) was developed using prismatic meshes - Code(s) employs a new edgebased prism for ease in modeling - Uses special preprocessors for mesh generation without a need to employ sophisticated meshing packages.Only the surface mesh Dist need be specified Pris - Many applications and Validations of the code were carried out during the 2nd year * New feed modeling techniques were developed in the context of FEM. Coaxial and aperture coupled feeds were examined torted m University of Michigan 4

Progress: New Archimedean Slot Spiral Slot Radome Modeling c..viy slot array placed 60mils below top surface of 90mils layer Tw n-wUT mJnI pirui nil m B 9Oamils absorber, Er=4. j10.4 I --- * Archimedean spirals were modeled and measured from 800 to 1200MHz. * A new microstrip feed was eveloped for slot antennas * Absorbing boundary conditions were examined for FEM mesh truncation conformal to the printed antenna * Several printed and slot antennas as well as slot radomes were examined (placed on planar and curved treated platforms) University of A. lichigan 5 */

Progress: Perfectly Matched Layers for Mesh Truncation Uniaxial mArtificial Absorber * A most important development was the Cavity -- - 1 l ~I 21 introduction, design and performance evaluation of a perfectly matched(PML T-i' ^ Iabsorber for truncating FEM meshes. T] 0.21cm 3.53c PML is ease to implement and can yielc almost any desired level of absorption. developed design guidelines for this absorber X. E.... * Developed an interface to integrate FEM and HF codes for pattern evaluation of printed antennas mountec on conformal platforms. Compared with measured data collected at China Lake(NAWC) and at MRC(Dayton) University of Michiga ) his We n 6 I

Progress: Frequency Extrapolation Techniques * Developed and Made Successful use of the Asympotic Waveform Evaluation(AWE) method for predicting the frequency response of antennas and microwave circuits using only a few calculated samples of the response. This method and its generalized version, Complex Frequency Hopping(CFH), has only been recently introduced in modeling the response of VLSI circuits and its successful application to electromagnetics is an important development. It can eb an important companion to hybrid FEM methods and can be used to generate several time domain responses from a few frequency domain points. University of Michigan 7

Demonstration of AWE With FEM-PML for (M)MIC Application 1.06 cm 1 - 2.38 cm Frtii of t - 1.... i, f' 2 2.5 frequency (GHz) 44.h1. ai:K. 8th..~ or AWE,.3 t.e t () l- lt.1.8 G..z si ig. l.. p:.).oi.:.t

Progress Summary * In summary, a new class of software has been developed for the characterization of slot and patch antennas on conformal platforms * The new codes (FEMA-PRISM and FEMA_TETRA) are capable for modeling a large class of antennas (dual patch, slot or microstrip fed, spirals, log-periodics) in the presence of platform treatments and bandpass radomes. * FEMA-PRISM is self-contained and does not require external meshing preprocessor * A new archimedean slot spirals was developed and patented for conformal installation on a automobiles and operation at broad frequencies from 800MHz up to 1900Mhz. * Measurements of the spiral antenna were carried out at the NASALangley anechoic chamber facility and these will continue during the coming fiscal year. * The frequency extrapolation techniques AWE and CFH were introduced as a means to recover broadband responses from a few frequency data. University of Michigan * 9

Current Emphasis Current emphasis is on the - design and characterization of broadband antennas such as the new archimedean spiral, log-periodics, modulated spirals etc. - Integration of antennas and microwave circuit networks - examination of ferrite, high contrast and anisotropic material loadings for miniaturation - further development of the AWE and CFH techniques - hybridization of FEM antenna codes, high frequency and fast integral methods for a characterization antenna-platform interactions for predicting automobile antenna patterns - Improved hybrid FEM algorithms for a new class of antennas involving a mix of small and large scale details. University of Michigan 10

A PLANAR SLOT SPIRAL FOR CONFORMAL VEHICLE APPLICATIONS M.W. Nurnberger*, J.L. Volakis*, D.T. Fralickt, F.B. Beckt * Radiation Laboratory Dept. of Elect. Engin. & Comp. Sci. University of Michigan, Ann Arbor, MI 48109-2122 t NASA-Langley Research Center Hampton, VA 23681 ABSTRACT A slot spiral antenna and its associated feed are presented for conformal mounting on a variety of land, air, and sea vehicles. The inherent broadband behavior and good pattern coverage of the spiral antenna is exploited for the integration of multiple frequencies, and thus multiple transmitting and receiving apertures, into one compact, planar antenna. The feasibility of the broadband slot spiral antenna relies on the use of an equally broadband, balanced, planar, and non-intrusive feed structure. The design of the slot spiral, its feeding structure, and the reflecting cavity are discussed with emphasis on the experimental validation and construction of the antenna. Keywords: Spiral Antennas, Planar, Conformal, Automotive, Antenna Measurements 1.0 INTRODUCTION Spiral antennas are particularly known for their ability to produce very wide-band, almost perfect circularly-polarized radiation over their full coverage region. Because of this polarization diveristy and broad spatial and frequency coverage, many applications exist, ranging from military surveilance, ECM, and ECCM to numerous commercial and private uses, including the consolidation of multiple low gain communications antennas on moving vehicles. For the typical wire spiral antenna, the performance advantages mentioned above come at the price of size and complexity. While the radiating elements of a wire spiral may be planar, the feed network and balun structures generally are not, and combine to add weight, depth, and significant 17

compelexity to the system. Furthermore, because the spiral antenna radiates bi-directionally, an absorbing cavity is generally used to eliminate the radiation in one direction, adding even more depth to the antenna. While some designs exist [1] that integrate the feed and balun into the cavity and reduce the complexity somewhat, the absorbing cavity is still at least a quarter-wavelength deep at the lowest frequency of operation, adding significant thickness to the antenna [2]. A slot spiral is not burdened with most of these difficulties. As was demonstrated previously [3], the balun and feed structure may be integrated into the planar radiating structure. This greatly simplifies the construction and increases the accuracy of the antenna by allowing standard printed circuit techniques to be used throughout the entire facrications process. Also, while the slot spiral also radiates bi-directionally, a deep, absorbing cavity is not necessary for uni-directional radiation. In fact, for conformal mounting, a very shallow reflecting cavity is sufficient, and also serves to increase the gain somewhat. 2.0 SPIRAL ANTENNA OPERATION The operation of a spiral antenna may most easily be understood by considering a two-wire transmission line that has been wrapped around itself to form a spiral. Clearly, prior to its deformation, the two-wire line did not radiate, as the currents in the adjacent wire were out of phase, each cancelling the other's radiation in the far field. Following the rearrangement, at the center of the spiral, the currents on the two adjacent wires are still 1800 out of phase. However, as the currents travel down the two wires, because they have traveled different distances, they reach a region where, instead of being out of phase with their neighbors, they are in phase. This annular region is called the active or radiating region, because here, instead of cancelling the effects of their neighbors, each current reinforces the others, creating appreciable radiation in the far field. It is interesting to note that, for each current element in this active region, there is another current element that is in both space and phase quadrature. Thus, the radiation from all the current elements in the active region combines in the far field to produce circularly polarized radiation. In order to be considered frequency independent, an antenna must obey both the "angle principle" and the "truncation principle," as defined by Rumsey [4]. Unfortunately, from its geometry, and as can be concluded by studying Figure 3, the Archimedean spiral antenna follows neither. Figure 3 shows the radiation pattern of a spiral antenna, operated at a frequency inside its "frequency independent" frequency range, measured with a spinning linear source. Specifically, the high axial ratio indicates that there is a significant amount of opposite sense circularly polarized energy combining with the desired sense circularly polarized radiation. This behavior, which is present over the entire operational frequency range, has been determined to be a result of the truncation of the spiral arms [2,5], and indicates that the currents are not dying out before reaching the end of the antenna. Rather, they are traveling to the end of the antenna, and reflecting back along the spiral arms into the active region, from which they re-radiate energy with the opposite sense circular polarization. Thus, a large portion of the effort in designing an Archimedean spiral antenna lies in minimizing this reflection. In essence, the structure must be made "frequency independent." Another important aspect of spiral antenna design is the necessity of a "frequency independent," balanced feed. Several different techniques have been developed [1,2], most notably Dyson's "infinite balun" [6], a form of which is used in the antenna under development. Besides providing a bal 18

anced feed from the unbalanced coaxial line, a properly designed balun prevents the pattern squint so often observed in early spiral antenna designs. Most spiral antenna designs also include a backing cavity to force uni-directional operation. It has been determined that, for wire spirals, an absorbing cavity offers the best trade-off in terms of frequency bandwidth and physical size [2,5]. For slot spirals, however, a reflecting cavity can be used to great effect, as the available bandwidth is much larger than that obtained when coupled with a wire spiral. From the above description it seems reasonable, and in the following discussion it will prove invaluable, to think of the spiral antenna in terms of three different regions of operation - the transmission line region, where very little radiation takes place, the active region, where the predominate amount of energy is radiated, and the termination region, where the unradiated energy is absorbed. This division of the antenna into smaller, almost independent pieces allows certain parameters and effects to be separated out and better understood by themselves before trying to understand the operation of the antenna as a whole. The above discussion has, for clarity's sake, implicitly assumed a wire spiral, and its associated electric currents. It may be extended to include the slot spiral by applying the theory of duality, in which the two-wire line is replaced with a slot line, and the electric currents with magnetic currents. 3.0 DESIGN INFORMATION The design procedure for the slot spiral antenna assembly is somewhat long and involved, and requires the careful choice of many inter-related design parameters. In this section we summarize the concerns in the design of the slot spiral, and discuss the more important design parameters briefly. Some guidelines are also given to assist the designer in developing a useful antenna. The design parameters are: * spiral growth rate, a * number of turns, N * microstrip line characteristic impedance(s), ZAf* * slot line characteristic impedance(s), Zs. * substrate characteristics: r, tan 6, metalization technique,... * thickness of the substrate, t * length and impedance of the tuning stub, IStub & ZStub * cavity depth, h * cavity wall material * slot line termination material and technique 19

* specific geometry of the microstrip-slot line transition. As with most antennas, the design procedure of the slot spiral is an iterative one. However, because the individual parameters are highly dependent on each other, the design procedure is particularly involved. It is perhaps easiest initially to choose the growth rate, and then to carefully inspect the initial design for inconsistencies. For the following discussion, please refer to Figures 1, 2, and 4-6 for clarification. In the 500 MHz to 4 GHz region, the authors have found a growth rate of 65.4 mils/rad to be a good initial choice. To provide a good match at the connector, the initial characteristic impedance of the microstrip line is chosen to be 50Q. The characteristic impedance of the microstrip inside the spiral will be decided later, based on the inter-slot spacing and construction feasibility. At the center of the antenna, the actual feed is contrived through the use of a microstrip to slot line transition. Here, the microstrip line views the slotline as a pair of shunt branches. For the transition to operate most effectively, the characteristic impedance of the microstrip line is chosen to be half that of the slot line, yielding a perfect match at the feed. Based on the dielectric constant and thickness of the substrate, and on fabrication concerns, a certain range of slot line impedances will be available. For a dielectric constant of 4.5 and a thickness of 15 mils, the authors have found that a 100Q slot line, which has a width of approximately 20 mils, is convenient. This yields a microstrip characteristic impedance of 50Q, and a microstrip groundplane width of approximately 185 mils. To minimize coupling between the mictrostrip feed line and the slotline, it is best to maximize the ratio of microstrip line to ground plane width. Many design issues can be traded off in this case - the most important are the growth rate, the substrate thickness and dielectric constant, and the slotline and microstrip characteristic impedances. In this iteration just the microstrip line width will be modified - the others may be adjusted in later iterations. To maximize the width ratio mentioned above, the impedance of the microstrip line is increased until fabrication difficulties determine the minimum width. For safety, the authors don't venture below a 5 mil feature size. Klopfenstein tapers [7] are used in both cases to ensure that these impedance changes do not hinder the broadband behavior of the antenna and feed. To terminate the microstrip line following the transition, for expediency the authors currently use a A9 /4 open-circuit stub of significantly lower impedance than the microstrip line. More wide-band techniques are currently being investigated. The termination of the slot line, as discussed above, is critical for making the antenna behave in a frequency independent manner, and is the subject of ongoing investigation. Initially, tapered applications of foam and ferrite-loaded rubber absorbers were used, but were found to be of little help. Currently under investigation is the use of "radar-absorbing paint," applied in a tapered fashion over the slots. The design of the reflecting cavity is driven by the upper and lower frequency extents of the antenna. The reflector must be close enough to the antenna so that, at the highest frequency, it is still less than A0/4 away, to avoid complete cancellation on boresight. However, it cannot be so close that it shorts the fields in the transmission line region. In practice, these fields are very tightly bound to the slots, but it is still quite possible to significantly disturb them. To help avoid resonant cavity effects, the walls of the reflecting cavity are also made of absorber. As was mentioned earlier, the quality of the slot line termination is very important for low axial ratio designs and consistent pattern quality over frequency. The growth rate is also important, and should be kept as small as possible, again to minimize the axial ratio. The microstrip to slot line transition is important as well, particularly in terms of the overall bandwidth and impedance match of 20

the antenna. Finally, choice of the substrate is critical to minimize losses, especially in the microstrip line, and to reduce unwanted coupling effects. 4.0 MEASUREMENTS A slot spiral antenna, with parameters as discussed above, and operating frequency range from 750 MHz to 1.8 GHz, was contructed on 15 mil Rogers TMM-4, with 1 oz. electrodeposited copper on each side. Figures 4-6 show the front, back, and disassembed views of the antenna itself. Calibrated pattern measurements were conducted at the NASA LaRC low frequency antenna chamber, using a spinning linear source antenna. The antenna was mounted in a 4.5' square ground plane attached to the roll-over-azimuth-over-elevation positioner in the chamber, as depicted in Figure 7. The edges of the ground plane were also treated with tapered R-card to reduce diffraction effects. Orthogonal plane cuts were taken through and at 90 to the connector in an effort to determine and minimize radiation effects from the connector. Figure 8 shows a typical pattern cut, orthogonal to the plane of the connector, with a plot of the axial ratio as well. This pattern was taken with the slot line termination constructed from ferrite-loaded rubber and absorbing foam. Figure 9 shows the radiation pattern for the same antenna, but with a tapered application of "radar-absorbing paint" as the termination. The decrease in axial ratio is significant, and demonstrates the necessity of a good termination. Patterns similar to these were obtained over the entire operating range of the antenna, each showing a similar reduction in axial ratio. Figures 3, 8, and 9 also show an unfortunate lack of gain - the peak gain is approximately -2.5 dBic. Preliminary calculations indicate that most of this loss is attributable to losses in the microstrip line, especially since the copper on the substrate was electrodeposited. Techniques to minimize these losses are currently being studied. 5.0 CONCLUSIONS The design of the slot spiral antenna was discussed, and experimental results were obtained to examine its performance. A useful broadband slot termination is now being investigated with promising results, implying that the Archimedean slot spiral can be made truly frequency independent. Current and future work includes further numerical and experimental studies of the slot line termination, synthesis of a more broadband microstrip to slot line transition, and pattern performance evaluation when mounted on automobiles or other moving vehicles. ACKNOWLEDGEMENTS The authors would like to thank Mr. Robert Holt, Mr. Carl A. Lipp, and Mr. Hunter Walden for performing the measurements on the slot spiral antenna, and deeply appreciate the donation of the TMM-4 substrate by the Rogers Corp. 21

REFERENCES 1. R. Bawer, J. J. Wolfe, "A Printed Circuit Balun for Use with Spiral Antennas," IEEE Trans. Microwave Theory Tech., Vol. MTT-8, May 1960, pp. 319-325. 2. R.G. Corzine, J.A. Mosko,"Four-Arm Spiral Antennas," Artech House, 1990. 3. M.W. Nurnberger, J.L. Volakis, "A New Planar Feed for Slot Spiral Antennas," Technical Report 030601-6-T, The University of Michigan Radiation Laboratory, Aug 1995 4. V.H. Rumsey, "Frequency Independent Antennas," Academic Press, New York, 1966. 5. R. Bawer, J. J. Wolfe, "The Spiral Antenna," IRE WESCON, Pt. 1, 1960, pp. 84-95. 6. J.D. Dyson, "The Equiangular Spiral Antenna," IRE Trans. Antennas Propagat., Vol AP-7, Apr 1959, pp. 181-187. 7. R.W. Klopfenstein, "A Transmission Line Taper of Improved Design," Proc. IEEE, Jan 1956, pp. 31-35. (cm) - - - Spiral Slotline - Spiral Microstrip Balun/Feed Figure 1: Schematic diagram of the slot spiral antenna and microstrip balun/feed. Figure 2: Enlarged view of the Section AA', Figure 1, showing the antenna cross-section and one slot termination technique. 22

0 I I I I I I I I I I I I I I I I I I I I I I I I I I -5 -10 -15 co -20 0 -25 -30 -35 I I I i I I i I I I i I I I I i I ' i I I I I. 1.. -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 1211 Og249:r 6 th 3 -..................................... -..... ^........ '.........................................................,....................... - I i! I I I I I i I I i I i I I I I I I I i Ii I -120 -110 -100 -9() -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 Azimuth (deg) Figure 3: Spinning linear plot of spiral antenna, showing effects of poor arm termination. Figure 5: Bottom view of Slot Spiral Antenna, showing absorber walls and reflecting cavity bottom. Figure 4: Top view of Slot Spiral Antenna, showing infinite balun microstrip feed. Figure 6: Disassembled pieces of the Slot Spiral Antenna. Figure 7: Slot Spiral Antenna mounted on ground plane, in anechoic chamber. 23

0 1 1 1 1 1 1 1 1 1 1 1 1 1 f I I I I I I I I -- T I - 5..... -10 -15 - ~-A5 -30 _40 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 Azimuth (deg) Figure 8: Spinning linear plot of spiral antenna, with ferrite-loaded rubber and absorbing foam arm termination. I I I I I I I I I II I I I........ 5 -1.................................. -40 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 12. - — I I:.... -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 Azimuth (deg) Figure 9: Spinning linear plot of spiral antenna, with "radar-absorbing paint" arm termination. 24

Triangular prisms for edge-based vector finite element analysis of conformal antennas1 T. Ozdemir and J. L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, Michigan 48109-2122 August 14, 1996 Abstract This paper deals with the derivation of the edge-based shape functions for the distorted triangular prism and their applications for the analysis of doubly curved conformal antennas in the context of the finite element method (FEM). Although the tetrahedron is often the element of choice for volume tessellation, mesh generation using tetrahedra is cumbersome and CPU intensive. On the other hand, the distorted triangular prism allows for meshes which are unstructured in two dimensions and structured in the third dimension. This leads to substantial simplifications in the meshing algorithm and many conformal printed antenna and microwave circuit geometries can be easily tessellated using such a mesh. The new edge-based shape functions are first validated by computing the eigenvalues of three different cavities (rectangular, cylindrical and pie-shell). We then proceed with their application to computing the input impedance of conformal patch antennas on planar, spherical, conical and other doubly curved (ogival) platforms, where the FEM mesh is terminated using an artifical absorber applied conformal to the platform. Use of artificial absorbers for mesh termination avoids introduction of Green's functions and, in contrast to absorbing boundary conditions, a knowledge of the principal radii of curvature of the closure's boundary is not required. 'This work was supported in part by the U.S. Air Force Rome Laboratory and the NASA Langley Research Center. 25

1 Introduction The brick and tetrahedron are popular elements for a finite element analysis of electromagnetic problems. The first is indeed attractive because of its simplicity in constructing volume meshes whereas the tetrahedron is a highly adaptable, fail-safe element. It is often the element of choice for three dimensional (3D) meshing but requires sophisticated and CPU-intensive meshing packages. The distorted prism (see Figure 1) is another volume element which provides a compromise between the adaptability of the tetrahedron and the simplicity of the brick. Basically, the distorted prism allows for unstructured meshing (free-meshing) on a surface and structured meshing in the third dimension. An approach for growing prismatic meshes is illustrated in Figure 2 and most volumetric regions in antenna and microwave circuit analysis can be readily tessellated using such a mesh. As seen, once the surface grid over the platform is constructed, the volume mesh is grown along the surface normal by repeating the same grid at multiple distances from the platform. This avoids use of CPU-intensive volumetric meshing packages and, in many cases, including some popular patch shapes, the surface grid can be easily generated without resorting to a surface gridding package. Examples include printed rectangular and circular patch antennas, and circuits comprised of rectangular shapes. Moreover, because of their triangular crosssection, the prisms overcome modeling difficulties associated with bricks at corners formed by planes or edges intersecting at small angles. A special case of the distorted prism is the right prism which is characterized by the right angles formed between the vertical arms and the triangular faces [1]. The top and bottom faces of the right prism are necessarily parallel and equal, restricting them to a limited range of applications, namely, geometries with planar surfaces. In contrast, the distorted triangular prism is almost as adaptable as the tetrahedron with the exception of cone-tips which are not likely to occur in printed antenna and microwave circuit configurations. In this paper, we introduce edge-based shape functions [2]-[4] for the most general distorted prism shown in Figure 1. These prisms have non-parallel triangular faces and each of their three vertical edges can be arbitrarily oriented. In the following, we first present the derivation and validation of the edge-based shape functions and then proceed with their application for a characterization of printed antennas on doubly curved surfaces. In this appli 26

z A 3 PI y1 x Figure 1: The distorted triangular prism shown with the direstions of the edge vectors cation, the finite element mesh is terminated conformal to the platform's surface using an artificial absorber rather than an absorbing boundary condition (ABC) [5]. The employed conformal mesh termination is easily implemented by using prisms, and in contrast to the ABCs, the artificial absorber does not require a priori knowledge of the closure's radii of curvature or the wave's propagation characteristics. The utility and versatility of the proposed finite element method (FEM) formulation is demonstrated by considering the analysis of several printed antennas on different platforms. Specifically, we include input impedance computations for rectangular and circular patches on planar, spherical and conical surfaces. The radiation from a patch on an ogive is also considered. 2 Vector edge-based basis functions Consider the distorted prism shown in Figure 3. The prism's top and bottom triangular faces are not necessarily parallel to each other and the three vertical arms are not perpendicular to the triangular faces. The first step toward specifying a set of shape functions for the prism is the identification of a unique cross-section containing the observation point (x, y, z) (see Figure 3). This is done by introducing a parametric representation for the nodes (xi, y., z.) of this cross-section as illustrated in Figure 4. These parametric 27

Surface normal Surface geometric fidelity Figure 2: Prismatic mesh: Unstructured over the antenna surface, structured along the surface normal. Perspective view Top view 'N - 71 7 Interior point (x,y,z) Triangular cross-section containing the interior point Figure 3: Variation of the triangular cross-section along the length of the prism 28

(x3, Y4 z, ) 3" - X Y * O s 1Z21t ) I AA (x1,,y,,.z) (-,- i e > d - of, a3' (< sy3 n) d2 A A,- -I (xsy~z) no T fl (Xsey,.Z) ^ --->-v-^ y- (2^); > * d, i, cross-section - 1 \* 2 a s(r t-rl ) \ a)('s.^. ). - / -- - \ / I e/ a (x, Yh zb) ---— 4 (xd, y^, zb) 3 2 (a) (b) Figure 4: (a) Nodal coordinates, (b) triangular cross-section with the local coordinates b and t. representations are given in the Appendix. They involve the parameter s such that (xi, yi, zi) reduce to the nodes of the bottom triangle when s = 0 and those of the top triangle when s = 1. Given a point (x, y, z) interior to the prism, a unique value for s can be computed as discussed in Appendix A. Having specified the cross-section through the point (x,y,z), we next proceed with the derivation of the basis functions. We chose to represent the field variation across the triangular cross-section using the Whitney-1 form [6]. A simple linear variation will be assumed along the length of the prism. Specifically, the vector basis functions for the top triangle edges can be expressed as N1 = d ( L2VL3-L3VL2)s N2 = d2 (L3VLi-LiVL3)s (1) N3 = d3(LVL2-L2VL,)s and correspondingly those for the bottom triangle edges will be Ml = d (L2VL3-L3VL2 )(1-s) M2 = d2 (L3VLi-LiVL3 )(1 -s) (2) M3 = d3 (L1VL2-L2VLi )( 1-s). The subscripts in these expressions identify the edge numbers as shown in Figure 1 and the distance parameters d- are equal to the side lengths of the 29

A1 A( XYZ V " \ -W ^ "1/. A A \ -- ' \X.-y. — (a) (b) (c) Figure 5: (a) Vector map of N2 or M2, (b) Vector map of K1, (c) Variation of v as a function of x, y, and z. triangular cross-section containing the observation point (see Figure 4(b)). Also, 1 Li(^,7) = hi L( cos a3 sina3 L r} = -_ _- ___-r7 (3) h2.h2 L3(, ) =cosa2 + sin a2 h3 h3 are the usual two dimensional scalar node-based basis functions [7] for the same triangle with ai denoting the interior vertex angles and hi being the node heights from the opposite side. The variables g and rq represent the local coordinates and are illustrated in Figure 4(b). As required with all edge-based shape functions, N e& and Mi ej have unity amplitude on the ith edge whereas Ni ej = Mi ej =0 for i 7 j. Their vector character is depicted in Figure 5(a) and they simply "curl" around the node opposite to the edge on which their tangential components become unity. It remains to define the shape functions for the three vertical edges and we 30

chose to express these by the representations (linear over the cross-section) Ki(I, q) = v(, q) Li() K2( )= (,1)L2((,) (4) K3(g, q) = ) L3(((,). As before, Li are the node-based shape functions defined in (3) and a pictorial description of K1 is found in Figure 5(b). Of particular importance in (4) is the unit vector V(4, i). It is a linear weighting of the unit vectors ai, v2 and V3 associated with the vertical arms (see Figure 5(c)), and is given by = I V= EZ=1 L )L (5) This particular choice of v is oriented parallel to the side faces of the prism when evaluated on those surfaces and minimizes tangential field discontinuity across the side faces. Another choice is V(r s) = VS (6) which is always oriented normal to the triangular cross-sections of the prism and ensures tangential field continuity across the top and bottom triangular faces. Both choices are equally useful. However, (5) is more computationally efficient and has been used in generating the results presented in sections 3 and 4. The shape functions derived above for the distorted prism reduce to the right prism shape functions presented by Nedelec [8]. However, it should be pointed out that the above shape functions do not guarantee tangential field continuity across the faces of neighboring prisms. The possible discontinuity is primarily an issue for highly distorted elements and is caused by the fact that the horizontal vector basis functions (N and M) may have small nonvanishing tangential conponents across the vertical faces of the prism (and vice versa for the vertical basis functions). The expressions given in Appendix B neglect contributions due to the tangential discontinuities across the inter-element boundaries. These contributions were ignored because the discontinuity depends on the distortion of the prism and in practice the elements are marginally distorted. For general applications, our basis functions can be safely employed (within the accuracy of the finite element method) as 31

long as the deflection of the vertical arms is less than 5~. In the results to be presented later, the deflection was always kept less than 3~. Furthermore, in the case of patch antennas, the discontinuity is nearly nonexistant because the horizontal prism faces (top and bottom) coincide with metallic surfaces and therefore any possible discontinuity on the vertical fields under the patch (which are most dominant) is eliminated. Shape functions ensuring tangential continuity across the distorted prism faces can be derived but will be more complicated in form and CPU intensive to implement. On the contrary, the ones presented here lead to a more useful code for design purposes. 3 Eigenvalue Computation In this section, we examine the validity of the presented edge-based functions. Specifically, we consider the eigenvalues of three different cavities using the distorted prism as the tessellation element. We begin by first deriving the matrix elements following Galerkin's testing. The weighted residuals of the vector wave equation are IJPNx.(VxVxE- k2E)dV = 0, i = 1,2,3, (7) J JJMi (V x V x E-k E)dV=0, i=1,2,3, (8) Ki (v x x E - k2E)dV = 0, i = 1,2,3, (9) in which N-, Mi and Ki comprise the nine edge-based vector basis functions defined in the previous section and E is the electric field vector. The matrix equations are generated by introducing the representation 3 E = E N [ENNi(r) + EiMMi(r) + EiKKi(r)] (10) i=1 where EiN, EiM and EiK are the expansion coefficients, and correspond to the average amplitudes of the tangential field vector along the prism's edges. Substituting (10) into (7)-(9), and invoking the divergence theorem yields the element equations 3 3 EjN[NNCij - k2NND23] + EjM[NMCi, - k2NMDij] + j=1 j=1 32

3 E EjK-[NKCj - k2NKDj] =0 (11) j=l 3 3 EN[MNCij - k2MNDij] + E EjM[MMCij - k2MMDij] + j=1 j=l 3 EEjK[MKCi j - k 2MKDj] = 0 (12) j=l 3 3 EjN[KNCj - k2KNDj] + E Ejm,[KMCi j - k2KMDj] + j=l j=l 3 E Ej.K[KKC j - k2KKDIDj] = 0, i = 1,2,3. (13) j=l where the quantities NNCij, NNDij, etc. are integrals (over the volume of the prism) involving vector basis functions. Explicit expressions for the integrals are given in Appendix B. Upon assembly of (11)-(13) and boundary condition enforcement, we obtain the generalized eigenvalue system [A]{x} = ko[B]{x} (14) in which A = ko2 represent the eigenvalues of the problem. The matrices [A] and [B] are real, symmetric and sparse ([B] is also positive definite). 33

1cm Mode k, cm1 % Error - * 0.5cm (Exact) Prism Brick Tetra. - <- - - 1 TE 1I 5.236 0.73 -1.36 0.44 I * 0 0 TMI 11 7.025 2.32 -2.23 0.70 0.75cm ----------- - -- v i TE 1 [ 7.531 0.53 -2.58 1.0() - - -- TE211 7.531 0.64 -3.13 -0.56 (Actual mesh) TM _ _ 8.179 0.22 2.09 2.29 (a) (b) Figure 6: (a) Rectangular cavity discretized using right triangular prisms, (b) eigenvalues for the air-filled cavity and comparison with the brick and tetrahedron results. Example 1: Rectangular Cavity The first example is the rectangular cavity shown in Figure 6(a). Results based on brick, tetrahedron [9] and prism discretizations are given in Figure 6(b), where they are compared to the exact values. We observe that, overall, the data based on the prism are better than those based on the brick and, at least, as good as those associated with the tetrahedron. Figure 6(a) displays the actual mesh used for both the brick and prism discretizations. The number of degrees of freedom for the prism, brick and tetrahedron computations were 382, 270 and 260, respectively. Example 2: Dielectric ring resonator Figure 7a shows a dielectric ring placed symmetrically inside a cylindrical cavity and of interest is the computation of the resonance frequency in the presence of the dielectric ring. The resonance frequency was measured [10] using a loop antenna connected to a directional coupler as shown and the cavity was excited by the same loop antenna. Since maximum coupling to the cavity occurs at resonance, minimum power is returned to the detector at this frequency. For computation, the cavity was discretized as shown in Figure 7b resulting in 1051 degrees of freedom. The computed frequency was 1257 MHz (based on the smallest eigenvalue) and this is within 2% of the measured value (1282 MHz). 34

Dielectric ring Directional -" Coupler ~-W Metal Cavity -—,dl %% I,H I I, X% -4 (b) - I r 1 i iI.J4 11 - ~~ns~r$~b. II (a) Cavity dimensions: Height = 2.2cm Inner diameter = 7.5cm Outer diameter = 11 cm ~ r=13.5 Dielectric ring info: Height = 13.84cm Diameter = 15.24 cm Results: Measured = 1282 MHz Computed = 1257 MHz ( Error= 2 %) Figure.7: Dielectric ring resonator, (a) geometry, (b) mesh. Example 3: Pie-shell Cavity The third example is a pie-shell sector as shown in Figure 8(a) modeled using distorted prisms. It is obtained by bending the rectangular cavity considered earlier. The computed and exact eigenvalues for the first five dominant modes are given in Figure 8(b), and these testify to the accuracy of the distorted prisms in modeling curved geometries. The number of degrees of freedom used for this computation was 382. 4 Application to antennas 4.1 Finite element-artificial absorber method Prisms facilitate the modeling of conformal antennas since the presence of curvature does not complicate the mesh generation. As is the case with all PDE methods, the mesh must be terminated using a boundary integral, some approximate boundary condition or an absorber. When a Green's function is available, the mesh termination can be achieved using the boundary integral method right on the antenna surface yielding exceptionally accurate results. However, there are two problems with the boundary integral mesh trunca 35

(Actual mesh) Mode k, cm-' %Error _ (Exact) (Computed) P- ---- - - ---- TM 11or4.693 1.55 0.75cm TM210 6.009 1.28 TE1ll 6.640 2.71 TE 211 7.513 -0.27 T011 7.579 2.73 0.5cm - 1cm (70~) (a) (b) Figure 8: (a) Pie-shell cavity discretized using distorted prisms, (b) eigenvalues for the air-filled cavity. Ea= a= 1-j2.7 I ~h= ~ ( PEC < Eb= ~c(-J2.7) 0. 15 [ I - lj2 ^= IJ2.7 Es 9b = IJ2 *Rb-,b <Vo.5 15 1-j2.7) 2.7 (a) (b) Figure 9: Finite element-artificial absorber technique: (a) Antenna geometry, (b) Mesh termination using a metal-backed lossy dielectric layer. 36

tion techniques. First, the Green's function is only available for canonical geometries, i.e., planar [11], cylindrical [12] and spherical [13]. Second, the presence of material overlays can prohibitively complicate the derivation and evaluation of the Green's function making it very difficult to implement the boundary integral termination method. Thus, for modeling non-canonical geometries with possible material overlays, it is preferable to avoid use of the Green's function. Instead, an artificial absorber will be used for truncating the mesh as illustrated in Figure 9. The resulting implementation will be referred to as the finite element-artificial absorber (FE-AA) method and has several attractive computational features. Among them are fast convergence [5] and a capability to provide accurate input impedance computations. Accurate radiation pattern data can also be obtained when the rest of the platform has little or no effect. In its simplest form, this is done by ignoring the extent of the platform and integrating the aperture fields using the free-space Green's function. We next discuss the artificial absorber and then proceed with the application of the FE-AA method for the analysis of a variety of patch antennas on different platforms, some of which are presented here for the first time. 4.2 Design of the artificial absorber The absorber consists of a lossy dielectric layer backed by a metal. Metal backing enables us to terminate the mesh while the lossy dielectric lining traps the incoming wave and absorbs it, thereby forming a trasparent boundary. The absorber material parameters are chosen to completely eliminate wave reflections at normal incidence. Away from normal, the absorber reflectivity increases monotonically but can be minimized by proper selection of the absorber thickness and material parameters. Clearly, a thicker absorber has a better performance but carries additional computational burden. In addition, higher attenuation is achieved by making the layer more lossy. However, in this case more samples are required along the thickness of the layer to accurately model the field's amplitude decay. For a given thickness and sampling, an optimization can however be carried out. Such a study was performed in [5] and it showed that a minimization of the reflectivity for a 3-layer, 0.15A0 (free-space) thick metal-backed absorber yields the constitutive parameters of 6r = y,= 1 - j2.7 (see Figure 9). Below, we make use of this absorbing layer for mesh truncation. A layer based on the recently 37

Input Impedance vs. Frequency r = 1.3cm (0.15 WL at resonance) 200, =2.17cm (0.25 WL at resonance) 150 r3 = 2.6cm (0.3 WL at resonance) 0 r4 = 1.3cm (0.15 WL at resonance) -1 (a) (b) 50 0 0.8cm -50 ---- --- -2 3 4 5 <Av: Frequency (GHz) <. r2 Finite Element - Boundary Integral ^ (using tetrahedra) r4 Finte Element - Metal Backed Absorber / (using prisms) Mel = r- -j 2.7 1.3cm Resonance Frequency = 3.55 GHz i s (disagreement = 1%).. [ O5cm r= 2.4 Figure 10: Cavity-backed circular patch on planar platform, (a) Configuration, (b) Input impedance vs frequency and comparison with the finite element-boundary integral method. introduced anisotropic perfectly matched absorber [14] could have been used as well, but this was not found necessary at this stage. 4.3 Input impedance computations Cavity-backed circular patch on planar platform We first consider a patch on a planar platform since validation data can be generated using the more rigorous finite element-boundary integral method. Figure lOa shows a circular patch backed by a circular cavity recessed in a planar ground plane. A probe feed has been placed between the patch and the ground plane and Figure lOb shows the variation of the input impedance as a function of the excitation frequency. As seen, the computed impedance is in excellent agreement with the reference data obtained by terminating the finite element mesh at the cavity aperture using the half-space Green's function [11]. 38

For this example, the computational domain was discretized using 7,440 prisms resulting in 12, 523 degrees of freedom. One frequency run took 3.5 minutes on an HP 715/64 workstation. Microstrip circular patch on sphere Figure Ila shows a microstrip circular patch placed on a sphere. This is an example which clearly shows the advantage of the finite element-artificial absorber technique. Once the antenna aperture is triangulated, the volume mesh is grown along surface normals (using prisms) and terminated with the artificial absorber. Figure 1lb shows the comparison with the moment method [13] and the effect of sphere radius on the resonance behavior. Resonance frequencies predicted by this method and the moment method are within 2% of each other. Although negligible, the difference must be evaluated in view of the fact that the reference data is based on the approximation that only the dominant TM1, mode exists under the patch. Such assumptions are absent in the finite element analysis. Clearly, the value of the input resistance is a strong function of the employed feed model and therefore, it is not surprising to see differences in the levels of the resistance as computed by the finite element method and the single mode moment method solution. Figure lib also shows the resonance behavior for different sphere radii, and it is observed that the resonance frequency decreases with increasing radius. This trend might be explained by noting that the shortest distance between the radiating edges of the patch is greater for a larger sphere radius. The patch radius (measured on the spherical surface) has been kept constant for this calculation. Figures 12c shows the input resistance variation with the excitation frequency for different substrate/superstrate material constants and thicknesses. Similarly to the planar patch, we observe the downward shift in the resonance frequency caused by the increasing relative permittivity of the substrate. Notice also that the shift is only half as much when the superstrate is present even though the increase in the refractive index of the superstrate layer is 1.6 times higher than that of the substrate. This is because the field is much stronger under the patch than above it. As expected, Figure Iic shows that the loss in the substrate material has no effect on the resonance frequency but it lowers the overall level of the input impedance. The typical bandwidth increase with increasing substrate thickness is also observed. 39

Feed location - z x 200 cn E 0 u Ce., U0 0. 0 - 100 F Moment 7!" Method r=oo = I I (r=3.16cm) JTmt Metallic ground y x plane 0 2.5 3.0 (GHz) 3.5 (b) (a) 100 ^-1 0 I 50 (A.ca -@ O-O 1: ~ri=2.47, tan&-, t1=0.16cm, no superstrate 2: cri=2.8, tan= —O, t —=0.16cm, no superstrate 3: ~ri=2.47, tan —=O, t=0.16cm, ~r2=, r2=, t2=0.2cm 4: ~ri =., tanS=., t =0.16cm, no superstrate 5: ~rl =2.47, tan6=., tI =0.2cm, no superstrate 2.5 3.0 (GHz) (c) 3.5 e, (1-jtan6) (d) Figure 11: Microstrip spherical-circular patch, (a) Geometry, (b) Comparison with reference data and effect of sphere radius on the resonance behavior, (c) resonance behavior as a function of substrate/superstrate thickness and material constants, (d) definition of the parameters. 40

Metallic Ground AZ h=O.ll4cm 200.rb 12.9cm I 14cm 3.0 3.1 3.2 3.3 _ m (GHz) x~ (a) (b) Figure 12: Microstrip sectoral patch on cone: (a) Configuration, (b) change in the resonance behavior as a function of the cone angle 0. For this example, the computation region has been discretized using 6,048 prisms resulting in 9,997 degrees of freedom. One frequency run took 4.5 minutes on an HP 715/64 workstation. Microstrip sectoral patch on cone This is an example which illustrates the effectivenenss of the new technique to model antennas on doubly curved platforms with varying radii of curvature. To our knowledge, this is the first analysis of patches on such a platform. Figure 12a displays a sectoral microstrip patch printed on a conical ground plane and Figure 12b shows the antenna resonance behavior as a function of the cone angle 9. We clearly observe that the resonance frequency 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 12, 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 HP 715/64 workstation. 41

6" 6" A e (Side view) (Top view) (a) Feed 6cm location |<. * cm 2cm 5cm Cavity atch 1. 125cm r = 2.17 3cm 0.08cm - 0 - Polarized Radiated Power 0 -10 -20 -30 -40 - -50 ----- -180 -90 0 90 180 Observation Angle, (Deg.) Computation ------ Measurement (b) (c) Figure 13: Cavity-backed rectangular patch on ogive: (a) Setup, (b) antenna dimensions, (c) Comparison with the measurement. 4.4 Radiation pattern computations Cavity-backed rectangular patch on ogive As another example, we present the FE-AA method's application to a more general doubly curved platform. Figure 13a shows the set-up, where a rectangular patch is placed on the aperture of a rectangular cavity recessed in the ogive's surface. Also, in Figure 13b we show the dimensions of the antenna. The radiated field is computed by integrating the tangential aperture fields using the free-space Green's function. Thus, the shown pattern accounts only for the antenna (and the curvature of the platform), but does not include interactions with the ogive's tips. Figure 13c shows the com 42

puted 0-polarized radiation as compared to the measurement [15]. Clearly, the agreement is very good for this polarization. However, predicting the 0 -polarized radiation (not shown) requires modeling of the entire ogive as this particular polarization has a vertical surface field component which is known to cause diffraction from the ogive's tips. A way to account for such secondary diffractions is to interface the FE-AA method with a high frequency technique and an encouraging study in this direction has been carried out in [16]. 5 Conclusions A new finite element technique was introduced for the analysis of printed antennas recessed in platforms of non-canonical shape. The distorted prism was introduced as the volume discretization element, and corresponding edgebased shape functions were derived and tested for eigenvalue computations. A major part of the paper though was devoted to characterizations of printed antennas using the new technique. Use of prismatic elements was found very attractive for this application and was essential in simplifying the modeling of antennas on doubly curved platforms. An artificial absorber was used to terminate the mesh conformal to the platform thereby minimizing the size of the problem while preserving the sparsity of the finite element matrix. Use of the absorber also avoids difficulties associated with the conformal application of the classical ABCs. Antennas on spherical, conical and ogival platforms were considered without using a Green's function and therefore, the superstrate/substrate metarials can be readily accounted for. A limitation of the technique is that it models the immediate neighborhood of the antenna, and therefore ignores the details associated with the rest of the platform and possible substructures around the radiator. However, the proposed method was shown to be a good approach for predicting the resonance behavior of antennas, and could evolve as an important tool for designing conformal antennas on doubly curved platforms. 43

Appendix A Referring to Figure 4a, a way to specify the nodes of the triangular crosssection through the point (x, y, z) is to use the parametric representation ri = rib + s(rit - rib), i = 1, 2, 3, (Al) where.ri are the nodal position vectors of the triangular cross-section (see Figure 4a) whereas rit and rib are associated with the top and bottom triangular faces, respectively, and 0 < s < 1. Thus, for s = 0 r, specify the nodes of the bottom triangle and for s = 1, ri reduce to the nodes of the top triangle. The cross-section sweeps the entire volume of the prism as s assumes values between these two limits. Given an interior point (x, yz), the task is to specify the cross-section of the prism that contains the point. In mathematical terms, a solution of s in terms of x, y, and z is required. On our way to finding this solution, we note that the normal distance of the point (x, y, z) to the cross-section must be zero, viz. (x -.xi) nX + (y - yi) ny + (z - Z1) nz = 0 (A2) where x1, yi and z1 are the coordinates for one of the nodes of the crosssection (see Figure 4a), and nx, ny and nz are the components of the crosssection's unit normal pointing toward the top surface and can be computed as ( ) (r2- rli) x (r3-r2) s) (r2-rl) x (r3-r2)|I where ri are as given in (Al). The left hand side of (A2) is a linear function of s for the right prism, and is nearly linear for a marginally distorted prism (which is the case in practice). Therefore, a very fast converging solution of (A2) is obtained by the Newton-Raphson method. In generating the results presented in sections 3 and 4, an average of three iterations sufficed to achieve a tolerance of 10-5 x (Prism Height). 44

Appendix B For the assembly of the finite element equations, the following quantities must be computed: NNCit = (V x Ni) (V x Ne)dV NMCe = J J I (v x N) (v x Me)dv NKC- = (VxNi) (V xKe)dV MMCie = I (V x Mi) (V x Me)dV MKC = 1 J (V x M.) (V x Ke)dV KKCe = JJJ(V x K) (V x Ke)dV NNDie = J Ni NedV NMDit = I J Ni MedV NKDie = J IJ Ni KedV MMDie = JI Mi MedV MKDie = IJ Mi KedV KKDie = Ki KedV. where the integrals are over the volume of the prism and they must be evaluated numerically. However, numerical integration over the distorted prism volume is cumbersome and therefore, the distorted prism is first mapped onto a right prism as shown in Figure B1. The integration over the new volume is then carried out by using the five point Gaussian formula. The relationship between the two integrations is given by i J f(x,y,z)dxdydz = J J f(x, y,z) det[J] ddr/dC, Yr J 77C- / ' i 45

where (x, y, z) are related to ((, /, () through the bilinear transformation x= 4 + C(alj + a2 + Xlt) (B1) y = r/ + ((blj + b27r + lt) (B2) Z = ((C1( + C27 + Zlt) (B3) with X2t - Xlt - X2b al = X2b X3t - Xlt - X3b - alX3b a2 = - r Y3b Y2t - Ylt bl = X2b Y3t - Y1t - y3b - blx3b b2 = Y3b Z2t - Z1t Cl - X2b Z3t - ClX3b - Zlt C2 - Y3b and (Xib, Yib) and (xit, yit, zit) are the coordinates of the prism's bottom and top surface nodes, respectively, as shown in Figure Bi. Also, [J] is the Jacobian of the bilinear transformation defined by (B1)-(B3) and is given by [J]- gg g 6 I -a ay a ax an az ax a: a Special Case: Right Prism For the right prism, closed form expressions can be derived for the above integrals. Referring to Figure 4b, we have C didt cOS 3kn cos jm - COS 3km n ENNCi- (,, Xjm + Xhkn - Xjnh c hkhn hjhm hkhm X cos/3 j 2 c2hid1 Thh~ Xkm + - hjhkhmh Sztn/jk SZmn ) n3i 0 ifl3~ 46

(x 2t ly2tl Z 2t ) (x 3b, Y3b' 1) (x It ly it, Z, t ) (0, 0,1) Bilinear transformation x (x~b 0, 0) (0, 0, 0) (x2b',O,) (0, 0, 0) (Distorted prism) (Right prism) Figure B 1: Mapping distorted prism onto a right prism ENMCie7 - j C -CS k Xj1m _ cs /jmXkn + hjhm CSkXjn + hkhm Cos/3j' 1 c2hid, h~~Xkm + - szn!%kszflI~mn ) 3 hh~~~ ENKC-iI EMMCieI EMKCje EKKC-tz ENNDie2 ENMDje2 hid, Cos/31e - 6 h -he COS /3ke hkhe = ENNCie - -ENKCje h1d, cOs/13ie = 2 hih1 dijCSO, COS j COS /3km dd o /k j + - Xkn - hkh 3Xn -C hkh /-3 hj-hm h 1 — ENNDje 2 c/3 Xkm) h3h 47

ENKDie = 0 EMMDie = ENNDie EMKDni = 0 EKKDie = cXie where c is the height of the prism, 0 O ifr=s 3rs r + Cs otherwise Xrs = hid {WrWs + - [(cotCt3 - COtc2)(,sWr + 7?rW) + 2(sWr + rWs)] + h2 1 [3(cota3 - cota 2)(Tq,r + iri3s) + 2rrs (Cot2 t2 - cota2cota3 + 12 cot2a3) + 6rs]}, r, s = 1, 2,3, the set {i, j, k} takes on the value {1, 2, 3}, {2, 3, 1} or {3, 1,2}, and the same is true for the set {/, m, n}, and sin ca3 W1 = 1 Cos= -e 12 = ~ w2=0 -- 2 —o 2= ~ = h2 = h W3 = 0 =sin a2 h3 h3 48

References [1] Sacks, Z., S. Mohan, N. Buris and J. F. Lee, "A prism finite element time domain method with automatic mesh generation for solving microwave cavities," IEEE APS Int. Symp. Digest, Vol. 3, pp. 2084-87, Seattle, Washington, June 19-24, 1994. [2] Nedelec, J. C.,"Mixed finite elements in R3," Numer. Math.,Vol. 35, pp. 315-341, 1980. [3] Bossavit, A.,"A rationale for 'edge-elements' in 3-D fields computations," IEEE Trans. Magn., Vol. 24, No. 1, pp. 74-79, January 1988. [4] Webb, J. P.," Edge elements and what they can do for you," IEEE Trans. Magn., vol. 29, no. 2, pp. 1460-1465, March 1993. [5] Ozdemir, T. and J. L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Sci., Vol. 29, No. 5, pp. 1255-1263, Sept.-Oct. 1994. [6] Bossavit, A.,"Whitney forms: a class of finite elements for threedimensional computations in electromagnetism," IEE Proc. Part A, Vol. 135, No. 8, November 1988. [7] Zienkiewics, 0. C. and R. L. Taylor, The Finite Element Method, 4th ed., Vol. 1, McGraw-Hill, New York, 1989. [8] Nedelec, J. C., "A New Family of Mixed Finite Elements in R3," Numer. Math., No. 50, pp. 57-81, 1986. [9] Chatterjee, A., J. M. Jin and J. L. Volakis, "Computation of Cavity Resonances Using Edge-Based Finite Elements," IEEE Trans. Microwave Theory Tech., Vol. 40, No. 11, pp. 2106-2108, Nov. 1992. [10] Muldavin, J. B., A. D. Krisch and M. Skalsey, Personal communication, University of Michigan, 1995. 49

[11] Gong, J. L. Volakis, A. C, Woo, and H. T. G. Wang, "A hybrid finite element-boundary integral method for the analysis of cavity-backed antennas of arbitrary shape," IEEE Trans. Antenna Propagat., Vol. 42, No. 9, pp. 1233-1242, Sept. 1994. [12] Kempel, J. L. Volakis and R. J. Sliva,"Radiation by cavity-backed antennas on a circular cylinder," IEE Proc.-Microw. Antennas Propag., Vol. 142, No. 3, pp. 233-239, June 1995. [13] Tam, W. Y., A. K. Y. Lai and K. M. Luk, "Input impedance of spherical microstrip antenna," IEE Proc.-Microw. Antennas Propag., Vol. 142, No. 3, pp. 285-288, June 1995. [14] Sacks, Z. S., 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, pp. 1460-1463, Dec. 1995. [15] Sliva, R. J., and H. T. G. Wang, Personal communication, Naval Air Warfare Center, Weapons Division, China Lake, California. [16] Ozdemir, T., M. W. Nurnberger, J. L. Volakis, R. Kipp and J. Berrie,"A hybridization of finite element and high frequency methods for pattern prediction of antennas on aircraft structures,"1EEE Antenna Propagat. Mag., Vol. 38, No. 3, pp. 28-38, June 1996. 50

Comrnparison 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 0(N'15) or 0(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. 51

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 computationally 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 NlogN [1],[2]. However, in general, the boundary integral is not convolutional and in that case the CPU requirements will be of O(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(N'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 I f Oil f Oin I [ [ ] [BP1] { I } } { {^nc} } (0) The vector and typical form of this system is given in Figure 2. {4} represents the magnetic field at the nodes within the groove and on the boundary. Also, {01} 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 [F1] is a full sub-matrix and thus 0(Nb) operations are needed to perform the product [P1]{'1} 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 [1Pi] { I }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 FMIM is a fast method for calculating the matrix-vector product. The computation of the matrix-vector product is illustrated with the boundary integral for TE 52

I A y Ad! ||-pec x W 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 H + EM (2) with Hs the scattered magnetic field, Hz the incident magnetic field and HFEM the total field which is employed to enforce tangential continuity and Hs =- MZ(P')I42o (koP - p) dl' (3) Mz = EX is the magnetic current and FI 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 tb groups are constructed. Next, by invoking the addition theorem for cylindrical wave functions, the Hankel function is expanded as Q/2 H(2) (ko pIW + p-P ) Jn(ko pl - )H( (kop (-p ) (4) n=-Q/2 where qil and O()pip are the angles pI' and p' - - make with the x-axis respectively and is valid for pl' > I P - 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 In(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 J(k ) = eko(p'-)-j'n(O-pp+/r2)d (6) 27r 27 53

0 ': ' h.h 60 " \ "'. \ 60 BI System 80 h 100: 120 L ' 1 0 20 40 60 80 100 120 column Figure 2: Typical form of the FE-BI system matrix arising from the scattering/radiation problem of a groove in a ground plane. in conjunction with (4) is used to write (3) as H ) k0Y0 V / -(o)TI()e-koPd (7) where the far-field pattern of the source group is given by V,()= j Mz(p')eC' Pdl' (8) and the translation operator Tw,(q) is given by Q/2 ()= E H(2)(kopIIt)e-(U-'+l/2) (9) n=-Q/2 with the propagation vector of the plane wave k0 = ko(x cos & + y sin (q) (10) The integral over q 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 o T4 r,(~)V1,()e5-iq'P (11) 4wr q=1 54

where AO = 27r/Q indicates the angular spacing between propagation vectors of plane waves emanating from a group and thus Aq = qAq q = 1... Q and kq = ko(x cos bq + y sin iq). 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 &. %#-o ' I Group P,, rop Group Group,'Group, I KNb d K! d i K -i -z x >d mi.' P Far group | P/ | < d min 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 V ('-') = M( I')e dl' (12) Evaluating this pattern for a single source group and in a single direction requires Mb 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, Al(Oq) = M^(0q)T111(0) (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,,Op ' = - AX Al(q)e-i (14) q=1 Thus, computing the redistributed field at a s in t requires Q operations, and for Lb groups each containing Mb unknowns the operation count is QLbMb. 55

The total operation count in the above sequence is CiQMbLb+C2QL2. By choosing, Q Mb N2 for convergence, the operation count is given by ClMbNb + C2 b-. and on setting Mb Nb, the final operation count is NJ5. 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 (lifferent 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 sinc 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 will (Oq) = f 1 (q ( - k-11'I < s) { 6c~(kqkj1,/i~ (I (k - q011i1 > iOs) (15) with s = sin- Q -m) (16) \2kopil) and ao is the taper factor. The discretized plane wave expansion is now H- = - A ~~ ~ (e-'T (17) q=1 The operation count now associated with (13) is now reduced to C3L NM2 with the cor2 2 Mb responding total operation count given by CMbNb + C4-2'. Grouping the unknowns with b 56

N1/3 per group, results in an operation count of 0(N4'13). 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 X.. -...I Groups 1 and 3 250 f- Q=54;P,1=5k (? — Q54; Pl5X - 2.00 ff0 00=1 ^ f; \ \ --— Q=54;P, /45?i F e 4: Te G oups 1 and 8 ~130 g v Groups 1 an o 27 1.00 - *^ AAI,' 1 'I VA.50I r 11 j jit I~A A -180 -120 -60 0 60 120 180 - Off(degrees) Figure we 4: The Translation operator for different groups on the boundary of a 50A wide groove; 750 BI unknowns; 27 groups 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 the addition theorem, in this algorithm the large argument approximation H' k) - p1 -jkopit1.7 2 6-jkoPjkop^,.j (18) w Vlkopp1 is used, where plil is the distance between the center of the test group I and the center of the source group 1'; pni' 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 57

-I ff=1I Group L Group Group./Group b d K+ I /' K-I nWmdow - - A y z x z x I I I I.' I I IP F> gdmin Far group:, p,; I P,< dmin# Near group 'ff 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 __ T/i 'I = t|-j kop~ I t lpnl vi/I = Z j16S j.= (19) 2. Since the above aggregation operation needs to be done for all source groups the operation count becomes O(NbMb) - O(Nb), where j- 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 N- test groups leading to a total operation count N 2 of O( b ) for aggregation. 3. The next step is a translation operation corresponding to Al/, = TIIVl where (20) (21) 2j ejkopili V 71 lc=pui Since this needs to be done only at the group level, it involves O( j ) operations for all b possible test and source group combinations and is the least computationally intensive step. 58

,1I Group, n[ / Group Group, /Group K+ 1 Kal K- I N.,,m, i 1 /e I 0 / I -it / I / \ 1 V^ \ / I I | I.I I I / I —. I- - " -4-1 - ^ IP >dmin#. m Far group Pij| < d mm > 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 ko y Nb/Mb H(p) = - 2 E All,-JkoP m (22) 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( Nb ) 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!^b test groups, the operation count would be of O( ). 5. The near field operation count being of O(NbMb) and the far field being O(^ ) gives a total operation count of NV2 Op.count - C1NbMb + C2 M (23) Typically, we can set the elements in each group, Mb = Nb and as a result the total operation count is O_ NVb 5. 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 59

and 50A. The depth of the groove was 0.35A with a material filling of Cr = 4 and Ur=1 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 (>1 dB). If the maximum tolerable RMS error is set at 1 dB [9], the FE-Windowed FMM 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 Total BI Width Unknowns Unknowns FE-BI (CG) FE-FAFFA FE-Exact FMM FE-WFMM 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) Is I -s - 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 compared 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. 60

10~ 2 \-2 10 -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, J. Gong, and A. Alexanian. A finite element boundary integral method for antenna RCS analysis. Electromagnetics, 14(1):83-85, 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, August 1990. 61

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 [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. M. 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. 62

30 A y A. pec x 20 o X-10 -20 -30 94 w= 10k (a) 0 30 60 90 120 Observation angle (deg.) (b) 150 180 RMS error (dB) Groove Width FE-FAFFA FE-Exact FMM FE-WFMM 10A 0.7627 0.1621 0.3291 (c) Figure 8: Scalability of the hybrid techniques to smaller problems (a) Problem geometry (b) Bistatic patterns (c) Error table 63

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. 64

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 Ez or Hz = e-ko(x cos+/sin ) (1) the reflection coefficient is RE(0) or RH(O), and the objective is to minimize these. If the layer is composed of a homogeneous isotropic material whose relative permittivity Cr and relative permeability pi are such that,r = /r = b = a - j/3 (say), then R/1 - b~ sn- 2 sin2 - j cos 4)tan(kobt /1 - b-2 sin2 d) RE (2) = /1 - b-2 sin2 q + j cos ~ tan(kobtl - b2 sin2 (2) (/1 - b-2 sin2 - j cos X cot(kobt/l - b2 sin2 ) (3) RHW (3) 1 - b-2 sin2 + j cos $ cot(kobt/l - b2 sin2 4) These differ because the presence of the metal backing has destroyed duality, and at grazing incidence (- = 7r/2), RE = RH = -1. If sinO > 1, i.e. - = 7r/2 + j8 with a > 0 so that sin c -cosh 8 and cos s =-jsin, IR differs from unity by only a small amount for both polarizations. The behavior of IRE,H(~)I as a function of sinkb is illustrated in Fig. 3. At normal incidence (4 -= 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 l3 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 oc, RE,H(q) - 0 only for normal incidence, but Sacks et al [4] have shown that a particular uniaxial anisotropic material has this property for all ( < 7r/2. The result is an example of a perfectly matched layer (PML), and if r = r = b - (b - )x (5) where I is the identity tensor and b a - j/3, then RE,H(S) = Te-2jkt(-j3)cos (6) which reduces to (4) in the particular case of normal incidence. If ko/3 >> 1 the reflection coefficients decay exponentially for all ( < 7/2, and since (6) is also valid for sink > 1, the choice a > 0 ensures an exponential decay for these angles as well. The behavior of IRE,H(k)| is illustrated in Fig. 3 for the same values of k0t, a and Q3 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. 65

3 Absorber Design The theoretical behavior of the reflection coefficient IRI for the metal-backed layers is relatively simple. In the case of the isotropic material, an increase in 3 and/or t decreases IR(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 = 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 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 IR\ 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)I 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 /t/oA0. 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 /t/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) and N versus /3t/Ao on the same figure, and these can be used to determine the optimum N for a given /t/Ao. 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 /t/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 66

the element size or the layer thickness. Thus, for elements 0.025Ao thick, we have t = 0.25Ao and y3 = 2.32. By increasing N we can improve the performance up to the limit provided by the theoretical value of IR(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.01061R| + 0.0433 (7) Ao N = 0.147exp [7.353/3t/Ao] (8) where IRI 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 3 = 2.6, and thus 13 = 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 (7 =: 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 67

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 l.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.Namely, 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. 68

Eb= c(l-J2.7) b= 1-j2.7 I ~a= a= 1-j2.7 0.15 X Es Eb= E(1-j2.7) b = 1-j2.7 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, 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. [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. 69

ty ti 1... N x x x Figure 2: Geometry ofl the metal backed absorber layer and its finite element discretization. 0 -10 -20 -30 -40 -50 - -60 -70 -80 0.0 0.2 0.4 0.6 sin o 0.8 1.0 1.2 Figure 3: 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. 70

0 -10 -20 -30 -40 -50 -60 -70 (.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 z 15 10 5 -80 ' 0.3 0.4 0.5 0.6 0.7 Pt/?o Figure N: ( 5: Absorber design curves. The -) exact IRI; (- — ) numerical straight lines give IRI IRI in dB, and the curved ones give 71

Metallic Ground Plane z t Patch I -// h= 0.114cm pEr = 2.32 rb = 12.9cm r = 10cm DO= 18.34~ 200 C/ 100 0 cl 3.0 3.1 3.2 3.3 (GHz).f — o 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 72

(a) RCS - NASA METALLIC ALMOND 0 20 40 60 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. 73

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 (FEMNI) 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 74

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 terminating 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 rnietal-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 compared 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. 75

'V 1y. 1... N x lo x 0 (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 E or Hz -jko(xcosq +y sin q) (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 (e: and relative permeability /Ur are such that -r -= 1t, = b = a - jf (say) where a and 3 are real, then RE() /- b-2 sin2 - j cos d tan(kobt\ - b- 2 sin2 4) y'1 - b-2 sin2 + jcos tan(kobt/1 - b2 sin2 () RH1 - b-2 sin2 -j cos 4cot(kobt - b2 sin2 ) 1-b2 si2 - b2 sin2 ),/ -6b-2 sin2 + j cos ( cot(kobt/l-b-2 sin2 +) (2) (3) These differ because the presence of the PEC backing has destroyed duality, and at grazing incidence (~ = 7r/2), RE = RH = -1. If sine > 1, i.e. ( = 7r/2 + jS with 8 > 0 so that sinq = cosh and cos q = -j sin, JRI differs from unity by only a small amount for 76

0 -10 -20 -30 -40 / -50 -60 -70 -80, 0.0 0.2 0.4 0.6 0.8 1.0 1.2 sin q 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(&)I as a function of sin $ is illustrated in Fig. 2. At normal incidence (- = 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 kot13 sufficiently large. As kot -- oc, RE,H(4) -- 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 Er =, = b- (b - b)xx (5) where I is the identity tensor, then RE,H(0) = e-23jko( - (6) which reduces to (4) in the particular case of normal incidence. If ko3 >> 1 the reflection coefficients decay exponentially for all q < 7r/2, and since (6) is also valid for sinkb > 1, the choice a > 0 ensures an exponential decay for these angles as well. The behavior of IRE,H() I is illustrated in Figure 2 for the same values of kot, a and / used for 77

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 -= {xy/(x)} a wave propagating into the material has the form e-jko {x/(z)cos O+y sin }(7 (7) and when (x)= 1 + (a - j/3 - 1) (t), (8) which tends to unity as x -+ 0+, the reflection coefficients of the layer are identical to those given in (6). With this expression for 7(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 \RI is relatively simple. In the case of the isotropic material, an increase in /3 and/or t decreases JR(0)l. 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, /3 and t, and our objective is to find the minimum number of sampling elements (or discrete layers) to achieve a specified IR[. It is anticipated that the errors introduced by the discretization will have a number of consequences. In particular, 78

0 -10 -20 -40 -50 -60 -70 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, (- - -) N3, ( — -) N=6 and (-) N=12. for a given number N of discrete layers and given t, increasing: will ultimately lead to an increase in IRI 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, 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 ciscretized 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 [R(0)I converge to the common theoretical value for both polarizations. 79

3.2 Dependence on ac and Q For a layer of constant thickness the theoretical value of IR(0)| is independent of a and polarization, but in the numerical implementation the behavior is much more complicated. Figure 4 shows IR(0)I plotted versus a and f3 for a layer of thickness t = 0.25Ao made up of 5 (=N) elements, where the darker tones indicate lower values. For small / the results are in close agreement with theory. As evident from the level lines, |R(0)| is almost independent of a and decreases exponentially with /3, leading to a linear decrease on a dB scale. For large f3, 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 / 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 /3 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 /3 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 /3 = 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. 80

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

0 -10 -20 -30 -40 -50 -60 0 1 2 3 4 5 Figure 5: |R(0)| with t = 0.25Ao and N = 5: ( —) Exact, ( -) a = 0 E pol., ( ---)a 0 Hpo., ( ---) a = 0.75 Epol. and (...) a = 0.75 H pol. 3.3 Dependence on /, N and t We now seek a connection between the values of /, N and t for which IR(0)I is minimized. To this end, we first examine IR(0)I as a function of / and N for constant t, and the resulting plot is shown in Figure 6 for E polarization with t = 0.25Ao 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 3. Consider now the behavior of IR(0)| for fixed N. As /3 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 3 at which [R(0)I 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 JR(0)| the curve of 3t/Ao versus N is virtually the same for all thin layers. The observation that 3t/Ao 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 82

I 1 0 10 4 - a6 | 1 1 -20 3N 1-30 IRI [dB] -40 0. The dashed curves are level lines and the solid curve gives the location of the minimum |R|. and N 3, then 2.13. If a smaller value of is chosen, -50 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. I.f 3 is set to a value larger than 2.13, (0). will still be 0 l.,............. -6 0 2 4 6 8 10 12 14 16 18 20 N Figure 6: JR1 as a function of 0 and N for E pol. with t = 0.25Ao and larger because the chosens aretoo smallevel linesto reproduce the rapolid curve gives the decay within the minimum layer. So far we ave considered only nor 2.13. If a smal incidence, but fo is chosen, IR() will be larglayer it is a simple Figumatter to extend the results to all real angles of inthat tcidence field reflected from the metal backing has not been absorptionenuated sufficieat normal incidence is set to a value larger thany angle if the layer thickness is be inver becausely proportional to coN is too small to reproduce thspeifying the decayer thickness the layer.wavelength along the normal (or x axis) to the considerial nterface. Sincidence, but for the anisotropic layeris nowt ndeis a simple matter to extend the scalablesults to all real angles of in-ormal icidence) become. As evident from the anisotropic layer, the absorption obtained at normal incidence is reproduced any angle if the layer thickness is by inversubstituting ely proportional to cos. This can be achieved by specifying theas a layer thickness t as a fraction of N duplicates the cur waveles shown in Figure 7. This notion can also be used to he maccunt for problems wherSince the outer cos not free space. In suc h cases, we have parameter t (at normal ncidence) incidence) becomes Ot/AX. For the anisotropic layer, all the results obtained at -normal incidence are made applicable for arbitrary 0 by substituting Ax for Ao. For example, plotting the optimum OtlAx as a function Of Al 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 -_oI/V (at normal incidence) 83

0.7 0.6 0.5 0.4 0.3 0.2 E. E I, -, 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 cr is the relative permittivity of the outer medium. Although the scaling property of /3t/A\ has been established only for oa = 0, it holds to a reasonable degree for small a -7 0, but as a increases, the 3t/A) versus N curves become increasingly dependent on a. The scalability also extends to the associated values of |R(0)[, 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)| are the same for layer thicknesses up to about 0.5A, at least, design curves can be obtained by plotting I|R() | and N versus 3t/A, 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 JR(0)[ curve intersects the -50 dB line at f3t/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)| which has been included in Figure 8. A good approximation to the short-dashed curves in Figure 8 obtained by linear 84

30 -a ' - 25 ---- //- 25 JRI -40.__ _ _ _2 / / 20 -50 '-....'. ~.. / / -15 -60 IN/ 10 -70 -80 0.3 0.4 0.5 0.6 0.7 0.8 Pt/Xx Figure 8: Absorber design curves. The straight lines give IR\ in dB, and the curved ones give N: ( — ) exact IRI, ( — -) homogeneous case and (- -) inhomogeneous case. regression is /t = -0.01061RI + 0.0433 (9) Ax N = 0.147exp [7.353/3/Ax] (10) where Ax = A0o/ cos A, 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 (, 0~ < q5 < 90~, for the anisotropic layer, and / = 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 -(x) = -j/3 (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 3t/Ax c_ 0.64 corresponding to N = 9. Approximations 85

to the long-dashed curves in Figure 8 are -O.01l91RI + 0.0451 (11) Ax N = 0.298exp [5.2633t/Ax] (12) where IRI and N are as before. Compared with the homogeneous material the decrease in the number of elements required becomes more pronounced as \R\ 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 Or = 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,r ir in the upper section and 3.2Er, /, in the lower section to match the substrate. The calculations were carried out at several frequencies using an FEMI 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 3 is shown in Figure 9. Recognizing that most of the field is confined to the substrate, we have Ax = A0o/ / = 0.559Ao. Using t = 1.789Ax and NI = 5 in the design formulas (9) and (10) gives /3 = 1.07 and 86

. Uniaxial absorber 3 095 cm 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 IR[ for the microstrip line at 4 GHz with c = 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 /3 - 1 and \R\ -- -50 dB. The fact that the minimum IRI 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 IRI, 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-j/3. 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 O3t/AX = 0.574 and IR] = -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 /t/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 ca of the material parameter b was set to 1, 87

-20 40c -30 4 cm -40 / 2.215 cm L Uniaxial Absorber -50 - ---------- ( --- 4.775 crnm 0 1 2 3 4 (a) Pt/x (b) Figure 10: (a) Geometry for the rectangular waveguide, (b) computed \RI 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 Z 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 88

(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. c =- 0, in the layer. Condition (iii) is a requirement for scalability, and though small values of ac 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. 89

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. 90

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 91

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:3 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 PMVL absorber's performance can be characterized in terms of the product Ag (a scalable quantity when a = 0) and the discretization rate. A two-dimensional analysis was carried out to determine the optimum values of 3 and N (the number of samples in the PML layer) for maximum absorption near normal incidence. It was determined that given a desired reflection coefficient IRI for the PML absorber, the optimum t and N values are approximately given by the expressions [7] /3t = -0.01061R| + 0.0433 Ag N = 0.147exp 7.353 -p1 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 IRI equal to -50dB, from the above formulae we have that 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 92

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 Aft _ 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 S1d 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. M[ittra, "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. 93

F F 1 1;.......r....:::I'I^^^^^^X II:"i 11lrb Prb 4r, Figure 1: Illustration of wave incidence upon a perfectly match interface (PML) with and without metal backing. Uniaxial '"'-Artificial Absorber Cavit y j 3.0953cm z4x 0.21cm x Strp Width = 1)/10 =0.71428cm Az=0.21 cm Ax=0.1 1905cm Ay=0.1 3636cm 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 absorber layer. 94

Figure 3: Illustration of the field decay pattern inside the PML layer. ol I I -5 -10 X-15 2 -20. -30 -E -35 -40 -45 -— f=4.0GHz --- f=4.5GHz - f=5.0GHz?-Y -50, I I 0 0.5 1 1.5 2 2.5 3 3.5 21t in Kg (a=1.0) 4 Figure 4: Reflection coefficient vs 2,3t/Ag with a=1, for the shielded microstrip line terminated by the perfectly matched uniaxial layer. 95

I I I 1.525 mm -..61 mm It Coi -bY1 Figure 5: Illustration of a meander line geometry used for comparison with measurement. -5 -10 -15 -20 -25 x — Measured x Calculated -30 -35 -401 6 8 10 12 14 16 18 20 22 24 26 frequency (GHz) Figure 6: Comparison of calculated and measured results for the meander line shown in Fig.5. 96

HYBRID FINITE ELEMENT METHODS FOR ELECTROMAGNETICS Applications to Antennas and Scattering JOHN L. VCLAKIS, JIAN GONG AND TAYFUN OZDEMIR University of Michigan Ann Arbor, Michigan U.S.A. Abstract This chapter provides an overview of the finite element method (FEM) as applied to electromagnetic scattering and radiation problems. It begins with a review of the methodology 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 applications are discussed. Specifically, new finite element modeling of antennas on doubly curved platforms, bandpass radomes and jet engine scattering are presented using a variety of mesh truncation schemes (boundary integral, absorbing boundary conditions and perfectly matched absorbers) are presented. Also, an entire section of the chapter is devoted to a reduced ordering method based on the Pade expansion. This reduce ordering method is used to obtain broadband responses from a few data points of the entire response. It is introduced in this chapter for the first time in connection with hybrid finite element systems and promises substantial computational savings. 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, 97

FIEM FOR ANTENNAS AND SCATTERING 2 microwave circuits and biomedical applications. For example, modern aircraft geometries contain large non-metallic 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 sub-wavelength size bodies of simple geometries, PDE methods and hybrid versions of these have shown a much greater promise for large scale simulations without placing restrictions on the geometry and material composition of the structure. This chapter 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 chapter 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. Section 6 describes a reduced order modeling approach for extrapolating broadband responses from a few data points. The method is based on the Asymptotic Evaluation Method(AWE) and is developed here for hybrid finite element systems. The final section of the paper gives some additional representative applications of the finite element 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 Vx[1 (V x E) -kr- E =-jkoZoJi-V x (1Mi) (1) 98

FEM FOR ANTENNAS AND SCATTERING 3 Conformal Patch Array / S/ / (b) So. Patches,- -— Z!.:!!, I -- Feed Microstrip line Absorber termination (a) (c) Figure 1: Computational domains for FEM analysis, (a) the various regions enclosed by So, (b) typical tetrahedral mesh, (c) computational regions for antenna analysis 99

FEM FOR ANTENNAS AND SCATTERING 4 in which E represents the total electric field, (er 7i4) 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 Mi 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 consider the weak form of (1), instead of solving it directly. This can be achieved by extremizing the pertinent functional [1] 1 -1 F(E) J ( x E) [- 1 (V x E)]-k(rE). E d + Ji E [ikoZoJi + V x (-1 ~Mi)] dv + jkoZo f E (H x n) ds O+Sf 1 r 1 d _E_ H+ (koZo { -(n x E). (n x E)ds + Z JJJ E EdV (2) 2 LJJSr R ZL JJJV 1) where R denotes the resistivity or impedance of the surface Sr, Sf stands for 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 nF. Also, V, is the source volume and Vi is the volume occupied by the load ZL, whose length and cross section are given by d and 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 resistive sheet condition [2] n x n x E = -Rn x (H+ - H-) where HA denote 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 required for the evaluation of the integral over Sr. 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 <R,T > / {(V xT) [.r (V x E)] - k.E. T} d + f/t T {j-koZ J + V x (-1. Mi)} dV + jkoZo f T. (H x n)dS So+Sf + koZo { (n x T) ( x E)dS + z-S JX T. 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. 100

FEM FOR ANTENNAS AND SCATTERING 5 The functional (2) 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 Ee = Eew% = [We]T {Ee} k=l (4) in which W6 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 Wk = 0 within the element) and normalized in such a manner so that the expansion coefficients Ee represent the average field value across the kth edge of the eth element. 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 101

FEM FOR ANTENNAS AND SCATTERING 6 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 funcNe tions, the discretization of (2) or (3) proceeds by introducing the expansion E = E, with e=l Ee as given by (4) 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) =, e = 1,2,..., N; k =1,2,...,M (5) kEk 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 Ne N8 Ne N, E[Ae]{Ee} + E[BS]{Es} + {IKe} + ETis} = 0 (6) e=l s=1 e=l s=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) or (3); Ns is equal to the aggregate of the surface elements (quadrilaterals for hexahedra or triangles for tetrahedra and prisms) used for the tessellation of SO, Sf and Sr; the column {K6} results from the discretization of the source integral over Vs; [B8] is the matrix associated with the surface integrals over So + Sf in (2) or (3) with Ns aggregated surface elements; and finally {Cs} results from the discretization of the corresponding surface integral in (2) or (3) involving the external incident field. Basically, {Cs} provides the excitations for scattering problems whereas {K6} 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 [Ae] = {[DW ][r]-l[DWe]T - k2[W ][er][We]T} dV (7) with _ T - [DWe { Wz}T We(8) [a w - {w}T {We and Ve denotes the volume of the discrete element while the subscripts (x, y, z) in (8) 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 )e), the diagonal entries of [A6] are simply modified with the addition of 102

FEM FOR ANTENNAS AND SCATTERING 7 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 -A f + V x [,]-I Mdv }l V (9) Ke}J [e]{ kozo Jpy }+ ([ is ) The specification of the matrix [Bs] and the entries of the corresponding excitation column {Cs} due to the incident field 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 (6), a process referred to as matrix assembly, the resulting matrix system will be of the general form [A0] [0] [0] f {EV} I b} (10) [ {E} [ [0] [0 ] {E} } { {b } } ( In this, {EV} denotes the field unknowns within the volume enclosed by So + Sf whereas {EB} represents the corresponding unknowns on the boundaries So and/or Sf. The excitation column {bV} results from the assembly of {Ke} and similarly {bB} is associated with {CS}. 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 O(N) and scalability can be attained on multiprocessor platforms. The boundary matrix [6] in (10) 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 103

FEM FOR ANTENNAS AND SCATTERING 8 0 20 40 60 80 IM:1:. I:.:1. *2: U. l2. 'I System 0 1oo00 -1. l:. I::., *::U OU 0:0 l:..:U 120L 0 I I I I I 20 40 60 columrn 80 100 120 x14 x 1 04 n U I m Mr E -. -.- -... - 0.51[ 1 1.5 2 2.5 wit -I dL-. -r J ip 54. '[ Jil I1*'i d IL' 2 *: 9 * miL'" Pr — j C 0 0.5 1 1.5 nz = 469151 2 2.5 3 X 104 Figure 3: Examples of matrix systems generated by the finite element-boundary integral (top figure) and finite element-ABC methods (lower figure) 104

FEM FOR ANTENNAS AND SCATTERING 9 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 [9] 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. 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 nx H(r) = x Hnc(r) n x Jj {[n'x H(r)] x V'Go(r, r') + jkoYo[' E(r)]Go(r, r') + j '. [ X E(r')]VG,(r, r') dS' (11) where r and r' denote the observation and integration points, respectively, and G,(r,r') = exp(-jkolr - r'l/(47rlr - r') 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]. 105

FENM FOR ANTENNAS AND SCATTERING 10 The discretization of (11) follows by introducing the expansion (over the sth surface element) n=# of edges E(r) - E S;(R) = [S]T {E*} (12) 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, n x Sk are simply equal to the traditional basis functions given by Rao et. al. [31]. However, as pointed out in [32], S' can be chosen independently of we provided care is exercised when (11) is substituted/combined with (2) or (3). Regardless, when (11) is substituted into (2) or (3) after discretization we will get the partly sparse, partly full system [33] given in (10) 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 1(b)), the integral equation (11) can be simplified substantially [36]. Specifically the integral relation between E and H on the aperture reduces to n x H = n x Hg~ + n x G(r, r') [E(r) x ni dS (13) 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,.2k 1 ) e- 1 jk-lr-'l G =-+-Jo (I+ _2VV) (14) \ k2 )47rlr- r'| In this, I is the unit dyad and the factor of 2 is due to image theory. Also, H9~ 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, H~c 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 137], [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 [9] {EB} within the iterative solver. For example, if the symmetric biconjugate gradient (BCG) method [45] is employed to solve (10) 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 [9] is Toeplitz or not depends on the shape of the surface SO and on the uniformity of 106

FEM FOR ANTENNAS AND SCATTERING 11 the mesh. It can be shown that [G] will be Toeplitz for planar and cylindrical surfaces provided the surface grid is uniform. The Toeplitz form of [S] 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 vorkstation. 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 4) 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]. I cuula >atds 1i:. 1 2 5 Figure 4: Overlay of a uniform grid on an unstructured mesh for implementation of the FFT to perform the matrix-vector products. 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 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 (N,2) down to 0 (NAI5) and reductions down to O (NA.33) have also been proposed. The main idea of the FMM [49] is to subdivide the surface So into groups, each containing approximately Mse, Nse 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 (N2I/Mse,) operations. The sum of these two operation counts yields a total operation count 107

FEM FOR ANTENNAS AND SCATTERING 12 of 0 (N5) for the boundary integral matrix when Mse V ANse-, 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 and 0.35A deep (all metallic) 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) RMS error (dB) Groove Width FE-BI (CG) FE-FMM Groove Width FE-FMM 25A 1072 277 25A 1.12 35A 2102 458 35A 1.2 50A 4291 784 50A 1.36 Table 1: CPU Times, storage requirement and error and filled with a dielectric having Or = 4 and /r = 1. By using a sampling rate of about 15 edges per wavelength, the system unknowns were Nee = 5256 and Nse = 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 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 (NA5) or even down to O (Nsl33). 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) 108

FEM FOR ANTENNAS AND SCATTERING 13 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 i) == a' Ett + - 3 [V x n(V x Escat),] + 7 Vt(V Ecat), (15) where a = alltlt + a2t2t2 P3 = (ltl + t2t2):/(jk0 + 3Km -g -k2K1)ol + 3Km a = ko + 3Km - - -2l)tl + (jko+ - 3-2/2)t2t2 jkO K"m j(ko Km a1 rT1[4n- Kg + D( jko - Kl) + 1 + KmAK] a2 = 7[4K - Kg + D(jko - K2) + 2- KmAK] 1 r/ D D - M2', AK = K1 - K,2 K 9 D = 2jko+5-Km-. tKm 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 K1, K2 are the corresponding curvatures, Km = (K1 + K2)/2 is the mean curvature and Kg = K12 is the Gaussian curvature. For K1 = K2 (spherical termination boundary), (15) 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 (15) into the surface integrals in (2) or 3), separating the scattered field from incident field, and making use of vector differential and integral identities [59], we have sascat ) A jkoZo ]Esat (H x n) dS = [(Ecat)2 + (Eat)2]dS + Jj r(V x Escat)2dS so - Jj (v. E t )[V (. Est)t]dS. (16) Next, expansion (4) or (12) is introduced and (16) is differentiated as in (5) to obtain the finite element equations (6). Since the ABC is a local condition, the sparsity of the finite element matrix is preserved and the resulting system is again given by (10) with [9] 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. 109

FEM FOR ANTENNAS AND SCATTERING 14 Metal Patch r 1-j2.7 Figure 5: Illustration of an antenna terminated by an isotropic absorbing material 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 5. 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) or (3). 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 5 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 6(a) displays a sectoral microstrip patch Metallic m Plandte fr xtens = t0.114cm etlatwih 3.1a 3. 33 t e f x, (a) (b) Figure 6: Sectoral microstrip patch on cone printed on a conical surface and Figure 6(b) shows the antenna resonance behavior for different cone angles. We clearly observe that the resonance frequency drops with the cone angle for the 110

FEM FOR ANTENNAS AND SCATTERING 15 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 6, 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 layer (PML) [20]. The intrinsic parameters of the PML are = a - j3 0 0 r =- 0 -j/3 0 a 0 0 1 0 a-jl3 J where a = (Z/Zo)2, Z is the intrinsic impedance of the medium being terminated by the absorber, whereas a, 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 7 [62] where it is shown that the anisotropic absorber retains its low reflectivity at; oblique incidences except near grazing ( = 0~). The performance shown in Figure 7 is based on a purely theoretical analysis but as shown in Figure 8, this performance can be realized in numerical simulations with careful choices of At and the sampling rate N within the absorber. It has been found [63] that the reflectivity curve shown in Figure 8 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/AX (A\ = Ao/ cos i) is a scalable parameter and that given the desired reflection coefficient [R], the formulas [62] =-0.01061RI + 0.0433 (17) x N = 0.147exp [7.353-1 (18) can be used to choose 3t/Ax and the minimum number of discrete samples to achieve the desired absorption. In these, IRI is specified in dB and Ao denotes the free space wavelength. These formulas are in agreement with the actual numerical results in Figure 8 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 (17)-(18) are given by Legault [62] for one such tapered anisotropic absorber. 111

FEM FOR ANTENNAS AND SCATTERING 16 y 0 x m cc: -10 -20 -30 -40 -50 -60 -70 Isotropic r = r =: b -80 - 0.0 0.2 0.4 0.6 sin 4 0.8 1.0 1.2 Anisotropic ~r= g r = b I - (b- i/b) x x (a) (b) Figure 7: Reflectivity of a A/4 thick planar metal-backed perfectly matched absorber (/3 = 1 - j2) as a function of incidence angle. (a) Geometry, (b) plane wave reflection coefficient vs. angle -25 "0 f 1.05 cm 0.21 cm t -. Uniaxial absorber Microstrip 3.095 cm - (a) (b) Figure 8: Reflection coefficient of the PML for terminating a microstrip line as extracted from a numerical implementation 112

FEM FOR ANTENNAS AND SCATTERING 17 4. Feed Modeling For scattering problems where the plane wave incidence is usually the 'source', the right-handside excitation has been explicitly given in (10) and will not be discussed further. However, for antenna parameter computations, the explicit form of {K6} in (9) will depend on the type of feeding scheme being employed. Below we discuss specific forms of {Ke} corresponding to different feeding choices. 4.1. SIMPLE PROBE FEED For thin substrates the coaxial cable feed may be simplified as a thin current filament of length I carrying an electric current I 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 Ill= jkoZoI I Wi(r), i = jl,j2,...jm (19) 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. 4.2. 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 V = E. d, where d is a vector whose magnitude is the gap width, and E the electric field across the gap, we have V that Ei = --—, where cosOi is equal to 1 if the ith edge is parallel to d. Numerically, this gap dcosOi voltage model can be realized by first setting the diagonal term Aii of the [A] matrix equal to unity and the off-diagonal terms Aj (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. 4.3. 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 J/ (E x H). dS (20) 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 9) [64] 113

FEM FOR ANTENNAS AND SCATTERING 18 E = -r, H= -, (21) r r with VC rc Io ho =- eo~-. (22) Zo 7r In these expressions, Crc is the relative permittivity inside the cable, E and H are the electric cavity patch patch r 2b 2a cavity-cable junction (a) (b) Figure 9: (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). and magnetic fields, respectively, measured at z = 0 and 1o 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 (22) 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 b AV = E:(b - a) = eoln-, i = Np(p = 1, 2,..., Nc) (23) a 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 OFc extend the mesh inside the cable or to employ a fictitious current probe. The derivation of OEZ and its incorporation into the system is 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. 114

FEM FOR ANTENNAS AND SCATTERING 19 4.4. OTHER FEED MODELS There are a few other commonly used feed models for simulating antennas and the associated network in the context of the finite element methods. In certain cases, the structures may contain detailed geometries which must be modeled with care to ensure the efficiency and accuracy of the simulation results. For instance, the configuration of an aperture coupled microstrip antenna may be efficiently modeled by applying the equi-potential continuity condition and the interested readers are referred to [65] for details. Also in modeling microwave circuits as antenna feed network, the excitation location along the network may have to be placed far from the antenna for probe models, and thus the modal excitation is an alternative to the probes described in section 4.1 and 4.2. This reduces the size of the computational domain without compromising accuracy. The modal field distribution is typically obtained using a simplified analysis model to truncate the 3D FEM domain. A 2D FEM code can be used as well for geometries having the same cross-section as the original feed network. In general, the antenna feed or feed network can be accurately modeled in the context of the FEM. Moreover, unlike the method of moments (MoM), the FEM provides the field distribution in the entire 3D computational space and this is particularly useful for visualization around the feed region and on the antenna. 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 0(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 CPU 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. 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: 115

FEM FOR ANTENNAS AND SCATTERING 20 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 5.1. 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 [66] 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 sections of long vectors and to achieve this for our type of matrices the jagged diagonal format [67] 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 M.Flops 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. 5.2. 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 116

FEM FOR ANTENNAS AND SCATTERING 21 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 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 10. This type of matrix reordering can deliver speed-ups as close to linear as possible. x 10 0.5 X 104 Figure 10: Reduced bandwidth of the FE-ABC system after application of the Gibbs-Pool-Stockmeyer profile reduction algorithm 117

FEM FOR ANTENNAS AND SCATTERING 22 5.3. 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 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 [68] 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. Reduced Order Modeling (ROM) of Frequency Responses Reduced Order Modeling (ROM) methods such as the Asymptotic Waveform Evaluation (AWE) have been successfully applied in VLSI and circuit analysis to approximate the transfer function associated with a given set of ports/variables in circuit networks [69, 70, 71]. The basic idea of the method is to develop an approximate transfer function of a given linear system from a limited set of spectral solutions. Typically, a Pade expansion of the transfer function is postulated whose coefficients are determined by matching the Pade representation to the available spectral solutions of the complete system. In the context of finite element systems, ROM can be employed to predict a frequency response of the antenna input impedance or the scattering cross section of a given structure from a few data points. That is, once a few frequency points have been obtained by solving the entire finite element system of equations, these solutions along with the associated matrices can be re-used to extrapolate a broadband response without a need to resolve the system at other frequency points. In this section we present the theoretical basis of ROM and demonstrate its validity for full wave simulations using the finite element method as the computational engine. In addition to using ROM for antenna impedance and radar cross section prediction as a function of frequency, the method can also be used to fill-in angular pattern data points, thus eliminating a need to recompute the entire solution at small angular intervals. Since typical partial differential ecluation (PDE) systems involve thousands of unknowns, ROM can indeed lead to dramatic reductions of CPU requirements in generating a response of antenna or scatterer without a need to resolve the system for the fields in the entire computational grid. However, it should be noted that the FEM matrix for the reference frequency points must be stored (in core or out of core) with the current development of ROM for frequency domain analysis and thus some trade-off between CPU and memory requirements is unavoidable. Nevertheless, in view of the large CPU saving afforded by ROM, this appears to be a very small price to pay. 118

FEM FOR ANTENNAS AND SCATTERING 23 6.1. THEORY OF REDUCED ORDER MODELING 6.1.1. FEM System When the functional (2) is discretized in connection with absorbing boundary conditions or artificial absorbers for mesh truncation, the resulting system can be decomposed into the form (Ao + kA1 + k2A2) {X} = {f} (24) where A- denote the usual square (sparse) matrices and k = 27r/A = w/c is the wavenumber of the medium. As usual, {f} is a column matrix describing the specific excitation. Clearly (24) can be solved using direct or iterative methods for a given value of the wavenumber as described earlier. Even though Ai is sparse, the solution of the system (24) is computationally intensive and must be repeated for each k to obtain a frequency response. Also, certain analyses and designs may require both temporal and frequency responses placing additional computational burdens and a repeated solution of (24) is not an efficient approach in generating these responses. An application of ROM to achieve an approximation to these responses is an attractive alternative. For these problems, the excitation column {f} is a linear function of the wavenumber and can therefore be stated as {I} = k {fi} (25) with {fi} being independent of frequency. This observation will be specifically used in our subsequent presentation. 6.1.2. Reduced Order Modeling To describe the basic idea of ROM in conjunction with the FEM, we begin by first expanding the solution {X} in a Taylor series about ka, the wavenumber at which the system solution is available. We have {X} == {Xa} + (k - ka) {X} + (k-ka)2 {X2} +... +(k - ka)1 {XI} + 0 {(k - ka)+1} (26) where {Xa} is the solution of (24) corresponding to the wavenumber ka. By introducing this expansion into (24) and equating equal powers of k in conjunction with (25), after some manipulations, we find that {Xa} = kaAo {fl} {X1} = Ao [{f} -A {Xa} - 2kaA2 {Xa}] {X2} = -A-1 [Al {X1} + A2({Xa} + 2ko {X1})] (27) {X1} = -Ao [A1 {X_i} + A2({X1_2} + 2ko {X-1})] with Ao = Ao + koA1 + kgoA2 (28) 119

FEM FOR ANTENNAS AND SCATTERING 24 Expressions (27) are referred to as the system moments whereas (28) is the system at the prescribed wavenumber (ka). Although an explicit inversion of A-1 may be needed as indicated in (27), this inversion is used repeatedly and can thus be stored out-of-core for the implementation of ROM. Also, given that for input impedance computations we are typically interested in the field value at one location of the computational domain, only a single entry of {X1(k)} need be considered, say (the pth entry) Xl(k). The above moments can then be reduced to scalar form and the expansions (27) become a scalar representation of Xf(k) about the corresponding solution at ka. To yield a more convergent expression, we can instead revert to a Pade expansion which is a conventional rational function. For transient analysis the Pade expansion can be cast by partial fraction decomposition [71, 74] into q Xzk k-= -k (29) Xq(k=XO +Zk -ka ki where Xqo is the limiting value as k tends to infinity. This is a qth order representation suitable for time/frequency domain transformation. As can be realized, the residues and poles (ri and ka + ki) in (29) correspond to those of the original physical system and play important roles in the accuracy of the approximation. As can be expected a higher order expansion with more zeros and poles can provide an improved approximation. The accuracy of ROM relies on the prediction of the dominant residues and poles located closest to ka in a complex plane. Its key advantage is that for many practical electromagnetic problems only a few poles and zeros are needed for a sufficiently accurate representation. For a hybrid finite element - boundary integral system, the implementation of ROM is more involved because the fully populated boundary integral sub-matrix of the system has a more complex dependence on frequency. In this case we may instead approximate the full sub-matrix with a spectral expansion of the exponential boundary integral kernel to facilitate the extraction of the system moments. This approach does increase the complexity in implementing ROM. However, ROM still remains far more efficient in terms of CPU requirements when compared to the conventional approach of repeating the system solution at each frequency. As an application of ROM to a full wave electromagnetic simulation, we consider the evaluation of the input impedance for a microstrip stub shielded in a metallic rectangular cavity as shown in figure 11. The stub's input impedance is a strong function of frequency from 1-3.2 GHz and this example is therefore a good demonstration of ROM's capability. The shielded cavity is 2.38cm x 6.00cm x 1.06cm in size and the microstrip stub resides on a 0.35cm thick substrate having a dielectric constant of 3.2. The stub is 0.79cm wide and A/2 long at 1.785 GHz and we note that the back wall of the cavity is terminated by a metal-backed artificial absorber having relative constants of 67 = (3.2, -3.2) and r = (1.0, -1.0). As a reference solution, the frequency response of the shielded stub was first computed from 1 to 3.2 GHz at 40 MHz intervals (50 points) using a full wave finite element solution. To demonstrate the efficacy and accuracy of ROM we chose a single input impedance solution at 1.78 GHz in conjunction with the 8th order ROM in (29) to approximate the system response. As seen in Figure 12, the 8th order ROM representation recovers the reference solution over the entire frequency band for both the real and reactive parts of the impedance. We conclude that the ROM technique is an extremely useful addition to electromagnetic 120

FEM FOR ANTENNAS AND SCATTERING 25 simulation codes and packages for computing wideband frequency responses, large sets of radar cross section(RCS) pattern signature, etc. using only a few samples of the system solution. 1.106 cm ~i -T 6 cm - 2.38 cm - Figure 11: Illustration of the shielded microstrip stub excited with a current probe. 1.5 2 2.5 frequency (GHz) Figure 12: Real and imaginary parts of input impedance computations based upon the 8th order ROM implementations using a single point expansion at 1.78 GHz. Solid lines: exact reference data; Dashed lines: 8th order ROM results. 7. Additional Applications We choose two 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 70) 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 121

FEM FOR ANTENNAS AND SCATTERING 26 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 [75]. 1A 4X 20. m co 10 4C. < b 10 -0 --10 --20 --30 - 4i I^,i '**n. I, I~ ^-I~0 BOR 44' * FEMATS -40- -,,.....,,,i,,,I,,,I, -90 -60 -30 0 30 60 9( Observation Angle O., deg. Figure 13: Backscatter pattern of a perfectly conducting conesphere for q$f and 90 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 [75]. Mesh termination surface is a rectangular box. Frequency Selective Surfaces (FSS): FSS structures [76] are arrays of tightly packed periodic elements which are typically sandwiched between dielectric layers. The periodic elements may be of printed form or slot configurations designed to resonate 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 122

FEM FOR ANTENNAS AND SCATTERING 27 slot array placed 60iils below top surface of 90mils layer slot array at 30mils below top of 90mils layer r = 4 /' 2 -90mils -1 m a3 -1 'o a 90Omils 15.85cm 2cm T Nabsorber, r 4. -jl0.4 6.45cm 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 through the FSS screen will retain its TEM character, the transmission coefficient of the FSS panel can be approximated as T() =101 og dB ae 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 T + 10 log [1 - (l -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 [77]. 123

FEM FOR ANTENNAS AND SCATTERING 28 8. Conclusion 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 requirements 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 [78], domain decomposition[79, 80], robustness, adaptive refinement[81], 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. Reduced order modeling techniques such as the AWVVE method are also very promising for reducing the computational requirements in generating broadband responses. Further development of AWE is certainly needed for its application in connection with hybrid finite element systems. 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. Reduced order modeling techniques such as the AWE method is another promising approach References [1] J. L. Volakis, A. Chatterjee, and L. C. Kempel, 'A review of the finite element method 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 124

FEM FOR ANTENNAS AND SCATTERING 29 [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 threedimensional 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. [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. 125

FEEM FOR ANTENNAS AND SCATTERING 30 [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. [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. 126

FEM FOR ANTENNAS AND SCATTERING 31 [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. ERoux, "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. [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. 127

FEM FOR ANTENNAS AND SCATTERING 32 [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. [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] J.L. Volakis, J. Gong and T. Ozdemir, "FEM Applications to Conformal Antennas, FINITE ELEMENT METHOD SOFTWARE IN MICROWAVE APPLICATIONS," Edited by Tatsuo Itoh, Giuseppe Pelosi and Peter Silvester, (Wiley, 1996) [66] 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. [67] 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. [68] 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. [69] S. Kumashiro, R. Rohrer and A. Strojwas, "Asymptotic waveform evaluation for transient analysis of 3-D interconnect structures," IEEE Trans. Computer-Aided Design of Integrated Circuits and Systems, Vol. 12, No. 7, pp 988-996, Jul, 1993 [70] L. Pillage and R. Rohrer, "AWE: Asymptotic Waveform Estimation," Technical Report, SRC-CMU Research Center for Computer-Aided Design, Carnegie-Mellon University, 1988 128

FEM FOR ANTENNAS AND SCATTERING 33 [71] E. Chiprout and M. Nakhla, "Asymptotic waveform evaluation and moment matching for interconnect analysis," Norwell, Kluwer Acad. Pubs, 1994 [72] J. L. Volakis, A Chatterjee and J. Gong, " A class of hybrid finite element methods for electromagnetics: a review," J. Electromagn. Waves Applications, Vol. 8, No. 9/10, pp. 1095-1124, 1994 [73] J. Gong, J.L. Volakis, A.C. Woo and H. G. Wang, "A Hybrid Finite Element-Boundary Integral Method for the Analysis of Cavity-backed Antennas of Arbitrary Shape," IEEE Trans Antenna and Propagat., Vol. 42, No. 9, pp 1233-1242, 1994 [74] Joseph Lehner, "Partial fraction decompositions and expansions of zero," Trans. Amer. Math. Soc., Vol. 87, pp 130-143, 1958 [75] 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 [76] E.L. Pelton and B.A. Munk, "Scattering from periodic arrays of crossed dipoles," IEEE Trans. Antennas Propagat., vol. AP-27, pp. 323-330. [77] H. Wang, Personal communication, China Lake, CA, 1995 [78] 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 [79] T. Cwik,"Parallel decomposition methods for the solution of electromagnetic scattering problems," Electromagnetics, Vol. 42 pp. 343-357, 1992 [80] 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 [81] 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 129