Bubble and multiscale stabilization of bilinear finite element methods for transient advection–diffusion equations on rectangular grids

Journal of Computational and Applied Mathematics (

)

–

Contents lists available at SciVerse ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Bubble and multiscale stabilization of bilinear finite element methods for transient advection–diffusion equations on rectangular grids Onur Baysal ∗ İzmir University, Mathematics and Computer Science Department, İzmir PO: 35350, Turkey

article

info

Article history: Received 21 December 2012 Received in revised form 12 March 2013 Keywords: Stabilized FEM Convection–diffusion equation Parabolic problem Multiscale function

abstract In this paper a Petrov–Galerkin type stabilization for a time dependent advection–diffusion equation is considered. We first enrich the bilinear test space with bubble functions and the bilinear trial space with a special combination of bubble and multiscale functions for the steady state advection–diffusion equation. It is known that solving the residual equation obtained by the bubble elimination procedure is as difficult as solving the steady case of the original problem, which makes the enriched methods quite costly for two-dimensional problems. In this study, we instead utilize their cheap and efficient approximations in each rectangular element. Then we suggest a recipe for a proper adaptation of these functions combined with the generalized Euler time integration to the unsteady problem. Some numerical experiments for typical model problems are presented to illustrate the accuracy of our method. © 2013 Elsevier B.V. All rights reserved.

1. Introduction It is known that the standard Galerkin finite element method (SGFEM) based on low order piecewise polynomials is unsuitable for the solution of singularly perturbed problems. In the case where the advection term dominates over the diffusion one, numerical solutions obtained by SGFEM suffer from nonphysical oscillations, unless appropriately designed meshes are used. Although a number of studies focus on the steady problems, little attention is given to the unsteady variants. For unsteady problems, authors generally consider methods which are based on separating spatial and temporal discretization. One of the most favored and efficient approaches consists of discretizing the problem first in space by using the streamline-upwind Petrov–Galerkin (SUPG) method proposed in [1], then in time by means of generalized Euler time integration (the θ -method). Important recent references concerning the analysis of the SUPG method for the unsteady equations are [2] for the pure advection equation, [3,4] for the advection–diffusion equation and [5–7] for the advection–diffusion–reaction equation. This method includes an additional consistent term in the streamline direction that provides an additional diffusion. This additional diffusion is controlled by a stabilization parameter which should be selected in a proper way, but its optimum value is not precisely known a priori. This is usually seen as a drawback of the method. Another popular strategy is based on enriched methods which consist of the Galerkin method with enhanced approximation spaces followed by the application of simple time integration to the resulting semidiscrete formulation [8–10]. In this paper we take up this issue. For this purpose we first concentrate on the steady equation with the stabilized finite element method recently introduced by Franca et al. in [11,12] and analyzed in [13]. More precisely, this method is related to the residual-free bubble (RFB) method introduced in [14] (see also [15–21]). The RFB method improves the

∗

Tel.: +90 5426630827. E-mail addresses: [email protected], [email protected].

0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.03.035

2

O. Baysal / Journal of Computational and Applied Mathematics (

)

–

Galerkin method by enriching the polynomial spaces with bubble functions which satisfy the steady differential operator condition inside each element with the homogeneous Dirichlet boundary condition. But a most widely known drawback of this method is that it produces nonphysical oscillations near the outflow boundaries. To remedy this problem, Franca et al. suggest the multiscale (MS) method for steady advection–diffusion–reaction equations in [11]. This method is achieved by enriching the test function space with bubble functions and the trial space with the so-called MS functions. Indeed, the MS functions satisfy the same differential equation as the RFB functions, but they, in contrast, satisfy a suitable restriction of the steady differential operator at the element boundary instead of vanishing. However, in contrast to the RFB method case, internal layers are not well captured by the MS algorithm if the mesh is not aligned with the advection field. Therefore, motivated by this observation, Franca and his co-workers presented the MIX method in [11]. According to this algorithm, the MS functions are used in elements connected to the outflow boundaries and the RFB ones in the rest of the domain. This nonconforming MIX method is successful for all regimes, but costly, as it requires the evaluation of the RFB and the MS functions. However, the authors offer in [22,23] pseudo-approximation techniques for these functions whereby for each element a few suitable single points are used in order to retain the stabilizing property of the enriching functions. It is reported in [23] that utilizing these techniques, the MIX algorithm becomes more practical and its outcomes are sufficiently close to the ones in [11] for the steady state equation. The main subject of this paper is how the stabilized finite element methods RFB, MS and MIX designed for steady problems can be properly combined with generalized Euler time integration for the numerical solution of the unsteady advection–diffusion equation. The efficiency of pseudo-approximation techniques is also discussed for the RFB and the MS methods. The outline of this paper is as follows. In Section 2, we present the steady advection–diffusion model and the enriched methods. In Section 3, pseudo-approximation techniques are summarized for the RFB and the MS functions. Adaptation of these functions combined with the generalized Euler time integration to the unsteady problem is considered in Section 4. Finally numerical results and conclusions are given in Section 5. 2. The steady state advection–diffusion problem Let Ω be a bounded open domain in R2 with polygonal boundary ∂ Ω . We consider the following linear advection–diffusion problem:

