NASA Contractor Report 175098

AIAA-86-1520

An LU Implicit Scheme for High

Speed Inlet Analysis

r

 (NftSa-CR-175098) iB: LU IMPLICIT! SCHEHE ?OR N86-23563

HIGH SPEED INLET ANALYSIS Final Beport

{Sverdrup Technology, Inc.) 17 p

HC AQ2/HF A01 CS.CL ;0;U Unclas

G3/02 _ 05983 __

Seokkwan Yoon

Sverdrup Technology, Inc.

Lewis Research Center

Cleveland, Ohio

and

Antony Jameson

Princeton University

Princeton, New Jersey

April 1986

Prepared for

Lewis Research Center

Under Contract NAS 3-24105

NASA

National Aeronautics and

Space. Administration
AN LU IMPLICIT SCHEME FOR HIGH SPEED INLET ANALYSIS

Seokkwan Yoon*

Sverdrup Technology, Inc.

Lewis Research Center

. Cleveland, Ohio

and

Antony Jameson*

Princeton University

; Princeton, New Jersey 08544

ABSTRACT

A numerical method 1s developed to analyze the 1nv1sc1d flowfield of a

high speed Inlet by the solution of the Euler equations. The LU implicit

scheme 1n conjunction with adaptive dissipation proves to be an efficient and

robust nonosdllatory shock capturing technique for high Mach number flows as

well as for transonic flows.

INTRODUCTION

Recent interest in the aerospace plane and other hypersonic vehicles

revitalized the research on high speed propulsion systems as well as hypersonic

aerodynamics. In the design of supersonic and hypersonic propulsion systems

the analysis of high speed flow past an Inlet plays a critical role. While

there are half a dozen propulsion study concepts for high speed flight, many

candidate concepts share a common Idea of combination of turboramjet engines

for sub and supersonic flights and scramjet (supersonic combustion ramjet)

engines for hypersonic flight. Resulting speeds of flows through the engines

range widely from subsonic to hypersonic regimes.

The analysis of hypersonic flows would require the full Navler-Stokes

equations with slip effects and chemical reaction. It 1s also Important to

*

Senior research engineer; member AIAA.

Professor of Mechanical and Aerospace Engineering; member AIAA.
understand the complex structure of shock waves. The turboramjet Inlet

flowfleld Includes the Incoming supersonic flow deflected by oblique shock

waves and the subsonic dlffuser flow after the terminal normal shock wave

while the scramjet Inlet flow 1s characterized by strong oblique shock waves.

The Euler equations which represent hyperbolic conservation law can be a useful

testbed for developing and evaluating a shock capturing numerical algorithm.

It has been a difficult task for computational aerodynamldsts to capture

nonosdllatory shock waves as a converged solution. Unbounded growth of

spurious oscillations often resulted 1n numerical Instability. It Is well

known that upwind difference schemes can eliminate oscillations In the

neighborhood of shock waves at the expense of a substantial Increase of

computational work. In parallel with the developments 1n upwind schemes It

has been found that steady aerodynamic flows containing moderately strong shock

waves can be quite well predicted by a central difference scheme augmented by

a carefully controlled blend of first and third order d1ss1pat1ve terms.

In this paper the performance of adaptive dissipation 1s demonstrated for

strong oblique shock waves 1n high Mach number flows on a near-uniform mesh.

Although a space marching method has been useful, 1t 1s not well-posed

when there 1s upstream Influence through subsonic portions of the flowfleld

such as a boundary layer. It also cannot handle the terminal shock and the

subsonic dlffuser flowfleld as well as the flow with streamwlse separation.

Early time-Integration codes for calculating the supersonic flow through an

2 3 Inlet used popular HacCormack schemes. ' A disadvantage of these schemes

for steady state calculation 1s that the computed steady state depends on the

time step. Another drawback of HacCormack1 s Implicit scheme ,1s the difficulty

in treating boundary conditions.

During the last decade, the Navier-Stokes equations have been the subject

of exploratory Investigations aimed at establishing the feasibility of their

2
solution, but the methods so far developed have been too expensive to permit

their use 1n a routine production mode. However, recent developments of modern

Implicit schemes ~ 1n conjunction with multlgrld methods are encouraging. The

authors developed an optimal dissipation model for the alternating direction

Implicit (ADI) scheme and proved that the Improved ADI scheme 1s Ideal for

