Report #028918-3-T Report Title: A Finite Element-Boundary Integral Method for Cavities in a Circular Cylinder Report Authors: Leo C. Kempel and J. L. Volakis Primary University Collaborator: John L. Volakis Volakis@um.cc.umich.edu Telephone: (313) 764-0500 Primary NASA-Ames Collaborator: Alex Woo woo@ra-next.arc.nasa.gov Telephone: (415) 604-6010 University Address: Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor MI 48109-2122 Date: December, 1992 Funds for the support of this study have been allocated in part by the NASA-Ames Research Center, Moffett Field, California, under interchange No. NCA2-653.

A Finite Element-Boundary Integral Method for Cavities in a Circular Cylinder Leo C. Kempel and John L. Volakis December 22, 1992 Abstract Conformal antenna arrays offer many cost and weight advantages over conventional antenna systems. However, due to a lack of rigorous mathematical, models for conformal antenna arrays, antenna designers result to measurements and planar antenna concepts for designing non-planar conformal antennas. Recently, we have found the finite element-boundary integral method to be very successful in modeling large planar arrays of arbitrary composition in a metallic plane. Herewith, we shall extend this formulation to conformal arrays on large metallic cylinders. In this report, we develop the mathematical formulation. In particular we discuss the shape functions, the resulting finite elements and the boundary integral equations, and the solution of the conformal finite element-boundary integral system. Some validation results are presented and we further show how this formulation can be applied with minimal computational and memory resources. 1

Contents 1 Introduction 3 2 FEM-BI Formulation 4 3 Vector Weight Functions 8 4 Dyadic Green's Function 10 5 Numerical Implementation 14 5.1 Matrix Assembly........................ 14 5.2 Excitation............................. 17 6 Engineering Quantities 20 6.1 Cavity Eigenvalues..................... 20 6.2 Far-field Evaluation........................ 21 7 Numerical Results 22 7.1 Resonant Wavenumbers..................... 23 7.2 Scattering......................... 24 8 Future Work 25 2

1 Introduction Conformal antenna arrays are attractive for aircraft, spacecraft, and land vehicle applications since these systems possess low weight, flexibility, and cost advantages over conventional antennas. The majority of previous developments in non-planar conformal antennas has been conducted experimentally due to a lack of rigorous analysis techniques. Various approximate analysis techniques are restricted in many respects, including accuracy and element shape, and are based on planar antenna models. Recently, we have found that the finite element-boundary integral (FEMBI) method can be successfully employed for the analysis of large planar arrays of arbitrary composition [1]. The resulting system is sparse due to the local nature of the finite element method whereas the boundary integral is convolutional, thus ensuring an O(N) memory demand for the entire system. In this report we will extend the FEM-BI formulation for aperture antennas conformal to a cylindrical metallic surface. Both the radiation and scattering problems will be developed in the context of the FEM-BI method. In contrast to the planar aperture array, the implementation of the cylindrically conformal array requires shell-shaped elements rather than bricks, and the required external Green's function is that of the circular perfectly conducting cylinder. In its exact form this Green's function is an infinite series which must be evaluated efficiently and must also be put in a convenient convolutional form for storage minimization. This report presents the FEM-BI formulation, the appropriate cylindrical shell elements, and the system evaluation strategy which will maintain low memory and computational load. The cylindrical elements will be chosen divergenceless while maintaining excellent geometrical fidelity. These elements will be analogous to the bricks used by Jin and Volakis [1]. Efficient asymptotic evaluations of the cylindrical Green's function will be discussed. In addition, the numerical implementation of the FEM-BI method will be presented along with some numerical results for validation purposes. These results will illustrate the accuracy of this formulation for the calculation of cavity eigenvalues and Radar Cross Section (RCS) calculations. Antenna pattern calculations will be presented in a future report. 3

