Ching Y. Loh and Ananda Himansu

Taitech, Inc., Beaver Creek, Ohio

Lennart S. Hultgren

Glenn Research Center, Cleveland, Ohio

A 3-D CE/SE Navier-Stokes Solver

With Unstructured Hexahedral Grid

for Computation of Near Field Jet

Screech Noise

NASA/TM—2003-212314

June 2003

AIAA–2003–3207

https://ntrs.nasa.gov/search.jsp?R=20030064175 2018-05-15T18:32:06+00:00Z
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA’s counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
Ching Y. Loh and Ananda Himansu

Taitech, Inc., Beaver Creek, Ohio

Lennart S. Hultgren

Glenn Research Center, Cleveland, Ohio

A 3-D CE/SE Navier-Stokes Solver

With Unstructured Hexahedral Grid

for Computation of Near Field Jet

Screech Noise

NASA/TM—2003-212314

June 2003

National Aeronautics and

Space Administration

Glenn Research Center

Prepared for the

Nineth Aeroacoustics Conference and Exhibit

cosponsored by the American Institute of Aeronautics and Astronautics

and the Confederation of European Aerospace Societies

Hilton Head, South Carolina, May 12–14, 2003

AIAA–2003–3207
Acknowledgments

This work received support from the Supersonic Propulsion Technology Project Office at

NASA Glenn Research Center.

Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

This report contains preliminary

findings, subject to revision as

analysis proceeds.

Available electronically at http://gltrs.grc.nasa.gov

The Propulsion and Power Program at

NASA Glenn Research Center sponsored this work.
Abstract

A 3-D space-time CE/SE Navier-Stokes solver

using an unstructured hexahedral grid is described and applied to a circular jet screech

noise computation. The present numerical results for an underexpanded jet, corresponding

to a fully expanded Mach number of 1.42, capture the dominant and nonaxisymmetric ‘B’

screech mode and are generally in good agreement with existing experiments.

1 Introduction

Jet noise is one of the most important topics in computational aeroacoustics. Many of its aspects are of

primary practical importance and the associated complicated physical phenomena are the topic of many experimental and theoretical investigations. Refs. [1–4]

provide a comprehensive discussion and further references. An over/under-expanded supersonic jet emits

mixing noise, broadband shock-associated noise, as well

as screech tones under certain conditions. The mixing

noise is directly associated with large-scale structures,

or instability waves, in the jet shear layer; whereas, the

broadband shock-associated noise and screech tones are

associated with the interaction of these waves with the

shock-cell structure in the jet core. The screech tones

arise due to a feedback loop, i.e., some of the acoustic waves generated by the wave/shock-cell interaction

propagate upstream and regenerate the instability waves

at, or in the vicinity of, the nozzle lip. The feedback loop

leading to screech tones is sensitive to small changes in

the system conditions, and the understanding of the phenomena is to date based mostly on experimental observations [5–8].

Although substantial progress in numerical computation of jet mixing noise has been made, reliable direct

numerical simulation of jet screech noise has up to quite

recently not been feasible. Shen and Tam [9, 10] obtained excellent results in direct numerical simulations

of screech for circular jets using the well-known DRP

scheme. In their 3-D computation [10], a spectral method

was adopted in the azimuthal direction, and by using only

a limited number of spectral functions substantial savings of computer memory and CPU time was achieved,

without deterioration of accuracy. Other recent computational work includes [11].

The present authors [12–14] successfully computed

axisymmetric near-field screech noise for round jets using the recent CE/SE (space-time conservation-element

and solution-element) method utilizing (unstructured)

triangulated grids. Because of the implementation (based

on flux balance) of the non-reflecting boundary conditions (NRBC), a much smaller near field computational domain can be used with this method. However,

when the CE/SE method is applied to a 3-D rectangular jet screech noise computation [15], a major challenge

emerges. Namely, the number of unstructured tetrahedral

cells required to achieve a certain resolution could be as

high as tens of millions, which is currently somewhat beyond the capability of common parallel computers such

as Linux PC clusters. An alternative is to employ an unstructured hexahedral grid. The number of cells may then

be greatly reduced without much loss of resolution. The

3-D CE/SE Navier-Stokes (N-S) solver needs to be modified to accommodate the unstructured hexahedral grid,

however.

A general description of the CE/SE Euler method with

a hexahedral grid was first given by Zhang et al [16]. In

the present paper, the current 3-D CE/SE N-S solver using an unstructured hexahedral grid is described in §2.

The parallel computation implementation of the scheme

