031307-5-T MODELING OF CONFORMAL ANTENNAS ON DOUBLY CURVED PLATFORMS AND THEIR INTERACTIONS WVITH AIRCRAFT PLATFORMS Annual Progress Report T. Ozdemir, J. Gong, S. Legault J. Volakis, T. Senior J. Berrie, R. Kipp, H. Wang Mission Research Corp. Rome Laboratories/ERPT 3975 Research Blvd. US Air Force Dayton, Ohio 45340 Griffiss AFB, New York 13441 October 1995

PROJECT INFORMATION PROJECI TITLE: Electromagnetic Modeling of Ultra Low Sidelobe Antennas in Presence of Complex Platforms REPORT TITLE: Finite Element Analysis of Conformal Antennas Using Prismatic Elements REPORT No.: 031307-5-T DATE: October 1995 SPONSOR: J. Berrie Mission Research Corp. (Prime: Rome Laboratories) 3975 Research Blvd Dayton, Ohio 45340 phone: (513) 429 9261 T. Blocker Rome Laboratory/ERPT U.S. Air Force Griffiss AFB New York 13441 SPONSOR CONTRACT No.: SC-0036-93-0002 (F30602-93-C-0036) U-M PRINCIPAL INVESTIGATOR: John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave Ann Arbor, MI 48109-2122 Phone: (313) 764-0500 FAX: (313) 747-2106 CONTRIBUTORS TO THIS REPORT: T. Ozdemir, J. Gong S. Legault, M. Nurnberger, J. Volakis, T. Senior R. Kipp and J. Berrie 1

TABLE OF CONTENTS Executive Summary Triangular prisms for edge-based vector finite element analysis of conformal antennas Describes the theory and validation of the new code FEMA_PRISM A hybridization of finite element and high frequency methods for pattern prediction of antennas on aircraft structures This document gives an overview of the finite element codes (FEMA_CYL, FEMA_TETRA and FEMA_PRISM) and the asymptotic code APATCH. The method of hybridization is discussed along and measured results are presented for validating this hybridization. Efficient finite element simulation of Slot antennas using prismatic elements Describes the performance of the prismatic code for modeling narrow slot antennas and bandpass radomes consisting of multiple slot and dielectric layers. Design of planar absorbing layers for domain truncation in FEM applications Provides the first design criteria for the new uniaxial artificial absorber. These were concluded after an extensive study and although based on two dimensional numerical computations, the design trends are expected to hold for three dimensions as well. 2

EXECUTIVE SUMMARY Project Goal The goal of this project is to develop techniques and codes for modeling antennas and arrays on doubly conformal platforms. Of particular interest is the integration of the developed antenna software with a high frequency code for computing the antenna interactions with the substructures. The development of the high frequency companion code is to be carried out by DEMACO, Champaign, IL. At U-M, the plan is to first develop a code for the characterization of antennas on a circular platform. The second year and subsequent effort will then be devoted on developing a similar code for printed antennas on doubly curved platforms and its integration with the companion high frequency code. Project Status This report describes most of our activity during the second year of the ultralow sidelobe project. During the first year we completed and validated the antenna code for cylindrical platforms. This code, referred to as FEMA_CYL was written by Leo Kempel as part of his graduate work. By all measures, this code was a great success and generated accurate results for scattering and radiation by uncoated and coated patch antennas and arrays on cylindrical platforms. The theory and validation of this code is described in the report 031172-2-T. User manuals (reports 031307-1-T and 031307-2-T) were also written at the end of the first year. The code FEMA_CYL was also interfaced with APATCH to obtain the results given in a section of this report. The emphasis during the past year has been on the development of a finite element code for the characterization of antennas on doubly curved platforms. This code, referred to as FEMA_PRISM, embraces several features which are particularly attractive for printed antenna analysis: * Employs prismatic elements rather than tetrahedrals. These elements lead to substantial simplification of the mesh generation task without compromising geometrical adaptability as was the case with bricks. * Narrow slot antennas and details in printed elements can be modeled with great ease. 3

*Mesh truncation is accomplished with the use of unique artificial absorbers. Unlike ABCs, artificial absorbers can be applied conformal to the surface and consequently antennas on doubly curved platforms can be analyzed with the same ease as those on planar and cylindrical platforms and without a need to employ a Green's function. Recent community emphasis on artificial absorbers for truncating finite element meshes confirms our early choice to follow this path. After 15 months or so under development, the FEMA_PRISM code has been completed and has been validated for input impedance and radiation pattern calculations. The list of validations include: * Eigenvalue computations for material loaded cavities * Rectangular patches on planar and cylindrical platforms * Circular patches on planar and spherical platforms ~ Patches on conical surfaces * Patches on ogival (doubly curved) surfaces * Slot antennas and radomes The latter validation is based on new measured data collected recently at the Naval Weapons Center (by R. Sliva and H.T. Wang) for the purpose of validating the FEMA_PRISM code. The model used for measurement was specified by the University of Michigan and was manufactured at the Naval Weapons Center in China, Lake. As of this moment, FEMA_PRISM is a stand-alone code and can be easily developed into a user oriented code since it is already interfaced with automated mesh generation routines. Our plan over the next few months is as follows: * Write users manual * Modify code to employ the new unixial artificial absorber for mesh truncation * Modify code to compute equivalent currents on a raised surface (not necessary conformal to the platform). These set of currents are intended to facilitate the interface with APATCH 4

* Complete interface with APATCH. * Develop an improved hybridization of antenna with platform. This is essential for computing low sidelobe patterns. * Check and validate improved hybridization with APATCH. This will involve the computation and integration of physical optics and asymptotic currents due to antenna radiation * Improve users interface * Add/Develop mesh generator libraries for additional antenna geometries * Incorporate improved feed modeling (aperture coupled, coax, etc. ) 5

Triangular prisms for edge-based vector finite element analysis of conformal antennas T. Ozdemir and J. L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, Michigan 48109-2122 October 16, 1995 Abstract This paper deals with the derivation and validation of edge-based shape functions for the distorted triangular prism. Although the tetrahedron is often the element of choice for volume tessellation, mesh generation using tetrahedra is cumbersome and CPU intensive. On the other hand, the distorted triangular prism allows for meshes which are unstructured in two dimensions and structured in the third dimension. This leads to substantial simplifications in the meshing algorithm and many printed (conformal) antenna and microwave circuit geometries can be easily tessellated using such a mesh. The new edge-based shape functions presented in this paper are validated by computing the eigenvalues of three different cavities (rectangular, cylindrical and pie-shell). A number of examples are given to illustrate the effectiveness of the new elements in analyzing conformal antennas on doubly-curved platforms. 6

1 Introduction The brick and tetrahedron are popular elements in finite element analysis of electromagnetic problems. The first is indeed attractive because of its simplicity in constructing volume meshes whereas the tetrahedron is a highly adaptable, fail-safe element. It is often the element of choice for three dimensional (3D) meshing but requires sophisticated and CPU-intensive meshing packages. The distorted prism (see Figure 1) is another volume element which provides a compromise between the adaptibility of the tetrahedron and the simplicity of the brick. Basically, the distorted prism allows for unstructured meshing (free-meshing) on a surface and structured meshing in the third dimension. An approach for growing prismatic meshes is illustrated in Figure 2 and most volumetric regions in antenna and microwave circuit analysis can be readily tessellated using such a mesh. As seen, once the surface mesh at the different surface levels is constructed, the prismatic mesh can be generated by simply connecting the nodes between the adjacent mesh surfaces. This avoids use of CPU-intensive volumetric meshing packages and in many cases the user/analyst can construct the mesh without even resorting to a surface meshing package. Examples include printed circuits and antennas on planar platforms. Moreover, because of their triangular cross-section, the prisms overcome modeling difficulties associated with bricks at corners formed by planes or edges intersecting at small angles. A special case of the distorted prism is the right prism which is characterized by the right angles formed between the vertical arms and the triangular faces [1]. The top and bottom faces of the right prism are necessarily parallel and equal, restricting them to a limited range of applications, namely, geometries with planar surfaces. In contrast, the distorted triangular prism is almost as adaptable as the tetrahedron with the exception of cone-tips which are not likely to occur in printed antenna and microwave circuit configurations. In this paper, we introduce edge-based basis functions [2]-[4] for the most general distorted prisms. These prisms have non-parallel triangular faces and each of their three vertical edges can be arbitrarily oriented. In the following, we first present the edge-based shape functions and then proceed with the derivation of the finite element matrix. Eigenvalue computations as well as antenna analysis results are included and used to validate the new edge-based basis functions for prismatic elements. 7

