Hyperbolic Problems

Back

____ Conservation Laws

We have successfully applied _nite element methods to elliptic and parabolic problems_

however_ hyperbolic problems will prove to be more di_cult_ We got an inkling of this

while studying convection_di_usion problems in Section __      _ Conventional Galerkin meth_

ods required the mesh spacing h to be on the order of the di_usivity _ to avoid spurious

oscillations_ The convection_di_usion equation __       ___ changes type from parabolic to hy_

perbolic in the limit as _ _ _ The boundary layer also leads to a jump discontinuity in

this limit_ Thus_ a vanishingly small mesh spacing will be required to avoid oscillations_

at least when discontinuities are present_ We_ll need to overcome this limitation for _nite

element methods to be successful with hyperbolic problems_

Instead of the customary second_order scalar di_erential equation_ let us consider

hyperbolic problems as _rst_order vector systems_ Let us con_ne our attention to con_

servation laws in one space dimension which typically have the form

ut _ f u_x _ bx_ t_ u__ _____a_

where

ux_ t_ _

_____

u_x_ t_

u_x_ t_

___

umx_ t_

_____

_ f u_ _

_____

f_u_

f_u_

___

fmu_

_____

_ bx_ t_ u_ _

_____

b_x_ t_ u_

b_x_ t_ u_

___

bmx_ t_ u_

_____

_____b_

are m_dimensional density_ _ux_ and load vectors_ respectively_ It_s also convenient to

write _____a_ as

ut _ Au_ux _ bx_ t_ u_ _____a_

_

_ Hyperbolic Problems

where the system Jacobian is the m _ m matrix

Au_ _ f u_u_ _____b_

Equation _____a_ is called the conservative form and _____a_ is called the convective

form of the partial di_erential system_

Conditions under which ______ and ______ are of hyperbolic type follow_

De_nition _______ If A has m real and distinct eigenvalues __ _ __ _ _ _ _ _ _m and_

hence_ m linearly independent eigenvectors p____ p____ _ _ _ _ p_m__ then _____a_ is said to

be hyperbolic_

Physical problems where dissipative e_ects can be neglected often lead to hyperbolic

systems_ Areas where these arise include acoustics_ dynamic elasticity_ electromagnetics_

and gas dynamics_ Here are some examples_

Example _______ The Euler equations for one_dimensional compressible inviscid _ows

satisfy

_t _ mx _ _ _____a_

mt _

m_

_

_ p_x _ _ _____b_

et _ _e _ p_

m

_

_x _ _ _____c_

Here __ m_ e_ and p are_ respectively_ the _uid_s density_ momentum_ internal energy_ and

pressure_ The _uid velocity u _ m__ and the pressure is determined by an equation of

state_ which_ for an ideal _uid is

p _ _ _ ___e _

m_

__

__ _____d_

where _ is a constant_ Equations _____a__ _____b__ and _____c_ express the facts that

the mass_ momentum_ and energy of the _uid are neither created nor destroyed and are_

hence_ conserved_ We readily see that the system ______ has the form of ______ with

u __

_

_

m

e

__

_ f u_ __

_

m

m___ _ p

e _ p_m__

__

_ bx_ t_ u_ __

_

 

 

 

__

_ ______

Example _______ The de_ection of a taut string has the form

utt _ a_uxx _ qx__ ____          a_

_____ Conservation Laws _

x = 0 x = L

T T

u(x,t)

Figure ______ Geometry of the taut string of Example ______

where a_ _ T__ with T being the tension and _ being the linear density of the string Fig_

ure _______ The lateral loading qx_ applied in the transverse direction could represent

the weight of the string_

This second_order partial di_erential equation can be written as a _rst_order system

of two equations in a variety of ways_ Perhaps the most common approach is to let

u_ _ ut_ u_ _ aux_ ____          b_

Physically_ u_x_ t_ is the velocity and u_x_ t_ is the stress at point x and time t in the

string_ Di_erentiating with respect to t while using ____           a_ and ____    b_ yields

u__t _ utt _ a_uxx _ qx_ _ au__x _ qx__ u__t _ auxt _ autx _ au__x_

Thus_ the one_dimensional wave equation has the form of ______ with

u _ _ u_

u_ __ f u_ _ _ _cu_

_cu_ __ bx_ t_ u_ _ _ qx_

 __ ____          c_

In the convective form _______ we have

A _ _  _a

_a  __ ____     d_

______ Characteristics

The behavior of the system ______ can be determined by diagonalizing the Jacobian

_____b__ This can be done for hyperbolic systems since Au_ has m distinct eigenvalues

De_nition _______ Thus_ let

P _ _p____ p____ _ _ _ _ p_m__ _____a_

and recall the eigenvalue_eigenvector relation

AP _ P__ _____b_

_ Hyperbolic Problems

where

_ _

_____

__

__

_ _ _

_m

_____

_____c_

Multiplying _____a_ by P__ and using _____b_ gives

P__ut _ P__Aux _ P__ut _ _P__ux _ P__b_

Let

w _ P__u ______

so that

wt _ _wx _ P__ut _ P___tu _ __P__ux _ P___xu__

Using ______

wt _ _wx _ Qw _ g_ _____a_

where

Q _ _P___t _ _P___x_P_ g _ P__b_ _____b_

In component form_ _____a_ is

wi_t _ _iwi_x _

m

Xj__

qi_jwj _ gi_ i _ __ __ _ _ _ _m_ _____c_

Thus_ the transformation ______ has uncoupled the di_erentiated terms of the original

system _____a__

Consider the directional derivative of each component wi_ i _ __ __ _ _ _ _m_ of w_

dwi

dt

_ wi_t _ wi_x

dx

dt

_ i _ __ __ _ _ _ _m_

in the directions

dx

dt

_ _i_ i _ __ __ _ _ _ _m_ _____a_

and use _____c_ to obtain

dwi

dt

_

m

Xj__

qi_jwj _ gi_ i _ __ __ _ _ _ _m_ _____b_

_____ Conservation Laws      

The curves _____a_ are called the characteristics of the system ______ _______ The

partial di_erential equations ______ may be solved by integrating the _m ordinary dif_

ferential equations _____a_ _____b__ This system is uncoupled through its di_erentiated

terms but coupled through Q and g_ This method of solution is_ quite naturally_ called

the method of characteristics_ While we could develop numerical methods based on the

method of characteristics_ they are generally not e_cient when m          __

De_nition _______ The set of all points that determine the solution at a point Px__ t__

is called the domain of dependence of P_

Consider the arbitrary point Px__ t__ and the characteristics passing through it as

shown in Figure ______ The solution ux__ t__ depends on the initial data on the interval

_A_B_ and on the values of b in the region APB_ bounded by _A_B_ and the characteristic

curves x_ _ __ and x_ _ _m_ Thus_ the region APB is the domain of dependence of P_

dx/dt =

x

1 dx/dt =

P(x ,t )

A B

t

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

m 

0 0

Figure ______ Domain of dependence of a point Px__ t___ The solution at P depends on

the initial data on the line _A_B_ and the values of b within the region APB bounded

by the characteristic curves dx_dt _ ___ _m_

Example _______ Consider an initial value problem for the forced wave equation

____    a_ with the initial data

ux_ _ _ u_x__ utx_ _ _ u_ _x__ __ _ x _ __

Transforming ____       a_ using ____  b_ yields the _rst_order system ______ with A and

b given by ____           __ Using ____ b__ The initial conditions become

u_x_ _ _ u_ _x__ u_x_ _ _ au_

xx__ __ _ x _ __

_ Hyperbolic Problems

With A given by ____  __ we _nd its eigenvalues as ____ _ _a_ Thus_ the character_

istics are

x_ _ _a_

and the eigenvectors are

P _

_

p_ _ _ _

_ __ __

Since P__ _ P_ we may use ______ to determine the canonical variables as

w_ _

u_ _ u_ p_

_ w_ _

u_ _ u_ p_

_

From _______ the canonical form of the problem is

w__t _ aw__x _

q

p_

_ w__t _ aw__x _

q