is outlined in §3. A 3-D circular jet screech problem

is described in §4. The numerical results are presented

and compared with the available experimental data in §5.

Concluding remarks are presented in §6.

NASA/TM—2003-212314 1

∗Member AIAA

†Associate Fellow AIAA

A 3-D CE/SE NAVIER-STOKES SOLVER WITH UNSTRUCTURED HEXAHEDRAL

GRID FOR COMPUTATION OF NEAR FIELD JET SCREECH NOISE

Ching Y. Loh* and Ananda Himansu

Taitech, Inc.

Beaver Creek, Ohio 45430

Lennart S. Hultgren†

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135
2 The Numerical Scheme

As stated in §1, for fully 3-D computations of complicated problems such as the jet noise problem, a large

number of computational cells are needed to provide sufficient numerical resolution. This imposes severe memory and speed requirements on the computer system. Of

course, the choice of unstructured grid type to be used in

the computation greatly affects these demands. It turns

out that the use of an unstructured hexahedral grid in the

CE/SE scheme rather than the heretofore standard tetrahedral one leads to significant reductions of the memory

and CPU time requirements. Its drawback, of course, is

a modest increase in dissipation in the scheme. Zhang et

al [16] first presented a 3-D CE/SE Euler method using

hexahedral cells instead of tetrahedral ones and provided

a detailed description of the implementation.

In general, the CE/SE method [17, 18] systematically

solves a set of integral equations derived directly from

the physical conservation laws and naturally captures

shocks and other discontinuities in the flow. In order to

have a compact cell stencil, both conservative variables

and their derivatives are computed simultaneously as unknowns. A brief sketch of the 3-D CE/SE scheme with

hexahedral cells is given below.

2.1 Conservation Form of the 3-D Unsteady

Compressible Navier-Stokes Equations

Consider a dimensionless conservation form of the unsteady 3-D Navier-Stokes equations of a perfect gas. Let

ρ, u, v, w, p, and γ be the density, streamwise transverse and spanwise velocities, static pressure, and constant specific heat ratio, respectively. The 3-D NavierStokes equations then can be written in the following

vector form:

Ut + Fx + Gy + Hz = 0, (1)

where x, y, and z are the streamwise, transverse, and

spanwise coordinates, and t is time. The five components

of the conservative flow variable vector U are given by

U1 = ρ, U2 = ρu, U3 = ρv, U4 = ρw,

U5 = p/(γ − 1) + ρ(u2 + v2 + w2)/2.

The flux vectors in the x, y, and z directions, F, G, and

H, respectively, are further split into inviscid and viscous fluxes,

F = Fi − Fv, G = Gi − Gv, H = Hi − Hv,

where the inviscid fluxes are

Fi1 = U2,

Fi2 = (γ−1)U5+



(3 − γ)U2

2 − (γ − 1)(U23 + U24 )



/2U1,

Fi3 = U2U3/U1, Fi4 = U2U4/U1,

Fi5 = γU2U5/U1 − (γ − 1)U2



U2

2 + U23 + U24



/2U2

1 ;

Gi1 = U3, Gi2 = U2U3/U1,

Gi3 = (γ−1)U5+



(3 − γ)U2

3 − (γ − 1)(U22 + U24 )



/2U1,

Gi4 = U3U4/U1,

Gi5 = γU3U5/U1 −(γ −1)U3



U2

2 + U23 + U24 )



/2U2

1 ;

Hi1 = U4, Hi2 = U2U4/U1, Hi3 = U3U4/U1,

Hi4 = (γ−1)U5+



(3 − γ)U2

4 − (γ − 1)(U22 + U23 )



/2U1,

Hi5 = γU4U5/U1 −(γ −1)U4



U2

2 + U23 + U24 )



/2U2

1 ;

and the viscous fluxes are

Fv1 = 0, Fv2 = μ(2ux − 2

3∆),

Fv3 = μ(uy + vx), Fv4 = μ(uz + wx),

Fv5 = μ[2uux + v(uy + vx) + w(wx + uz) − 2

3u∆+

γ

P r

∂

∂y (

U4

U1

− u2 + v2

2 )],

Gv1 = 0, Gv2 = μ(vx + uy),

Gv3 = μ(2vy − 2

3∆),

Gv4 = μ(vz + wy),

Gv5 = μ[2vvy + u(uy + vx) + w(wy + vx)w − 2

3 v∆+

γ

P r

∂

∂y (

U4

U1

− u2 + v2

2 )],

