035067-6-T HYBRID FINITE ELEMENT AND MOMENT METHOD SOFTWARE FOR THE SERAT ARRAY 5th Quarterly Report Sanders, A Lockheed Martin Co. 95 Canal Street NCA1-6268 P.O. Box 868 Nashua, NH 030601-0868 December 24, 1997 35067-6-T = RL-2486 1

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: U-M REPORT No.: CONTRACT START DATE: END DATE: DATE: SPONSOR: SPONSOR CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: Hybrid Finite Element Design Codes for the SERAT Array 5th Quarterly Report 035067-6-T October 1996 September 1998 Jan 5, 1997 (5th Quarterly Report) Roland Gilbert SANDERS, INC, A Lockheed Martin Co. MER 24-1583 PO Box 868 Nashua, NH 030601-0868 Phone: (603) 885-5861 Email: RGILBERT@mailgw.sanders.lockheed.com P.O. QP2047 John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave Ann Arbor, MI 48109-2122 Phone: (313) 764-0500 FAX: (313) 747-2106 volakis @umich.edu http://www-personal.engin.umich.edu/-volakis/ T. Eibert(UM), Y. Erdemli (UM), K. Sertel(UM), D. Jackson(UH), J. Volakis(UM) and D. Wilton(UH) CONTRIBUTORS TO THIS REPORT: 2

TABLE OF CONTENTS TABLE OF CONTENTS......................................................................................................................... 3 LIST OF FIGURE CAPTIONS IN SUMMARY SECTION.................................................................. 3 CHRONOLOGY of Events (Updated Every Quarter)..................................... 4 MEETINGS......................................................................................................................... UPDATED MILESTONE CHART (page 6).......................................................................... 6 Executive Summary and Project Status........................................ 7 Summary of Code Modules and Capabilities.......................................................................... 7 This Quarter's Accomplishments........................................ 8 FSS-BRICK Driver Description.................................................................................. 15 FSS-PRISM Driver -Update................................................................................................................. 16 FSS-PRISM speed improvements with the adaptive integral method (AIM)............................. 17 FSS-PRISM Update For Non-Commensurate SERAT Arrays.............................. 21 Curved SERAT(FSS-CURVE) Array Update (Sertel)................................. 25 LIST OF FIGURE CAPTIONS IN SUMMARY SECTION Figure 1. Validation of FSS-PRISM for a 5-layer non-commensurate FSS structure; details of the geom etry are given later................................................................................................................... 8 Figure 2. Illustration of a bandgap substrate layer modeled with the FSS PRISM (see later section for details)......................................................................9................................... 9 Figure 3. Calculations for a planar and a curved SERAT array of dipole elements.............................. 12 Figure 4. Illustration of the SERAT panel (Planar and Curved)............................ 14 Figure 5: General m ultilayered FSS structure...................................................................................... 15 LIST OF TABLES Table 1. CPU and Memory Reduction after Implementation of AIM into FSS-PRISM...................... 10 3

CHRONOLOGY of Events (Updated Every Quarter) * April 1996 Proposal Submission * July 1996 Answers to Proposal Questions * August 1996 Began Contract Negotiations * 20 Sept. 1996 Kickoff meeting at Ann Arbor (attended by Sanders, UM and UH) * October 1996 Contract Signed between U-M and Sanders in Mid October * October 1996 Subcontract to the Univ of Houston (formalized in early November) * 15 Nov. 1996 SERAT Review meeting (at Nashua) * 9 January 1997 Submitted First Quarterly Report Report Described Code Plan and Progress on the Moment Method FSS Code. Specifically, a new scheme was developed to accelerate the convergence of the periodic Green's function * 28 February 1997 Prepared viewgraphs on the project's progress review. Showed first validation results for the moment method FSS code with the new accelerated Green's function; showed results for a new algorithm to accelerate the boundary integral truncation of the planar and curved FSS hybrid FEM code using the Adaptive Integral Method(AIM) and CVSS, the new LU solver specialized to sparse matrices * 5 April 1997 Submission of Second Quarterly Report. Report included the first validation results for the stand alone small array FEM code (with dipole FSS elements and dipole antenna elements). A similar validation was done for the moment method FSS developed at Houston. The fast AIM algorithm was described for boundary truncation and the TRIANGLE surface mesher was introduced to generate the aperture mesh, subsequently grown down to the FSS. * 30 May 1997 Semi-annual review at the Univ. of Michigan (attended by all parties) Review covered progress up-to-date. At this meeting, emphasis was on the validation of the three codes which took 'shape and form' between March-May 1997 in accordance with the proposed schedule. Theory, validations and comparisons among the codes were presented. * 25 June 1997 SERAT review Meeting (Nashua) * 5 July 1997 Submission of Third Quarterly Report The major component of this report was the description and validation of the periodic hybrid FEM code, FSS-PRISM. Comparisons among the FSS-EIGER, FSS-BRICK and FSS-PRISM were given for the first time. Also, the geometry drivers for the FSS-BRICK and FSS-EIGER were given. 4

* 6 Oct 1997 * 24 Oct. 1997 * Dec 1997 Submission of Fourth Quarterly Report Delivered FSS-EIGER and FSS-BRICK, both with manuals and preprocessors. Primitives are used to specify antenna and FSS elements in each code. Report gives example calculations for non-commensurate FSS panels for the first time using the scaling approach implemented in FSS-PRISM. Another first is the hybridization of the fast multipole method with the finite element code PRISM as a step toward the curved FSS modeling. Review at the Univ. of Michigan Attended by U-M and Houston project people, Gibert, Pirrung and Asvestas. Presented status of FSS-BRICK, FSS-EIGER, FSS-PRISM and FSS-CURVE. Delivered manuals for FSS-EIGER and FSS-BRICK (and codes); Successes for non-commensurate array modeling were presented (5 layer example); Update on FSS_EIGER for non-commensurate was given; Initial implementation of finite curved array code was presented with the fast multipole method for mesh truncation. 5th Quarterly Report The major highlights of this progress report are * Implementation of the fast integral method into FSS-PRISM, making prism a practical analysis code, even when the sampling requirements are very high * Completion of FSS/Antenna geometry Driver for FSS-PRISM * First implementation of curved arrays with the FMM for mesh truncation. * Additional validations for non-commensurate FSS MEETINGS * A semi-annual review meeting was held on Oct. 24, 1997. * Detailed viewgraphs was handed out along with manuals for FSS-EIGER and FSS-BRICK. 5

