030395-1-T A new finite element formulation for modeling circular inlets with irregular terminations D. Ross, J. L. Volakis and T. Ozdemir Veda Inc. Suite 200 5200 Springfield Pike Dayton OH 45431 Wright Laboratory/AARA Wright-Patterson AFB OH 45433-7001 May 1993 30395-1-T = RL-2418

A new finite element formulation for modeling circular inlets with irregular terminations D. Ross, J. L. Volakis and T. Ozdemir Department Radiation Laboratory of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109-2122 Abstract This report describes a finite element formulation for computing the fields near the termination of a circular inlet. New node-based divergenceless tetrahedral elements are developed with 12 degrees of freedom and a new artificial absorber is described for terminating the finite element mesh only 0.3A away from the inlet termination. The developed formulation and code are validated by employing it to compute cavity resonances and the reflected mode coefficients in a circular cavity with a circular stub termination. 1

1 Introduction This report presents the current state of the jet engine cavity project that was undertaken in January 1993. It was proposed that a hybrid finite element/ modal technique be used with the finite element method (FEM) applied to the engine face and modal analysis used in the cavity (see Figure 1). It is intended, for practical application, that the incident field on the engine face (found by some other technique: ray tracing, GTD etc.) be decomposed into its constitutive modes. Each of these modes will act as an independent incident field, giving rise to a set of outgoing modes. The engine face can therefore be considered to be an N-port network where N is the number of travelling modes in the cavity. The FEM software will be used to generate the scattering parameters which will fully describe this network. The FEM analysis (which is computationally intensive) need only be done once for a given cavity, and the results stored as a scattering matrix. An important step in the validation of the method is the analysis of a geometry for which there is a rigorous solution. For this, a simple, relatively small inlet consisting of a concentric, circular, cylindrical termination (instead of an actual engine face) was chosen (see Figure 2). This configuration has been analyzed with a mode-matching technique for comparison with the FEM/modal results. In this report, the full, three-dimensional FEM formulation is given with an explanation of how the numerical shape functions are derived. Numerical shape functions allow the software to be written in a general way so that different approximation functions (basis) can be employed in each tetrahedral element. Among those tested was a new node-based, bilinear, divergencefree element. It will be shown how the FEM formulation and software have been validated by computing the resonant frequencies of metallic cavities (eigenvalue problem). The FEM software has been implemented for jet engine cavity analysis on a mode-by-mode basis. That is, for a given face and a given incident mode, the scattered field from the face is computed. The coefficients for the outgoing modes are then extracted. Results will be shown for the simple test inlet mentioned above. 2

2 Three-dimensional FEM formulation Consider the geometry in Figure 2. The volume enclosed by the metal walls of this geometry makes up the entire domain of the FEM solution. This volume, to be denoted by f, is enclosed by the surface F. The electric field in Q must satisfy the weak form of the wave equation {-iV x E *. V x W- k 2rEs * } dv Y r I -W n. x V x E'ds r r = W. {Vx -V x EinC - koerEi} dv /Itr + [I - ] W i x V x Eincds (1) ad /a2 /lj where W: weighting function E: scattered field Einc: incident field ko: free space wave number Fd: dielectric interface surface nd: normal to dielectric interface Special consideration must be given to the second term on the LHS of (1). Since we have chosen to use a metal-backed material absorber for terminating the computational volume (i.e., the FEM mesh), the scattered field on I must satisfy the boundary condition n x Es = -n x Ei on F Therefore, the additional unknowns V x E on F will uncouple from the system, and the second term on the LHS drops out. As an alternative to the material absorber, we could have used an absorbing boundary condition (ABC). An ABC is some differential operator D(.) that relates n x V x Es to Et (tangential component of E), i.e. n x V x Es = D(E) 3

