Continuous Interior Penalty Finite Element Method for Oseen's Equations

In this paper we present an extension of the continuous interior penalty method of Douglas and Dupont [Interior penalty procedures for elliptic and parabolic Galerkin methods, in Computing Methods in Applied Sciences, Lecture Notes in Phys. 58, Springer-Verlag, Berlin, 1976, pp. 207–216] to Oseen’s equations. The method consists of a stabilized Galerkin formulation using equal order interpolation for pressure and velocity. To counter instabilities due to the pressure/velocity coupling, or due to a high local Reynolds number, we add a stabilization term giving L2-control of the jump of the gradient over element faces (edges in two dimensions) to the standard Galerkin formulation. Boundary conditions are imposed in a weak sense using a consistent penalty formulation due to Nitsche. We prove energy-type a priori error estimates independent of the local Reynolds number and give some numerical examples recovering the theoretical results.

1. Introduction.The construction of finite element methods for the incompressible Navier-Stokes equations that are robust and accurate for a wide range of Reynolds numbers remains a challenging problem.The standard Galerkin method requires the fulfillment of the inf-sup or Babuska-Brezzi condition, which leads to the need for formulations using mixed interpolations (see [7,27]).From the computational point of view it is, however, more practical to use equal order interpolation for the velocity and pressure spaces, which requires that stability is imposed in some other fashion.One possibility is to construct stabilized finite element methods where some terms are added to the standard Galerkin formulation in order to enhance the stability properties of the method.To be useful the method must also be stable with respect to the convective terms and give sufficient control of the incompressibility condition.
A favored approach has been to stabilize both the velocities and the pressure using the streamline upwind Petrov-Galerkin (SUPG) method originally proposed by Brooks and Hughes in [9].This method was first analyzed for the Navier-Stokes equations in a velocity vorticity formulation by Johnson and Saranen in [32], and then in a pressure velocity formulation by Hansbo and Szepessy in [29], by Franca and Frey in [24], and by Tobiska and Verfürth in [38].The SUPG method owes its success to the unified treatment of velocities and pressures.It allows for a priori error estimates that are independent of the Reynolds number and has been used extensively in practice with good results.Nevertheless, the SUPG method has some undesirable features: • artificial boundary conditions on velocities and pressure are introduced; • artificial nonsymmetric terms are introduced; • the least squares term introduces nonphysical pressure/velocity couplings; • the least squares term makes mass lumping impossible and the choice of timestepping methods limited; most clear-cut from a theoretical point of view is a space-time finite element approach using discontinuous approximation in time; • it is not yet fully understood how to use mixed finite elements in combination with the SUPG method (for recent advances see [26]).To overcome these disadvantages, alternative stabilization techniques have been developed such as the projection method proposed by Codina [17] and Codina and Blasco [18], the subgrid viscosity method or local projection method proposed by Guermond [28] and Becker and Braack [1], the polynomial pressure projection method by Dohrmann and Bochev reported in [19], and the pressure-Poisson stabilization of the Stokes equations proposed by Bochev and Gunzburger in [2].
Recently, the continuous interior penalty method of Douglas and Dupont [20] was revived as an alternative.The idea is to add a least squares penalization on the gradient jump between neighboring elements as a unified treatment of all the abovementioned instabilities.It was shown in [13,15] that the method stabilizes both instabilities due to dominating convection and instabilities due to the velocity/pressure coupling.Moreover, it was shown in [10] how this method provides a natural link between conforming and nonconforming stabilized finite element methods.It was used in [14] to provide a Reynolds number independent stabilized formulation for the classical nonconforming P 1 Crouzeix-Raviart approximation for the velocities combined with elementwise constant pressures.
In this paper we extend the face oriented stabilization method to Oseen's equations, using equal order interpolation for velocities and pressure.For the case of similar stabilization strategies for element pairs satisfying the inf-sup condition we refer the reader to [11].We follow the framework proposed in [10] using weakly imposed boundary conditions as introduced by Nitsche (see [36,25]).Although the constants of our analysis inevitably depend on the parameters of the problem (since the solution depends on the physical parameters), the stabilization terms allow us to trade the need of coercivity in the H 1 -norm for coercivity in the weaker L 2 -norm plus the stabilization term, which vanishes at optimal rate under refinement.To exploit this in the analysis, we add a zero order term to Oseen's equations.With this additional term we obtain estimates that do not explode as the viscosity goes to zero, provided the exact solution is sufficiently regular.
With the proposed method, all the above-mentioned inconveniences of the SUPG method are alleviated.The formulation allows for general unstructured meshes and variable polynomial degree.The main inconveniences of the present method, however, are as follows: • Added couplings in the Jacobian matrix: the bandwidth of the system matrix doubles in two space dimensions and triples in three space dimensions.This may increase the computational cost of the linear system solution, for instance, if an incomplete LU factorization is used as preconditioner.
• The method requires stabilization terms to be evaluated on the faces of the elements and hence a table of nearest neighbors for computation of the jumps.Such features are present typically when using adaptivity with a posteriori error estimation or for discontinuous Galerkin (DG) methods.However, all stabilization terms have the same structure, allowing for one computation of one global stabilization matrix based on the gradient jumps of one component.The parameter values that change from time-step to time-step may then be updated using locally averaged weights.When time-stepping the Navier-Stokes equations this means that the stabilization matrix has to be constructed only once and for one component.This is in stark contrast to the SUPG method, where the stabilization terms have to be reconstructed at every time-step for consistency.
The outline of the paper is as follows: In the next section we introduce our model problem, Oseen's equations, and formulate the interior penalty finite element method.In section 3 we discuss the question of stability, we prove a lemma of fundamental importance for the stability of the method, and we show that the discrete problem has a unique solution.We then proceed and prove (quasi-) optimal a priori error estimates in section 4 with special focus on how to make the estimates independent of the local Reynolds number.Finally, in section 5, we study the performance of the numerical scheme on some linear model cases in three space dimensions.We make some concluding remarks in section 6 and some outlooks to future developments, with special emphasis on the relation between the present method and variational multiscale methods (VMS) in large eddy simulations (LES).

