# Continuous Interior Penalty Finite Element Method for Oseen's Equations

Erik Burman, Miguel Angel Fernández, Peter Hansbo

► To cite this version:

Erik Burman, Miguel Angel Fernández, Peter Hansbo. Continuous Interior Penalty Finite Element Method for Oseen's Equations. SIAM Journal on Numerical Analysis, 2006, 44 (3), pp.1248-1274. 10.1137/040617686 . hal-01345519

HAL Id: hal-01345519

https://hal.science/hal-01345519

Submitted on 13 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CONTINUOUS INTERIOR PENALTY FINITE ELEMENT METHOD FOR OSEEN'S EQUATIONS

ERIK BURMAN†, MIGUEL A. FERNÁNDEZ‡, AND PETER HANSBO§

**Abstract.** 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 $L^2$-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.

**Key words.** finite element methods, stabilized methods, continuous interior penalty, Oseen's equations

**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

<footnote>†Institut d'Analyse et Calcul Scientifique (CMCS/IACS), Ecole Polytechnique Federale de Lausanne, Switzerland (Erik.Burman@epfl.ch).</footnote>

<footnote>‡INRIA projet REO, Rocquencourt BP 105, F-78153 Le Chesnay Cedex, France (miguel.fernandez@inria.fr).</footnote>

<footnote>§Department of Applied Mechanics, Chalmers University of Technology, SE-41296 Göteborg, Sweden (hansbo@am.chalmers.se).</footnote>
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 time-stepping 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 above-mentioned 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 $\Omega$ be a Lipschitz-continuous domain in $\mathbb{R}^d$ ($d=2$ or 3) with a polyhedral boundary $\partial\Omega$ and outward pointing normal $\mathbf{n}$. We will consider the Sobolev spaces $W^{m,q}(\Omega)$, with norm $\|\cdot\|_{m,q,\Omega}$, $m \ge 0$, and $q \ge 1$. In particular, we have $L^q(\Omega) = W^{0,q}(\Omega)$. We use the standard notation $H^m(\Omega) \stackrel{\text{def}}{=} W^{m,2}(\Omega)$. The norm of $H^m(\Omega)$ is denoted by $\|\cdot\|_{m,\Omega}$ and its seminorm by $|\cdot|_{m,\Omega}$. The space of $L^2(\Omega)$ divergence free functions is denoted by $H_0(\text{div};\Omega)$. The scalar product in $L^2(\Omega)$ is denoted by $(\cdot, \cdot)$ and its norm by $\|\cdot\|_{0,\Omega}$. The closed subspaces $H_0^1(\Omega)$, consisting of functions in $H^1(\Omega)$ with zero trace on $\partial\Omega$, and $L_0^2(\Omega)$, consisting of functions in $L^2(\Omega)$ with zero mean in $\Omega$, will also be used. Oseen's equations take the form