One choice of D is [1] D(E') = aEt +,/V x [ni(V x ES)/,]Vt(V ES) where a = jko, 3 = 2- and this ABC has been employed in [1] for 3D scattering. The choice to use a material absorber instead of an ABC is based on the fact that a material absorber has been shown to behave more like free space than an ABC [2]. Moreover, the artificial absorber avoids numerical difficulties at the junction of the ABC with the metal surfaces. As part of this study a new artificial absorber was designed and illustrated in Figure 3. The effectiveness of this absorber was evaluated by placing it behind a threedimensional slot and observing how well the fields passing through the slot are absorbed. For a perfect absorber the fields computed by the configuration in Figure 3(b) must be identical to those of Figure 3(a). In Figure 4 we show the scattering by each of the three configurations in Figure 3, and it is clear that the proposed absorber is quite effective. These results were generated by using an existing three-dimensional finite element-boundary integral code described in [3, 4]. Note that the absorber is only 0.15A thick and its upper surface is only 0.15A away from the structure. Another possibility is to use a. boundary integral to exactly model the inside of the infinite cylinder. In this case, the Green's function for the interior must be found. It was decided that this approach would be impractical due to the complexity of the Green's function for the interior of the cavity. Also, in practice, the inlet cavity will not be circularly cylindrical, and the Green's function cannot in general be found. 3 Nodal-based numerical elements Standard node-based FEM element expansion functions are of the form N ES = Es = xE N(zx, y, z) E a=l N + Na (, y, z)Ea a=1 N + z Z Na(x, y, z) Ez a=1 4

where Et is the ith component of the scattered field at node a, and Na are referred to as the element shape functions. By choosing Wb = (x + y + )) Nb(x, y, z), i.e. weight and basis functions are the same, we can proceed with all of the same arguments normally used in variational methods to prove convergence, accuracy etc. [5]. Assuming that we are using tetrahedral elements (to facilitate mesh generation), we proceed with the derivation of the shape functions by approximating the field in the ith tetrahedron as E = al + a2x + a3y + a4z (2) where i denotes x, y or z. Next we assign Es the value E? at each of the four nodes (a = 1,2, 3, 4). In matrix form, El 1 1 Y ' 1 l 1 al E 1 X Y2 Z2 a2 EX 1 3 Zy3 a3 E4 1 X4 Y4 Z4 a4 or, {E} = [F]{a,} Solving for the a, coefficients yields {an}= [F]-1 {E} and upon substituting into (2) and grouping terms we obtain 4 E= Na(x,y,z) E a=1 {Na(, y, z)} [F]) {f(x,y, z)} (3) The superscript T denotes the transpose of the matrix [F]-1, fi(x,y,z)= 1, f2(x,y,z) = x f3(x, y, z) = y, f4(x, y, Z) = z 5

It is noted that the above procedure for deriving the shape functions can be generalized for any set of {fj(x, y, z)}. Specifically, if we set N E' - = aj fj(x, y, z) j=1 by following the same procedure we find El fi(Xi,yi,zi) f2(Xi, yI1, Z) f3(xi, yiZi) E] f (2f ( y2, z2) f2(X2, Y2, z2) f3(X2, Y2, Z2) E3 -= X fi(x3,3, z3) f2(x3, Y3, Z3) f3(3, Y3, Z3) EN fl(XN,YN, N) f2(xN, YN,N) f3(xN, YN,N) a1 a2 a3 aN or {E)} = [F]NxN {an}. The shape functions are then found from (3). One of the problems that have been found in three-dimensional FEM analysis is the occurrence of spurious modes due to the existence of nondivergenceless solutions to (1). If the shape functions can themselves be made divergence free, this difficulty is eliminated. This is the motivation behind edge-based elements [6]. Edge-based elements, however, are awkward to work with since they do not conform to standard FEM pre/post-processing techniques. Furthermore, their use doubles the unknown count. It is possible to construct nodal-based divergence-free elements by enforcing the divergenceless condition in the construction of the numerical shape functions as outlined above. This was first attempted by coupling the three field components together such that V. E' = 0 using standard linear elements. While this technique works at the element level, after assembling all elements and enforcing continuity, it was found that the elements lock. That is, the enforcement of the divergenceless condition completely (or almost completely) consumes all degrees of freedom associated with the wave equation. This phenomenon has also been observed in the field of incompressible fluid flow [7]. A set of divergence-free nodal-based shape functions that do not lock was created by adding an additional node at the center of the element and using the expansion Ex = a + a2+a3y+ a4z + a5yz 6