2.
A finite element method for Oseen's equations.Let Ω be a Lipschitzcontinuous domain in R d (d = 2 or 3) with a polyhedral boundary ∂Ω and outward pointing normal n.We will consider the Sobolev spaces W m,q (Ω), with norm • m,q,Ω , m ≥ 0, and q ≥ 1.In particular, we have L q (Ω) = W 0,q (Ω).We use the standard notation H m (Ω) def = W m,2 (Ω).The norm of H m (Ω) is denoted by • m,Ω and its seminorm by | • | m,Ω .The space of L 2 (Ω) divergence free functions is denoted by H 0 (div; Ω).The scalar product in L 2 (Ω) is denoted by (•, •) and its norm by • 0,Ω .The closed subspaces H 1 0 (Ω), consisting of functions in H 1 (Ω) with zero trace on ∂Ω, and L 2 0 (Ω), consisting of functions in L 2 (Ω) with zero mean in Ω, will also be used.Oseen's equations take the form where The well-posedness of this problem follows by the Lax-Milgram lemma applied in the space [H 1 0 (Ω)] d ∩ H 0 (div; Ω) (see, for instance, [27]).Let {T h } 0<h≤1 denote a family of triangulations of the domain Ω without hanging nodes.For each triangulation T h , the subscript h ∈ (0, 1] refers to the level of refinement of the triangulation, which is defined by with h K the diameter of K. We also define the elementwise constant function h| K = h K .The interior of a triangle K will be denoted by • K, and N (K) will stand for the set of elements sharing at least one node with the element K.Moreover, we will assume that the family {T h } 0<h≤1 has the following regularity properties: 1. Local shape regularity: for all K ∈ T h with h ∈ (0, 1] there holds where ρ K stands for the diameter of the largest ball contained in K, and c 0 is a fixed positive constant.2. Local quasi-uniformity: for all K ∈ T h with h ∈ (0, 1] there holds where ρ > 1 is a given parameter depending on the local uniformity of {T h } 0<h≤1 .We will also assume that the data are sufficiently well resolved in the sense that there exists ρ β > 1 such that Note that this is a hypothesis on the mesh and not on the data.Under the assumption that β ∈ W 1,∞ (N (K)) this can be ensured by for some constant c β > 0 small enough.
For the error analysis, we shall use the trace inequality where C T is a generic constant independent of h K (for a proof, see [37, p. 26]).
For a given piecewise continuous function ϕ, the jump [[ϕ]] over a face f is defined by where n f is a normal unit vector on f and x ∈ f .
In this paper we let V k h denote the standard space of continuous functions of piecewise polynomial order k ≥ 1, and H 2 (T h ) the space of piecewise H 2 functions For the velocities we will use the space [V k h ] d and for the pressure we will use In what follows, we let π h,k , i h,k , and C h,k denote (respectively) the L 2projection operator, the nodal interpolation operator, and the Clément interpolant onto the finite element spaces, and we make no notational difference between the projection onto the velocity and pressure spaces.We also introduce a piecewise linear approximated velocity Here and in the following C denotes a constant independent of the problem parameters and the local mesh size, but not necessarily of the local mesh geometry.Denoting the product space

