030601-1-T Radiation and scattering by cavity-backed antennas on a circular cylinder Leo C. Kempel and John L. Volakis July 1993 30601-1-T = RL-2422

030601-1-T TECHNICAL REPORT for NASA Grant NAG-1-1478 NASA Technical Monitor: Fred Beck GCD, Mail Stop 490 NASA Langley Research Center Hampton, VA 23681-0001 Telephone: (804) 864-1829 Grant Title: Analysis of Conformal Antennas in Composite Surfaces Radiation and Scattering by Cavity-Backed Antennas on a Circular Cylinder Report Title: Institution: Period Covered: Report Authors: Principal Investigator: Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor MI 48109-2122 February 1993 - August 1993 Leo C. Kempel and John L. Volakis John L. Volakis Telephone: (313) 764-0500 Volakis@um.cc.umich.edu

Radiation and Scattering by Cavity-backed Antennas on a Circular Cylinder Leo C. Kempel and John L. Volakis July 21, 1993 Abstract Conformal arrays are popular antennas for aircraft and missile platforms due to their inherent low weight and drag properties. However, to date there has been a dearth of rigorous analytical or numerical solutions to aid the designer. In fact, it has been common practice to use limited measurements and planar approximations in designing such non-planar antennas. In this report, we extend the finite element-boundary integral method to scattering and radiation by cavity-backed structures in an infinite,metallic cylinder. In particular, we discuss the formulation specifics such as weight functions, dyadic Green's function, implementation details and particular difficulties inherent to cylindrical structures. Special care is taken to ensure that the resulting computer program has low memory demand and minimal computational requirements. In this report, both scattering and radiation parameters are computed and validated as much as possible.

Contents 1 Introduction 4 2 FEM-BI for Circular Cylinders 5 2.1 Vector Weight Functions..................... 7 2.2 FEM Matrix.............................. 10 2.3 Boundary Integral Matrix..................... 12 3 Excitation 16 4 System Solution 21 5 Scattering and Radiation 23 5.1 Far-field Evaluation......................... 23 5.2 Input Impedance.......................... 25 6 Validation 26 6.1 Scattering............................. 26 6.2 Radiation................................ 33 7 Future Work 36 A BiCG Algorithm 40 B FFT-based Matrix-Vector Product 43 B.1 2-D Integral Equations...................... 43 B.2 3-D Integral Equations...................... 52 C Replication Rule (Fortran) 61 2

List of Figures 1 Cavity geometry............. 2 Cylindrical shell element........... 3 4 5 6 7 8 9 10 11 12 A-1 B-1 B-2 B-3 B-4 B-5 B-6 B-7 Geodesic paths................... Probe-fed conformal patch antenna....... Test body...................... Planar patch.................... Comparison with 2-D results............. Far-zone subtraction scheme........... Wraparound cavity scattering........... Antenna pattern for various cylinder radii.... Effect of cavity termination on antenna pattern. Input impedance................... BiCG convergence behavior............ Flat resistive strip geometry............ Operation count-2D case............. Circular arc geometry.............. 3-D geometries................... Rectangular patch discretization......... Local edge numbering scheme........... Operation count-3D case............. 6 9 15 20 27 29 30 31 32 34 35 37 41 44 49 51 53 56 57 60 3

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 protruding antennas. The majority of previous studies in non-planar conformal antennas has been conducted experimentally due to a lack of rigorous analysis techniques. Some approximate analyses have been considered but these are restricted in accuracy and element shape. Recently, the finite element-boundary integral (FEM-BI) method was successfully employed for the analysis of large cavity-backed planar arrays [1]. The resulting system is sparse due to the local nature of the finite element method whereas the boundary integral sub-matrix is fully populated. However, by resorting to an iterative solver such as the Biconjugate Gradient (BiCG) method, the boundary integral system may be cast in circulant form allowing the use of the Fast Fourier Transform (FFT) in performing the matrix-vector products. This BiCG-FFT solution scheme ensures O(N) memory demand for the entire FENI-BI system and minimizes the computational requirements. In this report we extend the FEM-BI formulation to 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 must satisfy the boundary conditions on the cylinder. In its exact form, this Green's function is an infinite series which imposes unacceptable computational burden on the method. However, for large radius cylinders, suitable asymptotic formula are available and herein used for an efficient evaluation of the Green's function. In addition, the resulting BI system is again cast in circulant form to ensure an 0(N) memory demand and take advantage of the FFT's speed when carrying out the matrix-vector product. At the end of this report, we also present computations of the relevant antenna and scattering parameters validated with available known results. One of the most interesting applications of the numerical laboratory developed in this project will be the ability to study curvature effects of real world antenna systems in an exact manner. Such studies will be pursued in the near future and the will be addressed in a future report. 4

2 FEM-BI for Circular Cylinders Consider a cavity recessed in an infinite,metallic cylinder, shown in figure 1. The cavity walls are assumed to coincide with constant p-, q- and z-surfaces and the cavity is assumed to be filled with inhomogeneous material. Interior sources and lumped loads may be present as well as surface metallization patterns and resistive cards. The FEM-BI technique [1] permits determination of the electric field within the cavity which are induced by either interior or exterior sources. The system of equations associated with the FEM formulation is sparse whereas the boundary integral sub-matrix is usually fully populated. However, with a judicious choice of boundary elements, the formulation will maintain O(N) memory demand when coupled with a Biconjugate Gradient-Fast Fourier Transform (BiCG-FFT) solver. Upon determination of the fields, the radiated and scattered patterns are readily calculated from the aperture fields while the input impedance, S-parameters and other pertinent antenna quantities may be computed from the appropriate interior fields. The FEM-BI system is developed directly from the inhomogeneous vector wave equation and a complete presentation of the derivation is given by Volakis et al. [2]. The resulting system is / fV x vj (p,q)V x Wz(p,, z) Vi, l /i (Pn > Z) -Or(p, Z)(p ),pz). ( )pd d d +(koa) 2a(j)6 a(i) J [fVfz(a )', )*p(a, z)x G2(a,, X ) x p(a, o z) * Wj(a, O', z ) do d' dz'd dz = t + f (1) where Wl;(p, c, z) are edge-based expansion functions where support is over the volume Vi and G2(a, X, z) is the pertinent dyadic Green's function. The free-space wavenumber is denoted by ko = 2, the cylinder has radius a whereas,. and Air are the relative constitutive parameters of the material filled cavity. Also V, represents the ith integration volume which extends over the support of WVi(p,, z) and in a similar fashion Si and Sj indicate surface integration over all aperture elements which extend over the support 5

0 Al *........ 'e R X -- XXR.,R.....z -.... -...',..........1.y'., OX'.E::':"''" ' ''.:::::.:....::..:::::::::-'::r': ~ ~ ~ ~ ~ ~ ~...::::::::':..:::...:::::::::::, ~~~~~~.-...:..:..........:.,::\..:':S:::.: ':'::':::'::'::'': '............... E':'::: ' - ~~~.''. -,..........::::...::::::i:,.:.:-'::::::::,:.::::: -:::.:'::-.:::.................................... 1:Ilstain.fte.aiy.em tr.iuae.n..eali yln e (a) Coordinat...... (b).yia.....ba ke......... 6

of the ith (test) and jth source basis functions, respectively. The function aL(i)ba(j) is a product of two Kronecker functions and it simply indicates that the boundary integral only contributes when both the test and source unknowns are on the aperture. We remark that the dyadic Green's function is convolutional (- =>- o, z = z- z ) and will be discussed later along with the functionals f/nt and fixt which depend on the interior and exterior excitations, respectively. Rather, we shall first proceed with a presentation of the appropriate weight functions and the evaluation of the FEM matrix entries which are represented by Aij in the FEM-BI matrix system [ + [9 [0]1 {E } (2) ~{Ej} - {f L} to be solved via the BiCG method. 2.1 Vector Weight Functions Traditional node-based finite elements associate the degrees of freedom (Ej) with the nodal fields. These elements impose both tangential and normal field continuity along inter-element boundaries even at material boundaries [3]. These elements also do not correctly model the null space of the curl operator and as a result spurious solutions are generated [3] which are typically suppressed by imposing a penalty function [4]. Edge-based elements, where the unknowns are associated with the free edges of the mesh, do not enforce normal field continuity (although tangential continuity is still maintained) and are therefore more suitable for electromagnetics applications. In addition, edge-based elements avoid explicit specification of the fields at geometry corners where the field may be singular [5]. Two popular volume elements are the brick [5] and the tetrahedron [6, 7] which are readily mated with rectangular and triangular surface elements, respectively. The brick is well-suited for geometries delimited by constant x-, y- and z-planes in the Cartesian system while the tetrahedra are more versatile (and consequently more complex). For cavities residing in a circular cylinder it is advantageous to employ the cylindrical shell element. Cylindrical shell elements, which are defined by constant p-, 0- and zsurfaces, are superbly suited for cavities recessed in a circular cylinder since 7

