Méthodes de Galerkin avec stabilisation symétrique de la pression pour l'équation de Stokes transitoire: analyse de stabilité et de convergence

Résumé : On considère l'analyse de stabilité et convergence de méthodes d'éléments finis stabilisées pour l'équation de Stokes transitoire. L'analyse est valable pour des stabilisations symétriques et faiblements consistantes de la pression. Si l'approximation de la vitesse initiale est donnée par un pro jection de Ritz spécifique, dépendante de la stabilisation, on montre stabilité inconditionnelle et convergence optimale pour la vitesse et la pression, dans des normes naturelles. Par contre, pour des approximations arbitraires de la vitesse initiale, la stabilité de la pression est soumise à une condition entre les param'etres de discrétisation en espace et en temps.

Mots-clés : Équation de Stokes transitoire, méthode d'éléments finis, stabilisation symétrique de la pression, discrétisation en temps, projection de Ritz

Introduction

In this paper we consider stabilized finite element methods for the transient Stokes' problem. For standard pressure stabilizations of Petrov-Galerkin (PSPG) or least squares (GaLS) type, the analysis of time-discretization schemes is a difficult issue, unless a space-time approach is applied with a discontinuous Galerkin discretization in time. Indeed, for standard finite difference type time discretizations, the finite difference term must be included in the stabilization operator to assure consistency (see e.g. [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF][START_REF] Franca | Stabilized finite element methods. II. The incompressible Navier-Stokes equations[END_REF]). It has been shown, in [START_REF] Bochev | On stabilized finite element methods for the Stokes problem in the small time step limit[END_REF], that even for first order backward difference schemes (where the consistency does not improve convergence) this perturbs the stability of the numerical scheme when the time step is small, unless a condition between the space mesh size and the time step is verified,

δt ≥ Ch 2 , (1)
where δt denotes the time-step and h the space discretization parameter. For higher order schemes, such as Crank-Nicholson or second order backward differencing, the strongly consistent scheme appears to be unstable (see e.g [START_REF] Badia | On a multiscale approach to the transient Stokes problem[END_REF]). Similar initial time-step instabilities were observed, in [START_REF] Codina | Time dependent subscales in the stabilized finite element approximation of incompressible flow problems[END_REF], for the algebraic (static) sub-scale stabilization scheme applied to the Navier-Stokes' equations, and they were cured by including time dependent sub-scales. Our aim in this work is to consider a fairly large class of pressure stabilization methods and show that convergence of velocities and pressures, for the transient Stokes' problem, can be obtained without conditions on the space and time discretization parameters (like (1)), provided the initial data is chosen as a specific (method dependent) Ritz-projection (see e.g. [START_REF] Suri | Mixed variational principles for time-dependent problems[END_REF][START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]) onto a space of discretely divergence free functions. Discretely divergence free should here be interpreted in the sense of the stabilized method. If, on the other hand, the initial data is chosen as some interpolant that does not conserve the discrete divergence free character, the following condition

δt ≥ Ch 2k , (2)
with k the polynomial degree of the velocity approximation space, has to be respected in order to avoid pressure oscillations in the transient solution for small times.

Although the stability conditions (1) and ( 2) are similar, their nature is different. As mentioned above, if [START_REF] Becker | A finite element pressure gradient stabilization for the Stokes equations based on local projections[END_REF] fails to be satisfied pressure instabilities appear when dealing with non-discrete divergence free initial velocity approximation, but they are not related to the structure of the pressure stabilization. For residual based stabilization methods (PSPG, GaLS, etc.), on the other hand, the finite difference/pressure coupling of the stabilization perturbs the coercivity of the discrete pressure operator (see [START_REF] Bochev | On stabilized finite element methods for the Stokes problem in the small time step limit[END_REF]) unless condition [START_REF] Badia | On a multiscale approach to the transient Stokes problem[END_REF] is satisfied (irrespectively of the divergence free character of the initial velocity approximation).

The analysis carried out in this paper is valid for pressure stabilization operators that are symmetric and weakly consistent, as well as for standard methods using inf-sup stable velocity/pressure pairs, but it does not apply to residual based pressure stabilizations (PSPG, GaLS, etc.). In particular, space and time discretization commute (i.e. lead to the same fully discrete scheme) for the methods we analyze.

We prove unconditional stability of velocities and pressures and optimal convergence (in natural norms) when the initial data is chosen as a certain Ritz-type projection. In case a standard interpolation of the initial data is applied, an inverse parabolic CFL-type condition must be respected in order to maintain pressure stability for small time-steps. We give the full analysis only for the backward difference formula of order one, and we indicate how the analysis changes in the case of second order approximations in time. Indeed any A-stable implicit scheme is expected to yield optimal performance.

The remainder of the paper is organized as follows. In the next section we introduce the problem under consideration and some useful notation. The space and time discretized formulations are introduced in section §3. In paragraph §3.1, the space discretization is formulated using a general framework, we also discuss how some known pressure stabilized finite element methods enter this setting. The time discretization is performed in paragraph §3.2 using the BDF1, Crank-Nicholson or BDF2 schemes. Section §4 is devoted to the stability analysis of the resulting fully discrete formulations. The convergence analysis for the BDF1 scheme is carried out in section §5. We illustrate the theoretical results with some numerical experiments in section §6, using interior penalty stabilization of the gradient jumps. Finally, some conclusions are given in section §7.

Problem setting

Let Ω be a domain in R d (d = 2 or 3) with a polyhedral boundary ∂Ω. For T > 0 we consider the problem of solving, for u : Ω × (0, T ) -→ R d and p : Ω × (0, T ) -→ R, the following time-dependent Stokes problem:

         ∂ t u -ν∆u + ∇p = f , in Ω × (0, T ), ∇ • u = 0, in Ω × (0, T ), u = 0, on ∂Ω × (0, T ), u(•, 0) = u 0 , in Ω. (3)
Here, f : Ω × (0, T ) -→ R stands for the source term, u 0 : Ω -→ R d for the initial velocity and ν > 0 for a given constant viscosity. In order to introduce a variational setting for (3) we consider the following standard velocity and pressure spaces

V def = [H 1 0 (Ω)] d , H def = [L 2 (Ω)] d , Q def = L 2 0 (Ω), normed with v H def = (v, v) 1 2 , v V def = ν 1 2 ∇v H , q Q def = ν -1 2 q H ,
where (•, •) denotes the standard L 2 -inner product in Ω. Problem (3) can be formulated in weak form as follows: For all t > 0, find

u(t) ∈ V and p(t) ∈ Q such that      (∂ t u, v) + a(u, v) + b(p, v) = (f , v), a.e. in (0, T ), b(q, u) = 0, a.e. in (0, T ), u(•, 0) = u 0 , a.e. in Ω, (4)
INRIA for all v ∈ V , q ∈ Q and with

a(u, v) def = (ν∇u, ∇v), b(p, v) def = -(p, ∇ • v).
From these definitions, the following classical coercivity and continuity estimates hold:

a(v, v) ≥ v 2 V , a(u, v) ≤ u V v V , b(v, q) ≤ v V q Q , (5)
for all u, v ∈ V and q ∈ Q. It is known (see e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF])

that if f ∈ C 0 ([0, T ]; H) and that u 0 ∈ V ∩ H 0 (div; Ω) problem (4) admits a unique solution (u, p) in L 2 (0, T ; V ) × L 2 (0, T ; Q) with ∂ t u ∈ L 2 (0, T ; V ′ ).
Thorough this paper, C stands for a generic positive constant independent of the physical and discretization parameters.

Space and time discretization

In this section we discretize problem (4) with respect to the space and time variables. Symmetric pressure stabilized finite elements are used for the space discretization (paragraph §3.1), and some known A-stable schemes for the time discretization (paragraph §3.2).

Space semi-discretization: symmetric pressure stabilized formulations

Let {T h } 0<h≤1 denote a shape-regular family of triangulations of the domain Ω. For each triangulation T h , the subscript h ∈ (0, 1] refers to the level of refinement of the triangulation, which is defined by

h def = max K∈T h h K ,
with h K the diameter of K. In order to simplify the analysis we assume that the family of triangulations {T h } 0<h≤1 is quasi-uniform. For more precise information on the constraint on the mesh we refer to the analysis of the various finite element methods in the steady case, see paragraph §3.1.1.

In this paper, we let X k h and M l h denote, respectively, the standard spaces of continuous and (possibly) discontinuous piecewise polynomial functions of degree k ≥ 1 and

l ≥ 0 (k -1 ≤ l ≤ k), X k h def = v h ∈ C 0 (Ω) : v h|K ∈ P k (K), ∀K ∈ T h , M l h def = q h ∈ L 2 (Ω) : q h|K ∈ P l (K), ∀K ∈ T h .
For the approximated velocities we will use the space

[V k h ] d def = [X k h ∩ H 1 0 (Ω)] d
and for the pressure we will use either

Q l h def = M l h ∩ L 2 0 (Ω) or Q l h def = M l h ∩ L 2 0 (Ω) ∩ C 0 (Ω).
In order to stabilize the pressure we introduce a bi-linear form j : Q h × Q h -→ R satisfying the following properties:

• Symmetry: j(p h , q h ) = j(q h , p h ), ∀p h , q h ∈ Q l h , (6)
• Continuity:

|j(p h , q h )| ≤ j(p h , p h ) 1 2 j(q h , q h ) 1 2 ≤ C p h Q q h Q , ∀p h , q h ∈ Q l h , (7)
• Weak consistency:

j(Π l h q, Π l h q) 1 2 ≤ C h sp ν q sp,Ω , ∀q ∈ H s (Ω), (8)
with s p def = min{s, l, l + 1}, l ≥ 1 denoting the order of weak consistency of the stabilization operator, and Π l h :

Q -→ Q l h a given projection operator such that q -Π l h q Q ≤ C ν 1 2 h l+1 q l+1,Ω , (9)
for all q ∈ H l+1 (Ω).

Finally, we assume that there exists a projection operator

I k h : V -→ V k h
satisfying the following approximation properties:

v -I k h v H + hν -1 2 v -I k h v V ≤ C I h ru v ru,Ω , (10)
|b(q h , v -I k h v)| ≤ Cj(q h , q h ) 1 2 ν 1 2 h -1 (v -I k h v) H + v -I k h v V , (11)
for all v ∈ [H r (Ω)] d , r u = min{r, k + 1} and (q h , v h ) ∈ Q l h × [V k h ] d .
Our space semi-discretized scheme reads: For all t ∈ (0, T ), find (u

h (t), p h (t)) ∈ [V k h ] d × Q l h such that (∂ t u h , v h ) + a(u h , v h ) + b(p h , v h ) -b(q h , u h ) + j(p h , q h ) = (f , v h ), u h (0) = u 0 h , (12)
for all

(v h , q h ) ∈ [V k h ] d × Q l h and with u 0 h a suitable approximation of u 0 in [V k h ] d .
The following modified inf-sup condition states the stability of the discrete pressures in [START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF]. Lemma 3.1 There exists two constants C, β > 0, independent of h and ν, such that [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]Corollary 2.4] and [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF]

