NASA ContractorReport172306

1C ASE 1NAsS40A;1CR2_872306

r,OT_ B_T; A,r.._IF,RO_.! "n_I3V,OO:.W

CONVERGENCE OF GENERALIZED MUSCL SCHEMES

Stanley Osher

Contract No. NASI-17070

February 1984

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

•:'._o'r<i_",4

National Aeronautics and

Space Administration L',!"_Lr:;' RE_E,r.RC:4,;ErJr_R

Lf.3P, ARy. NASA

Langley Research Center H.,_.:,q-on,V=R_r't._,

Hampton, Virginia 23665

https://ntrs.nasa.gov/search.jsp?R=19840011218 2018-11-11T00:03:29+00:00Z

'>- t, : 5 c.%:; .!;:\.z>I..s -i-. :. - % 7.32 - -- - , i'.. ... : ....... i. ..... l < ,-.-,..-, <-. . . . - . -. . - . . - . . - r, ./..< : .i --: !.>. < + i; i .>i 3-c. .< i r:; .:; i,..'. - ,. .. , 

-2 p:,; : *:y:> -2 !I.*. + 8, . ..-. T lz.!:; \F * <-... .2. <x::- T- : :,i .-. ,*% c- <\:-? .r ,-. :T.? * .-: ; . . . 2.. . , i z.:: . -. - .- . = -- :.. !"..A. 8 ".., . - i .F : - : z.7. .< ._. ! . -..-?..-I! ... ! ,z 

ac.Ts. hi&~-~-f~-$ 7.-:7:-.;: ,-% . . . - - - - ----- . . -. - - . - - i...! i 3: . ! 8:: !.L! . ..<! ... A < c.:.vj<> 

xChCC-C.:.~-?::l 2.i~ ^; .->c. <7'7.?p<= i-.T!.-<i '27 i:.-; ! .:: !.i = ,;<> . ! ,= >-.-. ~h!,~f. -.',-~;:>7;> hi~!-y -?7!:! h.:.*~ :nrt=-s~~.-;:::?:~*~.~: " I a-S .2 .-.!..!... !.'!!.-.?. ?.:.-*<..-> !.!!..-<A -I..' !I.-.! i'?: ; i:, 

-, > .-,!!! :'L-L.- " .*'.<-.?.-\ !.-*.--. -- -. - '2 ;:.,- .> -" ;'.. - gi$2c:Z<?: ..-,... ,;.. . .. , . ; , ., :> . .? .-. s. ..-. z. $ &!-.F 'TI: .=.-.-' : i:.::~; 5~=.:3. ; j ~r~ p,:-;rc e~ra -- - . ..... .. , . ,... .' . >_.; + ! - -._. ~:.-;:.-.i_ii:jir~; ; i ;‘ .:. -.k..:- .. - Tic== cj.-..-:: -~;..+-..-..la+ ..-,! !,-. <-.%->!~~.j~'~~~:i: :.-;! >..-! !..-! .> i 1 LC<-? y:;;:~i.-.~ z:.-.i r:~i.~~ ;-.>; . : ! i,r: i i....y!-.<->i ;- 

p.: ;, = 5 ., $7:~ ~22 .?. ti- -. .. . 7,';. ! .. i. =/ i-;.>; ]i-i.%. j .yi; .7..?..?.= p:!<<--:?. < &:- ..:= -T..{ -: !+-.<.: '$ 

: r:; :..c.-..>i : it.-;; <->ii: .: = 