p_

_

The characteristics integrate to

x _ x_ _ at_ x _ x_ _ at_

and along the characteristics_ we have

dwk

dt

_

q

p_

_ k _ __ __

Integrating_ we _nd

w_x_ t_ _ w_

_x__ _

_

p_ Z t

_

qx_ _ a_d

or

w_x_ t_ _ w_

_x__ _

_

ap_ Z x__at

x_

q__d__

It_s usual to eliminate x_ by using the characteristic equation to obtain

w_x_ t_ _ w_

_x _ at_ _

_

ap_ Z x

x_at

q__d__

Likewise

w_x_ t_ _ w_

_x _ at_ _

_

ap_ Z x

x_at

q__d__

The domain of dependence of a point Px__ t__ is shown in Figure ______ Using the

bounding characteristics_ it is the triangle connecting the points x__ t___ x_ _ at__ __

and x_ _at__ __ Actually_ with q being a function of x only_ the domain of dependence

only involves values of qx_ on the subinterval x_ _ at__ _ to x_ _ at__ ___

_____ Conservation Laws _

0

0

t

dx/dt = a dx/dt = -a

0 P(x ,t )

x

x - at x + at 0

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

____________________

0 0

Figure ______ The domain of dependence of a point Px__ t__ for Example _____ is the

triangle connecting the points P_ x_ _ at__ __ and x_ _ at__ __

Transforming back to the physical variables

u_x_ t_ _

_

p_

w_ _ w__ _

_

p_

_w_

_x _ at_ _ w_

_x _ at__ _

_

_a Z x_at

x_at

q__d__

u_x_ t_ _

_

p_

w_ _w__ _

_

p_

_w_

_x_at__w_

_x_at___

_

_a

_Z x

x_at

q__d_ _Z x

x_at

q__d___

Suppose_ for simplicity_ that u_ _x_ _ _ then

u_x_ _ _  _

_

p_

_w_

_x_ _ w_

_x___

u_x_ _ _ au_

xx_ _

_

p_

_w_

_x_ _ w_

_x___

Thus_

w_

_x_ _ _w_

_x_ _

au_

xx_

p_

_

and

u_x_ t_ _

a

_

_u_

xx _ at_ _ u_

xx _ at__ _

_

_a Z x_at

x_at

q__d__

u_x_ t_ _

a

_

_u_

xx _ at_ _ u_

xx _ at__ _

_

_a

_Z x

x_at

q__d_ _ Z x

x_at

q__d___

Since u_ _ aux_ we can integrate to _nd the solution in the original variables_ In order

to simplify the manipulations_ let_s do this with qx_ _ _ In this case_ we have

u_x_ t_ _

a

_

_u_

xx _ at_ _ u_

xx _ at___

_ Hyperbolic Problems

hence_

ux_ t_ _

_

_

_u_x _ at_ _ u_x _ at___

The solution for an initial value problem when

u_x_ _

_

x _ __ if _ _ _ x _

_ _ x_ if  _ x _ _

_ otherwise

is shown in Figure ______ The initial data splits into two waves having half the initial

amplitude and traveling in the positive and negative x directions with speeds a and _a_

respectively_

u(x,0)

-1 1 x

-1 1 x

-1 1

-1 1

x

x

u(x,1/2a)

u(x,1/a) u(x,3/2a)

Figure ______ Solution of Example _____ at t _  upper left__ ___a upper right__ __a

lower left__ and ___a lower right__

______ Rankine_Hugoniot Conditions

For simplicity_ let us neglect bx_ t_ u_ in _____a_ and consider the integral form of the

conservation law

d

dt Z _

_

udx _ _f u_j_

_ _ _f u__ t__ _ f u_ t___ ______

which states that the rate of change of u within the interval  _ x _ _ is equal to the

change in its _ux through the boundaries x _ _ __

If f and u are smooth functions_ then ______ can be written as

Z _

_

_ut _ f u_x_dx _ _

_____ Conservation Laws _

If this result is to hold for all _control volumes_ _ ___ the integrand must vanish_ and_

hence_ _____a_ and ______ are equivalent_

To further simplify matters_ let con_ne our attention to the scalar conservation law

ut _ fu_x _  ______a_

with

au_ _

dfu_

du

_ ______b_

and

ut _ au_ux _ _ ______c_

The characteristic equation is

dx

dt

_ _ _ au__ ______a_

The scalar equation ______c_ is already in the canonical form _____a__ We calculate

the directional derivative on the characteristic as

du

dt

_ utdt _ ux

dx

dt

_ ut _ au_ux _ _ ______b_

Thus_ in this homogeneous scalar case_ ux_ t_ is constant along the characteristic curve

_____a__

For an initial value problem for ______a_ on __ _ x _ __ t       _ the solution

would have to satisfy the initial condition

ux_ _ _ u_x__ __ _ x _ __ _______

Since u is constant along characteristic curves_ it must have the same value that it had

initially_ Thus_ u _ u_x__ _ u_

_

along the characteristic that passes through x__ __ From

______a__ we see that this characteristic satis_es the ordinary initial value problem

dx

dt

_ au_

_

__ t       _ x_ _ x__ _______

Integrating_ we determine that the characteristic is the straight line

x _ x_ _ au_

_

_t_ _____        _

This procedure can be repeated to trace other characteristics and thereby construct

the solution_

_ Hyperbolic Problems

0

0 x = x + at

at

u(x,0) = (x) u(x,t) = (x-at)

at

x

x

u(x,t)

t

a

x

1

Figure ____     _ Characteristic curves and solution of the initial value problem ______a_

_______ when a is a constant_

Example _______ The simplest case occurs when a is a constant and fu_ _ au_ All

of the characteristics are parallel straight lines with slope __a_ The solution of the initial

value problem ______a_ _______ is ux_ t_ _ u_x_at_ and is_ as shown in Figure ____        _

a wave that maintains its shape and travels with speed a_

Example _______ Setting au_ _ u and fu_ _ u___ in ______a_ ______b_ yields the

inviscid Burgers_ equation

ut _

_

_

u__x _ _ _______

Again_ consider an initial value problem having the initial condition ________ so the

characteristic is given by _____           _ with a_ _ ux__ _ _ u_x___ i_e__

x _ x_ _ u_x__t_ _______

The characteristics are straight lines with a slope that depends on the value of the

initial data_ thus_ the characteristic passing through the point x__ _ has slope __u_x___

_____ Conservation Laws __

The fact that the characteristics are not parallel introduces a di_culty that was not

present in the linear problem of Example ______ Consider characteristics passing through

x__ _ and x__ _ and suppose that u_x__          u_x__ for x_    x__ Since the slope of the

characteristic passing through x__ _ is less than the slope of the one passing through

x__ __ the two characteristics will intersect at a point_ say_ P as shown in Figure ______

The solution would appear to be multivalued at points such as P_

1

x

P

x x 0 1

1

0

x = x 1 + (x )t

x = x + (x )t 0 

1

0

1

t

Figure ______ Characteristic curves for two initial points x_ and x_ for Burgers_ equation

________ The characteristics intersect at a point P_

In order to clarify matters_ let_s examine the speci_c choice of u_ given by Lax ___

u_x_ _            

_

__ if x _

_ _ x_ if  _ x _ _

_ if _ _ x

_ _______

Using ________ we see that the characteristic passing through the point x__ _ satis_es

x _      

_

x_ _ t_ if x_ _

x_ _ _ _ x__t_ if  _ x _ _

x__ if _ _ x

_ _______

Several characteristics are shown in Figure ______ The characteristics _rst intersect at

t _ __ After that_ the solution would presumably be multivalued_ as shown in Figure

______

It_s_ of course_ quite possible for multivalued solutions to exist_ however_ i_ they

are not observed in physical situations and ii_ they do not satisfy ______a_ in any

classical sense_ Discontinuous solutions are often observed in nature once characteristics

of the corresponding conservation law model have intersected_ They also do not satisfy

__ Hyperbolic Problems

t

x

1

1

Figure ______ Characteristics for Burgers_ equation _______ with initial data given by

