- Email: [email protected]
A stabilized non-conforming finite element method for incompressible flow
Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899 www.elsevier.com/locate/cma
A stabilized non-conforming finite element method for incompressible flow Erik Burman a, Peter Hansbo a
b,*
Institut d’Analyse et de Calcul Scientifique, Ecole Polytechnique Fe´de´rale de Lausanne, CH-1015 Lausanne, Switzerland b Department of Applied Mechanics, Chalmers University of Technology, SE–41296 Go¨teborg, Sweden Received 31 January 2004; received in revised form 30 September 2004; accepted 23 November 2004
Abstract In this paper we extend the recently introduced edge stabilization method to the case of non-conforming finite element approximations of the linearized Navier–Stokes equation. To get stability also in the convective dominated regime we add a term giving L2-control of the jump in the gradient over element boundaries. An a priori error estimate that is uniform in the Reynolds number is proved and some numerical examples are presented. 2005 Elsevier B.V. All rights reserved. Keywords: Stabilized methods; Finite element; Crouzeix–Raviart; Non-conforming; Interior penalty; Incompressible flow
1. Introduction The solution of the Navier–Stokes equations for incompressible flow using finite element methods remains a challenging problem, in particular if the objective is to construct a method which remains robust and accurate for a wide range of Reynolds numbers. The discretization must assure not only satisfaction of the Babus˘ka–Brezzi condition but also stabilization of the convective terms and sufficient control of the incompressibility condition. Approximations using non-conforming Crouzeix–Raviart (CR) elements are attractive for the velocity approximation in combination with elementwise constant pressures, since they satisfy the Babus˘ka–Brezzi condition and have local conservation properties. This discretization was
*
Corresponding author. Fax: +31 772 3827. E-mail address: [email protected] (P. Hansbo).
0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.11.033
2882
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
proposed and analyzed in [14] and a stabilized version using the streamline diffusion stabilization was analyzed in [12]. In neither of these cases the high Reynolds number limit was treated. Moreover, in a recent paper, [4], the authors showed that a stabilized non-conforming finite element method using the Crouzeix– Raviart element remains uniformly stable in the vanishing viscosity limit for the generalized Stokes equation. However, the discretization of convection dominated problems using CR-elements are not stable without both a weak coupling on the element inflow boundary as in the discontinuous Galerkin method and streamline diffusion stabilization of the convective terms in the interior of the element (see [10,9,13]). In this respect, the element needs stabilization from both the streamline diffusion method and the discontinuous Galerkin method. Considering the fact that the method uses more degrees of freedom than the continuous Galerkin approximation this seems suboptimal. Moreover, the streamline diffusion stabilization has the drawback that it does not permit lumped mass for time stepping. In this paper we therefore propose to apply the recently introduced edge stabilization operator (see [2,3]) to the lowest order Crouzeix–Raviart element for the stabilization of the convective terms. We prove that this operator stabilizes exactly that part of the convective term which is not already included in the approximation space. In this sense this is the smallest perturbation needed to make the Crouzeix–Raviart element stable for convection– diffusion problems. This stabilization method has the advantage, as compared to other stabilized methods, that we may lump mass for efficient timestepping, we do not add any additional degrees of freedom, and we do not need any special structure of the mesh. For Oseens equation we prove an optimal a priori error estimate in the energy norm independent of the Reynolds number. Another attractive feature of the proposed stabilization is that, unlike SUPG, here the stabilization parameter is independent of the flow regime; we illustrate this by proving an L2 a priori error estimate for the velocities in the case of low local Reynolds number. Finally, we study the performance of the numerical scheme on some linear and nonlinear model cases.
2. A finite element method for the Oseen’s equation We consider, in X Rd with boundary oX, the problem of solving ru þ b ru þ rp 2lr eðuÞ ¼ f ;
in X;
r u ¼ 0; in X; u ¼ 0 on oX;
ð1Þ
d
where u, b 2 ½H 10 ðXÞ \ H 0 ðdiv; XÞ , b 2 W1,1(X), p 2 L20 ðXÞ, f is a given source term, r and l are bounded positive functions. By H0(div; X) we denote the functions in [L2(X)]d such that $ Æ u = 0, and by L20 ðXÞ the functions in L2(X) with zero mean value. In the following, we let (Æ, Æ) denote the L2-scalar product with the corresponding norm k Æ k = (Æ, Æ)1/2; the Hs(X) norm will be denoted by k Æ ks,X. d The weak form of (1) is to find ðu; pÞ 2 ½H 10 ðXÞ L20 ðXÞ such that (
aðu; vÞ þ bðp; vÞ ¼ ðf ; vÞ; bðq; uÞ ¼ 0;
8ðv; qÞ 2
½H 10 ðXÞd
ð2Þ
L20 ðXÞ;
where aðu; vÞ :¼ ðru; vÞ þ ðb ru; vÞ þ 2ðleðuÞ; eðvÞÞ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2883
and bðp; vÞ :¼ ðp; r vÞ. The well posedness of the above problem follows by the Lax-Milgram Lemma applied in the space ½H 10 ðXÞd \ H 0 ðdiv; XÞ. The finite element method consists of seeking a piecewise polynomial approximation uh 2 Vh, ph 2 Qh. We let Vh denote the space of the lowest order non-conforming Crouzeix–Raviart elements. Let Th S denote a shape regular triangulation of the domain X, EðKÞ the set of all faces of an element K 2 Th , E :¼ K2Th EðKÞ the set of all faces in Th , EoX :¼ fe 2 E : e oXg, and E0 :¼ E n EoX the set of the boundary and inner faces respectively. For a given piecewise continuous function u, the jump [u] and the average {u} on a face e 2 E are defined by 8 < limþ ðuðx snÞ uðx þ snÞÞ if e 6 oX; s!0 ½uðxÞ :¼ : lim ¼ uðx snÞ if e oX; s!0þ 81 < 2 limþ ðuðx þ snÞ þ uðx snÞÞ if e 6 oX; s!0 fugðxÞ :¼ 1 : 2 lim 2uðx snÞ if e oX; þ s!0
where n is a normal unit vector on e and x 2 e. If e oX we choose the orientation of n to be outward with respect to X otherwise n has an arbitrary but fixed orientation. For the non-conforming finite element functions, continuity across edges e will only be enforced with respect to Z je ðvh Þ :¼ ½vh ds. e
Using this definition our finite element space may be defined as V h :¼ fvh 2 ½L2 ðXÞd : vh jK 2 ½P 1 ðKÞd
8K 2 Th ;
je ðvh Þ ¼ 0; 8e 2 Eg.
Moreover, we introduce the space of piecewise constants with mean value zero, Qh :¼ fqh 2 L20 ðXÞ : qh jK 2 P 0 ðKÞg and the subspace Wh of Vh such that W h :¼ fwh 2 V h : ðr wh ; qh Þh ¼ 0 8qh 2 Qh g P where ðr wh ; qh Þh ¼ K ðr wh ; qh ÞK . Since the above spaces are H1-non-conforming we introduce the broken norm equivalent of the L2-norm X 2 2 kukh ¼ kukK K
and the broken H1-seminorm X 2 juj2h ¼ juj1;K . K
The local mesh size is defined by hK :¼ max hoK ; K
and we will assume that hK/hoK < C where C is a fixed constant. We will use C and c as generic constants taking different values at different occasions. To indicate their provenance or main dependence, a subscript may be added, e.g., cl. We introduce the interpolation operator rhu : [H1(X)]d ! Vh defined by
2884
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
rh uðxe Þ ¼
1 jej
Z u ds; e
where xe is the midpoint of the edge e. The L2-projections onto the spaces are also required for the analysis. Let p0,h : L2(K) ! P0(K), pd0;h : ½L2 ðKÞd ! ½P 0 ðKÞd denote the L2-projection onto the constant functions on K and p1,h : [L2(K)]d ! Vh the standard L2-projection onto the finite element space. For the above defined interpolation operator and projections we need some approximation and stability properties. These, and some inverse inequalities are collected in the following Lemmas. Lemma 1. For the interpolation operator rh there holds, if u 2 [H2(X)]d then krh u ukh þ hjrh u ujh 6 C r h2 .
ð3Þ
Moreover, if $ Æ u = 0 then rhu 2 Wh. Proof. The proof of the interpolation estimate is given by Crouzeix and Raviart [7]. The second claim is immediate noting that Z Z r rh u dx ¼ r u dx ¼ 0 K
K
by the definition of the interpolant.
h
Lemma 2. For the L2-projection the following H1 stability holds, jp1;h ujh 6 C s kuk1;X . For the error analysis, we shall use the following trace inequality. Lemma 3. For v 2 H1(K) there holds 2 2 kvk2oK 6 C t ðh1 K kvkK þ hK kvk1;K Þ
8v 2 H 1 ðKÞ;
ð4Þ
where Ct is a constant independent of hK. We also need the following local inverse inequality. Lemma 4. Let uh 2 Vh where Vh is defined on a shape regular mesh then kruh kK 6 h1 K C i kuh kK ; with Ci independent of K. Proof. For proofs of Lemmas 2–4, see, respectively, Carstensen [5], Thome´e [15], and Ciarlet [6]. Our finite element method reads: find uh 2 Vh such that ah ðuh ; vh Þ þ bh ðph ; vh Þ þ ju ðuh ; vh Þ ¼ ðf ; vh Þ bh ðqh ; uh Þ ¼ 0 8ðvh ; qh Þ 2 V h Qh ;
h
ð5Þ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2885
where ah ðuh ; vh Þ ¼ ðruh þ b ruh ; vh Þh
1X hb n½uh ; fvh gioK þ ð2leðuh Þ; eðvh ÞÞh ; 2 K
bh ðph ; vh Þ ¼ ðph ; r vh Þh ;
ð6Þ ð7Þ
and the jump terms take the form XZ cb 2 ju ðuh ; vÞ ¼ hoK bh ½ruh bh ½rvds oKnoX jbh j K XZ ca ðl þ jb njhoK ÞhoK ½t ruh ½t rvds þ K
þ
X K
oK
Z cb
hoK ½ðt ruh Þ n½ðt rvÞ nds.
ð8Þ
oK
Here bh is the interpolant of b on Wh and t is a unit vector perpendicular to n. The gradient jump term serves three purposes. It stabilizes the convective terms (the first sum in (8)), it assures that Korns inequality is satisfied (the cal part of the second sum in (8)), it gives additional control of the divergence inconsistency error (the third sum in (8)). The last two properties can be obtained by introducing a lower order penalizing term (see [4]), but the use of the jump of the gradient has the advantage of allowing for one point quadrature in the implementation and from the point of view of analysis it is practical. In the case of three space dimensions the tangent vector should be replaced by the tangent tensor $uh · n. In the following we will for simplicity only consider the two dimensional case for the tangent vectors. Remark 5. The analysis below holds with only minor modifications if b is replaced by bh also in the convective term. This is convenient when timestepping the Navier–Stokes equation: b may be taken as the solution uh of the previous timestep. In the analysis we will not distinguish between the different stabilization parameters, they will all be denoted c. We introduce the following shorthand notation A½ðuh ; ph Þ; ðwh ; qh Þ ¼ ah ðuh ; vh Þ þ bh ðph ; vh Þ bh ðqh ; uh Þ. To simplify the analysis we will assume that the exact solution (u, p) belongs to [H2(X)]d · H1(X); it then follows that the formulation (5) enjoys the following consistency property. Lemma 6. For u 2 H2(X) and p 2 H1(X) there holds A½ðu uh ; p ph Þ; ðvh ; qh Þ þ ju ðu uh ; vh Þ ¼ Rðu; p; vh Þ
ð9Þ
for all (vh, qh) 2 Vh · Qh. Where the consistency error due to the non-conforming approximation is given by 1X 1X Rðu; p; vh Þ ¼ h2leðuÞ n; ½vh ioK þ hp; ½vh nioK . 2 K 2 K Proof. This is an immediate consequence of the regularity hypothesis: if u 2 H2(X) then the trace of $u is well defined and hence j(u, vh) = 0. The consistency error is obtained by integration by parts, X X 1 ð2leðuÞ; eðvh ÞÞ ¼ ð2lr eðuÞ; vh Þ þ h2leðuÞ n; ½vh ioK 2 K K
2886
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
and
X
ðp; r vh ÞK ¼ ðrp; vh Þ
K
1X hp; ½vh nioK . 2 K
2.1. Preliminary lemmas In this section we will prove some preliminary results that will facilitate the analysis. The main result of this section is Lemma 7 where we show that that the jump term (8) controls the difference between the convective derivative and its quasi-interpolant on the finite element space. This lemma is the key ingredient to derive error bounds that are independent of the Peclet number. Then we prove the coercivity of the bilinear form. The triple norm that we will use is given by jjjðwh ; qh Þjjj2 ¼ kr1=2 wh k2h þ jl1=2 wh j2h þ kr wh k2h þ ju ðwh ; wh Þ þ cp kqh k2 ; where cp is a constant depending on the problem parameters r, l, b to be specified later. Note that although for the Crouzeix–Raviart element $ Æ uh = 0 on each triangle the formulation must include a stabilization of the jump of the normal velocity due to the H(div, X) consistency error. Therefore the triple norm must be d chosen as a discrete norm on ðH ðdiv; XÞ \ l1=2 ½H 10 ðXÞ Þ L20 ðXÞ. The triple norm is dominated, at low Rey1 nolds numbers, by the H (X) contribution, and at high Reynolds numbers by the part of the jump term controlling the inconsistency in the divergence. This latter term prohibits the decoupling of the velocities and the pressure and the order of the estimate can be no better that the approximation properties of the pressure space. We introduce the space of functions that are piecewise linear on each element d
Y h ¼ ½fy 2 L2 ðXÞ : yjK 2 P 1 ðKÞg . h : Y h ! V h . Let xi be the midpoint of the We now introduce a quasi-interpolant based on local averages p face shared by element K and element K 0 then fuðxi Þg for xi an interior node; h uðxi Þ ¼ p uðxi Þ for xi a boundary node. In the following Lemma we prove that the projection error is bounded by the jumps in the gradient. Lemma 7. Let bh 2 Vh and wh 2 Vh then 2
h ðbh rwh ÞÞk 6 jb ðwh ; wh Þ; kh1=2 ðbh rwh p where jb(wh, wh) is given by X Z jb ðwh ; wh Þ ¼ cb K
hK hoK ? ðbh ½rwh Þ2 ds
oKnoX
with hoK ? denoting the triangle size perpendicular to the side on oK and cb is a parameter depending only on the number of space dimensions. h ðbh rwh Þ makes sense. Now consider any Proof. First note that (Vh Æ $)Vh Yh so that the projection p d triangle K and note that using the Crouzeix–Raviart basis functions fui gi¼1 , with ui associated with node xi we may write bh rwh jK ¼
dþ1 X i¼1
bh rwh ðxi Þui ðxÞ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2887
and h ðbh rwh ÞjK ¼ p
dþ1 X
fbh rwh ðxi Þgui ðxÞ.
i¼1
Taking the difference of the two functions in a nodal point yields bh rwh ðxi ÞjK fbh rwh ðxi Þg ¼ bh ðxi Þ ½rwh ðxi Þ.
Note that if the node is on the boundary, the right hand side is zero. It follows that for any K 2 Th such that oK \ oX = ; kh
1=2
h bh ðbh rwh p
rwh Þk2K
¼
Z
dþ1 X ðbh rwh ðxi Þ fbh rwh gÞui
hK K
dx
i¼1
Z
dþ1 X
hK
6
!2
K
!2 bh ½rwh ui
dx.
ð10Þ
i¼1
We now evaluate the integral using nodal point quadrature and note that since the nodes of the Crouzeix– Raviart element are on the midpoints of the element sides this is exact for second degree polynomials in two space dimensions, hence 1=2
2
h bh rwh ÞkK ¼ khK ðbh rwh p
3 X k¼1
hK
3 X mK moK k 2 2 bh ðxk Þ ½ruh ¼ bh ðxk Þ ½ruh ðxk Þ ; hK h? oK 3 6 k¼1
ð11Þ
R where moK k ¼ oK k dx with oKk the face associated with quadrature point k. In three space dimensions the midpoints of the faces has to be supplemented with the six midpoints of the edges of the tetrahedron to yield an exact quadrature formula (with weights 1/15 for the midpoints of the faces and 3/20 for the midpoints of the edges). One may then easily show that 1=2
2
h bh rwh ÞkK 6 khK ðbh rwh p
3 5X moK k 2 hK h? ðbh ðxk Þ ½ruh ðxk ÞÞ . oK 6 k¼1 6
ð12Þ
It follows, using the Simpson quadrature formula in two dimensions and a quadrature taking the midpoint d , that of the face and the corner-points in three dimensions and noting that the weight for the midpoint is dþ1 Z 2 h ðbh rwh ÞÞk2K 6 cd kh1=2 ðbh rwh p hK hoK ? ðbh ½ruh Þ ds; ð13Þ oK
where cd = 1/4 in two dimensions and cd = 10/9 in three dimensions (the subscript d is here referring to ‘‘dimension’’). We conclude by taking the sum over all triangles K 2 T noting that all boundary contributions vanishes thanks to the definitions of the quasi-interpolant. Remark 8. A consequence of the above proof is that the jump term edge integral may be evaluated using midpoint quadrature on the faces. In fact the first part of the jump operator given in (8) can be substituted by the discrete operator given by (11) in two space dimensions and by (12) in three to get optimal values of the stabilization constants. The integral formulation however still remains practical from a theoretical viewpoint since it is consistent for H2-regular solutions. As was pointed out in the previous section we only stabilize using the jumps in the gradient. However we need to establish a result showing the equivalence between the jumps in the solution and the jumps in the tangential gradient.
2888
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Lemma 9. For the interior penalty term (8) there holds X 2 kl1=2 h1=2 ½wh koK 6 cju ðwh ; wh Þ; K
X
2
kh1=2 ½wh nkoK 6 cju ðwh ; wh Þ;
K
and ju ðwh ; wh Þ 6 cc jwh jh for all wh 2 Vh. Proof. By the midpoint continuity of the Crouzeix–Raviart element we note that we may write, with n a coordinate along e with midpoint ni, ½wh ðnÞje ¼ ½t rwh je ðn ni Þ and ½wh ðnÞ nje ¼ ½ðt rwh Þ nje ðn ni Þ. Hence we have Z Z Z 1 2 2 2 ½wh ðnÞ dn ¼ ð½t rwh ðn ni ÞÞ dn ¼ h2 ½t rwh dn; 12 e e e e which proves the first claim. The proof of the second claim is equivalent. The last claim finally is an immediate consequence of the trace inequality (4). h Lemma 10. For the consistency error the following upper bound holds jRðu; p; vh Þj 6 ðcl hkuk2;X þ Chkpk1;X Þju ðvh ; vh Þ. Proof. Using the zero mean value property of the Crouzeix–Raviart space followed by the Cauchy– Schwarz inequality we obtain
Rðu; p; vh Þ 6 c
X
!1=2 kh
1=2
ð2lÞ
1=2
ðeðuÞ n
pd0;h eðuÞ
2 nÞkoK
K
þc
X
!1=2 1=2 1=2
kl
h
2 ½vh koK
K
X
!1=2 hK kp
p0;h pk2oK
K
X
!1=2 h1 K k½vh
nk2oK
ð14Þ
.
K
The claim now follows using the trace inequality (4), standard interpolation and Lemma 9.
h
Let us now investigate the coercivity properties of the discretization of the convective terms. Lemma 11. For the convective terms there holds 1X ðb ruh ; uh Þh ¼ hb n½uh ; fuh gi. 2 K
ð15Þ
Proof. The result follows by integrating by parts elementwise in the left hand side and using the equality [xy] = [x]{y} + {x}[y]. The integration by parts yields
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Z X X 1X ðb ruh ; uh ÞK ¼ ½b nuh uh ds ðuh ; b ruh ÞK . 2 K oK K K
2889
ð16Þ
We rewrite the edge term in the following fashion ½b nuh uh ¼ 2b n½uh fuh g
ð17Þ h
and the proof is completed using (17) in (16).
To prove the coercivity of our operator we also need the following discrete Korns inequality. Lemma 12. For wh 2 Vh there holds C 2K jl1=2 wh j2h 6 kð2lÞ1=2 eðwh Þk2h þ ju ðwh ; wh Þ. Proof. See Brenner [1].
h
The coercivity of our formulation is an immediate consequence of Lemmas 11 and 12. Lemma 13. For all (wh, qh) 2 Vh · Qh there holds C 2K jl1=2 wh j2h þ kr1=2 wh k2 6 A½ðwh ; qh Þ; ðwh ; qh Þ. Proof. First of all notice that the terms bh(uh, ph) cancel. We may write 2
A½ðwh ; qh Þ; ðwh ; qh Þ ¼ kr1=2 wh kh þ kð2lÞ
1=2
2
eðwh Þkh þ ðb rwh ; wh Þh
1X hb n½wh ; fwh gi. 2 K
Using Lemma 11 for the convective term we get ðb rwh ; wh Þ
1X hb n½wh ; fwh gi ¼ 0. 2 K
The proof is then completed by applying Lemma 12.
h
3. Stability In this section we will prove an inf–sup condition for our discretization vital for the convergence analysis. Theorem 14 (Stability). If (uh, ph) 2 Vh · Qh then there holds cis jjjðuh ; ph Þjjj 6
sup ðwh ;qh Þ2V h Qh
A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ ; jjjðwh ; qh Þjjj
ð18Þ
where the constant cis depends only on the parameters l, r, b, c and remains bounded from below when l ! 0. Proof. We prove the above inf–sup condition in two steps. First we will prove that there exists (wh, qh) 2 Vh · Qh such that 2
jjjðuh ; ph Þjjj 6 A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ
ð19Þ
2890
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
and then we conclude by proving that jjjðwh ; qh Þjjj 6 Cjjjðuh ; ph Þjjj.
For the first step we note that by Lemma 13 we have choosing wh = uh, qh = ph 2
2
C 2K jl1=2 uh jh þ kr1=2 uh k þ ju ðuh ; uh Þ 6 A½ðuh ; ph Þ; ðuh ; ph Þ þ ju ðuh ; uh Þ.
ð20Þ
2
To control the L -norm of the pressure we note, following [8], that by the surjectivity of the divergence d operator there exists a function vp 2 ½H 10 ðXÞ such that $ Æ vp = ph and jvpj1,X 6 ckphk. Therefore we now choose wh = rhvp and qh = $ Æ uh. By the stability of the quasi-interpolation operator rh we have krh vp k þ jrh vp jh 6 ckph k.
ð21Þ
Moreover, using the definition of the quasi-interpolant rhvp we have X X 2 kph k ¼ ðph ; r vp Þ ¼ h½ph ; vp nie ¼ h½ph ; rh vp nie ¼ ðph ; r rh vp Þh . e
ð22Þ
e
As a consequence of (22) we may write 2
2
A½ðuh ; ph Þ; ðrh vp ; r uh Þ þ ju ðuh ; rh vp Þ ¼ kph k þ kr uh k þ ðruh ; rh vp Þ þ ðlruh ; rrh vp Þh 1X þ ðb ruh ; rh vp Þ hb n½uh ; frh vp gioK þ ju ðuh ; rh vp Þ. 2 K ð23Þ Using Cauchy–Schwarz inequality, Youngs inequality, and the stability (21), we readily deduce 2
2
ðruh ; rh vp Þ P cr kr1=2 uh k 18kph k ; ð2leðuh Þ; eðrh vp ÞÞh P cl jl
1=2
uh j2h
ð24Þ
1 kph k2 8
ð25Þ
and by applying the third inequality of Lemma 9 ju ðuh ; rh vp Þ P ju ðuh ; uh Þ
1=2
1=2
ju ðrh vp ; rh vp Þ
2
P cc ju ðuh ; uh Þ 18kph k .
ð26Þ
It now remains to bound the convective term. An integration by parts yields ðb ruh ; rh vp Þ
1X 1X hb n½uh ; frh vp gioK ¼ ðuh ; b rrh vp Þ hb nfuh g; ½rh vp ioK . 2 K 2 K
ð27Þ
For the first term we clearly have ðuh ; b rrh vp Þ P 2ckbkL1 ðXÞ kuh k kph k.
ð28Þ
For the second we use the H1-regularity of vp, the trace inequality (4) and the local inverse inequality to obtain X 1X hb nfuh g; ½rh vp vp ioK P 2 kbkL1 ðKÞ kuh koK krh vp vp koK 2 K K X 2 2 1=2 1=2 P 2 kbkL1 ðKÞ ðh1 hK kvp k1;K K kuh kK þ hK kuh k1;K Þ K
P 2C i ckbkL1 ðXÞ kuh kkph k.
ð29Þ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2891
Combining (28) and (29) we obtain 1X 1 ðb ruh ; rh vp Þ hb n½uh ; rh vp ioK P cb kuh k2X kph k2 . 2 K 8
ð30Þ
Using now (24)–(26) and (30) to find a lower bound for the expression (23) we get 2
2
2
2
A½ðuh ; ph Þ; ðrh vp ; r uh Þ þ ju ðuh ; rh vp Þ ¼ 12kph k þ kr uh k C lrbc ðkr1=2 uh k þ jl1=2 uh jh þ ju ðuh ; uh ÞÞ; ð31Þ where C lrbc
! cb cl ¼ max cr ; ; 2 ; cc . r CK
Combining (20) and (31) we conclude that (19) holds for the choice wh = uh + (2Clrbc)1rhvp, qh = ph + $ Æ uh. More precisely we have, with cp = (2Clrbc)1, 1 jjjðuh ; ph Þjjj2 2
6 A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ.
ð32Þ
It remains to prove that jjjðwh ; qh Þjjj 6 Cjjjðuh ; ph Þjjj; this is obtained simply by noting that 2
2
2
2
2
jl1=2 wh jh 6 jl1=2 uh jh þ jl1=2 cp rh vp jh 6 jl1=2 uh jh þ lc2p ckph k .
ð33Þ
Proceeding in the same fashion for the other terms in the triple norm yields jjjðwh ; qh Þjjj2 ¼ jl1=2 wh j2h þ kr1=2 wh k2 þ ju ðwh ; wh Þ þ kr wh k2 þ kqh k2 6 jl1=2 uh j2h þ kr1=2 uh k2 þ jðuh ; uh Þ þ 2kr uh k2 þ Ckph k2 6 Cjjjðuh ; ph Þjjj2
ð34Þ
and we conclude that 2
cis jjjðuh ; ph Þjjj jjjðwh ; qh Þjjj 6 jjjðuh ; ph Þjjj 6 A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ. Clearly the constant cis is independent of h and furthermore it does not vanish for vanishing l.
4. Error estimates We will now proceed to derive a priori error estimates in the triple norm. The estimate takes the form jjjðu uh ; p ph Þjjj 6 Ch. The energy norm estimate is independent of the Peclet number, indicating that the proposed method should be stable for a wide range of Reynolds numbers when applied to the full Navier–Stokes equations. Lemma 15 (Approximation). Consider the projection (p1,hu, p0,hp) of the exact solution d
ðu; pÞ 2 ½H 2 ðXÞ H 1 ðXÞ onto the finite element space Vh · Qh. For the projection error there holds jjjðp1;h u u; p0;h p pÞjjj 6 Ch; where C 6 c(1 + l1/2 + jbj1/2h1/2 + r1/2h + c1/2)kuk2,X + cpkpk1,X.
2892
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Proof. By the optimal approximation property of the L2-projection we have kr1=2 ðp1;h u uÞk 6 cr h2 and kp0;h p pk 6 Ch. We now consider rhu p1,hu. Noting that rhu p1,hu = p1,h(rhu u) we obtain using the H1-stability of the L2-projection for the Crouzeix–Raviart element on locally quasi-uniform meshes (see [5]), krðrh u p1;h uÞkh 6 C s krðrh u uÞkh 6 C s C r h
ð35Þ
and we conclude that kl1=2 rðu p1;h uÞkh 6 cl h and kr ðu p1;h uÞkh 6 Ch. Finally we estimate the penalizing term ju(p1,hu u, p1,hu u). We have that XZ XZ 2 2 ½rðp1;h u uÞ ds 6 2 ðrðp1;h u uÞÞ ds oK oK K K X 2 2 ðh1 62 K krðp1;h u uÞkK þ hK kuk2;K Þ 6 Chkuk2;X ; K
where we have used the trace inequality (4) in the second inequality. We thus conclude that ju ðp1;h u u; p1;h u uÞ
1=2
6 cc1=2 ð1 þ l1=2 þ jbj
1=2 1=2
h
Þhkuk2;X .
Lemma 16 (Continuity). Let g = p1,hu u and f = p0,hp p be the projection error of the velocity and the pressure respectively, then there holds A½ðg; fÞ; ðwh ; qh Þ þ ju ðg; wh Þ Rðu; p; wh Þ 6 Cðjjjðg; fÞjjj þ ðcl þ cb h1=2 Þhkuk2;X þ Chkpk1;X Þjjjðwh ; qh Þjjj. Proof. Clearly we have using Cauchy–Schwarz inequality ðrg; wh Þ 6 Ckr1=2 gk jjjðwh ; 0Þjjj 6 cr jjjðg; fÞjjj jjjðwh ; qh Þjjj; ð2leðgÞ; eðwh ÞÞ 6 Cjl1=2 gjh jjjðwh ; 0Þjjj 6 cl jjjðg; fÞjjj jjjðwh ; qh Þjjj; ju ðg; wh Þ 6 jjjðg; fÞjjj jjjðwh ; qh Þjjj. We consider now the terms expressing the pressure velocity coupling. By the orthogonality of the L2-projection p0,h we have bðf; wh Þ ¼ ðp0;h p p; r wh Þh ¼ 0.
Using once again the Cauchy–Schwartz inequality we readily obtain bðqh ; gÞ ¼ ðqh ; r gÞh 6 jjjðwh ; qh Þjjj jjjðg; fÞjjj. It remains to treat the convective term and the non-consistency term. Let us first consider the convective term, an integration by parts yield together with the addition and substraction of bh yields 1X hb n½g; fwh gi ¼ ðg; ðb bh Þ rwh Þ ðg; bh rwh Þ ðb rg; wh Þ 2 K 1X hb n; fgg½wh i ¼ I þ II þ III. ð36Þ þ 2 K
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2893
The first term is controlled using a local inverse inequality and a local approximation result for b bh. !1=2 X 2 2 2 I 6 kgk kbkW 1;1 ðKÞ hK krwh k 6 cb kgkC i kwh k 6 cb h2 kuk2;X jjjðwh ; qh Þjjj. K
Using Lemma 7 we have for the term II h ðbh rwh ÞÞÞ 6 kh1=2 gk kh1=2 ðbh rwh p h ðb rwh ÞÞk II ¼ ðg; ðbh rwh p 6 c1=2 kh1=2 gkju ðwh ; wh Þ
1=2
6 cc h3=2 kuk2;X jjjðwh ; qh Þjjj.
ð37Þ
For the term III we obtain using the trace inequality (4) and the approximation properties of the L2projection. !1=2 X 1=2 1=2 III 6 hjb nj ½wh ; ½wh ioK kb nkL1 ðKÞ kfggkoK K
6 ju ðwh ; wh Þ1=2
X
1=2
kb nkL1 ðKÞ kfggkoK 6 cb h3=2 kuk2;X jjjðwh ; qh Þjjj.
K
Only the non-consistent residual remains to be bounded and we conclude the proof by applying Lemma 10. Theorem 17. Let (u, p) 2 [H2(X)]d · H1(X) be the solution of (1) and let (uh, qh) 2 Vh · Qh be the finite element solution of (5). Then there holds jjjðu uh ; p ph Þjjj 6 ðcr h þ cl þ cb h1=2 þ cc Þhkuk2;X þ Chkph k1;X .
ð38Þ
Proof. We will consider the discrete error ehu ¼ p1;h u uh and ehp ¼ p0;h p ph since using Lemma 15 we have jjjðu uh ; p ph Þjjj 6 jjjðu p1;h u; p p0;h pÞjjj þ jjjðp1;h u uh ; p0;h p ph Þjj 6 Ch þ jjjðp1;h u uh ; p0;h p ph Þjjj; where C is the constant given in Lemma 15. By Lemma 14 we have cis jjjðehu ; ehp Þjjj 6
A½ðehu ; ehp Þ; ðwh ; qh Þ þ ju ðehu ; wh Þ . jjjðwh ; qh Þjjj
ð39Þ
Using now Galerkin orthogonality we may write cis jjjðehu ; ehp Þjjj 6
A½ðg; fÞ; ðwh ; qh Þ þ ju ðg; wh Þ Rðu; p; wh Þ ; jjjðwh ; qh Þjjj
where g = p1,hu u and f = p0,hp p. We conclude the proof by applying Lemmas 16 and 15.
ð40Þ h
Remark 18. Since the estimate is only first order, due to the low order approximation of the pressure and the inclusion of the divergence in the triple norm, the virtues of the streamline stabilization are not obvious. We could in fact proceed with an inverse inequality in term II of (36) and still have the same formal convergence order of the triple norm. It is however known that this would destroy stability of the velocities for problems with important gradients. This is illustrated in the numerical section. As we mentioned in Section 1 the interior penalty method is independent of the Reynolds number. Of course, for low local Reynolds number the numerical scheme is stable without stabilization (except for the Korns inequality), so that the stabilization may be eliminated. We will however show that this is unnecessary
2894
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
for our discretization by proving that even when keeping the stabilizing terms our discretization has optimal L2-convergence of the velocities in the local low Reynolds number regime in spite of the fact that cb is independent of l. Consider the dual continuous problem of seeking / 2 ½H 10 ðXÞ and r 2 L20 ðXÞ such that r/ b r/ 2lr eð/Þ þ rr ¼ e in X; r/¼0
ð41Þ
in X;
where e :¼ u uh, and assume that we have the regularity estimate k/kH 2 ðXÞ þ krkH 1 ðXÞ 6 kek.
ð42Þ
We recall that since uh 2 Wh and rh/h 2 Wh we have $ Æ e = 0 and $ Æ (/ rh/) = 0 elementwise. Multiply the first line of (41) by e and integrate by parts to obtain X X h2ln eð/Þ; ½eioK þ hr; ½n eioK ; kek2 ¼ ah ðe; /Þ þ ju ðe; /Þ K
K
where we have used the partial integration X X X 1X ðe; b r/Þ ¼ ðb re; /ÞK hb ne; /ioK ¼ ðb re; /ÞK hb n½e; f/gioK 2 K K K K and the divergence free property of the error and of /. Using Galerkin orthogonality, the divergence free property of the interpolant and the zero mean value property of the Crouzeix–Raviart element, we obtain X 2 kek ¼ ah ðe; / rh /Þ þ ju ðe; / rh /Þ h2lðn eð/Þ pd0;h n eð/ÞÞ; ½eioK
X
K
h2lðn eð/Þ
pd0;h n
eð/ÞÞ; ½/ rh /ioK
K
þ
X
hr p0;h r; ½n eioK þ
K
X hp p0;h p; ½n ð/ rh /ÞioK K
6 jjjðe; 0Þjjj jjjð/ rh /; 0Þjjj þ cb ju ðe; eÞ
1=2
X K
þ
X
hK k2l1=2 ðn eð/Þ pd0;h n
K
þ
X
þ
þ
!1=2 jjjðe; 0Þjjj
!1=2 hK kr p0;h rkL2 ðoKÞ
jjjðe; 0Þjjj !1=2
hK k2l1=2 ðn eð/Þ pd0;h n
K
X
k/
2 eð/ÞÞkL2 ðoKÞ
K
X
!1=2 2 rh /koK
2 eð/ÞÞkL2 ðoKÞ
jjjð/ rh /; 0Þjjj
!1=2 hK kp p0;h pkL2 ðoKÞ
jjjð/ rh /; 0Þjjj.
K
Using Lemma 15, the trace inequality (4) and error estimates for rh and for piecewise constant interpolation, we arrive at 2
kek 6 Chðjjjðe; 0Þjjj þ hkukH 2 ðXÞ þ hkpkH 1 ðXÞ Þðk/kH 2 ðXÞ þ krkH 1 ðXÞ Þ and thus we have
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2895
Theorem 19. Under the regularity assumption (42), the L2-error in the velocities can be estimated as kek 6 Ch2 ðkukH 2 ðXÞ þ kpkH 1 ðXÞ Þ.
ð43Þ
5. Numerical results 5.1. Convergence study in the case of small viscosity Let k = (l1 (l2 + 16p2)1/2)/2. Then the exact solution to (1) is given by u1 ðx1 ; x2 Þ ¼ 1 ekx1 cos 2px2 ; kx1 k e sin 2px2 ; u2 ðx1 ; x2 Þ ¼ 2p 1 p ¼ e2kx þ C; 2 cf. [11], with b = u, r = 0 and a right hand side matching the exact solution. In our examples, we also chose C to give zero mean pressure. We solved this problem approximatively on X = (1/2, 3/2) · (0, 2), using stability parameters ca = cb = 1/100 and cb = 1/4. Concerning the size of the stabilizing parameters, our numerical experience is that they should be of approximately this size; in particular the scheme becomes excessively diffusive with parameters around or above 1. The parameters in this method should thus be chosen smaller than, e.g., in the SUPG method. In Fig. 1 we show the convergence for l = 103, and in Fig. 2 for l = 105. Note that the absolute value of the pressure decreases linearly with the inverse of l in L2, which is why the absolute error in pressure is smaller in Fig. 2. 5.2. Stability in the Navier–Stokes case We show the influence of the different stabilizing terms for the lid-driven cavity flow in X = (0, 1) · (0, 1) and with l = 103. In Fig. 3 we show the numerical solution of Navier–Stokes equations (with b = u) after
-4 Velocity Pressure
log (L2 -error)
-4.5 -5 -5.5 -6 -6.5 -7 -4.5
-4
-3.5
-3
-2.5
-2
-1.5
log (h)
Fig. 1. Convergence for l = 103.
-1
-0.5
0
2896
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899 -3.5 -4 -4.5
log (L2-error)
-5
Velocity Pressure
-5.5 -6 -6.5 -7 -7.5 -8 -8.5 -5
-4
-3
-2
-1
0
1
log (h)
Fig. 2. Convergence for l = 105.
Fig. 3. Approximate solution of the velocities, (left) with fluxes only and (right) with fluxes plus jump of the convective derivative.
15 fixed point iterations, we present the solution using only upwind fluxes, as well as the fully stabilized solution. Note that we do not get a wiggle free solution without the jump in the convective derivative. On the other hand, in Fig. 4 we show the solution obtained using only the jump in convective derivative as stabilization. This solution is markedly less diffusive, yet still retains stability. 5.3. Navier–Stokes flow over a step Finally, we give an example of Navier–Stokes flow over a step using increasing Reynolds numbers (see Figs. 5–8). The computational domain is given by X = (0, 4) · (0, 1)n(1.2, 1.6) · (0, 0.4), and the boundary conditions are: u = (0, 0) at the upper and lower parts of oX; u = (4x2(1 x2), 0) at the inflow; natural boundary condition at outflow (not traction free: the viscous operator was written as lDu for this example). We give the velocities and pressures (shown L2-projected onto the continuous space fv 2 C 0 ðXÞ : vjK 2 P 1 ðKÞ; 8K 2 Th g for ease of presentation) computed without the edge fluxes, and compare in Fig. 8 with a computation with fluxes. Clearly, the fluxes introduce too much artificial viscosity into
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2897
Fig. 4. Approximate solution of the velocities with jump of the convective derivative only.
Fig. 5. Approximate solution of the velocities and pressures at Reynolds number 100.
the method (at least when combined with the gradient jumps). This could be improved by tuning the gradient jump parameter, but the conclusion is that the flux terms are indeed not necessary. We remark that a real flow at a Reynolds number 10 000 certainly will not be laminar; our method thus introduces numerical dissipation, and the example in Fig. 8 only serves to show the amount of numerical dissipation.
6. Conclusion We have studied a non-conforming stabilized finite element method for incompressible flow. The velocities were approximated using the Crouzeix–Raviart element and the pressures were chosen as piecewise constants. The numerical scheme proposed is of interior penalty type and remains stable in all flow regimes
2898
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Fig. 6. Approximate solution of the velocities and pressures at Reynolds number 1000.
Fig. 7. Approximate solution of the velocities and pressures at Reynolds number 10 000.
Fig. 8. Approximate solution of the velocities and pressures at Reynolds number 10 000, with fluxes.
without streamline diffusion type stabilization. Instead we stabilize the jump in the streamline derivative between adjacent elements. We prove that this stabilization controls the part of the streamline derivative which is not already in the approximating space, allowing an optimal order a priori error estimate in the
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2899
energy norm which is uniform in the Peclet (Reynolds) number. Moreover the stabilizing term has the right asymptotic behaviour in the low Peclet regime and optimal order L2 estimates for the velocities are proved in this case. We present numerical results for different Reynolds numbers showing the robustness of the method and indicating optimal convergence of the error in the L2-norm. We also test the stability of the scheme on the lid-driven cavity flow and observe that the jumps in the streamline gradient is the most important stabilizing term. We believe that this scheme offers an attractive alternative to the ones proposed in [14] and in [12]. We have stability for all Reynolds numbers and may still lump mass for efficient timestepping.
Acknowledgement This paper was written while the second author was a guest of the Bernoulli center at the Ecole Polytechnique Fe´de´rale de Lausanne, working within a program for the modeling of cardiovascular flow. The support of the Bernoulli center is gratefully acknowledged.
References [1] S. Brenner, Korns inequalities for piecewise H1 vector fields, Math. Comp. 73 (2004) 1067–1087. [2] E. Burman, A unified analysis for conforming and non-conforming stabilized finite element methods using interior penalty. Preprint available from:. [3] E. Burman, P. Hansbo, Edge stabilization for Galerkin approximations of convection–diffusion–reaction problems, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1437–1453. [4] E. Burman, P. Hansbo, Stabilized Crouzeix–Raviart element for the Darcy–Stokes problem, Numer. Methods Partial Differ. Equat. 21 (2005) 986–997. [5] C. Carstensen, Merging the Bramble–Pasciak–Steinbach and the Crouzeix–Thome´e criterion for H1-stability of the L2-projection onto finite element spaces, Math. Comp. 71 (2002) 157–163. [6] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [7] M. Crouzeix, P.-A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations I, Rev. Franc¸aise Automat. Inform. Recherche Ope´r. Se´r. Rouge 7 (R-3) (1973) 33–75. [8] V. Girault, P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations, Springer-Verlag, Berlin, 1986. [9] V. John, G. Matthies, F. Schieweck, L. Tobiska, A streamline-diffusion method for nonconforming finite element approximations applied to convection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 166 (1998) 85–97. [10] V. John, J.M. Maubach, L. Tobiska, Nonconforming streamline-diffusion-finite-element-methods for convection–diffusion problems, Numer. Math. 78 (1997) 165–188. [11] L.I.G. Kovasznay, Laminar flow behind two-dimensional grid, Proc. Cambridge Philos. Soc. 44 (1948) 58–62. [12] G. Lube, L. Tobiska, A nonconforming finite element method of streamline diffusion type for the incompressible Navier–Stokes equations, J. Comput. Math. 8 (1990) 147–158. [13] G. Matthies, L. Tobiska, The streamline-diffusion method for conforming and nonconforming finite elements of lowest order applied to convection–diffusion problems, Computing 66 (2001) 343–364. [14] R. Temam, Navier–Stokes Equations, North-Holland, Amsterdam, 1984. [15] V. Thome´e, Galerkin Finite Element Methods for Parabolic Problems, Springer-Verlag, Berlin, 1997.
A stabilized non-conforming finite element method for incompressible flow Erik Burman a, Peter Hansbo a
b,*
Institut d’Analyse et de Calcul Scientifique, Ecole Polytechnique Fe´de´rale de Lausanne, CH-1015 Lausanne, Switzerland b Department of Applied Mechanics, Chalmers University of Technology, SE–41296 Go¨teborg, Sweden Received 31 January 2004; received in revised form 30 September 2004; accepted 23 November 2004
Abstract In this paper we extend the recently introduced edge stabilization method to the case of non-conforming finite element approximations of the linearized Navier–Stokes equation. To get stability also in the convective dominated regime we add a term giving L2-control of the jump in the gradient over element boundaries. An a priori error estimate that is uniform in the Reynolds number is proved and some numerical examples are presented. 2005 Elsevier B.V. All rights reserved. Keywords: Stabilized methods; Finite element; Crouzeix–Raviart; Non-conforming; Interior penalty; Incompressible flow
1. Introduction The solution of the Navier–Stokes equations for incompressible flow using finite element methods remains a challenging problem, in particular if the objective is to construct a method which remains robust and accurate for a wide range of Reynolds numbers. The discretization must assure not only satisfaction of the Babus˘ka–Brezzi condition but also stabilization of the convective terms and sufficient control of the incompressibility condition. Approximations using non-conforming Crouzeix–Raviart (CR) elements are attractive for the velocity approximation in combination with elementwise constant pressures, since they satisfy the Babus˘ka–Brezzi condition and have local conservation properties. This discretization was
*
Corresponding author. Fax: +31 772 3827. E-mail address: [email protected] (P. Hansbo).
0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.11.033
2882
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
proposed and analyzed in [14] and a stabilized version using the streamline diffusion stabilization was analyzed in [12]. In neither of these cases the high Reynolds number limit was treated. Moreover, in a recent paper, [4], the authors showed that a stabilized non-conforming finite element method using the Crouzeix– Raviart element remains uniformly stable in the vanishing viscosity limit for the generalized Stokes equation. However, the discretization of convection dominated problems using CR-elements are not stable without both a weak coupling on the element inflow boundary as in the discontinuous Galerkin method and streamline diffusion stabilization of the convective terms in the interior of the element (see [10,9,13]). In this respect, the element needs stabilization from both the streamline diffusion method and the discontinuous Galerkin method. Considering the fact that the method uses more degrees of freedom than the continuous Galerkin approximation this seems suboptimal. Moreover, the streamline diffusion stabilization has the drawback that it does not permit lumped mass for time stepping. In this paper we therefore propose to apply the recently introduced edge stabilization operator (see [2,3]) to the lowest order Crouzeix–Raviart element for the stabilization of the convective terms. We prove that this operator stabilizes exactly that part of the convective term which is not already included in the approximation space. In this sense this is the smallest perturbation needed to make the Crouzeix–Raviart element stable for convection– diffusion problems. This stabilization method has the advantage, as compared to other stabilized methods, that we may lump mass for efficient timestepping, we do not add any additional degrees of freedom, and we do not need any special structure of the mesh. For Oseens equation we prove an optimal a priori error estimate in the energy norm independent of the Reynolds number. Another attractive feature of the proposed stabilization is that, unlike SUPG, here the stabilization parameter is independent of the flow regime; we illustrate this by proving an L2 a priori error estimate for the velocities in the case of low local Reynolds number. Finally, we study the performance of the numerical scheme on some linear and nonlinear model cases.
2. A finite element method for the Oseen’s equation We consider, in X Rd with boundary oX, the problem of solving ru þ b ru þ rp 2lr eðuÞ ¼ f ;
in X;
r u ¼ 0; in X; u ¼ 0 on oX;
ð1Þ
d
where u, b 2 ½H 10 ðXÞ \ H 0 ðdiv; XÞ , b 2 W1,1(X), p 2 L20 ðXÞ, f is a given source term, r and l are bounded positive functions. By H0(div; X) we denote the functions in [L2(X)]d such that $ Æ u = 0, and by L20 ðXÞ the functions in L2(X) with zero mean value. In the following, we let (Æ, Æ) denote the L2-scalar product with the corresponding norm k Æ k = (Æ, Æ)1/2; the Hs(X) norm will be denoted by k Æ ks,X. d The weak form of (1) is to find ðu; pÞ 2 ½H 10 ðXÞ L20 ðXÞ such that (
aðu; vÞ þ bðp; vÞ ¼ ðf ; vÞ; bðq; uÞ ¼ 0;
8ðv; qÞ 2
½H 10 ðXÞd
ð2Þ
L20 ðXÞ;
where aðu; vÞ :¼ ðru; vÞ þ ðb ru; vÞ þ 2ðleðuÞ; eðvÞÞ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2883
and bðp; vÞ :¼ ðp; r vÞ. The well posedness of the above problem follows by the Lax-Milgram Lemma applied in the space ½H 10 ðXÞd \ H 0 ðdiv; XÞ. The finite element method consists of seeking a piecewise polynomial approximation uh 2 Vh, ph 2 Qh. We let Vh denote the space of the lowest order non-conforming Crouzeix–Raviart elements. Let Th S denote a shape regular triangulation of the domain X, EðKÞ the set of all faces of an element K 2 Th , E :¼ K2Th EðKÞ the set of all faces in Th , EoX :¼ fe 2 E : e oXg, and E0 :¼ E n EoX the set of the boundary and inner faces respectively. For a given piecewise continuous function u, the jump [u] and the average {u} on a face e 2 E are defined by 8 < limþ ðuðx snÞ uðx þ snÞÞ if e 6 oX; s!0 ½uðxÞ :¼ : lim ¼ uðx snÞ if e oX; s!0þ 81 < 2 limþ ðuðx þ snÞ þ uðx snÞÞ if e 6 oX; s!0 fugðxÞ :¼ 1 : 2 lim 2uðx snÞ if e oX; þ s!0
where n is a normal unit vector on e and x 2 e. If e oX we choose the orientation of n to be outward with respect to X otherwise n has an arbitrary but fixed orientation. For the non-conforming finite element functions, continuity across edges e will only be enforced with respect to Z je ðvh Þ :¼ ½vh ds. e
Using this definition our finite element space may be defined as V h :¼ fvh 2 ½L2 ðXÞd : vh jK 2 ½P 1 ðKÞd
8K 2 Th ;
je ðvh Þ ¼ 0; 8e 2 Eg.
Moreover, we introduce the space of piecewise constants with mean value zero, Qh :¼ fqh 2 L20 ðXÞ : qh jK 2 P 0 ðKÞg and the subspace Wh of Vh such that W h :¼ fwh 2 V h : ðr wh ; qh Þh ¼ 0 8qh 2 Qh g P where ðr wh ; qh Þh ¼ K ðr wh ; qh ÞK . Since the above spaces are H1-non-conforming we introduce the broken norm equivalent of the L2-norm X 2 2 kukh ¼ kukK K
and the broken H1-seminorm X 2 juj2h ¼ juj1;K . K
The local mesh size is defined by hK :¼ max hoK ; K
and we will assume that hK/hoK < C where C is a fixed constant. We will use C and c as generic constants taking different values at different occasions. To indicate their provenance or main dependence, a subscript may be added, e.g., cl. We introduce the interpolation operator rhu : [H1(X)]d ! Vh defined by
2884
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
rh uðxe Þ ¼
1 jej
Z u ds; e
where xe is the midpoint of the edge e. The L2-projections onto the spaces are also required for the analysis. Let p0,h : L2(K) ! P0(K), pd0;h : ½L2 ðKÞd ! ½P 0 ðKÞd denote the L2-projection onto the constant functions on K and p1,h : [L2(K)]d ! Vh the standard L2-projection onto the finite element space. For the above defined interpolation operator and projections we need some approximation and stability properties. These, and some inverse inequalities are collected in the following Lemmas. Lemma 1. For the interpolation operator rh there holds, if u 2 [H2(X)]d then krh u ukh þ hjrh u ujh 6 C r h2 .
ð3Þ
Moreover, if $ Æ u = 0 then rhu 2 Wh. Proof. The proof of the interpolation estimate is given by Crouzeix and Raviart [7]. The second claim is immediate noting that Z Z r rh u dx ¼ r u dx ¼ 0 K
K
by the definition of the interpolant.
h
Lemma 2. For the L2-projection the following H1 stability holds, jp1;h ujh 6 C s kuk1;X . For the error analysis, we shall use the following trace inequality. Lemma 3. For v 2 H1(K) there holds 2 2 kvk2oK 6 C t ðh1 K kvkK þ hK kvk1;K Þ
8v 2 H 1 ðKÞ;
ð4Þ
where Ct is a constant independent of hK. We also need the following local inverse inequality. Lemma 4. Let uh 2 Vh where Vh is defined on a shape regular mesh then kruh kK 6 h1 K C i kuh kK ; with Ci independent of K. Proof. For proofs of Lemmas 2–4, see, respectively, Carstensen [5], Thome´e [15], and Ciarlet [6]. Our finite element method reads: find uh 2 Vh such that ah ðuh ; vh Þ þ bh ðph ; vh Þ þ ju ðuh ; vh Þ ¼ ðf ; vh Þ bh ðqh ; uh Þ ¼ 0 8ðvh ; qh Þ 2 V h Qh ;
h
ð5Þ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2885
where ah ðuh ; vh Þ ¼ ðruh þ b ruh ; vh Þh
1X hb n½uh ; fvh gioK þ ð2leðuh Þ; eðvh ÞÞh ; 2 K
bh ðph ; vh Þ ¼ ðph ; r vh Þh ;
ð6Þ ð7Þ
and the jump terms take the form XZ cb 2 ju ðuh ; vÞ ¼ hoK bh ½ruh bh ½rvds oKnoX jbh j K XZ ca ðl þ jb njhoK ÞhoK ½t ruh ½t rvds þ K
þ
X K
oK
Z cb
hoK ½ðt ruh Þ n½ðt rvÞ nds.
ð8Þ
oK
Here bh is the interpolant of b on Wh and t is a unit vector perpendicular to n. The gradient jump term serves three purposes. It stabilizes the convective terms (the first sum in (8)), it assures that Korns inequality is satisfied (the cal part of the second sum in (8)), it gives additional control of the divergence inconsistency error (the third sum in (8)). The last two properties can be obtained by introducing a lower order penalizing term (see [4]), but the use of the jump of the gradient has the advantage of allowing for one point quadrature in the implementation and from the point of view of analysis it is practical. In the case of three space dimensions the tangent vector should be replaced by the tangent tensor $uh · n. In the following we will for simplicity only consider the two dimensional case for the tangent vectors. Remark 5. The analysis below holds with only minor modifications if b is replaced by bh also in the convective term. This is convenient when timestepping the Navier–Stokes equation: b may be taken as the solution uh of the previous timestep. In the analysis we will not distinguish between the different stabilization parameters, they will all be denoted c. We introduce the following shorthand notation A½ðuh ; ph Þ; ðwh ; qh Þ ¼ ah ðuh ; vh Þ þ bh ðph ; vh Þ bh ðqh ; uh Þ. To simplify the analysis we will assume that the exact solution (u, p) belongs to [H2(X)]d · H1(X); it then follows that the formulation (5) enjoys the following consistency property. Lemma 6. For u 2 H2(X) and p 2 H1(X) there holds A½ðu uh ; p ph Þ; ðvh ; qh Þ þ ju ðu uh ; vh Þ ¼ Rðu; p; vh Þ
ð9Þ
for all (vh, qh) 2 Vh · Qh. Where the consistency error due to the non-conforming approximation is given by 1X 1X Rðu; p; vh Þ ¼ h2leðuÞ n; ½vh ioK þ hp; ½vh nioK . 2 K 2 K Proof. This is an immediate consequence of the regularity hypothesis: if u 2 H2(X) then the trace of $u is well defined and hence j(u, vh) = 0. The consistency error is obtained by integration by parts, X X 1 ð2leðuÞ; eðvh ÞÞ ¼ ð2lr eðuÞ; vh Þ þ h2leðuÞ n; ½vh ioK 2 K K
2886
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
and
X
ðp; r vh ÞK ¼ ðrp; vh Þ
K
1X hp; ½vh nioK . 2 K
2.1. Preliminary lemmas In this section we will prove some preliminary results that will facilitate the analysis. The main result of this section is Lemma 7 where we show that that the jump term (8) controls the difference between the convective derivative and its quasi-interpolant on the finite element space. This lemma is the key ingredient to derive error bounds that are independent of the Peclet number. Then we prove the coercivity of the bilinear form. The triple norm that we will use is given by jjjðwh ; qh Þjjj2 ¼ kr1=2 wh k2h þ jl1=2 wh j2h þ kr wh k2h þ ju ðwh ; wh Þ þ cp kqh k2 ; where cp is a constant depending on the problem parameters r, l, b to be specified later. Note that although for the Crouzeix–Raviart element $ Æ uh = 0 on each triangle the formulation must include a stabilization of the jump of the normal velocity due to the H(div, X) consistency error. Therefore the triple norm must be d chosen as a discrete norm on ðH ðdiv; XÞ \ l1=2 ½H 10 ðXÞ Þ L20 ðXÞ. The triple norm is dominated, at low Rey1 nolds numbers, by the H (X) contribution, and at high Reynolds numbers by the part of the jump term controlling the inconsistency in the divergence. This latter term prohibits the decoupling of the velocities and the pressure and the order of the estimate can be no better that the approximation properties of the pressure space. We introduce the space of functions that are piecewise linear on each element d
Y h ¼ ½fy 2 L2 ðXÞ : yjK 2 P 1 ðKÞg . h : Y h ! V h . Let xi be the midpoint of the We now introduce a quasi-interpolant based on local averages p face shared by element K and element K 0 then fuðxi Þg for xi an interior node; h uðxi Þ ¼ p uðxi Þ for xi a boundary node. In the following Lemma we prove that the projection error is bounded by the jumps in the gradient. Lemma 7. Let bh 2 Vh and wh 2 Vh then 2
h ðbh rwh ÞÞk 6 jb ðwh ; wh Þ; kh1=2 ðbh rwh p where jb(wh, wh) is given by X Z jb ðwh ; wh Þ ¼ cb K
hK hoK ? ðbh ½rwh Þ2 ds
oKnoX
with hoK ? denoting the triangle size perpendicular to the side on oK and cb is a parameter depending only on the number of space dimensions. h ðbh rwh Þ makes sense. Now consider any Proof. First note that (Vh Æ $)Vh Yh so that the projection p d triangle K and note that using the Crouzeix–Raviart basis functions fui gi¼1 , with ui associated with node xi we may write bh rwh jK ¼
dþ1 X i¼1
bh rwh ðxi Þui ðxÞ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2887
and h ðbh rwh ÞjK ¼ p
dþ1 X
fbh rwh ðxi Þgui ðxÞ.
i¼1
Taking the difference of the two functions in a nodal point yields bh rwh ðxi ÞjK fbh rwh ðxi Þg ¼ bh ðxi Þ ½rwh ðxi Þ.
Note that if the node is on the boundary, the right hand side is zero. It follows that for any K 2 Th such that oK \ oX = ; kh
1=2
h bh ðbh rwh p
rwh Þk2K
¼
Z
dþ1 X ðbh rwh ðxi Þ fbh rwh gÞui
hK K
dx
i¼1
Z
dþ1 X
hK
6
!2
K
!2 bh ½rwh ui
dx.
ð10Þ
i¼1
We now evaluate the integral using nodal point quadrature and note that since the nodes of the Crouzeix– Raviart element are on the midpoints of the element sides this is exact for second degree polynomials in two space dimensions, hence 1=2
2
h bh rwh ÞkK ¼ khK ðbh rwh p
3 X k¼1
hK
3 X mK moK k 2 2 bh ðxk Þ ½ruh ¼ bh ðxk Þ ½ruh ðxk Þ ; hK h? oK 3 6 k¼1
ð11Þ
R where moK k ¼ oK k dx with oKk the face associated with quadrature point k. In three space dimensions the midpoints of the faces has to be supplemented with the six midpoints of the edges of the tetrahedron to yield an exact quadrature formula (with weights 1/15 for the midpoints of the faces and 3/20 for the midpoints of the edges). One may then easily show that 1=2
2
h bh rwh ÞkK 6 khK ðbh rwh p
3 5X moK k 2 hK h? ðbh ðxk Þ ½ruh ðxk ÞÞ . oK 6 k¼1 6
ð12Þ
It follows, using the Simpson quadrature formula in two dimensions and a quadrature taking the midpoint d , that of the face and the corner-points in three dimensions and noting that the weight for the midpoint is dþ1 Z 2 h ðbh rwh ÞÞk2K 6 cd kh1=2 ðbh rwh p hK hoK ? ðbh ½ruh Þ ds; ð13Þ oK
where cd = 1/4 in two dimensions and cd = 10/9 in three dimensions (the subscript d is here referring to ‘‘dimension’’). We conclude by taking the sum over all triangles K 2 T noting that all boundary contributions vanishes thanks to the definitions of the quasi-interpolant. Remark 8. A consequence of the above proof is that the jump term edge integral may be evaluated using midpoint quadrature on the faces. In fact the first part of the jump operator given in (8) can be substituted by the discrete operator given by (11) in two space dimensions and by (12) in three to get optimal values of the stabilization constants. The integral formulation however still remains practical from a theoretical viewpoint since it is consistent for H2-regular solutions. As was pointed out in the previous section we only stabilize using the jumps in the gradient. However we need to establish a result showing the equivalence between the jumps in the solution and the jumps in the tangential gradient.
2888
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Lemma 9. For the interior penalty term (8) there holds X 2 kl1=2 h1=2 ½wh koK 6 cju ðwh ; wh Þ; K
X
2
kh1=2 ½wh nkoK 6 cju ðwh ; wh Þ;
K
and ju ðwh ; wh Þ 6 cc jwh jh for all wh 2 Vh. Proof. By the midpoint continuity of the Crouzeix–Raviart element we note that we may write, with n a coordinate along e with midpoint ni, ½wh ðnÞje ¼ ½t rwh je ðn ni Þ and ½wh ðnÞ nje ¼ ½ðt rwh Þ nje ðn ni Þ. Hence we have Z Z Z 1 2 2 2 ½wh ðnÞ dn ¼ ð½t rwh ðn ni ÞÞ dn ¼ h2 ½t rwh dn; 12 e e e e which proves the first claim. The proof of the second claim is equivalent. The last claim finally is an immediate consequence of the trace inequality (4). h Lemma 10. For the consistency error the following upper bound holds jRðu; p; vh Þj 6 ðcl hkuk2;X þ Chkpk1;X Þju ðvh ; vh Þ. Proof. Using the zero mean value property of the Crouzeix–Raviart space followed by the Cauchy– Schwarz inequality we obtain
Rðu; p; vh Þ 6 c
X
!1=2 kh
1=2
ð2lÞ
1=2
ðeðuÞ n
pd0;h eðuÞ
2 nÞkoK
K
þc
X
!1=2 1=2 1=2
kl
h
2 ½vh koK
K
X
!1=2 hK kp
p0;h pk2oK
K
X
!1=2 h1 K k½vh
nk2oK
ð14Þ
.
K
The claim now follows using the trace inequality (4), standard interpolation and Lemma 9.
h
Let us now investigate the coercivity properties of the discretization of the convective terms. Lemma 11. For the convective terms there holds 1X ðb ruh ; uh Þh ¼ hb n½uh ; fuh gi. 2 K
ð15Þ
Proof. The result follows by integrating by parts elementwise in the left hand side and using the equality [xy] = [x]{y} + {x}[y]. The integration by parts yields
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Z X X 1X ðb ruh ; uh ÞK ¼ ½b nuh uh ds ðuh ; b ruh ÞK . 2 K oK K K
2889
ð16Þ
We rewrite the edge term in the following fashion ½b nuh uh ¼ 2b n½uh fuh g
ð17Þ h
and the proof is completed using (17) in (16).
To prove the coercivity of our operator we also need the following discrete Korns inequality. Lemma 12. For wh 2 Vh there holds C 2K jl1=2 wh j2h 6 kð2lÞ1=2 eðwh Þk2h þ ju ðwh ; wh Þ. Proof. See Brenner [1].
h
The coercivity of our formulation is an immediate consequence of Lemmas 11 and 12. Lemma 13. For all (wh, qh) 2 Vh · Qh there holds C 2K jl1=2 wh j2h þ kr1=2 wh k2 6 A½ðwh ; qh Þ; ðwh ; qh Þ. Proof. First of all notice that the terms bh(uh, ph) cancel. We may write 2
A½ðwh ; qh Þ; ðwh ; qh Þ ¼ kr1=2 wh kh þ kð2lÞ
1=2
2
eðwh Þkh þ ðb rwh ; wh Þh
1X hb n½wh ; fwh gi. 2 K
Using Lemma 11 for the convective term we get ðb rwh ; wh Þ
1X hb n½wh ; fwh gi ¼ 0. 2 K
The proof is then completed by applying Lemma 12.
h
3. Stability In this section we will prove an inf–sup condition for our discretization vital for the convergence analysis. Theorem 14 (Stability). If (uh, ph) 2 Vh · Qh then there holds cis jjjðuh ; ph Þjjj 6
sup ðwh ;qh Þ2V h Qh
A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ ; jjjðwh ; qh Þjjj
ð18Þ
where the constant cis depends only on the parameters l, r, b, c and remains bounded from below when l ! 0. Proof. We prove the above inf–sup condition in two steps. First we will prove that there exists (wh, qh) 2 Vh · Qh such that 2
jjjðuh ; ph Þjjj 6 A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ
ð19Þ
2890
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
and then we conclude by proving that jjjðwh ; qh Þjjj 6 Cjjjðuh ; ph Þjjj.
For the first step we note that by Lemma 13 we have choosing wh = uh, qh = ph 2
2
C 2K jl1=2 uh jh þ kr1=2 uh k þ ju ðuh ; uh Þ 6 A½ðuh ; ph Þ; ðuh ; ph Þ þ ju ðuh ; uh Þ.
ð20Þ
2
To control the L -norm of the pressure we note, following [8], that by the surjectivity of the divergence d operator there exists a function vp 2 ½H 10 ðXÞ such that $ Æ vp = ph and jvpj1,X 6 ckphk. Therefore we now choose wh = rhvp and qh = $ Æ uh. By the stability of the quasi-interpolation operator rh we have krh vp k þ jrh vp jh 6 ckph k.
ð21Þ
Moreover, using the definition of the quasi-interpolant rhvp we have X X 2 kph k ¼ ðph ; r vp Þ ¼ h½ph ; vp nie ¼ h½ph ; rh vp nie ¼ ðph ; r rh vp Þh . e
ð22Þ
e
As a consequence of (22) we may write 2
2
A½ðuh ; ph Þ; ðrh vp ; r uh Þ þ ju ðuh ; rh vp Þ ¼ kph k þ kr uh k þ ðruh ; rh vp Þ þ ðlruh ; rrh vp Þh 1X þ ðb ruh ; rh vp Þ hb n½uh ; frh vp gioK þ ju ðuh ; rh vp Þ. 2 K ð23Þ Using Cauchy–Schwarz inequality, Youngs inequality, and the stability (21), we readily deduce 2
2
ðruh ; rh vp Þ P cr kr1=2 uh k 18kph k ; ð2leðuh Þ; eðrh vp ÞÞh P cl jl
1=2
uh j2h
ð24Þ
1 kph k2 8
ð25Þ
and by applying the third inequality of Lemma 9 ju ðuh ; rh vp Þ P ju ðuh ; uh Þ
1=2
1=2
ju ðrh vp ; rh vp Þ
2
P cc ju ðuh ; uh Þ 18kph k .
ð26Þ
It now remains to bound the convective term. An integration by parts yields ðb ruh ; rh vp Þ
1X 1X hb n½uh ; frh vp gioK ¼ ðuh ; b rrh vp Þ hb nfuh g; ½rh vp ioK . 2 K 2 K
ð27Þ
For the first term we clearly have ðuh ; b rrh vp Þ P 2ckbkL1 ðXÞ kuh k kph k.
ð28Þ
For the second we use the H1-regularity of vp, the trace inequality (4) and the local inverse inequality to obtain X 1X hb nfuh g; ½rh vp vp ioK P 2 kbkL1 ðKÞ kuh koK krh vp vp koK 2 K K X 2 2 1=2 1=2 P 2 kbkL1 ðKÞ ðh1 hK kvp k1;K K kuh kK þ hK kuh k1;K Þ K
P 2C i ckbkL1 ðXÞ kuh kkph k.
ð29Þ
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2891
Combining (28) and (29) we obtain 1X 1 ðb ruh ; rh vp Þ hb n½uh ; rh vp ioK P cb kuh k2X kph k2 . 2 K 8
ð30Þ
Using now (24)–(26) and (30) to find a lower bound for the expression (23) we get 2
2
2
2
A½ðuh ; ph Þ; ðrh vp ; r uh Þ þ ju ðuh ; rh vp Þ ¼ 12kph k þ kr uh k C lrbc ðkr1=2 uh k þ jl1=2 uh jh þ ju ðuh ; uh ÞÞ; ð31Þ where C lrbc
! cb cl ¼ max cr ; ; 2 ; cc . r CK
Combining (20) and (31) we conclude that (19) holds for the choice wh = uh + (2Clrbc)1rhvp, qh = ph + $ Æ uh. More precisely we have, with cp = (2Clrbc)1, 1 jjjðuh ; ph Þjjj2 2
6 A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ.
ð32Þ
It remains to prove that jjjðwh ; qh Þjjj 6 Cjjjðuh ; ph Þjjj; this is obtained simply by noting that 2
2
2
2
2
jl1=2 wh jh 6 jl1=2 uh jh þ jl1=2 cp rh vp jh 6 jl1=2 uh jh þ lc2p ckph k .
ð33Þ
Proceeding in the same fashion for the other terms in the triple norm yields jjjðwh ; qh Þjjj2 ¼ jl1=2 wh j2h þ kr1=2 wh k2 þ ju ðwh ; wh Þ þ kr wh k2 þ kqh k2 6 jl1=2 uh j2h þ kr1=2 uh k2 þ jðuh ; uh Þ þ 2kr uh k2 þ Ckph k2 6 Cjjjðuh ; ph Þjjj2
ð34Þ
and we conclude that 2
cis jjjðuh ; ph Þjjj jjjðwh ; qh Þjjj 6 jjjðuh ; ph Þjjj 6 A½ðuh ; ph Þ; ðwh ; qh Þ þ ju ðuh ; wh Þ. Clearly the constant cis is independent of h and furthermore it does not vanish for vanishing l.
4. Error estimates We will now proceed to derive a priori error estimates in the triple norm. The estimate takes the form jjjðu uh ; p ph Þjjj 6 Ch. The energy norm estimate is independent of the Peclet number, indicating that the proposed method should be stable for a wide range of Reynolds numbers when applied to the full Navier–Stokes equations. Lemma 15 (Approximation). Consider the projection (p1,hu, p0,hp) of the exact solution d
ðu; pÞ 2 ½H 2 ðXÞ H 1 ðXÞ onto the finite element space Vh · Qh. For the projection error there holds jjjðp1;h u u; p0;h p pÞjjj 6 Ch; where C 6 c(1 + l1/2 + jbj1/2h1/2 + r1/2h + c1/2)kuk2,X + cpkpk1,X.
2892
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Proof. By the optimal approximation property of the L2-projection we have kr1=2 ðp1;h u uÞk 6 cr h2 and kp0;h p pk 6 Ch. We now consider rhu p1,hu. Noting that rhu p1,hu = p1,h(rhu u) we obtain using the H1-stability of the L2-projection for the Crouzeix–Raviart element on locally quasi-uniform meshes (see [5]), krðrh u p1;h uÞkh 6 C s krðrh u uÞkh 6 C s C r h
ð35Þ
and we conclude that kl1=2 rðu p1;h uÞkh 6 cl h and kr ðu p1;h uÞkh 6 Ch. Finally we estimate the penalizing term ju(p1,hu u, p1,hu u). We have that XZ XZ 2 2 ½rðp1;h u uÞ ds 6 2 ðrðp1;h u uÞÞ ds oK oK K K X 2 2 ðh1 62 K krðp1;h u uÞkK þ hK kuk2;K Þ 6 Chkuk2;X ; K
where we have used the trace inequality (4) in the second inequality. We thus conclude that ju ðp1;h u u; p1;h u uÞ
1=2
6 cc1=2 ð1 þ l1=2 þ jbj
1=2 1=2
h
Þhkuk2;X .
Lemma 16 (Continuity). Let g = p1,hu u and f = p0,hp p be the projection error of the velocity and the pressure respectively, then there holds A½ðg; fÞ; ðwh ; qh Þ þ ju ðg; wh Þ Rðu; p; wh Þ 6 Cðjjjðg; fÞjjj þ ðcl þ cb h1=2 Þhkuk2;X þ Chkpk1;X Þjjjðwh ; qh Þjjj. Proof. Clearly we have using Cauchy–Schwarz inequality ðrg; wh Þ 6 Ckr1=2 gk jjjðwh ; 0Þjjj 6 cr jjjðg; fÞjjj jjjðwh ; qh Þjjj; ð2leðgÞ; eðwh ÞÞ 6 Cjl1=2 gjh jjjðwh ; 0Þjjj 6 cl jjjðg; fÞjjj jjjðwh ; qh Þjjj; ju ðg; wh Þ 6 jjjðg; fÞjjj jjjðwh ; qh Þjjj. We consider now the terms expressing the pressure velocity coupling. By the orthogonality of the L2-projection p0,h we have bðf; wh Þ ¼ ðp0;h p p; r wh Þh ¼ 0.
Using once again the Cauchy–Schwartz inequality we readily obtain bðqh ; gÞ ¼ ðqh ; r gÞh 6 jjjðwh ; qh Þjjj jjjðg; fÞjjj. It remains to treat the convective term and the non-consistency term. Let us first consider the convective term, an integration by parts yield together with the addition and substraction of bh yields 1X hb n½g; fwh gi ¼ ðg; ðb bh Þ rwh Þ ðg; bh rwh Þ ðb rg; wh Þ 2 K 1X hb n; fgg½wh i ¼ I þ II þ III. ð36Þ þ 2 K
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2893
The first term is controlled using a local inverse inequality and a local approximation result for b bh. !1=2 X 2 2 2 I 6 kgk kbkW 1;1 ðKÞ hK krwh k 6 cb kgkC i kwh k 6 cb h2 kuk2;X jjjðwh ; qh Þjjj. K
Using Lemma 7 we have for the term II h ðbh rwh ÞÞÞ 6 kh1=2 gk kh1=2 ðbh rwh p h ðb rwh ÞÞk II ¼ ðg; ðbh rwh p 6 c1=2 kh1=2 gkju ðwh ; wh Þ
1=2
6 cc h3=2 kuk2;X jjjðwh ; qh Þjjj.
ð37Þ
For the term III we obtain using the trace inequality (4) and the approximation properties of the L2projection. !1=2 X 1=2 1=2 III 6 hjb nj ½wh ; ½wh ioK kb nkL1 ðKÞ kfggkoK K
6 ju ðwh ; wh Þ1=2
X
1=2
kb nkL1 ðKÞ kfggkoK 6 cb h3=2 kuk2;X jjjðwh ; qh Þjjj.
K
Only the non-consistent residual remains to be bounded and we conclude the proof by applying Lemma 10. Theorem 17. Let (u, p) 2 [H2(X)]d · H1(X) be the solution of (1) and let (uh, qh) 2 Vh · Qh be the finite element solution of (5). Then there holds jjjðu uh ; p ph Þjjj 6 ðcr h þ cl þ cb h1=2 þ cc Þhkuk2;X þ Chkph k1;X .
ð38Þ
Proof. We will consider the discrete error ehu ¼ p1;h u uh and ehp ¼ p0;h p ph since using Lemma 15 we have jjjðu uh ; p ph Þjjj 6 jjjðu p1;h u; p p0;h pÞjjj þ jjjðp1;h u uh ; p0;h p ph Þjj 6 Ch þ jjjðp1;h u uh ; p0;h p ph Þjjj; where C is the constant given in Lemma 15. By Lemma 14 we have cis jjjðehu ; ehp Þjjj 6
A½ðehu ; ehp Þ; ðwh ; qh Þ þ ju ðehu ; wh Þ . jjjðwh ; qh Þjjj
ð39Þ
Using now Galerkin orthogonality we may write cis jjjðehu ; ehp Þjjj 6
A½ðg; fÞ; ðwh ; qh Þ þ ju ðg; wh Þ Rðu; p; wh Þ ; jjjðwh ; qh Þjjj
where g = p1,hu u and f = p0,hp p. We conclude the proof by applying Lemmas 16 and 15.
ð40Þ h
Remark 18. Since the estimate is only first order, due to the low order approximation of the pressure and the inclusion of the divergence in the triple norm, the virtues of the streamline stabilization are not obvious. We could in fact proceed with an inverse inequality in term II of (36) and still have the same formal convergence order of the triple norm. It is however known that this would destroy stability of the velocities for problems with important gradients. This is illustrated in the numerical section. As we mentioned in Section 1 the interior penalty method is independent of the Reynolds number. Of course, for low local Reynolds number the numerical scheme is stable without stabilization (except for the Korns inequality), so that the stabilization may be eliminated. We will however show that this is unnecessary
2894
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
for our discretization by proving that even when keeping the stabilizing terms our discretization has optimal L2-convergence of the velocities in the local low Reynolds number regime in spite of the fact that cb is independent of l. Consider the dual continuous problem of seeking / 2 ½H 10 ðXÞ and r 2 L20 ðXÞ such that r/ b r/ 2lr eð/Þ þ rr ¼ e in X; r/¼0
ð41Þ
in X;
where e :¼ u uh, and assume that we have the regularity estimate k/kH 2 ðXÞ þ krkH 1 ðXÞ 6 kek.
ð42Þ
We recall that since uh 2 Wh and rh/h 2 Wh we have $ Æ e = 0 and $ Æ (/ rh/) = 0 elementwise. Multiply the first line of (41) by e and integrate by parts to obtain X X h2ln eð/Þ; ½eioK þ hr; ½n eioK ; kek2 ¼ ah ðe; /Þ þ ju ðe; /Þ K
K
where we have used the partial integration X X X 1X ðe; b r/Þ ¼ ðb re; /ÞK hb ne; /ioK ¼ ðb re; /ÞK hb n½e; f/gioK 2 K K K K and the divergence free property of the error and of /. Using Galerkin orthogonality, the divergence free property of the interpolant and the zero mean value property of the Crouzeix–Raviart element, we obtain X 2 kek ¼ ah ðe; / rh /Þ þ ju ðe; / rh /Þ h2lðn eð/Þ pd0;h n eð/ÞÞ; ½eioK
X
K
h2lðn eð/Þ
pd0;h n
eð/ÞÞ; ½/ rh /ioK
K
þ
X
hr p0;h r; ½n eioK þ
K
X hp p0;h p; ½n ð/ rh /ÞioK K
6 jjjðe; 0Þjjj jjjð/ rh /; 0Þjjj þ cb ju ðe; eÞ
1=2
X K
þ
X
hK k2l1=2 ðn eð/Þ pd0;h n
K
þ
X
þ
þ
!1=2 jjjðe; 0Þjjj
!1=2 hK kr p0;h rkL2 ðoKÞ
jjjðe; 0Þjjj !1=2
hK k2l1=2 ðn eð/Þ pd0;h n
K
X
k/
2 eð/ÞÞkL2 ðoKÞ
K
X
!1=2 2 rh /koK
2 eð/ÞÞkL2 ðoKÞ
jjjð/ rh /; 0Þjjj
!1=2 hK kp p0;h pkL2 ðoKÞ
jjjð/ rh /; 0Þjjj.
K
Using Lemma 15, the trace inequality (4) and error estimates for rh and for piecewise constant interpolation, we arrive at 2
kek 6 Chðjjjðe; 0Þjjj þ hkukH 2 ðXÞ þ hkpkH 1 ðXÞ Þðk/kH 2 ðXÞ þ krkH 1 ðXÞ Þ and thus we have
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2895
Theorem 19. Under the regularity assumption (42), the L2-error in the velocities can be estimated as kek 6 Ch2 ðkukH 2 ðXÞ þ kpkH 1 ðXÞ Þ.
ð43Þ
5. Numerical results 5.1. Convergence study in the case of small viscosity Let k = (l1 (l2 + 16p2)1/2)/2. Then the exact solution to (1) is given by u1 ðx1 ; x2 Þ ¼ 1 ekx1 cos 2px2 ; kx1 k e sin 2px2 ; u2 ðx1 ; x2 Þ ¼ 2p 1 p ¼ e2kx þ C; 2 cf. [11], with b = u, r = 0 and a right hand side matching the exact solution. In our examples, we also chose C to give zero mean pressure. We solved this problem approximatively on X = (1/2, 3/2) · (0, 2), using stability parameters ca = cb = 1/100 and cb = 1/4. Concerning the size of the stabilizing parameters, our numerical experience is that they should be of approximately this size; in particular the scheme becomes excessively diffusive with parameters around or above 1. The parameters in this method should thus be chosen smaller than, e.g., in the SUPG method. In Fig. 1 we show the convergence for l = 103, and in Fig. 2 for l = 105. Note that the absolute value of the pressure decreases linearly with the inverse of l in L2, which is why the absolute error in pressure is smaller in Fig. 2. 5.2. Stability in the Navier–Stokes case We show the influence of the different stabilizing terms for the lid-driven cavity flow in X = (0, 1) · (0, 1) and with l = 103. In Fig. 3 we show the numerical solution of Navier–Stokes equations (with b = u) after
-4 Velocity Pressure
log (L2 -error)
-4.5 -5 -5.5 -6 -6.5 -7 -4.5
-4
-3.5
-3
-2.5
-2
-1.5
log (h)
Fig. 1. Convergence for l = 103.
-1
-0.5
0
2896
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899 -3.5 -4 -4.5
log (L2-error)
-5
Velocity Pressure
-5.5 -6 -6.5 -7 -7.5 -8 -8.5 -5
-4
-3
-2
-1
0
1
log (h)
Fig. 2. Convergence for l = 105.
Fig. 3. Approximate solution of the velocities, (left) with fluxes only and (right) with fluxes plus jump of the convective derivative.
15 fixed point iterations, we present the solution using only upwind fluxes, as well as the fully stabilized solution. Note that we do not get a wiggle free solution without the jump in the convective derivative. On the other hand, in Fig. 4 we show the solution obtained using only the jump in convective derivative as stabilization. This solution is markedly less diffusive, yet still retains stability. 5.3. Navier–Stokes flow over a step Finally, we give an example of Navier–Stokes flow over a step using increasing Reynolds numbers (see Figs. 5–8). The computational domain is given by X = (0, 4) · (0, 1)n(1.2, 1.6) · (0, 0.4), and the boundary conditions are: u = (0, 0) at the upper and lower parts of oX; u = (4x2(1 x2), 0) at the inflow; natural boundary condition at outflow (not traction free: the viscous operator was written as lDu for this example). We give the velocities and pressures (shown L2-projected onto the continuous space fv 2 C 0 ðXÞ : vjK 2 P 1 ðKÞ; 8K 2 Th g for ease of presentation) computed without the edge fluxes, and compare in Fig. 8 with a computation with fluxes. Clearly, the fluxes introduce too much artificial viscosity into
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2897
Fig. 4. Approximate solution of the velocities with jump of the convective derivative only.
Fig. 5. Approximate solution of the velocities and pressures at Reynolds number 100.
the method (at least when combined with the gradient jumps). This could be improved by tuning the gradient jump parameter, but the conclusion is that the flux terms are indeed not necessary. We remark that a real flow at a Reynolds number 10 000 certainly will not be laminar; our method thus introduces numerical dissipation, and the example in Fig. 8 only serves to show the amount of numerical dissipation.
6. Conclusion We have studied a non-conforming stabilized finite element method for incompressible flow. The velocities were approximated using the Crouzeix–Raviart element and the pressures were chosen as piecewise constants. The numerical scheme proposed is of interior penalty type and remains stable in all flow regimes
2898
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
Fig. 6. Approximate solution of the velocities and pressures at Reynolds number 1000.
Fig. 7. Approximate solution of the velocities and pressures at Reynolds number 10 000.
Fig. 8. Approximate solution of the velocities and pressures at Reynolds number 10 000, with fluxes.
without streamline diffusion type stabilization. Instead we stabilize the jump in the streamline derivative between adjacent elements. We prove that this stabilization controls the part of the streamline derivative which is not already in the approximating space, allowing an optimal order a priori error estimate in the
E. Burman, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 195 (2006) 2881–2899
2899
energy norm which is uniform in the Peclet (Reynolds) number. Moreover the stabilizing term has the right asymptotic behaviour in the low Peclet regime and optimal order L2 estimates for the velocities are proved in this case. We present numerical results for different Reynolds numbers showing the robustness of the method and indicating optimal convergence of the error in the L2-norm. We also test the stability of the scheme on the lid-driven cavity flow and observe that the jumps in the streamline gradient is the most important stabilizing term. We believe that this scheme offers an attractive alternative to the ones proposed in [14] and in [12]. We have stability for all Reynolds numbers and may still lump mass for efficient timestepping.
Acknowledgement This paper was written while the second author was a guest of the Bernoulli center at the Ecole Polytechnique Fe´de´rale de Lausanne, working within a program for the modeling of cardiovascular flow. The support of the Bernoulli center is gratefully acknowledged.
References [1] S. Brenner, Korns inequalities for piecewise H1 vector fields, Math. Comp. 73 (2004) 1067–1087. [2] E. Burman, A unified analysis for conforming and non-conforming stabilized finite element methods using interior penalty. Preprint available from: