Sin-Chung Chang

Glenn Research Center, Cleveland, Ohio

Chau-Lyan Chang

Langley Research Center, Hampton, Virginia

Joseph C. Yen

Jacobs Technology Inc., Tullahoma, Tennessee

Recent Development in the CESE Method for the 

Solution of the Navier-Stokes Equations Using

Unstructured Triangular or Tetrahedral Meshes

With High Aspect Ratio

NASA/TM—2013-218061

September 2013

AIAA–2013–3068
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings 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. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
Sin-Chung Chang

Glenn Research Center, Cleveland, Ohio

Chau-Lyan Chang

Langley Research Center, Hampton, Virginia

Joseph C. Yen

Jacobs Technology Inc., Tullahoma, Tennessee

Recent Development in the CESE Method for the 

Solution of the Navier-Stokes Equations Using

Unstructured Triangular or Tetrahedral Meshes

With High Aspect Ratio

NASA/TM—2013-218061

September 2013

AIAA–2013–3068

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

21st Computational Fluid Dynamics Conference

sponsored by the American Institute of Aeronautics and Astronautics

San Diego, California, June 24–27, 2013
Acknowledgments

Sin-Chung Chang would like to express his deep gratitude to Dr. Grant Cook, a friend and teacher,

for his encouragement during the course of the CESE development.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

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

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 
Abstract

In the multidimensional CESE development, triangles and tetrahedra turn out to be the most natural

building blocks for 2D and 3D spatial meshes. As such the CESE method is compatible with the simplest

unstructured meshes and thus can be easily applied to solve problems with complex geometries. However,

because the method uses space-time staggered stencils, solution decoupling may become a real nuisance in

applications involving unstructured meshes. In this paper we will describe a simple and general remedy

which, according to numerical experiments, has removed any possibility of solution decoupling. Moreover, in

a real-world viscous flow simulation near a solid wall, one often encounters a case where a boundary with high

curvature or sharp corner is surrounded by triangular/tetrahedral meshes of extremely high aspect ratio (up

to 106). For such an extreme case, the spatial projection of a space-time compounded conservation element

constructed using the original CESE design may become highly concave and thus its centroid (referred to as a

spatial solution point) may lie far outside of the spatial projection. It could even be embedded beyond a solid

wall boundary and causes serious numerical difficulties. In this paper we will also present a new procedure

for constructing conservation elements and solution elements which effectively overcomes the difficulties

associated with the original design. Another difficulty issue which was addressed more recently is the wellknown fact that accuracy of gradient computations involving triangular/tetrahedral grids deteriorates rapidly

as the aspect ratio of grid cells increases. The root cause of this difficulty was clearly identified and several

remedies to overcome it were found through a rigorous mathematical analysis. However, because of the

length of the current paper and the complexity of mathematics involved, this new work will be presented in

another paper.

1. Introduction

The space-time conservation element and solution element (CESE) method is a high-resolution and

genuinely multidimensional method for solving conservation laws [1–70]. 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 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

NASA/TM—2013-218061 1

Recent Development in the CESE Method for the Solution of the 

Navier-Stokes Equations Using Unstructured Triangular or 

Tetrahedral Meshes With High Aspect Ratio 

Sin-Chung Chang 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Chau-Lyan Chang 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, Virginia 23681 

Joseph C. Yen 

Jacobs Technology Inc. 

Tullahoma, Tennessee 37388 
any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability to capture shocks

without using Riemann solvers, and also partially contributes to the unusually simple CESE non-reflecting

boundary conditions); (iv) the requirement that each scheme be built from a nondissipative 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 mesh variables to be

solve for simultaneously. Note that CEs are non-overlapping 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 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 constructed from the solution information stored at a space-time solution point. In

general, a CE does not coincide with a SE.

Without using flux-splitting, dimensional-splitting, mesh-alignment or other special techniques, since its

inception in 1991 [1], the unstructured-mesh compatible 2nd-order CESE 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 [15]. The physical phenomena modeled include traveling and interacting shocks, acoustic waves,

vortex shedding, viscous flows, detonation waves, cavitation, flows in fluid film bearings, heat conduction

with melting and/or freezing, MHD vortex, hydraulic jump, crystal growth, chromatographic problems,

solar wind and stress waves in elastic solids [3–70]. In particular, its unexpected simple non-reflecting

boundary conditions [12,50] and rather unique capability to resolve both strong shocks and small disturbances

(e.g., acoustic waves) simultaneously [19,21,22,66,67] make the CESE method an effective tool for attacking

computational aeroacoustics (CAA) problems. Note that the fact that the second-order CESE method is

capable of solving CAA problems accurately is an exception to the commonly-held belief that a second-order

scheme is not adequate for modeling 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 nontraditional features and robust capabilities, the core ideas of the CESE method are

simple. In fact, all of its key features are the results of a pursuit driven by these simple ideas. The first

and foremost is the belief that the method must be solid in physics. As such, in the CESE development,

conservation laws are enforced locally and globally in their natural space-time unity forms for 1D, 2D and 3D

cases. Moreover, because direct physical interaction generally occurs only among the immediate neighbors,

use of the simplest stencil also becomes a common CESE feature. Obviously, this requirement has the effect

of simplifying boundary-condition implementation.

The second idea emerges from the realization that stability and accuracy are two competing issues in

time-accurate computations, i.e., too much numerical dissipation will degrade accuracy while too little of

it will cause instability. In other words, to meet both accuracy and stability requirements, computation

must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a

real dilemma in numerical method development. As an example, high order schemes generally have higher

accuracy and lower numerical dissipation. However, they are susceptible to computational instabilities. In

fact, in complicated real-world applications, not only they seldom live up to their nominal order of accuracy—

generally they possess only first-order accuracy when shocks are present, stability of these schemes often is

highly problem-dependent and difficulty to predict. Moreover, generally speaking, the higher the order of

the scheme, the more restritive its stability conditions become. To confront this issue head-on, in the CESE

development, generally it is required that a solver be built from a nondissipative (i.e., neutrally stable) core

scheme. By definition, computations involving a neutrally stable scheme are performed right on the edge of

instability and therefore the numerical results generated are nondissipative. As such numerical dissipation

can be controlled effectively if the deviation of a solver from its nondissipative core scheme can be adjusted

using some built-in parameters.

Moreover, because an accurate viscous flow simulation requires that the numerical dissipation be much

smaller than the physical dissipation which decreases as the Reynolds number increases, in principle, an

accurate and robust solver for high-Reynolds-number flows must be able to cut numerical dissipation as the

NASA/TM—2013-218061 2
Reynolds number increases. Obviously this requirement can only be met by a solver built from a

nondissipative core solver.

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner;

