Xiao-Yen Wang

Glenn Research Center, Cleveland, Ohio

A Summary of the Space-Time Conservation 

Element and Solution Element (CESE) Method 

NASA/TM—2015-218743

June 2015
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 Technical Report Server—Registered 

(NTRS Reg) and NASA Technical Report Server—

Public (NTRS) 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 counter-part 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 co-sponsored 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.

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 757-864-6500

• Telephone the NASA STI Information Desk at

 757-864-9658

• Write to:

NASA STI Program

 Mail Stop 148

 NASA Langley Research Center

 Hampton, VA 23681-2199
Xiao-Yen Wang

Glenn Research Center, Cleveland, Ohio

A Summary of the Space-Time Conservation 

Element and Solution Element (CESE) Method 

NASA/TM—2015-218743

June 2015

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Acknowledgments

The author acknowledges valuable input from Drs. Hung Huynh, Philip Jorgenson, Ken Loh, Wai-Ming To, 

Balaji Venkatachari, Chau-Lyan Chang, and Sin-Chung Chang. The support from Mrs. Mary Jo Davis, 

Dr. James DeBonis, and Dr. James Heidmann is greatly appreciated. This work was funded by the 

Transformative Tools and Technologies project under NASA’s Fundamental Aeronautics Program.

Available from

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. 

This report is a formal draft or working 

paper, intended to solicit comments and 

ideas from a technical peer group.

This report contains preliminary fi ndings, 

subject to revision as analysis proceeds.

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5285 Port Royal Road

Springfi eld, VA 22161

703-605-6000

This report is available in electronic form at http://www.sti.nasa.gov/ and http://ntrs.nasa.gov/
NASA/TM—2015-218743 1 

A Summary of the Space-Time Conservation Element and Solution 

Element (CESE) Method 

Xiao-Yen Wang 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Summary 

The space-time Conservation Element and Solution Element (CESE) method for solving conservation 

laws is examined for its development motivation and design requirements. The characteristics of the 

resulting scheme are discussed. The discretization of the Euler equations is presented to show readers how 

to construct a scheme based on the CESE method. The differences and similarities between the CESE 

method and other traditional methods are discussed. The strengths and weaknesses of the method are also 

addressed. 

1.0 Introduction 

The space-time Conservation Element and Solution Element (CESE) method was introduced by 

S.-C. Chang (Ref. 1) in 1995 as shown in the Figure 1 timeline. In 1999, the two-dimensional (2D) 

triangular and three-dimensional (3D) tetrahedral unstructured Euler solvers were published (Refs. 2 and 

3), which was followed by the extension to the 2D quadrilateral and 3D hexahedral meshes by Zhang 

(Ref. 4) in 2002. The Courant number (i.e., Courant-Friedrichs-Lewy (CFL) number) insensitive scheme 

was presented in 2002 (Ref. 5) and 2003 (Ref. 6) and gave users a better way to minimize the numerical 

dissipation. In the same year, Huynh published a paper to address the similarities and differences between 

CESE and other methods (Ref. 7) and gave an insightful view of the CESE method. Local time stepping 

was introduced in 2005 (Ref. 8) for time-accurate Euler solutions, whereby different timescales are used 

in the computational domain and the time step at each mesh point can vary during the time-marching to 

achieve higher CFL numbers at all locations and all the time. This results in better accuracy (less 

dissipation) and less computational time. The 2D/3D Navier-Stokes solver was published by C.-L. Chang 

in 2006 (Ref. 9) and 2007 (Ref. 10). In 2010, the high-order (fourth-order) CESE scheme was introduced 

by S.-C. Chang (Ref. 11) and was extended to the 2D Euler equations in 2013 (Ref. 12) and 2014 

(Ref. 13). The strategy for handling meshes with a high aspect ratio (a viscous mesh near a solid wall) is 

described in Reference 14. 

The CESE method can be categorized as a cell-centered finite volume method. It aims to solve 

