NASA/TM-1999-208793

Entropy Splitting and Numerical Dissipation

H. C. Yee, M. I,Tnokur, and M. J. Djomehri

August 1999

https://ntrs.nasa.gov/search.jsp?R=20010047483 2018-05-15T07:00:39+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 peer-reviewed 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-1320
NASA/TM-1999-208793

Entropy Splitting and Numerical Dissipation

H. C. Yee

Ames Research Center, Moffett Field, California

M. Vinokur

Ames Associate

Ames Research Center, Moffett Field, California

M. J. Djomehri

Sverdrup Technology Inc.

Ames Research Center, Moffett Field, California

National Aeronautics and

Space Administration

Ames Research Center

Moffett Field, California 94035-1000

August 1999
NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076-1320

(301 ) 621-0390

Available fi'om:

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

(703) 487-4650
ENTROPY SPLITTING AND NUMERICAL DISSIPATION

H.C. Yee 1, M. Vinokur 2 and M.J. DjomehrP

NASA Ames Research Center, Moffett Field, CA 94035, USA

Abstract

A rigorous stability estimate for arbitrary order of accuracy of spatial central difference schemes

for initial-boundary value problems of nonlinear symmetrizable systems of hyperbolic conservation

laws was established recently by Olsson and Oliger (1994) and Olsson (1995), and was applied

to the two-dimensional compressible Euler equations for a perfect gas by Gerritsen and Olsson

(1996) and Gerritsen (1996). The basic building block in developing the stability estimate is

a generalized energy approach based on a special splitting of the flux derivative via a convex

entropy function and certain homogeneous properties. Due to some of the unique properties of the

compressible Euler equations for a perfect gas, the splitting resulted in the sum of a conservative

portion and a non-conservative portion of the flux derivative, hereafter referred to as the "Entropy

• , Cy_, Sphttm_, . There are severaI potential desirable attributes and side benefits of the entropy splitting

for the compressible Euler equations that were not fully explored in Gerritsen and OIsson. The

paper has several objectives. The first is to investigate the choice of the arbitrary parameter that

determines the amount of splitting and its dependence on the type of physics of current interest to

computational fluid dynamics. The second is to investigate in what manner the splitting affects

the nonlinear stability of the central schemes for long time integrations of unsteady flows such as

in nonlinear aeroacoustics and turbulence dynamics. If numerical dissipation indeed is needed to

stabilize the central scheme, can the splitting help minimize the numerical dissipation compared to

its un-split cousin? Extensive numerical study on the vortex preservation capability of the splitting

in conjunction with central schemes for long time integrations will be presented. The third is

to study the effect of the non-conservative proportion of splitting in obtaining the correct shock

location for high speed complex shock-turbulence interactions. The fourth is to determine if this

method can be extended to other physical equations of state and other evolutionary equation sets.

If numerical dissipation is needed, the Yee, Sandham and Djomehri (1999)numerical dissipation

is employed. The Yee et al. schemes fit in the Olsson and Oliger framework.

t Senior StaR' Scientist

3 Ames Associate

SSenior Aerodynarnlcist_ Sverdrup Technology Inc. stationed at NASA Ames

Typeset by .AA4S-TEX
I. Introduction

The construction of efficient high order low dissipation numerical methods for nonlinear

conservation laws has been the subject of much research recently. For smooth flows, it is well

known that the standard high order non-dissipative central schemes generate spurious noise leading

to nonlinear instability, especially for long time integration applications such as in aeroacoustics,

rotorcraft dynamics and turbulence physics. On the other hand, central schemes in conjunction

with linear numerical dissipations are too diffusive for the physical problems in question. At the

same time the majority of the available high order high-resolution shock-capturing schemes are too

CPU intensive for practical computations. In spite of their high-resolution capability for rapidly

evolving flows and short term time integrations, for long time integrations these schemes often

exhibit undesirable amplitude errors for aeroacoustics problems. Current focus has been mainly

on algorithmic issues in constructing highly accurate methods away from boundaries. Rigorous

stability estimates for accurate and appropriate numerical boundary conditions associated with

fourth- or higher-order methods are equally important, and have been the major stumbling block

in the theoretical development of these schemes for nonlinear systems. Most of the existing theory

for nonlinear conservation laws and their finite discretizations is concerned with the initial value

problem (IVP). Standard practice in computational fluid dynamics (CFD) involving boundary

conditions relies on guidelines from theory for linear stability analysis of initial boundary value

problems (IBVPs) or IVP theories with the boundary conditions ignored. These iinearized stability

guidelines are only necessary but not sufficient for nonlinear stability. Spatial nonlinear stability of

IBVPs goes hand-in-hand with the appropriate amount of nonlinear numerical dissipation required

to stabilize the spatial scheme. The delicate balance of the numerical dissipation for stability

without the expense of excessive smearing of the flow physics after long time integrations poses a

severe challenge for unsteady flow simulations of this type. Employing a nonlinear stable form of

the governing equations in conjunction with the appropriate nonlinear stable scheme for IBVPs is

one way of minimizing the use of numerical dissipation.

Until recently it was not known how to derive the proper numerical boundary conditions for

a rigorous stability estimate for conventional spatial high order central difference schemes

for nonlinear hyperbolic IBVPs. Advances by Kreiss and Scherer (1977), Strand (1994) and

Olsson (1995a) led to the construction of high order boundary difference operators that enabled

the design of stable high order central schemes for linear hyperbolic systems. The major tool

used to overcome the stumbling block is a generalized energy method. The energy method for

deriving stability estimates for hyperbolic IBVPs was first applied to the nonlinear scalar case by

Gustafsson and Olsson (1994) for both higher-order central and compact symmetric schemes. It

was then generalized and extended to nonlinear systems of symmetrizable hyperbolic conservation

laws by Olsson and Oliger (1994) and Olsson (1995), and it was applied to the two-dimensional

(2-D) compressible Euler equations for a perfect gas by Gerritsen (1996) and Gerritsen and Oisson

(1996). With these recent developments, renewed interest has emerged in the use of spatial central

schemes where efficiency, simplicity and non-dissipative properties are their trademarks. See for

example, Yee et al. (1999) and Yee and Sandham (1998).
1.1. Relevance and Motivation

The basic building block in establishing a stable energy estimate for high order spatial central

schemes for nonlinear hyperbolic conservation laws consists of two parts. The first is a special

transformation of the conservation laws to an appropriate form for the application of the continuous

energy estimate for a stable IBVP of the governing equations. The second is a compatible numerical

boundary difference operator for the application of the discrete analogue of the continuous energy

estimate for a stable IBVP of the discretized counterparts. The special transformation relies on

the symmetrizability of the systems of nonlinear hyperbolic conservation laws, the possession

of a convex entropy function and a suitable splitting of the flux derivative vector with certain

homogeneous properties. The compatible boundary difference operator has to satisfy the discrete

analogue of the integration-by-parts procedure used in the continuous energy estimate. Olsson and

Oliger (1994) utilized the result of Harten (1983a) on symmetric forms for systems of conservation

laws as the backbone. Convexity of the flux functions is not required. These building blocks in

turn allow one to use a modified form of the energy estimate (or generalized energy estimate) in

deriving a compatible set of numerical boundary conditions that are stable for the higher-order

central differencing schemes.

Olsson proved that conservation is possible for second-order central schemes. To obtain a

rigorous estimate for higher-order central schemes, one must apply the scheme to the split form of

the flux derivative, written in non-conservative form, in terms of the transformed variables. For

the Euler equations, one can further simplify the final split form of the flux derivative leading to

a conservative portion and a non-conservative portion. The non-conservative portion as well as

the conservative portion can be recast in terms of the conservative variables, making it tractable

for practical applications in CFD. The proportion of the conservative to non-conservative parts

is dictated by a parameter which falls into two wide ranges of intervals. The resulting splitting

is hereafter referred to as the entropy splitting of the flux derivative or entropy splitting for

ease of reference. Here, the entropy splitting should not be confused with the traditional flux

vector splittings such as the Steger and Warming splitting (1981) or other variants. The traditional

flux vector splitting splits the flux function into different parts and most often into upwind and

downwind portions. However, the entropy spiitting splits the flux derivative using the properties

of the chosen entropy function and the symmetrizability of the conservation laws without reference

to any upwinding.

Harten showed that the viscous terms of the compressible Navier-Stokes equations can also be

symmetrized. In this case, only symmetry is needed in the derivation of the energy estimate. Due

to the parabolic nature of the boundary conditions, the homogeneity properties are not required for

the Navier-Stokes equations. For the numerical study involving the compressible Navier-Stokes

equations in the present study, we apply the entropy splitting to the inviscid fluxes and the

symmetric form of the viscous terms is not used. This is an attempt to examine if entropy splitting

of the inviscid flux derivatives alone will provide side benefits over the un-split approach.

Active research in the use of the symmetric form of the governing equations was carried out

by Hughes. Franca and Mallet (1986) and related recent work. Hughes et al. utilized only

the symmetric idea and employed the physical entropy as one of the entropy variables. Their
resulting inviscid flux vector and transformed state vector are not homogeneous in the entropy

variables. Unlike the entropy splitting, their transformed equations are in purely non-conservative

form. They have enjoyed improved results over the standard conservation law formulation. Their

approach, however, does not allow a rigorous stability estimate for IBVPs for nonlinear hyperbolic

conservation laws. In addition, due to their use of the purely non-conservative form, it is not

certain that a correct shock speed can be obtained in general.

The entropy splitting is not limited to smooth solutions. Olsson and Oliger (1994) also extended

their result to weak solutions (problems containing discontinuities) that are obtained as pointwise

limits of vanishing viscosity solutions. The entropy equality condition for the smooth solution

case now becomes an entropy inequality condition. In addition, appropriate numerical dissipation

is needed in conjunction with central schemes to pick out the physically relevant solutions if

weak solutions are present. Gustafsson and Olsson (1995) proposed a scalar filter as numerical

dissipation. Gerritsen and Olsson (1996) proposed the use of a slightly different nonlinear scalar

filter in conjunction with wavelets for sharp shock capturing and shock detection. The recently

developed high order low-dissipative shock capturing schemes using characteristic filters of Yee

et ai. (1999) fit in the entropy splitting framework. Instead of applying a scalar filter, they supply

nonlinear filters based on, locally, the different wave characteristics of the convective flux. For

complex shock waves, shear and turbulence interactions, one has better control of the amount

of dissipation associated with each wave. For efficiency, Yee et al. proposed a combination

of a narrow grid stencil higher-order non-dissipative classical spatial differencing schemes and

low order total variation diminishing (TVD), essentially non-oscillatory (ENO) or weighted ENO

(WENO) dissipations as nonlinear characteristic filters with an artificial compression method

(ACM) sensor. The ACM sensor is the same as that of Harten (1978) but applied in a slightly

different context.

The design principle of the Yee et ai. schemes consists of two steps. The first step is the high

order spatial and temporal base scheme. Unlike the hybrid schemes, the base scheme is always

activated. The high order central and compact symmetric schemes are their primary choice for

base schemes. In other words, the primary base schemes used by Yee et al. are exactly the central

schemes for which Olsson and Oliger (1994) have provided rigorous stability estimates. The second

step is the appropriate filter for stability and to capture shocks, shear-layers and fine scale flow

structure. Many of the TVD, positive, ENO and WENO dissipations, after a minor modification

by the ACM sensor, are suitable candidates as filters. The final grid stencil of these schemes

is five in each spatial direction if second-order TVD schemes are used as filters and seven if

second-order ENO schemes are used as filters for a fourth-order base scheme. There is only a 10%

increase in operations count over standard second-order TVD schemes for 2-D direct numerical

simulations (DNS). Based on their preliminary study for shock-turbulence flows, higher accuracy

was achieved with fewer grid points when compared with that of standard second-order TVD,

positive or ENO schemes. See Yee et al. (1999) or references cited therein. For the suppression of

unphysical high frequency oscillations due to insufficient grid resolution and nonlinear instability

away from discontinuities, higher-order spectral-like filters (Vichnevetsky (1974), Lele (1994),

Alpert (198 l), Visbal and Gaitonde (1998), Gaitonde and Visbal (1999)) might be needed where

the value of the ACM sensor is extremely small.
1.2. Objectives

Motivated by the aforementioned development of entropy splitting, Yee et al. (1999) proposed,

as a followup work, to apply their schemes to the entropy splitting form of the inviscid flux

derivatives. This paper is a sequel to Yee et al. Besides investigating some of the fundamental

issues described below, studies will be conducted to determine to what extent the entropy splitting

form of the flux derivative can help in minimizing numerical dissipation, or equivalently, in

improving nonlinear stability in conjunction with the Yee et al. (1999) schemes. Our main goal is

to explore the possible side benefits of the entropy splitting without considering the accompanying

stable numerical boundary difference operator as a complete package for stability requirements.

This is accomplished by choosing numerical examples with periodic boundary conditions, or

computational domains whose boundaries are far enough away so as not to affect the mainstream

flow activities, and/or by using lower order non-characteristic boundary schemes.

There are several potential desirable attributes of the entropy splitting for the compressible Euler

equations that were not fully explored in Gerritsen and Olsson. First, in regions of smooth flows,

additional numerical dissipation might not be required by the entropy splitting in conjunction

with non-dissipative spatial central difference schemes. Second, the splitting appears to improve

nonlinear stability over the un-split approach employing the same non-dissipative higher-order

central schemes even for periodic boundary conditions. Third, the non-conservative portion of

the flux derivative seems to have a small effect in obtaining correct shock speeds on the physical

problems that Gerritsen and Olsson considered. Fourth, the entropy splitting in conjunction with

higher-order central differencing could be a good candidate for nonlinear aeroacoustics, rotorcraft

dynamics and turbulence computations where simplicity, high accuracy and low numerical

dissipation are essential. But most of all, the splitting could possibly be extended to other physical

equations of state and other evolutionary equation sets. The objective of this paper is many

fold. The first is to investigate the choice of the arbitrary parameter that determines the amount

of splitting (conservative and non-conservative proportions) and its dependence on the type of

physics of current interest to CFD. The second is to investigate in what manner the splitting affects

the nonlinear stability of the central schemes for long time integrations of unsteady flows such as

in nonlinear aeroacoustics and turbulence dynamics. If numerical dissipation indeed is needed to

stabilize the central scheme, can the splitting help minimize the use of numerical dissipation over

its un-split cousin? Extensive numerical study on the vortex preservation capability of the splitting

in conjunction with the central schemes for long time integrations will be presented in section

4.1. At present, existing finite discretizations (applied to the un-split approach) that are capable

of preserving long distance vortex convection are ineffective and CPU intensive. Extensive grid

refinement and grid adaptation, and an unusually demanding time step reduction are necessary.

The third is to study the effect of the non-conservative proportion of the splitting in obtaining

the correct shock location for high speed complex shock-turbulence interactions. The fourth is to

determine if this method can be extended to other physical equations of state and evolutionary

equation sets.

Aside from stability considerations, as explained in Harten (1983a), another potential desirable

attribute in the use of the symmetric form of the governing equations is for the computation of the
steady-statesolutionof theconservationlaws. In solvingthesteadynonlinearconservationlaws,

thesymmetryof thematrix coefficientscouldpossiblyenhancethestructureof iterativematricesin

directNewton-iterationmethods.Fortime-marchingto steadystatesor time accuratesubiteration

proceduresusingimplicit methods,thesymmetricform in conjunctionwith theentropy splitting