and with
n the outward pointing normal to ∂Ω, and using the notation To keep down notation we have used a canonical stabilization parameter γ for all terms.In practice this term, however, can be chosen distinctly for different terms.The gradient jump terms serve three purposes: 1. stabilization of the convective terms (the first sum in (2.13)); 2. giving additional control of the incompressibility condition (the second sum in (2.13)); and 3. making the discretization inf-sup stable (the sum in (2.14)).We will see in the analysis that these three objectives are all obtained in the same fashion and that essentially the gradient jump operator can stabilize any instability provoked by a first order term.
Assuming sufficient regularity of the exact solution the above formulation is strongly consistent.More generally we have the following result.
Lemma 2.1 (modified Galerkin orthogonality).Assume that (u, p), the solution of (2.1), belongs to the space [H 3/2+ε (Ω)] d ×L 2 0 (Ω), with ε > 0, and let (u h , p h ) ∈ W k h be the solution of (2.10).Then Proof.This is an immediate consequence of the consistency of the standard Galerkin method and the fact that, under the regularity assumptions, j u (u, v h ) = 0 since [[∇u]] f = 0 for all interior faces f .

Stability of the method.
Stability in the face oriented stabilization method is based on the following lemma, which was proved for piecewise linear continuous approximation in [10] (for a similar result with applications to DG methods see [33]).
Here we extend this result to arbitrary polynomial degree.Note that we give a lower bound as well.This is not needed for the analysis but shows that in some sense the stabilizing terms are optimal.
Lemma 3.1.There exist an interpolation operator d and constants γ, γ low depending on the local mesh geometry and the polynomial degree, but not on the local mesh size, such that Proof.First note that, as pointed out in [10], For each node x i , let n i be the number of elements containing x i as a node.Then we define a quasi-interpolant π * h,k of degree k by For each element K ∈ T h consider the function Clearly, δ K (x j ) = 0 for each interior node x j ∈ • K, whereas on the element faces, i.e., for all nodes x j ∈ ∂K, we have where P (K, K ) stands for the set of faces between K and K (the shortest path).We now introduce the reference element K and, for each K ∈ T h , the affine mapping Therefore, by equivalence of norms on discrete spaces, using a standard scaling argument (see [27, p. 96]) and (3.1), it follows that where, in the two last inequalities, E(K) denotes the set of faces containing some node of K. On the other hand, the local quasi-regularity of T h implies that the maximum number of occurrences of a face in all the sets E(K) is bounded by a fixed constant independent of h K .Then, by summation on K, we get the upper bound The lower bound follows by considering the L 2 -norm of the discontinuous function δ over the reference patch Ĝ consisting of the reference element K and its nearest neighbors.Clearly if δ Ĝ = 0, then Hence by norm equivalence on discrete spaces we have The claim then follows in the same fashion as the first part of the proof by scaling and extension to all of T h .
Using the same technique we immediately have the following corollary where for simplicity the lower bounds are omitted.
Corollary 3.2.Under the same assumptions as Lemma 3.1 we have, with α > 0 and φ some function that is constant per element, for all (v h , q h ) ∈ W k h and with γ > 0 constants independent of h.We now introduce the following mesh-dependent norm for the velocity: The following lemma gives the coercivity of our discrete operator with respect to this mesh-dependent norm.