unsteady flow problems. The core ideas that motivated the development of the CESE method are 

(1) enforcing flux conservation over a discretized space-time domain, (2) constructing a nondissipative 

scheme as a baseline to fully control the numerical dissipation, (3) avoiding the Riemann solver and 

dimensional splitting by using a staggered mesh in space and time while maintaining robust handling of 

shock waves and discontinuities, (4) using a compact stencil (use immediate neighboring cells and/or 

points), and (5) using unstructured meshes for ease of dealing complex geometries. 
NASA/TM—2015-218743 2 

The CESE method has been used to solve computational aeroacoustic (CAA) benchmark problems 

involving fan noise and jet noise (Refs. 15 to 20). All problems were successfully solved as blind tests; 

numerical solutions were compared with analytical solutions and showed good agreement, demonstrating 

the ability of the CESE method to solve complicated flow physics (e.g., shocks + waves + eddies). In 

References 21 to 25, the capability of the CESE method for solving engineering problems including highspeed flows with shock waves are demonstrated. 

However, some issues remain, such as 

(1) The CESE solver becomes very dissipative when the CFL number is small (<<1); here, either a 

CFL-insensitive scheme might help (Refs. 5 and 6) or local time stepping for a time-accurate solution 

would be possible way to remedy this (Ref. 8). 

(2) A highly stretched viscous mesh near the wall might degrade the accuracy of the gradients of the 

marching variables. There are a few ways to get better accuracy of the gradients; these can be found in 

References 7 and 14. Recently, a higher order CESE scheme (fourth-order) has been developed (Refs. 11 

to 13), which might enhance accuracy for the viscous flow problems that are dominated by boundaries 

that are solid walls. 

2.0 Discretization 

In the following, the numerical scheme is constructed for the 2D and 3D Euler equations using the 

CESE method. The governing equations are introduced, then the control volume is described in 

comparison with the traditional methods, and the CESE 2D and 3D Euler solvers are presented. 
NASA/TM—2015-218743 3 

2.1 The Governing Equations 

The 3D Euler equations in conservation form are 

 5,4,3,2,1 ,0 





















m

z

g

y

f

x

e

t

u mmmm (1) 

where um are the conservative variables and em, fm, and gm are the x-, y-, and z-component inviscid fluxes 

that are defined as follows: 

  21     43   5 euwuvuuuu t ,,,, (2) 

