A Combined Finite Element-Boundary Element Formulation for Solution of Axially Symmetric Bodies Jeffery D. Collins John L. Votakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122

g,, fg A I L'\ b:2i

Abstract A new method is presented for the computation of electromagnetic scattering from axially symmetric bodies. To allow the simulation of inhomogeneous cross sections, the method combines the finite element and boundary element techniques. Interior to a fictitious surface enclosing the scattering body, the finite element method is used which results in a sparse submatrix, whereas along the enclosure the Stratton-Chu integral equation is enforced. By choosing the fictitious enclosure to be a right circular cylinder, most of the resulting boundary integrals are convolutional and may therefore be evaluated via the FFT with the system is iteratively solved. In view of the sparse matrix associated with the interior fields, this reduces the storage requirement of the entire system to O(N) making the method attractive for large scale computations. This report describes the details of the corresponding formulation and its numerical implementation.

Contents 1 Introduction 4 2 Analysis 6 2.1 Finite Element Formulation........... 8 2.1.1 Analytic CAP Formulation...................... 8 2.1.2 Discretization of the CAP Equations.. 11 2.2 Boundary Integral Formulation...................... 17 2.2.1 Derivation of the Modal Boundary Integral Equation........ 18 2.2.2 Discretization of the Modal Boundary Integral Equation... 23 3 Scattered Field Computation 28 4 Results 33 A Derivation of Modal Incident Field 43 B Maxwell's Equations for Axisymmetric Media 47 C Derivation of Boundary Conditions 48 C.1 Derivation of Axial Boundary Conditions.................. 48 1

C.2 Derivation of PEC Boundary Conditions................... 50 D Evaluation of Finite Element Contour Integral 53 D.1 Contour Integral Evaluation along Conducting Surfaces.......... 53 D.2 Contour Integral Evaluation along the Axis of Symmetry......... 54 D.3 Contour Integral Inter-element Connection Cancellation.......... 54 E Evaluation of the Finite Element Matrix Elements 55 F Boundary Integral Matrix Elements 60 F.1 Elements of PO.............................. 60 F.2 Elements of pt............................ 61 F.3 Elements of Q................................. 61 F.4 Elements of QO'................................. 62 F.5 Elements of Qt............................... 63 F.6 Self-Cell Evaluation....................... 63 2

List of Figures 2.1 General surface of revolution........................... 7 2.2 Cross section of a generating surface enclosed by the fictitious boundary ran....................................... 11 4.1 Mode 0 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 34 4.2 Mode 1 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 35 4.3 Modes 0+1 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 36 4.4 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for axial incidence....... 37 4.5 TM and TE bistatic scattering pattern from a perfectly conducting sphere of radius A for axial incidence......................... 38 4.6 TM and TE bistatic scattering pattern from a perfectly conducting ogive length 1A and maximum radius of 0.088A for axial incidence....... 39 3

Chapter 1 Introduction A restraining factor in the numerical simulation of three-dimensional structures for electromagnetic scattering computations is the storage requirement associated with the chosen method. For sub-wavelength structures traditional methods [1] have been found to work well. However, for structures spanning several wavelengths, the storage requirement limits the use of these methods. For the special case of axially symmetric structures or bodies of revolution (BOR), a reduction of the storage requirement is accomplished by reducing the three-dimensional problem to a set of two-dimensional ones. Several moment method codes have been developed for the solution of these ([2] - [7] and others). However, for large structures the required storage of O(N2), where N denotes the number of unknowns over the BOR cross section, limits their use. To further reduce the storage requirement, hybrid finite element methods ([8]-[12], etc.) may be used, since the storage associated with the finite element method is O(N) in contrast to the O(N2) requirement of moment methods. These methods differ from one another primarily by the application of the radiation condition. The most accurate 4

method enforces the radiation implicitly through an application of the boundary integral equation over a fictitious boundary enclosure [11], and in this case the storage is still O(N2 ), where Nb is the number of unknowns on the boundary. However, through a judicious choice of the enclosing boundary, the storage requirement can be reduced to O(N). This can be achieved by selecting the enclosing boundary to be rectangular or circular [15], [16], making some of the integrals convolutions which can then be evaluated via the FFT when an iterative solution scheme is employed. The proposed method combines the finite element and boundary element methods for the solution of inhomogeneous bodies of revolution. The coupled potential equations [10] are discretized via the usual finite element method, and the resulting system is augmented by a discrete representation (via the boundary element method [13]) of the Stratton-Chu equations [14]. By choosing a right circular cylinder to enclose the scatterer, some of the integrals become convolutions and their discrete counterparts are then evaluated via the FFT in conjunction with an iterative solution procedure as was done in the two-dimensional case [15]. With some care, the storage is reduced to O(N). In this report, we describe the formulation for the proposed finite element - boundary element method as applied to the body of revolution. Some preliminary results are shown to be in reasonable agreement with the method of moments (MOM). 5

Chapter 2 Analysis Consider the body of revolution (BOR) illustrated in fig. 2.1. To employ the proposed finite element - boundary element (FE/BE) method, the BOR is tightly enclosed in a fictitious finite length cylinder, which divides the entire space into two regions, i.e. the one enclosed by the cylinder and the other exterior to it. Since the interior region is generally inhomogeneous, the finite element method is suited for formulating the fields of that region, whereas the boundary element method is applicable for the exterior free space region. A usual approach [3] for treating BORs is to introduce a Fourier series (in the azimuthal coordinate b) representation of the fields, reducing the problem to a set of two-dimensional ones. The finite element - boundary element method is then used to compute each modal field and the final result is found by adding the modal fields. In the following, we present the finite element and boundary element formulations for each mode. First, the finite element formulation is developed.

ImlI (JO C)QC s-I~ CD~ 0'cC

2.1 Finite Element Formulation In this section, we derive the analytical coupled azimuth potential (CAP) equations [17] which are then discretized via the finite element method. 2.1.1 Analytic CAP Formulation Maxwell's equations in a source free region (a eiwt time dependence is assumed and suppressed) are given by V x E(r) = -jwpH (2.1) V x H(r) = jwaE (2.2) V. D(r() = O (2.3) V* B() = 0 (2.4) (2.5) For axially symmetric media, the fields may be represented as Fourier series in the cylindrical coordinate ~6 as 00 E(r) = E Zm(pz)eJiO (2.6) m=-oo,H() = E hm-(p,z)ejm (2.7) m=-oo and when these are substituted into Maxwell's equations (2.1) and (2.2), we obtain X [jmhmz - 1Z(Rhmo)] = jEremp (2.8) [xjmemz - o(Remo)] = -jlrhmp (2.9) R [ hmp - ORhmz] = jEr(Rem,) (2.10)