sup v h ∈[V k h ] d |b(q h , v h )| v h V + Cj(q h , q h ) 1 2 ≥ β q h Q , ( 13
)
for all q h ∈ Q l h . Proof. Let q h ∈ Q l h , from
there exists v q ∈ H 1 0 (Ω) such that ∇ • v q = ν -1 q h and I k h v q V ≤ C v q V ≤ C q h Q . ( 14
)
On the other hand, using [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], we have

q h 2 Q = b(q h , v q ) = b(q h , v q -I k h v q ) + b(q h , I k h v q ) ≤ Cj(q h , q h ) 1 2 ν 1 2 h -1 (v -I k h v) H + v -I k h v V + b(q h , I k h v q ) ≤ Cj(q h , q h ) 1 2 q h Q + b(q h , I k h v q ).

INRIA

We conclude the proof by dividing this last inequality by I k h v q V and using [START_REF] Burman | Edge stabilization for the generalized Stokes problem: a continuous interior penalty method[END_REF].

The above Lemma assures the well-posedness of problem [START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF]. This is stated in the following Theorem. Theorem 3.2 The discrete problem [START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF]

with u 0 h ∈ V div h,k def = {v h ∈ V k j : b(q h , v h ) = 0, ∀q h ∈ Q h ∩Ker j} has a unique solution (u h , p h ) ∈ C 1 (0, T ]; [V k h ] d × C 0 (0, T ]; Q k h .
To facilitate the analysis we introduce the following (mesh-dependent) seminorm, which is a norm for the velocity and a semi-norm for the pressure,

|||(v h , q h )||| 2 h def = v h 2 V + j(q h , q h ). ( 15
)
Remark 3.3 If the velocity/pressure finite-element pair V k h /Q k h is inf-sup stable, we can take j(•, •) def = 0 in (12), as usual. Obviously, this choice is compatible with hypothesis (6)-( 8) so that the results of this paper still apply. In particular, the relation (13) becomes the standard inf-sup condition between V k h and Q k h .

Examples

In this section we will review some of the most well-known pressure projection stabilization methods and discuss how they enter the abstract framework of the previous paragraph. For detailed results on analysis for the respective methods we refer the reader to the references considering the stationary case.

Recently, several different weakly consistent symmetric pressure stabilized finite element methods have been proposed. These methods take their origin in the work of Silvester [START_REF] Silvester | Optimal low order finite element methods for incompressible flow[END_REF] and the work by Codina and Blasco [START_REF] Codina | A finite element formulation for the Stokes problem allowing equal velocity-pressure interpolation[END_REF]. Further developments include the work by Becker and Braack [START_REF] Becker | A finite element pressure gradient stabilization for the Stokes equations based on local projections[END_REF] on local projection schemes, the extension of the interior penalty method, using penalization of gradient jumps, to the case of pressure stabilization by Burman and Hansbo [START_REF] Burman | Edge stabilization for the generalized Stokes problem: a continuous interior penalty method[END_REF] and the interpretation of these methods as minimal stabilization procedures by Brezzi and Fortin [START_REF] Brezzi | A minimal stabilisation procedure for mixed finite element methods[END_REF]. Similar approaches have been advocated in Bochev and Dohrmann in [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF] and a review of the analysis (with special focus on discontinuous pressure spaces and the Darcy problem) is given in [START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF].

The main idea underpinning all these methods is that, when using a velocitypressure space pair V h × Q h , the inf-sup stability constraint on the spaces may be relaxed by the addition of an operator penalizing the difference between the discrete pressure variable and its projection onto a subspace Qh ⊂ Q h , such that V h × Qh is inf-sup stable. The penalization may either act directly on the pressure, as in [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF][START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF], or on the gradient of the pressure, as in [START_REF] Becker | A finite element pressure gradient stabilization for the Stokes equations based on local projections[END_REF][START_REF] Codina | Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection[END_REF][START_REF] Burman | Edge stabilization for the generalized Stokes problem: a continuous interior penalty method[END_REF]. Generally speaking, the pressure approximation properties of the numerical scheme will be given by Qh , expressed in the weak consistency satisfied by the penalty operator. In some cases [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF][START_REF] Braack | Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method[END_REF][START_REF] Codina | Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales[END_REF], these stabilized methods may be extended to include high Reynolds number effects in the case of the Oseen's problem. The advantages and disadvantages of symmetric weakly consistent pressure stabilization methods compared to GaLS or PSPG approaches is discussed in a recent review paper [START_REF] Braack | Stabilized finite element methods for the generalized Oseen problem[END_REF].

The methods of Brezzi-Pitkäranta, Silvester and Bochev-Dohrmann

The original pressure stabilized finite element method was proposed by Brezzi and Pitkäranta in [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF]. Here, the velocity and pressure discrete spaces are chosen as the standard finite element space of piecewise affine continuous functions,

[V 1 h ] d × Q 1 h . The operator j(•, •) is given by j(p h , q h ) = h 2 ν ∇p h , ∇q h . (16)
A variant of this method was recently proposed by Dohrmann and Bochev, in [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF], using an equivalent stabilization operator, namely,

j(p h , q h ) = 1 ν (I -π 0 )p h , (I -π 0 )q h , (17)
where π 0 : Q -→ Q 0 h denotes the (element-wise) projection onto piecewise constants. Property [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] is verified after an integration by parts, with

I 1 h simply the Scott-Zhang interpolant onto [V 1 h ] d (see e.g. [31, 22]) b(q h , v -I 1 h v) = (∇q h , v -I 1 h v) ≤ j(q h , q h ) 1 2 h -1 ν 1 2 v -I 1 h v H + v -I 1 h v V .
One readily verifies that ( 7) and ( 8) hold. Moreover, in both cases, ( 16) and ( 17), the weak consistency property holds (with l = 1)

j(Π 1 h p, Π 1 h p) 1 2 ≤ C ν 1 2 h p 1,Ω ,
with Π1 h being, for instance, the L2 -projection onto Q 1 h (we could use instead the Clément [START_REF] Clément | Approximation by finite element functions using local regularization[END_REF] or Scott-Zhang interpolants). Indeed, for [START_REF] Codina | Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales[END_REF] we apply the H 1 -stability of the L 2 -projection (see e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF][START_REF] Crouzeix | The stability in L p and W 1 p of the L 2projection onto finite element function spaces[END_REF][START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF][START_REF] Boman | Estimates for the L 2 -projection onto continuous finite element spaces in a weighted L p -norm[END_REF]), whereas for [START_REF] Codina | A finite element formulation for the Stokes problem allowing equal velocity-pressure interpolation[END_REF] we add and subtract suitable terms (p and π 0 p) and use the approximation properties of π 0 and Π 1 h (see e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF]). As a result, our analysis for the time discretization is valid.

Another low order scheme, covered by the analysis, is the method which consists in using piecewise affine continuous velocities and element wise constants pressures, [START_REF] Hughes | A new finite element formulation for computational fluid dynamics. V. Circumventing the Babuška-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF]. Stability is obtained by the addition of the jump over element faces of the discontinuous pressure, namely,

[V 1 h ] d × Q 0 h , see e.g.
j(p h , q h ) = K∈T h ∂K\∂Ω h ν p h q h .
Here, q h denotes the jump of q h over the inter-element boundary, defined by,

q h (x) def = lim ǫ→0 (q h (x + ǫn) -q h (x -ǫn)) , ∀x ∈ F,
with n standing for a fixed, but arbitrary, normal to the internal face F . In this case, ( 11) is obtained after an integration by parts in the term b(q h , v -I 1 h v) and an element-wise trace inequality,

b(q h , v -I 1 h v) = - K ∂K\∂Ω q h (v -I 1 h v) • n ≤ j(q h , q h )
In addition, by taking, for instance, Π 0 h as the L 2 -projection onto Q 0 h , using an element-wise trace inequality and the approximation properties of Π 0 h (see e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF]), one also easily shows that the weak consistency property holds,

j(Π 0 h p, Π 0 h p) 1 2 = j (I -Π 0 h )p, (I -Π 0 h )p 1 2 ≤ C ν 1 2 h p 1,Ω , hence l = 1.
For details on the cases of stabilization of the pressure jumps only in macro elements, or the generalization to higher order finite element spaces of the Taylor-Hood family with discontinuous pressures, we refer to [START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF].

Orthogonal sub-scale stabilization

The orthogonal sub-scale stabilization was proposed by Codina and Blasco in [START_REF] Codina | A finite element formulation for the Stokes problem allowing equal velocity-pressure interpolation[END_REF]. Equal order (k = l ≥ 1) continuous approximation spaces are used for the velocities and the pressures.

Here the main idea is to penalize the difference between the pressure gradient and its projection onto the finite element space. This imposes the introduction of an auxiliary variable for the projection since it may not be localized and is given only implicitly. Hence, the stabilization operator is given by

j(p h , q h ) = h 2 ν (∇p h -π k h ∇p h ), ∇q h , where π k h : [L 2 (Ω)] d -→ [V k h ] d stands for the L 2 -projection onto [V k h ] d , which is given as the solution of the (global) problem (π k h ∇p h , ξ h ) = (∇p h , ξ h ), ∀ξ h ∈ [V k h ] d .
One may readily show that ( 7) and ( 8) hold. Disregarding for simplicity the boundary conditions, the projection operator 11) is here chosen also as the L 2 -projection onto [V k h ] d . This can be justified if boundary conditions are imposed weakly, for instance using Nitsche's method (see [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF]), and V k h includes the degrees of freedom on the boundary. Indeed, then we have

I k h = π k h of (
b(q h , v -I k h v) = (∇q h -π k h ∇q h , v -I k h v) ≤ j(q h , q h ) 1 2 h -1 ν 1 2 v -I k h v H + v -I k h v V . (18)
Finally, by taking

Π k h : Q -→ Q k h as the L 2 -projection onto Q k h
, adding and subtracting suitable terms (∇p and π k h ∇p) and using the approximation properties of π k h and Π k h (see e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF]), one readily verifies the weak consistency

j(Π k h p, Π k h p) 1 2 = h ν 1 2 (I -π k h )∇Π k h p 0,Ω ≤ C ν 1 2 h sp p sp,Ω ,
for all p ∈ H s (Ω) and with s p = min{k + 1, s}. In particular, l = l = k. The above analysis is hence valid also in this case (with some modifications of technical nature due to the weakly imposed boundary conditions).

Local projection stabilization

In the local projection stabilization, proposed in [START_REF] Becker | A finite element pressure gradient stabilization for the Stokes equations based on local projections[END_REF], stability is obtained by penalizing the projection of the gradient onto piecewise discontinuous functions defined on patches consisting of several elements, obtained by using hierarchic meshes, or by penalizing the gradient of the difference of the pressure and its projection on polynomials of lower polynomial order. The construction relies on the inf-sup stability of a velocity/pressure pair typically of mini-element character or of the Taylor-Hood family. Similar ideas were advocated in [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF].

The stabilization operator is written as

