RL-868 FINITE ELEMENT-BOUNDARY ELEMENT METHODS FOR ELECTROMAGNETIC SCATTERING Jianming Jin, Valdis V. Liepa, and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109 October 11, 1989 RL-868 = RL-868

ABSTRACT FINITE ELEMENT-BOUNDARY ELEMENT METHODS FOR ELECTROMAGNETIC SCATTERING 1 by Jianming Jin, Valdis V. Liepa, and John L. Volakis Recent advances in computer speed and storage have led to an increasing interest in developing new methodologies to satisfy a need for accurate and efficient numerical simulation of complex open region electromagnetic scattering problems. In this thesis, several hybrid techniques which employ finite element and boundary element methods are presented. A finite element-boundary element method is presented for computing the scattering by cylindrical objects. The method uses artificial boundaries that can follow the contour of the scatterer so that the discretization region is minimized. More im — portantly, the method results in a highly sparse or uniformly banded matrix which The material in this report is also the Ph.D. thesis submitted by Jianming Jin in partial fulfillment of the requirements for the Ph.D. degree in the University of Michigan. i

can be efficiently solved using special algorithms. Results are presented to demonstrate the accuracy of the method as well as its capability and versatility. The method is then extended to treat the scattering by coated wedges and halfplanes. In this case, the physical optics approximation is employed to approximate the surface fields away from the edge of the structure. Numerical results are presented to show the surface field and surface impedance behavior near the edge. A finite element-boundary element method is also developed to characterize the scattering and transmission properties of an inhomogeneously filled aperture in a thick conducting plane. Of particular interest in this formulation is the use of a fast Fourier transform to evaluate the boundary integrals and the use of a conjugate gradient method to solve the system of equations. As a result, the method is particularly efficient in dealing with large apertures. The use of isoparametric elements is subsequently presented for modeling arbitrarily shaped and curved geometries. This leads to an improved efficiency and accuracy not shared with traditional elements. Finally, formulations for several hybrid techniques which combine the finite element method with either a surface integral equation or an expansion of vector eigenfunctions are presented and discussed for three-dimensional scattering problems. 11

ACKNOWLEDGEMENTS The work described in this report has been partially supported by the National Science Foundation under Grant ECS-8319595, Mission Research Corporation, Santa Barbara, California under contract SC-5818-88-0001, and by the Radiation Laboratory and the Department of Electrical Engineering and Computer Science of the University of Michigan in the form of fellowships and graduate research assistantships. 111

TABLE OF CONTENTS ABSTRACT................................ i ACKNOWLEDGEMENTS.......................... iii LIST OF FIGURES................................ vii LIST OF TABLES.............................. xi CHAPTER I. INTRODUCTION............................ 1 1.1. Motivation 1.2. Background 1.3. Overview II. BASIC ELECTROMAGNETIC CONCEPTS............. 6 2.1. Maxwell's Equations 2.2. Wave Equations 2.3. Boundary Conditions 2.4. Radiation Conditions 2.5. Surface Integral Equations 2.6. Radar Cross Sections III. FINITE ELEMENT METHOD AND ITS APPLICATION TO TWODIMENSIONAL PROBLEMS...................... 14 3.1. Solution of Boundary-Value Problems 3.2. Basic Concepts of Finite Element Method 3.3. Finite Element Analysis of Two-Dimensional Problems IV. SCATTERING BY CYLINDERS.................... 33 4.1. General Review 4.2. Formulation 4.3. Numerical Results 4.4. Comparison with Other Methods v

4.5. Concluding Remarks V. SCATTERING BY COATED WEDGES AND HALF-PLANES... 67 5.1. Introduction 5.2. Formulation 5.3. Results and Discussion 5.4. Conclusion VI. SCATTERING BY INHOMOGENEOUSLY FILLED THICK APERTURES.................................... 84 6.1. Introduction 6.2. Theoretical Formulation 6.3. Numerical Discretization 6.4. CG-FFT Implementation 6.5. Numerical Results 6.5. Concluding Remarks VII. APPLICATION OF ISOPARAMETRIC ELEMENTS........ 121 7.1. Finite Element-Boundary Element Formulation 7.2. Isoparametric Elements 7.3. Conclusion VIII. FORMULATIONS FOR TIIREE-DIMENSIONAL SCATTERING. 1130 8.1. Problem Description 8.2. Finite Element Formulation 8.3. Hybrid Formulations IX. CONCLUSION AND FUTURE WORK................ 142 9.1. Conclusion 9.2. Future Work BIBLIOGRAPHY................................ 146 vi

LIST OF FIGURES Figure 4.1 Geometry of Wave Scattering by a Cylinder. (a) Illustration of Artificial Boundaries. (b) Example of Discretization............ 37 4.2 Illustration of a Line Segment on rB....................... 41 4.3 Error in Backscattered Far-Field Coefficient vs Distance between rA and................................... 48 4.4 Convergence of the FE-BE Solution for the Backscattered Field of a Coated Circular Cylinder for TM Case................... 50 4.5 Convergence of the FE-BE Solution for the Backscattered Field of a Coated Circular Cylinder for TE Case.................. 51 4.6 Convergence of the FE-BE Solution for the Backscattered Field of a Square Dielectric Cylinder....................... 53 4.7 Bistatic Scattering Pattern of a Coated Circular Cylinder. (a) Amplitude. (b) Phase............................. 55 4.8 Backscattering Pattern of a Coated Rectangular Cylinder. (a) TM Case. (b) TE Case............................. 56 4.9 Backscattering Patterns of a Bare, Single- and Double-Edge Coated Ogival Cylinder. (a) TM Case. (b) TE Case.............. 58 4.10 Bistatic Scattering Pattern of an Inhomogeneous Circular Cylinder. 59 4.11 Plane Wave Scattering by an Inhomogeneous Rectangular Cylinder. (a) Bistatic Scattering. (b) Backscattering................ 61 4.12 Plane Wave Scattering by an Inhomogeneous Triangular Cylinder. (a) Bistatic Scattering. (b) Backscattering.............. 62 5.1 Geometry of a Coated Perfectly Conducting Wedge........... 70 5.2 Geometry of a Coated Perfectly Conducting Plane.......... 74 vii

5.3 Surface Fields of a Half-Plane with the Upper Face Coated. (a) Electric Field at Upper Surface Sa. (b) Normal Derivative of Electric Field at Upper Surface S1. (c) Electric Field at S2. (d) Normal Derivative of Electric Field at S2. t = 0.1A, er = 2.4-jl.0, jr = 1.0, f = 135~.................................. 76 5.4 Surface Fields of a Coated Right-Angled Wedge. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.06A, e, =2.0 - j2.0, tPr = 2.0, - = 135~............... 78 5.5 Surface Fields of a Half-Plane with its Edge and Two Faces Coated. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.06A, r, = 2.0, Pr = 2.0, 4 = 135~....... 79 5.6 Surface Fields of a Half-Plane with its Two Faces Coated Inhomogeneously for Edge-On Incidence. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.1A, er = 2.0(1 + e-2), r = 1.0, = 180~................. 81 5.7 Surface Fields of a Half-Plane with Thick Coating on the Upper Face. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.15A, er = 4.0 - jl.0, Pfr = 1.0, q> = 90~.. 82 6.1 General Geometry of an Inhomogeneously Filled Thick Aperture.. 87 6.2 Equivalent Problems for That Shown in Figure 6.1. (a) For Region I. (b) For Region II. (c) For Region III.................... 88 6.3 Equivalent Magnetic Current at Normal Incidence and Backscattering RCS for a 1X Wide and 0.25A Deep Groove. (a) Magnitude. (b) Phase. (c) RCS. Solid Lines Correspond to r = 1 and fr = 1; Dashed Lines Correspond to er = 4 - jl and Pr = 1. Solid and Dashed Lines: This Formulation; Circles and Squares: Mode Matching Solution.............................. 103 6.4 Equivalent Magnetic Current at Normal Incidence and Transmission Coefficient for a 1.6A Wide and 0.25A Thick Slot. (a) Cr = 1, Pr = 1; Solid Line: Current on the Upper Aperture; Dashed Line: Current on the Lower Aperture. (b) er = 2.56(1 - jO.1), Pr = 1; Solid Line: Current on the Upper Aperture; Dashed Line: Current on the Lower Aperture. (c) Transmission Coefficient Corresponding to er = 1, Pr = 1 (Solid Line); er = 2.56(1 - j0.1), Pr = 1 (Shorter Dashed Line); and er = 2.56(1 -j0.1), Pr = 1.2 - j0.35 (Longer Dashed Line). Solid and Dashed Lines: This Formulation; Circles: Mode Matching Solution [Wang, 1984].................... 105 viii