Lu := −ϵ△u + β · ∇ u = f u=0

in Ω

(1)

on ∂ Ω

where ϵ > 0 is a constant diffusion coefficient. Let Th be a standard partition of the domain into rectangles K where h refers to the level of refinement of the discretization, and is defined by h := maxK ∈Th hK where hK is the diameter of K . Then we assume that the source term f and the convection field β = (β1 , β2 ) are piecewise constants with respect to the decomposition Th . We also assume that the components of β are positive in each element K . The outflow boundary is the part of ∂ Ω defined by ∂ Ω out := {(x, y) ∈ ∂ Ω | β · n(x, y) > 0} where n is the outflow normal to the boundary. We denote by Tout h the set of elements in Th each of which has at least one boundary contained by ∂ Ω out . The associated weak formulation is finding u ∈ H01 (Ω ) such that a(u, v) := ϵ(∇ u, ∇v) + (β · ∇ u, v) = (f , v) ∀v ∈ H01 (Ω ).

(2)

As usual, (·, ·)D denotes the L inner product on the open subset D of Ω . To simplify the notation, we just write (·, ·) when D = Ω . Due to the coercivity of the bounded bilinear form a(·, ·) over H01 (Ω ) and the Lax–Milgram theorem, the weak problem (2) is well-posed. To define the Galerkin approach, we introduce the finite dimensional subspace Vl (Ω ) of H01 (Ω ) as 2

Vl (Ω ) = {v ∈ H01 (Ω ) : v|K is a bilinear polynomial ∀K ∈ Th }. Then the Galerkin approach for (2) is given as follows: Find ul ∈ Vl ⊂ H01 (Ω ) such that a(ul , vl ) = (f , vl ) ∀vl ∈ Vl .

(3)

Here ul is the bilinear approximation of u and it can be represented by the linear combination of the standard nodal basis functions ψi with the coefficients ui . It is well known that the Galerkin method inherits the stability of the continuous problem and it yields spurious oscillations when the advection coefficient is larger than the diffusion one (ϵ ≪ |β|h). Since we are interested in finding a finite element discretization for (2) that is stable and coarse mesh accurate for all ϵ and β , we consider a Petrov–Galerkin type stabilization, so we respectively enrich the trial and test spaces as Uh (Ω ) := Vl (Ω ) ⊕ Eh (Ω ) and

Wh (Ω ) := Vl (Ω ) ⊕ B(Ω ).

Here B(Ω ) is the bubble space defined by B(Ω ) := {v ∈ H01 (Ω ) : v|K ∈ H01 (K ) ∀K ∈ Th }

(4)

O. Baysal / Journal of Computational and Applied Mathematics (

)

–

3

and we later define the enriching space Eh (Ω ). Now the Petrov–Galerkin problem (3) reads: Find uh ∈ U h (Ω ) such that a(uh , wh ) = (f , wh ) ∀wh ∈ Wh (Ω ).

(5)

Since a typical member uh of Uh (Ω ) can be split into a bilinear part ul ∈ Vl (Ω ) and an enriching part ue ∈ Eh (Ω ), solving (5) is equivalent to finding ul ∈ Vl (Ω ) such that a(ul + ue , vl ) = (f , vl ) ∀vl ∈ Vl (Ω )

(6)