multlgrld 1n two dimensions. Unfortunately, the ADI scheme 1n delta form to

ensure the time step Independent solution has stability and convergence

problems 1n three dimensions. Two new Implicit schemes which are

unconditionally stable in any number of space dimensions were successfully

developed by the authors recently. They are lower-upper (LU) Implicit

scheme and LU-SSOR (symmetric successive over-relaxation) scheme. The LU

implicit scheme has been tried on an H-mesh 1n this work for high speed flow

calculations:

GOVERNING EQUATIONS

The Euler equations are obtained from the Navler-Stokes equations by

neglecting viscous terms. Let p, u, v, E, H, and p be the density,

Cartesian velocity components, total energy, total enthalpy, and pressure, and

let x and y be Cartesian coordinates. Then for a two-dimensional flow

these equations can be written as

3W 3F 36 n ...

at * ax * ay

 = ° (1)

where W 1s the vector of dependent variables, and F and G are convectlve

flux vectors

W = (p,pu,pv,pE)T

2 T F = (pu.pu + p,pvu,u(pE * p)) (2)

G = (pv.puv.pv2 + p,v(PE + p))T

The pressure is obtained from the equation of state
,..,.. . p = p(Y - 1)|E.- 1 (u2 + y2)J . ,. . _. (3)

These equations are to be solved for a steady state aw/at = 0 where t

denotes time. •

1 . SEMI-DISCRETE FINITE VOLUME METHOD

A convenient way to assure a steady state solution Independent of the time

step is to separate the space and time discretization procedures. In

semi-discrete finite .volume method one begins by applying a semi-discretization

1n which only the spatial derivatives are approximated. The use of a finite

volume method for space discretization allows one to handle arbitrary-1

geometries and helps one. to avoid problems with metric singularities tha.t;are

usually associated with finite difference methods. The scheme reduces to a

central difference scheme on a Cartesian grid, and is second order accurate in

space provided that the mesh 1s smooth enough. It also has the property that

uniform flow is an exact solution of the difference equations.

;

 ;- NONLINEAR ADAPTIVE DISSIPATION

In typical calculation of flow with discontinuities by a central

difference scheme, wiggles appear in the neighborhood of shock waves where

pressure gradient 1s severe. In order.to suppress the tendency for spurious

odd and even point oscillations, and to prevent unsightly overshoots near

shock waves, the scheme 1s augmented by artificial disslpatlve terms. The

d1ss1pat1ve term, which 1s constructed so that 1t 1s of third order 1n smooth

' 'jf' '-.'' :

regions of the flow, 1s explicitly added to the residual. For the density

equation, for example, the dissipation has the form

dUl/2,j d1-l/2,j * d1,j+l/2 " d1J-l/2

where

