stabilized finite element method for the advection
Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018 www.elsevier.com/locate/cma
A multiscale/stabilized finite element method for the advection–diffusion equation A. Masud *, R.A. Khurram Department of Civil and Materials Engineering, The University of Illinois at Chicago, 842 West Taylor Street, 2095 Engineering Research Facility, Chicago, IL 60607-7023, USA Received 2 January 2003; received in revised form 9 September 2003; accepted 8 December 2003
Abstract This paper presents a multiscale method that yields a stabilized finite element formulation for the advection–diffusion equation. The multiscale method arises from a decomposition of the scalar field into coarse (resolved) scale and fine (unresolved) scale. The resulting stabilized formulation possesses superior properties like that of the SUPG and the GLS methods. A significant feature of the present method is that the definition of the stabilization term appears naturally, and therefore the formulation is free of any user-designed or user-defined parameters. Another important ingredient is that since the method is residual based, it satisfies consistency ab initio. Based on the proposed formulation, a family of 2-D elements comprising 3 and 6 node triangles and 4 and 9 node quadrilaterals has been developed. Numerical results show the good performance of the method on uniform, skewed as well as composite meshes and confirm convergence at optimal rates. 2004 Elsevier B.V. All rights reserved.
1. Introduction The advection–diffusion equation governs several important phenomena in physics and engineering, and also serves as a vehicle to study the more advanced Navier–Stokes equations. Therefore it has been the focus of intense research for quite some time (see e.g., [7,16,18,22] and references therein). For the advection dominated case this equation becomes hyperbolic and develops sharp features in the solution. Classical numerical methods for the advection–diffusion equation result in non-convergent elements. Specifically, methods based on the standard Galerkin finite element approach lack stability that manifests itself in terms of nonphysical oscillations. In order to correct the deficiencies in the standard Galerkin approach, Hughes and colleagues introduced streamline-upwind-Petrov–Galerkin (SUPG) technique [7,16]. This technique turned out to be the fore-runner of a new class of stabilization schemes, called the Galerkin/least-square (GLS) stabilization methods [18]. In the context of the advection–diffusion equation, the essential feature of this method is the stabilization of the advection operator without
*
Corresponding author. Fax: +1-312-996-2426. E-mail address: [email protected] (A. Masud).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.12.047
1998
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
upsetting consistency or compromising accuracy. The key idea in the GLS stabilization approach is to augment the Galerkin finite element formulation with a least-squares form of the residuals that are based on the corresponding Euler–Lagrange equations. These least-squares integrals contain stabilization parameters that are designed so that the method achieves exact solution in the case of one-dimensional model problems. A generalization of the nodaly exact stabilization parameters in 1-D to the case of nodaly exact parameters in multi-dimensions has been a challenging task and a literature review reveals various attempts in designing the optimal stabilization parameters [12,13]. GLS stabilization was followed by the ‘‘unusual stabilized methods’’ introduced by Franca and coworkers [8–12]. Concurrently, another class of stabilized methods that was based on the idea of augmenting the Galerkin method with virtual bubble functions was introduced by Brezzi and coworkers [1–6]. In [14,17], Hughes revisited the origins of the stabilization schemes from a variational multiscale approach. In the HughesÕ variational multiscale (HVM) method different stabilization techniques appear as special cases of the underlying sub-grid scale modelling concept. In this paper we follow this line of thought and present a multiscale/stabilized finite element formulation for the advection–diffusion equation. The new formulation, hereafter termed as the HVM formulation, has improved stability properties as compared to the Galerkin form of the problem. The main idea in the present method is a multiscale decomposition of the scalar field in the coarse and the fine scales. This decomposition is based on the assumption that the fine scale structures of the solution may not be captured by a given mesh. However, the influence of the fine scales on the coarse scales may not be negligible and therefore it must be accounted for. It is important to note that within the context of the Galerkin methods, the only possible way to capture the fine scale features is to successively refine the mesh, i.e., reduce the characteristic length scale of discretization. However this line of thought is extremely limited in its scope because it can quickly exhaust any given computational platform. In fact the question one needs to consider is as follows: Is it possible to capture all the essential features in the solution, both coarse and fine, while still using cruder meshes? In other words, is it possible to develop highly accurate numerical schemes on cruder discretizations? Taking this line of thought, which is based on the notion of existence of the fine scales in the problem, we proposed a stabilized formulation for the Darcy flow [19], the linearized incompressible Navier–Stokes equations [20], and for the convective–diffusive heat transfer equations [21]. The present paper also aims at taking this viewpoint and presents the HVM method for the advection–diffusion equation. However, in this paper we do not deal with either the discontinuity-capturing operators or the nonlinear methods that are at times employed to model sharp layers in the solution. An outline of the paper is as follows: Section 2 presents the governing equations and the standard Galerkin form. Emphasis in the paper is on the variational multiscale approach for the development of a new stabilized form, which is presented in Section 3. Section 4 presents the numerical results, and conclusions are drawn in Section 5.
2. The advection–diffusion equation Let X Rnsd be an open bounded region with piecewise smooth boundary C. We assume nsd P 2. The advection–diffusion equation is a rv jDv ¼ f v ¼ w on C;
in X;
ð1Þ ð2Þ
where v is the unknown scalar field, aðxÞ is the given flow velocity which is assumed solenoidal, i.e., r a ¼ 0 in X, j ¼ jðxÞ > 0 represents diffusivity, f ðxÞ is the prescribed source function, and w is the boundary condition.
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
1999
2.1. The standard weak form Let V H 1 ðXÞ \ C 0 ðXÞ denote the space of trial solutions and weighting functions for the unknown scalar field. The weak form is obtained by multiplying the governing equation with an admissible weighting function and integrating it over the domain, resulting in the following form. ðw; a rvÞ þ ðrw; jrvÞ ¼ ðw; f Þ; ð3Þ R where ð; Þ ¼ X ðÞdX is the L2 ðXÞ inner product and w is the weighting function for v. The standard Galerkin form is obtained by substituting the discrete counterparts for the weighting functions and the trial solutions, i.e., wh and vh , respectively. However, in the rest of the paper we have dropped the superposed h to keep the notation simple. Remark 1. The standard Galerkin form suffers from numerous technical drawbacks [7,16,18]. The most significant being the lack of stability for the advection dominated problems. In order to overcome this issue, stabilized methods such as the SUPG and the GLS methods have been developed. Remark 2. The trial solutions and the weighting functions are usually piecewise polynomials and can only capture variability at spatial scales that are larger than the characteristic length scale of the mesh. All spatial features that are smaller than the element size and represent the sub-grid variability are neglected in the standard Galerkin approach. The lack of stability of the Galerkin method for the advection dominated problems can be understood in this context. As the sub-grid scales are not captured adequately, rather they are completely ignored in the Galerkin method, their effect propagates to larger scales and deteriorate the coarse-scale calculations.
3. The variational multiscale method 3.1. Multiscale decomposition In this section we present the HughesÕ variational multiscale approach [14] for the development of a stabilized formulation for the advection–diffusion equation. A fundamental principle of multiscale approach is to acknowledge the presence of the fine scales that cannot be captured by a given spatial discretization. This is particularly important for the advection dominated problems where the solution develops sharp features that require an impractical grid resolution. Therefore, from the outset, we work with the premise of the existence of multiple scales in the problem. We consider the bounded domain X discretized into non-overlapping regions Xe (element domains) with boundaries Ce , e ¼ 1; 2; . . . ; numel , such that X¼
n[ umel
Xe :
ð4Þ
e¼1
We denote the union of element interiors and element boundaries by X0 and C0 , respectively. X0 ¼
n[ umel
ðintÞXe
ðelem: interiorsÞ;
ð5Þ
e¼1
C0 ¼
n[ umel e¼1
Ce
ðelem: boundariesÞ:
ð6Þ
2000
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
We assume an overlapping sum decomposition of the scalar field into the coarse scales or the resolvable scales and the fine scales or the sub-grid scales. vðxÞ ¼ vðxÞ þ v0 ðxÞ : |{z} |ffl{zffl} coarse scale
ð7Þ
fine scale
Likewise, we assume an overlapping sum decomposition of the weighting function into the coarse and the fine scale components indicated as w and w0 , respectively. wðxÞ ¼ wðxÞ þ w0 ðxÞ : |ffl{zffl} |fflffl{zfflffl} coarse scale
ð8Þ
fine scale
We further make an assumption that the sub-grid scales although non-zero within the elements, vanish identically over the element boundaries. v0 ¼ w0 ¼ 0 on C0 :
ð9Þ
Remark 3. Another viewpoint to look at this decomposition is as follows. Following (7) we can write v ¼ vh þ ve where vh is the solution that can be obtained with a good numerical scheme on a given mesh, and ve is the part of the solution that is lost because its scale is smaller than the characteristic length scale of the discretization. In other words ve represents the error in the solution. The method presented below can be considered as a systematic procedure whereby one can rebuild the error term in the weak form of the problem. This procedure then automatically yields a stabilized form of the problem that shows higher accuracy on cruder discretizations. We now introduce the appropriate spaces of functions for the coarse and the fine scale fields and specify a direct sum decomposition on these spaces. V ¼ V V0 ;
ð10Þ
where V in (10) is the space of trial solutions and weighting functions for the coarse scale velocity field and is identified with the standard finite element space. V ¼ fvjv 2 H01 ðXÞ; vðXe Þ ¼ Pk ðXe Þg; k
e
ð11Þ e
where P ðX Þ denotes complete polynomials of order k over X . On the other hand, various characterizations of V0 are possible, subject to the restriction imposed by the stability of the formulation that requires V and V0 to be linearly independent. Consequently, in the discrete case V0 can contain various finite dimensional approximations, e.g., bubble functions or p-refinements, that satisfy (9). V0 ¼ fv0 jv0 ¼ 0 on C0 g:
ð12Þ
Furthermore, because of (9), w0 2 V0 . 3.2. The multiscale variational problem We now substitute the trial solutions (7) and the weighting functions (8) in the standard variational form (3), and this becomes the point of departure from the conventional Galerkin formulations. ðw þ w0 ; a rðv þ v0 ÞÞ þ ðrðw þ w0 Þ; jrðv þ v0 ÞÞ ¼ ðw þ w0 ; f Þ:
ð13Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2001
With suitable assumptions on the fine scale field, as stipulated in (9), and employing the linearity of the weighting function slot, we can split the problem into the coarse and the fine scale parts, indicated as W and W0 , respectively. The coarse scale sub-problem W can be written as ðw; a rðv þ v0 ÞÞ þ ðrw; jrðv þ v0 ÞÞ ¼ ðw; f Þ:
ð14Þ
0
The fine scale sub-problem W can be written as ðw0 ; a rðv þ v0 ÞÞ þ ðrw0 ; jrðv þ v0 ÞÞ ¼ ðw0 ; f Þ:
ð15Þ
The general idea at this point is to solve the fine scale problem to obtain the fine scale solution v0 . This solution can then be substituted in the coarse scale problem (14), thereby eliminating the fine scales, yet retaining their effect. 3.3. The fine scale sub-problem (W0 ) Let us now consider the fine scale part of the weak form W0 , which, because of the assumption on the fine scale field, is defined over X0 . Exploiting linearity of the solution slot in (15) we have ðw0 ; a rvÞX0 þ ðw0 ; a rv0 ÞX0 þ ðrw0 ; jrvÞX0 þ ðrw0 ; jrv0 ÞX0 ¼ ðw0 ; f ÞX0 ; Pnumel R where ð; ÞX0 ¼ e¼1 ðÞ dX is the L2 ðXe Þ inner product. Xe Applying integration-by-parts to the third term on the left hand side and employing (9) ðw0 ; a rvÞX0 þ ðw0 ; a rv0 ÞX0 ðw0 ; jDvÞX0 þ ðrw0 ; jrv0 ÞX0 ¼ ðw0 ; f ÞX0 :
ð16Þ
ð17Þ
Rearranging terms, the fine scale problem reduces to ðw0 ; a rv0 ÞX0 þ ðrw0 ; jrv0 ÞX0 ¼ ðw0 ; f a rv þ jDvÞX0 :
ð18Þ
To crystallize ideas, and without loss of generality, we assume that the fine scales are represented via bubbles over element domains, i.e., v0 jXe ¼ be1 v0e on Xe ;
ð19Þ
w0 jXe ¼ be2 w0e on Xe ;
ð20Þ
where be1 and be2 represent the bubble shape functions for the fine scale trial solutions and the fine scale weighting functions, respectively. Furthermore, v0e and w0e represent the coefficients for the fine scale trial solutions and the weighting functions, respectively. Substituting (19) and (20) in the fine scale problem (18) we get ðbe2 w0e ; a rbe1 v0e Þ þ ðrbe2 w0e ; jrbe1 v0e Þ ¼ ðbe2 w0e ; f a rv þ jDvÞ: Taking the constant coefficients we get v0e ¼
w0e
and
v0e
ð21Þ
out of the integral expressions and employing arbitrariness of w0e ,
1 ðbe ; a rv jDv f Þ: ½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ 2
We now reconstruct the fine scale field via recourse to (19) 1 0 e e v ðxÞ ¼ b1 ðb ; a rv jDv f Þ : ½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ 2
ð22Þ
ð23Þ
In order to keep the presentation simple, and for the case where the residual of the coarse scales over element interiors can be considered constant, we can simplify Eq. (23) as follows:
2002
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
v0 ðxÞ ¼ s½a rv jDv f ; where s is defined as the stability function written as follows: Z 1 e s ¼ b1 be2 dX½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ :
ð24Þ
ð25Þ
Xe
Remark 4. Note that in the diffusion dominated case s is Oðh2 =jÞ, and in the advection dominated case s is Oðh=jajÞ. Consequently in both the limit cases s possesses the right order [10,18]. 3.4. The coarse scale problem (W) Employing the linearity of the solution slot in the coarse scale sub-problem (14) we get ðw; a rvÞ þ ðw; a rv0 Þ þ ðrw; jrvÞ þ ðrw; jrv0 Þ ¼ ðw; f Þ:
ð26Þ
Applying integration by parts to the second and the fourth terms, applying conditions (2) and (9), and then combining all the v0 terms we get ðw; a rvÞ þ ðrw; jrvÞ ða rw þ jDw; v0 Þ ¼ ðw; f Þ:
ð27Þ
Substituting v0 from (24) in (27) ðw; a rvÞ þ ðrw; jrvÞ þ ða rw þ jDw; sða rv jDvÞÞ ¼ ðw; f Þ þ ðsða rw þ jDwÞ; f Þ:
ð28Þ
3.5. The HVM form The stabilized form (28) is completely expressed in terms of the coarse/resolvable scales in the problem. Therefore, in order to keep the notation simple we drop the superposed bars and write the resulting HVM form as ðw; a rvÞ þ ðrw; jrvÞ þ ða rw þ jDw; sða rv jDvÞÞ ¼ ðw; f Þ þ ðsða rw þ jDwÞ; f Þ:
ð29Þ
Remark 5. It is remarkable that up to the definition of s, Eq. (28) derived above is identical to the modified GLS form first proposed in [10]. Remark 6. When compared with the standard Galerkin method, the multiscale approach involves additional integrals that are evaluated element wise. These additional terms represent the effects of the sub-grid scales that are being modeled in (29) in terms of the residuals of the coarse scales of the problem. Remark 7. The sub-grid scales are proportional to the residual of the coarse scales, i.e., it is a residual based method and therefore satisfies consistency ab initio. 3.6. The stability function The stability function is given in Eq. (25). We can write (25) as Z s ¼ be1 be2 dX^s; Xe
ð30Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2003
Fig. 1. 1-D bubbles for the weighting function and trial solution.
where ^s ¼ ½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ1 ;
ð31Þ
be1
be2
where is the bubble function for the trial solution and represents the bubble for the weighting function. In our numerical implementation we have used the absolute value of the advection velocity, i.e., ja1 j and ja2 j. Furthermore, in order to keep ^s positive, we take the absolute values of both the terms in the brackets. This restriction can however be relaxed by a better choice of the bubble for the fine scale weighting function w0 . The bubble functions for 1-D trial solutions v0 and weighting function w0 are presented in Fig. 1, and are given as follows. b1 ðxÞ ¼ ð1 x2 Þ; and
8 > > >
1 > > : ðx Xb Þ þ 1; 1 Xb
ð32Þ
1 6 x 6 Xb ; ð33Þ Xb 6 x 6 1;
where Xb is the location of the internal virtual node for the piecewise linear bubble for w0 . Remark 8. The bubble for the fine scale weighting function is piecewise continuous over Xe , but is not in C 1 [15]. Therefore this function is differentiable in a distributional sense only. We have employed this idea to develop piecewise bilinear functions for the fine scale weighting function for quadrilaterals and triangles. These functions are presented in Appendix A. 4. Numerical examples We have developed a family of 2-D elements comprising 3 and 6 node triangles and 4 and 9 node quadrilaterals (see Fig. 2). The bubble functions employed to compute ^s in (31) are given in Appendix A.
2004
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 2. A family of 2-D linear and quadratic elements.
In these functions Xb and Yb represent the location of the internal virtual node. In the numerical simulations presented here, we have used ðXb ; Yb Þ ¼ ð0:99; 0:99Þ for quadrilaterals, and ðXb ; Yb Þ ¼ ð0:34; 0:34Þ for triangles. 4.1. Rate of convergence study The first numerical simulation is a study of convergence rates on a model problem. The domain under consideration is a biunit square with origin at ðx; yÞ ¼ ð0; 0Þ. The exact solution is assumed of the following form. u ¼ sin
2px 2py sin : L L
The forcing function as evaluated from Eq. (1) is given below. 2p 2px 2py 2px 2py 4pj 2px 2py a1 cos sin cos þ sin sin f ðxÞ ¼ þ a2 sin : L L L L L L L L
ð34Þ
ð35Þ
In specifying the boundary-value problem, the forcing function is integrated over X while w ¼ ujC is prescribed nodally at the boundary. The rate of convergence study is divided into two parts. In the first part we show the rates for low Peclet number flows, and in the second part we present convergence rates for the high Peclet number flows. Throughout these convergence rate studies, we have taken j ¼ 1. Fig. 3a and b present the elevation plots of the forcing functions for low and high Peclet numbers, respectively. 4.1.1. Convergence study at low Peclet numbers This section presents the rate of convergence study for the (i) uniform, (ii) graded, and (iii) composite mesh cases, at low Peclet number. The convective velocity in this study is a1 ¼ a2 ¼ 10. The corresponding forcing function is presented in Fig. 3a. Case 1. Uniform meshes Fig. 4a and b show representative meshes for the linear quadrilateral and triangular elements, while Fig. 5a and b show the representative meshes for quadratic elements. For linear quadrilateral elements, the meshes employed consisted of 100, 400, 1600 and 6400 elements. The linear triangular element meshes consisted of exactly twice as many elements. For quadratic quadrilateral elements, the meshes employed consisted of 100, 400, 1600 and 6400 elements. Again, the quadratic triangular element meshes consisted of
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 3. Elevation plot of the forcing function for (a) low Pe (a1 ¼ a2 ¼ 10) and (b) high Pe (a1 ¼ a2 ¼ 106 ).
(a)
(b)
Fig. 4. Meshes of (a) 800 3-node triangles and (b) 400 4-node quadrilaterals.
(a)
(b)
Fig. 5. Meshes of (a) 200 6-node triangles and (b) 100 9-node quadrilaterals.
2005
2006
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 6. Convergence rates for linear elements (uniform meshes, low Pe).
twice as many elements. The element mesh parameter, h, is taken to be the edge length of the elements for the quadrilaterals, and the short-edge length for triangles. Convergence rates are presented in Figs. 6 and 7. Fig. 6 shows the L2 -norm and H 1 -seminorm of the scalar field for the linear quadrilateral and triangle. We see that optimal rates are achieved. Fig. 7 presents the corresponding rates for the quadratic elements. Once again we see optimal convergence rates. Fig. 8a and b present the elevation plots of the exact and the computed fields for a representative mesh of 4-node elements. We have also plotted the contours of the scalar field for 3-node triangles and the 9-node quadrilaterals in Fig. 9a and b, respectively. These plots present a qualitative comparison of the computed solution with the corresponding exact solution, and show a stable response of the elements. Similar plots were obtained for the 6-node triangles and are therefore not shown here. Case 2. Distorted structured meshes We have performed convergence study over structured, distorted, graded meshes composed of 3- and 6node triangles and 4- and 9-node quadrilaterals. The representative graded meshes for linear and quadratic
Fig. 7. Convergence rates for quadratic elements (uniform meshes, low Pe).
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2007
Fig. 8. (a) Elevation plot of the scalar field (exact solution). (b) Elevation plot of the scalar field for 4-node elements (20 · 20 element mesh).
Fig. 9. Contour plot of the scalar field for (a) 3-node triangles and (b) 9-node quadrilateral (low Pe).
elements are shown in Figs. 10a,b and 11a,b, respectively. The degree of distortion has been kept the same for all the element types. In this study the convective velocity a1 ¼ a2 ¼ 10, and the corresponding forcing function is presented in Fig. 3a. In all these cases convergence is attained at optimal rates and is presented in Figs. 12 and 13 for the linear and the quadratic elements, respectively. Fig. 14a and b presents the contours of the scalar field for the 4-node and the 9-node quadrilaterals, respectively. These plots present a qualitative comparison of the computed solution and show a stable response of the elements. Similar plots were obtained for the triangular elements and are therefore not shown here. Case 3. Composite meshes This study is interesting from a practical viewpoint where different element types are used in the same discretization. Fig. 15a presents a composite mesh composed of linear elements, i.e., combination of 3-node triangles and 4-node quadrilaterals in the same computational domain. Similarly, Fig. 15b presents the
2008
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
(a)
(b)
Fig. 10. Meshes of (a) 800 3-node triangles and (b) 400 4-node quadrilaterals (distorted meshes).
(a)
(b)
Fig. 11. Meshes of (a) 200 6-node triangles and (b) 100 9-node quadrilaterals (distorted meshes).
Fig. 12. Convergence rates for linear elements (distorted meshes, low Pe).
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2009
Fig. 13. Convergence rates for quadratic elements (distorted meshes, low Pe).
Fig. 14. Contour plot of the scalar field for (a) 4-node quadrilateral and (b) 9-node quadrilateral (distorted meshes, low Pe).
corresponding composite mesh for the quadratic elements. Fig. 16 presents the contour plot of the scalar field for the linear composite mesh. A similar plot was obtained for the composite mesh composed of quadratic elements and is therefore not shown here. The rate of convergence is a function of the order of interpolation polynomials and once again we get optimal rates as presented in Figs. 17 and 18. 4.1.2. Convergence study at high Peclet numbers In this section we repeat the rate of convergence study for the (i) uniform, (ii) graded, and (iii) composite mesh cases, at high Peclet number. The convective velocity in this study is a1 ¼ a2 ¼ 106 . The corresponding forcing function is presented in Fig. 3b. Case 1. Uniform meshes Meshes for this study are shown in Figs. 4 and 5. Convergence rates are presented in Figs. 19 and 20. Fig. 19 shows the L2 -norm and H 1 -seminorm of the scalar field for the linear quadrilateral and triangle. We see
2010
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
(a)
(b)
Fig. 15. Composite mesh of (a) linear elements and (b) quadratic elements.
Fig. 16. Contour plot of the scalar field for composite mesh of linear elements (low Pe).
Fig. 17. Convergence rates for composite mesh of linear elements (low Pe).
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2011
Fig. 18. Convergence rates for composite mesh of quadratic elements (low Pe).
Fig. 19. Convergence rates for linear elements (uniform meshes, high Pe).
that optimal rates are achieved. Fig. 20 presents the corresponding rates for the quadratic elements that show some degradation in the rate because of the high Peclet number. Case 2. Distorted structured meshes The representative graded meshes for linear and quadratic elements are shown in Figs. 10a,b and 11a,b, respectively. Figs. 21 and 22 present the convergence rates. Once again we obtain optimal convergence for the low order elements. Fig. 23a and b presents the contours of the scalar field for the 4-node and the 9node quadrilaterals, respectively, and show a stable response of the elements. Case 3. Composite meshes The composite linear- and quadratic meshes are shown in Fig. 15a and b. The rate of convergence is a function of the order of interpolation polynomials and once again we obtain optimal convergence rates for
2012
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 20. Convergence rates for quadratic elements (uniform meshes, high Pe).
Fig. 21. Convergence rates for linear elements (distorted meshes, high Pe).
the linear composite mesh composed of 3-node triangles and 4-node quadrilaterals (see Fig. 24). The convergence rates for the quadratic mesh are shown in Fig. 25. 4.2. Advection in a rotating flow field This numerical simulation is a test of very high Peclet number flows and is used to assess solutions which are essentially purely advective in nature. The problem is defined on a unit square of coordinates 0:5 6 x, y 6 þ 0:5, where the flow velocity components are given by a1 ¼ y;
a2 ¼ x:
ð36Þ
Along the external boundary v ¼ 0 and along the internal boundary (OA), v ¼ 12ðcosð4py þ pÞ þ 1Þ;
0:5 6 y 6 0:
ð37Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2013
Fig. 22. Convergence rates for quadratic elements (distorted meshes, high Pe).
Fig. 23. Contour plot of the scalar field for (a) 4-node quadrilateral and (b) 9-node quadrilateral (distorted meshes, high Pe).
A schematic diagram of the problem statement is shown in Fig. 26. The diffusivity is j ¼ 106 . A uniform mesh with 30 · 30 elements is employed. Elevation plots are shown in Fig. 27a and b for the HVM and the Galerkin methods, respectively. This problem has a smooth exact solution and therefore both methods perform well. There are however small amplitude oscillations in the Galerkin method that can be seen in Fig. 27b. 4.3. Thermal boundary layer problem This problem can be viewed as the simulation of a thermal boundary layer on a fully developed flow between two parallel plates, where the bottom plate is fixed and the top plate is moving with a unit velocity. The computational domain has dimensions Lx ¼ 1:0 and Ly ¼ 0:5. Schematic diagram of the problem is given in Fig. 28. Boundary conditions are given as follows.
2014
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 24. Convergence rates for composite mesh of linear elements (high Pe).
Fig. 25. Convergence rates for composite mesh of quadratic elements (high Pe).
u¼1
x ¼ 0; y ¼ 0:5;
u ¼ 0;
y ¼ 0;
u ¼ 2y;
x ¼ 1;
0 < y 6 0:5; 0 6 x 6 1;
ð38Þ
0 6 x 6 1:0;
ð39Þ
0 < y < 0:5:
ð40Þ
The flow components are a1 ¼ 2y;
a2 ¼ 0
in X
ð41Þ
and the diffusivity is j ¼ 7 104 . Taking the top plate velocity as the characteristic flow velocity, we have a Peclet number Pe ¼ uLy =j ¼ 714 for this flow.
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 26. Schematic diagram of the rotating flow field problem.
Fig. 27. Rotating flow problem. Elevation plot of the scalar field for the (a) HVM method and (b) Galerkin method.
Fig. 28. Schematic diagram of the thermal boundary layer problem.
2015
2016
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 29. Thermal boundary layer problem (a) HVM and (b) GLS.
We employ a mesh consisting of 21 equally spaced nodes in the x-direction, 11 nodes uniformly distributed in the interval 0 6 y 6 0:1 and the same number of nodes equally spaced on 0:1 6 y 6 0:5. This amounts to subdividing the domain into two different regions. Fig. 29a presents the elevation profile for the HVM method. For comparison we have plotted the elevation profile for the GLS method [10] in Fig. 29b. 5. Conclusions We have presented a new stabilized formulation for the advection–diffusion equation. The key point of the proposed formulation is a multiscale decomposition of the scalar field in the coarse and the fine scales. A significant feature of the present method is that it is free of user-defined stability parameters. The multiscale approach has several advantages over other stabilized formulations; (i) the stabilization term arises naturally, (ii) it is not restricted to a particular sub-grid model, and (iii) incorporates the effect of the fine (sub-grid) scales onto the coarse (grid) scales. When compared with standard Galerkin method, the multiscale method involves additional integrals that are evaluated element wise and incorporate the effect of the sub-grid scales on the coarse scales. The sub-grid scales are proportional to the residual of the coarse scales, i.e., the method is residual based and therefore, it is automatically consistent. We have developed a family of triangular and quadrilateral elements based on the proposed formulation. Numerical results show the good performance of the method and yield optimal convergence rates in the norms considered. Acknowledgements Support for this work was provided by the ONR grant N00014-02-1-0143. This support is gratefully acknowledged. The authors wish to thank the anonymous reviewers for helpful comments. Appendix A. Bubble functions for the 2-D cases A.1. Quadrilaterals The bubble for the fine scale trial solution is the usual quadratic bubble defined as b1 ðx; yÞ ¼ ð1 x2 Þð1 y 2 Þ:
ðA:1Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2017
Fig. 30. Schematic diagram showing the domains for the piecewise functions for (a) quadrilaterals and (b) triangles.
For the fine scale weighting function, the element domain is divided into four regions (see Fig. 30a) and the piecewise functions on these regions are given as follows. 8 x > > ð1 þ yÞ for x; y in region 1; > > > Xb ð1 þ Yb Þ > > > > > y > > > < Y ð1 X Þ ð1 xÞ for x; y in region 2; b b ðA:2Þ b2 ðx; yÞ ¼ > x > > ð1 yÞ for x; y in region 3; > > > Xb ð1 Yb Þ > > > > > y > > ð1 þ xÞ for x; y in region 4; : Yb ð1 þ Xb Þ where Xb and Yb represent the location of the internal virtual node in Xe . A.2. Triangles The bubble for the fine scale trial solution is a quadratic bubble defined as b1 ðx; yÞ ¼ 27xyð1 x yÞ:
ðA:3Þ
For the fine scale weighting function, the element domain is divided into three regions (see Fig. 30b) and the piecewise functions on these regions are given as follows. 8 xy > > for x; y in region 1; > > Xb Y b > > > < 1xy for x; y in region 2; ðA:4Þ b2 ðx; yÞ ¼ 1 Xb Y b > > > > > > xy > : for x; y in region 3: Xb Y b References [1] C. Baiocchi, F. Bezzi, L.P. Franca, Virtual bubbles and Galerkin-least-squares type methods (Ga.L.S), Comput. Methods Appl. Mech. Engrg. 105 (1993) 125–141.
2018
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
[2] F. Brezzi, M.O. 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 (1) (1992) 117–129. R [3] F. Brezzi, L.P. Franca, T.J.R. Hughes, A. Russo, b ¼ g, Comput. Methods Appl. Mech. Engrg. 145 (3–4) (1997) 329–339. [4] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residual-free bubbles for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 166 (1–2) (1998) 25–33. [5] F. Brezzi, P. Houston, D. Marini, E. Suli, Modeling subgrid viscosity for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1601–1610. [6] F. Brezzi, D. Marini, A. Russo, Application of the pseudo residual-free bubbles to the stabilization of convection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 166 (1998) 51–63. [7] 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 (1–3) (1982) 199–259. [8] L.P. Franca, C. Farhat, Bubble functions prompt unusual stabilized finite element methods, Comput. Methods Appl. Mech. Engrg. 123 (1–4) (1995) 299–308. [9] L.P. Franca, C. Farhat, M. Lesoinne, A. Russo, Unusual stabilized finite element methods and residual free bubbles, Int. J. Numer. Methods Fluids 27 (2) (1998) 159–168. [10] L.P. Franca, S.L. Frey, T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective–diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (2) (1992) 253–276. [11] L.P. Franca, A. Nesliturk, M. Stynes, On the stability of residual-free bubbles for convection–diffusion problems and their approximation by a two-level finite element method, Comput. Methods Appl. Mech. Engrg. 166 (1-2) (1998) 35–49. [12] L.P. Franca, F. Valentin, On an improved unusual stabilized finite element method for the advective–reactive–diffusive equation, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1785–1800. [13] I. Harari, L.P. Franca, S.P. Oliveira, Streamline design of stability parameters for advection–diffusion problems, J. Comput. Phys. 171 (1) (2001) 115–131. [14] T.J.R. Hughes, Multiscale 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) 387–401. [15] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewoods Cliffs, NJ, 1987, Dover edition, 2000. [16] T.J.R. Hughes, A.N. Brooks, Streamline-upwind/Petrov–Galerkin methods for advection dominated flows, in: Proceedings of the Third International Conference on Finite Element Methods in Fluid Flow, Banff, June 1980, pp. 283–292. [17] T.J.R. Hughes, G.R. Feijoo, M. Luca, Q. Jean-Baptiste, The variational multiscale method––a paradigm for computational mechanics, Comput. Methods Appl. Mech. Engrg. 166 (1–2) (1998) 3–24. [18] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (2) (1989) 173–189. [19] A. Masud, T.J.R. Hughes, A stabilized mixed finite element method for Darcy flow, Comput. Methods Appl. Mech. Engrg. 191/ 39–40 (2002) 4341–4370. [20] A. Masud, On a stabilized finite element formulation for incompressible Navier–Stokes equations, in: Proceedings of the Fourth US–Japan Conference on Computational Fluid Dynamics, Tokyo, Japan, May 28–30, 2002. [21] M. Ayub, A. Masud, A new stabilized formulation for convective–diffusive heat transfer, Numer. Heat Transfer Part B 44 (2003) 1–23. [22] T.E. Tezduyar, Y.J. Park, Discontinuity-capturing finite element formulations for nonlinear convection–diffusion–reaction equations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 307–325.
A multiscale/stabilized finite element method for the advection–diffusion equation A. Masud *, R.A. Khurram Department of Civil and Materials Engineering, The University of Illinois at Chicago, 842 West Taylor Street, 2095 Engineering Research Facility, Chicago, IL 60607-7023, USA Received 2 January 2003; received in revised form 9 September 2003; accepted 8 December 2003
Abstract This paper presents a multiscale method that yields a stabilized finite element formulation for the advection–diffusion equation. The multiscale method arises from a decomposition of the scalar field into coarse (resolved) scale and fine (unresolved) scale. The resulting stabilized formulation possesses superior properties like that of the SUPG and the GLS methods. A significant feature of the present method is that the definition of the stabilization term appears naturally, and therefore the formulation is free of any user-designed or user-defined parameters. Another important ingredient is that since the method is residual based, it satisfies consistency ab initio. Based on the proposed formulation, a family of 2-D elements comprising 3 and 6 node triangles and 4 and 9 node quadrilaterals has been developed. Numerical results show the good performance of the method on uniform, skewed as well as composite meshes and confirm convergence at optimal rates. 2004 Elsevier B.V. All rights reserved.
1. Introduction The advection–diffusion equation governs several important phenomena in physics and engineering, and also serves as a vehicle to study the more advanced Navier–Stokes equations. Therefore it has been the focus of intense research for quite some time (see e.g., [7,16,18,22] and references therein). For the advection dominated case this equation becomes hyperbolic and develops sharp features in the solution. Classical numerical methods for the advection–diffusion equation result in non-convergent elements. Specifically, methods based on the standard Galerkin finite element approach lack stability that manifests itself in terms of nonphysical oscillations. In order to correct the deficiencies in the standard Galerkin approach, Hughes and colleagues introduced streamline-upwind-Petrov–Galerkin (SUPG) technique [7,16]. This technique turned out to be the fore-runner of a new class of stabilization schemes, called the Galerkin/least-square (GLS) stabilization methods [18]. In the context of the advection–diffusion equation, the essential feature of this method is the stabilization of the advection operator without
*
Corresponding author. Fax: +1-312-996-2426. E-mail address: [email protected] (A. Masud).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.12.047
1998
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
upsetting consistency or compromising accuracy. The key idea in the GLS stabilization approach is to augment the Galerkin finite element formulation with a least-squares form of the residuals that are based on the corresponding Euler–Lagrange equations. These least-squares integrals contain stabilization parameters that are designed so that the method achieves exact solution in the case of one-dimensional model problems. A generalization of the nodaly exact stabilization parameters in 1-D to the case of nodaly exact parameters in multi-dimensions has been a challenging task and a literature review reveals various attempts in designing the optimal stabilization parameters [12,13]. GLS stabilization was followed by the ‘‘unusual stabilized methods’’ introduced by Franca and coworkers [8–12]. Concurrently, another class of stabilized methods that was based on the idea of augmenting the Galerkin method with virtual bubble functions was introduced by Brezzi and coworkers [1–6]. In [14,17], Hughes revisited the origins of the stabilization schemes from a variational multiscale approach. In the HughesÕ variational multiscale (HVM) method different stabilization techniques appear as special cases of the underlying sub-grid scale modelling concept. In this paper we follow this line of thought and present a multiscale/stabilized finite element formulation for the advection–diffusion equation. The new formulation, hereafter termed as the HVM formulation, has improved stability properties as compared to the Galerkin form of the problem. The main idea in the present method is a multiscale decomposition of the scalar field in the coarse and the fine scales. This decomposition is based on the assumption that the fine scale structures of the solution may not be captured by a given mesh. However, the influence of the fine scales on the coarse scales may not be negligible and therefore it must be accounted for. It is important to note that within the context of the Galerkin methods, the only possible way to capture the fine scale features is to successively refine the mesh, i.e., reduce the characteristic length scale of discretization. However this line of thought is extremely limited in its scope because it can quickly exhaust any given computational platform. In fact the question one needs to consider is as follows: Is it possible to capture all the essential features in the solution, both coarse and fine, while still using cruder meshes? In other words, is it possible to develop highly accurate numerical schemes on cruder discretizations? Taking this line of thought, which is based on the notion of existence of the fine scales in the problem, we proposed a stabilized formulation for the Darcy flow [19], the linearized incompressible Navier–Stokes equations [20], and for the convective–diffusive heat transfer equations [21]. The present paper also aims at taking this viewpoint and presents the HVM method for the advection–diffusion equation. However, in this paper we do not deal with either the discontinuity-capturing operators or the nonlinear methods that are at times employed to model sharp layers in the solution. An outline of the paper is as follows: Section 2 presents the governing equations and the standard Galerkin form. Emphasis in the paper is on the variational multiscale approach for the development of a new stabilized form, which is presented in Section 3. Section 4 presents the numerical results, and conclusions are drawn in Section 5.
2. The advection–diffusion equation Let X Rnsd be an open bounded region with piecewise smooth boundary C. We assume nsd P 2. The advection–diffusion equation is a rv jDv ¼ f v ¼ w on C;
in X;
ð1Þ ð2Þ
where v is the unknown scalar field, aðxÞ is the given flow velocity which is assumed solenoidal, i.e., r a ¼ 0 in X, j ¼ jðxÞ > 0 represents diffusivity, f ðxÞ is the prescribed source function, and w is the boundary condition.
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
1999
2.1. The standard weak form Let V H 1 ðXÞ \ C 0 ðXÞ denote the space of trial solutions and weighting functions for the unknown scalar field. The weak form is obtained by multiplying the governing equation with an admissible weighting function and integrating it over the domain, resulting in the following form. ðw; a rvÞ þ ðrw; jrvÞ ¼ ðw; f Þ; ð3Þ R where ð; Þ ¼ X ðÞdX is the L2 ðXÞ inner product and w is the weighting function for v. The standard Galerkin form is obtained by substituting the discrete counterparts for the weighting functions and the trial solutions, i.e., wh and vh , respectively. However, in the rest of the paper we have dropped the superposed h to keep the notation simple. Remark 1. The standard Galerkin form suffers from numerous technical drawbacks [7,16,18]. The most significant being the lack of stability for the advection dominated problems. In order to overcome this issue, stabilized methods such as the SUPG and the GLS methods have been developed. Remark 2. The trial solutions and the weighting functions are usually piecewise polynomials and can only capture variability at spatial scales that are larger than the characteristic length scale of the mesh. All spatial features that are smaller than the element size and represent the sub-grid variability are neglected in the standard Galerkin approach. The lack of stability of the Galerkin method for the advection dominated problems can be understood in this context. As the sub-grid scales are not captured adequately, rather they are completely ignored in the Galerkin method, their effect propagates to larger scales and deteriorate the coarse-scale calculations.
3. The variational multiscale method 3.1. Multiscale decomposition In this section we present the HughesÕ variational multiscale approach [14] for the development of a stabilized formulation for the advection–diffusion equation. A fundamental principle of multiscale approach is to acknowledge the presence of the fine scales that cannot be captured by a given spatial discretization. This is particularly important for the advection dominated problems where the solution develops sharp features that require an impractical grid resolution. Therefore, from the outset, we work with the premise of the existence of multiple scales in the problem. We consider the bounded domain X discretized into non-overlapping regions Xe (element domains) with boundaries Ce , e ¼ 1; 2; . . . ; numel , such that X¼
n[ umel
Xe :
ð4Þ
e¼1
We denote the union of element interiors and element boundaries by X0 and C0 , respectively. X0 ¼
n[ umel
ðintÞXe
ðelem: interiorsÞ;
ð5Þ
e¼1
C0 ¼
n[ umel e¼1
Ce
ðelem: boundariesÞ:
ð6Þ
2000
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
We assume an overlapping sum decomposition of the scalar field into the coarse scales or the resolvable scales and the fine scales or the sub-grid scales. vðxÞ ¼ vðxÞ þ v0 ðxÞ : |{z} |ffl{zffl} coarse scale
ð7Þ
fine scale
Likewise, we assume an overlapping sum decomposition of the weighting function into the coarse and the fine scale components indicated as w and w0 , respectively. wðxÞ ¼ wðxÞ þ w0 ðxÞ : |ffl{zffl} |fflffl{zfflffl} coarse scale
ð8Þ
fine scale
We further make an assumption that the sub-grid scales although non-zero within the elements, vanish identically over the element boundaries. v0 ¼ w0 ¼ 0 on C0 :
ð9Þ
Remark 3. Another viewpoint to look at this decomposition is as follows. Following (7) we can write v ¼ vh þ ve where vh is the solution that can be obtained with a good numerical scheme on a given mesh, and ve is the part of the solution that is lost because its scale is smaller than the characteristic length scale of the discretization. In other words ve represents the error in the solution. The method presented below can be considered as a systematic procedure whereby one can rebuild the error term in the weak form of the problem. This procedure then automatically yields a stabilized form of the problem that shows higher accuracy on cruder discretizations. We now introduce the appropriate spaces of functions for the coarse and the fine scale fields and specify a direct sum decomposition on these spaces. V ¼ V V0 ;
ð10Þ
where V in (10) is the space of trial solutions and weighting functions for the coarse scale velocity field and is identified with the standard finite element space. V ¼ fvjv 2 H01 ðXÞ; vðXe Þ ¼ Pk ðXe Þg; k
e
ð11Þ e
where P ðX Þ denotes complete polynomials of order k over X . On the other hand, various characterizations of V0 are possible, subject to the restriction imposed by the stability of the formulation that requires V and V0 to be linearly independent. Consequently, in the discrete case V0 can contain various finite dimensional approximations, e.g., bubble functions or p-refinements, that satisfy (9). V0 ¼ fv0 jv0 ¼ 0 on C0 g:
ð12Þ
Furthermore, because of (9), w0 2 V0 . 3.2. The multiscale variational problem We now substitute the trial solutions (7) and the weighting functions (8) in the standard variational form (3), and this becomes the point of departure from the conventional Galerkin formulations. ðw þ w0 ; a rðv þ v0 ÞÞ þ ðrðw þ w0 Þ; jrðv þ v0 ÞÞ ¼ ðw þ w0 ; f Þ:
ð13Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2001
With suitable assumptions on the fine scale field, as stipulated in (9), and employing the linearity of the weighting function slot, we can split the problem into the coarse and the fine scale parts, indicated as W and W0 , respectively. The coarse scale sub-problem W can be written as ðw; a rðv þ v0 ÞÞ þ ðrw; jrðv þ v0 ÞÞ ¼ ðw; f Þ:
ð14Þ
0
The fine scale sub-problem W can be written as ðw0 ; a rðv þ v0 ÞÞ þ ðrw0 ; jrðv þ v0 ÞÞ ¼ ðw0 ; f Þ:
ð15Þ
The general idea at this point is to solve the fine scale problem to obtain the fine scale solution v0 . This solution can then be substituted in the coarse scale problem (14), thereby eliminating the fine scales, yet retaining their effect. 3.3. The fine scale sub-problem (W0 ) Let us now consider the fine scale part of the weak form W0 , which, because of the assumption on the fine scale field, is defined over X0 . Exploiting linearity of the solution slot in (15) we have ðw0 ; a rvÞX0 þ ðw0 ; a rv0 ÞX0 þ ðrw0 ; jrvÞX0 þ ðrw0 ; jrv0 ÞX0 ¼ ðw0 ; f ÞX0 ; Pnumel R where ð; ÞX0 ¼ e¼1 ðÞ dX is the L2 ðXe Þ inner product. Xe Applying integration-by-parts to the third term on the left hand side and employing (9) ðw0 ; a rvÞX0 þ ðw0 ; a rv0 ÞX0 ðw0 ; jDvÞX0 þ ðrw0 ; jrv0 ÞX0 ¼ ðw0 ; f ÞX0 :
ð16Þ
ð17Þ
Rearranging terms, the fine scale problem reduces to ðw0 ; a rv0 ÞX0 þ ðrw0 ; jrv0 ÞX0 ¼ ðw0 ; f a rv þ jDvÞX0 :
ð18Þ
To crystallize ideas, and without loss of generality, we assume that the fine scales are represented via bubbles over element domains, i.e., v0 jXe ¼ be1 v0e on Xe ;
ð19Þ
w0 jXe ¼ be2 w0e on Xe ;
ð20Þ
where be1 and be2 represent the bubble shape functions for the fine scale trial solutions and the fine scale weighting functions, respectively. Furthermore, v0e and w0e represent the coefficients for the fine scale trial solutions and the weighting functions, respectively. Substituting (19) and (20) in the fine scale problem (18) we get ðbe2 w0e ; a rbe1 v0e Þ þ ðrbe2 w0e ; jrbe1 v0e Þ ¼ ðbe2 w0e ; f a rv þ jDvÞ: Taking the constant coefficients we get v0e ¼
w0e
and
v0e
ð21Þ
out of the integral expressions and employing arbitrariness of w0e ,
1 ðbe ; a rv jDv f Þ: ½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ 2
We now reconstruct the fine scale field via recourse to (19) 1 0 e e v ðxÞ ¼ b1 ðb ; a rv jDv f Þ : ½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ 2
ð22Þ
ð23Þ
In order to keep the presentation simple, and for the case where the residual of the coarse scales over element interiors can be considered constant, we can simplify Eq. (23) as follows:
2002
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
v0 ðxÞ ¼ s½a rv jDv f ; where s is defined as the stability function written as follows: Z 1 e s ¼ b1 be2 dX½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ :
ð24Þ
ð25Þ
Xe
Remark 4. Note that in the diffusion dominated case s is Oðh2 =jÞ, and in the advection dominated case s is Oðh=jajÞ. Consequently in both the limit cases s possesses the right order [10,18]. 3.4. The coarse scale problem (W) Employing the linearity of the solution slot in the coarse scale sub-problem (14) we get ðw; a rvÞ þ ðw; a rv0 Þ þ ðrw; jrvÞ þ ðrw; jrv0 Þ ¼ ðw; f Þ:
ð26Þ
Applying integration by parts to the second and the fourth terms, applying conditions (2) and (9), and then combining all the v0 terms we get ðw; a rvÞ þ ðrw; jrvÞ ða rw þ jDw; v0 Þ ¼ ðw; f Þ:
ð27Þ
Substituting v0 from (24) in (27) ðw; a rvÞ þ ðrw; jrvÞ þ ða rw þ jDw; sða rv jDvÞÞ ¼ ðw; f Þ þ ðsða rw þ jDwÞ; f Þ:
ð28Þ
3.5. The HVM form The stabilized form (28) is completely expressed in terms of the coarse/resolvable scales in the problem. Therefore, in order to keep the notation simple we drop the superposed bars and write the resulting HVM form as ðw; a rvÞ þ ðrw; jrvÞ þ ða rw þ jDw; sða rv jDvÞÞ ¼ ðw; f Þ þ ðsða rw þ jDwÞ; f Þ:
ð29Þ
Remark 5. It is remarkable that up to the definition of s, Eq. (28) derived above is identical to the modified GLS form first proposed in [10]. Remark 6. When compared with the standard Galerkin method, the multiscale approach involves additional integrals that are evaluated element wise. These additional terms represent the effects of the sub-grid scales that are being modeled in (29) in terms of the residuals of the coarse scales of the problem. Remark 7. The sub-grid scales are proportional to the residual of the coarse scales, i.e., it is a residual based method and therefore satisfies consistency ab initio. 3.6. The stability function The stability function is given in Eq. (25). We can write (25) as Z s ¼ be1 be2 dX^s; Xe
ð30Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2003
Fig. 1. 1-D bubbles for the weighting function and trial solution.
where ^s ¼ ½ðbe2 ; a rbe1 Þ þ ðrbe2 ; jrbe1 Þ1 ;
ð31Þ
be1
be2
where is the bubble function for the trial solution and represents the bubble for the weighting function. In our numerical implementation we have used the absolute value of the advection velocity, i.e., ja1 j and ja2 j. Furthermore, in order to keep ^s positive, we take the absolute values of both the terms in the brackets. This restriction can however be relaxed by a better choice of the bubble for the fine scale weighting function w0 . The bubble functions for 1-D trial solutions v0 and weighting function w0 are presented in Fig. 1, and are given as follows. b1 ðxÞ ¼ ð1 x2 Þ; and
8 > > >
1 > > : ðx Xb Þ þ 1; 1 Xb
ð32Þ
1 6 x 6 Xb ; ð33Þ Xb 6 x 6 1;
where Xb is the location of the internal virtual node for the piecewise linear bubble for w0 . Remark 8. The bubble for the fine scale weighting function is piecewise continuous over Xe , but is not in C 1 [15]. Therefore this function is differentiable in a distributional sense only. We have employed this idea to develop piecewise bilinear functions for the fine scale weighting function for quadrilaterals and triangles. These functions are presented in Appendix A. 4. Numerical examples We have developed a family of 2-D elements comprising 3 and 6 node triangles and 4 and 9 node quadrilaterals (see Fig. 2). The bubble functions employed to compute ^s in (31) are given in Appendix A.
2004
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 2. A family of 2-D linear and quadratic elements.
In these functions Xb and Yb represent the location of the internal virtual node. In the numerical simulations presented here, we have used ðXb ; Yb Þ ¼ ð0:99; 0:99Þ for quadrilaterals, and ðXb ; Yb Þ ¼ ð0:34; 0:34Þ for triangles. 4.1. Rate of convergence study The first numerical simulation is a study of convergence rates on a model problem. The domain under consideration is a biunit square with origin at ðx; yÞ ¼ ð0; 0Þ. The exact solution is assumed of the following form. u ¼ sin
2px 2py sin : L L
The forcing function as evaluated from Eq. (1) is given below. 2p 2px 2py 2px 2py 4pj 2px 2py a1 cos sin cos þ sin sin f ðxÞ ¼ þ a2 sin : L L L L L L L L
ð34Þ
ð35Þ
In specifying the boundary-value problem, the forcing function is integrated over X while w ¼ ujC is prescribed nodally at the boundary. The rate of convergence study is divided into two parts. In the first part we show the rates for low Peclet number flows, and in the second part we present convergence rates for the high Peclet number flows. Throughout these convergence rate studies, we have taken j ¼ 1. Fig. 3a and b present the elevation plots of the forcing functions for low and high Peclet numbers, respectively. 4.1.1. Convergence study at low Peclet numbers This section presents the rate of convergence study for the (i) uniform, (ii) graded, and (iii) composite mesh cases, at low Peclet number. The convective velocity in this study is a1 ¼ a2 ¼ 10. The corresponding forcing function is presented in Fig. 3a. Case 1. Uniform meshes Fig. 4a and b show representative meshes for the linear quadrilateral and triangular elements, while Fig. 5a and b show the representative meshes for quadratic elements. For linear quadrilateral elements, the meshes employed consisted of 100, 400, 1600 and 6400 elements. The linear triangular element meshes consisted of exactly twice as many elements. For quadratic quadrilateral elements, the meshes employed consisted of 100, 400, 1600 and 6400 elements. Again, the quadratic triangular element meshes consisted of
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 3. Elevation plot of the forcing function for (a) low Pe (a1 ¼ a2 ¼ 10) and (b) high Pe (a1 ¼ a2 ¼ 106 ).
(a)
(b)
Fig. 4. Meshes of (a) 800 3-node triangles and (b) 400 4-node quadrilaterals.
(a)
(b)
Fig. 5. Meshes of (a) 200 6-node triangles and (b) 100 9-node quadrilaterals.
2005
2006
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 6. Convergence rates for linear elements (uniform meshes, low Pe).
twice as many elements. The element mesh parameter, h, is taken to be the edge length of the elements for the quadrilaterals, and the short-edge length for triangles. Convergence rates are presented in Figs. 6 and 7. Fig. 6 shows the L2 -norm and H 1 -seminorm of the scalar field for the linear quadrilateral and triangle. We see that optimal rates are achieved. Fig. 7 presents the corresponding rates for the quadratic elements. Once again we see optimal convergence rates. Fig. 8a and b present the elevation plots of the exact and the computed fields for a representative mesh of 4-node elements. We have also plotted the contours of the scalar field for 3-node triangles and the 9-node quadrilaterals in Fig. 9a and b, respectively. These plots present a qualitative comparison of the computed solution with the corresponding exact solution, and show a stable response of the elements. Similar plots were obtained for the 6-node triangles and are therefore not shown here. Case 2. Distorted structured meshes We have performed convergence study over structured, distorted, graded meshes composed of 3- and 6node triangles and 4- and 9-node quadrilaterals. The representative graded meshes for linear and quadratic
Fig. 7. Convergence rates for quadratic elements (uniform meshes, low Pe).
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2007
Fig. 8. (a) Elevation plot of the scalar field (exact solution). (b) Elevation plot of the scalar field for 4-node elements (20 · 20 element mesh).
Fig. 9. Contour plot of the scalar field for (a) 3-node triangles and (b) 9-node quadrilateral (low Pe).
elements are shown in Figs. 10a,b and 11a,b, respectively. The degree of distortion has been kept the same for all the element types. In this study the convective velocity a1 ¼ a2 ¼ 10, and the corresponding forcing function is presented in Fig. 3a. In all these cases convergence is attained at optimal rates and is presented in Figs. 12 and 13 for the linear and the quadratic elements, respectively. Fig. 14a and b presents the contours of the scalar field for the 4-node and the 9-node quadrilaterals, respectively. These plots present a qualitative comparison of the computed solution and show a stable response of the elements. Similar plots were obtained for the triangular elements and are therefore not shown here. Case 3. Composite meshes This study is interesting from a practical viewpoint where different element types are used in the same discretization. Fig. 15a presents a composite mesh composed of linear elements, i.e., combination of 3-node triangles and 4-node quadrilaterals in the same computational domain. Similarly, Fig. 15b presents the
2008
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
(a)
(b)
Fig. 10. Meshes of (a) 800 3-node triangles and (b) 400 4-node quadrilaterals (distorted meshes).
(a)
(b)
Fig. 11. Meshes of (a) 200 6-node triangles and (b) 100 9-node quadrilaterals (distorted meshes).
Fig. 12. Convergence rates for linear elements (distorted meshes, low Pe).
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2009
Fig. 13. Convergence rates for quadratic elements (distorted meshes, low Pe).
Fig. 14. Contour plot of the scalar field for (a) 4-node quadrilateral and (b) 9-node quadrilateral (distorted meshes, low Pe).
corresponding composite mesh for the quadratic elements. Fig. 16 presents the contour plot of the scalar field for the linear composite mesh. A similar plot was obtained for the composite mesh composed of quadratic elements and is therefore not shown here. The rate of convergence is a function of the order of interpolation polynomials and once again we get optimal rates as presented in Figs. 17 and 18. 4.1.2. Convergence study at high Peclet numbers In this section we repeat the rate of convergence study for the (i) uniform, (ii) graded, and (iii) composite mesh cases, at high Peclet number. The convective velocity in this study is a1 ¼ a2 ¼ 106 . The corresponding forcing function is presented in Fig. 3b. Case 1. Uniform meshes Meshes for this study are shown in Figs. 4 and 5. Convergence rates are presented in Figs. 19 and 20. Fig. 19 shows the L2 -norm and H 1 -seminorm of the scalar field for the linear quadrilateral and triangle. We see
2010
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
(a)
(b)
Fig. 15. Composite mesh of (a) linear elements and (b) quadratic elements.
Fig. 16. Contour plot of the scalar field for composite mesh of linear elements (low Pe).
Fig. 17. Convergence rates for composite mesh of linear elements (low Pe).
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2011
Fig. 18. Convergence rates for composite mesh of quadratic elements (low Pe).
Fig. 19. Convergence rates for linear elements (uniform meshes, high Pe).
that optimal rates are achieved. Fig. 20 presents the corresponding rates for the quadratic elements that show some degradation in the rate because of the high Peclet number. Case 2. Distorted structured meshes The representative graded meshes for linear and quadratic elements are shown in Figs. 10a,b and 11a,b, respectively. Figs. 21 and 22 present the convergence rates. Once again we obtain optimal convergence for the low order elements. Fig. 23a and b presents the contours of the scalar field for the 4-node and the 9node quadrilaterals, respectively, and show a stable response of the elements. Case 3. Composite meshes The composite linear- and quadratic meshes are shown in Fig. 15a and b. The rate of convergence is a function of the order of interpolation polynomials and once again we obtain optimal convergence rates for
2012
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 20. Convergence rates for quadratic elements (uniform meshes, high Pe).
Fig. 21. Convergence rates for linear elements (distorted meshes, high Pe).
the linear composite mesh composed of 3-node triangles and 4-node quadrilaterals (see Fig. 24). The convergence rates for the quadratic mesh are shown in Fig. 25. 4.2. Advection in a rotating flow field This numerical simulation is a test of very high Peclet number flows and is used to assess solutions which are essentially purely advective in nature. The problem is defined on a unit square of coordinates 0:5 6 x, y 6 þ 0:5, where the flow velocity components are given by a1 ¼ y;
a2 ¼ x:
ð36Þ
Along the external boundary v ¼ 0 and along the internal boundary (OA), v ¼ 12ðcosð4py þ pÞ þ 1Þ;
0:5 6 y 6 0:
ð37Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2013
Fig. 22. Convergence rates for quadratic elements (distorted meshes, high Pe).
Fig. 23. Contour plot of the scalar field for (a) 4-node quadrilateral and (b) 9-node quadrilateral (distorted meshes, high Pe).
A schematic diagram of the problem statement is shown in Fig. 26. The diffusivity is j ¼ 106 . A uniform mesh with 30 · 30 elements is employed. Elevation plots are shown in Fig. 27a and b for the HVM and the Galerkin methods, respectively. This problem has a smooth exact solution and therefore both methods perform well. There are however small amplitude oscillations in the Galerkin method that can be seen in Fig. 27b. 4.3. Thermal boundary layer problem This problem can be viewed as the simulation of a thermal boundary layer on a fully developed flow between two parallel plates, where the bottom plate is fixed and the top plate is moving with a unit velocity. The computational domain has dimensions Lx ¼ 1:0 and Ly ¼ 0:5. Schematic diagram of the problem is given in Fig. 28. Boundary conditions are given as follows.
2014
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 24. Convergence rates for composite mesh of linear elements (high Pe).
Fig. 25. Convergence rates for composite mesh of quadratic elements (high Pe).
u¼1
x ¼ 0; y ¼ 0:5;
u ¼ 0;
y ¼ 0;
u ¼ 2y;
x ¼ 1;
0 < y 6 0:5; 0 6 x 6 1;
ð38Þ
0 6 x 6 1:0;
ð39Þ
0 < y < 0:5:
ð40Þ
The flow components are a1 ¼ 2y;
a2 ¼ 0
in X
ð41Þ
and the diffusivity is j ¼ 7 104 . Taking the top plate velocity as the characteristic flow velocity, we have a Peclet number Pe ¼ uLy =j ¼ 714 for this flow.
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 26. Schematic diagram of the rotating flow field problem.
Fig. 27. Rotating flow problem. Elevation plot of the scalar field for the (a) HVM method and (b) Galerkin method.
Fig. 28. Schematic diagram of the thermal boundary layer problem.
2015
2016
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
Fig. 29. Thermal boundary layer problem (a) HVM and (b) GLS.
We employ a mesh consisting of 21 equally spaced nodes in the x-direction, 11 nodes uniformly distributed in the interval 0 6 y 6 0:1 and the same number of nodes equally spaced on 0:1 6 y 6 0:5. This amounts to subdividing the domain into two different regions. Fig. 29a presents the elevation profile for the HVM method. For comparison we have plotted the elevation profile for the GLS method [10] in Fig. 29b. 5. Conclusions We have presented a new stabilized formulation for the advection–diffusion equation. The key point of the proposed formulation is a multiscale decomposition of the scalar field in the coarse and the fine scales. A significant feature of the present method is that it is free of user-defined stability parameters. The multiscale approach has several advantages over other stabilized formulations; (i) the stabilization term arises naturally, (ii) it is not restricted to a particular sub-grid model, and (iii) incorporates the effect of the fine (sub-grid) scales onto the coarse (grid) scales. When compared with standard Galerkin method, the multiscale method involves additional integrals that are evaluated element wise and incorporate the effect of the sub-grid scales on the coarse scales. The sub-grid scales are proportional to the residual of the coarse scales, i.e., the method is residual based and therefore, it is automatically consistent. We have developed a family of triangular and quadrilateral elements based on the proposed formulation. Numerical results show the good performance of the method and yield optimal convergence rates in the norms considered. Acknowledgements Support for this work was provided by the ONR grant N00014-02-1-0143. This support is gratefully acknowledged. The authors wish to thank the anonymous reviewers for helpful comments. Appendix A. Bubble functions for the 2-D cases A.1. Quadrilaterals The bubble for the fine scale trial solution is the usual quadratic bubble defined as b1 ðx; yÞ ¼ ð1 x2 Þð1 y 2 Þ:
ðA:1Þ
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
2017
Fig. 30. Schematic diagram showing the domains for the piecewise functions for (a) quadrilaterals and (b) triangles.
For the fine scale weighting function, the element domain is divided into four regions (see Fig. 30a) and the piecewise functions on these regions are given as follows. 8 x > > ð1 þ yÞ for x; y in region 1; > > > Xb ð1 þ Yb Þ > > > > > y > > > < Y ð1 X Þ ð1 xÞ for x; y in region 2; b b ðA:2Þ b2 ðx; yÞ ¼ > x > > ð1 yÞ for x; y in region 3; > > > Xb ð1 Yb Þ > > > > > y > > ð1 þ xÞ for x; y in region 4; : Yb ð1 þ Xb Þ where Xb and Yb represent the location of the internal virtual node in Xe . A.2. Triangles The bubble for the fine scale trial solution is a quadratic bubble defined as b1 ðx; yÞ ¼ 27xyð1 x yÞ:
ðA:3Þ
For the fine scale weighting function, the element domain is divided into three regions (see Fig. 30b) and the piecewise functions on these regions are given as follows. 8 xy > > for x; y in region 1; > > Xb Y b > > > < 1xy for x; y in region 2; ðA:4Þ b2 ðx; yÞ ¼ 1 Xb Y b > > > > > > xy > : for x; y in region 3: Xb Y b References [1] C. Baiocchi, F. Bezzi, L.P. Franca, Virtual bubbles and Galerkin-least-squares type methods (Ga.L.S), Comput. Methods Appl. Mech. Engrg. 105 (1993) 125–141.
2018
A. Masud, R.A. Khurram / Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018
[2] F. Brezzi, M.O. 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 (1) (1992) 117–129. R [3] F. Brezzi, L.P. Franca, T.J.R. Hughes, A. Russo, b ¼ g, Comput. Methods Appl. Mech. Engrg. 145 (3–4) (1997) 329–339. [4] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residual-free bubbles for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 166 (1–2) (1998) 25–33. [5] F. Brezzi, P. Houston, D. Marini, E. Suli, Modeling subgrid viscosity for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1601–1610. [6] F. Brezzi, D. Marini, A. Russo, Application of the pseudo residual-free bubbles to the stabilization of convection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 166 (1998) 51–63. [7] 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 (1–3) (1982) 199–259. [8] L.P. Franca, C. Farhat, Bubble functions prompt unusual stabilized finite element methods, Comput. Methods Appl. Mech. Engrg. 123 (1–4) (1995) 299–308. [9] L.P. Franca, C. Farhat, M. Lesoinne, A. Russo, Unusual stabilized finite element methods and residual free bubbles, Int. J. Numer. Methods Fluids 27 (2) (1998) 159–168. [10] L.P. Franca, S.L. Frey, T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective–diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (2) (1992) 253–276. [11] L.P. Franca, A. Nesliturk, M. Stynes, On the stability of residual-free bubbles for convection–diffusion problems and their approximation by a two-level finite element method, Comput. Methods Appl. Mech. Engrg. 166 (1-2) (1998) 35–49. [12] L.P. Franca, F. Valentin, On an improved unusual stabilized finite element method for the advective–reactive–diffusive equation, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1785–1800. [13] I. Harari, L.P. Franca, S.P. Oliveira, Streamline design of stability parameters for advection–diffusion problems, J. Comput. Phys. 171 (1) (2001) 115–131. [14] T.J.R. Hughes, Multiscale 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) 387–401. [15] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewoods Cliffs, NJ, 1987, Dover edition, 2000. [16] T.J.R. Hughes, A.N. Brooks, Streamline-upwind/Petrov–Galerkin methods for advection dominated flows, in: Proceedings of the Third International Conference on Finite Element Methods in Fluid Flow, Banff, June 1980, pp. 283–292. [17] T.J.R. Hughes, G.R. Feijoo, M. Luca, Q. Jean-Baptiste, The variational multiscale method––a paradigm for computational mechanics, Comput. Methods Appl. Mech. Engrg. 166 (1–2) (1998) 3–24. [18] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (2) (1989) 173–189. [19] A. Masud, T.J.R. Hughes, A stabilized mixed finite element method for Darcy flow, Comput. Methods Appl. Mech. Engrg. 191/ 39–40 (2002) 4341–4370. [20] A. Masud, On a stabilized finite element formulation for incompressible Navier–Stokes equations, in: Proceedings of the Fourth US–Japan Conference on Computational Fluid Dynamics, Tokyo, Japan, May 28–30, 2002. [21] M. Ayub, A. Masud, A new stabilized formulation for convective–diffusive heat transfer, Numer. Heat Transfer Part B 44 (2003) 1–23. [22] T.E. Tezduyar, Y.J. Park, Discontinuity-capturing finite element formulations for nonlinear convection–diffusion–reaction equations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 307–325.
