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

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 4th Quarterly Report 035067-4-T October 1996 September 1998 Oct. 5, 1997 (4th 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.................................................5......................................... Executive Summary and Project Status................................... 6 1. Small array code (FSS-BRICK).................................. 6 2. M om ent M ethod Code (FSS-E ER)................................................................. 7 3. Finite Element Planar Periodic Code (FSS-PRISM)............................................ 10 4. Doubly Curved SERAT code (FSS-CURVE)............................... 13 Project G oals....................................................................... 15 SUMMARY OF CODE FEATURES IN FSS-BRICK AND FSS-PRISM.................... 17 Non-Commensurate FSS Modeling in FSS-PRISM.......................................... 20 EIGER PROGRESS REPORT................................................................. 24 The Preprocessor............................................................. 24 EIGER and Postprocessor Functions.................................. 26 Modeling Non-Commensurate Structures.................................................. 28 Code User's Manual................................................................. 30 Hybrid FEM Code with the Fast Multipole Method for Curved SERAT Panels............. 31 LIST OF FIGURE CAPTIONS IN SUMMARY SECTION Figure 1. SERAT code development plan......................................... 8 Figure 2. Illustration of the scaling approach.......................................................... 11 Figure 3. Results from FSS-PRISM for non-commensurate......................................... 12 Figure 4. Grouping approach for implementing the FMM method in PRISM array........ 14 Figure 5. Matrix-fill CPU speed-up when the FMM is used in PRISM for array analysis.................................................................................................................................... 14 Figure 6. Illustration of the SERAT panel (Planar and Curved).................. 16 3

CHRONOLOGY of Events (Updated Every Quarter) * Anril 1996 Pronosal Siubmission ",t -X v vX Y v s %J v ss Y ss 4 - xA XL7L AI * July 1996 * August 1996 * 20 Sept. 1996 UH) * October 1996 * October 1996 * 15 Nov. 1996 * 9 January 1997 * 28 February 1997 * 5 April 1997 * 30 May 1997 Answers to Proposal Questions Began Contract Negotiations Kickoff meeting at Ann Arbor (attended by Sanders, UM and Contract Signed between U-M and Sanders in Mid October Subcontract to the Univ of Houston (formalized in early November) SERAT Review meeting (at Nashua) 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 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 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. 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. SERAT review Meeting (Nashua) Submission of Third Quarterly Report 0 * 25 June 1997 5 July 1997 4

* 6 Oct 1997 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. 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. MEETINGS No meetings were held during the 4th quarter (July 1- Sept 30, 1997) 5

Executive Summary and Project Status This is our 4th progress report and coincides with the completion of the first year on the project. At this time we are on schedule with all tasks and in some cases we are ahead of schedule as mentioned in our 3d quarterly report. In fact our third quarterly report provides a good assessment of the project status and will not be repeated here. To reiterate, the goal of this project is * to develop and deliver a set of codes from simple/faster to complex/slower capable of characterizing the SERAT array panels (with slot, strip and crossed dipole elements) on planar and doubly curved platforms. * the unique feature of our proposed effort was to employ well-established methods for the simpler planar arrays on FSS panels (low risk approach) and new more capable hybrid methods (higher risk) for finite/infinite planar and double curved arrays with possible non-planar inhomogeneities. The hybrid FEM method was proposed for this implementation using prismatic elements for ease in meshing and fast integral truncation schemes. As given in the simplified code development plan in Figure 1 (see also the Milestone chart in Figure 2), four different codes were proposed to address the speed and complexity requirements of the project. 1. Small array code (FSS-BRICK) FSS-BRICK is an upgrade from an existing code developed under Air Force sponsorship. It employs bricks for modeling the FSS and antenna volume in the context of the finite element method (see Progress reports 035067-2-T and 035067-3-T) and a fast boundary integral technique for mesh truncation. The intent is to use this code for fast finite array computations by compromising somewhat geometrical fidelity due to the use of rectangular bricks. Moreover, this code served as the reference to the future periodic codes FSS-EIGER and FSS-PRISM. FSS-BRICK was extended to handle multilayered FSS structures and a geometry driver was added for ease of use. In accordance with the milestone chart, this code was finished on schedule at the end of the first year. A users manual was prepared including sample I/O files with example results and validations for finite arrays on multilayered SERAT panels. Apart from a lack of full geometrical generality, FSS-BRICK permits modeling of * resistive sheets and cards at each dielectric layer * horizontal and vertical loads at each node/edge of the model * vertical and horizontal probe feeds * strip dipoles, slots and crossed elements separated by thin laminates. * vertical wire or strip feeds using a concatenation/combination of conducting posts at any edge/node location within the domain (both, vertical and horizontal loads can be handled) 6

Typically, any geometry of FSS or antenna element can be handled using FSS-BRICK provided the geometry fits within the rectangular grid specified by the discretization rate. However, for ease of use, direct use of the Driver permits a list of element choices (primitives), thus, bypassing a need for geometry meshing. However, this SERAT driver can be bypassed for a more general use of the code. FSS-BRICK permits scattering and antenna calculations for metal backed panels as well as transmission coefficient calculations for the FSS panels when no metal-backing is used Since FSS-BRICK handle only finite arrays, non-commensurate FSS panels can be handled with equal ease within its inherent gridding limitations 2. Moment Method Code fFSS-EIGER) This code was developed by the SERAT team members from the Univ. of Houston. It is based on the traditional moment method technique with the inclusion of the periodic Green's functions to handle the multilayer FSS and antenna panels. Basically, each periodic antenna element is modeled using triangular surface patches and the periodic Green's functions is used to radiated and couple the patch currents among the different periodic layers. By using separate periodic Green's functions for each FSS layer, different periodicities (non-commensurate) for each of the individual layers can be handled in a rather straightforward manner. Such an approach has been shown to be more effective than the supercell method which although more accurate, enlarges the computational domain and is furthermore restricted to layer periodicities that lead to rational ratios. Univ. of Houston decided to begin development of the SERAT moment method code starting with the existing EIGER code. EIGER was developed by the Univ. of Houston in conjunction with Sandia Labs. and has the unique feature of being quite modular (object oriented). The Green's function is provided through a separate object/module. For the SERAT application, it was modified to include the periodicity of the array and the FSS panel. For non-commensurate structures each layer will be associated with its own Green's function. FSS-EIGER for commensurate FSS is completed and delivered to the Univ of Michigan. This code package includes: * SERATBUILD: a preprocessor allowing the user to bypass the line by line raw input to EIGER. Using this package, the user specifies the element types (slots, strip dipoles etc.) for each of the antenna or FSS layers. As in the FSS-BRICK and FSSPRISM, meshing of these antenna/FSS elements is done internally without intervention by the user. Specification of the usual discretization rates along each direction is, of course, necessary. As noted above, by resorting to a set of 'primitive' shapes for the antenna and FSS elements we have eliminated a need for a complete integration of these packages through a single meshing package. SERATBUILD generates a input file which is then used as input to the FSS-EIGER code. 7

