Simulating transient phenomena via residual free bubbles

Comput. Methods Appl. Mech. Engrg. 200 (2011) 2127–2130

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Simulating transient phenomena via residual free bubbles A.L.G.A. Coutinho a, L.P. Franca b,⇑, F. Valentin c a

Center for Parallel Computing and Department of Civil Engineering, COPPE/Federal University of Rio de Janeiro, P.O. Box 68506, Rio de Janeiro, RJ 21945, Brazil University of Colorado Denver, Department of Mathematical and Statistical Sciences, P.O. Box 173364, Campus Box 170, Denver, CO 80217-3364, United States c National Laboratory for Scientific Computing – LNCC, Av. Getúlio Vargas, 333, 25651-070 Petrópolis, RJ, Brazil b

a r t i c l e

i n f o

Article history: Received 9 June 2010 Received in revised form 17 November 2010 Accepted 19 February 2011 Available online 27 February 2011

a b s t r a c t We derive two stabilized methods for transient equations using static condensation of residual-free bubbles. The methods enhance the stability of the Discontinuous Galerkin method. Ó 2011 Elsevier B.V. All rights reserved.

Keywords: Residual-free bubbles Transient phenomena Artificial viscosity in time Hyperbolic problem

1. Introduction Time dependent problems are generally discretized using finite elements in space and finite difference methods in time. This is termed a semidiscrete formulation. An exception to this approach is to fully discretize, i.e., use finite elements both in time and space. In order to make it feasible Discontinuous Galerkin method is used in time. In this paper we start with the Discontinuous Galerkin method in time for a couple of model problems using piecewise linear approximations in time and in space, then we enrich the trial and test functions using residual-free bubbles [1,2] and we use static condensation to derive stabilized methods. We would like to point out that this elimination is only possible because we add a temporal artificial diffusion that is taken to be zero at the end of the derivation of the methods. Temporal artificial diffusion (TAD) is not a common notion as spatial artificial diffusion (SAD). TAD has been introduced as an elliptic regularization in time [4,3]. In the next section we use this idea for an initial value problem and in the subsequent section we apply it to a purely advective problem.

Consider the IVP given by: find u(t) such that

⇑ Corresponding author. Tel.: +55 21 2535 1993. E-mail address: [email protected] (L.P. Franca). 0045-7825/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2011.02.016

uðtÞ ¼ u0 þ

ð1Þ

Z

t

f ðsÞds

ð2Þ

0

Herein we are interested in numerical approximations of (1). IVPs are generally discretized by finite difference methods such as Euler and Runge–Kutta methods (see a numerical analysis textbook for other finite difference methods). We now consider the Discontinuous Galerkin method for the IVP: Find uh(t), t 2 In such that

Bðuh ; v Þn ¼ Lðv Þn ;

n ¼ 0; 1; . . . ; N  1

ð3Þ

where

Bðuh ; v Þn ¼

Lðv Þn ¼

2. Initial value problem

u;t ¼ f ðtÞ in ð0; T uð0Þ ¼ u0 ;

where f(t) is a piecewise constant function defined in a partition of (0, T] into subintervals In = (tn, tn+1), n = 0, 1, . . . , N  1 with tN = T. The solution to (1) is straightforward and it follows from the Fundamental Theorem of Calculus as

Z

Z

t nþ1 tn

t nþ1

tn

v ðtÞuh;t ðtÞdt þ v

v ðtÞf ðtÞdt þ v

 þ h þ tn u tn

 þ h  tn u tn

ð4Þ

ð5Þ

  Here we used the notation v t n ¼ lim!0 v ðt n  Þ. The method implies satisfaction of the differential equation plus continuity weakly on each time interval. For any order of interpolation we can show that methods of this type are all unconditionally stable and high order accurate.

2128

A.L.G.A. Coutinho et al. / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2127–2130

Error estimates can be derived in the norm

1 2

jjjv jjj2 ¼

N1 X

1 ½v ðt n Þ2 þ ðv ðT  Þ2 þ v ð0þ Þ2 Þ 2 n¼1

ð6Þ

Lub ¼ ðLu1  f Þ in In

ð17Þ

ub ¼ 0 on @In

