Sin-Chung Chang

Glenn Research Center, Cleveland, Ohio

Courant Number and Mach Number

Insensitive CE/SE Euler Solvers

NASA/TM—2005-213868

September 2005

AIAA–2005–4355
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

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

has less stringent limitations on manuscript

length and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

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

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
Sin-Chung Chang

Glenn Research Center, Cleveland, Ohio

Courant Number and Mach Number

Insensitive CE/SE Euler Solvers

NASA/TM—2005-213868

September 2005

National Aeronautics and

Space Administration

Glenn Research Center

Prepared for the

41st Joint Propulsion Conference and Exhibit

cosponsored by the AIAA, ASME, SAE, and ASEE

Tucson, Arizona, July 10–13, 2005

AIAA–2005–4355
Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

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

1. Introdution

The space-time conservation element and solution element (CE/SE) method is a high-resolution and

genuinely multidimensional method for solving conservation laws [1–58]. Its nontraditional features include:

(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution

elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a novel time marching strategy

that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated

without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability

to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a

non-dissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v)

the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as

independent marching variables to be solve for simultaneously. Note that CEs are nonoverlapping space-time

subdomains introduced such that (i) the computational domain can be filled by these subdomains; and (ii)

flux conservation can be enforced over each of them and also over the union of any combination of them.

On the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can

be divided into several component parts with each of them belonging to a unique SE; and (ii) within a SE,

any physical flux vector is approximated using simple smooth functions. In general, a CE does not coincide

with a SE.

Without using flux-splitting or other special techniques, since its inception [1] the unstructured-mesh

compatible CE/SE method has been used to obtain numerous accurate 1D, 2D and 3D steady and unsteady

flow solutions with Mach numbers ranging from 0.0028 to 10 [42]. The phyical phenomena modeled include

traveling and interacting shocks, acoustic waves, shedding vortices, viscous flows, detonation waves, cavitation, flows in fluid film bearings, heat conduction with melting and/or freezing, electrodynamics, MHD

vortex, hydraulic jump, crystal growth, and chromatographic problems [2–58]. In particular, the rather

unique capability of the CE/SE method to resolve both strong shocks and small disturbances (e.g., acoustic

waves) simultaneously [11,13,14] makes it an effective tool for attacking computational aeroacoustics (CAA)

problems. Note that the fact that second-order CE/SE schemes can solve CAA problems accurately is an

exception to the commonly-held belief that a second-order scheme is not adequate for solving CAA problems.

Also note that, while numerical dissipation is needed for shock capturing, it may also result in annihilation of

small disturbances. Thus a solver that can handle both strong shocks and small disturbances simultaneously

must be able to overcome this difficulty.

In spite of its past successes, there is still room for improving the CE/SE method. An example is

the fact that, in a CE/SE simulation with a fixed total marching time, generally the numerical dissipation

NASA/TM—2005-213868 1

It has been known that the space-time CE/SE method can be used to obtain ID, 2D, and 3D steady and 

unsteady flow solutions with Mach numbers ranging from 0.0028 to 10. However, it is also known that a CE/SE 

solution may become overly dissipative when the Mach number is very small. As an initial attempt to remedy this 

weakness, new 1D Courant number and Mach number insensitive CE/SE Euler solvers are developed using several 

key concepts underlying the recent successful development of Courant number insensitive CE/SE schemes. 

Numerical results indicate that the new solvers are capable of resolving crisply a contact discontinuity embedded in 

a flow with the maximum Mach number = 0.01. 

Courant Number and Mach Number Insensitive 

CE/SE Euler Solvers 

Sin-Chung Chang 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 
increases as the value of the Courant number decreases from 1, its maximum stability bound. As such, in a

case with a large Courant number disparity (e.g., a simulation with a highly non-uniform spatial mesh and a

spatially independent time step), the performance sensitivity with respect to the Courant number can lead

to a solution that is highly dissipative in a region where the local Courant number  1. Another example

is the fact that a CE/SE solution may become overly dissipative when Mach number is very small.

By using the fact that each CE/SE scheme is built from a non-dissipative core scheme, robust onedimensional and multidimensional Courant number insensitive schemes were described in [50,52]. As a followup, new one-dimensional Courant number and Mach number insensitive CE/SE Euler solvers have been

developed using several key concepts underlying the development of Courant number insensitive schemes.

Numerical results indicate that the new solvers are capable of resolving crisply a contact discontinuity

embeded in a flow with the maxmum Mach number = 0.01. The rest of the paper is organized as follows.

A review of the existing CE/SE schemes is given in Secs. 2–4. The new Courant number and Mach number

insensitive schemes are described in Sec. 5. Numerical results are presented in Sec. 6. Conclusions and

discussions are given in Sec. 7.

2. Review of the basic 1D CE/SE method

For simplicity, we review the existing CE/SE schemes for the PDE

∂u

∂t

+ a

∂u

∂x

= 0 (2.1)

where a 6= 0 is a constant. Let x1 = x, and x2 = t be considered as the coordinates of a two-dimensional Euclidean space E2. Then, because Eq. (2.1) can be expressed as ∇·

~h = 0 with ~h

def = (au, u), Gauss’ divergence

theorem in the space-time E2 implies that Eq. (2.1) is the differential form of the integral conservation law

I

S(V )

~h · d~s = 0 (2.2)

As depicted in Fig. 1, here (i) S(V ) is the boundary of an arbitrary space-time region V in E2, and (ii)

d~s = dσ ~n with dσ and ~n, respectively, being the area and the unit outward normal of a surface element

on S(V ). Note that: (i) because ~h · d~s is the space-time flux of ~h leaving the region V through the surface

element d~s, Eq. (2.2) simply states that the total space-time flux of ~h leaving V through S(V ) vanishes;

(ii) in E2, dσ is the length of a line segment on the simple closed curve S(V ); and (iii) all mathematical

operations can be carried out as though E2 were an ordinary two-dimensional Euclidean space.

To proceed, let Ω denote the set of all space-time staggered mesh points in E2 (dots in Fig. 2(a)), where

n = 0, ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, j = n±1/2, n±3/2, n±5/2, . . .. Each (j, n) ∈ Ω is associated

with a solution element, i.e., SE(j, n). By definition, SE(j, n) is the interior of the space-time region bounded

by a dashed curve depicted in Fig. 2(b). It includes a horizontal line segment, a vertical line segment, and

their immediate neighborhood.

Let (x,t) ∈ SE(j, n). Then Eq. (2.2) will be simulated numerically assuming that u(x,t) and ~h(x,t),

respectively, are approximated by

u

∗

(x,t; j, n)

def = u

n

j + (ux)

n

j

(x − xj ) + (ut)

n

j

(t − t

n

) (2.3)

and

~h

∗

(x,t; j, n)

def =



au∗(x,t; j, n), u

∗

(x,t; j, n)



(2.4)

Note that (i) u

n

j

, (ux)

n

j

, and (ut)

n

j

are constants in SE(j, n), (ii) (xj ,t

n) are the coordinates of the mesh point

(j, n) with xj = j∆x and t

n = n∆t, and (iii) Eq. (2.4) is the numerical analogue of the definition ~h = (au, u).

Let u = u

∗

(x,t; j, n) satisfy Eq. (2.1) within SE(j, n). Then one has (ut)

n

j = −a (ux)

n

j

. As a result,

Eq. (2.3) reduces to

u

∗

(x,t; j, n) = u

n

j + (ux)

n

j



(x − xj ) − a (t − t

n

)



(2.5)

NASA/TM—2005-213868 2
i.e., u

n

j

and (ux)

n

j

are the only independent marching variables associated with (j, n).

Let E2 be divided into nonoverlapping rectangular regions (see Fig. 2(a)). As depicted in Figs. 2(c)–2(e),

(i) two such regions, i.e., CE−(j, n) and CE+(j, n), are associated with each interior mesh point (j, n) ∈ Ω;

and (ii) CE(j, n) is the union of CE−(j, n) and CE+(j, n).

Given the above preliminaries, we are ready to describe the existing CE/SE solvers for Eq. (2.1).

2.1. The a scheme

Note that, among the line segments forming the boundary of CE−(j, n), AB and AD belong to SE(j, n),

while CB and CD belong to SE(j − 1/2, n − 1/2). Similarly, the boundary of CE+(j, n) belongs to either

SE(j, n) or SE(j + 1/2, n − 1/2). As a result, by imposing two conservation conditions at each (j, n) ∈ Ω,

i.e.,

I

S(CE+(j,n))

~h

∗

· d~s = 0 and I

S(CE−(j,n))

~h

∗

· d~s = 0, (j, n) ∈ Ω (2.6)

and using Eqs. (2.4) and (2.5), one can obtain two equations for the two unknowns u

n

j

and (ux)

n

j

. In fact,

let (i) ν

def = a∆t/∆x, and (ii) for any (j, n) ∈ Ω,

(ux ̄)

n

j

def =

∆x

4

(ux)

n

j

(2.7)

then Eq. (2.6) implies that (i)

u

n

j =

1

2

n

(1 + ν)u

n−1/2

j−1/2 + (1 − ν)u

n−1/2

j+1/2 + (1 − ν

2

)

h

(ux ̄)

n−1/2

j−1/2 − (ux ̄)

n−1/2

j+1/2

io (2.8)

and, assuming |ν| 6= 1, (ii)

(ux ̄)

n

j = (u

a

x ̄

)

n

j

(|ν| 6= 1) (2.9)

with

(u

a

x ̄

)

n

j

def =

1

2

h

u

n−1/2

j+1/2 − u

n−1/2

j−1/2 − (1 + ν)(ux ̄)

n−1/2

j+1/2 − (1 − ν)(ux ̄)

n−1/2

j−1/2

i

(|ν| 6= 1) (2.10)

The a scheme, i.e., the inviscid case of the a-μ scheme [1,3,9], is formed by Eqs. (2.8) and (2.9). Note that,

because

∂u

∂x ̄

=

∆x

4

∂u

∂x

if x ̄

def = x/(∆x/4), the normalized parameter (ux ̄)

n

j may be interpreted as the value at (j, n) of the derivative

of u with respective to the normalized coordinate x ̄. Also note that the superscript symbol “a” in the

parameter (u

a

x ̄

)

n

j

is introduced to remind the reader that Eq. (2.9) is valid for the a scheme.

The review of the a scheme is concluded with the following remarks:

(a) Even though it is introduced to model a single PDE (i.e., Eq. (2.1)) with a single dependent variable

u, the a scheme is formed by two coupled discrete equations (i.e., Eqs. (2.8) and (2.9)) involving two

independent numerical variables u

n

j

and (ux)

n

j

. It is shown in [1] that Eqs. (2.8) and (2.9) are consistent

with a pair of PDEs with one of them being Eq. (2.1).

(b) The a scheme has the simplest stencil, i.e., a triangle with a vertex at the upper time level and the other

two vertices at the lower time level. Furthermore, the number of the independent marching variables

associated with a mesh point (j, n) ∈ Ω is equal to the number of the mesh points at the (n − 1/2)th

time level that are part of the stencil. Note that the same relation also holds for many 2D and 3D

CE/SE schemes [7,8,41].

(c) As shown in [3], the two amplification factors of the a scheme are identical to those of the leapfrog

scheme. As a result, the a scheme is non-dissipative and it is stable if |ν| < 1 (see the additional

discussions given in Sec. 2.2).

(d) Note that derivation of Eqs. (2.8) and (2.9) can be facilitated by the following observations: because

u

∗

(x,t; j, n) is linear in x and t, it can be shown that the total flux of ~h

∗

leaving CE−(j, n) or CE+(j, n)

NASA/TM—2005-213868 3
through any of the four line segments that form its boundary is equal to the scalar product of the vector

~h

∗

evaluated at the midpoint of the line segment and the “surface” vector (i.e., the unit outward normal

multiplied by the length) of the line segment.

(e) Because, for any (j, n) ∈ Ω, the total flux of ~h

∗

leaving each of CE−(j, n) and CE+(j, n) vanishes (see

Eq. (2.6)), CE−(j, n) and CE+(j, n), (j, n) ∈ Ω, will be referred to as the conservation elements (CEs)

of the a scheme. In addition, because (i) the vector ~h

∗ at any surface element lying on any interface

separating two neighboring CEs is evaluated using the information from a single SE, and (ii) the unit

outward normal vector on the surface element pointing outward from one of these two neighboring CEs

is exactly the negative of that pointing outward from another CE, one concludes that the flux leaving one

of these CEs through the interface is the negative of that leaving another CE through the same interface.

As a result, the local conservation relations Eq. (2.6) lead to a global flux conservation relation, i.e., the

total flux of ~h

∗

leaving the boundary of any space-time region that is the union of any combination of

CEs will also vanish. In particular, because CE(j, n) is the union of CE−(j, n) and CE+(j, n),

I

S(CE(j,n))

~h

∗

· d~s = 0, (j, n) ∈ Ω (2.11)

must follow from Eq. (2.6). In fact, it can be shown that Eq. (2.11) is equivalent to Eq. (2.8).

(f) In addition to the non-dissipative a scheme, as will be shown, there is a family of its dissipative extensions

in which only the less stringent conservation condition Eq. (2.11) is assumed [3]. Because Eq. (2.11) is

equivalent to Eq. (2.8), for each of these extensions, u

n

j

is still evaluated using Eq. (2.8) while (ux ̄)

n

j

is

evaluated using an equation different from Eq. (2.9).

2.2. The a- scheme and the c scheme

To proceed, consider any (j, n) ∈ Ω. Then (j ± 1/2, n − 1/2) ∈ Ω. Let

u

0 n

j±1/2

def = u

n−1/2

j±1/2 + (∆t/2)(ut)

n−1/2

j±1/2

(2.12)

With the aid of Eq. (2.7) and the fact that the Courant number ν

def = a∆t/∆x, a substitution of the relation

(ut)

n

j = −a(ux)

n

j

into Eq. (2.12) results in

u

0 n

j±1/2 = (u − 2ν ux ̄)

n−1/2

j±1/2

(2.13)

Note that, to simplify notation, in the above and hereafter we adopt a convention that can be explained

using the expression on the right side of Eq. (2.13) as an example, i.e.,

(u − 2ν ux ̄)

n−1/2

j±1/2 = u

n−1/2

j±1/2 − 2ν(ux ̄)

n−1/2

j±1/2

Also note that, by definition, (j ± 1/2, n) ∈/ Ω if (j, n) ∈ Ω. Thus u

0 n

j±1/2

is associated with a mesh point

∈/ Ω. The reader is warned that similar situations may occur in the rest of this paper.

According to Eq. (2.12), u

0 n

j±1/2

can be interpreted as a first-order Taylor’s approximation of u at

(j ± 1/2, n). Thus

(u

c

x ̄

)

n