*FSS-EIGER: This is the moment method code for simulating the FSS and the periodic antenna elements. It was upgraded from the original EIGER code with the inclusion of periodic Green's function. It is import to note that the periodic Green's function for the free-space layer employs an acceleration scheme for summing the modes which proved very crucial to developing a code that is rid of the usual convergence difficulties associated with most periodic Green's function. The theory of the formulation is described in the first quarterly report (Univ. of Michigan Radiation Lab report 035067-1-T). Both reflection and transmission coefficients of the FSS can be extracted from EIGER and in contrast to FSS-BRICK, the elements can be of any arbitrary orientation. However, FSS-EIGER does not allow modeling of vertical inhomogeneities as is the case with vertical wires, feeds, lumped loads, etc. However, EIGER does allow modeling of resistive cards in the plane of each antenna or FSS layer. The users manual for FSS-EIGER is currently being written and will be delivered along with the code at the end of October 1997. At the same time, Univ. of Houston is working on completing the non-commensurate implementation using the method noted above and described later in more detail. In accordance with the milestones chart, this code was to be developed at the end of the first year and work on upgrades during the 5th quarter. Thus, its development is still within our original schedule plans. Year 1 Year 1 Year 1 and 2 Year 2 Figure 1. SERAT code development plan. 8

Milestone Chart for EM Model Development (Tasks 1 and 2) Ouarterlv Progress 5th 6th 7th 8th Q. Q. Q. Q -- --- v --- —:_ I- I hi p= hi - -, - p- - Planar Periodic Array I oftware Integration ana 1/u Di I Software Support 9

3. Finite Element Planar Periodic Code (FSS-PRISM) FSS-PRISM refers to the hybrid FEM code. This code is the centerpiece of the proposed development and combines the geometrical adaptability and generality of the FEM for modeling materials and inhomogeneities with the rigor of the moment method for mesh truncation. This code is therefore the most general for planar SERAT panels. As can be understood this code is the slower of the FSS-BRICK and FSS-EIGER codes for planar SERAT panels. However, we plan to improve the CPU speed of FSS-PRISM with the inclusion of the adaptive integral method(AIM) for handling the matrix vector products in the solver. This will be done in the 5th quarter. The development of FSS-PRISM began in March 1997 as a modification of the nonperiodic finite element-boundary integral PRISM code, originally developed under Air Force sponsorship. It was chosen for the SERAT project because its prism elements permitted simplifications in specifying the geometric of the FSS panels and the antenna elements. Basically, the constant thickness layers used for building the SERAT panels eliminates a need to use the more general tetrahedral elements PRISM was upgraded to FSS-PRISM with the addition of the periodic Green's function obtained from FSS-EIGER and the implementation of periodic boundary conditions to reduce the computational domain to within a single periodic element. These modifications along with the inclusion of options for metal-backed or open FSS panels provided us with the first version of FSS-PRISM at the end of the 3rd quarter. Thus FSSPRISM was developed one quarter ahead of schedule. Validations and the theory relating to FSS-PRISM were described in the 3rd quarterly report (Univ. of Michigan Radiation Lab Report 035067-3-T). Having completed the first version of FSS-PRISM, during the 4th quarter we proceeded * to develop a geometry Driver for FSS-PRISM, where again the user will specify primitive shapes (slots, strips, crosses, etc.) to be internally used for developing the mesh without user interference. The user can specify such primitive shapes at each layer interface and distorted prisms are used for connect the nodes from one layer to the next. At this moment, the Driver can handle primitives which are oriented along the two principal directions and will be generally in October to handle more arbitrarily oriented elements. In addition, although FSS-PRISM can handle resistive cards, lumped loads, etc., these items have yet to be included in the Driver. * to generalize FSS-PRISMfor modeling non-commensurate FSS panels. Since no studies on this topic are available in the literature, it is necessary to pursue a research investigation. As described later in the report, several approaches were explored to find a method for modeling non-commensurate FSS. Among them, the supercell method, element scaling and layer periodicity decoupling were examined. As noted earlier, the supercell approach, although yielding exact results for a limited class of problems, it was not found attractive due to its large CPU requirements and 10

limited flexibility. The scaling approach is a theoretically sound method and was rather straightforward to implement. With this method, the periodic cell's vertical boundaries were adjusted (shrunk or expanded) to accommodate the periodicities of each layer. The scaling method has already been implemented into FSS-PRISM and a sample result is given in Figure 3. The comparison is made with data produced by Aroudaki and Hensen using the layer periodicity decoupling method. As seen FSSPRISM data predict the similar FSS behavior. At this moment, we are not as impressed with the CPU performance of the scaling approach for non-commensurate FSS. Substantial distortion of the prismatic FEM elements leads to slow convergence, a situation that also impact accuracy. Given that delivery of FSS-PRISM is not scheduled till the end of the 5th quarter, we are now looking at two possible ways to improve modeling of non-commensurate FSS. There are: ~ use of more robust iterative solvers Our intent is to use the preconditioned Flexible GMRES and work toward this direction is already under way. * revert to an implementation based on layer periodicity decoupling This implementation is based on the same principle to be used in FSS-EIGER and will likely lead to much faster convergence. However, its implementation is more complex and time consuming. In addition to the above tasks, during the 5th quarter we will also incorporate a fast integral method in FSS-PRISM for further CPU speed-up. L -- - Figure 2. Illustration of the scaling approach 11

h/p 9 9 9 9 9 27 m -........ 1t gtm g 41m 9 18 18 18 18 g 9r c/p m 6 14 14 14 14 mf r Er= 2.3 1- j.... v - - - - - - - - 1 0.9 t 0.8 0.7 0.6 0.5 0.4 0.3 E 0.2 I-.1 u 0 1000 2000 3000 4000 5000 6000 7000 8000 Frequency in GHz > Figure 3. Results from FSS-PRISM for non-commensurate 12