where for all K ∈ TK , ue is the weak solution of the following residual equation:

Lue = f − Lul in K ue = µ

(7)

on ∂ K .

In order to evaluate ue uniquely, we need to choose an appropriate boundary condition. It is known that for the RFB method we set µ = µb := 0. In this case, the enriching space Eh (Ω ) can be represented as a direct sum of the two-dimensional bubble spaces [24], such that Eh (Ω ) = Bh (Ω ) :=



Bh (K ).

K ∈Th

Here Bh (K ) = span{bK0 , bK2 } and bK0 , bK1 are the solutions of the following problems:



LbK0 = 1 bK0

in K on ∂ K

=0

 and

LbK1 = β1 (x − xK ) + β2 (y − yK )

bK1 = 0

in K on ∂ K

(8)

where (xK , yK ) is the barycenter of K . Then the enriching part of the solution can be written as ue = ub :=



(α0K bK0 + α1K bK1 )

(9)

K ∈Th

∂2u

where α0K = (f − β · ∇ ul )(xK , yK ) and α1K = − ∂ x∂ yl . For the MS method, we choose µ = µm satisfying the following ordinary differential equation:

L∂ K µm = −L∂ K ul on ∂ K µm = 0 at the vertices of K

(10)

where

L∂ K v := −ϵ

∂v ∂ 2v + P(β, s) ∂ s2 ∂s

with s is a variable that parameterizes ∂ K by arc-length and P(β, s) is the usual projection of the convection field onto ∂ K . In order to make the method more practical, the boundary condition in the original MS algorithm in [11] was modified as in (10), but numerical tests indicate that this difference can be negligible. Let I0 be the set of indexes of internal nodal points with respect to the discretization of Ω ; then the MS space can be given by Eh (Ω ) = Mh (Ω ) := span{φi , φf }i∈I0 where φi and φf are enriching basis functions, defined by the following auxiliary problems:



Lφi = −Lψi

in K on ∂ K

φi = νi

L∂ K νi = −L∂ K ψi νi = 0

 where

on ∂ K at the vertices of K

(11)

and

φf =



f |K bK0 .

(12)

K ∈Th

Then enriching part of the solution can be written in terms of the global MS bases functions ue = um :=



ui φi + φf .

(13)

i∈I0

It is reported in [11] that the MS method captures layers near the outflow boundary better than the RFB method does. On the other hand the RFB method produces more accurate results than the MS one in some parts of the domain close to the internal layer if the mesh is not aligned with the advection field. On the basis of this important observation, the MIX algorithm (RFB–MS) was proposed by Franca et al.. The idea of this method can be simply explained by the rule that MS

4

O. Baysal / Journal of Computational and Applied Mathematics (

)

–

Fig. 1. A typical element K (left) and the approximate bubble basis function bK0 (right).

functions are used in the elements adjacent to the outflow boundaries and the RFB functions are used in the rest of the domain. According to MIX method, we set the boundary condition of the residual equation (7) µ = µbm such that

µbm =



µm , µb ,

K ∈ Tout h otherwise.

In this case the enriching part of the solution can be represented as ue = ubm :=



ub | K +

K ∈Th /Tout h



um |K .

(14)

K ∈Tout h

3. Computing the enriching functions In order to obtain the enriching part of the solution ue , Eqs. (8) and (11) should be solved. This task is as difficult as solving the original problem, which makes the method impractical. So in this study we employ the practical and reliable approximations of enriching functions which are proposed in [22] for the bubble and in [23] for MS functions. We first consider the bubble equation in (8). Although Bh (K ) is spanned by two bubble basis functions bK0 and bK1 in each element, considering only the dominant one bK0 is sufficient for reducing the instability. The pseudo-RFB approximation is also utilized for bK0 , as illustrated in Fig. 1, right. Further information about the shape of this pyramid and some stability properties can be found in [22]. We secondly summarize the pseudo-MS approximation for the MS basis function φi presented in [23]. Without loss of generality, we may consider the sample element K which is shown in Fig. 1, left, and assume that (hx , hy ) is the ith inner node in the discretization. Let us rewrite the problem (11) in terms of the enriched basis function λi := φi + ψi :



Lλi = 0

λi = θi

in K on ∂ K

 where

