028918-2-T Electromagnetic scattering and radiation from microstrip patch antennas and spirals residing in a cavity J. L. Volakis, J. Gong, and A. Alexanian Northrop Corp. B-2 Division 8900 E. Washington vBlvd. Pico Rivera CA 90660

TABLE OF CONTENTS Page 1. Introduction 2 2. Formulation 3 2.1 Interior Region Formulation 4 2.2 Exterior Formulation 8 3. Computational and Modeling Considerations 1 1 3.1 Mesh Generation 1 1 3.2 Matrix Assembly and System Solution 1 3 4. Numerical Results 15 Example 1: RCS of a single residing in a rectangular cavity Example 2: RCS of a single patch near resonance Example 3: RCS of rectangular patch in a rectangular cavity Example 4: RCS of a loaded rectangular patch Example 5: RCS of a RAM coated patch Example 6: RCS of a patch with distributed resistive loading Example 7: RCS of a cavity-backed spiral Example 8: Input impedance of a rectangular patch Example 9: Input impedance of a circular patch 5. Conclusion 20 References 20 Figures 23-38 Appendix A: Matrix Elements for Rectangular Bricks 3 9 Appendix B Evaluation of Matrix Elements for Tetrahedrals 42 Appendix C: Biconjugate Gradient Algorithm... 45 Appendix D: Boundary Integral Implementation with Tetrahedral Elements 46

Electromagnetic scattering and radiation from microstrip patch antennas and spirals residing in a cavity J. L. Volakis, J. Gong, and A. Alexanian Abstract A new hybrid method is presented for the analysis of the scattering and radiation by conformal antennas and arrays comprised of circular or rectangular elements. In addition, calculations for cavitybacked spiral antennas are given. The method employs a finite element formulation within the cavity and the boundary integral (exact boundary condition) for terminating the mesh. By virtue of the finite element discretization, the method has no restrictions on the geometry and composition of the cavity or its termination. Furthermore, because of the convolutional nature of the boundary integral and the inherent sparseness of the finite element matrix, the storage requirement is kept very low at O(n). These unique features of the method have already been exploited in other scattering applications and have permitted the analysis of large-size structures with remarkable efficiency. In this report, we describe the method's formulation and implementation for circular and rectangular patch antennas in different superstrate and substrate configurations which may also include the presence of lumped loads and resistive sheets/cards. Also, various modelling approaches are investigated and implemented for characterizing a variety of feed structures to permit the computation of the input impedance and radiation pattern. Many computational examples for rectangular and circular patch configurations are presented which demonstrate the method's versatility, modeling capability and accuracy. 1

1 Introduction In designing low observable vehicles, the radar cross section (RCS) of conformal microstrip patch or cavity antennas and arrays is of crucial importance. An example of such an antenna configuration is illustrated in Figure 1 and in computing its RCS, it is necessary to consider the contributions of the radiating elements as well as those of the feed and cavity substructures. Conventional analyses based on the space domain and spectral domain integral equations [1]-[7] make use of the infinite substrate Green's function and are thus best suited for computing the RCS near resonance. In that case, the antenna size is relatively small and structural terminations are not as important. Outside resonance and particularly at higher frequencies, the moment method solutions lead to large systems which are not practical for including contributions due to the feed and substructure of the conformal antenna or array. In this report, we review a hybrid finite element method for computing the scattering of conformal antenna structures while taking into account the presence of the feed geometry and the cavity substructure. This method was originally presented by Jin and Volakis [8]-[11] and has recently been generalized to allow modeling of non-rectangular patches [12] and more intricate feed structures and distributed loads such as resistive sheets [13], often used for RCS reduction purposes. An important feature of the method is the use of the finite element method to discretize the cavity region occupied by the substrate, feed structures and radiating elements. This results in a sparse matrix for the solution of the cavity fields and is the primary advantage of the method because it leads to O(N) matrix system, where N denotes the number of edges used in constructing the volume mesh. At the aperture of the cavity, the exact boundary integral equation is enforced to terminate the mesh (i.e. not an absorbing boundary condition) and thus, in principle, the proposed hybrid finite element method is exact. To retain an O(N) storage requirement for the overall system, it is necessary to have a structured grid at the aperture allowing use of the FFT in conjunction with an iterative solver [14], [15], [16]. In the following, we first present the mathematical formulation of the method, allowing for the presence of lumped loads, resistive sheets and feeding cables or microstip lines. The resulting functional is discretized by introducing tetrahedral or rectangular elements. In particular, edge-based ele 2

ments are employed which have a number of inherent advantages over traditional node-based elements. Specifically, the edge-based elements are divergenceless within each element, thus, eliminating spurious solutions; they maintain tangential field continuity at dielectric interfaces without a need to introduce additional constraints; and are suitable for modeling corners, where certain field components are singular. Using codes based on the presented finite element boundary integral formulation, we compute the RCS of different configurations which demonstrate the accuracy of the method as well as its adaptability in modeling a large class of configurations. 2 Formulation Consider the cavity configuration illustrated in Figure 1. The cavity is recessed in a ground plane as is usually the case with conformal antennas on aircraft platforms. The actual antenna or array elements are placed either on the aperture or within the cavity, in which case they radiate in a substrate/superstrate configuration. It is assumed that the elements are fed in various ways, including via a stripline or a microstrip line, coaxial cable or alternatively, they may be aperture coupled to the feeding lines. Furthermore, this formulation allows for the cavity material and geometry to be arbitrary, and the same is true for the shape of the patches. The patches may be further loaded with resistive cards or lumped elements for bandwidth or frequency control. For the described cavity-backed antenna configuration, we are interested in computing its radar cross section over a wide range of frequencies. In accordance with the method introduced by Jin and Volakis [11], the problem described above can be formulated by separating the computational domain into two regions, to be referred to as regions I and II. In region II, which encompasses the volume occupied by the cavity, the finite element method will be employed for discretizing the fields. The primary reason for using the finite element method is its adaptability in modeling a variety of cavities and radiating elements whose geometries and material composition may vary widely. The fields in region I (exterior region) will be computed via the boundary integral method. This amounts to introducing equivalent sources on the cavity aperture which are then integrated to obtain the radi 3

ated fields. The final step of the formulation is to couple the fields in each region by enforcing tangential field continuity at the aperture. When the final equations are discretized, the resulting system involves a sparse section and a full submatrix. The first is associated with the interior fields and, thus, involves the bulk of the unknowns. The full submatrix is associated with the field components on the boundary and by resorting to a structured grid, this submatrix is Toeplitz which is important for retaining the O(N) storage requirement of the proposed hybrid finite element method. 2.1 Interior Region Formulation The goal with any partial differential equation (PDE) solution method, such as the finite element method, is to solve the vector wave equation V x (V x E) - korE -jkoZoJ' + V x (Mi) (1) where (Cr, r) denote the relative constitutive parameters of the material filling the cavity, ko is the free space wavenumber and Zo is the free space intrinsic impedance. In addition, Ji and Mi denote the impressed (if any) electric and equivalent magnetic sources, respectively. These are generally set to zero for scattering calculations but we shall retain their presence for the sake of generality. Instead of solving (1) directly, it is more appropriate for numerical purposes to construct its equivalent weak-form. This can be readily done via the weighted residual method. Accordingly, (1) is first weighted with the function T(r) giving the expression (R,T) =JI {VX (-V E) T-k6rE.T} d + T. jkoZoJ + V x (-M) dv = 0 (2);JJv [J (1\M)r ] in which V denotes the volume occupied by the cavity. It is desirable for numerical calculations to make the first volume integral of the weighted residual (R, T) symmetric with respect to E and T. To accomplish this we recall the identity JfJ V x (1 Vx E. Tdv = JJJI (V x E). V x T dv 4