$$ (2.1) \qquad \left\{ \begin{array}{ll} \sigma \mathbf{u} + \beta \cdot \nabla \mathbf{u} - 2\nabla \cdot (\nu \boldsymbol{\epsilon}(\mathbf{u})) + \nabla p = \mathbf{f} & \text{in } \Omega, \\ \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega, \\ \mathbf{u} = 0 & \text{on } \partial\Omega, \end{array} \right. $$

where $\mathbf{u} \in [H_0^1(\Omega)]^d \cap H_0(\text{div};\Omega)$, $\beta \in [W^{1,\infty}(\Omega)]^d \cap H_0(\text{div};\Omega)$, $p \in L_0^2(\Omega)$, $\boldsymbol{\epsilon}(\mathbf{u})$ stands for the strain rate tensor

$$ \boldsymbol{\epsilon}(\mathbf{u}) \stackrel{\text{def}}{=} \frac{1}{2} (\nabla \mathbf{u} + (\nabla \mathbf{u})^\mathrm{T}), $$

$\mathbf{f} \in [L^2(\Omega)]^d$ is a given source term, and $\sigma, \nu$ are positive constants.

The weak formulation of problem (2.1) reads as follows: find $(\mathbf{u}, p) \in [H_0^1(\Omega)]^d \times L_0^2(\Omega)$ such that

$$ (2.2) \qquad \left\{ \begin{array}{ll} a(\mathbf{u}, \mathbf{v}) + b(p, \mathbf{v}) = (\mathbf{f}, \mathbf{v}), & \forall (\mathbf{v}, q) \in [H_0^1(\Omega)]^d \times L_0^2(\Omega), \\ b(q, \mathbf{u}) = 0 & \end{array} \right. $$
where

$$
\begin{equation} \tag{2.3}
\begin{aligned}
a(\mathbf{u}, \mathbf{v}) &\stackrel{\text{def}}{=} (\sigma \mathbf{u}, \mathbf{v}) + (\boldsymbol{\beta} \cdot \nabla \mathbf{u}, \mathbf{v}) + 2(\nu \epsilon(\mathbf{u}), \epsilon(\mathbf{v})), \\
b(p, \mathbf{v}) &\stackrel{\text{def}}{=} -(p, \nabla \cdot \mathbf{v}).
\end{aligned}
\end{equation}
$$

The well-posedness of this problem follows by the Lax–Milgram lemma applied in the
space $[H_0^1(\Omega)]^d \cap H_0(\text{div}; \Omega)$ (see, for instance, [27]).

Let {$\mathcal{T}_h$}$_{0<h\le 1}$ denote a family of triangulations of the domain $\Omega$ without hanging nodes. For each triangulation $\mathcal{T}_h$, the subscript $h \in (0,1]$ refers to the level of refinement of the triangulation, which is defined by

$$
h \stackrel{\text{def}}{=} \max_{K \in \mathcal{T}_h} h_K,
$$

with $h_K$ the diameter of $K$. We also define the elementwise constant function $\tilde{h}|_K = h_K$. The interior of a triangle $K$ will be denoted by $\tilde{K}$, and $\mathcal{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 $\{\mathcal{T}_h\}_{0<h\le 1}$ has the following regularity properties:

1. Local shape regularity: for all $K \in \mathcal{T}_h$ with $h \in (0, 1]$ there holds

$$
(2.4) \qquad \frac{h_K}{\rho_K} < c_0,
$$

where $\rho_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 \in \mathcal{T}_h$ with $h \in (0, 1]$ there holds

$$
(2.5) \qquad \frac{1}{\rho} h_{K'} \le h_K \le \rho h_{K'} \quad \forall K' \in \mathcal{N}(K),
$$

where $\rho > 1$ is a given parameter depending on the local uniformity of
$\{\mathcal{T}_h\}_{0<h\le 1}$.

We will also assume that the data are sufficiently well resolved in the sense that there
exists $\rho\beta > 1$ such that

$$
(2.6) \quad \frac{1}{\rho\beta} \|\beta\|_{0,\infty,K'} \le \|\beta\|_{0,\infty,K} \le \rho\beta \|\beta\|_{0,\infty,K'} \quad \forall K' \in \mathcal{N}(K).
$$

Note that this is a hypothesis on the mesh and not on the data. Under the assumption
that $\boldsymbol{\beta} \in W^{1,\infty}(\mathcal{N}(K))$ this can be ensured by

$$
(2.7) \qquad |\beta|_{1,\infty,\mathcal{N}(K)} \le c_{\beta} h_K^{-1} \|\beta\|_{0,\infty,\mathcal{N}(K)}
$$

for some constant $c_{\beta} > 0$ small enough.

For the error analysis, we shall use the trace inequality

$$
(2.8) \quad \|v\|_{0,\partial K}^2 \le C_T \left( h_K^{-1} \|v\|_{0,K}^2 + h_K \|v\|_{1,K}^2 \right) \quad \forall v \in H^1(K),
$$

where $C_T$ is a generic constant independent of $h_K$ (for a proof, see [37, p. 26]).

For a given piecewise continuous function $\varphi$, the jump $\llbracket \varphi \rrbracket$ over a face $f$ is defined
by

$$
\llbracket \varphi \rrbracket(\mathbf{x}) \stackrel{\text{def}}{=} \begin{cases} \lim_{t \to 0^+} (\varphi(\mathbf{x} - t\mathbf{n}_f) - \varphi(\mathbf{x} + t\mathbf{n}_f)) & \text{if } f \not\subset \partial\Omega, \\ 0 & \text{if } f \subset \partial\Omega, \end{cases}
$$

where $\mathbf{n}_f$ is a normal unit vector on $f$ and $\mathbf{x} \in f$.
In this paper we let $V_h^k$ denote the standard space of continuous functions of piecewise polynomial order $k \ge 1$,

$$V_h^k \stackrel{\text{def}}{=} \{v \in H^1(\Omega) : v|_K \in P_k(K) \quad \forall K \in \mathcal{T}_h\},$$

and $H^2(\mathcal{T}_h)$ the space of piecewise $H^2$ functions

$$H^2(\mathcal{T}_h) \stackrel{\text{def}}{=} \{v : \Omega \to \mathbb{R} : v|_K \in H^2(K) \quad \forall K \in \mathcal{T}_h\}.$$

For the velocities we will use the space $[V_h^k]^d$ and for the pressure we will use $Q_h^k = V_h^k \cap L_0^2(\Omega)$. In what follows, we let $\pi_{h,k}, \tilde{h}_{h,k}$, and $C_{h,k}$ denote (respectively) the $L^2$-projection 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 $\beta_h \in [V_h^1]^d$ such that

$$(2.9) \qquad \|\beta - \beta_h\|_{0,\infty,K} \le C h_K |\beta|_{1,\infty,K} \quad \forall K \in \mathcal{T}_h.$$

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 $W_h^k \stackrel{\text{def}}{=} [V_h^k]^d \times Q_h^k$ our finite element method reads as follows: find $(\mathbf{u}_h, p_h) \in W_h^k$ such that

$$(2.10) \qquad a_h(\mathbf{u}_h, \mathbf{v}_h) + b_h(p_h, \mathbf{v}_h) - b_h(q_h, \mathbf{u}_h) + j_\mathbf{u}(\mathbf{u}_h, \mathbf{v}_h) + j_p(p_h, q_h) = (\mathbf{f}, \mathbf{v}_h)$$

for all $(\mathbf{v}_h, q_h) \in W_h^k$, and with

$$(2.11) \qquad a_h(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} a(\mathbf{u}_h, \mathbf{v}_h) - \langle 2\nu\epsilon(\mathbf{u}_h)\mathbf{n}, \mathbf{v}_h \rangle_{\partial\Omega} - \langle \mathbf{u}_h, 2\nu\epsilon(\mathbf{v}_h)\mathbf{n} \rangle_{\partial\Omega} \\ - \langle \boldsymbol{\beta} \cdot \mathbf{n}\mathbf{u}_h, \mathbf{v}_h \rangle_{\partial\Omega_{\text{in}}} + \langle \gamma(\nu/\tilde{h})\mathbf{u}_h, \mathbf{v}_h \rangle_{\partial\Omega} \\ + \langle \gamma \max\{|\boldsymbol{\beta}|, \nu/\tilde{h}\} \mathbf{u}_h \cdot \mathbf{n}, \mathbf{v}_h \cdot \mathbf{n} \rangle_{\partial\Omega},$$

$$(2.12) \qquad b_h(p_h, \mathbf{v}_h) \stackrel{\text{def}}{=} b(p_h, \mathbf{v}_h) + \langle p_h, \mathbf{v}_h \cdot \mathbf{n} \rangle_{\partial\Omega},$$

$$(2.13) \qquad j_\mathbf{u}(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \sum_{K \in \mathcal{T}_h} \gamma\xi(\text{Re}_K)h_K^2 \int_{\partial K} \|\boldsymbol{\beta} \cdot \mathbf{n}\|_{0,\infty,\partial K} \llbracket \mathbf{n} \cdot \nabla \mathbf{u}_h \rrbracket \cdot \llbracket \mathbf{n} \cdot \nabla \mathbf{v}_h \rrbracket \, \mathrm{d}s \\ + \sum_{K \in \mathcal{T}_h} \gamma\xi(\text{Re}_K)\|\boldsymbol{\beta}\|_{0,\infty,K}h_K^2 \int_{\partial K} \llbracket \nabla \cdot \mathbf{u}_h \rrbracket \cdot \llbracket \nabla \cdot \mathbf{v}_h \rrbracket \, \mathrm{d}s,$$

$$(2.14) \qquad j_p(p_h, q_h) \stackrel{\text{def}}{=} \sum_{K \in \mathcal{T}_h} \gamma\xi(\text{Re}_K) \frac{h_K^2}{\|\boldsymbol{\beta}\|_{0,\infty,K}} \int_{\partial K} \llbracket \nabla p_h \rrbracket \cdot \llbracket \nabla q_h \rrbracket \, \mathrm{d}s,$$

$\mathbf{n}$ the outward pointing normal to $\partial\Omega$, and using the notation

$$\text{Re}_K \stackrel{\text{def}}{=} \frac{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K}{\nu}, \quad \xi(\lambda) \stackrel{\text{def}}{=} \min\{1, \lambda\}, \quad \partial\Omega_{\text{in}} \stackrel{\text{def}}{=} \{\mathbf{x} \in \partial\Omega : (\boldsymbol{\beta} \cdot \mathbf{n})(\mathbf{x}) < 0\}, \\ \langle x, y \rangle_{\partial\Omega} \stackrel{\text{def}}{=} \sum_{K \in \mathcal{T}_h} \int_{\partial K \cap \partial\Omega} xy \, \mathrm{d}s, \quad \tilde{h} \in H^2(\mathcal{T}_h) \quad \text{with} \quad \tilde{h}|_K \stackrel{\text{def}}{=} h_K.$$

To keep down notation we have used a canonical stabilization parameter $\gamma$ 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 $(\mathbf{u}, p)$, the solution of (2.1), belongs to the space $[H^{3/2+\varepsilon}(\Omega)]^d \times L_0^2(\Omega)$, with $\varepsilon > 0$, and let $(\mathbf{u}_h, p_h) \in W_h^k$ be the solution of (2.10). Then

$$a_h(\mathbf{u} - \mathbf{u}_h, \mathbf{v}_h) + b_h(p - p_h, \mathbf{v}_h) - b_h(q_h, \mathbf{u} - \mathbf{u}_h) + j_{\mathbf{u}}(\mathbf{u} - \mathbf{u}_h, \mathbf{v}_h) - j_p(p_h, q_h) = 0$$

for all $(\mathbf{v}_h, q_h) \in W_h^k$.

*Proof.* This is an immediate consequence of the consistency of the standard Galerkin method and the fact that, under the regularity assumptions, $j_{\mathbf{u}}(\mathbf{u}, \mathbf{v}_h) = 0$ since $\|\nabla \mathbf{u}\|_f = 0$ for all interior faces $f$. $\square$

**3. 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 $\pi_{h,k}^* : [H^2(\mathcal{T}_h)]^d \to [V_h^k]^d$ and constants $\gamma$, $\gamma_{low}$ depending on the local mesh geometry and the polynomial degree, but not on the local mesh size, such that

$$\gamma_{low} j\beta(\mathbf{v}_h, \mathbf{v}_h) \le \|h^{\frac{1}{2}}(\beta_h \cdot \nabla \mathbf{v}_h - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{v}_h))\|_{0,\Omega}^2 \le j\beta(\mathbf{v}_h, \mathbf{v}_h)$$

for all $\mathbf{v}_h \in [V_h^k]^d$, where

$$j\beta(\mathbf{v}_h, \mathbf{v}_h) = \gamma \sum_{K \in \mathcal{T}_h} \int_{\partial K} h_K^2 |\beta_h \cdot \mathbf{n}|^2 \|\nabla \mathbf{v}_h \mathbf{n}\|^2 \mathrm{d}s.$$

*Proof.* First note that, as pointed out in [10], $\|\beta_h \cdot \nabla \mathbf{u}_h\| = \beta_h \cdot \mathbf{n} \|\mathbf{n} \cdot \nabla \mathbf{u}_h\|$. 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 $\pi_{h,k}^*$ of degree $k$ by

$$\pi_{h,k}^* \mathbf{v}(x_i) \stackrel{\text{def}}{=} \frac{1}{n_i} \sum_{\{K: x_i \in K\}} \mathbf{v}|_K(x_i) \quad \forall \mathbf{v} \in [H^2(\mathcal{T}_h)]^d.$$

For each element $K \in \mathcal{T}_h$ consider the function

$$\delta_K \stackrel{\text{def}}{=} h_K^{\frac{1}{2}} (\beta_h \cdot \nabla \mathbf{v}_{h|K} - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{v}_{h|K})) .$$
Clearly, $\delta_K(x_j) = 0$ for each interior node $x_j \in \hat{K}$, whereas on the element faces, i.e., for all nodes $x_j \in \partial K$, we have

$$ (3.1) \qquad \delta_K(x_j) = h_K^{\frac{1}{2}} \frac{1}{n_j} \sum_{\{K' : x_j \in K'\}} \beta_h \cdot (\nabla \mathbf{v}_{h|K}(x_j) - \nabla \mathbf{v}_{h|K'}(x_j)) \\ = h_K^{\frac{1}{2}} \frac{1}{n_j} \sum_{\{K' : x_j \in K'\}} \sum_{f \in P(K, K')} \beta_h(x_j) \cdot [\nabla \mathbf{v}_h]_f(x_j), $$

where $P(K, K')$ stands for the set of faces between $K$ and $K'$ (the shortest path). We now introduce the reference element $\hat{K}$ and, for each $K \in \mathcal{T}_h$, the affine mapping

$$ F_K(\hat{\mathbf{x}}) = B_K \hat{\mathbf{x}} + \mathbf{b}_K \quad \forall \hat{\mathbf{x}} \in \hat{K}, $$

such that $F_K(\hat{K}) = K$. Finally, let $\varphi_j^K$ for $j = 1, \dots, k$ be the basis functions on $K$. Since $\delta_K(x_j) = 0$ for each interior node $x_j \in \hat{K}$, $\|\delta_K \circ F_K\|_{0, \partial \hat{K}}^2 = 0$ implies that $\delta_K \circ F_K = 0$ in $\hat{K}$. Therefore, by equivalence of norms on discrete spaces, using a standard scaling argument (see [27, p. 96]) and (3.1), it follows that

$$ \begin{aligned} \|\delta_K\|_{0,K}^2 &= \det B_K \|\delta_K \circ F_K\|_{0,\hat{K}}^2 \\ &\le C \det B_K \|\delta_K \circ F_K\|_{0,\partial\hat{K}}^2 \\ &= \int_{\partial\hat{K}} \frac{1}{|B_K^{-\mathrm{T}}\hat{\mathbf{n}}|} |\delta_K \circ F_K|^2 \underbrace{\det B_K |B_K^{-\mathrm{T}}\hat{\mathbf{n}}| \mathrm{d}\hat{s}}_{\mathrm{ds}} \\ &\le C |B_K^{\mathrm{T}}| \int_{\partial K} |\delta_K|^2 \mathrm{d}s \\ &\le C h_K \int_{\partial K} |\delta_K|^2 \mathrm{d}s \\ &\le C h_K \int_{\partial K} \sum_{j=1}^k |\delta_K(x_j)|^2 (\varphi_j^K)^2 \mathrm{d}s \\ &\le C h_K^2 \int_{\partial K} \sum_{j=1}^k \frac{1}{n_j} \sum_{\{K' : x_j \in K'\} \, e \in P(K, K')} \sum_{f \in E(K)} |\beta_h(x_j) \cdot [\nabla \mathbf{v}_h]_f(x_j)|^2 (\varphi_j^K)^2 \mathrm{d}s \\ &\le C h_K^2 \int_{\partial K} \sum_{j=1}^k \frac{1}{n_j} \sum_{f \in E(K)} |\beta_h(x_j) \cdot [\nabla \mathbf{v}_h]_f(x_j)|^2 (\varphi_j^K)^2 \mathrm{d}s \\ &\le C h_K^2 \sum_{f \in E(K)} \int_f |\beta_h \cdot [\nabla \mathbf{v}_h]_f|^2 \mathrm{d}s, \end{aligned} $$

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 $\mathcal{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

$$ \begin{aligned} \|h^{\frac{1}{2}} (\beta_h \cdot \nabla \mathbf{v}_h - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{v}_h))\|_{0,\Omega}^2 &\le C \sum_{K \in \mathcal{T}_h} h_K^2 \sum_{f \in E(K)} \int_f |\beta_h \cdot [\nabla \mathbf{v}_h]_f|^2 \mathrm{d}s, \\ &\le C \sum_{K \in \mathcal{T}_h} \int_{\partial K} h_K^2 |\beta_h \cdot [\nabla \mathbf{v}_h]_f|^2 \mathrm{d}s. \end{aligned} $$
The lower bound follows by considering the $L^2$-norm of the discontinuous function $\delta$ over the reference patch $\tilde{G}$ consisting of the reference element $\tilde{K}$ and its nearest neighbors. Clearly if $\|\delta\|_{\tilde{G}} = 0$, then $\beta_h \cdot \nabla v_h = \pi_{h,k}^* \beta_h \cdot \nabla v_h$ in $\tilde{G}$. This means that $\beta_h \cdot \nabla v_h$ is continuous in $\tilde{G}$ and hence $\sum_{f \in E(K)} \int_f h_K [\beta_h \cdot \nabla v_h]^2 ds = 0$. Hence by norm equivalence on discrete spaces we have

$$ \sum_{f \in E(K)} \int_f [\beta_h \cdot \nabla v_h]^2 ds \le \|\delta_G\|_{0,G}^2. $$

The claim then follows in the same fashion as the first part of the proof by scaling and extension to all of $T_h$. $\square$

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 $\alpha > 0$ and $\phi$ some function that is constant per element,*

$$ \begin{gathered} \|\phi^{\frac{1}{2}}(\nabla \cdot \mathbf{v}_h - \pi_{h,k}^*(\nabla \cdot \mathbf{v}_h))\|_{0,\Omega}^2 \le \gamma \sum_{K \in \mathcal{T}_h} \int_{\partial K} \phi h_K [\nabla \cdot \mathbf{v}_h]^2 ds, \\ (3.2) \qquad \|\phi^{\frac{1}{2}}(\nabla q_h - \pi_{h,k}^*(\nabla q_h))\|_{0,\Omega}^2 \le \gamma \sum_{K \in \mathcal{T}_h} \int_{\partial K} \phi h_K [\nabla q_h]^2 ds \end{gathered} $$

for all $(\mathbf{v}_h, q_h) \in W_h^k$ and with $\gamma > 0$ constants independent of $h$.

We now introduce the following mesh-dependent norm for the velocity:

$$ \begin{aligned} \|\mathbf{v}_h\|^2 = \sup_{h \in \mathcal{T}_h} & \|\sigma^{\frac{1}{2}} \mathbf{v}_h\|_{0,\Omega}^2 + \|\nu^{\frac{1}{2}} \nabla \mathbf{v}_h\|_{0,\Omega}^2 + j_u(\mathbf{v}_h, \mathbf{v}_h) + \|\|\beta \cdot \mathbf{n}\|^{\frac{1}{2}} \mathbf{v}_h\|_{0,\partial\Omega}^2 \\ & + \|\gamma^{\frac{1}{2}}(\nu/\tilde{h})^{\frac{1}{2}} \mathbf{v}_h\|_{0,\partial\Omega}^2 + \|\gamma^{\frac{1}{2}} \max\{|\beta|, \nu/\tilde{h}\}^{\frac{1}{2}} \mathbf{v}_h \cdot \mathbf{n}\|_{0,\partial\Omega}^2 \end{aligned} $$

for all $\mathbf{v}_h \in [V_h^k]^d$.

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 $\Omega$ and $\gamma$, such that*

$$ a_h(\mathbf{v}_h, \mathbf{v}_h) + j_u(\mathbf{v}_h, \mathbf{v}_h) \ge C \|\mathbf{v}_h\|^2 $$

for all $\mathbf{v}_h \in [V_h^k]^d$.

*Proof.* From (2.10) we get

$$ \begin{aligned} a_h(\mathbf{v}_h, \mathbf{v}_h) + j_u(\mathbf{v}_h, \mathbf{v}_h) \ge & \|\sigma^{\frac{1}{2}} \mathbf{v}_h\|_{0,\Omega}^2 + 2\|\nu^{\frac{1}{2}} \epsilon(\mathbf{v}_h)\|_{0,\Omega}^2 + j_u(\mathbf{v}_h, \mathbf{v}_h) \\ & + \frac{1}{2}\|\|\beta \cdot \mathbf{n}\|^{\frac{1}{2}} \mathbf{v}_h\|_{0,\partial\Omega}^2 + \|\gamma^{\frac{1}{2}}(\nu/\tilde{h})^{\frac{1}{2}} \mathbf{v}_h\|_{0,\partial\Omega}^2 \\ & + \|\gamma^{\frac{1}{2}} \max\{|\beta|, \nu/\tilde{h}\}^{\frac{1}{2}} \mathbf{v}_h \cdot \mathbf{n}\|_{0,\partial\Omega}^2 \\ & - \langle 4\nu\epsilon(\mathbf{v}_h)\mathbf{n}, \mathbf{v}_h \rangle_{\partial\Omega}, \end{aligned} $$

where we used the fact that, after integration by parts and since $\nabla \cdot \beta = 0$,

$$ (\beta \cdot \nabla \mathbf{v}_h, \mathbf{v}_h) = \frac{1}{2} \langle \beta \cdot \mathbf{n} \mathbf{v}_h, \mathbf{v}_h \rangle_{\partial \Omega}. $$
The last term in (3.4) can be bounded using the Cauchy–Schwarz inequality followed by a trace inequality, to obtain

$$|\langle 4\nu\epsilon(\mathbf{v}_h)\mathbf{n}, \mathbf{v}_h \rangle_{\partial\Omega}| \leq 8\frac{C_T}{\gamma} \|\nu^{\frac{1}{2}}\epsilon(\mathbf{v}_h)\|_{0,\Omega}^2 + \frac{1}{2}\|\gamma^{\frac{1}{2}}(\nu/\tilde{h})^{\frac{1}{2}}\mathbf{v}_h\|_{0,\partial\Omega}^2$$

In what follows we will assume that

$$(3.5) \qquad \gamma > 4C_T > 0,$$

and therefore

$$\lambda(\gamma) \stackrel{\text{def}}{=} 2 - 8 \frac{C_T}{\gamma} > 0.$$

From (3.4), we then get

$$\begin{aligned} a_h(\mathbf{v}_h, \mathbf{v}_h) + j_u(\mathbf{v}_h, \mathbf{v}_h) &\geq \|\sigma^{\frac{1}{2}}\mathbf{v}_h\|_{0,\Omega}^2 + \lambda(\gamma)\|\nu^{\frac{1}{2}}\epsilon(\mathbf{v}_h)\|_{0,\Omega}^2 + j_u(\mathbf{v}_h, \mathbf{v}_h) \\ &\quad + \frac{1}{2}\|\boldsymbol{\beta} \cdot \mathbf{n}\|_{\frac{1}{2}}^2 \mathbf{v}_h\|_{0,\partial\Omega}^2 + \frac{1}{2}\|\gamma^{\frac{1}{2}}(\nu/\tilde{h})^{\frac{1}{2}}\mathbf{v}_h\|_{0,\partial\Omega}^2 \\ &\quad + \|\gamma^{\frac{1}{2}} \max\{|\boldsymbol{\beta}|, \nu/\tilde{h}\}^{\frac{1}{2}}\mathbf{v}_h \cdot \mathbf{n}\|_{0,\partial\Omega}^2, \end{aligned}$$

and consequently

$$\begin{aligned} a_h(\mathbf{v}_h, \mathbf{v}_h) + j_u(\mathbf{v}_h, \mathbf{v}_h) &\geq \|\sigma^{\frac{1}{2}}\mathbf{v}_h\|_{0,\Omega}^2 + j_u(\mathbf{v}_h, \mathbf{v}_h) \\ &\quad + \min_{\substack{K \in \mathcal{T}_h \\ K \cap \partial\Omega \neq \emptyset}} \left\{ \lambda(\gamma), \frac{\gamma}{4h_K} \right\} \left( \|\nu^{\frac{1}{2}}\epsilon(\mathbf{v}_h)\|_{0,\Omega}^2 + \|\nu^{\frac{1}{2}}\mathbf{v}_h\|_{0,\partial\Omega}^2 \right) \\ &\quad + \frac{1}{2}\|\boldsymbol{\beta} \cdot \mathbf{n}\|_{\frac{1}{2}}^2 \mathbf{v}_h\|_{0,\partial\Omega}^2 + \frac{1}{4}\|\gamma^{\frac{1}{2}}(\nu/\tilde{h})^{\frac{1}{2}}\mathbf{v}_h\|_{0,\partial\Omega}^2 \\ &\quad + \|\gamma^{\frac{1}{2}} \max\{|\boldsymbol{\beta}|, \nu/\tilde{h}\}^{\frac{1}{2}}\mathbf{v}_h \cdot \mathbf{n}\|_{0,\partial\Omega}^2. \end{aligned}$$

In particular, since $0 < h \leq 1$ and by choosing (accordingly with (3.5))

$$\gamma \stackrel{\text{def}}{=} \varepsilon + 4C_T,$$

with $\varepsilon > 0$ sufficiently small, one obtains

$$\lambda(\gamma) < \frac{\gamma}{4h_K} \quad \forall K \in \mathcal{T}_h, \quad K \cap \partial\Omega \neq \emptyset.$$

We then conclude the proof using Korn's inequality (see [6]). □

<figure><img src="image_2.png" /></figure>

In what follows, we shall make use of the following discrete pressure and velocity subspaces:

$$C_{h,k}^{1} \stackrel{\text{def}}{=} \{q_h \in Q_h^k : j_p(q_h, q_h) = 0\},$$

$$V_{h,k}^{\text{div}} \stackrel{\text{def}}{=} \{\mathbf{v}_h \in [V_h^k]^d : b_h(q_h, \mathbf{v}_h) = 0 \quad \forall q_h \in C_{h,k}^1\}.$$

In addition, $Q_h^k \setminus C_{h,k}^1$ will stand for the supplementary of $C_{h,k}^1$ in $Q_h^k$, i.e.,

$$Q_h^k = (Q_h^k \setminus C_{h,k}^1) \oplus C_{h,k}^1.$$
The following lemma ensures, in particular, that $V_{h,k}^{\text{div}}$ is not trivial.

**LEMMA 3.4.** There exists a constant $\beta > 0$, independent of $h$, such that

$$ \inf_{q_h \in C_{h,k}^1} \sup_{\mathbf{v}_h \in [V_h^k]^d} \frac{|b_h(q_h, \mathbf{v}_h)|}{\|q_h\|_{0,\Omega} \|\mathbf{v}_h\|_{1,\Omega}} \geq \beta. $$

*Proof.* Let $q_h \in C_{h,k}^1$. From [27, Corollary 2.4], there exists $\mathbf{v}_q \in [H_0^1(\Omega)]^d$ such that

$$ (3.6) \qquad \nabla \cdot \mathbf{v}_q = q_h, \quad \|\mathbf{v}_q\|_{1,\Omega} \leq C \|q_h\|_{0,\Omega}. $$

Thus, using integration by parts and (2.12), we have

$$ \begin{align} \|q_h\|_{0,\Omega}^2 &= (q_h, \nabla \cdot \mathbf{v}_q) \nonumber \\ &= (q_h, \nabla \cdot \mathbf{v}_q - \nabla \cdot \pi_{h,k} \mathbf{v}_q) + (q_h, \nabla \cdot \pi_{h,k} \mathbf{v}_q) \tag{3.7} \\ &= (\nabla q_h, \mathbf{v}_q - \pi_{h,k} \mathbf{v}_q) - \langle q_h, (\pi_{h,k} \mathbf{v}_q) \cdot \mathbf{n} \rangle_{\partial \Omega} \nonumber \\ &\quad + (q_h, \nabla \cdot \pi_{h,k} \mathbf{v}_q) \nonumber \\ &= (\nabla q_h, \mathbf{v}_q - \pi_{h,k} \mathbf{v}_q) - b_h(q_h, \pi_{h,k} \mathbf{v}_q). \nonumber \end{align} $$

Since $q_h \in C_{h,k}^1$, it follows that $\nabla q_h \in [V_h^k]^d$. Thus, using the orthogonality of the $L^2$-projection, we have

$$ (\nabla q_h, \mathbf{v}_q - \pi_{h,k} \mathbf{v}_q) = 0. $$

Thus, from (3.7), it follows that

$$ |b_h(q_h, \pi_{h,k} \mathbf{v}_q)| = \|q_h\|_{0,\Omega}^2. $$

In addition, using $H^1$-stability of the $L^2$-projection (see [5]) and (3.6), we have

$$ \begin{align*} \|\pi_{h,k} \mathbf{v}_q\|_{1,\Omega} &\leq C \|\mathbf{v}_q\|_{1,\Omega} \\ &\leq C \|q_h\|_{0,\Omega}, \end{align*} $$

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

$$ (3.8) \qquad \begin{aligned} A \mathbf{u}_h + B^T p_h &= M \mathbf{f} && \text{in } ([V_h^k]^d)', \\ B \mathbf{u}_h &= J p_h && \text{in } [Q_h^k]', \end{aligned} $$

with $A \in \mathcal{L}([V_h^k]^d, ([V_h^k]^d)'$), $M \in \mathcal{L}([L^2(\Omega)]^d, ([V_h^k]^d)'$), $B \in \mathcal{L}([V_h^k]^d, (Q_h^k)')$, and $J \in \mathcal{L}(Q_h^k, (Q_h^k)')$ defined by

$$ \begin{align*} \langle A \mathbf{u}_h, \mathbf{v}_h \rangle &\stackrel{\text{def}}{=} a_h(\mathbf{u}_h, \mathbf{v}_h) + j_u(\mathbf{u}_h, \mathbf{v}_h), \\ \langle M \mathbf{f}, \mathbf{v}_h \rangle &\stackrel{\text{def}}{=} (\mathbf{f}, \mathbf{v}_h), \\ \langle B \mathbf{v}_h, q_h \rangle &\stackrel{\text{def}}{=} b_h(q_h, \mathbf{v}_h), \\ \langle J p_h, q_h \rangle &\stackrel{\text{def}}{=} j(p_h, q_h). \end{align*} $$
We also introduce the operator $B^1 \in \mathcal{L}([V_h^k]^d, (C_{h,k}^1)')$ defined by

$$ \langle B^1 \mathbf{v}_h, q_h \rangle \stackrel{\text{def}}{=} b_h(q_n, \mathbf{v}_h) \quad \forall (\mathbf{v}_h, q_h) \in [V_h^k]^d \times C_{h,k}^1; $$

in other words,

$$ B^1 \mathbf{v}_h \stackrel{\text{def}}{=} (B \mathbf{v}_h)|_{C_{h,k}^1} \quad \forall \mathbf{v}_h \in [V_h^k]^d. $$

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_{h,k}^{\text{div}} \stackrel{\text{def}}{=} \text{Ker}(B^1) \neq \{0\}$.

Let us consider the following reduced formulation (derived from (2.10) with $(\mathbf{v}_h, q_h) \in V_{h,k}^{\text{div}} \times (Q_h^k \setminus C_{h,k}^1)$): find $(\mathbf{u}_h, \tilde{p}_h) \in V_{h,k}^{\text{div}} \times (Q_h^k \setminus C_{h,k}^1)$ such that

$$ \begin{array}{ll} (3.9) & A\mathbf{u}_h + B^T \tilde{p}_h = M\mathbf{f} \quad \text{in} \quad (V_{h,k}^{\text{div}})' \\ & B\mathbf{u}_h = J\tilde{p}_h \quad \text{in} \quad (Q_h^k \setminus C_{h,k}^1)' \end{array} $$

Since, by construction, $C_{h,k}^1 = \text{Ker}(J)$, we conclude that $J$ is invertible in $Q_h^k \setminus C_{h,k}^1$. Hence, from (3.9), we have

$$ (3.10) \qquad \tilde{p}_h = J_{|Q_h^k \setminus C_{h,k}^1}^{-1} B \mathbf{u}_h. $$

By plugging this expression into the first equation of (3.9), we obtain that $\mathbf{u}_h \in V_{h,k}^{\text{div}}$ solves

$$ (A + B^T J_{|Q_h^k \setminus C_{h,k}^1}^{-1} B) \mathbf{u}_h = M\mathbf{f} \quad \text{in} \quad (V_{h,k}^{\text{div}})' $$

Existence and uniqueness of $\mathbf{u}_h$ follow by the positivity of $A$ (Lemma 3.3) and the nonnegativity of $B^T J_{|Q_h^k \setminus C_{h,k}^1}^{-1} B$. We may then recover $\tilde{p}_h$ 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

$$ A\mathbf{u}_h + B^T \tilde{p}_h - M\mathbf{f} \in (\text{Ker}(B^1))^0, $$

with $(\text{Ker}(B^1))^0$ standing for the polar set of $\text{Ker}(B^1)$. From Lemma 3.4, it follows that $B^1$ is an isomorphism from $C_{h,k}^1$ onto $(\text{Ker}(B^1))^0$ (see [27, p. 58]). Thus, there exists a unique $p^1 \in C_{h,k}^1$ such that

$$ (3.11) \qquad A\mathbf{u}_h + B^T \tilde{p}_h - M\mathbf{f} = (B^1)^T p^1 \quad \text{in} \quad ([V_h^k]^d)' $$

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 $(\mathbf{u}_h, p_h \stackrel{\text{def}}{=} \tilde{p}_h - p^1)$. $\square$

**4. Convergence of the method.** The parameter for the pressure stabilization scales as $h_K^2/\|\boldsymbol{\beta}\|_{0,\infty,K}$ when the local Reynolds number $\text{Re}_K$ is big, and as $h_K^3/\nu$ when $\text{Re}_K$ is small. The stabilizing terms acting on the velocity scale as $\|\boldsymbol{\beta}\|_{0,\infty,K}h_K^2$ at a high local Reynolds number and as $\text{Re}_K \|\boldsymbol{\beta}\|_{0,\infty,K}h_K^2$ for a low Reynolds number. The factor $\text{Re}_K \|\boldsymbol{\beta}\|_{0,\infty,K}h_K^2$ 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, $(\mathbf{u}, p) \in [H^{k+1}(\Omega)]^{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 $\rho, \eta > 0$ sufficiently small and for all $\phi \in V_h^1$ satisfying

$$ \phi > 0, \quad |\nabla\phi(\mathbf{x})| \le \eta h_K^{-1} \phi(\mathbf{x}) \quad \forall \mathbf{x} \in K, \quad \forall K \in \mathcal{T}_h, $$

there holds

$$ \begin{aligned} \|\phi\pi_{h,k}u\|_{0,\Omega} &\le C\|\phi u\|_{0,\Omega} \quad \forall u \in L^2(\Omega), \\ \|\phi\nabla\pi_{h,k}u\|_{0,\Omega} &\le C\|\phi\nabla u\|_{0,\Omega} \quad \forall u \in H^1(\Omega). \end{aligned} $$

A direct consequence of this result is stated in the following corollary.

**COROLLARY 4.3.** Under the assumptions of the previous lemma, we have

$$ \left( \sum_{|\alpha| \le l} \|\phi \partial^{\alpha} (u - \pi_{h,k} u)\|_{0,\Omega}^2 \right)^{\frac{1}{2}} \le C \left( \sum_{K \in \mathcal{T}_h} \|\phi\|_{0,\infty,K}^2 h^{2(r_u-l)} \|u\|_{r_u,\Omega}^2 \right)^{\frac{1}{2}} $$

for all $u \in H^r(\Omega)$, with $r \ge 1$, $r_u \stackrel{\text{def}}{=} \min\{r, k+1\}$, $0 \le l \le r_u$, $\alpha \in \mathbb{N}^d$, and $\partial^\alpha$ 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.

**LEMMA 4.4.** Let $\phi_i \in H^2(\mathcal{T}_h)$, $i = 1, \dots, 5$, be piecewise constant functions defined by

$$ \begin{aligned} \phi_{1|K} &\stackrel{\text{def}}{=} \nu^{-\frac{1}{2}} \min \{\mathrm{Re}_K^{-\frac{1}{2}}, 1\}, & \phi_{2|K} &\stackrel{\text{def}}{=} \|\beta\|_{0,\infty,K}^{\frac{1}{2}} h_K^{-\frac{1}{2}}, \\ \phi_{3|K} &\stackrel{\text{def}}{=} h_K^{-\frac{1}{2}} \|\beta\|_{0,\infty,K}^{\frac{1}{2}} \xi(\mathrm{Re}_K)^{-\frac{1}{2}}, & \phi_{4|K} &\stackrel{\text{def}}{=} \phi_{3|K}^{-1}, & \phi_{5|K} &\stackrel{\text{def}}{=} h_K^{-r}, \end{aligned} $$

with $r \ge 1$, for all $K \in \mathcal{T}_h$, and let $\phi_i^* \stackrel{\text{def}}{=} \pi_{h,1}^* \phi_i$. Then, there holds

$$ \begin{aligned} \phi_i (\rho \beta \rho)^{-\frac{1}{2}} &\le \phi_i^* \le (\rho \beta \rho)^{\frac{1}{2}} \phi_i \\ |\nabla \phi_i^*| &\le c_0 (\rho \beta \rho - 1) h_K^{-1} \phi_i^* \\ \phi_i \rho \beta^{-\frac{1}{2}} \rho^{-1} &\le \phi_i^* \le \rho \beta^{-\frac{1}{2}} \phi_i \\ |\nabla \phi_i^*| &\le c_0 (\rho \beta \rho^2 - 1) h_K^{-1} \phi_i^* \\ \phi_5 \rho^{-r} &\le \phi_5^* \le \rho^r \phi_5, \\ |\nabla \phi_5^*| &\le c_0 (\rho^{2r} - 1) h_K^{-1} \phi_5^* \end{aligned} \quad \text{for } i=1,2, \\ \phi_i \rho \beta^{-\frac{1}{2}} \rho^{-1} &\le \phi_i^* \le \rho \beta^{-\frac{1}{2}} \phi_i \quad \text{for } i=3,4, \\ |\nabla \phi_i^*| &\le c_0 (\rho \beta \rho^2 - 1) h_K^{-1} \phi_i^* \\ \phi_5 \rho^{-r} &\le \phi_5^* \le \rho^r \phi_5, \\ |\nabla \phi_5^*| &\le c_0 (\rho^{2r} - 1) h_K^{-1} \phi_5^* \end{aligned} $$

in $K$, for all $K \in \mathcal{T}_h$, with $c_0 > 0$ the constant in (2.4).
*Proof.* We give the proof only for $\phi_1$; the argument for the rest is similar. First, note that for all $K \in \mathcal{T}_h$,

$$ (4.1) \qquad \begin{aligned} \max_{x \in K} \phi_1^* &\le \max_{K' \in \mathcal{N}(K)} \nu^{-\frac{1}{2}} \min\{\mathrm{Re}_{K'}^{-\frac{1}{2}}, 1\} \\ &= \max_{K' \in \mathcal{N}(K)} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K'}^{-\frac{1}{2}} h_{K'}^{-\frac{1}{2}}, \nu^{-\frac{1}{2}}\}. \end{aligned} $$

We now distinguish two cases. On one hand, if $\mathrm{Re}_K \le 1$, we have $\phi_{1|K} = \nu^{-\frac{1}{2}}$. Thus, from (4.1) and since $\rho\beta\rho > 1$, it follows that

$$ \max_{x \in K} \phi_1^* \le \nu^{-\frac{1}{2}} \\ \le (\rho \rho \beta)^{\frac{1}{2}} \phi_{1|K}. $$

On the other hand, if $\mathrm{Re}_K > 1$, we get $\phi_{1|K} = \|\boldsymbol{\beta}\|_{0,\infty,K}^{-\frac{1}{2}} h_K^{-\frac{1}{2}}$. Therefore, from (4.1) and the assumptions on the mesh (2.5) and (2.6), we have

$$ \begin{aligned} \max_{x \in K} \phi_1^* &\le \max_{K' \in \mathcal{N}(K)} \{\|\boldsymbol{\beta}\|_{0,\infty,K'}^{-\frac{1}{2}} h_{K'}^{-\frac{1}{2}}\} \\ &\le (\rho \rho \beta)^{\frac{1}{2}} \|\boldsymbol{\beta}\|_{0,\infty,K}^{-\frac{1}{2}} h_K^{-\frac{1}{2}} \\ &= (\rho \rho \beta)^{\frac{1}{2}} \phi_{1|K}. \end{aligned} $$

The lower bound follows in a similar fashion.

Finally, for the derivative, using the bounds on $\phi_1^*$ and the regularity of the mesh (2.4), and since $\rho\beta\rho > 1$, we obtain

$$ \begin{aligned} |\nabla \phi_1^*|_K | &\le \frac{\max_{x \in K} \phi_1^* - \min_{x \in K} \phi_1^*}{\rho_K} \\ &\le \frac{(\rho \beta \rho)^{\frac{1}{2}} - (\rho \beta \rho)^{-\frac{1}{2}}}{\rho_K} \phi_{1|K} \\ &\le c_0 (\rho \beta \rho - 1) h_K^{-1} \phi_1^*, \end{aligned} $$

which completes the proof. □

**REMARK 4.5.** It follows from Lemma 4.4 that for the weight functions $\phi_i^*$, $1 \le i \le 5$, the stability estimate of Lemma 4.2 holds, provided $\rho\beta$ and $\rho$ 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 $\|\cdot\|_K$.

**LEMMA 4.6 (velocity approximability).** Assume that $\rho\beta$ and $\rho$ are sufficiently close to 1. Then, there holds

$$ \|\mathbf{u} - \pi_{h,k} \mathbf{u}\|^2 \le C \sum_{K \in \mathcal{T}_h} \left( \sigma h_K^{2r_u} + \max\{\nu, \|\boldsymbol{\beta}\|_{0,\infty,K} h_K\} h_K^{2(r_u-1)} \right) \|\mathbf{u}\|_{r_u,K}^2 $$

for all $\mathbf{u} \in [H^r(\Omega)]^d$, with $r \ge 2$ and $r_u = \min\{k+1, r\}$.

*Proof.* First note that

$$ \|\mathbf{u} - \pi_{h,k} \mathbf{u}\|^2 \le \|i_{h,k} \mathbf{u} - \pi_{h,k} \mathbf{u}\|^2 + \|i_{h,k} \mathbf{u} - \mathbf{u}\|^2. $$
We give the proof for the first term only. The argument for the second term is similar. By the stability estimate for the L²-projection Lemma 4.2 we have

$$
\begin{align*}
\|\sigma^{\frac{1}{2}}(i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u})\|_{0,\Omega}^2 &\le C \|\sigma^{\frac{1}{2}}(\mathbf{u} - i_{h,k}\mathbf{u})\|_{0,\Omega}^2 \\
&\le C \sum_{K \in \mathcal{T}_h} \sigma h_K^{2r_u} \|\mathbf{u}\|_{r_u,K}^2.
\end{align*}
$$

Using now the $H^1$-stability of the $L^2$-projection on locally quasi-uniform meshes (see [5]) we get

$$
\begin{align*}
\|\nu^{\frac{1}{2}} \nabla (i_{h,k} \mathbf{u} - \pi_{h,k} \mathbf{u})\|_{0,\Omega}^2 &\le C \|\nu^{\frac{1}{2}} \nabla (i_{h,k} \mathbf{u} - \mathbf{u})\|_{0,\Omega}^2 \\
&\le C \sum_{K \in \mathcal{T}_h} \nu h^{2(r_u-1)} \|\mathbf{u}\|_{r_u,\Omega}.
\end{align*}
$$

We treat the boundary terms using the trace inequality (2.8) in combination with Lemma 4.2 and approximation, which yields

$$
\begin{align*}
(4.2) \quad & \| \max\{|\boldsymbol{\beta}|, \nu/\tilde{h}\}^{\frac{1}{2}} (i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}) \|_{0,\partial\Omega}^2 \\
& \le C \sum_{\substack{K \in \mathcal{T}_h \\ K \cap \partial\Omega \neq \emptyset}} \| \max\{\|\boldsymbol{\beta}\|_{0,\infty,K}, \nu h_K^{-1}\}^{\frac{1}{2}} (i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}) \|_{0,K \cap \partial\Omega}^2 \\
& \le C \sum_{K \in \mathcal{T}_h} \| h_K^{-\frac{1}{2}} \max\{\|\boldsymbol{\beta}\|_{0,\infty,K}, \nu h_K^{-1}\}^{\frac{1}{2}} (i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}) \|_{0,K}^2 \\
& \le C \sum_{K \in \mathcal{T}_h} \| \phi_3^*(i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}) \|_{0,K}^2 \\
& = C \sum_{K \in \mathcal{T}_h} \| \phi_3^* \pi_{h,k} (i_{h,k}\mathbf{u} - \mathbf{u}) \|_{0,K}^2 \\
& \le C \sum_K \max\{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K, \nu \} h_K^{2(r_u-1)} \|\mathbf{u}\|_{r_u,K}^2.
\end{align*}
$$

The interior penalty terms are treated in the same fashion as the boundary terms.
We have

$$
(4.3) \qquad
\begin{aligned}
j_{\mathbf{u}}(\mathbf{u} - \pi_{h,k}\mathbf{u}, \mathbf{u} - \pi_{h,k}\mathbf{u}) &\le j_{\mathbf{u}}(i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}, i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}) \\
&\quad + j_{\mathbf{u}}(\mathbf{u} - i_{h,k}\mathbf{u}, \mathbf{u} - i_{h,k}\mathbf{u}).
\end{aligned}
$$

The first term in this inequality can be estimated using that $\xi(\mathrm{Re}_K) \le 1$, the trace inequality (2.8), an inverse inequality, and the $H^1$-stability of the $L^2$-projection (see [5]), which yields

$$
\begin{align*}
(4.4) \quad & j_{\mathbf{u}}(i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}, i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}) \le C \sum_{K \in \mathcal{T}_h} \|\boldsymbol{\beta}\|_{0,\infty,K} \xi(\mathrm{Re}_K) h_K^2 \|\nabla(i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u})\|_{0,\partial K}^2 \\
& \le C \sum_{K \in \mathcal{T}_h} \|\boldsymbol{\beta}\|_{0,\infty,K} h_K^{-1} \|i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u}\|_{0,K}^2 \\
& \le C \|\phi_2^*(i_{h,k}\mathbf{u} - \pi_{h,k}\mathbf{u})\|_{0,\Omega}^2 \\
& \le C \sum_K \|\boldsymbol{\beta}\|_{0,\infty,K} h_K^{2r_u-1} \|\mathbf{u}\|_{r_u,K}^2,
\end{align*}
$$