Ey = a6 + a7x + asy + agz + aloxz Ez = all + al2x + al3y + al4z + al5xy The divergence-free condition for this expansion is a2+ as + a13 =0 (4) and it should be noted that the coefficients of the bilinear terms do not affect the divergence-free condition. To enforce (4) along with the nodal equations, one of the 15 degrees of freedom must be discarded. Arbitrarily, the z-component of the field at the center can be considered redundant with the divergence-free condition. Consequently, the equations that we solve for generating the shape functions are E \ /1 I1 yi zl (Ylzi) / al E2 1 x2 Y2 Z2 (Y2Z2) a2 E3 = 1 X3 Y3 Z3 (y3Z3) a3 (5) E4I 41 X4 Z4 (y4z4) a4 E5 \ 1 xs y5 Z5 (ys5) a5 Ey1\ 1 x, yJ zl (xlzl) \ a6 Ey2 1 X2 Y2 Z2 (X2Z2) a7 E = 1 X3 y3 Z3 (X3Z3) as (6) E 1 X4 Y4 Z4 (X4Z4) a9 Ey5 1 X5 y5 Z5 (x5z5) \ alo E \ 1 X1 yl z, (Xly) all Ez2 1 X2 Y2 Z2 (x2Y2) a2 7) 3 (7) E3 1 X3 Y3 Z3 (X3Y3) a13 \ Ez4 1 X4 Y4 Z4 (X4y4) a14 and a2 + as + a13 = (8) where I 4 5 = E xi i=1 1 4 Y = 4E Yi 4 i=1 7