ð18Þ

Here

where [] denotes the jump operator. Note that this norm is silent about the interior of the intervals. This led us to consider stabilization via enrichment functions.

Lu ¼ u;tt þ u;t

ð19Þ

Then

ub;tt þ ub;t ¼ ðu1;t  f Þ in In

ð20Þ

ub ¼ 0 on @In

ð21Þ

3. Enriching and deriving the stabilized method We start by adding an artificial viscosity in time written as

ðv ;t ; u;t ÞI

ð7Þ

n

where  > 0 is a small parameter and ð; ÞIn denotes integration over an interval In. This artificial diffusion parameter will be set to zero after we derive the stabilized method. We rewrite the Discontinuous Galerkin method as:



 v ;t ; uh;t



þ

In



v ; uh;t

 In

        þ v t þn uh t þn ¼ ðv ; f ÞIn þ v t þn uh t n

ð8Þ

Our functions uh and v are discretized by adding a bubble function to a piecewise linear function, i.e.,

v ¼ v1 þ vb

h

u ¼ u1 þ ub Taking

ðv b;t ; u1;t þ ub;t ÞI þ ðv b ; u1;t þ ub;t ÞI     ¼ ðv b ; f ÞI þ v b t þn ðu1 þ ub Þ t n n

n

    þ v b t þn ðu1 þ ub Þ t þn

  Using the zero boundary conditions for the bubble, i.e., ub tþ n ¼     þ ub t n ¼ v b ðt n Þ ¼ v b t nþ1 ¼ 0 and integration by parts we get for ub = cub and vb = ub

þ ðub ; u1;t ÞIn ¼ ðub ; f ÞIn

ub ¼ 0 on @ In

ð23Þ

Multiply by ub and integrate over the interval:

Z

In

ub ¼ ðub;t ; ub;t ÞIn þ ðub;t ; ub ÞIn

ð24Þ

The last term is zero by integration by parts, leaving us with

Z In

ub

ð25Þ

s¼

R

R

ub Þ2 ub ¼ In jIn j jIn jkub;t k2In ð

In

So we need

R

In

ð26Þ

ub to define the method. Use  ? 0 in (22)

~ b;t ¼ 1 in In u

ð27Þ

~ b ¼ 0 on t þn u

ð28Þ

ð10Þ

Solving for c:

c¼

ð22Þ

Substituting in the expression of s we get:

ð9Þ

n

n

ub;tt þ ub;t ¼ 1 in In

kub;t k2 ¼

v = vb we get

ckub;t k2I

For piecewise-constant f the solution is spanned by a single bubble basis function ub where

Then

1 2 b;t kIn

ku

ðf  u1;t ; ub ÞIn

We now take

ð11Þ

ð29Þ

and

v = v1 to get

ðv 1;t ; ðu1 þ ub Þ;t ÞI þ ðv 1 ; ðu1 þ ub Þ;t ÞI     ¼ ðv 1 ; f ÞI þ v 1 t þn ðu1 þ ub Þ t n n

n

R Dt

    þ v 1 tþn ðu1 þ ub Þ tþn

n

s¼ ð12Þ

Integrating by parts and using the zero boundary condition for the bubbles the formulation reduces to

ðv 1;t ; u1;t ÞI ¼ ðv 1 ; f ÞI

u~b ¼ t in ½0; Dt

    þ ðv 1 ; u1;t ÞIn  ðv 1;t ; ub ÞIn þ v 1 t þn u1 t þn     þ v 1 t þn u1 t n n

n

ð13Þ

0

tdt Dt ¼ 2 Dt

We can finally write the method (13) for

ð30Þ

 = 0 as follows:

    Dt ðv 1 ; u1;t ÞIn þ ðu1;t  f ; v 1;t ÞIn þ v 1 t þn u1 tþn 2     ¼ ðv 1 ; f ÞIn þ v 1 tþn u1 t n

ð31Þ

This is the Discontinuous Galerkin method for piecewise linears plus the additional term:

The method is equivalent to a stabilized method. There is no parameter to be selected. Its convergence now also depends on the term Dtkv 1;t k2In which is an interior term absent in the DG method.

 ðv 1;t ; ub ÞIn ¼ ðv 1;t ; cub ÞIn