and so the proof is finished. □

<figure><img src="image_3.png" /></figure>
For the pressure, we have the following result.

**LEMMA 4.7 (pressure approximability).** *Under the assumptions of Lemma 4.6, there holds*

$$ \begin{aligned} & \|\tilde{h}^{\frac{1}{2}} \phi_1^*(p - \pi_{h,k} p)\|_{0,\partial\Omega}^2 + \|\phi_1^*(p - \pi_{h,k} p)\|_{0,\Omega}^2 + j(\pi_{h,k} p, \pi_{h,k} p) \\ & \le C \sum_{K \in \mathcal{T}_h} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K}^{-1}, h_K/\nu\} h_K^{2s_p-1} \|p\|_{s_p,K}^2 \end{aligned} $$

for all $p \in H^s(\Omega)$ with $s \ge 1$ and $s_p \stackrel{\text{def}}{=} \min\{k+1, s\}$.

*Proof.* As $p$ may be only $H^1(\Omega)$, 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} \nabla p] = 0$ and using a trace inequality followed by an inverse inequality, we have

$$ \begin{aligned} j(\pi_{h,k} p, \pi_{h,k} p) &= \sum_{K \in \mathcal{T}_h} \xi(\mathrm{Re}_K) \frac{h_K^2}{\|\boldsymbol{\beta}\|_{0,\infty,K}} \|[\nabla \pi_{h,k} p - C_{h,k} \nabla p]\|^2 \\ &\le C \sum_{K \in \mathcal{T}_h} \xi(\mathrm{Re}_K) \frac{h_K}{\|\boldsymbol{\beta}\|_{0,\infty,K}} \|\nabla \pi_{h,k} p - C_{h,k} \nabla p\|_{0,K}^2 \\ &\le C (\|\phi_4^* \nabla \pi_{h,k} (p - C_{h,k} p)\|_{0,K}^2 + \|\phi_4^* \nabla (p - C_{h,k} p)\|_{0,K}^2 \\ &\quad + \|\phi_4^* (\nabla p - C_{h,k} \nabla p)\|_{0,K}^2) \\ &\le C \sum_{K \in \mathcal{T}_h} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K}^{-1} h_K, h_K^2/\nu\} h_K^{2(s_p-1)} \|p\|_{s_p,K}^2, \end{aligned} $$

