NASA Contractor Report 178053

leASE REPORT NO. 86-5

ICASE

APPLICATION OF A RUNGE-KUTTA SCHEME

FOR HIGH-SPEED INVISCID INTERNAL FLOWS

Anutosh Moitra

Eli Turkel

Ajay Kumar

Contract Nos. NASl-17D7D and NASl-18lD7

January 1986

NASA-CR-178053

19860011253

L'.'. ~l\_;' f..'c'3[t.!\C:-j CEnTER

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SJ\

National Aeronautics and

Space Administration

Langley Research center

Hampton.Virginia 23665

APPLICATIONOF A RUNGE-KUTTASCHEME

FOR HIGH-SPEEDINVISCIDINTERNALFLOWS

Anutosh Moitra

SASC Technologies, Inc., Hampton, Virginia 23665

Eli Turkel

Institute for Computer Applications in Science and Engineering

NASA Langley Research Center, Hampton, Virginia 23665-5225

and

Tel-Aviv University

Ajay Kumar

NASA Langley Research Center, Hampton, Virginia 23665-5225

ABSTRACT

• A multi-stage Runge-Kutta method is analyzed for solving the twodimensional Euler equations for external and internal flow problems.

Subsonic, supersonic and, highly supersonic flows are studied. Various

techniques for accelerating the convergence to a steady state are described

and analyzed. Effects of the grid aspect ratio on the rate of convergence are

evaluated. An enthalpy damping technique applicable to supersonic flows is

described in detail. Numerical results for supersonic flows containing both

oblique and normal shocks are presented confirming the efficiency of the

method.

Research was supported by the National Aeronautics and Space

Administration under NASA Contract Nos. NASI-17070 and NASI-18107 while the

second author was in residence at ICASE, NASA Langley Research Center,

Hampton, VA 23665-5225.

I° INTRODUCTION

Over the past few years, a number of numerical schemes have been devised

for solving the compressible Euler equations. An explicit finlte-volume scheme

with multi-stage Runge-Kutta integration in time, as proposed by Jameson,

Schmidt, and Turkel [i], has shown promise of being an efficient and accurate

Euler solver. This method has been found to be of great versatility in the

simulation of transonic flow, particularly when the method has been augmented

by such convergence accelerating techniques as local time-stepping, enthalpy

damping, and implicit residual averaging. The incorporation of these

techniques as described by Turkel [2] has produced remarkable increases in the

rate of convergence to a steady-state. However, the extent to which these

enhancements are necessary vary significantly from one problem to another.

Although the preceding method has been used successfully for transonic flow,

its use in high-speed flows involving shocks of substantial strength has not

been explored fully. Moitra [3] has presented some results for supersonic

flow demonstrating the possibility of its use in capturing complex shock

structures. The aim of this study is to investigate the applicability,

accuracy, and efficacy of the multi-stage Runge-Kutta method in solving the

compressible Euler equations for two-dimenslonal supersonic internal flow

problems involving complex shock discontinuities of substantial strength. The

qualitative and quantitative nature of the effects of various techniques for

accelerating the convergence to a steady-state are analyzed in detail. The

results from the present calculations are compared with experimental data and

results obtained by other numerical methods [4,5].
-2-

II. SOLUTION ALGORITHM

The Runge-Kutta time-stepping scheme is more efficient than many other

explicit schemes such as MacCormack's original scheme. In addition, the

steady-state solutions are independent of the time-step, and as mentioned

earlier, the scheme is amenable to a variety of techniques for accelerating

steady-state convergence. Before the results are presented, we present a

description of these accelerating techniques.

Local Time-Stepping

A variable time-step determined by the bound on the local Courant number

is used in each zone. The use of the largest time-step at any location allowed

by the local stability bound accelerates convergence to the steady-state at

the expense of tlme-accuracy of the solution.

Implicit Residual Averaging

This strategy for accelerating convergence consists of taking an average

of residuals in neighboring cells by the following expression

_Ri-I + (l-2s)Ri + gRi+l = Ri (I)

where Ri denotes the residual in the i-th cell and Ri is the residual

before the implicit averaging. This is done after every even stage of the

scheme. This results in a scheme which is stable for a CFL number _ if