________

u(x,3/2)

u(x,1/2)

u(x,1)

u(x,0)

0 1 x 0 1 2 x

0 1 2 x

2

0 1 2 x

Figure ______ Multivalued solution of Burgers_ equation _______ with initial data given

by ________ The solution ux_ t_ is shown as a function of x for t _ _ ____ __ and ____

______a__ but they might satisfy the integral form of the conservation law _______ We

examine the simplest case when two classical solutions satisfying ______a_ are separated

by a single smooth curve x _ _t_ across which ux_ t_ is discontinuous_ For each t        

we assume that  _ _t_ _ _ and let superscripts _ and _ denote conditions immediately

_____ Conservation Laws __

to the left and right_ respectively_ of x _ _t__ Then_ using _______ we have

d

dt Z _

_

udx _

d

dt

_Z __

_

udx _ Z _

__

udx_ _ _fu_j_

_

or_ di_erentiating the integrals

Z __

_

utdx _ u____ _ Z _

__

utdx _ u____ _ _fu_j_

__

The solution on either side of the discontinuity was assumed to be smooth_ so ______a_

holds in _ ___ and ___ __ and can be used to replace the integrals_ Additionally_ since

_ is smooth_ ___ _ ___ _ ___ Thus_ we have

_fu_j__

_ _ u___ _ fu_j_

__ _ u___ _ _fu_j_

__

or

__u_ _ u__ _ fu__ _ fu___ ______

Let

_q_ _ q_ _ q_ ______a_

denote the jump in a quantity q and write ______ as

_u___ _ _fu___ ______b_

Equation ______b_ is called the Rankine_Hugoniot jump condition and the discontinuity

is called a shock wave_ We can use the Rankine_Hugoniot condition to _nd a discontinuous

solution of Example ____         _

Example _____           _ For t _ __ the discontinuous solution of _______ _______ is as given

in Example ____          _ For t _ __ we hypothesize the existence of a single shock wave_ passing

through __ __ in the x_ t__plane_ As shown in Figure ______ the solution of Example

____    can be used to infer that u_ _ _ and u_ _ _ Thus_ fu__ _ u_____ _ ___ and

fu__ _ u_____ _ _ Using ______b__ the velocity of the shock wave is

__ _

_

_

_

Integrating_ we _nd the shock location as

_ _

_

_

t _ c_

__ Hyperbolic Problems

= (t + 1)/2

0 1

t

x

1

Figure ______ Characteristics and shock discontinuity for Example ______

u(x,0)

u(x,1) u(x,3/2)

1 2 x 0 2 x

0 1 x 0 1 2 x

0 1

2

u(x,1/2)

Figure ______ Solution ux_ t_ of Example _____ as a function of x at t _ _ ____ __ and

____ The solution is discontinuous for t             __

Since the shock passes through __ ___ the constant of integration c _ ____ and

_ _

_

_

t _ ___ _______

_____ Conservation Laws _   

The characteristics and shock wave are shown in Figure _____ and the solution ux_ t_

is shown as a function of x for several times in Figure ______

Let us consider another problem for Burgers_ equation with di_erent initial conditions

that will illustrate another structure that arises in the solution of nonlinear hyperbolic

systems_

Example ______ Consider Burgers_ equation _______ subject to the initial conditions

u_x_ _

_

_ if x _

x_ if  _ x _ _

__ if _ _ x

_ _______

Using _______ and ________ we see that the characteristic passing through x__ _ sat_

is_es

x _      

_

x__ if x _

x__ _ t__ if  _ x _ _

x_ _ t_ if _ _ x

_ _______

These characteristics_ shown in Figure _______ may be used to verify that the solution_

shown in Figure _______ is continuous_ Additional considerations and di_culties with

nonlinear hyperbolic systems are discussed in Lax ____

Example _______ A Riemann problem is an initial value Cauchy_ problem for ______

with piecewise_constant initial data_ Riemann problems play an important role in the

numerical solution of conservation laws using both _nite di_erence and _nite element

techniques_ In this introductory section_ let us illustrate a Riemann problem for the

inviscid Burgers_ equation ________ Thus_ we apply the initial data

ux_ _ _ _ uL_ if x _

uR_ if x _

_ _____           _

As in the previous two examples_ we have to distinguish between two cases when

uL         uR and uL _ uR_ The solution may be obtained by considering piecewise_linear

continuous initial conditions as in Examples _____ and ______ but with the _ramp_

extending from  to _ instead of from  to __ We could then take a limit as _ _ _ The

details are left to an exercise Problem _ at the end of this section__

When uL           uR_ the characteristics emanating from points x_ _  are the straight

lines x _ x_ _ uLt cf_ _________ Those emanating from points x_        are x _

x_ _ uRt_ The characteristics cross immediately and a shock forms_ Using _______ we

see that the shock moves with speed __ _ uL _uR____ The solution is constant along the

characteristics and_ hence_ is given by

ux_ t_ _ _ uL_ if x_t _ uL _ uR___

uR_ if x_t _ uL _ uR___

_ uL      uR_ ______a_

__ Hyperbolic Problems

t

0 1 x

1

Figure _______ Characteristics for Example ______

u(x,0) u(x,1/2)

1 2 x 0 2 x

0 1 2 x

0

1

1

1

2

1

1

1

u(x,1)

0 x

u(x,3/2)

Figure _______ Solution ux_ t_ of Example _____ as a function of x at t _ _ ____ __ and

____

_____ Conservation Laws __

Several characteristics and the location of the shock are shown in Figure _______

When uL _ uR_ the characteristics do not intersect_ There is a region between the

characteristic x _ uLt emanating from x_ _ _ and x _ uRt emanating from x_ _ _

where the initial conditions fail to determine the solution_ As determined by either

the limiting process suggested in Problem _ or thermodynamic arguments using entropy

considerations ____ no shock forms and the solution in this region is an expansion fan_

Several characteristics are shown in Figure ______ and the expansion solution is given

by

ux_ t_ _          

_

uL_ if x_t _ uL

x_t_ if uL _ x_t _ uR

uR_ if x_t _ uR

_ uL _ uR_ ______b_

t

1/u

1/u

R

L

x

t

1/u

1/uR

L

x

Figure _______ Shock left_ and expansion right_ wave characteristics of the Riemann

problem of Example ______

We conclude this example by examining the solution of the Riemann problem along

the line x _ _ Characteristics for several choices of initial data are shown in Figure

______ and_ by examining these and ________ we see that

u_ t_ _

            _

uL_ if uL_ uR    

uR_ if uL_ uR _

_ if uL _ _ uR    

uL_ if uL           _ uR _ _ uL _ uR___  

uR_ if uL           _ uR _ _ uL _ uR___ _

_

This data will be useful when constructing numerical schemes based on the solution of

Riemann problems_

Problems

__ Hyperbolic Problems

t

x

t

t

x

x

x x

t

t

t

x

Figure _______ Characteristics of Riemann problems for Burgers_ equation when uL_ uR    

 top__ uL_ uR _  center__ uL  _ uR _ _ uL _ uR___   bottom left__ and

uL _ _ uR          bottom right__

__ Show that the solution of the Riemann problem _______ _____     _ is given by

________ You may begin by solving a problem with continuous initial data_ e_g__

ux_ _ _           

_

uL_ if x _ __

uL

__ _ _ x_ _ uR

__ _ _ x__ if _ _ _ x _ _

uR_ if _ _ x

_

and take the limit as _ _ _

_____ Discontinuous Galerkin Methods __

____ Discontinuous Galerkin Methods

In Section ____ we examined the use of the discontinuous Galerkin method for time

integration_ We_ll now examine it as a way of performing spatial discretization of con_

servation laws _______ The method might have some advantages when solving problems

with discontinuous solutions_ The discontinuous Galerkin method was _rst used for to

solve an ordinary di_erential equation for neutron transport _____ At the moment_ it

is very popular and is being used to solve ordinary di_erential equations ____ ___ and

hyperbolic _     _ __ __ __ ___ ___ ___ ____ parabolic ____ _        __ and elliptic ___ __ ___ partial

