Electromagnetic Scattering from Jet Engine Inlets Using Analytical and Fast Integral Equation Methods by Hristos Thomas Anastassiu A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1997 Doctoral Committee: Professor John L. Volakis, Chair Professor Linda P.B. Katehi Professor Robert Krasny Professor Kamal Sarabandi RL-951 = RL-951

CO Hristos T. Anastassiu 1997 All Rights Reserved

ABSTRACT Electromagnetic Scattering from Jet Engine Inlets Using Analytical and Fast Integral Equation Methods by Hristos Thomas Anastassiu Chair: John L. Volakis Analytical and numerical techniques are employed for the evaluation of the Radar Cross-Section (RCS) of jet engine inlets. Initially, the Mode Matching technique is used for the extraction of reference data associated with canonical engine-like structures. Of most importance in this work is the utilization of the inherent cylindrical periodicity of the geometry to reduce the CPU time and memory requirements of rigorous integral equation methods for the characterization of jet engines. Finally, the Adaptive Integral Method (AIM) is applied to the RCS calculations of large, noncanonical geometries. The low computational complexity and storage requirements of AIM, combined with the exploitation of the structure periodicity, are especially attractive properties that allow the RCS analysis of very complex and electrically large jet engine inlets for the first time.

Z 'a VTOV~ 7FOt2 /iCVOVV/... 11

TABLE OF CONTENTS DEDICATION........ ii LIST OF FIGURES............................... v LIST OF TABLES....................................... viii LIST OF APPENDICES.. CHAPTER I. INTRODUCTION............ ix 1 1.1 1.2 1.3 Motivation and Previous Research......... Methods of Analysis................. Format of the Thesis................. 1 2 5 II. THE MODE MATCHING TECHNIQUE...... 2.1 Background of the Method............. 2.2 Cylindrical Hub Termination............ 2.2.1 2.2.2 2.3 Periodic 2.3.1 2.3.2 2.4 Periodic 2.4.1 2.4.2 Mode Matching Formulation............. Numerical Results................... Array of Grooves Termination............. Mode Matching Formulation............. Numerical Results................... Array of Curved Blades Termination......... Mode Matching Formulation............. Numerical Results................... 7 7 8 8 14 17 17 20 23 23 31 34 34 35 47 53 III. THE ADAPTIVE INTEGRAL METHOD (AIM)...... 3.1 3.2 3.3 3.4 Introduction........................... Mathematical background.................... Application to Relatively Flat Surfaces............. Numerical Results and Parallelization Issues.......... 111

IV. EXPLOITATION OF THE CYLINDRICAL PERIODICITY OF THE JET ENGINE GEOMETRY.............. 4.1 4.2 4.3 4.4 Introduction........................... Preliminary Setup for Perfectly Conducting Scatterers.... Decoupling of the Linear System of Equations......... Extension to Dielectric Scatterers................ 67 67 70 74 80 V. APPLICATION OF AIM TO THE JET ENGINE PROBLEM 84 5.1 5.2 5.3 5.4 Introduction.......................... The Cylindrical Waveguide Green's Function......... Application of the Slicing Scheme to AIM........... Numerical Results.......................... 84 85 89 97 VI. SUMMARY, CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK.............................. 108 6.1 Summary and Conclusions............... 6.2 Suggestions for Future Work.............. APPENDICES................................. 108.... 110....112 BIBLIOGRAPHY..................................148 iv

LIST OF FIGURES Figure 2.1 Cylindrical inlet terminated by a cylindrical hub................ 9 2.2 RCS (dB/A2) pattern calculated by the Mode Matching Technique (solid line) and the Hybrid Modal Technique (dashed line) for the geometry shown in Fig. 2.1. The specific dimensions for this calculation are: a = 1.5A,b = 3A,11 = 16.595A,12 = A a) 0 polarization, b) 0 polarization.................................... 15 2.3 Amplitudes of the elements of the generalized scattering matrix calculated by Mode Matching for a hub termination with a = 0.503A, b = 1.66A,11 = 16.595A,12 = 0.335A..................... 16 2.4 Cylindrical inlet terminated by an array of straight blades. 18 2.5 Tangential fields at both sides of the interface for the geometry corresponding to Fig. 2.4. The fields are plotted at p = 3A. a) p component, 0i = 0~, b polarization b) 4 component, 0i = 30~, 0 polarization....................................... 22 2.6 Data for an inlet terminated by an array of straight blades (6 and 8 G H z)................................... 24 2.7 Data for an inlet terminated by an array of straight blades (10 and 12 GHz)......................................... 25 2.8 Cylindrical inlet terminated by an array of curved blades....... 26 2.9 Tangential fields at both sides of the interface for the geometry corresponding to Fig. 2.8 (qb polarization). The fields are plotted at p = 3A. a) p component for Oi = 0~. b) q component for 0% = 30~.. 33 3.1 The Rao-Wilton-Glisson (RWG) basis functions........... 37 v

3.2 Definition of the far field region: fn lies in the near zone, whereas fp in the far zone of f..................................................... 40 3.3 Uniform AIM grid surrounding the scatterer...................... 41 3.4 Grid surrounding a Rao-Wilton-Glisson basis function........... 42 3.5 The local coordinate systems on triangles T~................... 44 3.6 Scatterer shape for various values of the flatness parameter y..... 50 3.7 Memory required for standard MoM and AIM (surface problems with expansion order M = 3, grid step h =: I and rthr = 101).......... 51 3.8 CPU time for standard MoM and AIM (surface problems with expansion order M = 3, grid step h = I and rthr = 101)............ 52 3.9 A 4A by 0.25A flat blade illuminated by a plane wave......... 54 3.10 Backscatter for the 4A by 0.25A flat blade of Fig. 3.9 (0 polarization, 0 = 0 cut). The blade lies on the xy plane with the long dimension along the x axis. The origin lies at the center of the blade. Grid step h = 0.05A......................................... 55 3.11 Backscatter for the 4A by 0.25A flat blade of Fig. 3.9 (0 polarization, = () cut). The blade lies on the xy plane with the long dimension along the x axis. The origin lies at the center of the blade. Grid step h = 0.05A................................. 56 3.12 A simplified aircraft structure..................... 59 3.13 Backscatter of the structure in Fig. 3.12 (~ polarization, 0 = 90~ cut).................................... 60 3.14 Backscatter of the structure in Fig. 3.12 (0 polarization, 0 = 90~ cut). 61 3.15 Profile on a CONVEX Exemplar for the scatterer analyzed in the previous figures............................ 62 3.16 The jagged plate geometry....................... 63 3.17 Monostatic RCS for the jagged plate geometry (0 polarization)... 65 3.18 Monostatic RCS for the jagged plate geometry (0 polarization)... 66 vi

