- Email: [email protected]
Stabilized solution of the multidimensional advection–diffusion–absorption equation using linear finite elements
Computers & Fluids 36 (2007) 92–112 www.elsevier.com/locate/compfluid
Stabilized solution of the multidimensional advection–diffusion–absorption equation using linear finite elements Eugenio On˜ate *, Juan Miquel, Francisco Za´rate International Center for Numerical Methods in Engineering (CIMNE), Universidad Polite´cnica de Catalun˜a, Edificio C1, Gran Capita´n s/n, 08034 Barcelona, Spain Received 23 June 2005; accepted 21 July 2005 Available online 1 December 2005
Abstract A stabilized finite element method (FEM) for the multidimensional steady state advection–diffusion–absorption equation is presented. The stabilized formulation is based on the modified governing differential equations derived via the finite calculus (FIC) method. For 1D problems the stabilization terms act as a nonlinear additional diffusion governed by a single stabilization parameter. It is shown that for multidimensional problems an orthotropic stabilizing diffusion must be added along the principal directions of curvature of the solution. A simple iterative algorithm yielding a stable and accurate solution for all the range of physical parameters and boundary conditions is described. Numerical results for 1D and 2D problems with sharp gradients are presented showing the effectiveness and accuracy of the new stabilized formulation. 2005 Elsevier Ltd. All rights reserved.
1. Introduction Considerable effort has been spent in recent years to derive finite element methods (FEM) [1] for the solution of the advection–diffusion–reaction equation. In this work, we will focus on the so called exponential regime originated by large absorptive (dissipative) reaction terms. Here the solutions are of the form of real exponential functions. Numerical schemes find difficulties to approximating the sharp gradients appearing in the neighborhood of boundary and internal layers due to high Peclet and/or Damko¨hler numbers. Non-physical oscillatory solution are found with the standard Galerkin FEM unless some stabilization procedure is used. Stabilized methods to tackle this problem have been based on streamline-upwind/Petrov–Galerkin (SUPG)
*
Corresponding author. Tel.: +34 93 2057016; fax: +34 93 4016517. E-mail address: [email protected] (E. On˜ate).
0045-7930/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2005.07.003
[2–7,16], Galerkin/least-squares [5–11], subgrid scale (SGS) [5,6,12,13] and residual free bubbles [14] finite element methods. While a single stabilization parameter suffices to yield stabilized (and even nodally exact results) for the one-dimensional (1D) advection–diffusion and the diffusion–reaction equations (Vol. 3 in [1] and [8–15]), this is not the case for the diffusion–advection–reaction equation. Here, in general, two stabilization parameters are needed in order to ensure a stabilized solution for all range of physical parameters and boundary conditions [4,10,14,16,17]. As reported in [12,13] the SUPG, GLS and SGS methods with a single stabilization parameter fail to obtain a stabilized solution for some specific boundary conditions in the exponential regime with negative (absorption) terms when there is a negative streamwise gradient of the solution. On˜ate et al. [18] have recently presented a stabilized FEM for the advection–diffusion absorption equation based on the use of a single stabilization parameter
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
which has the form of a diffusion term. In [18] the formulation is detailed for 1D problems and only a brief introduction to the multidimensional case is given. This paper extends the ideas presented in [18] and provides evidence of the effectiveness and accuracy of the new formulation to deal with multidimensional advection–diffusion–absorption problems with sharp gradients. The stabilized formulation is based on the standard Galerkin FEM solution of the modified governing differential equations derived via the finite calculus (FIC) method [19,20]. The FIC equations are obtained by expressing the balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitesimal theory which are a function of the balance domain dimensions. Although the FIC–FEM formulation here presented is general, we will restrict its application in this work to linear finite element approximations only. The Galerkin FIC–FEM formulation described here introduces naturally an additional nonlinear dissipation term into the discretized equations which is governed by a single stabilization parameter. In the absence of the absorption term the formulation simplifies to the standard Petrov–Galerkin approach for the advection–diffusion problem. For the diffusion–absorption case the diffusion-type stabilization term is identical to that recently obtained by Felippa and On˜ate using a variational FIC approach [15]. The general nonlinear form of the stabilization parameter is a function of the signs of the solution and its first and second derivatives. This introduces a nonlinearity in the solution scheme and a simple iterative algorithm is described. A simpler constant expression of the stabilization parameter is also presented. Details of the 1D formulation and its extension to deal with multidimensional problems are given. For the multidimensional case On˜ate et al. [30] have recently shown that a general form of the stabilization parameters can be found by writing the FIC equations along the principal curvature directions of the solution. The resulting FIC–FEM formulation is equivalent in this case (for linear elements) to adding a stabilizing diffusion matrix to the standard infinitesimal equation. The stabilizing diffusion matrix depends on the signs of the solution and its derivatives and on the velocities along the principal curvature directions of the solution. This introduces a nonlinearity in the solution process. We present a simple iterative scheme based in assuming that the main principal curvature direction at each point is coincident with the gradient vector direction. In the last part of the paper we present a collection of 1D and 2D examples showing the effectiveness and accuracy of the new FIC–FEM formulation for different values of the advective and absorptive terms.
93
2. FIC formulation of the 1D stationary advection– diffusion–absorption equation The governing equations for the 1D stationary advection–diffusion–absorption problem can be written in the FIC formulation as r
h dr ¼ 0 in x 2 ð0; LÞ 2 dx
u/ þ k
d/ h þ qp r ¼ 0 on Cq dx 2
/ /p ¼ 0
on C/
ð1Þ ð2Þ ð3Þ
where r ¼: u
d/ d d/ þ k s/ þ Q dx dx dx
ð4Þ
In above equations / is the state variable, x 2 [0, L] is the problem domain, L is the domain length, u is the velocity field, k P 0 is the diffusion, s P 0 is the absorption, dissipation or destruction source parameter, Q is a constant source term, qp and /p are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries Cq and C/, respectively and h is a characteristic length which plays the role of a stabilization parameter. In the 1D problem C/ and Cq consist of four combinations at the end points of the problem domain. Eqs. (1) and (2) are obtained by expressing the balance of fluxes in an arbitrary 1D domain of finite size within the problem domain and at the Neumann boundary, respectively. The variations of the transported variables within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitesimal theory [19,20]. The underlined stabilizing terms in Eqs. (1) and (2) emanate from these higher order expansions. Note that as the characteristic length parameter h tends to zero the FIC differential equations gradually recover the standard infinitesimal form. Successful applications of the FIC method to a variety of problems in computational mechanics can be found in [19–30,37].
3. Finite element formulation We will construct a standard finite element discretization {le} of the 1D analysis domain length L with index e ranging from 1 to the number of elements N [1]. The over the analysis state variable / is approximated by / is interpolated domain. The approximated variable / within each element with n nodes in the standard manner, i.e.
94
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
/’/
for
x 2 ½0; L
ð5aÞ
with ¼ /
n X
ð5bÞ
N i /i
i¼1
where Ni are the element shape functions and /i are Substitutnodal values of the approximate function /. ing Eq. (5a) into Eqs. (1) and (2) gives r
h dr ¼ rX 2 dx
þk u/
in x 2 ð0; LÞ
d/ h þ qp r ¼ rq dx 2
/p ¼ r / /
on C/
ð6Þ on Cq
ð7Þ ð8Þ
and rX, rq and r/ are the residuals of where r ¼ rð/Þ the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries Cq and C/, respectively. The weighted residual form of Eqs. (6)–(8) is written as Z h dr W i r dx 2 dx L d/ h ^ þ qp r þ W i u/ þ k ¼0 ð9Þ dx 2 Cq ^ i are test functions satisfying where Wi(x) and W ^ i ¼ 0 on C/. Wi ¼W Assuming smooth enough solutions and integrating by parts the term involving h in the first integral gives ^ i ¼ W i for W Z þ k d / þ qp W ir dx W i u/ dx L Cq X Z h dW i r dx ¼ 0 ð10Þ þ le 2 dx e The third term in Eq. (10) is computed as the sum of the integrals over the element interiors, therefore allowing for the space derivatives of r to be discontinuous. Also in Eq. (10) h has been assumed to be constant within each element, (i.e. dh ¼ 0 within le). dx The weak form is obtained by integrating by parts the advective and diffusive terms within r in the first integral of Eq. (10). This gives Z dW i dW i d / W i s/ þ W i Q dx ½W i qp Cq / k u dx dx dx L h dW i X Z dW i d / Q dx ¼ 0 ð11Þ bk dx dx 2 dx le e with
"
# 00 / s/ u þ 0 h b¼ 0 2k 2/ 2k /
ð12Þ
where a prime denotes the derivative with respect to the space coordinate. We see clearly that the last term of Eq. (11) introduces within each element an additional diffusion of value bk. Substituting expression (5b) into (11) and choosing a Galerkin method with Wi = Ni within each element gives the discrete system of FE equations written in the standard matrix form as ¼f K/ ð13Þ is the vector of nodal unknowns and the where / element contributions to matrix K and vector f are Z dN i dN i dN j K eij ¼ N j þ kð1 þ bÞ þ sN i N j dx u dx dx dx le Z h dN i Ni þ Q dx ðN i qp ÞCq fie ¼ 2 dx le
ð14Þ ð15Þ
The amount of balancing diffusion in Eq. (14) clearly depends on the (nonlinear) stabilization parameter b. The element and critical values of b are deduced in the next section for linear two node elements. i We note that the integral of the term h2 dN Q in Eq. dx (15) vanishes after assembly when h and Q are uniform over a patch of linear elements. 4. Computation of the stabilization parameter for linear elements The matrix Ke and the vector fe for two node linear elements are (for constant values of u, k, s and Q) 1 1 u 1 1 k sle 2 1 e e K ¼ þ þ e ð1 þ b Þ 2 1 l 6 1 2 1 1 1 ð16aÞ ( ) e Qle 1 h2 fe ¼ þ boundary term ð16bÞ e 2 1 þ h2
Fig. 1. Global axes (x, y) and principal curvature axes (n, g).
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
In Eqs. (16) index e denotes element values. Assuming Q = 0, a typical stencil for elements of equal size l can be written as / Þ ð1 þ bÞ/ þ 2ð1 þ bÞ/ cð/ iþ1 i1 i1 i w ð1 þ bÞ/iþ1 þ ð/i1 þ 4/i þ /iþ1 Þ ¼ 0 6
95
From Eq. (17) we deduce /iþ1 / w / i1 i1 þ 4/i þ /iþ1 þ 1 b¼c 2/ þ/ þ/ 6 / /iþ1 2/ i i1 iþ1 i iþ1 ð18Þ
ð17Þ
where for simplicity a constant value of b across the 2 ul and w ¼ sl4 mesh has been assumed. In Eq. (17) c ¼ 2k are the Peclet number and a velocity independent dimensionless number, respectively.
In the vecinity of a sharp gradient zone we can take / ’/ S1 / iþ1 i1 max 2/ þ/ ¼/ S2 / iþ1 i i1 max þ 4/ þ/ ¼/ S0 / ð19Þ i i iþ1 iþ1
Fig. 2. /p1 ¼ 8; /p9 ¼ 3; c ¼ 1 and x = 20. FIC results for a mesh of 8 linear elements obtained for b = 0 (Galerkin), be and bc. Comparison with the analytical solution.
Fig. 3. /p1 ¼ 8; /p9 ¼ 3; c ¼ 10 and x = 20. FIC results for a mesh of 8 linear elements obtained for b = 0 (Galerkin), be and bc. Comparison with the analytical solution.
96
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
Table 1 Comparison of 1D and 2D solutions for the advection–diffusion–absorption problem of Fig. 3 (cx = 10, w = 20) 1D
2D (nodes along line A–AÕ)
Fig. 3
4 node quads (Fig. 6)
3 node triangles (Fig. 7)
Node
¼ 0Þ /ðb
eÞ /ðb
Þ /ðb c
/ (exact)
FIC-1
FIC-2
FIC-1
FIC-2
1 2 3 4 5 6 7 8 9
8.00 2.94 1.32 1.80 0.599 0.633 1.16 1.83 3
8 3.06 1.17 0.447 0.172 0.0646 0.0264 0.0073 3
8 4 2 1 0.5 0.25 0.125 0.0625 3
8 3.08 1.19 0.457 0.176 0.0677 0.0261 0.01 3
8 3.99 2.00 1.00 0.49 0.248 0.125 0.0615 3
8 3.057 1.170 0.448 0.172 0.0648 0.0255 0.0101 3
8 4.0 2.0 1.0 0.499 0.2501 0.1250 0.0624 3
8 3.059 1.167 0.452 0.166 0.0681 0.0257 0.0072 3
Fig. 4. 2D advection–conduction–absorption problem over a square domain of size equal to 8 units. /p = 8 at x = 0, /p = 3 at x = 8, qn = 0 at y = 0 and y = 8. u = [2, 0]T, k = 1, s = 20, w = 20, cx = 1 and cy = 0. Galerkin and FIC solutions for a mesh of 8 · 8 four node square elements.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
where / max is the maximum value of the approximate in the patch of elements adjacent to the sharp function / gradient zone and 2 d/ d/ S 0 ¼ sign ð/Þ; S 1 ¼ sign ; S 2 ¼ sign dx2 dx ð20Þ denotes the sign of the magnitude within where sign ðÞ the brackets computed at the patch mid point. Substituting Eq. (19) into (18) leads to the following expression of the stabilization parameter: S0 w S1 b¼ þ c1 ð21Þ S2 6 S2 The element stabilization parameter be is now defined as be ¼ b e
b ¼0
for b > 0 for b 6 0
ð22Þ
97
where b is given by Eq. (21) and the signs S0, S1 and S2 are computed now at the element mid-point. It is clear from above that the computation of the stabilization parameter be requires the knowledge of the and that of the first sign of the numerical solution / and second derivatives of / within each element. This necessarily leads to an iterative scheme. A simple algorithm which provides a stabilized and accurate solution in just two steps is presented below. 4.1. Critical stabilization parameter and unstability conditions The following constant value of b over the mesh ensures a stabilized solution for all ranges of c and w: b 6 bc ¼ w6 þ jcj 1
Fig. 5. Solution of problem of Fig. 4 with a mesh of 8 · 8 · 2 linear triangles.
ð23Þ
98
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
where bc is the critical stabilization parameter. Note that bc corresponds to the maximum value of b in Eq. (21) for SS 02 ¼ SS 12 ¼ 1. A mathematical proof of Eq. (23) is given in [18]. Clearly the value of bc of Eq. (23) is meaningful only if bc > 0 and this can be taken as an indicator of an unstable solution. Conversely, a value of bc 6 0 indicates that no stabilization is needed.
ensures a stabilized, although sometimes slightly overdiffusive, solution. Step 2. Step 2.1. Compute the signs of the first and second 1 within each element. The second derivatives of / derivative field is obtained as follows: !e " !e !e # 1 ^1 ^1 d2 / 1 d/ d/ ¼ e ð24Þ l dx2 dx dx 2
4.2. Iterative solution scheme The following two steps iterative scheme is proposed in order to obtain a stabilized and accurate solution: 1 using Step 1. Compute a first stabilized solution / e the critical value b = bc given by Eq. (23). This
1
ð^Þei
denotes averaged values of the first where derivative at node i of element e. At the boundary nodes the constant value of the derivative of / within the element is taken in Eq. (24); i.e. ðeÞ e ð^Þi ¼ ðddx/ Þ ¼ /2 le /1 . Step 2.2. Compute the enhanced stabilized solution /2 using the element value of be given by Eq. (22).
Fig. 6. 2D advection–conduction–absorption problem over a square domain of size equal to 8 units. /p = 8 at x = 0, /p = 3 at x = 8, qn = 0 at y = 0 and y = 8. u = [20, 0]T, k = 1, s = 20, w = 20, cx = 10 and cy = 0. Galerkin and FIC solutions for a mesh of 8 · 8 four node square elements.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
In all the 1D examples solved the above two steps have sufficed to obtain a converged stabilized and accurate solution [18]. The reason of this is that the signs of the first and second derivative fields typically do not change any further after the second step over the elements adjacent to high gradient zones.
5. Extension to multi-dimensional problems Consider the general steady-state advection–diffusion–reaction equation written for the zero constant source case (Q = 0) as rð/Þ :¼ uT $/ þ $T D$/ s/ ¼ 0
u ¼ ½u; v ;
o o ; $¼ ox oy
T
;
D¼k
1
0
0
1
are respectively the velocity vector, the gradient vector and the diffusivity matrix, respectively. For simplicity we have assumed an isotropic diffusion matrix. The FIC form of Eq. (25) is written as 1 r hT $r ¼ 0 2
ð27Þ
where r is the original infinitesimal differential equation as defined in Eq. (25). The Dirichlet and boundary conditions of the FIC formulation are / /p ¼ 0
on C/
ð28Þ
ð25Þ
1 uT n/ þ nT D$/ þ qp hT nr ¼ 0 on Cq 2
ð26Þ
where n is the normal vector to the boundary where the normal flux is prescribed. As usual superindex p denotes the prescribed values.
For 2D problems T
99
Fig. 7. Solution of problem of Fig. 5 with a mesh of 8 · 8 · 2 linear triangles.
ð29Þ
100
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
In Eqs. (27) and (29) h = [hx, hy]T is the characteristic vector of the 2D FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs. (27) and (29) introduce the necessary stability in the numerical solution [19– 21]. If vector h is taken to be parallel to the velocity u the standard SUPG method is recovered [18–23]. The more general form of h allows to obtain stabilized finite element solutions in the presence of strong gradients of / near the boundaries (boundary layers) and within the analysis domain (internal layers) [30]. The FIC formulation therefore reproduces the best features of the so called shock-capturing or transverse-dissipation schemes [2,31–36]. Let us write down the FIC balance equation in the principal curvature axes of the solution n, g (Fig. 1). The FIC balance equation is
2 o/ o/ o / o2 / hn o un ug þk þ 2 s/ 2 on og og 2 on on 2 2 o/ o/ o/ o/ un ug þk þ s/ on og on2 og2 2 hg o o/ o/ o / o2 / un ug þk þ s/ ¼0 on og 2 og on2 og2 ð30Þ where un, ug are the velocities along the principal axes of curvature n and g, respectively. As n and g are the principal curvature axes of the solution then o2 / o2 / ¼ ¼0 on og og on
ð31Þ
Introducing this simplification into Eq. (30) we can rewrite this equation as
Fig. 8. Solution of problem of Fig. 4 with an unstructured mesh of 209 four node bi-linear quadrilaterals.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
101
1 ! 2 un hn shn o/ o2 / o/ kþ þ 2 2 2 on on on2 1 ! 2 ug hg shg o/ o2 / o/ þ þ kþ s/ 2 og2 2 2 og og hn o3 / hg o3 / k þ ¼0 ð32aÞ 2 on3 2 og3
values bnk and bgk along the n and g axes, respectively with 1 un hn shn o/ o2 / bn ¼ þ ; 2k 2k on on2 1 un hn shg o/ o2 / þ ð33Þ bg ¼ 2k 2k og og2
o/ o/ o2 / o2 / ug þ kð1 þ bn Þ 2 þ kð1 þ bg Þ 2 on og og on hn o3 / hg o 3 / s/ k þ ¼0 2 on3 2 og3
0 Þ$0 / s/ ¼ 0 u0T $0 / þ $0T ðD þ D
o/ o/ ug þ un on og
or un
ð32bÞ
We can see clearly from Eq. (33) that the FIC governing equations introduce orthotropic diffusion parameters of
Also note that the last term of Eq. (32b) will vanish after discretization for linear elements. Eq. (32b) can be rewritten in matrix form (neglecting the last term) as
T ½ono ; ogo , 0
ð34Þ
D is the ‘‘physical’’ where u 0 = [un, ug]T, $0 ¼ is the balancing diffuisotropic diffusion matrix and D sion matrix in the local axes n and g. The form of this matrix is
Fig. 9. Solution of problem of Fig. 4 with an unstructured mesh of 176 three node linear triangles.
102
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
" 0 ¼ D
bn k 0
0 bg k
#
ð35Þ
The velocities along the principal curvature axes un and ug can be obtained by projecting the cartesian velocities into the principal curvature axes n and g as un c a sa u u0 ¼ ¼ Tu with T ¼ ; u¼ ug sa ca v ð36Þ where ca = cos a, sa = sin a and a is the angle which the n axis forms with the x axis (Fig. 1). Note that as the solution is continuous the principal curvature directions n and g are orthogonal. The values of bn and bg are computed by considering the solution of two uncoupled 1D problems along the n and g directions. This gives from Eqs. (21) and (22)
bn ¼
sl2n S 0 wn S n1 u n ln þ ; wn ¼ cn 1 ; cn ¼ S n2 6 S n2 2k k ð37aÞ
bg ¼
2
slg S g1 S 0 wg u g lg þ ; wg ¼ cg 1 ; cg ¼ S g2 6 S g2 2k k ð37bÞ
where 2 o/ o/ S n1 ¼ sign ; S n2 ¼ sign on on2 2 o/ o/ ; S g2 ¼ sign S g1 ¼ sign ð38Þ og2 og is as usual the approximate solution. and / The lengths ln and lg are taken as the maximum projection of the velocities un and ug along the element sides S 0 ¼ sign ð/Þ;
Fig. 10. Solution of problem of Fig. 6 with an unstructured mesh of 209 four node bi-linear quadrilaterals.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
(for triangles) and the element diagonals (for quadrilaterals), i.e. li ¼ maxðdTj ui Þ;
i ¼ n; g
ð39aÞ
103
where the global diffusion matrix DG is given by DG ¼ D þ D
ð41aÞ
is where the global balancing diffusion matrix D
with j ¼ 1; 2; 3 ðfor trianglesÞand j ¼ 1; 2 ðfor quadrilateralsÞ
¼ TT D 0T D ð39bÞ
In Eq. (39a) un and ug contain the global components of the velocity vectors ~ un and ~ ug , respectively. For triangles dj are the element side vectors, whereas for quadrilaterals dj are the element diagonal vectors [30]. The next step is to transform Eq. (34) to global axes x, y. The resulting equation is uT $/ þ $T DG $/ s/ ¼ 0
ð41bÞ
ð40Þ
Remark. Similarly as for the 1D problem, a negative value of the parameters bn and bg indicates that no stabilization is needed along the directions n and g, respectively. A zero value of the corresponding stabilization parameter is chosen in this case. Remark. The expressions of bn and bg in Eqs. (37) can be simplified to
Fig. 11. Solution of problem of Fig. 6 with an unstructured mesh of 176 three node triangles.
104
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
hw i n þ jcn j 1 hw6 i ð42Þ g bg ’ bgc ¼ þ jcg j 1 6 This avoids the computation of the sign of the solution and of its first and second derivatives. The expressions of bnc and bgc in Eq. (42) are equivalent to that of the 1D critical stabilization parameter bc of Eq. (23). The main difference is that the computation of the local directions n and g is still mandatory in the multidimensional case and, therefore, the nonlinearity of the process can not be avoided here. bn ’ bnc ¼
5.1. Computation of the principal curvature axes for linear elements Excellent results have been obtained in our work by approximating the main curvature direction ~ n by the direction of the gradient vector $/.
This simplification allows us to estimate the direction ~ n in a very economical manner as the gradient vector $/ can be directly computed at any point of a linear element. Direction ~ g is taken orthogonal to that of ~ n in an anti-clockwise sense. is constant within the element. For linear triangles $/ varies linearly. We For four node quadrilaterals $/ have assumed in this case that the direction of ~ n is constant within the element and equal to the direction of computed at the element center. vector $/ The computation of the signs of the second deriva in Eq. (38) involves the following steps: (1) tive of / recovery of the cartesian first derivative field at the nodes using a nodal averaging procedure; (2) computation of the second derivative tensor at the element center and (3) transformation of this tensor to the local axes n and g. We note that in problems where positive values of / are prescribed at the Dirichlet boundary, the signs of
Fig. 12. 2D advection–diffusion–absorption problem over a square domain of size equal to 10 units. /p = 50 along x = 0 and y = 10 and /p = 100 along x = 10 and y = 0, u = [0, 15]T, k = 1, s = 12, w = 12, cx = 0 and cy = 7.5. Galerkin and FIC solutions for a mesh of 10 · 10 four node bi-linear square elements.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
S n2 , S g2 are positive due to the convexity of the numerical solution. As mentioned above the dependence of the balancing with the principal curvature direcdiffusion matrix D ~ tions n and ~ g introduces a nonlinearity in the solution process. A simple and effective iterative algorithm is described next. 5.2. General iterative scheme A stabilized numerical solution can be found by the following algorithm: Step 1. For elements in the interior of the domain choose 1n = u, i.e. the gradient direction is taken coincident with the velocity direction. If u = 0 then 1 n is taken coincident with the global x axis. In elements adjacent to a boundary choose 1n = n where n is the normal to the boundary. 0, 1D and 1DG using the expressions Compute 1 g; 1 D of bn and bg from Eq. (42).
105
that the solution is stable. This Solve for 1 /.Verify implies that there are not undershoots or overshoots in the numerical results with respect to the expected physical values. If the numerical solution is unstable, then go to step 2. Step 2. For all elements, compute at the element cen Then compute 2 g; 2 D 0, 2D and 2DG ter 2 n ¼ $1 /. using the expressions of bn and bg from Eqs. (37). Solve for 2 /. Step 3. Estimate the convergence of the process. We have chosen the following convergence norm: 1 k/k ¼ N/
max
" n X
i
/j i1 / j
2
#1=2 6e
ð43Þ
j¼1
where N is the total number of nodes in the mesh and /max is the maximum prescribed value at the Dirichlet boundary (if / max ¼ 0 then /max ¼ 1). In above steps the left upper indices denote the iteration number.In the examples shown in the next section e = 103 has
Fig. 13. Solution of the problem of Fig. 12 with an unstructured mesh of 432 four node bi-linear quadrilaterals.
106
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
been taken.If condition (43) is not satisfied, repeat steps 2 and 3 until convergence. Remark. For the advective–diffusive problems (i.e. s = 0) the expression of the balancing diffusion matrix in the interior elements for step 1 coincides with the standard (linear) SUPG form [30]. Remark. An alternative solution scheme is to use a time relaxation technique based in the solution of a pseudo transient problem with a forward Euler scheme and a diagonal mass matrix. 6. 1D numerical examples The examples presented in this section are solved in a 1D domain discretized with eight two-node linear elements of unit size. The examples are solved with the
same Dirichlet conditions /p1 ¼ 8 and /p9 ¼ 3 and two different values of c and w (c = 1, w = 20 and c = 10, w = 20). We note that for both cases the instability condition (bc > 0) is violated and, hence, the Galerkin solution should yield non-physical results. Figs. 2 and 3 show the numerical results obtained with the standard Galerkin method (b = 0) and using the element (be) and critical (bc) values of the stabilization parameter given by Eqs. (22) and (23), respectively. The exact analytical solution is also shown for comparison. Table 1 shows the nodal values of the results of the example of Fig. 3 for comparison with the 2D solutions presented in the next section. The following conclusions are drawn from the 1D results: 1. The Galerkin solution (b = 0) is unstable for both problems, as expected.
Fig. 14. Solution of problem of Fig. 12 with an structured mesh of 10 · 10 · 2 three node linear triangles.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
2. The solution obtained with the critical value bc is always stable, although it yields slightly overdiffusive results in both cases. 3. The results obtained with be are less diffusive and more accurate than those obtained with bc. The explanation is that the sign of the ratio S1/S2 is negative in the region close to the left end point of the domain. This naturally reduces the value of the stabilizing diffusion parameter b in Eq. (21) with respect to that of bc in Eq. (23) where the sign effect is not relevant. 4. The FIC algorithm provides a stabilized solution for Dirichlet problems when there is a negative streamwise gradient of the solution. This is an advantage versus SUPG, GLS and SGS methods using a single stabilization parameter which fail in some cases for these type of problems [12,13].
107
Above conclusions have been confirmed in the solution of a wider collection of 1D problems presented in [18].
7. 2D examples The analysis domain in the first two 2D examples presented is a square of size 8 units. The problems are solved first with relatively coarse meshes of 8 · 8 four node bi-linear square elements and 8 · 8 · 2 linear triangles. The boundary conditions for both examples are /p = 8 and /p = 3 at the boundaries x = 0 and x = 8, respectively and zero normal flux at y = 0 and y = 8. This reproduces the condition of the two 1D examples solved in the previous section. The first example is analyzed for
Fig. 15. Solution of problem of Fig. 12 with an unstructured mesh of 780 three node triangles.
108
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
u = [2, 0]T, k = 1 and s = 20 giving w = 20, cx = 1 and cy = 0 which corresponds to the first 1D example (Fig. 2). The correct solution for this problem has a boundary layer in the vecinity of the two sides at x = 0 and x = 8 where / is prescribed (Fig. 4). The numerical results obtained with the standard Galerkin solution are oscillatory as expected. The stabilized FIC formulation eliminates the oscillations and yields the correct physical solution. Good results are obtained for both meshes of linear rectangles and triangles (Figs. 4 and 5). Results labelled as FIC-1 and FIC-2 in the figures correspond to those obtained in the first and second iteration of the algorithm presented in Section 5.2, respectively. We note that the FIC-1 results agree precisely with those obtained in the 1D case for b = bc, whereas the FIC-2 results agree with the more accurate 1D values obtained with the element stabilization parameter be (see Fig. 2). The second example is similar to the first one with u = [20, 0]T, k = 1 and s = 20 giving w = 20, cx = 10
and cy = 0. These values correspond to the second 1D problem of the previous section (Fig. 3). The Galerkin solution is again oscillatory, whereas the FIC results are physically sound (Figs. 6 and 7). Once more the FIC-1 and FIC-2 results are in good agreement with the 1D values shown in Fig. 3 for bc and be, respectively for both meshes of square and triangular elements. The coincidence of the 1D and 2D results for this problem can be clearly seen in Table 1. Note that, similarly to the 1D case, the FIC-2 results are more accurate (less diffusive) than those obtained in the first iteration (FIC-1). This is due to the more precise evaluation of bn and bg in Eqs. (37) accounting for the correct sign of all the terms. Figs. 8–11 show results for the two 2D problems above described solved now with relatively coarse unstructured meshes of linear triangles and quadrilaterals. The effectiveness and accuracy of the FIC iterative scheme is again noticeable in all cases. Note the agreement of the FIC-2 results of Figs. 10 and 11 with
Fig. 16. Diffusive–absorptive problem over a square domain of size equal to 10 units. /p = 50 over x = 0 and y = 10 and /p = 100 over x = 10 and y = 0, u = [0, 0]T, k = 1, s = 20, w = 20, cx = 0 and cy = 0. Galerkin and FIC solutions obtained with structured meshes of four node quadrilaterals and linear triangles.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
the exact solution for the equivalent 1D problem of Fig. 3. Fig. 12 presents the solution of a similar problem where the values of / are prescribed at the four boundaries. The solution domain has now 10 units and the problem is solved first with a mesh of 10 · 10 four node square elements. Details of the physical parameters are given in Fig. 12. Excellent results are again obtained with the FIC scheme. Similar good results are obtained with a structured mesh of linear triangles (Fig. 13) as well as with non-structured meshes of linear quadrilateral and triangles (Figs. 14 and 15). The effectiveness of the FIC scheme for a diffusive– absorptive problem with Dirichlet boundary conditions is shown in Fig. 16. The results shown have been obtained with structured meshes of linear quadrilateral and triangles. Note that the four boundary layers are well captured in the first step of the iterative solution.
109
Similar good results have also been obtained with unstructured meshes not shown here. The final example is a standard benchmark problem of advection–diffusion where sharp layers appear at both the boundary and the interior of the domain. The problem is the advective–diffusive transport of / in a square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source terms (i.e. Q = 0 and s = 0). Fig. 17 displays the SUPG solution and FIC results obtained after two iterations using a structured mesh of 20 · 20 linear four node square elements. It is remarkable that the FIC results capture the sharp gradient zones at the boundaries where / is prescribed to zero and at the interior of the domain and eliminate all the spurious oscillations present in the SUPG method. Similar good results obtained with the FIC method for a wide range of advective–diffusive problems are presented in [30]. Recent applications of the FIC
Fig. 16 (continued)
110
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
Fig. 17. Square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 20 · 20 linear four node square elements.
method to incompressible fluid flow problems are reported in [37].
8. Conclusions The FIC–FEM formulation presented allows to obtain a stabilized and accurate solution for the advection–diffusion–absorption equation. For the 1D problem the formulation is equivalent to adding a nonlinear diffusion term to the standard discretized equations which is governed by a single stabilization parameter. The use of the constant critical value of the 1D stabilization parameter provides a stabilized numerical solution in a single step. A more accurate (less diffusive) solution can be obtained using the two step iterative scheme proposed. The equivalence of the FIC method with a nonlinear stabilizing diffusion term extends naturally to multidi-
mensional problems using structured and unstructured meshes. The key step is to express the governing equations of the FIC formulation in the principal curvature directions of the solution. The resulting FIC equation is equivalent to adding a nonlinear diffusion matrix to the infinitesimal governing equations. The solution process becomes nonlinear and a simple iterative algorithm has been presented. The approximation of the main principal curvature direction by that of the gradient vector simplifies the computations in the iterative scheme. Excellent results have been obtained for all the 2D problems solved in just two iterations for structured and nonstructured meshes. It is remarkable that, similarly to the 1D case, good stabilized results are obtained in the first iteration of the scheme proposed (FIC-1 results) and this may be sufficient for many practical cases. More accurate (less diffusive) results are obtained by performing a second iteration at a relatively small additional computational cost.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
Acknowledgements The authors also thank Profs. C. Felippa and S.R. Idelsohn for many useful discussions. This work has been sponsored by the Ministerio de Educacio´n y Ciencia of Spain. Plan Nacional, Project numbers: DPI2001-2193, BIA2003-09078-C02-01, and DPI2004-07410-C03-02.
References [1] Zienkiewicz OC, Taylor RL. fifth ed. The finite element method, vol. 3. Butterworth-Heinemann; 2001. [2] Tezduyar TE, Park YJ. Discontinuity-capturing finite element formulations for nonlinear convection–diffusion–reaction equations. Comput Methods Appl Mech Engrg 1986;59:307–25. [3] Tezduyar TE, Park YJ, Deans HA. Finite element procedures for time-dependent convection–diffusion–reaction systems. Int J Numer Meth Fluids 1987;7:1013–33. [4] Idelsohn S, Nigro N, Storti M, Buscaglia G. A Petrov–Galerkin formulation for advection–reaction–diffusion problems. Comput. Methods Appl Mech Engrg 1996;136:27–46. [5] Codina R. Comparison of some finite element methods for solving the diffusion–convection–reaction equation. Comput Methods Appl Mech Engrg 1998;156:186–210. [6] Codina R. On stabilized finite element methods for linear systems of convection–diffusion–reaction equations. Comput Meth Appl Mech Engrg 2000;188:61–82. [7] Burman E, Ern A. Nonlinear diffusion and discrete maxium principle for stabilized Galerkin approximations of the convection–diffusion–reaction equation. Comput Meth Appl Mech Engrg 2002;191:3833–55. [8] Harari I, Hughes TJR. Finite element methods for the Helmholtz equation in an exterior domain: model problems. Comput Meth Appl Mech Engrg 1991;87:59–96. [9] Harari I, Hughes TJR. Galerkin/least squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains. Comput Meth Appl Mech Engrg 1992;98:411–54. [10] Harari I, Hughes TJR. Stabilized finite element method for steady advection–diffusion with production. Comput Meth Appl Mech Engrg 1994;115:165–91. [11] Harari I, Grosh K, Hughes TJR, Malhotra M, Pinsky PM, Stewart JR, et al. Recent development in finite element methods for structural acoustics. Arch Comput Mech Eng 1996;3(2-3): 131–309. [12] Hauke G, Garcia-Olivares A. Variational subgrid scale formulations for the advection–diffusion–reaction equation. Comput Methods Appl Mech Engrg 2000;190:6847–65. [13] Hauke G. A simple subgrid scale stabilized method for the advection–diffusion reaction equation. Comput. Methods Appl Mech Engrg 2002;191:2925–47. [14] Brezzi F, Hauke G, Marin LD, Sangalli S. Link-cutting bubbles for the stabilization of convection–diffusion–reaction problems. Math Models Methods Appl Sci. World Scientific Publishing Company; 2002. [15] Felippa CA, On˜ate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Comput Mech, in print, CM-050155.
111
[16] Idelsohn SR, Heinrich JC, On˜ate E. Petrov–Galerkin methods for the transient advective–diffusive equation with sharp gradients. Int J Numer Meth Engng 1996;39:1455–73. [17] Harari I. Stability of semidiscrete advection–diffusion in transient computation. In: Yao ZH, Yuan MW, Zhong WX, editors. Proceedings of 6th World Congress on Computational Mechanics, Beijing, September. Tsinghua Univ. Press-Springer; 2004. [18] On˜ate E, Miquel J, Hauke G. Stabilized formulation for the advection–diffusion–reaction equations using finite calculus and linear finite elements. Comput Methods Appl Mech Engrg, submitted for publication. [19] On˜ate E. Derivation of stabilized equations for advective–diffusive transport and fluid flow problems. Comput Methods Appl Mech Engrg 1998;151(1-2):233–67. [20] On˜ate E. Possibilities of finite calculus in computational mechanics. Int J Numer Meth Engng 2004;60:255–81. [21] On˜ate E, Manzan M. A general procedure for deriving stabilized space–time finite element methods for advective–diffusive problems. Int J Numer Meth Fluids 1999;31:203–21. [22] On˜ate E, Idelsohn SR, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics. Applications to convective transport and fluid flow. Int J Numer Meth Engng 1996;39: 3839–66. [23] On˜ate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comput Methods Appl Mech Engrg 2000;182(1-2): 355–70. [24] On˜ate E, Garcı´a J. A finite element method for fluid–structure interaction with surface waves using a finite calculus formulation. Comput Methods Appl Mech Engrg 2001;191(6–7): 635–60. [25] Idelsohn SR, On˜ate E, Del Pin F. A Lagrangian meshless finite element method applied to fluid–structure interaction problems. Comput Struct 2003;81:655–71. [26] On˜ate E, Idelsohn SR, Del Pin F, Aubry R. The particle finite element method. An overview (PFEM). Int J Comput Methods 2004;1(2):267–307. [27] On˜ate E, Garcı´a J, Idelsohn SR. Ship hydrodynamics. In: Hughes T, de Borst R, Stein E, editors. Encyclopedia of computational mechanics, Vol. 3. J. Wiley; 2004. p. 579–607. Chapter 18. [28] On˜ate E, Taylor RL, Zienkiewicz OC, Rojek J. A residual correction method based on finite calculus. Eng Comput 2003;20(5/6):629–58. [29] On˜ate E, Rojek J, Taylor RL, Zienkiewicz OC. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. Int J Numer Meth Engng 2004; 59(11):1473–500. [30] On˜ate E, Za´rate F, Idelsohn SR. Finite element formulation for convection–diffusion problems with sharp gradients using finite calculus. Comput Meth Appl Mech Engng, submitted for publication. [31] Hughes TJR, Mallet M. A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective–diffusive system. Comput Methods Appl Mech Engrg 1986;58:329–36. [32] Codina R. A discontinuity-capturing crosswind dissipation for the finite element solution of the convection–diffusion equation. Comput Methods Appl Mech Engrg 1993;110: 325– 42. [33] Hughes TJR, Mallet M, Mizukami A. A new finite element formulation for comptutational fluid dynamics: II Beyond SUPG. Comput Methods Appl Mech Engrg 1986;54:341–55.
112
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
[34] Galea˜o AC, Dutra do Carmo EG. A consistent approximate upwind Petrov–Galerkin method for convection-dominated problems. Comput Methods Appl Mech Engrg 1988;68:83–95. [35] Tezduyar TE. Adaptive determination of the finite element stabilization parameters. Proceedings of the ECCOMAS Computational Fluid Dynamics Conference 2001, Swansea, Wales, UK, CD-Rom, 2001.
[36] Tezduyar TE. Computation of moving boundaries and interfaces and stabilization parameters. Int J Numer Methods Fluid 2003;43:555–75. [37] On˜ate E, Garcı´a J, Idelsohn S, Del Pin F. FIC formulation for finite element analysis of incompressible flows. Eulerian, Lagrangian and ALE approaches. Comput Methods Appl Mech Eng, in press.
Stabilized solution of the multidimensional advection–diffusion–absorption equation using linear finite elements Eugenio On˜ate *, Juan Miquel, Francisco Za´rate International Center for Numerical Methods in Engineering (CIMNE), Universidad Polite´cnica de Catalun˜a, Edificio C1, Gran Capita´n s/n, 08034 Barcelona, Spain Received 23 June 2005; accepted 21 July 2005 Available online 1 December 2005
Abstract A stabilized finite element method (FEM) for the multidimensional steady state advection–diffusion–absorption equation is presented. The stabilized formulation is based on the modified governing differential equations derived via the finite calculus (FIC) method. For 1D problems the stabilization terms act as a nonlinear additional diffusion governed by a single stabilization parameter. It is shown that for multidimensional problems an orthotropic stabilizing diffusion must be added along the principal directions of curvature of the solution. A simple iterative algorithm yielding a stable and accurate solution for all the range of physical parameters and boundary conditions is described. Numerical results for 1D and 2D problems with sharp gradients are presented showing the effectiveness and accuracy of the new stabilized formulation. 2005 Elsevier Ltd. All rights reserved.
1. Introduction Considerable effort has been spent in recent years to derive finite element methods (FEM) [1] for the solution of the advection–diffusion–reaction equation. In this work, we will focus on the so called exponential regime originated by large absorptive (dissipative) reaction terms. Here the solutions are of the form of real exponential functions. Numerical schemes find difficulties to approximating the sharp gradients appearing in the neighborhood of boundary and internal layers due to high Peclet and/or Damko¨hler numbers. Non-physical oscillatory solution are found with the standard Galerkin FEM unless some stabilization procedure is used. Stabilized methods to tackle this problem have been based on streamline-upwind/Petrov–Galerkin (SUPG)
*
Corresponding author. Tel.: +34 93 2057016; fax: +34 93 4016517. E-mail address: [email protected] (E. On˜ate).
0045-7930/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2005.07.003
[2–7,16], Galerkin/least-squares [5–11], subgrid scale (SGS) [5,6,12,13] and residual free bubbles [14] finite element methods. While a single stabilization parameter suffices to yield stabilized (and even nodally exact results) for the one-dimensional (1D) advection–diffusion and the diffusion–reaction equations (Vol. 3 in [1] and [8–15]), this is not the case for the diffusion–advection–reaction equation. Here, in general, two stabilization parameters are needed in order to ensure a stabilized solution for all range of physical parameters and boundary conditions [4,10,14,16,17]. As reported in [12,13] the SUPG, GLS and SGS methods with a single stabilization parameter fail to obtain a stabilized solution for some specific boundary conditions in the exponential regime with negative (absorption) terms when there is a negative streamwise gradient of the solution. On˜ate et al. [18] have recently presented a stabilized FEM for the advection–diffusion absorption equation based on the use of a single stabilization parameter
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
which has the form of a diffusion term. In [18] the formulation is detailed for 1D problems and only a brief introduction to the multidimensional case is given. This paper extends the ideas presented in [18] and provides evidence of the effectiveness and accuracy of the new formulation to deal with multidimensional advection–diffusion–absorption problems with sharp gradients. The stabilized formulation is based on the standard Galerkin FEM solution of the modified governing differential equations derived via the finite calculus (FIC) method [19,20]. The FIC equations are obtained by expressing the balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitesimal theory which are a function of the balance domain dimensions. Although the FIC–FEM formulation here presented is general, we will restrict its application in this work to linear finite element approximations only. The Galerkin FIC–FEM formulation described here introduces naturally an additional nonlinear dissipation term into the discretized equations which is governed by a single stabilization parameter. In the absence of the absorption term the formulation simplifies to the standard Petrov–Galerkin approach for the advection–diffusion problem. For the diffusion–absorption case the diffusion-type stabilization term is identical to that recently obtained by Felippa and On˜ate using a variational FIC approach [15]. The general nonlinear form of the stabilization parameter is a function of the signs of the solution and its first and second derivatives. This introduces a nonlinearity in the solution scheme and a simple iterative algorithm is described. A simpler constant expression of the stabilization parameter is also presented. Details of the 1D formulation and its extension to deal with multidimensional problems are given. For the multidimensional case On˜ate et al. [30] have recently shown that a general form of the stabilization parameters can be found by writing the FIC equations along the principal curvature directions of the solution. The resulting FIC–FEM formulation is equivalent in this case (for linear elements) to adding a stabilizing diffusion matrix to the standard infinitesimal equation. The stabilizing diffusion matrix depends on the signs of the solution and its derivatives and on the velocities along the principal curvature directions of the solution. This introduces a nonlinearity in the solution process. We present a simple iterative scheme based in assuming that the main principal curvature direction at each point is coincident with the gradient vector direction. In the last part of the paper we present a collection of 1D and 2D examples showing the effectiveness and accuracy of the new FIC–FEM formulation for different values of the advective and absorptive terms.
93
2. FIC formulation of the 1D stationary advection– diffusion–absorption equation The governing equations for the 1D stationary advection–diffusion–absorption problem can be written in the FIC formulation as r
h dr ¼ 0 in x 2 ð0; LÞ 2 dx
u/ þ k
d/ h þ qp r ¼ 0 on Cq dx 2
/ /p ¼ 0
on C/
ð1Þ ð2Þ ð3Þ
where r ¼: u
d/ d d/ þ k s/ þ Q dx dx dx
ð4Þ
In above equations / is the state variable, x 2 [0, L] is the problem domain, L is the domain length, u is the velocity field, k P 0 is the diffusion, s P 0 is the absorption, dissipation or destruction source parameter, Q is a constant source term, qp and /p are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries Cq and C/, respectively and h is a characteristic length which plays the role of a stabilization parameter. In the 1D problem C/ and Cq consist of four combinations at the end points of the problem domain. Eqs. (1) and (2) are obtained by expressing the balance of fluxes in an arbitrary 1D domain of finite size within the problem domain and at the Neumann boundary, respectively. The variations of the transported variables within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitesimal theory [19,20]. The underlined stabilizing terms in Eqs. (1) and (2) emanate from these higher order expansions. Note that as the characteristic length parameter h tends to zero the FIC differential equations gradually recover the standard infinitesimal form. Successful applications of the FIC method to a variety of problems in computational mechanics can be found in [19–30,37].
3. Finite element formulation We will construct a standard finite element discretization {le} of the 1D analysis domain length L with index e ranging from 1 to the number of elements N [1]. The over the analysis state variable / is approximated by / is interpolated domain. The approximated variable / within each element with n nodes in the standard manner, i.e.
94
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
/’/
for
x 2 ½0; L
ð5aÞ
with ¼ /
n X
ð5bÞ
N i /i
i¼1
where Ni are the element shape functions and /i are Substitutnodal values of the approximate function /. ing Eq. (5a) into Eqs. (1) and (2) gives r
h dr ¼ rX 2 dx
þk u/
in x 2 ð0; LÞ
d/ h þ qp r ¼ rq dx 2
/p ¼ r / /
on C/
ð6Þ on Cq
ð7Þ ð8Þ
and rX, rq and r/ are the residuals of where r ¼ rð/Þ the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries Cq and C/, respectively. The weighted residual form of Eqs. (6)–(8) is written as Z h dr W i r dx 2 dx L d/ h ^ þ qp r þ W i u/ þ k ¼0 ð9Þ dx 2 Cq ^ i are test functions satisfying where Wi(x) and W ^ i ¼ 0 on C/. Wi ¼W Assuming smooth enough solutions and integrating by parts the term involving h in the first integral gives ^ i ¼ W i for W Z þ k d / þ qp W ir dx W i u/ dx L Cq X Z h dW i r dx ¼ 0 ð10Þ þ le 2 dx e The third term in Eq. (10) is computed as the sum of the integrals over the element interiors, therefore allowing for the space derivatives of r to be discontinuous. Also in Eq. (10) h has been assumed to be constant within each element, (i.e. dh ¼ 0 within le). dx The weak form is obtained by integrating by parts the advective and diffusive terms within r in the first integral of Eq. (10). This gives Z dW i dW i d / W i s/ þ W i Q dx ½W i qp Cq / k u dx dx dx L h dW i X Z dW i d / Q dx ¼ 0 ð11Þ bk dx dx 2 dx le e with
"
# 00 / s/ u þ 0 h b¼ 0 2k 2/ 2k /
ð12Þ
where a prime denotes the derivative with respect to the space coordinate. We see clearly that the last term of Eq. (11) introduces within each element an additional diffusion of value bk. Substituting expression (5b) into (11) and choosing a Galerkin method with Wi = Ni within each element gives the discrete system of FE equations written in the standard matrix form as ¼f K/ ð13Þ is the vector of nodal unknowns and the where / element contributions to matrix K and vector f are Z dN i dN i dN j K eij ¼ N j þ kð1 þ bÞ þ sN i N j dx u dx dx dx le Z h dN i Ni þ Q dx ðN i qp ÞCq fie ¼ 2 dx le
ð14Þ ð15Þ
The amount of balancing diffusion in Eq. (14) clearly depends on the (nonlinear) stabilization parameter b. The element and critical values of b are deduced in the next section for linear two node elements. i We note that the integral of the term h2 dN Q in Eq. dx (15) vanishes after assembly when h and Q are uniform over a patch of linear elements. 4. Computation of the stabilization parameter for linear elements The matrix Ke and the vector fe for two node linear elements are (for constant values of u, k, s and Q) 1 1 u 1 1 k sle 2 1 e e K ¼ þ þ e ð1 þ b Þ 2 1 l 6 1 2 1 1 1 ð16aÞ ( ) e Qle 1 h2 fe ¼ þ boundary term ð16bÞ e 2 1 þ h2
Fig. 1. Global axes (x, y) and principal curvature axes (n, g).
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
In Eqs. (16) index e denotes element values. Assuming Q = 0, a typical stencil for elements of equal size l can be written as / Þ ð1 þ bÞ/ þ 2ð1 þ bÞ/ cð/ iþ1 i1 i1 i w ð1 þ bÞ/iþ1 þ ð/i1 þ 4/i þ /iþ1 Þ ¼ 0 6
95
From Eq. (17) we deduce /iþ1 / w / i1 i1 þ 4/i þ /iþ1 þ 1 b¼c 2/ þ/ þ/ 6 / /iþ1 2/ i i1 iþ1 i iþ1 ð18Þ
ð17Þ
where for simplicity a constant value of b across the 2 ul and w ¼ sl4 mesh has been assumed. In Eq. (17) c ¼ 2k are the Peclet number and a velocity independent dimensionless number, respectively.
In the vecinity of a sharp gradient zone we can take / ’/ S1 / iþ1 i1 max 2/ þ/ ¼/ S2 / iþ1 i i1 max þ 4/ þ/ ¼/ S0 / ð19Þ i i iþ1 iþ1
Fig. 2. /p1 ¼ 8; /p9 ¼ 3; c ¼ 1 and x = 20. FIC results for a mesh of 8 linear elements obtained for b = 0 (Galerkin), be and bc. Comparison with the analytical solution.
Fig. 3. /p1 ¼ 8; /p9 ¼ 3; c ¼ 10 and x = 20. FIC results for a mesh of 8 linear elements obtained for b = 0 (Galerkin), be and bc. Comparison with the analytical solution.
96
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
Table 1 Comparison of 1D and 2D solutions for the advection–diffusion–absorption problem of Fig. 3 (cx = 10, w = 20) 1D
2D (nodes along line A–AÕ)
Fig. 3
4 node quads (Fig. 6)
3 node triangles (Fig. 7)
Node
¼ 0Þ /ðb
eÞ /ðb
Þ /ðb c
/ (exact)
FIC-1
FIC-2
FIC-1
FIC-2
1 2 3 4 5 6 7 8 9
8.00 2.94 1.32 1.80 0.599 0.633 1.16 1.83 3
8 3.06 1.17 0.447 0.172 0.0646 0.0264 0.0073 3
8 4 2 1 0.5 0.25 0.125 0.0625 3
8 3.08 1.19 0.457 0.176 0.0677 0.0261 0.01 3
8 3.99 2.00 1.00 0.49 0.248 0.125 0.0615 3
8 3.057 1.170 0.448 0.172 0.0648 0.0255 0.0101 3
8 4.0 2.0 1.0 0.499 0.2501 0.1250 0.0624 3
8 3.059 1.167 0.452 0.166 0.0681 0.0257 0.0072 3
Fig. 4. 2D advection–conduction–absorption problem over a square domain of size equal to 8 units. /p = 8 at x = 0, /p = 3 at x = 8, qn = 0 at y = 0 and y = 8. u = [2, 0]T, k = 1, s = 20, w = 20, cx = 1 and cy = 0. Galerkin and FIC solutions for a mesh of 8 · 8 four node square elements.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
where / max is the maximum value of the approximate in the patch of elements adjacent to the sharp function / gradient zone and 2 d/ d/ S 0 ¼ sign ð/Þ; S 1 ¼ sign ; S 2 ¼ sign dx2 dx ð20Þ denotes the sign of the magnitude within where sign ðÞ the brackets computed at the patch mid point. Substituting Eq. (19) into (18) leads to the following expression of the stabilization parameter: S0 w S1 b¼ þ c1 ð21Þ S2 6 S2 The element stabilization parameter be is now defined as be ¼ b e
b ¼0
for b > 0 for b 6 0
ð22Þ
97
where b is given by Eq. (21) and the signs S0, S1 and S2 are computed now at the element mid-point. It is clear from above that the computation of the stabilization parameter be requires the knowledge of the and that of the first sign of the numerical solution / and second derivatives of / within each element. This necessarily leads to an iterative scheme. A simple algorithm which provides a stabilized and accurate solution in just two steps is presented below. 4.1. Critical stabilization parameter and unstability conditions The following constant value of b over the mesh ensures a stabilized solution for all ranges of c and w: b 6 bc ¼ w6 þ jcj 1
Fig. 5. Solution of problem of Fig. 4 with a mesh of 8 · 8 · 2 linear triangles.
ð23Þ
98
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
where bc is the critical stabilization parameter. Note that bc corresponds to the maximum value of b in Eq. (21) for SS 02 ¼ SS 12 ¼ 1. A mathematical proof of Eq. (23) is given in [18]. Clearly the value of bc of Eq. (23) is meaningful only if bc > 0 and this can be taken as an indicator of an unstable solution. Conversely, a value of bc 6 0 indicates that no stabilization is needed.
ensures a stabilized, although sometimes slightly overdiffusive, solution. Step 2. Step 2.1. Compute the signs of the first and second 1 within each element. The second derivatives of / derivative field is obtained as follows: !e " !e !e # 1 ^1 ^1 d2 / 1 d/ d/ ¼ e ð24Þ l dx2 dx dx 2
4.2. Iterative solution scheme The following two steps iterative scheme is proposed in order to obtain a stabilized and accurate solution: 1 using Step 1. Compute a first stabilized solution / e the critical value b = bc given by Eq. (23). This
1
ð^Þei
denotes averaged values of the first where derivative at node i of element e. At the boundary nodes the constant value of the derivative of / within the element is taken in Eq. (24); i.e. ðeÞ e ð^Þi ¼ ðddx/ Þ ¼ /2 le /1 . Step 2.2. Compute the enhanced stabilized solution /2 using the element value of be given by Eq. (22).
Fig. 6. 2D advection–conduction–absorption problem over a square domain of size equal to 8 units. /p = 8 at x = 0, /p = 3 at x = 8, qn = 0 at y = 0 and y = 8. u = [20, 0]T, k = 1, s = 20, w = 20, cx = 10 and cy = 0. Galerkin and FIC solutions for a mesh of 8 · 8 four node square elements.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
In all the 1D examples solved the above two steps have sufficed to obtain a converged stabilized and accurate solution [18]. The reason of this is that the signs of the first and second derivative fields typically do not change any further after the second step over the elements adjacent to high gradient zones.
5. Extension to multi-dimensional problems Consider the general steady-state advection–diffusion–reaction equation written for the zero constant source case (Q = 0) as rð/Þ :¼ uT $/ þ $T D$/ s/ ¼ 0
u ¼ ½u; v ;
o o ; $¼ ox oy
T
;
D¼k
1
0
0
1
are respectively the velocity vector, the gradient vector and the diffusivity matrix, respectively. For simplicity we have assumed an isotropic diffusion matrix. The FIC form of Eq. (25) is written as 1 r hT $r ¼ 0 2
ð27Þ
where r is the original infinitesimal differential equation as defined in Eq. (25). The Dirichlet and boundary conditions of the FIC formulation are / /p ¼ 0
on C/
ð28Þ
ð25Þ
1 uT n/ þ nT D$/ þ qp hT nr ¼ 0 on Cq 2
ð26Þ
where n is the normal vector to the boundary where the normal flux is prescribed. As usual superindex p denotes the prescribed values.
For 2D problems T
99
Fig. 7. Solution of problem of Fig. 5 with a mesh of 8 · 8 · 2 linear triangles.
ð29Þ
100
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
In Eqs. (27) and (29) h = [hx, hy]T is the characteristic vector of the 2D FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs. (27) and (29) introduce the necessary stability in the numerical solution [19– 21]. If vector h is taken to be parallel to the velocity u the standard SUPG method is recovered [18–23]. The more general form of h allows to obtain stabilized finite element solutions in the presence of strong gradients of / near the boundaries (boundary layers) and within the analysis domain (internal layers) [30]. The FIC formulation therefore reproduces the best features of the so called shock-capturing or transverse-dissipation schemes [2,31–36]. Let us write down the FIC balance equation in the principal curvature axes of the solution n, g (Fig. 1). The FIC balance equation is
2 o/ o/ o / o2 / hn o un ug þk þ 2 s/ 2 on og og 2 on on 2 2 o/ o/ o/ o/ un ug þk þ s/ on og on2 og2 2 hg o o/ o/ o / o2 / un ug þk þ s/ ¼0 on og 2 og on2 og2 ð30Þ where un, ug are the velocities along the principal axes of curvature n and g, respectively. As n and g are the principal curvature axes of the solution then o2 / o2 / ¼ ¼0 on og og on
ð31Þ
Introducing this simplification into Eq. (30) we can rewrite this equation as
Fig. 8. Solution of problem of Fig. 4 with an unstructured mesh of 209 four node bi-linear quadrilaterals.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
101
1 ! 2 un hn shn o/ o2 / o/ kþ þ 2 2 2 on on on2 1 ! 2 ug hg shg o/ o2 / o/ þ þ kþ s/ 2 og2 2 2 og og hn o3 / hg o3 / k þ ¼0 ð32aÞ 2 on3 2 og3
values bnk and bgk along the n and g axes, respectively with 1 un hn shn o/ o2 / bn ¼ þ ; 2k 2k on on2 1 un hn shg o/ o2 / þ ð33Þ bg ¼ 2k 2k og og2
o/ o/ o2 / o2 / ug þ kð1 þ bn Þ 2 þ kð1 þ bg Þ 2 on og og on hn o3 / hg o 3 / s/ k þ ¼0 2 on3 2 og3
0 Þ$0 / s/ ¼ 0 u0T $0 / þ $0T ðD þ D
o/ o/ ug þ un on og
or un
ð32bÞ
We can see clearly from Eq. (33) that the FIC governing equations introduce orthotropic diffusion parameters of
Also note that the last term of Eq. (32b) will vanish after discretization for linear elements. Eq. (32b) can be rewritten in matrix form (neglecting the last term) as
T ½ono ; ogo , 0
ð34Þ
D is the ‘‘physical’’ where u 0 = [un, ug]T, $0 ¼ is the balancing diffuisotropic diffusion matrix and D sion matrix in the local axes n and g. The form of this matrix is
Fig. 9. Solution of problem of Fig. 4 with an unstructured mesh of 176 three node linear triangles.
102
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
" 0 ¼ D
bn k 0
0 bg k
#
ð35Þ
The velocities along the principal curvature axes un and ug can be obtained by projecting the cartesian velocities into the principal curvature axes n and g as un c a sa u u0 ¼ ¼ Tu with T ¼ ; u¼ ug sa ca v ð36Þ where ca = cos a, sa = sin a and a is the angle which the n axis forms with the x axis (Fig. 1). Note that as the solution is continuous the principal curvature directions n and g are orthogonal. The values of bn and bg are computed by considering the solution of two uncoupled 1D problems along the n and g directions. This gives from Eqs. (21) and (22)
bn ¼
sl2n S 0 wn S n1 u n ln þ ; wn ¼ cn 1 ; cn ¼ S n2 6 S n2 2k k ð37aÞ
bg ¼
2
slg S g1 S 0 wg u g lg þ ; wg ¼ cg 1 ; cg ¼ S g2 6 S g2 2k k ð37bÞ
where 2 o/ o/ S n1 ¼ sign ; S n2 ¼ sign on on2 2 o/ o/ ; S g2 ¼ sign S g1 ¼ sign ð38Þ og2 og is as usual the approximate solution. and / The lengths ln and lg are taken as the maximum projection of the velocities un and ug along the element sides S 0 ¼ sign ð/Þ;
Fig. 10. Solution of problem of Fig. 6 with an unstructured mesh of 209 four node bi-linear quadrilaterals.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
(for triangles) and the element diagonals (for quadrilaterals), i.e. li ¼ maxðdTj ui Þ;
i ¼ n; g
ð39aÞ
103
where the global diffusion matrix DG is given by DG ¼ D þ D
ð41aÞ
is where the global balancing diffusion matrix D
with j ¼ 1; 2; 3 ðfor trianglesÞand j ¼ 1; 2 ðfor quadrilateralsÞ
¼ TT D 0T D ð39bÞ
In Eq. (39a) un and ug contain the global components of the velocity vectors ~ un and ~ ug , respectively. For triangles dj are the element side vectors, whereas for quadrilaterals dj are the element diagonal vectors [30]. The next step is to transform Eq. (34) to global axes x, y. The resulting equation is uT $/ þ $T DG $/ s/ ¼ 0
ð41bÞ
ð40Þ
Remark. Similarly as for the 1D problem, a negative value of the parameters bn and bg indicates that no stabilization is needed along the directions n and g, respectively. A zero value of the corresponding stabilization parameter is chosen in this case. Remark. The expressions of bn and bg in Eqs. (37) can be simplified to
Fig. 11. Solution of problem of Fig. 6 with an unstructured mesh of 176 three node triangles.
104
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
hw i n þ jcn j 1 hw6 i ð42Þ g bg ’ bgc ¼ þ jcg j 1 6 This avoids the computation of the sign of the solution and of its first and second derivatives. The expressions of bnc and bgc in Eq. (42) are equivalent to that of the 1D critical stabilization parameter bc of Eq. (23). The main difference is that the computation of the local directions n and g is still mandatory in the multidimensional case and, therefore, the nonlinearity of the process can not be avoided here. bn ’ bnc ¼
5.1. Computation of the principal curvature axes for linear elements Excellent results have been obtained in our work by approximating the main curvature direction ~ n by the direction of the gradient vector $/.
This simplification allows us to estimate the direction ~ n in a very economical manner as the gradient vector $/ can be directly computed at any point of a linear element. Direction ~ g is taken orthogonal to that of ~ n in an anti-clockwise sense. is constant within the element. For linear triangles $/ varies linearly. We For four node quadrilaterals $/ have assumed in this case that the direction of ~ n is constant within the element and equal to the direction of computed at the element center. vector $/ The computation of the signs of the second deriva in Eq. (38) involves the following steps: (1) tive of / recovery of the cartesian first derivative field at the nodes using a nodal averaging procedure; (2) computation of the second derivative tensor at the element center and (3) transformation of this tensor to the local axes n and g. We note that in problems where positive values of / are prescribed at the Dirichlet boundary, the signs of
Fig. 12. 2D advection–diffusion–absorption problem over a square domain of size equal to 10 units. /p = 50 along x = 0 and y = 10 and /p = 100 along x = 10 and y = 0, u = [0, 15]T, k = 1, s = 12, w = 12, cx = 0 and cy = 7.5. Galerkin and FIC solutions for a mesh of 10 · 10 four node bi-linear square elements.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
S n2 , S g2 are positive due to the convexity of the numerical solution. As mentioned above the dependence of the balancing with the principal curvature direcdiffusion matrix D ~ tions n and ~ g introduces a nonlinearity in the solution process. A simple and effective iterative algorithm is described next. 5.2. General iterative scheme A stabilized numerical solution can be found by the following algorithm: Step 1. For elements in the interior of the domain choose 1n = u, i.e. the gradient direction is taken coincident with the velocity direction. If u = 0 then 1 n is taken coincident with the global x axis. In elements adjacent to a boundary choose 1n = n where n is the normal to the boundary. 0, 1D and 1DG using the expressions Compute 1 g; 1 D of bn and bg from Eq. (42).
105
that the solution is stable. This Solve for 1 /.Verify implies that there are not undershoots or overshoots in the numerical results with respect to the expected physical values. If the numerical solution is unstable, then go to step 2. Step 2. For all elements, compute at the element cen Then compute 2 g; 2 D 0, 2D and 2DG ter 2 n ¼ $1 /. using the expressions of bn and bg from Eqs. (37). Solve for 2 /. Step 3. Estimate the convergence of the process. We have chosen the following convergence norm: 1 k/k ¼ N/
max
" n X
i
/j i1 / j
2
#1=2 6e
ð43Þ
j¼1
where N is the total number of nodes in the mesh and /max is the maximum prescribed value at the Dirichlet boundary (if / max ¼ 0 then /max ¼ 1). In above steps the left upper indices denote the iteration number.In the examples shown in the next section e = 103 has
Fig. 13. Solution of the problem of Fig. 12 with an unstructured mesh of 432 four node bi-linear quadrilaterals.
106
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
been taken.If condition (43) is not satisfied, repeat steps 2 and 3 until convergence. Remark. For the advective–diffusive problems (i.e. s = 0) the expression of the balancing diffusion matrix in the interior elements for step 1 coincides with the standard (linear) SUPG form [30]. Remark. An alternative solution scheme is to use a time relaxation technique based in the solution of a pseudo transient problem with a forward Euler scheme and a diagonal mass matrix. 6. 1D numerical examples The examples presented in this section are solved in a 1D domain discretized with eight two-node linear elements of unit size. The examples are solved with the
same Dirichlet conditions /p1 ¼ 8 and /p9 ¼ 3 and two different values of c and w (c = 1, w = 20 and c = 10, w = 20). We note that for both cases the instability condition (bc > 0) is violated and, hence, the Galerkin solution should yield non-physical results. Figs. 2 and 3 show the numerical results obtained with the standard Galerkin method (b = 0) and using the element (be) and critical (bc) values of the stabilization parameter given by Eqs. (22) and (23), respectively. The exact analytical solution is also shown for comparison. Table 1 shows the nodal values of the results of the example of Fig. 3 for comparison with the 2D solutions presented in the next section. The following conclusions are drawn from the 1D results: 1. The Galerkin solution (b = 0) is unstable for both problems, as expected.
Fig. 14. Solution of problem of Fig. 12 with an structured mesh of 10 · 10 · 2 three node linear triangles.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
2. The solution obtained with the critical value bc is always stable, although it yields slightly overdiffusive results in both cases. 3. The results obtained with be are less diffusive and more accurate than those obtained with bc. The explanation is that the sign of the ratio S1/S2 is negative in the region close to the left end point of the domain. This naturally reduces the value of the stabilizing diffusion parameter b in Eq. (21) with respect to that of bc in Eq. (23) where the sign effect is not relevant. 4. The FIC algorithm provides a stabilized solution for Dirichlet problems when there is a negative streamwise gradient of the solution. This is an advantage versus SUPG, GLS and SGS methods using a single stabilization parameter which fail in some cases for these type of problems [12,13].
107
Above conclusions have been confirmed in the solution of a wider collection of 1D problems presented in [18].
7. 2D examples The analysis domain in the first two 2D examples presented is a square of size 8 units. The problems are solved first with relatively coarse meshes of 8 · 8 four node bi-linear square elements and 8 · 8 · 2 linear triangles. The boundary conditions for both examples are /p = 8 and /p = 3 at the boundaries x = 0 and x = 8, respectively and zero normal flux at y = 0 and y = 8. This reproduces the condition of the two 1D examples solved in the previous section. The first example is analyzed for
Fig. 15. Solution of problem of Fig. 12 with an unstructured mesh of 780 three node triangles.
108
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
u = [2, 0]T, k = 1 and s = 20 giving w = 20, cx = 1 and cy = 0 which corresponds to the first 1D example (Fig. 2). The correct solution for this problem has a boundary layer in the vecinity of the two sides at x = 0 and x = 8 where / is prescribed (Fig. 4). The numerical results obtained with the standard Galerkin solution are oscillatory as expected. The stabilized FIC formulation eliminates the oscillations and yields the correct physical solution. Good results are obtained for both meshes of linear rectangles and triangles (Figs. 4 and 5). Results labelled as FIC-1 and FIC-2 in the figures correspond to those obtained in the first and second iteration of the algorithm presented in Section 5.2, respectively. We note that the FIC-1 results agree precisely with those obtained in the 1D case for b = bc, whereas the FIC-2 results agree with the more accurate 1D values obtained with the element stabilization parameter be (see Fig. 2). The second example is similar to the first one with u = [20, 0]T, k = 1 and s = 20 giving w = 20, cx = 10
and cy = 0. These values correspond to the second 1D problem of the previous section (Fig. 3). The Galerkin solution is again oscillatory, whereas the FIC results are physically sound (Figs. 6 and 7). Once more the FIC-1 and FIC-2 results are in good agreement with the 1D values shown in Fig. 3 for bc and be, respectively for both meshes of square and triangular elements. The coincidence of the 1D and 2D results for this problem can be clearly seen in Table 1. Note that, similarly to the 1D case, the FIC-2 results are more accurate (less diffusive) than those obtained in the first iteration (FIC-1). This is due to the more precise evaluation of bn and bg in Eqs. (37) accounting for the correct sign of all the terms. Figs. 8–11 show results for the two 2D problems above described solved now with relatively coarse unstructured meshes of linear triangles and quadrilaterals. The effectiveness and accuracy of the FIC iterative scheme is again noticeable in all cases. Note the agreement of the FIC-2 results of Figs. 10 and 11 with
Fig. 16. Diffusive–absorptive problem over a square domain of size equal to 10 units. /p = 50 over x = 0 and y = 10 and /p = 100 over x = 10 and y = 0, u = [0, 0]T, k = 1, s = 20, w = 20, cx = 0 and cy = 0. Galerkin and FIC solutions obtained with structured meshes of four node quadrilaterals and linear triangles.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
the exact solution for the equivalent 1D problem of Fig. 3. Fig. 12 presents the solution of a similar problem where the values of / are prescribed at the four boundaries. The solution domain has now 10 units and the problem is solved first with a mesh of 10 · 10 four node square elements. Details of the physical parameters are given in Fig. 12. Excellent results are again obtained with the FIC scheme. Similar good results are obtained with a structured mesh of linear triangles (Fig. 13) as well as with non-structured meshes of linear quadrilateral and triangles (Figs. 14 and 15). The effectiveness of the FIC scheme for a diffusive– absorptive problem with Dirichlet boundary conditions is shown in Fig. 16. The results shown have been obtained with structured meshes of linear quadrilateral and triangles. Note that the four boundary layers are well captured in the first step of the iterative solution.
109
Similar good results have also been obtained with unstructured meshes not shown here. The final example is a standard benchmark problem of advection–diffusion where sharp layers appear at both the boundary and the interior of the domain. The problem is the advective–diffusive transport of / in a square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source terms (i.e. Q = 0 and s = 0). Fig. 17 displays the SUPG solution and FIC results obtained after two iterations using a structured mesh of 20 · 20 linear four node square elements. It is remarkable that the FIC results capture the sharp gradient zones at the boundaries where / is prescribed to zero and at the interior of the domain and eliminate all the spurious oscillations present in the SUPG method. Similar good results obtained with the FIC method for a wide range of advective–diffusive problems are presented in [30]. Recent applications of the FIC
Fig. 16 (continued)
110
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
Fig. 17. Square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 20 · 20 linear four node square elements.
method to incompressible fluid flow problems are reported in [37].
8. Conclusions The FIC–FEM formulation presented allows to obtain a stabilized and accurate solution for the advection–diffusion–absorption equation. For the 1D problem the formulation is equivalent to adding a nonlinear diffusion term to the standard discretized equations which is governed by a single stabilization parameter. The use of the constant critical value of the 1D stabilization parameter provides a stabilized numerical solution in a single step. A more accurate (less diffusive) solution can be obtained using the two step iterative scheme proposed. The equivalence of the FIC method with a nonlinear stabilizing diffusion term extends naturally to multidi-
mensional problems using structured and unstructured meshes. The key step is to express the governing equations of the FIC formulation in the principal curvature directions of the solution. The resulting FIC equation is equivalent to adding a nonlinear diffusion matrix to the infinitesimal governing equations. The solution process becomes nonlinear and a simple iterative algorithm has been presented. The approximation of the main principal curvature direction by that of the gradient vector simplifies the computations in the iterative scheme. Excellent results have been obtained for all the 2D problems solved in just two iterations for structured and nonstructured meshes. It is remarkable that, similarly to the 1D case, good stabilized results are obtained in the first iteration of the scheme proposed (FIC-1 results) and this may be sufficient for many practical cases. More accurate (less diffusive) results are obtained by performing a second iteration at a relatively small additional computational cost.
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
Acknowledgements The authors also thank Profs. C. Felippa and S.R. Idelsohn for many useful discussions. This work has been sponsored by the Ministerio de Educacio´n y Ciencia of Spain. Plan Nacional, Project numbers: DPI2001-2193, BIA2003-09078-C02-01, and DPI2004-07410-C03-02.
References [1] Zienkiewicz OC, Taylor RL. fifth ed. The finite element method, vol. 3. Butterworth-Heinemann; 2001. [2] Tezduyar TE, Park YJ. Discontinuity-capturing finite element formulations for nonlinear convection–diffusion–reaction equations. Comput Methods Appl Mech Engrg 1986;59:307–25. [3] Tezduyar TE, Park YJ, Deans HA. Finite element procedures for time-dependent convection–diffusion–reaction systems. Int J Numer Meth Fluids 1987;7:1013–33. [4] Idelsohn S, Nigro N, Storti M, Buscaglia G. A Petrov–Galerkin formulation for advection–reaction–diffusion problems. Comput. Methods Appl Mech Engrg 1996;136:27–46. [5] Codina R. Comparison of some finite element methods for solving the diffusion–convection–reaction equation. Comput Methods Appl Mech Engrg 1998;156:186–210. [6] Codina R. On stabilized finite element methods for linear systems of convection–diffusion–reaction equations. Comput Meth Appl Mech Engrg 2000;188:61–82. [7] Burman E, Ern A. Nonlinear diffusion and discrete maxium principle for stabilized Galerkin approximations of the convection–diffusion–reaction equation. Comput Meth Appl Mech Engrg 2002;191:3833–55. [8] Harari I, Hughes TJR. Finite element methods for the Helmholtz equation in an exterior domain: model problems. Comput Meth Appl Mech Engrg 1991;87:59–96. [9] Harari I, Hughes TJR. Galerkin/least squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains. Comput Meth Appl Mech Engrg 1992;98:411–54. [10] Harari I, Hughes TJR. Stabilized finite element method for steady advection–diffusion with production. Comput Meth Appl Mech Engrg 1994;115:165–91. [11] Harari I, Grosh K, Hughes TJR, Malhotra M, Pinsky PM, Stewart JR, et al. Recent development in finite element methods for structural acoustics. Arch Comput Mech Eng 1996;3(2-3): 131–309. [12] Hauke G, Garcia-Olivares A. Variational subgrid scale formulations for the advection–diffusion–reaction equation. Comput Methods Appl Mech Engrg 2000;190:6847–65. [13] Hauke G. A simple subgrid scale stabilized method for the advection–diffusion reaction equation. Comput. Methods Appl Mech Engrg 2002;191:2925–47. [14] Brezzi F, Hauke G, Marin LD, Sangalli S. Link-cutting bubbles for the stabilization of convection–diffusion–reaction problems. Math Models Methods Appl Sci. World Scientific Publishing Company; 2002. [15] Felippa CA, On˜ate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Comput Mech, in print, CM-050155.
111
[16] Idelsohn SR, Heinrich JC, On˜ate E. Petrov–Galerkin methods for the transient advective–diffusive equation with sharp gradients. Int J Numer Meth Engng 1996;39:1455–73. [17] Harari I. Stability of semidiscrete advection–diffusion in transient computation. In: Yao ZH, Yuan MW, Zhong WX, editors. Proceedings of 6th World Congress on Computational Mechanics, Beijing, September. Tsinghua Univ. Press-Springer; 2004. [18] On˜ate E, Miquel J, Hauke G. Stabilized formulation for the advection–diffusion–reaction equations using finite calculus and linear finite elements. Comput Methods Appl Mech Engrg, submitted for publication. [19] On˜ate E. Derivation of stabilized equations for advective–diffusive transport and fluid flow problems. Comput Methods Appl Mech Engrg 1998;151(1-2):233–67. [20] On˜ate E. Possibilities of finite calculus in computational mechanics. Int J Numer Meth Engng 2004;60:255–81. [21] On˜ate E, Manzan M. A general procedure for deriving stabilized space–time finite element methods for advective–diffusive problems. Int J Numer Meth Fluids 1999;31:203–21. [22] On˜ate E, Idelsohn SR, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics. Applications to convective transport and fluid flow. Int J Numer Meth Engng 1996;39: 3839–66. [23] On˜ate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comput Methods Appl Mech Engrg 2000;182(1-2): 355–70. [24] On˜ate E, Garcı´a J. A finite element method for fluid–structure interaction with surface waves using a finite calculus formulation. Comput Methods Appl Mech Engrg 2001;191(6–7): 635–60. [25] Idelsohn SR, On˜ate E, Del Pin F. A Lagrangian meshless finite element method applied to fluid–structure interaction problems. Comput Struct 2003;81:655–71. [26] On˜ate E, Idelsohn SR, Del Pin F, Aubry R. The particle finite element method. An overview (PFEM). Int J Comput Methods 2004;1(2):267–307. [27] On˜ate E, Garcı´a J, Idelsohn SR. Ship hydrodynamics. In: Hughes T, de Borst R, Stein E, editors. Encyclopedia of computational mechanics, Vol. 3. J. Wiley; 2004. p. 579–607. Chapter 18. [28] On˜ate E, Taylor RL, Zienkiewicz OC, Rojek J. A residual correction method based on finite calculus. Eng Comput 2003;20(5/6):629–58. [29] On˜ate E, Rojek J, Taylor RL, Zienkiewicz OC. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. Int J Numer Meth Engng 2004; 59(11):1473–500. [30] On˜ate E, Za´rate F, Idelsohn SR. Finite element formulation for convection–diffusion problems with sharp gradients using finite calculus. Comput Meth Appl Mech Engng, submitted for publication. [31] Hughes TJR, Mallet M. A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective–diffusive system. Comput Methods Appl Mech Engrg 1986;58:329–36. [32] Codina R. A discontinuity-capturing crosswind dissipation for the finite element solution of the convection–diffusion equation. Comput Methods Appl Mech Engrg 1993;110: 325– 42. [33] Hughes TJR, Mallet M, Mizukami A. A new finite element formulation for comptutational fluid dynamics: II Beyond SUPG. Comput Methods Appl Mech Engrg 1986;54:341–55.
112
E. On˜ate et al. / Computers & Fluids 36 (2007) 92–112
[34] Galea˜o AC, Dutra do Carmo EG. A consistent approximate upwind Petrov–Galerkin method for convection-dominated problems. Comput Methods Appl Mech Engrg 1988;68:83–95. [35] Tezduyar TE. Adaptive determination of the finite element stabilization parameters. Proceedings of the ECCOMAS Computational Fluid Dynamics Conference 2001, Swansea, Wales, UK, CD-Rom, 2001.
[36] Tezduyar TE. Computation of moving boundaries and interfaces and stabilization parameters. Int J Numer Methods Fluid 2003;43:555–75. [37] On˜ate E, Garcı´a J, Idelsohn S, Del Pin F. FIC formulation for finite element analysis of incompressible flows. Eulerian, Lagrangian and ALE approaches. Comput Methods Appl Mech Eng, in press.