.. NASA Contractor Report 178235 

ICASE REPORT NO. 87-6 

1 

ICASE 

c 

ERROR ESTIMATES FOR CELL-VERTEX SOLUTIONS 

OF THE COMPRESSIBLE EULER EQUATIONS 

(SASA-Ca-1782.35) ERBCIi ES'11HAlES FOB 

L€LL-VERTEX 3CLt'llCYS CP 'ItE CCBEhESSIBLE 

EDLER EkUATPC!BS Final hepcrt (BASA) 43 p 

CSCL 01A 

~3/il2 

P. L. Roe 

I Contract No. NASl-18107 

January 1987 

-_ -- ---___ - -- - ~- ~ 

-__ __ 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERIN 

NASA Langley Research Center, Hampton, Virginia 

Operated by the Universities Space Research Association 

23665 

NASA 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 

N87- 18534 

Unclas 

43560 
ERROR ESTENATES FOR CELL-VERTEX SOLUTIONS OF TEE 

COMPRESSIBLE EULER EQUATIONS 

P. L. Roe 

Cranfield Institute of Technology, United Kingdom 

I ABSTRACT 

I The cell-vertex schemes due to Ni and Jameson, et al., have been 

subjected to a theoretical analysis of their truncation error. The analysis 

confirms the authors’ claims for second-order accuracy on smooth grids, but 

shows that the same accuracy cannot be obtained on arbitrary grids. It is 

shown that the schemes have a unique generalization to axisymmetric flow that 

I 

i 

I 

I 

I preserves the second-order accuracy. 

This research was supported under the National Aeronautics and Space 

Administration under NASA Contract No. NAS1-18107 while the author was in 

residence at the Institute €or Computer Applications in Science and 

Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665-5225. 

i 
I 1. Introduction 

i 

I Most Euler codes written for computing inviscid aerodynamic flows have 

made use of the finite-volume formulation. An early example was described by 

I 

I Rizzi and Inouye 111. Essentially, they are generalizations of the secondI 

order accurate Lax-Wendroff [2] algorithm to nonrectangular mesh geometries. 

I Conservation (which is a necessary ingredient toward ensuring convergence to 

correct weak solutions [3]) comes from regarding the numerical representation 

I 

I of the flow as comprising average states within each computational cell. The 

I rules for updating these values employ interface fluxes which transfer the 

I mass, momentum, and energy between neighboring cells, leaving the total quanI tities present unchanged except for boundary efEects. There seems to be 

, little analysis to indicate how accurate these schemes are on nonuniform I 

I 

grids. They are clearly second-order accurate on uniform, rectangular grids, 

I 

l and it has been generally assumed that second-order accuracy will. still be 

I obtained on a "sufficiently smooth" irregular grid. 

I Finite-volume schemes based on central-differencing algorithms have 