4. Doubly Curved SERAT code (FSS-CURVE) We began development of the curved SERAT code during the 4th quarter in accordance to the milestone chart and project plan in our proposal. As stated in our proposal, two approaches will be pursued. * One is an approximate method and will rely on results from the planar FSS code (FSS-PRISM, FSS-EIGER, FSS-BRICK). This implementation is really a modification of these planar FSS codes and will be carried out during the 6th quarter as a last resort. Basically, this is an approximate method based on the relocation of the unit cell currents (with the appropriate phase terms) as computed from the planar FSS codes. * The mainstream approach that was proposed is based on an exact implementation of the curved panels. Our approach combines the finite element method for interior modeling and fast integral methods for mesh truncations. Practically, this approach combines the geometrical adaptability of PRISM with the speed of fast integral methods such as the fast multipole method(FMM) During the 4th quarter we began the integration of the non-periodic PRISM with the fast multipole method. Our goals were * to examine the speed-ups attained with this method * to find the appropriate region (array size) for applying fast integral method * to examine the practical CPU requirements for typical SERAT problems * to provide a validation for the approximate curved SERAT code As of this moment, we have completed a first version of the PRISM code incorporating the FMM and have obtained several validations for large planar arrays. An example 6x6 antenna element array with the employed FMM grouping approach is displayed in Figure 4. The corresponding CPU time (fill time) is given in Figure 5 (this example corresponds to 2500 boundary integral unknowns) and shows a factor of 3 for this example. Since PRISM is capable of curved FSS modeling the generalization of this code to curved FSS panels is straightforward. This code is scheduled for completion at the end of the 6th quarter and we will report on it with more details in the 5th and 6th quarterly reports. 13

... " '......... ".".. lei W c t \Lf I L/l.el s; L t le.e s;-__ CFU Ti me Consumed for BI Matrlx Construction 1000 800 700 r1 400.....- -- - -.......L.. 200 ___, --- '" _._. --- — -, * —. ~ $,- t,, t $, t,.,t. Ot '.'.................... '....i 0 I 5000 I 500 2000 2500f Nurber of 9I1 Unknowns \anl yis/I/" I Figure. Mroix -fil lo iCPUlspeed-up i whe th F isu d in PRISM fr array aInalysis.r ot lf 0 - ---- /D t VIVI VIA "VIA VIAzi Vi, 7 0 0 - I FM 700 E 50 0 - |L 400 - 3l 300 - 200 - - =.J 2 - - 0 500 1000 1500 2000 2500 Number of Bl Unknowns Figure 5. Matrix-fill CPU speed-up when the FMM is used in PRISM for array analysis 14

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, crosseddipole, 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 FSSPRISM, 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 non-commensurate 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 15

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 A itilil -tlicll' *a ' Id.;l ---'I ~t tl ) Flfl tcflfl::*':i-m $: i.i. a t........ LI.E It.f An *|W 3^ 1^ *^'liijg^.. -.. I.-.. _~)"" "n s "ix~ i I... NI"r~p~.......~ + x Figure 6. Illustration of the SERAT panel (Planar and Curved) 16

Summary of Code Features in FSS-BRICK and FSS-PRISM In the project proposal, the following input/output (I/O) operations are described for the SERAT codes to be developed: 1) Automatic meshing of the required FSS and antenna elements in multilayered media should be available in the codes. Bricks and distorted prisms can be used for automated meshing. 2) The user can enter the number of layers above and below the antenna elements; the number of FSS layers, each layer's physical and material characteristics; the element type, its location and dimensions; periodicity; resistive sheet specifications; feed type and location; lumped load values and location; excitation information, and so on. 3) Strip dipoles, single slot elements, dual polarized (crossed) dipole and slot elements should be available as FSS and antenna elements. The crossed-dipole antenna element can be specified as two dipoles separated by a thin laminate. These elements can be oriented in any arbitrary direction at any layer interface. In addition, the automatic mesher can include square patches as an option element. 4) The user can have options of placing vertical and horizontal complex lumped loads between any two node location in the structure. 5) Complex resistive sheets can be placed at any location on the specified interface, between the elements and on the elements. 6) Horizontal as well as vertical feeding can be specified at any node location using current probes or potential specifications. 7) Modeling of a metal-backed or a free standing FSS panel should be available in the code. An option to model the bottom of the FSS as a perfectly matched absorber may also be included. Also, impedance boundary conditions may be allowed at all walls, in the vicinity of platform surface and on the base of the SERAT structure. 17

8) The user will have the option of defining the geometrical features of the feed/balun model along with its electrical composition. 9) The above input options are supplied by a driver in the code. Before proceeding with the final analysis, there will be a display of the input geometry and mesh for verification by the user. The output file should list all input parameters and computed results and will be available in graphical and tabulated form (on screen and postscript). The graphical output should include the capability to plot/display: a) Antenna gain as a function of angle, polarization, and frequency; b) Input impedance (VSWR) vs. frequency; c) Mutual coupling between any two elements in the array vs. frequency; d) Scattering as a function of incidence and scattered angles, polarization and frequency; e) Antenna gain or scattering at specific angles and frequency vs. variations in design parameters; f) The capability of computing the reflection (or transmission) coefficient at each FSS layer. The geometry driver development status for FSS-BRICK and FSS-PRISM codes are summarized in Table 1. Most of the I/O operations described above are included in FSS-BRICK Driver. After addition of the plotting package (XMGR), FSS-13RICK will be delivered with its user's manual in October 1997. FSS-PRISM Driver is still in progress, and the manual will be ready in November 1997. 18

