Multiscale finite element modeling of diffusion-reaction equation using bubble functions with bilinear and triangular elements

Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107 www.elsevier.com/locate/cma

Multiscale finite element modeling of diffusion-reaction equation using bubble functions with bilinear and triangular elements M. Parvazinia, V. Nassehi

*

Department of Chemical Engineering, Loughborough University, Loughborough, Leicestershire, LE11 3TU, UK Received 27 March 2006; received in revised form 8 August 2006; accepted 8 August 2006

Abstract A multiscale finite element method based on the bubble functions is used to generate stable and accurate solutions for the diffusionreaction equation for both dissipation and production cases. This equation shows multiscale behaviour when the process becomes reaction dominated. To validate the approach, the numerical results obtained for a benchmark problem are compared with the analytical solution. Further numerical experiments are carried out to investigate the performance of the developed technique over a range Damkohler number. The numerical experiments are performed using bilinear and triangular elements.  2006 Elsevier B.V. All rights reserved. Keywords: Diffusion-reaction equation; Bubble function; Multiscale; Finite element; Damkohler number

1. Introduction The numerical solution of diffusion-reaction equation becomes unstable using traditional methods where the reaction dominant case is considered. Therefore numerical solution of the equation is not a trivial matter. Well known examples of such difficulties are encountered in cases where the diffusion-reaction equation cannot be satisfied by a globally smooth function. In particular, if the field variables vary rapidly within thin layers, standard numerical schemes lead to inaccurate and unstable results [1]. However, almost all of these cases can be regarded as representation of multiscale phenomena in which both fine and coarse scale variations of the field variables need to be taken into account in the numerical solution scheme used. In theory any basically sound scheme should generate accurate numerical results if sufficiently fine computational grids are used. However in practice the simultaneous representation of all physical phenomena requires a very high level of discretisation which is not computationally cost *

Corresponding author. Tel.: +44 1509 222522; fax: +44 1509 223923. E-mail address: [email protected] (V. Nassehi).

0045-7825/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.08.011

effective [2]. Using a coarse discretisation the fine scale information is ignored and the solution becomes unstable and inaccurate. The development of algorithms which enable implementation of the theoretical concepts in the practical problems is not a trivial matter [3]. A theoretically rigorous numerical solution of the diffusion-reaction equation can be developed using the variational multiscale method proposed by Hughes [4] and Hughes and Stewart [5]. To develop a combined fine and coarse scale scheme, the field unknown (T) is divided into two parts as T = T1 + Tb, where Tb represents the fine scale variations of T which may be derived analytically, while T1 the coarse scale variations of T is approximated by the standard polynomial finite element discretisations. The use of bubble enriched elements offers a systematic approach to develop a combined discretisation. Bubble functions are, typically, high order polynomials which vanish on the element boundaries [6–9]. A systematic approach to derive bubble functions is the residual free bubble function (RFB) method [10–13]. The essential idea of this method is the use of the analytical solution of the governing differential equation within each element subject to homogeneous boundary conditions. This implies that the

1096

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

residual free bubble functions should satisfy, strongly, the governing equation within each element. However, in multi-dimensional problems, the analytical solution of governing equations, which are partial differential equations (PDE) can cause major difficulties. To overcome this difficulty a simple semi-discrete method has been proposed by Parvazinia et al. [14] in which the solution of a PDE is replaced by the analytical solution of ordinary differential equations. In this technique the exact solutions obtained from the ODE is expanded using the Taylor series and the multi-dimensional bubble functions are derived as the tensor products of one-dimensional functions for bilinear elements. This method has been used to model the flow in porous media and the convection–diffusion equation [14,15]. In this paper, we have extended the semi-discrete method to develop a practical scheme for the solution of diffusion-reaction equation for both dissipation and production cases. Residual free bubble (RFB) method is used in conjunction with bilinear elements and static condensation method (STC) is used for both bilinear and triangular elements to incorporate the elemental bubble functions and Lagrangian shape functions in a Galerkin finite element scheme. The elemental bubble function for bilinear elements is derived by the RFB method and an equivalent elemental bubble function is used for the triangular elements. Details of the procedure for deriving bubble enriched Lagrangian shape functions which can generate an appropriate scheme for the solution of the diffusion-reaction equation, together with their implementation are described. The performance of the developed method is evaluated through the comparison of the model results with the analytical solution of a bench mark problem. 2. Governing equation The steady state diffusion-reaction equation in domain X  Rd can be written as follows: kr  rT  sT ¼ f ;

ð1Þ

where k is the diffusion (conduction) coefficient and s is a source/sink term (s < 0 represents production and s > 0 stands for dissipation), T is the field unknown and f is a given source term. Using the following dimensionless forms: 8 T > : x ¼ x ; h where T1 is a reference value of the field variable, h is a characteristic length (e.g. width of the domain) and x represents position vector in the selected coordinate system. After substitution from Eq. (2) the governing equation is written in a dimensionless form as r  rT   Da T  ¼ f 

ð3Þ