mightresultin animprovedconvergencerateovertheun-splitapproach.This will beasubjectof

futureresearch.

Outline:

Section II reviews the entropy splitting and the numerical schemes for the 2-D compressible

Euler equations for a perfect gas. The choice of the entropy splitting parameter is discussed in

Section 2.2. Section III describes the extension of the entropy splitting to other physical equations

of states and evolutionary equation sets. Section IV illustrates the performance of the entropy

splitting for a variety of unsteady flows and compares the results with those obtained using the

un-split conservative approach (i.e., without splitting the inviscid flux derivatives). The study

concentrates only on perfect gases. The comparison uses the same spatial and time discretizations

for both the split form and un-split form of the inviscid flux derivatives. If numerical dissipation

is required, the Yee et al. (1999) filters are used as the numerical dissipation for nonlinear stability

and/or to insure physically relevant numerical solutions and to suppress spurious oscillations across

discontinuities. If spurious high frequency oscillations are present, the use of spectral-like filters

(Vichnevetsky (1974), Lele (1994), Alpert (1981), Visbal and Gaitonde (1998), Gaitonde and

VisbaI (1999)) in conjunction with the Yee et al. filters might be needed at the locations where the

value of the ACM sensor is very small. See Section 2.6 for a discussion.

In this paper, unless indicated, Euler or Navier-Stokes equations pertain to compressible fluids.

High order central difference schemes refer to fourth or higher-order spatial central difference

schemes (compact or non-compact methods) without numerical dissipation added. Compatible

time discretizations (in terms of stability and accuracy) should be used, but these are not the subject

of this paper. The terms "split" and "un-split" mean the application of the same discretization to

the "entropy splitting of the inviscid flux derivative" and the standard "inviscid flux derivative

in terms of the conservative variables without splitting".

II. Entropy Splitting for a Perfect Gas

This section reviews the basic building blocks for the entropy splitting for the 2-D compressible

Euler equations for a perfect gas. The mathematical theory is quite involved. Interested readers are

referred to references cited. The Yee et al. (1999) numerical methods used in conjunction with the

entropy splitting are also summarized.

2.1. Summary of Entropy Splitting for a Perfect Gas

In vector notation the 2-D compressible time-dependent Euler Equations in conservation form

for an equilibrium real gas can be written, in Cartesian coordinates, as
Ut + F. + Gy = O, (2.1.Ia)

where Ut = ou OF za -_-, F. = _- and G_ ="_d and the U, F, G, are vectors given by

U = [!1; F = I +1 pu p puv

m'" ; C = |_, + P/" (2.1.1b)

The dependent variable U is the vector of conservative variables, and (p, u, v,p) T is the vector of

primitive variables. Here p is the density, u and v are the velocity components, pu and pv are the

• - and V-components of the momentum per unit volume, p is the pressure, e = pie + (u s + v_)/2]

is the total energy per unit volume, and e is the specific internal energy.

For a thermally perfect gas, the equation of state is

p = pRT, (2.1.2)

where R is the specific gas constant, and T is the temperature with _ = t(T). For constant specific

heats (calorically perfect gas)

e - c,,T, (2.1.3)

where c,, is the specific heat at constant volume.

