RL866 A UNIFIED APPROACH TO INTEGRAL FORMULATIONS FOR ELECTROMAGNETIC SCATTERING Authors: Jian-Ming Jin and John L. Volakis May 1991 RL-866 = RL-866

A UNIFIED APPROACH TO INTEGRAL FORMULATIONS FOR ELECTROMAGNETIC SCATTERING (Tutorial Paper) Jian-Ming Jin and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109-2122 ABSTRACT A unified approach is presented to formulate various boundary and domain-boundary integral methods for electromagnetic scattering. Beginning with the weighted residual form of the wave equation, a variety of integral methods, including some well-known ones, are obtained by choosing different weighting and trial functions. Formulating different methods under such a common framework is important for instructional purposes and enables us to better understand their relationships. In fact, all known integral methods can be derived from the proposed general representation by choosing appropriate weighting and trial functions. To illustrate the approach, while maintaining simplicity, this paper deals only with the two-dimensional case.

I. INTRODUCTION Many numerical methods have been developed in the past for scattering computations. Interestingly, although the intent with all methods is to solve Maxwell's equations subject to given boundary conditions, their development has followed many different approaches. For example, the volume and/or surface integral equation methods are developed by invoking the equivalence principle [1], [2]. In contrast, the finite element methods have traditionally been based on variational principles which are closely related to the concept of conservation of energy [3], [4]. Clearly, although one should expect that all methods be inter-related, this is not apparent. In this paper, we employ a unified approach to formulate various integral and finite element methods. The approach begins with the weighted residual equations for the fields inside and outside the scatterer(s). By choosing appropriate weighting and trial field expansion functions we then derive a variety of numerical formulations including some of the most well-known. This type of derivation is important for instructional purposes because it places all popular numerical formulations under a common framework. As a result, one can compare the similarities and differences of the methods making possible an evaluation of their effectiveness for a given application. In the following, we illustrate the unified approach by restricting our attention to the two-dimensional case for the only purpose of retaining the simplicity of the presentation. Extensions to the threedimensional case are possible but are not considered. In the next section we consider the scattering by a homogeneous cylinder whose formulation involves only boundary integrals. This is followed by the formulation for inhomogeneous cylinders which, as expected, involves both domain and boundary inte 1