j(p h , q h ) = K ( h 2 ν κ∇p h , κ∇q h )
where κ is the so called fluctuation operator defined as κ def = I -πh where πh denotes a local projection operator onto either a polynomial of order k on a macro-patch consisting of three triangles (or four quadrilaterals) or a polynomial of order k -1 on the element. One may show that ( 11), ( 7) and ( 8) hold (for details on the construction of I k h see [START_REF] Becker | A finite element pressure gradient stabilization for the Stokes equations based on local projections[END_REF][START_REF] Braack | Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method[END_REF] and for general conditions on the finite element spaces and stabilization operators [START_REF] Matthies | A unified analysis for local projection stabilisations applied to the Oseen problem[END_REF]). In case we consider the projection πh onto polynomials of order k -1 the stabilization operator may be written as

j(p h , q h ) = K h 2 ν ∇(κp h ), ∇(κq h ) , (19)
or, equivalently, following [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF],

j(p h , q h ) = K 1 ν κ p h , κ q h .
In these latter cases, the condition ( 11) is obtained by choosing I k h as the Fortin interpolation operator associated to [V k h ] d × Qh where Qh is the space of continuous piecewise polynomial functions of order k -1. Clearly we then have [START_REF] Codina | Time dependent subscales in the stabilized finite element approximation of incompressible flow problems[END_REF] is treated in a similar fashion after an integration by parts. On the other hand, by taking Π l h : Q -→ Q l h as the L2 -projection operator onto Q l h and using approximation properties of Π l h and πh , we have

b(q h , v -I k h v) = b(κq h , v -I k h v) ≤ j(q h , q h ) since b(q h , v -I k h v) = 0 for all qh ∈ Qh . The form
j(Π l h p, Π l h p) 1 2 ≤ j (I -Π l h )p, (I -Π l h )p 1 2 + j(p, p) 1 2 ≤ C ν 1 2 h sp p sp,Ω , ∀p ∈ H s (Ω),
where s p = min{ l, s, l + 1} and l -1 denotes the polynomial order of the space on which the local projection is taken. Clearly if we project on polynomials of order k -1 the stabilization operator looses one order in the weak consistency, however the estimates remain optimal since we expect the velocities to be one order more regular than the pressure.

INRIA

Continuous interior penalty stabilization

The continuous interior penalty stabilization for the stationary Stokes' problem was proposed in [START_REF] Burman | Edge stabilization for the generalized Stokes problem: a continuous interior penalty method[END_REF] and generalized to the Oseen's problem in [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF]. It uses equal order continuous approximation spaces for velocities and pressures (k = l ≥ 1) and relies on the fact that the component of the pressure gradient orthogonal to the finite element space may be controlled by the gradient jumps using an interpolation estimate between discrete spaces. Indeed it was shown on [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF] that the following inequality holds

h(∇p h -ĩ∇p h ) 2 H ≤ K∈T h ∂K\∂Ω h 3 K ∇p h • n 2 , ( 20
)
for a certain Clément-type quasi-interpolation operator ĩ. This motivates the use of the pressure stabilization operator

j(p h , q h ) = K∈T h ∂K\∂Ω h 3 ν ∇p h • n ∇q h • n .
Clearly ( 7) and ( 8) are verified in this case. Moreover, [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] may be shown to hold if I k h is chosen to be the L 2 -projection onto [V k h ] d and boundary conditions are imposed weakly [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF]. To show the inequality we combine [START_REF] Codina | Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection[END_REF] with [START_REF] Crouzeix | The stability in L p and W 1 p of the L 2projection onto finite element function spaces[END_REF]. Finally, by taking Π k h as the

L 2 -projection onto Q k h , since C k h ∇p = 0 (with C k h the Clément interpolant onto [X k h ] d
), using an element-wise trace inequality, adding and subracting ∇p, and using the approximation properties of C k h and Π k h , one readily verifies (see [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF]Lemma 4.7]) the weak consistency

j(Π k h p, Π k h p) 1 2 ≤ C h ν 1 2 ∇Π k h p -C k h ∇p 0,Ω ≤ C ν 1 2 h sp p sp,Ω , ∀p ∈ H s (Ω),
with s p = min{k + 1, s}, so that l = l = k. We refer to [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF] for the details on the technical issue related to the weak imposition of the boundary conditions using Nitsche's method.

The Ritz projection operator

For the purpose of the stability and convergence analysis below we introduce the Ritz projection operator

S k,l h : [H 1 (Ω)] d × L 2 (Ω) -→ V k h × Q l h . For each (u, p) ∈ [H 1 (Ω)] d ×L 2 (Ω), the projection S k,l h (u, p) def = P k h (u, p), R l h (u, p) ∈ [V k h ] d × Q l h is defined as the unique solution of a(P k h (u, p), v h ) + b(R l h (u, p), v h ) =a(u, v h ) + b(p, v h ), -b(q h , P k h (u, p)) + j(R l h (u, p), q h ) =0, (21)
for all (v h , q h ) ∈ [V k h ] d × Q l h .
Problem ( 21) is well-posed thanks to the inf-sup condition [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF], in particular, we have the following a priori stability estimate:

||| P k h (u, p), R l h (u, p) ||| 2 h ≤ C u 2 V + p 2 Q , (22)
with C > 0 a constant independent of h and ν.

Finally, we have the following approximation result.

Lemma 3.4 Let (u, p) ∈ C 1 ([0, T ], [H r (Ω) ∩ H 1 0 (Ω))] d ∩ H 0 (div; Ω) × H s (Ω)) with r ≥ 2 and s ≥ 1.
The following error estimate for the projection S k,l h holds with α = 0, 1:

||| ∂ α t (u -P k h (u, p)), ∂ α t R l h (u, p) ||| h ≤ C ν 1 2 h ru-1 ∂ α t u ru,Ω + ν -1 2 h sp ∂ α t p sp,Ω , p -R l h (u, p) Q ≤ C ν 1 2 h ru-1 u ru,Ω + ν -1 2 h sp p sp,Ω ,
for all t ∈ [0, T ] and with r u def = min{r, k + 1}, s p def = min{s, l, l + 1} and C > 0 independent of ν and h. Moreover, provided the domain Ω is sufficiently smooth and if l ≥ 1, there also holds

∂ α t (u -P k h (u, p)) H ≤ Ch|||(∂ α t (u -P k h (u, p)), ∂ α t R l h p)||| h . (23)
Proof. For simplicity we here use the notation 15), the V -coercivity of a(•, •) (5), and the orthogonality provided by [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF] we have

u h def = P k h (u, p) and p h def = R l h (u, p). From (
||| u h -I k h u, p h -Π l h p |||
Finally, using ( 5) and ( 11), we have that

||| u h -I k h u, p h -Π l h p ||| 2 h ≤ u -I k h u V + p -Π l h p Q u h -I k h u V +C ν 1 2 h -1 (u -I k h u) H + u -I k h u V + j(Π l h p, Π l h p) 1 2 j(p h -Π l h p, p h -Π k l p) 1 2 .
We obtain the estimation for the velocity (α = 0) using the approximation properties of I k h and Π l h (( 10) and ( 9)), and the weak consistency (8) of the stabilizing term j(•, •). The convergence for the time derivative (α = 1) is obtained in a similar fashion after time derivation of [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF].

For the pressure estimate we use the generalized inf-sup condition ( 13) and the orthogonality provided by [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF]. We then have

β Π l h p -p h Q ≤ sup v h ∈[V k h ] d b(Π l h p -p h , v h ) v h V + Cj(Π l h p -p h , Π l h p -p h ) 1 2 ≤ sup v h ∈[V k h ] d b(Π k h p -p, v h ) -a(u -u h , v h ) v h V + Cj(Π l h p -p h , Π l h p -p h ) 1 2 .
We conclude using the continuity of a(

), p(t n )) in [V k h ] d × Q l h .

Backward Euler scheme (BDF1)

By introducing the first order backward difference quotient

Du n+1 h def = u n+1 h -u n h δt ,
our first fully discrete scheme reads as follows: For 0

≤ n ≤ N -1, find (u n+1 h , p n+1 h ) ∈ [V k h ] d × Q l h such that ( Du n+1 h , v h ) + a(u n+1 h , v h ) + b(p n+1 h , v h ) -b(q h , u n+1 h ) + j(p n+1 h , q h ) = (f (t n+1 ), v h ), (24)
for all

(v h , q h ) ∈ V k h × Q l h and with u 0 h a suitable approximation of u 0 in [V k h ] d .

Crank-Nicholson scheme