R [ O emp- 8 emz] = -jAr(Rhmo) (2.11) [jmh,, - (Rhm)] = jremz (2.12) U jmemp - 8 (Rem+)] = ihz (2.13) with R = kop (2.14) Z = koz (2.15) to be referred to as normalized coordinates. Substituting hmz of (2.13) into (2.8) gives emp = jfm [m (Remk) + R/pr jZ(Rhmo)] (2.16) where f = R -2 m2]_ - (2.17) K2 = /r r (2.18) Substituting emp of (2.8) into (2.13) we obtain hmz = jfm [m [m (Rhmqs) + REr &(Remr)] (2.19) Substituting hmp,, of (2.9) into (2.12) yields emz = jfm [m,Z(Remo) - Rir I (Rhmo,)] (2.20) Substituting emz of (2.12) into (2.9) yields hmp = jfm [m ~ (Rhmo)-Rer a (Rems)] (2.21)

Equations (2.16) through (2.21) may be written in compact form as x'(R, Z) = jfm [m$ X Vthee- prRVt/'h] (2.22) x h(R, Z) = ifm [mx x Vth + ErRVte] (2.23).*e(R,Z) = be/R (2.24) *. h(R, Z) = IOh/R (2.25) where Vt=A = a+ Z (2.26) Rewriting (2.10) and (2.11) as RVt (~ x hm) = -jrlIe (2.27) RVt (4 x em) = jPr.h (2.28) and then substituting (2.22) and (2.23) into them, we obtain the CAP equations Vt [fm(ErRVte + mq x Vth] + = 0 (2.29) Vt [fm(sryRVtemh - mO X Vt'e] + = 0. (2.30) This system may be written in operator form as Lb = 0 (2.31) where Vt. [fmrRVt] + R mvt[fmq x Vt] L = (2.32) -mVt[fmb x Vt] Vt - [fmrRVt] + j and - = [Oe Ih]T (2.33) 10

K^ R A Aaa2 Ra2 I I Z-3 Zal Figure 2.2: Cross section of a generating surface enclosed by the fictitious boundary ra. 2.1.2 Discretization of the CAP Equations To discretize (2.31), we first enclose the generating contours of BOR in a fictitious boundary Par, and the axis of symmetry as shown in Fig. 2.2. The contour ra is chosen to be rectangular in shape thus generating a right circular cylinder. The region interior to ra is divided into Ne linear triangular elements and within each element the corresponding 11

weighted residual expression is written IJRNe(R, Z)RQedSe = 0 (2.34) Se where (RNje) is the weighting function and RQ is the residual. Further, Nie is the usual linear shape function found in any finite element book [18]. Using this definition, (2.29) and (2.30) may be written JJRNe {Vt * [fm(ErRVtbe + mq X Vt rh] +'OR dS = 0 (2.35) JJRNie {Vt [fm(prRVth - mO X Vt'e] + R } dS = 0 (2.36) Se and upon using the identity (RNe)Vt AT = Vt (RNTeA) - * Vt(RNe) (2.37) we obtain JJRNf [V. {RNr fm (erRVt'be + mb x Vt'h)} Se +Er,eN -i f (ErRVti + mq x V ) Vtbh) Vt(RNi)] dSe = 0 (2.38) JJRNe [Vt {RNtfm (rLRVtlkh - mb x bVt' ) e Se +lrIhNi - fm (rRVt'h - mO x Vtoe) Vt(RNi)] dSe = O (2.39) Further, by invoking the divergence theorem (2.38) and (2.39) may be written as JJ { [-fm (rRVtoP + m x Vt4h)] * Vt(RNie ) + EroleNt} dSe Se + f [RNm (ErRVte + m X Vtih)] dI = (2.40) JJ { [-fm (PrRVt -Oh - mO x Vt.e) Vt(RNe) + 1?hNi} dSe Se +~ in* [RNfm ( rRVt/)hm X Vt'e)] dI = 0 (2.41) 12

where ni is the outward normal along the boundary, Ce, of the eth element. Finally, these may be simplified by making use of (2.22) and (2.23), yielding If { [-fm (ErRVtoe + mq X Vt'h)]. Vt(RN ) + EreNi } dSe Se - fIc RNi(jhmt)dle = 0 (2.42) JJ { [ —f (rRVth - m x Vt'e)] *Vt(RNe) + LrVhNi} dSe Se +/ RNje(jem)dle = 0 (2.43) where emt = t.Em (2.44) hint = t hm (2.45) with t= n x (2.46) To form a system of equations over the eth element, the fields are represented as a linear basis expansion as 3?e(R, Z) = E )b ejRNj(R, Z) (2.47) j= — 3 tIh(R, Z) = > tkejRNe(R, Z) (2.48) j=1 Substituting these into (2.42) and (2.43) yields j=LI {R [fErVt(RN). V(RNe) + ENrNe] -mfme x V(RN;) V,(RNfl)'hj} dSe] 13

L-. ~VRNe(jh t)dle = ~ (2.49) [-fmr RVt (RN). Vt (RNje) +,Lr NieN]' hj j=l C +mq x Vt(RNje) Vt(RNe )be } dSe] + f RNe(jemt)dle = 0 (2.50) To proceed further, it is necessary that we evaluate the integrals over the surface of the element. Assuming Er and p, are constant within a given element, (2.49) and (2.50) 3 ea" [ aieJ -+ bi ] + i, RNi(jhmt)dle = 0 (2.52) j=1 where as; = |J R [-fmVt(RNi). Vt(RNje) + NNj] dSe (2.53) Se b. = J/mfm5 x Vt(RNe). Vt(RNe)dSe (2.54) Se Summing over all elements to obtain a solution for the entire problem region, our system becomes Ne 3 Na E E - -a xejtj hj Cis.hmts a= 0 (2.55) e=l j=l s=1 Ne 3 Na Z -, [ e +,Aas Lo] + 5 clemts = 0 (2.56) e=l j=l s=1 where c= R Ni-9 Pj-dl8 (2.57) 14

and P1' is the pulse function equal to unity in the sth element. Note that in (2.55) and (2.56) the contour integral contribution canceled out except along the boundary ra as shown in Appendix D. In block matrix form (2.55) and (2.56) may be written emma AC, Afx 0 -Baa -Bai -Caa em,/ 0 ACa Ai 0 -BIa -BII 0 jemt 0 (2.58) Baa BaI Caa A~a AA 0 hma 0 BIa BII 0 ACa AA 0 hmqI 0 jhmt in which we have substituted'te and bh with emO and hm,, respectively, and Ne AC = E'alj (2.59) e=l Ne AU = Z,4aj (2.60) e=l Ne B = E a5 (2.61) Na C=Z EC (2.62) s=1 The subscripts on A'C,, -B and C refer to the various regions of SI and its boundary. The elements of (2.59) - (2.62) are derived in Appendix E and are listed as follows 15

aj = [- ajQ1o- (P a + PjLat)Qjj - 2(QyW + 7j 4)Q2o - 2(/3ye + Pj )-)Q21 -I P Q12 - (4'v77 + 3 /je)Q3o +cala Po + (ofag + Oa?)Pii + (,ea + 3e4)P20 + (I3ie7 + 37r )P21 + /3eP12 +'Yt' yP3O] ( 2.63) and =bj = (2e)2 [( -ai- )Qio + 2(/7 - i7j)Q2o] (2.64) where the Ps and Qs are defined in Appendix E. The elements of C are Cs = Cll C8+1,o = C1 (2.65) where Csl = Cs1 = - a2~ for r.2 Cs = (A ~ Rs )(R + R) F(R + MRsr + R s2) upper sign for ral (2.66) cs = {:[IR9(R~ + R1)-jI(R92 + RsRs + Rq2)] lower sign for ra3 To form a complete system, (2.58) must -be appended by a discrete version of the boundary integral equation to be discussed next. 16

2.2 Boundary Integral Formulation The electric and magnetic fields are represented in the unbounded region by E(r) = r(T) + (Tr) (2.67) H(ri) = H?(r) + (r) (2.68) where tE(f) and H7(r) are the incident fields and the scattered fields are given by the Stratton-Chu equations [14] = fJ) -jwp { [n' x r7(')] g(r;) + [n'E(Ti)] V'g(f,rT) S' where r' and T are the source and observation points, respectively and e-JkO Ir-F'l g(r i;') = if -r'I (2.71) is the free space Green's function. It is convenient for computational purposes to eliminate the presence of the normal field components and after some manipulation we obtain ES(r) = ]& {-j ko [n x r/oHf(r)] g(r, T) + 1 [n * V x iHo(TI)] V'g(rT') + [n' x E(T')] x V'g(T,l 7)} dS' (2.72) ~o~ (~)=fIjk [s' x E(i')] g(T,.) 1 i~o (Tk ) o n Vx ) g( - x (TVg(, ) + ['x oH( x 1( V'g(,T)} dS' (2.73) 17

For r = r', the integrals in (2.72) and (2.73) are singular and by removing these singularity, they may be rewritten in terms of principal integrals as IE(r;= (r) + 4 {-jko [n' x ioH(Y')g] g(Q,f') S' fjio I0 ( )] Vg(,) + [n x E(t) x Vg(r) dS' (2.74) 177oH() = 1o() +; {jko [ x (')] g(s') S' k [in'. V x qE(")] V'g(,')+ [i' x 7H(f')] x VGg(, )} dS' (2.75) where we have also made use of (2.67) and (2.68). These must now be enforced on the boundary so that they can be coupled with the FEM equations. Initially, we will allow S' to be a general surface of revolution and will then specialize it to the case of a right circular cylinder. In the next section, we derive the modal boundary integral equations by expressing the fields and the Green's functions as a Fourier series in the cylindrical coordinate 4. The resulting modal equations are then discretized and the resulting subsystem is augmented to the finite element system previously derived. 2.2.1 Derivation of the Modal Boundary Integral Equation Consider the general surface of revolution indicated in fig. 2.1 whose tangential unit vectors are denoted by q and t. The angle v is that between the i and the z-axis and is negative when t points toward the z-axis. Referring to the figure, we may represent the various unit vectors as n = cos v cos + Yos v sin sinv (2.76) = -x sin q + + cos q (2.77) 18

-t- =:n sin v cos + sinvsin + cos v (2.78) x = tsinvcos + f cosvcos - sin (2.79) y = isinvsinq + + cosvsin +cosq (2.80) z = tcosv- nsinv (2.81) Expressing the primed unit vectors in terms of the unprimed unit vectors results in t =' [sin v' sin v cos(4 - O') + cos v cos v'] +n' [cos v' sin v cos(O - O') - cos v sin vi + ~' [sin v sin(O - 4')] (2.82) n = i' [sin v' cos v cos(O -') - sin v cos v'] +n' [cos v' cos v cos(q - q') + sin v sin v'] + ~' [cos v sin(4 -')] (2.83) = t-' [sin v' sin(4 -')] - ii' [cos v' sin(4 - 4')] +' [cos(4 - 4')] -= p' sin(O - O') + 4' cos(o - 4) (2.84) Taking the 4 component of (2.74) and noting the identities, *. (n' x 1oH0) = - 0oHt sin v' sin(O - 4') - 10oHt cos(4 - 4') (2.85) *V'g = -. vg (2.86) n(V'' x 1oH) = [- -,T(p''toH) + -~,(t/oHt)] (2.87). [(n' x ) x V'g] [z'Et sin(4 - 4') + n'E, cos(4 - 4') + /'EO cos v' sin(o - 4')] - V'g (2.88) 19

we may rewrite (2.74) as 4E(7) = E'(r) J f {jko [Lo[t4 sin v' sin(b -') + +7oIt cos(4 - 4')] g(r, f') + -'1 oH4) +, ((i)oH)] s( Vg(',)?) + [Z'Et sin(4 - 4') + W'E cos( + 4/E[' cos v' sin(4 - )] * V'g} p'do'dr (2.89) Further, by carrying out the derivatives of the Green's functions, we have 2E () = E+(r) J fr {jko [77oH sin v' sin(q - 4') + ro&t cos(O - 4/)] g(, T') +C-.1 [-(p )+ (X7oHt)] sin(q O) dg + (Et sin(4 -') + Eq[p' cos v' cos( - 4/) -p cos v' + (z - z') sin v' cos( - 1')]) - (2.90) Ro dRo pd dJI' (2.90) in which Ro = /p2 + p2 _ 2pp' cos( - 4) + (z - z')2 (2.91) To generate the corresponding integral equations for the modal components, the fields and Green's function may be expanded as 00 E() = m(p,z)ejm4 (2.92) m=-oo 7oH() = Z hm(p,z)ejmO (2.93) m=-00 9 ()(, ) = g()(p, pI, z, z)en(O-') (2.94) n=-oo 20 20

where m (P, z) 2= E(P, u, z)e- )'udu (2.95) hm(P, z) =ro 2 (P u, z)e-imUdu (2.96) 1 /r e-jkoR g()(p, p', z, z') = n(p, P z) = 7o 4 cos(nu)du (2.97) g(9 )(, p',, Z') = [9nI (P, p', Z, Z') + 9n+1 (p,', z, Z')] 1 vr e-jkoR =-f cos U cos(nu)du (2.98) (2)(p, p', z, zt)= -- gn (P p Z) =- [9n-1 (P P' Z Z') - 9n+1 (P,, z, ZIA)] =*i e-jkoR - r sin u ~l sin(nu)du (2.99) g(0) (P, P', Z, Z') = g~n(p P" Z, Z') = 2 - 0 R cos(nu)du (2.100) 7rk Jo dR 9n (p, P, Z, z') = 2 [9gn- (P P Z, Z,') + gn+1 ( P,, Z,')] 1 ~P 1 dg cos u-= --- cos(nu)du (2.101) rk2Jo R dR gn'(P, p, z, z') = -2 [n-l (P, p, z, z') -- n+1 (p, p, z, z')] j'. ldg 7k2 f sin U-= - sin(nu)du (2.102) kJo RRdR with R= 2pp cos u(z z)2 (2.103) Substituting these into (2.90) yields noo Ioo 00 00 i' + ejnof f7i {jko [hmq sin v' (2) + hmtg(l)] n=-oo m=-oo +~k [- j-T(p hmo + jmhmt] ktg(2) 21

.-emt(z - z' )k2(2) + ek(p' cosg(1) -p cos v'g + (z - z') sin v'gl)')} ei(m-n)'lp'd'ldr (2.104) and by multiplying each side by e-jP" and integrating over (0, 27r) to extract the mth modal equation results in -emo (p, z) = e (p, z) +2ir f {jhm, sin v'g 2)-j k(p'hmkg$,2)' + (jhmt) [g) + jmgm2)'] -emt(z - z')kog(2)' + emqsko(p' cos v'g(1) -p cos v'gI + (z - z') sin v'g$(V') kop'dr (2.105) after combining terms and where we have used 2,r 2r m = n j ej(m-n)k'dq' = (2.106) 0 otherwise For the case in which r(= ra) is the generating cross section of a right circular cylinder (indicated in fig. 2.2) the integral in (2.105) may be written as a sum of three integrals, one over each side of ra(= 3=1 r) as em,0(p, z) = em(p, z) +27r {ijhmg$,2)-j; (p'hmo )g(2)' + (jhmt) [g() + jMg(2) -emt(z - z3)kog()' + emko(z - z3)gl) } kop'dp' +271 j {-j t(P2m)gm + (jhmt) [gj) + jmg$2)'] Z3 -emt(- + emko(29gl) - pM M)} kop2dz' +27( m2 (2)1 +2r' J {-jhmg (2) 8 (p/hm )g2)' + (jhmt) [gW)+ jmgm2)] 0 - r -O o )' M -emt(Z - z')kog (2)' - emko(z - z')g)' } kop'dp' (2.107) 22

Introducing the normalized coordinates R = kop R' = kop' Z = koz Z' = koz', (2.108) a = ko 9T (2.107) becomes em,(R, Z) = em(R, Z) +R2 {jhmg2) - j (R'hm)g(2)' + (jhmt) [g+ jmg[ )'] +j(jemt)(Z - Z3)g(2) + emq(Z - Z3)g()'} R'dR' + JZ {-j T, (R2hmo)g$2) + (jhmt) [g$) + jmg() +j(jemt)(Z - Z')g(2) + em(R29(g)' - Rg)} R2dZ' + JR2 {jhmg -j (R'hm)g)' + (jhmt) [9g) + jmg(2)'] +j(jemt)(Z - Z')g(2)' - em,(Z- Z')g(1)' R'dR' (2.109) This equation and its dual are discretized in the next section. 2.2.2 Discretization of the Modal Boundary Integral Equation Consider the fig. 2.2 where the rectangular boundary is divided into Na boundary elements and are equal in length along the ra2. Along the boundary, the fields are expanded into pulse basis functions as Na U(R', Z') = E Uj+i P(R +, - R', Z - Z') (2.110) where U represents any one of the components em,, hm~, emt or hmt and P(R R-, Zj+i - Z') - 2 (2.111) 0 otherwise 23

and Z+ = +l+ Rj+ = + (2.112) j+2 2 2 Substituting the pulse basis expansion into (2.109) and simplifying yields.em~(R, Z) = ei (R, Z) +N { {hm}j+ 2i jg(2)R'dR' + E 21 j=Nal +Na2+1 Rj+ i(Rh )j+ jg)' R'dR'+{jhm [gW + jhmgj 2)] R'dR +{j)mt}j j(z - Z3)g9m)'R'dR' + {em}j+ Rj+ (Z - Z3)9(1)'R'dR' Nal +Na2 Zj ~+ S j {(R2hm)'},i+ Jj - jg2) R2dZ' j=Nal+l +{jhmt}j J Z+ [g$) + jmg 2)'] R2dZ' +{jemt}j J j(Z)- Z)g2)'R2dZ' + {em,}j+ ( - )R2dZ'} Zj+(R ) = e( + Na {hm}i+ J - jg$(2)R'dR' +(a= +a) +1 Rj g Ri+1(Ri, = em(Ri,, ZRJ) +{(R'hms)'}i+2jm R- dR + {fjhmt}j f [g) + jmg)2eR] R'dR 2 24R+ 24

~+ j(Z, - Zs)g(2)'R'dR' + {eo}j+ R1+Rd3 +{jemt} j J (Z1 ) + {m}g + J (ZI -Z3)m Rj+ 3 Nal +Na2 Z + f {(R2hmo)'I}j+ - jg(2) R2dZ' j=Nal +1l Z,+a +{jhmt}j Je [() + jmg'(2)] R2dZ' +{emt}j,J j(Zi - + {em}+ f j (R2g(ml) - Ri+' )R2dZ'} /+{(R hm) 3 2)' jg)'d' + {jhmt JI [l + R dR +{jemt}j j (Z Z) Zi). eRR'dm}.114) for the field points on contours 1al and ra3 and em(R2 Z+) = (, z+ ) fNa 0m ]. jg(R2)R'dR' + Z {{hmj~},+ I -Rg+m SR) 2 f~Rj ig R'aR' tfRj~ IsR'dR' +jemt}j, Rii+ Z)g ) {'R'dR' + {em } j+ L (Z Z3)9W)' R'dR' -} 3 2 NR1 +Na2 Z + {(R2m) }R+ -m (R2, Zi+ + 2hmo}j+i 3- jg(2)R'dR' j=Noal, N2+1 RZ+l } Ri - (2)D [1) A jMg (2)' + Y Jg4'L(2)'dR' + {jhmt}j + [ + mm]R'dR' 25 R, Rj~~~~~~~~~~~~~~~~~~~~~~~~~~~.

+{jmt}j IRj+( (2 RZR - {em} 1Rj+Z Zi)g(1YR'dF1 115) +{iemt}2 t i(Zi+ -Z,)g() R dR - 0em}j+l Rj (Zi+ M Rj Ri 2 2 R ) for the field points on contour ra2. The above set of equations may be written more compactly as P2 f-P]c {emo}a -[Pt] {jemt}a- [4QOC + Q' D] {hm}, [Qt] {jihmt}a = {e m (2.116) where the matrix D arises from the derivative {R'hm.}'j+i = o (R'hm=)I R=R Rj{ hm}j -Rj+i{hm }j+ (2.117) 2 8 ~ += IRj+1 - Rj[ and C is a matrix comprised of'l's along the diagonal and superdiagonal. Also, the subscript 2 represents evaluation at the boundary element midpoints. In a parallel fashion, the dual of corresponding (2.116) may be written 4 [-I- P'] C {hm}a- [P ] {ijhmt} + [4QC + Q"D] {emO}a + [Qt] {je,,t} {hm'} (2.118) MO ay The matrices in (2.116) or (2.118) are 3 x 3 in size, each element of which is a matrix corresponding a particular integration and observation (field point) contours. Each element of the submatrices is in Appendix F. For non-self-cell terms, the integrals are evaluated via open formula numerical integration schemes. The self-cell terms are given in detail in Appendix F. The integrals involving gm are computed via Romberg integration with a specified convergence criterion to ensure accurate evaluation for any mode. 26

Finally, augmenting finite element system with that formed by (2.116) and (2.118), we derive the system Aea A a 0 -Baa -BaI -Caa emoa 0 A}~ AtI 0 -Bia -BIj 0 emqI 0 P 0 p -Q 0 _Qt jemt {em}aL e)2 (2.119) Baa BaI Caa AOa AA 0 hmoa 0 BIa BII 0 Aa AA 0 hmOI 0 Q 0 Qt P 0 -Pt jhmt {hM,}a. which is to be solved iteratively and where Q =Q'C + Q''D (2.120) P =-,[~r- P~]C (2.121) 27

Chapter 3 Scattered Field Computation In the far field the scattered fields are given by E(r) = EO(r)q + r0oH4o(f) (3.1) 7)ioH () = 77OH(f)O- -E(-)O (3.2) We wish to compute the radar cross section given by [19] lim 4sr2 IF~(f)12 I 1'(f)l2 = lim 4irr 2(r)- = lim 47rr22 (3.3) r-.~ -o ()2 roo i 2 -00 1o(fr)12 For TMz polarization we have M E,(r, 6, r ) = 2j ( e~,(r, 6) sin(mr) (3.4) m=1l M i7oH,(r, 69, ) = h)(r, 6) + 2 ~ h,O,(r, 6) cos(mq) (3.5) m=1l and for TE, polarization we have M E(r,, ) = e ) + 2 (r, 6) + 2 e(r, 6) cos(m4) (3.6) m=1 M 70oHa(r,, ) = 2j ~ hs ~(r, 6) sin(mk) (3.7) m=l 28

These combined with a unit amplitude incident field implies that (3.3) becomes 0TM. = lim 4rrr 2 2 E ee,(r,) sin(mo) + h| + 2 h'm(r, 8) cos(mw) I (3.8) L m=1 m=1 aTE. = lim 4rr2 [e + 2 E e8(r, 0) cos(mq) 2+ 2 E hO(r, 0) sin(mq) (3.9) m —-- m —1 We had previously discretized the Stratton-Chu integral equation for field points on the integration contour as given in (2.113). Eliminating the principle value factor for observation not on the integration contour, the corresponding scattered field equation may be written es(R, Z) = (2 {jhmog2) - j (R/hm)gm2)' + (jhmt) [g() + jmg(m2)'] +j(jemt)(z- Z3)g2)' + em,(Z- Z3)g(')l} R'dR' + J {j;T (R2hmo)gm2) + (jhmt) [g) + jmg$l)'] +j(jert)(Z- z')g($2)'+ emO,(R29g)' - RgI)) R2dZ' + j jh mog()-j -*(R'hmO )g )' + (jhmt) [g) +mm] +j(jemt)(Z - Z)g(2)' _ em(Z - Z)g()')} R'dR' (3.10) We wish to evaluate this expression for large kor = vR2 + T 2. For large r JR2+R2-2RR' cos u + (Z- Z)2 kr- ZZ RR' cos(3.11) Thus, we may write (2.97) 1 e-jk0r -zzI oir*' gm(RI R, Z Z) -- j 2kre e O cos(mu)du (3.12) Noting that the integral is related to the Bessel function of the first kind, we may write 29

(3.12) as e -jkor gm,(R, R', Z, Z') e eJ Z'csfm(R, 6) (3.13) g 2kor where we have used jmJm(3) = -j ei3~o cos(mx)dx (3.14) fm(R', 6) = jmJm(R' sin 6) (3.15) and the fact that R=korsin6 Z=korcos6 (3.16) Likewise (2.98) - (2.102) become e-jkor g()(R, R', Z, Z') _ 2kor ek z' Co~s fm(R', 6) (3.17) g j2kor eJZ' Of (R ) (3.19) s()(R, R', Z, Z')= ~2ko ejz cosof m(R'1) (3.18) gm2or _ tj e-j-.,. I g (R, R', Z, Z') _- ej z cos a fm (R', 8) (3.19) r 2kor jko e-Jkor r 2k-Z ejZ' cosg fcm(R',6) (3.20) gm I r 2kor g~)'(R, Rj, Z, Z') _ Jko e-jk~r ejz' Cos 0fsm(R', 8) (3.21) g RZ r 2kor where fm (R', 8) = jm- J (R' sin 8) (3.22) fm(R', 6) = - mjm Jm(R sin ) (3.23) R' sin 6 where the prime on J indicates differentiation with respect to the argument. Substituting these expressions into (3.10) results in the expression e-jkor em(R, Z) 2korfe( 0) (3.24) 30k

where f(m, 9) = ejZ3 cos| {hmqh(jfam) + (jhmt)fcm + (jemt) cos fsm -em,0 cos 0(jifcm)} R'dR' + eZ' cos {(jhmt)fc, + (jemt) cos fam + em, sin 6(jfm)} R2dZ' +eZ cos {-hmr(jfm) + (jhmt)fcm + (jemt) cos Of.m +em, cos O(jfcm)} R'dR' (3.25) Using a midpoint integration to compute the integrals, (3.25) becomes Na fc(m,6) = eiZ3 coS 0 Z {{hm,}.j+' [fam(Rj+L,')] + {jiht}jfcm j=N4,1 +Na2+1 +{jemt}j cos 9fam - {6em,}j+i cos 0(jfcm)} Rj+.Aj Nal +Na2 +? e.+ Co {{jhmt}j fcm(R2, 0) + (jemt}j cos Ofm, j=Nai +1 +{em,}j+ 1 sin (jifm)} R2Aj Nal1 +ezl cosS >E {-{hmcj},+lifsm(Rj+i,9)] + {jhmt}fm + {jet}j cos fm 5=12 +{fem,}j+3 cos(jfm)} Rj+lAj (3.26) Letting fh(m, 0) be the dual of (3.25) we may write (3.8) as TrMA(o,) 2 f[(m, 0) sin(m)2 Ao - m= + fh(0, 9) + 2 E fh(m, 9) cos(m4) 2] (3.27) aTB~(0, 4Ir I fr(O,9)+2 Z fe(m,9)cos(mq) + i. m=1 + 2 1j fh(m, 0)sin(mq) (3.28) m=1 -

where Aj is the length of the jth boundary element. 32

Chapter 4 Results The scattering patterns for three structures are shown in figs. 4.1 - 4.6. In each case the bistatic pattern is computed as a function of Os with Xi = 8. The first body is a conducting right circular cylinder of length 1A and radius o Fig. 4.1 shows both the TE and TM scattering patterns for broadside incidence (92 = 900) for mode 0 and as seen these are in good agreement with the corresponding pattern obtained from the MOM code CICERO [7]. The results for mode 1 are shown in fig. 4.2, while fig. 4.3 shows the sum of modes 0 and 1. Fig. 4.4 depicts the bistatic scattering pattern for the same geometry with axial incidence (0t = 0). Only mode 1 yields a non-zero solution in this case and the results are in good agreement with CICERO. The TE and TM bistatic scattering patterns of a perfectly conducting sphere of radius - are shown in fig. 4.5. Again, these are in good agreement with the data obtained from CICERO. The axial incidence scattering patterns for an ogive of length 1 A and radius.088 A are depicted in fig. 4.6 and agree with those generated by CICERO. 33

p=O.l1, 1=1.0 pc cylinder (m=O) 10.0.....,..'-'/ —' * * *''"'-" 10.0. -10.0 -. -30.0 -50.0. 0.0 30.0 60.0 90.0 120.0 150.0 180.0 e, [deg] (-=0, _=90 __ FE/BE (TE) 0 CICERO (TE) FE/BE (M) E CICERO (TM) Figure 4.1: Mode 0 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 34

p=O.l, 1=1.0X pc cylinder (m=l) 10.0,,.,,,,,, W'.... I......, i 10.0 -10.0 -30.0 -40.0 -50.0.. j -.-. * 0.0 30.0 60.0 90.0 120.0 150.0 180.0 08 [deg] (O —0, Oi=90 )FBE (E) 0 CICERO (TE)......... E/BE (TM) E[ CICERO (TM) Figure 4.2: Mode 1 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 35

p=O.l1, 1=1.0X pc cylinder (m=O-1) 10.0.... 0.0 -10.0' -20.0 -30.0 -40.0 -50.0..... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 o8 [deg] (, —0O, e0=90 ) FWdlE (M) CICERO (TE)......... FE/BE (TM) e CICERO (TM) Figure 4.3: Modes 0+1 TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 36

p=O.l 1, 1=1.OX pc cylinder 0.0.... -10.0 e CIcERo C)....FEBE -30.0 037 0.0 30.0 60.0 90.0 120.0 150.0 180.0 0 CICERO (TE) -----— FE/BE (TM) El CICERO (TM) Figure 4.4: TM and TE bistatic scattering pattern from a perfectly conducting circular cylinder of length 1A and radius 0.lA for axial incidence. 37