6.5 Illustration of an Aperture Filled with Three Dielectric Layers and a Strip Grating. Top and Bottom Layers: Er = 2.56, [r, = 1.0, Thickness = 0.2 cm. Middle Layer: c, = 4.0,,r = 1.0, Thickness = 0.2 cm................................... 106 6.6 Field Distribution Corresponding to the Geometry in Figure 6.5 at 10 GHz and Normal Incidence. (a) Upper Aperture Field. (b) Field at the Center Line Containing the Strip Grating. (c) Lower Aperture Field. L1 = L2 = 186, N = 1309, Tolerance=0.01, 139 Iterations.. 108 6.7 Bistatic Scattering Patterns for the Geometry in Figure 6.5 at 10 GHz. (a) itnc = 90~ (Normal Incidence), 139 Iterations. (b) inc = 60~, 176 Iterations. (c) (pinc = 30~, 161 Iterations. L1 = L2 = 186, N = 1309, Tolerance=0.01...................... 109 6.8 Plots of the Transmission Coefficient (T) Corresponding to the Geometry in Figure 6.5. (a) T versus Angle of Incidence (p) at 10 GHz; L1 = L2 =186, N = 1309, Tolerance=0.01. (b) T versus Frequency (f) at Normal Incidence; 0-10 GHz: L1 = L2 = 124, N = 625, Tolerance=0.01; 10-20 GHz: L1 = L = 186, N = 1309, Tolerance=0.01; 20-30 GHz: L1 = L2 = 248, N = 1747, Tolerance=0.01.......... 110 6.9 Magnitude and Phase of Ml (Circles) and M2 (Squares) for a 1A Wide and 0.25A Thick Slot Filled with a Dielectric Having c, = 1-j. L = L2 = 20, N = 126. Solid and Dashed Lines: This Formulation; Circles and Squares: Data from [Aucland and Harrington, 1978].. 111 6.10 Illustration of an Aperture Filled with (a) a Single Layered Strip Grating and (b) a Double Layered Strip Grating. d Denotes the Thickness of the Associated Dielectric Layer or Conducting Strip.. 113 6.11 Magnitude of M1 (Solid Line) and M2 (Dashed Line) for the Geometry in Figure 6.10(a) at 30 GIIz and Normal Incidence. L1 = L2 = 252, N = 1771, Tolerance=0.005, 389 Iterations............ 114 6.12 Magnitude of M1 (Solid Line) and AM2 (Dashed Line) for the Geometry in Figure 6.10(b) at 9 GHz and Normal Incidence. L1 = L2 = 248, N = 1743, Tolerance=0.01, 159 Iterations............ 115 6.13 Bistatic Scattering Patterns for the Geometry in Figure 6.10(a) at Normal Incidence. (a) 3 GHz, L1 = L2 = 126, N = 508, Tolerance=0.005, 97 Iterations. (b) 10 GHz, L1 = L2 = 189, N = 1140, Tolerance=0.005, 118 Iterations. (c) 30 GHz, L1 = L2 = 252, N = 1771, Tolerance=0.005, 389 Iterations.............. 117 6.14 Bistatic Scattering Patterns for the Geometry in Figure 6.10(b) at 9 GHz. L1 = L2 = 248, N = 1743, Tolerance=0.01. (a) cpinc = 30~, 102 Iterations. (b) inc" = 60~, 201 Iterations. (c) cp'n = 90~, 159 Iterations.................................. 119 ix

6.15 Plots of the Transmission Coefficients versus Frequency at Normal Incidence. The Circles Represent Actual Computation Points. (a) For the Geometry in Figure 6.10(a). (b) For the Geometry in Figure 6.10(b)................................... 120 7.1 Two Examples for Breaking a Circular Region. (a) A Five Quadrilateral Element Model with 20 Nodes. (b) A 12 Quadrilateral Element Model with 45 Nodes........................... 124 7.2 A Quadratically Curved Segment in the xy-Plane Can Be Transformed into a Straight Segment Lying on the (-Axis.............. 125 7.3 A Quadrilateral Element with Quadratically Curved Sides in the xy-Plane Can Be Transformed into a Square in the hq-Plane..... 128 8.1 Geometry of Wave Scattering by a Permeable Body.......... 132 x

LIST OF TABLES Comparison of the Two Approaches.................. 44 Comparison of the Three Methods.................... 64 Table 4.1 4.2 xi

CHAPTER I INTRODUCTION This dissertation deals with the development of new numerical methodologies to analyze and evaluate open region electromagnetic wave scattering. 1.1. Motivation In the frequency domain the problem of open region electromagnetic scattering is usually treated using integral equation methods (see, for example, Mei and Van Bladel [1963], Richmond [1965, 1966], Harrington [1968], and Wu and Tsai [1977a]), partial differential equation methods (see, for example, Mason [1982]), and hybrid techniques which combine the partial differential equation methods with a surface integral equation or an eigenfunction expansion (see, for example, Chang and Mei [1976], Marin [1982], and Jeng and Chen [1984]). The integral equation methods have the advantage of a simple numerical implementation with a minimum discretization region, since their formulations incorporate the radiation condition by employing the free-space Green's function. However, they have the disadvantage of a rather difficult formulation when dealing with complex media and usually result in large full matrices whose treatment is usually time consuming and costly. 1

2 In contrast, the partial differential equation methods have the advantage of a simple formulation even in the case of complex media and result in a simple numerical implementation. In addition, they produce highly sparse matrices which can be efficiently solved using special algorithms [George and Liu, 1981]. Their major disadvantage, though, is the need of extending the discretization to the far field region in order to enforce the radiation condition and this usually leads to an extremely large number of unknowns.1 The hybrid techniques represent a combination of the above two methods. They eliminate their disadvantages while maintaining their advantages. As a result, the hybrid techniques are more powerful especially for treating large and complex scatterers. However, the formulations and, especially, the numerical implementation of the hybrid techniques usually require more effort than for the integral equation and partial differential equation methods. Therefore, there exists a need to improve the existing hybrid techniques and to develop new ones that are more efficient for practical application. This is the goal of this dissertation. 1.2. Background There are basically two kinds of partial differential equation methods: finite element and finite difference methods, and both of them can be employed in the hybrid techniques. The hybrid techniques presented in this dissertation employ the finite element method and combine it with boundary element method. Therefore, they are referred to as finite element-boundary element methods. The following is a brief review of the background of the finite element-boundary element methods in electro1 Currently, much attention is aimed at employing more complicated radiation conditions or absorbing boundary conditions to reduce the discretization region (see, for example, McCartin et al. [1988], Mittra [1988], and Wilton and Richards [1988]).

3 magnetics, while a detailed review for each problem considered in this dissertation is given in the related chapter. The basic idea of the finite element-boundary element methods was first developed in mechanical engineering and introduced later in electromagnetics by Silvester and Hsieh [1971] and McDonald and Wexler [1972] when they were considering the finite element solution of exterior or unbounded field problems. The idea is to introduce an artificial boundary enclosing the scatterer, and express its exterior field by either an expansion of eigenfunctions or a boundary integral equation involving the freespace Green's function. Such expressions are imposed on the variational equation for the interior field as boundary constraints. The resultant variational equation is then solved via the finite element method. This idea has been further developed and improved upon in order to solve two-dimensional antenna problems [Washisu et al., 1979; Orikasa et al., 1983] and scattering problems [Marin, 1982; Jeng and Chen, 1984; Guo, 1985]. These methods, however, were not well-known in the mid 1980's since their capability and computational efficiency were not obvious. A better known method is the one called unimoment method and was developed by Mei and his students [Mei, 1974; Chang and Mei, 1976; Morgan and Mei, 1979]. In this method, the finite element method is used to generate a set of trial functions inside a mathematical circle encompassing the scatterer. The field outside the circle is expanded in terms of Hankel functions and then coupled with interior solutions on the circle. Recent advances in computer speed and storage and an increasing need to study more complicated scattering problems have led to a growing interest in using the finite element-boundary element methods. As a result, an increasing amount of research on the subject, including that of the Radiation Laboratory of The University of Michigan, has been reported in the past two years.

4 1.3. Overview This dissertation includes the part of my research conducted during the past four years on developing new finite element-boundary element methods for a number of scattering problems. The dissertation is organized in such a way that each chapter is self-contained and independent, and can be read without intensive reference to the other chapters. Chapter 2 presents a brief review of some basic concepts of electromagnetic theory and also contains pertinent equations that are used in subsequent chapters. Chapter 3 first briefly reviews several methods for solving boundary-value problems, then presents some basic concepts of the finite element method, and finally describes its application to two-dimensional problems. The chapter is intended for a reader who is not very familiar with the finite element method. Thus, in the later chapters we can concentrate our attention on the principles and formulations of the hybrid techniques, rather than the finite element formulation. Chapter 4 deals with the scattering by cylinders. A general review of various methods is given and the major shortcomings that exist in each are pointed out. A new formulation combining the finite element method with a surface integral equation is then presented which overcomes those shortcomings, and numerical results are given to verify the method as well as to demonstrate its capability. Finally, the method is compared with other commonly used numerical methods and its advantages are clearly shown. The work contained in this chapter was published in 1988 [Jin and Liepa, 1988a and 1988b]. Chapter 5 describes a technique for the scattering by a coated wedge and halfplane. This technique is an extension of that presented in Chapter 4 by incorporating the physical optics approximation for the fields far away from the edge. A brief review for the problem is given and then followed by description of the formulation.

5 Numerical results are subsequently presented to show the surface field and surface impedance behavior around the edge. The work reported here was published in 1989 [Jin and Liepa, 1989a]. Chapter 6 presents a very efficient technique for the electromagnetic characterization of the scattering and transmission properties of an inhomogeneously filled aperture in a thick conducting plane. The new technique again combines the finite element method with boundary integral equations as was done in Chapter 4. Of particular interest is the use of the conjugate gradient method [Hestenes and Steifel, 1952] and fast Fourier transform for the solution of the system. The technique is successfully used to solve some problems which are formidable when other techniques are used, and this is shown by numerical results. The work in this chapter was published in 1989 [Jin and Volakis, 1989c]. Chapter 7 describes the application of isoparametric elements in the finite elementboundary element methods. The purpose for using the isoparametric elements is to improve the computational efficiency and accuracy of the previously developed methods which usually employ the linear triangular elements. Chapter 8 presents several hybrid formulations involving the finite element method for the three-dimensional scattering. Numerical implementation of these formulations is a subject for further investigation. Chapter 9 concludes this dissertation with a brief summary and recommendations for future work.

CHAPTER II BASIC ELECTROMAGNETIC CONCEPTS The analysis of electromagnetic scattering and radiation is actually a problem of solving Maxwell's equations subject to given boundary conditions. In this chapter, we will briefly review some basic concepts and equations of electromagnetic theory, which will be used later on in our formulations. It should be noted that our emphasis is laid upon their mathematical expressions rather than their physical interpretations. 2.1. Maxwell's Equations Maxwell's equations are a set of fundamental equations governing all macroscopic electromagnetic phenomena. In differential form they are written as VxE — at VxH=J+ vx7=7++ a at V.D=p V.B=O (Faraday's law) (Maxwell-Ampere's law) (Gauss' law) (Gauss' law- magnetic) (2.1) (2.2) (2.3) (2.4) where E = electric field intensity (volts/meter) 6

7 D = electric flux density (coulombs/meter2) H = magnetic field intensity (amperes/meter) B = magnetic flux density (webers/meter2) J = electric current density (amperes/meter2) p = electric charge density (coulombs/meter3) Another fundamental equation, which is known as the equation of continuity, can be written as V.=- p (2.5) which states the conservation of charge. Among the above five equations, only three are independent, and thus called independent equations. Either the first three equations, (2.1)-(2.3), or the first two equations, (2.1) and (2.2), with (2.5) can be chosen as such independent equations. The other two equations, (2.4) and (2.5) or (2.4) and (2.3), can be derived from the independent equations, and thus called auxiliary or dependent equations. The three independent equations are in indefinite form, since the number of equations is less than the number of the unknowns. Maxwell's equations become definite when the constitutive relations between the field quantities are specified. The constitutive relations describe the macroscopic properties of the medium being considered. For a simple medium, they are D = eE (2.6) B = H (2.7) J = aE (2.8) where the constitutive parameters e, to and o denote, respectively, the permittivity (farads/meter), permeability (henrys/meter) and conductivity (mhos/meter) of the medium. These parameters are tensors for anisotropic media and scalars for

8 isotropic media. For inhomogeneous media, they are functions of position; while for homogeneous media they are not. When the field quantities in Maxwell's equations are harmonically oscillating functions with a single frequency, using complex phasor notations, (2.]1)-(2.5) can be written in a simplified form V x E=-jwsH (2.9) V x H = +jjeE (2.10) V. (eE) = p (2.11) V. (1H) = 0 (2.12) V J= -jwp (2.13) where the time convention e3wt is used and suppressed, and w is angular frequency. 2.2. Wave Equations To solve Maxwell's equations, one usually first converts the first-order differential equations involving two vector variables into the second-order differential equations involving only one vector variable. By eliminating H or E in (2.9) and (2.10), one obtains V x -V x E-w 2E = -jwJ (2.14) V x -V x H - o2iH = V x ( ) (2.1.5) where Ji is an impressed or source current, and sc (=e - joa/w) represents a combination of induced current and displacement current; however, for simplicity we next use e to denote ec. Equations (2.14) and (2.15) are called inhomogeneous vector wave equations.

9 For two-dimensional electromagnetic fields (the fields have no variation with respect to one Cartesian coordinate, say the z-coordinate), it can be shown that the z-components of (2.14) and (2.15) become (x —x+ y-O + k e Ez = j koZoJ (2.16) (a a + aI + Hz =- (J) + ( (2.17) \9x er ax ayr ay e ax u - 9y \a e where er (= e/Eo) and p,r (= p/,lo) denote the relative permittivity and relative permeability, respectively, and they are assumed here to be complex scalar functions of position; ko (= wV/o0j) is the wavenumber in free-space; Zo (= \fo/o) is the intrinsic impedance of free-space; and the electric constant c0 (= 8.854 x 10-12 farads/meter) and the magnetic constant lto (= 4w7 x 10-7 henrys/meter) are the permittivity and permeability of free-space. Equations (2.16) and (2.17) are called inhomogeneous scalar wave equations. 2.3. Boundary Conditions Solving the inhomogeneous vector or scalar wave equations given above in a domain of interest, one may obtain many solutions; however, only one of then is the real solution to the problem. To find out this real solution, one should know the boundary conditions associated with the domain. In other words, a complete description of an electromagnetic problem should include complete information about both differential equations and boundary conditions. This section presents some boundary conditions which apply to many practical problems. Boundary Conditions at the Interface between Two Media. At the interface between two media, say medium 1 and medium 2, the boundary conditions can be mathematically expressed as n x (E1 -E2) = 0 (2.18)

10 n. (D -D2) =0 (2.19) for electric fields, and similarly n x (H -H2) = 0 (2.20) n (B - B2) = 0 (2.21) for magnetic fields, where n is the unit vector normal to the interface, pointing from medium 2 into medium 1. It is assumed in (2.19) and (2.20) that neither surface currents nor surface charges exist at the interface. It is also understood that among the four boundary conditions only two of them are independent: one from (2.18) and (2.21), and the other from (2.19) and (2.20). Boundary Conditions at the Surfaces of Perfect Conductors. This is a special case when one of the media, say medium 2, becomes a perfect conductor. Since a perfect conductor cannot sustain a field inside, (2.18) reduces to n x E 0 (2.22) and (2.21) reduces to n - B=0 (2.23) Note that in this case the boundary can support a surface current (J,, = n x H) and surface charge (p, = i D). 2.4. Radiation Conditions When the domain of interest is unbounded (its boundary is at infinity), a condition can be specified at infinity. Assuming that the sources and scatterers are located within a finite distance from an origin, the fields are required to satisfy lim R V x +jkR x =0 (2.24) L-+o -f/' L \ / \ / J

11 where R = /x2 + y2 + Z2. This is referred to as the radiation condition for general three-dimensional fields. For two-dimensional fields, the radiation condition becomes lim E + jk 0 (2.25) where r = V/x2 + y2. 2.5. Surface Integral Equations This section presents the general solutions, in the form of integral equations, to the electromagnetic field problems within a homogeneous isotropic medium. Consider an electric current J radiating in a domain V bounded by the surface S. The electric field inside can be written as [Tai, 1971] E(R) = -j II| J )Go(,R') dV- { V x Go(R, ) [n' x E(R')] - j Go(RR') [n x H(R')]} dS' (2.26) where n' is the outward unit vector normal to the surface S, and Go(R,R') is the dyadic Green's function defined by Go(R, = (I + -VV) Go(R) with e-jklR-~R' Go(R, R) = - IR4' \R-R'I and I = xx + yy + zz Similarly, one can find the magnetic field as H(R) = Go(R R) V' x J(R )dV' - {V x Go(IRR') * [ x H(R7')] + ji o(R,') * [ x E(R)] } dS (2.27)

12 If we consider the domain V to be of infinite extent, the outer surface, denoted by Soo, recedes to infinity. Because of the radiation condition described by (2.24), the surface integral over SO vanishes, and thus, (2.26) and (2.27) are still valid. In scattering and radiation problems, S usually represents the surfaces of obstacles or scatterers. The first term on the right hand side of (2.26) and (2.27) represents the field radiated by the electric current in the absence of any scatterers, and thus can be referred to as the incident field. The second term comes from the integrals over the surfaces of scatterers, and thus represents the scattered field. Equations (2.26) and (2.27) can be uniformly written as F(R) = F'(R) -J { V x Go(R,R') [ni x F(R')] + Go(R,R') [n x x F(R)] dS (2.28) where F'(R) denotes the incident electric or magnetic field, and F denotes the corresponding total electric or magnetic field. Equation (2.28) is referred to as the surface integral equation for the three-dimensional case. For the two-dimensional case, in an unbounded domain with homogeneous isotropic medium containing scatterers, the solution to (2.16) and (2.17) is found to be F(r) = F'(r) - [Go(, ') a( - F(') G ]dC (2.29) where F = Ez for E-polarization, F = Hz for H-polarization, and C is the boundary of all scatterers. Also, Go(r, r') is the two-dimensional Green's function given by Go(r, ') =- 2)(kIr-' r) and F(r) is the incident field given by Ez (r) = -jW J GO(r, r')Jz(r')ds' for E-polarization, and Hz) Go( i aJ() a('y ds'

13 for H-polarization. Equation (2.29) is referred to as the two-dimensional surface integral equation or boundary integral equation. 2.6. Radar Cross Sections The radar cross section (RCS) is a quantity characterizing the scattering of an obstacle. It is defined for plane wave incidence, and is also known as the scattering cross section. In the three-dimensional case, the scattering cross section is defined by a(O, b)= lim 41rR2 (2.30) R —oo JF'12 where FS denotes the scattered field at the observation point (R, 0, ~). In the twodimensional case, the scattering cross section per unit length is defined by c(=) = lim 27r i(2.31) r —oo FP 12 where Fs is the scattered field at the observation point (r, q).

CHAPTER III FINITE ELEMENT METHOD AND ITS APPLICATION TO TWO-DIMENSIONAL PROBLEMS The finite element method (FEM) is a computer-aided mathematical technique for obtaining approximate numerical solutions to boundary-value problems of mathematical physics. The method has a history of about forty years. It was first proposed in the 1940's and its use began in the 1950's for aircraft design. Thereafter, the FEM has been developed and applied extensively to problems of structural analysis and less extensively, but increasingly, to problems in other fields. Today, the FEM has become recognized as a general method of wide applicability to engineering and mathematical problems. Many books have been written on the subject, of which Zienkiewicz's [1977] is probably one of the most popular reference texts. In this chapter, we will first briefly review several methods for solving boundary-value problems. Some basic concepts of the FEM will then be presented and followed by a FEM analysis for two-dimensional field problems. 3.1. Solution of Boundary-Value Problems Boundary-value problems have long been a major topic in mathematical physics and nowdays one has many methods available for finding their solutions,. A boundaryvalue problem can be mathematically defined by a governing differential equation in 14

15 a domain Q LQb = f (3.1) together with the associated boundary conditions on the boundary F of the domain. In (3.1), L is usually a differential or integral operator, f is the excitation or forcing function and q5 the unknown quantity. It is, of course, desirable to solve (3.1) analytically, whenever possible. However, this is generally the exception since an analytical solution can only be attained in a few cases for given F. To overcome this difficulty, numerical methods have been and continued to be developed for obtaining an approximate solution by means of computing q at discrete points. These methods can be categorized into three groups: the finite difference method, the variational method, and the method of weighted residuals. 3.1.1. Finite Difference Method The finite difference method approximates the differential operator and is best suited for problems with regular boundary domains and homogeneous media. Using this method, the differential operator in (3.1) is represented by a difference operator and thus (3.1) is transformed into a set of difference equations. This method has been mostly used for solving field problems with regular boundaries, such as rectangular waveguides in electromagnetics. However, it is not easily adapted to problems involving curved or irregular boundaries because of the difficulty of enforcing boundary conditions. The method is also not as effective if the problem involves different media since the boundary condition at the interface between media requires special treatment. As a result, it is difficult to develop general purpose computer programs using the finite difference method.

16 3.1.2. Variational Method The variational method formulates the boundary-value problems in terms of variational expressions, referred to as functionals, with (3.1) as their Euler equations under the given boundary conditions. The approximate solutions are then obtained by minimizing the functional with respect to its variables. Specifically, consider (3.1). If the operator L is self-adjoint, i.e. < L+1, 0 >=<,,Lo > (3.2) and positive-definite (>0 <#0 < L+, > (3.3) then the solution to (3.1) can be obtained by the minimization of the functional [Mikhlin, 1964] F(q') =< L+', ' > - < ', f > - < f, ' > (3.4) with respect to q'. The inner product, denoted by the angular bracket, is defined by < <,i >=- / *dQ (3.5) where the asterisk denotes the complex conjugate. Once the functional is found, the solution can be obtained by the procedure described below. For simplicity, let us assume that the problem is real-valued and suppose that q' in (3.4) can be approximated by the expression N E' = cv- CT = vTC (3.6) i=1 where vi are the chosen expansion functions and ci are constant coefficients to be determined. Also, the boldface letters denote column vectors and the superscript T denotes the transpose of the vector. Substituting (3.6) into (3.4), we obtain F('Y) = cT vLvTdQc - 2cT vfdQ (3.7)

17 To minimize F(q') we let its derivatives with respect to ci vanish. This yields the matrix equation [S]c = b (3.8) whose elements in [S] are given by S = I vLvjdf (3.9) and those in b by bi = fdQ (3.10) As a consequence of the self-adjoint property of the operator, the matrix [S] is symmetric. Solving (3.8), one obtains an approximate solution for (3.1). The variational method described above has been the basis for many finite element formulations since the beginning of the FEM history. It does not apply, though, to cases where the operator L is not self-adjoint. In that case, a functional can be constructed leading to an adjoint operator before applying the variational method procedure [McDonald and Wexler, 1980]. Finding and minimizing a functional is, however, not the only way to formulate the finite element equations. The weighted residual methods can provide an alternative approach. 3.1.3. Weighted Residual Methods Based on its definition, the weighted residual method weights a residual. Assume that q' is an approximate solution to (3.1). Substitution of q' for k in (3.1) would then result in a non-zero residual R= Lq'-f (3.11) The best approximation for q will be one that reduces the residual, R, to the least value at all points of Q. The weighted residual methods enforce the condition RwidQ=0 (3.12)

18 where w, is a chosen weighting function. Different choices of the weighting functions can result in different formulations as discussed next. Collocation Method. This is also known as the point matching method. Impulse functions are selected as the weighting functions (wi = 1 at some point i and zero everywhere else). This is equivalent to satisfying (3.1) at specific points. Subdomain Method. The weighting functions here are set equal to unity over a specific subdomain and zero elsewhere. This is equivalent to making the integral of the residual vanish over a number of subdomains. Least Squares Method. The least squares method uses the residual as the weighting function and minimizes a new error term defined by e = j R2dQ (3.13) The minimization is with respect to the unknown coefficients in the approximate solution. Galerkin's Method. In this method, the weighting functions are selected to be the same as those used for the expansion of the approximate solution. This usually leads to the most accurate solution and is, therefore, a popular approach in developing the finite element equations. To more explicitly illustrate the method, let us assume that the solution is represented as in (3.6). The weighting functions are then selected as wi = vi i = 1,2, 3,..., N (3.14) and (3.12) becomes j(vLvTc- vif)dQ = 0 i = 1,2,3,-.,N (3.15) This again leads to a matrix system such as that given in (3.8), although now the matrix [S] is not necessarily symmetric unless the operator L is self-adjoint. In that case, the variational and Galerkin's methods result in identical equations.

19 3.2. Basic Concepts of Finite Element Method As stated at the beginning of this chapter, the finite element method (FEM) is a numerical procedure for obtaining solutions to boundary-value problems of mathematical physics. The principle of the method is to replace a continuous domain by a number of subdomains, which are usually referred to as elements whose behavior is modelled by interpolation functions containing a few unknown values at the nodes of the elements. Thus, the original boundary-value problem with an infinite number of degrees of freedom is converted into a problem with a finite number of degrees of freedom, or in other words, the behavior of the whole system is approximately represented by a finite number of nodal values. Then, a set of algebraic equations or a system of equations is obtained by assembling all elements and the solution of the boundary-value problem is achieved by solving the system. The most important feature of the system of equations is its sparse form due to the highly localized interaction of nodes. In summary, the usual steps in a finite element solution of a boundary-value problem are as follows: (i) discretization of the domain, (ii) derivation of the element equations, (iii) assembly of the elements, and (iv) solution of the system of equations. Below we describe each step in some detail. 3.2.1. Discretization of Domains In a finite element analysis, the discretization of the domain is the first and most important step. The manner in which the domain is discretized will affect the computer storage requirement, the processing time, as well as accuracy of the numerical results. The discretization involves these essential tasks: (i) subdivision of the domain into elements, (ii) numbering of the elements and nodes, and (iii) selection of the expansion or interpolation function to model the behavior of the field within each elements.

20 Domain Subdivision. The boundary-value problem represented by (3.1) may be solved in a one-dimensional, two-dimensional or three-dimensional space. For a onedimensional domain which is actually a straight or curved line, the elements are shorter line segments interconnected to form the orginal line. For a two-dimensional domain, the elements are usually triangular and rectangular surface areas. The rectangular element is, of course, best suited for discretizing rectangular regions while the triangular one can be used for irregular regions. In a three-dimensional solution, the domain may be subdivided into tetrahedra, rectangular prisms or triangular prisms, among which the tetrahedron is the simplest and best suited for arbitrary volume domains. We note that the line segments, triangles and tetrahedra are the basic one-, two- and three-dimensional elements which model the curved lines or surfaces by straight line segments or planar patches. A better modelling can be achieved by introducing isoparametric elements, which will be discussed in Chapter 7. Node and Element Numbering. In a finite element analysis, it is necessary to label the elements and associated nodes for implementation purpose. A complete description of a node contains its coordinate values, local number and global number. The local number of the node indicates its position in the element, whereas the global number gives its position in the entire system. While specifying the coordinate values is a rather straightforward job, numbering the nodes and elements requires some strategy. Since the maximum difference between the global numbers of two adjacent nodes in an element determines the bandwidth of the system matrix, this bandwidth can then be minimized by properly numbering the nodes. Thus, if a banded matrix method is used to solve the final system of equations, the computer storage and processing cost can be reduced significantly. However, in the case that the bandwidth minimization is unnecessary, the numbering scheme can be arbitrary and is usually chosen to simplify the programming.

21 Interpolation Function Selection. The interpolation function provides an approximation of the unknown solution within an element. It is usually selected to be a polynomial of first (linear), second (quadratic) or higher order. High-order functions, although more accurate, usually result in a more complicated formulation. Hence, a simple and basic linear interpolation is still widely used. Once the order of the polynomial is selected, one can derive an expression for the unknown solution in n le = Nj (3.16) j=l where n is the number of nodes in the element, q$ the field value at node j, and Ne the expansion function which is also known as basis function or shape function. 3.2.2. Derivation of Element Equations The second step in a finite element analysis is to derive the element equation. Two popular approaches exist for this derivation. One employs the variational principle and the other uses Galerkin's method. The least squares method can also be employed, but will not be discussed here. Derivation Using Variational Method. The element equation can be derived by applying the variational method, described in Section 3.1.2. The functional of the system can be expressed as a summation of integrations over each element. Thus, for a real-valued problem one has M F(q') = Fe (l) (3.17) e=l where M is the total number of the elements and F ( ) = | cL'edQ - 2 fq'edQ (3.18) J~e "

22 Substituting (3.16) into (3.18) and differentiating Fe with respect to the nodal field values, we obtain I1F e n 2 T=E Ne LNedQN - fNdfQ i= 1,2,3,..-,n (3.19) Z j=l * e ~ This can be written in matrix form as 2 K e - H} (3.20) where [Ke] is an n x n matrix and {He} an n x 1 column vector with their elements given by K = Ie NLNJedQ (3.21) and H =Ie f NifdQ (3.22) We note that the element matrix [Ke] is symmetric since L is self-adjoint. Derivation Using Galerkin's Method. Galerkin's method is an alternative for the derivation of the element equation. When using this method, the weighted residual equation for the eth element is Qe ijN (LO - f)dQ =1,2,3,.., n (3.23) Substituting (3.16) into (3.23) then yields n NL N dQ I fN= i=1,2,3, *,n (3.24) which can again be written in a matrix form as [K ]{} - {H} = 0 (3.25) In this the matrix elements IKs and Hi are of the same form as (3.21) and (3.22), respectively. Since the operator L is not required to be self-adjoint, here the element matrix [Ke] is not necessarily symmetric. However, when the operator L is selfadjoint, the element matrices resulting from the variational and Galerkin's methods are identical.

23 3.2.3. Assembly of the System of Equations Once the element equation or subsystem has been derived, the entire system of equations can be obtained by assembling each subsystem in conjunction with the boundary conditions. Either assembling (3.20) for all elements and imposing the stationary requirement or assembling (3.25) yields M Ei([K[]{e} - {He}) = 0 (3.26) e=l By changing the local number of each node in (3.26) to the corresponding global number, (3.26) can be written more compactly as [K]{q} - {H} = 0 (3.27) where [K] is an N x N matrix, {H} an N x 1 known vector, {)} an N x 1 unknown vector whose elements are the nodal field values, and N denotes the total number of nodes in the system. Before solving (3.27), one has to apply the boundary conditions associated with the problem. There are two kinds of boundary conditions which are often encountered: the homogeneous Dirichlet and Neumann boundary conditions. Since the latter is usually satisfied automatically as a natural boundary condition, one only needs to impose the Dirichlet boundary condition by setting the fields at the appropriate nodes equal to zero. The resultant system is then solved via matrix inversion, LU decomposition, or iterative approach. 3.2.4. Solution of the System of Equations Solving the system of equations is the final step in a finite element analysis. The resultant system has one of the two following forms: ([A] - A[B]) {X} = 0 (3.28)

24 or [A] {X}= {Y} (3.29) Equation (3.28) is of the eigenvalue type and in electromagnetics it is usually associated with waveguide and cavity problems. In that case, A is the eigenvalue, {X} is the corresponding eigenvector, both of which are unknowns. Equation (3.29) is of the deterministic type and in electromagnetics it usually associated with scattering and radiation problems. Since we are interested in the latter in this thesis, only the solution to (3.29) will be discussed here. To choose an algorithm for the solution of (3.29), one needs to consider the properties of the matrix [A]. Theoretically, there are many algorithms available to solve this matrix equation. However, since [A] here is usually large but sparse and banded provided the nodes are properly numbered, it is instructive to look at algorithms that exploit these properties. We describe three of them here. One solution approach is the banded matrix method with a variable bandwidth storage scheme. Using this method, the matrix [A] is stored as a vector rather than as a two-dimensional array. Only the elements from the first non-zero to the last non-zero (or to the diagonal element for a symmetric matrix) in each row are stored. The solution proceeds by first decomposing the matrix into an upper triangular form using Gaussian elimination (however, if the matrix is symmetric and positive-definite, an LDLT decomposition will be more efficient). The solution is then obtained by back substitution. Another solution method for sparse systems is the frontal method which was developed in the early seventies [Irons, 1970]. This method applies Gaussian elimination during the assembling process of the matrix [A]. Such progressing continues until every element and every node has been processed. The final solution is then obtained by back substitution. Because the assembly and elimination are intertwined, the matrix [A] is never formed explicitly. As a result, the required computer storage

25 is reduced to minimum. Also, the operations on the zero elements of [A] are avoided. The Lanczos method is a third solution method which is also referred to as the "bi-conjugate gradient method" [Sewell, 1985]. If the matrix [A] is symmetric, the method is actually the usual conjugate gradient mthod.d is actuallt to the abovusual conjugate two direct methods involving Gaussian elimination, this one solves the system iteratively. The assembly p lace in each iteration and the calculation of the matrix products [A]{P} are performed by executing the summation ZEMl[te]{Pe} where [Ae] is the element matrix and {P} denotes a vector. Thus, like the frontal method, the Lanczos method does not actually form the matrix [A] and, as a result, it avoids the storage of the matric [A] as well as the operations on the zero elements which comprise most of the matrix. Among the three methods, the banded matrix method is the most preferable for matrices with dimensions up to a thousand or so. However, it requires much larger storage than the other two methods since the elements within the band of [A] must be stored explicitly. The frontal and Lanczos methods are capable of solving very large systems with many thousands of unknowns, but at the expanse of computing time. 3.3. Finite Element Analysis of Two-Dimensional Problems In the previous section we briefly introduced the basic concepts of the FEM and in this section we describe the application of the FEM to two-dimensional problems. The pertinent differential equation is of the form - ( [( )as(x)] - a ( )a(xY) + (x, y)(x, y) = f(x, y) (3.30) where +(x, y) is the unknown field, a(x, y), t/(x, y) and -y(x, y) are known parameters associated with the physical properties of the geometry, and f(x, y) is the source or excitation function. The ordinary two-dimensional Laplace equation, Poisson

26 equation and Helmholtz equation are only special forms of (3.30). To solve (3.30) via the FEM we will employ linear triangular elements and the system is derived using the variational and Galerkin's methods. The results obtained in this section will be frequently used in subsequent chapters. 3.3.1. Shape Functions for Linear Triangular Elements To apply the FEM to (3.30), one first divides the domain of interest, denoted by S, into a number of two-dimensional elements. If the linear triangular elements are used, the field within each element is approximated as 'le(x, y) = a + bx + cy (3.31) where a, b and c are constant coefficients to be determined and e is the element number. For a linear triangular element, there are three nodes located at the vertices of the triangle. Assume that the nodes are numbered counterclockwise by numerals 1, 2 and 3 with the corresponding field values denoted by $e^, q and O|, respectively. Solving for the constant coefficients a, b and c in terms of qie yields 3,e(, y) = EN(. (xY)i (3.32) i=1 where Ne(x, y) are referred to as the shape functions given by Ni(x,y) = 2 (ai + bix + ciy) i = 1,2,3 (3.33) in which al = x2y3- Y2x3; bl= Y2 - y3; C = x3- x2 a2 = x3y1 - y3x1; b2 = Y3 - Y1; c2 = x1 - x3 a3 = xly2 - ylx2; b3 = Y1 - Y2; C3 = X2- X1

27 and 1 X1 yl = 1 X2 Y2 = area of the eth element 1 X3 y3 It can be easily shown that the shape functions have the property 1 i-j Nie(xj, y)= -Sj (3.34) and as a result at the nodes the field expression (3.32) reduces to the corresponding nodal fields q7. 3.3.2. Derivation of Finite Element Equations Using Variational Method As stated earlier, the variational method is a traditional technique for deriving the finite element equations. Let us denote the contour enclosing S as C and assume that the field on a portion of C, say C1, satisfies the Dirichlet boundary condition Olat c1 = 0 (3.35) and the field on the remaining portion, say C2, satisfies the Neumann boundary condition (a + a ) n =0 (3.36) 0 0 9Jy at c2 where fn is the outward unit normal to C. It can be shown that the solution to the differential equation (3.30) subject to the above boundary conditions is equivalent to the variational problem ( SF() - 0 (3.37) I lat c, = 0 where F(Js') / a, + 0 + 2 ds - J fq'ds (3.38)

28 The functional F(q') here is obtained by first substituting (3.30) into (3.4), then invoking the divergence theorem, and enforcing the Neumann boundary condition (3.36). Next assume that the region S is divided into M triangular elements of surface area Se (e = 1,2,3,- *, M). The integral over S in (3.38) can then be expressed as a summation of M sub-integrals over each element. This in turn allows the subdivision of the functional F in (3.38) into M sub-functionals Fe which are identical to F except that S is replaced by Se. Introducing the expression (3.32) for /e and differentiating each Fe with respect to Xl yields { F } =[ + [AB]e } + [B{ He} (39) where {P FFe OF= OFe1T e {-}e =5~ [ ae; {>e} = [$ 2 >3] The elements of the matrices [Ae] and [Be] are given by aA = lax j + y ay ) dxdy (3.40) 13 e a\ x 9x ay ay BeJ = j e NJNdxdy i,j = 1,2,3 (3.41) and those of the vector {He} by Hi = Jj f Nedxdy i = 1,2,3 (3.42) It is noted that both [Ae] and [Be] are symmetric here. Assuming now that the coefficients a, P/, - and the source f are constant within each element and equal to ae, 3e, 7e and fe, respectively, (3.40)-(3.42) can then be analytically integrated. A basic result in this process is J js (N)' (N2e)m (N3)r dxdy= 3I!m!n! (1 + m + n + 2)!2A(

29 We find AX. = 4A (aebbj + /ec c,) (3.44) Bj = (1 + j) (3.45) Hi = f e (3.46) If ca,,, a and f are maintained arbitrary within each element, then the matrix elements must be evaluated numerically. Assembling all M elements and imposing the stationary requirement on F yields the complete system M E([Ae]{e} + [Be]{oe} - {He}) = 0 (3.47) e=l Finally, altering the local node numbers to the global ones, (3.47) can be written in a matrix form of the form the same as (3.27), where [K] = [A] + [B] and is also symmetric. 3.3.3. Derivation of Finite Element Equations Using Galerkin's Method In this section, the finite element equation (3.47) is rederived using Galerkin's method and compared with the results derived via the variational method. The residual for (3.30) is R -x ( - y -+7-S-f (3.48) and thus the weighted residual equation for the element e is / s N:Rdxdy = 0 i = 1,2,3 (3.49) Substituting (3.32) into (3.48) and then into (3.49) yields J= 1 JN ' [ (x a 9x a Ny + ay ]ddy = J Ntfdxdy i =1,2,3 (3.50)

30 Invoking the identities Ne -O( ON)\ A0 (1. ON ON )N ax Oax / Ox Ox 'x 9 axxy +a a a oyN y).N aNv i- 1,23 (3.53) a9y ay a S y ay a y ay and the divergence theorem [A e]{ + [B]{} dd= -( + V) {G} (3.54) G e a 9x ax Oy y ' ^ d ( ove r all elements - IINefdxdy~+ ZI l^ (aO+/3O! 9 iedl i=1,2,3 (3.53) where Ce is the boundary of Se and nA is the outward unit vector normal to Ce. In matrix form, the above equations become [Ae]{Je} + [Be]{e }) = ({H} + {G}) (3.54) where the elements in [Ae], [B.] and [He] are the same as those given by (3.40)-(3.42), respectively. The elements in {Ge} are given by Gc= Ni (n l (3.55) The complete system of finite element equations is then obtained by summing (3.54) over all elements M M EZ([Ae]{Oe} + [Be]{qe}) = Z({He} + {Gc}) (3.56) e=l e=l This system is obviously not the same as that in (3.47), since it contains an additional term (the second term on the right-hand side) arising from an integration around the sides enclosing the elements. However, since each of the internal sides,

31 those lying inside the domain, is shared by two adjacent elements, their contribution cancels each other provided that no source of a Dirac delta function form exists there. Therefore, the only contribution comes from the integrals along the external sides lying on C. If the boundary condition is of the Dirichlet type (the values of f are specified on C), the additional term will drop out once the boundary condition is imposed. Alternatively, if the boundary condition is of the Neumann type defined by (3.36), the integrand in the line integral along C becomes zero. Similar conclusion can also be reached when the boundary condition is of the Dirichlet type on portion of C and of the Neumann type on the remaining portion. Therefore, {Ge} drops out when all elements are assembled and the boundary conditions are enforced. The finite element equation derived via Galerkin's method is then identical to that derived via the variational method. 3.3.4. Concluding Remarks In this section, the FEM was applied to a two-dimensional boundary-value problem with a general governing differential equation in conjunction with the Dirichlet and/or Neumann boundary conditions. The finite element equation was derived using the variational and Galerkin's methods with the first order triangular elements. The derivation using high-order elements can, of course, be developed in a parallel manner, but is more involved. The equations obtained here can be used directly to find solution of bounded twodimensional problems. In electromagnetics, such problems include various waveguides filled with inhomogeneous and anisotropic media. One advantage of using the FEM for waveguide problems is the easy treatment of the boundary conditions. Therefore, one can easily deal with a waveguide of arbitrary cross-section. Another advantage is the easy treatment of inhomogeneity and anisotropy of the media. All one needs to do is to specify different physical properties for different elements. The

32 advantages in computational aspects are also notable. The utilization of the properties, such as sparseness, bandedness and symmetry, are key factors in reducing the storage demand and increasing the efficiency of the solution. Due to these advantages, the FEM has gained wide applications in bounded problems. However, if the problem is unbounded (its exterior boundary is at infinity), the FEM alone is no longer applicable, simply because the subdivision of an infinite region will result in an infinite number of elements and nodes. To extend the FEM to unbounded problems, various kinds of techniques, such as infinite elements, boundary elements, etc., have been developed. Though these techniques may appear in quite different form, the principle is the same, that is, to truncate an infinite domain into a finite one where the FEM can be applied. The electromagnetic scattering and radiation in free-space is a typical unbounded problem because of the requirement to enforce the radiation condition at infinity. To apply the FEM to such a problem, several techniques have been developed. In the next chapter, we will describe a hybrid technique which combines the FEM and a boundary element method for a surface integral equation in dealing with twodimensional electromagnetic scattering problems.

CHAPTER IV SCATTERING BY CYLINDERS In this chapter, we consider the problem of electromagnetic scattering by infinite cylinders for transverse magnetic (TM) and transverse electric (TE) wave incidence. The cylinders of particular interest are those having arbitrary cross sections and consisting of inhomogeneous materials. This chapter is organized as follows: first, we give a general review of the problem; then we introduce a finite element-boundary element method with a detailed description of the formulation; following this are the numerical results, which verify the formulation as well as demonstrate the capability of the method; and finally, we compare this hybrid method with other commonly used numerical methods. 4.1. General Review The problem of electromagnetic scattering from infinite cylinders has been studied intensively for many years using exact, asymptotic and numerical techniques. Among various solutions, the exact, which are usually in the form of infinite series, are most desirable, but are only available for very few geometries, such as perfectly conducting and impedance circular and elliptical cylinders, dielectric circular cylinders, and dielectric coated circular cylinders [Bowman et al., 1969; Ruck et al., 1970]. 33

34 Many geometries of practical interest are not solvable analytically. For electrically very small and very large cylinders, the asymptotic techniques, such as the Rayleigh scattering technique for low frequencies, and for high frequencies the geometrical optics (GO), the physical optics (PO), and a more accurate technique, the geometrical theory of diffraction (GTD), can be applied to obtain approximate solutions. However, there are still many problems which are practically important, but difficult to solve using either exact or asymptotic techniques. In such cases, numerical techniques can be used to solve the problems. As discussed in the first chapter, all numerical techniques can be categorized into three fundamental groups: the integral equation (IE) methods, the partial differential equation (PDE) methods, and the hybrid techniques. Among the IE methods, the method of moments (MM) is the best-known [Harrington, 1968] for solving scattering problems. This method is efficient for problems involving perfectly and imperfectly conducting cylinders and homogeneous dielectric cylinders. For such problems the method utilizes surface integral equations (SIE) which reduce the problems into one dimension. The method, however, is not efficient for solving problems involving inhomogeneous media. This is because the Green's functions in inhomogeneous media are usually unknown and thus volume integral equations (VIE) have to be used. As a consequence, the whole cross sections of cylinders need to be discretized, giving rise to a large number of unknowns, which may make the problem practically unsolvable due to limited computer storage and processing time. Recently, various iterative methods have been developed to complement the MM. Among them, the conjugate gradient-fast Fourier transform (CG-FFT) method appears to be a good candidate. The main advantage of the iterative methods is to avoid using a large amount of storage, which is a major drawback of the MM. The PDE methods alone cannot be efficiently applied to finding solutions of

35 scattering problems, since the regions of solution are unbounded or infinite. To make the PDE methods applicable, various hybrid techniques have been developed to truncate infinite regions into finite ones. The best-known of such techniques is the unimoment method, developed by Mei and his students [Mei, 1974; Chang and Mei, 1976; Morgan and Mei, 1979]. In the unimoment method, a PDE method, either the finite difference method (FDM) or the finite element method (FEM), is used to generate a set of trial functions inside a mathematical boundary encompassing the scatterer by solving a sparse or uniformly banded matrix. The field outside the boundary is expanded in the Hankel functions and then coupled with interior solutions on the boundary. Another example of the hybrid techniques is the so called hybrid element method presented by Jeng and Chen [1984], in which the exterior field is again expressed by an expansion of the Hankel functions. Such an expression is imposed on the variational equation for the interior field, which is then solved using the finite element and boundary element methods. However, this method is less efficient than the unimoment method, since it results in a partly full and partly sparse system matrix. In addition to the aforementioned two methods, we also note several other methods, including the one developed by Marin [1982] and the finite element-extended boundary condition method (FEM-EBCM) by Morgan et al. [1984]. In this chapter, we present a finite element-boundary element (FEl-BE) method which is a special form of the hybrid PDE-IE method using the FEM for the partial differential equations and boundary element method (BEM) for the SIE. The basic idea of this method was introduced by McDonald and Wexler [1972] in the early 1970's. The method formulates the exterior field using the SIE and the interior field using the FEM, and then couples them by applying the boundary cond nitions on the artificial boundary separating the exterior and interior regions. This method has been further improved and used to solve two-dimensional antenna and scattering

36 problems [Washisu et al., 1979; Orikasa et al., 1983; Guo, 1985]. In the previous works on scattering, two major shortcomings exist. First, a circle is often used as an artificial boundary to separate the exterior and interior regions. Although such an approach simplifies the formulation, it is not efficient when dealing with non-circular cylindrical problems, especially those with thin geometries. Secondly, the system matrices resulting from the timethod of Jeng and Cen [1984] and the previous FE-BE formulations are partly full and partly sparse, and such matrices are not easily amenable to banded matrix algorithms. These two shortcomings severely limit the capability of these methods. In the method we present here, these two shortcomings no longer exist. The first is overcome by developing a formulation which allows the artificial boundary to be of arbitrary shape. The second is overcome by modifying the coupling technique, resulting in a sparse or uniformly banded system matrix, as in the unimoment method. As a result, the method presented here is more efficient and more capable of dealing with complicated two-dimensional scattering problems. A detailed -formulation is given in the next section. 4.2. Formulation The problem under consideration is illustrated in Figure 4.1, where an electromagnetic wave is incident on a cylinder consisting of dielectric and magnetic materials and perfect conductors. The constitutive parameters of the materials are denoted by a complex permittivity c and a complex permeability A, which may be functions of position. The axis of the cylinder is assumed to be along the z-direction; hence, there is no variation of field quantities and material parameters with respect to the z-coordinate. We will consider TM and TE wave incidence separately. The solution of an arbitrary wave incidence can be obtained from TM and TE solutions, since an arbitrary

37 y r nA "A (a) 5 4 3 2 1 rFA 18 6 4 2 ~- FB 3 / 3 4 33 2.3 1 3 F 70 668 6e6 64 626 69 67 65 63 | 61 65 64 63 62 61: 130 128/ 126/ 124/122/ /1 2 9 1/ 2 7 |/ 2 5 /1 2 3 1 2 1 95 94 93 92 91 C (b) Figure 4.1: Geometry of Wave Scattering by a Cylinder. (a) Illustration of Artificial Boundaries. (b) Example of Discretization.

38 wave can be expressed as the sum of TM and TE waves. Note that the formulation we describe below is for a single cylinder; for multiple cylinders the formulation would be similar. 4.2.1. Interior and Exterior Decoupling As the first step in the FE-BE method, the infinite solution region is divided into two subregions: one is an infinite and homogeneous region and the other is a finite region where the FEM is applied (see Figure 4.1). For this, let the boundary of the cylinder be F. We first draw an artificial boundary FA to enclose F and then draw another artificial boundary FB lying between FA and F (some criteria for drawing FA and FB are discussed in Section 4.3). The region inside FA is treated using the FEM, and thus is called interior region or FEM region and is denoted by RI. The region outside FB is treated using the SIE, and thus is called exterior region or SIE region. The fields in the FEM and SIE regions are coupled by imposing the continuity conditions over the overlapping region between FA and FB. The displacement of FB from FA avoids the singularity problem typically encountered in the Green's function integration. In general, to minimize the FEM region, FA is drawn in such a way that only one layer of elements is needed between FA and F. 4.2.2. Finite Element Analysis In the interior region denoted by RI, the axial field component satisfies the partial differential (wave) equation V. (VEz + ko2r(r)E = 0 (4.1) for the TM case, and 12 V. VI7tI + koir(r)Hz = 0 (4.2)

39 for the TE case, where V = -x8/9x + yd/9y is the two-dimensional operator, ko = 2r/A is the free-space wavenumber, and A is the free-space wavelength. Note that (4.1) and (4.2) are special forms of (2.16) and (2.17) for the source-free case. It can be shown that the solutions to (4.1) and (4.2) can be obtained by minimizing the functional F = / 1 VEz VEz- krr E d - dl (4.3) JIJi P/r(r) Z Or z z AJ d nA for the TM case, and F = I {( VHZ VHzkO r(r)HZ Hz} dsI - Ha zdl (4.4) I Cr(r) J 'WA cl nA for the TE case, where hA denotes the outward unit vector normal to rA. Following the FEM analysis described in Section 3.3, we subdivide RI into M triangular elements. Using the linear interpolation for each element, we can express the field q (= E, for the TM case, and Iz for the TE case) in RI as M 3 (zxy) = E N(, y)c) (4.5) e=l i=1 where Ne(x, y) (i = 1,2,3) is the interpolating function for the eth element and he is the field at the ith node in the eth element, all well-defined in Section 3.3.1. Application of the FEM analysis to (4.3) or (4.4) results in the matrix equation [K]{A } = {b} (4.6) where the element expression for the matrix [K] can be obtained from Section 3.3.2, {( } is a column representing the discretized field, and {i} is a column resulting from the discretization of the line integrals in (4.3) and (4.4). In more detailed form, (4.6) can be written as KAA KAI O~A O A (4.7) L KIA KII Li 0 L J

40 where the subscript A denotes the nodes on FA, I the nodes interior to rA. The above equation is unsolvable since {4A} is unknown. However, such a problem can be bypassed by discarding the first row of (4.7) if we can find another relationship between {4A} and {I,}. The surface integral equation (2.29) can provide such a relationship. 4.2.3. Surface Integral Discretization In the exterior region, the field can be represented by the surface integral equation (2.29) which is rewritten here for convenience as +(r) = IN()- J [Go(B)r 0) - (B) Go(, B) c)( ~ ~~~___r)-cICr-f G('B0~s -4(B0~frB)l dlB (4.8) where rB is on FB, hB is an outward unit vector normal to FB, and G0o(r, B) is the two-dimensional free-space Green's function. Using the BEM, we can discretize (4.8) into a matrix equation so that it can be solved together with (4.7). Assuming that rB passes through MB elements, we can rewrite (4.8) as (er) - N (nC() L [Gt(e qr(rB) a- aGo(r rB] dlB where the quantities in the integral can be expressed as where the quantities in the integral can be expressed as (4.9) Go(r, rB) aG0(r, rB) 9nB rB) 9^(rB OnB = -H (kolr-rBl) koH2)(koIr- rBI) (I -X B cos sin ae) 4 |- + I\r- rB i 3 3 = E t (XB,B)i = E1 Z(ai + bixB + CiyB)Oi i=l =l rNi (O( xBYB) os e + i(XB, YB)sin a i= OXB aYB J 1 3 2,e Z(bi cos ae + ci sin ace)q i=l

41 n. B2 5 B1 Figure 4.2: Illustration of a Line Segment on FB. with sinae = (XB2 - XB)/L, cosae = (YB2 - YB)/L, L = |rB2 - rBlI and (XB1,YB1) and (XB2,YB2) are defined in Figure 4.2. In the above expressions, Nie(x, y) is the same as that in (4.5), which along with as, bi, ci and Ae is defined by (3.33). If we use the simplest one point integration, we can write (4.9) as MB ( 3 2e3r ()) = INC(+) _+ y tHo (kolr- rl)(bicosa + csin, a) (4.10) e=l 8Ati=1 - koH2 )(koIrB-r|) I f cos a + sin a ) (ai + bixB + ciYB)] s } |r- rB| +rB - rl si n where (XB, YB) is taken to be the mid-point of le. Applying (4.10) to calculate the fields at the nodes on rA, we obtain the matrix equation { A} = I{ INC - [PAA]{(A} - [PAI]] { 4I } (4.11)

42 Matrix equation (4.11) provides another relationship between {(JA} and {4I}. We can think it as a boundary condition or boundary constraint on FA. Since (4.8) satisfies the radiation condition at infinity, the solution of (4.11) should also satisfy the radiation condition. Thus, we have converted the radiation condition at infinity to a boundary constraint on the artificial boundary. 4.2.4. Interior and Exterior Coupling In Sections 4.2.2 and 4.2.3, we obtained a finite element equation (4.7) which can be written as [KIA]{$A} + [KI]{1,i} = 0 (4.12) where the first row in (4.7) is discarded here, and a boundary constraint (4.11) which can be written as [PAA]{)A} + [PA1]{I} = A{NC} (4.13) where [PAA] = [I] + [PAA] with [I] as the unit matrix. A joint solution of above two equations gives the numerical result of the fields at nodal points. There are several approaches to the solution of (4.12) and (4.13). The simplest is to solve the (NA + NI)x(NA + NI) matrix equation P'A PA INC PAA PAI A1 A (4.1.4) L KIA Kn j I 0 L where NA is the number of nodal points on rA and NJ is the number of nodal points interior to rA. However, this approach is not efficient. In the following, we discuss two alternative approaches. The first, which is the one commonly used, imposes (4.13) on (4.12) and gives the system of equations [K'I]{,I}l = {i} (4.15)

43 where [KI] = [Kll] - [KIA][PAA] [PA,], {f} = -=[KiA][PAA] {A} In this approach, two matrices have to be solved: one is the complex and full matrix [PAA] having size of NA x NA, and the other is the partly full and partly sparse complex matrix [K'I] having size of NI x NI. Usually, NI is much larger than NA, and hence the size of the scatterer to be treated by this method is mostly limited by the magnitude of NI. The second approach, the one proposed in [Jin and Liepa, 1988b], substitutes (4.12) into (4.13) and gives the equation [PAA]]{qA} = {IANC} (4.16) where [PAA] = PA] - [PAA I][KI] -[KA] Mathematically, this second approach is equivalent to the first one; however, computationally it is much more efficient. Here, we also need to solve two matrices: one is the complex and full matrix [PAA], but now the other is the symmetric and sparse matrix [KII] which becomes real-valued for lossless scatterers and can be narrowly banded if one numbers the nodes properly. A more obvious comparison is given in Table 4.1. The difference between the first approach and the second approach is in the properties of the matrices [KNI] and [AII]. Solving the symmetric, sparse or uniformly banded matrix [KII] is of course much easier and more efficient than solving the nonsymmetric matrix [K'I] with nonuniform block submatrix structures. Three algorithms, namely the banded matrix algorithm, the frontal algorithm and the Lanczos algorithm, which are discussed in Section 3.2.4, can be well applied to solving [K,,].

44 Table 4.1: Comparison of the Two Approaches. Matrix to Matrix Approach be solved order Properties of matrices be solved order [PAA] NA Complex and full 1st approach Complex, partly full and 1 partly sparse [PAA] NA Complex and full 2nd Symmetric, sparse or banded. approach [K] N Real for lossless materials and complex for lossy materials Before ending this section, we wish to point out that the above described interior and exterior coupling is actually a result of the continuity conditions. The continuity of the field and its normal derivative on the artificial boundary rB is satisfied by using the same field expansion (4.5) and solving (4.12) and (4.13) jointly. 4.2.5. An Alternative Formulation As we have seen, the general principle of the FE-BE method is first to use an artificial boundary to split the whole solution region into an interior and an exterior region, then to apply the FEM to formulating the interior field and the SIE to formulating the exterior field, and finally to couple these fields by enforcing the field continuity conditions. Any formulation which complies with the above principle will be able to solve the problem. The formulation described above is just one version of this general formulation. In the following we propose an alternative approach, in which the two artificial boundaries rA and rB will not be used. Let us use R to denote the region occupied by the cylinder and use r to denote

45 the boundary of the cylinder. In this case the functional to be minimized is F J= I VEz- VEz - k2cr()Ez Ez}ds j I () E zdl (4.17) O an for the TM case, and F = / { -)VHzV *H - ko2,lr(r)Hz * H ds- I()Hz dl (4.18) r () r r( and for the TE case, where ni denotes the outward unit vector normal to IF, and Ez and Hz denote, respectively, the electric and magnetic fields on r when the field point approaches F from the inside. Applying the FEM analysis to (4.17) and (4.18), we obtain the matrix equation [BTrp]{r} + [KTT]{qT} = 0 (4.19) where we use {I$T} to denote the discretized nodal fields on and inside r, and use {fir} to denote the discretized quantities (1/Lr)9Ez /9n for the TM case and (1/Er)Hz /lan for the TE case on F. In (4.19), the element expression for [KTT] is the same as [KII] in (4.7), while the element expression for [BTr] can be easily derived. The matrix [BTr] comes from the line integral on F, and hence is a boundary element matrix. In addition to (4.19), another relation between {Nr} and {q$T} is needed to provide a solvable system. The field outside F can be formulated using the surface integral ecuation (-) = 'INC(r)- [GO(r)() - Go( ] d' (4.20) where 0+ denotes the field on F when the field point approaches F from the outside. The discretization of (4.20) using the BEM provides the matrix equation [Prr]{r} + [PrT]{OT} = -{N} (4.21) where {qT} is the same as that in (4.19), and {4r} here denotes the discretized quantity 0o+/9n.

46 The boundary conditions on F require that E- = E+, H-H+ 1 2Ez_ OEZ 1 OH; _ aH 1r An An ' r On On As a result, the discretized field {(} and the quantity {?k} on r in (4.19) should be equal to the corresponding ones in (4.21). Therefore, by solving (4.19) and (4.21) jointly, we can obtain the solution to the problem. The two approaches discussed in Section 4.2.4 for solving the matrix problem are also applicable here. Three differences exist between this alternative formulation and the previous one. First, the previous formulation uses two artificial boundaries while this one does not use any artificial boundary. Second, the previous formulation does not use the line integral in the functional while this one does. Third, the previous formulation avoids the singularity of the Green's function while this one encounters the singularity in deriving (4.21); however, this singularity problem can be easily handled. These three differences are not independent; they are actually intimately related. From computational aspect, the two formulations have about the same efficiency. However, the previous one is easier to implement. The results presented in the next section were obtained using that formulation. 4.3. Numerical Results Based on the formulation described in Sections 4.2.1-4.2.4, a computer program was written for computing the nodal fields. Once the nodal fields are found, the far fields are calculated from the surface integral equation (4.8) using the large-argument approximation for the Hankel function. The program can be used to compute the bistatic and backscattered fields from cylinders having arbitrary cross sections and consisting of inhomogeneous materials for both TM and TE wave incidence. In this

47 section, we will first use this program to investigate two problems pertaining to the subdivision of the FEM region and then present six numerical examples. 4.3.1. General Considerations Before the program is used to obtain final results, we need to resolve the following two questions: first, what is the optimum distance between IA and F, and second, how fast does the numerical result converge with respect to the density of nodes? Let us consider the first question. Since we use linear interpolating functions for the elements between FA and F, the normal derivative of the field on FB, q/OlnB in (4.8), is approximated by the difference of the fields on FA and F. From this consideration, to make this approximation accurately, the distance between rA and I should be as small as possible. However, if this distance is very small, it will result in very narrow elements. Such narrow elements will introduce errors in the FEM analysis. Also, because of the finite (digit) resolution of computers, a very small distance can produce a so small difference between the fields on FA and r that it is beyond the last effective digit, resulting in a meaningless value. Therefore, the distance between FA and F should be properly chosen to best approximate the normal derivative of the field, and thus to provide the most accurate result. To find this optimum distance, we made computations for a coated circular conducting cylinder with a TM plane wave incidence. The radius of the conducting cylinder is 0.4A, and the coating thickness is 0.06A. The coating is a lossy material with r = 2 - j2 and r = 2 - j2. For computation, we divided the coating into two layers and the circumference into 52 segments. Such a subdivision can provide a sufficiently accurate result for the TM incidence. We computed the error in both amplitude and phase of the far-field coefficient (P) in the backward direction for different distances (d). The result is shown in Figure 4.3, where we see that a good

48 (a) Amplitude 4 >0m AI O 3 2 1 0o0.00 0.01 0.02 0.03 0.04 0.05 d/lambda (b) Phase 0.06 3 Q) 'O l0) A.lb, OM am. O) 2 1 0 -- 0.00 0.01 0.02 0.03 0.04 0.05 0.' d/lambda Figure 4.3: Error in Backscattered Far-Field Coefficient vs Distance between rA and F. 06

49 result (error in amplitude within 1% and error in phase within 1~) is obtained for distances of d between 0.005A and 0.03A. Therefore, a criterion for drawing IA is such that the distance between VA and F is within the range 0.005A to 0.03A. This criterion is generally applicable and is independent of the polarization of the incident field and the complexity of scatterers. Considering the second question regarding the convergence of the result, we find it difficult to give a clear answer. Though we know that generally the numerical result converges to the exact solution as the density of nodes increases or the area of each element decreases, we cannot find a general criterion to relate the accuracy of the result to the density of nodes. This is because such a relation depends on the polarization of the incident wave, the property of the material, and the geometry of the scatterer. To show this, let us consider the same problem as that used to generate data for Figure 4.3. The far-field coefficient in the backward direction was computed as a function of the number of nodes (N) along the circumference and the number of layers (L) in the coating. The results are shown in Figure 4.4 for the TM incidence and in Figure 4.5 for the TE incidence. In both figures, we see that as N increases, the far-field coefficient converges to a constant, and as L increases, this constant converges to the exact solution, which is P = 0.3737/14.2~ for the TM case and P = 0.3000L - 142.8~ for the TE case. However, we note that the convergence rates for the TM and TE case are different. For the TM case, the result converges to the exact solution rapidly, while for the TE case, such a convergence is rather slow. This is because in the TE case the coating can support a surface wave which is circulating around the cylinder, resulting in a drastically changing field distribution, while for the TM case the field changes slowly since surface waves do not exist. Therefore, we conclude that the convergence rate of the result depends on the individual case and one needs to check this convergence for each case. For a problem whose exact solution is unknown, perhaps the best way to assess the

50 (a) Amplitude 0.38 -0.37 - 0.36 -0.35 -0.34 5 m- L=1 ~-L-2 OM:u I 1 0 I 30 I 50 20 30 40 50 60 N (b) Phase 18 I 0 0L C: 17 - 16 - 15 - 14 - 13 - 12 - -X- L=1 * L=2 11 I I 10 20 30I I 40 50 20 30 40 50. 60 N Figure 4.4: Convergence of the FE-BE Solution for the Backscattered Field of a Coated Circular Cylinder for TM Case.

51 (a) Amplitude 0.40 0.35 X 0.30 0.25 0 L=2 0.20-.,,. 10 20 30 40 50 60 N (b) Phase -130 -140 CL - -150 L=3 -160 -160 - -- i -- i -- -- i --- * -- i --- * -- i --- -- 10 20 30 40 50 60 N Figure 4.5: Convergence of the FE-BE Solution for the Backscattered Field of a Coated Circular Cylinder for TE Case.

52 convergence of the numerical result is to perform computations for different density of nodes. To give an example, consider a plane wave scattering by a square cylinder with a size of 0.5A x 0.5A, a permittivity of Or = 2.5 and a permeability of,, = 1.5. If we divide each side of the cylinder into N segments, the resulting number of nodes is (N + 1)2. We computed the backward scattered field in the broadside direction for different N's. The results are shown in Figure 4.6 for both TM and TE cases, where we see that the results converge as N increases. The result for the TM case converges faster than that for the TE case. Such a difference in convergence rate is due to the difference in the values of c, and r,. Generally, the result for the TM case converges:fast for a small value of itr and slowly for a large value of #, while the result for the TE case converges fast for a small value of Or and slowly for a large value of c,. Though no general criteria can be obtained, less general criteria can be still developed for certain class of problems. For example, for dielectric cylinders involving no sharp edges, a subdivision of 16 points per material wavelength can provide an accuracy of ~10~ in phase and 5% in amplitude (corresponding to ~0.4 dB in RCS) or better for backward far-field computation. As experience is gained by performing various computations, one can develop a skill for predicting the accuracy of result based on his or her understanding of the field behaviour within the scatterer. 4.3.2. Examples Six examples, selected from numerous computations, are presented here to show the validity, versatility and capability of the FE-BE method (also called as FEMBEM). The first three have unknowns (or nodes) of about one hundred and were computed using the first coupling approach discussed in Section 4.2.4. The last three have unknowns as many as about two thousand and were computed using the

53 (a) Amplitude 1.2 1.0 I3 0.6 -0.4 -0.2 -0.0 - - E-pol. * H-pol. Il -. v * - 6 8 10 12 14 16 18 20 N (b) Phase 0) -0 0) Q. O% C1 6 8 10 12 14 16 18 20 N Figure 4.6: Convergence of the FE-BE Solution for the Backscattered Field of a Square Dielectric Cylinder.

54 second coupling approach. In these examples, the incident wave is assumed to be a plane wave; however, the program can handle non-plane wave incidence as well. A. A Coated Circular Cylinder This is the case for which an exact solution is available and hence is used to verify the program. Figure 4.7 shows the results of our numerical (FEM-BEM) and exact (series) computations for both TM (E-polarization) and TE (H-polarization) incidence. The dimension of the cylinder, the coating thickness and the coating electric properties, given in Figure 4.7, are the same as that used for investigation in Figures 4.3-4.5. For the numerical FEM-BEM computation we divided the coating into two layers and the circumference into 32 segments. Therefore, the dimension of the resulting matrix equation is 64 for E-polarization and 96 for H-polarization. The FEM-BEMN results agree with the exact solutions quite well. Better results can be obtained by using smaller elements and more nodal points as shown in Figures 4.4 and 4.5. Also given in Figure 4.7 are the results computed by using the moment method (MM) with impedance boundary condition (IBC() applied at the cylinder's surface [Knott and Senior, 1974]. In this computation the normal incidence impedance value (ry) which is appropriate for coated conducting plane was used. In this case, the MMIBC also provides good results. B. A Coated Rectangular Cylinder The second example concerns the backscattering from a coated rectangular cylinder. The backscattering cross section, denoted by oa, was computed using the FEMBEM and MM-IBC, and the results are shown in Figure 4.8, where the dimension of the cylinder, the coating thickness and the coating electric properties are also given. From the figure, we see a large difference between the FEM-BEM and MM-IBC results. To assure the validity of the FEM-BEM solution, we computed the same problem using a rigorous moment method with the volume-surface integral formu

55 4.0 3.0 IPI 2.0 1.0 0.0 180 90 arg(P) 0 -90 -180 0 20 40 60 80 100 q, degrees (a) 120 140 160 180 0 20 40 60 80 100 q, degrees (b) 120 140 160 180 Figure 4.7: Bistatic Scattering Pattern of a Coated Circular Cylinder. (a) Amplitude. (b) Phase.

56 8 4 h I4I t 0 0 U 0* I 1* 1 * cr/A (dB) W = 0.92A h = 0.17A t = 0.06A Er = 4 -j4 Pr= 1 77 = 0. 1102 + jO.4368 000 FEM-BEM *.. VSIEM **MM-IBC -8 -16 10 0 -10 -20 -30 -40 3 * I* * I3* 0000000000000000 * I 3** I - I -- - I I I I 0 -10 20 30 40 50 60 70 80 9( q$, degrees (a) cr/A (dB) *0 **0 * 4~~W ~0*** * *""g* * I 00 00 wz00s a - 00 0 0 0 00 0 D ~~I I I I I I I I I 0 0 20 30 40 50 60 70 80 90 q$, degrees (b) Figure 4.8: Backscattering Pattern of a Coated Rectangular Cylinder. (a) TM Case. (b) TE Case.

57 lation (VSIEM) [Jin and Liepa, 1989b]. As seen, the VSIEM results show a very good agreement with the FEM-BEM results. Therefore, the large difference between the FEM-BEM and MM-IBC results is due to the failure of the MM-IBC. Recalling the good accuracy of the MM-IBC results for the coated circular cylinder case, we conclude that the impedance value 71 used for computation, which is appropriate for the coated conducting plane, is an improper approximation near edges. Techniques such as the FE-BE method which is capable of finding accurate near fields can be used to compute this impedance value, from which a better guideline to predict the impedance near edges could be developed. C. An Ogival Cylinder with Coated Edges An ogival cylinder shape is often used to study the scattering behaviour from bare and coated edges. Here we present a sample computation using the FE-BE method for bare, single- and double-edge coated ogival cylinders. In this case the coating is homogeneous lossy material, but of varying thickness as shown by a sketch in Figure 4.9, where the computed backscattering is presented. The ogival cylinder is 1.OA wide and the maximum coating thickness is 0.1A at the edge. From the results we see that the single- and double-edge coating has very little effect on the return at broadside incidence, but at edge-on the double-edge coating decreased the return by some 5 and 16 dB for E- and H-polarizations, respectively. By examining the return at edge-on from the single-edge coated case, it is evident that for E-polarization the leading edge is the dominant scatterer, since coating the leading edge reduces the return (q = 0~), while coating the rear edge the return (q = 180~) is unchanged. In the case of H-polarization, usually the dominant scatterer is the rear edge (due to travelling wave effects), but in this example such is not obvious, since the cylinder is too narrow to sustain and launch appreciable travelling waves. D. An Inhomogeneous Circular Cylinder Starting with this example, we now consider the problems of having thousands of

58 8 0 ac/A (dB) -8 -16 0 20 40 60 80 100 X, degrees (a) 120 140 160 180 10 0 -10 r/A (dB) -- 0 ~;~: 0o ' e~~I@' 'o+* "** ~ * L* o ** n* ~ o** *. 0* 0 0* 0o 01.2 o 0o 0* o* o 0.0 * o - 0* 0 0 0 o -20 -30 -40 I I I I I I I I -W 0 20 40 60 80 100 X, degrees (b) 120 140 160 180 Double(b) TE Figure 4.9: Backscattering Patterns of a Bare, Single- and Edge Coated Ogival Cylinder. (a) TM Case. Case.

59 20. y Er=(2.4X-r)/r \ Ei'Hi L Pr1.0 -X 10. \2.4 X —~,, % \ IA -10.. -pol. 00.4 X, ' 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 4 (degrees) Figure 4.10: Bistatic Scattering Pattern of an Inhomogeneous Circular Cylinder.,,,, ##!,',/i % ## %% r r! o~~~ - ~-, / r'. ~. ss i, ~ 1! -10.-E-po., ":I pol/ Cylinder.

60 unknowns. Here is an inhomogeneous dielectric circular cylinder having radius 1.2A and the permittivity of the dielectric varying in the radial direction according to Er = (2.4A - r)/r. Since the permittivity tends to infinity at the center, a conducting cylinder of radius 0.2A is used to replace the dielectric there. The results for bistatic scattering are shown in Figure 4.10. For computation, the cross section of the dielectric cylinder is subdivided into 3948 triangular elements, resulting in NA = 94 and NI = 2068. E. An Inhomogeneous Rectangular Cylinder This example is an inhomogeneous rectangular cylinder of size 3.8Axl.2A. The permittivity of the cylinder varies along the y-direction as Cr = 1 + cos(y7r/a), with a = 3.8Ai being the height of the cylinder; the permeability is a constant, having Pr 1.5. Both bistatic and backscattering results are given in Figure 4.11. In this case, the number of triangular elements is 3648, resulting in NA = 200, NI = 1925. F. An Inhomogeneous Triangular Cylinder The last example is an inhomogeneous triangular cylinder with the length of each side equal to 2.77A. The permittivity and permeability vary along the x-direction as Er = 1.5 + 0.5x/h and pr = 2.0 -0.5x/h, respectively, where h (= 2.4A) is the height of the triangle. The results are presented in Figure 4.12, again for both bistatic and backscattering. In this computation, the number of triangular elements is 3844, resulting in NA = 186, NI = 2016. All above computations were performed on an Apollo Domain workstation (using Domain series 3000 and 4000 machines). The symmetry property of some examples was not utilized, though it could have been easily adopted to reduce the number of unknowns. Although there is no codes available at current to check the results of the last three examples, some computations were performed using the VSIEM for similar cylinders with reduced size and excellent agreement was observed between the FEM-BEM and VSIEM solutions [Jin et al., 1988c].

61 30. 20. 10. l "C0 IzC --- "e 0. -10. -20. -30. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. ~ (degrees) (a) 0. -10. -20. -30. -40. -50. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 4 (degrees) (b) Figure 4.11: Plane Wave Scattering by an Inhomogeneous Rectangular Cylinder. (a) Bistatic Scattering. (b) Backscattering.

62 20. 10. P7 0. -10. -20. 60. 80. 100. 120. 140. 160. 180. ) (degrees) (a) 10. 0. -10. -20. -300. 20. 40 '. 80. 100. 120. 140. 160. 180. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. ) (degrees) (b) Figure 4.12: Plane Wave Scattering by an Inhomogeneous Triangular Cylinder. (a) Bistatic Scattering. (b) Backscattering.

63 It should be noted that in the FE-BE method there exists a problem of singular frequencies due to the employment of the SIE. This problem occurs when the frequency of the incident field is at the resonant frequency of the interior of the boundary which is FA for the formulation in Sections 4.2.1-4.2.4 and F for the formulation in Section 4.2.5. At such frequencies the system matrix becomes singular and, thus, near and at these frequencies the solutions may be erroneous. Since the resonant frequencies depend on the shape of the boundary, the problem can be avoided by deforming the artificial boundary so that the resonant frequency is shifted away from that of the incident field. 4.4. Comparison with Other Methods In order to give an objective evaluation of our FEM-BEM, we compare it with two other methods: the unimoment method [Chang and Mei, 1976] and the conventional volume integral equation (VIE) method [Harrington, 1968], both of which have been widely used in scattering computations of inhomogeneous cylinders having both e and u different from their free-space values. In Table 4.2 we list the number of nonzero matrix elements that need be generated and the size of the matrices to be solved in each of the three methods. In the FEMBEM, [KII] has about 7NI nonzero elements, since a node is usually connected to six adjacent nodes, and due to the symmetry of the matrix, only 4N1 elements need be generated. 2NA is the number of nonzero elements in [KIA], and 2NA the number of nonzero elements in [PAA] and [PAI]. As previously shown, two matrices need to be solved: one is [PhA] of size NA x NA, and another one is the symmetric and sparse matrix [KII], which in banded form is of size NJ x B with B being the half bandwidth (including the diagonal) of the matrix. In the unimoment method, as originally developed, an artificial circle is drawn

64 Table 4.2: Comparison of the Three Methods. Method. FE-BE method Unimoment method VIE method Elements to generate 4N,+2NA+2NA| 4N+2Nc+4Q2 9N, Matrices to solve N x B, NAX NA NxB', 2Qx2Q 3Nix3Ni to enclose the entire scatterer. Assume N is the total number of the interior nodes and NC the total number of boundary nodes. One needs to generate 4N nonzero elements for the matrix [Q] and 2Nc for the matrix [T] (see [Chang and Mei, 1976], eq. 21). 4Q2 is the number of elements in a matrix of size 2Q x 2Q, resulting from the continuity conditions on the circle, where Q denotes the number of harmonics used to expand the field in the exterior region. In this method, also two matrices need to be solved: one is that with the size of 2Q x 2Q, and the other is the symmetric and sparse matrix [Q], which in banded form has the size N x B', where B' is its half bandwidth. In the VIE method, one needs to generate and solve a matrix of the size 3Ni x 3NI, since at each node there are three unknown equivalent polarized current components. In order to compare the present FEM-BEM with the unimoment method we need to consider three cases. For a scatterer having its cross section close to a circle, the magnitudes of NI and N are about the same, and so are the magnitudes of B and B'. Since the number 2Q can be smaller than NA for far field calculation, the unimoment method in this case is more efficient than the FEM-BEM. However, for a slender scatterer whose cross section substantially deviates from a circle, N is much larger than NI, and B' is also larger than B. Even though the number 2Q is possibly smaller than NA, the FEM-BEM in this case is expected to be competitive with the unimoment method. Finally, for the case of multiple scatterers, and especially when the scatterers are far apart, the FEM-BEM may prove to be more efficient. This is because in the FEM-BEM the artificial boundary rA can be split into several

65 contours, each enclosing one scatterer, while in the unimoment method a single circle is used to enclose all scatterers. However, a modification to the unimoment method, which uses a boundary conforming to the surface of the scatterer, would enhance its efficiency. One such modification, the FEM-EBCM, was presented by Morgan et al. [1984] for the case of dielectric bodies of revolution, and if applied to two-dimensional problems it would result in an efficient numerical technique. Comparing the present FEM-BEM with the conventional VIE method, we find that the FEM-3EM, like the unimoment method and the FEM-EBCM, is far more efficient than the VIE method. For the problems of size comparable to those of last three examples considered in the previous section, using the VIE method, one would have to generate many more matrix elements, about four hundred times the number needed in the FEM-BEM. The computational differences in solving the matrices is expected to be even larger, though no accurate number can be stated, since the computing time depends on the algorithm chosen for solving the matrix. We note that for the case of dielectric (i,. = 1) cylinders, the number of unknowns in the VIE method reduces to 2NI for H-polarization and to NI for E-polarization. We also note that the volume integral equation involving three equivalent polarized current components can be transformed to an integral equation having both volume and surface integrals but only involving one axial field component [Jin et al. 1988c]. Such a volume-surface integral equation provides a technique (VSIEM) much more efficient than the VIE method; however, it still would not be competitive with either the FEM-BEM or the unimoment method. 4.5. Concluding Remarks In this chapter, we described a technique which combines the finite element method and a surface integral equation for solving the scattering by infinitely long

66 cylinders. This technique relaxes the restriction that the artificial boundaries be circles to the condition that the boundaries follow the shape of the scatterer, and thus solves the problem numerically more efficiently. Numerical results were presented to validate the technique as well as to demonstrate its versatility and capability. We note that the corresponding computer program has been successfully used to generate solutions to verify other methods [Jin et al., 1988c; Jin and Liepa, 1989b; Ricoy and Volakis, 1989]. We also note that since 1987 when we presented the technique at the IEEE International AP-S Symposium/URSI Radio Science Meeting, we have seen more researchers attracted into this technique. As a result, more papers were presented on this subject at 1988 and 1989 AP-S International Symposia/URSI Radio Science Meetings [Wu et al., 1988; Collins and Volakis, 1989; Yuan et al., 1989; Boyes and Seidl, 1989]. Though they appeared in different forms and different names, the basic principle behind them is the same as presented in this chapter.

CHAPTER V SCATTERING BY COATED WEDGES AND HALF-PLANES In the previous chapter, we presented a hybrid method to deal with scattering by cylinders whose cross sections are of finite extent. In this chapter, we will consider scatterers whose cross sections are of infinite extent. Two such typical scatterers are the wedge and half-plane which still attract considerable interest in the field of electromagnetics. We will show that the method presented in the previous chapter can be extended to deal with such problems by incorporating the physical optics approximation. The particular problem we consider here is the scattering by conducting wedges and half-planes coated with dielectric and/or magnetic materials. We will first give a brief introduction to the problem and describe a hybrid formulation for E-polarization. Numerical results are then presented and discussed, followed by a short conclusion. 5.1. Introduction Electromagnetic scattering by a perfectly conducting wedge and a half-plane are two canonical problems in electromagnetics. Their analytical solutions are wellknown and can be found in many reference books [Bowman et al., 1969]. By coating 67

68 the wedge and the half-plane with dielectric and/or magnetic materials, their scattering patterns can be changed. This concept is often used for the purpose of radar cross section reduction. For these problems, no exact analytical solutions exist. An approach to the problems is to model the coated wedge (or the half-plane) by a wedge (or a half-plane) of the same geometry, but with an impedance boundary condition (IBC) applied to its surface, and then to obtain the solution using readily available methods [Senior, 1952; Maliuzhnets, 1958; Tiberio, 1985; Volakis, 1986]. However, due to the edge effects the surface impedance may not be uniform near the edge, and, even worse, the impedance value cannot be predicted there accurately. Hence, the applicability of the above approach remains questionable. One way to obtain valid results is to go to numerical techniques. For many years, numerical techniques have been developed for solving various scattering and antenna problems. Among them, the method of moments (MM) is best known [Harrington, 1968]. However, the finite element method (FEM) has advantages for dealing with problems involving inhomogeneous media [Silvester and Hsieh, 1971; McDonald and Wexler, 1972]. Nevertheless, both the MM and FEM are only applicable to electrically small or medium bodies. For dealing with electrically large bodies, such as wedges and half-planes, various hybrid techniques have been developed. Following is a brief review of the previous works. In 1971, Morita [1971] applied the MM to TM scattering from a conducting half-plane by subtracting out the physical optics (PO) current. In 1975, Burnside et al. [1975] combined the geometrical theory of diffraction (GTD) and the MM to treat a conducting wedge by adding a diffraction term to the PO current. Soon after, Wu and Tsai [1977b] applied the same technique to a dielectric wedge; however, in their formulation, only the PO current was used, since the diffraction current is unknown for dielectric wedges. More recently, Newman [1985; 1986] presented a MM/Green's function solution to the

69 scattering by a dielectric and/or magnetic cylinder in the presence of a half-plane. In his formulation, a half-plane Green's function was used as the kernel of integral equations. Note that all these efforts employed the MM. In this chapter, a numerical technique is presented for investigating the scattering by a coated wedge and half-plane. The basic concept employed here is similar to that used by Morita [1971], Burnside et al. [1975] and Wu and Tsai [1977b], except that a finite element-boundary element (FE-BE) method [Jin and Liepa, 1988a], instead of the MM, is used to compute fields near the edges. The use of the FEM allows one to efficiently handle more complicated problems involving inhomogeneous media. In Section 5.2, a procedure of combining the FEM with surface integral equation (SIE) and PO formulations is described for a coated wedge for E-polarization. The analysis for a coated half-plane and for H-polarization is similar. In Section 5.3, numerical results for the surface fields are presented and discussed to show the validity and applicability of this technique. Also presented are the corresponding surface impedances calculated from the surface fields. 5.2. Formulation Consider the scattering by a coated wedge with a TM (E-polarized) plane wave incidence. The coating is characterized by complex permittivity e and permeability p. To solve such an unbounded problem, the whole solution region is divided into three subregions (see Figure 5.1). The first region, where SIE applies, lies outside of abode. The second region is enclosed by abAA1 and edEE1. The coating in this region is assumed to be homogeneous and of uniform thickness, and thus the PO approximation can be used. The interior of the contour ABCDEFA forms the third region, in which the FEM is used to compute the field. In this region, the coating can be inhomogeneous and the geometries can be arbitrary. The problem is finally

70 Incident wave i Ez,Hz FEM region P.O. and S.W. region C A' C' walls P.O. and S.W. region E e Figure 5.1: Geometry of a Coated Perfectly Conducting Wedge.

71 solved by coupling these three subregions. Following is a detailed description of the technique. In the FEM region defined by the contour ABCDEFA, the electric field Ez satisfies the wave equation (4.1) which is rewritten here as V Ir( )VEz + kr(F)Ez = 0 (5.1) where Er(r) and Pr(r) denote the relative permittivity and permeability, respectively. The solution to (5.1) can be obtained by solving the variational problem 6F = O (5.2) with F = / {E VEz - koCr(r)Ez Ez} ds - XEV - E E, ndl (5.3) ABCDE r,( r) where RI denotes the FEM region, ko is the free-space wavenumber, and ni is the outward normal unit vector. Note the line integral over EFA is neglected since we assume the wedge to be a perfect conductor. The solution to (5.2) can be obtained only if one has complete boundary information on ABCDE. If the interfaces between the FEM region and the PO region, AB and DE, are far enough from the edge, the fields there can be approximately represented by the PO fields. As a consequence, one can obtain an approximate boundary condition VE, * n = UE, (5.4) where the unknown coefficient U will be derived later. The line integral in (5.3) can then be written as L E-zVEz * ndl ABCDE r (r) = iIE,.Ezdl + -i Ez Edl + LE/VEz. ndl (5.5) AB,rl JDE tr2 BCD

72 where the subscripts 1 and 2 refer to the corresponding parameters on AB and DE, respectively. The first two integrals in (5.5) can be discretized while the third one remains unknown. If we use subscript B to refer to nodes on BCD and I to refer to other nodes in the FEM region, application of the finite element analysis to (5.2) together with (5.3) and (5.5) results in the matrix equations [ABB]{)B} +[ABI]{4I} = {JB} (5.6) [AIB] {fB}+ [AII]{qI} = 0 (5.7) where {q} denotes the discretized fields Ez at nodes. {JB} in (5.6) comes from the line integral on BCD, which is unknown. However, the problem can still be solved if another relation between {q(B} and {fI} can be found to replace (5.6). Such a sidestepping of (5.6) was justified in the previous chapter. The required alternative relation between {JB} and {qfi} can be found using the SIE together with the PO approximation. In the SIE region, which is the outside of abcde, the field can be computed using a surface integral equation E(- E (o(r() -= (- {G0(, r')VE(i'). i* - Ez(r')VGo(r, r') * '} dl' (5.8) abcde where r' refers to integration point on abode and n' the outward normal unit vector. Go(r, r') is the well-known two-dimensional free-space Green's function. Using the PO fields to approximate the fields on ab and de and defining E()- EINC(r)- { }dl'- i{ }dl (5.9) one can write (5.8) as Ez(r) =E() - { }dl (5.10) (5.10

73 In accordance with the boundary element method, applying (5.10) to compute the fields at nodes on BCD, one then obtains the matrix equation {JB} = {qB$} + [PBB]{qB} + [PBI]{q$I} (5.11) where the subscripts B and I imply the same as in (5.6) and (5.7). Thus the matrix equation (5.11) provides the second relation between {qJB} and {Iri}. By jointly solving (5.7) and (5.11), one obtains the fields on BCD {+B} = [PBB]-1 {~ (5.12) where [PBB] = [I] - [PBB] + [PBI] [A] 1[AIB]. The fields inside can then be computed from (5.7) as {fi} = -[A,]- [AIB]{qB} (5.13) Knowing the fields {4B} and {l,}, one can easily find the scattered field. To implement the above procedure on a computer, one needs to derive the PO field expression to be used in (5.9) and the coefficient U in (5.4). To do this, consider a TM plane wave incident upon a coated conducting plane at an angle 0 from the normal (see Figure 5.2). The total field above the coating is the superposition of the incident and reflected fields EP~ r) = EINC(r) (1 + R,-j2koycos) (5.14) where R is the reflection coefficient given by A+ jB R= J with A = Z cos 0 tan(kt cos 0'), B = Zo cos 0', k = wf/i Z=- -, 0 '-cos - 1 - -- E C 6rIr

74 I < - x -n o t ~, 6 = so 1I Figure 5.2: Geometry of a Coated Perfectly Conducting Plane. The coefficient U in this case is then U = -jko sin (5.15) and this along with (5.4) is valid everywhere (above and within the coating) on the plane perpendicular to the x-axis. 5.3. Results and Discussion To show the validity and effectiveness of the technique described above, a computer code was developed to calculate the surface fields on coated wedges and halfplanes for unit amplitude plane wave incidence. Before the final results can be obtained, the following questions must be first resolved. The first one concerns the integration over the PO regions ab and de in (5.9). The PO regions extend to infinity; however, in computations one needs only to integrate over a finite section, since the Green's function and its derivative provide an integrand

75 whose magnitude monotonically decreases with an increasing argument. The length of the section to be integrated also depends on the direction of the incident wave. For a desired accuracy this length can vary from several wavelengths to hundreds of wavelengths for different directions of incidence. Hence one has to judiciously choose a criterion to truncate the integration. In the half-plane case, when the direction of incidence is near edge-on, the range of integration will be minimum, i.e., the integration converges rapidly. However, when the incident angle approaches 0~ or 360~, the integration range will be maximum. Exploratory computations show that for incident angles within 135~ < ( < 225~, a 10A (free-space wavelengths) integration range is sufficient. This problem is analogous to that of near-field calculations for a travelling wave antenna. The second question is how far from the edge the PO approximation becomes valid. To approach this, one can make an estimate by performing sample computations. The answer depends on many factors, such as the coating thickness, its constitutive parameters, the incident field direction, and the desired accuracy of the results. To give an example, computations were made on a half-plane with one face coated with a dielectric material and the other left bare. The problem is illustrated in Figure 5.3(c) where the coating thickness, the dielectric property and the incident angle q are shown. The FEM region is defined over 0 < x < I and the PO region over x > 1. Computations were made for four different I's and three of them are shown in Figure 5.3. The fourth computation, which is not shown in the figure, for I = 2.0A shows that if the observation point is farther than 1.5A away from the edge the surface fields deviate from the PO fields by less than 2%. It is interesting to observe in Figure 5.3 that even for an inadequately small I the results in the FEM region are still valid. However, this only occurs when 90~ < q < 270". If q < 900 or X > 270~, the approximate boundary condition (5.4) at x = i is not appropriate for small 1, since there are two waves propagating in opposite directions.

76 1.1 1 -C.E: E El 1.0 0.9 0.8 (a) 0.0 0.5 1.0 1.5 8,);,: E 7 - 6 - (b) I 5 i I I I v 0.0 0.5 1.0 1.5 C: 0 (C) 0.0 0.5 1.0 1.5 "0:0 4 3 2 0 -0- 1=0.5wavelength -*- 1=1.0wavelength ' — 1=1.5wavelength - - P.O. field.T^^^-gg^A ~ (d) 0.0 0.5 1.0 x (wavelength) 1.5 Figure 5.3: Surface Fields of a Half-Plane with the Upper Face Coated. (a) Electric Field at Upper Surface S1. (b) Normal Derivative of Electric Field at Upper Surface S1. (c) Electric Field at S2. (d) Normal Derivative of Electric Field at S2. t = 0.1A, er = 2.4-jl1.0, r = 1.0, > = 135~.

77 Since there are no results, to our knowledge, available for a coated wedge or halfplane to validate our results, surface currents of a perfectly conducting half-plane (a special case of a coated half-plane with Er = /r = 1) were computed for various angles of incidence. Very good agreement was observed between the numerical results and the exact solutions. To show the feasibility of the technique, four examples are given below. The first is a coated right-angled wedge, and is illustrated in Figure 5.4(a). The surface electric field and its normal derivative are given in Figures 5.4(a) and (b), respectively. The fields on the two sides of the wedge are identical because of the symmetric incidence (q$ = 135~). From the results we see that the fields change greatly near the edge, while far from edge they tend to constants which, as one would expect, are the PO fields. To check the validity of IBC on such a structure, the surface impedance (Etan/Htan) normalized to free-space intrinsic impedance Zo was calculated and is given in Figure 5.4(c). One finds that at points farther than 0.1A away from the edge, the IBC provides a very good approximation. However, close to the edge, the IBC would be invalid since the computed impedance value deviates drastically from that predicted by PO. It should be noted that in the computations, the phase reference point was usually shifted to the FEM/PO interface, and as a result the phase data for the field and its normal derivative are shifted by a constant phase. However, this has no effect on the surface impedance values. The next example is a half-plane with its edge and two faces coated by a lossless material, as depicted in Figure 5.5(c). The results are presented in Figure 5.5, where the edge effect on the surface fields and on the surface impedance is clearly seen. It is of interest to note that even for a lossless coating the impedance around the edge has a resistive component, indicating energy flowing into the coating. At the lower (shadow) surface the impedance has a negative resistance, which indicates that the energy is flowing out from the coating.

78 1 100 -0- magnitude -.- phase a) 0, E 1.1 1.0 0.0) a) V) Cn -100 m aL (a) o.c O.E 7 6 3 3 -200 0.0 0.2 0.4 0.6 0.8 1.0 0 E. 5 4 2 4 3 3 0.0 0.2 0.4 0.6 0.8 1.0 - 1- (b) C) (0 C 0 a0 (C) n. 0 - a) *'0 23 cE 2 1 0 0.0 0.2 0.4 0.6 d (wavelength) 0.8 1.0 Figure 5.4: Surface Fields of a Coated Right-Angled Wedge. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.06A, Er = 2.0 - j2.0, /tr = 2.0, ~ = 135~.

79 -200 "- magnitude -- phase a 3 0) rC E 100 0 (D ") ra) Q. (a) - -100 - -200 1.0 -r200 0.0 1 i I I -1.0 -0.5 0.0 0.5 8 1 6 0 E 4 2 o t 0-i-1.0 3 - -0.5 0.0 0.5 1.C 100 0) (, "O a 0) -100 -200 140 120 ~I) 100 -C a) 8, a) 80 (b) 2 0) - 0) E (C) 1 0o-1.0 +60 1.0 -0.5 0.0 0.5 d (wavelength) Figure 5.5: Surface Fields of a Half-Plane with its Edge and Two Faces Coated. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.06A, r = 2.0, lr = 2.0, - = 1350.

80 The third example shows an inhomogeneously coated half-plane with edge-on incidence, as sketched in Figure 5.6(c). The relative permittivity of the coating varies along the surface from er = 4 at the edge to r = 2 at infinity according to Or= 2.0(1 + e-x /A ). The results for the upper surface are given in Figure 5.6 and, due to symmetry, are the same for the lower surface. One may note that the examples presented above were made for thin coatings which are not able to support surface waves. As a final example, a half-plane with its upper face coated by a thicker lossy dielectric was considered. For top incidence the results are shown in Figure e 5.7, and there one clearly sees the oscillatory behavior of the surface fields due to the interaction of the incident and reflected fields with the surface wave excited by the edge. We note that even though the present analysis does not include the surface waves in the PO region formulation, the results show such waves and their correct interference patterns in the FEM region. This indicates that the present technique is also applicable to the thick and lossy coating cases. It should be noted that for lossless coatings which are sufficiently thick (t > 0.25A/(ErIr - ])1/2) to support surface waves, the present technique may not give accurate results. The difficulty will be more pronounced when one extends the technique to the transverse electric (TE) case, since there surface waves are always present, no matter how thin the coating is. To overcome such a difficulty, one would need to incorporate expressions for surface waves in the PO formulation. 5.4. Conclusion This chapter presents a numerical technique for computing TM scattering by perfectly conducting wedges and half-planes coated with dielectric and/or magnetic materials. The technique combines finite element method with surface integral equation and physical optics approximation. Numerical results are presented and discussed

81 1.0 0.8' -- magnitude -*- phase 0) E 0.6 0 a).0) U) QC (a) 0.4.).; O.C 4 2 3 0.0 0.5 1.0 1.5 2.0 - 3 3 0) E.E 1-: 0) 'a. (b) Q 2 1 0.0 2.0, 0.5 1.0 1.5 2.0 80 1 0) 3 rE a) (D () 03 CL (C) 0.0 0.5 1.0 1.5 2.0 x (wavelength) Figure 5.6: Surface Fields of a Half-Plane with its Two Faces Coated Inhomogeneously for Edge-On Incidence. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.1A, e, = 2.0(1 +e-x2/A2), r = 1.0, = 180~.

82 a) ~.. 0) E 1.5 1.4 1.3 1.2 1.1 -0- magnitude -- phase Q) (a) (0 0, Mr M= 0.0 0.5 1.0 1.5 2.0 -0 C 0) E 8 7 6 5 4 174 172 170 168 a (b) (I) 166 X 164 162 0.0 0.5 1.0 1.5 2.0 3 C 0) E (C) 0 Q) 1.0 x (wavelength) Figure 5.7: Surface Fields of a Half-Plane with Thick Coating on the Upper Face. (a) Electric Field. (b) Normal Derivative of Electric Field. (c) Surface Impedance. t = 0.15A, e, = 4.0 - jl.O, lr = 1.0, > = 90~.

83 for coated half-planes and a coated right-angled wedge. The technique allows one to study the scattering behavior of coated wedges and half-planes as well as the surface impedance behavior around the edges of such structures. From computed data one can extract diffraction coefficients, which then may be used to extend the GTD to coated structures. The shortcomings of this technique are also pointed out. Incorporation of surface waves in the formulation may allow to treat thick coatings and TE incidence case.

CHAPTER VI SCATTERING BY INHOMOGENEOUSLY FILLED THICK APERTURES The hybrid technique that combines the finite element method and the surface integral equation has been well-established in Chapter 4 for scattering by cylinders and in Chapter 5 it has been extended to solving for scattering by wedges and halfplanes by incorporating the physical optics approximation. In this chapter, we will deal with the electromagnetic characterization of the transmission and scattering properties of an aperture in a thick conducting plane filled with an inhomogeneous composite material, a problem which is of considerable interest in electromagnetic engineering. The principle of the technique for this analysis is basically the same as that in Chapter 4; however, to allow the treatment of large apertures, the conjugate gradient method and fast Fourier transform will be employed for the solution of the resulting system. Numerical examples will be presented which demonstrate the validity, versatility and capability of the technique. 6.1. Introduction The problem of electromagnetic diffraction by a slot in a perfectly conducting plane of finite thickness has been studied extensively. Lehman [1970] and Kashyap 84

85 and Hamid [1971] employed asymptotic techniques; Neerhoff and Mur [1973] considered a formulation of coupled integral equations; and Hongo and Ishii [1978] utilized the Weber-Schafheitlin integral technique. Also, Auckland and Harrington [1978; 1980] considered modal and nonmodal formulations where the method of moments was employed to solve the coupled integral equations. These solution approaches are usually restricted to rectangular slots or slots of arbitrary cross section filled with a homogeneous or at most partially homogeneous material. Also, except for the asymptotic solutions presented by Lehman [1970] and Kashyap and Hamid [1971], all others, especially the one by Auckland and Harrington [1980], require excessive storage to invert the resulting matrix. In this chapter we present a method that alleviates both of those limitations. The method combines the finite element method (FEM) and the boundary integral equation formulation. The FEM [Zienkiewicz, 1977] is known to be well-suited for dealing with material inhomogeneity and geometry irregularity and results in a very sparse matrix. It has recently been used by Jeng [1988] to compute the admittance matrix of an air-filled cavity-backed aperture for transverse electric (TE) incidence. In this chapter the FEM is employed to formulate the fields within the slot and establish a relationship with those at the aperture. The fields external to the slot are expressed as an integral over the aperture and a system of integral equations is then obtained by enforcing field continuity across the aperture. These are solved via the conjugate gradient method (CGM) and in the process the boundary integrals are efficiently evaluated via the fast Fourier transform (FFT). The capability of the FEM to treat structures of arbitrary geometry as well as material composition and that of the CG-FFT method to deal with large systems are demonstrated herein. Although this combined methodology is applicable to a wide class of configurations, in this chapter we restrict its application to that of a twodimensional slot in a thick perfectly conducting plane with TE or TMV (transverse

86 magnetic) incidence. To show the versatility of the combined FEM and CG-FFT method, computations are presented for slots filled with multilayered dielectrics that may also support embedded conducting strips. 6.2. Theoretical Formulation Consider the geometry illustrated in Figure 6.1. This specific configuration consists of an inhomogeneous material slab inserted between two metallic half-planes of the same thickness. The relative permittivity and permeability of the material will be denoted by c,.(r) and Ip(r(), respectively. For further reference, we will denote the upper half-space (y > 0) as region I, and the lower half-space (y < -t) as region II. Also, the cross section of the slot (O > y > -t) will be referred to as region III. In the following, we first present a general formulation for the scattering by such a structure and then give detailed formulations for TM and TE cases, respectively. 6.2.1. General Formulation For the problem of scattering by an aperture in a thick, infinite ground plane, it is customary to decouple the fields in the three regions by closing the upper and lower apertures of the slot with a perfect conductor and introducing the equivalent magnetic currents M1 and M2 over the extent of the apertures, as illustrated in Figure 6.2. Based on the equivalence principle [Harrington, 1961], M = E x, M2 = E2 x (-) (6.1) where E1 is the field at the upper aperture (y = 0) and likewise E2 is the field at the lower aperture (y = -t). Both currents reside on the surface of a ground plane which can be removed by application of image theory.

87 I y Region I t I --- — Region I-IEgg Region III x 2 Region II Figure 6.1: General Geometry of an Inhomogeneously Filled Thick Aperture.

88 i Region I conducting plane (a) conducting plane ------— ^ — ------ M2 Region II y=O y=-t (b) (c) Figure 6.2: Equivalent Problems for That Shown in Figure 6.1. (a) For Region I. (b) For Region II. (c) For Region III.

89 The field in region I can now be expressed as the radiation caused by M1 and the impressed sources (7, M'). We will denote the tangential (to the aperture) magnetic field at y = 0 as H(M1,7IMi). Similarly, the tangential magnetic field at y = -t in region II will be denoted by H7 (M2) as generated by the equivalent magnetic current M2 over the lower aperture. Enforcing continuity of the tangential electric fields at the upper and lower apertures of the slot, the fields within region III can be expressed as the radiation of the equivalent magnetic currents -M1 and -M2, as illustrated in Figure 6.2(c). The fields in each region are coupled by requiring continuity of the tangential magnetic fields across the slot opening. We then have H,(Ml, J ) = t (-M1, -M ) (6.2) Ht (M2) = Ht (-M1,-M2) (6.3) where Hll and HIII denote the tangential magnetic fields in region III at y = 0 and y = -t, respectively. Traditionally, (6.2) and (6.3) are employed in constructing integral equations for the solution of the magnetic currents M1 and 2. In this procedure, the fields in each region are expressed as an integral over the aperture currents utilizing the appropriate Green's function applicable to each region. In regions I and II the freespace Green's function is required, whereas in region III the appropriate Green's function must be derived to satisfy the pertinent boundary conditions. Clearly, for a slot of arbitrary cross section and filled with inhomogeneous material the derivation of the Green's function for region III may not be possible and this is the major limitation of traditional approaches. In the proposed formulation the variational equation F = 0 (6.4)

90 with F = J/ [ (V x E") (V x EI, ) - k2rEI E I] dxdy -2jkoYo M1 tl dz -H2jkoYo d (6.5) i J2 or F = Jj [-(V x H"'). (V x H 1) - koPr H"I H'] dxdy -2jkoYoj Mi~Htl dx - 2jkoYo M2 Ht2 d (6.6) 1 2 is employed in region III, where ko = 27x/A is the free space wavenumber and Zo = 1/Yo is the free space intrinsic impedance. Also, Q denotes the cross sectional area of region III, F1 is the line segment specifying the upper aperture and F2 is the corresponding line segment specifying the lower aperture. Equation (6.4) can be discretized in the standard manner employed in finite element techniques to generate a matrix relation between the aperture and internal fields. The incorporation of the boundary conditions to be satisfied by the internal fields is rather straightforward after a discretization of (6.4). Essentially, the introduction of (6.4) eliminates a need to find the Green's function associated with the internal structure of the slot. More importantly, the application of (6.4) is independent of the slot's cross sectional geometry and material filling once the appropriate finite element mesh has been generated. The boundary conditions (6.2) and (6.3) in conjunction with (6.4) and (6.5) or (6.6) imply a system of equations for the solution of the magnetic currents M1 and M2 to be employed for the computation of the scattered and transmitted fields. 6.2.2. Formulation for TM Case For TM (E-polarization) incidence, the impressed or incident electric field EinC

91 is z-directed and the equivalent magnetic currents M1 and M2 may be written as M1= xM(x) =-xEI- = — E on (6.7) M2 = M2() = E = ixE on r2 (6.8) The electric field in region I is now given by E(x y)= Enc(x, y) + efl(,y) -2 Ml( (2) (k z-')2+y2 dx' (6.9) 20y J O ) where EZrefl is the electric field reflected from the infinite ground plane in the absence of the slot and HQ2) is the zeroth order Hankel function of the second kind. The corresponding tangential magnetic field at the aperture is given by H(x, O) = 2H (x,0 ) 2 1 +- k 2 J Ml(x')Ho(kolx - x'dx' (6.10) where HnC denotes the x-component of the incident magnetic field. Similarly, the electric field in region II is given by E '(xy) 2 = y - M(x )H2) (ko (+y dx (6.11) and the corresponding tangential magnetic field at the aperture is found as H2l(x 0) = - o ( + f M2(x)H 2)(ko - x')dx' (6.12) 2 ko Ox2 a 2 The functional F to be employed in (6.4) becomes J { r [(Oa )2 (O aI)2]-k (E i)2} dxdy -2jkoZo I MlHIdx - 2jkoZo 2 M2HIdx (6.13) By further incorporating (6.2), (6.3) and (6.7), (6.8) along with (6.10) and (6.12), F can be more explicitly written as F = ilIOEZ)2 +,OE 2 ] k2yE2 P Q fr a x a 9y 0

92 +j J Ez(x) k[( + a2) j E(') (kolx- xI)dx' dx +j 2 Ez(x) [(ko + a) E(x' )H /(kolx - x'l)dx'] dx +4jkoZo j Ez(x)Hrlc(xr,O)dx (6.14) where the superscript III on Ez has been omitted for convenience. A numerical solution for Ez can then be obtained by solving the variational equation (6.4). Once the aperture electric fields are obtained, the scattered field can be computed from (6.9) and (6.11) and the scattering radar cross section (RCS) can be evaluated from [Bowman et al., 1969] a= k- P(>)2 (6.15) where p is the observation angle measured from the x-axis and P(p) is the far field coefficient P(o) given by P() = 2 j Ez(X)jk~ocs sin pdx 0 < p < ir (6.16) P(p) = -ko f Ez(x)ejk~~ sin dx ir < p < 2ir (6.17) The time average power per unit length (in the z-direction) transmitted into region II through the slot is given by Ptrans = -Re fr H(x)E (x)dx (6.18) 2 where the asterisk denotes the complex conjugate. From this, the transmission coefficient T of the slot, defined as the ratio of Ptrans to the time average incident power P.c intercepted by Ir per unit length, can be calculated in a straightforward manner.

93 6.2.3. Formulation for TE Case For TE (H-polarization) incidence, the impressed or incident magnetic field tfnC has only a z-component implying that the scattered magnetic field will also be zdirected. The equivalent magnetic currents M1 and M2 may then be written as M = zM(x) = E= zE on 1 (6.19) z = =-zEI on F2 (6.20) The magnetic field in region I due to M1 is given by ~HI() = zHi(r) = zHzc(x, y) + zHnc(x, -y) - A M(x)H(kolr - xx'j) dx' (6.21) 2 1 and similarly, the magnetic field in region II due to M2 can be expressed as = ~H(r) =) = 2 ~kY j M2(x')Ho2 (ko - I x'l)dx' (6.22) The functional F to be employed in (6.4) becomes F = J {i [( 6x + (OIy ) - k2[ (HZH) dxdy -2jko MlH,,dx -2jkoYo j M2HaIIdx (6.23) In this case, it is impossible to find an analytical expression for the functional F which couples the fields in the three regions by enforcing the boundary conditions (6.2)-(6.3), as is in the TM case. The coupling of the fields or the enforcement of the boundary conditions can only be carried out after we discretize (6.21)-(6.23). This is the main difference between the formulations for the TM and TE cases. Once the equivalent magnetic currents AllM and M2 are computed, the RCS of the structure is again given by (6.15) with P(p) given by P() = -k j Al (x)eik ~ dx 0 < < 7r (6.24) P() = k < < (6.25) P(GO) koYo -2 M~(X)ejkx~o~dx 7 < p < 27r (6.25) 2 JV>

94 and the time average power per unit length transmitted into region II through the aperture is given by Ptrans = -Re H(x)M2(x)dx (6.26) 6.3. Numerical Discretization In this section the functional expressions and integral equations will be discretized using the finite element and boundary element methods and systems of linear equations will be derived. For discretization, the cross sectional region Q is subdivided into M small triangular or rectangular elements. Also, the line segments F1 and F2 are broken into L1 and L2 short segments, respectively. On the assumption of a linear field distribution within each element, the field within the eth element having n nodes can be expanded as n Ez (or Hz) = yNi(x,y)O- = {N }T{qe} = {qe}T{Ne} (6.27) i=l where Ne(x,y) are the shape functions and at represent the nodal fields. In the following we will discuss the TM and TE cases separately since their formulations are different. 6.3.1. Discretization for TM Case For the TM case, we need to derive a system of equations by discretizing and minimizing (6.14). For the discretization of the line integrals in (6.14), we use the mid-point field value to represent the field on the segment. Thus, the field on the sth segment becomes Ez = (S1 + s2) = os (6.28)

95 where ~q1 and q82 denote the electric fields at the ends of the segment. Substituting (6.27) and (6.28) into (6.14), we obtain M L1 L1 F = Z({~ }T O[I(]{,} + Z}1 + E( s2)s E qPt e=l s=1 t=l 2L2 L1 + E(l + ) + 'Pt + 2jkoZo (' + +s2)Hnc(x,, 0)8 (6.29) s=1 t=1 s=1 where 6, denotes the length of the sth segment and xs is its mid-point. The n x n matrix [Ke] is given by [Kf] e a[ da + dy dy 5 2] e } T - kO{e } {N}N' } dxdy (6.30) in which Qe denotes the area of the eth element with (e,ese) being the relatively permittivity and permeability of that element. Finally, PSt = 2 [( ~ 2 tJ-t/2 42 (kox - x')dxl (6.31) 2 L\ d J=,c -i,/2 For constant er and e,, the elements of [Ke] can be evaluated analytically, resulting in the simple expressions given in Section 3.3.2 for linear triangular elements. The integral in Pt can also be evaluated analytically to find Pst = k2 [1- log(0.1638kos,)] 6 -jkoH(2)(koS/2) s=t (6.32) 2 k Ll7k PSt = 2ko H()(kolx - xt) + ko [~H2(kolxs - xt-,/21) T H(2)(kolx, - Xt + t/21)] x,>(xt ~ 8t/2) s t (6.33) where H42) is the first order Hankel function of the second kind. In accordance with the Rayleigh-Ritz procedure, differentiating F in (6.29) with respect to each nodal field and equating the resulting expression to zero yields the system of linear equations [A]{} = {b} (6. 34)

96 where [A] is an N x N square matrix with N being the total number of nodes within Q and on Fr and r2, {q} is a column vector representing the electric field at the N nodes and {b} is a known column vector representing the excitation. It is not difficult to show that [A] is a partly sparse and partly full symmetric matrix. 6.3.2. Discretization for TE Case For the TE case, we need to discretize and minimize (6.23) with respect to HIII to derive a set of linear equations which, together with another two sets of equations derived by discretizing (6.21) and (6.22), forms a system of equations for the solution of the problem. Substituting (6.27) into (6.23) and employing a pulse basis expansion for the magnetic currents M1 and M2, we have M L1 L2 F -= y{e}T[IKe]{qe} -j (3 sl + s -2)s. j (sl + sq2)s Ls (6-35) e=l s=1 s=1 where {(te} represents the nodal magnetic fields within the eth element, qo8l and q>s2 denote the magnetic fields at the ends of the sth segment, and,1 or i/2 denotes the quantity (koYoM1) or (koYoM2) on the sth segment. The matrix [IKe] is given by (6.30) provided e' and,4 are interchanged. Minimizing (6.35) with respect to the nodal fields yields the matrix equation [KI]{} ++ [B1]{41} + [B2]{V2} = 0 (6.36) where {(} is a column vector representing the magnetic field at the N nodes within region III and on r1 and r2. Also, {-1} = koYo{M1} and {f2} = koYo{M2 (6.37) are column vectors having L1 and L2 elements, respectively, equal to the number of segments employed for the discretization of rF and r2, respectively. The quantities

97 in the brackets represent sparse matrices assembled as M Li L2 [K] = E[Ke] [B1] = E[B], 1[B2]= Z[Bs] (6.38) e=1 s=1 s=1 where the elements of the column vector [B8] are given by Note that the submatrices [Ke] and [Bs] are properly augmented when assembled using the global nodal numbers. We also note that [K] is an N x N square matrix whereas [Bi] and [B2] are rectangular matrices having N x L1 and N x L2 elements, respectively. We now proceed with the discretization of (6.21) and (6.22) to derive another two sets of linear equations. Multiplying both sides of (6.21) and (6.22) by ( —j) and applying Galerkin's technique to the resultant equations, we obtain two matrix equations [C1]{l}= {inc} - [P1] {41} (6.39) [C2]{ 2} = -[P2]{Mk2} (6.40) where {l,} and {12} are column vectors having L1 + 1 and L2 + 1 elements, respectively, representing the nodal magnetic fields on rl and r2. The sparse matrices [Cl] and [C2] are assembled from L L2 [C1] = Z[CS], [C2] = Y[Cs] (6.41) s=l s=1 where [CS] is a row vector whose elements are given by C1 - C2 - -9s Again, [C8] is properly augmented when assembled using global numbering. We note that [C1] and [C2] are rectangular matrices having L1 x (L1 + 1) and L2 x (L2 + 1) elements, respectively. The column vector {bfinc} has L1 elements given by C'- = -2jI.nC(x,, 0)6( (6.42)

98 where x, denotes the mid-point coordinate of the sth segment on Fl. Also in (6.39) and (6.40), [P1] and [P2] are square matrices having L1 x L1 and L2 x L2 elements, respectively, given by Ps = - [1- 2 log (0.1638k0s)] s =t (6.43) Pst = -ZH(2)(kol.-,)6,t S t (6.44) The final system of equations is obtained by enforcing (6.2) and (6.3) on (6.36), (6.39) and (6.40). In doing so, we have K B1 B2 0 BT Pi 0 i c = < n (6.45) BT 0 P2 ~2 0 for a numerical solution of the magnetic fields and currents. Note that the superscript T in (6.45) denotes the transpose of the associated submatrix and, thus, the operator in (6.45) is symmetric. 6.4. CG-FFT Implementation While the systems of equations (6.34) and (6.45) can be solved via a direct matrix inversion, this approach requires a memory of O(N2) for (6.34) and 0((N + L1 + L2)2) for (6.45). If the sparsity property of the matrices is efficiently utilized (not necessarily a simple task), task), the memory demand can be reduced to 0(N) + O(L2) + O(L2) which still may exceed the capacity of the available storage. An approach in dealing with this difficulty is to employ an iterative solution process such as the CGM and exploit the convolutional nature of the boundary integrals to reduce the memory demand to O(N). In this section, we describe such an approach. Rewriting (6.34) and (6.45) as Az = b (6.46)

99 and assuming an initial solution vector z1, the CG algorithm [Hestenes and Steifel, 1952] for the solution of (6.46) can be outlined as follows. Initialize the residual and search vector: ri = b-Azl (6.47) o0 = (6.48) (Aar,Aari) (6. Pi =- oAarl (6.49) Iterate for k = 1,2,..., N for (6.34) or N + Li + L2 for (6.45): ak = (, k) (6.50) {ApA, Apk) Zk+l = Zk+ akPk (6.51) rk+l = rk- akApk (6.52) A k (6.53) (Aark+l, Aark+l) (6.53) Pk+l = Pk +k+lAark+l (6.54) Terminate at k = N for (6.34) or N + L1 + L2 for (6.45) or when k+1 12 < tolerance (6.55) As seen, in the implementation of the CG algorithm the coefficient matrix A is only involved in the computation of the two column vectors Ax and Aax, where x is a known column vector and the superscript a implies the adjoint operator. This computation can be accomplished rather efficiently without generating matrix A explicitly and to illustrate this process, let us consider the computation of Az for the TM and TE cases separately. TM Case The contribution of the M area elements to Az is Af E[KIe]{)} (6.56) e=l

100 and since the expression for [Ke] is quite simple, the above summation can be carried out very fast. The contribution of the L1 and L2 line segments to Az can also be computed efficiently via the FFT. For example, if we assume that each segment has the same length 8 and define fs as +[( -):62 (ko Ix dx s -1,2,..., Li (6.57) then, {f} = FFT-1 {FFT{g} o FFT{G}} (6.58) in which the symbol o denotes the Hadamard product. The elements of the column vector {g} are defined by t t =l 1,2,..,Li 9t = (6.59) 0 t-=L +l,L + 2,...,Np and those of {G} are given by k2 [1 - (2j/7r)log(0.1638koS)] 6 - 2koH(2) (ko/2) t = 1 koH2 ((t - 1)ko5) + ko ((t - 3/2)ko) (0) Gt = (6.60) - Hn2)((t - 1/2)ko)] t = 2,3,..., Np/2 + 1 GNp-t+2 t = N/2+ 2, N/2 + 3,..., Np where Np > 2L - 1. TE Case The first N elements of Az are attributed to (6.36) and can be computed as M L1 L2 [I[e]{~Oe} + E[B]{,,} + E[B`]{O0} (6.61) e=1 s=1 s=1 The next L1 elements of Az are attributed to (6.39) and can be computed as L1 >[Cs]{q>O} + FFT-1 {FFT{,1} o FFT{G}} (6.62) s=1

101 where the elements in {G} are given by -(j/2) [1- (2j/Tr)log(0.1638ko6)] 62 t = 1 Gt = -(j/2)H2)(ko(t - 1))6)2 t = 2,3,..., Np/2 + 1 (6.63) GN,-t+2 t = Np/2 + 2, Np/2 + 3,..., Np and 1t = 0 for t = Li + 1, L1 + 2,..., N. Similarly, the last L2 elements in Az are attributed to (6.40) and can be computed as L2 E[C']{d+} + FFT-1 FFTFT2} o FFT{G}} (6.64) 1=1 The computation of Aaz can be accomplished in a similar manner. But, since A is symmetric, then Aaz=A*z, where A*z can be computed using the above formulas with [Ke], [B8], [Cs] and {G} replaced by their complex conjugates. 6.5. Numerical Results Based on the formulation described above, a general computer code was developed and in this section we present some numerical results which validate the formulation and demonstrate the code's versatility. In all cases a plane wave illumination is assumed. 6.5.1. Results for TM Case For the sake of verification, Figure 6.3 shows the equivalent magnetic currents (normal incidence) and backscattering RCS for a rectangular groove, depicted in Figure 6.3(a), and Figure 6.4 shows the equivalent magnetic currents (normal incidence) and transmission coefficient for a rectangular slot, depicted in Figure 6.4(c), both filled with different materials. The results presented are those computed by this method and the method of moments in conjunction with the modal Green's function

102 3.0 2.0 0 1.0 0.0 I - -0.50 -0.25 0.00 0.25 Xx/ (a) 0.50 180. 90. 0. -90. -180. 00 0 I - - -I. -. ---. I. - A I -0.50 -0.25 0.00 x/X 0.25 0.50 (b)

103 20. 10. 0. -10. -20. -30. -40. 90. 120. 150. 180. 4p (degrees) (c) Figure 6.3: Equivalent Magnetic Current at Normal Incidence and Backscattering RCS for a 1A Wide and 0.25A Deep Groove. (a) Magnitude. (b) Phase. (c) RCS. Solid Lines Correspond to r, = 1 and Ptr = 1; Dashed Lines Correspond to Cr = 4-jl and,r = 1. Solid and Dashed Lines: This Formulation; Circles and Squares: Mode Matching Solution.

104 2.0 1.5 z 0 1.0 0.5 0.0 '-0.80 -0.40 0.00 0.40 x/Xl (a) 0.80 2.0 1.5 z 0 1.0 0.5 0.0 '-0.80 -0.40 0.00 0.40 x/X% (b) 0.80

105 9 -~, * V% rA 1.2 1.0 0.8 0.6 0.4 0.2 0.0 90. 100. 110. 120. 130. 140. 150. 160. 170. 180. qp (degrees) (c) Figure 6.4: Equivalent Magnetic Current at Normal Incidence and Transmission Coefficient for a 1.6A Wide and 0.25) Thick Slot. (a) Er = 1, iir = 1; Solid Line: Current on the Upper Aperture; Dashed Line: Current on the Lower Aperture. (b) Er = 2.56(1 - j0.1), ~r = 1; Solid Line: Current on the Upper Aperture; Dashed Line: Current on the Lower Aperture. (c) Transmission Coefficient Corresponding to Er = 1, tr- = 1 (Solid Line); Er = 2.56(1 - j0.1), ltr = 1 (Shorter Dashed Line); and Er = 2.56(1 - jO.1), i, = 1.2 - j0.35 (Longer Dashed Line). Solid and Dashed Lines: This Formulation; Circles: Mode Matching Solution [Wang, 1984].

106 1cm 2cm /////// - "'" - [ ///. 31cm c 4 — *1 Figure 6.5: Illustration of an Aperture Filled with Three Dielectric Layers and a Strip Grating. Top and Bottom Layers: Cr = 2.56, Pir = 1.0, Thickness = 0.2 cm. Middle Layer: r = 4.0, O, = 1.0, Thickness = 0.2 cm. of the cavity [Wang, 1984]. It is clearly seen that the both methods give identical results. To demonstrate the versatility and capability of the combined FE,-BE and CGFFT formulation as well as the pertinent computer code, we consider the computation of the plane wave scattering by the geometry shown in Figure 6.5. The illustrated configuration depicts a wide slot in a thick conducting plane filled with three dielectric layers with a perfectly conducting strip grating embedded in the center layer. The specific geometrical and material properties associated with the dielectric and pertinent strip grating are given in the figure. Figure 6.6 shows the electric fields at the upper aperture, the center line containing the thstrip grating and the lower aperture computed at 10 GHz and for normal incidence. We note that the field distribution within each period of the strip grating is very similar and as a result, the approximate scattering solution for the truncated periodic array proposed in [Jin

107 and Volakis, 1990] will be of reasonable accuracy for large apertures. Some bistatic scattering patterns are given in Figure 6.7 for three different incidence angles. The higher order Bragg diffraction effects are clearly seen in Figures 6.7(b) and (c). Finally, Figure 6.8 shows the transmission coefficient as a function of the incidence angle at 10 GHz and as a function of frequency at normal incidence. Note that in all Figures 6.6-6.8, we have included information relating to the discretization and convergence of the iterative solution. 6.5.2. Results for TE Case For the sake of verification, Figure 6.9 shows the equivalent magnetic currents M1 and M2 as computed by this method and the method of moments in conjuction with the slot's Green's function [Aucland and Harrington, 1978]. It is clearly seen that both results are practically identical. Next, we consider the computation of the plane wave scattering by the geornetries shown in Figure 6.10. The illustrated configurations depict wide slots in a thick conducting plane containing a (truncated) perfectly conducting strip grating embedded in a multilayer dielectric that fills the entire slot. Specifically, Figure 6.10(a) shows a single strip grating centered in the slot and residing within a 0.1 cm gap that separates two dielectric layers. The slot shown in Figure 6.10(b) contains two strip gratings, each embedded in the center of a dielectric layer. The specific geometrical and material properties associated with the dielectric and pertinent strip grating are given in the figure. Figures 6.11 and 6.12 show the equivalent magnetic currents computed at 30 GHz for Figure 6.10(a) and 9 GHz for Figure 6.10(b). We again observe that the current distribution within each period of the strip grating is very similar. Some bistatic scattering patterns for the geometry in Figure 6.10(a) are given in Figure 6.13 at three different frequencies, all for normal incidence. Also,

108 2.0 1.5 N w 1.0 5 0.5 -15.5 0.0 x (cm) (a) 15.5 3.0 2.0 N wz 1.0 ARA ffWAI 0.0 L-15.5 2.0 1.5 0.0 x (cm) (b) 15.5 N 1.0 0.5 I -15.5 0.0 15.5 x (cm) (c) Figure 6.6: Field Distribution Corresponding to the Geometry in Figure 6.5 at 10 GHz and Normal Incidence. (a) Upper Aperture Field. (b) Field at the Center Line Containing the Strip Grating. (c) Lower Aperture Field. L1 = L2 = 186, N = 1309, Tolerance=0.01, 139 Iterations.

109 30. 10. -10. -30. 0. 90. 180. 270. 360. p (degrees) (a) 30. 10. e -10. -30. 0. 90. 180. 270. 360. wp (degrees) (b) 30. -10. '3 -10. -30. 0. 90. 180. 270. 360. qp (degrees) (c) Figure 6.7: Bistatic Scattering Patterns for the Geometry in Figure 6.5 at 10 GHz. (a) 'i nc 90~ (Normal Incidence), 139 Iterations. (b) c'in = 60~, 176 Iterations. (c) W'nc = 30~, 161 Iterations. Li = L2 = 186, N = 1309, Tolerance=0.01.

110 1.2 1.0 0.8 0.6 0.4 0.2 0.0 90. 100. 110. 120. 130. 140. 150. 160. 170. 180. (p (degrees) (a) 1.2 1.0 0.8 H 0.6 0.4 0.2 0.0 0. 5. 10. 15. 20. 25. 30. f (GHz) (b) Figure 6.8: Plots of the Transmission Coefficient (T) Corresponding to the Geometry in Figure 6.5. (a) T versus Angle of Incidence (o) at 10 GHz; L1 = L2 = 186, N = 1309, Tolerance=0.01. (b) T versus Frequency (f) at Normal Incidence; 0-10 GHz: L1 = L2 = 124, N = 625, Tolerance=0.01; 10-20 GHz: L1 = L2 = 186, N = 1309, Tolerance=0.01; 20-30 GHz: L1 = L2 = 248, N = 1747, Tolerance=0.01.

111 2.0 1.5 1.0 0.5 0.0 -0.50 -0.25 0.00 0.25 0.50 (a) 180. 90. -90. -180. L'-0.50 -0.25 0.00 0.25 0.50 (b) Figure 6.9: Magnitude and Phase of M1 (Circles) and M2 (Squares) for a 1A Wide and 0.25X Thick Slot Filled with a Dielectric Having, = -1j. L1 -L2 20, N = 126. Solid and Dashed Lines: This Formulation; Circles and Squares: Data from [Aucland and Harrington, 1978].

112 Figure 6.14 includes three bistatic scattering patterns corresponding to the geometry in Figure 6.10(b) at 9 GHz for different incidence angles. Finally, Figure 6.15 shows the transmission coefficients for the two geometries as a function of frequency at normal incidence. 6.5. Concluding Remarks A hybrid technique was developed for a numerical characterization of the scattering and transmission properties of an inhomogeneously filled aperture in a thick conducting plane with TM and TE incidence. The technique combines the finite element method and boundary integral formulation. This resulted in a system of equations whose solution was obtained via the conjugate gradient (CG) method employing the fast Fourier transform (FFT) for the evaluation of the convolution integrals over the aperture. The use of the CG-FFT eliminated a need for an explicit generation of the pertinent matrices and thus reduced the memory demand to O(N) for the implementation of the entire system. Differing from other techniques, the proposed methodology can efficiently handle any slit cross section regardless of its material properties and geometry. New geometries can be simply accommodated by changing the finite element mesh without affecting the solution process. It is, therefore, a promising approach for treating large and complex radiation and scattering problems that are often formidable with traditional techniques.

113 l1cmj 2cm I" " i"' //////w////'/......... *.@ - b //////l/ 21cm I,I ~r=2.56, d=0.2cm - ~r=4, d=0.2cm - strip, d=0.lcm (a) 1cm 2cm I-. —....aI.. - - ' r-..........- -. *... b.............../ - w — ------ ------------- AM f //////~// 31cm..I PI PI 1-I I r=2, d=0.4cm -.~r=2.56, d=0.4cm - strip, d=0.1cm (b) Figure 6.10: Illustration of an Aperture Filled with (a) a Single Layered Strip Grating and (b) a Double Layered Strip Grating. d Denotes the Thickness of the Associated Dielectric Layer or Conducting Strip.

114 2.5 2.0 1.5 -10.5 10.5 Incidence. = L = 252, N= 1771, Tolerance=0.005, 389 Iterations. ~ ' 5, is 0.0 -10.5 10.5

115 2.0 a a 1.01 ["A I I I i, I i i \ t h i G i F g r I I a1 b I "' ', I I,,,!: -15.5 15.5 Incidence. L = L = 248, N = 1743,,Tolerance=0.01, 159 Iterations.

116 't3 20. 10. 0. -10. -20. 90. 120. 150. 180. 210. p (degrees) (a) 240. 270. 'R 30. 20. 10. 0. -10. -20. -30. 90. 120. 150. 180. 210. (p (degrees) (b) 240. 270.

117 40. 30. 20. 10. 0. -10. -20. Figure 6.13: 90. 120. 150. 180. 210. 240. 270. p (degrees) (c) Bistatic Scattering Patterns for the Geometry in Figure 6.10(a) at Normal Incidence. (a) 3 GHz, L1 = L2 = 126, N = 508, Tolerance=0.005, 97 Iterations. (b) 10 GHz, L1 = L2 = 189, N = 1140, Tolerance=0.005, 118 Iterations. (c) 30 GHz, L1 = L2 = 252, N = 1771, Tolerance=0.005, 389 Iterations.

118 m A/ 30. 20. 10. 0. -10. -20. 0. 60. 120. 180. 240. 300.:360. pq (degrees) (a) V1-N %0 I, 30. 20. 10. 0. -10. -20. 0. 60. 120. 180. 240. 300. qp (degrees) (b):360.

119 30. 20. 10. 0. -10. -20. 0. 60. 120. 180. 240. 300. 360. (p (degrees) (c) Bistatic Scattering Patterns for the Geometry in Figure 6.10(b) at 9 GIIz. L1 = L2 = 248, N = 1743, Tolerance=0.01. (a) Vinc = 30o, 102 Iterations. (b) Sin' = 60~, 201 Iterations. (c) pinc = 90o, 159 Iterations. Figure 6.14:

120 0.8 0.6 F 0.4 0.2 0.0 0. 5. 10. 15. 20. 25. 30. f (GHz) (a) -- 1.0 0.8 0.6 0.4 0.2 0.0 0. 3. 6. 9. 12. 15. f (GHz) (b) Figure 6.15: Plots of the Transmission Coefficients versus Frequency at Normal Incidence. The Circles Represent Actual Computation Points. (a) For the Geometry in Figure 6.10(a). (b) For the Geometry in Figure 6.10(b).

CHAPTER VII APPLICATION OF ISOPARAMETRIC ELEMENTS There are basically two approaches to achieve more efficient and accurate computation of electromagnetic scattering. One approach is to develop new methods. The development of several hybrid methods discussed in the previous three chapters are some examples in this direction. An alternative approach is to improve the existing methods, and here we will consider this approach. Specifically, we will describe the application of isoparametric elements in the finite element-boundary element method presented in Chapter 4 to improve the efficiency and accuracy of the solution for the scattering by cylinders. 7.1. Finite Element-Boundary Element Formulation To illustrate the use of isoparametric elements in finite element-boundary element methods, we choose the formulation presented in Section 4.2.5 as an example. For convenience, we introduce the following notations = E(r), () = ) a(, U(r), v(r) = E(r) for TM incidence and 1(r) = Hz1(r), )=, U(r), (r) = rr) 121

122 for TE incidence. Then, (4.17) and (4.18) can be uniformly written as F = R {u(i)V()V) (- + kOv(r)q(r) (f)} ds - jr) fir)dl (7.1) and (4.20) as INC f= [IN0G-(, I (r) = - (r) - [Go(r, )(r-) - () n ] dl' (7.2) Note that in (7.1) and (7.2) the boundary conditions on F have already been enforced by neglecting the superscripts in (4.17), (4.18), and (4.20). In accordance with the finite element-boundary element method, we need to transform (7.1) into the form of matrix equation (4.19) and (7.2) into that of (4.21) to form a linear system. For this, we first subdivide R into M elements and assume that L is the number of segments on F resulting from the subdivision. Then, we expand 6 in each element as m ()= EVI ()i {v )} } = {T =,I{v} e= 1,2,...,M (7.3) i=l and expand q and b on each segment as +(q( r) = EU}(r) =us = {ps}T{uJe} - =1,2,...,L (7.4) i=1 s( ) = EU(r){ = {us)T{W {'}T{ue} ) = 1,2,...,L (7.5) i=1 where VTe and U7' are the known expansion functions and q$, k? and I' are the corresponding unknown coefficients. Substituting (7.3)-(7.5) into (7.1) yields M L F = E{e}T[Ie]{ f} - E{()T[Bs]{ 7) (7.6) e=l s=1 in which the elements in [JIK] are given by j= J {u()VV() * VV() - v(r)V r)Vij l,2,...,m (7.7)

123 and the elements in [Be] are given by Bi = fjr U( (r)U;()dl i,j = 1,2,., 1 (7.8) where Qe denotes the area of the eth element and Ps the length of the sth segment. Application of the Rayleigh-Ritz procedure to (7.6) will then produce a matrix equation in the form of (4.19). Now substituting (7.4) and (7.5) into (7.2), we have I =1 Jr = ) INC()- [GO(r, I){U}T{s} -{U{)T{q'}V'Go(r, r'). '] dl' for r on r (7.9) where f denotes the Cauchy principle value integral. Application of the boundary element method to (7.9) then provides a matrix equation in the form of (4.21). From the formulation above, it is clear that the accuracy of the solution depends on how accurate (7.3)-(7.5) can represent the fields and how accurate the integrals in (7.7)-(7.9) can be evaluated. The second problem is closely related to the modelling of geometry. 7.2. Isoparametric Elements As we have seen, to obtain an accurate numerical solution one is required to model both the geometry and the fields accurately. Isoparametric elements were developed to provide such an accurate modelling. These elements belong to a special group in the family of parametric elements and particularly refer to the discretization elements which have the same order shape functions as the field representation within them. When using isoparametric elements, the region R can be broken into a number of elements whose sides can be curved. Figure 7.1 shows two examples of breaking a circular region into five and 12 quadrilateral elements, respectively, with

124 (a) (b) Figure 7.1: Two Examples for Breaking a Circular Region. (a) A Five Quadrilateral Element Model with 20 Nodes. (b) A 12 Quadrilateral Element Model with 45 Nodes.

125 I y 2 3 1 c Mo 1 3 2 - I I0 BO.? -1 o +1 x w 0 (a) (b) Figure 7.2: A Quadratically Curved Segment in the xy-Plane Can Be Transformed into a Straight Segment Lying on the - Axis. quadratic shape and expansion functions. In this section we will describe the formulation for quadratically curved line elements and that for quadrilateral elements with quadratically curved sides. 7.2.1. Curved Line Elements It can be shown that a quadratically curved segment in the xy-plane can be transformed into a straight segment lying on the (-axis (Figure 7.2) using the transformation 3 3 x = L()x, y =E L?(~)y i=l i=l with the quadratic shape functions L? (i = 1,2, 3) given by =- 1 22 2 (7.10) (7.11)

126 It is observed that LP() is unity at the ith node and thus by choosing (7.11) for the expansion functions in (7.4) and (7.5), i.e., U? = Li, then Oq will coincide with the nodal values of the field (r) and so will x,,. It can also be shown that the field is a quadratic expansion of x and y. From (7.10), we find d - + IJJ'l d (7.12) and the normal unit vector n is given by = ( - y ) IJ1 (7.13) With these relations the line integrals in (7.8) and (7.9) can be evaluated numerically. Specifically, B = - U (r)U(r)dl = L t(r)L( )|J'(s)7dg (7.14) and the line integral in (7.9) can be written as Is= B(r,r')U2(r')dl'= B(r,r')LS( ')1Js(')1d' (7.15) If a four point Gaussian integration formula is employed, then 4 B,~ = E WkL,. (~)Ls (G)IJs(G) (7.16) k=1 4 Is =- WkB(r, r () (7.17) k=1 where j = -4 - -0.8611363116, ~2 = - =3 - -0.3399810436, W1 = W4 = 0.3478548451, W2 = W3 = 0.6521451549, and r- = x(k)x + Y(~k)Y in which X((k) and Y(Gk) are given by (7.10).

127 7.2.2. Quadrilateral Elements with Curved Sides Similarly, it can be shown that a quadrilateral element with quadratically curved sides in the xy-plane can be transformed into a square in the br-plane (Figure 7.3) using the transformation 8 8 x = ENr(,r)xi, y = ENe(~,?)yi (7.18) i=1 i=1 with the shape functions Ne (i = 1,2,...,8) given by Ne =~(1 - )(1 - r)( + 77 + 1), Ne = 1(1 + )(1 - )( - - 1) 4 4 1 1 N = 4(1 + ~)(1 + r)(~ + 77 1), Ne = (1 1 + (1)(- + V/- 1) N = (1 - 2)(1 -), N6 -(1 + )(1 - 2) 1 1 N = (1 - 2)(1 + 7), N = (1 - 0( - 2) (7.19) Again it can be shown that Nje(~, a) is unity at the ith node and thus by choosing them for the expansion functions in (7.3), i.e., Vie = Nit, then qt will coincide with the nodal values of the field Oe(r). Also, the field expansion can be shown to be a quadratic function of x and y. From (7.18) the area element ds can be expressed as ds =(x ayd a d~d I JeldFd? (7.20) where IJel is the determinant of the Jacobian transformation matrix. The area integrals in (7.7) can then be evaluated numerically. Specifically, - J-1 U( ) [ 9xx 9' 9y Jy ] - k2v()Ne, r7)N(, )} IJe(E, )Idd (7.21) where ___e _1 V ye "y (1rj r(7.22) 9xf J 9 xy -IJeI a ax a ay aNj I 1 ( adONSj +x Na,j (7.23) 17 o 0OO~ O? /

128 A 1n (0,1) 7 (-1,1) 4 2 1 3 ( 8 (1,1) 3 (1,0) 2 (1,-1 ) 4,- v 0 1 5 (-1,-i) (0,-i) (a) (b) Figure 7.3: A Quadrilateral Element with Quadratically Curved Sides in the xy-Plane Can Be Transformed into a Square in the br/-Plane.

129 By using a 3x3 point Gaussian integration formula, (7.21) can be evaluated according to the formula 1 1 3 3 1 _ A(, 7)dd1 = EE WWjA(i, j) (7.24) I I ~i=1 j=l where 41 = - = -0.7745966692, 72 = /2 = 0.0, '3 = q3 = 0.7745966692, W1 = W3 = 0.5555555556, and W2 = 0.8888888889. We note that the numerical implementation of the above formulation is rather straightforward. The quadrilateral elements are just one type of two-dimensional isoparametric elements. Another type of elements is trilateral element which is also often used and its formulation can be found in the book by Zienkiewicz [1977]. 7.3. Conclusion In this chapter we have described two types of isoparametric elements: a quadratically curved line element and a quadrilateral element with quadratically curved sides. The use of these isoparametric elements in the finite element-boundary element methods can lead to an improved efficiency and accuracy not shared with the traditional elements. Since such an improvement is well-quantified in mechanical and structural engineering, it is expected that similar improvement can be obtained in electromagnetic scattering problems.

CHAPTER VIII FORMULATIONS FOR THREE-DIMENSIONAL SCATTERING For long cylindrical geometries with axial excitation the electromagnetic wave scattering can be well-approximated by two-dimensional models where one deals with only scalar fields. However, all practical geometries are finite. For the situation where there is a symmetry, such as a body of revolution with axial or circumferential field excitation, the problem can be still reduced into a two-dimensional one. But, for the geometries having no such a symmetry, a description by three independent coordinates and consequently a full three-dimensional treatment is required. The numerical analysis of three-dimensional scattering has been investigated using integral equation approaches with the method of moments (see, for example, Schaubert et al. [1984]). The rapid increase of the number of unknowns in the three-dimensional case, however, prohibits such a method for scatterers having a size comparable to the wavelength of the excitation field. Recently, the high sparsity of the matrices produced by the finite element method has attracted more people, and an increasing amount of effort has been concentrated on developing various finite element solutions. The vector property of the electromagnetic fields, however, makes such a development very difficult. In this chapter we will propose several techniques where the finite element method will be employed in conjunction with either a surface integral equation or an eigenfunction expansion. 130

131 8.1. Problem Description Consider a three-dimensional scattering problem illustrated in Figure 8.1, where the volume V is occupied by an electrically and magnetically permeable body with constitutive parameters e and,p which are complex if the material is lossy. The surrounding medium is assumed to be free-space with constitutive parameters co and po. An electromagnetic wave excited by a source with harmonically oscillating angular frequency w is impinging upon the body and then scattered. Let us denote the infinite volume exterior to the volume V as Vo and the closed surface separating V and V, as S. The electromagnetic fields (El, H1) inside the volume V satisfy the following vector partial differential (wave) equations 1 -- V x -V x E1- w2E1 = O (8.1) /2 1 V x -V x 7H1- w2JH1 = 0 (8.2) which actually are the special forms of (2.14) and (2.15) for the source-free case. The fields (E2, H2) outside the volume V satisfy the free-space wave equation V x V x E2 - kE2 = 0 (8.3) V x V x 2-k- kH2 = O (8.4) where ko(= w/Eo/t7) is the wavenumber of free-space. The general principle of the various methods for solving such a problem is as follows: a three-dimensional finite element method is used to formulate the fields inside the volume V while a surface integral equation or an eigenfunction expansion, which is a solution for the fields in VO, is used to provide a necessary boundary constraint on the surface S for the finite element solution. We will address the finite element formulation first in the next section.

132 V0 n - (E2,H2) Incident Wave / \ I I. \ (0o t0o ) Figure 8.1: Geometry of Wave Scattering by a Permeable Body.

133 8.2. Finite Element Formulation In this section we describe the formulation of the finite element method for the fields inside the volume V. 8.2.1. Stationary Functionals of Fields To solve (8.1) and (8.2) using the finite element method, we first find the corresponding functionals in the volume V. In the entire domain V + VOO the functional for (8.1) is [Silvester and Ferrari, 1983] F =, (V E) ).(V x El)-eEl - E1 dV +2 / vo [1 —(V x E1) ~ (Vx E1)-E1 *. El dV (8.5) Substituting the vector identity 1 1 2(V x E1) (V x E) = E Vx( Vx E1) -V (V x El-) x E (8.6) into the second integral of (8.5) and applying the divergence theorem, we obtain F = 1 (V x E1)-(V x El)-cEl- El] dV F = 2 Isw2 [(V x El) x El * n dS (8.7) where n is the normal unit vector pointing from V to Voo. Equation (8.7) is the functional expressed in terms of the electric field. A similar expression can be derived in terms of the magnetic field, which is F = l// [- (vx Hl) (v x Hl)-iH 1 HI] dV +2 /J s' - [(V x H1) x 7] * n dS (8.8) 9 t,2,

134 Note that the surface integrals in (8.7) and (8.8) are the same and they can be written as 2w (Ei X Hi). n dS (8.9) The scalar forms for (8.7) and (8.8) are respectively given by F = 1 i 1j [(OEz aEy2 (Ex OEz)2 2 W2 IL\ 9 ZJ \a z axJ +( a aEd )] E+E]}E dV + (EHm - EmHi) dS (8.10) and Fv= 1 J { (% z 0.)2 (a.H O9z2 F = 1 y 2 J Jv [ y z ) z x ) + y - - + + Hz2] d +2 Jj(EIHm- EmHi)dS (8.11) where the subscript "1" on both E and H has been dropped for convenience. Also, I and m are two orthogonal unit vectors tangential to the surface S and they are so oriented that (1,, mn) form a right-handed system. 8.2.2. Finite Element Discretization Let the volume V be subdivided into a number of three-dimensional elements which can be either tetrahedra, rectangular prisms, or triangular prisms, or even better, isoparametric elements. Within each element, the fields can be expressed as n n TE N= E(xy,z He = N(xy, (8.12) i=l i=l

135 tions, and (ET,H ) represent the fields at the ith node. Substituting (8.12) into of linear equations For example, application of the three-dimensional finite element discretization to (8.10) results in the matrix equations [K{] {E} + [IKs] * {E} = 0 (8.13) and [Ksi] {E} [Kis] {Es}t + [s] {sTi} = (8.14) where {Ei} denotes the electric field at the nodes interior to the surface S, {ES} denotes the tangential electric field at the nodes on S, and {77s} is the tangential magnetic field at the nodes on S. Also, [K] denotes a matrix having three dimensions. It is obvious that the finite element equations, such as (8.13) and (8.14), do not form a complete system and another matrix equation relating {Es} to {Hs} is required. In the next section we will discuss several formulations to provide such a relationship. 8.3. Hybrid Formulations In this section we present several formulations which combine the finite element method with either a surface integral equation or an eigenfunction expansion for solving the three-dimensional scattering problems.

136 8.3.1. Finite Element-Surface Integral Formulation In this formulation a surface integral equation, either (2.26) or (2.27), is employed to provide a matrix equation relating {Es} to {Hs}. For example, (2.26) can be written as E(R) = EZ (R) -( { V x Go(R, R [n' x E(R )] -jo o(R, R') * [n x H(R')] } ds' (8.-5) The discretization of the above equation on the surface S provides a matrix equation in the form of [6s].{ }+ [ss] { s} = {NC (8.16) where -{EINC is the tangential incident electric field at the nodes on 5, and again, B] denotes a matrix of three dimensions. Equations (8.13), (8.14) and (8.16) then form a complete system for the solution of the nodal fields. Clearly, this formulation is the same as that presented in Section 4.'2.5 for the twodimensional case. Another formulation corresponding to the one in Section 4.2 can also be developed for the three-dimensional case by introducing two artificial surfaces to enclose the scatterer and, as a result, the singularity problem in discretizing the surface integral equation can be avoided. 8.3.2. Finite Element-Eigenfunction Expansion Formulation It can be shown [Stratton, 1941] that the scattered field outside V can be expanded in the spherical vector harmonic functions M4 and N4 as E(WR) = E D [fMT/(koR) + gN(koR)] (8.17),,

137 where fV and g, are the unknown expansion coefficients, D, is a normalization constant, and v is a combined index incorporating the three spherical harmonic indices. The scattered magnetic field is then given by 7s(R) = V x E((R) = jo E, [fiNv(koR) + gJM(koR)] (8.18) and the total fields in the exterior region by E(R) = ENC(R) + E(R) (8.19) H(R) = HI (R) + H8(R) (8.20) Next, using (8.19) and (8.20) to evaluate the tangential electric and magnetic fields at each node (point-matching) on the surface S, we obtain the two matrix equations {Ts} = {E [PE] { } (8.21) 9V {Hs} {7HINC} + [P] } (8.22) gv where {ff,} and {g,} are column vectors representing the unknown expansion coefficients properly truncated from the infinite series. From (8.21) and (8.22) a matrix equation relating {Es} to {Hs} can be found as {Es} - [P1] [72]-. {H} = {-EC} - [P1] [72]-. {H } (8.23) which together with (8.13) and (8.14) also forms a complete system for the solution of the nodal fields. It should be noted that in the above formulation the number of harmonics used in (8.17) must be equal to the number of nodes on the surface S. However, the necessary number of harmonics to expand the exterior fields can be much smaller. Obviously, the above formulation does not exploit this advantage. In the following we propose an alternative formulation.

138 Let us first expand the tangential fields on the surface S as {Es} = c, {Ts,}; {Hs} = Ec {Hs,} (8.24) 11 A where {Es}4 and {Hs,} are the coupled basis functions, and cA is the expansion coefficient. The number of expansion terms in (8.24) is equal to the number of the harmonics in (8.17). {Es,} and {Hs,} can be generated by selecting a set of orthogonal known bases for {Es} in (8.13) and (8.14) and then solving (8.13) and (8.14) for {Hs}. Now let the surface S be divided into some smaller surfaces Si whose number also equals the number of the harmonics. Substituting (8.24) into (8.19) and (8.20) and integrating over Si (subdomain method) yields I Z ccLEsdS = J j EINcdS + j E9dS (8.25) c cHsdS = HINCdS + HdS i = 1,2,3,... (8.26) These two equations, after employing (8.17) and (8.18), can be transformed into the two matrix equations [1] {c} = {einc} + P1] v } (8.27) K2] {c} = {hinc} + [H2] fv (8.28) from which one can eliminate {c,) and then solve for {f,) and {g,}. 8.3.3. Finite Element-Extended Boundary Condition Formulation This formulation was proposed by Morgan et al. [1984] and applied to the case of axisymmetric inhomogeneous penetrable scatterers. In the formulation the surface integral equation embodied in the extended boundary condition method (EBCM)

139 [Waterman, 1965 and 1979] is combined with the finite element method, as shown below. According to Stratton [1941], the incident field can be expanded in the spherical vector harmonic functions M1 and N as EINC (I) = D [a Ml(koR) + byNl(koR)] (8.29) where a, and b, are the known expansion coefficients. The scattered field is again expanded as (8.17). In accordance with the EBCM, one can find the following four basic integral equations ko2 s (N Es - jYoM4 Hs) dS = j7ra, (8.30) ko2 f/ (M4 Es - jYoN. Hs) dS = j7rb, (8.31) and J7rf = kJo j(N Es - jYoM Hs)dS (8.32) j~rg = ko f (AMl Es-YoN l Hs)dS (8.33) Now assume that 1ES + E d S'1 (8.34) { -cd Hs A HA S." Sit where {Ts5} and {Hs}4 are again the coupled pairs of basis functions, and c, and d, are the expansion coefficients. {Es,} and {Hs} can be generated in the same manner as discussed in the previous subsection for (8.24). The expansion in (8.34) is actually the same as that in (8.24) except that it is now split into two parts for description convenience. Substituting (8.34) into (8.30)-(8.33), one obtains the two matrix equations 1 i4 j4( 83 a I^ ^ l = |(8.35) jr I' 4K L4 d,, J 1/11 I /I11 b I

140 and { fg 9v I1 1 c 1 IIJA JVJA A.7r I K' Ll d IIJA via A (8.36) where Jv k f(N.Es -jYoMt Hs). ko Es( - jYoM. Hs) = kJ j(. E * -JYoN H4) with t = 4 in (8.35) and t = 1 in (8.36). Eliminating {ca} and (8.36), one finds 1-1( {J}i 1p J [1 I4 J4 v| Iv L1 Kt.K L 1 bv and then the scattered field can be computed from (8.17). dS dS dS dS (8.37) (8.38) (8.39) (8.40) {d,} from (8.35) and I (8.41) 8.3.4. Discussion Among the above three three-dimensional hybrid formulations, the one in Section 8.3.1 has been verified in the two-dimensional case in this thesis and the one in Section 8.3.3 has been verified in the case of a body of revolution [Morgan et al., 1984]. Only the one proposed in Section 8.3.2 has not yet been verified. The major shortcoming in the finite element-surface integral formulation is the need to generate two full matrices [Bss] and [Bs], whose dimensions are proportional to the number of nodes on the surface S. These two matrices can be so large that their memory demand can easily exceed the capacity of a computing facility. A remedy for this is to use a rectangular box as the surface S, then compute the surface integrals via the fast Fourier transform and solve the system of equations using

141 the conjugate gradient method, as demonstrated by Collins and Volakis [1989] for the two-dimensional case. One obvious advantage of this formulation is that there is no need to compute the spherical vector harmonic functions, as in the other two formulations. The finite element-eigenfunction expansion formulation with point-matching appears very simple; however, it is not efficient since the matrices [P1] and [P2] are also full and have the same dimensions as [Bss] and Bss]. The alternative formulation with subdomain method is more efficient though it involves one more matrix inversion. Instead of using delta functions (point-matching) or subdomain functions (subdomain method) as the testing functions as in the previous formulation, the finite element-extended boundary condition formulation employs the spherical vector harmonic functions as the global testing functions and utilizes the orthogonality of these harmonic functions. As a result, this approach can be more efficient than the previous one. Its internal field computation, however, may be less accurate than that of the finite element-surface integral method. In all three formulations, one encounters the problem of dealing with the matrix equations (8.13) and (8.14). Generating these matrices and then solving them are the two major aspects of the problem. Also, selecting a set of known bases for {Es,} to compute the coupled bases {Hs,} is another important issue in the second (with subdomain method) and third approaches. These could be the subjects for future studies.

CHAPTER IX CONCLUSION AND FUTURE WORK In this thesis several finite element-boundary element methods were developed for solving the problem of open region electromagnetic scattering. 9.1. Conclusion A finite element-boundary element method was developed in Chapter 4 for the numerical analysis of the scattering by cylindrical objects. The method overcomes two major shortcomings which existed in the previous methods. First, it relaxes the restriction that the artificial boundaries be circles to the condition that they follow the shape of the scatterer and thus substantially reduces the solution domain. Second, with a modification on the coupling procedure it results in a highly sparse or uniformly banded system matrix which can be efficiently solved using special algorithms. Numerical results were presented, showing the high accuracy of the method as well as its notable capability and versatility. Through a comparison with other numerical methods, it is concluded that this method provides one of the most accurate and efficient means to solve complex two-dimensional scattering problems involving inhomogeneous media. We need to note that, though the computer code was developed based on the formulation using two artificial boundaries, the alternative formulation proposed 142

143 in Section 4.2.5 may turn out to be theoretically more attractive since it can be interpreted using the equivalence principle [Harrington, 1961]. Nevertheless, from the computational aspect both formulations have the same efficiency. The finite element-boundary element method was then extended in Chapter 5, by incorporating the physical optics approximation, to treat the scattering by a coated wedge and half-plane. Numerical results were presented for the transverse magnetic case to show the surface field and surface impedance behavior around the edge. It is demonstrated that the impedance boundary condition is invalid at a point close to the edge, though it provides a very good approximation when the point is about a tenth of a wavelength away from the edge. In Chapter 6, a finite element-boundary element method was developed and applied to characterizing the scattering and transmission properties of an inhomogeneously filled aperture in a thick conducting plane. Of particular interest in this chapter is the use of the fast Fourier transform to evaluate the boundary integrals and the use of the conjugate gradient method to solve the system of equations. The high accuracy, capability, and versatility were demonstrated by the numerical results. Some of the results, such as those for the truncated strip gratings embedded in dielectric layers, would be impossible to obtain using the conventional techniques. It is felt that this technique is a promising approach for treating large and complex radiation and scattering problems. In Chapter 7, the application of isoparametric elements in the finite elementboundary element method was described in order to more accurately model the curved and arbitrarily shaped scatterers. The use of the isoparametric elements can further enhance the efficiency and accuracy of the method. Chapter 8 presented and discussed several hybrid techniques which combine the finite element method with either a surface integral equation or an expansion of vector eigenfunctions for three-dimensional scattering. Some are just the extensions

144 of the previously developed methods for the two-dimensional case or for the case of a body of revolution. 9.2. Future Work An obvious, interesting, challenging, and very useful as well as very difficult problem is three-dimensional scattering from arbitrarily shaped electrically and magnetically permeable bodies which may contain conducting surfaces. This is a subject which is currently attracting more and more people and will become a popular field in the near future. The finite element-surface integral equation method, proposed in Section 8.3.1, is felt to be a promising approach for attacking the problem when it incorporates the conjugate gradient method and the fast Fourier transform as described in Chapter 6. While the theoretical formulation is available, its numerical implementation, which needs a tremendous effort, is the most difficult. The implementation of the finite element-extended boundary condition method, described in Section 8.3.3, is also very worthwhile since it could be more efficient for the far field computation. As a special three-dimensional problem, the scattering by an inhomogeneously filled finite aperture in a thick conducting plane, in contrast to the two-dimensional one treated in Chapter 6, is also a problem of practical interest. The formulation given in Section 6.2.1 is applicable; however, a numerical implementation is necessary. The finite element-boundary element method combined with the physical optics approximation, as described in Chapter 5 for coated wedges and half-planes, can be applied to treat the scattering by very large two-dimensional geometries. While the physical optics approximation (and surface wave formulation, if necessary) can be employed to approximate the fields on the smooth portion, the finite element method and boundary integral equation must be used where there is a surface discontinuity or

145 coating inhomogeneity. The incorporation of the physical optics approximation and surface wave formulation can remove the problem of singular frequencies introduced by the employment of the surface or boundary integral equation.

BIBLIOGRAPHY 146

147 BIBLIOGRAPHY [1] Auckland, D. T. and R. F. Harrington, "Electromagnetic transmission through a filled slit in a conducting plane of finite thickness, TE case," IEEE Trans. Microwave Theory and Tech., vol. MTT-26, pp. 499-505, July 1978. [2] Auckland, D. T. and R. F. Harrington, "A nonmodal formulation for electromagnetic transmission through a filled slot of arbitrary cross section in a thick conducting screen," IEEE Trans. Microwave Theory and Tech., vol. MTT-28, pp. 548-555, June 1980. [3] Axelsson, 0. and V. A. Barker, Finite Element Solution of Boundary Value Problems - Theory and Computation. London: Academic Press, 1984. [4] Bowman, J. J., T. B. A. Senior, and P. L. E. Uslenghi ed., Electromagnetic and Acoustic Scattering by Simple Shapes. North-Holland Publishing Company-Amsterdam, 1969. [5] Boyes, W. E. and A. A. Seidl, "An efficient hybrid finite element method for densely coated conductors," Program and Abstracts of 1989 URSI Radio Science Meeting, p. 333, San Jose, CA, June 1989. [6] Bucci, D. M., and G. Franceschetti, "Electromagnetic scattering by half plane with two face impedances," Radio Sci., vol. 11, pp. 49-59, 1976. [7] Burnside, W. D., C. L. Yu, and R. J. Marhefka, "A technique to combine the geometrical theory of diffraction and the moment method," IEEE Trans. Antennas Propagat., vol. AP-23, pp. 551-558, July 1975. [8] Chang, S. K. and K. K. Mei, "Application of the unimoment method to electromagnetic scattering of dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-24, pp. 35-42, Jan. 1976. [9] Collins, J. D. and J. L. Volakis, "A boundary integral conjugate gradient fast Fourier transform method for solving two-dimensional scattering problems," Program and Abstracts of 1989 URSI Radio Science Meeting, p. 217, San Jose, CA, June 1989. [10] George, A. and J. Liu, Computer Solutions of Large Sparse Positive Definite Systems. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1981.

148 [11] Guo, Y.-N., "Finite-element solution of electromagnetic scattering from inhomogeneous dielectric cylinder with arbitrary cross section," M.S. thesis, Nanjing University, 1985. [12] Harrington, R. F., Time-Harmonic Electromagnetic Fields. New York: McGraw-Hill, 1961. [13] Harrington, R. F., Field Computation by Moment Methods. New York: Macmillan, 1968. [14] Hestenes, M. R. and E. Steifel, "Method of conjugate gradients for solving linear systems," J. Res. Nat. Bur. Standard, vol. 49, pp. 409-436, Dec. 1952. [15] Hongo, K. and G. Ishii, "Diffraction of an electromagnetic plane wave by a thick slit," IEEE Trans. Antennas Propagat., vol. AP-26, pp. 494-499, May 1978. [16] Irons, B. M., "A frontal solution program for finite element analysis," lnt. J. Num. Meth. Eng., vol. 2, pp. 5-32, 1970. [17] Jeng, S. K. and C. H. Chen, "On variational electromagnetics: theory and application," IEEE Trans. Antennas Propagat., vol. AP-32, pp. 902-907, Sept. 1984. [18] Jeng, S. K., "Aperture admittance matrix by finite element method for scattering from a cavity-backed aperture," 1988 IEEE AP-S International Symposium Digest, vol. 3, pp. 1134-1137, Syracuse, NY, June 1988. [19] Jin, J.-M. and V. 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. 1988a. [20] Jin, J.-M. and V. V. Liepa, "A note on hybrid finite element method for solving scattering problems," IEEE Trans. Antennas Propagat., vol. AP-36, no. 10, pp. 1486-1490, Oct. 1988b. [21] Jin, J.-M., V. V. Liepa, and C.-T. Tai, "A volume-surface integral equation for electromagnetic scattering by inhomogeneous cylinders," The Journal of Electromagnetic Waves and Applications, vol. 2, no. 5/6, pp. 573-588, 1988c. [22] Jin, J.-M. and V. V. Liepa, "A numerical technique for computing I'TM scattering by coated wedges and half-planes," Electromagnetics, vol. 9, no. 2, pp. 201-214, 1989a. [23] Jin, J.-M. and V. V. Liepa, "Simple moment method program for computing scattering from complex cylindrical obstacles," Proc. Inst. Elec. Eng., part H, vol. 136, no. 4, pp. 321-329, 1989b.

149 [24] Jin J.-M. and J. L. Volakis, "New technique for characterizing diffraction by inhomogeneously filled slot of arbitrary cross section in thick conducting plane," Electronics Letters, vol. 25, no. 17, pp. 1121-1123, 17th August 1989c. [25] Jin J.-M. and J. L. Volakis, "Electromagnetic scattering by a perfectly conducting patch array on a dielectric slab," IEEE Trans. Antennas Propagat., vol. AP-38, 1990. [26] Kashyap, S. C. and M. A. K. Hamid, "Diffraction characteristics of a slit in a thick conducting screen," IEEE Trans. Antennas Propagat., vol. AP-19, pp. 499-507, July 1971. [27] Knott, E. F. and T. B. A. Senior, "Non-specular radar cross section study," The University of Michigan Radiation Laboratory Report No. 011764-1 -T(AFAL-TR-73-422), Jan. 1974. [28] Lehman, G. W., "Diffraction of electromagnetic waves by planar dielectric structures. I. Transverse electric excitation," Journal of Math. Physics, vol. 11, pp. 1522-1535, May 1970. [29] Maliuzhnets, G. D., "Excitation, reflection and emission of surface waves from a wedge with given face impedance," Sov. Phys. Dokl., vol. 3, pp. 752-755, 1958. [30] Marin, S. P., "Computing scattering amplitudes for arbitrary cylinders under incident plane waves," IEEE Trans. Antennas Propagat., vol. AP30, pp. 1045-1049, Nov. 1982. [31] Mason, J. L., "Finite element solution for electromagnetic scattering from two-dimensional bodies," Ph.D. dissertation, The University of Michigan, Ann Arbor, 1982. [32] McCartin, B. J., G. Meltz, R. Mittra, and L. J. Bahrmasel, "Application of the control region approximation in conjunction with absorbing boundary conditions to the direct solution of electromagnetic scattering problems," Program and Abstracts of 1988 URSI Radio Science Meeting, p. 14, Syracuse, NY, June 1988. [33] McDonald, B. H. and A. Wexler, "Finite-element solution of unbounded field problems," IEEE Trans. Microwave Theory Tech., vol. MTT-20, pp. 841-847, Dec. 1972. [34] McDonald, B. H. and A. Wexler, Finite Element in Electrical and Magnetic Field Problem (Eds. M. V. K. Chari and P. P. Silvester). Wiley, London and New York, 1980, chap. 9.

150 [35] Mei, K. K. and J. Van Bladel, "Scattering by perfectly conducting rectangular cylinder," IEEE Trans. Antennas Propagat., vol. AP-11, pp. 185-192, 1963. [36] Mei, K. K., "Unimoment method of solving antenna and scattering problems," IEEE Trans. Antennas Propagat., vol. AP-22, pp. 760-766, Nov. 1974. [37] Mikhlin, S. G., Variational Methods in Mathematical Physics. Macmillan, New York, 1964. [38] Mittra, R., "A critical look at absorbing boundary conditions for partial differential equations arising in electromagnetic scattering problems," Program and Abstracts of 1988 URSI Radio Science Meeting, p. 15, Syracuse, NY, June 1988. [39] Morgan, M. A. and K. K. Mei, "Finite-element computation of scattering by inhomogeneous penetrable bodies of revolution," IEEE Trans. Antennas Propagat., vol. AP-27, pp. 202-214, March 1979. [40] Morita, N., "Diffraction by arbitrary cross-sectional semi-infinite conductor," IEEE Trans. Antennas Propagat., vol. AP-19, pp. 358-364, 1971. [41] Morgan, M. A., C. H. Chen, S. C. Hill, and P. W. Barber, "Finite elementboundary integral formulation for electromagnetic scattering," Wave Motion, vol. 6, no. 1, pp. 91-103, Jan. 1984. [42] Neerhoff, F. L. and G. Mur, "Diffraction of a plane electromagnetic wave by a slit in a thick screen placed between two different media," Appl. Sci. Res., vol. 28, pp. 73-88, July 1973. [43] Newman, E. H., "TM scattering by a dielectric cylinder in the presence of a half-plane," IEEE Trans. Antennas Propagat., vol. AP-33, pp. 773-782, 1985. [44] Newman, E. H., "TM and TE scattering by a dielectric/ferrite cylinder in the presence of a half-plane," IEEE Trans. Antennas Propagat., vol. AP-34, pp. 804-813, 1986. [45] Orikasa, T., S. Washisu, T. Honma, and I. Fukai, "Finite element method for unbounded field problems and application to two-dimensional taper," Int. J. Num. Meth. Eng., Vol. 19, pp. 157-168, 1983. [46] Richmond, J. H., "Scattering by a dielectric cylinder of arbitrary crosssection shape," IEEE Trans. Antennas Propagat., vol. AP-13, pp. 334 -341, 1965. [47] Richmond, J. H., "TE-wave scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Antennas Propagat., vol. AP-14, pp. 460-464, July 1966.

151 [48] Ricoy, M. A. 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, March 1989. [49] Ruck, G. T., D. E. Barrick, W. D. Stuart, and C. K. Krichbaum, Radar Cross Section Handbook. (Ed. G. T. Ruck) Plenum Press, New York, 1970. [50] Schaubert, D. H., D. R. Wilton, and A. W. Glisson, "A tetrahedral modelling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies," IEEE Trans. Antennas Propagat., vol. AP-32, pp. 77-85, Jan. 1984. [51] Senior, T. B. A., "Diffraction by a semi-infinite metallic sheet," Proc. Roy. Soc., vol. 213A, pp. 436-458, 1952. [52] Sewell, G., Analysis of Finite Element Method: PDE/PROTRAN. SpringerVerlag New York Inc., 1985. [53] Silvester, P. P. and M. S. Hsieh, "Finite-element solution of 2-dimensional exterior field problems," Proc. Inst. Elec. Eng., vol. 118, pp. 1743-1747, Dec. 1971. [54] Silvester, P. P. and R. L. Ferrari, Finite Elements for Electrical Engineers. Cambridge University Press, Cambridge, 1983. [55] Stratton, J. A., Electromagnetic Theory. (Chapter 7) McGraw-Hill, New York, 1941. [56] Tai, C. T., Dyadic Green's Functions in Electromagnetic Theory, International Textbook Company, Scranton, Pennsylvania, 1971. [57] Tang, C. H., "Backscattering from dielectric-coated infinite cylindrical obstacles," J. of Appl. Physics, vol. 28, no. 5, pp. 628-633, May 1957. [58] Tiberio, R., G. Pelosi, and G. Manara, "A uniform GTD formulation for the diffraction by a wedge with impedance faces," IEEE Trans. Antennas Propagat., vol. AP-33, pp. 867-873, 1985. [59] Volakis, J. L., "A uniform geometrical theory of diffraction for an imperfectly conducting half-plane," IEEE Trans. Antennas Propagat., vol. AP-34, pp. 172-180, 1986. [60] Wang, N., "Electromagnetic scattering from a filled slit in a thick conducting screen, TM case," The Ohio State Univ. ElectroScience Lab. Tech. Report 714349-8, May 1984. [61] Washisu, S., I. Fukai, and M. Suzuki, "Extension of finite-element method to unbounded field problems," Electron. Letts., vol. 15, pp. 772-774, 1979.

152 [62] Waterman, P. C., "Matrix formulation of electromagnetic scattering," Proc. IEEE, vol. 53, pp. 805-812, 1965. [63] Waterman, P. C., "Matrix methods in potential theory and electromagnetic scattering," J. Appl. Phys., vol. 50, pp. 4550-4566, 1979. [64] Wilton, D. R. and W. F. Richards, "Application of near-field radiation conditions to the solution of electromagnetic scattering problems," Program and Abstracts of 1988 URSI Radio Science Meeting, p. 13, Syracuse, NY, June, 1988. [65] Wu, K.-L., G. Y. Delisle, and D.-G. Fang, "Electromagnetic scattering from an arbitrary dielectric coated conducting cylinder by coupled finite - boundary element method," 1988 IEEE AP-S International Symposium Digest, vol. I, pp. 380-383, Syracuse, NY, June 1988. [66] Wu, T. K. and L. L. Tsai, "Scattering by arbitrarily cross sectioned layered lossy dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-25, pp. 518-524, 1977a. [67] Wu, T. K. and L. L. Tsai, "Scattering by a dielectric wedge: a numerical solution," IEEE Trans. Antennas Propagat., vol. AP-25, pp. 570-571, 1977b. [68] Yuan, X., D. R. Lynch, and J. W. Strohbehn, "A new hybrid moment/finite element method for electromagnetic scattering from arbitrarily shaped inhomogeneous dielectric objects," Program and Abstracts of 1989 URSI Radio Science Meeting, p. 332, San Jose, CA, June 1989. [69] Zienkiewicz, 0. C., The Finite Element Method, 3rd ed. New York: McGraw-Hill, 1977.