di_erential equations_ A recent proceedings contains a complete and current survey of

the method and its applications ____

The discontinuous Galerkin method has a number of advantages relative to traditional

_nite element methods when used to discretize hyperbolic problems_ We have already

noted that it has the potential of sharply representing discontinuities_ The piecewise

continuous trial and test spaces make it unnecessary to impose interelement continuity_

There is also a simple communication pattern between elements that makes it useful for

parallel computation_

We_ll begin by describing the method for conservation laws ______ in one spatial

dimension_ In doing this_ we present a simple construction due to Cockburn and Shu ____

rather than the more standard_ approach ____ used in Section ___ for time integra_

tion_ Using a method of lines formulation_ let us divide the spatial region into elements

xj___ xj__ j _ __ __ _ _ _ _N_ and construct a local Galerkin problem on Element xj___ xj_

in the usual manner by multiplying _____a_ by a test function v and integrating to

obtain

Z xj

xj__

vT _ut _ f u_x_dx _ _ _____a_

The loading term bx_ t_ u_ in _____a_ causes no conceptual or practical di_culties and

we have neglected it to simplify the presentation_

Following the usual procedure_ let us map xj___ xj_ to the canonical element ___ __

using the linear transformation

x _

_ _ _

_

xj__ _

_ _ _

_

xj _ _____b_

Then_ after integrating the _ux term in _____a_ by parts_ we obtain

hj

_ Z _

__

vTutd_ _ vT f u_j_

__ _ Z _

__

vT

_ f u_d_ _____c_

where

hj _ xj _ xj___ _____d_

_ Hyperbolic Problems

Without a need to maintain interelement continuity_ there are several options available

for selecting a _nite element basis_ Let us choose one based on Legendre polynomials_

As we shall see_ this will produce a diagonal mass matrix without a need to use lumping_

Thus_ we select the approximation Ujx_ t_ of ux_ t_ on the mapping of xj___ xj_ to the

canonical element as

Uj__ t_ _

p

Xk__

ckjt_Pk__ _____a_

where ckjt_ is an m_vector and Pk__ is the Legendre polynomial of degree k in __ Recall

cf_ Section __  __ that the Legendre polynomials satisfy the orthogonality relation

Z _

__

Pi__Pj__d_ _

__ij

_i _ _

_ i_ j _  _____b_

are normalized as

Pi__ _ __ i _ _ _____c_

and satisfy the symmetry relation

Pi__ _ ___iPi____ i _ _ _____d_

The _rst six Legendre polynomials are

P___ _ __ P___ _ __

P___ _

___ _ _

_

_ P___ _

            __ _ __

_

_

P          __ _

_          _          _ ___ _ _

_

_ P__ _

___ _ ___ _ _  _

_

_ ______

These polynomials are illustrated in Figure _______ Additional information appears in

Section __       and Abromowitz and Stegen ____

Substituting _____a_ into _____c__ testing against Pi___ and using _____b_d_ yields

hj _cij

_i _ _

_ f Uxj_ t__ _ ___if Uxj___ t__ _ Z _

__

dPi__

d_

f Uj__ t__d__

i _ __ __ _ _ _ _ p_ _____a_

where __ _ d __dt_

Neighboring elements must communicate information to each other and_ in this form

of the discontinuous Galerkin element method_ this is done through the boundary _ux

_____ Discontinuous Galerkin Methods __

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Figure ______ Legendre polynomials of degrees p _ _ __ _ _ _ _       _

terms_ The usual practice is to replace the boundary _ux terms f Uxk_ t__ k _ j _ __ j_

by a numerical _ux function

f Uxk_ t_          FUkxk_ t___Uk__xk_ t__ _____b_

that depends on the approximate solutions Uk and Uk__ on the two elements sharing the

vertex at xk_ Cockburn and Shu ____ present several possible numerical _ux functions_

Perhaps_ the simplest is the average

FUkxk_ t___Uk__xk_ t__ _

f Ukxk_ t__ _ f Uk__xk_ t__

_

_ ____ a_

Based on our work with convection_di_usion problems in Section __  _ we might expect

that some upwind considerations might be worthwhile_ This happens to be somewhat

involved for nonlinear vector systems_ We_ll postpone it and_ instead_ note that an

upwind _ux for a scalar problem is

FUkxk_ t___ Uk__xk_ t__ _ _ fUkxk_ t___ if aUkxk_ t__ _ aUk__xk_ t__             

fUk__xk_ t___ if aUkxk_ t__ _ aUk__xk_ t__ _

____    b_

where

au_ _ fuu__ ____        c_

__ Hyperbolic Problems

A simple numerical _ux that is relatively easy to apply to vector systems and employs

upwind information is the Lax_Friedrichs function ____

FUkxk_ t__Uk__xk_ t__ _

_

_

_f Ukxk_ t__ _ f Uk__xk_ t__

__maxUk__xk_ t_ _ Ukxk_ t____ ____        d_

where _max is the maximum absolute eigenvalue of the Jacobian matrix fuu__ u

_Ukxk_ t___Uk__xk_ t___

Example _______ The simplest discontinuous Galerkin scheme uses piecewise_constant

p _ _ solutions

Uj__ t_ _ c_jt_P___ _ c_j _

In this case_ _____a_ becomes

hj _c_j _ f Uxj_ t__ _ f Uxj___ t__ _ __

In this initial example_ let_s choose a scalar problem and evaluate the _ux using the

average ____   a_

FUkxk_ t___ Uk__xk_ t__ _

fUkxk_ t__ _ fUk__xk_ t__

_

_

fc__k_ _ fc__k___

_

and upwind ____         b_

FUkxk_ t___ Uk__xk_ t__ _ _ fc__k__ if ac__k_ _ ac__k___          

fc__k____ if ac__k_ _ ac__k___ _

numerical _uxes_ With these _ux choices_ we have the ordinary di_erential systems

_ c_j _

fc__j___ _ fc__j___

_hj

_

and

_ c_j _

_ _ _j_fc__j___ _ ___j_fc__j_ _ _ _ _j___fc__j_ _ _ _ _j___fc__j___

_hj

_

where

_j _ sgnac__j_ _ ac__j_____

In the simplest_ case when fu_ _ au with a a positive constant_ we have the two

schemes

_ c_j _

ac__j__ _ c__j___

_hj

_  j _ _ __ _ _ _ _ J_

and

_ c_j _

ac__j _ c__j___

hj

_ _ j _ _ __ _ _ _ _ J_

_____ Discontinuous Galerkin Methods __

Initial conditions for c_j_ may be speci_ed by interpolating the initial data at the center

of each interval_ i_e__ c__j_ _ u_xj _ hj____ j _ __ __ _ _ _ _ J_

We use these two techniques to solve an initial value problem with a _ _ and

u_x_ t_ _ sin__x_

Thus_ the exact solution is

ux_ t_ _ sin__x _ t__

Piecewise_constant discontinuous Galerkin solutions with upwind and centered _uxes

are shown at t _ _ in Figure ______ A ___element uniform mesh was used and time inte_

gration was performed using the MATLAB Runge_Kutta procedure ode_      _ The solution

with the upwind _ux has greatly dissipated the solution after one period in time_ The

maximum error at cell centers

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

x

U

Figure ______ Exact and piecewise_constant discontinuous solutions of a linear kinematic

wave equation with sinusoidal initial data at t _ __ Solutions with upwind and centered

_uxes are shown_ The solution using the upwind _ux exhibits the most dissipation_

je__ t_j_ __ max

__j_J juxj _ hj___ t_ _ Uxj _ hj___ t_j

at t _ _ is shown in Table _____ on meshes with J _ ___ ___ and __ elements_ Since

the errors are decreasing by a factor of two for each mesh doubling_ it appears that the

__ Hyperbolic Problems

upwind__ux solution is converging at a linear rate_ Using similar reasoning_ the centered

solution appears to converge at a quadratic rate_ The errors appear to be smallest at the

downwind right_ end of each element_ This superconvergence result has been known for

