The Taylor-Galerkin discontinuous finite element method—An explicit scheme for nonlinear hyperbolic conservation laws

Finite Elements in Analysis and Design 10 (1991) 243-265

243

Elsevier

The Taylor-Galerkin discontinuous finite element method An explicit scheme for nonlinear hyperbolic conservation laws Kyu Y. Choe and Keith A. Holsapple Department of Aeronautics and Astronautics, University of Washington FS-IO, Seattle, WA 98195, USA

Received March 1991 Abstract. A two-dimensional Taylor-Galerkin discontinuous finite element method for the computation of

weak solutions of nonlinear hyperbolic conservation laws is presented. Linear stability is investigated numerically. The resulting scheme is explicit and solves the problem element by element. Extensive numerical tests show second-order convergence of the solution to the entropy one. The numerical solutions imply the total variation diminishing property when the extended flux limiter is applied.

Introduction

Conventional finite element methods use continuous functions for their base and 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. On the other hand, the discontinuous finite element (Galerkin) method, pioneered by Reed and Hill [1], Lesaint and Raviart [2], Johnson and Pitk~iranta [3] and Chavent et al. [4], uses only piecewise continuous functions for base and weighting functions and as a consequence is able to isolate the inter-coupling of the nodal variables within one element. The resulting local mass matrix is easily inverted exactly to give an equivalent explicit expression for the unknowns. The resulting algorithm 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 discontinuous finite element method ( D F E M ) needs to solve only the local Riemann problem to determine the inflow interface (Godunov) values. This makes the scheme practical even for higher order basis functions. Since the discontinuous finite element method is explicit, it could be easily combined with adaptive finite element methods. This scheme has an excellent discontinuity capturing ability due to the discontinuous base functions. The use of the solution to the local Riemann problem makes the scheme appropriate for shock front tracking. Unfortunately, the D F E M is stable only under the restrictive condition C o u r a n t Friedrichs-Lewy (CFL) number 3' ~ O(AvC~--x-).To enforce stability in non-linear problems a flux limiter must be combined with the scheme to give the T V D M (total variation diminishing in the mean) property under a maximum C F L number near 1/2, which guarantees the convergence to a weak solution. But, even with the flux limiter, exentsive numerical experi0168-874X/91/$03.50 © 1991 - Elsevier Science Publishers B.V. All rights reserved

244

K.Y.. Choe, K.A. Holsapple / Taylor- Galerk& discontinuous FEM

Lh I

I

I

I

Xi'2

xi'l

I

1--

~ 1

--I+

H.

I

vi

xi

I

vi÷lt

hi

Xi+l

~

I

xi+2

Fig. 1. Discontinuity of the approximate solution u h.

ments show the convergence of the scheme to a weak, but not necessarily an entropy solution. Also, in non-convex flux problems, the scheme is still anti-diffusive [5]. Recently a R u n g e - K u t t a discontinuous Galerkin method was introduced by Cockburn and Shu to improve linear stability [6] and Choe and Holsapple developed the Taylor-Galerkin discontinuous finite element method [7] for the same purpose. This paper extends the Taylor-Galerkin discontinuous finite element method ( T G - D F E M ) to multi-dimensional problems with a brief description of the fundamental idea introduced in [7].

Theoretical f o u n d a t i o n s in one d i m e n s i o n

The discontinuous finite element method

To illustrate the idea of the Taylor-Galerkin discontinuous finite element method, we first describe the discontinuous Galerkin method for nonlinear scalar conservation laws in one dimension:

ut+ [ f ( U ) ] x = 0 ,

(1)

with appropriate initial and boundary conditions. In contrast to the conventional finite element methods which use C o functions, the discontinuous Galerkin method employs piecewise continuous polynomials for its base functions. Consequently, at each element interface discontinuity of the solution is allowed (see Fig. 1) and the scheme is capable of capturing a shock perfectly without spurious oscillation when it falls on the element boundary. On the other hand the solution is not defined uniquely at each nodal point. Therefore, some appropriate technique to define the interface value is required. Chavent et al. introduced the Godunov value [4] which is essentially the exact solution of the Riemann problem constructed with the left and right limiting values of the numerical solution at the element interface. The use of the Godunov value at each element interface yields an upwind effect [8] to the numerical scheme. In problems where the exact solution cannot be easily determined, simple approximate Riemann solvers such as those of Roe [9] or of Engquist-Osher [10] can be used to determine an alternative to the Godunov value.

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

245

The discontinuous Galerkin method employs the forward Euler time stepping for the time discretization of the unsteady term to yield the weighted residual statement n+l

n