in which f* is the dimensionless source term and Da is the Damkohler number, respectively, defined as 8 sh > > < Da ¼ ; k ð4Þ h > > :f ¼ f: kT 1 In a two-dimensional system (x, y) Eq. (3) can be written as  2   oT o2 T  þ 2  Da T  ¼ f  : ð5Þ ox2 oy Corresponding dimensionless boundary conditions for the rectangular domain are: (a) Dissipation: T  ¼ 0 for y  ¼ 0; 0 6 x 6 1 and x ¼ 0; 0 6 y  < 1; T  ¼ 1 for x ¼ 1; 0 6 y  < 1 and y  ¼ 1; 0 6 x 6 1: ð6Þ (b) Production: T ¼ 0 

for x ¼ 0;

0 6 y  6 1;



T ¼1 for x ¼ 1; 0 6 y  6 1;   oT oT ¼  ¼ 0 for y  ¼ 0 and y  ¼ 1 at 0 < x < 1:  ox oy ð7Þ 3. Residual free bubble functions The derivation of the bubble functions is based on the analytical solution of the model differential equation within each element using homogeneous boundary conditions. More specifically, to derive the appropriate bubble function for the present diffusion-reaction equation we follow the method described by Franca and Russo [12] and Brezzi et al. [16]. Let us consider a boundary value problem defined in X  R2 as  LT ¼ f in X; ð8Þ T ¼0 on C; where L is a linear differential operator and f is a given source function defined on X. The standard Galerkin method is formulated in a subspace Vh  V, where V is the space of functions for which a solution of the continuous problem is sought. The Galerkin method aims to find uh 2 Vh such that aðT h ; vÞ ¼ ðLT h ; vÞ ¼ ðf ; vÞ; ð9Þ where a(Æ, Æ) is a bilinear form and (Æ, Æ) represents the scalar product of its arguments. In a two-scale method, the unknowns are divided into two parts T h ¼ T 1 þ T b;

ð10Þ

where Tb is the fine scale and T1 represents a standard finite element approximation polynomial (interpolation function). The fact that bubble functions disappear on element boundaries makes it possible to remove the equations that

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

correspond to these functions from the set of elemental equations using the static condensation procedure [17]. In the static condensation procedure we set v = vb in Eq. (9) to obtain aðT 1 þ T b ; vb ÞXe ¼ ðf ; vb ÞXe :

ð11Þ

The subscript Xe indicates that formulation is restricted to individual elemental domains. In the residual free method it is assumed that the bubble functions satisfy the differential equation in each element, i.e. LT b ¼ ðLT 1  f Þ

in Xe :

ð12Þ