some time ____ but other more general results were recently discovered ____

J Upwind Centered

jej_ jej_

__ ____ __      __

__ __   __ __

__ ___ _ ____

Table ______ Maximum errors for solutions of a linear kinematic wave equation with

sinusoidal initial data at t _ _ using meshes with J _ ___ ___ and __ uniform elements_

Solutions were obtained using upwind and centered _uxes_

As a second calculation_ let_s consider discontinuous initial data

u_x_ t_ _ _ __ if  _ x _ ___

___ if ___ _ x _ _

_

This data is extended periodically to the whole real line_ Piecewise_constant discontin_

uous Galerkin solutions with upwind and centered _uxes are shown at t _ _ in Fig_

ure ______ The upwind solution has_ once again_ dissipated the initial square pulse_

This time_ however_ the centered solution is exhibiting spurious oscillations_ As with

convection_dominated convection_di_usion equations_ some upwinding will be necessary

to eliminate spurious oscillations near discontinuities_

______ High_Order Discontinuous Galerkin Methods

The results of Example _____ are extremely discouraging_ It would appear that we have

to contend with either excessive di_usion or spurious oscillations_ To overcome these

choices_ we investigate the use of the higher_order techniques o_ered by _______ With

cij being an m_vector and i ranging from  to p_ we have p__ vector and mp___ scalar

unknowns on each element_

We will focus on the four major tasks_ i_ evaluating the integral on the right side

of _____a__ ii_ performing the time integration iii_ de_ning the initial conditions_

and iv_ evaluating the _uxes_ The integral in _____a_ will typically require numerical

integration and the obvious choice is Gaussian quadrature as described in Chapter __

This works _ne and there is no need to discuss it further_

Time integration can be performed by either explicit or implicit techniques_ The

choice usually depends on the spread of the eigenvalues _i_ i _ __ __ _ _ _ _m_ of the Jaco_

bian Au__ If the eigenvalues are close to each other_ explicit integration is _ne_ Stability

_____ Discontinuous Galerkin Methods _      

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

x

U

Figure ______ Exact and piecewise_constant discontinuous solutions of a linear kinematic

wave equation with discontinuous initial data at t _ __ Solutions with upwind and

centered _uxes are shown_ The solution using the upwind _ux is dissipative_ The solution

using the centered _ux exhibits spurious oscillations_

is usually not a problem_ An implicit scheme might be necessary when the eigenvalues are

widely separated or when integrating ______ to a steady state_ For explicit integration_

Cockburn and Shu ____ recommend a total variation diminishing TVD_ Runge_Kutta

scheme_ However_ Biswas et al_ ___ found that classical Runge_Kutta formulas gave sim_

ilar results_ Second_ and third_order and fourth_ and _fth_order classical Runge_Kutta

software was used for time integration of Example ______ If forward Euler integration of

_____a_ were used_ we would have to solve the explicit system

hj

_i _ _

cn__

ij _ cn

ij

_t

_ _f Unxj__ _ ___if Unxj____ _ Z _

__

dPi__

d_

f Unj

___d__

i _ __ __ _ _ _ _ p_

The notation is identical to that used in Chapter __ thus_ Unx_ and cn

ij are the approx_

imations of Ux_ tn_ and cijtn__ respectively_ produced by the time integration software

and _t is the time step_ The forward Euler method is used for illustration because of its

simplicity_ The order of the temporal integration method should be comparable to p_

__ Hyperbolic Problems

Initial conditions may be determined by L_ projection as

Z _

__

Pi___Uj__ _ _ u____d_ _ __ i _ _ __ _ _ _ _ p_ j _ __ __ _ _ _ _ J_ ______

One more di_culty emerges_ Higher_order schemes for hyperbolic problems oscillate

near discontinuities_ This is a fundamental result that may be established by theoretical

means cf__ e_g__ Sod __       ___ One technique for reduced these oscillations involves limiting

the computed solution_ Many limiting algorithms have been suggested but none are

totally successful_ We describe a procedure for limiting the slope _Ujx_ t___x of the

solution that is widely used_ With this approach_ _Ujx_ t___x is modi_ed so that_

__ the solution _____a_ does not take on values outside of the adjacent grid averages

Figure ______ upper left__

__ local extrema are set to zero Figure ______ upper right__ and

__ the gradient is replaced by zero if its sign is not consistent with its neighbors Figure

______ lower center__

Figure _____ illustrates these situations when the solution is a piecewise_linear p _ __

function relative to the mesh_

A formula for accomplishing this limiting can be summarized concisely using the

minimum modulus function as

_Uj_modxj_ t_

_x

_ minmod

_Ujxj_ t_

_x

_rUjxj_____ t___Ujxj_____ t__ _____a_

_Uj_modxj___ t_

_x

_ minmod

_Ujxj___ t_

_x

_rUjxj_____ t___Ujxj_____ t__ _____b_

where

minmoda_ b_ c_ _ _ sgna_minjaj_ jbj_ jcj__ if sgna_ _ sgnb_ _ sgnc_

_ otherwise

_____c_

and r and _ are the backward and forward di_erence operators

rUjxj_____ t_ _ Ujxj_____ t_ _ Ujxj_____ t__ _____d_

and

_Ujxj_____ t_ _ Ujxj_____ t_ _ Ujxj_____ t__ _____e_

With _Uj_modxj___ t___x and _Uj_modxj_ t___x_ determined_ _____a_b_ are used to re_

computed the coe_cients in _____a_ to reduce the oscillations_ However_ _____a_b_

_____ Discontinuous Galerkin Methods __

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

______

______

______

______

______

______

______

______

______

______

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

______

______

______

______

______

______

______

______

______

______

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

______

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

_____

j _____

j

j

Figure ______ Solution limiting_ reduce slopes to be within neighboring averages upper

left__ set local extrema to zero upper right__ and set slopes to zero if they disagree with

neighboring trends_

only provide two vector equations for modifying the p vector coe_cients cij_modt__ i _

__ __ _ _ _ _ p_ in _Ujx_ t___x_ When p _ __ _____a_b_ are identical and c_j_modt_ is

uniquely determined_ Likewise_ when p _ __ the two conditions _____a_b_ su_ce to

uniquely determine the modi_ed coe_cients c_j_modt_ and c_j_modt__ Equations _____a_b_

are insu_cient to determine the modi_ed coe_cients when p      _ and Cockburn and

Shu ____ suggested setting the higher_order coe_cients cij_modt__ i _ __ __ _ _ _ _ p_ to zero_

This has the disturbing characteristic of __attening_ the solution near smooth extrema

and reducing the order of accuracy_ Biswas et al_ ___ developed an adaptive limiter which

__ Hyperbolic Problems

applied the minimum modulus function _____c_ to higher derivatives of Uj_ They began

by limiting the p th derivative of Uj and worked downwards until either a derivative was

not changed by the limiting or they modi_ed all of the coe_cients_ Their procedure_

called _moment limiting__ is described further in their paper ____

Example _______ Biswas et al_ ___ solve the inviscid Burgers_ equation _______ with

the initial data

ux_ _ _

_ _ sin x

_

_

This initial data steepens to form a shock which propagates in the positive x direction_

Biswas et al_ ___ use an upwind numerical _ux ____  b_ and solve problems on uniform

meshes with h _ ____ with p _ _ __ __ Time integration was done using classical Runge_

Kutta methods of orders ____ respectively_ for p _ _ __ __ Exact and computed solutions

are shown in Figure ____         _ The piecewise polynomial functions used to represent the

solution are plotted at eleven points on each subinterval_

The _rst_order solution p _ _ shown at the upper left of Figure ____   is character_

istically di_usive_ The second_order solution p _ __ shown at the upper right of Figure

____    has greatly reduced the di_usion while not introducing any spurious oscillations_

The minimum modulus limiter ______ has _attened the solution near the shock as seen

with the third_order solution p _ __ shown at the lower left of Figure ____      _ There is

a loss of local_ monotonicity near the shocks_ Average solution values are monotone

and this is all that the limiter ______ was designed to produce__ The adaptive moment

