Report RL 935 HYBRID FINITE ELEMENT MODAL ANALYSIS FOR JET ENGINE INLET SCATTERING Dan Ross John Volakis Wright Laboratories/AARA Wright Patterson AFB OHIO 45433-7001 March 1996 RL-935 = RL-935

FORWARD This report describes portion of the jet engine analysis work funded by Wright Laboratory/AARA, Wright Patterson Air Force Base. The project was under the Direction of Maj. Dennis Andersh and was funded from 1994 to 1996. This report describes theory and development relating to a hybrid finite element method for characterization of jet engines with periodic blade terminations. This report is also submitted by the first author (DCR) as a dissertation in partial fulfillment of the requirements for the Ph.D. degree at the University of Michigan

MISSING

TABLE OF CONTENTS DEDICATION.......................................... ii ACKNOWLEDGMENTS................................... iii LIST OF FIGURES.................................... vii LIST OF TABLES......................................... x LIST OF APPENDICES.................................... xi CHAPTER I. INTRODUCTION................................... 1 1.1 M otivation............................................ 1 1.2 Jet Engines.............................................. 3 1.3 Available Data and Modeling Techniques...................... 6 1.4 Objective.............................................. 7 II. FULL THREE DIMENSIONAL HYBRID FEM MODAL IMPLEMENTATION............................... 9 2.1 Introduction........................................... 9 2.2 Finite Element-Modal Formulation........................... 12 2.2.1 Convergence of the Modal Solution....................... 18 2.2.2 Finite Element Analysis.............................. 20 2.2.3 Implementation of the FEM Solution..................... 26 2.3 Examples............................................ 30 2.3.1 Example 1....................................... 30 2.3.2 Example 2.......................................... 33 2.3.3 Example 3.......................................... 36 2.4 Sum m ary............................................ 39 v

III. OVERLAPPING MODAL AND GEOMETRIC SYMMETRIES..................................... 42 3.1 The Need for Drastic Computational Scaling.................... 42 3.2 Overlapping Modal and Geometric Symmetries........................... 45 3.3 Phase Boundary Conditions for 3-D FEM with Node-Based Elements............................................ 51 3.4 Example................................................ 55 3.5 Summary..............................................56 IV. THREE DIMENSIONAL, EDGE-BASED, FINITE ELEMENT ANALYSIS FOR DISCRETE BODIES OF REVOLUTION................................... 59 4.1 Introduction to Edge-Based Discrete Body of Revolution Analysis 6.. 0 4.2 Phase Boundary Conditions for Wave Electromagnetic Problems Using Edge-Based Elements............................... 64 4.3 Validation and Results.................................... 75 4.4 Summary................................................ 83 V. EFFICIENT COMPUTATION OF JET ENGINE MODULATION.................................. 88 5.1 Jet Engine Modulation.................................... 88 5.2 Modal Scattering Matrix.................................. 91 5.3 Efficient Computation of Modulation.......................... 94 5.4 Validation............................................. 99 5.5 Summary................................................ 101 VI. SUMMARY AND SUGGESTIONS FOR FUTURE WORK........................................ 105 6.1 Summary............................................. 105 6.2 Suggestions for Future W ork................................ C106 6.2.1 Task 1............................................. 107 6.2.2 Task 2............................................. 108 APPENDICES.............................................. 109 BIBLIOGRAPHY......................................... 133 vi

LIST OF FIGURES Figure 2.1 Hybrid jet engine inlet analysis................................... 2.2 Amplitudes of the elements of the generalized scattering matrix calculated by mode-matching for a hub termination............................... 2.3 Maximum mode index (n) of traveling modes vs. electrical radius.......... 2.4 Profile of BOR cavity used to test convergence of modal solution......... 11 16 20 21 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 Convergence of far field with increasing mode index n (2 GHz).......... 22 Convergence of far field with increasing mode index n (4 GHz)........... 23 Convergence of far field with increasing mode index n (6 GHz)........... 24 Transmission line model of cylindrical waveguide absorber.............31 Optimum cylindrical waveguide absorber............................. 32 Performance of optimum absorber for the 77 traveling modes in a guide with a 2 wavelength radius......................................... 33 Results for example 1: a shorted cylindrical cavity...................... 34 Results for example 2: a stub terminated inlet......................... 35 2.13 Results for example 3: a ridge terminated structure.. 37 3.1 Illustration of the hybrid finite element-modal analysis of jet engine inlet scattering................................................... 3.2 A unique slice of an engine-like termination.......................... 3.3 A periodic signal with progressive phase advance...................... 3.4 Reference signal used to construct S( )........................... 43 46 49 49 vii

3.5 Calculated scattered field Re(E ) for TM01 excitation of a shorted inlet using a four degree slice.......................................... 56 3.6 Calculated scattered field Re(E ) for TEBl excitation of a shorted inlet using a four degree slice.......................................... 57 4.1 A cylindrical guide terminated by a fan-like structure................... 65 4.2 A guide terminated in a set of straight blades (left). For this geometry a symmetric mesh for one lobe (right) can be generated with coincident nodes and elements on face 1 and face 2................................... 66 4.3 A guide terminated in a set of curved blades (left). For this geometry a symmetric mesh for one lobe (right) having coincident nodes and elements on face 1 and face 2 cannot be generated............................. 66 4.4 A unique slice of the FEM domain for a fan-like structure................ 69 4.5 Diagram used for physically deriving the phase boundary conditions....... 72 4.6 The mesh for an 8 degree lobe of a shorted cylindrical guide with an electrical radius of 6 wavelengths................................... 76 4.7 Field distributions for a higher order mode (TE 136 LH) showing the intricate field distributions being modeled correctly by the phase boundary conditions. The z directed scattered field is zero................................. 77 4.8 Field distributions for a higher order mode (TM45 RH) showing the intricate field distributions being modeled correctly by the phase boundary conditions...................................................... 78 4.9 An example of the erroneous fields generated when applying the phase boundary condition on the E field only. Scattered field due to TEl6 RH excitation is shown. Note large z-directed field, especially near the axis, which should be zero............................................. 80 4.10 Cross section of a small cylindrical guide, terminated by a 3 element thick anisotropic layer. Es for TEB I RH excitation is shown.................. 1 4.11 A cylindrical guide with a fan termination.......................... 33 4.12 Solid model and mesh used to analyze the fan termination................ 84 4.13 Backscatter RCS from fan-terminated guide (4.348 GHz)................ 85 viii... Vill