Let us consider now the scheme given by: For 0 ≤ n ≤ N -1, find (u n+1 h , p [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] uniquely determines u n+1 h , since u 0 h is given. For the pressure however neither p n+1 h nor p n h are used in [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. Therefore, by working with p

n+ 1 2 h ) ∈ [V k h ] d × Q l h such that ( Du n+1 h , v h ) + a(u n+ 1 2 h , v h ) + b(p n+ 1 2 h , v h ) -b(q h , u n+ 1 2 h ) + j(p n+ 1 2 h , q h ) = (f n+ 1 2 , v h ), ( 25
) for all (v h , q h ) ∈ [V k h ] d × Q l h , where u n+ 1 2 h def = 1 2 (u n+1 h + u n h ) and u 0 h is a suitable approximation of u 0 in [V k h ] d . Remark 3.5 Note that
n+ 1 2 h
as pressure variable, we do not need to provide an initial condition for the pressure. On the other hand, we do not have an approximation of p n+1 h unless one is constructed by extrapolation.

Second order backward difference (BDF2)

Finally, by considering the second order backward difference quotient

Du n+1 def = 1 2δt (3u n+1 h -4u n h + u n-1 h ),
we obtain the BDF2 scheme: For 1

≤ n ≤ N -1, find (u n+1 h , p n+1 h ) ∈ [V k h ] d ×Q l h such that ( Du n+1 h , v h ) + a(u n+1 h , v h ) + b(p n+1 h , v h ) -b(q h , u n+1 h ) + j(p n+1 h , q h ) = (f n+1 , v h ), ( 26
) RR n°6321 for all (v h , q h ) ∈ [V k h ] d × Q l h , and (u 1 h , p 1 h ) ∈ [V k h ] d × Q l h
given by the first step of backward Euler scheme [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF].

Stability

In this section we analyze the stability properties of the fully discrete schemes introduced in §3.2. For the sake of simplicity, full details will be only given for the backward scheme [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF]. Nevertheless, in §4.2, we will discuss how the results extend to the second order time stepping schemes: Crank-Nicholson and BDF2.

First order A-stable scheme

The next result provides the unconditional stability of the velocity. It also provides a uniform estimate for the pressure, in terms of the discrete velocity time derivative. Theorem 4.2 points out the role of the initial velocity approximation on the stability of the velocity time derivative approximations. Finally, Corollary 4.3 states the (conditional or unconditional) stability of the pressure depending on the choice of the initial velocity approximation.

Theorem 4.1 Let u 0 h be a given H-stable approximation of u 0 in [V k h ] d and {(u n h , p n h )} N n=1
be the solution of the fully discrete problem [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF]. The following estimate holds for 1 ≤ n ≤ N : Proof. Taking v h = u n+1 h and q h = p n+1 h in [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF], using the coercivity of the bi-linear form, the Cauchy-Schwarz's inequality and the Poincaré's inequality, we have

u n h 2 H + n-1 m=0 δt|||(u m+1 h , p m+1 h )||| 2 h ≤ C u 0 2 H + C 2 P ν n-1 m=0 δt f (t m+1 ) 2 H , n-1 m=0 δt p m+1 h 2 Q ≤ C β 2 n-1 m=0 δt |||(u m+1 h , p m+1 h )||| 2 h + ν -1 Du m+1 h 2 H + ν -1 f (t m+1 ) 2 H . (27
( Du n+1 h , u n+1 h ) + 1 2 |||(u n+1 h , p n+1 h )||| 2 h ≤ C 2 P 2ν f (t n+1 ) 2 H . (28)
Now, recalling that

( Du n+1 h , u n+1 h ) = 1 2 D u n+1 h 2 H + 1 2δt u n+1 h -u n h 2 H , (29)
we have

D u n+1 h 2 H + |||(u n+1 h , p n+1 h )||| 2 h ≤ C 2 P ν f (t n+1 ) 2 H ,
leading to, after summation over 0 ≤ m ≤ n -1,

u n h 2 H + n-1 m=0 δt|||(u m+1 h , p m+1 h )||| 2 h ≤ u 0 h 2 H + C 2 P ν n-1 m=0 δt f (t m+1 ) 2 H .

INRIA

For the pressure estimate, from ( 13), ( 24) (with q h = 0) and the Poincaré's inequality, we have

β p n+1 h Q ≤ C |||(u n+1 h , p n+1 h )||| h + ν -1 2 Du n+1 h H + ν -1 2 f (t n+1 ) H ,
which completes the proof.

The next theorem states some a priori estimates of the approximations of the velocity time derivative.

Theorem 4.2 Let {(u n h , p n h )} N n=1
be the solution of the fully discrete problem [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF].

• If u 0 ∈ [H 1 (Ω)] d and u 0 h = P k h (u 0 , 0), the following estimate holds for 1 ≤ n ≤ N : n-1 m=0 δt Du m+1 h 2 H + |||(u n h , p n h )||| 2 h ≤ C u 0 2 V + n-1 m=0 δt f (t m+1 ) 2 H . ( 30
)
• If u 0 ∈ [H r (Ω) ∩ H 1 0 (Ω)] d ∩ H 0 (div; Ω), r ≥ 2, and u 0 h = I k h u 0 , the following estimate holds for 1 ≤ n ≤ N : n-1 m=0 δt Du m+1 h 2 H + |||(u n h , p n h )||| 2 h ≤ C u 0 2 V + νh 2(ru-1) p 1 h 2 Q + u 0 2 ru,Ω + n-1 m=0 δt f (t m+1 ) 2 H , (31)
with r u def = min{k + 1, r}.

Proof. For 0 ≤ n ≤ N -1, by taking v h = Du n+1 h and q h = 0 in (24) and using the Cauchy-Schwarz inequality, we have

1 2 Du n+1 h 2 H + a(u n+1 h , Du n+1 h ) + b(p n+1 h , Du n+1 h ) = 1 2 f (t n+1 ) 2 H , (32)
On the other hand, for 1 ≤ n ≤ N -1, testing (24) at the time levels n and n + 1 with v h = 0 and

q h = p n+1 h we have b(p n+1 h , u n+1 h ) = j(p n+1 h , p n+1 h ), b(p n+1 h , u n h ) = j(p n h , p n+1 h ). ( 33
)
Therefore, by subtracting these equalities and using the bi-linearity of j(•,

•), we obtain b(p n+1 h , Du n+1 h ) = j( Dp n+1 h , p n+1 h ), (34)
for 1 ≤ n ≤ N -1. It then follows, from [START_REF] Silvester | Optimal low order finite element methods for incompressible flow[END_REF], that

1 2 Du n+1 h 2 H + a(u n+1 h , Du n+1 h ) + j(p n+1 h , Dp n+1 h ) ≤ 1 2 f (t n+1 ) 2 H . (35)
On the other hand, using the symmetry and bi-linearity of a(•, •) and j(•, •), we have

a(u n+1 h , Du n+1 h ) = 1 2 Da(u n+1 h , u n+1 h ) + δt 2 a( Du n+1 h , Du n+1 h ), j(p n+1 h , Dp n+1 h ) = 1 2 Dj(p n+1 h , p n+1 h ) + δt 2 j( Dp n+1 h , Dp n+1 h ).
Hence,

Du n+1 h 2 H + D a(u n+1 h , u n+1 h ) + j(p n+1 h , p n+1 h ) ≤ f (t n+1 ) 2 H , for 1 ≤ n ≤ N -1.
After multiplication by δt and summation over 1

≤ m ≤ n-1, it follows that n-1 m=1 δt Du m+1 h 2 H + |||(u n h , p n h )||| 2 h ≤ |||(u 1 h , p 1 h )||| 2 h + n-1 m=1 δt f (t m+1 ) 2 H . ( 36
)
In order to highlight the impact of the initial velocity approximation on the stability of the time derivative, we consider now the first time level (n = 0) of [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF]. By testing with v h = Du 1 h , q h = 0, after multiplication by 2δt and using the symmetry and bi-linearity of a(•, •), we get

δt Du 1 h 2 H + a(u 1 h , u 1 h ) -a(u 0 h , u 0 h ) + 2δtb(p 1 h , Du 1 h ) ≤ δt f (t 1 ) 2 H . (37)
If the initial velocity approximation is given in terms of the Ritz-projection,

u 0 h = P k h (u 0 , 0) with u 0 ∈ [H 1 (Ω)] d , by setting p 0 h def = R l h (u 0 , 0) it follows that (33) also holds for n = 0. Therefore, b(p 1 h , Du 1 h ) = j( Dp 1 h , p 1 h ), (38)
Thus, from the symmetry and bi-linearity of j(•, •) and (37), we have

δt Du 1 h 2 H + |||(u 1 h , p 1 h )||| 2 h ≤ |||(u 0 h , p 0 h )||| 2 h + δt f (t 1 ) 2 H . (39)
Estimate ( 30) is obtained by adding this last inequality to (36) and using the stability of the Ritz-projection [START_REF] Ern | Theory and practice of finite elements[END_REF],

|||(u 0 h , p 0 h )||| 2 h ≤ C u 0 2 V .

If the initial velocity approximation is given in terms of a general interpolant,

u 0 h = I k h u 0 with u 0 ∈ [H r (Ω) ∩ H 1 0 (Ω)] d ∩ H 0 (div; Ω), equality (
(p 1 h , Du 1 h ) = 1 δt j(p 1 h , p 1 h ) -p 1 h , ∇ • (I k h u 0 -u 0 ) ≥ 1 δt j(p 1 h , p 1 h ) - C I δt νh 2(ru-1) p 1 h 2 Q + u 0 2 ru,Ω , (40)
with r u def = min{k + 1, r}. As a result, from (37), it follows that

δt Du 1 h 2 H + |||(u 1 h , p 1 h )||| 2 h ≤ a(u 0 h , u 0 h ) + C I νh 2(ru-1) p 1 h 2 Q + u 0 2 ru,Ω + δt f (t 1 ) 2
H . We conclude the proof by adding this equality to (36) and using the stability of the Ritz-projection.

The next Corollary closes the problem of the stability of the pressures by combining the results of Theorems 4.1 and 4.2.

INRIA Corollary 4.3 Let {(u n h , p n h )} N n=1
be the solution of the fully discrete problem [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF]. Then:

• If u 0 ∈ [H 1 (Ω)] d and u 0 h = P k h (u 0 , 0), the following estimate holds for 1 ≤ n ≤ N : n-1 m=0 δt p m+1 h 2 Q ≤ C β 2 ν u 0 2 V + C β 2 n-1 m=0 δt |||(u m+1 h , p m+1 h )||| 2 h + ν -1 f (t m+1 ) 2 H . (41) • If u 0 ∈ [H r (Ω) ∩ H 1 0 (Ω)] d ∩ H 0 (div; Ω), r ≥ 2, u 0 h = I k h u 0 and 2C I β 2 h 2(ru-1) ≤ δt, (42)
the following estimate holds for 1 ≤ n ≤ N :

n-1 m=0 δt p m+1 h 2 Q ≤ C β 2 ν u 0 2 V + u 0 2 ru,Ω + C β 2 n-1 m=0 δt |||(u m+1 h , p m+1 h )||| 2 h + ν -1 f (t m+1 ) 2 H , (43)
with r u def = min{k + 1, r}.

Proof. Estimate (43) is a direct consequence of Theorem 4.1 and estimate [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]. On the other hand, from Theorem 4.1 and estimate (31), we have

β 2 δt -C I h 2(ru-1) p 1 h 2 Q + β 2 n-1 m=1 δt p m+1 h 2 Q ≤ C ν u 0 2 V + u 0 2 ru,Ω +C n-1 m=0 δt |||(u m+1 h , p m+1 h )||| 2 h + ν -1 f (t m+1 ) 2 H ,
which combined with the stability condition (42) leads to (43).

A few observations are now in order. Corollary 4.3 states the unconditional stability of the pressure provided the initial velocity approximation u 0 h is given in terms of the Ritz-projection operator [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF]. In the general case, i.e. whenever u 0 h does not satisfy a discrete divergence free condition (as u 1 h does), only conditional stability can be guaranteed. As a matter of fact, from the stability condition (42), pressure instabilities are expected for very small time steps. This issue will be illustrated by numerical experiments in section §6.

Finally, let us mention that residual based stabilization methods, such as PSPG and GaLS, combined with finite differences time discretization schemes, are known to give rise to pressure instabilities in the small time step limit, see [START_REF] Bochev | On stabilized finite element methods for the Stokes problem in the small time step limit[END_REF][START_REF] Codina | Time dependent subscales in the stabilized finite element approximation of incompressible flow problems[END_REF]. Indeed, it has been shown in [START_REF] Bochev | On stabilized finite element methods for the Stokes problem in the small time step limit[END_REF] that the finite difference/pressure coupling of the stabilization perturbs the coercivity of the discrete pressure operator unless a condition of the type

Ch 2 ≤ δt, (44)
is satisfied. It is worth emphasizing that, although the stability conditions (44) and (42) are somehow similar, their nature is different. Actually, the instabilities anticipated by Corollary 4.3 are related to the discrete divergence free character of the initial velocity approximation, but not to the structure of the pressure stabilization j(•, •).

Second order A-stable schemes

In this paragraph we discuss how the results of Theorems 4.1 and 4.2 and Corollary 4.3 extend to the second order time stepping schemes Crank-Nicholson and BDF2.

Crank-Nicholson

The following theorem summarizes the resulting stability estimates.

Theorem 4.4 Let u 0 h be a given H-stable approximation of u 0 in [V k h ] d and {(u n h , p n h )} N n=1
be the solution of the discrete scheme [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. Then the following estimate holds for 1 ≤ n ≤ N :

u n h 2 H + n-1 m=0 δt|||(u m+ 1 2 h , p m+ 1 2 h )||| 2 h ≤ C u 0 2 H + C 2 P ν n-1 m=0 δt f (t m+ 1 2 ) 2 H . Moreover, if u 0 ∈ [H 1 (Ω)] d and u 0 h = P k h (u 0 , 0), the following estimate holds for 1 ≤ n ≤ N : n-1 m=0 δt p m+ 1 2 h 2 Q ≤ C β 2 ν u 0 2 V + C β 2 n-1 m=0 δt |||(u m+ 1 2 h , p m+ 1 2 h )||| 2 h + ν -1 f (t m+ 1 2 ) 2 H .
On the other hand, if

u 0 ∈ [H r (Ω) ∩ H 1 0 (Ω)] d ∩ H 0 (div; Ω), r ≥ 2, u 0 h = I k
h u 0 and the stability condition (42) is satisfied, the following estimate holds for

1 ≤ n ≤ N : n-1 m=0 δt p m+ 1 2 h 2 Q ≤ C β 2 ν u 0 2 V + u 0 2 ru,Ω + C β 2 n-1 m=0 δt |||(u m+ 1 2 h , p m+ 1 2 h )||| 2 h + ν -1 f (t m+ 1 2 ) 2 H .
Proof. The first estimate, corresponding to Theorem 4.1, holds by taking

v h = u n+ 1 2 h
and q h = p . We first note that, since

u n+1 h -u n h = 2 u n+ 1 2 h -u n h , u n h = 2u n-1+ 1 2 h -u n-1 h , INRIA
we have

u n+1 h -u n h = 2u n+ 1 2 h + 4 n l=1 (-1) l u n-l+ 1 2 h -(-1) n 2u 0 h , for 0 ≤ n ≤ N -1.
Therefore, from ( 25) and using the bi-linearity of j(•, •), we get b(p

n+ 1 2 h , u n+1 h -u n h ) =j 2p n+ 1 2 h + 4 n l=1 (-1) l p n-l+ 1 2 h , p n+ 1 2 h -2(-1) n b p n+ 1 2 h , u 0 h . ( 45
)
On the other hand, for 0 ≤ n ≤ N -1, we introduce the following change of variables (or extrapolation)

1 2 p n+1 h + p n h def = p n+ 1 2 h
, with p 0 h ∈ Q l h to be specified later on. By inserting this expression in (45) we obtain b(p

n+ 1 2 h , u n+1 h -u n h ) =j p n+1 h -p n h , p n+ 1 2 h + 2(-1) n j p 0 h , p n+ 1 2 h -b p n+ 1 2 h , u 0 h . (46)
If u 0 h = P k h (u 0 , 0) and we choose p 0 21) 2 it follows that the last term in (46) cancels. Thus, we have b(p

h def = R l h (u 0 , 0), from (
n+ 1 2 h , Du n+1 h ) = j Dp n+1 h , p n+ 1 2 h = 1 2 Dj p n+1 h , p n+1 h ,
which corresponds to the Crank-Nicholson's counterpart of [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]. Finally, when the initial velocity approximation is given in terms of a general interpolant, u 0 h = I k h u 0 with u 0 ∈ [H r (Ω) ∩ H 1 0 (Ω)] d ∩ H 0 (div; Ω), we take p 0 h def = 0. Therefore, from (46) and using an approximation argument (as in (40)), we get b(p

n+ 1 2 h , Du n+1 h ) = 1 2 Dj p n+1 h , p n+1 h - 2 δt (-1) n b p n+ 1 2 h , u 0 h ≥ 1 2 Dj p n+1 h , p n+1 h - 2C I δt νh 2(ru-1) p n+ 1 2 h 2 Q + u 0 2 ru,Ω , (47)
which leads to the stability condition (42). The rest of the proof follows with minor modifications.

Remark 4.5 By comparing the proofs of Corollary 4.3 and the previous Theorem we can notice that, if the initial velocity approximation is not discretely divergence free, the stability condition (42) has to be satisfied at each time level when using the Crank-Nicholson scheme (due to (47)), whereas for the Backward Euler scheme that condition is only needed at the first time step (thanks to [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF] and (40)).

Second order backward difference

The following theorem summarizes the resulting stability estimates.

Theorem 4.6 Let u 0 h be a given H-stable approximation of u 0 in [V k h ] d , (u 1 h , p 1
h ) be the corresponding first time step of the Backward Euler scheme [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF] and {(u n h , p n h )} N n=2 be the solution of the discrete scheme [START_REF] Hecht | FreeFem++ v. 2.11. User's Manual[END_REF]. Then, the following estimate holds for 2 ≤ n ≤ N :

u n h 2 H +2 n-1 m=1 δt|||(u m+1 h , p m+1 h )||| 2 h ≤ C u 0 2 H + u 1 h 2 H + 2C 2 P ν n-1 m=1 δt f (t m+1 ) 2 H . Moreover, if u 0 ∈ [H 1 (Ω)] d and u 0 h = P k h (u 0 , 0), the following estimate holds for 2 ≤ n ≤ N : n-1 m=1 δt p m+1 h 2 Q ≤ C β 2 ν u 0 2 V + |||(u 1 h , p 1 h )||| 2 h + C β 2 n-1 m=1 δt |||(u m+1 h , p m+1 h )||| 2 h + ν -1 f (t m+1 ) 2 H .
On the other hand, if

u 0 ∈ [H r (Ω) ∩ H 1 0 (Ω)] d ∩ H 0 (div; Ω), r ≥ 2, u 0 h = I k h u 0 and the stability condition (42) is satisfied, the following estimate holds for 2 ≤ n ≤ N : n-1 m=1 δt p m+1 h 2 Q ≤ C β 2 ν u 0 2 V + u 0 2 ru,Ω + |||(u 1 h , p 1 h )||| 2 h + C β 2 n-1 m=1 δt |||(u m+1 h , p m+1 h )||| 2 h + ν -1 f (t m+1 ) 2 H .
Proof. The first estimate, corresponding to Theorem 4.1, holds by taking v h = u n+1 h and q h = p n+1 h in [START_REF] Hecht | FreeFem++ v. 2.11. User's Manual[END_REF] and applying the standard identity

(3a -4b + c)a = 1 2 a 2 -b 2 + (2a -b) 2 -(2b -c) 2 + (a -2b + c) 2 . ( 48
)
providing the numerical dissipation of the BDF2 scheme.

Since the pressure estimate is based here on the control of the time derivative, Du n+1 h , we take v h = Du n+1 h and q h = 0 in [START_REF] Hecht | FreeFem++ v. 2.11. User's Manual[END_REF]. In particular, for the coupling term b p n+1 h , Du n+1 h , using ( 26) and (48), we have

b p n+1 h , Du n+1 h = j Dp n+1 h , p n+1 h ≥ 1 4 D j(p n+1 h , p n+1 h ) + j(2p n+1 h -p n h , 2p n+1 h -p n h ) , (49)
for 2 ≤ n ≤ N -1, which corresponds to the BDF2 counterpart of [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]. On the other hand for n = 1, from ( 26) and ( 24), we obtain

b p 2 h , Du 2 h = 1 2δt 3j(p 2 h , p 2 h ) -4j(p 1 h , p 2 h ) + b(p 2 h , u 0 h ) . (50)

INRIA

If the initial velocity approximation is given in terms of the Ritz-projection,

u 0 h = P k h (u 0 , 0), it follows that b(p 2 h , u 0 h ) = j(p 0 h , p 2 h ), with p 0 h def = R l h (u 0 , 0). Thus, (50) reduces to b p 2 h , Du 2 h = j( Dp 2 h , p 2 h ),
so that (49) holds true also for n = 1. Finally, if the initial velocity approximation is given in terms of a general interpolant, u 0 h = I k h u 0 , we apply an approximation argument (as in (40)). Hence, from (50)

b p 2 h , Du 2 h ≥ 1 2δt 3j(p 2 h , p 2 h ) -4j(p 1 h , p 1 h ) -C I νh 2(ru-1) p 2 h 2 Q + u 0 2 ru,Ω (51)
which is the BDF2 counterpart of (40) and leads to the stability condition (42). The rest of the proof follows with minor modifications. Remark 4.8 When the initial velocity approximation is not discretely divergence free, the stability condition (42) has to be satisfied twice when using BDF2: At the first time step (according to (51)) and at the Backward Euler initialization (see Theorem 4.2).

Convergence

In this paragraph we provide optimal convergence error estimates for the discrete formulation [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF], the Backward Euler scheme.

Theorem 5.2 concerns the convergence for the velocity and gives an estimate for the pressure in terms of the error in the velocity time derivative. Theorem 5.3 closes the question of optimal convergence of the pressure by providing an optimal error estimate for the time derivative, provided the exact pressure is smooth. Finally, Theorem 5.4 provides an improved L ∞ ((0, T ), H) estimate that justifies the initialization of the BDF2 scheme with a Backward Euler step.

The following result expresses the modified Galerkin orthogonality in terms of the consistency error in space and time. Lemma 5.1 (Consistency error) Let (u, p) be the solution of (3) and {(u n h , p n h )} 0≤n≤N the solution of [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF].

Assume that u ∈ C 0 [0, T ]; V ∩ C 1 (0, T ]; H and p ∈ C 0 (0, T ]; Q . Then, for 0 ≤ n ≤ N -1, there holds Du(t n+1 ) -Du n+1 h , v h + a u(t n+1 ) -u n+1 h , v h + b(p(t n+1 ) -p n+1 h , v h ) -b(q h , u(t n+1 ) -u n+1 h ) = j p n+1 h , q h + Du(t n+1 ) -∂ t u(t n+1 ), v h , for all (v h , q h ) ∈ [V k h ] d × Q l h . Theorem 5.2 Assume that u ∈ H 1 (0, T ; [H r (Ω)] d ) ∩ H 2 (0, T ; [L 2 (Ω)] d ) and p ∈ C 0 ((0, T ]; H s (Ω)) with r ≥ 2 and s ≥ 1, and set u 0 h ∈ [V k h ] d a given
approximation of u 0 . Then the following estimate holds for 1 ≤ n ≤ N :

u n h -u(t n ) 2 H + n-1 m=0 δt|||(u m+1 h -u(t m+1 ), p m+1 h )||| 2 h ≤ I k h u 0 -u 0 h 2 H + Ch 2ru u 2 C 0 ([t1,tn];H ru (Ω)) + ν -1 ∂ t u 2 L 2 (0,tn;H ru (Ω)) +C δt 2 ν ∂ tt u 2 L 2 (0,tn;H) + h 2sp ν t n p 2 C 0 ([t1,tn];H sp (Ω)) +νh 2(ru-1) t n u 2 C 0 ([t1,tn];H ru (Ω)) n-1 m=0 δt p m+1 h -p(t m+1 ) 2 Q ≤ C 1 + 1 β 2 h 2sp ν t n p 2 C 0 ([t1,tn];H sp (Ω)) + C β 2 n-1 m=0 δt |||(u m+1 h -u(t m+1 ), p m+1 h )||| 2 h + ν -1 ∂ t u(t m+1 ) -Du m+1 h 2 H ,
with C > 0 a positive constant independent of h, δt and ν.

Proof. The error estimate for the velocity follows standard energy arguments, and for the pressure we use the modified inf-sup condition [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF]. We start by decomposing the velocity and pressure error using, respectively, the projection

I k h and Π h . This yields, u(t n+1 ) -u n+1 h = u(t n+1 ) -I k h u(t n+1 ) θ n+1 π + I k h u(t n+1 ) -u n+1 h θ n+1 h = θ n+1 π + θ n+1 h , p(t n+1 ) -p n+1 h = p(t n+1 ) -Π l h p(t n+1 ) y n+1 π + Π l h p(t n+1 ) -p n+1 h y n+1 h = y n+1 π + y n+1 h .
(52) The first term θ n+1 π can bounded using approximation [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF]. In order to estimate θ n+1 h we first note, using [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] and the coercivity of the bi-linear form a(•,

•) + j(•, •) 1 2 D θ n+1 h 2 H + |||(θ n+1 h , y n+1 h )||| 2 h ≤ ( Dθ n+1 h , θ n+1 h ) + |||(θ n+1 h , y n+1 h )||| 2 h ≤ ( Dθ n+1 h , θ n+1 h ) + a(θ n+1 h , θ n+1 h ) + b(y n+1 h , θ n+1 h ) -b(y n+1 h , θ n+1 h ) + j(y n+1 h , y n+1 h ) T n+1 1 . (53)
In addition, using (52) we have

T n+1 1 = -( Dθ n+1 π , θ n+1 h ) -a(θ n+1 π , θ n+1 h ) + j(Π l h p(t n+1 ), y n+1 h ) -b(y n+1 π , θ n+1 h ) + b(y n+1 h , θ n+1 π ) + ( Du(t n+1 ) -Du n+1 h , θ n+1 h ) + a(u(t n+1 ) -u n+1 h , θ n+1 h ) + b(p(t n+1 ) -p n+1 h , θ n+1 h ) -b(y n+1 h , u(t n+1 ) -u n+1 h ) -j(p n+1 h , y n+1 h ).
By the modified Galerkin orthogonality (Lemma 5.1) this expression reduces to

T n+1 1 = -( Dθ n+1 π , θ n+1 h ) + Du(t n+1 ) -∂ t u(t n+1 ), θ n+1 h -a(θ n+1 π , θ n+1 h ) + j(Π l h p(t n+1 ), y n+1 h ) -b(y n+1 π , θ n+1 ) + b(y n+1 h , θ n+1 π ). ( 54
)

INRIA

Now, using the Cauchy-Schwarz and the Poincaré's inequalities and [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], we have

T n+1 1 ≤ Du(t n+1 ) -∂ t u(t n+1 ) H + Dθ n+1 π H T n+1 2 C p ν 1 2 |||(θ n+1 h , y n+1 h )||| + θ n+1 π V + y n+1 π Q + ν 1 2 h -1 θ n+1 π H + j Π l h p(t n+1 ), Π l h p(t n+1 ) 1 2 |||(θ n+1 h , y n+1 h )|||. ( 55
)
The term T n+1 2 can be treated, in a standard way (see e.g [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]), using a Taylor expansion and the Cauchy-Schwarz inequality, which yields

T n+1 2 ≤ 1 δt tn+1 tn (δt ∂ tt u(s) H + ∂ t θ π (s) H ) ds ≤ δt 1 2 ∂ tt u(s) L 2 ((tn,tn+1);H) + δt -1 2 ∂ t θ π L 2 ((tn,tn+1);H) . (56)
Thus, from (55), using Young's inequality, it follows that

T n+1 1 ≤ 1 2 |||(θ n+1 h , y n+1 h )||| 2 +C C 2 P ν δt ∂ tt u 2 L 2 ((tn,tn+1);H) + δt -1 ∂ t θ π 2 L 2 ((tn,tn+1);H) + θ n+1 π 2 V + y n+1 π 2 Q + ν h -1 θ n+1 π 2 H + j Π l h p(t n+1 ), Π l h p(t n+1 )
By inserting this expression in (53), multiplying the resulting expression by 2δt and after summation over 0 ≤ m ≤ n -1 we obtain

θ n h 2 H + n-1 m=0 δt|||(θ m+1 h , y m+1 h )||| 2 h ≤ θ 0 h 2 H +C δt 2 ν -1 ∂ tt u 2 L 2 (0,tn;H) +ν -1 ∂ t θ π 2 L 2 (0,tn;H) + n-1 m=0 δt θ m+1 π 2 V + y m+1 π 2 Q +ν h -1 θ m+1 π 2 H +j Π l h p(t m+1 ), Π l h p(t m+1 ) .
Finally, the velocity error estimate is obtained using approximation [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF] and the consistency of the pressure stabilization [START_REF] Bramble | On the stability of the L 2 projection in H 1 (Ω)[END_REF], which yields

θ n h 2 H + n-1 m=0 δt|||(θ m+1 h , y m+1 h )||| 2 h ≤ θ 0 h 2 H + C δt 2 ν ∂ tt u 2 L 2 (0,tn;H) + h 2ru ν ∂ t u 2 L 2 (0,tn;H ru (Ω)) + νh 2(ru-1) n-1 m=0 δt u(t m+1 ) 2 ru,Ω + h 2sp ν n-1 m=0 δt p(t m+1 ) 2 sp,Ω ,
For the pressure error estimate we first note that, from (52), it suffices to control y n+1 h 0,Ω . To this aim, we use the modified inf-sup condition [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF]:

β y n+1 h Q ≤ sup v h ∈[V k h ] d |b(y n+1 h , v h )| v h V + Cj(y n+1 h , y n+1 h ) 1 2 . ( 57
)

RR n°6321

From (52) we get

b(y n+1 h , v h ) = -b(y n+1 π , v h ) + b(p(t n+1 ) -p n+1 h , v h ).
The first term can be bounded, using the continuity of b(•, •), [START_REF] Boman | Estimates for the L 2 -projection onto continuous finite element spaces in a weighted L p -norm[END_REF], this yields

b(y n+1 π , v h ) ≤ y n+1 π Q v h V .
On the other hand, using the modified Galerkin orthogonality (Lemma 5.1 with

q h = 0) we have b(p(t n+1 ) -p n+1 h , v h ) = -a u(t n+1 ) -u n+1 h , v h -∂ t u(t n+1 ) -Du n+1 h , v h ≤C|||(u(t n+1 ) -u n+1 h , 0)||| h v h V + ∂ t u(t n+1 ) -Du n+1 h H v h H
As a result, from the above estimations we have:

β y n+1 h Q ≤ C y n+1 π Q + |||(u(t n+1 ) -u n+1 h , y n+1 h )||| h + C P ν 1 2 ∂ t u(t n+1 )-Du n+1 h H .
Therefore,

β 2 n-1 m=0 δt y m+1 h 2 Q ≤C n-1 m=0 δt y m+1 π 2 Q + |||(u(t m+1 ) -u m+1 h , y m+1 h )||| 2 h + ν -1 ∂ t u(t m+1 ) -Du m+1 h 2 H ,
and we conclude using approximation and the error estimate for the velocity.

We close the problem of the pressure convergence by providing an error estimate for the time derivative of the velocity.

Theorem 5.3 Under the assumptions of Theorem 5.2, assuming that p ∈ C 0 ([0, T ]; H s (Ω)),

u 0 ∈ V ∩ H 0 (div; Ω) and if u 0 h def = P k h (u 0 , 0), for 1 ≤ n ≤ N , we have n-1 m=0 δt Du m+1 h -∂ t u(t m+1 ) 2 H +||| P k h u(t n ), p(t n ) -u n h , R l h u(t n ), p(t n ) -p n h ||| 2 h ≤ C δt 2 ∂ tt u 2 L 2 (0,T ;H) + h 2ru ∂ t u 2 L 2 (0,T ;H ru (Ω)) + C h 2sp ν p(0) 2 sp,Ω .
Proof. In order to provide an optimal error estimate, we decompose the error in terms of the Ritz-projection operator ( 21)

u(t n+1 ) -u n+1 h = u(t n+1 ) -P k h (u(t n+1 ), p(t n+1 )) θ n+1 π + P k h (u(t n+1 ), p(t n+1 )) -u n+1 h θ n+1 h = θ n+1 π + θ n+1 h , p(t n+1 ) -p n+1 h = p(t n+1 ) -R l h (u(t n+1 ), p(t n+1 )) y n+1 π + R l h (u(t n+1 ), p(t n+1 )) -p n+1 h y n+1 h = y n+1 π + y n+1 h . (58)

INRIA

Using the triangle inequality we then have

n-1 m=0 δt ∂ t u(t m+1 ) -Du m+1 h 2 H ≤ C n-1 m=0 δt ∂ t u(t m+1 ) -Du(t m+1 ) 2 H + Dθ m+1 π 2 H + Dθ m+1 h 2 H . ( 59
)
For the first term we proceed as in (56) using a Taylor expansion, which yields

∂ t u(t n+1 ) -Du(t n+1 ) 2 H ≤ δt 1 2 ∂ tt u(s) L 2 ((tn,tn+1);H) .
For the second term we have,

Dθ n+1 π H = 1 δt tn+1 tn ∂ t θ π (s) H ds ≤ δt -1 2 ∂ t θ π L 2 ((tn,tn+1);H). (60)
Finally, for the third term we use the modified Galerkin orthogonality (Lemma 5.1 with q h = 0) and the definition of the Ritz-projection ( 21) to obtain

Dθ n+1 h 2 H + a(θ n+1 h , Dθ n+1 h ) + b(y n+1 h , Dθ n+1 h ) = -( Dθ n+1 π , Dθ n+1 h )-a(θ n+1 π , Dθ n+1 h )-b(y n+1 π , Dθ n+1 h )+( Du(t n+1 )-∂ t u(t n+1 ), Dθ n+1 h ) = -( Dθ n+1 π , Dθ n+1 h ) + ( Du(t n+1 ) -∂ t u(t n+1 ), Dθ n+1 h ).
Young's inequality yields:

1 2 Dθ n+1 h 2 H + a(θ n+1 h , Dθ n+1 h ) + b(y n+1 h , Dθ n+1 h ) ≤ C Dθ n+1 π 2 H + Du(t n+1 ) -∂ t u(t n+1 ) 2
H . In addition, for 0 ≤ n ≤ N , testing [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF] at the time level n with v h = 0 we have

b q h , P k h (u(t n ), p(t n )) = j R l h (u(t n ), p(t n )), q h . (61)
On the other hand, for 1 ≤ n ≤ N , testing [START_REF] Freund | On weakly imposed boundary conditions for second order problems[END_REF] at the time level n with v h = 0 and since, by definition, u 0

h def = P k h (u 0 , 0) we have b(q h , u n h ) = j(p n h , q h ), (62)
for all q h ∈ Q l h and 0 ≤ n ≤ N and where we have defined p 0 h def = R l h (u 0 , 0). As a result, from (61)-(62), we have b(q h , θ n h ) = j(y n h , q h ), for all q h ∈ Q l h and 0 ≤ n ≤ N . We therefore have, for 0

≤ n ≤ N -1 b(y n+1 h , Dθ n+1 h ) = j( Dy n+1 h , y n+1 h ).
On the other hand, using the symmetry of a and j, we have

a(θ n+1 h , Dθ n+1 h ) = 1 2 Da(θ n+1 h , θ n+1 h ) + δt 2 a( Dθ n+1 h , Dθ n+1 h ), j(y n+1 h , Dy n+1 h ) = 1 2 Dj(y n+1 h , y n+1 h ) + δt 2 j( Dy n+1 h , Dy n+1 h
).

RR n°6321

So that

1 2 Dθ n+1 h 2 H + 1 2 D a(θ n+1 h , θ n+1 h ) + j(y n+1 h , y n+1 h ) ≤ Dθ n+1 π 2 H + Du(t n+1 ) -∂ t u(t n+1 ) 2 H .
So that, after multiplication by 2δt and summation over 0

≤ n ≤ N -1 we have n-1 m=0 δt Dθ m+1 h 2 H + |||(θ n h , y n h )||| 2 h ≤ |||(θ 0 h , y 0 h )||| 2 h + C n-1 m=0 δt Dθ m+1 π 2 H + Du(t m+1 ) -∂ t u(t m+1 ) 2 H . ( 63
)
For the initial terms we use the linearity of the Ritz-projection and its approximation properties (Lemma 3.4) to obtain

|||(θ 0 h , y 0 h )||| 2 h = |||(P k h (0, p(0)), R l h (0, p(0)))||| 2 h ≤ C ν h 2sp p(0) 2 sp,Ω .
Therefore, using (60) and (56), we have

n-1 m=0 δt Dθ m+1 h 2 H + |||(θ n h , y n h )||| 2 h ≤ C h 2sp ν p(0) 2 sp,Ω + δt 2 ∂ tt u 2 L 2 (0,tn;H) , for 1 ≤ n ≤ N .
Finally, for completeness, we here give a result of optimal convergence in the L ∞ ((0, T ), H)-norm. For this we assume that the domain Ω is such that the optimal convergence in the H-norm holds for the Ritz-projection (see Lemma 3.4). This result is of importance since it shows that the initialization of the BDF2 method using one BDF1 step is justified (i.e. we keep error optimality in time). Theorem 5.4 Assume that the domain Ω is sufficiently smooth so that the Hestimate (23) holds. Assume also that u ∈ H

1 (0, T ; [H r u (Ω)] d )∩H 2 (0, T ; [L 2 (Ω)] d ), p ∈ C 0 ([0, T ]; H sp (Ω)) with r u ≥ 2, s p ≥ 1, u 0 ∈ V ∩ H 0 (div; Ω) and u 0 h def = P k h (u 0 , 0). Then the following estimate holds for 1 ≤ n ≤ N : u(t n ) -u n h H ≤ C ν 1 2 h ru u 0 ru,Ω + h sp+1 p(0) sp,Ω + h ru ∂ t u L 1 (0,tn;H ru (Ω)) + δt ∂ tt u L 1 (0,tn;H) .
Proof. Since the proof is similar to that of Theorem 5.3 and we will only give the outline. Let θ n+1 h and y n+1 h be defined as in (58). From (53) and (54), it follows that

( Dθ n+1 h , θ n+1 h ) ≤ -( Dθ n+1 π , θ n+1 h ) + ( Du(t n+1 ) -∂ t u(t n+1 ), θ n+1 h
).

INRIA

Applying now the Cauchy-Schwarz inequality, we have

θ n+1 h H ≤ θ n h H + δt Dθ n+1 π H + Du(t n+1 ) -∂ t u(t n+1 ) H ,
and by summation over n, we get

θ n h H ≤ θ 0 h H + n-1 m=0 δt Dθ n+1 π H + Du(t n+1 ) -∂ t u(t n+1 ) H , for 1 ≤ n ≤ N .
The first term in the right hand side can be estimated using Lemma 3.4 since, by definition,

θ 0 h = P k h (u(0), p(0)) -u 0 h = P k h (u 0 , p(0)) -P k h (u 0 , 0) = P k h (0, p(0)). (64)
Finally, for the finite difference consistency terms we use a standard argument (see e.g. [34, Theorem 1.5, Page 14]).

Remark 5.5 From (64), one could pretend to initialize the time-stepping procedure with u 0 h = P k h (u 0 , p(0)) (as in [START_REF] Suri | Mixed variational principles for time-dependent problems[END_REF], for instance). In practice, however, the initial pressure is unknown, so that the choice u 0 h = P k h (u 0 , 0) is more convenient. Lemma 3.4 shows that we can preserve optimality while keeping this choice (see also [START_REF] Boffi | Analysis of finite element approximation of evolution problems in mixed form[END_REF]).

Remark 5.6 Note that the above convergence proofs only use stability, Galerkin orthogonality and the truncation error of the finite difference time approximation scheme. Hence the extension to the second order Crank-Nicholson or BDF2 schemes is straightforward. In particular we recall that the estimate of Lemma 5.4 shows that the initialization using one BDF1 step does not make the convergence deteriorate, provided the solution is sufficiently smooth under the first timestep. Indeed, for smooth solutions we expect ∂ tt u L 1 (0,δt;H) to be O(δt) and hence the global convergence will be second order in spite of the intial low order perturbation.

Numerical experiments

In this section we will consider some numerical examples using the continuous interior penalty stabilization (CIP), described in §3.1.1. We present computations demonstrating the optimal convergence using finite element spaces consisting of quadratic functions, for the space discretization, and first and second order backward difference formulas and the Crank-Nicholson scheme for the time discretization. We also verify numerically that, for small time steps, the pressure is unstable for initial data that are not discretely divergence free. All computations have been performed using FreeFem++ [START_REF] Hecht | FreeFem++ v. 2.11. User's Manual[END_REF].

Convergence rate in time

We consider problem (3) in two dimensions, Ω = [0, 1] × [0, 1] and T = 1, with non-homogeneous boundary conditions. The right hand side f and the We choose the initial data either as the Lagrange interpolant, u 0 h = I k h u 0 , or as the Ritz-projection, u 0 h = P k h (u 0 , 0). In Figure 2 we have reported the convergence history (in space) of the pressure error, at the first time step, using P 1 /P 1 finite elements for different time step sizes. The pressure instability for small time steps is illustrated in Figure 4(a), where the initial velocity approximation is given in terms of the Lagrange interpolant. Indeed, we can observe that the pressure error has the right convergence rate in space, but it grows when the time step is decreased. On the other hand, as shown in Figure 4(b), the instability is eliminated when the initial velocity approximation is provided by the Ritz-projection, as stated in Corollary 4.3. In this case the error remains bounded (dominated by the space discretization) while reducing the time step size. Similar results are found with P 2 /P 2 finite elements, as shown in Figure 3. In particular, we can notice, from Figures 2(a) and 3(a), that for quadratic approximations the pressure instability shows up only for very small time steps. As a matter of fact, condition (42) is less restrictive for quadratic than for affine velocity approximations of smooth initial data. Finally, some pressure contours are reported in Figure 4 for the Lagrange interpolation, and in Figure 5 for the Ritz-projection. The pressure degradation is clearly visible in Figure 4, whereas with the Ritz-projection initialization (Figure 5) the pressure remains unconditionally stable.

Conclusion

In this paper we have proved unconditional stability and optimal error estimates, in natural norms, for pressure stabilized finite element approximations of the transient Stokes' problem. It should be noted that the extension of the present results to mixed formulations of the Poisson problem is straightforward. We have shown that for small initial time-steps the use of a pressure stabilization dependent Ritz-projection, for the initial data, is essential to avoid pressure instabilities, unless a condition between time and space discretization parameters is satisfied. From the analysis, we also conclude that a second order scheme (e.g. BDF2) can be initialized (without optimality loss) using a first step with BDF1, provided that the Ritz-projection [START_REF] Dohrmann | A stabilized finite element method for the Stokes problem based on polynomial pressure projections[END_REF] is used for the initial data.

It is interesting to note that for low order elements the weakly consistent stabilization operators still yield optimal convergence in time when used with a second order scheme. However, in case SUPG type stabilization is used for the convective term the convergence order in time will be lost unless full consistency is guaranteed in the stabilization term. This is why SUPG type stabilizations prompt space time finite element formulations with discontinuous approximation in time.

Somme the methods described in paragraph §3.1.1, on the other hand, may be extended to the case of the Oseen's equations, handling all Reynolds numbers, by applying the same type of stabilizing term for the convection (see [START_REF] Burman | Continuous interior penalty finite element method for Oseen's equations[END_REF][START_REF] Braack | Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method[END_REF][START_REF] Codina | Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales[END_REF][START_REF] Braack | Stabilized finite element methods for the generalized Oseen problem[END_REF], for details).

  ) with C P > 0 the Poincaré's constant.

  38) does not hold in general. Instead, we can use an approximation argument to obtain b