Therefore, in the variational form of Eq. (12), which is similar to Eq. (11), static condensation is satisfied automatically [16]. Eq. (12) can be solved by considering that T b ¼ T 0b þ T fb  T 0b and T fb are, respectively, solutions of the following equations [16]: ( LT 0b ¼ LT 1 in Xe ; ð13Þ T 0b ¼ 0 on Ce ; and ( LT fb ¼ f T fb ¼ 0

in Xe ; on Ce :

ð14Þ

Assuming that / is a bubble shape function and w is a polynomial shape function, then Eqs. (13) and (14) can be rewritten as:  L/i ¼ Lwi in Xe ; ð15Þ /i ¼ 0 on Ce ; where wi and /i are functions associated with node i. Xe is the element domain and Ce is the element boundary. Hence ( L/f ¼ f in Xe ; ð16Þ /f ¼ 0 on Ce and Th ¼ T1 þ Tb ¼

n X

T i ðwi þ /i Þ þ /f ;

ð17Þ

i¼1

where n is the number of nodes per element. To solve Eqs. (15) and (16) it is assumed that [13]: N i ¼ wi þ /i :

ð18Þ

For the diffusion-reaction equation with respect to Eq. (3), operator L is defined as L ¼ D  Da

ð19Þ

Substituting Eq. (18) into Eq. (15), for a linear element on each node we have 8 > d2 N 1 > >  Da N 1 ¼ 0 for x 2 ½0  l; < dx2  ð20Þ N 1 ð0Þ ¼ 1; > > > ¼ w ) N : 1 1 N 1 ðlÞ ¼ 0;

8 > d2 N 2 > >  Da N 2 ¼ 0 for x 2 ½0  l; < dx2  N 2 ð0Þ ¼ 0; > > > N ¼ w ) : 2 2 N 2 ðlÞ ¼ 1;

1097

ð21Þ

where l is the characteristic element length and wi is a linear shape function. The solution of the above equations gives bubble shape functions expressed in a local elemental coordinate system as: 8 pffiffiffiffiffiffi sinh Da ðl  xÞ > > > pffiffiffiffiffiffi ; N ¼ 1 < sinh Da l ð22Þ pffiffiffiffiffiffi > > ax > N 2 ¼ sinh pD ffiffiffiffiffiffi : sinh Da l by solving Eq. (16) /b, elemental bubble function [9], can be derived as: f / Da b /b ¼ /1 þ /2 : /f ¼

ð23Þ

3.1. Polynomial bubble functions-coefficients derived directly by the RFB method Hyperbolic functions (2) can only be directly used if the integrals in the elemental equations are evaluated manually. Using Taylor expansions and truncating after a selected term these functions are expressed as polynomials sinhðxÞ ¼ x þ

x3 x5 þ þ  3! 5!

For example, truncating after the second term third order polynomial bubble functions are derived as   ðl  x2 Þ ðl  xÞ 1 þ Da l  x xðl  xÞð2l  xÞ 6    ;   N1 ¼ ¼ Da 6 l l l 1 þ h2 þ l2 Da 6 



Da x 1 þ x2 x xðl  xÞðl þ xÞ 6 ¼    : N2 ¼  Da 2 6 l 2 l l 1þ l þl Da 6

ð24Þ

ð25Þ

In Eqs. (24) and (25) the second parts represent third order bubble functions: 8 xðl  xÞð2l  xÞ > >  ; /1 ¼  > > 6 > 2 > > l þl < Da ð26Þ > xðl  xÞðl þ xÞ > >   ¼ : / > 1 > 6 > > : þ l2 l Da Using a local coordinate system of n(1, +1) the bubble functions are written as:

1098

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

8 ð1  n2 Þð3  nÞ > >  ¼ bð1  n2 Þð3  nÞ; > /1 ¼  > > 6 > > 8 1þ < D a l2 2 > ð1  n Þð3 þ nÞ > >  ¼ bð1  n2 Þð3 þ nÞ; /2 ¼  > > 6 > > : 8 1þ D a l2 where n ¼ 1  2xl and b ¼ 

1

8 1þ

6 Da l2

For ninth order bubble the coefficients are: ð27Þ

.

Using a similar procedure fifth order bubble enriched bilinear element can also be derived as: 8 2 / ¼ B½b2 ð1  n2 Þ þ b3 ð1  n2 Þð1  nÞ þ b4 ð1  n2 Þ > > < 1 2 þb5 ð1  n2 Þ ð1  nÞ; 2 2 2 2 > > /2 ¼ B½b2 ð1  n Þ þ b3 ð1  n Þð1 þ nÞ þ b4 ð1  n Þ : 2 2 þb5 ð1  n Þ ð1 þ nÞ ð28Þ in which   1 Da l D2a l3 l2 B¼   ;  ; b2 ¼  3! 5! 4 Da l2 D2a l4 l 1þ þ 3! 5!   Da 3Da l2 l3 D 2 l5 D2 l5 b3 ¼   ; b4 ¼ a ; b5 ¼ a : 3! 5! 8 5! 8 5! 32 For seventh order bubble the coefficients are:   1 Da l D2a l3 D3a l5 l2 B¼  ¼  ; b þ ; þ  2 5! 7! 4 3! Da l2 D2a l4 D3a l6 l 1þ þ þ 7!  3! 2 25!  2  Da 3Da l 5D3a l4 l3 2Da l 4D3a l3 l4 þ þ ; b4 ¼ ; þ b3 ¼  5! 7! 8 7! 16 5!   23! Da 6D3a l3 l5 3D3a l l6 b5 ¼ þ ; b6 ¼  ; 7! 64 5!  7! 32  D 3 l7 b7 ¼  a : 7! 128

1 B¼  ; Da l2 D2a l4 D3a l6 D4a l8 l 1þ þ þ þ 3! 5! 7! 9!  2 3 3 5 4 7 2 Da l Da l Dl Dl l þ þ a þ a ; b2 ¼  3! 5! 7! 9! 4   Da 3D2a l2 5D3a l4 7D4a l6 l3 þ þ þ ; b3 ¼  3! 5! 7! 9! 8  2  2Da l 4D3a l3 6D4a l5 l4 þ þ ; b4 ¼ 5! 7! 9! 16  2  Da 6D3a l3 15D4a l4 l5 þ þ ; b5 ¼ 5! 7! 9! 32  3  3Da l 10D4a l3 l6 þ ; b6 ¼  7! 9! 64  4  8  3  Da 10D4a l2 l7 4Da l l þ ; b8 ¼ ; b7 ¼  9! 256 7! 9! 128  4 9 Da l : b9 ¼ 9! 512 In which bi is coefficients of ith order term in the polynomial bubble function (see Eq. (28)). 3.2. Polynomial bubble functions-coefficients derived by the static condensation (STC) method For the bilinear elements the static condensation method can be carried out in one dimension for the linear element and then, the bilinear bubble enriched element can be obtained by the tensor product of one-dimensional shape functions. For the triangular elements the STC method must be directly applied in two dimensions to incorporate the bubble functions with the Lagrangian shape functions. Furthermore, it must be noted that solu-

T*=1

1

1

T*=0 T*=1

y

T*=0 T*=1

y

x x

0

1 T*=0

Fig. 1. Domain and the boundary conditions for the dissipation case (see Eq. (6)).

0

1

Fig. 2. Domain and the boundary conditions for the production case (see Eq. (7)).

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

1099

procedure. However, it is possible to incorporate any given bubble function with Lagrangian shape functions by means of the static condensation procedure. With respect to the derived polynomial of the residual free method and with respect to Eq. (23), the elemental polynomial bubble functions can be derived as:

1 0.9 0.8

Second order elemental bubble function: /b ¼ ð1  n2 Þ:

0.5

ð29Þ

Forth order elemental bubble function: 2

/b ¼ ð1  n2 Þ þ ð1  n2 Þ :

ð30Þ

Sixth order elemental bubble function:

y

2

3

/b ¼ ð1  n2 Þ þ ð1  n2 Þ þ ð1  n2 Þ :

ð31Þ

According to above polynomials it can be concluded that a nth order elemental bubble function may be written as: 0

1

0.5

x

2

3

n

/b ¼ ð1  n2 Þ þ ð1  n2 Þ þ ð1  n2 Þ þ    þ ð1  n2 Þ n X q ¼ ð1  n2 Þ : ð32Þ

Fig. 3. 10 · 10 equal density mesh with rectangular elements, BLNR1.

q¼1

To incorporate the bubble functions with ordinary shape functions we have

3.2.1. Bilinear elements In the RFB method, condensation takes place automatically for the derived bubble functions. Therefore the bubble coefficients are calculated as a part of the condensation

where wi in this work is the Lagrangian linear shape function and /b is the polynomial bubble function. Using the static condensation procedure the bubble enriched onedimensional shape functions can be generally derived as:

uh ¼ w1 u1 þ w2 u2 þ /b ub ;

Da=200, y=0.5

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

Da=200, y=0.5

1 0.9 0.8

Exact Galerkin 3rd, RFB 5th, RFB

T*

T*

tion of a PDE in a triangle is not a trivial matter and hence, for the triangular elements it is not possible to derive the bubble functions directly from the RFB method.

Exact Galerkin 2nd, STC 4th,STC

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

X* Da=200, y=0.8

1 0.9

1

0.8

1

Exact Galerkin 2nd, STC 3rd, RFB

0.8 0.7 0.6

T*

T*

0.8

Da=200, y=0.9

1 0.9

Exact Galerkin 2nd, STC 3rd, RFB

0.8 0.7 0.6

0.6

X*

0.5 0.4 0.3 0.2 0.1

0.5 0.4 0.3 0.2 0.1 0

0

0

0.2

0.4

0.6

X*

0.8

1

0

0.2

0.4

0.6

X*

Fig. 4. Results for the Da = 200 at different cross sections, BLNR1 mesh scheme, dissipation.

1100

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

N i ¼ wi þ b/b in which b is the bubble coefficient and is derived during the implementation of the static condensation method (see Eq. (11)). The bubble coefficient derived by the static condensation method for different bubble orders are as follows: Second order bubble (see Eq. (29)): b¼

1 16 1:60 þ D a l2

:

3.2.2. Triangular elements The bubble function in Eq. (32) is composed of a second order bubble, /b = (1  n2), which is the minimum possible order of a bubble function for a linear element. This bubble can also be derived by the product of the linear Lagrangian shape functions. If we do so for a triangular element the minimum order bubble is achieved as /b ¼ ngð1  n  gÞ:

Fourth order bubble (see Eq. (30)): b¼

1 : 31:238 3:09 þ Da l2 Da=400, y=0.5

1.2

Exact Galerkin 3rd, RFB 5th, RFB 7th, RFB

1

T*

0.8

Higher order bubbles with the same configuration as Eq. (32) for triangular element can be written as 2

0.6

q¼1

0.2

After application of the static condensation method the bubble coefficient can be derived as:

0 0

0.2

0.4

0.6

0.8

1

X*

Second order bubble function: 8 /b ¼ ngð1  n  gÞ; > > > < 1 ! b¼ : 2 1 1 > > þ þ 0:071429 > : Da l2x l2y

Da=400, y=0.5

1.2

Exact Galerkin 2nd, STC 4th, STC 6th, STC 8th, STC

1 0.8

T*

3

/b ¼ ngð1  n  gÞ þ ðngð1  n  gÞÞ þ ðngð1  n  gÞÞ þ  n X q ¼ ðngð1  n  gÞÞ : ð33Þ

0.4

0.2 0 0

0.2

0.4

X*

0.6

0.8

1

Da=400

1 0.9 0.8 0.7 0.6

ð34Þ

Fourth order bubble function: 8 > / ¼ ngð1  n  gÞ þ ðngð1  n  gÞÞ2 ; > > b > < 1 ! b¼ ; > 1 1 1 > > þ þ 0:124368 > : Da l2x l2y

0.6 0.4

3rd, RFB, y=0.5

ð35Þ

Da=600, y=0.5

1.2

Exact Galerkin 3rd, RFB 6th, STC

1

3rd, RFB, x=0.5

0.8

T*

T*

Sixth order bubble (see Eq. (31)): 1  b¼ 49:643 4:52 þ Da l2 and for the eighth order bubble we have: 1 : b¼ 71:554 5:92 þ Da l2

0.5 0.4 0.3 0.2 0.1 0

0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

X* or Y* Fig. 5. Results for the Da = 400, BLNR1 mesh scheme, dissipation.

1

0

0.2

0.4

0.6

0.8

X* Fig. 6. Results for the Da = 600, BLNR1 mesh scheme, dissipation.

1

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107 Da=1000, y=0.5

1

Exact Galerkin 3rd, RFB 6th, STC

1 0.8 0.6

T*

T*

0.6

Da=1000, y*=0.8

1.2

Exact Galerkin 3rd, RFB 6th, STC

0.8

0.4

1101

0.4

0.2

0.2

0 -0.2

0

0

0.2

0.4

0.6

0.8

1 -0.2

0

0.2

0.4

Da =1000, y*=0.9

1.2

1

3rd, RFB 6th, STC

0.8

1

0.8

1

3rd, RFB, y*=0.5 3rd, RFB, x*=0.5

0.8

T*

0.6

T*

0.8

Da=1000

1.2

Exact Galerkin

1

0.6

X*

X*

0.4 0.2

0.6 0.4

0 -0.2 0

0.2

0.4

0.6

0.8

1

0.2 0

-0.4

0

0.2

0.4

X*

0.6

x* or y*

Fig. 7. Results for the Da = 1000, BLNR1 mesh scheme, dissipation.

where lx and ly are characteristic element length in the x- and y-directions in a Cartesian coordinate system, respectively. 3.3. Two-dimensional bubble functions for bilinear elements

8 N1 > > > > N3 > > > : N4

¼ 14 ð1  nÞð1  gÞ þ bð1  n2 Þð1  g2 Þð3  nÞð3  gÞ; ¼ 14 ð1 þ nÞð1  gÞ þ bð1  n2 Þð1  g2 Þð3 þ nÞð3  gÞ; ¼ 14 ð1 þ nÞð1 þ gÞ þ bð1  n2 Þð1  g2 Þð3 þ nÞð3 þ gÞ; ¼ 14 ð1  nÞð1 þ gÞ þ bð1  n2 Þð1  g2 Þð3  nÞð3 þ gÞ; ð36Þ

Two-dimensional bubble functions can be derived using tensor products of one-dimensional functions. The derived bubble functions are then incorporated into normal interpolation functions of bilinear Lagrangian elements to obtain shape functions of a bubble enriched bilinear element. For example for a third order bubble we have:

where

1  b¼  6 8 1þ D a l2 0.9

Da=2000, y*=0.5

1.2

0.8

Exact Galerkin 3rd, RFB 6th, STC

1 0.8

I II

0.7 0.5

III

T*

0.6 0.4 0.2

y

0 -0.2 0 -0.4

0.2

0.4

0.6

0.8

1

X*

Fig. 8. Results for the Da = 2000, BLNR1 mesh scheme, dissipation.

0

1 x

Fig. 9. BLNR2 Mesh scheme with quadrilateral elements.

1102

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

in which l is a characteristic element length. For the triangular elements the bubble functions are added to the Lagrangian shape functions as: 8 > < N 1 ¼ ð1  n  gÞ þ bngð1  n  gÞ; N 2 ¼ n þ bngð1  n  gÞ; > : N 3 ¼ g þ bngð1  n  gÞ;

Z

2

0 Xn 1 Xn o2 j¼1 N j T j o2 j¼1 N j T j 4W i @ A þ ox2 oy 2 Xe 3 n X   Da W i N j T j  W i f 5 dx dy ¼ 0: j¼1

To implement the multiscale procedure using bubble functions, the weights (Wi) are similar to the ordinary Lagrangian shape functions while the shape functions (Ni) are

where b is the bubble coefficient defined by Eq. (34). Similarly the Lagrangian shape functions can be enriched by higher order bubble functions. The weighted residual statement of the governing equation (5) can be written as

Da=-200 Exact 8th, STC

Galerkin 9th, RFB

6th, STC 3RD, RFB

2 1.5 1

T*

0.5 Da=1000, line III

1.2

0

-0.5 0

0.2

0.4

Exact

1 0.8

Galerkin

-1

3rd, RFB

-1.5

6th, STC

X*

0.6

0.8

1

-2

T*

0.6

Fig. 11. Results for the Da = 200, BLNR1 mesh scheme, production.

0.4 0.2

Da=-400

0 -0.2

0

0.2

0.4

0.6

0.8

Exact 9th, RFB

1

Galerkin 3rd, RFB

8th, STC

1.5

X* 1 Da=1000, line II

1

0.5

0.8

T*

Exact Galerkin 3rd, RFB 6th, STC

T*

0 0

0.6

0.2

0.4

0.4

0.8

1

-1

0.2

-1.5

X*

0 0

0.2

0.4

-0.2

0.6

0.8

1

Fig. 12. Results for the Da = 400, BLNR1 mesh scheme, production.

X* Da=-600

Da=1000, line I

1

Exact 8th, STC

Exact

T*

0.6

-0.5

0.8

Galerkin

2

0.6

3rd, RFB 6th, STC

1.5

Galerkin 9th, RFB

1

0.4

0.5

T*

0.2

-0.2

0

-0.5 0

0 0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

-1 -1.5

-0.4

X* Fig. 10. Results for the Da = 1000 at different cross sections BLNR2 mesh scheme, dissipation.

-2

X* Fig. 13. Results for the Da = 600, BLNR1 mesh scheme, production.

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

eliminated during the assembly of elemental equations. This problem does not become apparent in the one-dimensional case as the boundary integrals are reduced to simple nodal flux terms. The variational formulation for the diffusion-reaction equation, after application of Green’s theorem is

Da=-1000 Exact 9th, RFB

2

Galerkin 8th, STC

1.5 1

T*

0.5

ðrT h ; rvÞ þ ðDa T h ; vÞ ¼ ðf ; vÞ:

0 -0.5 0

1103

0.2

0.4

0.6

0.8

1

ð37Þ

Substitution from Eq. (17) gives: ðrT 1 ; rvÞ þ ðrT b ; rvÞ þ ðDa T h ; vÞ ¼ ðf ; vÞ:

-1 -1.5

ð38Þ

If v is a linear test function (weight function) according to Green’s theorem [18] we have:

-2

X* Fig. 14. Results for the Da = 1000, BLNR1 mesh scheme, production.

Da=-2000 2

ð39Þ

where / is bubble function. Therefore Eq. (38) is reduced to: ðrT 1 ; rvÞ þ ðDa T h ; vÞ ¼ ðf ; vÞ:

Galerkin 8th, STC

Exact 9th, RFB

ðrv; r/ÞXe ¼ ðDv; /ÞXe þ ðrv; /ÞCe ¼ 0;

ð40Þ

1.5

1

1

T*

0.5 0 -0.5 0

0.2

0.4

0.6

0.8

1

-1 -1.5 -2

0.5

X* Fig. 15. Results for the Da = 2000, BLNR1 mesh scheme, production.

identical to the bubble enriched shape functions (for example similar to Eq. (36)).

0.3 y

3.4. Elimination of the boundary integrals

0

When the discretisation involves bubble functions the inter element boundary integrals are not automatically

1 x

Fig. 16. Mesh scheme, BLNR3, with quadrilateral elements.

Table 1 Comparison between exact and numerical solution at the sampling points for some extreme cases of production Sampling points

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Da = 1000, BLNR1 mesh scheme

Da = 1000, BLNR3 mesh scheme

Da = 2000, BLNR1 mesh scheme

Exact

9th RFB

8th STC

Sampling points

Exact

9th RFB y* = 0.5

Sampling points

Exact

9th RFB

0.00000 0.11101 0.22136 0.33002 0.43632 0.53997 0.64153 0.73908 0.83201 0.91974 1.00000

0.00000 0.12717 0.24189 0.36174 0.46439 0.55766 0.67951 0.76812 0.87191 0.94954 1.00000

0.00000 0.19160 0.37664 0.54881 0.70222 0.83163 0.93261 1.00170 1.03660 1.03600 1.00000

0.00 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 1.00

0.00000 1.40430 1.38690 1.36113 1.32705 1.28459 1.23387 1.17544 1.109663 1.03695 1.00000

0.00000 1.20430 1.19850 1.33730 1.42060 1.35270 1.36590 1.47820 1.45200 1.34950 1.00000

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.00000 1.37510 0.64185 1.07522 1.14467 0.53948 1.39687 0.11067 1.34536 0.73686 1.00000

0.00000 1.36170 0.63305 1.06740 1.12930 0.54236 1.38140 0.09986 1.33500 0.72050 1.00000

1104

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

As can be seen the bubble function does not affect the Laplacian term in the diffusion-reaction equation and therefore no boundary integral due to the bubble function exists.

Da=-1000, y*=0.3 Exact 8th, STC

8

Galerkin 9th, RFB

6 4 2

T*

4. Analytical solution To validate the numerical solutions, the analytical solution of the dimensionless diffusion-reaction equation is presented. The solution for the boundary conditions represented in Eq. (6) is based on the method of separation of variables and can be written as: n 1 pffiffiffi 2X 1  ð1Þ T  ðx ; y  Þ ¼ bn pffiffiffi fsinh cy  sinðnpx Þ p n¼1 n sinh c pffiffiffi þ sinh cx sinðnpy  Þg;

0 -2 0

0.2

0.4

0.6

0.8

1

-4 -6 -8

X* Da=-1000, y*=0.5

ð41Þ

Galerkin 9th, RFB

Exact 8th, STC

10 8 6 4

where pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ n2 p 2 þ D a : T*

2 0 -2 0 -4 -6

For the boundary conditions represented in Eq. (7), a one-dimensional solution can be used

0.2

0.4

0.6

0.8

1

-8 -10

X*

Da=-400, y*=0.3 Galerkin 9th, RFB

Exact 8th, STC

2

Fig. 18. Results for the Da = 1000 at different cross sections, BLNR3 mesh scheme, production.

1.5 1 1

T*

0.5

0.9

0 -0.5 0

0.2

0.4

0.6

0.8

1

0.8

-1 -1.5 -2

X*

0.5

Da=-400, y*=0.5 Galerkin 9th, RFB

Exact 8th, STC

2.5 2

y

T*

1.5 1 0.5 0 -0.5 0 -1 -1.5

0.2

0.4

0.6

0.8

1

1

0 x

-2 -2.5

Fig. 19. Mesh scheme TRI1 with triangular elements.

X* Fig. 17. Results for the Da = 400 at different cross sections, BLNR3 mesh scheme, production.

T  ðx Þ ¼

pffiffiffiffiffiffi sinð Da x Þ pffiffiffiffiffiffi : sinð Da Þ

ð42Þ

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

The main objective of the present work is the construction of a multiscale finite element scheme for the diffusionreaction problem. To validate the scheme its results are compared with the analytical solutions obtained for a bench mark problem. To carry out the simulations, an in house computer code in FORTRAN has been developed. The results are presented for bilinear and triangular elements with different mesh schemes for both dissipation and production cases. In all figures second, third, fourth, etc. represent the order of the bubble functions, RFB is symbol of residual free bubble and STC represents static condensation.

Da=1000, y*=0.5, TRI1

1

Exact 4th, STC

0.8

2nd, STC Galerkin

0.6 0.4 0.2

Da= 200, y*=0.5, TRI1

1 0.9 0.8 0.7 0.6

0

Exact 4th, STC

-0.2

2nd, STC Galerkin

0

0.2

0.4

Exact 4th, STC 2nd, STC Galerkin

0.8 0.6

0.4

X*

0.6

1

Da=1000, y*=0.9, TRI1

1

0.2

0.8

X*

0.5 0.4 0.3 0.2 0.1 0 0

0.6

0.8

1

0.4

T*

T*

Figs. 1 and 2 show the domain and its boundaries for the dissipation and production cases, respectively. Figs. 1–8 show the results for the dissipation case with bilinear elements for BLNR1 mesh scheme shown in Fig. 3 for Da up to 400, bubble functions with different orders. For higher Da values third and sixth order bubbles are compared, where stable-accurate results are generated.

T*

5. Results and discussion

1105

0.2

Da =200, y*=0.8, TRI1

1

0

Exact 4th, STC 2nd, STC Galerkin

0.9 0.8 0.7

-0.2

0

0.2

0.4

-0.4

T*

0.8

1

X*

0.6 0.5

0.6

Fig. 21. Results for the Da = 1000 at different cross sections, TRI1 mesh scheme, dissipation.

0.4 0.3 0.2

1

0.1 0

0.9

0

0.2

0.4

0.6

0.8

1

X* Da= 200, y*=0.9, TRI1

1 0.9

Exact 4th, STC 2nd, STC Galerkin

T*

0.8 0.7 0.6

y

0.5 0.4 0.3 0.2

0.4

0.1 0 0

0.2

0.4

X*

0.6

0.8

1

Fig. 20. Results for the Da = 200 at different cross sections, TRI1 mesh scheme, dissipation.

0

1 x

Fig. 22. Mesh scheme, TRI2, with triangular elements.

1106

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107

The numerical results show that the third order RFB method has the same results as the sixth order STC method. Fig. 10 shows the results using BLNR2 mesh scheme, represented in Fig. 9 at Da = 1000, as an extreme case. At the line I which is close to the boundary the third order shows a better performance than the sixth order while has some degree of instability. The results with the BLNR2 mesh scheme show that the mapping error has no deteriorating effect on the results and only close to the boundaries slight deviation from exact solution is observed at Da = 1000. Figs. 11–15 show the results with BLNR1 mesh scheme for production case where the boundary conditions of Fig. 2 is used. At Da = 200 Galerkin method can not match with the analytical solution while the ninth order

Da=1000, y*=0.4, TRI2

1

Exact 4th, STC

0.8

2nd, STC Galerkin

T*

0.6 0.4 0.2

bubble can cope with the analytical solution. At Da = 400 the third order bubble fails to predict correct solution while in the dissipation case gives stable-accurate result even for the Da = 2000 (see Fig. 8). In the dissipation case the multiscale behaviour happens only in the diffusion boundary layer while in the production case the multiscale behaviour is repeated over the entire domain. Figs 13–15 show that the ninth order bubble copes very well with the analytical solution even at the Da = 2000. Table 1 compares the numerical values at the sampling points. To test the performance of the bubble functions in the production case with non-rectangular elements the mesh scheme shown in Fig. 16, BLNR3, was used. As the results in Figs. 17 and 18 show, the eighth order bubble function which is incorporated to the bilinear element by the STC method fails to predict correct results at both Da = 400 and 1000 while the ninth order bubble of RFB method gives accurate results according to the exact solution. Figs. 20 and 21 represent the results with triangular elements with TRI1 mesh scheme shown in Fig. 19. At Da = 200 standard Galerkin and bubble functions give stable results while at Da = 1000 solution is unstable although, the bubble enriched elements give more accurate results according to the exact solution. If the results in Fig. 21 is compared with those in Fig. 23 it becomes clear

0 0

0.2

0.4

0.6

0.8

1

-0.2

k=-200, y*=0.4, TRI1

X*

Exact

Da =1000, y*=0.9, TRI2

1

Exact 4th, STC 2nd, STC Galerkin

0.8

T*

T*

0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

3 2.5 2 1.5 1 0.5 0 -0.5 0 -1 -1.5 -2 -2.5

4th, STC

0.2

0.4

0.6

Galerkin

0.8

1

X*

-0.2

X* Da=-200, y*=0.4, TRI2 Da=1000, y=0.4

1

Exact

Galerkin, TRI1

0.8

Galerkin, TRI2

0.4

T*

T*

0.6

0.2 0 0

0.2

0.4

0.6

0.8

1

-0.2

X* Fig. 23. Results for the Da = 1000 at different cross sections, TRI2 mesh scheme, dissipation.

3 2.5 2 1.5 1 0.5 0 -0.5 0 -1 -1.5 -2 -2.5

0.2

4th, STC

0.4

0.6

Galerkin

0.8

1

X* Fig. 24. Results for the Da = 200, TRI1 and TRI2 mesh schemes, production.

M. Parvazinia, V. Nassehi / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1095–1107 Da=-400, y*=0.4, TRI2 Exact

4th, STC

Galerkin

2.5 2 1.5 1

T*

0.5 0 -0.5 0

0.2

0.4

0.6

0.8

1

-1

1107

those derived by the STC method. The third order RFB bubble performs better than the sixth order STC bubble function in the dissipation case. In the production case the ninth order RFB bubble shows excellent results. For the triangular elements it seems that the mesh structure have a significant effect on the accuracy of the results. The fourth order bubble yields accurate results in the production case for a relatively wide range of Damkohler number. Using higher order bubble functions accurate solutions can be obtained for Da < 400.

-1.5

References

-2 -2.5

X* Fig. 25. Results for the Da = 400, TRI2 mesh scheme, production.

Da=-600, y*=0.4, TRI2 Exact

4th, STC

Galerkin

8 6 4

T*

2 0

-2 0

0.2

0.4

0.6

0.8

1

-4 -6 -8

X* Fig. 26. Results for the Da = 600, TRI2 mesh scheme, production.

that TRI2 mesh scheme, shown in Fig. 22, improves the results. According to the last graph in Fig. 23 which compares the Galerkin solutions using TRI1 and TRI2 mesh schemes at the same cross sections, it can be concluded that the improvement is likely to be the result of mapping error reduction in the TRI2. In the production case, as Fig. 24 shows, using TRI1 mesh scheme the accurate results can not be achieved by the fourth order bubble. The fourth order bubble can produce acceptable results for Da = 200 and 400 with TRI2 mesh scheme while it fails for Da = 600. The accurate results can be obtained using a higher order bubble function which is out of the scope of the present paper (see Figs. 25 and 26). 6. Conclusions A multiscale finite element model based on bubble functions has been developed which yields stable-accurate solutions for the diffusion-reaction equation using the Galerkin scheme for a range of Damkohler number. Both bilinear and triangular elements are used with different bubble functions for production and dissipation cases of diffusion reaction equation. The results show that the bubble functions derived by the RFB method yields better results than

[1] C. Johnson, U. Navert, J. Pitkarenta, Finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 285–312. [2] M. Parvazinia, V. Nassehi, R.J. Wakeman, M.H.R. Ghoreishy, Finite element modelling of flow through a porous medium between two parallel plates using the Brinkman equation, Transport Porous Med. 63 (2006) 71–90. [3] V. Gravimier, W.A. Well, E. Ramm, A three-level finite element method for the instationary incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1323–1366. [4] T.J.R. Hughes, Multi scale phenomena, Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 381–401. [5] T.J.R. Hughes, A space time formulation for multiscale phenomena, Comput. Methods Appl. Mech. Engrg. 74 (1996) 217–229. [6] F. Brezzi, M. Bristeau, L.P. Franca, M. Mallet, G. Roge, A relationship between stabilized finite element methods and the Galerkin method with bubble functions, Comput. Methods Appl. Mech. Engrg. 96 (1992) 117–129. [7] C. Baiocchi, F. Brezzi, L.P. Franca, Virtual bubbles and Galerkinleast-squares type methods, Comput. Methods Appl. Mech. Engrg. 105 (1993) 125–141. [8] L.P. Franca, T.J.R. Hughes, R. Stenberg, Stabilized finite element methods for the Stokes problem, in: R.A. Nicolaides, M.D. Gunzberger (Eds.), Incompressible Fluid Dynamics – Trends and Advances, Cambridge University Press, Cambridge, 1993. [9] F. Brezzi, L.P. Franca, T.J.R. Hughes, A. Russo, Comput. Methods Appl. Mech. Engrg. 145 (1997) 339–392. [10] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residual free bubbles for Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 40 (1997) 4003–4009. [11] L.P. Franca, A. Russo, Deriving upwinding, mass lumping and selective reduced integration by residual free bubbles, Appl. Math. Lett. 9 (1996) 83–88. [12] L.P. Franca, A. Russo, Unlocking with residual-free bubbles, Comput. Methods Appl. Mech. Engrg. 142 (1997) 361–364. [13] L.P. Franca, A. Russo, Mass lumping emanating from residual free bubbles, Comput. Methods Appl. Mech. Engrg. 142 (1997) 353–360. [14] M. Parvazinia, V. Nassehi, R.J. Wakeman, Multi-scale finite element modelling of laminar steady flow through highly permeable porous media, Chem. Engrg. Sci. 61 (2006) 586–596. [15] M. Parvazinia, V. Nassehi, R.J. Wakeman, Multi-scale finite element modelling using bubble function method for a convection–diffusion problem, Chem. Engrg. Sci. 61 (2006) 2742–2751. [16] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residualfree bubbles for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 166 (1998) 25–33. [17] J. Donea, A. Huerta, Finite Element Methods for Flow Problems, Wiley, Chichester, 2003. [18] L.P. Franca, C. Farhat, Bubble functions prompt unusual stabilized finite element methods, Comput. Methods Appl. Mech. Engrg. 123 (1995) 299–308.