2 FEM-BI Formulation Consider the configuration illustrated in figure 1 where a cavity is recessed in an infinite,circular metallic cylinder. The cavity walls are assumed to coincide with either constant p-,O- or z-planes and we allow the possibility of radiating elements on the surface of or within the inhomogeneous substrate which fills the cavity. The FEM-BI formulation [1] allows the determination of the electric field present in the cavity due to interior or exterior sources. This method utilizes the finite element method to formulate the interior field and yield a sparse system of equations. This system is coupled to a second set of equations generated by enforcing a boundary integral equation on the cavity aperture. Although in general the boundary integral system of equations if full, a judicious choice of boundary elements allows 0(N) storage to be maintained. The scattered or radiated field is readily calculated from the aperture fields. The development of the FEM-BI system begins with the vector wave equation valid in the interior of the cavity where we also allow the possibility of interior electric (Jint) and/or magnetic (Mint) sources. Specifically, we have V x x(p,, z) - koe(p, i, Z) (p, q, ) = -j koZoJint (p, q, z) a + x, (1) rP, (r[P, Z) relative permittivity and permeability of the substrate, ko is the free-space system of equations from (1) we apply the method of weighted residuals J, {V x [ E(P, < z) -k2r(p, >, Z)E(p, q, z). Wi(p, q, z)}pdp d dz = 4

' Vx { V X p -,j, zo int,, z, z)p dp do dz Jvi y Lpr(P r,Z) Z (2) where Wi(p, X, z) are subdomain vector weight functions to be specified. The substrate is discretized using cylindrical-shell volume elements such as the one shown in figure 2. The vector weight function coefficients represent the field at the edges of these elements and the integration volume (Vi) corresponds to the elements which possess the ith edge. The forcing function due to the interior impressed sources (Jint int) is given by ~ (p, 0, z ) ' Ain=t L V l [ "nt (P ps z) ] j kZ nt p A TV(p, $, z)pdp dq dz (3) We recognize (2) as the weak form of the wave equation and upon application of a standard vector identity and the divergence theorem [2] (2) becomes V fVxE( p,<, z) V x Wi(p, q, z) ~Jvi, ur (p,, z) -koer(p, f, z)E(p, O, z) Wi(p, q, z)}pdp dq dz -jkozo j n(p,,, z) x H(p, q, z) Wi(p, X, z)dS = ft (4) with nh(p, (, z) indicating the outward directed unit normal for the surfaces of the elements associated with the ith edge, Si is the total surface of those elements, and H(p, 4, z) is the total magnetic field. It can be shown that the surface integral in (4) vanishes for all surfaces except where the surface corresponds to the aperture. Thus, Si is the subdivision of the aperture surface associated with the elements which possess the ith edge. Unfortunately, (4) requires the determination of H(a, X, z) over the aperture which increases the required number of unknowns. To eliminate H(a, ), z) from (4) we use the integral representation of the magnetic field due to surface fields H(a, q, z) = H(a,,z) +H(a, z)+ 5

jYaka J p(a, )', z) x E(a,, )z) Go(a, - ), i)d(' dz +a j p(a, Q', z) x H(a, $, z') V x Go(a, ~, z)d$', dz (5) with q = -, z - z and S encompassing the entire surface of the cylinder. The first two terms are the incident field and the scattered field from the cylinder with the aperture removed. We shall denote these two terms by the symbol HCYl(a, (, z). The integral equation (5) cannot be readily coupled with (4) because it involves both the tangential electric and magnetic fields at the aperture. However, we are free to choose a dyadic Green's function which satisfies the radiation condition and the additional Neumann surface boundary condition V x Gm2(a, X, z)= 0. (6) With this condition, the second integral in (5) vanishes and the integration surface need only extend over the aperture surface. The dyadic Green's function which satisfies-(6) is usually termed the magnetic dyadic Green's function of the second kind [3]. We may now write (5) as H(a,,z) = H'(a,,, z) +jYokoa j p(a, O', z') x E(a, ', z') * Gm2(a, 4, z)dq' dz (7) a where Sa indicates integration over the entire aperture. Upon inserting (7) into (4), we obtain [V x E(p,, z) * V x Wi(p,, z) JE(1' ir(p, <, Z) -kor(p, ), z)E(p, A, z) * Wi(p, (, z)}p dp d dz +(ka)2 j [Wi(a, q, z) * p(a,, z) x Gm2(a, (, z) x p(a, ', z')E(a, >', z')] dq' dz'dq dz = ft + fez (8) in which a standard dyadic identity has been used and the forcing function due to exterior sources is fet = jZokoa j,(a,, z) p(a,, z) x 'HCY(a, ', z')d5' dz (9) i S9 6

We recognize that (8) has unknowns corresponding to the electric field within the cavity and on the aperture. Following the principles of Galerkin's method, the electric field in (8) is expanded in terms of the same vector subdomain weight functions as was used for testing, i.e., N E(p,, z) = EjWj(p,, z). (10) j=1 In this expansion (interior + aperture edges) N is the total number of unknowns or edge fields. Since the vector wave equation (1) requires the electric field to be divergenceless, the weighting functions used in (10) must also be divergenceless to avoid the use of a penalty term [1] which is required if this condition is not met (e.g. node-based elements [4]). By necessity the aperture field takes the form N E(a,, z) = Eja(j)W j(a, q, z) (11) j=1 where 6a(j) = 1 if j E aperture edge = 0 else (12) and the vector weighting function is the same as that used in (10) when evaluated on the surface. Combining (8),(10) and (11) we obtain the FEMBI system i JV{ <x j(p,i, z).Vx Wi(p,$, z) tr(P, O, Z) -k~o2 r(p,, z)l;Vj(p,, z). W1i(p, c, z) p dp do dz +(koa) a(j)6a(') / / [v&(a,, z) p(a, z) x Gm2(a,, Z) x p(a,, z) z 1j(p, z)] dQ dz dqdz =' f ent + ft(13) in which Si and Sj represent the surface elements bordered by the ith and jth edges, respectively. Below we discuss the specific vector weighting functions for the shell element shown in figure 2. 7

3 Vector Weight Functions Traditional finite elements associate the unknown field coefficients (Ej) with the node field. However, it has been found that such elements enforce nonphysical field continuity at element boundaries since both the tangential and normal field components are continuous which leads to spurious modes that must be suppressed with a penalty function [4]. However, edge-based elements have only tangential continuity and are therefore better suited for electromagnetics applications where jump discontinuities in the normal component is permitted. In addition, edge-based functions avoid an explicit specification of the fields at corners where the field may be singular [5]. Various volume element shapes have been used for finite element applications. Two of the most popular are: bricks [5] and tetrahedra [6, 7]. Although the later element is quite versatile, we shall develop a cylindrical shell edge-based elements since we require geometrical fidelity which is not available with bricks or tetrahedra. Bricks are the natural element for discretization of a volume which is defined by portions of the Cartesian planes. Likewise, cylindrical shell elements are superbly suited for cavities where the six faces are defined by portions of the constant p-, q-, and z-planes. Thus, the cylindrical shall elements will be analogous to brick elements used by Jin and Volakis [5]. By inspection, the shell element in figure 2 is comprised of twelve edges which will constitute the vector weight functions. These weight functions should satisfy the following properties: 1. subdomain (finite element) 2. |1 l'(p,,z) 11= 1 if (p,, z) E jth edge 3. j| Wj(p,, z) 11= 0 if (p,, z) E any edge I jth edge 4. V. Wj(p, ), z)=0 lWj(p, X, z) will only satisfy this requirement within the volume of the element. These weighting functions introduce artificial charges on the faces of the element and are not divergenceless at element interfaces. This is required since these elements do not guarantee normal field continuity across the element faces. 8

These conditions are met by the following set of weight vectors W12(p, >, Z) = W,(p, 0, z;, Or, Zt, +), W56(p, ), Z) = Wp(p, 4O, Z;, zr, Zb, -), W43(P,, Z) = Wp(p, ), z;., 1i, Zt, -) W87(p, 4), Z) = Wp(p, (, z;, z; l, Zb, +) Wl4(p, q, z) W58(p, ), z) = w(p, o, z; pb,, Zt,+), = W,(p, ), z; pb,, Zb,-), W23(p, 4, Z) = W (p, 4, Z; pa,., Zt, -) W67(p, 4, Z) = Wo(p, ~, Z; Pa,. Zb, +) W15(p, O, Z) =;lz(Pi, O, Z; Pb, Or,, +), W26(P, O4, Z) = /z(p, 4, z; Pa, )r, ', -) W48(p, 0, z) = Iz(p, (, Z; Ph, i, -,-), 1V37(p, O, Z) = Hz(p,, z; pa, 01,,,+) (14) where Wik denote the edge which is defined by local nodes (l,k) which is shown in figure 2. Each vector weight is represented by one of three fundamental vectors (one for each orthogonal direction) ~vp(pl Ol z; 7i,4),~l.)= VV, (p, ), Z;b,q5, i, 9) = Spa ( )- ()(Z -- ) ^ ah p (P - )( - ) (p - p)(z - z):A tc ^ (15) with t = Pb - Pa, a = q - )01, and h = t - Zb. The element parameters (Pa, Pb, )I, Or, Zb, Zt) are shown in figure 2. We observe that all four of the requirements listed above are met by (14) with the understanding that these weight vectors are non-zero only within the element. Interestingly, in order to satisfy the divergenceless requirement by forcing each component of the element field to be divergenceless, we need the p term in the p-directed weight (15). Of course, for very large radius cylinders and small elements, the curvature of these cylindrical shell elements decreases resulting in weight vectors which are functionally similar to the bricks used by Jin and Volakis [5]. We shall now specify the appropriate dyadic Green's function which will allow us, in conjunction with the vector weight functions given above, to fully specify the FEM-BI system (13). 9