Input Features [FSS-BRICK (Finite Array) FSS-PRISM (Infinite Array) Antenna Elements Dipole V/ Slot / V/ Crossed-dipole V 3rd week of Oct'97 FSS Elements Strip / / Slot V/ V Crossed-strip / 3rd week of Oct'97 Crossed-slot V 3rd week of Oct'97 Arbitrary Orientation N/A End of Oct'97 Resistive Cards (surface) V/ 2nd week of Oct'97 Lumped Loads (between any two grid location) Vertical (z-directed) V 2nd week of Oct'97 Horizontal (x- or y-directed) V/ 2nd week of Oct'97 Feed Wires (or Short-ckt. Pins) v/ 2nd week of Oct'97 Sources (at any grid location) x-directed Probes / 2nd week of Oct'97 y-directed Probes / 2nd week of Oct'97 z-directed Probes V 2nd week of Oct'97 Output Features Radiation Pattern (with the closed backing) Antenna Gain V 3rd week of Oct'97 Input Impedance V 3rd week of Oct'97 (as a function frequency, and scan angles) Scattering Pattern (with the open backing) Transmission Coefficient V/ Monostatic RCS V 2nd week of Oct'97 (as a function frequency, and scan angles) Output Display Format Tabulation l/ V/ Plotting (XMGR) 1st week of Oct'97 3rd week of Oct'97 Manual Almost done 1st week of Nov'97 Table 1: Geometry Driver development status.

Non —Commensurate FSS Modeling in FSS-PRISM T. F. Eibert In the hybrid finite element (FE) method the actual periodic structure is modeled by finite elements and the solution domain is terminated using a boundary integral formulation for the top and bottom half spaces, respectively. If unique periodicities can be defined for the whole structure, it is sufficient to only model one unit cell for an exact analysis. For non-commensurate structures, one has to account for different layers each having different periodicities. If the individual periodicities of the layers are integer multiples of each other the so-called supercell approach can be employed. This approach combines groups of several FSS elements in a single periodic cell which runs vertically through the FSS structure. However, for irrational ratios the periods of the supercell tend to infinity and an exact solution is no longer practical. Nevertheless, it is possible to construct approximate supercell periodicities of finite size which can give reasonably accurate results. In general, the supercell approach is very time consuming. So, our goal is to work with FE meshes that can employ different periodicities for each layer as it is illustrated in Fig. 1. This approach in general introduces some approximation,// X M L[ X X X\ I/ / i l \ \ \ \I Figure 1: Scaled PRISM mesh for a periodic cell involving non-commensurate FSS layers into the model, but it is very flexible and can easily be adapted to meet arbitrary requirments on the performed analysis. Within the PRISM-BI code the surface mesh structure in the individual layers is fixed. However, a variation of the periods can be relatively easily introduced by scaling the surface meshes in each layer as be seen in Fig. 1. A difficulty with this approach is that the periodic boundary 20

condition for the scaled mesh must be implemented for vertical edges which are at skew angles with each other. Thus, some distortion of the fields near the boundaries will be introduced by scaling the mesh. The differences between the scaling factors in the various layers should not be very large to maintain accuracy. However, as mentioned above, the amount of scaling can easily be controlled by grouping more or less unit cells within the various layers of the periodic structure. Another very interesting effect of the scaling approach is that the resulting fields must be also scaled. This can be illustrated by the following reasoning. If we assume a plane wave in the upper half-space is travelling towards the periodic structure a certain amount of power enters the FE mesh on its top aperture. The amount of this power is determined by the power intensity as well as the cross sectional area (parallel to the layers) of the FE mesh. A scaling of the FE mesh changes the cross sectional area of the mesh dependent on the mesh height. Thus, a fixed amount of power is maintained over the cross sectional areas of the periodic cell as the wave travels through the mesh. This causes changes to the power intensity as well as the associated fields at each layer. For a correct modeling of the interactions between neighboring layers of the mesh this is not an issue and for an evaluation of reflection and transmission coefficients the scaling of the fields h/grn g/gm c/gm 9II 9 6 9 9 9 9 -4 - Li 18 14 18 14 18 14 12 8 l E r= 2.31 ij..... F m Figure 2: Geometry of 5-layer FSS can easily be removed. In the numerical solution of the discretized field problem the scaling of the FE mesh can result in higher numbers of iterations to achieve convergence when an iterative solver is employed. This is due to the broader eigenvalue spectrum that is typical for irregular meshes. Some compensation of this effect can be obtained by the introduction of small losses in the dielectrics of the layers, a situation that is usually true for practical problems. As an example a 5-layer FSS structure with non-commensurate periodicities was 21

investigated as illustrated in Fig. 2. The structure was presented in [1] and it was used as a low-pass filter for radioastronomical measurements. The mesh used for the calculations with the FSS-PRISM code consisted of fourty prism layers and comprised a total of more than 50 000 unknowns. In the uppermost layer four unit cells were meshed so that periods of 18 um resulted for this layer equal to the periods of the layers below. In the lowest layer also four unit cells were meshed, but this resulted in periods of 24 jm for this layer. Thus, for all layers of the mesh scaling factors of 1.0 could be employed except for the lowest layer where a scaling factor of 1.333 was needed. In the transition region between the lowest layer and the immediate layer above the scaling factor was adapted gradually. First, as a sanity check all layers were assumed to be air and the metallic patches were not present. Hence, we had a situation as illustrated in Fig. 3 for a commensurate structure which transmission coefficient should be exactly 1. The results Air BI/PGF * I I I PBC a, Air, PBC I I * I I.~ Air \ BI/PGF Figure 3: Infinite air layer for the scaled mesh are given in Fig. 4. It can be seen that the transmission coeffecient exhibits errors of some percent due to the vertical edges at skew angles with each other at the periodic boundaries. As can be expected, the error decreases for higher frequencies. In Fig. 5 results are given for the 5-layer low-pass filter. One curve shows the results from [1] which where obtained using a MoM code based on spectral domain Green's functions. These results were also validated by measurements. A second cuve in the figure shows results of the PRISM code being calculated with a scaling factor of 1.0 for all layers. So, the periods of the lowest layer were 18 gm as in all layers above. The transmission factor for this commensurate filter shows similar frequency behavior as the reference curve from [1]. However, the transmission factor in the passband is considerably larger than in the reference curve. At this point it has to be mentioned that the results in [1] were obtained with relatively high dielectric losses for the layers. The values of the losses are not given in the paper. The calculations with the PRISM code were performed with tan b = 0.02. The third curve in the figure gives the results for the scaled PRISM mesh. It can be seen that its transmission factor is closer to the results of [1], but the irregu 22

lar behavior of the transmission factor with respect to frequency also reveals that there are some convergence problems in this case. 1 i I I I I 0.95.L-. -.....'... —.. -.. - - — ' 0.85 ------- -I ---- ------ ------ ------ ------ 0.8 --—. ----. -----. —. -. --- —. -------- 0.75 -- ------------- e 0.7 ------- ------ ------ ------ ------ r --- —-— T F e i G F r, T c ' cek 0.6 0. — i --- —----------— L — 4 0.8 " _ — -— _ 0.5 ------- ------ - ---— L — -- ---------— 1 --- —-. 0. ------------ --- ---: PRIS,,com 0...... a)04 --- —-------— ^-, --- —-— _-^ —I_ 0.6-.. — --—....-.......-..-............ 0.1 i Figure 5: Transmission coefficient (5-layer FSS) References [1] Aroudaki, H., Hansen, V., "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, pp. 1486-1491, Dec. 1995. Dec:*.l 1995 23

EIGER PROGRESS REPORT This report summarizes the features of FSS/EIGER, as currently implemented. Some changes and improvements are to be expected, particularly regarding input and output data formats, but the basic computational engine for periodic structures and the preprocessor for setting up the problem description are essentially complete. Additional features that have not been implemented, but which may be considered for possible future extensions to the code, are also mentioned here. The Preprocessor Preprocessor organization FSS/EIGER currently comprises two codes: SERATBUILD, a preprocessor for setting up problem description for periodic structures, and EIGER, a general purpose computational engine with integrated postprocessor for deriving figures-of-merit for periodic FSS and array structures. The preprocessor has been designed so that it may be integrated with EIGER by changing only a few lines of code. However, the codes have been left separate at this time for the following reasons: * It is anticipated that a number of changes or additions to the preprocessor may be required to improve the user interface. Furthermore, during the testing phase a number of different users or testers may need access to the user interface source code. Since the preprocessor has no licensing or access restrictions, it was separated to facilitate such access. * The process of generating the geometry and problem specification file is often an iterative process which the user must go through before obtaining a correct data file for running EIGER. This is particularly true during the learning phase in running the code(s). Separating the preprocessor therefore quite naturally separates this problem definition phase. As a user becomes more familiar with the code, it may be desirable to integrate the preprocessor with EIGER, and this can be done by changing only a few lines of code. In most computational electromagnetics applications, it is desirable to also separate the postprocessing software from the preprocessor and computational engine. Primarily this is because once one has solved for the currents on a structure (electric and equivalent magnetic) at a given frequency, one may then examine afterwards any number of quantities, such as pattern cuts, near fields, surface fields, impedances, etc. Keeping this process separate also separates the field observation location data from the problem geometry and excitation data. 24

Indeed, this approach was originally taken with FSS/EIGER since some of the needed postprocessing software already existed. However, it was later decided to integrate the postprocessing phase with the computational phase for the following reasons: * For periodic structures, the most important observables are simply the amplitudes and phases of a few of the excited Floquet harmonics (and, for radiation problems, the input impedance); practically any other desired data or figures-of-merit may be easily computed from this data. Furthermore, all the data needed to specify and determine such quantities - such as the wavenumbers of the Floquet modes - is immediately derivable from the excitation data. Hence, FSS/EIGER always assumes that these quantities will be of interest and automatically computes them without the need for user specification. * From a software point of view, 95% of the software needed to perform the postprocessing function resides within EIGER. Indeed, most of the postprocessing can be performed with the addition of only one or two subroutines. Hence, from a resources point of view, it also makes sense to integrate this function within EIGER. Preprocessor functions The functions and features of the current preprocessor, SERATBUILD, may be summarized as follows: * SERATBUILD constructs planar, rectangular or cross-shaped elements which are assumed to be metallic - unless they lie in a ground plane, in which case they are assumed to be apertures (slots). Although EIGER is not generally restricted to planar elements, the current version for periodic structures is so limited because Green's function acceleration procedures have not yet been developed for out-of-plane geometries. Most of the other modeling restrictions which follow are not due to limitations of EIGER, but of the preprocessor. * All FSS or antenna elements are constructed from triangles. This choice was initially dictated by the desire for flexibility in geometry modeling, as well as the perceived need to generate numerical Green's functions for Michigan's PRISM code. For the strictly rectangular geometries proposed for SERAT, however, it would probably be computationally more efficient to use rectangular patches. Indeed, EIGER accepts rectangular patches as input, and this feature has been successfully tested, but unfortunately the input data to EIGER must be generated by other means; SERATBUILD only generates a triangular mesh. * The element arrangement may be arbitrary, i.e., the periodic lattice may be composed of rectangular or parallelogram unit cells. Specifying a rectangular geometry simplifies data input. * While the maximum number of layers is currently dimensioned within SERATBUILD to 15, a simple change in the dimensions allows an essentially arbitrary number of material layers to be considered. Layers may have complex 25

relative permittivities and/or permeabilities (choosing unit permeabilities simplifies data input). The boundaries between any of the layers may be specified as ground screens and any ground screen may contain a slot or aperture of rectangular or cross shape. * For antenna modeling, one may specify delta-gap voltage sources whose terminals establish a potential difference of arbitrary complex amplitude between any number of pairs of adjacent triangles of the mesh. * For FSS modeling, incident plane waves of arbitrary polarization (TE or TM to the direction normal to the layer boundaries) and arbitrary angle of incidence may be specified. The amplitudes of the incident waves are automatically chosen to produce unit-amplitude transverse electric field components, so that similar components of the scattered Floquet modes may be directly interpreted as scattering coefficients. * A desirable feature which is currently missing in FSS/EIGER is an option to compute a structure's response to both incident polarizations (TE and TM) at the same frequency and incidence angles, without having to recompute the impedance matrix. Also missing is the ability to specify arbitrary reference planes for measuring reflection and transmission coefficients (i.e., Floquet modal coefficients). Currently, the z=O plane is always taken as the reference plane. EIGER and Postprocessor Functions This section summarizes the current features and functions of the EIGER computational engine and the integrated postprocessor. They are as follows: * The unknown quantities computed by EIGER are the electric and magnetic currents induced on metallic and aperture (slot) elements, respectively. * A sophisticated acceleration procedure is used to evaluate the layered-media Green's functions for the scalar and vector potentials. The Green's functions are initially formulated in spectral form as a doubly infinite series. When appropriate, the direct and (up to) three quasi-static periodic image contributions are extracted from the spectral representation. These contributions are efficiently evaluated by the Ewald method. Under a variety of conditions, however, the remaining spectral series is still slowly convergent and thus dominates the computation. * The active element input impedance and admittance are computed whenever only one voltage source is specified for a problem. * Both polarizations of the specular component of the scattered or radiated fields are computed from the induced electric and magnetic currents. In the plane-wave illumination case, because of the normalization of the incident field, the amplitudes of the scattered fields are equivalent to the generalized reflection coefficients (corresponding to the field scattered into the TE and TM specular waves from a TE or 26

TM incident wave). In the antenna problem, the fields radiated by the induced currents may be directly related to the active element pattern. One additional output item which would be desirable is the calculation of reflection and transmission coefficients for higher order (including cutoff) Floquet modes. This information would provide some assessment of effects of grating lobes whenever several modes are propagating, as well as provide data for analyzing cascaded FSS structures by the generalized scattering-matrix method. Other features which might be considered in future development plans are the following: * The specification of resistance cards. There are two problems which precluded developing a treatment of resistance cards during the current effort: 1) It is extremely difficult to specify and mesh the geometry to the detail needed. This is particularly true if the resistance cards overlap the conductors so that the mesh must conform to the boundaries of the overlap region. This modeling would also be simplified if the cards and conductors still maintain a rectangular geometry so that rectangular subdomains can be used. 2) The modeling procedure may depend on knowledge of the thicknesses and conductivity properties of the resistance cards, especially if it is deemed necessary to model the region of overlap region between cards and conductors. * The modeling of feed lines. Only EIGER's current restriction to handling planar equivalent currents in the periodic case prevents the modeling of feed lines which are non-planar. An extension of the mixed-potential integral equation approach to nonplanar periodic sources must be developed to implement this feature. Wire-to-surface junction modeling features are already available within EIGER. * Further acceleration of the remaining spectral series in the Green's functions. Several approaches could be considered here, among which are the following: * Since singular contributions have been removed, the resulting residual Green's function should be relatively smooth and hence could possibly be efficiently evaluated by interpolation from precomputed tables. * Since exact image representations (in complex space) for current sources in half spaces are available, these additional image contributions could also be extracted. This would result in Green's functions that always decay exponentially, even when the FSS elements are at an interface (the most typical case). The complex half-space image theory could be combined with the Ewald acceleration scheme to efficiently evaluate the resulting periodic 27