Hv1 = 0, Hv2 = μ(wx + uz),

Hv3 = μ(vz + wy), Hv4 = μ(2wz − 2

3∆),

Hv5 = μ[2wwz + u(uz + wx) + v(wy + vx) − 2

3 v∆+

γ

P r

∂

∂y (

U4

U1

− u2 + v2

2 )],

where u, v, w, ux, uy, uz, etc. are the flow velocities

and their spatial derivatives, which can be expressed in

terms of the conservative variables U1, U2, etc. along

with their gradients. P r = 0.72 is the Prandtl number, μ

the kinematic viscosity, and

∆ = ux + vy + wz

is the velocity divergence.

By considering (x, y, z, t) as coordinates of a fourdimensional Euclidean space, E4, and using Gauss’ divergence theorem, it follows that Eq. (1) is equivalent to

the following integral conservation law:

I

S(V )

Im · dS = 0, m = 1, 2, 3, 4, 5, (2)

where S(V ) denotes the surface around a volume V in

E4 and Im = (Fm, Gm, Hm, Um).

NASA/TM—2003-212314 2




 





















































































































  







 



 













































1

2

3

4

5

6

7

8

A

B

C

Figure 1: A hexahedral cell with its six neighboring average cell centers (ACs). For clarity, only three of them,

A, B, C are shown. Together, they form a polyhedron of

24 faces or the base of the CV .

2.2 Unstructured Hexahedral Grid

Fig. 1 illustrates a hexahedron cell 1234 − 5678 in the

unstructured grid. Each of the six quadrilateral surfaces

(e.g. 1265, 2376, etc.) is associated with a neighboring hexahedral cell. Note that the four nodes of such a

quadrilateral may not lie on the same plane. Assume A,

B, C,..., are the average centers (ACs) of the respective

neighboring cells. The coordinates of the average center

of a hexahedron are the simple average of the coordinates

of its eight vertices. The average center may not necessarily coincide with the geometrical (gravity) center. The

eight vertices of each hexahedron cell and the six average

centers (A, B, C..., only three are shown in Fig. 1) of the

neighboring cells form a polyhedron of 24 faces, which,

when incorporating the time t direction, is the finite volume V in E4 in Eq. (2) or the control volume CV .

2.3 Compact Updating

For any explicit time-marching scheme, flow data at the

neighboring nodes at the previous time step are required

to update the flow data at the current node to the present

time level. In a finite difference method, a difference

scheme is utilized to do the updating. In most finite volume schemes, the computational cell (e.g., 1234-5678 in

Fig. 1) is the control volume (CV ). First, the fluxes on

the (hyper) surfaces S(V ) (e.g., 1265 in Fig. 1) need

to be updated. Conventionally, flow data is extrapolated from several neighboring nodes to the center of

the hyper-surface. In the CE/SE algorithm, a compact

updating is employed [17, 18]. The compact updating

achieves high resolution by using a cell stencil consisting only the immediate neighboring cells. For example, the CV is expanded from the hyper-hexahedral cell

1234 − 5678 to a hyper-polyhedron of 24 faces (Fig. 1),

with all the six neighboring cell centers (A, B, C...) included as vertices. When applying Eq. (2) to this hyperpolyhedron in E4, the flux associated with the face 1265

is now replaced by the fluxes associated with the four

triangles ∆A12, ∆A26, ∆A65 and ∆A51 for the new

CV . The advantage of this procedure is that flow data

at the cell center A is known and since A lies in each of

these triangles, no extrapolation through the interior of

the CV is needed when evaluating the fluxes associated

with these triangles. The updating is thus conservative.

In the CE/SE scheme, in order to achieve higher order

accuracy, flow data are extrapolated by Taylor expansion from A to the triangle centers along the triangle surface plane. Consequently, not only the conservative flow

variables U but their spatial derivatives also are considered as the unknowns (totally 20 scalar unknowns) Note

that Ut is obtained by evaluating Eq. (1). The hyper24-face-polyhedron is the control volume where Eq. (2)

is applied for conservative updating. Any of these 24

surfaces (segments of solution elements or SEs) is associated with one of the six neighboring average centers (ACs,e.g. A, B, C, ...), where the solution U and

its gradients are already given at the previous time step

n. Then both the inviscid and viscous fluxes on these

24 faces (hyperplane surfaces in E4 space) can be calculated by evaluating the flux functions F, G and H at the

face geometrical centers through Taylor expansion. Consequently, the unknown Un+1 at the new time level n+1

is evaluated at the geometrical center of the polyhedron

