FINITE ELEMENT SOLUTION FOR ELECTROMAGNETIC SCATTERING FROM TWO-DIMENSIONAL BODIES by John Lawrance Mason A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1982 Doctoral Commi ttee: Professor Thomas B.A. Senior (Co-Chairman) Professor William J. Anderson(Co-chairman) Associate Professor Valdis V. Liepa Professor Ronald J. Lomax Professor Herschel Weil (Aerospace Engineering) RL-734 = RL-734

For Linda, Jane, and David -11 -

ACKNOWLEDGEMENTS The author wishes to acknowledge the financial support provided in the early stages of the research by the Radiation Laboratory of the Department of Electrical and Computer Engineering at The University of Michigan (a stipend) and by Western Michigan University, his employer (a sabbatical leave). At Western, Dean W. Chester Fitch and Assistant Dean Robert E. Boughner of the College of Applied Sciences, and Professor Cassius A. Hesselberth, Chairman of the Electrical Engineering Department were particularly helpful in obtaining support. The author is grateful for the aid and advice given him by the Doctoral Committee. He also received advice from the following individuals: Professor C. Bruce Sharpe of The University of Michigan, and Professor Paul J. Eenigenburg and Assistant Professor Dean R. Johnson, both of Western Michigan University. Assistance with the computer program was received from the staff at the Western Michigan University Computer Center, particularly from William Granet, George R. Kohrman, and Brian T. Mitchell. The author wishes to thank Linda Mason and Wanita Rasey for typing the first and final versions of the dissertation respectively and for their editorial assistance. Also, the author is grateful for the overnight accommodations provided by Leroy and Sara Kottke on numerous occasions. Finally, the author wishes to thank his wife and children for their patience and understanding during the time this work was in progress. -iii

TABLE OF CONTENTS Page DEDICATION i ACKNOWLEDGEMENTS - LIST OF ILLUSTRATIONS vi LIST OF TABLES ix LIST OF APPENDICES x CHAPTER I. INTRODUCTION 1 1.1 Historical Background 1 1.2 Formulation of Scattering Problems 2 1.3 Previous Work 4 1.4 Objectives of the Present Study 7 CHAPTER II. FORMULATION OF THE SCATTERING PROBLEM 10 2.1 Statement of the Problem 10 2.2 Integral Equation Formulation 14 2.3 Differential Equation Formulation 14 2.4 Far Scattered Field and the Radar Cross Section 19 CHAPTER III. SOLUTION BY THE FINITE ELEMENT METHOD 24 3.1 Summary of the Problem 24 3.2 Derivation of an Integral Statement of the Problem 25 3.3 Interpolation Functions and Expansion of the Unknown 28 3.4 Conversion to a Set of Algebraic Equations 29 3.5 The Element Mesh 34 3.6 Quadrilateral Finite Elements 36 3.7 Triangular Elements at the Strip Edge 41 CHAPTER IV. SOLUTION PROCEDURE 46 4.1 The Computer Program 46 4.2 Element Matrix Coefficients for Quadrilateral Elements 47 4.3 Element Matrix Coefficients for Triangular Elements 55 -1 v

Page 4.4 Incident Field Terms 64 4.5 Induced Current 66 4.6 Radar Cross Section 71 CHAPTER V. RESULTS OF THE NUMERICAL SOLUTION 80 5.1 Error Criteria for Numerical Results 80 5.2 Selection of Element Mesh Parameters 80 5.3 Effects of Element Density 84 5.4 Effects of Strip Resistivity 89 5.5 Solution Cost Considerations 98 CHAPTER VI. VOLUME SCATTERING EXAMPLE 105 6.1 Statement and Formulation of Example 105 6.2 Selection of an Incident Field 109 6.3 Implementation of the Solution Procedure 111 6.4 Numerical Results 112 CHAPTER VII. CONCLUSIONS AND RECOMMENDATIONS 117 7.1 Overview of the Work 117 7.2 Limitations of the Method 119 7.3 Recommendations for Further Work 120 APPENDICES 121 BIBLIOGRAPHY 1 54

