- Email: [email protected]
A new upwind function in stabilized finite element formulations, using linear and quadratic elements for scalar convection–diffusion problems
Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402 www.elsevier.com/locate/cma
A new upwind function in stabilized finite element formulations, using linear and quadratic elements for scalar convection–diffusion problems Eduardo Gomes Dutra do Carmo
a,*
, Gustavo Benitez Alvarez
b
a
b
Department of Nuclear Engineering, COPPE––Federal University of Rio de Janeiro, Ilha do Fund~ao, P.B. 68509, Rio de Janeiro, RJ 21945-970, Brazil LNCC––National Laboratory of Scientific Computation, Av. Getulio Vargas 333, P.B. 95113, Petropolis, RJ 25651-070, Brazil Received 15 April 2003; received in revised form 7 November 2003; accepted 7 January 2004
Abstract This article presents a study on the function of the Peclet number that appears in the upwind function of the GLS, CAU and SAUPG methods, using linear and quadratic triangular elements, as well as bilinear and quadratic quadrilateral elements. The numerical experiments indicate that in the case of elements with faces on the ‘‘outflow boundary’’, the upwind function is strongly dependent on the geometry and on the degree of the interpolate polynomial. A new function of the Peclet number is therefore proposed, to take these effects into consideration, creating a new method of stabilization for higher order elements. This new stabilized and accurate finite element formulation for convection dominated problems, based on the idea of the SAUPG method, uses the GLS formulation as starting point. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Stabilized FEM; Convection–diffusion; Stabilization
1. Introduction In general, the solution of diffusive-convective problems possesses boundary layers, that are small subregions where derivatives of the solution are very large. In this case, the classic finite element method or Galerkin method is inappropriate, because its numerical solution presents spurious oscillations. During the last two decades, various attempts to obtain new methods have been carried out to avoid such oscillations. These new methods are nowadays known as stabilized finite element methods. References about some of them can be found in [14]. Among them, the SUPG [1] (Streamline Upwind/Petrov–Galerkin) and later the GLS [2] (Galerkin least-squares) are two methods that stand out. Both methods add stabilization terms to the Galerkin formulation, providing good stability and accuracy to the numerical solution, when the exact
*
Corresponding author. E-mail addresses: [email protected] (E.G.D. do Carmo), [email protected] (G.B. Alvarez).
0045-7825/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.01.015
2384
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
solution of the problem is smooth [3]. However, the spurious oscillations remain when the exact solution possesses boundary layers. Many attempts have used the SUPG or GLS methods as a starting point for new stabilized formulations. Also, other attempts have added discontinuity capturing terms to the SUPG and GLS formulations [4,5]. In particular, the CAU method [6] (Consistent Approximate Upwind) maintains the term of SUPG and adds a nonlinear term providing an extra control concerning the function’s derivative in the direction of the approximate gradient, and eliminates the spurious oscillations. Nevertheless, when the exact solution of the problem is smooth, the CAU’s approximate solution presents an undesirable crosswind diffusion [7]. Building up a stable method that captures discontinuities and also preserves the accuracy of the SUPG or GLS methods for problems with smooth solution has been the greatest challenge. Recently, the authors developed the SAUPG method [14], which takes advantage of the best qualities of the SUPG and CAU methods, solving in part the difficulty placed above. We say ‘in part’, because SAUPG continues to be a nonlinear method, even when the physical problem is linear. In this work, based on the idea of the SAUPG method, we developed a new stabilized method using the GLS method as a starting point with a new upwind function. This new method introduces a new upwind function dependency on the geometry of the element and on the degree of the shape-functions’ polinomial. We studied the upwind function dependent on the geometry of the element (triangular and quadrilateral elements) and on the degree of the polynomial in the shape functions. It is well known that to solve the incompressible Navier–Stokes equations using a mixed formulation in velocities and pressures, the interpolation spaces of these two fields must satisfy the Babuska–Brezzi conditions [8,9]. The simplest elements that verify these conditions are quadratic elements, in terms of velocity, and the linear ones, in terms of pressure. The numerical results of GLS, CAU and the new method are shown to facilitate comparison. The GLS and CAU methods were chosen because they are representative of the two great difficulties that should be faced when solving this kind of problem: the lack of stability and the loss of accuracy. This paper is organized as follows. First, present the continuous problem and a review of the GLS and CAU methods. In the next section we develop the GLSAU method (Galerkin least-squares and approximate upwind). In the following section, several numerical examples for the three methods are shown, as a way to evaluate the new method’s performance. In Section 5, we finally present our conclusions.
2. The stationary scalar diffusion–advection equation 2.1. The boundary value problem Let X Rn be an open bounded domain, whose boundary C is a piecewise smooth boundary. The unit outward normal vector ^ n to C is defined almost everywhere. We shall consider the problem: r ðKr/Þ þ u r/ ¼ f /¼g
in X;
ð1Þ ð2Þ
on C;
together with the corresponding reduced problem with K ¼ 0, u r/ ¼ f
in X;
/ ¼ g on C ;
C ¼ fx 2 C : n u 6 0g;
ð1aÞ ð2aÞ
where / denotes the unknown quantity of the problem, K is the diffusivity tensor, u ¼ ðu1 ; . . . ; un Þ is the transport advective field, f is the volume source term and g is the boundary value prescribed for / on C or C .
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2385
2.2. Variational formulation Consider the set of all the kinematically admissible functions S and the space of the admissible variations V defined as: S ¼ fw 2 H 1 ðXÞ : w ¼ g on Cg, V ¼ fg 2 H 1 ðXÞ : g ¼ 0 on Cg, where H 1 ðXÞ is the standard Sobolev space. The variational problem associated to the boundary value problem, defined as (1) and (2), involves finding / 2 S that satisfies the variational equation: Z Z ½ðKr/Þ rg þ ðu r/Þg dX ¼ f g dX 8g 2 V : ð3Þ X
X
2.3. Finite element formulations S Consider M h ¼ fX1 ; . . . ; Xne g a finite element partition of X, such that: X ¼ ne e¼1 Xe , X ¼ X [ C, Xe ¼ Xe [ Ce and Xe \ Xe0 6¼ 0 if e 6¼ e0 . The respective finite element set and space of S and V are defined as: S h;k ¼ f/h 2 H 1 ðXÞ : /he 2 P k ðXe Þ and /h ¼ g on Cg, V h;k ¼ fgh 2 H 1 ðXÞ : ghe 2 P k ðXe Þ and gh ¼ 0 on Cg, where P k ðXe Þ is the space of degree k polynomials, gh denotes the interpolation of g and /he denotes the restriction of /h to Xe . The approach to the problem (3) consists of: finding /h 2 S h;k that satisfies 8gh 2 V h;k , AG ð/h ; gh Þ ¼ FG ðgh Þ
ð4Þ
Galerkin method ;
AG ð/h ; gh Þ þ ALS ð/h ; gh Þ ¼ FG ðgh Þ þ FLS ðgh Þ
ð5Þ
GLS method ;
AG ð/h ; gh Þ þ ASUPG ð/h ; gh Þ þ ACAU ð/h ; gh Þ ¼ FG ðgh Þ þ FSUPG ðgh Þ AG ¼
ne Z X
¼
½r ðKe r/he Þ þ ðue r/he Þpeh dXe ;
FLS ¼
Xe
s e he ½r ðKe rghe Þ þ ue rghe 8Xe 2 M h ; 2jue j
n X oni be;i ¼ ue;j ; ox j j¼1
ASUPG ¼
ne Z X
ne Z X Xe
e¼1
ACAU ¼
ne Z X e¼1
he jue j ; Pe ¼ 2K
Xe
f ghe dXe ;
ð6Þ ð7Þ
Xe
ne Z X e¼1
Xe
he ¼ 2
jue j ; jbe j
fe peh dXe ;
ð8Þ
ð9Þ
C1 se ¼ max 0; 1 ; C 2 Pe
½r ðKe r/he Þ þ ðue r/he Þ
Xe
e¼1
FSUPG ¼
ne Z X e¼1
ne Z X e¼1
peh
FG ¼
Xe
e¼1
ALS ¼
½ðKr/he Þ rghe þ ðu r/he Þghe dXe ;
CAU method ;
he s e ðue rghe ÞdXe ; 2jue j
ð10Þ
ð11Þ
he s e ðue rghe Þ dXe ; 2jue j
ð12Þ
Ce ð/he Þr/he rghe dXe ;
ð13Þ
fe
2386
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Ce ð/he Þ ¼
8 > > > 0; > <
if
" # h c c > s h s h jReð/ Þj jReð/he Þj > e e e e e > > ; : 2 sh h jue jjr/e j jr/he j e e
Reð/he Þ ¼ ue r/he r ðKe r/he Þ fe hce ¼ 2
jue vhe j ; jbce j
bce;i ¼
jReð/he Þj s c hc P e e; h s e he jue jjr/e j
jReð/he Þj s c hc if < e e; h jue jjr/e j se he
8Xe 2 M h ;
n X oni ðue;j vhe;j Þ; ox j j¼1
hce jue vhe j C1 ; sce ¼ max 0; 1 ; 2jKe j C2 Pec 8 > if jr/he j ¼ 0; < ue ; h Reð/e Þ vhe ¼ r/he ; if jr/he j 6¼ 0; > : ue 2 jr/he j Pec ¼
ð14Þ
ð15Þ ð16Þ
ð17Þ
ð18Þ
where ni ðXe Þði ¼ 1; . . . ; nÞ is the local coordinates system that map Xe in the usual standard elements Xe ; ue;j is the ‘‘j’’ component of ue ; and C1 and C2 are constants that depend of the element type, Xe and k (C1 ¼ C2 ¼ 1 for bilinear quadrilaterals elements). Remark 1. When the advective term is dominant, the Galerkin method can present spurious oscillations. In e 8Xe 2 M h , where he denote this case, the method is totally unstable, unless the mesh is very fined (he 6 2K jue j the diameter of Xe ). The instability of the Galerkin method is due to the lack of control on the gradient of the solution. The term ALS ð/h ; gh Þ introduces control in the streamline direction. Even so, when the problem possesses boundary layers the GLS method presents spurious oscillations. The approximate upwind direction Ueh ¼ ue vhe is in the direction of r/he . This allows the CAU method to possess greater stability, in regions near to internal and/or external boundary layers, than the GLS or SUPG methods.
3. The Galerkin least-squares and approximate upwind method (GLSAU) The fundamental idea that support the development of SAUPG method [14] is the divergence in the performance of both SUPG and CAU methods. The SUPG method lacks stability when applied to problems with boundary layers, while the CAU method loses accuracy in the solution for smooth problems. Here the same idea is used and we substituted the stabilization term of SUPG by the GLS term. Therefore, the new formulation involves finding /h 2 S h;k that satisfies 8gh 2 V h;k : AG ð/h ; gh Þ þ ALS ð/h ; gh Þ þ AAU ð/h ; gh Þ ¼ FG ðgh Þ þ FLS ðgh Þ; AAU ð/h ; ghe Þ ¼
ne Z X e¼1
Be ð/he Þr/he rghe dXe ;
ð19Þ ð20Þ
Xe
Be ð/he Þ ¼ B1 ½B2 c1 ½Ce ð/he Þc ¼ B1 ½B2 Ce ð/he Þc1 Ce ð/he Þ;
ð21Þ
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2387
where B1 , B2 and c are functions depending on the regularity of the approximate solution and were h
e Þj determined for the SAUPG method. In [14] the nondimensional variable ae ¼ jujReð/ , in order to evaluate jjr/h j e
e
the regularity of the approximate solution, and the functions B1 and B2 were defined as: B1 ðae Þ ¼ ½b1 c1 ðae Þ ¼
b1 ¼
ðc1Þðc1 þ1Þ
q2 ðae Þ 0
8 <1
: b2 ¼
b3 ¼ ½ he
q1 ðae Þ
;
;
B2 ¼
2 ; jue jhe se
ð22Þ
if ae < 1; if ae P 1;
ð23Þ if ae P 1;
b3 if b3 < ae 1 ½a þ b3 if b3 P ae 2 e he ¼ minf1; adimðhe Þg
ð24Þ
if ae < 1; 8Xe 2 M h ;
ð25Þ
where adimðÞ is a function that turns its argument dimensionless. Following the procedure described in [14], the function c is given by: 8 1 if ae P 1; > # <" P 1 1 if c 2 c¼ ð26Þ 1c2 if ae < 1; > : 1þc2 if c2 < 1 ½2 c2 8 1 1 > > if ae < ; < c4 2 4 ð27Þ c2 ¼ c3 ¼ c ; 1 1 > ½b2 3 if ae P ; > : þ ae 4 4 ( c4 ¼
½ae
q3 ðae Þ
½Reae
q3 ðae Þ
if ½Re
q3 ðae Þ
6 1;
if ½Re
q3 ðae Þ
> 1;
ð28Þ
where Re ¼ adimðjue jjr/he jÞ. At the same time, the functions q1 , q2 and q3 dependent on the parameter ae are obtained with the procedure described in [14]. 2
q1 ðae Þ ¼ q2 ðae Þ ¼ 1 ðae Þ ;
2
q3 ðae Þ ¼ 3 þ 12ae þ ðae Þ :
ð29Þ
Remark 2. When the solution is very smooth, c tends to 2 and the ALS ð/h ; gh Þ term prevails on AAU ð/h ; gh Þ in (20). In the case where the solution is not smooth, the contribution of the AAU ð/h ; gh Þ term is then significant. Remark 3. The GLSAU method was constructed to have a greater stability than the GLS method and a greater accuracy than the CAU method. When the solution of the problem is smooth, the accuracy of the GLSAU method is smaller or equal to the GLS method’s accuracy, but never greater. The GLS method’s solution is the limit, when referring to the GLSAU method’s accuracy. When the problem possesses boundary layers, the GLSAU method’s stability is smaller or equal to the CAU method’s stability. The CAU method’s solution is the limit, when referring to the GLSAU method’s stability [14]. This means that the GLSAU method degenerates into the GLS and CAU methods, when in extreme cases. This performance is
2388
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
possible with the help of he and Re numerical parameters, which reduce the upwind function depending on the regularity of the problem. Referring to Eqs. (25) and (28), the method is conservative, since it activates the CAU method at the slightest flaw indication in the GLS method’s solution. Remark 4. For the new method, the parameter se defined by (10), as well as, the sce given by (17) are redefined as: se ¼ so maxf0; 1 P1e g and sce ¼ so maxf0; 1 P1ec g, where Pe and Pec are respectively given by (10) and (17). The parameter so will be determined later on through the numerical experiments and was determined for bilinear quadrilateral elements in [14] as so ¼ 1.
4. Numerical results In this section we present the numerical results obtained from several standard numerical tests, which are divided into two types of problems: with internal and/or external boundary layers (Examples 1 and 2) and with smooth solutions (Example 3). In all cases, the medium is assumed homogeneous and isotropic with K ¼ 1010 . The considered domain is a square of unitary sides ð0; 1Þ ð0; 1Þ. In Fig. 1 we illustrate the two types of meshes used in the numerical tests. Five iterations were accomplished for the CAU and GLSAU methods. Example 1. Advection skew to the mesh. In this problem we have f ðx; yÞ ¼ 0 and the Dirichlet boundary conditions: /ðx; 0Þ ¼ 0;
/ðx; 1Þ ¼ 1;
/ð1; yÞ ¼ 0;
8y 2 ½0; 1;
/ð0; yÞ ¼ 0;
8y 2 ½0; 0:6;
/ð0; yÞ ¼ y 0:6;
8y 2 ½0:6; 0:65;
/ð0; yÞ ¼ 18ðy 0:65Þ þ 0:05;
8y 2 ½0:65; 0:70;
Fig. 1. The meshes used in the numerical tests.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
/ð0; yÞ ¼ ðy 0:70Þ þ 0:95; /ð0; yÞ ¼ 1;
2389
8y 2 ½0:70; 0:75;
8y 2 ½0:75; 1;
and three differents advective fields skew to the mesh were considered. Example 2. External boundary layers and semicircular internal layer. In this example, we have f ðx; yÞ ¼ 0, where the Dirichlet boundary conditions are the same as in the latter example and the advective field is u ¼ ðy; xÞ. Example 3. Advection of a sine hill in a rotating flow field. The problem statement is shown in Fig. 2. The flow’s rotation is determined by the advective field u ¼ ðy; xÞ. Along the external boundary / ¼ 0, and on the internal ‘boundary’ OA a sine hill function /ð0; yÞ ¼ sinð2pyÞ 8y 2 ½0:5; 0 is prescribed. The objective of following subsections is to verify, numerically, if the functions qi ðae Þ (i ¼ 1; 2; 3) (29), determined through the numerical experiments for bilinear quadrilateral elements in [14], continue being valid for linear or quadratic triangular elements, as well as, for quadratic quadrilateral elements. The results following confirm that the functions qi ðae Þ do not depend on the geometry of the element nor the degree of the polynomial in the shape functions. Also, it was verified in all the examples that the three methods introduce diffusivity in excess for the external boundary layers, when the so defined for bilinear quadrilateral elements is used. The numerical experiments led us to redefine so 8Xe 2 M h , as follows:
1 if Ce \ Cþ ¼ ; seo if Ce \ Cþ 6¼ ; ; Cþ ¼ fx 2 C : ^ n u > 0g
so ¼
Fig. 2. Advection in a rotating flow field: problem statement.
ð30Þ
2390
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
At present, not a perfect theoretical support has been obtained and not much has been studied about the parameter se dependence on the element Peclet number for the higher order quadratic rectangular or triangular elements. Some few attempts can be found in [10–13]. Many of them are based on analyses carried out for the one-dimensional case. In these attempts, the parameter is optimally determined, in the sense that it leads to nodally exact solutions, as they also try to extend the analysis for the two-dimensional case. Other attempts are based on heuristic, intuitive and experimental criterion. We propose (30), based on a numerical analyses of the stability and on the fact that the internal and external boundary layers are of different physical nature [3]. These are our arguments to justify the discontinuity of so in (30). We have defined the criterion of ‘optimality’ in the sense that diffusive or oscillatory results are not obtained, i.e., that the upwind function is not overestimated or underestimated. All figures coming afterward show the numerical results for each example. In each figure, cases (a) correspond to each one of the three methods, using so ¼ seo for the reduced problem [1a–2a]. Cases (b) refer to the solutions of the methods obtained with so given by (30). Cases (c) were obtained with so ¼ 1, that is to say, the upwind functions corresponding to the value determined for the bilinear quadrilateral elements. 4.1. Approximate solutions using quadratic quadrilateral elements (9-noded) For Examples 1 and 2, we can observe in Figs. 3–6 that the approximate solutions of the GLSAU and CAU methods are very close. In the cases (a) overshoots are obtained for the internal sharp layer, whereas in the cases (c) undershoots occur for the external boundary layers. For the Example 3, the GLS and GLSAU methods yield to very good results (see Fig. 7). 4.2. Approximate solution using linear triangular elements In this subsection, the previous subsection’s comments are still valid. In the cases (b) of the Figs. 8–12 we obtained the best numerical solution for each method. It should be highlighted, that the internal sharp layer of Example 1 (see Fig. 9) is very well approached by the GLS, CAU and GLSAU solutions. This is something expected, since in this case the advective field u is aligned to the mesh. Fig. 1b illustrates the mesh of triangular elements used in this subsection. Similar results are obtained for quadratic triangular elements, so they are not shown here. Here, again we should highlight the overshoots obtained for the internal sharp layer in the cases (a). 4.3. Some comments about the convergence of CAU and GLSAU methods The CAU and GLSAU methods are nonlinear methods because they need an approximate solution of the problem to estimate the residual of Eq. (1). Thus, each method generates a sequence of approximate solutions that we can denote as f/CAU g and f/GLSAU g. n n The results presented in this subsection were obtained using three meshes with bilinear quadrilateral elements (see Fig. 1a). The distance between the elements of each sequence was calculated using the Euclidean norm. Tables 1 and 2 show the convergence results, requiring that j/in1 /in j 6 e, where i ¼ fcau; glsaug. We can observe in these tables that both, CAU and GLSAU sequences are Cauchy sequences. We can also observe that, in a general, the sequence generated by GLSAU possesses a faster convergence than the sequence of the CAU method. The Example 3 deserves an outstanding place. It is well known that, for this case, the exact solution is smooth and, therefore, the GLS method produces a very good approximate solution. Thus, the residual of the Eq. (1) for the GLS or SUPG solutions is small. In spite of this, we can observe in Tables 1 and 2 that,
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2391
Fig. 3. Example 1, u ¼ ð2; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
while the CAU method needs to increase the number of iterations when the mesh is refined, the GLSAU method does not need to.
2392
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 4. Example 1, u ¼ ð1; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 5. Example 1, u ¼ ð1; 2Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2393
2394
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 6. Example 2, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 7. Example 3, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2395
2396
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 8. Example 1, u ¼ ð2; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 9. Example 1, u ¼ ð1; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2397
2398
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 10. Example 1, u ¼ ð1; 2Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 11. Example 2, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2399
2400
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 12. Example 3, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2401
Table 1 Iteration number of the GLS and GLSAU methods for e ¼ 102 Examples
n-Iterations he ¼ 0:1
Example Example Example Example Example
1, 1, 1, 2, 3
u ¼ ð2; 1Þ u ¼ ð1; 1Þ u ¼ ð1; 2Þ u ¼ ðy; xÞ
he ¼ 0:05
he ¼ 0:025
CAU
GLSAU
CAU
GLSAU
CAU
GLSAU
4 4 4 5 4
4 3 4 4 3
6 4 5 5 4
5 4 5 4 2
10 5 7 10 7
12 5 7 5 2
Table 2 Iteration number of the GLS and GLSAU methods for e ¼ 103 Examples
n-Iterations he ¼ 0:1
Example Example Example Example Example
1, 1, 1, 2, 3
u ¼ ð2; 1Þ u ¼ ð1; 1Þ u ¼ ð1; 2Þ u ¼ ðy; xÞ
he ¼ 0:05
he ¼ 0:025
CAU
GLSAU
CAU
GLSAU
CAU
GLSAU
6 5 5 7 5
5 5 5 6 3
10 5 7 10 7
9 6 7 5 3
25 8 13 26 16
23 9 13 7 3
5. Conclusions In this paper, we have studied and compared the GLS, CAU and a new stabilized finite element formulation (GLSAU) to solve scalar convection–diffusion problems. The new formulation preserves the best qualities of the GLS and CAU methods, in terms of stability and accuracy. The numerical experiments showed that the functions q1 ðae Þ, q2 ðae Þ and q3 ðae Þ do not depend on the geometry of the element nor on the degree of the polynomial in the shape functions. We also studied the upwind function dependence on the geometry of the element and the degree of the polynomial in the shape functions. The numerical results led us to define so as in (30), i.e., value 1 is necessary for the internal sharp layer, while seo is the optimum value to avoid overshoots and undershoots for the external boundary layer. This modification on stabilizing parameter near the outflow boundaries increases the quality of the numerical solution for the GLS, CAU and GLSAU methods. Although we lack a theoretical error estimate, the simplicity of the new scheme, as well as the parameter s must be considered for practical purposes. Acknowledgements The authors wish to thank the Brazilian research funding agencies CNPq and FAPERJ for its support to this work. References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259.
2402
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
[2] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VII. The Galerkin-least-squares method for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173–189. [3] C. Johnson, U. N€avert, J. Pitkaranta, Finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 285–312. [4] T.J.R. Hughes, M. Mallet, A. Mizukami, A new finite element formulation for computational fluid dynamics: II Beyond SUPG, Comput. Methods Appl. Mech. Engrg. 54 (1986) 341–355. [5] T.J.R. Hughes, M. Mallet, A new finite element formulation for computational fluid dynamics: IV A discontinuity capturing operator for multidimensional advective–diffusive systems, Comput. Methods Appl. Mech. Engrg. 58 (1986) 329–336. [6] A.C. Gale~ ao, E.G. Do Carmo, A consistent approximate upwind Petrov–Galerkin method for convection-dominated problems, Comput. Methods Appl. Mech. Engrg. 68 (1988) 83–95. [7] E.G. Do Carmo, A.C. Gale~ao, Feedback Petrov-Galerkin methods for convection-dominated problems, Comput. Methods Appl. Mech. Engrg. 88 (1991) 1–16. [8] I. Babuska, Error bounds for finite element method, Numer. Math. 16 (1971) 322–333. [9] F. Brezzi, On the existence, uniqueness and approximation of saddle point problems arising from Lagrange multipliers, RAIRO, Ser. Rouge Anal. Numer. 8 (1976) R-2. [10] I. Christie, A.R. Mitchell, Upwinding of high order Galerkin methods in conduction–convection problems, Int. J. Numer. Methods Engrg. 14 (1978) 1764–1771. [11] J.C. Heinrich, On quadratic elements in finite element solution of steady-state convection–diffusion equation, Int. J. Numer. Methods Engrg. 15 (1980) 1041–1052. [12] A. Mizukami, An implementation of the streamline-upwind/Petrov–Galerkin method for linear triangular elements, Comput. Methods Appl. Mech. Engrg. 49 (1985) 357–364. [13] R. Codina, E. O~ nate, M. Cervera, The intrinsic time for the streamline upwind Petrov–Galerkin formulation using quadratic elements, Comput. Methods Appl. Mech. Engrg. 94 (1992) 239–262. [14] E.G.D. do Carmo, G.B. Alvarez, A new stabilized finite element formulation for scalar convection–diffusion problems: the sreamline and approximate upwind/Petrov–Galerkin method, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3379–3396.
A new upwind function in stabilized finite element formulations, using linear and quadratic elements for scalar convection–diffusion problems Eduardo Gomes Dutra do Carmo
a,*
, Gustavo Benitez Alvarez
b
a
b
Department of Nuclear Engineering, COPPE––Federal University of Rio de Janeiro, Ilha do Fund~ao, P.B. 68509, Rio de Janeiro, RJ 21945-970, Brazil LNCC––National Laboratory of Scientific Computation, Av. Getulio Vargas 333, P.B. 95113, Petropolis, RJ 25651-070, Brazil Received 15 April 2003; received in revised form 7 November 2003; accepted 7 January 2004
Abstract This article presents a study on the function of the Peclet number that appears in the upwind function of the GLS, CAU and SAUPG methods, using linear and quadratic triangular elements, as well as bilinear and quadratic quadrilateral elements. The numerical experiments indicate that in the case of elements with faces on the ‘‘outflow boundary’’, the upwind function is strongly dependent on the geometry and on the degree of the interpolate polynomial. A new function of the Peclet number is therefore proposed, to take these effects into consideration, creating a new method of stabilization for higher order elements. This new stabilized and accurate finite element formulation for convection dominated problems, based on the idea of the SAUPG method, uses the GLS formulation as starting point. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Stabilized FEM; Convection–diffusion; Stabilization
1. Introduction In general, the solution of diffusive-convective problems possesses boundary layers, that are small subregions where derivatives of the solution are very large. In this case, the classic finite element method or Galerkin method is inappropriate, because its numerical solution presents spurious oscillations. During the last two decades, various attempts to obtain new methods have been carried out to avoid such oscillations. These new methods are nowadays known as stabilized finite element methods. References about some of them can be found in [14]. Among them, the SUPG [1] (Streamline Upwind/Petrov–Galerkin) and later the GLS [2] (Galerkin least-squares) are two methods that stand out. Both methods add stabilization terms to the Galerkin formulation, providing good stability and accuracy to the numerical solution, when the exact
*
Corresponding author. E-mail addresses: [email protected] (E.G.D. do Carmo), [email protected] (G.B. Alvarez).
0045-7825/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.01.015
2384
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
solution of the problem is smooth [3]. However, the spurious oscillations remain when the exact solution possesses boundary layers. Many attempts have used the SUPG or GLS methods as a starting point for new stabilized formulations. Also, other attempts have added discontinuity capturing terms to the SUPG and GLS formulations [4,5]. In particular, the CAU method [6] (Consistent Approximate Upwind) maintains the term of SUPG and adds a nonlinear term providing an extra control concerning the function’s derivative in the direction of the approximate gradient, and eliminates the spurious oscillations. Nevertheless, when the exact solution of the problem is smooth, the CAU’s approximate solution presents an undesirable crosswind diffusion [7]. Building up a stable method that captures discontinuities and also preserves the accuracy of the SUPG or GLS methods for problems with smooth solution has been the greatest challenge. Recently, the authors developed the SAUPG method [14], which takes advantage of the best qualities of the SUPG and CAU methods, solving in part the difficulty placed above. We say ‘in part’, because SAUPG continues to be a nonlinear method, even when the physical problem is linear. In this work, based on the idea of the SAUPG method, we developed a new stabilized method using the GLS method as a starting point with a new upwind function. This new method introduces a new upwind function dependency on the geometry of the element and on the degree of the shape-functions’ polinomial. We studied the upwind function dependent on the geometry of the element (triangular and quadrilateral elements) and on the degree of the polynomial in the shape functions. It is well known that to solve the incompressible Navier–Stokes equations using a mixed formulation in velocities and pressures, the interpolation spaces of these two fields must satisfy the Babuska–Brezzi conditions [8,9]. The simplest elements that verify these conditions are quadratic elements, in terms of velocity, and the linear ones, in terms of pressure. The numerical results of GLS, CAU and the new method are shown to facilitate comparison. The GLS and CAU methods were chosen because they are representative of the two great difficulties that should be faced when solving this kind of problem: the lack of stability and the loss of accuracy. This paper is organized as follows. First, present the continuous problem and a review of the GLS and CAU methods. In the next section we develop the GLSAU method (Galerkin least-squares and approximate upwind). In the following section, several numerical examples for the three methods are shown, as a way to evaluate the new method’s performance. In Section 5, we finally present our conclusions.
2. The stationary scalar diffusion–advection equation 2.1. The boundary value problem Let X Rn be an open bounded domain, whose boundary C is a piecewise smooth boundary. The unit outward normal vector ^ n to C is defined almost everywhere. We shall consider the problem: r ðKr/Þ þ u r/ ¼ f /¼g
in X;
ð1Þ ð2Þ
on C;
together with the corresponding reduced problem with K ¼ 0, u r/ ¼ f
in X;
/ ¼ g on C ;
C ¼ fx 2 C : n u 6 0g;
ð1aÞ ð2aÞ
where / denotes the unknown quantity of the problem, K is the diffusivity tensor, u ¼ ðu1 ; . . . ; un Þ is the transport advective field, f is the volume source term and g is the boundary value prescribed for / on C or C .
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2385
2.2. Variational formulation Consider the set of all the kinematically admissible functions S and the space of the admissible variations V defined as: S ¼ fw 2 H 1 ðXÞ : w ¼ g on Cg, V ¼ fg 2 H 1 ðXÞ : g ¼ 0 on Cg, where H 1 ðXÞ is the standard Sobolev space. The variational problem associated to the boundary value problem, defined as (1) and (2), involves finding / 2 S that satisfies the variational equation: Z Z ½ðKr/Þ rg þ ðu r/Þg dX ¼ f g dX 8g 2 V : ð3Þ X
X
2.3. Finite element formulations S Consider M h ¼ fX1 ; . . . ; Xne g a finite element partition of X, such that: X ¼ ne e¼1 Xe , X ¼ X [ C, Xe ¼ Xe [ Ce and Xe \ Xe0 6¼ 0 if e 6¼ e0 . The respective finite element set and space of S and V are defined as: S h;k ¼ f/h 2 H 1 ðXÞ : /he 2 P k ðXe Þ and /h ¼ g on Cg, V h;k ¼ fgh 2 H 1 ðXÞ : ghe 2 P k ðXe Þ and gh ¼ 0 on Cg, where P k ðXe Þ is the space of degree k polynomials, gh denotes the interpolation of g and /he denotes the restriction of /h to Xe . The approach to the problem (3) consists of: finding /h 2 S h;k that satisfies 8gh 2 V h;k , AG ð/h ; gh Þ ¼ FG ðgh Þ
ð4Þ
Galerkin method ;
AG ð/h ; gh Þ þ ALS ð/h ; gh Þ ¼ FG ðgh Þ þ FLS ðgh Þ
ð5Þ
GLS method ;
AG ð/h ; gh Þ þ ASUPG ð/h ; gh Þ þ ACAU ð/h ; gh Þ ¼ FG ðgh Þ þ FSUPG ðgh Þ AG ¼
ne Z X
¼
½r ðKe r/he Þ þ ðue r/he Þpeh dXe ;
FLS ¼
Xe
s e he ½r ðKe rghe Þ þ ue rghe 8Xe 2 M h ; 2jue j
n X oni be;i ¼ ue;j ; ox j j¼1
ASUPG ¼
ne Z X
ne Z X Xe
e¼1
ACAU ¼
ne Z X e¼1
he jue j ; Pe ¼ 2K
Xe
f ghe dXe ;
ð6Þ ð7Þ
Xe
ne Z X e¼1
Xe
he ¼ 2
jue j ; jbe j
fe peh dXe ;
ð8Þ
ð9Þ
C1 se ¼ max 0; 1 ; C 2 Pe
½r ðKe r/he Þ þ ðue r/he Þ
Xe
e¼1
FSUPG ¼
ne Z X e¼1
ne Z X e¼1
peh
FG ¼
Xe
e¼1
ALS ¼
½ðKr/he Þ rghe þ ðu r/he Þghe dXe ;
CAU method ;
he s e ðue rghe ÞdXe ; 2jue j
ð10Þ
ð11Þ
he s e ðue rghe Þ dXe ; 2jue j
ð12Þ
Ce ð/he Þr/he rghe dXe ;
ð13Þ
fe
2386
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Ce ð/he Þ ¼
8 > > > 0; > <
if
" # h c c > s h s h jReð/ Þj jReð/he Þj > e e e e e > > ; : 2 sh h jue jjr/e j jr/he j e e
Reð/he Þ ¼ ue r/he r ðKe r/he Þ fe hce ¼ 2
jue vhe j ; jbce j
bce;i ¼
jReð/he Þj s c hc P e e; h s e he jue jjr/e j
jReð/he Þj s c hc if < e e; h jue jjr/e j se he
8Xe 2 M h ;
n X oni ðue;j vhe;j Þ; ox j j¼1
hce jue vhe j C1 ; sce ¼ max 0; 1 ; 2jKe j C2 Pec 8 > if jr/he j ¼ 0; < ue ; h Reð/e Þ vhe ¼ r/he ; if jr/he j 6¼ 0; > : ue 2 jr/he j Pec ¼
ð14Þ
ð15Þ ð16Þ
ð17Þ
ð18Þ
where ni ðXe Þði ¼ 1; . . . ; nÞ is the local coordinates system that map Xe in the usual standard elements Xe ; ue;j is the ‘‘j’’ component of ue ; and C1 and C2 are constants that depend of the element type, Xe and k (C1 ¼ C2 ¼ 1 for bilinear quadrilaterals elements). Remark 1. When the advective term is dominant, the Galerkin method can present spurious oscillations. In e 8Xe 2 M h , where he denote this case, the method is totally unstable, unless the mesh is very fined (he 6 2K jue j the diameter of Xe ). The instability of the Galerkin method is due to the lack of control on the gradient of the solution. The term ALS ð/h ; gh Þ introduces control in the streamline direction. Even so, when the problem possesses boundary layers the GLS method presents spurious oscillations. The approximate upwind direction Ueh ¼ ue vhe is in the direction of r/he . This allows the CAU method to possess greater stability, in regions near to internal and/or external boundary layers, than the GLS or SUPG methods.
3. The Galerkin least-squares and approximate upwind method (GLSAU) The fundamental idea that support the development of SAUPG method [14] is the divergence in the performance of both SUPG and CAU methods. The SUPG method lacks stability when applied to problems with boundary layers, while the CAU method loses accuracy in the solution for smooth problems. Here the same idea is used and we substituted the stabilization term of SUPG by the GLS term. Therefore, the new formulation involves finding /h 2 S h;k that satisfies 8gh 2 V h;k : AG ð/h ; gh Þ þ ALS ð/h ; gh Þ þ AAU ð/h ; gh Þ ¼ FG ðgh Þ þ FLS ðgh Þ; AAU ð/h ; ghe Þ ¼
ne Z X e¼1
Be ð/he Þr/he rghe dXe ;
ð19Þ ð20Þ
Xe
Be ð/he Þ ¼ B1 ½B2 c1 ½Ce ð/he Þc ¼ B1 ½B2 Ce ð/he Þc1 Ce ð/he Þ;
ð21Þ
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2387
where B1 , B2 and c are functions depending on the regularity of the approximate solution and were h
e Þj determined for the SAUPG method. In [14] the nondimensional variable ae ¼ jujReð/ , in order to evaluate jjr/h j e
e
the regularity of the approximate solution, and the functions B1 and B2 were defined as: B1 ðae Þ ¼ ½b1 c1 ðae Þ ¼
b1 ¼
ðc1Þðc1 þ1Þ
q2 ðae Þ 0
8 <1
: b2 ¼
b3 ¼ ½ he
q1 ðae Þ
;
;
B2 ¼
2 ; jue jhe se
ð22Þ
if ae < 1; if ae P 1;
ð23Þ if ae P 1;
b3 if b3 < ae 1 ½a þ b3 if b3 P ae 2 e he ¼ minf1; adimðhe Þg
ð24Þ
if ae < 1; 8Xe 2 M h ;
ð25Þ
where adimðÞ is a function that turns its argument dimensionless. Following the procedure described in [14], the function c is given by: 8 1 if ae P 1; > # <" P 1 1 if c 2 c¼ ð26Þ 1c2 if ae < 1; > : 1þc2 if c2 < 1 ½2 c2 8 1 1 > > if ae < ; < c4 2 4 ð27Þ c2 ¼ c3 ¼ c ; 1 1 > ½b2 3 if ae P ; > : þ ae 4 4 ( c4 ¼
½ae
q3 ðae Þ
½Reae
q3 ðae Þ
if ½Re
q3 ðae Þ
6 1;
if ½Re
q3 ðae Þ
> 1;
ð28Þ
where Re ¼ adimðjue jjr/he jÞ. At the same time, the functions q1 , q2 and q3 dependent on the parameter ae are obtained with the procedure described in [14]. 2
q1 ðae Þ ¼ q2 ðae Þ ¼ 1 ðae Þ ;
2
q3 ðae Þ ¼ 3 þ 12ae þ ðae Þ :
ð29Þ
Remark 2. When the solution is very smooth, c tends to 2 and the ALS ð/h ; gh Þ term prevails on AAU ð/h ; gh Þ in (20). In the case where the solution is not smooth, the contribution of the AAU ð/h ; gh Þ term is then significant. Remark 3. The GLSAU method was constructed to have a greater stability than the GLS method and a greater accuracy than the CAU method. When the solution of the problem is smooth, the accuracy of the GLSAU method is smaller or equal to the GLS method’s accuracy, but never greater. The GLS method’s solution is the limit, when referring to the GLSAU method’s accuracy. When the problem possesses boundary layers, the GLSAU method’s stability is smaller or equal to the CAU method’s stability. The CAU method’s solution is the limit, when referring to the GLSAU method’s stability [14]. This means that the GLSAU method degenerates into the GLS and CAU methods, when in extreme cases. This performance is
2388
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
possible with the help of he and Re numerical parameters, which reduce the upwind function depending on the regularity of the problem. Referring to Eqs. (25) and (28), the method is conservative, since it activates the CAU method at the slightest flaw indication in the GLS method’s solution. Remark 4. For the new method, the parameter se defined by (10), as well as, the sce given by (17) are redefined as: se ¼ so maxf0; 1 P1e g and sce ¼ so maxf0; 1 P1ec g, where Pe and Pec are respectively given by (10) and (17). The parameter so will be determined later on through the numerical experiments and was determined for bilinear quadrilateral elements in [14] as so ¼ 1.
4. Numerical results In this section we present the numerical results obtained from several standard numerical tests, which are divided into two types of problems: with internal and/or external boundary layers (Examples 1 and 2) and with smooth solutions (Example 3). In all cases, the medium is assumed homogeneous and isotropic with K ¼ 1010 . The considered domain is a square of unitary sides ð0; 1Þ ð0; 1Þ. In Fig. 1 we illustrate the two types of meshes used in the numerical tests. Five iterations were accomplished for the CAU and GLSAU methods. Example 1. Advection skew to the mesh. In this problem we have f ðx; yÞ ¼ 0 and the Dirichlet boundary conditions: /ðx; 0Þ ¼ 0;
/ðx; 1Þ ¼ 1;
/ð1; yÞ ¼ 0;
8y 2 ½0; 1;
/ð0; yÞ ¼ 0;
8y 2 ½0; 0:6;
/ð0; yÞ ¼ y 0:6;
8y 2 ½0:6; 0:65;
/ð0; yÞ ¼ 18ðy 0:65Þ þ 0:05;
8y 2 ½0:65; 0:70;
Fig. 1. The meshes used in the numerical tests.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
/ð0; yÞ ¼ ðy 0:70Þ þ 0:95; /ð0; yÞ ¼ 1;
2389
8y 2 ½0:70; 0:75;
8y 2 ½0:75; 1;
and three differents advective fields skew to the mesh were considered. Example 2. External boundary layers and semicircular internal layer. In this example, we have f ðx; yÞ ¼ 0, where the Dirichlet boundary conditions are the same as in the latter example and the advective field is u ¼ ðy; xÞ. Example 3. Advection of a sine hill in a rotating flow field. The problem statement is shown in Fig. 2. The flow’s rotation is determined by the advective field u ¼ ðy; xÞ. Along the external boundary / ¼ 0, and on the internal ‘boundary’ OA a sine hill function /ð0; yÞ ¼ sinð2pyÞ 8y 2 ½0:5; 0 is prescribed. The objective of following subsections is to verify, numerically, if the functions qi ðae Þ (i ¼ 1; 2; 3) (29), determined through the numerical experiments for bilinear quadrilateral elements in [14], continue being valid for linear or quadratic triangular elements, as well as, for quadratic quadrilateral elements. The results following confirm that the functions qi ðae Þ do not depend on the geometry of the element nor the degree of the polynomial in the shape functions. Also, it was verified in all the examples that the three methods introduce diffusivity in excess for the external boundary layers, when the so defined for bilinear quadrilateral elements is used. The numerical experiments led us to redefine so 8Xe 2 M h , as follows:
1 if Ce \ Cþ ¼ ; seo if Ce \ Cþ 6¼ ; ; Cþ ¼ fx 2 C : ^ n u > 0g
so ¼
Fig. 2. Advection in a rotating flow field: problem statement.
ð30Þ
2390
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
At present, not a perfect theoretical support has been obtained and not much has been studied about the parameter se dependence on the element Peclet number for the higher order quadratic rectangular or triangular elements. Some few attempts can be found in [10–13]. Many of them are based on analyses carried out for the one-dimensional case. In these attempts, the parameter is optimally determined, in the sense that it leads to nodally exact solutions, as they also try to extend the analysis for the two-dimensional case. Other attempts are based on heuristic, intuitive and experimental criterion. We propose (30), based on a numerical analyses of the stability and on the fact that the internal and external boundary layers are of different physical nature [3]. These are our arguments to justify the discontinuity of so in (30). We have defined the criterion of ‘optimality’ in the sense that diffusive or oscillatory results are not obtained, i.e., that the upwind function is not overestimated or underestimated. All figures coming afterward show the numerical results for each example. In each figure, cases (a) correspond to each one of the three methods, using so ¼ seo for the reduced problem [1a–2a]. Cases (b) refer to the solutions of the methods obtained with so given by (30). Cases (c) were obtained with so ¼ 1, that is to say, the upwind functions corresponding to the value determined for the bilinear quadrilateral elements. 4.1. Approximate solutions using quadratic quadrilateral elements (9-noded) For Examples 1 and 2, we can observe in Figs. 3–6 that the approximate solutions of the GLSAU and CAU methods are very close. In the cases (a) overshoots are obtained for the internal sharp layer, whereas in the cases (c) undershoots occur for the external boundary layers. For the Example 3, the GLS and GLSAU methods yield to very good results (see Fig. 7). 4.2. Approximate solution using linear triangular elements In this subsection, the previous subsection’s comments are still valid. In the cases (b) of the Figs. 8–12 we obtained the best numerical solution for each method. It should be highlighted, that the internal sharp layer of Example 1 (see Fig. 9) is very well approached by the GLS, CAU and GLSAU solutions. This is something expected, since in this case the advective field u is aligned to the mesh. Fig. 1b illustrates the mesh of triangular elements used in this subsection. Similar results are obtained for quadratic triangular elements, so they are not shown here. Here, again we should highlight the overshoots obtained for the internal sharp layer in the cases (a). 4.3. Some comments about the convergence of CAU and GLSAU methods The CAU and GLSAU methods are nonlinear methods because they need an approximate solution of the problem to estimate the residual of Eq. (1). Thus, each method generates a sequence of approximate solutions that we can denote as f/CAU g and f/GLSAU g. n n The results presented in this subsection were obtained using three meshes with bilinear quadrilateral elements (see Fig. 1a). The distance between the elements of each sequence was calculated using the Euclidean norm. Tables 1 and 2 show the convergence results, requiring that j/in1 /in j 6 e, where i ¼ fcau; glsaug. We can observe in these tables that both, CAU and GLSAU sequences are Cauchy sequences. We can also observe that, in a general, the sequence generated by GLSAU possesses a faster convergence than the sequence of the CAU method. The Example 3 deserves an outstanding place. It is well known that, for this case, the exact solution is smooth and, therefore, the GLS method produces a very good approximate solution. Thus, the residual of the Eq. (1) for the GLS or SUPG solutions is small. In spite of this, we can observe in Tables 1 and 2 that,
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2391
Fig. 3. Example 1, u ¼ ð2; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
while the CAU method needs to increase the number of iterations when the mesh is refined, the GLSAU method does not need to.
2392
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 4. Example 1, u ¼ ð1; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 5. Example 1, u ¼ ð1; 2Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2393
2394
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 6. Example 2, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 7. Example 3, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2395
2396
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 8. Example 1, u ¼ ð2; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 9. Example 1, u ¼ ð1; 1Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2397
2398
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 10. Example 1, u ¼ ð1; 2Þ, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 11. Example 2, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
2399
2400
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
Fig. 12. Example 3, (a) so ¼ seo ¼ 12, (b) seo ¼ 12, (c) so ¼ 1.
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
2401
Table 1 Iteration number of the GLS and GLSAU methods for e ¼ 102 Examples
n-Iterations he ¼ 0:1
Example Example Example Example Example
1, 1, 1, 2, 3
u ¼ ð2; 1Þ u ¼ ð1; 1Þ u ¼ ð1; 2Þ u ¼ ðy; xÞ
he ¼ 0:05
he ¼ 0:025
CAU
GLSAU
CAU
GLSAU
CAU
GLSAU
4 4 4 5 4
4 3 4 4 3
6 4 5 5 4
5 4 5 4 2
10 5 7 10 7
12 5 7 5 2
Table 2 Iteration number of the GLS and GLSAU methods for e ¼ 103 Examples
n-Iterations he ¼ 0:1
Example Example Example Example Example
1, 1, 1, 2, 3
u ¼ ð2; 1Þ u ¼ ð1; 1Þ u ¼ ð1; 2Þ u ¼ ðy; xÞ
he ¼ 0:05
he ¼ 0:025
CAU
GLSAU
CAU
GLSAU
CAU
GLSAU
6 5 5 7 5
5 5 5 6 3
10 5 7 10 7
9 6 7 5 3
25 8 13 26 16
23 9 13 7 3
5. Conclusions In this paper, we have studied and compared the GLS, CAU and a new stabilized finite element formulation (GLSAU) to solve scalar convection–diffusion problems. The new formulation preserves the best qualities of the GLS and CAU methods, in terms of stability and accuracy. The numerical experiments showed that the functions q1 ðae Þ, q2 ðae Þ and q3 ðae Þ do not depend on the geometry of the element nor on the degree of the polynomial in the shape functions. We also studied the upwind function dependence on the geometry of the element and the degree of the polynomial in the shape functions. The numerical results led us to define so as in (30), i.e., value 1 is necessary for the internal sharp layer, while seo is the optimum value to avoid overshoots and undershoots for the external boundary layer. This modification on stabilizing parameter near the outflow boundaries increases the quality of the numerical solution for the GLS, CAU and GLSAU methods. Although we lack a theoretical error estimate, the simplicity of the new scheme, as well as the parameter s must be considered for practical purposes. Acknowledgements The authors wish to thank the Brazilian research funding agencies CNPq and FAPERJ for its support to this work. References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259.
2402
E.G.D. do Carmo, G.B. Alvarez / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2383–2402
[2] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VII. The Galerkin-least-squares method for advective–diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (1989) 173–189. [3] C. Johnson, U. N€avert, J. Pitkaranta, Finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 285–312. [4] T.J.R. Hughes, M. Mallet, A. Mizukami, A new finite element formulation for computational fluid dynamics: II Beyond SUPG, Comput. Methods Appl. Mech. Engrg. 54 (1986) 341–355. [5] T.J.R. Hughes, M. Mallet, A new finite element formulation for computational fluid dynamics: IV A discontinuity capturing operator for multidimensional advective–diffusive systems, Comput. Methods Appl. Mech. Engrg. 58 (1986) 329–336. [6] A.C. Gale~ ao, E.G. Do Carmo, A consistent approximate upwind Petrov–Galerkin method for convection-dominated problems, Comput. Methods Appl. Mech. Engrg. 68 (1988) 83–95. [7] E.G. Do Carmo, A.C. Gale~ao, Feedback Petrov-Galerkin methods for convection-dominated problems, Comput. Methods Appl. Mech. Engrg. 88 (1991) 1–16. [8] I. Babuska, Error bounds for finite element method, Numer. Math. 16 (1971) 322–333. [9] F. Brezzi, On the existence, uniqueness and approximation of saddle point problems arising from Lagrange multipliers, RAIRO, Ser. Rouge Anal. Numer. 8 (1976) R-2. [10] I. Christie, A.R. Mitchell, Upwinding of high order Galerkin methods in conduction–convection problems, Int. J. Numer. Methods Engrg. 14 (1978) 1764–1771. [11] J.C. Heinrich, On quadratic elements in finite element solution of steady-state convection–diffusion equation, Int. J. Numer. Methods Engrg. 15 (1980) 1041–1052. [12] A. Mizukami, An implementation of the streamline-upwind/Petrov–Galerkin method for linear triangular elements, Comput. Methods Appl. Mech. Engrg. 49 (1985) 357–364. [13] R. Codina, E. O~ nate, M. Cervera, The intrinsic time for the streamline upwind Petrov–Galerkin formulation using quadratic elements, Comput. Methods Appl. Mech. Engrg. 94 (1992) 239–262. [14] E.G.D. do Carmo, G.B. Alvarez, A new stabilized finite element formulation for scalar convection–diffusion problems: the sreamline and approximate upwind/Petrov–Galerkin method, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3379–3396.