they possess both geometrical fidelity as well as simplicity. These shell elements are analogous to the bricks used by Jin and Volakis [5] and belong in the general class of curvilinear brick elements. Figure 2 illustrates a typical shell element which has eight nodes and twelve edges: four edges aligned along each of the three orthogonal directions of the cylindrical coordinate system. The edge-based weight functions for these elements must satisfy the following properties: 1. subdomain (finite element) 2. WI lj(p, 0, ) 11= 0 if (p, 4, z) E any edge ~ j and I jth edge 3. V. II(P,, ~)= 01 Such expansion/weight functions can be represented as WIV2(pi, z) = W'lp(p, O, z;, 0, t, +), WV43(p, O, z) = W,(p, (, z;.,, t, -) W.56(,, p,) = Pl' (p, 0,; *, r, zb,-), V87(p, Z) = Wp(p, ),;, Or, Zb, +) WI14(p,, Z) = I(p,,; Pb, 'I,t, +), W123(p, 0, ) = WO(p, ),Z; pa,, Zt,-) 158(p, =, Z) = Wl/"(p, ), Z; Pb,, Zb, -), V67(p,, ) = /(p, i), Z; pa,, Zb, +) IV1.5(P, O, Z) = L-'V(p, ), Z; Pb, Or,, +), /V26(p, q, z) = Wz(p, O, Z; Pa, r, r, -) W148(P, z) lV(p,, z; Pb, )i, ', -), W37(p 0X, z) = /zV(p, r, z; Pa, l1,, +) (3) where WIk denote the edge which is defined by local nodes (l,k) as shown in figure 2. That is three fundamental vector functions are required for the complete representation of the edge-based cylindrical shell element. They are given by SPb (O - ~ )(Z- ), ah P lWj (p,, 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

z z I II I I I Zb-0 y x Figure 2: Cylindrical shell element. 9