j

def =

∆x

4

 

u

0 n

j+1/2 − u

0 n

j−1/2

∆x

!

(2.14)

is a central-difference approximation of ∂u/∂x ̄ at (j, n). Note that: (i) the superscript “c” is used to remind

the reader of the central-difference nature of the term (u

c

x ̄

)

n

j

; and (ii) by using Eqs. (2.13), (2.14) and (2.10),

one has

(u

c

x ̄

)

n

j =

1

4

h

(u − 2νux ̄)

n−1/2

j+1/2 − (u − 2νux ̄)

n−1/2

j−1/2

i

(2.15)

NASA/TM—2005-213868 4
and

(u

c

x ̄ − u

a

x ̄

)

n

j =

1

2

h

(ux ̄)

n−1/2

j+1/2 + (ux ̄)

n−1/2

j−1/2

i

−

1

4



u

n−1/2

j+1/2 − u

n−1/2

j−1/2



(2.16)

The a- scheme is formed by Eq. (2.8) and

(ux ̄)

n

j = (u

a

x ̄

)

n

j + 2(u

c

x ̄ − u

a

x ̄

)

n

j

(2.17)

where  is a real number. Obviously the a- scheme reduces to the a scheme when  = 0. Also, for the case

 = 1/2, Eq. (2.17) reduces to

(ux ̄)

n

j = (u

c

x ̄

)

n

j

(2.18)

Because (u

c

x ̄

)

n

j

represents a central-difference approximation, hereafter, to simplify its frequent references,

the special a- scheme with  = 1/2 will be referred to as the c scheme.

To proceed, several key remarks about the a- scheme are presented:

(a) At each mesh point (j, n) ∈ Ω, Eqs. (2.8) and (2.9) are the results of Eq. (2.6). Because Eq. (2.17) does

not reduce to Eq. (2.9) except in the special case  = 0, at each mesh point (j, n) ∈ Ω, generally the

a- scheme satisfies only the single conservation condition Eq. (2.11) (which is equivalent to Eq. (2.8))

rather than the two consevation conditions Eq. (2.6). However, because (u

a

x ̄

)

n

j

generally is present on

the right side of Eq. (2.17), the a- scheme generally will still be burdened with the cost of solving two

conservation conditions at each mesh point. The exception occurs only for the special case  = 1/2

(i.e., the c scheme) in which Eq. (2.17) reduces to Eq. (2.18). As it turns out, implementation of

a multidimensional Euler version of the c scheme does not require inverting any system of equations

while a similar implementation involving a version of any other a- scheme ( 6= 1/2) generally requires

inverting, per mesh point and per time step, a system of several linear equations (to be exact, a system

of eight and fifteen equations, respectively, for 2D and 3D Euler equations) [7,8,41]. As such, it is much

more cost effective to use a multidimensional Euler version of the c scheme than using that of any other

a- scheme. Partly for this reason, extensions of the c scheme have been used extensively.

(b) For the a- scheme, it is shown in [3] that the principal and spurious amplification factors per ∆t,

respectively, are (λ+)

2 and (λ−)2 with

λ±(, ν, θ)

def =  cos(θ/2) − iν sin(θ/2) ±

q

(1 − )



(1 − )cos2(θ/2) + (1 − ν

2)sin2

(θ/2)(2.19)

Here (i) i

def =

√

−1, and (ii) θ, −π < θ ≤ π, is the phase angle variation per ∆x. In addition, it is shown

that (i) the necessary and sufficient conditions for the stability of the a- scheme are

0 ≤  ≤ 1, and |ν| < 1 (2.20)

and (ii) the a- scheme becomes progressively diffusive as the value of  increases from 0 to 1. Note

that, unless specified otherwise, in the remainder of the paper the ranges of , ν and θ, respectively, are

defined by Eq. (2.20) and −π < θ ≤ π.

(c) Let k be a constant. Then u = e

ik(x−at)

represents a plane wave solution to Eq. (2.1). For this solution

the exact amplification factor per ∆t

def =

e

ik[x−a(t+∆t)]

e

ik(x−at)

= e

−ika∆t = e−iνθ

(2.21)

where θ = k∆x.

(d) According to Eq. (2.19), [λ±(0, ν, θ)]2, the amplification factors of the a scheme (which corresponds to

the case  = 0) per ∆t, have the following properties:

|[λ±(0, ν, θ)]2| = 1 (2.22)

NASA/TM—2005-213868 5
lim ν→±1

[λ+(0, ν, θ)]2 = e

∓iθ (2.23)

lim ν→±1

[λ−(0, ν, θ)]2 = e

±iθ (2.24)

and

[λ±(0, 0, θ)]2 = 1 (2.25)

On the other hand, e

−iνθ

, the exact amplification factor per ∆t, has the following properties:

|e

−iνθ

| = 1 (2.26)

lim

ν→±1

e

−iνθ = e∓iθ (2.27)

and

e

−iνθ = 1 if ν = 0 (2.28)

For the a scheme, Eqs. (2.22)–(2.28) imply that: (i) the two amplification factor of the scheme, and

the exact amplification factor all have the same constant absolute value (= 1) and, thus, the scheme is

non-dissipative; (ii) in the limit of |ν| → 1 (i.e., ν → 1 or ν → −1), the principal amplification factor is

identical to the exact amplification factor and, thus, the former has no dissipative or dispersive error in

this limit; (iii) also in the limit of |ν| → 1, the phase angle associated with the spurious amplification

factor is exactly the negative of that associated with the exact amplification factor and, thus, the

spurious amplification factor has a large dispersive error in this limit except when |θ|  1(i.e., when

the wavelengths of the errors  1); and (iv) when ν = 0, the two amplification factors of the scheme,

and the exact amplification factor are all equal to 1 and, thus, the two amplification factors of the

scheme have no dissipative or dispersive error if ν = 0. Because the accuracy of a scheme is essentially

hinged on the behaviors of the principal amplification factor [1], according to the facts stated above,

the a scheme tends to become very accurate when |ν| approaches 1 or 0. However, the short-wavelength

errors associated with the spurious amplification factor (which could be introduced at t = 0 as a result of

an inaccurate initial-value specification [1]) may appear in a solution as persistent (i.e., non-dissipative)

numerical wiggles when |ν| approaches 1 [1,9].

(e) According to Eq. (2.19), [λ±(1/2, ν, θ)]2, the amplification factors of the c scheme (which corresponds

to the case  = 1/2) per ∆t, have the following properties:

lim ν→±1[λ+(1/2, ν, θ)]2 = e

∓iθ (2.29)

lim

ν→±1

[λ−(1/2, ν, θ)]2 = − sin2(θ/2) (2.30)

and

[λ±(1/2, 0, θ)]2 =

1

2

h

1 ± cos(θ/2)p2 − cos2(θ/2)i

(2.31)

For the c scheme, Eqs. (2.27)–(2.31) imply that: (i) in the limit of |ν| → 1, the principal amplification

factor is identical to the exact amplification factor and, thus, the former has no dissipative or dispersive

errors in this limit; (ii) also in the limit of |ν| → 1, the spurious amplification factor has large dissipative

and dispersive errors; and (iii) when ν = 0, both the principal and spurious amplification factors

generally have large dissipative errors but no dispersive errors. According to the facts stated above, like

the a scheme, the c scheme also tends to become very accurate when |ν| approaches 1. However, unlike

the a scheme, the errors associated with the spurious amplification factor of the c scheme generally do

die out rapidly when |ν| approaches 1. Also, in sharp contrast to the a scheme, the c scheme becomes

highly dissipative when ν approaches 0.

From the above discussions, one concludes that:

NASA/TM—2005-213868 6
(a) The advantages of the a scheme include: (i) it is non-dissipative throughout the range of  and ν defined

in Eq. (2.20); and (ii) when the value of |ν| is close to 0 or 1, the scheme is very accurate. On the other

hand, its disadvantages include: (i) because it is non-dissipative, its extensions for nonlinear equations

generally are unstable; (ii) when the value of |ν| is close to 1, the short-wavelength errors associated

with the spurious amplification factors will not die out rapidly and, therefore, appear in a solution as

persistent numerical wiggles; and (iii) comparing with the c scheme, it costs more to implement.

(b) The advantages of the c scheme include: (i) when the value of |ν| is close to 1, it is very accurate and

the short-wavelength errors associated with the spurious amplification factor also die out rapidly; (ii)

because it is dissipative, its extensions for nonlinear equations generally are stable; and (iii) in terms of

ease of implementation and computer cost, it is much more superior than any other a- scheme. On the

other hand, the c scheme has a serious disadvantage, i.e., it is very dissipative when ν approaches 0.

In Sec. 3, it will be shown that new solvers of Eq. (2.1) can indeed be constructed such that they possess all

the advantages but none of the disadvantages listed above. Specifically, each of these solvers will be formed

by Eq. (2.8) and a new equation in which (ux ̄)

n

j

is evaluted using a simple central-differencing procedure

similar to that used to obtain (u

c

x ̄

)

n

j

. In addition, (ux ̄)

n

j

so obtained will be (i) identical to (u

c

x ̄

)

n

j

in the limit

of |ν| → 1, and (ii) identical to (u

a

x ̄

)

n

j when ν = 0. As such, each of these new solvers (i) is comparable to

the c scheme in ease of implementation; (ii) becomes the c scheme in the limit of |ν| → 1; and (iii) becomes

the a scheme when ν = 0.

2.3. The w-α scheme—a special wiggle-suppressing scheme

If discontinuities are present in a numerical solution, any a- scheme such as the c scheme is not equipped

to suppress numerical wiggles that generally appear near these discontinuities. To serve as a preliminary

for future development, here we shall briefly review an extension of the c scheme which was introduced as a

remedy for this deficiency [3,35].

To proceed, let

(ux ̄−)

n

j

def =

∆x

4

 

u

n

j − u

0 n

j−1/2

∆x/2

!

(2.32)

and

(ux ̄+)

n

j

def =

∆x

4

 

u

0 n

j+1/2 − u

n

j

∆x/2

!

(2.33)

i.e., (ux ̄−)

n

j

and (ux ̄+)

n

j

are numerical analogues of ∂u/∂x ̄ at (j, n) evaluated from the left and the right,

respectively. It can be shown that

(u

c

x ̄

)

n

j =

1

2

(ux ̄− + ux ̄+)

n

j

(2.34)

i.e., (u

c

x ̄

)

n

j

is the simple average of (ux ̄−)

n

j

and (ux ̄+)

n

j

. As such, the c scheme can be extended by replacing

(u

c

x ̄

)

n

j

in Eq. (2.18) with an weighted average of (ux ̄−)

n

j

and (ux ̄+)

n

j

. In other words, the resulting extension

is formed by Eq. (2.8) and

(ux ̄)

n

j = (w−)

n

j

(ux ̄−)

n

j + (w+)

n

j

(ux ̄+)

n

j

(2.35)

where (w−)

n

j

and (w+)

n

j

, the weight factors associated with (ux ̄−)

n

j

and (ux ̄+)

n

j

respectively, must satisfy the

condition

(w−)

n

j + (w+)

n

j = 1 (2.36)

at all (j, n) ∈ Ω. In addition, the expression on the right side of Eq. (2.35) represents an interpolation (rather

than an extrapolation) of (ux ̄−)

n

j

and (ux ̄+)

n

j

if and only if

(w−)

n

j ≥ 0 and (w+)

n

j ≥ 0 (2.37)

For real variables x−, x+, and α ≥ 0, let W− and W+ be the functions defined by: (i) W−(x−, x+; α) =

W+(x−, x+; α) = 1/2 if x− = x+ = 0; and (ii)

W−(x−, x+; α) =

|x+|

α

|x−|

α + |x+|α

and W+(x−, x+; α) =

|x−|

α

|x−|

α + |x+|α

(2.38)

NASA/TM—2005-213868 7
if either x− 6= 0 or x+ 6= 0. Furthermore, let

(w±)

n

j = W±((ux ̄−)

n

j

,(ux ̄+)

n

j

, α) (2.39)

Then (w−)

n

j

and (w+)

n

j

so defined satisfy Eqs. (2.36) and (2.37) and have the property that

(w−)

n

j = (w+)

n

j = 1/2 if α = 0 or |(ux ̄−)

n

j

| = |(ux ̄+)

n

j

| (2.40)

Note that: (i) to avoid dividing by zero, in practice a small positive number such as 10−20 is added to each

of the denominators in Eq. (2.38); and (ii) the special cases of Eq. (2.38) with α = 1 and α = 2 are used in

the slope-limiter proposed by van Leer [59], and van Albada et al. [60].

An extension of the c scheme is formed by Eqs. (2.8) and (2.35) with (w−)

n

j

and (w+)

n

j being defined

by Eq. (2.39). Because it involves an weighted average which is dependent on a parameter α, hereafter the

scheme is referred to as the w-α scheme. Let α > 0 and |(ux ̄−)

n

j

| 6= |(ux ̄+)

n

j

|. Then Eq. (2.38) implies that,

of (ux ̄−)

n

j

and (ux ̄+)

n

j

, the one with smaller absolute value is associated with an weight factor > 1/2. This

observation coupled with Eqs. (2.34)–(2.37) leads to the conclusion that, of (ux ̄−)

n

j

and (ux ̄+)

n

j

, (ux ̄)

n

j will

have an algebraic value closer to the one with smaller absolute value if (ux ̄)

n

j

is evaluated as an weighted

average of (ux ̄−)

n

j

and (ux ̄+)

n

j

according to Eq. (2.35). As a result, (ux ̄)

n

j

so evaluated has a smaller absolute

value than that evaluated using Eq. (2.18). In turn, numerical wiggles or overshoots can be annihilated by

the additional numerical dissipation introduced as a result of this local “flattening” of (ux ̄)

n

j

. It has been

shown numerically that the extension is stable if |ν| < 1 and α ≥ 0. Moreover, as a result of Eqs. (2.18),

(2.34), (2.35) and (2.40), (i) the extension reduces to the c scheme when α = 0; and (ii) even if α > 0,

the extension behaves very much like the c scheme in any smooth solution region (where the condition

(ux ̄−)

n

j = (ux ̄+)

n

j more or less prevails) or at a solution extremum (where the condition (ux ̄−)

n

j = −(ux ̄+)

n

j

more or less prevails). As such, the wiggle-suppressing power of the extension takes effect only if α > 0 and

only in a solution region where |(ux ̄−)

n

j

| and |(ux ̄+)

n

j

| differ substantially.

3. The c-τ and c-τ

∗

schemes

In this section, the ideal solvers of Eq. (2.1) mentioned at the end of Sec. 2.2 will be constructed. As a

preliminary, we shall show that (u

a

x ̄

)

n

j

can also be cast into a central-difference form when ν = 0.