i 4 z5 = z The shape functions are now vectors and can be written as N = xNj + yNY + zNZ for a = 1,2,3,4,5 except that Nz = 0. They are obtained as before by inverting the system of equations listed in (5)-(8). Specifically, 1N:) {f;} {NeA} = ([F-l]) {f{ } {N:}J {f} J in which [F] denotes the matrix of the system implied by (5)-(8), fl = 1, f2 =x, f =y, f =z, f = yz f = 1i,f = f3y = Y, f4 =, f5z =xz and f = 1, f2=x, f3 =, f4 =z, f = xy In practice, the degrees of freedom associated with the node at the center of an element can be statically condensed at the element level so that only 12 degrees of freedom remain. These elements have been implemented, and tested successfully for the eigenvalue problem (see next section). They do, however, suffer from an inherent ill-conditioning when two of the elements nodes lie on the x, y or z axis. If, for example, we consider the x-axis, then the x-component of the field at these two nodes will exactly specify the factor ONa ax everywhere on the element. Consequently, the enforcement of the divergencefree condition at the central node leads to a singular system. It is then necessary to use a different coordinate system for each element so as to avoid this phenomenon. 8

The judicious choice of a coordinate system for each element can become a very intensive task, and no good scheme for doing this has been developed yet for large-scale problems (hundreds of thousands of elements). An alternate technique for eliminating the non-divergenceless modes is to add the term I (V. W)(V E) dv to the LHS of (1). This has the effect of shifting the frequency of the nondivergence-free modes up beyond the sampling rate of the mesh. Consequently, the spurious modes do not affect the solution [8]. This technique has also been successfully implemented using linear elements. 4 Validation/resonance computations The finite element method can be used to calculate the resonant frequencies of closed metal cavities by writing (1) as ([K] - A[M]){} = 0 where [K] = { V xN, V x Nb d P r [M] = rNa N.Nbd Qt: element domain A = k2 = W2pLe0 The term fAe 1(V ' Na)(V Nb)dv can be added to [K] so that the eigenvalue (A) corresponding to a non-divergence-free eigenvector {4} will be shifted up beyond the sampling rate of the mesh. In this way, conventional, node-based elements can be used. The above was implemented for bi-linear divergence-free elements and conventional linear elements. It was found that with the added expense of finding an appropriate coordinate system for each bilinear element, conventional elements are the preferred choice for large practical problems. 9

The eigenvalues for a rectangular cavity with dimensions 1A x 0.5A x 0.75A were found (122 degrees of freedom) (0.2A mesh size) with no spurious modes present. hi-linear (divergenceless) linear Mode analytic calculated calculated TElo 5.236 5.26 5.6 TMllo 7.025 7.35 7.8 TEoi 7.531 7.43 8.7 TE201 7.45 8.7 TM111 8.179 8.10 9.6 TElll 8.12 9.6 TM210 8.886 8.96 10.6 TE102 8.947 10.49 10.6 Another example run with our preliminary code pertains to the inlet termination displayed in Figure 2. As seen, the termination is a circular stub of radius b = 0.503A. The outer radius of the inlet is a = 1.66A, and the artificial absorber's face was placed A/3 away from the stub. This circular inlet termination was excited with the TE~eve mode, and the finite element code was used to compute the total fields everywhere within the "cavity". Given the total field, the reflected mode coefficients were computed by evaluating the inner-product of the cross-sectional fields with the circular inlet eigenmodes. The table below compares the mode coefficient amplitudes as computed by the finite element code and the (exact) mode-matching code described in [9]. Amplitudes of the reflected mode coefficients for the inlet in Fig. 2 with b = 0.503A and a = 1.60A Mode Finite element Mode-matching TE11 0.632 0.625 TE12 0.512 0.49 TE13 0.22 0.20 TE14 0.105 0.05 TE15 0.04 0.01 As seen, the mode coefficients are in good agreement, and this provides a validation of the finite element code. 10

5 Future work The goal of this project is the development of a general-purpose finite element code for computing the fields reflected by jet engine inlet terminations. At this stage, we have developed and provided a basic validation of the code. Over the next few months our intent is to provide a more extensive validation of the code and to further enhance its user-oriented features. Specifically, we will exercise the code in computing the mode scattering coefficients associated with more complicated inlet termination which better represent a typical engine face. A reference code will be developed for this purpose to validate the finite element code. Perhaps even more important tasks are: (a) an assessment of the code's computational demands for modeling large size inlets; (b) a study on the (minimum) parameters required for a complete characterization of the field reflected by the jet engine termination; and (c) development of techniques for coupling the finite element fields to different ray techniques. References [1] A. Chatterjee, J.M. Jin and J.L. Volakis, "Edge-based finite elements and vector ABC's applied to 3-D scattering," IEEE Trans. Antennas Propagat., vol. 41, pp. 221-226, Feb. 1993. [2] J.L. Volakis and T. Ozdemir, "Field calculations for finite slots using artificial absorbers for terminating the finite element grid," in preparation. [3] J.M. Jin and J.L. Volakis, "Electromagnetic scattering by an aperture formed by a rectangular cavity recessed in a ground plane," J. of Electromagnetic Waves and Application, vol. 5, no. 7, pp. 715-734, 1991. [4] J.M. Jin and J.L. Volakis, "A hybrid finite element method scattering and radiation from microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas and Propagat., vol. AP-39, pp. 1598-1604, Nov. 1991. [5] G. Strang and G. Fix, An analysis of the finite element method, Englewood Cliffs, N.J.: Prentice-Hall, 1973. 11

[6] M.L. Barton and Z.L. Cerdes, "New vector finite elements for threedimensional magnetic field computation," J. Appl. Phys., vol. 61, no. 8, pp. 3919-3921, April 1987. [7] D.S. Malkus and T.J.R. Hughes, "Mixed finite element methods-reduced and selective integration techniques: a unification of concepts," Computer Methods in Applied Mechanics and Engineering, vol. 15, no. 1, pp. 63-81, 1978. [8] K.D. Paulsen and D.R. Lynch, "Elimination of vector parasites in FE Maxwell solutions," IEEE MTT Transactions, vol. 39, no. 3, pp. 395 -404, March 1991. [9] H.T. Anastassiu and J.L. Volakis, "A mode matching technique for electromagnetic scattering by circular inlets with canonical terminations," Univ. of Michigan Radiation Laboratory Report #030395-2-T. 12

Connectivity Boundaries SBR engine Finite element region Connectivity bounds Figure 1. Illustration of the hybrid finite element-high frequency method for jet engine inlet analysis.

Metal Artificial absorber Figure 2. Illustration of the artificial absorber termination scheme.

(a) FE-BI Solution (Slot Code Used) ___ w - - ---- PEC t D | PEC L (b) FE-Artificial Absorber (Cavity Code Used) ' W w Absorber PEC (c) Cavity (Cavity Code Used) Hw-H PEC D D PEC PEC Figure 3. Geometries referred to in Figure 4. NOTE: All the slots mentioned in this report have square apertures.

20 10 I r -— I vc L ---i C14 -E C) a:) t) 0 -10 -20 -30 -40 20 10 0 -10 -20 -30 -40 recslotW0.525_D0.075, -=0 FE-BI Solution ' —. - Eo FE-Artificial Absorber ------ Cavity - I, I.! I. I. I. (a) 0 10 20 30 40 50 60 70 80 90, 0 [deg] recslot_W1.0_D0.2, (=0 r-. (D 0> (b) I 0 10 20 30 40 50 60 70 80 90 0 [deg] Figure 5. Comparison of the RCS data generated by the finite element slot code using the boundary integral (exact) and the artificial absorber for truncating the mesh cavity data (slot with PEC bottom) have been included to illustrate the amount of energy absorbed by the artificial absorber. (a) Backscattering by a rectangular slot 0.525k x 0.525k and 0.075X thick; (b) backscattering by a rectangular slot 1X x 1I and 0.2k thick.