where we concluded using the stability lemma, Lemma 4.2, with weight function $\phi_4^*$, and the optimal approximation properties of the Clément interpolant (see [16, 21]). □

**4.1. 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 $\mathbf{v}_h \in [V_h^k]^d$, there holds*

$$ \begin{aligned} \sum_{K \in \mathcal{T}_h} h_K^2 \int_{\partial K} \|\boldsymbol{\beta}_h \cdot \mathbf{n}\|_{0,\infty,\partial K} \|[\mathbf{n} \cdot \nabla \mathbf{v}_h]\|^2 \mathrm{d}s &\le C \left( j_{\mathbf{u}}(\mathbf{v}_h, \mathbf{v}_h) + \|\nu^{\frac{1}{2}} \nabla \mathbf{v}_h\|_{0,\Omega}^2 \right), \\ \sum_{K \in \mathcal{T}_h} \phi_{1|K}^{-1} h_K \int_{\partial K} \|[\nabla \cdot \mathbf{v}_h]\|^2 \mathrm{d}s &\le C \left( j_{\mathbf{u}}(\mathbf{v}_h, \mathbf{v}_h) + \|\nu^{\frac{1}{2}} \nabla \mathbf{v}_h\|_{0,\Omega}^2 \right). \end{aligned} $$

*Proof.* Let $A_1$ denote the set of elements $K \in \mathcal{T}_h$ such that $\xi(\mathrm{Re}_K) \ge 1$, and $A_2$ the set of elements such that $\xi(\mathrm{Re}_K) < 1$. It then follows that $|\boldsymbol{\beta}|_{\infty,0,K} h_K < \nu$ for $K \in A_2$, and we may write
$$
\begin{align*}
& \sum_{K \in \mathcal{T}_h} h_K^2 \int_{\partial K} \|\boldsymbol{\beta}_h \cdot \mathbf{n}\|_{0,\infty,\partial K} |\llbracket \mathbf{n} \cdot \nabla \mathbf{v}_h \rrbracket|^2 \, ds \\
& \le \sum_{K \in A_1} h_K^2 \int_{\partial K} \|\boldsymbol{\beta}_h \cdot \mathbf{n}\|_{0,\infty,\partial K} |\llbracket \mathbf{n} \cdot \nabla \mathbf{v}_h \rrbracket|^2 \, ds + \sum_{K \in A_2} h_K \nu \int_{\partial K} |\llbracket \mathbf{n} \cdot \nabla \mathbf{v}_h \rrbracket|^2 \, ds \\
& \le \sum_{K \in \mathcal{T}_h} h_K^2 \xi(\mathrm{Re}_K) \|\boldsymbol{\beta}_h \cdot \mathbf{n}\|_{0,\infty,\partial K} |\llbracket \mathbf{n} \cdot \nabla \mathbf{v}_h \rrbracket|^2 \, ds + C \|\nu^{\frac{1}{2}} \nabla \mathbf{v}_h\|_{0,\Omega}^2,
\end{align*}
$$

