A Finite Element - Boundary Integral Method for Electromagnetic Scattering by Jeffery David Collins A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1992 Doctoral Committee: Associate Professor John L. Volakis, Chairperson Associate Professor Linda B. Katehi Professor Richard A. Scott Assistant Research Scientist Jian-Ming Jin Research Scientist Valdis V. Liepa

i~;"'o]I 6c- oI P'

ABSTRACT A Finite Element - Boundary Integral Method for Electromagnetic Scattering by Jeffery David Collins Chairperson: J.L. Volakis A method that combines the finite element and boundary integral techniques for the numerical solution of electromagnetic scattering problems is presented. The finite element method is well known for requiring a low order storage and for its capability to model inhomogeneous structures. Of particular emphasis in this work is the reduction of the storage requirement by terminating the finite element mesh on a boundary in a fashion which renders the boundary integrals in convolutional form. The fast Fourier transform is then used to evaluate these integrals in a conjugate gradient solver, without a need to generate the actual matrix. This method has a marked advantage over traditional integral equation approaches with respect to the storage requirement of highly inhomogeneous structures. Rectangular, circular and ogival mesh termination boundaries are examined for two-dimensional scattering. In the case of axially symmetric structures, the boundary integral matrix storage is reduced by exploiting matrix symmetries and solving

the resulting system via the conjugate gradient method. In each case, several results are presented for various scatterers aimed at validating the method and providing an assessment of its capabilities. Important in methods incorporating boundary integral equations is the issue of internal resonance. A method is implemented for their removal, and is shown to be effective in the two-dimensional and three-dimensional applications.

To Sheila 11

ACKN OWLED GEMENTS I would like to express my appreciation to everyone who helped me in completing this dissertation. First, I would like to thank Dr. Timothy J. Peters for inspiring this work. I would also like to thank Drs. Jian-Ming Jin and Kamal Sarabandi for their insightful technical discussions, invaluable to the completion of this research. Many other members of the Radiation Laboratory, both present and past, have been a source of encouragement and friendship throughout the course of my research. I wish to extend a special thanks to Drs. Kasra Barkeshli, Kyle McDonald, Leland Pierce, Mark Ricoy, Hasnain Syed, Michael Whitt, Norman Vandenburg and Emilie VanDeventer and Messrs. Richard Austin, George Eleftheriades, Andrew Engel Jr., Leo Kempel, Mark Portelli and James Stiles. I would also like to thank Mr. Daniel Ross and Ms. Pamela Haddad for their constant source of inspiration; and Mr. Douglas Craig for his support, encouragement and friendship throughout my entire academic career. I would like to extend a special thanks to my parents Clarence and Sharon, my grandmother Isabelle and my siblings Claudia, Angela, Rhonda and Ray for having the faith in me to complete this work. I wish to also thank to Dr. William Kent for his financial support and encouragement through the final stages of my research; to my committee members for their careful review of this document and helpful comments; and to Prof. John Volakis

for providing the opportunity, financial support, and commitment necessary for the successful completion of this highest degree.

TABLE OF CONTENTS DEDICATION.................................. ii ACKNOWLEDGEMENTS.......................... iii LIST OF FIGURES............................. viii LIST OF TABLES................................ xii CHAPTER I. Introduction.............................. 1 II. Fundamental Concepts........................ 6 2.1 Basic Electromagnetic Theory................. 6 2.2 The FE-BI Approach...................... 8 2.3 The Conjugate Gradient Method............ 13 III. A Finite Element - Boundary Integral Formulation for Twodimensional Scattering with a Rectangular Termination Boundary.................................... 15 3.1 Introduction........................... 15 3.2 Analysis.............................. 15 3.2.1 Discretization of the Object and Field Quantities.. 18 3.2.2 Derivation of the Finite Element Matrix....... 19 3.2.3 Evaluation of Boundary Integral........... 23 3.2.4 A CGFFT Implementation.............. 32 3.3 Computational Considerations................. 34 3.3.1 Storage Efficiency................... 34 3.3.2 Computational Efficiency............... 40 3.4 Far Field Computation.................... 40 3.5 Code Validation......................... 43 3.6 Summary........................... 54

IV. A Finite Element - Boundary Integral Method for Twodimensional Scattering Using Circular and Ogival Termination Boundaries............................ 55 4.1 Introduction........................... 55 4.2 Analysis........................ 56 4.2.1 Circular Enclosure................... 57 4.2.2 Ogival Enclosure................. 66 4.3 A CGFFT Implementation................... 73 4.4 Scattered Field Computation.................. 74 4.4.1 Circular Boundary................... 75 4.4.2 Ogival Boundary................... 76 4.5 Results.............................. 79 4.6 Summary............................. 86 V. The Elimination of Boundary Integral Interior Resonances in Two-dimensional Finite Element - Boundary Integral Formulations................................ 87 5.1 Introduction........................... 87 5.2 Method........................ 88 5.3 Results................................ 91 5.4 Summary............................. 96 VI. A Finite Element - Boundary Integral Formulation for Axially Symmetric Structures................... 97 6.1 Introduction........................... 97 6.2 Finite Element Formulation.. 97 6.2.1 Analytic CAP Formulation.............. 98 6.2.2 Discretization of the CAP Equations........ 101 6.2.3 An Improvement to the Finite Element Formulation 110 6.3 Boundary Integral Formulation................. 114 6.3.1 Derivation of the Modal Boundary Integral Equation 115 6.3.2 Discretization of the Modal Boundary Integral Equation........................... 119 6.4 Sources of Electromagnetic Radiation............. 124 6.4.1 A Plane Wave Excitation............... 125 6.4.2 An Electric Dipole Source............. 129 6.5 Scattered Field Computation.............. 130 6.6 Code Validation and Performance.............. 136 6.6.1 Storage Efficiency................... 136 6.6.2 Computational Efficiency and Accuracy....... 139 6.6.3 Elimination of the Internal Resonance Corruption of the Boundary Integral................ 142

6.6.4 Scattering from Various Test Bodies......... 142 6.7 Conclusions............................ 149 VII. Conclusion............................... 150 7.1 Summary............................. 150 7.2 Future Work......................... 152 BIBLIOGRAPHY............................... 154

LIST OF FIGURES Figure 3.1 Geometry of the scatterer.................... 16 3.2 Partially discretized body.................... 16 3.3 Example of the mesh of a penetrable structure............ 38 3.4 Example of the mesh of an impenetrable structure.......... 38 3.5 Ez backscatter from a.25 x 2A body................. 45 3.6 Hz backscatter from a.25 x 2A body............... 46 3.7 Ez backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties er = 5. - j.5, 1r =... 47 3.8 Hz backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5, p/ = I... 48 3.9 Ez backscatter from a.25 x 1I perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5, tp = 1.5 - j.5 49 3.10 Hz backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5, i, = 1.5 - j.5 50 3.11 Ez backscatter from a.25 x 1A perfect conductor with two A/20 thick top material coatings. The lower layer has the properties Er = 2. - j.5, pr = 1.5 - j.2, and the upper layer has properties r = 4.- j.5, pr = 1.5 - j.2....................... 51 3.12 Hz backscatter from a.25 x 1A perfect conductor with two A/20 thick top material coatings. The lower layer has the properties er = 2. - j.5, pr = 1.5 - j.2, and the upper layer has properties Er = 4. - j.5, r = 1.5 -j.2...................... 52 viii

3.13 Hz scattering from a composite body. Both the perfect conductor and the dielectric body are A/2 in each dimension. The material properties are er = 5. - j.5, lr = 1.5 - j.5............... 53 4.1 Geometry of the scatterer.................... 56 4.2 Partially discretized body in a circular enclosure......... 58 4.3 Partially discretized body with an ogival enclosure.......... 66 4.4 Ez and Hz bistatic echowidth of a perfectly conducting circular cylinder of radius 0.5.......................... 80 4.5 E, and Hz bistatic echowidth of a perfectly conducting circular cylinder with a conductor radius of.5A and a coating thickness of.05A containing material properties er = 5 - j5, ir = 1.5 - jO.5..... 81 4.6 Ez and Hz bistatic echowidth of a coated circular cylinder with a conductor radius of 3A and coating thickness of 0.05A with material properties c, = 5- j5, fL, = 1.5- jO.5................ 82 4.7 Ez and Hz backscatter echowidth of a 0.5 x 1A perfectly conducting ogive................................... 83 4.8 Ez and Hz backscatter echowidth of a.5 x 1A perfectly conducting ogive with a 0.05A thick material coating containing the properties Er = 3 -j5, r = 1.5 -jO.5.................... 84 4.9 Ez and Hz backscatter echowidth of a 1 x 4A perfectly conducting ogive with a.05A thick material coating containing the properties er = 3 -j5, pr= 1.5 -j0.5.................... 85 5.1 Comparison of the far zone and near zone fields for TM plane wave incidence on a circular metallic cylinder as computed by the total field FE-BI method and the eigenfunction series. (a) backscatter echo width vs. kao - the lines over the horizontal axis correspond to the eigenvalues of a circular conducting waveguide. (b) magnitude of the TM scattered field on the enclosure at the resonant frequency kao = 23.586............................. 92 5.2 Far and near zone fields for TM plane wave incidence on a circular metallic cylinder as computed by the total field FE-BI method with k replaced by ka in the BI equation. (a) backscatter echo width vs. ka,. (b) magnitude of the TM scattered field on the enclosure at the resonant frequency ka, = 23.586.................... 93

5.3 Far and near zone fields for TM plane wave incidence on a circular metallic cylinder as computed by the scattered field FE-BI method with k replaced by ka in the BI equation. (a) backscatter echo width vs. kao. (b) magnitude of the TM scattered field on the enclosure at the resonant frequency kao = 23.586................ 94 5.4 Bistatic echo width for TM plane wave incidence on a 5.256A x 5.256A metallic square cylinder with the boundary integral place at a radius ao = 3.754A............................... 95 6.1 The general axially symmetric surface with source (primed) and observation (unprimed) points and the corresponding unit vectors... 98 6.2 Cross section of a generating surface enclosed by the fictitious boundary.................................... 102 6.3 Typical triangular elements in the vicinity of the line singularity at Ro =! for real K.............. 106 6.4 A typical pair of source (e') and observation (e) elements associated with the discretization of the boundary integral equation...... 120 6.5 The configuration of a typical boundary element e' passing through triangular element e...................... 134 6.6 The mesh used for the frequency sweep of a conducting sphere is shown in (a) and the expanded region shows the values of ka at which the line singularity intersects the nodes. The frequency sweep in (b) demonstrates the inaccuracies associated with the line singularities. 140 6.7 A sphere coated with E, =, = 1 - j5 is swept in frequency.... 141 6.8 The axial backscatter cross section is displayed as a function of normalized radius k0oa. The structure is a conducting sphere of radius a coated with a dielectric with c, = 2 - j4 has an outer radius of ao = 1.02. Employing a complex wavenumber k = k0a, the resonance behavior is removed for the complex values in comparison to the c = 1 case.............................. 143

6.9 The bistatic scattering pattern for a coated sphere. (a) Shows the summation of modes converging to the correct solution for an incidence of 60 degrees. The converged solution for modes 0-10 at 60 degrees (plotted on a scale shifted by -60 degrees to coincide with the pattern due to axial incidence) shown with the axial incident result and a comparison to moment method results......... 145 6.10 The bistatic scattering pattern for a perfectly conducting sphere of radius 1A. Both E-plane (TE) and H-plane (TM) are shown along with data generated via the moment method............ 146 6.11 The bistatic scattering pattern for a 2A x 0. 176A perfectly conducting with a.0.05A coating with Er = 2 - j2 for an axially incident plane wave. Both E-plane (TE) and H-plane (TM) are shown along with data generated via the moment method............... 147 6.12 The mesh of the sphere-capped cone frustrum............ 147 6.13 The bistatic scattering pattern is shown for a plane wave axially incident on the coated sphere-capped cone frustum. Part (a) shows the E-plane pattern, while (b) demonstrates the H-plane pattern.. 148

LIST OF TABLES Table 3.1 Definitions for the finite element mesh................. 18 3.2 Definitions of the field vectors..................... 19 3.3 List of major memory-consuming variables.............. 35 3.4 List of major memory-consuming variables (continued)...... 36 3.5 Summary of major memory consumption............... 37 3.6 Summary of major memory consumption: filled mesh....... 37 3.7 Summary of major memory consumption: open mesh.. 40 3.8 A comparision of computation efficiency of the FE-CGFFT with the CGFFT method............................. 41 4.1 Definition of various quantities associated with the finite element mesh 58 4.2 Definition of various field vectors associated with the finite element mesh and its boundary....................... 59 4.3 Geometric quantities in reference to the figure above......... 67 4.4 Definitions for the finite element mesh with an ogival enclosure. 67 4.5 Definition of the field vectors for the ogival boundary........ 67 6.1 An algorithm for the boundary integral system assembly for a generic element interaction matrix [A].................. 123

CHAPTER I Introduction Over the years, integral equation (IE) techniques have been employed for evaluating the electromagnetic fields scattered from various structures. Special cases of the IE methods are surface integral equation (SIE) methods and volume integral equation (VIE) methods, which have been the workhorses in computational electromagnetics for the past two decades. They have been used to solve the scattering from simple structures, such as smooth perfectly conducting bodies and homogeneous material structures to more complicated ones such as multi-layered coated conductors. However, the classical SIE and VIE techniques result in fully populated matrices whose required storage is O(N2), where N is the number of unknowns in the model. Consequently, as the size of the structure increases, the memory demand may become much greater than the storage available on existing computing resources. To circumvent this, the IE may be cast in convolutional form by a proper choice of the discretization model. The resulting system can then solved via the conjugate gradient method, using the fast Fourier transform (FFT) to evaluate the discrete integral operator. The use of the FFT requires storage on the same order as the dimensionality of the scattering body (i.e., 0(Nld) for three-dimensional scatterers using VIE, where Nd denotes the number of cells per dimension), but is still lower than that required

by the fully populated version. This method, referred to as the CGFFT method has been recently employed by several authors to solve two-dimensional (2-D) and three-dimensional (3-D) inhomogeneous scattering problems [1, 2, 3]. Nevertheless, for 3-D applications the storage requirement remains problematic and furthermore, rectangular grid meshing results in a stair-case approximation of otherwise smooth scatterers. In an effort to further reduce the storage requirement, partial differential equation (PDE) methods were considered, since these methods provide for an O(N) storage requirement. Among them are the finite difference (FD) methods [4] and finite element (FE) methods [5]. Finite difference techniques typically suffer from the same staircasing problems as the CGFFT methods and may have difficulties in modeling numerical dispersion when solved in the time domain (TD). Frequency domain solutions (of primary interest in this study) allow for accurate (conformal) modeling of the scattering structure since, for instance, triangles and tetrahedrals can be used for discretizing surfaces and volumes, respectively. However, when the finite element method (or any other PDE method) is employed, the mesh must be truncated at a distance from the structure on which an approximate or exact boundary condition is applied. Various methods for truncating the FE mesh have been employed. Among them are the simple enforcement of an approximate absorbing boundary conditions (ABC) at the mesh boundary and the unimoment method. The ABC's are popular because they result in a banded sub-matrix structure. However, they require additional unknowns since the enclosure must be placed at a distance approximating the far field region. This distance may be large for structures with sharp edges and as a result a large number of unknowns is required, especially for 3-D applications. In the

unimoment method [6], the scattered field in the unbounded region is represented by an eigenfunction expansion. The coefficients of the expansion are then determined by enforcing continuity at the circular mesh termination boundary. This method produces a dense square sub-matrix whose dimension is proportional to the number of modes. It also requires the truncation of an infinite series which may be slowly convergent for irregular structures, thus resulting in a large storage requirement. An exact termination method is the boundary integral equation, first introduced by Hsieh [7] and McDonald and Wexler [8] in the early 1970's. The method is heretofore referred to as the finite element - boundary integral (FE-BI) technique and is the subject of this thesis. The boundary integral equation employs the freespace Green's function which implicitly satisfies the radiation condition at infinity. Since the integral equation is exact, the mesh termination boundary may be brought very close to the scatterer to minimize the meshing of the free-space regions enclosed by the termination boundary. The main drawback of the method is again the dense submatrix structure associated with the boundary integral equation. However, by carefully choosing the shape of the boundary, some of the boundary integral terms are cast into convolution form. Thus, on employing a conjugate gradient solver, these integral operators may be evaluated via the FFT as was done in the CGFFT method. The required storage of the boundary integral is then reduced toward O(N) and allows for the accurate solution of large inhomogeneous scatterers as well as small targets. In this thesis, the FE-BI technique (presented in general terms in chapter II) is developed for two-dimensional and axial-symmetric structures. The two-dimensional case is based on [9]. Chapter III contains a FE-BI formulation for rectangular enclosures, leading to simple boundary integrals some of which have convolutional form,

thus leading to a reduction in the memory requirement. In chapter IV, circular and ogival enclosures are considered. For the circular boundary, the boundary integral is entirely convolutional in form and an O(N) storage requirement is thus achieved at all times. Circular enclosures are consequently preferred for storage efficiency if the structure's outer boundary does not substantially deviate from a circle. In chapter VI, the FE-BI formulation is developed for axially symmetric structures. The finite element method for this problem was originally developed by Morgan and Mei using the coupled azimuthal potential (CAP) equations for generating the FE matrix system. In that formulation, the unimoment method was employed for terminating the mesh, whereas in our implementation the boundary integral equation is instead employed for this purpose. It is related to Flemings implementation [10] for circular boundaries, but the present formulation allows for an arbitrarily shaped enclosure. Consequently the boundary integrals are no longer convolutional and the storage reduction is in this case achieved by exploiting certain symmetry properties of the boundary integral subsystem. A problem with the CAP equations is the presence of a singularity in the finite elements which tends to corrupt the solution, and a method is presented for circumventing this difficulty. A difficulty with most boundary integral formulations is the appearance of internal resonances which corrupt the solution. These resonances correspond to the eigenvalues of the integral operator, and also occur at the same location as the cutoff frequencies of a resonator with conducting walls of the same shape as the mesh termination boundary and filled with the same material as the unbounded medium. The resonances are pare ticularly problematic when the structure becomes electrically large, in which case the eigenvalues are very closely spaced. Without any correction, computations for large structures are unreliable, and a method for correcting the

resonance problem is developed in chapter V. Specific contributions of this work include the implementation of the FE-BI for rectangular, circular and ogival boundaries as described in chapters III and IV in a manner taking advantage of storage reduction schemes. A new method is also presented for suppressing the resonance corruption problem existing in almost all implementations employing some form of a boundary integral equation over a closed surface or contour. This method involves the introduction of a complex wavenumber and is demonstrated for both 2-D and axial-symmetric bodies. The implementation of the FE-BI method for axially symmetric structures as described in chapter VI is for the most part new and incorporates a scheme for treating the line singularity associated with the CAP equations. This resulted in real matrices for lossless scatterers, a property consistent with 2-D and 3-D implementations of the finite element method.

CHAPTER II Fundamental Concepts In this chapter, some basic concepts are presented as applied to electromagnetic scattering and its computation via the FE-BI method. At the end of the chapter, the conjugate gradient matrix solver is also appended. 2.1 Basic Electromagnetic Theory The overall goal in any electromagnetic problem is the solution of Maxwell's equations subject to given boundary conditions. In the frequency domain, Maxwell's equations in free space are V x E = -jwiH (2.1) V x H = jwcE (2.2) V (EE) = 0 (2.3) V. (pH) = 0 (2.4) where E(r) = electric field strength [Volts/meter] (2.5) H(r) = magnetic field strength [Amperes/meter] (2.6) E(r) = electric permittivity [Henrys/meter] (2.7)

p(T) = magnetic permeability [Farads/meter] (2.8) Note that the last two equations, which are based on Gauss' law, can be obtained by taking the divergence of the first two equations. Thus for time harmonic fields, only the first two of Maxwell's equations are needed for a unique solution of E and H. By substituting (2.2) into (2.1) and vice versa, we obtain 1 2 V X -V X E - w2cE = 0 (2.9) V x -V x H - w2pH = 0 (2.10) which are commonly referred to as the vector wave equations. They must be solved and are subject to the boundary conditions on the particular scatterer as well as the radiation condition. Typically when solving (2.9) and (2.10) for a piecewise homogeneous scatterer, the following conditions must be explicitly imposed: * Continuity of tangential E and H * Continuity of normal pH and eE * Tangential E = 0 on metallic surfaces For radiation and scattering problems, the radiation condition must be also satisfied. This condition assumes that the wave must be outwardly propagating and must attenuate no slower than the inverse of the distance far from the source. Expressed in mathematical terms, lim [rH-r x E]-0 (2.11) lim [E-r x rH] = 0 (2.12) where r- = is the impedance of unbounded medium and i is the unit radial vector in spherical coordinates.

The media considered is this work is linear and isotropic and the material parameters y and e are often normalized to their free space (vacuum) values of Eo, 8.85 x 10-12 H/m and gso = 4ir x 10-7 F/m. The relative constitutive parameters of the medium will be denoted by er = (2.13) pr -- (2.14) /lo A fundamental solution of the wave equation is the plane wave E(r) e-jk= (2.15) where k = w/V and w = 2irf is the frequency in rad/sec. For lossy media, p and e are complex and thus k = kr + jki, implying an exponentially decaying field. Note that the required condition, ki < 0, ensures that the wave does not grow exponentially and, consequently, the imaginary part of p and e must likewise be less than zero. 2.2 The FE-BI Approach To solve PDE's (2.9) and (2.10), the scattering body is first enclosed in an artificial surface S bounding a volume V and for this application all sources will be assumed to reside in the unbounded region outside S. To satisfy Maxwell's equations in this volume, the vector wave equation for the electric field in (2.9) is discretized via the method of weighted residuals. Though not given here, a similar approach in terms of the magnetic field may developed in a parallel fashion. The volume V is first divided into Ne elements. Within each element, the weighted residual expression is given by W| w -(V x-V x E-w2eE)dV = 0 (2.16) Ve P

where W[ is the ith weighting function over the eth element. First, we temporarily define A as A=-V x E (2.17) and recall the identity W *.V x A = V. (A x W) + A. V x We (2.18) Using this identity, the first term in (2.16) can be rewritten as JfJWe V x AdVe = JJ (A x We)dVe + 7A V x WedVe Ve Ve Ve n=. (A x W~)dSe + A. x wedVe Se Ve (2.19) where the last equality was obtained by invoking the divergence theorem and thus, Se is the surface enclosing Ve and nfi is its outward unit normal. Substituting (2.19) into (2.16) yields fff k(-V x E) X (V x We) dV = e (1 V x 7) x WedSe (2.20) and upon rearranging this we have /1 [(!V.' x E) (V x W) )- W2-E] dV -= f jwW. (ni x H)dSe (2.21) Ve Se For scatterers comprised in part of conducting material, the surface boundary of the volume V is S = S, + S, where the subscript a denotes the exterior boundary and the subscript c denotes the conducting boundary. In the assembly of the finite element equations, the fully discretized form of (2.21) is summed over all elements. Since the tangential fields are continuous across material interfaces, the surface integral in (2.21) will cancel everywhere except on Sa and S,. Thus, unless Se intersects Sa or Sc, the surface integral may be removed.

To construct a system of equations from (2.21), the field in element e may be represented as n.= Nj(x,y,Z)Ej (2.22) j=1 m x x Nj(x, y, )E (2.23) j=1 m n x He = zn x Nj(x y z)Hte (2.24) j=1 where the superscript t denotes the tangent to the boundary of element e. (These vector shape functions must also be chosen to satisfy the divergence condition given by (2.3).) Substituting this expansion into (2.21) and using Galerkin's approach (i.e., W' = Ni) gives n ///1 -- ZE,JJJ [(V x iN) (V x N) - kN iNj] dVe 3=1 Ve IL = H f" jwNi (Zn x Nj)dSe + Ht fJ jwN. (n x )dSe (2.25) j=l Se E Sa j=l Se E Sc where k = wx/pe. After summing over all elements, our system takes the form {En} {El} Itn Ktt Kat Katt IKtn Baa 0 ~aa If aa I al ac If ac {E,) K1a Kta KI KIn I= 0 - 0 (2.26) IJ nn It In Knn int o I I I I o eec }ca }cl cc cc Ifn tIKtt K t Ktn Ktu 0 BcH~t ca ca I c In these the subscript a refers to unknowns on Sa, c refers to the set of unknowns on SI and I refers to the remaining unknowns. Furthermore, the superscripts n and t refer to the normal and tangential components, respectively, of the weighting and/or

11 testing functions. The application of the boundary condition on S, results in the following system Ifa aKn t I K'ani Kanc 0 0 0 {E'} Ictna Ktt I<, Kl a Ba a 0 | I Io {El) nla JEIa Ken KIc ~ ~ | | {E} |~ (2.27) KIfl, I'Ku Ki K ~I 0 0 0 Ia ca Ic{ cc 0 0 0 0 I 0 0 (HJt} Clearly, another equation is necessary relating E and Ht on Sa. The FE-BI approach employs the Stratton-Chu integral equation for the purpose. This equation relates the electric and magnetic fields via the integral equation E(r) = r(o) + 5I {-jkO [' X 3o0H()] g(n,) +,V' x )]V) Si +[' x E()] x )dS (2.28) where y/o = -Vo/eO and g is the free space Green's function. Evaluating this equation on the surface Sa and discretizing it, we obtain the subsystem Laa{Eat} + Maa{Ha} (2.29) which relates the tangential fields on Sa. Combining (2.29) and (2.27) we obtain the

