b 

NASA Technical Memorandum 4020b9' 

,9 11 I! I 

ICOMP-89-13 

I 

A Genuinely 

Upwind Cellfor the Euler Equ 

_. 

[HA SA-TU-1 0 202 9) A GEHUINELY 

MULTI-DII!lENSIO#BL UPWIND CELL-YBRTEX SCHEME 

POR THE EULEB EQUATIONS 

22 P CSCL 12A Dnclas 

N 89- 2 48 72 

[Hichigan Uaiv,) 

G3/64 0217234 

Kenneth G. Powell 'and Bram 

The University of Michigan 

Ann Arbor, Michi& 

and Institute for Cdrnputational 

Lewis Research Center 

Cleveland, Ohio 

May 1989 

I 

B 
A Genuinely Multi-Dimensional Upwind Cell-Vertex 

Scheme for the Euler Equations 

Kenneth G. Powell’ and Bram van Leer’ 

Department of Aerospace Engineering 

The University of Michigan, Ann Arbor, Michigan 48109 

and Institute for Computational Mechanics in Propulsion 

Lewis Research Center, Cleveland. Ohio 44135 I 

Abstract 

A scheme for solving the twedimensional Euler equations is developed. It is based on a new scheme for 

the two-dimensional linear convection equation, and the 

Euler-equation decomposition developed by Hirsch et 

al. [l]. The scheme is genuinely two-dimensional. At 

each iteration, the data are locally decomposed into 

four variables, allowing convection in appropriate directions. This is done via a cell-vertex scheme with a 

downwind-weighted distribution step. The scheme is 

conservative, and third-order accurate in space. The 

derivation and stability analysis of the scheme for the 

convection equation, and the derivation of the extension 

to the Euler equations are given. Preconditioning techniques based on local values of the convection speeds 

are discussed. The scheme for the Euler equations is 

applied to two channel-flow problems. It is shown to 

converge rapidly to a solution t,hat agrees well with that 

of a third-order upwind solver. 

In t rod u c t ion 

Much of the understanding of modern upwind schemes 

for the Euler equations has come from designing algorithms for the onedimensional linear convection equat ion 

au au 

- + c- = 0. 

at az 

As a consequence of this, problems in two or three dimensions are typically solved in a direction-split manner, with the upwinding directions normal to the faces 

of the computational cell. This leads to schemes that 

are strongly coupled to the grid on which they are applied. Discontinuities that lie along grid lines are represented properly when treated in this manner, but ones 

that are oblique to the grid are interpreted incorrectly 

by the built-in “Riemann solver” [2]. This suggests the 

need for designing an upwind-differencing scheme for 

the Euler equations that is truly multi-dimensional, and 

therefore less strongly coupled to the grid. The design 

‘Work funded under Space Act Agreement C99066G. 

of an algorithm of this type should be motivated by the 

two-dimensional linear convection equation 

and by an understanding of the wave-like character of 

the two-dimensional Euler equations. 

Characteristic information has been used in the past 

to formulate schemes. hloretti’s X-scheme [3] and 

the QAZlD algorithm of Verhoff and O’Neil [4] are 

two examples of non-conservative, characteristic- based 

schemes that use grid-decoupled stencils. Conservative 

schemes that are decoupled from the grid are more rare, 

however. Davis [5] has formulated an upwind scheme 

in which the Riemann problem is not solved normal to 

cell faces, but normal to shock waves. Levy et ai. [6] 

have extended this work, including other possible upwinding directions. Hirsch et al. [I] have developed a 

method of decomposing the Euler equations into a set 

of convection equations. They have formulated a firstorder scheme based on this decomposition. Roe [2] has 

developed a different decomposition method, based on 

locally decomposing the data into waves. He has formulated a first-order scheme that makes use of his decomposition. All of these conservative algorithms are 

extremely nonlinear. Differences are not taken in gridcontravariant directions, but in directions determined 

by local values of the flow variables. In general, the directions are actually based on derivatives of flow quantities. For this reason, these schemes are inherently less 

robust than schemes that use the grid-contravariant direc t ions. 

Most of the upwind schemes used to date are cellcentered schemes. While cell-vertex schemes have advantages in terms of accuracy [7], the ones that have 

been developed for the Euler equations thus far are 

based on central differencing [8,9] or on the LaxU’endroff scheme [10,11]. In the central-differenciilg 

version of a cell-vertex scheme, the residual for the cell 

is distributed equally to the four nodes of the cell. In 

the Lax-Wendroff version, this distribution is altered 

by the higher-order terms, so that the nodes receive 

unequal portions of the residual. This can be generalized so that the nodes receive some weighted fraction 

1 
of the residual, where the weight is determined from 

the stability analysis of the scheme. For a convection 

problem, the weights should be such that the residual 

is "pushed" downwind. For the Euler equations, there 

is the added difficulty of determining what variables 

should be convected, and in what directions. 

For the design of a genuinely multi-dimensional upwind cell-vertex scheme, then, the following compoand are approximated here by the third-order accurate 

one-parameter formulas 

Ui,j + Ui+l,j 

2 ui+1 . = a 13 

6: -- ((1 + 6,) ui+l,j + (1 - 0,) ui,J } (3) 24 

Uij + Ui,j+l 

2 nents are necessary: Ui,j++ = 

1. a cell-vertex scheme with a downwind-biased distribution for a scalar convection equation in which 

the grid components of the convection speed are 

known; 

62 -6 ((1 + @y) ui,j+l + (1 - ey) Uij} t (4) 

where 62 and by' are the centered second-difference operators 

2. a method of locally decomposing the Euler equa- 6zuipj = ui+l,J - 2uilj + ~i-l,j tions into a set of convection equations; 

62uij = ui+l,j - 2uij + ui-l,j . 3. an extension of the scalar scheme to a system, such 

that mm, momentum and energy are conserved. Using these formulas for the cell-face averages, the 

Fourier footprint of the flux integration for a cell is 