half-space Green 's functions. If adjacent layers are thin, however, the presence of nearby higher-order images may partially defeat this approach. * The calculation of self impedance and mutual coupling between isolated antenna elements in the presence of the FSS layers. This is an extremely computationally-intensive calculation, involving integration over all wavenumbers in a unit cell of the reciprocal (wavevector) lattice, in addition to the double summation in the spectral series. Considerable effort would be required to implement this feature. Modeling Non-Commensurate Structures The analysis of noncommensurate structures has been considered, where the ratios of the periods of the FSS layers are not rational numbers. The presence of noncommensurate FSS layers introduces significant challenges into the analysis that are not met in the commensurate case. For a commensurate structure, where the ratios of all periods are rational numbers, it is always possible to define a superperiod that contains an integer number of elements of each FSS layer. Therefore, the composite structure can, in principle, be modeled in the usual way as a periodic structure with a period equal to the superperiod. A structure consisting of noncommensurate FSS layers, on the other hand, is in actuality a nonperiodic structure, which never repeats. That is, there is no well-defined unit cell for the composite structure. Therefore, strictly speaking, the concepts of "reflection coefficient" and "active scan impedance" are not well-defined. Ultimately, however, the solution of any FSS structure by a periodic moment-method code (such as FSS/EIGER) must rely on the successful decomposition of the problem into one or more periodic structure solutions. Several different approaches to this have been considered. These are outlined below. 1) The simplest method is an approximate one, similar to methods that have been used in the literature by Vardaxoglou and others [1], [2]. In the proposed method, the currents on each FSS layer are assumed to be periodic, with a uniform progressive phase shift corresponding to the incident wave (or voltage excitation). This is an approximation, since the composite structure is not really periodic. The field produced by each periodic FSS layer is then found using a periodic Green's function, using the periodicity appropriate to the particular layer. The field at any observation point is then a superposition of the fields radiated by the different FSS equivalent currents, each one having the form of a periodic summation of Floquet modes (with a different Floquet mode set for each source layer). The resulting superposition of fields would not be a periodic function with a uniform progressive phase shift. Therefore, different results for the solution of the equivalent currents will be obtained depending on which unit cell in each layer the testing procedure is performed in the moment method solution. 2) A perturbational method has been considered. This method can be viewed as an extension to the previous method. In this method, a zero-order perturbational solution for the FSS currents is found first. This is the solution that is obtained by assuming a uniform progressive phase shift on all FSS layers. In this zero-order all Floquet 28

