035675-1-T USER'S MANUAL FOR AIMJET Hristos T. Anastassiu, Mikhail Smelyanskiy, and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan 1301 Beal Avenue Ann Arbor MI 48109-2122 FAX: (313) 647-2106 hristosa@engin.umich.edu msmelyan@engin.umich.edu volakis@engin.umich.edu 35675-1-T = RL-2494

1 Introduction Aimjet is a FORTRAN77 code which uses the Adaptive Integral Method (AIM) [1] to model electromagnetic scattering from perfectly conducting jet engine models. It is based on a free space Moment Method code written by Pamela Haddad at the Radiation Laboratory in 1991, on a free space AII code developed by Hristos Anastassiu in 1995 and a jet engine MoM code written by Hristos Anastassiu in 1997. Section 2 of this manual gives a brief description of AIM applied to the engine problem. Section 3 describes the code operation, and section 4 presents a few examples of Radar Cross Section calculations. 2 AIM applied to the jet engine problem A major reason for an emphasis on the EM characterization of jet engines is their critical importance in determining the airplanes's reflectivity and for target identification. Even though the engine is a small component of the overall structure, its complexity and visibility in the front sector of the aircraft make its contribution a large component of the overall airframe scattering. Moreover, the rotating engine blades cause unavoidable doppler shifts which are important for target identification purposes. The latter can be an important discriminator for both civilian and military applications and more reliable than other methods. As applied to electromagnetic scattering problems, the Moment Method (MoM) has proven to be a very accurate integral equation technique. However, its brute force application to the jet engine problem is not feasible due to excessive electrical size of the structure. To reduce the CPU time and storage requirements down to manageable levels, the inherent blade periodicity of the jet engine was exploited to show that the computational domain can be reduced down to a single blade or engine "slice" [2, 3]. In free space we have demonstrated [2] that for large scatterers the Adaptive Integral Method (AIM) is much more efficient than MoM in terms of memory and CPU time, yet it retains MoM accuracy. Therefore, it is expected that application of AIM to the jet engine 2

0 Wrong axis alignment I - -, ~..... / Correct axis alignment \ / Figure 1: Cylindrical grid (front view). problem will demonstrate similar merits. Although the application of the slicing scheme greatly enhances the efficiency of MoM, the benefits for AIM are not that impressive without substantial modifications of the method for its application to this structure. The main reason is that the Toeplitz property of the Green's function matrix can only be retained by extending the rectangular grid over the whole scatterer, and not over a single slice. The details of this problem are explained in [2]. To avoid using a global rectangular grid and, hence, large dimensions for the FFT, the only possible solution is to construct a grid that still retains the Toeplitz property for the Green's function matrix [G]. Such a grid is not rectangular, but cylindrical, as depicted in Fig. 1. In cylindrical coordinates the Green's tensor is not Toeplitz along the p direction, therefore the FFT can only be used along the 4 and z directions, meaning that matrixvector product calculations are slower than in the case of the rectangular grid. This drawback of the cylindrical grid is balanced by the smaller size of the FFT dimensions, since the computational domain is reduced down to a single slice. 3

Tha mathematical details of the implementation are out of the scope of this manual. but all essential information can be found in [2]. 3 Code Operation 3.1 Code structure Aimjet models the scattering from metallic cylindrical jet engine inlets terminated by an arbitrary, PEC, cylindrically periodic structure. A simplified geometry is illustrated in Fig. 2, where the termination is a cylindrical hub. A plane wave is assumed to illuminate the open end and couple to waveguide modes in Region 1, propagating towards the termination (Region 2). Each of these modes is used as an excitation for AIM applied to one periodic sector of the termination (see Figs. 3, 4). Solving the linear system via AIM yields the equivalent currents on the termination, and in turn the scattered field within the inlet. To extract the coefficients of the scattered modes propagating towards the open end an integration [2, 4] is performed at a hand-off plane, which is perpendicular to the inlet axis and lies at a fraction of a wavelength in front of the termination (see Fig. 5). Each of the scattered modes propagates to the open end and radiates in free space. The overall modal contribution yields the Radar Cross Section (RCS) of the engine. To exploit the inherent periodicity of the engine, the waveguide modes must have exponential, and not trigonometric dependence on q. Immediate consequence of this is that the modal orders may be negative, positive or zero, and there is no distinction between "odd" and "even" behavior, as opposed to the Mode Matching code that has been developed in the past [6]. 4

, L.761V.U Is a Region 2 - L -1 1- 1 ------- z=0. Figure 2: Cylindrical inlet terminated by a cylindrical hub 5

C44 '1 CD CD Si 4

7 A'-tr~~:b '\* 4L# ":r Figue 4 quaterof ablad geoetr 7