L∂ K θi = 0

θi = ψi

on ∂ K at the vertices of K .

(15)

By using the separation of variables method, the solution of (15) can be written as

λi (x, y) = λxi (x)λyi (y) where λxi (x) := θi |z1 and λyi (y) := θi |z4

(16)

y i

and the single-variable functions λxi and λ are the solutions of the following two-point boundary value problems:

  

L∂ K |z1 λxi = −ϵ

λxi (0) = 0

d2 λxi dx2

+ β1

dλxi

dx and λxi (hx ) = 1

= 0 in (0, hx )

(17)

and

 y 2 y L | λy = −ϵ d λi + β dλi = 0 in (0, h ) ∂ K z4 i 2 y dy2 dy  y y λi (0) = 0 and λi (hy ) = 1.

(18)

These homogeneous equations, (17) and (18), with nonhomogeneous boundary conditions, can be easily transformed to the form of bubble equations (nonhomogeneous equations with homogeneous boundary conditions); then a sub-grid point for y each λxi and λi can be chosen to construct their approximations. In order to keep their stabilizing property, the selection of these points is very important. In the literature, two different approaches are presented, specifically in [25,26]. The following

O. Baysal / Journal of Computational and Applied Mathematics (

0.1

0.1

0.05

0.05

0

0

−0.05

−0.05

−0.1

−0.1

0.1 0.08 0.06 0.04 0.02

0.02

0 0

0.04

0.06

0.08

0.1

)

0.1 0.08 0.06 0.04 0.02 0 0

–

0.02

5

0.04

0.06

0.08

0.1

Fig. 2. φi (left) and φ˜ i (right) for ϵ = 0.1, β = (1, 2) and hx = hy = 0.05.

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

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

0.1 0.08 0.06 0.04 0.02

0.02

0 0

0.04

0.06

0.08

0.1

0.1 0.08 0.06 0.04 0.02 0 0

0.02

0.04

0.06

0.08

0.1

Fig. 3. φi (left) and φ˜ i (right) for ϵ = 0.001, β = (1, 2) and hx = hy = 0.05. y

˜ xi , is obtained by direct application of the idea in [25] to the function λxi (λ˜ i (y) can be constructed in a similar approximation, λ way):  x β1 ξx x   + ,  2h ϵ h x x λ˜ xi (x) = β1 ηx (hx − x) x    + , 2hx ϵ hx

x ≤ ηx x > ηx

 2ϵ   hx − , β1 where ηx =  h  x, 2

ϵ≤

β1 hx 4

(19)

otherwise.

and ξx = hx − ηx . Thus, the approximate enriching basis function φ˜ i is evaluated by using

φ˜ i (x, y) = λ˜ i (x, y) − ψi (x, y) = λ˜ xi (x)λ˜ yi (y) − ψi (x, y).

(20)

In Figs. 2 and 3, a comparison of φi and its approximation φ˜ i are presented on a patch of four rectangular elements for the diffusion dominated case and the advection dominated case, respectively. As we see in these figures, although we have used just a single additional node in each element, the approximations are very close to the exact forms. Hence we utilize the approximate enriching functions φ˜ i instead of φi in the related algorithms. 4. An evolutionary problem In this part we consider the following initial–boundary value problem:

Lt u := ut + Lu = f u=0

in Ωt := Ω × (0, T )

on ∂ Ω × [0, T ],

u(·, 0) = u0 (·) in Ω

(21)

6

O. Baysal / Journal of Computational and Applied Mathematics (

)

–

where L is taken as in the previous section, u0 ∈ L2 (Ω ) is an initial datum and the right hand side function f (t ) is chosen from L2 (Ω ) for each t ∈ [0, T ]. Then the weak formulation associated with (21) reads: For all t ∈ (0, T ], find u(t ) ∈ H01 (Ω ), u(0) = u(·, 0) satisfying d dt

(u(t ), v) + a(u(t ), v) = (f (t ), v) ∀v ∈ H01 (Ω ).

(22)

In order to construct an approximation to (21), we first discretize the space variable. This process leads to a system of ordinary differential equations whose solution uh (t ) is an approximation of the exact solution u at t ∈ [0, T ]. Then the semidiscrete Petrov–Galerkin approximation of (21) can be given as follows: For all t ∈ (0, T ], find uh (t ) ∈ Uh such that d dt