12 block matrix IKnn Knt Jn Kt 0 {En} 0 aa ac I'tn Ktttt K a c Baa {Et} 0 K a IIa K I 0 {EI} 0 (2 {C~E,'} 0(2.30) IJnn IJnt Kn K~n 0 O {Ec } o 0 O O 0 I 0 {Et} O o La O O 0 Mat {H } { Ui which, after applying the condition on S, can be solved uniquely for the interior and boundary fields. The resulting system is then solved via the conjugate gradient method. After the fields in the solution region have been determined, the scattered fields are computed by evaluating the integral equation in the far zone. In particular a quantity called the echowidth is computed for 2-D problems and is given by [11] a= lim 27rp 2 (2.31) p-.oo2 IJi12 where p is the usual cylindrical coordinate and X denotes either E, or Hz. For threedimensional applications, the quantity of interest is the radar cross section (RCS) and its expression is r = lim 4rr (2.32) r-oo Ji(-)l2 where r is the radial spherical coordinate. Expressed in component form we have p = lim 47rr2 IEP( (2.33) =D (2.33) where p and q represent the polarization of the scattered and incident fields, respectively.

13 2.3 The Conjugate Gradient Method The CG algorithm is employed throughout this work and is given as follows for solving the system Ax = b. Initialize the residual and search vectors nib =)11 [ Oinc O O O T 112-1 b 112 s = A-(~) r(l)= b-s s = Aar(1)?s- =11 S 112 B(~) = )-1 Iterate for k = 1,..., Ng s = Ap(k) 78 = II' JII a(k) = 7;1 0(k+l) -= (k) + a(k)p(k) r(k+1) = r(k) _ a(k)p(k) 7r -11 r(k+l) 112 S = Alr(k+l) pl(k) = Y-1 p(k+l) = p(k) + p(k)S(k)

Terminate when k = N. or \/ < tolerance. The given tolerance must be used with caution. In each of the FE-BI techniques developed, the excitation b is only partially full. If the simple preconditioner of dividing by the diagonal value is employed, then those rows associated with an excitation component of zero will consequently not scale the excitation. In other words for those rows, b will not change in value while the matrix itself does. As a result, the residual error will be scaled differently and, consequently, the tolerance criterion will be scaled differently as well. This may lead to artificially small normalized residual errors and cause the convergence criterion to be satisfied before convergence is actually achieved. One way to avoid this is to compute the far field in terms of o at one or more angles every few iterations and compare the current value to the previous one. When the difference is sufficiently small for 15 or 20 consecutive iterations, convergence is assumed to have been achieved.

CHAPTER III A Finite Element - Boundary Integral Formulation for Two-dimensional Scattering with a Rectangular Termination Boundary 3.1 Introduction In this section an FE-BI method is developed for two-dimensional scattering in which the exterior boundary is rectangular in shape. This enclosure ensures some of the boundary integral terms are convolutional and are therefore amenable to evaluation via the FFT in the conjugate gradient solver. Results are presented for several structures. 3.2 Analysis Consider a cylindrical body of arbitrary cross-section and composition illuminated by the plane wave ~, (p) Z= i9Z/C(T) -.ejko,P cos(O- o) (3.1) as indicated in Fig. 3.1. To evaluate the fields scattered from this object, two boundaries are placed tightly around the body as shown in Fig. 3.2. Inside the outer boundary, the Finite Element Method is applied to solve the Helmholtz equation

16 y x Figure 3.1: Geometry of the scatterer 2 i4 3 T Figure 3.2: Partially discretized body

given by V. [v(p)VO(q)] + k2u(p)q(p) = O (3.2) where 0(p) - E,() (3.3) v = t() (3.4) u(i) = C,(P) (3.5) for E-polarization and =(p) Hz( ) (3.6) v() = (3.7) u(=) - r,(T) (3.8) for H-polarization. Also, ko = wr/i is the wave number, and,r and Er are the relative permeablility and permittivity, respectively. The appropriate boundary condition is enforced on the surface of the impenetrable boundary, while the radiation condition is satisfied implicitly by evaluating the integral equation (P) = -'"(n) - {G( P ) [')]- ) G( P )] } dl' (3.9) on the boundary rF, where G(p,) = -~H2)(ko- 1) (3.10) is the 2-D Green's function in which Ho2)(.) denotes the zeroth order Hankel function of the second kind. Furthermore, p and A' are the usual observation and source position vectors, respectively, and ~P - 61 = J~C l.............

Symbol Description Ng number of nodes in the finite element mesh Ne number of elements in the finite element mesh Nx number of nodes on ra or rb along the x-direction Ny number of nodes on r, or rb along the y-direction Na total number of nodes on ra Nb total number of nodes on rb Nab Na + Nb Pa Et=l rai rb El41 rbi r Ei l rc, (Xa,,yai) coordinates of a point on contour F,a (xb,, Yb,) coordinates of a point on contour rb, Table 3.1: Definitions for the finite element mesh The normal derivatives are taken in the direction of the outward normal of r,. 3.2.1 Discretization of the Object and Field Quantities In Fig. 3.2, Pa is the field/observation point boundary, and rc is the integration contour, which is placed midway between ra and rb. Also, rd denotes the contour enclosing the impenetrable portion of the scatterer. Herewith, each side of Pa, rb and rc are numbered counterclockwise starting from the top side, as indicated in Fig. 3.2. The fields in the region between ra and rd satisfy (3.2) in conjunction with the required boundary condition on rd. The boundary integral equation (3.9) will be enforced on ra,.

Symbol Description ai nodal fields on Fat bti nodal fields on Fb,, gnodal fields on Fb, qd nodal fields on rd ~1 sregion enclosed by rb and rd, exclusive Table 3.2: Definitions of the field vectors To numerically solve (3.2), it is required that the region within Fa be discretized. This is done in a traditional manner when employing the finite element method. The global node numbering starts from the right endpoint of contour ral and proceeds counterclockwise. The numbering continues beginning at the right endpoint of contour b1, and proceeds counterclockwise. Within rb, the nodes are numbered arbitrarily. The definitions pertaining to the FE mesh are given in Table 3.1. Each node corresponds to an unknown field value to be determined. It is important to associate the unknown field values corresponding to the various node groups on contours ra and rb by using different variables. The labeling scheme is given in Table 3.2, and this discrimination of the nodal fields is required, since they are treated differently in the analysis. 3.2.2 Derivation of the Finite Element Matrix One of several approaches may be used to generate the finite element matrix, such as the variational approach or the method of weighted residuals. In this development, we will utilize Galerkin's method, which is a special case of the method of weighted