Conducting Sphere (p=0.5X) 15.0... 10.0 5.0 - 0.0..'0 -5.0 -10.0 -15.0['''...0.0 30.0 60.0 90.0 120.0 150.0 180.0 0. [deg] ( —-O, — ) — 0E (E) 0 CICERO (TE)......... FE/BE (TM) CICERO (TM) Figure 4.5: TM and TE bistatic scattering pattern from a perfectly conducting sphere of radius A for axial incidence. 38

p=.088X, 1=1.OX pc Ogive 0.0..... -10.0 m' -20.0 - -0 -30.0 - -40.0 - -50.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 0, [deg] ( —=0-0, = FE/BE E) e CICERO CTE)......... E/BE (TM) 3 CICERO C(M) Figure 4.6: TM and TE bistatic scattering pattern from a perfectly conducting ogive length 1A and maximum radius of 0.088A for axial incidence. 39

Bibliography [1] S.M. Rao, D.R. Wilton and A.W. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 409-418, 1982. [2] M.G. Andreasen, "Scattering from bodies of revolution," IEEE Trans. Antennas Propagat., pp. 303-310, March 1965. [3] J.R. Mautz and R.F. Harrington, "Radiation and scattering from bodies of revolution," App. Sci. Res., pp. 405-435, June 1969. [4] T.K. Wu and L.L. Tsai, "Scattering from arbitrarily-shaped dielectric bodies of revolution," Wave Motion, vol. 12, pp. 709-718, Sep.-Oct. 1977. [5] A.W. Glisson and D.R. Wilton, "Simple and efficient numerical techniques for treating bodies of revolution," The University of Mississippi Tech. Rept. No. 105, March 1979. [6] A.W. Glisson, D. Kajfez and J. James, "Evaluation of modes in dielectric resonators using a surface integral equation formulation," IEEE Trans. Microwave Theory Tech., vol. MTT-31, pp. 1023-1029, Dec. 1983. 40