(uh (t ), vl ) + a(uh (t ), vl ) = (f (t ), vl ) ∀vl ∈ Vl (Ω ).

(23)

To simplify the algorithm we enrich only the trial space, as depicted in (4), for each t ∈ (0, T ]. However, the test space will be taken as Wh = Vl and the enriched part of the solution of the steady equation will be adapted for the unsteady case by adopting the following recipe. RFB method: ue (t ) = ub (t ) :=

 (f (t ) − Lul (t ))|(xK ,yK ) bK0 .  K ∈Th MS method: ue (t ) = um (t ) := i∈I0 ui (t )φi + φf (t ).   MIX method: ue (t ) = ubm (t ) := ub (t )|K + K ∈Tout ue (t )|K . K ∈Th /Tout h h  K Here φf (t ) = K ∈Th f (t )|K b0 and ui (t ) is the unknown coefficient of the ul (t ) such that  ui (t )ψi (x, y). ul (t ) = i∈I0

In order to obtain full discretization of (21), we take uniform subintervals for the time variable with t n = nτ , n = 0, . . . , N. Here τ is the time step and N = T /τ . Next we replace the time derivative with suitable difference quotients in (23); therefore we construct a sequence unh (x, y) which is the approximation of u(x, y, t n ). For simplicity, we restrict ourselves to the one-step-scheme θ method [27]. Applying this method to Eq. (23), we get the following discrete problem: Given u0h as some suitable approximation of u(0), for n ≥ 1, find unh ∈ Uh such that



unh − uhn−1

τ

 , vl + a(unh−θ , vl ) = (f n−θ , vl ) ∀vl ∈ Vl

(24)

where unh−θ = θ unh−1 + (1 − θ )unh . The values θ = 0, θ = 1/2, and θ = 1 refer to the backward Euler, Crank–Nicolson and forward Euler methods, respectively. Since unh = unl + une , the full discrete algorithm (24) can be given: For all t ∈ (0, T ], find ul (t ) ∈ Vl such that



unl − uln−1

τ

+

une − une −1

τ

 , vl + a(unl −θ + une −θ , vl ) = (f n−θ , vl ) ∀vl ∈ Vl .

(25)

Here une denotes the discrete version of the enriched function ue (t ).

 n K n K ∈Th (f − Lul )|(xK ,yK ) b0 .  n n n n MS method: ue = um := i∈I0 ui φi + φf .   n n MIX method: une = unbm := K ∈Th /Tout ub |K + K ∈Tout ue |K .

RFB method: une = unb :=

h

h

We also employ the functions φ˜ i and b˜ K0 instead of their exact counterparts in the full discrete equation (25). 5. Numerical experiments and the conclusion We illustrate the accuracy of the method by means of two different test problems on Ω = (0, 1)2 with uniform 20 ∗ 20 rectangular discretization. For all numerical simulations, ϵ = 10−6 and τ = 0.0025 are used. For the first test problem, we chose f = 0, β = (1, 1), θ = 0 with L-shape discrete type initial condition (see Fig. 4, left and center). The unstable result for the standard Galerkin solution at T = 1/8 has been shown in Fig. 4, right. We also compared the enriching finite element method RFB, MS and MIX algorithms in Figs. 5 and 6 at the final times T = 1/8 and T = 1/4, respectively. It is known that the MS method captures the outflow layer and the RFB methods capture the internal layer well for the steady problem. In these figures, similar results can be seen for the unsteady problem, verifying the accuracy of the MIX algorithm.

O. Baysal / Journal of Computational and Applied Mathematics (

)

–

7

Initial Data

Initial Data

1

Std. Gal.

1

1

0.8

0.8

0.6

0.6

0.8

0.4 0.2

1

0 0 0.5 x

y

y

0.6

0.4

0.4

0.2

0.2

0.5 y

1 0

0

0

0.2

0.4

x

0.6

0.8

0

1

0

0.2

0.4

x

0.6

0.8

1

Fig. 4. Initial condition (left), its counter-plot (center) and the Galerkin solution at T = 1/4 (right).

1

1

1

0.5

0.5

0.5 1

1 0 0 0.5 x

0.5 y

1 0

