Comparison of several finite difference schemes for time marching methods as applied to one dimensional nozzle flow

van KARMAN

INSTITUTE

FOR FLUID DYNAMICS

TECHNICAL NOTE 132

COMPARISON OF SEVERAL FINITE VIFFERENCE SCHEMES

FOR

'

TIME MARCHING METHOVS AS APPLIEV TO

ONE VIMENSIONAL NOZZLE FLOW

W. VAN HOVE & A. ARTS

JUNE 1979

TECHNISCHE UNIVERSITEIT DELFT LUCHTVAART· EN RUIMTEVAARTTECHNIEK

BIBUOTHEEK

Kluyverweg 1 - 2629 HS DELFT

1 JAN.

1988

~A~

-~O~-

RHODE SAINT GENESE BELGIUM

~VW

~

CHAUSSEE DE WATERL00₇ 72

B - 1640 RHODE SAINT GENESE, BELGIUM

TECHNICAL NOTE 132

COMPARISON OF SEVERAL TIME VIFFERENCE SCHEMES FOR TIME MARCHING METHOVS AS APPLIEV

TO ONE VIMENSIONAL NOZZLE FLOW

W. VAN HOVE & A. ARTS

JUNE 1979

ABSTRACT

The one dimensional, inviscid, transonic flow with shock through a convergent divergent nozzle is used as a test case to compare several explicit finite difference schemes for time marching methods. The comparison was carried out both theoretically and numerically, and accounted for stabi-lity conditions and convergence speed. Additional means to improve the convergence of time marching methods are also discussed.

TABLE OF CONTENTS

Abstract . . . List of symbols List of figures

1. INTRODUCTION

2. ONE DIMENSIONAL FLOW THROUGH A CONVERGENT DIVERGENT NOZZLE . . i i i i v 1 3

2.1 Partial differential problem . . . 3 2.2 Stability of the discretized system 7 2.2.1 Courant-Friedrichs-Lewy stability condition 7 2.2.2 von Neumann stability analysis . . . . 9 2.3 Practical elaboration of the discretization 11 2.4 Euler scheme . . . 12 2.5 Lax scheme . . . 13

2.6 Corrected viscosity scheme . 15

2.7 Lax-Wendroff scheme 18

2.8 Corrected damping method 21

2.9 Influence of the initial conditions 24

2.10 Variable time step 24

2.11 Mesh refinement 26

References Figures . .

28 31

a a· J A

df