become very popular following the work of Jameson et al. [41. Because these 

schemes update the solution via a five-cell rather than nine-cell stencil, I 

1 

I their accuracy is easier to analyze and has been treated theoretically by 

I Turkel [51, whose results have led to modified schemes offering better 

I 

accuracy. Numerical experiments confirming this are reported by Turkel et al. 

[61 

Yet another possibility is the scheme recently introduced by Ni [7]. 

Here the numerical values represent states found at the corners (vertices) of 

the computational cells. Ni originally described his scheme in finitedifference terms, with the values approximating point samples of the continuum 
-2- I 

solution. Subsequently, it has been described in finite-element terms with 

the vertex values implying a bilinear interpolation function within each cell I 

(Davis et al. [8], Lohner et al. [91). This latter view enables a precise 

I 

I 

interpretation of conservation to be given. At a time when various schools of 1 

thought are learning to use each other's language and ideas, we prefer to 

I classify Ni's method under the neutral term of "cell vertex'' scheme. One 

attraction of Ni's scheme (which has been very significantly refined by Hall 

[lo]) is that it seems to handle arbitrary, even completely unsmooth, meshes 

I 

in a very natural way. However, we will show that claims of uniformly secondorder accuracy are based on an incorrect argument. The object of this paper 

is to provide a correct analysis of the accuracy of cell-vertex schemes. 

This is made easier by the fact that these schemes, in a sense, factorize 

the process of updating the solution. Each cell is examined to see if the 

fluxes around it are in balance; if not, some changes are made. The updating 

process at each vertex uses nine vertex states so that, overall, the scheme is 

as complex as Lax-Wendroff. However, the test for balance involves only four 

vertices, making the first stage very amenable to analysis. This first stage 

is the only part to be analyzed here, but this is the only analysis needed to 

determine the accuracy of the steady state which is reached when all cells are 

in equilibrium. That is to say, the integral 

evaluated by the trapezium rule, vanishes around every cell (Figure la). Note 

that there are generally more vertices than cells so that this condition 

leaves the solution undetermined. Boundary conditions need to be added, but 
-3- 

~ carefully. If too many are prescribed, then the integral cannot vanish everyet al. [13] have proposed a scheme which also makes use of flow variables 

vertices outnumber cells by roughly 2:l; so the problem would be overconstrained even without boundary conditions. Instead, Jameson's procedure is 

1 such that the solution ceases to change when the flux integral vanishes around 

I all control volumes such as the one shown in Figure lb. 

I 

The analysis contained in this paper applies equally to Ni's and to 

Jameson's schemes because it evaluates the accuracy of the approximate flux 

integral in terms of the geometry o€ the polygon around which it is 

I evaluated. The claim, mentioned earlier, of second-order accuracy arising 

, from the trapezium rule is fallacious because, as we refine the grid, the conI tour around which we integrate does not remain fixed but shrinks with the grid I 

l and always contains just a few intervals. In fact, it is easy enough to make 

I 

1 

I 

an analysis based on a local two-dimensional Taylor expansion which demonstrates that the integrals are only evaluated to second-order accuracy if the 

polygons meet certain geometric criteria. These criteria are not met by any 
-4- 

triangle or by any quadrilateral except parallelograms. They are satisfied by 

all regular n-gons (n > 3) and in a seemingly sporadic manner by various 

irregular n-gons (n > 4). In practice, the errors will be acceptable .if the 

I polygons have sufficiently small error constants, the formulae €or which are I 

I 

given. 

There are two offshoots to our analysis. One is that we derive a con- I 

sistent finite-element interpretation of the cell-vertex schemes. We show 

that it gives results which differ from those of the finite-difference interI 

pretation by terms of the same magnitude as the truncation error inherent in I 

I 

either. Finally, we consider the application of these schemes to problems 

which are pseudo-two-dimensional, €or example, axially symmetric. We show I 

that out of a family of plausible extensions, only one member is capable of 

retaining second-order accuracy near the axis. 

2. Description of Cell-Vertex Schemes 

We consider steady inviscid two-dimensional compressible flow governed by 

the Euler equations 

where 

where 

aF aG 

-++- -0 - - 

ax aY 

P =- (y-l) p[h - 1 (u2 + v 2 11 

Y 
-5- 

and h is total specific enthalpy, assumed constant. We intend to satisfy 

these equations in an approximation to their integral form 

where the integral is around any arbitrary closed contour. To this end, the 

plane is divided into cells that are either regularly arranged quadrilaterals 

(Figure la) or unstructured triangles (Figure lb). In either case, the 

vertices are numbered, and a fluid state xi, where u = (p, pu, pv), is 

assigned to the ith vertex. Ni’s scheme for quadrilaterals approximates (2.4) 

by trapezium rule integration around each cell. This is 

- 

where the summation is over all sides of the cells, and all differences are 

taken anticlockwise. Note that D has dimensions of the time derivatives of 

the actual conserved quantities (not their densities), so that indicates 

how masses, etc., associated with the points A, B, C, D are to be changed. 

Different time-marching schemes result from specifying different weights with 

which the changes are distributed over A, B, C, D. Ni’s original recipe [7] 

reproduced a type of Lax-Wendroff algorithm, and Hall [lo] has presented a 

Runge-Kutta version. Here, however, we are interested only in the steady 

- D 

state which we assume is characterized by the vanishing of D for every 

cell. In practice, this will only happen if the boundary conditions are 

specified carefully. (For example, on a grid of n x n cells there will be 

n2 equations and (n+1)2 unknowns; hence 2n+l boundary conditions must be 

- 
-6- 

prescribed.) Also, the residuals do not in general vanish near a captured 

shock ill]. This implies that our interest here is chiefly in the accuracy 

obtained in smooth parts of the flow. 

In Jameson’s triangular grid scheme [13], the average flux across each 

line joining two vertices is evaluated as 

1 1 E =- (F + EB) (J’, - yA> - - (G cB) (xB - XA)* -AB 2 -A 2 -A 

Then masses, etc., proportional to E At are added at X -AB and subtracted at Y. When this has been done for all edges in the mesh, t.ien every 

point has been changed by a net amount proportional to the flux across the 

perimeter of the polygon Formed by all triangles meeting at that point. At 

the steady state, we again assume that all these flux integrals vanish - due 

care being taken at the boundary. 

It is the ability to characterize steady state solutions in this rather 

simple way that makes the accuracy of cell-vertex schemes much easier to 

analyze than Lax-Wendroff schemes where the statement of local equilibrium 

typically involves nine fluid states and thirty-two mesh coordinates. The 

statement of equilibrium for a Jameson-type scheme on quadrilaterals involves 

only five fluid states, but analysis of its accuracy needs twenty-four mesh 

coordinates. 

3. Finite-Difference Analysis 

We wish to know how accurately the formula 
-7- 

approximates the integral 

-F I =$Fdy - (3.2) 

taken around the same polygonal contour. At this stage, the values of F 

are thought of as point samples, in the spirit of finite difference methods. 

This is most easily done by assuming that - F is representable by a Taylor 

expansion about some arbitrary nearby origin. We consider each side of the 

polygon in turn. The numerical approximation to the integral is 

- 

(3.3) 

4 + 0 (h 1. 

The leading terms in the exact integral are found by introducing a local 

coordinate 5 varying linearly from 0 to 1 over AB, so that a typical point 

on AB is defined by 

and 
-8- 

F 

-XY 

) F ] + O(h’) . 1 2 2 

+ (YA + YAYB + YB -yy 

Subtracting (3.5) from (3.3) gives the error on AB as 

I 

5 I 

I 

(3.5) 

Hence the total error may be arrived at simply by summing the contribution from all sides. Clearly we can find the error in the - G integral 

similarly. The final result is 

/Fdy=lFAy-E F +2E F -E F - - 2 -xx 3 -xy 4 -yy 

/G~X=~GAX-E~~~+~E - G -E G 

- - 2 -xy 3 -YY 

(3.6a) 

(3.6b) 
-9- 

I 

where i El = iij- 1 1 (Ax) 3 

= - 1 1 (AX) 2 Ay E2 12 

= - 1 1 AX (Ay) 2 E3 12 

E4 = iij- 1 1 (AYl3- 

(3.7) 

It is pleasing that the error constants take these simple forms and useful that there are only four of them. Observe that the quantities we wish to 

I approximate are O(h2) and that, superficially, the error constants are 

O(h 1, making the method first-order accurate. However, some cancellation 

I 

I 

I 

I 3 

I takes place inside the summations in (3.7), which holds out the prospect of 

I better accuracy, depending on geometric properties of the grid which we now 

I investigate. 

I 4. Polygons with Vanishing Error Constants 

i It is easy to see that any polygon with central symmetry will have El 

I 

to E4 identically zero because contributions from opposite sides will 

vanish. Also, given one error-free polygon, we can form others by permuting 

the sides and forming mirror images. It is straightforward to show that an 

error-free polygon remains error-free under rigid rotations. Certain other 

results are also obtainable. 
-10- 

Theorem I : A quadrilateral is error-free iff it ,s a parallelogram. 

Proof: For a quadrilateral ABCD, the error constants can be manipulated I 

where 

and S is the area of the quadrilateral. The only terms which may vanish in 

these expressions are ’X +x ’X >, (y -y +y -y >. These quantities, 

divided by two, are the components of the vector which joins the midpoints of 

the two diagonals. Thus El to E4 vanish iff the diagonals bisect each other, 

which only happens if the quadrilateral is a parallelogram. 

(xA B C D ABCD 

Corollary: No triangle is error-free. 

Proof: Since triangles are degenerate quadrilaterals, they are excluded 

by Theorem 1. 
i -1 1- 

Theorem 2: Any regular n-gon (n > 3) is error-free. 

I 

Proof: For n even the result is trivial by symmetry. A proof, which 

includes n odd, follows by considering the polygon in the complex plane. I 

I 

Without loss of generality, one side can be made parallel to the x-axis, and 

each side is then represented by one of the nth roots of unity, say 

zk (k = 0, ... n-1) where 

2nk 2n k z = COS - + i sin - 

k n n 

then 

2sk sin2 2ak 3 2.nk s-n - - 3 cos - -- i sin - 

3 2nk 2 2nk 2nk 3 = 1 cos - + 3- cos - 

k k n n n n n n 

1 'k 

2nk 2 2nk 

- 3 cos - (1 - cos -) 2nk sin - 

2 2nk = 1 cos3 = + 31 cos - 

n n n n n k 

- i sin - 2nk (1 - cos 2 -) 2nk 

n n 

= 4(E1 + iE2) 

= (by a similar argument) - 4(E t iE4) - 3 

But 

3k k=n- 1 3 

k=O k=O 

k=n- 1 

1 1 z= 

2;" - 1 

=3 

z -1 1 
-12- 

The numerator vanishes because zn = 1. So the expression equals zero 1 

unless the denominator also vanishes; that is, unless n = 3 (when the expression equals 3). So unless n = 3, E = E = E = E = 0. 1234 

For n > 4, we can find a variety of asymmetric error-free polygons. The 

general procedure is to prescribe (n-2) vertices. Then requlring that El 

to E4 vanish gives four equations for two unknown coordinate pairs. Unfortunately, these are simultaneous cubic equations for which no systematic 

solution is apparent. To give a flavor of the possibilities, we consider a 

pentagon with prescribed vertices at (O,O), (O,l), and (1,O). We maintain 

symmetry by supposing that an error-f ree pentagon exists with additional 

vertices at (l+a,b), (b,l+a). Then we require 

3 33 1 AX)^ = 1 (AY)~ = 1 - (l+a-b) + a - b = 0 

and 

2 1 (Ax)~A~ = 1 (Ay2)Ax = ab2 - (l+a-bI3 - ba = 0. 

Condition (4.2) simplifies to 

3(a-b)(a+l)(b-l) = 0; (4.4) 

if (a-b) = 0, there are no solutions to (4.3). If b = 1, then (4.3) yields 

2 a3+a -a=~ 
-13- 

whereas, if a = -1 then 

2 b3-b -b=O 

1 

2 b = 0, - (1 i 47). 

Alt..ough there appear to be six solutions, Ley actually appear in 

pairs. In each pair, the sides are permutations of the same set of vectors. 

The solution to the problem consists of three independent sets of vectors, 

each of which, by permutation, generates 24 different pentagons, if we allow 

the gtven vectors to permute also. However, out of a given set of vectors, 

only two orderings will produce a convex polygon, that is, those which are 

monotonically ordered by orientation. Those two are mirror images of each 

other. The independent convex pentagons derivable from (4.2) and (4.3) are 

shown in Figure 2. The square is a degenerate case for which one side 

vanishes. 

Error-free polygons can also be constructed by uniting one or more errorfree polygons that have edges in common; but, though the construction of 

error-f ree polygons has some mathematical interest, it seems unlikely that any 

practical mesh could be constructed from them, except for completely regular 

meshes. We turn, therefore, to the more practical question of finding 

polygons for which the error constants are acceptably small. 

5. Polygons with Small Error Constants 

The only complete results that have been obtained relate to quadrilaterals, for which the error constants are given by (4.1). The local trunca- 
-14- 

tion error will be O(h2) provided these expressions are O(h4). We can I 

clearly identify factors which are O(h2), and the only possibility for I 

second-order accuracy is that 

1 (5.la) - XB + xc - XD) = 0 (h 2 1, (xA 

and 

(5.lb) I 

Although these conditions are not met by arbitrary quadrilaterals, they 

are met by quadrilaterals which are produced as part of an analytic grid. 

Suppose we have grid generating functions I 

and the mesh points are the intersections of for integer 

values of m, n. Then, if the mesh is systematically refined, AS, An tend 

to zero, and 

5 = mAS, n = nAn, 

provided the cross-derivat ives are not large. 

Note that since the error terms involve only the local intrinsic geometry 

of each cell, neighboring cells may be of markedly different shape, size, or 

orientation. That statement cannot be made for finite-volume schemes 
-1 5- 

[5,11]. Indeed, the local quality of cell-vertex schemes makes them secondorder accurate on grids generated more straightforwardly than those described 

above. Suppose a coarse grid is determined arbitrarily, and fine grids are 

constructed from it by repeatedly bisecting the cell edges. Within each of 

the origlnal cells, PQRS, there is a mapping defined by 

(5.3a) 

(5.3b) 

which meets the conditions (5.2). 

Related results have been found by Paisley [ll], who finds that in the 

cell-vertex method, it is allowable to perturb a rectangular mesh by O(h2) 

3 whereas with Jameson's [4] scheme, the allowable perturbation is only O(h ). 

6. A Finite-Element Approach 

As mentioned in the introduction, the schemes under study have been 

described elsewhere as finite element schemes. Whether this is justifiable or 

not, there is one objection that can be removed. This objection is that the 

formula 

is obtained by integrating over each triangular (quadrilateral) element 

assuming a linear (bilinear) distribution €or both - F and - G, which is incon- 
-16- 

sistent because there is a nonlinear functional relationship between - F and 

G. A consistent finite-element model can be derived by assuming that the 

vector w varies linearly (bilinearly) over the elements, where 

- 

- w=p "2 (i) . 

These variables were introduced by the author to facilitate the construction of an approximate Riemann solver [14]. They have the property that 

every component of 2, - F, - G is bilinear in the components of - W. Thus 

F= c 

"lW2 '+ hw 2 + - y+l 2 y-1 2 

1 2y w2 - 7g w3 

w2w3 

where h is a constant for steady flow, equal to the specific total 

enthalpy. In computing a typical term of (6.1) such as 

IAB = I, B dy 

there will be subterms of the form 

IAB i,j JA Bw i w j dy . 

The consistent evaluation of this term is 
-17- 

+w w +w,w +2w w 1 = (yB 'A) [2WiBWjB iA jB iB jA iA jA 

1 

- (yB - y )(F 2 A -B The first term contributes the 'naive' approximation +EA); 

the second term is part of a correction that makes the approximation selfconsistent. The leading contributions to it are due to first derivatives of 

- W. The correction term in (6.5) is, to leading order, 

When we sum such terms over every side of the polygon, expressions arise 

proportional to the constants El, E*, E3, E4 given by equations (3.7). 

Where these error constants are small enough, the error due to inconsistency 

is negligible, and both the consistent and inconsistent formulae approximate 

the flux integrals to O(h2). 

7. An Integral Formulation of Axisymmetric Flow 

Flow which is identical in every plane through some axis of symmetry is 

governed by differential equations that can be written in either of the forms 
-1 8- 

u + F + G = (v/r) u - t -x -r - (7.1) 

or 

(r2lt + (rFIx - + (r - GIr = - p (7.2) 

I 

where E, - F, - G are the usual conserved quantities in two-dimensional flow, 

except that r is now a radial coordinate, and v a radial velocity. Here 

- p is the vector (0, 0, p, O)T. 

As a basis for numerical work, (7.1) is tempting because it requires only 

slight modification to a two-dimensional code. However, (7.2) is a better 

I 

basis for shock-capturing because it expresses conservation more exactly. To 

see this, imagine that we are performing a truly three-dimensional computa- I 

tion, using the same grid in every axial plane, so that each 3-D cell is a 2-D 

cell rotated through some angle 9 with respect to the axis (Figure 3). 

The exact angle does not matter, so long as we avoid the extremes Q scd 27. 

If we make a flux balance on such a control volume, the contributions 

from the four faces which do not lie in the radial planes can be evaluated 

analytically, by integrating with respect to the angle. The result is 

a -/{ r u dA - { (rx dr - r - Gdx) - I/ - p dA = 0 

an n 

- 

at n 

(7.3) 

where n is the interior of the 2-D cell, and an its boundary. 

In order to check that a proposed solution was conservative (and hence a 

correct weak solution), we would check that all such integral equations were 

satisfied. A numerical method which imposes this as a constraint will be a 

correct shock-capturing method. In this way, we can give an interpretation to 

conservation of "radial momentum" via the third equation in (7.3). 
-19- 

Equation (7.3) can be the basis for extending to axisymmetric flow either 

of the cell-vertex methods under discussion. In Ni‘s scheme, one integrates 

the quantities rF, -- rG around quadrilaterals; in Jameson’s scheme, the 

quantity r(FAy - - ~AX) is evaluated for each side and used to update the 

adjoining vertices. The problem we now address is how to evaluate these 

things numerically. 

8. Error Analvsis for the Axisvmmetric Methods 

In either method, the central ingredient is an expression of the form 

L Br F dr or lABr - G dx. - 

Consider the first of these. Two plausible computational formulae are 

1 

I AB 

= - 1 (rB - rA) (rA +rB) [F~ + 

‘AB 4 

As in the two-dimensional case, we find what the exact value of the 

integral should be in terms of a Taylor expansion about an arbitrary nearby 

point. The result is 
-20- 

+3(r:+r 1 r +r 2 )F AB B -r 

1 

24 A A AAB BA BAB 

+ - (3r x + 2r x x + rA% + r x 2+ 2r x x 

+ - 1 (3rA 2 xA + 2rArBxA + rB 2 xA + rg 2 XA + rA 2 rB 

12 

2 3 ) F ] + ... . 13 2 + (rA + rA rB + r r + rB -rr AB 

(8.3) 

+ 3r x 2)~ B B -XX 

2r r x + 3rB 2 xB) xxr 

ABB 

I If we make similar expansions for the approximations (8.11, (8.21, and I 

subtract from them (8.3), we find error estimates as follows. 

1 -1 1 3 IAi - IAB = (rB-rA) L#rB-rA)(x -x )F + :(r - rA)-F- B A -X 6 R -L 

- z(~B-~A)(3rA~A 1 + r x - r x - 3rBxB) F AB BA -xx 

- -(r 1 -1- ) (3rAxA - r x + r n - 3rBxB) F 12 B A AB BA -xr 

+ 1 (rB-rA)(rB 22 - r ) F ] + ... A -rr 

xAB 2 - = (rB-rA> [- -(r 1 -r >(x -x F - - 1 (rB-rAl2 F 12 B A B A -X 12 -r 

+ - 1 (xB-xA> (rAxB - r x ) F 

12 B A -XX 

1 

12 B A -xr 

+ - (rB-rA) (rAXB - r X ) F ] + 
-21- 

Comb€ning terms for a particular polygonal control volume, we obtain the 

following equations relating the exact integral to the approximations (8.1), 

(8.2). 

wherein we again €ind the error constants defined in (3.7). At first sight, 

this is a satisfactory result showing that the errors are O(h ) on a smooth 

grid. However, the quantity being approximated is I,rr%dA, which is 

O(rh ). The problem is that, near the axis, we find cells for which r is 

O(h), and the error terms are O(h) relative to the true values. On a nonsmooth mesh, they are even O(1). 

4 

2 

We observe from (8.6) and (8.7) that these particular errors could be 

eliminated by choosing the linear combination 

Note that this is the combination obtained if we assume that K, - F, both 

vary linearly along AB, and then integrate exactly. The error analysis €or a 

single side now gives 

1 L + - (rB - rA) F 1. 24 -r r 
-22- 

Therefore, the total error around a particular polygon is i 

I 

3 1 2 1 2 ~ 

I 

1 IAB - I r Fds = - F 1 ;(Ax) Ar + - F 1 r (Ax) (Ar) - 12 -xx 6 -xr 

(8.10) 

1 3 

12 -rr 

+ - F 1 r (Ar) + ... . 

The analysis is now identical to the two-dimensional case except for the 

appearance of the mean radii (r) in every term. These €actors, however, 

complicate matters considerably. For example, it is no longer true that a 

regular polygon has vanishing error constants because contributions from 

opposite sides no longer cancel in general. To put the error constants into 

forms which make their order of magnitude self-evident needs considerable 

manipulation, which has only been carried out for quadrilaterals. The easiest 

term to analyze then is 

3 4 3 3 4 

A 

- ri) = r - 2r r + 2rBrA - r c (ri + r.1+1) (ri+1 B BA 

43 3 4 

B 

+ r - 2r r + 2rCrB - r C CB 

4 3 3 4 

C 

+ r - 2r r + 2rDrC - r D DC 

4 3 3 4 

D 

+ r - 2r r + 2rArD - r 

A AD 

= 2rA 3 (rB - r ) + 2rg 3 (rC - r + 2rC 3 (rD - r ) + 2rD 3 (rA - rc) D A B 
-23- 

2 2 -rr -r 2 ) 

A AC C + rDrB + rB - r = 2(rA - rc) (rD - rB) (rD 2 

3 

2 B = - (rA - rc) (rB - rD) (rA + r + rC + rD) (rA - rB + rC - rD) 

(8.11) + - 1 (r - rc) (rB - rD> (rA + rB - r - rD> (rA - rB - rC + rD>- 

2A C 

The first term is O(rh4) on grids which are smooth in the sense of 

the previous discussion; it is O(rh3) on non-smooth grids. The second 

term appears to be O(h4). That it shall be O(rh 4 ) imposes a natural condition on the grid near the axis. If the cell nearest the axis has one side 

on the axis, then it can be shown that one of the factors (rA + rB - rC - 

rD), (rA - rB - rC + r,) is O(h); the other is O(rh)--which one is which 

depends on the labelling. It is assumed that this cell is also smooth in the 

same sense as the others. 

We will omit details of the calculations that show the other terms in 

(8.10) to be O(rh 4 ), simply exhibiting the final rearrangements. We find 

1 ;(AxI2Ar = (xB- x,,) (xA - xc> (rA 2 - r + r 2 - r 2, B C D 

(8.12) 

2 2 

+ (xB - xD> (rA - rc (xA - x + xc -xD). B 
-24- 

1 

Here the last two terms are O(rh4) on any smooth grid, and the first I 

I 

one is so if the axis is included. Also, with great patience, we can obtain: ~ 

21 1 r(Ax)(Ar) = (rA + r + rc + rD)(rB - rD>(rA - rC>(XA - xB+ xc - xD) B 

+ s (rA + rB + rc + rD>(rA - rB + rC - rD) 

1 

2 A + - s (r + rB - rc - rD)(rA - rB - rc + rD> 

(8.13) 

where S is the cell area. 

On smooth grids, the first three lines are O(rh4); the fourth is 

O(h6). 

Now we turn to the terms involving - G and - p in equation (8.3). Both 

of these .terms are O(h2), but near the axis they almost cancel so the 

requirement for accuracy is the same as for the other terms. A straightforward evaluation of /Ip dA cannot avoid errors O(h4) so some subtlety is 

needed. 

Integration by parts shows that 

(8.14) 

and we can use this to rewrite the terms in G and 2 as - 

I r - G dx - - P dA = I r (E - p) dx + r 2 dA (8.15) ar 

an n an n 
-25- 

I 

where both terms on the RHS are of order (rh2). Now we can evaluate the 

contour integral by the methods just discussed; the error terms will be 

(8.16) 

The summations multtplying Sr, --rr G have already been shown to be 

O(rh4). The first summation can be evaluated as 

1 + -2 S[3(XA - XB + xc - XD> 2 + (XA - x - x + XD> 2 BC 

+ (x* + XB - xc - XD) 2 3 (8.17) 

Every term in this expression can be justified as on a smooth mesh, 

except for the third line. Here S is O(h2) and the first term inside the 

bracket is and 

the other O(h2). Thus, the overall order of the expression must be h4. 

O(rh4) 

O(h 4 ), but of the two remaining terms, one will be O(r2h2) 

However, this does not destroy the second-order accuracy of the scheme near 

the axis, because in that region (2 - p>xx is O(r) so that all contributions to (8.16) are actually O(rh4). 
-26- 

Lastly, we have to evaluate 

r d~. 

n 

(8.18) , 

A formula for this can be found using finite element techniques, if the I 

gradient is rewritten as 

where J is the Jacobian of the transformation (E ,TI) + (x, y). Since we 

also have dA = J dS dy, (8.18) can be found straightforwardly. The result 

is 

The error term turns out to be 

(8.19) 
-27- 1 

I 

+ - 1 [(rA + r + r + rD)IZs(xA - xB + xC - + 24 BC 

+ (XB - XJ (XA - xc> (rA - rB + r - rD)l C 

(8.20) 

All terms in this expression are O(rh4) provided the grid is smooth and 

contains the axis. The scheme which has emerged, although more cumbersome 

than one would hope for, is the unique cell vertex scheme retaining secondorder accuracy in axisymmetric flow. It may be noted that it has also one 

property commonly demanded of a numerical code. It is satisfied identically 

by uniform flow parallel to the axis. For such a flow, - F is a constant, (2 - 

- p) vanishes , and ap/a r vanishes. Under these circumstances, all contributions to the balance equations vanish identically. 

The expressions which have been found as approximations, in effect, to 

aF/ax - and aG/ar - have been derived by other authors from other considerations. Margolin and Adams [15], in the context of a Lagrangian scheme, sought 

I 

numerical approximations that would make the equation 

--- dV - div u (8.21) V dt - 

an identity at the discrete level. Here V is the volume of a moving parcel 

of fluid. The unique solution to this problem gave the differentiation 

operators derived here. Margolin and Adams report greatly improved accuracy 

in their computations. 
-28- 

The same operators were also obtained by Holm et al. [18] as the unique 

operators to preserve in the discrete formulation certain Hamiltonian properties of the differential equations. 

If the scheme were to be used in practice, it may be noted that some 

simplifications of the expressions are possible, at least for quadrilateral 

control volumes. The terms involving F can be made more compact by 

rearrangement. The terms.involving p in (8.19) can be used partly to cancel 

the terms in (2 - p). What emerges from this is the following, which we 

claim is the simplest version of the unique formula 

- 

- 

(8.22) 
-29- 

A striking Eeature is that the approximations to I r - G dx and ll 2 dA 

are precisely those that would have been obtained by f inite-element integrations assuming bilinear variations for G and p. However, neither of those 

terms individually is second-order accurate. The combination is, given suitable conditions on the grid. 

- - 

9. Conclusions 

The cell vertex schemes studied in this paper have been claimed as 

second-order accurate on arbitrary grids. Certainly the results show an 

impressive convergence as the grids are refined [lo]. However, we have shown 

that on arbitrary grids, the local truncation error may be first-order. 

Explicit error constants associated with grid geometry have been derived. 

These indicate what sort of cell shapes to avoid and could also be used as 

part of a mesh refinement strategy. 

The extension to axisymmetric flow has been made, and it seems harder to 

achieve accuracy then. The unique formula which retains second-order accuracy 

near the axis has been derived. However, the extension to fully threedimensional schemes is an extremely difficult piece of algebra which has not 

yet been attempted. 
-30- 

REFERENCES 

[l] A. Rizzi, M. Inouye, "A time-split finite-volume technique for three- , 

I 

: dimensional blunt-body flow." AIAA Journal, No. 11, 1973, pp. 1478- 

1485. 9 

[2] P. D. Lax, B. Wendroff, "Difference schemes for hyperbolic equations 

with high order of accuracy." Comm. Pure and Appl. Math., Vol. 17, 

1964, pp. 381-398. 

[3] P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their I 

numerical computation." Comm. Pure and Appl. Math., Vol. 7, 1954, pp. 

159-193. 

[4] A. Jameson, W. Schmidt, E. Turkel, "Numerical solutions o€ the Euler 

equations by finite volume methods using Runge-Kutta time-marching 

schemes." AIAA paper 81-1259, 1981. 

[5] E. Turkel, "Accuracy of schemes with nonuniform meshes for compressible 

fluid flows." Appl. Numer. Math., in press. 

[6] E. Turkel, S. Yaniv, U. Landau, "Accuracy of schemes for the Euler 

Equations with nonuniform meshes .I1 ICASE Report No. 85-59, 1985. Also 

AIAA paper 86-0341. 
-31- 

I 

! 

I 

j 

I 

I 

I 

1 

I 

I 

~ 

I 

I 

I 

I 

I 

1 

I 

I 

I 

I 

R-H. Ni, "A multiple-grid scheme for solving the Euler equations." - AIAA 

Journal, Vol. 20, No. 11, 1982, pp. 1565-1571. 

R. L. Davis, R-H. Ni, W. W. Bowley, "Prediction of compressible laminar 

viscous flows using time-marching control volume and multiple grid 

technique." AIAA Journal, Vol. 22, No. 11, 1984, pp. 1573-1581. 

R. Lohner, K. Morgan, J. Peraire, 0. C. Zienkiewiez, "Finite-element 

methods for high-speed flows." AIAA paper 85-1531-CP, 1985. 

M. G. Hall, "Cell vertex multigrid schemes for solution of the Euler 

equations." In Numerical Methods for Fluid Dynamics 11, K. W. Morton, M. 

J. Banes (eds.), Oxford University Press, 1986, pp. 303-345. 

M. F. Paisley, "A comparison of cell centre and cell vertex schemes for 

the steady Euler equations." Oxford University Computing Laboratory, 

Numerical Analysis Group, Report 86/ 1, 1986. 

K. W. Morton, M. F. Paisley, "On the cell-centre and cell-vertex 

approaches to the steady Euler equations and the use of shock 

fitting." 10th Int. Conf. Num. Meth. Fluid Dyn., Beijing, 1986. 

A. Jameson, T. J. Baker, N. P. Weatherill, "Calculation of inviscid 

transonic flow over a comlete aircraft." AIAA DaDer 86-0103. 1986. 
-32- 

[14] P. L. Roe, "Approximate Riemann solvers, parameter vectors and 

difference schemes." J. Comput. Phys., Vol. 43, 1981, pp. 357-372. 

[15] L. F. Margolin, T. F. Adams, "Spatial differencing for finite diEference 

codes." Report LA-10249, Los Alamos National Laboratory, January 1985. 

(161 D. D. Holm, R. A. Kupershmidt, C. D. Levermore, "Hamiltonian 

differencing of fluid dynamics." Advances in Applied Mathematics, Vol. 

6, 1985, pp. 52-84. 
-33- I 

I 

APPENDIX 

Some ProDerties of the Jameson Formula 

At convergence, Jameson’s scheme for triangles [13] satisfies 

around the perimeter of every polygon formed from the set of triangles meeting 

at a common vertex. This has been presented as a set of flux balance 

equations around an overlapping set of control volumes. A straightforward 

rearrangement of (A.l) is 

G(x -X )=o. yi-1) - 7 -i i+l i-1 

Now because if the polygon has an odd number of sides, 

we travel twice round its perimeter to rearrive at our starting vertex, and 

with even sides there are two closed polygonal paths. Therefore, we can add 

arbitrary constants to the 3, 5. Let us choose 3, GJ) which are the 

values at the common vertex (Figure lb). 

c (Yi+1- Yi-1 ) = 0, 

Rewrite (A.2) as 

I 

I The reason for the factors is shown in Figure 4, where P, Q are 

the controids of OABS, OBC. Then I 

I 
-34- 

So (A.3) can be interpreted as the flux balance around non-overlapping 

control volumes whose edges join the centroids. This interpretation leads to 

a dual form of the algorithm as follows. 

Scheme II 

For all edges MN in the mesh, compute 

where S, T are the other vertices of the triangles adjacent to MN, and use 

this flux to update M, N. 

This has an identical outcome to the standard scheme, which we repeat for 

contrast. 

Scheme I 

For all edges MN in the mesh, compute 

and use this flux to update S, T. 

Indeed, there are two more possibilities, both of which will also yield 

the same results at convergence. 

1 

I 

I 

i 

I 

I 

I 

i 

! I

1 

I 

I 

I 

I 

I 

I 

1 

Scheme 1x1 

For all edges MN in the mesh, compute 
-35- 

and use this flux to update My N. 

Scheme IV 

For all edges MN in the mesh, compute 

and use this flux to update S, T. 

Given the interpretation that u is piecewise linear within each triangle, one is tempted to look for control volumes allowing a more accurate 

integration. For example (Figure 5), we may consider the control volume 

formed by joining the centroids of each triangle to the center of each side. 

Along PI the value of F, say, varies linearly from - (F + EB) to 

(In view of 1 

- (F + EA + xB) with average value 34 

the discussion in Section 6, we neglect the nonlinear variations of Some 

algebra will show that this does nothing more than reproduce equations (A.l) 

or (A.3). This result is obvious if we think of integrating over the 

area of the control volumes. Within OAB, is a constant and the area of 

OJPI is one-third that of OAB. 

- 

1 

- 2-0 

1'5% + 2F + 5 F ). 12 -A -B 

- F.) 
-36- 

! 

" (b) 

Figure 1. (a) Part of a structured quadrilateral mesh. 

(b) Part of an unstructured triangular mesh. 
-37- 

Figure 2. Examples of polygons with zero error constants. 
-38- 

r 

Figure 3. A computational cell partly rotated about the axis. 
-39- 

C 

Figure 4. Alternative control volume for Jameson’s scheme. 
-40- 

C 

Figure 5. Control volume using cell bisectors. 
Standard Bibliographic Page 

17. Key Words (Suggested by Authors(s)) 

numerical approximation, compressible 

flow, error estimates 

. Report No. NASA CR-178235 12. Government Accession No. 

18. Distribution Statement 

02 - Aerodynamics 

64 - Numerical Analysis 

Unclassified - unlimited 

ICASE Report No. 87-6 

ERROR ESTIMATES FOR CELL-VERTEX SOLUTIONS OF THE 

COMPRESSIBLE EULER EQUATIONS 

. Title and Subtitle 

19. Security Classif.(of this report) 

Unclassified 

’. Author(s) 

P. L. Roe 

20. Security Classif.(of this page) 21. No. of Pages 22. Price 

Unclassified 42 A0 3 

1. yfo mi rga izatio Name and nsiiWe tor eomputer qepdrperications in Science 

and Engineering 

Mail Stop 132C, NASA Langley Research Center 

Hamton. VA 23665-5225 

2. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Washington, D.C. 20546 

5. Supplementary Notes 

Langley Technical Monitor: 

J. C. South 

Final Report 

Submitted 

3. Recipient’s Catalog NO. 

5. Report Date 

January 1987 

6. Performing Organization Code 

8. Performing Organization Report No. 

87-6 

10. Work Unit No. 

11. Contract or Grant No. 

NASl-18107 

13. Type of Report and Period Covered 

CouUactor Report 

14. Sponsoring Agency Code 

505-90-21-01 

to AIAA J. 

6. Abstract 

The cell-vertex schemes due to Ni and Jameson, et al. have been subjected to 

a theoretical analysis of their truncation error. The analysis confirms the 

authors’ claims for second-order accuracy on smooth grids, but shows that the 

same accuracy cannot be obtained on arbitrary grids. It is shown that the 

schemes have a unique generalization to axisymmetric flow that preserves the 

second-order accuracy. 

For sale by the National Technical Information Service, Springfield, Virginia 22161 

NASA-Langley, 1987