NASA Contractor Report 178038

lCASE REPORT NO. 85-59

ICASE

ACCURACY OF SCHEMES FOR THE EULER EQUATIONS

WITH NON-UNIFORM MESHES

NASA-CR-178038

19860009595

E. Turkel

S. Yaniv

U. Landau = ~- '-. ~- c. i_~ "

Contract Nos. NASI-17070 and NASI-18107

December 1985

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

Natlonal Aeronautics and

Space Administration

Langley Research center

Hampton, Virginia 23665

ACCURACY OF SCHEMES FOR THE EULER EQUATIONS NITH NON-UNIFORH MESHES

E. Turkel

Institute for Computer Applications in Science and Engineering

and

Tel-Aviv University

S. Yaniv and U. Landau

IAI, Israel

ABSTRACT

The effect of non-uniform grids on the solution of the Euler equations is

analyzed. We consider a Runge-Kutta type scheme based on a finite volume

formulation. We show that for arbitrary grids the scheme can be inconsistent

even though it is second-order accurate for uniform grids. An improvement is

suggested which leads to at least first-order accuracy for general grids.

Test cases are presented in both two- and three-space dimensions.

Applications to finite difference and implicit algorithms are also given.

Research for the first author was supported in part by the National

Aeronautics and Space Administration under NASA Contract Nos. NASI-17070 and

NAS1-18107 while he was in residence at ICASE, NASA Langley Research Center,

Hampton, VA 23665-5225.

I. INTRODUCTION

In recent years much progress has been made in the solution of the steadystate Euler equations in both two and three dimensions. For complex shapes

these calculations are usually based on a body-fitted curvilinear grid. For

general three-dimensional bodies it is very difficult to construct a bodyfitted coordinate system, see e.g., [3]. In particular, the grid system used

frequently is not smooth. In some cases there may exist discontinuities in

the gradients of the grid while in other cases the gradients vary sharply but

not discontinuously.

Figure1
-2-

In Figure 1 we present a typical grid constructed by FL057 for a wing-alone

configuration. This is an expanded view of the trailing edge for one slice in

the spanwise direction. We see that even in this simple case there are large

variations between neighboring cells.

2. SPACE DISCRETIZATION

We consider the Euler equations for an inviscid fluid written in

conservation form. For simplicity of representation we shall consider twodimensional flow, however, all the results generalize to three dimensions in a

straightforward manner. We thus consider

wt + fx + gy 0. (i)

where w = (p, pu, pv, E)t. Rewriting (i) in integral form, we get

ff w dv + f F.nds = 0 (2)

V S

where F = (f,g)t and n is the outward normal. Letting _ be the cellaveraged values of w, we can rewrite (2) as

V -_ + f + f + f + f (fdy - gdx) = 0 (3) AB BC CD DA

for the zone shown in Figure 2. At this stage, (3) is still exact. In order

to solve (3), we now replace the integrals in (3) by an integration rule.
--3 --

C F

Figure 2

Midpoint Integration Rule

Using the formulation presented in [2] we assume that the metrics are

defined at the mesh nodes while the dependent variables are cell averages or

else are located in the center of the cell. Replacing the integrals in (3) by

the midpoint rule we get

V _t + QE + QF + _ + QH = 0 (4)

where Q = fay - gAx. We examine this more carefully by considering the

point G = (i+I/2, j), then
-4-

QG = fG(Yc - YB ) - gG(Xc - XB)' (5)

We are then left with the job of evaluating fG, gG" In the original

formulation [2], w was calculated at G by

wG = _- wi,j + wi+ l,j , (6)

and then fG = f(wG), gG = g(wG)" Thus, both cells contribute equally to the

value of w along the face independent of the geometry of the grid. Thus, we

do not account for non-uniform volumes and stretchlngs and deformations in the

various directions. Another possibility is to average f and g directly.

fG = _IIfi+l,j fi,j )" 0a