Lemma 3.3 (coercivity).
There exists a constant C > 0, depending only on Ω and γ, such that where we used the fact that, after integration by parts and since ∇ • β = 0, The last term in (3.4) can be bounded using the Cauchy-Schwarz inequality followed by a trace inequality, to obtain In what follows we will assume that γ > 4C T > 0, (3.5) and therefore From (3.4), we then get and consequently In particular, since 0 < h ≤ 1 and by choosing (accordingly with (3.5)) with ε > 0 sufficiently small, one obtains We then conclude the proof using Korn's inequality (see [6]).
In what follows, we shall make use of the following discrete pressure and velocity subspaces: In addition, Q k h \C 1 h,k will stand for the supplementary of The following lemma ensures, in particular, that V div h,k is not trivial.Lemma 3.4.There exists a constant β > 0, independent of h, such that Thus, using integration by parts and (2.12), we have Thus, using the orthogonality of the L 2 -projection, we have Thus, from (3.7), it follows that In addition, using H 1 -stability of the L 2 -projection (see [5]) and (3.6), we have which completes the proof.We now state the well-posedness of the discrete problem.Theorem 3.5.The discrete problem (2.10) has a unique solution.
Proof.Problem (2.10) can be written, in operator form, as with We also introduce the operator in other words, From Lemma 3.4, it follows that B 1 is surjective and (B 1 ) T is injective (see [27, p. 58]).We then deduce that V div h,k def = Ker(B 1 ) = {0}.Let us consider the following reduced formulation (derived from (2.10) with Since, by construction, C 1 h,k = Ker(J), we conclude that J is invertible in By plugging this expression into the first equation of (3.9), we obtain that .
Existence and uniqueness of u h follow by the positivity of A (Lemma 3.3) and the nonnegativity of B. We may then recover ph uniquely from (3.10).Therefore, the reduced problem (3.9) has a unique solution.On the other hand, from the first equation of (3.9), it follows that with Ker(B 1 ) 0 standing for the polar set of Ker(B 1 ).From Lemma 3.4, it follows that B 1 is an isomorphism from C 1 h,k onto Ker(B 1 ) 0 (see [27, p. 58]).Thus, there exists a unique Therefore, from (3.11) and (3.9), and by noticing that (B 1 ) T p 1 = B T p 1 and Jp 1 = 0, it follows that problem (2.10) has a unique solution, given by (u h , p h def = ph −p 1 ).

4.
Convergence of the method.The parameter for the pressure stabilization scales as h 2 K / β 0,∞,K when the local Reynolds number Re K is big, and as h 3 K /ν when Re K is small.The stabilizing terms acting on the velocity scale as β 0,∞,K h 2 K at a high local Reynolds number and as Re K β 0,∞,K h 2 K for a low Reynolds number.The factor Re K β 0,∞,K h 2 K in the velocity stabilization may be omitted in the low Reynolds regime without perturbing the convergence.We will now show that this scaling gives optimal a priori error estimates in the high (local) Reynolds number regime when the solution is smooth, (u, p) ∈ [H k+1 (Ω)] d+1 , and in the low (local) Reynolds number regime under standard regularity assumptions.We then prove, using the Aubin-Nitsche duality technique (see, e.g., [21]), that the velocities have optimal convergence order also in the L 2 -norm, when the local Reynolds number is low, without any modification of the stabilization.
First, we summarize some stability properties of the L 2 -projection with weighted norms and show an approximability result for the triple norm (3.3).
Remark 4.1.In the remainder of this section, C > 0 stands for a generic constant independent of h and the physical parameters.
To prove approximability for the L 2 -projection on locally quasi-uniform meshes we need some additional stability for the L 2 -projection from [3] that we state here without proof.
Lemma 4.2.For ρ, η > 0 sufficiently small and for all φ A direct consequence of this result is stated in the following corollary.Corollary 4.3.Under the assumptions of the previous lemma, we have for all u ∈ H r (Ω), with r ≥ 1, r u def = min{r, k + 1}, 0 ≤ l ≤ r u , α ∈ N d , and ∂ α the standard multi-index notation for high order derivatives.
In order to obtain localized estimates we now show that the weights appearing in our stabilization allow for L 2 -stability.
Proof.We give the proof only for φ 1 ; the argument for the rest is similar.First, note that for all K ∈ T h , We now distinguish two cases.On one hand, if Re K ≤ 1, we have φ 1|K = ν − 1 2 .Thus, from (4.1) and since ρ β ρ > 1, it follows that On the other hand, if Re K > 1, we get φ 1|K = β K .Therefore, from (4.1) and the assumptions on the mesh (2.5) and (2.6), we have The lower bound follows in a similar fashion.Finally, for the derivative, using the bounds on φ * 1 and the regularity of the mesh (2.4), and since ρ β ρ > 1, we obtain which completes the proof.Remark 4.5.It follows from Lemma 4.4 that for the weight functions φ * i , 1 ≤ i ≤ 5, the stability estimate of Lemma 4.2 holds, provided ρ β and ρ are sufficiently close to 1. From now on we assume that this is the case.
The following lemma states the approximation properties of the L 2 -projection in the triple norm ||| • |||.
Lemma 4.6 (velocity approximability).Assume that ρ β and ρ are sufficiently close to 1.Then, there holds We give the proof for the first term only.The argument for the second term is similar.By the stability estimate for the L 2 -projection Lemma 4.2 we have Using now the H 1 -stability of the L 2 -projection on locally quasi-uniform meshes (see [5]) we get We treat the boundary terms using the trace inequality (2.8) in combination with Lemma 4.2 and approximation, which yields The interior penalty terms are treated in the same fashion as the boundary terms.We have The first term in this inequality can be estimated using that ξ(Re K ) ≤ 1, the trace inequality (2.8), an inverse inequality, and the H 1 -stability of the L 2 -projection (see [5]), which yields and so the proof is finished.
For the pressure, we have the following result.Lemma 4.7 (pressure approximability).Under the assumptions of Lemma 4.6, there holds for all p ∈ H s (Ω) with s ≥ 1 and s p def = min{k + 1, s}.Proof.As p may be only H 1 (Ω), we must replace the nodal interpolant by the Clément interpolant in the analysis.The proof for the first term is similar to (4.2), by replacing i h,k by C h,k .The estimate for the second term follows from Corollary 4.3.Finally, for the interior penalty term, since [[C h,k ∇p]] = 0 and using a trace inequality followed by an inverse inequality, we have where we concluded using the stability lemma, Lemma 4.2, with weight function φ * 4 , and the optimal approximation properties of the Clément interpolant (see [16,21]).

Energy norm error estimate.
In this section we prove convergence in the triple norm.These results are optimal independently of the local Reynolds number when the exact solution is sufficiently smooth.
We start by proving a technical lemma.Lemma 4.8.For all Proof.Let A 1 denote the set of elements K ∈ T h such that ξ(Re K ) ≥ 1, and A 2 the set of elements such that ξ(Re K ) < 1.It then follows that |β| ∞,0,K h K < ν for K ∈ A 2 , and we may write where the last inequality follows by a trace inequality, an inverse inequality in the second term, and extending the sums over all T h .The second inequality follows in a similar fashion, noting that and so the proof is completed.
The main result of this paragraph is stated in the following theorem.Theorem 4.9.Assume (u, p) ∈ [H r (Ω)] d × H s (Ω), with r ≥ 2 and s ≥ 1, is the solution of (2.1) and (u h , p h ) ∈ W k h is the solution of (2.10).Then, under the assumptions of Lemma 4.6, there holds +C max , with r u = min{k + 1, r} and s p = min{k + 1, s}.
Proof.Let us decompose the error u − u h in two parts: We also consider the discrete pressure error For the second boundary term we use the Cauchy-Schwarz inequality followed by a trace inequality and an approximation argument, similar to (4.3)-(4.4), to obtain The convective term is controlled using a local inverse inequality, Lemma 4.4, Corollary 4.3, and the orthogonality of the L 2 -projection, after having replaced the continuous velocity field β by its piecewise linear interpolant β h , where we used Corollary 4.3 and Lemma 4.8 in the last inequality.
Collecting terms we have For the second term in (4.6), using the orthogonality of the L 2 -projection, Lemmas 4.7 and 4.8, and replacing u with p in (4.2), we have In a similar fashion, after integration by parts in the third term, one obtains Finally, using Lemma 4.7, for the interior penalty terms we have .
We conclude the proof by collecting the results of (4.9)-(4.12) in (4.6) and applying the approximation lemma, Lemma 4.6.
The following corollary follows from (4.6) in combination with (4.5) and Lemma 4.7.
Remark 4.11.In a physical situation the velocity gradient on boundaries with noslip conditions is known to scale as |β| 1,∞,∂Ω ∼ ν − 1 2 .If in the boundary layer h K ∼ ν, that is, a low local Reynolds number on the boundary, then the estimate is dominated by the H 1 (Ω) contribution from the boundary that converges at optimal rate since the layer is resolved.The condition (2.7) is satisfied with c β ∼ ν 1 2 showing that the strongest constraint on the mesh is not that of (2.7), but that of the h K ν contribution on the boundary.In laminar free-flow we can expect |β| 1,∞,K ≤ c β 0,∞,K to hold, and hence the convergence in the L 2 -norm in this regime is of the quasi-optimal rate h k+ 1 2 for a sufficiently regular solution.
4.2.Recovering the pressure.In this section, we provide an estimate of the L 2 -norm of the pressure error.This is the aim of the following theorem, which ensures that the pressure converges at the rate of the velocity.
Theorem 4.12.Assume (u, p) ∈ [H r (Ω)] d × H s (Ω), with r ≥ 2 and s ≥ 1, is the solution of (2.1) and (u h , p h ) ∈ W k h is the solution of (2.10).Then, under the assumptions of Lemma 4.6, there holds with C u the convergence rate of |||u − u h ||| given by Theorem 4.9, and C L a positive constant depending only on Ω.
Proof.Following [27, Corollary 2.4], there exists with C L > 0 a constant, depending on Ω, which scales as a distance.Thus, using the modified Galerkin orthogonality (Lemma 2.1), we readily obtain Thus, after integrating by parts, we get For the first term, using the orthogonality of the L 2 -projection, the Cauchy-Schwarz inequality, Corollary 3.2, (4.13), and Corollary 4.3, we get Using the definition (2.11) of the bilinear form a h , and after integration by parts in the convective term, we have For the convective term we have, using the H 1 -stability of the L 2 -projection (see [5]) and (4.13), The boundary terms are controlled in the following fashion: In addition, as in (4.7) and (4.8), we have 2 (e π ) 0,∂Ω + (ν h) In the same fashion, we obtain Finally, from (4.13), it follows that  Theorem 4.14.Assume that the solution (u, p) of (2.1) belongs to [H 2 (Ω)] d × H 1 (Ω) and let (u h , p h ) ∈ W k h be the solution of (2.10).Assume also that and that the solution (ϕ, ψ) of the adjoint problem ] and satisfies Then, there holds u − u h 0,Ω ≤ ch 2 u 2,Ω + p 1,Ω , with constant c > 0 independent of h, but depending on the physical parameters.
Proof.Multiplying the first equation of (4.23) by u − u h and the second by −(p−p h ), integrating by parts, and using the modified Galerkin orthogonality (Lemma 2.1), it follows that Following the argument of the proofs of Theorems 4.9 and 4.12, and using Lemma 4.6 and (4.22), we get Using Cauchy-Schwarz, Lemma 4.6, and (4.22), one obtains Finally, from Lemma 4.7 and (4.22), for the last term we have The proof concludes by combining the above estimations with Theorems 4.9 and 4.12, Corollary 4.10, (4.22), and the assumed regularizing behavior (4.24).Let us sum up the results provided by Theorems 4.9 and 4.12.When the local Reynolds number is high and the solution is regular, we enjoy an optimal O(h k+ 1 2 ) convergence order of the error in the L 2 -norm for the velocity and the pressure.For less regular solutions, for instance, when the pressure is in H 1 (Ω) and the velocity is in [H 2 (Ω)] d , we get an optimal O(h) estimate in the energy norm, when the local Reynolds number is low, but a suboptimal estimate of O(h 1 2 ) when the local Reynolds number is high.This is due to the fact that the inconsistencies in the pressure stabilization pollute the energy norm estimate for the velocities.
Remark 4.15.Note that by adding the L 2 -coercivity, we can use the stabilization term to control the convective term without using the H 1 -coercivity; this leads to a quasi-optimal estimate in the weaker L 2 -norm, with a ν-weighted H 1 contribution showing that the stabilization handles the numerical instability induced by treating nonsymmetric terms using the standard Galerkin method.In case σ = 0 the H 1 estimate obtained by a standard energy argument will scale as ν − 1 2 , reflecting the physical instability of the problem.