from the divergence theorem (2) in E4 space. Note that

the ‘geometrical center’ or centroid of the polyhedron is

in general different from its average center (AC). After

evaluating the new gradients (see next subsection) Ux,

Uy, Uz, it takes only one more step of Taylor expansion

to extrapolate Un+1 to the average center (AC).

The present compact updating avoids the uncertainty

of dimensional-splitting and extrapolation, and hence

yields better accuracy. In addition, no Riemann solver is

needed at these hyper-surfaces. Incorporation of the gradients of U further enhances the accuracy of the scheme

and makes it possible to use a compact cell stencil.

After all the Un+1 at the polyhedron centroids are updated, the next step is to calculate their corresponding

spatial gradients.

2.4 Evaluation of Spatial Gradients

As in the standard CE/SE procedure [17], unknown spatial derivatives/gradients Ux, Uy, and Uz are evaluated

using the weighted average technique or an extended van

Albada limiter [19]. The six neighboring polyhedral centers A, B, C, D, E, F around a hexahedron cell (Fig. 2)

form an octahedron. The vertices of each of the triangular face of the octahedron and the polyhedral center O

NASA/TM—2003-212314 3
A

B

C

D

E

F

O

Figure 2: An octahedron formed by the six neighboring

polyhedron centroids A, B, C, D, E, F, with the current

polyhedron cell centroid O in the middle. With O always

as a vertex, eight tetrahedra are formed and hence eight

sets of gradients exist.

associated with the current hexahedral cell form a tetrahedron (e.g. ABEO in Fig. 2). The Un+1 are already

calculated at all the vertices. For each tetrahedron, a set

of Un+1

x , Un+1y , and Un+1z can be directly calculated by

solving a linear equation system. Totally, there are eight

different tetrahedra, and so also eight sets of the spatial

derivative data. The final values of the spatial derivatives

are obtained by applying the extended van Albada limiter

(weighted average) to these eight sets of data. Let Ux

(i)

,

Uy

(i)

, Uz

(i)

, i = 1, 2, ..., 8 be respectively the gradients

from the ith set data. Let

τi = [(Ux

(i)

)

2 + (Uy

(i)

)

2 + (Uz

(i)

)

2]α,

where α is any real index number. Let ti = τ1

−1; if

τi = 0, ti = 0 and that

t = t1 + t2 + t3 + ... + t8,

then, with the van Albada limiter,

Ux = ΣUx

(i)

ti

t , Uy = ΣUy

(i)

ti

t , Uz = ΣUz

(i)

ti

t .

As is well-known, the van Albada limiter is less diffusive. In practice, it was found that if the power index number α is chosen to be slightly negative, say,

α = −0.2, the numerical dissipation resulting from averaging could be further reduced.

Once the new gradients are evaluated, Un+1 at the AC

of the hexahedron is obtained by Taylor expansion from

the centroid of the polyhedron. One marching step is thus

completed.

mesh details at the center

Figure 3: 2-D quadrilateral mesh in a cross section in a

3-D hexahedral grid

3 Hexahedral Grid Generation and Parallel

Computation

In order to substantially reduce the required number

of cells and to increase numerical stability in the computation, a hexahedral grid is adopted for the 3-D CE/SE

Navier-Stokes (N-S) solver. For the simple geometry of

a circular jet, with the x axis assumed in the stream direction, the hexahedral grid can be generated as follows:

1. a 2-D quadrilateral unstructured mesh is generated

in a circular domain on the y − z plane (Fig. 3);

2. the 2-D mesh is translated step by step in the x direction to form the 3-D hexahedral mesh. At locations that are occupied by the jet nozzle body, no

cells are generated.

In the present 3-D screech noise computation, employing a hexahedral grid instead of a tetrahedral one helps to

NASA/TM—2003-212314 4
Figure 4: Schematic diagram of parallel computation.

reduce the number of cells from tens of millions to a few

millions. Still, the number of computational cells is very

large. Due to the large number of computational cells

(about 3.67 millions in this case), parallel computation

with multiple processors becomes necessary in view of

computation turn-around time and memory size.

The parallelization procedure is similar to the one described in [15] and is sketched in Fig. 4:

1. the unstructured hexahedral grid is generated for the

entire computational domain;

2. the domain is decomposed into subdomains according to the assigned number of processes (usually

one-to-one with CPU’s, here 60 are used), using the

METIS code. METIS is an efficient mesh partitioning code and is freely available from the University

of Minnesota [20]—for example, Fig. 5 illustrates