UPDATED MILESTONE CHART (page 6) 6

Executive Summary and Project Status In our previous quarterly report we gave an extensive description of the project status which was followed by the semi —annual project review held on 24 October 1997. Our last report along with the viewgraph slides and code manuals given at the review provide a clear picture of the project's status. Basically, as of this moment we have four operational codes whose development is either on schedule or ahead of schedule (refer to Milestone chart given on page 6). Summary of Code Modules and Capabilities CODE FSS-BRICK FSS-EIGER DESCRIPTION AND STATUS Small (finite) Array Code * Finite element-Boundary Integral (FEBI) code using brick elements * Completed and Validated * Geometry Driver described in Report 035067-5-T Moment Method code * Uses multilayered Green's function and triangular boundary elements (RaoWilton-Glisson formulation) * Uses Ewald acceleration for periodic Green's function * Validated (slot and dipole elements) and delivered for commensurate FSS * Geometry Driver Manual delivered in October 1997 * Non-commensurate FSS modeling still in progress Finite Element-Boundary Integral Code * Combines flexibility of finite elements for volume and boundary element for robust mesh truncation * Uses Ewald acceleration for periodic free space Green's function * Incorporates Adaptive Integral Method for fast Boundary element evaluation * Uses layer de-coupling to handle noncommensurate FSS. * Validated for commensurate and noncommensurate FSS and slot and printed antenna elements * Geometry Driver is now available completed for commensurate and noncommensurate FSS FSS-PRISM 7

FSS-CURVE Curved SERAT array code * Based on the non-periodic version of PRISM * Incorporates the fast multipole method (FMM) for fast non-planar boundary integral evaluation * Tested for planar large array simulations * Extended to non-planar simulations * Under testing for curved array simulations * FSS Geometry driver to be developed in the 6th quarter. This Quarter's Accomplishments As of this time we continue to be ahead of schedule in developing FSS-PRISM and FSS-CURVE as well as their associated geometry drivers and validations. This quarter's summary of accomplishments are as follows (refer to the attached document for more technical details): 1. After examining several approaches (see previous quarterly report), we finally hpim 9 -9! 9 9 9 27 r glpim - 9 18 I ---- 9 ------- is 18 18 12 E=231-j.... L -I rl 6 14 14 14 c 8 _ I U. I 0.8 - - - - 0.7 — t- - ---- -- I II 0.6 t tc t -;9l t - ---- -— e- - 0.5 ----------—; ----. --- —---- I I iI I I 0.4 ------— i --- —7 —l ---E ---J ---0.4 i i L I I I 0.3 Aroudaki -— t ---- PRISM E=2.3 -j0.12 I -...PRISM Er=2.31-j.l5 - - --- 0.1 --— I ---- --- -- 0 1000 2000 3000 4000;5000 6000 7000 8000 Frequencyin GHz: Figure 1. Validation of FSS-PRISM for a 5-layer non-commensurate FSS structure; details of the geometry are given later. resorted to the method of layer de-coupling to handle non-commensurate FSS. This approach was successfully implemented and tested using reference data for a 5-layer FSS structure as depicted for example in Figure 1. Clearly, this is a major highlight of the project since it readily permits the 8

implementation of non-commensurate FSS structures without any geometrical or material restrictions and with good control on the accuracy of the implementation. Any antenna radiator (slot or printed) or FSS element can be accommodated provided the geometry and mesh is appropriate discretized and specified using the supplied geometry Driver or an external mesher. 2. Substantial speed-up was achieved with the implementation of the adaptive integral method (AIM) for calculating the boundary integral portion of the matrix-vector products. This implementation was considered as optional in the proposal due to its recent introduction and uncertain level of success. Also, our intend was to only include a fast integral method with curved array only. Nevertheless, its implementation into our most capable general code FSS-PRISM will allow to deliver a software package that incorporates the best and fastest computational approaches for SERAT modeling. We have already performed (see later section) several CPU and memory comparisons to evaluate the savings when AIM is invoked in the iterative solver. More specifically, for a 32x32 elements periodic cell (see Table 1) * memory was reduced by a factor of 10 (an order of magnitude) * CPU was reduced by a factor of 2.5 allowing the calculation of each frequency point within 50 minute on a conventional workstation for a system which included over 40,000 unknowns, 6000 of which were part of the dense boundary integral subsystem For a large 44x44 elements periodic cell grid, use of AIM resulted in solving a system over 70,000 unknowns (with 12,000 as part of the dense boundary integral subsystem) using only 70MB. Without using AIM, the corresponding system would require 850MB of memory and would therefore be beyond the capability of desktop workstations. The CPU time was only 2 hours for this very large problem demonstrating that the developed FSS-PRISM code is also a very practical code for the analysis (and possibly design) of very complex multilayered FSS and antenna geometries. 3. Validated and extended the capability of FSS-PRISM for modeling bandgap antennas such as those shown in Figure 2. Figure 2. Illustration of a bandgap substrate layer modeled with the FSS PRISM (see later section for details). 9