Numerical results.
In this section we report several numerical experiments that show the good convergence properties of our stabilized finite element method.In particular, we recover the convergence rates obtained in section 4. We consider problem (2.1) in three dimensions with nonhomogeneous boundary conditions.The right-hand side f and the boundary data are chosen in order to ensure that the exact solution of (2.1) is given by the following expression [22]: The resulting continuous problem was solved approximately using the stabilized discrete formulation (2.10); however, the boundary conditions were strongly enforced.All numerical tests have been performed using conforming linear and quadratic finite elements for velocity and pressure, namely, P 1 /P 1 and P 2 /P 2 (implemented in a threedimensional research code [23]).The stabilization parameter involved in the jumps terms (2.13) and (2.14) were chosen as In Figure 5.1 we show, respectively, the velocity and pressure convergence histories for k = 1 and k = 2.Note that, in both cases, the numerical solution exhibits optimal convergence order and is hence in agreement with Theorems 4.9 and 4.12.We show in Figure 5.2 the pressure contours in two different meshes (which are depicted in Figure 5.3) using linear elements.No spurious pressure oscillations are observed.We report in Figure 5.4 the contours of the second component of u h , u h2 , in the left plot with full stabilization and in the right plot setting the stabilization parameter for the term associated with the streamline derivative to zero, on the cutting plane x = 0.5.Although the exact solution is smooth, the plot of the unstabilized solution (right) exhibits spurious oscillations.Note that the spurious velocity oscillations (right) are completely controlled by the streamline-derivative jumps (left).In what follows we will replace, in (5.1), the expression for the pressure by Clearly, this function satisfies p ∈ H 1 (Ω) but does not belong to H 2 (Ω). Figure 5.5  shows the pressure contours in a cut of a coarse and a fine mesh.Once more no spurious pressure oscillations are observed.Figure 5.6 shows the velocity and pressure convergence histories using linear elements.We get the suboptimal O(h 1 2 ) order for the velocity in the H 1 -norm in the case of high local Reynolds numbers, in agreement with Theorem 4.9.The L 2 -norm of the velocities, on the other hand, is still not far from the quasi-optimal O(h left graphic.In addition, as predicted in Theorem 4.14, we notice that the convergence order for the velocity in the L 2 -norm is O(h 2 ).

Conclusion and outlook.
In this paper, we have extended the results reported in [15,13] to Oseen's equations using equal order interpolation and finite element spaces of arbitrary polynomial order.The stability properties of the method are based on an interior penalty term giving L 2 -control of the jump of the gradient over interior element faces.We have shown that such a stabilization operator may be used to control all the nonsymmetric first order terms of Oseen's equations and that they give control only of the part of the operator that is not in the finite element space.In this sense the proposed method is a minimal stabilized method (see [8]).
The convergence analysis shows that the method has (quasi-) optimal convergence properties both in the L 2 -norm and in the energy norm when the solution is sufficiently regular or the local Reynolds number is low.When physically realistic regularities are considered (p ∈ H 1 (Ω)) and the local Reynolds number is high, the convergence may become suboptimal O(h 1 2 ) due to the inconsistencies in the pressure stabilization.In some numerical examples we illustrated the theoretical results.The method shows very good performance in all regimes.In particular, we observe that in the high Reynolds number regime the scheme degenerates to the theoretical O(h k+ 1 2 ) convergence in the L 2 -norm predicted by the theory only in the case where the pressure is only H 1 and where the theoretical prediction is O(h 1 2 ).The method presented here has some common features with VMS for LES as introduced in [30].However, unlike the VMS, where two scales V h and V H defined by hierarchic meshes are considered (see, e.g., [35,31,4]), in our case the finite element space V h represents the only resolved scale and the "turbulent" viscosity acts only on the gradient component that is not resolved on V h .Recently, John and Kaya [31] proposed a VMS using a projection method framework which essentially takes the form of a standard Galerkin formulation for u h supplemented with the turbulent viscosity acting only on the fine scales in the form of an additional term (ν T (I − P H ) (u h ), (I − P H ) (v h )), (6.1)where P H is some map from fine scales to coarse scales.Comparing this now with the face oriented stabilization method, we would choose H = h and thus make the turbulent viscosity act only on the scales that are not resolved on the space V h .Applying Lemma 3.1 we immediately get an interior penalty interpretation of the term (6.1), with ] ds, and we conclude that a possible subgrid modeling term would be where the choice of ν T now is a modeling issue.It should be noted that the choice ν T = γh K gives us a term which is asymptotically equivalent to the face penalty operator using the whole gradient.However, other choices of ν T based on modeling considerations are possible.
For sufficiently high polynomial degree there exists a C 1 subspace of V h with approximation properties.It follows that the solution may be decomposed into one C 1 part which is untouched by the stabilizing terms and another C 0 part which is penalized.We conclude that the method enjoys the scale separation property characteristic for VMS as proposed in [30] by polynomial order rather than by hierarchic meshes.Future work will focus on the extension of the present method to the Navier-Stokes equations both from a numerical and a theoretical standpoint.
Finally, we remark that the Nitsche-type weak boundary conditions used in this paper, while nonstandard, have the benefit of acting as slip boundary conditions in the high Reynolds number regime and as conditions when the boundary layers are resolved; this may be favorable in LES (see Layton [34]).

1 2 p
− p h 0,Ω .(4.21)We conclude the proof by collecting the estimations (4.15)-(4.20) in (4.14), using (4.21) and Theorem 4.9.Remark 4.13.Let us notice that the three terms appearing in the error estimate of the previous theorem scale with the right dimensions.

4. 3 .
Low Reynolds number optimality.The following theorem gives an optimal L 2 -error estimate for velocity when the local Reynolds number is low.