[7] J.M. Putnam and L.N. Medgyesi-Mitschang, "Combined field integral equation formulation for axially inhomogeneous bodies of revolution," McDonnell Douglas Research Laboratories report QAO03, Dec 1987. [8] B.H. McDonald and A. Wexler, "Finite-element solution of unbounded field problems," IEEE Trans. Microwave Theory Tech., vol. MTT-20, pp. 841-847, Dec. 1972. [9] K.K. Mei, "Unimoment method of solving antenna and scattering problems," IEEE Trans. Antennas Prop., vol AP-22, pp. 760-766, Nov. 1974. [10] M.A. Morgan, S.K. Chang and K.K. Mei, "Coupled azimuthal potentials for electromagnetic field problems in inhomogeneous axially symmetric media," IEEE Trans. Antennas Propagat., pp. 413-417, May 1977. [11] Jian-Ming Jin and Valdis V. Liepa, "Application of hybrid finite element method to electromagnetic scattering from coated cylinders," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 50-54, Jan. 1988. [12] S.-Y. Lee, "Highly storage-efficient hybrid finite element/boundary element method for electromagnetic scattering," Electronic Letters, vol. 25, pp. 1273-1274. [13] Brebbia, C.A. and Dominguez, J., Boundary Elements: An Introductory Course, New York: McGraw-Hill, 1989. [14] J.A. Stratton, Electromagnetic Theory, New York: McGraw-Hill, 1941. 41