These components are described below for the case of 

two-dimensions. F(AtRes)=-L [,,sin$ (13~0s~ P 

6 2 

Scheme for the Convection Equation 

The heart of the new scheme for the Euler equations is 

a cell-vertex scheme for a two-dimensional convection 

equation 

aU aU au 

at ax - = -cz- - 

This scheme is analyzed below. On a uniform Cartesian 

grid, the residual for the convection equation is given 

by 

I where the semi-integer index denotes a average over a 

cell-face, i.e. 

These cell-face averages, to fourth order, may be written as 

3Qz 

3sin $)I , 

where the p's are the Fourier variables and the u's are 

the Courant numbers 

c, At 

u, = - 

Ax 

cy At 

vy = -. 

AY 

To update the nodes, the cell-centered residual, given 

in Equation 2, multiplied by At, will be sent to the 

nodes (i.j), (i + lj), (i + 1,j + 1) and (i,j + 1) with 

weights wIwr w,~, wne, and w,, respectively (see Figure 1). The Fourier footprint of this distribution step 

is given by 

+ (wnW - wse) sin - -"I . (5) 2 

2 
If a simple forward-Euler time-stepping scheme is 

used, the net amplification factor for the entire scheme 

is 

1 

G(V,,~~,~,,~~) = 1 + 3(At&s)F(Dist) (6) 

The appropriate values for the 0’s and thew’s remain 

to be determined from the stability analysis. The w’s 

correspond to convection directions, and should therefore be determined by enforcing st,ability for the long 

waves (/3z,,b’y + 0). The 0’s control the high-order difference terms of the scheme, and should therefore be 

determined by enforcing stability for the short waves 

Taking the limit of Equation 6 as &,,,By -.+ 0 yields 

(&, Py - .). 

the constraint 

k2 (bPz + VyPy) = 

where k is some real constant. An added constraint 

(conservation) is that the entire residual must be distributed, 

 ne 4- wsw + wnw +use = 1 

Also, by symmetry, if the v, and vy are such that convection is directly towards one node, all of the residual 

is distributed to that node, i.e. 

w,, = 1 if v, = vy > 0 ; 

wIw = 1 if v, = vy < 0; 

wnw = 1 if -v, = vy > 0; 

w,, = 1 if -u, = vy < 0 

Combining these conditions gives the distribution coefficients 

wne = max (0, ) (7) (vr + VY) 

I4 + vyl + Ivz - vyl 

These formulas state that the residual is sent only to 

the nodes that define the downwind face, and is distributed in a weighted manner between the two nodes 

on that face. For a plane wave moving in one of the 

coordinate directions, the two downwind weights are 

equal, and the scheme reduces to the standard onedimensional upwind scheme. 

i,j+l i + l,j + 1 

.. 1 

Figure 1: Cell-Vertex Distribution 

A short-wave analysis shows that a necessary condition for stability is 

(ivz+vyi- IV, -vy~)(vzey+vy~,) > 0. 

r=3 e 

0, v, 

For this to hold for all values of v, and vy , the constraint 

must be met, so that the 0’s must be given by 

e, = 

ey = atv, 

For steady solutions that are independent of At, the 6’s 

must be independent of At. This leads to the choice 

(12) avy ey = max(lvzI1 bYl) 

where a is a positive parameter of order one. It is 

interesting to note that this says that the 0’s must be 

downwinded, i.e., in the cell-face average calculation of 

Equations 3 and 4, the 0’s must be chosen so as to 

give more influence to the second-difference about the 

downwind node of the face. 

The only parameter in the scheme that remains to 

be determined is the value of a. Figures 3-5 show the 

effect of different values of cr on the amplification factor C of the scheme. The maximum amplitude of t he 

amplification factor over the high-frequency region (see 

Figure 2) 

3 
ur = 0.90 uy = 0.30 

Damping characteristics 

0 75- 

Cmsr 

0 50- 

- 

0 25- 

- 

0 00 

High Requcncy Domain 

II1111111 

orlr 

Figure 2: High-Frequency Region Used to Determine a 

ut = 1.00 Yy = 0.00 

Figure 3: Effect of Q on Stability - 0" \Vave 

a 

Figure 4: Effect of cr on Stability - 20' Wave 

u, = 0.ao uy = 0.80 

DamDinR characteristics 

U 

Figure 5: Effect of Q on Stability -- 45' \Vavr 

4 
New Cell-Vertex Scheme 

St ability Boundary 

us = 1.00 and uY = 0.00 

Figure 6: Stability Boundary for Scheme 

c 

is plotted as a function of a for waves traveling at O", 

20" and 45". Based on these results, the value Q = 1.0 

was chosen for the scheme. The stability boundary for 

the scheme, with a = 1.0, is shown in Figure 6. 

The locus of the scheme (i.e. the Fourier footprint 

of F(AtRes)F(Dist)) is shown for the Oo, 20" and 45" 

waves in Figures 7-9. The plots are generated by varying p, continuously and p,, discretely, which leads to a 

mesh of points within the continuous footprint of the 

locus. The actual locus of the scheme consists of the 

lines in the plot, and all the space between the lines. 

The circular stablity boundary of forward-Euler timestepping is circumscribed about the loci for reference. 

The loci are very different from those of first-order upwind or central-difference schemes. It is the wave that 

is convected at 45' that is damped the best, while waves 

at 0" (or 90") are not damped well. This can be seen 

clearly in the contours of the amplification factor for 

the scheme, shown in Figures 10-12. 

Some numerical results for a convected Gaussian on 

a 32 x 32 grid are given in Figures 13-15. In each case 

the Gaussian propagates across the grid virtually undistorted. The onset of a zebra instability can be seen in 

the 0" case, as predicted in the stability analysis. The 