5.1 Location of hand-off plane where modal scattering matrix is defined....... 89 5.2 Diagram of 8-dimensional sparse matrix indices........................ 93 5.3 Instantaneous rest-frame rotating along with blades..................... 96 5.4 Simplified engine section model................................. 100 5.5 RCS of simplified engine section model for blade position 1 and position 2.. 101 5.6 Comparison of RCS for blade position 2 by mode matching and the technique described in this chapter............................... 102 5.7 Far field amplitude modulation pattern due to rotating blades (90 polarization)................................................... 103 5.8 Far field amplitude modulation pattern due to rotating blades ((p(p polarization)................................................... 104 A. 1 Mesh storage data structure.................................... 112 A.2 Operation used to fill the mesh tree with all of the elements.............. 113 A.3 Operation to efficiently find all free surfaces in the mesh................. 115 A.4 The correct direction for the normal on a free face of an element........... 116 A.5 Operation to divide free surface into unconnected pieces................ 118 A.6 Recursive function to remove connected free elements from the mesh tree.. 119 A.7 Integration contour around a coated conductor with a crack............... 121 ix

LIST OF TABLES Table 2.1 Number of traveling modes for cylindrical guides with different radii....... 17 2.2 Mode indices corresponding to the results in Figure 2.2................. 17 2.3 Impedance of the eight different traveling modes in a cylindrical guide with radius of 0.66k............................................. 29 3.1 Behavior of modes on axis...................................... 53 3.2 Error in the calculated reflection coefficients for an inlet terminated in a short, using a 4 degree slice..................................... 55 3.3 Storage requirement on [S] for lm diameter inlet...................... 58 4.1 Some modal reflection coefficients calculated for the large shorted guide shown in Figure 4.6............................................ 76 4.2 Electrical parameters for the geometry shown in Figure 4.12. Bottom two rows are projected............................................ 86 4.3 Time requirements for each stage of analysis (per frequency) @24 MFlops for the geometry shown in Figure 4.12. Bottom two rows are projected...... 7 4.4 Memory requirements for each stage of analysis for the geometry shown in Figure 4.12. Bottom 2 rows are projected 7...........................:7 x

LIST OF APPENDICES Appendix A. O(N) Graph Algorithms for Mesh Preprocessing............. 110 B. Quasi-Minimum Residual (QMR) Solver........................... 122 xi

CHAPTER I INTRODUCTION 1.1 Motivation Electromagnetic signatures can give clues for the identification and classification of unknown flying targets. Since World War II it has been known that radar echoes from propeller-driven aircraft have a superimposed modulation which make these flying targets distinguishable from other radar reflecting objects, such as storm clouds or birds. However, it wasn't until later that this modulation was exploited as a feature for distinguishing between different aircraft. The effect of rotating jet engines on the radar signatures of aircraft remains an important problem for developing effective automatic target identification capabilities. Another effect of jet engines on the radar signature of aircraft is that they cause a washing out or blurring of high resolution radar images. This blurring is due to the presence of the engine inlets, or ducts which draw air into the front of the engine. Jet engine inlets are resonant cavities, spreading the radar return out over time. Since the very nature of an image is to assume that energy comes from points, by using standard image processing techniques, the spreading out of energy in time due to the engine results in a blurry image. A new technique described in [1] and applied to simple cavities on aircraft I

2 in [2] allows the separation of mechanisms localized in time from mechanisms spread out in time. Naturally, the mechanisms that are spread out in time are localized in frequency, and these frequencies are related in complex ways to the dimensions and materials of the inlet cavity and the face of the jet engine as seen by the radar. Radar signatures that can be measured in practice are extremely complex. Since the target is not stationary with respect to the radar making the measurement, the signature is influenced by the relative motion of the aircraft with respect to the radar. It is indicated in [2], however, that the scattering mechanisms that are localized in frequency due to the inlet cavity and engine face may not be sensitive to the orientation of the aircraft. Additionally, long range radars must have a relatively low (1 KHz) pulse repetition frequency (PRF) to avoid range ambiguity. These long times mean that the aircraft may have changed aspect from pulse to pulse. Jet Engine Modulation (JEM) however, is a fast mechanism since jet engines rotate at approximately 20,000 RPM. JEM therefore occurs in a small time window, on the order of the length of a single radar pulse and the effect has been observed out to sixty degrees off nose-on aspect [3]. Unfortunately, the modulation effect is extremely under-sampled by a radar operating with a low PRF, making the modulation effect appear as random noise. An investigation into the signal processing of such an under-sampled signal is given in [4], where it is shown that since the modulation is periodic, the modulation pattern can be partially reconstructed from the under-sampled data. What is lacking is a robust simulation technique capable of accurately modeling the jet engine cavity and engine face from a radar point of view. The need for such a model

3 provides the motivation for the work within this thesis. 1.2 Jet Engines Jet engines operate on a few simple principles. Their operation is basically the same as that of a piston engine in a car, except that the combustion cycle in a jet is continuous as opposed to the four stroke, or two stroke combustion process in automobile engines. Early jet engines were axial flow engines. This name is derived from the fact that all of the air that enters the engine travels down the main axis of the engine and is used in the combustion process. In an axial flow engine, there are guide vanes which change pitch, acting as a throttle. These guide vanes do not rotate, but are interspersed with rotating compressor blades, which do not change pitch. The series of compressor blades andior guide vanes serves to compress the air many times over. The compressor section may include dozens of compressor blade sets and at least a few guide vanes located in the front section of the compressor series. Once the highly compressed air leaves the compressor section it enters the combustion chamber. Here jet fuel is squirted into the air and is automatically ignited due to the high compression. To start the engine, a spark plug is used to initiate the process. The explosive force of combustion causes hot gas to be propelled back towards a turbine which then rotates and keeps the entire process going. The temperature is the hottest in the turbine section. It is in this section that sophisticated metallurgical processes are needed to manufacture turbine blades which can withstand extremely high heat. Once the hot gas leaves the turbine section, which may have one or more turbines, it is expelled out of the back of the engine. A relatively crude technique (after-burning) for increasing the thrust is

4 to squirt more fuel directly into the exhaust gas. While this technique allowed some aircraft to briefly reach super-sonic speeds, or was used during take off, it tended to waste considerable amounts of fuel. The material used to manufacture axial flow engines was nickel alloy. In modern aircraft, a series of incremental improvements allowed the weight of the compressor section to be reduced by making use of titanium instead of nickel alloy. Nickel alloy is still used in the tortuous environment of the turbine. Of the many incremental improvements made to jet engines over the past half of a century, the development of the turbo fan w as probably the most significant. In a turbo fan engine, a large fan in front of the engine bypasses some of the air around the compressor. The bypassed air is used to cool the engine and decrease the turbulent flow, thus reducing engine noise. Turbo fan engines allowed aircraft to sustain supersonic speeds without the need for fuel wasting after burners. An important parameter used in the specification of a turbo fan engine is its bypass ratio. This number indicates the ratio of bypassed air to the amount of air actually used in combustion. Early turbo fan engines had high bypass ratios while the newest engines are being redesigned with lower bypass ratios as more optimizations are being introduced. All of these improvements have the effect of increasing the overall performance of these amazing, reliable machines. From an electromagnetic point of view, we are only concerned with the part of the engine visible to a radar. From the back, a radar will see the final turbine section only. The turbines have very dense blades and the radar returns are dominated by the scattering froIn

5 the last turbine blade set. Of course, the exhaust flaps (sometimes referred to as "turkey feathers") of the engine are visible to the radar as well. These exhaust flaps are used to redirect the thrust, thus increasing maneuverability, and to reverse the direction of the thrust when the aircraft is braking. From the front, the face of the engine will be visible to the interrogating radar along with the air inlet. This front face may be a titanium turbo fan, with a cone to cover the hub (as in most commercial carriers), or a front frame. A front frame is a stationary set of guide vanes that holds the center shaft fixed in the engine cavity for the rest of the engine to rotate about. An important fact about jet engines from a radar point of view is that they rotate at a fixed speed, usually around 20,000 RPM. This is because the mechanical turbine engine has a natural, optimum speed at which it operates. There is a similar situation in a car where the engine has an optimum speed and a mechanical transmission is used to supply power to the drive shaft at a wide range of speeds. In a jet, however, one method of controlling the amount of air/fuel mixture is by adjusting the pitch of the guide vanes. The change in pitch alters the amount of thrust while keeping the speed of the engine relatively constant. In modern turbo fan engines the fan may be fixed to a separate shaft which allows the fan speed to be slower than the rest of the engine. Additionally, some modern engines have alternating compressor stages which rotate in opposite directions. From an electromagnetic point of view, when the radar is looking at the front of an aircraft, only the turbo fan needs to be modeled since it has a very dense blade set. Alternatively, a combination of a stationary front frame and one compressor section would

6 need to be modeled for engines so equipped. Even with this assumption, we will see that the modeling of aircraft engines from an electromagnetic point of view is a tremendously challenging task. The work presented in this thesis however shows that the problem has at least become feasible for perhaps the first time. 1.3 Available Data and Modeling Techniques Measured data on radar interactions with fan blades is available from only a few sources. Unfortunately, this data often lacks a complete description of the measurement technique used or the geometry of the target and is of little use for reference purposes. Additionally, the targets modeled are often complete engines, or scale models of complete engines. This data cannot be used to validate a modeling technique incrementally. Incremental validation is needed due the complexity of the problem which requires that computational electromagnetic techniques specifically suited to this problem be created. One relatively early source of measured data for a simplified engine model is found in [5]. This work serves as an excellent review of the techniques that were available at the time, or that were being developed for mathematically modeling the engine face. The measurements made in [5] were done to eliminate all of the current candidates for jet engine modeling. For example, the proposal to use a random, statistical function for tracing rays in and out of the duct was discredited with a careful series of measurements showing that the total exposed blade area does not correlate to the amount of power reflected from the blades. A random model is still often proposed today by newcomers to the problem and therefore the measurements contained in [5] remain timely. Also in [5] it was demonstrated that diffraction effects from the blades are essential

7 and that these diffraction effects must include blade interactions. This fact eliminates most high frequency ray optics or physical optics techniques. It is concluded in [5] that only rigorous, integral equation techniques could accurately account for the physics inherent in engine scattering. Given the limited power of the computational resources available at the time and the limited development of exact techniques, the task of doing rigorous analysis on jet engines was considered to be quite daunting. During the 70's, many researchers spent their time developing new mathematically exact models rather than investigating large scale computational efforts. A very good summary of these efforts is contained in [6]. While these mathematically exact techniques were capable of giving good answers to canonical problems, such as simple open-ended waveguides and even a cone-terminated waveguide, they were limited to electrically small problems and could not begin to predict modulation. Still, insights gained into plane wave to mode conversions contained in these works persist as an important contribution to the problem, and some of the work contained in this thesis relies on those advancements. More of the important insights were later distilled into modern works and only these modern works will be referenced in the remainder of this thesis. 1.4 Objective In this thesis, we apply modern finite element methods (FEM) to the jet engine scattering problem with the goal of developing software to be used in generating radar signatures of existing and future aircraft. A tool capable of predicting signatures from aircraft at the design cycle would be a tremendous asset to the aircraft industry, which is currently embracing the technology of rapid computer aided-prototyping.

8 Since full aircraft simulation is the ultimate goal, the FEM analysis must be coupled to a high frequency simulation of the remainder of the aircraft. Thus, the objective is to develop a hybrid finite element modal method where the connection to the high frequency model will be done using a modal scattering matrix.

CHAPTER II FULL THREE DIMENSIONAL HYBRID FEM MODAL IMPLEMENTATION A hybrid finite element modal formulation is presented for the analysis of cavities with complex terminations, with the goal of characterizing jet engine inlets. The finite element method is used to find the generalized scattering matrix for an N-port representation of the complex termination, where N is the number of traveling modes in the cavity. The cavity is assumed to be circular at the termination (engine), but the remainder of the cavity can be of arbitrary cross section. The scattered fields are obtained by tracing the fields back out of the cavity via a high frequency or modal technique with the generalized scattering matrix used in determining the fields at an aperture near the irregular cavity termination. "Proof of concept" results are presented and several issues relating to the implementation of the FEM are addressed. Among these, a new artificial absorber is developed for terminating the FEM mesh and the suitability of edge or nodebased elements is examined. 2.1 Introduction The simulation of radar scattering from jet engine inlets is an important step toward the characterization of aircraft structures. While high frequency techniques can accurately simulate many scattering mechanisms on a typical aircraft, these techniques are 9

10 not suitable for resonant or guiding structures such as antennas/radomes and jet engine inlets. The engine face is an intricate target, possessing complex geometrical features at the wavelength level or less, and is therefore inappropriate for high frequency analysis. By comparison, the finite element method (FEM) [7] [8] is well suited for the analysis of geometrically complex, inhomogeneous volumetric targets such as the engine face. However, in spite of its inherent 0(n) storage demand, the FEM analysis would still require prohibitive computational resources were it to be also used for modeling the volume enclosed by the inlet leading to the engine. To overcome this difficulty, in this chapter we describe a new hybrid finite element method for the analysis of the inlet-engine configuration. The proposed hybrid FEM was originally proposed in [9] and combines ray techniques (for propagating the field to and from the engine face) and the FEM (for computing the fields scattered from the engine face). Figure 2.1 shows the regions to be characterized either by the finite element method or some high frequency technique. Obviously, the high frequency method is best suited for modeling the fields in the large cavity region between the inlet mouth and the engine face. Any of the well known ray methods such as the shooting and bouncing ray (SBR) [10], the generalized ray expansion (GRE) [11], or even a modal decomposition technique can be used for coupling the fields into the inlet region and guiding them from the inlet mouth to the engine face. The same ray or modal technique can also be used for propagating the fields scattered by the engine back out of the inlet. Of importance in this analysis is the interface between the ray/modal and FEM

11 MArn n c lv -I Li FEM [s]/ REGION ENGINE FACE \ CONNECTIVITY BOUNDARY F' RAY OR MODE FIELDS MESH TERMINATION BOUNDARY F' Figure 2.1: Hybrid jet engine inlet analysis. methods, and the truncation of the finite element mesh. Given that the inlet cross section near and at the engine location is circular, and our desire to propose an efficient and flexible coupling scheme, the coupling of the FEM and ray fields in this thesis is accomplished via the generalized modal scattering matrix. That is, the FEM analysis generates the modal scattering matrix which can then be interfaced with any high frequency technique for computing the engine scattered fields without reference to the geometry of the jet engine. Regarding the truncation of the FEM mesh, several schemes are considered including absorbing boundary conditions, the unimoment method (integral equation), and a new artificial broadband absorber, with the latter found most effective for this application. The chapter begins with a section describing the proposed hybrid FEM method, termed as the FEM modal method because it generates the modal scattering matrix. This section describes the role of the modal scattering matrix for interfacing with the FEM and

12 ray techniques used for propagating the fields to and from the engine face. The next sections discuss implementation issues of the FEM technique including elements, application of boundary conditions, and truncation of the FEM mesh. Finally, results are presented for three different inlet terminations. These terminations are rather simple and serve to validate the proposed hybrid FEM modal method since reference calculations using different techniques are available. Given that the emphasis of the thesis is on the development of the FEM-modal technique and not on the inlet field propagator, calculations refer to straight circular inlets. However, the inherent flexibility of the generalized scattering matrix allows the characterization of the same inlet terminations connected to different inlet geometries. 2.2 Finite Element-Modal Formulation Consider the three-dimensional cavity configuration shown in Figure 2.1. The cavity is excited by an arbitrary field (typically a plane wave) through its opening at the left side and is assumed to have a complex geometrical configuration (an engine) at its right end. We are interested in computing the field scattered by this complex cavity termination due to a given excitation which is assumed to be specified at the left opening of the cavity. In our analysis, the cross section of the cavity is assumed to be arbitrary and of diameter greater than a free space wavelength up to the connectivity boundary as shown in Figure 2.1. Beyond this connectivity boundary, the cavity's outer perimeter is assumed to be circular but may enclose complex geometrical configurations such as an engine. In accordance with the proposed hybrid FEM modal formulation, the fields entering the left opening of the inlet will be modeled and propagated up to the

13 connectivity boundary using some ray technique (such as the SBR [10] and the GRE [11] or modal method [12]). These techniques are well understood and this thesis is not concerned with their description and implementation. Instead, our emphasis is on the proposed techniques for modeling the complex cavity termination and coupling of the fields associated with the different techniques at the connectivity boundary. As noted earlier, the FEM will be used to model the fields in the vicinity of the complex cavity termination (i.e. to the right of the connectivity boundary). We could indeed use the modal or ray techniques to generate the excitation to the FEM system of equations. However, this approach would require the solution of the FEM system for each field excitation, a rather inefficient way of characterizing the interior cavity scattering at all incidence angles of the impinging plane wave. Instead, given that the cavity is circular at the connectivity boundary, a convenient way to characterize the termination is by determining its generalized modal scattering matrix. Since each field distribution at the connectivity boundary can be expressed as a sum of incoming (or outgoing) cylindrical waveguide modes, the modal scattering matrix provides us with a unique method for characterizing the cavity termination without consideration of the technique used for modeling the fields to the left or right of the connectivity boundary. The generalized scattering matrix [S] of a given termination relates the coefficients of the incoming modes to the coefficients of the corresponding outgoing modes. That is [S]{a} = {b} (2.1)

14 where the elements of the vectors {a} and {b} are simply the coefficients of the incoming and outgoing modes respectively. They are defined by Je C(x, ) y - " t (x y. zo)ds a = r (2.2).J ( X, Y. Zo) '*m(X, y, Zo)ds -* e (x, y, Zo) T m(x, y, zo)ds bm I - (2.3) IM (x y, z7) -T m (x, y, zo)ds F for m = 1, 2, 3,...N, in which e and e are the incident and scattered transverse electric fields, respectively, on the connectivity boundary surface F, and -m denotes the mode functions of the circular waveguide of a given radius. A single mode index is used to compactly represent the totality of the even/odd, TE or TM modes, each of which is associated with indices n and p. Additionally, each 4T is an odd/even pair, having either sin [np] or cos[n(p] angular dependence. In accordance with the proposed FEM-modal method, the entries of the generalized scattering matrix are computed from the solution of the FEM system. More specifically, the FEM solution proceeds as follows: Given that the cavity accommodates N traveling modes, (a) Use the q'th mode as the excitation

15 (b) Find aq at the connectivity boundary F (c) Solve the finite element system to find the scattered field (d) Calculate the inner product of the scattered field with each out-going traveling mode on the connectivity boundary F to find bp (p = 1, 2, 3,...N) (e) Calculate the q'th column of the scattering matrix as S = Pq aq (f) Repeat for all q = 1, 2, 3,...N It was observed that if the connectivity boundary is placed )/4 from the -2 termination, the coefficients of the outgoing evanescent modes were quite small (<10 ) and were not therefore included in the final calculation of the scattering pattern. We note that the generalized scattering matrix has certain distinct properties: * Its size is N x N, where N is the number of traveling modes (See Table 2.1). * It is symmetric and unitary since all modes are defined to have unit power. It has so far been observed to be sparse since the incoming modes tend to couple more strongly to those outgoing modes havin indices similar to those of the incoming mode. As an example, we illustrate in Figure 2.2 the scattering matrix for a stub terminated inlet as a function of the mode index using the ordering given in Table 2.2. The results shown were generated via the theoretically exact mode-matching technique [13], and the characteristic sparseness of the scattering matrix has been observed for all terminations investigated so far. Given the generalized scattering matrix of the termination, the evaluation of the

16.../ Figure 2.2: Amplitudes of the elements of the generalized scattering matrix calculated by mode-matching for a hub termination. scattered field from the structure is feasible. Since the incident field is known, the coefficients { a} of the incoming modes are also known [15], therefore the coefficients {b} of the outoin modes can be readily evaluated by (2.1). In the context of the modal technique, all outgoing modes are tracked to the open end, where an aperture integration is

17 radius I I l () 0.5 0.66 1.66 3 4 5 10 N 6 10 60 186 328 508 2008 Table 2.1: Number of traveling modes for cylindrical guides with different radii. m Mode.. e -10 TMo 1- TMo l 11-20 TM 1 - TM 10 21-30 TM - TM, 10 31-40 TE, - TE, 10 41-50 TE 1 - TE~, 10 51-60 TE2 - TE; 10 Table 2.2: Mode indices corresponding to the results in Figure 2.2. performed to calculate the radiated field. Since aperture integration is based on the physical optics approximation, it is often necessary to correct the result by including the effect of the rim, i.e. by considering the contribution of the fringe-wave currents [16]. The total far field can be evaluated in closed form for rectangular or cylindrical inlet crosssections [17]. Although there can be higher order scattering mechanisms involving diffraction of the outgoing modes from the mouth which scatter back into the inlet, since the rim effects have been observed to be rather small these higher order mechanisms are not included in the analysis. If the inlet is electrically very wide, one can employ high frequency techniques to model the propagation through the inlet body. The most commonly used techniques are the Shooting and Bouncing Ray method (SBR) [10], and the Generalized Ray Expansion (GRE) [11]. Neither is as accurate as the exact modal method, but they are much simpler

18 from a computational point of view. Moreover, unlike the modal technique, they do not require that the geometry be of canonical shape. In the context of the SBR method, the incident field on the open end is decomposed into a set of parallel ray tubes which are tracked into the cavity. When using the GRE method, the open end aperture is divided into a number of subapertures, and the incident field is decomposed into a set of rays emanating from each subaperture. Unlike with SBR, the rays are not necessarily parallel to each other, and GRE is thus capable of tracing non-planar wavefronts. Regardless of which method is used (SBR or GRE), as soon as the rays reach the connectivity boundary the incoming field is transformed into a superposition of modes using [15] (it is assumed that modes can be defined in the vicinity of the connectivity boundary). The generalized scattering matrix is then computed via the FEM analysis and the amplitudes of the out-going modes are computed from (2.1). The scattered field can be evaluated by means of the reciprocity integral method [18] thus eliminating the need to track the rays back to the open end. 2.2.1 Convergence of the Modal Solution For an inlet with a given radius, we would like to know how many modes for which we will actually have to perform the rigorous hybrid FEM modal analysis, not just the theoretical number of traveling modes that can exist. Given the large number of modes that will be present in realistic size inlets, we wish to know if the solution can be expected to converge before all of the modes have been used. We will investigate the convergence of the modal solution by considering the variation in the far field patterns for a canonical problem for which a body of revolution

19 (BOR) moment method solution exists [14]. By increasing the number of modes used in the body of revolution analysis until convergence is reached, we can determine the highest mode index (n) that will be required in the hybrid FEM modal analysis. That is, each time we add a mode to the body of revolution analysis, it is equivalent, in the context of the hybrid FEM modal technique, to adding another mode group containing the next higher mode index n. In Figure 2.3 we plot the theoretical, highest mode index (n) for the highest order traveling mode as a function of increasing inlet radius. Note that the TE modes always have the highest index. A linear interpolation to the curve (TE) gives a useful engineering formula for the highest index present in an inlet with electrical radius a) nmax- 5.9aX - 1.5 (2.4) Next we consider a representative body of revolution problem whose BOR profile is given in Figure 2.4. The BOR analysis was performed at 2, 4 and 6 GHz where the corresponding electrical radii were 1, 2 and 3 X. By increasing the number of BOR modes up to the maximum value of n as given in Figure 2.3, we observe the convergence of the far field solution with respect to n. We observe in Figure 2.5 that for the 2 GHz case, the solution converges rapidly at the lower observation angles, but converges at the highest observation angles only once all of the modes are used. For the 4 GHz results (Figure 2.6), we see that the solution converges rapidly at the low observation angles, but again almost all of the modes must be used to achieve convergence at the higher observation angles. Except for the theta-theta polarized case, we see that the solution has mostly converged

20 Maximum n of Traveling Mode vs. Radius (n: angular variation of mode) 24.0 ----- TE: --- TM T........................ 2.0.................................................................. 16.0 1' / | E: / /,: I * * *................................................................. 0......... I 4.0. 0.01 0.0 1.0 2.0 3.0 4.0 5.0 Radius (wavelengths) Figure 2.3: Maximum mode index (n) of traveling modes vs. electrical radius. out to about 45 degrees when the n=7 and 8 modes were not included. Of most importance is the 6 GHz result (Figure 2.7) which shows that the solution has converged long before the highest traveling mode index has been included. The solution has in fact converged once the n=1 1 mode is included, even though there are traveling modes up to n=18. This is a promising result which indicates that although there will be a very large number of traveling modes, in practice the higher order traveling modes may be unimportant, thus saving a lot of needless computation. 2.2.2 Finite Element Analysis The traditional FEM analysis involves the solution of the time harmonic weak form of the wave equation in a bounded volumetric region of space. Being a partial

21 30.3175 cm 7.5cm |15 cm BOR - I axis 20 cm Figure 2.4: Profile of BOR cavity used to test convergence of modal solution. differential equation method, the FEM analysis leads to sparse matrices (with about 10 to 50 nonzero entries per row) and permits the modeling of complex inhomogeneous regions without need for special care and considerations. For the problem at hand, the FEM analysis region is shown in Figure 2.1 and is seen to extend a bit beyond the connectivity boundary F to a mesh termination boundary F'. On F' it is necessary to enforce an absorbing boundary condition (ABC) or some other mesh termination scheme which ensures the outgoing nature of the waves. That is, F' must be a non-reflecting boundary. This will be discussed later in more detail. We remark that the cross section between r and r' is again assumed circular for this analysis, but this assumption has no bearing on the actual cavity of Figure 2.1. On the basis of the FEM-modal formulation, we are interested in computing the fields scattered by the cavity termination due to the modal excitation E ic q E (2.5)

22 I phi-phi polarization - N- I2 4.0 _ -.._N-5= | -4.0 -12.0 j i \ -16.0 - - 4.0: theta-theta polarization dBSM -12.0 -\ _ \ / -\, -16.0 - \, -20.0 -20.0r 0.0 10.0 20.0 3 0.0 4 0.0 50.0 6 0.0 70.0 th eta \ theta-theta polarization - _-,~: o _/ ',. _\ -4.0 - - 0.O - 1 C 2 0 4 -12.0 ----- theta Figure 2.5: Convergence of far field with increasing mode index n (2 GHz).

23 10.0 -, { 6.0 -1 2.0 -2.0 ' -' dBSM I -6.0A A -10.0 -j -1B. 0 phi-phi polarization I N ', - - - N 2 ------- % ,%, ---- N =4 ---- N =', N — I- I --- N =7 i - - - - - - 1% =47 1 i i II i i i I 1. i N I1 I Ie., v I I I I I I I 0.0 10.0 20.0 30.0 40.0 50.0 60.0 7 0. 0 theta 24.0 -- 16.0 -0.0 - dBSM -3. 0 A -16.0 -24.0 2 -32.0j theta-theta polarization - N=2 - N= - 1% 4 N .- I - N =6 - N =7 - N=6 I I I I O.C I II I I I i i I. i I 0 1 0.0 20.0 13 0. 0 40.0 50.0 60.0 7 0. 0 theta Figure 2.6: Convergence of far field with increasing, mode index n (4 GHz).

24 24.0 i i phi-phi polarization N — 16.0 ---- N=1 d 32V N=14 -3. -— 'M-', --- —-------------- 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70C.0 theta 16.0 - _ I theta-theta polarization — N — - 24.o --- \, ".!f -32.0 -N I I 0,.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 th eta 8. ht-ht oaithetation dBSM 10. 2. 0. 5._ 6.0 --— ~=~ t"/... N= Figure 2.7: Convergence of far field with increasing mode index n (6 GHz).

25 The total field in the FEM region is then — inc -s E E + (2.6) It is well known that the unknown scattered electric field satisfies the weak wave equation (see for example [7] [8] and [20]) - 2 F' (2.7) F 1 V -JJJq(VX E )(VXT)-kOrT *EjdV -2 f T -r. [nx (VxE f-2 f f (Vx~T )(VxT)-E T where T is referred to as a testing function and must be at least square integrable. For Galerkin implementations, T is set equal to one of the expansion basis functions and for each basis function, a system equation is constructed. However, before proceeding with -s the system construction, it is necessary to first introduce the magnetic field HF = VE and rewrite the above weak equation as J J -(VxES VxT) -T- koE dv - j{jk T - (nx H)ds + lt r i r (2.8) f {I (. )T(vT)dv = -2f T Vvx P - Vk OE cdv in which Qd is the volume occupied by the dielectric not including the fictitious absorber (which will be discussed in the next section). Also, the third term on the left hand side of

26 (2.8) is the penalty term which ensures that the divergence condition be satisfied in the Galerkin least squares sense and is also discussed in the next section. On our way to discretizing the weighted residual equation given above, it is necessary to: (a) tessellate the volume Q2 into smaller elements, (b) choose an expansion for the field in each element, and (c) relate H and E on the outer boundary of Q including F'. These issues are addressed in the next section. 2.2.3 Implementation of the FEM Solution Since the focus of the chapter is to present the hybrid FEM modal analysis, only those aspects of the implementation that are unique to this application will be discussed. Node-based elements were used in this initial implementation because it was found that standard edge elements [20] broke down when the field was purely TE (no z component of the scattered electric field) for large problems. The breakdown of standard edge elements was highly mesh dependent and was manifested in the form of a poorly conditioned global system. In Chapter 4 we will show how the edge elements can be modified to make them more suitable for hybrid FEM modal analysis. Although edge elements are attractive for imposing the boundary conditions on metallic surfaces and corners, we nevertheless resort to standard node elements for this initial implementation because of the aforementioned breakdown. It is well known that node-based elements can have difficulties when applied to electromagnetic problems. The source of these problems has been reported from different points of view [19] [21] [22] [23], and remedies have been proposed. Among them, the

27 penalty-term method has been widely used, but a recent formulation using potentials rather than fields [22] appears promising. However, difficulties again appear with the imposition of boundary conditions on the potentials when dealing with disjointed structures. Since the engine face is a large target with many disconnected components we resorted to the traditional node-based formulation where the penalty term [23] was used for enforcing the divergenceless condition in a least squares sense. Although not rigorous, the traditional node-based implementation [23] gave good results provided care was exercised when enforcing the boundary conditions at metallic boundaries and edges. For a node lying on a conductor, the boundary conditions are st, -t, z1 (2.9) -S - inc t2 E =-t E where t 1. denote the orthonormal unit vectors tangent to the metallic surface. There are many possible ways of coupling the global equations to enforce the boundary conditions (2.9), but there is always a best way which preserves the system's condition. If the boundary conditions are enforced arbitrarily, it is then possible to completely destroy the condition of the system and to generate wrong results. The following procedure will guarantee a well-conditioned final system: Given the metal surface normal n at the node, find the two tangent vectors %t and t? as follows *if y x h > 0. 15 then set 1 = x n Ly '! Yxn\

28 x x n else set tl = I x txxh | X 7| ti Xn *~ t = i X nGiven the three global equations for Ex, E1 and E at the node - find the largest component of 7t (x, y or z) and replace the corresponding global equation with t l E -t i E - find the largest component of 2 (x, y or z) and replace the corresponding gJo-,-S -pinc bal equation with t2 E =-t2 E At the open end of the mesh, the fields radiate into an infinite, cylindrical waveguide. The boundary condition at this open end must absorb all traveling modes in the guide (see Figure 2.1), and the most obvious choice for terminating the FEM mesh is to expand the scattered magnetic field as a sum of all traveling modes and to couple this expansion directly to the FEM equations through the second term on the LHS of (2.8). While this scheme worked well for shorted inlets, it was found to be unstable for the stub terminated inlet. Even if it was stable, the full submatrix resulting from this scheme would make it impractical for large radius cavities. The full submatrix has a size of M x Ar where M is the number of modes in the cavity (including some evanescent modes) and Nr is the number of degrees of freedom on the connectivity boundary F. Both M and Nr increase as (approximately) the square of the radius of the guide. A more efficient scheme is to use a fictitious material absorber designed especially for cylindrical, traveling modes. The use of material absorbers in FEM computation also preserves the sparsity of

29 the overall system and results in a better conditioned system than other mesh termination schemes such as numerical absorbing boundary conditions. The design of a fictitious material absorber to absorb all traveling modes is complicated because the impedance and propagation constants of these modes vary greatly. The impedance of the first eight traveling modes for a guide of radius 0.66?. is given in Table 2.3. These eight modes were used to design an optimized absorber. Clearly, if the absorber was allowed to be very long, it could then be designed to have nearly perfect absorption. However, since the absorber will be part of the FEM mesh, as a design constraint it should be no more than about a half of a free space wavelength long to be practical. Mode TM eo1 TM ll TE~11 TE 21 TM~ I TEe01 TEe 11 TE 21 Z(Q) 307 144 421 557 144 986 421 557 Table 2.3: Impedance of the eight different traveling modes in a cylindrical guide with radius of 0.66X. A Monte Carlo optimization scheme based on [24] which has the desirable properties of finding not only the best performing but also the most stable design was used. A five section, metal backed absorber was used as a starting point with e- and aJr set to unity in each section. The length of each section was set initially to 0.1X. The 15 parameters (material constants and lengths of each section) were varied randomly, independent of one another. Transmission line theory was used to calculate the reflection coefficient for each mode (see Figure 2.8), and the algebraic sum of reflection coefficients

30 (the sum of their absolute values) was used as the global optimization parameter. A pass/ fail criterion on the global parameter was set to give approximately a 50% yield. For each random design, a pass or fail was noted for each parameter, and after a number of designs have been sampled a pattern begins to emerge. Some parameters indicate a tightening up is in order since most of the 'good' designs were centered about a certain value. Other parameters indicate a don't care condition as all values worked equally well. A new set of initial values and ranges was chosen and the process was restarted. Again, the global parameter threshold was chosen to give about a 50% yield. This global threshold continues to drop each iteration until a stable design was found. The final design and its performance is shown in Figure 2.9. Note that the reflection coefficient for all modes is less than 0.1 (-10 dB). The performance of this absorber was tested for a larger guide of radius 22 having 77 traveling modes, and from Figure 2.10 it is seen to have good broad-band modal performance. 2.3 Examples 2.3.1 Example 1 The simplest configuration that could be analyzed is a cylindrical cavity terminated with a flat plate. To terminate the FEM mesh the absorber was placed 0.33X away from the cavity base. We remark that no mode coupling is present in this configuration and that the mesh was generated using the commercial CAD system SDRC I-DEAS with a global element size of 12. For this calculation, the connectivity boundary was placed 0.15; from the base and the scattering matrix calculated using the FEM-modal formulation, which was in turn used to calculate the outgoing mode coefficients. These modes were

31 SH~T 13 '2 '1 - ENGINE TE Ziki zi= z i Zi - ri 1 ri TM Zikzi z. = lv ki = k0 orJrir Figure 2.8: Transmission line model of cylindrical waveguide absorber.

32 0.1 o ~, Io 0.21, 0.18,o l-, I I 1 I =~I - I -c I - I SHORT CYLINDRICAL GUIDE 0r~fi, i, r = 0.95-j1.22 ~rl = 1.33 -jl.5,r2 = 1.04- j0.08 ~rl = 1.15-j0.16 tr3 = 1.25-j0.03 sr3 = 1.34 - j.04.1 IRI.05 - 01 Modal Reflection Coefficients from Optimum Absorber Figure 2.9: Optimum cylindrical waveguide absorber.

33 1.00 - 0.80 - 0.60 iRI 0.40 - -i 0.20 i II I 0.00i 0. 20. 40. 60. 80. MODE INDEX Figure 2.10: Performance of optimum absorber for the 77 traveling modes in a guide with a 2 wavelength radius. then traced back out of the inlet to find the radiated field. The results are shown in Figure 2.11 for a 1 X-long, 1.32X -wide cylindrical cavity. Good agreement is seen as compared to a mode matching solution [13]. Note that the penalty parameter Tc in (2.8) was set to unity for this example. 2.3.2 Example 2 The next simplest cavity termination is a circular stub. This geometry can also be analyzed via a mode matching solution [13] since all of the conductors are on curves of constant coordinates. The absorber was placed 0.33h from the stub and the scattering

34 MONOSTATIC RCS (dx polarization 25.0 -,.I O0 polarization on n I dB/ 5.0 2 -5.0 -i I - i -15.0 - -I -25.0 K 0. dB -12.0 MODE MATCHING -20.0 --- FEMMODAL 1 -20.0 ------------------ 0. 20. 40. 60. 0 (deg.) M ODE MATCHIN( -- FEM/MODAL I I I I 20. 40. 0 (deg.) i I. 60. —....... -1.i32 1, Figure 2.11: Results for example 1: a shorted cylindrical cavity. matrix was computed at a distance of 0.25X in front of the stub. The co-polarized backscatter patterns are given in Figure 2.12. It was found that the proper handling of the boundary conditions at the inner edge (the rim of the stub) was critical. If the total field were set to zero at this inner edge, the

35 MONOSTATIC RCS <to polarization 20.0 20.0 polarization 20.0 - 12.0 - 4.0 - dB/X2 -4.0 -12.0 - -20.0 - I 1 2.0 - - 11-1 \I N \ MODE MATCHING --- FEM:E = 0 on inside edge ------ FEM: E0 = 0 on inside edgei 9 4.0 dB/ X2 -4.0 -12.0 -20.0 -i __ MODE MATCHING --- FEM: E = 0 on inside edge i ------ FEM: E = 0 on inside edgel 9 ~1 I1I ( I I i I I 0. 20. 40. — 4 60. I - I I I I I I 0. 20. 40. 60. 0 (deg.) 0 (deg.) 1 X, inside edge K (.) 01.33X 0.6X 0.33X Figure 2.12: Results for example 2: a stub terminated inlet. results were not as good. However, if only the component tangential to the edge (Ed) was set to zero, whereas E- and E were allowed to float, the results were much improved as is depicted in Figure 2.12.

36 While the results shown are in good agreement for horizontal polarization, there is some discrepancy in the vertical polarization. Finer sampling around the rim does improve the comparison slightly but with very slow convergence. The results shown were obtained using a global element size of k and an element size of Ak around the rim. Again, the penalty term was set to unity. 2.3.3 Example 3 The third example is for an inlet with a ridge termination as shown in Figure 2.13. The termination consists of four ridges (grooves) each with an angular span of 45 degrees. The absorber was placed 0.5X from the ridge interface and the scattering matrix was evaluated at a distance of 0.25X from the discontinuity. For this problem, it was found that the implementation of the penalty term had modest effects on the results. The penalty term contains an arbitrary parameter Tr which in the previous examples was set to unity. However, by using an analytical solution as a reference, or by fitting the formulation to a higher order scheme, it is possible to find an expression for Tl as a function of the element size and local material parameters that minimizes the net error in the FEM solution. This type of approach has been used with success in the field of applied mechanics [25] [26]. For this example we have set (j) V = 842.1 - (2.10) where V is the volume of the element and this value where V is the volume of the element and this value of t1 was determined numerically by

37 MONOSTATIC RCS p$ polarization 99 polarization 20.0 - 20.0 -j -- 1 i i t 12.0I 12.0 -PI dB/: 2 dB/?2 - -4.0 -i -4.0 - -120 -- MODE MATCHING -120 - MODE MATCHING -- FEM/MODAL --- FEM/MODAL -20,0 i-20.0 1 i 0. 20. 40. 60. 0. 20. 40. 60. 0 (deg.) 0 (deg.), o.... j.33.. 0.6. 0.33X side view front view Figure 2.13: Results for example 3: a ridge terminated structure.

38 minimizing the error for the shorted inlet of the same radius. It is noted that for a global mesh size of about 12, the value of T1 is unity. Again, boundary conditions at the internal edges and corners needed to be handled properly in order to achieve physically meaningful results. From example 2 it was found that only the tangential field to the edges should be specified at the internal edges while the other components remained as degrees of freedom. For this example, the additional complexity of the ridge geometry gives rise to corners as well, and it was found that the boundary conditions at these corners could modestly affect the scattering patterns. The boundary conditions imposed at the corners must of course be consistent with the known field behavior and discontinuities [27]. With this in mind, the total field components at the corners formed by the ridges and the outer cylindrical guide were all set to zero. At the inner corners, where the ridge meets the inner cylinder, the geometry is discontinuous in the z direction only, and thus the z directed electric field was a degree of freedom while the other components of the total field were set to zero. Figure 2.13 shows the principal plane, co-polarized RCS patterns computed via the modal-FEM scheme and a mode matching solution [13]. It is seen that the RCS patterns compare adequately with the mode matching solution especially around normal incidence. For this example a global mesh size of 'R was used with a local mesh size of 1 around the discontinuities. It is suspected that the errors in the modal-FEM results are due to the simplistic penalty term and not due to the material absorber (based on the good results achieved for the shorted inlet).

39 2.4 Summary The benchmark tests show the validity of the overall FEM-modal scheme and the utility of the broadband cylindrical mode absorber for terminating the mesh. Also, the results indicate that node-based elements can be used for scattering analysis if boundary conditions are handled properly so as not to disturb the condition of the FEM system. Additional consideration must however be given to the proper implementation of the penalty term by fitting the numerical scheme to known, analytical solutions and selecting a penalty function to minimize the net error. Since the edge elements broke down for the pure TE case, no comparison between edge and node-based elements was possible for this problem. Several different methods and schemes were tried before the above formulation was chosen. A mode matching technique for terminating the mesh was tried but was found to be unstable in addition to its large storage requirements. A mathematical absorbing boundary condition was considered, but has been found to produce poorly conditioned systems in comparison to fictitious material absorbers. Thus, we resorted to a specially designed material absorber for terminating the mesh, and this was shown to have good modal absorption for inlets larger than was initially considered. The validation of the method for non-trivial cavity terminations and the proper implementation of the boundary conditions near edges and discontinuities were indeed challenging tasks. With regards to the validation, a mode-matching code was written which was crucial to the validation of the FEM-modal formulation since it is difficult to isolate the termination scattering by measurement techniques. As seen, the comparison

40 between the results from the two methods were quite good given the geometric complexity and diversity of the solution procedures. However, with the goal of characterizing jet engines at radar frequencies, we clearly have a long way to go. First of all, the electrical size of the problems examined in this chapter are small (radius about one wavelength). These small electrical sizes correspond to engine inlets at about 100-300 MHz. Although these frequencies may be of interest in practice for very long range radars, our primary interest is for microwave frequencies. Secondly, we found that when increasing the electrical size of the problem up to about 2 wavelengths in radius, the nodal based formulation given in this chapter began to yield highly erroneous results, regardless of the choice of penalty parameter. The failure of node-based elements provided the motivation for finding the cause of the breakdown of the edge elements since edge-based elements are described in the literature as being free of the theoretical inconsistencies inherent in nodal elements. Thirdly, the material absorber used here is a resonant design and requires a different mesh for analysis at different frequencies. Also, this design was seen to give measurable errors in the far field, significant enough to warrant the search for an improved mesh truncation scheme. Finally, nothing as yet has been done to address the fact that we must model rotating blades. Solution of the entire problem over and over again for different blade positions will require prohibitive computational expenses.

41 Each of these four issues: (1) enormous electrical size of the problem, (2) the need for edge-based analysis, (3) improved mesh truncation scheme, and (4) the modeling of rotating blades, is addressed in the following chapters and represents the most important contributions contained in this thesis. In the next chapter, the theoretical basis for limited mode coupling is introduced. This limited mode coupling explains the sparse structure of the modal scattering matrices encountered so far. The limited mode coupling also allows the solution of very large problems by considering only one unique slice of the blade termination. The sliced FEM method is given next for small problems, still in terms of node-based elements, but in Chapter 4 a small change to the standard definition of edge elements is given which is found to make them well suited to the new sliced FEM formulation. Also, a broad-band, anisotropic absorber layer is introduced as a candidate for mesh truncation within the sliced FEM problem. In Chapter 5, the most important contribution contained in this work is given. Here we reveal a technique for computing the time-varying return from the rotating engine blades as a simple post-processing step, based on a single solution for the stationary blades.

CHAPTER III OVERLAPPING MODAL AND GEOMETRIC SYMMETRIES By examining the scattering from a cylindrical inlet terminated by a fan-like structure possessing discrete angular symmetry, it is found that only very limited intermodal coupling is possible. The previous chapter considered the solution of the blade structure without any simplifications due to the obvious periodicity of the engine blades. In this chapter, the limited mode coupling from the periodic engine blades is exploited in a node-based hybrid finite element modal scheme. We develop a very efficient solution, in which only one slice of the geometry need be modeled. It is shown that a phase boundary condition at the interior walls of the mesh is sufficient for the complete solution of tilhe problem. The implementation of the phase boundary condition is detailed for the full three-dimensional node-based case, including special considerations for enforcing the boundary conditions along the axis. A simple example is shown for the sliced hybrid finite element modal scheme to validate the method. 3.1 The Need for Drastic Computational Scaling The use of a hybrid finite element modal technique to model the radar scattering from jet engine inlets is depicted in Figure 3.1 [28] [29]. Briefly, the finite element method (FEM) is employed to generate a modal scattering matrix for the engine face while some 42

43 ---- 3D FEM Region - Incident wave is decomposed into modes where modal scattering matrix is delned. Figure 3.1: Illustration of the hybrid finite element-modal analysis of jet engine inlet scattering. generate the scattering matrix. That is, the incoming field is decomposed into waveuirber modes prior to the application of the FEM and the FEM is used only netion boundary where modal scatterin matrix is dened. Fiure 3.1:llri the hybrid finite element m tecniqe as avalydate in Cnhater of inet scattering. himoh frequency or modal technique is used to trace the fields in and out of the inlet. It is necessary to perform the FEM analysis once for each traveling mode in the circular inlet to generate the scatterin mtrix. Tht iarge thructures would involve computational cnto waveghich aruide modes prior to the application of the FEM, and the FEM is used only to henerate the modal scatterinow as matrix. While the hybrid finite element modal technique was validated in Chapter 2 for an engine-like termination consisting of straight blades, the electrical sizes considered were small (approximately I radius), whereas typically the inlet can have a radius of 10 or more. Because the number of degrees of freedom grows as the square of the radius, to apply the method directly to large structures would involve computational costs which are indeed staggering. Also, the number of traveling modes and the size of the scattering matrix grow as the radius squared (note that the analysis must be repeated for each mode).

44 For the inlet configurations considered in Chapter 2 approximately 50,000 elements were needed with about 20,000 degrees of freedom. The analysis was repeated approximately 10 times, once for each mode. Accordingly, an inlet which is 102 in radius would require 100 times the computational resources (5,000,000 elements) and the analysis would need to be repeated 100 times over (2,000,000 degrees of freedom, 1,000 times), thus increasing the total computational cost by about 10,000. In effect, the computational cost increases as the electrical radius to the fourth power. Obviously, some physically derived simplification is needed to scale the problem to a workable size. By exploiting the cyclic geometric symmetry which exists in an engine face, it is shown that the entire problem can be reduced to a single unique slice of the geometry. For example, if the engine face has 40 blades, it is sufficient to model, and carry out the analysis for only a single angular period of the geometry, which encompasses onefortieth of the total computational volume. To achieve this computational scaling, it is necessary to work with modal field excitations and not plane wave excitations. By exploiting the modal (excitation) and geometric symmetries, it can further be shown that a very limited set of scattered modes are possible. This leads to a very sparse scattering matrix (i.e. very few observables are needed for characterizing the angularly periodic engine face) as has been observed in the previous chapter and is mathematically proven in this chapter. It is also demonstrated that all of the scattered modes have a constant phase shift from one slice to another. From a computational point of view, since all scattered modes have equal phase shift across the slice, a phase boundary condition can be imposed at the two interior faces of the FEM mesh to bound the problem. This technique has been used successfully in [30] and [31] and is extended to three-dimensions in this thesis.

45 Section 3.2 of this chapter justifies that the problem can be scaled down to only one slice by showing that the overlapping modal and geometric symmetries give rise to a limited set of scattered/reflected modes. It is shown that these modes can be modeled with a phase boundary condition on the interior walls of the mesh separating the periodic sectors of the engine face. The theory presented is essentially a special extension of Floquet theory to cylindrical structures. Section 3.3 describes the implementation of the phase boundary conditions for the full three-dimensional problem including some important numerical considerations that if overlooked can lead to ill-conditioned systems. Also, enforcement of the boundary conditions along the axis is discussed in light of the fact that the phase boundary conditions cannot be defined there, and thus a different, physically derived set of boundary conditions must be used. Section 3.4 shows results for the hybrid finite element modal analysis making full use of the overlapping modal and geometric symmetries. 3.2 Overlapping Modal and Geometric Symmetries Within a cylindrical waveguide, electromagnetic radiation propagates as modes. These modes are classified as being either TE (transverse electric) or TM (transverse magnetic). Each of these two mode classes are associated with four parameters: the integers n, p, ot, and z. The index n is associated with the angular dependence of the modes, while the index p is associated with the radial dependence of the modes. Although it is often convenient to use cylindrical waveguide modes whose angular dependence is either cos(nO) or sin(no), to exploit the symmetry of the engine face, the modes must be defined as having e dependence. It is the simple angular phase shift of the modes

46 across a sector that makes possible the exploitation of the geometric symmetry [30]. The index cy determines the sense of rotation as the mode propagates down the guide, while the index C7 determines the direction of propagation. Both of these indices take on the values of ~ 1. All modes have a z dependence of e:- where the propagation constant is a function of the mode indices n and p and the radius of the guide. A mode having 6 = +1 corresponds to an incoming (incident) mode while a mode having.- = -1 corresponds to an outgoing (scattered) mode. The modal electric field will be denoted by II p for incident modes and s p for scattered modes, each actually corresponding to a TE/TM pair. Consider a unique slice of an engine-like termination as shown in Figure 3.2. Let Os be the angular extent of the unique slice of the geometry. For any fan-like structure, blades_ llUU -\ symmetry face 2 symmetry face 1 r Figure 3.2: A unique slice of an engine-like termination. ZD z

47 O~ = 2 where N~ is the symmetry number (number of blades). Since the incident field will be a cylindrical mode with an angular dependence of the form e, the FEM system resulting from a solution of the entire problem would take the form I K" K3 Es KN.kEX f ~jlz. /o fe ~2yins ( fe (3.1) ~(Ns - 1)ninOs 'v where Ek denotes the components of the unknown scattered electric field in slice k. Because the geometry is the same in each slice, K = K. = K3 and by linearity, the unknown scattered fields must all be equal to within a phase factor. That is 2 1 +jninos Es = Ese 3 1 +j2n iQs E -1 E e E. = Ee (3.2), 1 ~j(N,- 1 )nis Es = E in and consequently, the scattered field is a periodic function in r with period Os and a progressive phase advance of ej in 'in each period (slice). This restricts the possible scattered field modes, and to find which scattered modes can make up such a field, we

48 consider the modal decomposition of the scattered field. Specifically, on a cross section F of the guide located somewhere within the termination region (say at z = z ) Es (X, y, ) = 7 b, Px yZ0) (3.3),,, = -oo p = 0 where the coefficients b, are found by taking the inner product of the transverse scattered field (es) with each of the outgoing (scattered) modes sY n S -S. e (x, y, zo),,, p(X, y, zo)ds b (3.4),,J, VI p(x, y, 7o) ~ A,,, p((x, y, o) ds F Considering the angular dependence only, the above expansion is seen to reduce to a Fourier transform with respect to 0 since the scattered modes have an eJ ' dependence. That is, the variables o and n7. are a Fourier transform pair. It is thus important to investigate the spectral content of a signal which is periodic in O, but has a progressive phase advance e from one period to the next. In so doing, the set of possible scattered modes will be found. Consider a signal S(4) which is periodic in o but has a progressive phase advance e as illustrated in Figure 3.3. This signal S(O) can be thought of as a sum of A' signals, each with an angular offset and a phase shift from a reference signal SO() which is shown in Figure 3.4.

49 ZO Zejn'n. Zej2n,. Z eJ3n., Zej(N, - 1)nO. /O S(o) 5 --- —-------- I ure 3.3: A periodic sinal with proressive phase advance. Figure 3.3: A periodic signal with progressive phase advance. Z0 i i rI L /0 Li | F (b 110 Figure 3.4: Reference signal used to construct S() ). Mathematically, N,. -1 S(O)= = SO(O-kos)eji O. k= (3.5) and the Fourier transform (o - n ou) of such a signal can be readily found using the out Z —) Z I-nI~~IV II properties of shifting and linearity to give N. - I S(no,,) = SO(no,,), ei k=O ejN, ( in,-.,), 1 o - =|~ j — - "'u) 1S(O) nin n )O.,' SO(n ) = =- e - I (3.6) noul: nin~mNs NsSO(nout) nout = nin+mNs

50 where m is any integer. That is, only if not = ni + mNs m:any integer (3.7) is S(n0,) nonzero. Thus the spectrum of S(O) is made up of a limited set of modes, whose index is given by (3.7). The implication of (3.7) is that for a given Ns (symmetry number or number of blades), only certain scattered modes can exist for a given incident mode, and this of course implies a very sparse scattering matrix. For example, if there are 17 blades and the incident mode is associated with ni = 1. then only scattered modes which have nout = 1,18,-16,35,-33... can exist. For nn = 2 only scattered modes with nout = 2, 19, -15,36,-32... will be present in the scattered field. For an inlet termination with axial symmetry (a body of revolution, such as a circular stub, a cone, a bulb etc.) the symmetry number N, goes to infinity. Consequently, on the basis of the above discussion coupling can only occur to modes which carry the same angular dependence (i.e. nou n= ). This is consistent with what is known from classical body of revolution theory [32]. An interesting observation from the above analysis is that coupling does not occur to modes having n0^ = n~ + 1 since a symmetry number of 1 (N, = 1 ) implies that the termination has no angular symmetry. This fact will turn out to be of great importance for establishing boundary conditions on the axis of the FEM solution. To consider further the

51 impact of the limited mode property on the FEM solution, note that all of the possible scattered modes share a common phase shift from symmetry face 1 to symmetry face 2. Thus, all scattered modes are related, on a cut of constant z, from face 2 to face 1 by 2 ^ ~jninOs Ep = Epe 2 1 +ins (3.8) Eo = Ece 1 ~-J in's E2 Ee E = E e 9 1 where E2 is the scattered field on face 2 and E is the scattered field on face 1. This fact was first exploited in [31] to efficiently compute eigenmodes within a cyclotron using FEM. The implementation of the phase boundary condition for node-based threedimensional FEM analysis of the engine face is discussed next. 3.3 Phase Boundary Conditions for 3-D FEM with Node-Based Elements In the previous section it was shown that the only condition needed to relate the fields on the two symmetry faces is a phase boundary condition (3.8). However, the implementation of this condition is complicated because it must be enforced over a surface which includes the cylinder axis. The phase conditions cannot be defined on the axis and must therefore be replaced. For the initial implementation of the FEM, the mesh which is generated to model the slice must have coincident nodes and elements on the two symmetry faces. In this manner, for each degree of freedom on face 2 (E,EE, E ), there is a corresponding degree

52 of freedom on face 1 (E1,Ey, El ). During assembly of the FEM system, degrees of freedom on face 2 can be discarded in favor of degrees of freedom on face 1 by enforcing the phase boundary condition as,2,, 2 I2 1 (xEx = p (xE El )e.2 2? 1 1nos.. I - - in O ((xEx+yE~) = d) (x:E + yEy)e^ +jE in 4 s E = E1e After some algebra (3.9) becomes 1 2 e 2jn, 1 J1 PX - p e E + PvP - 2 P e+Jn E X ( 2\ 2 O2fxA pxp+ y-Ox — 5 Ey yX - IK <Px kPx) +2 ~jn;2,. 1 E = e E, 1 1 2 2 where P(/ y)' y) and P(x y)' ~(x) are the components of the polar unit vectors at face 1 and 2, respectively. Expression (3.10) can be used directly to assemble degrees of freedom on face 1 in favor of degrees of freedom on face 2. Since the hybrid FEM modal formulation as shown in Figure 3.1 makes use of an

53 absorbing layer to truncate the mesh, the space between the engine face and the absorber must include the axis of the guide. Consequently, some degrees of freedom are located on the axis (since part of the axis is non-metallic) and this introduces a complication in applying the slicing scheme not previously encountered in the literature. While the transverse polar field components cannot be defined on the axis, the transverse cartesian components have specific values. This fact is a motivation for formulating the finite element solution in terms of cartesian rather than polar components. In practice, boundary conditions on the axis must be imposed differently for each modal excitation. By considering the field behavior along the axis and the limited mode coupling effect, a consistent set of boundary conditions can be applied along the axis. First, consider the behavior of the modes on the axis as shown in Table 3.1. MODE n=0 n=1 n>l1 TE Ex, E= 0 ~0 = 0 TE E =0 =0 =0 TM Ex, Ey = 0 0 =0 TM E. 0 =0 =0 Table 3.1: Behavior of modes on axis. Since it was previously noted that mode coupling does not occur to modes having not, = nin ~ 1, there will never be a mode from column 1 and column 2 that exist concurrently. If a mode included in column 1 is present, i.e. in modulo N = 0, then the boundary conditions to be enforced on the axis are Ex = E,. = 0 which holds for all possible scattered modes. If a mode from column 2 is present, i.e. (ni, + 1 ) modulo

54 N = 0 then the boundary condition E_ = 0 is enforced. If all scattered modes are such that nou > 1, then the conditions Ex = E = E = 0 are imposed. Careful implementation of the boundary conditions on conductors that cross the symmetry faces was also necessary to preserve the condition of the system. In this case, conditions must be enforced so that during assembly, both the phase and the metal boundary conditions are enforced without causing deterioration of the system's condition. This can be accomplished using the following procedure which is performed at the element level. For a degree of freedom located on a conductor surface with normal h at symmetry face 1: *if I( x h > 0.15, then set l = else set = '= 2x1 ' n = X nZX t * Given the three local equations for E,E and Es - find the largest component of tl (x, y or z) and replace the corresponding equation in the element system with ^ ~ Es = - l E - find the largest component of t2 (x, y or z) and replace the corresponding -s T-inc equation in the element system with t2 ' E = - t E For a decree of freedom that is located on a conductor at symmetry face 2, the procedure is the same except that the largest component of the tangent vectors (_t and?t,) must be found for the corresponding degree of freedom on face 1 (not face 2 ) and the local

55 equations replaced accordingly. In this manner, when (3.10) is used for assembly, both boundary conditions (on face 1 and face 2) will be enforced correctly and the condition of the system will be preserved. 3.4 Example As an example, an inlet terminated in a short with radius of 0.662. is analyzed by using only a 4 degree slice of the original problem. While this geometry is in fact a body of revolution, it is nevertheless an important benchmark for validating the correct implementation of the phase boundary conditions. The absorber is placed 0.52. from the short and the connectivity boundary (where the scattering matrix is calculated) is located 0.252 from the short. In this case, any size slice is sufficient, but a balance was achieved by considering that if the slice is too thin, more elements will be needed as the element size must shrink to fill the narrow slice. The calculated scattered fields on the slice boundary for two modal excitations are depicted in Figure 3.5 and Figure 3.6 and are seen to have the correct behavior everywhere including the axis. When calculating the scattering matrix, this termination should simply generate a diagonal matrix with each mode having a reflection coefficient of -1. The errors in our calculations using a 4 degree slice are given in Table 3.2 for the first five modes. These errors are within acceptable ranges for finite element implementations and are mainly attributable to the finite reflections from the absorber used to terminate the mesh. Mode TMol TMI l TEoi TEI TE21 % error -0.99 -5.77 +.12 +2.35 -.06 Table 3.2: Error in the calculated reflection coefficients for an inlet terminated in a short, using a 4 degree slice.

56 Figure 3.5: Calculated scattered field Re(Ep) for TM0o excitation of a shorted inlet using a four degree slice. 3.5 Summary In this chapter we introduced the concept of limited mode coupling for simplifying and scaling the hybrid finite element modal solution of jet engine inlet scattering applications. It was shown that angularly periodic inlet/waveguide terminations lead to very sparse scattering matrices and mode-to-mode coupling can be predicted a priori to reduce storage and computing requirements. Most importantly, it was shown that periodic phase boundary conditions can be used to restrict the computational volume to a single periodic cell, as is typically done with ordinary planar periodic antenna arrays. The theory presented here is essentially a special extension of Floquet theory to cylindrical, periodic structures. Since the analysis of jet engine inlet scattering by 3-dimensional FEM requires VLI VCH~ VU YIIVV I~I ULIIJ ~U V JV to

57 Figure 3.6: Calculated scattered field Re (E) for TEB 1 excitation of a shorted inlet using a four degree slice. that electric field degrees of freedom be defined along the axis, a consistent set of physically derived boundary conditions for the fields along the axis was given. The phase boundary condition alone is not sufficient to solve the problem, as the conditions for the fields along the axis must be included. Also, the extension of the periodic phase boundary conditions to a full node-based 3-dimensional FEM solution required the development of practical procedures for enforcing the boundary conditions without destroying the system condition. It was shown that the above analysis could accurately predict the modal fields scattered from a shorted inlet using only a 4 degree slice of the actual cylindrical region. In this chapter, we have addressed one of the four practical issues brought up at the end of Chapter 2, namely the enormous electrical size of the problem. The computational

58 cost associated with the solution of this enormous problem has been scaled by many orders of magnitude in two ways. First, the storage requirements for the sparse, modal scattering matrices as predicted by the limited mode theory presented in this chapter are drastically lower than for the full modal matrices. In Table 3.3 the storage requirements on the modal scattering matrix for a realistic size (Im diameter) inlet with either 20 or 87 blades at 2, 10 and 18 GHz are given. Note that the storage savings realized with the sparse scattering matrices over the full matrices is of 2-3 orders of magnitude for the 10 GHz case. f GHz Modes Blades M Bytes M Bytes (sparse) (full) 2 228 20.01.8 10 5550 20 6.2 500 18 7972 20 12.7 1000 2 228 87.002.8 10 5550 87 1.4 500 18 7972 87 2.9 1000 Table 3.3: Storage requirement on [S] for Im diameter inlet. Another computational cost saver developed in this chapter is the scaling of the FEM solution down to a single slice of the blade termination. The slicing scheme will be more fully developed in the next chapter where we introduce a special edge-based formulation for discrete bodies of revolution. This formulation and the resulting software discussed next are a general contribution to computational electromagnetics even though they are developed here specifically for the jet engine scattering problem.

CHAPTER IV THREE DIMENSIONAL, EDGE-BASED, FINITE ELEMENT ANALYSIS FOR DISCRETE BODIES OF REVOLUTION A simple modification to the standard edge element definitions is found to allow edge elements to be used within the context of the jet engine analysis problem. This modification tends to make the overall FEM system less sensitive to element distortions which always occur in large meshes. An extension to three dimensional, edge-based finite element analysis for modeling electrically large fan-like structures as discrete bodies of revolution is given. By exploiting the overlapping symmetries between a fan-like structure and a modal expansion of the electromagnetic fields, only one lobe of the problem need be solved by the edge-based finite element method without introducing any approximations. This computational scaling makes it possible to solve electrically large structures much more efficiently. A periodic phase boundary condition must be applied to the faces of the mesh describing a single slice of the body, and it is found that the phase boundary condition must be applied to both the electric and magnetic fields for a robust solution with edge-based analysis. Details of the implementation of the phase boundary conditions with edge elements are given along with results to validate the overall technique which involves the concurrent enforcement of phase boundary conditions on both electric and 59

60 magnetic fields. Although developed for the jet engine, the implementation of phase boundary conditions is applicable to any other periodic problem in electromagnetics. 4.1 Introduction to Edge-Based Discrete Body of Revolution Analysis For electrically large objects, brute force application of the Finite Element Method (FEM) for scattering and radiation problems leads to very large systems. However, for some classes of large problems, the domain of the FEM solution can be scaled down by taking advantage of symmetries. Ideally, this scaling should be done without introducing any approximations as is the case for the new technique presented here. Some large problems can be scaled down to workable size without introducing any approximations into the solution but other problems can only be solved approximately. For example, a finite antenna array can be analyzed rigorously element by element with the assumption that each element is a member of an infinite array of elements [33]. However, the infinite array approximation is valid only for elements that do not lie near the edges of the array. Also, scattering from three dimensional Bodies of Revolution (BOR) can be solved by scaling the problem down to a two dimensional one without introducing any approximations into the solution [32]. The power of the exact BOR formulation for electromagnetic scattering computations has allowed modeling of practical vehicle communications problems that would otherwise be prohibitively expensive to simulate due to their large electrical size [34]. In this chapter we extend the application of the edge-based Finite Element Method

61 to Discrete Bodies of Revolution (DBOR). No approximations are introduced when using the DBOR formulation to solve the scattering from fan-like structures within a cylindrical guide. A finite number of discrete modes exist within a cylindrical guide, making the exact solution available by solving the problem on a mode-by-mode basis. For open DBOR structures, such as a fan in free space, the same technique for scaling the problem can be used, but an approximation is introduced since in this case a continuous spread of modes exists. For solution, the continuous spread of modes for the free space case must be made discrete to solve the DBOR problem and thus an approximation is introduced. Resolving the free space fields into a discrete set of modes involves converting a plane wave spectrum into an equivalent set of cylindrical mode functions [35] which can in turn be used to solve the DBOR scattering problem mode by mode. For a fan-like scatterer within a cylindrical guide, the edge-based FEM solution need only be done for one lobe (a unique slice) of the scatterer and then solved for each mode supported by the guide. If the fan-like structure is situated in free space, we can also solve the FEM solution for only one lobe of the scatterer, but we must then use as many modes as are necessary to represent the plane wave excitations of interest. In any event, to scale the domain of the problem down to one lobe of the DBOR structure, we must again introduce a local periodic phase boundary condition on the sides of the mesh used to represent a single slice of the scatterer. Justification for the use of local phase boundary conditions for the analysis of fanlike structures using hybrid FEM modal analysis was given in [36]. The implementation of phase boundary conditions was also discussed in Chapter 3 for node-based elements. In

62 this chapter, phase boundary conditions are implemented with edge-based elements which are known to be more robust than node-based elements and are free of spurious solutions [20] [39]. A fresh perspective on the reasons node-based elements tend to fail for electromagnetics problems while edge-based elements work is given in [37] from a fluid dynamics point of view. In [37], the authors rigorously demonstrate that when the strong statement of the boundary value problem (Maxwell's equations and boundary conditions), is converted to the weak form, the divergence equation from the strong statement becomes an additional boundary condition on the outer edge of the domain. That is, the divergencefree condition which is enforced throughout the domain in the strong statement of the problem is shifted to the boundary of the FEM domain in the weak statement of the problem. Therefore, in the absence of impressed charges, the zero divergence condition must be enforced on the boundary of the FEM domain. This indicates that any element of sufficiently low order (so as to meet the criterion of zero divergence) and having tangential electric field continuity will be appropriate. This result is consistent with first order tetrahedral edge-based elements which have zero divergence, and is consistent with the simplest type of element known to work for electromagnetics, the brick element [38]. However, it appears that only those elements on the boundary of the mesh must enforce the divergence condition. A small change to the standard edge element definition was required in order to avoid the problems discussed in Chapter 2 (pg. 26). Although lengthy analysis went into determining the source of the numerical problems in the edge element definitions, a

63 simple remedy was found by redefining the edge-based scalar degrees of freedom to be Ei = E ei instead of Ei E e, where E is the unknown scattered electric field and ei is the i'th edge unit vector, while e- has magnitude equal to the length of the edge. That is, the edge-based element degrees of freedom are the projection of the field onto the i'th edge, not simply the electric field in the direction of the i'th edge. This redefinition has the effect of reducing the sensitivity of the system's condition number to element distortions. Fairly severe element distortions inevitably occur in large meshes. While this simple remedy was sufficient to allow solution of the large FEM systems presented here, it should be noted that modifications to edge-based approaches are continually evolving to address the fact that large problems often give rise to poorly conditioned systems when edge-based analysis is used [40]. It has been found that the phase boundary condition must be applied to both the electric field and the magnetic field for a robust solution when using edge-based analysis. This fact is in contrast to what was shown in the previous chapter using node-based elements and the available literature for low frequency motor problems [31] and twodimensional edge-based eigenfunction/eigenvalue (resonant field) problems [30]. The need to enforce the phase boundary condition on both the electric and magnetic field with edge elements will be shown algebraically as well as with examples. Enforcement of the phase boundary condition on the magnetic field adds a constraint on the FEM system which must be condensed out at the global level. An efficient procedure for assembling the FEM system with the correct enforcement of the phase boundary condition (on both E and H) will also be given to aid the code developer.

64 This chapter is structured as follows. In Section 4.2, the implementation of phase boundary conditions is given for three-dimensional, edge-based elements on a mesh with internal symmetry (coincident nodes and elements on the faces of the slice). In Section 4.3, results are shown to validate all aspects of the technique. Since our implementation is one of the first to make use of the phase boundary condition for wave applications, validation of the correct implementation of the phase boundary condition is stressed, especially the need to enforce the phase boundary condition on both the electric and magnetic fields. In the summary we indicate how the overall computational scaling yields immense reduction in computational costs over brute force methods. It is stressed that this DBOR method will be able to solve problems which would otherwise be prohibitively expensive to solve. 4.2 Phase Boundary Conditions for Wave Electromagnetic Problems Using Edge-Based Elements A cylindrical guide, terminated by a fan-like structure (see Figure 4.1) is the geometry most amenable to the Discrete Body of Revolution (DBOR) analysis technique presented here. For the type of structure shown in Figure 4.1, the DBOR formulation does not introduce any approximation since the cylindrical guide supports only a finite set of modes. For this work, we will consider only the internal problem of a fan within a cylindrical guide as a means of validating the DBOR formulation. This problem is of practical interest in determining the radar scattering from jet engines. For a cylindrical guide or inlet terminated by a fan-like structure, where each lobe of the structure is itself symmetric, it is possible to generate a mesh for a single lobe having coincident nodes and elements on the two symmetry faces. For example, a guide

65 Figure 4.1: A cylindrical guide terminated by a tan-like structure. terminated in a set of symmetric straight blades can be analyzed by considering only a single lobe and employing a phase boundary condition to connect the fields on the two faces of the slice (see Figure 4.2). For this case, a mesh can be generated having coincident nodes and elements on the two slices by exploiting the internal symmetry within the slice itself. In the case where there is no internal symmetry, such as a set of curved blades, we cannot generate a mesh having coincident nodes and elements on the slice faces (see Figure 4.3). Thus, for curved blade structures, the phase boundary condition must be interpolated from one side to the other. This interpolation is a standard FEM technique commonly used for joining dissimilar mesh areas, or for joining FEM meshes to boundary elements [41]. However here the implementation of the phase boundary condition is made considerably more complex by interpolation. The details of the implementation will be saved for future work. As noted earlier, it has been found that the phase boundary condition must be

66 Figure 4.2: A guide terminated in a set of straight blades (left). For this geometry a symmetric mesh for one lobe (right) can be generated with coincident nodes and elements on face 1 and face 2. face.......... Figure 4.3: A guide terminated in a set of curved blades (left). For this geometry a symmetric mesh for one lobe (right) having coincident nodes and elements on face 1 and face 2 cannot be generated. applied to the manetic as well as the electric field for a robust solution with edge

67 Maxwell's equations (see for example [7]) leads to a global FEM system of the form FK K1 ~EE _ rFoi il 12K ot l (4.1) LK21 K229 lE iI 0o where we have partitioned the FEM system by splitting the degrees of freedom into two groups. In particular, Ei represents the field values within the volume and E ^ is the vector containing the field values on the mesh boundary. Note that since we are using edge based FEM analysis, the sub vectors Eout and Ei, contain scalar projections of the unknown electric field onto the external or internal edges, respectively. Also, the elements of the subvector Tou are given as boundary integrals involving the unknown magnetic field. = -n 'x ( xE)ds = NhxH}ds (4.2) r r J r and rF is the portion of the outer contour of the mesh over an element face containing the edge-based shape function NT. Typically, the Dirichlet (E) or Neumann (H) boundary condition may be imposed on certain portions of the outer boundary (F). A generalized impedance boundary condition [42] may also be used which relates the H field to the E field and the local derivatives of the E field. In each of these cases, enforcement of one and only one of the boundary conditions on any part of the outer boundary (T) is necessary and sufficient to

68 close the domain in a manner that yields a unique solution. However, we will find that for the phase boundary condition it is necessary to enforce a condition on the electric and magnetic fields separately. Consider a unique slice of the DBOR termination (see Figure 4.4). Here the symmetry faces that define the slice can have coincident nodes and elements due to the internal symmetry of the slice. Once each of the edges on the symmetry faces are found, along with a map giving corresponding edge pairs on the two faces, (4. 1) can be written as 11 K12 K1 Eut T out K2 K22 K23 E = 2 (4.3) 21 22 K23 Eout Tout 3K32 K33 Ein 0 where the subvectors E out and Eout contain edge-based scalar electric field unknowns on symmetry face 1 and 2, respectively. Also, the entries in the subvectors T out and p 2out are given by fou = jC)lo n { nxH ds 2J _ (4.4) -, i - 2 i where H and H' denote the magnetic fields on face 1 and 2, respectively. By making use of the overlapping modal and geometric symmetries as discussed in [36], we also know that under single mode excitation, where the incoming mode has an

69. Symmetry face I Symmetry face 2 O\ s Figure 4.4: A unique slice of the FEM domain for a fan-like structure. exponential angular dependence of e, the scattered field must be periodic, slice by jnin,. slice, with the exception of a progressive phase advance e in each slice, where Os is the width of the slice. Using signal analysis, it can be then be shown [36] that the scattered field will consist of a linear combination of modes satisfying the condition nout = ni ~ kNs (4.5) where Nl is the number of blades and k is any integer. From the above discussion 2 1 Jl, io Eolt = Eoute (4.6) and upon inserting (4.6) into (4.3) we get H~~U L~rVI JIIUVI1~1 \ IVj ~~CV I~Z ~ ~

70 (Kll +K12e ')Eout+K13E, = out jn.i,, j 2 (K21 + K22e )Eout + K23Ein = o ( jni., 1 (K31 +K32e )Eout + K33Ein = 0 1 2 We observe that T1 and t, are not decoupled from the unknown electric fields, and it is clear that (4.7) involves four unknowns, but only three equations. Thus, an 1 2 additional constraint is needed to eliminate P,t or Decoupling of the boundary integral terms (involving the unknown scattered magnetic field) does occur for any portion of the outer boundary on a perfect conductor. This includes the portion of the outer boundary along the walls of the guide. If we consider a partitioning of the FEM system as K1 2 K3 EphaseBC TphaseBC KI1 K12 K13 out out K21 K22 K23 EPEC = PEC (4.8) out out K32 K32 K33 E. 0 L J where we have separated the fields on the outer contour into either phase boundary condition terms (unknowns lying on either of the two symmetry faces), or perfect electric conductor (PEC) terms (unknowns lying on metal surfaces). Next, when we enforce the PEC boundary conditions on the PEC terms as Eout = 0 we get

71 KI EhaseBC +KE = phaseBC II+ut K13Ein -out ~.phaseBC + PEC 4.9) K21Eout + K23in = out K phaseBC K31out + K33Ein = 0 Inspection of (4.9) reveals that once the boundary integral terms for the phase boundary portion have been specified, we are left with only two unknowns but three equations. However, the second and third equation are redundant. That is, we only need to use the first and third equations in (4.9) to find the unknown electric fields EhseBC and Ein. This decoupling would have occurred if, instead of setting the fields to zero along the PEC, we had set them equal to an impressed field distribution (as in a scattered field formulation, which is what must be used in practice). Decoupling of the boundary integral terms involving the unknown magnetic field does not occur in (4.7) for the phase boundary condition. We therefore must make a similar argument [36] to relate the magnetic fields on each side of the slice as was done for the electric fields. An additional condition giving the needed relation between the two unknown boundary integral terms in (4.7) can be found physically. To make a physically intuitive argument, we must first introduce a little more notation to keep track of the electric fields and the boundary integral terms in each slice (see Figure 4.5). Along with superscripts 1 and 2, which indicate that the field values lie on the first or second symmetry face defining each slice, respectively, we must also introduce an additional subscript to indicate to which slice the fields belong. Since there is no physical boundary at the interface between two slices not already

72 Slice 2 Slice 1 El, HI 1 1 n2 nl Figure 4.5: Diagram used for physically deriving the phase boundary conditions. taken into account by a traditional type of boundary condition (Dirichlet. Neumann or impedance), the boundary integrals on either side of a slice boundary must exactly cancel. This implies (using the notation shown in Figure 4.5) that I =- -xT2 (4.10) However, since we actually want to relate PTT to T I, we first note that E2 = Ele (4.11) which immediately implies

73 VxE2 = VxE 1 VxE2 = VxEI~ (4.1.2) As a result jH n =O. H2 = Hle (4.13) and accordingly the magnetic field is also periodic, with a progressive phase shift of e in each slice. Thus, the boundary integral terms as defined by (4.4) must satisfy the relation 1 I jein,., '2 = 1e (4.14) and from (4.10) 2 1 jn,,.. = (4.15) which is the missing constraint. Making use of the above result in (4.7) gives Z: V V IV V VIVVI 1 IDl jni,,. 1I (Kl l + K2e )Eo 4 + K13Ei- = out (K21 + K22e O)Elot + K23Ein - jn1 in. I (K31 + K32e )EouE + K33Ein = 0 (4.16) and by eliminating t t we get

74 j K n jns,in.v 2jQin,, ) I j;nO, (K21 + K22e n., +Ke K12e )Eut + (K13e +K23)Ein = 0 (4.17) jniO., (K31 + K32e )Eout + K33Ein = 0 which can be uniquely solved for E anE and E i. Recall that the subscripts on K in (4.17) refer to the partitioning given in (4.3). Although the enforcement of the phase boundary conditions as given by (4.1"7) requires a specialized assembly procedure, the added complexity is nevertheless needed for a robust solution. Formulations based on the enforcement of the phase boundary condition on E only yield erroneous results [30] [31] for some modes and tend to give inaccurate field values near the axis of the cylindrical guide. The formulation based on the global system (4.17) is free of these inconsistencies. Note that the resulting FEM system based on (4.17) will be asymmetric. This asymmetry is a natural consequence of truncating an FEM domain and imposing an asymmetric boundary condition. Of course, the final linear system could always be preconditioned to yield a symmetric system, but the condition of the system would be compromised. In practice, we must retain an asymmetric system and use an appropriate iterative linear solver. For this work, we have had success with the transpose-free QuasiMinimum Residual (QMR) solver which is described in [43] and in Appendix B. Fundamental mesh preprocessing algorithms for FEM analysis are given in Appendix A. In addition to these, specialized algorithms for handling the sliced FEM problem must also be used. These specialized algorithms find the edges on the two symmetry planes, the edges which lie on a conductor, and establish a map between the

75 edges on symmetry plane 2 to symmetry plane 1. This map can then be used in a specialized assembly procedure to enforce the phase boundary conditions as given by (4.17). A detailed description of the specialized mesh preprocessing algorithms along with the assembly procedure for enforcing the phase boundary condition is given in [44]. 4.3 Validation and Results To show the validity of the DBOR formulation, we will first apply the method to an electrically large, shorted, cylindrical guide. In this case, there are no additional approximations in the method and thus we are able to assess the robustness of the edgebased DBOR FEM formulation itself. Most importantly, we do not need to use an approximate boundary condition at the open end of the mesh. We can simply use the exact modal solution for a shorted guide to terminate the mesh. The exact solution is simply a single outgoing mode with a reflection coefficient of -1 or +1 (for TE or TM modes respectively) with the phase referenced to the location of the short. The analysis is done using the same hybrid FEM modal technique given in Chapter 2, except that edge-based elements are used along with the new DBOR formulation for scaling down the size of the FEM domain to a single lobe. The DBOR FEM analysis was applied to an eight degree slice of a shorted guide having an electrical radius of 6 wavelengths (see Figure 4.6). The mesh had a depth of 1/2 free space wavelengths down the guide with the exact boundary condition applied at the open end of the mesh. The modal reflection coefficients were extracted from the mesh on a contour 0.165 free space wavelengths in front of the short by taking a numerical inner Z-n Z- -VC ~V rHV ~IIIV I

76 product. The mesh contained 37808 elements, 7903 nodes and the resulting sparse FEM system based on (4.17) had 48640 degrees of freedom, containing 737626 nonzero entries. The iterative QMR solver (Appendix B) required about 20 to rative solver (Appendix required about 20 to 200 iterations and about 1 to 5 minutes per mode to reach a tolerance of 10 (Euclidean norm of the residual error vector) on an SGI Indigo rated at 24 MFlops..5; o Figure 4.6: The mesh for an 8 degree lobe of a shorted cylindrical guide with an electrical radius of 6 wavelengths. In Table 4.1 we give the reflection coefficients for a representative sampling of the traveling modes in the shorted guide. Since a guide of this electrical size supports 7-32 traveling modes, all of the modes could not be shown. Mode TMol TE11 TM45 TE57 TM63 TE13,6 TM19,2 Reflection 1.0062 -0.9997 0.9926 -0.9998 1.0021 -0.9952 0.99824 Coefficient +j 0.0021 +j 0.0398 -j 0.0017 +j 0.0002 +j 0.0511 +j 0.00011 +j 0.0623 Table 4.1: Some modal reflection coefficients calculated for the large shorted guide shown in Figure 4.6.

77 The traveling modes, including high order traveling modes, are modeled correctly with the DBOR formulation. The intricate scattered modal field distributions are shown in Figure 4.7 and Figure 4.8 for a few of the modes. Note that the modal fields are modeled correctly, even near the axis of the guide. 1.76 1 2.92 -10.C54 -4 O.-73 _: j -o93e8942 Re(Ep).'Re(E -8.08 -4.80I -10.54 -6.77 Figure 4.7: Field distributions for a higher order mode (TE136 LH) showing the intricate field distributions being modeled correctly by the phase boundary conditions. The z directed scattered field is zero. If the phase boundary conditions are applied to the E field only, we find that modal reflection coefficients from the short are incorrect for many modes. This is especially true for the higher order modes. The erroneous reflection coefficients are due to large errors in the fields near the axis of the guide. To show how this error is manifested, we will use the same mesh for the shorted guide (Figure 4.6) but apply the phase boundary conditions to the E field only. The resultant fields are shown in Figure 4.9 and are highly erroneous and ---- -~)~- VIIJ ~ II ~~H II)~~CIIVUU IV IIV 1J~ Il Z give rise to scattering coefficients that do not conserve energy. This same modal field is

cIQ 60 C) ~ C) C C) CD -0 u *1 f.... V 00 A> I I I, I I...... I.......8

79 modeled correctly when the phase boundary condition is applied to both the electric and magnetic fields. To extend the formulation to a fan-like obstacle within a cylindrical guide, we introduce a new boundary condition at the open end of the mesh. Truncation of an FEM mesh for simulating open problems is a topic well represented in the literature, and there is a large selection of techniques for enforcing the out-going nature of fields on a truncated open boundary. Most of the standard techniques, however, have drawbacks when applied to the DBOR FEM problem. For example, the 3-layer absorber considered in Chapter 2 is a resonant design and a new mesh must be generated for different frequencies. A more broad-band, recently introduced technique is to use an anisotropic material absorber, referred to in the literature as a Perfectly Matched Layer (PML) [46]. We have successfully applied the PML absorber to the DBOR FEM formulation by using a metal backed layer of anisotropic material with parameters (4.2 - j4.2) 0 0 l ~= r= 0 (4.2 -j4.2)5Y 0 (4. 0 0 1 (4.2 - j4.2)~. A three element layer of anisotropic material having the parameters given in (4.18) can successfully absorb outgoing modes in the cylindrical guide with very little error introduced. For example, in Figure 4.10 we show the scattered field due to the dominant mode excitation on a cross section of a shorted guide terminated in a 3 element layer of anisotropic material having the parameters given by (4.18). The guide has an electrical radius of 0.33 wavelengths and an electrical length of 3 wavelengths. With this length, we

CD t?14 -.CD....... M 0 N Q 0* CD CD 0: - 0 CD )_____________________________ CD ND 0CD oC-D....... CD.......

81 can assess the accuracy of the anisotropic absorber by observing the standing wave pattern generated in the guide. For the case shown in Figure 4.10, we observe a reflection coefficient of about -50 dB. This reflection coefficient is computed by taking the modal inner product at a "worst case" position in the guide. This is one practical way of defining a reflection coefficient for a non-TEM mode where a voltage cannot be uniquely defined. Figure 4.10: Cross section of a small cylindrical guide, terminated by a 3 element thick anisotropic layer. Es for TE1 RH excitation is shown. Fortunately, the phase boundary condition is unchanged by the presence of the PML. As long as the optic axis of the crystal, or in this case the PML, is along the axial direction, the implementation of the phase boundary condition does not need to be

82 modified. We note that for a non-diagonal anisotropic material or a material whose optic axis is not in the axial direction, the implementation of the phase boundary condition would need to be modified extensively. To demonstrate the compatibility of the PML with the phase boundary conditions, we will model the electrically large fan termination depicted in Figure 4.11. This canonical termination has a reference mode-matching solution [46] which we will use to demonstrate the accuracy of the overall scheme. The FEM mesh used for this analysis is shown in Figure 4.12. Using the procedure discussed in Chapter 2, the modal scattering matrix was computed at a handoff plane 6.25 cm in front of the blades. Backscatter RCS data was then generated by computing the plane wave to modal coupling for a given angle and polarization, propagating the modes down the guide calculating the scattered modes from the handoff plane, propagating the scattered modes back up the guide and computing the modal to plane wave coupling. We note that the plane wave to modal and modal to plane wave coupling is based on physical optics with fringe-wave currents added to account for the effects of the rim [16] with no higher order diffraction effects. A plot of the backscatter RCS patterns at a frequency of 4.348 GHz is shown in Figure 4.13 and a comparison is made with the mode-matching solution [47]. We note that for this model there were 64,050 elements, with 87,020 edge-based degrees of freedom. Although the results were in very good agreement, convergence of the QMR solver was very slow for some modes, requiring as many as 10,000 iterations to reach convergence. As a result, the analysis required over a week of running on an SGI Indigo. The slow convergence can be attributed to the PML layer. We note that although the PML -is

83, *:: Inlet radius: 36.5 cm ',,:' ~ 'is:'i:i "i? iii~i{'~':":~:.~,'i ~,'~::,~ 75 j::j2:-~::~:~:,:~: Inlet length: 50 cm Ei~i:i,:iii~ii.......................... (19.685"):iii~iiiiiiiiiisiii~* Inner hub radius: 23.495 cm (9.25") ~- ~ ~ ~ i...........blad angle. Blade::,iii::::::::!:::s~~!::i. ii:I:~!i.!:ii~i!i::i~'~~:: 2 0 btlades. 2dge........... Inade ne.r h Blradiue: depth 6.985 cm "....(2.7 5") Figure 4.11' A cylindrical guide with a fan termination. incorporate the PML layer and to include a free-space betw een the scatterer and the absorber. Toour m.:Ak theindBORa iEMwt fruation temintorepatclnw.ilrqir pca modralv meshatunctos ee thisacuayatd fullyt expl oit thpesymentryos of the prblmLand ldoe not destroy the condition of the resulting system. 4.4 Summary Substanetial rductons in cormputationa costs hrciave bee realzedqfore thecFinit Element Method (FEM) modelincg of Discrete Bodies of Revolution (DBOR) The key feature of our formulation is the correct implementation of the phase boundary conditions

84,-:'.,........ ' --. blade '-:::.*.7.. ' * f'' ^:'^ ^ /' '"; -: ', -. -.-,... — i -'. 4' ' E Figure 4.12: Solid model and mesh used to analyze the fan termination. and this is one of the first implementations of phase boundary conditions for electromagnetic wave problems. It was found that the phase boundary conditions must be applied to both the electric and the magnetic fields for a robust solution and this was shown mathematically and with three-dimensional examples. The computational scaling made possible by the DBOR formulation reduces the overall size of the FEM system by a factor of NS, where NA is equal to the number of.. blades in the object. Since the required time to solve an FEM system of size N tends to be

85 E..., m.. i".:,.'/ '.. "o...................:................ 5.00 -!K 0.0::'. -5.001. l - -- -15.00.......... -20.00 -60.00 -36.00 -12.00 12.00 36.00 60.00 Theta (deg.) Figure 4.13: Backscatter RCS from fan-terminated guide (4.348 GHz). 2 of O(N-), the solution time per mode is scaled by N. Although the solution must be repeated for each mode, each of these runs can be performed in parallel since they do not require communications. Scattering matrix entries computed by each processor can simply be concatenated to form the complete modal scattering matrix. Scaling a problem for parallel implementation with minimum communications is at the heart of parallel programming, and this formulation therefore represents an advancement in the applicability of parallel computing to electromagnetic FEM analysis. Computational costs for the overall method are shown using three tables. Table 4.2

86 gives the electrical parameters for this geometry at different frequencies (octaves) and indicates the number of modes and size of the mesh. Table 4.3 indicates the time required electrical degrees of Octave freq. (MHz) radius Modes Elements freedom radius freedom 1 274-548.33 - 0.66 2 - 10 5096 7k 2 548-1096.66 - 1.33 10- 36 15902 16k 3 1096-2192 1.33- 2.66 36- 152 35010 45k 4 2192-4348 2.66-5.33 152-578 64050 87k 5 4348-8696 5.33 - 10.66 578 - 2288 125000 200k 6 8696-17392 10.66- 21.33 2288 - 9058 250000 300k J Table 4.2: Electrical parameters for the geometry shown in Figure 4.12. Bottom two rows are projected. per mode, again at different octaves. These times are for an SGI Indigo rated at 24 MFlops. The times are split up into the three phases of analysis: mesh preprocessing, FEM solution and RCS computation. The RCS computation code also fully exploits symmetry by using a sparse modal scattering matrix and by only using those travelling modes that are physically possible when computing the radiated fields. Table 4.4 indicates the storage requirements for each phase of analysis, again for different frequency bands. At this point then, three out of the four practical issues discussed at the end of Chapter 3 have been addressed. The enormous electrical size of the problem has been scaled down as much as possible. Edge-based elements have been introduced successfully and a new broad-band absorber for mesh truncation was shown to be a promising candidate for mesh truncation. The DBOR analysis, along with the exploitation of the sparse modal matrices brings this enormous problem at least into the realm of possibility. The last and most significant issue, that of modeling rotating blades will be discussed in

87 RCS FEM (per (200 angles) Octave preproc. t mode) (per blade angle) 1 10 sec 1 min. 20 sec. 2 30 sec. 5 min. 30 sec. 3 1.5 min. 10 min. 50 sec. 4 5 min. 20 min. 3 min. 5 10 min. 45 min. 10 min. 6 20 min. 180 min. 30 min. Table 4.3: Time requirements for each stage of analysis (per frequency) @24 MFlops for the geometry shown in Figure 4.12. Bottom two rows are projected. Octave preproc. FEM RCS 1 4 MBytes 3 MBytes 1 KByte 2 9 MBytes 8 MBytes 200 KBytes 3 22 MBytes 18 MBytes 1 MBytes 4 45 MBytes 35 MBytes 3 MBytes 5 120 MBytes 100 MBytes 7 MBytes 6 250 MBytes 220 MBytes 13 MBytes Table 4.4: Memory requirements for each stage of analysis for the geometry shown in Figure 4.12. Bottom 2 rows are projected. the next chapter. Since Jet Engine Modulation (JEM) is an important aspect of aircraft target identification, the technique revealed next is probably the most significant contribution in this thesis.

CHAPTER V EFFICIENT COMPUTATION OF JET ENGINE MODULATION By making use of the overlapping symmetries between the modal inlet fields and the engine face, it is found that the radar modulation due to a set of rotating blades can be computed using only one solution for any blade position. Frequency domain solutions for the stationary engine face can immediately produce modulation data with a simple postprocessing step. The validity of the modulation scheme is demonstrated with the aid of a simplified engine model possessing a semi-analytical mode-matching solution. 5.1 Jet Engine Modulation Jet engines are important contributors to the radar scattering cross section (RCS) of aircraft structures, and it is therefore necessary to provide an accurate characterization of their scattering contributions. One of the most important effect the engines have on the signature is that they introduce a modulation into the radar echo. Computation of jet engine modulation is a rather challenging task because of the intricate details associated with the engine blade configuration and the engine's large electrical size (up to sixty wavelengths in diameter) in the frequency range of interest. Taking into consideration the engine face complexity and the length of the inlet Z:llV ZD Z=rlVI L 88

89 (which may span hundreds of wavelengths), a hybrid finite element method modal technique was described in Chapter 2 for the analysis of jet engine inlets. In accordance with this hybrid approach, a ray method such as the SBR [10] or the GRE [11] can be used to propagate the field from the inlet opening to the engine face and from there on the FEM is employed for modeling the intricate compressor and engine blades. The interface of the two computational regions is accomplished using the generalized modal scattering matrix, and this relies on the fact that the inner cross section of the inlet is circular at the engine location. Basically, at some handoff surface (see Figure 5.1), the incoming ray field is decomposed into cylindrical modes and each of these modes is used as the excitation to the FEM computational section. The FEM or some other numerical technique must be used to return all mode scattering coefficients associated with a given incoming mode. By repeating the process for each incoming mode, all entries of the modal scattering matrix are computed. hand-off plane Figure 5.1: Location of hand-off plane where modal scattering matrix is defined.

90 The scattering matrix is unique to each engine configuration, and as a result it does not need to be recomputed for each excitation. Instead, for different inlets and excitations, it only needs to be weighted with the mode coefficients of the incoming wave at the handoff surface. In this manner the computational burden is substantially reduced. Also, it was pointed out in Chapter 4 that a major portion of the scattering is null due to mode decoupling (a consequence of the engine's circular periodicity) which can be predicted a priori. As a result storage requirements for the scattering matrix are drastically reduced and additional computational savings can be realized. However, even with the above domain decompositions and computational improvements, the CPU requirements for a characterization of the jet engine modulation is still daunting. It is quite impractical to recompute the scattering pattern of the jet engine using a rigorous numerical method at each blade position in order to synthesize the modulation pattern. The subject of this chapter is to show that the modal scattering matrix of the engine at different rotation positions can be generated by modifying the available scattering matrix from a previous engine position. That is, the modulation pattern of a rotating engine can be computed from the numerical data available for a single stationary position of the engine. This chapter reviews the structure and characteristics of the modal scattering matrix. By considering the blade periodicity in conjunction with the structure of the modal scattering matrix, it is shown that the scattering matrices at other engine positions can be obtained as a postprocessing step without the need to repeat the computationally intensive full-wave electromagnetic analysis. Considerations relating to relativistic effects and

91 engine rotation frequency versus radar frequencies are also discussed to support the validity of the model. Finally, numerical results are presented, and rigorous validation of the technique is made via comparison with results obtained from simulations which do not ma: —' use of the proposed simplification. 5.2 Modal Scattering Matrix Since approximate methods, such as physical optics, geometric optics. or shooting and bouncing rays, are inappropriate for characterizing an aircraft engine face, more rigorous techniques such as the mode-matching method or the finite element method must be used. The solution for the scattering from the engine face must be connected to the solution for the remainder of the aircraft. By defining the modal scattering matrix [S] at some hand-off plane sufficiently far in front of the engine face where the inlet is still circular, only propagating modes need be considered (see Figure 5.1). Multiple engine sections can be modeled by cascading scattering matrices together for each separate section of the engine. In this way, very complex modulation patterns can be computed for engine faces having multiple sections rotating at different speeds. This complex modulation can can still be computed as a simple post-processing step using the techniques outlined here. Mathematically, the modal scattering matrix [S] is defined by bi s1j = (5.1) a1 where

92 (X, ) ' t (X, y, Zo)ds J(x y o) J (x y, zo)ds a = r b. - r (5.2). * JS* Ji(x, y, Zo) ' ti (, y, zo)d, y, zo) ' (x, y, zo)ds r r inc s for i, j = 1, 2, 3,...N, in which e and e are the incident and scattered transverse electric fields, respectively, at the hand off plane. In this chapter we do not discuss how the elements of the scattering matrix are generated, but rather how once found they can be used to compute the radar modulation based on one solution for the stationary engine. We simply remark that the elements of the scattering matrix Si can be found rigorously for canonical terminations by the mode matching method [47] and for more realistic structures using the hybrid finite element modal method (Chapter 2), provided that thate discrete angular symmetry of the geometry is exploited to scale the problem down to a workable size (Chapter 4). It was proven in Chapter 3 that the modal scattering [S] matrix has a very sparse structure which can be determined in advance. It is fortunate that the modal scattering matrix is very sparse since many thousands of travelling modes can exist in a real jet engine inlet. Consider for example that at X-band, a typical engine on a commercial carrier may have a radius of about 30?X which would allow 17,878 distinct travelling modes. Only by exploiting the sparsity of the modal scattering matrix can we manage the volume of data that results from multiple frequency runs. The sparse modal scattering matrix that represents each section of the engine is in

93 fact an 8-dimensional sparse matrix whose indices can be conceptualized using the diagram shown in Figure 5.2. Each outgoing mode is either TE or TM, indexed by p (radial variation), and n (angular variation), and with either sin(n(p) (odd) or cos(n(p) (even) (with respect to the radial directed electric field) angular dependence. The sparse structure of the scattering matrix is a result of the fact that the relation (3.7) must hold. I 1 TE even TM even 3 X 2 I even p 1 2 3. TE TM even Incoming Mode Index Outgoing Mode Index Figure 5.2: Diagram of 8-dimensional sparse matrix indices.

94 5.3 Efficient Computation of Modulation Once found, the sparse modal scattering matrix for the stationary blade set can be used for computing the radar modulation due to the rotating blades. If we assume that the blades are turning sufficiently slowly so that: (1) relativistic effects are unimportant, (2) effects due to free charges are unimportant, and (3) the frequency of revolution of the blades is much less than the radar frequency, then we can make use of an instantaneous rest-frame that rotates along with the blades. With an instantaneous rest-frame we assume that the steady-state solution is valid at each incremental rotation of the instantaneous restframe, and that the modulation can be gleaned by computing the far fields for each incremental blade rotation until one period has been found. The three assumptions made above require some justification and can be seen to be valid once we introduce some realistic numbers. First, the assumption that relativistic effects are not present is confirmed in [48] and [49]. It was noted in [48] and [49] that relativistic effects are observed only if the blade tip travels at speeds above 10- c where c is the velocity of light in free space. A conservative estimate of realistic blade tip velocity -4 is 10 c if we assume a 1 meter radius fan turning at 100,000 RPM. In short, no relativistic effects are required to predict the modulation from jet engines. Second, the assumption that no effects due to free charges exist must be explored. Normally, free charge within conductors are only taken into consideration for very fast pulses (nanosecond pulses) through lossy materials where the relaxation time is of the same order as the pulse duration. Conversely, if the engine were to spin fast enough, the blade tips would be subjected to an excitation by an incident field varying so rapidly that the free

95 charges would not be able to "keep up." To see if this is in fact the case, assume conservatively that the relaxation time for the material in the engine blades (titanium alloys) has a relaxation time of 10- sec. Therefore, the boundary conditions at a given -9 location on the blade tip would also have to change on the order of 10 sec as the blade rotates through the incident field. Now consider the velocity required for a point on the tip of the blade to see a 180~ phase change in the boundary conditions for a 30 GHz -9 (. = 1cm) incident plane wave. The blade tip velocity would need to be 0.5cm/109 sec or 0.5 x 10-7m/sec = 10-2c. But our conservative estimate of linear blade tip velocity was only 104c. Thus, no free charge effects are needed to compute the radar modulation. The third assumption stating that the frequency of revolution is much less than the radar frequency is most important and also more difficult to justify. Using the same conservative estimate for blade tip velocity based on a im radius and 100,000 RPM, consider that at 10 GHz the blades rotate only 105 degrees during one cycle of the interrogating radar. Thus, even after 100 cycles, the blade would have only rotated through 103 degrees. Indeed we have observed experimentally [50] that although small angular — 3 displacements (0.1~- 1~) of blades within a circular cavity can produce drastic changes in the scattering pattern, angular displacements on the order of only 10-3 degrees do not produce any observable changes. That is, very small rotations of the blades occurring within the 100 cycles of a 10 GHz radar have always produced indistinguishable modal scattering matrices, and therefore identical far field patterns. Finally, the condition that the

96 frequency of the blade rotations must be much lower than the radar frequency also insures that there is no mode which is "slow" enough that the blades would have moved appreciably during the time required for the mode to propagate in and out of them. The speed of the "slowest" mode in the guide can be determined experimentally by observing the ringing in a time domain response of an inlet with a bladed termination. It has been observed [50] that the ringing decays within 10 ns from the initial return at the front of the -2 blades. Therefore the slowest mode will see less than a 10 degree change in the angular position of the blades during a round trip through the blade section. Again this angular displacement is too small to produce observable differences in the scattered field. Based on the above, we can algebraically update the scattering matrix for a series of incremental angular displacements (to account for the rotation of the blades) using data extracted from the instantaneous rest-frame analysis. That is, once the scattering matrix is found for position 1 (see Figure 5.3), the scattering matrix for position 2, rotated by A(p, Yrt-fe Yrest-frame - rest-frame Xrest-frame Position 1 Position 2 (b} = [S]{a} {bA } = [SA]{a} Figure 5.3: Instantaneous rest-frame rotating along with blades.

97 can be found without recomputing the complex scattering problem. However, to do this, we must first convert the modal scattering matrix [S] to an equivalent matrix [S'] for modes having exponential angular dependence e-J2(P. To convert the scattering matrix vwe define two new mode functions with either e n or e -2c angular dependence as t+ = e even + jd dd =y- = even _ j.dd The elements of the new modal scattering matrix [S'] can be found in terms of elements of [S] by considering both systems to be made up of 2 x 2 subsystems of the form beven) Lc dJ a } bodd (5.4) le fa} I F" Lg h a-J bIn these subsystems, a refer to the coefficients of the modes e v en, odd, whereas a+' refer to the corresponding coefficients of the modes P+'-. From (5.3), it is easily found that (a +d) (b-c) (a-d) (b + c) 2 2j 2 (2j) (5.5) (a - d) (b+c) h (a+d) (b-c) - + h2j - 2 2j 2 2j To account for the rotation of the blades we simply transform the original problem

98 into the instantaneous rest frame which rotates along with the blades. The updated system for blade position 2 can then be found by introducing a simple phase shift e ' into the incoming modes to simulate the blade rotation, and then an opposite phase shift -jnz(Azp) e ) of the outgoing modes to rotate the outgoing modes back to the original reference frame which is stationary. The corresponding scattering matrix of the rotated blades is obtained by modifying each submatrix of [S] as (e)ej(n - n)AO (g)e J(nj + ni;)Ao jf j(-ni-ni)Ao e ' (f)e ieA fA (h )ej(-ni' A + hA) (5.6) The submatrices of the odd/even subsystem are now aA bA CA dAj (5.7) where (eA + fA + gA + hA) a, = 2 (e + f - g-h) CAo= - 2j b = (e- fA + gA- hA) 'A 2 d = (e -fA-g + hA) dA -2 (5.8) The updated scattering matrix can now be used to compute the outgoing mode coefficients for the blades in position 2 using the same set of incoming modal coefficients. For a set of blades rotating with angular velocity co we have A(p = co(t - to). where lo corresponds to the time when the blades are at position 1. By computing the far fields

99 based on the updated outgoing modal coefficients for a sufficiently large number of times t, the modulated radar return is found. This modulation must be periodic with a period of mod = ) (59) The radar return will be both amplitude and phase modulated and can be analyzed in a variety of ways, but not simply with Doppler analysis which only makes use of the phase modulation. For example the product Ns ) is inversely related to the fundamental period of the amplitude modulation. We will only consider the efficient determination of the complex modulation pattern itself and leave the signal analysis for a future work. 5.4 Validation To validate the modulation scheme presented above, we first see if the exact solution for two different blade positions can be recovered from the scattering matrix at a single position using the technique presented here. A simplified engine section was constructed (see Figure 5.4) and analyzed via the semi-analytical mode matching method [47] as a reference solution. The implementation of this mode matching solution has been experimentally validated in [50] and has been shown to be accurate for this particular model over a frequency range of 2 - 12GHz. The 6 GHz, monostatic RCS patterns, for two different blade positions are shown in Figure 5.5 with the cone removed. The blades are aligned with the measurement plane for position 1 and rotated 20~ counter-clockwise for position 2 (as shown in Figure 5.3 with A(p = 20~). When the scattering matrix for position 1 is modified using (5.4)-(5.8) with Acp = 20~ to account for the blade rotation,

100 we see that the scattered field pattern is essentially identical to the mode matching solution for position 2 (see Figure 5.6). In fact, the scattering matrices obtained using this method are nearly identical term by term with those generated using the direct analytical method but the additional computational cost for updating the scattering matrix using this new technique is insignificant as compared to the cost of redoing the analysis. Although not shown, the other polarization can be recovered in the same manner. *Overall length of inlet: 30 cm. * Inlet open at front, closed at back. * Length of blades: 10 cm. * Inner Diameter of inlet: 30 cm. Diameter of inner hub: 15 cm. Number of blades: 8..*..Blades are straight and 40 wide........ Cone interior angle: 60 * Blades are in contact with inlet * All components made using SLA and are painted copper amplitude modulation pattern f.r two different azimuth angles and for both polarizations is shown in Figure.7 ad F. The moulation is plotted iasa function of the blade mangle ather han as a funcltion of tiern to srepes ated e ofaular velocit of the bldes angle;p rather than as a function of time to stress that the angular velocity of the blades is b" `~ L~~Clr~CL VIVIL~JVI ~1V VCC~ZU n

101 RCS of inlet with straight blades: 6 GHz, (06) 10 Blade position 1.;;.k**'.::- '.".. "::'......... X."~~'~~ j ~ I-~,, 0 -5 1,& Blade position 2 dBSM -10 - -15 - -20 - -25 - I I I I -30 0 10 20 0 40 50 60 70 0 (degrees) Figure 5.5: RCS of simplified engine section model for blade position 1 and position 2. tzn Z: -. ~Y unimportant for the modulation pattern. With different angular velocities, the modulation pattern would be the same, only the period would change as given by (5.9). Note that both the overall backscatter level and the dynamic range of the modulation are higher for the 00 co-polarized case than for the (p co-polarized case. 5.5 Summary Computing the scattering from jet engines is a difficult and expensive task. The scheme presented here allows the radar modulation pattern to be found as an inexpensive post-processing step based on a single solution for stationary blade sets. Doing rigorous

102 RCS of inlet with straight blades: 6 GHz, (00) *.... 10 5 o 0 dBSM This method, Mode Matching - %. kt,:a OF 1, A... -5 F 10 F -1 5 I 1. I I I 0 10 20 30 40 50 60 70 0 (degrees) Figure 5.6: Comparison of RCS for blade position 2 by mode matching and the technique described in this chapter. analysis for different blade positions is prohibitively expensive in practice and therefore the technique presented here will be essential for carrying out scattering simulations of jet engines including modulation. Only a single, sparse modal scattering matrix is required per frequency and scattering mating rices for each engine section can be modified term by term and cascaded to simulate the different rotating or stationary engine sections. t ~rV Z:VIWI V1

103 RCS amplitude modulation of inlet with straight blades: 6 GHz, (00),,,,,,! 5 Azimuth(0): 20~ 30~ 45~ <'.A -10 / I dBSM - 15 l / 11,i' I 11 1~r -30 F e.:.1.....11 " —o -35 - I- -- I - I I --- I -- I I 5 10 15 20 25 A(p (degrees) 30 35 40 45 Figure 5.7: Far field amplitude polarization). modulation pattern due to rotating blades (00

104 RCS amplitude modulation of inlet with straight blades: 6 GHz, ((pp) Azimuth(O): 20~ 30 45 - B-. /. 1 0-. dBSM\ '/ 1 -15 -25 -30 0 5 10 15 20 25 30 35 40 45 A(p (degrees) Figure 5.8: Far field amplitude modulation pattern due to rotating blades ((p4 polarization).

CHAPTER VI SUMMARY AND SUGGESTIONS FOR FUTURE WORK 6.1 Summary In this thesis a hybrid finite element modal solution for the jet engine scattering problem was given in a manner such that its interface with a high frequency or other ray optics solution for the remainder of the aircraft can be accomplished. The connection between the FEM solution and the high frequency method is done using a modal scattering matrix to represent the engine termination within a cylindrical inlet. Therefore, the goal of this work was to create a tool capable of generating the modal scattering matrices to replace the engine face in a high frequency model. Initially, a 3-dimensional hybrid finite element modal solution was implemented using standard node-based elements and the far field scattering patterns were compared to an exact mode-matching solution as a reference. Only small electrical sizes could be considered, corresponding to real jet engines at about 100-300 MHz. Since the goal was to create a tool for simulating engines at radar frequencies, drastic reductions in computing costs were sought. To this end, a theory of overlapping modal and geometric symmetries was given, uncovering the mechanism by which limited mode coupling occurs. This theory is 105

106 essentially a special extension of Floquet theory to angularly periodic structures. Limited mode coupling gives rise to very sparse modal matrices for each engine section, thus reducing the computer storage space required for the modal scattering matrices by at least 2 orders of magnitude. A full three dimensional discrete body of revolution formulation was given which takes full advantage of the symmetry to reduce the problem down to a single slice of the termination. A small change in the standard edge-based element definition was made which allowed edge elements to be used successfully in this context. The discrete body of revolution formulation allowed much larger electrical sizes to be modeled, and with the use of edge-based elements, robust, accurate FEM solutions result. Although the computational costs are still high, the problem has at least been made feasible. The new formulation and the resulting software are considered to be a significant contribution to computational electromagnetics which grew out of this research. The single most important contribution made in this work is the technique presented in Chapter 5 which allows the modulation to be computed as a post-processing step. By using only a single solution for the stationary engine blade, the modulated return due to the rotating engine can be found immediately without the need to recompute the scattering for every blade angle. 6.2 Suggestions for Future Work There are two clear engineering tasks that could be done to make the analysis software a more practical tool for real engine simulation.

107 6.2.1 Task 1 In practice, the mesh that will be used to model a single lobe (unique slice) of the fan-like scatterer will not generally have coincident nodes and elements on the faces of the slice. Since the phase boundary condition relates the fields on the two slices, an interpolation scheme must be introduced so that the phase boundary condition can be imposed on meshes which do not have coincident nodes and elements on the two symmetry faces defining the slice. This interpolation technique will be of practical importance since symmetric meshes cannot be generated for realistic structures except for the special case where each lobe of the structure is itself symmetric about the mid-line of the slice. In this thesis, we considered only those problems where a symmetric mesh can be generated, yielding coincident nodes and elements on the two faces of the slice. While extension to unsymmetrical meshes involves application of standard interpolation techniques (conventionally used for joining dissimilar mesh areas), due to the complexity of the phase boundary condition this interpolation scheme will be more involved. When simply joining dissimilar mesh areas together, only continuity of the tangential electric field need be enforced, and in fact the magnetic field must not be made continuous or the FEM solution will be overconstrained. The implementation is not as straightforward for the phase boundary conditions as for interpolation between dissimilar mesh areas since both the electric and magnetic fields must be interpolated from face to face. This is a clear engineering task however, and would require modification of the assembly operator only. An additional data structure would be needed to hold the map from edges on symmetry face 2 to element faces on symmetry face 1 and an additional 0(N) pre-processing step would be needed to fill this structure.

108 6.2.2 Task 2 The use of a Perfectly Matched Layer (PML) for truncating the open end of the mesh was discussed and results were shown to demonstrate the compatibility of the PML with the phase boundary conditions. However, the PML yielded poorly conditioned FEIM systems. An additional and final use of symmetry must therefore be exploited to create a special, sparse modal integral mesh truncation scheme. It appears that the source of the numerical instability found in the modal, integral mesh truncation scheme described in Chapter 2 (pg. 28) was due to the inclusion of physically impossible modes and that the integral equation was enforced over redundant sections of the open end. If so, this would be a new manifestation of the relative convergence phenomenon. By including only those modes which satisfy the periodicity conditions as described in Chapter 3, and by imposing the integral equation on a single slice of the open end, the domain of the integral equation will be profoundly scaled and the formulation will make the fullest exploitation of symmetry. Also, in contrast with the PML, a sparse modal boundary integral mesh truncation scheme will not require an additional volume added to the mesh. With the sparse modal mesh truncation scheme, the DBOR FEM formulation presented here will allow practical simulation of jet engines and open up an entirely new class of problems for which rigorous, robust electromagnetic analysis will be available.

APPENDICES 109

110 APPENDIX A 0(N) Graph Algorithms for Mesh Preprocessing Electromagnetic analysis with the finite element method can be thought of as a four step process: mesh generation, preprocessing, analysis, and postprocessing. The preprocessing step is often incorporated into the analysis code or considered to be part of the mesh generation step. For the large, complex objects (especially in three dimensions) that arise naturally in electromagnetic analysis, this task can be very cumbersome and should be considered as a separate process. Preprocessing involves the extraction of information which is required for imposing boundary conditions on metal surfaces and mesh truncation surfaces where some absorbing boundary condition is enforced, as well as material discontinuities and outward normals to these surfaces. For field calculations outside the region of the mesh an integration contour or surface must also be found which completely encloses the object while passing through element centers where field derivatives can be evaluated accurately for calculating the scattered field. Without an accurate evaluation of the field derivatives on such a surface, substantial error (3dB) in far field calculations and even more inaccuracy in near field calculations, is likely to result. Although all of the above information could be specified during mesh generation, this would in general involve a

111 great deal of unnecessary work. The goal of the preprocessor should be to minimize the amount of interaction with the user. The following algorithms have been developed for preprocess a fin inite element mesh for electromagnetic analysis. The mesh is only required to contain a minimum amount of geometric information. The algorithms will extract from the mesh metal boundaries (any number of metal objects), outer boundary, material boundaries, outward pointing normals, element edges lying on conductors (for edge based analysis) and an integration contour passing through element centers. Although the algorithms are appropriate for many different elements in two and three space dimensions. the code fragments shown are for the particular case of tetrahedral elements. The code is written in C since recursion and dynamic memory allocation are needed. Data Structures The mesh is specified by filling the simple data structures shown in Figure A..1 with the data from the mesh file which contains the node locations and element specifications as generated by the mesh generator. This is enough information to find metal surfaces, the outer surface of the mesh, material discontinuities, outward pointing surface normals, element edges lying on metal surfaces and an integration contour that encloses the target but passes through element centers. Nodal grouping is required only to identify infinitely thin conductors or edges and corners that are to be treated differently by the finite element code (singular elements). No other information is needed from the user.

112 int nnodes; int nelements; struct node { float x,y,z; } *nodes; struct element { int nodes [4 ], material; } *elements; struct material{ float ReE,ImE,ReM,ImM; } *materials; Number of nodes in mesh Number of elements in mesh Array of node coordinates Array of tetrahedral elements (4 nodes and a material index number) Array of materials (Complex ~ and g) Figure A. 1: Mesh storage data structure. Once the mesh data structures are filled, the mesh tree must be created. The mesh tree is a two dimensional array that stores the element numbers of elements which share a common node. The structure is cheap to fill and allows for the efficient extraction of all the necessary information. The operation shown in Figure A.2 fills the mesh tree with all of the elements in the mesh. Operations *Find free surface This algorithm operates on the mesh tree filled with all elements. Let nface be the number of coplanar nodes defining an element face. (nface = 3 for a tetrahedral element.) The algorithm shown in Figure A.3 efficiently scans the mesh tree for groups of

113 C 0 0 #define MAXCONNECTIONS 50 int **meshtree int n,m,nl,i; 0 0 Maximum number of edges sharing a node. Two dimensional array to store mesh tree. Some integers. 0 meshtree= (int **) calloc (nnodes+l, Dynamically allocate memory for sizeof ( *meshtree)) the mesh tree and initialize to zero. for (n=0;n<=nvertices+l;n++) mesh tree[n]=(int *) caloc( MAXCONNEC - TIONS,sizeof(**mesh-tree)); 0 0 0 for (n=l;n<=nelements;n++) for (m=0;m<4;m++) nl=elements[n].vertices[m]; i=0; while(mesh_tree[nl] [i]!=0) i=i+l; mesh_tree[nl][i]=n; mesh_tree[nl] [i+l]=0; } Fill mesh tree with all elements. Figure A.2: Operation used to fill the mesh tree with all of the elements. nface coplanar nodes (faces) shared by more than one element. A face which is not shared by more than one element is a free face. Upon completion of this operation there is one free surface. This will be divided into pieces (outer surface, metal surfaces) by filling the mesh tree with the free elements only and doing a divide free surface operation. *Find material discontinuities and normals Inside the find free surface loop, check the material numbers of each element pair

114 0 0 0 #define MAXCONNECTIONS 50 struct freeElement{ int element; int nodes[3]; int surface; } *freeElements; int nFreeElements; int e,i,j,k,n,m,foundflag; int nFreeElements=0; int ns[4]; int ecount [3*MAXCONNECTIONSi [2]; 0 0 0 Maximum number of edges sharing a node. Structure for storing free elements (triangular faces). It contains: the element that the face belongs to, the three nodes that define the face,surface number. Some integers Common element counter for (n=l;<=nelements; n++) Every element for (i=1;i<=4;i++) Every face switch (i) { Get nodes for each face case 1; ns[0]=elements [n].node[0]; ns[1]=elements[n].node[1]; ns[2]=elements[n].node[2]; break; case 2; ns[0]=elements[n].node[0]; ns[1]=elements En].node[1]; ns[2]=elements[n].node[3]; break; case 3; ns[0]=elements [n].node[1]; ns[1 ]=elements[n].node[2]; ns[2]=elements n].node[3]; break; case 4; ns[0 =elements [n].node 0]; ns[] =elements [n].node[2]; ns[2]=elements[n].node[3]; break; } for (j=0;j<3*MAXCONNECTIONS;j++) { ecount [ j ] [ 0] =0; Zero out common element ecount[j] [1]=0; counter } for (j =0; j <=2; j ++) { Fill common element counter m=0; while (meshtree[nsj]] [m]!=0) { e=mesh_tree[ns [E]] [m]; k=0; while (ecount[k] [O]!=e&& ecount[k] [0]!=0) ++k);

115 if (ecount[k] [0]==O) ecount[k] [O]=e; ecount[k] [1]=1; I { else ++ecount[k] [I]; ++m; } k=O; foundflag=O; while (foundflag==O) { while (ecount[k] [1]!3 && ecount[k] [1]!=0) ++k; if (ecount[k] [0]==n) ++k; else foundflag=l; } Add new element to common element counter Increment common element counter Determine if face is free or shared Don't count self shared face if (ecount[k] [ 1 <3) { Test if free face ++nFreeElements; Free face freeElements [nFreeElements.Update freeElements (face) data element=n; structures. freeElements[nFreeElements]. nodes[0]=ns[0]; freeElements[nFreeElements]. nodes[1] =ns[1]; freeElements[nFreeElements]. nodes[2]=ns[2]; The following operations can be put here if needed. 0 0 Figure A.3 Operation to efficiently find all free surfaces in the mesh. Figure A.3: Operation to efficiently find all free surfaces in the mesh. sharing a face. If they are different, store the shared face as a free surface and calculate the normal in the same manner as above. Consistent normal directions are found by ensuring that the normal points towards the lower material index. if (elements[n].material < elements[ecount[k] [].material) { 0 0 0.3

116 *Find outward normals to free surface To ensure that all free face normals have the correct sign, the following inexpensive test will guarantee that all normals are pointing outward. Given an element with a free surface defined by the four nodal points (xi, yi, zi), i=1,2,3,4 (as shown in Figure A.4) the normal (to within a sign) is given by (r2 - r3) X (T4 - r3) (?r2-3) X (r4-r3)i (A.1) where ri = xxi + yyi + ZZi (A.2) node 4 ^ -Free face - node 2 node 1 ---- node 3 Figure A.4: The correct direction for the normal on a free face of an element. To find the correct sign, compare the distance from the node not on the free surface (node 1 for the example shown in Figure A.4) to a point that is a unit length along the normal in either direction. The direction giving the larger distance is the correct one.

117 *Divide free surface (find metal surfaces and outer surface of mesh) This algorithm divides the entire free surface found by the find free surface operation into disconnected pieces. To accomplish this, the mesh tree is filled with all of the free elements found by the find free surface operation and is recursively scanned to find groups of connected elements. Each surface element is tagged with a surface number and the volume of each surface is compared. The surface with the largest volume is then tagged as the outer surface and all other surfaces are tagged as metal surfaces. There can be multiple, unconnected metal bodies in the mesh. This operation is shown in Figure A.5 and Figure A.6. Now that the free surface has been split into pieces, the surface with the largest volume must be the outer surface and all others must be metal surfaces. To calculate the volumes of the surfaces, add up the contributions from each free face as N1ur] V = ( x v^) o (A.3) e= 1 where v - = x (x3- x2) + (Y.- Y2) + 2.(3 z2) 2 = (x4 - x2) + y(y4 - Y2) + (4 - 2) (A.4) an = x(x2) + y(Yc) + Z(z2) and (x., Y, Z) are the nodal coordinates of the element having a free face (see Figure A.4) and N/sur is the number of free faces in the surface.

118 for(n=l;n<=NFreeElements;n++) for(m=0;m<3;m++) { nl=FreeElements[n]. 4; Fill mesh tree with free surface elements only. i=0; while( mesh_tree[nl] [i]!=0) i=i+l; mesh_tree[nl] [i]=n; mesh_tree[nll[i+l]=0; } surface=1; Find connected free elements by recurwhile (surface!=0) { sively scanning the mesh tree. i=l; j=0; foundflag=0; while(foundflag==) { Drop out when mesh tree is empty while(mesh_tree[i][j]== -1) j=j+l; if (meshtree[i] [j]==0) if (i <= nvertices) { i=i+l; j=0; } else { foundflag=l; surface=0; } else { foundflag=l; e=mesh_tree[i] [j]; } if (surface!=O) { Recursively update mesh tree.. error=update(e,mesh_tree, FreeElements,surface); Keep count of distinct surfaces. nsurfaces=surface; surface=surface+l; Figure A.5: Operation to divide free surface into unconnected pieces. *Find edges lying on a conductor Now that all of the free surfaces have been found, each free face belonging to a metal surface (all surfaces except the surface with the largest volume) can be split into three edges. This information is necessary for edge-based finite element analysis. Z -- Z:V~I CV I~V~VVUrJ ~I VtV VHVLI~~V IIIVI CICI ~ ~

119 int update(e,meshtree,FreeElements, surface) struct FreeElement{ int element; int FreeNodes[l]; int FreeVertices[3]; int surface;; int e,surface; struct FreeElement *FreeElements; int **meshtree; { int n[3],j,m,f,i,error; n[O]=FreeElements[e].FreeVertices[O]; n[l]=FreeElements[e].FreeVertices[l]; n[2]=FreeElements[e].FreeVertices[2]; for (i=0;i<3;i++) { j=o; m=n[i]; while(meshtree[m][j]!=) { f=mesh_tree[m][j]; if (f!= -1) { FreeElements[f].surface= surface; mesh_tree[m][j] = -1; error=update (f,mesh_tree, FreeElements,surface); } j=j+l; } } return(!); Recursive function to remove connected free elements from mesh tree. Recursive call Figure A.6: Recursive function to remove connected free elements from the mesh tree. *Find enclosing integration contour This operation finds a surface (three-dimensions) or contour (two-dimensions) that completely encloses the target and passes through element centers. This routine relies on the previous ones since it requires the divided free surfaces and material interfaces. The mesh tree must be loaded with all elements (not just free surface elements).

120 The tree is then scanned for elements having the following properties: 1. The element's material is air (or other surrounding media) and has one, two or three nodes lying on a metal surface or a dielectric interface. 2. The element does not have more than one face on either a metal surface or a dielectric interface. Requirement 1 ensures that all of the scatterer (metal and dielectric) is enclosed. Requirement 2 ensures that the surface/contour will be closed and relatively smooth. Any element meeting these requirements is divided into two pieces by extracting all vertices which have one and only one node on a metal surface or dielectric interface. A new node is created at the midpoint of each vertex and a new face is created from these points. Taken collectively, these new faces make up a closed surface/contour which completely encloses the scatterer and passes through element centers where field derivatives can be calculated accurately. In Figure A.7, part of the integration contour for a two-dimensional coated conductor with a crack is shown.

121 y Integration contour Material coating Figure A.7: Integration contour around a coated conductor with a crack.

122 APPENDIX B Quasi-Minimum Residual (QMR) Solver The solutions of large, sparse linear systems which result in the practice of the finite element method for electromagnetics is obtained most appropriately by an iterative, minimum residual method. Large sparse systems require prohibitive memory demands when using a direct solver, and therefore the use of an iterative solver is a matter of necessity and not convenience. Iterative solvers based on the Conjugate Gradient (CG) method are not well suited to the linear systems resulting from electromagnetic analysis since these systems are not Hermitian. A pre-processing step is required for CG based methods in order to make the system Hermitian. This pre-processing step simply involves multiplying the linear system Ax = b (B.1) by its adjoint A (complex conjugate transpose) The CG method is then applied to the resulting system (A A)x = A b = A'x = b' (B.2) Unfortunately, this operation squares the condition number of the system, and in practice,

123 these poorly conditioned systems are prohibitively expensive to solve, and the solutions are numerically unreliable. An alternative to CG based methods is the Bi-Conjugate Gradient (BCG) method. While BCG algorithms can be used to solve the non-Hermitian systems that naturally occur in electromagnetics without squaring the condition number, BCG methods are prone to breakdowns. These breakdowns are manifested in a division by zero which occurs during the algorithm when the search vector becomes orthogonal to the excitation vector. In practice the breakdowns, or near breakdowns, result in wild oscillations of the norm of the residual are more than a nuisance. The following Quasi-Minimum Residual (QMR) algorithm (algorithm 5.1 in [43]) is well suited for the complex, non-Hermitian, unsymmetrical systems which result in electromagnetic analysis such as the application of FEM to problems presented in this thesis. Note that this algorithm has approximately the same convergence characteristics as the BCG method but is free from the numerical breakdowns that plague the BCCG method. The following implementation assumes that the sparse system is stored in SLAP (Sparse Linear Algebra Package / Lawrence Livermore National Laboratory) column format. The sparse matrix/vector product routine (DSMV) is not shown but is a standard routine in the SLAP software package. The linear system should be diagonally scaled (diagonal preconditioner) before calling this routine. The algorithm template is first given followed by an efficient implementation where memory optimizations have been introduced.

124 QMR Algorithm template: Solve Ax = b (1) Start: (a) Choose initial guess xO e C; (b) Set w1 =yl = ro = b - Axo, vO = Ayl, do = 0; T0 = Irol, 0o = 0, nio = 0; - -H (c) Choose random ro such that po = ro ro r 0 (orthogonality check). (2)Forn = 1,2,... do: -H Pn-1 (a) Set n - I = ro,, 1, n - 1 - 1; Y2n = Y2n- 1 n- I Vn- 1I; (b) Form = 2n- 1,2n do: Set wm + I =wm -a (,!AYm; m= lim+Ill/'Cm-,m,n = 1/( +0); Tm = Tm- 1mCm' lm = Cmgn-1; dm = Ym + (m-i m-I/, i)dm-; Xm = Xm- I+ rlmdm; If x, has converged: stop; -H (c) Set Pm = rOw2, + i, H,~ = P,,/P,,-!; y2, + i = I v,+ + nZyn; V, = Ay2,1 + + 3,(AY2, + /3v,,7 - )

125 C ' C * BLOCK: sprslpqmr Dan Ross C * C * DESCRIPTION: This routine solves the linear system Ax=b C * for x given A and b where A is complex/non-hermitian/ C * unsymmetric sparse stored in SLAP format. C* C * Based on: Roland W. Freund,"A Transpose-Free Quasi-Minimal C * Residual Algorithm for Non-Hermitian Linear Systems",SIAM C * vol. 13 1992 C * C* C * GLOSSARY: C A,IA,JA: SLAP format sparse system. C * b: Excitation vector C * x: Solution vector C * n: Size of system (nXn) C * tol: QMR Terminates when residual error<tol C * maxits: Maxinum # of iterations (the number of residual C * (error values will be twice this.) C * nnz: Number of nnz entries in sparse system. C reset: (1) Set initial guess to zero (0) use x C w,y 1,y2,r,v,d,temp,temp2: Temporary vectors C* Subroutine sprslpqmr(A,IA,JA,b,x,n,tol,maxits,nnz,reset, & w,y 1,y2,r,v,d,temp,temp2) IMPLICIT NONE Include 'bf3d.inc' Integer n,nnz,maxits,reset,IA(nnz),JA(nnz) Complex* 16 A(nnz),b(n),x(n),w(n),y 1 (n),y2(n),r(n),v(n), & d(n),temp(n),temp2(n) Complex* 16 tau,theta l,theta2,c,bnorm,rnorm,norm,ubound Complex* 16 eta 1,eta2,sigma,rho 1,rho2,alpha,beta,diprod Real*8 tol,error Integer success.seed,i,j,m Real randu Write (6,*) 'BEGIN QMR solver: nnz,n = ',nnz,n

126 C Initialize QMR success = 0 bnorm = norm(b,n) If (reset.EQ. 1) Then C No initial guess. Do 10 i=l,n x(i)=(0.0,0.0) 10 Continue End if call DSMV(n,x,temp,nnz,IA,JA,a,0) call vecsub(temp,b,temp,n) Do 20 i=l,n y l(i)=temp(i) w(i)=temp(i) d(i)=0 20 Continue call DSMV(n,y 1,v,nnz,IA,JA,a,0) tau = norm(temp,n) C Check if initial guess is exact solution. If (tau.EQ.0) Then error=0 Write (6,*) ':**WARNING:Exact Guess or Trivial solution' Return End If thetal=0 eta 1=0 seed=7 Write (6, ) 'QMR random seed = ',seed r( 1 )=Randu(seed) Do 30 i=l,n r(i) = 2.*Randu(seed)-l. 30 Continue rho 1 =diprod(temp,r,n) C Iterate untill Convergence Criteria has been met Do 40 i= 1,maxits C (a) sigma=diprod(v,r,n)

127 alpha=rho 1/sigma Call vecsac(y2,y 1,-1.*alpha,v,n) C (b) Do 50 m=2*i-1,2*i Call DSMV(n,y 1,temp2,nnz,IA,JA,a,0) Call vecsac(w,w,- 1.*alpha,temp2,n) theta2= norm(w,n)/tau c = 1./DSqrt( 1. + Dreal(theta2)**2) tau=tau*theta2*c eta2=c**2 * alpha Call vecsac(d,yl,(theta 1**2*etal/alpha),d,n) Call vecsac(x,x,eta2,d,n) If (m.EQ.2*i) Then ubound=DSqrt(DFloat(m)+ 1.)*tau error = Dreal(ubound)/Dreal(bnorm) If ( error.LT. 10000*tol ) Then call DSMV(n,x,temp,nnz,IA,JA,a,0) call vecsub(temp,b,temp,n) rnorm = norm(temp,n) error = Dreal(rnorm)/Dreal(bnorm) Write (6,*) 'Iteration: ',i,' Error = ',error If (error.LT.tol) Then success= 1 Goto 500 End If Else Write (6,*) 'Iteration: 'i, & 'Error bound = ',error End If If (QMRconvergence.EQ..TRUE.) Then Write (55,*) i,error End If End If theta =theta2 eta 1=eta2 Do 60 j=,n yi(j)=y2(j) 60 Continue 50 Continue C (c) rho2=diprod(w,r,n) beta=rho2/rho 1 Call vecsac(y2,w,beta,y 1,n) Call DSMV(n,y2,temp,nnz,IA,JA,a,0) Call vecsac(temp,temp,beta,temp2,n)

128 Call vecsac(v,temp,beta*betav,n) rho 1 =rho2 Do 65 j=,n y (j)=y2(j) 65 Continue 40 Continue 500 If (success.EQ.0.AND.tol.NE.0) then write (6,*) k*** ****** WARNING * > * ********' write (6,*) '***** Tolerance level not acheived ****' write (6,*) '\; *; * *;*;*:: 8 R * k k;:.********' End if Return End C * C BLOCK: vecsac C* C * DESCRIPTION: This routine scale a vector by a complex constant C * and adds it to another vector. C * result=vec l+scale*vec2 C* C ' INPUTS: (complex) vecl C * (complex) vec2 C (complex) scale C * (integer) n C OUTPUTS: (complex) result C* Subroutine vecsac(result,vec 1,scale,vec2,n) Complex* 16 result,vec 1,vec2 Complex* 16 scale Integer n Dimension result( 1 ),vec 1(1 ),vec2( 1) Do 10 i=l,n result(i)=vec 1 (i)+scaie*vec2(i) 10 Continue

129 Return End C* C ' BLOCK: norm C* C * DESCRIPTION: This function returns the C * Euclidean norm of a complex vector. C* C * INPUTS: (complex) vec C * (integer) n C * OUTPUTS: (real) norm C* Complex* 16 Function norm(vec,n) Complex* 16 vec Integer n Dimension vec(n) Complex* 16 xr,xi,temp,sum 1 sum1=0 DO 100 I=1,n xr=dreal(vec(I)) xi=dimag(vec(I)) temp = xr*xr + xi*xi suml = sumi + temp 100 CONTINUE norm=DSqrt( DReal(suml)) Return End ction-randu --- —-ix) --- -------------------------------- c --- —---------------------------------------- real function randu(ix) c

130 c Declare seed an integer, and init to 7. c then just call randi(seed) every time that you c want a real number between 0 and 1. c It is a uniform distribution, every number c between 0 and I being equally likely. c purpose: c This function calculates a uniform random number, c between 0 and 1. It uses a multiplicative congruential generator. c ie: x(i+l)=(x(i)*(7**5)) mod(2**31 -1). c See Kelton and Law for the algorithm used (p.227). c c input parameters: c c ix i*4 the latest random integer generated. c or the initial value input by user. c c output parameters: C *:' *:. i*:p. ~*:' *-. e -.* * *: ** * * * *. c The uniform random number is contained in RANDU c c commons used: c * * ** * * * * * * NONE c c files used: c* s*********** NONE c c subroutines called: c NONE c c c C --- —--------------------------------- -------------- c c ALPHABETICAL LISTING OF ALL VARIABLES USED c --- ------------ ----------- ---------- c c seed i'4 initial value for uniform random number c generator. in /randnum/ c c iy iY4 temp place to save product c

131 C --- —------------------------------ integrer a,p,ix,b15,b 1 6,xhi,xalo,leftlo,fhi,k data a116807/,b15/32768/,bl6/65536/,p/2 147483647/ c write(6,1I0) 1 0 format( 'entering, randu') xhi=ix/bl16 xalo=(ix-xhi*bl16)*a leftlo=xalo/b 16 fhi=xhP~a+leftlo k=fhi/b 15 ix=(((xalo-leftlo*b 1 6)-p)+(fhi-k*b 1 5)Th 1 6)+k if(ix.lt. 0) ix=ix+p randu=float(ix) *4.65 66 12875e- 1 0 c write(6,30)randu c 30 format('randu=',f 17. 10) c write(6,20) 20 format('leaving, randu') return end C CC BLOCK: diprod C DESCRIPTION: This routine dots two complex vectors. C result=vecI. vec2 C INPUTS: (com-plex) vecl~vec2 C (integer) n C OUTPUTS: (complex) diprod

132 Complex* 16 Function diprod(vec 1,vec2,n) Complex* 16 vec 1 vec2 Integer n Dimension vec 1 (1),vec2(1) diprod=0.0 Do 10 i=l,n diprod=diprod+vec 1 (i) * DCMPLX( DReal(vec2(i)),- 1. *Dlmag(vec2(i))) 10 Continue Return End C: C* C * BLOCK: vecsub C* C * DESCRIPTION: This routine subtracts two vectors. C * result=vec 1-vec2 C * C INPUTS: (complex) vecl,vec2 C * (integer) n C * OUTPUTS: (complex) result C Subroutine vecsub(result,vec 1,vec2,n) Complex* 16 result,vec 1,vec2 Integer n Dimension result( 1),vec 1(1 ),vec2(1) Do 10 i=l,n result(i)=vec 1 (i)-vec2(i) 10 Continue Return End

BIBLIOGRAPHY 133

134 BIBLIOGRAPHY [1] S. Qian and D. Chen, "Signal Representation Using Adaptive Normalized Gaussian Functions." Signal Processing, vol. 36, pp. 1-11, Mar. 1994. [2] L.C. Trintinalia and H. Ling, "Extraction of Waveguide Scattering Features Using Joint Time-Frequency ISAR," IEEE Microwave and Guided Wave Letters, vol. 6, no. 1,pp. 10-12, Jan. 1996. [3] R.E. Gardner, Report 5656, Naval Research Laboratory, Washington DC, Aug. 1961. [4] M.R. Bell and R.A. Grubbs, "JEM Modeling and Measurement for Radar Target Identification," IEEE Transactions on Aerospace and Electronic Systems, vol. 29, no. 1, pp. 73-87, Jan. 1993. [5] J.R. Coleman, "Investigation of Radar Scattering by Simulated Aircraft Duct/Engine Combinations," Northrop Corporation, Technical Report AFAL-TR-73-361, Nov. 1973. [6] T.W. Johnson, "Electromagnetic Scattering by Open Circular Waveguides," Ph.D. Dissertation, Ohio State University, Columbus, Ohio, 1980. [7] P.P. Silvester and R.L. Ferrari, Finite Elements for Electrical Engineers. Cambride: Cambridge University Press, 1986. [8] J.L. Volakis, A. Chatterjee and L.C. Kempel, "A Review of the Finite Element Method for Three-Dimensional Electromagnetic Scattering," Journal of the Optical Society of America - A, pp. 1422-1433, Apr. 1994. [9] D.C. Ross, J.L. Volakis and T. Ozdemir, "A New Finite Element Formulation for Modeling Circular Inlets with Irregular Terminations," Technical Report 030395-1 -T, The University of Michigan, EECS Radiation Laboratory, Ann Arbor MI, May 1993. [10] H. Ling, R. Chou and S.W. Lee, "Shooting and Bouncing Rays: Calculating the RCS of an Arbitrary Shaped Cavity," IEEE Transactions on Antennas and Propagation, vol. AP-37, no. 2, pp. 194-205. Feb. 1989. [11] P.H. Pathak and R.J. Burkholder, "High Frequency EM Scattering by Open-Ended Waveguide Cavities," Radio Science, vol. 26, no. 1, pp. 211-218, Jan.-Feb. 1991. [12] 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 Transactions on Antennas and Propagation, vol. AP-36, no. 1, pp. 84-96, Jan. 1988. [13] H.T. Anastassiu and J.L. Volakis, "The Mode Matching Technique for Electromagnetic Scattering by Inlets with Complex Terminations," Technical Report 030395-3 -T, The University of Michigan, EECS Radiation Laboratory, Ann Arbor MI, October 1993.

135 [14] D.C. Craig and J.L. Volakis, "Body of Revolution Code (DBR) Programming Guide," Technical Report 389604-8-T, The University of Michigan, EECS Radiation Laboratory, Ann Arbor MI, May 1989. [15] C. Huang, "Ray Analysis of EM Backscatter from a Cavity Configuration," Ph.D. Dissertation, The Ohio State University, ElectroScience Laboratory, Columbus OH, 1982. [16] P.Y. Ufimtsev, "Method of Edge Waves in the Physical Theory of Diffraction," (in Russian) Izd-Vo Sov. Radio, pp. 1-243. 1962. (English Translation: U.S. Foreign Tech. Office Div., Doc. ID FTD-HC-23-259-71, Wright-Patterson AFB, Dayton OH, 1971). [17] P.H. Pathak, C.W. Chuang, M.C. Liang, "Inlet Modeling Studies," Technical Report 717674-1, The Ohio State University, ElectroScience Laboratory, Columbus OH, Oct. 1980. [18] P.H. Pathak and R.J. Burkholder, "A Reciprocity Formulation for the EM Scattering by an Obstacle Within a Large Open Cavity," IEEE Transactions on Microwave Theory and Techniques, vol. MTT-41, no. 4, pp. 702-707, Apr. 1993. [19] D.R. Lynch and K.D. Paulsen, "Origin of Vector Parasites in Numerical Maxwell Solutions," IEEE Transactions on Microwave Theory and Techniques, vol. MTT-39, no. 3, pp. 383-394, Mar. 1991. [20] M.L. Barton and Z.J. Cendes, "New Vector Finite Elements for Three-Dimensional Magnetic Field Computation," Journal of Applied Physics., vol. 61, no. 8, pp. 3919 -3921, Apr. 1987. [21] W. Schroder and I. Wolff, "The Origin of Spurious Modes in Numerical Solutions of Electromagnetic Field Eigenvalue Problems," IEEE Transactions on Microwave Theory and Techniques, vol. MTT-42, no. 4, pp. 644-653, Mar. 1994. [22] K.D. Paulsen, W.E. Boyse and D.R. Lynch, "Continuous Potential Maxwell Solutions on Nodal-Based Finite Elements," IEEE Transactions on Antennas and Propagation, vol. AP-40, no. 10, pp. 1192-1200, Oct. 1992. [23] B.M. Azizur-Rahman and J.B. Davies, "Penalty Function Improvement of Waveguide Solutions by Finite Elements," IEEE Transactions on Microwave Theory and Techniques, vol. MTT-32, no. 8, pp. 922-928, Aug. 1984. [24] A.B. MacFarland, "Parameter Centering and Tolerancing," Master's Thesis, Electrical Engineering Graduate School, University of Idaho, Aug. 1986. [25] F. Brezzi, M. Bristeau, L.P. Franca, M. Mallet, G. Roge, "A Relationship Between Stabilized Finite Element Methods and the Galerkin Method with Bubble Functions," Computer Methods in Applied Mechanics and Engineering, vol. 96, pp. 11 7 -129, 1992. [26] R. Pierre, "Simple C~ Approximations for the Computation of Incompressible Flows," Computer Methods in Applied Mechanics and Engineering, vol. 68, pp. 205 -227, 1988.

136 [27] J. Van Bladel, Singularities of the EM Field. Oxford: Oxford University Press, 1990. [28] D.C. Ross, J.L. Volakis and H.T. Anastassiu, "Hybrid Finite Element-Modal Analysis of Jet Engine Inlet Scattering," IEEE Transactions on Antennas and Propagation, vol. AP-43, no. 3, pp. 277-285, Mar. 1995. [29] R. Lee and T. Chia, "Analysis of Electromagnetic Scattering from a Cavity with a Complex Termination by Means of a Hybrid Ray-FDTD Method," IEEE Transactions on Antennas and Propagation, vol. AP-41, no. 11, pp. 1560-1569, Nov. 1993. [30] E.M. Nelson, "High Accuracy Electromagnetic Field Solvers for Cylindrical Waveguides and Axisymmetric Structures Using the Finite Element Method," Ph.D. Dissertation, Stanford University, Stanford CA, 1993; also SLAC-431. [31] A. Frenkel, J.R. Brauer and M.A. Gockel, "Complex Periodic Boundary Conditions for AC Finite Element Models," Proceedings of the 10th Annual ACES Conference, pp. 195-202, Mar. 1994. [32] J.R. Mautz and R.F. Harrington, "Radiation and Scattering from Bodies of Revolution," Applied Science Research, vol. 20, pp. 405-435, Jun. 1969. [33] S.D. Gedney, J. Lee and R. Mittra, "A Combined FEM/MoM Approach to Analyze the Plane Wave Diffraction by Arbitrary Gratings," IEEE Transactions on Microwave Theory and Techniques, vol. MTT-40, no. 2, pp. 363-370, Feb. 1992. [34] A.W. Glisson and C.M. Butler, "Analysis of a Wire Antenna in the Presence of a Body of Revolution," IEEE Transactions on Antennas and Propagation, vol. AP-28, no. 9, pp. 604-609, Sep. 1980. [35] A.E. Siegman, "Quasi Fast Hankel Transform," Optics Letters. vol. 1. no. 1, pp. 13 -15, July 1977. [36] D.C. Ross, J.L. Volakis and H.T. Anastassiu, "Overlapping Modal and Geometric Symmetries for Computing Jet Engine Inlet Scattering," IEEE Transactions on Antennas and Propagation, vol. AP-43, no. 10, pp. 1159-1163, Oct. 1995. [37] B. Jian, J. Wu and L.A. Povinelli, "The Origin of Spurious Solutions in Computational Electromagnetics," NASA Technical Memorandum 106921, COMP-95-8, May 1995. [38] J.M. Jin and J.L. Volakis, "Electromagnetic Scattering by and Transmission through [38] J.M. Jin and J.L Volakis Z 1) Z-t-) a Three-Dimensional Slot in a Thick Conducting Plane" IEEE Transactions on Antennas and Propagation, vol. AP-39, no. 4, pp. 543-550, Apr. 1991. [39] A. Chatterjee, J.M. Jin and J.L. Volakis, "Edge-Based Finite Elements and Vector ABCs Applied to 3D Scattering," IEEE Transactions on Antennas and Propagation, vol. AP-41, no. 2, pp. 221-226, Feb. 1993. [40] R. Dyczij-Edlinger and 0. Biro, "A Joint Vector and Scalar Potential Formulation for Driven High Frequency Problems Using Hybrid Edge and Nodal Finite Elements," IEEE Transactions on Microwave Theory and Techniques, ol. MTT-44, no. 1, pp. 15-23, Jan 1996.

137 [41] J.Gon, J.L.Volakis, A.C. Woo and H.T.G Wan, "A Hybrid Finite Element-Boundary Integral Method for the Analysis of Cavity-Backed Antennas of Arbitrary Shape," IEEE Transactions on Antennas and Propagation, vol. AP-42, no. 9, pp. 1233-1242, Sep. 1994. [42] T.B.A. Senior and J.L. Volakis, Approximate Boundary Conditions in Electromagnetics. London: Institute of Electrical Engineers, 1995. [43] R.W. Freund, "A Transpose-Free Quasi-Minimal Residual Algorithm for Non-Hermitian Linear Systems," RIACS Technical Report 91.18, NASA Ames Research Center, Sep. 1991. [44] D.C. Ross, J.L. Volakis, "Users Manual for FEMENG," Technical Report 033279-5 -T, The University of Michigan, EECS Radiation Laboratory, Ann Arbor MI, Feb. 1996. [45] D.C. Ross, "Some Finite Element Preprocessing Algorithms for Electromagnetic Scattering," IEEE Antennas and Propagation Magazine, vol. AP-35, no. 3, pp. 68 -72, Jun. 1993. [46] Z.S. Sacks, D.M. Kingsland, R. Lee, and J. Lee. "A Perfectly Matched Anisotropic Absorber for Use as an Absorbing Boundary Condition," IEEE Transactions on Antennas and Propagation, vol. AP-43, no. 12, pp. 1460-1463, Dec. 1995. [47] 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, vol. 9, no. 11/12, pp. 1363-1391, Nov.-Dec. 1995. [48] R.D. Graglia, A. Freni, and G. Pelosi, "A Finite Element Approach to the Electromagnetic Interaction with Rotating Penetrable Cylinders of Arbitrary Cross Section," IEEE Transactions on Antennas and Propagation, vol. AP-41, no. 5, pp. 635 -650, May 1993. [49] J. Van Bladel, "Electromagnetic Fields in the Presence of Rotating Bodies," Pro[49] J Van Bladell, In 1-n ceedings of the IEEE, vol. 64, no. 3, pp. 301-318, Mar. 1976. [50] G. Crabtree, W. Huegle, D. Salisbury, "RCS Compact Range Test Results for a Set of Simplified Engine Face Models," Report TM94226, GE Aircraft Engines, One Neumann Way, Cincinnati OH, Nov. 1994.