(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of 1D

schemes so that multidimensional shocks can be captured without using any mesh alignment technique;

(iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the method is

compatible with the simplest unstructured meshes and thus can be easily used to solve problems with

complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible,

and special techniques that has only limited applicability and may cause undesirable side effects be avoided.

Fortunately for the CESE development, realization of the above lesser ideas (i)–(iv) follows naturally from

the first two core ideas.

Building on its initial successes, efforts to refine and improve the CESE method have continued in the

past few years. As an example, it was shown in [5] that the numerical dissipation of a dissipative extension

of a CESE core scheme may increase to an intolerable level as the value of the CF L number decreases from

its maximum stability bound. As such, in a case with a large CF L 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 CF L number can lead to a solution that is highly dissipative in a region where the

local CF L number  1. Even though a remedy was suggested in [5], a simple and robust solution to this

problem had not arrived until a family of new Courant number insensitive schemes [25,31,46–49,51,54,55,58]

was developed with fresh insights. Note that these new schemes have one important advantage, i.e., all

variables at each mesh point can be evaluated explicitly without resorting to solving a system

of linear/nonlinear algebraic equations involving these local mesh variables even in applications

where systems of multidimensional nonlinear PDEs are solved.

As another example, even though they are accurate enough to solve CAA problems, second-order CESE

solvers are not capable of resolving fine flow structures within a boundary layer without using relatively fine

meshes and/or meshes with large aspect ratios. To overcome this limitation, two neutrally stable 4th-order

schemes, referred to as the a(3) and a-4 schemes [59–61], were developed in 2006 and intended to serve

as the core schemes of other high order CESE schemes. Unfortunately, the dissipative extensions of these

new schemes turned out to have the same intractable stability problem that afflicts traditional high order

schemes. With the aid of a conceptual leap, this difficulty was finally overcome through the development

of a new approach [63] by which one can build from a given 2nd-order CESE scheme its 4th-, 6th-, 8th-,...

order versions which have the same stencil and same stability conditions of the 2nd-order scheme,

and also retain all other advantages of the latter scheme. As such, these new high order schemes

can avoid the common shortcomings of traditional high order schemes including: (a) susceptibility to

computational instabilities; (b) computational inefficiency due to their local implicit nature

(i.e., at each mesh points, need to solve a system of linear/nonlinear equations involving all the

mesh variables associated with this mesh point); (c) use of large and elaborate stencils which

complicates boundary treatments and also makes efficient parallel computing much harder;

(d) difficulties in applications involving complex geometries; and (e) use of problem-specific

techniques which are needed to overcome stability problems but often cause undesirable side

effects.

In the current paper, we will describe other recent developments which successfully address some of

the few remaining CESE issues. Some of the new techniques devised have been implemented in the NASA

code ez4d (developed by Dr. Chau-Lyan Chang, the second author) and the Jacobs Technology Inc. code

JUSTUS (developed by Dr. Joseph C. Yen, the third author). Their effectiveness has been verified by

applications involving hypersonic viscous flow over rough solid surfaces [62,64,65] and also by applications

involving studies of the WICS (weapons Internal Carriage and separation) test data produced at Arnold

Engineering Development Center [66,67]. Note that the current paper would represent the first systematic

documentation of the new techniques even though the improved numerical results were presented in [62,64–

67]. To explain clearly these new techniques and their significance, we will begin our discussion by providing

some basic information about the CESE method.

NASA/TM—2013-218061 3
For simplicity, first we review existing 2nd-order CESE solvers for the simple partial differential equation

(PDE)

∂u

∂t + a

∂u

∂x = 0 (1.1)

where the advection speed a = 0 is a constant. Let 6 x1 = x, and x2 = t be the coordinates of a two-dimensional

Euclidean space E2. Then, because Eq. (1.1) can be expressed as ∇ · ~ ~h = 0 with ~h

def = (au, u) (i.e., au and u

are the x1- and x2- components of ~h, respectively), Gauss’ divergence theorem in the space-time E2 implies

that Eq. (1.1) is the differential form of the integral conservation law

I

S(V )

~h · d~s = 0 (1.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 vector 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. (1.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 differential 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 Ω1 denote the set of all space-time staggered mesh points (j, n) 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, . . .. Thus (j + n) is an

half integer for each (j, n) ∈ Ω1. Each (j, n) ∈ Ω1 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. By definition,

the end points of the line segments referred to above are excluded from SE(j, n) so that two SEs will not

overlap.

Eq. (1.2) will be simulated numerically assuming that, for any (x, t) ∈ SE(j, n), 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

) (j, n) ∈ Ω1 (1.3)

and

~h

∗

(x, t; j, n)

def =



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



(x, t) ∈ SE(j, n) and (j, n) ∈ Ω1 (1.4)

Note that: (i) u

n

j

, (ux)

n

j

, and (ut)

n

j

are constants in SE(j, n), and (in a rough sense) they can be considered

as the numerical analogues of the values of u, ∂u/∂x, and ∂u/∂t at the mesh point (j, n), respectively, (ii)

(xj , tn) are the coordinates of the mesh point (j, n) with xj = j∆x and t

n = n∆t, and (iii) Eq. (1.4) is the

numerical analogue of the definition ~h = (au, u).

Let u = u

∗

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

(ut)

n

j = −a (ux)

n

j

(j, n) ∈ Ω1 (1.5)

As a result, Eq. (1.3) reduces to

u

∗

(x, t; j, n) = u

n

j + (ux)

n

j



(x − xj ) − a (t − t

n

)



(j, n) ∈ Ω1 (1.6)

i.e., u

n

j

and (ux)

n

j

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

Let E2 be divided into non-overlapping rectangular regions (see Fig. 2(a)) referred to as basic conservation elements (BCEs). As depicted in Figs. 2(c)–2(e), (i) each (j, n) ∈ Ω1 is assigned with two BCEs,

i.e., CE−(j, n) and CE+(j, n); (ii) each BCE has one and only one pair of diagonally opposite vertices which

belong to Ω1; (iii) the space-time E2 can be filled by CE−(j, n) and CE+(j, n), (j, n) ∈ Ω1; and (iv) CE(j, n),

which is the union of CE−(j, n) and CE+(j, n), is referred to as a compounded conservation element (CCE).

NASA/TM—2013-218061 4
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, among the line segments forming the boundary

of CE+(j, n), AF and AD belong to SE(j, n), while EF and ED belong to SE(j + 1/2, n − 1/2). As will be

shown, by imposing the two local flux conservation conditions at each (j, n) ∈ Ω1, i.e.,

I

S(CE+(j,n))

~h

∗

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

and

I

S(CE−(j,n))

~h

∗

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

one can obtain two equations for two unknowns u

n

j

and (ux)

n

j

.

Let

ν

def = a

∆t

∆x

and (ux ̄)

n

j

def =

∆x

4

(ux)

n

j

(j, n) ∈ Ω1 (1.9)

where ν is the Courant number, and (ux ̄)

n

j

is the normalized version of (ux)

n

j

. Then, with the aid of Eqs. (1.4),

(1.6) and (1.9), it can be shown that Eqs. (1.7) and (1.8) are equivalent to

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

n

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

n−1/2

j+1/2

(j, n) ∈ Ω1 (1.10)

and

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

n

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

n−1/2

j−1/2

(j, n) ∈ Ω1 (1.11)

respectively. To simplify notation, in the above and hereafter we adopt a convention that can be explained

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

[u + (1 + ν)ux ̄]

n

j = u

n

j + (1 + ν)(ux ̄)

n

j

By adding Eqs. (1.10) and (1.11) together, one has

u

n

j =

1

2

n

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

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω1 (1.12)

Let 1 − ν

2 = 0, i.e., 1 6 − ν 6= 0 and 1 + ν = 0. Then Eqs. (1.10) and (1.11) can be divided by (1 6 − ν) and

(1 + ν), respectively. By subtracting the resulting equations from each other, one has

(ux ̄)

n

j =

1

2

n

[u − (1 + ν)ux ̄]

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω1 (1.13)

Hereafter, by definition the a scheme is formed by Eq. (1.12) and (1.13) for any ν. According the above

derivation, a solution to the a scheme will satisfy Eqs. (1.10) and (1.11) and therefore the two conservation

conditions Eqs. (1.7) and (1.8) for any ν. However, because Eqs. (1.7) and (1.8) ⇒ (i.e., “imply”) Eq. (1.12)

for any ν, but they ⇒ Eq. (1.13) only if ν

2 = 1, one concludes that a solution to Eqs. (1.7) and (1.8) may 6

not be a solution to the a scheme if ν

2 = 1. As such the a scheme represents a stronger condition than that

of Eq. (1.7) and (1.8).

The a scheme is an explicit and neutrally stable (i.e., non-dissipative) solver of Eq. (1.1) in its stability

domain |ν| < 1 [63]. Also, according to Eqs. (1.12) and (1.13), it has a space-time staggered stencil, i.e.,

any (j, n) ∈ Ω1 is associated with a stencil formed by the three mesh points (j, n), (j − 1/2, n − 1/2) and

(j + 1/2, n − 1/2).

Eqs. (1.7) and (1.8), i.e., the flux conservation conditions over the two BCEs CE−(j, n) and CE+(j, n)

for all (j, n) ∈ Ω1, are enforced by any solution to the a scheme. Because (i) the vector ~h

∗ at any surface

element lying on any interface separating two neighboring BCE is evaluated using the information from a

NASA/TM—2013-218061 5
single SE, and (ii) the unit outward normal vector on the surface element pointing outward from one of these

two neighboring BCEs is exactly the negative of that pointing outward from another BCE, one concludes that

the flux leaving one of these BCEs through the interface is the negative of that leaving another BCE through

the same interface. As a result of this flux cancellation, the local flux conservation relations Eqs. (1.7) and

(1.8) lead to a global flux conservation relation, i.e., the total flux of ~h

∗

leaving the boundary of any spacetime region that is the union of any combination of BCEs 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) ∈ Ω1 (1.14)

must follow from Eqs. (1.7) and(1.8). In fact, it can be shown that Eq. (1.14) is equivalent to Eq. (1.12).

There is a family of the dissipative extensions of the a scheme in which only the less stringent conservation

condition Eq. (1.14) is assumed [5,26,31,63]. Because Eq. (1.14) is equivalent to Eq. (1.12), for each of these

extensions, u

n

j

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

n

j

is evaluated in terms of the mesh variables at

(j − 1/2, n − 1/2) and (j + 1/2, n − 1/2) using an equation different from Eq. (1.13). In fact, (i) the a-

scheme is formed by Eq. (1.12) and

(ux ̄)

n

j =

1

2

n

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

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω1 (1.15)

where  is a real parameter, and (ii) the c-τ scheme is formed by Eq. (1.12) and

(ux ̄)

n

j =

1

2(1 + τ )

n

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

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω1 (1.16)

where τ is a real parameter 6= −1.

Note that:

(a) By comparing Eqs. (1.13) and (1.15), it is seen that the a scheme is the special case of the a- scheme

when  = 0.

(b) By using Jordan canonical form theorem [71], it can be shown that the a- scheme is von Neumann

stable ⇔ (i.e., “if and only if”) either

0 ≤  ≤ 1 and |ν| < 1 (1.17)

or

 = |ν| = 1 (1.18)

Also one can show that the a- scheme becomes more dissipative as the value of  increases from 0 to 1

[5].

(c) By using Jordan canonical form theorem [71], it is shown rigorously in [46] that the c-τ scheme is von

Neumann stable ⇔

ν

2 ≤ 1, τ ≥ τo(ν2

), and (ν

2

, τ ) = (1 6 , 1) (1.19)

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

(1.20)

It can be shown that: (i) τo(s) is continuous at s = 0, i.e.,

lim

s→0+

τo(s) = τo(0) = 0 (1.21)

NASA/TM—2013-218061 6
(ii) τo(s) is consistently defined at s = 3/11, i.e., both the second and third expressions on the right

side of Eq. (1.20) implies that τo(3/11) = 1/3; (iii) τo(s) is continuously differentiable at s = 3/11, i.e.,

lim

s→ 3

11

−

τ

0

o

(s) = lim

s→ 3

11

+

τ

0

o

(s) = 121/90 (1.22)

where τ

0

o

(s)

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

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

s < τo(s) <

√

s, 0 < s < 1 (1.23)

Moreover, (i) the stability conditions Eq. (1.19) agree completely with those generated through numerical

experiments [31] (Note: the reader is warned that the condition τ ≥ τo(ν

2

) that appears in Eq. (1.19)

here is expressed as τ ≥ τo(|ν|) in Eq. (3.12) of [31], i.e., the function τo defined here is different from

the function τo introduced in [31]), and (ii) for any given fixed value of |ν| < 1, the c-τ scheme tends to

become more dissipative as the value of τ increases from its minimum stability bound τo(ν

2

).

The a scheme and its extensions were defined over Ω1, the set of space-time staggered mesh points marked

by dots in Fig. 2(a). Obviously, the a scheme can also be defined over Ω2, the set of all space-time staggered

mesh points (j, n) in E2 where n = 0, ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, j = n, n ± 1, n ± 2, n ± 3, . . ..

By definition, each (j, n) ∈ Ω2 is associated with an whole-integer value of (j + n), and represented by a

mesh point not marked by dots in Fig. 2(a). Because one can carry out the time marching (using the a

scheme or any one of its 1D extensions) exclusively over either Ω1 or Ω2, there is no need to carry out a

marching over both Ω1 and Ω2. Note that this conclusion generally remains valid even if an 1D nonuniform

spatial mesh is used [31]. The exception may occur for a marching scheme where the time-step size varies

from one spatial region to another such as the local time stepping procedure described in [42].

Note that the two disjoint sets Ω1 and Ω2 have the property that, if a space-time mesh point (j, n)

belongs to Ω1 (Ω2), then its four immediate space-time neighbors (j ±1/2, n) and (j, n±1/2) must belong to

Ω2 (Ω1). As a result, for a scheme with a space-time staggered stencil such as the a scheme or the classical

leapfrog scheme, this implies that a time marching carried out over Ω1 is completely decoupled from that

carried out over Ω2. As such one only needs to carry out time marching over one of these two sets. For

this reason, the problem of “solution decoupling” (which is famously associated with the leapfrog scheme)

usually can be avoided in the 1D case.

However, as will be shown shortly, for a space-time mesh built from a multidimensional unstructured

spatial mesh, one generally cannot define two disjoint sets of space-time mesh points which possess the same

property of Ω1 and Ω2 discussed above. As such, the problem of solution decoupling may become a real

nuisance in multidimensional CESE simulations involving unstructured spatial meshes.

Is was shown earlier that, for an 1D CESE scheme, (i) each mesh point is assigned with two BCEs; and

(ii) for each dependent physical variable u, each mesh point (j, n) is assigned with two independent mesh

variables u

n

j

and (ux)

n

j

. Similarly [12,13], for a 2D (3D) CESE scheme, (i) each space-time solution point

(to be defined shortly) is assigned with three (four) BCEs and one CCE, with the CCE being the union

of the three (four) BCEs; and (ii) for each dependent physical variable, each space-time solution point is

assigned with three (four) independent mesh variables—note that there are two (three) independent spatial

derivatives in a problem with two (three) spatial dimension. Because a triangle and a tetrahedron have three

and four sides, respectively, triangles and tetrahedra, respectively, become the simplest and most natural

building blocks of the spatial meshes for 2D and 3D CESE schemes [10–13,16]. As such, in the following

discussion about solution decoupling in a CESE simulation involving an unstructured 2D spatial domain, we

consider an unstructured domain built from triangles.

Consider a spatial domain formed by triangles of arbitrary shapes (see Fig. 3). Here (i) any two

neighboring triangles share a common side, and (ii) the vertices and centroids of triangles are marked by

dots and circles, respectively. As shown in Fig. 3, the triangle 4F BD has three neighboring triangles. Let

(i) point G be the centroid of 4F BD, and (ii) points A, C and E, respectively, be the centroids of the three

neighbors of 4F BD. Then points A, B, C, D, E, and F form a hexagon. The centroid of the hexagon

NASA/TM—2013-218061 7
is denoted by G0 and marked by a cross. Hereafter point G0

is referred to as the spatial solution point of

4F BD. In a similar manner, we define spatial solution points A0, C

0

, E0,... of other triangles.

A space-time solution point is a mesh point which resides at a time level n = 0, 1/2, 1, 3/2, 2, ... (t =

t

n = n∆t at the nth time level) and has a spatial solution point being its spatial projection. Given any

spatial solution point G0 and any time level n, let (G0, n) denote the space-time solution point which resides

at the nth time level and has point G0 being its spatial projection. As explained in [13], the three BCEs

associated with (G0, n) are constructed such that their spatial projections are the quadrilaterals GF AB,

GBCD, and GDEF, respectively. Also, for a reason explained in [13], independent mesh variables are

stored at space-time solution points.

As shown in [13], in a 2D CESE scheme using a triangular spatial mesh, the mesh variables at (G0, n) are

evaluated in terms of those at (A0, n−1/2), (C

0

, n−1/2) and (E0, n−1/2), i.e., the scheme has a space-time

staggered stencil formed by one point at the nth time level and three points at the (n − 1/2)th time level.

Also it was shown in [13] that, given an unstructured triangular spatial mesh, the time marching cannot be

carried out independently over two disjoint sets of space-time solution points unless the triangles forming

the spatial mesh can be divided into two disjoint sets Ψ1 and Ψ2 such that two triangles sharing a common

side always belong to two different sets. As will be shown, because point F is the common vertex of an odd

number (5) of triangles, the triangles depicted in Fig. 3 cannot be divided into Ψ1 and Ψ2.

To prove the last proposition by contradiction, let 4F BD ∈ Ψ1. Then because 4F HB and 4F BD

share a common side, by definition, 4F HB ∈ Ψ2. Similarly, because (i) 4F NH and 4F HB share a

common side, (ii) 4FMN and 4F NH share a common side, and (iii) 4F DM and 4FMN share a

common side, one concludes that 4F NH ∈ Ψ1, 4FMN ∈ Ψ2, and 4F DM ∈ Ψ1. As such, we have arrive

at the contradiction that both 4F BD and 4F DM belong to Ψ1 in spite of the fact that they share a

common side. Because we will reach a similar contradiction by assuming 4F BD ∈ Ψ2, the proposition has

been proved.

The triangles depicted in Fig. 3 cannot be divided into two disjoint sets Ψ1 and Ψ2. As such, according

to another proposition stated earlier, the CESE time marching cannot be carried out independently over two

disjoint sets of space-time solution points constructed from the triangular spatial mesh depicted in Fig. 3.

Moreover, because the incidence that an odd number of triangles sharing a common vertex occurs regularly

in an unstructured triangular mesh, the CESE time marching generally cannot be carried out independently

over two disjoint sets of space-time solution points if they are constructed from an unstructured triangular

spatial mesh.

Similarly, for a 3D CESE scheme using an unstructured spatial mesh built from tetrahedra, it can

be shown that generally the time marching cannot be carried out independently over two disjoint sets of

space-time staggered solution points.

According to the above discussions, a CESE multidimensional time marching using an unstructured

triangular or tetrahedral spatial mesh generally will involve all space-time solution points. This coupled

with the fact that a CESE scheme has a space-time staggered stencil implies that solution decoupling may

become a real nuisance. To deal with this problem once and for all, recently new CESE schemes have been

developed such that possibility of solution decoupling can be removed completely. The description of these

new schemes is one of the topics covered in the current paper. As will be shown, (i) these new schemes are

constructed without compromising the basic CESE ideas, and (ii) their effectiveness has been verified by

real-world applications. [62, 64–67].

To describe another topic of the current paper, note that the spatial solution point G0 depicted in Fig. 3

is the centroid of the hexagon ABCDEF where points A, C and E, respectively, are rather arbitrarily

designated as the centroids of the three neighbors of 4F BD. In fact, to insure that CEs will not overlap in

space and thus no loss of local and global flux conservation will occur, it only requires that points such as

G, A, C and E lie in the interiors of 4F BD, 4F HB, 4BLD, and 4F DM, respectively.

For the rather benign triangular domain depicted in Fig. 3, the hexagon ABCDEF happens to be

convex, i.e., a line segment joining any two points in this hexagon always lies entirely in it. Because the

centroid of a convex polygon always lies within its interior, the centroid (i.e., point G0) of the hexagon

ABCDEF lies within its interior. However, this may not be true for a more pathological triangular domain

NASA/TM—2013-218061 8
formed from triangles with very large aspect ratios and sharply different orientations. As an example,

consider Fig. 4 where 4F BD is again surrounded by three neighboring triangles 4F HB, 4BLD and

4F DM. However, in Fig. 4, the aspect ratio of each of the four triangles depicted is much larger than

that of its counterpart depicted in Fig. 3. Moreover, in Fig. 4, the orientations of the sides BD and F D of

4F BD are sharply different from those of the sides BH and F H of 4F HB, respectively. Let points G, A,

C and E depicted in Fig. 4 again be the centroids of 4F BD, 4F HB 4BLD, and 4F DM, respectively.

Then the resulting hexagon ABCDEF may become highly concave and thus its centroid G0(not shown in

Fig. 4) may no longer lie in its interior.

Similarly, a CESE tetrahedral mesh is built from tetrahedra where any two neighboring tetrahedra

share a common triangular face. Given any two neighboring tetrahedra, in the original CESE tetrahedralmesh construction, the 5-vertex polyhedron with its vertices being (i) the centroids of these two neighbors

and (ii) the three vertices of the triangular face shared by them, becomes the spatial projection of a BCE.

Because each tetrahedron is associated with four neighbors, it is associated with four 5-vertex polyhedra and

their union—a 8-vertex polyhedron. The 8-vertex polyhedron and its centroid, respectively, are the spatial

projections of a CCE and a spatial solution point. For some pathological case where the 8-vertex polyhedron

becomes highly concave, its centroid (i.e., the spatial solution point associated with this polyhedron) may

lie outside of its interior.

In a real-world CFD simulation near a solid wall, one often encounters a case where a boundary with

high curvature or sharp corner is surrounded by triangular/tetrahedral meshes of extremely high aspect

ratio (up to 106). For such an extreme case, a spatial solution point may lie far outside of the associated

hexagon/8-vertex polyhedron. It could even be embedded beyond a solid-wall boundary and causes serious

numerical difficulties. To overcome this problem, a new way of building spatial meshes from triangles and

tetrahedra has been developed. Specifically, the centroids of triangles depicted in Fig. 3 have been replaced

by the corresponding in-centers or other specially-defined interior points in the new construction. However,

the solution point is still the centroid of the new hexagon formed by vertices and in-centers (or other

specially-defined interior points). Similarly, centroids of tetrahedra also are replaced by corresponding incenters or other specially-defined interior points in the new construction. According to the numerical results

obtained [62,64,65], the new construction has the following advantages: (i) it can properly model curved

boundaries—without it, computed surface gradient such as heat flux will be much less accurate (and noisier

in many cases); (ii) it can reduce aspect ratios of triangular or tetrahedral cells used in gradient calculations

in a highly stretched mesh and, as a result, can improve the overall accuracy for triangular and tetrahedral

meshes; (iii) for vortex-shedding problems, better resolution is observed by using the in-center formulation.

In this paper we will describe this new development and also provide the theoretical justification for the new

construction.

As mentioned earlier, mesh values of the physical dependent variables and their spatial derivatives are

treated as independent mesh variables in a CESE scheme. As it turns out, for a case where a triangular or

tetrahedral mesh is used, numerical evaluation of the spatial derivatives using the established CESE approach

tends to become less accurate as the mesh cell aspect ratio becomes larger. In the latest development, the

root cause behind this tendency was clearly identified and several remedies to overcome it were found through

a rigorous mathematical analysis. However, because of the length of the current paper and the complexity

of the mathematics involved, this latest work will be presented in another paper.

Because it is rather simple and straightforward, in the following, first we will describe the new way of

building spatial meshes from triangles and tetrahedra. The new flux-conserving solution-coupling procedures

will be presented in Sec. 3.

2. New construction of triangular and tetrahedral meshes

2.1. Triangular meshes

Again consider the triangle 4F BD and the hexagon ABCDEF depicted in Fig. 4. As noted earlier,

(i) because it is the centroid of the hexagon ABCDEF, the spatial solution point G0(not shown in Fig. 4)

always lies in the interior of the hexagon if it is convex but may lie outside of it if it is concave; and (ii)

serious numerical difficulties may arise if point G0lies outside of the hexagon ABCDEF. As such, to avoid

NASA/TM—2013-218061 9
numerical difficulties, it is essential to ensure that ABCDEF is a convex hexagon.

Because any triangle such as 4F BD is always convex, the hexagon ABCDEF will also be convex if its

boundary more or less coincides with that of 4F BD. In other words, the hexagon is convex if points A, C,

and E are located close enough to the sides BF, BD and F D, respectively. Unfortunately, because (i) A is

the centroid of 4F HB, and (ii)

|BF|  |BH| and |BF|  |F H| (2.1)

(hereafter, as an example, the length of BF is denoted by |BF|), point A generally is far away from the

side BF. As such the hexagon ABCDEF may become concave. For this reason and the previously stated

requirement that point A, C, and E must lie in the interiors of 4F HB, 4BLD, and 4F DM, respectively,

in the following, we will describe a new approach in which the locations of A, C, and E are chosen so that:

(i) each of them again is an interior point of the associated triangle, and (ii) a point such as A will remain

close to the side BF even for the pathological case Eq. (2.1).

Consider Fig. 5. Here (i) points P1, P2, P3 are the vertices of the triangle 4P1P2P3 which lies on the

x-y plane and has its area

A(4P1P2P3) > 0 (2.2)

(ii) P0 is a point on the x-y plane to be defined shortly; (iii) for each k = 0, 1, 2, 3, xk and yk are the x- and

y- coordinates of point Pk, respectively; (iv) point Q1 is the point on P2P3 (or its extension) which meets

the condition

P0Q1 ⊥ P2P3 (2.3)

i.e., Q1 is the perpendicular projection of P0 on the straight line which contains P2P3; (v) point Q2 is the

point on P3P1 (or its extension) which meets the condition

P0Q2 ⊥ P3P1 (2.4)

i.e., Q2 is the perpendicular projection of P0 on the straight line which contains P3P1; (vi) point Q3 is the

point on P1P2 (or its extension) which meets the condition

P0Q3 ⊥ P1P2 (2.5)

i.e., Q3 is the perpendicular projection of P0 on the straight line which contains P1P2; (vii) α1, α2, and α3

are the internal angles of 4P1P2P3 facing the sides P2P3, P3P1, and P1P2, respectively; (viii)

`1

def = |P2P3|, `2

def = |P3P1| and `3def = |P1P2| (2.6)

(ix)

hk

def = |P0Qk| k = 1, 2, 3 (2.7)

and (x) the x-y-z coordinate system depicted in Fig. 5 is a right-handed system, i.e.,

~ex × ~ey = ~ez, ~ey × ~ez = ~ex and ~ez × ~ex = ~ey (2.8)

where ~ex, ~ey, and ~ez are the unit vectors in the x-, y-, and z- directions, respectively.

Note that, according to Eq. (2.6) and Fig. 5, for each k = 1, 2, 3, `k denotes the length of the side of

4P1P2P3 facing the vertex Pk. Also from the definition of points Q1, Q2, and Q3, and that of hk, k = 1, 2, 3,

one concludes that h1, h2, and h3, respectively, are the (shortest) distances from P0 to the straight lines that

contain the sides P2P3, P3P1, and P1P2, respectively.

For each k = 0, 1, 2, 3, let −→P k be the position vector of point Pk, i.e.,

−→P k

def = xk ~ex + yk ~ey k = 0, 1, 2, 3 (2.9)

NASA/TM—2013-218061 10
Also let −→P 0 be an weighted average of −→P k, k = 1, 2, 3, i.e.,

−→P 0 = w1

−→P 1 + w2

−→P 2 + w3

−→P 3 (2.10)

where w1, w2, and w3 are real weight factors with

w1 + w2 + w3 = 1 (2.11)

Moreover, given a pair of integers k1 and k2 with k1 6= k2 and k1, k2 = 0, 1, 2, 3, let −−−−→ Pk1 Pk2 be the

displacement vector joining points Pk1

and Pk2

and pointing in the direction from Pk1

to Pk2

, i.e.,

−−−−→ Pk1 Pk2

def =

−→P k2 −

−→P k1

, k1 6= k2 and k1, k2 = 0, 1, 2, 3 (2.12)

Then, by using Eq. (2.12) and other well-known properties of the vector product such as

−→A ×−→A = 0 and −→A ×−→B = −−→B ×−→A (2.13)

for any vectors −→A and −→B, one can show that

−−−→ P1P2 ×−−−→ P1P3 =

−→P 2 −

−→P 1



×

−→P 3 −

−→P 1



=

−→P 1 ×−→P 2 +−→P 2 ×−→P 3 +−→P 3 ×−→P 1 (2.14)

Similarly, with the aid of Eqs. (2.10)–(2.14), one can show that

−−−→ P0P2 ×

−−−→ P0P3 = w1

−→P 1 ×

−→P 2 +−→P 2 ×−→P 3 +−→P 3 ×−→P 1



= w1

−−−→ P1P2 ×

−−−→ P1P3



(2.15)

−−−→ P0P3 ×−−−→ P0P1 = w2

−→P 1 ×

−→P 2 +−→P 2 ×−→P 3 +−→P 3 ×−→P 1



= w2

−−−→ P1P2 ×

−−−→ P1P3



(2.16)

−−−→ P0P1 ×−−−→ P0P2 = w3

−→P 1 ×

−→P 2 +−→P 2 ×−→P 3 +−→P 3 ×−→P 1



= w3

−−−→ P1P2 ×

−−−→ P1P3



(2.17)

Next, by using (i) a well-known property of the vector product, (ii) the fact that A(4P1P2P3) denotes

the area of 4P1P2P3, and (iii) Fig. 5, one has

A(4P1P2P3) = (|

−−−→ P1P2||−−−→ P1P3|sin α1)/2 = |−−−→ P1P2 ×−−−→ P1P3|/2 (2.18)

Similarly, with the aid of Eqs. (2.15)–(2.18), one can show that

A(4P0P2P3) = |

−−−→ P0P2 ×−−−→ P0P3|/2 = |w1| × A(4P1P2P3) (2.19)

A(4P0P3P1) = |

−−−→ P0P3 ×−−−→ P0P1|/2 = |w2| × A(4P1P2P3) (2.20)

and

A(4P0P1P2) = |

−−−→ P0P1 ×−−−→ P0P2|/2 = |w3| × A(4P1P2P3) (2.21)

where A(4P0P2P3), A(4P0P3P1), and A(4P0P1P2) denote the areas of 4P0P2P3, 4P0P3P1, and 4P0P1P2,

respectively. In turn, Eqs. (2.19)–(2.21) ⇒

A(4P0P2P3) + A(4P0P3P1) + A(4P0P1P2) = (|w1| + |w2| + |w3|) × A(4P1P2P3) (2.22)

Because a triangle such as 4P1P2P3 is intrinsically convex, according to a theorem on convex sets, point

P0 defined by Eqs. (2.10) and (2.11) lies on the boundary or in the interior of 4P1P2P3 ⇔

0 ≤ wk ≤ 1 k = 1, 2, 3 (2.23)

NASA/TM—2013-218061 11
Moreover, because Eq. (2.2) implies that the interior of 4P1P2P3 is not empty, P0 lies in the interior of

4P1P2P3 ⇔

0 < wk < 1 k = 1, 2, 3 (2.24)

For the case where P0 lies on the boundary or in the interior of 4P1P2P3, Eqs. (2.11), (2.22), and (2.23)

now ⇒ the following obvious result:

A(4P0P2P3) + A(4P0P3P1) + A(4P0P1P2) = A(4P1P2P3) (2.25)

On the other hand, for the case where P0 lies outside of 4P1P2P3, Eq. (2.23) is no longer true, i.e., some

of the weight factors must be negative. Because Eq. (2.11) requires that some of the weight factors must be

positive, we arrive at the conclusion that

|w1| + |w2| + |w3| > |w1 + w2 + w3| = 1 (2.26)

if P0 lies outside of 4P1P2P3. In turn, Eqs. (2.22) and (2.26) ⇒

A(4P0P2P3) + A(4P0P3P1) + A(4P0P1P2) > A(4P1P2P3) (2.27)

if P0 lies outside of 4P1P2P3.

Note that, for any real parameter ζ, Eqs. (2.2) and (2.6) ⇒

`k > 0 and (`k)

ζ > 0 k = 1, 2, 3; −∞ < ζ < +∞ (2.28)

As such the weight factors

wk(ζ)

def =

(`k)

ζ

(`1)

ζ + (`2)ζ + (`3)ζ

k = 1, 2, 3; −∞ < ζ < +∞ (2.29)

are well defined and satisfy the conditions

0 < wk(ζ) < 1 k = 1, 2, 3, −∞ < ζ < +∞ (2.30)

and

w1(ζ) + w2(ζ) + w3(ζ) = 1 − ∞ < ζ < +∞ (2.31)

In this subsection hereafter we consider only the special case where

wk = wk(ζ) k = 1, 2, 3; −∞ < ζ < +∞ (2.32)

For this special case, Eqs. (2.30)–(2.32) imply that both Eqs. (2.11) and (2.24) are satisfied for any real

number ζ. Thus, for the case Eq. (2.32), (i) point P0 defined by Eq. (2.10) must always lie in the interior of

4P1P2P3, and (ii) Eqs. (2.19)–(2.22) reduce to

A(4P0P2P3) = w1 × A(4P1P2P3) (2.33)

A(4P0P3P1) = w2 × A(4P1P2P3) (2.34)

A(4P0P1P2) = w3 × A(4P1P2P3) (2.35)

and Eq. (2.25), respectively. Moreover, for the special case ζ = 0, we have

w1 = w2 = w3 = 1/3 and −→P 0 =

−→P 1 +

−→P 2 +−→P 3



/3 (ζ = 0) (2.36)

NASA/TM—2013-218061 12
Because the position vector of the centroid of a triangle is the simple average of those of its three vertices,

Eq. (2.36) ⇒ point P0 is the centroid of 4P1P2P3 if ζ = 0.

Note that, for a polygon other than a triangle (which is the simplest polygon), generally it is not true

that the position vector of its centroid is the simple average of those of its vertices. As an example, the

position vector of the centroid of a quadrilateral may not be the simple average of those of its four vertices

unless the quadrilateral is a parallelogram.

Next, with the aid of Eqs. (2.6) and (2.7), and Fig. 5, one can see that

A(4P0P2P3) = h1`1

2

, A(4P0P3P1) = h2`2

2

and A(4P0P1P2) = h3`3

2

(2.37)

In turn, by substituting Eq. (2.37) into Eqs. (2.33)–(2.35) and using Eqs. (2.29) and (2.32), one has

hk`k

2

= wkA(4P1P2P3) = (`k)

ζ × A(4P1P2P3)

(`1)

ζ + (`2)ζ + (`3)ζ

k = 1, 2, 3; −∞ < ζ < +∞ (2.38)

and thus

h1

(`1)

ζ−1

=

h2

(`2)

ζ−1

=

h3

(`3)

ζ−1

=

2 × A(4P1P2P3)

(`1)

ζ + (`2)ζ + (`3)ζ

− ∞ < ζ < +∞ (2.39)

Because (i) h1, h2 and h3 are the distances from point P0 to the sides P2P3, P3P1, and P1P2, respectively,

and (ii) `1, `2, and `3 are the lengths of P2P3, P3P1, and P1P2, respectively, Eq. (2.39) states that, for a

given 4P1P2P3 and a given ζ, the distance from point P0 to any side of 4P1P2P3 is proportional to the

(ζ − 1)th power of the length of this side.

In fact, for the case ζ = 0 (i.e., when P0 is the centroid of 4P1P2P3), Eq. (2.39) reduces to

hk =



2 × A(4P1P2P3)

3



1

`k

k = 1, 2, 3 if ζ = 0 (2.40)

i.e., the distance between the centroid P0 and a side of 4P1P2P3 is inversely proportional to the side length.

Thus the centroid is farthest away from the shortest side and closest to the longest side—a result consistent

with what we observe from the triangle 4F HB depicted in Fig. 4.

On the other hand, for the case ζ = 1, Eq. (2.39) reduces to

h1 = h2 = h3 =

2 × A(4P1P2P3)

`1 + `2 + `3

if ζ = 1 (2.41)

i.e., point P0 has the same distance to the three sides of 4P1P2P3 if ζ = 1. By definition, this implies that

P0 is the in-center of 4P1P2P3 for the case ζ = 1.

In general one can infer from Eq. (2.39) that:

(a) For the case ζ > 1, point P0 is farthest away from the longest side of 4P1P2P3 and closest to its shortest

side. The trend becomes stronger as the value of ζ increases from 1.

(b) For the case ζ < 1, point P0 is farthest away from the shortest side of 4P1P2P3 and closest to its longest

side. The trend becomes stronger as the value of ζ decreases from 1.

Let the interior point P0 of any triangle 4P1P2P3 defined by Eqs. (2.10) and (2.32) be referred to as

the interior point of this triangle associated with the parameter ζ. Then, from the above observations and

the discussions given in Sec. 1, it becomes clear that a spatial solution point such as the point G0 depicted

in Fig. 3 will be much less likely to lie outside of the associated hexagon ABCDEF and causes numerical

difficulties if each of the centroidal mesh points A, C, G and E depicted in Figs. 3 and 4 is replaced by a an

interior point associated with a parameter ζ ≥ 1. In fact, it has been shown through numerical experiments

that the numerical difficulties associated with spatial solution points being embedded beyond a solid-wall

boundary can be overcome by replacing centroids of triangles with in-centers (i.e., interior points associated

with ζ = 1) in the construction of a spatial triangular mesh.

NASA/TM—2013-218061 13
2.2. Tetrahedral meshes

Consider a 3D spatial domain formed by tetrahedra with the understanding that any two neighboring

tetrahedra share a common triangular face. As an example, consider the tetrahedron P1P2P3P4 depicted in

Fig. 6(a) where P1, P2, P3 and P4 are the vertices of the tetrahedron. The tetrahedron has four triangular

faces, i.e., 4P1P2P3, 4P2P3P4, 4P3P4P1, and 4P4P1P2. As shown in Fig. 6(b), the tetrahedron P1P2P3P4

shares these faces with four neighboring tetrahedra, i.e., P1P2P3D, P2P3P4A, P3P4P1B, and P4P1P2C,

respectively. Note that, to simplify notation, the same symbols used in Sec. 2.1 may also be used in this

subsection albeit that those used here generally denote objects in the 3D x-y-z space.

Because a tetrahedron is intrinsically convex in a 3D x-y-z space, its centroid must lie in its interior if

its volume > 0. Hereafter, only tetrahedra of volume > 0 will be considered. In particular, we assume that

V (P1P2P3P4)

def = the volume of the tetrahedron P1P2P3P4 > 0 (2.42)

Let the points P0, A0, B0, C0, and D0 depicted in Fig. 6(b) be the centroids of the tetrahedra, P1P2P3P4,

P2P3P4A, P3P4P1B, P4P1P2C, and P1P2P3D, respectively. Then each of these five centroids lies in the

interior of the associated tetrahedron. As such the four polyhedra P0P2P3P4A0, P0P3P4P1B0, P0P4P1P2C0,

and P0P1P2P3D0 (each of them has five vertices, six triangular faces, and nine edges) do not overlap. In

the original CESE 3D tetrahedral mesh construction, (i) each of these 5-vertex polyhedra is the spatial

projection of a BCE, (ii) the union of the above four 5-vertex polyhedra is a polyhedron with 8 vertices,

12 triangular faces, and 18 edges, and also is the spatial projection of a CCE, and (iii) the centroid of the

8-vertex polyhedron referred to above is a spatial solution point, i.e., the spatial projection of a space-time

solution point.

As in the 2D case, for some pathological case where an 8-vertex polyhedron become highly concave, its

centroid (i.e., the associated spatial solution point) may lie outside of the polyhedron and causes numerical

difficulties. In the following, we explain how these difficulties can be overcome by a new construction of 3D

tetrahedral meshes.

Consider Fig. 6(a) again. Here (i) P0 is a point in the x-y-z space to be defined shortly and not necessary

the centroid of the tetrahedron P1P2P3P4; (ii) for each k = 0, 1, 2, 3, 4, xk, yk, and zk are the x-, y- and zcoordinates of point Pk, respectively; (iii) point Q1 (not shown in Fig. 6(a)) is the point that lies on the

plane containing the face 4P2P3P4, and also meets the condition

P0Q1 ⊥ 4P2P3P4 (2.43)

i.e., Q1 is the perpendicular projection of P0 on the plane containing 4P2P3P4; (iv) point Q2 (not shown)

is the point that lies on the plane containing the face 4P3P4P1, and also meets the condition

P0Q2 ⊥ 4P3P4P1 (2.44)

i.e., Q2 is the perpendicular projection of P0 on the plane containing 4P3P4P1; (v) point Q3 (not shown) is

the point that lies on the plane containing the face 4P4P1P2, and also meets the condition

P0Q3 ⊥ 4P4P1P2 (2.45)

i.e., Q3 is the perpendicular projection of P0 on the plane containing 4P4P1P2; (vi) point Q4 (not shown)

is the point that lies on the plane containing the face 4P1P2P3, and also meets the condition

P0Q4 ⊥ 4P1P2P3 (2.46)

i.e., Q4 is the perpendicular projection of P0 on the plane containing 4P1P2P3; (vii)

A1

def = A(4P2P3P4), A2

def = A(4P3P4P1), A3

def = A(4P4P1P2), and A4

def = A(4P1P2P3) (2.47)

NASA/TM—2013-218061 14
where A(4P2P3P4), A(4P3P4P1), A(4P4P1P2), and A(4P1P2P3) are the areas of the four triangular faces

4P2P3P4, 4P3P4P1, 4P4P1P2, and 4P1P2P3, respectively; (viii)

hk

def = |P0Qk| k = 1, 2, 3, 4 (2.48)

and (ix) the x-y-z coordinate system depicted in Fig. 6(a) is also a right-handed system which satisfies

Eq. (2.8).

Note that, according to Eq. (2.47) and Fig. 6(a), for each k = 1, 2, 3, 4, Ak denotes the area of the

face of the tetrahedron P1P2P3P4 which is opposite to the vertex Pk. Also from the definition of points

Q1, Q2, Q3, and Q4 and that of hk, k = 1, 2, 3, 4, one concludes that h1, h2, h3, and h4, respectively, are

the (shortest) distances from P0 to the planes that contain the faces 4P2P3P4, 4P3P4P1, 4P4P1P2, and

4P1P2P3, respectively.

For each k = 0, 1, 2, 3, 4, let −→P k be the position vector of point Pk, i.e.,

−→P k

def = xk ~ex + yk ~ey + zk ~ez k = 0, 1, 2, 3, 4 (2.49)

Also let −→P 0 be an weighted average of −→P k, k = 1, 2, 3, 4, i.e.,

−→P 0 = w1−→P 1 + w2

−→P 2 + w3

−→P 3 + w4

−→P 4 (2.50)

where w1, w2, w3, and w4 are real weight factors with

w1 + w2 + w3 + w4 = 1 (2.51)

Moreover, given a pair of integers k1 and k2 with k1 6= k2 and k1, k2 = 0, 1, 2, 3, 4, let −−−−→ Pk1 Pk2 be the

displacement vector joining points Pk1

and Pk2

and pointing in the direction from Pk1

to Pk2

, i.e.,

−−−−→ Pk1 Pk2

def =

−→P k2 −

−→P k1

, k1 6= k2 and k1, k2 = 0, 1, 2, 3, 4 (2.52)

Then, by using Eq. (2.52) and other well-known properties of the vector product such as Eq. (2.13), one can

show that Eq. (2.14) is still valid in the current 3D case. In turn, with the aid of Eqs. (2.14) and (2.52), and

other well-known properties of the scalar triple product such as

(

−→A ×−→B) ·−→C = (−→B ×−→C ) ·−→A = (−→C ×−→A) ·−→B = −(−→A ×−→C ) ·−→B (2.53)

and

(

−→A ×−→A) ·−→B = (−→A ×−→B) ·−→A = (−→B ×−→A) ·−→A = 0 (2.54)

for any vectors −→A,

−→B, and −→C , one can show that

∆

def = (−−−→ P1P2 ×

−−−→ P1P3)·−−−→ P1P4 = (−→P 1 ×−→P 2)·−→P 4 + (−→P 2 ×−→P 3)·−→P 4 + (−→P 3 ×−→P 1)·−→P 4 + (−→P 3 ×−→P 2)·−→P 1 (2.55)

Moreover, by using Eqs. (2.13), (2.50) and (2.52), one has

−−−→ P0P1 ×−−−→ P0P2 = (1 − w1 − w2)−→P 1 ×−→P 2 +



w3

−→P 3 + w4

−→P 4



×

−→P 1 −

−→P 2



(2.56)

Next, by using Eqs. (2.50)–(2.56), one can show that

(

−−−→ P0P1 ×−−−→ P0P2) ·−−−→ P0P3

= w4

h

(

−→P 1 ×−→P 2) ·−→P 3 + (−→P 2 ×−→P 1) ·−→P 4 + (−→P 3 ×−→P 4) ·−→P 1 + (−→P 4 ×−→P 3) ·−→P 2

i

= −w4 ∆

(2.57)

NASA/TM—2013-218061 15
In turn, by carrying out the index mapping

(1, 2, 3, 4) → (2, 3, 4, 1)

of Eq. (2.57) repeatedly and using Eqs. (2.53)–(2.55), one obtains

(

−−−→ P0P2 ×−−−→ P0P3) ·−−−→ P0P4

= w1

h

(

−→P 2 ×−→P 3) ·−→P 4 + (−→P 3 ×−→P 2) ·−→P 1 + (−→P 4 ×−→P 1) ·−→P 2 + (−→P 1 ×−→P 4) ·−→P 3

i

= w1 ∆

(2.58)

(

−−−→ P0P3 ×−−−→ P0P4) ·−−−→ P0P1

= w2

h

(

−→P 3 ×−→P 4) ·−→P 1 + (−→P 4 ×−→P 3) ·−→P 2 + (−→P 1 ×−→P 2) ·−→P 3 + (−→P 2 ×−→P 1) ·−→P 4

i

= −w2 ∆

(2.59)

and

(

−−−→ P0P4 ×−−−→ P0P1) ·−−−→ P0P2

= w3

h

(

−→P 4 ×−→P 1) ·−→P 2 + (−→P 1 ×−→P 4) ·−→P 3 + (−→P 2 ×−→P 3) ·−→P 4 + (−→P 3 ×−→P 2) ·−→P 1

i

= w3 ∆

(2.60)

For any four different integers k1, k2, k3, k4 = 0, 1, 2, 3, 4, let

V (Pk1 Pk2 Pk3 Pk4

)

def = the volume of the tetrahedron with the vertices Pk1

, Pk2

, Pk3

, and Pk4

(2.61)

Then by using the definition of the scalar triple product along with Eqs. (2.55) and (2.57)–(2.60), one has

V (P1P2P3P4) = |(

−−−→ P1P2 ×−−−→ P1P3) ·−−−→ P1P4|/6 = |∆|/6 (2.62)

V (P0P2P3P4) = |(

−−−→ P0P2 ×−−−→ P0P3) ·−−−→ P0P4|/6 = |w1 ∆|/6 = |w1| × V (P1P2P3P4) (2.63)

V (P0P3P4P1) = |(

−−−→ P0P3 ×−−−→ P0P4) ·−−−→ P0P1|/6 = |w2 ∆|/6 = |w2| × V (P1P2P3P4) (2.64)

V (P0P4P1P2) = |(

−−−→ P0P4 ×−−−→ P0P1) ·−−−→ P0P2|/6 = |w3 ∆|/6 = |w3| × V (P1P2P3P4) (2.65)

and

V (P0P1P2P3) = |(

−−−→ P0P1 ×−−−→ P0P2) ·−−−→ P0P3|/6 = |w4 ∆|/6 = |w4| × V (P1P2P3P4) (2.66)

Eqs. (2.63)–(2.66) ⇒

V (P0P2P3P4) + V (P0P3P4P1) + V (P0P4P1P2) + V (P0P1P2P3)

= (|w1| + |w2| + |w3| + |w4|) × V (P1P2P3P4)

(2.67)

Because the tetrahedron P1P2P3P4 is intrinsically convex in the x-y-z space, point P0 defined by

Eqs. (2.50) and (2.51) lies on the boundary or in the interior of the tetrahedron P1P2P3P4 ⇔

0 ≤ wk ≤ 1 k = 1, 2, 3, 4 (2.68)

Moreover, because Eq. (2.42) implies that the interior of the tetrahedron P1P2P3P4 is not empty, P0 lies in

the interior of the tetrahedron P1P2P3P4 ⇔

0 < wk < 1 k = 1, 2, 3, 4 (2.69)

For the case where P0 lies on the boundary or in the interior of the tetrahedron P1P2P3P4, Eqs. (2.51),

(2.67), and (2.68) now ⇒ the following obvious result:

V (P0P2P3P4) + V (P0P3P4P1) + V (P0P4P1P2) + V (P0P1P2P3) = V (P1P2P3P4) (2.70)

NASA/TM—2013-218061 16
On the other hand, for the case where P0 lies outside of the tetrahedron P1P2P3P4, Eq. (2.68) is no longer

true, i.e., some of the weight factors must be negative. Because Eq. (2.51) requires that some of the weight

factors must be positive, we arrive at the conclusion that

|w1| + |w2| + |w3| + |w4| > |w1 + w2 + w3 + w4| = 1 (2.71)

if P0 lies outside of the tetrahedron P1P2P3P4. In turn, Eqs. (2.67) and (2.71) ⇒

V (P0P2P3P4) + V (P0P3P4P1) + V (P0P4P1P2) + V (P0P1P2P3) > V (P1P2P3P4) (2.72)

if P0 lies outside of the tetrahedron P1P2P3P4.

Note that, for any real parameter ζ, Eqs. (2.42) and (2.47) ⇒

Ak > 0 and (Ak)

ζ > 0 k = 1, 2, 3, 4; −∞ < ζ < +∞ (2.73)

As such the weight factors

wk(ζ)

def =

(Ak)

ζ

(A1)

ζ + (A2)ζ + (A3)ζ + (A4)ζ

k = 1, 2, 3, 4; −∞ < ζ < +∞ (2.74)

are well defined and satisfy the conditions

0 < wk(ζ) < 1 k = 1, 2, 3, 4 − ∞ < ζ < +∞ (2.75)

and

w1(ζ) + w2(ζ) + w3(ζ) + w4(ζ) = 1 − ∞ < ζ < +∞ (2.76)

Note that, even though the same notation is used, the weight factors wk(ζ), k = 1, 2, 3, 4, defined in Eq. (2.74)

are different from those defined in Eq. (2.29). The definition Eq. (2.74) will be used in this subsection.

Moreover, hereafter we consider only the special case where

wk = wk(ζ) k = 1, 2, 3, 4; −∞ < ζ < +∞ (2.77)

For this special case, Eqs. (2.75)–(2.77) imply that both Eqs. (2.51) and are satisfied for any real

number ζ. Thus, for the case Eq. (2.77), (i) point P0 defined by Eq. (2.50) must always lie in the interior of

the tetrahedron P1P2P3P4, and (ii) Eqs. (2.63)–(2.67) reduce to

V (P0P2P3P4) = w1 × V (P1P2P3P4) (2.78)

V (P0P3P4P1) = w2 × V (P1P2P3P4) (2.79)

V (P0P4P1P2) = w3 × V (P1P2P3P4) (2.80)

V (P0P1P2P3) = w4 × V (P1P2P3P4) (2.81)

and Eq. (2.70), respectively. Moreover, for the special case ζ = 0, we have

w1 = w2 = w3 = w4 = 1/4 and −→P 0 =

−→P 1 +

−→P 2 +−→P 3 +−→P 4



/4 (ζ = 0) (2.82)

Because the position vector of the centroid of a tetrahedron is the simple average of those of its four vertices,

Eq. (2.82) ⇒ point P0 is the centroid of the tetrahedron P1P2P3P4 if ζ = 0.

Note that, for a polyhedron other than a tetrahedron (which is the simplest polyhedron), generally it is

not true that the position vector of its centroid is the simple average of those of its vertices. As an example,

NASA/TM—2013-218061 17

(2.69)
the position vector of the centroid of a hexahedron may not be the simple average of those of its eight vertices

unless the hexahedron is a parallelepiped.

Next, with the aid of Eq. (2.47) and (2.48), and Fig. 6(a), one can see that

V (P0P2P3P4) = h1A1

3

, V (P0P3P4P1) = h2A2

3

, V (P0P4P1P2) = h3A3

3

and V (P0P1P2P3) = h4A4

3

(2.83)

In turn, by substituting Eq. (2.83) into Eqs. (2.78)–(2.81) and using Eqs. (2.74) and (2.77), one has

hkAk

3

= wkV (P1P2P3P4) = (Ak)

ζ × V (P1P2P3P4)

(A1)

ζ + (A2)ζ + (A3)ζ + (A4)ζ

k = 1, 2, 3, 4; −∞ < ζ < +∞ (2.84)

and thus

h1

(A1)

ζ−1

=

h2

(A2)

ζ−1

=

h3

(A3)

ζ−1

=

h4

(A4)

ζ−1

=

3 × V (P1P2P3P4)

(A1)

ζ + (A2)ζ + (A3)ζ + (A4)ζ

− ∞ < ζ < +∞ (2.85)

Because (i) h1, h2, h3, and h4 are the distances from point P0 to the triangular faces 4P2P3P4, 4P3P4P1,

4P4P1P2, and 4P1P2P3 respectively, and (ii) A1, A2, A3, and A4 are the areas of 4P2P3P4, 4P3P4P1,

4P4P1P2, and 4P1P2P3, respectively, Eq. (2.85) states that, for a given tetrahedron P1P2P3P4 and a given

ζ, the distance from point P0 to any face of this tetrahedron is proportional to the (ζ − 1)th power of the

area of the face.

In fact, for the case ζ = 0 (i.e., when P0 is the centroid of the tetrahedron P1P2P3P4), Eq. (2.85) reduces

to

hk =



3 × V (P1P2P3P4)

4



1

Ak

k = 1, 2, 3, 4, if ζ = 0 (2.86)

i.e., the distance between the centroid P0 and a face of the tetrahedron P1P2P3P4 is inversely proportional

to the area of this face. Thus the centroid is farthest away from the smallest face and closest to the largest

face.

On the other hand, for the case ζ = 1, Eq. (2.85) reduces to

h1 = h2 = h3 = h4 =

3 × V (P1P2P3P4)

A1 + A2 + A3 + A4

if ζ = 1 (2.87)

i.e., point P0 has the same distance to the four faces of the tetrahedron P1P2P3P4 if ζ = 1. By definition,

this implies that P0 is the in-center of the tetrahedron P1P2P3P4 for the case ζ = 1.

In general one can infer from Eq. (2.85) that:

(a) For the case ζ > 1, point P0 is farthest away from the largest face of the tetrahedron P1P2P3P4 and

closest to its smallest face. The trend becomes stronger as the value of ζ increases from 1.

(b) For the case ζ < 1, point P0 is farthest away from the smallest face of the tetrahedron P1P2P3P4 and

closest to its largest face. The trend becomes stronger as the value of ζ decreases from 1.

Let the interior point P0 of any tetrahedron P1P2P3P4 defined by Eqs. (2.50) and (2.77) be referred to

as the interior point of this tetrahedron associated with the parameter ζ. Then, from the above observations

and the discussions given in Sec. 1, it becomes clear that a spatial solution point which is the centroid of

the 8-vertex polyhedron referred to earlier will be much less likely to lie outside of the associated 8-vertex

polyhedron and causes numerical difficulties if each of the centroidal mesh points P0, A0, B0, C0 and D0

depicted in Fig. 6(b) is replaced by an interior point associated with a parameter ζ ≥ 1. In fact, it has been

shown through numerical experiments that the numerical difficulties associated with spatial solution points

being embedded beyond a solid-wall boundary can be overcome by replacing centroids of tetrahedra with

in-centers (i.e., interior points associated with ζ = 1) in the construction of a spatial tetrahedral mesh.

NASA/TM—2013-218061 18
In this section, it will be described how an 1D CESE scheme can be modified so that the two decoupled solutions of the original scheme can become coupled in the new scheme while still preserving a flux

conservation condition to be defined shortly. This description is then followed by a sketch of the 2D and 3D

extensions.

3.1. 1D CESE solution-coupling schemes

As a preliminary, note that the CEs, SEs and the 1D schemes, which were defined over Ω1 in Sec. 1,

hereafter will be defined over

Ω

def = {(j, n)|j, n = 0, ±1/2, ±1, ±3/2, . . .} = Ω1 ∪ Ω2 (3.1)

In particular, it should be understood that the symbol Ω1 from now on is replaced by Ω in each of Eqs. (1.3)–

(1.16). Also, as an example, the extension (defined over Ω) of the a- scheme will be referred to as the dual

a- scheme.

To clarify the significance of the extension, we offer the following observations:

(a) The solution of a dual scheme is formed by two decoupled solutions which are defined over Ω1 and Ω2,

respectively.

(b) The same space-time region can be designated as a BCE of a mesh point ∈ Ω1 and also as a BCE of

another mesh point ∈ Ω2. For examples, the space-time region ABCD depicted in Fig. 2(c) can be

designated as CE−(j, n) and CE+(j−1/2, n), respectively, while the region AF ED depicted in Fig. 2(d)

as CE+(j, n) and CE−(j + 1/2, n), respectively. Hereafter, as an example, points (j, n) and (j − 1/2, n)

(i.e., points A and B) will be referred to as the hosts of the region ABCD.

(c) Let the region ABCD be designated as CE−(j, n). Then, on its boundary, the line segments AB and

AD belong to SE(j, n) while CB and CD belong to SE(j − 1/2, n − 1/2), i.e., over AB and AD, the

vector ~h

∗

should be evaluated using the solution data stored at point A while, over CB and CD, it

should be evaluated using those stored at point C.

Alternately, let the region ABCD be designated as CE+(j − 1/2, n). Then, on its boundary, the line

segments BA and BC belong to SE(j −1/2, n) while DA and DC belong to SE(j, n−1/2), i.e., over BA

and BC, the vector ~h

∗

should be evaluated using the solution data stored at point B while, over DA

and DC, it should be evaluated using those stored at point D. Note that, in the above and hereafter,

as an example, a line segment joining points A and B is denoted by AB if it belongs to SE(A) (i.e., the

SE centered at point A) while the same line segment is denoted by BA if it belongs SE(B) (i.e., the SE

centered at point B).

(d) From the discussions given in the above item (c), it becomes clear that

I

S(CE+(j,n))

~h

∗

· d~s = 0 and I

S(CE−(j+1/2,n))

~h

∗

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

represent two different flux conservation conditions even though CE+(j, n) and CE−(j + 1/2, n) occupy

the same space-time region for any (j, n) ∈ Ω.

As will be shown, each marching step of the 1D solution-coupling scheme is formed by two sub-steps.

In the first sub-step, referred to as the dual-scheme marching step, the intermediate mesh values at the

nth time level are evaluated from the given mesh values at the (n − 1/2)th time level using a dual-scheme

(see Fig. 7). In the second sub-step, referred to as the solution-coupling step, the final mesh values at the

nth time level are updated from the intermediate mesh values through a flux-conserving weighted-averaging

procedure. Let the intermediate mesh values of u and its normalized spatial derivative ux ̄ at any (j, n) ∈ Ω

be denoted by ˆu

n

j

and (ˆux ̄)

n

j

, respectively. Then, as an example, ˆu

n

j

and (ˆux ̄)

n

j

can be evaluated in terms of

the given mesh values at the (n − 1/2)th time level by using the dual a- scheme defined by Eqs. (1.12) and

(1.15), i.e.,

uˆ

n

j =

1

2

n

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

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω (3.3)

3. Flux-conserving solution-coupling procedures

NASA/TM—2013-218061 19
and

(ˆux ̄)

n

j =

1

2

n

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

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω (3.4)

Note that, by using a relation similar to the second equation in Eq. (1.9), one can define the “unnormalized”

version of (ˆux ̄)

n

j

, i.e.,

(ˆux)

n

j

def =

4

∆x

(ˆux ̄)

n

j

(j, n) ∈ Ω (3.5)

To avoid introducing unnecessary new notations, in the following description of the solution-coupling

step, at each (j, n) ∈ Ω, the final mesh values after updating from the intermediate values are simply denoted

as u

n

j

and (ux ̄)

n

j

, respectively.

Note that the line segment AF depicted in Fig. 7 belongs to SE(j, n). Thus, by using (i) Eqs. (1.4),

(1.6) and (1.9), and (ii) Fig. 2(d), one concludes that the flux of ~h

∗

leaving CE+(j, n) through AF evaluated

using the final solution data stored at point A is given by

F(AF)

def =

Z xj+(∆x/2)

xj

u

∗

(x, tn; j, n) dx =

∆x

2

(u+)

n

j

(j, n) ∈ Ω (3.6)

where

(u+)

n

j

def = (u + ux ̄)

n

j

(j, n) ∈ Ω (3.7)

Let Fˆ(AF) be the intermediate value of F(AF), i.e., the value of F(AF) when u

n

j

and (ux ̄)

n

j

are replaced

by ˆu

n

j

and (ˆux ̄)

n

j

, respectively. Then Eqs. (3.5)–(3.7) ⇒

Fˆ(AF) = ∆x

2

(ˆu+)

n

j

(j, n) ∈ Ω (3.8)

where

(ˆu+)

n

j

def = (ˆu + ˆux ̄)

n

j

(j, n) ∈ Ω (3.9)

At this juncture, note that: (i) at the midpoint M+ of AF, x = xj + (∆x/4) and t = t

n; and (ii)

Eqs. (1.6), (1.9), and (3.7) ⇒

u

∗

(xj + (∆x/4), tn; j, n) = (u+)

n

j

(j, n) ∈ Ω (3.10)

Thus, (u+)

n

j

is the value of u at the midpoint M+ of AF evaluated using the values u

n

j

and (ux ̄)

n

j with

the aid of Eqs. (1.6) and (1.9). Similarly (ˆu+)

n

j

is the value of u at the midpoint M+ evaluated using the

intermediate values uˆ

n

j

and (ˆux ̄)

n

j with the aid of Eqs. (1.6), (1.9), and (3.5).

Because (u+)

n

j

and (ˆu+)

n

j

are evaluated using the solution data stored at point (j, n) (i.e., point A),

hereafter, respectively, they will be referred to as the final and intermediate values of u at the midpoint

M+ associated with the host point A of the region AF ED. Similarly, F(AF) and Fˆ(AF), respectively, are

referred to as the final and intermediate values of the flux leaving the top face of the same region associated

with the same host.

Moreover, with the aid of Eqs. (3.6) and (3.8), and Fig. 7, one concludes that F(AF) (Fˆ(AF)) is equal

to the length of AF multiplied by the approximated value (u+)

n

j

((ˆu+)

n

j

) of u at the midpoint M+. In other

words, the value of u at the midpoint M+ is the average value of u over AF—a result stemmed from the

fact that u

∗

(x, t; j, n) is a first-order Taylor’s expansion in x when t = t

n (see Eqs. (1.3) and (1.6)).

Similarly, with the aid of Fig. 7, one concludes that:

(a) The flux of ~h

∗

leaving CE−(j, n) through the line segment AB evaluated using the final solution data

stored at point A is given by

F(AB)

def =

Z xj−(∆x/2)

xj

u

∗

(x, tn; j, n) dx =

∆x

2

(u−)

n

j

(j, n) ∈ Ω (3.11)

NASA/TM—2013-218061 20
where

(u−)

n

j

def = (u − ux ̄)

n

j

(j, n) ∈ Ω (3.12)

is the value of u at the midpoint M− of AB evaluated using the values u

n

j

and (ux ̄)

n

j with the aid of

Eqs. (1.6) and (1.9).

Also Fˆ(AB), the intermediate value of F(AB), is given by

Fˆ(AB) = ∆x

2

(ˆu−)

n

j

(j, n) ∈ Ω (3.13)

where

(ˆu−)

n

j

def = (ˆu − uˆx ̄)

n

j

(j, n) ∈ Ω (3.14)

is the value of u at the midpoint M− evaluated using the intermediate values uˆ

n

j

and (ˆux ̄)

n

j with the aid

of Eqs. (1.6), (1.9), and (3.5).

Because (u−)

n

j

and (ˆu−)

n

j

are evaluated using the solution data stored at point (j, n) (i.e., point A),

hereafter, respectively, they will be referred to as the final and intermediate values of u at the midpoint

M− associated with the host point A of the region ABCD. Similarly, F(AB) and Fˆ(AB), respectively,

are referred to as the final and intermediate values of the flux leaving the top face of the same region

associated with the same host.

Moreover, with the aid of Eqs. (3.11) and (3.13), one concludes that F(AB) (Fˆ(AB)) is equal to the

length of AB multiplied by the approximated value (u−)

n

j

((ˆu−)

n

j

) of u at the midpoint M−.

(b) The flux of ~h

∗

leaving CE−(j + 1/2, n) through the line segment F A evaluated using the final solution

data stored at point F is given by

F(F A)

def =

Z xj+(∆x/2)

xj

u

∗

(x, tn; j + 1/2, n) dx =

∆x

2

(u−)

n

j+1/2

(j, n) ∈ Ω (3.15)

where

(u−)

n

j+1/2

def = (u − ux ̄)

n

j+1/2

(j, n) ∈ Ω (3.16)

is the value of u at the midpoint M+ of F A evaluated using the values u

n

j+1/2

and (ux ̄)

n

j+1/2 with the aid

of Eqs. (1.6) and (1.9). Note that: (i) because (j, n) ∈ Ω ⇔ (j + 1/2, n) ∈ Ω, Eq. (3.12) ⇔ Eq. (3.16);

and (ii) according to Eqs. (3.6) and (3.15), F(AF) 6= F(F A) even though AF and F A denote the same

line segment. In fact, F(AF) is a flux evaluated using the solution data stored at point A while F(F A)

is another flux evaluated using those stored at point F. The reader is warned that other similar pairs

of notations will appear throughout the paper.

Also Fˆ(F A), the intermediate value of F(F A), is given by

Fˆ(F A) = ∆x

2

(ˆu−)

n

j+1/2

(j, n) ∈ Ω (3.17)

where

(ˆu−)

n

j+1/2

def = (ˆu − uˆx ̄)

n

j+1/2

(j, n) ∈ Ω (3.18)

is the value of u at the midpoint M+ evaluated using the intermediate values uˆ

n

j+1/2

and (ˆux ̄)

n

j+1/2 with

the aid of Eqs. (1.6), (1.9), and (3.5). Obviously Eq. (3.14) ⇔ Eq. (3.18).

Because (u−)

n

j+1/2

and (ˆu−)

n

j+1/2

are evaluated using the solution data stored at point (j + 1/2, n) (i.e.,

point F), hereafter, respectively, they will be referred to as the final and intermediate values of u at the

midpoint M+ associated with the host point F of the region AF ED. Similarly, F(F A) and Fˆ(F A),

respectively, are referred to as the final and intermediate values of the flux leaving the top face of the

same region associated with the same host.

NASA/TM—2013-218061 21
Moreover, with the aid of Eqs. (3.15) and (3.17), one concludes that F(F A) (Fˆ(F A)) is equal to the

length of F A multiplied by the approximated value (u−)

n

j+1/2

((ˆu−)

n

j+1/2

) of u at the midpoint M+.

(c) The flux of ~h

∗

leaving CE+(j − 1/2, n) through the line segment BA evaluated using the final solution

data stored at point B is given by

F(BA)

def =

Z xj

xj−(∆x/2)

u

∗

(x, tn; j − 1/2, n) dx =

∆x

2

(u+)

n

j−1/2

(j, n) ∈ Ω (3.19)

where

(u+)

n

j−1/2

def = (u + ux ̄)

n

j−1/2

(j, n) ∈ Ω (3.20)

is the value of u at the midpoint M− of BA evaluated using the values u

n

j−1/2

and (ux ̄)

n

j−1/2 with the

aid of Eqs. (1.6) and (1.9). Obviously, Eq. (3.7) ⇔ Eq. (3.20).

Also Fˆ(BA), the intermediate value of F(BA), is given by

Fˆ(BA) = ∆x

2

(ˆu+)

n

j−1/2

(j, n) ∈ Ω (3.21)

where

(ˆu+)

n

j−1/2

def = (ˆu + ˆux ̄)

n

j−1/2

(j, n) ∈ Ω (3.22)

is the value of u at the midpoint M− evaluated using the intermediate values uˆ

n

j−1/2

and (ˆux ̄)

n

j−1/2 with

the aid of Eqs. (1.6), (1.9), and (3.5). Obviously Eq. (3.9) ⇔ Eq. (3.22).

Because (u+)

n

j−1/2

and (ˆu+)

n

j−1/2

are evaluated using the solution data stored at point (j −1/2, n) (i.e.,

point B), hereafter, respectively, they will be referred to as the final and intermediate values of u at the

midpoint M− associated with the host point B of the region ABCD. Similarly, F(BA) and Fˆ(BA),

respectively, are referred to as the final and intermediate values of the flux leaving the top face of the

same region associated with the same host.

Moreover, with the aid of Eqs. (3.19) and (3.21), one concludes that F(BA) (Fˆ(BA)) is equal to the

length of BA multiplied by the approximated value (u+)

n

j−1/2

((ˆu+)

n

j−1/2

) of u at the midpoint M−.

Next recall that two different final fluxes F(AF) and F(F A), respectively associated with the hosts

points A and F of the region AF ED, are defined over the top face of this region. For a reason explained in

[42], the dual flux (denoted by F(A; F) or F(F; A)) associated with these two fluxes is defined as the simple

average of them, i.e.,

F(A; F) ≡ F(F; A)

def =



F(AF) + F(F A)



/2 (3.23)

Similarly, we have

Fˆ(A; F) ≡ Fˆ(F; A)

def =

h

Fˆ(AF) + Fˆ(F A)

i

/2 (3.24)

The solution-coupling step will be constructed so that, at the top face of any BCE (such as the space-time

region AF ED) sandwiched between the nth and (n − 1/2)th time levels, the dual flux evaluated using u

n

j

and (ux ̄)

n

j

, (j, n) ∈ Ω, is equal to that evaluated using ˆu

n

j

and (ˆux ̄)

n

j

, (j, n) ∈ Ω, i.e.,

F(A; F) = Fˆ(A; F) (j, n) ∈ Ω (3.25)

By substituting Eqs. (3.6), (3.8), (3.15), and (3.17) into Eqs. (3.23) and (3.24), and then imposing the

condition Eq. (3.25), one has

(u+)

n

j + (u−)

n

j+1/2 = (ˆu+)

n

j + (ˆu−)

n

j+1/2

(j, n) ∈ Ω (3.26)

Recall that: (i) (u+)

n

j

and (u−)

n

j+1/2

, respectively, are the final values of u at the midpoint M+ evaluated

using the solution data stored at the two hosts (points A and F) of the region AF ED; and (ii) (ˆu+)

n

j

and

NASA/TM—2013-218061 22
(ˆu−)

n

j+1/2

, respectively, are the intermediate-value counterparts of (u+)

n

j

and (u−)

n

j+1/2

. As such, Eq. (3.26)

states that the average of the two final values of u at the midpoint M+ evaluated, respectively, using the

solution data stored at the two hosts of the region AF ED, is identical to its intermediate-value counterpart.

Let

(u+)

n

j = (1 − η)(ˆu+)

n

j + η (ˆu−)

n

j+1/2

(j, n) ∈ Ω (3.27)

where η is a real adjustable parameter. Then (u+)

n

j

is an weighted average of (ˆu+)

n

j

and (ˆu−)

n

j+1/2

. In turn,

by substituting Eq. (3.27) into Eq. (3.26), one has

(u−)

n

j+1/2 = η (ˆu+)

n

j + (1 − η)(ˆu−)

n

j+1/2

(j, n) ∈ Ω (3.28)

i.e., (u−)

n

j+1/2

is also an weighted average of (ˆu+)

n

j

and (ˆu−)

n

j+1/2

. In fact, Eqs. (3.27) and (3.28) ⇒ (i)

(u+)

n

j = (ˆu+)

n

j and (u−)

n

j+1/2 = (ˆu−)

n

j+1/2

if η = 0 (j, n) ∈ Ω (3.29)

(ii)

(u+)

n

j = (u−)

n

j+1/2 =

1

2

h

(ˆu+)

n

j + (ˆu−)

n

j+1/2

i

if η = 1/2 (j, n) ∈ Ω (3.30)

and (iii)

(u+)

n

j = (ˆu−)

n

j+1/2

and (u−)

n

j+1/2 = (ˆu+)

n

j

if η = 1 (j, n) ∈ Ω (3.31)

Because (j, n) ∈ Ω ⇔ (j − 1/2, n) ∈ Ω, Eq. (3.27)–(3.31) remain valid if each index j in these equations

is replaced by j − 1/2. In particular, Eq. (3.28) ⇔

(u−)

n

j = η (ˆu+)

n

j−1/2 + (1 − η)(ˆu−)

n

j

(j, n) ∈ Ω (3.32)

In turn, by substituting Eqs. (3.7), (3.9), (3.12), (3.14), (3.18), and (3.22) into Eqs.(3.27) and (3.32), one

has

(u + ux ̄)

n

j = (1 − η)(ˆu + ˆux ̄)

n

j + η(ˆu − uˆx ̄)

n

j+1/2

(j, n) ∈ Ω (3.33)

(u − ux ̄)

n

j = (1 − η)(ˆu − uˆx ̄)

n

j + η(ˆu + ˆux ̄)

n

j−1/2

(j, n) ∈ Ω (3.34)

By adding Eqs. (3.33) and (3.34) together, one has

u

n

j = (1 − η)ˆu

n

j + (η/2) h

(ˆu − uˆx ̄)

n

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

n

j−1/2

i

(j, n) ∈ Ω (3.35)

On the other hand, by subtracting Eq. (3.34) from Eq. (3.33), one has

(ux ̄)

n

j = (1 − η)(ˆux ̄)

n

j + (η/2) h

(ˆu − uˆx ̄)

n

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

n

j−1/2

i

(j, n) ∈ Ω (3.36)

Eqs. (3.35) and (3.36) form the solution coupling step, i.e., with the aid of these two equations, at each mesh

point (j, n) ∈ Ω, u

n

j

and (ux ̄)

n

j

can be determined in terms of the intermediate mesh values at the mesh

points (j − 1/2, n), (j, n), and (j + 1/2, n).

As expected from Eq. (3.29), we have

u

n

j = ˆu

n

j and (ux ̄)

n

j = (ˆux ̄)

n

j

if η = 0 (j, n) ∈ Ω (3.37)

i.e., solution coupling does not occur when η = 0. On the other hand, because of Eq. (3.30), the case with

η = 1/2 is referred as the full coupling case.

NASA/TM—2013-218061 23
Let the dual-scheme marching step be executed using the a- scheme, i.e., Eqs. (3.3) and (3.4) are

assumed. Then, by substituting Eqs. (3.3) and (3.4) into Eqs. (3.35) and (3.36), one obtains

u

n

j =

η

2

[(2 − )u − νux ̄]

n−1/2

j

+

(1 − η)

2

n

(1 + ν)u + (1 − ν

2

)ux ̄

n−1/2

j−1/2

+



(1 − ν)u − (1 − ν

2

)ux ̄

n−1/2

j+1/2

o

+

η

4

n

( + ν)u + (2 + ν − ν

2

)ux ̄

n−1/2

j−1

+



( − ν)u + (ν + ν

2 − 2)ux ̄

n−1/2

j+1 o

(j, n) ∈ Ω (3.38)

and

(ux ̄)

n

j =

η

2



νu + (2 − 2 − ν

2

)ux ̄

n−1/2

j

+

(1 − η)

2

n

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

n−1/2

j−1/2

+



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

n−1/2

j+1/2

o

−

η

4

n

( + ν)u + (2 + ν − ν

2

)ux ̄

n−1/2

j−1

−



( − ν)u + (ν + ν

2 − 2)ux ̄

n−1/2

j+1 o

(j, n) ∈ Ω (3.39)

Because Eqs. (3.38) and (3.39) are derived assuming that (i) the dual-mesh marching step is carried out using

the dual a- scheme, and (ii) the solution-coupling step is formed using an weighted-averaging procedure

involving the parameter η (see Eq. (3.27)), hereafter the full solution-coupling scheme formed by Eqs. (3.38)

and (3.39) will be referred to as the a--η scheme.

Alternatively, the dual marching step can be executed using the dual c-τ scheme which are defined

by Eqs. (1.12) and (1.16). For this case, the intermediate values ˆu

n

j

and (ˆux ̄)

n

j

should be evaluated using

Eq. (3.3) and

(ˆux ̄)

n

j =

1

2(1 + τ )

n

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

n−1/2

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

n−1/2

j−1/2

o

(j, n) ∈ Ω (3.40)

By substituting Eqs. (3.3) and (3.40) into Eqs. (3.35) and (3.36), one obtains the c-τ -η scheme, i.e.,

u

n

j =

η

2(1 + τ )

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

n−1/2

j

+

(1 − η)

2

n

(1 + ν)u + (1 − ν

2

)ux ̄

n−1/2

j−1/2

+



(1 − ν)u − (1 − ν

2

)ux ̄

n−1/2

j+1/2

o

+

η

4(1 + τ )

n

(τ + ν + τν)u + [2(ν + τ ) − (1 + τ )ν

2

]ux ̄

n−1/2

j−1

+



(τ − ν − τν)u + [2(ν − τ ) + (1 + τ )ν

2

]ux ̄

n−1/2

j+1 o

(j, n) ∈ Ω (3.41)

and

(ux ̄)

n

j =

η

2(1 + τ )



(1 + τ )νu + [2 − (1 + τ )ν

2

]ux ̄

n−1/2

j

−

(1 − η)

2(1 + τ )

n

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

n−1/2

j−1/2

−



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

n−1/2

j+1/2

o

−

η

4(1 + τ )

n

(τ + ν + τν)u + [2(ν + τ ) − (1 + τ )ν

2

]ux ̄

n−1/2

j−1

−



(τ − ν − τν)u + [2(ν − τ ) + (1 + τ )ν

2

]ux ̄

n−1/2

j+1 o

(j, n) ∈ Ω (3.42)

NASA/TM—2013-218061 24
The von Neumann stability analysis was carried out for both a--η and c-τ -η scheme. In particular,

the amplification matrices and the associated amplification factors were obtained for both schemes involving

arbitrary values of , τ , η, and the phase angle θ [72]. Because the a--η and c-τ -η schemes reduce to

the dual a- and c-τ schemes, respectively when η = 0, as expected, the stability conditions of these two

schemes also reduce to those of the a- and c-τ schemes, respectively when η = 0. For the general case

involving an arbitrary value of η, it turns out that the stability conditions of these schemes generally become

so complicated that they can only be expressed numerically. As such, the discussion of their stability

conditions for the general case is beyond the scope of the current paper.

However, for the full coupling case, i.e., when η = 1/2, it can be shown that [72] the a--η scheme is von

Neumann stable ⇔

ν

4 ≤  ≤ 1 (the a--η scheme with η = 1/2) (3.43)

Let

 = β |ν| β > 0 (3.44)

where β > 0 is an adjustable parameter. Then Eq. (3.43) becomes

ν

4 ≤ β |ν| ≤ 1 (3.45)

Excluding the meaningless case ν = 0, Eq. (3.45) ⇔ (|ν|

3 − β) ≤ 0 and |ν| ≤ 1/β, i.e.,

|ν| ≤ min{β

1/3

, 1/β} =







1 if β = 1

β

1/3 < 1 if 0 < β < 1

1/β < 1 if β > 1

(ν 6= 0) (3.46)

In turn, Eqs. (3.44) and (3.46) ⇒

 = β |ν| =

(1 if β ≥ 1

β

4/3 < 1 if 0 < β < 1

(ν 6= 0) (3.47)

Also, for the case η = 1/2, it has been shown that the stability conditions of the c-τ -η scheme is

f(ν

2

) ≤ τ (ν

2 ≤ 2) (3.48)

if ν

2 ≤ 2, where

f(s)

def =

s(s + 2)

1 + 2s − s

2

(0 ≤ s ≤ 2) (3.49)

Note that (i)

1 + 2s − s

2 = −(s +

√

2 − 1)(s −

√

2 − 1) > 0 if 0 ≤ s ≤ 2 (3.50)

i.e., the denominator of the expression on the right side of Eq. (3.49) is not zero throughout the domain of

f(s) and therefore the function is well defined over its domain; and (ii)

df(s)

ds =

2(2s

2 + s + 1)

(1 + 2s − s

2)2

> 0 if 0 ≤ s ≤ 2 (3.51)

i.e., f(s) is monotonically increasing in its domain and therefore

0 = f(0) ≤ f(s) ≤ f(2) = 8 (0 ≤ s ≤ 2) (3.52)

NASA/TM—2013-218061 25
The stability conditions of the c-τ -η scheme for the case η = 1/2 and ν

2 > 2 can also be cast into an

analytical form. However, the form is much more complicated and will not be presented here.

3.2. 2D and 3D extensions

For a 2D case using triangular meshes, or a 3D case using tetrahedral meshes, there is a natural extension

of each of the 1D a--η and c-τ -η schemes. As an example, consider Fig. 8 where the triangle 4BDF and

its three neighbors 4BF H, 4BID, and 4DJF (which all lie on the x-y plane) are shown. As explained

in Sec. 2, each of the centroid mesh points G, A, C, and E shown in Fig. 3 is replaced in Fig. 8 by another

associated interior point which is (i) defined by Eqs. (2.10) and (2.29) with ζ ≥ 1; and (ii) denoted by the

same symbol and also marked by an open circle. The spatial solution point G0(also marked by a cross)

associated with 4BDF again is defined as the centroid of the hexagon ABCDEF. Other spatial solution

points A0, C

0

, and E0 are defined similarly. In addition, points K, L, and M (marked by triangles) depicted

in Fig. 8, by definition, are the centroids of the quadrilaterals GF AB, GBCD, and GDEF, respectively.

Let each spatial solution point be assigned an unique identification index (an integer). As an example,

let points G0, A0, C

0

, and E0 be identified by the indices j, j1, j2, and j3, respectively. As such, as an

example, the x- and y- coordinates of point G0 will be denoted by xj and yj , respectively. Moreover, the

space-time solution points (G0, n), (A0, n), (C

0

, n), and (E0, n), which lie on the nth time level and have

the points G0, A0, C

0

, and E0 being their spatial projections, respectively, will be denoted by (j, n), (j1, n),

(j2, n), and (j3, n), respectively. Also each space-time solution point (j, n) is associated with [13] (i) one SE

denoted by SE(j, n); and (ii) three BCEs denoted by CE(k)(j, n), k = 1, 2, 3, respectively, and a CCE which

is the union of CE(k)(j, n), k = 1, 2, 3, and denoted by CE(j, n). By definition, CE(1)(j, n) is the space-time

cylinder sandwiched between the nth and (n − 1/2)th time levels with the quadrilateral GF AB being its

spatial projection. CE(2)(j, n) and CE(3)(j, n) are defined identically except that the quadrilateral GF AB

are replaced by the quadrilaterals GBCD and GDEF, respectively. It is seen that (i) the solution points

(j, n) and (j1, n) are the hosts of the space-time region occupied by CE(1)(j, n), (ii) the solution points (j, n)

and (j2, n) are the hosts of the space-time region occupied by CE(2)(j, n), and (iii) the solution points (j, n)

and (j3, n) are the hosts of the space-time region occupied by CE(3)(j, n),

Let the 2D counterpart of Eq. (1.1) be

∂u

∂t + ax

∂u

∂x + ay

∂u

∂y = 0 (3.53)

where ax, and ay are constants. Then the 2D counterparts of Eqs. (1.4) and (1.6) are

~h

∗

(x, y, t; j, n)

def =



axu

∗

(x, y, t; j, n), ayu

∗

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



(3.54)

and

u

∗

(x, y, t; j, n) = u

n

j + (ux)

n

j

[(x − xj ) − ax(t − t

n

)] + (uy)

n

j

[(y − yj ) − ay(t − t

n

)] (3.55)

respectively. Thus there are three independent marching variables, i.e., u

n

j

, (ux)

n

j

, and (uy)

n

j

associated with

a space-time solution point (j, n) (i.e., point (G0, n)).

Next let (i) (K, n), (L, n), and (M, n) denote the points on the nth time level with their spatial projections being points K, L, and M, respectively; (ii) the intermediate solution data stored at (j, n), (j1, n),

(j2, n), and (j3, n) be evaluated from the given solution data at the (n − 1/2)th time level using the 2D

versions of the a- or c-τ scheme [7,13,31]; (iii) [ˆu(K, n)]n

j

and [ˆu(K, n)]n

j1

denote the intermediate values of u

at (K, n) evaluated using the intermediate solution data stored at the hosts (j, n) and (j1, n) of CE(1)(j, n),

respectively; (iv) [ˆu(L, n)]n

j

and [ˆu(L, n)]n

j2

denote the intermediate values of u at (L, n) evaluated using the

intermediate solution data stored at the hosts (j, n) and (j2, n) of CE(2)(j, n), respectively; (v) [ˆu(M, n)]n

j

and [ˆu(M, n)]n

j3

denote the intermediate values of u at (L, n) evaluated using the intermediate solution data

stored at the hosts (j, n) and (j3, n) of CE(3)(j, n), respectively; (vi) [u(K, n)]n

j

and [u(K, n)]n

j1

denote the

final values of u at (K, n) evaluated using the final solution data stored at the hosts (j, n) and (j1, n) of

CE(1)(j, n), respectively; (vii) [u(L, n)]n

j

and [u(L, n)]n

j2

denote the final values of u at (L, n) evaluated using

NASA/TM—2013-218061 26
the final solution data stored at the hosts (j, n) and (j2, n) of CE(2)(j, n), respectively; and (viii) [u(M, n)]n

j

and [u(M, n)]n

j3

denote the final values of u at (M, n) evaluated using the final solution data stored at the

hosts (j, n) and (j3, n) of CE(3)(j, n), respectively.

Moreover, for each k = 1, 2, 3, let (i) Ak(j, n) denote the area of the spatial projection of CE(k)(j, n);

(ii) Fˆ

k(j, n) denote the intermediate value of the flux of ~h

∗

leaving CE(k)(j, n) through its top face evaluated

using the intermediate solution data stored at the space-time solution point (j, n); (iii) Fˆ

k(jk, n) denote the

intermediate value of the flux of ~h

∗

leaving CE(k)(j, n) through its top face evaluated using the intermediate

solution data stored at the space-time solution point (jk, n); (iv) Fk(j, n) denote the final value of the flux

of ~h

∗

leaving CE(k)(j, n) through its top face evaluated using the final solution data stored at the space-time

solution point (j, n); and (v) Fk(jk, n) denote the final value of the flux of ~h

∗

leaving CE(k)(j, n) through its

top face evaluated using the final solution data stored at the space-time solution point (jk, n). According

to the above definitions, A1(j, n), A2(j, n), and A3(j, n) are the areas of the quadrilaterals GF AB, GBCD,

and GDEF, respectively. Because these quadrilaterals are the spatial projections of CE(k)(j, n), k = 1, 2, 3,

respectively, for each k = 1, 2, 3, Ak(j, n) is the area of the top face of CE(k)(j, n).

With the aid of (i) the above definitions, (ii) the fact that point K depicted in Fig. 8 is the centroid

of the quadrilateral GF AB (i.e., the spatial projection of CE(1)(j, n)), and (iii) Eq. (3.55) is a first-order

Taylor’s expansion, by using Eqs. (3.54) and (3.55), one can establish the following relations for several fluxes

of ~h

∗

leaving CE(1)(j, n) through its top face:

Fˆ

1(j, n) = [ˆu(K, n)]n

j A1(j, n) (3.56)

Fˆ

1(j1, n) = [ˆu(K, n)]n

j1A1(j, n) (3.57)

F1(j, n) = [u(K, n)]n

j A1(j, n) (3.58)

and

F1(j1, n) = [u(K, n)]n

j1A1(j, n) (3.59)

The current solution-coupling step will be constructed assuming the following dual-flux conservation relation:

F1(j, n) + F1(j1, n) = Fˆ

1(j, n) + Fˆ1(j1, n) (3.60)

By substituting Eqs. (3.56)–(3.59) into Eq. (3.60), one obtains

[u(K, n)]n

j + [u(K, n)]nj1 = [ˆu(K, n)]nj + [ˆu(K, n)]nj1

(3.61)

Similarly, by considering the fluxes at the top faces of CE(2)(j, n) and CE(3)(j, n) and imposing the dual-flux

conservation relations similar to Eq. (3.60), one has

[u(L, n)]n

j + [u(L, n)]n

j2 = [ˆu(L, n)]n

j + [ˆu(L, n)]n

j2

(3.62)

[u(M, n)]n

j + [u(M, n)]n

j3 = [ˆu(M, n)]n

j + [ˆu(M, n)]n

j3

(3.63)

We assume the following weighted-averaging relations:

[u(K, n)]n

j = (1 − η)[ˆu(K, n)]n

j + η [ˆu(K, n)]n

j1

(3.64)

[u(L, n)]n

j = (1 − η)[ˆu(L, n)]n

j + η [ˆu(L, n)]n

j2

(3.65)

and

[u(M, n)]n

j = (1 − η)[ˆu(M, n)]n

j + η [ˆu(M, n)]n

j3

(3.66)

NASA/TM—2013-218061 27
where η is a real adjustable parameter. In turn, by substituting Eqs. (3.64)–(3.66) into Eqs. (3.61)–(3.63),

one has

[u(K, n)]n

j1 = η [ˆu(K, n)]nj + (1 − η)[ˆu(K, n)]nj1

(3.67)

[u(L, n)]n

j2 = η [ˆu(L, n)]nj + (1 − η)[ˆu(L, n)]nj2

(3.68)

and

[u(M, n)]n

j3 = η [ˆu(M, n)]nj + (1 − η)[ˆu(M, n)]nj3

(3.69)

respectively.

According to Eqs. (3.64)–(3.66), [u(K, n)]n

j

, [u(L, n)]n

j

and [u(M, n)]n

j

can be determined in terms of

[ˆu(K, n)]n

j

, [ˆu(K, n)]n

j1

, [ˆu(L, n)]n

j

, [ˆu(L, n)]n

j2

, [ˆu(M, n)]n

j

, [ˆu(M, n)]n

j3

, and a given value of η. On the other

hand, because t = t

n at points (j, n), (K, n), (L, n), and (M, n), Eq. (3.55) along with the definitions of

[u(K, n)]n

j

, [u(L, n)]n

j

, and [u(M, n)]n

j

, ⇒

u

n

j + (ux)

n

j

[x(K) − xj ] + (uy)

n

j

[y(K) − yj ] = [u(K, n)]n

j

(3.70)

u

n

j + (ux)

n

j

[x(L) − xj ] + (uy)

n

j

[y(L) − yj ] = [u(L, n)]n

j

(3.71)

and

u

n

j + (ux)

n

j

[x(M) − xj ] + (uy)

n

j

[y(M) − yj ] = [u(M, n)]n

j

(3.72)

Here (i) x(K) and y(K) are the x- and y- coordinates of point K, (ii) x(L) and y(L) are the x- and ycoordinates of point L, and (iii) x(M) and y(M) are the x- and y- coordinates of point M.

Let

det





1 x(K) − xj y(K) − yj

1 x(L) − xj y(L) − yj

1 x(M) − xj y(M) − yj





≡ det





x(L) − x(K) y(L) − y(K)

x(M) − x(K) y(M) − y(K)



 = 0 (3 6 .73)

i.e., points K, L, and M do not lie on a straight line. Then Eqs. (3.70)–(3.72) can be inverted and, as a

result, u

n

j

, (ux)

n

j

, and (uy)

n

j

can be explicitly determined in terms of [u(K, n)]n

j

, [u(L, n)]n

j

, and [u(M, n)]n

j

if the spatial coordinates of points (j, n), K, L, and M are given. As such, given all the spatial geometric

data, the final solution values at any (j, n) can be determined in terms of the given solution values at the

(n − 1/2) time level by using the 2D solution-coupling procedure describe above.

In a CESE 3D scheme using a spatial mesh constructed from tetrahedra, generally a solution point is

associated with four independent mesh variables, four BCEs and one CCE. Also each BCE has two hosts. As

such, construction of the 3D solution-coupling procedure is a straightforward extension of its 2D counterpart.

Finally, it should be emphasized that the use of a solution-coupling procedure generally

will introduce numerical dissipation and, as a result, degrade accuracy. As such it is strongly

recommended that, as long as it serves the purpose of preventing solution decoupling, this

procedure should be applied as sparingly as possible, i.e., it should be used with a small value

of η and/or applied once after a preset number of marching steps.

4. Conclusions and discussions

As a preliminary for the later development, a review of basic CESE ideas and schemes along with a

discussion of several issues related to the application of the CESE method involving unstructured triangular/tetrahedral meshes was presented in Sec. 1. Among the issues discussed are: (A) geometric difficulties

associated with a boundary which has high curvature or sharp corner and is also surrounded by triangular/tetrahedral meshes of extremely high aspect ratio, and (B) possible occurrences of solution decoupling

in CESE simulations involving unstructured triangular/tetrahedral meshes.

To tackle the issue (A), Sec. 2 began by identifying the root cause of the geometry difficulty for the

triangular-mesh case, i.e., the hexagon ABCDEF depicted in Fig. 4 may become highly concave if (i) the

NASA/TM—2013-218061 28
mesh points A, C, E, and G are designated, respectively, as the centroids of the associated triangles, and

(ii) two neighboring triangles with high aspect ratio also have sharply different orientations. A rigorous

mathematical analysis was then presented to validate the observation. According to this analysis, the

centroid and in-center of a triangle are but two special cases with ζ = 0 and ζ = 1, respectively, of the

family of the interior points defined by Eqs. (2.10), (2.29), and (2.32). By using the results presented in

Eqs. (2.39)–(2.41), it becomes clear that, by replacing the centroid with another interior point with ζ ≥ 1,

the hexagon ABCDEF referred to above will become convex or much less concave. This conclusion has been

verified by numerous numerical experiments. In fact, the geometric difficulty can be overcome by simply

designating points A, C, E, and G (see Fig. 4), respectively, as the in-centers (ζ = 1) of the associated

triangles. Moreover, it was shown in Sec. 2 that the tetrahedral-mesh case can be dealt with in a completely

parallel manner.

The new solution-coupling procedures described in Sec. 3 have become an integral part of the CESE

code development. In fact, the nuisance of solution decoupling is avoided completely with the use of these

new procedures.

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, On An Origin of Numerical Diffusion: Violation of Invariance under Space-Time Inversion,

in Proceedings, 23rd Conference on Modeling and simulation, April 30-May 1, 1992, Pittsburgh, PA,

USA, edited by W.G. Volt and M.H. Mickle, Part 5, p. 2727. Also published as NASA TM 105776.

3. 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.

4. S.C. Chang, New Developments in the Method of Space-Time Conservation Element and Solution

Element—Applications to the Euler and Navier-Stokes Equations, NASA TM 106226.

5. 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).

6. 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).

7. 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.

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

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

9. 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.

10. 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. I. The Two-Dimensional Time Marching Schemes, NASA TM—1998-208843, December, 1998.

11. 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. II. Numerical Simulation of Shock Waves and Contact Discontinuities, NASA TM—1998-208844,

December, 1998.

12. 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,

NASA/TM—2013-218061 29

J. Comput. Phys., 156, 89 (1999).
13. 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).

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

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

15. N.S. Liu and K.H. Chen, Flux: An Alternative Flow Solver for the National Combustion Code, AIAA

Paper 99-1079.

16. 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.

17. 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).

18. J. Qin, S.T. Yu, Z.C. Zhang, and M.C. Lai, Direct Calculations of Cavitating Flows by the Space-Time

CESE Method, J. Fuels & Lubricants, SAE Transc., 108(4), 1720 (2000).

19. 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).

20. 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).

21. 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).

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

23. 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).

24. S. Cioc and T.G. Keith, Application of the CESE Method to One-Dimensional Flow in Fluid Film

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

25. S.C. Chang, Courant Number Insensitive CESE Schemes, AIAA Paper 2002-3890 (2002).

26. S. Cioc and T.G. Keith, Application of the CESE Method to Two-Dimensional Flow in Fluid Film

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

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

the CESE Method, STLE Tribology Transactions, 46(1), 128 (2003).

28. S. Cioc, F. Dimofte, and T.G. Keith, Application of the CESE Method to Wave Journal Bearings, STLE

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

29. 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).

30. I.S. Chang, Unsteady Underexpanded Jet Flows, AIAA Paper 2003-3885.

31. S.C. Chang and X.Y. Wang, Multidimensional Courant Number Insensitive CESE Euler Solvers for

Applications Involving Highly Nonuniform Meshes, AIAA Paper 2003-5280.

32. 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).

33. 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).

34. 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 (CESE) Method, Computers and

Chemical Engineering, 28, 1309 (2004)

NASA/TM—2013-218061 30
35. 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).

36. Y.I. Lim, An Optimization Strategy for Nonlinear Simulated Moving Bed Chromatography: Multi-level

Optimization Procedure (MLOP), Korea J. Chem. Eng., 21(4), 836 (2004).

37. 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).

38. 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).

39. 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).

40. S. Jerez, J.V. Romero, and M.D. Rosello, A Semi-Implicit Space-Time CESE Method to Improve Mass

Conservation through Tapered Ducts in Internal Combustion Engines, Math. and Computer Modeling,

40, 941 (2004).

41. 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.

42. 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, International J. Comput. Fluid Dynamics, 19(5),

359 (2005).

43. Y.I. Lim, S.B. Jorgensen, and I.H. Kim, Computer-Aided Model Analysis for Ionic Strength-Dependent

Effective Charge of Protein in Ion-Exchange Chromatography, Bio-Chem. Eng. J., 25(2), 125 (2005).

44. T.I. Tseng and R.J. Yang, Simulation of the Mach Reflection in Supersonic Flows by the CESE Method,

Shock Waves, 14(4), 307 (2005).

45. B. Wang, H. He, and S.-T.J. Yu, Direct Calculation of Wave Implosion for Detonation Initiation, AIAA

J., 43(10), 2157 (2005).

46. S.C. Chang, Explicit von Neumann Stability Conditions for the c-τ Scheme—A Basic Scheme in the

Development of the CESE Courant Number Insensitive Schemes, NASA TM 2005-213627, April 2005.

47. J.C. Yen and D.A. Wagner, Computational Aeroacoustics Using a Simplified Courant Number Insensitive CESE Method, AIAA Paper 2005-2820.

48. I.S. Chang, C.L. Chang, and S.C. Chang, Unsteady Navier-Stokes Rocket Nozzle Flows, AIAA Paper

2005-4353.

49. S.C. Chang, Courant Number and Mach Number Insensitive CESE Euler Solvers, AIAA Paper 2005-

4355

50. S.C. Chang, A. Himansu, C.Y. Loh, X.Y. Wang, and S.T. Yu, Robust and Simple Non-Reflecting Boundary Conditions for the Euler Equations–A New Approach Based on the Space-Time CESE Method, in

Proceedings, NSF-CBMS Regional Research Conference on Mathematical Methods in Nonlinear Wave

Propagation, North Carolina A&T State University, Greensboro, North Carolina, May 15-19, 2002,

edited by D.P. Clemence and G. Tang, p. 155-190, Vol. 379 in Contemporary Mathematics, American

Mathematical Society (2005). Also published as NASA/TM-2003-212495/Rev1.

51. B.S. Venkatachari, G.C. Cheng, and S.C. Chang, Courant Number Insensitive Transient Viscous Flow

Solver Based on CESE Framework, AIAA Paper 2005-00931.

52. T.I. Tseng and R.J. Yang, Numerical Simulation of Vorticity Production in Shock Diffraction, AIAA

J., 44(5), 1040 (2006).

53. M. Zhang, S.-T. Yu, S.C.H. Lin, S.C. Chang, and I. Blankson, Solving the MHD Equations By the

Space-Time Conservation Element and Solution Element Method, J. Comput. Phys., 214, 599 (2006).

54. I.S. Chang, C.L. Chang, and S.C. Chang, 3D Unsteady Navier-Stokes Rocket Nozzle Flows, AIAA Paper

2006-4775.

55. C.L. Chang, Time-accurate, Unstructured-Mesh Navier-Stokes Computations with the Space-time CESE

Method, AIAA Paper 2006-4780.

56. S.C. Chang, On Space-Time Inversion Invariance and Its Relation to Non-Dissipativeness of a CESE

Core Scheme, AIAA Paper 2006-4779. Also published as NASA TM 2006-214385.

NASA/TM—2013-218061 31
57. X. Feng, Y. Zhou, and S.T. Wu, A Novel Numerical Implementation for Solar Wind Modeling by the

Modified Conservation Element/Solution Element Method, Astrophys. J., 655(1), 1110 (2007).

58. C.L. Chang, Three-Dimensional Navier-Stokes Calculations Using the Modified Space-Time CESE

Method, AIAA Paper 2007-5818.

59. S.C. Chang, The a(3) Scheme—A Fourth-Order Neutrally Stable CESE Solver, AIAA Paper 2007-4321.

60. S.C. Chang, The a(4) Scheme—A High-Order Neutrally Stable CESE Solver, AIAA Paper 2007-5820.

61. S.C. Chang, The a(3) Scheme—A Fourth-Order Space-Time Flux-Conserving and Neutrally Atable

CESE Solver (A revised and expanded version of AIAA Paper 2007-4321), NASA TM 2008-215138,

April 2008.

62. C.-L. Chang and M. Choudhari, “Hypersonic Viscous Flow over Large Roughness Elements,” AIAA

Paper 2009-0173, to aapear in Theoretical and Computational Fluid Dynamics Journal, special issue on

Global Instability of Complex Flows.

63. S.C Chang, “A New Approach for Constructing Highly Stable High Order CESE Schemes,” AIAA Paper

2010-543. Also appear as NASA/TM–2010-216777.

64. C.-L. Chang, M. Choudhari and F. Li, “Numerical Computationals of Hypersonic Boundary Layers over

Surface Irregularities,” AIAA Paper 2010-1572, 2010.

65. G.C. Cheng, B. Venkatachari, C.-L. Chang and S.C. Chang, “The Space-Time CESE Method for HighFidelity Flow Simulations,” presented at Parallel CFD 2010 Conference, Kaohsiung, Taiwan, May 2010.

66. J.C. Yen, A.B. Kimbrell and C.H. Connor, “Study of WICS Data Using an Emerging Lower-Order CAA

Method,” AIAA Paper 2010-1742.

67. A.B. Kimbrell, J.C. Yen and C.H. Connor, “CAA Study of WICS Data Using JUSTUS,” presented in

2010 Aircraft-Stores Compatibility Symposium sponsored by International Test and Evaluation Association.

68. Y.Y. Chen, L. Yang and S.T. John Yu, Simulation of Waves in Elastic solids of Cubic Symmetry by the

Conservation Element and Solution Element Method, Wave Motion, in press.

69. L. Yang, Y.Y. Chen and S.T. John Yu, Velocity-Stress Equations for Waves in Solids of Hexagonal

Symmetry Solved by the Space-Time CESE Method, ASME J. of Vibration and Acoustics, in press.

70. Other CESE references are posted on: http://www.grc.nasa.gov/www/microbus.

71. R.A. Horn and C.R. Johnson, Matrix Analysis, Cambridge University Press, 1985.

72. To be presented in a NASA Technical Memorandum.

NASA/TM—2013-218061 32
Figure 1.—A surface element on the boundary

 S(V) of an arbitrary space-time volume V.

x

ds

→

→

→

dr

t

V

S(V)

r = (x,t)

NASA/TM—2013-218061 33
t

∆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

x

j = –3/2

j = –1

j = –1/2

(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).

Figure 2.—The SEs and CEs. 

2(b) SE(j,n).

∆t

∆t

∆x ∆x

2(a)—The space-time mesh.

NASA/TM—2013-218061 34
D

L

B H

M N

A'

A

E

G

G'

E'

P

P'

F

Q

Q'

C

C'

x

x

x

x

x

x

x

y

x

Figure 3.—A spatial domain formed by triangles.

NASA/TM—2013-218061 35
D

L

B

H

M

A

E

G

F

C

Figure 4.—Neighboring triangles with high aspect ratios and sharply different orientations.

NASA/TM—2013-218061 36
Q3 P1

P2

P0

(x2, y2)

(x1, y1)

(x3, y3)

(x0, y0)

P3

Q1

α1

α3

α2

Q2

Figure 5.—The perpendicular projections of paint P0 on the

 three sides (or their extensions) of a triangle ∆P1, P2, P3.

y

z

x

NASA/TM—2013-218061 37
P1

P1

P2

P2

P0

(x2, y2, z2)

(x1, y1, z1)

(x3, y3, z3)

(x0, y0, z0)

P4

P4

(x4, y4, z4)

P3

P3

Figure 6.—Five neighboring tetrahedra in the x-y-z space.

6(a).—A tetrahedron P1 P2 P3 P4 and a point P0.

6(b).—A tetrahedron and four neighboring tetrahedra.

y

z

x

P0

C0

B0

B

C

D0

D

A0

A

NASA/TM—2013-218061 38
C

(j–1/2, n–1/2)

(j–1/2, n)

B

(j–1/4, n)

M–

D

(j, n–1/2)

(j, n)

A

E

(j+1/2, n–1/2)

n–1/2

(j+1/2, n)

F

Figure 7.—Construction of the a–ε–η and c–τ–η schemes.

∆t

2

∆x

4

×

(j+1/4, n)

M+

× n

∆x

4

∆x

4

∆x

4

NASA/TM—2013-218061 39
D

I

B

L K

G

H

M F

J

A'

A

E

G'

E'

C'

C

x

x

x

x

Figure 8.—Construction of the 2D a–ε–η and c–τ–η schemes.

NASA/TM—2013-218061 40

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The 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 Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-09-2013 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Recent Development in the CESE Method for the Solution of the Navier-Stokes Equations 

Using Unstructured Triangular or Tetrahedral Meshes With High Aspect Ratio 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Chang, Sin-Chung; Chang, Chau-Lyan; Yen, Joseph, C. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 599489.02.07.03.04.13.01 

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

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18704 

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

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2013-218061 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 34 and 64 

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

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

In the multidimensional CESE development, triangles and tetrahedra turn out to be the most natural building blocks for 2D and 3D spatial 

meshes. As such the CESE method is compatible with the simplest unstructured meshes and thus can be easily applied to solve problems 

with complex geometries. However, because the method uses space-time staggered stencils, solution decoupling may become a real 

nuisance in applications involving unstructured meshes. In this paper we will describe a simple and general remedy which, according to 

numerical experiments, has removed any possibility of solution decoupling. Moreover, in a real-world viscous flow simulation near a solid 

wall, one often encounters a case where a boundary with high curvature or sharp corner is surrounded by triangular/tetrahedral meshes of 

extremely high aspect ratio (up to 106). For such an extreme case, the spatial projection of a space-time compounded conservation element 

constructed using the original CESE design may become highly concave and thus its centroid (referred to as a spatial solution point) may lie 

far outside of the spatial projection. It could even be embedded beyond a solid wall boundary and causes serious numerical difficulties. In 

this paper we will also present a new procedure for constructing conservation elements and solution elements which effectively overcomes 

the difficulties associated with the original design. Another difficulty issue which was addressed more recently is the well- known fact that 

accuracy of gradient computations involving triangular/tetrahedral grids deteriorates rapidly as the aspect ratio of grid cells increases. The 

root cause of this difficulty was clearly identified and several remedies to overcome it were found through a rigorous mathematical analysis. 

However, because of the length of the current paper and the complexity of mathematics involved, this new work will be presented in another 

paper. 

15. SUBJECT TERMS

The space-time CESE method; Computational fluid dynamics; Numerical analysis 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

48 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18