4 Dyadic Green's Function The boundary integral present in (13) requires numerical computation and an efficient evaluation of the dyadic Green's function, Gm2(a, ), z-), is crucial. This dyadic Green's function may be derived in its eigenfunction form from a set of cylindrical vector wave functions using the procedure espoused by Tai [8]. The resulting expression contains an infinite sum of angular eigenvalues and an infinite integration over the axial eigenvalues. It is well known that the angular eigenvalue expansion for cylinders exhibits extremely poor convergence characteristics for large radius cylinders [9]. We must therefore resort to approximate evaluations of the Green's function which do not possess the prohibitive computational demand of the exact modal expressions. As described above, one may readily derive a dyadic Green's function which satisfies both the radiation condition (by virtue of the chosen wave functions) and the Neumann boundary condition on the surface of the cylinder (6). The resulting modal expression is Gm2(a, X, z) = q+'G'(a, z) + [qz + zX] GQZ(a, q, z) + zz (a, X, z) (16) where unprimed unit vectors are functions of the observation point (a, q, z) while primed unit vectors are functions of the source point (a, q', z ) and z z(' = (2)2Z L| ( ) - 1H(2) (')) 1 [0 (nk / H$2)(Q) j(nq-kz,f)~, (27r)2 no -. ( 2) (2) I-c (2Hn)( n=-oo) kako ) H(2)() k 00 kn~a — 2 oo (17) (17) with kp = k -k, 7y = ka, a is the radius of the cylinder, n denotes the angular mode number, and kz is the axial eigenvalue. As stated previously, these expressions exact a high computational demand especially when the source and observation point coalesce. We therefore want to recast (17) 10

I i in some other approximate form which improve as the radius of the cylinder increases and the distance between the source and observation points becomes large. Several research groups have developed in the past such evaluations of (17) and the most successful are attributed to Pathak [10], Boersma and Lee [11], and Bird [12, 13]. All three approximations utilize Watson's transformation to convert the slowly converging angular eigenvalue representation (17) into a rapidly converging series over the radial eigenvalue. The Hankel function and its derivative are approximated using either the Debye expansions or the uniform asymptotic expansion [14] and finally the axial eigenvalue integral is evaluated via the steepest decent method. These three formulas differ only in the level of approximation offered. The three formula are listed in tabular form by Bird [12] where he compared the accuracy of each. In developing these formula, the following quantities are necessary in describing the on-surface ray contributions: 1. Geodesic path length (s = /(aI) + 2) 2. Geodesic trajectory (0 = tan-1 [1a]) 3. Approximation order (q = __) and in these CD = 6 or I = 2rr - 6 depending on which path as illustrated in figure 3 (the short or long) is used. Boersma and Lee [11] employed approximate expressions from a Debye approximation of the Hankel functions and a first-order evaluation of the axial wavenumber integral to obtain GZZ(a, Z, ) -j qe-ko (cos2 + q( - q)(2 - 3co )) v/3) 27T ^ v / -q/3v (/3) [cos 20q + q + tan - cos 12 6 3 64cs 6O*(a ~,S) ~ jko GIz (a,q$, ) q 2e-Jke sin 0 cos 0 (1 -3q(1 -q)) v(/) +q/3v'(/3) [ta20 - + q (1 -sec0)]} -34 ~ 64 4 11