U-:??.>-. f :-:-.: - <?j 2 .; . 1 ,-<, ri,;,fATI ?-iTTC I::... ."""- !i? (in': ,t.Ac C-.J;>+ f !:!!>'-=!..! I.-. 

.-::.--- . -. .. . -. ,.!! . ! !.-. r::.-,.:,; :.ii r5:.-; ;".-'iJS: ::av. /<+~~"R$->a.~~"'L"fi"\ '+'r:>w""f-zT?p:": ; C.:!ip.~~<~CfiA\i:!~~.$T~~.%y~:eE~~..~T~.<->r:'s.;;<<+~~~:-~~ . . .... .. . i...:.--; .... zz$=;-':\ ~LT:PL?.! ! !! 1 !....-?fi.ii'iF.! i<->y.i; : ,!.>>i-c..::zi i=.-;i.,j .-F.ij<.-*, ..-., -,!,!.< ,-. 8 ".. .; ; ;-fij F.i-j.T ! ! 5 :.:i >;)i-;j L-5 y< ?:e ... ::.t:- 2 ri !, .l.. i/ $!-GO), I TI-;?:;! *Fi-ii-j i 3 t;/NAM si/ I NEQ-!G(: .,. !,- 1 - !T + CEI. . -.-. , .-;:.-._.nini. 

.:. CC <'i': 5E.C - - :.\ . hi . " ". :-.. ..;. -. -- !!.-.!!. r!:,-?i-i?.2i . . ~RII.= 12.~.:fi5-x-:i~~--~-~.~~. ~s~~-:~~2:~--.-=r..~..y. .-.:+ ;~...*Eie~;c~Fic~ :-.:.-!-;C.F ~:..!?.=..~..:-i.r..+-.. -.i <-.r-.~!i:~-..:-..*.:*r- ! !.-8 .-. .>..-!!! ! . . . . " .. . - .-.- >..-!!..-! .-5: ~~%~i-~;-;ii~, ;?! ..-> ! ..-*..- ?. ;-Ti .ys :.>i : i-: ;><-;:.-+:.-!: i\-; .+ .? 

=:,-.;-..:-..:<:.-.. I< r ...- ..st< +-. >= *?..& ?:A: it=<-.! <- ..?. + .....-*..-. -.:s .*., ;- .:.->F:c * +*! !:.-. * .-....- $ c- >;=- < +-.:ys ,+% < At-.: -:?... .-..; !:xi:ic! r..i;:.;;i;z! . - . . . , -.. - 2 .-!!..- yi:.>.-s,-.x- c,=.-.: icciic! .2.{ ; -...-,! !.> ,-: .-a.,L .-T<L! .> ,-,d... ,- ;.:.:= ?!.! ;. i-:. ,. ,:.i.;y 

;I:. ,id? <?ass ;is ~rs,j%n to 55 SZ -: ~~;f'~-~:~. appr~im~tj~~ t~ f$-.. .-!:c -7.. :._.<-;! -..-. M>- c:.:;- - .-...-.+ 

:-. r..:-..z. > :-. L= .-A .-. =-. + :-. k... :-., y;; :-., r. < :-. 7*-: -= ..-i .-. . <. ? :;.: :: .. - . -. ,, !.! <-! >.,-! ! .- :-,;-.i - 1-;;t-.-. . - . . . - - ! z, !;!,>:.>? = 

3 117 62

CONVERGENCEOF GENERALIZEDM_SCL SCHEMES

Stanley Osher

University of California

Los Angeles

Abstract

Semi-discrete generalizations of the second order extension of Godunov's

scheme, known as the MUSCL scheme, are constructed, starting with any three

point "E" scheme. They are used to approximate scalar conservation laws in

one space dimension. For convex conservation laws, each member of a wide

class is proven to be a convergent approximation to the correct physical

solution. Comparison with another class of high resolution convergent schemes

is made.

Research supported in part by NSF Grant No. MCS 82-00788, ARO Grant No.

DAAG 29-82-k-0090, and NASA Grant No. NAG-I-270. Part of the research was

carried out while the author was a visitor at the Institute for Computer

Applications in Science and Engineering (ICASE), NASA Langely Research Center,

Hampton, VA., which is operated under NASA Contract No. NASI-17070.

i

I. Introduction and Preliminaries. Recently, there has been an

enormous amount of activity related to the construction and analysis of

"high resolution" schemes approximating hyperbolic systems of conservation laws. Some examples of the successful consequences of this

activity can be found in the proceedings of the latest (sixth) AIAA Computational Fluid Dynamics Conference, [I], [IO], [18]. Extensive

bibliographies can also be found in these papers.

Our aim here is merely to construct and prove the convergence of a

subclass of these schemes approximating scalar convex conservation laws.

This subclass is based on an idea of van Leer [17]. He christened his

algorithms MUSCL (Monotonic Upstream Centered Schemes for Conservation

Laws) schemes, and with mixed emotions, we shall use his acronym here.
In future work with S. Chakravarthy, we shall extend this construction to systems in multi-dimenslons, using trlangle-based algorithms.

This work in progress will stress the computational aspects of the algorithms, especially as they relate to the Euler equations of compressible

gas dynamics. In earlier work with the same author, [12], we constructed, and proved convergence of, a class of high resolution schemes

approximating scalar convex conservation laws. We also showed for certain high resolution approximations to systems, that limit solutions

satisfy an entropy inequality.

We shall consider numerical aproximatlons to the initial value

problem for a single conservation law in one space dimension

(1.1)(a) wt + f(w)x : O, t > O, -I < x < I,

with a periodic boundary condition:

(b) w(x + 1,t) _-w(x,t)

and initial condition:

(o) w(x,O) : Wo(X).

It is well-known that solutions of (1.1) may develop discontinuities in finite time, even when the initial data are smooth. Because of

this, we seek a weak solution of (1.1), i.e. a bounded measurable func-

-2-
tion w, such that, for all g <- Co(R x R+),

(1.2)(a) _ (w_t + f(w)_x) dx dt = 0

RxR+

(b) t-->O lim _lw(x,t) - Wo(X)11L1 : O.

Solutions of (1.2) are not necessarily unique. For physical reasons, the limit solution of the viscous equation, as viscosity tends to

zero, is sought. In the scalar case, this solution must satisfy, for

co

_ CO (R x R+), 6 [ O, and all real constants c:

(1.3)(a) -_ (lw - cl_t + sgn(w - c)(f(w) - f(c))dx) dx dt _ 0

This is equivalent to the statement:

(1.3)(b) _-{ lw -cl + ((f(w) - f(c)) sgn(w - c)) < 0,

in the sense of distributions.

Such solutions are called entropy solutions. Kruzkov has shown in

[9], that two entropy solutions satisfy

(1.4) llw(x'tl) - u(x'tl)llL1 ! llw(x'to) - u(x'to)llL1

for all tI _ to. Hence, (1.3) guaranteesthe uniquenessof solutions

to the scalarversion of (1.2). Existencewas also obtainedthere.

-3-
Equations (1.3) may be viewed as the statement

(1.5) _-_V(w) + F(w) < 0

for each of the convex "entropy" functions of solutions w to (1.1):

V(w): lw- c[

and their associated flux functions

F(w) = (f(w) - f(c)) sgn(w - c).

For scalar convex conservation laws whose solutions lie in the space

BV, a single entropy inequality (1.5) for any strictly convex V(w)

and its associated entropy flux F(w), satisfying

(1.6) F'(w) : V'(w)f'(w)

is enough to imply existence and uniqueness of solutions to (1.1). This

follows from the results of DiPerna [3].

Next we consider a semi-discrete, method of lines, approximation to

(1.1). We break the interval (-I,1) into subintervals:

lj: {xl(j- I/2)A_Ax < (j.I/2)_}

j = O, ±I, ..., ±N, with (2N + I)/\: 2.

-4-
Let xj : JA, be the center of each interval Ij, with end

points xj_1/2,x.3+1/2.

Define the step function for each t > 0, as

U/_(x,t)= uj(t),

for x _ Ij.

The initial data is discretlzed via the averaging operator T/__,

I

___wO(x)= _ _Ij w0(s) ds = uj(0) for x _ Ij.

For any step function, we define the difference operators

/___±uj: i(uj± 1 - uj)

D+_ uj : 1__±uj /\

A method of lines, conservation form, discretization of (1.1), is a system of differential equations

(1.8) _-_uj + D+hj_i/2 = 0, j = 0, ±I, ..., ±N

U/._(x,O) : _.._Wo(X), for x + Ij.

Here, the numerical flux defined by:

-5-
(1.9) hj_i/2: h(uj+k_l,...,uj_k).

for k > I, is a Lipschitzcontinuousfunctionof its arguments,satlslies the consistencycondition:

h(w,w,..,w) : f(w).

It is well known that bounded a.e. limits,as /__k-_0, of approximate solutionsconvergeto weak solutionsof (1.1), i.e. (1.2)(a)is

satisfied. However,this does not also imply that the limit solutions

will satisfyany of the entropyconditions(1.5),let alone the general

condition(1.3). Some restrictionson h are required.

A simpleclass of flux functions h, for which (1.8)converges,

for all f, to the unique entropysolutionin L°°(LI(R);[0,T]),as

/i-_ 0, for any T > O, is the class of "E" schemesintroducedin

[11]. Such schemes satisfythe following:

A consistentschemewhose numericalflux satisfies

(1.10) sgn(uj - uj_1)[hj_i/2- f(u)] ! 0

for all u between uj_I and uj is said to be an E scheme. This is

at present, the most general class of schemes known to converge in the

nonconvex case.

It is clear that this class includes the widely known class of

-6-
three point monotone schemes, i.e. those for which

hi_I/2= h(uj,uj_1),

with hj_i/2 nonincreasing in its first argument, nondecreasing in its

second.

We denote partial derivatives of a numerical flux via

_-_j+y h(uj+k,...,uj_k+I) : hy.

Thus a three point scheme with a differentiable flux function is monotone iff

h1< 0 < ho.

Any numerical flux can be written

(1.11) hj.I/2: I/2(f(uj) + f(uj.1)) - I/2Qj.I/2(A.uj)

where Qj+I/2 can be viewed as the viscosity of the scheme [15].

One particular three point scheme is due to Godunov, [5], and has a

special significance in this theory. The flux for Godunov's scalar

scheme can be defined by

(1.12) hG_i/2= hG(Uj,Uj_I) : mln f(u), if uj_I < uj

uj_1<_u<_uj

-7-
= max f(u), if uj_ I > uj.

u.3- 1>u>u. - - 3

One can thus, [11], characterize E schemes as precisely those for

which:

G

(1.13)(a) hj_i/2_ hi_I/2, if uj < uj_I

G

(b) hj_i/2_ hj_i/2, if uj > uj_1,

or, as Tadmor [15] pointed out, those which have at least as much

viscosity as Godunov's scheme

G

(1.14) Qj_I/2_ Qj_I/2.

It follows from [11], Lemma (2.1), that these approximations are,

at most, first order accurate.

Together with an entropy inequality, a key estimate involved in

many convergence proofs, is a bound on the variation. For any fixed

t _ 0, the x variation of UA(x,t) is

B(U/k)= _- I/k+uj(t)_. J

If we can write

(1.15)(a) _ +hi_l/2 : -Cj+I/2 /k +uj + Dj_I/2 /k uj

-8-
(b) cj+I/Z2o

(c) Dj_I/2>_ 0

then it is easy to show, [12], using an argument of [19], that, for

t1_> t2 >_0

(I.16) B(U/__(",tI)) __B(U/__(",t2)).

Marten in [7], pointed out for explicit methods, that this decomposition could be obtained for schemes which are higher order accurate. In

our present method of lines context, it involves a five point, consistent, approximation:

_t I (1.17)(a) : -_ /._..hj_i/:2Cj.I/2A +uj - Dj_I/2A_u j

with

(b) Cj+I/2= C(uj+2,uj+I,Uj,Uj_I) _>0

(c) Dj_I/2: D(uj+I,Uj,Uj_I,Uj_2) _>0,

both Lipschitz continuous functions of their arguments. (See also van

Leer [16].)

In addition to (1.16) we have a maximum principle for (1.17),

[12],:

--9--
(1.18) min Uk(O) < uj(t) < max Uk(O) k k

for each j and all t > O.

Morover, in [12] we also showed a limit on the possible accuracy of

approximations of the type (I.17). Any scheme of the type (I.17) is at

most first order accurate at nonsonic critical points of u. (A sonic

point u is one such that f'(u) = 0). Thus, although schemes of this

type can be made to be as high as third order accurate, Lipschitz continuity implies a local degeneracy to first order accuracy at smooth

maxima and minima. This local degeneracy, together with some results on

initial boundary value problems in [6], indicate strongly that overall

second order accuracy is the best possible.

Following Harten [7], we call algorithms of the type (1.17), total

variation diminishing, or TVD schemes.

It is known, that although (1.17) generates a compact family of

solutions, UA(x,t), in L°°(LI(R),[0,T]), (if w0 <-LI _ L°° N BV) ),

limit solutions may not satisfy any of the entropy conditions (1.15),

and hence may not even be unique (1.12) See e.g. [4], [1i], and [15].

One way of avoiding this difficulty is to use an inequality obtained in

[11].

Let V(w) be any convex function. We showed in [11], Section III,

that, for any solution o£ any scheme (1.8):

-i0-
(1.19) /_(_-_ d V(uj) + D+FA(Uj)) = _ujuj+1 dw V"(w)[hj+i/2_ f(w)],

where the approximate entropy flux is defined through

(1.20) FA(Uj) = F(uj) + V'(uj)[hj_i/2- f(uj)].

Thus, a sufficient condition that any limit solution satisfy (1.5)

for a fixed convex V, is that

_uj+1_

(1.21) _u. aw V"(w)[hj+i/^-'_f(w)] < 0 3

In order that the above inequality be valid for all convex V, it

is necessary and sufficient that hj+I/2 correspond to an E scheme,

which implies again that the approximation be at most first order accurate. Thus, for the schemes to be constructed in the following section

we shall only obtain our entropy inequality (1.21) for a single V(w),

say V(w) = I/2w2. The following Theorem summarizes the technical

hypotheses needed for convergence.

THEOREM 1.1. The sequence of approximate solutions converses _._.

a__sA-_ 0, to the unique solution of the scalar, convex conservation law

(_._) provided that th___elnitialdat____i_ss ai__nnBVandtha a _ inequalities

(1.17) and (1.21) are valid for a sinsle convex V(w).

The proof is analogous to that of Theorem (4.1) in [12].

-ii-
II. Construction of TVD MUSCL Schemes. The MUSCL scheme as discussed e.g., in [2] and [8] is a second order accurate extension of

Godunov's method [5] that is based on ideas expressed by van Leer in

[16].

We assume that uj(t) is known. The first step in MUSCL is the

reconstruction of a piecewise linear description of the solution. In

the interval I. the result of this operation is

3

(2.1) z(x,t) = uj(t) + (x - xj)sj(t), for (x - xj) < A/2.

Here sj(t) are slopes that satisfy

(2.2) sj = Wx(Xj,t) + 0(A)

subject to some monotonicity constraints discussed below.

Godunov's scheme, in our semi-discrete context, is the following.

numerical flux hG(uj+1,uj) is computed by solving the Riemann

The

problem, i.e. the initial value problem (1.1) with initial data

w _-uj(t), for x < xj+I/2

w _ uj+1(t), for x _ xj+i/2.

The resulting unique entropy condition satisfying solution is a

function of the type

-12-
tj+I/2J

A closed form for w was recently obtained in [11].

Then Godunov's flux is defined through

(2.31 hi+I/2= f(w(O))

which is the same as (1.12).

A fully discrete, explicit in time, Godunov scheme has the same

numerical flux, after we impose a CFL restriction which prevents the

interaction of solutions of adjacent Riemann problems The literal

second order accurate extension to this fully discrete approach would be

to compute the solution to the initial value problem (1.11, with initial

data

w _ uj(t) + (x - xj)sj(t), for x < xj+i/2

uj+1(t) + (x - xj)sj+1(t), for x _ xj+i/2.

Find f(w(xj+i/2,s)) for t + X/_ s _ t, With _ > 0 satisfying a CFL

restriction, then compute

(2.41 1_/__ ft+___ "t f(w(xj+I/2,s))ds.

This is the MUSCL numerical flux.

-13-
Unfortunately, obtaining the exact solution to this nonlinear initial value problem with piecewise linear initial data is a non-trivial

business, even in the scalar case. However, at s = t+ , x = xj+i/2, the

solution is the same as for the Riemann problem with initial data:

W _ Uj(t) +_sj(t), for x < Xj+I/2

W w Uj+ I - _sj+1(t), 2 for t h xj+I/2

which is easily calculated, e.g., from the formula in [11]. Thus, the

semi-discrete MUSCL extension of Godunov's algorithm comes through

_t I - hG(uj j+1 j _ "

We generalize this still further as follows. Let

h(uj+1,u j) = hj+i/2 be an arbitrary first order accurate numerical flux

function. Then our generalized MUSCL algorithm is merely:

_u. 1 A A

(2.6) "_t : -mA /--_h(uj+1 --_-sj.1,u _ j +-_'sj).

We have:

LEMMA 2.1. At points w(xj,t) : uj(t) in a neishborhood of which

h(uj+1,uj) i__sC2 with Lipschitz continuous partial derivatives, and

: u ' the al_orithm (2.6) is at least second order

-14-
accurate for smooth functions w.

Proof. Let

(2.7) uj+I/2= I/2(uj+I + uj).

We shall show:

(2.8) /_ - (h(uj+I - _s. 2 3+1'uj + 2_sj) - f(uj+I/2)) = o(/S2)

at these points.

Using consistency of h(uj+1,uj), the left side of (2.8) can be

written as

/\ h I uj+i/2+ Y-_ J I - + dY I - -- 2 _ uj_,uj _

The result is immediate.

Next we make an observation about the viscosity of this scheme,

Qj+l,,.12as defined in (1.11).

LEMMA 2.2. If h is monotone, then the viscosity of algorithm (2.6)

is a deereasin_ function of sj+1

_j an_AC_uj d

-15-
The proof comes directly from (1.11).

sj and sj+1 are both restricted to be

Thus, as expected, if /_ uj /._uj

positive, the most viscous, hence the least accurate, case occurs when

these values are both zero. The scheme then degenerates to first order

accuracy. We can compress or smear ot]tthe solution locally by increasing or decreasing these ratios.

We wish this scheme to be TVD, which will give a natural restrictions on these above ratios.

LEMMA 2.3. If h is a flux correspondin$ to an E scheme, then the

scheme is TVD if

/__sj _ sj+I

(2.9) 0 _</.,_.uj ' _ uj -< I

for each J.

Proof. We have

(2.10) /_/_h(uj - 2_sj,uj_1 + 2_sj_1 )

: /i_ [h(uj+I - ,uj+ _'sj) - f(uj+ _sj)

.r(uj

-16-
sj /_Isj

+ f(uj - -T) - h(uj - T, Uj_1 + _sj_1) ]

Thus, in (1.15) we can write:

(2.11)(a) Cj+I/2= A _ uj) [(uj+1 - /_puj + 2_sj) - f(uj + _sj)

(b) Dj_I/2_ /_ (_ uj) [h(uj - 2 J +_sj) - f(uj - _sj)

._°_-_-._o_- _,u___._.___

Both Cj+I/2,Dj_I/2 are nonnegative because of (1.10) and (2.9).

The restrictions (2.9) can be relaxed somewhat if h is monotone;

we have:

LEMMA 2.4. If h is a flux correspondin5 to a monotone scheme, then

the scheme is TVD if

(2.12)(a) 1 > (s_+I - sj) A

- 2_ uj near points where hI _ 0

-17-
(b) I

--

> (-sJ+1 2A +

u.

sj) _. near points where h0 _ 0

for each j.

Proof. We have

(2.13) _/__ h(uj -_ sj,uj_l + _sj_ I)

: /_ [h(uj+1 ' + 3 - 2 j + s

1 _ )]

+ _ [h(uj - ,Uj + --_-_)- h(uj - 2 'uj-1 + "

Thus, in (1.15), we can write:

I (Gj _2-_)___uj -2_sj+1- sj)] (2.14)(a) Cj+I/2= A _ uj) hl ,uj +

(b) Dj-I/2- /__(_uj)

I

ho (uj - %,_j)[/_._. uj + _2 (sj - sj_1)]

Both Cj_i/2 and Dj_I/2 are nonnegative because of monotonicity

and (2.12).

Remark 2.5. These last two Lemmas can be made "local,"i.e. in

regions of monotonicitywe may relax the restrictionsfrom (2.9)to

(2.12).

-18-
III. Convergence of MUSCL Schemes. In order to prove convergence

for convex f, we need only verify the discrete entropy inequality

(1.21) for the schemes constructed in the previous section (Theorem

(1.1)). We shall do this for the entropy

2

V(w) I, = w

and for schemessatisfyingthe TVD hypothesisof Lemma(2.3),in additionto someotherrestrictionsnearsonlcshockpoints.

Suppose uj < u.j+l , i.e we are at a rarefaction. We then have,

letting uj+i/2 : I/2(uj + uj+ I) as _bove

(3.1) j. j+1 dw[h(uj+ I _ ,uj + ) - f(w)] -uj

/ksj+ I /ks.

= _ uj)(h(uj+1 - 2 , uj +_/) - f(uj+i/2))

+ _uj+1 uj dw[f(uj+i/2)- f(w)]

= [I]+ [II].

Now, using the hypotheses of Lemma (2.3), we have:

(3.2) [I] < 0

-19-
Next we have

:_uj uJ+1 dw (uj+i/2-w) f'(uj+I/2)

_u_.x_wu_,__s_ __._.f,.,c_d_ "U . "W

3

:-;j+1"U. dw _ j+I/2(s-w) f"(s)ds

3

: _ ___/2_J*1 u -u.ds f"(s) (minls - ui*1o_Is, - ujl)2 < O. 3

Suppose uj > uj+I, i.e. we are at a shock. Let _j+I/2 be chosen

so that

_j+l f'(w)(w dw 0

j - aj+i/2) :

ioeo

_j+1 wf'(w) dw

F(uj)

C3.4) gj+'n,,/2 : _ :

I\ f(uj) "uj+1 f'(W) --+ U.

3

-20-
where F(u) is the entropy flux corresponding to our entropy,

V(u) = I/u2. 2 Of course we need the denominator in (3.4) to be nonvanishing. In fact we restrict ourselves to the situation where

(3.5) uj+ I _ _j+I/2_ uj.

We note

_]j+1 J (w- uj+I/2)f'(w) dw

(3.6) Gj+I/2- uj+i/2: /_ f(uj)

uj+1 (w - uj+i/2) (f'(w) - f'(uj+i/2))dw

/__ f(uj)

Thus,we have

sup f"(u)

, I 2 uj+1_u!uj

(3.7) :_j+I/2- uj+i/2,< _ _ uj) :sj+i/2: ,

where sj+i/2 /__ f(uj)

We further restrict the sj so that

(3.8) 0 _< -Asj _<2 max(min((uj - gj+I/2),(_j_l/_2 uj)),0).

In view of (3.7), this restriction does not affect the second order

-21-
accuracy of the scheme, except at sonic points.

We now write

J

lU /'-ksj+1 :_ uj[h j+1 _, Uj + - f(Gj+i/2)]

uj

dw [f(w) - f(Gj+I/2)] +_uj+ I

=[III] + [IV].

Now, the assumptions (3.5), (3.8) imply

(3.1o) [izz] a< o.

An argumentanalagousto that of (3.3) sufficesto show that

[zv]<0.

Finally,suppose uj+I < uj and (3.5) is violated. The inequalities of Lemma (2.1),and (3.8) implythat sj and sj+I vanish. Then

inequality(1.21)is immediate,since h is an E flux.

We may summarizeall thisin

-22-
THEOREM 3.1. Th___esequence of approximate solutions satisfyin_

(2.6) converges a.e. to the unique solution of the scalar convex conservation law (1.1) provided that the initial data is in B.V and that for

eachj

_< I

° </__%uj -

o < _sj < I

-_uj -

and, in addition if uj > uj+1,

-/Ssj _ 2 max (min ((uj - Qj+I/2) , _j_i/2 - uj)), 0).

IV. Steady Discrete Shocks. We now check for the existence of

discrete, steady, shock solutions to the convergent MUSCL schemes constructed in the previous section, based on the Engquist-Osher flux, [4].

Since the scheme satisfies a maximum principle (Lemma (2.1) of [12]),

any profile must be monotone.

Let f"(u) > O, and for simplicity, we take f'(O) : 0 : f(O). In

this case, the monotone E-O scheme becomes

(4.1)(a) hEO (uj+1,uj) : f_(uj+I) + f+(uj)

where

-23-
(b) f (u)--f(u) if u < 0

f (u)-_ 0 if u > 0

(c) f+(u)-_f(u) if u > 0

f+(u)-_ 0 if u < 0

Let uL, uR be the left and right states for a physically correct,

steadyshock, i.e.

f(uL) : f(uR), uL >O>u R.

A steady discrete shock {uj} will satisfy

L

(4.2)(a) lim u_ : u

J

J_oo

R

(b) lim u. : u

j-_oo 3

(c) hEOIu _sj+1 uj AsJl . f(uL) lJ+1- 2 ' + 21-

for all J.

-24-
We shall seek steady, discrete, shocks of the same general form

obtained in [4] for the E-O scheme, e.f. also [12], section V.

These are

(4.3)(a) u. _-uL, j < -I j

uR

(b) uj _- , j _> 2

(c) u L _> u0 > 0 > u I > uR.

For the first order scheme, u0 can be viewed as a smooth function

of uI satisfying the above, and, in addition

f(uO) + f(uI) = f(uL)

For the present scheme, we shall get a different one parameter family of

intermediate states.

It follows for all j _ O, that (4.3) implies (4.2)(c).

For j : 0, we have the following equation:

One special solution is, again, u I : u R, u0 : O.

For (Uo,UI) close to (O,uR) and satisfying (4.3)(c), we have,

-25-
from the hypotheses of Theorem (3.1):

(uR _ uI)

(4.5) Sl = al/2 /k

So:O

Here 0 _ al/2 is assumed, in addition, to be a C2 function of

R _a_

UI near uI : u , with u_1: O.

A straightforward calculation gives us

(4.6)(a) _ (u 0) : f'(uR) I + < 0

_G

(b) _ (uR,o) : 0

(c) _)2G _ (uR,0) = f"_(UR) I + > 0

5u I

_2G

(d) _ (uR,o) = 0

(e) _2G = f".(0)> 0

 Uo2

A simple application of the implicit function theorem gives us:

-26-
THEOREM 4.1. Given the hypotheses of Theorem (3.1) and the mild

technical assumption concernin5 ai/2, there exists a family of sharp,

discrete, shock solutions to (2.6) of the type (4.3), with u I a smooth

function of Uo, and 0 < Uo, small enough. By symmetry, the same i__ss

true, with u 0 a smooth function of u I < O, and -u I small enough.

V. Com arison With Other Convergent High Resolution Schemes. In

[12] we constructed a class of convergent high resolution schemes, using

flux limiters and second order accurate upwinding. Our first step was

to use a first order E scheme flux, h(uj+1,uj) , to construct a second

order accurate TVD scheme. The resulting scheme is

_u.

(5.1)

:- H(uj.2,ujI,. uj,uj_1),

with

(5.2) H(uj+2, uj+ I, uj, uj_I) = h(uj+I, uj) - I/2_ (Rj+1)(h(uj+1,-uj) - f(uj))

+ I/2 _(R_)(f(uj+ I) - h(uj+ I, uj)).

Here

(5.3) R_ _ f(ul) - h(uj, uj_1)

J f(uj+I) - h(uj+I, uj)

R_ _ h(uj+1,u j) - f(uj)

h(uj, uj_1) - f(uj_1)

-27-
The slope limiter @ (R) satisfies _(I) : I and the following

inequalities

]

-4 - ) >o

(5.4)(a) I + I/2 R- (Rj+I - J

['_ (R+) 1

1/21"-2-+-1_(R o

'+ L5 >--

Various slope limiters have been developed. See Sweby [13], for a

numerical and theoretical analysis of their properties. Perhaps the

simplest example is the so called min-mod, [14]:

(5.5) _(R) : max (O,min(R,1)).

We have the following analogue of Lemma (2.2):

LEMMA 5.1. If h is an E flux, then the viscosity of (5.1) is a

decreasin_ function of 9.

Next we compare the schemes (2.6) and (5.1). If f'(w) g 0 and w

is smooth at xj, then it is easy to see that

_ uj _2--_] O( uj)2 uj)2) H(uj+2, uj+1, uj, uj_I) : h Uj+ I - 2 ' + + _ +

if

-28-
(5.6) A sj : @(R_)Z_" uj + 0((_L uj)2 + (_ uj)2), if h1(w,w) _ 0

and/or

sj=_(R])_ uj.0(_ uj)2. (_uj)2),if ho(w,w) _0

Dropping the quadratic terms, inequality (2.9) gives

_-_ R_) < I if hI

(5.7)(a) 0 _ R_ ' g( - (w,w) _ 0

(b) 0 <_(-_

-- . , g(R_) _ I if h0(w,w) _ 0

Rj

These inequalities are a bit more restrictive than those needed in

the monotone case. There we have

(5.8)(a) I _ I/2 I-_ + _(Rj+I) 1 if h1(w,w) _ 0

1 0 0

which are the same as inequalities (5.4)(a),(b).

For the entropy condition to be proven valid, the flux was modified. We replaced H by Hac (where a.c. stands for artificial

29
compression), as follows.

First we modify the flux differences in (5.2):

(5.9)(a) (h(uj+1, uj) - f(uj))M

: (h(uj+1, uj) -f(uj)) I + "h(uj+1,hICu uj)j,uj)) -f(uj) uj)1

(b) (f(uj+I) - h(uj+I, uj))M

a- I,(Ah^(u , u ))A u

_-,/:_"--+ u J J

= (f(uj.I) -h(uj+1, uj)) I - _uj+ 11 -h(uj+1, uj)3

.

Here a_+1_ are both positive, chosen first so that the quantities

in the brackets in (5.9) are between 0 and 2, (inequalities (5.13)

below).

We next let:

M (f(u) -h(u , u ))M

+ j j j-1

Rj :

(f(uj+I) - h(uj+ I, uj))M

M_ (h(uj+1, uj) - f(uj))M R. :

3 (h(uj, uj_I) - f(uj_1))M

-30-
Then we define our scheme:

(5.10) _t : - i/-_ Hac (uj+2, uj+I, uj, uj_I)

with

(5.11) HaC(u.+^, 3 _ uj+1, uj, uj_I)

M

: h(uj+1, uj) - I/2_(Rj+I) (h(uj+1, uj) - f(uj))M

M

+ I/2$(Rj+) (f(uj+I) -h(uj+ I, uj))M.

The resulting scheme is convergent if the following inequalities

are valid:

(5.12) _uj+1 f"(w)dw ((I/2_ uj))2 - (w - I/2(uj+I + uj))2) uj

- a;+l_ (_/_ uj)2 [h1(uj+I, uj+I) - h1(uj, uj)]

- a+j+_I (A-_uj)2 [ho(Uj+1,uj+1) - ho(Uj, uj)] < 0

(h(uj+I, uj) - f(u.))

(5.13)(a) la;+I/2(h1(uj+I, uj+ I) - h1(uj, uj))l _ - /_ uj J

+ f(uj+I) - h(uj+ I, uj)

(b) Iaj+I/2(ho(uj+I, uj+I) - ho(Uj, uj)): ! /_ uj '

-31-
and if _(R) is definedby (5.5).

Examplesare given in [12].

It is interestingthat this case allowscompressionto be added at

+

shocks, i.e. the restrictionson the a3+i/2allow negativeviscosityto

be added to (5.2) if uj > uj+1, while, again for the proof of convergence,positiveviscositymust be added for rarefactions,uj < uj+I. In

the presentMUSCL ease the situationis (unnaturally)reversed. A comparison of the convergentMUSCL schemewith this one yields two main

points:

(I) Away from sonic points,the fluxesdiffer only in 0((__ uj)3)

terms

(2) At sonic points,the fluxesdiffer in 0(__ uj)2) terms.

BIBLIOGRAPHY

[I] S. R. CHAKRAVARTHYand S. OSHER, High resolutionapplicationsof

the Osher ugwind scheme for the Euler equations,Prec. AIAA

Comp. Fluid DynamicsConf., Danvers,Mass. (1983),pp. 363-372.

[2] P. COLELLA,A directEulerianMUSCL schemefor gas d_namics,

LawrenceBerkeleyLab. Report# LBL-14104(1982).

[3] R. J. DIPERNA,Convergenceof approximatesolutionsto conservation laws, Arch. Rat. Mech. and Analysis,82 (1983),pp. 27-70.

[4] B. ENGQUISTand S. OSHER, Stableand entropyconditionsatisf[in5

approximationsfor transonicflow calculations,Math. temp.,34

(1980),pp. 45-75.

-32-
[5] S. K. GODUNOV, _ finite difference method for the numerical computation of discontinuous solutions of the equations of fluid

dynamics, Mat. Sb., 47 (1959). pp. 271-290.

[6] B. GUSTAFSSON, The convergence rate for difference approximations

to mixed initial boundary value problems, Math. Comp., 29

(1975), pp. 396-406.

[7] A. HARTEN, High resolutlon schemes for hyperbolic conservation

laws, J. Comp. Phys., 49 (1983), pp. 357-393.

[8] A. HARTEN, On second order accurate Godunov-ty__peschemes, NASA

_{ES Report # NCA2-ORS25-201.

[9] S. N. KRUZKOV, First order quasi-linear equations in several

independent variables, Math. USSR Sb., 10 (1970), pp. 217-2q3.

[I0] W. A. MULDER and B. VAN LEER, Implicit u_wind computations for the

Euler equations, AIAA Comp. Fluid Dynamics Conf., Danvers,

Mass., (1983), pp. 303-310.

[11] S. OSHER, Rlemann solvers, the entropy condition, and difference

approximations, SINUM, (to appear).

[12] S. OSHER and S. R. CHAKRAVARTHY, High resolution schemes and the

entropy condition, SINUM, (to appear).

[13] P. K. SWEBY, High resolution schemes using flux limiters for

hyperbolic conservation law____s, SINUM (to appear).

[14] P. K. SWEBY and M. J. BAINES, Convergence of Roe's scheme for the

general non-linear scalar wave equation, Reading University

Numerical Analysis Report 8/81.

-33-
[15] E. TADMOR, Numerical viscosity and the entropy condition for conservative difference schemes, NASA Contractor Report 172141,

(1983), NASA Langley.

[16] B. VAN LEER, Towards the ultimate conservative difference scheme,

II. Monotonieity and conservation combined in a second order

scheme, J. Comp. Phys., 14 (1974), PP. 361-376.

[17] B. VAN LEER, Towards the ultimate conservative difference scheme.

V. A. second-order sequel t_ooGodunov's method, J. Comp. Phys.,

32 (1979) PP. 101-136.

[18] H. C. YEE, R. F. WARMING, and A. HARTEN, Implicit total variation

diminishing (TVD) schemes for steady state calculations, Proc.

AIAA Comp. Fluid Dynamics Conf., Danvers, Mass., (1983), PP.

110-127.

[19] R. SANDERS, On convergence of monotone finite difference schemes

with variable spatial differencing, Math. Comp., 40 (1983), PP.

91-106.

-34-



1. ReportNo. 2. GovernmentAccessionNo. 3. Recipient'sCatalogNo.

NASA CR- 172306

4. TitleandSubtitle 5. ReportDate

February 1984

Convergence of Generalized MUSCL Schemes 6. PerformingOrganizationCode

7. Author(s) 8. Performing Organization Report No.

Stanley Osher 84-10

10. WorkUnit No.

9. Performing Organization NameandAddress

Institute for Computer Applications in Science

and Engineering '11. Contract or Grant No.

Mail Stop 132C, NASA Langely Research Center NASI-17070

Hampton, VA 23665 13. Type of Report and Period Covered

12. Sponsoring Agency Name and Address

contractor report

National Aeronautics and Space Administration 14. Sponsoring Agency Code

Washington, D.C. 20546

5.SupplementaNryotes Additional Support: NSF Grant No. MCS 82-00788, ARO Grant No. DAAG

29-82-k-0090, NASA Grant No. NAG-I-270.

Langley Technical Monitor: Robert H. Tolson

Final Report

16. Abstract

Semi-discrete generalizations of the second order extension of Godunov's scheme,

known as the MUSCL scheme, are constructed, starting with any three point "E" scheme.

They are used to approximate scalar conservation laws in one space dimension. For

convex conservation laws, each member of a wide class is proven to be a convergent

approximation to the correct physical solution. Comparison with another class of

high resolution convergent schemes is made.

17. Key Words (Sugg_ted by Author(s)) 18. Distribution Statement

high resolution 02 Aerodynamics

MUSCL approximation 64 Numerical Analysis

conservation laws

Unclassified-Unlimited

19. Security_a_if.(ofthisreport] 20. S_urit¥Classif.(ofthis_ge) 21. No. ofpages 22. Dice

Unclassified Unclassified 36 A03

Forsale by the NationalTechnical InformationService,Springfield,Virginia 22161

NASA-Langley, 1984

NASA Technical Library

3 1176 01416 9362