limiter of Biswas et al_ ___ reduces the _attening and does a better job of preserving local

monotonicity near discontinuities_ The solution with p _ _ using this limiter is shown in

the lower portion of Figure ____          _

Example _______ Adjerid et al_ ___ solve the nonlinear wave equation

utt _ uxx _ u_u_ _ __ _____a_

which can be written in the form _____a_ as

u__t _ u__x _ u__ u__t _ u__x _ u__u_

_ _ __ _____b_

with u_ _ u_ The initial and boundary conditions are such that the exact solution of

_____a_ is the solitary wave

ux_ t_ _ sechx cosh

_

_

_ t sinh

_

_

_ _____c_

cf_ Figure _______

Adjerid et al_ ___ solved problems on ____ _ x _ ____  _ t _ _ by the discontin_

uous Galerkin method using polynomials of degrees p _  to __ The solution at t _ _

_____ Discontinuous Galerkin Methods __

Figure ____     _ Exact line_ and discontinuous Galerkin solutions of Example _____ for

p _ _ __ __ and h _ _____ Solutions with the minmod limiter ______ and an adaptive

moment limiter of Biswas et al_ ___ are shown for p _ __

performed with p _ _ and J _ __ is shown in Figure ______ The entire solitary wave is

shown_ however_ the computation was performed on the center region ____ _ x _ ____

_ Hyperbolic Problems

Discretization errors in the L_ norm

ke__ t_k _

J

Xj__

Z xj

xj__ jUx_ t_ _ Ujx_ t_jdx

are presented for the solution u for various combinations of h and p in Table ______

Solutions of this nonlinear wave propagation problem appear to be converging as Ohp___

in the L_ norm_ This can be proven correct for smooth solutions of discontinuous Galerkin

methods ___ ___ ____

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-10 -8 -6 -4 -2 0 2 4 6 8 10

Figure ______ Solution of Example _____ at t _ _ obtained by the discontinuous Galerkin

method with p _ _ and N _ ___

J p _  p _ _ p _ _ p _ _ p _ _

_ ____e__       ___e__ ____e__ ____e__ ____e__

__ ____e__ ____e__ ____e_ ____e__ ____e__

__ ____e__ ____e__ ___e__ ___e__ __                  e__

__ ____e__ ___e_      ____e__ ____e__ ___e___

___ ____e__ ____e_  __       _e__ ___e__ ____e___

_          _ __     _e__ ____e__             ___e__

Table ______ Discretization errors at t _ _ as functions J and p for Example ______

Evaluating numerical _uxes and using limiting for vector systems is more complicated

than indicated by the previous scalar example_ Cockburn and Shu ____ reported problems

when applying limiting component_wise_ At the price of additional computation_ they

applied limiting to the characteristic _elds obtained by diagonalizing the Jacobian fu_

Biswas et al_ ___ proceeded in a similar manner_ _Flux_vector splitting_ may provide a

compromise between the two extremes_ As an example_ consider the solution and _ux

vectors for the one_dimensional Euler equations of compressible _ow _______ For this

_____ Discontinuous Galerkin Methods __

and related di_erential systems_ the _ux vector is a homogeneous function that may be

expressed as

f u_ _ Au _ fuu_u_ _____a_

Since the system is hyperbolic_ the Jacobian A may be diagonalized as described in

Section ___ to yield

f u_ _ P___Pu _____b_

where the diagonal matrix _ contains the eigenvalues of A

_ _

_____

__

__

_ _ _

_m

_____

__

_

u _ c

u

u _ c

__

_ _____c_

The variable c _ p_p___ is the speed of sound in the _uid_ The matrix _ can be

decomposed into components

_ _ __ _ __ _____a_

where __ and __ are_ respectively_ composed of the non_negative and non_positive com_

ponents of _

__

i _

_i _ j_ij

_

_ i _ __ __ _ _ _ _m_ _____b_

Writing the _ux vector in similar fashion using ______

f u_ _ P____ _ ___Pu _ f u__ _ f u___ _____c_

Split _uxes for the Euler equations were presented by Steger and Warming _____ Van

Leer ____ found an improvement that provided better performance near sonic and stag_

nation points of the _ow_ The split _uxes are evaluated by upwind techniques_ Thus_ at

an interface x _ xj _ f_ is evaluated using Ujxj_ t_ and f_ is evaluated using Uj__xj_ t__

Calculating _uxes based on the solution of Riemann problems is another popular

way of specifying numerical _uxes for vector systems_ To this end_ let wx_t_ uL_ uR_

be the solution of a Riemann problem for _____a_ with the peicewise_constant initial

data _____      __ The solution of a Riemann problem _breaking_ at xj_ tn_ would be

wx_xj__t_tn__Ujxj_ tn__Ujxj___ tn___ Using this_ we would calculate the numerical

_ux at xj_ t__ t             tn_ as

FUjxj_ tn__Uj__xj_ tn__ _ f w_Ujxj_ tn__Uj__xj_ tn___ _______

__ Hyperbolic Problems

Example _______ Let us calculate the numerical _ux based on the solution of a Rie_

mann problem for Burgers_ equation ________ Using the results of Example ______ we

know that the solution of the appropriate Riemann problem is

w_ Uj_ Uj___ _

            _

Uj _ if Uj_ Uj__           

Uj___ if Uj_ Uj__ _

_ if Uj _ _ Uj__            

Uj _ if Uj          _ Uj__ _ _ Uj _ Uj_____        

Uj___ if Uj       _ Uj__ _ _ Uj _ Uj_____ _

_

The arguments of Uj and Uj__ are all xj_ tn__ These have been omitted for clarity__

With fu_ _ u___ for Burgers_ equation_ we _nd the numerical _ux

FUj_ Uj___ _

            _

U_

j ___ if Uj_ Uj__          

U_

j_____ if Uj_ Uj__ _

_ if Uj _ _ Uj__            

U_

j ___ if Uj         _ Uj__ _ _ Uj _ Uj_____        

U_

j_____ if Uj      _ Uj__ _ _ Uj _ Uj_____ _

_

Letting

u_ _ maxu_ __ u_ _ minu_ __

we can write the numerical _ux more concisely as

FUj_ Uj___ _ max_U_

j _____ U_

j________

When used with a piecewise_constant basis and forward Euler time integration_ the result_

ing discontinuous Galerkin scheme is identical to Godunov_s _nite di_erence scheme _____

This was the _rst di_erence scheme to be based on the solution of a Riemann problem_

This early work and a subsequent work of Glimm ____ and Chorin ___ stimulated a great

deal of interest in using Riemann problems to construct numerical _ux functions_ A

summary of a large number of choices appears in Cockburn and Shu _____

____ Multidimensional Discontinuous Galerkin Meth_

ods

Let us extend the discontinuous Galerkin method to multidimensional conservation laws

of the form

ut _ r _ fu_x _ bx_ y_ z_ t_ u__ x_ y_ z_   _ t              _ _____a_

where

fu_ _ _f u__ gu__ hu__ _____b_

_____ Multidimensional Discontinuous Galerkin Methods __

and

r _ fu_ _ f u_x _ gu_y _ hu_z_ _____c_

The solution ux_ y_ z_ t__ componenets of the _ux vector f u__ gu__ and hu__ and the

loading bx_ y_ z_ t_ u_ are m_vectors and  is a bounde region of ___ Boundary conditions

must be prescribed on _ along characteristics that enter the region_ We_ll see what this

means by example_ Initial condtions prescribe

ux_ y_ z_ _ _ __ x_ y_ z_     _ _ _____d_

Following our analysis of Section ____ we partition  into a set of _nite elements  j _

j _ __ __ _ _ _ _N__ and construct a weak form of the problem on an element_ This is done_

as usual_ by multiplying _____a_ by a test function v  L_ j__ integrating over  j _ and

applying the divergence theorem to the _ux to obtain

v_ ut_j_ _ v_ f _ n       j _rv_ f_j _ v_ b_j_ _v  L_ j__ _____a_

where

v_ u_j _ Z_j

vTudxdydz_ _____b_