Table 1. CPU and Memory Reduction after Implementation of AIM into FSS-PRISM CPU Times 32x32 Grid Run Data DOFs = 35,432(vol) + 2x3136(Surface) Memory: with AIM: 2183k elements (26MB) w/o AIM: 20,000k elements (240MB) CPU(IBM RS6000): FILL= 38 sec(FEM)+ 1930sec(BI)+95 sec(BI) BCG Solver= 415 Iterations (1081 sec) Total CPU=3049sec (~50min) FILL=38sec(FEM)+7330sec(BI) Conventional BCG Solver=364 Iter. (1446sec) Total CPU=8814sec(~2.5 hours) AIM reduced memory by a factor of 9 CPU by a factor of 3 10

4. Completed the Geometry Driver of FSS-PRISM for commensurate and non-commensurate elements. With this, the periodic planar codes were completed on schedule (end of fifth quarter). A manual is currently being written which will parallel the manual written delivered during this quarter for the FSS-BRICK (see UM Radiation Lab report 035067-5-T) 5. Developed a first version of the curved finite periodic array code FSS-CURVE which incorporates the fast multipole method (FMM) to speed-up boundary integral calculations. The FMM was employed rather than AIM because FMM is easier to implement for non-planar geometries. Nevertheless, AIM is much more attractive for planar geometries and is the reason for using it in connection with the planar code FSS-PRISM. At this stage we have already extended the implementation to curved arrays as shown in Figure 3. Again the fast integral method, FMM, plays a crucial role in reducing the CPU by a factor of 5 or more. We therefore anticipate that the resulting curved SERAT code will also be a practical useful tool. At this point, the development of the curved SERAT code is ahead of schedule. Over the next quarter, our goal 'will be to * validate and examine capabilities of the curved SERAT code for several practical geometries, including the five layer FSS modeled by planar SERAT code FSS-PRISM * incorporate the PRISM geometry Driver into the curved SERAT array code. 11

80 cm, l I I II I I 30 - 20... 10 - 40,.- *. *............ 40 1200 f000 *-r SOO IS 0 E 600 -, o 10 I I I - 10 I I -2 - 10 - 40 -so I i i i I I II I i................... -........................ ~..i...................". ".. —..-... - I" ':u.:,a..............11... 0 1000 200 GO' 000 4000 S000 6000 H 4 rm br or B I u kn own:; -$0 -O6 -40 -20 0 20 40 60 so i Figure 3. Calculations for a planar and a curved SERAT array of dipole elements. Upper left: Planar array geometry Upper right: Curved array geometry Lower left: CPU comparisons with and without FMM Lower right: C-omparison of radiation patterns of the shown curved and planar array 12

Project Goals The goal of the SERAT project at the University of Michigan (with subcontract to Univ. Of Houston) is to develop a suite of software for the analysis of strip and slot dipoles on multilayered substrates backed by a frequency selective surface. The dipoles are equipped with photonic switches permitting variable electrical dipole lengths for broadband performance and the FSS is suitably designed to simulate a variable substrate thickness for optimal operation. A general view of the geometry is given in Figure 1. The UM/UH team proposed to construct a code which combines various computational modules interfaced with appropriate pre-processors and post-processors. The computational modules include: * Stand-alone moment method simulation of the FSS with up to 10 layers with commensurate and non-commensurate periodicities. * Simple moment method simulation of the antenna elements on the FSS panels * Hybrid FEM simulation modules for small arrays, planar periodic arrays and curved arrays on FSS panels. Various options for modeling the FSS and for mesh truncation were proposed to provide a compromise between speed and accuracy. These are outlined in the proposal and summarized in the included milestone chart (repeated from the proposal). As called for in the milestone chart, we are proceeding in accordance with the schedule in our proposal. Specifically: 1. The commensurate moment method FSS and moment method periodic array code was completed and tested by Houston and delivered to UM. This code is referred to as FSS-EIGER and will be delivered in October 1997 as scheduled. A formulation for the non-commensurate version of the code has been developed and is under implementation at this time. 2. A geometry driver has been developed for FSS-EIGER to handle the dipole, crossed-dipole, slot and crossed-slot elements. 3. UM completed ahead of schedule the small array FSS code. This code will be delivered in October 1997 along with a geometry Driver for antenna and FSS elements 4. A first version of the centerpiece of this Code-Suite, the hybrid FEM code FSS-PRISM, is completed for commensurate arrays and combines the powerful FEM for material modeling with the robust boundary integral method. A Driver is now under development and a noncommensurate version of the FSS-PRISM is being tested in accordance to the schedule. 5. Several code validations have been given for slot/dipole periodic arrays on FSS structures and example data for non-commensurate panels will be shown in October 1997 6. We have already began the incorporation of fast algorithms into the SERAT hybrid code as scheduled for an accurate model of the curved SERAT panels 13

-Manct leAl S. I AI 0 F'l................................................................................................................................................................................................................................................................................................................................................................................................................................................... 1-..........................w.............................................................................................. -............................................................................-............................................ 11...................................I - I...............-.......................,,, ''.,,..........................................................I.....................................................................................................................................:................................................................................................................................................................................................................................................................................................................................................................................................ I.................................................. -,............................................................................................................................................................................................... I........................................................................ I...................................................................................... -.............................................................................................. I........................................................................................................................................................................................ I................................................................................. Figure 4. Illustration of the SERAT panel (Planar and Curved) 14

FSS-BRICK Driver Description Refer FSS-BRICK Manual (UM Radiation Lab report 035067-5-T) FSS-BRICK is used to analyze electromagnetic scattering and radiation characteristics of planar finite antenna arrays and/or frequency selective surfaces (FSS) residing in a three dimensional metal-backed cavity. The general multilayered FSS geometry is shown in Figure 1. /Vf /Vf A/ /V A/ /Vf // L/ ZL./ Y IV / -.. r -------- -IV-e s ---........ -5 vv Figure 5: General multilayered FSS structure. The driver is an I/O tool which enables the users to enter the necessary geometry and excitation information for antenna or FSS analysis. It is also employed to display the computed gain and input impedance or transmission coefficient characteristics of the corresponding structure. Using the geometrical input data, the FSS-BRICK driver generates the antenna mesh consisting of rectangular brick elements without a need for an external meshing package. The details of the Driver features are described in FSSBRICK User's Manual (UM Radiation Lab report 035067-5-T). 15