modes are accounted for when calculating the reaction between the field of a particular FSS layer and the currents on the same FSS layer. However, only a single Floquet mode, corresponding to the incident wavenumber, is used to calculate the reaction between different FSS layers. After the zero-order approximation is obtained, a first-order perturbation for the FSS currents is found by including the interaction of the Floquet modes launched by the zero-order currents, accounting for all Floquet modes except the ones corresponding to the incident wavenumber (since this wavenumber has already been used to find the zero-order currents). The process can then in principle be repeated to obtain higher-order perturbational corrections. This method has the advantage that it would be somewhat compatible with the existing EIGER code, and would probably not require extensive code modifications, at least if the solution were restricted to the first-order solution only. If the solution process was automated to allow for an arbitrary number of perturbations, the coding could become much more involved. The disadvantage of this method is that it is not exact, unless the number of perturbations were allowed to increase indefinitely. However, this method would be an improvement over the first method discussed above, and even the first-order perturbation may be sufficiently accurate for FSS layers that are not too closely spaced. 3) A "ray-bounce" method has been considered. In this method the Floquet modes (plane waves) that result from multiple internal reflections and transmissions from the different FSS layers, starting with the incident plane wave, are tracked and summed to yield a total solution in terms of a mixture of Floquet modes. Each time a Floquet mode impinges on an FSS layer in the solution process, the mode spawns a new set of Floquet modes that have a periodicity corresponding to the FSS layer. Thus, the total number of Floquet modes could potentially rise very quickly in the numerical solution. However, only Floquet modes with an amplitude above a specified threshold would be kept in the tracking process. This method has the potential of being accurate, with a user-controlled convergence, defined by the threshold value. This method would probably be the most efficient for FSS layers that are relatively far apart, and for periodicities that are small enough so that only the incident wavenumber corresponds to a propagating mode in each region. 4) A "discrete spectral propagator" method has been considered. In this method the entire two-dimensional transverse wavenumber space is discretized into a set of plane waves, both propagating and evanescent. A generalized scattering matrix is then found for each FSS layer. This involves treating each plane wave in the discrete wavenumber set as an incident plane wave, and then finding the amplitude of the scattered plane waves in the discrete set by solving a periodic-structure scattering problem. Once the generalized scattering matrix has been obtained for each FSS layer, the matrices may be cascaded to obtain a generalized scattering matrix for the composite structure. This method is quite general and does not assume anything about the periodicities of the FSS layers. In fact, the same methodology may be used to solve for the scattering from a composite structure consisting of arbitrary scattering bodies, not necessarily periodic FSS layers. In principle, the method is limited in accuracy only by the density of the sample wavenumber space. In practice, the computation time may become prohibitive for even moderate densities. The 29