4. Hyperbolic problem

¼ v 1;t jIn

1 2 b;t kIn

Z ðf  u1;t Þð ub Þ2

ku In 2 u Þ b In ðu1;t  u1;tt  f ; v 1;t ÞIn ¼ jIn jkub;t k2In ð

R

ð14Þ Consider the purely advective equation given by: find u(x, t) such that

ð15Þ

This is a perturbation term of a stabilized method with stabilization parameter given by:

s¼

R

ub Þ2 jIn jkub;t k2In ð

In

ð16Þ

This derivation works for any single bubble added to an interval In. We now select residual free bubbles to be these bubbles. The equations of residual-free bubbles are given by

u;t þ au;x ¼ f ðtÞ in Q ¼ ½0; L  ð0; T

ð32Þ

uð0; tÞ ¼ gðtÞ on ð0; T

ð33Þ

uðx; 0Þ ¼ hðxÞ on ½0; L

ð34Þ

where f(t) is a piecewise constant function defined in a partition Q into elements K = (xi1, xi)  (tn1, tn). The Discontinuous Galerkin (DG) method in time reads: Find uh(x, t) such that

2129

A.L.G.A. Coutinho et al. / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2127–2130

Bðuh ; v Þn ¼ Lðv Þn ;

n ¼ 0; 1; . . . ; N  1

ð35Þ

where

  Bðuh ; v Þn ¼ uh;t þ auh;x ; v

Qn

Lðv Þn ¼ ðf ; v ÞQ n þ

Z

þ

Z 0

L

    uh x; tþn v x; t þn dx

s¼ ð36Þ

L h

u

0

ðx; t n Þ

v

ðx; t þn Þdx

ð37Þ

where Qn = In  [0, L]. As in the IVP case, DG implies satisfaction of the differential equation plus continuity weakly on each time interval. We now introduce small viscosities in space and time to derive our method (SAD) and (TAD). At the end we set these viscosities to zero. The modified PDE reads:

u;t þ au;x ¼ u;tt þ ju;xx þ f

Thus we have a stabilization parameter s given by:

ð38Þ

R ð K ub Þ2

ð46Þ

ðkub;t k2K þ jkub;x k2K ÞjKj

This method holds for any single bubble added per element. Let us now examine what happens if we use residual-free bubbles. The residual-free bubbles equations are:

Lub ¼ ðLu1  f Þ in K

ð47Þ

ub ¼ 0 on@K Here

Lu ¼ u;tt  ju;xx þ u;t þ au;x

ð48Þ

Then

 ub;tt  jub;xx þ ub;t þ aub;x ¼ ðu1;t þ au1;x  f Þ in K

The DG method becomes:

ub ¼ 0 on @K

ðu;t ; v ;t ÞQ

Since the right-hand-side is piecewise constant then ub = cub

þ

Z

0

L

n

þ jðu;x ; v ;x ÞQ n þ ðu;t þ au;x ; v ÞQ n

    u x; t þn v x; tþn dx

¼ ðf ; v ÞQ n þ

Z

L 0

    u x; t n v x; tþn dx

ð39Þ

Take v = vb and use the zero boundary condition for the bubbles on each element boundary @K:

ðub;t ; v b;t ÞK þ jðub;x ; v b;x ÞK þ ðu1;t þ ub;t þ au1;x þ aub;x ; v b ÞK ¼ ðf ; v b ÞK Take

c¼

1

kub;t k2K þ jkub;x k2K

Now take

¼ ðf ; ub ÞK

ðf  u1;t  au1;x ; ub ÞK

Z

K

ub ¼ ðub;t ; ub;t ÞK þ jðub;x ; ub;x ÞK ¼ kub;t k2K þ jkub;x k2K R

R ð K ub Þ2 ðkub;t k2K þ jkub;x k2K ÞjKj

ð52Þ

¼

K

ub

jKj

ð53Þ

 ? 0 and j ? 0 in (50) we have

For

ð42Þ

~ b;x ¼ 1 in K ~ b;t þ au u ~ b ¼ 0 on @K  u

Z

Z

n

ð54Þ

K

L