a typical partition of the computational domain for

the current circular jet noise problem;

3. the N-S flow solver is modified to use MPI library

calls and applied to each subdomain to conduct

computations, with neighboring domains exchanging pertinent results.

MPI, or message passing interface, is an interprocessor

communication protocol standard. The software library

is prepared by the Argonne National Laboratory [21].

4 The 3-D Jet Screech Noise Problem

The circular jet in 3-D space is sketched in Fig. 6. The

flow at the nozzle exit is choked, i.e., the nozzle exit

jet

nozzle

Figure 5: A typical partition of the computational domain by METIS, with different shading or color indicating the subdomains.

Mach number, Me, is unity, and the ambient air is stationary. The case of jet Mach number Mj = 1.42 is

considered. These conditions correspond to the experimental conditions of Panda [5, 6]. In these experiments,

it was shown that for Mj = 1.42, the jet noise field exhibits truly 3-D phenomena, e.g., the flapping ‘B’ mode.

Hence, a 3-D Navier-Stokes solver described above is required.

In the investigation, our attention is focussed on the

near field of the nozzle since this is the source region

of the noise. The inner diameter, D, of the jet nozzle

is chosen as the length scale. The density, ρ0, speed of

sound, a0, and temperature T0 in the ambient flow are

taken as scales for the dependent flow variables.

In order to clearly display the upstream propagating

screech waves, the computational domain was extended

2D upstream of the nozzle exit. The full computational

domain is a circular cylinder of 10D axial length and

7.5D radius. At the nozzle exit, the inflow plane is recessed by two cells so as not to numerically restrict or

influence the feed-back loop. A straight nozzle lip of

0.25D in thickness is adopted in the computations. Note

that in the experimental setup of Panda [5, 6], the actual

nozzle diameter D = 25.4 mm and the nozzle exit lip is

beveled.

The unstructured hexahedral grid currently used is

generated as described in §2. The hexahedral cells are

non-uniform since good grid resolution is needed in the

jet shear layer and in the screech feedback loop paths.

Cell numbers in the x and radial directions are typically

NASA/TM—2003-212314 5
L

2R

S

inner

diam. D

Figure 6: Geometry of the computational domain: the

inner diameter of the nozzle is D; the nozzle extends

into the computational domain by S = 2D, the domain

length L = 8D and the radius is R = 5.5D (sponge

zones are excluded)

400 and 160 respectively. In the azimuthal direction,

there are 64 cells, corresponding to totally about 3.67

million hexahedral cells. The last 6 cells in the streamwise and radial directions have exponentially growing

size, respectively, and serve as sponge zones to substantially eliminate any small remaining numerical reflection from the outflow and circumferential non-reflecting

boundaries.

4.1 Initial Conditions

Initially, the flow of the entire domain is set at the ambient flow conditions (using nondimensional variables),

i.e.,

ρ0 = 1, p0 = 1

γ ,

u0 = 0, v0 = 0, w0 = 0.

4.2 Boundary Conditions

At the inlet boundary, the conservative flow variables and

their spatial derivatives are specified to be the same as the

ambient flow, except at the nozzle exit, where an elevated

pressure is imposed, i.e., the jet is under-expanded, as in

the physical experiments. By using the ideal gas isentropic relations, it follows that the nondimensional flow

variables at the nozzle exit, with Me = 1, are given by

ρe = γ(γ + 1)pe

2Tr

,

pe = 1

γ

"

2 + (γ − 1)M2

j

γ + 1 # γ

γ−1

,

ue =

 2Tr

γ + 11/2

, ve = 0, we = 0,

where Tr is the reservoir (plenum) temperature. We will

also follow the experimental cold-flow condition where

the reservoir temperature equals the ambient one, i.e.,

Tr = 1.

At the circumferential and outflow boundaries, the

Type II and Type I CE/SE non-reflecting boundary conditions as described in the next subsection are imposed,

respectively. The no-slip boundary condition is applied

on all the nozzle walls.

4.3 Non-Reflecting Boundary Conditions

In the CE/SE scheme, non-reflecting boundary conditions (NRBC) can be easily constructed based on planewave propagation theory for hyperbolic conservation

laws [22]. There are various implementations of the

non-reflecting boundary condition (NRBC) and in general they have proven to be well suited for aeroacoustic

problems [14, 23]. The following 3-D NRBCs are employed in this paper.

For a ghost grid node (j, n) lying at the outer radius of

the domain the non-reflective boundary condition (Type

II) requires that