4.1 Illustration of the computational decomposition. 4.2 Simplified model of a jet engine......... 69 4.3 Symmetry in a DBOR body............ 5.1 Decomposition of the computational domain. 5.2 Local rectangular grid for the basic slice.. 5.3 Cylindrical grid (front view)........... 5.4 Rotated cartesian coordinate systems...... 5.5 A quarter of a hub geometry........... 5.6 Backscatter for the hub geometry of Fig. 5.5 ( cut)....................... 5.7 Backscatter for the hub geometry of Fig. 5.5 ( cut)....................... 5.8 A quarter of a blade geometry.......... 5.9 Backscatter for the blade geometry of Fig. 5.8 ( cut)....................... 5.10 Backscatter for the blade geometry of Fig. 5.8............. 70............. 71............. 88........... 91........... 91............. 96............ 100 * - q) polarization, 0 polarization,....... (O polarization, polarization.. (0 polarization, 101 102 103 104 0-=0 0=0 cut).................................... 105 5.11 Backscatter of zero order propagating modes in a shorted inlet of 2A diameter (O polarization, 4 = 0 cut)................. 5.12 Backscatter of zero order propagating modes in a shorted inlet of 2A diameter (0 polarization, 0 = 0 cut).................. B.1 Transformation to a right isosceles triangle.............. 106 107 139 vii

LIST OF TABLES Table 3.1 Speedup due to FFT parallelization for the scatterer of Fig. 3.12. 59 3.2 CPU time comparisons for the scatterer of Fig. 3.16 with zero thickness................................. 64 3.3 Memory comparisons for the scatterer of Fig. 3.16 with zero thickness. 64 3.4 Memory comparisons for the scatterer of Fig. 3.16 with non-zero thickness.................................. 64 viii.. Vlll

LIST OF APPENDICES Appendix A. MATHEMATICAL DETAILS OF THE MODE MATCHING TECHN IQ UIU E..................................... 113 A.1 Explicit expressions of the modes in the cavity........ 113 A.2 Closed Form Expressions for the Coefficients of the Mode Matching System of Equations (Cylindrical Hub Termination) 118 A.3 Compact Expressions for the Coefficients of the Mode Matching System of Equations (Groove Array Termination)..... 121 A.4 The Eigenvalue Symmetry Property of the G.T.E. Matrix.. 124 A.5 Application of the Mode Matching Technique to the Generalized Modes within the Compressor............... 126 A.6 The Exponential Matrix Solution of the Generalized Telegraphist's Equations........................ 130 A.7 Explicit Radiation Coefficients for Circular Inlets....... 13-1 B. MATHEMATICAL DETAILS OF THE ADAPTIVE INTEGRAL METHOD (A IM ).................................... I 3z4 B.1 The AIM Impedance Matrix.................. 134 B.2 The Moments of the Basis Functions.............. 137 C. MATHEMATICAL DETAILS OF THE PERIODICITY EXPLOITAT IO N..................................... 143 C.1 The Cylindrical Waveguide Dyadic Green's Function..... 143 C.2 Elements of the Impedance Matrix for PEC Scatterers.... 145 C.3 Scattered Field from PEC DBOR's............... 146 ix

CHAPTER I INTRODUCTION 1.1 Motivation and Previous Research One of the most challenging problems in modern computational electromagnetics is the scattering characterization of aircraft configurations. Our interest in such an analysis is due to the fact that reliable target identification, a very important task in both military and civilian applications, is based on accurate Radar Cross-Section (RCS) calculations. Among the various contributors to the RCS, measurements have shown that the jet engine inlets have very significant impact to the overall signature of the airplane structrure. As a result, over the years numerous studies have been initiated for the analysis of jet engine inlets and their components. Numerous studies have been published on the analysis of propagation through ducts of rectangular and cylindrical cross-sections. Among them, Johnson and Moffat [1] presented a rigorous Wiener-Hopf analysis for the circular straight duct termi — nated by a short, and Ching-Chao Huang [2, 3] considered a GTD/modal analysis of circular ducts terminated by a simple blade configuration. The case of coated circular waveguides was investigated in [4] while the relative significance of the modes within the cavity was studied in [5]. In the case of large size, arbitrary cross-section ducts, high frequency (asymptotics) methods have been introduced. Among them, 1

2 the Shooting and Bouncing Ray (SBR) method [6] has been succesfully employed to track the ray fields within the duct. The SBR decomposes the aperture fields into a set of parallel rays which are then tracked into the duct in accordance with Geometrical Optics ((GO0) laws. An important advantage of the SBR is its simplicity and ease of interfacing with solid modeling packages [7]. However, although the SBR is attractive for general-purpose calculations, it lacks accuracy because it neglects rim diffraction and fails at caustics. Some of the disadvantages of the SBR method can be overcome by subdividing the aperture into smaller subapertures, and this is the basic concept of the Generalized Ray Expansion (GRE) method [8, 9]. Unlike the SBR, the rays are not necessarily parallel to each other, and consequently the GRE is capable of tracing non-planar wavefronts. Finally, Lu and Bansal [10] proposed a Boundary Relaxation Method (BRM) to numerically evaluate the model reflection coefficients at the duct termination. Although the above ray methods have been quite successful in modeling the field propagation through large cross section ducts, they are not suited for characterizing ray reflections from complex jet engine terminations which have irregular boundaries and edges. With this in mind, our motivation in our research is to present methods for the characterization of irregular engine-like terminations placed at the back end of the duct. 1.2 Methods of Analysis The jet engine problem is considered extremely difficult and challenging, mainly due to two reasons: large size and geometrical complexity. A typical jet engine inlet has a diameter of about 50 wavelengths in the X band, thus ruling out the use of any brute-force standard numerical techniques, which are typically suited for

3 3 - 5A structures on non-distributed platforms. On the other hand, high frequency techniques may accurately model the field propagation within the hollow duct, but fail to characterize the termination region which includes the electrically small details of the engine blades. Modal methods can treat some canonical blade configurations, but cannot be applied to realistic engines with more complex termination geometries. So far no single method has been adopted for the characterization of the entire engine, therefore hybridization of several techniques is absolutely essential for a complete analysis. One of the most serious drawbacks encountered during our investigation of the problem was lack of reference data. All previous analytical or numerical results were associated with either overly simplistic or too small geometries. Specifically, no data had ever been published regarding the presence of engine blades inside the cylindrical inlet. Since complex computer codes based on the Finite Element Method (FEM) were the original main focus of our research [11, 12], our need for code validation urged us to actually develop reference results on our own. The first objective of the thesis is to develop a rigorous method applicable to scattering from engine-like canonical structures.. The purpose of this effort was to develop an understanding of the scattering mechanisms within the inlet (coupling patterns among modes, dependence of fields on blade geometry, e.t.c.), as well as to calculate RCS results of some canonical configurations that can be used as reference data. Since rigor and best possible accuracy was our goal, modal techniques was our choice. Specifically, the Mode Matching (MM) method [13, 14] was applied to three different cavity terminations (cylindrical hub, array of straight blades, array of curved blades). The results were found to be in excellent agreement with the pertinent measurements, and therefore the developed MM code has since been an

4 indispensable tool for the validation of other, more versatile numerical techniques. Our calculated RCS patterns for inlets terminated with engine-like structures were published for the very first time in the literature [15, 16]. Due to its limitations, the MM method can only be applied to canonical geolnetries, where analytical expansions of the fields in terms of modes are possible for some critical sections of the termination. On the contrary, the geometry of realistic engines is irregular, and not amenable to modal representations. Furthermore, even among various canonical geometries, the MM technique must be reformulated for each particular case. Also, memory and CPU time considerations for modern computer facilities restrict its performance to inlets smaller than 20 wavelengths in diameter. To address the limitations of the MM code, a Finite Element (FEM) code was developed [11, 12], but because of its volumetric nature, this approach also has exorbitant memory requirements. Finally, standard integral equation methods, such as the Moment Method (MoM), is completely out of question, for the same reasons. Our need for versatile numerical techniques with low computational complexity urged us to investigate a few of the so-called "fast integral equation methods". Specifically, the Adaptive Integral Method (AIM) [17, 18] is given particular attention, due to its very low computational requirements for flat surfaces, such as jet engine blades. The application of AIM to the problem of scattering from jet engine inlets was extensively studied leading to very successful simulations of the jet engine. Originally, a free space AIM analysis was developed and validated, and the method was subsequently applied to field calculations inside the engine. Several fundamental modifications of the method were necessary before the technique could be used within a circular waveguide geometry. It is shown that very significant reduction of computational requirements can be achieved over the expensive Moment Method (MoM)

5 without accuracy compromise. Crucially important to this work is the exploitation of the jet engine cylindrical periodicity to reduce the computational domain down to a single periodic sector ("slice") of the structure. The versatility and inherently low computational complexity of AIM, combined with the "slicing scheme" enables us to analyze electrically large, realistic engines for the first time. 1.3 Format of the Thesis Chapter II presents the application of the Mode Matching Technique to a number of canonical engine-like structures. Three different types of termination are considered, and in the case of curved blades a novel analysis based on generalized modes is performed. Several comparisons with measurements are shown, and the excellent agreement establishes the method as a reliable source of reference data. The basics of the Adaptive Integral Method (AIM) are discussed in Chapter III. The method is presented in a mathematically rigorous way, along with considerations regarding its application to relatively flat surfaces. Numerical results and comparisons with the MoM show the impressive accuracy of the method. Parallelization schemes are also developed and the significant advantages of AIM over the MoM in terms of storage requirements and computational complexity are clearly demonstrated. In Chapter IV the circular periodicity of the jet engine geometry is exploited to reduce the effective size of the engine by a factor equal to the number of blades. Analytical results are derived for both perfectly conducting and dielectric materials. The immediate consequence of this analysis is a very significant reduction in memory and CPU time requirements of any integral equation method applied to the jet engine geometry.

6 In Chapter V AIM is directly applied to field calculations inside the cylindrical duct. The basic AIM concepts are substantially modified and a cylindrical AIM grid is introduced in conjunction with the "slicing scheme" for further storage and CPU time savings. Various resulting RCS patterns are shown, proving the validity and computational efficiency of the method. Several appendices describe pertinent mathematical details not given in the main text.

CHAPTER II THE MODE MATCHING TECHNIQUE 2.1 Background of the Method The Mode Matching (MM) method [13, 14] is an analytical technique applicable to the modeling of electromagnetic propagation through waveguide junctions. The fundamental requirement of the method is that the geometries of the waveguide sections be separable. In the associated coordinate systems the Helmoltz equation can therefore be solved analytically and the fields can be expressed as superpositions of modes (eigenfunctions). The unknown coefficients of the field expansions are determined by invoking the appropriate boundary conditions at the interface between the waveguide regions, i.e. continuity of the tangential components of both the electric and magnetic fields. Application of the latter conditions, in conjunction with the orthogonality properties of the eigenfunctions leads to an infinite linear system of equations, the solution of which yields the undetermined coefficients. A closed form solution is possible in a few simple cases [13], but in general, approximate, numerical solvers must be used, after careful truncation of the system. The method has been successfully applied to several microwave waveguide problems [19, 20, 21, 22, 23, 24], but usually the analyzed structures have been fairly small, supporting only few propagating modes. In this chapter we apply the method 7

8 to jet engine inlets, where a much larger number of propagating modes exists, requiring much more intense and challenging calculations. 2.2 Cylindrical Hub Termination 2.2.1 Mode Matching Formulation In this section we apply the MM technique to the geometry of Fig. 2.1. We assume that a plane wave is incident on the open end of the cylindrical duct and we are interested in evaluating the scattered field as a function of the incidence angle. In accordance with the MM technique, the fields in regions 1 and 2 are first expressed as a weighted sum of waveguide modes propagating either along the positive or negative z direction. Of course, the coefficients of the expansion are then determined by enforcing tangential field continuity across the interface separating regions 1 and 2. Before the application of the method it is necessary to evaluate the coupling of the incident plane wave with the forward propagating modes of region 1. The geometries of regions 1 and 2 are separable, and therefore Helmholtz's equa — tion can be solved analytically to yield the appropriate mode representations in each region [25]. These representations are given explicitly in Appendix A. As in Appendix A, the transverse fields of each mode are denoted by ei, h.,n where e stands for the electric field, and likewise h stands for the magnetic field. The subscript i may be either 1 or 2, corresponding to regions 1 and 2 respectively, and n is the characteristic integer of the mode. In region 1 the modes can be TM or TE, whereas in region 2 they are TM, TE or TEM. Using the mode functions given in Appendix A, the transverse components of the incident field are represented as a superposition of modes propagating along the -z direction, viz.

9 -. l,...:: 4 - 1z - O z=O. Figure 2.1: Cylindrical inlet terminated by a cylindrical hub 00 E E = anei exp{j31,nz} (2.1) n=l oo Ht = ahi,nexp {j/,31,z} (2.2) n=l can be found in terms of the incident plane wave as [2, 3]: where the an coefficients 7wrej kr an = Einc (0) Eadn (r) (2.3) jWI/t where r is the spherical far zone distance from the center of the open end of the cavity, Einc (0) denotes the incident plane wave at the center of the open end before it couples into the cylinder and Erad,n (r) is the radiated field of the nth outgoing mode of unit power impinging at the inlet mouth from the interior of the waveguide (see Appendix A). Eq. ( 2.3) can be proven via the reciprocity theorem [26], by using

10 a delta current test source at a far field location r. The transverse fields reflected at the interface take similar forms, namely 00 Er = bnel, exp {-jPz} (2.4) n=l oo H1r = - b,hl,exp{-j/l,nZ} (2.5) n=l where bn are unknown coefficients to be determined. The transverse fields in region 2 can be expressed as 00 E2 = cme2,m exp {j2,mZ} + m=l oo + E dme2,mexp {-j2,mZ} (2.6) m=l oo Ht = E cmh2,m exp {j:2,mz}m=l oo - S dmh2,mexp{-j32,mz} (2.7) m=l The latter are sums of +z and -z propagating modes, where the coefficients cm, dmi, bi must be determined by enforcing continuity of the tangential fields at the interface (z = 0) and the required boundary condition at the back of the guide (z = -/2). Mathematically, field continuity at z = 0 implies Et/l=o + Etrlzo = E=lzo (2.8) H ~tIz=o + Hlr I = H Iz=o (2.9) Vp c [a, b], V E [0, 27] and

11 Etolz=o + Ejlz=o = 0 (2.10) vp c [0, a], v C [0, 2r] To facilitate the analysis, we introduce the vectors {a} [al,...,aM]T (2.11) {b} - [bl,...,bM ]T (2.12) {c} - [cl,...,cM2]T (2.13) {d} - [dl,..,d(M2] (2.14) where Mi is the number of modes in the ith region. Evidently, for an exact modal representation, an infinite number of modes has to be used, i.e. for both regions M1,2 -+ o. Substituting ( 2.1), ( 2.4) and ( 2.6) into ( 2.8) and ( 2.10), taking the dot product of the resulting equations with elm, and then integrating over the cross-section of the interface yields [U] ({a} + {b}) = []T ({c} + {d}) (2.15) where the superscript T denotes the transpose of the matrix. The elements of the square matrices [P] and [U] are given by rb r27r Pmn -- e2,m e,,npdpdq (2.16) Ub 2,7r Umn n m / el,m, el,mpdpdoi (2.17) Jo Jo

12 in which 8mn is the Kronecker delta. Next, on substituting ( 2.2), ( 2.5) and ( 2.7) into ( 2.9), upon taking the dot product of ( 2.9) with h2,m and integrating over the portion of the interface which excludes the perfectly conducting disk defined by p < a yields [R] ({a} - {b}) = [V] ({c} - {d}) (2.18) where Rmn - j h2,m hl,npdpdoq (2.19) Vm n 6mn j h2,m h2,mpdpdq$ (2.20) Finally, upon enforcing the boundary condition at the back end (z =-12) we obtain the system [D]1 {d} = [Send] [D] {c} (2.21) where Dmn = mmn exp {-j32,m 12} (2.22) and [Send] denotes the scattering matrix at the back end of the termination. It can be shown that [P] [U]-1 = [V]-1 [R]- =[L] (2.23) provided the modes in both regions are normalized to unity power according to ( A.10) of Appendix A. Using this identity, it follows from ( 2.15), ( 2.18), ( 2.21) and ( 2.23) that

13 {b} = ([W] + [I])-: ([W]- [I]) {a}- [S] {a} (2.24) where [W] = [L] ([I] + [D] [Send] [D)([D ] [D] [Send] [D] ) [L] (2.25) and [I] is the unit matrix. The matrix [S] = ([W] + [I])1 ([W]- [I]) (2.26) is the generalized scattering matrix of the termination and is independent of the excitation. It is worth mentioning that the elements of [L] for this geometry can be evaluated in closed form (see Appendix A for the relevant analytical expressions). It is evident from ( 2.24) that the infinite modal scattering matrix must be truncated for a numerical evaluation of the vector {b} and extreme care is necessary to achieve an appropriately convergent result. Mittra and Lee [13] and Shih [14] discuss these difficulties, and point out that reliable convergence can be achieved by choosing an appropriate ratio M1/M2 of the number of modes used in the two regions. This ratio depends on the geometry, and it can be shown [27] that once the relative convergence problem is resolved by appropriately choosing the ratio MA1/M2, good conditioning of the system is guaranteed. Once {b} is determined, the far zone radiated fields from the open end are obtained by integrating the modal fields over the aperture (Kirchhoff approximation) and the appropriate expressions are given in Appendix A. Summing over all the modes, using the bm coefficients as weights, yields the total scattered field of the structure.

14 2.2.2 Numerical Results Figure 2.2 shows some numerical results for the geometry shown in Fig. 2.1. In generating the calculated data, 420 modes were used to represent the fields in each region, i.e. the order of the cylindrical functions ranged from n = 0 to n = 14 and 7 modes for each n were taken into account (i.e. we used 15 7 = 105 TM modes, plus 105 TE modes, multiplied by 2 for even and odd angular behavior. For region 2, the TEM mode substituted one of the less important evanescent TM modes). Clearly, the MM calculations are in good agreement with data based on the Hybrid Modal Technique [28], where the termination is modeled by the method of moments. It is worth noting that the MM results converge rather slowly, and a smaller number of modes may not yield correct values for the RCS. On the other hand, the inclusion of more modes requires the inversion of larger matrices which may be possibly associated with worse matrix condition numbers. Figure 2.3 depicts the amplitudes of the elements of the generalized scattering matrix [S] corresponding to the above geometry. It is evident that coupling exists only among modes of the same order (of the cylindrical function), and this can be proven analytically. Also, TM and TE modes do not couple to TM and TE modes (respectively) of different angular variation, e.g. a TM mode of cos no variation does not couple to a TM of sin no variation. However, intercoupling among TM and TE modes does occur only for modes having different angular variation (i.e. one of them having a cos no and the other a sin no variation). Moreover, the coupling between any two given modes is invariant if the angular variations of both modes are switched. Finally, all modes couple significantly only to their closest neighbors, resulting to almost banded scattering matrices. Large off-diagonal elements (whenever they occur) correspond to evanescent modes.

15 40.0 0o0 - 21,.1 O.l - U.o - 0.0 30.0 20.0 0.000 dIO.0 0 I 30.00 (dca) (a) -.......-.. — - -- I T~ -1k - -------- - r - 1-. 20.0 0.I -. -- - --- -. -... --— 0- -. X) IO.0 20.00 3x0 40,X 5.00 e, (demyes) (b) Figure 2.2: RCS (dB/A2) pattern calculated by the Mode Matching Technique (solid line) and the Hybrid Modal Technique (dashed line) for the geometry shown in Fig. 2.1. The specific dimensions for this calculation are: a = 1.5A, b = 3A, 11 = 16.595A, 12 = A a) 4 polarization, b) 0 polarization.

16 O. 75 6 0.5 0.25 40 020 40 Modes 1-10: TMo,1 to TMo,lo (even) Modes 11- 20: TM1,1 to TM1,10 (even) 60' 0 Modes 21- 30: TM2,1 to TM2,10 (even) Modes 31- 40: TEo,i to TEo,io (odd) Modes 41- 50: TE1,1 to TE,1io (odd) Modes 51- 60: TE2,1 to TE2,10 (odd) Figure 2.3: Amplitudes of the elements of the generalized scattering matrix calculated by Mode Matching for a hub termination with a = 0.503A, b 1.66A, 1' = 16.595A, 12 = 0.335A

17 2.3 Periodic Array of Grooves Termination 2.3.1 Mode Matching Formulation In this section we apply the Mode Matching (MM) technique to the more cornplicated geometry of a hollow cylinder terminated by a cylindrical array of grooves as shown in Fig. 2.4. Evidently, the modes in region 1 are identical to those of the previous geometry. To provide a modal representation of the fields in region 2, we will assume that the array consists of J grooves, that are essentially cylindrical (pie-shell) sectors of angular extent,. In each of the pie-shell sectors the wave equation can be solved analytically, leading to the modal solution given in Appendix A. The modes in Appendix A refer to the 'th groove whose rightmost edge makes an angle =- 1 + ( - 1) (2.27) with the x-axis. Note that no TEM mode exists within the grooves, since there is no seperation among the conducting surfaces of each groove. The fields in each region can be again represented by ( 2.1)- ( 2.7), with the understanding that e2,n and h2,n are different. To find the coefficients cm, dm, bm we proceed with the enforcement of the continuity conditions of the tangential fields at the interface (z = 0), viz. Ei=o + El=o= Etz=o (2.28) Hti =o + Htr =Htl o (2.29) Vp V [a, b], V< E [0, on + Ow],+ e {1,..., J} and to these we must add

18 Side view i... X k b ". a — (P K - ip N Front view Front view )w. I11,_ 12 Side view Figure 2.4: Cylindrical inlet terminated by an array of straight blades.

19 Etilz=o + E[z o = 0 elsewhere (2.30) i.e. on the metallic portion of the interface. Following the same mode matching procedure used for the previous geometry, ( 2.28)-( 2.30) can be used to determine the scattering matrix of the structure. As before, upon applying the MM technique at the interface we obtain a system of equations for the mode coefficients. Then, on eliminating the coefficients of the modes in region 2 (groove region) we deduce that ( 2.26) is again valid except that the matrix [W] must be replaced by [W] =: [Lk]T ([I] + [D] [Sn,end] [D] )( [I] - [D] [Sn,end] [D] ) [L,] (2.31) x,=1 in which [I] is the unit matrix and [S,end] is the modal scattering matrix of the back end of each groove. The elements of the [L,] and [D] matrices are given by L - PKrmn (2.32) Un b 0K +Ow PI,mn e2,^,m el,npdpdO) (2.33) fb '27V Un = ei,n e-l,npdpdc (2.34) Dmn = mnexp {-j/2,m12} (2.35) Clearly, the expressions derived for the scattering matrix of the periodic groove termination are quite similar to the corresponding ones for the cylindrical hub. However, apart from the sum in ( 2.31) over the total number of grooves, there are a few additional important differences. The order of the cylindrical functions in region 2 is, in general, non-integer and, of course, the +-integration is performed only over the angular extent of the groove. Consequently, in general, all modes between regions 1

20 and 2 couple to each other, leading to a fully populated matrix [L,]. In the case of the stub termination, orthogonality among the trigonometric functions resulted in coupling only among those modes of the same order (i.e. same order of the Bessel functions), leading to banded matrices. Nevertheless, the scattering matrix [S] must be banded [12], and this indeed occurs after carrying out the matrix operations speci — fied in ( 2.26),( 2.31). Finally, we note that as opposed to the single hub termination, the elements of the [L,] matrix must be evaluated by numerical integration. 2.3.2 Numerical Results When our method was first developed (Fall 1993), no reference data existed for this geometry, and to validate the computer code we computed and compared the transverse electric fields at both sides of the interface separating regions 1 and 2. Fig. 2.5 depicts the p (a) and X component (b) of the total electric field at both sides of the interface. The dimensions for this calculation are: b = 3.5A, a = 2.5A, 11 = A, 12 =, 4 blades,, = 45~, 1 = 0~. The waveguide field is plotted at p = 3A, as a function of A, for an incidence angle of O- = 0~ and q polarization (a) and for an incidence angle of Os = 30~ and 0 polarization (b). It is apparent that the transverse fields are continuous across the interface between regions 1 and 2, implying that the boundary conditions ( 2.28) and ( 2.29) are satisfied. Moreover, the total tangential field on the metallic surface of the blades is 20 or so dB below the field values elsewhere. The transverse field values on the top surface of the hub were also found to be very low, implying that the boundary condition ( 2.30) is also satisfied fairly well. Consequently, since all boundary conditions are satisfied, we were convinced that any RCS results that we had obtained so far must have been accurate within the limits of the Kirchhoff approximation employed at the aperture of the cylinder's

21 open end. A year later, actual measurements were collected for a cylindrical cavity ter — minated with an array of straight blades [29], finally enabling us to compare our calculations with the measured data. The model shown in Fig. 2.4 was constructed and measured. Following the notation used in the figure, the specifications of the fabricated structure were b = 15 cm, a = 7.5 cm, 11 = 20 cm (length of hollow section), 12 = 10 cm (length of termination), &~ = 40~, q1i = 2.5~, and the inlet was backed by a perfectly conducting plate. We investigated the backscatter in the b = 0 plane. The solid and (lotted lines shown in the following figures correspond to calculated and measured data, respectively. To give more details about the computations using the MM method, we define the following parameters: eN1,2: number of orders of cylindrical functions used in the hollow section (subscript 1) and the termination section (subscript 2). eM1,2: number of modes per order used in the hollow section (subscript 1) and the termination section (subscript 2) In generating the calculated curves shown in Figs. 2.6, 2.7, the above parameters were selected as follows: 6 GHz: N1 = 21, M1 = 10, N2 = 6, M2 = 8. 8 GHz: N1 = 25, M1 = 10, N2 = 7, M2= 8. 10 GHz: N == 30, M1 = 12, N2 = 7, M2 = 12. 12 GHz: N1 = 33, M1 = 13, N2 = 10, M2 = 8. In spite of the model's and code's complexity, as well as the high sensitivity of the patterns on the geometrical parameters, the agreement between measured and calculated data is, in general, very good. Most disrepancies are likely due to diffraction from the exterior back section of the cylinder (not included in the calculated

22 U. 0. fz 1-11 m -.1 -1 lu - Oi -10. -20. -30. -~, -40. -50. 20. 10. * (deg) (a) n - cidmt FPiid.. Tota Fid (Region I) -—... ToW Fed (Raei 2) Ca i 10:2 0. -10. -20. -30. -40. 0. 45. 90. 135. 180. 225. 270. 315. 360. * (deg) (b) Figure 2.5: Tangential fields at both sides of the interface for the geometry corresponding to Fig. 2.4. The fields are plotted at p = 3A. a) p component, 0i = 0~, $ polarization b) q component, Qi = 30~, 0 polarization.

23 data) or the Kirchhoff approximation of the field radiated through the aperture. The reliability of the Kirchhoff approximation increases with the aperture size, resulting in better agreement at higher frequencies. Finally, some inaccuracies stem from truncation of the infinite modal expansions. It should be pointed out that the memory and CPU time required by the MM method increase drastically with the diameter of the structure. Inlets wider than 20A in diameter require the storage and inversion of several full matrices, larger than 1500 x 1500, and therefore must be analyzed on a multiprocessor computing facility. 2.4 Periodic Array of Curved Blades Termination 2.4.1 Mode Matching Formulation None of the previously considered geometries represents a realistic jet inlet termination. With the goal of characterizing the scattering of more realistic jet inlet terminations, in this section we present a mode matching solution for an inlet termination consisting of a special class of curved blades (see Fig. 2.8) Basically, each blade is allowed to be curved with the restriction that its base still remains perpendicular to the hub at any point of intersection (see Fig. 2.8, front view). However, the angle formed by the rightmost edge of the boundary of the blade and the x-axis is allowed to vary with z. Specifically the blade or groove boundary is now allowed to make an angle Or (z) -= (0) + ( - 1) 7- + F(z) (2.36) with the x axis, where F is an arbitrary function of the longitudinal coordinate z2' and is measured in radians. Of course, at the face of the fan (z = 0), the condition

24 6 GHz Theta Pol. 20.0 - -- -" "'.C. —100 -- C 0 --- —-,0 ^.. -30,0 -40.0 - 6 GHz Phi Pol. 20.0 --, 1-1 IE C) col 10.0. --- --- --------- ---- ' ------ -20.0 - - -- -20.0...Af.l..........................tff 'l............ — I-.00 10.00 0;0 30.00 40.00 50.00 60.00 % (dgms).00 10.00 20.00 30.00 4,000 50.00 60.00 ei (dgrS) 8 GHz Theta Pol. 8 GHz Phi Pol. *5j S i s3.00 1O.00 20,00 30.0 40.00 50.00 60.00 0, (dogmt) e (degn s) Figure 2.6: Data for an inlet terminated by an array of straight blades (6 and 8 GHz).

25 100GHz Theta Pol. 10 G~lz Phi Pol. Li.00 20.00 20.00 30.00 40.00 50.00 60.00 0, (degree).00 I10.00 20.00 30.00 40.00 50.00 60ff 0j (depfts) 12 GHz Phi Pol, 120GHz Theta Pol C4 s CA la %..I' 6 04 'V U,) CZ.D0 [0.00 20.00 3D.DO 40.00 $0.00 60.00 9, (degrees).00 i., 10,00 20-00 30.00 40.00 50.00 M0.0 Si (degsr=e) Figure 2.7: Data for an inlet terminated by an array of straight blades (10 and 12 GHz).

26 Side view Perspective view b ' i!:! r.:... -~. pw, I1 12 (, Front view Side view Figure 2.8: Cylindrical inlet terminated by an array of curved blades. F() = 0 (2.37) should be satisfied, so that ( 2.36) is self-consistent. Note that ( 2.36) is valid provided all blades have the same z-dependence, resulting in grooves of constant width for any z; that is each pair of blades is assumed to form an identical guiding region, as is typically the case with realistic engine compressor configurations. To analytically characterize the fields within each groove formed by a pair of these curved blades, we note its similarities with the straight groove. We observe that for any given z the cross-section of the curved blade pair is identical to that of a straight blade pair, the only difference being the rotation about the z axis. Consequently, the

27 modes expressed in Appendix A still satisfy the appropriate boundary conditions, provided, is a function of z satisfying ( 2.36). It is therefore reasonable to assume that the transverse fields within the Kth groove (i.e. the waveguide region between the (K- l)th and the Kth blades) can be expressed as a superposition of the aforementioned modes, viz. 00 E; = E V, (z)e., (p, X; z) (2.38) i=1 00 Ht, = E I,,i(z)h i(p,; z) (2.39) i=1 In the latter expressions, the modes are normalized in such a way that Jij(z) e< e d2S = - 8 (2.40) f, hz ed2S = (2.41) I IS^{z) where S, (z) is the cross-section of the Kth groove at z. It can be shown [30] that the coefficients V, (z). I, i(z) satisfy the following infinite system of differential equations: V,1 (Z) T11 (z) T12 (z)... -jZ.. V,1 (z) V,2 (z) T21 (z) T22 (z)... 0 -j2Z2... V,2 (z) dz f,1^) Ij 0 -T1 (Z) -T21 (Z)... In, (z) In,2 (z) 0... -T12 (z) -T22 (z)... I,2 (z) (2.42) where

