An elemental scale adjustment method for the finite element and finite difference solutions of diffusion-reaction and convection–diffusion equations

Finite Elements in Analysis and Design 58 (2012) 1–9

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

An elemental scale adjustment method for the finite element and finite difference solutions of diffusion-reaction and convection–diffusion equations M. Parvazinia n Iran Polymer and Petrochemical Institute, Tehran, PO Box 14965-115, Iran

a r t i c l e i n f o

abstract

Article history: Received 14 December 2011 Received in revised form 4 April 2012 Accepted 5 April 2012 Available online 27 April 2012

A method for the stable-accurate solution of convection–diffusion and diffusion-reaction equations is proposed to produce a solution similar to highly refined meshes for any level of discretisation. The method is applied on finite element and finite difference methods. The idea is based on the analytical or numerical solutions of the governing equation in a test domain and then determining the level of adjustment in an element with length l by solving the global system of equations in the test domain. The adjustment can be done either on the coefficients of the equation or enrichment of normal shape functions in the Galerkin finite element scheme. The numerical experiments are performed in one and two dimensional cases. Different mesh schemes and boundary conditions are used. To validate the approach, the numerical results obtained for a benchmark problem are compared with the analytical ¨ solution in a wide range of Damkohler and Peclet numbers. & 2012 Elsevier B.V. All rights reserved.

Keywords: Finite element Finite difference Enrichment Diffusion-reaction equation Convection–diffusion equation

1. Introduction In the case of sharp gradients in which the field variables vary rapidly within thin layers, inaccurate and unstable solution is observed in the classical numerical solution of such singularly perturbed problems. If the variations take place in a length scale smaller than the discretisation level the unstable-oscillatory behaviour is observed. It means that the diffusion is limited into a thin boundary layer or interlayer and the process is reaction or convection dominant. Using sufficiently fine computational grids the stable-accurate solution can be achieved but, it is not computationally cost effective particularly for practical problems [1]. Many multiscale finite element schemes have been developed for different types of problems [2]. Using bubble functions to enrich the usual polynomial shape functions can make the solution stable and accurate for the above mentioned situations. Bubble functions are nonzero within the element and vanish on the element boundaries [3–5]. Therefore any desired bubble function can be included in the usual polynomial shape functions by using static condensation (STC) method. Using the STC method the resulting shape function is divided into two parts as Ni ¼ ci þbfi in which ci is usual polynomial shape function associated with node i, fi is the bubble function associated with node i and b is the bubble coefficient. The essential question is how to calculate b to satisfy the perfect solution at least in one dimensional problem?

n

Tel.: þ98 21 48662428. E-mail addresses: [email protected], [email protected]

0168-874X/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.finel.2012.04.001

During the use of STC method b is calculated but at sharp gradients (highly reaction or convection dominant cases) the solution becomes over-diffusive. By increasing the order of the bubble function the numerical approximation approaches the exact solution at nodal points [6]. Using the residual free bubble function (RFB) method makes it possible to have perfect solution at nodal points for any level of discretisation. In this method the analytical solution of the governing differential equation within each element subject to homogeneous boundary conditions is derived and the resulting shape function is an analytical function [7–11]. Derivation of analytical solution in multidimensional problems and implementation in the finite element scheme is not a simple task and needs further considerations to apply in the practical problems [12]. The variational multiscale (VMS) method, proposed by Hughes [13], is based on the dividing of field unknown (T) into two parts. The field variable is written as T¼T1 þ Tb.T1 is approximated by the standard polynomial finite element discretisations. Although the VMS provides a firm theoretical framework to improve our understanding of the stabilised methods and derivation of subgrid scale models but it has to be noted that the derivation of the fine scale Tb based on the element Green’s function is not a trivial matter. In this paper to determine the level of adjustment a global solution in a test domain with just one non-boundary node, is considered. At the 1st step the exact solution in the test domain can be derived using an analytical or a numerical solution. In the next step, the level of adjustment in the element of the length l is derived via the global system of equations by having the exact solution from the 1st step. The adjustment can be done by

2

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