(Ux)

n

j = (Uy)

n

j = (Uz)

n

j = 0,

while Un

j is kept fixed at the initial steady boundary

value. At the downstream boundary, where there are

substantial gradients in the radial direction, the nonreflective boundary condition (Type I) requires that

(Ux)

n

j = 0,

while Un

j , (Uy)n

j and (Uz)nj are now defined by simple

extrapolation from the nearest interior node j0

, i.e.,

Un

j = Unj0 ,

(Uy)

n

j = (Uy)

n

j0 , (Uz)

n

j = (Uz)

n

j0 .

As will be observed later, these NRBCs, when combined

with the above buffer/sponge zones, are robust enough

to allow a near field computation without disturbing or

distorting the flow and acoustic fields.

5 Numerical Results

In this section, 3-D numerical results for the underexpanded circular jet described above are presented and

compared to experimental results [5, 6]. Computations

are conducted for the jet Mach number Mj = 1.42. At

this moderate jet Mach number, the dominant unsteady

motion in the experiments, see [5], is truly three dimensional. With an appropriate time step size, a large number (510,000) of time steps are performed in order to

achieve sufficient accuracy in the Fourier analysis of time

series data. Note that no harmonic forcing is imposed

in the numerical simulation. The initial impact of the

boundary condition at the nozzle exit stimulates the jet

shear layer and triggers the feedback loop that generates

the (then) self-sustained screech waves.

NASA/TM—2003-212314 6
Figure 7: Time-averaged experimental (top) and instantaneous numerical (bottom) Schlieren pictures, showing

shock cell structure; Mj = 1.42.

5.1 Shock-Cell Structure

Experimental results for jets are often documented in

terms of Schlieren pictures. It is straightforward to construct Schlieren plots (density-gradient modulus) from

the numerical results. For the case of Mj = 1.42,

Fig. 7 shows the experimental [5, Fig. 4(a)] and numerical Schlieren pictures. Good agreement in shockcell structure is shown. For example, the shock cell

width (spacing) is about 1.28D in the streamwise direction. Note that the experimental Schlieren plot is a timeaveraged results, while the numerical one is a snapshot at

a relatively early stage (time step 110,000). In the experimental Schlieren plot, it is observed that the first shock

cell appears to be sharp and clear since the shear-layer instability wave is too weak at this location to significantly

affect the shock cell. However, once the instability wave

has gained a sufficient amplitude through its streamwise

growth, it interacts strongly with the shock cells. This

is evident from the deformations from the second shock

cell and onwards and the additional blurring due to timeaveraging in top panel of Fig. 7.

5.2 Near-Field Radiating Screech Waves

Figs. 8- 12, which represent a series of instantaneous

snapshots of pressure iso-surfaces in the flow field, illustrate the generation and propagation of screech waves.

Since no forcing is applied in the numerical simulation,

these waves are a clear indication of a self-sustained oscillation. For Mj = 1.42, the jet screech is in the flapping B mode, which is a truly three dimensional one. The

screech waves not only propagate in the upstream direction but also swirl and flap around the jet core.

Figure 8: Pressure iso-surfaces at time step 120,000; a

x − z cross section slice is also shown.

After the computation has run for about 110,000 time

steps, truly three dimensional asymmetric wave patterns

begin to appear. Fig. 8 and Fig. 9 demonstrate the pressure iso-surfaces at time steps 120,000 and 130,000 respectively and from different angle. The iso-surfaces

clearly take on nonaxisymmetric forms. However, at

time steps of 200,000 and 320,000, the pressure isosurfaces again turn out to look more axisymmetric, as

demonstrated in Fig. 10 and Fig. 11 respectively. At

time step 510,000, when the computation was stopped, a

highly asymmetric pattern can be observed, see Fig. 12.

5.3 Screech Frequency and Sound Pressure

Level (SPL)

The numerical time history is recorded for the location

(2.0,0.6,0) at the nozzle lip in the flow field and later

post-processed to obtain spectral information using Fast

Fourier Transform (FFT) techniques. The recording begins after an initial time period has elapsed (180, 000 actual time steps) ensuring that any start-up transients have

left the computational domain. Figure 13 displays the

SPL for the case of Mj = 1.42. The SPL plot shows

the spikes at the fundamental ‘B’ mode frequency and

its harmonics. The fundamental frequency spike correspond to screech frequency of 4500 Hz (SPL=122 dB),

which agrees quite well with the experimental value of

4350 Hz.

5.4 Performance of the Parallel Computation

