Space and time error estimates for a first order, pressure stabilized finite element method for the incompressible Navier-Stokes equations
Abstract
In this paper we analyse a pressure stabilized, finite element method for the unsteady, incompressible Navier-Stokes equations in primitive variables; for the time discretization we focus on a fully implicit, monolithic scheme. We provide some error estimates for the fully discrete solution which show that the velocity is first order accurate in the time step and attains optimal order accuracy in the mesh size for the given spatial interpolation, both in the spaces L^2(W) and H(W) the pressure solution is shown to be order 1/2 accurate in the time step and also optimal in the mesh size. These estimates are proved assuming only a weak compatibility condition on the approximating spaces of velocity and pressure, which is satisfied by equal order interpolations.
Full text
Space and time error estimates for a rst order, pressure stabilized nite element metho d for the incompressible Navier{Stokes equations J. Blasco (1) , R. Co dina (2) (1) Departament de Matematica Aplicada I, Universitat Politecnica de Catalunya, Edici H, Campus Sud, Avgda. Diagonal 647, 08028 Barcelona. E-mail: [email protected] c.es (2) Departament de Resistencia de Materials i Estructures a l'Enginyeria, Universitat Politecnica de Catalunya, c/ Jordi Girona 1-3, Edici C1, 08034 Barcelona. E-mail: co [email protected] c.es Abstract In this pap er we analyse a pressure stabilized, nite element metho d for the unsteady, incompressible Navier{Stokes equations in primitivevariables for the time discretization we fo cus on a fully implicit, monolithic scheme. We provide some error estimates for the fully discrete solution which show that the velo city is rst order accurate in the time step and attains optimal order accuracy in the mesh size for the given spatial interp olation, b oth in the spaces L 2 () and H 1 0 () the pressure solution is shown to b e order 1 = 2 accurate in the time step and also optimal in the mesh size. These estimates are proved assuming only a weak compatibility condition on the approximating spaces of velo city and pressure, which is satised by equal order interp olations. key words: Finite elements. Incompressible ow. Pressure instability. Navier- Stokes equations. 1Intro duction The purp ose of this pap er is to provide some error estimates for a pressure stabilized, nite element metho d for the numerical solution of the unsteady, incompressible Navier{Stokes equations in the primitivevariables velo cityandpressure. The metho d was intro duced in 7] as an extention to the transient case of a technique initially develop ed for the Stokes problem (5]) and then extended to the steady, incompressible Navier{Stokes equations (6]). The stabilization of the pressure in incompressible ow problems has received much attention in the last decades. Numerical schemes have b een develop ed whichbypass the 1
need for the approximating spaces of velo city and pressure to satisfy the compatibility condition met when using standard Galerkin metho ds. Stabilized formulations were rst intro duced under the idea of Petrov{Galerkin metho ds (17]), whichthen led to Galerkin Least Squares (GLS) techniques. These were rst develop ed in the context of advection{diusion equations (18]), and then extended to the linearized, steady incompressible Navier{Stokes equation in 9] (see also 10] and the references therein). More recently, the GLS technique has evolved into the idea of subgrid{scale models (see 16], 4]). All these techniques have b een analysed in the literature for steady problems using arbitrary nite element interp olations. Error analysis b oth in space and time for stabilized formulations of transient problems have been given in 19], for advection-diusion problems, and 13], for the incompressible Navier{Stokes equations. In this last reference, the analysis was based on the assumption that the time step t is of the same order as the mesh size h : t ' h . Moreover, it was restricted to the case of piecewise linear elements. On the other hand, some combinations of nite element spaces which satisfy the discrete compatibility condition have b een analysed for the Stokes problem and proven to b e stable and yield optimal order accuracy of the solution (see, for instance, 3], 22]). Assuming a (mixed) nite element pair which satises the discrete compatibility condition, some analysis of metho ds for the unsteady problem have b een given: Heywood and Rannacher (14], 15]) proved second order error estimates in the time step and optimal order accuracy in the mesh size for a mixed metho d using a Crank-Nicholson time integration scheme Boukir et al. (1]) also proved second order estimates in time and optimal order in space for a characteristic-based metho d under a stability restriction on the time step of the form t h d= 6 , where d is the dimension of space nally, Guermond and Quartap elle (12]) have recently given an analysis of the classical fractional-step pro jection metho d of A.J.Chorin and R.Temam in its incremental form, which yields a rst order scheme also optimal in space (a second order metho d can also be develop ed). Their analysis is based on the satisfaction of the LBB condition, which has traditionally b een considered unnecessary in pro jection metho ds based on a Poisson equation for the pressure. This condition can be avoided assuming t h l +1 , where l is the order of the spatial interp olation, in the stability analysis of the non-incremental form of the metho d, but not in the convergence one. We analyse here a stabilized formulation of the unsteady problem which employs a nite element, pressure gradient pro jection technique (5]) and a fully implicit, backward Euler scheme for the time integration. Weshow that rst order accuracy in time is maintained in the fully discrete metho d, which attains optimal order accuracy in space for the given interp olation. The analysis is carried out assuming only a weak compatibility condition on the approximating spaces of velo cityand pressure, which was proven to b e satised by simplicial equal order nite elementinterp olations in 5]. The error estimates obtained are given in terms of a certain norm of the velo cityin L 2 () and H 1 0 () and the pressure and its gradient in L 2 (). We rst analyse the temp oral error by considering a semidiscrete approximation of the problem, and then study the fully discrete metho d, with b oth a linearized and a nonlinear approximation of the convective term. It has to b e remarked that the purp ose of the technique employed here is to stabilize 2
the pressure solution the instabilities due to the convective term at high cell Reynolds numbers are not addressed at with this formulation. Moreover, the interest here relies on showing how the technique that we use to stabilize the pressure, with resp ect to the spatial interp olation, can b e analysed in transient problems, regardless of the particular time integration metho d employed. We concentrate on a fully implicit, backward Euler scheme, which, although being only rst order accurate, is unconditionally stable however, other metho ds could also be considered (see 7]). The resulting scheme is computationally feasable (see 7]), and also suitable as an iterative metho d to reach steady states. Our presentation is split into two Sections. In Section 2 we state the problem to solve, recall some known prop erties of its solution and intro duce some notation we then present the semidiscrete approximation considered and nally the fully discrete, stabilized nite element metho d. In Section 3 we state and prove our error estimates, rst for the semidiscrete and then for the fully discrete problems. We rst recall a stability estimate which was proven in 7] under weak assumptions on the continuous solution then weprove some optimal order error estimates for the velo city, from which we obtain an improved stability estimate as aside pro duct. We nally analyse the pressure solution, for which we also obtain optimal order error estimates. 2 Description of the metho d 2.1 Problem statement The evolution of viscous, incompressible uid ow in a b ounded domain IR d ( d = 2 3) is governed, in the primitivevariable formulation, by the unsteady, incompressible Navier{Stokes equations: @ u @t +( u r ) u ; u + r p = f in (0 T ) (1) r u = 0 in (0 T ) (2) on (0 T ) (with T > 0agiven nal time), where u ( x t ) 2 IR d is the uid velo cityat p osition x 2 and time t 2 (0 T ), p ( x t ) 2 IR is the uid kinematic pressure, > 0 is the kinematic viscosity, f ( x t )is an external force, r is the gradient op erator, r is the divergence op erator and is the Laplacian op erator (here, and in what follows, b oldface characters denote vector quantities). Boundary conditions have to be given to complete the equation system (1){(2). For the sake of simplicity, only homogeneous Dirichlet typ e boundary conditions are considered here: u = 0 on ; (0 T ) (3) where ;= @ . An initial condition must also b e sp ecied for the velo city: u ( x 0) = u 0 ( x )in (4) The treatment of the ab ove equations of motion requieres of the usual Sob olev spaces H m () m 0, consisting of functions with distributional derivatives up to 3
order m b elonging to L 2 (). The scalar pro duct in H m () is denoted by ( u v ) m (the subscript m may be omitted when it equals 0) and its norm by jj u jj m . The closed subspaces H 1 0 (), consisting of functions in H 1 () with zero trace on ;, and L 2 0 (), made up with functions in L 2 () with zero mean on , will also be needed. Also, let H ; 1 () denote the dual space of H 1 0 (), the duality between these two spaces b eing denoted by h i , and let: W = f u 2 H 1 0 () = r u =0 g Assuming u 0 2 H 1 0 () and f 2 L 2 (0 T H ; 1 ()) and if is b ounded and Lipschitz continuous, problem (1){(2){(3){(4) has at least one solution u 2 L 1 (0 T L 2 ()) \ L 2 (0 T H 1 0 ()) (see 24]). Uniqueness and more regularity of the solution can be achieved by assuming more regularity on f , u 0 and . In particular, we assume hereafter that the continuous solution ( u p ) of (1){(2){(3){(4) is unique and satises: R1) u 2 L 1 (0 T H 2 ()) \C 0 (0 T W ) p 2 L 1 (0 T H 1 ()) \C 0 (0 T L 2 0 ()) R2) u t 2 L 2 (0 T L 2 ()) R3) R T 0 t jj u tt ( t ) jj 2 ; 1 dt C R4) R T 0 jj u tt ( t ) jj 2 W 0 dt C Here, and in what follows, the subscript t is employed for @ @t ,andweuse C as a generic constant dep ending of f , u 0 , and ,but not on the time step t nor on the mesh size h also, W 0 is the dual space of W . Sucient conditions for R1 , R2 and R3 to hold can be found 14] for R4 ,see 20] and 21]. In particular, it is required that f 2 L 2 (0 T L 2 ()), which we assume from now on. Let us call V = H 1 0 () and Q = L 2 0 (). In what follows the following notation will be used for the weak form of the dierent terms in equations (1){(2): a ( u v ) = ( r u r v ) u v 2 V b ( q v ) = ; ( q r v ) v 2 V q 2 Q c ( u v w ) = (( u r ) v w ) + 1 2 ( r u ) v w u v w 2 V All these forms are continuous on the sp ecied spaces, and the expression taken for the trilinear form c arising from the convective term in (1) is skew{symmetric in its last two arguments (see 24]) under the incompressibility condition (2), this expression is equivalent to that obtained from the original convectiveterm in (1). Besides, a is co ercive as a consequence of the Poincare{Friedrics inequality, that is, there exists a constant K a > 0such that: a ( u u ) = jjr u jj 2 0 K a jj u jj 2 1 8 u 2 V: and b satises the (continuous) inf{sup condition, that is, there exists a constant K b > 0 such that: inf q 2 Q sup v 2 V b ( q v ) jj v jj 1 jj q jj 0 K b > 0 : (5) 4
(inma and suprema are always taken with resp ect to non zero functions). Condition (5) is usually refered to as the inf{sup or LBB condition, after the work of O.A. Ladyzhenskaya, I. Babu^ska and F. Brezzi. Finally, c satises other continuity prop erties, some of which are (see, for instance, 8]): c ( u v w ) 8 > > > < > > > : C jj u jj 1 jj v jj 1 jj w jj 1 C jj u jj 0 jj v jj 2 jj w jj 1 C jj u jj 2 jj v jj 1 jj w jj 0 C jj u jj 0 jj v jj 1 jj w jj L 1 () 2.2 Finite element approximation The numerical approximation of problem (1){(2) that we analyse here was intro duced in 7] as an extention to the transient case of a nite element metho d originally develop ed for steady problems. It is well known that discrete approximations of incompressible ow problems in primitive variables are restricted by the discrete inf{sup condition, that is, the discrete counterpart of condition (5) this prevents the use of many simple nite element combinations for the discrete spaces of the velo city and the pressure, such as equal order ones. The metho ds based on a pressure gradient pro jection circumvent this restriction by intro ducing the pro jection of the gradient of the discrete pressure onto the space of discrete velo cities as a new variable of the problem this allows, in particular, the use of equal order interp olations. In the transient case, this metho dology can b e applied together with dierent time integration schemes we concentrate here on an implicit, monolithic scheme using the trap ezoidal rule, but extentions to other schemes such as fractional{step or multistep metho ds can be derived in a similar way (see 7] for a description of some of them). 2.2.1 Semidiscrete problem We consider a parameter 2 (0 1] and discretize equations (1){(2) in time rst by the following implicit scheme, which we write in variational form: given a time step size t > 0, let N = T=t ] ; 1 for n 2 f 0 :::N g ,let t n = nt given u n 2 V and p n 2 Q , approximations of u ( t n )and p ( t n ), resp ectively, nd u n +1 2 V and p n +1 2 Q such that: ( u n +1 ; u n t v ) + c ( u n + u n + v ) + a ( u n + v ) + ( r p n + v ) = ( f n + v ) (6) b ( q u n +1 ) = 0 (7) for all ( v q ) 2 V Q , where for a given function g the notation g n + stands for: g n + = g n +1 + (1 ; ) g n The parameter which app ears in the approximation of the nonlinear term in (6) can take the values 0 and 1, corresp onding to a linearized and a nonlinear approximation of convection, resp ectively. The rst option is suitable in the rst order, backward 5
Euler case = 1, since the approximation it provides is also rst order accurate and it results in a lower computational cost of the fully discrete problem (which is then linear in each time step) the second option, however, enhances stability for highly convective ows and is compulsory in the Crank{Nicholson case = 1 = 2 to maintain second order accuracy. In this sense, we use the expression c ( u n + u n + v )which difers from ( c ( u u v )) n + (the form to which strict application of the trap ezoidal rule would lead) by a second order term and is computationally simpler. Moreover, we assume that the semidiscrete pressures satisfy r p n +1 2 L 2 () conditions on the data for this condition to hold can be found, for instance, in 11]. 2.2.2 Fully discrete metho d Wenow pro ceed to intro duce a spatial approximation of the semidiscrete problem (6){ (7). Let h denote a nite element partition of the domain of diameter h . We assume that all the element domains K 2 h are the image of a reference element ^ K through p olynomial mappings F K , ane for simplicial elements, bilinear for quadrilaterals and trilinear for hexahedra. On ^ K we dene the p olynomial spaces R k ( ^ K ) where, as usual, R k = P k for simplicial elements and R k = Q k for quadrilaterals and hexahedra. The nite element spaces we need are: Q h = n q h 2 C 0 () \ L 2 0 () j q h j K = ^ q F ; 1 K ^ q 2 R k q ( ^ K ) K 2 h o V h = n v h 2 ( C 0 ()) d j v h j K =^ v F ; 1 K ^ v 2 ( R k v ( ^ K )) d K 2 h o V h 0 = n v h 2 V h j v h j ; = 0 o Notice that b oth the velo city and pressure nite element spaces V h 0 and Q h are referred to the same partition and b oth are made up with continuous functions. These nite element spaces satisfy the following approximating prop erties (see 23], for instance): given v 2 H r () r 2 and q 2 H s () s 1, there exist ! h 1 ( v ) 2 V h 0 , ! h 2 ( q ) 2 Q h and ! h 3 ( r q ) 2 V h such that: k v ; ! h 1 ( v ) k m 1 C 1 h k 1 ; m 1 k v k k 1 k q ; ! h 2 ( q ) k m 2 C 2 h k 2 ; m 2 k q k k 2 kr q ; ! h 3 ( r q ) k m 3 C 3 h k 3 ; m 3 kr q k k 3 for 0 m i k i ( i =1 2 3), where: k 1 =min f rk v +1 g k 2 =min f s k q +1 g k 3 =min f s ; 1 k v +1 g Let now > 0 b e a given parameter. Given ( u n h p n h n h ) 2 V h 0 Q h V h , approximations of ( u n p n r p n ), we discretize (6){(7) in space by nding ( u n +1 h p n +1 h n +1 h ) 2 V h 0 Q h V h such that: 6
( u n +1 h ; u n h t v h )+ c ( u n + h u n + h v h )+ a ( u n + h v h )+( r p n + h v h ) = ( f n + v h ) (8) ; b ( q h u n +1 h )+ ( r p n +1 h r q h ) ; ( n + h r q h ) = 0 (9) ; ( r p n +1 h h )+( n +1 h h ) = 0 (10) for all ( v h q h h ) 2 V h 0 Q h V h , where again either =0 or =1. Equation (10) says that n +1 h is the L 2 -pro jection of r p n +1 h onto the space V h thus, the cases =0 and = 1 corresp ond to an explicit and an implicit approximation of the pressure gradient pro jection in the mo died continuity equation (9), resp ectively (see 7]). In the formulation (8){(9){(10) wehave used a 'global' parameter , with the same value on all the element domains the numerical analysis of this metho d then requires of some regularity prop erties of the nite element mesh suchas its quasi-uniformity . However, this restriction can be relaxed by considering a set of elemental parameters K , K 2 h and replacing the L 2 -scalar pro ducts app earing in (9){(10) by a sum of pro ducts weighed in each element by K .This extention to lo cal parameters was analysed in 6] for the steady, incompressible Navier{Stokes equations, and the analysis given there can b e readily applied to the transient case. We restrict our attention here to the global parameter case to simplify the presentation. 3 Stability and error analysis We now present a numerical analysis of the nite element metho d (8){(9){(10). For the time approximation, we restrict to the fully implicit, backward Euler case = 1, which is rst order accurate in the time step. We split the errors of the metho d into atemp oral error, due to the semidiscretization (6){(7), and a spatial error, due to the stabilized, fully discrete metho d (8){(9){(10). In the case of study = 1, rst order accuracy in the time step for the semidiscrete velo city solution can b e shown by standard arguments we include a pro of of this result for completeness. We consider both the linearized metho d = 0 and the fully nonlinear scheme = 1. We then concentrate on the spatial approximation in the implicit pressure gradientcase =1. 3.1 Error estimates for the semidiscrete solution Let us dene the continuous errors (as for the spatial variables) as: e n +1 c = u ( t n +1 ) ; u n +1 r n +1 c = p ( t n +1 ) ; p n +1 g n +1 c = r r n +1 c We then have: Theorem 1: assume R1 , R2 and R4 hold. Then, thereisaconstant C independent of t such that: 7
jj e N +1 c jj 2 0 + t N X n =0 jj e n +1 c jj 2 1 Ct 2 (11) If =1 , (11) holds for suciently smal l t . proof: We call R n the truncation error dened by: 1 t ( u ( t n +1 ) ; u ( t n )) ; u ( t n +1 ) + ( u ( t n +1 ) r ) u ( t n +1 ) + r p ( t n +1 ) = f ( t n +1 ) + R n (12) so that: R n = 1 t Z t n +1 t n ( t ; t n ) u tt ( t ) dt Multiplying (12) by v 2 V and (2) (at t = t n +1 ) by q 2 Q , and subtracting (6) (with =1) and (7) from them, resp ectively, we nd: ( e n +1 c ; e n c t v )+ ( r e n +1 c r v )+( r r n +1 c v ) = h R n v i + c ( u n + u n +1 v ) ; c ( u ( t n +1 ) u ( t n +1 ) v ) (13) b ( q e n +1 c ) = 0 (14) Taking v =2 t e n +1 c in (13) and q = r n +1 c in (14), and using the identity( a ; b 2 a )= j a j 2 ;j b j 2 + j a ; b j 2 ,we get: jj e n +1 c jj 2 0 ;jj e n c jj 2 0 + jj e n +1 c ; e n c jj 2 0 + 2 t jjr e n +1 c jj 2 0 = 2 t h R n e n +1 c i + 2 t NLT where NLT stands for: NLT = c ( u n + u n +1 e n +1 c ) ; c ( u ( t n +1 ) u ( t n +1 ) e n +1 c ) For the Taylor residual term, one has (see, for instance, 20]): 2 t h R n e n +1 c i t 3 jjr e n +1 c jj 2 0 + Ct 2 Z t n +1 t n jj u tt jj 2 W 0 dt The treatment of the NLT is dierentinthe cases = 0 and =1: Linearized case: when =0, wehave: 2 t NLT=2 t c ( u n u n +1 e n +1 c ) ; c ( u ( t n +1 ) u ( t n +1 ) e n +1 c ) = 2 t ; c ( u n e n +1 c e n +1 c ) ; c ( e n c u ( t n +1 ) e n +1 c ) ; c ( u ( t n +1 ) ; u ( t n ) u ( t n +1 ) e n +1 c ) = T 1 + T 2 + T 3 8
where T 1 = 0 due to the skew{symmetry of the trilinear form c , and, due to its continuity prop erties and the regularity prop erty R1 of u : T 2 = ; 2 t c ( e n c u ( t n +1 ) e n +1 c ) Ct jj e n c jj 0 jj u ( t n +1 ) jj 2 jj e n +1 c jj 1 t 3 jjr e n +1 c jj 2 0 + Ct jj e n c jj 2 0 T 3 = ; 2 t c ( u ( t n +1 ) ; u ( t n ) u ( t n +1 ) e n +1 c ) Ct jj u ( t n +1 ) ; u ( t n ) jj 0 jj u ( t n +1 ) jj 2 jj e n +1 c jj 1 t 3 jjr e n +1 c jj 2 0 + Ct 2 Z t n +1 t n jj u t jj 2 0 dt Therefore: jj e n +1 c jj 2 0 ; jj e n c jj 2 0 + jj e n +1 c ; e n c jj 2 0 + t jj e n +1 c jj 2 1 (15) Ct 2 Z t n +1 t n jj u tt jj 2 W 0 dt + Ct 2 Z t n +1 t n jj u t jj 2 0 dt + Ct jj e n c jj 2 0 Adding up (15) for n = 0 :::N , and using the regularity prop erties R2 and R4 of the continuous solution, we get: jj e N +1 c jj 2 0 + N X n =0 jj e n +1 c ; e n c jj 2 0 + t N X n =0 jj e n +1 c jj 2 1 Ct 2 + Ct N X n =0 jj e n c jj 2 0 Applying the discrete Gronwall inequality,thisimplies: jj e N +1 c jj 2 0 + N X n =0 jj e n +1 c ; e n c jj 2 0 + t N X n =0 jj e n +1 c jj 2 1 Ct 2 (16) and (11) follows. Nonlinear case: when =1 we have: 2 t NLT = 2 t c ( u n +1 u n +1 e n +1 c ) ; c ( u ( t n +1 ) u ( t n +1 ) e n +1 c ) = 2 t ; c ( u n +1 e n +1 c e n +1 c ) ; c ( e n +1 c u ( t n +1 ) e n +1 c ) = T 1 + T 2 where again T 1 =0 due to the skew{symmetry of the trilinear form c ,and: T 2 = ; 2 t c ( e n +1 c u ( t n +1 ) e n +1 c ) Ct jj e n +1 c jj 0 jj u ( t n +1 ) jj 2 jj e n +1 c jj 1 t 3 jjr e n +1 c jj 2 0 + Ct jj e n +1 c jj 2 0 Therefore: jj e n +1 c jj 2 0 ; jj e n c jj 2 0 + jj e n +1 c ; e n c jj 2 0 + t jj e n +1 c jj 2 1 (17) Ct 2 Z t n +1 t n jj u tt jj 2 W 0 dt + Ct 2 Z t n +1 t n jj u t jj 2 0 dt + Ct jj e n +1 c jj 2 0 9
C h jj e n + d jj 1 jj e n +1 d jj 1 jj P h 1 ( r q h ;r p n +1 h ) jj 0 C h ( jj e n + d jj 2 1 + jj e n +1 d jj 2 1 ) jj P h 1 ( r q h ;r p n +1 h ) jj 0 c ( e n + c e n +1 d P h 1 ( r q h ;r p n +1 h )) C jj e n + c jj 1 jj e n +1 d jj 1 jj P h 1 ( r q h ;r p n +1 h ) jj 1 C h jj e n + c jj 1 jj e n +1 d jj 1 jj P h 1 ( r q h ;r p n +1 h ) jj 0 C h ( jj e n + c jj 2 1 + jj e n +1 d jj 2 1 ) jj P h 1 ( r q h ;r p n +1 h ) jj 0 ; c ( u ( t n + ) e n +1 d P h 1 ( r q h ;r p n +1 h )) C jj u ( t n + ) jj 2 jj e n +1 d jj 1 jj P h 1 ( r q h ;r p n +1 h ) jj 0 C jj e n +1 d jj 1 jj P h 1 ( r q h ;r p n +1 h ) jj 0 due to Remark 1 and the regularity of the continuous velo city. Assuming h C 1 = 2 in the last term, we get: T 2 C I 1 ( p n +1 q h ) + 1 t jj e n +1 d ; e n d jj 0 + 1 = 2 h ( jj e n +1 d jj 1 + jj e n d jj 1 ) + h ( jj e n +1 d jj 2 1 + jj e n d jj 2 1 + jj e n +1 c jj 2 1 + jj e n c jj 1 ) Moreover, due to condition (23) and since P h 3 = Id ; P h 12 and n +1 h = P h 12 ( r p n +1 h ), we have: T 3 = jj P h 2 ( r q h ) ; P h 2 ( r p n +1 h ) jj 0 C jj P h 1 ( r q h ) ; P h 1 ( r p n +1 h ) jj 0 + jj P h 3 ( r q h ) ; P h 3 ( r p n +1 h ) jj 0 C T 2 + jj P h 3 ( r q h ) jj 0 + jj P h 3 ( r p n +1 h ) jj 0 C T 2 + jjr q h ;r p n +1 jj 0 + jjr p n +1 ; P h 12 ( r q h ) jj 0 + G n +1 = C ( T 2 + I 1 ( p n +1 q h ) + T 1 + G n +1 ) Finally: T 4 = jj P h 3 ( r p n +1 h ) jj 0 = G n +1 and (30) follows. In our convergence analysis we will also need the following assumption: H4: There exists C > 0 indep endent of h and t such that: t Ch 2 (31) This condition do es not imp ose an upp er bound on the time step, so that the metho d remains unconditionally stable (see also Remark 5). Our main result of this Section is the following: Theorem 3: assume R1, R2, R4, H1, H2, H3 and H4 hold then, there exists a constant C > 0 independent of t and h such that, for smal l enough h and, if = 1 , 16
smal l enough t : jj e N +1 d jj 2 0 + t N X n =0 jj e n +1 d jj 2 1 C ( E ( h )) 2 + E ( h ) t 2 (32) proof: Let us call: A = e n +1 d ; e n d t e n +1 d + ( r e n +1 d r e n +1 d ) + ( r r n +1 d e n +1 d ) + ( r e n +1 d r n +1 d ) + ( r r n +1 d r r n +1 d ) ; ( g n +1 d r r n +1 d ) ; ( g n +1 d r r n +1 d ) + ( g n +1 d g n +1 d ) = 1 2 t jj e n +1 d jj 2 0 ;jj e n d jj 2 0 + jj e n +1 d ; e n d jj 2 0 + jj e n +1 d jj 2 1 + jj n +1 h ;r p n +1 h jj 2 0 Given ( v h q h h ) 2 V h 0 Q h V h arbitrary,wetake( v h ; u n +1 h q h ; p n +1 h h ; n +1 h ) as test functions in (26)-(27)-(28), resp ectively, to get: A = e n +1 d ; e n d t u n +1 ; v h + ( r e n +1 d r ( u n +1 ; v h )) + ( r r n +1 d u n +1 ; v h ) + c ( u n + h u n +1 h u n +1 h ; v h ) ; c ( u n + u n +1 u n +1 h ; v h ) + ( r e n +1 d p n +1 ; q h ) + ( r r n +1 d ; g n +1 d h ;r q h ) We bound each term as follows: e n +1 d ; e n d t u n +1 ; v h 1 t jj e n +1 d ; e n d jj 0 I 0 ( u n +1 v h ) 1 4 t jj e n +1 d ; e n d jj 2 0 + C t I 0 ( u n +1 v h ) 2 ( r e n +1 d r ( u n +1 ; v h )) C jj e n +1 d jj 1 I 1 ( u n +1 v h ) 10 jj e n +1 d jj 2 1 + CI 1 ( u n +1 v h ) 2 ( r e n +1 d p n +1 ; q h ) 10 jj e n +1 d jj 2 1 + CI 0 ( p n +1 q h ) 2 ( r r n +1 d ; g n +1 d h ;r q h ) = ( r p n +1 h ; n +1 h h ;r q h ) G n +1 ( I 0 ( r p n +1 h )+ I 1 ( p n +1 q h )) Ch 2 G n +1 ( I 0 ( r p n +1 h )+ I 1 ( p n +1 q h )) ; h 2 3 G 2 n +1 + Ch 2 I 0 ( r p n +1 h ) 2 + Ch 2 I 1 ( p n +1 q h ) 2 where ; was dened in (20). Moreover, due to Lemma 1 we have: ( r r n +1 d u n +1 ; v h ) jjr r n +1 d jj 0 I 0 ( u n +1 v h ) C I 0 ( r p n +1 h ) + I 1 ( p n +1 q h ) + 1 t jj e n +1 d ; e n d jj 0 + G n +1 17
+ 1 = 2 h ( jj e n +1 d jj 1 + jj e n d jj 1 ) + h ( jj e n +1 d jj 2 1 + jj e n d jj 2 1 + jj e n +1 c jj 2 1 + jj e n c jj 1 ) I 0 ( u n +1 v h ) C h 2 I 0 ( r p n +1 h ) 2 + h 2 I 1 ( p n +1 q h ) 2 + 1 h 2 I 0 ( u n +1 v h ) 2 + 1 4 t jj e n +1 d ; e n d jj 2 0 + C t I 0 ( u n +1 v h ) 2 + ; h 2 3 G 2 n +1 + 10 jj e n +1 d jj 2 1 + C 1 = 2 h jj e n d jj 1 I 0 ( u n +1 v h ) + ( jj e n +1 d jj 2 1 + jj e n d jj 2 1 + jj e n +1 c jj 2 1 + jj e n c jj 2 1 ) C h I 0 ( u n +1 v h ) We split the convective terms the following way: c ( u n + h u n +1 h u n +1 h ; v h ) ; c ( u n + u n +1 u n +1 h ; v h ) = c ( e n + d u n +1 u n +1 h ; v h ) + c ( u n + h e n +1 d u n +1 h ; v h ) = ; c ( e n + d e n +1 c u n +1 h ; v h ) + c ( e n + d u ( t n +1 ) u n +1 h ; v h ) ; c ( e n + d e n +1 d u n +1 h ; v h ) ; c ( e n + c e n +1 d u n +1 h ; v h ) + c ( u ( t n + ) e n +1 d u n +1 h ; v h ) = c ( e n + d e n +1 c e n +1 d ) ; c ( e n + d e n +1 c u n +1 ; v h ) ; c ( e n + d u ( t n +1 ) e n +1 d ) + c ( e n + d u ( t n +1 ) u n +1 ; v h ) + c ( e n + d e n +1 d e n +1 d ) ; c ( e n + d e n +1 d u n +1 ; v h ) + c ( e n + c e n +1 d e n +1 d ) ; c ( e n + c e n +1 d u n +1 ; v h ) ; c ( u ( t n + ) e n +1 d e n +1 d ) + c ( u ( t n + ) e n +1 d u n +1 ; v h ) Due to the continuity prop erties of the trilinear form c ,its skew symmetry in its last two arguments, the results of Theorem 1, the regularity assumed for the continuous solution u and Young's inequality,we have: c ( e n + d e n +1 c e n +1 d ) C jj e n + d jj 1 jj e n +1 c jj 1 jj e n +1 d jj 1 Ct 1 = 2 jj e n + d jj 1 jj e n +1 d jj 1 10 jj e n +1 d jj 2 1 + Ct jj e n + d jj 2 1 ; c ( e n + d e n +1 c u n +1 ; v h ) C jj e n + d jj 1 jj e n +1 c jj 1 I 1 ( u n +1 v h ) Ct 1 = 2 jj e n + d jj 1 I 1 ( u n +1 v h ) t jj e n + d jj 2 1 + CI 1 ( u n +1 v h ) 2 ; c ( e n + d u ( t n +1 ) e n +1 d ) C jj e n + d jj 0 jj u ( t n +1 ) jj 2 jj e n +1 d jj 1 C jj e n + d jj 2 0 + 10 jj e n +1 d jj 2 1 c ( e n + d u ( t n +1 ) u n +1 ; v h ) C jj e n + d jj 0 jj u ( t n +1 ) jj 2 I 1 ( u n +1 v h ) C jj e n + d jj 2 0 + I 1 ( u n +1 v h ) 2 c ( e n + d e n +1 d e n +1 d ) = 0 18
; c ( e n + d e n +1 d u n +1 ; v h ) C jj e n + d jj 1 jj e n +1 d jj 1 I 1 ( u n +1 v h ) 10 jj e n +1 d jj 2 1 + C jj e n + d jj 2 1 I 1 ( u n +1 v h ) 2 c ( e n + c e n +1 d e n +1 d ) = 0 ; c ( e n + c e n +1 d u n +1 ; v h ) C jj e n + c jj 1 jj e n +1 d jj 1 I 1 ( u n +1 v h ) 10 jj e n +1 d jj 2 1 + C jj e n + c jj 2 1 I 1 ( u n +1 v h ) 2 ; c ( u ( t n + ) e n +1 d e n +1 d ) = 0 c ( u ( t n + ) e n +1 d u n +1 ; v h ) C jj u ( t n + ) jj 2 jj e n +1 d jj 1 I 0 ( u n +1 v h ) 10 jj e n +1 d jj 2 1 + CI 0 ( u n +1 v h ) 2 Taking all the previous inequalities into account, and using (20), we nd: jj e n +1 d jj 2 0 ; jj e n d jj 2 0 + jj e n +1 d ; e n d jj 2 0 + t jj e n +1 d jj 2 1 + t h 2 G 2 n +1 Ct I 1 ( u n +1 v h ) 2 + h 2 I 1 ( p n +1 q h ) 2 + h 2 I 0 ( r p n +1 h ) 2 + 1 h 2 I 0 ( u n +1 v h ) 2 + I 0 ( p n +1 q h ) 2 + CI 0 ( u n +1 v h ) 2 + t ( jj e n +1 d jj 2 1 + jj e n d jj 2 1 ) C h I 0 ( u n +1 v h ) + t ( jj e n +1 c jj 2 1 + jj e n c jj 2 1 ) C h I 0 ( u n +1 v h ) + Ct 2 jj e n + d jj 2 1 + Ct jj e n + d jj 2 0 + C t 1 = 2 h jj e n d jj 1 I 0 ( u n +1 v h ) + Ct jj e n + d jj 2 1 I 1 ( u n +1 v h ) 2 + Ct jj e n + c jj 2 1 I 1 ( u n +1 v h ) 2 Taking the inmum with resp ect to ( v h q h h ) 2 V h 0 Q h V h , we get: jj e n +1 d jj 2 0 ; jj e n d jj 2 0 + jj e n +1 d ; e n d jj 2 0 + t jj e n +1 d jj 2 1 + t h 2 G 2 n +1 Ct ( E n ( h )) 2 + Ch 2 ( E n ( h )) 2 + Ct ( jj e n +1 d jj 2 1 + jj e n d jj 2 1 + jj e n +1 c jj 2 1 + jj e n c jj 2 1 ) E n ( h ) + Ct 2 jj e n + d jj 2 1 + Ct jj e n + d jj 2 0 + Ct 1 = 2 jj e n d jj 1 E n ( h ) + Ct jj e n + d jj 2 1 + jj e n + c jj 2 1 ( E n ( h )) 2 (33) Adding up (33) from n =0to N , using assumption H4 , the denition of E ( h ) and the estimates of Theorem 1, we get: jj e N +1 d jj 2 0 + N X n =0 jj e n +1 d ; e n d jj 2 0 + t N X n =0 jj e n +1 d jj 2 1 + t h 2 N X n =0 G 2 n +1 C ( E ( h )) 2 + C ( t N X n =0 jj e n +1 d jj 2 1 ) E ( h ) + C ( t N X n =0 jj e n +1 c jj 2 1 ) E ( h ) + Ct 2 N X n =0 jj e n + d jj 2 1 + Ct N X n =0 jj e n + d jj 2 0 + C ( t 1 = 2 N X n =0 jj e n d jj 1 ) E ( h ) 19
+ C t N X n =0 jj e n + d jj 2 1 + t N X n =0 jj e n + c jj 2 1 ( E ( h )) 2 C ( E ( h )) 2 + C ( t N X n =0 jj e n +1 d jj 2 1 ) E ( h ) + Ct 2 E ( h ) + Ct 2 N X n =0 jj e n + d jj 2 1 + Ct N X n =0 jj e n + d jj 2 0 + C ( t N X n =0 jj e n d jj 2 1 ) 1 = 2 E ( h ) C ( E ( h )) 2 + C ( t N X n =0 jj e n +1 d jj 2 1 ) E ( h ) + Ct 2 E ( h ) + Ct 2 N X n =0 jj e n + d jj 2 1 + Ct N X n =0 jj e n + d jj 2 0 + (1 = 2) ( t N X n =0 jj e n +1 d jj 2 1 ) since E ( h ) Ch , ( E ( h )) 2 E ( h )for h small enough and j Z j l 1 ( X ) C j Z j l 2 ( X ) for any Z and X (see Remark 1). For suciently small h ,the second term in the right hand side can b e passed over to the left hand side, since E ( h ) tends to 0 as h tends to 0. By the discrete Gronwall inequality, this implies, for suciently small t in the case =1: jj e N +1 d jj 2 0 + N X n =0 jj e n +1 d ; e n d jj 2 0 + t N X n =0 jj e n +1 d jj 2 1 + t h 2 N X n =0 G 2 n +1 C ( E ( h )) 2 + Ct 2 E ( h ) (34) and (32) follows. Remark 3: for equal order interp olations of degree k , the spatial error function E ( h ) b ehaves like h k , the worst case b eing that of linear ( P 1 ) and multilinear ( Q 1 ) elements. In general, one always has E ( h ) Ch due to assumption (31), this result proves in particular that the discrete velo cities are b ounded in l 1 ( H 1 0 ()) by a constant indep endent of t and h , since: jj u n +1 h jj 1 jj u ( t n +1 ) jj 1 + jj e n +1 c jj 1 + jj e n +1 d jj 1 C 1+ t 1 = 2 +( ( E ( h )) 2 t ) 1 = 2 C 1+( h 2 t ) 1 = 2 C This is the key p ointtoobtainimproved stability estimates in the next Section. Remark 4: the last term in the estimate (32) for the discrete velo city is due to the presence of the convective term in the equations (it is not present in an analysis of the linear Stokes case) and arises from the estimates of the semidiscrete problem. Again, since E ( h ) Ch , this extra term is always smaller than t 2 ,and the metho d remains rst order accurate in time for the velo city. 3.4 Improved stability estimate As a consequence of the convergence analysis of the previous Section, the stability results of Section 3.2 can be improved as follows: 20
Prop osition 2: assume R1, R2, R4, H1, H2, H3 and H4 hold then, there exists a constant C > 0 independent of t and h such that, for smal l enough h and, if =1 , smal l enough t : t h 2 N X n =0 jjr p n +1 h jj 2 0 C (35) proof: in a similar way to 7], taking v = u n +1 h in (8) (with = 1), q h = p n +1 h in (9) and h = n +1 h in (10), and adding them up, we get: ( u n +1 h ; u n h t u n +1 h ) + jjr u n +1 h jj 2 0 + jjr p n +1 h ; n +1 h jj 2 0 = ( f n +1 u n +1 h ) (36) From (36), it is found that: jj u N +1 h jj 2 0 + N X n =0 jj u n +1 h ; u n h jj 2 0 + t N X n =0 jj u n +1 h jj 2 1 + t N X n =0 jjr p n +1 h ; n +1 h jj 2 0 C ( t N X n =0 jj f n +1 jj 2 0 + 1) C ( Z T 0 jj f ( t ) jj 2 0 dt +1) = C Thus, the third comp onent P h 3 ( r p n +1 h )= r p n +1 h ; n +1 h in the decomp osition of r p n +1 h in E h is b ounded due to assumption (23), it only remains to b ound P h 1 ( r p n +1 h ), which b elongs to V h 0 . Using the continuity of the forms a and c ,the inverse estimate (22) and the result of Remark 2, we have: jj P h 1 ( r p n +1 h ) jj 2 0 = ( r p n +1 h P h 1 ( r p n +1 h )) = ; ( u n +1 h ; u n h t P h 1 ( r p n +1 h )) ; a ( u n +1 h P h 1 ( r p n +1 h )) ; c ( u n + h u n +1 h P h 1 ( r p n +1 h )) + ( f n +1 P h 1 ( r p n +1 h )) jj P h 1 ( r p n +1 h ) jj 0 1 t jj u n +1 h ; u n h jj 0 + jj f n +1 jj 0 + C h jj u n +1 h jj 1 + C h jj u n + h jj 1 jj u n +1 h jj 1 jj P h 1 ( r p n +1 h ) jj 0 1 t jj u n +1 h ; u n h jj 0 + jj f n +1 jj 0 + C h Dividing this estimate by jj P h 1 ( r p n +1 h ) jj 0 , squaring the result, multiplying by t h 2 and adding up for n =0 :::N , we nd: t h 2 N X n =0 jjr p n +1 h jj 2 0 C h 2 t N X n =0 jj u n +1 h ; u n h jj 0 + 1) C due to the assumed b ehaviour (31) on the time step size. 21
3.5 Error estimates for the pressure We b egin this Section with an estimate for the discrete pressure gradient: Prop osition 3: assume R1, R2, R3, R4, H1, H2, H3 and H4 hold then, there exists a constant C > 0 independent of t and h such that, for smal l enough h and, if =1 , smal l enough t : t h 2 N X n =0 jjr r n +1 d jj 2 0 C ( E ( h )) 2 + E ( h ) t 2 (37) proof: from Lemma 1, we have: jjr r n +1 d jj 2 0 C ( ( I 0 ( r p n +1 h )) 2 + ( I 1 ( p n +1 q h )) 2 + 1 t 2 jj e n +1 d ; e n d jj 2 0 + h 2 jj e n +1 d jj 2 1 + ( G n +1 ) 2 ) Thus: t h 2 N X n =0 jjr r n +1 d jj 2 0 C ( t h 2 N X n =0 ( I 0 ( r p n +1 h )) 2 + t h 2 N X n =0 ( I 1 ( p n +1 q h )) 2 + h 2 t jj e n +1 d ; e n d jj 2 0 + t N X n =0 jj e n +1 d jj 2 1 + t h 2 N X n =0 ( G n +1 ) 2 ) Taking the inmum with resp ect to h and q h and using (31), this implies: t h 2 N X n =0 jjr r n +1 d jj 2 0 C ( t N X n =0 ( E n ( h )) 2 + N X n =0 jj e n +1 d ; e n d jj 2 0 + t N X n =0 jj e n +1 d jj 2 1 + t h 2 N X n =0 ( G n +1 ) 2 ) and (37) follows from (34) and the denition of E ( h )(29) . Since we have obtained error estimates for the fully discrete pressure gradient and the semidiscrete pressure itself, we now present some estimates for the fully discrete pressure solution, which are based on a classical duality argument: Prop osition 4: assume R1, R2, R3, R4, H1, H2, H3 and H4 hold then, there exists a constant C > 0 independent of t and h such that, for smal l enough h and, if =1 , smal l enough t : t 2 N X n =0 jj r n +1 d jj 2 0 C ( E ( h )) 2 + t 2 (38) 22
proof: Let z 2 H 1 0 () and 2 L 2 0 () b e the solution of the following Stokes problem: ; z + r = 0 in (39) r z = r n +1 d in z = 0 on ; Standard results for this problem yield: jj z jj 1 C jj r n +1 d jj 0 (40) jj jj 0 C jj r n +1 d jj 0 If z h 2 V h 0 now satises: jj z ; z h jj m Ch 1 ; m jj z jj 1 (41) for m =0 1, we have: jj r n +1 d jj 2 0 = ( r n +1 d r n +1 d ) = ( r z r n +1 d ) = ; ( z r r n +1 d ) = ; ( z ; z h r r n +1 d ) ; ( z h r r n +1 d ) = ; ( z ; z h r r n +1 d ) + ( e n +1 d ; e n d t z h ) + ( r e n +1 d r z h ) ; c ( u n + h u n +1 h z h ) + c ( u n + u n +1 z h ) Furthermore: ; ( z ; z h r r n +1 d ) jj z ; z h jj 0 jjr r n +1 d jj 0 Ch jjr r n +1 d jj 0 jj z jj 1 Ch jjr r n +1 d jj 0 jj r n +1 d jj 0 ( e n +1 d ; e n d t z h ) 1 t jj e n +1 d ; e n d jj 0 jj z h jj 0 1 t jj e n +1 d ; e n d jj 0 ( jj z ; z h jj 0 + jj z jj 0 ) 1 t jj e n +1 d ; e n d jj 0 ( Ch jj z jj 1 + C jj z jj 1 ) 1 t jj e n +1 d ; e n d jj 0 jj r n +1 d jj 0 ( r e n +1 d r z h ) C 1 = 2 jj e n +1 d jj 1 jj z h jj 1 C 1 = 2 jj e n +1 d jj 1 jj r n +1 d jj 0 ; c ( u n + h u n +1 h z h ) + c ( u n + u n +1 z h ) = c ( u n + e n +1 d z h ) + c ( e n + d u n +1 h z h ) C jj u n + jj 1 jj e n +1 d jj 1 + jj e n + d jj 1 jj u n +1 h jj 1 jj z h jj 1 C jj e n +1 d jj 1 + jj e n + d jj 1 ( jj z ; z h jj 1 + jj z jj 1 ) C 1 = 2 jj e n +1 d jj 1 + jj e n + d jj 1 jj r n +1 d jj 0 Estimate (38) is obtained dividing by jj r n +1 d jj 0 throughout, squaring the result, multiplying by t 2 and adding up from n =0 to N , due to (34) and (37). 23
3.6 Global error behaviour As a consequence of the previous results, we have: Corolary 1: assume R1, R2, R3, R4, H1, H2, H3 and H4 hold assume also that, for n =0 :::N , u n +1 2 H r () , r 2 and p n +1 2 H s () , s 1 ,and that they are uniformly bounded in these spaces. Then, there exists a constant C > 0 independent of t and h such that, for smal l enough h and, if =1 , smal l enough t : jj e N +1 jj 2 0 + t N X n =0 jj e n +1 jj 2 1 + t 2 N X n =0 jj r n +1 jj 2 0 C t 2 + h 2 k (42) where k = min( r ; 1 sk v k q +1) . proof: this estimate follows from Theorems 1 and 3, Prop ositions 1 and 4, assumption (31), the regularity assumed of the semidiscrete solution ( u n +1 p n +1 ) and the approximating prop erties of the nite element spaces considered. Remark 5 :The condition t Ch 2 arises due to the pro of technique employed, which deals with the temp oral error rst and then the spatial error. However, according to the results of Corolary 1, accuracy considerations indicate that, when equal order interp olation of degree k is used, t should be of order h k for linear ( P 1 ) and bilinear ( Q 1 ) elements, one has k =1, so that assumption H4 is fullled. Even for quadratic ( P 2 )andbiquadratic ( Q 2 ) elements, one still has k = 2, making H4 acceptable. References 1] K. Boukir, Y. Maday, B. Metivet, E. Razandrakoto, A high-order characteristic/nite element metho d for the incompressible Navier-Stokes equations, International Journal for Numerical Metho ds in Fluids, 25 (1997) 1421{1454. 2] S.C. Brenner, L.R. Scott, The mathematical theory of nite element metho ds, Springer{Verlag, 1994. 3] F. Brezzi, R.S. Falk, Stability of higher{order Ho o d{Taylor metho ds, SIAM Journal on Numerical Analysis, 28 (1991) 581{590. 4] R. Co dina, A stabilized nite element metho d for generalized stationary incompressible ows, Computer Metho ds in Applied Mechanics and Engineering, (1999) submitted. 5] R. Co dina, J. Blasco, A nite element formulation for the Stokes pro-blem allowing equal velo city{pressure interp olation, Computer Metho ds in Applied Mechanics and Engineering, 143 , (1997) 373{391. 24
6] R. Co dina, J. Blasco, Analysis of a pressure stabilized nite elementapproximation of the stationary Navier{Stokes equations, Numerische Mathematik, (1999) to app ear. 7] R. Co dina, J. Blasco, Stabilized nite element metho d for the transient Navier{ Stokes equations based on apressure gradient pro jection, Computer Metho ds in Applied Mechanics and Engineering, (1999) to app ear. 8] Constantin, P.Foias, Navier{Stokes equations, Chicago Lectures in Mathematics, The University of Chicago Press, Chicago and London, 1988. 9] L.P. Franca, S.L. Frey, Stabilized nite element metho ds: II. The incompressible Navier{Stokes equations, Computer Metho ds in Applied Mechanics and Engineering, 99 (1992) 209{233. 10] L.P.Franca, T.J.R. Hughes, Convergence analysis of Galerkin least{squares methods for symmetric advective{diusive forms of the Stokes and incompressible Navier{Stokes equations, Computer Metho ds in Applied Mechanics and Engineering, 105 (1993) 285{298. 11] V. Girault, P.A. Raviart, Finite Element Approximation of the Navier{Stokes Equation, Springer{Verlag, New York, 1986. 12] J-L. Guermond, L- Quartap elle, On stability and convergence of pro jection meth- o ds based on pressure Poisson equation, International Journal for Numerical Methods in Fluids , 26 (1998) 1039{1053. 13] P. Hansb o, A. Szep essy, A velo city{pressure streamline diusion nite element metho d for the incompressible Navier{Stokes equations, Computer Metho ds in Applied Mechanics and Engineering, 84 (1990) 175{192. 14] J.G. Heywo o d, R. Rannacher , Finite element approximation of the nonstationary Navier{Stokes problem: Part 1: Regularity of solutions and second order error estimates for spatial discretization, SIAM Journal on Numerical Analysis, 19 (1982) 275{311. 15] J.G. Heywo o d, R. Rannacher , Finite element approximation of the nonstationary Navier{Stokes problem: Part 4: Error analysis for second{order time discretization, SIAM Journal on Numerical Analysis, 27 (1990) 353{384. 16] T.J.R. Hughes, Multiscale phenomena: Greens functions, subgrid scale mo dels, bubbles and the origins of stabilized metho ds, Computer Metho ds in Applied Mechanics and Engineering, 127 (1995) 387{401. 17] T.J.R. Hughes, L.P. Franca, M. Balestra, A new nite element formulation for computational uid dynamics: V. Circumventing the Babu^ska{Brezzi condition: a stable Petrov{Galerkin formulation of the Stokes problem accomo dating equal{ order interp olations, Computer Metho ds in Applied Mechanics and Engineering, 59 (1986) 85{99. 25