where the last inequality follows by a trace inequality, an inverse inequality in the
second term, and extending the sums over all $\mathcal{T}_h$.

The second inequality follows in a similar fashion, noting that

$$
\begin{align*}
\sum_{K \in \mathcal{T}_h} \phi_{1|K}^{-1} h_K \int_{\partial K} |\nabla \cdot \mathbf{v}_h|^2 ds &\le C \sum_{K \in \mathcal{T}_h} h_K \max\{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K, \nu\} \int_{\partial K} |\nabla \cdot \mathbf{v}_h|^2 ds \\
&\le \sum_{K \in A_1} h_K^2 \|\boldsymbol{\beta}\|_{0,\infty,K} \xi(\mathrm{Re}_K) \int_{\partial K} |\nabla \cdot \mathbf{v}_h|^2 ds + \|\nu^{\frac{1}{2}} \nabla \mathbf{v}_h\|_{0,A_2}^2,
\end{align*}
$$

and so the proof is completed. $\square$

The main result of this paragraph is stated in the following theorem.

THEOREM 4.9. Assume $(\mathbf{u}, p) \in [H^r(\Omega)]^d \times H^s(\Omega)$, with $r \ge 2$ and $s \ge 1$, is the solution of (2.1) and $(\mathbf{u}_h, p_h) \in W_h^k$ is the solution of (2.10). Then, under the assumptions of Lemma 4.6, there holds

$$
\|\mathbf{u} - \mathbf{u}_h\| \le C \left[ \sum_{K \in \mathcal{K}} \left( \sigma h_K^{2r_u} + \max\{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K, \nu\} h_K^{2(r_u-1)} \right) \|\mathbf{u}\|_{r_u, K}^2 \right]^{\frac{1}{2}} \\
+ C \max_{K \in \mathcal{T}_h} \left\{ \sigma^{-\frac{1}{2}} |\boldsymbol{\beta}|_{1,\infty,K} h_K^{r_u} \right\} \|\mathbf{u}\|_{r_u, \Omega} + C \left( \sum_{K \in \mathcal{T}_h} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K}^{-1}, h_K/\nu \} h_K^{2s_p-1} \|p\|_{s_p, K}^2 \right)^{\frac{1}{2}},
$$

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:

$$
\mathbf{u} - \mathbf{u}_h = \underbrace{\mathbf{u} - \pi_{h,k} \mathbf{u}}_{\mathbf{e}^{\pi}} + \underbrace{\pi_{h,k} \mathbf{u} - \mathbf{u}_h}_{-\mathbf{e}_h} = \mathbf{e}^{\pi} - \mathbf{e}_h.
$$

We also consider the discrete pressure error

$$
(4.5) \qquad y_h \stackrel{\text{def}}{=} p_h - \pi_{h,k} p.
$$

It follows then that

$$
\| \mathbf{u} - \mathbf{u}_h \| \le \| \mathbf{e}^\pi \| + \| \mathbf{e}_h \|.
$$

Lemma 4.6 gives an estimate for $\|\mathbf{e}^\pi\|$. Hence, it suffices to estimate $\|\mathbf{e}_h\|$.

Using coercivity and orthogonality, namely, Lemmas 3.3 and 2.1, we get

$$
\begin{equation} \tag{4.6}
\begin{aligned}
C \|\mathbf{e}_h\|^2 + j_p(y_h, y_h) &\le a_h(\mathbf{e}_h, \mathbf{e}_h) + j_\mathbf{u}(\mathbf{e}_h, \mathbf{e}_h) + j_p(y_h, y_h) \\
&= a_h(\mathbf{e}^\pi, \mathbf{e}_h) + b_h(y^\pi, \mathbf{e}_h) - b_h(y_h, \mathbf{e}^\pi) \\
&\quad + j_\mathbf{u}(\mathbf{e}^\pi, \mathbf{e}_h) - j_p(\pi_{h,k}p, y_h).
\end{aligned}
\end{equation}
$$
By an application of the Cauchy–Schwarz inequality in the symmetric part of the discrete elliptic operator and integrating by parts in the convective term, we obtain