rv_ f_j _ Z_j

_vT

x f u_ _ vT

y gu_ _ vT

z hu__dxdydz_ _____c_

f _ n _ fn _ f u_n_ _ gu_n_ _ hu_n__ _____d_

and

_ v_ f _ n         j_ Z__j

vTfndS_ _____e_

The vector n _ _n__ n__ n__T is the unit outward normal vector to _ and dS is a surface

in_nitessimal on _ j _

Only the normal component of the _ux is involved in _______ hence_ its approxi_

mation on _ j is the same as the one_dimensional problems of Section ____ Thus_ the

numerical normal _ux function can be taken as a one_dimensional numerical _ux using

solution values on each side of _ j _ In order to specify this more precisely_ let nbj_k_

k _ __ __ _ _ _ _NE_ denote the indices of the NE elements sharing the bounding faces of

 j and let _ j_k_ k _ __ __ _ _ _ _NE_ be the faces of  j Figure _______ Then_ we write

_____a_ in the more explicit form

v_ ut_j _

NE

Xk__

_ v_FnUj_Unbj_k _    j_k _rv_ f_j _ v_ b_j_ _v  L_ j__ ______

__ Hyperbolic Problems



j

nb

nb

j,1

nb

j,3 j,2

Figure ______ Element j and its neighboring elements indicating that the segments _ j_k_

k _ __ __ _ _ _ _NE_ _

Without the need to maintain inter_element continuity_ virtually any polynomial basis

can be used for the approximate solutionUjx_ y_ z_ t_ on  j _ Tensor products of Legendre

polynomials can provide a basis on square or cubic canonical elements_ but these are

unavailable for triangles and tetrahedra_ Approximations on triangles and tetrahedra can

use a basis of monomial terms_ Focusing on two_dimensional problems on the canonical

right _  __ triangle_ we write the _nite element solution in the usual form

Ujx_ y_ t_ _

np

Xk__

ckjNk__ ___ ______

where np _ p _ __p _ ____ is the number of monomial terms in a complete polynomial

of degree p_ A basis of monomial terms would set

N_ _ __ N_ _ __ N_ _ __ _ _ _ _ Nnp _ _p_ ____   _

_____ Multidimensional Discontinuous Galerkin Methods _    

All terms in the mass matrix can be evaluated by exact integration on the canonical

triangle cf_ Problem _ at the end of this section_ as long as it has straight sides_

however_ without orthogonality_ the mass matrix will not be diagonal_ This is not a

severe restriction since the mass matrix is independent of time and_ thus_ need only be

inverted factored_ once_ The ill_conditioning of the mass matrix at high p is a more

important concern with the monomial basis ____         __

Ill_conditioning can be reduced and the mass matrix diagonalized by extracting an

orthogonal basis from the monomial basis ____           __ This can be done by the Gram_

Schmidt orthogonalization process shown in Figure ______ The inner product and norm

are de_ned in L_ on the canonical element as

procedure gramN_

!N

_ __ N__kN_k___

for k __ _ to np do

t __ Nk _Pk__

i__ Nk_ !N

i__ !N

i

!N

k __ t_ktk___

bf end for

return !N

Figure ______ Gram_Schmidt process to construct an orthogonal basis !N

k k _ __ __ _ _ _ _ np

from a basis of monomials Nk_ k _ __ __ _ _ _ _ np _

u_ v__ _ Z _

_ Z ___

_

uvd_d__ kuk___ _ u_ u____

_ _ _____a_

The result of the Gram_Schmidt process is a basis !N

k_ k _ __ __ _ _ _ _ np that satis_es the

orthogonality condition

 !N

i_ !N

k_ _ _i_k_ i_ k _ __ __ _ _ _ _ np_ _____b_

The actual process can be done using symbolic computation using a computer algebra

system such as MAPLE or MATHEMATICA cf_ Remacle et al_ ____ and Problem _

at the end of this section__

Example _______ We will illustrate some results using the discontinuous Galerkin

method to solve two_ and three_dimensional compressible _ow problems involving the

__ Hyperbolic Problems

Euler equations_ This complex nonlinear system has the form of _____a_ with

u _

______

_

m

n

l

e

______

_ fu_ _ _f u__ gu__ hu__ _

______

m n l

m___ _ p nm__ lm__

mn__ n___ _ p ln__

ml__ nl__ l___ _ p

e _ p_m__ e _ p_n__ e _ p_b__

______

_

bx_ t_ u_ _

______

 

 

 

 

 

______

_ _____a_

Here_ _ is the _uid density_ m_ n_ and l are the Cartesian components of the momentum

vector per unit volume_ e is the total energy per unit volume_ and p is the pressure_ which

must satisfy an equation of state of the form

p _ _ _ ___e _ m_ _ n_ _ l_______ _____b_

This equation of state assumes an ideal _uid with gas constant __

Let us consider a classical Rayleigh_Taylor instability which has a heavy _ _ __ _uid

above a light _ _ __ _uid Figure _______ This hydrostatic con_guration is unstable and

any slight perturbation will cause the heavier _uid to fall and the lighter one to rise_ The

_uid motion is quite complex and Remacle et al_ ____ simulated it using discontinuous

Galerkin methods_ They considered two_dimensional motion l _ _ ___z _  in _______

with the initial perturbation

_ _ _ __ if  _ y _ ___

__ if ___ _ y _ _

_ p _ _ ___ _ y_ if  _ y _ ___

__ _ y__ if ___ _ y _ _

u _ _x sin __x cos _y sin___ _y_ v _ __y cos _ sin_ _y_

Here u_ v_ and w are the Cartesian velocity components and _ _        ___  _ __ and _x and

_y were chosen to be small_ The boundary conditions specify that u _  on the sides and

top and v _  on the bottom_

Solutions for the density _ at t _ ___ are shown in Figure ___ for computations

with p _  to __ The mesh used for all values of p is shown in Figure ____ The total

number of vector degrees of freedom for two_dimensional discontinuous Galerkin methods

is N_np_ Since there are four unknowns per element __ m_ n_ and e_ for two_dimensional

_ows_ there are ____ ____ _____ and ___ unknowns for degrees p _ _ __ __ and

__ respectively_ Fluxes were evaluated using Roe_s linearized _ux approximation _____

No limiting was used for this computation_ A high_frequency _ltering ____ was used to

suppress oscillations in the vicinity of the interface separating the two _uids_

_____ Multidimensional Discontinuous Galerkin Methods __

1/4

1/2

1/2



Figure ______ Con_guration for the Rayleigh_Taylor instability of Example ______ There

are solid walls on the bottom and sides and open _ow at the top_

The results with p _  show very little structure of the solution_ Those with p _ _

show more_and_more detail of the _ow_ There is no exact solution of this problem_ so

it is not possible to appraise the e_ects of using higher degree polynomials_ however_

solutions with more detail are assumed to be more correct_

Remacle et al_ ____ also did computations using adaptive p_re_nement_ There is no

error estimate available for the Euler equations_ so they used an error indicator Ej on

element j consisting of

Ej _ Z_j r_ _ r_dV _

_

Xk__

_____

Z__j

_j _ _nbj_k _dS_____

_

This can be shown ____ to be the length of the interface that separates the two _uids

on  j _ Remacle et al_ ____ increased the degree on elements where Ej was above the

median of all error indicators_ Results using this adaptive p_re_nement strategy with p

ranging from _ to _ are shown in Figure ____ The mesh used for these computations was

__ Hyperbolic Problems

Figure ______ Densities for the Rayleigh_Taylor instability of Example _____ at t _ ___

and p _  to __ The mesh used for all computations is shown at the left_

a uniform bisection of each element of the mesh shown in Figure ___ into four elements_

Successive frames in Figure ___ show the selected values of p and the density _ at

t _ __   _ ____ and __ _ The computations show the complex series of bifurcations that

occur at the interface between the two _uids__

Example _______ Flaherty et al_ ____ solve a _ow problem for the three_dimensional Eu_

ler equations ______ in a tube containing a vent Figure ____ using a piecewise_constant

