The discontinuous finite element method with the Taylor-Galerkin approach for nonlinear hyperbolic conservation laws

Computer Methods in Applied Mechanics and Engineering 95 (1992) 141-167 North-Holland

The discontinuous finite element method with the Taylor-Galerkin approach for nonlinear hyperbolic conservation laws Kyu Y. Choe and Keith A. Hoisapple Department of Aeronautics and Astronautics, University of Washington, Seattle, WA 98195, USA Received 2 April 19911

A one step, explicit finite element scheme which is second order accurate in both time and space is developed for the computation of weak solutions of nonlinear hyperbolic conservation laws in one dimension. The scheme is an improved version of the discontinuous finite element method using the Taylor-Galerkin procedure. It is linearly stable under a fixed CFL number up to nearly 0.4 (or 1/2 for the alternative two step explicit scheme also developed here) and TVDM which guarantees convergence to a weak solution in nonlinear problems when the flux limiter is applied with the CFL number up to 0.366 (or 1/2 for the two step scheme). Numerical experimentation demonstrates the convergence of the solution to the entropy one even for nonconvex flux cases. The scheme captures stationary discontinuities perfectly.

1. Introduction

The finite difference method is the most widely used technique for the computation of hyperbolic equations involving conservation laws in fluid dynamics. Schemes in this method can be grouped largely into upwind and centered difference types. Among the centered difference schemes, the Lax-Wendroff [1] scheme is a typical one. It utilizes the Taylor series expansion in temporal derivative terms for second order accurate time integration and then replaces the appropriate temporal derivative terms with the spatial ones using the governing equation. The viscous-like terms resulting from this approach not only give the scheme the higher order accuracy but also make the scheme stable. However, the Taylor series expansion assumes the solution field to be smooth, and, due to the small amount of numerical viscosity, the scheme shows spurious oscillations near discontinuities of the solution. In nonlinear cases, it may even fail to converge to a weak solution. Explicit artificial viscosities [2-4] can be used to remedy this problem. The upwind scheme was pioneered by Courant et al. [51 who used one sided differencing in space in the characteristic direction to account for the direction in which the signal propagates. This idea has recently been extended to popular schemes for systems of nonlinear equations such as flux limiter schemes [6]. The Godunov scheme [7], which Harten et al. [8] put into the category of upstream methods in their review of finite difference methods, is unique in its approach. It is neither a finite difference scheme nor a finite element scheme. The basic building block of the method is the exact solution of the local Riemann problem resulting from the projection of the initial condition onto the space of piecewise constant functions. In the special case of a linear flux function, the resulting scheme is identical to the C-I-R scheme. This scheme was extended to second and third order by Van Leer [9] who used piecewise continuous linear and quadratic functions for the approximate solution. However, these 0045-7825/92/$05.00 ~) 1992 Elsevier Science Publishers B.V. All rights reserved

142

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

schemes are not monotone and, in nonlinear cases, one must solve more than a Riemann problem at each discontinuity due to the nonconstant base functions for the approximate solution. The first shortcoming could be remedied by a flux limiter devised to enforce the monotonicity of the scheme, but the second problem makes the scheme impractical because, in general, there does not exist a simple solution for a nonlinear problem with a discontinuous initial condition. The finite element method has demonstrated its power in a variety of solid and structural mechanics problems [10]. The ability to treat boundary conditions naturally, and the ease of handling complex geometries has made the method attractive. The natural treatment of the boundary terms from the weak formulation is especially suited to elliptic problems and the freedom of arbitrary element shapes gives the method the ability to handle any kind of complex geometry. Furthermore, the centered difference form natural to the Galerkin formulation results in higher order accuracy than it does in corresponding finite difference schemes [11, 12]. Unfortunately, the finite element method loses the first of these merits in hyperbolic equations, where both the boundary and the initial conditions are required, and no conditions are specified at the final stage of the solution. However, the second merit remains true for all kind of problems. That could be the most attractive feature of the method in application to computational fluid dynamics problems, due to the difficulty of mesh generation in complex geometries and the resulting effects o,! the solutions. As is well known in finite difference methods, the centered difference form in space combined with forward Euler time steps results in an unconditionally unstable algorithm. This is again due to the role of characteristics as preferred directions of propagation. To fix this problem, various approaches have been presented. Among them, and corresponding to the upwind scheme in the finite difference method, are those with 'biased' weighting functions [13-15]. Morton and Parrott [16] also used modified weighting functions devised to satisfy the 'unit CFL property' for undistorted wave propagation when CFL = 1 to produce an upwind effect. Hughes and Brooks [17] introduced upwinding only in the streamline direction to devise the streamline upwind Petrov-Galerkin method. That method has been analyzed by Johnson et al. [18-20]. These Petrov-Galerkin methods choose their weighting functions from a function space different from that of the base functions to produce an upwind effect. The latter two schemes are known to be particularly effective, but they invariably involve a free parameter to determine the optimal form of the weighting function. On the other hand the Taylor-Galerkin method introduced by Donea [21] preserves the Gal~rkin approach to keep higher order accuracy both in space and time, but also adds stability to the scheme by following Lax-Wendroff's idea. Therefore, the Taylor-Galerkin method corresponds to the centered difference scheme in the finite difference method. In this approach the time discretization is performed prior to the spatial discretization by a Taylor series expansion to achieve higher order accuracy, and then the governing equation is used to replace the augmented temporal derivative terms with the spatial ones which gives damping in the resulting scheme. This scheme does not involve a free parameter and results in a natural form with the optimal value of the free parameters appearing in the Petrov-Galerkin methods. It also shows a particularly high phase accuracy but, due to the inherent small amount of numerical viscosity and the underlying assumption for the Taylor series expansion that the solution field be smooth, the scheme results in a wild oscillation near discontinuities of the solution and, in multidimensional nonlinear problems, may even fail to converge to a weak solution. This kind of problem could be improved by the Special treatment of the, discontinuities like the pseudo-Galerkin method [22] or, as usual, by adding explicit artificial viscosity to the scheme [23]. All the finite element methods described above use continuous functions for their bases and

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

143