+ jkoZo s T (H x )ds (3) in which S encloses the volume V, h denotes the outward unit normal to S and H = V x E/(jkoZo) is the magnetic field. When (3) is employed into (2), we obtain (R, T) = 1 (f(V x E) (V x T) - krE T) dv + jkoZo T [(Hi + H + H) x n] ds + JJT. jkoZoJ + V x (IMi) dv = 0 (4) i//v [J (ur )] where the closed surface in (3) has been replaced by the cavity aperture surface Sap since E x n vanishes on the interior cavity walls. Consequently, we have set H = Hi + Hr + Hs, where Hi denotes the incident plane wave excitation and H' is the reflected field from the ground plane with the aperture removed. The remaining field Hs is that scattered by the aperture and must be related to the electric field on Sap either through a local or a global boundary condition. A local boundary condition, such as an absorbing boundary condition, will only provide an approximate relationship between E and H on Sap and is not therefore considered in the context of this formulation. In contrast, a global boundary condition such as the Stratton-Chu integral equation is an exact relation between E and H on Sap. The treatment of the surface integral over Sap will be discussed later in more detail. For the moment we return to the overall numerical solution of (4) for the field E. To generate from (4) an appropriate linear system, it is necessary to first discretize the volume V in smaller volume elements. The field in each element can then be approximated by introducing some linear or higher order expansion. Subsequently, by using different testing or weighting functions T = Wj (j = 1,2,..., N) a linear set of equations can be generated for the solution of the field expansion coefficients. Some typical volume elements for discretizing V are illustrated in Figure 2 and any one of these can be employed for meshing the cavity volume. Obviously, elements such as rectangular bricks are best suited for rectangular patches, whereas tetrahedrals can be used for discretizing cavities which house antenna elements of more arbitrary shape. An advantage of the rectangular bricks is their geometrical 5