28 TJ (Z)- Z e,,j. d2S (2.43) 3i is the propagation constant and Zi is the impedance of the ith mode. Given F(z), analytical evaluation of a is feasible. Also, in ( 2.43) the integration over ( can be performed in closed form, but the integration over p must be done numerically. It is possible to solve ( 2.42) analytically, provided the geometry obeys certain restrictions. First we note that the cross-section of all grooves remains invariant in shape along the axis, and thus the propagation constants and the mode impedances do not depend on the longitudinal coordinate z. Therefore, ( 2.42) can be rewritten more compactly as dz {UK (z)} = [M (z)] {U, (z)} (2.44) where {U, (z)} - [M (z)] - [D] V,1 (Z) V,,2 (Z) Zl/I,1 (z) Z2I,,2 (z).) [D] -_ [p]T F(z) (2.45) (2.46 ) [D] = diag[-ji,3, -j32...] (2.47)

29 Pij - e.,. d2S (2.48) J (z) <93 Oo and the prime on F(z) denotes differentation with respect to z. From the explicit expressions of e, i in Appendix A, we conclude that Pij does not depend on z. Most importantly, the elements of [M] carry a dependence on z only through the presence of F'(z). That is if F'(z) = const., [M] becomes independent of z and this can be exploited to obtain a closed form solution of {UN (z)}. Specifically, it follows that if F(z) is linear in z, the explicit solution for {U, (z)} is {U, (z)} = exp (z [M]) {U, (0)} (2.49) where {U, (0)} is the value of {U, (z)} at z = 0. Although this is an explicit expression, the numerical evaluation of the matrix exponential is not necessarily an easy task. To evaluate it, the standard approach is to rewrite it as [31] ezM] = [K] diag [ez, zA2,..., ez] [K]-1 (2.50) where Ai denote the eigenvalues of [M] and [K] is the matrix of the corresponding eigenvectors [31]. The latter expression is valid provided the eigenvalues are distinct, which is expected to hold in this situation, unless degenerate generalized modes exist.. In theory, the matrix [M] is infinite, but for practical purposes it must be truncated, taking into account only the traveling and most significant evanescent modes. It is important to note that, as shown in Appendix A, if A is an eigenvalue of [M], then -A is also an eigenvalue of [M]. This is an expected property since each guiding region supports identical pairs of modes propagating along the +z and -z directions. By making use of this property, ( 2.49) can be rewritten as

30 V, (ZK) l 0 [Kl1] [K22] 0 [K] V (0) J (Z) [K21] [K22] 0 e-[A]z [K21] [K22] Ji (0) (2.51) where [A] - diag [A,..., A,] (2.52) Rei(A) > 0 Oor (2.53) Re(A ) = 0 Im(A) > 0 (2.54) and JZ,i = Z,<,i (2.55) We remark that ( 2.51) is a much more simplified expression and its form is crucial in making reliable numerical calculations. In the following analysis it turns out that only e-[A]l 1 > 0 appears in the calculations. Therefore, since the real parts of the elements of [A] are chosen to be nonnegative, no numerical instabilities occur. Using ( 2.51) and following the analysis in Appendix A, the scattering matrix of the array of curved grooves at its interface is derived to be [S] = ([W] + [I])-1 ([W]- I]) (2.56) where J [W] -= [L]NT [FK] [G,]-1 [L,] (2.57) n=1 and [I] is the unit matrix. Also