The eigenvalues associated with the flux Jacobian matrices of F and (7 are (u, u, u + e) and

(v,v,v -t-c), where c is the sound speed. The two u,u and v,v characteristics are linearly

degenerate. Hereafter, we refer to the fields associated with the u + c and v + c characteristics as

the nonlinear fields and the fields associated with the u, u and v, v characteristics as the linear

fields.

Here we outline the entropy splitting for a perfect gas in 2-D Cartesian coordinates. Formulas for

the corresponding 3-D case can be found in Appendix B. Gerritsen and Olsson (1996) extended the

summation-by-parts idea of Strand (1994), and the entropy splitting of Olsson (1995) and Olsson

and Oliger (1994) to the 2-D Euler equations for an ideal gas in conjunction with high order central

schemes. The first step in deriving the entropy splitting for the compressible Euler equations for a

perfect gas is to introduce a symmetry transformation from the vector of conservative variables U to

a new vector of symmetry variables W, referred to as the "entropy variables". The transformation

is chosen so that the flux Jacobian matrix with respect to W is symmetric, and the Jacobian matrix

of the transformation is symmetric and positive definite. A family of symmetry transformations,

based on a scalar convex function ,7, referred to as an "entropy function", derived for the Euler

equations for a perfect gas by Harten (1983a), was employed by Gerritsen and Olsson. It has the

form
h

r1 -" p_(S). (2.1.4)

The function _(S) is an arbitrary but differentiable function of a dimensionless physical entropy

S = log(pp -'r) (2.1.5)

where S has been non-dimensionalize by c,. The entropy variables W are then given by W = ea--_u.

The entropy function r/is not to be confused with the "physical entropy" S or the entropy variables

vector W. The next step is to restrict the transformations to those that allow a special splitting of

the flux derivatives. This requires that both the flux and conservative variables are homogeneous

functions of the new variables. Under these conditions one can rigorously establish a bound on

the rate of growth of the energy norm in terms of the absolute eigenvalues corresponding to the

in-coming characteristic variables at the boundary of the domain. The relation between entropy

functions and symmetrizable variables can be found in Mock (1980).

In other words, one introduces an entropy variable transformation W = W(U) such that

Fw = or' and Ow = _oa #tr are symmetric, and Uw = _ is symmetric and positive definite.

The entropy variable W is chosen such that F(U(W)), G(U(W)) and U(W) are homogeneous

functions of W of degree/_, i.e., there is a constant _ such that for all r

The homogeneity property implies that

FwW =/_F(U(W))

UwW = _U.

Then the splitting of F. results in

F. - _3 +1 l ( Fw W ) " + -_----_ 1 Fw W, - t3 /3+ -_F,. + _-_-_FwW,, t3 # -1.

A similar splitting can be written for G u and Us.

For a perfect gas, the required entropy function is obtained by letting t_(S) = KerZ-_

K and a are constants. The corresponding W can be written as

(2.1.10)

, where
a-I w = [,.1 ,,,, ,,,, w,I_"= P-[. ÷ _-___p-m

p

and the upper triangular part of the symmetric matrix Uw is

(2.1.11a)

UW = --1 [ap

p" [

apu apv

apu 2 - p apuv

apv 2 - p

_p(u _ + v 2) - 7_-_y__p

u[_p(,,,+_,) - bp]

_[_p(,,,+ _,)- bp]

_ - t,r,(,,'+ ,,') + _p(,,' + ,,_)' 3,--1 p

(2.1.11b)

The constants a and b are a = (1 - a - 7)/a and b = 7/(7 - 1). Here, p, and p are related

through

with X = -_ < 0. In the authors' opinion, the simplest choice is to set K =/3. The parameter

is given by

/3'= +--------_, a (2.1.11d)

1-7

Using (2.1.1b), (2.1.2), (2.1.3), (2.1.1 la) and (2.1.11d), we can show (see Section 3.1) that U, F

and G are homogeneous functions of W of degree/3. The positive definite condition on Uw (see

Section 3.1) restricts a to the two ranges a > 0 ora < -7.

The flux vectors, expressed in the W variables are given by

p [_,_,_ +p" _ _(=1 + .,_,pJ;, F(tr(W))p" ,,,. .,,.

(2.1.11e)

a(u(w))= p'--P [_,,,,_,..-__e.. +P" __,,(=_ +____?p.)]T.(2.z.110

The upper triangular parts of the symmetric matrices F(U(W))w and G(U(W))w, expressed in

the U variables, are given by

apu apu 2 - p apuv

,,(a_, 2 - 3r,) ,,(,,p,,_ - r')

"('_' - I")

,,[_p(,_, + ,,2)_ t,r,]

Ii

it'(" + v2)

._[_p+ _p(,,' + ,,')]

,,[_ +_,,(u,+ ,,,)+ _,,(,,,+,,')']

(2.1.11g)
10

apv apuv apv 2 - p

Cw = 1 "("P"2 - P) "("_' - P)

p-"_ v(apv 2 - 3p)

where _ = (1 - 23')/(7 - 1).

3 I 2 101

(2.z.zth)

2.2. Choice of the Entropy Splitting Parameter

From the structure of (2.1.10), the entropy splitting divides the flux derivative into a conservative

and a non-conservative part. The ratio between the conservative and non-conservative parts depends

on the choice of the parameter/3. Both Harten (1983a) and Gerritsen and OIsson (1996) introduced

the parameter at. In the authors" opinion, the introduction of at is not necessary. However, to adhere

to the discussion when referencing their work, we retain the use of at in the perfect gas case. The

convexity condition on the entropy function r/restricts the value at to two possible ranges; namely,

at > 0 orat < -7 (or equivalently, /3 < _ or/3 > 0). Although Gerritsen and Olsson considered

the at < -7 range which Harten (1983a) overlooked, they set at = I - 23' (/9 = 1) in conjunction

with high order central differencing schemes in all of their numerical examples. This choice of

at corresponds to the splitting of the flux derivative into equal conservative and non-conservative

proportions. They did not give any guidelines or examples of the effect of the choice of at on the

quality of the numerical solutions for different flow physics. In addition, all of their examples deal

with at most simple shock waves, if present. Wavelets are used as shock detectors and to guide

the grid adaptation. Due to the type of problem they addressed and the dense clustering of the grid

points near the shocks using very small time steps, it is not certain that correct shock speeds were

really obtained with a reasonably practical time step and grid distribution. It is the purpose of this

section to discuss the choice of the at parameter value. We discuss at > 0 and at < -7 separately.

Theat > 0 (or_ < I_--_) case:

This is the only case that Harten considered. This corresponds to a negative/3 which results in

a conservative proportion fro = -2- > I and a non-conservative proportion fr,,c -0.- V_* <

*+_

As at _ 0 ÷, f,= _ 3'- and f_',,= _ (1 - 3') +. Here, the superscripts "+" and "-" indicate

the values approach the limit from above and below respectively. Thus, it appears that at > 0

is "non-standard" or "nonphysical" in CFD practice in the sense that a larger than 100% of

the conservative proportion and a negative non-conservative proportion is used. As ct _ 0 +, the

proportion becomes extremely unphysical. As at _ oo, f"e _ 1 + and f_',¢ _ 0-. Therefore as

at _ oo, the proportion becomes more physical.

The at <-7 (orfl>O) case:

The at < -7 case corresponds to a positive/3 and consequently, .fr¢ < 1 and .f,,.,= < 1.

This range of a conforms more closely to standard CFD practices. We have the following five

situations.
i) Asa_-%_--,0 +,fee _0 +andfr.c _lii) Fora = 1 - 27,/9 = 1, fr, = 1/2 and fr,c = 1/2

iii) For I - 23, < a < -3', fr¢ < 1/2 and frn, > 1/2

iv) For a < 1 - 27, fr, > 1/2 and fr,,, < 1/2

v) As a _ -oo, fr_ _ 1- and fr,.,,. _ 0 +

Section IV gives a parameter study of a for three different types of flow physics.

2.3. Numerical Methods

ll

The spatial discretizations for weak solutions proposed in Yee et al. (1999) consist of two

parts, namely, a base scheme and a filter. When numerical dissipations or filters are not used,

the scheme consists of only the base scheme. This section discusses the base scheme and Section

2.4 discusses the form of the filter (numerical dissipation) for complex shock waves, shear and

turbulence interactions. Section 2,6 discusses the blending of the Yee et al. (1999) filters with

other filters for the suppression of spurious high frequency oscillations.

2.3.1. Spatial Base Schemes for the Convection Terms

Denote F./,k as the discrete approximation of the convection flux F at (jAz, kay), where A;e

and Ay are the grid spacing in the z- and y-directions and j and tt are the corresponding spatial

indices. Possible non-dissipative high order base schemes for F, (similarly for Gy) can be of the

following two types.

Central Differencings: (fourth and sixth-order)

F, _ _ 1(Fi+2,_ - 8Fj+1,k + 8Fj_,,_ - Fj-,,h ), (2.3.1)

1 / \

[Fj+s,/,- 9Fj+2,t, + 45Fj+Ij,- 45Fj_1,& + 9Fj-2,k- Fj-s,k}.

k /

(2.3.2)

Compact Central Differencings: (fourth and sixth-order)

F, ,-_\ / ./,I,

where for a fourth-order approximation

(A,,F)jj, = -_I (Fj+L. + 4Fj,. + Fj-I,.)

= 1(f - ),

(2.3.3a)

(2.3.3b)

(2.3.3c)
12

and for a sixth-order approximation

(A.F)j,_, = -_I(Fj+, ,h + 3Fj,_ + Fj__,h)

(B.F)j,h = _-_l (Fi+,j, + 28Fi+,j, - 28Fi-I,_,- Fj-2,k ) •

(2.3.3d)

(2.3.3e)

2.3.2. Spatial Schemes for Viscous Terms

For simplicity let V.. be a viscous term in one dimension. The possible high order base schemes

for V.. can be

Central Differencings: (fourth and sixth-order)

If.._ 12Az------ 1(_ 5+, - 165+I + 305 - 165-* + 5-, ),

(2.3.4)

1 / \

/2_+, - 27_+,+270v.,-490_+270vj_,-27v__,+2v__,/.(2.3.5)

180Az 2 \ /

Compact Central Differencings: (fourth and sixth-order)

V.. _ _ C2'D.V ,

Az2 i

where for a fourth-order approximation

(c.vb = i_1(D+, + lO5 + 5-, ),

(D.V)j = 5+* - 25 + 5-*,

and for a sixth-order approximation

(C-Vb = 5+* + a*5 + 5-*,

(D.Vb = bo 5+,- 25 + 5-, + ¥ 5+_- 2v_+ 5-_

a0 = 5.5,

bo = 4(ao - 1)/3,

co = (10 - ao)/3.

(2.3.6a)

(2.3.6b)

(2.3.6c)

(2.3.6d)

(2.3.6e)

(2.3.6t")

(2.3.6g)

(2.3.6h)
2.4. Filters

13

In this section we first review the procedure for applying the characteristic filter to multistage

Runge-Kutta type and linear multistep method (LMM) types of time discretizations (Yee et al.

(1999)). Examples of explicit LMMS are forward Euler and Adams-Bashforth. Examples of

implicit LMMs are backward Euler, trapezoidal rule, and three-point backward differentiation.

The one-leg formulation of the LMMs of Dahlquist (1979) is also applicable. We then discuss

forms of the characteristic filter.

2.4.1. Procedure to Apply the Filter Step

If a multistage time discretization such as the Runge-Kutta method is desired, the spatial

differencing base scheme discussed in the previous section is applied at every stage of the RungeKutta method. If viscous terms are present, we use the same order and type of base scheme for the

viscous terms as for the convection terms.

There are two methods for applying the characteristic filter. Method I is to apply the filter

at every stage of the Runge-Kutta step. Method 2 is to apply the filter at the end of the full

Runge-Kutta step. For inviscid and strong shock interactions, method 1 might be more stable.

If one desires a time discretization that belongs to the class of LMMs, then the filter can be

applied as a numerical dissipation vector in conjunction with the base scheme. The filter in this

case is evaluated at U" for explicit LMMs. For implicit LMMs additional similar filters evaluated

at the n + 1 time level are involved. Alternatively, method 2 can be applied to LMMs as well.

In this case, we apply the filter after the completion of the implicit time step. One can minimize

flux evaluations by using the one-leg formulation of the LMMs of Dahlquist (1979). The only

non-dissipative (in time) second-order, two-time level one-leg method is the mid-point implicit

method. Note that the noniterative linearized form of the midpoint implicit formula reduces to the

regular noniterative linearized trapezoidal formula.

For time marching to steady states using implicit LMMs, certain flow physics only requires an

explicit dissipation term. Also, the implicit operator can be different from the explicit operator.

See Yee (1985, 1986, 1989), Yee et al. (1990) for some efficient conservative linearized implicit

forms.

2.4.2. Nonlinear Characteristic Filters

There are many possible candidates for the filter operator in conjunction with high order base

schemes. Here, we propose using filter operators that have similar width of grid stencils as the base

scheme for efficiency and ease of numerical boundary treatment. Higher than third-order filter

operators are of course applicable, but they are more CPU intensive and require special treatment

near boundary points for stability and accuracy.

We use nonlinear dissipation terms in conjunction with the Harten ACM sensor applied to each

characteristic wave as the filter vector. [n essence, the nonlinear dissipation terms act as second or
14

third-order ACM-like operators instead of Harten's first-order ACM (Harten (1978)). The sensor

is used to signal the amount of nonlinear dissipation to be added to the high order non-dissipative

scheme, one wave at a time. Thus, the current approach is also different in spirit from using ACM

to sharpen the contact discontinuities in the original Harten (1983b) second-order TVD scheme.

Let the filter vector in the z-direction be of the form

(2.4.1)

Fj+ l,t' is the modified form of the nonlinear dissipation portion of the standard numerical flux.

For characteristic based methods, the quantity Rj+ ] is the right eigenvector matrix of _-UsFusing,

for example, Roe's approximate average state. Note that the eigenvector Ri+ 1⁄2 should not be

confused with the R in (2.1.2). We cast the G;,h+ 1⁄2 in the same manner. The elements of ,_*./+I

i

denoted by (4,./+ ] )" are

, _;0.,, ,_., , 1⁄2. (2.4.2) i $t" i It

_b_+, in (2.4.2) are the elements of _./+ 1⁄2 - the dissipative vector of the high resolution schemes

resulting from using a TVD, ENO or WENO scheme. Hereafter, we refer (2.4.2) as the ACM filter.

The function ,_0_+ 1⁄2 is the key mechanism for achieving high accuracy of the fine scale flow

structure as well as shock waves in a stable manner. In other words, the elements of ,I,"./+ [ are the

same as the nonlinear dissipation term of the TVD, ENO or WENO scheme with the exception

of premultiplying by t¢0_+1⁄2. The parameter _¢ is problem dependent. For smooth flows, _¢ is

used to improve nonlinear stability and can be very small. Different physical problems require a

different value of t¢ because of the large variation in flow properties. The _¢ value may vary from

one characteristic wave to another, and from one region of the flow field to another region with

different flow structure. The range of _ for our present numerical experiments is 0.0 < s _< 0.7.

The function 0./+l _ is the Harten ACM sensor. For a general 2m + 1 points base scheme, Harten

recommended

Otj+l = max (_-,.,+1, ..., _+_), (2.4.3)

[ %+It -

(2.4.4)

The a: -1 Uj). j+ 1⁄2 are elements of R./+ 1⁄2(U./+l -

Instead of varying _ for the particular flow problem, one can vary p. For largerp, less numerical

dissipation is added. Note that by varying p > 1 in (2.4.4), one can essentially increase the order

of accuracy of the dissipation term. The order of the dissipation depends on the value ofp (Bjorn
15

Sjogreen, private communication, 1998). One can switch from p = 1 near shock locations to p > 1

at smooth regions. For all of the numerical examples, we use p = 1 and

0_+1⁄2 = ma.x:: (Oj,Oj+,). (2.4.5)

The shock-turbulence interaction problems appear to favor this form of 0 l

j+1⁄2

!

Formulae for _b./+ i are well known and can be found in the literature. In most of the numerical

computations in Section IV, we use the Harten and Yee (1985) second-order upwind TVD

numerical dissipation. Computations using the symmetric TVD dissipation (Yee, 1985) will also

be presented. See Yee et al. (1999) for details and for a discussion of other possible filters.

2.5. Computer Implementation

To avoid some conditional statements in the actual computer code and to promote vectorization,

several of the functions inside the filter with the potential of dividing by zero are modified. See

Yee et ai. (1999) for details. In particular, the sensor (2.4.4), with p = 1, is replaced by

I' j÷,/21 t (2.4.6) -la _Zl21

Oj = l !

+ 1%_,/21 + e"

In all of the computations, we take e = 10 -r. (Actually, e should have the same dimension as

a_+xl2)

2.6. Blending ACM Filters with Other Filters

The ACM filter (2.4.2) might not be sufficient for (a) time-marching to steady state and (b)

spurious high frequency oscillations due to insufficient grid resolution and nonlinear instability

away from discontinuities, especially for turbulent and large-eddy simulations. This Section

discusses the blending of other filters with ACM filters.

(a) Time-Marching to Steady State:

For time-marching to steady state one usually needs to add fourth-order linear dissipation to a

second-order spatial differencing scheme (Beam and Warming (1976)). For the present schemes

using characteristic filters, in addition to the ACM filter, one might need to add a sixth-order

linear dissipation to a fourth-order spatial base scheme and an eighth-order linear dissipation to a

sixth-order spatial base scheme in regions away from shocks for stability and convergence. Let/_d

be such an additional filter operator. Take the case of a Runge-Kutta time discretization. There

are again two ways of incorporating the La operator. One way is to incorporate the La operator

at every stage of the Runge-Kutta method. Another way is to include the Ld operator after the

completion of the Runge-Kutta full step. The best way of applying the Ld operator is most likely

problem dependent and time integrator dependent. For LMM type of time integrators, L,t is used

in conjunction with the ACM filter step as an additional dissipation.
To minimizetheamountof dissipation due to L,/in the vicinity of shock waves, there should be

a switching mechanism '_d (different from _ in (2.4.2)) to turn off the Ld operator in the vicinity

of shock waves. The Ld operator can be applied to the conservative, primitive or characteristic

variables. The simplest form is to apply Ld to the conservative variables. Alternatively, since

all of the work in computing the average states and the characteristic variables is done for the

shock-capturing filter operator, one can apply the La operator to the characteristic variables. In

this case, the switching mechanism k,_ can be a vector so that it is more in tune with the ACM

shock detector using the approximate Riemann solver. For example, one can set _ = 0 for the

linearly degenerate fields and blend a small amount of _d to remove spurious noise generated by

the lack of ACM filters. This blending of the ACM filter with the Ld operator can be applied to

time-accurate computations as well.

(b) Suppression of Spurious High Frequency Oscillations:

The ACM filters might not be able to remove spurious high frequency oscillations effectively

unless sufficient grid points are used. For the suppression of unphysical high frequency oscillations

due to insufficient grid resolution and nonlinear instability away from discontinuities, higher-order

spectral-like filters (Vichnevetsky (1974), Lele (1994), Alpert ( 198 l), Visbal and Gaitonde (1998),

Gaitonde and Visbal (1999)) might be needed at the locations where the value of the ACM

sensor is very small or zero. If spectral-like filters are needed, a proper blending of ACM filters

with spectral-like filters should be applied. In this case, we can use the same procedures as the

time-marching to the steady state except the La operator should be replaced with the spectral-like

filters (for compact central schemes).

IlL Extension to Other Equations of State and Equation Sets

In this section we discuss the extension of entropy splitting to more general cases. The method

originally was developed for the 2-D Euler equations in Cartesian coordinates for a perfect gas.

We show here how it can be extended to flow of a gas that is only thermally perfect. For

maximum generality, the analysis is presented for arbitrary three-dimensional, time varying grids.

For compactness, we employ the vector approach of Vinokur (1989). Here the word vector

refers to a physical vector such as velocity or momentum, as distinguished from an algebraic

vector representing a set of variables. For completeness, the Roe Riemann solver for a thermally

perfect gas is included. This is motivated by the fact that if the characteristic type of nonlinear

dissipation or filter is desired, Roe's Riemann solver is normally employed. It is noted that

both Steger-Warming flux-vector splitting and Roe's approximate Riemann solver have an exact

extension for this case. For the readers unfamiliar with the vector approach of Vinokur (1989), the

results for the 2-D and 3-D Euler equations in Cartesian coordinates are given in Appendices A

and B. The caloric equation for an ideal diatomic gas is given in Appendix C. We also examine

the possible extensions to a nonequilibrium mixture of species, magnetohydrodynamics, and the

artificial compressibility method for incompressible flow.
17

3.1. Entropy Splitting for a Thermally Perfect Gas

In this subsection we derive entropy splitting for a gas that is only thermally perfect, with

the internal energy being an arbitrary function of temperature. This law is valid for a dilute gas

consisting of a single chemical species, and is also a very good approximation for air below the

temperature when oxygen starts to dissociate (approximately 2000* K).

The following development describes the derivation leading to a final form of the entropy

splitting of the flux derivative for a thermally perfect gas. It has the same form as the perfect gas

case, but, the corresponding ranges of the/3 parameter are different, and W, Uw and Fw have

different expressions. In fact, the derivation of entropy splitting for a perfect gas has to be modified.

Certain parameters that are constant for a perfect gas are no longer constant for a thermally perfect

gas. In particular, Harten and Gerritsen normalized their entropy S by c,,, which is now a variable

for a thermally perfect gas. We therefore normalize our S by the gas constant R. This results in

differences in our results from theirs involving the quantity 3' - 1. The positive definite condition

(or equivalently, the convexity condition) on Uw again restricts _ to be in two possible ranges. As

mentioned previously, Harten (1983a) overlooked one of the more physical ranges, and needlessly

introduced a parameter a in his solution. Such a parameter has no analogue for the more general

thermally perfect gas case, and only serves to complicate the derivation. Harten also introduced

an arbitrary constant K, which he then set equal to a particular value in order to simplify the

final expressions. We chose a particular value from the beginning, and avoided introducing an

un-needed constant.

We repeat the equation of state for a thermally perfect gas (2.1.2)

p = paT, (3.1.1)

where p, p, T, and R are the pressure, density, temperature, and gas constant, respectively. The

entropy S and internal energy e are then related to p and T through the first and second laws of

A

dS = de.., tip (3.1.2)

T P'

thermodynamics by

where we have introduced the normalized variables

S and f = RT. (3.1.3)

Equation (3.1.1) can then be rewritten as

(3.1.4)

While S is dimensionless, T has the same dimension as e. It follows from (3.1.2) that e -- _(T)

only. All real gases satisfy the conditions _ > 0 and _"> 0, where _ = de/dT.

Equation (3.1.2) can be integrated to obtain

p/= (3.1.5)
18

where

f(T) = exp(- / _dT). (3.1.6)

The arbitraryconstantin the integralof (3.1.6)can be absorbed inthe definitionof S.

Since conservation laws are expressed in terms of of conserved quantitiesper unitvolume, it

is convenient to introduce the internal energy per unit volume _" = p_. If u is the fluid velocity

vector, then the set of conservative variables U can be represented by the vector

G--[p m e] T, (3.1.7)

where m = pu is the momentum vector per unit volume, and e = _'+ ]pu. u is the total energy

per unit volume. Note the algebraic vector U has three elements, of which the second element is

the physical vector m. The temperature T(U) is obtained by solving the equation

e lm-m (3.1.8) "(T)-p 2 p

Equation (3.1.8) has a unique solution since i > 0. Let n be the unit normal vector in a positive

direction to a cell surface in a finite-volume grid, or a coordinate surface in a finite-difference

grid. If v, is the normal component of the velocity of a time-varying surface, and u, = u • n, one

can define the normal relative velocity component u' = u, - v,. The set of inviscid normal flux

components F, is given by the vector

F,=[pu' mu'+pn eu'+pu,] r. (3.1.9)

Following Harten (1983a), we obtain the transformed variable W from

o. (3.1.1o) W'=

where the function r/(U) is given by

. = (3.1.11)

Here again the second element of the algebraic vector Fn is a physical vector quantity, t_(._) is

an arbitrary function of S, and ,.q can be expressed as a function of U using (3.1.5), (3.1.6) and

(3.1.8). Using (3.1.2) and (3.1.8), we can express the differential of (3.1.11) as

_ {[e- 2"_-p(1+t3)]dp-m.dm+ pde}, drl= p

(3.1.12)

where _ = d_/dg and

(3.1.13)
Ig

Note that/3 is in

multiplication as

general a function of 3. We can rewrite (3.1.12) in

P ["1de

the form of a matrix

(3.1.14)

where the vector dot product is implied in multiplying the second element of the row vector by the

second element of the column vector. In the rest of the paper, we will use the convention that in

forming the product of two matrix elements, a dot product is implied if each element is either a

physicaI vector or a tensor. From (3.1.10) we can express dr/as dr/= WTdU. We thus obtain

W=[_" w _]r=p[e-27-p(l+/3) -m pit. (3.1.15)

In order to investigate the homogeneity condition, it is useful to introduce the function

_" W'W u(w) - - = - +/3). (3.1.16)

w 2(9)2

Note that ,u involves only thermodynamic variables, and is also a homogeneous function of W of

degree 0. Since _ is in general a nonlinear function of T, the homogeneity condition can only be

satisfied if T is a homogeneous function of W of degree 0. In view of (3.1.5) and (3.1.16) this can

only be accomplished if

/3(g) = constant. (3.1.17)

We can now solve for ,_(S) from (3.1.13). The sign of _ determines the positive definite condition

on Uw; the scale of _ does not affect any numerical calculations. Anticipating the positive definite

condition which will be derived below, we find that the simplest solution of (3.1.13) is

=/3e (3.1.18)

which then gives

= -e -_/_. (3.1.19a)

Substituting (3.1.5), this can be expressed as

= _(pf)l/a. (3.1.19b)

The expression for p(W) obtained using (3.1.4), (3.1.7), (3.1.15) and (3.1.19b) is

TO

p- f(_)(-w) _. (3.1.20)

Therefore p is a homogeneous function of W of degree/3. From (3. l. 15) and the definition of m

it follows that u, u,_ and u' are all homogeneous functions of W of degree 0. We conclude from
2O

definitions (3.1.7) and (3.1.9) that the vectors U and N,., are both homogeneous functions of W of

degree/3.

While we cannot obtain an explicit formula for U(W), we can derive an explicit expression for

the symmetric matrix Uw as functions of U. The upper triangular part can be written as

1 Iap apu ae + bp

Uw : _ L apuu - pI u[ae + (b - 1)p] , (3.1.21a) "" --_- -- U • u) + + 3)

where

1

"(T'_) = _ + 1 +/3 1, (3.1.21b)

= (1 -e1 + ,8, (3.1.21c)

and I is the identity tensor.

Let A = Uw and u = qe, where e is a unit vector in the direction of u and q is the magnitude

of u. The positive definite condition for A is determined by the quadratic form XTAX, where the

vector X has the general form X = [ zl z2e z3 IT, and the dot product is implied in a product

involving the unit vector e. Since e • ee • e = e • I • e = e • e = 1, the quadratic form can be

written in terms of ordinary scalar quantities as x'TA'x ', where X' = [ zl z2 z_ ]T, and the

upper triangular part of A ' becomes

A, = -_i lap apq'apq- p q[ae ae+ (b+ bp- 1)p] ] . (3.1.22a)

-7- ae 2 + lbe 1 )+ .

Therefore the positive definite condition for Uw is obtained by calculating the signs of the leading

principal minors of (3.1.22a). Since the elementary operation of subtracting a multiple of one row

from another leaves a determinant unchanged, we can reduce (3.1.22a) by a series of elementary

operations to the matrix

- -p - . (3.1.22b)

0

This matrix is positive definite if the leading principal minors are positive. Since p > 0 and p > O,

we obtain the conditions

a/_ > 0, (3.1.23a)

-a/_ 2 > 0, (3.1.23b)

and

bjS/_' > O. (3.1.23¢)

From (3.1.23a) and (3.1.23b) it follows that _ < O, a condition already satisfied by (3.1.19a). It

then follows from (3.1.23c) that b_ < O, which can be reduced to

1 -1

> 1 + t" (3.1.24)
21

Condition (3.1.24) is satisfied if j9 > 0 or,8 < -(1 + _). It is easy to show that a value of,8 in

these two ranges satisfies a < 0, as required by (3.1.23a). Since _ > 0, the maximum value of

occurs at T,,,,.. Therefore, for fl < 0,

< -[i+ (3.1.25)

A sufficiency condition, independent of the flow problem, is obtained by replacing _(T,.,,,.) by

We can also derive an explicit expression for the symmetric flux Jacobian with respect to the

transformed variables, (F,)w, as functions of U. The upper triangular part can be written as

u'ap u' apu - pn

u'(apuu - pI) - _un + nu)

,,'(ae + br,)- ,,,,p"

a23

(113

(3.1.26a)

where

and

= +(b- - - -Pp(e + p)n,

a3z u'[_- p(2be u. u)+---_(1 + ;9)]-2ur, P(e+p). : + P

(3.1.26b)

(3.1.26c)

3.2. Roe Riemann Solver for a Thermally Perfect Gas

The extension of Roe's approximate Riemann solver to a thermally perfect gas has been given

by Abgrall (1990) and also Spekreijse and Hagmeijer (1990). (They actually considered the more

general case of a mixture of thermally perfect gases, valid for non-equilibrium flow.) We present

the results here for arbitrary three-dimensional grids, using our compact vector notation. The

Riemann solver is based on properties of the ordinary flux Jacobian matrix A = _ From (3.1,7) 8U'

and (3.1,9) it follows that we need the pressure differential

dp = x dp + ,_ arg, (3.2.1)

where

^ 1

,c(T) = = and x(T) = f- _e. (3.2.2)

The matrix A can then be written as

A

--v. n 0

K,n - u.u un - snu + u'I sn

(K, - H)u. Hn - _u.u u' + _u.

, (3.2.3)

1

where K, = _,_u. u + X, H = h + _-u- u is the total enthalpy per unit mass, and h = e + T is

the specific enthalpy. The three distinct eigenvalues of A are

A' = u' , A _ = u' +c, and A 3 =u ' -c, (3.2.4)
22

where the speed of sound c is given by

c2 --- X + _h. (3.2.5)

The multiple eigenvalue A x is associated with those conservative variables whose flux is

purely convective. These are p and the tangential component of m. In order to construct the

corresponding linearly independent eigenvectors associated with the multiple eigenvalue, we span

the plane normal to n by an arbitrary set of two basis vectors bi and the set of reciprocal

basis vectors b j, satisfying hi" b j = 6_/, where 6{ is the Kronecker delta. It follows that

bi • n = b "/• n = 0. The right eigenvector matrix R can then be written as

R = cbi u + en u- en

K2 cbi • u H + cun H -cun

, (3.2.6)

where K2 = _u. x u -- X/t¢.

Among the various approximate Riemann solvers, the most common one uses the Roe average

because of its simplicity and its ability to satisfy the jump conditions across discontinuities exactly.

In those solvers based on local linearization, the flux at a point separating two states ULand U R is

based on the eigenvalues and eigenvectors of some average A. The optimum choice for A is one

satisfying

AF,, = a _V, (3.2.7")

where A(.) = (.)R _ (.)z. This choice of Acaptures discontinuities exactly. One way of obtaining

is to seek an average state U, which is a function of U z and U R, such that

A = A(U). (3.2.8)

Such a state is known as a Roe-averaged state. One can easily show that

U = au L + (1 - a)u a (3.2.9)

and

where

= ,_HL+ (1 - ,_)HR, (3.2.10)

1

a = (3.2.11)

1 + v,Tr/p _

From the definition of H one then obtains

_=g- _u.u. (3.2.12)

The discrete form of (3.2.1) yields the relation

Ap + _ A_" = Ap. (3.2.13)
23

The Roe-averaged sound speed is given by (3.2.5) as

_2 _- _ + _. (3.2.14)

Equation (3.2.13) provides only one relation to determine _ and _. Since X and _ are functions of

only, the simplest assumption is that _ and _ depend only on Tnand _z. Eliminating p, using

(3.1.4), we rewrite (3.2.13) as

(3.2.15)

Equating the coefficients of,o R and pL on the two sides of (3.2.15) we obtain the relations

= -- (3.2.16a) Ae

and

_ = . (3.2.16b)

Ae

Equations (3.2.16a,b) are replaced by (3.2.2) when AE ---, 0. An important quantity in the

approximate Riemann solver is the column vector R-XAU ". Its components are the jumps in the

characteristic variables. It is given by

where

R-xAU =

Ap- Ap/ 

_bJ. Au/_

 (zxp/ ÷ Au/ )

_ Au/ )

(3.2.17)

= (3.2.18)

3.3. Entropy Splitting for Other Equation Sets

In this subsection we examine the possibility of applying entropy splitting to other equation

sets. We first consider non-equilibrium flow, consisting of a mixture of different species, each

obeying a thermally perfect gas law. The motivation is again the fact that this case can also be

treated exactly with Steger-Warming flux-vector splitting and Roe's approximate Riemann solver.

We next consider the equations of magnetohydrodynamics, since there is much current interest in

their solution. Finally, we investigate the artificial compressibility method applied to the solution

of the incompressible equations.

3.3.1. Non-equilibrium Flow

In non-equilibrium flow, we consider a mixture of species, each obeying a thermally perfect

gas law. The conservation law now takes the form

Ut + V. F = S, (3.3.1)
24

where S is a vector consisting of source terms for each species. The equation of state for species i

is

pi = piRiT, (3.3.2)

where pi, pi, and R i are the pressure, density, and gas constant for species i, respectively, and T

is the temperature of the mixture, assumed to be in thermal equilibrium. If there are ns species in

the mixture, the index i takes on values from

then related to pl and T by

dS' =

l to ns. The entropy fi'i and internal energy e i are

dei dPi (3.3.3)

RiT pi '

where we have introduced the normalized variable

,,,. S i

s'= (3.3.4)

Note that we have not introduced a normalized T, since the R _ are different for each species. It

again follows from (3.3.3) that ei = el(T) only. All real species satisfy the conditions _i > 0 and

_ > 0, where _i = de_/dT. Equation (3.3.3) can be integrated to obtain

p_fl = e -g', (3.3.5)

where

The arbitrary constant in the integral of (3.3.6) can be absorbed in the definition of ,_i

We can now define the density of the mixture, p, as

(3.3.6)

,,= (3.3.7)

i

the pressure p as

where

and the entropy per unit volume, ,oS, as

p = Ep'= pRT, (3.3.8)

i

,oR-- _'_.,piRi, (3.3.9)

pS -- E pi Si" (3.3.10a)

Using (3.3.4), the last equation can be written as

peg= _" pi Ri si,

i

(3.3.10b)
25

where

A

S= S. (3.3.11)

R

Note that since R is no longer a constant, S is not proportional to S. Finally, we define the internal

energy of the mixture per unit volume as

(3.3.12)

The set of conservative variables U can be represented by the vector

U [R m e] T

where the vector _ is defined as

(3.3.13a)

R=(p' p' ... p"'). (3.3.13b)

The temperature T(U) is obtained by solving implicitly the equation

Z Im. m (3.3.14) piei(T) = e 2 p '

i

where p is given by (3.3.7). Equation (3.3.14) has a unique solution since pi > 0 and _i > 0. The

set of inviscid normal flux components F,, is given by the vector

F,_ = [Ru' mu' + pn eu' + pu,_ IT. (3.3.15)

The procedure to obtain the transformed variable W follows that of subsection 3.1. Equations (3.1.I0)-(3.1.13), and (3.1.17)-(3.1.19a) in that subsection are still valid. With the aid of

(3.3.3), (3.3.7), (3.3.9), (3.3.10b), and (3.3.14) we obtain

W=[W w _]r, (3.3.16a)

where

W = (w' w' ... w "°), (3.3.16b)

^. lm.m

w i=--R-_ [-RiT(I+S-S')-_RT-e I+2 p-i ]

i = 1, ..., ns, (3.3.16c)

_rn (3.3.16d)

w - pRT'
2O

and = _ (3.3.16¢)

w - RT"

It is again useful to define a set of functions of W that are homogeneous of degree 0, involving

only thermodynamic variables, as

wi w. w _RIT(1 + g_ ._i) _/3RT_ ei i = 1,...,ns. (3.3.17)

_'(w)- _ 2(_) _ -

Equations (3.3.16e) and (3. I. 19a) can be combined to yield

(-w )a (RT) _3 = e -'2. (3.3.18)

In order to prove homogeneity, we let

p_= ai(-_)a. (3.3.19)

Substituting (3.3.5), (3.3.18) and (3.3.19) into (3.3.17), we obtain

RIT + -_- + 1 - t31n(RT) + ln(aif ') = O. (3.3.20)

Combining (3.3.4), (3.3.9), and (3.3.10b) to eliminate p', yields

R - Ei aiai

_i ai (3.3.21)

Substituting (3.3.5), (3.3.18), and (3.3.19) into (3.3.10b), we obtain

1_In(RT) = _-,i a' R i ln( a' f' ) (3.3.22)

EiaiR i

Equation (3.3.20), (3.3.21), and (3.3.22) comprise a set of coupled non-linear equations for R, T,

and a i as functions of #i. Since the pi are homogeneous functions of W of degree 0, it follows

that R, T, and a i are all homogeneous functions of W of degree 0. It is then easy to show that U

and F,., are homogeneous functions of W of degree/3.

In order to obtain an expression for Uw, we can combine the differentials of (3.3.5), (3.3.17) and

(3.3.18) to express dp i as a linear combination ofdw i, dw, d-_, dR, and dT. From the differentials_

of (3.3.4), (3.3.5), (3.3.10b),and (3.3.18) we find that dT equals a linear combination of d-_, dR,

and all the dp i. By eliminating d o from the differentials of (3.3.4) and (3.3.5) we see that dR is also

equal to a linear combination of all the dp i. Therefore, obtaining Uw requires inverting a dense

linear system. It would thus be difficult to establish the positive definite condition. Therefore the

extension of the method to non-equilibrium flow is not practically feasible. If the homogeneity

condition is not required, then one can use symmetry variables based on the physical entropy, as

was shown by Chalot et al. (1990).
27

3.3.2. Magnetohydrodynamics

The set of conservative variables U for magnetohydrodynamics is represented by the vector

U:[p m e B] r, (3.3.23a)

where B is the magnetic field vector, and

lm.nl 1

B.B. (3.3.23b)

The set of inviscid normal flux components F,_ is given by the vector

I PUI

pu'u + n(p+ _B-B) - B.B (3.3.24) F,, = ,,,, +(p+ iB-' B),,. - B.(,,. B '

u.B - B.u

where B,, = B • n.

1

Equation (3.3.23a) shows that B is an element of U, while (3.3.23b) shows that [B • B is a

term in e, which is also an element of U. It follows that both B and e cannot have the same degree

of homogeneity, and therefore U cannot be a homogeneous function. This result is not surprising,

since the magnetic field vector is not a true conservative variable, and is not expected to behave

the same way as the physical conservative variables.

3.3.3. Artificial Compressibility Method for Incompressible Flow

For the artificial compressibility method for incompressible flows, the set of conservative

variables U for this case is given by

U = [p u] r, (3.3.25)

while the set of inviscid normal flux components Fn is given by

Fn = [o'u' u'u + pnl r, (3.3.26)

where o- is the artificial compressibility. Since the second element of U is u, while one of the

terms in the second element of F,_ is uu. n it follows that U and Fn cannot have the same

degree of homogeneity. This is also not surprising, since p is not a true conservative variable.

In the compressible case, the velocity has homogeneity of degree 0, and it is the density that is

homogeneous of degree _. When the density is no longer a variable, the homogeneity property

disappears.

IV. Numerical Examples

The numerical experiments will be limited to a perfect gas. Three test cases are considered.

The first is inviscid and the last two are compressible mixing layer computations. These test cases
28

were also considered in Yee et al. (1999). The three test cases are: (l) a horizontally convecting

vortex, (2) a vortex pairing in a time-developing mixing layer with shock waves formed around the

vortices, and (3) a shock wave impinging on a spatially evolving mixing layer where the evolving

vortices must pass through a shock wave, which in turn is deformed by the vortex passage. For

the two mixing layer computations, the study wilt be limited to the choice of the/3 that determines

the amount of splitting in obtaining the same shock location as the un-split approach. For the

Navier-Stokes computations involving entropy splitting, the splitting is applied to the inviscid flux

terms, and the symmetric form of the viscous flux is not used (see Section I)

In all of the computations the classical fourth-order Runge-Kutta time discretization is employed.

For the purposes of this paper we concentrate on the non-compact central schemes (2.3.l) and

(2.3.2) with the same order of accuracy and type of base scheme for the convection and viscous

terms (if viscosities are present). Compact symmetric schemes (2.3.3) are also applicable,

but require nearly twice the CPU time that the non-compact central schemes require for 2-D

compressible mixing layer computations and will not be addressed in this paper.

If numerical dissipation is added, the filters (2.4.2), (2.4.5) and (2.4.6) are used at the end of

the full Runge-Kutta time step. Roe's (1981) average states are used in (2.4.1). For most of the

computations, the Harten and Yee (see Yee and Harten (1985a) and Yee (1989)) second-order

upwind TVD dissipation for _bt i in (2.4.2) is used. These will be notated as ACM with the

J+_

following numbers indicating the order of the base scheme for the convection and viscous terms. For

example, ACM44 means the use of fourth-order central as the base scheme for both the convection

and viscous terms. In order not to introduce additional notation, inviscid flow simulations are

designated by the same notation, with the viscous terms not activated. Computations using

symmetric TVD dissipation (Yee, 1985, limiter (2.7b)) are indicated by adding the letter "S", as in

ACM44S. Computations using entropy splitting are indicated by adding the letters "ENT" at the

end as in ACM44-ENT. To examine the performance of the entropy splitting schemes where shock

waves are absent, the computations also employ only the non-dissipative central base schemes

(without the ACM filters) designated as CEN22, CEN44 and CEN66 for second, fourth and sixth

order respectively.

The inviscid case uses a uniform Cartesian grid. The two compressible mixing layer test

cases use a uniform Cartesian grid in the z-direction and a mildly stretched Cartesian grid in the

y-direction. [n order to assess the true performance of the algorithm, no attempt is made to enhance

the resolution using appropriate adaptive grid procedures. The code used for the Yee et al.(1999)

study is employed for the present study. For non-periodic boundary conditions (BCs), the code

reduces to lower order central base schemes near the boundary points. For the current study, we

employ the same numerical BCs treatment in order to have a one-to-one comparison with the

results obtained in Yee el al. (1999). Appropriate stable boundary difference operators developed

by Strand (1994) should be used but are not yet implemented for the present study. The global

accuracy of the scheme related to intermediate BC treatment for the multi-stage Runge-Kutta

method (Carpenter et al. (1995)) is not addressed here. Except for the vortex convection problem,

all computations impose intermediate BC updates. Special treatment of the temporally (Carpenter

et al.) and spatially dependent physical BCs related to the Runge-Kutta method is not considered
29

for the two compressible mixing layer cases and such treatment is beyond the scope of this

paper. Nonreflecting BCs or characteristic inflow and outflow boundary treatments are also not

implemented. As indicated in the objective section, we explore the possible side benefits of

the entropy splitting without considering the accompanying stable numerical boundary difference

operator as a complete package for stability requirements. The three numerical examples were

chosen to consist of periodic BCs, or computational domains whose boundaries are far enough

away so as to not affect the mainstream flow activities. For the non-periodic cases, lower order

non-characteristic boundary schemes are used. Evaluation of the performance of these schemes for

the two compressible mixing layer test cases should take the above assumption into consideration.

4.1. Isentropic Vortex Evolution

The first test case is the evolution of a 2-D inviscid isentropic vortex in a free stream with

periodic BCs in both spatial directions. The free-stream flow velocity, u,,. and v_, pressure,

p_, and density, poo are (uoo,voo) = (1,0) and p,= = p_. =1. An isentropic vortex with no

perturbation in entropy (6,5' = 0) is added to the free-stream flow field as initial conditions. The

perturbation values are given by

^

= 1_1,,..! (4.1.1)

where _ is the vortex strength and 7 - 1.4. Note that the vortex strength _ should not be confused

with the/3 in section 2.1.2. Here T = _, T_ = 1.0, (_', _) = (_ - _**, y - Yvo ), where z,,. and

Y,,o are the initial coordinates of the center of the vortex, and r 2 = _-2 + if2. The entire flow field

is required to be isentropic. Thus, for a perfect gas, _ = 1.

From the relations, u = uoo + 6u, v = v,,. + 6v, T = Too + 6T, and the isentropic relation, the

resulting initial state for the conservative variables is given by

[ ]

p =T _'--_ = (Too + 6T) _'-_-_= 1 (7 - 1)/32e,-.'

87r

=o(-= + 6-) : p i-2 j

t_p_ I

p =p'r

p 1

=-C=--f-i + _P(_'_ + _)"

1

"I'-I

(4.1.3)

(4.1.5)

(4.i.s)
3O

The exact solution with given initial states is just a passive convection of the vortex with the

free-stream velocity and thus provides a good measure of the accuracy of the schemes for smooth

solutions of the nonlinear Euler equations. The vortex strength _ = 5 is fixed for all runs. Let k

be the order of the central scheme; then the initial vortex covers a domain 0 < z < 10 + 0.125k

and -5 _< y < 5 and its center is placed at (Zoo,Yo,) = (5,0). A uniform grid spacing of

Az = 0.125 and Ay = 10/(79 + k) is used. Although the actual grid size is 80 × 79, regardless

of the order of the scheme, the grid size including ghost cells to accommodate the periodic BCs

is (80 + k) × (79 + k). The reason for using an odd number of grid points in the ?/-direction is

due to the compressible mixing layer structure of the code to accommodate fluctuations added to

the inflow. The vortex is convected to the right by the mean flow velocity. Since there are no

shock waves or steep gradient regions for this flow, the filter is used only to stabilize the nonlinear

governing equations. For this reason, the filter coefficient _ (2.4.2), if needed, should be kept very

small. We use 0.001 _< _; < 0.07 for the computations. Due to the isentropic flow property, one

can set p" (2.1.1 lc) to be a constant, p* = 1 is used for this test case.

Density profiles at the centerline, y = 0, cutting through the center of the initial vortex are used

for comparing the various schemes. Due to the time and spatial discretization numerical errors,

the vortex, after long time integrations, can drift away from the centerline. The amount of drift

depends on the scheme, grid size and the time step. If the computed vortex drifts away from the

centerline but still preserves the vortex shape and strength, the centerline, 9' = 0, density profiles

do not convey the full information and can be misleading. We complement the comparison with

snap shots of density contours at different times up to 1 I0 spatial periods. Here, one period is

defined as the length of the periodic computational domain. The time required for one spatial

period is t = 10. In all of the computations 6 = 0.1, where 6 is defined by (2.26a) of Yee et

al. (1999). The limiter used is that given by Eqs. (2.25f) of Yee et al. (1999). Here, 6 is the

entropy satisfying parameter of Harten and Hyman (1983) for TVD schemes. As recommended by

Carpenter et al. (1995), no intermediate BC update is imposed to improve the time accuracy of the

multistage Runge-Kutta methods. The computations using intermediate BC updates do not have a

drastic effect, but tend to diverge a little earlier than the case where there are no intermediate BC

updates.

We present results for sixth-order schemes. (Comparisons with second and fourth-order results

are made in the following discussions.) Centerline distributions for CEN66, CEN66-ENT, ACM66,

ACM66-ENT and ACM66S-ENT are shown in Figs. 4.1.1 - 4.1.7. Snap shots of selected density

contours for these schemes are shown in Figs. 4.1.8 - 4.1.16.

The performance of the central schemes (with or without ACM) using entropy splitting and

their un-split cousins is evaluated based primarily on vortex preservation capability after long time

integrations of up to 130 periods (t = 1300). The discussion of numerical results is divided into

the following:

(a) Effect of the order of accuracy of the base scheme (At = 0.04):

We ran some computations with second and fourth-order accurate base schemes for At = 0.04

(not shown), and compared them with the sixth-order results. CEN22 diverged after 3.5 periods,
31

CEN44 diverged after 6.8 periods, and CEN66 diverged after 5.17 periods. With no dissipation,

the more accurate CEN66 computations had non-linear instability which caused it to diverge

slightly earlier than the CEN44 case. Even so, the CEN66 density distribution was very accurate

up to 5 periods, as seen from Fig. 4. I. lc. It did show oscillations at period 5, indicating incipient

instability at period 5.17.

With the entropy splitting, the CEN22-ENT case diverged at 4.4 periods, the CEN44-ENT at

13 periods, and the CEN66-ENT at 17 periods. The latter distribution is shown in Fig. 4.1.2c.

Note that the entropy splitting allowed the calculation to proceed many periods further before it

diverged, but, more significantly, the stability was improved as the order of accuracy was increased.

This demonstrates the stabilizing effect of the entropy splitting. For a smaller At = 0.01, the

CEN66-ENT is stable up to 53 periods. See the later discussions.

The addition of dissipation improved the performance markedly. Using the same At = 0.04,

the ACM22 case with _; = 0.07 diverged after 10 periods, the ACM44 with _ = 0.06 became

very distorted after 40 periods (with _ = 0.04 diverged after 32 periods), and the ACM66 with

= 0.06 remained stable for as long as we ran. (We stopped computing after 120 periods.) The

density contours for this last case are shown in Fig. 4.1.12. Note that the vortex center starts to

drift vertically and horizontally, and undergoes gradual smearing and distortion starting at period

40 (see later discussion). The centerline distribution in Fig. 4.1.3c is therefore shown only up to

period 30. The agreement with the exact solution is excellent.

The addition of the entropy splitting yields further improvement. The ACM22-ENT with

= 0.07 remained stable, but the solution became very distorted at 30 periods. The ACM44-ENT

with _ = 0.04 similarly became severely distorted, and drifted vertically and horizontally beyond

40 periods. The ACM66-ENT solution with _; = 0.04 remained very good up to period 120, as

shown in Fig. 4.1.13, although it started to drift to the right and upward even at period 30. This

drift is evident from the centerline distributions in Fig. 4.1.4d.

(b) Effect of ACM dissipation:

We first will compare sixth-order computational results without dissipation (CEN66) with those

due to added second-order upwind dissipation (ACM66). Results without entropy splitting for

At = 0.04, 0.02, and 0.01 are shown in Figs. 4.1. I, 4.1.8 and 4.1.9 for CEN66, and Figs. 4. 1.3 and

4.1.12 for ACM66. The CEN66 computations diverged shortly after five periods for all three time

steps. The ACM66 computations remained stable for as long as we ran (120 periods) for all three

time steps. Due to the drift of the vortex center, the centerline distributions in Fig. 4.1.3 are shown

only up to 30 periods. Note that up to that time, the larger time step, At = 0.04 experienced less

drift than the smaller time steps. The density contours for At = 0.04 show smearing and distortion

as we increase the duration of the time integration.

Results with entropy splitting are shown in Figs. 4.1.2, 4. l. 10 and 4.1.11 for CEN66-ENT, and

Figs. 4.1.4, 4.1.6, and 4.1.13-4.1.15 for ACM66-ENT. The entropy splitting for the non-dissipative

computation (CEN66-ENT) allows a much longer time integration than that of CEN66 before it

diverged, but even for the smallest time step, At = 0.01, the solution for CEN66-ENT deteriorated
32

after 30 periods. Since ACM66-ENT allows a stable computation with less dissipation, the

entropy splitting (_ = .01 and At = 0.02, 0.01) eliminated the distortion and smearing found in

the ACM66 computations. The center of the vortex still drifted as we increased the time of the

computation, but for the smallest time step, At = 0.01, there was only a very small drift to the

right, even at 130 periods. Otherwise the vortex remained undistorted. Figures 4.1.5 and 4.1.6

show the centerline density distribution for 10 - 130 periods with a 10 period increment. One can

see the drifting effect as a function of the At. Aside from the drifting, the vortex is still quite

accurate. This is evident from the density contours Figs. 4.1.14 and 4.1.15.

We would like to point out that the vertical drifting of the vortex away from the centerline

y = 0 and horizontal drifting (or rather shifting) are quite common for all schemes beyond 30

periods. Depending on the scheme, the amount of numerical dissipation and the time step, drifting

can occur as early as 5 periods. We believe that the vertical drifting is due largely to the spatial

numerical dissipation of the scheme. This is evident from the n refinement study. See Figs. 4.1.13

and 4.1.15. The horizontal drifting is due largely to the phase error of the time integrator. This

is evident from the time step refinement study on ACM66-ENT using n = 0.01. See Figs. 4.1.5,

4.1.6, 4.1.14 and 4.1.15.

(c) Effect of the adjustable ACM parameter n:

Although we experimented with various values of _; (0.01 __<_ < 0.07) to find the optimum

value, we show some results only for the case giving the best solution, namely ACM66-ENT.

For the larger time step, At = 0.04, Fig. 4.1.4 shows a negligible effect up to 30 periods when

increasing the value of _ from 0.01 (Fig. 4.1.4c) to 0.04 (Fig. 4.1.4d). The computation _ = 0.04

and At = 0.04 remained stable, although with some distortion and smearing, as shown by the

density contours in Fig. 4.1.13. Actually, this case gives slightly better results than the one using

= 0.01 and the same time step size (not shown). The computation is not stable for _; < 0.01 and

At > 0.01 for the ACM66-ENT scheme.

(d) Effect of entropy splitting:

The cases discussed above clearly show the advantages of using entropy splitting. For the cases

without dissipation, CEN66-ENT, the splitting allowed the computation to proceed for a much

longer time before it became unstable. The real advantage came when used in conjunction with

the upwind filter, ACM66-ENT, where excellent solutions with just a very small amount of drift

were obtained after long time integrations of 130 periods.

(e) Effect of the time step:

In general, decreasing the time step At gave a better solution. For the base scheme CEN66,

the improvement in decreasing the time step from 0.04 to 0.01 was negligible. However, with

entropy splitting, decreasing the time step had a significant effect, as seen from the results for

CEN66-ENT in Figs. 4.1.2, 4.1.ll and 4.1.10. The improvement with decreasing time step was

not as marked when the ACM filter was applied. ( Fig. 4.1.3 showing the density distributions for

ACM66 actually appears to show an improvement with increasing time step. This is misleading,
33

since the slight upward drift of the vortex center produced a shift in the centerline distributions for

the smaller time steps. The density contours illustrate similar vortex preservation capability (not

shown)) The improvement in decreasing the time step from 0.02 to 0.01 is clearly shown in the

density distribution and density contours for ACM66-ENT of Figs. 4.1.5, 4.1.6, 4.1.14 and 4. l. 15.

The significant downward drift for At = 0.02 has been totally eliminated for At = 0.01. The

drift to the right has also been virtually eliminated. A quantitative evaluation of the solution for

At = 0.01 can be obtained from the centerline distribution for ACM66-ENT in Fig. 4.1.6. Even

at 130 periods, the profile is undistorted, with a shift due to the small drift to the right.

(f) Effect of symmetric vs. upwind ACM:

We confine the comparison to the best solution, namely ACM66-ENT for At = 0.01. The

centerline distributions up to period 30 for the symmetric ACM case, ACM66S-ENT in Fig. 4.1.7

are as good as for the upwind case, ACM66-ENT in Fig. 4.1.4a. The symmetric TVD dissipation of

Yee (1985), limiter 2.7b was used. Note that the value of _ in the former case has been decreased

to 0.005. A more meaningful comparison can be made by examining the density contours up to

110 periods for the two cases in Figs. 4.1.16 and 4.1.15. Note that the symmetric ACM solution

undergoes some distortion and upward drift as the period is increased. Actually ACM66S-ENT

with _ = 0.005 and At = 0.01 produce better results than ACM66-ENT with _ = 0.04 and

At = 0.01 or _ = 0.01 and At = 0.04 The drifting behavior also occurs with ACM66-ENT if

>__0.04 or for _ = 0.01 and At = 0.04 If we had used a larger value of _ = 0.01 for the

symmetric case, the solution would have become very inaccurate. On the other hand, a smaller

value of _ = .001 produced an unstable solution due to insufficient dissipation. It appears that

the use of symmetric ACM filter in conjunction with entropy splitting is also computationally

attractive. In all of the computations using ACM, the solution is quite sensitive to the value of

and the time step size, although, the upwind ACM appears to be a bit less sensitive.

(g) Effect of the splitting parameter/3:

All of the calculations shown for the entropy splitting have been for a value of/3 = 1

(a = -I.8). This produces an equal amount of conservative and non-conservative splitting.

We have also run some cases for/3 = 0.5 (a = -1.6). This gives a splitting that is one-third

conservative and two-thirds non-conservative. The results are slightly worse than for the/3 = 1

case. Increasing the conservative proportion beyond 80% will defeat the purpose of using the

splitting for this particular example since the gain in stability is diminished by the expense of the

added CPU computation required by the splitting.

In summary, the use of entropy splitting in conjunction with an upwind TVD ACM filter has

preserved a horizontally convecting vortex with great accuracy after long time integration of 130

periods. The splitting helps minimize the use of numerical dissipation. To the authors' knowledge,

highly accurate finite discretization computations previously reported in the literature were only

carried out up to 10 periods of integration.
34

4.2. Vortex Pairing in a Time-Developing Mixing Layer

This test case involved vortex growth and pairing in a temporal mixing layer at a convective

Mach number equal to 0.8. At this Mach number there are shock waves (shocklets) that form

around the vortices and the problem is to compute accurately the vortex evolution while avoiding

oscillations around the shocks. Previous calculations of the problem can be found in Sandbam

and Reynolds (1989), Lumpp (1996), Fu and Ma (1997), Sandham and Yee (1998) and Yee et al.

(1999). Figure 4.2.1 shows a schematic of the physical problem. Here we set up a base flow as in

Sandham and Yee (1989)

u = 0.5 tanh(2y), (4.2.1)

with velocities normalized by the velocity jump ua - u2 across the shear layer and distances

normalized by the vorticity thickness,

ul - uz (4.2.2)

= (au/ay),,,.."

Subscripts i and 2 refer to the upper (y > 0) and lower (y < 0) streams of fluid respectively. The

normalized temperature and hence local sound speed squared is determined from an assumption of

constant stagnation enthalpy

c'= cl + ,,,). (4.2.3)

Equal pressure through the mixing layer is assumed. Therefore, for this configuration of u2 = -ul

both fluid streams have the same density and temperature for y ---, +e_. The Reynolds number

defined by the velocity jump, vorticity thickness and kinematic viscosity at the free-stream

temperature is set equal to 1000. The Prandtl number is set to 0.72, the ratio of specific heats

is taken as 3' = 1.4, and Sutherland's law with reference temperature TR = 300"K is used for

the viscosity variation with temperature. The reference sound speed squared, c_, is taken as the

average of c2 over the two free streams.

Disturbances are added to the velocity components in the form of simple waves. For the normal

component of velocity we have the perturbation

2

v' = _ ap, cos(2a'kz/L_, + q_)exp(-y2/b), (4.2.4)

/o=1

where L, = 30 is the box length in the z-direction and b = 10 is the y-modulation. In our test

case we simulate pairing in the center of the computational box, by choosing the initially most

unstable wave k = 2 to have amplitude a2 = 0.05 and phase ,;b2 = -a-/2, and the subharmonic

wave k = 1 with ax = 0.01 and _bl = -a'/2. The u-velocity perturbations are found by assuming

that the total perturbation is divergence free. Even though these fluctuations correspond only
35

approximately to eigenfunctions of the linear stability problem for a compressible mixing layer,

they serve the purpose of initiating the instability of the mixing layer and have the advantage as a

test case in that they can be easily coded.

Numerically the grid is equally spaced and periodic in the z-direction and stretched in the

y-direction, using the mapping

Ly sinh(byrl) (4.2.5)

Y= Y si (by) '

where we take the box size in the v-direction Ly = 100, and the stretching factor by = 3.4. The

mapped coordinate r/ is equally spaced and runs from -1 to +1. The boundaries at +Ly/2 are

taken to be slip walls. For example, at the lower boundary

Pl : P2

=

= o,

(e), = [4(e)2-

(4.2.6a)

(4.2.65)

(4.2.6¢)

(4.2.6d)

where subscripts here refer to the grid point and e is the total energy.

We compute this test case on a 101 × 101 grid. A gird refinement study was performed in Yee

et al. (1999) Figure 4.2.2, a reference solution taken from Yee et al. (1999) using ACM44, shows

a snapshot of the temperature contours at t = 40, 80, 120 and 160 using ACM44, illustrating the

roll-up of the primary vortices followed by vortex merging. Shock waves and shears form around

the vortices with a peak Mach number ahead of the vortex of approximately 1.55 at t = 120. The

grid is 201 × 201.

For this vortex pairing in a time-developing mixing layer, we study only the effect of the choice

of the arbitrary splitting parameter ,8 (i.e., the proportion of conservative and non-conservative

parts of the splitting) in obtaining the same shock location as the un-split approach with scheme

ACM66-ENT using At = 0.1. In all of the computations for the vortex pairing case, limiter

(2.25h), and 6 = 0.25 (2.26a) of Yee et ai. (1999) are used and _ = 0.7 (2.4.2) is used for the

nonlinear fields for the ACM methods. Intermediate BC update is imposed in order to have a

one-to-one comparison with the Yee et al. (1999) results.

We consider a = -100,-10,-5,-3,-2,-1.8 ('3 = 1), -1.6, 0.1, 1 (,3 = -6),

2, 5, 10, 100 with '3 = (a + 3,)/(1 - 7)- The scheme diverges for a = 0.1. This corresponds to 136.36% of the conservative proportion and -36.36% the non-conservative proportion.

By monitoring the left A-shock location, studies indicate that for a > 151, the same shock

location and shock strength of the )_-shock are obtained as in the un-split approach. With the

exception of a small increase in spurious noise in the vicinity of the shock (not shown), it is

surprising to see that a slightly over 100% conservative proportion (a > 0) and the corresponding
36

negative non-conservative proportion would give the correct shock strength and location. Away

from the A-shock area, the solution is less sensitive to the choice of a > [5].

For the physical choice of a < -% we obtain the opposite effect as compared to that of a > 0,

in terms of spurious noise. As a + 1 - 27, the entropy splitting has a spurious noise reduction

capability when compared with the un-split approach. For example, the shock strength and location

that are a bit away from the A-shock are almost the same for a = -3 (/3 = 4) as for a = -5

(/3 = 9), except that a = -3 has a bigger smoothing effect on the spurious noise generated

by the scheme (especially when a more compressive flux limiter is employed). In addition, for

-5 < a < -7, a bigger negative a in that range results in more shift of the A-shock location away

from the un-sp[it approach location. For example, for a = -3, there is a shift of approximately

1 - 1 1/2 grid points. The shock strength reduction at the A-shock location is very small. A

reduction in At might improve the accuracy of the shift in the shock location, based on the vortex

convection case. Here, we only compare the results with the un-split approach using the same time

step and BCs as reported in Yee et al. (1999). Figure 4.2.3 illustrates the normalized temperature

at t = 160 using a = -5 and a = -3 on a 101 x 101 grid. Except for a slight noise reduction in

the vicinity of the shock, the a = -5 solution is almost identical to the un-split computation. To

illustrate the noise reduction capability of the splitting, Fig. 4.2.4 shows a comparison of the split

and un-split forms with the ACM filter turned off for the linearly degenerate fields (u, u and v, v

characteristic fields) using a = -3 (/3 = 4). One can see the noise reduction effect of the entropy

splitting on the scheme which, at the same time, maintains the accuracy of the shock and shears

away from the A-shock location as in the un-split approach. Since the entropy splitting requires the

same amount of filter as the un-split approach for this type of rapidly developing shock-turbulence

interactions, its stabilizing effect is only on the spurious noise reduction, and the benefit is not

as pronounced as for the smooth flow case. However, for turbulent flows involving long time

integrations that contain weak or no shock waves, the entropy splitting could help minimize the

use of numerical dissipation due to its unique nonlinear stability property. A separate investigation

is in progress.

4.3. Shock Wave Impingement on a Spatially-Evolving Mixing Layer

The third test case has been developed to test the behavior of the schemes for shock waves

interacting with shear layers where the vortices arising from shear layer instability are forced to

pass through a shock wave. Figure 4.3.1 shows the schematic of the physical problem. An oblique

shock is made to impact on a spatially-developing mixing layer at an initial convective Mach

number of 0.6. The shear layer vortices pass through the shock system and later through another

shock, imposed by reflection from a (slip) wall at the lower boundary. The problem has been

arranged with the Mach number at the outflow boundary everywhere supersonic so that no explicit

outflow boundary conditions are required. This allows us to focus on properties of the numerical

schemes rather than on the boundary treatment.

Figure 4.3:2, a reference solution taken from Yee et a1.(1999) using ACM44, shows the nature of

the flow on a 641 x 161 grid illustrating the pressure, density and temperature fields using the

ACM44 method and the same limiter as the pairing case with _: = 0.35 for nonlinear characteristic
37

fields and _ = 0.175 for linear characteristic fields. The time step is At = 0.12. The entropy

satisfying parameter 6 of Harten and Hyman (1983) is set to 0.25. An oblique shock originates

from the top left hand corner and this impacts on the shear layer at around z = 90. The shear layer

is deflected by the interaction. Afterwards we have a shock wave below the shear layer and an

expansion fan above it. The shock wave reflects from the lower solid wall and passes back through

the shear layer. The lower wall uses a slip condition so no viscous boundary layer forms and we

focus on the shock-wave interaction with the unstable shear layer.

The inflow is specified again with a hyperbolic tangent profile, this time as

, = 2.5+ 0.5t h(2y), (4.3.1)

giving a mixing layer with upper velocity Ux = 3, lower velocity u2 = 2, and hence a velocity

ratio of 1.5. Equal pressures and stagnation enthalpies are assumed for the two streams, with

convective Mach number, defined by

M, - ul - u2 (4.3.2)

ci + c_

where cl and c2 are the free stream sound speeds, equal to 0.6. The reference density is taken as

the average of the two free streams and a reference pressure as (pa + p2)(ua - u2)2/2. This allows

one to compute the inflow parameters as given in the first two columns of Table 4.3.1. Inflow

sound speed squared is found from the relation for constant stagnation enthalpy (4.2.3). The 0 in

Table 4.3.1 is the flow inclination angle with respect to the z-direction.

The upper boundary condition given in column 3 of Table 4.3.1 is taken from the flow properties

behind an oblique shock with angle/3 = 12 °, The table also gives the properties behind the

expansion fan (column 4) and after the oblique shock on the lower stream of fluid (column 5)

computed by standard gasdynamics methods with/3 = 38.118". In practice, the conditions in

regions 4 and 5 do not correspond exactly to the simulations due to the finite thickness of the

shear layer. The Mach number of the lower stream after this shock is approximately Ms = 1.6335

and remains supersonic through all the successive shocks and expansion fans up to the outflow

boundary. The resulting shock waves are not strong, but tests showed that they could not be

computed without using shock-capturing techniques. The lower boundary was specified with the

same slip condition used for the pairing case (Equation (4.2.6)).

The Prandtl number and ratio of specific heats were taken to be the same as for the vortex pairing

test case. The Reynolds number was chosen to be 500.

Fluctuations are added to the inflow as

2

v' = Z ah eos(27rkt/T + Ct,) exp(-y2/b),

/e=l

(4.3.3)
38

with period T = _/uc, wavelength _ = 30, convective velocity uc = 2.68 (defined by

u, = (ute2 + u2e_)/(el + c2)) and b = 10. For k = 1 we take al = 0.05 and @ = 0, and for k = 2

we take a2 = 0.05 and _b = a-/2. No perturbations are added to the u-component of velocity.

The grid is taken to be uniform in • and stretched in y according to equation (4.2.5) with bu = 1.

This stretching is much milder than for the pairing problem, as we have to resolve the shear layer

even when it deflects away from y = 0. The box lengths were taken to be/3, = 200 and/;y = 40.

The reference solution indicates that vortex cores are located by low pressure regions and the

stagnation zones between vortices by high pressure regions. The shock waves are seen to be

deformed by the passage of the vortices. Another interesting observation is the way the core of

the vortex at _ = 148 has been split into two by its passage through the reflected shock wave.

In spite of the relatively high amplitude chosen for the subharmonic inflow perturbation, there is

no pairing of vortices within the computational box. We do, however, begin to see eddy shock

waves around the vortices near the end of the computational box where the local convective Mach

number has increased to around 0.66. The oscillations seen near the upper boundary for z > 120

occur where the small Mach waves from the initial perturbations arrive at the upper boundary. The

use of characteristic boundary conditions should remove this problem. Practically, the amplitude

of oscillations is not sufficient to cause numerical instability or affect the remainder of the flow.

For this shock wave impingement on a spatially-evolving mixing layer, again, we study

only the effect of the choice of the arbitrary splitting parameter /3 in obtaining the same

shock location as the un-split approach. The study is limited to ACM66-ENT and ACM66

using the fifth limiter of equation (2.25h) of Yee et al. (1999), and the same _ value and

time step size (At = 0.12) as the reference solution. Intermediate BC updates are imposed

in order to have a one-to-one comparison with the Yee et al. (1999) result. We consider

= -I00,-I0,-5,-3,-2,-1.8 (/3 = I),-1.6, 0.i, I (/3 = -6), 2, 3, 4, 5, I0, 20, I00.

Studies indicate that fora = -2, -1.8, -1.6, 0.1, 1, 2, 3, 4, 5 the solution diverges. The rest of

the ct values for [a[ > 10 produce almost identical results as the un-split case. This example poses

a more stringent requirement on the a range than the vortex pairing case. Figure 4.3.3 compares

the un-split pressure contours with the split case for a = -t-10 at t = 120 on a 321 × 81 grid with

At = 0.12. The value of a = +10 (t3 = -28.5) corresponds to a 103.6% conservative proportion

and a small negative non-conservative proportion. The _ = -10 (13 = 21.5) corresponds to a

less than 100% conservative proportion. Note that the c_ = -10 solution produces spurious noise

reduction, while the a = +10 solution actually induces more spurious noise than the un-split

approach. This opposite effect of the spurious noise phenomena for o_ > 0 and o_ < -7 is shared

with the vortex pairing example. Again, for turbulent flows involving long time integrations that

contain weak or no shock waves, we believe the entropy splitting could help minimize the use of

numerical dissipation due to its unique nonlinear stability property.

4.4. Computational Costs

For thecompressible mixing layer computations using the fourth-order Runge-Kutta method,

the central base schemes with the ACM filter are only around 25% more expensive than the same

base schemes without ACM filter. This has been achieved by only requiring one application of
39

the ACM filter per full time step for the convection terms. For LMM time discretizations, the

central base scheme with the ACM filter is only 10% more expensive than standard second-order

TVD schemes. The entropy splitting is approximately 20% more expensive than the un-split

conservative form for the 2-D mixing layer computations in conjunction with the fourth-order

Runge-Kutta method. The extra CPU time is mainly due to the fact that, for each direction, four

entropy splittings are required. If two to three time level LMM type of time discretizations are

used, less CPU time can be realized.

Summary

Our study shows that the entropy splitting can be formally extended to a thermally perfect gas,

with the internal energy being an arbitrary function of temperature. For non-equilibrium flows

consist of a mixture of different species, each obeying a thermally perfect gas law, extension of the

splitting is problematic. While we were able to prove the symmetry and homogeneity properties,

the degree of homogeneity can only be obtained by solving a system of nonlinear equations.

In addition, to obtain the 3acobian of the transformation required inverting a non-sparse linear

system. It would therefore be difficult to establish the positive definite condition. Consequently,

the extension of the method to non-equilibrium flows is not practically feasible. If the homogeneity

condition is not required, then one can use symmetry variables based on the physical entropy, as

was shown by Chalot et al. (1990). In this case, the resulting PDEs are in pure non-conservative

form and entropy splitting is no longer possible. For magnetohydrodynamics, the magnetic field

has to be added as a "conservative" variable. But the square of the magnetic field is one of the

terms in the definition of the total energy. Thus, from dimensional arguments it is clear that one

cannot obtain the homogeneity condition. A similar situation exists for the artificial compressibility

method of treating incompressible flow.

Using the same high order central schemes, numerical experiments with a 2-D vortex convection

Euler computation consisting of periodic BCs indicate that entropy splitting is more stable than the

un-split (purely conservative) approach. With an appropriate time step, numerical dissipation is

not required for up to 30 spatial periods with nearly perfect vortex preservation as opposed to only

5 periods for its un-split cousin. For even longer time integration, although numerical dissipation

is needed to stabilize the schemes, the amount required is much less than for its un-split cousin. A

nearly perfect vortex preservation of up to 120 periods was achieved.

For the mixing layer study, in order to obtain the same )_-shock strength and shock location

as the un-split approach using the same scheme, the range of the arbitrary splitting parameter/3

has to be confined to the use of at least 90% of the conservative proportion of the flux derivative.

For problems without ,_-shocks, a wider range of/3 can be used. Only a slight advantage of

the entropy splitting over the un-split approach was observed for this type of flow physics. The

advantage is in terms of noise reduction and improved nonlinear stability. There is no reduction in

the use of the ACM filter. This might largely be due to the rapidly developing flow and the high

percentage of conservative proportion required. Unlike the vortex convection with periodic BCs

in all spatial directions, these two more complicated cases consist of rapidly developing flows.

Not all of the physical BCs are periodic. In addition, the BCs consist of spatially or temporally
40

sinusoidaldisturbances.Thus,theperformanceof theschemesforthemixing layercasesispartially

cloudedby the spatially and temporallydependentphysicalBCsthat requiredspecialtreatment

in conjunctionwith the Runge-Kuttatime integrator,andalsoby the factthat wedid not impose

the morestableandappropriateboundarydifferenceoperator.Theuseof thesymmetricform of

Harten(1983a)for the viscousterm might bea sourceof improvement.In addition,a time step

reductionwill alsobenefitthe useof entropysplitting as indicated in the vortex convection study.

Without additional study, the benefit of using the entropy splitting is inconclusive for compressible

turbulence mixing applications. Since the entropy splitting requires the same amount of filter as

the un-split approach for rapidly developing shock-turbulence interactions, its stabilizing effect

is not as pronounced as for the smooth flow case. However, for turbulent flows involving long

time integrations that contain weak or no shock waves, the entropy splitting could help minimize

the use of numerical dissipation due to its unique nonlinear stability property. More rigorous

implementation of the BCs and extensive study are needed. These are the subjects of future

research.

Overall, the three numerical examples indicate a positive side benefit of the entropy splitting.

The splitting can stabilize spurious noise generated by the non-dissipative or low dissipative spatial

discretizations which are a major cause of nonlinear instability. Modern high-resolution numerical

dissipation has been the major player in improving nonlinear instabilities for short or moderate

time integrations (unsteady). Most often, added numerical dissipation is necessary for longer

time integration at the expense of excess smearing of the flow physics without resorting to finer

grids and extremely small time steps. The use of the entropy splitting form of the flux derivative

in conjunction with high-resolution filters can minimize the use of numerical dissipation. We

believe that the use of the entropy splitting is not limited to spatial central schemes (compact or

non-compact), but to spectral and spectral-like spatial schemes as well. This and the blending of

ACM filters with other filters on the possible suppression of spurious high frequency oscillations

as discussed in section 2.6 are the subjects of our near future research.

Acknowledgment

Special thanks to Tom Coakley and Dennis .l'espersen for their critical review of the manuscript.

Thanks also to Cetin Kiris for his help in the graphics illustrations.

APPENDIX A

Entropy Splitting for a Thermally Perfect Gas in 2-D Cartesian Coordinates

In this appendix we specialize the equations derived in Sections 3.1 and 3.2 to 2-D, stationary

Cartesian coordinates. The variables in (2.1.1) are defined as

U= F= rn2/p + P G= nu (A.I)

e + p)m/p + p)n/pj

where m = pu, n = pv. u and v are the Cartesian velocity components, and the other variables
41

are as defined in Section 3.1. The temperature T(U) is obtained by solving implicitly the equation

e(T) - e 1 (m z + nZ)_. (A.2)

p 2 p_

The transformed variable W takes the form

W=[u,1 _,_ w3 w,l_'=_[e-2_-r'(1+_) -_-P_' pit (A.3)

P

where _ is given by (3.1.19b). Note that w, corresponds to the component _, w2 and ws correspond

to the components w, and w4 corresponds to the component w in Section 3.1. The upper triangular

part of Uw is now written as

ap aptt apv ae + bp ]

1 -_,' - p ,_' "["' + (b- 1)r,l

Uw = _ apv z - p v[ae + (b- 1)p] '

*,' ,4 zs, _ a2_ + b__(1 + .8)

(A.4)

where qZ = uz + v 2 and a and b are given by (3.1.21b) and (3.1.21c). The upper triangular parts

of Fw and Gw can be written as

apu

(A.5)

where

and

GW=_ 1 [apv

u [ae2 + p{2(b - 1) e - qZ} + p2 {b(1 +/3) - 2}]

where

apuv

_,(._' - r,)

apv 2 - p

'.'('w"' - r')

v(apv _ - 3p)

,,[,_.+ (b- t)r,]

.,,[a. + (b- 2).1

,,_[.. + (b- 2)r,]- _(. + r,)

944

g44= ,,["e--_-_ + p{2(b- 1)_-- q'} + _-{b(1+ _) - 9-}1.

p P

The matrix A, =Ftr can be written as

0 2 1 0 A. = K1 - u (2-,_)u -_v

[( KI -,,,, v _, - H)u I:[ n,u2 -_uv

0] 0

(1 + ,,),,

(A.6)

(A.7)

(A.8)

(A.9)
42

whereKl = _gqZ2 + x, H = h + i qz2isthetotalenthalpyperunitmass, h=_+Tisthespecific

enthalpy, and X and _; are given by (3.2.2). The three distinct eigenvalues of A,. are

A,. 1 =u, A, _ =u+c, and A, s=u-c, (A.10)

where c is given by (3.2.5). The right eigenvector matrix R,. can be written as

R,:

°

0 u+c1 u c (A.11)

K2 re H+ue H-uc.I

where K2 = iq2 t - X/_, and the first two columns correspond to A,, x

R,-IAU is given by

The column vector

R.-1AO " =

Ap- Ap/_ _

_(a_,/-e_+ _a./-_)

1⁄2(ar,/_' - -_a.l-_)

(A.12)

The Roe-averaged states are defined in Section 3.2.

The matrix A v = Gtr can be written as

I 0 0

--UI._ V

Av = K1 - v 2 -_u

I.(Kz -//') v -_uv

1

(2- _)_

0

0

(I+ ,,),,

(A.13)

The three distinct eigenvalues of Ay are

Ay 1 =v, Av 2 =v+c, and Av s=v-c. (A.14)

The right eigenvector matrix Ru can be written as

1 1j C B B

D

""v = 0 v +c v -c

Q

Ks uc H + vc H vc

(A.15)

The column vector Ru-_AU is given by

Ry-a AU =

Ap- ,,,p/_, ]

_(ar,/_' + -_Av/-_)I" (A.16)
43

If we set _ = 1/(7 - 1), and make further simplifications, (A.3 - A. 16) reduce to the perfect

gas case.

APPENDIX B

Entropy Splitting for a Thermally Perfect Gas in 3D Cartesian Coordinates

In this appendix we specialize the equations derived in Sections 3.1 and 3.2 to 3D, stationary

Cartesian coordinates. The conservation law in this case can be written in the form

Ut + E, + Fy +G, =0, (B.I)

where U_ = ._g.i_,E.oU= oEo__i_,F_= _°P and Gz = oao__7. The dependent variable U is the vector of

conservative variables, while E, F and G are the flux vector components. The variables in (B.I)

are defined as

P 13/ p mu nu

U= , E= , F= m2/p+ p , G= nv , (B.2)

L Iw J rnw rt2/p+p

e + p)I/p( e + p)rn/p (e + p)n/p

where l = pu, m = pe. n = pw. u, v and w are the Cartesian velocity components, and the other

variables are as defined in Section 3.1. The temperature T(U) is obtained by solving implicitly the

equation

,(_)= e_ l(l2+m 3+n 2) (B.3)

p 2 p2

The transformed variable W takes the form

= p[e-2_'-p(l+8) -m -pe -p_, p]_', (B.4)

where _ is given by (3.1.19b). Note that w2, ws and w4 correspond to the components of w in

Section 3.1 and w5 corresponds to _ in Section 3.1. The upper triangular part of U_ is now

written as

"ap apu ape apw ae + bp

apu 2 - p apuv apuw u[ae + (b - l)p]

I apv 2 - p apvw v[ae + (b - l)p] , (B.5)

trw = _ .,,._ - p _[°. + (b - 1)p]

,,____3 + o(2b, _q2) + bp__(1 +/3) p "" p

where q2 = u 2 + v 2 + w 2 and a and b are given by (3.1.21b) and (3.1.21c). The upper triangular

parts of Ew,Fw" and Gw, can be written as

.pu ._, _ _.._ .puw .[._ + (b- 1)p] E.. = _ .(ape _ - p) o._,_ ._[._ + (b- 2)pl , (B.6)

u(apw 2 - p) uw[ae + (b - 2)p]

e55
44

where

and

e2s = u2[ae + (b- 2)p] - P- p(e + p),

ess= .joe' + p{2(b - l)e- - q2} + P2 {b(l +/9) - 2}];

P P P

(s.7)

(B.8)

apv apuv apv 2 - p apvw

u(apv 2 - p) a_vw

1 _(apu' - p) v(._' - 3p) w(._' - p)

v[,,_+ (b- x)p]

,,v[,,e+ (b- 2)p]

_,.[ae + (b- 2)p]

, (B.9)

where

f,,= _'[..+ (b- 2)p]- ;(_+p), (B.10)

and

f55 --_[ _e2

P

.2

-- + p{2(b - 1) e- - qZ} + L{b(1 +/9) - 2}];

P P

(B.11)

and

apw apuw apvw apw 2 - p

apuvw u(aFw 2 - p)

1 '1'(""' -P) ,.,,(,:,p,,, _ p) ,,(,,p_,_ p)

_"" = _ u,(,,p='- 3p)

_[.,+ (b- 1)p]

_.[..+ (b- 2)p]

g45

g55

(B.12)

where

g,, = w'[,,_+ (b- 2)p]- e(. + p),

p

(B.13)

and ,,2

+ p{2(b - 1) e-- - q'} + r---{b(1 +/3) - 2}].

P P

(B.14)

The matrix A, = Eu can be written as

0 1 0 0 i ]

KI--' (2-,¢)u -_v -,¢w

--UI_ I; U 0

--UlD _ 0 U

, (B.15)

where K1 =itcq* 2 +x,H=h+i q* _ is the totalenthalpy per unit mass, h = e + Tis the specific

enthaipy, and X and _ are given by (3.2.2). The three distinct eigenvalues of A= are

A=,* =u, A,,' =u+c, and A, 3=u-c, (B.16)
45

where e is given by (3.2.5). The right eigenvector matrix R. can be written as

I 0 0 1 I 1

I

u 0 0 u+c u-c

v c 0 v v

W 0 C

K2 ve wc H + uc H - ucj

where K2 = ,

iq2 - X/_, and the first three columns correspond to A.*

R,-ZAU is given by

R,-IAU =

Ap_ ap/_'

_(Ap/_' + -_,_,,/_)

The Roe-averaged states are defined in Section 3.2.

The matrix Ay = Fu can be written as

0 0 1 0 0

-uv v u 0 0

A_. = Kz - v 2 -_¢u (2- _)v -_w _;

( -wv 0 w v 0 K1- H) v -_uv H-Icy 2 -_wv (l+_;)v

The three distinct eigenvalues of Av are

Av I =v, Av 2=v+c, and Ay s=v-c.

The right eigenvector matrix Ry can be written as

oo1 11 C 0 U U

R_ = 0 0 v+c v-c .

0 C tO tO

LK2 uc we H+vc H vc

The column vector Ry-*AU is given by

R.-*AU =

Ap - Ap/_2 l

-_A_,/_ I

-_Aw/-e I

_( Ar, l_ 2 - _AvI-_) J

(B._?)

The column vector

(B._8)

(B._0)

(B.20)

(B.21)

(B.22)
48

The matrix A, = Gtr can be written as

AZ 0 0 0 1

--Uto tO 0 U

--Vto 0 tO V

KI - to2 -,_u -,., (2 - ,_)to

L(Kx - H) to -_uto -_vw It - _to2

The three distinct eigenvalues of Az are

_ 1 =to, A 2=w+c, and A 3

The fight eigenvector matrix R, can be written as

RZ

1 0 0 1

B C 0 B

U 0 C 17

to 0 0 to + e

K2 uc vc 11+toe

The column vector R,-aAU is given by

R,-IAU =

(B.23)

= to-c. (B.24)

W--C

-- toC

Ap- Ap/-6 2

-_Au/_

_Av/e

_(z_p/-_2 - -pzxto/-e)

(B.25)

(B.26)

APPENDIX C

Caloric Equation for an Ideal Diatomic Gas

In this appendix we present the internal energy functions for an ideal diatomic gas. This gas

is one for which rotation is fully excited, the vibrating molecule is approximated as a harmonic

oscillator, and electronic contribution is neglected. The harmonic oscillator model implies that

the temperature is small compared to the dissociation temperature. Under these assumptions, the

internal energy is given by

5 0

_(T) = _ + ---:-_,,, (c.1)

where 0 is the characteristic temperature for vibration, and

,. : 0/9. (c.2)

The quantity _ = de/dT is needed to solve Eq. (3.1.8) iteratively, and to evaluate the functions a,

b, and X in (3.1.21b), (3.1.2 It), and (3.2.2). It is given by

5 r2e" (C.3)

e(T) : 2 + (e" - 1) 2.
The function f(T) defined in (3.1.6) is given by

[ re. ]

f(ff') = T-S/'(e" _ 1)exp - [e-V-L-__ lj

47

(C.4)

References

Abgrall, R. (1990) "Extension of Roe's Approximate Riemann Solver to Equilibrium and

Nonequilibrium Flows," Notes on Numerical Fluid Mechanics, Vol. 29, pp. 1-10.

Alpert, P. (1981), "Implicit Filtering in Conjunction with Explicit Filtering," J. Comput. Phys.,

Vol. 44, pp. 212-219.

Beam, R.M. and Warming, R.F. (1976), "An Implicit Finite-Difference Algorithm for Hyperbolic Systems in Conservation Law Form," J. Comput. Phys., Vol. 22, pp. 87-110.

Carpenter, M.H., Gottlieb, D., Abarbanel, S. and Don, W.-S. (1995), "The Theoretical Accuracy

of Runge-Kutta Time Discretizations for Initial value Problem: A Study of the Boundary Error,"

SIAM J. Sci. Comput., Vol. 16, No. 6, pp. 1241-1252.

Chalot, F., Hughes, T.J.R., and Shakib, F. (1990), "Symmetrization of Conservation Laws with

Entropy for High-Temperature Hypersonic computations," Computing Systems in Engineering,

Vol. 1, pp. 465-521.

Dahlquist, G. (1979), "Some Properties of Linear Multistep and One-Leg Methods for Ordinary

Differential Equations," Conference Proceeding, 1979 SIGNUM Meeting on Numerical ODE's,

Champaign, I11.

Fu, D. and Ma, Y. (1997), "A High Order Accurate Difference Scheme for Complex Flow

Fields," J. Comput. Phys., Vol. 134, pp. 1-15.

Gaitone, D.V. and Visbal, M.R. (1999), "Further Development of a Navier-Stokes Solution

Procedure Based on Higher-Order Formulas," AIAA Paper 99-0557, Reno, NV.

Gear, C.W. (1971), Numerical Initial Value Problems in Ordinary Differential Equations,

Prentice-Hall, Englewood Cliffs.

Gottlieb D. and Turkel, E. (1976), "Dissipative Two-Four Methods for Time Dependent

Problems", Math. Comp., Vol. 30, pp. 703-723.

Gerritsen, M. (1996), "Designing an Efficient Solution Strategy for Fluid Flows," PhD Thesis,

SCCM Program, Stanford University, Oct. 1996.

Gerritsen, M. and Oisson, P. (1996), "Designing an Efficient Solution Strategy for Fluid

Flows: I. A Stable High Order Finite Difference Scheme and Sharp Shock Resolution for the Euler

Equations," J. Comput. Phys., Vol. 129, pp. 245-262.
48

Gustafsson, B. and Olsson, O. (1995), "Fourth-Order Difference Method for Hyperbolic

IBVPs," J. Comput. Phys., Vol. 1177, pp. 300-317.

Harten, A. (1978), "The Artificial Compression Method for Computation of Shocks and

Contact Discontinuities: III. Self-Adjusting Hybrid Schemes," Math. Comp.. Vol. 32, No. 142,

pp. 363-389.

Harten, A. (1983a), "On the Symmetric Form of Systems for Conservation Laws with Entropy,"

ICASE Report 81-34, NASA Langley Research Center, 1981 also, J. Comput. Phys., Vol. 49, pp.

151-164.

Harten, A. (1983b), "A High Resolution Scheme for Computation of Weak Solutions of

Hyperbolic Conservation Laws," J. Comput. Phys., Vol. 49, pp. 35-393.

Harten, A. and Hyman, J.M. (1983), "A Self-Adjusting Grid for the Computation of Weak

Solutions of Hyperbolic Conservation Laws," J. Comput. Phys., Vol. 50, pp. 235-269.

Hughes, T.J.R., Franca, L.P. and Mallet, M. (1986), "A New Finite Element Formulation for

Computational Fluid Dynamics: I. Symmetric Forms of the Compressible Euler and Navier-Stokes

Equations and the Second Law of Thermodynamics," Computer Methods in Applied Mechanics

and Engineering, North-Holland, Vol. 54, pp. 223-234.

Kreiss, H.-O. and Scherer, G. (1977), "On the existence of Energy Estimates for Difference

Approximations for Hyperbolic Systems, Technical Report, Department of Scientific Computing,

Uppsala University.

Lele, S.A, (1992), "Compact Finite Difference Schemes with Spectral-like Resolution," J.

Comput. Phys., Vol. 103, pp. 16-42.

Lumpp, T. (1996), "Compressible Mixing Layer Computations with High-Order ENO

Schemes," 15th Intl. Conf. on Num. Meth. in Fluid Dynamics, Monterey, June 1996.

Olsson, P. (1995a), "Summation by Parts, Projections and Stability," Math. Comp., Vol. 64,

pp. 1035-1065.

Olsson, P. (1995b), "Summation by Parts, Projections and Stability II.," Math. Comp., Vol.

64, p212.

Olsson, P.

Report 95-06,

Olsson, P.

Conservation

(1995), "Summation by Parts, Projections and Stability III," RIACS Technical

NASA Ames Research Center.

and Oliger, J. (1994), "Energy and Maximum Norm Estimates for Nonlinear

Laws," RIACS Technical Report 94-01, NASA Ames Research Center.

(1981), "Approximate Riemann Solvers, Parameter Vectors, and Difference

Comput. Phys., Vol. 43, pp. 357-372.

Sandham, N.D. and Reynolds, W.C. (1989), "A Numerical Investigation of the the Compressible

Mixing Layer," Report TF-45, Department of Mechanical Engineering, Stanford University.
4g

Sandham,N.D. and Yee, H.C. (1989),"A NumericalStudyof a Classof TVD Schemesfor

CompressibleMixing Layers," NASA TM- [02194, May 1989.

Sandham,N.D. and Yee, H.C. (1998)"Performance of Low Dissipative High OrderTVD

SchemesforShock-TurbulenceInteractions,"RIACSTechnicalReport98.10,April 1998,NASA

AmesResearchCenter.

Spekreijse,R. and Hagmeijer, R. (1989) "Derivation of a Roe Schemefor an N-Species

ChemicallyReactingGasin Thermal Equilibrium,"Notes on Numerical Fluid Mechanics, Vol.

29, pp. 522-532.

Steger, J.L. and Warming, R.F. (1981), "Flux Vector Splitting of Inviscid Gas Dynamics

Equations with Applications to Finite Difference Methods," J. Comput. Phys., Vol. 40, pp.

263-293.

Strand, B. (1994), "Summation by Parts for Finite Difference Approximations for d/dz," J.

Comput. Phys., Vol. 110, pp. 47-67.

Vinokur, M. (1989) "An Analysis of Finite-Difference and Finite-Volume Formulations of

Conservation Laws," J'. Comput. Phys., Vol. 81, pp. 1-52.

Visbal, M.R., and Gaitonde, D.V. (1998), "High-Order Accurate Methods for Unsteady

Vortical Flows on Curvilinear Meshes," AIAA Paper 98-0131, Reno, NV.

Vichnevetsky, R. (1974), "Numerical filtering for Partial Differencing Equations," Numerical

Applications Memorandum, Rutgers University, NAM 156.

Yee, H.C. and Harten, A. (1985a), "Implicit TVD Schemes for Hyperbolic Conservation Laws

in Curvilinear Coordinates," AIAA Paper No. 85-1513, Proceedings of the AIAA 7th CFD

Conference, Cinn., Ohio, July, 1985, also AIAA J., Vol. 25, No. 2, 1987, pp. 266-274.

Yee, H.C. (1985b), "On Symmetric and Upwind TVD Schemes," Proceedings of the 6th

GAMM-Conference on Numerical Methods in Fluid Mechanics, Sept. 1985, Gottingen, Germany.

Yee, H.C. (1986) "Linearized Form of Implicit TVD Schemes for Multidimensional Euler

and Navier-Stokes Equations," Computers and Mathematics with Applications, Vol. 12A, pp.

413-432.

Yee. H.C. (1987) "Construction of Explicit and Implicit Symmetric TVD Schemes and Their

Applications," J. Comput. Phys., .Vol. 68, pp. 151-179; also NASA TM-86775 July 1985.

Yee. H.C. (1989), "A Class of High-Resolution Explicit and Implicit Shock-Capturing

Methods," VKI Lecture Series 1989-04 March 6-10, 1989, also NASA TM-101088, Feb. 1989.

Yee. H.C., Klopfer, G.H. and Montagne, J.-L. (1990), "High-Resolution Shock-Capturing

Schemes for Inviscid and Viscous Hypersonic Flows," J. Comput. Phys., Vol. 88, pp. 31-61.

Yee. H.C., Sandham, N.D. and Djomehri, M.J. (1999), "Low Dissipative High Order ShockCapturing Methods Using Characteristic-Based Filters," RIACS Technical Report 98.1l, May

1998, also, J. Comput, Physics, Vol. 150, pp. 199-238.
Property

u-velocity

v-velocity

8 (degrees)

density p

pressure p

sound speed c

Mach number IMI

(1)

3.0000

0.0000

0.0000

1.6374

0.3327

0.5333

5.6250

(2)

2.0000

0.0000

0.0000

0.3626

0.3327

1.1333

1.7647

(3)

2.9709

-0.1367

2.6343

2.1101

0.4754

0.5616

5.2956

(4)

2.9792

-0.1996

3.8330

1.8823

0.4051

0.5489

5.4396

(5)

1.9001

-0.1273

3.8330

0.4173

0.4051

1.1658

1.6335

Table4.3.1.Flow propertiesfortheshock-wave/shear-layer testcaseinvariousregionsofthe

flow:(I)upperstreaminflow,(2)lowerstreaminflow,(3)upperstreamafteroblique

shock,(4)upperstreamafterexpansionfan,(5)lowerstreamaftershockwave.

50
\

\

i

\

\

\m

\

\

\

.

\
!

/
II

<I

o

- !

I U

"_'0

il

I_.-. ............... (_............ :.......... " _. _.......................... E.......... _......i I _ "_

l li ili_ _| I,_ I

i ,', i ',!!

d

J

S3
r'--

.............................................. .................... _ _. _................ ._.............................. ,....................

! ............... .... _ •_,_ ,.,.,,_,_

-_

........ _ .................. .-- .,_

ClI_ i ! ..........._ u . 11

....... / i"° 1_ _ " _ _ t.- . ,.,..,,

k_ ]: ] i ,r,. "

tl ,_i

°_"

II -_ .............. "... II

_4.
_ .-_ _,._ ! II___ li __-× •

Ii !11 i tif _.,., °,-,=_ o,,

_P

rL

ss
m

_ _._

'N6_

_.- _ a _

_-_

_",_ _ _

b_

56
- !

u t'4

m

_oo

I I

l a

_x

......................... ".................._..............._"............. _ =..0°-_ o

;_.-_

E

57
c: ×

U..,

58
w v • i

I I I I

I I I I

I I I I I I

I ! 1 I I I

,...-',,,

_m

I I I I I I !

!

m

m

m

J

I I I I I I I

I I ! I 1 I I

u

..__

_×

! I

59
I!

I!

!!

I1I | II!

I!!I!!l

!I H

!!

II

II I I 1 I I I I I

II@u

d

6O
D

I I ! | I I | | ! I

H

F I I

[ 1 I I I I I

0

61
I

!

62
\

,.%,.

I I I I

63
!

I

I !

I I I

,oo

=__.

64
H

_1 ! I i ! ! ! I I I

i I ! I I I i I l I

H

T I ! l ! ! I I I

| I I I I I I

D

\

_I i ! i I I I I I !

,,.11

65
!! I

!I||IIIIIIIII!I

il

IIIIIIIIIIIIIIlI II I i I I I I I I x Uc_

_eN

u ,-

_2

66
I

I

n

! I I I I I ! I I

I I I I I I i I I

n g

I I I ! i ! !

U

°

U..

67
II

I

II

I I I I I I I

I I I I I | I

• s...l

68
H II

,=n

T I I I ! I I I I I

I I I I I I I I I !

g

H

| | | ! | I I I 1 I

[ l I I I I I I ! I

69
I I I ! I I I !

I l I l I I ! I

H

I I

H

I !

! I I I I I ! I

! I I I I ( l I

_... II

0 II

_._'_

._=_

o ,_ ,,,4

_ ,-,

,d

'70
I

g

I I

g

I I I I I I |

| I I I I I !

I I I 1 I I I

il

; ! • w |

71
)4

I i

q

! 1 I 1 I d\" e ,

t I M

I

q

q

q

q

H

.c

0

2-

e4

72
0

I---

0 0 0 0 0

I I

x

0

"1o

0

0

o,l

x

0 0 0 0 0

I I

o

II

I--

0

?3
II II

._ _

_,u o_

_.o

"1:3 iti II

• e" w U

• _1_,

e4

74
\

Fig. 4.2.4. Vortex pairing: Comparison of normalized tcmperaturt contours for a = -3 with the

un-split approach at time f = 160 on • 101 x 101 grid with s = 0.7 for the nonlinear

fields and _ = 0 for the linear fields using ACM66.

75
_o

L_

b_

c_

11t

76
o

0

0 0 0

tN

0

0

0 0

tN

! I

0 0 0 0 0

I I

0

0

o

o°

0'_

d

o

0

I--

o

0

lllllilllJllll

0 0 0 0 0

l I

0

0

_ ^ II

_. _._.

_ x

if,

0

77

f

f

f

_ k

J
I I

0

C

78
Form Approved

REPORT DOCUMENTATION PAGE oMaNo ozo..o,oe

PuOhc tet0orhng burden for this Collection of _nformahOn _g a$1dr_aled to average I hour oar response, qncfud=ng the hme tor reviewing _nstrucIions. searching exist*rig data sources,

qathenng an(::l ma=nla_n_nq Ihe data needed, and completing and reviewing Ihe colleCt=on of _nformahon Send commenl$ regard=ng this burden estimate or any other asloeCt of th*s

ColleChon of _nformaIlon, _ncluding suggesbon$ Jot reducing th_$ burden, to Washington Headquarters Serwces. Dlrecforate for Jnlormalion Operalions and Reports, t215 Jefferson

Daws H_qhway. Su*le t204 Arlin(;iton. VA 22202-4302, and to Ihe Office of Management and Budget. Paperwork Reduchon Profecl (0704-01881. WaShington. OC 20503.

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE

August 1999

4. TITLE AND SUBTITLE

Entropy Splitting and Numerical Dissipation

S. REPORT TechnicalTYPEMemorandum ANDDATESCOVERED

5. FUNDINGNUMBERS

S. AUTHOR(S)

H. C. Yee, M. Vinokur.* and M. J. Djomehri'

*Ames Associate, Ames Research Center. 'Sverdrup Technology Inc.. Ames Research Center.

PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

Ames Research Center

Moffett Field, CA 94035-1000

7.

9. SPONSORING/MONITORINGAGENCY NAME(S)AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-000l

522-31 -12

8. PERFORMING ORGANIZATION

REPORT NUMBER

A-99V0035

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA/TM-1999-208793

11. SUPPLEMENTARY NOTES

Point of Contact: H.C. Yee, Ames Research Center, MS T27B- 1, Moffett Field, CA 94035- I000;

To be presented at the 8th International Symposium on CFD, Sept. 5-10, 1999, Bremen, Germany.

Manuscript submitted to J. Comput. Phys., July 30, 1999.

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified -- Unlimited

Subject Category 02, 34, 64 Distribution: Standard

Availability: NASA CASI (301) 621-0390

(650) 604-4769

2b. DIS'TRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

A rigorous stability estimate for arbitrary order of accuracy of spatial central difference schemes for initial-boundary value

problems of nonlinear symmetrizable systems of hyperbolic conservation laws was established recently by Olsson and Oliger

(19o4) and Olsson (1995), and was applied to the two-dimensional compressible Euler equations for a perfect gas by Gerritsen and

Olsson (1996) and Gerritsen (1996). The basic building block in developing the stability estimate is a generalized energy approach

based on a special splitting of the flux derivative via a convex entropy function and certain homogeneous properties. Due to some

of the unique properties of the compressible Euler equations for a perfect gas. the splitting resulted in the sum of a conservative

portion and a non-conservative portion of the flux derivative, hereafter referred to as the "Entropy Splitting." There are several

potential desirable attributes and side benefits of the entropy splitting for the compressible Euler equations that were not fully

explored in Gerritsen and Olsson. The paper has several objectives. The first is to investigate the choice of the arbitrary parameter

that determines the amount of splitting and its dependence on the type of physics of current interest to computational fluid

dynamics. The second is to investigate in what manner the splitting affects the nonlinear stability of the central schemes for long

time integrations of unsteady flows such as in nonlinear aeroacoustics and turbulence dynamics, tf numerical dissipation indeed is

needed to stabilize the central scheme, can the splitting help minimize the numerical dissipation compared to its un-split cousin?

Extensive numerical study on the vortex preservation capability of the splitting in conjunction with central schemes for long time

integrations will be presented. The third is to study the effect of the non-conservative proportion of splitting in obtaining the

correct shock location for high speed complex shock-turbulence interactions. The fourth is to determine if this method can be

extended to other physical equations of state and other evolutionary equation sets. If numerical dissipation is needed, the Yee.

Sandham. and Djomehri (1999) numerical dissipation is employed. The Yee et al, schemes fit in the Olsson and Oliger framework.

14. SUBJECT TERMS

Entropy variables, Symmetrizable nonlinear conservation laws, Entropy splitting, Thermally perfect

gas. [.,ow dissipation, High order methods. Compact schemes, High order central differencings, Shockcapturing methods. Nonlinear filters. Conservative differencing, TVD schemes. Positive schemes,

END schemes

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION

OF REPORT OF THIS PAGE OF ABSTRACT

Uncla._sitied Unclassified

NSN 7540-01-280-5500

15. NUMBER OF PAGES

83

i

16. PRICE CODE

A05

20. LIMITATION OF ABSTRACT

Standard Form 298(Rev 2-89)

Pt_ctlbo=J b v ANSI _fd Z_l_" I a