Hl/2J(pUlJ -

 P1J> - eU1/2,j(pU2J ' 3pHl,J * 3p1,j

(4)
Let S be the cell area which 1s equivalent to the Inverse of the determinant

of transformation Jacoblan. Both coefficients Include a normalizing factor

S (2)

1 ,/2 . , t ^proportional to the length of the cell side, and e^/2 < 1s

also made proportional to the normalized second difference of the pressure

(5)

1n the adjacent cells. The third order terms provide background damping of

high frequency modes. The first order terms are needed to control oscillations

1n the neighborhood of shock waves, and are turned on by sensing strong

pressure gradients in the flow. The dlsslpatlve terms for the other equations

are constructed from similar formulas with the exception of the energy equation

where.'the differences are of pH rather than pE. The purpose of this 1s to

allow.a:steady state solution for which H remains constant. Increasing the .

amount of artificial viscosity Improves the rate of convergence although too

much dissipation can hurt 1t. However, 1t 1s desirable to make the amount be

as small as possible 1n order not to degrade the accuracy of solution. Typical

amount of the third order terms 1s almost negligible when compared to the

physical viscosity.

LU IMPLICJT SCHEME

Let the Jacoblan matrices be

- .

 A = aw ' B = aw

and let the correction be

,n 6W = W" ' - W1

here n denotes the time level.

The linearized Implicit scheme for a system of nonlinear hyperbolic

equations such as the Euler equations can be formulated as

(l +• 0 At (D A + D B)}«W * At R = 0 (7)

J

"

; 5 '
where R 1s the residual "

, . -,. ..... . , R = DXF(WV+ DyG(Wn) --..:-. ; ..-;., , ,

Here D and D are central difference operators that approximate a/ax

x y

and a/ay'.' • ' """ '

If B = 1/2 the scheme, remains second order accurate 1n time, for other

values of 6, the time accuracy drops to first order. The unfactored Implicit

scheme (Eq. (7)) produces a large block banded matrix which 1s very" costly "to"

Invert and requires hugh storage. One can solve the system Indirectly using a

relaxation algorithm. Then It is desirable that the matrix should be

diagonally dominant ""'to meet a convergence criterion of a relaxation method.

This cari be achieved by flux splitting at the expense of a substantial increase

in the computational work. Moreover, it seems that second order flux splitting

methods 1ri conjunction with relaxation algorithms are either conditionally -

stable or slow. • <

The operation count can be reduced by factorlzlng the operator of Eq. (7)

approximately 1n various ways. The first way 1s known as the ADI scheme. ':(

(I + B At 0 A)(I + B At D B)SW 4- At R = 0 (8)

x y

Although the Introduction of optimal artificial dissipation makes the

scheme be very desirable 1n two dimensions, the scheme 1n delta form 1s only

conditionally stable 1n three dimensions. The ADI scheme Introduces error

'(... 3 " :

terms of order (At) 1n three dimensions which reduce the convergence rate.

If one concerns about memory requirement, each factor can be split into two

Q

subfactors. If B = 1. the scheme becomes

(I + At D~A*)(I t At D*A~)(I -i- At 0~B*)(I + At' D*B''~)«W•+ At R =

;0 (9)

where D~ and D~ are backward difference operators and D and D

x y • o-•-. - • • . ..x •• .. y .

are forward difference operators. Each factor can be constructed using the

g

diagonally dominant ADI factorization. This scheme has six factors in three

6
dimensions and Introduces error terms of order (At) which reduce the

convergence rate further.-;

While the ADI scheme has been valuable 1n two dimensions, Its Inherent

limitations 1n three dimensions suggest an alternative approach. An

unconditionally stable Implicit scheme which has error terms at most of order

2

(At) 1n any number of space dimensions can be derived by the LU

factorization.

(l t (3 At" (0~A* + D~B*)}(l * 13 At (DV * oV))«W * At R = 0 (10)

x y x y J

Here, A , A~, B , and B~ are constructed so that the eigenvalues of "*"

matrices are nonnegative and those of "-" matrices are nonposltlve.

•A* = \ (A * rAI), A" = \ (A - rAI)

* • j • ,• ' ...

B* = \ (B * rBI), . B- = \ (B - rBI)

(11)

where •, -

rA > max(|XA | ), r B >max(|x B | ) (12)

Here, X. and X_ represent eigenvalues of Jacob.lan matrices. Equation (10)

A D -

can be Inverted 1n .two..steps i. The LU Implicit scheme needs the Inversion of

sparse triangular matrices which can be done efficiently without using large

storage. This scheme has only two factors 1n three dimensions. Other forms

of factorization 1n conjunction with flux splitting can be found 1n Ref. 10.

For example, ,

(l t 0 At (D~A+ + DyB)}'(l * B At DV 6W + At R = 0} (13)

This scheme requires a block-tr.1d1agonal Inversion 1n one direction. If one

wants to Include thin-layer viscous terms 1n the Implicit operator, this scheme

may be useful. However, 1t does not seem to be necessary to Insert viscous

terms Into the operator when only the steady-state solution 1s desired.
--.- - ----- - - RESULTS- .-.----, . -- . : .

In order to test the performance of the LU Implicit scheme for high speed

flow calculations, a two-dimensional model problem was selected: Figure 1

shows a typical hypersonic Inlet and F1g. 2 shows a 54 by 32 H-mesh for a

schematic high speed Inlet. This mesh was used 1n all cases except the

terminal shock wave problem where a 104 by 32 mesh was used for better

resolution of the normal shock wave. However, all figures of convergence

history are the results on the 54 by 32 mesh for comparison. The ramp angle

is 9° and the shoulder angle 1s set to 0.5° for the terminal shock wave

problem.

At the inflow boundary all the flow quantities are specified, and they

are extrapolated from the Interior at the outflow boundary for supersonic

outf-low.- For the terminaL shock wave^problem the^pressure.was prescribed at

the outflow boundary.

Four plots including Mach number contours, Mach number and pressure along

the centerline, and the convergence history are shown 1n each set of Figs. 3

to 6. Terminal shock wave problem with freestream Mach number 2 1s shown 1n

Fig. 3. Figures 4 to 6 are for supersonic throughflows with freestream Mach

numbers 5, 10, and 20, respectively. The pressure plots show the values of

the pressure normalized by freestream static pressure so that the strengths of

shock waves can be compared. Two Indicators 1n the convergence histories are

the maximum and the average density residual 1n logarithm scale.

As the figures show the wave structures 1n high Mach number flows

Including oblique shock waves, reflected shock waves, expansion fans, and the

Interaction of shock waves with expansion fans are successfully captured.

These results clearly demonstrate the capability of the present numerical

method for high speed flows. Figures 5 and 6 show that the location of the

shock waves 1s hardly changed as the Mach number Increases from 10 to 20.

8
However, the pressure plots show that the strengths of shock waves are quite

different.

The convergence histories show that the residuals drop linearly and

continuously. These prove the efficiency of the present numerical method.

However, as the Mach number Increases the convergence rate 1s slowed down.

This problem will be Investigated 1n the future. Another difficulty

encountered 1n high Mach number flows 1s the sensitivity to the way of starting

the solution procedure. Sudden Introduction of the boundary condition to the

freestream uniform flow 1s likely to cause numerical instability 1n high Mach

number flows. Gradual increase of the time step 1s found to be effective to

fix this problem.

CONCLUSION

The LU implicit scheme combined with the nonlinear adaptive dissipation

is successfully developed as a robust and efficient shock capturing method for

high Mach number Inlet flows. It seems to be possible to Improve the

resolution and the accuracy of shock waves by using a total variation

diminishing scheme. Extension to the Navier-Stokes equations 1s desirable for

more-accurate "simulation of the flow through Integrated high speed propulsion

system.

ACKNOWLEDGMENT

The first author would like to thank Dr. Peter Sockol 1n NASA Lewis

Research Center for helpful suggestions. Thanks also go to Drs. Rod Chima,

E11 Turkel, and'Al Bishop for their help 1n graphics.

REFERENCES

1. Jameson, A., "Steady State Solution of the Euler Equations for Transonic

Flow," Transonic. Shock, and Multidimensional Flows, edited by R.E.

Meyer, Academic Press, New York, 1982, pp. 37-70.
2. MacCormack, R...W., "The Effect of Viscosity In Hype rye locltyImpact

CrateMng," AIAA Paper 69-354, Apr. 1969.

3. MacCormack, R.W., "A Numerical Method for Solving the Equations of

Compressible Viscous Flows," AIAA Paper 81-0110, Jan. 1981.

4. Jameson, A., and Baker, T.J., "Solution of the Euler Equations for

Complex Configurations," AIAA Paper 83-1929, July 1983.

5. Jameson, A., and Yoon, S., "Multigrld Solution of the Euler Equations

. using Implicit Schemes," AIAA Paper 85-0293, Jan. 1985.

6. Jameson, A., and Yoon, S.," LU Implicit Schemes with Multiple Grids for

.the Euler Equations," AIAA Paper 86-0105, Jan. 1986.

7. Beam, R., and Warming, R.F., "An Implicit Factored Scheme for the

Compressible Navler-Stokes Equations," AIAA Journal. Vol. 16, No. 4,

_.Ap_r._ 1978_, pp. 393-402.

8. Obayashl, S., Kuwahara, K., Fujli, K., and Matushima, K., "Improvements

1n Efficiency and Reliability for Navler-Stokes Computations using the

LU-AOI Factorization Algorithm," AIAA Paper 86-0338, Jan. 1986.

9. Lombard, C.K., Bardina, J., Venkatapathy, E., and Oliger, J.,

"Multidimensional Formulation of CSCM - An Upwind Flux Difference

Eigenvector Split Method for the Compressible Navier-Stokes Equations,"

AIAA Paper 83-1895, July 1983.

10. Steger, J.L... "Implicit Finite Difference Simulation of Invisdd and

Viscous Compressible Flow," Transonic. Shock, and Multidimensional Flows.

edited by R.E. Meyer, Academic Press, New York, 1982, pp. 181-217.

11. Pulllam, T.H., "Euler and Thin Layer Navler-Stokes Codes: ARC2D, ARC3D,"

.. CFD User's Workshop Notebook. Univ. of Tennessee, 1984, pp. 15-1 to 15-85.

10
ORIGINAL PAGE IS

OF POOR QUALfTY

Figure 1. -A typical hypersonic inlet.

Figure 2. - 54x32 H-mesh for a schematic high speed inlet. Ramp angle

9° and shoulder angle OP or 0.5°.
(a) Mach number contours.

.2 .4 .6 .8

(b) Mach number along the centerline.

1.0

4.5

3.5

2.5

1.5

.5

.2 .4 .6 .1

(cl Pressure along the centerline.

1.0

-2

o

3 -4

100 200 300 400 500

(d) Convergence history.

Figure 3. -Mach 2 inlet with the terminal Shockwave.
(a) Mach number contours.

<c

ZJ

o

o

o

5.2

4.8

4.4

4.0

-2

-8

-10

' 0 .2 .4 .6 .8 1.0

(b) Mach number along the centerline.

3 i—

1

 0 .2 .4 .6 • .8

- (c) Pressure along the centerline.

2 r~

100 200 300

(d) Convergence history.

Figure 4. -Mach 5 inlet.

400

1.0

500
(a) Mach number contours.

ce.

LLJ

00

o

<t

<t

3

o

o

p

10.5

9.5

8.5

7.5

6.5

8 r—

2 |—

-2

-10

.2 .4 .6 .8

(b) Mach number along the centerline.

.2 .4 .6 .8

(c) Pressure along the centerline.

1.0

_]

1.0

0 . 100 : - .200 300 400

(d) Convergence history. •

Figures. -Mach 10inlet.

500
241—

(a) Mach number contours.

20

16

12

o.

J/5

Q25.0

20.0

15.0

10.0

5.0

.2 .4 .6 .8

(b) Mach number along the centerline.

T'"T""T

.2 .4 .6 .8

(c) Pressure along the centerline.

1.0

1.0

n

O

O

-2

-4

I . I

100 200 300 400

(d) Convergence history.

Figure 6. -Mach 20 inlet.

500
1. Report NO. NASA CR_-| 75098

AIAA 86-1520

2. Government Accession No. 3. Recipient's Catalog No.

4. Title and Subtitle

An LU Implicit Scheme for High Speed Inle.t Analysis

5. Report Date

April 1986

6. Performing Organization Code

7. Authors)

Seokkwan Yoon and Antony Jameson

8. Performing Organization Report No.

None

10. Work Unit No.

9. Performing Organization Name and Address

Sverdrup Technology, Inc. and Princeton University

Lewis Research Center Princeton, New Jersey

Cleveland, Ohio 44135 08544

11. Contract or Grant No.

NAS 3-24105

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Washington, D.C. 20546

13. Type of Report and Period Covered

Contractor Report

14. Sponsoring Agency Code

505-62-21 (E-3022)

15. Supplementary Notes

Final report. Project Manager, Peter M. Sockol, Internal Fluid Mechanics

Division, NASA Lewis Research Center, Cleveland, Ohio 44135. Prepared for the

22nd Joint Propulsion Conference, cosponsored by the AIAA, ASME, SAE, and ASEE,

Huntsvllle, Alabama, June 16-18, 1986.

16. Abstract

A numerical method 1s developed to analyze the 1nv1sc1d flowfleld of a high speed

Inlet by the solution of the Euler equations. The LU Implicit scheme 1n conjunc

tlon with adaptive dissipation proves to be an efficient and robust nonosdllatory shock capturing technique for high Mach number flows as well as for transonic flows.

17. Key Words (Suggested by Author(s))

High speed flow

Computational aerodynamics

18. Distribution Statement

Unclassified - unlimited

STAR Category 02

9. Security Classif. (of this report)

Unclassified

20. Security Classif. (of this page)

Unclassified

21. No. of pages 22. Price'

*For sale by the National Technical Information Service, Springfield, Virginia 22161
National Aeronautics and

Space Administration

Lewis Research Center

Cleveland. Ohio 44135

Official Business

Penalty for Private Use $300

SECOND CLASS MAIL

ADDRESS CORRECTION REQUESTED

Postage and Fees Paid

National Aeronautics and

Space Administration

NASA-451

NASA