s _'_- - 1 (2)
-3-

where % is the CFL limit of the original scheme. Since % can be much

higher than X*, the rate of convergence is accelerated. However, choosing

too large a CFL limit, _, slows down the convergence. An optimal choice is

X/X _ 2. See reference [2] for further details.

Artificial Dissipation

The method requires artificial dissipative terms in order to alleviate

two difficulties associated with centrally differenced schemes. The first one

is that of odd-even point decoupling that adversely affects convergence.

Addition of a fourth-difference artificial viscosity reduces this difficulty.

The second difficulty manifests itself in oscillations in the vicinity of

shocks. An additional viscous term is introduced, which is an approximation to

a second derivative in order to overcome this problem. To maintain secondorder accuracy this term contains a coefficient that itself depends on the

second difference of the pressure. Since the fourth-order difference can

create oscillations near the shock, it is turned off in the vicinity of

shocks.

Enthalpy Damping

The enthalpy damping procedure has been modified from the original

version [i] of FLO52 in order to enable it to be used in supersonic flow.

Hence, we shall describe the new enthalpy damping in more detail.

In solving the full potential equation it is well known that adding

artificial terms such as #t can accelerate the convergence to a steady

state. For a first-order system of equations such as the Euler equations,

this can not be done. Instead, it is suggested [I] that one can add forcing
-4-

functions that depend on h-h0, h = (E+P)/o; where h and h 0 denote the

local and the free stream values of the total specific enthalpy. In the

steady state, the total specific enthalpy h is constant along each

streamline in an inviscid flow. Assuming that the flow comes from a common

reservoir, in the steady state h = h0 everywhere, where h0 is determined

by the inflow boundary conditions. Based on the analogy with the full

potential equations and the fact that h - h0 = _t' it was suggested that the

forcing function be added to the density equation. Using density O and

velocity components (u,v) the modified equations are

Ot + L 1 = -ap(h - h0) (2a)

u t + L2 = 0 (2b)

v t + L 3 = O, (2c)

where Lj are the usual space derivatives of the Euler equations. When one

solves for the momentum variables pu and pv, then the forcing function

also occurs in these equations so that

(pu) t + L2m = -apu(h - h0) (3a)

(PV) t + L3m = -apv(h - h0) (3b)

where the subscript m denotes the space derivative terms in the momentum

equations. The finite difference approximation to the density equation is
-5-

n+l n

- p = -At L 1 - sat pn+llhn',- h0_

P (4)

or

n

n+l P - At L 1

p = . (5)

i + _t(h n - ho)

Freezing coefficients and changing variables [2], one can reduce (2) to a

second-order equation for p,

P_ - (c2 - u2)pqq + 2KpT = 0 (6)

where q and T are new independent variables and the positive constant K

depends on i - M 2 and linearly on a. Thus, we obtain the potential-like

equation with the artificial term Kp T. To see the effect of K, we assume

that p has the simple form

p = e8T ei_q, -_ < n < _. (7)

Assuming subsonic flow, letting q2 = c2 _ u 2 and substituting (7) into (6),

we find

B2 + 2K8 + q2 _2 = 0

or

8 = -K + / K2

2 2

- q _ . (8)

In order to reach a steady state as fast as possible we wish the real part of

8, i.e, Re 8 to be as negative as possible. Note for K = 0, i.e., no
-6-

enthalpy damping, Re B = 0 and we have no term driving the system to a steady

state. We see from (8) that Re B is most negative when 0 < K < q_, and

K is as large as possible. Thus, we wish to choose K slightly less than

q_.

We first consider the case in which we want to damp the low frequencies,

since these are the slowest to decay. In this case, we have _ = 0(i) and

so K = 0(I). We stress that the constant HM in the FL052 code is equal to

cat. Hence, if _ = 0(i), this implies that HM should vary linearly with At

and therefore, is very small. We next consider the case where we wish to damp

the high frequencies, e.g., multigrid applications. Now _ = T/An and so

K varies as i/An or equivalently I/At. In this case, HM = cat should be

a constant.

This analysis implies that the enthalpy damping suggested in [I] is

useful if:

(I) the flow is locally subsonic;

(2a) HM is proportional to At to damp low frequencies;