discontinuous Galerkin method_ A van Leer _ux vector splitting _____ _ ______ ____

was used to evaluate _uxes_ No limiting is necessary with a _rst_order method_ The main

tube initially had a supersonic _ow at a Mach number ratio of the speed of the _uid to

the speed of sound_ of _____ There was no _ow in the vent_ At time t _  a hypothetical

diaphragm between the main and vent cylinders is ruptured and the _ow expands into

the vent_ Flaherty et al_ citeFLS__ solve this problem using an adaptive h_re_nement

procedure_ They used the magnitude of density jumps across element boundaries as a

re_nement indicator_ Solutions for the Mach number at t _  and ___ are shown on the

left of Figure ___ for a portion of the problem domain_ The mesh used in each each case

_____ Multidimensional Discontinuous Galerkin Methods __

Figure ____     _ Density for the Rayleigh_Taylor instability of Example _____ at t _ __          _

____ and __    left to right_ obtained by adaptive p_re_nement_ The values of p used on

each element are shown in the _rst_ third_ and _fth frames with blue denoting p _ _ and

red denoting p _ __

is shown on the right of the _gure_

A shock forms on the downwind end of the vent tube and expansion forms on the

upwind end_ The mesh is largely concentrated in these regions where the rapid solution

changes occur_ The initial mesh consisted of ______ elements_ This rose to more than

__ elements during the adaptive enrichment_ This computation was done on __

processors of a parallel computer_ The coloring of the images on the right of Figure ___

indicates processor assignments_

The discontinuous Galerkin method is still evolving and many questions regarding _ux

evaluation_ limiting_ a posteriori error estimation_ the treatment of di_usive problems_

and its e_ciency relative to standard _nite element methods remain unanswered_

Problems

__ Construct a typical term in the mass matrix on the canonical element by integrating

Z _

_ Z ___

_

Nm__ __Nn__ __d_d_

using the basis of monomials ____       __

__ Use the monomial basis ____         _ and the Gram_Schmidt process of Figure _____

to construct an orthogonal basis on the canonical right triangle for polynomials of

_ Hyperbolic Problems

Figure ______ Mach contours left_ and adaptive meshes right_ used to solve the com_

pressible _ow problem of Example _____ at t _  top_ and t _ ___ bottom__

degree p _ _ or less_

Bibliography

___ M_ Abromowitz and I_A_ Stegun_ Handbook of Mathematical Functions_ volume                   

of Applied Mathematics Series_ National Bureau of Standards_ Gathersburg_ _____

___ S_ Adjerid_ K_D_ Devine_ J_E_ Flaherty_ and L_ Krivodonova_ A posteriori error esti_

mation for discontinuous Galerkin solutions of hyperbolic problems_ In preparation_

__

___ F_ Bassi and S_ Rebay_ A high_order accurate discontinuous _nite element method

for the numerical solution of the compressible navier_stokes equations_ Journal of

Computational Physics_ _______"____ _____

___ C_E_ Baumann and J_T_ Oden_ A discontinuous hp _nite element method for

convection_di_usion problems_ to appear_ _____

_          _ K_S_ Bey and J_T_ Oden_ hp_version discontinuous galerkin method for hyper_

bolic conservation laws_ Computer Methods in Applied Mechanics and Engineering_

_____  _"____ _____

___ K_S_ Bey_ J_T_ Oden_ and A_ Patra_ hp_version discontinuous galerkin method for hy_

perbolic conservation laws_ A parallel strategy_ International Journal of Numerical

Methods in Engineering_ _______"____ ___  _

___ K_S_ Bey_ J_T_ Oden_ and A_ Patra_ A parallel hp_adaptive discontinuous galerkin

method for hyperbolic conservation laws_ Applied Numerical Mathematics_ _____"

____ _____

___ R_ Biswas_ K_D_ Devine_ and J_E_ Flaherty_ Parallel adaptive _nite element methods

for conservation laws_ Applied Numerical Mathematics_ ____                       "____ _____

___ A_J_ Chorin_ Random choice solution of hyperbolic systems_ Journal of Computa_

tional Physics_ _          _          __"       ___ _____

__

__ Hyperbolic Problems

___ B_ Cockburn_ G_ Karniadakis_ and C__W_ Shu_ editors_ Discontinous Galerkin Meth_

ods Theory Computation and Applications_ volume __ of Lecture Notes in Compu_

tational Science and Engineering_ Berlin_ __ Springer_

____ B_ Cockburn_ S__Y_ Lin_ and C__W_ Shu_ TVB Runge_Kutta local projection discon_

tinuous _nite element method for conservation laws III_ One_dimensional systems_

Journal of Computational Physics_ ____"____ _____

____ B_ Cockburn and C__W_ Shu_ TVB Runge_Kutta local projection discontinuous

_nite element method for conservation laws II_ General framework_ Mathematics of

Computation_ _____"__        _ _____

____ K_ Devine and J_E_ Flaherty_ Parallel adaptive hp_re_nement techniques for conser_

vation laws_ Applied Numerical Mathematics_ _____"____ _____

____ K_ Ericksson and C_ Johnson_ Adaptive _nite element methods for parabolic prob_

lems I_ A linear model problem_ SIAM Journal on Numerical Analysis_ _____"___

_____

__        _ K_ Ericksson and C_ Johnson_ Adaptive _nite element methods for parabolic prob_

lems II_ Optimal error estimates in l_l_ and l_l__ SIAM Journal on Numerical

Analysis_ _____"___ ___       _

____ J_E_ Flaherty_ R_ Loy_ M_S_ Shephard_ B_K_ Szymanski_ J_ Teresco_ and L_ Ziantz_

Adaptive local re_nement with octree load_balancing for the parallel solution of

three_dimensional conservation laws_ Journal of Parallel and Distributed Computing_

______"_        __ _____

____ J_ Glimm_ Solutions in the large for nonlinear hyperbolic systems of equations_

Communications on Pure and Applied Mathematics_ ______"__        _ ___   _

____ S_K_ Godunov_ A _nite di_erence method for the numerical computation of dis_

continuous solutions of the equations of _uid dynamics_ Mat_ Sbornik__ ______"___

__        __

____ C_ Johnson_ Error estimates and adaptive time step control for a class of one step

methods for sti_ ordinary di_erential equations_ SIAM Journal on Numerical Anal_

ysis_ _ ___"____ _____

___ P_D_ Lax_ Hyperbolic Systems of Conservation Laws and the Mathematical Theory

of Shock Waves_ Regional Conference Series in Applied Mathematics_ No_ ___ SIAM_

Philadelphia_ _____

_____ Multidimensional Discontinuous Galerkin Methods __

____ W_H_ Reed and T_R_ Hill_ Triangular mesh methods for the neutron transport

equation_ Technical Report LA_UR________ Los Alamos Scienti_c Laboratory_ Los

Alamos_ _____

____ J__F_ Remacle_ J_E_ Flaherty_ and M_S_ Shephard_ Adaptive order discontinuous

galerkin methods_ In preparation_ __

____ P_L_ Roe_ Approximate Riemann solvers_ parameter vectors_ and di_erence schemes_

Journal of Computational Physics_ ____          _"____ _____

____ P_ Le Saint and P_ Raviart_ On a _nite element method for solving the newtron

transport equations_ In C_ de Boor_ editor_ Mathematical Aspects of Finite Elements

in Partial Di_erential Equations_ pages __"__   _ New York_ _____ Academic Press_

__        _ G_A_ Sod_ Numerical Methods in Fluid Dynamic_ Cambridge University Press_

Cambridge_ ___          _

____ J_L Steger and R_F_ Warming_ Flux vector splitting of the inviscid gasdynamic

equations with applications to _nite di_erence methods_ Journal of Computational

Physics_ _____"____ _____

____ B_ van Leer_ Flux_vector splitting gor the Euler equations_ Lecture Notes in Physics_

___      _"         ___ _____

____ M_F_ Wheeler_ An elliptic collocation__nite element method with interior penalties_

SIAM Journal on Numerical Analysis_ _         __        _"____ _____

Навигация