$$a_h(\mathbf{e}^{\pi}, \mathbf{e}_h) \le \|\mathbf{e}^{\pi}\| \|\mathbf{e}_h\| + |(\mathbf{e}^{\pi}, \boldsymbol{\beta} \cdot \nabla \mathbf{e}_h)| \\ - \langle 2\nu\boldsymbol{\epsilon}(\mathbf{e}^{\pi})\mathbf{n}, \mathbf{e}_h \rangle_{\partial\Omega} - \langle \mathbf{e}^{\pi}, 2\nu\boldsymbol{\epsilon}(\mathbf{e}_h)\mathbf{n} \rangle_{\partial\Omega},$$

where, for simplicity, the boundary term from the integration by parts has been included in the first term on the right-hand side. We note that in the same way we have, using the Cauchy–Schwarz inequality, a trace inequality and a local inverse inequality,

$$(4.7) \qquad \langle 2\nu\boldsymbol{\epsilon}(\mathbf{e}_h)\mathbf{n}, \mathbf{e}^{\pi} \rangle_{\partial\Omega} \le C \|\mathbf{e}_h\| \|\mathbf{e}^{\pi}\|.$$

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

$$(4.8) \qquad \langle 2\nu\boldsymbol{\epsilon}(\mathbf{e}^{\pi})\mathbf{n}, \mathbf{e}_h \rangle_{\partial\Omega} \le C \left( \sum_{K \in \mathcal{T}_h} \nu h_K^{2(r_u-1)} \|\mathbf{u}\|_{r_u, K}^2 \right)^{\frac{1}{2}} \|\mathbf{e}_h\|.$$

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 $\boldsymbol{\beta}$ by its piecewise linear interpolant $\beta_h$,

$$|(\mathbf{e}^{\pi}, \boldsymbol{\beta} \cdot \nabla \mathbf{e}_h)| \le |(\mathbf{e}^{\pi}, (\boldsymbol{\beta} - \beta_h) \cdot \nabla \mathbf{e}_h)| + |(\mathbf{e}^{\pi}, \beta_h \cdot \nabla \mathbf{e}_h)| \\ \le C \sum_{K \in \mathcal{T}_h} |\boldsymbol{\beta}|_{1, \infty, K} \|\mathbf{e}^{\pi}\|_{0, K} h_K \|\nabla \mathbf{e}_h\|_{0, K} \\ + |(\mathbf{e}^{\pi}, \beta_h \cdot \nabla \mathbf{e}_h - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{e}_h))| \\ \le C \sum_{K \in \mathcal{T}_h} \sigma^{-\frac{1}{2}} |\boldsymbol{\beta}|_{1, \infty, K} h_K^{r_u} \|\phi_5^* \mathbf{e}^{\pi}\|_{0, K} \|\sigma^{\frac{1}{2}} \mathbf{e}_h\|_{0, K} \\ + |(\mathbf{e}^{\pi}, \beta_h \cdot \nabla \mathbf{e}_h - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{e}_h))| \\ \le C \max_{K \in \mathcal{T}_h} \{\sigma^{-\frac{1}{2}} |\boldsymbol{\beta}|_{1, \infty, K} h_K^{r_u}\} \|\mathbf{u}\|_{r_u, \Omega} \|\mathbf{e}_h\| \\ + \|\phi_2 \mathbf{e}^{\pi}\|_{0, \Omega} \|\phi_2^{-1} (\beta_h \cdot \nabla \mathbf{e}_h - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{e}_h))\|_{0, \Omega}.$$

Now we apply Lemma 3.1 to obtain

$$\|\phi_2 \mathbf{e}^{\pi}\|_{0, \Omega} \|\phi_2^{-1} (\beta_h \cdot \nabla \mathbf{e}_h - \pi_{h,k}^*(\beta_h \cdot \nabla \mathbf{e}_h))\|_{0, \Omega} \\ \le C \|\phi_2^* \mathbf{e}^{\pi}\|_{0, \Omega} \left( \sum_{K \in \mathcal{T}_h} \int_{\partial K} h_K^2 \|\boldsymbol{\beta} \cdot \mathbf{n}\|_{0, \infty, \partial K} \|[\mathbf{n} \cdot \nabla \mathbf{u}_h]\|^2 ds \right)^{\frac{1}{2}} \\ \le C \left( \sum_{K \in \mathcal{T}_h} \|\boldsymbol{\beta}\|_{0, \infty, K} h_K^{2(r_u-1)} \|\mathbf{u}\|_{r_u, K}^2 \right)^{\frac{1}{2}} \|\mathbf{e}_h\|,$$

where we used Corollary 4.3 and Lemma 4.8 in the last inequality.

Collecting terms we have

$$(4.9) \qquad a_h(\mathbf{e}^{\pi}, \mathbf{e}_h) \le C \|\mathbf{e}^{\pi}\| \|\mathbf{e}_h\| + C \max_{K \in \mathcal{T}_h} \{\sigma^{-\frac{1}{2}} |\boldsymbol{\beta}|_{1, \infty, K} h_K^{r_u}\} \|\mathbf{e}_h\| \\ + C \left( \sum_{K \in \mathcal{T}_h} \max\{\|\boldsymbol{\beta}\|_{0, \infty, K} h_K, \nu\} h_K^{2(r_u-1)} \|\mathbf{u}\|_{r_u, K}^2 \right)^{\frac{1}{2}} \|\mathbf{e}_h\|.$$
For the second term in (4.6), using the orthogonality of the L²-projection, Lemmas 4.7 and 4.8, and replacing **u** with **p** in (4.2), we have

$$
\begin{equation} \tag{4.10}
\begin{split}
b_h(y^\pi, \mathbf{e}_h) &= -(y^\pi, \nabla \cdot \mathbf{e}_h - \pi_{h,k}^*(\nabla \cdot \mathbf{e}_h)) + \langle y^\pi, \mathbf{e}_h \cdot \mathbf{n} \rangle_{\partial\Omega} \\
&\leq \|\phi_1^* y^\pi\|_{0,\Omega} \|\phi_1^{-1}(\nabla \cdot \mathbf{e}_h - \pi_{h,k}^*(\nabla \cdot \mathbf{e}_h))\|_{0,\Omega} + \|\phi_1^* \tilde{h}^{\frac{1}{2}} y^\pi\|_{0,\partial\Omega} \|\mathbf{e}_h\| \\
&\leq C \left( \sum_{K \in \mathcal{T}_h} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K}^{-1}, h_K/\nu\} h_K^{2s_p-1} \|p\|_{s_p,K}^2 \right)^{\frac{1}{2}} \|\mathbf{e}_h\|.
\end{split}
\end{equation}
$$

In a similar fashion, after integration by parts in the third term, one obtains

$$
\begin{equation} \tag{4.11}
\begin{split}
b_h(y_h, \mathbf{e}^\pi) &= -(y_h, \nabla \cdot \mathbf{e}^\pi) + \langle y_h, \mathbf{e}^\pi \cdot \mathbf{n} \rangle_{\partial\Omega} \\
&= (\nabla y_h, \mathbf{e}^\pi) \\
&= (\nabla y_h - \pi_{h,k}^*(\nabla y_h), \mathbf{e}^\pi) \\
&\le C \|\phi_3^{-1}(\nabla y_h - \pi_{h,k}^*(\nabla y_h))\|_{0,\Omega} \|\phi_3^* \mathbf{e}^\pi\|_{0,\Omega} \\
&\le C j_p(y_h, y_h)^{\frac{1}{2}} \left( \sum_{K \in \mathcal{T}_h} \max\{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K, \nu\} h_K^{2(r_u-1)} \|\mathbf{u}\|_{r_u,K}^2 \right)^{\frac{1}{2}}.
\end{split}
\end{equation}
$$

Finally, using Lemma 4.7, for the interior penalty terms we have

$$
\begin{equation} \tag{4.12}
\begin{split}
j_{\mathbf{u}}(\mathbf{e}^{\pi}, \mathbf{e}_h) + j_p(\pi_{h,k}p, y_h) &\le C \|\mathbf{e}^{\pi}\| \|\mathbf{e}_h\| + j_p(\pi_{h,k}p, \pi_{h,k}p)^{\frac{1}{2}} j_p(y_h, y_h)^{\frac{1}{2}} \\
&\le C \|\mathbf{e}^{\pi}\| \|\mathbf{e}_h\| + C j_p(y_h, y_h)^{\frac{1}{2}} \left( \sum_{K \in \mathcal{T}_h} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K}^{-1}, h_K/\nu\} h_K^{2s_p-1} \|p\|_{s_p,K}^2 \right)^{\frac{1}{2}}.
\end{split}
\end{equation}
$$

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.

COROLLARY 4.10. Under the assumptions of Theorem 4.9, there holds

$$
\begin{align*}
j_p(p_h, p_h) &\le C \max_{K \in \mathcal{T}_h} \left\{ \sigma^{-\frac{1}{2}} |\boldsymbol{\beta}|_{1,\infty,K} h_K^{r_u} \right\} \|\mathbf{u}\|_{r_u, \Omega} \\
&\quad + C \left[ \sum_{K \in \mathcal{K}} \left( \sigma h_K^{2r_u} + \max\{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K, \nu\} h^{2(r_u-1)} \right) \|\mathbf{u}\|_{r_u, K}^2 \right]^{\frac{1}{2}} \\
&\quad + C \left( \sum_{K \in \mathcal{T}_h} \min\{\|\boldsymbol{\beta}\|_{0,\infty,K}^{-1}, h_K/\nu\} h_K^{2s_p-1} \|p\|_{s_p, K}^2 \right)^{\frac{1}{2}}.
\end{align*}
$$