[15] J.D. Collins, J.M. Jin and J.L. Volakis, "A combined finite element - boundary integral formulation for solution via CGFFT of two-dimensional scattering problems", IEEE Antennas Propagat., AP-38, pp. 1852-1858, Nov. 1990. [16] J.D. Collins, J.M. Jin and J.L. Volakis, "A combined finite element - boundary element formulation for solution of two-dimensional problems via CGFFT," to appear in Electromagnetics, 1991. [17] M.A. Morgan, "Numerical computation of electromagnetic scattering by inhomogeneous dielectric bodies of revolution," Ph.D. Dissertation, June 1976. [18] O.C. Zienkiewicz and K. Morgan, Finite Elements and Approximation, New York: John Wiley and Sons, 1983. [19] C.A. Balanis, Advanced Engineering Electromagnetics, New York: John Wiley and Sons, 1989. 42

Appendix A Derivation of Modal Incident Field Consider a field incident at a point T = (r, 4, z) at an angle (0i, Oi) (as indicated in fig. 2.1 ) of the form (o(,; p, 0, Z) = e-jko f (A.1) 1(9i, i; pi e, z) = -6ie-jko~' (A.2) where the q direction is perpendicular to the plane of incidence and 9i direction is in the plane of incidence. Using = x sin t cos 4 + + sin 8 sin q + + cos 8 (A.3) z' = sin 8O cos Oi + y sin 8' sin 4' + i cos Oi (A.4) the argument of the exponential becomes ko * = kor(-z'. ) (A.5) = -ko [p sin 6' cos(, -') + z cos 8'] (A.6) 43

in cylindrical coordinate system. Using these and the fact that 5i' = -x sin 9i + y cos Oi (A.7) 0i = ~ cos 0i cos 0i + y cos'i sin q - z sin 8 (A.8) x = r sin 0 cos 4 + + cos 8 cos - sin b (A.9) y = r sin e sin 4 + 0 cos e sin + cos (A.10) z = rcosO - Osin (A.11) (A.1) and (A.2) become (i; p, - -, z) = [ sin(9 - qbi) + kcos(4- _i)] ejko[PsinSi cos(&-i')+zcosOi] (A.12) (?;p 4 _ Ot z) = - [ COs t cos(') _ )- cos Ot sin(q- _ ) - sin O'] ejko[psine' cos(O-'i )+zcos8i] (A.13) The previously derived fields may be expanded into a Fourier series in the parameter -_ 5i by first writing (A.1) and (A.2) as 00 Z(6i;p,!- i, Z) = E Zm(6i; p, z)ejm(-i) (A.14) m=-oo oo (ti; p, - ti, z) = IE,m(Oi;p,z)ejm(4-' ) (A.15) m=-oo and then making the definitions f(ei; p, _- 4i) = ejkopsin'i cos(0-~'i) (A.16) f(t; p, q t- = cos(COO - 4)f(0t;, - ) (A.17) fs(OG; p, 4 - qt) = sin(4 - _q)f(O9; p, 4 - _ i) (A.18) 44

Expanding each of these into a Fourier series in ( -'i) and using the fact that f(~q) = f(- q) - fm(u) = f-m(U) (A.19) f(qb) = -f(-q) fm(u) = -f-m(u) (A.20) we have 00oo f(oi;p, 4,- i) = fo(si,p)+ 2 E fm(Oi, p) cos[m( -i)] (A.21) m —=l oo fc(9?; P I - 5i) = fo(9i, p) + 2 E fCm(i', p) cos[m( - q)] (A.22) 00 f8(6i; p, ( 4_ -i) = 2j Z fm(Oi, p) sin[m(4 - ~i)] (A.23) m= —1 where fm(Oip) = 1 eikopsin~' cosu cos(mu)du (A.24) f m(oi p)= os ue,,kopin' cos cos(mu)du (A.25) fm(Oi, P) = -J j sin uejkoPin' cosu sin(mu)du (A.26) Noting the identities jmJm(p) =- e=p~ Cj cos(mx)dx (A.27) jm-1J (,) = - cos ze/p co"s cos(mx)dx (A.28) - jmJ() = - fL sin xejpCOS sin(mx)dx (A.29) where the last two are derived from the first by differentiation with respect to 3 and integration by parts respectively, we may write (A.24)-(A.26) as fm(Oi, p) = jmJm(kop sin i) (A.30) 45

fcm (ei, p) = jm-l Jm(kop sin 9i) (A.31) fm( e',P) = kopsine fm(' p) (A.32) With these, we may proceed to rewrite (A.14) as (Oi; p - i k z)- eJko z cosi [,f m(Gip)+ +fCM(6,p)] ejm(g'i) (A.33) m=-oo Z(9i;p, - i,z)= - _ejkozcos' Z [Pfm(Oi, p) cosO'i - cfm(Oi, p) cosOi -: fm(Oi, p) sin Oi e'm(4-') (A.34) =-o00 or, using (A.19) and (A.20), we have 00 (Oi; p, - qi, z) = eJkozosoi p2j E fmm(Oi,p) sin[m( _qi)] m=1 +f eo(i, p) + 42 j frn (0i, p) cos[m(4 - _i)] (A.35) m=1l Z(oi; p, _ i, Z) = _ejkoz cos' cos O fCo (O', p) + 2 1 fm (', p) cos[m(o - B)I)] - cos 0i 2j Z fm (0i, p) sin[m(4 - )] -z sin i fo(Oi, p) + 2 fm(Oi p) cos[m( - ~i)]], (A.36) In this work, we will use the q components of each of these equations. 46

Appendix B Maxwell's Equations for Axisymmetric Media The usual Maxwell's equations in a source free region are given by V x E() = -jwH (B.1) V x H(r) = jWcE (B.2) V. D(r) = O (B.3) v. B(r;) = O (B.4) In cylindrical coordinates the electric and magnetic fields may be expanded into a Fourier series in S as E(r)= em(p,z)eijm (B.5) m= —oo H1(r)= E hm(p,z)ejmk (B.6) m=-oo Substituting these into Maxwell's equations, we obtain V x em - im x J = -jahm (B.7) P 77 47M

Appendix C Derivation of Boundary Conditions In this appendix, the axial and perfectly conducting boundary conditions are derived. C.1 Derivation of Axial Boundary Conditions Substituting the Fourier series representation of the electric field into the divergence condition we obtain in the normalized cylindrical coordinate system V. ('emejmO) = koeijm [~Eremp + a(Eremp) + jErTemk + 9(Eremz)] = 0 (C.1) Thus, Cremp + R remp)+ Eremz)] = -jmem (C.2) as Morgan had previously derived. Taking the limit of this expression as R -, 0+, we obtain emp + jmem4, = 0 (C.3) Expanding the derivative w.r.t. R in ~ [jmemp- (Remo)] = jiprhmz (C.4) 48~

and taking the limit as R —. 0+, we obtain emp + emqo = 0 (C.5) Combining (C.3) and (C.5) and solving for em, we have (m2 - 1)emlR.-o+ = 0 (C.6) In a similar manner we obtain the dual expression (m2 - 1)hm.IR-0++ = 0 (C.7) For m $ 1, the following axial condition holds emq1iR=O+ = hml1R=+ = o (m 4 1) (C.8) To derive the condition for m = 1, lets first consider emz = Jfm [m j-z(Rem )- Rpr -R(Rhmo)] (C.9) As R -+ 0+, emz - 0 for m $ O. Differentiating (C.9) with respect to Z we have a emz = jfm {mR e R (hmo + R rhmo) +Pr ( 8zhm + R 82 hm0)]} (C.10) Clearly, as R O 0+ e,,mz = 0 for m # 0. Differentiating (C.2) with respect to R after dividing by Er, we obtain A-emp + [ A(fremp) + f (Eremz)] IR=O+ = -jm R-lemtiR=0+ (C.11) Accounting for the behavior of emz and s emz (C.11) becomes 2 -aemp + emp -(ln E,) + jm snemm = 0 (C.12) 49

To find another equation in terms of hemp and -emp, we multiply (C.4) by R and differentiate it with respect to R to obtain R B(Prhmz) + Prrhmz = j [2 +em + R emp -jm hemp] (C.13) Letting R -, O+ we obtain 2j 1- e,nm + m 9- emk = 0 (C.14) Substituting (C.14) and (C.5) into (C.12) we obtain (4- m2) e, + e R(ln r) = 0 R = O, m # 0 (C.15) In an analogous fashion, the dual of (C.15) is given by (4- m2) hm + hm, a+(ln r)= 0 R = O, m # 0 (C.16) For r and jr constant in R at the axis of symmetry and for m = 1, (C.15) and (C.16) reduce to aR emu = 0 (C.17) hm 0 (C.18) C.2 Derivation of PEC Boundary Conditions On a perfect conductor the condition n x E = 0 (c.19) Substituting the Fourier series expansion for the field into this boundary condition yields the following condition on each mode n x m = 0 (C.20) 50e

The second Maxwell's equation for the mth mode is given by (see an appendix) V x hm + 3hm x - jWE Tm (C.21) p Crossing this equation with n' and noting that ni hm = O on the conducting surface, we obtain n x (V x h) = 0 (C.22) Carrying out the curl in cylindrical coordates yields n x {p [pWhmz - hmO]+k[ 8hmp iphmz] +Z [ 8(phm) - Ahmp] } = 0 (C.23) Noting the identities n x/ p = (n *z-) (C.24) nx =t (C.25) n x Z = -q(i. p) (C.26) we find that the middle term of (C.23) implies #&hmp =hmz (C.27) and the first and third terms may be written 2n z) [1 9hmz- 8hmo]- q(n *) [ f (phmo) - hmp] = 0 (C.28) Rearranging terms, we have n * [Vt(phmm)- t hmt] = O (C.29) 51

af;~h = 0 (C.30) and we have used the following Vt = A im + A a (C.31) -n = n * Vt (C.32)'Ph = kophm., (C.33) 52

Appendix D Evaluation of Finite Element Contour Integral D.1 Contour Integral Evaluation along Conducting Surfaces It is shown in the appendix that along perfectly conducting surfaces the conditions,e = O (D.1) O.nWh = O (D.2) must hold. During the assembly of the finite element equations (i.e., when the summation over all elements is performed), those rows and columns of the finite element matrix corresponding to nodes on the conducting boundary are eliminated. As a result, the corresponding contour integral vanishes along a conducting boundary. Imposing the condition (D.2) results in the elimination of the associated contour integral since on the conducting surface n ( x Vte) = t lVte, = 0 (D.3) (check this stuff) 53

D.2 Contour Integral Evaluation along the Axis of Symmetry In the appendix, the axial boundary conditions are derived em, = 0 (D.4) hmq = 0 (D.5).emo, = 0 (D.6) iRhm = O (D.7) Conditions (D.4) and (D.5) results in the elimination of the rows and columns of the assembled finite element matrix associated with nodes on the axis. Alternatively, since R -+ 0 all terms in the contour integral are zero by virtue of the chosen weighting function. (may explore the possiblity of a different weighting function which does not guarantee this) D.3 Contour Integral Inter-element Connection Cancellation Since the argument of the contour integrals are tangential fields at the element boundary, they will be continuous between adjacent elements. As a result, the contour integrations along the element intersection cancel. 54

Appendix E Evaluation of the Finite Element Matrix Elements In the evaluation of a! and bij, integrals of the form Pab JJRaZdRdZ (E.1) Se and Qab =J R22 m2 dRdZ (E.2) se Clearly, Qab exhibits singularities for real K. To evaluate this integral, consider an integral of the form I = JJg(R)ZbdRdZ (E.3) Se To evaluate the integral, first transform it to an integration along the element boundary via Zb+1 F(R, Z)= -g(R) b (E.4) Using Stokes' theorem JJ(Vx)CdS= - di (E.5) 55

and V x l = -X g(R)Zb (E.6) dS = q dRdZ (E.7) Inserting these into (E.5) yields JJg(R)ZbdRdZ= - g(R)Zb+ldR (E.8) Via (E.8), (E.1) and (E.2) become respectively Pab = b1 J RaZb+ldR (E.9) and Qab = +1 R2 m2 dR (E.10) T+-1reR2nV2 -m2d where the contour integration is taken in a counterclockwise fashion. For linear triangular elements, re is represented by a summation of three contours, one for each side of the triangle. The variable Z may be thus expressed as Z(R) = vR + ul (E.11) where U -- RZ+ - Rl (E.12) =1+1 - RI vl = Z.- ulRj (E.13) Then Zb+1 may be expressed as Zb+l= (u + Rv)b+l = b+l(l+ R)b+ p= V b+p y) (E.14) 56

where =m) m!(n-m)! (E.15) Thus, by writing the integral in (E.10) as a sum of an integral along each side of the triangular element, it be rewritten as Qab =+ E b+1E ( p )' IbRl - dR (E.16) l=lp1=O P ~, R2JR 2 - rM Clearly, integrals of the form I(n, m)= JdRR+i R (E.17) JR, R2K2 m2dR for n = 0,1,..., a+ b + 1 must be solved. For n = 0 RI+, R1 m=R I(O, m)= < 1 R)+] (E.18) 2l [n(m- RK) - ln(m + R)] m > 0 For n = 1 it is easily shown that (1, m) = ln(m - Rv) + ln(m + RE)] m > 0 (E.19) Using the definition of the principle branch of the natural logarithm in the equations above guarantees that the singularity is properly handled. For values of n > 1, the recursive formula [17] I(n, m) = I(n, O) + — I(n - 2, m) (E.20) is used. Thus, (E.16) may be written in terms of I as 1 3 b+1 1 b 5 Qab 1 ul)Pvb+lI(a + p, m) (E.21) 57

In a similar fashion, p-0 b+1 b+1 1 P ~ b+1 Ra+p+l RI+, Pa b = b+ 1 E VI Ia + p v (E.22)+ 1=1 p=O RI The shape function is written in expanded form as NA,(R, Z)= 2. + / +/Z + -yieR) (E.23) where e= = (/ e?e7e (E.24) e= ZeRe -Zke (E.25) i = M -Rk (E.26) ~e = k - e (E.27) We had derived in section 2.1 a= JJR [-fm, ErVt(RNe)* Vt (RNj) + ErNieNY dSe (E.28) Se Noting that'Vt(RNi')-p [+2Qe ]ie S + z 2e (E.29) NNj = ( [Lc + Z(PWa + jeac) + R(7ay + yjWi) +RZ(, yje +,By,) + Z2p_3je + R2Y1.e] (E.30) Substituting these into (E.28) and reducing we obtain the desired result e = [-4aWQjo - (,O + ejei )Qj -2(Qyc + Yea )Q20 -2(3/; + j/y)Q2 - /3.P,3Qt 1 - (4v7r;' + 3 l~)Qo 58

+esa PeP1 + (f3W + /3ea)P1j + (7a + ajea )P2o + (At iY + / at )P21 +t Pj P12 + at.j P30] (E.31) e (2Qe)2 In a similar manner, we may write bis = JJmf'm x Vt(RNje) Vt(RN,)dSe (E.32) Se as biM = 2e m + e + dSe (E.33) 2.7 Me f e J 7 2fe -I + 2fe and likewise as b= (2)[(l3eaC - /Xea)Qjo + 2(e,/, -Be7)Q2o] (E.34) t3 (2- S )2 t t - X 59 59

Appendix F Boundary Integral Matrix Elements In this appendix, the elements for the discrete boundary integral system are presented. F.1 Elements of PO [P.]ij = o (F.1) ] = [R2g1)'-R+ g' (Ri+, R2, Z1, Z')] R2dZ' (F.2) [l =3 (Z1 - Z3) J ms' (Ri+, R', Z, Z3)R'dR' (F.3) 137 2+2-~,fRj Rj+p P21J.. = (Zi+, - Zi) g()',(R2, R', Zi+,Z)R'dR' (F.4) [P2 j= J [R2g(ml) - R2g (R2,R2, Zi+,, Z)] R2dZ' (F.5) =[23] ij= (Z+ -Z3~)f+ _ g(,)'(R2, R', Z+, Z3)R'dR' (F.6) [3] = (Z3 -Z1) g()'(Ri+ R. Z3, Z))RdR, (F.7) 60e

[p] = I2 [R2g()' -Ri + (Ri+7R2,Z3,Z')] RdZ' (F.8) [Pli = o (F.9) F.2 Elements of pt [po]j = o (F.10) [2 I= z 9 j(Zl - Z')g(2)' (Ri+, R2, Z1, Z')R2dZ' (F.1) [3]ij = (Z1 -Z)j jg )'(Ri+, R',Z1,z3)RdR (F.12) [ i = (Z,1 - Z3)fl jgl2)' (R2, R', Z1i+ )R'dR' (F.13)'3 JRF 2 [3]i = (Zi+ 2) RZ3),/ 1 m)(2 R', Zi+, Z3)R'dR' (F.15) [B32]ii = ZP i(Z3-Z,)9m2)'(Ri+ R', Z3, Z)R'dR' (F.1) 61 F.3 Elements of QP Pill 0 (F.10) [Qf3],j -/+ P113 (Z _ Z-R ijg()(R+~, R', Z1, Z)R'dR' (F.12) 61

[].]ij = fR/+- - jgm)(R2, R', Zi+, Z1)R'dR' (F.22) [Q22], o (F.23) [Q23] i =f Ri )(R21 R', Zi+1 Z3)R dR (F.24) []31i = J+R jg)(Ri+ R', Z3, Z1)R'dR' (F.25) [Q32] = o (F.26) [Q]33J = j+ jg )(R.i+ R', Z3, Z3)R'dR' (F.27) F.4 Elements of QO' I[o];' = / Rgim)-,)(R,+, R Z1, Z)R'dR' (F.28) [Q0]i, = -j9()'(Ri+, R2, Z1, Z')R'dZ' (F.29) [Q12]2 = - j)'(Ri+ R, Z, Z3)R'dR (F.30) Zj+ 2 [Q21 =JR - jg)(R2, R', Z0, Zi)R'dR (F.31) ii= J )(R2 R Z3)RdR (F.33) [Q30] = j j - jg()'(R2', R', Z3, Zi)R'dR' (F.34) [Q3,2] = j - jgm2)(Ri+, R2, Z3, Z')R'dZ' (F.35) zjf+i [ = f (Ri+ R', Z3, Z)R'dR' (F.36) 62 2

F.5 Elements of Qt [Q]li f RjM [gm + jmg R)'i+ 1 RR', Z, ZI)] R'dR' (F.37) [Qj2].= f' [g() + jmg(2,)'(R+, R2, Z1, Z')] R'dZ' (F.38) [3 = IJR-+I [g(l) + jmg(2)'(Ri+i, R',, Z3)] R'dR' (F.39) [Q2] i =j jR [9() + jmg( ) (R2, R', Z,+ l )] R'dR' (F.40) [Q22],j = fz+ [g') + jmg$()'(R2, R2, Zi+, Z')] R'dZ' (F.41) [Q3] j = Rj+ [g) + jmg(2)'(R2, R', Zi+, Z3)] R'dR' (F.42) R [g) + jmg, (Ri+, R',Z3, Z)] R'dR' (F.43) F.6 Self-2Cell Evaluation [312] - z [)+ j)(Ri+, Z3,Z,>] R'dZ' (F.43) The integrals in the matrix elements [ [Qk]ii and [Qkk]ii contain integrable singularities. They could be integrated numerically without modification as long as the singularity point is avoided, but costs excessive computation time. To avoid the resulting excessive computation time and innacuracies, the integrals are evaluated as in [5]. 2 gm(Ri+2,' Zi+,Z')R'dl' 63

= 2 iJ f [ R cos(mu)R'- +2] dudl' 2n J rl 0 Ro Ro +I(Ri+i li+, 1l, /12) (F.46) where Ri+ 12 1 2 i In I(R,,11,12) =. [K 22 R i2 +) i+ dl' T'92 IZ21(F.47) and where K is the complete elliptical integral of the first kind, I may be either Z or R and' J1 = (Ri+ - R')2 + (Zi+1 - Z')2 (F.48) z2 = ( + R')2 + (Zi+, - Z')2 (F.49) Also, Ro = R2 + R'2 - 2R+ R'cos u+ (Z+ Z)2 (F.50) (F.51) The first and second integrals of (F.46) may be computed using an open interval numerical scheme that also avoids the midpoint. The integral expression for the self-cell of P22 may be rewritten as [Zi [P2e2] = - cosu ( ) 2 cos(mu)du - /1 +r (1 j+io) e-jR~ 2 o 2 2 cos(mu)du R2dZ' 2r 2 fo 2Ro z, 1'r (1 + jR+) e-j) 21u2dZ =...+ [ (i: _ sin2() cos(mu)du R2dZ' (F.52) 64~

where we have used the identity 1 - cos u = 2 sin2(2) (F.53) The solution to (F.52) is [ ii2 -2 J f [ Ro eR cos(mu) sin2(~) Zi+,; L Ro Ro 4 [R2U2 + (z+ - ZI,)2]3 /2 dZ' R2 + 2 I'(R2, Zi+ I Zi+, Zi) (F.54) where ['(Rj 1 11, 12) = j2 f dudl' 11 0 [R2U2 + (I - 1,)2]3/2 R3 {(- l )ln [R + (-11)2+R22] +(12- l)ln [Rr + (12 - 1)2+ R22] -(I- 1l)ln(l-n1 -)- (12 - 1)ln(12 - 1)} (F.55) In the same manner we have [Qa] -j f12 f (1 +jo) sin(mu) sin(u)R' 21 0 RO 2R0 mu2Ri+ 1 r2 2 3/2 duR'dl' 2 [R 2 u2 + (Ri+1 - R')2] mRi+ 1 + 2m +t(Ri+1 R Ri,Ri+R ) (F.56) 67r 2 2

where 11 = Ri 12= Ri+l li+. = Ri+ 1 = R' for a = 1 11 = Zi+l 12 = Zi li+. = Zi+ I' = Z' for a = 2 (F.57) 11 = Ri+i 12= Ri li+, = Ri+ l' = R' for a = 3 Finally, we treat each term in [Qta]ii seperately and obtain [Q=a]ii j2: [1g) + jmg (2) (Rj+, R', Zj+, Z')] R'dR' 1f2 [ [Ro cosucos(mu)R', - i+ dudl' 27rT, jo RO Ro +I(R, 1,11, 12) {+jM f e2 - sin(mu)sinuR' jmR2 2RI + 2r 2 (R.+,li+-,ll 2)j (F.58) where (F.57) is used to determine the expression for each value of a. The self cells involved in the other matrices contain non-singular integrands and may thus be integrated numerically without modification. 66

UNIVERSITY OF MICHIGAN 3IIIII Ill 0 1111111 11 i llllll 0252 792 3 9015 02527 7925