(.,u~ --uhv dx + f~,+,Of(u~) Jx, -At x, --Vox dx + {f(u +') -f(~)}v(xi)

- {f(ui-+t ) -f(,~+m)}V(X,+l) = 0,

(2)

where uT, = Uh(t,) and t o = 0, t l, t 2. . . . . t N = t is the partition of the finite time interval [0, t]. The resulting discrete equation is very simple to code and purely explicit. The problem can be solved element by element and does not require any matrix inversion process which can introduce convergence problems. However the scheme has a restrictive stability limit [11] and fails to converge to the entropy solution for nonconvex flux functions [5].

Taylor-Galerkin discontinuous finite element formulation The instability and the non-convergence of the discontinuous Galerkin method necessitated addition of artificial diffusion to the scheme. Any popular technique such as Lapidus, Jameson or MacCormack damping in finite difference methods or the streamline upwind Petrov-Galerkin method [12] or Taylor-Galerkin method [13] in finite element methods could be used for this purpose. To obtain higher-order time accuracy, here we will employ the Taylor-Galerkin method for linear stability. Near a shock in nonlinear problems, the flux limiter is the most appropriate and natural way for the discontinuous Galerkin method to add further numerical dissipation. To introduce Lax-WendroWs idea, take the forward Taylor series expansion of the temporal derivative term about time level n. Replacing the higher-order time derivative terms with the spatial ones using the governing equation, we approximate the governing equation (1) by

u "+l - u "

At 0 I[A 20u ~ n

At

2 Ox~

~x)

Of + O-x=0'

(3)

up to second order. To avoid losing the second-order derivative terms using piecewise linear discontinuous base functions, the additional diffusion term is replaced by the mixed variable defined by At 0u r = - --A z2 0x"

(4)

With continuous piecewise linear functions for r and the trapezoidal quadrature rule for integration, the weighted residual form of the Taylor-Galerkin discontinuous finite element method is written as fxi+lU~+l

x,

-At

U hn

fXi+lar n

v dx + ix,

O---~-vd x -

+f(¢n+l)U(Xi+l) --f(~?)V(Xi) = O,

rXi+l

Jx,

Ov

f(Uh)~xdXn (5)

with At

fb

2 ()U

fabrWdx=-5 Joa wdx

(6)

where a and b are the boundaries of the domain. Due to the higher-order accuracy of the Taylor-Galerkin approach, the resulting scheme is formally uniformly second-order accurate in both time and space.

K Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM

246

A Fourier series analysis gives the theoretical maximum CFL number (y) of 0.392647 with the one-step method and 1 / 2 with a two-step alternative. Also, TVD stability is guaranteed under the CFL condition 3' ~ 0.3660 with the one-step method and 3' ~< 1 / 2 with the two-step method when the flux limiter is applied. For the details of the theoretical analysis of the Taylor-Galerkin discontinuous finite element method, refer to [7].

TG-DFEM in multi-dimensions

Mathematical formulation Hyperbolic conservation laws in multi-dimensions for a scalar quantity u(x, t) are written

ut+V'f(u

)=0

in $2,

(7)

with the boundary condition

u(x,t)]rm=Ur(S,t)

on Fin,

(8)

and the initial condition

u ( x , O) =Uo(X )

in g2,

(9)

where g2 is the spatial domain of interest, F its boundary, Fin the inflow boundary for which the inequality n 'A < 0

(10)

holds, s the arc length along Fi, from a fixed point, n is the outward normal to the boundary and A = f ' = d f / d u is the flux gradient. To apply the Taylor-Galerkin methodology, eqn. (7) is rewritten in the form

u, = - V . f ( u )

= - A .Vu

(11)

and differentiated with respect to time to obtain

utt = - [ V ' f ] t

= - V . [ f t ] = - V . [Aut] = V . [ A ( u ) A ( u ) . V u ] .

(12)

The Taylor series expansion of the temporal derivative term is U n+l

ut -

--U

At

n

At 2

ut~t + O ( A t 2 ) .

(13)

Therefore, up to second-order accuracy in time, eqn. (7) is approximated by

U n+l __ Un At

At + 17 . f " - - - V . 2

[anA n ' V u ]

= 0,

(14)

where abbreviations u n = u(x, t,), f " = f ( u " ) and A" = d f ( u " ) / d u = A(u") are used. The weighted residual statement for the solution of eqn. (14) with (8) is to find a function u n+l in 5 ~ such that for all v ~ ~" U n+l

ffl

--

At

Un

At ud~'-~q'-f~2~'fnl) d~'-~---ff2~'[AnAn'~un]ud~'~2

+ f_ [ f ( u ~ , ) - f ( u " ) ] -/

in

.nv OS = O;

247

K . I ( Choe, K.A. Holsapple / T a y l o r - G a l e r k i n discontinuous F E M

or, using integration by parts, U n+l _ U n

fo

At

v dO-fs

n.Vv d O - At 2 f V. [A~A n.vu~]v dO

+ frif(u~)_ "nv d S + fv,,.f(un)"nv d S = 0

Vv ~ ~',

(15)

where ~ ' c HI is the space of base functions, ~ c H 1 is the space of weighting functions and H 1 is the Sobolev space of functions with square-integrable first derivatives and u r is the prescribed boundary value on the inflow boundary Fin.

Finite element approximation To approximate the solution of eqn. (7), we introduce the space of discontinuous piecewise polynomials of degree k > 0,

M~={v~L2(g2)'vlr~Pk(T)

VT~Ch},

(16)

for the finite element space where C h = {T} is a finite element partitioning of the domain /2 indexed by the parameter h representing the maximum diameter of T ~ C h and Pc(T) is the space of polynomials of degree ~
(17)

and approximate equation (15) by the discontinuous mixed finite element formulation: Find a function u n÷l ~M~ such that u n + 1 __ l.ln

fT

?,7 vdD-f/"'VvdD-fr(V'd")vdg2+fa/(,)'nvdS=O (18)

for all v ~ Mg where d~' ~ Wh1 (i = 1, 2) is defined by

fadnw dO = ½AtfAn(V'f)w dO =½Atfa (Az4 ~'vun)w

dO

Vw~ W~h*, (19)

where Wh~ is the space of continuous piecewise linear polynomial functions of degree k, i.e.,

W~={w~Hl(O):Wlr~&(T)

VT~Ch},

(20)

and sc is the Godunov value defined on the element interface ~T similar to the one-dimensional case. Denoting bases of the function spaces M 2 and W~ as {N/} and {wi} respectively (see Fig. 2), we approximate the solutions by u " + 1 = Y'.uT+ 1Ns,

v =Nj,

(21)

d = ( d l, d2),

d k = £d~w i,

w = Wi

(22)

/(Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

248

D Fig. 2. Discontinuous and continuous piecewise linear functions t, e M~I and w e W~ for basis functions in the mixed formulation.

Then we obtain from (18)

<'+'

.I

0~

ON,

fr " AT UTNiN, ds2- jri fl-~x +f277-y ) + ~T(fi(~)n

./d, 0N,

i ON•

d.Q (23)

, + f 2 ( ~ ) n i ) N j dS = 0

and from (19)

. r

aN,

oN,. 1

dl fowiwi d O = }Atuij'ail A~(u)"Tx + A,(u) A=(u)v iw, d,O, d~2fwiwj dO = J S~

½Atu i

fo[ A l ( U ) A 2 ( u

)

~Ni +AZ(u) --£--,Wj oy j dO

Vw~ e w~. (24)

wi

Note that N, = in an element T. Equation (23) contains the integration of Godunov fluxes fi(sc) along the element interface. In contrast to the one-dimensional case, the general solution to the two-dimensional Riemann problem is not available. Therefore we approximate the local solution of the two-dimensional Riemann problem at each point on the interface by that of the corresponding one-dimensional problem in the direction normal to the boundary. Thus for the approximation of the boundary integral term, numerical quadrature rules such as Gaussian integration should be employed. Then it is convenient to use the one-dimensional Riemann problem solution at each integration point for the approximation of the solution to the two-dimensional Riemann problem (see Fig. 3).

Flux limiter extension

It is proved in one dimension that the T G - D F E M combined with van Leer's flux limiter [14] is total variation diminishing (TVD) [7]. The consequence of the TVD property is the convergence of the numerical solution to a weak solution even for nonlinear problems, and non-oscillatory solutions near a shock. There are many ways to define the flux limiter in two dimensions, depending on the criterion introduced to define the flux-limited solution. One example is the direct extension of

K.Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM

// /

,

,-

- . . \ - \

/ / . /

/

249

//

n~ \

I I\

ff : Outward Normal

Fig. 3. Integrationpoints for the approximation of the generalizedRiemann problem.

van Leer's flux limiter with distance from the original solution for the criterion. Another example is the one introduced by Cockburn et al. [15] using the maximum principle for the criterion. In this paper, the direct extension of van Leer's flux limiter is used for its simplicity. The idea is to find a function as close as possible to the original solution derived by the algorithm described in the previous sections, under some restricton, by minimizing some measure of deviation between them. The restrictions are directly extended from those of the one-dimensional van Leer's flux limiter, i.e., the numerical solution on an element should lie between the minimum and the maximum of the average solution values on the neighboring elements. The 2-norm of the difference vector f i - u is introduced as the measure of the difference to the original solution. Then the problem reduces to one of quadratic programming which is convex and separable. Dualizing the problem, we can find a unique analytical solution. For the details of the flux limiter, see Chavent and Jaffre [5] or Choe [16]. It has not been proved that the discontinuous Galerkin method with the flux limiter is T V D in multi-dimensions. Nevertheless, the numerical tests here strongly suggest that the T G - D F E M with flux limiter is T V D in two dimensions. In the presence of a strong shock, the flux limiter implemented in this p a p e r did capture the shock without spurious oscillation. Unfortunately, any T V D scheme in more than one dimension is at most first-order accurate [17].

Numerical stability considerations An attempt to evaluate the stability limit of the proposed scheme in two dimensions has been made. Again, the Fourier method is used for the analysis of linear stability. Consider the discrete version of eqns. (23) and (24) with a linear flux vector function fi(u) = aeu, where ae are constants (i = 1, 2) for the Taylor-Galerkin discontinuous finite element method in the rectangular grid of element size Ax X Ay. Take a typical term in the Fourier series expansion of the solution

u( x, y, t) = ~ Zmn( t ) e i~m* e i~,,y, n,m

where Zmn(t) is the amplification coefficient function of time, i = ~/- 1 is the unit imaginary

250

K.Y. Choe, K.A. Holsapple / T a y l o r - Galerkin discontinuous F E M

number, m, n are integers, /3m = m ' r r / L x, [3n = n r r / L `' are the wave numbers in each direction and L , and L,. arc the lengths of the rectangular domain in the x- and y-direction respectively. Then Hi ± l,j + m = e + ilOI e + im°2tli'J k ~'k ,

for l, m = 0, 1 where 0~ =/3,, Ax and 0 2 = ~ Ay. Using this in the discrete equation, after some rearrangement we obtain the following set of algebraic equations

{, ~. } = [G(3/1, 3/2, 01, 02)]{u~:J}, (25) where G(V~, 3'2, 0~, 0=) is the 4 × 4 magnification matrix with 4 parameters, t,~'j is the updated solution at time level n + 1 at node k of the element Ti, i and y~ = a~ A t / A x , Y2 = a 2 A t / A y are the CFL numbers in the x- and y-directions respectively. The stability requirement of the numerical scheme is [I G II ~ 1 + c At, where C is a constant independent of Ax, Ay, At and O's, and I1" II is a norm of a matrix. Since II G II is a function of Yl, Y2, Oi and 02, IJGII 4:p(G) in general where p(G) = max, t Ai] is the spectral radius and A/(TI , 3/2, 01, 02) is an eigenvalue of the magnification matrix G. But the equality holds only if A's are distinct from each other. Following the one-dimensional stability analysis, we need an explicit expression for A, to study further. The characteristic equation of a 4 × 4 matrix is a polynomial equation of degree 4. It is always possible to solve a polynomial equation of degree 4 or less. But in practice, it does not appear to be possible to obtain an explicit solution of this problem due to the complexity of the expression for each element of the magnification matrix G. Thus, instead of the analytical study of linear stability, a numerical estimation is used as an alternative. The numerical result is given in Fig. 4. One can easily find that this numerical evaluation of linear stability gives an overestimate to the

(/.5

1.0

"""

""

""" " "

'"

""

0.4

Yy o.3

0.2

0.1 1.0

0.0 0.0

0. l

0.2

m 0.3

r 0.4

0.5

0.6

Tx Fig. 4. M a x i m u m absolute e i g e n v a l u e s of the m a g n i f i c a t i o n matrix G as a function of y ' s - - C o n t o u r plot.

K.Y. Choe, K.A. Holsapple / Taylor - Galerkin discontinuous FEM

251

CFL limit by observing that at y~ = 0 or Y2 = 0, which includes the one-dimensional case, the two-dimensional CFL limit cannot be larger than that of the corresponding one-dimensional problem which has been analyzed in [7]. This overestimate seems to be caused by the numerical errors introduced during the eigenvalue calculations. Further numerical experiments show that 0.31 is the maximum CFL number for linear stability of two-dimensional problems. Even though this approach gives a rough numerical estimation for linear stability, it reveals at least one interesting behavior of the CFL limit in two dimensions as a function of Yl and Y2. In contrast to the linear stability result for the conventional Taylor-Galerkin finite element method presented by Selmin et al. [18], where it is shown that the CFL limit deteriorates considerably as the flow direction deviates from one of the grid axes, the stability limit for the suggested scheme in each direction is preserved almost independent of the flow direction. The T V D stability analysis for the discontinuous finite element method in two dimensions is not available yet to the author's knowledge. Therefore the T V D stability analysis for the two-dimensional T a y l o r - G a l e r k i n discontinuous finite element method is left unsolved along with the analytical solution of yon Neumann analysis for the linear stability.

Numerical tests

Various examples of flux functions ranging from linear to nonlinear fluxes have been selected and tested. Four-node bilinear isoparametric quadrilateral elements were used in all cases. These examples were computed on the San Diego Supercomputer Center C R A Y Y-MP 8/864. The results were directly compared to those of the original discontinuous finite element method on the same discretization. Exact solutions were used for the computation of the numerical errors whenever possible.

Linear hyperbolic equations Linear equations are the most interesting examples for the numerical tests because they have a direct relationship to the linear stability analysis and easily demonstrate the effectiveness of a numerical method. In the following, linear hyperbolic conservation laws,

u, + A~( x )ux + A2( x)u~ = O, are investigated with various flow directions and initial conditions.

Translating sine wave R e m e m b e r i n g that the T a y l o r - G a l e r k i n procedure is based upon the existence of the higher-order derivatives in space, we can expect this scheme to work best for a smooth solution field. To see the performance of the suggested method, a translating wall of sine wave profile is chosen for the test problem. The flow direction is normal to the initial solution contour line. To study the mesh effect on the numerical solutions, various flow directions including 0 = 0 ° , arctan 1/2, arctan 2 / 3 and 45 ° are investigated with the flow velocity A 1 = c o s 0, ~z~ 2 = sin 0 in a uniformly discretized (Ax = Ay) rectangular domain. To exaggerate the differences of the results, a rather coarse mesh size of Ax = 0.2 (10 elements for the whole sine wave length) is chosen. Periodic boundary conditions are used in the flow direction. To minimize the effect of the boundary condition on the solution algorithm, the solution is assumed to be a portion of the solution in an infinite domain cut along the boundary. Even though the numerical estimation of the two-dimensional linear stability limit (see Section on Numerical Stability Considerations) predicts the maximum CFL number of over 0.4, the numerical tests showed that a CFL number of about 0.3 is the practical limit for

252

K.Y. Choe, K.A. Holsapple / Taylor -- Galerkin discontinuous FEM

0

.g

B

?

¢.,

,,.J e-

? ~7 o

-,4

E 0 i1--,

a

,_

o

©

.g

[-,

Ii

....~

[,A ~Z

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

253

the extreme case of 7x = 7y- At this CFL number, all solutions were stable. For the sake of safety, a CFL number of 0.27 is chosen for the rest of the numerical simulations, including the nonlinear ones, unless indicated otherwise. Among the numerical tests, the case of Cx = Cy (45 ° flow angle) is the least diffusive and gives the best solution whereas it gives the most restrictive CFL condition (7 ~ 0.3). The T G - D F E M solution at t = 2.0 is shown in Fig. 5 with its profile compared with the exact solution. The D F E M solutions are not shown since all of them were unstable at this C F L number.

Rotating Gaussian hills A standard test problem for a linear convection equation is the Gaussian hills in a rotating flow field. This problem is insensitive to the mesh orientation. The problem is described in Fig. 6 along with the initial condition. Note that the exact solutions at time t = nw coincides with the initial condition (Fig. 7). The maximum CFL number is restricted to 0.35.

Y 7-

[] ~

VI=-(y-y0) V2

=

x - x0

~ CenterofRotation 1,i5 FlowDirection ~ yo) ~ d

FlowDirection 1.,

\

Af

.......

~

q

A7

,-

~i I

1

//~

0.8 0.20.40"6/

i ~

-d 2

U=exp[-~5-] 135)

115 2.0 2.5 3 3.5 4 4.5 SectionA-A Fig.6.InitialconditionandproblemdescriptionforthetwinGaussianhillsinarotatingflowfield.

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

254

Fig. 7. Exact solution at t = nnz; n = 0 is the initial condition. In Fig. 8 numerical results of D F E M and T G - D F E M at t = v are compared side by side. Solutions are shown from the top both without and with the flux limiter. Solutions of the D F E M are shown on the left and T G - D F E M results are shown on the right. The numerical solution of the D F E M without the flux limiter is still unstable. The results at later times are shown in Fig. 9. The amplitude of the Gaussian hills became much larger than 1, indicating instability. Behind the Gaussian hill, following wakes are developed and will grow due to the instability. We notice that the T G - D F E M improved the result considerably. But for this smooth solution case, the use of the flux limiter just distorted the numerical solution of T G - D F E M without any advantage, whereas it stabilized the numerical solution of D F E M at the price of severe distortion. We conclude that for a smooth solution field, the flux limiter is not recommended for the T G - D F E M scheme. The order of convergences of L~, L 2 and L~ error norms are plotted together in Fig. 10 for the T G - D F E M results at t = "rr. The order of convergence rapidly converges to 2 as the mesh is refined showing that the T G - D F E M is a second-order scheme.

Burgers' equation As a classical example of nonlinear convex flux functions, the inviscid Burgers' equation

is tested. We have seen that the scheme captures a stationary shock perfectly in one dimension. Therefore it is interesting to see how it works in two dimensions with different mesh effects.

Sine wave initial condition The first example for Burgers' equation is the retest of the one-dimensional problem with sine wave of period 2 in two dimensions. As in the linear case, a CFL number 0.27 is used. In Fig. 11, numerical solutions with two different element sizes are shown at t = 0.5. The

K.Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM

255

6[.-., z~ k~

E~

h.

E ~ E

,2

Z k~

K.Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM

256

i

3

g)

f, c

3-

Fig. 9. N u m e r i c a l solutions for the twin G a u s s i a n hills in a r o t a t i n g flow field at t i m e t = 5.005.

one-dimensional profiles of the solutions at t = 0.5 are compared with the exact solution, within the resolution of the numerical solution. In the coarse mesh (20 elements along the x-axis), the solution deviates slightly from the exact one near the shock and at the right boundary. The inaccuracy of the solution near the right boundary is due to the rough treatment of the periodic boundary condition which caused the change in the periodic length by Ax. All the errors vanished in the fine mesh (40 elements along the x-axis) and the numerical solution is almost identical to the exact one. As in the one-dimensional case, the method captured the stationary discontinuity perfectly. The same problem has been computed with a rotated axis. The angle between the contour line and x-axis is 45 ° in a counterclockwise direction. Figure 12 demonstrates how the solutions at t = 0.5 converge to the exact one. The mesh effect resulted in a smeared shock front of the numerical solution, but spurious oscillations near the shock were not observed. The order of convergence is plotted in Figs. 13 and 14 for the fixed time t = 0.5 when the shock strength has reached its peak. Despite the existence of the strong shock, the conver-

K.Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM 2

T

--

I

I

257

I

1.5

g

///~/

1

/

13

L2 Norm

- - ~ - - L~ Norrn

/ I

0.5

15- 30 elemts

I

I

I

30-60 elemts 60-120 elemts No. of Elements

120-240 elemts

Fig, 10. Order of convergence of the T G - D F E M scheme for the rotating flow. The convergence rate increases uniformly toward 2.0

lllllllllllllIlllllJ J

s

,/

8.

0 ©

t.i

Fig. 11. T G - D F E M solutions for Burgers' equation with sine wave initial condition (contour lines are parallel to the y-axis). From the left, solutions with 20 elements and solutions with 40 elements are shown. In the second row the solution profiles are compared with the exact solution at t = 0.5.

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

258

O.D

~0

LO

O0

l0

LO

Fig. 12. T G - D F E M solutions for Burgers' equation with 45 ° sine wave initial condition. From the left, solutions with 20 elements and solutions with 40 elements are shown. The solution profiles are compared with the exact solution at t = 0.5. 4 :

0 I ~

//

i

\.

i

L I Order L2 O r d e r

L ~

I -

72

0

.......................

!../

........................

....................

..................

......................... -I

I

10 -

. . . . . . . . . .

I 20

20

- 40

.....

i I 40

- 80

80

__

- 160

No. Elements Fig. 13. Order of convergence for Burgers' equation. The shock line is parallel to the y-axis.

K.Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM

259

1.5

1

...................

o

L 1 Order

[]

L20rOer

.................. i............................. ~

--o-- L~ Order

"i

J~

'

!

.............................. :[~ :

r-,

e, =~ 0.5 o

~5

0

-0.5

..............................................

-1

I

10 - 2 0

~ . . . . . . Z _._.._

I

20 - 40

40 - 80

No. Elements

Fig. 14. Order of convergence for Burgers' equation. The shock angle is 45 o to the y-axis. gence test shows the scheme is still second-order accurate when the shock is aligned with the element interfaces. On the other hand, when the shock line is oblique to the boundaries, the order of convergence drops to 1. This confirms the second-order accuracy of the scheme for a smooth solution but possible deteriorated convergence under the presence of a shock misaligned to the mesh. Adaptive finite element methods or front tracking methods could make the scheme uniformly second order irrespective of the presence of the shock. Propagating shock wave

As the second example for Burgers' equation, a propagating discontinuity with the Dirichlet boundary condition is simulated. In this case, the solution surface is also discontinuous along the sides as shown in Fig. 15. The solution surface of the test is extremely distorted both at the front and the sides of the shock. Therefore we can easily observe the cross-wind diffusion effect near the sides of the shock, as well as the spurious oscillation near the shock front at the same time. Any numerical scheme would have difficulty in reproducing accurate results of this kind of pedagogical problem. The D F E M and T G - D F E M solutions are compared in Fig. 16. Apparently the T G - D F E M scheme reduced the magnitude of the oscillations which appear in the solutions of the D F E M . In Fig. 17, the flux limiter was applied. The flux limiter slightly smeared the shock front in the T G - D F E M solution compared to the D F E M solution. Nevertheless, the shock position is very accurate for all cases. The flux limiter eliminated the oscillations near the shock showing that the resulting scheme is virtually TVD. Rarefaction wave

The same problem has been computed with the inverse initial condition. In this case the initial discontinuity must be smoothed out for an entropy solution as time passes. Figure 18 shows the initial condition and the exact solution at t = 2.0.

260

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

cO e-

8 e-

©

K.Y. Choe, K.A. Holsapple / Taylor - Galerkin discontinuous F E M

~i

///VV V / l / Y / ~

{

261

IIIVVVVI~I~

Exact

Exact

o

!

ko

i o.o

1',~

~OW A ~ S

~OW-AII$

Fig. 16. D F E M solution and T G - D F E M solution with the Dirichlet boundary condition. Shock angle = 45 °, t = 3.0

The numerical solutions are presented and compared in Fig. 19. A "wiggle" near the stationary point is noticed in the D F E M solution. This is due to the numerical oscillation introduced by the jump discontinuity at the very beginning of the computation. Also the D F E M solution has adverse errors along the smooth zone, but bends back to give a relatively accurate wave front position. On the other hand, the T G - D F E M solution does not show the wiggles near the stationary point and is more accurate except at the wave front. The wave front in the T G - D F E M solution is rather smeared out due to the added damping. The flux limiter eliminated the wiggles in the D F E M solution completely in Figure 20. But it aggravated the solution accuracy elsewhere. The T G - D F E M solution does not show much change but became slightly better than without the flux limiter.

Conclusions

A multi-dimensional extension of the Taylor-Galerkin discontinuous finite element method is studied. The method has been developed by the authors [7,16] based on the discontinuous finite element (Galerkin) method and the Taylor-Galerkin method. The scheme demon-

K.Y. Choe, K.A. Holsapple

262

/ Taylor- Galerkin discontinuous FEM

o

O-DFEM Exact

i

:1 U

U

1.o

3.O

Z°Flow - Axis

1,0

~.e

0.0

10

2.0

3.0

I0

Flow - Axis

Fig. 17. Flux limited solutions with the D i r i c h l e t b o u n d a r y condition. Shock angle = 45 o, t = 3.0.

8

? ,..< ~ "

D 5 X -A~<]S

~

Fig. 18. Initial c o n d i t i o n and exact solution at t = 2.0 for B u r g e r s ' e q u a t i o n (rarefaction case). T h e r e c t a n g u l a r d o m a i n is discretized into 25 × 25 elements.

K.Y. Choe, K.A. Holsapple / Taylor- Galerkin discontinuous FEM

263

60

Exact O.DFEM ~ : ~ - . . . -

8-

U

Exact

8

F E M ~ TG-D

U d

o.o

Io

zo

~o

Flow - Axis

40

~o

oo

1o

~o

3.0

Flow - Axis

4.0

~o

Fig. 19. D F E M solution (left) and T G - D F E M solution (right) for Burgers' equation (rarefaction case) without flux limiter. The solution profiles are compared with the exact solution at t = 2.0, CFL = 0.27.

strates superior results compared to the discontinuous finite element method, with higherorder time accuracy and better stability. Numerical tests confirm that the scheme is second order accurate for a smooth solution without a discontinuity. Near a discontinuity, the order of convergence dropped down to 1. In two-dimensions, numerical estimation shows the maximum CFL number for stability is approximately 0.3 which is comparable to the LaxWendroff scheme and much better than the discontinuous finite element method which requires that CFL ~ ~ for stability. The scheme has negligible cross-wind diffusion effect, and thus captures side discontinuities crisply irrespective of the mesh orientation. It is easily combined with a flux limiter, resulting in a superior discontinuity capturing capability even for a moving shock. Unlike the conventional (continuous) finite element method, which always gives an implicit form due to the coupling of the unknowns of neighboring elements, the TG-DFEM is one-step and purely explicit, and therefore does not involve time consuming iterative procedures which may introduce additional convergence problems during the inversion of the consistent mass matrix. This advantage may become particularly significant in multi-dimensional problems where the band-width of the consistent mass matrix grows rapidly with the spatial dimension. The handling of discontinuities and the implementation of the adaptive

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

264

0t

!I/ i

~J

~o

lo

zo

30 Flow - Axis

o

~o

o.o

[!xact

i.o

zo

30

40

50

Flow ~ Axis

Fig. 20. D F E M solution (left) and T G - D F E M solution (right) for Burgers' equation (rarefaction case) with flux limiter. The solution profiles are compared with the exact solution at t = 2.0, CFL = 0.27.

finite element method is easier in the current approach than in continuous finite element schemes. Despite the application of the Taylor series expansion for the time derivative terms, the TG-DFEM does not lose the ability to perfectly capture a stationary discontinuity when the element interfaces are aligned with the discontinuity. This suggests the introduction of the front tracking method [19] or moving finite element method [20] to the present scheme in future research for perfect capturing of a moving shock.

References [1] W.H. REED and T.R.HILL, Triangular mesh methods for the neutron transport equation, Rept. No. LA-UR-73479, Los Alamos Scientific Laboratory, Los Alamos, NM, 1973. [2] P. LESAINT, and P.A. RAVIART, "On a finite element method for solving the neutron transport equations", in: Mathematical Aspects of Finite Elements in Partial Differential Equations, edited by C. DE BOOR, Academic Press, New York, pp. 89-123, 1974. [3] C. JOHNSON, and J. PITK~RANTA, "An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation", Math. Comput. 46, pp. 1-26, 1986.

K.Y. Choe, K.A. Holsapple / Taylor-Galerkin discontinuous FEM

265

[4] G. CHAVENT, and G. SALZANO,"A finite-element method for the 1-D water flooding problem with gravity", J. Comput. Phys. 45, pp. 307-344, 1982. [5] G. CHAVENT,and J. JAFFRE, Mathematical Methods and Finite Elements for Reservoir Simulation, North-Holland, Amsterdam, 1986. [6] B. COCrmURN, and C.-W. SHu, "TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework", Math. Comput. 52 (186), pp. 411-435, 1989. [7] K.Y. CHOE, and K.A. HOLSAPPLE, "The discontinuous finite element method with the Taylor-Galerkin approach for nonlinear hyperbolic conservation laws", Comput. Methods Appl. Mech. Eng. 94, to appear, 1992. [8] S.K. GODUNOV, "Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics", Mat. Sb. 47 (3), pp. 271-306, 1959. [9] P.L. ROE, "Approximate Riemann solvers, parameter vectors and difference schemes", J. Comput. Phys. 43, pp. 357-372, 1981. [10] B. ENGOUIST, and S. OSHER, "Stable and entropy satisfying approximations for transonic flow calculations", Math. Comput. 34, pp. 45-75, 1980. [11] G. CHAVENT, and B. COCKBURN, Consistance et stabilit~ des sch6mas LRG pour les lois de conservation scalaires, Rept. No. 710, INRIA, France, 1987. [12] T.J.R. HUGHES, and A. BROOKS, "A multidimensional upwind scheme with no cross wind diffusion", in: Finite Element Methods for Convection Dominated Flows, edited by T.J.R. HUC;HES, ASME, AMD-Vol. 34, New York, pp. 19-36, December 1979. [13] J. DONEA, "A Taylor-Galerkin method for convective transport problems", Int. J. Numer. Methods Eng. 20, pp. 101-119, 1984. [14] B. VAN LEER, "Towards the ultimate conservation scheme IV: A new approach to numerical convection", J. Comput. Phys. 23 (3), pp. 276-299, 1977. [15] B. COCKaURN, S. Hou and C.-W. SHU, "The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case", Math. Comput., to appear. [16] K.Y. CHOE, The discontinuous finite element method with the Taylor-Galerkin approach for nonlinear hyperbolic conservation laws, Thesis, University of Washington, 1991. [17] J. GOODMAN and R. LE VEQUE, "On the accuracy of stable schemes for 2D scalar conservation laws", Math. Comput. 45, pp. 15-21, 1985. [18] V. SELMIN, J. DONEA and L. QUARTAPELLE, "Finite element methods for nonlinear advection", Comput. Methods Appl. Mech. Eng. 52, pp. 817-845, 1985. [19] J. GLIMM, C. KLINGENBERG, O. McBRYAN, B. PLOHR, D. SHARP and S. YANIV, "Front tracking and two-dimensional Riemann problems", Adv. Appl. Math. 6, pp. 259-290, 1985. [20] K. MILLER,"Recent results on finite element methods with moving nodes", in: Accuracy Estimates and Adaptive Refinements in Finite Element Computations, edited by I. BAnUSKA, O.C. ZIENKIEWlCZ,J. GAGO and E.R. de A. OLIVEIRA, Wiley, New York, 1986.