REMARK 4.11. In a physical situation the velocity gradient on boundaries with no-slip conditions is known to scale as $|\boldsymbol{\beta}|_{1,\infty,\partial\Omega} \sim \nu^{-\frac{1}{2}}$. If in the boundary layer $h_K \sim \nu$, that is, a low local Reynolds number on the boundary, then the estimate is dominated by the $H^1(\Omega)$ contribution from the boundary that converges at optimal rate since the layer is resolved. The condition (2.7) is satisfied with $c\boldsymbol{\beta} \sim \nu^{\frac{1}{2}}$ showing that the strongest constraint on the mesh is not that of (2.7), but that of the $h_K^{-1}$ contribution on the boundary. In laminar free-flow we can expect $|\boldsymbol{\beta}|_{1,\infty,K} \le c\|\boldsymbol{\beta}\|_{0,\infty,K}$ to hold, and hence the convergence in the $L^2$-norm in this regime is of the quasi-optimal rate $h^{k+\frac{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 $(\mathbf{u}, p) \in [H^r(\Omega)]^d \times H^s(\Omega)$, with $r \ge 2$ and $s \ge 1$, is the solution of (2.1) and $(\mathbf{u}_h, p_h) \in W_h^k$ is the solution of (2.10). Then, under the assumptions of Lemma 4.6, there holds

$$ \|p - p_h\|_{0,\Omega} \le C \left( C_L \sigma^{\frac{1}{2}} + \max_{K \in \mathcal{T}_h} \{\|\boldsymbol{\beta}\|_{0,\infty,K} h_K, \nu\}^{\frac{1}{2}} + \sigma^{-\frac{1}{2}} \|\boldsymbol{\beta}\|_{0,\infty,\Omega} \right) C_{\mathbf{u}}, $$

with $C_{\mathbf{u}}$ the convergence rate of $\|\mathbf{u} - \mathbf{u}_h\|$ given by Theorem 4.9, and $C_L$ a positive constant depending only on $\Omega$.

*Proof.* Following [27, Corollary 2.4], there exists $\mathbf{v}_p \in [H_0^1(\Omega)]^d$ such that

$$ (4.13) \quad \nabla \cdot \mathbf{v}_p = p - p_h, \quad \|\mathbf{v}_p\|_{0,\Omega} \le C_L \|p - p_h\|_{0,\Omega}, \quad |\mathbf{v}_p|_{1,\Omega} \le C \|p - p_h\|_{0,\Omega}, $$

with $C_L > 0$ a constant, depending on $\Omega$, which scales as a distance. Thus, using the modified Galerkin orthogonality (Lemma 2.1), we readily obtain

$$ \begin{aligned} \|p - p_h\|_{0,\Omega}^2 &= (p - p_h, \nabla \cdot \mathbf{v}_p) \\ &= (p - p_h, \nabla \cdot (\mathbf{v}_p - \pi_{h,k} \mathbf{v}_p)) + \langle p - p_h, \pi_{h,k} \mathbf{v}_p \cdot \mathbf{n} \rangle_{\partial \Omega} \\ &\quad + a_h(\mathbf{u} - \mathbf{u}_h, \pi_{h,k} \mathbf{v}_p) + j_{\mathbf{u}}(\mathbf{u} - \mathbf{u}_h, \pi_{h,k} \mathbf{v}_p). \end{aligned} $$

Thus, after integrating by parts, we get

$$ (4.14) \quad \|p - p_h\|_{0,\Omega}^2 = \underbrace{(\nabla(p - p_h), \mathbf{v}_p - \pi_{h,k} \mathbf{v}_p)}_{T_1} + \underbrace{a_h(\mathbf{u} - \mathbf{u}_h, \pi_{h,k} \mathbf{v}_p) + j_{\mathbf{u}}(\mathbf{u} - \mathbf{u}_h, \pi_{h,k} \mathbf{v}_p)}_{T_2}. $$

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

$$ \begin{aligned} T_1 &= (\nabla(p - p_h) - \pi_{h,k} \nabla p + \pi_{h,k}^* \nabla p_h, \mathbf{v}_p - \pi_{h,k} \mathbf{v}_p) \\ &\le \|\tilde{h}(\nabla p - \pi_{h,k} \nabla p)\|_{0,\Omega} \|\tilde{h}^{-1}(\mathbf{v}_p - \pi_{h,k} \mathbf{v}_p)\|_{0,\Omega} \\ &\quad + C \|\phi_3^{-1}(\nabla p_h - \pi_{h,k}^* \nabla p_h)\|_{0,\Omega} \|\phi_3^*(\mathbf{v}_p - \pi_{h,k} \mathbf{v}_p)\|_{0,\Omega} \\ &\le C \left[ \left( \sum_{K \in \mathcal{T}_h} h_K^{2s_p-1} \|p\|_{s_p,K}^2 \right)^{\frac{1}{2}} \right. \\ &\qquad \left. + \max_{K \in \mathcal{T}_h} \{\nu, \|\boldsymbol{\beta}\|_{0,\infty,K} h_K\}^{\frac{1}{2}} j_p(p_h, p_h)^{\frac{1}{2}} \right] \|p - p_h\|_{0,\Omega}. \end{aligned} \quad (4.15) $$

Using the definition (2.11) of the bilinear form $a_h$, and after integration by parts in the convective term, we have

$$ (4.16) \quad T_2 \le \|\mathbf{u} - \mathbf{u}_h\|\|\pi_{h,k} \mathbf{v}_p\| + (\mathbf{u} - \mathbf{u}_h, \boldsymbol{\beta} \cdot \nabla \pi_{h,k} \mathbf{v}_p) \\ \qquad - \langle 2\nu\epsilon(\mathbf{u} - \mathbf{u}_h)\mathbf{n}, \pi_{h,k} \mathbf{v}_p \rangle_{\partial\Omega} - \langle \mathbf{u} - \mathbf{u}_h, 2\nu\epsilon(\pi_{h,k} \mathbf{v}_p)\mathbf{n} \rangle_{\partial\Omega}. $$

For the convective term we have, using the $H^1$-stability of the $L^2$-projection (see [5]) and (4.13),

$$ (4.17) \quad (\mathbf{u} - \mathbf{u}_h, \boldsymbol{\beta} \cdot \nabla \pi_{h,k} \mathbf{v}_p) \le \sigma^{-\frac{1}{2}} \|\boldsymbol{\beta}\|_{0,\infty,\Omega} \|\sigma^{\frac{1}{2}}(\mathbf{u} - \mathbf{u}_h)\|_{0,\Omega} \|\nabla \pi_{h,k} \mathbf{v}_p\|_{0,\Omega} \\ \le C \sigma^{-\frac{1}{2}} \|\boldsymbol{\beta}\|_{0,\infty,\Omega} \|\mathbf{u} - \mathbf{u}_h\| \|p - p_h\|_{0,\Omega}. $$
The boundary terms are controlled in the following fashion:

$$
(4.18) \quad \begin{aligned}
& \langle 2\nu\epsilon(\mathbf{u}-\mathbf{u}_h)\mathbf{n}, \pi_{h,k}\mathbf{v}_p \rangle_{\partial\Omega} + \langle \mathbf{u}-\mathbf{u}_h, 2\nu\epsilon(\pi_{h,k}\mathbf{v}_p)\mathbf{n} \rangle_{\partial\Omega} \\
& \le C\|(\nu\tilde{h})^{\frac{1}{2}}\epsilon(\mathbf{u}-\mathbf{u}_h)\|_{0,\partial\Omega}\|\pi_{h,k}\mathbf{v}_p\| + C\|(\nu\tilde{h})^{\frac{1}{2}}\epsilon(\pi_{h,k}\mathbf{v}_p)\|_{0,\partial\Omega}\|\mathbf{u}-\mathbf{u}_h\|.
\end{aligned}
$$

In addition, as in (4.7) and (4.8), we have

$$
(4.19) \qquad \begin{aligned}
\|(\nu\tilde{h})^{\frac{1}{2}}\epsilon(\mathbf{u}-\mathbf{u}_h)\|_{0,\partial\Omega} &\le \|(\nu\tilde{h})^{\frac{1}{2}}\epsilon(\mathbf{e}^\pi)\|_{0,\partial\Omega} + \|(\nu\tilde{h})^{\frac{1}{2}}\epsilon(\mathbf{e}_h)\|_{0,\partial\Omega}, \\
&\le C \left[ \left( \sum_{K \in \mathcal{T}_h} \nu h_K^{2(r_u-1)} \|\mathbf{u}\|_{r_u,K}^2 \right)^{\frac{1}{2}} + \|\mathbf{e}_h\| \right].
\end{aligned}
$$

In the same fashion, we obtain

$$
(4.20) \qquad \|(2\nu\tilde{h})^{\frac{1}{2}}\epsilon(\pi_{h,k}\mathbf{v}_p)\|_{0,\partial\Omega} \le C\|\pi_{h,k}\mathbf{v}_p\|_{0,\infty}.
$$

Finally, from (4.13), it follows that

$$
(4.21) \quad \| \pi_{h,k} \mathbf{v}_p \| \le C \left( C_L \sigma^{\frac{1}{2}} + \max_{K \in \mathcal{T}_h} \left\{ \| \boldsymbol{\beta} \|_{0,\infty,\kappa} h_K, \nu \right\}^{\frac{1}{2}} \right) \| p - p_h \|_{0,\Omega}.
$$

We conclude the proof by collecting the estimations (4.15)–(4.20) in (4.14), using
(4.21) and Theorem 4.9. $\square$

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.

THEOREM 4.14. Assume that the solution $(\mathbf{u}, p)$ of (2.1) belongs to $[H^2(\Omega)]^d \times H^1(\Omega)$ and let $(\mathbf{u}_h, p_h) \in W_h^k$ be the solution of (2.10). Assume also that

$$
(4.22) \qquad \|\boldsymbol{\beta}\|_{0,\infty,\kappa} h_K \le \nu \quad \forall \mathbf{K} \in \mathcal{T}_h,
$$

and that the solution $(\varphi, \psi)$ of the adjoint problem

$$
(4.23) \qquad \left\{
\begin{array}{ll}
\sigma\varphi - \boldsymbol{\beta} \cdot \nabla\varphi - 2\nu\nabla\cdot\boldsymbol{\epsilon}(\varphi) - \nabla\psi = \mathbf{u} - \mathbf{u}_h & \text{in } \Omega, \\
\nabla\cdot\varphi = 0 & \text{in } \Omega, \\
\varphi = 0 & \text{on } \partial\Omega
\end{array}
\right.
$$

belongs to $[H^2(\Omega)]^d \times [H^1(\Omega)]$ and satisfies

$$
(4.24) \qquad \|\varphi\|_{2,\Omega} + \|\psi\|_{1,\Omega} \le C \|\mathbf{u} - \mathbf{u}_h\|_{0,\Omega}.
$$

Then, there holds

$$
\| \mathbf{u} - \mathbf{u}_h \|_{0,\Omega} \le c h^2 ( \| \mathbf{u} \|_{2,\Omega} + \| p \|_{1,\Omega} ),
$$

with constant $c > 0$ independent of $h$, but depending on the physical parameters.
*Proof.* Multiplying the first equation of (4.23) by $\mathbf{u} - \mathbf{u}_h$ and the second by $-(p-p_h)$, integrating by parts, and using the modified Galerkin orthogonality (Lemma 2.1), it follows that

$$ \begin{aligned} \|\mathbf{u} - \mathbf{u}_h\|_{0,\Omega}^2 &= a_h(\mathbf{u} - \mathbf{u}_h, \varphi) + b_h(p - p_h, \varphi) - b_h(\psi, \mathbf{u} - \mathbf{u}_h) \\ &= \underbrace{a_h(\mathbf{u} - \mathbf{u}_h, \varphi - \pi_{h,k}\varphi) + b_h(p - p_h, \varphi - \pi_{h,k}\varphi) - b_h(\psi - \pi_{h,k}\psi, \mathbf{u} - \mathbf{u}_h)}_{T_1} \\ &\quad + \underbrace{j_u(\mathbf{u} - \mathbf{u}_h, \varphi - \pi_{h,k}\varphi)}_{T_2} + \underbrace{j_p(p_h, \pi_{h,k}\psi)}_{T_3}. \end{aligned} $$

Following the argument of the proofs of Theorems 4.9 and 4.12, and using Lemma 4.6 and (4.22), we get

$$ \begin{aligned} T_1 &\le \|\mathbf{u} - \mathbf{u}_h\|\|\varphi - \pi_{h,k}\varphi\| + |(\mathbf{u} - \mathbf{u}_h, \beta \cdot \nabla(\varphi - \pi_{h,k}\varphi))| \\ &\quad - \langle 2\nu\epsilon(\mathbf{u} - \mathbf{u}_h)\mathbf{n}, \varphi - \pi_{h,k}\varphi \rangle_{\partial\Omega} - \langle \mathbf{u} - \mathbf{u}_h, 2\nu\epsilon(\varphi - \pi_{h,k}\varphi)\mathbf{n} \rangle_{\partial\Omega} \\ &\le Ch (\|\mathbf{u} - \mathbf{u}_h\|\|\varphi\|_{2,\Omega} + \|p - p_h\|_{0,\Omega}\|\varphi\|_{2,\Omega} + \|\mathbf{u} - \mathbf{u}_h\|\|\psi\|_{1,\Omega}). \end{aligned} $$

Using Cauchy–Schwarz, Lemma 4.6, and (4.22), one obtains

$$ \begin{aligned} T_2 &\le j_u(\mathbf{u} - \mathbf{u}_h, \mathbf{u} - \mathbf{u}_h)^{\frac{1}{2}} j_u(\varphi - \pi_{h,k}\varphi, \varphi - \pi_{h,k}\varphi)^{\frac{1}{2}} \\ &\le Ch^{\frac{3}{2}} \|\mathbf{u} - \mathbf{u}_h\|\|\varphi\|_{2,\Omega}. \end{aligned} $$

Finally, from Lemma 4.7 and (4.22), for the last term we have

$$ \begin{aligned} T_3 &\le j_p(p_h, p_h)^{\frac{1}{2}} j_p(\pi_{h,k}\psi, \pi_{h,k}\psi)^{\frac{1}{2}} \\ &\le Ch j_p(p_h, p_h)^{\frac{1}{2}} \|\psi\|_{1,\Omega}. \end{aligned} $$

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+\frac{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(\Omega)$ and the velocity is in $[H^2(\Omega)]^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^{\frac{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 $\nu$-weighted $H^1$ contribution showing that the stabilization handles the numerical instability induced by treating nonsymmetric terms using the standard Galerkin method. In case $\sigma = 0$ the $H^1$ estimate obtained by a standard energy argument will scale as $\nu^{-\frac{1}{2}}$, reflecting the physical instability of the problem.

**5. 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.
<figure><img src="image_4.png" /><figcaption>FIG. 5.1. Convergence history: Linear elements (k = 1) (left) and quadratic elements (k = 2) (right).</figcaption></figure>

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]:

$$
\begin{equation}
\begin{aligned}
u_1(x_1, x_2, x_3) &= be^{a(x_1-x_3)+b(x_2-x_3)} - ae^{a(x_3-x_2)+b(x_1-x_2)}, \\
u_2(x_1, x_2, x_3) &= be^{a(x_2-x_1)+b(x_3-x_2)} - ae^{a(x_1-x_3)+b(x_2-x_3)}, \\
u_3(x_1, x_2, x_3) &= be^{a(x_3-x_2)+b(x_1-x_2)} - ae^{a(x_2-x_1)+b(x_3-x_1)}, \\
(5.1) \qquad \qquad \qquad \qquad \qquad \qquad p(x_1, x_2, x_3) &= (a^2 + b^2 + ab) \left[ e^{a(x_1-x_2)+b(x_1-x_3)} + e^{a(x_2-x_3)+b(x_2-x_1)} \right. \\
& \qquad \qquad \qquad \qquad \qquad \left. + e^{a(x_3-x_1)+b(x_3-x_2)} \right]
\end{aligned}
\end{equation}
$$

with $\beta = \mathbf{u}$, $\sigma = 1$, $\nu = 10^{-4}$, $a = b = 0.75$, and $\Omega = (0, 1)^3$ the unit cube.

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 three-dimensional research code [23]). The stabilization parameter involved in the jumps terms (2.13) and (2.14) were chosen as

$$
\gamma = \begin{cases} \frac{1}{8} & \text{if } k = 1, \\ \frac{1}{32} & \text{if } k = 2. \end{cases}
$$

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 $\mathbf{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).
<figure><img src="image_5.png" /><figcaption>FIG. 5.2. Pressure contours: Coarse mesh (left) and fine mesh (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5.3. Coarse mesh (2929 tetrahedra) and fine mesh (196608 tetrahedra).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5.4. Velocity ($u_{h2}$) contours on a cutting plane: Stabilized (left), with $\gamma\beta = 0$ (right).</figcaption></figure>

In what follows we will replace, in (5.1), the expression for the pressure by

$$p(x_1, x_2, x_3) = \begin{cases} 2x_2 & \text{if } 0 \le x_2 \le \frac{1}{2}, \\ 2(1 - x_2) & \text{if } \frac{1}{2} \le x_2 \le 1. \end{cases}$$

Clearly, this function satisfies $p \in H^1(\Omega)$ but does not belong to $H^2(\Omega)$. Figure 5.5
<figure><img src="image_8.png" /><figcaption>FIG. 5.5. Cutting plane pressure: Coarse mesh and fine mesh.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5.6. Convergence history: Linear elements, nonsmooth pressure, stabilization, parameters chosen as in (2.13) and (2.14).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 5.7. Convergence history: Linear elements, nonsmooth pressure, parameters chosen as in (2.13) and (2.14). Left: High Reynolds number. Right: Low Reynolds number.</figcaption></figure>

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^{\frac{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^{\frac{3}{2}})$ convergence order. As expected, when the local Reynolds is low (for instance $\nu = 0.1$), we recover the optimal $O(h)$; see Figure 5.7,
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)$.

**6. 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 \in H^1(\Omega)$) and the local Reynolds number is high, the convergence may become suboptimal $O(h^{\frac{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+\frac{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^{\frac{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

$$ (6.1) \qquad (\nu_T(I - P_H)\epsilon(\mathbf{u}_h), (I - P_H)\epsilon(\mathbf{v}_h)), $$

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 $P_H \stackrel{\text{def}}{=} \pi_{h,k}^*$,

$$ \|\nu_T^{\frac{1}{2}}(I - \pi_{h,k}^*)\epsilon(\mathbf{u}_h)\|_{\Omega}^2 \le \sum_{K \in \mathcal{T}_h} \int_{\partial K} \nu_T h_K [[\epsilon(\mathbf{u}_h)]] : [[\epsilon(\mathbf{u}_h)]] \, ds, $$

and we conclude that a possible subgrid modeling term would be

$$ j_T(\mathbf{u}_h, \mathbf{v}_h) = \sum_{K \in \mathcal{T}_h} \int_{\partial K \setminus \partial \Omega} \nu_T h_K [[\epsilon(\mathbf{u}_h)]] : [[\epsilon(\mathbf{v}_h)]] \, ds, $$

where the choice of $\nu_T$ now is a modeling issue. It should be noted that the choice $\nu_T = \gamma h_K$ gives us a term which is asymptotically equivalent to the face penalty operator using the whole gradient. However, other choices of $\nu_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 no-slip conditions when the boundary layers are resolved; this may be favorable in LES (see Layton [34]).

**Acknowledgment.** The authors wish to thank the anonymous reviewers for their careful reading of the manuscript and their many constructive comments.

## REFERENCES

[1] R. BECKER AND M. BRAACK, *A finite element pressure gradient stabilization for the Stokes equations based on local projections*, Calcolo, 38 (2001), pp. 173–199.

[2] P. BOCHEV AND M. GUNZBURGER, *An absolutely stable pressure-Poisson stabilized finite element method for the Stokes equations*, SIAM J. Numer. Anal., 42 (2004), pp. 1189–1207.

[3] M. BOMAN, *Estimates for the $L_2$-projection onto continuous finite element spaces in a weighted $L_p$-norm*, BIT, 46 (2006), to appear.

[4] M. BRAACK AND E. BURMAN, *Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method*, SIAM J. Numer. Anal., 43 (2006), pp. 2544–2566.

[5] J. BRAMBLE, J. PASCIAK, AND O. STEINBACH, *On the stability of the $L^2$ projection in $H^1(\Omega)$*, Math. Comp., 71 (2002), pp. 147–156.

[6] P. G. CIARLET, *Mathematical Elasticity*, Volume I: Three-Dimensional Elasticity, North-Holland, Amsterdam, 1988.

[7] F. BREZZI AND M. FORTIN, *Mixed and Hybrid Finite Element Methods*, Springer Ser. Comput. Math. 15, Springer-Verlag, New York, 1991.

[8] F. BREZZI AND M. FORTIN, *A minimal stabilisation procedure for mixed finite element methods*, Numer. Math., 89 (2001), pp. 457–491.

[9] A. N. BROOKS AND T. J. R. HUGHES, *Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations*, Comput. Methods Appl. Mech. Engrg., 32 (1982), pp. 199–259.

[10] E. BURMAN, *A unified analysis for conforming and nonconforming stabilized finite element methods using interior penalty*, SIAM J. Numer. Anal., 43 (2005), pp. 2012–2033.

[11] E. BURMAN AND M. A. FERNÁNDEZ, *Stabilized finite element schemes for incompressible flow using velocity/pressure spaces satisfying the LBB-condition*, in Proceedings of the 6th World Congress in Computational Mechanics (WCCM VI), Beijing, China, 2004.

[12] E. BURMAN, M. A. FERNÁNDEZ, AND P. HANSBO, *Edge stabilization: An interior penalty method for the incompressible Navier-Stokes equation*, in Proceedings of the Fourth European Congress on Computational Methods in Applied Sciences and Engineering, Vol. I, Jyväskylä, Finland, 2004.

[13] E. BURMAN AND P. HANSBO, *Edge stabilization for the generalized Stokes problem: A continuous interior penalty method*, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 2393–2410.

[14] E. BURMAN AND P. HANSBO, *A stabilized nonconforming finite element method for incompressible flow*, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 2881–2899.

[15] E. BURMAN AND P. HANSBO, *Edge stabilization for Galerkin approximations of convection-diffusion-reaction problems*, Comput. Methods Appl. Mech. Engrg., 193 (2004), pp. 1437–1453.

[16] P. CLÉMENT, *Approximation by finite element functions using local regularization*, Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge Anal. Numér., 9 (1975), pp. 77–84.

[17] R. CODINA, *Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales*, Publication CIMNE 289, International Center for Numerical Methods
in Engineering (CIME), Barcelona, Spain, 2006.

[18] R. CODINA AND J. BLASCO, *Analysis of a pressure-stabilized finite element approximation of the stationary Navier-Stokes equations*, Numer. Math., 87 (2000), pp. 59-81.

[19] C. R. DOHRMANN AND P. B. BOCHEV, *A stabilized finite element method for the Stokes problem based on polynomial pressure projections*, Internat. J. Numer. Methods Fluids, 46 (2004), pp. 183-201.

[20] J. DOUGLAS AND T. DUPONT, *Interior penalty procedures for elliptic and parabolic Galerkin methods*, in Computing Methods in Applied Sciences (Second International Symposium, Versailles, 1975), Lecture Notes in Phys. 58, Springer-Verlag, Berlin, 1976, pp. 207-216.

[21] A. ERN AND J.-L. GUERMOND, *Theory and Practice of Finite Elements*, Appl. Math. Sci. 159, Springer-Verlag, New York, 2004.

[22] C. ETHIER AND D. STEINMAN, *Exact fully 3-D Navier-Stokes solutions for benchmarking*, Internat. J. Numer. Methods Fluids, 19 (1994), pp. 369-375.

[23] M. A. FERNÁNDEZ, L. FORMAGGIA, A. GAUTHIER, J. GERBEAU, C. PRUD'HOMME, AND A. VENEZIANI, *The LifeV Project*, http://www.lifev.org.

[24] L. P. FRANCA AND S. L. FREY, *Stabilized finite element methods. II. The incompressible Navier-Stokes equations*, Comput. Methods Appl. Mech. Engrg., 99 (1992), pp. 209-233.

[25] J. FREUND AND R. STENBERG, *On weakly imposed boundary conditions for second order problems*, in Proceedings of the Ninth International Conference on Finite Elements in Fluids, Venice, Italy, 1995, pp. 327-336.

[26] T. GELHARD, G. LUBE, M. A. OLSHANSKII, AND J.-H. STARCKE, *Stabilized finite element schemes with LBB-stable elements for incompressible flows*, J. Comput. Appl. Math., 177 (2005), pp. 243-267.

[27] V. GIRAULT AND P. A. RAVIART, *Finite Element Methods for Navier-Stokes Equations*, Springer Series in Comput. Math. 5, Springer-Verlag, Berlin, 1986.

[28] J. L. GUERMOND, *Stabilization of Galerkin approximations of transport equations by subgrid modeling*, M2AN Math. Model. Numer. Anal., 33 (1999), pp. 1293-1316.

[29] P. HANSBO AND A. SZEPESSY, *A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations*, Comput. Methods Appl. Mech. Engrg., 84 (1990), pp. 175-192.

[30] T. HUGHES, L. MAZZEI, AND K. JANSEN, *Large eddy simulation and the variational multiscale method*, Comput. Vis. Sci., 3 (2000), pp. 47-59.

[31] V. JOHN AND S. KAYA, *A finite element variational multiscale method for the Navier-Stokes equations*, SIAM J. Sci. Comput., 26 (2005), pp. 1485-1503.

[32] C. JOHNSON AND J. SARANEN, *Streamline diffusion methods for the incompressible Euler and Navier-Stokes equations*, Math. Comp., 47 (1986), pp. 1-18.

[33] O. KARAKASHIAN AND F. PASCAL, *A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems*, SIAM J. Numer. Anal., 41 (2003), pp. 2374-2399.

[34] W. LAYTON, *Weak imposition of “no-slip” conditions in finite element methods*, Comput. Math. Appl., 38 (1999), pp. 129-142.

[35] W. LAYTON, *A connection between subgrid scale eddy viscosity and mixed methods*, Appl. Math. Comput., 133 (2002), pp. 147-157.

[36] J. NITSCHE, *Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind*, Abh. Math. Sem. Univ. Hamburg, 36 (1971), pp. 9-15.

[37] V. THOMÉE, *Galerkin Finite Element Methods for Parabolic Problems*, Springer Series in Comput. Math. 25, Springer-Verlag, Berlin, 1997.

[38] L. TOBISKA AND R. VERFÜRTH, *Analysis of a streamline diffusion finite element method for the Stokes and Navier-Stokes equations*, SIAM J. Numer. Anal., 33 (1996), pp. 107-127.