cp Cv CFL e E 9 . J § i IC [/ k L M MRS N Nd Nv O(n) p R Re t T u U ,..' VC x X LIST OF SYMBOLS speed of sound eigenvalues ofc7t

section of the one dimensional nozzle

coefficient matrix of the partial differential system specific heat at constant pressure

specific heat at constant density

Courant-Friedrichs-Lewy stability condition interna 1 energy

total energy

eigenvalues of§

von Neumann amplification matrix imaginary constant (i 2

₌

_-1) initial conditions

unity matrix

damping coefficient in the corrected damping metn..nd

-length of the one dimensional nozzle Mach number

Mach number ratio of the shock

number of time steps needed to reach convergence

number of time steps with constant damping correction number of time steps with constant viscosity correction order of magnitude n pressure gas constant Reynolds number time temperature velocity Fourier amplitude of u viscosity coefficient

one dimensional coordinate along the nozzle axis number of discret;zation points

y fit fix E: À W p

e

0 1 2 th j k kin t tin m -+

ratio of specific heats (y

₌

_cp/cv)

time step discretization s te p convergence error wave 1ength ( À

=

wave number density

e

=

wflX Subscripts tota1 in1et outlet throat j

₌

1 ,

...

,

n; n at the point at the point at time 1 ev e 1 at time level maximum va1ue Superscripts dimensiona1 vector with with t

=

t

=

for 27T/W)

=

number of equations coordinate x = k6X coordinate x

=

(k±n)flx t6t (t±n) fit all J

LIST OF FIGURES

Fig. No. Title page No.

1 Computational domain 31

2 Characteristic directions in a boundary point 32

3 CFL stability condition 33

4 Correction factor on the CFL number for

CORstant total enthalpy 33

5 Nozzle geometry 34

6 Calculation result with the Lax scheme 35 7 Convergence of the Lax scheme and the

corrected viscosity scheme 36

8 Stability of the corrected viscosity scheme 37 9 Calculation result with the corrected

viscosity scheme 38

10 Calculation result with the corrected

viscosity scheme 39

11 Influence of initial conditions

initial data 40

12 Influence of initial conditions

comparison of convergence error 41 13 Influence of a variable time step 42 14 Influence of reducing the grid spacing

on the convergence error 43

15 Influence of the mesh refinement on

the convergence error 44

16 Influence of the grid spacing on the

1. INTRODUCTION

Time marching methods become more and more popular for the calculation of inviscid, transonic flows. The main reason for this evolution is of course that the partial differential equa-tions describing a time dependent flow are hyperbolic with respect to time in both subsonic and supersonic flow regions, so that the same finite difference operator can be used over the complete flow field. With the shock capturing technique and appropriate sbheme, shocks will appear automatically during the instationary evolution and no complicated shock fitting is needed.

With a time marching method the stationary flow is cal-calculated as the asymptotic limit af ter a sufficient long time of the solution of the instationary equations, starting from given initial conditions. These initial conditions can be considered as a large perturbation of the final solution, giving rise to pertur-bation waves. Convergence is reached when these perturpertur-bation waves have died out. As the time dependent equations of inviscid flow (the Euler equations) do not have any dissipative effect, the only damping is introduced by the discretization and must remain very small to obtain a sufficient accuracy. Therefore the convergence is very slowand requires many iterations and a long computational time. This demand of computer time is the main disadvantage of time marching methods and all the development work is therefore directed

towards an increase of the convergence speed. However ~ the popu-larity of time marching methods will remain closely linked to the availability of new generations of faster computers.

In order to evaluate their potential, several finite difference sche~.s for time marching methods were applied on the calculation of the one dimensional, inviscid, transonic flow with shock through a convergent-divergent channel. Due to the simplicity of the problem and the fact tbat the correct analytical solution is known, it is possible to make extensive tests on accuracy and convergence speed and to compare theoretical prediction methods with numerical results. This work was carried out preliminary

to the development of two and three dimensional calculation methods for the inviscid, transonic flow through turbomachinery blade rows. Therefore, only explicit finite difference schemes which are easy to program and require a minimum of computer memory storage are investigated. More specifically, we are

interested in the performance of the corrected viscosity scheme, which is in our experimence very adequate for time marching

calculations. Other means to improve the convergence speed of a particular finite difference scheme are discussed as well.

2. ONE DIMENSIONAL FLOW THROUGH A CONVERGENT DIVERGENT NOZZLE 2.1 Partial differential problem

The conservation laws for mass, momentum and energy for the instationary, inviscid, one dimensional flow through a con-vergent-divergent nozzle with section A can be written in dimen-sional for as follows

apÄ

"

_{+ apüA}

=

0

at

ax

apüÄ

a(p+;;ü

2

_)Ä

_-

_dÄ

+

_-

=

P dx

at

ax

( 1 )

dPEÄ

₊

aü{p+pË}Ä

₌

₀

at

ax

p

=

pRT

=

-

e + -

ü

2 = 2

These equations are normalized by the total inlet conditions, the length [ of the nozzle and the throat section A

th

-

ii x x

=

A

₌

1 \ [ _Ä th

t

J

t =

R

TOl [

-

_T p

₌

- L

_-

p

₌

-p- T

=

( 2 ) TOl POl POl

u,a

₌

u,a e,E

₌

e

,f

IR

TOl R TOl

P

=

pT e

=

T E

=

e + u2 a2

₌

_yT

y-1 2

au _{+ af} +

_h

=

at ax

rA

[PUA

1

0 +

puA

1

(p+pu2_{)A ~}

_dAl

_{( 3 )}

u

=

=

=

p _Ox

pEA u(p+pE)A 0

(y-1) p (E- u2

P

=

- )

2

The use of this conservative form is necessary for a correct shock calculation but for theoretical investigations the quasi linear form + + au _+c:!J~

₌

at ax

û

=

[

~

1

is more adequate

h

u

c;If=

(y -1 )-e _p 0 p

~Y-l)l

~

=

u (y-1)e

o

pu dA TOx eu dA -(y-1) T Ox ( 4 ) The energy equation can be replaced by an algebraic relation expressing that total enthalpy remains constant. But this is only true for stationary flow and is not compatible with th~ instationary equation5~ so that the time marching calculation will 'not follow the physics of the flow. It will nevertheless converge to the correct stationary solution (as constant total enthalpy is compatible with this solution) which is our only interest. In this way one can gain approximately 1/3 (1/4 for two dimensional problems) of the calculation time. However~ care should be taken when applying the constant total enthalpy rela-tion, beçause the stability conditions of the discretized equa-tions changes (see § 2.2). With constant total enthalpy the sys-tem to be solved becomes, in conservative form :

ait

₊ at + af +

=

h ax

-+

[::A]

-+

[PUA

]

-+

[: d

u = f = _h = _{( 5 )} '(p+pu2_)A

P

= p(l - -y-l u2 ) 2y

and in quasi linear form

-+ -+ _-+ au _+c:II-~ h = at ax

[:]

_:IY]

[- P

u

dAl

[U

-+ A ë1x -+

c:II--

h ( 6 ) u = = a2_{jyp 0}

Both systems (4) and (6) are hyperbolic and posses respectively 3 and 2 families of characteristics. The extreme characteristic directions in a certain point limit the domain

of influence and the domain of dependenee for that point (Fig. 1). The slope of the characteristics is given by

dx

dt

=

a· J

(7)

with a_j the eigenvalues of the matrixc7f. The expressions for a_j are, for system (4)

a·

=

u, u±a

J ( 8 )

and for system ( 6 )

~ u± j(Y-l) 2 2 1

_a

2

a·

=

u +

-J _{2y 2y y} (9 )

A mathematically well posed problem (the solution exists, is unique and depends continuouslyon the initial and boundary

conditions) for the above described system is obtained by adding appropriate conditions on the limits of the computational domain (Fig. 1). This we" posed character of the problem is a necessary

condition for a stable and convergent solution.

At t

=

0, the complete flow field must be defined. These initial conditions give rise to perturbation waves and convergence is reached when these waves have disappeared. The form of the

initial conditions will be given in § 2.3 and their influence on the convergence speed will be discussed in § 2.9.

For t > 0 certain boundary conditions must be specified on the upstream or inlet (x

=

0) and on the downstream or outlet

(x

=

1) boundary in order to obtain the correct stationary solution. The number of conditions to impose on each boundary can be obtained by counting the number of characteristics coming from outside the computational domain in a point on the boundary (Refs. 7,8). The system (4) for example has three families of characteristic direc-tions with slope given by (7) and (8). If the velocity is subsonic (u

<

a), there are two characteristics with positive and one with negative slope. Figure 2a shows these characteristics in a point on the inlet boundary. The two characteristics with positive slope come from outside the domain and the information they carry must be replaced by two boundary conditions. If the inlet velocity is supersonic, the three characteristics have a positive slope and three boundary conditions must be imposed. At the outlet the same considerations can be made but this time the characteristics with positive slope come from inside the computational domain. The system (6) has only two families (Fig. 2b) of characteristic di-rections, one with positive and one with negative slope if u

<

a and two with positive slope if u > a. The following table summa-rizes the number of boundary conditions to impose in each case.

SYSTEM (4) SYSTEM (6)

INLET OUTLET INLET OUTLET

M<1 2 1 M < 1 1 1

The nature of the boundary conditions is not a problem if the flow at the boundary is supersonic : one has to impose all flow variables if the boundary is the inlet boundary and no flow variables if the boundary is the outlet boundary. Such a super-sonic outlet boundary is an advantage from the point of view of convergence speed as no boundary condition can be imposed and the perturbation waves can leave freely the computational domain. All other boundaries are constrained and reflect the perturbation waves back into the domain. This reflection can be neutral or can damp or even amplify the waves. In reference 7, Essers investigated the interaction between a boundary and a perturbation wave for all kinds of boundary conditions. From this study we recall that the optimum combination of boundary conditions, provoking the most important damping of the perturbation waves, is obtained by fixing total pressure (and also total temperature for system 4) at the inlet together with static pressure or static density at the out-let. But even for this combination of boundary conditions, the damping remains very small.

2.2 Stability of the discretized system

Even when neither the partial differential equations nor the boundary conditions amplify certain perturbation waves, it is possible that the solution of the discretized equations diverges explosively. This instability is due to the discretiza-tion and can therefore only be studied on the discretized equa-tions. Because of the importance of the stability condition for explicit discretiz~tion schemes, we will briefly repeat the two most important methods to analyze stability.

The CFL condition is a necessary but not sufficient condition for the stability of any explicit discretization of a hyperbolic system. Figure 3 shows the discretization grid for a system with two independent variables (x and t). The point Q on time level t is calculated from the points A, B, C on time level

t-2~t, etc. Clearly the numerical domain of dependence of the point Q is limited by QAD and QCH. The physical domain of depen-dence of Q is limited by the extreme characteristics a and c through Q. To obtain a convergent and thus stable calculation, it is necessary that any perturbation that can influence Q phy~ sically must also be able to do this numerically (the discreti-zation may not cut off physically possible influences). This requires the physical domain of dependence to be contained com-pletely in the numerical domain of dependence. In view of the relation (7) for the slope of the characteristics, this leads to the CFL condition

la.1 Ä!~ 1

J m ~x .

(10)

with lajlm the maximum absolute value of the eigenvalues a j . As the eigenvalues a j represent the propagation speed of the perturbation waves along the characteristics, relation

(10)

also states that any perturbation originating in a certain grid point, may not have passed the neighbouring grid point one time step later.

Because the CFL condition

(10)

ts a necessary condi-tion for stability, it gives the maximum allowable time step for any explicit discretization scheme. Since the matrixc7f0f the systems

(4)

and

(6)

depends on the unknowns,

(10)

is alocal condition to be fulfilled in every point. In practice one cal-culates the limiting time step in every point using the local values for lajlm and ~x (for a non uniform grid) and the most.

restrictive time step is selected.

For the system (4) _{the eigeavalues a j are given by}

(8) and the maximum value is

Ipjl

=u+a

m

( 11 )

leading to the classical form of the CFL condition. But using system (6) with constant total enthalpy gives

I

a

.1

=

y+1 u + j(Y-l)'

J

+ 1

i

=

(u+a).f(M) J m 2y 2y y (12a) with y+l M-JtY-l) \ , _{+ -}1 f(M)

=

2y 2y y ~

-

1 M + 1

..;y

(12b)

wh ere the function f(M) is represented in figure 4 and can be approximated by the value l/vr for Mach numbers of practical use. This means that the maximum permissible time step for system (6) is about 15% higher than for system (4), what does

not mean, however, that the calculation will converge faster. As was mentioned in § 2.1, keeping total enthalpy constant is not compatible with the original Euler equations and system (6)

describes in fact another kind of flow with other characteristic directions but the same stationary solution as system

(4).

The perturbation waves propagate at lower speed for system (6)

thus increasing the total time needed to reach convergence. This compensates the increased time step and the total number of

iterations remains the same (provided the appropriate form of lajlm is used for each system).

Performing a Fourier analysis of the solution of a discretized equation

+ t +t iwkllx

uk

=

~ U e (13 )

and reintroducing this in the discretized equations leads to the following relation for the amplitude of every Fourier component

e

=

wllX ( 14 )

If the partial differential system is linear with constant

coefficients, if the computational domain is infinite in extent in x direction or has periodic boundary conditions, if the grid

spacing is uniform (llX

=

ct) and if gis a normal matrix, the van

Neumann necessary and sufficient stability condition is

Ig·1

<

1 + O(llt)

J m for all w

(15a)

where the term of order llt must remain bounded if llt + 0 and

llX

=

f(llt) + O. This term permits an exponentially growing

solu-tion that could be a correct solusolu-tion for certain partial dif-ferential equations. But the systems (4) and (6) do not have such a solution and therefore we can use the following condition

for all w (l5b)

In practice, however, most of the conditions are not fulfilled. When the computational domain is bounded, the von Neumann condition, like the CFL condition, is an internal tta-bility condition that gives no information on the influence of the boundary conditions. When the system of partial differential equat;ons is linear with variable coefficients or non linear or

when llx is not constant

(15)

is alocal stability condition that

is neither necessary nor sufficient. But experience shows that at least for quasi linear systems theresults are reliable.

In reference 7 it is shown that a non derivative term

+

(like the term h in (3)) which does not introduce an amplifica-tion of the perturbaamplifica-tion waves in the partial differential system, also will not introduce instability in the discretized system if the non derivative term is approximated by the mean value of the surrounding points. In this case only the derivative terms must be taken into account when performing a stability analysis.

+

This can be applied to the term h in the systems (3) to (6) since these systems are derived from the Euler equations which do not possess any amplification or damping of the perturbation waves.

2.3 Practical elaboration of the discrètization

Figure 5 presents the geometry of the two nozzles that are used in the numerical investigation. Nozzle 1 is designed to give a linear Mach number distribution before the shock (where the flow is isentropic). Nozzle 2 has a linear variation of the sec-tion. Both nozzles are discretized with 46 points.

Both the system (3) with energy equation and the

system (4) with constant total enthalpy are applied. The x deri-vitive is discretized with a second order cen~ral difference scheme and the non derivative term is approximated by the mean values of the surrounding points.

Two sets of initial conditions are used. The first one is the correct subsonic solution with sonic throat section

(ICI).

The second set is calculated with the isentropic relations from a parabolic pressure distribution

(IC2).

This second set of initial conditions does not satisfy the system of partial differential equations, but has a smooth variation of all flow variables where-as the first set hwhere-as a pressure step in the outlet section. The influence of the initial conditions will be discussed further in § 2.9.

At the inlet total pressure (and also total temperature for system 3) are imposed, but due to the normalization of the equations this reduces to the use of isentropic relations. At the outlet, the normalized outlet pressure P2

=

P2/POI is imposed, fixing shock position and strength.

But thesé boundary conditions which we will refer to as physical boundary conditions are not sufficient to determine all flow variables on the boundaries. Additional conditions are needed to calculate the missing flow variables. These numerical boundary conditions have an important influence on the stability and the accuracy of the calculation and general rules to choose these conditions cannot be given. For the calculation of the

one dimensional flow through a convergent divergent nozzle, good results are obtained setting the inlet mass flow equal to the throat mass flow and applying momentum conservation (and also energy conservation for system 2) between the outlet section and the preceeding section.

In the subsequent chapters all discretization schemes are written in quasi linear and in conservative form. Although the quasi linear form is used for theoretical investigations, the conservative form must be used for practical calculations to

obtain a correct shock capturing.

At every time step (every Nv time steps for the corrected viscosity scheme) a convergence error E is calculated

E

=

MAX

k=l,X

( 16)

The convergence criterium used to compare different calculations is

E

<

0,2%

2.4 Euler scheme

The most straightforward discretization of the time deri-vative is a forward difference operator, giving the Euler scheme. conservative form :

( 17a)

(17b)

A von Neumann stability analysis (taking only the deri-vative terms into account) shows immediately that this scheme is unconditionnally unstable. The amplitude relation over one time step is

-+i+l -+i

U =§.U

where the amplification matrix §is given by (17is the unity matrix) §= [j _ ~b.t sine b.x with eigenvalues g. = 1 - i a. b.t sine J J b.X 2

I

9

·1

J m 2 2

=

1 +

la·1

(~)

J m b.X

and the stability relation (15) will never be fulfilled. 2.5 Lax scheme

The simple Euler scheme can be stabilized with a spatial smoothing applied on the first term on the right hand side of

e q u a t ion (117) . conservative form = b.t (-+ft -+fi ) b.t (-+h i -+hi ) (18 ) k+l- k-l + k+l+ k-l a 2b.x 2

(18 b)

With the von Neumann analysis one obtains the amplitude relation

-+R-+1

-+R-U =§.U

where the amplification matrix§is this time given by

§= [lcose -

ic:4~

sine

!1x with eigenvalues 9 . = cose

-

i a . !1 t sine J J _!1X 2 2

I

9

·1

= cos2e +

I

a

·1

(~) J m J m !1x 2 2 sin e

and the stability condition (15) is fulfilled for the CFL

condi-tion (10).

The discretized equation (IS) can also be written as follows

and a comparison with equation(17)shows already that the smoothing

• . -+R.+1

-+R-term takes the form a second derlvatlve. At convergence uk _{=u k}

and by introducing the Taylor series expansion for

~~+I'

-+R. -+R.

-+R-uk-I' hk+I' hk- 1 in the above equation (to investigate the

system converges to the solution of the following stationary system

The Lax scheme introduces a second derivative with a coefficient of O(~x) (~x and ~t are of the same order of magnitude because of the CFL condition) and is thus only first order accurate. By analogy with the Navier Stokes equations, the second derivative is a viscous term, the viscosity of which can be associated with a Reynolds number Re

=

O(x) (based on the length X·~x of the nozzle and the velocity

IWj

I

_m).

The influence of this viscosity on the solution is

clearly visible in figure 6 where the shock is completely smoothed. The favourable effect of the viscosity on the convergence speed can be seen in figure 7, where the convergence of the Lax scheme is compared with the convergence of the corrected viscosity scheme.

2.6 Corrected viscosity scheme

A second order accurate scheme can be constructed by periodically correcting the viscosity term in the Lax scheme. conservative form

,11.+1

Û k

quaSi linear form

-( 19a)

where the last term or correction term is updated every Nv

itera-tions (t*

=

0, Nv, 2Nv, 3Nv, ... ). This scheme is also known in

the literature as the damping surface technique and was first pro-posed by McDonald (Ref. 3).

During the Nv time steps with constant correction term the scheme has the properties of the simple Lax scheme, namely stable for the CFL condition and a high convergence speed but a poor accuracy. Every Nv iteration, however, the solution is cor-rected for the viscosity effects to obtain finally at convergence a second order accurate solution.

The stability over Nv time steps can be investigated with a von Neumann analysis. The amplitude relation for one time step is ~t+l U where

!B

=[/ cos a - i

df

~

sin a I1X d

=

a(l-cosa)

A recursion formula over Nv time steps gives the relation ~t+Nv

a

~t

U ='::7° U

with amplification matrix

with eigenvalues Nv Nv b.

-

1 gj

=

bj J d b .

-

1 J

where

tlt

b.

=

COS6 - i a· -- sin6

J J tlx

One can prove that inequality (15) will never be satisfied if 6 can take all values. Both Couston (Ref. 2) and Lehthaus (Ref. 16) calculated numerically the maximum eigenvalue Igj1m for the discrete spectrum of wave numbers that can exist in a discretized domain (due to the discretization only waves with a wavelength which is a multiple of tlx can exist; the Fourier spectrum is then discrete and given by w = k~x or À = 2ktlx, k = 1, 2, ... ,X/2).

They found that if the CFL condition is satisfied, the scheme is stable when X/Nv<1. Smolderen (Ref. 3) came to the same con-clusion applying the corrected viscosity scheme on a simple model equation. Rather than trying to solve the inequality (15) for the above expression for gj' we investigated the stability of the corrected viscosity scheme experimentally. The results of these numeri cal tests are presented in figure 8, showing that values of 2 and more can be used for the ratio X/Nv before instability occurs. This important conclusion is confirmed by the practical results of Couston (Refs. 1, 2).

If too low values of Nv can cause instability, too high values of Nv on the contrary reduce the convergence speed. The following table shows an example of the influence of Nv on the total number of iterations N to reach convergence (for com-partson: the Lax scheme converges in 80 iterations).

NOZZLE 1, IC1, X = 31, CFL = 1, MRS = 1.48, VC = 0.98

Nv 10 14 17 20 25 30 40 50 60 80

N unstable 520 460 680 675 960 1240 1500 1740 2320 Clearly the value of Nv should be chosen close to the stability

limit. But the choice is not very critical as within a certain range of Nv values the convergence speed does not change very much This confirms the theoretical investigations of Smolderen (Ref. 3) and Essers (Ref. 7) on the subject of convergence speed.

The coefficient a in (19) permits to retain a certain

amount of viscosity to obtain a correct shock capturing. At

con-9,+1 R, R,*

vergence uk

=

uk

=

uk and like for the Lax scheme we find that (19) converges to the solution of the following stationary system

+ + AX2 ~2u+

c:/f~

=

h + (I-a) Ll 0 + 0(r:,x2 )

ax

2r:,t

ax

2

(20)

If I-a

=

O(~x) the scheme is second order accurate and the

remaining viscosity corresponds with a Re

=

0

(~).

In practical

l-a

application a is a function of the density gradient :

a

=

VC

[1 -

~ IP::1-2;:+p~:111

where VC is a given viscosity coefficient. In this way a is more important in the shock region and affects less the accuracy

further away. Optimum values of VC can be given in function of the expected shock strength (MRS

=

Mach number up~tream /Mach number downstream of the shock).

MRS VC shock free solution 0,99 1,50 1,80 0,98 0,96 1,95 2,10 2,40 0,95 0,94 0,93

Figures 9 and 10 finally show two results which permit to judge the accuracy one can obtain.

2.7 Lax Wendroff scheme

A classical, second order accurate, explicit discreti-zation scheme for time marching calculations is the Lax Wendroff scheme. There exist several versions in one or two steps of this scheme, but a very interesting and frequently used form is the two steps scheme proposed by Richtmyer.

= 1

2

quasi linear farm

(21a.l)

(21a.2)

(21b.1)

(21b.2)

This scheme is in fact a sequence of a Lax discretiza-tion (18) ahd a midpoint leapfrog discretizadiscretiza-tion (which is a centered time centered space scheme). The Lax scheme is stable for the CFL condition, is highly dissipative and converges very fast but is only first order accurate. The midpoint leapfrog

scheme (Ref. 11) is second order accurate but is only marginally stable if the CFL condition is fulfilled and contains na dissi-pation at all. A van Neumann analysis far bath steps tagether gives where with eigenvalues 2 gj = 1 - 2a. (~t) J ~x 2

I

9 _J

·1

= 1 -

41

a j

l

2 2

sin2_{e -} 2ic7t~t _sinecase

~x

sin2_e

-

_2ia·

!!

_sinecase

J ~x

2

[1

_,a.,2(~t)21

(!!)

sin4_e

showing that the two steps scheme is also stable for the CFL condition.

As the Lax Wendroff scheme is derived from a second order accurate Taylor series expansion in time, the dissipation introduced by this scheme is very small [corresponds with a fourth order derivative with a coefficient of 0(~x3)]. This limited dissipation can ensure a stable and convergent calcu-lation but not a correct shock capturing. Additionally, the scheme is only marginally stable on a sonic line (if one of the ~igenvalues a_j = 0, what is the case when u =

a,

then

Igjl

=

1 and (15) is marginally satisfied), which could give

m

rise to non linear instabilities. For these two reasons an artificial viscosity must be introduced in the scheme. Usually this is done by adding a third step or viscosity step (Refs. 3, 7, 9) of the form

+ t+2

u

k,sm (21.3)

with

I-a

=

O(~x) to maintain a second order accuracy. Like for the corrected viscosity scheme the viscosity coefficient a can

be a function of the density gradient for example to increase the viscosity in the shock region.

Couston (Ref. 3) applied both the Lax Wendroff scheme (21) and the corrected viscosity scheme (19) to calculate the one dimensional flow with shock through the nozzles of figure 5. From his results he concluded that both schemes give the same final accuracy with a slower convergence for the Lax Wendroff scheme. Essers (Ref. 7) also performed calculations with the Lax Wendroff scheme on the one dimensional flow with shock

through a convergent divergent nozzle. Using the same geometry, the same number of discretization points, the same initial con-ditions and the same convergence test, we found that the cor-rected viscosity scheme (with optimum values for Nv and VC)

converged in 820 iterations' (CFL

=

1) against 650 iterations (CFL

=

0.9) for the Lax Wendroff scheme. However, every itera-tion with the two step Lax Wendroff scheme consists of three equations and is therefore more time consuming. Our experiences confirm that both schemes give the same final accuracy with a faster convergence but a higher calculation time for the Lax Wendroff scheme. The corrected viscosity scheme is also easier to program. Moreover, it is less problematic to obtain a good shock capturing. Accurate shock calculations are easily obtained by using a correct value for the viscosity coefficient in the corrected viscosity scheme, whereas great care must be taken wh en introducing artificial viscosity in the Lax Wendroff scheme to avoid oscillations around the shock.

One has to keep in mind, however, that the corrected viscosity scheme only reaches second order accuracy at convergence whereas the Lax Wendroff scheme is also second order accurate in time and can therefore, unlike the corrected viscosity scheme, calculate correctly truly instationary flows.

2.8 Corrected damping method

An attempt was made to increase the convergence speed by introducing a damping in the partial differential system as proposed by Smolderen and Essers (Ref. 7). We already indicated earl ier that the Euler equations do not attenuate the perturbation waves. Howeyer, a non derivative term of the form ku on the left hand side of one of the systems (3) to (6) introduces adamping of the perturbation waves if k >

o.

Th avoid that the calculation converges to the solution of an erroneous stationary system, the non derivative term is periodically corrected in the same way as the viscosity term in the corrected viscosity scheme. The partial differential system to be solved becomes, in

conservative form

+ +f +

au + ~ + k(u-u*)

=

h (k

>

0) (22a)

quasi linear form

aû

-flau

+ +* +

+C7T-- + k(u-u )

=

h (k > 0) (22b)

at

ax

+

wherec/fand f take the form given by one of the equations (3) to ( 6 ) .

As the artificial term is a non derivative term, it does neither change the type of the system nor the characteristic direc-tions and therefore the physical boundary conditions and the CFL stability condition remain the same. The system (22) c~ now be

+

discretized by any scheme and, as for the non derivative term h, the damping term will not introduce instability provided it is approximated by a mean value of the surrounding points.

Unfortunately, the convergence speed is very sensitive to the choice of the two parameters k and Nd (number of iterations without correcting the damping term) as is shown theoretically

in reference 7. Applying the corrected viscosity scheme (19) on the system (22) gives the following results for the number of iterations N to reach convergence:

NOZZLE 1, IC2, X

=

46, CFL

=

1, Nv = 15 MRS = 1.58, VC = 0.97

~

0 0.4 0.8 1.2 1.6 40 50 871 796 736 691 721 60 70 1.8 2 2.2 2.4 706 721 721 646 646 676 691 661 661 661 691

MRS

=

1.79, VC

=

0.96

-A

0 0.2 0.4 0.5 0.6 0.8 1.2 1.6 2 30 1126 40 1126 50 1096 1111 1051 1051 1081 1096 1081 1126 1216 60 1066 70 1051 80 11066 MRS

=

2.0, VC

=

0.95

A

0 0.4 0.6 .· 0.7 0.8 0.85 . 0.9 ₁ 1.2 1.6 2 2.5 20 976 40 1021 50 991 S91 1021 1036 976 991 991 1021 1051 1126 1186 1246 60 1006 70 1036 MRS

=

2.45, VC

=

0.93

A

0 0.2 0.3 0.4 0.6 20 796 841 826 871 916 50 796 841 871 916 80 946

quite laboriously, but also that the optimum values vary from case to case, giving ca 30% gain at the best and no gain at the worst. Clearly the corrected damping method is not suited for use wtth the corrected viscosity method and further investiga-tion was discontinued.

2.9 Influence of the initial conditions

Intuitively one would expect that a better initial

guess of the flow field automatically gives a faster convergence. However, increasing the accuracy of the initial conditions has no proportional effect on the number of iterations to reach convergence.

The convergence process can roughly be devided in two parts. The first part corresponds to the elimination of the large initial errors and depends strongly on the initial conditions. But when the flow is more or less settled, the convergence is determined by the dissipative effect of the scheme and this dissipation must remain small to achieve a good accuracy. This second part is characterized by an exponential decay of the

error and deRends on the scheme but not on the initial conditions. The influence of the initial conditions was tested

numerically with three different sets of initial data, repre-sented in figure 11. IC1 and IC2 were already described in § 2.3,

IC3 are uniform conditions with a chosen Mach number . . Figure 12 gives for each IC the evolution of the convergence error. From these (and other) results we can conclude that it is not

worth-while spending computer time for the calculation of a sophisticated first guess but that it is more important to have smooth initial conditions.

2.10 Variàble time step

The Euler equations (1) are linear equations with variable coefficients (quasi linear equations) and in this case both the CFL stability condition and the von Neumann stability

analysis give alocal condition, to be fulfilled in every point. Therefore at every iteration the most restrictive time step over the whole spatial domain must be selecte4. But this means that in only 1 point the time step equals the maximum permissible one and in all other points the time step is of ten considerably

smaller than allowed. In these points the stability of the

scheme is underestimated, the perturbation waves do not propagate as far as they might and the convergence is slowed down.

The problem can be avoided by solving the system -+ -+

au af -+

I

a

·1

- +

=

h

J m _{at ax} (23a)

or in quasi 1 i nea r form

-+ -+ _-+

I

a

·

1

au +c:;If~

=

h

J m _{at ax} (23b)

which has a different transient evolution but the same statio-nary solution as the systems (3) to (6). The system (23b) is equivalent with "-+u -+-+ _0 _ + d(~

=

hl at ax where -+ 1-+ hl

=

•

h

la ·1

J m

and the CFL condition reduces to

I

as the maximum eigenvalue of

c:;If

is 1. The time step is now independent of x (if the grid is uniform) and is thus equal to the maximum permissible one in every point. In practice one uses

the condition (10) to calculate a time step that is function of x. This variable ~t is introduced in the schemes described in the previous chapters, leading to an equivalent result. Figure 13 presents a comparison of the convergence error for two cal-culations using system (5) and (23a), showing a gain of ca. 20% in number of iterations. The extension for a non uniform grid is straightforward.

2.11 Mesh refinement

The total time T which a time marching calculation

uses to converge is determined by the time the perturbation waves need to die out or to convect out of the domain. This

time depends on the partial differential system,on the boundary conditions and on the discretization scheme, but not on the time step (if the stability condition is fulfilled). The number of iterations N to reach convergence is therefore inversely proportional to the time step (T

=

N~t) and using the CFL con-dition also inversely proportional to ~x or directly proportional to the number of points X. The calculation time per iteration is also directly proportional with X so that the total computer time changes quadratically with the number of points X.

Relations like (20) show that to increase the accuracy of the final result, the discretiz~tion step must be reduced leading to aquadratic increase of the computer time (cubic for a two dimensional and with the 4th power for a th ree dimensional problem).

This dependence can be relaxed by starting the calcu-lation with a coarse grid (one usually takes half of the dis-cretization points) and change then to the fine grid to increase the accuracy. Qne can in fact consider the results of the coarse grid as improved initial conditions for the calculation with the fine grid. In principle, this approach could be used several

time during one calculation, but unless a very high final accuracy is needed (requiring a very dense grid), it is not worthwhile

This method was used with considerable success by Couston (Ref. 2), although he somewhat overestimated the gain in his theoretical analysis. Figure 14 gives the convergence error for two calculations performed on the same nozzle under the same conditions, but one with twice the number of discre-tization points than the other showing clearly the proportion-ality between N and ~x. Figure 15 gives the convergence error when changing from coarse to fine grid at 800, 1300 and 1800

iterations with the coarse grid, corresponding with a not completely converged, just converged and fully converged

solution. Clearly the mesh refinement introduces perturbation waves which need a certain number of iterations to disappear. It is of no advantage to wait for complete convergence of the solution with the coarse grid. The gain is not so much due to a gain in number of iterations but rather to the fact that approximately half of the iterations are performed with only half the computational time per iteration. Figure 16 gives the result obtained with the fine grid. Comparison with the result of figure 9 shows a better accuracy for the shock.

REFERENCES

1. COUSTON, M.: Time marching finite area method.

VKI LS 84 "Transonic Flows in Axial Turbomachinery", February 1976.

2. COUSTON, M.: Méthode de calcul de l'écoulement inter-aube pseudo-tridimensionnel en régime transsonique. Doctoral Thesis, Université Libre de Bruxelles, September 1976.

3. COUSTON, M.; McDONALD, P.W.; SMOLDEREN, J.J.: The damping surface technique for time dependent solutions to fluid dynamic problems.

VKI TN 109, March 1975.

4. DENTON, J.D.: Extension of the finite area time marching method to three dimensions.

V KIL S 84 11 Tra n s 0 n i c F 1 0 w sin A x i alT u rb 0 m ach i ne r y 11 ,

February 1976.

5. McDONALD, P.W.: The computation of transonic flow through two dimensional gas turbine cascades.

ASME P 71 GT 89, March 1971.

6. ESSERS, J-A.: Time dependent methods for mixed and hybrid steady flows.

VKI LS 1978-4 "Computational Fluid Dynamics", March 1978.

7. ESSERS, J-A.: Méthodes nouvelles pour le calcul numériQue

d'écou-lements stationnaires dé 'fluides ~arfaits compressibles.

Doctoral Thesis, Université de Liège, 1977.

8. GOPALAKRISHNAN, S. & BOZZOLA, R.: Numerical representation of

inlet and exit boudnary conditions in transient cascade flows.

ASME P 73 GT 55, Apri 1 1973.

9. GOPALAKRISHNAN, S.: Fundamentals of time-marching methods. Application of the time dependent technique to turbomachinery cascades.

VKI LS 59 "Transoni c Flows in Turbomachi nery", May 1973. 10. RICHARDSON, S.M.: Stability theory of time-marching methods with

extensions to the theory of ill-conditioning of steady state methods.

ARC CP 1383, September 1976.

11. ROACHE, P.J.: Computational fluid dynamics. Hermosa Publ., Albuquerque, 1976.

12. SMITH, G.O.: Numerical solution of partial differential equations: finite difference methods.

13. SMOLDEREN, J.J.: Partial differential equations. VKI CN 102.

14. LEHTHAUS, F.: Anwendung eines Zeit-Schritt-Verfahrens zur Berechnung der transsonischen Durchströmung ebener Turbinengitter und experimentelle ÜberprUfung. DFVLR-AVA Bericht 251 77 A 01, February 1977.

15. SIEVERDING, C.: The turbine blade definition. Experimental data on two transonic turbine blade sections and comparison with various theoretical methods.

VKI LS 59 "Transonic Flows in Turbomachinery", May 1973. 16. VEUILLOT, J.P.

&

VIVIAND, H.: Méthodes pseudo-instationnaires

pour le calcul d'écoulements transsoniques. ONERA Publ. 1978~4, 1978.

)(=0 INLET CONDITIONS t _INFLU -OF THE POINT

---

---x

=

1 OUTLET CONDITIONS o~---~-.-t = 0 INITIAL CONDITION

,

x QJbJc: CHARACTERISTICS IN P.

t /j,x

a

.

lNLET t SLOPE (M<l) 1. U+a 2. U 3. U-a /j,x b.OUTLET

FIG. 2 - CHARACTERISTIC DIRECTIONS IN A BOUNDARY POINT.

x

/ C F

,

"

x

,

"

_"

,

_,

,

_,

Q,b,c CHARACTERISTIC DIRECTIONS

FIG. 3 - C.F.L. STABILITY CONDITION

1 f (M) 1

"V

0.75 ~ 0.5 L - -_ _ ...J.L--_ _ -&' _ _ _ ....L' _ _ _ " - - - I _ _ ...J

o

1 2 3 4 M 5 FIG. 4- CORRECTION FACTOR ON THE C.F.L. NUMBER

FOR CONSTANT TOTAL ENTHALPY.

1.6 A Ath 1.4 1.2 1

!---10 NOZZLE 1 i

I

I

/

~

/

----

19 28 37

x

(NUMBER OF POINTS IN AXIAL DIRECTION )

/

46 1.6~---~---~---~---~---~ A Ath 1.4 ~---+---+- ---NOZZLE 2 1.2~----~~---~---~---+---~ 1 10 19 28 37

x

(NUMBER OF POINTS IN AXIAL DIRECTION )

FIG. 5 - NOZZLE GEQMETRY

1.4 1.2 1 .

.8

. 6 .4 MACH

v

/

/

/

••••••••

_{• •}

••

••

_•

_•

••

e e_

"

4~ • " e. e

_...

~

1 10 19 CORRECT SOLUTION • CALCULATED SOLUTION NOZZLE 1 INITIAL CONDITION 1 MRS =2. CF L =1. 28 37

_x

(NUMBER OF POINTS IN AXIAL DIRECTION )

FIG. 6 - CALCULATION RESULT WITH THE LAX SCHEME.

,.

~---~---~---~---~ 3. H-~---~---~---+_---~ 2.~---1~--+---~---+_---~ 1. H+---~--~----~---~---+_---~

O.

500 1000 1500 N 2000 (NUMBER OF ITERATIONS ) NOZZLE

,

INITIAL CONDITIONS 1 P2

=

0.800 CF L

=

1.0 LAX SCHEME

2 CORRECTED VISCOSITY SCHEME VC

=

0.9'0 Nv= 25

FIG. 7 - CONVERGENCE OF THE LAX SCHEME AND THE CORRECTED VISCOSITY SCHEME.

C.F.L. 1.5 1. .5

o.

.1 M.R.S. Cl) SHOCK FREE x 1.59

•

SHOCK FREE

•

2.0 STABLE 1. NOZZLE 1 INITIAL CONDITIONS P2 VC 0.868 1. 0.850 0.97 0.790 0.99 0.800 0.9' 1 X Nv X -46 '6 31 '6 (*) INITlAL CONDITIONS

=

FINAL SOLUTION

10.

6.t

CONSTANT (*) CONSTANT CONSTANT VARIABLE

1.6 MACH 1.L. 1.2 1 . . 8 .---r---~---_1--r_----_1---~ .6 ~---__+---_t_----_+_---.. - - - -- - I - -- - - ---i 1 10 19

CORR EeT SOLU TlON

•

CALCULATED SOLUTION NOZZLE 1 INITIAL CONDITIONS 1 28 37

_x

( NUMBER OF POINTS IN AXIAL DIRECTION ) CFL = 1. Nv = 25 VC

₌

0.9L.0 MRS ₌ 2.0 N = 1350

FIG. 9 - CALCULATION RESULT WITH THE CORRECTED VISCOSITY SCHEME.

1.6 MACH 1.4 r - - - t - - - 1 r - - - t - - - t - - - j 1.2 1 . . 8 ~---r---+_---1---~~---~---~ .6 ~---~~~----~~---r---~---~ .4 ~---r---1r---r---t---1 1 10 19 CORRECT SOLUTION

•

CALCULATED SOLUTION NOZZLE 2 INITIAL CONDITIONS 1

37

28

x

(NUMBER OF POINTS IN AXIAL DIRECTION ) CFL

₌

1 . Nv

₌

25 VC

₌

0.960 M RS

₌

1.69 N

₌

1300

FIG. 10 - CALCULATION RESULT WITH THE CORRECTED VISCOSJTY SCHEME

1.6 MACH 1.2 1 .

.8

.6

.

,

/-1 , / I , / , / I / ' , / I , / / ' _I , / / ' I , / / ' I , / / ' I / ' / ' _I . / / '

~

I

V

I

~

I I I

~

I

~

IC3 _r--. ...

~

... 1 10 19 ---CORRECT SOLUTION 28 37 X (NUMBER OF POINTS IN AXIAL DIRECTION )

FIG. 11 - INFLUENCE OF INITIAL CONDITIONS : INITIAL DATA

500 1000 1500 N 2000 (NUMBER OF ITERATIONS)

1. INITIAL CONDITIONS 1 (EXACT SUBSONIC SOWTION)

2 INITlAL CONDITIONS 2 (PARABOLIC DISTRIBUTION OF PRESSURE )

3 INITIAL CONDITIONS 3 (CONSTANT DISTRIBUTION )

FIG. 12 - INFLUENCE OF INITlAL CONDITIONS: COMPARISON OF CONVERGENCE ERROR.

3 . _

-I

I

, I

I 2 _I 1

o

500 1000 1500 2000 (NUMBER OF ITERATIONS) 1 CONSTANT /:,. t 2 VARIABLE /:,. t

3. 2. 1.

o.

1000 2000 3000 N (NUMBER OF ITERATIONS) 1. X

=

'6

2. X

=

91 ( X

=

NUMBER OF POINTS IN AX lAL DIRECTION )

FIG. l' - INFLUENCE OF REDUCING THE GRID SPACING ON THE CONVERGENCE ERROR.

E (%) 3. 2. 1.

O.

1000 2000 3000 N 1.000 (NUMBER OF ITERATIONS)

1. 800 ITERATIONS BEFORE MESH REFINEMENT

2. 1300 ITERATIONS BEFORE MESH REFINEMENT

3. 1800 ITERATIONS BEFORE MESH REFINEMENT

FIG. 15 - INFLUENCE OF THE MESH REFINEMENT ON THE

MACH 1.6 , - - - , - - - , - - - , - - - r - - - , 1.4 t - - - _ _ +- - - -- + - - - -... > - t - - - + - - - i 1.2 I - - - +- - -- _ _ . . . . ' - - - + - - - + - - - + - - - i

1 .

. 8~----~---__+---_r_+---~----~ .6 I---+---+---+---~~~---i .4 ~---4---_+---~---~---~ 1 19

37

CORRECT SOLUTION

•

CALCULATED SOLUTION NOZZLE

,

INITlAL CONDITIONS 1 55

73

x

(NUMBER OF POINTS IN AXIAL DIR ECTION )

C F L

₌

1 .

Nv

₌

50

VC

₌

0.9'0

MRS

₌

2.0

N

₌

2650

FIG. 16 - INFLUENCE OF THE GRID SPACING ON

THE FINAL ACCURACY.