residuals with the distinction that the testing functions are the same as the basis functions. Proceeding with the finite element analysis, we may rewrite (3.2) as ax [(x Y) (x, Y)] + 4 y [v(xy)y (x, y)] + kou(z,y) (x,y) = O (3.12) the residual of which is given by R= + [V( xY)xz(xy) y (,y) - k2u(x,y Y)q(x,y) (3.13) The field within r, may be represented as a summation of piecewise continuous functions and, thus, may be written as Ne (x*, y) = e(x, y) (3.14) e=1 where Oe(x, y) is the field within element e. It is expanded as n e'(X y) _ Nje(x, y) dj (3.15):=l where Nie's are the standard shape functions (found in any standard finite element book), he's are the fields at the nodes, and n is the number of nodes per element (n = 3 for linear elements considered here). The weighted residual equation for the eth element is defined by JJ NR dxdy = 0 i = 1,..., n (3.16) Se where Se denotes the surface area of the eth element. Inserting (3.15) and (3.13) into (3.16) yields rff N? -[_ v ( ]ne - (y v O ) k2uNc qj dxy -d O i 1,2,.., n (3.17)

21 Further, by using the identities ta ( N\a ) a ( ae a\aNtaN Nt yy / y y y Y ay ay (3.19) and the divergence theorem x + dxdy =- n(u + v) dl (3.20) where Ce is the contour enclosing the eth element, (3.17) becomes / N a a aN v aNeN j +V k2uNee)i qe dxdy sx ax ay ay 0 3/3 j= i 3 v v ) n i 1 (3.21) This may be written in matrix form as Aee = be (3.22) where aNP aM aNe aNe [AeLi= If ( S 3+v 3- L dxdy (3.23) avx ax -ay ay-y 3 y and lb} = Nv' h dl i -, 12,..., n (3.24) For linear triangular elements, Ne are given by 1 xE yl Qe 1 y ( b1 +..... (3.26) 1 x (bc- bjc) 1 4 and~~

22 e = e e - e e (3.27) ai Xjk - kyj (3.27) be = y_ - y (3.28) c = Xe - x (3.29) where (xe, ye) is the coordinate of the ith node of the eth element. From (3.25) aNye_ bi (3.30) &x 2Me 9Ne c? ay 2Qe Substituting these and the formula [5] (N1 )P(N2)q dxdy = 2Qp P!ql (3.32) se \~r2 I C~L UY - J (p + q + 2)! se into (3.23), we find Ve e e +"e [Ae]ij = 4 (bb + cc) - ko2e (1 + j) (3 33) where 1 if ij iij = (3.34) O otherwise In (3.33) we have assumed that u and v (the material constitutive parameters) are constant within each element and are given by Ue and Ve, respectively. By summing over all elements as implied by (3.14), we may write the overall system in block form as Aaa Aab 0 0 Oa ba Aba Abb AbI 0 | b (3.5) O Alb AII Aid kI O 0 0 AdI Add'Ad 0

23 In this, the values of the elements in the submatrix Apq are the contributions associated with the nodes in group p which are connected directly to the nodes in group q. One can easily show that the line integral contribution (3.24) of those elements vanishes everywhere, unless the element has a side in common with,Fa. As a result, be contributes only from the boundary Fa of the finite element region, as indicated by the presence of the vector ba in (3.35). Without a priori knowledge of the total field on that boundary, ba cannot be determined. We may, however, provide the appropriate condition on this boundary by utilizing the integral equation (3.9). This amounts to replacing the first block-row of the matrix (associated with Oa on Fa) with a discrete form of this integral equation. 3.2.3 Evaluation of Boundary Integral The boundary integral in (3.9) may be written as a summation of four integrals, one for each side of the contour F, as = G) - p$ [( G(T, p' dl -fr, [(, 7) —(') + p- (p) P( C dany 3 GI Tj 1) 0 (1) + G01;51 ) dl C -L rc.da, a 1 C2 anx z -~,[ I\IG(;\ q$(')-G(T5 d JrG4 a - Ox' dl 4 (3.36) C4 L ~~nx O(P O~g) axl C4 where the derivatives along the x and y directions, denoted by a a re ian., (ny, e spectively, have been left in this form for the later convenience of determining them numerically. More explicitly, we have q~(X, yj) = ~tlc(X, ya i) 3 3

- r [G(x - X',ya,,y ) Yc (, YX) - ( x',yc ) G(x - x',Ya Ycl)] dx' - [G(, c2,ya,Y )n ((c2,Y ) + b(zx2 y')YYG(x, c,,ya.i y')] dy' C2 3 3n, OX' _-Jr G(X - x,7 Ya YC3) O(X',YC3) + q(x'YC3)2 G(x - x', Ya, c3) dx' C3 3 )3 - - [G( C4 )Ya, yy') q(xc4,') - ((xC4,) y') G(x, C4, Ya I, y)] dy' (3.37) and O(Xa,2, Y) ='c(Xaa2 Y) 4 4 -r [G(a, X, y, yc,) (X, YI) - (xx YcI) G(xa2, y ycl dx 4 G4xan, x4 2 y - y' dy a a r G(Xa 2' y - ) xn y'(xc2 Y) + q$(x 2 y'> G(Xa2, yXyc' ) dy' - f G(xa,x, y, YC3) (x YC3) + ( YC3) G(a,,2 Y C3) d'C3 [4 Of31 a a - [G(Xa2 X c4, Y - ) X(C4,yY') y') - (Zxc4, y') G(Xa2 Z C4, y-y') dy C4 4n, 4 (3.38) where the first subscript on x or y refers to the contour (a, b or c), and the second refers to the contour number (see Fig. 3.2 and Table 3.1). It is noted that the arguments of the Green's functions have been modified to imply a convolution when appropriate, and this representation will be used throughout. The normal derivatives of g may be evaluated via the central difference formulas O(Xc y') = [b(xa, y') - ((xb, y')] + O(A2) (3.39) an((x, yc) = [q(x', y)- q(x', Yb)] + O(L2) (3.40) where A is the displacement of rFa from rb (A is usually less than one tenth of a

25 wavelength). Substituting (3.39) and (3.40) into (3.37) and (3.38), we obtain O(X, ya ) = OinC(X, Ya ) 3 3 - frK [7G(- X',Ya Yc1)q(X', Ya - IG(z - X', Ya,, YC, )(x', Yb )] dx' C i 3 3 1 3 3 -r [f I G(x- X', ya,yc3)+(x,ya3 )- K-G(x -',YaI,yc3)1(x', yb3)] dx' - KG [I ( xG(, ZC4,Yai, y')(Xay y') -KG(xXc4ya, by') dy' ~' L'C4 3 (3.41) and (Xa2, y) Y= inc(Xa2, y) 4 4 -r I [KI+G(Xza2 X,', y, YCI )y(X', Yal ) - K G(, )('y,ybl)] dx' -/r [I' G(Xa 2 xC2 Y - Y') ((Xza2, y') -I +G(Xa 2 c2 x Y-Y ) (xb2 aY')] dy' -IJrC3 [I 4G(Xa2 X y, Yc3 ) i(x', Ya3) - I(G(Xa2 y YYC3)q$(X'Yb3)] dx'' f 1( G(Xa 2, Xc4, y y )(Y')O(Xa4 Y')- YG(a2 Iv C4 y y )(Xb4Iy) dy' 4 4 (3.42) in which 1 i0 K = ~ 2 ax' (3.43) 1 10.+_ - 1 0 (3.44) v A 2 0y' Assuming a pulse basis expansion for the nodal fields (i.e., piecewise constant functions centered at nodes on contours rP and rb), a midpoint integration may be performed for the evaluation of the integrals in (3.41) and (3.42), to obtain (xi, Ya i) = /fc (Xi, Ya ) 3 3

26 - z/i' [Ks G(i,, yX, ya,Yj)~(Xa2, yj)- ~ 6(i, xyc2,, ya j)1 (~, yj -nf~ ~ G(x - xj, ya, i, y,Y)d(xj, ya) - I(',jG(x - X,,Ya i,Yc3)(X3j,YbY)] j=1 3 3 - j -1 [KI G(xiXC,yYai,yj) (Xa,yj) - KxYiG(Xi, x c4 Yai Y,y)(xb4,Yj)] j~~~l ~3 3 (3.45) and (Xa2, Yi) = fC (Xa2, Yij) A- [I -G(xa2 zxj, yi, Yc)(Xj, I Ya1) - K+G(xa2 Xj, Yi, Ycl )q$(xj, Ybl)] - ~\ [Iii iG(Xa2z Xc2, Yi -yj)q(a2, yj)- IjG(Xa2,Xc2 - yi-yj)q5(Xb2, Yi)] 3= 1 - Ajj [IxG(x r2 Zi Yi, Yc3)(j, Ya3)-IG(a2 Xj, Yi, yc3)(X Yb3)] - ~ [I G(X~a2 Xc4 7Yi -Iyj) O(Xa4,,yj) - KG(Xa2Xc4, Yi Yi)k(Xb4, yi)] (3.46) In these xi and y, denote the ith matching/testing points corresponding to the nodal locations on Fa, while xj and yj denote locations on rb. We recognize some of the terms in (3.45) and (3.46) as discrete convolutions amenable to numerical evaluation via FFT. The subsystem (3.45) and (3.46) may be written more concisely as O(al 2 7 y oi'c (SXb T 2 S b3R Tb4 al 4a 4 | K|Tinbc T Sb2 T+b3 S Rb4Y | a2 a3 St 2bl TaXb2 Sab3 Tab4 Xa3'ea4 S ftTR S S2b Ta ab4 -raa3b4 ab ab3

27 Sl bl Tal b2 Salb3R Ta b rb Ta+2bl a2b2 Ta2b3 Sa2b4 Y b2 S3blR Ta3 b2 S3b3 T+ b4b3 a b Sab2 R Ta b3 a+ (3.47) with the various parameters to be given explicitly later. The matrices R,,y simply reverse the order of the unknown vector so that the convolutions may be performed properly. This is required solely because of the employed counterclockwise nodal numbering scheme. Since ") (') /~X(3.48) (bi) last O= (bj first j = 2, 3, 4, 1 the vector T |(bl'k b, 2'ib3 rbk4 ] (3.49) can be related to the actual unknowns on the contour Fb via a transformation Db as )b- =Dbqb (3.50) and (3.47) may then be written as (I + Laa)qa - LabDb qb =- di (3.51) or [AB,] | = | $fC (3.52)

28 where I + Sab, Tb2 Sb3 R Tal b4 T~,l I+ r + T + Sa I + Laa Tb I+S 2 Tb3 S2 (353) Sa blR T3 b2 I+ b3 Ta b4 Tabl Sab2 RY T4 I + Sab4 and T+b Sab2 T2b3 Sa2b4 Ry Lab = S+ Rx S- T+ SbR Ta3b2 a3Sb3 aTb4 T+b Sa b2 Ry Tab3 Sai b4 (3.54) After replacing the first block row of (3.35) with (3.51), the complete system may be finally written as I + Laa -LabDb 0 0 Oka acifC Aba Abb AbI 0 ||b 0 (3.55) O AIb AII Aid'I 0 0 0 AdI Add'kd 0 to be solved via the CG algorithm. The elements of ABI defined above may be evaluated via the discrete Fourier transform. Specifically, we have Sb Obl = DFT-i {DFT[G(xYa, Yb,, Gya(xy, y ab )Y]DFT[bl ]} (3.56) 3 3 3 3 3 Sa3b, qbl = DFT' IDFT[G(,Y, ybi) ~ G(x, yaY )]DFT[+b l] } (3.57) 3 3 3 33 S b -b2 DFT' {DFT[G(x a2 xb:2, ):)+ G(xa2, b, y)]DFT[q>b2]} (3.58)

29 S=4b2 Ob2 = DFT- {DFT[G(xz ZXb2, y) + G(Xza4, Xb2, y)]DFT[b2]} (3.59) in which DFT denotes the discrete Fourier transform operator. Also 0 4 G(z, x', y, y') = 4 H2)(kolp-p ) = H (ko (x-x)2 + (Y y)2) (360) Gx~xi 0 G(x, I y, y') Gx(x,x',y,y') = -2-G (x'xyy) 2 ax' kH2) (ko/(x- x')2 + (y - x') (3.61) 8/ ( - x')2 + (y - y' ) Gy(x, x', y,y') = A G(x, x', y, y') - H o(2) ( ( -x')2 + (y- y)2 (3.62 -3,Ako,- __1 (Y - y') (3.62) - x')2 (y - y, )2 Special cases of the convolution operators for the chosen mesh are given as Gx(Xa2, Xb2,Y y=) 4 4 I(2)(koZXa2 - Xb2 + (y - y3)2)2) k f(Xa2 -Xb2 )2 + (y _ y,)2 4 4 (a2 4 14 Gy(x - X', y,,ybl) = 3 3 H(2) (ko(X - X)2 + (Y - b )2): o8 (x - x')2 + (ya, - yb )2 3 3 Y 3 3 and the corresponding expressions for G are implied by the arguments in the previous two equations. Additionally, the upper coordinate in the braces corresponds to the upper sign and likewise for the lower one. The cross-term element submatrices are given by [Tdl1b] = G(Xi, b2, YlaYi) + GX(Xi, Xb2, yal, Yj) (3.65) 3: ij 4 3 4 3 [TAbl] = G(zxa2,XZj, Yi, Ybl ) Gy(xa2,,Zj, Yi, Ybl ) (3.66) 3~ ij 3 a 3

30 with Gx(xi, Xb2, Ya I, Yj) = 4 3 H:(2) (ko.Xi - Xb2)2 + (Ya yj)2 t3 kok; -?i X- IA (3.67) 8 ~ Xi -- Xb2 )2 + (Ya - Y-)2 i: 36 4 3 and Gy(xa2, j, Yi, yb) )' 4 3 i(2)(ko(2 - a )2 + (Yi - Yb)2) 4 3 -jIyi/- yb, IA yb3 (3.68) 8 H(k(a2 - Xj)2 + (Yi — YbY 3 Ybl where again the corresponding expressions for G are implied by the arguments of the previous two equations. Making the substitutions (xi - Xb2)2 = (i - 1)2A2 (3.69) 4 (Ya - yj)2 = j2A2 (3.70) 3 3 Xi Xb2 )2 + (Yai - yj)2 = A2 1)2 + j2 (3.71) and (Xa2 - xj)2 = j2A2 (3.72) 4 (Yi - Yb )2 = (i _ 1)2A2 (3.73) 3 I(a2 -Xj)2 + (Yi-Ybi)2 = j2+(ij1)2 (3.74) 4 3 we may write Gx and Gy as (2) 12 +j2 1' 1' 2 7j,2 4G; 8x(-i) + X (3.75) Xb H Ya iY) = T3(~' ) 1i {b Gy(xa2,, xji, Y yb1 ) =:k (k 2 +(i - IY- (3.76)

31 to be used in the actual implementation of the system. Since each of the above relations are similar, we are required to store only one of them and alter the signs accordingly. It should be noted that, however, this is not the most efficient method of storage. Storing only a few of the cross term values and using an interpolation scheme will reduce the storage considerably. Of course, an interpolation table of (3.75) and (3.76) will lead to a substantial reduction in memory at the expense of some computational efficiency. Assuming that the positive sign is chosen in equations (3.75) and (3.76), we have T+ a2bl -'a2 abl 4 4 alb2 a Ib2 3 3 T+ -T+ a 2 b3 a 2 b3 4 4 ab4 - ba Ib4 3 3 (3.77) Choosing the positive sign for the (3.63) and (3.64), we also find s+ -s albl albl1 3 3 a2b2 Sa2b2 4 4 a3bi a3b 3 3 s+ -s+ a4b2 = a4b2 4 4 (3.78)

32 Thus, the elements of ABI may be written as I + S a1 bl b2 Salb3 R fal b4 T_ + I+S S R I+Laa - a2b3 a2b4 Y I + LaaSa3bRz Tb3b2 I + Sa3b3 a3b4'T__ S+ Ry T+ a4 bl ab2R a4 b3 I + Sb4 S bi al b2 Sa b3RX al bs Taabl Sa+2b2 a2b3 Sa2b, Ry Lab Sa3bR fa3b2 Sa3b3 a3b4 Ta bl Sa4 b2 RY a b3 Sa4 b4 (3.79) The elements of the adjoint of ABI required in the implementation of the CG algorithm are I+(Srbl,)a (2b ) R(Sa3bl)a (7 bl )a (I +Laa) = | (Salb2)a I + (Sab2) (a3b2) R (Sab2) (Taj b+ )a Ry (S 4)a (bT+)a I + (S- )a [ 1 0(Salb b ) (<a~bl) Rr:(S,;b, )a (S - )a (LabDb) = |(,b2)a (S2b2)a (7Lb2)a Ry (S)Ob2)a RT (Sa b3)b3 (7a;b3 ) (Sa3b3) ( b3)a (7+) b )a RT jS2 T)a (+)a (S- )a ( alb4( Ry (Sa+ b a3b4 I (Sab b4 (3.80) 3.2.4 A CGFFT Implementation The conjugate gradient algorithm presented in section 2.3 is employed for solving the FE-BI system (3.55). The required numerical computation of Az and Aaz are

33 performed by first decomposing the system into a summation of two matrices; one involving operators associated with the boundary integral and another involving the elements of the finite element matrix. Then system matrix-vector products may then be written as {S} = {S}BI + {S}FE (3.81) where I + Laa LabDb 0 0 Z1 0 0 0 0 Z2 {S}BI = (3.82) 0 0 0 0 z3 0 0 0 z4 and O 0 0 0 Z1 {S}FE = Aba Abb AbI 0 Z2(3.83) 0 AIb All Aid z3 0 0 AdI Add z4 For the adjoint operations, we have I+ La 0 0 0 Z1 DTLab 0 0 0 Z2 {S}BI- (3.84) o 0 0 0 z3 0 0 0 0 z4 O O O O Z4 and 0 Aba 0 0 Z2 {S}FE = 0 Ab Ab 0 (3.85) O Ab, A,, Ad, z3 0 0~ A Ad Ad z4

In each case, the operation is performed such that only the resulting vector {s} need be stored. 3.3 Computational Considerations The FE-BI method for rectangular enclosures is efficient in terms of memory usage and computation time, and each of these aspects is discussed in detail below. 3.3.1 Storage Efficiency The fundamental advantage of this method is the reduction of storage requirements, so that the scattering by electrically large bodies may be evaluated. To show that the low storage requirement of O(Ng) is assured, we refer to Tables 3.3 and 3.4. These contain a list of all major memory consuming variables. A summarized list is also given in Table 3.5. Specifically, Table 3.5 includes the memory requirements pertaining to the finite element matrix (FE), fast Fourier transforms (FT), boundary integral cross terms (Cross) and conjugate gradient variables (CG). N, is one less than the number of elements connected to a particular node, and a typical value of 5 is used here. To put the quantities of Table 3.5 in terms of N., the total number of nodes, we consider two special geometries. The mesh in Fig. 3.3 corresponds to a penetrable body, while that of Fig. 3.4 corresponds to an impenetrable structure tightly enclosed by the picture frame. Within each special case two extremes are considered; a mesh corresponding to a square object (N, = Ny) and a long strip (N, >> Ny). In each case, Ng is assumed to be large. Alluding to Table 3.6 the total storage is O(Ng) for the square region, but is somewhere between O(Ng) and O(N2) for the (N, >> Ny) case. This is based

35 Memory Consumption variable typeI count ] comment Mesh Xg R N9 X coordinate of global nodes Yg R N9 Y coordinate of global nodes Nglob I 3Ne Node-element connectivity Pointers ABint I Nab Observation and integration points Pnodes I Pnum Nodes on conducting bodies dmatl I N - Nab Element material properties Finite Element Matrix (FE) Ar C (N+i )(Ng - Na) Non-zero values of FE matrix col I (N )(Ng - N,) Column pointer of FE matrix row I N - Na Pointer to rows of FE matrix Conjugate Gradient (CG) Phiz C. Ng Unknown vector CG1 C Ng Residual vector CG2 C Ng Search vector CG3 C Ng Temporary vector q C MAX(Nx, Ny) Temporary vector phiinc C Na, Incident field vector Table 3.3: List of major memory-consuming variables

36 Memory Consumption (continued) code variable type count comment Fourier Transforms (FT) FTHxl C 2Nx Fourier transform along x-direction FTHx2 C 2Nx FTHx3 C 2Nx FTHx4 C 2Nx FTHyl1 C 2N, Fourier transform along y-direction FTHy2 C 2Ny FTHy3 C 2N, FTHy4 C 2N, FT C 2MAX(Nx, Ny) FT of unknown sub-vector WR R 2MAX(Nx, Ny) Temporary array WI R 2MAX(Nx, Ny) Temporary array Cross-Term Matrices (Cross) PQP C | MAX(N, Ny ) PQm C |MAX(NN, Ny) Legend R = REAL*4 C = COMPLEX I = INTEGER*4 Table 3.4: List of major memory-consuming variables (continued)

37 Major Memory Consumption (N, > Ny) DItem [Type J Count FE COMPLEX (N2 )[Ng - 2(NS + Ny)] FT COMPLEX 12Nx + 8Ny Cross COMPLEX 2Nx2 CG COMPLEX 4N9 Table 3.5: Summary of major memory consumption Major Memory Consumption: Penetrable Item N=.Ny Nx >> Ny FE (Ne+1)(Ng -44) (Nc+I )Ng(l - 2) FT 20/Ng 12Ng/(Ny + 2) Cross 2Ng 2(Ng/Ny ) CG 4Ng 4Ng total I9N9 2() N2)2 + 62Nfi + 7Ng Table 3.6: Summary of major memory consumption: filled mesl

38 Figure 3.3: Example of the mesh of a penetrable structure Figure 3.4: Example of the mesh of an impenetrable structure

39 on the assumption that all cross terms are individually stored, but by using an interpolation table, the O(N9) memory requirement can be assured regardless of the value of N, with respect to Ny. In Table 3.7, more dramatic results for the storage of the cross term are listed. In this case, all of the unknowns are on the outer two boundaries, so it appears that the storage is O(Ng) for the square case. One must note, however, that the factor multiplying the N9 term may be quite small. The strip case, on the other hand, requires an O(N2) storage. This case would be an unlikely candidate for the use of this method, since that structure would be handled much more efficiently via a direct implementation of the CGFFT method. As noted above, the storage of the cross terms may be brought down to O(N9) for all cases by using an interpolation table, and this will certainly be necessary in a 3-D implementation.

40 Major Memory Consumption: Impenetrable Item_ Nx = Ny N~ >> Ny FE ( +)Ng/2 ( )Ng/2 FFT 5N /2 3Ng Cross N2/32 N2/8 l CG 4Ng 4Ng total N72/32 + 8Ng K N2/8 + 17Ng/2 Table 3.7: Summary of major memory consumption: open mesh 3.3.2 Computational Efficiency Since the primary importance of the FE-BI method is storage reduction, a comparable level of efficiency with alternative methods is a bonus. A method for which a fair comparison may be made is the standard CGFFT. This requires a 2-D FFT, which is slower than using multiple l-D FFTs for large bodies. We compared the two methods for a specific penetrable scatterer using an Apollo 3500 without code optimization. The pertinent CPU times are compared in Table 3.8. The comparison provides only a relative measure of the speed difference. 3.4 Far Field Computation The scattered fields may be computed as =-r {G(~) [n )]- ) [dG(P')] } dl (3.86) Y~,=-J~. ~~(p~ P ) [ dn' (~8'

41 Body Properties FE-CGFFT CGFFT Composition Dimensions T/I (s)]I li Total T/I (s) I Total dielectric 2A x 2A 8.6 155 1333 170 33 5610 er -4 -j.1 Legend T/I = time/iteration I = number of iterations Table 3.8: A comparision of computation efficiency of the FE-CGFFT with the CGFFT method Using the discretization scheme developed earlier, we have 5(x, y) = -{| [I;cG(x,x', y,yc,)k(X', yai) - IKG(x, x', y, yci )(x Yb )] dx' + [if+ IfrG(, x2,yy)q+(xa2,y')-K- G(x, x2, y, y')(xb2,y')] dy' + j| [IfK+G(x,x',yyc3)q(xya3) - Ky G(xz, x',yc, (3)(xyb3) dx' + [KG(x, X4, y, y')q(Xa4, Y') - K+G(x, x, y, y')q$(Xb4,')] dy'} (3.87) where the definitions for I~~ and Is are as specified previously. Letting (x, y(, YYc) = ~ J_ - G(x, x', y, y)dx' (3.88) fY(x, xc, y) = G(Z,,, y, y')dy' (3.89)

42 -X(zXlyYC)= G(x, x', y, y,)dx' (3.90) 1 Y' + A -9 (,ZY)1= - G(x, xc, y, y')dy' (3.91) (3.87) becomes $3(x, y) = NT +- ([/ x'(x,yY) +)- (x,YYcy)] {,al2}j - [/ (xYyc) + 7Y(x,Yyc)] {kbl2}j) j=1 Ny + Z ([ (X,, yC) + YX(, C, y)] { ) - a2} - [(X, XC,21Y) y y (X IY)] {b2b}j) j=1 N= + E ([/ (X,,yYc3s) +'Y(X, y,ycs)] {fba., }j - [/2(XYy,Yc3s) -- X(X,y,Yc3)] {~bb3}j) j=l + 1 ([13 (X, Xc4,Y) - 7Y(X, Xc,,y)] { Oa, }j - [P3Y(X Xc, x y) + -Y (X, Xc 4, )] Ob, }j) j=l (3.92) valid for all observation points (x, y). To specialize (3.92) to far zone computations, we must introduce the appropriate asymptotic expansion for the Hankel functions implied in (3.88)-(3.91). In doing so, we have X(X, Y, Yc ) = jfo(p)f (0, yC )ejko~x cos (3.93) 3 3,Y(x, xC2,Y) = jfo(P)f2(O, xC2 )ejk~oyj sin (3.94) 4 4 7(x, y, Yci ) = -fo(P)fl (0, Yc, )koA sin Oejkojx cos (3.95) 3 3 ^Y(Zx, xz, y) = -fo(P)f2(o, X2 )ko,~ cos OeJkoyj sine (3.96) 4 4 where fo(p) = l 4 e2-jkop (3.97) jkyc 1 sin G fi (O, Ycl) = -e 3 (3.98) f2(O,x 2) = -e korc (3.99) 4

43 in which (p, 0) imply the usual cylindrical coordinates of the observation point. Substituting expressions (3.93)-(3.96) into (3.92), we obtain q$ (x, y) -- fo(P) (Li + koA sin 0] {aj - [j - koA sin 0] { b l ij) (0(, Yc )ejkox cos~ Ny + E ([j - koA cos 0] { }j - [j + koA cos 0] {b2 }j) f2 (0, X2 ) k in j=1 Nz + ([j - koALsin0] {a3})j - [i + koALsin0] {Ib3j)fi(08, y3)ejk~2jc~sO j=1 + ([j + koA cos 0] {(q.,ai} - L -koA cos ] {b4, }j) f2(O XC)ejkoyj sin j=1 (3.100) From (3.100) the echowidth becomes a = 1k: (L[ + koA sin 0] { (a, }j - [j - ko sin 0] { bb, }j) f, (0, Yc, )ejko~x cos8 4ko j=1,Ny + ([j - koA cos 0] {a2 }i- Lj + koA cos 0] {b2}j) Xj)f(2)ekyj csin j=1 + Z ([j - koA sin 0] {Oa3.- [j + koA sin 0] {b3 }j) fl (0,yc3)eiko~j Cos -1 NNy 2 +. ([j + koA cos 0] 0 U)a. i - koA cos 0] {4b, }j) f2 (0, xc ) eikoyj sin j=1 (3.101) 3.5 Code Validation The scattering patterns from several rectangular structures are presented. The echowidth is computed for each structure and compared to the results of the moment method. The bodies are discretized at a sampling rate of 20 samples/free-space wavelength. Results are presented for the following cases:

* perfectly conducting bodies (Figs. 3.5 and 3.6) * partially and fully coated perfectly conducting cylinders (Figs. 3.7 - 3.12) * composite body (Fig. 3.13) In each figure, the discretized geometry is shown along with the corresponding results. As seen in all cases, the generated patterns using the FE-CGFFT formulation are agree with the moment method data.

45.25 x 2. X Perfect Conductor (E-pol) FE-GCFFr 20.0 OM 10.0 - 0.0 cZ -10.0 - CZ CP -20.0 -30.0. I....... 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 3.5: E_ backscatter from a.25 x 2A body

46.25 x 2. X Perfect Conductor (H-pol) 30.0. I. I.. I I I I I I I FE-GCFFr 20.0 10.0 0.0 -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 3.6: Hz backscatter from a.25 x 2A body

47.25 x 1. X Coated Perfect Conductor (E-pol) 30.0 FE-GCFFT 20.0 --------- MOM 10.0 0.0;CZ E -10.0 m -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure3.7: EP backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties, = 5. - j.5, p, = 1

48.25 x 1. X Coated Perfect Conductor (H-pol) 30.0 FE-GCFFT 20.0 30.0....,..... MOM 10.0 0.0 co -10.0 mP -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure3.8: Hm backscatter from a.25 x 1l perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5, lr- = 1

49.25 x 1. X Coated Perfect Conductor (E-pol) 30.0 FE-GCFFT ~ m 20.0 MOM 10.0 - 0.0 - cur -10.0 -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) material coating containing the properties Er = 5. - j.5, hr = 1.5 - j.5

50.25 x 1. X Coated Perfect Conductor (H-pol) 30.0 FE-GCFT 20.0 - MOM 10.0 0.0 cu - - -10.0 m -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure3.10: H, backscatter from a.25 x 1l perfect conductor with a A/20 thick material coating containing the properties -. = 5.- j.5, = 1.5- j.5

51.25 x 1. X Double Topped Conductor (E-pol) 30.0.... FE-GCFFT,- 20.0 20.0 o MOM 10.0 0.0 -10.0 -2o.o I - I -30.0 -90.0 -60.0 -30.0 0.0 30.0 60.0 90.0 Angle (deg) Figure 3.11: E, backscatter from a.25x IA perfect conductor with two A/20 thick top material coatings. The lower layer has the properties c, = 2.- j.5, #u = 1.5 - j.2, and the upper layer has properties Er = 4. - j.5, /tr = 1.5 - j.2

52.25 x 1. X Double Topped Conductor (H-pol) 30.0.... FE-GCFFT 20.0 - 10.0 0.0 -10.0 -20.0 -30.0..... -90.0 -60.0 -30.0 0.0 30.0 60.0 90.0 Angle (deg) material coatings. The lower layer has the properties E I = 2.- j.5, L, 1.5 - j.2, and the upper layer has properties = 4. - j.5,L = 1.5- j.2

53.5 x 1. A Composite Body (H-pol) 30.0.... I I.. I I I.. FE-GCFFT 20.0 10.0 ~O 0.0 CZ -10.0 C: -20.0 -30.0 i 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle (deg) Figure 3.13: H, scattering from a composite body. Both the perfect conductor and the dielectric body are A/2 in each dimension. The material properties are T = 5. -j.5, L, = 1.5 - j.5

54 3.6 Summary A procedure was developed for computing the scattering by 2-D structures. This procedure combined the boundary integral and finite element methods, and the resulting system was solved via CGFFT. The main advantage of the new approach was a reduction in memory demand to O(N) compared to O(N2) required with traditional solution techniques. A detailed map of the storage requirements was presented, and the principle memory consuming arrays were discussed. Also, the computational efficiency of the technique was examined and found favorable. To validate the proposed solution approach, several backscatter patterns were presented and compared with results based on traditional solution methods.

CHAPTER IV A Finite Element - Boundary Integral M1ethod for Two-dimensional Scattering Using Circular and Ogival Termination Boundaries 4.1 Introduction It is possible to choose other boundaries that result in convolutional integrals, and in this chapter we consider circular and ogival enclosures. Clearly, a circular enclosure would be attractive for circular scatterers whereas an ogival boundary will be more attractive for those structures conforming to this boundary. In the case of the circular boundary the entire integral is convolutional ensuring the O(N) memory demand of the system provided an iterative solver is used. When an ogival enclosure is used the integral becomes convolutional only if the observation and source points are on the same arc, but an efficient storage scheme is again required for the remaining "cross-terms" 1 In the following sections, the pertinent FE-BI formulations are developed for the circular and ogival boundaries. Results for several circular and ogival structures are presented and shown to be in excellent agreement with that obtained by traditional methods. 1 across terms" refer to integrals for which the source and observation points are not on the same arc

56 y x Figure 4.1: Geometry of the scatterer 4.2 Analysis Consider the plane wave ~fnc(_) =;qtinC() = ejkoPcos(Oo) (4.1) illuminating a composite cylinder as shown in Fig. 4.1 and we are interested in computing the scattered field. For the application of the Finite Element - Boundary Element Method the target is enclosed in a fictitious circular or ogival boundary as shown in Figs. 4.2 and 4.3. Within the boundary La, the finite element method is used to solve the Helmholtz equation V. [v(p)Vq(p)] + ko2v(p)q(p) = 0 (4.2) where +(p) = Ez(p), u(i)= ( v(j) = er(0) (4.3) Jiii~~i r(P)ii~~~iii

for E-polarization and O(p) = Hz(T), u(p) -= V(p) = Ur(J) (4.4) for H-polarization. The free-space wave number is k, = wafj7i and,r and Er are the relative permeablility and permittivity, respectively. On the boundary r, the Helmholtz integral equation O(~) = oinc(T)-f {G(7 pa) [$(Ta)] -(Pa) [OG(ia)] }dia (4.5) a ] c( [ 0G, )] d/ (4.5n provides the required boundary constraint, implicitly satisfying the radiation condition. In (4.5) G(,P a) - H(2) (ko I;- al) (4.6) is the 2-D free space Green's function where H42)(.) denotes the zeroth order Hankel function of the second kind. Also, a- denotes differentiation with respect to the outward normal, whereas p and Pa are the usual source and observation points, respectively, and IPPaI-x)2 + (Y - y a)2 (4.7) 4.2.1 Circular Enclosure Discretization of the Scatterer and Field Quantities The region enclosed by ra, denoted as Ra, is discretized into N, finite elements as illustrated in Fig. 4.2. In the figure, Pa is the radius of the circle and cxa is the integration angle along this boundary (Further definitions for the finite element mesh are indicated in Table 4.1, while the definitions of the field vectors are indicated in Table 4.2.). We note that nodes along ra are equispaced with angular displacement A.

58 -. ~ iiiiiiiira... ra t Figure 4.2: Partially discretized body in a circular enclosure Symbol Description Nn number of nodes in the finite element mesh Ng number of unknowns Ne number of elements in the finite element mesh Na number of nodes or elements on ra Nd number of nodes or elements on Fd Table 4.1: Definition of various quantities associated with the finite element mesh

59 Symbol Description qOa field at the nodes on ra /)a normal derivative of the field at the nodes on Fa qoI field at the nodes region I enclosed by ra and Fd qd field at the nodes on I'd Table 4.2: Definition of various field vectors associated with the finite element mesh and its boundary Derivation of the Finite Element Matrix The weighted residual expression over each element may be written as [5] JJReWiedSe = 0 i =1,2,3 (4.8) Se where Re = z _- u(x, y)) x('., - yy [(x, y)J y - V(x, y))-( (X, y) (4.9) In (4.9), Wie is the ith weighting function associated with the eth element. Substituting (4.9) into (4.8) and invoking the divergence theorem yields a~e awei a__ ae -J { x4 + a xi + k2vOe We} dQe. + j| Wiye'defe - 0 (4.10) where re denotes the contour enclosing the eth element. Additionally, e a U (4.11) is zero outside element e. Summing over Ne elements we obtain e=l- 1 [e aO e + e ax a a2yXe JO y Na Nd ~-zj WiS sd + zf/r' WSSdI 0= (4.12)

60 where the summations over s refer to the elements with sides adjacent to the fictitious (P,) and conducting (rd) boundaries. The integral over the conducting boundary in (4.12) vanishes all cases. This is obvious when no conductor is present. When X = Hz, the normal derivative of the field is zero on the conductor and the field unknowns on the boundary are allowed to "float" (i.e., they are not fixed to a specific value but are determined by solving the system). Finally, when q = E,, imposing the Dirichlet condition after assembling the finite element system results in the elimination of those equations associated with the integral over rd. Proceeding with the discretization, the field and its derivative within each element may be expanded into a linear combination of shape functions 3 =O E Nje (4.13) 3 f'E Ne ~] (4.14) k=l Substituting (4.13) and (4.14) into (4.10) and choosing Wie = Nie (Galerkin's method), we obtain 3 3 3 ZE ajq- E b!(k)Oj/e = 0 (4.15) j=1 k=1 j=1 where fra= u[a id +.. i I k2vN. T de (4.16) ex ax d ay ay and bY(k) = Nie die (4.17) For linear triangular elements, N?" are given by Ne = 2e(ae + be + cy) (4.18) r 2f " v

61 with Qe = -(b ec - bjce) (4.19) ae= Xk -k Xkyj (4.20) be = yej- _k (4.21) ci = Xek - x (4.22) and (xz, ye) being the coordinates of the ith node of the eth element. From (4.18) =N_ be (4.23) ax 2f~e OyN= _ (4.24) Using (4.23), (4.24) and the identity /(Nl)p(2)dxdy = Qe p q! (4.25) (N)P(N)dxdy 2(p + q + 2)! aij in (4.16) reduces to e ue Ue a3j - 4Qe (bb; + crcT ) - k2ve 12 (1 + &j) (4.26) where 1 if ij =ij = (4.27) O otherwise We note that in deriving (4.26) we have assumed that u and v (the reciprocal of the material constitutive parameters) are constant within each element and are given by ue and ve, respectively. To find an algebraic expression for b, we may reparametize the integral in (4.17) as bik = Pi8Psrada (4.28) 1

62 where Pi' and P2 are given by P )- = a P1(8) - 1-,, -,, (4.29) P;((a) = a a (4.30) Integrating, we have bk r 6 (6ik + 1) (4.31) Substituting the previous equations into (4.12) a sparse matrix is obtained for the nodal fields that has the form Aaa Aai 0 -Baa Oa 0 AIa Ali Aid 0 | I 0 (4.32) 0 AdI Add 0 kd 0 0 0 0 0 V-'a 0 In this, the values of the elements in the submatrix Apq are the contributions associated with the nodes in group (region or boundary) p which are connected directly to the nodes in group q. Also, [Baa]ik = b s k = 6 (6i-l,k + 46ik + 5i+1,k) (4.33) 8=1 The last row in (4.32) has been intentionally left blank to imply a need for another set of equations relating the fields and its derivatives on Fa. This additional set of equations is produced by discretizing the boundary integral equation. Evaluation of the Boundary Integral The boundary integral in (4.5) may be rewritten in cylindrical coordinates via the transformations IP - i: = (p cos a - p cos aa) - y(p sin a - pa sin a.)| _ p2 + pa2pp cos(a-a) (434)

63 where (p, a) and (pa, a) are the usual source and observation points in cylindrical coordinates. For Ip = IPal, IT- Pa = 2p Isin (a-,)1 (4.35) and the Green's function and its normal derivative may be written as G(p, Pa) =- H(2) (2kopa sin (~ —a)) (4.36) 40 2 G( jk H2) (2kopa sin (~-Ha))sin aa) (4.37) dfa 4 2 We may now write (4.5) as ~1 (p, a) = qinC(p, a) - fo(p, a) + fi(p, a) (4.38) where as a result of (4.36) and (4.37) fo(p, a) = - ap, a (p,,c a)H2) (2kop, sin (~a-)) daa (4.39) fi (p, a) = a f(pa, aa)Hf2) (2kopa sin ( ~- ))sin (a -n) da, (4.40) 4 02o with (pa, aa) = I a)(paaa) (4.41) The factor of 2 in (4.38) accounts for the singularity associated with H (2)() and the f (4.40) denotes principal value. We may now discretize (4.39) by expanding the field using pulse basis functions as Na p(p, apa) _ Z PA (aa - aj) pj (4.42) j=1 where )a(a j O otherwise a)=-(4.43) 0 otherwise

64 and A is the angular width of the integration cell as indicated in Fig. 4.2. Thus, the discrete version of (4.39) may be written as fo(pa) =-, jb J ) H(2) (2ko pa sin (- aa)) dae (4.44) j=1 J3W Performing point collocation and letting u' = a - a,, we have fo(p, ai) - jp A N [('i-3)+2 H2(2) (2kopa sin ()) du' (4.45) 4.5 which may be written in compact form as fo(p, ai) = jPa, /jho(ai- aj) (4.46) j=1 where ho(ai - j) = ( ) H(2) (2k0oPa sin (u)) du (4.47) (-)i-cJ)- 9 It is clear that (4.46) is in the form of a discrete convolution and can thus be written as fo(p, a) = DFT-' {DFT(b) * DFT(ho)} (4.48) where the elements of ho are given by a| {1-j2 [ln (ykop a) 1]} p=0 ho(pA) = (4.49) (P+)a2 H (2 (2kopasin (U)) du' p =l1.Na -1 (p-)A 2 where in developing the p = 0 term, the small argument approximation of the Hankel function was used and is given by [12] J -j;1n(.P) n-0 lim H(2)(kp)=n 2 (4.50) ~p *0 1 (kp)+Jn) n > (1)

65 in which -y ~ 1.781. Through a similar analysis, the field may be approximated by the expansion Na ~)(pa, Qa) _ E PA(aa - aj)Oqj (4.51) j=1 and by substituting this into (4.40), we obtain fi(p, ai) -= j4P jhl(ai-j) (4.52) j=1 where /(a_"j)+ " hi (ai - aj ( = J, (2(kopa si s ( ) du' (4.53) J (hi- )j-) 2 Clearly, (4.52) may again be written in operator form as f,(p, a) = DFT-1 {DFT(q) * DFT(hl)} (4.54) where kopa (A- sin A)+j hI(pzA) = (4.55) I r(P+') H?2H) (2kopa sin (u))sin (u)du' p=l,..., Na 1 where again (4.50) was employed. Point matching (4.38) at each node results in the system =21i = rinc - fo(p, ai) + fi(, ail) (4.56) which may be written in operator form as M..4a - La. O,. - 4inc (4.57) Maaqa Laa2I)a = a(457) where [Laa]ij _ jpa ho(ai - aj) (4.58) [Maa]ij = I- h(-a) (4.59) -~ 2 ~1 4

66 aa2 Paa 2 Pal a Figure 4.3: Partially discretized body with an ogival enclosure A final system is obtained by combining (4.57) with (4.32) to yield Aaa AaI 0 -Baa Ota 0 Ala AII AId 0 I 0 (4.60) 0 AdI Add 0 ~ d 0 Maa 0 0 -Laa Ia LOaC which can be solved via the conjugate gradient algorithm to obtain the nodal fields. 4.2.2 Ogival Enclosure Discretization of the Scatterer and Field Quantities The region within ra, denoted Ra, is discretized into Ne finite elements and a partial discretization is shown in Fig. 4.3 for the circular case. With respect to deA final system ios obtahe finite element mesh are ining (4.57) wi th (4.32) to yield vector definitions are indicated in Table 4.5.

Symbol Description Ap angular displacement between nodes on Fap Pap radii of Fap aap angular integration variable along Fap t distance between centers of curvature of rap ycp y-coordinate of the center of curvature of rap Table 4.3: Geometric quantities in reference to the figure above Symbol Description Nn number of nodes in the finite element mesh Ng number of unknowns Ne number of elements in the finite element mesh Na number of nodes (=Nal + Na2) on Fa Table 4.4: Definitions for the finite element mesh with an ogival enclosure Symbol Description'kap fields corresponding to the nodes on rap, p = 1,2 0'ap fields corresponding to the midpoints of the nodes on rap I %p+~ fields at the nodal midpoints on rap XI fields corresponding to region I enclosed by Fa and Fd'd fields corresponding to the nodes on the rd Table 4.5: Definition of the field vectors for the ogival boundary

Derivation of the Finite Element Matrix The derivation of the finite element matrix follows that described in section 4.2.1 with the exception of the matrix B,,aa. Consider the ogival boundary as indicated in Fig. 4.3. The boundary contour r, is comprised of two arcs labeled Fal and Fa2, which form the vertices of the ogive where they meet. At the vertices the unknown normal field is discontinuous and will therefore be evaluated at the midpoint. Also, in evaluating the contour integral, the field derivative are expanded in terms of pulse basis functions, rather than linear functions. This results in a different Baa matrix and involves the replacement of Pee in (4.28) by the pulse basis function expansion 1 if 0 < la - aj < P~a (a - a- ) =t-i< 2 (4.61) 0 otherwise By integrating in cylindrical coordinates we then obtain b= -(6+ i,j+l), j = 1, i= 1,2 (4.62) where 1e is the length of the eth boundary element along Fa and is equal to Pa,,pap for rLp, p = 1,2. Performing a summation over all boundary elements then yields Ne l j [Baa]ij = bi 2 (6ij + i,j+1) (4.63) e=1 where ii is the length of the jth element since the jth "node" (associated with the unknown $bj) is at the center of the jth boundary element. The remainder of finite element analysis for this case proceeds exactly as in section 4.2.1. Evaluation of the Boundary Integral The evaluation of the boundary integral along an ogival contour is similar to that described for the circular boundary. For integration and observation points

69 on the same arc, the integrals become convolutions. On the other hand, when the integration and observation points reside on different contours, the integrals have no special form and must be discretized and stored in memory as efficiently as possible. The distance between the source and observation points in terms of cylindrical coordinates for points on the same arc is given by P- PapI = /p2 + p p -2pp cos(a - crp) p = 1, 2 (4.64) When the source and observation points are along different arcs, (4.64) becomes IPq - Pap = /(p ccos aap)2 + (p sin aq - pa sinaa + YCq - YCp) p, q =- 1, 2 (4.65) in which the subscript ap refers to the integration coordinates along contour p and the subscript q refers to the observation coordinates. Also, yp is the y-coordinate of the center of curvature for contour p for p = 1,2. For further reference we note that (4.65) may be also rewritten as 2 1 (PI + P - 2p1 pa2 cos(ca - a2 ) + t2: 2t(pi sin ct - pa2 sin a2 ) (4.66) V 2 1 2 1 2 2 2 1 1 where t = yc2 - c, and the upper sign corresponds to the upper set of subscripts. To discretize (4.5), the fields are expanded as Nal Na q(pa, aa) _ E Pa((aa - aj)6j+~ + >j Pa(Qa - 9)j+~i (4.67) j=1 j=Nal +1 Nal Na #(p, aa) E Pa(aa - aji)/j + E PA((a - aj)Oj (4.68) j=l j=Nal +1 where as before Px(e~ - c |) = - tajil (4.69) 1 0 otherwise

70 and &j+I = ](Sj + qj+') (4.70) Substituting (4.67), (4.68) and (4.69) into (4.5) then yields! q(p (,a) = i"'(pl,al) E- Z4 fJ Go(pl, pal,I - al,)paldal j=l crj + / j+a f O Go(pl, Pala - aal)Paldaal Na aCt3+A2 - ctj + A Go(pl,Pa2 Ol la2 )pa2da2 j=Na, +1 -P-Go0(p1,Pa2,I,ola2 )pa2daa2 (4.71) j=Na, +1 apa2 when the observation point is on ral and,(P2, a2) = oinc(p2, a2) Nal aj+A1 + fj + / 0Go (2, Pa Go (P2, aal)Pa) draal, j=Nal +1 /xi + E +.j a Go(p2, Pa2a, a2 - aa2)Pa2dCaa2 (4.72) j=Nal +1 aj pa for observation on rF2. Performing point collocation at the nodal midpoints, (4.71) and (4.72) further reduce to: (Pl, ai+~ (m, ai+~) -- b E,\oi-otj q Go(PlPal, U)Padu j=l 2

71 Na ca j +A2 /~,-' Ij+ | Go(pl, pa2, i+u, aa2)Pa2dala2 j=Nal +1 l Na rj+A2 + Z j+ f a Go(pl, Pa2, ti+, aa2)pa2daa2 (4.73) j=Na1+1 + 1Pa for observation on Fa1 and -1 (p2, a+L) -= iC(p2,Q i+) Nal.aj+Al - ~>+ fX' Go(p2 pal ai+, aa)Paldaal Nal a+i.+Al + Z Oj+1 Go(P2, pal ai+ 1, aa,)pal daal Na ai-aj+ - - — Go (P, Pa2, a U)pa2 d (474) j=Na,+1 c1 - 2,+7 + ig djt 2 i~.i4, a +&aGo(P2,Pa,, U)Pa2 dL (4.74) where the'' in the subscript refers to the fictitious "node" midway between the actual nodes. A system of equations can now be obtained by testing (4.73) and (4.74) at a sequence of points on the contours. This yields -CI ial =. - + {ll4al + P12a2,- 1MI Cl al - Q1 2C2a2 _C2(a2 _ I- {P21 Oal + L22ia2- Q21C al-C1 M22C2a2 } (4.75) 4C q~a2 = a2+~ 2 P212al + L224'a2 which can be alternatively written as aal M 0 Q1202 D LD P12 | |a2 |. LA+ (4.76) 21 C1 M/22C2 IP21 L22 4ai -.JL 4 a:+~Q

72 In this, the nonzero values of the upper bidiagonal matrix C are 1, and $incl the value of the incident field evaluated at the nodal midpoints. The matrix D accounts for the double use of the nodes at the endpoints and the remaining elements are given by JMPP -2 I( - MPP) (4.77) Qpq - Qpq (4.78) (4.79) for p = 1,2 in which [MPP]ij 2= Go(pp, Pap u)papdU (4.80) [Lpp]ij = f/ + Go(ppPapU)PapdU (4.81) [Qpq]ij = ] Go(pp, paq, a i+~ cIaq)paq,daaq (4.82) and [Ppq]ii =G J p —Go(pp, pa, i+, Caq)paqdcaq (4.83) j, j =Paq 2 More explicitly, upon evaluation of the integrals f(' koa|, j [koPap( - sin 2 ) + jap [Mpplij = a (4.84) 1~* i'-i 4'+' H )(2kopap sin U!)sin'du i y j 4 fce j +'&4 2-11} [Lpp]ij = (4.85) -j'Pa - f + H(2) (2kopap sin )du i j [] =~ pa, S~i ~C-j 24 [Q ] jkoPa2 |j+A2 H )(ko/:: ) IQk121 [Pa-2 P cos(aa- Qa 2 )+t sincra2] dsa2 (4.86) 1Pj] = i-: jf'+| HnO2)(ko.T)d.a2 (4.87) -jq, a

where the upper sign corresponds to the upper set of subscripts on P or Q, while the lower sign corresponds to the lower set of subscripts. Introducing the definitions 11QC1 Q12C2 K1- IfI| M 1~ C1 O2 ] t ](4.88) Q21C1 M22C2 Li, P12 IK2 = (4.89) P21 L22 the system (4.76) may be combined with that derived via the finite element method to obtain Aaa AaI O -Baa O 0 AIa Ail Aid 0 qI 0 (4.90) 0 AdI Add 0 5d 0 1(1 0 0 K2 Oa jai We note that (4.90) can be solved via the CG algorithm to take advantage of the convolution operators M and L in reducing the memory requirements. This algorithm is given next. 4.3 A CGFFT Implementation In a manner similar to the previous chapter, the matrix vector multiplications Az and Aaz is represented as a summation of matrices, one corresponding to the finite element portion of the system and the other corresponding to the boundary integral portion. Thus, a typical product may be represented as {s} = {S)BI + {S}FE (4.91)

74 where O O O O z 0 00 0 z2 {S}BI = (4.92) 0 0 0 0 Z4 K1 0 0 IfzJL4 and Aaa AaI 0 Baa Z1 AIa AII Aid 0 Z2 {S}FE = (4.93) 0 0 Add 0 Z3 L0 0 0 0 z4 For the adjoint operations, we have 0 0 0 K1 z1 O O O O z2 {S}BI = (4.94) 0 0 0 0 Z3 0 0 0 K2 z4 and Ai Aa: Aq AI, Aa I AId 0 z2 {S}FE = (4.95) 0 AdI Add 0 z3 Baa 0 0 0 Z4 Again, the operation is performed such that only the resulting vector {s} needs to be stored. 4.4 Scattered Field Computation In this section the expressions for the scattered field and radar cross-section are developed for both the circular and ogival boundaries. The scattered field is com

75 puted from the identity $() G(=j { a [ )] _ + (Ta) G(, } dia (4.96) from which the echowidth is calculated. 4.4.1 Circular Boundary The scattered field expression (4.96) may be written as q$(p, a) = -fo(p,) + fi (p, a) (4.97) where [21r H'(2 2 fo(p,) = Pa j o(Paa)H2) (kop + P- Pa cos(a - aa)) da (4.98) and fi(pv)=,Pakoj (p 4 2) (ko\ p2+Pp-pPacos(e - aa)) fl (p, a~) Pa= o.a pC ao 4 ~ ~/p2 + pa -2pp cos ( - a) [Pa - pcos(a - aa)] doa (4.99) To evaluate the integrals in (4.98) and (4.99) we invoke the field expansions (4.42) and (4.51). We have fo(p a) = I-pZ J' | H42) (ko~/p2 + Pa-2ppa cos(a —a))) dca (4.100) and j Na j+4 H(2) (ko/p2 + p - 2PPa cos(-a)) fi(p, ac) = kopa j=I aj2 + I22 _ ppa COS(-a) [pa - pcos( - Ca)] d a~ (4.101) where the remaining integrals over the subsections must be evaluated numerically for arbitrary observation. However, for far-field computations (p - oo), the Hankel

76 functions may be approximated as Hn(2)(kp) 2ne-jkp (4.102) and since + - cos(a - ) for amplitude (4.103) /p2 + p2 _ 2ppa cos(a- aa) _- P(4.103) p - Pa cos(a - aa) for phase terms (4.100) and (4.101) become Afo(p a)=2j e-jkop i ejkoP cos(a-_i) (4.104) 4 irk0 p 3=1 and Pa koA'2j-o~Na f1(p, a) = 2j -koP E Oj cos(a - aj)edkoP. cos(a-") (4.105) r4 7kop j=1 Substituting (4.104) and (4.105) into (4.97) we obtain,(p, a) =PaA J jko p [, N jko COS(-a (p a) 4 rk j eikopaco +klo E Oj cos(a - aj)ejkoPacos(-aij) (4.106) j=1 and from (2.31) the echowidth is given by yields the echowidth a)2 = N) Na 2 = (p'a/ 5 j Z eikopcos(c) + k0o Oj, cos(a - aj)ekopos(a-( ) (4.107) 4.4.2 Ogival Boundary Following the same discretization scheme used in Section 4.2.2, (4.96) may be written as N01 Na q'(p, a) = -{ (?/fil(p, a, aj) + ~ /f12(p,, aj) j=1 j=Na1 +1 N,1 Na - 7 cf21(p, a, ai) - Z qf22(P,aaj)} (4.108) j=l j=Nal +1

77 where fip(p,a,, c) = Go(p, pap, a, aap) Papdcap (4.109) f2p(P, a, i) = | ap Go(p, Pap aap)papdaap (4.110) (4.111) in which Go(p, Pap, a, rap) =- H2) (kop2+p 2PPap cos(a -ap) +y -2ycp(p sin a-Pap sin aap)) (4.112) and pap Go (p, Pap c, aa, p) jko Hi21 (ko/p + Pp- 2ppap cos(a- aap) + y2 _ 2yc(psin a- Pap sin aap)) 4 /p2 + p - 2ppap cos(a - aap) + y2 - 2yc(p sina - Pap sin aap) [pap - p cs(a - aap) + ycp sin aap] (4.113) and ycp are the corresponding y-coordinates of the arc Fap. Using the large argument approximation for the Hankel function and the approximation J/p + pa - 2ppap cos( - caap) + yc2 + 2ycp(psin - Pap sin caap) p for amplitude terms (4.114) P - pap cos(a - caap) + ycp sin a for phase terms for p - oo, the Hankel function simplifies to 2j e-jkope-jko[- apos(o- saap)-ycsin] (4.115) YCJ (4.1p

78 Similarly, 1| } (\/;;7i[ cp j kop e-jko [-Pap cos(a-Caap)-ycp sin a] COS(Q - aap) IY eC P cos~a - a"'p) (4.116) Substutiting these into (4.109) and (4.110) and performing midpoint integration yields LApPap 2j ejkop -jko[-Pap cos(ca-Oj-P)-ycpsina] (4 117) ~fpp 2. =ejkope-jko[-Papcos(aJt-c)-Ycpsina] f2p(p, a, aj) = k -j4rkp.pe r'kop Cos(aj + - aap) (4.118) Thus, from (4.108) 1 2/2-jkop q(p, a) - 4 rko p{i~p, ~ p.+~al e- jko~[-Pai COS(~(0'- -Al )-Ycl sinc] Na1 JAliPal Y1 Oj+ I ek2)-y sina] +Ja2Pa2 E 4j+i e-jko[ Pa2 COS(a- - A )/cY2 sina] i=Nal +1 1 k~p! E j+'' -jko[-pal co(a-aj )- sina COS( - ) + koLAPal E Oj+ e jko Pal cos(aa3 2) )-Y' sin a] aCO j)] j= 2 2 Na (4.119) the echowidth becomes a r|1 PalE t e-jk~[-pOl cos(a-a-i')-ycl sinca] + j2Pa2 1 4j+, -Jiko~[-Pa cos(a-aJ-.2-)-Yc2 sina] j=Nal +1 Na ~1 + ko pa 1Z O j+ e-jko[-Pa l cos((-aj)-yc- sn] cos(a, + - ) 4 —1

kNPa e2koKPa2cs(ca - sin] cos(a + 2) 3=N0 +12 (4.120) 4.5 Results The scattering patterns of several circular and ogival cylinders for both E- and H- polarization are shown in the figures to follow. Figs. 4.4-4.6 contain circular geometries both coated and uncoated, while Figs. 4.7-4.9 contain coated and uncoated ogival structures. The echowidth is computed for each structure and compared to the results of the series solution for the circular bodies and moment method [13, 14] for the ogival structures. As seen in all cases, the generated patterns via the FE-BI formulation are in excellent agreement with the corresponding data based on the Mie Series and Moment Method Solutions.

80 Perfectly Conducting Cylinder R=.5 X (E-pol) 20.0. g I i * X X 15.0 FE!{BE' 10.0 1 —----— ~......... Solufion -, 10.0 _e 5.0 0.0 LmL -5.O. -10.0 -15.0 -20.0..... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Perfectly Conducting Cylinder R=.5 X (H-pol) 20.0. I I I I I 15.0BE 00 —. —-- SeriesSolution -U 10.0.5 5.0 0o.o -5.0 -10.0 -15.0 -20.0. 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Figure 4.4: E, and Hz bistatic echowidth of a perfectly conducting circular cylinder of radius 0.5A

81 Coated Conducting Cylinder R=.5 X (E-pol) 30.0 20.0 m320.0. [..FE/BE _0 10.0 -—.. —-~ Sries Solution 10.0 o -10.0 -20.0 -30.0.... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Coated Conducting Cylinder R=.5 X (H-pol) 30.0. --—......... SriesSolution 0.0 -10.0 -20.0 -30.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Figure 4.5: E, and H, bistatic echowidth of a perfectly conducting circular cylinder with a conductor radius of.5A and a coating thickness of.05A containing material properties er = 5- j5, I, = 1.5- j0.5

82 Coated Conducting Cylinder R=3 x (E-pol) 30.0'. 25.0 ~.X 20.0 - PE/BE -......... Sris Solution 15.0 o.o 10.0 -5.0 -10.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Coated Conducting Cylinder R=3 A (H-pol) 30.0.. I I 25.0 M. 20.0 | FE/BE........ Series Solution 15.0 10.0 5.00.0 -5.0 -10.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Figure 4.6: Ez and H, bistatic echowidth of a coated circular cylinder with a conductor radius of 3A and coating thickness of 0.05A with material properties E = 5 - j5, /lr = 1.5 - jO.5

83.5 x 1.X PEC Ogive (E-pol) 20.0. I I.. 15.0 FE/BE 10.0 5.0 o.o MeO -5.0 -10.0 -15.0 -....,....,.... [.... X.... -20.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle [deg].5 x 1.f PEC Ogive (H-pol) 20.0.... I 15.0 FE/E!0.0 5.0 -10.0 -15.0 -20.0.... 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle [deg] Figure 4.7: E, and H, backscatter echowidth of a 0.5 x IA perfectly conducting ogive

84.5 x 1X Coated Conducting Ogive (E-pol) 20.0...... 15.0 mBE.. —----- MoM 10.0 5.0 0.0 -5.0 -10.0 -15.0 -20.0....................... 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle [deg].5 x 1X Coated Conducting Ogive (H-pol) 20.0... I I I I * 15.0 PEBE......... MoM 10.0 5.0 -t 0.0 -5.o.10.0 -15.0 -20.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle [deg] Figure 4.8: Eh and H, backscatter echowidth of a.5 x 1A perfectly conducting ogive with a 0.05A thick material coating containing the properties Er = 3j5, Ir = 1.5 - jO.5

85 1 x 4X Coated Conducting Ogive (E-pol) 20.0.... 10.0 0.0 o -10.0 -20.0 -30.0. -30.0......... i....,.... 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle [deg] 1 x 4X Coated Conducting Ogive (H-pol) 20.0... 10.0 - 0.0 o -10.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle [deg] Figure 4.9: E, and Hz backscatter echowidth of a 1 x 4A perfectly conducting ogive with a.05A thick material coating containing the properties c, = 3 - j5, /#r = 1.5 - jO.5

4.6 Summary The scattering from targets surrounded by ogival and circular boundaries has been presented. The finite element method produces the usual sparse sub-matrix, while a discrete version of the boundary integral results in a dense sub-matrix. The mathematical boundary enclosing the scattering structure may be judiciously chosen such that the boundary integrals are convolutional. As a result, they become amenable to evaluation via the FFT and leads to an O(N) storage requirement. Among the circular and ogival boundaries considered, the circular boundary satisfies the above requirements. The ogival boundary results in convolutions only when the source and observation points are along the same arc, while the non-convolutional cross-terms must be stored efficiently to guarantee the required storage requirement. When considering circular and ogival structures, the associated circular and ogival boundaries are usually conformal to the structure, thus providing an additional reduction in the number of unknowns. To validate the method and associated computer code, scattering patterns of several circular and ogival structures were given and compared with data generated by proven methods.

CHAPTER V The Elimination of Boundary Integral Interior Resonances in Two-dimensional Finite Element - Boundary Integral Formulations 5.1 Introduction The interior resonance corruption of boundary integral solutions for scattering computations is well known, and its treatment has been a subject of research for the last two decades. Methods based on the "complexification" of the wavenumber [15], the overspecification of the boundary conditions [16], [17], and the linear combinations of integral equations [18], [19] have been proposed, while others have focused on the solution technique rather than the system formulation [20]. Not surprisingly, when the boundary integral equation is used to terminate the finite element mesh, the interior resonance corruption persists and this has restricted the application of an otherwise promising method for large scale computations of highly inhomogeneous structures. To demonstrate the seriousness of the problem, Fig. 5.1a shows the backscatter echo width of a circular conducting cylinder of radius ao for TM plane wave incidence. The mesh is terminated on a circle of radius a0 = 1.01a, on which the boundary integral equation is applied. The results, displayed as a function of ka0, were obtained

via the FE-BI described in chapter IV. The unknown in this implementation is the total field, and as seen the solution is severely corrupted, a difficulty which persists for other scatterers as well. (For reference, the locations of the interior resonant frequencies are displayed by the vertical lines at the bottom of the figure.) The corruption is further evidenced in the near zone scattered field plotted in Fig. 5.lb, obtained by subtracting the incident from the total field generated via the aforementioned FE-BI method at kao = 23.586. To render the FE-BI method robust at all frequencies, it is thus essential to remove the problem associated with the interior resonances. However, employing traditional methods such as those described in [2] - [6] requires either significant modifications to the original FE-BI formulation or substantial computing time, thus affecting the efficiency and performance of the method. Further, for the "complexification" scheme proposed in [15] to be effective, we verified that three different computations are required for each incident angle when combined with the total field FE-BI system in (4.60). That is, the total field FE-BI solution must be repeated three times with different complex propagation constants, all slightly perturbed from the free space wavenumber. Fig. 5.2, the counterpart to Fig. 5.1, demonstrates how the amplitude of the field quantity varies for ca = - jO.001 and ca = 1-jO.005. In the following section we present a simple modification to the FE-BI formulation which renders it relatively immune to the resonance problem without the need to repeat the solution. 5.2 Method The proposed approach for eliminating the failure of the FE-BI method at interior resonant frequencies is based on the observation that the specification of the

89 tangential electric or magnetic field excitation alone at the integration boundary is not sufficient to yield a unique solution [19]. Therefore, we move the excitation away from the integration boundary by employing the scattered field as the working variable. This is accomplished by first writing the total fields everywhere in space as a superposition of the incident field Xi and the scattered field O$s as = qY + Os (5.1) The boundary integral in (4.5) may be expressed entirely in terms of the scattered field as Os(7)=_ if {G(pPa) [a Sn a (Pa (Ta) [ aa G(;71 }a)dla (5.2) which differs in form from (4.5) by the incident field term. The excitation is instead associated with the finite element portion of the system on conduction surfaces and material interfaces. This becomes clear after substituting (5.1) into (4.10) to obtain II {_ [8+? 9~ie + Hi a~j ] + ko i d + Ire i dr -u ax aax x axy + yy 5 3vqS'W df+.d -- - U axx ax y +- k~vqbiWf dQ (5.3) where the quantity x, has been left in terms of the total field to ensure tangential field continuity between adjacent elements which may contain different materials. The expressions for TM (q = Ef) and TE (Os = Hz) differ by the application of the boundary condition on perfectly conducting surfaces. For TE incidence, b = 0 and the contour integral contribution on this portion of the path rF disappears. For TM incidence, the condition Os = -qY is applied after assembly, and this results in the elimination of the associated contour integral. Thus, in following the discretization procedure outlined in section 4.2.1, we obtain the system 3 3 3 a,{}~ - _ E~ be,(k){/,8} e -Ue (5.4) j=l k=l j=l

90 where e { [( o~ + + rO aN + k~v2 ei Note that a.. and bVj are unchanged. Assembling the final system and applying the appropriate boundary conditions, we have Aaa AaI -Baa 0 -Ua AIa Ail AId 0 |I -UI (5.6) o AdI Add 0 ~ -Ud Maa 0 0 -Laa a O for TE and Aaa AaI -Baa 0 Ua AIa Aii Aid 0 |1| -U57 (5.7) 0 0 I o0 I d Maa 0 0 -Laa a 0 for the TM case. Clearly, now, the excitation is present entirely in the FE region of the system, as opposed to the BI portion as was the case in the total field formulation of chapter IV. The system is then solved with the introduction of a small loss in the propagation constant appearing in the Green's function in (5.2) accomplished by replacing k with ak, where aC = 1 - j6. This substantially improved the convergence of the employed biconjugate gradient solver for the cases considered. However, in contrast to the "complexification" scheme employed with the total field formulation, the scattered field solution in (5.6) and (5.7) is relatively insensitive to 6 (provided, of course, S is very small).

5.3 Results To demonstrate the effectiveness, efficiency and accuracy of the proposed method, let us reconsider the problem of scattering by a circular cylinder via the scattered field FE-BI formulation. As seen in Fig. 5.3a, the far field is no longer corrupted by the fictitious interior resonances and the same holds for the near zone field as displayed in Fig. 5.3b. The results, shown for a = 1-jO.OO1 and a = 1 -jO.005, are also seen not to deviate from the series solution, although the convergence rate of the solver varied significantly as a function of a. For example, the unperturbed (no complexification, i.e., aC = 1) scattered field formulation converged in approximately 0.15N iterations over the frequency band considered in Fig. 5.3a, where N denotes the unknown count. For a = 1 - jO.001, the solution converged in 0.13N iterations whereas for ca = 1 - jO.005, convergence was achieved in approximately O.08N iterations. Note that this fast convergence rate is due in part to the fact that the discrete boundary integral occupies a substantial portion of the total FE-BI system. For complex structures, this may not be the case and the affect is expected to be less pronounced. Also, we consider the TM illumination of a 5.256A x 5.256A metallic square cylinder. For the implementation of the scattered field FE-BI formulation, the boundary integral was enforced on a circle of radius ao = 3.754A so that kao is the fifth zero of the Bessel function of order 6. With the incident field direction normal to one of rectangle's faces, Fig. 5.4 depicts the corresponding bistatic scattered field obtained via the total and scattered field FE-BI formulation with ca = 1 - jO.005. Clearly, the pattern based on the scattered field FE-BI formulation agrees everywhere with the moment method data. In contrast, the results based on the total field FE-BI formulation are substantially in error.

92 20.0... 15.0 o.0 5.0 - FE-BI......... Series Solution 0.0 11111111 18.0 24.0 30.0 36.0 42.0 48.0 54.0 60.0 66.0 kao (a) 2.0 00 1.5 0.5 0.0 0.0 60.0 120.0 180.0 240.0 300.0 360.0 Angle [Deg] (b) Figure 5.1: Comparison of the far zone and near zone fields for TM plane wave incidence on a circular metallic cylinder as computed by the total field FE-BI method and the eigenfunction series. (a) backscatter echo width vs. kao - the lines over the horizontal axis correspond to the eigenvalues of a circular conducting waveguide. (b) magnitude of the TM scattered field on the enclosure at the resonant frequency kao = 23.586

93 20.0 15.0. —---------- 10.0 a=1-jO.001 5.0 5.0 -.. a=1-j0.005 ------ Series Solution 0.0,.,, 18.0 21.0 24.0 27.0 30.0 33.0 kao (a) 1.5.....,...............,....., 1.0 co 0.5 -a=1-j0.001 -------- a=1-j0.005 ------ Series Solution 0.0......... 0.0 60.0 120.0 180.0 240.0 300.0 360.0 Angle [Deg] (b) Figure 5.2: Far and near zone fields for TM plane wave incidence on a circular metallic cylinder as computed by the total field FE-BI method with k replaced by kca in the BI equation. (a) backscatter echo width vs. kao. (b) magnitude of the TM scattered field on the enclosure at the resonant frequency kao = 23.586

94 20.0-,....., 15.0 *' 10.0 t 5.0 - a=1 -jO.001 5.0 50 —- =1-jO.005 ------ Series Solution 0.0 18.0 24.0 30.0 36.0 42.0 48.0 54.0 60.0 66.0 kao (a) r 1.0 0.0 0.0 60.0 120.0 180.0 240.0 300.0 360.0 Angle [Deg] (b) Figure 5.3: Far and near zone fields for TM plane wave incidence on a circular metallic cylinder as computed by the scattered field FE-BI method with k replaced by kac in the BI equation. (a) backscatter echo width vs. kao. (b) magnitude of the TM scattered field on the enclosure at the resonant frequency kao = 23.586

95 30.0 Scattered field FE-BI (a=l-jO.005) 25.0 m-" 2.0 -- - Total field FE-BI 6- 20.0 15.0 10.0 10.0 -- -0 5.256X -!0.0...... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle [deg] Figure 5.4: Bistatic echo width for TM plane wave incidence on a 5.256A x 5.256A metallic square cylinder with the boundary integral place at a radius a. = 3.754A ri~ 111 a a a.~~~~~~~~~~~~~~~~~~~~~~a a a aa aaa aaaa a

5.4 Summary A method was presented for the elimination of the interior resonance corruption of the boundary integral in the FE-BI formulation. By expressing the system in terms of the scattered field and employing a complex wavenumber in the boundary integral, the effect of resonances were removed. This implementation of the method was shown to be superior to that associated with the total field formulation presented in chapter IV in that only one sample per frequency was needed. Though the development was implemented for the two-dimensional FE-BI formulation, it is applicable to the three-dimensional one as well, as seen in chapter VI.

CHAPTER VI A Finite Element - Boundary Integral Formulation for Axially Symmetric Structures 6.1 Introduction A finite element - boundary integral approach is applied to the case of axially symmetric structures. The method follows the same procedure outlined in section 2.2. In this implementation, the coupled potential equations [21] are discretized via the usual finite element method, and the resulting system is augmented by a discrete form of the Stratton-Chu equations [22]. The storage reduction associated with the boundary integral is achieved by exploiting matrix symmetries and the final system is computed by employing a conjugate gradient solver. In this chapter, the formulation for the FE-BI is described for axially symmetric scatterers. The results presented demonstrates the accuracy of the method along with showing its limitations. 6.2 Finite Element Formulation In this section, we derive the analytical coupled azimuth potential (CAP) equations [23] which are then discretized via the finite element method. A consequence of the formulation is a line singularity, which tends to corrupt the computed fields

98 erXialys trcsfewth y general axially symmetric surface with source (prime vation (unprimed) points and the corresponding unit vectors for scattering domains incorporating lossless media. A regularization approach is suggested for its removal. 6.2.1 Analytic CAP Formulation Maxwell's equations in a source free region are given by [12] V x E() = -jwp#H (6.1) V x H(F) = jweE (6.2) V *-D(F) = 0 (6.3) X7 B() -= 0 (6.4)

99 For axially symmetric media such as that indicated in Fig. 6.1, the fields may be represented by a Fourier series in the cylindrical coordinate q as 00 E(r)= Z m(p,z)edm( (6.5) m= -oo oo rH(r) = E hm(p,z)ejm, (6.6) m=- o0 where m (Pz) = j1 E(2)e-Jm dq (6.7) hm(p, z)= 2- | t(=)e-Jd (6.8) Up substituting (6.5) and (6.6) into (6.1) and (6.2), we obtain [jmhmz - (Rhm)] = jEremp (6.9) [jmemZ (Rem,)] = -jrhmp (6.10) R [zhmp- a hmz] = jer(Rems) (6.11) R [h emp - eemz] = -jlr(Rhmo) (6.12) [jmhmp- aR (Rhm-)] = -j Eemz (6.13) R [jmemp - (Reom)] = jurhmz (6.14) with R-= kop, Z = koz (6.15) to be referred to as normalized coordinates. Substituting hmz of (6.14) into (6.9) gives emp = jfm [ma 8(Rem6) + RIr a(Rhmk)] (6.16) where f= [R272-m2]-1 (6.17)

100 and 2K 2- rr (6.18) Analogously, substituting emp of (6.9) into (6.14) yields hmz = jfm [ma (Rhmq) + Rer a(Rem6)] (6.19) while a similar procedure for combining (6.10) and (6.13) yields the pair emz = jfm [maa (Remk) - Rra -(Rhme)] (6.20) hmp = jfm [m _..(Rhm) - Rcr (Rem,)] (6.21) Equations (6.16) through (6.21) may be expressed in compact form as q x em(R, Z) = jfm [mq x Vt,e - prRVth] (6.22) x h-m(R, Z) = jfm [mq x Vtbh + ErRVtbe] (6.23) q. em(R, Z) =,eR (6.24) *hm(R,Z) = bh/R (6.25) where Vt = P8Ri + z- (6.26) and Obe and Oh are herein referred to as the azimuthal potentials. Rewriting (6.11) and (6.12) as RV, (q x hm) = -jer/e (6.27) RVt. (~ x m) - jurfIh (6.28) and then substituting (6.22) and (6.23) into them produces the CAP equations Vt* [fm(erRVtfpe + m Xvh)] + E4R = 0 (6.29) Vt' [fm(,RVt'h-m x Vtb)] + /r'/h = 0. (6.30) Li - gm) R

101 This system may be written in operator form as Lo = 0 (6.31) where Vt' [fmerRVt] + M mvt [frm x Vt] L= (6.32) [-mVt [frR X Vt] Vt [fmtirRVt]+ rJ and = [Oe Ohh]T (6.33) The three dimensional axially symmetric problem has, thus, been cast into a two-dimensional one and its discrete representation is formulated in the following section. 6.2.2 Discretization of the CAP Equations To discretize (6.31), we first enclose the generating contours of BOR in a fictitious boundary r (= ra + r + r,) as shown in Fig. 6.2. The region interior to r is divided into Ne triangular elements and within each element the corresponding weighted residual expression is JfNfe(R,Z) L dS = =O (6.34) Se where Nie is the usual shape function [5], so chosen to satisfy the Dirichlet boundary condition of 0,e and /h on r,. Substituting the first of (6.31) into (6.34) gives JJNi Vt,. [fm(ErRVt + mO x Vth] + } dS = 0 (6.35) and upon employing the identity NeVt. As = Vt (NeAS) - As. VtNe (6.36)

102 R n Z t.....e.......tp PEC ft [v.{'f ERtb q b)rz}... rz Furthermore, by invoking the divergence theorem, (6.37) may be written If f [-fm (ErRVtke + mq x Vt?,bh)]. VtN + ER }dj gvs - th RV~e e.R xpressj IiR i:] on(X7 + [n [Niefm (erRVt7e + m x Vth)] di = O (6.38) where ni is the outward normal along the boundary Ce of the eth element. Finally, on using (6.22) these may be simplified, yielding Ji|| { [i-fm (ErRVte + m x Vth)] VN~iR + rIN'. } NdS- Nri(jhmt)dl=0 ~(6.39) Furthermore, by invoking the divergence theorem, (6.37) may be writt........ $, /~~~~~~............ on using (6.22) these may be simplified........................ $, ~ ~ PE -~~~~~~~~-c -----------------------------— 0 ][HZ z~~~~~~~~~~~~~~(.9

103 where hit = t. hm (6.40) with t=i x (6.41) Equation (6.39) and its dual constitute a weak form of (6.31) over the eth element. The development thus far has employed the total potential as the working variable, but may also be expressed in terms of the scattered potential as well. To this end, any function i$ satisfying Maxwell's equations may be written as a superposition of the incident and scattered potentials, i.e., be = Oe + pe (6.42) Obh = O, + 0bh (6.43) where the superscript s denotes the scattered potential and i denotes the incident potential (i.e., that potential present in the absence of the scatter). Substituting these into (6.39), the corresponding expression in terms of the scattered potential is J/ { [-fmr (erRV'tt: + VN + N x } de -Ice Ne(jhmt)dl = JJ/ { [m (rRVt' + m VtN + eN } dSe (6.44) where the contour integral has been left in terms of the total field until the final system assembly (of equations) is performed. To form a discrete system of equations, the potentials in the e element are expressed as ik(R, Z) =' N (R, Z) {i,}e (6.45) j=l

104 bh(R, Z) = Nj (R, Z) {jb+} (6.46) j=1 and upon substituting these into (6.44) gives L' { [.~-fmerRVtNie VtNje + ENr (b] ) - ICe Nf(jhmt)dle = {u} (6.47) (6.47) where {U} is given below. Assuming Er and Pr are constant within the element, (6.47) may be written as 3 Z [fe [a],j {;e}j- [b]j {-j] = js Ni (jhmt)dle + {U}) (6.48) where [aI~j = I -Ifm RVtNte.VtNlN + ] dS (6.49) [b]iei = Jmfm$ x VtNe VtNe dS' (6.50) Se {U} = -JJ {[-fm (r R Vt/ \+ m Xvt)]. VtNi+ ER/N }dSe(6 51) These constitute the equations for the eth element, and the final system is assembled by summing over all elements in the discrete model. Before pursuing the step of system assembly, explicit expressions for (6.49) and (6.50) may be developed by first writing them as [a~I1 =+]JN jdSe (6.52) t3 (2Qe)2 s R S* and [b] - = m ( ) lo (6.53)

105 where the linear shape function and its coefficients are given by ayn, z) "? +z +I ve Nje(RI Z) = te + 3Z + R (6.54) 2Me e =- Zi+lRi%+2-Z i+2Ri1 (6.55) 3 - = Re+1 -RR2 (6.56) = Z,+2- Zl1 (6.57) and I1 R -dS (6.58) 10 fJ 1 (6.59) Se;2R2 m2d Evaluating (6.58) and (6.59) over a triangular element yields.1 ='-'[Ie (m) + Ti (-m)] (6.60) o = ( -m) - (-m)] (6.61) with 2~?(m) = (KR? + m) log(KRt + m) (6.62) /3,e+l f+2 The second integral in [a]ij may be computed numerically via gaussian integration. The sums appearing in Io and I,, however, experience problems when an element edge is parallel or nearly parallel to the axis of revolution (i.e., element d in Fig. 6.3). In this case, two of the three terms have denominators which are nearly zero, resulting in cancellation errors upon their summation. To avoid such errors, the sum of two consecutive terms, i.e., 2k + I+, is carried out by first expanding the numerator of either term Z and then performing the addition analytically. In this way the offending term is canceled and the final result is well behaved.

106 R m \a / /b R C k Z Figure 6.3: Typical triangular elements in the vicinity of the line singularity at Ro = m for real /c Proceeding then, we note that the Taylor series expansion of x Iln x about ox0 is xInx = xlnxo + (x-xo) [1 + SX Xo)] (6.63) where S(y)- ( —1) +lyn (6.64) n=1 n(n + 1) and the sum converges when Ix - xol < R, where 1? is the radius of a circle centered at x0 in the complex x-plane for which x ln x is analytic [24]. The series may be truncated in N terms when the error, given by the ratio of the Nth term to the first term in the series is less than some tolerance c, or more explicitly when 2yN-1 N(N~+1) <E (6.65) For e = 10-7, a fit to the nonlinear function in (6.65) for IYl E [0.1, 0.7] is N(y) = Round(eIyl(9.97633elyl - 5.94387)) (6.66) Once N is known, the series is evaluated efficiently using Horner's rule. Employing the expansion (6.63) to the partial sum Ii(m) + Ii+ (m) yields [ [(^Rs+2+m)ln(KR,+l+m) 1-S(K K/+2 )] form1 zi(m) + Ii+l(m) = +1 [(, P Ri+,+m'(6.67) l1 r[ (Ri+2+m)ln(KRi+ + 1 + S \I] form 2 Kf3~~~~ LI ~1~a+,+

107 where form 1 was obtained by expanding the numerator of ZI(m) about Ri+1, while form 2 results by expanding the numerator of Zi+l(m) about Ri. For elements well away from the line singularity, either form is valid. However, as KcR + m - 0, the radius of convergence 1R is reduced to zero due the branch point of (KR+ m)ln(KR+ m). The method for choosing the appropriate form may be done numerically by computing the value argument of S in (6.67) for each form. The one which gives the smallest value y E [0.1,0.7] in S(y) is chosen for computation, since the series will converge the fastest for it. If y falls in the range 0.7 < y < 1, the series (6.64) will still converge, but very slowly as y approaches 1. If IyI > 1, the series will not converge and the associated form (1 or 2, or both) will be invalid. Note that neither of them is valid if both points lie along the line singularity (in which an infinite result is obtained) or if one lies above it and the other below (as in element b of Fig. 6.3). In the latter case, (6.60) and (6.61) must be evaluated directly, but a new mesh may be necessary for the former. An approach for bypassing this difficulty is discussed in the following section. Summing over all elements to obtain a solution for the entire computational domain Q, our system becomes Ne 3.Ne Ne E E r adj {e I j {a Sh} - N. (i =t)dl+ I N (jhmt)dl e=l j=l e=l a e=l z ~~Ne Ne + E Ne (jhmt)dl + j{(U} (6.68) e=1 c e —l Note that since Ne and hit are individually continuous between adjacent elements (and the same is true for their product), the contour integrals cancel everywhere except on the domain boundary r. Accounting for Ne = 0 along the z-axis, this

108 system combined with its dual may be written in block matrix form as ACa 0 A", -Baa 0 -B 1i {/e}a {U}a o A A A 0 -Bcc -Bc| {~}c {U}c AIa AIc A B - -Bc -BII {;e}I _ {U}I Baa 0 BaI A~a 0 AA i {;/ }a {V}a 0 Bc Bc 0 AA AlC {}}C { Bia Blc Bjl ALa AAL AL {ib,}I {V}i + [ VeI jNte hmtd dl ~~ -IC N h<, j V, em T +0 e frc jN hmt dl 0 0 - Ne=l frc jN emtdl 0 0 (6.69) where Ne [A']kl = er [a]el (6.70) e=l Ne [A]kl = >,,. [a]e (6.71) e=l Ne [B]kl = [b]kl (6.72) e=l Nc Uk = ] {U}k (6.73) e=l where the ijth member of the eth local element matrix is related to the klth member of the global matrix by a node-element connectivity transformation function p as k = p(e, i) (6.74) 1 = p(e,j) (6.75) In (6.69) the subscript I refers to group of nodes in region the Q excluding r, and the subscripts on A"', B, U and V refer to the various regions of Q and its boundary.

109 For example, AI, refers to the matrix resulting from the interactions between nodes in region Q and on rF, whereas AC, refers to the matrix resulting from the interaction of the nodes on r,. Each matrix is sparse, having nonzero elements when the nodes share a common element. Note also that V is the dual of U, which was defined in (6.51). Two options exist at this point for evaluating the contour integrals. They may either remain on the excitation side and the tangential modal fields expressed in terms of a condition on rF, or the tangential fields may be expressed as unknown scattered potentials and moved into the matrix. The second is required for the FE-BI formulation and is given by {4e}a ACa 0 A, -Baa 0 -BI -C { {U}a 0 I O O O 0 Caa 0 {0e}I -{e}c Aa A A" -BIa -Bc -BI 0 0 { h}a {U} Baa 0 BaI APa 0 AAi 0 0 {0 } c V}a 0 BCC BcI 0 AAC AA 0 0 {Ifi}I {V}c BIa BIc BlI AA ACA 0 0 {/OSt}a {V}I I Oht IC +[ el fr. Nie(jhm )dl 0 0 -S I fra N(i emt)dl 0 0O] (6.76) where the boundary conditions Fe, = 0 and emt = 0 on rc were also enforced. Additionally, C[C Ckt =Ia'ie'R die (6.77)

110 and 3 =et = jRemt = E NJ {4} (6.78) j=1 3,t jRh = jRh NJ {fbOt}; (6.79) j=l This form is now suited to augmentation by either a discrete form of a boundary integral equation explored in section 6.3. 6.2.3 An Improvement to the Finite Element Formulation The integrand of each of [a]TJ and [b]s- in (6.49) and (6.50), respectively, becomes singular when R = j~m for real /. Since we are only concerned about solutions for which m > 0 (those for m < 0 are found via symmetry considerations), only the positive sign is considered, or Ro = m. The location of the singularity is independent of Z and is hereafter termed "line singularity." The line singularity intersects any element e containing the radius Ro (as seen in Fig. 6.3 for a homogeneous medium), or is near an element if Min(lRo - Rs I) is small (not defined now) for the normalized radii of element e, R. for i E [1,3]. The origin of this singularity has been discussed previously in [23] for the CAP equations. For elements containing Ro, [a]~, and [b]s, gain an additional residue contribution, automatically included by the use of the logarithm of (6.62) in the sums (6.60) and (6.61). However, for elements oriented nearly parallel to the line singularity, as in element b in Fig. 6.3, the contribution to (6.60) and (6.61), and consequently [a]3 and [b]~,, become large. In the case of element a in Fig. 6.3, they become infinite. Large matrix elements are associated with the slow convergence of the conjugate gradient solver and inaccuracies in the resulting solution. Presented here is a way to avoid this problem.

111 To understand the nature of the problem, we recall (6.23) d x hm(R, Z) = jfm [mq X Vtbh + -rRVtbe] (6.80) the right hand side of which appears in (6.44) and ultimately in [a]- and [b]- after discretization. Since q x hm must remain bounded everywhere in space (barring edges of perfect conductors), the right hand side of (6.80) must also. Clearly then, the bracketed terms involving,e and 1Ckh must combine to cancel with the singularity in fm = 1 - The finite element discretization, however, separates these (KR+m)(K.R-m) ~ two terms and each will individually become large when an element edge is nearly coincident with the line singularity. To regularize [a],e and [b]Ti, we first define the integral I = Jffm(rRVt-) VtNidS + fm (X V ) VtNidSe (6.81) Se Se which appears in a portion of (6.44). Each of the integrals in (6.81) must be regularized, and in doing this we must subtract and add a term from each of the form h(Z,R) (6.82o) KR- m where h(Z, Ro) is the factor of the integrand which is well behaved at R = Ro. Applying this technique to (6.81) gives I = || 1 r rERVts. VTNe _ (Z, Ro)RoVtg'b((Z, R). VN(Z,Ro)] dS1'fJ OKRIm LKR+m K/cRO+mt IJ /1 KR-rm VKRt + V N KRo0 +r VtN(Z Ro)] dS 1 _ ( Zi Ro)RoVttbS(Zp Ro) + mR x V+tb (Z, Ro) $sJJ tcR-mm L~Rm + m VtN i (v e +JJ K KZ, Kor.VNi(Z, Ro)dSe where h(Z, Ro) for each integral is clear by the argument (Z, R0). Where not explicitly shown, all functions exhibit a (Z, R) dependency.

112 The key to implementing this method is the elimination of the third integral in (6.83). Note that the numerator in the integral of this terms is simply the bracketed factor in (6.80) evaluated at R = Ro. As mentioned before, X x hm(Ro, Z) must be finite, the bracketed term is zero at R = Ro. Thus, mO X Vth(Z, Ro) = -r,(Z, Ro)RoVt4e(Z, Ro) (6.84) The third integral of (6.83) vanishes when the expression is substituted into its integrand. After following the usual discretization of the azimuthal potentials, (6.83) may be written in discrete form as r / 1 i ERVtNje R. VeN 1{esi i;R-m LK R+m r:(Z, Ro)RoVtN (Z, to) VtNe(Z, 1o)] dS + y.{- m x x V,N. + h} R-m [ cR + m ~j=1 eRRo+m] +x 7vtNje(Z, R~) V Ne(Z, R) e (6.85) For linear shape functions, the gradients are constant (as well as Er) and are factored out to yield j=1 SJ 2rK(KR+ m) 3 I+ = I}jJ 2,m(KpR + mj X VtNj VtNetdSe j=1 se (6.86) which are clearly regular when R = m. Thus, the matrix elements in (6.49) and

113 (6.50), may now be written [a/e /[= II V21( + VNje VtNe + Nfj]dS (6.87) 2i;(K;R + m) R se i ff 1 I m x VtNe. VtNiedS' (6.88) Ltij $I 2(iR + m) These were implemented, but did not work very well. Note that since the potentials 0/ (Z, Ro) and 0'h(Z, Ro) were expanded in terms of linear functions, their respective gradients are are constant everywhere in the element independent of R and consequently independent of Ro. Since the derivatives V1e and Vzh are constant over the element, this approach to regularization is of low order. To improve the accuracy, the functional dependency of this on R must be increased and for the existing discretization scheme, this may be done by refining the mesh in the vicinity of the singularity. This, however, is not a realistic option for frequency sweeps for scatterers containing lossless materials. Using higher order shape functions or employing the modal field as the working variable (i.e., Vb, = V(Re,g,) = RVemr + e,sVR) would alleviate this problem. It is interesting to note that given Ie = Rem,, a linear variation in Z of em, corresponds to the same in e,. However, a linear variation of e,m in R corresponds to a quadratic variation of the same in e,. But since Ike is expressed in terms of linear shape functions here, the radial behavior is lost. In summary, then, we have developed a method by which to alleviate the problem associated the line singularity at Ro =, for real ic. This not only eliminates the need for a residue contribution to [a]s. and [b]-, but also prevents them from ever becoming infinite. In addition, a direct consequence is that these element matrices are now purely real for scatterers containing lossless materials, which is not the case for the traditional formulation presented in section 6.2.2.

114 6.3 Boundary Integral Formulation In the previous section, a weak form of the wave equation was developed for scattering bodies in Q and bounded by r =,a + rc + r + The boundary conditions on rc and rF have already been included in the system (6.76). The resulting discrete system remains incomplete, however, since the boundary conditions (which provide a relationship between the tangential E and H fields) on the exterior boundary F, remain unspecified. The Stratton-Chu integral equation (S-C) provides the necessary relationship for source and integration points along ra, when the integral equation is expressed in terms of fields tangent to the surface of revolution. By employing a Fourier series expansion of the S-C and thereafter discretizing it on Pa, the resulting system provides the needed equations for solving (6.76). As a preliminary step, the electric and magnetic fields in the unbounded region are represented by E(=) = ET8() + 7 (T) (6.89) 7H(-) = H (Y) + H'() (6.90) where 7(T) and W(F) are the incident fields and the scattered fields are given by the Stratton-Chu equations [22] ES(r) = } {-jko [n' x yoH (i)]g(T, ) + k [- [ V' x q o H (T)] Vg(T g ) S7it(T) = fi{jko[h'x 7S( )]g(T9) - i [kvn V x ES()]V'g(rt) +[h' x Eo'(F')] x V'g(,') }dS' (6.92) I F(;4 (f y4 VXES(4)Vg -l14 77077.9~~~~~~~~~~ j(6.02) nr

115 where A' and T are the source and observation points, respectively and e-jko Ir-A' I g(,rt) e 4 1- l (6.93) is the free space Green's function. In this form, all field quantities are tangent to S'. Since in the development of the discrete boundary integral system the source and field points reside on the surface S', it is convenient to remove the singularity at S = T' by expressing the integrals in (6.91) and (6.92) in terms of their respective principal values as () f{jko[n' x qo'()] 9(()+ [n V x oH (T)] V'g(r,) + [ x Es(t)] x V'g(r,7)}dS' (6.94) HS7r) = 4{Jk0 [f' x E (j)] g(r,') - ko [ + [n' x 7BH()] x V'g(rf,)}dS'(6.95) Looking back to (6.76), it is clear that two additional equations relating the unknown potentials O?/, tk9, et9 and /t/ are necessary to form a complete system. This may be achieved, for example, by using the ti component of the modal form of (6.94) and the q component of the modal form of (6.95). In fact, many combinations are possible but for simplicity and symmetry, the q component of (6.94) and (6.95) is considered below. 6.3.1 Derivation of the Modal Boundary Integral Equation In a fashion analogous to the finite element formulation, the development is provided explicitly for (6.94) and duality is employed to obtain the corresponding expression for (6.95). First, consider the general surface of revolution indicated in Fig. 6.1 whose tangential unit vectors are denoted by X and t. The unit vector i subtends

116 an angle v with the z-axis, and t' subtends an angle v' with the z-axis. With reference to Fig. 6.1, the various unit vectors are given by n = cos v cos 0 + y cos v sin - sin v (6.96) = - sin + y cos (6.97) t = sin v cos + ysinvsin + z cosv (6.98) = tsinvcos + cosvcos - sino (6.99) y = tsinvsin + ncosvsin + 0cosq (6.100) z= icosv-nhsinv (6.101) Expressing the primed unit vectors in terms of the ones results in i = t' [sin v' sin v cos(4 - q') + cos v cos v'] + ni' [cos v' sin v cos(4 - q') - cos v sin v'] + q' [sin v sin(q - q')] (6.102) n = t' [sin v' cos v cos(q5 - 0') - sin v cos v'] + n' [cos v' cos v cos(4 - 0') + sin v sin v'] + 4' [cos v sin(0 - 0')] (6.103) 0 = -P sin(o - 0') + O cos(q - 4') (6.104) which clearly reveal the dependency of the unit vectors on the azimuthal variation - q'. Further, the following identities can be shown to hold: 0 V'g -= Vg (6.105) * (n' x toH7T) = -7roH; sinv'sin( - q) - toHt' cos(0 - q') (6.106) n''. (V' x o0Hi) = I [-(p'H) + (H] (6.107) q. [(*' x E ) x V'g] = [,'Et sin(q - q') + n'Eg cos(q - q') + q'E cos v'sin( - q')] V'g (6.108)

117 Using these in (6.94) gives'E(=) I fri j ko[7oH; sinv'sin(q - q') + 7oH: cos( -')]g(L,') jkop' [-I PoH, + a~(r0H')]. Vg(r,') + [)z'Et' sin(q- -') jkcop' [-~r~rg~frr~rX rAV + n'E; cos(q - q') + q'E cos v'sin(q -')] V'} p'dd'dt' (6.109) and by carrying out the derivatives of the Green's functions, we find E2's) = r I o {jko [7loH; sin v' sin(q - O') + yoHt' cos(O - O')]g(v, T') 1 ( dg - jko [-t (P 7oHs) + a~ (770Hts)] sin(q;- O')Ro dRo + (-E(z-z)sincos(')[(z z')sin v' - ( - p')cos v'] - 2pcos v' sin2( 2 0)}HR d-g p'd'dt' (6.110) in which Ro Jp2 + p2 2pp' cos(q- q_') + (z - z')2 (6.111) The integrand of (6.110) is expressed explicitly in terms of q and 4' and is now suitable for harmonic decomposition. To generate the corresponding modal integral equations in terms of the Fourier coefficients, the fields and Green's function may be expanded as 00 E~()= Z e -(p,z)ejmk (6.112) m=-oo 00 1oH'(r)= E h1m(p,z)ejm4O (6.113) m=-oo oo00 g(k) k)(p, p, z, z')ejn(O-O') (6.114) n=-00 where e(p, z) =- EE(p, u, z)e-J"mudu (6.115) 2'/!_

118 hm(pz) = / -J H(p,,U z)e-jmudU (6.116) (p, Z) 27r (pU,_ 1 fr e-jkoR gn$o)(Pp', ZZ) = - JkRcos(nu)du (6.117) ~ o 47'R 1 ir e-jkoR gn)( p', z, z') = - cosu u cos(nu)du (6.118) ir -jkoR gn)(p,', zz) = sin u R sin(nu)du (6.119) i'zr 1 dg g'n)(pp' P, Z') = oJ o f dR cos(nu)du (6.120) 1 f~ 1 dg g() (pp z) = k2 COS UdR cs(nu)du (6.121) (2' i, = 1 dg g(2)(P, p', z,z')= sinu sin(nu)du (6.122) r ko Jo R dR in which R = p + p2 - 2pp' cos u + (z - z')2 (6.123) Substituting these into (6.110) yields eZ (p, z)e ndn = j en] + {jko [hm sin v'g(2) + hamtn 1)] n=-oo n=-oo m=-oo ra J - [- [-(p'hs ) + jmh t] kog(') -et(Z - z')kog2 + esmko(p' cos - pcos v'g + (z - z') sin v'g)')} eJ(m-n)'p'd'dt' (6.124) where t' is the parameter along the contour ra and increases in the clockwise direction as shown in Fig. 6.2. Further, upon multiplying each side by e3P'0 and integrating over (0, 27r) to extract the mth harmonic equation gives emk(p, z) = 27r Jr0 {jhh sin v'g$ j( )-(P ) + (jh')[g () + jmg +- e~,~ko(p' cos v'g$)' - p cos v'gm + (z- z')sin v'g$) ) -emt(z- z')kog$) )}kop'dt' (6.125)

119 after also combining similar terms and where we have used the identity 27r ej(m-")' =d 2r m = n (6.126) 0 otherwise Introducing the normalized coordinates R = kop R'= kop' Z= koz Z' = koz' a 0 A9 t ko/ 5,-,(6.127) Or' = k Or' in (6.125) and replacing the field quantities with the azimuthal potentials yields 2Re =ko f r{ h O sin v gm ] ( gt] - h) [j Rg ]2)' + Ot [g$l) + jmg2)'] ko +?r [Rcosv$(1) + (Z()'] Os [R' c'- R cos v'g'c + (Z - Z') sin v g,)' + t [j(Z - Z')g()'] }dr' (6.128) This equation and its dual form a pair of integral equations to be imposed at the mesh termination boundary. Their discretization is considered next. 6.3.2 Discretization of the Modal Boundary Integral Equation In discretizing the modal boundary integral in (6.128), the contour rF is first divided into Na elements each of length Ae. Within the eth observation element, the parametric representation R = RR + r sin ve (6.129) Z = Ze + r Cos v (6.130) is adopted as shown in Fig. 6.4 and is consistent with the required counterclockwise path traversal of Pa. Likewise, within the source element e', the parametric

120 element e element e' R /N2() R A A A n ite representatione' Na 2 e'=1 j l N1(2>~~~~N 1 T (6.134)Ve N N e(T 2 ^ employed in the finite element discretization of the boundary integral equationof weighted IN 2 el=lej=l are linear parameterized functions of a'. (Note that the linear basis functions are

121 residuals over the eth observation element and Galerkin's technique to (6.128) gives Na 2 A-e NN 2~r Na 2 -- jl Z {~ R e ko, RNdT + p Z Z f {fbh}" N -N'J sin v'g()] e'=l j=1 e'=1 j=1 -{i N,(8riNf) [jR'gm']2 + {IVt}' N.Ne' [9m) + jmg,)] O's Nie Ne' R tg i t + {sk}' NeNj [R' cosv'g (1)' - Rcosv + (Z Z') sin v'g)'] 1+ ({Cj}' NeNj_. [j(Z - Z')g,)'] dr'dr = 0 (6.136) In deriving the first term in equation (6.136), it is useful to note that e(r.) may be written ~b,(r) = 4Ir b(r')b(r - -r')dr' (6.137) where S is the Dirac delta function. Substituting (6.133) into (6.137) yields Na 2 =eT Ej ~{be}' Ir N (')S( - r')dr' e'=l j=1 a VNa 2 = Z Z{ib}je'Nfe'(rT) (6.138) e'=1 j=1 Taking the inner product of (6.138) with N yields Na 2 r Ne(r)Ne'(r)( 4sti(7r) dr = ~ b} Idr (6.139) R d = E f O{ j....... e'=1 j=l a as expected. The system (6.136) may be written in compact matrix notation Na 2r z z[ [L6]f {t:}j + [M;]j {h}j; + [Ltj {st}j' + [Mt], { };htj = 0(6.140) e'=l j=l where we easily see that [L]tj' = | dr + [R'cosvg$) - R cosv gm +(Z- Z')sinv'gO)'] dr'dr (6.141)

122 [M>] = - ~f {N. NJ N'[j sin v'g)] - Ni(Nf') [jR'g)']}dr'dr(6.142) [t]j = j N Ne' [j(Z - Z')g('2)] dr'dr (6.143) [Mt] f NeN' j[g() + jmg(2)'] dr'dr (6.144) Each matrix of the form [A],te' can be termed "element interaction matrix" (after the finite element term "element matrix" for the case e = e'), since it represents the interaction between the source element e' and the observation element e for i, j E 1,2 of the local weighting and expansion functions, respectively. The numerical evaluation of the integrals in (6.141) - (6.144) is performed by breaking the integral into two integrals with continuous integrands and performing midpoint integration for each of them. This applies to the integration in r' as well as for T. When the source and observation points coincide, an average value is computed by moving the observation point / away from the sub-cell center. This has been shown to work well in [25], and avoids the need for elaborate self-cell evaluation techniques. In all cases, Gaussian quadrature is employed for the q' integration in (6.117) - (6.122). The local "interaction element matrices' may be assembled to form a global boundary integral system in a fashion analogous to the finite element method by first summing (6.140) over all observation elements e as Na Na 2 z[ [L]{l }' + IMp]e' 0`' + ] +LJ= 0 e=1 e'=l j=I (6.145) where the local element subscripts i, j have been replaced by their global counterparts k, I and are related by the node-e lement connectivity transformation function q as k = q(e, i) (6.146) = q(e',j) (6.147)

123 A(k,l) = 0 k,l E {1,Na} Do e = 1 to Na Do e' = 1 to Na Do i=1 to 2 Do j=1 to 2 k = q(e, i) = q(e', j) Form [A]". A(k, 1) = A(k, 1) + [A]e' End Do End Do End Do End Do Table 6.1: An algorithm for the boundary integral system assembly for a generic element interaction matrix [A]~' In Table 6.1 is illustrated an algorithm where a generic matrix [AO]eel (where A represents any of L+, MO, Lt or Mt) is assembled to form a global system. This algorithm differs from that of the finite element method, in that the loop in e' is eliminated. The final FE-BI system is formed by augmenting finite element system with

124 (6.145) and its dual giving ACo 0 Al -Baa -B 0 -Caa } {U}a O I O O O O Caa O 0oe }c - C Aa AC A BIa - BBIa - BII O O {Se}I {U)i Baa O BaI A"a O AA 0 0 {?.I}a | {V}a o0 m= BcI 0 AC Al1 0 0 {h } { V}c BIa BIc BIjI Aa AIc AA O0 { h}I {V}I L,4 M, 0 0 0 0 Lt Mt {/et }a 0 -Me Lx 0 0 0 0 -Mt Lt { t}c + E[ fra Ne(jhmt)dl 0 0 - e=l fr N (j et)dl 0 0] (6.148) which can be solved by the conjugate gradient method. Before considering different solutions of this system, we must develop the harmonic coefficients of various sources of excitation. This is considered next. 6.4 Sources of Electromagnetic Radiation The modal forms for a plane wave excitation and an electric dipole source are developed in this section. The former is employed for scattering problems while the latter may be used for both scattering and radiation. For both cases being considered, the resulting expressions for the azimuthal potentials are used in the expressions developed in section 6.2.

125 6.4.1 A Plane Wave Excitation A plane wave incident at angles (9i, q0i) and observed at r = (r, 0, z) (see Fig. 6.1 ) has the form us,(Oti, qi; p, ~, z) = O'e-jko%' (6.149);U(Oi, 0i; p, q, z) = oie-jkor (6.150) where qY is perpendicular to the plane of incidence and O' is in the plane of incidence. Expressions (6.149) and (6.150) may be written as explicit functions of q by first writing the argument of the exponential function as lo * _ kor( * r =-ko [psin O' cos(- i) + zcos i] (6.151) since r = sin 0 cos + Y sin 0 sin q + z cos O (6.152)' = x sin O' cos qi + } sin Oi sin i + z cos Oi (6.153) Also, making use of the unit vector transformations i= -= sin O' + Y cos 0i (6.154) = i cos 0i cos qY + Y cos 8' sin qY _ Z sin O' (6.155) = p cos O - sinq (6.156) y-= sin~ + +cos q (6.157) =z z (6.158) (6.149) and (6.150) can be rewritten as T,(0'; p, q-',z) = [psin(- Y) + q$cos(q$ - 0) ] eiko["sin'i cos(d-')+zcos~']

126 (6.159) uie(Oi; p, - qi, z) = [ cos Oi cos(5- 0i) - cos O' sin(o - qi) -' sin 0'] eiko [P sin O' cos(k-, 0)+z cos 6i] (6.160) These may be expanded into a Fourier series in the parameter (q- qY) by first writing (6.149) and (6.150) as 00 vU(09; p, c Z - z) = m(; p, z)ejm("-') (6.161) m=-oo 00.oiie(; p, q - q9, Z) = umoe(O9; p, z)ejm(-'i) (6.162) m= -00 Each component of (6.159) and (6.160) may be expressed in terms of one of the functions f(0i; p, / - ~,) - ejkopsinG' cos(4-4') (6.163) fc(8O; p, q - bi) = cos(q$ - Y)f(90; p, q - q ) (6.164) f8(Oi; p, q _- i) = sin( - _i)f(0i; p, q -_ qi) (6.165) Expanding each of these into a Fourier series in (4 - i) and using the even/odd properties s(q) = s(-q) ) sm(u) = s_m(u) for f,fc (6.166) S() = -S(-q) ~ Sm(U) = -S-m(U) for f, of the functions f in (6.163) - (6.165) may be expanded as 00 f(O9; p,'- q ) = fo(Ot,p) +2 Z f,(O,p)cos[m(d - qY)] (6.167) m=l f (Oi; p, d~ - ~b) = fo(8i, p) + 2 ~ f~ (0i, p) cos[m(b - ~i)] (6.168) m=l o0 f(i; p, q - q6) = 2j ~ fsm(Oi, p) sin[m( - q )] (6.169) m=l

127 with the corresponding Fourier coefficients 1 oc fm (, p) = 1 j ekopsin-' cosU cos(mu)du (6.170) fcm (i, p) = 1, cos uejkoPsinG' cos(mu)du (6.171) fsm(Oi, p) - | sin uejk~opsinO cos u sin(mu)du (6.172) 7r Comparing these to the.Bessel function identities jmJm(/)= _I ejpc~s cos(mx)dx (6.173) jm-1 JI() = cos | xCos cos(mx)dx (6.174) _3 mJm(3) = -* j| sin xejPC~O sin(mx)dx (6.175) where the last two are derived from the first by differentiation with respect to, and integration by parts respectively, we conclude that fm(Oi, p) = jm Jm(kopsin 0i) (6.176) fm(Oi, p) = jm-1Jm(kop sin 0i) (6.177) f m(0i, P) = - kop ifm (0', p) (6.178) Using these, along with (6.163) - (6.165) in (6.159) and (6.160), the Fourier coefficients of (6.161) and (6.162) are found to be Um (0i; p, ) = ejkozcosO' [pfm(8i, p) + -fcm(9ip)] (6.179) ime(80; p, z) = eJkoCOB [ f_(oi p) cos 80 -Afsm (~ p) cos 0 - fm( ) sinj C (6.180) or, using (6.166), (8h;pp, _ 4,qz) - ejkozcos, [2j E fm(9',p)sin[m(4- 4i)] m=1

128 + gfco(Oi, p) + g2 z fcm(O, p) cos[m(q- qi)]] (6.181) m=1 u-e(oi; p, 4 - 4', z) = ejkozcos8i pcosOa [fo(Oi, p) +2 E fcm(Op)cos[m(-i)] -Acos0i [2j Z fm(Oi,p)sin[m(q- i)] (6.182) m =l - z sin 9 [fo(9', p) + 2 E fm (0, p) cos [m( - )]J j The modal coefficients u'im and ~UM can now be used for determining the harmonic coefficients of the incident and magnetic fields. For TEz incidence, we have'm — = Um) O(6.183) hm = Umo (6.184) and for TMz incidence em = -Umo (6.185) hm = UmO (6.186) More explicitly and in component form em = fj(i, p) ejko z cos (6.187) ht = fm(Oi, p) Cos iO ekozcosi (6.188) emt = fsm(O, p) sin v e jkoosO (6.189) ht = [fm (Oi, p) sin 8O cos v - fm(Oi, p) cos Oi sin v] ejko~coO' (6.190) for TEz polarization and e = - fsm, (Oi, p) cos 0i ejkozcos8' (6.191) him, = fcm(Oi;p) eskozcos~' (6.192)

129 emt = [fC (i, p) cos 9i sin v - fm (', p) sin 9i cos v] ejkoz cos' (6.193) ht = fsm (Oi, p)sin v e 3kzcs (6.194) for TM, polarization, where v is the angle between the z-axis and the vector t in the p - z plane. The corresponding potential forms are obtained by multiplying the previous expressions by R = kop. 6.4.2 An Electric Dipole Source An i-directed electric dipole of moment x4lre/k3 located at a point zo on the z-axis is given by [26],e-jo [kr(kr -kzO cos) ) 3 3i j 1] (6l95) =e-iRp4 (kkz q- + 1 ) P + 1 sin 0 cos - (6.195) e-j& [krkzo sin 20 3 3j 1 j 1)os0]cos (6 196) -ro [ 1 ] (6.19 i) E - Ro [R + -1 sinoq (6.197) yHe =- j ]k (krk - kzo cos 0) sinn (6.199).e-j~ [Ro+ jV1 ](kr cos 0 - kzo)cos (6.200) where Ro = kr2 sin 0 + (r cos 0 - zo)2 (6.201) When using this source in connection with the current implementation, the harmonic coefficients of the fields must first be determined. This may be done rather trivially by noting the identities sin b = - (6.202) cos = + (6.203)

130 Introducing these into (6.195) - (6.200) and comparing with (6.161) and (6.162) we find (on setting i' = 0) that e3 kr(k -1 3 3j 1 emr [r(r - zo cos 8) + - 1 sin 0 (6.204) e-j [krkzo sin2 3 3] (6.205) emo 2= [ (i - 1) -( - -cos0 (6.205) 2Ro Fj 1R1 e-iS [j~o I em = -j 2RRo + -1] (6.206) hmr 2~= L 2+ kzo sin 0 (6.207) hm 2Ro o + ](kr - kzo cos 0) (6.208) hm=j -oRo ]+ (kr cos0- kzo) (6.209) where m = ~1. For all other m, the coefficients are zero. Additional field components needed in the finite element portion of the system are emt and hmt and are given by emt = emr cos(O - v) - emo sin(O - v) (6.210) hmt = hmr cos(O - v) - hme sin(O - v) (6.211) where the identities =r tcos(O - v) + i sin(0 - v) (6.212) = -tsin(O - v) + tcos(0 - v) (6.213) have been involved in deriving (6.210) and (6.211). 6.5 Scattered Field Computation Once the modal coefficients are determined from the solution of (6.76), the goal is to proceed with the determination of the scattered or radiated field. In the following, the evaluation of these from a knowledge of the modal scattered field coefficients is discussed.

131 The scattered fields in the far zone are given by Es(r) = E+(r)+ + i7oHk(r)O (6.214) 77H0FS() = o70Ho(r)A - Eo(r)0 (6.215) These involve the q components of the electric and magnetic fields and these are related to the modal field components via the Fourier series as M E;(r,,,q) = 2j >E es(r, 0) sin(m) (6.216) m=l M o10H,(r,0,, qf) = hs(r, 0) + 2 > hS,(r, 0) cos(m0) (6.217) m=l and that for the 0-polarized receiver, and by the series M E;(r, 0, q ) = e)(r, 0) + 2 E e,(r, 0) cos(mq) (6.218) m=l M 770oH,(r, 0, ) = 2j E h, (r, 0) sin(m0) (6.219) m=1 for the +-polarized receiver. For a unity amplitude incident field, these imply the radar cross section formulae (2.33) r=lioo4n 1 e 2 E+ em (r,0)cos(mq) (6.220) m=1 M 2 O14 =libmceo7 ho + 2 E, h',(r, 0) cos(mq ) (6.221) m=1 a = lirn + 2 4r 0(r, 0)sin(mq) (6.222) M=1 1 M 2 -= lir -2 E (r, 0) sin(mo) (6.223) r —1oo 47r m=1 where es = 4irr es, (6.224) hea = 4rr hr (6.225)

132 The evaluation of e and hm, in the far zone using the near field values. The Stratton-Chu integral equation used for calculating the scattered fields was discretized in section 6.3.2. For the potential distribution computed by solving the system (6.76), and these values are known on any path rf passing through the solution region ~ and which encloses the scatterer. By subdividing rf into Nf contour (boundary) elements, the scattered field may be expressed in a form similar to (6.136) as e~6= p P Q ] Q] ht (6.226) ht where Nf 2 P, = Z -'{P,}e' (6.227) e'=l j=l Nf 2 Qo= Z EZ{QO}' (6.228) e'=l j=l Nf 2 Pt= Z E{Pt}j (6.229) e'=l j=l N1 2 Qt E= Z {Qt}X' (6.230) e'=1 j=l and in which Ip 2r(1)' {P};' = Irk j N' [R' cos vg( ) R cos v'gm + (Z - Z') sin v'gm ]dr' (6.231) {Q} = j {N [j sin ()]-( aN) [jR'g)'] }dr' (6.232)'.e' Ae' = N;' [j(Z- Z')g)']dr' (6.233) {Qt;= -- j N;'[g) + jmgg)' ]dr' (6.234) \~~sstlj Ic o I

133 which must be specialized for far zone (r - oo) computations. In deriving expressions for the far fields, the argument of the Green's function is replaced by its usual far zone approximation R = /R2 + R - 2RR' cos u + (Z Z;)2 kor - Z' cos 0 - R' cos u sin 0 (6.235) where Z = kor cos 0 (6.236) R = kor sin 0 (6.237) The Green's function and its derivative may then be written as ei r ej(z' cos +R' cos u sin ) (6.238) 4rr O'- 9 (6.239) R OR kor and when these are used in (6.118) - (6.122), we obtain g4) 4rr f m(0, R') eiZ'cosO (6.240) (2) ~e-jkor ( 6 21, J r-jkor m2) -4ir fsm(0, R') ejz"ose (6.241) m -~ 4r fe(k or) ejZ'Cs (6.242) kor 4irr gm kor 47rr fm (, R') ejz S (6.243) k2r e-jk~ s gm2) f( R') ec jZ S (6.244) kor 47rr where fm(0, R') = jmJm(R' sin 0) (6.245) f (0, R') = jm-lJ (R' sin 0) (6.246) fsm (8, R') = R' sin 0fm (0, R') (6.247)

134 element e boundary element e' {R)e {R }1/2 Figure 6.5: The configuration of a typical boundary element e' passing through triangular element e Furthermore, substituting (6.236) and (6.237) into (6.231) - (6.234) yields {pj - eIkoj | e Nf(r')[cos 0 sin v'f, - sin 0 cos v'fm] eiZ'c osdr/ (6.248) {P)''- 2kor Jo *'~j e2k or Nje (r')sin v'f,,m ejz'cos~drT (6.249) {Q} J.e-jkor f {Pt}' -i 2krO Nje ()cosfs,, ejZ' CoS drT (6.250) {QtI Nj (71)fcm ejZ'cs~dr' (6.251) If e't and hst are not known on rf, then they must be computed by numerically differentiating the potentials for each triangular element traversed by rFf. For convenience, it is assumed that the e'th boundary element passes through the eth triangular as indicated in Fig. 6.5 (i.e., the boundary element must begin at one edge of the triangle and end on another). In this case, ers is found by first taking the dot product of ii' (the outward normal of rf) with (6.22) to obtain jemt = -fm[m t. Vte - R. * V~h] (6.252)

135 Since in the eth triangular element 3 {/p}e = E Nje {} (6.253) j=1l we have 3 Vt {I} = VtNe{Ie} (6.254) j=l where =e + e + 7R(6.255) VtN / + 72e (6.256) Further, on noting the identities tP' = cos v' t'.' = sin v' (6.257) n' * =- sin v' n' * = cos v' the derivatives can be written more explicitly as ^'. Vt{~~}' =1 E~3 -#je sin v' + yj7 cos v' into = Z{} (6.258) j=1 {t f e = 1Z im_ (6.259) {R, {j,t MV = te * + n V (626) j=1 A similar procedure is used for. V and. Vt{ and substituting these into (6.252) gives tR~e' In~,J ~el'' }j+}A R} CV=s -{R}R [ml' Vt - {R}%+f' Vk](6.260) ~ m2 {rj+' = {j+ R}1 [ml'. Vt3 + {} 1' iV:] V9(6.261) [{2j+] 2- m2 h where {R;1} is introduced to express the left hand side as an azimuthal potential. case of F1 coinciding with Fa, then the e~n and hSt are computed during the system solution. Then the need to compute the i- components is unnecessary.

136 6.6 Code Validation and Performance In this section, the numerical implementation of the FE-BI method is evaluated in terms of storage and computational efficiency and accuracy. Also, the method for eliminating the boundary integral resonances presented in chapter V is tested here for spherical enclosures. Finally, the scattering from several structures is presented for validation purposes. 6.6.1 Storage Efficiency This implementation of the FE-BI method discussed in this chapter differs from the two-dimensional FE-BI formulation discussed earlier in that arbitrarily shaped mesh terminations are allowed. As a result, an increased storage efficiency is not obtained from the convolutional properties of the boundary integral operator, but is instead realized by exploiting matrix symmetries. It is, nevertheless, understood that the controlling factor in the storage of large FE-BI systems is the dense boundary integral subsystem, which grows as O(n2) for n unknowns on the boundary, while the interior FE system grows as O(N) for N unknowns in the interior region. It is possible, however, to choose an enclosure on which some of the integrals are convolutional. The Green's function in (6.93) appearing in the integrand of the integral equations (6.91) is a function of the distance between the source and field points, which when expressed in cylindrical coordinates takes the form - = p2 + p'2 - 2pp' cos( -') + (z -')2 (6.262) After employing the Fourier series expansion of the field quantities and Green's functions, the azimuthal dependency is removed and the result is given by (6.128). This expression is a convolution only on contours parallel to the z-axis. Consequently,

137 a suitable enclosure is rectangular in shape, generating a right circular cylinder in three-dimensional space. Though analogous to the two-dimensional rectangular enclosure explored in chapter III, the integrals are in the form of convolutions only when source and observation points lie on the portion of ra which is parallel to the z-axis. All remaining terms are considered "cross terms" and must be stored efficiently, relying on symmetry to achieve this goal. The discrete version of a typical boundary integral equation block matrix (of Lt, say) has the form all a12 a13 a21 a22 a23 (6.263) a31 a32 a33 w lote the observation on side p and source on side q of the rectangular boundary (the z-axis is not applicable to the boundary integral). As a consequence of the convolution, a22 is Toeplitz in form and, effectively, only the first row need be stored when an FFT is employed to evaluate the associated matrixvector products. For a22 to dominate the storage requirement of the BI system, the enclosure must be long with respect to its radius. For scatterers not satisfying this requirement, the reduction of a22 is not substantial. Consequently, the method was abandoned in favor of a general conformal boundary. To determine the storage requirement for a general boundary, it is first noted that each matrix in the BI subsystem of (6.148) need only be stored once. Additionally, after examining (6.141) - (6.144) Mt and Lt are the only symmetric matrices. Thus, for n nodes on the boundary the storage for the each of the symmetric matrices is n(n + 1)/2 and that for the unsymmetric matrices is n2. The overall storage requirement of the BI matrix is then n(3n - 1). Though this continues to be 0(n2),

138 the constant in front is reduced from 8 to approximately 3. For the particular case of a purely smooth metallic scatterer, the full storage requirement of the standard surface integral equation (SIE) for BOR structures is 4n2, slightly less than that of the FE-BI system. However, for structures such as a corrugated cylinder, n becomes very large for the SIE but in the case of the FE-BI method, the termination boundary may be chosen to occupy a minimum path and thereby minimizing the associated storage requirement. In passing we note that it is possible to optimize the storage by choosing the termination boundary far enough from the scatterer so as to reduce the sampling the requirement of the boundary integral (which grows as 0(n2)), but without substantially increasing the storage associated with the additional elements required in the finite element region (which grows as O(N)). This is particularly useful for structures containing very thin high contrast material layers, in which case a strictly conformal boundary would result in a large number of unknowns on the boundary. For a given amount of data storage space, S [bytes], the maximum allowable length I of the termination boundary depends on the uniform sampling rate A as I = nA (6.264) The single precision (8 bytes) storage requirement for the discrete boundary integral thus becomes S = 24() (6.265) and then solving for the length, I= A 4S (6.266) Thus, if S = 10 megabytes and A = A/20 l _ 32A (6.267)

139 This represents an upper limit on I given the memory and sampling rate without regard to the finite element storage nor the overhead associated with the geometry and solution storage, which mutually grow as O(N). Consequently in practice, n = must be reduced until the storage of the entire model is accomodated, a requirement which depends on the scattering body. 6.6.2 Computational Efficiency and Accuracy The FE-BI system given by (6.148) is solved via the conjugate gradient method. It is well known that the rate of convergence is proportional to K2, where rc is the condition number of the matrix. Clearly, those factors which influence C must be examined, since in part, cpu time of the solver and subsequently the accuracy of the solution depends on this. Since the FE-BI system is comprised of two incomplete subsystems, it is expected that the weakness from each portion contributes to a reduction in Kc. Among those due to the finite element subsystem, are the shape and size of the element, and also the lines of singularity at Ro = m for real K, discussed previously. Elongated elements yield larger matrix values than an equilateral triangle of the same area, and the line singularities produce matrix values roughly two orders of magnitude larger than the average matrix value. The combination of large and small matrix elements can undermine the condition of the matrix. Those factors influencing the BI portion are the element size and the total number of boundary elements. Elements which are very small (< 0.01A) result in subsequent rows which are nearly equal and give the appearance of a singular matrix. Large dense systems are known to have a smaller Kc simply because of the number of matrix elements.

140 2.16 2.21 2.44 2.49 2.81 4.21 (a) 6.0.............. 4.0 2.0 i, -2.0 - -4.0 - FE-BI........ exact -6.0 2.0 2.5 3.0 3.5 4.0 4.5 5.0 k0ao (b) Figure 6.6: The mesh used for the frequency sweep of a conducting sphere is shown in (a) and the expanded region shows the values of ka at which the line singularity intersects the nodes. The frequency sweep in (b) demonstrates the inaccuracies associated with the line singularities

141 0.000 -2.00 mc -4.00 - -6.00 - -8.00 -10.0.... 2.00 2.50 3.00 3.50 4.00 4.50 5.00 ka Figure 6.7: A sphere coated with, =, = 1 - j5 is swept in frequency The issue of accuracy, which is closely tied to the issue of system condition, is examined as a function of frequency for two spheres. Each is comprised of a conductor of radius a with a mesh termination at ao = 1.02a. The second sphere contains dielectric material with c, = - 1 - j5 in the finite element region. As seen in Fig. 6.6, there are locations where the solution has a large error with respect to the exact result [27], plotted as a function of the outer radius normalized to the free space wavenumber ko. The vertical lines denote the regions of greatest change, and from (a), it is clear that the edge of an element is nearly parallel to this line. This problem, discussed previously in section 6.2.3, is due to the fact the matrix elements are becoming large and the sensitive cancellation required is not occurring accurately. When a lossy coating is introduced, the line-singularity is no longer present, and as shown in Fig. 6.7, the results follow quite well with the moment method data.

142 6.6.3 Elimination of the Internal Resonance Corruption of the Boundary Integral The corruption of the solution due to interior resonances of the boundary integral perseveres in the three-dimensional case. The approach taken in chapter 3 is also employed here. It is again shown that for the coated sphere test structure, the "complexification" of the wavenumber in conjunction with the scattered field FE-BI formulation successfully removes the resonances. Since the theory surrounding the method was previously presented, it will not be repeated here. Consider the conducting sphere of radius a enclosed in a boundary of radius ao = 1.02a and f, = 2 - j4. Shown in Fig. 6.8 is a as a function of kao of the outer boundary for axial incidence and observation (thus requiring the solution for the m = 1 mode). Clearly, the expected resonances indicated by the vertical lines are present and are evident in the solution for a = 1 by the spiked behavior, except at koaO c 10.6 at which the resonance is shifted due to numerical errors. These are eliminated by setting k = ako and the results for a = 0.005, a =.007 and a =.01 are presented in the same figure. As seen by these it is clear that the solution is relatively insensitive to a. 6.6.4 Scattering from Various Test Bodies'lo validate the code, bistatic scattering patterns are considered in the four various sample targets. In each case the radar cross section in (2.33) is computed and the results are compared to a moment method solution [25, 19]. Most of the cases considered involve axial incidence (along the z-axis) and in this case TEz and TMz no longer have meaning. Planes defined by the electric field vector and the z-axis form the E-plane, and the magnetic field vector and the z-axis form the H-plane.

143 -.0.,..FE-BI(- l =l) 13.0 -------- BI (abr-jO.005) ------ FE-BI (l-j.007) ------ FE-BI (a=4I-jO.O1) 9.0 7.0 5.0 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 k0aO Figure 6.8: The axial backscatter cross section is displayed as a function of normalized radius koao. The structure is a conducting sphere of radius a coated with a dielectric with Cr = 2 - j4 has an outer radius of ao = 1.02. Employing a complex wavenumber k = koa, the resonance behavior is removed for the complex values in comparison to the ca = 1 case

144 The patterns shown are computed in these planes. First we consider the scattering from a coated conducting sphere of radius 0.098m with a 0.02m coating with c, = 2-j4. Fig. 6.9a demonstrates the convergence of the solution as more Fourier harmonics are summed together, and shown are the partial sums up to m=10. Fig. 6.9b shows the final converged result of the (a) (plotted on a scale shifted by -60 degrees to coincide with the pattern due to axial incidence) along with the axial incident result (requiring only m=1). Only the E-plane patterns were computed. A comparison to a method of moments result demonstrates agreement. Next we consider the scattering from a perfectly conducting sphere of radius 1A as indicated in Fig. 6.10. Both E-plane and H-plane patterns are included and the comparison with moment method results is favorable. In Fig. 6.11 is shown the patterns for a 2A x 0.176A perfectly conducting ogive with a dielectric coating of thickness 0.05 with = 2 - j2, shown for an axially incident plane wave. Both E-plane (TEs) and H-plane (TMz) are shown along with data generated via the moment method. Next we consider the case of a coated sphere-capped cone frustum, the mesh of which is indicated in Fig. 6.12. The structure is 4A in length, the spheres are of radii 1A and A and the coating thickness A. Additionally, the material is dielectric with cr = 4 - j5. The bistatic cross section computed for this structure for an axially incident plane wave is indicated for TEz (or Eplane) and TM, (or H-Plane) incidence in Fig. 6.13. Reference data provided indicates a good agreement.

sl[nsaI poiqlapi uawuow ol uos!.eduwoc d pue linsal quapTau!.T!xe aql qIAi umoqs (aauap!ou! Iv!xv ol anp usla$ed aiq tizM ap3puioo o saasdap 09- Xq pajj!ts jajis uo pal -joId) saaa2ap 09 IV 01-0 SapOU o10J uolnlios p;aDJAUOuo aqSL'sasaIap 09 JO a)uap!Lou us Joj uo!QnUos 3valjoD aiq o3 uUI!2aAUOD SapoWU JO UOIVWu -uns ail smoS (e) aLIayds palvo3 ioj uialled fuIualmes Dilisiq atqJ:6-9 aIn~r (q) [2ap] aliuv 0081 OO1 0OZI 0'06 0'09 0'0 0'0..........,......,........,....,..... 0'0 001 o,01 u010 o WO~w0'0 LW. ~ ~ —~ul.009 --------- 0.0~ O'O'OZ 0'01Z 0'08 0I 1 OS I 0' 06 0'09 OI=tU,OZL=tu ~~~~3tLI~~~~~~~~~OO....'....'....'............ fog

146 25.0..... 20.0 20.0 - FE-BI (CE)....... FE-BI rM) 15.0 - X mom C._M) _ -5.0 0.0 7 -5.0...... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 0 [deg] Figure 6.10: The bistatic scattering pattern for a perfectly conducting sphere of radius 1A. Both E-plane (TE) and H-plane (TM) are shown along with data generated via the moment method

147 10.0....,..,..,..,.., FE-BI (TE). 0.0........ FE-BI TM),' o MoM (TE): 9 MoM (FTM) ~' -10.0 d - O -20.0 -30.0 - -40,0........,............... -40.0..... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 O [deg] Figure 6.11: The bistatic scattering pattern for a 2A x 0.176A perfectly conducting with a.0.05A coating with Ec = 2 - j2 for an axially incident plane wave. Both E-plane (TE) and H-plane (TM) are shown along with data generated via the moment method Figure 6.12: The mesh of the sphere-capped cone frustrum

148 30.0 FE-BI (TE). MoM (TE) 20.0 10.0 0.0 -10.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 0 [deg] (4=0) (a) 30.0. - FE-BI (TM) 20. 0 —MoM (TM) 20.0 10.0 - / 0.0 -10.0..... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 O [deg] (Q=90) (b) Figure 6.13: The bistatic scattering pattern is shown for a plane wave axially incident on the coated sphere-capped cone frustum. Part (a) shows the E-plane pattern, while (b) demonstrates the H-plane pattern

149 6.7 Conclusions We have presented the finite element - boundary element formulation for axially symmetric bodies. The storage requirement associated with the boundary integral was reduced by symmetry considerations. Several patterns were generated to demonstrate the validity of the method, and the model shortcomings were presented as well. Also, a method for eliminating the problem of the line singularity inherent in the CAP finite element formulation was given and was shown to give real matrices for lossless structures. It should be noted that frequency sweeps and single mode radiation problems make the most efficient use of the employed CG solver. The method is not well suited, however, to problems incorporating multiple excitations or those requiring many Fourier modes.

CHAPTER VII Conclusion 7.1 Summary In this thesis, the FE-BI technique (presented in general terms in chapter II) was developed for two-dimensional and axially symmetric structures. The twodimensional case was based on the moment-method version developed by Jin [9]. An FE-BI formulation for rectangular enclosures was developed in Chapter III and lead to simple boundary integrals some of which had convolutional form which was exploited for reducing the memory requirement. In chapter IV, circular and ogival enclosures were considered. For the circular boundary, the boundary integral was entirely convolutional in form and an O(N) storage requirement was thus achieved at all times. It was shown that circular enclosures are consequently preferred for storage efficiency if the structure's outer boundary does not substantially deviate from a circle. In chapter VI, the FE-BI formulation was developed for axially symmetric structures. The finite element method for this problem was based on the CAP equations for generating the FE matrix system, and a boundary integral equation was used for the boundary condition on an arbitrarily shaped mesh termination boundary. Consequently the boundary integrals were no longer convolutional and a storage reduction 150

151 was, instead, achieved by exploiting certain symmetry properties of the boundary integral subsystem. Also, a new method was presented for circumventing the singularity problem associated with the CAP equations in the finite element portion of the system. A method was developed in chapter V for eliminating the internal resonance corruption of the solution, a difficulty found with most boundary integral formulations. These resonances correspond to the eigenvalues of the integral operator, and may be also viewed as the cut-off frequencies of a resonator with conducting walls of the same shape as the mesh termination boundary. Because the eigenvalues become very closely spaced for electrically large structures, any scattering computations become unreliable unless the resonances are suppressed. In summary, specific contributions of this work include the development and implementation of the FE-BI for rectangular, circular and ogival boundaries as described in chapters III and IV in a manner taking advantage of storage reduction schemes. A new method was also presented for suppressing resonance corruptions existing in almost all implementations employing some form of a boundary integral equation over a closed surface or contour. This method involved the introduction of a complex wavenumber and was demonstrated for both 2-D and axially symmetric bodies. The implementation of the FE-BI method for axially symmetric structures as described in chapter VI was for the most part new and incorporated a scheme for treating the line singularity associated with the CAP equations. This resulted in real matrices for lossless scatterers, a property consistent with 2-D and 3-D implementations of the finite element method.

152 7.2 Future Work The FE-BI technique presented in this work remains a promising technique for computing the scattering from ever larger structures. However, as the structure increases in size, so does the condition number of the resulting system due in part to its size. To reduce the number of unknowns, higher order basis functions must be used. For electromagnetic scattering, bases formed with hierarchical functions [5] show promise. Instead of standard functions, in which each undetermined coefficient corresponds to an unknown field quantity, those associated with the hierarchical function are comprised of both field quantities and field errors. For instance, a quadratic hierarchical function for a one-dimensional element takes the form 2 &= N;(x)q$' + af(x) (7.1) j=1 where f(z) is an arbitrary quadratic function, and the coefficient a is proportion to the deviation of the linear approximation to the quadratic one. Clearly, as higher order terms are added to the sum in (7.1), the smaller the coefficients become. Furthermore, the quadratic function f(x) may be chosen to be approximately orthogonal to polynomials of different order. Thus, in a finite element implementation, the off diagonal elements may become very small increasing the diagonal dominance of the matrix and leading to better matrix conditioning. The above idea can be extended to include functions other than polynomials. The physical optics solution could be employed as the fundamental solution and linear (or higher order polynomials) functions could then be employed as correcting functions. That is, the unknown coefficients in the solution would then be the deviation of the physical optics solution from the actual one. This could consequently aid in reducing the number of unknowns, particularly if the original approximation was reasonably

153 good. Note that this technique is applicable for the discretization of the boundary integral and finite element equations.

BIBLIOGRAPHY 154

155 BIBLIOGRAPHY [1] T.K. Sarkar, E. Arvas and S.M. Rao, "Application of FFT and the conjugate gradient method for the solution of electromagnetic radiation from electrically large and small conducting bodies," IEEE Trans. Antennas Propagat., vol. AP30, pp. 409-418, May 1986. [2] T.J. Peters, "Computation of the Scattering by Planar and Non-planar Plates using a Conjugate Gradient FFT Method," Ph.D. dissertation, Radiation Laboratory, The University of Michigan, 1988. [3] K. Barkeshli, "Applications of the CGFFT method in Scattering and Radiation Including Simulations with Impedance Boundary Conditions," Ph.D. Dissertation, The University of Michigan, 1991. [4] A. Taflove and K.R. Umashankar, "The finite-difference time-domain method for numerical modelling of electromagnetic wave interactions," Electromagnetics, vol. 10, pp. 105-126, Jan-June 1990. [5] O.C. Zienkiewicz and K. Morgan, Finite Elements and Approzimation, New York: John Wiley and Sons, 1983. [6] K.K. Mei, "Unimoment method of solving antenna and scattering problems," IEEE Trans. Antennas Prop., vol AP-22, pp. 760-766, Nov. 1974. [7] P. Silvester.and M.S. Hsieh, "Finite-element solution of 2-dimensional exterior field problems," Proc. Inst. Elec. Eng., vol. 118, pp. 1743-1747, Dec. 1971. [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] 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 [10] A.H.J. Fleming, "A finite element method for composite scatterers,", Pier 2: Progress in Electromagnetics Research, ed. M.A. Morgan, Ch. 2. [11] C.A. Balanis, Advanced Engineering Electromagnetics, New York: John Wiley and Sons, 1989.

156 [12] Harrington, R.F., Time-Harmonic Electromagnetic Fields, New York: McGrawHill, 1961, p. 203. [13] M.A. Ricoy and J.L. Volakis, "Integral equations with reduced unknowns for the simulation of two-dimensional composite structures," IEEE Trans. Antennas Propagat., Vol. AP-37, pp. 362-372, 1989. [14] J.M. Jin and V.V. Liepa, "Simple moment method program for computing scattering from complex cylindrical obstacles," IEE Proceedings, Vol. 136, part H, No. 4, pp. 321-329, Aug. 1989. [15] W.D. Murphy, V. Rokhlin and M.S. Vassiliou, "Solving electromagnetic scattering problems at resonance frequencies," J. Appl. Phys., vol. 67, pp. 6061-6065, May 1990. [16] A.D. Yaghjian, "Augmented electric- and magnetic-field equation," Radio Science, vol. 16, pp. 987-1001, Nov. 1981. [17] R. Mittra and C.A. Klein, "Stability and convergence of moment method solutions," in Numerical and Asymptotic Techniques in Electromagnetics, ed. R. Mittra. New York: Hemisphere, 1987 Reprint. [18] M.B. Woodworth and A.D. Yaghjian, "Derivation, application and conjugate gradient solution of the dual-surface integral equations for three-dimensional, multi-wavelength perfect conductors," in PIER 5: Progress in Electromagnetics Research, ed. Tapan K. Sarkar, Ch. 4. [19] J.R. Mautz and R.F. Harrington, "H-field, E-field and combined-field solutions for conducting bodies of revolution," AEU, pp. 157-164, 1978. [20] F.X. Canning, "Singular value decomposition of integral equations of EM and applications to the cavity resonance problem," IEEE Trans. Antennas Propagat., vol. AP-37, Sep. 1989. [21] 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. [22] J.A. Stratton, Electromagnetic Theory, New York: McGraw-Hill, 1941. [23] M.A. Morgan, "Numerical computation of electromagnetic scattering by inhomogeneous dielectric bodies of revolution," Ph.D. Dissertation, June 1976. [24] R.V. Churchill and J.W. Brown, Complez Variables and Applications, New York: McGraw-Hill, Inc., 1984. [25] J.M. Putnam and L.N. Medgyesi-Mitschang, "Combined field integral equation formulation for axially inhomogeneous bodies of revolution," McDonnell Douglas Research Laboratories report QAOO3, Dec 1987.

157 [26] J.J. Bowman, T.B.A. Senior and P.L.E. Uslenghi, Electromagnetic and Acoustic Scattering by Simple Shapes, Amsterdam: North-Holland Publishing Co., 1969. [27] H.L. Thal, Jr., "Exact circuit analysis of spherical waves," IEEE Trans. Antennas Propagat., pp. 282-287, March 1978.

3 9015 0514 8043 THE UNIVERSITY OF MICHIGAN DATE DUE 7(.oV 4,,