amplitude of oscillations in this case is very small (on 

the order of The convergence history for each of 

the cases is shown in Figures 16-18. The Gaussian convects at almost one cell per iteration. so that the slope 

of the residual curve changes drastically after approximately forty iterations. The 45' case, which has the 

best high-frequency damping, converges very quickly, 

while the 0' case converges very slowly. Table 1 shows 

the results of a grid-refinement study confirming the 

third-order accuracy of the scheme. 

Wave number locua 

-2.00 -1.60 -1.20 -0.80 -0.40 0.00 

Wz) 

Figure 7: Fourier Footprint of Scheme - 0" wave 

us = 0.90 and vY = 0.30 

-2.00 -1.60 -1.20 -0.80 -0 40 0 00 

WZ) 

Figure 8: Fourier Footprint of Scheme - 20" wave 

5 
= 0.80 Md vu = o.ao 

Wave number locus 

1.00 1 

-1.00: I , I 

-2.00 -1.60 -1.20 -0.80 -0.40 0.00 

Wz) 

Figure 9: Fourier Footprint of Scheme - 45" wave 

ur = 1.00 and uI = 0.00 

Amplification factor 

1.00 , 1 

1 

0.60 

0.20 

-0 20 

-0.60 

1 0 0000 

2 0 1000 

3 0 2000 

4 0 3000 

5 0 4000 

6 0 5000 

7 0 6000 

8 0 7ooil 

9 0 8000 

10 0 9000 

0 

Figure 11: Amplification Factor of Scheme - 20' wave 

-I 00 -0.60 -0.20 0 20 0 60 100 