The computation was mainly carried out on a Linux Pentium III PC cluster (NASA Glenn Research Center CW-7

cluster) using 60 processes, running on 20 CPUs. Normally, one would want a one-to-one relationship between

NASA/TM—2003-212314 7
Figure 9: Pressure iso-surfaces at time step 130,000.

isobars

200k steps

Figure 10: Pressure iso-surfaces at time step 200,000.

iso−surface levels: .7132, .71355, .7137

Figure 11: Pressure iso-surfaces at time step 320,000; a

x − z cross section slice is also shown.

the number of processes and processors. This less than

optimal situation lead to about 15 % increase in the wallclock time and was due to the computation having been

started on a different Linux PC cluster using 60 CPUs.

In the optimal situation, it takes about 8 seconds wallclock time to march one step on the CW-7 cluster. The

parallel computation was also tested on the SGI (Silicon Graphics ) Origin 3800 workstation cluster with 64

and 128 processors at NASA Ames Research Center.

From 64 to 128 processors, the clock time reduces linearly with the increasing number of processors, but the

reduction will begin to level off if more processors are

used. As a result of the explicit time-marching in the

scheme, when running with 64 processors, the number

of MFLOPS (megaflops) per processor remains between

170 and 183, exceeding 20 % (160) of the theoretical

peak MFLOPS. This performance is considered excellent by code-performance specialists at NASA Ames.

6 Concluding Remarks

In this paper, a 3-D CE/SE N-S solver using an unstructured hexahedral grid is briefly described and tested

in a 3-D circular jet screech problem. The use of a hexahedral grid rather than a tetrahedral one enhances the numerical stability of the scheme and significantly reduces

both memory size and CPU time, making the CE/SE

scheme a viable tool for near-field CAA simulation.

For the test case of Mj = 1.42, the jet screech frequency of the dominant nonaxisymmetric ‘B’ mode and

the shock cell structure agree well with the experimental data [5, 6]. Perhaps due to the relatively coarse grid

used, the computed SPL is somewhat lower than the experimental one. Further tests with refined grids are being

carried out and will be reported in the future.

NASA/TM—2003-212314 8
Figure 12: Pressure iso-surfaces at time step 510,000

viewed from two different angles; x − z cross section

slices are also shown. SPL (dB)

frequency (Hz)

CE/SE simulation: M_j = 1.49 -- simplified nozzle geometry, 3.7M cells

(J. Panda experiment: B-mode 4350 Hz 161 dB)

V

| 4500 Hz

2nd harmonic

3rd harmonic

 0. 5000. 10000. 15000. 20000. 

 60. 

 100. 

 140. 

 180. 

Figure 13: SPL at nozzle wall.

References

[1] Tam, C. K .W., “Supersonic Jet Noise,” Ann. Rev.

Fluid Mech. vol. 27, pp. 17-43 (1995).

[2] Seiner, J. M., “Advances in High Speed Jet Aeroacoustics,” AIAA Paper 84-2275 (1984).

[3] Tam, C. K. W., “Jet Noise Generated by Large

Scale Coherent Motion,” NASA RP-1258, pp. 311-

390 (1991).

[4] Raman G., “Advances in Understanding Supersonic Jet Screech: Review and Perspective,” Prog.

Aerosp. Sci. vol. 34, pp. 45-106 (1998).

[5] Panda, J., “Shock Oscillation in Underexpanded

Screeching Jets,” J. Fluid Mech., vol. 363, pp. 173-

198 (1998).

[6] Panda, J., “An Experimental Investigation of

Screech Noise Generation,” J. Fluid Mech., vol.

378, pp. 71-96 (1999).

[7] Panda, J., Raman, G. and Zaman, K. B. M. Q., “Underexpanded Screeching Jets from Circular, Rectangular and Elliptic Nozzles”, AIAA paper 97-

1623 (1997).

[8] Ponton, M. K., Seiner, J. M. and Brown, M. C., “

Near Field Pressure Fluctuations in the Exit Plane

of a Choked Axisymmetric Nozzle”, NASA TM

113137 (1997).

[9] Shen, H. and Tam, C. K. W., “Numerical Simulation of the Generation of Axisymmetric Mode Jet

Screech Tones,” AIAA Paper 98-0283 (1998).

[10] Shen, H. and Tam, C. K. W., “Three-Dimensional

Numerical Simulation of the Jet Screech Phenomenon” AIAA Paper 2001-0820 (2001).

[11] Imamoglu, B. and Balakumar, P. “Computation of