weighting functions and have the disadvantage that, even with the forward Euler time discretization, the final discrete equation is implicit due to the coupling of the nodal variables in the consistent mass matrix. A lesser known finite element scheme, called the Discontinuous Finite Element Method (DFEM), uses only piecewise continuous functions for base and weighting functions and hence has a potential ability to isolate the inter-coupling of the nodal variables within one element. The first appearance of this scheme is in Reed and Hill's report [24, 25] on ti~e linear neutron transport problem. Rigorous mathematical analyses were given by Lesaint and Raviart [26], and Johnson and Pitkiiranta [27, 28]. For these linear problems, the characteristic direction is fixed (constant) in space, and the inflow and outflow boundaries could be easily determined prior to the solution process. Therefore the problem could be successively solved element by element. This idea was further extended to nonlinear conservation law problems in petroleum engineering by Chavent et al. [29, 30] with the introduction of the Godunov value [31] for the predetermination of the inflow element boundary condition. The approach is similar to the hybrid method [I0] with additional interface variables, but with the difference that the interface variable is predetermined by the exact solution (Godunov value) of the local Riemann problem at each element interface where the solution is discontinuous. In this sense, the Discontinuous Finite Element Method is the finite element version of Van Leer's higher order extension of Godunov's scheme. But, unlike the Van Leer scheme, which requires the exact solution of the general discontinuous initial-boundary value problem, the DFEM needs to solve only the local Riemann problem to determine the inflow interface (Godunov) values. This makes the scheme practical even for the higher order basis functions. In this way the DFEM produces an explicit scheme which is first order accurate in time (with the forward Euler time discretization) and second order accurate in space. Unfortunately, the scheme is unconditionally unstable under a fixed CFL number and is stable only under the restrictive condition C F L - O ( ~ ) [32]. To enforce stability in nonlinear problems a flux limiter is combined with the scheme to give the TVDM (Total Variation Diminishing in the Mean) property under a maximum CFL number near 1/2 which guarantees the convergence to a weak solution [33]. But, even with the use of the flux limiter, extensive numerical experiments show the convergence of the scheme to a nonentropy (but weak) solution [33] and also, in nonconvex flux problems, the scheme is still anti-diffusive [34]. Recently a Runge-Kutta iteration based scheme was introduced for higher order time accuracy and improved linear stability giving general practical applicability [35]. This paper describes an alternative method to improve both' the time accuracy and stability based on the following interpretation on the instability of the DFEM: (1) Without the use of the Godunov value, the scheme takes the centered difference form in space: a;'+' = u,_. + A t ! f l U "' + ' )

