RANDOLPH E. BANKy AND BRUNO D. WELFERTz. Abstract. We derive and ..... ( ]A) and jump ( ]J) of a function v across an edge e are de ned edgewise by: 4 ...
A POSTERIORI ERROR ESTIMATES FOR THE STOKES PROBLEM
RANDOLPH E. BANKy AND BRUNO D. WELFERTz
Abstract. We derive and analyze an a posteriori error estimate for the mini-element discretization of the Stokes equations. The estimate is based on the solution of a local Stokes problem in each element of the nite element mesh, using spaces of quadratic bump functions for both velocity and pressure errors. This results in solving a 9 9 system which reduces to two 3 3 systems easily invertible. Comparisons with other estimates based on a Petrov-Galerkin solution are used in our analysis, which shows that it provides a reasonable approximation of the actual discretization error. Numerical experiments clearly show the eciency of such an estimate in the solution of self adaptive mesh re nement procedures. Key words. Mixed nite element methods, Stokes equations, a posteriori error estimates, mesh adaptation, mini-element formulation, Petrov-Galerkin formulation. AMS subject classi cations. 65F10, 65N20, 65N30.
1. Introduction. The need for accurate solutions of large scale problems (in particular) in Computational Fluid Dynamics has made the use of adaptive, automatic re-meshing very attractive for nite element computations of approximations to solution of partial di erential equations [2]. A posteriori error estimates/estimators were introduced in order to provide an information about the local and global quality of the computed nite element solution. They allow the automatic determination of the zones in the mesh which require some re nement or unre nement. In this paper we present and analyze an error estimate for the Stokes problem [13], which plays a center role in the solution of more complicated problems arising in particular in Computational Fluid Dynamics [12] [18] [13] [15] [14]. In section 2, we introduce the equations and the notations used. These equations were solved using a two-level iterative scheme applied to a mini-element discretization of the corresponding variational formulation [6]. Numerous a priori and a posteriori error estimates for elliptic problems [3], for inde nite problems like the Stokes equations [16] or more general problems [2] have already been derived. In particular our new estimate can be viewed as a simpli cation of the one presented in [16]. Because a direct analysis of this estimate was quite dicult, we compared the mini-element formulation [1] [7] [8] with the method proposed by T.J.R.Hughes et al. [10] for solving the Stokes problem, and performed our analysis on the Petrov-Galerkin scheme. Section 3 contains a presentation of the error estimate, which is based on the solution of a local Stokes problem in each element. This estimate is shown to be both a global upper and lower bound of the discretization error, by comparing it to an estimate derived from the Petrov-Galerkin approach. Received by the editors February 16, 1990; accepted for publication (in revised form) August 24, 1990. y Department of Mathematics, University of California at San Diego, La Jolla, California 92093. The work of this author was supported by the Oce of Naval Research under contract N00014-89J1440, by Dassault Aviation, 78 quai Marcel Dassault, 92214 St Cloud, France and by Direction des Recherches Etudes et Techniques, 26 boulevard Victor, 75996 Paris Armees, France. z Department of Mathematics, Arizona State University, Tempe, Arizona 85282. The work of this author was supported by Dassault Aviation, 78 quai Marcel Dassault, 92214 St Cloud, France and by Direction des Recherches Etudes et Techniques, 26 boulevard Victor, 75996 Paris Armees, France. 1
In section 4, we test this estimate on several classical problems and demonstrate its eciency in grid adaptation. 2. The Stokes Equations. In this section we consider a mixed nite element approximation of the following Stokes equations [9]: Find u (velocity eld, 2 components) 2 (H1 ( ))2 (the usual Sobolev space) and p (pressure eld) 2 L2 ( ) (the usual Lebesgue space) such that (2.1)
8 > < ? u + rp ru > : u
= f = 0 = g
in
in
on @
in a bounded domain R2 . The function f is aRsmooth function on R2 , g is linear and satis es the compatibility condition @ g nds = 0. Furthermore Rpiecewise p d
is assumed to be 0. The constant is a viscosity parameter.
We de ne the spaces (2.2) (2.3) (2.4)
Hg1 ( ) = fu 2 (H1 ( ))Z2 ; u = g on @ g L20 ( ) = fp 2 L2 ( ); p d = 0g
H g = Hg1 ( ) L20 ( )
and the two bilinear forms
Z
(2.5)
a(u; v) =
(2.6)
b(u; p) = ?
Z
rurvd
u; v 2 (H1 ( ))2
p r:ud
u 2 (H1 ( ))2 ; p 2 L2 ( )
(; ) will denote the L2 inner product associated with the norm k k. We de ne the energy norm kj(; )kj by kj(u; p)kj2 = a(u; u) + 1 kpk2 u 2 H1 ; p 2 L2 (2.7)
Then a classical variational formulation of the system of equations (2.1) reads Find (u; p) 2 Hg such that (2.8)
a(u; v) + b(v; p) b(u; q)
= (f; v) v 2 (H01 ( ))2 = 0 q 2 L20 ( )
It satis es the following \inf-sup" condition: (2.9)
inf2
p 2 L ( ) p 6= 0
b(u; p) > 0 sup kr uk:kpk 1 u 2 (H1 ( ))2 u 6= 0
which guarantees existence and uniqueness of a solution (u; p) of (2.8) (or (2.1)). Let T be a triangulation of such that any two triangles in T share at most a vertex or an edge. For 2 T let h be the diameter of and E the set of (three) 2
edges on @ . Let h = max h . The set E contains all interior edges and for e 2 E 2T we denote by he the length of e. We suppose also that the triangulation T satis es a minimal angle condition, i.e. the smallest angle in triangle 2 T is bounded away from zero by some constant independent of h. This is equivalent to
h? 1 emin h C1 > 0 and h? 1 emax h C2 2 T 2E e 2E e Furthermore, for ? = e, E , or some subset of E , we de ne the inner product and
(2.10)
associated norm
Z
hu; vi? = uvds = ?
XZ e2? e
uvds
kuk? = (hu; ui? )1=2 Let C 0 be the space of continuous functions over T . Let i = i ( ); i = 1; 3 be the barycentric coordinates (linear nodal basis functions) in the triangle . In order to facilitate the introduction of local function spaces and inner products we will need in our analysis, we consider the following spaces of piecewise H1 functions
HT =
(2.11) and the spaces (2.12)
L=
(2.13)
K=
(2.14)
B=
(2.15) (2.16) (2.17)
Y 2T
Y
2T
Y
Y
2T
L = K = B =
H1 ( ) = fu; uj 2 H1 ( ); 2 T g
Y 2T
Y
2T
Y
2T 2T X = HT \ L \ C 0 X = (X B)2 Y = L20 \ L \ C 0
spanf i( ); 1 i 3; 2 T g spanf i( ) j ( ); 1 i < j 3; 2 T g spanf 1( ) 2 ( ) 3 ( ); 2 T g
Q
and set Q = L K for 2 T and Q = 2T Q . L is the space of piecewise linear functions and Q the space of piecewise quadratic functions on T . The elements of B are referred to as bubble functions. For u; v 2 (HT )2 and p 2 L20 the forms a(; ) and b(; ) are interpreted as (2.18)
a(u; v) =
(2.19)
b(u; p) =
X
2T
X
a(u; v) = b(u; p) = ?
XZ
2T Z
X
rurvd p r ud
2T are contained in L2 the L2 -inner product on X and Y 2T
Note that since X and Y is the usual L2 -inner product. Problem (2.8) can be solved using a good choice of spaces for u and p. 3
For example, the mini-element discretization [1] of the system (2.1) is given by: Find (uh ; ph ) 2 Xg Y such that
a(u ; v) h
(2.20)
b(uh ; q)
+ b(v; ph ) = (f ; v) = 0
for all (v ; q) 2 X0 Y . This formulation also satis es an inequality of the type (2.9) [9], (2.21)
inf
(u; p) 2 X Y (u; p) 6= 0
a(u; v) + b(u; q) + b(v; p) 1 kj(u; p)kj:kj(v; q)kj (v; q) 2 X Y (v; q) = 6 0 sup
for a constant 1 > 0 which is independent of the mesh size h (both constants 1 in (2.21) and (2.9) are denoted with the same symbol since no confusion can arise). This condition implies the unique solvability of the system (2.20). The decomposition uh = uh;l + uh;b , with uh;l 2 (L \ C 0)2 and uh;b 2 B2, is unique. A static condensation of the bubble unknowns uh;b yields, for v 2 X and q 2 Y:
8 a(uh;l; v) + b(v; ph) < (2.22): b(uh;l ; q) ? X 1 (rph ; rq) 2T 3600
= (f ; v) X 1 = ? 60 (f ; b ) rq 2T
Here = j 1j (r b ( ); r b ( )) (j j represents the area of the triangle ) and 2 b ( ) = 1 ( ) 2 ( ) 3 ( ) in triangle (bubble function). is of order O(h? ) in the sense that there exist two positive constants C3 and C4 depending on the minimal angle in the triangulation such that
C3 h2 ?1 C4 h2 Using elementwise integration by parts on (2.22), a simple calculation yields, for (u; p) solution of the system (2.1) and (v ; p) in HT HT L20
8 @ (u ? uh;l) > ; [v]J i a ( u ? u ; v ) + b ( v ; p ? p ) = ( r ; v ) ? h > h < @@nu A +h(p ? ph )n; [v]J iE + h @ nh;l ; [v]A iE > > J : h;l
(2.23)
b(u ? uh;l ; q) = (s; q)
where the vector n = (nx ; ny ) is the unit normal to a triangle, and average values ([]A ) and jump ([]J ) of a function v across an edge e are de ned edgewise by: 4
CHHH HH C HH C : H C out in C n bb C bb C e bC
[v]A = 21 (vin + vout ) and [v]J = vout ? vin
Here in denotes the current triangle. The quantities r and s are the residuals of (2.1) de ned as
(r
= f ? rph s = r uh;l
(2.24)
The velocity term does not appear explicitly in the residual r since uh;l = 0. However, an integration by part of (r; v) yields : (r; v) =
X
2T
(r; v)
= (f ; v) ? a(uh;l ; v) ? b(v; ph ) ? h @@unh;l ? phn; [v]J iE A @u h;l ? h @ n ; [v ]A iE J which is the computational form we will use in the following. Note also that since the computed pressure is continuous, no jump of the quantity p ? ph appears in (2.23). A classical development of an a posteriori error estimate for this system is dicult, mainly because the mini-element discretization is not a member of a sequence of discretizations of varying degrees of approximation. However, we can take advantage of the similarity between the mini-element formulation and the Petrov-Galerkin method of T. J. R. Hughes et al [10] using continuous piecewise linear interpolation for both pressure and velocity terms. Let G = L \ C 0 or Q \ C 0 , and (u ; pG) (u ; pL) when G = L \ C 0 and (u ; pG ) (u ; pQ ) when G = Q \ C 0 . Then the Petrov-Galerkin formulation reads: (2.25)
G
(2.26)
L
G
Q
8 a(u ; v) + b(v; p ) = (f ; v) v 2 G2 > G G > X 1 > < b(uG; q) ? 3600 ((rpG; rq) ? (uG; rq)) 2T > X 1 > (f ; rq) q 2 G = ? > : 2T 3600
It is well known that the use of either piecewise linear or piecewise quadratic velocity and pressure terms yields a stable formulation, provided that the coecient (3600 )?1 is small enough. We now show that the matrix of the resulting system is 5
non singular (equivalent to the fact that the formulation satis es an inf-sup condition). To see this we will need the following lemma: Lemma 2.1.
kuk2 2160 kruk2 1944 105 2 kuk2 u 2 Q p Also kr(ru)k 2160 15 kuk u 2 K
For 2 T
Proof. Let f1; 2; 3g be a numbering of the vertices in the triangle and let i be the angle at vertex i; 1 i 3. We have Q = spanf1; 2; 3 ; 4 1 2 ; 4 2 3 ; 4 3 1 g. So let
u = u1 + u2 2 + u3 3 + 4u1b 2 3 + 4u2b 3 1 + 4u3b 1 2 and set xt = (u2 ; u3; u1b ; u2b ; u3b ). De ne the 5 5 matrices M = 00 D0 and N = BAt B C where the 2 2 matrix A, the 2 3 matrix B and the 3 3 matrices C and D are
respectively de ned by c + c ?c 1 A = 2 1?c1 3 c1 +1c2
0
c ?c3 ?c2 C = 34 @ ?c3 c ?c1 ?c2 ?c1 c
B = ? 32
?c c + c ?c 3 1 3 1 ?c2
?c1
c1 + c2
0c 1 1 2880 @ c1 A ? c c c t A D= 2 1 2 3 c c3
Here ci = cotg(i ) (1 i 3) and c = c1 + c2 + c3 . Note also that c = 180 j j. Since
kuk2 = xt M x and kruk2 = xt N x we have the inequality
kuk2 1 kruk2 with 1 being the largest (positive) root of det( 1 M ? N ) = 0. The matrix M is positive semi de nite and of rank 1 and N is positive de nite, so that this equation has root 0 of multiplicity 4 and a positive root, which is 1 . A direct calculation (add the last three rows in the determinant) shows that 1 = 2160, i.e. this value is independent of the geometry of the triangle . We have also Q = spanf 1; 2 ; 3 ; 4 1 2 ; 4 2 3 ; 4 3 1 g. So let u = u1 1 + u2 2 + u3 3 + 4u1b 2 3 + 4u2b 3 1 + 4u3b 1 2 and set yt = (u1 ; u2 ; u3; u1b ; u2b ; u3b ). Now we have kuk2 = yt S y and kruk2 = yt N 0 y 6
where
0 30 BB 15 B 15 j j S = 180 B BB 12 @ 24
1
15 15 12 24 24 30 15 24 12 24 C C zt 15 30 24 24 12 C 0 C N = z N 24 24 32 16 16 C C 12 24 16 32 16 A 24 24 12 16 16 32 and = c2 +2 c3 , zt = 61 (?3c3 ; ?3c2; ?4(c2 + c3 ); 4c3 ; 4c2 ). The matrix S is positive de nite while N 0 is now positive semi-de nite, which implies that the maximal eigenvalue of N 0 is smaller than the trace Tr(N 0 ) of N 0 . Thus yt N 0 y Tr(N 0 )yt y = 5c yt y p 901 j j ' 1:093 10?2j j > 10?2j j, 31 ? Since the minimal eigenvalue of S is = 90 we get: yt N 0 y 500c y t S y = 9 104 yt S y
j j
which was to be shown. Finally, for u 2 K , let u = 4u1b 2 3 + 4u2b 3 1 + 4u3b 1 2 . Then, in the two-dimensional plane (x; y): kr(ru)k2 = j j (u2xx + 2u2xy + u2yy ) = 32j j 2( +(
3 X
i=1
192j j
3 X
i=1
uib i+1;x i+2;x )2 + 2(
3 X
i=1
uib i+1;y i+2;y )2
uib ( i+1;x i+2;y + i+1;y i+2;x ))2
3 X
i=1
!
u2ib kr i+1 k2 kr i+2 k2
768 902j j 2
3 X
i=1
u2ib
192 902 45 2 kuk2 where the indices x; y; xx; xy and yy denote partial derivatives, and the indices i + 1 and i + 2 represent the modulo 3 values of i + 1 and i + 2 in the set f1; 2; 3g. These inequalities also apply to u 2 Qk (or u 2 K k ), k > 1, by simply adding
the results on each component. Let
S = f(u; p) 2 (Q \ C 0 )3 ; u = g on @ ; Introducing the (nonsymmetric) bilinear form
B ((u; p); (v ; q)) = a(u; v)+ b(v; p) ? b(u; q)+ 7
X
Z
pd = 0g
1 ((rp; rq) ? (u; rq) ) 2T 3600
and the semi-norm k k de ned by
k(u; p)k2 kruk2 +
X
1 krpk2 3600 2T
we have then
Theorem 2.2. For (u; p) 2 Q Q Q we have
B ((u; p); (u; p)) !k(u; p)k2 Proof. It is similar to the proof of the coercivity of the very same bilinear form in [8] for a special choice of the stability constant (3600 )?1 for 2 T . A direct application of Lemma 2.1 yields X ? B ((u; p); (u; p)) = a(u; u) + 36001 krpk2 ? (u; rp) 2T X 1 ? 2 p a(u; u) + 3600 krpk ? 2160 kruk krpk 2T
! p X 1 3 2 2 (1 ? p ) kruk + 3600 krpk 2 5 2T p (2.27) = (1 ? p3 )k(u; p)k2 2 5 p and the theorem holds with ! = 1 ? 0:5 0:6 ' 0:612. Then, since the semi-norm is a norm on subspaces of Q Q Q (like K K K) the form B (; ) is positive de nite on these same subspaces (and in particular on K K K). Note that kj kj and k k de ne two norms on K K K which are equivalent, i.e.
1 k(u; p)k kj(u; p)kj 2 k(u; p)k (u; p) 2 K K K for some positive constants 1 and 2 . The rst inequality results from a direct application of Lemma 2.1, and is more generally valid on Q Q Q ( 1 = 25), while the second one is a consequence of the minimal angle condition expressed by the inequality
(2.28)
k k pc5 kr k for 2 K
It is true in particular in K K K (but not all of Q Q Q). As an indication, c5 < 2 when the minimal angle in the triangulation is greater than 1 degree, which gives 2 < 120 (c5 < 0:4 for a minimal angle greater than 5 degrees, or 2 < 24). Another important case where the semi-norms are equivalent is the following Lemma 2.3. Let (u; p) 2 Q Q Q such that
a(u; v) + b(v; p) = 0 all v 2 L L 8
Then there exist two positive constants 3 and 4 such that
3 k(u; p)k kj(u; p)kj 4 k(u; p)k Proof. Our proof is based on the duality argument given in [7]. Consider the dual problem
8 > < ?w + r rw > : w
= 0 = p = 0
in
in
on @
so that the regularity inequality
krwk + kk c6 kpk holds. We have then
kpk2 = (p; r w) = (p; r v) + (p; r (w ? v)) where v is the linear interpolator of w at the vertices of the triangulation so that
kw ? vk pc7 krwk
and
krvk krwk + kr(v ? w)k c8 krwk for all triangles ( = O(h? 2 )).
Integrating the second term by parts and using the hypothesis on (u; p) for the rst term, we get
kpk2 = a(u; v) ? (rp; wX ? v) kruk krvk + krpk kw ? vk 2T
c8 kruk krwk + c7
X
2T
p1 krpk krwk
c8 kruk krwk + c7 krwk c6 c8 kruk kpk + c6 c7 kpk so that
kpk c6 c8 kruk + c6 c7 9
X krpk2 !1=2
2T
X krpk2 !1=2 2T
X krpk2 !1=2 2T
Consequently
?1 kpk2 2c26c28 kruk2 + 2c26 c27
X krpk2 2T
and nally
kj(u; p)kj 4 k(u; p)k Conversely
1 krpk2 25 X kpk2 = 25kpk2 2T 2T 3600
X
which implies
X
1 krpk2 2T 3600 2 kruk2 + 25 kpk 25kj(u; p)kj2
k(u; p)k2 = kruk2 +
Note that the second part of the proof remains valid for any (u; p) 2 Q Q Q. In order to show the next result (continuity property) on the bilinear form B , we shall use the inequality
kr uk2 2 kruk2 u 2 HT HT Lemma 2.4.
jB ((u; p); (v ; q))j 27kj(u; p)kj kj(v; q)kj for all (u; p) and (v ; q) in Q3 .
(2.29)
Proof. A straightforward calculation leads to
jB ((u; p); (v; q))j kruk krvk + kpk kr v k + kqk kr uk X pk krqk X kuk krqk + kr3600 + 2T 3600 2T p p kruk krvk + 2kruk kqk + 2krvk kpk X p X + 25 k p k k q k + 15 kruk kqk 2T 2T p p kruk krvk + 2kruk kqk + 2krvk kpk p 15kruk kqk k p k k q k + + 25 ? ? 18 kruk2 + 27 ?1kpk2 1=2 2 krvk2 + 27 ?1kqk2 1=2 27kj(u; p)kj kj(v; q)kj 10
The mini-element form after static condensation (2.22) and the Petrov-Galerkin form using piecewise linear velocities and pressure (G = L in (2.26)) are nearly identical: the sti ness matrices are identical and the right-hand sides di er by only a small quantity since 60(1; b) = j j so that j(f ; rq) ? 60(f ; b ) rqj = (f ? f ; (1 ? 60 b)rq) c9 h krf k krqk where f is the mean value of f in triangle and c9 is a constant independent of . Consequently, we have Theorem 2.5.
2 kj(uL ; pL) ? (uh;l ; ph )kj c10ph krf k
Proof. By Theorem 2.2 we have
!k(uL ; pL) ? (uh;l ; ph)k2 B ((uL ; pL) ? (uh;l ; ph ); (uL; pL ) ? (uh;l ; ph)) X ((f ; r(pL ? ph)) ? 60(f ; br(pL ? ph)) ) = 3600 2T c9 X h krf k :kr(p ? p )k 3600 L h 2T p 2X p 1 krf k :kr(pL ? ph )k cp9 C4 h 3600 2T 3600 p 2 cp9 C4 h krf k k(uL ; pL) ? (uh ; ph)k 3600 The linear part (uh;l ; ph ) of the computed solution (uh ; ph ) in the mini-element formulation is certainly a good approximation of the solution (u; p) of the initial problem, and that the bubble part of this solution is only introduced for stability reasons and does not improve the approximation of the velocity and pressure terms, as was pointed out in [16]. In the next section we take advantage of this special similarity between these two formulations in order to carry out the analysis of our error estimate.
3. An a posteriori error estimate based on the solution of a local Stokes problem. The goal of this section is to de ne an estimation of the discretization errors (e; ) = (u; p) ? (u ; pL) and (e0 ; 0 ) = (u; p) ? (u ; ph) for the Petrov-Galerkin L
h;l
and Mini-Element formulations respectively. The estimates are based on the solution of local Stokes problems (i.e. de ned in each element). The analysis and derivation of these small problems extend the work done by Bank and Weiser for elliptic problems [5]. 11
The main results presented here are Theorem 3.6 and Corollary 3.11, which show that the estimates provide a good global assessment of the discretization error, under reasonable assumptions. Let (e ; Q ) = (u ; pQ ) ? (u ; pL ). An error analysis on the system (2.26) similar to the one on (2.20) leads to the following system of equations: Q
Q
L
8 @e > a(e; v) + b(v; ) = (r; v) ? h @ n ; [v ]J i > > A@ u > > +hn; [v]J i + h @ n ; [v]A i > < J X (3.1) > 1 b(e; q) ? 3600 (r; rq) ? (e; rq) = (s; q) > 2T > X > > ? 36001 (r; rq) : L
2T
Equations (3.1) are valid for all (v ; q) 2 (H (T ))3 . This is equivalent to (3.2)
B((e; ); (v; q)) F ((v; q))
(v ; q) 2 (H (T ))3 :
= L((v ; q)) + J ((e; ); (v; q))
with the forms L() and J () de ned by
X L((v ; q)) = (r; v) ? (s; q) + h @@un ; [v ]A i + 36001 (r; rq) J 2T @e J ((e; ); (v ; q)) = ? h @ n ; [v]J i + hn; [v]J i L
A
and the residuals r and s have the same meaning as in (2.24) with (uh;l ; ph ) replaced by (u ; pL). Note that J ((v; q)) = 0 when v is continuous. From equations (3.2) we have for (v; q) 2 (L \ C 0 )3 L
(3.3)
B ((e ; Q ); (v ; q)) = B ((u ; pQ ); (v ; q)) ? B ((u ; pL); (v; q)) = 0 Q
Q
L
and in particular
a(eQ ; v) + b(v; Q ) = 0 which gives, using Lemma 2.3,
3 k(eQ ; Q )k kj(eQ ; Q )kj 4 k(eQ ; Q )k Note that the classical relation (3.4) also holds.
B ((e; ); (v; q)) = B ((u; p); (v; q)) ? B ((u ; pL ); (v; q)) = 0 L
12
For (v; q) 2 (Q \ C 0 )3 we also have
B ((e ; Q ); (v ; q)) = B ((u ; pQ ); (v ; q)) ? B ((u ; pL); (v ; q)) = B ((u; p); (v ; q)) ? B ((u ; pL); (v; q)) = B ((e; ); (v; q)) = L((v; q)) We now de ne the local error estimate (e; ) 2 K3 by (3.6) B ((e; ); (v; q)) = L((v; q)) (v ; q) 2 K3 ( ) Since the bilinear B (; ) form is positive de nite on K3 , this estimate is well de ned, on each triangle 2 T . The velocity components as well as the pressure component are described in each triangle by quadratic bump functions, which schematically
(3.5)
Q
Q
L
L
correspond to the following degrees of freedom:
y y y
A A A A A A
degrees of freedom for velocity and pressure errors Equation (3.6) is a 9 9 system to be solved in each triangle. Note that using an integration by parts on (r; v) similar to the one performed to get (2.25) the right-hand side L((v; q)) takes the form
L((v; q)) = (f ; v) ? a(uL; v) ? b(v ; pL) + b(uL; q) + h @@un ? pL n; vi@ + 36001 (f ? rpL ; rq) L
A
In the sequel of the paper we will need some notion of convergence of the nite element solutions (u ; pL ) and (u ; pQ ) to the weak solution (u; p) of (2.1). In particular, for 2 T , we make the following saturation assumption L
(3.7)
Q
@ (u ? u
2 1 = 2
kj(u; p) ? (u ; pQ )kj + he @ n
Q
Q
2 kj(u; p) ? (u ; pL)kj2
) 2 j(p ? p )nj2 !1=2
2 Q
+
A
E
L
where = o(1), which implies in particular (3.8)
(1 ? )kj(e; )kj kj(e ; Q )kj (1 + )kj(e; )kj Q
This is not a very strong condition since (u ; pQ ) is an approximation to (u; p) in a higher degree polynomial space than (u ; pL). It supposes however that the solution is more than H1 -regular. Q
L
13
We now present a few basic estimates, which will be used in the next part of this section. Their proofs make use of the following inequalities, valid for 2 T and e 2 E : ?1=2
@kvv
ke c11h kvk v 2 Q( ) Q( )
c12h? 1=2krvk v 2 Q( ) Q( ) @n e
(3.9) (3.10)
These inequalities are not satis ed for general functions v; they remain valid however for polynomial functions, the constants c11 and c12 depending then on the degree of the polynomial. Lemma 3.1.
kvk2 c13 h2 a(v; v)
v 2 K2 ( )
Proof. It is a simple reformulation of inequality (2.28) where = O(h? 2 ). Lemma 3.2.
p
h?1=2 v
c a(v; v)1=2
e @ 14
v 2 K2 ( )
Proof. for e 2 @ , using inequality (3.9) and Lemma 3.1, we have
kvk2e c211 h ?1 kvk2 c13 c211 h a(v; v) so that
2
h?e 1=2 v
@ 3C1?1 c211 c13 a(v; v)
Lemma 3.3. Assume (3.7) holds. Then
h1=2 @ e 2 + jnj2 !1=2
c kj(e; )kj
15 @n A
e E
Proof. The proof is similar to the proof of Lemma 3.1 in [5] with only a few di erences due to the pressure term. For e 2 E such that e = in \ out
2
h1=2 @ (u ? u ) 2 + j(pQ ? pL)nj2 !1=2
@n A
e
e
! 1=2
2 2 2
=
h1e=2 @@en + jQnj
A e Q
L
Q
14
(
2 k( n)j k2
@ e
2) @ e Q e
he @ n j + + @ n j
e e 2 c Q
in
in
Q
out
1 kreQ k2 ) + 11 h?1 kQ k2 he c212 (h? in1 kreQ k2 in + h? out out in in C2 (c212 + c211 )(kj(eQ ; Q )kj2 in + kj(eQ; Q )kj2 out )
Therefore, summing over all edges, we get
h1=2 @ (u ? u @n
e Q
(3.11)
L
!1=2
2 ) 2 + j(pQ ? pL )nj2
3C2(c212 + c211 )kj(eQ; Q)kj2 A
E 3C2 (c212 + c211 )(1 + )2 kj(e; )kj2
Finally, combining (3.7), (3.11) and using the triangular inequality we get
h1=2 @ e 2 + jnj2 !1=2
@n A
e
E
2 j(p ? p )nj2 j nj2 !1=2
2
@ ( u ? u ) @ e Q
h1e=2 2 @ n Q + 2 @ nQ + 2 + 2 Q
E A A
! p @ (u ? uQ ) 2 j(p ? pQ )nj2 1=2
1 = 2
he 2 @ n +
A E
!1=2
2 2 p
+
h1e=2 2 @@enQ + jQnj
E A q 2 2 p 2kj(e; )kj + (1 + ) 6C2 (c12 + c11 )kj(e; )kj Lemma 3.4. The \continuity" inequality
jB ((e; ); (v ; q))j c16 kj(e; )kj kj(v; q)kj holds for all (v; q) 2 K K K.
Proof. We bound successively each term in the de nition of B (; ). The rst three terms are easily bounded, namely:
ja(e; v)j krek krvk p jb(v; )j kk kr v k p2kk krvk and jb(e; q)j kr ek kqk 2krek kqk while the two remaining terms need to be integrated by parts; for 2 T we have: (r; rq) = ?(; q) + hn; rqi@ 15
On one hand we get
X p p (; q) 15 X kk kqk 15 kk kqk 2T 3600 2T
and on the other hand the boundary term gives
X hn; rqi@ = ? 1 Xh; 1 @q ie 2T 3600 3600 e2E @ n J
1 X kk
1 @q
3600 e2E e @ n J e
Now, for e = in \ out 2 E :
1 @q
2 2 2 ?1 ?2
@ n J e 2c12 (h krqk in
in
in
1 ?2 krq k2 ) + h? out out out
18 104 c212 C4 (h in kqk2 in + h out kqk2 out ) 18 104 c212 C4 C1?1 he (kqk2 in + kqk2 out ) Note that the same inequality holds with the jump of 1 @@qn replaced by its average value. Consequently
X c qC C ?1 X !1=2 hn; rqi@ 12 p4 1 2 he kke 2T 3600 6 2 e 2 E
jj
kqk c17
h1e=2 p
p E q ?1
3
X 2T
kqk2
!1=2
c12 C4 C
p 1 . where c17 = 2 6 Likewise, an integration by parts on the last term yields (e; rq) = ?(re; r(rq)) + h @@ne ; rqi@ and by Lemma 2.1:
p X (re; r(rq)) 3p 3 krek kqk 2T 3600 5
As for the boundary term involving the velocity, a calculation comparable to the one performed above for the pressure implies X @e h @n ; rqi@ = ? 1 X h @ e ; 1 rq i 2T 3600 3600 e2E @ n A J @ e 1 +h @ n ; rq i A J 16
h1e=2
@@ne
+
@@ne
(kqk2 in + kqk2 out )1=2 3 e2E A e J e X 1=2
@ e
2 c17 3he @ n (kqk in + kqk2 out )1=2
c17 p
X
@ e A
e kqk
p 3c17
h1e=2 @ n
p A E e2 E
Adding all the terms together nally yields:
p jB ((e; ))j kr e k kr v k + 2 kk krvk p p p 15 3 + 2 krek kqk + kk kqk + p 3 krek kqk
jj
kqk
p 5 @ e
kqk +c17
h1e=2 p
p + 3c17
h1e=2 @ n
p E A E 1 = 2 2 2 42 krek2 + 17 kk krv k2 + kqk 5
jj
2
p @ e
2 !!1=2 +9c217
h1e=2 p
+
h1e=2 @ n
E A E 2 2 kjq(v; q)kj (17kj(e; )kj + 9c17 c15 kj(e; )kj2 )1=2 17 + 9c217 c15 kj(e; )kj kj(v; q)kj
Lemma 3.5. Let (v; q)
constant < 1 such that
2 K3 and ( ; ) 2 L3 . Then there exists a positive
kj(v; q)kj p 1 2 kj(v + ; q + )kj 1? Proof. Let 2 T with angles i ; i = 1; 3. De ne d = cos2 1; +cos22; +cos2 3; . From [11] we know that
(rv; r ) krvk kr k and
p
r
(q; ) 415 kqk kk
d ? 43 . If the triangulation satis es a minimal angle condition, p 15 ) < 1. (For example, when the minimal then d < 3 and thus max(max
; 2T 4 angle is 5 degrees we have ' 0:99873). Next we get ( q; ) j((v; q); ( ; ))j (rv ; r ) + kj(v; q)kj kj( ; )kj
with 2 = 12 + 31
17
and nally kj(v + ; q + )kj2 = kj(v; q)kj2 + kj( ; )kj2 + 2((v; q); ( ; )) kj(v; q)kj2 + kj( ; )kj2 ? 2 kj(v; q)kj kj( ; )kj (1 ? 2 )kj(v; q)kj2 The following theorem proves that the estimate (e; ) is a global upper and lower bound of the error (e; ) in the kj kj norm. Theorem 3.6. there exist two positive constants c18 and c19 such that (3.12) c18 kj(e; )kj kj(e; )kj c19 kj(e; )kj Proof. Let (eLQ ; LQ ) and (eQQ ; QQ ) be respectively the linear and (continuous) quadratic parts of (eQ ; Q ), so that (eQ ; Q ) = (eLQ ; LQ ) + (eQQ ; QQ ) Recalling (3.5), (3.6) and using (3.3), we have B ((e; ); (eQQ ; QQ )) = B ((e ; Q ); (eQQ ; QQ )) = B ((e ; Q ); (eQ ; Q )) Now we have, from Theorem 2.2 and Lemmas 2.4 and 3.5: kj(eQ ; Q )kj2 42 k(eQ ; Q )k2 42 !?1 B ((eQ ; Q ); (eQ ; Q )) = 42 !?1 B ((e; ); (eQQ ; QQ )) 27 42!?1 kj(e; )kj kj(eQQ ; QQ )kj 2 ?1 p 4 ! 2 kj(e; )kj kj(eQ ; Q)kj 27 1? and thus 2 ?1 p 14?! 2 kj(e; )kj (1 ? ) kj(e; )kj kj(eQ ; Q )kj 27 Conversely we have kj(e; )kj2 22 k(e; )k2 22 !?1 B ((e; ); (e; )) = 22 !?1 (B ((e; ); (e; )) ? J ((e; ); (e; ))) The jump term on the right-hand side is bounded using Lemmas 3.2 and 3.3: Q
Q
@ e jJ ((e; ); (e; ))j = h @ n ? n; [e]J iE A
@ e 2 j nj2 !1=2
kh?e 1=2p [e]J kE
h1e=2 @ n +
E
A ! p X ?1=2 2 1=2 khe [e]J ke c15 kj(e; )kj e2E p 2c15 c14 kj(e; )kj kj(e; )kj 18
Finally, by Lemma 3.4, we have kj(e; )kj2 22 !?1 (27kj(eQ ; Q )kj kj(e; )kj + c16 kj(e; )kj kj(e; )kj) and hence kj(e; )kj 22 !?1 (27(1 + ) + c16 )kj(e; )kj Let e = u ? uh;l and = p ? ph . We now de ne the local error estimate (~e; ~) 2 K3 ( ) by 0
0
B ((~e; ~); (v; q)) = L ((v; q)) (v ; q) 2 K3 ( )
(3.13) with
0
L ((v ; q)) = (f ; v) ? a(uh;l ; v) + (ph ; r:v ) ? (r:uh;l ; q) + h @@un ? phn; vi@ + 36001 (f ? rph ; rq) 0
h;l
A
Lemma 3.7. For all (v; q) 2 Q3 we have
(3.14)
X L((v; q)) ? L ((v; q)) c20kj(uL; pL) ? (uh;l; ph)kj kj(v; q)kj 0
2T
Proof. We write the di erence between the right-hand sides of (3.6) and (3.13) in triangle 2 T as L((v; q)) ? L ((v ; q)) = ?a(uL ? u ; v) + (pL ? ph; r:v) ?(r:(uL ? u ); q) + h @ (uL@?nu ) ? (pL ? ph )n; vi@ A 1 ? 3600 (r(pL ? ph ); rq) 0
h;l
h;l
h;l
Each of these terms is now bounded in terms of the energy norm of the di erence (uL ? uh;l ; pL ? ph ) between the solution from the Petrov-Galerkin and mini-element formulations, as in the proof of Lemma 3.4: First ja(uL ? u ; v) j p kr(uL ? uh;l )k a(v ; v)1 =2 h;l
j(pL ? ph ; r:v) j kppL ? ph k kr:vk 2kpL ? ph k krvk and similarly
p j(r:(uL ? u ); q) j 2kr(uL ? uh;l )k kqk h;l
The bound on the boundary term is then derived using Lemma 3.2: @ (u ? u )
@ (u ? u )
h L
h1e=2 L h;l
h?e 1=2v
; v i @ @n @n @ A A @ h;l
19
for e 2 @ we have
@ (u ? u )
2 h (
@ (u ? u )
2
@ (u ? u )
h1e=2 L h;l
e
L h;l j
+
L h;l j @n @n @n A e 2 e 2 he c12 h?1 kr(u ? u )k2 + h?1 kr(u ? u )k2
out
L h;l 2 C2 c12 kr(u ? u )k2 L h;l [ out 2
2
L
out
2)
e
h;l out
Thus summing over edges in
0 1
@ (u ? u )
2 X
h1e=2 L h;l
C2c212 @ 3 kr(uL ? uh;l)k2 + 1 kr(uL ? uh;l)k2 A @n 2 2 A @ 0
0
where is a triangle neighbor of (total of 3 neighbors in the 2D case). Similarly 0
jh(pL ? ph)n; vie j kh1e=2(pL ? ph)ke kh?e 1=2v ke C21=2 c11 kpL p?ph k c14 a(v; v)1 =2 Finally, summing over all triangles, we get
X X L((v; q)) ? L ((v; q)) p kr(uL ? uh;l)k a(v; v)1 =2 2T 2T Xp Xp 25 X 0
+ + +
2T
X
2T
X
2T
2kpL ? ph k krvk +
2kr(uL ? uh;l )k kqk + kpL ? ph k kqk 2T
2T
11=2 0 X p C21=2 c12 @ 3 kr(uL ? uh;l )k2 + 1 kr(uL ? uh;l )k2 A c14 a(v; v)1 =2 2
2
0
0
3C21=2 c14 c11 kpL ? ph k a(v; v)1 =2
X 2T
!1=2
(a(v ; v) + ?1 kqk2 )
(3 + 2C2 c214 c212 )
X 2T
kr(uL ? uh;l )k2 + (627 + 9C2 c214 c211 )
c20 kj(v; q)kj kj(uL ; pL) ? (uh;l ; ph )kj
X kpL ? phk2 !1=2 2T
This estimate of the di erence between the two right-hand sides of the systems de ning the two error estimates (~e; ~) and (e; ) allows us to get a bound for the di erence between the estimators themselves, since both systems have the same lefthand side which is positive de nite on the space K3 . Indeed, for (v; q) = (~e; ~) ? (e; ) 2 Q3 in the previous lemma, we have:
kj(~e; ~) ? (e; )kj2 22 k(~e; ~) ? (e; )k2 20
22 !?1 = 22 !?1
X 2T
B ((~e; ~) ? (e; ); (~e; ~) ? (e; ))
X
L((~e; ~) ? (e; )) ? L ((~e; ~) ? (e; )) 0
2T 2 ? 1 2 ! c20 kj(~e; ~) ? (e; )kj kj(uL; pL ) ? (uh;l ; ph )kj
Noting that (uL ; pL) ? (uh;l ; ph ) = (e; ) ? (e ; ) and recalling Theorem 2.5 we thus have the 0
0
Theorem 3.8. 2
kj(~e; ~) ? (e; )kj c14 kj(e; ) ? (e ; )kj c14pc10 h krf k
(3.15)
0
0
Note that if neither (uL ; pL) nor (uh;l ; ph ) is the exact solution to the original Stokes system, then both terms (e; ) and (e ; ) are of order O(h) and their di erence is of order O(h2 ). However, when for example (uL ; pL) is the exact solution (which is possible for properly chosen functions f and g and a particular domain ), (e; ) is zero while (e ; ) is of order O(h2 ). Therefore there is no constant c21 such that 0
0
0
0
c?211 (h)kj(e; )kj kj(e ; )kj c21 (h)kj(e; )kj 0
0
Finally the best result we can get reads Theorem 3.9. There exist two positive constants c22 and c23 such that (3.16)
c22 kj(e ; )kj ? O(h2 ) kj(~e; ~)kj c23 kj(e ; )kj + O(h2 ) 0
0
0
0
Proof. It is a direct consequence of Theorems 3.6 and 3.8. In the last step of our analysis we simplify the system (3.13). We introduce the nal estimate (e ; ) 2 K3 de ned in triangle by 00
00
8 a(e ; v) + b(v; ) > > < (3.17) > > : b(e ; q) ? (r ; rq) 00
00
00
00
3600
= (f ; v) ? a(uh;l ; v) + (ph ; r:v) +h @@unh;l ? phn; vi@ A = (r:uh;l ; q) ? 36001 (f ? rph ; rq)
for (v; q) 2 K3 . Note that the right-hand side has not changed and that only the Laplace term has been removed from the left-hand side, when compared to the system (3.13). The corresponding 9 9 matrix M of this system takes the form (3.18)
0 A M = @ 0
1
t B ;x t A A B ;y B ;x B ;y ?S
0
21
where the 3 3 matrices A , B ;x and B ;y are respectively de ned by (3.19) (3.20) (3.21)
Z
r ib ( )r ib ( )d 1 i; j 3 ib ( ) d 1 i; j 3 (B ;x )i;j = ? jb ( ) @ @x Z ib ( ) d 1 i; j 3 (B ;y )i;j = ? jb ( ) @ @y (A )i;j =
Z
( 1b = 4 2 3 and cyclically) and S = 36001 2 A . In the case of a boundary triangle, the equations in the system (3.17) associated with the corresponding edge(s) are replaced by a scaled version of the Dirichlet boundary condition e00 = u ? uh;l , where uh;l represents the linearly interpolated value of uh;l at the midpoint of the edge. Thus all elements but the diagonal terms in the corresponding rows of M are zeroed out. Since K( ) does not contain the constant functions, the matrix A is symmetric positive de nite. If is an interior triangle, then the Schur complement C = C 0 + S t + B ;y A?1 B t , is therefore well de ned, and because of M , with C 0 B ;xA? 1 B ;x ;y S is positive de nite, C is also positive de nite. The matrix C 0 can in fact be shown to be independent of the geometry of the triangle , even though the matrices A , B ;x and B ;y are not (see [17]). In a boundary element, the matrix C 0 does now depend on the geometry of the element, but is still positive semi-de nite, so that C is non-singular. Hence the matrix M is non singular since 0 A 0 Bt 1 0 A 0 0 1 0 I 0 A?1Bt 1 ;x ;x t A=@ 0 t A A 0 A @ 0 I A? 1 B ;y (3.22)@ 0 A B ;y B ;x B ;y ?C B ;x B ;y ?S 0 0 I so that the system (3.17) has a unique solution (the non-singularity of the Schur complement C (and hence its positive de niteness), together with the non-singularity of A , is equivalent to the stability of the discretization of the error, i.e. the discretization will satisfy a local Babuska-Brezzi type condition). In our nal theorem we compare this last estimate with the discretization error resulting from the mini-element formulation. Theorem 3.10. There exist two positive constants c24 and c25 such that c24 kj(e ; )kj kj(~e; ~)kj c25 kj(e ; )kj 00
Proof.
00
00
kj(e ; )kj2 22 k(e ; )k2 00
00
00
00
00
X
a(e ; e ) + 36001 kr k2 2T = 22 L0 ((e ; )) = 22 B ((~e; ~); (e ; )) 27 22 kj(~e; ~)kj kj(e ; )kj = 22
00
00
00
00
00
00
00
00
22
00
!
Conversely,
kj(~e; ~)kj2 22 k(~e; ~)k2 22 !?1B ((~e; ~); (~e; ~)) = 22 !?1L0 ((~e; ~))
X
a(e ; e~) + 36001 (r ; r~) + b(~e; ) ? b(e ; ~) 2T p p 22 !?1 k(e ; )k k(~e; ~)k + 2k k kre~k + 2k~k kre k p 22 !?1( 2 + 1?2 )kj(e ; )kj kj(~e; ~)kj
22 !?1
00
00
00
00
00
00
00
00
!
00
00
This equivalence between the two estimate is true locally as well, since all inequalities in the proof remain valid in any triangle 2 T . Corollary 3.11. There exist positive constants c26 , c27 , c28 and c29 such that
c26 kj(e ; )kj ? c27 h2 kj(e ; )kj c28 kj(e ; )kj + c29 h2 0
0
00
00
0
0
This last result means that the error estimate (e ; ) is a reasonable global estimate of the discretization error resulting from the use of the mini-element. No local inequality is proved. However, as will be seen in section 4, this estimate seems to be also a good local indicator of the size of the error and becomes an ecient tool in the adaptation process of the mesh to the solution of the original problem (2.1). 00
00
4. Numerical Results. The system (3.17) is solved in each triangle using the decomposition (3.22). Thus we need to solve eight 3 3 systems with the matrix A t and A?1 B t ), and then another (among which six are necessary to compute A? 1 B ;x ;x 3 3 system with C to get the pressure error and nally the velocity error. We present three numerical examples which demonstrate the eciency of our error estimator in controlling the mesh adaptation process and also in estimating the discretization error. 4.1. Driven Cavity. We consider rst the standard case in CFD of a driven cavity in a unit square. In this example = 0:01. The right hand side f in (2.1) is set to zero and boundary conditions are zero except on the upper side of the square where the velocity is tangential and has value 1: 23
ut = (1; 0) -
ut = (0; 0)
ut = (0; 0) ut = (0; 0)
This domain is rst triangulated into NT = 8 triangles, to form the level 1 grid (Fig. 4.1(a)). After computing the solution of (2.1) on this coarse grid, the mesh is either uniformly (Fig. 4.1(b)) or adaptively (Fig. 4.1(c)) re ned. Each triangle is subdivided into 4 smaller triangles by joining the midpoints of the three sides (regular re nement). \Green" triangles are added to make the mesh an admissible triangulation. Re nement stops when the total number of vertices NV in the mesh reaches a preassigned target value. For further details on the re nement strategy the reader can refer to [4]. We can notice that the adaptive strategy created a lot of triangles (up to 21 levels) in the two upper corners of the cavity, where two singularities arise due to the discontinuities in the boundary conditions.
(a)
(b)
(c)
Fig. 4.1. (a) Initial grid with N V = 9 vertices and N T = 8 triangles; (b) uniformly re ned grid with N V = 289 vertices and N T = 512 triangles (4 levels) and (c) adaptively re ned grid with N V = 282 vertices and N T = 496 triangles (10 levels).
The resulting velocity eld is plotted on Figure 4.2. Each arrow represents the velocity at the inner vertices, with a streamline direction and a length proportional to the norm of the velocity. Although the picture seems to be nicer for the uniform re nement case because of the equal spacing of the vertices, the one corresponding to the adaptive strategy gives a better idea of the real eld, especially at the two top corners: on the left the uid is taken away and is forced back into the cavity on the 24
right side, hence the corresponding negative and positive pressures on Figures 4.3(a) and (b). The capture of the discontinuities is much more e ective on the adapted mesh than on the uniform mesh, and the level of the pressure is about 104 times the level reached on the uniform grid. Error estimates are shown on Figures 4.4(a) and (b). In order to draw continuous functions, errors at the vertices are computed by averaging the estimates on neighboring triangles. 4.2. Backward Facing Step. Another classical case in CFD is the backward facing step. Here the spikes in the pressure are not due to a discontinuity in the boundary conditions but only to a discontinuity in the geometry of the boundary. Boundary conditions are as follows:
ut = ( (y? 2 )(b a?y) ; 0) a
ppp p p pppp
ut = (0; 0) ?
& 6 -
6
ut = (0; 0)
p pp ppp pp p
-
ut = ( y(a8?b y) ; 0)
The initial grid (Fig. 4.5(a)) is uniformly re ned (Fig. 4.5(b), 4 levels of re nement) or adaptively re ned (Fig. 4.5(c), 9 levels of re nement) based on the error estimate computed from the solution ( = 1). In the latter case triangles were added near the re-entrant corner to resolve the discontinuity. Pro les of the velocity eld are very similar this time (Fig. 4.6,4.7) but the drop in pressure is better represented with the adaptive strategy (Fig. 4.8). Note that in Fig. 4.5(c) some re nement was done near the non homogeneous Dirichlet boundaries.
25
(a)
(b)
Fig. 4.2. Velocity eld on the uniform grid (a) and on the adapted grid
(a)
(b).
(b)
Fig. 4.3. Pressure elevation on the uniform grid (a) and on the adapted grid (b); it is negative at the left upper corner of the cavity (suction e ect) and positive at the right upper corner; at these points the level of the pressure is actually much higher on the adapted grid (order of magnitude 102 ) than on the uniform mesh (order of magnitude 100 ).
26
(a)
(b)
(a) and adapted (b) grids. Although the represen(b) is less smooth than (a), the errors are much more uniformly reparted in (b) than in
Fig. 4.4. Estimated error on the uniform
tation (a).
(a) Initial grid with
(b) Uniformly re ned grid with
(c) Adaptively re ned grid with
NV
= 36 vertices and
NT
= 44 triangles.
NV
= 1513 vertices and
NT
= 2816 triangles (4 levels).
NV
= 1517 vertices and
NT
= 2863 triangles (9 levels).
Fig. 4.5.
27
(a) Velocity eld on the uniform grid.
(b) Velocity eld on the adapted grid. Fig. 4.6.
Fig. 4.7. Detail of velocity eld near the re-entrant corner of the step.
(a)
(b)
Fig. 4.8. Pressure elevation on the uniform grid (a) and on the adapted grid (b); note that at the top corner of the step level of the pressure is higher on the adapted grid (order of magnitude 102 ) than on the uniform mesh (order of magnitude 101 ).
28
(a)
(b)
Fig. 4.9. Estimated error on the uniform
29
(a) and adapted (b) grids.
4.3. Disk with a Crack. Finally we test our error estimate on a Stokes ow in a disk of radius 1 with a crack joining the center to the boundary; the right-hand side f is 0 and the boundary conditions are 3 3 3 t u = 2 cos 2 ? cos 2 ; 3 sin 2 ? sin 2 where (r; ) is a polar representation of a point in the disk. The exact solution is then given by
pr 3 3 cos ? cos ; 3 sin ? sin
ut = 3
2
2
2
p = ? p6r cos 2
2
2
and is singular at the end of the crack, i.e. at the center of the disk.
(a)
(b) Fig. 4.10.
(c)
(a) Domain ; (b) velocity ; (c) pressure.
We rst solve our Stokes problem on a coarse grid consisting of NV = 15 vertices and NT = 16 triangles (Fig. 4.10(a)), then re ne either uniformly or adaptively, thus creating two sequences of meshes of increasing and comparable size (or degrees of freedom = d. of f.)(see Tables 4.1 and 4.2) ( = 1). Pressure elevations are plotted on the uniform and adaptive re ned grids (Figures 4.12(a)(b)). The solution was computed on each grid, along with error estimates. During the re nement process, these estimates were computed using an interpolation scheme (one Jacobi sweep) for the values at the new nodes. Consequently their accuracy deteriorates along with the number of re nement steps; however had we used intermediate recalculations to base the computation of the estimates on, we would have gotten a mesh with more levels of re nement around the singularity (hence giving a higher level for the pressure). in that regard the use of interpolated values instead of computed solution values had a grid smoothing e ect. 30
(a)
(b)
(c)
Fig. 4.11. (a) Initial grid with N V = 15 vertices and N T = 16 triangles; (b) uniform mesh with N V = 561 vertices and N T = 1024 triangles (4 levels) and (c) adapted mesh with N V = 563 vertices and N T = 1054 triangles (9 levels).
(a)
(b)
Fig. 4.12. Pressure elevation on the uniform (a) and adapted (b) meshes. The level of pressure is much higher on the adapted grid (order of magnitude 103 ) than on the uniform mesh (order of magnitude 102 ).
31
We can measure the eciency of our estimate both locally and globally: we perform a convergence analysis on the two sequences of (uniform and adapted) meshes and evaluate in particular the \e ectivity ratio" q de ned as the ratio of the exact global error kj(e; )kj to the estimated global error kj(e ; )kj, and report the results in Tables 4.1 and 4.2. The column labeled \digits" gives the estimated number of correct digits in the solution of the discrete problem. 00
00
Uniform re nement
levels NV NT d. of f. 1 15 16 53 2 45 64 215 3 153 256 875 4 561 1024 3539 5 2145 4096 14243
digits 0.341 0.455 0.586 0.732 0.883
q 0.74 0.99 1.10 1.13 1.14
Table 4.1
Convergence analysis and e ectivity ratio for a sequence of uniform meshes.
Adaptive re nement
levels NV NT d. of f. 1 15 16 53 3 41 64 219 5 143 250 861 8 556 1042 3616 12 2139 4124 14361
digits 0.341 0.606 0.772 1.077 1.413
q 0.74 0.50 0.73 0.93 1.11
Table 4.2
Convergence analysis and e ectivity ratio for a sequence of adapted meshes.
From the results in the column \digits" it is clear that the convergence is much better in the adaptive case than in the uniform one. we plot the estimated error vs the exact error to test the local behavior of the estimate, and its propensity to recognize the regions in the mesh needing some re nement or unre nement (Figures 4.13(a)(b) and 4.14(a)(b)). In Table 4.3 are listed the convergence rates for both uniform and adaptive cases, in the energy norm, as well as in the H01 norm (kr k) for the velocity and the L2 norm of the pressure, which both are regrouped in the energy norm. These numbers are such that 10?digits NV ? =2 h (in the least square sense). The row labeled \B-D" refers to a priori estimates results published by Brezzi and Douglas [8]. All convergence rates are based on a least square tting from the number of correct digits in the solution for each of the grids in the (uniform and adapted) sequences. In the uniform case these values are smaller than expected from the a priori estimates given by Brezzi and Douglas, except maybe for the H01 norm of the pressure, mainly because the solution does not have here the regularity required in the derivation of their estimates. Note the singularity has about the same e ect on the L2 and H01 norms for the velocity, this e ect being less obvious on the pressure. Note also that the di erence between the convergence in L2 and H01 norms is of the order of unity, 32
which corresponds to one power of h in the estimates in [8]. In the adaptive strategy the rates of convergence are increased back to the expected levels (velocity) or more (pressure), yielding a superlinear convergence in the energy norm (on a computational point of view, the adaptive re nement has somehow regularized the solution around the singularity). Finally we plot the (estimated and exact) errors in both cases. Since the estimates are given by triangles, they are transformed into errors on the vertices by averaging between all the triangles neighbor of a node, as in the previous examples. On the other hand exact errors are known at the vertices. However, in order to compare similar results, these are used to compute errors in each triangle based on the energy norm of interpolated errors at the midpoints of all interior edges. Then an averaging identical to the one above is performed to get an error at the vertices. On Figures 4.13(a)(b) (uniform case) and 4.14(a)(b) (adapted case) we can note that the estimate is in good agreement with the exact error, thus providing a nice tool for adapting grids, especially when discontinuities or steep variations in the solution occur.
(a)
(b)
Fig. 4.13. (a) Estimated error and (b) exact error on the uniform grid. The error ranges from 7 10?3 to 8 10?1 for both estimated and exact errors. :
5. Conclusion. In this paper we derived an a posteriori error estimate which can be used toward the solution of the Stokes equations on geometries of industrial interest. On the test problems of section 4 the computation of this estimate took about one fourth of the total time for the solution process. It remains to compare it with other estimates, in particular with those derived by R. Verfurth in [16].
33
(a)
(b)
Fig. 4.14. (a) Estimated error and (b) exact error on the adapted grid. The estimated error ranges from 2:10?3 to 10?1 , compared to 3:10?3 to 2:10?1 for the exact error.
REFERENCES [1] D. N. Arnold, F. Brezzi, and M. Fortin, A stable nite element for the Stokes equations, Calcolo, 21 (1984), pp. 337{344. [2] I. Babuska, O. C. Zienkiewicz, J. Gago, and E. R. de A. Oliveira, Accuracy Estimates and Adaptive Re nements in Finite Element Computations, John Wiley and Sons, 1986. [3] R. E. Bank, Analysis of a local a posteriori error estimate for elliptic equations, in Accuracy Estimates and Adaptive Re nements in Finite Element Computations, John Wiley and Sons, 1986, pp. 119{128. [4] , Pltmg users' guide - edition 5.0, tech. report, University of California, San Diego, 1988. [5] R. E. Bank and A. Weiser, Some a posteriori error estimators for partial di erential equations, Math. of Comp., 44 (April 1985), pp. 283{301. [6] R. E. Bank, B. D. Welfert, and H. Yserentant, A class of iterative methods for solving saddle point problems, Numer. Math., (to appear). [7] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, R.A.I.R.O., 8 (1974), pp. 129{151. [8] F. Brezzi and J. J. Douglas, Stabilized mixed methods for the Stokes problem, Numer. Math., 53 (1988), pp. 225{235. [9] V. Girault and P. Raviart, Finite Element Methods for Navier-Stokes Equations, SpringerVerlag, Berlin, 1986. [10] T. J. R. Hughes, L. P. Franca, and M. Balestra, A new nite element formulation for computational uid dynamics: V. circumventing the Babuska-Brezzi condition: A stable Petrov-Galerkin formulation for the Stokes problem accommodating equal-order interpolations, Comp. Meth. in Appl. Mech. and Engr., 59 (1986), pp. 85{99. [11] J. F. Maitre and F. Musy, The contraction number of a class of two-level methods ; an exact evaluation for some nite element subspaces and model problems. 1982. [12] J. F. Maitre, F. Musy, and P. Nigon, A fast solver for the Stokes equations using multigrid with a Uzawa smoother, in Notes on Numerical Fluid Mechanics, Volume 11, Vieweg, Braunschweig, 1985, pp. 77{83. [13] R. Verfurth, A combined conjugate gradient-multigrid algorithm for the numerical solution of the Stokes problem, IMA J. Numer. Anal., 4 (1984), pp. 441{455. , A multilevel algorithm for mixed problems, SIAM J. Numer. Anal., 21 (1984), pp. 264{ [14] 271. , Iterative methods for the numerical solution of mixed nite element approximations of [15] the Stokes problem, Tech. Report 379, Institut National de Recherche en Informatique et en Automatique, 1985. [16] , A posteriori error estimators for the Stokes equations, Numer. Math., 55 (1989), pp. 309{325. 34
[17] B. D. Welfert, A Posteriori Error Estimates and Adaptive Solution of Fluid Flow Problems, PhD thesis, University of California, San Diego, 1990. [18] G. Wittum, Multigrid methods for Stokes and Navier Stokes equations. Transforming Smoothers: Algorithms and Numerical Results, Tech. Report 446, Sonderforschungsbereich 123, Universitat Heidelberg, 1988.
35
type of L2 norm H01 norm L2 norm H01 norm energy norm re nement (velocity) (velocity) (pressure) (pressure) kj(; )kj uniform 1.51 0.59 0.96 0.09 0.87 adaptive 1.79 1.03 1.66 0.29 1.46 B-D 2.00 1.00 1.00 0.00 1.00 Table 4.3
Rates of convergence in L2 ; H10 and energy norms.
36