simplicity, leading to a rather simple geometry preprocessor. In contrast, the mesh generator associated with the tetrahedral elements must be far more sophisticated. Nevertheless, regardless of the element type used, the fields within each element (say the eth element) can be expressed as Nv E= E EeW (5) i=l In this, N, denotes the number of edges forming the element (twelve for bricks and four for tetrahedrals) and We are the element basis or expansion functions. For our analysis they have been chosen such that the expansion coefficients Ee represent the value of the field component along the ith edge of the eth element when evaluated at the location of the same edge. Substituting the expansion (5) into (4) and setting T = Wj (Galerkin's testing) gives Flv + Fs = FI + FF (6) where Ne Ne Nv Flv = EF ev = E eeAe- (7) e=l e=l i=1 Ne NT FF = F= Ej Ie (8) e=l e=l Nse Nse Fs = F, F, Z F Qi E Q? (9) e=l e=l eESap eESap with AJ = ( -(V W) ( V X W J ( X ) - k W: w) dv (10) Ke = - J/ W { jkozoJi-V x ( M)}dv (11) Fs = jkoZo S S ' [H x h] ds (12) 2jkZjse s [H x ] ds (13) Q = -2jkoZo/ S [Hi x n] ds (13) 6

The number of volume elements comprising V has been denoted by N, whereas the number of surface elements comprising Sap has been denoted by N8e. For this application, nh = and Sj is the same as Wj when the last is evaluated on the face of the volume element coinciding with Sap. If resistive cards, lumped loads or coaxial cables are considered, the lefthand side of (6) must be supplemented by other terms. Specifically, in the presence of a resistive card (of resistivity R, occupying the surface Sr = E Se), we must introduce the function Nsr Nr eESr i=1 r eESr i= = ELe enx S\ ~ (n x S i= ) do = E E e (14) in which Se denotes the surface that belongs to the eth element, has Nsr sides and coincides with Sr. Correspondingly, Se are the We basis functions when evaluated on Se. When a lumped load or conducting post is placed at the ith edge of the eth element, we must then introduce the function F3v = j ko~Et ~ *W WJ = E C!ej (15) ZL with Wj evaluated at the location of the load or conducting post. Usually, Cej = 0 for i:7 j. In (15) ZL denotes the impedance of the load and 1 is its length, equal to the length of the ith edge. Finally, when the antenna is fed with a coaxial cable, we must add to the left-hand side of (6) the term (see fig. 3) F = J Nk0 E EZ Z eI S> - e-~ p ds eEScab i=l cab eEscab i=1 As usual, Scb denotes the surface over the face of the eth element whose face, having NSc sides, coincides with Scab. The basis functions Sj are the same as Wj when the last is evaluated on Sb and Io is the strength of the current excitation on the center conductor. Also, p is the radial distance measured from the same center conductor and p is the associated unit vector along the radial direction. We note that the expression (16) was derived on the assumption that the coaxial cable supports the lowest order mode. 7

Using the aforementioned matrix elements, (6) can be rewritten in discrete form as Ne {[A] + [Bt] + [Cr] + [D]} {Ee} + Fs e=l Ne Nse = E IK)+ E {Q) E {Pe} (17) e=l e=l eEScab eESap After the element assembly, the above sum results in a sparse matrix which is the primary advantage of the method. The sparsity of this matrix is assured because the domain of the shape functions Wj is restricted within a single element. Obviously, a single edge can be shared by more than one element leading to more equations than unknowns which are reduced to the actual number of unknowns during the assembly process. The evaluation of the matrix elements Aj, Bj, etc. is straightforward once the basis functions We are specified. The edge-based expansion functions for rectangular bricks and tetrahedral elements are given in the Appendices along with the corresponding expressions for the matrix elements. 2.2 Exterior formulation As discussed earlier, to solve the system resulting from (4) or (6), it is necessary to replace the magnetic field appearing in the surface integral in terms of the electric field, i.e. the working variable. An exact relation between the tangential electric and magnetic fields at the aperture is.2ko 1 H- H+ HI- J oo (I + k-VV) Go(r,r'). [E(r') x z] ds' (18) where Go(r, r') is the free space Green's function -jko olr-r'l Go(r, r') = r' (19) with r and r' representing the observation and integration points, respectively. Also, I denotes the unit dyad. Substituting (18) into (4) or (12) and invoking the divergence theorem twice along with the identity V (GoE) = E. VGo + GoV E (20) 8

it follows that Fs = 2k2 lko ( x S(r)) * {s [E(r') x A] Go(r,r') ds ds (21) a'p ^ ap +2//s Vs* (S (r) x ) {/ V* [E(r') x ] Go(r, r') ds' ds ap p On the aperture, the expansion (5) becomes Ns E(r) = E;' so'(r) (22) i=l where Se (r) are the basis We (r) given in the Appendices when evaluated on Sap which is the surface element bordering the e'th volume element. When (22) is introduced in (21), we obtain Nse Ns Fs= E 'G e (23) e'=l i=l with G = 2klfs [ x Ss(r)] { JJ [Se'(r') x z] Go(r,r')ds'} ds (24) + 2Js Vs * (S;(r) x V) {ffs V* [St'(r') x A] Go(r, r') ds'} ds We observe that the matrix elements G~' will be non-zero even when Sj(r) does not overlap with S '(r). In contrast, the matrix element A, given by (10) will be zero whenever We does not overlap with the expansion function We. Consequently, the matrix [Gf.e] will be full, whereas the matrix resulting after the assembly of [Arj] will be very sparse and, thus, mostly populated by zeros. To proceed with the numerical evaluation of the elements Gfe, it is necessary to introduce the explicit expressions for the basis functions S (r). For rectangular bricks, these are equal to the basis functions Np = - Wj" (p = x or y) defined in appendix A with the local z variable set to zero. With these basis functions, the surface area Sap is a rectangle and therefore standard 9

numerical integration routines can be used except when e' = e. In that case, Go(r, r') is singular and the integrals must then be evaluated with care. In accordance with the above, when the eth and the e'th elements are not the same or adjacent, then midpoint integration can be used to find lGe % -2apSaSpGo (r, r) {k2S(X, y). S (x', y ) (25) + 2 SV [SS(xe c) x z] V * [S(xy) x z] Go(r, r')} where Sp and Sap are the areas of the associated elements, rc = xcx + Ycy is the midpoint of the eth surface element and likewise r' is the midpoint of the e'th surface element. When the eth and e'th elements are adjacent or coincident, the rectangular surface areas can be subdivided into 3 x 3 small rectangles. On replacing the expansions functions with their midpoint value at each of the nine sub-elements we get 3 3 eete 2 2Sm Se G/ j 2 Z Z gmnpq,{-kS(rmn) Si (rpq) m,n=l p,q=l + VS [Sj(rmn) X Z]VS * [S (rpq) x Z]} (26) where rmn and rpq denote the midpoints of each sub-element, and gmnp =JJ {IJ Go(r,r')ds' ds (27) rMn;; The kernel of gmnpq is singular when r = r', in which case it is appropriate to rewrite it as gmnpq = JjS [JjI (Go(r, r) -4- ) ds' ds +/ II --- ds'ds /S, JJs4r- r ' If uJlSe' 47Ir - r'I ds' ds The first of these integrals is now well behaved and can be evaluated numerically using a 2 x 2 Gaussian integration. The second integral, although singular, can be evaluated analytically as described in [17]. In the case of tetrahedral elements, the surface basis Sj(r) become the same as those given by Rao, etc. [18] when defined for a pair of triangles which share an edge. The evaluation of the matrix elements G_' for this case is more involved and is best accomplished using area coordinates as done in [19, 20]. 10

3 Computational and Modeling Considerations 3.1 Mesh generation A first step in any numerical analysis is the generation of an appropriate discrete model of the structure. For our case, this amounts to subdividing the cavity into bricks or tetrahedrals, a procedure referred to as mesh generation. For rectangular cavities and patches, the mesh generation using rectangular bricks is fairly simple. An approach which was employed is to subdivide the cavity volume into layers of thickness Az. Each layer is then readily subdivided into Nse x Nse rectangular bricks. If the cavity housing the antenna structure is nonrectangular, rectangular bricks can still be used, but in that case staircasing will be unavoidable. The discretization of the cavity using tetrahedrals requires a more sophisticated mesh generation scheme and several commercial packages are available for this purpose. We have used a package marketed'by SDRC IDEAS which can be run on many different workstation platforms. Regardless of which mesh generation scheme is employed, the final result is a file often called the Universal file. As a minimum, this file contains the following Tables: 1. A list of the global volume element Nos. and their material code. 2. A list of the nodes forming each volume element. 3. A list of the (x, y, z) coordinates for each node. 4. A list of global edge Nos. and a list of the pairs of global node Nos. forming an edge. The last table is not usually contained in the Universal file but can be readily constructed from the first three tables. To the above lists we must also add information pertaining to the antenna geometry, feeding lines, resistive cards, lumped loads and about the presence of a coax cable or a probe feed. This information is generated by a preprocessor based on information provided by the user. For example, in specifying the location of a rectangular patch or stripline, the user may enter the z location of the patch along 11

with the (x, y) coordinates for two opposite corners. For non-planar metallic surfaces, it is necessary to flag or group the appropriate surface elements bordering the metallic surface. This must, of course, be done when generating the mesh. A similar procedure is required for identifying the presence of resistive cards. The resistivity must be specified for each resistive surface element or patch. The specification of lumped loads or conducting posts is relatively simple. These are assumed to be placed at the edges of the volume elements and consequently the mesh must be generated with this in mind. Once the mesh is completed, the user specifies the location or edge number and the impedance value of each load or post. If the post is perfectly conducting, i.e. a short, it is flagged as such and during matrix assembly, the field coefficient associated with the edge occupied by the short is set to zero. The center conductor of the coaxial cable is defined in this manner. Based on the antenna and cavity structural and loading specifications, the following Tables are generated: 5. List of edge Nos. on a perfectly conducting surface (cavity wall, patch, stripline, wires, metallic posts, etc.). 6. List of volume elements along with their global edge Nos. which border the aperture and do not coincide with a patch or other metallic surface. 7. List of volume elements and global edge Nos. bordering a resistive surface along with the associated resistivity. 8. List of global edge Nos. along with subglobal edge Nos. for combinations. The preprocessor (see flow chart in fig. 4) which can be run interactively, with or without a graphical user interface, depending on which mesh generation package is employed. The matrix elements are computed directly from the data provided by these tables using the algorithms outlined in the previous section. The excitation vectors (to the right of (17)) are readily filled once the user specifies the location and type of source for antenna parameter calculations. For scattering calculations {KI } and {Pe} are set to zero with {Ql} being the only excitation vector. 12

3.2 Matrix Assembly and System Solution Once all input tables are read, the goal is to assemble the matrix system EA.{n } + 1[ [~] 1 {E } {bv (2 [ ] } + [[0] [] {E} (28) where [A] is a very sparse square symmetric matrix, [Q] is an Nse x Nse square full matrix, the column {EV} denotes the field values at the edges not bordering Sap, and {EJ } are the corresponding coefficients for the edges on Sap. For scattering computations, the column {bV} is set to zero and {bA} is equal to {Qme as defined in (13). All matrix coefficients in the above system are numbered in accordance with the unique global or their own subglobal numbering schemes. If the pth global edge lies on a metallic surface or coincides with a metallic wire or post, Ev or EB is zero and is, thus, not included in the unknown field column. Consequently the computation of Ap, is skipped. The non-zero matrix elements Amn are computed by summing A-j, Be-,' C, and Ds. Note that with this notation, m and n refer to global edge numbers, whereas i and j are local indices for the eth element. Obviously, an edge may belong to more than one volume element and in that case A is constructed by summing the element equations for that edge. It is important to note that only the non-zero elements of the [A] matrix are stored, thus, ensuring the O(N) memory demand of the system. The [Q] matrix is the same as the [Gfe] matrix, once the local indices of the last are transformed to global incides. Because this matrix is full, it is important to employ uniform gridding to permit the computation of the associated matrix-vector products via the FFT. That is, for rectangular bricks, all elements of Sap must be chosen to be rectangular and equal. This may place some minor restrictions on the separation of the patches but is well worth it, because it leads to substantial computational savings. In the case of tetrahedrals, we must choose all surface elements to be right triangles of equal size. Again, this is not difficult to achieve with the mesh generator, but may place some restrictions on the specific shape and separation of the elements. The [9] matrix is symmetric by its nature, and by employing uniform discretization as noted above it is necessary to store only one of the rows of each Toeplitz block comprising the matrix. Many solvers can be used for solving this system, including LU decomposition or iterative schemes. Because 13

of the sparsity of [A] and the special form of [a], iterative solvers are best suited as they require O(N) storage when the fast Fourier transform is used for computing the product [g]{Ef}. Specifically, we employed a biconjugate gradient algorithm for symmetric systems. This algorithm is given in Appendix C and its storage requirement is 6.25Nv + 10.5Ne complex numbers. We have observed rather impressive convergence rates using this algorithm (with diagonal preconditioning) for systems having even more than 100,000 unknowns. For example, for 180,000 unknowns, convergence is achieved in 57 iterations. The CPU time per iteration on a HP9000/750 for this system was 40 sec. Another system of 25,000 unknowns required 66 iterations per angle with an average CPU time of 2 sec per iteration. Generally, the sampling rate was 12 to 15 elements per linear wavelength. Other iterative algorithms such as the Conjugate Gradient (CG) or the conjugate gradient squared (CGS) could be used instead of the biconjugate gradient method. The CG and CGS guarantee convergence and are monotonic. They achieve this by making the original matrix positive definite, a procedure which squares the condition of the original system. Since the convergence rate of the algorithm is strongly dependent on the condition of the matrix (high condition numbers are associated with slower convergence), the convergence rate of the CG or CGS algorithm is slowed down substantially. In contrast the biconjugate gradient (BiCG) algorithm does not square the condition number of the system. Consequently, its convergence rate is faster but has the drawback of not always converging. Nevertheless, for the symmetric system associated with this application, we have not experienced difficulties in achieving convergence. In fact, as noted above, the employed BiCG algorithm converged with rather impressive rates. This is in spite of the erratic residual errors. Specifically, the residual error does not decrease monotonically and this is another drawback of the BiCG. If desirable, one could instead use the QMR algorithm which is very similar to the BiCG other than its monotonic convergence curve. That is, QMR does not again guarantee convergence and is associated with faster convergence rates than CG and CGS. It has been pointed out that for layer systems (over 200,000 unknowns), the QMR converges faster than the BiCG. It should be noted though that the QMR algorithm is more involved and requires more memory than the BiCG algorithm. Consequently, since we have found the convergence rate of the BiCG algorithm to be quite satisfactory, we did not explore that QMR algorithm for this application. More details on the different iteration 14

algorithms can be found in [21, 22, 23]. 4 Numerical Results In this section we present radar cross section computations for different cavity-backed antenna configurations. Once the fields in the cavity volume are computed by solving (28), only the tangential aperture fields need be retained for evaluating the antenna RCS (or radiated field). Specifically, the scattered field is given by H"~(r)= - -j koYo r - A. H'(r) = -jkoyo-2r /fS (+0~+).[E(r') xz]ekos in(x'cos+Y'sn) ds' (29) 27rr JJSap where (r, 0, q$) denote the usual spherical coordinates of the observation point. The RCS of the cavity-backed antenna is then given by o'pq - lim 4rrz2 [H (r)'.q2 pq = l47r (30) r —*oo IH (r)12 in which p and q are equal to either 0 or q. Below we present a number of radar cross section patterns for different cavity-backed antenna configurations. The first few computational examples are intended for validation purposes, and the remaining expose the reader to RCS reduction techniques. Input impedance computation examples are also included. Example 1: RCS of a single patch residing in a rectangular cavity The geometry of this configuration is shown in Figure 5 along with the computed backscatter cee RCS pattern. The patch is 1.448" x 1.038" and resides on a dielectric substrate having thickness t = 0.057" and relative permittivity,r E 4.0. The substrate is housed in a 2.89" x 2.10" x 0.057" rectangular cavity recessed in a ground plane. The excitation was a 0-polarized plane wave at 9.2 GHz in the xz plane, and the analysis was done using rectangular bricks. As seen the computed oes pattern is in good agreement with the measured data [11]. We note that for this example the cavity terminations play a more important role because the computation frequency is not near the resonant frequency of the patch. The measured data were collected by 15

placing the cavity-backed antenna configuration on a 5 foot long low cross section body [18]. Example 2: RCS of a circular patch near resonance The geometry of this patch configuration is shown in Figure 6 along with the calculated ae0 backscatter RCS as a function of frequency. The diameter of the circular patch is 1.3 cm and resides on a substrate 0.0787 cm thick, having er - 2.33. The excitation is a plane wave incoming at (i0 = 60~, qi = 180~), and the computations were done using tetrahedral elements with the substrate terminated at a diameter of 1.6148 cm. However, because the calculations are near the patch resonance (~ 8 GHz), the substrate terminations do not contribute significantly to the RCS of the antenna. As seen, our computations are in good agreement with the measured data [19] which were collected with the patch placed in a cavity having a different periphery. Example 3; RCS of a rectangular patch in a rectangular cavity This configuration consists of a rectangular patch 0.5 cm x 0.5 cm residing on a substrate housed by a rectangular cavity 1 cm x 1 cm x 0.5 cm. The substrate's relative permittivity is 2.33. In Figure 7, we show the ass backscatter RCS pattern of this configuration at 3 GHz as a function of angle. The incident plane wave is in the xz plane (principal plane) and the computations were done using rectangular bricks and tetrahedral elements for discretizing the cavity region. As seen, both discretization schemes give identical results. Example 4: RCS of a loaded rectangular patch A common method for reducing the RCS of a patch antenna is to apply resistive loading. Figure 8 shows the relative RCS reduction which can be achieved by placing four loads at the edges of the microstrip patch. It is seen that as the load's value is reduced to zero, the RCS as well as the gain of the antenna are also reduced. As is well known [5, 20], the dB RCS reduction at resonance is twice that of the gain, and this is observed in Figure 8 when the feed's presence is modeled with a matched load (the feed is at x = 2.5 cm, y = 1.7 cm). However, the 2 to 1 dB reduction is not precisely observed because the resonant frequency of the patch is a function of the load's value. 16

The calculation frequency of 1.97 GHz corresponds to the resonant frequency with the loads set to 300 Q each. Unfortunately, as shown in Figure 9, the predicted RCS reduction is not maintained outside the operating band of the antenna. Away from the resonant frequency, the resistive loads play a much lesser role since the RCS is dominated by the antenna structure, including the cavity terminations. In a later example, we describe a treatment which is better suited for RCS reduction outside the first resonance. Example 5: RCS of a RAM coated patch Another approach for reducing the RCS of the patch is to coat it with radar absorbing material (RAM) as shown in Figure 10. In parity with the study performed in [5], this lossy coating is only over the extent of the patch, thus minimizing radiation losses. It should also be noted that this specific configuration can only be analyzed using a volume formulation and illustrates the adaptability of the finite element method. The relative reduction achieved with the placement of a specific coating as a function of its thickness t is illustrated in Figure 10. These calculations were performed at the resonant frequency of the patch which varied from 1.944 GHz at t = 0 cm to 1.893 GHz at t = 0.3 cm. As a result, the typical relation between gain and RCS reduction is observed with much more fidelity. Unfortunately, as seen from Figure 11, the observed RCS reduction is again maintained only near resonance. At higher frequencies the thickness of the coating causes the overall RCS of the structure to be increased. This occurs primarily because of its larger electrical thickness. Example 6: RCS of a patch with distributed resistive loading With the motivation of reducing the patch RCS at frequencies outside its operational band without substantial compromise in gain, we consider the placement of a narrow resistive ribbon (R-skirt) around the periphery of the patch as illustrated in Figure 12. Because the ribbon's electrical width increases with increasing frequency, it is expected to have a rather noticeable effect on the RCS of the patch at higher frequencies. Thus, greater RCS reduction can be attained at higher frequencies without compromising the antenna's gain performance. If a certain level of RCS reduction is desired at resonance, lumped loads can still be used in conjunction with the R-skirt 17

loading. The ass RCS of the patch as a function of frequency depicted in Figure 12 is given in Figure 13 before and after loading with the R-skirt. The given backscatter RCS data correspond to a plane wave excitation incident in the xz plane and at an angle of 0 = 70~ with respect to the z axis. However, as shown in Figure 13, the RCS reduction is fairly uniform over all 0 angles. As seen, the resistive skirt loading reduces the RCS in and out of the band as opposed to the narrowband performance of the lumped loads. Not surprisingly, the placement of the R-skirt is seen to slightly change the resonant frequency of the patch from 1.9522 GHz down to 1.744 GHz. From the inset in Figure 12, the RCS reduction at resonance is about 9 dB for all angles of incidences and thus the corresponding gain reduction is approximately 4.5 dB. Example 7: RCS of a cavity-backed spiral Cavity-backed spirals are popular broadband antennas whose RCS is certainly of concern because of the complex nature of the antenna. Figure 14 shows the top and cross-sectional view of such an antenna. The spiral is rectangular in shape for modeling convenience and is housed in a metallic cavity 2.8125 cm x 2.8125 cm with and without absorber at the base of the cavity. From the given dimensions, the lower limit of the antenna's operational frequency is approximately 3 GHz. We shall therefore consider the RCS of the antenna at 4.5 GHz and our goal is to study the effect of the absorber at the cavity base and that of a resistive sheet placed uniformly over the cavity aperture. The chosen absorber consists of three uniform layers and was designed to be similar to a commercial absorber giving 18 dB reduction at 4.5 GHz. The principal plane ao0 backscatter RCS of the spiral antenna using different treatments is shown in Figure 15 as a function of angle. The corresponding radiation patterns are given in Figure 16. As seen the absorber alone reduces the RCS of the antenna by approximately 6 dB and this was to be expected (the gain reduction was 3 dB). Clearly, the R card has a more dramatic effect on the RCS of the spiral antenna. Some RCS reduction is still observed even when the resistivity of the sheet is 10Zo or 5Zo. However, the gain of the antenna (not shown) drops substantially more and in particular nearly a 10 dB gain reduction was observed with R = 5Zo. For lower val 18

ues of the resistivity the gain reduction is so substantial (as is also the case with the RCS of the spiral) making such a treatment impractical. Perhaps a better treatment (which reduces the RCS without substantial compromise in gain) would be to apply a lossy dielectric coating over the spiral/cavity aperture. This coating would be equivalent to a resistive sheet having a resistivity which is a function of frequency. By controlling the thickness and loss tangent of the coating, a frequency response may be achieved which is appropriate for a specific application. Example 8: Input impedance of a rectangular patch The specific antenna geometry is shown in Figure 17 consisting of a patch 3.4 cm x 5.0 cm residing at the aperture of a 7.5 cm x 5.1 cm cavity. The substrate is of thickness t = 0.08779 and has a relative permittivity Cr = 2.17 with a loss tangent of 0.0015. Figure 17 shows the complex input impedance as a function of frequency with the probe feed placed at (1.22 cm, 0.85 cm). These calculations were done in the presence of a 50 Q resistor (lumped load) placed between the patch and substrate at (-2.2 cm, -1.5 cm). Our calculations are seen to be in good agreement with the measured data [11] whereas corresponding calculations using traditional methods are not as accurate. We remark that the 50 Q load (50 Q resistor) serves to reduce the RCS of the patch at resonance. However, it also reduces the antenna gain. It has been observed though that at resonance the dB reduction in RCS is twice that of the gain (i.e. if the gain is reduced by 5 dB, the RCS will be reduced by approximately 10 dB). Unfortunately, as noted in a previous example the resistors do not have any appreciable effect on the RCS when the patch is illuminated with a frequency out of its resonant frequency. Example 9: Input impedance of a circular patch The specific antenna configuration is a circular patch having a radius of 1.3 cm and residing at the center of a substrate filled cavity of radius 2.115 cm. The cavity depth is t = 0.406 cm (-A/15) and the substrate's relative dielectric constant is er m 2.9. We computed the input impedance of this patch placed on the surface of the cavity for a probe feed located 0.8 cm from the center of the patch. The computed and measured impedance are depicted in a Smith chart shown in Figure 18 by varying the frequency from 3 GHz 19

to 3.8 GHz. The agreement between our calculations and reference data is indeed quite satisfactory. 5 Conclusion In this report we reviewed a hybrid finite element-boundary integral method for computing the RCS of planar antennas. The pertinent finite element and boundary integral formulations were presented along with the method's attributes. It was demonstrated that by virtue of the finite element method, the formulation is adaptable and rather efficient in modeling differing geometries and material compositions, including lumped and distributed loads. General-purpose computed codes were written on the basis of the presented formulation. These were used to generate the presented RCS computations for circular and rectangular patch configurations, spiral antennas and other non-rectangular configurations. Among the presented examples some served for validating the method whereas others exposed the reader to RCS reduction schemes. References [1] M.C. Bailey and M.D. Deshpande, "Integral Equation formulation of microstrip antennas," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 651-656, July 1982. [2] D.M. Pozar, "Input impedance and mutual coupling of rectangular microstrip antennas," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 1191-1196, Nov. 1982. [3] J.R. Mosig and F.E. Gardiol, "General integral equation formulation for microstrip antennas and scatterers," IEE Proc., vol. 132, part H, Dec. 1985. [4] D.M. Pozar, "Radiation and scattering from a microstrip patch on a uniaxial substrate," IEEE Trans. Antennas Propagat., vol. AP-35, pp. 613-621, June 1987. 20

[5] D.R. Jackson, "The RCS of a rectangular microstrip patch in a substrate/superstrate geometry," IEEE Trans. Antennas Propagat., vol. 38, pp. 2-8, Jan. 1990. [6] J.T. Aberle, D.M. Pozar and C.R. Birtcher, "Evaluation of input impedance and radar cross section of probe-fed microstrip patch elements using an accurate feed model," IEEE Trans. Antennas Propagat., vol. 39, pp. 1691-1696, Dec. 1991. [7] K.A. Michalski and D. Zheng, "Analysis of microstrip resonators of arbitrary shape," IEEE Trans. Microwave Theory Techn., vol. 40, pp. 112-119, Jan. 1992. [8] J.M. Jin and J.L. Volakis, "A finite element-boundary integral formulation for scattering by three-dimensional cavity-backed apertures," IEEE Trans. Antennas Propagat., vol. AP-39, pp. 97-104, Jan. 1991. [9] J.M. Jin and J.L. Volakis, "Electromagnetic scattering by and transmission through a three-dimensional slot in a thick ground plane," IEEE Trans. Antennas Propagat., vol. AP-39, pp. 543-550, April 1991. [10] J.M. Jin, J.L. Volakis and J.D. Collins, "A finite element-boundary element method for scattering by two and three dimensional structures," IEEE Antennas & Propagat. Magazine, vol. 33, no. 3, pp. 22-32, June 1991. [11] J.M. Jin and J.L. Volakis, "Scattering and radiation from microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas Propagat., vol. AP-39, pp. 1598-1604, Nov. 1991. [12] J. Gong, J.L. Volakis, A. Chatterjee and J.M. Jin, "Characterization of cavity-backed conformal antennas and arrays using a hybrid finite element method with tetrahedral elements," 1992 IEEE AP-S International Symposium, Chicago IL, Digest pp. 1629-1632. [13] J.M. Jin, J.L. Volakis, C.L. Yu and A. Woo, "Modeling of resistive sheets in the context of finite element solutions," IEEE Trans. Antennas Propagat., vol. 40, pp. 727-731, June 1992. 21

[14] 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. [15] K. Barkeshli and J.L. Volakis, "Electromagnetic scattering by an aperture formed by a rectangular cavity recessed in a ground plane," J. Electromagnetic Waves Applications, vol. 5, no. 7, pp. 715-734, 1991. [16] S.M. Rao, D.R. Wilton and A.W. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 409-418, May 1982. [17] D.R. Wilton, S.M. Rao, A.W. Glisson, D.H. Schaubert, O.M. Al-Bundak and C.M. Butler, "Potential integrals of uniform and line source distributions on polygonal and polyhedral domains," IEEE Trans. Antennas Propagat., vol. AP-32, pp. 276-281, March 1984. [18] J.M. Jin and J.L. Volakis, "Radiation and scattering from microstrip patch antennas via a hybrid finite element method," International Conference on Electromagnetics in Aerospace Appl., Politecnico di Torino, Italy, Sept. 1991, Conference Proceedings pp. 195-198. [19] D.G. Shively, M.C. Bailey, C.R. Cockrell, M.D. Deshpande, "Scattering from microstrip patch antennas using subdomain basis functions," Electromagnetics, 1993. [20] R.B. Green, "The general theory of antenna scattering," The Ohio State Univ., Ph.D. dissertation, Nov. 1963 (also Ohio State Univ. ElectroScience Laboratory Report #1223-17). [21] T.K. Sarka, "On the application of the Generalized Biconjugate Gradient method," J. Electromagnetic Waves Appl., vol. 1, no. 3, pp. 223-242, 1987. [22] G. Golub and C.F. Van Loan, Matrix Computations, Johns Hopkins Univ. Press, 1983. [23] J.L. Volakis and K. Barkeshli, "Applications of the Conjugate GradientFFT method to radiation and scattering," in Application of Iterative 22

Methods to Electromagnetic and Signal Processing, ed. T. Sarkar, Elsevier Publ. Co., pp. 159-240, 1991. List of Figures I It figure Al. Local geometry of the rectangular brick. 'igure B1. Edge and vertex numbering for the tetrahedron. Figure 1. Antenna geometry. Figure 2. Geometries for various finite element shapes. Figure 3. Illustration of Sap and Scab. Figure 4. Computer code flow chart. Figure 5. Comparison of measured and calculated RCS patterns for the shown rectangular patch [11]. Figure 6. Comparison of measured and calculated backscatter RCS data for the shown circular patch as a function of frequency (i0 = 60~, 4i = 180~). Figure 7. RCS of a 0.5 cm x 0.5 cm patch in a rectangular cavity 1 cm x 1 cm x 0.5 cm. Comparison of calculations using rectangular bricks and tetrahedra. Figure 8. Relative normal incidence RCS/Gain Reduction in dB of the shown rectangular patch at f = 1.97 GHz as a function of the four lumped loads placed at the edges of the patch. The cavity is 0.17558 cm deep and is filled by a substrate having a relative dielectric constant of 0.17. Figure 9. RCS of the patch and cavity shown in Figure 2 as a function of frequency. Figure 10. Relative normal incidence RCS/Gain of the unloaded antenna configuration in Figure 2 where the patch is coated with RAM having E, = 10.2 - j3.8 and O, = 2.12 - jl.5. Calculations are referred to the resonant frequency of the coated patch. I 23

Figure 11. Geometry of a patch in a cavity. The patch is loaded with a square resistive ribbon having a resistivity of R = 0.5Zo. Figure 12. aeg backscatter RCS of the patch in Figure 11 as a function of frequency with and without resistive skirt loading; incident wave direction Oi = 70~, Oi = 180~. Figure 13. a6e backscatter RCS of the patch in Figure 11 as a function of angle with and without resistive loading. Figure 14. Spiral configuration with small cavity aperture. Figure 15. RCS of the spiral configuration in Figure 14. Figure 16. Radiation power pattern of the spiral antenna shown in Figure 14. Figure 17. Comparison of computed and measured input impedance for a loaded microstrip patch antenna. (a) Real part. (b) Reactive part. Figure 18. Comparison of input impedance computations for a circular patch. The patch radius is 1.3 cm and resides in a cavity of radius 2.115 cm and 0.406 cm deep filled with a dielectric of c,. 2.9. The feed is 0.8 cm from the center of the patch and the frequency is swept from 3 GHz to 3.8 GHz. 24

Cavity (a) Plan view Patches / \ Region I stripline / (b) Cross-sectional (xz) view Figure 1. Antenna geometry. 25

I AI Skewed Triangular Prism Skewed/Distorted Brick/Six sided Prism Tetrahedron Right Angled Triangular Prism Right Angled Brick Figure 2. Geometries for various finite element shapes. Cavity 5ap I.- ground plane coaxial cable Scab Figure 3. Illustration of Sap and Scab. 26

IMPLEMENTATION FLOW CHART I I I I I I I I I I I I I Figure 4. Computer code flow chart. 27

y c I X T = 0.057 in r = 4 f = 9.2 GHz -I ---— i 2.89 in- -10. -20: -30. -40. -50. with patch 0. 15. 30. 45. 60. 75. 90. 0 (degrees) Figure 5. Comparison of measured and calculated RCS patterns for the shown rectangular patch [11]. 28

0.000 I -- I S I - ~ ~~~ - -7 C Ad -5.00 / - cm, inCrm X|5cm -10.0 n FO - 5 cm -- o0 Location of Loods -15.0 / I x -5 / --- Normalized Gain in dB, --- Normalized RCS in dB. No load at feed -20.0 j ----- Normalized Rcs in dB with 48.92+j 1.13 Ohms at feed v (matched load for 300 Ohms) -25.0 - < 0 300 600. 9(0) 1200 1500 Load Resistance in Ohms Figure 8. Relative normal incidence RCS/Gain Reduction in dB of the shown rectangular patch at f = 1.97 GHz as a function of the four lumped loads placed at the edges of the patch. The cavity is 0.17558 cm deep and is filled by a substrate having a relative dielectric constant of 0.17. -10.0 -20.0 -30.0 -40.0 -50.0 E.~5 U/ cr.1 I' T 1 r I I I I __ RJ Unloaded patch 1.80 Four 300 Ohms resistors at the patch edges Four 300 Ohms resistors at the patch edges and 48.92+j1.13 Ohms at feed (matched load at 1.97GHz) 2.00 2.20 -60.0 1 -70.0 - 1.00 ~,,1.... I.... I.... I...!........ -...-I 2.00 3.00 4.00 5.00 6.00 Frequency in GH2 7.00 8.00 9.00 10.00 Figure 9. R.CS of the patch and cavity shown in Figure 2 as a function of frequency. 30

0.000 -1.00 - \ - Normalized RCS in dB -2.00 -----—..- Normalized Gain in dB -3.00 \.. -4.00 \ -5.00 -.. -6.00 RAM Coating -7.00. 0. 17558cm -8.00 - ' _ -9.00 -10.0 0.00 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.30 Thickness of the RAM coating (t) in cm Figure 10. Relative normal incidence RCS/Gain of the unloaded antenna configuration in Figure 2 where the patch is coated with RAM having cr = 10.2 - j3.8 and,r = 2.12- jl.5. Calculations are referred to the resonant frequency of the coated patch. 31

9375cm 9.375cm Y (0,0) x Resistive skirt 0.15625cm wide, R=0.5 Dielectric er=2.17 pr=1.0 / Lower Left corner of patch at (x,y)=(2.1875cm,2.96875cm) Patch dimensions 5.00cm x 3.4375cm [j] Feed at (x,y) = (3.4375cm,3.75cm) Figure 11. Geometry of a patch in a cavity. The patch is loaded with a square resistive ribbon having a resistivit y of R = 0.5Zo. 32

-10.0 -30.0 ' U -40.0 -30.0 - '40. 0," / _ No R-skirt around patch -50.0 -50.0 ------- R-skirt around patch 0.15625 cm wide R = 0.5 -60.0.... 1.50 10 1.70 1.80 1.90 2.00 -60.0 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00 Frequency in GHz Figure 12. ase backscatter RCS of the patch in Figure 11 as a function of frequency with and without resistive skirt loading; incident wave direction Hi = 70~, i = 180~. 33

0.000 i............ |........|.... i.... -10.0 -20.0 - U -30.0 - 0 No Rskirt at resonant f=1.9522GHz —.... 0=O With Rskirt at resonant f=1.744GHz -40.0.... 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00 0 in degrees Figure 13. ais backscatter RCS of the patch in Figure 11 as a function of angle with and without resistive loading. 34

Frequency 4.5 GHz Coordinates of cross (X at lower right corner of conducting strip): (1.406 cm, 1.406 cm) z~: 0.09375 cm I I I 2.8125 cm y+ 0.09375 cm - - - I x. --- —--------- I I J L I w I x 2.8125 cm Spiral m______ mmm mm without absorber 9375 cm AIR Spiral AIR x e,= 1.1-j0.5 ~r= 1.5-jl.0, = =2.0-j4.0 \ \\\\\\\\\\\\, \ \, with absorber I + N 0.9375 cm 0.6cm 0.6 cm 0.6 cm 18 dB Absorber at 4.5 GHz 7 N Figure 14. Spiral configuration with small cavity aperture. 35

-20.0........ I.... -.. |.. |. 30.0 - E; --- —- ---- ----- _ --- —-------- -40.0 Spiral -50.0 -—. --- — Spiral+Absorber.. ---- Spiral+Absorber+Rcard over the aperture (R=5.0): ------ Spiral+Absorber+Rcard over the aperture (R=1.0) -60.0. 60.00 65.00 70.00 75.00 80.00 85.00 90.00 8 in degrees Figure 15. RCS of the spiral configuration in Figure 14. 50.0........... 40.0 'c, 230.0 f. --- —------— ". —......^-.- - | 20.0 - "' — 10.0. - O o.ooo ` --- —----- ------------- ~ = 0.000 - -- aa -10.0^ -20.0 Spiral Sp Zin=109-j340 ',, -20.0 Ce -------- Spiral+Absorber Zin=78-j437 -30.0 ------ Spiral+Absorber+Rcard (R=5.0) Zin=61-j389 -40.0 - --- Spiral+Absorber+Rcard (R=1.0) Zin=31-j395 50.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00 0.0000 0.00 20.0 30.00 40.00 50.00 60.00 70.00 80.00 90.00 0 in degrees Figure 16. Radiation power pattern of the spiral antenna shown in Figure 14.

. - tn E r - o0 C-^8 cr C. nr (xf,yf) = (0.85 cm,1.22 cm) (x,y) = (-1.5 cm,-2.2 cm) T=0.08779 cm F (GHz) (a),E c -4 o~ i - I - cr=2.17 tan6=0.0015 ZL = 50 ohms 0.0 I.1 If-~%.V t I ' ' ' 'I I I 'r' Ts I I - I I I ii I I T0 I I t5 2.0 2.5 3.0 3. F (GHz) (b) Figure 17. Comparison of computed and measured input impedance for a loaded microstrip patch antenna. (a) Real part. (b) Reactive part. 37

*ZHD 8' o~ ZHD ) moij Idaxs st:(uanbalj ail pue tplld all jo Jaluao aIll moj mU3 8'0 si pauj utll '6'S m j jo 3!.3Iajip 6e i!AM pallm daap Um3 90P*0 pue mU3 9'II- snipel jo XA!^)e3 e uI saplsaI pue ul 1g)' si snip1j i31ed atIL 'plsd Iejn3Jp e Joj suo!llndtuoD aiuepadum! Indui jo uosr!edumoD 8T ajnrj 0'I O 7 -0.1 --- 0 - O'l W t 90W0. -.............................:::: i::::::: i::

Appendix A: Matrix elements for rectangular bricks Referring to Figure Al, it is convenient to first introduce the definitions {Et} = {~L,;J,;J} {WT} = {Nj:yN } i 3 3 3 (Al) where j = 1,2,3,4. Thus, for example, Ee = 522, E- = y3, We = xNx2, W[ = 'Ny3, etc. The field expansions can then be written as Ex = Nxj(y, z) O; j=1 4 Ey= Nyj(z x)l y; j=1 4 Ez = E Nzj(x, y) Oj j=l (A2) where Ne = (b-y')(c-z) T y'(c-z'). Nae - (b-y')z. x1 be ' ~2 bc ' 3 be Nx4 =- c ' =4 - bc; N e - (c-z)(a-x). yl ~ ca I ATe _ (a-x')(b-y'). z1 - ab I Nye _ z'(a-x'). y2 - ca; Ne _ x(b-y') z2 - ab I re c-z(xt y re ZX. N 3 - ca 4 - ca; Ne _ (a-x')y Ne N z3 - ab ' z4 ab (A3) In these, (x',y',z') denote the local coordinates specifying a point within the eth element and from an examination of the expansion functions we observe that O<- represents an average of the field component EX along the edge segment (1,2). Likewise, O}2 is associated with the Ex component along the edge (3,4) and so on. The elements of the 12 x 12 [Ae] matrix defined in (10) can now be readily evaluated in a straightforward manner since Npj (p = x, y, z) are simple linear functions. To do so, for notational and programming convenience, we shall denote these elements as A(pj),(qk), where p = x,y,z; q = x,y, z; j = 1,2, 3, 4; and k = 1,2,3,4. We next define the matrices 'abc [K0] = -36 4 22 1 24 1 2 2 1 4 2 1 224 39

2 -2 1 -1 [/fA1 1- 1 2 -2 [IA] = 7 9-1 -1 1 -2 2 '2 1 -2 -1 r]- 1 2 -1 -2 IAz- — 2 1 2 1 -1 -2 1 2 2 1 -2 -1 -2 1 2 1 [i3] = 1 2 -1 -2 -1 -' 1 2 where a, b and c are the side lengths of the brick as defined in Figure Al. On using this notation we have [AXjzk] [Azj,yk] [Axz]zi [Azj,xk] [Azj,yk] [Azj,zk] I ac[1 b 2a [Azj,2k] = [Ir'] +$ 6c[2]-ko 2?r[IJ} [Axj,xk] = [ +{ 6a [ K2] kCr r [Ko] c b [Azj,yk] [K3= - ], [AKzk] = []T 6/Ir 6[,. c a [Ayj,3k]- [K3]T [Ayjzk] = [a3] 6/r 6/1-' 37{ 6a 21 + [I b a T [Azj,k] = - [K3], [Azj,yk] - [K3] G6L G6/, [A,j,zk],= -L f K1 + -K2] k- t[Ko]' If a resistive sheet of resistivity R borders a face of the eth brick, then Bi. will be non-zero. Denoting the field components at the edges of the subject rectangular face by E~, Ei2, E1, and E2, it follows from (14) that e~l~i = B2,i2 B~e~j1 j, ko Zo iil, 2 - i2,i2 Bjjl = 22- lo d B BR 40

BkoZo i1,i2 i2,il = B1,j2 2 = 6R didj where di and dj are the corresponding side lengths of the rectangle and R is the (unnormalized) resistivity of the sheet. All other elements of the [B] matrix are zero. 41

Appendix B: Evaluation of matrix elements for tetrahedrals Referring to fig. B1 and the associated table, the fields in the eth tetrahedron are expanded as 6 E= EEW i=1 where the basis functions We are given by We -(r) =- f-i + g7-1 X r r E Ve 7-i = O0 outside element b7-i f7- = - r; x ri, ri1, ri position vectors of vertices 6'~e?'~il and i2 (see Table) bib7-i g7- eV (ri9 - ri ) ei -- tbi bi = Ir2 - ri 1 = length of the ith edge (see Table) V, = element's volume Wle note that V * We = 0 V x W = 2gi indicating that W? are divergenceless. Furthermore, W () j - ij ~ 0 i j where rj has its tip on the jth edge of the tetrahedron. This last property ensures that the coefficients Ee = E ei represent the average field value at the ith edge of the tetrahedron. Using the above basis functions, we now proceed with the derivation of the matrix elements Ae. We have (v x wi) e(7 x We) = -gi.gj /Jf J J r 42

Also, er JW WWjdv = er f {(fi, fj)+(rD)+(gi xr).(gj xr)} dv = r(I1 + 2 + 13) where D = (fi x gj) + (fj x gi) and = fi 'fj'iv '2 = /vr D d 13 = /J (gi x r)(gj x r) dv Since f is a constanit vector, 11 reduces to II= f~ fjVe To evaluate 12 we first set 4 x = L — i=1 4 y = EZLY i=l 4 z = Lizi i=l where xi, yi, z: (i = 1,..., 4) denote the (x, y, z) coordinates of the tetrahedron's vertices and Li are the simplex coordinates or shape functions for the same element. That is, Li is the normalized area of the tetrahedron formed by its three corners other than the ith, and the point (x, y, z) located within the tetrahedron. Using the standard formula for volume integration within a tetrahedral element and simplifying, we have 12 Z 4 4 4 I2 = + D+ E+D y + D E i=1 i=1 i=1 43

where Dm is the mnth component of D. The evaluation of 13 can be simplified by the use of basic vector identities. We have, 13 = gi gjj rl2 dv - gv(i r)(g. r) dv = (giygjy + gizgjz) v x2 dv + (gixgjx + gizgjz) v y2 dv + (gixYgx + giygjy) v z2 dv - (gixgjy + gjx9iy) 'v xy dv - (gixjz + gjixgiz) zx dv - (iySgi + Sjysi.-) J yz dv where gi,, represents the 77th component of the vector gi. Each of the above integrals can be easily evaluated analytically and the result can be expressed in the general form 1K al,, dv = -0 aliami + al ami 20=l =l i=l for 1, 7 = 1,..., 3. The parameters al or am can represent any of the rectilinear variables x, y,. 44

Appendix C: Biconjugate gradient algorithm for solving the symmetric system Ax = b Initialize the residual and search vectors with an initial guess xo: po = o= b - Axo Iterate for k = 0, 1,2,...: ak = (rk, rk) (pk, Apk) Xk+l = Xk + akPk rk+l = rk - CkApk = (rk+l,rk+l) (rk,rk) Pk+l = rk+l + fkPk Terminate when err = 7rk+l 7 < tolerance. Ilbll In the algorithm, (x, y) = xTy, where T denotes the transpose of the column. 45

Appendix D: Implementation with tetrahedral elements In the case of tetrahedral elements, the surface basis functions Sj(r) become the same as those given by Rao, etc. [18] and are best defined by grouping the triangles in pairs. For computational efficiency (i.e. to allow use of the FFT) and without loss of generality, it is necessary that the aperture elements on Sap are chosen to be right triangles. This leads to the faceting shown in Figure D1, where the new indices (m, n) denote the grid element location. It is clear from Figure D1 that the new numbering scheme makes the three different vector surface basis Sj1'2'3 to be identical for each right triangle. Two of these are associated with triangle edges which are orthogonal and are directed along x and y, respectively. The third edge makes an angle depending on the size Ax and Ay with respect to the x or y direction. Because of the surface basis uniformity (which is a result of the geometrical uniformity), the boundary integral matrix resulting from the discretization of Fs is block Toeplitz, which is the primary advantage of the aperture faceting using right triangles. Obviously, all symmetry properties are with respect to the (m, n) indices and it is therefore instructive to redefine the surface basis functions in terms of (m, n). With this in mind, we introduce the definition SJ(r) = S(z, y), j =1,2,3 (1) where S12n'3 are the expansion functions associated with the three different classes of edges (triangle pairs) mentioned earlier. Each of these edges is shared by a (+) and a (-) triangle as illustrated in Figure D1. We shall denote these triangles as T,+ and T7 where (m, n) refers to the left lower node of the rectangle which is composed by these triangles. From [18] or Appendix B, we have that Sm(z, y) 1 (nAy - y)x + (x- mAx)y r E T+ = A j A [y - (n + l)Ay]x + [(m + 2)x - x]y r E Th+ln (2) S zA, A)ny sL,(^)t 46

I2 (nAy - y)x + [x- (m + l)Ax+]y r E T+ [y - (n + l)Ay]x + (mAx - x) r E Tr,, (3) 0 otherwise SMn(X, Y) I [ [(n + 2)Ay -y] + [x- (m + l)Ax]y r E TZ,.+ = 3AA (y - nAy) + (mnAx-) r E Tn (4) O otherwise where lj denote the length of the edge shared by the triangle pair. With the above definitions in place, the [G] matrix elements are given by G. Gmn,m' n =-2ko2 J/ [Smn(r) z]. {JJ [S (r') x ] Go(r r') ds'} ds (5) +2 /Js V [Sj(r) x z] {J/J V ' [Simn(r') x z] Go(r, r') ds'} ds in which Sm, denote the surface area encompassed by the triangle pair. This expression can be simplified further by invoking the special properties of Sm (x, y) given in (2)-(4). Specifically, we find that V *[SL(r') x] = [V x S'm(r)] 2= - ej(r) (6) ArAy ej(r) = { r E Tm (7) mn Also, since (S x ) * (Sn, X Zs SMin', it follows that G',,ni, can be rewritten as -mnm~ sJJStn(r) i cG~,mn = -- l J n Jjs ~ StuM(r) Smn(r) Go(r, r') ds' ds + (A^ y)2 Si. Jll e(r) E'(r') Go(r, r') ds' ds (8) The excitation vector Qj is given by Qm? - -2jkoZo J S (H x ) ds +2 O n X (' Hr)ds s +2joZo/ x; ~ Z H'(r) d (9) n 47

where j = 1,2,3 as before. An explicit form of the boundary integral subsystem which is suitable for the application of the FFT is as follows: [G m',n-nl] [Gm-m',n-n'] [G mInn] E } - {Qmn} m —m,n-n,n-n [ GGm m,n-n] i {E n'} {Qmn} 31 32 33 { a n} { }Q3 m [G m] -m'n] [G-n' [mm,n-n] [Gm m {nmlnl {Q}mn (10) In this, we have used the important property that Gmn,m'n' = - -mm,n-n' (11) which is a consequence of the chosen symmetry among the surface elements. The columns { Em,, } refer to the unknown fields at the jth class edges whose lower right node is (n',n') as illustrated in Figure D1. In practice, the {EjA'} and {QJn} columns are constructed by incrementing a single index p = (m, n), where n(M-1) + (m +1) j =1 elements p = nM + (m + 1) j = 2 elements (12) (n - 1)M + (m + 1) j = 3 elements with n = 0,1,2,...,N and m = 0,1,..., N. The most important aspect of the system (11) is the convolutional nature of each of the nine submatrices comprising the [G] matrix. Each of these submatrices can be actually decomposed into a set of M x M or N x N small submatrices which are Toeplitz. Thus, [Gm,, nn, ] are block Toeplitz matrices which is again a basic property of any two-dimensional convolutional integral. The convolutional nature of these nine matrices can be exploited for computing the matrix products using (N In N) x (M In M) operations rather than N2 x M2 operations. In particular, we are interested in evaluating the sum AM-1 N-1 EC~Gj E l -" (13) m-E E mn-n m'n' Qmn (13) m'=O n'=O where Qmn is an M x N matrix which contributes to the final entry Qnm as described in (10). By carrying out the sum operation (13) for each (m,n) pair we will generate the column {Qn}j such that 3 E LQm j -= {Qmn} (14) j=l 48

Although it appears that the sum (13) can be computed via that FFT in a rather straightforward application of the discrete convolution theorem, special care must be exercised before use of the FFT. The discrete FFT operation can only be applied to periodic functions, implying that [G_m,nin',] must be circulant before its application for computing {Q-n}j in connection with the sum (13). Otherwise, aliasing cannot be avoided. The [Gtj] blocks can be recast into circulant matrices by extending their size from N x M to (2N - 1) x (2M - 1). That is, the original matrix blocks must be doubled along each linear dimension (the memory size is actually increased by a factor of four). The additional entries must be properly filled to make it circulant. To achieve this, we first rewrite [Gij] in an expanded form as Gsj G... fij3(N U1 G0 -(N-2) G0 G... Gs(N3) =[Gmin-ni] (15) Gs2... =0 L _ N-1 N-2* G where each sub-block G^ is actually another M x M matrix which has the Toeplitz property noted earlier. Specifically [G<] = [G,,,,k] (16) We also remark that Gk - Gik, otherwise G- Z G k, i.e. the matrix (15) is non-symmetric. With these properties in mind, we can readily proceed to rewrite [G6_m,,nnI, ] in circulant form. We will basically expand (15), so that each row of [Gm,, n-ni] can be obtained from the first by a simple shift operation. Specifically, we replace [Grmi n,n'] by the new extended matrix Go G-1 *. G-(N-1) G(N-1) G(N-2) *. G1 G1 Go *. G-(N-2) G(N-2) G(N-3) ** G2 G(N-1) G(N-2) *.' Go G-1 G-2.*. G-(N-1) G-(N-1) G(N-I) *.. G1 Go G-1 *. G-(N-2) G-1 G-2 **" G-(N-2) G-(N-1) G(N-1) *. Go 49

which now has (2N - 1) x (2N - 1) entries Gk. In addition, each of the original Gk Toeplitz matrices is recast in a circulant form by a similar procedure, thus extending its size from M x M to (2M - 1) x (2M - 1). This new form of [G_,,i n-n, ] is now circulant, allowing us to proceed with the application of the FFT. It is important to note that the new circulant matrix [Gm,, n-n, ] can be generated by computing only a single row. Thus, the storage requirement of each of these matrices is only (2M - 1) x (2N - 1). In practice, we store the unique matrix entries using the arrangement {G-(N-1), G-(N-2), G-(N-3), '', GO,... GN-1} Of course, only the first row of each Gk matrix need be stored. Using the above extended form of the [G-mi, ] matrix, the sum (13) can be rewritten as 2AI-1 2N-1 G'E E Gm E = (17) m1=O n'=O where the column {E~ljA where the column {Enm,,} has been padded with zeroes to extend its size from N x M to (2N - 1) x (2M - 1). The column {Q,,} can now be generated via the FFT operation FFT-1 {FFT2D{Gn} FFT2D{E } = Q (18) This operation can be introduced in the system (??) discussed in Section 3.2 of the report. Specifically, the overall system can be rewritten as 1A]EVEA In2 (19) E {EV } + {Qbv} {bA} [r 4] {EA} + { n - {Q 2} (19) with {Qmn} computed from 3 {Qmn} = Z{Qmn} j=1 3 = rFFT {FFT2{ n} FFT2D{Em } (20) j=1 50

The given form of the system (19) assumes that the mesh be constructed so that the top layer of surface elements consists of equal right triangles as shown in Figure D1. This implies that the column {E } satisfies the relation EA) _ 1A { T )2A IT { aA }T {En} = {{Emn}, {Emn}, Emn (21) However, it is not necessary to force the original mesh to conform to the cell arrangement shown in Figure D1. If the surface elements of the given mesh are not right triangles, to take advantage of the computational and storage efficiencies associated with the structured mesh in Figure D1, we can overlay on the original mesh the rectangular grid shown in Figure D1. The edges of the rectangular grid will be associated with the field columns {Em}A but the relationship (21) will no longer hold. This relation must be replaced by a new one which can be obtained by linearly interpolating between {EA} and {EJA}, both defined on the same surface. Specifically, (21) must be replaced by a relation of the form r { Emn} {E} - {E2A}) (22) where [C] is referred to as the connectivity matrix. It is a very sparse and of narrow bandwidth matrix and does not therefore increase the memory requirements to a degree of concern. As can be expected, the [C] matrix is not square. 51

Z' 87 8 7 4 3 1 ^2 b. — a -J Figure Al. Local geometry of the rectangular brick. 1 I I Edge Number 2 Q~ud0 0x 0I D 3,3 Figure B1. Edge and vertex numbering for the tetrahedron.

N Lx < > //// --- — 4 Ly 1 n:O m:O 1 * x 2.... M (a) Aperture faceting using rectangular triangles (2 (1) (m,n) (3) Plus triangle (2) (1) 7 (m,n) Minus triangle (b) Illustration of the (+) and (-) triangles (m,n) (m+l,n) edge class 2 (m,n+l) H edge class 3 (m,n) (m,n) / (m+l,n) edge class 1 (c) Illustration of the three different triangle pairs Figure D1. Aperture faceting using right triangles

UNIVERSITY OF MICHIGAN 3 90111 111115 0265 111 0423 3 9015 02651 0423