RFB

1

0 0

0.5 y

0.5 x

1 0 0

MS

1 0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0.4

0.6

0.8

1

0 0

0.2

0.4

1 0 MIX

1

0.8

0 0

0.5 y

0.5 x

1 0

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

Fig. 5. Results obtained by RFB (left), MS (center) and MIX (right) algorithms at T = 1/8.

In the second test problem, we examine the enriching methods with the problem parameters f = 1, β = (1, 1/2), T = 1/2 and homogeneous initial condition. We again compared the RFB, MS and MIX methods in Figs. 7 and 8 for backward Euler (θ = 0) and Crank–Nicolson (θ = 0.5) time discretization, respectively. In both cases, the MIX algorithm shows more stable behaviors. These numerical simulations also validate the accuracy of the MIX algorithm combined with the implicit Euler time integration for the unsteady advection–diffusion problem. Therefore, this combination algorithm can be considered as not only reliable but also a very practical method. In conclusion, it should be emphasized that since the enriching basis functions φi and bK0 are obtained from the steady equation and applied directly to the unsteady variant, their shapes do not change at different time levels. This situation makes the methods quite practical, but the result is not entirely consistent with the steady case considered in [11,23]. For example Franca and his co-workers reported in [11] that if the mesh is aligned with the advection field, the MS algorithm gives a nearly exact solution of the problem. Similar results were found in the work [23] by using the approximate forms of the enriching functions. However, the computational experiments indicate that this idea is no longer valid for the unsteady problem. We believe that the construction of time dependent MS functions conserves this property of the MS method, so

8

O. Baysal / Journal of Computational and Applied Mathematics (

)

1

1

1

0.5

0.5

0.5 1

1 0 0

0 0

0.5 y

0.5 x RFB

1 0 0

0.5 y

0.5 x

1 0

1

–

1 0 MS

1 0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

1 0 MIX

1

0.8

0 0

0.5 y

0.5 x

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

Fig. 6. Results obtained by RFB (left), MS (center) and MIX (right) algorithms at T = 1/4.

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 0.4 y 0.2

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 1 0.8 0.4 0.6 0.4 y 0.2 0.2 x 0 0

0 0

RFB

1

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 1 0.8 0.4 0.6 y 0.2 0.4 0.2 0 0 x MS

1

0.9

0.9

0.9

0.8

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0 0

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0 0

0.6

0.8

MIX

1

0.8

0.4 0.2 x

0.2

Fig. 7. Results obtained by RFB (left), MS (center) and MIX (right) algorithms for θ = 0.

0.4

0.6

0.8

1

1

O. Baysal / Journal of Computational and Applied Mathematics (

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 0.4 y 0.2

0 0

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 1 0.8 0.4 0.6 0.4 y 0.2 0.2 x 0 0 RFB

1

–

0.9

0.9

0.9

0.8

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0 0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0

0.4 0.2 x

0.6

0.8

1

MIX

1

0.8

0 0

9

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 1 0.8 0.4 0.6 0.4 y 0.2 0.2 0 0 x MS

1

)

0.2

0.4

0.6

0.8

1

Fig. 8. Results obtained by RFB (left), MS (center) and MIX (right) algorithms at θ = 1/2.