FSS-PRISM Driver -Update The PRISM driver has an 110 operation similar to that of the BRICK driver. The Driver allows the interactive specification of resistive cards, horizontal and vertical loads, vertical short-circuit pins, horizontal and vertical probe feeds in a multilayered structure. The FSS elements can be either strip dipoles, slots and cross elements. The user can also specify FSS elements other than those primitives by bypassing the geometry Driver and making use of a commercial meshing tool. Commensurate as well as non-commensurate geometries can be handled by the PRISM driver. Using the input data, the Driver generates a surface mesh consisting of triangular elements at each interface. The volume mesh is then grown along the depth of the structure forming three-dimensional prismatic elements. The surface mesh for each FSS layer and a 3-D view of the whole geometry can be displayed using MatLab which is a platform independent: tool and reads the output file of the Driver/PRISM mesher. Each discrete element in the structure is assigned a different colors in the Matlab display. The Driver stores the necessary excitation and mesh information in data files to be used by FSS-PRISM. The development of FSS-PRISM User's Manual is now in progress, and the Driver features will be described in detail in an upcoming Manual. 16

Hybrid FE/BI Modeling of Commensurate and Non-Commensurate 3D Doubly Periodic Structures Thomas F. Eibert and John L. Volakis Radiation Laboratory, EECS Department The University of Michigan, Ann Arbor, MI 48109-2122 1 Introduction By applying appropriate periodicity conditions, the computational domain of infinite periodic structures can be reduced down to a single unit cell. Previous analyses have mainly been restricted to the application of Floquet's theorem to construct periodic Green's functions in the context of integral equation formulations. Fully three-dimensional modeling capability can be achieved by using hybrid finite element (FE) / boundary integral (BI) methods as presented in [1, 2]. Both of these methods are based on FE modeling employing tetrahedral meshes and a spectral domain Floquet mode expansion of the BI fields. In this paper we propose a hybrid FE/BI method employing distorted triangular prismatic elements for volume tesselation and a mixed potential integral equation (MPIE) formulation in spatial domain for mesh truncation. The triangular prismatic elements provide for full geometrical adaptability in the plane of the triangles and structured meshing along the depth of the cell. For the evaluation of the spatial periodic Green's function the Ewald acceleration technique is applied [3]. In general, an exact formulation of the periodic array problem is possible only, if the whole structure has a unique periodicity (commensurate case). However, if the periods change from layer to layer, referred to as non-commensurate periodicity, an exact analysis is also possible, if a super-period can be defined which is an integer multiple of all different layer periods. Although this approach does provide a solution option, in many cases it is an expensive alternative since the super-cell may include a large number of the different layer periods. In principle, arbitrarily accurate solutions for non-commensurate periodic lattices can be found by separately modeling each of the periodic layers and tracking the individual Floquet modes forward and backward from one layer to the other. This cascading procedure is also very time consuming prompting the introduction of approximate methods that can give good results with lower computational effort. Such an approach was proposed in [4] and [5] for planar integral equation methods and is based on decoupling the periods of the currents residing in different layers. In the case of the hybrid FE/BI method the situation is somewhat involved because the fields themselves are modeled within the layers and the meshes of neighboring layers must be connected to each other. Therefore, we propose an approach which is based on decoupling of the field periodicities in the single layers. By grouping together several unit cells of the individual layers, the accuracy of the method can be improved but with higher computational effort. 2 Formulation We consider the periodic structure illustrated in Fig. 1. For time harmonic electromagnetic fields (e-wt time convention) the pertinent finite element functional is

F(E, Ea = / [- (V x Ead) (V x E)- k E Ead + jkoZo Ead int dv V +jkoZo JEad (H x n) ds, (1) S where Ead is the solution of the adjoint field problem, Jint denotes an excitation current within the FE domain, S represents the planar surfaces on top and/or bottom of the FE domain, and n is the unit surface normal directed out of the FE domain. Also, ko and Zo are the free space wave number and wave impedance, respectively. Invoking the surface equivalence theorem, <i, - - ------ - /- --- - I /- r -, - - - - I- -I -, \ I, / I / I I I.,/ - Dxp Figure 1: Infinite periodic structure the magnetic field in the surface S of (1) can be expressed as ko 1 H =-2j O g(rr)(E x)ds+ V J(r, r)V' (E x f) ds + Hinc (2) where g(r, r') is the appropriate scalar Green's function and Hinc i is an incident wave in the presence of a metallic interface in S. Assuming an array periodic in the xy-plane with Dxp and Dyp as the periods in the x- and y-directions, respectively, the periodicty condition for E in a rectangular lattice is given by E(x + m Dxp, y + n Dyp) = E(x, y) e-jtxm DX e-J:lyn Dyp (3) where 3x = ko sin o cos cpo, /3y = ko sin90sin po (4) and (,9o, 9o) represent the scan direction of the array. Using the periodicity condition and the boundary integral representation (2) the computational domain can be reduced to a single unit cell of the array. In the implementation of the hybrid FE/BI method this requires the imposition of appropriate phase boundary conditions in the FE and in the BI portions of the method. Moreover, in (2) g(r, r') must be replaced by the Green's function for an infinite array of d-sources. In our implementation, the infinite series representation of the periodic Green's function is calculated using the Ewald transform [3] which subdivides the series into spatial and spectral domain components in a manner that achieves optimal convergence. In the case of non-commensurate arrays, the fields in the different layers are assumed to fulfill the periodicity conditions of that layer resulting in a decoupling of the different periodicities.