27r 3i v) G~+(a, ') -j ~ e- k t (sin20 + q(l-q)(2 - Sifl.)} +q [sec20 (u(3) -v(/3)) + 3 (v'(/) (tan20 11 si2) 3 \ 1 2 u (3)sec20 + v (/3)q ( - - 2tan20 + si20))] (18) Pathak's expressions are based on a uniform asymptotic expansions of the Hankel functions and a first-order evaluation of the axial wavenumber integral. His results are GZZ(a,,) ~ -jk2 e- jks{ (cosO2 + q( -q)(2 - 3cos20)) v() G'z (a, < Z) j~qe-Jk~ sin 0 cos 0 (1 - 3q(1-q)) v(3)} Jk0 -ks 2 2 G*(a,, z) - 2 qe- 7 (si2 + q(1 - q)(2 - 3sn2)) v(/3) +q [sec 0 (u(3)- ())] (19) Finally, Bird [12] used higher order terms from the uniform asymptotic formulas and a second-order evaluation of the integral to obtain GZZ(a,, ) -- e-jk (cos20 + q( - q)(2 -3cos20)) v(3) 31 5 (11 17.20 v() ko -jksi2o - - [(l324 ) v-l3) +' (60 -36i20) si120 + o )v2(1) + J' \(tl } GoZ(a, Z,) J qe-iko sin 0 cos 0 (1 - 3q(1- q))v(/3) +[q ( sec 20- v() + 36 - 45s8 20) l() 12

( 1 sec20 - ) V2(/) + 3s5ec2o() (3)] } \15 24 G(z(a, X 7) 2 qe-jk~ (sin2o + q(1 - q)(2 - 3sin20)) v(/) + q 2 o 8 20 _3i 20\V O) +q (u(3) - v($))sec20 + (9tan20- 7sin20) v(/3) + (sin20 - -tan20) V (/3) \36 45 + a - -ta 2 (si 20) v2(/3) +;5 tan20a(1) j) (20) 2 where 3 = ks [J~]5. In (18)-(20), u(3) and v n (/) represent the soft and hard surface Fock functions, respectively, while u'(/) and v'(/) denote the first derivative of the zeroth order surface Fock functions. These functions are characteristic of the creeping waves on a circular cylinder and are discussed in detail by Logan [15]. We will consider only the two direct path contributions from the creeping wave series which are shown in figure 3. Accordingly, the ray parameters given above depend on which direct ray contribution is being computed. Bird [12] has concluded that his expressions (20), which include terms up to O(q-2) are the most accurate throughout the full range of X and z except when s is very small or the observation point is in the paraxial region of the source. In this region, the formula given by Boersma and Lee (18) are most accurate since the Debye approximation used for the Hankel function is well suited for the paraxial region and O(q-3) terms have been retained. A notable exception to this conclusion occurs when the cylinder is sufficiently large while the aperture is small. In this situation, the aperture appears to be a hole in a planar metallic screen and the formula given by Pathak (19) turns out to be the most accurate. Inspection of (18)-(20) reveals all three approximate expressions include the free-space dyadic Green's function terms up to O(q-2) modulated by a surface Fock function whereas Bird's as well as Boersma and Lee's formula contain 'extra' terms. These terms are not present in the case of a metallic screen and thus provide reduced accuracy when the aperture looks planar. However, for cavity configurations where curvature is considered, Bird's conclusions should be valid. Therefore, 13

we shall use Boersma and Lee's formulas (18) if Icos l < 0.1 or s < 1A. Otherwise, Bird's formulas (20) will be used. The derivation of (18),(19) and (20) is described in [11],[10] and [13], respectively. 5 Numerical Implementation In the preceding sections of this report, we described the vector weight functions and dyadic Green's function appropriate for the solution of (13). Below we now explicitly develop the admittance matrix and excitation vector entries necessary to solve (13) for the electric field. Once the unknown field is determined, engineering quantities such as cavity eigenvalues, RCS and antenna radiation patterns may be determined. We may write the N equations implied by (13) in terms of N unknowns(Ej) as a matrix equation EP} [rap] fet [] J t [ [] ~] {E l f (21) where [A] is a very sparse N x N symmetric matrix, [9] is an Nap x Nap square full matrix, the Np x 1 column vector {E^P} denotes field values on the aperture edges, the N, x 1 column vector {Ent} is field values on interior edges, and the total number of unknowns (N) is the sum of the free aperture edges (Nap) and the free interior edges (N,). A free edge is any edge which is not tangential to and on a metallic wall since such an edge is known to be associated with a null field due to the natural boundary condition on the metallic wall. Hence, it is not included as a system unknown. The Nap x 1 column vector {f ext} is the forcing function attributed to an incident external field (9) or impressed sources on the aperture (3) while the N, x 1 column vector {fInt} is due to the impressed sources within the substrate (3). 5.1 Matrix Assembly The sparse nature of [A] is revealed upon inspection of (13). Since the edge vectors are non-zero only within the elements which possess that particular edge, the volume integral is identically zero unless both the source and testing 14

edges are within the same element. These matrix entries may be written as A.- 1 (1)ij -k2 (2)i (22) -3- -- st ~ ko r.st [Ir where we have assumed constant material properties within each element with Er and jIr denoting the relative permittivity and permeability of the element. The subscripts (ij) indicate the row and column of the matrix and are related to the global edge numbers. Since each edge may belong to more than one element, the matrix entry Aij is computed by adding the contributions for each element which have both the ith and jth edges. The auxiliary terms IJ)i" and ( )ij are defined by S(l)ij = j V x S(p(, z,; pj j) V x wt (p, a, z; pi, i) Zi, SZi)pdp dodz i(2)i= (t2)i =- J W,(p,, z; j, j, j, j) Wt(p,, z; Pi, i, i, i)p dp d dz (23) with (s,t) e {p, pq, z} indicating the direction of the source and test edges, respectively, while the weight vectors are given by (14). Since the fundamental weight vectors (15) have orthogonal orientations and (23) is symmetric with respect to the source and test vectors, (1) need be determined for six combinations of (s,t) and I2) need be determined for three such combinations. These integrals may be evaluated analytically upon performing the vector operations indicated in (23). In addition to [A], we must also evaluate the entries of [9] which includes all the boundary integral contribution to the FEM-BI system (13). The entries of [9] involve the integral Gij - (ka)2 / / Wt(a, q, z; pi, 4i, zi, si) [, z) x Gm(a,, )] p(a, q, z) x Gm2(a, qfz) x /(a, ', z) Wvs(a, O( z';pj, j, zj, Sj)d)' dzdd z d d(24) which unfortunately cannot be evaluated analytically due to the presence of the Green's function. From an inspection of the approximate Green's function formulas (18)-(20) we observe that Gm2(a,,z-) becomes singular 15

