c
TU Delft, The Netherlands, 2006
FULL WAVE ANALYSIS OF THE INFLUENCE OF THE JET
ENGINE AIR INTAKE ON THE RADAR SIGNATURE OF
MODERN FIGHTER AIRCRAFT
1Duncan R. van der Heul, Harmen van der Ven, and Jan-Willem van der Burg
National Aerospace Laboratory NLR Aerospace Vehicles Department Section Flight Physics and Loads
Anthony Fokkerweg 2, P.O. Box 90502, 1006 BM, Amsterdam, The Netherlands e-mail: vdheuldr@nlr.nl
web page: http://www.nlr.nl
Key words: Computational Electromagnetics, Radar cross section, cavity scattering, higher order edge elements
Abstract. Radar cross section prediction techniques are used to determine the radar signature of a military platform when the radar signature can not be determined experi-mentally, because the platform is not available or for reasons of time and cost. For classic jet aircraft the radar cross section for forward observation angles is dominated by the contribution of the open ended cavity formed by the jet engine air intake and compres-sor fan. This cavity is characterized by its large depth (L/d > 3), curved centerline and nonuniform cross section, for which the scattering characteristics can not by analyzed by approximate high frequency methods. Jin et al. have published a numerical method based on a higher order finite element discretisation of the Maxwell equations, where the re-sulting linear system is solved by means of a frontal solution method. The method takes full advantage of the topology of the cavity scattering problem and has been successfully applied for the analysis of cavities of intermediate size. In this paper an adaptation of their algorithm is discussed that can efficiently compute the electric field scattered by very large cavities, in particular the jet engine air intake cavity for X-band radar frequencies.
1 INTRODUCTION
The detectability of a weapon system by radar, often expressed in the quantity radar cross section, is a key parameter in assessing the possible threats imposed by hostile platforms and systems. The radar cross section can be determined experimentally, either for a scale model or for full scale and in-flight conditions. Full-scale and in-flight testing is very appropriate if the platform is available for measurements. However, often there is an
1
interest from a tactical point of view to gain insight in the detectability of weapon systems that are not available for measurements. In these situations use is made of theoretical and computational methods to determine the radar cross section.
A first order analysis of the radar cross section of a geometrically complicated object can be made by approximating the geometry by a set of simple geometrical shapes ( cylinders, balls, ogives and/or cylindrical cavities) for which the radar cross section has been assessed theoretically or experimentally. The contributions of all comprising parts are added to provide the radar cross section of the complicated object, often with a (surprisingly) high accuracy, because at relatively high frequencies the interaction between the comprising parts is small. When this analysis is performed for a jet fighter aircraft at X-band frequency (Lynch18), it reveals that the radar cross section for forward directions of observation is severely dominated by the contribution of the jet engine air intake. The deep cavity formed by the air intake and the compressor fan of the jet engine contributes not just for nose-on observation but for a large range of observation angles centered around this direction.
2 A MODEL FOR JET ENGINE AIR INTAKE SCATTERING
of much larger electrical size, i.e. a complete jet engine air intake at X-band frequen-cies on a relatively large computational platform. The nondimensional cross section of a generic fighter aircraft is of the order of 400λ2, whereas the depth to diameter ratio is approximately 10.
3 GOVERNING EQUATIONS
Starting from the time-harmonic Maxwell equations, a vector wave equation can be derived for the electric field (phasor) E. On the nearly perfectly conducting walls and bottom of the cavity a homogeneous Dirichlet boundary condition is imposed for the tangential electric field. On the aperture an integral radiation boundary condition is imposed. De derivation of the continuous equations is extensively discussed by Jin14 and Volakis et al.16 in the context of the analysis of cavity backed antennas.
The electric field inside the cavity E(r) is found by locating the stationary point of the functional (Liu et al.9):
F = 1 2 Z Z Z Vc 1 µr
(∇ × E(r))(∇ × E(r)) − ko2rE(r) · E(r) dV − (1) k2o Z Z S M(r) · Z Z Sa M(r)G(r, r’)dS0 dS + Z Z Sa ∇ · M (r) Z Z Sa G(r, r0)∇0· M (r0)dS0 dS − 2ik0Z0 Z Z Sa M(r) · Hi(r)dS, where µr and r denote the relative permeability and permittivity respectively, Vcand Sa the volume of the cavity and the surface of the aperture respectively. Furthermore, k0 and Z0 denote the wavenumber and the free space impedance. G(r, r0) is the free space Green’s function defined as:
G(r, r0) = e
−iko(r−r0)
4π(r − r0), (2)
and M (r) is the equivalent magnetic surface current on the aperture plane defined as: ˆ
n× E(r), (3)
where ˆn is the unit outward normal vector on the aperture surface. The electric field E(r) and the magnetic current are expanded in a series of test functions with compact support on tetrahedral and triangular elements, respectively:
E(r) = Ninternal X i=1 EiΩi(r) + Ninternal+Naperture X i=Ninternal+1 EiΩi(r)|r ∈ Vc, (4) M(r) = Naperture X i=1 Ei+NinternalΛi(r)|r ∈ Sa,
Substitution of the expansion (4) in the functional (1) and application of the Rayleigh-Ritz procedure leads to a system of linear equations for the unknown complex coefficients Ei in the expansion. The characteristic structure of the resulting finite-element/boundary integral method system matrix is the following:
Aii Aia Aaa Aaa+ M Ei Ea = 0 b , (5)
where Ei and Ea now denote the vector of degrees of freedom within the cavity and those on the aperture, respectively. Matrices Akk are the discretisation of the vector wave equation and sparsely populated, while M is the discretisation of the integral boundary condition and fully populated. Because the cavity is excited on the aperture surface, only the equations for the degrees of freedom on the aperture surface have a nonzero right hand side. In the next section the particulars of the test functions in the expansion and the computation of the element matrices are discussed.
4 HIGHER ORDER DISCRETISATION
4.1 Vector wave equation
As extensively discussed by Jin et al.7,8,9,11 it is essential to use a higher order discreti-sation of the vector wave equation, because for deep cavities the large dispersion error introduced by the standard zero order edge based elements makes these prohibitively inefficient to use.
The functions Ωi utilized for the expansion of the electric field in the weak form of the vector wave equations are higher order curl-conforming interpolatory vector test functions that were introduced by Graglia et al.3. In order to improve the efficiency of the algorithm, rectilinear elements have been chosen here. The zeroth-order curl conforming test functions are given by:
Ωγβ(r) = ξn∇ξm− ξm∇ξn, (6)
where {γ, β, m, n} are the six even permutations of {1, 2, 3, 4} such that γ < β. Ωγβ(r) is associated to the edge that connects face γ to face β of the tetrahedron. The curl conforming higher order basis functions of order p can be expressed as:
Ωγβijkl= Nγβ(p + 2)2ξγξβαˆijkl(ξ) iγiβ
Ωγβ(r), (7)
where iγ, β is taken to be i, j, k, l for γ, β = 1, 2, 3, 4 respectively. The Silvester-Lagrange interpolation polynomial ˆαijkl(ξ) is defined as:
ˆ
where the Shifted Sylvester polynomials ˆRi of order p are defined as ˆ Ri(p, ξ) = 1 (i−1)! Qi−1 k=1(pξ − k), 2 ≤ i ≤ p + 1 1, i = 1. limi→0 ˆ Ri(p,ξ) i ≡ 1 pξ. (9)
The normalization factor Nγβ is defined as: Nγβ = p + 2 p + 2 − iγ − iβ |lγβ|, (10) where lγβ = J (∇ξj× ∇ξk), J = li· lj× lk, li = ∂r ∂ξi . (11)
Because the elements are rectilinear the element matrix of the vector wave equation can be evaluated very efficiently. For a zeroth, first or second order discretisation there are 6, 21 or 45 test functions in each element.
One issue not discussed by Graglia et al.3 is how to ensure compatibility of the discreti-sation in neighboring elements while preserving linear independence of the basis functions. For each interpolation point on the faces of the element three degrees of freedom, each associated with one of the three edges are defined. However, the tangential field on the face has two independent components. Hence, one of the degrees of freedom has to be discarded. Obviously, in two neighboring elements identical degrees of freedom have to be discarded.
In the grid all edges of the elements have unique consecutive integer labels. The convention chosen here is that in each surface interpolation point the degree of freedom is discarded which is associated with the edge with the lowest label.
For the internal interpolation points six degrees of freedom, each associated with one of the six edges, are defined. Since they constitute a vector field with three independent components, three of them depend linearly on the other three, and can be discarded. Since the internal points are not connected to any neighboring element, these three degrees of freedom can be chosen arbitrarily.
4.2 Magnetic field integral equation
The boundary integral in (1) is discretized with test functions that conform to the test functions used to the discretisation of the vector wave equation according to (4). The zeroth order test functions are given by:
Λβ(r) = −ˆn× Ωβ = (12)
1
J (ξβ+1lβ−1− ξβ−1lβ+1) , β = 1, 2, 3,
where lβ are the edge vectors of the triangles. These test functions are equivalent to the standard Rao Glisson Wilton15 test functions commonly used for the discretisation of integral equations that describe scattering by (nearly) perfectly conducting bodies. The higher order interpolatory divergence-conforming test functions are given by:
Λβijk(r) = −ˆn× Ωβijk(r) = N β(p + 2)ξβαˆijk(ξ) iβ Λβ(r), (13) where ˆ αijk(ξ) = ˆRi(p + 2, ξ1) ˆRj(p + 2, ξ2) ˆRk(p + 2, ξ3), (14) and the normalization constant Nβ is defined as:
Nβ = p + 2 p + 2 − iβ
lβ. (15)
Using the expansion of the basis functions in the unitary basis vectors (Graglia et al.3), the entries to the linear system can be evaluated in a straightforward way. The evaluation of the regular integrals is performed using Gauß quadrature rules, while the singular integrals that result when the observation point and the source point coincide are evaluated using Duffy’s transformation.
5 SOLVING THE LINEAR SYSTEM
5.1 State of the art frontal solver
Jin et al.7,8,9,11 propose the use of a frontal solver, based on the classic algorithm by Irons20. To be able to factorize the system matrix in core, without having to store it completely, the algorithm uses a simultaneous factorization and assembly strategy. During this process element matrices are added to the system matrix while simultaneously degrees of freedom are eliminated from the instantaneous system matrix, also known as the frontal matrix. However, the algorithm can not be vectorized easily and does not exploit any sparsity in the frontal matrix.
In this paper, the frontal solvers as implemented in the HSL libraries19 are used. More specifically, the routine MA42, resp. ME42, is used to solve the resulting real, resp. complex, system. The routines are very easy to use, and being supplied as source, can be modified to suit specific purposes.
The computational complexity of a frontal solver is proportional to N · N2
f, where N is the number of degrees of freedom, and Nf is the number of degrees of freedom in the front. So it is of the utmost importance to keep the number of degrees of freedom in the front to a minimum. For cavities of nearly uniform cross sectional area, the number of degrees of freedom in the front is nearly constant, but dependent on the ordering strategy. This ordering process can be interpreted geometrically: for a given mesh, the first so many elements should fill the bottom of the cavity toward the aperture, in such a way that the number of exposed tetrahedral faces is minimal. How to achieve the optimal ordering is discussed in the next section.
Jin et al.7 discuss that straightforward application of the higher order tetrahedral elements would increase problem complexity with respect to a discretisation using mixed order prism elements. Their answer is their so-called element-based algorithm, where the internal degrees of freedom of the elements are eliminated before assembly of the element. The MA42 and ME42 algorithms accept the complete element matrices, and automatically apply static condensations to eliminate the internal degrees of freedom.
5.2 Adaptations
The frontal solver routines have been modified for the specific application in mind. For general applications the frontal solver requires the complete storage of the upper triangular matrix in order to solve for all equations. In the current application, however, we are only interested in the solution at the aperture. The degrees of freedom at the aperture are the last in the list of degrees of freedom, and to solve them only the lower square of the upper triangular matrix is needed. So the solver has been modified such that the upper triangular matrix is discarded but for the last part. This significantly reduces storage requirements.
as real multiplications, it is worthwhile to exploit the fact that the internal element matrices are real valued. Note that the equations at the aperture are always complex valued, through the appearance of the Green’s function. A solution is to reduce the internal real matrix up to the aperture, using the frontal solution algorithm. Then the assembled matrix, only containing equations for the aperture degrees of freedom, is added to the boundary integral part. The resulting full matrix is solved using a standard dense-matrix LU factorization algorithm. The HSL routines have been modified in such a way that the elimination process is stopped at the aperture and the matrix is dumped to file. Jet engine air intakes are fitted with radar absorbing linings, consisting of multiple layers of different materials with complex valued permittivity or permeability. Hence, the special case of a cavity with lossless material is only of interest for evaluating the effectiveness of such a coating.
6 GRID GENERATION
Use is made of the FASTFLO grid generation system, developed originally for the generation of hybrid grids for computational fluid dynamics, see Van der Burg12,13.
Figure 1: Element ordering by minimizing distance of element centers to bottom of cavity leads to a ’rugged’ surface with an unnecessary large number of degrees of freedom in the frontal matrix.
6.1 Ordering the elements
Figure 2: Element ordering by stacking sets of three elements into single prisms minimizes the number of degrees of freedom in the frontal matrix to nearly the number of degrees of freedom on the aperture surface.
the tetrahedral elements are added one-by-one to the instantaneous grid. All degrees of freedom not located on the exposed surface can be eliminated. Degrees of freedom on the circumferential surface and on the bottom of the cavity are eliminated from the element matrices prior to assembly. Obviously, the instantaneous number of degrees of freedom in the frontal matrix is directly proportional to the instantaneous number of triangular faces exposed in the front. The order of assembly of the element matrices follows the order of the labels of the elements in the grid. Hence, the order of the elements in the grid should be chosen in such a way the maximum number of instantaneous number of triangular faces in the front is minimized.
Different strategies have been used to minimize the maximum number of degrees of freedom in the front. A straight forward ordering strategy is to order the elements by ascending distances of the center of gravity of the elements to the bottom of the cavity. The ordering can be accomplished efficiently, when the elements are put in an octal tree data structure. This type of ordering will be referred to as the Tetris ordering, Additionally, the HSL library19 contains two re-ordering algorithms: a direct and indirect reordering algorithm . Newer versions of the HSL library also contain a spectral reordering algorithm. Application of the Tetris ordering algorithm, the direct or the indirect HSL reordering strategy did not result in a significant difference in minimizing the maximum instantaneous number of degrees of freedom in the front
lead to a maximum number of degrees of freedom in the front close to the number of degrees of freedom located on the aperture plane. Neither of the different strategies used was able to accomplish this. Visual inspection of the front shows a ’rugged’ instantaneous front surface.
To decrease the number of degrees of freedom in the front, a different grid generation strategy was followed. Instead of generating a grid of tetrahedral elements inside the cavity geometry, a grid was generated of layers of prisms parallel to the aperture plane. Next each prism is divided into three tetrahedral elements. The triangles on the cavity bottom are ordered starting from a triangle adjacent to the cavity wall. Next a two dimensional Tetris ordering is applied to the triangles in a arbitrary direction transversal to the centerline of the cavity. The ordering of the triangles is used to order the prism elements in each layer and the layers are ordered starting from the bottom of the cavity toward the aperture. Figure 1 and Figure 2 show snapshots of the instantaneous frontal surface for both meshing and ordering strategies for a rectangular cavity. For this specific case, the maximum occurring number of degrees of freedom in the front for the first strategy is 4100, but for the alternative approach it is reduced to 2300. This approach is not restricted to straight cavities. The grid generation algorithm is able to adapt the prismatic layers to conform to the cavity shape, as illustrated in Figure 3 and Figure 4 for a section of jet engine air intake. For more complicated geometries the divided prism elements can be combined with local clusters of tetrahedral elements. Furthermore, different prismatic grids can be merged in the case of cavities with a strong variation in cross sectional area or multiple disjunct apertures.
6.2 Radar absorbing linings
To reduce the electromagnetic power scattered by the jet engine air intake, a radar absorbing lining can be applied. Thin linings with a thickness tlining λ can be modelled by introducing impedance surface elements, as for instance discussed by Volakis et al.16. Coatings with a thickness tlining≈ λ can be modelled by discretising the volume occupied by the coating and prescribing the appropriate material properties of the coating to those elements.
X Z
Y
Figure 3: Layers of prismatic elements are stacked from the bottom of the cavity, before each prismatic element is subdivided in three tetrahedral elements. For this jet engine air intake geometry the prismatic layers smoothly conform to the continuously changing boundary of the cavity.
7 RESULTS
All computations are performed on a NEC SX-5 computer, on a single processor for the smaller test cases and on 4 processors for the larger test cases. We used a number of test cases presented by Liu et al.9 for validation of the algorithm. One example is the small cavity presented in Figure 8. For this small test case the elements are ordered by ascending distance from the bottom of the cavity. Figure 9 and Figure 10 show the dimensionless radar cross section for a fixed mesh and for zeroth, first and second order of accuracy of the discretisation. The results conform closely to those presented by Liu et al.9. For this cavity, the second order solution is converged for h = λ/2. In all following results a second order of accuracy discretisation is used. A second, larger example is shown in Figure 11. In Figure 12 the dimensionless radar cross section for a vertically polarized excitation computed for two different mesh widths is compared with a reference solution obtained using a modal expansion method. Next the algorithm was applied to the analysis of an even larger circular cavity of diameter D = 20λ and depth L = 40λ (Figure 13). Ling et al.6 use a combination of ray tracing/Physical Optics for the analysis of the same geometry.
X Y
Z
Figure 4: A snapshot of the front that shows how a near optimal ordering of the elements has been achieved.
computed for a mesh width of h = 0.6λ compared with a reference solution obtained using a modal expansion method. The accuracy of the computed solution is comparable to the results presented by Ling et al.6. The slight mislocation of the extrema is a result of the numerical dispersion error, that can be reduced by further decreasing the mesh width of the discretisation.
As a final example a section of a generic uncoated jet fighter engine air intake is analyzed. The geometry is shown in Figure 3. The cavity is approximately 70λ deep with a cross sectional area of 600λ2 for a frequency of 10 GHz. The computational mesh contains 430000 elements, resulting in 8000000 degrees of freedom of which 17500 are located on the aperture surface, The maximum number of degrees of freedom in the front matrix is approximately 18500. The radar cross section of the generic jet engine air intake section is shown in Figure 15.
advan-X Y
Z
Figure 5: Surface triangulation of a jet engine air intake with structured surface triangles. The coating is modelled by a two-layer prismatic grid.
tageous to further increase the order of the discretisation from second to third or fourth order.
8 CONCLUSIONS
X Y
Z
Figure 6: The transposed prismatic grid of Figure 5 having the same grid nodes.
REFERENCES
[1] O.C. Zienkiewicz and R.L. Taylor. The finite element method, McGraw Hill, Vol. I., (1989), Vol. II., (1991).
[2] S. Idelsohn and E. O˜nate. Finite element and finite volumes. Two good friends. Int. J. Num. Meth. Engng., 37, 3323–3341, (1994).
[3] R.D. Graglia, D.R. Wilton and A.F. Peterson, Higher Order Interpolatory Vector Bases for Computational Electromagnetics, IEEE Transactions on Antennas and Propagation, Ap-45(3),329-342, (1997).
[4] H.T. Anastassiu. A Review of Electromagnetic Scattering Analysis for Inlets, Cavities and Open Ducts,IEEE Transactions on Antennas and Propagation, Ap-45(6), 27-40,(2003).
[5] A. Altintas, P.H. Pathak and M.C. Liang, A Selective Modal Scheme for the Anal-ysis of EM Coupling Into or Radiation from Large Open-Ended Waveguides, IEEE Transactions on Antennas and Propagation, Ap-36(1), 84-96, (1988).
X Y
Z
Figure 7: Final tetrahedral grid of a jet engine air intake fitted with a radar absorbing lining with different material markers for the lining and the empty inner part
[7] J-M. Jin, J. Liu, Z. Lou and C.S.T. Liang, A Fully High-Order Finite-Element Sim-ulation of Scattering by Deep Cavities, IEEE Transactions on Antennas and Propa-gation, Ap-51(9), 2420-2429, (2003).
[8] J. Jin, Electromagnetic scattering from large, deep and arbitrarily shaped open cav-ities, Electromagnetics, 18, 3-34, (1998).
[9] J. Liu and J-M. Jin, A Special Higher Order Finite-Element Method for Scattering by Deep Cavities, IEEE Transactions on Antennas and Propagation, Ap-48(5), 694-703, (2000).
[10] A. Barka, P. Soudais and D. Volpert, Scattering from 3-D Cavities with a Plug and Play Numerical Scheme Combining IE, PDE, and Modal Techniques, IEEE Transactions on Antennas and Propagation, Ap-48(5), 704-712, (2000).
[11] J. Liu and J-M. Jin, Scattering Analysis of a Large Body With Deep Cavities, IEEE Transactions on Antennas and Propagation, Ap-51(6), 1157-1167, (2003).
[12] J.W. van den Burg, Tetrahedral grid optimisation: towards a structured tetrahedral grid, nlr Technical Report nlr-tp-2000-343, (2000).
x -0.75 -0.5 -0.25 0 0.25 0.5 0.75 y -0.75 -0.5 -0.25 0 0.25 0.5 0.75 z -0.6 -0.4 -0.2 0 X Y Z
Figure 8: Computational grid for the 1.5λ × 1.5λ × 0.6λ rectangular cavity analyzed by Liu et al.9
.
Conference on Numerical Grid Generation in Computational Field Simulations, San Jose, U.S.A.,(2005).
[14] J. Jin, The Finite Element Method in Electromagnetics, (1993).
[15] S.M. Rao, D.R. Wilton and A.W. Glisson, Electromagnetic scattering by surfaces of arbitrary shape, IEEE Transactions on antennas and propagation, AP-30(3), 409-418, (1982).
[16] J.L. Volakis and A. Chatterjee and L.C. Kempel, Finite Element Method for Elec-tromagnetics, Antennas, Microwave circuits and scattering applications, IEEE/OUP Series on Electromagnetic Wave Theory, IEEE Press, (1998).
[17] J.L. Karty and J.M. Roedder and S.D. Alspach, CAVERN: A Prediction Code for Cavity Electromagnetic Analysis, IEEE Antennas and Propagation Magazine, 37(3), 68-72, (1995).
[18] Lynch, David Jr., Introduction to RF Stealth, SciTech Publishing, (2004). [19] Harwell Software Library 2004, http://www.hsl-library.com.
θ σφφ / λ 2 0 20 40 60 80 -20 -15 -10 -5 0 5 10 15 20 σφφ0-order, 88 dof σφφ1-order 512 dof σφφ2-order 1542 dof Rectangular cavity l x b x d = 1.5λx 1.5λx 0.6λ h=λ/2
Figure 9: Dimensionless radar cross section of 1.5λ × 1.5λ × 0.6λ rectangular cavity for horizontally polarized excitation, using zeroth, first and second order of accuracy discretisation on a fixed mesh.
θ σθθ / λ 2 0 20 40 60 80 -20 -15 -10 -5 0 5 10 15 20 σθθ0-order 88 dof σθθ1-order 512 dof σθθ2-order 1542 dof Rectangular cavity l x b x d = 1.5λx 1.5λx 0.6λ h=λ/2
x -2.5 0 2.5 y -2.5 0 2.5 z -10 -5 0 X Y Z
Figure 11: 5λ × 5λ × 10λ rectangular cavity analyzed by Liu et al.9
. The elements are ordered with ascending distance from the bottom of the cavity.
θ σθθ / λ 2 d B 0 10 20 30 40 20 25 30 35 h=λ/2 h=λ/3 modal solution
z -10 -5 0 5 10 y -10 -5 0 5 10 x -40 -30 -20 -10 0 X Y Z
Figure 13: Large circular cavity with diameter of 20λ and depth 40λ.
θ σθθ / λ 2 d B 0 20 40 20 25 30 35 40 45 50 55 60 65 70 σθθComputed
σθθmodal solution from Ling et al.
φ σ d B m 2 -30 -20 -10 0 10 20 30 σff2-order 8000000 dof
Figure 15: Radar cross section of a generic jet engine air intake section of approximately 70λ deep and with a cross sectional area of 600λ2