PI/* 

us = 0.90 and uu = 0.30 

AmDliflcation factor 

u, = 0.M) and uY = 0.80 

Amplification factor 

1 0.0000 

2 0.1000 

3 0.2000 

4 0.3000 

5 0.4000 

6 0.5000 

7 0.6000 

9 0.8000 

10 0 9000 

a 0.7000 

Figure 10: Amplification Factor of Scheme - 0" wave Figure 12: Amplification Fact,or of Scheme - 45" wave 

6 
1.00 

1 I 

- 1 J' ' 

>-A 

0.80-, / - 

C / - 

- - 1- 

1 - 

1 

0.60- 

Y 

- 

0.20 - 1 0.0000 

2 0.1000 

3 0.2000 

4 0.3000 

5 0.4000 

6 0.5000 

7 0.6000 

8 0.7000 

9 0.8000 

10 0.9000 

10 0.00 IIIIIIIII 

0.00 0.20 0.40 0.60 0 80 100 z 

Figure 13: Gaussian Convected at 0' 

Y= = 0.80 vI = 0.80 C = 0.25 

Convected Gaussian 

10.0000 

2 0.1000 

3 0.2000 

4 0.3000 

5 0.4000 

6 0.5000 

7 0.6000 

8 0.7000 

9 0 8000 

10 0.9000 

0.40 

0.20 

0.00 

0.00 0.20 0.40 0.60 0.80 1.00 r 

Figure 15: Gaussian Convected at 45' 

vt = 0.90 vI = 0.30 C = 0.25 vI = 1.00 uV = 0.00 .$ = 0.25 

I iteration 

Figure 14: Gaussian Convected at 20' Figure 16: Convergence History for 0' Case 

10 

7 
Y, = 0.90 vY = 0.30 

Convernence History 

= 0.?5 

Iteration 

Figure 17: Convergence History for 20' Case 

11 Grid 

Table 1: Order of Accuracy Study - Convection of 

Unit-Ampilitude Gaussian, ( = 0.25, at 4.5' 

Scheme for the Euler Equations 

The Euler equations are 

- dU + - 8F + dG -- - 0, at az ay 

where U is the state vector of conserved variables 

and F and G are the flux vectors 

Y, = 0.80 Y# = 0.80 C = 0.25 

Convergence History 

1 -1.00 

-2.00 

Residual 

Figure 18: Convergence History for 45' Case 

To solve the Euler equations with a scheme analagous 

to the one above, the system must be decomposed into 

a set of two-dimensional convection equations, with or 

without source terms. Once the equations have been 

decomposed, each component can be treated with the 

convection scheme described above. The distribution 

step carries over in a very straightforward manner; the 

flux calculation (particularly the higher-order terms) 

must be treated carefully to ensure that the formulation 

for the system is consistent with the formulation for 

the scalar equation, and that, the resulting scheme is 

conservative. 

Decomposition of the Euler Equations 

Roe [2,12] has formulated a decomposition of the twodimensional Euler equations, based on the eigenvectors 

of the matrix 

dF dG . A = ---coso+ -smB dU dU 

The eigenvectors of A represent a shear wave, a contact 

discontinuity and two acoustic waves. Roe makes use 

of these eigenvectors, decomposing the flow into, for 

example, four acoustic waves, one shear wave and an 

entropy wave. He uses local values of the flow gradients 

to compute the strength and angle of inclination of each 

of the waves in his model. 

Hirsch et al. [I] have formulated a different decornp 

sition, which converts the Euler equations to the form 

8W' + Dx- +D -- - s, ax ay 

dW* 8W' 

at 

- 

where W' is a vector of convected quantities (entropy, a 

component of velocity, and two acoustic-like variables). 

8 
Dx and Dy are diagonal matrices of convection speeds, 

and S is a source term. Similarly to Roe's decomposition, the decomposition of Hirsch is based on the matrix 

dF aG 

-6, + - 

- au aU - 

where IC= and ny are the components of a unit propagation vector. Diagonalization of A gives 

A = PD(6)P-l 

where 

1U.K 0 0 O1 

0 0 u.n-c J lo 

and P and P-' are the left- and right-eigenvector matrices. The choice of the unit vector K is based on 

local flow gradients. Hirsh chooses K so as to minimize 

the components of the source term S. To minimize 

all components, one needs two different K vectors; 

for the velocity-component convection and K(~) for the 

acoustic-like convection. With two K vectors, 

and 

1. O 0 O1 

0 

0 

0 0 U+K?)C 

Dx= I:: 0 

r' O 0 O1 

v - npc 

0 

0 0 V+KY (2) c 

Dy = I: : 0 

and 

0 

These equations are general, holding for any choices of 

IC(') and Hirsch shows that, in order to minimze 

th source terms, one needs a IC(') that is aligned locally 

with the pressure gradient, and a that is related to 

the strain-rate tensor. That is, K(') is given by 

and 

the following way: if 

is computed from the velocity derivatives in 

auax 2 (E+$) -4avay 50 

then the propagation angle 

is calculated; otherwise, the possible propagation angles 

are given by 

2 (%+%)*d(%+%) -4ce 

tme = . (14) 2% 

The value of K(~) is then 

K(~) = cos ei + sin ej . 

The proper branch for Equation 14 is the one that maximizes the inner product IC(~) . da). This inner product 

appears in the denominator of entries of the transformation matrix P' (described below); the two vectors 

therefore must not be perpendicular. This is ensured 

by taking K(~) := n(') if the nominal value of the inner 

product is less than 1/10. 

Hirsch's decomposition was chosen for this study because the matrix P' is square (4 x 4)' as opposed to 

Roe's decomposition, which yields a 6 x 4 matrix. 

Extension of Convection Scheme 

Just as in the scheme for the convection equation, the 

scheme for the Euler equations is made up of two primary steps: 

1. a residual calculation based on a flux integral; 

9 
2. a residual distribution. 

Each of these steps is somewhat more complicated for 

the system, however. 

For the Euler equations, the residual for a cell is given 

by 

(15) 

%Si++,j++ =-;;ii,[Fdy-Gdz] 1 t 

where A is the cell area and the integral is taken along 

the cell's boundary 8A. The cells are now quadrilaterals, with facea lying along curvilinear coordinate lines 

< = & and = 9,. The boundary integral of Equation 15 is composed of contributions of fluxes normal 

to cell faces. To extend the approximation of Equations 3 and 4 90 that they apply to the above residual, 

it is necessary to convert these to flux approximations. 

Equation 4, for instance, when multiplied by c,Ay, becomes an expression for the total flux acrms the cellface centered at (i,j + 1/2): 

With regard to the Euler equations, this translates directly into 

1 +z (0, )i,j++ [(Fi,j+l - Fi,j)A,y 

- (Gi,j+l - Gij) As"] } 7 

in which the following notation is used: 

A7x = ~ij+l - ~;,j 

A7Y = !Ji,j+l - Yij 

AV = Vj+l-qj 

The quantity $' denotes the flux normal to a cell-face, 

scaled such that 

Fdq = Fdy - Gdz . 

The matrix 0, replaces the scalar Oy; it acts as a scalar 

in each of the convection equations generated by the 

transformation matrix P' 

0, = P'diag {e?)} Po-' , ( 16) 

with 

The contravariant Courant numbers vik) and vp' and 

the transformation matrix P' are defined further below. 

Note that P' and fly) are defined per cell, and must be 

averaged over neighboring cells to yied a cell-face value. 

The analog of Equation 3 is 

with 

and 

Gd< = Fdy - Gdx . 

The matrix 0~ is given by 

OE = P'diag {e?'} P*-' , (17) 

with 

The contravariant Courant numbers vlk) and v,(jk) are 

related to the wave speeds normal to the cell faces. 

Thus we have, for instance, 

with 

1 

(Aq~)i++,j++ = -(xij+l 2 - zij + zi+lj+l- Zi+l,j) 

(&~);++,j++ = 5 (Yi,j+l - Yi,j + Yi+l,j+l - Y~+I,]) 

The Cartesian wave speeds c, and cy are the diagonal 

elements of the matrices Dx and Dy introduced previously, and are evaluated at (i + 1/2. y + 1/2) by using 

cell-averaged state quantities and the cell-centered K. 

values. Note that the factor At in Equation 18 drops 

1 

10 
i 

out in the definition of By). Analagous to Equation 18 where L varies from one to four. In the above, 

is denotes the kth column of P’ corresponding to the kth 

wave, and btU(k) is the portion of time change in the 

= A (c=ACy - cyAez). A . A , (19) state vector caused by the kth wave. Each 6tU(k) may 

be divided between the two vertices of the downwind 

cell-face, according to the weights of Equations 7-10, 

with v, and vy replaced by vlk) and vf’, e.g. 

At 

lt aJt a 

with 

1 -- (Acz)j++j++ - (zi+l,j - xij + zi+l,j+l - Zi,+lj) 

1 

(ACY)~++~+* = 5 (gi+l,j - gi,j + gi+l,j+l - Yi,j+l) . 

The transformation matrix P’ is given, in columns, 

, 

- 

1 

U 

V 

2- 

Oi 

The distribution step requires, in each cell, projection of the residual onto the columns of the matrix P’, 

giving weights r(’)), and multiplication of each of the 

resulting vectors by an appropriate time-step: 

As indicated in Equation 21, it is not necessary to 

take the same value of At for each cell, or even for each 

wave. Spatial conservation, and, therefore, the ability 

to find steady weak solutions, is guaranteed by formulating the discrete residual on the basis of Equation 15. 

Using different time steps in different cells is a wellknown technique called “local time-stepping;” using different time steps for different waves is new, and will be 

called “characteristic time-stepping.” Mathematically 

speaking, the use of a non-constant At is equivalent 

to preconditioning the equations. Local time-stepping 

takes away the stiffness due to spatial variations, and 

may be called “spatial preconditioning;” characteristic 

time-stepping removes the stiffness due to the differences among the local wave speeds, and may be called 

“wave-preconditioning .” 

In local time-stepping, one chooaes, with some safety 

margin, the largest single time-step value that satisfies 

the stability criterion (see Figure 6) for the local values 

of all pairs v:’), vik’. In characteristic time-stepping, 

the time-step for each pair is maximized separately. 

The validity of this practice hinges on the assumption 

that, in the steady state, each residual (Equation 15) 

vanishes separately; for a cell-vertex scheme this, however, is not generally true. All one can assume is that 

the sum of all residual components sent to a particular 

vertex vanishes in the steady state. To prevent an imbalance among contributions 61 U(k) for a particular Ll 

arriving in a vertex from different cells, both local and 

characteristic time-step values need to be assigned to 

vertices (i, j) rather than cell-centers (i + 1/2, j + 1/2). 

The use of characteristic time steps requires special 

provisions near sonic lines, steady shocks and stagnation points, i.e. in regions where one of the convection speeds vanishes. The linear stability criterion then 

allows of arbitrarily large values of At; in practice, 

however, its value must be constrained by a solutiondependent upper limit. How to do this robustly in 

the multi-dimensional case is not yet known: for onedimensional flow, some progress has been reported [13]. 

Boundary conditions are imposed at the walls by enforcing a tangency condition at the vertices on the walls 

and zeroing the mass-flux for the faces on the walls. 

Boundary conditions at inflow and outflow are imposed 

by a non-reflecting condition described by Lindquist 

and Giles [14]. 

11 
Results 

The scheme described above was used to compute 

steady flows in a twedimensional channel, with cosineshaped walls yielding a 10% constriction at the throat. 

Two different inflow Mach numbers were taken; 

1. M, = 1.75; 

2. M, = 3.50. 

Both cases were run on the 64 x 32 grid shown in Figure 19. It is worth mentioning here that the full grid 

was used in the calculations: although the final steady 

state is symmetric about the channel axis, the transient 

states are not, due to the asymmetry of the This 

asymmetry of the Hirsch decomposition is suspected to 

have a negative effect on the convergence to a st,eady 

state; this remains to be investigated. In comparison, 

the decomposition of Roe preserves flow symmetry. 

The results for the first case are shown in Figures 20- 

30. Figure 20 shows the Mach number contours of the 

steady flow. The compression waves caused by the cosine bump, the coalescence into shock. waves, and the 

reflection of the shocks are clearly seen. There are 

some oscillations at the shocks, due to the fact that 

the scheme is third-order accurate everywhere. For 

comparison, results from a grid-biased cell-centered up 

wind scheme [15] are shown in Figures 21 and 22. This 

scheme is linearly third-order accurate but formally 

third-order accurate only for one-dimensional flow; the 

approximate Riemann solver used for upwinding is Osher’s. Figure 21 shows the results produced without 

limiting of higher-order terms, allowing a fair comparison with Figure 20. The results are very similar, the 

cell-centered scheme being a little less oscillatory, but 

yielding less defined reflected shocks. The results of 

Figure 22 were rendered monotone by the use of Van 

Albada’s limiter [16], leading to a clear loss of resolution 

for all waves, especially the reflected shocks. The remaining plots show the details of the shock-intersection 

region. Figure 23 shows the grid in this region; Figure 24 shows the Mach number contours. Comparison 

of the two shows that, although the shock is oblique 

to the grid, it is everywhere captured across two cells. 

The K(’) vectors (related to the pressure-gradient) and 

the vectors (related to the strain-rate tensor) are 

shown in Figures 25 and 26. These give rise to the convection directions shown in Figures 27-29. The first 

shows the convection direction for the shear variable 

and the entropy variable; this is simply the stream direction. The remaining two show the convection directions for the acoustic-like variables. The shocks are evident in the acoustic directions. Convergence histories 

for this case are shown in Figure 30. The two different 

convergence rat.es are for: 

1. constant At (no preconditioning); 

2. a different At in each cell (spatial preconditioning). 

M, = 1.75 10% bump 

-l.OO! , , , , I , , I I 

-1.00 -0.60 -0.20 0.20 0.60 1.00 I 

Figure 19: M, = 1.75 Case - Grid 

M, = 1.75 10% bump 

Mach Number Contoum 

1.00 

1 10000 

1 2 1 1000 I 3 1 2000 

0.60 

0.20 

Y 

-0.20 

-0.60 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

13000 

14000 

15000 

16000 

17000 

1 do00 

19000 

? 0000 

2 1000 

? 2000 

2 3000 

2 4000 

2 5000 

2 6000 

Id 2 7000 

19 2 8000 

-1 00 20 2 goon 

-1 00 -0 60 -0 20 0 20 0 60 100 z 

Figure 20: hi, = 1.75 Case - hlach Number Contours 

12 
M, = 1.75 10% bump 

Mach Number Contourn 

1 1.0000 

2 1.1000 

3 1.2000 

0.60 

0.20 

Y 

-0.20 

-0.60 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

1.3000 

1.4000 

1.5000 

1.6000 

1.7000 

1 .a000 

1.9000 

2.0000 

2.1000 

2.2000 

2.3000 

2.4000 

2.5000 

17 2.6000 

18 2.7000 

19 2.8000 

-1.00 20 2.9000 

-1.00 -0.60 -0.20 0.20 . 0.60 . 1.60 z 

Figure 21: M, = 1.75 Case - Grid-Biased Upwind 

Results (Limiter Off) 

1.0( 

0.6( 

0.2c 

Y 

-0.20 

-0.60 

-1.00 - 

M, = 1.75 10% bump 

lach Number Contours 

1 1.0000 

2 1.1000 

3 1.2000 

4 1.3000 

5 1.4000 

6 1.5000 

7 1.6000 

9 1.8000 

10 1.9000 

11 2.0000 

12 2.1000 

13 2.2000 

14 2.3000 

15 2.4000 

16 2.5000 

a 1.7000 

17 2.6000 

Id 2.7000 

19 2.8000 

20 2.9000 

-0.60 -0.20 0.20 0.60 1.00 

M, = 3.50 10% bump 

Grid 

0.50 

0.30 

0.10 

Y 

-0.10 

-0.30 

-0.50 

-0.50 -o-.30 -0.10 . o.io . 0.30 ’ 0.50 0 

Figure 23: M, = 1.75 Case - Grid 

M, = 1.75 10% bump 

Mach Number Contours 

0 50- 

0 30 

0 10 

Y 

-0 10 

-0 50 

-050 -030 -0 10 0 10 030 C z 

1 10000 

2 11000 

3 1?000 

4 13000 

5 14000 

6 15000 

7 16000 

8 17000 

9 1 8000 

10 19000 

11 2 0000 

12 2 1000 

13 2 2000 

14 2 3000 

15 2 4000 

16 2 5000 

17 2 6000 

18 2 7000 

19 2 (3000 

30 2 9000 

1 

Figure 22: M, = 1.75 Case - Grid-Biased Upwind 

Results (Limiter On) Figure 24: M, = 1.75 Case - Mach Number Contours 

13 
M, = 1.75 10% bump M, = 1.75 10% bump 

Ri VcctOn Entropy and Shear Convection 

0.50 

0.30 -l 

0 10 

Y 

-L I 

-0 

-0 

-0 

1.10 - 

- 

1.30 

-0.50 -0.30 -0.10 0.10 0.30 0.50 z 

Figure 27: M, = 1.75 Case - Convection Directions 

for Shear and Entropy Figure 25: M, = 1.75 Case - dl) Directions 

M, = 1.75 10% bump 

"2 Vectors 

-0.50 -0.30 -0.10 0.10 030 0.50 I 

M, = 175 10% bump 

First Acoustic Convection 

0 50 

0 30 

0.10 

Y 

-0 10 

-0 30 

-0.50 

-650 ' -0'30 . -0'10 . oio . 030 . 050 X 

Figure 28: M, = 1.75 Case - Convection Directions 

for First Acoust ic-like irariab le Figure 26: .Zf, = 1.75 Case - I&') Direct.ions 

14 
M, = 1.75 10% bump 

-0.50 . -0.30 ' -6.10 . o.io . 0.30 . 0.50 t 

Figure 29: M, = 1.75 Case - Convection Directions 

for Second Acoustic-like Variable 

M, = 1.75 10% bump 

Effect of preconditioning 

-1.00 

- none 

-2.00 1 A- h 

\ 

\ 

\ 

-6.00 IIIIIIII 

0. 200. 400. 600. 800. 1000 

Iteration 

' 

Wave-preconditioning was also tried, but the lack ofa 

precise control for the Adk) in regions where low wave 

speeds occur (inside steady shocks) made this calculation actually converge slightly more slowly than the 

calculation with spatial preconditioning only. 

The results for the second case are shown in Figures 31-41. Figure 31 shows the Mach number contours 

for the steady flow. The shock pattern is similar to that 

of the first case, but the shocks move at a shallower angle, so that they do not reflect from the walls before 

reaching the outflow boundary. Due to the strength 

of the shocks in this case, some extra damping was 

necessary to capture them without large oscillations; 

this was provided by smoothing the residuals after each 

time-step with a biharmonic operator. Comparison results, non-limited and limited, are shown in Figures 32 

and 33. Again, the new cell-vertex scheme gives virtually the same results as the non-limited cell-centered 

scheme. Figure 33 shows the substantial loss in resolution caused by turning on the limiter. The remaining 

plots again show the details of the shock-intersection region. Figure 34 shows the grid in this region; Figure 35 

shows the Mach number contours. The KC') and n(') 

vectors are shown in Figures 36 and 37, the convection 

directions in Figures 38-40. Convergence histories for 

this case are shown in Figure 41. The three different 

convergence rates are for: 

1. constant At (no preconditioning); 

2. a different At in each cell (spatial preconditioning); 

3. a different At for each convection equation (wavepreconditioning). 

Due to the relatively high Mach number, the wave 

speeds do not exhibit much of a spread, so that the 

preconditioning gives only about a 20% gain. In this 

case, the convergence history clearly resembles that of 

the scalar equation: after 150 iterations, the residual 

drops rapidly. 

Conclusions 

A genuinely multi-dimensional upwind-differencing 

cell-vertex scheme for a two-dimensional convection 

equation has been formulated, analyzed and tested. 

It has been extended to the Euler equations, giving a 

third-order accurate conservative scheme. The numerical results thus far are promising. The four essential 

elements of the Euler scheme are: 

Figure 30: M, = 1.75 Case - Effect of Local At 

(CFL=0.5) 

15 
M, = 3.50 10% bump 

Mach Number Contours 

1 2.6500 

2 2.7250 

3 2.8375 

M, = 3.50 10% bump 

Mach Number Contours 

1 2 6500 

2 2.7493 

3 2 3487 

0 60 4 2 9375 060 4 2 9480 

5 3 0474 

6 3 1467 

7 3 2461 

0 20 833375 0 20 8 3 3454 

10 3 5441 

11 3 6434 

13 3 8421 

14 3 9414 

15 4 0408 

-0 60- 16 4 1401 

5 3 0375 

6 3 1375 

7 3 2375 

Y 9 34375 y 9 3 4447 

-0 20 12 3 7375 -0 20 12 3 7423 

10 3 5375 

11 3 6375 

13 3 8375 

14 3 9375 

15 4 0375 

17 4.2375 

18 4.3375 

19 4.4375 

-1.00 -1.00 I -0.60 -0.20 0.20 0.60 1.00 20 4.5375 

17 4 ?395 

18 4.33% 

19 4.4382 

- 1.00 20 4 5373 

-1.00 -0.60 -0.20 0.20 0.60 1.00 I 

Figure 31: M, = 3.50 Case - Mach Number Contours Figure 33: M, = 3.50 Case - Grid-Biased Upwind 

Results (Limiter On) 

M, = 3.50 10% bump 

l,oo ,hh& Number Contours , 

1 

0.60 

0.20 

Y 

-0.20 

I -0.60 

1 2.6500 

2 2.7493 

3 2.d487 

4 2.94do 

5 3.0474 

6 3.1467 

7 3.2461 

9 3.4447 

10 3.5441 

11 3.6434 

12 3.7428 

13 3.8421 

14 3.9414 

15 4.04Od 

16 4.1401 

17 4.2395 

18 4.33- 

19 4.4382 

a 3.3454 

a0 4.5375 

-1 00 -060 -020 020 060 100 z 

Figure 32: M, = 3.50 Case - Grid-Biased Upwind 

Results (Limiter Off) 

M, = 3 50 10% bump 

Grid 

0.50 

0 30 

0 10 

Y 

-0 10 

-0 30 

-0 50 

-0 50 -030 -0 10 0 10 030 0 50 

I 

Figure 34: .ti, = 3.50 Case - Grid 

16 
M, = 3.50 10% bump 

Mach Number Contoura 

M, = 3.50 10% bump 

62 Vectors 

0.50 

1 2.6500 

2 2.7250 

3 2.8375 

0 30 4 2.9375 

5 3.0375 

6 3.1375 

7 3.2375 

0.10 a 3.3375 

Y 9 3.4375 

10 3.5375 

11 3.6375 

-0.10 12 3.7375 

13 3.8375 

14 3.9375 

15 4.0375 

-0.30 16 4.1375 

17 4.2375 

Id 4.3375 

19 4.4375 

-0.50 20 4.5375 

0.00 0.20 0.40 0.60 O.dO 1.00 E 

0 

Y 

- 0. 

Figure 35: M, = 3.50 Case - Mach Number Contours Figure 37: M, = 3.50 Case - d2) Directions 

M, = 3.50 10% bump 

61 Vectors 

E 

M, = 3.50 10% bump 

Entrow and Shear Convection 

0 

0 

Y 

-0 

0.50 

0.30 

0 

-0 

-0 

0 

.30 

.50 

0.00 . 010 . 040 . 060 . 0.20 . 1 bo 

Figure 38: M, = 3.50 

for Shear and Entropy Figure 36: irt, = 3.50 Case - Jl) Directions 

E 

Case - Convection Direct ions 

17 
M, = 3.50 10% bump 

0. 

Y 

-0. 

First Acoustic Convection 

0.50 

0.30 

0 

0 

-0.30 

-0.50 

0.00 0.20 0.40 0 60 O.dO 100 b -1.00 

-2.00 Figure 39: M,, = 3.50 Case - Convection Directions 

for First Acoustic-like Variable 

-3 oa 

LodReS) 

-4.00 

M, = 3 50 10% bump 

Second Acoustic Convection 

0 50 

0 30 

0 10 

Y 

-0 10 

-0 30 

-0 50 

0 00 0 20 040 0 60 0 c30 100 

I 

-5.00 

M, = 3.50 10% bump 

Effect of preconditioning 

- none -- space 

____ wave 

-6.00 

Figure 41: M, = 3.50 Case - Effect of Local At and 

Characteristic At (CFL=0.5) 

Figure 40: :VIm = 3.50 Case - Convection Directions 

for Second Acoustic-like Variable 

18 
0 the calculation of the convection directions; 

0 the residual calculation; 

0 the distribution step; 

0 the time-step calculation. 

Three of the four warrant further study to produce a 

practical scheme. 

The convection directions chosen here were the ones 

derived by Hirsch et al. [l]. As these are based on 

derivatives of the flow variables, the numerical values 

can be very noisy. One method used to minimize the effect of the noisiness of the Hirsch K’S was to freeze them 

after the residual had dropped two orders of magnitude. 

This improved convergence considerably. Experiments 

with other directions (e.g. the streamline direction) 

for the K’S showed that they could lead t,o faster convergence. Such alternative directions did not work in 

every case, however. Other wave models, such as that 

of Roe [2], remain to be investigated. 

The residual calculation derived here gives thirdorder accuracy everywhere, which is not an advantage 

near shocks. It is not yet clear how to modify the residual formulas in order to ensure monotonicity. In addition, the highest-order terms in this scheme are strongly 

coupled to the choice of the IC’s. These terms turned 

out to be destabilizing in regions where the K’S were 

highly oscillatory. In particular, subsonic cases (not 

shown here) converged only very slowly, and the final 

solutions were not smooth unless a more reliable fourthorder term (e.g. a biharmonic term) was added. 

Not all of the underdamped behavior of the scheme 

can be traced to its high nonlinearity. As can be seen 

from Figure 10, the basic convection scheme does not 

damp any combination of a high spatial frequency along 

one coordinate with a low frequency along the other 

coordinate, if the convection is precisely in one of the 

coordinate directions. This lack of damping is caused 

by the vanishing of either Bo or By. Improvement of 

the convection scheme in this respect requres the introduction of additional finite-difference terms; these may 

actually be formulated as a smoothing term following 

the distribution step. 

Finally, the time-step calculation, aimed at achieving 

optimal convergence, is far from robust. The technique 

of preconditioning by calculating a value of At for each 

convection equation, at each cell-vertex, can lead to 

large improvements in convergence [13,17]. When one 

of the convection speeds is very small, the potential 

benefit of wave-preconditioning is greatest, but so is the 

danger of taking the time step too large. A satisfactory 

analysis of this remains to be carried out. 

Acknowledgments 

The computer program that provided the cell-vertex 

reults reported above was initially developed during 

~ 

19 

a stay of the authors at the Institute for Computational Mechanics in Propulsion (ICOMP), NASA Lewis 

Research Center, Cleveland, OH. The work was also 

partially funded by NASA Langley Research Center, 

Hampton, VA, under grant NAG-1-869. The authors 

are indebted to Barry Koren and Piet Hemker from the 

Center for Mathematics and Computer Science (CWI) 

in Amsterdam, The Netherlands, for a copy of their 

cell-centered, grid-biased upwind code, the results of 

which were used for comparion with the cell-vertex results. The code was run locally by Eric Hsu, a graduate student in the Department of Aermpace Engineering. Last, but not least, the authors wish to thank Phil 

Roe for numerous tutorials on multi-dimensional upwind differencing, which enabled them to develop the 

present distribution scheme. 

References 

[l] C. Hirsch, C. Lacor, and H. Deconinck, “Convection algorithm based on a diagonalization procedure for the multidimensional Euler equations,” in 

AIAA 8th Computational Fluid Dynamics Conference, 1987. 

[2] P. L. Roe, “Discrete models for the numerical 

analysis of time-dependent multidimensional gasdynasrnics,” Journal of Computational Physics, 

vol. 63, 1986. 

[3] G. Moretti, “The X-scheme,” Computers and Fluids, vol. 7, 1979. 

[4] A. Verhoff and P. J. O’Neil, “A natural formulation for numerical solutions of the Euler equations,” Tech. Rep. MCAIR 83-031, McDonnell Aircraft Company, 1983. 

[5] S. F. Davis, “A rotationally-biased upwind difference scheme for the Euler equations,” Journal of 

Computational Physics, vol. 56, 1984. 

[6] D. Levy, I(. G. Powell, and B. van Leer, “Implementation of a grid-independent upwind scheme 

for the Euler equations.” Submitted for the AIAA 

9th Computational Fluid Dynamics Conference, 

1989. 

[7] P. L. Roe, “Error estimates for cell-vertex solutions of the compressible Euler equations.” ICASE 

Report 87-6. 1987. 

[8] A. Jameson, “A vertex based multigrid algorithm 

for three-dimensional compressible flow calculations,” 1086. Presented at t,he AME Symposium 

for Numerical Methods for Compressible Flow. 

[9] Ii. G. Powell and E. 11. Murman, “An embedded 

mesh procedure for leading-edge vortex flows.” in 

Proceedings of the Transonic Symposium, 1988. 
[lo] M. G. Hall, “Cell-vertex multigrid schemes for solution of the Euler equations,” in Numerical Methods for Fluid Dynamics II (K. W. Morton and 

M. J. Banes, eds.), Oxford University Press, 1986. 

[ll] R. H. Nil “A multiple-grid scheme for solving the 

Euler equations,” AIAA Journal, vol. 20, 1981. 

[12] P. L. Roe, “A basis for upwind differencing of 

the two-dimensional unsteady Euler equations,” in 

Numerical Methods for Fluid Dynamics lI, Oxford 

University Press, 1986. 

[13] B. van Leer, C. H. Tai, and E;. G. Powell, “Design of optimally-smoothing multi-stage schemes 

for the Euler equations.” Submitted for the AIAA 

9th Computational Fluid Dynamics Conference, 

1989. 

[14] D. R. Lindquist and M. B. Giles, “A comparison 

of numerical schemes on triangular and quadrilateral meshes.” Presented at the Eleventh International Conference on Numerical Methods in Fluid 

Dynamics, 1988. 

[15] P. W. Hemker and B. Koren, “Defect correction 

and nonlinear multigrid for the steady Euler equations,” Center for Mathematics and Computer Science Note NM-N8801, Amsterdam, 1988. 

[16] G. D. van Albada, B. van Leer, 

and J. W. W. Roberts, “A comparative study of 

computational methods in cosmic gas dynamics,” 

Astronomy and Astrophysics, voi. 108, 1982. 

[17] B. van Leer, W. T. Lee, and I<. G. Powell, “Sonicpoint capturing.” Submitted for the AIAA 9th 

Computational Fluid Dynamics Conference, 1989. 

20 
Naioonal Aeronautics and 

Space Admlnlstration 

87. Key Words (Suggested by Author(s)) 

Two-dimensional Euler equations 

Cell-vertex schemes 

Downwind weighting 

Report Documentation Page 

18. Distribution Statement 

Unclassified - Unlimited 

Subject Category 64 

2. Government Accession No. NASA TM- 102029 

ICOMP-89- 13 

1. Report No. 

9. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No of pages 

Unclassified Unclassified 21 

4. Title and Subtitle 

A Genuinely Multi-Dimensional Upwind Cell-Vertex 

Scheme for the Euler Equations 

22. Price' 

A03 

7. Author(s) 

Kenneth G. Powell and Bram van Leer 

9. Performing Organization Name and Address 

National Aeronautics and Space Administration 

Lewis Research Center 

Cleveland, Ohio 44135-3191 

2. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Washington, D.C. 20546-0001 

3. Recipient's Catalog No. 

5. Report Date 

May 1989 

6. Performing Organization Code 

8. Performing Organization Report No. 

E-4772 

10. Work Unit No. 

505-62-2 1 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 

5. Supplementary Notes 

This material was presented at the 27th Aerospace Sciences Meeting sponsored by the American Institute of Aeronautics and 

Astronautics, Reno, Nevada, January 9-12, 1989 (AIAA Paper No. 89-0095). Kenneth G. Powell and Bram van Leer, Dept. of 

Aerospace Engineering, The University of Michigan, Ann Arbor, Michigan 48 109 and Institute for Computational Mechanics in 

Propulsion, Lewis Research Center (work funded under Space Act Agreement C99066G). Space Act Monitor: Louis A. Povinelli. 

6. Abstract 

A scheme for solving the two-dimensional Euler equations is developed. It is based on a new scheme for the twodimensional linear convection equation, and the Euler-equation decomposition developed by Hirsch et al. [ 11. The 

scheme is genuinely two-dimensional. At each iteration, the data are locally decomposed into four variables, 

allowing convection in appropriate directions. This is done via a cell-vertex scheme with a downwind-weighted 

distribution step. The scheme is conservative, and third-order accurate in space. The derivation and stability 

analysis of the scheme for the convection equation, and the derivation of the extension to the Euler equations are 

given. Preconditioning techniques based on local values of the convection speeds are discussed. The scheme for 

the Euler equations is applied to two channel-flow problems. It is shown to converge rapidly to a solution that 

agrees well with that of a third-order upwind solver.