ZA 3 2 (II y x Figure 1: The distorted triangular prism shown with the directions of the edge vectors. Surface normal \S urface geometric Surface:~~ ~ n~fidelity Distorted triangular prisms stacked up Figure 2: Illustration of an approach for constructing prismatic meshes (unstructured surface mesh and structured mesh along the third dimension) 8

Perspective view Top view \. containing the interior point, K... " i- /' "". " Interior point (x,y,z) \/, - /'s. Triangular cross-section containing the interior point Figure 3: Variation of the triangular cross-section along the length of the prism 2 Vector Edge-based Basis functions Consider the distorted prism shown in Figure 3. The prism's top and bottom triangular faces are not necessarily parallel to each other and the three vertical arms are not perpendicular to the triangular faces. On our way to specifying a set of shape functions for the prism, we proceed with the identification of a triangular cross-section of the prism which can be uniquely defined given a point (x, y, z) interior to the prism. The cross-section contains this point as illustrated in Figure 3. A way to specify the nodes of such a triangular cross-section is to use the parametric representation ri = rib + s(rit - rib), i = 1,2,3, (1) where ri are the nodal position vectors of the triangular cross-section (see Figure 4(a)) whereas rit and rib are associated with the top and bottom triangular faces, respectively, and 0 < s < 1. Thus, for s = 0 ri specify the nodes of the bottom triangle and for s = 1 ri reduce to the nodes of the top triangle. Since s assumes values between these two limits, the cross-section sweeps the entire volume of the prism. Clearly, each (x, y, z) point belongs to one of the triangles specified by ri and consequently there is a unique correspondence between s and a given point in the prism's volume. 9

< —0 3 3^ d,.; _-3, Z 1.. - - (x, y,,z) ( * -* (xyi ^ 1, ^ —-1 _ $ {_L_ (x.y,zz ~' 0 —4* (xyz- * hp (X 1b' 1' Z Id * * (X. Y, Z 2. t 22 (a) (b) Figure 4: (a) Nodal coordinates, (b) triangular cross-section with the local coordinates 4 and i1. To determine s in terms of (x,y,z), we first note that (1) represents nine equations with ten unknowns (s and the coordinates of the three nodal points). Thus, additional equations are needed and these are the equations defining the plane coinciding with the triangle, viz. Xi Y" Z' xi^+ + z+1=0_O. i = 0 1,2,3 (2) a b c (a + )+ 1 = (3) a b c in which the constants a, b, and c correspond to the locations at which the plane inn tsects the xy, and z axes, respectively. Substituting in (2) the values of x, yi and zt obtained from (1) yields three equations involving four unknowns a, a, c and s. Solving these three linear equations for -, l and 1 polynomial — i+ + -- +1 03 i-1 2 a (2) a3s3 + as2 +as + a= (4) whose real root is the desired value of s. The expressions for the coefficients ap can be given in closed form in terms of x, y, and z, but are rather lengthy and bear no significance on the rest of the analysis. Therefore, they have been omitted but can easily be derived as outlined above. The appropriate solution of (4) in terms of ao,l,2,3 is [5] s = (s1 + s2) -( (5) 3a3 10

where S1 = [r + (q3 + r2)] (6) 82 = [r-(q + r2)]3 (7) r = 6(a2 - 3aa3) - -(a2/a3)3 (8) 1 12 q = a - -a2 (9) 3 9 ( and this completely specifies s in terms of x, y, and z. We next proceed with the derivation of the basis functions. We choose to represent the field variation across the triangular cross-section (defined by s) using the Whitney-1 form [6]. A simple linear variation will be assumed along the length of the prism. Specifically, the vector basis functions for the top triangle edges can be expressed as N1 = d1 (L2VL3- L3VL2) s N2 = d2 (L3VL1-L1VL3 ) (10) N3 = d3(L1VL2-L2VL1 s and correspondingly those for the bottom triangle edges will be M1 - dl (L2VL3- L3VL2)(1 -) M2 = d2(L3VL1 -L VL3)(1 -s) (11) M3 = d3 ( LVL2-L2VL ) (1 -). The subscripts in these expressions identify the edge numbers as shown in Figure 1 and the distance parameters di are equal to the side lengths of the triangular cross-section containing the observation point (see Figure 4(b)). Also, cos a3 sin a3 L2(, r7) -= rh- h 71 (12) h2 h2 cos ao2 sin &2i L3(7, r/) -= ( r1 11

2 3 (a) (b) (c) Figure 5: (a) Vector map of N 2 or M2 (b) Vector map of Kl (c) Variation of ) as a function of x, y, and z. are the usual two dimensional scalar node-based basis functions [7] for the same triangle with ai denoting the interior vertex angles and hi being the node heights from the opposite side. The variables - and r/ represent the local coordinates and are illustrated in Figure 4(b). As required with all edge-based shape functions, Ni ~* e and Mi * ei have unity amplitude on the ith edge whereas Ni ~ ej = Mi ~ ej = 0 for i - j. Their vector character is depicted in Figure 5(a) and as seen they simply "curl" around the node opposite to the edge on which their tangential components become unity. It remains to define the shape functions for the three vertical edges and we chose to express these by the linear representations K1 (, ) (, ) L1 (, V) K2(r, r) - v(, r) L2(r, V) (13) K3(S, V) = v(, v) L3(, r) As before, Li are the node-based shape functions defined in (12) and a pictorial description of K1 is found in Figure 5(b). Of particular importance in (13) is the unit vector v(i, r). It is a linear weighting of the unit vectors vi, v2 and V3 associated with the vertical arms (see Figure 5(c)), and is given 12

by V() = -11 3 A^ (14) This particular choice of vi is oriented parallel to the side faces of the prism when evaluated on those surfaces and minimizes tangential field discontinuity across the side faces. Another choice is v(, r) = Vs (15) which is always pointed normal to the triangular cross-sections of the prism and ensures tangential field continuity across the top and bottom triangular faces. Both choices are equally useful. However, while the one in (14) can be computed using a simple analytical formula, the evaluation of (15) has to be carried out numerically. Therefore, in terms of accuracy, (14) is a more efficient and accurate choice for v. Contrary to the rectangular brick, tetrahedron and right prism, the edgebased vector basis functions for the distorted prism do not ensure tangential field continuity across the faces, nor are they divergenceless. The same is true for the curvilinear bricks [8] and, in the case of the distorted prism, the field discontinuity and divergence increase with surface curvature (see Figure 2). This is certainly an undesirable feature but in most cases (particularly when sampling at 15 or so points per linear wavelength and the radii of curvature are greater than a wavelength), the angular deviation of the vertical arms is quite small. Consequently, for all practical purposes the field discontinuity and divergence are negligible. 3 Eigenvalue Computation In this section, we examine the validity of the presented edge-based functions. Specifically, we consider the eigenvalues of three different cavities using the edge-based distorted prism as the tessellation element. We begin by first deriving the matrix elements following Galerkin's testing. The weighted residuals of the vector wave equation are I I I Nj * (V x V xE - k2E)dV =0, i = 1, 2, 3, (16) S fMi * (V x V x E- k2E)dV =, i= 1,2,3, (17) 13

I y Ki v x v xE - kr)d = u2 = 1,2,3, (18) in which Ni, Mi and Ki comprise the nine edge-based vector basis functions defined in the previous section and E is the electric field vector. The matrix equations are generated by introducing the representation 3 E = E [EiNN(r) + EiMMi(r) + EKKi(r)] (19) i=1 where EiN, EiM and EiK are the expansion coefficients, and correspond to the average amplitudes of the field vector at the ith edge. Substituting (19) into (16)-(18), and invoking the divergence theorem yields the element equations 3 3 Z EjN[NNCij NND E NMC - NN ] + EM[NM - k2NMDij] + j=1 j=l Z EK[NKCij- k2NKDIj] 0 (20) j=1 3 3 EN[MNCij - k2MNDij] + E EM[MMCj3 - k2MMDj] + j=l j=1 3 EK[MKDi - k2MKDij] = 0 (21) j=1 3 3 EjN[KNCij - kiKNDij] + 5 EjM[KMCij - k2KMDij] + j=1 j=l 3 EjK[KKCij - k2KKDij] = 0, i = 1,2,3. (22) j=1 where NNCie = JJ (V x Ni) * (V x Ne)dV (23) NMCie = JJ f (V x Ni) * (V x Me)dV (24) NKCie = JJJ (V x Ni) (V x Ke)dV (25) 14

MMCi, = JJJ (v x Mi) e (V x Me)dV (26) MKCUt = / /J (V x Mi) * (V x Ke)dll (27) KKCit = l (V x Ki) * (V x K()dV (28) NNDit = Ni.NedV (29) NMDie = J| Ni MedV (30) NKDit = ll Ni.KedV (31) MMDje = IJ MiMedV (32) MKDJi = JJ Mi * KedV (33) KKDie = K i KedV. (34) The above integrals must be evaluated numerically and to integrate over the volume, the prism is divided into three regions (see Figure 6) of triangular or quadrilateral cross-section. In the case of the right prism, the integrals have closed form expressions and are given in the Appendix. Upon assembly of (20)-(22) and boundary condition enforcement, we obtain the generalized eigenvalue system [A]{x} = ko[B]{x} (35) in which A = k2 represent the eigenvalues of the problem. The matrices [A] and [B] are real, symmetric and sparse ([B] is also positive definite). Example 1: Rectangular Cavity The first example is the rectangular cavity shown in Figure 7(a). Results based on brick, tetrahedron [9] and prism discretizations are given in Figure 7(b), where they are compared to the exact values. We observe that, overall, the data based on the prism are better than those based on the brick and, at least, as good as those associated with the tetrahedron. Figure 7(a) displays the actual mesh used for both the brick and prism discretizations. Number of degrees of freedom for the prism, brick and tetrahedron computations were 382, 270 and 260, respectively. 15

Region 1 (0 < z < z ) Region 2 ( z< z < Z2) Region 3 ( z< z < z3) Triangular cross-section Quadrilateral cross-section Triangular cross-section 6z -Z3 - - - - - - - -- \/ A' /'' \ / I \ \Z __V Y:-.-^' Y:...;'/ Figure 6: The three regions of integration over the volume of the distorted prism. -—. Prism volume, - integration subregions, subregion cross-section. Icm — ^ Mode k, cm'-1 % Error 0. 5cm Mode(Exact) Prism Brick Tetra. TE l0 5.236 0.73 -1.36 0.44 TM o 7.025 2.32 -2.23 0.70 0.75cm TEo0, 7.531 0.53 -2.58 1.00 TE 0, 7.531 0.64 -3.13 -0.56 (Actual mesh) TM,1 8.179 0.22 -2.09 2.29 (a) (b) Figure 7: (a) Rectangular cavity discretized using right triangular prisms, (b) eigenvalues for the air-filled cavity. 16

Radius = 1cm (Actual mesh) (Actual mesh)_~_~~_~ _~; Mode k, cril % Error _ —- - _._.; i(Exact) (Computed) E r i___ _____^_ ":TMolo 2.405 1.29 II |i j O~ I TE Il 3.640 2.17 I'4.. —--.' -:. TM1loi 3.830 -2.90 ~ i ii r. TMoi i3.955 0.81' I....' =.... TE211 4380 -8.97 (a) (b) Figure 8: (a) Circular cavity discretized using right triangular prisms, (b) eigenvalues for the air-filled cavity Example 2: Circular Cavity The tessellation of the circular cavity (drum) is shown in Figure 8(a). The mesh was constructed by first creating a surface grid of triangles at the base of the drum. The prisms were then generated by growing the mesh along the axis of the cylinder. In this case, the resulting volume element is the right prism and the eigenvalues are given in Figure 8(b) along with the exact. Except for the higher order mode TE211, the remaining eigenvalues were computed to within three percent of the exact. The number of degrees of freedom for this computation was 544. Example 3: Dielectric ring resonator Figure 9 shows a dielectric ring placed symmetrically inside a cylindrical cavity and of interest is the computation of the resonance frequency in the presence of the dielectric ring. The resonance frequency was measured [10] using a loop antenna connected to a directional coupler as shown and the cavity was excited by the same loop antenna. Since maximum coupling to the cavity occurs at resonance, minimum power is returned to the detector at this frequency. For computation, the cavity was modeled using a mesh 17

( Height = 2.2cm Inner diameter = 7.5cm (Height = 13.84cm I Outer diameter = 11cm J 3 Diameter = 15.24 cm ~ r= 13.5 Metal Cavity Dielectric ring / Directional Coupler' --- Measured = 1282 MHz -- l..., Computed = 1257 MHz (Error = 2%) Figure 9: Dielectric Ring Resonator similar to that for the previous example resulting in 1051 degrees of freedom. The computed frequency was 1257 MHz (based on the smallest eigenvalue) and this is within 2% of the measured frequency (1282 MHz). Example 4: Pie-shell Cavity The fourth example is a pie-shell sector as shown in Figure 10(a). It is obtained by bending the rectangular cavity considered earlier and the resulting volume element is the distorted prism with a vertical arm angular deviation of five degrees. The computed and exact eigenvalues for the first five dominant modes are given in Figure 10(b), and these testify to the accuracy of the distorted prisms in modeling curved geometries. The number of degrees of freedom involved for this computation was 382. 4 Antenna results In this section, we will illustrate the effectiveness of the prisms in characterizing conformal antennas on doubly curved platforms. Figure 11 illustrates antenna modeling using the finite element method. It also shows how the 18

(Actual mesh) _Mode k, cm-' % Error:~ —.; —a (Exact) (Computed):'- r -^ "..TM ll0~ 4.693 1.55...-=.. — /..0.75cm TM210 6.009 1.28 TE l 1 6.640 2.71 cm'-....... —, TE1c 211 7.513 -0.27 0.5cm ( TE 01 7.579 2.73 (a) (b) Figure 10: (a) Pie-shell cavity discretized using distorted triangular prisms, (b) eigenvalues for the air-filled cavity mesh is terminated using metal-backed absorber. Note that the absorber works on the principle of matched wave impedances across the interfaces and therefore different absorber sections must be used to absorb surface waves supported by diffent superstrates and substrates. Accuracy of this termination scheme has already been demonstrated for scattering problems [11]. Examples given below establish its efectiveness for radiation problems. First the surface grid is created. Then the volume mesh is grown along surface normals as illustrated in Figure 2. Resulting volume element is the distorted prism. Process of assembling finite element equations is the same as for the eigenvalue computations since here also the termination boundary is entirely metallic. Example 1: Cavity-backed rectangular patch on planar platform Figure 12 shows a rectangular patch backed by a rectangular cavity recessed in a planar ground plane. Input impedance as a function of frequency has been computed. It is in excellent agreement with the more regirous finite element-boundary integral technique(FE-BI). In FE-BI technique, only the cavity volume is modeled by finite element technique and the mesh is terminated at the cavity aperture using the half-space Green's function, making it an accurate technique [12]. Note the dominant mode distribution of the field under the patch. 19

(a) Cavity-backed antenna b Cs(i2.7)Es = Ca= 1-j2.7 ~ % PEC,;\' O'\j /\ (Actual geometry) (FEM Simulation) (b) Microstrip antenna ~c C~a= a= 1j2.7,b Es (I-J2.7) Eb= EC(l- i2.7). 0. 15 Es b.-.2.7 PEC (Actual geometry) (FEM Simulation) Figure 11: Finite element - metal backed absorber technique for antenna modeling 20

Input Impedance vs Frequency I _00 _ —— _-__l__ l —-_l_ —-l-| llcm | 5cm 2cm Er = rl-j2.7 rin aI (a\l ) ~E( (b). t Abso r-ai niterElem ent-BoundaryIntegra 2.0 3.0 4.0 5.0 71 Tl.5cm (G H z c 1.5cm ----- Finite Element-Metal Backed Absorber I |10.08cm Er =2.17 Resonant frequency =3.17 GHz (disagreement < 1% ) Vertically directed electric field magnitude under the patch iesh termination boundary (PEG) Absorber-air interface Cavity wall... -of the patch. 21 of the patch 21dm.

Example 2: Cavity-backed circular patch on planar platform Figure 13 shows a circular patch backed by a circular cavity recessed in a planar ground plane. Again the input impedance computation is compared to the FE-BI solution and the disagreement in resonant frequency is about 1%. Example 3: Microstrip circular patch on sphere Figure 14 shows a microstrip circular patch placed on a sphere. Here the input impedance as a function of frequency is computed and compared to a moment method result [13]. The disagreement in resonance frequency is 1.7%. This is a very good result in terms of the accuracy of the finite element modeling considering that the referance solution only considers the current distribution of the dominant mode TM1l whereas the finite element solution models the entire physics of the problem. Note that the field distribution under the patch resembles the one for the cavity-backed pacth on planar platform implying that the surface curvature causes little difference in the field structure. Hoever, as will be illustrated later, the curvature does change the resonance frequency. Example 4: Microstrip sectoral patch on cone Having established the acccuracy of the new technique, we have chosen this example to illustrate the capabilities of the new technique as an analysis tool for characterizing conformal antennas on doubly curved polatforms. Figure 15 shows a microstrip sectoral patch placed on a cone. Input impedance as a function of frequency is computed. The predicted resonance frequency is 3.115 GHz. Cavity model predicts 3.214 GHz. Finite element result is the more reliable one as the cavity model is a simplistic one. Note that the field distribution under the patch is different from the rectangular patch case given in Example 1. This is because two radiating edges are not of the same length. Example 5: Cavity-backed rectangular patch on ogive This is another display of capability of the new technique to model antennas on doubly-curved platforms. Figure 16 shows the set-up, where a rectangular patch is placed on the aperture of a rectangular cavity recessed 22

Input Impedance vs. Frequency ^y200.,. r! = 1.3cm (0.15 WL at resonance) E 10 I = 2.17cm (0.25 WL at resonance) 15o 0 /I r ~ r3 = 2.6cm (0.3 WL at resonance) ~ 100 _o - / x (a) r4 = 1.3cm (0. 15 WL at resonance) (b) 50 ~ 50 0.8cm 2 3 4 5 \0 Frequency (GHz) -- Finite Element - Boundary Integral (using tetrahedra) -- Finte Element - Metal Backed Absorber (using prisms) Metal E = U = 1 —j 2.7 Resonance Frequency = 3.55 GHz.3c (disagreement = 1%) g. A5Scm Er= 2.4 Vertically directed electric field magnitude under the patch Metal Termination I iS Boundary ad m h t ( r Absorber-air interface C avity w all................ Patch " J'1~ Feed location j maqip Figure 13: Ractangular patch backed by a rectangular cavity recessed in planar platform,: (a) input impedance computations, (b) antenna geometry and mesh termination, (c) radiating sides of the patch. 23

Input Impedance vs Frequency 200 x X O 0 " I\ ^ ".- - c...; (S ide view, X in.z' in | Top view) Feedi (a) E (b) *- b.' Finit (Side vie- A 3GHz.17 2-.14G L,.:i Metallic ground plane ----- Moment Method (TM11 mode) Disagreement in Resonsonant frequency = 1.7% Vetically directed electric field magnitude under the patch Metal termination Pa( (Cross-sectio view) Fe (Top view) Feed location Figure 14: Microstrip circular patch on sphere: (a) Input impedance, (b) geometry, (c) field distribution of dominant mode. 24

Input impedance computation Metallic 200 \' Ground Plane h 014 h = 0.114cm 0 Patch rb =12.9cm r =10cm -100 r 3.0 3.1^ 1 3.2 3.3 4%= 18.340 (GHz) 0= 33.88~ fr= 3.115 GHz y Cavity model estimate =3.214 GHz X Vertically directed electric field strength under the patch Metal Termination Boundary Absorber-air interface Microstrip patch J!l' ~ Feed location SW. An,,,,,,,,.., _........ _.... Figure 15: Microstrip sectoral patch on cone: (a) Input impedance, (b) geometry, (c) field distribution of the dominant mode 25

t6_" 6" 6" (Side view) (Front view) Feed 6cm location A (Top view) Cavity Patch o/ - ^ A1.125cm ~ r = 2.17. Observation nY 1r 3cm path 0.08cm - Polarized Radiated Power ()- Polarized Radiated Power -10 I 0 o 0 -20: \ -0 \ 2' - -30 -30' -40 -40 -50 -50....... -180 -90 0 90 180 -180 -90 0 90 180 Observation Angle d (Deg.) Observation Angle j (Deg.) ------ Measurement Computation Figure 16: Cavity-backed rectangular patch antenna on ogive: (a) Geometry of antenna+platform, (b) geometry of the antenna, (c) radiation pattern 26

in an ogival ground plane. Antenna region has been again modeled as in Example 1. Upon obtaining the tangential electric fields on the cavity aperture, cavity aperture is sealed with metal and equivalent magnetic currents (M = E x n) are placed where the cavity aperture used to be. To produce the results given in Figure 16c, physical Optics approximation has been used, i.e., the magnetic currents have been doubled in strength and allowed to radiate in free-space. Computed radiated power is plotted against the measured data [14]. For 0-polarization, agreement is good. However, for Cpolarization, the computed pattern follows the same trend as the measured pattern but with substential disagreement. This is not surprising since we have not properly accounted for the platform interactions when radiating currents. Results also point out that 0-polarized radiation is not as sensitive to the shape of the platform as the q-polarization. To improve especially the +-polarization results, true shape of the platform must be taken into consideration when radiating magnetic currents. This can be accomplished in two ways: In addition to the direct radiation from the antenna, one must also consider creeping waves which radiate while creeping on the surface of the body. Another way is to use high frequency codes such as APATCH to include platform interactions [15]. Example 6: Effect of curvature on the resonant frequency of patch antennas One of the most important design parameters of antennas is the resonance frequency. For patch antennas, resonance frequency, to large extent, depends on the dimensions of the patch. However, it also depends on the curvature of the patch surface. Figure 17 shows the shift in the resonance frequency as a function of the platform curvature for rectangular as well as circular patches. It is shown that the shift in the resonance frequency is larger when the patch is curved such that the radiating sides get closer to each other. 5 Conclusions The distorted prism is an indispensable tool for the analysis of many doubly curved antenna and microwave geometries. It provides simple volume meshing without compromising accuracy (see the comparison with the tetrahedron in Figure 7(b)). In this paper, we introduced a set of edge-based shape 27

70 i....l 60.. i 50 /P /.:' \ ~.' Radiating (a) E 40 - / \ \ sidesi; 30 / X 20 - 10 1o,-K_ 3. i f2- f3 = 112MHz 3.10 3.12 3.14 3.16 3.18 3.20 f3 f =22MHz Frequency (GHz) / Y. —..~-" /(10/ 0 / / / - / I 0 - 0~~ p ~~~f2- fl = 70 MHz 3.0 3.2 3.4 3.6 3.8 4.0 f3 fl = 150MHz Frequency (GHz) Figure 17: Effect of platform curvature on resonance frequency: (a) rectangular patch, (b) circular patch 28

functions for the distorted prism and used them to generate the element equations. For validation purposes, the eigenvalues of three different cavities (rectangular, cylindrical and pie-shell) were computed and compared to the exact solutions. In almost all cases, the error was less than three percent. Validation has also been carried out for a number of antenna configurations. Eigenvalue and antenna computations together leave no doubt about the effectiveness of the prisms. Results pertaining to antennas on conic and ogival surfaces have been included in order to illustrate the capabilities of the new technique to model antennas on doubly-curved platforms. Appendix For the right prism, closed form expressions of the integrals in (23)-(34) are possible. Referring to Figure Al, we have didt cos /3kn cos 3jm COS /km ENNCu = -( Xjm + h Xkn- Xjn - c hkhn hjhm hkhm COS 3jn 2 C2hldl Xkm + sinfpjk sinfmn ) hj hn 3 hjhkhmhn did ( cos/ kn cos/ljm cos/ km ENMCie = ( Xjm - Xkn + Xjn + c hkhn hjhm hkhm COS /3j 1 c2 h d Xkm + L Sin3jk SinPmn ) hjhn 3X hjhkhmhn ENKCit = hldl dicos pje _ cos /3k ENMKdfe = ( 6 hjhe hkhe EMMCe = ENNCie EMKCie = -ENKCie EKKCr = hldl cos i EKKCii cd- c~ Oi^7 2 hihl 29

did cos!kn COS 3jm COS.km COS km) ENNDit =C ( Xm + Xkn - -- Xkm) 3 (63 hkh hjhm khm h hjhn ENMDie = -ENNDi 2 ENKDie = 0 EMMDie = ENNDie EMKDie = 0 EKKDie = cXie where r 0 ifr-s rs- ar + cs, otherwise hld1 hi Xrs = 2 {WrWs + — [(cota3 - COta2)(lsWr + tIrws) + 2(QsWr +,rWs)] + 2 3 -'[3(cota3 - cota2)(lsir + 7rrs) + 271rTrs(cot22 - cota2cota3 + 12 cot a3) + 6rgS]}, r, = 1,2, 3, the set {i,j, k} takes on the value {1,2,3}, {2,3,1} or {3,1,2}, and the same is true for the set {l, m, n}, W1 = 1 1 =-h 71=0 sin a3 W2=0 2= - -- ~h2 h2 W3 -= 0 cs =h 7/3 =sin a2 w3=0 3= -eh3 s3 References [1] Sacks, Z., S. Mohan, N. Buris and J. F. Lee, "A prism finite element time domain method with automatic mesh generation for solving microwave 30

3 N2 N i2r~~~ ~O- 3 K K N d h, - iK / i.... M2','-.~'M'1 d3 /. ). M, (1') d 3 (2) M3 (a) (b) Figure Al: Right triangular prism: (a) Indexing of nodes and edges, (b) parameters defining the triangular cross-section. cavities," IEEE APS Int. Symp. Digest, Vol. 3, pp. 2084-87, Seattle, Washington, June 19-24, 1994. [2] Nedelec, J. C.,"Mixed finite elements in R3," Numer. Math.,Vol. 35, pp. 315-341, 1980. [3] Bossavit, A.,"A rationale for'edge-elements' in 3-D fields computations," IEEE Trans. Magn., Vol. 24, No. 1, pp. 74-79, January 1988. [4] Webb, J. P.,"Edge elements and what they can do for you," IEEE Trans. Magn., vol. 29, no. 2, pp. 1460-1465, March 1993. [5] Abramowitz, M and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, Inc., New York, 1972. [6] Bossavit, A.,"Whitney forms: a class of finite elements for threedimensional computations in electromagnetism," IEE Proc. Part A, Vol. 135, No. 8, November 1988. [7] Zienkiewics, O. C. and R. L. Taylor, The Finite Element Method, 4th ed., Vol. 1, McGraw-Hill, New York, 1989. 31

[8] Antilla, G. E. and N. G. Alexopoulos, "Scattering from complex threedimensional geometries by a curvilinear hybrid finite-element-integral equation approach," J. Opt. Soc. Ame. A, Vol. 11, No. 4, pp. 1445-57, April 1994. [9] Chatterjee, A., J. M. Jin and J. L. Volakis, "Computation of Cavity Resonances Using Edge-Based Finite Elements," IEEE Trans. Microwave Theory Tech., Vol. 40, No. 11, pp. 2106-2108, Nov. 1992. [10] Muldavin, J. B., A. D. Krisch and M. Skalsey, Personal communication, University of Michigan, 1995. [11] Ozdemir, T. and J. L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Sci., Vol. 29, No. 5, pp. 1255-1263, Sept.-Oct. 1994. [12] Jin, J. M. and J. L. Volakis,"A hybrid finite element method for scattering and radiation from microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas and Propagat., Vol. 39, pp. 1598-1604, November 1991. [13] Tam, W. Y., A. K. Y. Lai and K. M. Luk, "Input impedance of spherical microstrip antenna," IEE Proc.-Microw. Antennas Propag., Vol. 142, No. 3, June 1995. [14] Sliva, R. J., Personal communication, Naval Air Warfare Center, Weapons Division, China Lake, California. [15] DEMACO, Inc., 100 Trade Centre Drive, Champaign IL 61820. 32

A hybridization of finite element and high frequency methods for pattern prediction of antennas on aircraft structuresl T. Ozdemir, M. W. Nurnberger, and J. L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, Michigan 48109-2122 volakis@um.cc.umich.edu R. Kipp DEMACO, Inc. 100 Trade Centre Drive Champaign, IL 61820 kipp@demaco.com J. Berrie Mission Research Corporation 3975 Research Blvd Dayton, Ohio 45430-2108 jberrie@mrcday.com 1 Abstract This paper considers the hybridization of the finite element and high frequency methods for predicting the radiation pattern of printed antennas mounted on aircraft platforms. The finite element method is used to model the cavity-backed antennas whereas the interactions between the radiators and the substructures are treated via a high frequency technique such as the GTD, PO/PTD or SBR. We present comparisons between measurements and calculations along 1 This work was supported in part by the Rome Laboratory. 33

with a qualitative description of the employed finite element and high frequency codes. 2 Introduction Placing antennas on aircraft inevitably introduces distortion (often significant distortion) in the resulting radiation pattern. This distortion has important implications both for on-board communications and remote sensing systems. For example, once the antenna is mounted, one may not get the desired coverage for effective communications. In the case of a radar antenna, pattern side lobes may be generated which reduce capability for isolating and locating targets (see Figure 1). An accurate prediction capability for pattern distortion is thus of high value to system designs. This application is also part of a wider class of military and commercial problems involving antennas mounted on any type of scattering platform, including satellites, buildings, and ground vehicles. One aspect of pattern prediction common to all of the above applications is the requirement for a characterization of fine geometrical details in the presence of large structures. At the small scale, simulating the behavior of the antenna is best done using a rigorous technique such as the finite element method (FEM) [1]-[6] or the method of moments (MOM) [7]-[9]. This may include features of the geometry immediately surrounding the antenna to predict loading, which will affect both the antenna current distribution and the terminal impedance. At the large scale, the platform on which the antenna is mounted often spans tens or hundreds of wavelengths in all three dimensions. It is not therefore computationally practical to apply rigorous methods to such large scatterers. Instead, one can use a high- frequency technique which takes advantage of the large scale geometry. Methodologies in this category include geometrical optics (GO) with geometrical theory of diffraction (GTD) [10] and physical optics (PO) with physical theory of diffraction (PTD) [11]. A total solution to the platform mounted antenna thus calls for a hybrid approach, one that uses both low-frequency and high-frequency methods. 34

This paper presents a straight-forward hybridization of the FEM and two high-frequency codes, one based on uniform theory of diffraction (UTD) and one based on PO and the shooting and bouncing rays (SBR) [12], SBR being a multiple bounce implementation of PO. The FEM code analyzes the antenna in the absence of large-scale structures such as wings, fins, engines, or surface details. The computed pattern or aperture currents for the "free-space" antenna are then used as the source(s) for the high-frequency code. The highfrequency code radiates these sources in the presence of the aircraft to produce the overall radiation pattern. In the following sections, we first give a qualitative description of the particular FEM, UTD, and SBR codes employed as well as the major features of those methodologies. We then present details of the hybridization procedure. Finally, the effectiveness of this hybridization is evaluated by comparing results of analysis with measurements for a patch antenna radiating on a finite cylinder (i.e., a fuselage) in the presence of an attached flat plate (i.e., a wing), as shown in Figure 2. 3 Finite Element Code Description Several finite element codes have been developed at the Univ. of Michigan for the analysis and design of printed antenna configurations. Typically, the printed antenna configuration is assumed to be recessed in some metallic or coated platform and the various codes differ in the element used for the tessellation of the antenna, the type of platform assumed in the analysis (planar, cylindrical or doubly curved) and the closure condition employed for terminating the finite element mesh as illustrated in Figure 3 (cylindrical and doubly curved platform). The following codes are available for antenna radiation and scattering analysis. In all of the above codes, the FEM is used to model the cavity region and the boundary integral or absorbing boundary conditions (ABC) or artificial absorbers are employed for truncating the mesh. FEMA-CYL: Code employs cylindrical shell elements (see Fig35

ure 4) and is therefore best suited for cavity-backed rectangular patch antennas recessed in metallic or coated cylindrical platforms(see Figure 3). Its main advantage is simplicity of mesh generation done automatically and quickly after the specification of the few geometrical parameters. FEMA-PRISM: Code employs distorted prisms for modeling the computational domain. As a result, doubly curved platforms (see Figure 5) can be simulated without much difficulty in creating the volume mesh. Once the mesh is created on the surface of the interface containing the patches, the volume mesh can be generated by growing the prismatic elements above and below the surface of the printed antenna. Any patch shape can be modeled using prisms since the faces of the prismatic elements bordering the patches are triangles. The mesh is terminated using artificial absorbers. FEMA-TETRA: Code employs the more general tetrahedral elements for modeling the cavity volume (see Figure 6). Use of this element makes the code adaptable to any antenna geometry recessed in a cavity residing in planar platforms. However, this necessitates use of sophisticated finite element mesh generation packages such as SDRC I-DEAS and PATRAN [13]-[14]. The theoretical background of FEMA-CYL is described in [15] and [16], and the code has been validated for antenna and scattering applications. An example calculation is illustrated in Figure 7. The measured data given in Figure 7 were obtained at NASA-Langley [17] with the patch antenna placed on the cylindrical section of the ogive+cylinder structure as illustrated in the Figure. For this calculation, the cylinder Green's function of the second kind was employed for truncating the finite element mesh to generate a combined finite element-boundary integral system. The latter is solved using the biconjugate gradient [18] or QMR [19],[20] iterative solver and the FFT is used to speed-up the matrix-vector product evaluation in the solver [21]. Consequently, in spite of the fact that the boundary integral matrix is fully populated, the code's computational and memory requirements remain at O(N). The system solution gener36

ates the fields within and on the surface of the antenna cavity. In all codes, for the radiation and scattering pattern calculations only the surface electric fields are used to generate equivalent magnetic currents. These are subsequently substituted in the radiation integral with the appropriate platform Green's function to generate the antenna parameters. For antenna excitation, a probe or a coaxial cable feed model is employed and in the case of scattering the excitation becomes the aperture fields due to the incoming field. When a probe feed model is used, the input impedance is the line integral of the electric field along the length of the probe divided by the magnitude of the probe current. When a coaxial feed is employed, the excitation is introduced by setting the fields coinciding with the edge elements bordering the opening of the coax cable equal to AV/AL, where AV is the given potential between the outer and inner surface of the cable and AL b- a, where b and a correspond to the outer and inner radius of the cable, respectively. The center conductor is modeled by setting to zero the field unknowns that are associated with the edges coinciding with the conductor. The improved results of this feed modeling approach are shown in Figure 8 and the details of the modeling are described in [22]. The codes can be used to include the effect of superstrate materials (see Figure 3) which may extend over the antenna platform. To avoid use of complicated Green's functions, in this case the mesh is extended a fraction of a wavelength over the cavity's aperture and an absorbing boundary condition (ABC) or an artificial absorber is employed for truncating the FEM mesh [11]. With this type of mesh truncation, the entire system is sparse and thus the memory and CPU requirement are again O(N) without a need to make use of the FFT. We have found that placement of the mesh at a distance of 0.3 wavelengths away from the coatings surface is sufficient to model the antenna with reasonable accuracy as illustrated in Figure 7(b) and in Figure 9. The data in the latter were obtained using FEMAPRISM and the data in the former was computed by introducing equivalent electric and magnetic currents placed at the surface of the dielectric coating. 37

Results generated by using the FEM codes refered to above will not, of course, include effects due to interactions between the antenna and the substructures residing on the same platform. To include these, the FEM codes are interfaced with the high frequency analysis packages described next. 4 High Frequency Code Descriptions GTD/UTD Code: The hybridization of GTD and low frequency codes such as the moment method were considered in the early 1980s [23],[24] for scattering and antenna radiation analysis. Typically, the GTD was employed to account for diffraction contributions from edges and corners near the radiating elements and a review of some of this work was given more recently in [25]. In the case of scattering, hybrid GTD and moment method codes have been used to model small details on larger structures. Specific applications include uses of the moment method to model a small wire or a crack in the presence of a large platform which is associated with several scattering centers caused by geometrical optics reflections and diffractions from edges and corners. For the most part these hybrid codes have been restricted to specialized implementations where the scattering substructure is hard-wired into the code and is characterized by a few scattering centers. Also, for radiation analysis, available GTD codes have been restricted to include a small class of antenna elements such as narrow slots, monopoles and infinitesimal current elements. A well developed UTD code, referred to as NEC-BSC was developed by Marhefka and Burnside [26] at the Ohio State University and was used in conjunction with the FEM codes for printed antenna pattern calculations on simple aircraft platforms. The UTD code computes the interactions between current elements and the aircraft substructure using ray diffraction methods. In the employed UTD code the fuselage is modeled as a truncated cylinder or an ellipsoid and the attached wing and fins are modeled as combination of flat plates. 38

The overall radiation pattern of the specified current element is then computed by adding the direct, geometrical optics, and diffraction (including surface waves) contribution from all possible primary and secondary sources which arrive at the observer. Shadowing is also performed as part of the pattern calculation. In our implementation, the UTD code is supplied with the location and amplitude of the magnetic surface current elements as computed by the FEM code in the absence of the fuselage appendages. Consequently, these currents include all coupling effects internal to the antenna cavity and among the antenna elements. However, they do not include coupling effects due to substructure contributions which return back to the antenna aperture. One way to include these effects is by re-executing the FEM code with the secondary external diffraction contribution as part of the excitation in addition to the feed. Alternatively, reciprocity can be used to first compute the secondary fields arriving at the antenna aperture due to plane wave incidence before executing the FEM code. APATCH Code: This code, developed at DEMACO, Inc., computes the radiation pattern of antennas mounted on electrically large scattering platforms. Like its predecessor RCS code, XPATCH [27][28], APATCH uses 3-D ray tracing on CAD models composed of triangular facets to implement a PO/SBR [29] scattering solution. There are several important differences between this approach and a GO/GTD approach, which can be explained with the aid of Figure 10. First, the PO/SBR approach in APATCH starts by adding the incident field at all observation points, whether blocked from the source or not, and then adding the computed scattered field to produce the total field, which is the net radiation pattern. In GO/GTD, the observed fields start at zero and only receive a direct contribution from the source if the path is unblocked. To rephrase it, APATCH generates field shadowing and specular reinforcement through wave interference, whereas GO/GTD generates these effects through ray optics. Second, in the PO/SBR approach, many rays are launched in all directions, and each ray contributes to each observation point by integrating the ray tube fields at the surface of the scatterer. The ray trace is independent of the observation points. In 39

GO/GTD, ray paths are found which directly or indirectly link the source to each observation point according to the rules of specular reflection and Keller cone diffraction. Hence, only a handful of rays comprise the computed field at each observation point. The practical consequence of this difference is that GO/GTD gives better solutions for geometries involving several canonical shapes, and will compute those solutions relatively quickly. The PO/SBR approach is better suited for complex scattering geometries described in CAD files and not readily or accurately reducible to canonical shapes and has correspondingly longer run times. Figure 10 provides a 2-D illustration of the differences between GO/GTD and PO/SBR for a point source and two observation points in the presence of a curved surface and a plate. In the GO/GTD formulation, only observation point # 1 gets a direct contribution; observation point # 2 is blocked. However, both observation points get a reflection contribution. For clarity, we leave out of the figure possible multi-bounce contributions between the curved surface and the plate, as well as GTD contributions from the two edges of the plate. In the PO/SBR formulation, the incident field is added at both observation points. Although many rays would be launched and traced until they escaped, only one is shown in Figure 10. The dashed tracks labelled PO and SBR indicate field/Green's function convolution integrals over the scattering surface and evaluated at the observation points. Note how the PO contribution from the first bounce interferes with the incident field at observation point # 2 to produce partial cancellation. The first bounce contribution in APATCH is a PO computation. A test ray is launched from the center of each facet back toward the point source. If the ray path to the source is not blocked by other surfaces, then the facet is directly illuminated, and the PO fields are integrated to yield the scattered field. Otherwise, the facet is shadowed, and no first bounce contribution is computed. Facets which are electrically large are first systematically subdivided into sub-facets, as shown in Figure 11. This addresses the problem of partially illuminated facets. Also, by setting the maximum sub40

facet size to a fraction of a wavelength (i.e. A/5), one can assume, for the purposes of far-field scattering, that the incident field local to each sub-facet is constant in magnitude and linear in phase. In this way, APATCH achieves a uniform first bounce ray density over the surface of the scatterer. Multiple bounce contributions are computed in APATCH using SBR. This part of the solution starts by launching ray tubes in the specular reflection direction from each illuminated sub-facet. In Figure 11, the shadowed sub-facets from the first bounce blockage check are blackened. Each ray is then traced, with specular reflection at each bounce surface, until either it has escaped or it has reached a set bounce limit, say 10 bounces. The ray tubes diverge spatially as they propagate and the ray tube fields are updated according to the Fresnel reflection coefficients of the surface. These coefficients are in turn based on the material properties of the surface coatings. At each bounce point, the ray tube is projected on to the surface and the tube fields are analytically integrated over this surface with free-space Green's function for evaluation at all observation points. The APATCH code accepts as input either the antenna aperture magnetic current elements computed by the FEM code or the corresponding far-field radiation pattern. The former is more accurate, while the latter is less computationally expensive and may be adequate if the major scattering fixtures are sufficiently far away. If the source information is given by the magnetic currents, the code generates rays with a short-dipole weighting. Further, scattering surfaces near these sources will be illuminated according to the near field representation of the short dipoles, including radial components of the dipole field. Each current source is treated individually according to superposition, and it is therefore appropriate to group nearby sources for CPU time reduction. An example of antenna analysis by APATCH is illustrated in Figure 12, which shows a Cassegrain antenna with a rectangular guide feed. Each reflector consists of 3600 facets treated as perfect electrical conductors. The rectangular guide is not part of the CAD model. Rather, a vertically polarized parametric cosO source is placed at the mouth of the feed with its main beam in the x-direction and no radiation in the x < 0 half-space. The setting for power q in the parametric pattern is 11.7 in both the E- and H-planes. In the first bounce 41

blockage check, APATCH attempts 20,252 sub-facets and finds that 11,880 are blocked. These are the sub-facets on the Cassegrain primary reflector, which is not directly illuminated. Of these, 8,372 rays go on to illuminate the secondary reflector. For this problem, the multi-bounce contribution is essential since that is the mechanism for including the effect of the primary reflector. The 8,372 rays incident rays reflect off the secondary, and all but 512 of these rays go on to hit the primary. The ray tubes projected onto the primary and secondary reflectors are integrated for all observation angles and added to the free-space incident field from the feed to give the radiation pattern of the Cassegrain antenna. A comparison of Apatch PO/SBR results with measurement is given in Figure 13. 5 Comparison of Measurements and Calculations To examine the effectiveness of the proposed hybridization of finite element and high frequency codes, the cylinder+wing structure shown in Figure 2 was constructed by Naval Air Warfare Center Weapons Division (China Lake, CA)2 and the University of Michigan. The patch antenna geometry is shown in Figure 14 and was placed on the cylinder's surface at a = 28.7~, 45~, and 90~ (see Figure 2). The radiation patterns of the patch antenna in the presence of the wing were obtained at the facilities of Mission Research Corporation (Dayton, Ohio) and are shown in Figure 15. It is seen that the patterns at a = 28.7~ and a = 45~ are affected quite substantially by the presence of the plate since the plate is visible to the antenna. However, when the patch is on the cylinder's top (a = 90~), the plate's effect diminishes and the overall pattern is nearly identical to that in the absence of the plate. The computations for the configuration in Figure 2 were done by first using the FEM code to obtain the surface electric fields over the non-metallic portion of the aperture housing the patch. These were then grouped in patches of 3x3 pixels in size and turned into equivalent magnetic current infinitesimal dipoles whose stength was set equal to (zMz + qM,) a S6 Sz, where Mz and MO are the surface 2R. Sliva and H. Wang provided the original cylinder model with the antenna cavities. 42

current densities over the pixel and a 6 6z is its area. The calculated radiation patterns for the case where a = 45~ are shown in Figure 16 and seen to be in good agreement with measurements. In particular, the pattern based on the UTD code is in better agreement with the measured data in the shadow region (below the plate) but both high frequency codes give nearly identical results in the lit region. This is primarily because APATCH employs the less accurate PO integration to obtain the fields in the shadow region below the plate. In contrast, the UTD code uses geometrical optics and diffraction theory which is known to be more accurate for flat plates with straight edges. This better performance of the UTD code is even more pronounced when the patch antenna is brought closer to the plate's surface as shown in Figure 17. Clearly, this is due to the proximity of the patch to the plate's surface (about 1A), resulting in stronger secondary currents on the plate's surface and thus the PO approximation is no longer accurate. The UTD pattern is in good agreement with the measured data except near the shadow boundary associated with the plate because the UTD diffracted fields account for the curvature of the near-zone antenna field. For practical aircraft configurations, the antenna will be placed at far distances from the wing and in that case the APATCH code is more attractive because of its geometrical adaptability and capability to handle materials in the context of the PO approximation. 6 Conclusions We presented a rather simple hybridization of finite element and high frequency codes for antenna pattern calculations in the presence of a complex structure such as an aircraft. Basically, the finite element code was employed to generate the aperture fields on the antennas surface and these were then turned into equivalent magnetic currents. These currents were subsequently used as the sources to the high frequency codes and the antenna pattern was calculated as the sum of the direct antenna radiation pattern and that generated by ray interactions with the substructure. The accuracy of this procedure was verified by comparing the calculated results with measured data. Since the measured set-up consisted of 43

canonical components, the hybridization with the UTD code leads to more accurate results in the shadow region. However, for general airframe configurations, the hybridization with the APATCH code which combines the PO, PTD and SBR methods is more attractive due to its geometrical adaptability and capability to handle material coatings. Although we did not consider antenna loading effects due to the substructure re-radiation, this can be easily incorporated into the analysis by executing the finite element code in the presence of the antenna feed excitation and the secondary excitation arriving at the antenna aperture after undergoing reflection and/or diffraction. References [1] J. M. Jin and J. L. Volakis,"A finite element- boundary integral formulation for scattering by three-dimensional cavity-backed apertures," IEEE Trans. Antennas and Propagat. Vol.39, pp.97104, January 1991. [2] J. L. Volakis, A. Chatterjee and L. C. Kempel, "A review of the finite element method for three dimensioanal scattering," J. Optical Society of America A, pp.1422-1433, April 1994. [3] J. L. Volakis, A. Chetterjee, and J. Gong, "A class of hybrid finite element methods for electromagnetics: A review," J. Electromagnetic Wave and Applications, Vol.8, No.9/10, September 1994. [4] J. Gong, J. L. Volakis, A. C. Woo, and H. T. G. Wang, "A hybrid finite element-boundary integral method for the analysis of cavity-backed antennas of arbitrary shape," IEEE Trans. Anten. Prop., Vol.42, No.9, pp.1233-1242, September 1994. [5] J. M. Jin and J. L. Volakis,"A hybrid finite element method for scattering and radiation from microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas and Propagat., Vol.39, pp.1598-1604, November 1991. 44

[6] J. L. Volakis, J. Gong, and A. Alexanian, "Electromagnetic scattering from microstrip patch antennas and spirals residing in a cavity," Electromagnetics, Vol.4, No.1, pp.63-85, 1994. [7] J. R. Mosig and F. E. Gardiol, "General integral equation formulation for microstrip antennas and scatterers," IEE Proceedings, Vol. 132, Pt. H, No. 7, 1985. [8] K.-L. Wu, J. Litva, R. Fralich, and C. Wu, "Full wave analysis of arbitrarily shaped line-fed microstrip antennas using triangular finite-element method," IEE Proceedings, Vol. 138, No. 5, October 1991. [9] R. F. Harrington, Field Computation with Moment Methods, Macmillan, New York, 1968. [10] R. S. Hansen, ed., Geometrical theory of diffraction, IEEE Press, 1981. [11] T. B. A. Senior and J. L. Volakis, Approximate Boundary Conditions in Electromagnetics, IEE Press, 1995. [12] H. Ling, R. Chou, and S. W. Lee," Shooting and Bouncing Rays: Calculating the RCS of an arbitrarily shaped cavity," IEEE Trans. Antennas Propagat., Vol.37 pp.194-205, 1988. [13] A. C. Woo, J. L. Volakis, G. Hulbert, and J. Case, "Survey of volumetric grid generators," IEEE Anten. Prop. Magaz. Vol.36, No.4, pp.72-75, August 1994. [14] A. C. Woo, T. Ozdemir, J. L. Volakis, G. Hulbert, and J. Case, "A survey of volumetric grid generators, Part II: Commercial grid generators," IEEE Anten. Prop. Magaz. Vol.37, No.1, pp.96-98, February 1995. [15] L. C. Kempel and J. L. Volakis,"Scattering by cavity-backed antennas on a circular cylinder," IEEE Trans. Antennas Propagat., Vol.42, No.9, pp.1268-1279, 1994. [16] L. C. Kempel, J. L. Volakis, and R. Sliva,"Radiation by cavitybacked antennas on a circular cylinder," IEE Proceedings-H, 1995. 45

[17] M. Gilreath, Personal communication, NASA Langley Research Center, September 1994. [18] J. M. Jin and J. L. Volakis,"A biconjugate gradient FFT solution for scattering by planar plates," Electromagnetics, Vol.12, No.1, pp.105-119, 1992. [19] R. W. Freund,"Conjugate gradient type methods for linear systems," SIAM J. Sci. Stat. Comput. Vol.13, pp.425-448, 1992. [20] R. W. Freund,"A transpose-free quasi-minimal residual algorithm for non-hermitian linear systems," SIAM J. Sci. Cornput., pp.470-482, 1993. [21] L. C. Kempel, "Radiation and scattering from cylindrically conformal printed antennas," Ph.D. Thesis, Dept. of Elect. Engr. and Comp. Sci., University of Michigan, Ann Arbor, MI 48105, April 1994. [22] J. Gong and J. L. Volakis, "An efficient and accurate model of the coax cable feeding structure for FEM simulations," to appear in IEEE Trans. Anten. Prop., 1995. [23] G. A. Thiele and T. H. Newhouse,"A hybrid technique for combining moment methods with the geometrical theory of diffraction," IEEE Trans. Antennas Propagat., Vol. AP-23, pp. 62-69, January 1975. [24] W. D. Burnside, C. L. Yu, and R. J. Marhefka,"A technique to combine the geometrical theory of diffraction and the moment method," IEEE Trans. Antennas Propagat., vol. AP-23, pp.551558, July 1975. [25] G. A. Thiele, "Overview of selected hybrid methods in radiating system analysis," IEEE Proc., Vol.80, No.1, January 1992. [26] R. J. Marhefka and W. D. Burnside,"Numerical Electromagnetic Code - Basic Scattering Code: NEC-BSC (Version 2), Part I: User's Manual," Tech. Rep. # 712242-14 (713742), ElectroScience Laboratory, Department of Electrical Engineering, The Ohio State University, Columbus, Ohio 43212. 46

[27] Users Manual for XPATCH, DEMACO. 100 Trade Centre Drive, Champaign IL, 1993. [28] D. J. Andersh, M. Hazlett, S. W. Lee, D. D. Reeves, D. P. Sullivan, and Y. Chu," XPATCH: A high-frequency electromagneticscattering prediction code and environment for complex threedimensional objects," IEEE Antenna and Propagation Magazine, Vol.36, No.l, February 1994. [29] J. Baldauf, S. W. Lee, L. Lin, S. K. Jeng, S. M. Scarborough, and C. L. Yu,"High frequency scattering from trihedral corner reflectors and other benchmark targets: SBR versus experiment," IEEE Trans. Antennas Propagat., pp.1345-1351, September 1991. 47

o / / Interference Lobe Figure 1. Illustration of antenna pattern illumination and interaction with platform. 48

1.219 m, 60.96 cm antenna aperture dia = 30.48 cm ~ 60.96 cm' cylinder / wing Figure 2. Measurement model for evaluating the hybrid FEM-High Frequency implementation. (a) _ _ _ —- -.. T Sabc'.' Patches Region I C,,,,,,,ompofi te skin...........'. "'.~""!'""" "'' * """:........... Composite skin.... Composite skin Metal Region 1 1I~ (b) Figure 3. (a) Antenna array on cylindrical platform, (b) FE-ABC modeling of an antenna element

Skewed Tetrahedron Triangular Prism Cylindrical shell Figure 4. Some commonly used tessellation elements for generating volume mesh. Surface normal \urface geometric fidelity Distorted prism x Figure 5. Illustration of the prismatic mesh and element. 50

Figure 6. Modeling of patch antennas using tetrahedral elements. 51

(a) (b) - 10 -, 10 -0 -20 - - - MEASUREMENT -20 MEASUREMENT I.0 __ __ __ _E -3 _ _I _ _ ~_ 0 - -30.. -l30. Z -180 -90 0 90 180 Z -180 -90 0 90 180 Angle (0) [deg] (O=90~) Angle (~) [deg] (8=90~) ABC mesh Y.... Yji ~ —— t te~_ termination (z=0) surface /0.64"`~~.., BI mesh / (,.1 termination surface (c) ~x x ~1.047" 12" Frequency = 3.52 GHz 0.69" 1.973" Figure 7. Comparison of measured and calculated radiation patterns for a rectangular patch on a cylindrical platform, (a) without dielectric covering, (b) with dielectric covering, (c) test body: circular ogive cylinder. 52

l ~o - P: cK X- ln Sao 0 - 3 1w Probe Fr'*' FEM IE./ \ (Voltage -. FEM /a cal Model) e | /7 -\~ Model), (Voltage, Model) o C'.- -Measured 0 ^ Probe 2a 5 2.5 2.5 Zee 2.7 275 2. 2.6 2.1 25 2S5 2. 260 2.7 275 2. 2.6 Frequency (GHz) Frequency (GHz) Figure 8. Input impedance calculations using an improved coaxial cable modeling. 53

Input Impedance vs. Frequency r = 1.3cm (0.15 WL at resonance) r = 2.17cm (0.25 WL at resonance) "00 -- -- a — -' —.-.T-' —. _ — 2o 0,',l\ r, = 2.6cm (0.3 WL at resonance) X 50 I', _ 150 a^ ^^r = 1.3cm (0. 15 WL at resonance) X. - 0.8cm 0 4 -50 ^ Q 2.0 2.5 3.0 3.5 4.0 4.5 5.0 Frequency (GHz) 3 -- Finite Element - Boundary Integral Metal (FEMA-TETRA) rMel r= 1-j2.7 -. — Finite Element - Metal Backed Absorber 1. (FEMA-PRISM) Resonant frequency = 3.55 GHz 0.4cm Error in resonant frequency estimation = 1% e = 2.4 Vertically directed field magnitude under the patch 500 V/cm 200 100 for input impedance computations. 54

GO PO-SBR obv 1 obv 1 reflected S \, SBR Integral \reflected X/ obv 2 \ PO iobv: direct incident ray raya escape O' ~incident sousource source Figure 10. Comparison of GO and PO/SBR methodologies. First Bounce Multi-Bounce toward source facet subfacet wavelength ray density Figure 11. Implementation of PO-SBR in Apatch for facet models. 55

fy... L" 92-mm 94 GHz \ g 28 mm <. WI igur 1. 5 mmrn << ----- 84mMk.ri eme.-, — Figure 12. Cassegrain geometry and ray trace. 56

Cassegrain H-plane Pattern 40 30 10 t % _ 20 - ~1 0 >*| Apatch E*| measurement -20 -15 -10 -5 0 5 10 15 Azimuth (deg) Figure 13. Comparison of Apatch with measurement for Cassegrain antenna. 57

/2.36" --- Probe feed 10.79" [9~j ~ _ iriTci 1.97" v_ - XFreq.= 3.3 GHz Patch 4^ —- 1.18" ~- Cavity 0.44", 0.031" Er= 2.17 Figure 14. Patch antenna geometry and material properties. -5 of plate w co -15 f E -20 - 0 — E~ — v N0~. 1 I at = )=28.7 h-"25/ |V/ ~ I =45 shown in Fig. 5. The shown pattern in the absence of the plate was -30.. calculated using FEMA-CYL 58

10 0IM -2 a45 -10 -3 -150 -100 -50 0 50 100 150 Observation Angle' (deg) Figure 16. Comparison of measured and calculated patterns when the pattch antenna is placed at a = 45~- using the hybridization of the finite element and high freguency codes. 1I~~~~ ~ -40 1 Meas.'" i —BSC -50. i^Y * ij | Apatch U -60 -150 -100 -50 0 50 100 150 Observation Angle 4 (deg) Figure 16. Comparison of measured and calculated patterns when the pattch antenna is placed at a = 4528.7~ using the hybridization of the finite element and high freguency codes. 10 =- -30-100 -40 Mea s. a-n a t the patch antenna is placed at ca = 28..70 using the hybridization of the finite element and high frequency codes.

Efficient Finite Element Simulation of Slot Antennas Using Prismatic Elements Jian Gong, John L. Volakis* and H.T.G. Wangt Abstract A hybrid finite element - boundary integral (FE-BI) simulation technique is discussed to treat narrow slot antennas etched on a planar platform. Specifically, the prismatic elements are used to reduce the redundant sampling rates and ease the mesh generation process. Numerical results for an antenna slot and frequency selective surfaces are presented to demonstrate the validity and capability of the technique. 1 Introduction It has been reported that a hybrid finite element-boundary integral technique [Jin, et. al. 1991, Silvester and Pelosi (1994)] can be employed for characterizing conformal antennas of arbitrary shape [Gong, et. al. (1994)]. Indeed, planar/non-planar, rectangular/non-rectangular designs, ring slot or spiral slot antennas with probe, coax cable or microstrip line feeds can be simulated with this formulation. This is because of the geometrical adaptability of tetrahedral elements used for the implementation. However, in practice, certain configurations require extremely high sampling rates due to the presence of fine geometrical details. Among them are a variety of slot antennas (spirals, rings, slot spirals, cross slots, log-periodic slots, etc.), where the slot width is much smaller than the other dimensions (cavity diameter or inter-distance of slots). In these cases, the mesh is extremely dense (with over 50, 100 or even higher samples per wavelength), whereas typical discretizations involve only 10-20 elements per wavelength. This *Jian Gong and John L. Volakis are with the radiation laboratory, University of Michigan, Ann Arbor, MI 48109-2122. tH.T.G. Wang is with the Naval Air Warfare Center, China Lake, CA 93555. 60

dense sampling rate is especially severe for 3-D tetrahedral meshes, where the geometrical details usually distort the tetrahedrals. The numerical system assembled from this type of mesh often leads to large system conditions due to the degraded mesh quality. Also, mesh generation is tedious and the solution CPU time is unacceptably large. In this paper, we propose a finite element-boundary integral formulation using edge-based triangular prism elements. It can be shown that this element choice is ideally suited for planar antenna configurations, including spirals, circular and triangular slots. Among the many advantages of the prismatic elements, the most important is the simplicity of mesh generation. Also, much smaller number of unknowns is required for an accurate and efficient modeling of complex geometries. Below, we begin by first outlining the finite element-boundary integral (FE-BI) formulation for slot antenna modeling. A new, physically meaningful, set of edge-based functions for prisms is then presented to generate the discrete system of equations. The final section of the paper gives results for antenna radiation and transmission through frequency selective surfaces. Comparisons with reference and measured data are given and the efficiency of the implementation is discussed. 2 Formulation Consider the cavity-backed slot antenna shown in Fig. 1 where the cavity is recessed in a ground plane. To solve for the E-field inside and on the aperture of the cavity, a standard approach is to extremize the functional F(E) = 2l J(V x E(V E (V x E E)- k -OE}dV + JJJ E (jkoZoJi + V x T M) dV + jikoZoff E(H x )dS (1) jko+Sf where e and p denote the relative tensor constitutive parameters of the cavity medium, Zo and ko are the free space impedance and propagation constant, respectively, So represents the aperture excluding the metallic portions and Sf denotes the junction opening to the guided feeding structures. Also, V, is the volume occupied by the source(s) and H is the corresponding magnetic field on So and Sf whose outer normal is given by n. The explicit knowledge of H in (1) is required over the surface So and Sf (referred to as mesh truncation surfaces) for a unique solution of E. Specifically, the 61

Planar ~ ^.. Ground $,: Plane.:..',Sf * -i::~:-s _ s~.'.A.". Figure 1: Geometry of cavity-backed microstrip antennas magnetic field H over So may be replaced in terms of E via a boundary integral (BI) or absorbing boundary condition (ABC), whereas H on Sf is determined on the basis of the given feeding structure. In this paper, we will employ the boundary integral method [Jin, et. al., 1991] for truncating the mesh, a technique commonly referred to as the finite element - boundary integral (FE-BI) method. In the context of the FE-BI, H is represented by the integral H = H90 + Jj [z x E(r')]. G(r, r') ds', (2) where G is the electric dyadic Green's function of the second kind [Tai (1994)] such that n' x (V x G) = 0 is satisfied on the (planar, spherical or cylindrical) metallic platform. For the antenna problem shown in Fig. 1 where the platform is a planar ground plane, G becomes the half space 1 e-j4co Ir-r'I G 2jkW(I+ VV) -r, (3) with r and r' being the observation and integration points, respectively, and I -= + yy + zz is the unit dyad. In connection with our problem, i.e. 62

that of a cavity recessed in a ground plane, H90 is equal to the sum of the incident and reflected fields for scattering computations, or zero for antenna parameter evaluations. To discretize the functional (1), we choose to subdivide the volume region using prismatic elements as shown in Fig. 2 and Fig. 3, The field in each of the prisms can be approximated using the linear edge-based expansion [Nedelec (1980), Webb (1993), Bossavit (1989)] 9 E = EeVj = [V]T{Ee}, (4) j=1 where [V]e = [{Vx}, {Vy}, {Vz}], and {Ee} = {E, Ee,..., E}T. The vectors {V,}, u = x, y, z, are of dimension m = 9 and they simply represent the x, y, z components of V associated with the jth edge of the eth element. Since Vj are chosen to be edge-based functions, the unknown coefficients Ej represent the average field along the jth edge of the eth element. A corresponding representation for the aperture fields is 3 E(r) = E S(r) = [S]T{E'}, (5) i=1 where [S], = [Sx,Sy], and V(r) reduces to S8(r) when the position vector is on the slot. To generate the discrete system for Eje, (4) and (5) are substituted into (1) and subsequently F(E) is differentiated with respect to each unknown EJ. With the understanding that the surface field coefficients EJ are a subset of Ee, we obtain OFV N N Nv N. {Ee = Z[A]{Ee} + E[B] {E9} + E{Ke} + E{L} = 0 (6) TE " e=l e=l e=l s where the sums are over the total number of volume or surface elements. In this, the matrix elements are given by [A]e = JfJ {[DV]e[ T'][DV] - kg[VIe[r] T[V]} dv (7) [ f Jix Mix {K}e M= JJJ [Vi e jkoZo {J + v x [ L1] Miy dv (8). I 63 63

Fieure 2: Illustration of tessellation using prisms A z 2 1 Z=Zc 6Z i=1,2,3 j=4,5,6 k=7,8,9 Figure 3: Right angled prism 64 9 00 8 ea s n 5 tot 4 NONE go Z=Z c 6 i=1,2,3 j=4,5,6 k=7,8,9 64n P

[OV]T I —}T v [DV]ze - {ZV}T (9) Li- ay [BS] = Jj j{-2ko]s[] +2 [{S } - [ {,5}-T _ at~.{$}T]} Go(r r')dS dS' () go\ Js8 \ -Hxg/ {L,} - jkoZos [$] _Hgo )dS (11) where {L,} is removed in case of radiation problems and the same holds for {Ke} when the scattering problem is considered. 3 Edge-Based Elements Consider the right angled prism shown in Fig. 3 whose vertical (z-directed) sides are parallel (right-angled prism). The height of the prism and the triangle area will be designated as h=.'(r-) = rei z x (r- ri) (12) where ri denotes the location of the ith node, ei is the unit vector along the ith triangular edge, li denotes the length of this edge and r is any position vector terminated inside the triangle. One way to obtain an edgebased field representation for the prism is to utilize the nodal basis functions [Zienkiewicz (1989)] and then apply the procedure discussed in [Nedelec (1980), Bossavit and Mayergoyz (1989)]. However, an alternative and more physically meaningful approach can be employed for the construction of the edge elements. Referring to Fig. 2, it is evident that if r is in the x-y plane, then Se in (12) gives the area of another triangle 12'3' such that the lengths of edges joining the nodes 2- 3 and 2' - 3' are equal. With this definition of r, the vector ii si = 2S; x (r-ri) (13) 65

has a magnitude which is equal to the ratio of the areas of the triangle 12'3' to that of 123. We observe that (13) is simply the edge-based expansion for the triangular elements [Rao, et. al. (1982)] and is the appropriate expansion to be used in (5). The corresponding volumetric basis functions can be obtained by inspection, viz. (z - Az) Vi =( Si i-=1,2,3 Az Vj =(Z+ Z) j = 4,5,6 (14) Az Vk = (k k=7,8,9 where (k is the triangle simplex coordinate associated with the kth prism vertex at (xk, Yk). As illustrated in Fig. 3, zc and h = Az represent the offset coordinate and the prism height, respectively. When (14) are substituted into (7), the resulting integrals can be evaluated in closed form as given in the Appendix. However, the integrals resulting from the substitution of (13) into (10) must be carried out numerically, except the self-cells which must be performed analytically as discussed by Wilton (1981). 4 Applications Radiation and scattering by an Annular Slot: To evaluate the accuracy and efficiency of the prismatic mesh and the aforementioned implementation, we first consider the analysis of the narrow annular slot (0.75cm wide) shown in Fig. 4. The slot is backed by a metallic circular cavity 24.7 cm in diameter and 3 cm deep. The FE-BI method is quite attractive for this geometry because the slot is very narrow and most of the computational requirements are shifted on the finite element portion of the system. The calculation shown in Fig. 5 were carried out using the prismatic and tetrahedral elements [Gong, et. al. (1994)]. As seen, they overlay each other. However, only 1024 prisms were needed for modeling the cavity, whereas the number of the tetrahedral elements for this homogeneously filled cavity were 2898 for acceptable element distortion. If a multi-layered structure was considered, or a similar system condition was used as a criterion for mesh generation, then much more tetrahedrals than prisms would be needed for modeling such a structure. Moreover, the prismatic mesh is trivially generated given the slot outline. In contrast, substantial time investment is required for generating and post-processing the tetrahedral mesh. Frequency Selective Surfaces (FSS): FSS structures [Pelton and Munk (1979), Mittra et.al. (1988)] are arrays of tightly packed periodic 66

,.' ^-', a = 12.35 cm b - 0.75 cm,t be, p0o= 7.7 cm ~S-.! / 0.7 < f< I GHz d 3 1m [ e - 1ee"35 Figure 4: Geometry of the annular slot antenna backed by a cavity 23.7 cm in diameter and 3 cm deep -.1~0,~~ ~ ^BIBBB'^ (~~~~~~-10-50 -.' _, JO!'''3 -a0.....-. -70 -70 0 10 20 30 40 5 0 70 6 0 8 0 100 100 -0 -60 -40 -2 0 0 20 40 60 80 100 tate (digmr) ~w tt(dqgrM) Figure 5: Scattering: Bistatic (co-pol) RCS patterns computed using the tetrahedral FE-BI code and the prismatic FE-BI code. The normally incident plane wave is polarized along the X = 0 plane and the observation cut is perpendicular to that plane. Radiation: X-pol and Co-pol radiation patterns in the 4 = 0 plane from the annular slot antenna shown in figure 4. The solid lines are computed using the tetrahedral FE-BI code whereas the dotted lines are computed using the prismatic FE-BI code. The excitation probe is placed at the point (y=0) marked in figure 4. 67

elements which are typically sandwiched between dielectric layers. The periodic elements may be of printed form or slot configurations designed to resonate at specific frequencies. As such, they are penetrable around the element resonances and become completely reflecting at other frequencies. To meet bandwidth design specifications, stacked element arrays may be used in conjunction with dielectric layer loading. Here we shall consider the analysis of FSS structures (with slot elements) via the FE-BI method. Because of the fine geometrical detail associated with the FSS surface, the finite element method has yet to be applied for the characterization of FSS structures, but use of prismatic elements makes this a much easier task. Of particular interest in FSS design is the determination of the transmission coefficient as a function of frequency, and since the array is periodic, it suffices to consider a single cell of the FSS. For computing the transmission coefficient T, the periodic cell is placed in a cavity as shown in Fig. 6 and the structure is excited by a plane wave impinging at normal incidence. Assuming that near resonance the wave transmitted through the FSS screen will retain its TEM character, the transmission line concept can be used to find the scattered field oaT2 Es = 1- aR where T is the transmission coefficient of the FSS, R = 1 - T and a is the reflection coefficient associated with the cavity base. To reduce the multiple interactions within the cavity, it is appropriate to terminate the cavity with some absorber, thus reducing the value of a to less than 0.1. Then, since R is also small near resonance, a good approximation for T is T() = 10 log -E dB - a and upon considering the next higher order cavity interactions, we have TdB T() +10 log [1 - a(1 -T(~)). A more direct and traditional computation of TdB would involve the placement of the FSS element in a thick slot [Jin and Volakis (1991)]. However, this requires enforcement of the boundary integral over the entire lower surface of the slot, leading to a much more computationally intensive implementation. The above FSS modeling approach was applied for a characterization of single layer and multi-layer FSS structures. In both cases, the periodic 68

1391 cm. Radome Cover (circular on e i4.5 substate) l 0.61~ cmA /0.0762cm: - 2 0.0762cm 1.2 cm Metal Absorber Figure 6: Illustration of the setup for computing the FSS transmission coefficient Upper figure: periodic element (top view); Lower figure: periodic element in cavity (cross-sectional view) 69

PMM Calc...on........ Measure-' /!.. * FE-BI -5 0 2 4 6 8................-... 6 1 - 25 * _................ Figure 7: Calculations and comparisons of transmission through the FSS structure shown in Fig. 6 element was a slot configuration. The geometry of the single layer periodic element is shown in Fig. 6 and consists of a planar slot array on a dielectric layer 0.0762 cm thick and having cE = 4.5. The FE-BI calculation using prismatic elements is given in Fig. 7. Clearly, our calculations are in good agreement with the measurements and data based on the more traditional PMM approach [Berrie (1995), Henderson (1983)]. The geometry of the multilayer radome considered in our study is given in Fig. 8. The total thickness of the FSS was 6.3072cm and is comprised of two slot arrays (of the same geometry) sandwiched within the dielectric layers. For modeling purpose, a 1.54cm thick absorber is placed below the FSS as shown in Fig 8. From the calculated results, it is seen that the results generated by the FE-BI method are in good agreement with the measurements. 5 Conclusion A hybrid finite element-boundary integral (FE-BI) formulation was presented for modeling narrow slots in metal backed cavities. Prismatic elements were used in connection with the FE-BI implementation, and in contrast to the tetrahedral elements, these offer several advantages. Among 70

slot array placed 60mils below top surface of 90mils layer slot array at 30mils below top of 90mils layer Er 4 O9 mils r 7m 1 < 5.m.85 ~.~ 2cm 90mils c. absorber,E r 4 - lOJ 104 6.45cm -10 C I — -20 / | measured calculated -25 -30 0- 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 frequency (GHz) Figure 8: Upper figure: geometry of the multilayer frequency selective surface (FSS) used for modeling; lower figure: measured and calculated transmission coefficient through the FSS structure 71

them,low sampling rates are needed for generating meshes and the mesh generation process is substantially simplified. Other advantages of the prismatic elements over the tetrahedral elements include better system conditions and faster pre/post data processing. The explicit expressions for FE-BI implementation of prismatic elements were tabulated and numerical results for slot antennas and frequency selective surfaces were presented to demonstrate the validity and capability of the technique. Appendix For FEM implementation, the following quantities are required Pn = Vx Vm V x VndV (15) Qmn = Vm n dV (16) where the curls are given by V x Vi = -2Se [(x - xi) + (Y - yi)( - z(z - i =12, 3 i4 VxVj = 2 [(x- j) +(y-yj)y+z(zc+ az-z)] j=4,5,6 (17) 2SeLz V Vk = [(k2 - kl) + (k2-ykl)Y] k = 7,8,9. To this end, we follow the notation defined in (13) and (14), where i, i'=1,2,3 represent the top triangle edges, j, j'=4,5,6 denote the bottom triangle edges and k, k'=7,8,9 stand for the vertical three edges. It is found that (16) and (17) can be analytically evaluated and we tabulate the results as follows Pii' = Cii' Diiz + Se(Az)3] (18) Pjj, = Cjj, [Djj,Az + 43e(AZ)3 (19) Pkk' = S ke (20) 4 S e ( 2 0 ) Pij = Pj, = -Cij [DijA - S(AZ)3] (21) Pik = Pki = -4(Se)2 [ * l(SX- iSe) + y.l (SY -yiSe)] (22) 72

?jk - Pkj-4(se)2 = 4(Se)2 [ -(SX- zijSe) + yl.kSY-yS)] (23) Qii = (-3 CiiDii (24) (AZ)3 Qjj, = ( —-Cjj3 j (25) Qkk' = AZSeTkk' (26) (AZ)3 Qij = Qji= CjDi (27) Qik = Qki = Qjk = Qkj = 0 (28) where Tkk' = 1/6 for k = k'; 1/12 for k $ k';ilj ci = 4(SeZ)2 (29) Dij = SXX - (xi + xj)SX + xixjSe + SYY - (yi + yj)SY + yiyjSe The remaining quantities in the above list of the expressions are defined as Se = / dxdy SX = x dxdy SY = s ydxdy SXX = x2 dxdy SYY = L y2dxdy SXY = xy dxdy These integrals can be expressed in terms of the global coordinates of the three nodes (Xi, Yi), (Xj, Yj), (Xm, Ym). Specifically, assuming that the three nodes i,j and m of a triangle are in counterclockwise rotation, we then have, - 1 Xi yi S =s dxdy 1 x y 1 Xm ym SX = xdxdy= S-(Xi + Xj + Xm) 73

SY = Jydxdy ( + Yj + Ym) SX = J xd2 ddy = (Xi +Xj+Xm(+ X + + )} SYY = y2 dzdy = {(+( j+Ym)2 (yi+ y y2)} SXY = zy ddy- = {(xi + xj + Xm) (Yi + Y + Ym) + (XiYi + XjYj + XmY,)} 6 References Berrie,J. (1995). Personal communication, Mission Research Corp. Bossavit, A. and Mayergoyz, I. (1989).'Edge-elements for scattering problems', IEEE Trans. magnetics, vol. 25, no. 4, pp. 2816-2821. Gong, J., Volakis, J.L., Woo, A., and Wang, H. (1994).'A hybrid finite element boundary integral method for analysis of cavity-backed antennas of arbitrary shape', IEEE Trans. Antennas Propagat., vol. 42, pp. 1233-1242. Henderson, L.W. (1983),'The scattering of planar array of arbitrary shaped slot and/or wire elements in stratified dielectric medium', Ph.D. dissertation ElectroScience Lab, Ohio State University Jin, J.M., Volakis, J.L. and Collins, J.D. (1991).'A finite element-boundary integral method for scattering and radiation by two- and three-dimensional structures', IEEE Antennas Propagat. Magazine, vol. 33, no. 3, pp. 22-32. Jin, J.M., Volakis,J.L. (1991).'Electromagnetic scattering by and transmission through a three dimensional slot in a thick conducting plane', IEEE Trans. Antennas Propagat., vol. 39, pp. 534-550. Mittra, R., Chen, C.H. and Cwik, T. (1988).'Techniques for analyzing frequency selective surfaces -A review', Proc. IEEE, vol. 76, pp. 1593-1615. Nedelec, J.C. (1980).'Mixed elements in R3', Numer. Math., vol. 35, pp. 315-341. 74

Pelton, E.L. and Munk, B.A. (1979).'Scattering from periodic arrays of crossed dipoles', IEEE Trans. Antennas Propagat., vol. AP-27, pp. 323330. Rao, S.M., Wilton, D.R. and Glisson, A.W. (1982),'Electromagnetic Scattering by surfaces of arbitrary shape', IEEE Trans. Antennas Propagat., Vol. AP-30, pp. 409-418. Silvester, P.P. and Pelosi, G. (1994),'Finite Elements for Wave Electromagnetics - Methods and Techniques', IEEE Press Tai, C.T. (1994),'Dyadic Green Functions in Electromagnetic Theory', (2nd Ed.) IEEE Press Webb, J.P. (1993).'Edge elements and what they can do for you', IEEE Trans. Magnetics, pp. 1460-1465. Wilton, D.R. and Butler, C.M. (1981),.'Effective method for solving integral and integro-differential equations', Electromagnetics, vol. 1, pp. 289-308. Zienkiewicz, O.E. and Taylor, R. (1989). The Finite Element Method, 4th ed., New York: McGraw-Hill. 75

Design of Planar Absorbing Layers for Domain Truncation in FEM Applications S.R. Legault, T.B.A. Senior and J.L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2122 Abstract A metal-backed layer of absorbing material offers a number of advantages for truncating the computational domain in a finite element simulation. In this paper we present design curves for the optimal selection of the parameters of the layer to achieve a specified reflection coefficient. The curves are based on one-dimensional finite element simulations of the absorbers, and the optimization is therefore a function of the sampling rate. Three types of material are considered, including the recently introduced perfectly matched uniaxial material, either homogeneous or with a quadratic material profile. A three-dimensional application to a microwave circuit is also presented and used to examine the validity of the design curves. 76

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

2 Analytical Study (onsider the metal-backed planar layer slown il Fig. l(a). The surface.r = 0 i: llhe ilc ('rfaclc between free space (in x < 0) and a lossy material (in.7 > 0) backeid by a P'C( at a = /. For an incident plane wave E: or H = -Iko(rcoso+ysino) (1) the reflection coefficient is RE(W) or RH(O). and the objective is to minimize these. If the layer is composed of a homogeneous isotropic material whose relative permittivitv or and relative permeability p, are such that er = = b = a - j3 (say). then 1 b- 2 sin2 - j cos 0 tan(kobt\ - b-2 sin2 ) RE(0) - / * (2) 1 - b- 2 2sin2 j + jcos tan(kobt/ - b2 sin2 () 1 - b-2 sin2 4- j cos ) cot(kobt 1 -b-2 sin2 4) RH() =- / (3) /1 - -2 sin2 + j cos 4 cot(kobt 1- b- 2 sin2 ) These differ because the presence of the PEC backing has destroyed duality, and at grazing incidence () = ir/2), RE = RH = -1. If sin 4 > 1, i.e. 4 = ir/2 + j' with 6 > 0 so that sin X = cosh 6 and cos 4 = -j sin 6, IRI differs from unity by only a small amount for both polarizations. The behavior of IRE,H()Il as a function sin q is illustrated in Fig. 2. At normal incidence (4 = 0), (2) and (3) give RE,H(O) = Fe-2jk~(- i) (4) whose magnitudes are independent of a and can be made as small as desired by choosing kot3 sufficiently large. As kot - oo, RE,H() -+ 0 only for normal incidence, but Sacks et al [4] have shown that a particular uniaxial anisotropic material has this property for all 4 < 7r/2. The result is an example of a perfectly matched layer (PML), and if =r, = bU - (b- b)X& (5) where I is the identity tensor, then RE,H() = Fe-2jkot(a-j()cos (6) which reduces to (4) in the particular case of normal incidence. If koh > 1 the reflection coefficients decay exponentially for all 4 < 7r/2, and since (6) is also valid for sin > 1, the choice a > 0 ensures an exponential decay for these angles as well. The behavior of IRE,H(4)) is illustrated in Figure 2 for the same values of kot, a and #f used for the isotropic layer. Clearly, a major advantage of the PML is that its reflection coefficient remains low for a wide range of angles of incidence. 78

Although the outer surface of the layer is reflect ionless for all o. the i l)rplll)t i (lalil i tlie material properties' at x = 0 may produce a contribution in an INFFI I solltion. \\ can eliminate the discontinuity by tapering the properties as a function of. to prod ucel anll inhomogeneous anisotropic laver. As shown by Legault and Senior [7]. if 1) =) wave propagating into the material has the form e-jko {rr(r) cos o+y sin} (7) and when y(x)= 1 + (a- j/-1) () (8) which tends to unity as x -- 0+, the reflection coefficients of the laver are identical to those given in (6). With this expression for y(x), the attenuation is less where the field is larger, i.e. close to the interface, and increases as the field is absorbed. A simplified version of (8) is employed in Section 3.4. 3 Numerical Study For all three types of layer the theoretical behavior of jRI is relatively simple. In the case of the isotropic material, an increase in / and/or t decreases IR(0)I. The uniaxial material has this behavior for all real angles of incidence, and while a plays little or no role, large values of a do produce higher absorption for complex angles. It follows that for a uniaxial layer of given thickness t, a and /3 can be chosen sufficiently large to produce high absorption over a broad angular spectrum, with angles near = 7r/2 providing the only exception. Unfortunately, the analytical results do not immediately translate into numerical performance. Because of the discretization inherent in an FEM implementation, the fields inside the layer are reproduced only approximately, and this is particularly true for a rapidly decaying field. To design a good absorber it is necessary to understand the impact of the sampling rate on the choice of a, /3 and t, and our objective is to find the minimum number of sampling elements ( or discrete layers) to achieve a specified JR\. It is anticipated that the errors introduced by the discretization will have a number of consequences. In particular, for a given number N of discrete layers and given t, increasing / will ultimately lead to an increase in JIR because of the inability to model the increasing attenuation, and an increase in a will likely produce a similar effect. To obtain some insight into the roles played by N, a, / and t, we now consider a simple FEM model of the layers. 3.1 Numerical Model A simple one-dimensional FEM code was used to examine the numerical performance of the absorbing layers. The computational domain was limited to the discretized layer structure shown in Fig. l(b), with the appropriate boundary conditions applied at the interface x = 0 and the PEC backing x = t. 79

\\e consider first a homogeneous isotropic layer at normal ilcildence tor wlill tllic tilt theoretical reflection coefficients are given in (4). In spite of the fact tliat thle Illiljit l( nd' ale the same for b)oth polarizations. a polarization dependence shows up il t lie M1 Illplementation. This is illustrated in Fig. 3. and we note that as.\ increases. tl(he IE.L valutes of iR(0)I converge to the common theoretical value for both polarizations. 3.2 Dependence on a and f3 For a layer of constant thickness the theoretical value of IR(0)I is independent of a and polarization, but in the numerical implementation the behavior is much more complicated. Figure 4 shows IR(0)\ plotted versus a and i3 for a layer of thickness t = 0.25Ao made up of 5 (=N) elements, where the darker tones indicate lower values. For small 3 the results are in close agreement with theory. As evident from the level lines, IR(O)j is almost independent of a and decreases exponentially with /3, leading to a linear decrease on a dB scale. For large /, however, the behavior is quite different, and the most striking feature is the series of deep minima whose spacing in a increases with increasing a and decreasing /. These are numerical artifacts which are common to both polarizations and may depend on the particular numerical code employed. The minima for the two polarization are interlaced, and for H polarization the first minimum occurs at a = 0, /3 = 1.6. Their locations also depend on t and N. If N is fixed, the spatial sampling is inversely proportional to t. Decreasing t results in better sampling, pushing the minima to higher values of / and producing agreement with the theoretical values for larger /3 than before. Increasing t has the opposite effect. On the other hand, if t is fixed, increasing N improves the accuracy, and shifts the minima to higher /3. Apart from the minima, the reflection coefficients for fixed / increase slightly with increasing a, and it is therefore sufficient to confine attention to the lower values of a. In Figure 5 the reflection coefficients are plotted as functions of P for the same layer with a = 0 and a = 0.75. The curves correspond to vertical cuts through the patterns in figure 4, and we also show the theoretical value obtained from (4). We observe that as /3 increases the reflection coefficients decrease initially at almost the same rate implied by (4), but beyond a certain point they begin to increase. The deep minimum at a = 0 and f = 1.6 in Figure 4(a) is clearly seen, but for design purposes it is logical to focus on the worst case, i.e. the polarization for which the reflection coefficient is larger. The upper curves in Figure 5 are almost identical and constitute this case. Since they correspond to two different values of a, either of them would suffice, but for reasons that will become evident later, we choose a=0. 3.3 Dependence on /3, N and t We now seek a connection between the values of /, N and t for which IR(0)1 is minimized. To this end, we first examine IR(0)I as a function of /3 and N for constant t, and the resulting plot is shown in Figure 6 for E polarization with t = 0.25Ao and a = 0 as before. For fixed 80

.3 tlhe reflection coefficient tends to its exact values as.\ increases.'1 i i> evitl'II(t fIomIl illi level curves and. as expected. the convergence is better for tlie smaller.'. (onisilder iio\\w tli b)ehavior of IR(O)I for fixed N'. As 3 increases from zero. the reflection (oefficient (ld(real(s to a minimum and then increases. The location of the miniina is indicated )\by t le solidl line. This is consistent with the behavior shown in Figure.5 and the upper curve is. in fact. just a vertical cut through Figure 6 at A' = 5. The solid line in Figure 6 therefore gives tlie value of 3 at which IR(O)I is a minimum as a function of the number of elements. If the process is repeated for other layer thicknesses, it is found that for minimum IR()0) the curve of ft/Ao versus N is virtually the same for all thin layers. The observation that O3t/Ao is a scalable parameter is an important conclusion of our study, and by choosing a constant layer thickness we can produce a universal curve for the optimal choice of AN and 3 in FEM simulations. Such a curve is shown in Figure 7 and can be interpreted as giving the value of fit/Ao for a specified N to minimize the reflection coefficient IR(O)I. For example. if t = 0.2Ao and N = 3, then 0 = 2.13. If a smaller value of 3 is chosen, IR(0)} will be larger (see Figure 5), and this can be attributed to the fact that the field reflected from the metal backing has not been attenuated sufficiently. If / is set to a value larger than 2.13, IR(O)I will still be larger because the chosen N is too small to reproduce the rapid field decay within the layer. Although the scaling property of ft/Ao has only been established for a = 0, it holds to a reasonable degree for small ca 0, but as c increases, the Pft/Ao versus N curves become increasingly dependent on a. The scalability also extends to the associated values of jR(O)I, and this enables us to provide a simple design prescription for an absorbing layer. 3.4 Design Curves Since the quantities lit/Ao and IR(0)| are the same for layer thicknesses up to about 0.5Ao at least, design curves can be obtained by plotting IR(0)| and N versus ft/Xo on the same figure as shown in Figure 8. To see how to use the figure, suppose that the desired reflection coefficient at normal incidence is -50 dB. In Figure 8 we observe that the IR(0)\ curve intersects the -50 dB line at it/Ao r 0.58, and referring now to the N curve, the number of elements required is N = 10. The value of / can then be found by specifying either the element size or the layer thickness. Thus, for elements 0.025A0 thick, we have t = 0.25Ao and 3 = 2.32. By increasing N we can improve the performance up to the limit provided by the theoretical value of IR(0)I which has been included in Figure 8. A good approximation to the long-dashed curves in Figure 8 obtained by linear regression is t = -0.01061Rj + 0.0433 (9) AO N = 0.147exp [7.353Pt/Ao] (10) where IRI is measured in dB and N is the smallest integer equal to or exceeding the right hand side of (10). 81

So far %we have considered only a homogeIneou isotropic laver at lorllal ilncidencllc. Inlt <, evident from (4). (6) an'd Figure 2. thle important feature of the anisot roic uiaial Ic mtterial is that it provides almost the same reflection coefficient for a range of aniles of inl cidenlc( about normal. Thus. for the homogeneous uniaxial material. tle (lesigli clurves ill li11ur1e are applicable for these angles of incidence as well. The performance can be iIllprove(l 1b making the anisotropic material inhomogeneous. and to illustrate this we consider tile case,(x) = -jf (t) for which the theoretical reflection coefficient is the same as before. Tlhe scalability is still preserved and the resulting curve is shown in Figure 8. The fact that the curve for the quadratically tapered layer lies below that of the homogeneous material confirms the improvement in performance, and we can now achieve a reflection coefficient of -50 dB by choosing ft/Ao " 0.64 corresponding to N = 9. Approximations to the short-dashed curves in Figure 8 are t = -0.01191R + 0.0451 (11) Ao N = 0.298exp [5.263/t/Ao] (12) where IRj and N are as before. Compared with the homogeneous material the decrease in the number of elements required becomes more pronounced as \RI is reduced. 4 Three-Dimensional Verification As noted earlier, a PML is particularly attractive for terminating a finite element mesh in the simulation of microwave circuits. For these applications a PML has an advantage over a traditional ABC because it does not require a priori knowledge of the guided wave propagation constant. It is therefore of interest to examine the performance of the layer when used to terminate a microwave transmission line (see Figure 9) and to see how this relates to the design curves generated on the basis of the one-dimensional simulations. The microstrip line has width w = 0.71428 cm, substrate thickness 0.12 cm and relative permittivity r, = 3.2, and is enclosed in a metallic cavity whose dimensions are shown in Figure 9. It should be noted that the height of the cavity from the microstrip line is sufficiently large to suppress any reflections from the cavity walls. As a result, the characteristic impedance of the line should be that same as if the line was in free space. The microstrip line was terminated using a two-section homogeneous uniaxial absorber having material parameters r,,,r in the upper section and 3.2=r, r. in the lower section to match the substrate. The calculations were carried out at several frequencies using an FEM code [8] and we show the results for 4.0 GHz. At this frequency the element width was 0.05Ao and a five layer absorber having a total thickness of t = 0.25A was used. With a = 0 the* computed reflection coefficient of the transmission line structure as a function of fi is shown in Figure 10. For comparison, with t = 0.25Ao and N = 5 the design formulas (9) and (10) give fi = 1.92 and IRI = -41 dB. However, these are based on matching free space, and 82

in the present instance the optinium design.i mlust (be scaied. itecoM4i i/ilz tl I11 11 Iu>1 f Iil( power is confined to the substrate. the.3 required for colmparison witl tl' (le coll)lln't (ldat i: 1.92//3.2 = 1.07. and this is in reasonable agreement with tile location of tl liIiilllllll ill Figure 10. The fact that the minimum \RI is lower than predicted is. perha)ps. not suIrprlisisin. We recall that the design curves are based on the worst case. i.e. tile polarization plrovidingi the largest minimum IRI. and the curve in Figure 10 resembles more thle H polarization curve in Figure 5 than the E polarization which constitutes the worst case. 5 Conclusion A uniaxial perfectly matched layer provides a powerful means for truncating finite element meshes close to the modeled structure. By properly selecting the material properties and sampling rate, almost any desired level of absorption can be attained, and typically very few samples (less than five) are needed to achieve a reflection coefficient of -40 dB over a wide range of incidence angles. In this paper we described a detailed study of three types of layer material including homogeneous and inhomogeneous uniaxial ones, and by identifying the scalable parameters of the layers, universal design curves and formulas were developed. The curves or formulas can be used to specify the numerical, geometrical and electrical parameters of the PML to achieve a desired absorption down to -60 dB or lower. As expected, a lower material loss requires a thicker absorber to produce the same reflection coefficient. On the other hand, a higher attenuation rate requires more samples to attain a lower reflection coefficient in a numerical implementation. As expected, an inhomogeneous (tapered) PML is better than a homogeneous one since the material loss can be increased to larger values close to the metal backing where the field is smallest. To test the applicability of the design criteria in a three-dimensional setting, a shielded microstrip line was considered. With the line terminated in a homogeneous PML, the results were in reasonable agreement with prediction, and the discrepancies were no more than could be expected in view of the conditions under which the criteria were established. These conditions are: (i) use of a particular one-dimensional FEM code (ii) based on the worst case polarization, i.e. the polarization for which the minimum reflection coefficient is largest (iii) restriction to a range of angles of incidence about normal (iv) assumption of a pure imaginary propagation constant, i.e. a = 0, in the layer. Condition (iv) is a requirement for scalability, and though small values of ax are still admissible, the condition is clearly inappropriate if there is substantial power at complex angles of incidence for which large a is required for absorption. If the polarization can be specified, (ii) is also inappropriate, and the design criteria may underestimate the performance that 83

can bie achieved. In an! given prob)lem where there is t he' luxury! of tlesti in' \ rlit t ot I8 \!(v specifications. it is probable that a performance cal ibe achieved wi\ch iv l)et tIer tllil tlt predicted by the criterion. but even then the design values are a logical p)lace to stalrt. In t he more likely situation where prior testing is not feasible. awe believe tiat tlie design criterlia provide a looical basis for specifying the parameters of the PNIL and its samplin,. Acknowledgment We are indebted to Ms Lin Zhang for her assistance with the calculations for the microstrip line. References [1] T.B.A. Senior and J.L. Volakis, Approximate Boundary Conditions in Electromagnetics. Stevenage, UK, IEE Press, 1995. [2] C. Rappaport and L. Bahrmasel,"An absorbing boundary condition based on anechoic absorbers for EM scattering computation," J. Electromagn. Waves Appl. 6, pp. 16211634, 1992. [3] T. Ozdemir and J.L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Sci. 29, pp. 12551263, 1994. [4] Z.S. Sacks, D.M. Kingsland, R. Lee and J.F. Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition," to appear in IEEE Trans. Antennas Propagat. [5] J.P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comp. Phys. 114, pp. 185-200, 1994. [6] D.S. Katz, E.T. Thiele and A. Taflove, "Validation and extension to three dimensions of the Berenger PML absorbing boundary conditions for FD-TD meshes," IEEE Microwaves and Guided Wave Lett. 4, pp. 268-270, 1994. [7] S.R. Legault and T.B.A. Senior, "Matched planar surfaces and layers," University of Michigan Radiation Laboratory Report, 1995. [8] J. Gong, J.L. Volakis, A.C. Woo and H.T.G. Wang, "A hybrid finite element-boundary integral method for the analysis of cavity-backed antennas of arbitrary shape," IEEE Trans. Antennas Propagat. 42, pp. 1233-1242, 1994. 84

Figure 1 (Gometry of the metal-backed absorber layer \ a aii its l [' l illl ll'lilt'ill tll 101. Figure 2 Analytical results for homogeneous isotropic and anisotropic absorblilIo layers with t = 0.25Ao and b = 1 - j2: (. ) isotropic E pol.. (- -) isotropic H pol. and (- ) anisotropic. Figure 3 Numerical results for a homogeneous isotropic layer with t = 0.15Ao and b = -j2.5. The top four curves are for E pol.. the bottom four for H pol.: ( ) exact, (- - -) N=3. (- - -) N=6 and (- -) N=12. Figure 4 Plot of JR(O)j in dB for (a) an H polarized and (b) an E polarized wave incident on a homogeneous isotropic layer with t = 0.25Ao and N =5. The solid curves are level lines. Figure 5 IR(0)I with t = 0.25Ao and N = 5: (-) Exact, (- -) a = 0 E pol., (- - -) a 0 H pol., (- - -) a = 0.75 E pol. and ( ) a = 0.75, H pol. Figure 6 IRI as a function of,3 and N for E pol. with t = 0.25Ao and a = 0. Figure 7 /t/Ao computed from the E pol. case with a = 0: (- ) t = 0.15Ao, (- -) t = 0.25Ao and (- — ) t = 0.5Ao. Figure 8 Absorber design curves. The straight lines give IRI in dB, and the curved ones give ( —) exact IRI, (- -) homogeneous case and (- - -) inhomogeneous case. Figure 9 Geometry for the microstrip line. Figure 10 Computed IRI for the microstrip line at 4 GHz with a = 0, t = 0.25Ao and N = 5. 85

/ -JOS2 98 (q) (R) N.. x A *

-10, — -20, -"'"" -60.; -40 -.4. 0.0 0.2 0.4 0.6 0.8 1.0 1 -60 - o. -80 I 0.0 0.2 0.4 0.6 0.8 1.0 1.2 sin ) 87 Ff6oec 2

0 -10 -20 m -30-50 "^ -60 - -70 -" —--- 0.0 0.2 0.4 0.6 0.8 1.0 1.2 88 S.oDr 7

(q) Volt 9, ~ L 0 -0 9 *C L 0 WR"'. W'VO* )BID; |~'=~=F/ _r~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

0 -10 -20 - _o -30.....III -40- \.-. \,-50 - II!re -60 ---— \II 0 1 2 3 4 5 90

5 2 4 6 8 10 12 14 16 18 20 N 91

0.7 0.6 0.5 0.4 0.3 0.2 1- - -- - I' I- i. I 2 4 6 8 10 12 14 16 18 20 N 92 F/sUCE- rl

IRI IdBI 00o 0 en v3 X 00 0 0 0 0 0 0 o0 o o o o o o e I' I 1 / I.o 0c ^ \ /*. I / /\ / * * / /' — O / / v "-"- 00 01 C- M CO N

Uniaxial /"Artificial Absorber Cavity, ii m _r[ 3.0953cm 0.21 cm 94

UNIVERSITY OF MICHIGAN II11111111 111: 1111 IIi l Illlll 3 9015 03527 1652 -15 -20 -25 -30 -?L-35 - -40 - -45 - -50 - -55 0 0.5 1 1.5 2 2.5 3 95