These two procedures are equivalent to within second-order accuracy. There

are some minor differences with respect to shock resolution. In order to

check the accuracy of (6) we first consider a one-dimenslonal example.

The one-dimensional equivalent of (4) - (5) (see Figure 3) is

d xj+ I/2

d-_ f wdx + fj+l/2- fj_l/2= 0. (7)

xj_ 1/2

wj_ 1 wj wj+ 1

I I I I I I I

xj_3/2 xj_ 1/2 xj+I/2 xj+3/2

hi_ 1 hj _ _ hj+l

Figure 3
-5-

Using the averaging on f given by (6a) we get

dwj + fj+l - fj-I = 0. (8)

dt 2h.

3

A straightforward Taylor-series expansion shows that (8) is inconsistent

unless

hj+1 - 2hj + hj_1 = 0(h) (9) h.

3

where h = max(hi). Let rj = hj+i/h j. In [6] we define a grid as algebraic J

if rj = 1 + 0(hp) and exponential if rj+i/r j = 1 + 0(h). Thus, (8) is

consistent only if the grid is algebraic, and is second-order accurate only if

we also have p = 2.

We next replace (6) by

h_ w$+ 1 + hi+ 1 wj (10)

+

wj+1/2 hj hi+1

and (8) by

dw. f(wj+ I_ ) - f(wj_ I_ ) --_ + = 0. (II)

dt h.

3

This scheme is now first-order accurate and furthermore, it is second-order

accurate for all algebraic grids. In [6] it is also shown how to construct a

second-order finite difference version for exponentially stretched meshes.

Until now we have discussed the accuracy of the flux differences. In the

Runge-Kutta schemes [2,5], it is also necessary to add an artificial
-6-

viscosity. This consists of two parts: one is an approximation to a second

difference with a nonlinear coefficient that is of higher order in smooth

regions. This essentially enforces an entropy-like condition and ensures the

correct shock conditions. The second portion of the viscosity is an approximation to a fourth difference and is a linear term. This term is needed in

smooth regions to prevent decoupling the odd and even points and to suppress

the highest frequencies [5]. These viscosities are now also used in implicit

A.D.I. algorithms [4] and also in Lax-Wendroff type methods [I]. Both these

viscosities add third-order terms to the truncation error when the grid is

sufficiently smooth. However, for rapidly varying grids these viscosities can

reduce the order of the method. In particular, we have found the lift

coefficient to be sensitive to the coefficients of these viscosity terms.

Implementation - Two-Dimensions

The formula given by (I0) is a one-dlmenslonal formula. In order to

generalize to two dimensions we use this formula in each direction. Using the

original mapping we calculate the position of the center of the cell and also

the center of each face. Let (_,n) be a plane for which the grid is uniform

and (x,y) be the physical plane. We assume that the mesh is given by a

mapping T, (xi,Yj) = T(_i,_j ). Then the center of a cell is defined by the

image of _i+i/2, nj+I/2) while the center of the faces are given by the