0

    u1 x; t n v 1 x; t þn dx

0

h

^ a

ð55Þ

1þa

^ a

ð43Þ

where h is the longest segment in the triangle K parallel to the vec^ and tor a

1 3

h

^ a

s ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi

ð56Þ

1þa

Thus the method (44) for

0

Z

1 3

~ b ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi jKj u 2

0

þ jðu1;x ; v 1;x ÞQ n þ ðu1;t þ au1;x ; v 1 ÞQ n  ðub ; v 1;t Z L     þ u1 x; t þn v 1 x; t þn dx

¼ ðf ; v 1 ÞQ n þ

Multiply by ub and integrate over the element K

ð41Þ

Integrating by parts this simplifies to:

þ av 1;x ÞQ n

ð51Þ

s¼

þ jððu1 þ ub Þ;x ; v 1;x ÞQ n þ ððu1 þ ub Þ;t Z L     þ aðu1 þ ub Þ;x ; v 1 ÞQ n þ u1 x; tþn v 1 x; tþn dx

n

ub ¼ 0 on @K

^ ¼ ð1; aÞ starting The solution is a linear function in the direction a ^  n < 0g. Thus from zero at the inflow boundary @K  ¼ fx 2 @Kja

ððu1 þ ub Þ;t ; v 1;t ÞQ

ðu1;t ; v 1;t ÞQ

ð50Þ

Substituting in the expression for s:

v = v1

¼ ðf ; v 1 ÞQ n þ

ub;tt  jub;xx þ ub;t þ aub;x ¼ 1 in K

ð40Þ

vb = ub, ub = cub then

ckub;t k2K þ jckub;x k2K þ ðu1;t þ au;x ; ub ÞK

ð49Þ

 = j = 0 becomes ^ a

L



u1 x; t n



v



þ 1 x; t n



dx

ð44Þ

This is the DG method for piecewise linears plus a perturbation term given by

ðub ; v 1;t þ av 1;x ÞQ n ¼ cðub ; v 1;t þ av 1;x ÞK v 1;t þ av 1;x ¼ kub;t k2K þ jkub;x k2K Z  ðf  u1;t  au1;x ÞjK ð ub Þ2 K R ð K ub Þ2 ¼ ðkub;t k2K þ jkub;x k2K ÞjKj

1 h pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu1;t þ au1;x  f ; v 1;t þ av 1;x ÞQ n 3 1 þ a2 Z L   þ     þ u1 x; t n  u1 x; tn v 1 x; t þn dx

ðu1;t þ au1;x ; v 1 ÞQ n þ

0

¼ ðf ; v 1 ÞQ n

ð57Þ

This is a stabilized space–time formulation derived from residualfree bubbles. The derivation was only possible due to a Temporal Diffusive Artificial term. The s parameter is derived and therefore is not left to be chosen as in standard stabilized methods [5]. Acknowledgments

 ðu1;t þ au1;x  u1;tt  ju1;xx  f ; v 1;t þ av 1;x ÞK ð45Þ

During the course of this work L.P. Franca was a Visiting Professor at the Federal University of Rio de Janeiro. The authors thank the support of CNPq and Petrobras.

2130

A.L.G.A. Coutinho et al. / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2127–2130

References [1] F. Brezzi, A. Russo, Choosing bubbles for advection–diffusion problems, Math. Models Methods Appl. Sci. 4 (1994) 571–587. [2] L. Franca, A. Russo, Deriving upwinding, mass lumping and selective reduced integration by residual-free bubbles, Appl. Math. Lett. 9 (1996) 83–88.

[3] T. J. R. Hughes, G. Scovazzi, L. P. Franca, Multiscale and stabilized methods, in: Encyclopedia of Computational Mechanics, E. Stein, R. D. Borst, T. Hughes (Eds.), John Wiley & Sons, 2004. [4] T.J.R. Hughes, J. Stewart, A space-time formulation for multi-scale phenomena, J. Comput. Appl. Math. 74 (1996) 217–229. [5] C. Johnson, U. Navert, J. Pitkaranta, Finite element methods for linear hyperbolic problem, Comput. Methods Appl. Mech. Eng. 45 (1984) 285–312.