This is achieved by meshing the structure according to the geometric periods of each layer as illustrated in Fig. 2. For the implementation of the phase boundary condition within the FE mesh we must carefully deal with the staircased side walls. By periodic continuation of the different decoupled layers the unknowns for the edges on the horizontal boundaries separating the layers can be related to the corresponding edges inside the volume mesh to bound the FE domain (see Fig. 2). In general, for the boundary integrals the method requires that Green's functions with different periodicities for the top and bottom layers must be utilized. 1 jvl I I I '\I I [A r-l I Figure 2: FE/BI modeling of non-commensurate structures (gridded portion is the computational domain which comprises a single period from each of the layers) 3 Results Bandgap Lattices: As a first example let us consider the dielectric slab in Fig. 3 with embedded periodic material blocks. These lattices are often referred to as photonic bandgap materials. The diagram in Fig. 3 shows the reflection coefficient of a normally incident plane 1.0 FE/BI 0.8 - - - Yang et al. 1997 (vol. integr. eq.) E 0.6 C 0.4., 0.2 z 0.0.......... 0.0 4 6 8 10 12 14 16 Frequency (GHz) Figure 3: Specular reflection coefficient of a plane wave ('0 = 0, pO = 0) from a dielectric slab (Cr = 4) with planarly embedded periodic material blocks (Er = 10); Slab height: 0.2 cm, period: 2 x 2 cm, block side length: 1 x 1 cm. wave onto the slab. The reflection coefficient curve exhibits the typical resonances of photonic bandgap materials. Compared to calculations obtained by a volume integral equation method [6] the first resonance is slightly shifted to a lower frequency whereas the frequency shift for the second resonance is larger due to the decreasing sampling rate per wavelength. Calculations for TM-waves with oblique incidence resulted in a shift of the resonances to higher frequencies and this is in agreement with [6]. Low-Pass FSS: Fig. 4 gives the geometry and results for a 5-layer low-pass FSS filter with non-commensurate layer periodicities (periods as illustrated in the graph). The given curves represent FE/BI results for different models based on the proposed method for non

commensurate periodicities and compare well to calculations performed using a spectral domain integral equation method [5]. 1.0 - FE/BI model 1 -- FE/BI model 2 0.8 / 5 FE/BI super cell -' -y Aroudaki et al. 1995 o d \(integral equation) O 0.6 h/[pm g/Pm c/mP 0.04 9,um, 1 patch with period 18 sum, 1 patch with period 12 m. Model 2: 2x2 patches with period 9 Atm, 1 patch with period 18 sum, 2x2 patches with period 12 Am. Super-cell: 4x4 patches with period 9 Urn, 2x2 patches with period 18 Atm, 3x3 patches with period 12 Atm. Acknowledgments This work was supported in part by the U.S. Navy and Sanders, A Lockheed Martin Corp.. The first author acknowledges the Fellowship support from the "German Academic Exchange Service (DAAD)". References [1] E. W. Lucas, T. W. Fontana, A 3-D Hybrid Finite Element/Boundary Element Method for the Unified Radiation and Scattering Analysis of General Infinite Periodic Arrays, IEEE Trans. AP, Vol. 43, No. 2, Feb. 1995, pp. 145-153. [2] D. T. McGrath, V. P. Pyati, Periodic Structure Analysis Using a Hybrid Finite Element Method, Radio Science, Vol. 31, No. 5, Sep/Oct. 1996, pp. 1173-1179. [3] D. R. Jackson, D. R. Wilton, Personal Communications, 1997. [4] J. C. Vardaxoglou, A. Hossainzadeh, A. Stylianou, Scattering from Two-Layer FSS with Dissimilar Lattice Geometries, lEE Proc. H, Vol. 140, No. 1, Feb. 1993, pp. 59-61. [5] H. Aroudaki, V. Hansen, H.-P. Gemund, E. Kreysa, Analysis of Low-Pass Filters Consisting of Multiple Stacked FSS's of Different Periodicities with Applications in the Submillimeter Radioastronomy, IEEE Trans. AP, Vol. 43, No. 12, Dec. 1995, pp. 1486-1491. [6] H.-Y. D. Yang, IR. Diaz, N. G. Alexopoulos, Reflection and Transmission of Waves from Multilayer Structures with Planar Implanted Periodic Material Blocks, J. of the Opt. Soc. of Am., B, Oct. 97.

Adaptive Integral Method for Hybrid FE/BI Modeling of 3D Doubly Periodic Structures Thomas F. Eibert and John L. Volakis Radiation Laboratory, EECS Department The University of Michigan, Ann Arbor, MI 48109-2122 1 Introduction The application of hybrid finite element (FE) / boundary integral (BI) methods to infinite periodic structures (antennas or frequency selective surfaces) [1, 2] is very attractive because it provides full 3D modeling flexibility. Basically, the finite element method is used to model the dielectric section of the unit cell, whereas the boundary integral provides for a rigorous mesh truncation at the upper and/or lower surfaces of the unit cell. Most practical problems can be analyzed by using planar BI surfaces and the corresponding half-space periodic Green's functions. However, for large unit cell apertures, the resulting BI matrix still becomes the most CPU intensive portion of the method. To alleviate the CPU and memory bottleneck, in this paper we present an acceleration and memory reduction scheme for the BI portion of a hybrid FE/BI method. The approach is based on the adaptive integral method (AIM) [3] and is adapted here to periodic structures. For the given problem, AIM results in low O(ns) storage and O(ns logns) CPU time requirements for the execution of the matrix vector products in the applied iterative solver (ns = number of surface unknowns). The paper focuses especially on AIM issues related to infinite periodic structures. Also, we present CPU time and storage comparisons with the more conventional implementation of the FE/BI method. 2 Formulation The conventional implementation of the hybrid FE/BI method for an infinite periodic structure (ei't time convention) results in a linear algebraic system of the form Aint Across Across F int F 0 0 1 F int f i A] Altop l~ Eto Z 0 Mt Lcross bound f JEbound 1 [ n Zo I ) rboundfb _ J. (1) A2,top Atop | \ btop + | ~ Ztop 0 | top bound ( * (1) cross bound Wbound L 2,bot ~ Abot. { Ebot ) L ~ ~ Zbot ' I Ebot ) The matrices A are sparse (20 to 40 non-zero elements per row) and are associated with the FE portion of the volume. The BI matrices Z are fully populated and are associated with the boundary edges on the top and bottom boundaries of the discretized periodic cell. In our implementation we employ triangular prismatic elements in E field formulation to generate the volume mesh giving rise to triangular surface meshes with Rao-Wilton-Glisson basis functions for the magnetic currents on the planar BI surfaces. For arrays with arbitrary scan angles, all matrices may be non —symmetric and this is due to the phase boundary conditions and the periodic Green's function. For the solution of the system the biconjugate gradient (BiCG) or generalized minimal residual (GMRES) solver is used, both of which rely on an efficient evaluation of the matrix-vector products associated with (1). The computational complexity per matrix-vector product of the total A matrix is of 0(nv) where nV denotes the number of volume edges. However, the complexity in executing the matrix-vector products [Z]{x} is of order 0(n2) and storage requirements are also of the same order. For a more efficient implementation of the method, it is therefore crucial to reduce the complexity of the latter

matrix-vector product.. Here we propose to accomplish this using AIM. The basic idea of AIM is to split the Z matrices as [Z] = [Znear] + [zfar], (2) where [Znear] contains the elements of [Z] that are near the self-cell and correspondingly [Zfar] contains the remaining "far-zone" elements of [Z]. AIM reduces the CPU time and memory of the iterative solver by exploiting the convolutional properties of the Green's function for the evaluation of the matrix-vector products associated with the mostly full matrix [Zfar] [3]. That is, the far-zone matrices are not explicitly generated and the matrix-vector products are performed in the discrete Fourier domain (DFT) utilizing appropriate 2D FFT algorithms [5]. However, the convolutional properties can only be exploited on a regular grid and therefore, the basis functions on the original and possibly irregular triangular mesh must first be mapped onto a regular grid. To ldo so, we introduce auxiliary basis functions on a regular rectangular grid which is coincident with the original triangular grid [6]. The auxiliary expansion for the magnetic surface currents is given by M N MUx = E E [Mm + M mn] 6(x -xo- mA) 6(y - yo - nAy), (3) m=l n=l where (xo, yo) is the lower left corner of the regular grid and Ax, Ay are the sample distances in the x and y directions, respectively. Since we work with a mixed potential integral equation (MPIE) formulation for the BI implementation, we introduce an additional expansion M N V. Ma = E qn (x - - mAx) (y - o - ny) (4) m=1 n=l1 for the divergence of the surface currents. In mapping the original basis functions onto the regular grid we relate each basis function to the auxiliary J-basis functions on the 9 closest grid points of the regular grid. This is done by equating a sufficiently large number of the corresponding moments of basis functions from the two coincident grids [3]. All the moments can be evaluated analytically and for each basis function of the original mesh three linear algebraic systems of order 9 must be evaluated. This leads to the construction of the sparse mapping matrices [A]x, [A]y, and [A]q with ns rows and 9 columns. So, the Z matrices for the rectangular grid can be expressed as [Z]AI= M = [h [A][ [G] [A] + [A ]q[G] [A, (5) which can be interpreted as the result of a numerical integration of the original Z matrices [3]. [G] is a Toeplitz matrix whose elements are the periodic Green's function evaluated on the grid points of the regular grid. For the implementation of AIM it is assumed that [Z]tAota is sufficiently accurate for [Zfar] but not so for the near-zone elements. Therefore, we decompose the AIM matrices as AIM = [Z]AIM + [Z]AIM (6) and when this is combined with (2), we can rewrite the original Z matrices as [z]apprx = ([Z]near - [z]near) + [total ( AIM] + [ZAIM (7) With this representation of [Z] the near-zone elements are evaluated without compromise in accuracy. However, since the great majority of [Z]approx includes the Toeplitz kernel [G], the

associated matrix-vector products can be performed using only 0{(ns) memory and 0({ns log ns) CPU time. In the final numerical implementation, a near-zone threshold is defined so that [Z]approx is a sufficiently accurate representation of [Z]. This threshold is mostly affected by the quasi-static singularities of the Green's function and can be reduced if the density of the regular grid is increased. In the case of an infinite periodic Green's function, we must also account for the singularities associated with the elements which neighbor the boundaries of the unit cell since the latter are adjacent to the next periodic cell. Therefore, test basis functions which supports are close to these singularties have to be considered to be in the near-zone of the source element. For the calculation of the matrix-vector products in the iterative solver [Z]fa is not computed explicitly. After mapping the actual source distributions onto the regular grid through the A matrices, the pertinent matrix-vector products are performed in the DFT domain using a FFT algorithm for the corresponding transformation. After transforming the results back into spatial domain, the fields on the original mesh are then obtained by reverse mapping between the auxiliary unknowns and the original grid unknowns. However, for the computation of [Z]nea, it is advantageous to collect the contributions of the regular grid in spatial domain before the Toeplitz Green's function is transformed into the DFT domain. 3 Results To show the speed improvements and the storage reduction due to AIM, we consider here a relatively simply strip dipole frequency selective surface (FSS) as shown in Fig. 1 (see also [4]). We calculated the reflection coefficient for three different FE meshes as described in Table 1 using triangular surface meshes for discretizing the surface of the unit cell and right-angled prisms for modeling the volume. As described above, the BI was applied on top and bottom of the volume meshes resulting into two fully populated identical BI matrices for the given problem. Therefore, 1.0 -0 — FE/BI 20x20.2 0.8 -0-0- FE/BI 32x32 0. — 6 FE/BI 44x44 3 A MOM 9x18 06 MOM 3x8 u 0.6 c o l ocm 0.4II. 0.2 0.0 12 14 16 18 20 22 24 Frequency (GHz) Figure 1: Resonance curves for strip dipole FSS only one of the matrices was computed. Table 1 shows the model parameters and the size of the system matrices with and without AIM for the BI matrix computations. The given numbers of matrix elements, which determine the storage requirements, in columns 5 and 6 of this table comprise both the top and bottom BI matrices. In the case of AIM, only the near-elements of the BI matrices have to be calculated explicitely. This results in memory savings of upto a factor of 15 for model 3. Moreover, the saving will increase further for surface meshes with a larger number of unknowns. The CPU times for the different models are summarized in Table 2. The

FE matrix fill time is, of course, the same for AIM and the conventional implementation of the hybrid method. For the investigated problem the largest CPU savings were obtained due to the reduced BI matrix fill time. In the solver, the necessary AIM overhead even resulted in a slightly higher CPU time for model 1, but for the larger systems, the lower complexity of AIM resulted in considerable CPU reduction (see CPU for model 2). For model 3 we only performed the calculations based on AIM since the storage requirements of more than 850 MByte without AIM are too large for desktop calculations (AIM reduced the storage down to about 70 MByte). The obtained resonance curves for the FSS structure in Fig. 1 show that the FE/BI results for lower resolution exhibit a small shift to a lower resonance frequency when compared to data calculated with a moment method code using multilayered media Green's function. The FE/BI curves for model 2 and 3 are identical implying that convergence was achieved with about 32 elements per linear wavelength (at resonance). However, we should note that our converged results still have a small frequency shift of about 1% between the FE/BI and the moment method results. Mod. Surf. mesh Surf. unk. Vol. unk. Matrix elements __ ________AIM Conv. 1 20x20 2*1240 7435 588 079 3 000 000 2 32x32 2*3136 35432 2 183 120 20 000 000 3 44x44 2*5896 67000 4 672 610 71 000 000 Table 1: Unknowns and storage for the different discretization models Mod. FE fill BIfill Solver Iterations Total l |AIM |Conv AIM Conv. AIM Conv. |~C AIM Conv. 1 3 736 1400 228 166 227 234 967 1 569 2 38; 1 930 7 330 1 081 1 446 415 364 3 049 8 814 3 70 33 821 5 073 676 8 964 Table 2: CPU times (in sec) for the different discretization models Acknowledgments The moment method code for planar periodic structures in multilayered media was provided by D. R. Jackson and D. R. Wilton from the University of Houston. This work was supported in part by the U.S. Navy and Sanders, A Lockheed Martin Corp.. The first author acknowledges the Fellowship support from the "German Academic Exchange Service (DAAD)". References [1] E. Lucas, T. W. Fontana, IEEE Trans. AP, Vol. 43, No. 2, pp. 145-153, Feb. 1995, [2] D. T. McGrath, V. P. Pyati, Radio Science, Vol. 31, No. 5, pp. 1173-1179, Sep/Oct. 1996. [3] E. Bleszynski, M. Bleszynski, T. Jaroszewicz, Radio Science, Vol. 31, No. 5, pp. 1225-1251, Sep./Oct. 1996. [4] R. Mittra, C. H. Chan, T. Cwik, Proc. IEEE, Vol. 76, No. 12, pp. 1593-1615, Dec. 1988. [5] J. L. Volakis and K. Barkeshli, Appl. of It. Meth. to Electromag. and Signal Proc., ed. T. Sarkar, Elsevier Pub. Co., pp. 159-240, 1991. [6] J. Gong, J. L. Volakis, A. Woo and H. Wang, IEEE Tians. AP, Vol. 42, No. 9, pp. 1233-1242, Sept. 1994.

FMM Solutions of Finite Element - Boundary Integral Implementations for Antenna Modeling Involving Doubly Curved Geometries K. Sertel and J. L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 1 Introduction Finite element - boundary integral formulations have been successfully applied to the analysis of cavity-backed antennas that possess fine details as is the case of with slot arrays, spirals, log-periodic and other broadband antennas. In those cases, the computational demand of the boundary-integral is large compared to that of the finite-element volume sub-system. Therefore, a fast integral method is attractive to speed-up the solution algorithm. For cavity-backed antennas recessed in an infinite ground plane, the boundaryintegral formulation is quite straightforward[l, 2, 3]. However, if the cavity resides in a curved platform periodic boundary conditions cannot be applied for array modeling. In this paper we present an approximate fast multipole method (FMM) solution to problems involving antennas recessed on doubly curved platforms. The approach is to use the halfspace Green's function as if the cavity resides in a flat ground plane. Clearly the limit of this approximation is applicable for geometries which are not highly curved. The region of validity of the approximation is investigated. 2 FE-BI Formulation for Cavity-Backed Antennas Consider a cavity-backed antenna recessed in a ground plane as depicted in Figure 1. This class of configurations have been modeled successfully using the finite element method [1, 2, 3]. The most rigorous of the implementations is to employ the finite element method to model the interior volume below the cavity and the boundary integral for truncating the finite element mesh on the antenna/cavity aperture [1, 2]. Extensive finite-element method formulations of the problem have been presented by various authors[l, 2, 3]. In the context of the FE-BI method, the mesh truncation condition between H and E is given by the boundary integral equation H = H9g~- 2ikoYo f G(r, r') (z x E(r')) dS' (1) so where G is the electric dyadic Green's function of the first kind such that n x G = 0 is satisfied on the metallic platform. For a cavity recessed in a ground plane, G becomes the

Plan....::: ii.......... -T"' — S Figure 1: Geometry of a cavity-backed slot antenna in a ground plane half space dyadic Green's function G= (- I VV 4RR (2) with R = Ir - r'i and I being the unit dyad. For this problem, H9~ is equal to the sum of the incident and ground plane reflected fields for scattering computations and zero for antenna analysis. The linear system generated by applying the standard FE-BI procedure is of the form [A] f {EV [0] [0] {E} {b(3) {E} [0] [L] {E } {b} Here, {EV} denotes the field unknowns within the volume enclosed by So, whereas {Es} represents the corresponding unknowns on the boundary So. The excitation column {bv} are due to internal sources and {bs} is associated with incident field excitations (for scattering). 3 The Fast Multipole Method The FE-BI system (3) is partly sparse and partly dense. More specifically [A] is sparse, whereas [B] is dense. Thus, although [B] is much smaller in rank than [A], it is usually responsible for most of the CPU and memory when an iterative algorithm is used for solving (3). To alleviate the CPU and memory requirements, fast integral methods have recently been introduced to perform the matrix-vector product [B]Es much faster and with much less memory. Two of these approaches are the adaptive integral method (AIM) and the fast multipole method (FMM). Both FMM and AIM reduce the CPU time and memory requirement from O(N2) down to O(Na) where a < 1.5. The main feature of AIM and FMM is the decomposition of the matrix as [B] = [BL3]ea + [B]far (4) based on some threshold distance referred to as the near-zone radius. The matrix [BI]near contains the interactions between elements separated less than the threshold distance, whereas [B1]f" contains the remaining interactions. The elements of [B3]nea' are evaluated without approximation. However, the product [B]farEs is evaluated in an approximate

manner leading to a much faster execution. Extensive formulations of FMM and its variations have been presented by [4, 5, 6]. To realize the FMM speed-up the computational domain is divided into M groups. The memory is then reduced to O(N2/M). Also, the CPU requirement of the FMM is O(NM)+O(N2/M)[5, 6]. This can be minimized by choosing M = VN and the result is an O(N1'5) algorithm. The memory required for the FMM also becomes O(N1'5). In practice, both the computational complexity and storage requirements of the FMM are less than those of standard FE-BI formulation for problem sizes larger than 1000 boundary unknowns. This makes FMM more suitable for problems involving large number of boundary unknowns. 4 Results and Comparisons We have applied the approximate approach to a 4x1 cavity-backed dipole array. Fig. 2 shows the layout of the curved array obtained by wrapping the flat array on a cylinderical surface. The geometry was run for several spacings of the dipoles and several different discretizations resulting in different numbers of boundary unknowns. The comparison of the radiation patterns for the flat and curved array is given in Fig. 3(a). The effect of curvature shows itself as an increase in the amplitudes of the side lobes and a decrease in the amplitude of the main lobe. Fig. 3(b) is the comparison of the convensional FEBI method with and without FMM. The feasibility of FMM is observable for all problem sizes of practical interest. Fig. 4(a) depicts the trend for the matrix filling time. Also Fig. 4(b) shows the time consumed per iteration in the iterative solver. We can conclude that the FMM implementation is both robust and reauires much less resources for large scale implementations. References [1] J.L. Volakis,T. Ozdemir and J. Gong "Hybrid Finite Element Methodologies for Antennas and Scattering", IEEE Trans. Antenna Propagat., pp. 493-507, March 1997 [2] J. Gong, J.L. Volakis and H.T.G. Wang, "Efficient Finite Element Simulation of Slot Antennas Using Prismatic Elements", Radio Sci., Vol. 31, No. 6, pp. 1837-1844, NovDec. 1996 [3] T. Ozdemir and J. L. Volakis, "Triangular Prisms for Edge-Based Vector Finite Elements Analysis" IEEE Trans.Antennas Propagat., pp. 788-797, May 1997 [4] R. Coifman, V. Rokhlin, and S. Wandzura, "The Fast Multipole Method for the Wave Equation: A Pedestrian Prescription", IEEE Antennas and Propagation Magazine, Vol. 35, pp. 7-12, June 1993 [5] W. C. Chew, J. Jin, C. Lu, E. Michielssen, and J. M. Song, "Fast Solution Methods in Electromagnetics", IEEE Trans. Antenna Propagat., Vol. 45, No. 3, pp. 533-543, March 1997 [6] S. S. Bindiganavale and J.L. Volakis, "Comparison of Three FMM Techniques for Solving Hybrid FE-BI Systems", IEEE Antennas and Propagation Magazine, Vol. 39, No. 4, pp. 47-59, August; 1993

40 20 10 0 -40 40 20 0 -20 0 -20 40 -40 Figure 2: 4 x 1 dipole array on a curved platform. 2 r.J.I 2000 1500 IE 1000 0 500 / FE-BI-FMM | | - FE-.81 &.................. / ',0,00 5020020 -80 -60 -40 -20 0 20 40 60 80 Iu 500 1000 1500 2000 2500 Number of BI Unknowns (b) 3000 3500 (a) Figure 3: (a) The radiation patterns of the planar and the curved array, (b) Comparison of the total CPU time consumed for the solution. - lZW,.-,1000'1000........... 800 - / FE-BI-FMM 0/ 0 — FE-BI 400 / 200 0 100 200 300 40 %$ 0 1000 2000 3000 4000 Number of BI Unknowns (a) 5000 6000 (b) Figure 4: (a) Comparison of the CPU time consumed for BI matrix filling, (b) Comparison of the CPU time consumed per iterarion.