as the source and observation point coalesce. When these points are close, a reasonable approach for large radius cylinders is to replace Gm2(a,,z) with the metallic screen Green's function whose singularity has been treated previously in the context of this integral [4]. Accordingly, we use the following boundary integral when the source and test edges are either in the same or some adjacent element where Gij ( oa2 j a zpi i i) W s(a z;pj z dz dz 2~A -2~i j v.[p(a, X, z) x GW.(a, ), z; pi, i,, )] oV* [p(a:) z ) x wS(a,) z; pIj, jj,,j) ] -R dZ dzd~ dz (27) This form of the boundary integral may be readily evaluated even as R vanishes by employing the regularization procedure used by Jin and Volakis [4]. The asymptotic dyadic Green's function used here consists of two direct interactions as shown in figure 3. We may use the planar approximation (26) for the short path contribution when the distance between the source and test elements is small. The asymptotic formulas are always used for the long path contributions, i.e. interactions between distant elements when curvature effects are likely to be prominent. 16

Explicit formulas for the computation of the admittance matrix entries have been given above which may be readily implemented. The [A] matrix may be considered the finite element portion of the system while [9] is the boundary integral contribution. The matrix [A] is very sparse and its entries may be computed from analytical evaluations of the integrals (23). Although the [9] submatrix is generally full and its entries must be computed numerically using for example two-point Gaussian quadrature, considerable storage and computation economy is possible by utilizing a CG-FFT type of solution procedure. The convolutional kernel of the dyadic Green's function makes this possible (17). The computation of the forcing functions in (21) must now be addressed. 5.2 Excitation Two principal types of excitation are found in the FEM-BI system (21). One is the excitation due to internal sources (possibly on the aperture) when the radiation is being considered. The entries for the resulting column vector {fAnt} is given by (3) and require the specification of the impressed currents (Jint, Mint). The second type of excitation is due to external sources which are non-zero for scattering problem computations. This excitation column vector {fXt)} are given by (9) which requires the surface magnetic field (HCY') when no aperture is present. This later excitation will now be examined assuming plane wave incidence. Circular cylinder scattering is one of the few exact solutions available in electromagnetics. Traditionally, the field scattered by a cylinder due to plane wave incidence is expressed in modal form. Assuming that the incident field is given by E = 2e e-jk0(k^') Hi = Y(i x e ke-o(kt) = Yo ['i sin a cos Oi- Ci cos a - z sin a sin Oi] ek3[p in cos (q-qi)+zcos.i] (28) where (Oi,oi) indicate the direction of incidence, a is the polarization angle where e = cos caOt + sin acit, and Yo is the free-space intrinsic admittance. The total surface field HCYl(a,, Z) = Hi(a,, z) + Hcy,(a, 5, z) (29) 17