Wo (p, o; p 0 ) (P th \4 (p, Iz';p, Izs) = -(P-P)(+-+)z (4) where t = Pb - pa, a = r- (fi, and h = 2t - Zb and the element parameters (pa pb, 01, ri, Zbr Zt) are shown in figure 2. We remark that the 1 term which appears in the definition of the '-directed weight is essential in satisfying the divergenceless property of the edge-based expansions. 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]. 2.2 FEM Matrix The inherent locality of a partial differential equation formulation such as the FEM suggests that [A] is sparse. This FEM matrix results from the discretization of the first two integrals of (1) and owing to the finite support of the basis functions its entries are identically zero unless both the test and source edges are within an element. These matrix entries can be represented as A2 =- 1(i)ij - k22 (2)3 (5) SLi- st 7 St st( Ailr where we have assumed constant material properties within an element (Cr and Ji,) and the subscripts (ij) refer to the row and column of the [A] matrix, respectively. The two auxiliary functions are defined by (ij = V x S(pt 0, ^ Pin ~win)* V x Vtjp, fzpi,(fi,,)pdpdodz 1(2)itj=J /5(<2)tJ = AV S(P,5 ^pjO -jj s ' V t(p,,z; i, i,- 9 s)pdpdfdZ (6) where (s,t) e {p,o,z} indicate the direction of the source and test edges, respectively. Since the fundamental weight vectors (4) are aligned along orthogonal directions, (6) is symmetric with respect to the source and test edges (e.g. (1)i = I)3), and thus, only six combinations of (s,t) for I(1) need be determined and only three such combinations are required for I(2). They are 10

I(1) pp S- S S t 2h In Pb Pb (ah)2. Pa a 9 (P - 1) p~ jOr - )(o - 4t)do + j (z - 3)( - t )d z i(l) p(i) pz P Sl [ I (Pb' ( Pba 1 ft - bIn (P + (I - J ( )(Z - -.)z Pb - r - t [ (Pb -P) + (Ps + Pt) (p3 -p2 + P2" ( 2 2)) (t h.)2. 4 3t (2 ( - Pa) - _t (Ps + Pt) + PsPt In (-)d = - \t P (P- Ps)(p- Pt)dp t J pa,J P(1) j(l) Iz FZ -zt-Z~z z(z - )(z- It)dz Zb + (Pa) S.,,.z t h (to )2 (cah) SsStP I[ (- (pb - p) -t (s + pt) + t In t Pb- Pa) j - s)(o - 2 / \ O>0r f zt b In (b - )}(4t- )dI (Z 2 Pa} Jb [ ( P) + (I +p.i) (p - P) + (z - z)(z -tdz b - Z')(z - Zt)dz I(2) J(2) zz ')Ps Pt P 2 x - Pa - [ ( - ) -~ P ~ P ( P P ) (o- (s) (o- qt)dq l + IPsPt (Pb -Pa) ] x (7) where each of the unevaluated integrals are of the form j ( - ) (g -) dg =! (L2 - U2) ( + ) + (U3- L3) + g.,t (U - L) 2 3 - (8) 11

The [A] matrix is assembled by evaluating (5) for each element combination which contain the ith and jth unknowns. Since the integrals (6) are symmetric, only the upper triangle or lower triangle of the FEM matrix need be stored. 2.3 Boundary Integral Matrix The entries of the boundary integral sub-matrix, [5], couple the interior fields of the FEM system and any external impressed fields. In addition, the boundary integral provides an exact boundary condition for mesh closure. The entries of [9] are Gj (kL= a o )2 J t(a pi i ) p(a, f, ) x G(a,(, ) x p(a, ~ ts(a. X, z; pj, j, )j, sj)d' dz d.dz (9) and we note that the global nature of the Green's function implies a fully populated sub-matrix (e.g. all boundary unknowns couple to all the other boundary unknowns). The surface weight functions in (9) are the volume weight functions (3) collapsed down to a cylindrical-rectangular patch on the surface of the cylinder, p = a. The dyadic Green's function G2(a, 5, z) is of the second kind and it satisfies the Neumann boundary condition x x G2(a, -) =0. (10) on the cylinder's surface as well as the Sommerfeld radiation condition. It can be expressed component-wise as G2(a,,) = 'GO(a,X, ) + [~ + i'] G'z(a,,X) + izzG(a,, ) (11) where the unprimed unit vectors are functions of the test point (a,4,z) and the primed unit vectors are functions of the source point (a,,z ). The components may be exactly expressed as an angular eigenvalue series [8] 1,s) - (2 2 1 H (2) GZz (a, (, jr (P) 2,2) - n=-e ( kz )d 12

1 ~ n )_ H (2(-) i(n kzidk Coz(a, 6 =) - -- -F (^\Hw dk, (27r)2 J 2: y () ( (I r ) = 1 I~ [ ffH(2)() ( nkz 2 H(2) ei(n-kz)dk G`0(a, -J y_ -^ [ k ^ e' dkz (12) However, for large radius cylinders, (12) is computationally prohibitive. Instead, it is advantageous to employ one of the available asymptotic approximations of G2(a,, ) [9, 10, 11, 12]. These approximations improve as the cylinder radius and the geodesic distance between the source and observation point increase. Since each approximation invokes Watson's transformation, thus converting the angular eigenvalue series of (12) into a creeping wave series, the geodesic or on-surface ray parameters are necessary in describing the interaction between surface fields. The formula due to Pathak and Wang [9] have proven most accurate for our applications since the regions where Bird's formula [12] and Boersma and Lee's expressions [10] are superior (paraxial and far-zone) correspond to situations where coupling is small. In addition, we shall find it necessary to employ a regularization process to correct for errors in comnputing the near region (s -O 0) and to also allow evaluation of the singular integrand as (<,,) -) (O', z') which mitigates any advantage that the more complex formula have with respect to Pathak's expressions in the near-zone. Using a uniform expansion of the Hankel functions in (12) and a first-order evaluation of the axial wavenumber integral, Pathak found that Gzz (a,~) I - qeiks (cos2 + q(1 - q)(2 - 3cos2O)) v(3)} G~z(ao) jqqeikssinOS0cos0{(l-3q(1 -q))v(3)} 27r 1 J Gq(a,e~{) -s 2 + q(1 - q)(2 - 3sin2)) v(13) +q [sec20 (u(3)-v(3))] } (13) where 3 = ks [C ] 3 and the geodesic parameters are given by 13

1. Geodesic path length (s = (a) + 2) 2. Geodesic trajectory (0 = tan-1 [i]) 3. Approximation order (q = -) where - = O or ( = 27r - depending on which of the two direct paths illustrated in figure 3 (the short or long) is used. In (13), u(f) and vm(fd) represent the soft and hard surface Fock functions, respectively. These functions are characteristic of the creeping waves on a circular cylinder and are discussed in detail by Logan [14] along with their evaluation. As stated previously, a regularization process is necessary for the accurate evaluation of (9) especially when the source point and test point coalesce. Bird has shown [12] that Pathak's approximation (13) reduces to the metallic screen Green's function, 2Go(a, G, - ), modulated by v(3) within the available order (O(q-2)). This suggests a regularization procedure where the metallic screen Green's function is subtracted from the cylindrical Green's function, thus rendering it non-singular, followed by an addition of the planar contribution. The non-singular Green's function is given by G(a, -), j _COS Gzz(a, 1) - -o + q(1 - )( - 3cos20))[v(/)-1] G(a, ) qejsincos0 (l-:3q(1 - q))[v(3) - 1] T} Gk ( O,- 2 jk2 -K 2 0 G- -) { qe - sin2 + q(1 -q)(2-3sin20)) [v()- 1 +q [sec0 (u(3)-v(3))] } (14) and this asymptotic form of the Green's function is used for all short path contributions whereas (13) is used for the long path contribution in which case curvature effects are expected to be prominent. When the regularized integrand is used, the contribution of the planar Green's function GP = 2(ka)2 1 Vt(a,, XZ; pin Oii,, i). zJ J 5, JSJ 14

I geodesic -- -; - paths S (...... z). - I Fgr. Goei. pat.h. o a c... Figure 3: Geodesic paths on a circular cylinder. 15

where G0(a,(> ) [7 71~o e)jkpa,,R ) G,,)= [ +,2 4 (16) is added to (14) Gj= Gj +GP (17) which is used instead of (9) with R = s and I = xx + yy + zz. Upon use of a common vector identity and the divergence theorem [4], we obtain from (15) G (kta2 [ Wt(a, @ z; P.jzi, i)* 27 iS, is, 11W5(a0, (f;p;,s p,;,;)~sdo cz d4(l) -ai I, [ v *p(a pz) x (ta, 2 Pii, i Si)\ e'kR v [pa z )x W.(a, z pj o,S\ R dodz'di6dz (18) and 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]. We note that v(3) -+ 1 as a -> oo and hence the regularized integrand (14) vanishes leaving only the planar contribution (18). With the specification of both the FEM matrix, [A], and the boundary integral sub-matrix, [9], the unknown fields can be determined for a given excitation. In the next section, explicit evaluations of the excitation functions appropriate for plane wave scattering and probe-fed antenna applications are given. 3 Excitation Two sources of electromagnetic fields are considered: external sources (plane waves) for scattering analysis and internal sources (probe-feed) for antenna 16

parameter calculations. The use of the exact boundary condition in (1) allows the coupling of an exterior excitation field into the cavity while the FEI formulation itself [2] readily permits modeling of interior sources. In this sectiol we describe the form of the source functionals fer' and ftfl along with their numerical implementation. The forcing functional, feLt due to exterior sources is given by 'f, =jZ~ [a 0 za) (a, 2')z ) x HCy(a,z )do' dz (19) where I'-(p, O, z) is the testing vector for the ith row of the matrix and Hc1 represents is the magnetic field on the cylinder's surface in the absence of the cavity. For plane wave incidence i - ( ie-jko(k'r) H =,(k x eC)6c-Sj() = Yo [p' sin 7 cos 0, - ' cos 7 - sin i sin 0] e3jk[Psn co(-)+zcos s] (20) where (0i,Oi) indicate the direction of incidence, 7 is the polarization angle and ei = 0cos y + O'sin y. The total surface field is given by Hc,((a, z ) = H(a, z, z) + Hl(a,, z) (21) with H'-( a, O, )= j21'o sin 7ej koo, 0,= z. e(2 ] 7koa,=-Ho 2)(koasini) H\ = -2Y eoco sO* [ cos, + Hkoa'~n a~ =- I( koasin ki) n sin cos0 1 e(+-) (22) 3,a sin i H2)( koa sin O i) These expressions may be computed by summing only a few terms of the series if koasin O is small. However, as this parameter becomes large (e.g. for large a and i0 -- 90~), (22) may be replaced with equivalent asymptotic representations similar to those considered earlier. Utilizing Watson's 17

transformation and Fock theory [13] in connection with (22), we find that 2 Yo sin'cosz - asin s iekcsz ejkoip g( ) (0() H"- -Yosinh sinie^^^^ e-^0'"110^ [0 p=l 2 2 H~,CYI 7 j2Ykoos 7 l eJ-koasinop [f(o)(mI)]l j 2o cos', eo,a sin p=l -Yo sin y cos ieJkocosO'z Z(-l)pe-J kasin0P aop(mp) p=1 -i m g(1)(7l(I) p)] (23) -ka sin 0i (23) in which 3 = (-i), 2 = - ) - m = [koa i i, and the- '' symbol indicates complex conjugation. The appropriate Fock functions are [14]2 = j dt g'(E)= T JF W{ (t) f ( = l( dt (24) where wi(t) and its derivative wl(t) denote the Airy functions and the integration contour is given by Logan [14]. The asymptotic formulas (23) are quite accurate except when q O i. In this region, Goriainov's [15] expressions HCy " -Yo sinacsin Oieiocos'z e-{jkOsinll'1 [g(~)(mil)]x +jkoa sin 0 cos( ) (- s ( - ))]' } H ~u "j2Y kcos a.-ka sin 0f{m Y +ekoa sin cos( I-,) [F(-rn cos (o -,))]' } 2Logan [14] uses the e-jt time dependency in the definition of these functions requiring the complex conjugation in (23) 18

+1'O ncsina c kCos iekoc0 z{e-jkoa si" [g(~O)( ) -J- -— 1)1)] ka sin Oi - jkoa sin, cos (0-,) [G(-r cos (-Ji m G(1)(-mcos ( - i)]} (25) with [14] c-( ) _ = g )()e3 F () f()()ei (26) have been found to be more accurate and can be used instead of (23). These surface field expressions may now be used to efficiently calculate the entries of the column vector {fext} via a numerical evaluation of (19). In particular, the modal series (22) is used when koa sin 0 < 10 and either (23) or (25) for kAasinOi > 10 as appropriate. Conformal antenna patches are typically fed by a microstrip line printed along with the radiator on the surface of the substrate The microstrip lines are in turn fed by a coaxial probe which originates behind the cavity as shown in figure 4. The patch may also be fed directly by the coax feed or through some form of aperture coupling. Nevertheless, the principal excitation for the system is given by fil = - x [ + )] +jkoZoint(p,z) VV(p,,z)pdpddz (27) where Jfnt and MInt are the impressed electric or magnetic current densities representing the sources. For a radially (p) directed probe feed, the impressed monopole current located at (0, zs) is given by Jfnt = po-(X- s)(Z-ZS) (28) P

Patch or microstrip line / Cavity Metallic cylinder Coax feed Figure 4: Probe-fed conformal patch antenna. 20

which results in an excitation function (27) fji = -jkoZ Io In (-) [( -ci) (z-i))] (29) ci hi \ Pa With both types of excitation and the FEM-BI matrix now specified, the BiCG method may be used to determine the unknown electric fields within the cavity. 4 System Solution The FEM-BI system (2) may be solved using one of the popular direct methods such as Gauss-Jordan elimination, Gaussian elimination and LU decomposition. Alternatively an iterative methods such as Hotelling's method, conjugate gradient method or the biconjugate gradient method may be used. We have chosen to use the symmetric form of the BiCG solver because it requires much less memory than a direct method and more importantly it can be implemented in a manner which is computationally efficient (utilizing only one matrix-vector product per iteration). Although use of an iterative method such as the BiCG method can require more wallclock time than a direct solver when multiple right-hand sides are considered (such as is the case with backscatter calculations), memory demand was deemed the most critical and expensive resource. The BiCG algorithm is given in Appendix A and the efficient FFT-based calculation of the boundary integral matrixvector product is discussed in Appendix B. In this section, we will present details of this iterative solver specific to this application. The BiCG algorithm requires one matrix-vector product per iteration as shown in Appendix A. This operation represents the bulk of the computational demand of the method and requires O(N2) complex operations in the case of non-symmetric matrices. The matrix-vector product is carried out by executing the sum N y[n] = [A]{x} = A[n,n']x[n'] n = 1,2,3,...,N (30) n'=1 and if the matrix is sparse, a storage scheme such as the Compressed Sparse Row (CSR) format may be used to reduce both the memory demand and 21

computational load. Using the CSR scheme, (30) can be rewritten as r[n] y[n] = [A] {x} = ~ A[e[n, t]]rx[n'] n = 1,2, 3,...,N (31) n = 1 where r[n] is the number of non-zero entries per row of the matrix and e[n, n'] indicates which entry of the long data vector, A, is associated with the matrix entry A[n, n']. The FEM matrix, [A] in (2) is such a sparse matrix. Although additional memory savings are possible since the FEM matrix is symmetric (thus only upper triangle of the matrix need be stored), experience has shown that use of a symmetric matrix-vector product leads to a severe performance degradation on vector computers (such as a CRAY) due to the resulting short vector lengths. Therefore, the entire sparse FEM matrix is stored and used in the product so that the code's performance is maximized when executed on vector architectures. The boundary integral matrix-vector product involves the fully populated matrix, [g]. If uniform surface elements are used in the discretization, this matrix-vector product may be expressed as a truncated,discrete,linear convolution and thus amenable to efficient calculation using the Fast Fourier Transform (FFT). Although uniform zoning imposes restrictions on the geometries which can be analyzed by this FEM-BI technique, the resulting memory and computational efficiencies have proven to be well worth the sacrifice. The boundary integral product is implemented as described in Appendix B with the following problem specific exception. The cross-term arrays do not possess the property: G,.[m - mn n-n'] = Gz [n - m, n - n] and hence the periodic replication rule (B-38) cannot be used here. Instead, a different replication rule was developed in performing the discrete convolution. The resulting code has proven robust even for large arrays. In addition, the cross-term replication rule must be changed for wraparound arrays since the cavity discretization is fundamentally different than is the case for cylindrical-rectangular cavities. We note that these replication rules are not unique and are implementation dependent. The resulting replication rules used in this project are given as Fortran code in Appendix C. 22

5 Scattering and Radiation Once the cavity aperture and volume electric fields have been determined for either an external excitation (scattering) or an internal excitation (radiation), several engineering quantities may be calculated. The aperture fields may be used to determine the Radar Cross Section (RCS) for scattering or the radiated fields for antennas. This entails the convolution of the aperture fields with an appropriate Green's function. In addition, the input impedance may be calculated by using the interior cavity fields. In this section we will present the relevant formula for calculating the far-zone radiated/scattered fields and the input impedance from the electric fields. 5.1 Far-field Evaluation Two of the most important applications of the presented formulation deal with the calculation of the cavity's RCS and the radiation pattern due to sources placed within or on the aperture of the cavity. We begin with the integral representation of the scattered magnetic field in terms of the aperture fields. We have, Hs(r, O) = jYoka G2(r, 0, 0; a, ). p(a, O, z) x E~a a, z) do dz (32) with (r,0,4) indicating the observation point in spherical coordinates. When the observation point is very far from the cylinder, the dyadic Green's function in (32) can be replaced by its far-zone representation G2(r, 0, [; a, X ) G ' + Gz0 + G (33) where the unprimed unit vectors are functions of the observation position and primed ones are functions of the integration point in (32). The components of this far-zone Green's function are determined by a mode matching procedure giving Gj ~ j 2ko0 cos z kC n ejn{({+(_,)) (27r)2 (koa sin) 0) 23

OZ k Co 0' in(- (O_ G 7- C E C. 2C (2w)2 a H'(2) (ka sin 0) 2 k, 1 jn(,+(-')) ('34) (" 2 7r)2a sin 0C z Hn2(kasn 3n)2 a siI1$n= - oo H n(')(k,aa sin O) As one might expect, these series converge rather slowly for large kIa sin 0. They must therefore be recast in another form by employing WYatson's transformation and Fock theory as described before. We have, 41, Cos 0 kcs0'k, i 2 Gei cs ' ( l)Pe-jkosiOP [Y(O)(m.Dp) j k 9 (1) (m(DP)] GG k=sin0 jkscosiz 2 4w, zsijnain0~ p=1 G~~z Co Sn$kj, cos OZ CkasnOD 90 wp 47r = p=1 where the appropriate Fock functions are given by (214). As expected, when e ~, the formulas attributed to Goriainov [15] -k Cos0 ek, osOZ' eJkoasinOl8D [(o)("12~ ) r - 42w sin 0sn are mor uste foul. Th fatr-zoecateed t orino raitd il cnb cmue nuerca by usin cos r (3)(), cos (36) in (32). h/ 4rr (P I k X,a sin 0 k, sin 0 j k,,, Cos 0 Z' jk,,,asinO(Dj (0)(771(DI) * G ~ Sila z Ccs -. j~~i~~~I" 9 ~~l' 4 7 0 72 jko cos Oz' k,,a sin 04D1. l,,,,,' G e e3fO(?(j '-)an sin 0 e,, k,,, a sin cos F mro)(-m cos (36) are more: Useful. The far-zone scattered or radiated field can be computed numerically by using either (34),(35), or (36) in (32). )24

For the scattering problem, the RCS is most often the quantity of interest. This is given by oa(, ) = lim 47r2I( )I (37) Alternatively, the antenna gain may be computed from the far-zone fields as GdB(0,) = 101ogo [4 ( )2E r(0) 2] + 101ogo [ ] (38) + \cm o Zo,, (8) where Am is the wavelength in centimeters, Rin is the input resistance which is given in the next section and E' is the radiated electric field as r -- oo. For comparison with other techniques which define the RCS in terms of the electric field instead of the magnetic field in (37), these two fields may be related by E, = -Zo H Ee = ZoH, (39) in the far-zone. This must be considered when referring to the polarization of an incident or scattered field. 5.2 Input Impedance In addition to the antenna. pattern and gain, designers are concerned with the input impedance of an antenna for feedline matching purposes. The FEM approach allows the calculation of the input impedance of the radiating structure in a rather eloquent manner. The input impedance is comprised of two contributions [19] Z = ZP + ZD (40) where the first term is the probe's self-impedance which is the impedance ( e.g. the probe's impedance in the absence of the patch) and the second term is the contribution of the patch current to the total input impedance. The probe self-impedance accounts for the finite radius of the probe and hence is 25

omitted when a zero-thickness probe is assumed. Ignoring the probe-feed's self impedance, we have [19] J 1J E(p, -) J. (p, 0 z)pdpdod- (41) 0 where the impressed current is given by (28), Vi refers to the volume elements containing the probe-feed, the electric field is the interior field associated with the feed edge and Io is the constant current impressed on the probe. Utilizing (4) and (28) into (41) yields In h (Pa ) f ()[(,) (.-i)] (42) which must be summed over the four radial edges of the element which contains the feed. 6 Validation Having developed the FEM-BI formulation for cavities recessed in an infinite,metallic cylinder and having implemented the technique in a manner which has low computational and storage demand, a essential task is the validation of the written computer code. Unfortunately, although their is much interest in the scattering and radiation characteristics of conformal patch arrays, a survey of the literature indicates a dearth of published data. We are currently awaiting measured data pertaining to a four-patch wraparound antenna mounted on a low observable test structure which is illustrated in figure 5. Until that data is available, we must resort to limited validation which either considers quasi-planar configurations or curve patches excited by a normally incident plane wave. We shall now look at such preliminary validation of the code. 6.1 Scattering The first validation effort for scattering by cavity-backed patch antennas relies on the fact that a small patch on a very large radius cylinder is quasiplanar and approximates rather well an equal sized planar patch. For our 26

,vrap-arourC-, Datch Cgr.e section r*-130i a.,I ~,i* 3' 1 i 01! I F..." 3" 1 f4t I047' elemners lOl /1001 Top-view K {0 69" 1.047, Figure 5: Wraparound test body. 27

test we chose as a reference a planar 1.448" x 1.083" patch residing on a 2.89" x 2.10" x 0.057" cavity filled with a dielectric having or = 4. The equivalent patch on a 10A cylinder is 6.46~ x 1.083 residing on a 12.90~ x 2.10" x 0.057" cavity. Figure 6 compares this work with the planar result computed using a similar FEM-BI code for cavity-backed antennas recessed in a groundplane. Although figure 6 illustrates only monostatic scattering in the 0 = 90~ plane, additional runs for normally incident monostatic scattering and various bistatic situations yield similar agreement. Comparisons may also be made for elongated cavities and 2-D MoM results.,ong narrow cavities have very little axial interaction for principal plane (0 = 90~) excitation and therefore results based on this formulation should compare well with corresponding 2-D data. It is well known, that the RCS of the 3-D scatterer of length L > Ao is related to the corresponding 2-D scattering of the same cross section via the relation '3D 2 ( A 020 (43) Such a comparison is shown in figure 7 for monostatic scattering by a 45~ x 5A x 0.1A cavity for both principal polarizations. Once again the agreement between the two results is excellent, thus providing a partial validation of the new code for highly curved geometries. We remark that similar agreement has been observed for bistatic scattering in the 0 = 90~ plane. The planar approximation eliminates the effects of curvature, which is a primary interest in this work, and the 2-D comparisons done above are only valid for normal incidence. To consider oblique incidence on a highly curved structure, we resort to comparisons with a Body of Revolution (BOR) code for wraparound cavities. Since the BOR code can only model finite structures. we simulate an infinite cylinder by coherently subtracting the farzone fields of the finite structure without a cavity from similar data which includes the cavity. This process is illustrated in figure 8. This procedure is suitable for near normal incidence and was found acceptable for near grazing incidence as well in the case of H-polarization (a = 90~) case shown in figure 9 where the data is taken for bistatic scattering in the b = 0~ plane due to plane wave incidence at Oi = 90~ and 4i = 0. Although the agreements with the results presented above give us confidence that our implementation yields accurate data, agreement with measured data of the geometry shown in figure 5 would provide a more secure 28

20.0 10.0,, -0.0 -10.0 -20.0,/ ' Cylinder: E-pol -30.0 -------- Cylinder: H-pol 0 " Planar: E-pol -40.0 a Planar: H-pol 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Angle (0) [deg] Figure 6: Comparison of a planar patch (1.488" x 1.083") residing on a 2.89" x 2.10" x 0.059" cavity which is filled with cr = 4 dielectric and a quasi-planar patch on a lO,\o cylinder. 29

"0$ Q-. 01 rj.0) 0n 30.0 20.0 10.0 0.0 -10.0 -20.0 -30.0 -40.0 -50.0 -60.0 -70.0 -80.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle (() [deg] Figure 7: Comparison of 2-D MoM results and this work for a 45~ x 5A0 x 0.1 Ao air-filled cavity. 30

Procedure: Wrap-around Cavity End-cap minus End-cap ~IIL approximately I I Figure 8: Far-zone subtraction scheme for the simulation of an infinite cylinder. 31

20.0 t5 r-" C) * -4 Co < -i c~ 0/3 5A e v.", 10.0 0.0 -10.0 -20.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 9: Comparison of this work and processed BOR data for a 3A0 x O.1lo air-filled wraparound cavity excited by a normally incident H-polarized (a = 90~) plane wave. 32

validation. The test body is a twelve inch long metallic cylinder with a radius of six inches. The cylinder is terminated with 53~ ogival end-caps which minimizes the scattering by the terminations (tips). The wraparound cavity is 2" x 0.16" and four identical metallic patches are symmetrically place around the cylinder where each patch is 1.047" x 0.69". 6.2 Radiation As was the case for scattering, the amount of published data suitable for validating our formulation and associated code is quite scarce. There are various approximate cavity models available for coated cylinders and even some measured data [16, 17, 18]. However, there is no cavity-backed patch antenna data yet available. Nonetheless, we may compare our results with planar results for an initial validation. Figure 10 compares the co-polarized antenna pattern for a 4cm x 3cm patch antenna in a 6cm x 5cm x 0.0795cm cavity recessed in a metallic cylinder whose radtius is allowed to vary. The substrate's permittivity is Er = 2.32 and the feed point is at (als = 2cm, zs = 1cm) relative to the center of the patch. The small glitches in the curves associated with the small radius cylinder is due to the far-zone formulas (34) and not due to the FEM formulation. The series (34) rely on delicate cancellations which are difficult to achieve numerically. Nevertheless, the agreement between the planar result and the case where the cylinder of radius 100 cm (~ 10\ at 3.17 GHz) is excellent. Clearly, curvature has no effect on the forward direction (d = 0~) but the curvature effect is pronounced at ~ = 90~. We believe that this broadening of the pattern as the radius decreased may be explained by considering the apparent size of the patch when observed at (O = 0~,0 = 90~). As the radius of the cylinder decreases, the patch's planar projection becomes smaller resulting in a broader pattern. Another concern for an antenna designer is the effect of the cavity termination on the radiation patterns. Figure 11 shows the gain pattern for the same patch placed on cavities whose azimuthal extend varies from 65~ to 360~ (wraparound). Obviously, near the resonant frequency, the effect of the cavity termination edges is only apparent near 0 = 1800 provided the cavity is small. For this test, the patch occupies 45~ of the cavity and hence the 65~ cavity has only 7.5~ (or 6.5 mm) to spare on either side of the patch. Larger cavities are associated with higher gain in the rear direction (q = 180~) possibly due to creeping wave effects. The 33

-0 c0 v 0.0 -10.0 -20.0 -30.0 -40.0 -50.0 -60.0 a = 5 cm --------- a= 10cm ------ a = 20 cm ------ planar o a= 100 cm. ' I..... I.... I..... ~ 0 \ \ - I-.. i - I-...... I...... I..... I..... I.... %N - % 1, \ ^ %% x '' 14 %,e -,%%...... I..... I..... I................. 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle (~) [deg] Figure 10: Co-polarized antenna patterns for a 4cm x 3cm cm patch antenna on a 6cm x 5cm x 0.0795cm cavity fed at (as = 2cm,z = 1cm) relative to the center of the patch. The cavity is filled with cr = 2.32 dielectric and the radius of the cylinder is allowed to vary. 34

0.0 -10.0 -20.0 _________________ 65~ x 5 cm - 90~ x 5 cm 90"":*.c -30.0 ---- 180~ x 5 cm...... 360~ x 5 cm -40.0..... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle (4) [deg] Figure 11: Effect of cavity termination on the radiation pattern generated by the antenna used in Figure 10. 35

backward gain is of interest to low observable antenna designs and we note that planar approximations are of no use for such antenna designs since there is no backward radiation. In addition to radiation patterns, designers are interested in input impedance calculations. Figure 12 compares our result for a quasi-planar patch and validated data for the planar patch antenna considered previously. Also shown is the input impedance calculated as a function of frequency for a highly curved cylinder (a = 5 cm). The decrease in the input resistance is typical of curved patches (see for example [16]) and we note that the resonant frequency has not moved consistent with T' 1ol mode excitation. We are currently pursuing the validation of our code with reference measured data for curved patches. 7 Future Work In this report, we have presented a FEM-BI formulation suitable for cavitybacked antennas which are recessed in a circular,metallic cylinder. Along with the formulation, we presented an implementation strategy which minimizes both computational and storage demand. Key to this goal was the use of asymptotic formulas for the dyadic Green's function and the use of the FFT-based matrix-vector product in the BiCG solver. The resulting computer program has been validated, to the extent possible, for both scattering and radiation problems. In the near future, we want pursue further validation of the written FEMBI code by comparison with measured data for curved cavity-backed antennas. This will be the first such comparison available to the electromagnetics community and hence will form an important resource for future work. Upon establishing full confidence in our implementation, we will then be able to study effects due to patch and cylinder of curvature on antenna parameter and scattering properties. We believe that the numerical laboratory developed herein will provide a powerful tool for the analysis and design of conformal antenna arrays on curved surfaces. Having produced a code suitable for antennas mounted in metallic surfaces, we want to continue this study for coated structures. The FEM-BI method will allow accurate computation of the input impedance, radiation patterns, scattering for large complex antennas mounted on inhomogeneous substrates. Another possible implementation could use Absorbing Bound 36

4 cm x 3 cm Patch Antenna - Planar — B — a= 100cm — &- a=5 cm Figure 12: Input impedance for the patch considered in Figure:10. 37

ary Conditions (ABCs) for mesh termination thus retaining a sparse system. However, the accuracy of FEM-ABC formulations for input impedance calculations have to be established and this will be an important part of this investigation. 38

Appendices 39

A BiCG Algorithm The biconjugate gradient (BiCG) algorithm is one member of a family of iterative solvers which have proven useful in computational electromagnetics [20]. The BiCG unlike the conjugate gradient (CCG) method does not guarantee convergence but does have the advantage of utilizing only one matrix-vector product in its symmetric implementation. Although the convergence characteristics of the BiCG algorithm is erratic (see for example figure A-i), it often converges in fewer iterations than the CG algorithm. This appendix lists the BiCG algorithm appropriate for use with symmetric matrices [21]. Consider the system [A] {x} = {b} (A-1) where [A] is a symmetric matrix, {x} is the unknown data vector and {b} is the excitation data vector. The BiCG pseudocode follows (assuming an initial guess {x}1 which may be {0}: Initialize: {r}1 = {b} - [A] {x}, {p}1 = {b }-[A] {x}1 Iterate: < {r}. { r}n > < [A] {P}n {P)n > {}n+l= {x}x, + an {p}n {}n+1 = {r} n- an [A] {p} c < {- }n+ {r}+ > c < {r} { r}n > {P} = {r}n+l + Cn {P}n (A-2) where the Euclidean norm is given by M < {x}, {b} > = L x[m,]b*[m] (A-3) m=1 40

2.00 ( 1.75 - " ' i rs %* 1.50 L 0.2',-, '.-r,:, BCGM 1.00. ~ ~ 75 I a I: 0.25 - r' I *- o 0.0 L,-so 0 25 50 75 100 125 150 175 200 Number of iterations Figure A-1: Convergence behavior of a typical system solve using the BiCG algorithm. algorithm. 41

In (A-2) the subscripts refer to the iteration and the symbol '*' denotes complex conjugation. Many termination criteria have been used in the past. One of the most popular is < {r}+, {r} > < n+ > < (A-4) < {b}, {b} > where e < 0.01 is a typical acceptable tolerance. As can be expected, the tolerance may need to be tightened or relaxed depending on the problem at hand and the desired accuracy of the result. Note that e should be kept small for antenna impedance computations but can 1)e relaxed for scattering and radiation pattern calculations. 42

B FFT-based Matrix-Vector Product The numerical solution of integral equations (IE) (or boundary integrals) is often performed by converting the IE into a linear system of equations using the method of moments (MoM) procedure. The MoM solution often requires the generation and storage of O(N2) matrix entries and, if a direct solver such as LU decomposition is utilized, O(N3) operations are required, where N is the number of unknowns. If an iterative solver such as the biconjugate gradient (BiCG) method is used, the solution can be found in O(r i. N2) operations where r is the number of right-hand sides and i is the number of iterations per right-hand side. However, for circulant or block circulant matrices, the solution may be reached in O(r. i NlogN) operations. This requires the use of FFTs for computing the matrix-vector products in the BiCG algorithm and consequently the resulting solution scheme is often referred to as the BiCG-FFT method. In this appendix, we will present specific examples of this efficient technique for 2-D and 3-D geometries. We will first look at the relatively straightforward 2-D problem followed by the necessarily more involved 3-D case. B.1 2-D Integral Equations Suppose a flat resistive strip centered at the origin of the y=0 plane is excited by an E-polarized (TMI) plane wave as shown in figure B-1. The appropriate Electric Field Integral Equation (EFIE) may be formed by enforcing the resistive transition condition on the strip giving E~(x) = R(x)J,(x) + 4 J( x)H2 (kolx - x') dx (B-1) where EW(x) is the incident electric field, R(x) is the normalized resistivity as a function of lateral position, J.(x) is the equivalent electric current on the strip induced by the incident field and ko = 2 is the free-space wavenumber. The time dependency ej3t is assumed and suppressed. In (B-l), the primed coordinate denote the source point while the unprimed ones indicate the test point. Once the current has been determined by solving (B-l), the scattered 43

y p x w/2 w/2 Figure B-l: Flat resistive strip geometry. 44

field in the far-zone is given by -Es( p — );Zo [ JZ(x )eikox cosj (B-2) where Zo is the free-space intrinsic impedance. We proceed with the numerical solution of (B-1) by expanding the unknown current in terms of subdomain basis functions as N-1 J,(x) = E J[n]1 n'(x) (B-3) n=0 where lV '(x) = Wl(x) n Ax- - < x < (n + 1)Ax - = 0 else (B-4) and Ax = N. The weight functions W(.x) may assume various forms, one of the simplest being a pulse (W(x) = 1). Substituting (B-3) into (B-1) and performing Galerkin's testing, we obtain w N-I w J VWn(x)E (x)dx = E J[n} 2 j Vn(x)R(x)dx + 2 n =O 2 n=n J J1 ( -W n(x)x' x')H (I-x'-I\dx'dx} 2 T (B-5) which is the discrete form of the integral equation. When pulse basis functions are used, (B-5) becomes (n+l)Ar- _1 f (n+l)Ax- w 2 E(x)dx = J[n '] {6[n - n 2 R(x)dx + 2 n =0 2 4 / 2 H2) (klx - x'l)dx'dx 4 nx- W Ax-B-6) (B-6) 45

where the Kronecker delta function 6[n - n] = 1 n = n 0 7n n' (B-7) has been introduced. Assuming that the resistivity is constant within a segment (e.g. R[n] = R(x) for nAx < x < (n. + 1)Ax) and making the change of variables x = x + ' = x' + ~-, we have I (n+l)Ax N-1 (n+l)A Ej U -,') d = fJ[ln' ] {/?[A].[nl - n'] + t i(n+l ).x j(n +1)x H2 - I)d d 4nl n Ar (B-8) We now observe that the double integral is in convolutional form and since the segments are of uniform length (Ax), we may introduce the discrete function Z, j(n+l)Ax j( +rn\x H - ( g[n -_ '], = k Lx i2)(x H (x-xI)dd- (B-9) 4 Anx JnA x and rewrite (B-8) as N(n+IAx N-1 N-1 ] E()d = Ax J[n']R[n]6[nt-n'] + Z J[n']g[n- n'] n. no n=O (B-10) The first sum in (B-10) is recognized as the product of a diagonal matrix (AxR[n]6[n - n7]) and the unknown data vector (J[n'], '= 0,1,2,..., N - 1) where as the second sum is a truncated,discrete,linear convolution. The discrete form of the IE (B-10) may be written as a matrix equation [Z]{J} = {f} (B-l1) where the excitation vector entries are given by f[n] = / E -- ) dg (B-12) Jn~r \ ^ 46

and the impedance matrix is given by Z[n, n'] = AxR[n]6[n - ] + g[n - n'] (B-13) Alternatively, in preparing to take advantage of the convolution property, we may rewrite (B-10) as [A xR[n]6[n - n'] {J} + [g[n- n']] {J} = {f} (B-14) Obviously, the first matrix-vector product can be trivially computed in N operations y[n] = JP [n]A'XR[n] n = 0, 1,2,..., N - 1 (B-15) Ibut the second matrix-vector product involves a fully populated matrix and would thus normally require O(N2) operations for its execution. However, since g[n - 7'] is a discrete convolution operator, the product may be computed in only (N\logN) operations by invoking the discrete convolution theorem {JP} [g[- -n']] = H[N]'F' {FD {JP} e.FD {I}} (B-16) where * denotes a Hadamard product and the discrete Fourier transform (DFT) pair is given by I -1 FD{I[]} = f[}]etl e 31 {F[q] } =,- E F[q]7mq (B-17) Q=0 The pulse operator H[N] indicates that only N of a possible 2N - 1 values of the discrete convolution are retained. Also, the DFT (and its inverse) operate on periodic sequences of period Al (all sequences with a tilde are periodic in this appendix). In (B-16), the DFTs must be of length lM > 2N - 1 since a minimum of 2N - 1 entries are needed for the complete specification of g[n - n'] as n and n' vary from 0 to N - 1. The unknown iterate (JP) is given by JP[n] = J[]] 7. = 0,1,2,..., N-1 0 7 = N,N + 1, V+ 2,...,M- 1 (B-18) 47

whereas the Green's function sequence is given by [n] - g[n] n = 0,1,2,...,^- 1 = 0 nN=N NV+1,+ 2,..., -1 = g[AM-] n = + 1, + 2,..., AI-1 (B-19) and we have assumed that g[n n'] g[n' - n]. Note that in (B-19), if Al > 2N- 1 it is necessary to pad g[n] with zeros in the middle of the sequence as shown. Also, only the first row(columnn) of the Green's function matrix ( [g[n - n']] ) need be explicitly computed since this is a Toeplitz matrix thus ensuring low O(N) memory demand. In generating the Green's function sequence, typically the first segment (left most) is used as a source while the non-periodic sequence (g[n]) computed by testing at all segments (e.g. let n' = 0, n = 0,1,..., N - 1). Thus all test segments are to the right of the source segment. If this sequence did not have the symmetry property (g[n - ] = g[n' - n]), tlhe interactions with test segments which are to the left of the source would need to be computed explicitly. Alternatively, if the sequence possessed anti-symmetry (g[n - n'] = -g[n -n]), periodic replication is still possible with the following modification to (B-19) [,] = g[], = 0, 1 2,..., N - 1 = 0 n = N,N + 1, AI + 2,..., - -1 = -g[AI -n] n - +1, + 2,..., M-1 (B-20) Combining (B-15) and (B-16), the matrix vector product in the iteration cycle of the BiCG-FFT method is efficiently computed as [Z]{JJ} = Ax{J[n]R[n]} + H[NV] 1 {FD {J} FD {g}}(B-21) This computation requires O(N + (A/) log2(MI)) operations provided radix-2 FFTs are used in (B-21) and should be compared to the standard matrix vector product N-1 [Z] { JP} = JP[n']Z[n, n] n = 0, 1,2,..., N -1 (B-22) n =0 48

*4 -S (d E 0 C) 10000.0 9000.0 8000.0 7000.0 6000.0 5000.0 4000.0 3000.0 2000.0 1000.0 0.0 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 Number of Unknowns (N) Figure B-2: Comparison of operation count for a full matrix-vector product verses a FFT-based matrix-vector product. 49

Figure B-2 illustrates the comparison between (B-22) and discrete convolution (B-21) computed using radix-2 FFTs. Clearly, for N > 30, (B-16) is more efficient than (B-22). Another 2-D geometry which yields systems which may be converted into circulant matrices is the circular strip such as the one shown in figure B-3. The EFIE for a resistive circular strip is given by E4(0) R( )Jz() + J)J ( )H2) (2ka sin( ) ) do'd 2 (B-23) where a is the radius of the arc and a is the subtended angle. Note that if the arc is closed, the a = 3600. A discrete form of (B-23) may be obtained by using pulse basis and uniform zoning (Ao = ). We have J n..'n =0 -ka (+) A,., (n',x O, 22 (,+ H+o ( 2k sin (I- - d4) di] 4 Kkoa J 2 / / (B-24) which is similar to (B-8). The entries of the Green's function matrix are given by [n = -a n/ H2 ) 2ka sin /( }) )d) do'd 4 AO J',xo 2 d \d (B-25) Indeed, (B-24) may be solved in exactly the same manner as (B-8) using FFTs. The only difference between the flat strip and the circular arc is a possible additional symmetry present in the Green's function sequence. If the arc is closed (a = 360~), we find that g[N-n] = gin] n = 1,2,3,..., - 1 (B-26) which indicates that 'V + 1 entries of the sequence need be computed rather than N. Therefore, the solution of (B-24) for a = 360~ requires roughly 50

. / x Figure B-3: Circular arc geometry. 51

half the number of matrix entry evaluations as compared to an equal length flat strip. For a < 3600, a sinilar symmetry exists to a lesser extent as long as ao > 180~. In this case, although the Green's function sequence (B25) is periodic, it is incorrect to assume that (B-24) now involves a circular convolution. The convolution is still a t.runcated,linear,discrete convolution and in practice the additional symmetry of (B-26) is not exploited unless matrix build time is excessive. B.2 3-D Integral Equations As was the case for 2-D geometries, there are certain 3-D geometries which adlmit efficient solution methodologies by making use of the FFT-based implementation of matrix-vector products. Three of these geometries are the flat plate, an impedance insert in a ground plane and an impedance insert in an infinite metallic circular cylinder. These geometries are shown in figure B4. We shall now proceed with the formulation for the later two problems (e.g. the planar and cylindrical inserts) using a general coordinate system. In the following, the general coordinates (u, v) may be considered as (u = x, v = y) for the planar case and (iL = ap,v = z) for the cylindrical case. An integral equation may be obtained by enforcing the standard impedance boundary cond(ition (SIBC)?, x Ih x E(u,v) = -i1Zoh x H (uv) (B-27) where i is an outward directed unit normal and the total magnetic field is given by H(l, v) = Hg(uv, ) + jkl, I G2(u - U', - I), [n x E(u, v')] dudv(B-28) where S denotes the surface of the insert. In (B-28), the geometrical optics fields Hg~(u, )) is comprised of the incident field H(u, v) and the reflected field Hr(u, v). For the planar geometry, the second kind Green's function is given by G2(u - ', - v) = 2G-o(u t v -, ) = 2 + )2 2 ] c- 2(u-u)2+(v-v')2 (B-29) 52

tz:::: O -::::: -41::.m:.: ^ ~ ~.:.~:..~............. — 7L 2b- / x z y YI I I I I I Icz Figure B-4: 3-D geometries: (a) flat impedance insert in a groundplane (or impedance plate in free-space) (b) impedance insert in a circular cylinder (2a = A0). 5:3

where Go is the free-space ldyadic Green's function and I = ~x.+ yy+z+ is the i(lem factor. For a cylinder, the appropriate Green's function is available as an eigenvalue series [8] for small radius cylinders or as a creeping wave series [9] for large radius cylinders. Irrespective of the form of G2, upon making use of (13-28) into (B-27), we have -Z = X x x E(u,,v) jko h x G2(u -, v - v) [ n' x E( t', v ) du dv' (B-30) We note that for planar geometries, an appropriate right-handed system is (U.vt.) whereas for cylindrical problems, the right-handed system is given by (7n,,bv). i;pon expanding the fields (EH90) and the dyadic Green's function (G2) in terms of tangential comp)ollents (~t,b), (B-30) yields the following coupled set of integral equations Ej(u,v) a: Z., HI0(ut ) = - - + jko [E(u u', v)Gu(u -,v - ') - EU(,v )CGlv(u - au, v - v )]du dv E (it, v) b: -ZoH ~(u, ) = - + jko J [E (l, v )Gu(uu, - L ) - Evut',V )GV (u - u1,v- v )du dvB-31) To convert (B-31) into a discrete set of equations, we expand the unknown electric field components in terms of subdomain basis functions N, -3 Nu -2 Eu(u, v) = y y Eu[t, s]WVMs(u,v) t=0 s=0 Nv-2 Nu-3 Ev(t, v) = E E [ts]Wts(u,v) (B-32) t=0 s=0 54

where ANu denotes the number of nodes in the u-direction and NV is the -inumber of nodes in the /,-direction. As shown in figure B-5. t denotes the row number of the edge discretization and s is the column number. For this example, tile unknowns are associated with the free-edges and the associated basis functions may be represented as VHts(u, u) = (- Vb) (Vt - Vb) (Vt - V) (.t - Vb) 0 1 S 7( u,) = ( r - l ) ( Ur -.t ) (i - il) (ur - 'Ll) 0 if local edge = 1 if local edge = 2 el.se if local edge = 3 if local edge = 4 else (B-33) where the local edge numbering is illustrated inl figure B-6. Note that a free edge is one that is not tangential to any metallic walls and that use of (B-33) requires a finite element type assembly. For the discretization shown in figure B-6, there are (A,,- 1)(NA - 1) elements and a total of 2N,Nv - 3( N + Nv) unknowns. Of these, (N, - 1)(N, - 2) unknowns are associated with the l-dlirectel field while (NI - 2)(N1 - 1) unknowns represent the v-directed field. Substituting the field expansions (B-32) into the coupled integral equations (B-31) and employing the Galerkin's testing procedure, we have Fu[ts] = NV-3 Nu-2 E E Eu[t', s']g,[t - t', s - S'] + t'=0 s'=0 N,-2 Nu-3 t E,[t,.s ]gu,[t - t,.s- ] t = 0,1,2,...,N,- 3 t'=0 s'=0 s = 0, 1,,..., Nu-2 N,,-3 Nu-2 F,[ts] = E E l[t', s']g2,v[t - t',s -.']+ t'=0 S'=0

N - V -I3-4- -+ —I t 11 li i I 4.- - I ~ITJ _JFTi _ _ 5i U t = u t I I s = 0 12 3 N 2 (a) Nv- 2 3 2 1 t =0............................. I -... 1.. ............. I I -.............. I- ----- -----— I........................ I -........ -.- -........... I.........................I - -........ —.1................. -1............................................ 1-............. I I............................... -.1-.......................... --..................................-............. - - -........................................................I...................................................................I U S= 0 12 3 N 3 (b) Figure B-5: Rectangular patch discretization: (a) u-directed edges (b) vdirected edges. 56

2 V -, 3 4 Vb --- I I I Il I I I I I I I I U1 Ur U Figure B-6: Local edge numbering scheme. 57

N,,2 Nu-3 1: 1 E, t'j s']g,, [t - t', s - -s']t = 0,1,2I..., Nv - 2 l=0 s =0 (B-34) where gilL [t - 1', I S - Si] - SK,~1 J U', (U', v')IV (u, v) du dv = jk0 [ I1/V,(u', v)lI7" ( a, v) Gul(it - t', v -v')du'dv'du dv = ko J J ~ V(it', v')VIt/>u, v) GVV (u- ut', u - v')du'dv'dudz) gvt- t',,. - S'] LIu [ t - t', s - s2']1 FU~l, s = Z0 H7(,i, v) H9'(u, v)dudv d Ft,.]=-Z0 J / WLu, v) H90(u,zy dcu dv (B-35) and c refers to the test element while e' denotes the source element. We note that each of the double sums in (B-34) is a truncated,discrete,linear convolution and hence amenable to the BiCG-FFT method. To proceed, we define the 2 ---D DFT pair (analogous to (B-17)) All -I A12 -1_j 2 ptq -F2D {f[tjs] = Z: Z f[(,.5]C7J-AfI2(s~~ t=0 s=0 I M-I -I M2 -1 2_____ tq 2A{F[q, p]} M =:D A1M E Z P[q,p]e ~25Pt) (B-36) II2t = 0 S= Using (B-36), the convolutions in (B-34) can then be rewritten as All- 1 A,12 -I Z: Z E[t', s']g[t - t', s - s']: = -DI Ij {F2 D{E[t, s]}eF 2D 1g[t, S]}} tI0 S =0 58

(B-37) The order of the relevant DFTs must be AI, > 2(number of rows) - 1 and A112 > 2(nIumiber of columns) - 1 where the number of rows and columns of the discretization may vary with each convolution (see figure B5). For example, the first convolution in (B-34) is associated with Li-testing and /-source edges and hence the number of rows and columns is (N, - 2) and (N~ - 1), respectively. The field sequences are loaded into a Al1 x M12 array in row/column order of the field discretization anid the remaining entries form a zero pa(d. The Green's function seqlluence must be loaded into a similar array (in the same manner) and periodic replication must be I)erformed to provide the necessary "negative lags". If the sequence has the property, y[t -,.s - '] = g[t' - t,.s -.s], then tills replication process takes the form Y<, = <M o<t<I - o<s<^2-l = [ + 2 -, s] < t < -1 0.s < - 1 ' - - = (g[t, 12 + 2-s] 0 < t < -1 < s < M - 1 = y[I, + 2 -t, 12 + 2 - s] < t< I-1 <l < < 2 - 1 (B-38) For the presented example, gut[t, s] and gyv[t, s] possess this property while the cross-terms, uv[t, s] and g,,u[t, s], do not. If anti-symmetry is present then a more complex replication scheme similar to (B-20) may be used. Otherwise, all possible lags must be computed requiring longer matrix build time since more than the first ti-directed and,-directed edges need be used as sources. Once the periodic arrays are loaded, the required matrix-vector product for the zti-interactions may be performed in 0((A/1logAl1)(M2logA12)) operations rather than the 0(((Nu - 2)(NA - 3))2) operations required for a standard matrix-vector product. The comparison is shown in figure B-7 with Al1 = 2(N, - 3), A2 = ')2(N - 2) and NA = NU = N. Clearly, when the number of nodes per side exceeds 10-15, the FFT-based matrix-vector product is more efficient than a. conventional matrix-vector product. In practice, the FFT-based product is more efficient than a standard product in terms 59

.60E+06 |.55E+06.50E+06 2.5E+06 --- O(((N-2)(N-3))2).45E+06 -------- O( (MlogM1)(M210og2M2) ).40E+06.35E+06.30E+06 E o.25E+06.20E+06.15E+06.1OE+06.50E+05.......... -.00E+00 0.0 5.0 10.0 15.0 20.0 25.0 30.0 Number of Rows/Columns in Grid (N) Figure B-7: Comparison of oleration count for a full matrix-vector product verses a FFT-based matrix-vector product. N is the number of nodes in each direction of the grid, All = 2(N-3) and A'2 = 2(N-2). 60

of wallclock time for N < 10 since in order to exploit the memory savings afforded by uniform zoning of a convolutional kernel without using FFTs, ad(litional overhead is incurred to match the appropriate matrix entry with the correct vector entries. Similar results are obtain-ed for the other convolutions il (B-34). C Replication Rule (Fortran) This appendix contains the Fortran source code which performs the periodic replication of the four components of tile lbouindary integral sub-matrix. It is included to illustrate the special care required for the cross-terms (- z and z- ) as opposed to tlhe application of (B-38) for the like-terms ( - and z - z). c c Augment the arrays C C c c c Provide negative lags for like-terms c Do row = 1,rowFFT Do column = l,colFFT If((row.LE.rowFFT/2).AND.(column.,E.colFFT/2)) Then c First quadrant, do nothing... ElseIf((row.GT.rowFFT/2).AND.(column.LE.colFFT/2)) Then c Second (iuadrant gUU2D(row,column) = gUUt2D(rowFFT+2-row,column) gVV2D(row,column) = gVV2D(rowFFT+2-row,column) ElseIf((row. LE.rowFFT/2).AND.(column.GT.colFFT/2)) Then c Third quadrant gUU2D(row,column) = gUU2D( rowcolFFT+2-column) gVV2D(row,column) = gVV2D(row,colFFT+2-column) Else c Fourth quadrant gUU2D(row,column) = gUJU2D(rowFFT-+ 2-row,colFFT+ 2-column) gVV2D(row,column) = 61

gVV2D( rowF FT+2-row,colF FT+2-column) Endlf EndDo EndDo c c Special treatment for cross-terms c If(.NOT. wrapAround) Then Do row = 1,rowFFT/2 Do coluImni = I,colFFT/2 C c Replicate for IUV C c c Fourth quadrant from first qluadrant gUV2d(rowFFT+ 1-row,colFFT+1-column) = gUTV2d( row,col umn n) c Second qucadrant from first quadrant gUV2d( roFT+ -owFFT +column) = -guvv2d(row,column) c c Replicate for VU c c Fourth quadrant from first quadrant gVU 2d( rowFFT+ 1-row,colFFT+ 1-column) = gVU2d( row,column) c Third quadrant from first quadrant gVU2(( row,colFFT+ 1-column) = -gVU2d(row,column) EndDo EndDo c Do row = l,rowFFT/2 Do column = 1,colFFT/2 c c Now mirror into first and third quadrants for UV c c First quadrant from second quadrant gUV2d(row,column) = (O.ODO,O.ODO) gUV2d(row,column) = -gUV2d(rowFFT- l-row,column) 62

c Third quadrant from fourth quadrant glTV2d(row,colFFT+1-column) = -gUV2d(rowFFT-1-row,colFFT+ 1-column) c c Now mirror into first and second quadrants for VU c c First quadrant from third quadrant gVU2d(row,column) = (O.ODO,O.ODO) gVU2d(row,columni) = -gVU2d(row,colFFT-1-colurInn) c Second quadrant from fourth quladrant g\V IU2d(rowFFT+l -row,col u ) = -gVU 2d( row FFT+ I-row,colFFT- 1-column) End Do EndDo Else c WRAPAROUND CAVITY Do row = 1,rowFFT/2 Do column = 1,colFFT/2 c c Replicate for JUV c c Fourth quadrant from first quadrant gUJV2d(rowFFT+ 1-row,colFFT+ 1-column) = gU V2d( row,column) c Second quiadrant from first quadrant gUV2d(rowFFT+1-row,column) = gUV2d( row,guvColMax+ 1-column) c c Replicate for VU c c Fourth quadrant from first quadrant gVU2d(rowFFT+ 1-row,colFFT+ 1-column) = gVU2d(row,column) c Third quadrant from first quadrant gVU 2d(row,colFFT+ 1-column) = -gVU2d(row,column) EndDo EndDo 63

c Do row = lrowFFT/2 Do columnn = 1,colFFT/2 c c Now mirror into first and third quadrants for UV c c First quadrant from second quadrant gUV2d(row,column) = (O.ODO,O.ODO) gUV2d(row,column) = -gUV2d(rowFFT-1-row,column) c Third clladrant from fourth lquadrant gU V2d(row,colFFT+ 1-column) -gl TV2d ( rowFFT-l -row.colFFT+1 -column) c c c Now mirror into first and second quadrants for VU c c First quadrant from third quadrant gVU2d(row,column) = (O.ODO,O.ODO) g\VU2d(row,column) = gVU2d(row,colFFT-gvuColMax + column) c Second quadrant from fourth quadrant gVU2d(rowFFT+ -row,column) = gVU2d(rowFFT+l-row,colFFT-gvuColMax + column) End Do EndDo EndIf

References [1] J-M Jin and J.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] J.L. Volakis, A. Chatterjee and J. Gong, "A class of hybrid finite element methods for electromagnetics: A review," JEWA, 1993. [3] Z.J. Cendes, "Vector finite elements for electromagnetic field computation"' IEEE Trans Magnetics, Vol. 27, No. 5, pp. 3958-3966, Sept. 1991. [4] J-M.Jin and J.L. Volakis, "A finite element-boundary integral formulation for scattering by three-dimensional cavity-backed apertures," IEEE Trans..Antennlas and Propagat., Vol. 39, No. 1, pp. 97-104, Jan. 1991. [5] J-M Jin and J.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] C-T Tai, Dyadic Green's Functions in Electromagnetic Theory. International Textbook Co., Scranton, 1971. [9] 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. 65

[10] J. Boersma and S.W. Lee, "Surface field due to a magnetic dipole on a cylinder: Asymptotic expansions of exact solution," Univ. Illinois Electromagtnetics Lab., Report 78-17, 1978. [11] T.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-124-1, Nov. 1984. [12] T.S. Bird, "Accurate asymptotic solution for the surface field due to apertures in a conducting cylinder," IEEE Trans. Antennas and Propagat., \ol. 33, No. 10, pp. 1108-1117, Oct. 1985. [13] 0. Eillarsson, R.E. Kleiiima.l, P. Laurin, and P.L.E. Uslenghi, "Studies in raldar cross sections L - Diffraction alnd scattering by regular bodies IV: The circular cylinder," University of.lichigan Technical Report No. 7133-3-T, 1966. [14] N. A. Logan, "General research in diffraction theory," Lockheed Aircraft Corp., Mlissiles and Space Div., vol. 1 and 2, Report LMSD-288088, Dec. 1959. [15] 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 Radiotekhlnica i Elektronica, Vol. 3) [16] J.S. Dahele, R.J. Mitchell, K.M. Luk and K.F. Lee, "Effect of curvature on characteristics of rectangular patch antenna," Electronics Letters, Vol 23, pp. 748-749, 2 July 1987. [17] K-M Luk, K-F Lee and J.S. Dahele, "Analysis of the cylindricalrectangular patch antenna," IEEE Trans. Antennas and Propagat., Vol. 37, No. 2, pp. 143-147, Feb. 1989. [18] K-L Wong and S-Y Ke, "Cylindrical-rectangular microstrip patch antenna for circular polarization," IEEE Trans. Antennas and Propagat., Vol. 41, No. 2, pp. 246-249, Feb. 1993. [19] T.M. Hashaby, S.M. Ali and J.A. Kong, "Input impedance and radiation pattern of cylindrical-rectangular and wraparound microstrip antennas," 66

IEEE Trans. Antennas and Propagat., Vol. 38, No. 5, pp. 722-731, May 1990. [20] T.K. Sarkar, "On the application of the generalized biconjugate gradient method," JEWA, Vol. 1, No. 3, pp. 223-242, 1987. [21] C.F. Smith, A.F. Peterson and R. Mittra, "The biconjugate gradient method for electromagnetic scattering," IEEE Trans. Antennas and Propagat., Vol. 38, No. 6, pp. 938-940, June 1990. 67