Modal or Ray region Fields to Ray/Modes Numerical or Fids on Surface Rigorous Modeling (Hand-off Surface) EEngine F omposi he computational domain. Figure 5. Decomp has thare tesform wall up to the hand-off lan ra based format m1 g o The current version of Mornet requires an e fla-tPlate (see Fg ecebti the filee formt wie goometr ris not typical for th plane and centered at the ormng format of the input file. The geomet simplicity is exp has the form 8

98.1.1.0.1 -.1.0 -.1.1.0 -.1 -.1.0.1.0.0.0.1.0 -.1.0.0.0 -.1.0.0.0.0 16 8 78 56 59 69 39 79 89 29 48 47 16 15 37 36 28 25 1 9 10 2 1112 234 5 6 13 4 5 14 167 8 15 7 16 8 3 (N nodes = 9, Ntriangles = 8) (the next ANlnodes entries are the cartesian coordinates of the nodes) (total number of edges Nedges ) (number of exterior edges Next) (next Nedges - Next entries give the pair of nodes defining the edges) (the remaining entries give the triangle connectivity. Each line entry contains 3 numbers corresponding to the edge numbers forming the triangle in a clockwise order) To obtain this format: 9

* Create the mesh geometry using a meshing package such as IDEAS and export it to a universal file called, for instance, geo-file.unv. * Convert geo-file.unv to a converter file, called cnv, by running program u2c-new2.f. * Convert cnv to the desired format geo-file by running program c2p-fast.f. Both u2cnew2.f and c2p-fast.f are available in the same directory with the main code. 3.3 Preprocessing * In module dim.inc set maxnod greater or equal to the number of nodes in geo-file. * In module dim.inc set nx, ny, nz in accordance to the cylindrical grid structure. nx is the number of grid points along the p direction, ny is the number of grid points along the 4 direction and nz is the number of grid points along the z direction. The number of grid points is determined by the grid steps and the dimensions of the structure. A typical maximum grid step is chosen to be about 0.05A. * Compile the preprocessing code by typing make -f Makestartjet. * Run the preprocessing code by typing startj. * While running startj enter the mesh file name (geo-file) and the AIM near field threshold in wavelengths. A typical value of the threshold is 0.6 wavelengths, but for most geometries decent accuracy can be retained for thresholds as low as 0.3 wavelengths. * The output of startj is the number of near field element pairs (nonzero) and also the FFT orders orx2, ory2, orz2 along the three axes p, q0, z. * It module dim.inc set nonzero, orx2, ory2, orz2 according to the output of startj. * In module param.inc set * in=number of nonnegative modal orders in hollow section (Region 1). *iinumber of modes per order. Chosen to include at least all propagating modes. * ngas=number of points per subinterval for hand-off plane integration. * nsym=number of periodic sectors (slices) of engine termination. 10

nangle=number of different incidence angles (0z,,c). bgas=Gaussian integration tolerance (at the hand-off plane). * exrho=2intrho=-number of integration subintervals along the p direction. * exphi=29ntPhi-=number of integration subintervals along the 4 direction. * zhandoff=z location of the hand-off plane (the tip of the termination is assumed at z =0. *str=grid step along p direction (in wavelengths). Usually set to 0.05A. *stp=grid step along q direction (in radians). Chosen so that the maximum distance between any two grid points is less than a given threshold (e.g. 0.05A). *stz=grid step along z direction (in wavelengths). Usually set to 0.05A. 3.4 Compilation and running * Compile the code by typing make -f Makeaimjet. * Run the code by typing aimj. The interaction is self-explanatory. See section 4 for further details. 3.5 List of subroutines 1. Initialization codes: u2c-new2.f: converts geo-file.unv to a converter file, called cnv. c2p-fast: converts cnv to the desired format geo-file. 2. Subroutines in Makefile: Aimjet.f: main code. ac-intcyl2.f: performs a double surface integration over two triangles of a function in local area coordinates. an-intl.f: performs an analytical double surface integration over the singular part of 11

the integrand. edge-len.f: updates the edge length table given the edge numbers of the source and observation triangles. excitcyl.f: computes elements of excitation vector (right hand side of linear system). prepfill.f: Includes several subroutines associated with AIM calculations. Specifically it includes prepaimcyl.f which calculates the mesh to grid mapping coefficients, green.f which computes the Green's function matrix and its 2D FFT, fillzaimcyl.f which calculates elements of the impedance matrix, nearaim.f which calculates the near field interactions, farraim.f which calculates the far field interactions, and several auxiliary subroutines, i.e. multil, multi2, discrete, regmom, regmomdiv, prodmix, vemult and FOURN. See code listing for further information on the function of each subroutine. FOURN is the Numerical Recipes FFT subroutine, which must be replaced by a parallel scheme for improvement of the code efficiency. bcgcyl.f: Solves the AIM linear system using the Biconjugate Gradient Method (BiCG). The currenrt BiCG subroutine is valid only for symmetric systems, which is not the general case when the periodicity "slicing" scheme is exploited. For general terminations this subroutine must be replaced by a non-symmetric iterative scheme. get-R.f: reads in triangle resistivity values from file specified by user. According to Pamela Haddad, author of the basic MoM code, the algorithm has not been validated for non-metallic scatterers. get-meshl.f: reads in mesh data from file specified by user. get-verts.f: returns the vertices of a triangle. normal.f: constructs normal vector to a triangle. ops.f: various small functions and subroutines performing auxiliary mathematical calculations. pot-ints-ss.f: computes the potential integrals for uniform and linearly varying surface 12

sources distributed on a planar polygon S. Dependent subroutine in the same module: i-params. pwcyl5.f: returns the value of (r-rm) Ei, where Ei is the incident cylindrical mode and rm is the position vector to the mth vertex of the observation triangle. Dependent subroutines in the same module: integration, discreteap, fieldap, vectintcyl, ffargcyl, dyadfree, cGsldx, GAULEGX, cGsldy, GAULEGY. See code listing for more information on the function of each subroutine. resfnc.f: returns the value of the integrand of the resistive term (r-rm) (r-rn), where rm, rn are the position vectors to the mth/nth vertices of the given observation triangle. xsintcyl.f: performs a single surface integration of a function over a triangle in local area coordinates. solve2.f: linear system solvers. update-signs.f: updates the triangle sign and edge use tables for an input observation interior edge and an input source interior edge. usr-datacyl.f: prompts user for, and reads in, data pertaining to file names, incident field, rcs observation cut via standard I/O. vectint.f: performs a single surface integration of vector function over a triangle. inverd.f, invers.f: sets of subroutines for the inversion of linear systems. cylinsubr2.f: includes several subroutines and functions related to the modal analysis of the hollow section (duct). Subroutines included: excitvec, radcrosec, vemult, modesl. See code listing for more information on the function of each subroutine. blockdata.f: initializes common block which contains the local area coordinate integration points and weights used in all integration routines. 4. Auxiliary files: const.inc: includes various constants. 13

dim.inc: includes uniform array dimensioning. param.inc: includes modal and integration parameters. 14

N3 N6 N1 E14 Ell. '.. T5 T2 -\ E5 I\ E2 E13 "\ E4 \ E12 T4 T3 N7 N9 N5....6.,......E3 i \ T6 T8 '. E1 E8 E10 \ E7. E16 T1.. I.T7 N4 N8 \ E N2.E9 -E15 N Figure 6: Mesh with node, edge and triangle numbering. 15

4 Examples This example is of a flat plate termination of radius A. To apply the periodicity slicing scheme only half of the termination is meshed. The length of the hollow region is A and the hand-off plane is located 0.1A in front of the termination. The maximum radial, azimuth or axial distance among grid points is 0.05A, the expansion order is Mi = 3 and the near field threshold is set equal to rthr = 0.7A. For testing purposes, only propagating modes of order 0 are used in the waveguide. To account for a 0.05A grid step along all threee directions, in dim.inc set nx=26, ny=70 and nz=7 (grid point numbers in the p, 4 and z directions respectively). Running startj with a 0.7 wavelength threshold yields nonzero=191355, orx2=64, ory2=256, orz2=16. These values must be set in module dim.inc before code operation. The param.inc file should read integer ii,in,iis,ii2,ii4 integer ngas,nsym,nangle integer intrho,intphi,exrho,exphi real*8 bgas,zhandoff parameter(in== 1,ii=3) parameter(iis=in*ii) parameter(ii2=2*iis) parameter(ii4:=2*ii2) parameter(ngas=3) Parameter (nsym=2) Parameter (nangle=121) parameter(bgas=.Od-2) parameter(intrho=2) 16

parameter(intphi =2) parameter(exrho=92**intrho) parameter(exphi=2**intphi) parameter(zhandoff=O. idO) parameter(pifix=3.141592654d0) parameter(str=0.05d0,strh=0.5d0*str) parameter(stp=pifix/31.OdO,stph=0.5d0*stp) parameter(stz=0.05d0,stzh=0.5d0*stz) We emphasize that in general in and ii are chosen so that all propagating modes are included. In our case we set in and ii very low only for testing purposes. The code output yields the propagation characteristics of each mode, hence running the code itself for a few test cases eventually provides the correct values of in and ii. For this geometry the input reads: Enter inlet radius in wavelengths 1.OdO Enter mesh file name: disk-lw-2 Enter length of hollow section in wavelengths 1.OdO Enter 1 for inclusion of rim contribution and 2 for non-inclusion 2!(not including Ufimtsef currents on the rim) Enter 1 for inclusion of evanescent modes 17

and 2 for exclusion 2 Enter measurement frequency in GHz 0.3d0!(to obtain RCS in dB/,2) Enter near field threshold in wavelengths 0.7d0 The output reads: TM modes (n,m,krho,prop./evan.) 0 1 2.40482565831768 1 0 2 5.52007838590681 1 0 3 8.65372785443314 -1 TE modes (n,m,krho,prop./evan.) mode indices, eigenvalues in A-1 1 for propagation and -1 for evanescence 0 1 3.83170597025498 1 0 2 7.01558679227743 -1 0 3 10.17346810959348 -1 Incident mode order= 0!(working with mode order=O) Order of system = 575 No. of triangular elements = 401 Slice 0!(working with slice 0) 18

Phi-phi Polarization 1.0000000000000E-02.51 -52.7482157063922 1.01 -40.8889918329867 1.51 -33.9207128371362 2.01 -28.9770105579115 2.51 -25.150081887174 3.01 -22.0336614785254 3.51 -19.4104360759395 4.01 -17.150449664679 4.51 -15.1697659720945 5.01 -13.4110150873241 -121.047319669989 RCS vs. 0inc in degrees Theta-theta Polarization 1.OOOOOOOOOOOOOOOE-02 -114.692693414436.51 -46.3946089608898 1.01 -34.5383648650336 1.51 -27.5750249452714 2.01 -22.6382194275045 2.51 -18.8201429361788 3.01 -15.7145271891473 3.51 -13.1040553080986 4.01 -10.8587669470665 4.51 -8.89472076651106 19

5.01 -7.15454102186346 References [1] E. Bleszynski, M. Bleszynski and T. Jaroszewicz, "AIM-Adaptive integral method compression algorithm for solving large scale electromagnetic scattering and radiation problems" Radio Science, vol. 31, no. 5, pp. 1225-1251, Sept.-Oct. 1996. [2] H.T. Anastassiu, Electromagnetic Scattering from Jet Engine Inlets Using Analytical and Fast Integral Equation Methods, Ph. D. Thesis, University of Michigan, 1997. [3] H.T. Anastassiu, User's Manual for MOMJET, University of Michigan, 1997. [4] D. C. Ross, J. L. Volakis and H. T. Anastassiu, "Hybrid Finite Element-Modal analysis of jet engine inlet scattering" IEEE Trans. Antennas and Propagation, vol. 43, no. 3, pp. 277-285, March 1995. [5] D. C. Ross, J. L. Volakis and H. T. Anastassiu, "Overlapping modal and geometric symmetries for computing jet engine inlet scattering" IEEE Trans. Antennas and Propagation, vol. 43, no. 10, pp. 1159-1163, Oct. 1995. [6] 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. 20

USER'S MANUAL FOR MOMJET Departr Hristos T. Anastassiu, Mikhail Smelyanskiy, and John L. Volakis Radiation Laboratory nent of Electrical Engineering and Computer University of Michigan 1301 Beal Avenue Ann Arbor MI 48109-2122 FAX: (313) 647-2106 hristosa@engin.umich.edu msmelyan@engin.umich.edu volakis@engin.umich.edu Science June 18, 1997 1

1 Introduction MIomjet is a FORTRAN7 code which uses the Moment Method (MIoNI) to model electromagnetic scattering from perfectly conducting jet engine models. It is based on a free space Moment Method code written by Pamela Haddad at the Radiation Laboratory in 1991. Section 2 of this manual gives a brief description of MoM applied to the engine problem. Section 3 describes the code operation, and section 4 presents a few examples of Radar Cross Section calculations. 2 MoM applied to the jet engine problem A major reason for an emphasis on the EM characterization of jet engines is their critical importance in determining the airplanes's reflectivity and for target identification. Even though the engine is a small component of the overall structure, its complexity and visibility in the front sector of the aircraft make its contribution a large component of the overall airframe scattering. Moreover, the rotating engine blades cause unavoidable doppler shifts which are important for target identification purposes. The latter can be an important discriminator for both civilian and military applications and more reliable than other methods. As applied to electromagnetic scattering problems, the Moment Method (MoM) has proven to be a very accurate integral equation technique. However, its brute force application to the jet engine problem is not feasible due to excessive electrical size of the structure. To reduce the CPU time and storage requirements down to manageable levels, the inherent blade periodicity of the jet engine was exploited to show that the computational domain can be reduced down to a single blade or engine "slice". This discrete body of revolution (DBOR) approach reduces the number of unknowns and storage requirements by a factor equal to the number of blades. In [1] we applied the DBOR concept in the context of integral equation methods for modeling the complex jet engine configurations. By considering mode by mode excitation [1, 2, 3] it is shown that the analysis over 2

Figure 1: Simplified model of a jet engine. the entire engine can be reduced down to a surface integral equation over a single blade. It is further demonstrated that the resulting modal scattering matrix is sparse, leading to additional storage and CPU time reductions. To demonstrate the reduction of the computational domain we consider a cylindrically periodic scatterer (Fig. 1). We define the current column vector {J}(mS) of the mth slice by {J}(m9) [i(" ), M),.. I()]T (1) where I(m) are the elementtay currents on the mth slice. Brute force MoM yields the linear system 3

[A](00~~){J}(~) + [A](~l){J}(1) +... + [A](~'N-I){J})(-1 ) = {b}(~) [A](1"){J}(O) + [A](11){J}(1) +... + [A](l N' -1){J}('I,-l) = {b}(1 [A] (N-IO){J}(O) +... + [A] (Ns-1Ns-1){J}(NA-1) = {b}' — (2) where A(nq m jkZ J fp(r) R GR.r, R r') R fq(r')d Sd2S' (3) are the entries of the submatrices [A](nm's) representing the interactions between the fp and fq elementary currents located on the mh and the n' slices, G is the cylindrical waveguide dyadic Green's function and R is the rotation dyadic [1]. Also, {b}(n) is the excitation column vector of the th slice defined by {b}(n) - [V(nSv), V2(n..V(ns)T (4) where () = E. (R".r) R fp (r) dS (5) In the latter equation, E' is the incident field, So is the outer surface of the slice within the volume Vo. Clearly, the index ms indicates the slice where the source point is located, whereas ns is the index representing the slice containing the observation point. Both indices run from 0 to Ns - 1. Following the analysis in [1] it can be shown that it is sufficient to determine only the currents on the basic slice by solving [[A](~~) + [A](01)e i'm4 +... + [A]O-(N-1)e3(N-1)n f ] {J}(o) {b}(~) (6) If the cylindrical waveguide dyadic Green's function is to be used, ( 6) can be further simplified, but that procedure is mainly of theoretical importance; the cylindrical waveguide 4

dvadic Green's function is very difficult to handle computationally. and therefore the free space Green's function is used throughout the code. Additional elementary currents are placed on the waveguide walls to account for the appropriate boundary conditions. An additional result of [1] is that for a given order of the incident mode. only a limited set of scattered modes is excited. Namely, only the scattered modes with orders n that satisfy n = ni, + vNs, v E Z are reflected back, and this is in agreement with the result given by the FEM analysis of a similar problem [3]. The importance of the periodicity analysis cannot be exaggerated. Since the problem essentially reduces to modeling only one slice of the scatterer, the number of equations or unknowns is reduced by a factor of NV,. For a typical Ns, = 40, this implies a CPU time and memory reduction each by a factor of 1600. For large scatterers with large periodicity numbers N., such as jet engines, the problem can thus be scaled to a tractable size. Moreover, the limited coupling among the scattered modes results in sparse modal scattering matrices which are much easier to store and handle. 3 Code Operation 3.1 Code structure Momjet models the scattering from metallic cylindrical jet engine inlets terminated by an arbitrary, PEC, cylindrically periodic structure. A simplified geometry is illustrated in Fig. 2, where the termination is a cylindrical hub. A plane wave is assumed to illuminate the open end and couple to waveguide modes in Region 1, propagating towards the termination (Region 2). Each of these modes is used as an excitation for the MoM applied to one periodic sector of the termination (see Figs. 3, 4). Solving the MoM linear system yields the equivalent currents on the termination, and in turn the scattered field within the inlet. To extract the coefficients of the scattered modes propagating towards the open end an integration [1, 2] is performed at a hand-off plane, which is perpendicular to the inlet axis and lies at a fraction of a wavelength in front of the termination (see Fig. 5). Each of the scattered modes propagates to the open end and radiates in free space. The overall modal contribution yields the Radar Cross Section (RCS) of the engine. 5

A A Region I ^ ------- \ --------- 2 --- z=O. Figure 2: Cylindrical inlet terminated by a cylindrical hub To exploit the inherent periodicity of the engine, the waveguide modes must have exponential, and not trigonometric dependence on q. Immediate consequence of this is that the modal orders may be negative, positive or zero, and there is no distinction between "odd" and "even" behavior, as opposed to the Mode Matching code that has been developed in the past [4]. Additionally, this representation enables calculation of RCS modulation (caused by blades rotation) in a simple manner as a postprocessing job [5]. 6

4 - /;: w.. kr:?? *I,:::t % F!- 'igure '.3.' A quateri.ofhugeoetr!~~~~~~~~~~~~~~~ ~ ~-i,,,,-. —...,......*....../,!~ i~' i~~~~~~~~~w' '~~~~~,':,,:,Z,.~L..,~ ~~-*',.~,? '. 1 7

i \\ I ~-,-?yc-t>:iS~* Y \ '7 Figure 4 A quarter of a blade geometry 8

I EMNOWWW" - X2 I I I I I I Engine Modal or Ray region Fields to Ray/Modes Integration Surface (Hand-off Surface) I I - i I Numerical or Rigorous Modeling Figure 5: Decomposition of the computational domain. 3.2 Input file.gd The input file comprises of a discretized sector of the termination, including the waveguide wall up to the hand-off plane. Typical inut geometries are iustrated in Figs. 3, 4. The.. o edge-based forme n te r current version of Momjet requires an edg ae rat the geometry fie (nter nodes and corresponding edges are listed first). To illustrat e the file format e onsider a 0.2A x 0.2A flat plate (see Fig. 6), lying on the xy plane and centered at the origin. This geometry is not typical for the code, but its simplicity is exploited to clarify the meshing format of the input file. The geometry file has the form 9

9 8.1.1.0.1 -.1.0 -.1.1.0 -.1 -.1.0.1.0.0.0.1.0 -.1.0.0.0 -.1.0.0.0.0 16 8 7 8 5 6 5 9 6 9 3 9 7 9 8 9 2 9 4 8 4 7 1 6 1 5 3 7 36 2 8 2 5 1 9 10 2 1112 234 5 6 13 4 5 14 1 6 7 8 15 7 16 8 3 (N\nodes - 9. Ntriangles - 8) (the next Nnodes entries are the cartesian coordinates of the nodesinner nodes first followed by outer (total number of edges Nedges) (number of exterior edges Next) (next Nedges - Next entries give the pair of nodes defining the edges) (the remaining entries give the triangle connectivity. Each line entry contains 3 numbers corresponding to the edge numbers forming the triangle in a clockwise order) 10

N3 N6 N1 E14 Ell N3 N6 Ni T5 T2 E5 E2 E13 E4 E12 T4 T3 N7 N9 N5.................................................... E 6..........................................................E 3............. T6 T8 \ E8 "\ E1 E10 E7 E16 T1 T7 N4 N8 N2 FQ F-l1........... -................................. i............................................................................ _.. Figure 6: Mesh with node, edge and triangle numbering. 11

To obtain this format: * Create the mesh geometry using a meshing package such as IDEAS and export it to a universal file called, for instance, geo-file.unv. * Convert geo-file.unv to a converter file, called cnv, by running program u2c-new2.f. * Convert cnv to the desired format geo-file by running program c2p-fast.f. Both u2cnew2.f and c2p-fast.f are available in the same directory with the main code. Momjet also prompts user for outer radius (in terms of A), inlet length up to the handoff plane (in terms of A), measurement frequency (choose 0.3dO if want RCS in dB/(A2)). Inclusion of Ufimtsef currents and evanescent modes is optional. All input parameters may be specified in the input file. 3.3 Preprocessing Prior to running the code, several parameters must be set appropriately: * In module dim.inc set maxnod greater or equal to the number of nodes in geo-file. * In module param.inc set *in=number of nonnegative modal orders in hollow section (Region 1). *ii=number of modes per order. Chosen to include at least all propagating modes. * ngas=number of points per subinterval for hand-off plane integration. * nsym=number of periodic sectors (slices) of engine termination. * nangle=number of different incidence angles (Oi,c). * bgas=Gaussian integration tolerance (at the hand-off plane). * exrho=2intrh~onumber of integration subintervals along the p direction. * exphi=2i=tphi=number of integration subintervals along the q$ direction. * zhandoff=z location of the hand-off plane (the tip of the termination is assumed at z= 0. *irotmax=max rotational angle 4 (sugested integer value). 12

"irotst=rotation step (sugested integer value). *irott=number of rotation steps. 'kstepp=O step for RCS calculation. *kfica=number of sempling points on the hand-off plane for unique engine ID representation. In order to define accurate number of propagating modes (in and ii), two runs of the code have to be (lone, according to following steps: 1) inicialize in and ii sufficiently high. 2) uncomment line in the Makefile for compilation without optimisation, and comment the line with optimization includded. 3) compile and run the code. 4) check max in and ii when modes are propagating, stop execution, and change parameters in the param.inc file. 5) uncomment line with included optimization and comment line without optimization. 6) compile and run the code. 3.4 Compilation and running * Compile the code by typing make (there is only one Makefile in the directory). * Run the code by typing momj. The interaction is self-explanatory. See section 4 for further details. If you want to use input parameter redirection, then type momj < input.dat. 3.5 List of subroutines 1. Initialization codes: u2c-new2.f: converts geo-file.unv to a converter file, called cnv. 13

c2p-fast: converts cnv to the desired format geo-file. 2. Subroutines in Makefile: Momjet.f: main code. ac-intcyl2.f: performs a double surface integration over two triangles of a function in local area coordinates. an-intl.f: performs an analytical double surface integration over the singular part of the integrand. back-subst.f: solves the MoM linear system using LU decomposition. edge-len.f: updates the edge length table given the edge numbers of the source and observation triangles. excitcyl.f: computes elements of excitation vector (right hand side of linear system). fillzcyl2.f: calculates elements of impedance matrix Z corresponding to incident modal order nine. get-R.f: reads in triangle resistivity values from file specified by user. According to Pamela Haddad, author of the basic MoM code, the algorithm has not been validated for non-metallic scatterers. get-meshl.f: reads in mesh data from file specified by user. get-verts.f: returns the vertices of a triangle. normal.f: constructs normal vector to a triangle. ops.f: various small functions and subroutines performing auxiliary mathematical calculations. pot-ints-ss.f: computes the potential integrals for uniform and linearly varying surface sources distributed on a planar polygon S. Dependent subroutine in the same module: i-params. pwcyl5.f: returns the value of (r-rm) Ei, where Ei is the incident cylindrical mode and 14

rm is the position vector to the nth vertex of the observation triangle. Dependent subroutines in the same module: integration, discreteap, fieldap, vectintcyl. ffargcyl. dyadfree, cGsldx, GAULEGX, cGsldy, GAULEGY. See code listing for more information on the function of each subroutine. resfnc.f: returns the value of the integrand of the resistive term (r-rm) (r-rn). where rm, rn are the position vectors to the mth/nth vertices of the given observation triangle. xsintcyl.f: performs a single surface integration of a function over a triangle in local area coordinates. solve2.f: linear system solvers. update-signs.f: updates the triangle sign and edge use tables for an input observation interior edge and an input source interior edge. usr-datacyl.f: prompts user for, and reads in. data pertaining to file names, incident field, rcs observation cut via standard I/O. vectint.f: performs a single surface integration of vector function over a triangle. cylinsubr2.f: includes several subroutines and functions related to the modal analysis of the hollow section (duct). Subroutines included: excitvec, radcrosec, vemult, modesl. See code listing for more information on the function of each subroutine. blockdata.f: initializes common block which contains the local area coordinate integration points and weights used in all integration routines. 4. Auxiliary files: const.inc: includes various constants. dim.inc: includes uniform array dimensioning. param.inc: includes modal and integration parameters. 5. Output files: rcsphi.txt: RCS = RCS(Ofin) for q polarization. rcsth.txt: RCS = RCS(Oinc) for 0 polarization. 15

modul.txt: RCS = RCS( oot )le,,,=const modulation file. etotal.txt: total unweighted scattered field on the hand-off plane (ID for every engine termination). Need to be weighted in the future. Number and position of sempling points is determined by the parameters ngas, exrho and exphi (see param.inc file and figure 7). 16

Figure 7: Sampling points on the hand-off plane. 17

4 Examples This example is of the 4 blade termination described in Fig. 4. The outer radius is 2A, the hub radius is A, the termination is 0.5A deep, the length of the hollow region is 3A and the hand-off plane is located 0.1A in front of the termination. Since the number of nodes for this geometry is 656, the dim.inc file should read Integer nmax,maxnod,maxedg,maxtri,maxz,maxrcs Parameter (nmax=4) Parameter (maxnod=660) Parameter (maxtri=2*maxnod) Parameter (maxedg=3*maxnod) Parameter (maxz=3*maxnod) Parameter (maxrcs=200) and the param.inc file should read integer ii,in,iis,ii2,ii4 integer ngas,nsym,nangle,kfica,irott,irotmax,itotst integer intrho,intphi,exrho,exphi,kstepp real*8 bgas,zhandoff parameter(in= 11,ii=4) parameter(iis=in*ii) parameter(ii2=2*iis) parameter(ii4=2*ii2) parameter(ngas=3) 18

parameter (nsym=4) parameter (nangle=121) parameter( bgas= 1.Od-2) parameter(intrho=2) parameter( intphi =2) parameter(exrho=2**intrho) parameter(exphi=2P**intphi) parameter(zhandoff=O. idO) parameter(irotmax=91 ) parameter(irotst=l) parameter(irott=irotmax/irotst + 1) parameter(kstepp= 1) parameter(kfica=exrho*exphi*ngas*ngas) We emphasize that in and ii are chosen so that all propagating modes are included. The code output yields the propagation characteristics of each mode, hence running the code itself for a few test cases eventually provides the correct values of in and ii. For this geometry the input reads: Enter inlet radius in wavelengths 2.0d0 Enter mesh file name: blade_2w_4 Enter length of hollow section in wavelengths 3.0d0 19

Enter 1 for inclusion of rim contribution and 2 for non-inclusion 2!(not including Ufimtsef currents on the rim) Enter 1 for inclusion of evanescent modes and 2 for exclusion 2 Enter measurement frequency in GHz 0.3d0!(to obtain RCS in dB/A2) For redirection running, the input file is: 2.0d0 blade_2w_4 3.0d0 2 2 0.3d0 The output reads: TM modes (n,m,krho,prop./evan.) 01 1.20241288350811 1 0 2 2.760039056538 1 20

0 0 1 1 1 1 2 3 4 1 2 3 4 1 4.32686395109136 1 5.89576721702858 1 1.91585298512173 1 3.50779348764358 1 5.08673402355608 1 6.66184574789692 -1 2.56781132204712 1 mode indices. eigenvalues in A-1' 1 for propagation and -1 for evanescence TE modes (n,m,krho,prop./evan.) 0 1 1.91585298512173 1 0 2 3.50779348764358 1 Incident mode order= -10 Order of system = 1777 No. of triangular elements = 1216 Slice 0!(working with slice Output file rcsphi.txt is: Phi-phi Polarization 1.000000000000OOOOOOOOOOOOOOOE-02 29.631303099879.51 29.61263172062007 1.01 29.55830862102566 1.51 29.4690832676185!(working with mode order=-10) 0) }38 RCS vs. 0inc in degrees 21

2.01 29.34619755314 773 2.51 29.19138128112484 3.01 29.00683892922878 3.51 28.7952212610019 4.01 28.55957452047508 4.51 28.30326037655948 5.01 28.02984213788905 5.51 27.74293755636333 Output file rcsth.txt is: Theta-theta Polarization 1.OOOOOOOOOOOOOOOE-02 29.63128937395119.51 29.57655924653133 1.01 29.41605365137127 1.51 29.14821273587283 2.01 28.77040147133607 2.51 28.27883999812705 3.01 27.66851579327062 3.51 26.93309717096325 4.01 26.06489475433391 4.51 25.0549787304352 5.01 23.89369678494705 5.51 22.57214244083259 Output file rnodul.txt is: Phi-phi Polarization 22

Radar Cross Section I~ 35 30 25 201 C) 0 cr 151 \\..................... x...............:............,....:....................... phi-phi polarization theta-theta polarization \\.. 10l 5[ OL ) 10 20 30 Incident angle theta 40 50 60 Figure 8: RCS for both polarizations. 23

theta= 1.000000000000000E-02 1 29.63130309269711 2 29.67204477665604 3 29.672235844,46394 4 29.67223669880159 5 29.67223664691986 Theta-theta Polarization theta= 1.OOOOOOOOOOOOOOOE-02 1 29.6312893811335 2 29.67203110417029 3 29.67222223582765 4 29.67222317841145 5 29.67222323874154 Output file etotal.txt is: X,Y,Z,Unwaited E-Field at the Handoff Plane 5.641662856629827E-02 1.015140670462751E-02.1 81.08614154147517 4.049268616987463E-02 4.057375264402420E-02.1 80.9713012643719 1.003855320833504E-02 5.643681853501174E-02.1 80.26172131340697 -1.015140671619875E-02 5.641662856421617E-02.1 79.68034874500992 -4.057375265232938E-02 4.049268616155283E-02.1 78.92981225914005 24

Phi-phi modulation *J I I --- —I ----III ----I- - --- - - 28 271.............................. ~~~~~~~~~~~~~~~~.o............................. 26' 25 - ~ () 0 n 24 1 - 23 \:. \ \ - ~- I~'~ 221-................. /: / /. ____ I ____! \ 21 - -\ \..... 20 -... -- theta=5 - theta=15. theta=25 I..:........... \. I * * * * * f ~/ /. i. '........ J.. '. Z.. 190 0 * I 10 20 30 40 50 60 70 80 90 Rotation angle 100 Figure 9: Modulation pattern for rotating blades-phi polarization. 25

Theta-theta modulation 30 / -.............. 28/.,' '......:.: 26. 24.... 1 8 -......................................................... 22eta / \.:. / '.. 18 \: / - - theta=15 | X, \ / "\ A - theta=15 \ ' l L - 1 theta=25 16 0 10 20 30 40 50 60 70 80 90 100 Rotation angle Figure 10: Modulation pattern for rotating blades-theta polarization. 26

References [1] H.T. Anastassiu, Electromagnetic Scattering from Jet Engine Inlets Using Analytical and Fast Integral Equation Methods, Ph. D. Thesis, University of Michigan, 1997. [2] D. C. Ross, J. L. Volakis and H. T. Anastassiu, "Hybrid Finite Element-Modal analysis of jet engine inlet scattering" IEEE Trans. Antennas and Propagation, vol. 43, no. 3, pp. 277-285, March 1995. [3] D. C. Ross, J. L. Volakis and H. T. Anastassiu, "Overlapping modal and geometric symmetries for computing jet engine inlet scattering" IEEE Trans. Antennas and Propagation, vol. 43, no. 10, pp. 1159-1163, Oct, 1995. [4] 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. [5] C. Ross, J. L. Volakis and H. T. Anastassiu, "Efficient Computation of Radar Scattering Modulation from Jet Engines" Radio Science, vol. 31, no. 4, pp. 991-997, July-Aug., 1996. 27