possibly more stable results are obtained. But in this case the applicability of the methods may be reduced. This issue requires additional advanced research and will be a subject of a further publication. References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259. [2] E. Burman, Consistent SUPG-method for transient transport problems: stability and convergence, Comput. Methods Appl. Mech. Engrg. 199 (2010) 1114–1123. [3] E. Burman, G. Smith, Analysis of the space semi-discretized SUPG method for transient convection–diffusion equations, Math. Models Methods Appl. Sci. 21 (10) (2011) 2049–2068. [4] P.B. Bochev, M.D. Gunzburger, J.N. Shadid, Stability of the SUPG finite element method for transient advection diffusion problems, Comput. Methods Appl. Mech. Engrg. 193 (2004) 2301–2323. [5] G. Lube, D. Weiss, Stabilized finite element methods for singularly perturbed parabolic problems, Appl. Numer. Math. 17 (1995) 431–459. [6] G. Hauke, M.H. Doweidar, Fourier analysis of semidiscrete and space–time stabilized methods for the advective–diffusive–reactive equation. I. SUPG, Comput. Methods Appl. Mech. Engrg. 194 (1) (2005) 45–81. [7] V. John, J. Novo, Error analysis of the SUPG finite element discretization of evolutionary convection–diffusion–reaction equations, SIAM J. Numer. Anal. 49 (3) (2011) 1149–1176. [8] L.P. Franca, J.V.A. Ramalho, F. Valentin, Enriched finite element methods for unsteady reaction–diffusion problems, Comm. Numer. Methods Engrg. 22 (6) (2006) 619–625. [9] M.I. Asensio, B. Ayuso, G. Sangalli, Coupling stabilized finite element methods with finite difference time integration for advection–diffusion–reaction problems, Comput. Methods Appl. Mech. Engrg. 196 (2007) 3475–3491. [10] J. Frutos, J. Novo, Bubble stabilization of linear finite element methods for nonlinear evolutionary convection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 197 (2008) 3988–3999. [11] L.P. Franca, J.V.A. Ramalho, F. Valentin, Multiscale and residual-free bubble functions for reaction–advection–diffusion problems, Int. J. Multiscale Comput. Eng. 3 (3) (2005) 297–312. [12] L.P. Franca, A.L. Madureira, F. Valentin, Towards multiscale functions: enriching finite element spaces with local but not bubble-like functions, Comput. Methods Appl. Mech. Engrg. 194 (2005) 3006–3021. [13] L.P. Franca, A.L. Madureira, L. Tobiska, F. Valentin, Convergence analysis of a multiscale finite element method for singularly perturbed problems, Multiscale Model. Simul. 4 (3) (2005) 839–866. [14] F. Brezzi, A. Russo, Choosing bubbles for advection–diffusion problems, Math. Models Methods Appl. Sci. 4 (1994) 571–587. [15] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residual-free bubbles for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 166 (1998) 25–33. [16] F. Brezzi, T.J.D. Hughes, L.D. Marini, A. Russo, E. Süli, A priori error analysis of residual-free bubbles for advection–diffusion problems, SIAM J. Numer. Anal. 36 (1999) 1933–1948. [17] F. Brezzi, D. Marini, A. Russo, Applications of the pseudo residual-free bubbles to the stabilization of convection diffusion problems, Comput. Methods Appl. Mech. Engrg. 166 (1998) 51–63. [18] F. Brezzi, D. Marini, E. Süli, Residual-free bubbles for advection–diffusion problems: the general error analysis, Numer. Math. 85 (2000) 31–47.

10

O. Baysal / Journal of Computational and Applied Mathematics (

)

–

[19] L.P. Franca, A.I. Neslitürk, M. Stynes, On the stability of residual-free bubbles for convection–diffusion problems and their approximation by two-level finite element method, Comput. Methods Appl. Mech. Engrg. 166 (1998) 35–49. [20] G. Sangalli, Global and local error analysis for the residual-free bubbles method applied to advection-dominated problems, SIAM J. Numer. Anal. 38 (2000) 1496–1522. [21] G. Sangalli, Capturing small scales in elliptic problems using a residual-free bubbles finite element method, Multiscale Model. Simul. 1 (2003) 485–503. [22] A.I. Neslitürk, On the choice of stabilizing subgrid for convection–diffusion problem on rectangular grids, Comput. Math. Appl. 59 (2010) 3687–3699. [23] A.I. Neslitürk, O. Baysal, Pseudo multiscale functions for the stabilization of convection–diffusion equations on rectangular grids, Int. J. Multiscale Comput. Eng. (2012), forthcoming (http://dx.doi.org/10.1615/IntJMultCompEng.2012004234). [24] L.P. Franca, L. Tobiska, Stability of the residual free bubble method for bilinear finite elements on rectangular grids, IMA J. Numer. Anal. 22 (2002) 73–87. [25] F. Brezzi, G. Hauke, L.D. Marini, G. Sangalli, Link-cutting bubbles for the stabilization of convection–diffusion–reaction problems, Math. Models Methods Appl. Sci. 13 (2003) 445–461. [26] A.I. Neslitürk, A stabilizing subgrid for convection–diffusion problem, Math. Models Methods Appl. Sci. 16 (2) (2006) 211–231. [27] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer, 1996.