To proceed, note that by assumption a 6= 0. Thus ν = 0 if and only if ∆t = 0. Because |EF| = |AD| =

|CB| = 0 (see Figs. 2(c,d)) when ∆t = 0, the two conservation conditions given in Eq. (2.6) for the case

ν = 0, respectively reduce to the following conditions: (i) the flux leaving CE+(j, n) through the top face

AF is equal to that entering the same CE through the bottom face ED; and (ii) the flux leaving CE−(j, n)

through the top face AB is equal to that entering the same CE through the bottom face CD. According

to Remark (d) given at the end of Sec. 2.1, the flux leaving CE+(j, n) through the top face AF is equal to

the value of u

∗ at the midpoint of AF (evaluated using the marching variables at point A) multiplied by

|AF|, while that entering it through the bottom face ED is equal to the value of u

∗ at the midpoint of ED

(evaluated using the marching variables at point E) multiplied by |ED|. With the aid of these observations

and the fact that |AF| = |ED|, the above condition (i) implies that, when ν = 0, the value of u

∗ at the

midpoint of AF evaluated using the marching variables at point A is equal to that at the midpoint of ED

evaluated using the marching variables at point E. As such, the first conservation condition in Eq. (2.6) is

equivalent to

(u + ux ̄)

n

j = (u − ux ̄)

n−1/2

j+1/2

(ν = 0) (3.1)

if ν = 0. Similarly, by using the above condition (ii), it can be shown that the second conservation condition

in Eq. (2.6) is equivalent to

(u − ux ̄)

n

j = (u + ux ̄)

n−1/2

j−1/2

(ν = 0) (3.2)

if ν = 0. Because Eqs. (2.8) and (2.9) (which form the a scheme) are equivalent to Eq. (2.6) if |ν| 6= 1, they

must be equivalent to Eqs. (3.1) and (3.2) when ν = 0. In fact, by substracting Eq. (3.2) from Eq. (3.1), one

obtains Eq. (2.9) where (u

a

x ̄

)

n

j

is the reduced form of Eq. (2.10) for the case ν = 0, i.e.,

(u

a

x ̄

)

n

j =

1

2

h

(u − ux ̄)

n−1/2

j+1/2 − (u + ux ̄)

n−1/2

j−1/2

i

(ν = 0) (3.3)

NASA/TM—2005-213868 8
Moreover, by summing over Eqs. (3.1) and (3.2), one has the reduced form of Eq. (2.8) for the case ν = 0.

With the aid of Eq. (3.3) and the facts that: (i) (u−ux ̄)

n−1/2

j+1/2

and (u+ux ̄)

n−1/2

j−1/2

, respectively, represent

an approximation of u at the midpoint of ED and that at the midpoint of CD (see Fig. 2(c,d)); and (ii)

the distance between the two midpoints referred to above is ∆x/2, it becomes obvious that, for the special

case ν = 0, (u

a

x ̄

)

n

j

is indeed a central-difference approximation of ∂u/∂x ̄ at (j, n − 1/2) (which coincides with

(j, n) when ν = 0). QED

According to the above discussions, construction of the ideal solvers defined at the end of Sec. 2.2 is

hinged on finding central-difference approximations for (ux ̄)

n

j

such that each approximation (i) becomes

(u

c

x ̄

)

n

j

in the limit of |ν| → 1, and (ii) reduces to the expression on the right side of Eq. (3.3) when ν = 0.

As a result of these observations, these new solvers can easily be constructed as the subschemes of the c-τ

scheme, a new class of CE/SE solvers for Eq. (2.1) to be described immediately.

3.1. The c-τ scheme

To proceed, refer to Fig. 3. Here M+ and M− denote the midpoints of AF and AB, respectively. Also

P

+ and P− are two points on BF that satisfy the following conditions: (i) P

+ coincides with M+ if and

only if P

− coincides with M−; (ii) P+ is to the right (left) of M+ if and only if P− is to the left (right)

of M−; and (iii) |M+P +| = |M−P −|, i.e., M+P + and M−P − have the same length. In addition, let the

parameter τ be defined by: (i) τ = 0 if P

+ coincides with M+; (ii) τ ∆x/4 = |M+P +| if P

+ is to the

right of M+; and (iii) τ ∆x/4 = −|M+P +| if P

+ is to the left of M+. Obviously, it follows from the above

definitions that (i) τ = 0 if P

− coincides with M−; (ii) τ ∆x/4 = |M−P −| if P

− is to the left of M−; and

(iii) τ ∆x/4 = −|M−P −| if P

− is to the right of M−.

Moreover, let

u(P

+)

def = [u + (∆t/2)ut − (1 − τ )ux ̄]

n−1/2

j+1/2

(3.4)

and

u(P

−)

def = [u + (∆t/2)ut + (1 − τ )ux ̄]

n−1/2

j−1/2

(3.5)

With the aid of Eq. (2.7), it is seen that u(P

+) is a first-order Taylor’s approximation of u at P+ evaluated

using the marching variables at point E, while u(P

−) is a first-order Taylor’s approximation of u at P−

evaluated using the marching variables at point C. Also note that, by using the relation (ut)

n

j = −a(ux)

n

j

,

Eqs. (3.4) and (3.5), respectively, can be simplified as

u(P

+) = [u − (2ν + 1 − τ )ux ̄]

n−1/2

j+1/2

(3.6)

and

u(P

−) = [u − (2ν − 1 + τ )ux ̄]

n−1/2

j−1/2

(3.7)

At this juncture, note that P

+ and P− generally lie outside of SE(j+1/2, n−1/2) and SE(j−1/2, n−1/2),

respectively. Yet here, by definition, u(P

+) and u(P−) are evaluated as though P+ and P− lie within

(j +1/2, n−1/2) and (j −1/2, n−1/2), respectively. At first glance, the current practice is inconsistent with

a previously established rule. However, as explained by the reasons given below, the definition of u(P

+) and

u(P

−) is perfectly legitimate:

(a) Recall that solution elements were introduced such that the boundary of a CE can be divided into

several component parts with each of them belonging to a unique solution element. As such, the flux

over a component part that belongs to a special solution element, say SE(j, n), can be unambiguously

determined in terms of the marching variables at the mesh point (j, n). In other words, in related to

evaluation of any flux conservation condition over any CE, Eqs. (2.3)–(2.5) can be applied only to a

point (x,t) ∈ SE(j, n).

(b) On the other hand, u(P

+) and u(P−) introduced here have nothing to do with flux evaluation. In fact,

they will be used only in the construction of some numerical analogues of ∂u/∂x at (j, n).

NASA/TM—2005-213868 9
To proceed, note that: (i) the mesh point (j, n) (i.e., point A depicted in Fig. 3) is the midpoint of

P −P +, and (ii) |P −P +| = (1 + τ )∆x/2. Thus

(uˆx ̄)

n

j

def =

∆x

4



u(P

+) − u(P−)

(1 + τ )∆x/2



(τ 6= −1) (3.8)

represents a central-difference approximation of ∂u/∂x ̄ at the mesh point (j, n). Thus the new scheme formed

by Eq. (2.8) and

(ux ̄)

n

j = (uˆx ̄)

n

j

(3.9)

represents a solver for Eq. (2.1). Because (i) (uˆx ̄)

n

j

represents a central-difference approximation of (ux ̄)

n

j

,

and (ii) the approximation is associated with the parameter τ , hereafter the new scheme will be referred to

as the c-τ scheme.

To explore the c-τ scheme, note that Eqs. (3.6)–(3.8) can be combined to yield

(uˆx ̄)

n

j =

1

2(1 + τ )

n

[u − (2ν + 1 − τ )ux ̄]

n−1/2

j+1/2 − [u − (2ν − 1 + τ )ux ̄]

n−1/2

j−1/2

o

(τ 6= −1) (3.10)

Moreover, by using Eqs. (2.10), (2.16) and (3.10), one has

(uˆx ̄)

n

j = (u

a

x ̄

)

n

j +

2τ

1 + τ

(u

c

x ̄ − u

a

x ̄

)

n

j −

ν(1 − τ )

2(1 + τ )

h

(ux ̄)

n−1/2

j+1/2 − (ux ̄)

n−1/2

j−1/2

i

(3.11)

By comparing Eq. (3.11) with (2.17), one concludes that the c-τ scheme generally is different from the a-

scheme. In fact, a special case of the c-τ scheme can be turned into that of the a- scheme and vice versa if

and only if either (i) τ = 1 or (ii) ν = 0. For the case τ = 1, Eq. (3.11) implies that (uˆx ̄)

n

j = (u

c

x ̄

)

n

j

. In other

words, the c scheme is the special case of the c-τ scheme with τ = 1, a fact that can also be deduced from

the observation that the points P

+ and P− depicted in Fig. 3, respectively, coincide with points F and B

(i.e., the mesh points (j + 1/2, n) and (j − 1/2, n)) if τ = 1. On the other hand, it is seen that, when ν = 0,

the c-τ scheme become the a- scheme with  = 2τ/(1 + τ ). In fact one can further deduce that c-τ scheme

reduces to the a scheme if and only if ν = τ = 0.

Because the c-τ scheme is formed by two rather complicated equations (i.e., Eqs. (2.8) and (3.9))

involving two parameters ν and τ , it were not expected that its von Neumann stability conditions could

be cast into an explicit analytical form. However, it is shown rigorously in [55] that the c-τ scheme is von

Neumann stable if and only if

ν

2 ≤ 1, τ ≥ τo(ν2

), and (ν

2

, τ ) 6= (1, 1) (3.12)

where

τo(s)

def =







0 if s = 0

4 − s − 2

p

2(2 − s − s

2)

s

if 0 < s ≤

3

11

s − 1 +

√

1 − 2s + 5s

2

2s

if 3

11 ≤ s ≤ 1

(3.13)

Note that:

(a) It can be shown that [55]: (i) τo(s) is continuous at s = 0; (ii) τo(s) is consistently defined at s = 3/11;

(iii)

lim

s→ 3

11

−

τ

0

o

(s) = lim

s→ 3

11

+

τ

0

o

(s) = 121/90

where τ

0

o

(s)

def = dτo(s)/ds; (iv) τo(s) is strictly montonically increasing in the interval 0 < s < 1; (v)

τo(1) = 1; and (vi)

s < τo(s) <

√

s, 0 < s < 1 (3.14)

NASA/TM—2005-213868 10
(b) For any given fixed value of |ν| < 1, the c-τ scheme tends to become more dissipative as the value of τ

increases from τo(ν

2

), its minimum stable value.

(c) The function τo defined here is different from the function τo defined in [52]. In fact, for any ν ∈ [0, 1],

the value of τo(ν

2

) defined here is identical to that of τo(|ν|) defined in [52]. Because the analytical

ν-τo(ν

2

) relation depicted in Fig. 4 here is virtually identical to the numerical ν-τ (|ν|) relation depicted

in Fig. 4 of [52], the analytical stability conditions given here are in complete agreement with those

generated numerically and stated in Eq. (3.12) of [52].

With the above preliminaries, the ideal solvers of Eq. (2.1) defined at the end of Sec. 2.2 will be

constructed in Sec. 3.2.

3.2. The c-τ

∗

schemes

The value of τ used in the c-τ scheme generally can be chosen independent of ν. Here we will introduce

a subset of the c-τ scheme in which τ is a function of ν

2

for each member of this subset. As a preliminary,

note that, by using the properties of τo(s) presented earlier, it can be shown that there exist infinitely many

choices of a strictly monotonically increasing smooth function h(s), 0 ≤ s < 1, which possesses the following

properties:

h(0) = 0; lim

s→1−

h(s) = 1; and h(s) ≥ τo(s) if 0 < s < 1 (3.15)

For each function h(s) that satisfies the above conditions, each member of the subset referred to earlier is

defined using the relation

τ = h(ν

2

) (ν

2 < 1) (3.16)

Note that, using the definition of h and Eq. (3.16), one can easily infer from Fig. 3 a simple relation between

the value of ν

2 and the locations of P+ and P−, i.e., as the value of ν2

increases from 0 to 1, P

+ will move

away from M+ and edge toward the mesh point (j + 1/2, n) while P

− will move away from point M− and

edge toward the mesh point (j − 1/2, n). Also note that, by using the above preliminaries, one can show

that: (i)

τ = 0 if ν = 0 (3.17)

(ii)

lim

ν

2→1−

τ = 1 (3.18)

and (iii)

τ ≥ τo(ν

2

) (ν

2 < 1) (3.19)

Recall that (i) (uˆx ̄)

n

j = (u

a

x ̄

)

n

j

if τ = ν = 0; and (ii) (uˆx ̄)

n

j = (u

c

x ̄

)

n

j

if τ = 1. As such, Eqs. (3.17) and

(3.18) imply that, for each member in the subset, (i) (uˆx ̄)

n

j = (u

a

x ̄

)

n

j

if ν = 0; and (ii) (uˆx ̄)

n

j = (u

c

x ̄

)

n

j

in

the limit of ν

2 → 1−. In other words, all members in the subset are ideal solvers in the domain ν2 < 1.

Moreover, by using Eq. (3.12) and (3.19), one can also show that these ideal solvers are also stable in the

same domain. Hereafter, each of these ideal solvers will be referred to as an ideal c-τ

∗

scheme.

Corresponding to infinitely many choices of h, there are infinitely many different ideal c-τ

∗

schemes. In

particular, because

h(s) =

√

s (0 ≤ s < 1) (3.20)

is a legitimate choice (see Eq. (3.14)), Eq. (3.16) implies that an ideal c-τ

∗

scheme can be defined using the

relation

τ = |ν| (ν

2 < 1) (3.21)

Any ideal c-τ

∗

scheme described above meets all the requirements of an ideal solver defined in Sec. 2.2.

However, because of stability problem, an ideal c-τ

∗

scheme and its multidimensional extensions may not

be robust enough for some complicated real-world applications. To overcome this difficulty, we introduce a

special c-τ scheme which is defined by Eqs. (2.8), (3.9), and (3.10) with

τ = β |ν| (β ≥ 1; ν

2 < 1) (3.22)

NASA/TM—2005-213868 11
where β ≥ 1 is an adjustable parameter. For this scheme, Eq. (3.17) is also valid. Thus it becomes the

nondissipative a scheme when ν = 0 and, therefore, it will not become overly dissipative as ν → 0 (i.e., the

scheme’s numerical dissipation is Courant number insensitive). Moreover, because β ≥ 1, Eqs. (3.12) and

(3.14) along with the comment (b) given following Eq. (3.14) imply that the scheme is stable if ν

2 < 1 and

it becomes more dissipative (i.e., more stable) as the value of β increases. However, also because β ≥ 1,

lim

|ν|→1−

τ = β ≥ 1 (3.23)

As such, except for the case β = 1, the scheme does not become the c scheme in the limit of |ν| → 1