is given by sina ej oo o n(s Oz ) H'Y(,<M = zT)0 '' -- ------ ~jko cos Oi z COS HY< (a ), z) -2Yoron(2) (koa sin Oi) ka sin a cosi ] n(+ (3 koa sin i Hn(2)(koa sin i)e0 These expressions may be computed by summing only a few terms of the series if ka sin Oi is small. However, as this parameter becomes large (e.g. for large a and 0, -~ 90~), asymptotic evaluations similar to those used for the dyadic Green's function must be applied to (30) to maintain computational economy. Utilizing Watson's transformation and Fock theory [9] in connection with (30) we find 2 H? o -Yo sin a sin Oiejkocos'z e-jkoasinOi<p [g(~)(m p)]* p=l m 2 2 Y j2Yo cos c0 kza sin ic z [f(O)(m()]P ko-,a sin 0I:p=l p=i J ko a sin Oi (m )] (31) 1 * -j (') (m * (31) k Ia sin 0' in which 1 = 2 -(- ), 1 = (- r m = [ in and the '*' indicates complex conjugation. The appropriate Fock functions are [15]2 = J' J eddt 9 '[' = Jr wi (t) J fr WIet f(l)() = w(t) dt (32) 2Logan uses the e-jw time dependency in the definition of these functions requiring the complex conjugation in (31) 18

where the Fock form of the Airy functions may be written in the more familiar Miller notation wl (i) = v/ [Bi(t) + jAi(t)] and the integration contour is given by Logan [15]. The asymptotic formulas are quite accurate compared to (30) except when X,;i. In this region, Goriainov [16] derived the following more accurate expressions from Fock theory Hi1 -Yo sin a sin OiJ k cos { e-jkoasin Oi ( [g(o)(m i) )]* + ejka sin i cos (-i ) [G(-m cos ( i)- 4i))]* } H yl j2Yo cosa, e i22 oei~ cos~s iz e-jkasinibl [f(O)(m)l)]* " -jkO a sinicos(q —,i) [F(-_ cos ($ -- i))]* } +Yo sin a cos Oiejkcos iz{e-jkoasinOi )l[g(O)(ml) mJ, g(1)1 ((l) -J koa sin I _jkoasini cos (-i) [G(-m cos( - i)) -J m i G(1)(-m cos (z - i)] a } (33), t oa sin ei s ( 3) F() - f(')()e3 (34) These surface field expressions may now be used to efficiently calculate the entries of the column vector {f t} via a numerical evaluation of (9). In particular, the modal series (30) is used when koasinOi < 10 and either (31) or (33) for koa sin Oi > 10 as appropriate. We note that the unit vectors in (9) and (28) are not the same since the former depend on the integration variable (5') while the latter depend on the incidence angle (qj). Thus, the 19

excitation integral (9) is f?' = jZo(ka ) cos (' -,Oi)i(a, ', z') [iHY (a, ', z )-) ]H'(a, z)] dq' dz (35) 6 Engineering Quantities The FEM-BI system (21) has now been fully specified and either a direct or an iterative solution of (21) for the system unknowns {Ej} may be applied. Once the expansion coefficients have been determined, we may calculate useful engineering quantities such as cavity eigenvalues, RCS and radiation patterns. 6.1 Cavity Eigenvalues One of the uses of the Finite Element Method in electromagnetics which has become popular recently is the calculation of cavity resonances numerically. Often analytical determination of the cavity modes is not possible if the cavity has an odd shape or is filled with an inhomogeneous material. Chatterjee et al.[7] have shown that edge-based tetrahedral finite elements are well-suited for computing cavity eigenvalues. These elements are able to model a wide range of cavity shapes and material filling. Likewise, the cylindrical shell elements introduced in this report are ideal for the calculation of modes present in any cylindrical cavity which does not intersect the z-axis. For eigenvalue calculations, the cavity aperture is covered with metal which closes the finite element mesh, implying [9] = [0]. Therefore, the eigenvalues may be computed by setting the interior portion of the system (22) equal to zero and identifying the eigenvalues as A = k2 [~ i] {E} = A [(2)j] {E} (36) We recognize (36) as a generalized eigenvalue problem where A denote the eigenvalues of the cavity. These eigenvalues may be determined directly or by exploiting the inherent matrix sparsity, Lanczos method may be used. Such a solution of (36) yield the resonant fields {Ej} (eigenvectors) and corresponding wavenumbers (eigenvalues). 20

6.2 Far-field Evaluation Two of the most important applications of the formulation presented in this report is the calculation of the cavity's radar cross section (RCS) and the radiation pattern due to sources placed within or on the aperture of the cavity. In this section, we will describe far-field calculation due to equivalent sources on the aperture. This will entail a similar analysis as was performed for the excitation fields. We begin with the integral representation of the scattered magnetic field in terms of the aperture fields. We have, Hs(r, O, =) ji Yoka jGm2r,, a, a, z) [p(a,, z ) x E(a,:', z')] d ' do (37) with (r,0,0) indicating the observation point in spherical coordinates. When the observation point is very far from the cylinder, the dyadic Green's function has the asymptotic form -jkr[ m2(r, 0,; a,, z) k o [G'+ O o + Gz + G+ ] (38) where unprimed unit vectors are functions of the observation position while primed unit vectors are functions of the integration point in (37). The components of this asymptotic Green's function are determined by a mode matching procedure GB 3 roo 2kOkcos Jko sz' n n (+( (27r)2 (koa sin 0)2 ~ Hn(2)(ka sin A ) (27)2 a E=- H(2)( (koa sin 0) *. <-) 0 I ji jka cos O z (39) '(22sn H ^ ~ }{) G~z (27r)2 a 2) As one might expect, these series converge rather slowly for large koa sin 0. They must therefore be recast in another form by employing Watson's transformation and Fock theory as described before. We have, G t ]o cos kcos (_l) e-k asinp () )(mp)-j zl ( Goo cs47 Z Z()P( L k0a sin 0 ( (iP)J p —1 21

Gz - eiko os Oz' -Cjkoasin Op [9(o)(m )]* 47r p=] 2 2 G2ar sin2 jkocos' e-jkkoasinO<p [f(o)(m0p)] (4 2a7r sin0 P= where the appropriate Fock functions are given by (32). As expected in the region )' - Sq, the formulas attributed to Goriainov [16] k asicos(-^) eG(( COS (4 - 4)) -. m * 7C(- 1 4 jkanLo cos -o gO)( ) -i jka sin 0 J G^ - _ sn ekocos'- { ekoasino~1 [g(o)(m4i)]* + ekoka sincos(- ) [G() (-m cos ( - ') ) ] }m Go4 _e g e + 2aw sin [ ma~ jin;eek~~~5~Z ~-jkoasinOIl [f(O)(mI>l)]* + ejkoasin cos( [r(0)(-m cos ( - ' ))] } (41) are most useful. The far-zone scattered field can be computed numerically by using either (39),(40), or (41) in (37). For the scattering problem, RCS is most often the quantity of interest. Once the far-zone magnetic field (37) is computed, the RCS is given (0, ) - lim 4rr2 H(r, ) (42) r- oo IHi(r, 0, q) 7 Numerical Results The preceding section of this report presented several engineering quantities of interest: cavity eigenvalues (resonant modes), RCS and antenna pattern calculations. We will now present some preliminary results as a means of validating the numerical implementation. In particular, we shall compare 22

the numerically computed mode wavenumbers to the analytical solution for air-filled cavities. We shall also present some preliminary comparisons of RCS calculations with known results. Antenna pattern calculations will be deferred to a later report. 7.1 Resonant Wavenumbers An important cavity characteristic is its resonant frequencies. These may be determined analytically from the eigenvalue series solution of the wave equation for simple geometries. For example, if an air-filled rectangular cavity is considered, the resonant wavenumbers may be determined from the characteristic equation. Such a cavity is shown in figure 4a. Likewise, the resonant wavenumbers of a cylindrical sector cavity (see figure 4b) may be determined from its characteristic equation. The dimensions of the cylindrical sector cavity may be chosen so that this cavity is equivalent to a rectangular one which then requires that the resonant wavenumbers be identical. This was done for a equivalent rectangular cavity with dimensions: a = 0.25 cm, b = 0.5 cm, c = 1.0 cm. The following table illustrates the capability of the presented formulation to numerically compute the first five resonant wavenumbers mode exact computed error TEl11 5.236 5.306 1.34 TM1lo 7.025 7.441 5.92 TEoll 7.531 7.908 5.01 TE201 7.531 7.928 5.27 TM111 8.179 8.562 4.68 The mode number corresponds to the equivalent rectangular cavity and the same is true for the exact wavenumbers. The the computed wavenumbers were determined by solving the generalized eigenvalue problem (36) and we note that the percent error is comparable to the error reported using bricks to discretize a rectangular cavity [7]. A cylindrical sector cavity is shown in figure 4b. We have compared the exact resonant wavenumbers with computed results for a cavity where pb = 5 cm, pa = 4.75 cm, a = r, - 01 = 5~, zt = 0.5 cm, and Zb = 0.0 cm in the following table of the first five modes. 23

mode 11 exact computed error TEl11 9.695 9.856 1.66 TEl01 14.051 14.283 1.65 TE112 14.575 14.816 1.65 TM11o 14.575 15.539 6.61 TM111 15.872 16.134 1.65 Again the percent error is within the expected bounds. The exact resonant wavenumbers of a cylindrical shell (see figure 4c) may also be found if the cavity is filled with a homogeneous material. For a shell where pb = 1.0 cm, Pb = 0.9 cm, Zt = 1.0 cm and zb = 0.0 cm the comparison between the calculated and exact wavenumbers are given in the table below. mode | exact I computed J error TM10oo6.246 6.450 3.27 TM11o 6.393 6.557 2.57 TE111 6.428 6.608 2.80 TM120 6.814 6.713 1.51 TE121 6.831 6.733 1.46 7.2 Scattering The FEM-BI formulation has been shown to provide an excellent tool for computing the scattered field attributed to a cavity recessed in a ground plane [4]. This report presents an analogous method when the cavity is recessed in an infinite,metallic circular cylinder. We have shown that the finite element method may be used to numerically calculate interior resonances for closed cavities. We now wish to validate the proposed formulation for scattering calculations from an open cavity due to plane wave incidence. Unfortunately, we currently do not have a suitable experimental or numerical data to determine the accuracy of this method. However, for small cavities embedded in a large radius cylinder, we expect our results to compare quite favorably with a similar method which was formulated for planar structures [4]. Below, we perform such comparisons for air-filled cavities. Comparisons have been made between a cavity recessed in a cylinder such as the one shown in figure 1 and an equivalent cavity in a metallic plane. First we consider the cavity whose dimensions are w x I x d = 2A x 2A x 0.25A 24

cavity. The azimuthal width is given by w = aAZ, where a is the cylinder radius in wavelengths and Aq is the azimuthal angle subtended by the cavity. The axial length (1) and cavity depth (d) are normalized by the wavelength. The incident wave is given by (28) and we will consider only TM polarization where a = 0, and consequently E- = 0. Figure 5 illustrates the comparison of the scattered fields for such a cavity in a cylinder with a = 20A and a similar cavity in a ground plane for normal incidence (0i = 90~, 5i = 0~). Additional bistatic patterns are shown in figures 6 and 7 for oblique incidence with (0i = 45~, Oi = 0~) and (Oi = 20~, qs = 0~), respectively. In all of these three figures, the azimuth observation angle is q = 0~ and the elevation angle (0) is allowed to vary. Figure 8 is a backscatter calculation of the same cavity for fixed X = 0~ and varying 0. Finally, we repeat these calculations for a w x I x d = 2A x 0.5A x 0.25A cavity in figures 9-12. Clearly the two solutions are in excellent agreement. 8 Future Work In this report, we have presented the formulation and preliminary validation of a FEM-BI method suitable for cavities recessed in a circular cylinder. Some validation of this formulation was presented for the special case of a large radius cylinder. However, we are actively seeking validation data where the curvature of the structure has a prominent effect on the scattering. Such validation data will either be found by measurement or comparison with another numerical solution. Such a solution will most likely be for an groove in a cylinder which will necessitate our consideration of a long narrow cavity and normal incidence only. Comparisons with curved cavities will help in determining which asymptotic Green's function formula (18)-(20) is best suited for these calculations. It will also aid in the identification of the source causing the difficulties with the other principal polarization. Inaccuracies using Pathak's formula appears to be the cause of this difficulty. After suitable scattering validation has been accomplished and a final decision as to which Green's function formula is most accurate, we shall implement a host of features which are useful for antenna design. Metallic patches and resistive cards will be placed on the aperture and different feeds will be included. Since the mutual coupling between antenna array elements is expected to be crucial for conformal arrays, input impedance and mutual 25

impedance calculations will be also considered. Currently, we have implemented the formulation without exploiting the convolutional kernel of the boundary integral. In fact, we are using a sparse (direct) matrix solver requiring the full boundary integral matrix be stored in its entirety. Such a direct solver allows moderately large geometries without the need to resolve the system for each new incidence angle. However, for large problems, the boundary integral must be evaluated using a Fast Fourier Transform (FFT) in conjunction with an iterative solver such as the BiConjugate Gradient (BiCG) method. The analogous FEM-BI formulation for planar structures [4] uses such a solver with impressive results both in terms of accuracy and efficiency. 26

References [1] Jian-Ming Jin and John L. Volakis, "A Hybrid Finite Element Method for Scattering and Radiation by Microstrip Patch Antennas and Arrays Residing in a Cavity," IEEE Trans. Antennas and Propagat., Vol. 39, No. 11, pp. 1598-1604, Nov. 1991. [2] Jian-Ming Jin, John L. Volakis, and Jeffrey D. Collins, "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, June 1991. [3] Chen-To Tai, "On the Eigen-function Expansion of Dyadic Green's Functions," University of Michigan Technical Report No. 011136-1-T, 1973. [4] Jian-Ming Jin and John L. Volakis, "A Finite Element-Boundary Integral Formulation for Scattering by Three-Dimensional Cavity-Backed Apertures," IEEE Trans. Antennas and Propagat., Vol. 39, No. 1, pp. 97-104, Jan. 1991. [5] Jian-Ming Jin and John L. Volakis, "Electromagnetic Scattering by and Transmission Through a Three-Dimensional Slot in a Thick Conducting Plane," IEEE Trans. Antennas and Propagat., Vol. 39, No. 4, pp. 543 -550, Apr. 1991. [6] M.L. Barton and Z.J. Cendes, "New vector finite elements for threedimensional magnetic field computation," J. Appl. Phys., Vol. 61, No. 8, pp. 3919-3921, Apr. 1987. [7] A. Chatterjee, 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. [8] Chen-To Tai, Dyadic Green's Functions in Electromagnetic Theory. International Textbook Co., Scranton, 1971. [9] 0. Einarsson, R.E. Kleinman, P. Laurin, and P.L.E. Uslenghi, "Studies in Radar Cross Sections L - Diffraction and Scattering by Regular Bodies 27

IV: The Circular Cylinder," University of Michigan Technical Report No. 7133-3-T, 1966. [10] P.H. Pathak and N.N. Wang, "An analysis of the mutual coupling between antennas on a smooth convex surface," Ohio State Univ. ElectroScience Lab., Report 784583-7, Oct. 1978. [11] J. Boersma and S.W. Lee, "Surface field due to a magnetic dipole on a cylinder: Asymptotic expansions of exact solution," Univ. Illinois Electromagnetics Lab., Report 78-17, 1978. [12] Trevor S. Bird, "Comparison of Asymptotic Solutions for the Surface Field Excited by a Magnetic Dipole on a Cylinder," IEEE Trans. Antennas and Propagat., Vol. 32, No. 11, pp. 1237-1244, Nov. 1984. [13] Trevor S. Bird, "Accurate Asymptotic Solution for the Surface Field Due to Apertures in a Conducting Cylinder," IEEE Trans. Antennas and Propagat., Vol. 33, No. 10, pp. 1108-1117, Oct. 1985. [14] J.J. Bowman, T.B.A. Senior, and P.L.E. Uslenghi, Electromagnetic and Acoustic Scattering by Simple Shapes. Hemisphere Publishing, New York, 1987. [15] N. A. Logan, "General research in diffraction theory," Lockheed Aircraft Corp., Missiles and Space Div., vol. 1 and 2, Report LMSD-288088, Dec. 1959. [16] A.S. Goriainov, "An Asymptotic Solution of the Problem of Diffraction of a Plane Electromagnetic Wave by a Conducting Cylinder," Radio Engr. and Electr. Phys., Vol. 3, pp. 23-39, 1958. (English translation of Radiotekhnica i Elektronica, Vol. 3) 28

List of Figures Figure 1. Illustration of the cavity geometry situated on a metallic cylinder. Figure 2. Cylindrical shell element. Figure 3. Geodesic paths on a circular cylinder. Figure 4. Cavity geometries for eigenvalue computation. a) Rectangular cavity. b) Circular cylinder sector cavity. c) Concentric circular cylinder cavity. Figure 5. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder for Eo polarization and Oi = 90~, qi = = = 0~ compared with an equivalent cavity recessed in a ground plane. Figure 6. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder for TM polarization and Oi = 45~, qi = ( = 0~ compared with an equivalent cavity recessed in a ground plane. Figure 7. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder for TM polarization and Oi = 20~, Oi = ) = 0~ compared with an equivalent cavity recessed in a ground plane. Figure 8. Backscatter scattering by a 2A x 2A x 0.25A cavity in a cylinder for TM polarization and 0 = 0~ compared with an equivalent cavity recessed in a ground plane. Figure 9. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and Oi = 90~, qi = < = 0~ compared with an equivalent cavity recessed in a ground plane. 1

Figure 10. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and Oi = 45~, qi = ) = 0~ compared with an equivalent cavity recessed in a ground plane. Figure 11. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and Oi = 20~, Oi = s = 0~ compared with an equivalent cavity recessed in a ground plane. Figure 12. Backscatter scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and q = 0~ compared with an equivalent cavity recessed in a ground plane. 2