Nomenclature b Function coefficient c Adjustment parameter Da Damkohler number f A given source term h A characteristic length (e.g., width of the domain) K Diffusion coefficient l1 and l2 Element length Pe Peclet number q1 and q2 Boundary integrals S Source/sink term T Field variable T1 Numerically resolved part of the solution Tb Numerically unresolved part of the solution T0 and T1Reference values for independent variable in Eqs. (7) and (9) T1, T2, T3 Nodal values of unknown in the global system of equations in the test domain

modification of the equation parameters or enrichment of the approximation functions. The method is applied in finite element and finite difference methods. The study is done in a wide range of Damkohler and Peclet numbers and different mesh schemes in one and two dimensional problems.

T1 v v w x

A reference value of the field variable Velocity vector Test function in the Eqs. (26) and (29)–(31) Weight function Position vector in the selected coordinate system

Superscript n

Dimensionless variable

Greek letters

a

Value of unknown inside the test domain Usual polynomial shape function Enriching shape function Variable in local coordinate system

c f x

where T is field unknown v is the velocity vector, k is diffusivity and f is a source term. r denotes the spatial gradient operator. Using the below mentioned dimensionless parameter [12]: (

T ¼ T 0 þ T n ðT 1 T 0 Þ xn ¼

ð7Þ

x h

2. Governing equation The steady state diffusion-reaction (DR) equation in domain O CRd can be written as follows kr:rTsT ¼ f

ð1Þ

In which k is the diffusion (conduction) coefficient and s is a source/sink term (so0 represents production and s40 stands for dissipation), T is the field unknown and f is a given source term. We prefer dimensionless form of the governing equations to work with dimensionless groups with clear physical meanings. Using the following dimensionless forms 8 < T n ¼ TT 1 ð2Þ : xn ¼ hx 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 n Da T n ¼ f n

ð3Þ

in which fn is the dimensionless source term and Da is the Damkohler number, respectively, defined as 8 < Da ¼ n

:f ¼

sh k h kT 1

where T0 and T1 are reference values for independent variable (e.g., temperature), and h is a characteristic length (e.g., width of the domain) dimensionless governing equation becomes: Pe rT n r:rT n ¼ f

n

ð8Þ

In which Pe is the Peclet number defined by: Pe ¼

vh k 2

n

f ¼f

h kðT 1 T 0 Þ

ð9Þ

Eq. (8) can be written in two dimensions as: !   @T n @T n @2 T n @2 T n n P ex n þ P ey n  þ ¼f @x @y @xn2 @yn2

ð10Þ

pex and pey correspond to x and y components of the velocity vector. It is assumed that Pe is equal in both directions. Corresponding dimensionless boundary conditions for the rectangular domain are: a) BC1:

f

ð4Þ

Tn ¼ 0 n

T ¼1 In a two-dimensional system (x,y) Eq. (3) can be written as ! @2 T n @2 T n n þ ð5Þ Da T n ¼ f @xn2 @yn2 The steady-state homogeneous convection-diffusion (CD) equation can be described by: v:rTkr:rT ¼ f

ð6Þ

for yn ¼ 0, 0 r xn r1 and xn ¼ 0, 0 r yn o 1 for xn ¼ 1, 0 r yn o 1 and yn ¼ 1, 0 r xn r 1

ð11Þ

b) BC2: for xn ¼ 0, 0 r yn r 1

Tn ¼ 0

for xn ¼ 1, 0 ryn r 1

n

T ¼1 @T n @xn

¼

@T n @yn

¼0

for yn ¼ 0 and yn ¼ 1 at 0 o xn o1

ð12Þ

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

3. The proposed method The fundamental idea of the proposed method is to construct global system of equations in a test domain in which just one node lays within the domain and other nodes are boundary nodes as shown in Fig. 1. Since the field variable is just unknown at one node (inner node) therefore the global system of equations reduces to one equation. By having the solution in the computational test domain via an analytical or numerical solution, the exact value of the unknown is determined at the inner node within the test domain. Therefore what remains unknown in the global set of equations in the test domain (which is now just one equation) is the level of adjustment that appears as a parameter. In this work the level of adjustment is calculated by adjusting the model parameter in each element in finite element or finite difference methods or by enriching the shape function in the finite element method. For practical means a one dimensional test domain can be used. In one dimension, the domain consists of just two elements as shown in Fig. 2. It is assumed that at the middle node the exact value is known by an analytical or numerical solution. To extend the results to multidimensional problems instead of the element length l, we have a characteristic length of the element. Other parameters in calculating the level of enrichment are the equation parameters like Peclet and Damkohler numbers which remain unchanged. In the Galerkin finite element method in the case of shape function enrichment, extension to multidimensional cases are done similar to normal lagrangian shape functions by tensor product of one dimensional functions. 3.1. Adjustment of the parameters of equations As mentioned earlier in Section 3 it is possible to modify the parameters of equations such that the stable-accurate solution can be achieved. The upwind methods [14] proposed using the artificial coefficients for stable solution in the expense of having over-diffusive solution particularly in the case of sharp gradients but, in the proposed method the basis of calculating the coefficients is quite different and the solution at the nodal points is perfect.