image of (_i±I/2, nj) and (_i' nj±I/2)" We note that using the explicit

mapping of [3] the mapping must be adjusted so that the center of the outermost cell is defined and is not at infinity. Having found these positions we

calculate the distance from the center of the zone to the center of each face.
-7-

This gives the distance _ and m shown in Figure 4. Formula (i0) is then

replaced by

_fi+1,j+ mfi,j (12)

fG = _+m

= _ fi+l,j + (l-e)fi,y

with _ = _/(_+m).

Figure 4

We see from (12) that only the ratio _/(_+m) need be known for each zone.

Hence, once the cell centers and surface centers are found, then only one

value need be stored per cell per direction. Thus, in three dimensions we

need to store three three-dimensional vectors and then the positions of the

cell and face centers can be discarded.

In the codes used, FLO5Z, FL053, FL057, FL059 the metrics are given

explicitly by analytic formulae. Hence, it is easy to find the cell and face
-8-

centers. For other mesh generation techniques it is not as easy to explicitly

find these other positions. One possibility is to construct a grid with twice

as many grid lines. Then the nodal positions are found by keeping every other

point. The other grid points on the fine grid give the cell and face centers

on the final grid. Since we do not need to store all this information, this

is not a large cost. A more approximate procedure is to calculate the cell

and face centers by simple averaging from the nodal positions. When the grid

is highly stretched in the outer field one must be careful calculating the

averages in the far field.

We note that the formulas given by (12) are one-dimensional formulas given

along each coordinate line given by the mapping and not along x and y

coordinate lines. Hence, this improvement accounts for non-uniform

stretchings but does not necessarily account for shearing effects between

cells. Furthermore, if we wish a constant solution to be preserved in the far

field, it is necessary to average w by (12) and then calculate f.

Averaging f will destroy constant solutions since the metrics will be

averaged differently in each direction.

Trapezoidal Integration Rule

In the previous sections we replaced each line integral by a midpoint

rule. Thus, (see Figure 2)

C

f fdy- gdx --(fAy - gAx) G. (13) B

For a general mesh this formula is only first- order accurate since G is
-9-

defined as the image of (_i+1/2, _j) and is not necessarily the center of

the line BC. Choosing G to be the center of the line BC would cause trouble

with the interpolation formulas, e.g., (12). An alternative is to replace the

midpoint rule (13) by the trapezoidal rule

f fdy - gdx -- Ay(f(B) + f(C)) - Ax(g(B) + g(C)) • (14)

B

This formula is now second-order accurate for non-uniform meshes. In order to

evaluate f at node, e.g., B, we use billnear interpolation. This interpolation is done in a standard square. Thus, we know the value of w at the cell

centers (xi,j) , (Xl+l,j, Yi+l,y), (xi,j+l, Yi,j+l ) and (Xi+l,j+ I,

Yi+l,j+l ). We now use an isoparametric billnear mapping to map the

quadrilateral given by these points into the unit square in ($,_) space.

Hence,

x = A I + A2 _ + A 3 n + A 4 _n (15a)

Y = B 1 + B2 _ + B3 n + B4 _n. (15b)

We then assume that w satisfies the same mapping, so

w = C1 + C2 $ + C3 n + C4 _n- (15c)

Knowing x, y, w at the four cell centers corresponding to _,_ = 0,I we can

evaluate the coefficients Ai, Bi, C i. We next calculate the value of _ and

that gives the values of x and y for node B. Given that value of
-I0-

and n we calculate w at the node B by (15c). As is standard in finite

element theory these mappingspreserve the continuityof x, y and w between

zones and the interpolation is first-order accurate.

This method is more accurate than the midpoint rule previously discussed

but is also more time consuming. Further, the bilinear transformation needs

to be modified near boundaries.

In the case of a uniform grid the midpoint rule and the trapezoidal rule

give rise to different schemes for advancing w in time. For simplicity of

notation, we assume that the fluxes are averaged rather than w. For a

uniform grid the midpoint rule reduces to

8wij fi+l,j - fi-l,j + gi,j+l - gi,j-i + = O. (16) 8t 2Ax 2Ay

For the uniform grid the trapezoidal rule reduces to

_w 1

__ZY+

_t _ fi+l,j+l + 2fi+l,j + fi+l,j-i

- (fi-l,j+l + 2fi-l,j + fi-l,j-l) 1

+ 8--_ gi+l,j+l + 2gi,j+l + gi-l,j+l

- (gi+l,j-I + 2gi,j-I + gi-l,j-l) 1 = 0. (17)

We now analyze the stability of both (16) and (17). Linearizing the equation

and freezing coefficients, we obtain
-iiwt + Awx + Bwy = 0. (18)

Using a Leapfrog or Runge-Kutta scheme the stability depends on the eigenvalues of the matrix

D = aA + 8B (19)

where

= sin 8, 8 = sin _ for (16),

a = sin e(I + cos _) ( )

2 , 8 = sin _ I + cos 2 8 (for (17)

We consider a general coordinate system (_,_) and denote the Cartesian

coordinates by (x,y). To simplify the analysis we symmetrize D. We then

find (see [5]) that

l w ac bc 0 1

TDT- 1 = ac w 0 0

DO = bc 0 w 0 '

(20)

0 0 0 w

where c is the speed of sound and

a = _Yn - BYe, b = 8x_ - axn, w = aq + Br

with (21)

q = Yn u - x v and r = x v - y u.
-12-

Hence, the eigenvalues of D are given by

d = W, W ± a2 + b2 c . (22)

It follows that the time step restriction for (16) is of the form

At < K (23)

V lq[ + [r[ + L 1 c

where K is a scheme-dependent Constant, and

L I = x_ + x2n + Y_ + Y_ + (x + Y_ n - yn) + f(x_ x + y_ yn).

For an orthogonal mesh, L reduces to

L I = 2.max(x_ r 2 + y_,2 x 2 + y_).

The time restriction for the trapezoidal rule (17) is of the form

At < K (24)

V

]q[ + [r[ + max([q[, 2 [r[) + L2 c

and

27

L2 - 32 LI"

Thus, the time restriction for the trapezoidal rule is slightly less stringent

than that for the midpoint rule. The coupling of modes between mesh points

also differs between the midpoint and trapezoidal rules.
-13-

3. COMPUTATIONAL RESULTS

We first consider the two-dimensional code FL052. This uses a Runge-Kutta

algorithm for an 0 mesh with acceleration by local time step, enthalpy

damping, and residual smoothing [2,5] but without any multigrid techniques.

All calculations were done on a coarse 64x16 mesh. As the mesh is refined the

accuracy of all the methods gets better.

The first case that we consider is flow about a NACA 0012 airfoil with

M = 0.3 and _ = I0°. The resultant flow is subsonic everywhere and so the

solution could be computed using the full potential equation. The lift

coefficient is then found to be about 1.27. In Table 1 we present some

results for this case. In the first column we indicate which case is being

presented. In the second column we present the coefficient of the fourthorder viscosity. In the third and fourth columns we present CL and CD

respectively. In the last column we give the RMS of the density residual

after i000 time steps. All calculations were done with a four-step RungeKutta method with _I = 1/4, _2 = I/3, _3 = 1/2, _4 = i. This gives fourthorder accuracy in time when the coefficients of the differential equation are

independent of time. Since the Euler equations are nonlinear there is only

second-order time accuracy. In any case, the time accuracy is not of

importance since we are considering only steady flows.

Since the flow is subsonic the second-order viscosity does not play an

important role. We see from Table i that accounting for the non-uniform mesh

increases the accuracy of both the lift and drag coefficients. We also see

that for this coarse mesh that the value of VIS4 also greatly affects the

accuracy. All these calculations were done with a distribution of points in

the normal direction as given in the original version of FL052. Other runs
-14-

not shown were done with an exponential-like stretching in the normal

direction. This provides a smoother mesh and less severe stretching in the

far field. As expected the use of (I0) does not improve things as much when

the mesh is smoother. In fact, for a uniform mesh (6) and (I0) coincide. In

all cases computed the use of (i0) improved the accuracy of the solution.

The next case that we consider is flow about a NACA 0012 airfoil with M =

0.85 and _ = 1°. In this case there is a strong shock. Comparisons with

other codes indicates that the CL ~ 0.36 - 0.40 and CD ~ 0.55. As before,

the lift coefficient is underpredicted by the coarse mesh. Increasing VIS4

f

increases CL and also slightly increases CD. As before we find that the

use of the weighting formula (I0) increases the accuracy of both CL and

CD. In this case, we see that the weighted formula also substantially

increases the convergence rate in some cases. As before, using a smoother

stretching in the normal direction increases the accuracy of all the cases,

though using (I0) still marginally increases the accuracy. In this case, we

again see that the value of VIS4 affects accuracy while other computations

show that VIS2 can also have a strong effect on the accuracy of the code.

As the final two-dimenslonal inviscid calculation we consider flow about

an RAE 2822 airfoil with M = 0.75 and _ = 3° . This is also a transonic

case which needs good shock resolution. Runs with finer meshes indicate that

CL ~ I.I0 and CD N .042. In this case the CL is again underpredicted by

10-20% and the use of the weighting formula (I0) again increases the accuracy

as does increasing VIS4. It is to be noted that the fine mesh calculations

are still done with the Euler equations and we are not considering viscous

effects. In this case the use of a smoother exponential-like stretching in

the normal direction again dramatically increases the accuracy but at the
-15-

expense of a slower convergence rate.

We next tested the weighting formulas in three dimensions using FL057 for

a wing-alone configuration. We found that the code still converged for

extreme cases as M = 0.60 and _ = 24° . The code was run with a coarse

mesh of 64xi0×I0 and a finer (though still coarse) mesh of 96x16x20. The

grid generation used the standard geometry routines contained in FLO57.

Typical results are shown in Table 4. The use of small value for VIS4, less

than 0.5, introduced oscillations in the solution and sometimes the method

would not even converge. The use of (10) in eachocoordinate direction allowed

for faster convergence and fewer oscillations in the converged solution. As

VlS4 was increased to 1.0 the oscillations became less important. Also, as

the mesh was refined the oscillations diminished. For the case M = 0.60,

= 24 and the 96x16x20 mesh, the standard method would not converge

unless we set VlS2 = 1.0. Using the weighting formula (i0) the code coverged

with VIS2 = 0.5. Recalculating the midpoints of the cells and cell surfaces

at each time step adds about 30% to the running time. If the distance ratios

in each direction are stored then the increase in running time is negligible

but three additional three-dimensional vectors need to be stored.

4. CONCLUSIONS

Using a Taylor-series expansion one can show that the standard weighting

of the fluxes (6) at interfaces for a finite volume scheme can be even

inconsistent for general meshes. A new weighting formula (10) guarantees that

the scheme is at least first-order accurate for all meshes. As expected,

computations demonstrate that for smooth meshes the differences between (6)
-16-

and (i0) are minimal. When the stretchings are more severe then (I0) can

offer significant advantages, especially on coarser meshes. In addition, in

three dimensions we found that the weighting (I0) can also improve the

convergence properties of the scheme. We also found that the constants

multiplying the artificial viscosity can have a large effect on the accuracy

of the solution, especially on the lift coefficient. Weighting formulas for

nodal schemes and MacCormack schemes are given in [6].
-17-

_FE_N_S

[i] Davis, R. L., Ni, R. H., and Carter, J. E., "Cascade Viscous Flow

Analysis Using the Navier-Stokes Equations," AIAA Paper 86-0033, 1986.

[2] Jameson, A., Schmidt, W. and Turkel_ E., "Numerical Solutions of the

Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepplng

Schemes," AIAA Paper 81-1259, 1981

[3] Jameson, A. and Baker, T. J., "Solution of the Euler Equations for

Complex Configurations," AIAA 6th Comput. Fluid Dynamics Conference,

Paper 83-1919, 1983.

[4] Pulliam, T. H. "Artificial Dissipation Models for the Euler Equations,

AIAA Paper 85-0438, 1985.

[5] Turkel, E., "Acceleration to a Steady State for the Euler Equations,"

Numerical Equations for the Euler Equations," Numerical Methods for the

Euler Equations of Fluid Dynamics, SIAM, 1985, pp. 281-311.

[6] Turkel, E., "Accuracy of Schemes with Non-uniform Meshes for

Compressible Fluid Flows," to appear J. Numerical Mathematics.

[7] Woan , C. J. and Bonnet, E., "An Investigation of Gridding on the

Accuracy of FLO57MM and FLO57MG Euler Wing Solutions, AIAA Paper 85-

4090, 1985.
-18-

TABLE I: Flow about NACA 0012 with M = 0.3, _ = I0°.

0 mesh with 64 x 16 grid.

Weighting VIS4 CL CD Residual

(6) 0.2 1.227 .008 9 x 10-9

0.5 1.240 .007 6 x 10-9

(i0) 0.2 1.264 .002 2 x 10-8

0.5 1.288 -.001 1 x 10-8

finer mesh ~1.274 0

o

TABLE 2: Flow about NACA 0012 with M = 0.85, _ = 1

0 mesh with 64 x 16 grid.

Weighting VIS4 CL CD Residual

(6) 0.2 .2631 .0487 1 × 10-9

0'5 .2738 .0499 1 x 10-5

(i0) 0.2 .2772 .0492 I × I0-I0

0.5 .3006 .0513 2 x 10-8

finer mesh N.36-.40 ~.054
-19-

TABLE 3: Flow about RAE 2822 with M = 0.75, _ = 3 ° .

Weighting VIS4 CL CD Residual

(6) 0.2 .9586 .0395 1 x I0-II

0.5 .9910 .0407 5 x 10-12

(I0) 0.2 1.003 .0414 3 × I0-II

0.5 1.059 .0431 5 × 10-12

finer mesh NI.10 ~.042

TABLE 4: Three-Dimensional Case with M = 0.60, a = 24° .

NO.

Weighting VIS4 CL CD Iterations Residual

(6) 0.25 1.0248 .4170 300 8 x 10-5

1.00 1.1359 .4611 150 7 x 10-6

-5

(I0) 0.25 1.0155 .4132 300 4 x i0

1.00 1.1265 .4571 150 7 x 10-6


1. Report 1110 NASA CR-178038

ICASE Report No. 85-59

I2. Govern~t Accaslon 1110.

4. TItle .nd SubtItle

ACCURACY OF SCHEMES FOR THE EULER EQUATIONS

WITH NON-UNIFORM MESHES

7. Author(s)

E. Turkel, S. Yaniv and U. Landau

5. Report o.t.

December 1985

6. PerformIng Orpanliitlon Code

8: Performing Orpanilltion Report tIIo.

1-------------------------------., 10. Work Unit No.

9. Performing Organization Name and Address

Institute for Computer Applications in Science

and Engineering

Mail Stop 132C, NASA Langley Research Center TT UA ?':lhhC;_C;')')C;

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Washington, D.C. 20546

11. Contract or Grant No.

13. Type of Report lind Period Covered

Contractor Report

14.Spomoring Agency Code

C;f\C; .':l1_R1_nl

15. Supplementary Notes

Langley Technical Monitor:

C. South Jr.

Final Report

16. Abstract'

Submitted to SIAM J. Sci. Stat. J.

Comput.

The effect of non-uniform grids on the solution of the Euler equations is

analyzed. We consider a Runge-Kutta type scheme based on a finite volume

formulation. We show that 'for arbitrary grids the scheme can be inconsistent

even though it is second-order accurate for uniform grids. An improvement is

suggested which leads to at least first-order accuracy for' general grids.

Test cases are presented in both two- and three-space dimensions.

Applications to finite difference and implicit algorithms are also given.

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

accuracy

Euler equations

Runge-Kutta scheme

non-uniform grids

18. Distribution Statement

64 - Numerical Analysis

02 - Aerodynamics

Unclassified - Unlimited

19. Security Oassif. (of this report)

Unclassified

20. ~urity Classif. (of this page)

Unclassified

21. No. of Pages

21

22. Price

AOZ

For sale by theNational Technical Information Service. Springfield. Virginia 22161