−, i.e.,

it does not meet all the requirements of an ideal solver.

Note that, unless specified otherwise, in Sec. 4 we consider only the cases in which either Eq. (3.16) or

Eq. (3.22) is assumed.

4. Extensions of the c-τ

∗

scheme and related weighted averagings

To proceed, let

(uˆx ̄−)

n

j

def =

∆x

4

 

u

n

j − u(P

−)

(1 + τ )∆x/4

!

(4.1)

and

(uˆx ̄+)

n

j

def =

∆x

4

 

u(P

+) − un

j

(1 + τ )∆x/4

!

(4.2)

Because |AP −| = |AP +| = (1 + τ )∆x/4 (see Fig. 3), it is easy to see that (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j

are two

one-sided difference approximations of ∂u/∂x ̄ at the mesh point (j, n) with one being evaluated from the left

and another from the right. Also, it follows immediately from Eqs. (3.8), (4.1) and (4.2) that

(uˆx ̄)

n

j =

1

2



(uˆx ̄−)

n

j + (uˆx ̄+)

n

j



(4.3)

Moreover, by (i) substituting Eqs. (2.8), (3.6) and (3.7) into Eqs.(4.1) and (4.2), and (ii) using Eqs. (3.10)

and (3.17), one arrives at the conclusion that

(uˆx ̄−)

n

j = (uˆx ̄+)

n

j = (uˆx ̄)

n

j

(ν = 0) (4.4)

when ν = 0.

With the above preliminaries, several extensions of the c-τ

∗

scheme will be constructed in the following

subsections.

4.1. Scheme w-1

A comparison of Eqs. (4.1)–(4.3) with Eqs. (2.32)–(2.34) reveals that an obvious extension of the c-τ

∗

scheme can be obtained by replacing (ux ̄−)

n

j

and (ux ̄+)

n

j

in Eqs. (2.35) and (2.39) with (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j

,

respectively. In other words, the new extension is formed by Eq. (2.8) and

(ux ̄)

n

j = (w−)

n

j

(uˆx ̄−)

n

j + (w+)

n

j

(uˆx ̄+)

n

j

(4.5)

with

(w±)

n

j = W±((uˆx ̄−)

n

j

,(uˆx ̄+)

n

j

, α) (4.6)

Because the scheme is the first extension of the c-τ

∗

scheme in which (ux ̄)

n

j

is expressed as an weighted

average of (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j

, for simplicity, hereafter it will be referred to as Scheme w-1. It has been

shown numerically that Scheme w-1 is stable if |ν| < 1 and α ≥ 0.

Note that, as a result of Eqs. (2.38), (4.4) and (4.6), one concludes that, for any given α ≥ 0, (w−)

n

j =

(w+)

n

j = 1/2 if ν = 0. In other words, for Scheme w-1, the “weighted” average on the right side of Eq. (4.5)

NASA/TM—2005-213868 12
becomes a simple average if ν = 0. According to an explanation given in the last paragraph of Sec. 2, this

implies that Scheme w-1 will lose its capability to suppress wiggles or overshoots when ν becomes small.

For this reason, even though the Euler version of Scheme w-1 performs much better than that of the special

scheme referred to in Sec. 2.3 in its ability to resolve shocks and contact discontinuities crisply in a wide

range (from 1 to less than 0.001) of the global Courant number (i.e., the maximal value of local Courant

numbers), it has a serious shortcoming, i.e., wiggles or overshoots can appear near a discontinuity in a

generated solution when the local Courant number there becomes extremely small. In the following, it will

be shown that this weakness can be overcome by simple modifications of Eq. (4.6).

4.2. Scheme w-2

A new scheme, referred to as Scheme w-2 is formed by Eqs. (2.8) and (4.5) with (w±)

n

j being given

by Eq. (2.39). In other words, although (ux ̄)

n

j

is still constructed as an weighted average of (uˆx ̄−)

n

j

and

(uˆx ̄+)

n

j

, the associated weight factors (w±)

n

j

are evaluated using (ux ̄−)

n

j

and (ux ̄+)

n

j

. Because the last two

parameters, respectively, are identical to the special cases of (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j with τ

def = 1 (see Eqs. (2.12),

(2.32), (2.33), (3.4), (3.5), (4.1), and (4.2)), their values do not vary with ν. As such, (w±)

n

j

6= 1/2 and

therefore the weighted average on the right side of Eq. (4.5) will not turn into a simple average when ν = 0.

In other words, Scheme w-2 is still capable of annihilating the numerical wiggles near a discontinuity even if

ν becomes small. It has been shown numerically that Scheme w-2 again is stable if |ν| < 1 and α ≥ 0.

Note that a possible drawback of Scheme w-2 is that the relation |(ux ̄−)

n

j

| < |(ux ̄+)

n

j

| (|(ux ̄−)

n

j

| >

|(ux ̄+)

n

j

|) does not automatically follow from |(uˆx ̄−)

n

j

| < |(uˆx ̄+)

n

j

| (|(uˆx ̄−)

n

j

| > |(uˆx ̄+)

n

j

|) and vice versa.

As a result, at some local mesh points, it may happen that, of (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j

, the one with smaller

absolute value may not be associated with a weight factor > 1/2. According to a discussion given in the

last paragraph of Sec. 2, this implies that there is no guarantee that, at all localities, the weighted-averaging

induced numerical dissipation will be available to suppress wiggles or overshoots. Despite this possible failing,

fortunately it has been demonstrated numerically that, not only are they capable of suppressing wiggles or

overshoots robustly, Scheme w-2 and its Euler extensions are also highly accurate.

In the following, schemes that overcome the weakness of Scheme w-1 and also avoid the theoretically

possible failing associated with Scheme w-2 will be constructed using new weighted-averaging formulae more

advanced than that given in Eq. (2.38).

4.3. New weighted-averaging techniques

To pave the way, first we shall discuss a limitation of Eq. (2.38) as a generator of weight factors.

Let x± 6= 0. Then, for a given α > 0, obviously W− → 1/2 and W+ → 1/2 as |x+/x−| → 1. As such,

when |x+/x−| is very close to 1, then both W− and W+ will be very close to 1/2 unless α  1. As a result,

in case that (i) (uˆx ̄±)

n

j

6= 0, (ii) |(uˆx ̄+)

n

j

/(uˆx ̄−)

n

j

| is very close to 1; and (iii) Eqs. (4.6) is assumed, then the

only way to prevent the weighted average that appears on the right side of Eq. (4.5) from becoming almost

a simple average is to increase the value of α used. However, this approach may be impracticable because

numerical evaluation of a quantity such as x

α for any real number x generally is hampered by round-off

errors and thus becomes highly inaccurate if the value of α becomes too large, say 100. It is the purpose

of this subsection to introduce new weighted-averaging techniques that do not have the limitaton discussed

above.

For motivation, note that Eqs. (4.5) and (4.6) can be expressed as

(ux ̄)

n

j = w1x1 + w2x2 (4.7)

and

w1 =

s1

s1 + s2

and w2 =

s2

s1 + s2

(s1 + s2 > 0) (4.8)

respectively if

x1

def = (uˆx ̄−)

n

j and x2

def = (uˆx ̄+)

n

j

(4.9)

w1

def = (w−)

n

j

and w2

def = (w+)

n

j

(4.10)

NASA/TM—2005-213868 13
s1

def = |(uˆx ̄+)

n

j

|

α

and s2

def = |(uˆx ̄−)

n

j

|

α

(α ≥ 0) (4.11)

Equation. (4.7) represents an weighted average of only two values x1 and x2. However, for the sake of

generality, weighted averages of two or more values will be considered in the following development.

To proceed, let (i) N be an integer ≥ 2, (ii) s`, ` = 1, 2, . . . , N, be given positive numbers, and (iii)

w`

def =

s`

S

, ` = 1, 2, . . . , N (4.12)

where

S

def =

 X

N

`=1

s`

!

> 0 (4.13)

(Note: to streamline the following development, here we assume that s` > 0, ` = 1, 2, . . . N, instead of s` ≥ 0,

` = 1, 2, . . . , N, as could be inferred from Eq. (4.11). However, without causing any practical harm, one can

add a very small positive number, such as 10−20, to each member of a set of nonnegative numbers and turn

all of them into positive numbers). It follows from Eqs. (4.12) and (4.13) that

X

N

`=1

w` = 1, and 1 > w` > 0, ` = 1, 2, . . . , N (4.14)

As such,

W

def =

X

N

`=1

w` x` (4.15)

is an “interpolated” weighted average of the real numbers x`. Note that, unless specified otherwise, hereafter

` = 1, 2, . . . , N is assumed.

Let

δ`

def = w` −

1

N

(4.16)

Then

w` =

1

N

+ δ` (4.17)

Also, with the aid of Eq. (4.14), Eq. (4.16) implies that

X

N

`=1

δ` = 0 (4.18)

Note that W becomes the simple average of x`, ` = 1, 2, . . . , N, if all δ` = 0. Thus the set {δ1, δ2, . . . , δN }

provides a measure of how far the weighted average is deviated from the simple average. In the following, a

simple way to adjust this deviation will be introduced.

Let

δmin

def = min{δ`} and δmax

def = max{δ`} (4.19)

Then Eq. (4.17) and the fact that 1 > w` > 0 for all ` imply that

1 >

1

N

+ δmax and 1

N

+ δmin > 0 (4.20)

Let some δ` 6= 0 (i.e., the case with all w` = 1/N is excluded). Then Eq. (4.18) implies that δmax > 0 > δmin.

The last inequality and Eq. (4.20) can be combined to yield

1 −

1

N

> δmax > 0 > δmin > −

1

N

(some δ` 6= 0) (4.21)

NASA/TM—2005-213868 14
Note that an immediate result of Eq. (4.21) is

σmax

def = min 

1

δmax

(1 −

1

N

), −

1

Nδmin



> 1

(some δ` 6= 0) (4.22)

Given any adjustable real parameter σ > 0, let

δ

0

`

def = σ δ` (4.23)

Then Eq. (4.18) implies that

X

N

`=1

δ

0

` = 0 (4.24)

In turn Eq. (4.24) and

w

0

`

def =

1

N

+ δ

0

`

(4.25)

imply that

X

N

`=1

w

0

` = 1 (4.26)

As such, w

0

`

, ` = 1, 2, . . . , N, form a new set of weight factors. From Eqs. (4.23) and (4.25) one also concludes

that the disparity of the weight factors (i.e., the deviation of the values of the weight factors from 1/N) will

be amplified (reduced) if σ > 1 (σ < 1).

The condition that

1 ≥ w

0

` ≥ 0 (4.27)

will be imposed in the current development. With the aid of Eqs. (4.19), (4.23), and (4.25), and the original

assumption that σ > 0, one concludes that Eq. (4.27) is true if and only if

1 ≥

1

N

+ σδmax and 1

N

+ σδmin ≥ 0 (4.28)

For the special case that in which all δ` = 0, one has δmin = δmax = 0 and thus Eq. (4.28) is true always.

On the other hand, for the case that some δ` 6= 0, Eqs. (4.21) and (4.22) imply that Eq. (4.28) is equivalent

to

σmax ≥ σ (some δ` 6= 0) (4.29)

Let some δ` 6= 0. Then, according to Eq. (4.22), σmax > 1. Moreover, σmax increases as |δmax| and

|δmin| decrease. In fact, σmax → +∞ as δmax → 0

+ and δmin → 0−. Thus the range of the values of σ

allowed becomes larger when |δmax| and |δmin| become smaller. Note that, when W defined in Eq. (4.15)

almost becomes a simple average (i.e., when |δmax|  1 and |δmin|  1), the disparity of the weight factors

must be sharply amplified such that the weight average

W0 def =

X

N

`=1

w

0

` x` (4.30)

will deviate substantially from the simple average. In this case, the larger range of the values of σ allowed

meets the need to use a larger value of σ. In practice, the value of σ used can be generated using a preset

formula as long as the generated value is less than or equal to σmax. For the case that the value generated

using the preset formula is larger than σmax, σ = σmax is assumed.

NASA/TM—2005-213868 15
As an example, consider the N = 2 case in which x` and s`, ` = 1, 2, are defined by Eqs. (4.9) and

(4.11). It was explained earlier that, for this case, w1 → 1/2 and w2 → 1/2 as ν → 0. In other words,

the weighted average w1x1 + w2x2 almost becomes a simple average when |ν|  1. To prevent this from

happening, the weight factors w1 and w2, respectively, are replaced by the new weighted factors w

0

1

and w

0

2

generated assuming

σ = min σmax,

σo

|ν|



(σo > 0) (4.31)

where σo > 0 is a preset parameter in the order of 1. According to Eqs. (4.16), (4.19), and (4.22), the fact

that w1 → 1/2 and w2 → 1/2 as ν → 0 implies that σmax → ∞ as ν → 0. As such Eq. (4.31) implies that

σ  1 when |ν|  1.

Note that, for any N = 2 case, one of δ1 and δ2 is δmax while another is δmin. As a result, Eqs. (4.18),

(4.20), and (4.22) imply that

0 < δmax = −δmin < 1/2 (some δ` 6= 0) (4.32)

and

σmax =

1

2 δmax

(some δ` 6= 0) (4.33)

Also for any case with N = 2, δmax > 0 and σ = σmax, Eqs. (4.23), (4.25), (4.32), and (4.33) imply that:

(i) w

0

1 = 1 and w

0

2 = 0 if δ1 = δmax, and (ii) w

0

2 = 1 and w

0

1 = 0 if δ2 = δmax.

This completes the description of a new approach by which the weight factors w

0

`

, ` = 1, 2, . . . , N, are

generated from the given weight factors w`, ` = 1, 2, . . . , N. In the following, Another approach will be

described.

To proceed, the indices of s`, ` = 1, 2, . . . , N, will be reshuffled such that

sN ≥ sN−1 ≥ . . . ≥ s1 > 0 (4.34)

As such, Eqs. (4.12) and (4.13) imply that Eq. (4.14) can be replaced by a set of stronger conditions, i.e.,

X

N

`=1

w` = 1 and 1 > wN ≥ wN−1 ≥ . . . ≥ w1 > 0 (4.35)

Next let

η`

def =

s`+1

s`

− 1, ` = 1, . . . , N −1 (4.36)

Then (i) η` ≥ 0, ` = 1, . . . , N −1, and (ii)

s`+1 =

"Y

`

`

0=1

(1 + η`

0 )

#

s1, ` = 1, . . . , N −1 (4.37)

Given any adjustable real parameter σ > 0, let (i) s ̃1 = s1 and

s ̃`+1 =

"Y

`

`

0=1

(1 + ση`

0 )

#

s ̃1, ` = 1, . . . , N −1 (4.38)

and (ii)

w ̃`

def =

s ̃`

S ̃

, ` = 1, 2, . . . , N (4.39)

NASA/TM—2005-213868 16
where

S ̃

def =

 X

N

`=1

s ̃`

!

> 0 (4.40)

Because σ > 0 and η` ≥ 0, ` = 1, . . . , N −1, Eq. (4.38) implies that

s ̃N ≥ s ̃N−1 ≥ . . . ≥ s ̃1 > 0 (4.41)

Also, as a result of Eqs. (4.39)–(4.41), one has

X

N

`=1

w ̃` = 1 and 1 > w ̃N ≥ w ̃N−1 ≥ . . . ≥ w ̃1 > 0 (4.42)

As such, w ̃`, ` = 1, 2, . . . , N, form a new set of weight factors and

W ̃

def =

X

N

`=1

w ̃` x` (4.43)

is an “interpolated” weighted average of the real numbers x`. Note that, for the special case that sN = sN−1 =

. . . = s1 > 0, it is easy to see that (i) w` = w ̃` = 1/N, ` = 1, 2, . . . , N, and (ii) η` = 0, ` = 1, . . . , N −1.

Let `1 and `2 be any pair of integers with 1 ≤ `1 < `2 ≤ N. Then Eqs. (4.12) and (4.37)–(4.39) imply

that

w`2

w`1

=

`

Y2−1

`

0=`1

(1 + η`

0 ) (4.44)

and

w ̃`2

w ̃`1

=

`

Y2−1

`

0=`1

(1 + ση`

0 ) (4.45)

Because σ > 0 and η` ≥ 0, ` = 1, . . . , N −1, a comparison of Eqs. (4.44) and (4.45) reveals that w`2 /w`1 =

w ̃`2 /w ̃`1 = 1 if η` = 0 for all ` with `1 ≤ ` ≤ (`2 − 1). However, in case that η` 6= 0 for at least one ` with

`1 ≤ ` ≤ (`2 − 1), one has

w ̃`2

w ̃`1







>

w`2

w`1

if σ > 1

<

w`2

w`1

if σ < 1

=

w`2

w`1

if σ = 1

(4.46)

¿From the above discussions, one concludes that, except for the special case in which sN = sN−1 = . . . =

s1, the disparity of w ̃` is greater (less) than that of w` if σ > 1 (σ < 1). Note that the current approach

for amplifying the weight factors has one advantage over the approach described earlier, i.e., in the current

approach, there is no upper bound for the value of σ one could use. Thus, in the current approach, Eq. (4.31)

can be simplified as

σ =

σo

|ν|

(4.47)

where σo > 0 again is a preset parameter in the order of 1.

Note that, after sorting through the differences in the notations used in [49] and here, and using the

fact that the value of the smaller of the parameters (s−)

n

j

and (s+)

n

j defined in Eq. (3.23) of [50] is zero, it

can be shown that the weighted-averaging technique introduced in Eqs. (3.23), (3.26), and (3.27) of [50] is

equivalent to the special case N = 2 and σo = 1/2 of the second approach just described above.

NASA/TM—2005-213868 17
Consider the N = 2 case in which x` and s`, ` = 1, 2, are defined by Eqs. (4.9) and (4.11). Let (w

0

−)

n

j

and (w

0

+)

n

j

, respectively, be the weight factors associated with (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j

generated using the first

approach described in Sec. 4.3. Then, by definition, Scheme w-3 is formed by Eq. (2.8) and

(ux ̄)

n

j = (w

0

−)

n

j

(uˆx ̄−)

n

j + (w

0

+)

n

j

(uˆx ̄+)

n

j

(4.48)

On the other hand, let (w ̃−)

n

j

and (w ̃+)

n

j

, respectively, be the weight factors associated with (uˆx ̄−)

n

j

and (uˆx ̄+)

n

j

generated using the second approach described in Sec. 4.3. Then, by definition, Scheme w-4 is

formed by Eq. (2.8) and

(ux ̄)

n

j = (w ̃−)

n

j

(uˆx ̄−)

n

j + (w ̃+)

n

j

(uˆx ̄+)

n

j

(4.49)

5. Courant number and Mach number insensitive solvers

We consider a dimensionless form of the 1-D unsteady Euler equations of a perfect gas. Let ρ, v, p, and

γ be the mass density, velocity, static pressure, and constant specific heat ratio, respectively. Let

u1 = ρ, u2 = ρv, u3 = p/(γ − 1) + (1/2)ρv2(5.1)

f1 = u2 (5.2)

f2 = (γ − 1)u3 + (1/2)(3 − γ)(u2)

2

/u1 (5.3)

and

f3 = γu2u3/u1 − (1/2)(γ − 1)(u2)

3

/(u1)

2

(5.4)

Then the Euler equations can be expressed as

∂um

∂t

+

∂fm

∂x

= 0, m = 1, 2, 3 (5.5)

The integral form of Eq. (5.5) in space-time E2 is

I

S(V )

~hm · d~s = 0, m = 1, 2, 3 (5.6)

where ~hm = (fm, um), m = 1, 2, 3, are the space-time mass, momentum, and energy current density vectors,

respectively.

As a preliminary, let

fm,k

def = ∂fm/∂uk, m, k = 1, 2, 3 (5.7)

Let F be the Jocobian matrix formed by fm,k, m, k = 1, 2, 3. Then Eqs. (5.2)–(5.4) imply that

F =





0 1 0

γ − 3

2



u2

u1

2

(3 − γ)

u2

u1

γ − 1

(γ − 1) 

u2

u1

3

− γ

u2u3

(u1)

2

γ

u3

u1

−

3(γ − 1)

2



u2

u1

2

γ

u2

u1





(5.8)

Let c be the sonic speed. Then

c =

vuutγ(γ − 1) "

u3

u1

−

1

2



u2

u1

2

#

(5.9)

NASA/TM—2005-213868 18

4.4. Schemes w-3 and w-4
Let G be the 3 × 3 matrix defined by

G

def =





1

u1

√

2c

u1

√

2c

u2

u1

u2

√

2c

−

u1

√

2

u2

√

2c

+

u1

√

2

1

2



u2

u1

2

(u2)

2

2

√

2cu1

−

u2

√

2

+

u1c

√

2(γ − 1)

(u2)

2

2

√

2cu1

+

u2

√

2

+

u1c

√

2(γ − 1)





(5.10)

Then the inverse of G is given by

G

−1 =





1 +

1 − γ

2c

2



u2

u1

2

(γ − 1)u2

c

2u1

1 − γ

c

2

(γ − 1)(u2)

2

2

√

2c(u1)

3

+

u2

√

2(u1)

2

(1 − γ)u2

√

2c(u1)

2

−

1

√

2u1

γ − 1

√

2cu1

(γ − 1)(u2)

2

2

√

2c(u1)

3

−

u2

√

2(u1)

2

(1 − γ)u2

√

2c(u1)

2

+

1

√

2u1

γ − 1

√

2cu1





(5.11)

Moreover, for any numbers a1, a2, . . ., an, let diag(a1, a2, . . . , an) denote the diagonal matrix with a1, a2,

. . ., an being the diagonal elements on the first, second, . . ., and n-th rows, respectively. Then, by using

Eqs. (5.8)–(5.11) and v = u2/u1, one has

G

−1 F G = diag(v, v − c, v + c) (5.12)

For any (j, n) ∈ Ω and any (x,t) ∈ SE(j, n), um(x,t), fm(x,t), and ~hm(x,t) are approximated by

u

∗

m(x,t; j, n), f

∗

m(x,t; j, n), and ~h

∗

m(x,t; j, n), respectively. They will be defined shortly. Let

u

∗

m(x,t; j, n)

def = (um)

n

j + (umx)

n

j

(x − xj ) + (umt)

n

j

(t − t

n

) (5.13)

where (um)

n

j

, (umx)

n

j

, and (umt)

n

j

are constants in SE(j, n).

By definition, for each m = 1, 2, 3 and each k = 1, 2, 3, fm and fm,k are functions of u1, u2 and u3. Let

(fm)

n

j

and (fm,k)

n

j denote the values of fm and fm,k, respectively, when the independent variables u1, u2,

and u3, respectively, assume the values of (u1)

n

j

, (u2)

n

j

, and (u3)

n

j

. Let

(fmx)

n

j

def =

X

3

k=1

(fm,k)

n

j

(ukx)

n

j

(5.14)

and

(fmt)

n

j

def =

X

3

k=1

(fm,k)

n

j

(ukt)

n

j

(5.15)

Because

∂fm

∂x

=

X

3

k=1

fm,k

∂uk

∂x

(5.16)

and

∂fm

∂t

=

X

3

k=1

fm,k

∂uk

∂t

(5.17)

NASA/TM—2005-213868 19
(fmx)

n

j

and (fmt)

n

j

can be considered as the numerical analogues of the values of ∂fm/∂x and ∂fm/∂t at

(xj ,t

n), respectively. As a result, we assume that

f

∗

m(x,t; j, n)

def = (fm)

n

j + (fmx)

n

j

(x − xj ) + (fmt)

n

j

(t − t

n

) (5.18)

Because ~hm = (fm, um), we also assume that

~h

∗

m(x,t; j, n)

def = (f

∗

m(x,t; j, n), u

∗

m(x,t; j, n)) (5.19)

Note that, by their definitions, (i) (fm)

n

j

and (fm,k)

n

j

, m = 1, 2, 3, are functions of (um)

n

j

, m = 1, 2, 3, (ii)

(fmx)

n

j

, m = 1, 2, 3, are functions of (um)

n

j

and (umx)

n

j

, m = 1, 2, 3, and (iii) (fmt)

n

j

are functions of (um)

n

j

and (umt)

n

j

, m = 1, 2, 3.

Moreover we assume that, for any (x,t) ∈ SE(j, n), um = u

∗

m(x,t; j, n) and fm = f

∗

m(x,t; j, n) satisfy

Eq. (5.5), i.e.,

∂u

∗

m(x,t; j, n)

∂t

+

∂f

∗

m(x,t; j, n)

∂x

= 0 (5.20)

According to Eqs. (5.13) and (5.18), Eq. (5.20) is equivalent to

(umt)

n

j = −(fmx)

n

j

(5.21)

Because (fmx)

n

j

are functions of (um)

n

j

and (umx)

n

j

, Eq. (5.21) implies that (umt)

n

j

are also functions of

(um)

n

j

and (umx)

n

j

. From this result and the facts stated following Eq. (5.19), one concludes that the only

independent discrete variables needed to be solved in the current marching scheme are (um)

n

j

and (umx)

n

j

.

In the current development, the Euler counterpart of Eq. (2.11), i.e.,

I

S(CE(j,n))

~h

∗

m · d~s = 0 (5.22)

is assumed. For any (j, n) ∈ Ω, let

(umx ̄)

n

j

def =

∆x

4

(umx)

n

j

(5.23)

and

(sm)

n

j

def = (umx ̄)

n

j +

∆t

∆x

(fm)

n

j +

(∆t)

2

4∆x

(fmt)

n

j

(5.24)

Then, with the aid of Eqs. (5.13)–(5.15), (5.18), (5.19), and (5.21), Eq. (5.22) implies that

(um)

n

j =

1

2

h

(um)

n−1/2

j−1/2 + (um)

n−1/2

j+1/2 + (sm)

n−1/2

j−1/2 − (sm)

n−1/2

j+1/2

i

(5.25)

Eq. (5.25) forms the first component of each of the Euler schemes to be constructed here. As will be shown,

the second component which evaluates (umx ̄)

n

j

is scheme dependent.

5.1. The Euler c-τ scheme

To proceed, consider any (j, n) ∈ Ω. Let ~u, ~u

n

j

, (~ux)

n

j

, (~ux ̄)

n

j

, and (~ut)

n

j

, respectively, denote the column

matrices formed by um, (um)

n

j

, (umx)

n

j

, (umx ̄)

n

j

, (umt)

n

j

, m = 1, 2, 3. Moreover, recall that the fluid velocity

v, the sonic speed c and the matrices F, G, and G−1 are functions of um, m = 1, 2, 3. Thus we will define

v

n

j

, c

n

j

, F

n

j

, Gn

j

, and (G−1)

n

j

, respectively, to be the values of v, c, F, G, and G−1 when the independent

variables um, m = 1, 2, 3, respectively, assume the values of (um)

n

j

, m = 1, 2, 3. Given the above definitions,

let

(ν1)

n

j

def =

v

n

j ∆t

∆x

, (ν2)

n

j

def =

(v − c)

n

j ∆t

∆x

, and (ν3)

n

j

def =

(v + c)

n

j ∆t

∆x

(5.26)

NASA/TM—2005-213868 20
and

(νmax)

n

j

def = (|v

n

j

| + |c

n

j

|)

∆t

∆x

(5.27)

Because

(νmax)

n

j = max{|(ν1)

n

j

|, |(ν2)

n

j

|, |(ν3)

n

j

|} (5.28)

(νmax)

n

j

can be interpreted as the Courant number at the mesh point (j, n).

Next, for each m = 1, 2, 3, let the points P

+

m and P

−

m, and the parameter (τm)

n

j

shown in Fig. 5 be

defined in the exact same manner by which the points P

+ and P−, and the parameter τ were defined (see

Fig. 3). Moreover, let (i) [

~b ]m denote the mth component of any column matrix ~b, (ii)



(G

−1

)

n

j ~u



m

(P

+

m)

def =

h

(G

−1

)

n

j ~u

n−1/2

j+1/2

i

m

+ (∆t/2) h(G

−1

)

n

j

(~ut)

n−1/2

j+1/2

i

m

−



1 − (τm)

n

j



h

(G

−1

)

n

j

(~ux ̄)

n−1/2

j+1/2

i

m

(5.29)

and (iii)



(G

−1

)

n

j

~u



m

(P

−

m)

def =

h

(G

−1

)

n

j

~u

n−1/2

j−1/2

i

m

+ (∆t/2) h(G

−1

)

n

j

(~ut)

n−1/2

j−1/2

i

m

+



1 − (τm)

n

j



h

(G

−1

)

n

j

(~ux ̄)

n−1/2

j−1/2

i

m

(5.30)

Note that, in the current development of the procedures for evaluating (~ux ̄)

n

j

, (G−1)

n

j

is treated as a fixed

constant square matrix for any given fixed (j, n) ∈ Ω while ~u is treated as a variable column matrix (this

practice, in spirit, is similar to the definition of u

∗

(x,t; j, n) given in Eq. (2.3) where u

n

j

, (ux)

n

j

, and (ut)

n

j

are treated as constants while x and t are treated as variables). Thus

∂



(G−1)

n

j

~u



m

∂t

=



(G

−1

)

n

j

∂~u

∂t



m

and

∂



(G−1)

n

j

~u



m

∂x ̄

=



(G

−1

)

n

j

∂~u

∂x ̄



m

(5.31)

It follows that h(G−1)

n

j

(~ut)

n−1/2

j±1/2

i

m

and h(G−1)

n

j

(~ux ̄)

n−1/2

j±1/2

i

m

, respectively, are the numerical analogues of

∂



(G−1)

n

j

~u



m

/∂t and ∂



(G−1)

n

j

~u



m

/∂x ̄ at the mesh point (j ± 1/2, n − 1/2). With the aid of this interpretation, Eqs. (5.23), (5.29), and (5.30) imply that 

(G−1)

n

j

~u



m

(P

+

m) is a first-order Taylor’s approximation of



(G−1)

n

j

~u



m

at point P

+

m evaluated using the marching variables at the mesh point (j + 1/2, n − 1/2) while



(G−1)

n

j

~u



m

(P

−

m) is a first-order Taylor’s approximation of 

(G−1)

n

j

~u



m

at point P

−

m evaluated using the

marching variables at the mesh point (j −1/2, n−1/2). As such, Eqs. (5.29) and (5.30) can be considered as

the Euler versions of Eqs. (3.4) and (3.5), respectively. In the following, we will construct the Euler versions

of Eqs. (3.6)–(3.9).

By using Eqs. (5.14), (5.21), and (5.23), one has

(∆t/2)(G

−1

)

n

j

(~ut)

n−1/2

j±1/2 = −(2∆t/∆x)(G

−1

)

n

j F

n−1/2

j±1/2 G

n

j

(G

−1

)

n

j

(~ux ̄)

n−1/2

j±1/2

(5.32)

Let F

n−1/2

j±1/2 ≈ F

n

j

. Then Eqs. (5.12) and (5.26) imply that

∆t

∆x

(G

−1

)

n

j F

n−1/2

j±1/2 G

n

j ≈

∆t

∆x

(G

−1

)

n

j F

n

j G

n

j = diag((ν1)

n

j

,(ν2)

n

j

,(ν3)

n

j

) (5.33)

Combining Eqs. (5.32) and (5.33), one has

(∆t/2) h(G

−1

)

n

j

(~ut)

n−1/2

j±1/2

i

m

≈ −2(νm)

n

j

h

(G

−1

)

n

j

(~ux ̄)

n−1/2

j±1/2

i

m

(5.34)

Substituting Eq. (5.34) into Eqs. (5.29) and (5.30), one arrives at the current Euler versions of Eqs. (3.6)

and (3.7), i.e.,



(G

−1

)

n

j

~u



m

(P

+

m) ≈

h

(G

−1

)

n

j

~u

n−1/2

j+1/2

i

m

−



2(νm)

n

j + 1 − (τm)

n

j



h

(G

−1

)

n

j

(~ux ̄)

n−1/2

j+1/2

i

m

(5.35)

NASA/TM—2005-213868 21
and



(G

−1

)

n

j ~u



m

(P

−

m ) ≈

h

(G

−1

)

n

j ~u

n−1/2

j−1/2

i

m

−



2(νm)

n

j − 1 + (τm)

n

j



h

(G

−1

)

n

j

(~ux ̄)

n−1/2

j−1/2

i

m

(5.36)

Next we introduce the Euler version of Eq. (3.8), i.e.,



(G

−1

)

n

j

~u



mx ̄

(j, n)

def =

∆x

4



(G−1)

n

j

~u



m

(P

+

m) −



(G−1)

n

j

~u



m

(P

−

m)



1 + (τm)

n

j



∆x/2



(τm)

n

j

6= −1



, m = 1, 2, 3 (5.37)

where (G−1)

n

j

~u



m

(P

+

m) and 

(G−1)

n

j

~u



m

(P

−

m), respectively, are to be evaluated using Eqs. (5.29) and (5.30)

instead of the approximated formulae Eqs. (5.35) and (5.36). With the aid of Fig. 5, Eq. (5.37) implies that



(G−1)

n

j

~u



mx ̄

(j, n) is a central-difference approximation of ∂



(G−1)

n

j

~u



m

/∂x ̄ at the mesh point (j, n).

To proceed, note that

∂



(G−1)

n

j

~u



∂x ̄

= (G

−1

)

n

j

∂~u

∂x ̄

(5.38)

Thus

∂~u

∂x ̄

= G

n

j

∂



(G−1)

n

j

~u



∂x ̄

(5.39)

Let (G−1)

n

j

~u



x ̄

(j, n) be the column matrix formed by



(G−1)

n

j

~u



mx ̄

(j, n), m = 1, 2, 3 (which are defined in

Eq. (5.37)). Then, with the aid of Eq. (5.39), one concludes that

(

~uˆx ̄)

n

j

def = G

n

j

n

(G

−1

)

n

j ~u



x ̄

(j, n)

o

(5.40)

represents a central difference approximation of ∂~u/∂x ̄ at the mesh point (j, n). Thus

(~ux ̄)

n

j = (

~uˆx ̄)

n

j

(5.41)

is an Euler version of Eq. (3.9). The Euler c-τ scheme is defined by Eqs. (5.25), (5.29), (5.30), (5.37), (5.40)

and (5.41). It has been shown by numerical experiments that stability of this scheme generally requires that

(νmax)

n

j < 1; and (τm)

n

j ≥ τo



((νm)

n

j

)

2



((j, n) ∈ Ω), m = 1, 2, 3 (5.42)

Obviously, the form of Eq. (5.42) is very similar to that of Eq. (3.12).

Consider the special case in which (τ1)

n

j = (τ2)

n

j = (τ3)

n

j = τ

n

j

. For this case, (i) points P

+

1

, P

+

2

and P

+

3

coincide (the resulting common point will be denoted by P

+); and (ii) points P

−

1

, P

−

2

, and P

−

3

also coincide

(the resulring common point will be denoted by P

−). As such it can be shown easily that Eqs. (5.29) and

(5.30) reduce to

~u(P

+) =



~u + (∆t/2)~ut − (1 − τ

n

j

)~ux ̄

n−1/2

j+1/2

(5.43)

and

~u(P

−) =



~u + (∆t/2)~ut + (1 − τ

n

j

)~ux ̄

n−1/2

j−1/2

(5.44)

respectively. Also Eqs. (5.37) and (5.40) can be used to conclude that

(

~uˆx ̄)

n

j

def =

~u(P

+) − ~u(P−)

2(1 + τ

n

j

)

(τ

n

j

6= −1) (5.45)

Note that: (i) with the aid of Eq. (5.28) and the fact that τo(s) is strictly monotonically increasing in the

interval 0 < s < 1, one concludes that Eq. (5.42) reduces to

(νmax)

n

j < 1; and τ

n

j ≥ τo



((νmax)

n

j

)

2



((j, n) ∈ Ω) (5.46)

NASA/TM—2005-213868 22
for this special case; and (ii) Eqs. (3.4), (3.5), and (3.8), respectively, are very similar to their Euler versions

Eqs. (5.43)–(5.45).

The above simplified scheme has the advantage in its structural simplicity. However, it suffers from the

loss of the capability to adjust the value of individual (τm)

n

j

. This is a serious disadvantage because the

theory presented in Sec. 3 along with the fact that Eqs. (3.6) and (3.7), respectively, are very much similar

to Eqs. (5.35) and (5.36) in algebraic form strongly suggests that the Euler c-τ scheme will perform better

if, for each m, the value of (τm)

n

j

can be adjusted to match the value of τo



((νm)

n

j

)

2



. This is particularly

true if, for some (j, n) ∈ Ω, there is a large disparity among the values of |(νm)

n

j

|, m = 1, 2, 3, i.e.,

(νmax)

n

j  min{|(ν1)

n

j

|, |(ν2)

n

j

|, |(ν3)

n

j

|} (5.47)

As an example, consider an Euler flow solution with an extremely small Mach number everywhere (i.e.,

|v

n

j

|  |c

n

j

|, (j, n) ∈ Ω), For such a case (see Eqs. (5.26)–(5.28)),

|(ν1)

n

j

|  (νmax)

n

j ≈ |(ν2)

n

j

| ≈ |(ν3)

n

j

| ((j, n) ∈ Ω) (5.48)

As such, Eqs. (3.13) and (5.46) can be used to show that

τo



((ν1)

n

j

)

2



 τ

n

j

(5.49)

Because the c-τ scheme tends to become more dissipative as the value of τ increases from τo(ν

2

) (see comment

(b) given following Eq. (3.14)), the fact that Eqs. (5.35) and (5.36), respectively, are very much similar to

Eqs. (3.6) and (3.7) in algebraic form strongly suggests that the component



(G−1)

n

j

~u



1

will become highly

dissipative under the condition Eq. (5.49). As such, one expects that a solution to the simplfied Euler c-τ

scheme would become highly dissipative at a region with a small Mach number.

5.2. The Euler c-τ

∗

schemes

Let hm(s) (0 ≤ s < 1), m = 1, 2, 3, be strictly monotonically increasing smooth functions which satisfy

the current version of Eq. (3.15), i.e.,

hm(0) = 0; lim

s→1−

hm(s) = 1; and hm(s) ≥ τo(s) if 0 < s < 1, m = 1, 2, 3 (5.50)

Then, based on the similarity in form that existed between Eqs. (3.6)–(3.8) and Eqs. (5.35)–(5.37), an ideal

Euler c-τ

∗

scheme can be formed as a special Euler c-τ scheme in which

(τm)

n

j = hm



(νm)

n

j

2



(|(νm)

n

j

| < 1), m = 1, 2, 3 (5.51)

Obviously, Eq. (5.51) is the Euler version of Eq. (3.16).

For the same reason that justifies the use of the relation Eq. (3.22), a special Euler c-τ

∗

scheme can be

defined as a special Euler c-τ scheme in which

(τm)

n

j = (βm)

n

j

|(νm)

n

j

| ((βm)

n

j ≥ 1; |(νm)

n

j

| < 1), m = 1, 2, 3 (5.52)

Here (βm)

n

j ≥ 1, m = 1, 2, 3, are adjustable parameters which may vary from one mesh point to another.

Corresponding to the simplified Euler c-τ scheme defined by Eqs. (5.25), (5.41), and (5.43)–(5.45), one

can construct an Euler c-τ

∗

scheme in which

τ

n

j = h



(νmax)

n

j

2



(|(νmax)

n

j

| < 1) (5.53)

where h(s) (0 ≤ s < 1) is any strictly monotonically increasing smooth function which satisfies Eq. (3.15).

Obviously, we can also construct another scheme in which

τ

n

j = β

n

j

|(νmax)

n

j

| (β

n

j ≥ 1; |(νmax)

n

j

| < 1) (5.54)

NASA/TM—2005-213868 23
Here β

n

j ≥ 1 is an adjustable parameter which may vary from one mesh point to another.

5.3. Other Euler extensions

Let



(G

−1

)

n

j ~u



mx ̄−

(j, n)

def =

∆x

4



(G−1)

n

j

~u

n

j



m

−



(G−1)

n

j

~u



m

(P

−

m)



1 + (τm)

n

j



∆x/4

(5.55)

and



(G

−1

)

n

j

~u



mx ̄+

(j, n)

def =

∆x

4



(G−1)

n

j

~u

n

j



m

(P

+

m) −



(G−1)

n

j

~u

n

j



m

1 + (τm)

n

j



∆x/4

(5.56)

By definition, (G−1)

n

j

~u



mx ̄−

(j, n) and (G−1)

n

j

~u



mx ̄+

(j, n) are two one-sided difference approximations

of ∂



(G−1)

n

j

~u



m

/∂x ̄ at the mesh point (j, n) with one being evaluated from the left and another from

the right. Let (τm)

n

j be defined using either Eq. (5.51) or Eq. (5.52). Then, for each m, one can define



(G−1)

n

j

~u

w

mx ̄

(j, n) to be an weighted average of (G−1)

n

j

~u



mx ̄−

(j, n) and (G−1)

n

j

~u



mx ̄+

(j, n) using any

weighted-averaging technique described in Sec. 4. As an example, an weighted average constructed using the

second approach described in sec. 4.3 with N = 2 is given by



(G

−1

)

n

j

~u

w

mx ̄

(j, n)

def =

[1 + σmηm−]

n

j



(G−1)

n

j

~u



mx ̄+

(j, n) + [1 + σmηm+]

n

j



(G−1)

n

j

~u



mx ̄−

(j, n)

[2 + σm(ηm− + ηm+)]n

j

(5.57)

Here

(ηm±)

n

j

def =









(G−1)

n

j

~u



mx ̄±

(j, n)







min n







(G−1)

n

j

~u



mx ̄+

(j, n)





 ,









(G−1)

n

j

~u



mx ̄−

(j, n)







o − 1 ≥ 0 (5.58)

and

(σm)

n

j

def =

σo

|(νm)

n

j

|

(5.59)

with σo > 0 being a preset number in the order of 1. Note that: (i) for each (j, n) ∈ Ω, the value of the

smaller of (ηm+)

n

j

and (ηm−)

n

j

is zero; and (ii) to avoid dividing by zero, in practice, a small positive number

such that 10−60 should be added to each of the denominators that appear in Eqs. (5.58) and (5.59).

Let (G−1)

n

j

~u

w

x ̄

(j, n) be the column matrix formed by



(G−1)

n

j

~u

w

mx ̄

(j, n), m = 1, 2, 3 (which are defined

in Eq. (5.57)). Then, with the aid of Eq. (5.39), one concludes that

(~u

w

x ̄

)

n

j

def = G

n

j

n

(G

−1

)

n

j

~u

w

x ̄

(j, n)

o

(5.60)

represents an weighted-averaging approximation of ∂~u/∂x ̄ at the mesh point (j, n). Thus an Euler weightedaveraging c-τ

∗

scheme can be formed by Eq. (5.25) and

(~ux ̄)

n

j = (~u

w

x ̄

)

n

j

(5.61)

For the simplified case in which (τ1)

n

j = (τ2)

n

j = (τ3)

n

j = τ

n

j

, Eqs. (5.55) and (5.56) are replace by

(~ux ̄−)

n

j

def =

~u

n

j − ~u(P

−)

1 + τ

n

j

(5.62)

and

(~ux ̄+)

n

j

def =

~u(P

+) − ~un

j

1 + τ

n

j

(5.63)

NASA/TM—2005-213868 24
respectively. Let (i) τ

n

j be defined using either Eq. (5.53) or Eq. (5.54); and (ii) (umx ̄±)

n

j denotes the mth

component of (~ux ̄±)

n

j

. Then, for each m, an weighted average of (umx ̄−)

n

j

and (umx ̄+)

n

j

constructed using

the second approach described in Sec. 4.3 is given by

(u

w

mx ̄

)

n

j

def =



1 + σ

n

j

(ηm−)

n

j



(umx ̄+)

n

j +



1 + σ

n

j

(ηm+)

n

j



(umx ̄−)

n

j

[2 + σ(ηm− + ηm+)]n

j

(5.64)

Here

(ηm±)

n

j

def =

|(umx ̄±)

n

j

|

min{|(umx ̄+)

n

j

|, |(umx ̄−)

n

j

|} − 1 ≥ 0 (5.65)

and

σ

n

j

def =

σo

(νmax)

n

j

(5.66)

with σo being a preset number in the order of 1. Obviously a simplified Euler weighted-averaging c-τ

∗

scheme

can again be formed using Eqs. (5.25) and (5.61) if (~u

w

x ̄

)

n

j denotes the column matrix formed by (u

w

mx ̄

)

n

j

,

m = 1, 2, 3.

Accuracy of two special Euler c-τ

∗

schemes will be evaluated in Sec. 6. The first is a simplified scheme

in which we assume that

|(νmax)

n

j

| < 1; τ

n

j = β|νmax|; and σ

n

j =

σo

(νmax)

n

j

(j, n) ∈ Ω (5.67)

Note that: (i) β ≥ 1 and σo > 0 are preset numbers in the order of 1; and (ii) the mesh point dependent

parameter β

n

j which appears in Eq. (5.54) is replaced by the mesh point independent parameter β here.

Moreover, for a reason given immediately following Eq. (5.49), one would expect that a solution to such a

simplified scheme becomes highly dissipative at a region with a small Mach number, i.e., the scheme is Mach

number sensitive. Therefore, the first scheme will be referred to as the “Mach number sensitive scheme” in

Sec. 6. Nevertheless this scheme generally is still Courant number insensitive. In fact the so called “new ”

Courant number insensitive solutions presented in Figs. 4–7 of [50] and Figs. 9 and 10 of [52] are generated

using the current simplified scheme with β = 1.0 and σo = 0.5.

In the second scheme to be evaluated in Sec. 6, we assume that

|(νmax)

n

j

| < 1; (τm)

n

j = βm|(νm)

n

j

|; and (σm)

n

j

def =

σo

|(νm)

n

j

|

((j, n) ∈ Ω), m = 1, 2, 3 (5.68)

Note that: (i) βm ≥ 1, m = 1, 2, 3, and σo > 0 are preset numbers in the order of 1; and (ii) the mesh

point dependent parameter (βm)

n

j which appears in Eq. (5.52) is replaced by the mesh point independent

parameter βm here. As will be shown in Sec. 6, not only is it Courant number insensitive, the second scheme

is also Mach number insensitive. As such it will be referred to as the Mach number insensitive scheme in

Sec. 6.

6. Numerical results

In this section, accuracy of the Mach number sensitive and insensitive schemes defined at the end

of Sec. 5 will be evaluated by comparing their numerical solutions with the known analytical solutions of

several shock tube problems. Without exception, (i) the spatial computational domaim which is defined by

−0.505 ≤ x ≤ 0.505, is divided into 101 uniform intervals with ∆x = 0.01; and (ii) the specific heat ratio

γ = 1.4.

At any time t ≥ 0, the exact solutions of a set of shock tube problems considered here are given by

(ρ, p, v) =







(0.125, 1.0,

√

0.14 ξ) if x <

√

0.14 ξt

(10.0, 1.0,

√

0.14 ξ) if x >

√

0.14 ξt

(6.1)

NASA/TM—2005-213868 25
Here, for each problem in the set, ξ is a defining constant parameter. As such, for each problem, (i) p and

v do not vary with x and t; and (ii) there is a contact (density) discontinuity which moves with the velocity

= v. Obviously the discontinuity occurs at x = 0 when t = 0. Moreover, because γ = 1.4, Eq. (6.1) implies

that

c =

rγp

ρ

=







√

11.2

.= 3.3466 if x <

√

0.14 ξt

√

0.14 .= 0.3742 if x >

√

0.14 ξt

(6.2)

and

M

def =

v

c

=







√

0.0125 ξ

.= 0.1118ξ if x <

√

0.14 ξt

ξ if x >

√

0.14 ξt

(6.3)

First consider the case with ξ = 0.01. According to Eqs. (6.1) and (6.2), the flow field is characterized by

extremely small values of M (0.001118 ≤ M ≤ 0.01) and a contact discontinuity moving with an extremely

small velocity (

.= 0.00374). This problem, designated as Problem no. 1, is solved by the Mach number

sensitive scheme assuming ∆t = 0.0024, σo = 2.0, and β = 1.0. It is also solved by the Mach number

insensitive scheme assuming ∆t = 0.0024 and σo = β1 = β2 = β3 = 2.0. Based on Eqs. (6.1) and (6.2), and

the given values of ξ, ∆x, and ∆t, it is estimated that the maximum local Courant number ecountered in

each simulation is 0.804. At t = 60.0 = 25000∆t (i.e., √

0.14 ξt

.= 0.2245), numerical values of pressure and

velocity obtained from both simulations match the constant exact solution values to at least seven significant

digits. As such no graphical comparisons of these numerical variables with their exact solution values are

given. On the other hand, numerical values of density and Mach number at t = 60.0 are compared with

the exact solution values in Fig. 6. It is seen that, for this problem characterized by extremely small values

of M, the contact discontinuity is resolved much more crisply by the numerical values generated using the

Mach number insensitive scheme than those generated using the Mach number sensitive scheme. Note that

this difference in the schemes’ capability to resolve the contact discontinuity could become more pronounced

if the chosen value of β is raised from 1.0 to 2.0, i.e., the value shared by β1, β2, and β3. However, for the

current case with very small values of (ν1)

n

j

, reducing each value of σo and β1 from 2.0 to 1.0 will result in

computational instabilty—obviously, at some (j, n), the assigned values of (σ1)

n

j

and (τ1)

n

j become too small

to sustain stabilty.

Next consider the case with ξ = 10.0. For this case, the flow field is characterized by relatively large

values of M (1.118 ≤ M ≤ 10.0) and a contact discontinuity moving with a relatively large velocity (

.= 3.74).

This problem, designated as Problem no. 2, is solved by the Mach number sensitive scheme assuming

∆t = 0.0012, σo = 2.0, and β = 1.0. It is also solved by the Mach number insensitive scheme assuming

∆t = 0.0012 and σo = β1 = β2 = β3 = 2.0. It is estimated that the maximum local Courant number

ecountered in each simulation is 0.851. At t = 0.06 = 50∆t (i.e., √

0.14 ξt

.= 0.2245), numerical values of

pressure and velocity obtained from both simulations again match the constant exact solution values to at

least seven significant digits. Moreover, as shown in Fig. 7, for this problem with relatively large values of

M, the contact discontinuity is resolved crisply by both schemes. A comparion of these results and those

shown in Fig. 6 reveals that the Mach number sensitive scheme is indeed Mach number sensitive while the

Mach number insensitive scheme is indeed Mach number insensitive.

The last shock tube problem to be considered is Sod’s problem [61]. For this problem, (i) at t = 0,

(ρ, p, v) =







(1.0, 1.0, 0.0) if x < 0

(0.125, 0.1, 0.0) if x > 0

(6.4)

and (ii) 0.0 ≤ M ≤ 0.929 for all x and t ≥ 0. This problem, designated as Problem no. 3, is solved by the

Mach number sensitive scheme assuming ∆t = 4 × 10−6 and σo = β = 1.0. It is also solved by the Mach

number insensitive scheme assuming ∆t = 4 × 10−6 and σo = β1 = β2 = β3 = 1.0. The estimated maximal

local Courant number encounted in each simulation is 0.00088. At t = 0.2 = 50000∆t, the computed solutions

are compared with the exact solution in Fig. 8. In spite of the extremely small maximal local Courant number

NASA/TM—2005-213868 26
encountered, it is seen that the numerical results generated by both schemes match very well with the exact

solution.

Problem no. 3 is also solved by (i) the Mach number sensitive scheme assuming ∆t = 4×10−3 and σo =

β = 1.0; and (ii) the Mach number insensitive scheme assuming ∆t = 4×10−3 and σo = β1 = β2 = β3 = 1.0.

The estimated maximal local Courant number encountered for each simulation is 0.88. At t = 0.2 = 50∆t,

the computed solutions are compared with the exact solution in Fig. 9. It is seen that the numerical solutions

shown in Fig. 8 do not deteriorated much from the current results generated with a much larger maximal

local Courant number. As such, both schemes are indeed Courant number insensitive.

7. Conclusions and discussions

Generally speaking, a stable numerical marching for a non-linear problem requires the presence of

a sufficient amount of numerical dissipation. However, accuracy of the numerical results, especially for

an unsteady problem, will suffer if too much numerical dissipation is present. As such, a careful control

of numerical dissipation is a must for an accurate and stable non-linear unsteady numerical simulation.

However, a proper control of numerical dissipation is a very difficult task. Although one can increase the

numerical dissipation rather easily, it is much harder to reduce it when accuracy consideration requires it.

The CE/SE method is developed from a set of non-dissipative solvers. As such each CE/SE solver is

an extension of a core non-dissipative scheme. It is this unique feature that make it much easier to reduce

numerical dissipation in a CE/SE simulation. It is also the key reason behind the successful construction of

the Courant number and Mach number insensitive Euler solvers described in this paper.

For the 2D and 3D unsteady Euler equations, there are two and three associated Jacobian matrices,

respectively. In general, it is impossible to diagonalize these associated Jacobian matrices simultaneously

using the same diagonalization matrix. Thus extension of the current work to a space of higher dimension

is by no means trivial.

References

1. S.C. Chang and W.M. To, A New Numerical Framework for Solving Conservation Laws–The Method

of Space-Time Conservation Element and Solution Element, NASA TM 104495, August 1991.

2. S.C. Chang and W.M. To, A brief description of a new numerical framework for solving conservation

laws—The method of space-time conservation element and solution element, in Proceedings of the Thirteenth International Conference on Numerical Methods in Fluid Dynamics, Rome, Italy, 1992, edited

by M. Napolitano and F. Sabetta, Lecture Notes in Physics 414, (Springer-Verlag, New York/Berlin,

1992), p. 396.

3. S.C. Chang, The method of space-time conservation element and solution Element—A new approach

for solving the Navier-Stokes and Euler equations, J. Comput. Phys., 119, 295 (1995).

4. S.C. Chang, S.T. Yu, A. Himansu, X.Y. Wang, C.Y. Chow, and C.Y. Loh, The method of space-time

conservation element and solution element—A new paradigm for numerical solution of conservation

laws, in Computational Fluid Dynamics Review 1998 edited by M.M. Hafez and K. Oshima (World

Scientific, Singapore), Vol. 1, p. 206.

5. T. Molls and F. Molls, “Space-Time Conservation Method Applied to Saint Venant Equations,” J. of

Hydraulic Engr., 124(5), 501 (1998).

6. C. Zoppou and S. Roberts, “Space-Time Conservation Method Applied to Saint Venant Equations: A

Discussion,” J. of Hydraulic Engr., 125(8), 891 (1999).

7. S.C. Chang, X.Y. Wang, and C.Y. Chow, The space-time conservation element and solution element

method: A new high-resolution and genuinely multidimensional paradigm for solving conservation laws,”

J. Comput. Phys., 156, 89 (1999).

8. X.Y. Wang, and S.C. Chang, A 2D non-splitting unstructured triangular mesh Euler solver based on the

space-time conservation element and solution element method, Computational Fluid Dynamics Journal,

8(2), 309 (1999).

9. S.C. Chang, X.Y. Wang and W.M. To, Application of the space-time conservation element and solution

element method to one-dimensional convection-diffusion problems, J. Comput. Phys., 165, 189 (2000).

NASA/TM—2005-213868 27
10. J. Qin, S.T. Yu, Z.C. Zhang, and M.C. Lai, Direct Calculations of Cavitating Flows by the Space-Time

CE/SE Method, J. Fuels & Lubricants, SAE Transc., 108(4), 1720 (2000).

11. C.Y. Loh, L.S. Hultgren and S.C. Chang, Wave computation in compressible flow using the space-time

conservation element and solution element method, AIAA J., 39(5), 794 (2001).

12. Z.C. Zhang, S.T. Yu, and S.C. Chang, A Space-Time Conservation Element and Solution Element

Method for Solving the Two- and Three-Dimensional Unsteady Euler Equations Using Quadrilateral

and Hexahedral Meshes, J. Comput. Phys., 175, 168 (2002).

13. K.B.M.Q. Zaman, M.D. Dahl, T.J. Bencic, and C.Y. Loh, Investigation of A ‘Transonic Resonance’

with Convergent-Divergent Nozzles, J. Fluid Mech., 463, 313 (2002).

14. C.Y. Loh and K.B.M.Q. Zaman, Numerical Investigation of ‘Transonic Resonance’ with A ConvergentDivergent Nozzle, AIAA J., 40(12), 2393 (2002).

15. S. Motz, A. Mitrovic, and E.-D. Gilles, Comparison of Numerical Methods for the Simulation of Dispersed Phase Systems, Chemical Engineering Science, 57, 4329 (2002).

16. S. Cioc and T.G. Keith, Application of the CE/SE Method to One-Dimensional Flow in Fluid Film

Bearings, STLE Tribology Transactions, 45, 167 (2002).

17. S. Cioc and T.G. Keith, Application of the CE/SE Method to Two-Dimensional Flow in Fluid Film

Bearings, International J. of Numer. Methods for Heat & Fluid Flow, 13(2), 216 (2003).

18. S. Cioc, F. Dimofte, T.G. Keith, and D.P. Fleming, Computation of Pressurized Gas Bearings Using

the CE/SE Method, STLE Tribology Transactions, 46(1), 128 (2003).

19. S. Cioc, F. Dimofte, and T.G. Keith, Application of the CE/SE Method to Wave Journal Bearings,

STLE Tribology Transactions, 46(2), 179 (2003).

20. A. Ayasoufi and T.G. Keith, Application of the Conservation Element and Solution Element Method

in Numerical Modeling of Heat Conduction with Melting and/or Freezing, International J. of Numer.

Methods for Heat & Fluid Flow, 13(4), 448 (2003).

21. A. Ayasoufi and T.G. Keith, Application of the Conservation Element and Solution Element Method in

Numerical Modeling of Axisymmetric Heat Conduction with Melting and/or Freezing, JSME International J. Series B, 47(1), 115 (2004).

22. A. Ayasoufi and T.G. Keith, Application of the Conservation Element and Solution Element Method in

Numerical Modeling of Three-dimensional Heat Conduction with Melting and/or Freezing, Transactions

of the ASME, J. of Heat Transfer, 126(6), 937 (2004).

23. Y.I. Lim, S.C. Chang, and S.B. Jorgensen, A Novel Partial Differential Algebraic Equation (PDAE)

Solver: Iterative Space-Time Conservation Element/Solution Element (CE/SE) Method, Computers

and Chemical Engineering, 28, 1309 (2004)

24. Y.I. Lim and S.B. Jorgensen, A Fast and Accurate Numerical Method for Solving Simulated Moving

Bed (SMB) Chromatographic Separation Problems, Chemical Engineering Science, 59, 1931 (2004).

25. C.K. Kim, S.T. John Yu, and Z.C. Zhang, Cavity Flow in Scramjet Engine by the Space-Time Conservation Element and Solution Element Method, AIAA J., 42(5), 912 (2004).

26. M. Zhang, S.T. John Yu, S.C. Lin, S.C. Chang, and I. Blankson, Solving Magnetohydrodynamic Equations Without Special Treatment for Divergence-Free Magnetic Field, AIAA J., 42(12), 2605 (2004).

27. K.S. Im, M.C. Lai, S.T. John Yu, and Robert R. Matheson, Jr., Simulation of Spray Transfer Process

in Electrostatic Rotary Bell Sprayer , ASME J. of Fluid Engineering, 126(3), 449 (2004).

28. S.C. Chang, Y. Wu, V. Yang, and X.Y. Wang, Local Time Stepping Procedures for the Space-Time

Conservation Element and Solution Element Method, accepted for publication in International J. of

Computational Fluid Dynamics.

29. X.Y. Wang, C.Y. Chow, and S.C.Chang, Numerical Simulation of Flows Caused by Shock-Body Interaction, AIAA Paper 96-2004 (1996).

30. S.C. Chang, A. Himansu, C.Y. Loh, X.Y., Wang, S.T., Yu, and P. Jorgenson, Robust and Simple

Non-Reflecting Boundary Conditions for the Space-Time Conservation Element and Solution Element

Method, in A Collection of Technical Papers, 13th AIAA CFD Conference, June 29-July 2, 1997,

Snowmass, Colorado, AIAA Paper 97-2077.

NASA/TM—2005-213868 28
31. C.Y. Loh, L.S. Hultgren and S.C. Chang, Vortex Dynamics Simulation in Aeroacoustics by the SpaceTime Conservation Element and Solution Element Method, AIAA Paper 99-0359 (1999)

32. X.Y. Wang, S.C. Chang and P.C.E. Jorgenson, Accuracy Study of the Space-Time CE/SE Method for

Computational Aeroacoustics Problems Involving Shock Waves, AIAA Paper 2000-0474 (2000).

33. C.Y. Loh, L.S. Hultgren, S.C. Chang and P.C.E. Jorgenson, Noise Computation of a Supersonic ShockContaining Axisymmetric Jet by the CE/SE Method, AIAA Paper 2000-0475 (2000).

34. C.Y. Loh, X.Y. Wang, S.C. Chang, and P.C.E. Jorgenson, Computation of Feedback Aeroacoustic System by the CE/SE Method, in Proceedings of the First International Conference on Computational Fluid

Dynamics, Kyoto, Japan, 10-14 July, 2000, edited by N. Satofuka, (Springer-Verlag Berlin Heiderberg

2001), p. 555.

35. C.Y. Loh, L.S. Hultgren and P.C.E. Jorgenson, Near Field Screech Noise Computation for An Underexpanded Supersonic Jet by the CE/SE Method, AIAA Paper 2001-2252 (2001).

36. X.Y. Wang, S.C. Chang, and P.C.E. Jorgenson, Numerical Simulation of Aeroacoustic Field in a 2D

Cascade Involving a Downstream Moving Grid Using the Space-Time CE/SE method, in Proceedings of

the First International Conference on Computational Fluid Dynamics, Kyoto, Japan, 10-14 July, 2000,

edited by N. Satofuka, (Springer-Verlag Berlin Heiderberg 2001), p. 543.

37. X.Y. Wang, S.C. Chang, A. Himansu, and P.C.E. Jorgenson, Gust Acoustic Response of A Single Airfoil

Using the Space-Time CE/SE Method, AIAA Paper 2002-0801 (2002).

38. S.T. Yu and S.C. Chang, Treatments of Stiff Source Terms in Conservation Laws by the Method of

Space-Time Conservation Element and Solution Element, AIAA Paper 97-0435 (1997).

39. S.T. Yu and S.C. Chang, Applications of the Space-Time Conservation Element / Solution Element

Method to Unsteady Chemically Reactive Flows,” AIAA Paper 97-2099, in A Collection of Technical

Papers, 13th AIAA CFD Conference, June 29-July 2, 1997, Snowmass, CO.

40. S.T. Yu, S.C. Chang, P.C.E. Jorgenson, S.J. Park and M.C. Lai, “Treating Stiff Source Terms in Conservation Laws by the Space-Time Conservation Element and Solution Element Method,” in Proceedings

of the 16th International Conference on Numerical Method in Fluid Dynamics, Arcachon, France, 6-10

July, 1998, edited by C.H. Bruneau, (Springer-Verlag Berlin Heidelberg 1998), p. 433.

41. X.Y. Wang and S.C. Chang, A 3D structured/unstructured Euler solver based on the space-time conservation element and solution element method, in A Collection of Technical Papers, 14th AIAA CFD

Conference, June 28–July 1, 1999, Norfolk, Virginia, AIAA Paper 99-3278.

42. N.S. Liu and K.H. Chen, Flux: An Alternative Flow Solver for the National Combustion Code, AIAA

Paper 99-1079.

43. G. Cook, High Accuracy Capture of Curved Shock Front Using the Method of Conservation Element

and Solution Element, AIAA Paper 99-1008.

44. S.C. Chang, Y. Wu, X.Y. Wang, and V. Yang, Local Mesh Refinement in the Space-Time CE/SE

Method, in Proceedings of the First International Conference on Computational Fluid Dynamics, Kyoto,

Japan, 10-14 July, 2000, edited by N. Satofuka, (Springer-Verlag Berlin Heiderberg 2001), p. 61.

45. S.C. Chang, Z.C. Zhang, S.T. John Yu, and P.C.E. Jorgenson, A Unified Wall Boundary Treatment for

Viscous and Inviscid Flows in the CE/SE Method, in Proceedings of the First International Conference

on Computational Fluid Dynamics, Kyoto, Japan, 10-14 July, 2000, edited by N. Satofuka, (SpringerVerlag Berlin Heiderberg 2001), p. 671.

46. Z.C. Zhang, S.T. John Yu, S.C. Chang, and P.C.E. Jorgenson, Calculations of Low-Mach-Number

Viscous Flows without Preconditioning by the Space-Time CE/SE method, in Proceedings of the First

International Conference on Computational Fluid Dynamics, Kyoto, Japan, 10-14 July, 2000, edited by

N. Satofuka, (Springer-Verlag Berlin Heiderberg 2001), p. 127.

47. A. Himansu, P.C.E. Jorgenson, X.Y. Wang, and S.C. Chang, Parallel CE/SE Computational via Domain

Decomposition, in Proceedings of the First International Conference on Computational Fluid Dynamics,

Kyoto, Japan, 10-14 July, 2000, edited by N. Satofuka, (Springer-Verlag Berlin Heiderberg 2001), p.

423.

48. Y. Wu, V. Yang, and S.C. Chang, Space-Time Method for Chemically Reacting Flows with Detailed

Kinetics, in Proceedings of the First International Conference on Computational Fluid Dynamics, Kyoto,

NASA/TM—2005-213868 29
Japan, 10-14 July, 2000, edited by N. Satofuka, (Springer-Verlag Berlin Heiderberg 2001), p. 207.

49. I.S. Chang, Unsteady Rocket Nozzle Flows, AIAA Paper 2002-3884.

50. S.C. Chang, Courant Number Insensitive CE/SE Schemes, AIAA Paper 2002-3890 (2002).

51. I.S. Chang, Unsteday Underexpanded Jet Flows, AIAA Paper 2003-3885.

52. S.C. Chang and X.Y. Wang, Multidimensional Courant Number Insensitive CE/SE Euler Solvers for

Applications Involving Highly Nonuniform Meshes, AIAA Paper 2003-5280.

53. B.S. Venkatachari, G.C. Cheng, and S.C. Chang, Development of A Transient Viscous Flow Solver

Based on Conservation Element-Solution Element Framework, AIAA Paper 2004-3413.

54. B.S. Venkatachari, G.C. Cheng, and S.C. Chang, Courant Number Insensitive Transient Viscous Flow

Solver Based on CE/SE Framework, AIAA Paper 2005-00931.

55. S.C. Chang, Explicit von Neumann Stability Conditions for the c-τ Scheme—A Basic Scheme in the

Development of the CE-SE Courant Number Insensitive Schemes, NASA TM 2005-213627, April 2005.

56. J.C. Yen and D.A. Wagner, Computational Aeroacoustics Using a Simplified Courant Number Insensitive CE/SE Method, AIAA Paper 2005-2820.

57. I.S. Chang, C.L. Chang, and S.C. Chang, Unsteady Navier-Stokes Rocket Nozzle Flows, AIAA Paper

2005-4353

58. Other CE/SE references are posted on: http://www.grc.nasa.gov/www/microbus.

59. B. van Leer, Toward the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical

Convection, J. Comput. Phys., 23, 276 (1977).

60. G.D. van Albada, B. van Leer, and W.W. Robert, A Comparative Study of Computational Methods in

Cosmic Gas Dynamics, Astronom. and Astrophys., 108, 76 (1982).

61. G.A. Sod, A Survey of several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws, J. Comput. Phys., 27, 1 (1978).

NASA/TM—2005-213868 30
x

ds

dr

t

r=(x,t)

V S(V)

NASA/TM—2005-213868 31

Figure 1.—A surface element on the boundary 

S(V) of an arbitrary space-time volume V. 
t

∆x/2 ∆x/2

∆x/2 ∆x/2

∆t/2

n=0

n=1/2

n=1

n=3/2

n=2

j=0

j=1/2

j=1

j=3/2

∆t/2

∆t/2

∆t/2

x

j=–3/2

j=–1

j=–1/2

2(a)—The space-time mesh.

(j,n)

2(b).—SE(j,n)

A

A

B

B

C

C

D

A

D

D

E

E

F

F

(j-1/2,n-1/2) (j+1/2,n-1/2)

(j,n)

(j–1/2,n–1/2) (j+1/2,n–1/2)

(j,n) (j,n)

2(e).—CE(j,n)

2(c).—CE(j,n) 2(d).—CE(j,n)

NASA/TM—2005-213868 32

Figure 2.—The SEs and CEs. 
∆x

2

∆x

4

(1–τ)∆x

4

1

2

(1+τ)∆x

4

∆x

4

P

– M–

∆x

2

∆t

2

∆x

4

(1+τ)∆x

4

(1–τ)∆x

4

∆x

4

P

(j– , n) M+ +

(j, n)

1

2

1

2

(j– , n– )

1

2

(j+ , n)

1

2

1

2

(j+ , n– )

C E

B F

A

D

0.0 0.2 0.4 0.6 0.8 1.0

ν

0.0

0.2

0.4

0.6

0.8

1.0

τ0

(ν

2

)

Figure 4. The function τo(ν

2

).

NASA/TM—2005-213868 33

Figure 3.—Definition of points P−

, M−, M+, and P+. 
∆x

2

[1–(τm)

j

 ]∆x

4

1

2

∆x

4

P

– M–

∆x

2

∆t

2

P (j– , n) M+ +

(j, n)

1

2

1

2

(j– , n– )

1

2

(j+ , n)

1

2

1

2

(j+ , n– )

C E

B F

A

D

m m

n

(τm)

j ∆x

4

n

[1–(τm)

j

 ]∆x

4

∆x

4

n

(τm)

j ∆x

4

n

NASA/TM—2005-213868 34

Figure 5.—Definition of points P−

m, M−, M+, and P+m. 
−0.5 0.0 0.5

x

0.0

2.0

4.0

6.0

8.0

10.0

rho

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

−0.5 0.0 0.5

x

0.000

0.002

0.004

0.006

0.008

0.010

Mach Number

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

Figure 6. Solution comparison for Problem No. 1.

NASA/TM—2005-213868 35
−0.5 0.0 0.5

x

0.0

2.0

4.0

6.0

8.0

10.0

rho

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

−0.5 0.0 0.5

x

0.0

2.0

4.0

6.0

8.0

10.0

Mach Number

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

Figure 7. Solution comparison for Problem No. 2.

NASA/TM—2005-213868 36
−0.5 0.0 0.5

x

0.0

0.2

0.4

0.6

0.8

1.0

rho

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

−0.5 0.0 0.5

x

0.0

0.2

0.4

0.6

0.8

1.0

u

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

−0.5 0.0 0.5

x

0.0

0.2

0.4

0.6

0.8

1.0

p

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

Figure 8. Solution comparison for Problem No. 3 (CFL=0.00088).

NASA/TM—2005-213868 37
−0.5 0.0 0.5

x

0.0

0.2

0.4

0.6

0.8

1.0

rho

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

−0.5 0.0 0.5

x

0.0

0.2

0.4

0.6

0.8

1.0

u

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

−0.5 0.0 0.5

x

0.0

0.2

0.4

0.6

0.8

1.0

p

Mach number sensitive scheme

Mach number insensitive scheme

Exact solution

Figure 9. Solution comparison for Problem No. 3 (CFL=0.88).

NASA/TM—2005-213868 38
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

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

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

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

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

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

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

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

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Technical Memorandum

Unclassified

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field

Cleveland, Ohio 44135–3191

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

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

National Aeronautics and Space Administration

Washington, DC 20546–0001

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

September 2005

NASA TM—2005–213868

AIAA–2005–4355

E–15238

WBS–22–781–30–69

44

Courant Number and Mach Number Insensitive CE/SE Euler Solvers

Sin-Chung Chang

Space-time; The CE/SE method

Unclassified -Unlimited

Subject Categories: 34, 64, and 71

Prepared for the 41st Joint Propulsion Conference and Exhibit cosponsored by the AIAA, ASME, SAE, and ASEE,

Tucson, Arizona, July 10–13, 2005. Responsible person, Sin-Chung Chang, organization code RTS, 216–433–5874.

It has been known that the space-time CE/SE method can be used to obtain 1D, 2D, and 3D steady and unsteady flow

solutions with Mach numbers ranging from 0.0028 to 10. However, it is also known that a CE/SE solution may become

overly dissipative when the Mach number is very small. As an initial attempt to remedy this weakness, new 1D Courant

number and Mach number insensitive CE/SE Euler solvers are developed using several key concepts underlying the

recent successful development of Courant number insensitive CE/SE schemes. Numerical results indicate that the new

solvers are capable of resolving crisply a contact discontinuity embedded in a flow with the maximum Mach number

= 0.01.