3

In this section we apply the proposed method for both diffusion-reaction and convection–diffusion equations using finite difference and finite element methods. 3.1.1. Application to Galerkin finite element scheme 3.1.1.1. DR Equation. As mentioned earlier for practical means a one dimensional test domain can be used (Fig. 2). By constructing the global stiffness matrix the level of adjustment is determined by calculating the parameter c in the governing equation in the test domain as (see Fig. 2)

r:rT n cDa T n ¼ 0 in Ot The boundary conditions are: n

T ¼0

at x ¼ 0

Tn ¼ 1

at x ¼ l1 þ l2

In which Ot is the test domain and c is adjustment parameter. The global system of equations is constructed as [12]: 2

1 l1

þ cD3a l1

6 6  1 þ cDa l1 6 l1 6 4 0

 l11 þ cD6a l1 1 l1

2 3 2 3 q1 3 T1 6 7 6 7 7 6 6 76 7 6 7 7 6 7  l12 þ cD6a l2 7 T2 7 76 7¼60 7 56 6 7 6 7 1 þ cD3a l2 4 5 4 5 l2 q2 T3 0

þ cD3a l1 þ l12 þ cD3a l2  l12 þ cD6a l2

ð13Þ where, l1 and l2 are elements length in the test domain, T1 ¼0,T3 ¼1(boundary conditions for the test domain), T2 ¼ a and c is the adjustment parameter multiplied by Da. By solving the middle row equation c is calculated as: c¼

ðð1aÞ=l2 Þða=l1 Þ Da ððaðl1 þl2 Þ=3Þ þ ðl2 =6ÞÞ

ð14Þ

In this work a is calculated analytically: pffiffiffiffiffiffi sinð Da l1 Þ pffiffiffiffiffiffi a¼ sinð Da hÞ

ð15Þ

3.1.1.2. CD Equation. Similar to DR equation, the level of adjustment is determined by calculating the parameter c in the governing equation in the test domain. cP e rT n r:rT n ¼ 0 in Ot

The boundary conditions are: n

T ¼0

at x ¼ 0

Tn ¼ 1

at x ¼ l1 þ l2

Fig. 1. Test domain in two dimensions.

The global system of equations for the CD equation in the test domain can be written as [12]: 2

1 l1

 12 cP e

6 6  1  1 cP e 6 l1 2 4 0 Fig. 2. Test domain in one dimension.

 l11 þ 12 cPe 1 þ l12 l1  l12  12 cP e

 l12  l12

2 3 2 3 q1 3 T1 6 7 6 7 0 7 6 7 76 6 7 6 7 þ 12 cP e 7 0 7 T2 7 ¼6 76 6 6 7 56 7 7 6 7 1 þ 2 cP e 4 5 4 5 q2 T3

ð16Þ

4

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

T1 ¼0,T3 ¼1 (boundary conditions for the test domain), T2 ¼ a and c is adjustment parameter multiplied by Pe. By solving the middle row equation c is calculated as

where a is the analytical solution using the BC on the test domain as mentioned above.

a¼ 1=l2 ð1=l1 þ1=l2 Þa c¼ 0:5Pe

1 þexpðl1 P e Þ expðP e hÞ1

ð24Þ

ð17Þ 3.2. Adjustment of approximation function

In this work a is calculated analytically: 1 þexpðl1 P e Þ expðPe hÞ1

a¼

3.1.2. Application to finite difference scheme 3.1.2.1. DR Equation. The discrete form of the governing equation in the finite difference method in one dimension can be written as T k þ 1 2T k þT k1