n+ 1 2 h 2 h

 22 in[START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. The pressure estimate requires an a priori bound of the discrete velocity time derivative. As in Theorem 4.2, such an estimate can be obtained by taking v h = Du n+1 h and q h = 0 in[START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], for 0 ≤ n ≤ N -1. The main difference, with respect to the proof of Theorem 4.2, arises in the treatment of the coupling term b, in the Crank-Nicholson scheme incompressibility is enforced on u n+ 1 instead of u n+1 h

Remark 4 . 7 A

 47 bound for the Backward Euler initialization terms u 1 h H and |||(u 1 h , p 1 h )||| h , appearing in the above estimates, is provided by Theorems 4.1 and 4.2 with n = 1.

Figure 1 :

 1 Figure 1: Convergence history in time: P 2 / P 2 CIP stabilized finite elements

Figure 2 :Figure 3 :

 23 Figure 2: Convergence history: P 1 /P 1 finite elements

(

  a) δt = 10 -1 (b) δt = 10 -6

Figure 4 :

 4 Figure 4: Pressure contour lines with P 2 /P 2 finite elements in a 40 × 40 mesh: u 0 h = I 2 h u 0

INRIA

  (a) δt = 10 -1 (b) δt = 10 -6

Figure 5 :

 5 Figure 5: Pressure contour lines with P 2 /P 2 finite elements in a 40 × 40 mesh: u 0 h = P 2 h (u 0 , 0)

1 1

 11 First order A-stable scheme . . . . . . . . . . . . . . . . . . . . . 14 4.2 Second order A-stable schemes . . . . . . . . . . . . . . . . . . Convergence rate in time . . . . . . . . . . . . . . . . . . . . . . 27 6.2 Behavior in the small time step limit . . . . . . . . . . . . . . . . 28 7 Conclusion

  In this paragraph we discretize[START_REF] Burman | Pressure projection stabilizations for Galerkin approximations of Stokes' and Darcy's problem[END_REF] with respect to the time variable. To this aim, we will use some known A-stable time discretization schemes for ODE's.Let N ∈ N * be given, we consider a uniform partition {[t n , t n+1 ]} 0≤n≤N -1 ,

	with t n	def = nδt, of the time interval of interest [0, T ] with time-step size δt	def =
	T /N . The discrete pair (u n h , p n h ) stands for an approximation of (u(t n

•, •) and b(•, •), approximability, the weak consistency of j(•, •) and the previous error estimate. For a proof of the optimality in the H-norm we refer e.g. to [13, Theorem 4.14]. INRIA 3.2 Fully discrete formulation: time discretization

RR n°6321

h -1 ν 1 2 v -I 1 h v H + v -I 1 h v V .INRIA

h -1 ν 1 2 v -I k h v H + v -I k h v V ,

h =a(u -I k h u, u h -I k h u) + b(p -Π l h p, u h -I k h u) + b(p h -Π l h p, u -I k h u) + j(Π l h p, p h -Π l h p).

boundary and initial data are chosen in order to ensure that the exact solution is given by: u(x, y, t) = g(t)

sin(πx -0.7) sin(πy + 0.2) cos(πx -0.7) cos(πy + 0.2) , p(x, y, t) = g(t) sin(x) cos(y) + (cos(1) -1) sin(1) , with g(t) = 1 + t 5 + e -t 10 + sin(t). In order to illustrate the convergence rate in time of the discrete solution, we have used quadratic approximations in space and a mesh parameter h = 0.01. In this case, the stability condition (42) is always satisfied for the range of time steps considered. Thus, the choice of the Lagrange interpolant or of the Ritz-projection as approximation of the initial velocity give similar results.

In figures 1(a)-(c) we report the convergences of the errors for the velocities ( • L ∞ (0,T ;L 2 (Ω)) ) and the pressures ( • L 2 (0,T ;L 2 (Ω)) ) for the first order backward difference formula (BDF1), Crank-Nicholson and the second order backward difference formulas (BDF2). In all the numerical examples both the velocities and the pressures converge at the optimal rate (O(δt) for BDF1 and O(δt 2 ) for Crank-Nicholson and BDF2). The BDF2 scheme was initialized using one step of BDF1.

Behavior in the small time step limit

In this paragraph we illustrate the impact of the initial velocity approximation on the approximate pressures for small time steps. For non discrete divergence free initial approximations, a pressure instability is predicted by Corollary 4.3 unless condition (42) is satisfied. In other words, pressure instabilities are expected for very small time steps.

We consider problem (3) in two dimensions and with non-homogeneous boundary conditions. We set Ω = [0, 1] × [0, 1], and the right hand side f and the boundary data are chosen in order to ensure that the exact (steady) solution is given by: u(x, y, t) = sin(πx -0.7) sin(πy + 0.2) cos(πx -0.7) cos(πy + 0.2) , p(x, y, t) = sin x cos y + (cos(1) -1) sin(1).

This numerical experiment is, in some degree, motivated by the work reported in [START_REF] Bochev | On stabilized finite element methods for the Stokes problem in the small time step limit[END_REF] (see also [START_REF] Codina | Time dependent subscales in the stabilized finite element approximation of incompressible flow problems[END_REF]), where pressure instabilities, of a different nature, are illustrated for pressure stabilizations involving residuals of the partial differential equation (e.g PSPG and GaLS). Indeed, the time derivative involved in the residual perturbs the coercivity of the space semi-discrete operator, which leads to pressure instabilities for (sufficiently) small time steps (see [START_REF] Bochev | On stabilized finite element methods for the Stokes problem in the small time step limit[END_REF]). Let us emphasize that, according to section §4, such kind of instabilities do not appear here, in particular since the CIP pressure stabilization (and the other examples of paragraph §3.1.1) is consistent without introducing the time derivative.

For different initial velocity approximations, we compare the behavior of the error in the pressure after one time step of the Backward Euler scheme, i.e.,


Éric Lubat

Silvano Dal

A Short Overview on Diagnosability of Patterns in Timed Petri Net

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A Short Overview on Diagnosability of Patterns in Timed Petri Net Éric Lubat and Silvano Dal Zilio LAAS-CNRS, Université de Toulouse, CNRS, INSA, Toulouse, France eric.lubat@laas.fr

Diagnosability is a basic property of Discrete Event Systems (DES) that relates to the observability of concealed events. Basically, it means that every failure (a distinct instance of unobservable event) can be eventually detected after a nite number of observations. In this work, we are interested by the diagnosability of systems modelled using labelled Time Petri nets (TPN), an extension of Petri nets in which we can associate timing constraints to transitions [START_REF] Merlin | A study of the recoverability of computing systems[END_REF]. This means that we take into account the date at which events are observed and that we want to detect failures in a bounded time. We are also interested by the detection of patterns of events (sequence of observable events that are part of some given regular language) instead of the occurrence of a single fault.

The problem of fault diagnosis was introduced by Sampath et al [START_REF] Sampath | Diagnosability of discrete-event systems[END_REF] and is well studied in the context of discrete event systems [START_REF] Cassandras | Introduction to Discrete Event Systems[END_REF]10]. The problem has also been studied on timed models, see for instance the work of Tripakis [START_REF] Tripakis | Fault diagnosis for Timed Automata[END_REF] with Time Automata. In these works, diagnosability is often reduced to properties on the trace language of systems and their composition. These problems are much more dicult when we need to take timing constraints into account. For instance, one of the main problem when considering Time Petri nets is that it is not possible to build the synchronous composition of two transitions when they have non-trivial time constraints. Therefore it is not possible to syntactically dene the intersection of the trace languages of two TPN.

We tackle the diagnosability problem using model-checking techniques, which means that we reduce diagnosability to the problem of checking the validity of a temporal logic formula on a nite-state model. This approach relies on a new class of TPN, called Product TPN (PTPN), introduced in [START_REF] Lubat | A state class construction for computing the intersection of time petri nets languages[END_REF]. The idea is to enrich TPN with a new synchronization operator (×) that can be used to compute the (language) intersection of two or more TPN. We show in [START_REF] Lubat | A state class construction for computing the intersection of time petri nets languages[END_REF] that it is possible to extend the classical State Class Graph construction to this model.

This work makes two main contributions. First, we describe how to use PTPN to decide the diagnosability of a system. Our approach extends the classical techniques based on the twin-plant construction, which relies on the composition of two copies of (the model of ) a system. Next, we show that we can also dene a notion of pattern diagnosability as a straightforward extension of the single fault problem. We conclude with a short discussion on future work.

Diagnosability of Single Faults. There are two main problems when checking timed extensions of Petri nets. First, the state space of a net is typically innite when we use a dense time model, that is when time delays can be arbitrarily small. Therefore we need to work on an abstraction of the transition system.

Second, there is no natural way to dene the (structural) composition of two transitions that have non-trivial time constraints.

A solution to the rst problem was proposed by Berthomieu and Menasche in [START_REF] Berthomieu | An enumerative approach for analyzing time Petri nets[END_REF], where the authors dene a state space abstraction based on state classes.

A state class captures a convex set of constraints on the time at which transitions can re. This approach is used in several model-checking tools, such as Tina [START_REF] Berthomieu | The tool TINAconstruction of abstract state spaces for Petri nets and time Petri nets[END_REF] for instance, and makes it possible to check LTL properties on bounded TPN. To solve the second problem, we propose to use Product PTPN, a new model where we can dene groups of transitions, in a net, that should re concurrently, meaning together, and at the same date. The main idea behind this product operator is to mimic the behaviour of a synchronous product between transitions.

Given two TPN, N and N , we can dene a PTPN N × N that has exactly the (observable) behaviours that are in both N and N .

Armed with PTPN, we can easily dene the twin-plant construction of a net N (see e.g. [START_REF] Yoo | Polynomial-time verication of diagnosability of partially observed discrete-event systems[END_REF]), as the composition of two copies of the net, say N.1 × N.2, on the set of their observable labels. In the following, we consider that failures are transitions on a common unobservable label, say f . As a consequence, a system is diagnosable when we cannot nd a pair of execution traces such that: [START_REF] Berthomieu | An enumerative approach for analyzing time Petri nets[END_REF] they are compatible (observable events occur in the same orders and at the same date); and (2) only one of them exhibits a failure. This would entail that we have an execution, with the occurrence of a failure, that is indistinguishable from a nominal one. This case (often referred to as the existence of a critical pair) can be directly translated into traces in the twin-plant. In particular, system N is diagnosable if all the (maximal) executions of the twin-plant N.1 × N.2 satises the LTL formula (♦f.1) ⇔ (♦f.2). This support the idea that diagnosability can be reduced to checking some temporal properties on the twin-plant model, an idea that was already exploited in [START_REF] Tripakis | Fault diagnosis for Timed Automata[END_REF][START_REF] Yoo | Polynomial-time verication of diagnosability of partially observed discrete-event systems[END_REF].

In this context, we use the general assumption that systems are ultimately observable, meaning that they do not block and that, on every execution, we will always eventually nd an observable event (after a bounded number of transitions). Hence a system is diagnosable when pairs of compatible traces have the same failures, or if one of them block, meaning that it is not possible to extend them with some compatible event. We can further simplify this property by using the inherent symmetry of the problem. Theorem 1. a TPN N is diagnosable if and only if all the maximal executions of the PTPN N.1 × N.2 satisfy the LTL formula (f.1 ⇒ ♦(dead ∨ f.2)).

Theorem 1 provides a constructive method for checking the diagnosability of a TPN, as long as the net is bounded. We have implemented this method inside a tool called TWINA. Diagnosability of Patterns. It is possible to extend this method to the observability of patterns of events [START_REF] Gougam | Diagnosability analysis of patterns on bounded labeled prioritized Petri nets[END_REF]. We dene a pattern as a labelled Petri net, P , representing the set of behaviours (sequences of labels) that we want to detect. We also identify a distinguished place in P , say found , that corresponds to the detection of the pattern. We say that the pattern is detected in the net N when the place found is marked in N × P . Drawing a parallel with the case of single failure, we say that a TPN N is diagnosable for the pattern P if we cannot nd a pair of compatible (and maximal) executions in N × P such that the pattern is detected in one and not the other. Again, we propose to reduce this diagnosability problem to a model-checking problem on the twin-plant, this time on the product of the system and the pattern. Theorem 2. Assuming some well-formedness conditions on the pattern P , TPN N is diagnosable for pattern P i all maximal executions of the PTPN (N.1 × P.1) × (N.2 × P.2) satisfy the LTL formula (found .1 ⇒ ♦(dead ∨ found .2)).

This time, we need to impose some extra restrictions. Intuitively, patterns are very much like observers in observer-based verication methods. Like with observers, we want to make sure that P does not interfere with the system it interacts with. For example, it should not prevents the execution of the system.

To this end, we impose three well-formedness conditions, which are sucient conditions for Theorem 2 to hold: (1) the pattern must to be total (it should never block or delay a transition on one of its label, so there is always at least one transition for any given label); (2) the pattern must be deterministic (the same observations should lead to the same states, so there is at most one possible transition for any given label); and (3) the net P is bounded and place found is a trap (once it is marked, it stays marked). Future Work. We give a very broad overview of a constructive method to decide the diagnosability of single failures and patterns of failures in TPN. Our presentation emphasize the parallel structure between these two properties. In the future, we hope to benet from this interconnection in order to dene more elaborate observability properties. For instance we would like to study the link with the notion of opacity and see if we can come up with an equivalent formulation with patterns.