y 0i x Figure 1. Illustration of the cavity geometry situated on a metallic cylinder.

4 I e I -4 -C! (a. E v) Q3 c"3 0) v *C 7= u C) To 0) Cz bO 414 N N

geodesic paths Figure 3. Geodesic paths on a circular cylinder.

z ~ —a —*/ /1 x a) z i abI -- _ t l b) ff _. j /^^^^^ ZX ~C) Figure 4. Cavity geometries for eigenvalue computation. a) Rectangular cavity. b) Circular cylinder sector cavity. c) Concentric circular cylinder cavity.

30.0 20.0 Cq 10.0 -10.0........ 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 5. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder for Eg polarization and Oi = 90~, - = - = 0~ compared with an equivalent cavity recessed in a ground plane.

3 0.0 ~..................I.........I............................................. Cylinder: a = 20 o Planar 20.0 10.0 "0 1I ~0.0 -10.0 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 6. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder for TM polarization and Oi = 45~, qi = ( = 0~ compared with an equivalent cavity recessed in a ground plane.

30.0 20.0 10.0 -10.0.................. |......... |......... L....... I.................. 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 7. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder for TM polarization and 0, = 20~, i- = 0 = 0~ compared with an equivalent cavity recessed in a ground plane.

3 0.0..................................................................... 20.0 C'q F 10.0 "U 0.0 Cylinder: a = 20X 0 Planar -1 0.0........................I.........I.................................... 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 8. Backscatter scattering by a 2A x 2A x 0.25A cavity in a cylinder for TM polarization and - = 0~ compared with an equivalent cavity recessed in a ground plane.