computation time should not depend on the periodicities or the spacing between the FSS layers, however. The present plan is to implement method (1) above, which is the simplest method, but also the most approximate. The implementation of the other methods, which are potentially more accurate, but also more computationally expensive, could be done in the future. In the implementation of method (1), the location of the unit cell on each FSS layer that will be used for the moment-method testing procedure can be specified by the user. The resulting solution for the equivalent currents on the structure will depend on the choice of the unit cell locations. The user will thus have the ability to study how the solution varies with the testing cell locations. The variation of the reflection coefficient with unit cell locations may give an approximate appraisal of how the true reflection coefficient varies with position on the structure. Code User's Manual A user manual for the SERATBUILD preprocessor and the FSS/EIGER computational engine and postprocessor is underway, and should be finished shortly. This manual will discuss the use of SERATBUILD to construct the FSS geometries, providing examples and illustrating with sample input and output files. Sample output data from the EIGER postprocessor will also be presented to allow for easy validation checks by the user. REFERENCES [1] J. C. Vardaxoglou, A. Hossainzadeh, and A. Stylianou, "Scattering from two-layer FSS with dissimilar lattice geometries," IEE Proc., Vol. 140, no. 1, pp. 59-61, Feb. 1993. [2] H. Aroudaki, V. Hansen, H.-P. Gemiind, and E. Kreysa, "Analysis of Low-Pass Filter Consisting of Multiple Stacked FSS's of Different Periodicities with Applications in the Submillimeter Radioastronomy," IEEE Trans. Antennas Propagat., vol. AP-43, pp. 1486-1491, Dec. 1995. 30