, 3 4 5 t )(,,, upeeuweuvepueue

2

21  (3) 

vpefvwfpvfvufvf ,, 4 5 t )(,, 2

21 3  (4) 

,,, 5 t )(, wpegpwgwvgwugwg

2

1 2 3 4  (5) 

Here  is density, (u,v,w) are velocity components, p is pressure, and et is the total energy per unit mass. 

Its integral form for a finite control volume V in space (x,y,z) is 

 0d



 





























 V

z

g

y

f

x

e

t

u

V

mmmm (6)

Defining  gfeh mmmm ),,(  in (x,y,z), Equation (6) then becomes 

d    0d





 V Vh

t

u

V m V

m  (7) 

with 

























 

zyx ,, . Using Gauss’s theorem, the volume integration in Equation (7) can be converted 

into surface integration: 

 





)(

d 0d

Vs m V

m shV

t

u   (8)

where s(V) is the boundary of the finite volume V. 

Equation (8) can be discretized using different traditional numerical methods. The volume integration 

of the time derivative allows people to use a different discretization approach for the time derivatives, 

which may be desirable for saving computational time, obtaining higher accuracy for steady-state 

solutions, or maintaining flux conservation in both space and time for transient problems. 

Contrary to the standard approach outlined in Equations (6), (7), and (8), with the CESE method we 

define  ugfeh mmmmm ),,,(  in space and time (x,y,z,t), and Equation (6) becomes 

   

V m Vh 0d  (9) 
NASA/TM—2015-218743 4 

where , ,,































 

tzyx

 and V is a finite control volume in space and time. Next, using Gauss’s 

theorem and converting the volume integration in Equation (9) into surface integration: 

0 d

)(   Vs m sh

  (10) 

The CESE method is about the discretization of Equation (10). Note that the integral form of the 

Euler equation in Equation (10) is the original form of conservation laws. Equations (1) and (7) can be 

derived from Equation (10) only if the flow is continuous. For the time derivative term, CESE does 

surface integration instead of volume integration, which results in no ambiguity for the discretization of 

the time derivative term. Also, the scheme is second-order accurate in both space and time. The flux 

conservation in both space and time is guaranteed. 

In the following section, the discretization of Equation (10) will be described using a staggered 

(CESE method) versus nonstaggered (traditional method) control volume. A 2D unstructured triangular 

mesh is used as an example to explain how to construct a scheme for Equation (10) based on the CESE 

method in comparison with traditional methods. 

2.2 Control Volume (Traditional Method Versus CESE Method) 

Figure 2 shows an unstructured triangular mesh as an example. Consider a triangle, ∆BDF, which has 

three neighboring triangles. Let points A, C, and E be the centroids of those three neighboring triangles. 

With the traditional method, each triangle will be the control volume where the flux is assumed to be 

conserved. The solution point is at the center of each triangle. Then an interface exists between 

neighboring control volumes. Each interface is associated with multiple solution points. Therefore, the 

uniqueness and smoothness of the interface flux has to be addressed. Most of the time, the interface flux 

is computed using the solution point on the upwind side, a mixture of both sides, or an average of both 

sides, which results in the various classes of upwind and central differencing schemes. 

In the CESE method, for the triangle ∆BDF, the control volume is the hexagon ABCDEF formed by 

connecting the three neighboring centroids to the vertices of the triangle ∆BDF. One of the neighboring 

triangles of ∆BDF is ∆DIF, and the corresponding control volume is the hexagon FGDHIJ. Note that 

there is no interface between the control volumes of two neighboring triangles, but the two control 

volumes overlapped on GDEF. The solution point will be the geometric center of the hexagon, not the 

geometric center of the triangle. Each triangle is associated with a hexagon and a solution point, which 

again is the centroid of the hexagon. 

The conservation element (CE) is the same as the traditional control volume but in both space and 

time. The solution element (SE) is the boundary of the CE for each solution point, also in space and time. 

For an example, the CE for point G is the hexagonal cylinder with the base of the hexagon ABCDEF, and 

the SE will be the boundary surfaces of the hexagonal cylinder. 
NASA/TM—2015-218743 5 

Within the SE associated with a solution point (j,n), the numerical approximation of 

 ugfeh mmmmm ),,,(  is denoted as ),,,( *****  ugfeh mmmmm

 with 

)()()()()()()()()(),;,,,( * n m jmtjjmzjjmyjjmxjm  ttuzzuyyuxxuunjtzyxu (11) 

 )()()()()()()()()(),;,,,( * n m jmtjjmzjjmyjjmxjm  ttezzeyyexxeenjtzyxe (12) 

)()()()()()()()()(),;,,,(* n m jmtjjmzjjmyjjmxjm  ttfzzfyyfxxffnjtzyxf (13) 

)()()()()()()()()(),;,,,( * n m jmtjjmzjjmyjjmxjm  ttgzzgyygxxggnjtzyxg (14) 

For the CE at (j,n), Equation (10) becomes 

0 d

)(

*   CEs m sh

  (15) 

Details on solving Equation (15) will be given in Section 2.5, Three-Dimensional Formulation. Here 

the difference between the traditional and CESE methods is summarized in Table I, using ∆BDF as an 

example. 

TABLE I.— COMPARISON BETWEEN TRADITIONAL (UPWIND) AND CONSERVATION ELEMENT 

AND SOLUTION ELEMENT (CESE) METHODS 

[Points A through G are illustrated in Figure 2.] 

 Definition Upwind scheme CESE scheme Comments on CESE 

Control volume (CV) ∆BDF (triangle, 

nonoverlap) 

ABCDEF (hexagon, 

partially overlap) 

Inflated CV, but the same stencil 

Boundary of CV BD, DF, and FB (3) AB, BC, CD, DE, EF, and 

FA (6) 

Doubled number of boundaries 

Number of interfaces between 

CVs of two neighboring 

triangles 

Three Zero (no interface) No need to consider interface 

fluxes 

Number of solution points each 

boundary of CV is associated 

with 

Two One Solution point is on the boundary 

How flux at center of boundary 

is computed 

Taylor series expansion 

from left and right and 

then upwind flux 

Taylor series expansion 

from the solution point 

Uniquely defined 

Example For BD, use solution at 

either C or G or both 

based on upwind 

For BC (and CD), use the 

solution at C 

Sum fluxes through all boundaries 

to update solution at G 
NASA/TM—2015-218743 6 

2.3 Time-Marching for the Primary Unknowns 2/1 )(

n u jm With CESE Method 

Figure 3 shows the control volume associated with the triangle ∆BDF in space and time (x,y,t). Given 

the solution at points A, C, and E at time level n, the solution at point G’ at t = n+1/2 is updated. Each 

solution point at t = n is associated with two surfaces of the control volume along the time direction and 

one surface at a constant time level. The flux is computed at the geometric center of each surface from the 

solution point based on a Taylor series expansion. For the surface on the top (time level n+1/2), the flux is 

simply the solution at the geometric center of the hexagon multiplied by the surface area. Some details are 

given as follows. 

(1) Considering the neighboring point A(j1,n) 

Compute ),,( ****  ufeh mmmm

 at the centroid of surface A’B’BA (A’F’FA ) using the Taylor series 

expansion from the solution point A as shown in Equations (11) to (14), then multiply the surface 

area, sj 1,1 ( sj 2,1 ), and its normal vector in the form of  nnn yx )0,,(  . For the surface of ABGF, 

where time is constant, again compute ),,(

****  ufeh mmmm

 at the centroid of the surface ABGF that has 

the normal vector in the form of n  )1,0,0(  and an area of sj 0,1 , and then we have the flux 

 0,10,1 * 2

1

,1,1 ** 1 ( ) )( jjm

k

kjkjymxm nj   susnfneR



 (16) 

(2) Repeat the same procedure for the second neighboring solution point C(j2,n) to get the flux 2)( n Rj

through the surface of B’C’CB, C’D’DC, and CDGB. 

(3) Repeat the same procedure for the third neighboring point E(j3,n) to get the flux 3)( n Rj through 

E’D’DE, E’F’FE, and EFGD. 
NASA/TM—2015-218743 7 

(4) At t = n+1/2, the surface A’B’C’D’E’F’ whose area is denoted as sj, the flux is simply 

2/1 )(   n us jmj because the solution point is located at the geometric center of the hexagon 

ABCDEFG. 

(5) Summing all fluxes through the CE for solution point G’(j,n+1/2) based on Equation (15) will 

give 

 ()( 321 ) 2/1 n

j nj nj n  us jmj  RRR  (17) 

which gives the explicit time-marching scheme to update 2/1 )(

n u jm . The stability condition of the 

scheme is CFL < 1. 

2.4 Update the Gradients 

A weighted-average limiter based on the gradient amplitude is used to compute umx and umy at G using 

the three neighboring solution points (A, C, and E) (Fig. 4). 

Three sets of gradients, such as  )()( , l

my

l

mx uu , where l = 1,2,3, can be computed using the three planes 

GAC, GCE, and GEA using a finite-difference approach. Define 2)(2)( , )()( l

my

l

lm mx  uuq , and 

   m  mm 3,2,1, m  mm 3,1,2, m  qqwqqwqqw mm 2,1,3, )(,)(,)( , where  = 0,1,2. The weighted-average limiter is 

3,2,1,

)3(

3,

)2(

2,

)1(

1,

3,2,1,

)3(

3,

)2(

2,

)1(

1, , mmm

mymmymmym

my

mmm

mxmmxmmxm mx www

uwuwuw

u

www

uwuwuw

u



  

  (18) 

Equation (18) is the original way of computing the gradients proposed by S.C. Chang (Ref. 2). In 

Reference 6, S.C. Chang describes a modified version of Equation (18) to provide better accuracy of 

gradients of the time-marching variables. 
NASA/TM—2015-218743 8 

2.5 Three-Dimensional Formulation 

Consider a tetrahedron ABCD in Figure 5, in which one of its four neighboring tetrahedrons ABCP is 

also plotted. The tetrahedron shares the face ABC. Points G and H are the centroids of ABCD and ABCP, 

respectively. The polyhedron GABCH is the subconservation element (sub-CE) associated with solution 

point G. 

Similarly, three additional sub-CEs associated with point G can be constructed by considering in turn 

the three tetrahedrons that share one of its other three surfaces with ABCD. The summation of the four 

sub-CEs will be the CE, the control volume where flux is assumed to be conserved. It includes the 

tetrahedron itself (ABCD) and one-fourth of each of its neighboring tetrahedrons. The solution point is 

located at the geometric center of the CE, which is not necessarily at the center of the tetrahedron. 

The control volume in space and time (x,y,z,t) cannot be plotted. The solution at point G’(j,n+1/2) is 

updated using four neighboring solution points at time level n. Each solution point at t = n is associated 

with three volumes of the control volume along the time direction and one volume at constant time. The 

flux is computed at the geometric center of each volume from the solution point based on a Taylor series 

expansion. At the volume on the top (time level n+1/2), the flux is simply the solution at the geometric 

center of the CE multiplied by its volume. The fluxes are then summed over the control volume based on 

Equation (15). Some details are given as follows. 

1. Considering the neighboring point H(j1,n) 

Compute ),,,(

*****  ugfeh mmmmm

 at the centroid of volume HABH’A’B’ (HACH’A’C’ and 

HBCH’B’C’) using a Taylor series expansion from the solution point H, then multiply the volume 

v ,1 kj , where k = 1,2,3 and the normal vector of the surface HAB (HAC, HCB) that is in the form of 

 nnnn zyx )0,,,(  . (H’,A’,B’,C’ is the extrusion of the point H,A,B,C in the time direction, not 

plotted in the figure.) For the volume GABCH, where time is constant (t = n), again compute 
NASA/TM—2015-218743 9 

),,,(

*****  ugfeh mmmmm

 at the centroid of the volume GABCH that has the normal vector in the form of 

n  )1,0,0,0(  and volume of v j 0,1 . We have the flux 

 0,10,1 * 3

1

,1,1 *** 1 ( ) )( jjm

k

kjkjzmymxm nj   vuvngnfneR



 (19) 

2. Repeat the same procedure for the rest of the volumes associated with the other three neighboring 

solution points to get the fluxes n

j nj nj 32 and ,, RRR 4 .

3. At t = n + 1/2, the volume of the CE is v j , and the flux is simply 2/1 )(

n u jm  v j because the 

solution point is located at the geometric center of the CE. 

By summing all fluxes through the CE for solution point G’(j, n+1/2), Equation (15) yields 

()( 4321 ) 2/1 n

j nj nj nj n  uv jmj  RRRR  (20) 

which results in the explicit time marching scheme for 2/1 )(

n u jm . The gradients can be easily extended to 

3D formulation based on Equation (18) and will not be repeated here. 

The location where the flux is computed at the boundary of the CE is summarized in Table II to give 

readers a better idea how it is defined. 

TABLE II.—LOCATION OF THE FLUX IS COMPUTED IN ONE, TWO, AND THREE DIMENSIONS 

Dimension Boundary of conservation element (CE) along the time direction Boundary of CE at t = n

3D Three surfaces, geometry center of each surface, time is at t = t

n + dt/4 One volume, center of the volume, t = tn

2D Two lines, geometry center of each line, time is at t = t

n + dt/4 One surface, center of the surface, t = tn

1D One point, at that point, time is at t = t

n + dt/4 One line, center of the line, t = tn

2.6 CESE Nondissipative Scheme 

The schemes in Equation (17) for 2D and Equation (20) for 3D are used in most computations of 

Euler flow. With the weighted-average limiter, the numerical dissipation is adjusted from the 

nondissipative scheme with added numerical dissipation as needed. The nondissipative scheme is the 

foundation of the CESE scheme and guides the CESE scheme construction under a set of rigorous rules 

(Ref. 1). When solving real problems, how much dissipation is needed depends on the physics of the 

problem. In general, numerical dissipation has to be less than the physical viscosity in order to obtain the 

correct physical solution. Sometimes the physical viscosity is extremely small, and then it requires the 

scheme to be highly accurate yet stable. 

For the 2D case described in Section 2.2, Control Volume (Traditional Method Versus CESE 

Method), the conservation element in space is defined as the hexagon ABCDEF. It can be seen that the 

hexagon is formed by three quadrilaterals, ABGF, BCDG, and DGEF, that can be defined as three subCEs. Assuming Equation (15) is true for the three sub-CEs, there will be three equations for three 

unknowns 2/1 2/1 2/1 )(and ,)(,)(   n

jmy njmx njm uu u for each m. For the 3D case, each conservation element has 

four sub-CEs that yield four equations for four unknowns 2/1 2/1 2/1 2/1 )(and ,)(,)(,)(    n

jmz njmy njmx njm uuu u

for each m. The marching variable 2/1 )(

n u jm and its spatial derivatives are solved using the conservation 

laws for the nondissipative CESE scheme. 
NASA/TM—2015-218743 10 

2.7 CESE Fourth-Order Scheme Compared to Traditional Methods 

The fourth-order CESE Euler solver has been developed recently. Validation and assessment is under 

way (Refs. 12 and 13). The stencil and control volume used in the fourth-order scheme is the same as that 

used in the second-order scheme. The time-marching scheme is the same as that of the second-order 

scheme, but it is used twice. First, the differentialized governing equation for the higher order derivatives 

(uxx and uxxx) is solved; then the original governing equation for u and ux is solved. Higher order Taylor 

series expansion is used to compute the flux on the boundary of control volume. The stability condition is 

CFL < 1, the same as for the second-order scheme. The shock-capturing limiter is the same as that for the 

second-order scheme. Upwinding is not needed since the flux is uniquely computed at the interface. The 

traditional high-order finite-volume methods use a high-order polynomial to represent the solution inside 

each cell. They (1) divide each cell into the subgrid based on the order of accuracy, (2) decide the 

polynomial coefficients using the subgrid solution, and (3) reconstruct the flux on the interface of the 

control volume. Note this has to be carefully handled since the accuracy of the scheme is highly 

dependent on this. The higher the order of accuracy, the more restricted the stability conditions become; 

that is, the smaller the time step size will be. 

3.0 Discussion 

The CESE scheme has some concepts that are similar to those that were previously described in the 

literature. 

Space-time formulation was mentioned in Roe’s paper in 1983 (Ref. 26), but the importance of the 

space-time formulation for the scheme’s accuracy was not discussed. CESE emphasizes space-time flux 

conservation and applies it to the discretization scheme. Given that flux conservation is applied in space 

and time, there is no ambiguity on how to discretize the time derivative terms with the CESE method. 

A staggered mesh was used by Nessyahu and Tadmor (N-T) on a 2D quadrilateral mesh in 1990 (Ref. 

27); but it was not extended to a triangular mesh for 2D or any 3D formulation. The CESE method is 

consistent in its formulation for one, two (triangle), and three (tetrahedron) dimensions. Also, the N-T 

scheme does not have a nondissipative scheme. The CESE nondissipative scheme is the core of the 

scheme development. 

The weighted-average limiter was proposed by van Albada in 1982 (Ref. 28), but it is for a 1D 

formula; the triangular mesh version of the corresponding limit was not developed. The CESE method 

has a coherent formulation of the weighted-average limiter from one, two, and three dimensions, which 

works robustly for capturing discontinuities, such as shock waves and contact surfaces. 

Staggered versus nonstaggered control volumes is the major difference between the CESE and 

traditional methods. Although a staggered control volume avoids the interface so there is no need of 

upwinding for the flux at the interface, its effect on the overall performance of the scheme compared with 

that of a nonstaggered control volume needs to be further evaluated using real engineering problems; that 

is, 3D unsteady viscous turbulent flows. The performance of the scheme includes the cost (computation 

time), accuracy, and robustness. These efforts are underway. 

4.0 Concluding Remarks 

The Conservation Element and Solution Element (CESE) method for solving conservation laws has 

been described and compared with other traditional methods. The numerical schemes for both two- and 

three-dimensional Euler equations were illustrated based on the CESE method. The major differences and 

similarities between the CESE method and other traditional methods have been highlighted and addressed 

from different perspectives. 
NASA/TM—2015-218743 11 

References 

1. Chang, Sin-Chung: The Method of Space-Time Conservation Element and Solution Element—A 

New Approach for Solving the Navier-Stokes and Euler Equations. J. Comput. Phys., vol. 119, no. 2, 

1995, pp. 295–324. 

2. Chang, S.-C.; and Wang, X.-Y.: A 2D Non-Splitting Unstructured Triangular Mesh Euler Solver 

Based on the Space-Time Conservation Element and Solution Element Method. Computational Fluid 

Dynamics Journal, vol. 8, no. 2, 1999, p. 309. 

3. Wang, Xiao-Yen; and Chang, Sin-Chung: A 3D Non-Splitting Structured/Unstructured Euler Solver 

Based on the Space-Time Conservation Element and Solution Element Method. AIAA–99–3278, 

1999. 

4. Zhang, Zeng-Chan; Yu, S.T. John; and Chang, Sin-Chung: 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., vol. 175, 2002, pp. 168–199. 

5. Chang, Sin-Chung: Courant Number Insensitive CE/SE Schemes. AIAA 2002–3890, 2002. 

6. Chang, S.C.; and Wang, X.Y.: Multi-dimensional Courant Number Insensitive CE/SE Euler Solvers 

for Applications Involving Highly Nonuniform Meshes. AIAA 2003–5285, 2003. 

7. Huynh, H.T.: Analysis and Improvement of Upwind and Centered Schemes on Quadrilateral and 

Triangular Meshes. AIAA 2003–3541, 2003. 

8. Chang, Sin-Chung, et al.: Local Time-Stepping Procedures for the Space-Time Conservation Element 

and Solution Element Method. International J. Comput. Fluid Dynamics, vol. 19, no. 5, 2005, 

pp. 359–380. 

9. Chang, Chau-Lyan: Time-Accurate Unstructured-Mesh Navier-Stokes Computations With the SpaceTime CESE Method. AIAA 2006–4780, 2006. 

10. Chang, Chau-Lyan: Three-Dimensional Navier-Stokes Calculations Using the Modified Space-Time 

CESE Method. AIAA 2007–5818, 2007. 

11. Chang, Sin-Chung: A New Approach for Constructing Highly Stable High Order CESE Schemes. 

AIAA–2010–543, 2010. 

12. Chang, Chau-Lyan; Venkatachari, Balaji; and Cheng, Gary: Time-Accurate Local Time Stepping and 

High-Order Space-Time CESE Methods for Multi-Dimensional Flows With Unstructured Meshes. 

AIAA 2013–3069, 2013. 

13. Bilyeu, David, et al.: A Two-Dimensional Fourth-Order Unstructured-Meshed Euler Solver Based on 

the CESE Method. J. Comput. Phys., vol. 257, pt. A, 2014, pp. 981–999. 

14. Chang, Sin-Chung; Chang, Chau-Lyan; and Yen, Joseph C.: Recent Developments in the CESE 

Method for the Solution of the Navier-Stokes Equations Using Unstructured Triangular or 

Tetrahedral Meshes With High Aspect Ratio. AIAA 2013–3068 (NASA/TM—2013-218061), 2013. 

15. Wang, Xiao-Yen; Chow, Chuen-Yen; and Chang, Sin-Chung: Numerical Simulations of Gust 

Generated Aeroacoustics in a Cascade Using the Space-Time Conservation Element and Solution 

Element Method. AIAA 98–0178, 1998. 

16. Wang, Xiao-Yen; Chang, Sin-Chung; and Jorgenson, Philip C.E.: Prediction of Sound Waves 

Propagating Through a Nozzle Without/With a Shock Wave Using the Space-Time CE/SE Method. 

AIAA 2000–0222 (NASA/TM—2000-209937), 2000. 

17. Wang, X.Y., et al.: Gust Acoustic Response of a Swept Rectilinear Cascade Using the Space-Time 

CE/SE Method. NASA/TM—2001-210815 (FEDSM−2001–18134), 2001. 

18. Wang, X.Y., et al.: Gust Acoustics Response of a Single Airfoil Using the Space-Time CE/SE 

Method. AIAA–2002–0801 (NASA/TM—2003-211513), 2002. 

19. Loh, Ching Y.; Hultgren, Lennart S.; and Chang, Sin-Chung: Wave Computation in Compressible 

Flow Using the Space-Time Conservation Element and Solution Element Method. AIAA J., vol. 39, 

no. 5, 2001, p. 794. 

20. Loh, Ching-Y.; and Zaman, K.B.M.Q.: Numerical Investigation of Transonic Resonance With a 

Convergent-Divergent Nozzle. AIAA J., vol. 40, no. 12, 2002, p. 2393. 
NASA/TM—2015-218743 12 

21. Chang, Chau-Lyan; Choudhari, Meelan M.; and Li, Fei: Numerical Computations of Hypersonic 

Boundary-Layer Over Surface Irregularities. AIAA 2010–1572, 2010. 

22. Chang, Chau-L., et al.: Effects of Cavities and Protuberances on Transition Over Hypersonic 

Vehicles. AIAA 2011–3245, 2011. 

23. Wang, Bao; He, Hao; and Yu, S.-T. John: Direct Calculation of Wave Implosion for Detonation 

Initiation. AIAA J., vol. 43, no. 10, 2005, pp. 2157–2169. 

24. Tseng, T.I.; and Yang, R.J.: Simulation of the Mach Reflection in Supersonic Flows by the CE/SE 

Method. Shock Waves, vol. 14, no. 4, 2005, pp. 307−311. 

25. Liu, Kaixin, et al.: A Review on the CE/SE Method. Adv. in Mech., vol. 41, no. 4, 2011, 

pp. 447−461. 

26. Roe, P.L.: An Introduction to Numerical Methods Suitable for the Euler Equations. VKI for Fluid 

Dynamics Lecture Series 1983−01, von Karman Institute for Fluid Dynamics, Belgium, 1983. 

27. Nessyahu, Haim; and Tadmor, Eitan: Non-oscillatory Central Differencing for Hypersonic 

Conservation Laws. J. Comput. Phys., vol. 87, no. 2, 1990, pp. 408−463. 

28. van Albada, G.D.; van Leer, B.; and Roberts, W.W., Jr.: A Comparative Study of Computational 

Methods in Cosmic Gas Dynamics. Astron. Astrophys., vol. 108, no. 1, 1982, pp. 76−84.