20.0......... I.........I.....I.........I.....I.........I.............I.....I......... 1 -~~~~~~~~~~~~~ r ^ -^ ^ ^ ^ ~"`l""`l"" * ^^ k-" Cl 04 rA. V. 15.0 10.0 - 5.0 I...1......I....1...I...1.....I......1. I I1,,,,!....... I.........I... 0.c I 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 9. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and Oi = 900, Oi = - = 00 compared with an equivalent cavity recessed in a ground plane.

20.0......I.........I 11.........I.........I......... I......... I......... I.... r-l Cl cn I:, Q 04 15.0 10.0 5.0 ~~~~~~~~~~~ ~~ ~ ~~~~~~~~~~~~~~~~= Cylinder: a = 20X o Planar I........ I....I.........I.....I.... I.a..I......... 0.0 I 0.0 10.0 20.0.30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 10. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and 0i = 45~, i - = = 0 compared with an equivalent cavity recessed in a ground plane.

20.0 "l. I.........I.........I....... IIIIIIII '1..........................I "0 N%00 i (A.T 15.0 10.0 5.0 -111 0.0 I ~.lz...... I...........I.........I.........I_ I.......1.I....1...... I 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 11. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and Oi = 20~, qi = - = 00 compared with an equivalent cavity recessed in a ground plane.

30.0 lt3 t-4 cr: +4-4 m u 033 20.0 I 10.0 0.0 Cylinder: a = 20X 0 Planar I i I.... I...... I....I.......... -10.0 I. I a...1.~ t......I.......m..llll llllllllI Illllllill1~~1~~l~ l~~~ ~ ~lll i l............. Bunti i B.................................. 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Observation Angle (0) [deg] Figure 12. Backscatter scattering by a 2A x 0.5A x 0.25A cavity in a cylinder for TM polarization and q = 00 compared with an equivalent cavity recessed in a ground plane.

UNIVERSITY OF MICHIGAN 3 9015 02651 0415