31 [F,] [K1l] + [K12] [Ra] (2.58) [G,] - [K21] + [K22] [R] (2.59) [R] - e-[]2 [K12]-1 [Sn,end] [Kll] e-[A]2 (2.60) [Si,end] is the scattering matrix at the back end (z = -12) of the ith groove, /2 is the length of the compressor and [L,] is identical to the [L,] matrix defined in section 2.3. 2.4.2 Numerical Results Since there were no reference solutions when the method was developed, for validation purposes we initially relied on the examination of the boundary conditions satisfied by the transverse electric fields. Fig. 2.9 depicts the p component (a) and the ( component (b) of the total electric field at both sides of the interface. The dimensions for this calculation are: b = 3.5A, a == 2.5A, 11 = A, 12 = A, 4 blades, = — 45~, 01 = ~,F'(z) = 30~/A. The field is plotted at p = 3A, as a function of 4, for an incidence angle of 0i = 0~ (a) and 0- = 30~ (b). The plots correspond to ~3 polarization. As in the case of the straight blades, the incident field is also plotted for reference purposes. It is evident that the fields on each side of the interface are almost equal and the tangential fields on all metallic surfaces are very low. Consequently, the boundary conditions are satisfied, implying that pertinent RCS plots are correct within the limits f the Kirchhoff approximation employed at the aperture. Indeed, our RCS results proved to be in very good agreement with actual measurements [29]. Unfortunately, the terms of the contract do not allow us to publish the comparisons. It must be pointed out that in the case of curved blades, as opposed to simpler geometries, a rigorous convergence check is a very hard task, because the condition

32 number of the involved matrices and the individual eigenvalues tends to deteriorate severely as the number of modes increases.

crq 0 0o0 CD 0C Ci) Normalized field (dB) Normalized field (dB) C 5- 8 v tA I 0 e 5 w (. b 9 9 9 S 9 I A -0 -4.4 9-% 0 -teaV,%go 9 "I. aW.'-"pw-'.,. - fl? A... 1 ts..4 i LA r- -f cn o t14 9 p Uti r --- — -, -- ----- - -r. - - --- ^ - W ew -^.,.. ^ t. -. I.t.-. -$..l.. OA! I I. I". I a: 1 -.9 T mq

CHAPTER III THE ADAPTIVE INTEGRAL METHOD (AIM) 3.1 Introduction In Chapter II we applied the Mode Matching (MM) method to the problem of electromagnetic scattering from jet engine inlets with three different terminations. The method is rigorous, accurate and powerful, since it yields excellent results for even fairly large structures (we showed data for inlets up to 12A in diameter). Nevertheless, the method faces serious limitations. Its most significant disadvantage is that it can handle only canonical geometries, where modes (eigenfunctions) can be defined. In reality, the termination of a jet engine is geometrically very complex, and consequently the Helmoltz equation cannot be solved analytically in that region. Moreover, even for canonical terminations, analysis of a realistic problem, which is about 50A wide in diameter, yields extremely large, full linear systems, almost intractable by modern computing facilities. In addition to the MM analysis, we also investigated other efficient numerical algorithms applicable to the same problem. A Finite Element Method (FEM) code was written and validated for a number of engine-like geometries [11]. The FEM method proved very versatile since it can handle arbitrary terminations. However, it also suffered from storage and CPU time restrictions that may become severe for 34

35 inlet diameters larger that 2A. Our need for flexible, numerical methods, applicable to very large, irregular geometries led us to the investigation of fast integral equation techniques. Among these, the Adaptive Integral Method (AIM) [17, 18] was given special attention due to its promising properties. The main idea of the method is to replace the full Moment Method (MoM) impedance matrix with a sum of special format matrices aimed at reducing the MoM computational complexity and memory requirements. This improvement can be achieved by approximating the interactions among far zone elementary currents using delta current sources located at the nodes of a rectangular grid. It can be shown that the original impedance matrix reduces to a sum of two terms: the first is a sparse matrix and the second is a sum of products of sparse and Toeplitz matrices. In this chapter the mathematical development of AIM is presented. Results are also given using a computer code based on this formulation. 3.2 Mathematical background Consider an arbitrary, perfectly conducting (PEC) surface S illuminated by some incident field EP. The scattered field E' is given by Es = — jwA - (3.1) where A is the magnetic potential defined as A (r) = 4 ( Js (r ) R' d 2S (3.2) and the scalar potential q is given by

36 = s (r') e -jkR}d2S' (3.3) In these equations R - Ir - r', Js and Ps denote the surface current and charge densities, respectively. The continuity equation Vs J, -jiwps (3.4) provides the relation between the two densities. For a numerical solution of the pertinent integral equation we discretize the surface into small triangular patches and expand the unknown current J, using a suitable set of basis functions fn(r). We let N J5(r) _ S Inf (r) (3.5) n=l where In are unknown coefficients (elementary currents). A popular choice for fn(r) are the Rao-Wilton-Glisson (RWG) basis functions defined by n+ P+ ifr T+ IAn - fn) = n if r T- (3.6) 0 otherwise where In is the length of the nth edge and A' is the area of triangle T$ (see Fig. 3.1). It is crucial to point out that fn(r) is tangential to the nth patch and vanishes outside it. On our way to construct an integral equation for the solution of the elementary currents, we enforce the boundary condition (El + Es). t = 0 on S (3.7)'

37 In Tn 0 Figure 3.1: The Rao-Wilton-Glisson (RWG) basis functions

38 where t is any tangential vector on the surface, or, approximately, tangential to the triangular patches. Finally, to test the integral equation per Galerkin we write, using ( 3.1) and ( 3.7) E2. fmd2S j jwA fmd2S + J V fmd2S or, by invoking an appropriate vector identity (3.8) J Ei fmd2S = jwA fmd2S — JV, V fmd2S JY d JwA Jy,. vs- m m = 1,..., N (3.9) which is equivalent, after use of the expansion ( 3.5), to the matrix equation [Z]{I} = {V} (3.10) The elements of the impedance matrix [Z] are given by Zmn = 8''n, exp I — jk — - [rJ I A exp{-jAti} p*d2Sd2Sl - T mMn exPP kR} d2Sd2$ (3.11) where +1 Cm = - I-1 if r T+ ifr C Tm (3.12) and I+1 Cn = - -1 if r' E T+ if r' C T (3. 13)

39 whereas?]= //c and Vm= EI fmd2S (3.14) Since the impedance matrix in ( 3.10) is fully populated, serious memory limitations arise when the geometry is electrically large. Also, the required numerical integration to evaluate the [Z] matrix entries and its direct inversion are both extremely time consuming tasks. These difficulties can be overcome by utilizing the AIM technique. A basic AIMI concept is to partition all pairs of matrix entries into "far field" and "near field" ones. A far field threshold rthr is set according to the desired accuracy. If the centers of two interacting current elements mth and nth lie at a distance larger than rthr, they are considered to be in the far field region of each other, (see Fig. 3.2), and their interaction can be modeled in an alternative way. To describe the far field interactions we envision the entire scatterer being submerged into a rectangular grid (Fig. 3.3). We introduce an auxiliary set of basis functions, namely 4bm, represented by clusters of delta current sources located at the nodes of the grid, namely A13 m(r) = E ( - mq)S(y - Ymq)(Z - Zmq) [A x + A qy + Azq] (3.15) q=1 where M is equal to the expansion order, and rmq are points on the grid surrounding the mth edge. In Figs. 3.3, 3.4 the relative position of the nth RWG patch in a M = 2 order grid is depicted. The A coefficients are determined to ensure that the two sets of basis functions 4'm and fm are equivalent. However, since 7Im are defined in the whole volume, whereas

40, ' -" I,.-.I.:;;4 '' l ->Ji ' ~?. t''' I;< I r \: -~ s..-., I. I,,~ A, s. i i "I". t:"..'.' I:~<::;~....!UH IY. -U *A< " I '. rthri i. 'y Figure 3.2: Definition of the far field region: fn lies in the near zone, whereas fp in the far zone of fm.

41 Figure 3.3: Uniform AIM grid surrounding the scatterer.

42!o ------------------------- ------ II I I I I I I I * I I I I * I I I *I I I f n I I I I I I I I I a I+\ I I I I +. I I I I I I I \ / I I * I \I I I I I-...X.... I I I I I I I I I I I I I I/ \ / \ I I / I I / I I I II I I r 0 Figure 3.4: Grid surrounding a Rao-Wilton-Glisson basis function

43 each fm is nonzero only on the mth patch, formal equivalence can be achieved only after a slight modification of the RWG basis function definition. Specifically, let fm(r) -fmS(z)r (3.16) where 6(z~) is a ID delta function and z~ is the local coordinate along the normal to the + or - triangles which comprise the RWG diehedral mth patch (see Fig. 3.5). This is an important step to ensure the mathematical rigor of AIM, and was omitted from the original presentation given by E. Bleszynski et al. [17, 18]. The two sets 4mh and fm can become equivalent, by imposing equality of their moments up to order M -X oc with respect to the midpoint ra of the mth edge. These moments are defined by 0t *OO 00 Mq 3 / / (r ) - x,)'(y - ya)2(Z- za,)3dxdydz = -oo -oo J -oo M3 (= (x - T)l (yq - Ya)m2(Zmq - Za)3 [AnqX + Ay q + Aq] (3.17) q=1 00oo oo roo - Mlq = fm(r - )( xa)ql(y a)q2( - Z a)3dxdydz (3.18) -00 -00 -00 Apart from fm itself, the surface divergence of fm must be similarly approximated. We define M3 X(r) = E (x - Xmq)(y - Yq)S(z - Zmq)Amq (3.19) q=1 The pertinent moments are Doo roo r00 D = (. J- - a(r)(x -a)ql(y-ya)q2(z- Za 3dxdydz -- -00 00 J -00

44 Vertex 2~ Vertex 3~ y z z Vertex 1+ Vertex 1 - Figure 3.5: The local coordinate systems on triangles T~

45 M3 = (mq - Xa) 1(Ymq - Ya) 2 (Zrq - Za )q3Adq (3.20) q=l DqM23 = / J Vs. ( - xa()l()) y- ya)l 2( - za)q3 dxdydz (3.21) -OC 00 -00 By equating Mq M to Mm and Dm to Dm we obtain four M3 x M3 q l q2 q3 ql q2 q3 ql q2 q3 qlq2 q3 systems that yield the A coefficients as a solution. A considerable portion of our effort was focused on the efficient calculation of the moments in ( 3.18) and ( 3.21). A very efficient, closed form scheme was developed, the mathematical details of which are described in Appendix B. Since the problem is now formulated in the whole volume, and not on surface Z, integral equation ( 3.8) must be modified. At any point in the space an equation of the form Ei(r) + E (r) = h(r) (3.22) is valid, where h(r) is an unknown function with the following property: h(r) t = 0 on S (3.23) t being any tangential vector on E. After the surface discretization, ( 3.23) holds approximately for the assembly of the triangular patches. Furthermore, ( 3.22) yields v [E(r) + Es(r)] f (r)d3v = h(r) fm(r)d3v = J h(r). f(r)8(zm~)(r)d3 = 0 (3.24) The equality to zero is due to ( 3.23), the definition of the delta function and the

46 fact that fm is tangential to the mth diehedral. Therefore, due to ( 3.24) the testing scheme reads Ei fmd3v = jiwA fmd3v + V& fmd3v (3.25) together with ( 3.5). For M - oo the sets mb, and fm are equivalent, hence the testing in ( 3.25) may be performed by either of them. For practical applications, though, M is finite and only interactions among elements that lie in the far field of each other can be modeled through the set bm. However, since the whole geometry is submerged in the grid, it is inefficient to separate near field and far field interactions a priori. It is preferable to calculate all possible interactions using the auxiliary basis functions and afterwards replace the near field interactions by their exact values. In Appendix 13 it is shown that the approximate impedance matrix for the whole geometry calculated via the auxiliary basis functions is given by 4 AIM = [(i)][G][L(i)] (3.26) i=1 where [L(i)] are sparse matrices, containing scaled versions of the equivalence coefficients A in ( 3.15), ( 3.19) and [G] is a matrix with Toeplitz properties. Since the [L(i)] matrices map the original RWG mesh onto the rectangular grid, we will refer to them as "mapping matrices", whereas [G] will be referred to as the "Green's function matrix". Explicit expressions for the elements of these matrices are given in Appendix B. Obviously, [Z]AIM can be split into two parts, corresponding to the near and far field interactions, namely [z]rtotal = rnear rl[Z f(ar LAIM [ ZAIM + [L ]AIM (3.27)

47 Furthermore, the exact impedance matrix in (3.10) is likewise split as [z] = [z]A + [z]far (3.28) Since for the far field interactions [Z]far c [Z|faj it follows that the impedance matrix in (3.10) can be approximated by [Z] - [Z] -[Z] + [Z] = [S] + [L()][G][L()] (3.29) i=1 where [S] is a sparse matrix corresponding to the difference between the exact and the AIM-modeled near field interactions. Of most importance in the formulation of the impedance matrix is that [G] is Toeplitz and thus the FFT algorithm can be used to significantly accelerate the calculation of matrix-vector products required for iterative solution of the linear system. Details on the utilization of the Toeplitz. property are given in Appendix B. 3.3 Application to Relatively Flat Surfaces Computational cost estimates show that the memory requirements and complexity of the algorithm depend both on the original number of unknowns N, and the total number of grid points N9. For large N, our double precision implementation has a memory requirement of Mem - (368 + 32M3)N + 314Ng 4- 16Nnear bytes (3.30) and the number of complex multiplications per iteration, assuming a symmetric Biconjugate Gradient (BiCG) algorithm and using a radix-2 FFT is Nmuit " 8M3N + 540N9 + 120N9 log2 N9 + Nnear (3.31)

48 where Nnear is the number of non-zero entries in [S], and depends on the geometry of the scatterer. For rectangular surfaces it has been shown [18] that Ng is asymptotically proportional to N3/2. However, Ng is in general highly dependent on the geometry. To observe this, let us consider a scatterer of rectangular form whose sides are denoted as a, b, c. If the grid step is chosen as h, then abc g - h3 (3.32) Assume that the surface of the scatterer is discretized by equilateral triangles of edge length 1. The area of each triangle is therefore ST = V12/4 and the number of triangles over the whole surface is NT 32 (ab + be + ca) (3.33) 312 Since the number N of RWG elements is approximately N z 3NT/2 it follows that," abc 3 Ng N3/2 (3.34) [43 (ab+bc+ca)] ' 3/2 A similar estimate is possible for Nnear. Assuming a locally smooth surface, the circular disk around a particular RWGC element lying within the near field of the element has area equal to Trr2h. Assuming triangles with edge length I the number of RWG elements in the disk is equal to NRWG 2V2 (rj) (3.35) and therefore Nnear = NRWGN ( N (3.36) 2 K i

49 The factor 1/2 is due to the symmetry of the matrix. It is pointed out that these estimates are rough and not very accurate for general geometries. They only serve as a qualitative measure for comparisons between AIM and other relevant techniques.. The above analysis, though approximate, shows that AIM applicable to surface has O(N3/2) memory requirements, and O(N3/21og2N3!2) complexity. Therefore, for a large number of unknowns N it always performs more favorably than the MoM (with an iterative solver), which is a O(N2) method. To demonstrate the dependence of VNg on the relative dimensions of the rectangular scatterer, let c = ya = yb, where -y is defined as the "flatness" parameter of the rectangle (see Fig. 3.6). As shown in Fig. 3.6 for small 7 the rectangle is relatively flat and becomes more like a plate as 7 - 0. The constant in front of N3/2 in ( 3.34) becomes smaller as 7y decreases and thus it is important to examine the performance of AIM for flat geometries. Performance comparisons between AIM and iterative MoM for various flatness parameters, h = I and rthr = 101 are given in Figs. 3.7 and 3.8. For both MoM and AIM the BiCG solver was assumed and it is clear that AIM is much more efficient than the MoM for large N. As expected from ( 3.34) its efficiency increases with the flatness of the scatterer, and for the chosen parameters AIM is preferrable to MoM in terms of memory and complexity for about N > 1500. As N becomes larger the improvement of AIM over MoM cannot be overemphasized. Moreover, for relatively flat surfaces, the AIM grid is also flat, permitting the use of much faster, nearly two-dimensional FFT's, resulting in a dramatic reduction of the computational cost.

50 7=1 C -- -----------, -b a A l0.5 -L — -- ---- ~ -- a 70.1 C b - a- - - a ~ ---b Figure 3.6: Scatterer shape for various values of the flatness parameter y.

51 Memory requirements for M=3 0 >. 0 (U) 1200. 1000. 800. 600. 400. 200. 0. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100. N x 1000 Figure 3.7: Memory required for standard MoM and AIM (surface problems with expansion order M = 3, grid step h = I and rthr = 101).

52 CPU Al.h time per 100 iterations (C90 at 275 Mflops).............................................................,T+-J. 360. 300. 240. 180. I""''. VE E +-j u I............................... - MoM AIM with y / AIM with =O.1 - - AIM with y=0.01 7, t,t. 7..1T T 1 I. I, I., /'-:.; -;L'.:,:,r,'7,7,3,,,,.,,-,,,.:,t,:,,:;I.,-,::.,1.,,',.I,,,,, 120. 60. 0...t~l.. B-. - -t- v.. <.......................................................................................................................................................................................... 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 100. N x 1000 Figure 3.8: CPU time for standard MoM and AIM (surface problems with expansion order M = 3, grid step h = I and rthr = 101).

53 3.4 Numerical Results and Parallelization Issues Based on the analytical results described in this chapter, an AIM-based code was developed and validated. An already existing 3D MoM code [32] was modified so that the far field interactions among the current elements were modeled according to AIM principles. The scatterer was a flat PEC plate with dimensions 4A by 0.25A (see Fig. 3.9), corresponding to an original 317 by 317 MoM system. The grid spacing was h = 0.05A, and the total number of grid points was 4500. The edge length of each elementary triangle was I = 0.1A and the expansion order was M = 3. Numerical results for the Radar Cross Section (RCS) were obtained and compared to standard MoM data, showing excellent agreement (see Figs. 3.10 and 3.11). The near field threshold rthr was set equal to A, but it could be decreased down to 0.3A without significant deterioration of the RCS pattern. To further test the code for more complex geometries, we investigated a simplified airplane-like structure (see Fig. 3.12). The RCS results are given in Figs. 3.13 and 3.14, and again they are in excellent agreement with the reference MoM. The original system was 690 by 690, the grid spacing was h = 0.05A, the expansion order was M = 3 and the near field threshold rthr was set equal to 0.8A. Profiling the AIM code showed that a very large portion (up to 75%) of the total CPU time was consumed by the FFT. Although the AIM system solution is faster than the MoM, the major bottleneck of the AIM algorithm is filling [S] in ( 3.29), a process which may also utilize the FFT. These observations were not unexpected and point to a need for parallelization. For general geometries the algorithm has been parallelized via domain decomposition [18], but this approach is fairly complicated and difficult to implement. In our case, due to the geometrical characteristics of

54 X Einc / -x Figure 3.9: A 4A by 0.25A flat blade illuminated by a plane wave. the scatterer (i.e. nearly flat), the FFT grid is almost two-dimensional, and also no excessive free space is wasted by the rectangular grid enclosing the scatterer. Therefore, speed-up via parallelization can be easily achieved, by simply parallelizing the FFT subroutines employed at each iteration step of the solver. Our parallelization steps for the FFT are [45]: * The array with dimensions n, nyI nz is distributed into planes, equally divided among all available processors of total number p. * One-dimensional FFT's are performed along each of the two dimensions x and y, and this is done nz/p times for each node concurrently, since each node contains nz /p planes. * The array is transposed, while still residing on all processors. * One-dimensional FFT's are then performed concurrently on all nodes along the third dimension.

55 4X by 0.25X plate (p Pol.) 20. 10. Cl U4 0. -10. -20. -30. 0. 30. 60. 90. 0 (deg) Figure 3.10: Backscatter for the 4A by 0.25A flat blade of Fig. 3.9 (4 polarization, =- 0 cut). The blade lies on the xy plane with the long dimension along the x axis. The origin lies at the center of the blade. Grid step h = 0.05A.

56 4X by 0.25X plate (0 Pol.) 20. 10. Ct C --- 7d 0. -10. -20. -30. 0. 30. 60. 90. 0 (deg) Figure 3.11: Backscatter for the 4A by 0.25A flat blade of Fig. 3.9 (0 polarization, - = 0 cut). The blade lies on the xy plane with the long dimension along the x axis. The origin lies at the center of the blade. Grid step h = 0.05A.

57 * At this point the transpose of the Fourier transform is available, hence an additional transposition is required to complete the FFT routine. Communication time consumption for a transposition across several processors is significant, and therefore in an optimized code only one absolutely necessary transposition may be performed. If the FFT output and input data must be in the same format, though, the data must be redistributed after the FFT calculation, resulting in a total of two transpositions. The parallelized version of the code was initially tested on the IBM SP2 for the airplane geometry of Fig. 3.12 amd the speed-up results are shown in Table 3.1; the outcome is evidently very favorable, since the overall computation time decreases dramatically with the number of processors. Moreover, to demonstrate the AIM efficiency for relatively flat structures, the parallelized version was also tested for the jagged plate geometry of Fig. 3.16. The dimensions of the plate were 3A x 3A x 0.2A, and therefore the FFT grid was quite flat. We chose a rectangular grid spacing of 0.05A for all directions, and the moment expansion order as M = 3, whereas the original number of unknowns was N = 3108. The monostatic RCS for this configuration is plotted at a X = 0~ cut, for 0 < 0 < 90~ with a 2~ step, and both polarizations (Figs. 3.17, 3.18). Two sets of AIM results are shown, one where the near field threshold is rth, = 0.8A, and one for rthr = 0.3A. From ( 3.30), ( 3.31) and ( 3.36) both memory and complexity increase with rthr and thus a smaller rthr is desirable provided the accuracy of the solution is maintained. As seen from Figs. 317 3.18 the smaller threshold of 0.3A leads to sufficiently accurate results and can be used to evaluate the performance of the AIM code. Tables 3.2 and 3.3 sum up the CPU time and storage comparisons between AIM and the MoM, when the jagged plate is assumed to have zero thickness. In terms of memory the AIM savings

58 are dramatic, especially for mall thresholds rthe. In terms of CPU time, for the zero thickness plate each matrix vector product of the serial AIM algorithm is comparable to the MoM, but significant speed-up is gained for a larger number of processors. For small, but non-zero plate thickness, though, the advantages of AIM over the MoM become much more impressive. For a similar jagged plate with non-zero thickness, the number of RWG unknowns is approximately doubled, resulting in quadruple memory and CPU time requirements for the MoM. On the other hand, since the number of the AIM grid points Ng remains the same, the AIM CPU time, which is dominated by Ng, is virtually unaffected, and even the serial AIM code becomes four times faster than the MoM. Furthermore, the AIM storage requirements are only slightly increased (see Table 3.4 for the thicker plate), whereas those of the MoM would not allow its application on conventional platforms.

59 # Processors Total time (h) Serial section (h) Parallel section (h) 1 (serial) 10.3 10.3 (100%) 0.0 (0%) 1 17.7 1.2 (7%) 16.5 (93%) 2 10.0 1.2 (12%) 8.8 (88%) 4 6.0 1.2 (20%) 4.8 (80%) 8 4.3 1.3 (30%) 3.0 (69%) 16 3.5 1.3 (37%) 2.2 (63%) Table 3.1: Speedup due to FFT parallelization for the scatterer of Fig. 3.12 ---- 0.5k --- A- 0.32k --- 0.2X 4 0.06k -— k 1.16k z Figure 3.12: A simplified aircraft structure.

60 Phi Polarization (0=90~ cut) 10..............'.............. o........."..... '1 '',-............I....I............. MoM 0. __ / -5. -10............... 1 1..................... I 0. 15. 30. 45. 60. 75. 90. f (deg) Figure 3.13: Backscatter of the structure in Fig. 3.12 (4 polarization, 0 = 90~ cut).

61 Theta Polarization (0=90~ cut) 10......I.......................I......................... MoM --------- AIM 5-. --- —------- -------- -5. -10. '........................ I.............. I I 0. 15. 30. 45. 60. 7', (deg) 5. 90. Figure 3.14: Backscatter of the structure in Fig. 3.12 (0 polarization, 0 = 90~ cut).

62 Figure 3.15: Profile on a CONVEX Exemplar for the scatterer analyzed in the previous figures.

63 3X 1*% 4** 4*11 I %4. z4 4*x 4* AL*......1~.2k 4*1.4 3), -10 Figure 3.16: The jagged plate geometry.

64 # Processors Total time (h) Mat.-vec. prod. time (s) FFT / total time 1 (serial MoM) 22.3 5.42 0% 1 (serial AIM) 29.2 5.36 70% 2 30.4 5.64 69% 4 18.7 3.12 55% 8 13.8 2.04 45% 16 11.2 1.48 29% Table 3.2: CPU time comparisons for the scatterer of Fig. 3.16 with zero thickness. Method Memory (Mbytes) MoM 77.3 AIM (rthr = 0.8A) 28.94 AIM (rthr = 0-.3) 17.19 Table 3.3: Memory comparisons for the scatterer of Fig. 3.16 with zero thickness. Method Memory (Mbytes) MoM 309.2 AIM (rthr = 0.8A,) 74.41 AIM (rthr = 0.3,) 27.41 Table 3.4: Memory comparisons for the scatterer of Fig. 3.16 with non-zero thickness.

65 Jagged plate (Phi Polarization) 30................................................ 30. MoM -.,,,f ------- AIM (r,,r=0.8X) - - - - ------ AIM (rthr=O.3X) CII 10. -' 10. 0. -10. '........................................... 0. 30. 60. 90. Oi (deg) Figure 3.17: Monostatic RCS for the jagged plate geometry (4 polarization).

66 Jagged plate (Theta Polarization) 30.............................................. MoM AIM (rthr=0.8X) 20 --- AIM (rthr=0.3X) inV I,r // \r N u 0. -10.......................................... 0. 30. 60. 90. 0i (deg) Figure 3.18: Monostatic RCS for the jagged plate geometry (0 polarization).

CHAPTER IV EXPLOITATION OF THE CYLINDRICAL PERIODICITY OF THE JET ENGINE GEOMETRY 4.1 Introduction It has already been discussed that for an efficient solution of the realistic jet engine problem, decomposition of the computational domain into a numerical and a high frequency region is the most appropriate approach. As shown in Fig. 4.1, the region surrounding the complex engine termination is modeled numerically, whereas ray or modal methods are used to propagate the field through the duct. In this manner, significant amount of memory is saved by not discretizing the duct region. One way to couple the two computational regions is through the modal scattering matrix. However, even with this decomposition, a direct application of numerical methods such as the FDTD [33] and FEM [11] leads to intractable problems when dealing with realistic jet engine sizes spanning 50 wavelengths in diameter. In Chapter III a novel numerical technique, the Adaptive Integral Method (AIM) was described and shown to be particularly suitable for the analysis of complex blade terminations. However, direct application of AIM to the entire scatterer would not be efficient either. Multispectral characterization requiring radar cross-section (RCS) computation at many 67

68 frequencies place additional CPU time problems. To reduce the CPU time and storage requirements down to manageable levels, the inherent blade periodicity of the jet engine was exploited to show that the computational domain can be reduced down to a single blade or engine "slice". This discrete body of revolution (DBOR) approach reduces the number of unknowns and storage requirements by a factor equal to the number of blades. However, it has so far been implemented only in the context of differential equation methods [12]. In this chapter we apply the DBOR concept in the context of integral equation methods for modeling the complex jet engine configurations. By considering mode by mode excitation [11, 12] it is shown that the analysis over the entire engine can be reduced down to a surface integral equation over a single blade. It is further demonstrated that the resulting modal scattering matrix is sparse, leading to additional storage and CPU time reductions. Further computational reductions can be achieved by incorporating fast integral equation algorithms such as AIM, as discussed in Chapter III. We have seen that the latter redistributes the currents on the blade onto a canonical grid such that the scattered field due to the fictitious sources on the nodes of the grid remain the same as that of the original sources up to a certain order. Combining blade periodicity and the AIM concept results in drastic complexity reductions. The DBOR procedure (which takes advantage of blade periodicity) leads to a compact integral equation whose domain is confined over the single engine blade. Initially we consider perfectly conducting scatterers, but the last section of the chapter gives an extension of the method to periodic angular sectors containing non-metallic materials. Appendix C provides certain mathematical details ommited from the main text, regarding the dyadic Green's function for cylindrical waveguides and the

69 I \ ^-C / I a Engine Modal or Ray region I I I ; allow - - -- -11 I ' __ ---0 — Fields to Ray/Modes Integration Surface (Hand-off Surface) Numerical or Rigorous Modeling I raion of the computational decomposition. Fipuree 4.1'. ll1tration of th comp

70 Figure 4.2: Simplified model of a jet engine. explicit expressions of the elements of the pertinent MoM impedance matrix. 4.2 Preliminary Setup for Perfectly Conducting Scatterers Consider the perfectly conducting (PEC) scatterer residing inside a cylindrical metallic waveguide and occupying the volume V. as shown in Fig. 4.2. This structure is used to terminate a cylindrical waveguide and consists of N, sectors or "slices". The slices are identical and rotated around the z axis with respect to each other by integer multiples of the angular period 27r 2 = (4.1) Ns where 0s is the angular opening of each slice (see Fig. 4.3). We will refer to this termination as a 'Discrete Body of Revolution" (DBOR) since the geometry can be

71 V2 V1 Vo ` - f VNs -1 = 27 p 2s 2ic/Ns - I a a z Figure 4.3: Symmetry in a DBOR body. generated by revolving the basic slice, occupying the volume Vo a discrete multiple of azimuth angles s,. Because of this inherent periodicity, the vector r of a point in the mth periodic slice occupying the volume Vms can be written as r Rm ro ro E Vo (4.2) where the superscript m, denotes the power of the dyadic R. The latter is the rotation dyadic, defined in cartesian coordinates by the matrix

72 cos 4s -sin ) 0 [R]= sin O cos Os 0 (4.3) 0 0 1 It can be readily shown that R-1 RI (4.4) R- (s) = R(ms) (4.5) which are important properties to be exploited later in the analysis. For the problem of scattering, an incident field E2 is assumed to impinge upon the periodic structure, where E2 is a sum of cylindrical waveguide modes. This excitation induces the surface current J giving rise to the scattered field [34] E'(r) = -jkZ j G(r, r') J(r')d2S' (4.6) where k = y, Z = t//cl, E is the outer surface of the scatterer and G is the dyadic Green's function for the Helmholtz equation inside a cylinder (electric type of the first kind) stated in Appendix C. An integral equation for J can be derived by invoking the boundary condition satisfied on E, namely (Es + E) 0 t =- 0 (4.7) where t is the tangential vector on Z, demanding that the tangential electric field vanish on the surface of the periodic scatterer. A numerical solution for J can be obtained by casting ( 4.7) into a discrete system of equations. To (lo so, we first approximate the surface current J by the expansion

73 J(r') = f (4.8) q where fq (r') are the chosen basis functions and Iq are coefficients to be determined. A standard set of linear basis functions are those given by Rao et.al. [35] defined as q + + if r T+ + 2A- q i q (r)= if r T- (4.9) 2A q ifqqE T 0 otherwise where A' denote the area of the two triangles forming the dihedral patch, as illustrated in Fig. 3.1. The basis functions represent the current flowing through the qth edge from one triangle to the other, where p' is measured from the vertices opposite to the qth edge. A useful property of the basis functions is fq (RmS r) = Rm. fq(r) (4.10) and this is a consequence of linearity. By introducing ( 4.8) into ( 4.7) upon Galerkin's testing we obtain the linear system [A]{I} = {V} (4.11) where Vp= j E(r) fp(r)d2 (4.12) are the elements of the excitation vector {V} and {I} is the column of unknown current coefficients Iq. As usual, [A] is the impedance matrix whose entries are given by Apq = jkZ / j fp(r) G(r, r') fq(r')d2SdS' (4.3) J (4.13)

74 Note that ( 4.11) involves the currents over the whole scatterer, covering the entire angular domain 0 < 4 < 27r. Thus, for large diameter terminations, the size of the impedance matrix quickly becomes unmanageable. In the next section it will be shown that the scattering from the slices can be isolated from each other, yielding an equivalent system much smaller than ( 4.11). 4.3 Decoupling of the Linear System of Equations To exploit the periodicity of the geometry, we proceed to establish a relationship between the currents among the identical slices. Assuming that each periodic slice is identically discretized into Q patches, the current expansion ( 4.8) can be rewritten as N-l1 Q J(r)= E Er), rI(mo c Vo (4.14) ms=O q=l This is an expression of the total current in terms of local currents on the individual slices with index m,. To facilitate further manipulation we define the current column vector {J}(ms) of the mth slice by J}()- [I, 1 )..., I (4.15) Using these notations, (4.4), ( 4.5) and ( 4.10) can be used to more explicitly rewrite the system ( 4.11) as

75 [A](0){J}(~) + [A](01){J(:) +... + [A](ON,-1){J}(N -1) = {b}() [A](1O){J}(0) + [A](11)J(1) + [A](N-){J}(N-) = {b}(1) [A] (N-,){J} (0) +.. [A](N —1{J}(N-1) - {b}-(N1-) (4.16) where n,,~m,,) n,, f M, 2(R T Apq jkZ (R( R- ~. (r, R ~ r.fq(r')d Sd2S' (4 17) are the entries of the tsubmatrices [A](nms) representing the interactions between the mh and the nth slice currents. Also, {b}(n,) is the excitation column vector of the nth slice, defined by {b}('~) - [V(T'), V(n,... V(n)] (4.18) where V ) = E (R. r) R'. fp (r) d2S (4.19) In the latter equation, So is the outer surface of the slice within the volume Vo. Clearly, the index m, indicates the slice where the source point is located, whereas nr is the index representing the slice containing the observation point. Both indices run from 0 to N, - 1. To take advantage of the blade periodicity we proceed to establish a relationship among submatrices [A](nsms). The principal part of the entries of [A](nsms) (see Appendix C) can be written as

76 oo ( o A(nsm) = jkZ e-n(m-n,)s m=l n=-oo cE fprMnm(r; ~TE TE 2 {CT7 j fp(r) Mnm(r; ~3)Mnm(r'; 13) fq(r')d Sd2S'+.So nno CM j fp(r) Nnm(r; ~;m )Nnm(r'; T/I3m ) fq(r')dSd2S'} (4.20) from which we readily deduce that A(q7ns) depends only on the difference m,-n~, and not on the individual values of m, and ns. This is expected, since physical intuition dictates that any interaction between two given slices should depend on their relative, and not absolute location. Hence, we introduce the superscript, = ms- nr to rewrite the system ( 4.16) as [A](~){J}() + [A](1){J}() +... + [A](N-1){J}(N-1) = {b}(~) [A](J-1)J} + [A](0){J}(1) +... + [A](Ns-2){J}(N-l) = {b}(1) [A]( —N+I)({J}() +... + [A](~){J}(N-1) = {b}(N-1) (4.21) The lone superscript of the impedance submatrices now indicates the relative location of the interacting slices. Since the waveguide fields are identical modulo 27r, it follows that [A]) (N -'-,) (4.22) This "modulo Ns"' property will prove extremely important because of its essential role in decoupling the subsystems comprising the overall system ( 4.21). To proceed further, we assume that the incident field is a single cylindrical waveguide mode of order ni, as described in Appendix C. To permit decoupling of the

77 subsystems comprising ( 4.21), it is important that the angular dependence be in the form of an exponential, and not trigonometric, i.e. the fields must be proportional to expf{jn4n}. The excitation column vectors on the Kth and (K + l)th slices are hence related via {b}-('+) = eJ'Is{b}() > = - e n's{b}(0) K = 0,..., Ns - 1 (4.231) We will refer to the block rows of ( 4.21) by the superscript of the right hand side, i.e. we will use the terms '0th block row", "1st block row, (N,- )t block row". Multiplying now the Kth block row of ( 4.21) by e-i"'i s, for all K = 1,..., Ns - 1, and subtracting the 0th block row from all others yields the equivalent system [A](~){J}(~) + [A](1){J}(1) +... + [A](N-l){J}(Ns-) = {b}(~) [A](O) ({J}(1)e-'Jni - {J}(o)) + + [A](1) ({J}(2)e-ini - {J}(1)) +... + + [A](N-l ) ({J}(O)e-jni {J}(Ns-l)) = {0} [A](0) ({j}(N-l)e-(N-1)ni -- {J}(o)) + + [A](1) ({J}(O)eJ(Ns-l)nis - {J}(1)) +... + + [A](Ne-1) ({J}(N-2)e-i(Ar-l)s - {J}(Ns-)) - {0} (4.24) where we have made use of ( 4.22). Apart from the top one, we can satisfy all block rows of ( 4.24) by setting

78 {J}(1) {J}(2) - ein' {J}(O) = ej2ni, j}(O) {J}(N-1) - ei(N,-l)ni s {j}(o) (4.25) and thus the top block row of ( 4.24) can be written as [[A](O) + [A](1)eJni"s +... + [A](N —1)ej(N-1)n'" ] {J}() = {b}(~) (4.26) Since ( 4.21) and ( 4.24) have the same solution, on the basis of uniqueness ( 4.25) is the only possible choice. A more compact expression for ( 4.26) can be obtained by using the properties of the Dyadic Green's function. First we define the left hand side of ( 4.26) as a single matrix [K] [K] [ A]) + [A](l)ej + +... 4 [A] (NS-1 )e(N'-1)"i'] (4.27) Next, by using ( 4.20) and the geometric sum N_-1 _ N, AN, (n,n,-) E el(ni-n)}S = S =0 O if n = ni + vN,, v E Z else (4.28) it follows that the principal part of the matrix entries Kpq reduces to

79 oo oo KIp = jkZ A A ANv (n,n) - m=1 n=-oo TE M_ - {C { j X fp(r). Mnm(r; ~/3m)M._nm( T/m) f(r')dSd2S + C 'M /j fp(r) Nnm(r;+ TM )N-nm(r'; ) f(r)dSd2S} (4.29) The implication of this result is that for a given order of the incident mode, only a limited set of scattered modes is excited. Namely, only the scattered modes with orders n that satisfy n = n, + vNs, v C Z are reflected back, and this is in agreement with the result given by the FEM analysis of a similar problem [12]. For a body of revolution (BOR), Ns - oo and in this case ( 4.28) implies that n = n, only, which is consistent with classical BOR theory [36]. The most important consequence of this analysis is that, for a given incident mode, it suffices to solve the integral equation over only one slice, using a modified version of the dyadic Green's function. Indeed, if we define this modified Green's function by the periodicity dyadic rn 5(r,- r')z' + 00 00 +: AN, (n, ni) CTMnm(r; ~T)Mn(r'; i E) m=l n=-oo oo oo + E E AN (n,n) Cm Nm(r;,m )N-nm(r'; qT. ) (4.30) rr=l nI=-oo the entries of [K] are given in compact form as Kpq =jkZ j fp(r) rn,(rr') fq(r')d2Sd2S' (4.31) The currents on the reference slice So are simply the solutions to the system

80 [K]{J}() = {bt)}(0) (4.32) and those on the other slices are obtained from ( 4.25). The importance of ( 4.32) cannot be exaggerated. Since the problem essentially reduces to modeling only one slice of the scatterer, the number of equations or unknowns is reduced by a factor of Ns. For a typical Ns = 40, this implies a CPU time and memory reduction each by a factor of 1600. For large scatterers with large periodicity numbers Ns, such as jet engines, the problem can thus be scaled to a tractable size. Moreover, the limited coupling among the scattered modes results in sparse modal scattering matrices which are much easier to store and handle. Also, calculations involving the Green's function in ( 4.30) can be performed more efficiently, since a number of terms corresponding to significant modes of low order are now absent. Finally, it is important to observe that the above formulation does not demand any restriction on the shape of the periodic slice. Therefore the shape of the exterior surface of the slices is irrelevant to the integral equation method and the technique is applicable to arbitrary DBOR's such as realistic jet engines. 4.4 Extension to Dielectric Scatterers It is straightforward to extend the above analysis to a DBOR which consists of non-metallic sections. To construct the integral equation we begin by introducing the scattered field expression ES(r) = -jw / G(r, r') J,(r')d3v' (4.33)

81 where J, now represents the equivalent volume current density replacing the dielectric region with permittivity ed(r), given by Jv(r') = j [d(r') - e] E(r') (4.34) The current density is again expanded into volume basis functions 1bq via Jv (r') Iq lb(r') (4. 35 ) q A choice for Pq is [37] 3 b q (r) = Pq,(r)u (4.36) K=1 where u, 1, 2, 3 denotes unit vectors spanning 7Z3 and P 1 if r E AVq Pq~ (r) = (4.37) 0 otherwise Using the volume equivalence principle it can be shown that the integral equation for the current can be written as [37] Ei(r) = jWu G (r, r') Jv(r')d3v' (4.38) where ~d r S'r - r') Z' 8(r - r')I (r r') ( Cd (r) W) [Ced (r') -] + 00 00 + Z C Mnm (r; +/?T)M-nm(r'; nm m=1 n=-oo + CnmNnm(r; ~/3m )N-nm(r';:TnM) (4.39) m=1 n=-oo

82 in which I is the unit dyadic and Mn (r), Nnm(r) are defined in Appendix C. On the basis of a DBOR scatterer, Cd(R m r)= Ed(r) Vm, C Z (4.40) and by invoking a similar procedure as in the case of metallic scatterers, we can express the entries of the impedance submatrix [A](nsms) as A(n Sin Jvo 1n [r p (r). z] [q (r). ] d3 - p- msnv -r (r) 'q (r) d3v + jWI e-jn( s )s 0 dr) m=1 n=-oo {CnmJ L J (r). Mnm(r; ~/3m)Mnm(r'; T/nm ) fq(r)d vd3v + nm nm nm )q(r')d vdv (4.41) where &mn is the Kronecker delta. Also, the entries of the excitation column vector are given by (cf. ( 4.19)) V (n) - ] R (r) d3v (4.42) Apart from slight modification, the expressions for the impedance matrix and the excitation vector are nearly identical to those for the perfectly conducting DBOR. For reference, the periodicity dyadic in ( 4.30) for the dielectric DBOR is

E- H- H; Q -H I -' .4L) -H -H & III < < 22 II -I- -I

CHAPTER V APPLICATION OF AIM TO THE JET ENGINE PROBLEM 5.1 Introduction In Chapter III we described the Adaptive Integral Method (AIM) and its application to scattering from arbitrary, perfectly conducting objects in free space. It was shown that the technique is particularly suited for electrically large structures since its memory and CPU time requirements are much lower than the Moment Method (MoM) as the problem size increases. Also, it was shown that the storage requirements and complexity are reduced further for relatively flat surfaces, since the FFT in such cases is virtually two-dimensional. In Chapter IV we also showed how the cylindrical periodicity of the jet engine can be exploited to reduce the computational domain of an integral equation method, such as AIM, down to a single periodic sector of the structure. We developed a "slicing scheme" which dramatically reduced memory requirements and the computational complexity of any integral equation method. The main objective of our analysis is to combine the concepts in Chapters III and IV in applying AIM to the jet engine inlet structure. We have already demonstrated that the most appropriate way to analyze this problem is to decompose the engine 84

85 into two computational domains: the hollow duct, where a modal or high frequency method is applicable, and the termination structure. For the latter, a numerical technique is necessary due to its complex geometrical features. The fields in the two computational domains are matched at the interface through a modal scattering matrix. In this chapter we apply AIM to the field computation inside the termination of a jet engine inlet. Since the target lies in the interior of a cylindrical waveguide, and not in free space, several AIM concepts must be modified accordingly. The cylindrical periodicity of the geometry (see Chapter IV) will also be exploited. Finally, several numerical results will be presented, to demonstrate the validity of our method. 5.2 The Cylindrical Waveguide Green's Function In Chapter IV we used the properties of the cylindrical waveguide modes and the pertinent dyadic Green's function to show that an integral equation modeling of the cylindrically periodic termination reduces to the analysis of only one periodic "slice". Although the analysis is theoretically sound, the dyadic Green's function for this geometry is difficult to treat computationally. The most serious difficulty of the dyadic Green's function [34] is that the double modal series (see Appendix C) converges very slowly if the source and observation points lie close to each other, and it actually diverges if they coincide. A similar situation occurs for the free space Green's function, which has a pole at r = r'. In the case of free space it is easy to show that the singularity is integrable. A standard technique which circumvents numerical problems is the extraction of the singularity from the original kernel, and its separate analytical integration. In the case of the waveguide Green's function, though, it is not clear as to how the singularity can be

86 extracted. For the cylindrical waveguide the analysis is especially cumbersome due to the presence of Bessel functions, their derivatives, their zeroes and the zeroes of their derivatives, the latter two not being given analytically. This particular difficulty can be avoided by transforming the summation over the zeroes to an infinite integral [38].. It can then be shown that the integral diverges again if r = r'. A temporary remedy to this situation is to subtract the asymptotic expression of the integrand, and then add it back after integrating it analytically. The integration can be carried out by using the Poisson's formula and other relevant lemmas. Unfortunately it turns out that, even so, the second series (over the orders of the Bessel functions) diverges. A second drawback of the standard expression of the dyadic Green's function, is the odd behavior of the delta function term along the zzc direction (see Appendix C) for surface problems. If the scatterer is perfectly conducting, the Green's function is integrated over the surface, and not the volume of the target. However, the delta function is three-dimensional, meaning that its one-dimensional component along the normal to the surface is not intted, and hence the result of the surface integration is equal to infinity on the scatterer. Apparently, this infinity can be cancelled by the modal double sum which also diverges in this case. However, the mathematical details of this cancellation are not known. Since use of the dyadic Green's function for the cylindrical waveguide leads to serious difficulties, the most efficient way to perform numerical calculations in the termination region is to instead use the free space Green's function, subject to sev — eral modifications of the initial geometry. Indeed, although the excitation is still a cylindrical waveguide mode, the free space Green's function does not account for the boundary conditions on the waveguide walls. Thus auxiliary electric currents are also placed on the walls, computed by enforcing the vanishing tangential field boundary

87 condition on these walls. This is essentially equivalent to considering the wall as part of the scatterer. Of course, meshing the whole wall surface is not efficient, since this results in an unnecessary increase of the number of unknowns. Domain decomposition is a much more attractive approach (see Fig. 5.1). The termination and a small portion of the walls beyond the hand-off plane is discretized for analysis via the integral equation method. The incident plane wave is coupled to waveguide modes in the hollow duct, as in Chapter II. The modal field at the hand-off plane (see Fig. 5.1) is the actual excitation of the termination analysis. After numerical solution the scattered field is calculated at the hand-off plane and the modal coefficients are extracted by projecting the fields at the hand-off plane onto the waveguide modes [11]. The rest of the procedure is identical to the one described in Chapter II. Although use of the free space Green's function inside the cylinder is not expected to yield exact results, it has been shown [39] that the numerical errors are not significant, especially for large inlet diameters. Use of the free space Green's function implies that the results of Chapter IV have to be modified accordingly. Specifically, although ( 4.26) is still valid, the sum cannot be given in closed form, meaning that submatrices [A](k) must be calculated explicitly. This procedure is evidently time consuming, but is still NS faster than the regular MoM (without invoking the slicing scheme), whereas memory is still reduced by a factor of N2. The results of Chapter IV concerning the coupling patterns of the modes are still taken into account, retaining the sparsity of the modal scattering matrix.

88:S -1...... Modal or Ray region I I I^^ I Engine I -r. Fields to Ray/Modes Integration Surface (Hand-off Surface) Numerical or Rigorous Modeling Figure 5.1: Decomposition of the computational domain

89 5.3 Application of the Slicing Scheme to AIM Although the application of the slicing scheme greatly enhances the efficiency of MoM, the benefits for AIM are not that impressive without substantial modifications of the method for its application to this structure. The main reason is that the Toeplitz property of [G] in ( 3.26) can only be retained by extending the rectangular grid over the whole scatterer, and not over a single slice. This can be seen by examining ( 4.26). Suppose we need to calculate the second term, i.e. the interaction between slices 0 and 1. Let slice 0 be embedded into the local rectangular grid, as shown in Fig. 5.2. Also, let the mapping from the original basis functions on slice o to the uniform grid be performed via matrices [As], [Ag], [Az], [Ag] as explained in Chapter III. To construct a local grid for slice 1 we cannot simply rotate the respective grid of slice 0 without destroying the Toeplitz property of [G]. The only possible solution is to extend the local grid of slice 0 so that it also encapsulates slice 1. By doing so, the mapping matrices [Af], [AY], [Af], [Ad] for slice 1 are, in general, different from those of slice 0. Therefore, the interaction between slices 0 and 1, i.e. the second term in ( 4.26), according to ( 3.26) is equal to ejtnSs Z[L(i)][Go1][Lo)]T (5.1) i-=l where [Go,] is the Green's function matrix defined on the grid covering both slices 0 and 1. Similarly, the interaction between slices 0 and 2, i.e. the third term in ( 4.26) is equal to i [L([Go[L) (5.2) i=l

90 e.t.c. To efficiently compute the sum in ( 4.26), it is necessary that the domains of all matrices [Goi] be extended over a grid which encapsulates all slices. Of course, the domains of the mapping matrices [L] must also be extended with the additional entries set to zero. If [G] is the extended matrix, the sum in ( 4.26) is readily given by 4 E[L()] [G][L(()]T (5.3) i=l where [L()] [[L(0)] + [Li)]e/n? +.+ [LNs-]e ) ] (5.4) Although the faormulation is compact, the result is not favorable from a computational point of view. The reduction in memory by employing the slicing scheme is still important, since there is no need to separately store the mapping matrices for any but two slices (slice 0 and a collective "slice" corresponding to ( 5.4)). However, since the extended matrices are defined on the global grid, which encapsulates the entire geometry, the FFT dimensions are quite large and result in serious deceleration of the algorithm. To avoid using a global grid and large dimensions for the FFT, the only possible solution is to construct a grid that still retains the Toeplitz property for the Green's function matrix [G]. Such a grid is not rectangular, but cylindrical, as depicted in Fig. 5.3. The numbering of a grid point at location r,, can be characterized by three integers, namely pi, p2, p3. If hp, hq, hz are the grid steps along the three directions p, X, Z, the cylindrical coordinates of the grid point at rp are (plhhp p2h,p3hz). Similarly, the cylindrical coordinates of the grid point at rq are (q1hp, q2h, q3hZ). It follows

91 ' Slice 2 / Slice 1 \, *,//,./'.I *, * *S * * \ 0 //, v * * * * I /' ' 0 Slice 0 Local rectangular grid of slice 0 Figure 5.2: Local rectangular grid for the basic slice. Wrong axis alignm < / x \ ^ - i - gx lent Correct axis alignment z. 0 Figure 5.3: Cylindrical grid (front view).

92 that the distance Irp - rq between the points is Irp-r l = \/p h2 + qh2 - 2pqgh2 cos [(P2 - q2) ho] + (p3 - q3)2 h (5.5) and the entries of the Green's function matrix [G] are equal to Gpq =exp {-jk-r — rq[} (5.6) 47r rp - rqJ 3 As in Appendix B, a tensor can be defined by 3 Gq q23 = G (5.7) associated with the aforementioned [G] matrix. From this definition it is evident that the elements of the tensor are actually functions of the differences p2 - q2, p3 - q3, but not of pi - ql, as opposed to the rectangular grid case. It follows that tensor products of the form GPl,P2,P3,q,q2,3 (5.8) ql,q2,q3 (5.8) can be written as E E E Gqip2 -q2 p3-q13 q23 (5.9) ql q2 q3 This is a two-dimensional discrete convolution followed by a regular matrix-vector product, rather than a three-dimensional convolution as was the case with rectangular grids. Since the Green's tensor is not Toeplitz along the p direction, the FFT can only be used along the k and z directions, meaning that matrix-vector product calculations

93 are slower than in the case of the rectangular grid. This drawback of the cylindrical grid is balanced by the smaller size of the FFT dimensions, since the computational domain is reduced down to a single slice. To illustrate the mathematical details associated with the cylindrical grid, we focus again on equation ( 4.26). The auxiliary AIM basis functions defined in slice 0 are (cf. ( 3.15), ( 3.19)) M3 (r) = Z (x - ~)&)(y - (o))(z- Z()) r -- 5x mq -x )5 y- mq mq q=1 [Ax q(O) + Ay(~0) + Az (0)] (5.10) M3 Od(~)(r) =- x())S(y — ym~))(z- z())Amq (5.11) /1"M'mq Ymq mq mq q=1 As before, superscript (0) refers to the local cartesian coordinates associated with slice 0 (see Fig. 5.4). According to Appendix B the AIM approximation to the first term in ( 4.26) is [A](0) = jw {[Ax][G](o)[Ax]T + [AY][G](o)[AY]T + [AZ][G](o)[A]T} + 1 +..[A] [G](O) [Ad]T (5.12) where 40 r(0) G(n)- - Irp (513) andpq and r(0o) = h2 + q2h- 2plqlh os [(p2 - q2) h] + (p3 - q3)2 hz2 (5.14)

94 (see above discussion on the definitions of the various parameters). It is pointed out that p, q span the local cylindrical grid over slice 0 only. Similarly, the AIM auxiliary basis functions located in the K, slice are M3 ^S)(r) = x-x())(y - yy ))6(z - Z ) q=l [Ax A('s) + A ^( + AzA(})] (5.15) M3 I s)(r) = 8x ))6(y -- ( )( y )8(z -,))Ad (5.16) q=l The A coefficients are the same as in slice 0, since the moments of the rotated basis function calculated in the rotated cartesian coordinate system are invariant. Using the transformation formulae x(^) = cos(,,)(0) + sin(3,)yr(~) (5.17) y(-) = -sin(s, )kx(~0) 4-cos( &$)Y(~) (5.18) z() = (O) (5.19) for the Ks term in ( 4.26) we obtain [A]) = jw {[Ax][G]( )[A]T cos(ss) - [Ax][G](K)[AY]T sin(is,,) + [AY][G](s)[A ]T sin(ic/) + [AY][G]('K)[Al]T cos(so)S) + [AZ][G](K)[AZ]T} + 1 [A] [G](K )[Ad]T (5.20) JW( where G( e{ ) } (5.21) pq4 (Kqs) 47Frpq

95 and r() = 2q2h + q2h2 - 2pqlh cos [(P2 - q2) h - Kso,] + (P3 - q3)2 h2 pq h' Pil o 2 _ 2 pl q (5.22) Again, it is crucial to point out that p), q span the grid of slice 0 only, although the Ks slice is being considered. The effect of the rotation is taken into account by the term -Ks.s in the cosine argument of ( 5.22). To rewrite ( 4.27) in compact form, we define the following matrices: [rc]- E [G](s)ej"nis cos(isOb) (5.23) Ks=0 N -1 [rn],- E [G]( S)e~nz' sin(,O,+) (5.24) Ks=0 Ns -1 [r]- E [G]( )e'n^ (5.25) we =0 we obtain the AIM approximation to the matrix [K] in By using ( 5.23)-( 5.25) ( 4.27) as [K]AIM = jwlA {[A'] [r, ][Ax] [Ax][F]- [A Y]T + + [AY] [rL.] [A ]T + [Ay] [rF] [Ay]T + + [Az][r~L][Az]T} + I[Ad] [f] [Ad]T n, ]WE n, (5.26) The matrix products in ( 5.26) are calculated via a two-dimensional FFT, as opposed to the three-dimensional FFT of the free space AIM. An important issue concerning the nature of the cylindrical grid is its orientation in space. No points of the grid may lie on the cartesian axes x and y, because in that case, in accordance with the properties of Vandermonde matrices [18], the systems

96 (1) Y (~) (2) --- y \ \ (2) \ \ Slice 2 ( 1 \Slice 1 Slice 0 X Z / Figure 5.4: Rotated cartesian coordinate systems. for obtaining the mapping A coefficients through moment equality become singular. Therefore, care must be taken so that the cylindrical grid is slightly misaligned with respect to the cartesian coordinate system (see Fig. 5.3). For the same reason the origin should be excluded from the grid. The efficiency and generality of the grid is not affected by this choice. Having developed all necessary tools, the method of analysis can be summarized a follows: Given an engine inlet and a plane wave illuminating its open end, we calculate the coefficients of the incoming cylindrical waveguide modes that couple to the plane wave, as described in Chapter II. Using one incident mode at a time as an excitation at the hand-off plane (see Fig. 5.1), the scattering in the termination

97 region is modeled via AIM. A cylindrical grid is used, in conjunction with the slicing scheme described in Chapter IV. The scattered field is calculated at the handoff plane, where the coefficients of the outgoing modes are extracted via numerical integration. Finally, the outgoing modes propagate along the hollow duct and radiate in free space through the aperture (see Chapter II) yielding the RCS of the structure. 5.4 Numerical Results To check the validity of the formulation, a MoM code was initially written for RCS computations of the jet engine inlets. The free space, instead of the cylindrical waveguide Green's function was used for reasons explained earlier in this chapter, and the cylindrical periodicity of the geometry was taken into account, as shown in Chapter IV. The code was initially tested for a cylindrical hub without blades. The dimensions of the geometry were (using the same notation as in Chapter II) b = 2A, a = A, AI = 3A, 2 = 0.5A. Only one quarter of the geometry was modeled (see Fig. 5.5) and the number of unknowns was N = 1866. A total of 176 modes was used (11 Bessel function orders times 4 modes per order, times 4 for all combinations of TM and TE modes) and the scattered field was calculated at a hand-off plane located at a distance of 0.1A from the hub top. The rim contribution was not taken into account. The RCS results are plotted and compared with Mode Matching data in Figs. 5.6 and 5.7 showing very good agreement. Any disrepancies are due to the series truncation of the Mode Matching solution or the use of the free space Green's function rather than that of the cylindrical duct. The formulation was also tested for a similar geometry with four straight blades, placed perpendicular to the inlet axis. Again, only a quarter of the geometry was modeled, by invoking the periodicity of the structure (see Fig. 5.8). The dimensions

98 of the inlet and the hub were the same as in the previous case, and the number of unknowns was N = 1777. The blade was 45~ wide and was centered around the x-axis. RCS results for this configuration are plotted and compared with Mode Matching data in Figs. 5.9 and 5.10 showing again very good agreement. Finally, an AIM code for the inlet problem was developed, based on the aforementioned MoM code and the mathematical analysis described in this chapter. To prove the validity of the AIM analysis within the cylindrical inlet, the RCS pattern of a limited set of modes was investigated for a shorted waveguide having 2A diameter and an overall length equal to A. Only half the cross section was modeled in this case by invoking the slicing scheme combined with AIM. The maximum distance among grid points was 0.05A, the expansion order was M = 3 and the near field threshold was set equal to rthr = 0.7A. For testing purposes, only propagating modes of order 0 were used in the waveguide, and the relevant RCS results are plotted in Figs. 5.11 and 5.12. The rim contribution was not taken into account. The comparison with MoM results is excellent, demonstrating the validity of the cylindrical AIM algorithm. Memory savings of the AIM compared to the MoM are very significant and analogous to the free space case. However, the true merits of the AIM with respect to speed become evident only for a very large number of unknowns, and especially after code parallelization. This is even more important for inlet problems, because the pertinent impedance matrix is usually ill-conditioned, and a large number of iterations is necessary before convergence is achieved. It is therefore absolutely imperative that an optimized, parallel FFT be used. Considering our experience with the free space case, it is speculated that this parallelization/optimization scheme will permit the present formulation to handle jet engine problems with size and geometrical

99 complexity not amenable to any other known method, numerical or analytical.

100 In,. las'^^ 9 —*.i.' ~~~;-~rr-C ~1;~6 k,? ~~'. s.ts.:;- vX -`;t-~.;;r:: - a:,.l~~ -..s8 s-, A. -. {^y^S ~*; m/ I:# *^T-M '-' ~'*^* ' ' - 'a, ''':' * \ ' *-;..l,,, Fi;gur 5.5 A qurefau mtr Figure 5.5: A quarter of a hub geometry.

101 Phi Polarization 40. v:$ Ou c/) u 30. 20. -. i 10. 0. 0. 10. 20. 30. 40. 50. &O (deg) Figure 5.6: Backscatter for the hub geometry of Fig. cut). 5.5 (0 polarization, 0 = 0

102 Theta Polarization 40. 30. 20. 0 20. \ 0. 10. Mode Matcl. ----- MoM 0. 10. 20. 30. Oi (deg) 40. 50. Figure 5.7: Backscatter for the hub geometry of Fig. 5.5 (0 polarization, 4 = 0 cut).

C7. CD CD C#k Ilk CDA i

104 Phi Polarization 40. 30. 20. 10. 0. [111l l [ 1 111111.. I. I I I I I I I I I I I I I I I I I I I I I I I I I I. I................ I 5I, I.............................. I.....ll l l 0. 10. 20. 30. 40. 50. &i (deg) Figure 5.9: Backscatter for the blade geometry of Fig. cut). 5.8 (0 polarization, - = 0

105 Theta Polarization 40. 30. 11-\ C14 C-I 7 cn u P4 20. 10. 0. 0. 10. 20. 30. 40. 50. O1 (deg) Figure 5.10: Backscatter for the blade geometry of Fig. 5.8 (0 polarization, -b = 0 cut).

106 Phi Polarization 20. 10. "0 (. V~ 0. -10. -20. 0. 15. 30. 45. 60. Oi (deg) Figure 5.11: Backscatter of zero order propagating modes in a shorted inlet of 2A diameter (q$ polarization, 6- = 0 cut).

107 Theta Polarization 20. 10. CVU ^4 0. -10. -20. 0. 15. 30. 45. 60. Oi (deg) Figure 5.12: Backscatter of zero order propagating modes in a shorted inlet of 2A diameter (0 polarization, - = 0 cut).

CHAPTER VI SUMMARY, CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 6.1 Summary and Conclusions We investigated the problem of electromagnetic scattering from jet engine inlets, and proposed various analytical and numerical methods for an efficient calculation of their Radar Cross-Section (RCS). In Chapter II the Mode Matching Technique, a semi-analytical method, was applied to characterize a cylindrical waveguide terminated by a number of non-trivial engine-like configurations. The motivation of this study was to generate reference data for validating more general numerical techniques. Three different termination geometries were considered, including a termination consisting of an array of curved grooves. The latter has a close resemblance to the jet engine face whereas the stub and straight groove terminations served to validate the mode-matching implementation for this application. Scattering results were presented for each configuration and a serious effort was devoted toward their validation either by comparing them to other reference solutions (when available) or by examining the satisfaction of the boundary conditions at the interface of the termination. Comparisons with actual measured data were performed, demonstrating excellent agreement. 108

109 Since analytical techniques, such as the Mode Matching method, are applicable only to canonical geometries, we investigated the use of more versatile, numerical methods, capable of modeling arbitrary structures. In Chapter III the Adaptive Integral Method (AIM), an integral equation technique especially developed for large scale problems, was described. The relevant mathematical background was given in a rigorous manner, novel analytical calculations were performed and an efficient computer code was written for the free space case. Based on these simulations it was shown that memory requirements and computational complexity of AIM is much smaller than the Moment Method (MoM), especially for relatively flat surfaces, as in the case with engine blades, whereas accuracy is retained at a desired level. Several numerical results were validated for relatively flat and corrugated plates and were found to be in excellent agreement with MoM. Also, parallelization of the FFT algorithm demonstrated the impressive computational efficiency of the method in multi-processor environments. To further reduce the complexity of AIM for jet engine modeling, we exploited the inherent symmetry of the blade structure. In Chapter IV it was shown that substantial computational efficiency can be achieved by taking advantage of the blade periodic properties in much the same way done for Bodies of Revolution (BOR's). The jet engine geometry was characterized as a "Discrete Body of Revolution" (DBOR). Specifically, the following simplifications were achieved: (a) the modal scattering matrix was very sparse and its non-zero entries could be predicted a priori, leading to large storage savings; (b) the computational domain was reduced down to a single periodic sector of the scatterer (e.g. a single blade for the jet engine). This was done regardless of the geometry of the periodic sector or slice and as can be realized the computational savings are dramatic. Specifically, the number of unknowns is reduced

110 by a factor equal to the number of slices N8, with the corresponding computational savings being equal to N2. The above CPU time simplifications for a DBOR parallel those already known for regular BOR's and therefore the presented analysis can be considered as a generalization of the latter. Similar simplifications can be carried out for PDE simulations, but the one presented here for the integral equation formulation is more attractive for metallic DBOR structures, since the computational domain is restricted to the scatterer's surface and no cumbersome phase boundary conditions are needed to take advantage of the periodicity. Instead, for integral equations, the only modification is the introduction of the periodic Green's function. In Chapter V we applied AIM inside the jet engine inlet geometry. We showed that the "slicing scheme" developed in Chapter IV can be efficiently combined with AIM only if fundamental modifications of the pertinent grid are carried out. Specifically, it was shown that the Toeplitz property, the cornerstone of AIM efficiency, can be retained only if a cylindrical, rather than a rectangular grid is utilized. Fast Fourier Transform (FFT) can be used only along two dimensions, as opposed to the free space case where three-dimensional FFT is employed, but the important advantage is that only a single sector ("slice") of the engine needs to be modeled, requiring much smaller FFT dimensions. Several numerical results were presented and compared with reference solutions, showing the validity of the method. 6.2 Suggestions for Future Work Although the methods presented in this thesis were shown to be quite accurate. further improvements must be carried out before they yield RCS results for realistic engines of diameters up to 50 wavelengths. As far as the Mode Matching method is concerned, very little futher improvement seems possible, since non-canonical geome

111 tries are not amenable to analytical solutions. On the other hand, AIM is probably the method of choice, provided a speed-up of the algorithm is attained. Parallelization of the code is absolutely essential, since sequential multi-dimensional FFT's are extremely time consuming. Two different options can be considered: parallelization of the FFT only, or domain decomposition of the geometry and parallelization of the entire code. FFT parallelization has already been tested for the free space case, and the results were very satisfactory, but the second option should also be investigated, especially for arbitrary, not necessarily flat surfaces. As soon as parallelization and optimization of the code are complete, it is believed that AIM will be capable of performing RCS calculations for realistically large jet engine inlets.

APPENDICES 112

113 APPENDIX A MATHEMATICAL DETAILS OF THE MODE MATCHING TECHNIQUE A.1 Explicit expressions of the modes in the cavity The transverse fields corresponding to each mode are given by: (the longitudinal dependence has been suppressed; e stands for the electric and h stands for the magnetic field) Region 1. eTlnpM 1, cos( u - eTM = T1pl,np TM / in n} up p 7l,np p n - l, npP - sin nb U l,npJn (l1,npP) { sin b} (A.2) uTE NTE [ (JE4 cos no e -= n Nn~p Jn -l,npPP { - sinun f ~ + rTY j/ & rLp) { sin - } u] ( -+ 7,(TE, COsin nlj 9 (A.3)

114 T{E N sTE n \ TE r} I TE S+o up hlnp - 1,np, [^ npJn 7l'npP) C O co Up + " in TE1 cosnp n A. 4) + p - sinn } (A.4) where the TM, p are found from JT(y b) = 0 (A.5) J(^Eb) 0 (A.6) The normalization factors are arbitrary. They are usually defined by [5] [I -1 NTM t W l(,np- Y/),n p bcn] (A.7) NTE - n1 TE T EE bE2 - n2 (A.8) J1,np 2 1,ipb kZI,npZ3 b)2 (A.8) (A.9) 2 n=-0 Cn { 1 n: 0 (A.9) and on the basis of this normalization eIp hel, np h d2S=1 (A.10) Region 2 (Cylindrical Hub Termination). TM TM \y(TM f _TM _ \ j _YTM \y j TM M emq Yl2,mq [Y m (2,mqa) J 2m (\2,mqP)- m (2,mqa)m / ( 2,mq)]!TM t TM 2,mq/ J sin Ir * Nfqee { cosmn }up + + - [Ym (q 2 mqa) Jm ("2,mqP) - Jm (Y2,mMqa) Ym ("2,mqp)] * N2'mq;e -sinmw uc (A.11) 2,mq -sin MO

115 hMmq = p [Ymi (72,mqa) Jm ("Y2,mqP) -Jm (72,m'qa) Ym (72,mqP)] 12,mq mq [m mq) m Y2,mqP) - a2,mq 2,mqP)\ N2,mq (AU1) 2 [Ym (Mqa) Jm (Tp) - (^ qa) m ( TEp)] NTEM sinm U(A.12) e2,mq cosmj, +2,mq 7 [m (72,mqa) Jm (72,mqP) -JIM ('T a) Y ( TE )] N2q/~mq { cos } up +(A TEM TE M T h +, -, [Ym (2qa) Jm (72TqP) - Jm (7y2q a) Ym (y72mqP)] N2,mq {, (A.14) Ios m(pO eTEM = -NTuT (A.15) where the JM 7 are found from Ym'(yqa)Jm Y2qb) - Jm (<qa)Yi(2,mqb) = 0 (A. 17) Th -J abn= 0u (A.14) The normalization factors N N NEM are fouagain arbitrary, but are usuallyom defined by [40]

116 N (3 cm) I__ IA.9 Nmq oT(Cm{[iM 2 T2 ~-(y;qb2 rT2 n=O Reio 2, Mgrooe (goov arryeminto) e2,,mq 2~2m [' 7/2m] I2mr TEV\12ma bv~2mV b 2Mn, inv(n5- q 2,mq * ~ (A.23) TEa h~~mq = U[YL (72Lmq Ji'(y'q)-v(yq)Y(y'q] NTAI * fsi2( q)]~(A.24) TE kq T2,mqJ Jv~ 2,mqP tJ~ 2,mq'JV 2mqJ *NIM 2qqsin [v u, ]U

117 2,~q 2m [7 1TEm 1TEm hTE TE q (yr 2,mqa) IT j2qP) - J'1 (j72qa)Y 1~ Q2,mqP) *N2,E Cosq IoV [v (- On) UO) u TE,2mq 2, JCmqXq )Y(W/ ) JV Tq)YE Tqb - J j(fvb)i 7T a)= (A.25) (A.26) (A.27) (A.-'28') and m~r LI M = 0,1,72,. 2,mq N TE - /3~- I WE W (3 Em,Jj&L1m (A.29) (A. 30) 7- 2(3Em-)f 2 o5w [F (q 2] - [ (y~qa2] }f 2 m2 Em Z-j M#0O (A.31) (A.32) 0,, is the angular extent of each groove.

118 A.2 Closed Form Expressions for the Coefficients of the Mode Matching System of Equations (Cylindrical Hub Termination) The integrals defined in eq. ( 2.16), ( 2.17), ( 2.19), ( 2.20) can be evaluated in closed form through use pro of properties oBessel functions. Particularly useful in the evaluation is eq. (11.3.29), p. 484 of [41]. In the text we used only one index to characterize the modes, and also we dropped the specifications TM, TE, TEM, for clarity purposes. Here we define the integrals in a more detailed manner. We set rb 27T Pmqnp e2,mq ' el,nppdpdo Ja JO b r27' Unp - j - e l,np elnppdpdo b 27r Rmq,np h2,mq h hl,nppdpdo rb h 2, h h2,mq Vmq - J h2,mq - h2,mqpdpd~o Ja JO (A.3a3) (A.34) (A.35a) (A.36) In order to write the latter integrals in closed form, it is convenient to define the following functions: Q11) (x; c) - 2 1) [aJn ('ax) Jn- (x) - Jn-1 (ax) Jn (x)] + + 2 (a2 -1) [7Jn+2 (ax) Jn+l (x) - Jn+l (ax) Jn+2 (x)], 2 (a2 - 1} if n 0 0(11) (x; ) - 2 - [aJ2 (ax) Jl (x) - J1 (ax) J2 (a)] if n = 0 (A.37) (A.38)

119 Q(l2) (X; a) ax -2(a2 - 1) [aJn (aX) Yn- (X) - Jn1-l (aX) Yn (X)] ~ + 2 ax [aJri+2 (aX) Yri+l (X) - J17+l (ax) Yn+2 (X)], if n#O~ (A. 39)' Q(Il2) (X a) 2ax if n =0 Then it can be shown that: eRegion 1: TMXP modes, region 2: TMmq modes Pmq,np - [\TTM TM/31nP/32 ' - l'rp 2,nqWE C (A.40') rT T M._ TM y(TM a' IQ(il) T~bi,"Y '\ / (1 TM a,__l _ LYh \-2 nq n 2,nq} II TLM - ~12,niqaI TMJ L \jY 2nnq Y2,niq/ K1- TM TM\KTa TM TM (J'a (i Q1 2),nb TMb)'ln - Q(i2) TMq a, 71,n (A.41') L 7Y2,nq 72n Region 1: TMV~n modes, region 2: TEmq modes Pmcj~np nU 2nT%7Fmjn WCh -1 }n,~ TM~a 2 (A.42)I * Region 1: TMV~n modes, region 2: TEM mode FOO,n -= -NU Nf EM/14P 2wr Jo ([a) * Region 1: TEnp modes, region 2: TMmq modes (A.43) Pmq rnp = 0 (A.474)

120 * Region 1: TEnp modes, region 2: TEmq modes Pm - VTE VTE T' mq,n p l,np 2,nq 7mn Y' (Ta ) TE b1 'rp - Q(11) TE a, 71np 7n I2,nqa 2,nq T n 2,nqnqa; TE -, 72, nq 72,nq J [ / TE \ / TE \ | _ j,.,TEna Q (12) TE b lnp (12) TE a np \ A Jn 2,nq n 72,nqb; TE -n /2,nqa; TE (A45), 72,nq/ \ Y2,nq * Region 1: TEnp modes, region 2: TEM mode Pmq,np = 0 (A.46) Finally, on the basis of ( A.10) for all cases, Unp 2Znp (A.47) In the latter expressions, 8mn is the Kronecker delta and Znp is the characteristic impedance of the np mode. It is interesting to observe that only modes of the same order (i.e. order of the Bessel function) couple to each other, Also, TE modes in the first region couple only to TE mode of the second region, in agreement with the results in [42]. Similar closed form formulas hold for ( A.35) and ( A.36), but the explict expressions are not important, due to property ( 2.23).

121 A.3 Compact Expressions for the Coefficients of the Mode Matching System of Equations (Groove Array Termination) As opposed to the cylindrical hub termination, the coefficients of the matrices involved in the mode-matching system of equations cannot be evaluated in closed form. The ( integration can be performed analytically, but the p integration has to be carried out numerically. Again, we define the pertinent coefficients as P - 'P,mq,np Unp fb 0+Obw / e2,,E,mq ' el,nppdpdo Ja K K /b 2 e7r J J elnp * elnppdpdo o Jo (A.48) (A.49) To expres ies the coefficients in a copact form, it is convenient to define the folowing quantities and functions (t, m, n E No, v = mn/O,)): I(1), ~n m K,3n j(3) tK,nrn -- J cos n,, cos [v (o -;,)] do = -= sinc 2- w cos ( 2w + n + 2 2 i 2 + — sinc Ip Icos~p> +n~ o. (~w (n - v \ (n - v \ + wsinc w cos ( - w + n,;3,) 2 2 / 2 OK+OW = sin nq cos [, (4 -- q,)] doq = J^K ~>w. n + vo. n + v 0 = sinc w J sin p 2w + n, + 2 2 2 + -sinc ( - w sin ( I- w + n 2 2 -- J cos noq sin [v (43 -- q,)] do$ = K4 (A.50) (A.51)

122 - sine (f -Lqw) sin c/$, + n, - -7slnc (n Ow sin (n/qw +no,) (A.52) Irrn ] sin n i v(,]d - 2 ~sinc 2l qw Cos( 2 Oi$+n,, + ~ 2 sinc ~~O o O o (A.53) M(" (Xi, X2; a) j J71 at J, ' ~c (A.54) 42Qm (XI, X2; a) j J ()J,()d A.) AI$~ (Xr, X2; a) j J, (a~) J, ( -)d~ (A.50) The oeficintsdefnedimq.8- 4)cnb rte s r rf()2___ ( __ jY m (XyfXa) a) imP Yjll d2,ma<A2m59 f(3)2)l2,m I(') TM __ w ~T f, n )11 (TM TM bJ 1irip I V12mq~ (2) J nrm ~2,mqa< 2,mq 'I TMj'

123 - J (2,mqa) { -.(^TMa) { I (4) K,nm J(3) K,nm i TTM TM \ 1,np M(12)_ TM TM b1 np TM rIm \12,mq aY2,mq I TM,Y2,mq Y2,mq }^nm n 72,mqa 72,mqb; TM J 2,mq I(1) n,nm _JI(2) K,nm (A.62) * Region 1: TMl/np modes, region 2: TEmq modes Pnmqnp:: Km q,nip -- OTM NTMN TE l,np l,np 2,mq * { (ITEa) + y (<TE) { f(4) TM TM Knm I 7Yl,np a(11) TE TE b; l,np + J I(3) vTE ~ -nm I Y2,mq: 72,mq I TE I Knm J/ i2,mq \Y2,mq/ /(I ) m TM \ ),nm r(im) TE TE 1,rip I (2) r nm 7Y2,mq /Y2,mq TE I -r,nm 2) 2,mq/ (4) A TM TM\ n,nm 7Hl,np (12) TE TE 7. l,np _ I(3) j LTE nMrm 72,mqa q2,mq 7 TE K,nm 2'2,mq \t2,mq I() TM\ r()m E EV_(12) TEa TE b. "l,_np _63 __) nm 72,mqaK 7y2,mq TE f (A.3,) rc,nm \Y2, mq - J' (y4qa) - J, (mqa) { { * Region 1: TEnp modes, region 2: TMmq modes PKrnq,np = 0 * Region 1: TEnp modes, region 2: 'TEmq mod? (A.64) PKNmq,np =: NTE TE _ {YlN2,Nqa) npY (2,mqa) { *~ TE Y 1~2~{ { 1(2) K,nm IMl) -_I(3) Knm I(4) K,nm TE TE T El,np (1 ) TE T E 71,np / TE }nvA, TE a qTE b- - 4-rp TE, TE TE nm "12,mqa '2,mq I TE j 722, mq \ 2,mq/ } Anm 2qa Yqb; } J J ( '2, mq ) - J (.2,qa) { I(2) K rinm I(i) Knm } J- (2,q a){ _I(3) K,nm 1(4) K,rnm A.65)

124 Of course, on the basis of ( A.10) for all cases, Unp = 2Znp (A.66) Like in section A.2, Smn is the Kronecker delta and Znp is the characteristic impedance of the np mode. Again, TE modes in the first region couple only to TE mode of the second region, in agreement with the results in [42]. A.4 The Eigenvalue Symmetry Property of the G.T.E. Matrix The following property holds for matrix [M] of eq. ( 2.46): If A is an eigenvalue of [M], -A is also an eigenvalue of [M]. The proof is based on the three following Lemmata: Lemma 1: For any matrix [A], det [A] = det [AT] (Proof well-known) Lemma 2: Let [A/j] be n x n complex matrices. Then [All] [A12] [All] - [A12] (A.67 [A21] [A22] - [A21] [A22] Proof: Let [I] be the n x n unit matrix. Evidently E [All] - [A2] 1 - [I] [0] [ [Al,] [A12 ] [ [I] []. - [A21] [A2] ] I] [A21 [A22]] [ [0] -] (A Given that [I] [0] -1 (A.69) [0] — [I]

125 and that the determinant of a product of matrices is equal to the product of their determinants, the proof of the lemma is obvious. Lemma 3: Let [A], [B], [C] be n x n complex matrices. Then [A] [C] [B] [C] [C] [B] I- [C] [A] ( Proof: Similar to Lemma 2. Note that [B] [C] [ [o] [I] [A] [C] [ [0] [I] (A.71 [C] [Al [A] [B L J L [] [0] [C] [B] [I] [0] [] J and consequently ( A.70) holds since [0] [I] 1 (A.72) [I] [0] and since the determinant of a product of matrices is equal to the product of their determinants. Using the above three Lemmata, we can develop the proof of the eigenvalue symmetry property as follows: Proof of the Eigenvalue Symmetry Property: Let [A] [I)] (A73) [M]- [D] -[A]T (A.3) where [D] is a diagonal matrix. If A is an eigenvalue of [M], then by definition [A] - A [I] [D] ( [D] -[A]T -A[ [ =I] 0 -

126 Further, by using Lemma 1, and since [D] is diagonal, we have [A]T - A [I] [D] [D] = 0 -[A] -- A [I] (A.75) or, equivalently - [A]T + X [I] -[D] - [D] [A] + A [I] A.76) yields Next, from Lernmma 2, ( (A.76) (A.77) - [A]T + A [I] [D] [D] - = [A] + A [I] and finally Lemma 3 gives [A] + A [I] [D] [D] -[A]T +A[I] (A7 This result shows that -A is an eigenvalue of [M], Q.E.D. A.5 Application of the Mode Matching Technique to the Generalized Modes within the Compressor In this section we apply the Mode Matching Technique to the geometry of Fig. 2.8 in order to evaluate the scattering matrix of the curved blades array. In Region 1 the transverse fields are given as a superposition of cylindrical waveguide modes: 00 Et = [ape1,p exp {j/31,pz} + bpe1,p exp {-j,p}] (A.79) p=1

127 o00 Ht = [aphl,p exp {j1,pz} + bph1,pxpxp {-1,pZ}] (A.80) p=l where ap,bp are the coefficients of the incident and the reflected fields due to the termination. Similarly, according to ( 2.39) the transverse fields in the Kth groove of region 2 are expressed as superpositions of the generalized modes: 00 E2K = EV(Z)e2,,(pz) (A.81 ) i=1 00 2, = Ir,,(z)h2,,(p,; Z) (A.82) i=1 From ( 2.51) we deduce that in the Kth groove {V,} = [Kll] ez[] { + K12] e-z[A] {d} (A.83) {JK} [K21] e {c,} + [K22] e-z[A] {d (A.84) where {c,}, {dj} are vectors of arbitrary constants. To evaluate {b} {b1, b2..}T in terms of {a} -{a,a2,a...}T we must impose the continuity of the tangential fields at the interface. Furthermore, the tangential electric fields must vanish on any PEC surfaces. Let us locate the origin of the coordinate system at the interface between the two regions. Enforcing continuity of the tangential electric field at the interface yields Eitz=o Et z=o (A.85) Ht1 lz=o = Haz=o (A.86) \Vp C [a, b], V' C [OK (0), K (0) + W], Vn { 1,..., J}

128 On combining ( A.85) with ( A.79) and ( A.81), taking the dot product with the modes el,q and integrating over the whole interface, we obtain the system J [U] ({a} + {b}) = [M,] ([K1l] {c,} + [K12] {d,}) (A.87) where rb r(b^K+w Mrqp - q e2,K,q ~ el,ppdpdo (A.88),b r27r Umn Smn ej el,p epdpdqc (A.89) and Smn is the Kronecker delta. In a similar way, on combining ( A.86) with ( A.80) and ( A.82), taking the dot product with the modes h2,K,q and integrating over the portion of the interface common to both regions, we obtain the system [QK] ({a} -- {b}) = [V,] ([K21] {c} + [K22] {d)K) Vr/ 1,.., J (A.90) where fb ykK+0w Qqp = J h2,,, hl,ppdpd) (A.91) J a /o rb r27r Vmn mn h2,m, h2,npdpd) (A.92) Ja JO and 8mn is the Kronecker delta. If [Srend] is the scattering matrix of the termination corresponding to the,th groove, then, by definition,

129 [K12] e12[A] {d} = [Sn,end] [Kll] e-12[A] {c,} (A.93) where 12 is the length of the second region. Expression ( A.93) provides a relation between {ci} and {d,}. To find {b,} we must solve systems ( A.87), ( A.90) and ( A.93). Obviously, {be} will be a function of {a,}, and to avoid repeating the solution for every excitation, it is customary to instead compute the scattering matrix [S] defined by {b} = [S] {a,} (A.94) To help define the scattering matrix, we introduce some auxiliary definitions: [R,] - e-[A]2 [K12]- [S,,,nd] [Kll] e-[A]12 (A.95) [F,] - [K11] + [K12] [Rn] (A.96) [G,] [K21] + [K22] [Rn] (A.97) [LK] - [V,]-1 [Q] (A.98) Then, from ( A.87), ( A.90), ( A.93) and ( A.94) it follows that [S] ([W] + [I])1 ([W]- [I]) (A.99) where [W] E [LK] [FF ] [K,]-' [L] (A.100) r-=1

130 and [I] is the unit matrix. A.6 The Exponential Matrix Solution of the Generalized Telegraphist's Equations. In this section we give an explicit form of the matrix exponential appearing in ( 2.49) when the blades are straight. It is shown that the well-known solution of the conventional Telegraphist's Equations is recovered. In the case of straight blades, F'(z) = 0. Hence [M - [ ] [D] (A.101) [D] [0] One observes that [M] = [O [D]2 = [Q (A.102) [[0] [D]2 Jwhere [Q] = [ [D] J (A.103) Thus exp [M] = [I] + [M]2 + [M] + [M]3 +... = {[I]+ [M]2 + [M]4 +...}+ {M][M] 3 + [ [M]...} = {[I]+ [Q]2 + [Q]4 + } + [M] [Q]+ ] + [ Q] + } [Q ] +=cosh [Q] + [M] sinh [Q] [Q]

131 cosh[D] sinh [D] (A. 104) sinh [D] cosh[D] (A.104) The latter equation, in conjunction with ( 2.45) and ( 2.49) yields V,,i = c,,ie-J/iz + dn,ieJiz (A.105) Zol,i c,,ie-j3iz -- d,ieie3i (A.106) where ci- [V,, (0) + ZoI,, (0)] (A.107) d,- [VI (0) - Zol (0)] (A.108) This is the well-known solution of the conventional Telegraphist's Equations. A.7 Explicit Radiation Coefficients for Circular Inlets In the near zone, there is no simple closed form expression for the radiated field due to a single outgoing waveguide mode. However, it is possible to derive closed form expressions for far zone observations. From [5], the far zone field due to each mode is given by Erad = -E1 ncos ~ E sin no} u E - (A.109) ad -sinn4 J no cos no J r Eo = Eok +Eou (A.110) EO = Ek+E Eu (A.1Ill)

132 where Eok,Eok are associated with the contribution from the Kirchhoff integral and Eo,, EsU correspond to the contribution from the equivalent fringe rim current (Ufimtsef type). Explicit expressions of Eok, EOk, Eou, Eou are as follows: TMnp modes: Eek = jikNTM nl b. sin 0 2 (cos - cs) JJn (kb sin 0) (A. 1 12) 2 [COSJ ( b - COS 0/2,np c E;k = 0 (A.113):n T OS (0/-N) _ COS (ITnM/2) — j lnp Jn (l,npb) ~ lnpOn lnp cos - COS 0 [/2TM COt 0 kTM 0 * [n2 dTIM sin 2,kb Jn (kbsin ) + + k71,Pbsin J" (kbsin0)] (A.114) Eu = nNTM I (TMb cos (0/2) - cos (,PT/2) -- 1,np -n b,1,np CO VjTh - C 0 rTM 1* {/T n sin 1 (kb sin0) + + k^TMb cObs]in2 kb [Jc J(kbsin0) kbsin f } (A.115) TEnp modes: 1 I,np - E,= rjkZN[l TEn 1 + cos 0 cosn (kb (in \) 2 sin JJn (1kb sin 0) (A.116)

133 E,(k = jnkZN nrp b. snp j (,pb Jn (kbsin 0) (A.117) 2 (cos 4pTE, - cos o) v n E jZI nc Cos (0/2) - cos (Ts 2) ^ = 3,np n^ b) - ^_- COS *T o - Cnp OS TE 0 / Ilnp csin J" (kb sin 0) -- - k'y 'npbsin 2 kbn -) (A.118) TETE\cos (0/2) -n ZN 'MTE Tbln Zu I on hn pCS T l,npbTE - COS 0 kTE2,np sinl (kbsin0) - - 2/31,psin - [J (kb sin 0) - } (A.119) where, by definition oTM,TE TM TE _ ip (A.120) COS YFln~ --

134 APPENDIX B MATHEMATICAL DETAILS OF THE ADAPTIVE INTEGRAL METHOD (AIM) B.1 The AIM Impedance Matrix In this section explicit expressions for the elements of the AIM impedance matrix are given, and we also explain how the Toeplitz property can be utilized, so that fast calculations of matrix-vector products can be performed. In the text it was proven that the integral equation is written as Ei fmdv = jwA -fvd - oV, * fmd (B.1) By using ( 3.2) and by substituting the basis functions by the auxiliary set Im for the far field interactions, the vector potential integral in ( B.1) is equal to N M3 M3 J = jw In E E (AnpA + ApA + A Aj) G(rp, rm) (B.2) n:=l p=l q=l where rmp and rnq are the locations of the grid nodes and

135 G(rnp,rmq) P {-Ij (B-.3) r47rlr - rmq[ We need to express ( B.2) in matrix form. To do so, we extend the summation over p and q to the whole grid that encapsulates the geometry, thus eliminating indices m and n from the arguments of the Green's function. Furthermore, matrices [Ar], [AY], [AZ] are defined in such a way that their elements Amq are zero unless q is one of the M3 grid nodes around the m element (see fig. 3.4). If Ng is the total number of grid nodes, we obtain N Ng Ng JA =- jW ^ ^^I+A^A^ ~AZYq} +A ^rq) jA n In A Anp mq + AnpAmq G(rp, r) n=l p=l q=1 j [EA] {I} (B.4) where {I} - [I1,I2..., N]T (B.5) and [E ] =[Ax][G][A ]T + [A] [G][A ]T + [Az][G][Az]T (B.6) where exp {-pjc k,-) Gpq -e {4 pq (B.7) 4Wrrpq with

136 rpq == rp-rq (B.8) Similarly, the scalar potential integral in ( B.1) is equal to J = [E]] {I} (B.9) jwe where [E ] = [Ad][G][Ad] (B.10) The [L] matrices in ( 3.26) are defined by [L(1)] j- w/,[A] (B.11) [L(2)] - Wl[A~] (B.12) [L(3)] - jwl[AZ] (B.13) [L(4)] - [Ad] (B.14) Matrix [G] has the Toeplitz property, provided the element numbering is based on six indices, instead of two. Indeed, suppose that h is the grid step and that p = h (pi +p2y +p3z) (B.15) rq = h (qlX + q2Y + q3) (B.16) where pi, qi are integers. Then a 3 tensor can be defined by Z\ Z 3

137 GPl,P2 73 - Gpq (B.17 ql,q2,q3 - 1 ) for the aforementioned [G] matrix. It is obvious from the definition that the elements of this tensor are actually functions of the differences pi - ql,p2 - q2,P3 - q3, which means that tensor products of the form GPl,P2,P3,2 q3 (B.18 ql,q2,q3 ($ 8 can be written as E GPl -qlP2 - -2 p3 — q3 q1'q2'q3 (B.19) ql q2 q3 which is evidently a three-dimensional discrete convolution. By defining the periodic extension of the three-dimensional sequence G and after suitable zero padding, all matrix-vector products involving the [G] matrix can be efficiently calculated via three-dimensional FFT subroutines. An FFT-based computer routine for the rapid calculation of the tensor products in ( B.18) has been developed, validated and included in the AIM code. B.2 The Moments of the Basis Functions In this section it will be demonstrated how the moments defined in ( 3.18) can be efficiently calculated by closed form expressions, thus avoiding time-consuming numerical integrations. Assume that the global coordinate system is denoted by (x, y, z), while an addi

138 tional local coodinate system (x~, y~, z~) is defined for each one of the triangles T~ that make up a tetrahedral patch (see Fig. 3.5). The origins of the local coordinate systems are located at the free vertices of the triangles and their z~ axes are perpendicular to the surface of the triangles. Given the definitions of the moments in ( 3.18) and the basis functions in ( 3.6) and ( 3.16), it is clear that the calculation reduces to integrations of the form M q = 2 / (r- r+)(z+)(x -- xza)q(y- ya)2(z- za) 3dxdydz + + 2 J- J(r - rl)6(z-)(x -- xa)q(y - ya)2( - za)q3dxdydz (B.20) and due to the properties of the delta function it is obvious that the initial volume integral is decoupled into two separate surface integrals over T+ and T-. Let us con — centrate on any one of them, dropping the superscripts ~. From now on we will use primes to denote the local coordinate system on any triangle. To calculate any of the integrals in ( B.20) we need to express all quantities in the local coordinate system. Let [C] be the directional cosines matrix that transforms (x', y', z') to (x, y, z). Then r- r = (x - x1)x + (y - Yl)Y + (Z - 1)z = x - xZ1 X = [x ] - Y =- [,, z] [C] - y' (B.21) Z -Z1 Z! Next, we transform the arbitrary triangle T into a right isosceles To, by use of the transformation (see Fig. B.1)

139 (X'2,y'2) 1 1 (x'3,y'3) Figure B.1: Transformation to a right isosceles triangle

140 Y [2 Y3 ] (B.22) y{x [ ] { } The new orthogonalization coordinates E, 77 simplify the integrations over the triangle, since it can be shown that 1 - P!q /& / 7 /1q rf 'dIdrd =- (B.23) T 0Pqdd Jo 1Jo ~Pdd - (p + q - 2)! (B.23) where p,q are nonnegative integers, while de( 2('l) i) (B.24) where A is the area of he triangle. Using the fact that on the triangle surface z' = Z = = Y 0, ( B.22) can be written as ( x' X1 X<2 X3 0 Iy' - y 3 yy ~ 'J (B.25) ~I ZZ1 2 23 1 After some simple algebraic manipulation, ( B.21) and ( B.25) yield r-r1 = [(x2 - x1)( + (23 - 1)7] + [(Y2 - Y1) + (Y:3 -yl)r7] y + + [(Z2- Z1) + (z3-z l)]] Z (B.26) Similarly we obtain

141 r- ra = r- rl + r -- ra = (X21 + X317 + Xla) + + (Y21 + y31l + Yla) + + (Z21 + Z311 +~ Zla) (B.27) where we have defined Xrs = Xr - Xs Yrs y= r ~- Ys, Zrs ~ - to simplify the notation. To calculate the moments, it is sufficient to evaluate the integral J(Q) = (oA + A0) (I ( + ~oq) + q4) d dq (B.28) q —0 where Q is the sumn of the moment indices. For example, to calculate the x component of M210 we set Q = 2 + 1 + 0 = 3 and A(O)= AX21, = X31 (0) - n (o) = (o) I f10 1 - 0 0 0 -1 (1) (1) (1) /10 )= X21, /01 = X31, / 00= Xla (2) (2) (2) 1X10 X21,//01 X31 /00 - Xla P1) = Y21,i (1) = y31, oo) = Yla (B.29) The integral J(Q) can be evaluated recursively in the following way: For Q = 0 we have N~O) +A(o) j() + (B.30) 6

142 The integrand of J(Q) is the polynomial Q+1 Q+1 p(Q) >3 Y>3 A8?)) iuJ i=0 3=O and the integral is (see ( B.23)) (B.31) Q+1 Q+1 j(Q) = E E i=0 j=:0 i!jI I](I +j-F 2)! (B.32) The integrand of j(Q+') is the polynomial where Q+2 Q+2 p(Q+l) - p(Q) ( (Q+1)~ + /(Q+1)r + (Q+1))AQ+) 10i=0 j=0 — QljQ)Q+l)1) + + if i> O0 (B.33) (B.34) and A =+' 0 if I < 0 or j < 0 (. 237 - (B-35)

143 APPENDIX C MATHEMATICAL DETAILS OF THE PERIODICITY EXPLOITATION C.1 The Cylindrical Waveguide Dyadic Green's Function Following [34] we define two sets of cylindrical vector wavefunctions, namely M m (r; /3E) Nnrn(r; /3, ) - V x [ Jill( p)ee- M -V x V x [JInIn] (C.1) (C.2) The functions M,,,m correspond to TE modes and the functions Nnm are associated with TM modes. If the waveguide radius is a, then JTM,'TE and QTIMTE are defined by J TMa = 0 Jn (Nm a) o Jn('YMa) = 0 ( TM,TE)2 + (3M,TE2 = k2 )nm +,nm, -- (C.3) (C.4) (C.5)

144 Carrying out the vector operations, Mnm and Nnm can be written explicitly as TE jn rlTE Mm(r; ) = Jlnl(pTE e-jn jzp - TEJ l(,TEYmp)eJ e-J z (C.6) TMTM M "/ri m i hm it TM ^\,_j~ TMZ + n J (7 mp)ein e-J +z (C. 7) The dyadic Green's function of the electric type satisfying Dirichlet boundary conditions (first kind) is given by G (r, r') =- (r-r') + 00 00 + S S Cnj Mn7n(r; ~!3rn)M-nm(r'; }fl3m) 77m=1 n=-oo + 5C C Nnm (r; ~ )Nnm(r; T3m) (C.8) m=1 n=-oo where the upper sign is used for z > z' and the lower sign for z < z'. The constants are given by C -j{2r [(yn7a) 2-n2] J(7m a):l} (C.9) C = -j{27(r m ( a)2l (YTN a)] n-}/ (C.10) Notice that exponential, and not trigonometric ependence on is used. Notice that exponential, and not trigonometric dependence on 4>is used.

145 C.2 Elements of the Impedance Matrix for PEC Scatterers In the text we showed that the impedance matrix entries for PEC terminations are given by A;mfr) -= kZ f(-r) G(R- r, Rm r') R f(r') dd2S' (C.11) where R is the rotation dyadic. To eliminate the rotation dyadics we can use the properties of G as defined in Appendix I. First, we express Mnm(r) as 3 Mnm (r)= E M( (r)ei(r) (C.12) i=l where {&1, &2,&3} = {3p, t>, z}. It then follows that 3 Mnm(R. R r= E M (R. r)i(R. r) = i=1 3 3 = eJss > Mi(r)R- s ei(r) = eJ E M3 ei.(r) R- = i=l i=l = eJnKssMnm(r) R (C.13) Similarly Nnm(R r) = enKs'N m(r) - R (C.14) Next, we define the principal part A^3 of ( C.11) as the integral over the modes only (not including the delta function term in ( C.8)). By invoking ( C.13) and ( C.14), expression ( C.11) yields

146 oo oo A(nms) = jkZ > E3 e-jn(ns-ns)s. pq m=1 n=-oo { f(r) Mnm(r; ~+3m)M_, (r'; F/3TE). fq(r')d2Sd2S' + So o Cnm (r- n n( f; (r'; ) )d Sd2S } (C.15) ~ 0 0 C.3 Scattered Field from PEC DBOR's It is possible to express the scattered field from a perfectly conducting DBOR in a very compact form involving the currents of the reference slice o0. Given an incident mode of order ni, the scattered field is generally given by ( 4.6). Using the properties of the Dyadic Green's function and the expansion ( 4.14) together with ( 4.25), after some algebra similar to Appendix II we deduce that for a point r away from the scatterer, such that z > z' for any source point r', 00 oo Enr) = -jkZ E [Bnmnnm(r;L 3T) + B Tnm.Nnm(r;/[ i )] (C.16) m=l n=-oo where Q Bnmn AN (n, ni)CTnE / Iq~) M- nm(r;-) fq(r)dS (C 17) g m,n — nm nm q=1 0 Q TM _ECTM L f i)Crn3 I2S Bnmnni -- /N(11 Ti)Cnm I0)i N- nm(r'; - fq(r )d (C.18) rin nm -i )( q=1 0 Evidently, ( C.16) yields the scattered field in terms of TE and TM modes, whereas the coefficients of the expansion (i.e. the elements of the modal scattering matrix) are conveniently given by ( C.17) and ( C.18) in terms of the currents (~) on slice q Q,ni

147 So only. A similar expression holds for dielectric scatterers. We emphasize again that due to the term AN,(n, ni) only a limited set of scattered modes are returned.

BIBLIOGRAPHY 148

149 BIBLIOGRAPHY [1] T. W. Johnson and D. L. Moffat. Electromagnetic scattering by open circular waveguides. Technical Report 710816-9, The Ohio State University, ElectroScience Laboratory, Columbus OH, Dec. 1980. [2] Ching-Chao Huang. Ray analysis of EM backscatter from a cavity configuration. PhD thesis, The Ohio State University, ElectroScience Laboratory, Columbus OH, 1982. [3] Ching-Chao Huang. Simple formula for the RCS of a finite hollow circular cylinder. Electronics Letters, 19(20):854-856, Sept. 1983. [4] C. S. Lee and S. W. Lee. RCS of a coated circular waveguide terminated by a perfect conductor. IEEE Trans. Antennas and Propagation, 35(4):391-398, April 1987. [5] A. Altintas, P. H. Pathak, and M. C. Liang. A selective modal scheme for the analysis of EM coupling into or radiation from large open-ended waveguides. IEEE Trans. Antennas and Propagation, 36(1):84-96, Jan. 1988. [6] H. Ling, R. C. Chou, and S. W. Lee. Shooting and Bouncing Rays: Calculating the RCS of an arbitrarily shaped cavity. IEEE Trans. Antennas and Propagation, 37(2):194-205, Feb. 1989. [7] D. J. Andersh, M. Hazlett, S. W. Lee, D. D. Reeves, D. P. Sullivan, and Y. Chu. XPATCH: A high-frequency electromagnetic scattering prediction code and environment for complex three-dimensional objects. IEEE Trans. Antennas and Propagation Magazine, 36(1):65-69, Feb. 1994. [8] R. J. Burkholder and P. H. Pathak. High frequency asymptotic methods for analyzing the EM scattering by open-ended waveguide cavities. Technical Report 719630-3, The Ohio State University, ElectroScience Laboratory, Columbus OH, Sept. 1989. [9] P. H. Pathak and R. J. Burkholder. High-frequency EM scattering by openended waveguide cavities. Radio Science, 26(1):211-218, Jan.-Feb. 1991. [10] Z. N. Lu and R. Bansal. Boundary relaxation method for the analysis of terminations in ducts. J. of Electromagnetic Waves and Applications, 6(12):1689-1703, 1992.

150 [11] D. C. Ross, J. L. Volakis, and H. T. Anastassiu. Hybrid finite element-modal analysis of jet engine inlet scattering. IEEE' Trans. Antennas and Propagation, 43(3):277-285, March 1995. [12] D. C. Ross, J. L. Volakis, and H. T. Anastassiu. Overlapping modal and geometric symmetries for computing jet engine inlet scattering. IEEE Trans. Antennas Propagation, 43(10):1159-1163, Oct. 1995. [13] R. Mittra and S.W. Lee. Analytical Techniques in the Theory of Guided Waves. Mac Millan, 1_971. [14] Y. C. Shih. The modal-matching method. Chapter 9 in Numerical Techniques for Microwave and Millimeter-Wave Passive Structures, T. Itoh, ed., Wiley, 1989. [15] H. T. Anastassiu, J. L. Volakis, and D. C. Ross. The mode matching technique for electromagnetic scattering by cylindrical waveguides with canonical terminations. Journal of Electromagnetic Waves and Applications, 9(11/12):1363-1391, Nov./Dec. 1995. [16] H. T. Anastassiu, J. L. Volakis, D. C. Ross, and D. Andersh. Electromagnetic scattering from simple jet engine models. IEEE Trans. Antennas and Propagation, 44(3):420-421, March 1996. [17] E. Bleszynski, M. Bleszynski, and T'. Jaroszewicz. A fast integral equation solver for electromagnetic scattering problems. In IEEE Int. Conf. AP-S Digest, pages 416-419, Seattle WA 1994. [18] E. Bleszyrski, M. Bleszynski, and T. Jaroszewicz. AIM-adaptive integral method compression algorithm for solving large scale electromagnetic scattering and radiation problems. Radio Science, 31(5):1225-1251, Sept.-Oct. 1996. [19] T. S. Chu, T. Itoh, and Y. C. Shih. Comparative study of Mode Matching formulations for microstrip discontinuity problems. IEEE Trans. Microwave Theory and Techniques, 33(10):1018-1023, Oct. 1985. [20] R. Safavi-Naini and R. H. MacPhie. Scattering at rectangular-to-rectangular waveguide junctions. IEEE Trans. Microwave Theory and Techniques, 30(11):2060-2063, Nov. 1982. [21] R. R. Mansour and R. H. MacPhie. Scattering at an N-furcated parallel-plate waveguide junction. IEEE Trans. Microwave Theory and Techniques, 33(9):830 — 835, Sept. 1985. [22] J. D. Wade and R. H. MacPhie. Scattering at circular-to-rectangular waveguide junctions. IEEE Trans. Microwave Theory and Techniques, 34(11):1085-1091, Nov. 1986.

151 [23] R. H. MacPhie and K.. L. Wu. Scattering at a junction of a rectangular waveguide and a larger circular waveguide. IEEE Trans. Microwave Theory and Techniques, 43(9):2041-2045, Sept. 1995. [24] G. V. Eleftheriades. Analysis and design of integrated-circuit horn antennas for millimeter and submillimeter-wave applications. PhD thesis, University of MIchigan, Radiation Laboratory, Ann Arbor MI, 1993. [25] R. F. Harrington. Time-harmonic electromagnetic fields. Mc Graw-Hill, 1961. [26] P. H. Pathak and R. J. Burkholder. A reciprocity formulation for the EM scattering by an obstacle within a large open cavity. IEEE Trans. Microwave Theory and Techniques, 41(4):702-707, April 1993. [27] M. Leroy. On the convergence of numerical results in mode analysis. IEEE Trans. Antennas Propagation, 31(7):655-659, July 1983. [28] R. J. Burkholder. The Ohio State University, ElectroScience Laboratory (personal communication). [29] G. Crabtree, W. Huegle, and D. Salisbury. Rcs compact range test results for a set of simplified engine face models. Technical Report TM94226, GE Aircraft Engines, One Neumann Way, Cincinnati OH, Nov. 1994. [30] G. Reiter. Generalized Telegraphist's Equation for waveguides of varying crosssection. Proc. I.E.E., 106B(suppl. 13):54-57, 1959. [31] Gilbert Strang. Linear algebra and its applications. HBJ, Third Edition, 1988. [32] Pamela Haddad. Scattering by non-planar surfaces. Directed Study, University of Michigan, Radiation Laboratory. [33] R. Lee and T.T. Chia. Analysis of electromagnetic scattering from a cavity with a complex termination by means of a hybrid ray-FDTD method. IEEE Trans. Antennas and Propagation, 41(11):1560-1669, Nov. 1993. [34] Chen-To Tai. Dyadic Green Functions in Electromagnetic Theory. IEEE Press, 1994. [35] S.M. Rao, D.R. Glisson, and A.W. Wilton. Electromagnetic scattering by surfaces of arbitrary shape. IEEE Trans. Antennas and Propagation, 30(3):409 — 418, May 1982. [36] J.R. Mautz and R. F. Harrington. Radiation and scattering from bodies of revolution. Appl. Sci. Res., 20:405-435, June 1969. [37] J.J.H. Wang. Analysis of a three-dimensional arbitrary shaped dielectric or biological body inside a rectangular waveguide. IEEE Trans. Microwave Theory and Techniques, 26(7):457-461, July 1978.

152 [38] R. E. Collin. Field Theory of Guided Waves. IEEE Press, 1991. [39] P. R Rousseau and R. J. Burkholder. A hybrid approach for calculating the scattering from obstacles within large, open cavities. IEEE Trans. Antennas and Propagation, 43(10):1068-1075, October 1995. [40] N. Marcuvitz. Waveguide Handbook. Krieger, IEE Press, P. Peregrinus Ltd., 1986. [41] M. Abramowitz and I. E. Stegun. Handbook of Mathematical Functions. Dover, 1970. [42] G. G. Gentili. Properties of TE-TM Mode-Matching techniques. IEEE Trans. Microwave Theory and Techniques, 39(9):1669-1673, Sept. 1991. [43] E. Bleszynski. (personal communication). [44] J.L. Karty, J.M. Roedder, and S. D. Alspach. CAVERN: A prediction code for cavity electromagnetic analysis. IEEE Trans. Antennas and Propagation Magazine, 37(3):68-72, June 1995. [45] Mikhail Srnelyanskiy. EECS, University of Michigan (personal communication).