Hybrid FEM Code with the Fast Multipole Method for Curved SERAT Panels Kubilay SERTEL, RadLab, UM. In the simulation of doubly-curved SERAT problems, the periodicity properties of the geometry cannot be utilized to reduce the size of the problem as is done in the planar periodic SERAT problems. The size of the resulting matrix equation by the discretization of the reallistic doubly curved SERAT geometry is usually large. The conventional FE-BI can no longer be utilized due to its large memory requirement and computational complexities. 1 Fast Algorithms The current FE-BI version of FEMA-PRISM code is being modified to be capable of solving larger problems involving curved geometries in a shorter time period using less memory. To be able to solve large problems on a given hardware in a shorter time, the large storage and computational complexities in the current version version should be decreased. This can be accomplished by employing the Fast Multipole Method (FMM) in the boundary-integral formulation of the problem. Direct implemantation of the FE-BI method requires the computation of NBA double surface integrals appearing as the elements of the resulting BI matrix. Iterative solvers dominantly require 0(NB,) operations per iteration. The memory requirement of the direct formulation is also O(NB2). This large order for storage limits the size of the problem that can be solved on a given hardware, and the high operation cost poses a limit to the size of problems that can be solved in a practically acceptable period of time. On the other hand the FMM reduces the computational complexity and the memory requirement complexity to O(NB'M) without sacrificing solution accuracy. For the simulation of large, doubly-curved, non-periodic (finite) SERAT configurations, the FMM is being implemented and tested in terms of solution time, memory requirement and solution accuracy. At this moment, we have completed a first version of the PRISM code incorporated with the FMM and obtained several results for validation and comparison with the direct 31

method (without the FMM). We are examining the speed-ups attained with this method, the reduction in the memory demand and the practical CPU requirements for typical SERAT problems. Generalization of the code to be capable of curved FSS modeling is straightforward. The lower storage and computational complexities of the FMM code makes it very attractive in the solution of large realistic problems such as large antenna arrays on multiple FSS layers. The next section outlines the formulations used in the FMM and the following section present some results and comparisons between the FMM and the convensional method in terms of execution times. 32

2 FMM Implementation In FEM formuations the weighted error in the computational volume V of Figure 1 is forced to be zero in an average sense. V x (V x E) - k2rE] Wjdv = -iZw~o J Wjdv (1) Jv [ Ir /Jv After a series of algebraic manupulations Eq. (1) can be converted into ( VxE) V x Wj) dv-k2 (E) Wdv -i o j (H x j). nds = -iwo J Wjdv. (2) The unknowns E field is expanded in terms of known basis functions as E = E iN/ (3) i where Xi's are the edge unknowns (average tangential electric field along the i'th edge). Hence Eq. (1) is converted into a matrix equation provided that H appering in the surface integral is expressed in terms of E. Making use of the surface equivalence theorem and after some algebraic manipulations, this surface integral can be converted into a form Zi = - 2 fds [, ds'Go(r,r')(Ni x )] (W, x i) (4) + V2 fds [, ds'Go(r, r')V' (Ni x )] V (Wj x ). (5) By defining a new set of testing and basis functions, this integral can be converted into the form r1 r') r')V] eikR r ji 2dstj(r)., ds' [ji(r) + 2V j() R. (6) When the basis functions of Tayfun and Volakis[l] are used in the FE formulations, the new set of basis and testing functions (ji's and tj's) appearing in Eqn. (6) are exactly the basis functions of Rao, Wilton and Glisson[2] defined for the approximation of the surface currents. 33

Boundary S FSS Layers Cavity V Figure 1: The geometry layout of the problem. The FMM is developed using two elementary identities. The first is the expansion of the scalar Green's function appearing in Eq. (6) as iklr+dl oo A - = ik (-1)'(2/1+ 1)j (kd) h (kr) P,(d r), (7) [r + d I =0 which is a form of Gegenbauer's addition theorem[3]. Here ji is the spherical Bessel function, h(l) is the spherical Hankel function of the first kind, Pi is the Legendre polynolmial, and d < r is the condition for the validy of the expansion. In the FMM formulations of scattering problems, where the source point is denoted by x' and the observation point by x, r will be chosen to be close to x -- x' so that d will be small as depicted in Fig. 2. The second identity is the expansion of jlPt product appearing in Eq. (7) as a sum of propagating plane waves: 47ri'jl(kd)Pl(d. ) = f d2keiP'(k r). (8) The Green's function in Eq. (7) can be rewritten using Eq. (8) as e iklr+d 2 id A oo eikr+d i d2ked (21 + i)hl)(kr)P( r), (9) [r + di 47r =o 34

x X r Figure 2: The basic geometry illustrating the relationship between x,x', r, and d. where the orders of summation and integration are interchanged. The idea of the FMM is that the function L TL(kr, k. r) -= (21 + 1)h(1)(kr)P(k ) (10) 1=0 can be computed for various values of kr which is independent of kd[4]. The series is truncated at the Lth term in numerical practice. The number of terms kept, L + 1, depends on the maximum allowed value of kd, as well as the desired accuracy. The choice of L will be mentioned later. Using Eq. 9, Eq. (7) becomes eiklr+dl ik fd2 ik.d I+ dl 4 de TL(kr, k r. (11) The direct path from a source point to the field point can be decomposed into three parts as in Fig. 3, where rji = rjm + rmm - rim,. (12) The unknowns are geometrically grouped into M clusters. The idea to be noted is that the same path will be used for all source point in cluster m' to translate their field to all observation points in cluster m. Equation (9) can be rewritten as i I d 2keik(rjrim)TL(krmm, k * rmm') (13) 7 ji 47r 35

source r.. field *- -- - - center m rImmi center m Figure 3: The geometry construction used in FMM formulations, illustrating the relation between source point, field point and the group centers. and the Green's function becomes G(rj,ri) -= I- VV — r d k [I - -VV'] ek'(r-*mrm)TL(krmm,, k irmm') = d2k [- kk] eik(r-rim)TL (krmm k. rm'). (14) Using the above equations, a matrix element as in Eq. (6) is approximated by Zmn 2 dstm(r) ds' [jn(r/) + -V' jn(r')V] R 2 Jd2:kVfmj(k) TL(krmm, k r'mmI)V;m-(k) (15) where Vsmli(k) = ds eikim [I-kk].ji(rim,), V/fj(k) = f dseik'm [I-kk] tj(rjm) (16) are the Fourier transforms of the basis and testing functions, respectively, and the superscript * denotes complex conjugation. The memory required for the FMM can be considered in two parts, the sparse-matrix storage and the FMM elements' storage. The storage of the sparse Z' near-field matrix requires O(N2/M) memory locations where M 36

is the number of clusters generated. The FMM aggregations need O(KN) memory locations, and the FMM translations need O(KLM2) memory locations. Hence the total memory storage needed is O(N2/M) + O(KN) + O(KLM2). Using the proportionalities K ca L2, D2 oc N/M, and L oc D, this expression can be simplified to C1(N2/M) + C2(NM N/M), where C1 and C2 are machine- and implementation-dependent constants. The coefficient C2 is so small compared to C1 for all problem sizes that can be solved with the FMM that the memory required is dominated by the O(N2/M) term. The computational complexity of the FMM can be determined by counting the number of floating-point operations required at each step of the algorithm. The aggregation step requires MIN/M = KIN operations. The translation step requires KM2 operations with the precomputed KM2 values of the translation function given in Eq. (10). The disaggregations require MKN/M = KN operations, and finally the sparse matrix-vector product requires N2/M operations. Using the proportionalities IK c L2, D2 oc N/M, and L oc D, the total cost of the matrix-vector product is found as O(NM) + O(N2/M). 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). Both the operation cost and the memory requirement of the FMM is less than those of standart MoM formulation for problem sizes larger than 1000, which makes the FMM more suitable for the solution of large problems. 37

3 Examples and Comparisons Y feed o x d Figure 4: The patch antenna used as the elements the arrays. The antenna is a 3.5cm. x2.5cm. patch antenna residing in a 4.5cm. x3.5cm. x 16cm. cavity. The cavity is filled with a dielectric with cr = 2.17. The feed is located at (0.75cm,0.25cm) away from the origin as depicted. The antenna is operated at 2.64GHz. The FE-BI discretization results in 104(56 for BI) unknowns. 90 1 120 60 150/. -.. ' \ 30 180......... I 0 210\ -: /330 240 300 270 Figure 5: Normalized E- and H-plane radiation patterns of the single patch antenna. (-) E-plane, ( —) H-plane. 38

nn i On O Li Figure 6: 3x3 element patch antenna array constructed from the patch antenna of Figure 3. The elements are uniformly excited.The FE-BI discretization results in 1120(600 for BI) unknowns. 90 1 150/ I-.-...30.. 30 18|...... 0 210 \.. /330 240 300 270 Figure 7: Normalized E- and H-plane radiation patterns of the 3 x3 uniformly excited patch antenna array. (-) E-plane, ( —) H-plane. 39

LII I LIZ LI LII LIZ LI] LI I LI I LI L] I IZ0 LII] [IZ LIZ LI L]LI] I LII]LI LI] LIZ LIIZ LIZ LIII] LI LI LI]I Figure 8: 6x6 element patch antenna array constructed from the patch antenna of Figure 3. The elements are uniformly excited.The FE-BI discretization results in 4699(2496 for BI) unknowns. 90 1 21 3300 1 240. 300 180......................... 0 270 Figure 9: Normalized E- and H-plane radiation patterns of the 6 x 6 uniformly excited patch antenna array. (-) E-plane, ( —) H-plane. 40

0 a. 0 1000 1500 Number of BI Unknowns 2500 Figure 10: Comparison of the CPU times consumed per one iteration in the iterative solution of the resulting matrix equation. CPU Time Consumed for BI Matrix Construction 0 0 E Q. 0 1000 1500 Number of BI Unknowns 2500 Figure 11: Comparison of the CPU times consumed for the filling of the full BI-matrix. 41

References [1] T. Ozdemir and J. L. Volakis, "Triangular Prisms for Edge-Based vector finite elements analysis" IEEE Trans.Antennas Propagat., pp. 788-797, May 1997 [2] S. M. Rao, D. R. Wilton, and A. W. Glisson, "Electromagnetic Scattering by Surfaces of Arbitrary Shape", IEEE Trans.Antennas Propagat., pp. 409-418, May 1982 [3] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical Functions", National Bureau of Standards, 1972 [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 42