LIST OF ILLUSTRATIONS Figure Page 2.1 A strip lying in the x-z plane. 11 2.2 An axial view of the strip. 13 2.3 A diagram for the far field computation. 22 3.1 The solution region. 26 3.2 The element mesh. 36 3.3 A quadrilateral finite element. 38 3.4 A square element in normalized coordinates. 39 3.5 Triangular elements in mesh at strip edge. 42 3.6 A general triangular element. 43 3.7 A triangular element in normalized coordinates. 44 4.1 The quadrilateral element adjacent to the strip. 50 4.2 A quadrilateral element adjacent to the contour rB. 53 4.3 A triangular element adjacent to the strip at x = w/2. 58 4.4 A triangular element adjacent to the strip at x = -w/2. 60 5.1 Error in IP(4,P0 )I vs. element density, (M )(Xo/w) = (M)(x /Pp ), R =. = 1.0, w/Xo = 0.5 = 05 backscattering (q = p ). 85 5.2 Error in arg P((,qo) vs. element density, (M )(x /w) = (M )(X /p ), R = 1.0, w/X = 0.5, p /x = 0.5 y 0 1 0 10 backscattering (8 = o). 86 0 5.3 Error in |P(,(o) )I vs. f for R = 1.0, w/X = 0.5, p /X 0.5 backscattering (q = ). 87 -vi -

Figure Page 5.4 Error in arg P({,( ) vs. ( for R = 1.0, w/x = 0.5, p /x 0.5, backscattering (q = O ). 88 5.5 Error in IP(P,(o)I vs. ( for R= 0.0, w/x = 0.5, p /X A 0.5, backscattering ( = () ). 90 0 ' 5.6 Error in arg P(,(O) o) vs. ( for R = 0.0, w/xo = 0.5, p /X - 0.5, backscattering (q = o) 91 5.7 Error in lP({,b o)I vs. (1/2)k w cos ( for R = 1.0 (M )(x /w) 16, (M )(x /p ) 16, backscattering 0 y o ( = ) 92 5.8 Error in arg P(,( 0o) vs. (1/2)kow cos ( for R = 1.0, (Mx)(xo/w) 16, (My)(x /p1) 16, backscattering ( = ). 93 5.9 Error in IP((, q)I vs. ( for R = 0.1, w/X0 = 0.5, p / 0.5, (M )(x /p ) 17, backscattering 1 o y 0 1 ( =o ). 96 5.10 Error in arg P(q(, ) vs. ( for R = 0.1, w/X = 0.5, p /xo 0.5 (M )(X /p ) 17, backscattering ( = P)o. 97 6.1 The extended solution region. 106 6.2 c/xo vs. ( for bistatic scattering, o = 0~ R = 1.0, w/x = 0.25. 113 6.3 Arg P((,( 0) vs. ( for bistatic scattering, qO = 0~, R = 1.0, w/X = 0.25. 114 6.4 o/X vs. (P for bistatic scattering, R = 1.0, w/x = 0.25, dielectric cylinder (R = 2) radius/x = 0.25. 115 0 -vii

Figure Page 6.5 Arg P(G,) o) vs. for bistatic scattering R = 1.0, w/Xo = 0.25, dielectric cylinder (ER = 2) radius/x = 0.25. 116 A.1 M vs. (EL)M. 122 y x A.2 A vs. SMEC. 123 B.1 Scattered field magnitude vs. distance outward from the strip center ( 0 = 90~, w/x = 0.5, M = 6, x = 0). 1 027 B.2 Scattered field phase vs. distance outward from the strip center (o = 90~, w/X = 0.5, = 1.0, M = 6, x = 0). 129 X C.1 Finite elements. 134 C.2 Transmission line. 143 C.3 Symmetrical T network. 145 C.4 Terminated transmission line. 147 C.5 IESCI vs. y/x for resistive sheet (R = 1.0) Z 0 j 0 = 1.0 and "free space termination" (KII given by C.33). 150 C.6 IEsci vs. y/X for resistive sheet (R = 1.0), = 1.0 and "characteristic impedance termination" o (KI given by C.57). 151 C.7 Arg Esc vs. y/x for resistive sheet (R = 1.0), zo ' A = 1.0 and "free space termination" (KII given by C.33). 153 -viii

LIST OF TABLES Table Page 5.1 Error in P(,(0o) for Backscattering (q = 0o) (Mx)(o/w) = 18, (My)(o/p1) - 17, w/ = 0.5, p /A 0 5 94 5.2 Variable Memory Requirement for a Resistive Strip (M)(o/w) = (My)(xo/p) 101 5.3 Global Matrix Memory for a Resistive Strip for (Mx)(o /w) = 16 102 5.4 Computation Time for a Resistive Sheet (Mx /w = Myx /p = 16) 104 xo yo 1 -ix

LIST OF APPENDICES Page APPENDIX A. APPENDIX B. APPROXIMATE RELATIONSHIP BETWEEN ELEMENT MESH PARAMETERS TRANSMISSION LINE INTERPRETATION FOR NUMERICAL RESULTS B.1 The Finite Element Mesh as a Transmission Line B.2 Effect of Standing Waves APPENDIX C. COMPUTATION OF SCATTERED FIELD FROM PLANE WAVE NORMAL INCIDENCE ON INFINITE RESISTIVE SHEET C.1 Introduction C.2 Finite Element Formulation C.3 Element Matrix Coefficients C.4 Set of Linear Algebraic Equations C.5 Transmission Line Interpretation 121 125 125 128 131 131 132 138 140 141

November 1982 ABSTRACT FINITE ELEMENT SOLUTION FOR ELECTROMAGNETIC SCATTERING FROM TWO-DIMENSIONAL BODIES by John Lawrance Mason Co-Chairmen: Thomas B.A. Senior, William J. Anderson The problem of two-dimensional electromagnetic scattering from bodies in free space has been formulated using a differential equation finite element method. An incident plane wave polarized so as to excite only E-wave fields has been assumed. A weak integral statement was obtained using an approach equivalent to the weighted residuals method, the weighting functions being selected by Galerkin's method. The finite element mesh of general first-order quadrilaterals extended outward into the far field region of the scattering body, where the boundary condition at the outer boundary was obtained from the asymptotic expression for the scattered field. In order to compare the finite element method with an integral equation moment method, a special case was solved numerically, namely that of a thin, infinitely long resistive strip. Accurate numerical results from a moment method solution were available to serve as a check on the accuracy of the finite element method results. For a perfectly conducting strip (the limiting case of a resistive strip) the finite element mesh contained special triangular elements at the edges of the strip to model the singulari ties in the magnetic field -8 -

which occur there. Two important disadvantages of the finite element solution became apparent. First, the boundary condition used at the outer boundary leads to the presence of standing waves in the computed scattered field, a non-physical result. Second, the computer memory required increased much more rapidly with strip width than in the moment method solution. The finite element method used was thus unable to compete with the moment method solution for the thin strip problem. However, to demonstrate a problem for which a moment method solution may be less attractive because of its increased memory requirements, numerical results were obtained using the finite element method for scattering from a body consisting of a thin resistive strip embedded in a dielectric cylinder. -9 -

CHAPTER I. INTRODUCTION 1.1 Historical Background The scattering of electromagnetic waves from bodies consisting of a simple shape of homogeneous material has been studied for years. Much of the early work focussed on bodies whose geometry coincided with coordinate surfaces of systems for which the wave equation is separable. Typically, modal expansions of the fields were written, followed by application of boundary conditions. The solution was in the form of an infinite series. Numerical values of the fields were then obtained by summing a finite number of terms of the series, which was generally a time-consuming process. The solution of scattering problems using integral equation formulations was not accomplished except in two cases. In these cases the body surfaces were complete coordinate surfaces of a system where the wave equation is separable. An integral transform was found to invert the equation for a solution. The techniques mentioned above are usually called analytical techniques, even though they are often followed by a numerical procedure to evaluate the resulting expressions at particular points. For bodies having complicated shapes and/or inhomogeneous material properties, these techniques fail. Instead, so-called numerical techniques are used. Until about 1970, most electromagnetic scattering problems which were solved numerically involved finite -1 -

-2 - difference, variational, or moment methods (Wexler, 1969; Jones, 1956; Harrington, 1968). Since 1970 the finite element method, utilized earlier in applied mechanics, has been applied to scattering problems with some success. All of these numerical methods provide a procedure for converting a scattering problem, formulated in terms of either a differential equation or an integral equation, into a problem of solving a system of linear algebraic equations. 1.2 Formulation of Scattering Problems For numerical solutions of scattering problems involving bodies in unbounded space, integral equation formulations have been preferred to those of differential equation formulations. The solution region for the integral equation is confined to the surface of the scattering body with boundary and radiation conditions automatically taken into account. Integral equation formulations are well suited for scattering by thin sheets and by bodies into which the field does not penetrate. However, it is more difficult to obtain such formulations for scattering problems involving thick, inhomogeneous bodies into which the field does penetrate. The most commonly used numerical technique for the solution of integral equations arising from scattering problems is the method of moments. The resultant system of linear algebraic equations has a dense matrix of coefficients. In view of the current computer hardware available, this method can be applied to scattering bodies with surface areas up to about one square wavelength or, in the case of two-dimensional problems, to lengths of about five

-3 - wavelengths. A characteristic of the method is that the solution becomes non-unique at frequencies for which the interior of the scattering body is resonant (Klein and Mittra, 1975). Differential equation formulations have been used primarily for interior scattering problems when numerical solutions have been made. For exterior problems, the solution region for the differential equation is unbounded, requiring application of the radiation condition at infinity. The solution region is thus larger than that for an integral equation formulation and is of higher dimension. However, a differential equation formulation is convenient for inhomogeneous media. Interior problems formulated in terms of a differential equation are usually solved numerically with finite difference, moment, or finite element methods. The finite difference algorithms are usually the simplest. However, the finite element method has advantages over both the other methods in terms of accuracy and ability to handle bodies having complex shapes. All of these methods lead to systems of linear algebraic equations having sparse matrices. Proper arrangement of the equations results in a banded matrix allowing use of efficient equation solvers so that computer time and memory requirements can be reduced. These comments apply to exterior problems for which, however, all the methods require application of the radiation condition and, in general, require a large system of equations due to the large solution region. Difficulties associated with these two requirements for an

-4 - exterior problem may easily overshadow the advantages of a differential equation formulation over an integral equation formulation, except for the ability to handle inhomogeneous media. The work reported here is an effort to determine if, and under what conditions, a finite element differential equation solution might be more economical than a moment-method integral equation solution for exterior scattering by electrically large bodies (dimensions of one half wavelength or greater). Also a finite element solution is demonstrated for a case involving a thick body for which a moment method integral equation solution may not be practical. 1.3 Previous Work The finite element method has seen limited application in solution of integral equations (for example, Jeng and Wexler, 1977, and Sankar and Tong, 1975) because of the efficiency of the method of moments. The finite element method was first applied successfully to interior boundary value problems because of the difficulty in handling the radiation condition in exterior problems. For example, Arlett et al., (1968) solved the Helmholtz equation in a waveguide using a variational formulation combined with a finite element discretization. Silvester and Hsieh (1971) reported solving Laplace's equation for an exterior problem using a variational formulation. They divided the solution region into two parts, one part being a local region containing, for example, sources, conductors and dielectrics. The local region was treated using finite elements.

-5 - The other region, extending to infinity was treated as a single exterior element. In the portion of the formulation related to this element it was necessary to require that the stored energy in the exterior region be finite, a requirement which is met by solutions of Laplace's equation. Determination of the element matrix involved the solution of an integral equation. McDonald and Wexler (.1972) also divided the solution region into a local region and external region. An integral equation using a free-space Green's function was used to relate the field at points within the local regions to points on the boundary of the external region. Examples of cases satisfying the Laplace and Helmholtz equations were solved using a variational formulation. As a result of the use of the integral equation, the global matrix obtained was more dense than the typical banded matrix of a finite element solution. The works of Silvester and Hsieh, and McDonald and Wexler are examples of the combined use of finite element and boundary integral methods. Another example is the study by Berkhoff (1975) on linear water wave propagation. A survey of work using this combined approach was given by Zienkiewicz, Kelly and Bettess (1977). The technique was also discussed by Brebbia and Walker (1980), Zienkiewicz (1977), and McDonald and Wexler (.1980). In another approach known as the unimoment method and developed by Mei (.1974), the fields in the exterior region are expressed in a modal expansion and the interior region is handled using finite element or finite difference techniques. The interior

-6 - and exterior problems are then coupled together by enforcing boundary conditions for the fields at the boundary between the regions. The resulting system of equations has a banded matrix. The unimoment method has been applied to the scattering of dielectric cylinders (Chang and Mei, 1976) and to inhomogeneous penetrable bodies of revolution (Morgan and Mei, 1979). Another way of handling the outer boundary condition for an exterior problem, proposed by Ungless (1973) and applied by Bettess (1977), involves the use of elements extending to infinity. The shape functions for such infinite elements must realistically model the behavior of the unknown as the distance increases from the local region of interest and should lead to integrals over the element domain which are finite. Bettess and Zienkiewicz (.1977) used a combination of finite and infinite elements in a study of water waves. Included in the infinite element shape functions was an exponential factor with a decay length. Although this method does not require detailed knowledge of the form of the solution in the outer region, some knowledge of the exact solution is required to be able to select the decay length. This approach does not give an accurate indication of the behavior of the unknown towards infinity; however, the effect of the far region on the local region of interest is introduced. Thatcher (1978) used an infinite number of elements defined in a systematic way to solve Laplace's equation in an unbounded region. This method will also handle singularities occurring at

-7 - boundaries of the solution region. Application of this method to the solution of the Helmholtz equation in unbounded regions has apparently not been made. Other work has been reported on absorbing or non-reflecting boundary conditions for use with finite difference or finite element solution of wave propagation problems, whereby the exterior problem is converted to a pseudo-interior problem. Smith (1974) proposed the superposition of the solutions of Dirichlet and Neumann problems to completely eliminate reflections at a plane boundary of a half space. However, this approach does not seem applicable to other geometries. Orlanski (1976) used the Sommerfeld radiation condition and a numerical evaluation of phase velocity in the vicinity of the boundary. In a mechanics problem involving determination of strain, Hanson and Petschek (1976) used a spring-dashpot system as a terminating network on a one-dimensional system to absorb a single wavelength. They reported successfully using a similar technique on a two-dimensional problem. Engquist and Majda (1977) developed "perfectly absorbing" boundary conditions for general classes of wave equations. They then derived local boundary conditions, suitable for numerical calculations, which approximate the theoretical boundary conditions. 1.4 Objectives of the Present Study The first objective of the present study is to develop a finite element formulation for an exterior scattering problem

-8 - Specifically, a differential equation formulation is used. The variational method, typical of finite element work, is replaced by a Green's identity and Galerkin's method. This amounts to using the method of weighted residuals. Also, the handling of the boundary condition at the edge of the finite element region is performed in such a manner that: (1) the influence of the boundary condition is made through individual element contributions and (2) the band nature of the global system matrix is preserved. This rules out the use of boundary integral techniques. The second objective is to apply the finite element formulation to the solution of a specific scattering problem to determine if it has an advantage over a moment method integral equation solution for electrically large bodies. Scattering from an infinitely long, finite width, infinitesimally thin, flat resistive strip is considered since an integral equation formulation and accurate moment method numerical solution (Knott, Liepa and Senior, 1973) are available. Briefly a resistive sheet has the properties that: (1) the component of the electric field tangent to the surface has no discontinuity at the surface, and (2) the component of the magnetic field tangent to the surface has a discontinuity at the surface which is directly proportional to the tangential electric field. A perfectly conducting sheet can be considered as a resistive strip with zero resistivity. Both formulations are capable of handling a strip whose resistivity varies in the direction perpendicular to the axis of the strip. Further,

-9 - since a singularity in the magnetic field occurs at the edge of a perfectly conducting strip, the finite element discretization includes special elements to accurately model this behavior. The practical requirement of each solution method is to compute the far scattered field and the radar cross section. The effect of the element density on the accuracy of these results for various strip widths and resistivities, and for different angles of incidence is investigated. In order to accomplish this, an automatic element mesh generator is included in the computer program so that the element density, the length-to-width ratio of the elements, and the distance outward to which the finite element region extends can be easily changed. The effect of the element density on the computer memory requirements and CPU time is also examined. The third objective is to apply the finite element formulation to the solution of scattering from a thick body into which the field penetrates so as to demonstrate a case where the moment method integral equation approach may not be practical. The specific case of a resistive strip embedded in a dielectric cylinder is considered.

CHAPTER II. FORMULATION OF THE SCATTERING PROBLEM 2.1 Statement of the Problem An infinitely long, infinitesimally thin strip of width w lying parallel to the z-axis and in the x-z plane is immersed in an infinite, uniform medium with permeability P0 and permittivity Eo (Fig. 2.1). The strip is a resistive sheet whose resistivity in ohms per square is RZ, where Z = v /7E. The sheet is defined by its boundary conditions [n ^ E] = 0 (2.1) and RZ n ^[H] = -n A (n A E), (2.2) where n is a unit normal to the strip in the positive y direction and the notation [ ]+ represents a discontinuity across the strip. E and H are, respectively, the total electric and magnetic field strengths. The normalized resistivity R may be a function of x, i.e., R = R(x). A plane wave, whose direction of propagation is perpendicular to the z-axis, is incident on the strip. The polarization is such that i ^ A -ik (x cos cp + y sin f ) E = Ezz = z e (2.3) -10 -

-1l y I I 0 x = -w/2 ff IL x =7 / x z x = w/2 Fig. 2.1: A strip lying in the x-z plane.

-12 - where o0 is the angle that the direction of propagation makes with the plane of the strip (Fig. 2.2). An e -Wt time dependence is assumed for all time-varying quantities. The incident fields E and Hi are the fields which would exist everywhere if the strip was not present. In the presence of the strip, the total fields are E and H which are the sums of the incident and scattered fields: E = Ei + Esc (2.4) and H = H1 + Hsc (2.5) The scattered fields are associated with an outward radiation of energy and must satisfy a radiation condition at infinity. A surface current induced on the strip is given by K = n [H], (2.6) but since the incident field has no discontinuity at the strip, K = n [. (2.7) If a knowledge of the surface current density on the scattering body can be obtained, the far scattered field and the radar cross section can be determined using an integral representation such as that of Franz (1948, 1949).

-13 - y -0*10 - ~.x = -w/2 0 x= w/2 Fig. 2.2: An axial view of the strip.

-14 - Of interest here, is the prediction of the radar cross section of the strip for various angles of incidence, resistivities, and electrical strip widths. 2.2 Integral Equation Formulation Formulation of a scattering problem in terms of an integral equation requires an integral representation (e.g., Stratton-Chu, 1939, or Franz) of the scattered field in terms of the surface current density on the scattering body. Expressions for the total fields in the space near the body can then be written. An integral equation for the surface current density is obtained if the observation point is allowed to approach the surface and the appropriate boundary conditions are applied (Maue, 1949 and Poggio and Miller, 1973). A numerical solution of the integral equations of scattering problems is normally carried out using the method of moments (Harrington, 1968). Discretization of the equation and the unknown by this method results in a system of linear, algebraic equations which may be solved simultaneously using a digital computer. 2.3 Differential Equation Formulation Maxwell's equations describing the electromagnetic field at any point in the region of free space surrounding the strip are V A E iwnojH (2.8) and v, H = -icwc E 0 - (2.9)

-15 - The strip is assumed to be longitudinally invariant along z such that all field components are independent of z. Expanding (2.8) and (2.9) in component form yields DE DE A E E x.... ay - x z =y iP (H x + Hyy + H z) (2.10) ay DX y x y 0 x z and AH ^ ) H 3H A z x+ ( - X) z = -iW (E x + E y + E z) (2.11) ay -— y x -: x y z Decomposing (2.10) and (2.11) into scalar component equations yields the following two independent systems of equations: DE ZE- i H ay = i,UoH x DE Z-i = i y > (2.12) - X - 0-iwE aH aH ax ay o z and Hz -iwc E aH x ax = Ey (2.13) DE DE Ey iW H ax Dy o z

-16 - Examination of the preceding equations reveals that (2.12) couple field components Ez, Hx, and H describing the field configuration ' x' y of an E-wave, while (2.13) couple field components Hz, Ex and Ey describing the field configuration of an H-wave. Since the incident plane wave has its electric field vector directed in the z-direction, only E-wave fields will be present. Eliminating H and Hy from (2.12) yields the following two-dimensional wave equation for Ez in which k2 = o2: ' 0 0 0 a2E 3'E Z + z + k2E = 0. (2.14) ax2 ay2 The incident plane wave field E1 satisfies the same wave equation, and thus it follows that the scattered field Esc = E - E satisfies D2ESC 2ESC z + - +k2Esc = 0. (2.15) ax2 ay2 z Since the scattered field is an outwardly propagating wave, it satisfies the radiation condition (Bowman et al., 1969, p. 5): 1im likoEP () =, (2.16) where p is a position vector in the x-y plane. At the strip the boundary conditions may be re-expressed as [Esc]+ = 0 (2.17) =.

-17 - and A R(x)Z0n, [H SCI+ = -n, (n E) (2.18) Since there is no discontinuity in the incident field and n = y, (2.18) becomes -R(x)Z[HSC ]+ = E = E +E Z Z Z (2.19) Further, because (2.12) applies to the scattered fields, aESC HSC i z x WI ay and (2.20) za [ay = E + E Z Z (2.21) In the case of a perfectly conducting strip, R(x) = 0, implying that ESC E Z i = -E z (2.22) at the strip. Some symmetry properties of the scattered field evident from an integral representation (e.g., 2.38) should be noted at this point: Esc(xy) z Hsc(xy) x = ES (x,-y) = -Hx(x,-y), 5 (2.23) (2.24)

-18 - This means that (2.19) can be written as -2R(x)ZH + = E + ES. (2.25) Using (2.20), this becomes ESC 2iR(x) Z = E + ESC (2.26) z.+ z on the strip. Further, from symmetry considerations, HSC equals zero off the strip in the x-z plane, implying that aEsc = 0 (2.27) ay for Ixl > w/2, y = 0. These symmetry properties allow the solution of the problem by solving the wave equation only in the half-space y > 0. The boundary conditions at the edge of the strip are not as easily deduced as they were elsewhere on the strip. At the edge there is some ambiguity concerning the definition of normal and tangent vectors. Meixner (1949) proposed an additional boundary condition which, when imposed, makes the solution to the boundary value problem unique. His "edge condition" states that the electromagnetic energy density must be integrable over any finite domain even if the domain contains singularities of the electromagnetic field. Bouwkamp (1946) showed that in the case of a perfectly conducting infinitesimally thin

-19 - plane surface with an edge, near the edge the singular components of the electric and magnetic field vectors are of the order of r-1/2, where r is the distance from the edge, while components of the fields parallel to the edge are always finite. It should be noted here that the singularities which do occur are in the scattered field components and not in the incident field. Also, the induced surface current density parallel to the edge of a perfect conductor is infinite at the edge. Extending the work of Meixner to deal with imperfectly conducting edges, Fawzi and Burke (1974) concluded that there are no singularities in any of the field components at the edge of a nonmagnetic, finite conductivity wedge. Senior (1979) found that the current density along the edge of a resistive sheet is always finite. A numerical solution for the scattered field E c using the differential equation formulation is carried out by the finite element method as described in the following chapters. 2.4 Far Scattered Field and the Radar Cross-Section In order to obtain the far scattered field, the induced surface current density K on the strip is calculated. Using (2.7) and (2.20) Sc - - rE. ( +28) Kz -[Hc SC - k 1(2.28) x - k z a

-20 - For a resistive sheet, (2.19) yields K 1 (Ei z R(x)Z ( 0 +E ) z (2.29) where R(x) f 0. Using Franz's representation, the scattered field may be expressed as ES(r) = v ^ v ^, (r), (2.30) where - is a Hertz vector of the electric type given by -n(r) i WE 0 S 5 K(r')G (r,r')ds' - 0 - (2.31) in which G (r,r') = e - 4TI r - r' I is the free space Green's function. This field satisfies the radiation condition at infinity. Since VAV II, = v(v iiR) - V211l and V2'1 + k21n = 0 - 0 - (2.32) in source-free regions,

-21 - ESC = v(v * ) + k2II (2.33) o0 -at any point off the sheet S. Since the problem being considered in the present work is two-dimensional (no variation in the z-direction) and since K = K z, then l = n z and v * - = O. Thus: E z s = Ez = k z, (2.34) where Kz = H )(kol - p')dp'. (2.35) o P Here r is the strip and (i/4)HOl)(kolp - p'l) is the free space Green's function in two dimensions. In the far field region kolp - p_' >> 1 and p >> p'. Thus (1) 2 ik p-i(1/4) -ik p*p' Ho (o - P') rrk e e. (2.36) ^ O In terms of source coordinates, p' = x'x and pIp' = x' cos p. -=Figure 2.3 shows the relationship between these various quantities. Therefore, in the far field region Ez (P) = JT p e ip-i() (2.37) 2 o-i (40

-22 - y / / / / P / I 'I I x=-w/2 0 x = w/2 Fig. 2.3: A diagram for the far field computation.

-23 - where k w/2 -ik x' cos k P(,o) - 4 ZoK (x') e dx'. (2.38) x =-w/2 The angle of incidence 0o enters the integral by way of the current density. For a two-dimensional problem the radar cross section is defined (Bowman, 1969) as ESC 2 Co(,0) = lim 2rp E1 (2.39) p + o0 which becomes c(- o) = P(, (o)2 (2.40) ko in the present case.

CHAPTER III. SOLUTION BY THE FINITE ELEMENT METHOD 3.1 Summary of the Problem The scattering problem formulated is essentially a scalar problem in two dimensions. The unknown field ESC is denoted by the scalar u here and in the following chapters. The scalar u satisfies scalar u here and in the following chapters. The scalar u satisfies au2 -- + ax2 u+ k2 = 0 ayZ ~ (3.1) On the strip located at y = 0 and Ixl < w/2, 2iR(x) ko - u+ u1 (3.2) Here ui = exp[-ik (x cos o + addition on y = O, Ixl > w/2, y sin 0)], the incident field. In ua= 0 Ty (3.3) The field at the edge of the strip must satisfy Meixner's edge condition if the strip is a perfect conductor. Finally, the radiation condi ti on au(p ) lim ua — )- ik u(p ) p ->+ o = (3.4) must be satisfied. -24 -

-25 - 3.2 Derivation of an Integral Statement of the Problem The solution for the scalar field u is to be obtained in the region Q bounded by the contour r as shown in Fig. 3.1. The region lies in the upper half-space y > 0. The contour r is composed of three sections: 1. rp is the line y = 0, -w/2 < x < w/2. 2. r0 are the lines y = 0, |x| > w/2. 3. rB is a semicircle with a radius p which is sufficiently large to place rB in the far field region (see Section 5.2). Green's first identity (Stratton, 1941, p. 165) indicates that for scalar functions u and v vv2u d = - vtv * vtu do + v - dr, (3.5) ^ I r where n is an outward normal to r and vt = x(2/Dx) + y(D/3y). The function u is the scattered field which is continuous and has continuous first and second derivatives in Q expect at the edges of the strip (x = ~w/2) if the strip is a perfect conductor. In this case singularities occur at the edges in the first derivatives which are integrable in accordance with Meixner's edge condition. The sum of the second derivative terms, v2u, has no singularity anywhere since u has no singularity and u satisfies (3.1). The function v is an arbitrary function at this stage except that it must be sufficiently continuous so that any singularities in the first derivatives are integrable.

-26 - y / I I I Q -p I I T - - -EMMON4 x=-w/2 7 F p x=W/2 F 0 0 Fig. 3.1: The solution region.

-27 - Use of (3.1) through (3.3) results in 0 2 v(u + f (tv v vtu - k2vu)dQ - i 2 S U dr Frp P - v -u dr = 0. (3.6) rB In the far field region u is approximately given by (2.37) which satisfies (3.4). Then on r, L1 - P=P -au =au = - 2-1 + ik uJ (3.7) Using this, (3.6) becomes k (vt u - kvu)dQ - i (u ) dr Q rp + P - ik f [vu]= dr = 0 rB 1 (3.8) This is known as a "weak form" integral statement (Zienkiewicz, 1977) since the order of the highest derivative is one lower than the order of the highest derivative in the original differential equation (3.1). In the case of a perfectly conducting strip, R(x) = 0, and (3.8) cannot be used since R(x) is in the denominator. Returning to (3.5), and using (3.1) and (3.7) yields

-28 - (vtv * v u - k2vu)da - f udr t t d n 2 r P + i] ik [vu dr = 0. (3.9) For the perfect conductor a prior knowledge of au/an on the strip (rp) does not exist, except that the derivative is proportional to the induced current on the strip, but the boundary condition, u + u1 = 0, on rp is yet to be imposed. Equation (3.9) is the integral statement for the perfectly conducting case. The derivations of (3.8) and (3.9) are essentially equivalent to the weighted residual method described in the finite element literature (Zienkiewicz, 1977, Chapter 3). 3.3 Interpolation Functions and Expansion of the Unknown In order to solve for u, it is approximated by I u: u = Nn(x,y)a (3.10) n=i and v is restricted to be v = w (m = 1,...,I) one of a set of weight functions yet to be prescribed. The domain Q is divided into a number of small domains Qe called elements with nodes n on the element boundaries. As indicated in the summation above, there is a total of I nodes. For each node n there is an interpolation or shape function Nn which is unity at that node and zero at all other nodes. The value of the unknown u at node n is then approximately

-29 - an, a complex constant. Thus, the determination of the scattered field amounts to determining the value of each of the an constants. Selection of interpolation functions and weight functions must be made in accordance with the necessity that the integrals of the weak form statement (3.8) or (3.9) be finite. In the weak form statement, the appearance of first derivatives requires that u and wm must be continuous throughout Q and on r, but allows discontinuities in their first derivatives. Further discussion on the selection of interpolation functions appears in Section 3.5. Another name is often attached to the finite element weighted residual procedure depending on the choice of the weighting functions. If Wm N (3.11) m the method of Galerkin is being used. In energy conserving problems, such as in a perfect dielectric medium, this method leads to a symmetric matrix of coefficients for the set of algebraic equations from which the an's are determined. Also, for this case, the method carries a guarantee of convergence to the proper numerical result, which other choices of weighting functions do not have (Zienkiewicz, 1977, Chapter 3). 3.4 Conversion to a Set of Algebraic Equations Using the approximation u for u and Galerkin's method in selecting the weighting functions wm one obtains for the weak form statement (3.8):

-30 - f Nm S^ I [z NnanA n= n a Ln=i + 0- I ay ' n=i nan 1Na dQ Nnan dr n=l I - k0 f Nm Nnan + 1 k r]f Ln=l ik dQ - -0 2 Nnan] dr f Nm R( x) ik 0 -2 rp N ui mR(x) dr (3.12) for m = 1,...,I. Interchanging the order of integration and summation, {JN9 Q Iff a UL" a 3N n Dx + k2N N dQ an ay -9y 0 a n n I i k 0 2 n= rp i'P Nm Nn R(Ix) I an + -n-= 2p ik) L ' 'B NmN dr} m n j a n ik 2J 1p N mu R(x) dr

-31 - for m = b.I This is written as ~ n= K a mn n f u nu = 5 (3.1 3) where KD nun DNn + 31m 3Nn 3x + y ~y i k pNmNn __ dr r p i k Im: -- 0 - k2NN] do + \21 - i-ko NmNn dr (3.14) and fp Nm m- d r (3.1 5)

-32 - Equations (3.13) represent a set of linear algebraic equations to be solved for the unknown a 's, the values of the scattered field n at the nodes. Calculation of the coefficients Kn and f. in the finite element method is done by an element-by-element procedure. The integrals over the entire region 2 are replaced by sums of integrals over individual elements making up the domain. Thus J K = Kn (3.16) and e=1 and J f = \ fe, fm j fm (3.17) e=1 where J is the number of elements, 9N DN aN aN n = f m n + n k2N N dQ in = D lx ax ay ay o m n e Q rp B and ik f N u fe = 2 J R0x dr (3.19) mn fm 2 Rx) dr. re P For a perfectly conducting strip, the set of algebraic equations should be modified by removing equations corresponding to nodes where

-33 - the field strength an is known. The values of an for nodes on the strip is known due to the boundary condition u + u = 0 on rp, Thus on rp an = -u1 (xn ) = -u n' n (3.20) In this case (3.10) can be written as u I I' - Nn(X Y)u + n=I+1n n=I'+l n=l Nn(.x,y)an, (3.21) where I - I' is the number of nodes on the strip. Proceeding as with the resistive strip using u for u and Galerkin's method in selecting wm, the weak form (3.9) becomes aNx ax F __ I n=I'+1 I ' N u1 + nn n=l Nnan) DN + m Dy I a n=I'+ n=I ' +1 I' N ui + < n nn n=l N a d2 - k2 Nm n n -0 Q I E N ui n n n=I'+1 I' n=+ n~l Nnan) dQ - rP N ( + n=I'+1 I ' N u1 + nn nJ n=l Na) dr n ^-ikj f (2p o 1 B p B m ( I I' n N ui + n n n n=I'+1 n=l Nnan d = 0 (3.22)

-34 - for mn = 1,...,I'. After interchanging the order of integration and summation, and transferring the terms corresponding to the nodes on the strip to the right-hand side, the resulting set of linear algebraic equations is I ' Kinnan g m 1,...I' (3.23) n=l where F "N aN 3N ONn Kn = m m - k2NmN do Kmn x yx my nyon B and I 9 = * (3.25) n=I'+l 3.5 The Element Mesh The selection of a particular type of element and its interpolation functions for a particular problem depends on, in addition to the nature of the integral statement, the geometry of the problem, the nature of the medium, the desired accuracy of the solution, and the computer facilities available for the computation. In the present study, the geometry of the problem and the desire for computerized mesh generation lead to the selection of the elements.

-35 - Most of the elements are general quadrilaterals having a node at each corner. In order that the node locations be generated automatically in the computer program, the nodes are placed at the intersections of the ellipses and hyperbolas of an elliptic coordinate system (Stratton, 1941, p. 52). The sides of the elements are straight lines between the nodes (Fig. 3.2). Triangular elements are used at the edges of the strip. This is done because of the need to properly model the behavior of the fields at the edge of a perfectly conducting strip. Since some components of the fields become infinite at the edge, an element is needed with shape functions whose derivatives have appropriate singularities. Such elements have been developed by W. S. Blackburn (1973) and a number of others (Zienkiewicz, 1977, Chapter 23). It can be stated as a general rule that use of higher-order interpolation functions requires fewer elements to achieve a given degree of accuracy of solution. Thus savings in computation time and computer memory requirements can be obtained. On the other hand, the higher-order functions lead to more complicated computations for the element matrix coefficients (Zienkiewicz, 1977, Chapter 7). In this work the quadrilaterals used are of higher order than the more commonly used linear triangular elements and offer the potential of greater accuracy. Linear interpolation is used in the quadrilaterals rather than higher order interpolation functions in the interest of keeping the element matrix computations simple.

-36 - Y X 0 Fig. 3.2: The element mesh.

-37 - 3.6 Quadrilateral Finite Elements As can be seen in Fig. 3.2, the quadrilateral elements are of various sizes and shapes. Rather than write interpolation functions for each of these elements, a transformation of the elements to a set of normalized coordinates allows the writing of the interpolation functions for all the quadrilaterals in a single concise form. A general quadrilateral in the x-y plane is shown in Fig. 3.3. The x- and y-coordinates of the nodes are known values. A transformation of coordinates is used to map this quadrilateral into a square element in the a-3 coordinate system as shown in Fig. 3.4. That is, the x-coordinate of any point in the element is a function of a and S. The same is true for the y-coordinate. Thus, x = x(a,3) (3.26) and y = y(a,3). (3.27) Before the exact expressions for this coordinate transformation are stated, a discussion of the interpolation functions Nn for the quadrilateral elements will be given. The simplest interpolation function is a linear function. It provides for continuity of u at element boundaries but will allow discontinuities in the first derivatives of u. The function Nn must satisfy

-38 - y (x3Y3) (x4,y4) 2 (x2,Y2),yl ) x Fig. 3.3: A quadrilateral finite element.

-39 - B (-1,1 ) 4 (1,1) 3 C (-,-1 ) (1,- ) 1 2 Fig. 3.4: A square element in normalized coordinates.

-40 - N = 0 at node n (3.28) at all other nodes A function meeting this requirement and having a linear variation with position n = (1 + an)( + Sn ) (3.29) is selected where, n n = value of a at node n, = value of B at node n. Thus in the a-3 coordinates the value of u at some point (a,c) in the square element is 4 u(a,) n = n=l Nnan (3.30) Linear interpolation is also used for the coordinate transformation so that (3.26) and (3.27) become x (.a,f) = 4 Nn (a, O)xn n=l 4 and (3.31) (3.32) y(a,3) = n=

-41 - where xn and Yn are the x- and y-coordinates of node n. Because the same interpolation functions are used for the function u and for the coordinates x and y, these elements are said to be "isoparametric". In order to carry out the integrations for the element matrix coefficients Kn (.3.18) or (3.22) and the incident field term fe (3.19), the integrals are transformed to the a-3 coordinates, simplifying the specification of the limits of the integrals. The details of these calculations are presented in Chapter IV. 3.7 Triangular Elements at the Strip Edge At both edges of the strip are regions which may be conveniently divided into two triangular elements rather than one quadrilateral element (Fig. 3.5). Integrations over these elements are most easily handled if the elements are transformed in a manner similar to that used for the quadrilateral elements. A general triangular element is shown in Fig. 3.6. This element is transformed into an a-f coordinate system as shown in Fig. 3.7 using linear interpolation functions: x(a,:) = (1 - a - )x + ax + x 1 2 3 and. (3.33) y(A,) = (1 - - a )y + ay + y 1 2 3

-42 - y x 0 +w/2 Fig. 3.5: Triangular elements in mesh at strip edge.

-43 - y (x3, 3) 3 2 (x2,Y2) 1 (x1 y1) X - Fig. 3.6: A general triangular element.

-44 - 3 h (0,1) 1 \(1,0) (0o,) - 2 Fig. 3.7: A triangular element in normalized coordinates.

-45 - The interpolation functions Nn are N = 1 - a -, 1 N = a, (3.34) 2 N =. 3 The field is interpolated with the same shape functions 3 ua,) => Nan, (3.35) n=l and the element is isoparametric. This approximation for (3.35) is valid for a resistive strip, but for a perfectly conducting strip a special set of interpolation functions N = 1 - a + N = a/VJ~ + (3.36) 2 N = e//a + \ 3 are used to model the behavior of the fields near the edge. Following the method of Blackburn (1973) this set of functions not only models the behavior of u near the edge, but also has a linear variation along the side of the triangle opposite the node at the edge of the strip to assure continuity of u at the element's boundary with the adjacent quadrilateral element. Since the interpolation functions (3.34) and (3.36) are not the same, this special element is not isoparametric.

CHAPTER IV. SOLUTION PROCEDURE 4.1 The Computer Program The computer program for the numerical solution generates the mesh of nodes and elements by calculating x and y coordinates of every node. The elements and nodes are then automatically numbered. Each node is given a unique global node number and a local node number associated with each element of which it is a corner. Next, the element matrix coefficients Ke and incident field terms fe are mn m calculated for one element and inserted into the proper location of the global matrix for the system of equations (3.13) or (3.23). This process is repeated for all the elements. The resulting global matrix has all its nonzero coefficients in a relatively narrow band centered on the principal diagonal. Computer memory requirements are kept to a minimum by storing only the band of nonzero coefficients. The system of equations is solved by a procedure beginning with a triangular decomposition of the matrix followed by a variation of Gauss elimination. Solution of the system of algebraic equations yields the value of the scattered field at every node. The induced current on the strip is then calculated, and from this the radar cross section is obtained. Details of the calculation of element matrix coefficients, the incident field terms, induced currents, and the -46 -

-47 - radar cross section are given in the remaining sections of this chapter. 4.2 Element Matrix Coefficients for Quadrilateral Elements The integrals for the element matrix coefficients Ke (3.18) or mn (3.24) are transformed to the a-s coordinates using (3.31) and (3.32)., The first portion of Ke mn m N +N aN kN m n + m n - k2NmNn dx dy 0, L 1 =-l N N aN aN ay m n k N d ds Dx 9y 3y omn d dcidI (4.1 ) where IJ| is the determinant of the Jacobian matrix J. Dy J = I~ (4.2) and I I = det J = a y ax a y Da 3' 3^ 3 - B a, (4.3) Using (3.31) and (3.32), the partial derivatives of J are

-48 - ax 1 aj = 1 (D + 3D ) (D + aD ) 3 2 (4.4) (D + BD ) a- 4 5 a = (D + D ) 4 6 5 where D = - + x +x -x 1 1 2 3 4 D = x -x +x -x 2 1 2 3 4 D = -x -x +x +x 3 1 2 3 4 (4.5) D = -Y + y +y -y 4 1 2 3 4 D = y -y +y -y 5 1 2 3 4 D = -y -y +y +y 6 -1 2 +3 + 4 The terms in (4.5) are the coordinates of the nodes. The partial derivatives in the integral of (4.1) are ann_ N aN N n 1 Naan n _a\ (4.6) DX Till Da D6 6 a /4.") and 9Nn 1 ( 3n ax N n ax) a 1 N ny Tn D3 a a _ (.y 8 Since N is given by (3.29), then 3N Da 1n(l a + '-) (4.8) = 4~n(~ aS

-49 - and aN as 4 = ( + aa) * (4.9) Although the transformation of the integral (4.1) to the a-8 variables simplifies the limits of integration, the integral is not easily evaluated due to the complexity of the integrand. Numerical integration (Gaussian quadrature) is used in the computer program. The next portion of Kn given in (3.18) is the integral across the resistive strip. This integral has a nonzero result for only those elements adjacent to the strip. For a perfectly conducting strip, Kn has no such integral on the strip as can be seen in (3.24). Figure 4.1 shows a typical element adjacent to the strip. From Fig. 4.1 it is apparent that ik N NN ik x4 [NmNn ]= dx = m - n o dx (4.10) e x rp 1 Transforming the integral to the a-s variables and using =x d x dx = da + (4.11) yields ik o NmN iko m n d: 2 RLx) 2YT R - =_e _=-1 a=-I 'p

-50 - (x3,y3) 3 r (x2 Y2) X (x,y1) x = w/2 Fig. 4.1: The quadrilateral element adjacent to the strip.

-51 - since the line between nodes 1 and 4 corresponds to the line a=-l in the a-B coordinate system. To allow for a number of different resistivity profiles, R at node n is assumed to be of the form Rn = A + BX + CX,(4.13) where A, B and C are constants which are part of the input data for the computer program. If A, B and C are zero, the computer program interprets the data to indicate a perfectly conducting strip. The expression for R(a,B) along the line a = -1 between nodes 1 and 4 is a linear interpolation in terms of R and R 1 4 [R(a,B)] = N R + N R X=-1 1 1 4 4 or [R(a,p)] - (R + R ) - B(R - R ). (4.14) 1 4 1 4 Using (4.4) and (4.5) yields - (x - x ) (4.15) and using (3.29) yields [NmNn 1 6 (1 - am)( 1 - n )(l1 + * )(4.16) a =- 1

-52 - Thus (4.12) becomes iko NN iko(1 - m)(l - n)( - 4 ) m n - J Rx dx = - - --- 32x e 1 p (1 + 6Bm)(1 + Bn )d6 (R + R ) - e(R - R ) (4.17) - 1 4 1 4 using (4.14), (4.15) and (4.16). This integral is also evaluated by numerical integration. The final portion of Kmn given in (3.18) or (3.24) is the integral along the contour rB in the far field region. This integral has a nonzero result only for the elements adjacent to this contour. Figure 4.2 shows a typical element adjacent to the strip. The integral is transformed to the a-6 coordinates. If r is the distance along rB from node 2, the integral is distance along d (21 - ik) f NN dr = ( - ik) Nf [N n d. (4.18) 2p O m n 2p o m n ' 1e 1 3=-c =1 B The expression for r is r (x - )2+ (y y )2 (1 + B).(4.19) 2 3 2 3 3 However, x = p cos (, y = p sin q, x = p cos 2 1 1 2 1 1 3 1 2 and y = p sin. 3 1 2

-53 - I / f / / / g I /1/ /0 // 01 -[ 0 Fig. 4.2: A quadrilateral element adjacent to the contour r.

-54 - Thus (4.19) becomes P r = -L /I cos(G - ) /-T 2 1 (1 + f) (4.20) Then dr d i=1 P = - /1 - cos ( - p ) dB 2 1 (4.21) Using (3.29) yields [NmNn] 0:= 1 1 6 (1 + a)(1 + )(1 + ( )( + n) (4.22) The integral in (4.18) becomes f e B N N dr m n p 1 - cos( 2 - )7 1 6 2 1 (1 + am)(l + an) (1 + MB )(1 + 3B )dg (+ m n (4.23) when (4.21) and (4.22) are used with (4.18). The integral in (4.23) is evaluated exactly so that 1 ik) 13 N N dr m n 2p 1 i k) p /1 - COS (2 - ) 8 /2 * (1 + am)(1 + n)m + mn) (4.24)

-55 - In summary, the element matrix coefficient K8n (3.18) for a quadrilateral element is the sum of three terms given by (4.1), (4.17) and (4.24) for a resistive strip. For a perfectly conducting strip Ke (3.24) is the sum of (4.1) and (4.24). mn 4.3 Element Matrix Coefficients for Triangular Elements The integrals for the element matrix coefficients Ke for mn the triangular elements are handled in much the same manner as those for the quadrilateral elements. The first portion of Ke is transmn formed to the a-6 coordinates N 3N 3N aN m n m n 2NN dxdy x x y y - km d e 1 l-ao:J J + m n k2NmN JJj da d (4.25) a=0 8:0 L where integration is over the triangle shown in Fig. 3.7. Using (3.33) yields ax X y = x = x -y 3a 2 1 yc 2 1 (4.26) x = x - x = y -y ~ 3 1 aP 3 1 so that

-56 - II = (x -x )(y - y )- (x 2 1 3 1 3 - x )(y - y ). (4.27) 1 2 1 The interpolation functions (3.34) are written in a general form Nn = an n - an - n, (4.28) where a = 1 1 a = -1 2 a = 0 3 = 1 1 B = 0 2 = -1 3 The partial derivatives in (4.25) Using (4.26) and (4.28) yields are given by (4.6) and (4.7). aN ax n I 3 n 2 1 aN n 1 a X t1 (x - x ) + ~ (x - x )] 3y |JT n n 1 2 (4.29) (4.30) and Then using (4.28) through (4.30), the integral (4.25) becomes

-57 - f [aN aN aN aN m n m n - k2 N dx d ax ax ay y om x dy Qe L[N = f [( - X )2 + (y - y )2] - (a + m a ) = 1( 3 1 3 1 m mm) n m n 2 3 2 1 3 1 2 1 3 1 mn 2 1 + (y - y )2] - 2:1 J k2 IJI a l ( m + a a ) 0 k n n - 6 Wenm+ma n'nam - 1 (a m n + an1 ) + I (a + anm) + (aa + n 6(am in n n (2 m n n ' m m n (4.31) The next portion of Kmn given in (3.18) is the integral on the resistive strip. This integral is nonzero only for the triangles which are adjacent to the strip. Figure 4.3 shows a triangular element adjacent to the strip at x = w/2. From the figure it is apparent that ik r 0 -I re p X3 NN i xk R(x) = -2 j x=x1 [NmNn ]yo dx R(x) (4.32) Transforming the integral to the a-B coordinates yields ik N N ik 2 \ R (x dx = 2-J p.,e B=0 'p N N m n R(a, _) L- -o=O (-X) d, (4.33)

-58 - 2 (x,y ) 2 2 (x3, y3) X w = w/2 Fig. 4.3: A triangular element adjacent to the strip at x = w/2.

-59 - since the line between nodes 1 and 3 corresponds to the line a = 0 in the a-B coordinate system. The expression for R(a,B) along the line a = 0 is a linear interpolation from R and R: 1 3 [R(a,3)] a=O = R + 3(R - R ) 1 3 1 (4.34) Using (4.26), (4.28) and (4.34), the integral (4.33) becomes iko NmNn ik0(x -X3) 2J R( — dx 0 - re =o PP (m - B)Bm( - n )d, R + B(R - R ) 1 3 1 (4.35) A change of variable in the integral, B = (a + 1)/2, yields ik N NN d er p ik (x - x ) 8 1 (2am - 1 - )Bm(2n - 1 - n)n — (R + R ) + (R - R )- d 1 3 3 1 (4.36) a form which is convenient for numerical integration. Figure 4.4 shows the triangular element adjacent to the strip at x = -w/2. The equivalent of (4.36) for this triangle is

-60 - y (x3 Y3) 3 ) x x = -w/2 2 Fig. 4.4: A triangular element adjacent to the strip at x = -w/2.

-61 - ik N N - rr m ndx e rp P 1 iko(x2 - l) f (2m - 1 - o)am(2n - 1 - o)an = ( R + R ) + (R - R d- (4.37) C-i 1 2 2 1 which is evaluated in the computer program by numerical integration. In the case of a perfectly conducting strip, the element matrix coefficient (3.24) for the triangular elements becomes aN aN aN aN k2N N dxdy. (4.38) n D ax ax ay ay dx (4 Qe - This integral is transformed to a-8 coordinates using (3.33): 1 i-ao e aFN aN 3N aN n = J J m n + m n k2N NIJ da de (4.39) mn DJax ax ay ay o m n a=O 3=0 where IJI is given by (4.27). The interpolation function N is given by (3.36) or can be written in a general form a a + Bn n = a -, (4.40) n nn Va + where an and 3n are given by

-62 - c = 1 = 1 1 1 a = -1 = 0 2 2 a = 0 = -1 3 3 The partial derivatives in (4.39) are given by (4.6) and (4.7). From (4.40), aN a a + (-2a + + ) n _ n n n (4.4n) (4.41) 2(c + )3 and aNn (an - 2Bn)a )- Bn (4.42) 2(a + )3/2 Then using (4.26), the partial derivatives (4.6) and (4.7) become n 1 n n n ax - T-JT — l [ +n (-2c + B3) B (an - 2 Y (4.43) Ke (4.39) is evaluated by exact integration and using (4.40), (4.43) 2(( + 3B)2 2 and n 1 ((4n - 2Dn)a - Bn L 2 1 ar " TT; 2(,m. )3/2 --a an + (-2an n ) 2- + B-s/2n (x - x ). (4.44) e (4.39) is evaluated by exact integration and using (4.40), (4.43) and (4.44). Defining

-63 -F = (x - x )2 + (y - y )2 1 3 1 3 1 F = (x - x )(x - x ) + (y - y )(y - y ) (4.45) 2 2 1 3 1 2 1 3 1 F = (x - x )2 + (y y )2 3 2 1 2 1 the coefficients Kn for the triangular element are Ke 1 (F -2F + F ) K1 = 4-TJ ( F - F)1I 2 3J -3 = 41T 21 + 2 3 30 Ke = k^rlJ! +F +2F -J Ke = - 1 F + F + F 12 = 1 2 2 3 3 0 Ke -F + F 3F) 0 13 2 2F 2 3 30 Ke 3 ( F + F + 1 F k J 21 2 23] 30 eK7 4 1 \ 2 k!JI Ke F + - F + 4 0 (4.46) 22 - 3 231 9 Ke -- (5F + 7 F + 5 F3 Ke F +F 3 F k 31 2 1 2 3 30 33 1 32 33 9 /

-64 - In summary, the element matrix coefficient Kn (3.18) for a triangular element is the sum of two terms (4.31) and (4.36) or (4.37) for a resistive strip. For a perfectly conducting strip Ken (3.24) is given by (4.46). 4.4 Incident Field Terms The integral for the incident field terms for the element adjacent to the resistive strip. (Fig. 4.1), (3.19) becomes (3.19) is evaluated For a quadrilateral N u1 dx mR R(x), fe nl ik - 0o 2 e P i x N ui ik 4 m dx = ~ x 1 or fe ni 1 iko 2 B=-1 N mui (,) (x) R(a,B) J' =-1 r Lu(aB (4.47) upon transformation to the the a-B coordinates is a-B coordinates. The incident field in u (c,f ) = exp[-ik x(ca,B)cos 01, (4.48) where x(a,B) is given by (3.31). Then u1 ( )jI Wx = -1 = exp- iko(cos o) 1 [(x + x ) - 3(x - x )] 1 4 I 4 (4.49) Using (3.29), (4.14), (4.15) and (4.49), the incident field term (4.47) becomes

-65 - e iko( - am)(xl - x4) x. /l f = - 8 4- exp -o-k 2 1xi - 1 (1 + e m)exp ik ( 2 -(R + R ) - (R - B=-1 1 4 1 which is evaluated using numerical integration. ) cos 1i (4.50) For the triangular element adjacent to the resistive strip at x = w/2 (Fig. 4.3), fe m ik 0 2 X 3 X=X, 1 N ul dx m R(x) or e m ik o Nu' (aB)1 2 J LR(a, B) 3=0o ( =O ( x\ (4.51) In this case the incident field ul(a,B) for a = 0 is obtained using (3.33) ul (oa,6) a=O = exp -ik cos (, [(1 - )x + ex ] 1 3 (4.52) Using (4.26), (4.28), (4.34) and (4.52), the integral (4.51) becomes fe m ik (x - X ) 1 - exp[-ik x 2 0 1 cos o] 1 6=0 (m - )m exp[iko(xl - x3)(cos o)B] R + (R - R ) d1 3: (4.53)

-66 - or, if B = (o + 1)/2, ik (x- x x) + fp = -~- exp -ik - 3 cos m 4 2 (R ) + (R + R)+ (R - R do:=-1 1 3 3 1 (4.54) which is evaluated using numerical integration. For the triangular element adjacent to the strip at x = -w/2 (Fig. 4.4), a similar calculation yields e =ik (x2 - exp l [-1 x1 x 1X /X x \ f = exp -ik 2 ) cos 2)(COs m(2am - 1 - o)exp iko 2 ) 2co - L l — -— oj do, (4.55) G=-1 (R + R ) + a(R - R ) which is numerically evaluated. In summary, the incident field terms are given by (4.50) for a quadrilateral, (4.54) for the triangle at x = w/2, and (4.55) for the triangle at x = -w/2. 4.5 Induced Current The induced current in the resistive strip is obtained using (2.29) or K - R (ui + u) (4.56) z R7x)Z 0

-67 - and is evaluated at the center of the side of each element adjacent to the strip. If the element is a quadrilateral (Fig. 4.1) where R(x) = N R + N R 1 1 4 4 x - x4 N -L 1 X - X 1 4 (4.57) (4.58) (4.59) and x - x N = 1 4 X - x 1 4 and Rn is given by (4.13). The scalar u on the strip is given by u = N a + N a 1 1 4 4 (4.60) and the incident field on the strip is given by u = exp[-iko x cos ) I (4.61) Since the x-coordinate of is along the x-axis is the center of the side of the element which x + x X = 1 4 c 2 (4.62) the current (4.56) at that point is written as

-68 - X1 L G a + aX exp -ikoi + o 2 K + (4. z z0R(Ri + R4I by using (4.57), (4.60) and (4.61). Similarly, the induced current at the center of the triangular element at x = w/2 (Fig. 4.3) is.63) 1 Kz - Z 10 exp -iko 2 cos )Ioj Ri + R3 2 (4.64) and at the center of the triangular element is I 2 2 1z K z Zo )(R+R2 at x = -w/2 (Fig. 4.4) a + a ( 2 I +R + R V 1 (4.65) If the strip is a perfect conductor, the expressions for K having R(x) in the denominator cannot be used. Instead (2.28) is used. Due to symmetry (2.24) in the scattered field, (2.28) is rewritten using (2.20) to obtain 2i au z Z- o y (4.66)

-69 - For the quadrilateral elements (Fig. 4.1) the derivative in (4.66) is now evaluated. The scalar u is given by (3.10) or by 4 u z u = N a (4.67) n n n=i in a single element. Thus the derivative can be expressed by 4 3u 3Nn yN ay an a (4.68) n=i Here N is given by (3.29) and (aNn/ay) is given by (4.7). Since the element is adjacent to the x-axis, then y = 0 and y = 0. Noting that a = -1, where y = 0, (4.68) is evaluated to be l+u ) L( + x2 + 3 -x ) +(x -x + - )] y x - 2 3 4 1 2 3 4 -1 4 + (a - a )(1 - B) + (a - a )(1 + y + y ) + (y - y 2 1 3 4 3 2 3 2 (4.69) At the center of the side of the quadrilateral, B = 0 and the induced current is a1 a - - 4 (-x + x + - ) + (a -a ) + (a - a ) 2i x X4 1 2 3 3 4 2 1 3 z kZ -y +y 0 r)3 2 (4.70) For the triangular element (Fig. 4.3) adjacent to the strip at x = w/2, the derivative in (4.66) is now evaluated. The scalar u is given by

-70 - 3 u L u nan (4.71) n=1 where N is given by (4.40). Then n 3 ~ n V a (4.72) by + a Y y=o n ' n= where aN /ay is given by (4.44). Since nodes 1 and 3 are on the x-axis, y = 0 and y = O. Noting again that a = -1 where y = o, 1 3 (4.72) is evaluated to be au a (x2 - X3) + 2a2(xl - x3) + [(x2 - x1) + (x3 - xl)]a3 y 21/2 y (x -x) 2 1 3 (4.73) At the center of the side of the triangle, B = 1/2 and the induced current is 2i -(x2 - x3)al + 2(x1 - X )a2 + [(x2 - xl) + (x3 - x)]a K kzZ oo L (x - x )y 1 3 2 (4.74) For the triangular element (Fig. 4.4) adjacent to the strip at x = -w/2, the derivative in (4.66) is now evaluated. Since nodes 1 and 2 are on the x-axis, y and y equal zero. Equation (4.72) 1 2 becomes

-71 - au -a (x2 - 3) - a2[(x2 - x) + (x - xl ) + 2a3(X - 1) 2 + 2/2y (x -x ) 3 2 1 (4.75) At the center of the side of the triangle, a= 1/2 and i -( - x 3)a - (x2 - x )]a + 2(x - 1 )a3 Kz o- 2y(x -x) 3 2 1 (4.76) In summary, for a resistive strip the induced current Kz is given by (4.63) for a quadrilateral, (4.64) for the triangle at x = w/2, and (4.65) for the triangle at x = -w/2. For a perfectly conducting strip K is given by (4.70) for a quadrilateral, (4.74) for the triangle at x = w/2, and (4.76) for the triangle at x = -w/2. 4.6 Radar Cross Section The radar cross section given by (2.40) requires evaluation of (2.38), the expression for P(k,o0). Using (4.56) yields w/2 exp[-ik x'(cos * + cos q) P(') { - - ---- R(X') dx' x'=-w/2 w/2 u(x',0O)exp[-ik x' cos] p] + R(x') dx' (4.77) x'=-w/2

-72 - for a resistive strip. These integrals can be converted to a series of integrals over individual element edges along the strip. J P(q,M) = Pe(,mo) (4.78) e=1 where k = - exp[-ik0x'(cos o + cos,)] Pe (c^ ) = -4 {J R- ) +cos- dx' P u(x',0) exp[-ikox' cos ], + R(x') dx' (4.79) re P and J is the total number of elements. In (4.78) only those elements adjacent to rp make a nonzero contribution to P(k, o). For a quadrilateral element (Fig. 4.1), the integrals in (4.79) are transformed to the a-s coordinates using (3.31) and (3.32). k S - exp[-ik x(U, )(cos c + cos fo)] p o) = _ 0 + x e 4 R(c,) a =[ ru(a,B)exp[-ik x(ca,)cos (~] 1 I +R(a,) d. (4.80) L=-1 aJC=

-73 - Since u(a,f) is approximately given by (3.30), then u(ca,3) u c(t -i 2 (1 - l )a + I (1 + )a4 (4.81) Using (4.14), (4.15), (4.49) and (4.81) yields e - '4 FCik0 + A4C Pe(,o) = ~ (x - x)exp Li k( 2 (cos + cos o) 1 exp [iko ( 2 4) (cos ( + cos <o) f r -_- ________ d_ (R + R ) - 6(R - R ) 1 1- 1 4 1 4 + 2 xp -iko i_2 ) cos 1 (a1 + a4) - S(a1 - a4) *I ----- -- - -exp =-_ (R1 + R4) - B(R1 - R4) iko ( 2 4 )(cos 4)j dCJ (4.82) For the triangular element adjacent to the strip at x = w/2 (Fig. 4.3)

-74 - p e(~ 9,q 0) e0 * -1 - - T x- x3)expL[ ik (XI + X3)(COS q + COS ~X1 3X exp i k0 2 (Cos ~p + Cos d (JR + JR )- - ~( R ) 1 31 3 + (X2 XL F + expL-ik IX1 2~ Cos 1 (a1 + a 3) - f3(a1 - a3 ) (JR1 + JR 3) - (R1- R 3) exp [iko (.X - X) (cos fl 4) ( 4.83 ) and for the triangle adjacent to the strip at x = -w/2 (Fig. 4.4) pe '~~0) 1 ~=- 1. k0 = - T (x 2 - x 1) exp Fik L 1 2 (Cos + cost 0 eXP 2 - Xi) (Cos ~p + Cos ~) (JR1I+ R2) - ~(JR2- R1) + X - X) exp [ i k (1 2) Cos 1 =3- 1 (a1I+ a2) - ~3(a2 - a1) (R + JR ) - (JR - JR ) 1 2 2 1 exp [iko X - (2 ) (Cos f)l] d (4.84)

-75 - For the case of a perfectly conducting strip, evaluation of (2.38) requires the use of (4.66) rather than (4.56). Thus P(cp,p ) k 0 w/2, w/2 exp -ikx' cos ~ dx" xI w/2 oy + 0 x':-w/2 - (4.85) As with the resistive strip, this integral is converted to a series of integrals along individual elements. Using (4.78) again, Pe(,qo ) e 0 ko p 2i mU* 1 dx' i -u exp -ikox' cos dx k0 y Y + (4.86) For a quadrilateral element adjacent to the strip (Fig. 4.1), transformation of (4.86) to a-s coordinates yields k P (4',4o) - 40 o W0 1 "X J 2-i expk [)y B: —I (4.87) Since x(a,B) is given by (3.31), then exp [-ikox(a,,),cosl 0o —1 ep X + X4) (1 ) (s (4.88). exp ik 2 (cos l)). (4.88)

-76 - Using (4.15), (4.69) and (4.88) yields i(x,- x) x + x Pe( o) 1= - 4 exp -iko(1 )cos 1 - I exp Lik 2 (cos ~) 6=- (Y + y ) + (y - y ) 1 3 2 3 2 r a-a * - - (D + BD )+ (-a + a + a a L X- x 2 1 2 3 4 1 4 + B(a - a + a - a )Jd, (4.89) 1 2 3, where D and D are given by (4.5). 1 2 For the triangular element adjacent to the strip at x = w/2 (Fig. 4.3), transformation of (4.86) to a-B coordinates yields 1 (k ) = 2i u x c a=oO (4.90) Using (3.31), (4.26) and (4.73) yields =-i(x -x) exp [-,ik x cos e - 2 xp 1 r,-al(x2_ 3) + 2a2(xI - x3) + [(x2 - Xl) + (x3 - xl)]a3 ] 2y (x -x ) 2 3 )( exp- k (x - x3)(cos q:)J J1o dB (4.91) p1/2

-77 - or, if B = (a + 1)/2, i(xl x3) e x + x p ( = — T 3 exp ~iok ( ' 2 3) (cos ) a Ix X3) + 2a (X - ^X + [(X - x) + (x3 )l { -a 2 y(x - x1 2 1 3 r X1 X3 exp [iko 1 2 (cos ))o] J /v' + 1 0i -1 Using (4.74), P (),2 ) may be written as (4.92) k (xK - X3) [ (XI + X3 0 0 Pe 5o) = - 4) z 0 Kz 2 iCos [ x I - x ) ] 1 exp [ik0 ( 3) (cos ) cb+ = - -- ~ - -- da + ~ /a + 1 C =-1, (4.93) where Z K is the induced current at the center of the element. 0 Z For the triangular element adjacent to the strip at x = -w/2 (Fig. 4.4), the expressions for P(p,O ) are similar to (4.92) and (4.93), that is

-78 - i(x2 - xl) [ 1 + X P (4~4) = -..-4r- exp -ik (cos e 4o 2 ( cos f a(x2 - x3) - a2[(2 - x1) + (x3 - xl)] + 2a3(x2 - t 2 y (x - x ) 3 2 1 1 exp ik0 1 X (cos 2)] - du (4.94) c=-1 c + 1 or k x2 - x x + X2 \ P ): --- ZoKz- 2 -)exp[ -iko - cos r x x exp [iko (2 (cos ])a] C7=-l where Z K is the induced current at the center of the element. O z The radar cross section is frequently given as a dimensionless quantity c(qp,<o)/x and is given in decibels. Since k = 27/x, (2.40) may be written as o(,o ) 2 = 2 I |P( |)I2 (4.96) - ''o'

-79 - Expressed in decibels, 10 log X ~ = 10 lo + 20 log IP(,O )1, (4.97) where P(, 0o) is given by (4.78). In summary, for a resistive strip Pe(,q o) is given by (4.82) for a quadrilateral, (4.83) for the triangle adjacent to the strip at x = w/2 and (4.84) for the triangle at x = -w/2. For a perfectly conducting strip, P e(,o ) is given by (4.89) for a quadrilateral, (4.93) for the triangle at x = w/2 and (4.95) for the triangle at x = -w/2. Numerical integration is used in the computer program to evaluate the integrals in these expressions for Pe(p,q o).

CHAPTER V. RESULTS OF THE NUMERICAL SOLUTION 5.1 Error Criteria for Numerical Results The numerical results of the finite element method (FEM) were obtained using the DEC-10 computer at Western Michigan University. They are compared with results obtained from an integral equation formulation of the same problem. The integral equation was solved by the method of moments (MOM) at the Radiation Laboratory of the University of Michigan (Knott, Liepa and Senior, 1973). The moment method results for the radar cross section of the resistive strip are known to agree quite well with measurements. The error in magnitude of any quantity Q calculated by the finite element method is defined as follows: IQIFEM - IQIMOM percent error in IQI = I — QMOM x 100. (5.1) I^IMOM The phase error is given by error in arg Q = (arg Q)FEM - (arg Q)MOM (5.2) 5.2 Selection of Element Mesh Parameters The objective to be met in selection of the mesh parameters is to provide a sufficient density of elements so that the incident -80 -

-81 - field is adequately sampled at the strip, and so that magnitude and phase errors in the computed scattered field are not large, allowing a certain accuracy of solution to be obtained. Of course, the cost of the solution is also a consideration. The semicircle rB is to be in the far field region of the scattered field. Selection of the radius p of rB is made using the criterion 2d2/x0 for the near field radius, where d is the largest dimension of the scattering body (Kouyoumjian and Peters, 1965). There is no reason to extend the finite element region further outward because the accuracy of solution is not improved and the cost of solution is increased. Thus - 2W2 2w2 (5.3) 1 A0 Along a radial line extending outward from the strip, the scattered field is essentially a decaying sinusoid of wavelength X The results of the one-dimensional problem (scattering from an infinite sheet) in Appendix C indicate that the maximum error in the magnitude of the scattered field is less than one percent if the density of the elements is 16 or more elements per wavelength. This density of elements is generally viewed as adequate for sampling a sinusoid. Thus if My is the number of elements in the radial direction, the selection is My > 16 x (5.4) 0

-82 - Using (5.3) yields 32 w2 M > 32 (5.5) Y - x2 0 In order to adequately sample the incident field on the strip, the density of elements there should also be about 16 elements per wavelength. If M is the number of elements adjacent to the strip, the selection is Mx > 16. (5.6) o A lower density may be used for cases where the incidence is nearly broadsi de. The element mesh, discussed in Section 3.5, has a nonuniform distribution of elements. Thus the element densities implied by (5.4) and (5.6) must be considered as average densities. Near the strip the density of elements is greater than near the semicircle Fr. When p, M and M have been selected, computation of the x x y parameters actually used in the computer program is begun. The integer number Mx is called NUMELS in the program. The semicircle of radius p bounding the finite element region coincides approximately with an ellipse of eccentricity g-1. For -1 < 0.5, the ellipse is essentially a circle. Taking p as the average of the semi-major and semi-minor axes of the ellipse yields. (5.7) 1 2'

-83 - Using (5.3) yields -1 = _ — (5.8) 4 o which is the required eccentricity. This value is the computer program parameter SMEC: SMEC = (5.9) 4 M 0 In order that the largest ellipse inside the circle p = p be nearly a circle, SMEC must be 0.5 or less. A value of SMEC equal to 0.5 corresponds to w/x = 0.5 according to (5.9) and p = w according to (5.3). Consequently, for w/x < 0.5 the parameter SMEC must remain at 0.5 and p must equal the strip width. Another of the computer program parameters is an aspect ratio EL. It is the ratio of the dimension of a quadrilateral element in the direction parallel to the strip to the dimension in the direction perpendicular to the strip. Assuming values for My, NUMELS, and SMEC have been tentatively selected, an approximate value for EL can be obtained from M - 1 M - 1 EL = Y - Y (5 10) - (NUMEL S)A (510) where A = -0.83 log SMEC + 0.17. This expression for EL is due to 10 the properties of the mesh generation portion of the computer program as explained in Appendix A.

-84 - 5.3 Effects of Element Density In order to examine the effect of element density on the accuracy of the solution, numerical values of the magnitude and phase of P(q5,po) for backscattering (q = 0 ) were obtained for a strip of width w/xo = 0.5. The radius in wavelengths of the finite element region p /X is 0.5 as specified by the selection rule (5.3). The errors in the numerical values are defined by (5.1) and (5.2). The numerical results indicate that an increase in accuracy in both magnitude and phase of P(4, o) for all angles of incidence is most economically accomplished, in terms of computer memory, by making element densities (M x)x /w and (M )x /p equal and increasing x o y o 1 them both. It is possible to make one density substantially larger than the other to reduce the errors at certain angles of incidence. However, the effects on the errors at other angles of incidence are not always predictable if this is done. Figures 5.1 and 5.2 show the errors in P(p,p0o) versus element density for R = 1. For 50 degree and 90 degree angles of incidence, element densities of 16 elements per wavelength in both directions appear to be adequate to reduce the error in the magnitude of P(q, o) to be less than one percent, and to reduce the phase error to 1.0 degree or less. The errors for incidence at zero degrees are larger, apparently due to the presence of a null in the radiation pattern near this angle. The effect of the null can be seen in Figs. 5.3 and 5.4 which show the errors in P((,co ) versus q for backscattering. The errors at the null which occur at about ~ = 25~ can be reduced by

2 I 900 0.-. c H-e cr LUJ 0ij x x 00 Io 00 S - - / I- -. -. N". N / / / 6- 0 — - - = 00 I -'3 - I N1 -4 71 -, -, -, 10 15 201 25 ELEMENTS PER 4A.VELENGTH Error in Jp(4,p0)I vs. element density, (M~ )(x/w) = (Mi) (x0/p ), R = 1.0, W/A p /A= u.,backscattering (p p1 Fig. 52 1: = 0. 5,

5 ____ ___500 Kx _______, —. A A - - -M 900 4) = 0 - -- 'o -- LI) LUJ LUI LU LUJ 0 - -5k~ 1- 0 0 - A Il. 0 - I co CY) I / -10 1 - / / -15 I 10 15 20 25 ELEMENTS PER WAVELENGTH Fig. 5.2: Error in arg P (4),4) ) vs. element density, (M )(x 0/w) = (M )(x 0/p ), R = 1.*0, wh = 0.5, = 0.5, backscattering (4' = 4) ~). 0 v1 0 0

-87 -4 o-'o (Mx (A 0/w) = 16, (M~ )(/p1 17 ~~~~(M )x 1w) = 20, (N )(x p)= 20 (Mx)(x /w) = 24, (M ( /p 24 0 -2/ 04 CD L.N LO F —6 LUJ LU0l-8 -10 -12 I 0 30 60 90 Fi g. 5. 3: Error in Plp(,p0)I vs. cp for R =1.0, wlx =0. 5, 0 0 p /x0 =0.5, backscattering 0)= 0

-88 - 0 — (M )(X/w) = 16, (M)(/p ) = 17 x )(Mx)(x /w) = 20, (My)(x/p ) = 20 4 -. --- (M )(x /w) = 24, (M)( /p ) = 24 x o ' y o 1 2 -LI. / 2 - \ CL: r,,. -4 - l -6 0 30 60 90 Fig. 5.4: Error in arg P(~,o ) vs. < for R = 1.0 w/x = 0.5, p l/Xo 0.5, backscattering (p = o0).

-89 - increasing the element densities. The same general results for R = 0 are shown in Fig. 5.5 and 5.6 except that the errors are larger than for R = 1. The errors in P(p, o) are plotted versus (1/2)kow cos p for three different strip widths in Figs. 5.7 and 5.8. Here both the element densities are 16 elements per wavelength and R = 1. Data for strips wider than 0.71 wavelength could not be obtained due to large amount of computer memory required. These two figures show that for (1/2)kow cos ( < 0.3 a the magnitude error is less than one percent and the phase error is less than 1.5 degrees for all strip widths. The wider strips have nulls in their radiation patterns for (1/2)kow cos ( 1 0.5 a and consequently have larger errors in P(~,5o) at certain angles due to the rapid changes in phase of the scattered field in the vicinity of a null. For the widest strip considered (w/x = 0.71), the errors for (1/2)kow cos ~ > 0.6 Tr are in approximately the same range as those for (1/2)kow cos ) < 0.3 '. Increasing the element density above 16 elements per wavelength will reduce errors mainly in the vicinity of a null. 5.4 Effects of Strip Resistivity Table 5.1 shows the effect of changing the strip resistivity on the error in P(p,po) for backscattering. The table shows that as R decreases below 0.1 the error grows significantly. In the range 0 < R < 0.05 the numerical results are highly inaccurate, except for near-broadside incidence. Two effects lead to the deterioration of the numerical results for low R. One can be explained by considering the finite element

-90 - 4 2 0 -2 0 H-e LUJ LUJ -4 -6 -8 -10 -12 -14 O 0 300 600 900 Fig. 5.5:. Error in P(, Ivs. q for R = 0.0~, w/x0 0.5,3 rp /x z 0.5, backscattering (p = cp0).

-91 - 1 0 8 6 4 LUL LU LUJ (-, LU 2 0 -2 -4 O 0 300 600 900 Fig. 5.6;: Error in arg P(p,p ) vs. ~ for R = 0.0, w/x~ = 0.5, 0 pk z 0.5,, backscatteri~ng (~=~) pI/ 0

-92 - I 61. t-0 rr —l Q'1-1 0 C-< M: >-) 'n Q0 r I 10 0 LO) CD CD C\j LO) ro 0 0 c< ~-< 3 -x o 0 I I I?.11 LO C) V) t7 0 I~z- U CD ) (\J C) 0 CD -eU) 0 u C~J * 0 U) -e-e- C: Q) -4-) (-) s* --- J 9. J6X-e I 4 x LO l:J- cl -- J C D C\j C:I- LO co CD I I I I 0I 1( 0 Vfdl NI d1Od3 3OVINJDHdd C\j

-93 — -co Q CC) CC *u 0 x /CD 0 coJ

-94 - Table 5.1 Error in P(c,p 0o) for Backscattering (( = o0) xx w p (,) k~ = 18, (M ) 0 — 17, - 0.5 05 1 o o -0.673% ^/ ^-4.118~

-95 - mesh as a radial transmission line. The boundary condition imposed at rB acts as an imperfect termination for the transmission line. Multiple reflections occur at rB and at the strip rp so that standing waves appear in the numerical data. Lowering the value of R increases the magnitude of the wave reflected from the strip and thus increases the standing-wave ratio. The transmission line interpretation and the nature of the error caused by the standing waves are discussed in Appendix B. The second effect leading to the poor numerical results for low R is due to the coefficients Kmn (3.14), corresponding to nodes on the strip, becoming much larger as R decreases. This results in the system of equations (3.13) becoming increasingly ill-conditioned. The effect can be partially offset by controlling the size of the term in (3.18) which contains R. This term is proportional to [R(Mx)xo/w]-1, where (Mx)xo/w is the density of elements along the strip. The data of Table 5.1 indicates that if R(Mx)x /w < 1, there will be large errors in P(U, o). A selection rule of R(Mx)xo/w > 5 (5.11) is appropriate. Figures 5.9 and 5.10 show the errors in P(),O ) versus the angle of incidence for backscattering with R = 0.1. The curves for constant (Mx)(x /w)R show how the accuracy improves as this parameter increases.

-96 - 8 -& / X3.6 -4/ -8 1. 8 Lu Lu -1 2 -1 6 -20, 00300 600 90C Fig. 5.9: Error in IP(t~ 0) vs. for R = 0.1, w =X 0.53, p /X0 z 0.51, (M )(xU/p ) 17, backscattering ( 1 Y Q

-91 - 2 0 -2 LUJ LU -4 -6 -8 -10 -12 -14 Fi g. 5. 10: Error p AI0 1 JU30 600 900 in arg P(p,~ ) vs. ~ for R = 0.1, w = 0.5,1 -. 0.5, (M y)(x /p1 ) z17, backscattering (~ = ~ )

-98 - 5.5 Solution Cost Considerations The cost of the numerical solution is a function of the number of elements used in the discretization of the space surrounding the scattering body. The number of elements is J = MMy + 2. (5.12) The actual number of unknowns, however, is equal to the number of nodes at which the scattered field is not known a priori. The number of unknowns for a resistive strip is Mu = (M + 1)(M +1). (5.13) For a perfect conductor the number is less, Mu (Mx (M My, (5.14) since the scattered field at the nodes on the strip is specified by the boundary condition there. The computer memory requirements for the finite element solution may be classified either as being fixed or as being variable with the number of elements. In the variable category, the largest block of memory is for band storage of the global matrix for the system of equations (3.13). The variable memory requirements (in words) are:

-99 - Resistive Strip Global Matrix Memory (band storage) = 6M M2 + 6M2 + 20 M M + 14 M + 20 M xy y xy x y (5.15) Matrices x y x y Perfectly Conductin= 14 M + 10 MStrip 10 M Perfectly Conducting Strip (5.16) Memory = 6 M M2 (band storage) x y + 6M2 + 8M M + 8 M y xy y (5.17) (5.18) = 2 M2M + 18 M M + 5 M + 11 M xy xy x y If the element densities equal in view of the results of terms of Mx using (5.3) if w/xo Mxxo/w and M x /p are to be kept Section 5.3, M may be written in 0.5. > 0.5. M = 2M w y xx 0 (5.19) For w/x < 0.5, it is required that p = w instead of (5.3) as discussed in Section 5.2. In this case M =M y x (5.20)

-100 - The number of elements (5.12) becomes J = 2M2 w + 2 x x 02 = M2 + 2 5 x w/X > 0.5 0 - w/X < 0.5 0 (5.21) Then for equal element densities, My may be eliminated from (5.15) through (5.18). For example, (5.15) and (5.16) may be written in terms of the element density M x /w and the strip width in wavelengths w/x0: Resistive Strip (w/x > 0.5) 0 E i Global Matrix Memory (band storage) = 24(M >) x w ( 5 W X 0 = / 2 + M w } 3 X + 40 ) 4 o 2 (W 7) ) 2t + 40 (Mx ' w + 14 (Mx 0 (5.22) Memory, Other Matrices 28M 2 3 x 2 28Mx) ( ) + 20(Mx w3j ) 0 0 +10 (M x0 (5.23) Table 5.2 shows the number of elements, number of unknowns and variable memory for different element densities and strip widths. Table 5.3 compares the memory requirement for the band-storage

-101 - Table 5.2 Variable Memory Requirement for a Resistive Strip \ x (M ~ = M ~ ) x w y 1 w Number of Number of Total Variable x w X Elements, J Unknowns, M Memory (words)* 0.25 18 25 1,240 0.50 66 81 6,064 16 0.71 178 204 25,160 1.00 514 561 123,200 0.25 27 36 2,020 20 0.50 102 121 10,540 0.71 282 315 47,762 1.00 802 861 230,480 Includes band-storage global matrix.

-102 -Table 5.3 Global Matrix Memory for a Resistive Strip for M 0~ = 16 x w Global Matrix Memory (words) w F.E.M* F.E.!.* M.O.M. X Square Matrix Banded Matrix Square Matrix 0.25 1,250 950 32 0.50 13,122 5,022 128 0.71 83,232 22,440 242 1.00 629,442 115,566 512 * X X M ~ = M ~ x w yp 1

-103 -global matrix to a full square matrix and to the matrix for a method of moments integral equation solution. Table 5.4 gives the actual computation time for three selected examples on the DEC-10 computer.

-104 - Table 5.4 Computation Time for a Resistive Sheet (Mxo/W = MyX/p1 ~ 16) Y1 w/x 0.25 0.50 0.71 o Number of elements, J 18 66 178 Number of unknowns, Mu 25 81 204 Matrix fill time (sec) 0.84 3.61 10.04 Solution of equations (sec) 0.19 1.10 6.22 Total CPU time* (sec) 1.41 5.10 16.60.-,........ Backscattering for one angle of incidence. Backscattering for one angle of incidence.

CHAPTER VI. VOLUME SCATTERING EXAMPLE 6.1 Statement and Formulation of Example One of the advantages of the finite element method lies in its ability to solve problems in which inhomogeneous media are present or multiple regions of differing homogeneous media are present. The distinction between these two types of problems is not important for a finite element solution. Only a minor difference in the computer program is needed in assigning material properties to particular finite elements. An example is considered here in which the resistive strip of the preceding chapters is embedded in a dielectric cylinder, radius w, of uniform relative permittivity cr and uniform permeability po as shown in Fig. 6.1. The cylinder is, in turn, immersed in free space. The radar cross section of the combined dielectric cylinder and resistive strip is computed. As before, the incident plane wave (2.3) will excite only E-wave fields. In this case the total field satisfies a2E D2E z + + z + 2 s(x,y)E = 0. (6.1) 3x2 y2 ~ The incident field E1 which is the field present in the absence of the scattering body, satisfies a2E1 92E1 - z + z + m2 o E1 = 0. (6.2) ax2 ay2 ~ 0 z -105 -

-106 - Y BF X Fig. 6.1: The extended solution region.

-107 - SC i The scattered field E = E- E then must satisfy D2ESc a2ESC z + z + 2 (x,y)E -k x(x,) (6.3) ax2 ay2 where X -1 = R-1 (6.4) 0 is the electric susceptibility for the dielectric. The derivation of an integral statement of the problem parallels that of Section 3.2. In this case however, the scattered field does not, in general, have the symmetry properties (2.24) and (2.27) which exist for the strip alone. Thus the region o must extend into the lower half space y < 0 as shown in Fig. 6.1. The region is bounded by the contour r which is composed of two sections: 1. rp is the line y = 0, -w/2 < x < w/2. 2. rB is a circle centered at the origin with a radius p which is sufficiently large to place rB in the far field region of the dielectric cylinder with the embedded strip. Combining the Green's identity (3.5) and differential equation (6.3) yields [vv vt u k2(l + x (x,y))vu]dQ = v a dr + k2 xe(x,y)vu dQ. (6.5) r 52

-108 - i where u and u are the scattered and the contour integral, the portion on rp, the boundary condition is iR(x) ul+ The contour integral on rp is incident fields as before. Of rB is the same as before. On + u + U (6.6) Iv d dr rp w/2 x=-w/2 u dx W-i (6.7) Combining (3.7),(6.6) and (6.7) with (6.5) yields the integral statement upon which the finite element solution is based. The result is [vtv vtu - ko(l + Xe(x,y))vu]d + f ( p iko) vu dr o r 1 B w/2 2 - - iko f x=-w/2 L y=o w/2 Fi dx = ikoy= x=-w/2 y=d dx + ko xe(x,y)vui dQ Qh (6.8) The remainder of the solution of this example closely parallels the procedure of Chapters III and IV. However, in order to calculate the far scattered field and the radar cross section, an integral such as (2.35) cannot be used because of the presence of the dielectric

-109 - cylinder. Instead, an integration of the scattered field values on a closed contour outside and surrounding the dielectric cylinder must be evaluated. If the contour is a circle of radius p = p as shown in Fig. 6.1, the far scattered field is ESC(p = 2pexp[ikp - i(T/4)]P(,o). (6.9) where, from the integral representation of the field, kP( o i aEz 2 ES (P,) ) P(,O) 4 4 ko Dp z 2 J exp[-ikop cos(+ - +)]P2 d' (6.10) 2 In the computer program, the radius p was taken to be slightly larger 2 than w, or more precisely, one element width out from the dielectric cylinder. Although the closed contour could be on the surface of the cylinder, in general, for an inhomogeneous body with a less distinct boundary the contour must be outside the body. Further, for a body with a noncircular cross section, it is desirable to place a circular contour outside the body for ease of computation of the far scattered field. 6.2 Selection of an Incident Field In order to demonstrate the use of this formulation without extensive modification of the computer program used to solve the problem of Chapters II and III, an incident wave is selected so that

-110 - the scattered field has the symmetry properties (2.24) and (2.27). This allows the solution to be obtained in the upper half space as before. In general, a scattered field can be interpreted as resulting from electric and magnetic current sources induced by the incident wave in the scattering body. In the present case p = V everywhere so that no magnetic current sources will be induced. For an incident field which excites only E-wave fields as discussed in Chapter II, the only possible electric currents must be in the z-direction. In order to maintain Hsc(x,O) equal to zero for |x| > w/2 as implied by (2.27), the induced electric current distribution must be an even function of y. For example, the electric current element I(x,y )z produces an x-component of the magnetic field on the x-axis which is cancelled by an x-component of the magnetic field of a current element I(x,-y )z. The incident field must then have the property that Ei(x,y) = Ei(x, - y) in order that the symmetry properties (2.24) and (2.27) hold. Here the incident field chosen with this property consists of two plane waves Ez = - exp[-iko(x cos % + y sin o)] + - exp[-iko(x cos q ~ - y sin p )], (6.11) whose angles of incidence are o0 and -p0, respectively.

-111 - 6.3 Implementation of the Solution Procedure The selection of the particular incident field (6.11) leads to a symmetric scattered field which is obtained by solution in the upper half space only. The weak form integral statement (6.8) becomes [vtv vu - k2( + x (xy))vu]dQ+ j ( iko vu dr Q rB B i k w/2 0 7 vu 2 R(x)y x=-w/2 y=o jk w/2 r V i dx = 2 / y=o x=-w/2 y=o dx + k2 Xe(x,y)vui dQ 0 e (6.12) where Q refers to the region in the upper half plane bounded by rB and the line y = 0. Following the procedure of Chapter III, the element matrix coefficients are given by n KNm JN, -N ( n n ax ax ay ay (1 + x (x,y))N N N d m^n --- x ay ay ko Xeen eL ik0 NN 1 2 R(X)dr d+ 2 iko) re r and the incident field ter and the incident field terms are given by NmNn dr m n (6.13) fe m iko 2 e rp ' ~^J Nmu R(x) dr + k2 '2e,n xe(x'Y) Nmui d (6.14)

-112 - 6.4 Numerical Results With these modifications of the computer program the radar cross section and the phase of the scattered field are computed for w/xo = 0.25 and ER = 2. Figures 6.2 and 6.3 show that the bistatic scattering from the quarter wavelength width strip without the dielectric cylinder to be relatively independent of scattering angle. The magnitude and phase of the induced currents on the strip are relatively uniform. However, when the dielectric cylinder is present the scattering action becomes more complicated. On the one hand, there is scattering off the surface of the cylinder. Also the incoming field which does penetrate the dielectric, interacts with the strip and the cylindrical surface setting up standing waves inside the cylinder and producing a second scattered field which combines with that scattered off the cylindrical surface. In this process, the orientation of the strip relative to the scattering direction becomes more important. Figures 6.4 and 6.5 show bistatic scattering for three angles of incidence. From these figures it is clear that the scattering is primarily characterized by radiation from the leading and trailing edges of the strip. For a scattering angle of z 90 degrees interference between leading and trailing edge radiation produces a null in the radar cross section.

1 6 8 0 0 (dB) -8 __j -16 -24 180 (DEGREES) Fi g. 6. 2: a/A vs. ~ for bistatic scattering, = 00, R =!.0, w/xA 0. 25. 0 0 0

210 X X x --- x -— x X - x- x - - 180 - - -1 - RESISTIVE STRIP ONLY 150 120,) LU UJ 0 -30 -30 ki 0 30 60 90 120 150 180 + (DEGREES) Fig. 6.3: Arg P(,pO ) vs. ( for bistatic scattering, <j = 0~, R = 1.0, w/X = 0.25.

12 8 4 0 (dB) 0 -4 -8 -12 Fig. 6.4: In 0 30 60 90 120 150 180 + (DEGREES) G/A vs. f for bistatic scattering, R = 1.0, w/xn = 0.25, dielectric cylinder (Ep = 2) radius/xo = 0.25.

-116 - C,) LUJ CD Lu C-) Q) 4-) L.) Q~J II C) 4-) Q) (I) C) (0 4 -0 V) S — co LO ) N om CO i-I) I1) (Y) 0)

CHAPTER VII. CONCLUSIONS AND RECOMMENDATIONS 7.1 Overview of the Work The finite element method was applied to the solution of two-dimensional electromagnetic scattering from bodies in unbounded space. The technique combined Green's first identity and the differential equation to form a weak integral statement equivalent to that obtained by the weighted residual method. Galerkin's method was used in the selection of the weighting functions. The finite element solution region was extended outward so that its outer boundary lay in the far field region of the scattering body. The boundary condition on the outer boundary was evaluated lby using the asymptotic expression for the scattered field in the fair-field region. The finite element mesh consisted primarily of quadrilaterals with linear shape functions and was generated by a computer program. Two special cases were solved using this approach: a thin resistive strip and a resistive strip embedded in a dielectric cylinder. Special triangular elements were used at the edge of the thin strip for the case of zero resistivity to model the singularity occurring in the magnetic field. Numerical results for scattering by the resistive strip were compared with results of a moment method integral equation solution known to produce highly accurate values of radar cross section and scattered field in the far field region. Results -11 7 -

-118 - indicate the element densities in the radial direction and adjacent to the strip should be kept equal. For an element density of 16 elements per wavelength in both directions, a strip width of a half wavelength, and a normalized resistivity R = 1.0, the magnitude of the far scattered field differed by less than four percent from the moment method result at all angles of incidence, except at deep nulls in the radiation pattern. The phase of the scattered field differed by less than six degrees at all angles of incidence. For low resistivity strips, pronounced standing waves appeared in the computed values of the scattered field. These standing waves were apparently due to multiple reflections between the strip and the outer boundary. The standing waves led to large errors in the radar cross section and the far scattered field. Their presence in the results prevented an evaluation of the effectiveness of the special triangular elements in modeling the field at the edge of a perfectly conducting strip. Another effect occurring for low (non-zero) resistivity was the ill-conditioning of the system of algebraic equations resulting in highly inaccurate results for most angles of incidence. This effect can be partially offset by increasing the element density. The finite element method was applied to the case of the resistive strip embedded in the dielectric cylinder. No moment method results are available for comparison, but the numerical results obtained here appear to be reasonable.

-119 - 7.2 Limitations of the Method The formulations and the computer program used here have some limitations. First, all work was done assuming an incident plane wave which excited only an E-wave. The formulation for H-wave excitation has not been developed. Second, the boundary condition on the outer boundary (rB) is not accurate enough to prevent standing waves from occurring in the computed scattered field, particularly for low resistivity. Increasing the element density will improve the boundary condition and reduce the standing waves but at increased computation cost (see Appendix B). Third, the memory requirement for the computer system increases rapidly with the size of the scattering body due to the need to extend the finite elements into the far field region. As an example, for the resistive strip in free space, if a certain element density is to be maintained, the band storage requirement for the global matrix in the finite element method increases as the fifth power of the strip width for w/x > 1 but only as the second power in the moment method. Fourth, the system of linear algebraic equations becomes ill-conditioned for low (non-zero) strip resistivities. A partial solution to this problem involves increasing the element density. The second, third and fourth items above are either not limitations for a moment method solution or are less of a limitation. Because of this the finite element method is not competitive with the method of moments for scattering from a thin body in a homogeneous unbounded space. On the other hand, for a thick

-120 - penetrable body or a body in an inhomogeneous medium, the moment method may be less attractive. This may be due either to increased memory requirements or to difficulties with an integral equation formulation in an inhomogeneous medium. The finite element method has been demonstrated to be a feasible approach for these cases. 7.3 Recommendations for Further Work 1. In order that scattering due to arbitrarily polarized incident plane waves may be determined, the formulation for H-wave excitation should be developed. 2. The use of quadrilateral elements with higher order shape functions should be investigated. It is expected that the same accuracy of solution could be achieved with significantly fewer elements of a higher order. 3. Element and global matrix coefficients should be interpreted in energy and circuit terms. Such interpretations could lead to a nonreflecting boundary condition for the outer boundary of the finite element region as was done for a one-dimensional case in Appendix C. The work of Engquist and Majda (1977) on absorbing boundary conditions should also be considered in improving the boundary conditions. 4. Three dimensional scattering problems could be formulated using methods similar to those described in this work.

APPENDIX A. APPROXIMATE RELATIONSHIP BETWEEN ELEMENT MESH PARAMETERS In selecting the values of the mesh parameters that are to be used for a particular computer run it is useful to know approximately the relationship between the parameters. The number of elements in the radial direction M is directly proportional to the product of EL and Mx where Mx (called NUMELS in the computer program) is the number of elements adjacent to the strip and EL is an aspect ratio for the elements. The parameter EL is the ratio of the dimension of a quadrilateral element in the direction parallel to the strip to the dimension in the direction perpendicular to the strip. The relationship in more precise form is M = A(EL)(M ) + B, (A.1) y x where A and B are constants yet to be determined. Figure A.1 shows the results of a number of computer runs. From this graph it can be seen that B is approximately unity and A is actually a function of SMEC. The parameter SMEC is approximately the eccentricity of the largest ellipse of the mesh, which forms the outer boundary at p = p. Figure A.2 is a plot of the slopes of the lines in Fig. A.1 vs. SMEC. This shows that A - -0.83 log SMEC + 0.17. (A.2) 10 -121 -

-122 - 25 20 1,5 l0 5 0 10 20 30 40 (EL) Mx 50 Fig. A.: M vs. (EL)M. y x

1.4 1.2 1.0 A 0.8 0.6 0.4 M(4 0.05 0.07 0.1 0.15 0.2 log lo SrVEC 03 0.4' 0.5 0.6 Fig. A.2: AVvs. SMEC.

-124 -The desired approximate relation is then M y = A(EL)(M ) + 1 x (A.3)

-125 - APPENDIX B. TRANSMISSION LINE INTERPRETATION FOR NUMERICAL RESULTS B.1 The Finite Element Mesh as a Transmission Line The resemblance of the set of linear algebraic equations for the nodal values in the finite element procedure to a set of node voltage equations for a two-dimensional electrical transmission line structure suggests that the results of the finite element solution be interpreted in terms of a transmission line. The use of electrical networks for the solution of electromagnetic field problems is well established. A great deal of work has been done in the past on the subject of wave propagation through two-dimensional periodic structures (Brillouin, 1953). Although the finite element mesh in the present work is not uniform and not perfectly periodic, the transmission line characteristics are evident in the numerical results. As an aid in the transmission line interpretation, the finite element solution for the problem of a plane wave normally incident on a plane resistive sheet of infinite extent is provided in Appendix C. This problem is analogous to that of a uniform one-dimensional transmission line and the method of solution generally parallels that used for the two-dimensional problem which is the subject of Chapters III through V. It is found that the characteristic impedance zc of the transmission line is always less than the characteristic impedance of free space Z0. This means that the boundary condition

-126 - (C.5) acts as an imperfect termination for the line and that standing waves appear in the numerical results. Further, the one-dimensional problem shows the phase constant k for the line is always less than the phase constant of free space ko so that the phase of the computed scattered field is always less than that of the exact value, except at the resistive sheet. If the density of the finite elements is increased, Zc approaches Zo and k approaches ko. This reduces the standing waves and the phase errors. These results for the one-dimensional uniform line apply qualitatively to the two-dimensional line extending outward from the resistive strip. However, since the two-dimensional line is not uniform in the direction of wave propagation, the impedance mismatch will depend on the local density of elements at the edge of the finite element region. An additional consequence of the non-uniformity of the line is the dependence of the standing-wave ratio on the impedance of the strip. A larger standing-wave ratio is observed for strips of low resistivity. For the two-dimensional problem, backscattering for normal incidence most closely resembles the one-dimensional problem of Appendix C. The transmission line characteristics are clearly apparent for this case. For example, Fig. B.1 is a plot of the magnitude of the scattered field as a function of the distance from the strip. The results of a method of moments (MOM) calculation are plotted for comparison. Figure B.1 shows that increasing the number of elements (My) results in a lower standing-wave ratio. The element mesh

-127 - ESC 1E | \ \\ ^' \ x. 0.4 0. 3 0.5- 2 -0.1 0.2 0.4 0.6 0.8 1.0 1.2 0 Fig. B.1: Scattered field magnitude vs. distance outward from the strip center (o = 90~, w/x = 0.5, M = 6, x = 0). 0 o x

-128 - parameters for the curves of Fig. B.1 were selected so that the standing waves could be easily seen in the data; the selection rules given in Section 5.2 were not used. Figure B.2 is a plot of the phase of the scattered field versus distance from the strip with MOM curves for comparison. It is evident that the error in phase for the FEM results increases with distance. Increasing (My) provides a small improvement in the phase accuracy. B.2 Effect of Standing Waves The value of the radar cross section of the strip may be obtained directly from the scattered field values computed in the far field region such as at p = p, or it may be obtained by means of a computation of the induced current on the strip followed by an integration over the strip (2.38). The numerical results obtained using the first method are generally less accurate than those having the second method when standing waves are present. To examine the effect of the standing waves on the current calculation in the second method, assume that near the strip Esc = [A(x) + B(x)y]exp[i(koy + oe(x))] (B.1) where A, B and eO are respectively the magnitude of IEzCI, the slope of IEzCI, and the phase of Esc at the strip (y = 0). Using (4.66)

-129 - 150~ ARG ESC -0 -50~ -100~ / -]5o~ /y - 150~ -200 _ 0 Fig. B.2: Scattered i center ( 0 0.2 0.4 0.6 0.8 y/A 1.0 Field phase vs. distance outward from the strip = 90~, w/x = 0.5, X = 1.0, M = 6, x = 0). 90 ' 0 ' X

-130 - 1 + exp i + ta -Ako. (B.2) 0 O z Z0 Ak0 -Ak 0 To apply this to the current calculation at the center (x = O) of the strip for R = 0, A0 = 1, and normal incidence, the following values are obtained from Figs. B.1 and B.2: MOM: A = 1.0, B z -0.2, e = -180 degrees FEM: A = 1.0, B B 0.6, e = -180 degrees Then (B.2) yields MOM: ZoK = 2.001 ei(1*820) FEM: ZoK = 2.009 ei(5420) It can be seen from this that the error in the slope of IEzCI at the strip caused by the standing wave results primarily in a phase error in ZoKz. Such phase errors for current in each element adjacent to the strip enter into the integration of ZoKz (2.38) to obtain P(p,4 o). 0

-131 - APPENDIX C. COMPUTATION OF THE SCATTERED FIELD FROM PLANE WAVE NORMAL INCIDENCE ON AN INFINITE RESISTIVE SHEET C.1 Introduction In this appendix a one-dimensional problem closely related to the problem stated in Chapter II is discussed. An infinitesimally thin, infinite sheet lying in the y = 0 plane is immersed in free space. The sheet is a resistive sheet having boundary conditions given by (2.1) and (2.2). An incident plane wave is polarized so that its electric field is i i' -ik y^ = Ez = e z, (C.1) _- z where an e iwt time dependence is assumed. The direction or propagation of Ei is in the negative y direction. The scattered field Ezc propagates in the positive y direction and satisfies d2ESC z + k2E = (C.2) dy2 ~ Z At the sheet the boundary condition for a resistive sheet is, using (2.26) iR d ESC 2iR z E= + Esc (C.3) k0 dy + z z

-132 - and for a perfectly conducting sheet is E + ES = 0 Ez z (C.4) Since the scattered satisifed for any y field is an outgoing wave, a condition must be > 0: dE sc z *. ESC dy l z (C.5) C.2 Finite Element Formulation From this point on, let u = EsC(y). The solution for u is to be obtained in the region 0 < y < L. A one-dimensional version of Green's first identity for scalars u and v is L y=o v du dy = - dy2 L J dv du dy y dydy y=o F L V dul dy= (C.6) Use of (C.2), (C.3) and (C.5) yields L fl dv du dy - k2 dy dy o y=o L ik fvu dy - iko[vu]yL - 2R [v(u + u )] y=o y=o (C.7) for a resistive sheet.

-133 - For a perfectly conducting sheet (C.7) cannot be used since R is in the denominator. Instead, the use of (C.2) and (C.5) in (C.6) yields L y=o dv du dy - k2 dy dy o L y=o vu dy - iko[vuy=L - v i y=o =. (C.8) The value of (du/dy) that at y = 0 at y = 0 is unknown, however (C.4) indicates i u = -u (C.9) This condition is imposed later. In order to solve for u, the region 0 < y < L is divided into J equal length elements e with nodes n at the element boundaries (Fig. C.1). For each node there is an interpolation function which is unity at that node and zero at all other nodes. The scalar u is approximated by u(y) - u(y) I n=l N ()an n n (C.10) where a = u(y ) and I = J + 1 The function v is restricted to be v(y) = wm(y), m = 1,...,I (C.ll)

-1 34 - n=lI n= I I I 1- 1 -- I -- -k 1 1 I1 I-. --- e::l1 e=2 FL e=J y=o y=L Fig. C. 1: Finite elements.

-135 - where wm(y) is one of a set of weight functions yet to be specified. Using Galerkin's method, wm(y) N (y) m (C.12) Using (C.10), (C.11) and (C.12) in (C.7) yields L y=o dNm(y) dy d t dy { Nn (y) a n=l dy - k2 0 L y=o I Nm(Y) l> n=l Nn(y)an dy -ik Nm(Y) I '1= 11 =. Nn (Y)an y=L ik F 2R Lm 6.~ I Nn Y an n=l y=o ik 2R [Nm(y)u]ly=o 2R m m = 1,...,I (C.13) Interchanging the order of integration and summation, (C.13) is written as I n=l K a mn n f m m = 1,...,I, (C.14) where Kmn L =y=o LdN m dy dN n - k2N N dy o m n.] dy - ik [NmN ]yL o m n y=L ik -R [NmNnly=O 2Ryo (C.15) and

-136 - ik i R u [Nm]y (C.16) m 2R 1 my=o Equation (C.14) represents a set of linear algebraic equations to be solved for the unknown a 's. Calculation of the coefficients K in the finite element mn method is done by an element-by-element procedure. The global integration over the entire region is replaced by a sum of integrals over individual elements lying within the interval. Thus J Kmn = Ken ' (C.17) e=t where Kdm n k2N N dy - ik [N N I - [ N N ] mn =J dy dy o m n o mI n y=L 2R m y=o e L J (C.18) For a perfectly conducting sheet, the set of algebraic equations should be modified by removing the equation corresponding to the node at the sheet where the field strength a is known. Using (C.9), a = -u (0) = -u. (C.19) 1 1 In this case (C.10) is written as I u(y) z u(y) -N (y)ul + N (y)a (C.20) n= n=2

-137 - Using (C.20), (C.11) and (C.12) in (C.8) yields I L N dy y -N (y) +)an dy - N (y)u dy dy 1 n n k m I y=o - n=2 y=o n n m 1 1 n n=2 n=2 y=L - N I (- N (y)ui + Nn(y)an) n=2 y=o = 0, m = 2,...,I (C.21) After interchanging the order of integration and summation, and transferring the term corresponding to the node on the sheet to the right-hand side, the resulting set of linear algebraic equations is I n=2 K a = m ' mn n m (C.22) where Kmn J -= e=l Ke mn (C.23) dN dN Ke = r m n - k2N N mn J dy dy o m n e dy - iko [NmNn] y=L (C.24) and gm = K u "m m1 (C.25)

-138 - C.3 Element Matrix Coefficient Ke mn Linear interpolation functions are used for the elements: Y - Yn Nn(y) 1 - Yn]n - Yn Y Y n * (C.26) n n-i In a typical element between y = Yk and Yk+i the interpolation functions can be written in a general form 2 where an = -1 for n = k and an = 1 for n = k + 1. The integral in Kmn (C.18) or (C.24) is now evaluated. dN dN ond = a k2 dm n k2N N dy = mn L 1 + 1 dy dy o m n L/( I - 1 3 mn ' eL (C.28) where m and n can each be "k" or "k + 1". Thus (C.18) becomes oilk2 Ke = ~m n o L (+ ) - ik [NmN mn LI - 1 3mn o mn y= ik o2R [NmNn]y=o (C.29) 2R y o

-139 - and Ke for the perfectly conducting strip case is given by mn k2 Ke = mOn ( L )(1 + 1 GtmOn [ y= mn L/(I ) 4 I - 1 4 1 -iko[NmNny=L (C.30) For a general element between nodes k and k + 1 (element not adjacent to y = 0 or y = L) element matrix coefficients are k2 e I 1 o L k,k - L 3 I - 1 k2 e - I 1 o L k,k+li L - 6 I - 1 (C.31) k2 e I - 1 o L k+,k L 6 I- 1 k2 e - I 1 o L k+l,k+i L 3 - 1 For the element adjacent to y = 0 (e = 1) the coefficients are k2 ik, I- 1 o L o k2 K -I- L k2 K - I -1 o L 21 - L 6 I- 1 k2 iK I 1 o L 22 L 3 I- 1 22^

-140 - For the element adjacent to y = L (e = J) the coefficients are k2 J I - 1 L I-iI-1 L - 3 I - 1 k2 K I- 1 o L KI-lI L - 6 I - 2. (C.33) J I - 1 o L KI-1 L - 6 I - 1 k2 J I - 1 L -k KII IL - 3 I - o In the case of a perfectly conducting strip, the coefficients for the elements are the same as (C.31), (C.32) and (C.33) except that K1 in (C.32) does not have the term with R in it. 11 C.4 Set of Linear Algebraic Equations Using (C.17), the coefficients for the set of equations (C.14) are written. If m f 1 and m f I, the mth equation is Km-a + + Km )a+ Km a 0 (C.34) m,-1 m-1 mm mm m mm+1 m+1 If m = 1, the equation is ik K1 a + K1 a 2 u (C.35) 11 1 12 2 2R 1 and if m = I, the equation is

-141 - KI-1 aIII-l I!-1 + KIaI = 0 (C.36) As an example, if there are six nodes (I = 6) and five elements (J = 5) the set of equations in matrix form is K 1 K1 21 0 0 0 0 O K1 12 (KI+K2 ) 22 22 K2 32 0 0 0 0 K2 23 (K2 +K3 ) 33 33 K3 43 0 0 0 0 K3 34 (K3 +K4 ) 44 44 K4 54 0 0 0 0 K4 45 (K4 +K5 ) 55 55 K5 65 0 0 0 0 K5 56 K5 66 0 — a 1 a 2 a 3 a 4 a 5 a 6 ik 2~ ui 2R 1 0 0 0 0 0 (C.37) for a resistive sheet. For a perfect conductor (C.22) and example, with I = 6 and J = 5, the system (C.23) are used. For of equations is (K1 +K2 ) 22 22 K2 32 0 0 0 O K2 23 (K2 +K3 ) 3 3 33 K3 43 0 0 0 K3 34 (K3 +K4 ) 44 44 K4 54 0 0 0 K4 45 (K4 +K5 ) 55 55 K5 65 0 0 0 K5 5s6 K5 66 a 2 a 3 a 4 a 5 a 6 ul K1 21 0 0 0 0 (C. 38) C.5 Transmission Line Interpretation The set of linear algebraic equations (C.14) may be interpreted in terms of a transmission line made up of discrete

-142 - components. A typical equation such as (C.34) resembles a node voltage equation if the coefficients Kmn are divided by -col. Each coefficient is then an admittance. The equation for node n is i n-1 1 [Kn-i n, n i K a + [K+Kna + a + i K a 1o n,n-i n-i L -o nn Knn n ( o n,n+i an+i = 0 (C.39) The node voltage equation for node n written in terms of the circuit components (Fig. C.2) is -Y1 a + (Y + (Y + Y3) a - Y3a+ -Y n-1 i 2 3 n 3 n+l = 0. (C.40) Upon comparing (C.39) and (C.40), it is apparent that Y = i n-1 1 ~ to n,n-1 i y i n+i Y - K 3 (i)0 n,n+l ' y + y + = (Kn + K ) 1 2 3 cO p0 nn nn (C.41) (C.42) (C.43) and From these Y = (Kn-i + Kn + Kn2 Io nn nn n,n-i 0on + Kn+l ) n,n+i (C.44) Using (C.31) yields

-143 - Fig. C.2: Transmission line.

-144 - k2 ] Y =Y i[ Ll+k2 - (C.45) 1 3 o Lt 6 1- I-1 and Y2 - C( -C i (C.46) The characteristic impedance of a transmission line made up of T sections (Fig. C.3) is (Johnson, 1950, p. 42) Z2 Zc =+ (C.47) Z 1 2 4 Using (C.45) and (C.46) yields Z = - L (C.48) 1 L and Z (C.49) 2 -1) Also k2. Z 1 - 1B (C.50) 0 ~ 1 1 - 4 1 + 0 +6 (I - 1 6 I L 1 J for which lim I+ Zc = - (C.51) 0

-145 - Fig. C.3: Symmetrical T network.

-146 - Thus the characteristic impedance of the finite element transmission line is always less than the free space value, but approaches 377 ohms as the density of elements is increased. Continuing the transmission line analogy, the proper termination of the line to prevent reflections and standing waves may now be determined. The transmission line (Fig. C.2) is terminated in its characteristic impedance Zc (Fig. C.4). The node voltage equation for node I is -Y a, I + +Y + aI = 0. (C.52) i I-i L 1 2 Zi I 2 c The finite element equation for the last node is,-K a + iK a = 0. (C.53) PO III -1 o II O Comparing (C.52) and (C.53) indicates that Y = - K( (C.54) 1 1,1 -o KP II-1 and y + Y + 1 = i K * (C.55) 1 2 Z' II 2 c

Fig. C.4: Terminated transmission line.

-148 - Then using (C.45) and (C.46) in (C.55) yields I = LI - 1 L 2 L o (C.56) 1+ Using (C.48) and (C.50) in (C.56) yields J 1 o L o L L 3 I 2 I - 1 L/(I -2 1 ) 77 1 70] 0 (L 0 - 4[ L )2]2 L/( I -1) |4 1 +(C.57) Use of K11 given by (C.57) corresponds to a transmission line terminated in its characteristic impedance. This expression should be compared with that given in (C.33) for K. The value K I of (C.57) approaches K I of (C.33) as I Ha. Thus use of KII in (C.33) corresponds to having an impedance mismatch at the termination of the transmission line. The greater the density of the elements, the less the mismatch is. The mismatch will cause standing waves in the numerical results for the scattered field. On the other hand, use of K( of (C.57) always leads to numerical results without standing waves. The propagation constant k for the transmission line is given by (Watkins, 1958, p. 28)

-149 - LY cos k - 1/ = 1 + 2Y (C.58) Use of (C.45) and (C.46) yields 2 cos k I - 1 2 or k k ~ (C.59) for (L/I-1) << x0. Thus the propagation constant is always less than that of free space k0, and the phase angle of numerical results will always be less than that of the actual scattered plane wave. The greater the density of elements, the less the phase error will be. Sample numerical results are given in Fig. C.5 through C.7 for various values of the element density (I-l)x o/L in elements per wavelength. The magnitude of the scattered field is plotted versus distance from the resistive sheet in Figs. C.5 and C.6. Figure C.5 shows the standing waves resulting from the imperfect termination of the transmission line and Fig. C.6 shows the result when the line is terminated in its characteristic impedance. The phase of the

-150 - 0.42 0.40 0.38 0.36 D sc i Z 2 0.34 0.32 0.30 0 0.2 0.4 0.6 0.8 1.0 Fig. C.5: j ESC vs. y/x for resistive sheet (R = 1.0), x = 1.0 and "free space termination" (Kij given by C.33).

-151 - 0. 42 0.40 0. 38 & lw 46 0. 36 0. 34 0. 32 4 I - - - & - - -.0 & - - — v & - - -..w & - - - 0 - - - 0 - -.- - '119 >( X >< --- X - X K - x x - -0 3 S & * p 0- I -13. 33 0- - 6. 67 x -—,X9. 33 0- 0 16. 68 EXACT I I I- I f'A 11 n U. ju 0. 2 0. 4 0.6 0. 8 1.0 YX0 Fi g. C. 6: I slvs. y/X for resistive sheet (R = 1.0), AC = 1.0 z ~ 0O and "characteristic impedance termination" (K 1 given by C.57).

-152 -scattered field versus distance from the sheet is shown in Fig. C.7 for the imperfect termination. The location of the exact values, if plotted, would essentially coincide with the line for 16.68 elements per wavelength. A set of phase curves for the characteristic impedance termination differs only slightly from Fig. C.7.

-153 - 300 200 li 100 a) w 0 I cSI Q -100 -200 0 0.2 0.4 0.6 0.8 y/I 1.0 Fig. C. 7: Arg ESC vs. y/X for resistive sheet (R = 1.0), X = 1.0 and "free space termination ( and "free space termination" (KII given by C.33).

BIBLIOGRAPHY Arlett, P. L., Bahrani, A. K., and Zienkiewicz, 0. C. (1968), "Application of Finite Elements to the Solution of Helmholtz's Equation," Proc. IEE 115, 1762-1766. Bettess, P. (1977), "Infinite Elements," Int. Jour. Num. Meth. Eng. 11, 53-64. Bettess, P. and Zienkiewicz, 0. C. (1977), "Diffraction and Refraction of Surface Waves Using Finite and Infinite Elements," Int. Jour. Num. Meth. Eng. 11, 1271-1290. Blackburn, W. S. (1973), "Calculation of Stress Intensity Factors at Crack Tips Using Special Finite Elements," in the Mathematics of Finite Elements and Applications (Whiteman, J. R., ed.), Academic Press. Bouwkamp, C. J. (1946), "A Note on Singularities Occurring at Sharp Edges in Electromagnetic Diffraction Theory," Physica 12, 467-474. Bowman, J. J., Senior, T.B.A., and Uslenghi, P.L.E., editors (1969), Electromagnetic and Acoustic Scattering by Simple Shapes, North-Holland Publishing Company. Brebbia, C. A. and Walker, S. (1980), "Boundary Element Techniques in Engineering," Newnes-Butterworths, London. Brillouin, L. (1953), Wave Propagation in Periodic Structures, Dover Publications, New York. Chang, S. K. and Mei, K. K. (1976), "Application of the Unimoment Method to Electromagnetic Scattering of Dielectric Cylinders," IEEE Trans. on Antennas and Propagation AP-24, 35-42. Engquist, B. and Majda, A. (1977), "Absorbing Boundary Conditions for the Numerical Simulation of Waves," Mathematics of Computation 31, 629-651. Fawzi, T. H. and Burke, P. E. (1974), "Edge Effects in Induction Problems," IEEE Trans. on Magnetics MAG-10, 429-430. -154 -

-155 - Franz, W. (1948), "Zur Formulierung des Huygensschen Prinzips," Z. Naturforschung, Part 3a, 500-506. Franz, W. (1949), "Zur Theorie der Beugung," Z. Physik 125, 563-596. Hanson, M. E. and Petschek, A. G. (1976), "A Boundary Condition for Significantly Reducing Boundary Reflections with a Lagrangian Mesh," Jour. Comp. Phys. 21, 333-339. Harrington, R. F. (1968), Field Computation by Moment Methods, MacMillan Co., New York. Jeng, G. and Wexler, A. (1977), "Isoparametric, Finite Element, Variational Solution of Integral Equations for Three-Dimensional Fields," Int. Jour. Num. Meth. Eng. 11, 1455-1471. Johnson, W. C. (1950), Transmission Lines and Networks, McGraw-Hill, New York. Jones, D. S. (1956), "A Critique of the Variational Method in Scattering Problems," IRE Trans. on Antennas and Propagation AP-4, 297-301. Klein, C. A. and Mittra, R. (1975), "An Application of the 'Condition Number' Concept to the Solution of Scattering Problems in the Presence of the Interior Resonant Frequencies," IEEE Trans. on Antennas and Propagation AP-23, 431-435. Knott, E. F., Liepa, V. V., Senior, T.B.A., (1973), "Non-Specular Radar Cross-Section Study," AFAL-TR-73-70, University of Michigan Radiation Laboratory Report 011062-1-F. Kouyoumjian, R. G. and Peters, L. (1965), "Range Requirements in Radar Cross-Section Measurements," Proc. IEEE 53, 920-928. Maue, A. W. (1949), "Zur Formulierung eines allgemeinen Beugungsproblems durch eine Integralgleichung," Z. Physik 126, 601-618. McDonald, B. H. and Wexler, A. (1972), "Finite-Element Solution of Unbounded Field Problems," IEEE Trans. on Microwave Theory and Techniques MTT-20, 841-847. McDonald, B. H. and Wexler, A. (1980), "Mutually Constrained Partial Differential and Integral Equation Field Formulations," in Finite Elements in Electrical and Magnetic Field Problems (Chari, M.V.K. and Silvester, P. P., editors), John Wiley and Sons, New York.

-156 - Mei, K. K. (1974), "Unimoment Method of Solving Antenna and Scattering Problems," IEEE Trans. on Antennas and Propagation AP-22, 760-766. Meixner, Josef (1949), "Die Kantenbedingung in der Theorie der Beugung Electromagnetischer Wellen an Vollkommen Leitenden Ebenen Schirmen," Ann. Physik 6, 2-9. Meixner, Josef (1972), "The Behaviour of Electromagnetic Fields at Edges,"' IEEE Trans. on Antennas and Propagation AP-20, 442-446. Morgan, M. A. and Mei, K. K. (1979), "Finite-Element Computation of Scattering by Inhomogeneous Penetrable Bodies of Revolution," IEEE Trans. on Antennas and Propagation AP-27, 202-214. Orlanski, I. (1976), "A Simple Boundary Condition for Unbounded Hyperbolic Flows," Jour. Comp. Phys. 21, 251-269. Poggio, A. J. and Miller, E. K. (1973), "Integral Equation Solutions of Three-Dimensional Scattering Problems," in Computer Techniques for Electromagnetics (Mittra, R., ed.), Pergamon Press. Sankar, A. and Tong, T. C. (1975), "Current Computation on Complex Structures by Finite-Element Method," Electronics Letters 11, 481-482. Senior, T.B.A. (1979), "Backscattering from Resistive Strips," IEEE Trans. on Antennas and Propagation AP-27, 808-813. Silvester, P. and Hsieh, M. S. (1971), "Finite-Element Solution of 2-Dimensional Exterior-Field Problems," Proc. IEE 118, 1743-1747. Smith, W. D. (1974), "A Nonreflecting Plane Boundary for Wave Propagation Problems," Jour. Comp. Phys. 15, 492-503. Stratton, J. A. and Chu, L. J. (1939), "Diffraction Theory of Electromagnetic Waves," Physical Review 56, 2nd Series, 99-107. Stratton, J. A. (1941), Electromagnetic Theory, McGraw-Hill, New York. Thatcher, R. W. (1978), "On the Finite Element Method for Unbounded Regions," SIAM Jour. Num. Anal. 15, 466-477.

-157 - Ungless, R. F. (1973), "An Infinite Finite Element," M.Sc. Thesis, Department of Civil Engineering, University of British Columbia, Canada. Watkins, D. A. (1958), Topics in Electromagnetic Theory, John Wiley and Sons, New York. Wexler, A. (1969), "Computation of Electromagnetic Fields," IEEE Trans. on Microwave Theory and Techniques MTT-17, 416-439. Zienkiewicz, (). C. (1977), The Finite Element Method, 3rd edition, McGraw —Hill, London. Zienkiewicz, 0. C., Kelly, D. W., and Bettess. P. (1977), "The Coupling of the Finite Element Method and Boundary Solution Procedures," Int. Jour. Num. Meth. Eng. 11, 355-375.