-/(u;')}.

U[n

-W

111

lit

Ill

up

7,".

u.~l

I

I

I|1 Ui+l

i

I

I!l!

I

~'

~,._ . . . . . . . . . . . . . . node i

uZ'.

£1emen[ i

~1111 - ~ --- X node i+ I

Fig. 1. Stencils used in D F E M without the O o d u n o v values defined.

144

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method EP-t = u;"= ur

/

, #'

D

[l~ ; illU~

[ i i I J

II[I II

&

/~"

node i-I Element

l-I

41

__

U

Ui+! - Ui+--~l

:

i i I I

I I,. +. /lui+l I11 I I~x

"J

n o d e i+ l

node J

Element

i

Fig. 2. Stencils used in DFEM with the Godunov values for f'(u) > 0 when the piecewise constant base functions are used.

(2)

I n the case

f'(u) >

0 with piecewise constant base functions the use of the Godunov value

at each interface gives the scheme an upwind effect because u 7 = 6i_ t -n+

1

At

= a7 + At

= t/7 + ~ (3)

-n

{i(u,+,)-

f(u?")}

6" {f(t/~')-f( i-,)}.

But the use of higher order base functions puts the scheme back to the centered difference form because u~- ~ 6i_ i.

a;'+'= a;'+

At

-n

{/(u,+,)-i(u?")}.

Therefore the resulting scheme combined with the forward Euler time discretization becomes unconditionally unstable. This observation naturally leads to the consideration of the application of the Taylor-Galerkin approach to the scheme for both time accuracy and stability. The resulting scheme, which will be called Taylor-Galerkin Discontinuou~ Finite Element Method ( T G - D F E M ) , is a purely explicit, one step method (in contrast to the R u n g e - K u t t a D F E M which is iterative) with an acceptable maximum CFL number near 0.4 for linear stability and preserves the ability to perfectly capture stationary discontinuities of the solution. Also like the D F E M or RKD F E M , flux limiters can be easily applied to give a T V D M property to the scheme. Numerical experiments show the improved overall performance of the suggested scheme.

u;"

= ur ,.-,f-"

u[~l

~r

I

II

II

I

l

I

\tC ,I

node i

:

Ui+ !

"+~ ui+1

=-x

Element i

node i + !

Fig. 3. Stencils used in DFEM with the Godunov values for f'(u) > 0 when the discontinuous piecewise linear base functions are used.

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

145

2. Mathematical formulation of the discontinuous finite element method

Consider a nonlinear scalar conservation law

ou Ot

of(u)

+

=o,

Ox

(1)

with the initial condition u(x, O) = Uo(X)

Vx ~ O ,

(2)

and the boundary condition

u(a, t) = ua(t )

and/or

u(b, t)

-

-

Ub(t ) or none

(3)

depending on the sign of f'(u) at the boundaries of ~ , x = a and x = b. For an approximate solution to this problem, discretize the space domain ~2 = (a, b) into I intervals (x~,x~+l), i = 1, 2 , . . . , I. Define a function space, V h = { v h E?. L2(/~

)" Vh[t~,.~,., ) E ~k, i = 1, 2 , . . . ,

(4)

I~,

and the approximation of the initial condition u 0, Uoh ~ Vh ,

(5)

where ~ k is the space of polynomials of degree k. For simplicity and practicality, choose k = 1 in (4), i.e. choose discontinuous piecewise linear polynomials for base functions as follows: ÷ Vi=

f(x,+, - x)/h,, [ O,

_ f (x - xi)/h,, Vi+! : ][ 0 ,

if x E [x i, xi+ 1], if x ~'[x,, x,. I] ;

(6)

if x • [ x i , x i . l ] , if x ~'[x,, x i* I],

where h~ = x~. ! - x, is the element size. With these base functions, the approximate solution can be written as I Uh" = ~

•tu,+n

(7)

v,+ + u/+~v,+l],

i=l

Uh

I I I.

.I

I I

B

I I

I I

0

~X Xi-2

Xi-I

xi

hi

Xi+l

xi+2

Fig. 4. Shape functions for the Discontinuous Finite Element Method.

146

K . Y . Chae, K . A . Holsapple, Discontinuous Taylor-Galerkin finite element method Uh

l

U.+

I

I

I

I /"

I

Uhll

ui IS.~

I

l u¢+i

-

I I

ui'-! I

~..+ '+ZLU'+2 I~

I

I I I I

I I I

Xi-2

Xi-1

I I I I X

I I I

Xi+l

~x

Xi+2

Fig. 5. Discontinuity of the approximate solution u h. II

-

11

II

II

where u [ " = Uh(X ~ + 0 ) and u~ =Uh(X~--0 ). Note that the approximate solution u h is discontinuous at each element interface x~ where it has two limiting values u~" and u~-" as shown in Fig. 5. In addition to the use of discontinuous base functions, a distinguishing feature of the discontinuous finite element method is the introduction of the so-called Godunov value at each element interface where the approximation function u h is discontinuous. The Godunov value ~'~' is defined by •

/

slgnkui

+ n

-

- ui

n

n

•

/

)" f ( ~ i ) <-slgn(ui

VkE(u~-",u~")

+ n

il

- u:, ). f ( k ) ,

(8)

Vi=l,2,...,l+l,

and can be shown to be the exact solution of the Riemann problem at x = x i [29]. The use of the Godunov value at each element interface yields an upwind effect [7, 36] to the numerical scheme. In problems where the exact solution cannot be easily determined, simple approximate Riemann solvers such as those of Roe [37] or of Engquist and Osher [38] can be used to determine an alternative to the Godunov value. With the definition of the Godunov value ~:i, the weighted residual statement for the semi-discretized problem determines the function u h E Vh by u h ( O ) = Uoh , "'+' OUh

J~,

Ot v d x + - - ?!

fxri + I

n

Of(Uh)

~

,..

+.

v d x + { .ttui

Ox

t!

-{f(u,+,)-f(~,+,)}v(x,+,)=O

Vv=•*

,,

) - f ( ¢i )}v(xi) Vi=l,2,...,l,

(9)

which is in turn fully discretized by the forward Euler time discretization approximating the temporal derivative term in (9) to yield [29]

( xi +I Uh"+' J~~

-- ii h"

At - - I!

-- { f ( U i + l ) - -

["+' o d x + J.'i

0f(u;:) v d x Ox

+ { j ~" t u , +" ) -

n

f(~i+i)}V(Xi+l)--O,

(10)

where u'/, = Uh(t,) and t 0 = 0, t 1, t 2 , . . . , t N = T is the partition of the finite time interval

[0, T].

In the approximation procedure, the quadrature rules used for the integration of the first term in (10) affects the stability of the resulting numerical scheme considerably, as will be

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

147

shown later. With the definitions in (6) and (7) and with Simpson's quadrature rule for the integration of the second (flux) term of the equation, the weighted residual formulation (10) leads to the following fully discretized difference equations under uniform space discretization: u+n+ l

= u~ n - Y[+/z{f(u~ n) + 4f(uT+,/2) + f(u,-+~,)} - ( 3 / z - 1)f( ~:7+i) - (3/~ + 1)f(¢7)1,

-n+!

Ui+i

(11) --?!

--n

n

fl

= ui+ , - T[-/~ {f(u +') + 4f(ui+ 112) "[-f(ui +, )} + (3p + 1 )f( sei+, ) +(3~ -

1)f(~7)l,

where /~ =

1 1/3,

if Simpson's quadrature rule is used, if the trapezoidal quadrature rule is used,

for the integration of the first (temporal) term, y = At/h where h is the uniform element size = + -.+, and Ui+ll 2 u i + | )/2. The major advantage of the above scheme over the conventional (continuous) finite element method is seen from (11): the scheme is purely explicit and simple to code. But unfortunately the CFL number for linear stability varies as ~ , which is too restrictive for practical applications [32]. To avoid this difficulty, Chavont and Cockburn [34] employed Van Leer's [9] flux limiter algorithm as follows: let

~i=(uT,+,-u:)/2

and

ti:=(u~-+,+u

+

i)/2.

Then the flux limited solution of the discontinuous finite element scheme is defined as (see Fig. 6)

ui-*=ui,

ui~*=Ominm°d(A-ui, u~, A+ui) , (12)

Ui

_, _ a, ~ U i

a,+,

,

=

-, + U i

,

whore A+ and A_ are the usual forward and backward difference operators, respectively, 0 is a positive real number (usually 0 ~< 0 ~< 1) and min mod is the usual rain mod function defined as man mod(a, b, c)

fs min(lal, Ibl, Icl), [ 0,

if sign(a)= sign(b)= sign(c)= s , otherwise.

Uh

X xi-2

xi-I

xi

Fig. 6. Flux limiter.

Xi+l

Xi+2

148

K.Y. Choe, K.A. Holsapple, Discontinuous Tayior-Galerkin finite element method

It is proved in [33] that the resulting scheme is TVDM (Total Variation Diminishing in the Mean) under the condition 0 ~< 1 and -/~< 1 / 2 which guarantees the convergence of the scheme to a weak solution even for nonlinear problems. In the limit 0 = 0, the scheme is identical to the Godunov [7] scheme. Despite the introduction of the flux limiter, the solution of a nonlinear problem with a nonconvex flux function f(u) is still anti-diffusive [34] so that the numerical shock front lags behind that of the exact solution, and even the smooth part is affected giving inaccurate results. There have been attempts [34] to supply enough damping to the numerical scheme by adjusting the flux limiter parameter 0 to some optimal value between 0 and 1, but the 'optimal value' of 0 is highly problem dependent and there is no general rule to predict it except by a numerical simulation for each problem.

3. Taylor-Galerkin discontinuous finite element method

The discontinuous finite element technique is unconditionally unstable. Even though it is convergent to a weak solution with the flux limiter (with 0 = 1), the scheme does not have enough numerical dissipation to produce an entropy solution [33]. Whereas the definition of the Godunov value at each element interface provides an upwind effect to the scheme, the use of higher order (k I> 1) basis functions brings the scheme back to the form of centered differencing in space. As is well known, the centered spatial differencing with forward Euler time stepping results in unconditionally unstable algorithms in the finite difference context. A Von Neumann stability analysis shows [32] that the eigenvalues of the magnification matrix of the discontinuous finite element scheme are larger than 1 in the vicinity of/3 = 0 irrespective of the CFL number (/3 is the wave number multiplied by Ax). To improve this stability problem, Lax-Wendroff's idea is introduced into the discontinuous finite element scheme for additional dissipation as well as higher order time integration resulting from a Taylor series expansion of the temporal derivative term. Following the approach of Donea [21], take the forward Taylor expansion of the temporal derivative term to obtain ,

U" + I - u "

u, =

At

At n At2 n 2 u , , - ~ u,,, + O(At3).

(13)

For simplicity, terms through second order are included, with discontinuous piecewise linear base functions. With the aid of the governing equation, (13) can be rewritten as un÷I

u',' =

-- u n

At

At 0 ( Ou)l" 2 Ox A' -~x + O ( A t 2 ) '

(14)

Therefore (1) is written in the semi-discretized form

u"+~-u"At

-

-

At ( 2-Ox -O A2 Ou)["

-

Of

(15)

and the corresponding weighted residual statement for the discontinuous finite element method is to find a function u h E V~ which satisfies

f x,+, u,, +'At-Uh v dx - -~At rJ~,x,+, OxO ( A2 -~xOU)1"v d x + jx,f ~'+' Of(US',) Ox v dx jx~

j~

{

~,,I

jtu,

"4-t!

--

H

II

.--

o

°

(16)

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkinfinite element method

149

The difficulty in handling the second derivative terms occurring in the additional diffusion term can be easily overcome by the use of the mixed finite element method. Introduce an additional variable

At A 2 0 u r

~

~

2

(17)

Ox

or equivalently f f rw d x = - - -At 2 ;

A2 ~X OU w d x

(18)

Vw~_ ~ l o,

where a and b are the boundaries of the domain. With the choice of continuous piecewise linear functions for w and the trapezoidal quadrature rule for the integration, the Galerkin approximation of (18) gives the difference formulas for r: r~ = - ~ T A 2 ( u ~ ) ( u ~ - u ~ ) T [A2(u;)(u~. _ u i + ) + A2(u~.)(u~.+, - u j )+] , rj = -- -~

j=2,3,...,I,

(19)

+

rt+ i = ---~T A2(U-~+I)(U-~+I_ U! ) Note that the right-hand side of the second equation represents the arithmetic average of the left and right limiting values of - ½ A t A 2 ( u ) u x at the interface. With this introduction of the new variable r in (17) the weighted residual form of the Taylor-Galerkin DFEM is written as

f~'+ ' unn+' J~, At

Uh

f~'+' Or n

v dx + j~,

.... Ox v dx -

+ f ( ~;'+, ) o ( x , +, ) - f ( ¢ 7 ) v ( x , )

ff+' oo , f(Uh) ~x dX (20)

= o.

Proceeding as in Section 2, the following fully discretized difference scheme is obtained: n

n

--n

u +n+l = u ; n - Tl{r,+i - r i ) + Id,{f(u; n) "t" 4f(uT+t,2) + f(u,+l))

-(3p,-

1)f(s~7+~) - (3p, + 1)f(~7) ] , (21)

-n+l

Ui+I

-n

n

n

-n

= ui+ l - T[{ri+l - rT} - / z { f ( u +n) + 4f(ui+l,2) + f ( u i + l ) }

+(3tL + I)f(~7+,) + (3tL - 1)f(67)1, with n T [A2(U~n)(u~n__ U,_l) +n -n + Ae(uTn)(ui+l -- U,+ n )] ri = ----4

where /z = 1 for Simpson's rule and # = 1 [3 for the trapezoidal rule in the numerical integration of the temporal term. In the linear case where f = u, a Von Neumann analysis shows that the eigenvalues of the amplification matrix of the scheme are •

;

A__(%/3) = 1 - ~T [(1 + 3p,) + (3/z - 1) e -i~] + ~,~/A(~,,/3), where

1 A(~,,/3)= ~ [(1 + 3/z) + (3# - 1) e-~]2 - 3 / z ( 2 - i~, s i n / 3 ) ( 1 - e - * )

150

K . Y . Choe, K . A . Holsapple, Discontinuous Taylor-Galerkin finite element method

and only A+ is relevant to the stability and accuracy. (A_ represents the 'Stegosaur bias' error mode, see [9]). As shown in Fig. 7, the only possible source of instability lies near the point/3 ~ 0 where the asymptotic expansion of A+ can be shown to be

A+m[l

(')/3)21-3 2 + 2-4

+ lOy-- 3y2 y3

(,y/3)4

] +

[ 1 / 6 - y/4 ] i (-y/3) + Y2 (yfl)3

(22)

as/3--->0. Equation (22) shows that the CFL condition for linear stability can be estimated from the relation

1

+

- 3 + lOy 3y 2 _ 1/3 -- y/2) 12y3 y2 (y/3)4 ~< 1 -

or

,)t3 + y2 + 2 y -

1 ~<0,

(23)

which gives y <~0.392647. Also an extremum analysis shows that [A+],= o = 1 , = ( I x ÷ l = V ,~=0=(]A+ ' I=V/~=o " =o

(24)

,

(IA+ 2V'"J~=o= - 6 y + 12y 2 + 6y 3 + 6y 4

and the condition for maximum (IA+[2) ''' ~<0 gives the same criterion as (23). It=I0 OriginalDI=EM L Locus

.....'....:.. i ,,,%

°'"' .... ]a=l/3, Original DFEM '~.' ' ' 'Locus

''i'''''''''"'''

i

..~o, p- 1/ 3, raylor-GalerklnDFEM ~.Locus

."

~'~

'k\

//

/

! : :

t t

:

-I

I

-15

:

-05

,

I

,

05

\ \

q Unii~C i r c l e

\ , : :

•

%',°

° .°'" • ,°°

..,° • °,, ,,° . . . . . . . . . ,° °°°,°

Fig. 7. Locus of the eigenvalues A~ (y = 0.39) on the complex plane.

151

K . Y . C h o e , K . A . Holsapple, Discontinuous T a y l o r - G a l e r k i n finite e l e m e n t m e t h o d

Figure 8 compares the modulus and the relative phase error of the amplification factor 4+ with those of other popular finite difference schemes. The modulus of the present scheme is closest to the unit circle which means the scheme is the least diffusive one. The relative phase error lies between those of Van Leer's second order scheme [9] and Lax-Wendroff scheme. The lagging phase error (inside the unit circle in contrast to the leading phase error of Van Leer's scheme which lies slightly outside the unit circle) shows much improved phase accuracy of the scheme compared to the Lax-Wendroff method. R E M A R K . Alternatively, a two step method described below can be used to increase the CFL limit of the scheme to 1/2.

For hyperbolic equations the characteristic speed is finite and therefore it is natural to use at, explicit scheme, whereas the characteristic speed of parabolic equations is infinite and it is natural to use an implicit scheme. However, for efficiency, the scheme is split into two explicit steps as follows: Step 1. Let ~,=(uT,+ ! - u ; ) / 2 --, u i -

--n

n

ylf(~:,+,) -f(~')}

u i ~' n

"¢ ~

and ~ , = ( u L ,

Then

+ u;)/2.

(25a)

. l|

u i = u, - / z y [ -

n

I|

(25b)

{ f(ui-2,) + 4f(u,+ ,,.) + f(u~+")} + 3 I f ( 6 i + , ) + f ( s¢` )}].

which is the first order DFEM. Step 2. Apply the diffusion term by ~n+ ! i --U* ~ -

"'+'

,/(r

?+

-ri), (26)

a*

Wi

=

,

where r ,*

-

=

~

[A"(u~*)(uT,

* -

u , _+* ,)

+

A2( u+*

) ( u , +-,*

-

u, +

* )l

=---Y [A2(u~-*)lT*, + A 2(u,+ *)a*] 2

.~-|.

UnitCircle TG-DFEM 1,,~ /'~"

..'b(/ /

=1

van Leer

/

x

~

/ ~,,,Lax-Wendroff

,.¢~"

,.'¢ To_o," T

\

¢'

///-,./

i I~

,,,,..... "" " " " ~

iS

1.

/

/

I t ~ ""

"~-.~..

_ ~ ~ . . ,

/

,'°.'] I

,'..,.o.','l

.

,,%,"~

'\

. .

~Lax-Wendroff

'~.

_

\~,nLoer

0.2

-0.5

,; XUnitCircle 0.5

(a)

I

,I -I

/

, -0.5

,

02 l "" - /

' 0.5

i

l

Co)

Fig. 8. Polar plot,,, of (a) the modulus [A+I and (b) the relative phase speed Arg(h+)/(-/t At) at y = 0.39.

152

K . Y . Choe, K . A . HoisapFle, Discontinuous Taylor-Galerkin finite element method

The amplification factor for one combined step is A±(/3) = 1 + ~T [(1 - e-'t3)(-1 + i3/z3' 2 s i n / 3 ) - 3/~(1 + e-iO)] +__' y ~ / 3 ) , where

1

A ( / 3 ) = -~ {(1 - e-~a)(-1 + i3/.~3'2 sin/3) - 3p,(1 + e-i°)} 2

- 30,(1 - e - i a ) ( 2 - i3" sin/3), which becomes A+ "" 1 + - i T / 3

. , /22 / 3 +2 i3'a/33 ( 121 q4 16.,/2 _.,/)

+ .,/,4/3.( 1

1 + 5 8.,/3 12.,/2

4

3.,/)

+ o ( ~ ~) as /3--*0. Therefore the scheme is second order accurate. An extremum analysis and an asymptotic expansion about fl ~ 0 shows that [A+ [ is maximum at/3 = 0 irrespective of the value of y but [A+ [ near/3 = ~r restricts the CFL number to y < 0.5. For nonlinear problems, the CFL condition for linear stability is only a necessary condition. To guarantee the convergence of a nonlinear problem, a flux limiter introduced in the previous section is again a convenient tool to ensure the T V D M property. To obtain the CFL-like condition for TVDM, rewrite (25a) in the form -n+l

ui

-n

-n

+n

n

= u, - 3"A_h(ui+t, u i + l ) - TA_ri+, ,

(27)

where A _ h ( tl i+nl ,

and

+n

n

(28)

n

u,÷,)=f(¢~+,)-f(¢,)

1 -","u,+, ÷. _ z . ÷. _. A_ r,,n+i = - g ~ ,t*~", ,% q. ,+ J~ -" - u, ) + , ~ t u , + , ) ( u , + , -- { ~

~ i

J[ iA21u'-n~'[u-n

+n

2

+n

2 Y['A2+'A+ui - "42A , a7]

+n

- u,

)}1

'

_ fl

= --

-n

- u,_,) + a (u,)(u,+,

+.

Ui+ 1)}

(29)

with r

,~/2 = ~ [ A e ( U ; n ) • =

[ a2(u7

-n

Ui

+n -- Ui-1

-,, a_u~

+

n ) /'~iLl ~ +1-1A2,' q,u i+n ) ~ /'~7

A_ui

-n

--2-. +n A tu i ) ]

A_uiJ

+n 1

Ui+ I -- Ui

a_a"

.

(30)

Following [39] or [40], rewrite (27) in the form -n+!

Ui

2

-" = u, + v[c,+,a+a"-

-,,

= ui +3"

"/

D,a_aT] + T

)_n(

3' .,~2+ a+u,

+l + ~

l

-

[,~L,a+a"-,~

D,

+~

-"

, a_u,]

,a_a7

--n

-- u, + ~ [ c , ÷ , a ÷ a 7 - fi, a _ a T l ,

(31)

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

where q + , = -h~,(1

a.aT

and

A"~+ti~']

D i = h,(1 + A_a~./

153

(32)

with the Lipschitz coefficients of h, h

--n

h I ---

-n

--n

, ui+n)

~n

//i+I -- Ui h

h2 --

+fi

(Ui+ I' Ui ) _ h(ui

-n

+n

-n

(33)

+n

(ui+,, ui+ t) - h(u i+,, ui ) +n

+n

Ui+I--U i

and

~ C+! ~-C+I + 23'/~2+, ,

/), = D, + ~3,/~2, •

(34)

For the Lipschitz continuous monotone flux h(., .) such as in Godunov's, Roe's or EngquistOsher's scheme, h~ ~>0 and h 2 ~ 0 [39, 41]. Now the TVDM property requires (i) C i, D i >~0, (ii) 3,(C~ ++ DD,)~ , ) ~ 1.1. (ii) 3,(C~ ~'n --n Under the condition [A +u~ [A+u~[<~ 1 which can be ensured by the flux limiter (with O = 1), it is easy to show that C~, D~ >~0 and, since ,A2 I> 0 by (30) with the flux limiter, it follows that C~, Di~>0. To satisfy requirement (ii), consider (30). Since 0 ~< ~i/A±a i <~1 under the flux limiter,

[

a,"_, + A2(u7 ") a7-n ] --~-[A2(u. ") + A2(u~'")] ~ 2A~,,,,

t~ i2 = A 2 ( u 7 n) ~

a_u~

where

Amax = umax[A(u)[,

a_u~J

,,~ = (m!n(uTn), max(uy')). I

(35) (36)

!

Therefore Ci+! + fii+l

= Ci+ t + Di+ ! + "~ (24/2+1+ 2~/2)~ Ci+ 1 + Oi+ l + 23,A2ax

(a+a;') --_--'a

----- h 2 1

A+U i

( a+aT~+ 27A,na,,

+ h~ 1 + ~

A+U i ]

<~2(ht - h2) + 23,A2ax

2

•

(37)

for the monotone flux h(., .). Since h I - h 2 ~-1f'(12i)] for the monotone flux schemes considered here [42,43], an approximate estimation of the CFL limit can be obtained from the TVDM requirement (ii) and (37) with h ~ - h 2 replaced by Areax, Y(Ci+l

+

"

2 2

Di+l)~2yAmax + 2T ~4max~ 1

or 2 2

2y 71max+ 27Area x - 1 ~ 0 ,

.

which gives the CFL condition 0 <~ 3,Areax ~<0.3660 . . . . REMARK.

If it is desired to recover the CFL limit to 1/2, the two step method described in

154

K.Y. O~oe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

the previous section can be combined with the flux limiter as follows: Step 1. Flux limiter (12)+ the DFEM solution (25). Step 2. Flux limiter (12)+ the diffus'~t~.a step (26).

For the first step above, it is proved that the scheme is TVDM under the CFL-like condition YAmax ~ 1/2 [35]. In the second step, rewrite the equation in the form

-,,+, --Ui- , _ v A - q+ 1 Ui 2 - * + T31 [ ,'~ 2+1 a + -U,* -- A~Aia*i] __ Ui - u- *i

[c7+ l a + n *

+

(38)

i - D ,*A _ a T ] ,

where u* is the original discontinuous FEM solution defined in (25), ~" = [ AE(uT, *) ~_~-~ U *i-I + A2(u +*) A~

U i'~* ]

~" ! "*2 and C*i + l = !2 TA/2+ 3 0 , D *i = ½7A, 3 0 . Then the TVDM requirement is

~l(Ci* ! + Di*t)=

/2+ "tk" 2~* "

T

I

i+

i]~/-,~

Area x ~l

or • 1 3tAm"x <~ V ~ '

(39)

where A,,,..,~ =

maxJA(u*)l

Vu* ~ (mjn(u?*), max(u ?*)). ! I

Since the first step above is TVDM, it follov's that (m!n(u? *), max(u? *)) c_ (mjn(u,7-"), max(u?")), I

I

i

I

so that A m a x -
1

1

Vmmax ~< 3,A~a~ ~ ~ ~ ~

'

i.e., the second step automatically satisfies the TVDM requirement. Therefore the maximum CFL number 1/2 for the TVDM property is maintained. Since the flux limiter alters the intermediate solution by the order of h 2 [35], the accuracy of the scheme is preserved.

4. Numerical experiments

Three different flux models including a nonconvex one are selected for numerical tests and, for each model, both continuous and discontinuous initial conditions are used to show the

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

155

general performance of the scheme. In all the examples tested, the one step TG-DFEM scheme is directly compared with the DFEM scheme with the trapezoidal quadrature rule for the integration of the time discretization terms and with Simpson's quadrature rule for the other spatial discretization terms. For comparison purposes, in the numerical tests without the flux limiter, the CFL condition of the DFEM y--- CV~--x with C = 1 is used to determine the time step increment. Otherwise, some fixed CFL number near the largest one is used whenever possible. In all the cases where the flux limiter is applied, 0 is fixed to 1 for the consistent (parameter free) application of the method. Since the flux limiter degrades the order of accuracy of the scheme to O(At) near the extremum, the modified flux limiter developed in [44] and [35], ui with

=

min mod(A_a~, ui, A+ ai) =

2 Ki = ~ (3 + 10M) M

f t min mod(A_ai, ffi, A+ff/),

if I ,1 ~< Ki h2 , otherwise,

h2

h "~+

I/t+ail

+ la_a,I

and M = maxlux~(x, 0)l, is used to preserve second order accuracy throughout the solution. One should remember that even though At is chosen to satisfy the linear stability of the DFEM scheme for comparison purposes, it is not reasonable to use such a restrictive time step in practice. 4.1. Linear equation The governing equation (1) with f ( u ) = u is solved with the smooth initial condition 5I

u=

0,

sin(x

-

½)~r+ ½ ,

for 0 ~
and the boundary condition u(0, t ) = 0. The interval is uniformly discretized into 50 elements (Ax = 1 ] 10) and a CFL number y = C ~ with C = 1 is chosen for the time step advance. The exact solution is the initial condition shifted to the right by t. Figure 9 shows the comparison of the DFEM and the TG-DFEM solutions without the flux limiter. It is noted that the first order scheme is still unstable in spite of the restrictive CFL number used and has a wild oscillation at the trailing edge due to the large phase error. The TG-DFEM solution is almost indistinguishable from the exact one except at the leading and trailing edges of the wave and shows a definite improvement in solution quality. To compare the effects of the flux limiters, both Van Leer's flux limiter and the modified flux limiter were applied to the same problem. Figure 10 is the solution with Van Leer's flux limiter applied. Again, the TG-DFEM solution approximates the exact solution much more closely than the first order scheme, but as was pointed out, degradation of accuracy to first order near the extremum due to the flux limiter is noticed. When the modified flux limiter is applied to the TG-DFEM, a result identical to the solution without the flux limiter is obtained, in contrast to the DFEM solution which was improved (Fig. 11). Since the Taylor-Galerkin process is based on the assumption of a continuous solution field, an extreme example for the scheme is a case of a propagating discontinuity. Among such problems, the contact discontinuity is the most difficult one to simulate numerically, because there is no mechanism to recover the discontinuity once it is spread by the numerical diffusion of the scheme. To see the worst performance of the TG-DFEM scheme, the governing

156

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method 1.5

I

U[_~p.I/3

TG-DFEM \ ~

{y=c wi . c-i l(At - 0.03162)

I st Order DFEM

/

.

I.O

0.5

0,0

I

I

I

I

0

I

I

1

I

I

1

Step

Huelber

=

I

I

I

I

2

3

64

Time

Fig. 9. Comparison of the D F E M and T G - D F E M

I

I

I

I

I

I

4

=

0, 2023857E+01

solutions for the linear equation.

Exact Solution

Ist Order DFEM

TG-DrEM Scheme Linear Equation p =I13 Flux Llmlter with 0=I dx = I/I0

0

6

0

4

0

2

I

Step

Hueber

2

-

64

3

T ice

=

4

O. 2 0 2 3 6 6 1 E + 0 1

Fig. lO. Solutions of the linear equation with the Van Leer's flux limiter.

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

/F- \,1 TG-DFEM Scheme

Exact Solution

Ist orderDFEM

157

I I

Linear Equation

l

Modified p=I/3

I

Flux Limiter with 0:I dx : I110

l

I

InitialCondition

1 1 1

\ \ \

\,I/ \

0

I Step

Number

v,

Ix =

I

3

2 fl

6~IL'v4

Time

4 =

0.2023681E+01

Fig. 11. Solutions of the linear equation with the Shu's modified flux limiter.

equation was solved with the discontinuous initial condition u=

1, O,

for O~
and the boundary conditions u(0, t) = 1, u(1, t) = 0 with the interval 0 ~< x ~< 1 divided into 50 elements and At determined by tLe CFL limit y = C ~ to give linear stability to the DFEM scheme. Figure 12 compares the numerical solutions of the two schemes. Clearly the overshoots and the undershoots are greatly reduced using the proposed scheme. When the (modified) flux iimiter is applied, both numerical solutions represents the correct shock position but the TG-DFEM solution is a little bit more smeared than the DFEM solution due to the added diffusion (Fig. 13).

4.2. Inviscid Burgers' s equation For an example of nonlinear problem with a convex flux function, the inviscid Burgers' equation f ( u ) = ½u is tested with several initial conditions. The first problem is Burgers' equation with the smooth initial condition u = sinOrx),

for 0 <~ x ~< 2 ,

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

158 I .5

K

p:l13 No Flux Limiter C=I for ¥ = C ~ dx = 1150

|

TG-DFEM Scheme

I ,0

~ I

1

/

InitialCondition. . ~

Exact Solution

0.5

0.0

I

0.0

0.2 Step

I

0.4

Number

=

I77

I

0.6 Time

-

0.8

1 .0

0,5006312E+00

Fig. 12. Propagation of the contact discontinuity without flux limiter. Ist Order DFEM I

0

TG-DFEM Solution

m

0

8

~olution I

0

6

0

4

m

m

0

2 m

0

I

0

I I

0.0

0.2 Step

Number

0.4 -

175

0.6 T i me

-

t

0.8 0 . 5020458E+00

Fig. 13. Propagation of the contact discontinuity with the modified flux limiter.

1.0

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerk#l finite element method

159

and the boundary conditions u(0, t) = u(2, t) = 0. The exact solution is smooth to t = 1/~r and then starts to form a stationary shock at x = 1. That shock grows gradually to full strength at t = 0 . 5 and then decays out. For discretization of the domain, Ax = 1/20 and A t = 0 . 0 1 ( y ~ X/-A-x) have been chosen for stability of the D F E M scheme. Figure 14 is a comparison to the exact solution of both solutions without the flux limiter at t = 0.3, just before a shock starts to form. Both solutions approximate the exact one very closely, except at the steep wave front where the numerical shock has already started. Nevertheless, the T G - D F E M scheme resulted in a slightly better approximation to the exact ~olution than did the D F E M . Figure 15 shows the solutions without the flux limiter at t = 0.5 when the maximum strength shock has just developed. The T G - D F E M scheme much improved the solution, and is almost identical to the exact one. That is in contrast to the D F E M solution which shows an overshoot near the shock. Also notice that the discontinuity is perfectly captured despite the Taylor-Galerkin procedure applied to the T G - D F E M scheme. This would always be possible whenever the discontinuity of the exact solution is on an element interface. The second problem is to solve the equation with the initial condition u=-2x+

l.,

for 0 ~ < x ~ < l ,

and the boundary condition u(0, t ) = 1, u(1, t) = - 1 . The exact solution is continuous initially but gradually forms a steep transition between the two constant states, that transition later becomes a stationary shock. Figures 16 and 17 show the shock formation process at several time steps in the D F E M and T G - D F E M solution. The space interval is discretized into 50 elements and At is determined by ), = C ~ with C = 1. The T G - D F E M scheme greatly reduced the overshoot and undershoot behind the steep wave front compared to the D F E M solution and shows a sharper shock front at t = 0.5. Notice that the T G - D F E M scheme

-

u

I . 0

TG-DFEM -

DFEM

Solution

5olution

I

-

Exact

Solution

w

I

0,5

I

I

u

0.0

I

-0.5

w

m

-I

.0

--

:

I

I

I

0,0

I

]

I

I

I

0,5

Step

Number

I

I

I

I

I

I

1 .0 =

30

I

I

I

1,5 Time

-

I

I

2,0

0.3000000E+O0

Fig. 14. Solutions without the flux limiter at t = 0.3 just before a shock starts to f o r m at x = |.

K.Y. (?hoe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

160

m

u

TG-DFEM Solution ! .0

DFEM Solution Exact

5olution

O .5

I

0.0

-0

.5 I

m

-I

m

.0 I

-

I

I

I

I

I

0,0

I

I

I

I

I

0.5 Step

I

I

I

I .~

Number

-

50

I

I

I

I

I

I , 5 Time

•

I 2.0

0,4999998E+00

Fig. 15. Solutions without the flux limiter at t = 0.5 for which a full strength shock has just developed at x = 1

!

5

I,_~_.°I__86_ _ ~ I

0

i

ODFEbl solution

[~ = 1/3

I

-0

INo Flux Limiter

I

l~--c~, c = Idx = 1 / 5 0 -I I

I

-!

-., 0.0

I

I 0.2

I

I 0.4

I

I 0.6

I

I

I

0.8

Fig. 16. Stationary shock formation of the D F E M scheme without flux limiter.

I I .0

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method I

161

5

~I

I~o;o.~oo~ 1 182 ~

l

0

0

5

0

0

1~ . . [

B

Initial

t = 0.7006858 1 n= 256

o n d i t i o n ~

I

w

o0

I TG-DFEM solution i.= 1/3 INo Flux Limiter

_-

5

_ -

lT=e~,

-

-!

0

c = 1

I dx = 1/50

lh-.~f

I

m

-I

-

5

I

I

0,0

I

I

0.2

I

0.4

i

I

I

0.6

I

I

0.8

I .0

Fig. 17. Stationary shock formation of the TG-DFEM scheme without flux limiter.

I .0 Ist

Order DFEM

TG-DFEM Solution i

0.5

i

M

m

i

0.0

m

I

Modified Flux Llmlter p - 1/3 7= 112 dx = 1150

i

m

0.5

I

Exact Solution

m

m

I

!.0

I

0.0

I

I

1

0.2 Step

Number

0.4 =

30

I

I

\~--,~ I

0.6 Time

I

I

0.8 -

0.3000000E+O0

Fig. 18. Solutions of the Burgers' equation with the modified flux limiter at t = 0.3.

I ! .0

K.Y. Choe, K.A. Hoisapple, Discontinuous Taylor-Galerkin finite element method

162

captured the discontinuity perfectly without any smearing of the shock front despite the use of Taylor-Galerkin approach. Figure 18 shows a comparison of the two solutions with the modified flux limiter at t = 0.3 for the same problem. The flux limiter did not affect the quality of the final solution (perfect discontinuity) hence that final solution is not shown. As shown in the figure, the T G - D F E M gives a better approximation to the exact solution. The final example for Burgers' equation is a test of the algorithm for a rarefaction wave. This example shows the convergence of the scheme to the entropy solution (expansion wave rather than a rarefaction shock which is also a weak solution of the problem). The initial condition is given by u=

0,

for0~
1,

for0.3~
(40)

and the boundary condition is u(0, t) = 0, u(1, t) = 1. The space interval is discretized into 25 elements and the CFL number 7 = 1 /2 is used with the flux iimiter. In Fig. 19, both solutions approximate the exact solutions closely but, in this special case, the DFEM solution shows a slightly better result.

4.3. Buckley-Leverett equation As the final numerical test for the scheme, the Buckley-Leverett equation 2

U

f(u) = u 2 + (1

u)212

-

E x a c t 5olutlon

o Ist O r d e r D F E M

8

m

F l u x L t m l t e r O .-- I -/=1/2

T_G-DFEM Scheme~

dx = I / 2 5 6

m

4

m

2

m

o

I

I

o.o

JI

I

0.2 Step

Number

0.4 =

213

I

I

I

0.6 Time

I 0.8

=

0.6024597E+00

Fig. 19. Convergence of the solutions to an expansion wave.

I

I ! .0

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

163

0

Ist Order D F E M --

0

B --

0

--

6

Exact

_

No F l u x L i m i t e r

Solution

DFEM S c h e m e

dt = 0 . 0 0 3 3 3 5 __

0

4

0

2

0

0

dx

=

1125

I

I

I 0.0

I

I

0.2 Step

I

I

0.4 Number

-

I

0.6

90

0.8

Time

•

I

I

I .0

0.3001499E÷00

Fig. 20. Solutions of the two schemes without flux limiter at t = 0.3.

I

m

0

Ist Order D F E M

0

B

0

6

Exact Solution at t = 0 3

m

Modified Flux Limiter 0

4

TG-DFEM Scheme

dt = 0.00667 ( y ~ 0.35 ) dx = 1/25 m

P

2

\ 0

I

I

0

I

I

0.2

0.0

Step

I

0

0,4

Number

"

l

qs

6

Time

I

I

0.8

=

0.3001499E+00

Fig. 21. Solutions of the two schemes with the modified flux iimiter at t = 0.3.

I I .0

164

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finiw element method

TG-DFEMScheme ExactSolution/

_

0

'

'

~

.~l_

0.0

I

I

O.Z Step

I

1

O,dl Number

•

75

I

0.6 Time

[

I

0,~ •

I

.0

0.5002499E*00

Fig. 22. Solutions of the two schemes with the modified flux limiter at t = 0.5.

in [29] is chosen for the model equation for the non-convex flux function. The initial condition is u= l / l l ,

forO~
and the boundary conditions are u(0, t)= 1 and u(1, t)= l / l l . In this example, 25 uniform elements are used and At = 0.003335 (~
5. Conclusions

A finite element method based on the Discontinuous Finite Element Method and the Tayior-Galerkin approach is developed. The method demonstrates superior results compared to the Discontinuous Finite Element Method. The major improvements are higher order time

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

165

accuracy and stability under a fixed CFL number up to nearly 0.4, in contrast to the unconditional instability of the DFEM. The scheme preserves the TVDM property under a C F L number up to over 0.366 thus guaranteeing the convergence to a weak solution when combined with the flux limiters described. Furthermore, numerical tests show the convergence of the solutions to the entropy ones without any difficulty, even for nonconvex flux functions in nonlinear problems, and with a correct shock speed. Despite the application of the Taylor series expansion for the time derivative terms, the T G - D F E M does not lose the ability to perfectly capture a stationary discontinuity. An alternative two step method increases the CFL limit up to 1/2. The scheme can be considered a finite element version of Van Leer's higher order extension of the Godunov scheme in the finite difference method, since it takes into account the slopes of the approximate solution projected inside the elements through the Taylor-Galerkin process for the correction of the element interface boundary conditions. But, because this scheme needs only the solutions of the Riemann problem, it is more practical than Van Leer's scheme, which requires the general solution of the nonlinear discontinuous initial value problems at each boundary interface. Unlike the conventional (continuous) finite element method, which always gives an implicit form due to the coupling of the unknowns of neighboring elements, the T G - D F E M is one-step and purely explicit, and therefore does not involve the time consuming iterative procedures which may introduce additional convergence problems during the inversion of the consistent mass matrix. This advantage would become particularly significant in multidimensional problems where the bandwidth of the consistent mass matrix grows rapidly with the spatial dimension. The handling of discontinuities is generally easier in the current approach than in continuous finite element schemes. The investigation of multidimensional cases is in progress and will be presented in subsequent papers [45, 46].

Acknowledgment We thank Professor Bernardo Cockburn for helpful discussions and information.

References [1] P.D. Lax and B. Wendroff, Systems of conservation laws, Comm. Pure Appl. Math. 13 (1960) 217-237. [2] A. Lapidus, A detached shock calculation by second-order finite differences, J. Comput. Phys. 2 (1967) 154-177. [3] R.W. MacCormack and B.S. Baldwin, A numerical method for solving the Navier-Stokes equations with application to shock-boundary layer interactions, AIAA Paper (1975) 75-1. [4] A. Jameson, W. Schmidt and E. Turkel, Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes, AIAA Paper (1981) 81-1259.

[5] R. Courant, E. Isaacson and H. Rees, On the solution of nonlinear hyperbolic differential equations, Comm. Pure Appl. Math. 5 (1952) 243-255. [6] J.T. Steger and R.F. Warming, Flux vector splitting of the inviscid gas-dynamicequations with application to finite-difference methods, J. Comput. Phys. 40 (1981) 263-293. [7] S.K. Godunov, Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics, Math. USSR-Sb. 47 (1959) 271-306. [8] A. Harten, P. Lax and B. Van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev. 25 (1983). [9] B. Van Leer, Towards the ultimate conservation scheme IV: A new approach to numerical convection, J. Comput. Phys. 23 (1977) 276-299.

166

K.Y. Choe, K.A. Hoisapple, Discontinuous Taylor-Galerkin finite element method

[lO] O.C. Zienl~iewicz and R.L. Taylor, The Finite Element Method, 4th Edition, Voi. 1 (McGraw-Hill, London, 1989). [ll] P.M. Gresho, R.L. Lee and R. Sani, Advection-dominated flows with emphasis on the consequences of mass lumping, in: R.H. Gallagher, O.C. Zienkiewicz, J.T. Oden, M. Morandi Cecchi and C. Taylor, eds., Finite Elements in Fluids, Vol. 3 (Wiley, New York, 1978) 335-350. [12] A.J. Baker and M.O. Soliman, Utility of a finite element solution algorithm for initial-value problems, J. Comput. Phys. 32 (1979) 289-324. [13] J.C. Heinrich, P.S. Huyakorn, O.C. Zienkiewicz and A.R. Mitchell, An 'upwind' finite element scheme for two-dimensional convective transport equation, Internat. J. Numer. Methods Engrg. 11 (1977) 131-143. [14] J.C. Heinrich and O.C. Zienkiewicz, The finite element method and 'up-winding' techniques in the numerical solution of convection dominated flow problems, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows, AMD, 34 (ASME, New York, 1979) 105-136. [15] D.F. Griffiths and A.R. Mitchell, On generating upwind finite element methods, in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows, AMD, 34 (ADME, New York, 1979) 91-104. [16] K.W. Morton and A.K. Parrott, Generalized Galerkin methods for first order hyperbolic equations, J. Comput. Phys. 36 (1980) 249-270. [17] T.J.R. Hughes and A. Brooks, A multidimensional upwind scheme with no cross wind diffusion, in: T.J.R. Hugh~:s ed., Finite Element Methods for Convection Dominated Flows, AMD, 34 (ASME, New York, 1979) 19-36. [18] C. Johnson, U. N~ivert and J. Pitk~iranta, Finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 285-312. [191 C. Johnson and J. Saranen, Streamline diffusion methods for incompressible Euler and Navier-Stokes equations, Math. Comput. 47 (1986) 1-18. [20] C. Johnson and A. Szepessy, On the convergence of a finite element method for a non-linear hyperbolic conservation laws, Math. Comput. 49 (1987) 427-444. [21] J. Donea, A Taylor-Galerkin method for convective transport problems, lnternat. J. Numer. Methods Engrg. 20 (1984) 101-119. [22] V. Selmin, J. Donea and L. Quartapelle, Finite element methods for non-linear advection, Comput. Methods Appl. Mech. Engrg. 52 (1985) 817-845. [23] R. L6hner, K. Morgan and O.C. Zienkiewicz, The solution of non-linear hyperbolic equation systems by the finite element method, Internat. J. Numer. Methods Fluids 4 (1984) 1043-1063. [24] W.H. Reed and T.R. Hill, Triangular mesh methods for the neutron transport equation, Los Alamos Scientific Laboratory, Los Alamos, NM, LA-UR-73-479, 1973. [25] W.H. Reed, T.R. Hill, F.W. Brinkley, Jr. and K.D. Lathrop, TRIDENT, a two-dimensional, multigroup, triangular mesh, discrete ordinates, explicit neutron transport code, Los Alamos Scientific Laboratory, Los Alamos, NM, LA-6735-MS, 1977. [26] E Lesaint and EA. Raviart, On a finite element method for solving the neutron transport equations, in: C. de Boor, ed., Mathematical Aspects of Finite Elements in Partial Differential Equations (Academic Press, New York, 1974) 89-123. [27l C. Johnson and J. Pitk/iranta, Convergence of a fully discrete scheme for two-dimensional neutron transport, SIAM J. Numer. Anal. 20 (1983) 951-966. [28] C. Johnson and J. Pitk~iranta, An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation, Math. Comput. 46 (1986) 1-26. [29] G. Chavent and G. Salzano, A finite-element method for the 1-D water flooding problem with gravity, J. Comput. Phys. 45 (1982) 307-344. [30] G. Chavent, G. Cohen, J. Jaffre, M. Dupuy and I. Ribera, Simulation of two-dimensional water flooding by using mixed finite elements, Soc. Pet. Eng. J. 24 (1984) 382-390. [31l A.Y. Le Roux, A numerical conception of entropy for quasi-linear equations, Math. Comput. 31 (1977) 848-872. [32] G. Chavent and B. Cockburn, Consistance et stabilite des schemas LRG pour ies lois de conservation scalaires, INRIA Report no. 710, 1987. [33] G. Chavent and B. Cockburn, The local projection P0 Pl-discontinuous-Galerkin finite element method for scalar conservation laws, IMA Preprint Series # 341, University of Minnesota, 1987, to appear in M2AN. [341 G. Chavent and J. Jaffre, Mathematical Methods and Finite Elements for Reservoir Simulation (NorthHolland, Amsterdam, 1986). [35] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Math. Comput. 52 (1989) 411-435. [36] A. Harten and P.D. Lax, A random choice finite difference scheme for hyperbolic conservation laws, SIAM J. Numer. Anal. 18 (1981) 289.

K.Y. Choe, K.A. Holsapple, Discontinuous Taylor-Galerkin finite element method

167

[37] EL. Roe, Approximate Riemann solvers, parameter vectors and difference schemes, J. Comput. Phys. 43 (1981) 357-372. [38] B. Engquist and S. Osher, Stable and entropy satisfying approximations for transonic flow calculations, Math. Comput. 34 (1980) 45-75. [39] S. Osher, Convergence of generalized MUSCL schemes, SIAM J. Numer. Anal. 22 (1985) 947-961. [401 A. Harten, On a class of high resolution total-variation-stable finite-difference schemes, SIAM J. Numer. Anal. 21 (1984) 1-23. [41] B. Cockburn and C.-W. Shu, The Runge-Kutta local projection Pl-discontinuous-Galerkin finite element method for scalar conservation laws, IMA Preprint Series #388, University of Minnesota, 1988. [42] A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357-393. [43] E. Tadmore, The large time behavior of the scalar, genuinely nonlinear Lax-Friedrichs scheme, Math. Comput. 43 (1984) 353-368. [44] C.-W. Shu, TVB uniformly high-order schemes for conservation laws, Math. Comput. 49 (1987) 105-121. [45] K. Choe and K. Holsapple, The Taylor-Galerkin discontinuous finite element method- An explicit scheme for nonlinear hyperbolic conservation laws, Finite Elements Anal. Des. 10 (1991) 243-265. [46] K. Choe, The discontinuous finite element method with the Taylor-Galerkin approach for nonlinear hyperbolic conservation laws, Ph.D. Thesis, University of Washington, 1991.