grals. II. BOUNDARY INTEGRAL FORMULATIONS Consider the scattering problem illustrated in Figure 1, where a harmonic electromagnetic wave, having an angular frequency A, impinges upon a homogeneous, infinitely long cylinder immersed in an unbounded homogeneous medium. We shall assume that the cylinder has constitutive constants E2 and A2 and its axis is parallel to the z-axis. The constitutive constants of the surrounding medium will be denoted as El1 and [ul. The electromagnetic field outside the cylinder, to be denoted as 1, is comprised of the incident and scattered fields generated by the induced currents and charges in the volume and surface of the cylinder. It can be expressed as 1 = in + in o (1) in which /inc and ~>s denote the incident and scattered fields, respectively, and Qo denotes the infinite region exterior to the cylinder. For E-polarization ( = Ez and for H-polarization s = H1. Since oinC satisfies the wave equation, the scattered field satisfies the same equation V2 S+ k24> = O in QO (2) subject to the Sommerfeld radiation condition lim V(% + jkl5) = (3) where kl = wv/jT and p = /x2+ y2. Similarly, the field inside the cylinder, 22, satisfies the Helmholtz equation V202 + kb22 = 0 in Q (4) 2

where k2 = WV/'22A and Q denotes the region occupied by the cylinder. The fields inside and outside the cylinder are then coupled by the boundary conditions b1 = q2 = 0, U1-n U2 = q on (5) where u = 1/,rt for E-.polarization and u = 1/Er for H-polarization with r, and pr being the relative permittivity and permeability, respectively. Also, r denotes the boundary of the cylinder, separating Q and O,O and?n is the associated unit normal vector pointing from Q to Qo,. Equations (1)-(5) now completely define the boundary-value problem and we are interested in a solution of q> and q2. Let us consider the scattered field first. By denoting the weighting function as W1, the weighted residual form for (2) can be written as /j W1 (V2 + k 2, )ds = 0 in Qoo (6) O0o where f>3 now becomes the trial function for the scattered field. Invoking Green's second identity lJJ(WV2> - $V2W)ds = f (W' - ^ ) dl (7) (6) can then be written as /j + krds (w1i ~ - OsW1 dl J0o Jr\ n 9n ) d + (W -s ) dl = 0 in QO (8) (roo Op OPl) where rFO denotes a circle of radius approaching infinity. Since '>$ satisfies the Sommerfeld radiation condition (3), the integral over rF vanishes provided we also choose Wi 3

to satisfy the same condition. Thus, we have s(v2W1+ k ds- f w' - nO ) dl= - in Q (9) In a similar manner, we obtain the weighted residual form for (4) as 2J 2(VW2 + k2W2)s + (W2 - 2 a) dl = 0 in Q (10) where W2 is another weighting function and 02 now becomes the trial function for the interior field. Equations (9) and (10), together with the boundary conditions (5), form a complete set of integral equations for Os and 02. It remains to choose the weighting functions W1 and W2 and expansions for the trial functions $8 and <2 to obtain a system of algebraic equations for a numerical solution. It will be seen below that different choices of the weighting and trial functions lead to different formulations, some of which are well-known and commonly used. Note that since QOo and QS are homogeneous it would be advantageous to choose the weighting functions W1 and W2 so that the area integrals in (9) and (10) are eliminated and this is the principle followed in the formulations described below. Method 1: Green's Functions for Weighting and Polynomials for Expansion In this approach, we choose W1 = Gl(p, p') and W2 = G2(p, p') where G1 and G2 are the two-dimensional unbounded space Green's functions satisfying the partial differential equation V2G1,2 + kl,2G, - (p- P') (11) 4

where b(p - p') is the usual Dirac delta function. Substituting (11) into (9) and (10) yields 0 for p E Q [Gl(p qp') -S(p)] -=-d (1-c)qS_(p) for pon F (12) I^(p) for p E Q and <2(P) for p e Q JX [G(p, 'j) P' - 2(P')] dl' = c2(p) for p on F (13) 0 for p C QOO where f- denotes the Cauchy principle value integral with the singularities removed, and c = 1/2 if p is on the smooth portion of F or c = a/27r if p is at a sharp corner having an internal angle of a (see Figure 2). From (12) it can be further shown that 4 pC(p) [Gl ( p) - () G (P' P) ] dl' r0 for p EQ (1 (l-c)>1 (p) for p on F (14) 561 (p) for p e foo upon making use of (1) and J inc(p) for p E Q [ p Qinc(p) P, )t - inc(p)Gl(p P ] d ' c Xcinc(p) for pon (15) 0 for p E Qoo Equations (13) and (14) are recognized as the well-known boundary integral equations. To solve for $3 and 02 we may invoke the boundary conditions (5) in (13) and 5

(14) to obtain the coupled set of integral equations (1- c)q(p) + f [-Gl (P,p')q(p') - (pt)] dl = n q () (16a) cq(p) - f [12(p, p')q(p') - G(p,p'P)] dl'= 0 (16b) valid for p on r. Equation (16) has been widely used in moment method implementations [5], [6] and to illustrate this procedure, we expand b and q as N N = Lj j, q= L'qj (17) j=1 j=1 with (j and qj being the unknown expansion coefficients whereas Lj and L'- are either entire domain or subdomain polynomial expansion functions defined on r. Substituting (17) into (16) and applying Galerkin's technique yields the matrix system A B \ b A ( } {l (18) C D q 0 where the elements of the submatrices [A], [B], [C] and [D] are given by Aij = (1 - c) fLi(p)Lj(p)dl - [Li(p)f Lj(p)OG1 P)dl'] dl (19) U1 B 3 u1 f [\i(P) LjL(p')G(ppp')dl'] dl (20) Cj = c Li(p)Lj(p)dl+ j [Li() Lj(p' dl (21) D, -i-i [ i(P) TL(P )G2(p, p')dldl (22) U2J' L and the elements of {b} are given by b, = j Li(p)0(p)(dl (23) (r 6

If we replace the weighting functions Lj(p) in (19)-(23) by 8(p - pi), (18) becomes the well-known point-matching system [2]. Method 2: Eigenfunctions for Weighting and Polynomials for Expansion In this approach, we choose W1 = E1i and W2 = T21i where 1i, and T2i satisfy the Helmholtz equations V21i, + k2li = 0 for p in QOO (24a) and V2t2i + k222i = 0 for p in Q (24b) respectively. Provided we choose E1i to satisfy the Sommerfeld radiation condition (3), from (9) and (10) we then obtain I Fi - ps1i) dl =0 in Qf (25) nr n on on ) and T ()2i 2 an ) dl = 0 in Q (26) Further, by invoking the boundary conditions (5), (25) and (26) become - A dl l= II i - - -- dl) (27a) U1 i q n an an (2q- n dl O (27b) for p on r, representing two coupled integral equations which can also be solved via the numerical procedure illustrated in Method 1. In particular, upon substitution of the expansion of > and q given in (17), we obtain a system of algebraic equations of the same 7

form as in (18). However, the elements of the submatrices are different and in this case we find Aij = 'tl1ijdl (28) Bi= -- j SiL'dl (29) Ul 3 Cij = I diLjdl (30) Dij = - F2iLndl (31) U2 3 and b i i -nc dl (32) It remains to specify Eij and T2i and this may depend on the geometry of the scatterer; however, satisfactory results can be obtained for smooth convex scatterers having nearly circular cross sections by setting TIi and T2i equal to the cylindrical eigenfunctions T(1) and (4i) respectively, defined later in (37) and (39). The approach with this choice of weighting functions is designated as Method 2. We remark that in comparison with (19)-(22), the integrals in (28)-(31) do not have singular integrands and can, thus, be evaluated numerically without difficulty. In addition, this method has the added advantage of not suffering from the internal resonances which are associated with (16). It has recently been employed for a solution of the eigenvalue problem pertaining to optical waveguides [7]. However, the method is not popular and possible reasons for this are (i) it requires the computation of high order Bessel and Hankel functions, (ii) the procedure for treating the integrand singularities in the evaluation of (19)-(22) has been well-established and (iii) Method 1 has been successfully applied to a variety of 8

geometries. In passing, we note that the eigenfunctions are not the only choice for 4li and T2i and particularly, for cylinders whose cross sections deviate substantially from a circle, they may result in ill-conditioned systems. In that case a different choice may be more appropriate. Observing that (25) and (26) are valid for any j1 and T2i provided they satisfy (24), we may choose E1i to be H(2)(kilp- Pi) an1d T2ito be Jo(k2lp-pil) where HO 2)() denotes the zeroth order Hankel function of the second kind, Jo(@) denotes the zeroth order Bessel function, and pi denotes the location of the ith testing point residing on or inside F. The only concern with such a choice is that pi (i = 1,2,3,...) must be judiciously chosen so that H 2)(kilp - Pil) form a complete set for representing S8 and so do Jo(k2\p - pil) for $2. We may also choose 42i to be H 2)(k2\p - Pil) instead of Jo(k2Ip - pil) but Pi must now reside on or outside r so that T2i is non-singular in 2. Finally, we remark that if we choose -li to be H (2)(kilp - Pi) and I2i to be H(2)(k2lp - Pil) with pi residing on r, then (27) becomes (16) and the resulting system given by (28)-(32) becomes identical to that given by (19)-(23) with point-matching. Method 3: Eigenfunctions for both Weighting and Expansion In this method, we again choose the eigenfunction Eij satisfying (24a) as the weighting function TW1. With this choice, (9) becomes (27a) and can be written as rn (41 On On assuming again that nli satisfies the radiation condition (3). We observe that all quantities associated with the right hand side integral of (33) are known and to facilitate its evaluation we may introduce a fictitious circular boundary Fl, sufficiently large to 9

enclose the entire scatterer (see Figure 3). By employing Green's second identity it can be shown that F _ in dl = I q T- lin d (34) Jr a n a ( n an provided no sources of the incident field exist in the region enclosed by r and FI and thus, (33) can be rewritten as tPl dl = f ( iS- tn, qn~n dl. (35) ir On nn on The evaluation of the right hand side integral can now be substantially simplified if /inc is expanded in terms of basis which are orthogonal to jli. We may choose to expand inc as 00 00 nc = in aiti4) = iaiJm(klp)sinmO; m = int(i/2) (36) i=1 i=1 where sin mq is used for even i and cos m$ is used for odd i. Also, Jm(e) denote the mth order Bessel functions and the expansion constants ai can be found by invoking the orthogonality of the expansion functions. Substituting (36) into (35) and choosing i = (1) = H()(kip)csm (37) with H( 2)() denoting the mth order Hankel functions of the second kind, we obtain / idl =: ai (38) 2 On On.... upon invoking the orthogonality of the sine and cosine functions and making use of the Wronskian for the Bessel functions. This provides a global relation between the field and its normal derivative on the boundary. 10

Equation (38) is sufficient for the solution of q1 or 09q(/On provided either of these vanishes on r as is the case for a perfectly conducting boundary. For the general case, however, a solution for both 1i and 0Q0/1n requires that (38) be supplemented by an additional equation which can be obtained from the interior field formulation. An appropriate eigenfunction representation for the interior field $2 is 00 00 2 = gi() = EgiZJm(k2p)sinOm (39) i=l i=l where again rn = int(i/2) and gi are unknown expansion constants. Strictly speaking, this expansion is only valid for p residing on or inside a circle which is completely enclosed by r. However, Waterman [11] showed that the expansion functions?(4) also form a closed complete set capable of representing the unknown surface field on r. He further showed that the expansion (39) converges in the mean on r and is differentiable there. This was proven through the identity ((W 2 (9W2 2 q52OW2 W -) = ( - ) dl (40) On on 2 J n On valid for W2 being non-singular within the region enclosed by r and F2, where F2 is a circular contour completely enclosed by r as shown in Figure 3, and by choosing W2 to be T() and (4). Thus, (39) can be used to approximately represent the boundary field and its normal derivative on ' and by invoking the boundary conditions (5) and substituting (39) into (38) we obtain [Q]{g} = {a} (41) 11

where Q ij = U — ) 2Il _IF?4) dl (42) Q 2 Ul j an 2j n(42) This matrix system can be solved for the expansion coefficients gi which can then be used in (39) to obtain the boundary field and its normal derivative. Elsewhere, the scattered field can be found by evaluating the pertinent boundary integral (see for example (14)). An alternative way to evaluate the scattered field is to expand it as oo oo s = Efi(l) = EfiH 2)(kip)smq p > pi (43) i=l i=1 where fi are the unknown constant coefficients to be determined and p\ is the radius of the circle Fl. Substituting (43) into the integral identity (Waln ) dl = (w1 n - oinc an dl + (vl O Os ) dl (44) J'nn - n "- )n dn (which can be proved using Green's second identity) and setting W1 = I) we obtain f2 ( 4i A <1 i )dl (45) since the first right hand side integral in (44) vanishes when the expansion (36) is introduced. Further, by invoking the boundary conditions (5) and making use of the expansion (39) in (45), we find {f} = [P]{g} (46) where the elements of [P] are given by X (4) (4)(4 P i U2 I (4)a J 4) r - - _ dl (47) 2 U1 l On -2j 12

When (46) is combined with (41) we then obtain the system {f} = [T]{a}, [T] = [P][Q]-' (48) which can be solved by properly truncating the infinite series in (36), (39) and (43). The formulation described herein is the two-dimensional version of the so-called extended boundary condition method which is also known as the T-matrix method or the null field method. The method was originally developed by expanding the Green's dyadic and employing the concept of analytic continuation [8]. It can be also developed by invoking Huygens' principle [9] or, alternatively, by using Schelkunoff's equivalence principle [10]. The formulation presented here is, however, more parallel to the one described in [11] which employs the concept of flux conservation. Similarly to Method 2, this formulation does not exhibit difficulties associated with internal resonances. Furthermore, for scatterers having nearly circular cross sections the number of unknowns required for a given accuracy can be smaller than that for the previous two methods. However, for elongated cylinders, the field representions (39) and (43) may be inadequate or slowly converging (consequently the system may be ill-conditioned) making the method much less attractive. It was seen, though, that (43) can be avoided if (41) is solved independent of (46) and (14) is then used to compute the scattered field. Nevertheless, the approximate representation of (39) for the boundary field and its normal derivative is essential to the method and its use is the major source of error when simulating scatterers with elongated cross sections. The aforementioned three different methods for the solution of the boundary-value problem defined by (1)-(5) are summarized in Table 1. These methods share a common 13

feature that the same type of weighting or trial field expansion functions are employed for both fields inside and outside the scatterer. This is, of course, not required and in fact we may choose to formulate the interior field using Method 1 and the exterior field in accordance with Method 2 or 3. Such combinations lead to a variety of formulations (not necessarily more advantageous than those just presented). In passing, we should also note that the methods described herein are equally applicable to multilayered, conducting and impedance scatterers provided the boundary conditions are accordingly modified. In particular for perfectly conducting and impedance cylinders the internal field vanishes leading to substantial simplication in formulations and a reduction in unknowns. Table 1: Different Choices for Weighting and Expansion Functions Method Weighting Functions Expansion Functions Eqn. # Method 1 Green's Functions Polynomials (18), (19)-(23) Method 2 Eigenfunctions Polynomials (18), (28)-(32) Method 3 Eigenfunctions Eigenfunctions (42), (47), (48) III. DOMAIN-BOUNDARY INTEGRAL FORMULATIONS The solution methods presented in the previous section are restricted to homogeneous scatterers, and whereas the expressions involving the exterior field <1 are still valid for inhomogeneous scatterers, those associated with the interior field S2 must be modified to account for the inhomrogeneity. When the scatterer is inhomogeneous, the appropriate wave equation for the interior field is V7 (u2Vq2) + ko2V2(2 = 0 in Q (49) 14

where X and u are the same as those defined earlier, ko is the free-space wavenumber, and v = Cr for E-polarization and v = Mtr for H-polarization.. The boundary-value problem is then completely defined by (1)-(3), (5) and (49). To derive integral equations for q2, we can follow a similar procedure to that outlined for the homogeneous scatterers. The integral of the weighted residual of (49) can be written as f/ W2[V (u2V42) + kov24}2]ds = 0 in Q (50) and by invoking the identity V * (W2u2V2) = W2[V (u2Vq2)] + u2Vq * VW2 (51) and the divergence theorem we obtain Jj(u2VW2 V<2 - kov2W22ds - u2W2 - dl = 0 in Q (52) As before, different choices for W2 lead to different integral formulations and some of these are discussed below. Method 1: Green's Function for Weighting The Green's function G2 associated with (U2,V2) is usually not available. However, we may choose G1, defined in (11), for W2 to rewrite (52) as ff [U2(p)VGi(p,p'). V+2(p)- k~2v2(p)G(p. p')2(p)] ds 1- u2(p)Gi (p,p') O2(P)dl = 0 (53) and through the repeated use of the identity (51) and the divergence theorem we find 14 [kV2(P)i(p'p, + u2(p)V2G(p,p') + VU2(p) * VGi(p,p')] q2(p)ds 15

+ u(P) [G(p,P) a - (P) q$2(P GlP) dl = 0 (54) a~n an Further, by making use of (11), (54) becomes J { [k2(p) )- k] U2(p')Gl(p, p') + V'2 (p') V'Gl(p, p')} +2(P')ds' + - 2 (P ) [Gi(p,p) a)' (P) - (P') G1(P )] dl'.,nr an' u2(p'b2(p) for p E P = c2(p)q2(p) for ponr (55) 0 for p E Qoo where we have also interchanged the primed and unprimed coordinates. Combining this with (14), which is still valid for the exterior field, through the boundary conditions (5), we finally obtain J/ { [k2(p) -- ] U2(p')Gll(p p') + V'U2(p') V'Gl(p, p')} 42(p')ds' +-/ [Ul - u2(p')] +(P')0G1(P'P )dl' + UQinC(p) J u2(p)q2(p) for p E [(1 - c)ui + cu2(p)] +(p) for p on r (56) U11(p) for p E Qoo which is the integral equation derived in [12]. We note that this integral equation involves only a single unknown field component within the inhomogeneous scatterer and along the contours defining an abrupt change in material constants. It has been shown, however, to be equivalent to the volume integral equations involving the vector polarization currents given in [13]- [15] which require three unknowns for a simulation of general inhomogeneous scatterers and are also associated with higher singular kernels. A solution of (56) for 16

general scatterers has been given in [16] using pulse basis and point matching and has also been implemented with isoparametric elements [17] to yield a more accurate simulation in the case of high contrast dielectrics. Method 2: Eigenfu.nctions for Weighting The eigenfunctions t2i associated with (u2,v2) are usually not available. However, if we choose W2 to be the eigenfunctions Ti satisfying the Helmoholtz equation V2I; + k2'i = 0 for p in Q (57) (52) becomes (u2V.T Vq2 - kov2 2) ds - u2i -2 dl = 0 in Q (58) In (57), k is an arbitrary constant and could be chosen as the average value of k2 for best result. Following the same procedure as described above under Method 1, (58) can be rewritten as Jj [(k2 - k2 )U2i + Vu2 * V[i] b2dS + fu P2( -2 ) n dl =0 in Q (59) which parallels that in (55). This together with (25) for the external fields and the boundary conditions (5) yield the two coupled equations I 4liq - p aic- linc I( l- =) IF i - - - (60a) UJr O n r n an L (k2 - 2k2 U22i + Vu2 * V7i] 02ds + f (Siq - U2 an ) dl = 0 (60b) 2o O r n 17

These integral equations can now be discretized by expanding the interior field, the boundary field and its normal derivative. Clearly, the integrands of (60) have no singularity and their evaluation can therefore be carried out numerically without difficulty. It should be noted however that in most cases (56) yields a more efficient solution because it does not involve the normal derivative of the boundary field. Similarly to that remarked earlier in Method 2 of Section II, we may also choose other solutions to (24a) and (57) for Tij and T' (other than the eigenfunctions). Particularly, if we choose both TIi and Pi to be H(2 (ki\p - pil) with pi residing on r, then (60a) and (60b) can be combined to yield an integral equation identical to (56). Method 3: Subdomain Basis for Weighting and Expansion If region Q is subdivided into a number of small elements and within each element a set of polynomials is used for the weighting and field expansion functions, we then obtain the well-knowrL finite element method. Specifically, the field in Q and its normal derivative on 1F are expressed as NT q0$2 NB 02=ENgj, q =u2 =E Lqj (61) j=1 j=1 where NT denotes the total number of nodes in 2 (including boundary nodes) and NB denotes the number of boundary nodes. Also, Nq and L are known expansion or shape functions chosen so that qj and qj represent, respectively, the unknown 42 and q at the jth node. Substituting (61) into (52) and choosing W2 to be Nf, we obtain the finite element matrix equation rr Kr 1 Or Hrr 0 q I J (62) LKir Kti 4 OI L 0 0 0 18

where the subscripts Fr and I are used to denote the fields at the nodes on r and those interior of r, respectively. The elements of the matrices [K] and [H] are given by i= JJ (=U2VN9 * VNJ - k 2v2N Nq) ds (63) Hij = r Nf Lq dl (64) *r 3 Clearly, the system (62) provides a relation between the boundary field and its normal derivative in a discrete form and to solve this we must introduce an additional relation between the two quantities at the boundary. This must be derived from the exterior field formulation and can take the form of an integral relation based on one of the methods described earlier. Below we discuss various approaches. A. Finite element/boundary element method If the exterior field is formulated in accordance with Methods 1 and 2 described in Section II, the resulting relation between the exterior boundary field and its normal derivative is the first row of (18). When this is combined with (62), we obtain the system A 0 B r b Krr Kfri -Hrr =O 0 (65) Ki r KII 0 q 0 where the elements of [A], [B] and {b} are given by (19)-(23) or (28)-(32), depending on the method used to formulate the exterior field. We observe that the [K] submatrices are sparse and banded whereas [A] and [B] are fully populated matrices. However, substantial memory reduction can be achieved if we first solve (62) to obtain {q r} = [E]{q} and then solve the equation ([A][E] + [B]) {q} = {b} for {q}. This finite element/boundary element method has been recently applied to a variety of two-dimensional scattering com 19

putations [18]-[21] and has also been found to be particularly attractive for scattering by grooves and slots in thick conducting planes [22]-[24]. B. Finite element/extended boundary condition method This method was proposed by Morgan et al. [25] to formulate the scattering by a body of revolution. The exterior fields are formulated following Method 3 of Section II. In particular, from (38) we have df ()q -- dl = ai (66) 2 U1 )ii On where ai correspond to the expansion coefficients associated with the incident field. Equations (62) and (66) now provide a coupled pair of relations between the field and its normal derivative at the boundary. To solve this coupled set of equations we can first expand q and Or in a form compatible with (61) and this results in the same finite element/boundary element method described above with the exterior fields formulated in accordance with Method 2 of Section II. An alternative is to generate from (62) a set of coupled expansion bases for q and Or to be used in (66) for a solution of the expansion coefficients. Specifically, we first choose a set of known bases for {q}, denoted by q^j, and for each of these we can use (62) to compute the corresponding coupled basis for {fr}, denoted as 4j. Since the system (62) is banded and sparse, the computational demands associated with this repetitive process are quite managable. Thus, we can expand q and O>r using a single set of the expansion coefficients gj as N N q=Egjqj, qr = gjiqj (67) j=l j=l where (bj, qj) are the known numerical bases. When these expansions are substituted 20

into (66) we obtain the system [Q]{g} = {a} (68) which is of the same form as in that (41) with the elements of [Q] given by Qi 2 j- ( j ) (69) It remains to choose ij and this can be done rather arbitrarily as long as a set of qj are capable of representing the true solution of q on F. For example, we may choose the entire domain basis qj = sin 2- (for odd j) and qj = cos y~ (for even j) for certain geometries. Once the expansion coefficients gj are solved from (68), the scattered field can be computed by evaluating the integrals involving the boundary field and its normal derivative given in (67). Alternatively, we can follow the procedure illustrated in Method 3 of Section II by first expanding the scattered field as in (43) and then solving for the expansion coefficients fi from (46), with Pij now given by (4).dl (70) C. Unimoment method The unimoment method was first proposed by Mei [26], [27] and is very similar to the finite element/extended boundary condition method. In fact, the latter is an extension of the former. Based on this formulation the boundary r is chosen to be a circle, say F1, of radius pi. As a result, the scattered field can be represented by (43) and through the same procedure described above, a set of coupled bases similar to those given in (67) can 21

be generated from (62) for the fields inside FI. The boundary conditions (5) are then imposed through their weighted residual forms i W (8 + >inc) dl = j Wqrdl (71a) l Wu n\ dl = Wqdl (71b) Jri \ - -n an J / where W denotes the weighting function. Substituting (43) and (67) into the above and choosing a set of weighting functions for W yields a system of equations for a solution of fi and gi. It is then seen that the unimoment method results in a 2Nx2N matrix whereas the finite element/extended boundary condition method leads to an NxN matrix. Further, this approach introduces additional nodal fields to be computed in the region between the structure's boundary and the fictitious circular boundary Ii. In passing, we note that the finite element/extended boundary condition method and the unimoment method usually result in a smaller size matrix than the finite element/boundary element method. This has often been attributed to the expansion of the scattered field in the form of (43), but is not true. In the formulation of the finite element/extended boundary condition method it is clearly seen that the expansion (43) is not an essential element since a solution can be obtained directly from (68). The smaller size matrix should actually be attributed to the use of coupled bases. A matrix of the same size can also be obtained when the boundary integral formulation (16a) or (27a) is used in conjunction with the coupled bases and this becomes obvious when the coupled basis expansion (67) is substituted into (16a) or (27a). The only concern with the use of the coupled bases is how to choose the known bases qj which should be capable of representing the true solution of q on the boundary to be used to compute the coupled 22

bases qOj. The safest choice is to choose a set of orthogonal pulse basis, i.e., qj = 1 and qi = 0 for i $ j. However, this choice offers no advantage over the standard approach adopted in the finite element/boundary element method since the resulting full matrix will have the same size. The coupled bases will have an advantage only if the number of the coupled bases is smaller than the number of the boundary nodes. Table 2: Different Choices for Weighting Function W2 Method Weighting Function W2 Resulting Formulation Eqn. # Method 1 Green's Function Singular Integral Equation (56) Method 2 Eigenfunctions Non-singular Integral Equation (60) Method 3 Subdomain Basis Finite Element Method (62) The three methods presented in this section are summarized in Table 2. In Method 1, subsectional bases are usually chosen to expand the boundary and interior fields, and similar expansion functions could be used in Method 2. For Method 3, subsectional expansion and weighting functions must be chosen and three variations of this method were presented depending on how the exterior fields are formulated. IV. CONCLUDING REMARKS In this paper, we presented a unified approach to formulate various boundary and domain-boundary integral methods, including the finite element method, for electromagnetic scattering. First, we introduced weighting and trial functions to formulate two weighted residual equations for the fields inside and outside the scatterer. It was shown 23

that different choices for the weighting and trial functions result in different integral equations or methods for the solution of the exterior and interior fields. Among the various methods presented in this paper, those employing the unbounded space Green's functions in conjunction with subsectional field expansions (Method 1 of Sections II and III) have been most frequently used primarily because they are capable of treating arbitrarily shaped geometries. Method 2 of Sections II and III, where eigenfunctions are employed for weighting, was presented to merely show that there are other possible ways to formulate the problem. The extended boundary condition method (Method 3 of Section II) could be efficient for far field computations, but it is not likely to provide accurate result for the induced field particularly for those structures which deviate substantially from a circular cross-section. This is primarily due to the approximate representation of the field and its normal derivative at the boundary in terms of cylindrical eigenfunctions. However, by employing the finite element method to formulate the interior fields, the accuracy of the method can be restored. The resulting formulation would be the finite element/extended boundary condition method (Method 3B of Section III). Nevertheless, to date, among the formulations which incorporate the finite element method those discussed under Method 3A of Section III has been most widely used. REFERENCES [1] R. F. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGraw-Hill, 1961. [2] R. F. Harrington, Field Computation by Moment Methods. New York: Macmillan, 24

1968. [3] 0. C. Zienkiewicz, The Finite Element Method, 3rd ed. New York: McGraw-Hill, 1977. [4] P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers. Cambridge University Press, 1983. [5] K. K. Mei and J. Van Bladel, "Scattering by perfectly conducting rectangular cylinder," IEEE Trans. Antennas Propagat., vol. AP-11, pp. 185-192, March 1963. [6] T. K. Wu and L. L. Tsai, "Scattering by arbitrarily cross sectioned layered lossy dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-25, pp. 518-524, July 1977. [7] J.-G. Ma, "Full modes analysis of rectangular optical waveguide using new boundary element method.," International Journal of Infrared and Millimeter Waves, Vol. 11, pp. 425-433, March 1990. [8] P. C. Waterman, "Matrix formulation of electromagnetic scattering," Proc. IEEE, vol. 53, pp. 805-812, Aug. 1965. [9] P. C. Waterman, "Symmetry, unitarity, and geometry in electromagnetic scattering," Phys. Rev. D, series 3, vol. 3, pp. 825-839, Feb. 1971. [10] P. W. Barber and C. Yeh, "Scattering of electromagnetic waves by arbitrarily shaped dielectric bodies," Appl. Opt., vol. 14, pp. 2864-2872, Dec. 1975. 25

[11] P. C. Waterman, "Matrix methods in potential theory and electromagnetic scattering," J. Appl. Phys., vol. 50, pp. 4550-4566, July 1979. [12] J. M. Jin, V. V. Liepa, and C. T. Tai, "A volume-surface integral equation for electromagnetic scattering by inhomogeneous cylinders," The Journal of Electromagnetic Waves and Applications, vol. 2, pp. 573-588, 1988. [13] J. H. Richmond, "Scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Antennas Propagat., vol. AP —13, pp. 334-341, March 1965. [14] J. H. Richmond., "TE-wave scattering by a dielectric cylinder of arbitrary crosssection shape," IEEE Trans. Antennas Propagat., vol. AP-14, pp. 460-464, July 1966. [15] E. H. Newman, "TM and TE scattering by a dielectric/ferrite cylinder in the presence of a half-plane," IEEE Trans. Antennas Propagat., vol. AP-34, pp. 804-813, June 1986. [16] J. M. Jin and V. V. Liepa, "Simple moment method program for computing scattering from complex cylindrical obstacles," Proc. Inst. Elec. Eng., part H, vol. 136, pp. 321-329, 1989. [17] J. M. Jin, J. L. Volakis, and V. V. Liepa, "A moment method solution of the volumesurface integral equation using isoparametric elements and point-matching," IEEE Trans. Microwave Theory Tech., vol. MTT-37, pp. 1641-1645, Oct. 1989. [18] S. P. Marin, "Computing scattering amplitudes for arbitrary cylinders under incident plane waves," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 1045-1049, 26

Nov. 1982. [19] J. M. Jin, "Finite element - boundary element methods for electromagnetic scattering," Ph.D. Thesis, The University of Michigan, 1989. [20] Z. Gong and A. WV. Glisson, "A hybrid equation approach for the solution of electromagnetic scattering problems involving two-dimensional inhomogeneous dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-38, pp. 60-68, Jan. 1990. [21] K. L. Wu, G. Y. Delisle, D. G. Fang, and M. Lecours, "Coupled finite element and boundary element methods in electromagnetic scattering," in PIER 2: Finite Element and Finite Difference Methods in Electromagnetic Scattering, M. A. Morgan, Ed. Elsevier, 1990, ch. 3. [22] S. K. Jeng, "Aperture admittance matrix by finite element method for scattering from a cavity-backed aperture," 1988 IEEE AP-S International Symposium Digest, vol. 3, pp. 1134-1137, June 1988. [23] J. M. Jin and J. L. Volakis, "TE scattering by an inhomogeneously filled aperture in a thick conducting plane," IEEE Trans. Antennas Propagat., vol. AP-38, pp. 1280-1286, Aug. 1990. [24] J. M. Jin and J. L. Volakis, "TM scattering by an inhomogeneously filled aperture in a thick conducting plane," Proc. Inst. Elec. Eng., part H, vol. 137, pp. 153-159, June 1990. [25] M. A. Morgan, C. H. Chen, S. C. Hill, and P. W. Barber, "Finite element-boundary 27

integral formulation for electromagnetic scattering," Wave Motion, vol. 6, pp. 91 -103, Jan. 1984. [26] K. K. Mei, "Unirnoment method of solving antenna and scattering problems," IEEE Trans. Antennas Propagat., vol. AP-22, pp. 760-766, Nov. 1974. [27] S. K. Chang and K. K. Mei, "Application of the unimoment method to electromagnetic scattering of dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-24, pp. 35-42, Jan. 1976. 28

FIGURE CAPTIONS Fig. 1 Geometry of the scattering problem. Fig. 2 Illustration of the internal angle. Fig. 3 Illustration of IF and F2. 29

.4.4 I I QOO I I I4 I4 F I I I I I I I I. I. I..4.4.4 4..4 Fig. 1 17" d W -.4a".4t 00.4 op % QOO F if I I, If I" )I I I I I.4 Fig.- 2 Fig. 3 0? 2,,