Shock Induced Noise in Imperfectly Expanded Supersonic jets”, AIAA Paper 2002-2527 (2002).

[12] Loh, C. Y., Hultgren, L. S., and Jorgenson, P. C.

E., “Near Field Screech Noise Computation for an

Underexpanded jet by the CE/SE Method,” AIAA

Paper 2001-2252 (2001).

[13] Loh, C. Y. and Hultgren, L. S., “Computing Jet

Screech - A Complex Aeroacoustic Feedback System”, presented at the 2nd Int’l Conference on CFD

(ICCFD2), also NASA/TM-2002-211807 (2002).

[14] Jorgenson, P. C. E. and Loh, C. Y., “Computing Axisymmetric Jet Screech Tones Using Unstructured

Grids”, AIAA Paper 2002-3889 (2002).

NASA/TM—2003-212314 9
[15] Loh, C. Y., Himansu, A. Wang, X.-Y. and Jorgenson, P. C. E., “Computation of an Underexpanded

3-D Rectangular jet by the CE/SE Method”, AIAA

Paper 2001-0986 (2001).

[16] Zhang, Z.-C., Yu, S. T. and Chang, S.-C., “A SpaceTime Conservation Element and Solution Element

method for Solving the Two and Three Dimensional Unsteady Euler Equations Using Quadrilateral and Hexahedral Meshes”, J. Comp. Phys. vol.

175, pp. 168-199 (2002).

[17] Chang, S.-C., Wang, X.-Y. and Chow, C.-Y., “The

Space-Time Conservation Element and Solution

Element Method—A New High Resolution and

Genuinely Multidimensional Paradigm for Solving

Conservation Laws,” J. Comp. Phys. vol. 159, pp.

89-136 (1999).

[18] Wang, X.-Y. and Chang S.-C., “ A 2-D Nonsplitting Unstructured Triangular Mesh Euler

Solver Based on the Space-Time Conservation Element and Solution Element Method” C.F.D. J. vol.

8, pp. 309-325 (1999).

[19] van Albada, G. D., van Leer, B. and Roberts, W. W.,

“A Comparative Study of Computational Methods

in Cosmic Gas Dynamics,” Astronomy and Astrophys., vol. 108, pp. 76-84 (1982).

[20] Karypis, G. and Kumar, V. “Multilevel k-way Partitioning Scheme for Irregular Graphs”, Univ. of

Minnesota Dept. of Comp. Sc./Army HPC Research Center Tech. Report 95-064 (1995).

[21] website: http://www-unix.mcs.anl.gov/mpi

[22] Loh, C. Y., “On a Nonreflecting Boundary Condition for Hyperbolic Conservation Laws,” AIAA Paper 2003-3975 (2003).

[23] Loh, C. Y., Hultgren, L. S. and Chang S.-C.,

“Computing Waves in Compressible Flow Using

the Space-Time Conservation Element Solution Element Method,” AIAA J., Vol. 39, pp. 794-801

(2001).

NASA/TM—2003-212314 10
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Technical Memorandum

Unclassified

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field

Cleveland, Ohio 44135–3191

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546–0001

Available electronically at http://gltrs.grc.nasa.gov

June 2003

NASA TM—2003-212314

AIAA–2003–3207

E–13903

WBS–22–708–90–25

16

A 3-D CE/SE Navier-Stokes Solver With Unstructured Hexahedral Grid

for Computation of Near Field Jet Screech Noise

Ching Y. Loh, Ananda Himansu, and Lennart S. Hultgren

CE/SE method; Unstructured hexahedral grid; Jet screech noise

Unclassified -Unlimited

Subject Categories: 02, 71, and 61 Distribution: Nonstandard

Prepared for the Nineth Aeroacoustics Conference and Exhibit cosponsored by the American Institute of Aeronautics

and Astronautics and the Confederation of European Aerospace Societies, Hilton Head, South Carolina, May 12–14,

2003. Ching Y. Loh and Ananda Himansu, Taitech, Inc., Beaver Creek, Ohio 45430; Lennart S. Hultgren, NASA Glenn

Research Center. Responsible person, Ching Y. Loh, organization code 5880, 216–433–3981.

A 3-D space-time CE/SE Navier-Stokes solver using an unstructured hexahedral grid is described and applied to a

circular jet screech noise computation. The present numerical results for an underexpanded jet, corresponding to a fully

expanded Mach number of 1.42, capture the dominant and nonaxisymmetric 'B' screech mode and are generally in

good agreement with existing experiments.