cDa T k ¼ 0 in Ot 2 l The boundary conditions are T k1 ¼ 0 at node 1 ðk1Þ T K þ 1 ¼ 1 at node 3 ðk þ1Þ

ð18Þ

where, c is the adjustment parameter and the test domain is shown in Fig. 2. By imposing boundary conditions and replacing Tk þ 1 ¼1,Tk  1 ¼0 the known value for Tk ¼ a the coefficient c is calculated as c¼

12a

ð19Þ

2

Da l a

where a is calculated analytically as:

2

The FD discretization is     d T k þ 1 T k T T cPe k k1 ¼ 0 dx l2 l1 T k þ 1 T k 2

l2



  T k T k1 T T cP e k k1 ¼ 0 l1 l2 l1

r:rT n Da T n ¼ f n in O Tn ¼ a

where O CR2 is a bounded domain with boundary G and fn is a given source function. The variational formulation of the above problem is finding T nh , the approximation of Tn in Galerkin method, such that

ð20Þ

ð21Þ

ð22Þ

T k þ 1 ¼ 1,T k1 ¼ 0,T k ¼ a

ð26Þ

(.,.) denotes inner product and v is the test function. Any function which is non-zero within the element can be regarded as a bubble function. We use this definition for the representation of the enriching function as ( f ¼ 1 in Oe ð27Þ f ¼ 0 on Ge

T nh ¼

nod X

ðc þ bfÞT i

in Oe

ð28Þ

i¼1

In which c denotes the usual polynomial shape function, f is the enriching shape function represented in (27), b is the function coefficient which determines the level of enrichment and ‘‘nod’’ denotes the number of nodes in each element. In this method a test domain is used in which just one node lays within the domain and other nodes are boundary nodes as shown in Fig. 1. To determine the level of enrichment the governing equation is solved in the test domain. The above variational formulation is written in each element as ! ! nod nod X X  r ðc þ bfÞT i , rv þ Da ðc þbfÞT i ,v ¼ ðf ,vÞOk i¼1

i¼1

Ok

Ok

ð29Þ

Replacing in Eq. (22)   1a a a  cP ¼0 e 2 l l l 1 2 1 l2

Summation of above formula over all elements in the test domain gives 8 ! ! 9 en < nod nod en = X X X X ¼  r ðc þ bfÞT i , rv þ Da ðc þ bfÞT i ,v ðf ,vÞOi : ;

c is derived as 2

l1 ð1aÞ=l2 a=l2 Pe a

v A H10 ðOÞ

Therefore in each element T nh can be written as

Based on the proposed method we have

c¼

ð25Þ

on G

n

3.1.2.2. CD Equation. The CD equation in the test domain is written as d T dT ¼ 0 in Ot þ cP e dx dx2 The boundary conditions are T k1 ¼ 0 at x ¼ 0 T K þ 1 ¼ 1 at x ¼ l1 þ l2

3.2.1. DR Equation The model problem is represented as

ðrT nh , rvÞ þ ðDa T nh ,vÞ ¼ ðf ,vÞ

pffiffiffiffiffiffi sinð Da l1 Þ pffiffiffiffiffiffi a¼ sinð Da hÞ



In this section the proposed method is applied to enrich the approximation function. The overall procedure is similar to what presented in Section 3.1 but here instead of calculating the adjusted coefficients, the parameter c, the enrichment function is calculated. For this reason the Galerkin finite element method is used and the normal shape functions are enriched. The enrichment function is added to the normal shape functions apparently similar to the addition of bubble functions but, the method of calculating the level of enrichment is different. The enrichment function is vanished at the element boundaries and is continuous inside the element. In below the method is applied on the DR and CD equations.

i¼1

ð23Þ

i¼

Oi

i¼1

Oi

i¼1

ð30Þ

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

In which, en is the number of elements in the test domain. Pen P Pnod fðr nod i ¼ 1 cT i , rvÞOi ðDa i ¼ 1 cT i ,vÞOi þðf ,vÞOi g ð31Þ b¼ i¼1 Pnod Pnod ðDa i ¼ 1 fT i ,vÞðr i ¼ 1 cT i , rvÞ For practical means a one dimensional test domain can be used. In one dimension, the domain consists of just two elements as shown in Fig. 2. It is assumed that at the middle node the exact value is known by an analytical or numerical solution as explained in Section 3.1. By having the value of unknown inside the test domain what remains unknown in the global set of equations in the test domain is the level of enrichment, b. The homogeneous governing equation and boundary conditions in the test domain are (see Fig. 2):

r:rT n Da T n ¼ 0 in Ot Tn ¼ 0

at x ¼ 0

Tn ¼ 1

at x ¼ l1 þ l2

ð32Þ

The elemental stiffness matrix for the above equation is   3 2R  R l dc2 dw1 l dc1 dw1 þ Da w1 N1 dx þDa w1 N2 dx 0 0 dx dx dx dx 6    7 ð33Þ 4 R l dc dw 5 R l dc2 dw2 1 2 þ Da w2 N1 dx þDa w2 N2 dx 0 0 dx dx dx dx In which ci represents the linear Lagrangian shape functions (

In which f is the enriching function and b is its coefficient. In the proposed method the enrichment function is ( f ¼ 1 in Oe ð36Þ f ¼ 0 at Ge where, Oe indicates the element interior and Ge element boundary. To calculate b the global system of equations is written as AT ¼ B

ð37Þ

We can decompose A into two parts, ordinary and enriched part. A ¼ GþE

ð38Þ

In which G represents ordinary Galerkin part and E is the enriched part. Therefore we have ET ¼ BGT

ð39Þ

Since we have just two elements the global matrix is 3  3, Eq.(39) can be written on the 2nd row as:

j¼1

E2j T j ¼ Bj 

3 X

G2j T j

ð40Þ

j¼1

E21 T 1 þE22 T 2 þ E23 T 3 ¼ B2 ðG21 T 1 þ G22 T 2 þ G23 T 3 Þ

Eij ¼ beij

ð45Þ

Thus b¼

G22 a þ G23 ae22 þ e23

ð46Þ

To calculate b, the global system of equations has to be constructed. As mentioned above, the global matrix can be written in two parts as A¼G þbe therefore, we have 2 1 Dl 3 þ 3a 1  l11 þ D6a l1 0 l 6 1 7 D l D l 7 6 1 Dl 1 Dl 1 1 A ¼ 6  l1 þ 6a 1 l1 þ 3a 1 þ l2 þ 3a 2  l2 þ 6a 2 7 4 5 1 0  l12 þ D6a l2 þ D3a l2 l2 2 Da l1 3 Da l1 0 2 2 6 Da l1 Da Da l2 7 7 ð47Þ þb6 2 ðl1 þ l2 Þ 2 5 4 2 Da l2 Da l2 0 2 2 Solving the global set of equations in the test domain we have. ððl1 þl2 Þa=l1 l2 Þð1=l2 Þ þððl2 Da þ 2ðl1 þ l2 ÞDa aÞ=6Þ 0:5Da ðl2 þ ðl1 þ l2 ÞaÞ

ð48Þ

3.2.2. CD Equation The same procedure can be applied for the CD equation. The elemental stiffness matrix based on the local coordinate x( 1,þ1) is:    3 2R R þ 1 2 dc2 dw1 þ 1 2 dc1 dw1 1 2 dx dx þ P e w1 dN þ P e w1 dN 1 1 dx dx l dx dx l dx dx 6   7 4 R þ 1  dc dw 5 R þ 1 2 dc2 dw2 dN 1 dN2 2 1 2 þ P e w2 dx dx þ P e w2 dx dx 1 1 l dx dx l dx dx ð49Þ The relation between x( 1, þ1)and x (0, l) is x ¼(l/2)(1þ x). Since in this method the enriching function is a constant number, its derivative is zero. It is proved that the level of enrichment of a symmetric function and its derivative within the element is equal [15]. Z þ1 Z þ1 Z þ1 W i fi dx ¼ W i fxi dx ¼ fi dx ð50Þ 1

1

1

With respect to the above equation the level of enrichment for the convection term can be calculated. 2 3 1 2 3  12 Pe  l11 þ 12 P e 0 0 2pe 2pe l1 6 7 1 6 1 1 6 7 þ l12  l12 þ 12 P e 7 A ¼ 6  l1  2 P e 7 þ b4 2pe 4pe 2pe 5 l1 4 5 1 1 1 1 0 2pe 2pe 0   Pe þ Pe l2

2

l2

2

ð51Þ We have: T 1 ¼ 0,T 2 ¼ a,T 3 ¼ 1

ð52Þ

Therefore: ð41Þ

Nodes 1 and 3 are boundary nodes and the value of the field variable at node 2 is a known value i.e., T 1 ¼ 0,T 2 ¼ a,T 3 ¼ 1

ð44Þ

Similar to bubble functions the enrichment does not apply on the laplacian term [6].Eij can be written as

ð34Þ

The enriched shape functions consist of two parts (similar to bubble function), an ordinary Lagrangian shape function and enrichment function as ( N 1 ¼ c1 þ bf ð35Þ N 2 ¼ c2 þ bf

3 X

Since Ni ¼ ci þ bf

b¼

c1 ¼ w1 ¼ 1 xl c2 ¼ w2 ¼ xl

5

ð42Þ

b¼

ðða1Þ=l2 Þ þ ða=l1 Þ þ 0:5pe pe ð4a þ2Þ

ð53Þ

4. Results and discussion

Therefore Eq. (41) can be written as: E22 a þ E23 ¼ ðG22 a þG23 Þ

ð43Þ

In this work the proposed method is evaluated for DR and CD equations. To validate the method its results are initially compared

6

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

1

Da=1000

1

Exact

0.8

0.9

FD

0.6

0.8

MFD

0.4 0.2

T*

0 -0.2 -0.4

0.5

-0.6 -0.8 -1

0

0.2

0.4

0.6

0.8

1

X*

y

Fig. 6. Comparison of normal FD and Modified FD for DR equation in one dimension Da ¼1000.Non-equal mesh size. Dissipation case.

Da=-500

4

0

0.5

x

1

Exact

3

FD

MFD

2 Fig. 3. Mesh scheme BLNR1.

T*

1

1

0 -1 -2

0.8

-3 -4

0

0.2

0.4

0.6

0.8

1

X* Fig. 7. Comparison of normal FD and Modified FD for DR equation in one dimension Da ¼  500. Production case.

y

0

0.8

x

Da=-2500

5

1

Exact

4

FD

MFD

3

Fig. 4. Mesh scheme BLNR2.

2

Da=1000

1 0.9

-2 -3

FD

0.7

-4

MFD

0.6

-5

0.5

0

0.2

0.4

0.6

0.8

1

X*

0.4

Fig. 8. Comparison of normal FD and Modified FD for DR equation in one dimension Da ¼  2500. Production case.

0.3 0.2 0.1 0

0 -1

Exact

0.8

T*

T*

1

0

0.2

0.4

0.6

0.8

1

X* Fig. 5. Comparison of normal FD and Modified FD for DR equation in one dimension Da ¼1000. Dissipation case.

with the analytical solutions in one dimension and then two dimensional examples are presented. This comparison shows the ability of the method to generate theoretically expected simulations. In all figures the analytical solution is indicated by ‘‘Exact’’. The Convection-Diffusion equation represented by CD, the DiffusionReaction equation by DR, Finite Difference method by FD, Finite

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

Pe=20

Da=500

1

1

0.8

Exact

0.8

0.6

MFEM

T*

T*

1.2

0.6 Exact

0.4

0

FEM ENFEM

0.4 0.2

FD MFD

0.2

7

0

0

0.2

0.6

0.4

0.8

1

-0.2

0

0.2

0.4

0.6

0.8

X*

1

X*

Fig. 12. Comparison of normal FEM and Modified vrsions for DR equation in one dimension Da ¼500. Dissipation case.

Fig. 9. Comparison of normal FD and Modified FD for CD equation in one dimension Pe ¼20. Non-equal mesh size.

Pe=100

1.2

T*

T* Exact

0.4

-0.2

MFD

RFB

0

0.2

0.4

0.6

0.8

1

X*

-0.4

0

0.2

0.4

0.6

0.8

1

X*

Fig. 13. Comparison of normal FEM and Modified versions with RFB method for DR equation in one dimension Da ¼ 1000. Dissipation case. Non-equal mesh size.

Fig. 10. Comparison of normal FD and Modified FD for CD equation in one dimension Pe ¼100. Non-equal mesh size.

EXACT

3

FEM

0.8 MFEM

0.6

ENFEM

T*

1 0.4

0.6

0.8

FEM ENFEM

RFB

0.2 0 -0.2

0.2

Exact MFEM

0.4

2

0

Pe=100

1

Da=-500

4

T*

ENFEM

0.2 0

FD

0.2

1

0

0.2

0.4

0.6

0.8

1

-0.4 -0.6

-2 -3

MFEM

0.4

0.6

-1

FEM

0.6

0.8

0

Exact

0.8

1

0

Da=1000

1

X*

-0.8 X*

-4 Fig. 11. Comparison of normal FEM and Modified vrsions for DR equation in one dimension Da ¼  500. Production case.

element method by FEM, the modified coeficient or enriched versions of the proposed method is indicated by prefix ‘‘M’’ or ‘‘EN’’, respectively. Figs. 3 and 4 The results for the FD method in one dimension are represented in Figs. 5–10. Figs. 5 and 6 show the results for the DR equation using FD method in dissipation case. In both figures the

Fig. 14. Comparison of normal FEM and Modified versions with RFB method for CD equation in one dimension Pe ¼100.

normal FD is somewhat over-diffusive while the modified FD is quite accurate. The modification of the parameters of DR equation is capable of producing accurate solution in equal and non-equal size meshes. In production case where the solution is sinusoidal the normal FD has failed while the modified version can produce accurate solutions (Figs. 7 and 8). Figs. 9 and 10 show the results for CD equation. In Pe ¼20 and Pe ¼100 the proposed method has generated accurate results with non-equal meshes while the normal FD is over-diffusive.

8

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

Da=1000

1

0.4

0.2

T*

T*

0.4 0

0.2

-0.2

-0.2

0

0.2

0.4

0.6

0.8

0.4

-0.6

FEM

ENFEM

1

Da=2000

0.9

Exact

0.8

FD MFD

0.7

0.5 0

0.2

0.4

0.6

0.8

1

T*

0.6 0.5 0.4 0.3

-1

0.2

-1.5

0.1

-2

0

X*

0

0.2

0.4

0.6

1

FEM

0.8

MFEM

4 2

0.2

0

0

-2

0.4

0.6

0.8

1

X*

-0.6 Fig. 17. Comparison of normal FEM and Modified version for CD equation in two dimensions Pe ¼50. BLNR1BC1, yn ¼ 0.5.

Figs. 11–14 show the results for the FEM method. Fig. 11 shows the inability of normal FEM in predicting accurate solution for DR equation while modified FEM and Enriched FEM are capable of producing accurate solution for production case. The same result is observed in Fig. 12 for dissipation case. In Fig. 13 the result of Modified FEM and Enriched FEM are compared with residual free method for the DR equation and all of them generate accurate solutions while normal FEM is not capable to produce stable results. For the CD equation as Fig. 14 shows the Modified FEM and Enriched FEM can produce accurate solutions similar to RFB method while normal FEM has failed.

T*

0.4

0.2

Da=-200

6

0.6

-0.2 0

1

Fig. 19. Comparison of normal FD and Modified version for DR equation in two dimensions,BLNR1BC1, Da¼ 2000.yn ¼ 0.5.Dissipation case.

PE=50

1.2

0.8

X*

Fig. 16. Comparison of normal FEM and Modified version for DR equation in two dimensions. Production case. Da ¼  1000. BLNR1BC2.

T*

0.8

X*

1

1

-0.4

0.6

Fig. 18. Comparison of normal FEM and Modified version for CD equat,ion in two dimensions Pe ¼50. BLNR2BC1, yn ¼0.4.

Da=-1000

1.5

T*

0.2

X*

2

-0.5

0

-0.4

1

Fig. 15. Comparison of normal FEM and Modified version for DR equation in two dimensions. Dissipation case. Da ¼ 1000. BLNR1BC1, yn ¼ 0.5.

0

ENFEM

0.6

ENFEM

0.6

0

FEM

0.8

FEM

0.8

Pe=50

1

0

0.2

0.4

0.6

0.8

1

-4 -6

FD MFD

-8 -10

X*

Fig. 20. Comparison of normal FD and Modified version for DR equation in two dimensions BLNR1BC1 Da ¼  200. yn ¼0.5. Production case.

Figs. 15–18 show the results for FEM method. With equal size meshes the proposed method generates stable solution for DR equation with BC1 boundary conditions (see Fig. 15). In the production case the proposed method produces accurate solution with BC2 boundary conditions (see Fig. 16). Figs. 17 and 18 show the results for CD equation in two dimensions for the BC1 boundary conditions. For both mesh schemes the normal FEM produces unstable solutions.

M. Parvazinia / Finite Elements in Analysis and Design 58 (2012) 1–9

Pe=50

6 FD MFD

4 2

T*

0 -2

0

0.2

0.4

0.6

0.8

1

-4 -6

X*

-8

9

convection-diffusion equations with any level of discretisation while the solution is stable and accurate. The simulation results show that the method is capable of producing stable solutions with coarse meshes and non-equal mesh sizes (Figs. 6,9,10,13) in convection or reaction dominant cases. The results of the proposed method were compared with residual free bubble function method and performed the same but the proposed method is flexible enough to be applied on different numerical schemes. Here in this work it is applied on finite element and finite difference schemes. The method applied in two dimensional cases and therefore, it can possibly be applied on practical problems.

-10 Fig. 21. Comparison of normal FD and Modified version for CD equation in two dimensions BLNR1BC1 Pe¼ 50. yn ¼ 0.5.

PE=50

1.5 FD MFD

1

T*

0.5 0 -0.5

0

0.2

0.4

0.6

0.8

1

-1 -1.5

X*

Fig. 22. Comparison of normal FD and Modified version for CD equation in two dimensions BLNR1BC1 Pe¼ 50.yn ¼0.9.

Figs. 19–22 show the comparison between normal and modified FD methods in two dimensions using BC1 boundary conditions. For DR equation the normal FD is over-diffusive with respect to Fig. 19. In production case the normal FD is unstable (see Fig. 20). Figs. 21 and 22 show that, the normal FD is quite unstable to solve CD equation in both cross sections.

5. Conclusions A method based on an elemental adjustment technique was proposed. The method presented to solve Diffusion-reaction and

References [1] 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, Transp. Porous Media 63 (2006) 71–90. [2] Y. Effendiev, T.Y. Hou, Multiscale Finite Element Methods Theory and Applications, Springer, 2009. [3] F. Brezzi, M. Bristeau, L.P. Franca, M. Mallet, G. Roge, A relationship between stablized finite element methods and the Galerkin method with bubble functions, Comput. Meth. Appl. Mech. Eng. 96 (1992) 117–129. [4] C. Baiocchi, F. Brezzi, L.P. Franca, Virtual bubbles and Galerkin-least-squares type methods, Comput. Meth. Appl. Mech. Eng. 105R (1993) 125–141. [5] F. Brezzi, L.P. Franca, T.J.R. Hughes, A. Russo, b ¼ g, Comput. Meth. Appl. Mech. Eng. 145 (1997) 339–392. [6] M. Parvazinia, V. Nassehi, Multiscale finite element modelling of diffusionreaction equation using bubble function with bilinear and triangular elements, Comput. Meth. Appl. Mech. Eng. 196 (2007) 1095–1107. [7] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residual free bubbles for Helemholtz equation, Comput. Meth. Appl. Mech. Eng. 40 (1997) 4003–4009. [8] L.P. Franca, A. Russo, Deriving upwinding, mass lumping and selective reduced integration by residual free bubbles, Appl. Math. Lett. 9 (1996) 83–88. [9] L.P. Franca, A. Russo, Unlocking with residual-free bubbles, Comput. Meth. Appl. Mech. Eng. 142 (1997) 361–364. [10] L.P. Franca, A. Russo, Mass lumping emanating from residual free bubbles, Comput. Meth. Appl. Mech. Eng. 142 (1997) 353–360. [11] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residual-free bubbles for advective–diffusive equations, Comput. Meth. Appl. Mech. Eng. 166 (1998) 25–33. [12] V. Nassehi, M. Parvazinia, Finite Element Modelling of Multiscale Transport Phenomena, Imperial College Press, 2010. [13] T.J.R. Hughes, Multi scale phenomena, Green’s functions, the Dirichlet -toNeumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Meth. Appl. Mech. Eng. 127 (1995) 381–401. [14] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Meths. Appl. Mech. Eng. 32 (1982) 199–259. [15] M. Parvazinia, A multiscale finite element method for the solution of transport equations, Finite Elem. Anal. Des. 47 (2011) 211–219.