(2b) HM is constant to damp high frequencies.

In order to explain the effectiveness of the enthalpy damping for nonmultigrid applications and to extend its usefulness to supersonic applications

we must discuss the energy equation which was ignored in (2). Since the

enthalpy damping was originally based on an analogy with the full potential

equation there is no obvious way to treat the energy equation. The original

idea was to append to (2) the equation

_s + = 0, (9) _--t-L4s
-7-

i.e., nothing special is done to the entropy equation. Rewriting (9) for the

energy equation we obtain

_-_ BE + L4E = - a(h - h0)h.

However, now the forcing function contains a quadratic term in h and was

found to be destabilizing rather than stabilizing. To avoid this difficulty

an ad hoc correction was used and the energy equation became

BE + = -a(h . B--t L4E - ho) (2d)

Instead of the ad hoc procedure, we shall replace (9) by

Bh = -8(h . -_ + L4h - h0) (9")

Thus, instead of ignoring the entropy we try and drive the total enthalpy to

its steady state value. We note that if B = = we recover the isoenergetic

approximation in which the energy equation is replaced by Bernoulli's law.

Thus, our new system is a compromise between the original Euler equations and

the isoenergetic equations. Of course, in the steady state all of these

variations are equivalent. Eliminating h we can rewrite (9") as

BE = - (eE + B0 (h . (2d')

+ L4E _--) - ho)

We thus, conclude that enthalpy damping helps in the non-multigrid version of

FLO52 because of the energy equation and not because of the continuity
-8-

equation, i.e., we can set s = 0 in (2) and the coefficient B assumes the

important role. This improved enthalpy damping can be used equally well in

supersonic regions. Furthermore, one can choose B as large as one desires

without causing instabilities, while choosing too large an sat causes

instabilities. However, choosing B too large may slow down the convergence

rate. In all the cases presented in the results section we chose s = 0 and

BAt to be constant.

Combining all the previous discussions we prescribe for the system (2a),

(3a), (3b), and (2d') that:

(I) For locally supersonic flow we must have s = O. For locally subsonic

flow we advise choosing s to depend on 1 - M 2. The coefficient B in

(2d') can be used for both subsonic and supersonic flows and does not depend

on the Mach number.

(2) To dampen the low frequencies, HM = sat should depend linearly on

At while to damp the high frequencies, HM = sat should be chosen to be

constant. The choice of B is frequency independent, but BAt should be

constant. Hence, even for a multigrid version one should choose sat and

BAt independently, especially since s should depend on 1 - M 2.

III. RESULTS

A wide variety of test cases involving different Mach numbers and shock

strengths were investigated in order to assess the ability of the numerical

scheme to simulate flows ranging from subsonic to high supersonic Mach

numbers. Convergence histories for the different test cases are presented

with a view to establishing the nature and extent of the effects of the
-9-

different convergence accelerating techniques as well as different grid

sizes. The CPU-seconds execution times presented are those for scalar

processing on a CDC Cyber-175 processor. Computed results have been compared

with experimental data and results from other schemes whenever available.

Test Case !

The first case investigated is the inviscid shock reflection problem in

which an incident shock wave of prescribed strength reflects off a flat

plate. A schematic of the configuration is presented in Figure I. Uniform

grids with various numbers of points in the x-dlrectlon (0 < x < 4.0) and

the y-direction (0 < y < 1.5) were considered. At the inflow boundary for

y > i and along the upper boundary all dependent variables are specified

based on a prescribed incident shock-angle _. All variables are extrapolated

from the interior at the supersonic outflow boundary.

The first set of computations for this configuration was performed with a

grid of 51 points in the x-dlrection and 31 points in the y-dlrectlon. The

shock-lncidence angle _ was 29° . A contour plot of the computed pressures

appears in Figure 2. Both the incident and the reflected shocks are seen to

be reasonably sharp. In Figure 3, a comparison of the pressure profile at the

y = 0.25 line computed by the present method with results obtained by the

original MacCormack scheme [4] and the Osher [5] upwind scheme is presented

and a good agreement of the results is noticed. Next, a set of convergece

histories corresponding to the use of the different convergence accelerating

techniques, is presented in Figure 4. The use of a locally maximum time step

is seen to substantially increase the rate of convergence over that obtained

by using a global minimum time step. The use of enthalpy damping does not
-i0-

create any advantage for this particular configuration; however, the

beneficial aspects of enthalpy damping will be apparent in the results

presented for the next test case. Finally, the use of the implicit residual

averaging technique results in a slight further increase in the convergence

rate. In all these plots the logarithm of the residual, i.e., the average

value of the time derivative of density, is plotted against actual CPU seconds

of computation.

The next series of computations for this configuration were done in order

to ascertain the effect of the grid aspect ratios on the rate of convergence.

Seven uniform grids of different sizes were used to compute the flow with

M = 2.9 and _ = 29° and the correspondinormack scheme [4] and the Osher

[5] upwind scheme is presented and a good agreement of the results is

noticed. Next, a set of convergece histories corresponding to the use of the

different convergence accelerating techniques, is presented in Figure 4. The

use of a locally maximum time step is seen to substantially increase the rate

of convergence over that obtained by using a global minimum time step. The

use of enthalpy damping does not create any advantage for this particular

configuration; however, the beneficial aspects of enthalpy damping will be

apparent in the results presented for the next test case. Finally, the use of

the implicit residual averaging technique results in a slight further increase

in the convergence rate. In all these plots the logarithm of the residual,

i.e., the average value of the time derivative of density, is plotted against

actual CPU seconds of computation.

The next series of computations for this configuration were done in order

to ascertain the effect of the grid aspect ratios on the rate of convergence.

Seven uniform grids of different sizes were used to compute the flow with
-IiM = 2.9 and _ = 29° and the corresponding convergence histories are

presented in Figure 5. The different aspect ratios of the cells are seen to

have pronounced effects on the convergence rates. The fastest rates of

convergence were generated by three particular grids, namely those with

41×25, 51×31, and 61×37 points. These three grids have a cell aspect ratio

of approximately 1.62. A departure from this value of the cell aspect ratio

is seen to substantially slow the rate of convergence. This indicates that

the convergence of the finite-volume Euler scheme is strongly influenced by

the aspect ratio of the computational cells. We next consider the possibility

that the angle of the shock may have some influence on the value of the cell

aspect ratio for which the fastest rate of convergence is attained. To test

this hypothesis an additional series of computations were performed for the

same configuration with the shock angle _ changed to 23° , all other

parameters remaining the same. The resulting convergence histories are

plotted in Figure 6. As is immediately obvious, the grid sizes with cell

aspect ratio of 1.62 again give rise to the fastest convergence rates. A

departure from this value, as in the case of the 51x43 grid, is again seen

to result in substantial worsening of the rate of convergence. This leads to

the conclusion that moderate changes in the shock angle do not substantially

influence the value of the optimal cell aspect ratio for which the fastest

convergence rate is achieved.

Test Case 2

The next test cases involved internal flow in a single throat nozzle [6].

Flow regimes ranging from subsonic to highly supersonic in nature have been

investigated. A variety of complex shock-structures have been successfully
-12-

captured. At supersonic inflow boundaries all variables are specified and at

supersonic outflow all the variables are extrapolated from the interior. For

subsonic inflow boundaries characterlstlc-type boundary conditions are

applied. At outflow boundaries characteristic extrapolation is supplemented

by a radiation boundary condition. A uniform gird of 61x31 points is used

in all of these computations.

A fully supersonic internal flow is considered as the first test problem

with this configuration. In this problem Mach 3.5 flow enters the twodimensional nozzle which is symmetric about the centerline. A contour plot of

the computed pressures is presented in Figure 7. A complex shock structure

consisting of multiple-reflected shocks is noticed. Interaction of the

expansion fan emanating from the nozzle wall with the shock reflected off the

centerline is clearly visible in the throat region. This demonstrates the

ability of the present Euler method to capture repeatedly reflected shock

structures at moderately high supersonic Mach numbers. The effects of the

various convergence accelerating strategies are apparent in the convergence

histories plotted in Figure 8. Enthalpy damping, in this case substantially

increases the rate of convergence. The use of implicit residual averaging is

seen to remarkably accelerate convergence and causes an approximately 200

seconds reduction in the computational time required to reach the same level

of the residual as that in the case of enthalpy damping alone.

The next problem provides a more severe test of the solution procedure's

ability to capture shocks and to handle subsonic and supersonic flow regions

in the configuration. The incoming flow is subsonic at a Mach number of

0.46. The flow accelerates as it travels down the nozzle and becomes supersonic at the throat creating a supersonic region in the middle section of the
-13-

nozzle. This supersonic region ultimately terminates at a normal shock beyond

which the flow is again subsonic. The terminal shock is clearly visible in

the pressure-contour plot presented in Figure 9. The presence of these

different flow regimes do not necessitate any special treatments as far as the

algorithm is concerned. A comparative plot of surface pressure computed by

the present algorithm and available experimental data [6] is presented in

Figure i0. The agreement of the results is noticeably good and confirms the

ability of the method to predict complex internal flows accurately. Some

experiments at hypersonic Mach numbers have been performed by the authors and

the results obtained generated encouragement as to the applicability of the

method to computation of hypersonic flow.

IV. CONCLUSIONS

Enthalpy damping, heretofore restricted in its applicability to subsonic

flows, has been demonstrated to be successfully applicable to supersonic

flows. Significant savings in computational time have been achieved by the

use of the present enthalpy damping technique. Problems involving supersonic

and subsonic flow regions have been solved successfully demonstrating the

versatility of the Runge-Kutta scheme. A particular value of the cell aspect

ratio is seen to cause the fastest rate of convergence; however, moderate

changes of the shock angle are not seen to have significant influence on this

value of the aspect ratio.
-14-

_FE_N_S

[I] Jameson, A.; Schmidt, W.; and Turkel, E.: Numerical solutions of the

Euler equations by finite volume methods using Runge-Kutta time-stepping

schemes. AIAA Paper No. 81-1259, 1981.

[2] Turkel, E.: Acceleration to a steady state for the Euler equations.

Numerical Methods for the Euler Equations of Fluid Dynamics, SLAM, 1985,

pp. 281-311.

[3] Moitra, A.: Numerical solution of the Euler equations for hlgh-speed,

blended wlng-body configurations. AIAA Paper No. 85-0123, 1985.

[4] Kumar, A.: Numerical analysis of the scramjet inlet flow field by using

two-dlmensional Navier-Stokes equations. NASA TP-1940, 1981.

[5] Chakravarty, S. R. and Osher, S.: High resolution applications of the

Osher upwind scheme for the Euler equations. AIAA Paper No. 83-1943,

1983.

[6] Boger, T. J.; Sajben, M.; and Kroutil, J. C.: Characteristic frequency

and length scales in transonic diffuser flow oscillations. AIAA Paper

81-1291, 1981.
-15-

M

2.9 _ _9° 1.5

Y- 0.25

Y

4.0

Figure I. Configuration of 2-D shock-reflection problem.
-16-

Figure 2. Contour plot of computed pressure for 2-D shock-reflection.

M = 2.9, c_= 29° .
-17-

5- Grid Size

51x 21

-- A Al

4 - 61x 21 _'-_- " -

31x 21 _,I

3

/i

P _=--_,,_/_/ Scheme

...... . ./J • • • OsherScheme

1 _ --- MacCormack

I 1 l [ ! I

.67 1.33 2.00 2.67 3.33 4.00

Figure 3. Plot of computed pressure at y = 0.25 for 2-D shock

reflection. M = 2.9, = = 29° .
-18-

Grid size51×31

0 _ATmin

o o ATlocal

• [] [] EnthalpyDamping

-I.0 _ ¢ ResidualAveraging

-1.5 "

20 '_

-2.5 _,*,'3. _

-3.0 o_o "JL.,

o

-3.5

-4.0 \\

-4.5

-.5 0 ,,I......,,,I,,,,,......... ,,,,I I......... I......... I

0 I00 200300400500600700800q00I000

CPUSEC

Figure 4. Convergence histories for different accelerating techniques for 2-D

shock-reflectlon. M = 2.9, _ = 29° , Grid: 51 x 31.
-19-

Grid size

31x 31

o o 41x31

[] [] 51x 31•

0 O <> 61x31

/x A 51X43

-.5 r,, _ 51X52

-I. 0 _ E361X37 _

-1..5 C_ C_41x 25_

-2.0 ,- t/_ r,

_,_ , .

-4.0 '_

-4..5

-5.0

0 100200300400500600700800900I000

CPU SEC

Flgure 5. Convergence histories for different grid sizes for 2-D shock

reflection. M = 2.9, _ = 29° .

* denotes cell aspect ratio = 1.62.
-20-

0 Grid size

-, 5 -- 41× 25*

o o 51x 31*

-I. 0 [] [] 61x 37*

-I. 5 _ _ 51x 43

-2,0 _

_: -2.5

 -3.o

-3.5 o

-4.0 -4.5

-5.0

0 I00200300400500600700800gO0I000

CPU SEC

Figure 6. Convergence histories for different grid sizes for 2-D shockreflection. M = 2.9, = = 23° .
-21-

Figure 7. Contour plot of computed pressure for 2-D nozzle. Min = 3.5.
-22-

400 600

CPU Sec

o

o

o

o

o

o

o

o

o

<)

o

-- ~T Min

o 0 ~T local

o 0 Entha Ipy Da mpi ng

a <> <> Residual Averaging

.5 ~~~~-­

1.0

1.5

V)

. 2. a

Q.)

0::: 2.5 .

.33. a

3.5

4.0

4.5

5. 0~w.uw.u.Uw.w.u.L.1w.w.&.w.w...w..u.L..u.u.~~~u.whw.wW

a 200

Figure 8. Convergence histories for different accelerating techniques for

2-D nozzle. Min = 3.5.
-23-

Figure 9. Contour plot of computed pressure for 2-D nozzle. Min = 0.46.
-24-

182

• Experment i

-- Euler

1.36

P 91 " =

,,, ,,*., |l,.,,.,ll, I**. ,, I i , *|,,,, tii ,, |, ,,,

".17 .04 .08 .20 .32

Figure I0. Comparison of computed and experimental surface pressure for 2-D

nozzle. Min = 0.46.

1. Report No. NASA CR-178053 12. Government Accession No. 3. Recipient's Catalog No.

ICASE Report No. 86-5

4. Title and Subtitle 5. Report Date

APPLICATION OF A RUNGE-KUTTA SCHEME T lQH~

FOR HIGH-SPEED INVISCID INTERNAL FLOWS 6. Performing Organization Code

7. Author(s} 8. Performing Organization Report No.

Anutosh Moitra, Eli Turkel, and Ajay Kumar A/;._"

9. Performing Organization Name and Address 10. Work Unit No.

Institute for Computer Applications in Science

and Engineering 11. Contract or Grant No.

Mail Stop 132C, NASA Langley Research Center NASl-17070j NASl-18l07 TT UI1 ?'l/;./;.<;_"??<;

12. Sponsoring'Agency Name and Address 13. Type of Report and Period Covered

Contractor Report

National Aeronautics and Space Administration 14. Spo"nsoring Agency Con"

tJ:lqh;nat-nn n r. ?()"L..~ 505-31-8 1 -01

15. Suppleme'iitary Notes .-

Langley Technical Monitor: AIM Paper No. 86-0104

J. C. South AIM 24th Aerospace Sciences

Meeting, 1/6-9/86, Reno, NV

J;';n:ll

16. Abstract .

A multi-stage Runge-Kutta method is analyzed for solving the twodimensional Euler equations for external and internal flow problems.

Subsonic, supersonic and, highly supersonic flows are studied. Various

techniques for accelerating the convergence to a steady state are described

and analyzed. Effects of the grid aspect ratio on the rate of convergence are

evaluated. An enthalpy damping technique applicable to supersonic flows is

described in detail. Numerical results for supersonic flows containing both

oblique and normal shocks are presented confirming the efficiency of the

method.

17. Key Words (Suggested by Authors(s» 18. Distribution Statement

Euler, 34 - Fluid Mechanics & Heat Transfer

Runge-Kutta,

high-speed,

internal flow Unclassified - unlimited

19. Security Classif.(of this report} 120. Security Classif.(of this page} 21. No. of Pages 122. Price

Unclassified Unclassified 26 A03

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

NASA-Langley, 1986