A locally conservative LDG method for the incompressible Navier-Stokes equations

Abstract

In this paper a new local discontinuous Galerkin method for the incompressible stationary Navier-Stokes equations is proposed and analyzed. Four important features render this method unique: its stability, its local conservativity, its high-order accuracy, and the exact satisfaction of the incompressibility constraint. Although the method uses completely discontinuous approximations, a globally divergence-free approximate velocity in $H(\mathrm{div};\Omega )$ is obtained by simple, element-by-element post-processing. Optimal error estimates are proven and an iterative procedure used to compute the approximate solution is shown to converge. This procedure is nothing but a discrete version of the classical fixed point iteration used to obtain existence and uniqueness of solutions to the incompressible Navier-Stokes equations by solving a sequence of Oseen problems. Numerical results are shown which verify the theoretical rates of convergence. They also confirm the independence of the number of fixed point iterations with respect to the discretization parameters. Finally, they show that the method works well for a wide range of Reynolds numbers.

1. Introduction

In this paper we propose and analyze a local discontinuous Galerkin (LDG) method for the stationary incompressible Navier-Stokes equations

$$\begin{alignat}{2} -\nu \Delta \mathbf{u}+\nabla \cdot (\mathbf{u}\otimes \mathbf{u})+\nabla p&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},{\tag{1.1}}\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma =\partial \Omega $}. \end{alignat}$$

Here $\nu$ is the kinematic viscosity, $\mathbf{u}$ the velocity, $p$ the pressure, and $\mathbf{f}$ the external body force. For the sake of simplicity, we take $\Omega$ to be a polygonal domain of $\mathbb{R}^2$.

This paper is the fourth in a series (9, 8, and 7) devoted to the study of the LDG method as applied to incompressible fluid flow problems. In 9, we considered the Stokes equations

$$\begin{alignat}{2} -\nu \Delta \mathbf{u}+\nabla p&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},{\tag{1.2}}\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma $}, \end{alignat}$$

and focused on the problem of how to deal with the incompressibility condition. Later, in 8, we considered the Oseen equations

$$\begin{alignat}{2} -\nu \Delta \mathbf{u}+(\mathbf{w}\cdot \nabla )\,\mathbf{u}+\nabla p&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},{\tag{1.3}}\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma $}, \end{alignat}$$

where the convective velocity $\mathbf{w}$ was taken to be a smooth function, and focused on the problem of how to incorporate the linear convective term. The resulting method was shown to be optimally convergent and robust for a wide range of Reynolds numbers. A succinct review of this work can be found in 7.

In this paper, we continue our study of LDG methods for incompressible flows and consider their application to the Navier-Stokes equations 1.1. Our main concern is to devise LDG methods that are locally conservative and can be proven to be stable. The local conservativity, a property highly valued by practitioners of computational fluid dynamics, is a discrete version of the identity

$$\begin{equation} \int _{\partial K} \left( -\nu \, \nabla \mathbf{u}\cdot \mathbf{n}_K+ (\mathbf{u}\cdot \mathbf{n}_K)\,\mathbf{u} + p\,\mathbf{n}_K \right)\,ds = \int _K \mathbf{f}\,d\mathbf{x}, {\tag{1.4}} \end{equation}$$

where $K$ is an arbitrary subdomain of $\Omega$ with outward normal unit vector $\mathbf{n}_K$. This property can be easily enforced by LDG methods as soon as the equations are written in divergence form. However, to devise an LDG method that can be proven to be stable (that is, that satisfies a discrete version of the stability estimate for the continuous case,

$$\begin{equation} \|\mathbf{u}\|_1 \leq \frac{ C_P \|\mathbf{f}\|_0}{\nu }, {\tag{1.5}} \end{equation}$$

where the constant $C_P$ is the Poincaré constant) is extremely difficult.

The reason for this is that, in order to obtain the stability estimate 1.5, the incompressibility condition must be used. It is well known that, for many numerical methods for the Stokes and the Oseen equations, a weakly enforced incompressibility is enough to guarantee stability. However, this is not so for the Navier-Stokes equations because of the presence of nonlinear convection. Moreover, the now standard solution to this problem, 21, 22, which is based on a suitable modification of the nonlinearity of the Navier-Stokes equations, cannot be used. This happens because such a modification does not have divergence form and hence prevents LDG methods from being locally conservative.

The main contribution of this paper is to show how to overcome this difficulty. In fact, we show that this can be done in two ways. The first one focuses on the convective nonlinearity and is based on a new modification in divergence form of the nonlinearity; it will be explored elsewhere in detail. The second one, which constitutes the main subject of this paper, focuses on the incompressibility constraint and is based on discretizing the Oseen equations 1.3 where the convective velocity $\mathbf{w}$ is taken to be a projection of the approximate velocity $\mathbf{u}_h$,

$$\begin{equation} \mathbf{w}=\mathbb{P} \mathbf{u}_h, {\tag{1.6}} \end{equation}$$

into the space of globally divergence-free functions. This projection is a slight modification of well-known projections $\Pi _h$ with the property

$$\begin{equation*} P_h \nabla \cdot =\nabla \cdot \Pi _h, \end{equation*}$$

where $P_h$ is an $L^2$-projection; see 5. Its implementation is very efficient as it can be computed in an element-by-element fashion.

Thus, given the convective velocity $\mathbf{w}=\mathbb{P}\mathbf{u}_h$, the resulting scheme is nothing but the LDG method 8 already studied for the Oseen equations 1.3. Since that method is stable, high-order accurate and locally conservative, so is the LDG method under consideration. Moreover, the approximation to the velocity given by $\mathbf{w}$ has continuous normal components across elements and is globally divergence free in $H(\operatorname {div};\Omega )=\{\mathbf{v}\in L^2(\Omega )^2: \nabla \cdot \mathbf{v}\in L^2(\Omega )\}$. To the knowledge of the authors, no other numerical scheme for the incompressible Navier-Stokes equations has all these properties.

Of course, the convective velocity $\mathbf{w}$ depends on the discrete velocity field $\mathbf{u}_h$ through 1.6 and, hence, we need to use an iterative method to compute it. To do so, we note that if $S({\overline{\mathbf{u}}}_h)$ is the LDG approximate velocity of the Oseen problem with convective velocity $\mathbf{w}=\mathbb{P} {\overline{\mathbf{u}}}_h$, then the approximate velocity $\mathbf{u}_h$ of the LDG method under consideration is a fixed point of $S$. If $S$ is proven to be a contraction, to compute the approximate solution $\mathbf{u}_h$, we can use the fixed point iteration

$$\begin{equation*} \mathbf{u}_h^{\ell +1} := S({\mathbf{u}}_h^\ell ). \end{equation*}$$

This is nothing but a discrete version of the argument used to prove the existence and uniqueness of the exact solution of 1.1. It ensures the existence and uniqueness of the exact solution $(\mathbf{u},p)\in H^1_0(\Omega )^2\times L^2(\Omega )/\mathbb{R}$ of 1.1 under a smallness condition of the type

$$\begin{equation} \frac{ C_\Omega C_P \|\mathbf{f}\|_0}{\nu ^2} < 1, {\tag{1.7}} \end{equation}$$

where $C_\Omega >0$ only depends on $\Omega$; see 17, Theorem 10.1.1 and the references therein. We mimic this argument to show that the approximate solution of the LDG method exists and is unique under a similar condition.

Let us point out that exact incompressibility can be achieved trivially if the LDG method has a velocity space that is div-conforming, i.e., it is included in $H(\operatorname {div};\Omega )$. In this particular case, weak incompressibility implies exact incompressibility, provided that the discrete spaces are matched correctly. As will be discussed below, this approach can be viewed as a particular LDG method for which the operator ${\mathbb{P}}$ is chosen to be the identity. Consequently, all the results of this paper hold true verbatim for methods that are based on div-conforming velocity spaces. Furthermore, we note that, although we have used the LDG method to discretize the terms associated with the viscosity effects, any other DG discretization whose primal form is both coercive and continuous could have been used to that effect; see the discussions in 2 and 19.

The organization of the paper is as follows. In Section 2 we discuss the ideas that motivate the devising of the LDG method we propose in this paper. In Section 3 we present the LDG discretization in detail and verify its local conservativity. In Section 4 we state and discuss the main results, namely, the stability of the method, the convergence of the fixed point iteration, and the a priori error estimates, and in Section 5 we present their proofs. In Section 6 we present numerical experiments verifying the theoretical results. We end in Section 7 with some concluding remarks.

2. Devising the LDG method

In this section we discuss the ideas that led us to the devising of an LDG method that is both stable and locally conservative. To keep the discussion as simple and clear as possible, we do not work with the numerical method. Instead, we work directly with the equations 1.1 and infer, from their structure, the properties of the corresponding LDG method.

2.1. A locally conservative LDG method

Since the incompressible Navier-Stokes equations 1.1 are written in divergence form, a locally conservative LDG method can be easily constructed. However, it is very difficult to prove its stability. Let us illustrate this difficulty by using the equations for the exact solution.

If we multiply the first equation of 1.1 by $\mathbf{u}$, integrate by parts and use the boundary conditions, we get

$$\begin{equation*} \nu \,\int _{\Omega } \nabla \mathbf{u}: \nabla \mathbf{u}\, d\mathbf{x} +\frac{1}{2} \int _{\Omega } |\mathbf{u}\,|^2\, \nabla \cdot \mathbf{u}\, d\mathbf{x} - \int _{\Omega } p\, \nabla \cdot \mathbf{u}\, d\mathbf{x} = \int _{\Omega } \mathbf{f}\cdot \mathbf{u}\, d\mathbf{x}. \end{equation*}$$

We see that we must use the incompressibility condition to obtain the equation

$$\begin{equation*} \nu \,\int _{\Omega } \nabla \mathbf{u}: \nabla \mathbf{u}\, d\mathbf{x} = \int _{\Omega } \mathbf{f}\cdot \mathbf{u}\, d\mathbf{x}, \end{equation*}$$

from which the stability estimate 1.5 immediately follows.

In general, since exact incompressibility is very difficult to achieve after discretization, it is usually only enforced weakly. This weak incompressibility is enough, in a wide variety of cases, to guarantee that the discrete version of the term

$$\begin{equation*} \int _{\Omega } p\, \nabla \cdot \mathbf{u}\, d\mathbf{x} \end{equation*}$$

is exactly zero, as for most mixed methods for the Stokes and Navier-Stokes equations, or nonnegative, as for the LDG methods considered for the Stokes 9 and Oseen 8 problems. Unfortunately, this is not true for the discrete version of the term

$$\begin{equation*} \frac{1}{2} \int _{\Omega } |\,\mathbf{u}\,|^2\, \nabla \cdot \mathbf{u}\, d\mathbf{x}, \end{equation*}$$

because the square of the modulus of the approximate velocity does not necessarily belong to the space of the approximate pressure.

2.2. The classical modification of the nonlinearity

A solution to this impasse can be obtained by using a now classical technique proposed back in the 1960s; see 21 and 22. From our perspective, it consists in modifying the nonlinearity of the equations as follows:

$$\begin{alignat*}{2} -\nu \Delta \mathbf{u}+\nabla \cdot (\mathbf{u}\otimes \mathbf{u}) -\frac{1}{2}(\nabla \cdot \mathbf{u})\,\mathbf{u}+ \nabla p&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma $}. \end{alignat*}$$

To see that this solves the problem, multiply the first equation by $\mathbf{u}$, integrate by parts, and use the boundary conditions to get

$$\begin{equation*} \nu \,\int _{\Omega } \nabla \mathbf{u}: \nabla \mathbf{u}\, d\mathbf{x} - \int _{\Omega } p\, \nabla \cdot \mathbf{u}\, d\mathbf{x} = \int _{\Omega } \mathbf{f}\cdot \mathbf{u}\, d\mathbf{x}. \end{equation*}$$

As a consequence, stability can follow from weak incompressibility. This is the case for most mixed methods for the Navier-Stokes equations. It is also the case for the first discontinuous Galerkin method for the incompressible Navier-Stokes equations 13, a method which uses locally divergence-free polynomial approximations of the velocity, and for the more recent discontinuous Galerkin method developed in 10.

The only problem with this approach is that local conservativity cannot be achieved because the first equation is not written in divergence form.

2.3. A new modification of the nonlinearity

To overcome this difficulty, it is enough to take another glance at the first equation in this section to realize that instead of working with the kinematic pressure $p$, we should work with the new variable

$$\begin{equation*} P=p-\frac{1}{2}\,|\,\mathbf{u}\,|^2. \end{equation*}$$

If we incorporate this unorthodox pressure into the Navier-Stokes equations, we get

$$\begin{alignat*}{2} -\nu \Delta \mathbf{u}+\nabla \cdot (\mathbf{u}\otimes \mathbf{u}) +\frac{1}{2}\nabla \,|\,\mathbf{u}\,|^2+ \nabla P&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma $}. \end{alignat*}$$

We now see that, since the above modification is in divergence form, locally conservative LDG methods can easily be constructed. Moreover, since we have

$$\begin{equation*} \nu \,\int _{\Omega } \nabla \mathbf{u}: \nabla \mathbf{u}\, d\mathbf{x} - \int _{\Omega } P\, \nabla \cdot \mathbf{u}\, d\mathbf{x} = \int _{\Omega } \mathbf{f}\cdot \mathbf{u}\, d\mathbf{x}, \end{equation*}$$

we also see that the stability of the LDG method can follow from weak incompressibility. The LDG method obtained with this approach can indeed be proven to have those properties; it is going to be studied thoroughly in a forthcoming paper.

2.4. Enforcing exact incompressibility

As we have seen, suitable modifications of the nonlinearity can be introduced which allow us to obtain stability by using only weakly incompressible approximations to the velocity. The approach we consider in this paper does not rely on a modification of that type. Instead, it is based on enforcing exact incompressibility in the space

$$H(\operatorname {div};\Omega ):=\{ \mathbf{v}\in L^2(\Omega )^2: \nabla \cdot \mathbf{v}\in L^2(\Omega )\}.$$

The idea that allows this to happen is based on two observations. The first is that we can rewrite the Navier-Stokes equations as the Oseen problem

$$\begin{alignat*}{2} -\nu \Delta \mathbf{u}+(\mathbf{w}\cdot \nabla )\,\mathbf{u} + \nabla p&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma $}, \end{alignat*}$$

where, of course, $\mathbf{w}=\mathbf{u}$. If we multiply the first equation by $\mathbf{u}$, integrate by parts, and use the boundary conditions, we get

$$\begin{equation*} \nu \,\int _{\Omega } \nabla \mathbf{u}: \nabla \mathbf{u}\, d\mathbf{x} -\frac{1}{2} \int _{\Omega } |\,\mathbf{u}\,|^2\, \nabla \cdot \mathbf{w}\, d\mathbf{x} - \int _{\Omega } p\, \nabla \cdot \mathbf{u}\, d\mathbf{x} = \int _{\Omega } \mathbf{f}\cdot \mathbf{u}\, d\mathbf{x}. \end{equation*}$$

This suggests that we consider an LDG method with two different (but strongly related) approximations to the velocity: one approximation for $\mathbf{u}$ and another for $\mathbf{w}$. The stability for the LDG method would then be achieved if the approximation to $\mathbf{u}$ is weakly incompressible and if the approximation to $\mathbf{w}$ is exactly incompressible. Further, local conservativity can be readily achieved for such an LDG method. Indeed, the fact that the equations are not written in conservative form can be compensated for by the fact that the approximation to $\mathbf{w}$ is globally divergence free.

The second observation is that, if the approximation to $\mathbf{u}$ given by an LDG method $\mathbf{u}_h$ is weakly incompressible, it is possible to compute, in an element-by-element fashion, another approximation, $\mathbf{w}=\mathbb{P}(\mathbf{u}_h)$, which is exactly incompressible. As we shall see, the post-processing operator $\mathbb{P}$ is a slight modification of the well-known Brezzi-Douglas-Marini (BDM) interpolation operator; see 4.

We note that the post-processing procedure can be omitted in the particular case where the velocity space is div-conforming. Indeed, in this case the fact that the approximation to the velocity is weakly incompressible does imply that it is exactly incompressible, provided the pressure space is chosen suitably. Hence, we can take $\mathbf{w}=\mathbf{u}_h$.

We are now ready to describe the LDG method in detail.

3. The LDG method

In this section, we introduce a locally conservative LDG discretization for the Navier-Stokes equations 1.1.

3.1. Meshes and trace operators

We begin by introducing some notation. We denote by ${\mathcal{T}}_h$ a regular and shape-regular triangulation of mesh-size $h$ of the domain $\Omega$ into triangles $\{K\}$. We further denote by ${\mathcal{E}}_h^{\mathcal{I}}$ the set of all interior edges of ${\mathcal{T}}_h$ and by ${\mathcal{E}}_h^{\mathcal{B}}$ the set of all boundary edges. We set ${\mathcal{E}}_h={\mathcal{E}}_h^{\mathcal{I}} \cup {\mathcal{E}}_h^{\mathcal{B}}$.

Next we introduce notation associated with traces. Let $K^+$ and $K^-$ be two adjacent elements of ${\mathcal{T}}_h$. Let $\mathbf{x}$ be an arbitrary point of the interior edge $e=\partial K^+\cap \partial K^-\in {\mathcal{E}}_h^{\mathcal{I}}$. Let $\varphi$ be a piecewise smooth scalar-, vector-, or matrix-valued function and let us denote by $\varphi ^\pm$ the traces of $\varphi$ on $e$ taken from within the interior of $K^\pm$. Then, we define the mean value $\{\!\!\{\cdot \}\!\!\}$ at $\mathbf{x}\in e$ as

$$\{\!\!\{\varphi \}\!\!\}:=\frac{1}{2}(\varphi ^++\varphi ^-).$$

Further, for a generic multiplication operator $\odot$, we define the jump $[\![\cdot ]\!]$ at $\mathbf{x}\in e$ as

$$[\![\varphi \odot \mathbf{n}]\!]:=\varphi ^+\odot \mathbf{n}_{K^+} + \varphi ^-\odot \mathbf{n}_{K^-}.$$

Here, $\mathbf{n}_K$ denotes an outward unit normal vector on the boundary $\partial K$ of element $K$. On boundary edges, we set accordingly $\{\!\!\{\varphi \}\!\!\}:=\varphi$, and $[\![\varphi \odot \mathbf{n}]\!]:=\varphi \odot \mathbf{n}$, with $\mathbf{n}$ denoting the outward unit normal vector on $\Gamma$.

3.2. The LDG method for the Oseen equations

We now recall the LDG method for the Oseen equations 1.3. We assume that the convective velocity field $\mathbf{w}$ is in the space

$$\begin{equation*} \mathbf{J}({\mathcal{T}}_h)=\{ \mathbf{v}\in L^2(\Omega )^2: \nabla \cdot \mathbf{v}\equiv 0 \text{ and }\mathbf{v}|_K \in H^1(K)^2,\ K\in {\mathcal{T}}_h \}. \end{equation*}$$

We begin by introducing the auxiliary variable $\underline{\sigma }=\nu \nabla \mathbf{u}$ and rewriting the Oseen equations as

$$\begin{alignat}{2} \underline{\sigma }& =\nu \nabla \mathbf{u} & \qquad \qquad &\text{in $\Omega $},\\ -\nabla \cdot \underline{\sigma }+(\mathbf{w}\cdot \nabla )\mathbf{u}+\nabla p&=\mathbf{f} & \qquad \qquad &\text{in $\Omega $},{\tag{3.1}}\\ \nabla \cdot \mathbf{u}&= 0 & &\text{in $\Omega $},\\ \mathbf{u}&=\mathbf{0} & &\text{on $\Gamma $}. \end{alignat}$$

Next, we introduce the space $\underline{\Sigma }_h\times \mathbf{V}_h\times Q_h$ where

$$\begin{equation*} \begin{split} \underline{\Sigma }_h &=\{\, \mathbf{v}\in L^2(\Omega )^{2\times 2} : \underline{\tau }|_K \in {P}_{k}(K)^{2\times 2}, \ K\in {\mathcal{T}}_h \,\},\\ \mathbf{V}_h&=\{\, \mathbf{v}\in L^2(\Omega )^2 : \mathbf{v}|_K \in {P}_{k}(K)^2, \ K\in {\mathcal{T}}_h \,\},\\ Q_h&=\{\, q\in L^2(\Omega ): q|_K \in {P}_{k-1}(K),\ K\in {\mathcal{T}}_h,\ {\textstyle \int _\Omega \, q\,d\mathbf{x}=0} \, \}, \end{split} \end{equation*}$$

for an approximation order $k\geq 1$. Here ${P}_{k}(K)$ denotes the space of polynomials of total degree at most $k$ on $K$. For simplicity we consider here only so-called mixed-order elements where the approximation degree in the pressure is of one order lower than the one in the velocity.

Finally, we define the approximate solution $(\underline{\sigma }_h,\mathbf{u},p_h)\in \underline{\Sigma }_h\times \mathbf{V}_h\times Q_h$ by requesting that for each $K\in {\mathcal{T}}_h$,

$$\begin{alignat}{1} &\int _K\, \underline{\sigma }_{h}:\underline{\tau }\, d\mathbf{x} = -\nu \int _K\, \mathbf{u}_{h}\cdot \nabla \cdot \underline{\tau }\, d\mathbf{x} +\nu \int _{\partial K}\, \widehat{\mathbf{u}}_{h}^\sigma \cdot \underline{\tau }\cdot \mathbf{n}_K\, ds,{\tag{3.2}}\\ &\int _K\, \big [\underline{\sigma }_{h}: \nabla \mathbf{v} - p_h \, \nabla \cdot \mathbf{v}\big ] \, d\mathbf{x} -\int _{\partial K}\, \big [ \widehat{\underline{\sigma }}_{h}:(\mathbf{v}\otimes \mathbf{n}_K) - \widehat{p}_h \, \mathbf{v}\cdot \mathbf{n}_{K} \big ] \, ds{\tag{3.3}}\\ & -\int _K\, \, \mathbf{u}_{h} \cdot \nabla \cdot (\mathbf{v}\otimes \mathbf{w})\, d\mathbf{x} +\int _{\partial K}\, \mathbf{w}\cdot \mathbf{n}_{K} \, \widehat{\mathbf{u}}^{\mathbf{w}}_{h}\cdot \mathbf{v} \, ds = \int _K \mathbf{f}\cdot \mathbf{v} \, d\mathbf{x}, \\ &-\int _K\, \mathbf{u}_h\cdot \nabla q\, d\mathbf{x}+ \int _{\partial K}\, \widehat{\mathbf{u}}^p_h\cdot \mathbf{n}_K q\, ds =0, {\tag{3.4}} \end{alignat}$$

for all test functions $(\underline{\tau },\mathbf{v},q)\in \underline{\Sigma }_h\times \mathbf{V}_h\times Q_h$. Each of the above equations is enforced locally, that is, element by element, due to the appearance of the so-called numerical fluxes $\widehat{\mathbf{u}}^\sigma _{h}$, $\widehat{\underline{\sigma }}_{h}$, $\widehat{p}_h$, $\widehat{\mathbf{u}}^{\mathbf{w}}_{h}$, and $\widehat{\mathbf{u}}^p_{h}$.

Thanks to this structure of the LDG method we immediately get that

$$\begin{equation*} \int _{\partial K} \left( -\widehat{\underline{\sigma }}_h\cdot \mathbf{n}_K +(\mathbf{w}\cdot \mathbf{n}_K)\,\widehat{\mathbf{u}}_h^{\mathbf{w}} + \widehat{p}_h\,\mathbf{n}_K \right)\,ds -\int _{K}\,\mathbf{u}_h\,\nabla \cdot {\mathbf{w}} \,d\mathbf{x} = \int _K \mathbf{f}\,d\mathbf{x}, \end{equation*}$$

and since $\mathbf{w}$ is globally divergence free, we obtain a discrete version of the property of local conservativity 1.4, namely,

$$\begin{equation*} \int _{\partial K} \left( -\widehat{\underline{\sigma }}_h\cdot \mathbf{n}_K +(\mathbf{w}\cdot \mathbf{n}_K)\,\widehat{\mathbf{u}}_h^{\mathbf{w}} + \widehat{p}_h\,\mathbf{n}_K \right)\,ds = \int _K \mathbf{f}\,d\mathbf{x}. \end{equation*}$$

In other words, the LDG method is locally conservative.

To ensure that the method is also stable (and high-order accurate), the numerical fluxes, which are nothing but discrete approximations to the traces on the boundary of the elements, must be defined carefully. As we shall prove, the numerical fluxes that define the LDG method for the Oseen equations 8 do ensure stability. For the sake of clarity we consider the fluxes in their simplest form.

The convective numerical flux

For the convective flux $\widehat{\mathbf{u}}_h^{\mathbf{w}}$ in 3.3, we take the standard upwind flux introduced in 1518. For an element $K\in {\mathcal{T}}_h$, we set

$$\begin{alignat}{1} \widehat{\mathbf{u}}^{\mathbf{w}}_h(\mathbf{x})= \begin{cases} \lim _{\epsilon \downarrow 0} \mathbf{u}_h\left(\mathbf{x}-\epsilon \mathbf{w}(\mathbf{x})\right), &\qquad \mathbf{x}\in \partial K\setminus \Gamma _-,\\ \mathbf{0},&\qquad \mathbf{x}\in \partial K\cap \Gamma _-, \end{cases} {\tag{3.5}} \end{alignat}$$

where $\Gamma _-$ is the inflow part of $\Gamma$ given by

$$\begin{equation*} \Gamma _-=\{\,\mathbf{x}\in \Gamma :\ \mathbf{w}(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})< 0\,\}. \end{equation*}$$

The diffusive numerical fluxes

If a face $e$ lies inside the domain $\Omega$, we take

$$\begin{equation} \widehat{\underline{\sigma }}_h=\{\!\!\{\underline{\sigma }_h\}\!\!\}- \kappa [\![\mathbf{u}_h\otimes \mathbf{n}]\!],\qquad \widehat{\mathbf{u}}^\sigma _h=\{\!\!\{\mathbf{u}_h\}\!\!\}, \tag{3.6} \end{equation}$$

and, if $e$ lies on the boundary, we take

$$\begin{equation} \widehat{\underline{\sigma }}_h=\underline{\sigma }_h -\kappa \mathbf{u}_h\otimes \mathbf{n},\qquad \widehat{\mathbf{u}}^\sigma _h=\mathbf{0}. \tag{3.7} \end{equation}$$

As will be shown later, the role of the parameter $\kappa$ is to ensure the stability of the method; see also 6.

The numerical fluxes related to the incompressibility constraint

The numerical fluxes associated with the incompressibility constraint, $\mathbf{\widehat{u}}^p_h$ and $\widehat{p}_h$, are defined by using an analogous recipe. If the face $e$ lies on the interior of $\Omega$, we take

$$\begin{equation} \mathbf{\widehat{u}}^p_h=\{\!\!\{\mathbf{u}_h\}\!\!\},\qquad \widehat{p}_h=\{\!\!\{p_h\}\!\!\}. {\tag{3.8}} \end{equation}$$

On the boundary, we set

$$\begin{equation} \mathbf{\widehat{u}}^p_h=\mathbf{0},\qquad \widehat{p}_h=p_h. {\tag{3.9}} \end{equation}$$

This completes the definition of the LDG method for the Oseen problem in 3.1.

Remark 3.1.

Notice that one can take the div-conforming velocity space

$$\begin{equation} \widetilde{\mathbf{V}}_h = \{ \mathbf{v}\in \mathbf{V}_h : \nabla \cdot \mathbf{v}\in L^2(\Omega )\}, {\tag{3.10}} \end{equation}$$

while keeping the other spaces and definitions above unchanged.

3.3. The post-processing operator

To complete the definition of the LDG method for the Navier-Stokes equations 1.1, it only remains to introduce what we refer to as the post-processing operator $\mathbb{P}$ and to set $\mathbf{w}={\mathbb{P}}\mathbf{u}_h$ in the approximation 3.2–3.4.

For a piecewise smooth velocity field $\mathbf{u}$, we define the operator ${\mathbb{P}}$ by

$${\mathbb{P}}\mathbf{u}|_K ={\mathbb{P}}_K \big (\mathbf{u}|_K, \widehat{\mathbf{u}}^p\big ), \qquad K\in {\mathcal{T}}_h,$$

where $\widehat{\mathbf{u}}^p$ is the numerical flux 3.8–3.9 related to the incompressibility constraint. For each element, the local operator ${\mathbb{P}}_K$ is given via the moments

$$\begin{alignat*}{2} \int _e \, {\mathbb{P}}_K\mathbf{u}\cdot \mathbf{n}_K\, \varphi \, ds & = \int _e \, \widehat{\mathbf{u}}^p\cdot \mathbf{n}_K\, \varphi \, ds & \qquad \qquad & \forall \varphi \in {P}_{k}(e), \ \text{for any edge $e\subset \partial K$},\\ \int _K\, {\mathbb{P}}_K\mathbf{u}\cdot \nabla \varphi \,d\mathbf{x} & = \int _K\, \mathbf{u}\cdot \nabla \varphi \, d\mathbf{x} & \qquad \qquad & \forall \varphi \in {P}_{k-1}(K),\\ \int _K \, {\mathbb{P}}_K\mathbf{u}\cdot \mathbf{\Psi }\, d\mathbf{x} & = \int _K \, \mathbf{u}\cdot \mathbf{\Psi } \, d\mathbf{x} & \qquad \qquad & \forall \mathbf{\Psi } \in \mathbf{\Psi }_{k}(K), \end{alignat*}$$

where

$$\mathbf{\Psi }_k(K) =\{\mathbf{\Psi }\in L^2(K)^2: DF_K^t \mathbf{\Psi }\circ F_K \in \mathbf{\Psi }_k(\widehat{K})\}.$$

Here, $F_K:\widehat{K}\to K$ denotes the elemental mapping and $DF_K$ its Jacobian. On the reference triangle $\widehat{K}=\{(\widehat{x}_1,\widehat{x}_2):\widehat{x}_1>0, \ \widehat{x}_1+\widehat{x}_2<1 \}$, the space $\mathbf{\Psi }_k(\widehat{K})$ is defined by

$$\begin{equation*} \mathbf{\Psi }_k(\widehat{K})= \{ \mathbf{\Psi } \in P_k(\widehat{K})^2: \ \text{$\nabla \cdot \mathbf{\Psi }=0$ in $\widehat{K}$}, \ \text{$\mathbf{\Psi }\cdot \mathbf{n}_{\widehat{K}}=0$ on $\partial \widehat{K}$} \}. \end{equation*}$$

The post-processing operator $\mathbb{P}$ is well defined and can be computed in an element-by-element fashion. Moreover, if $\mathbf{u}_h\in \mathbf{V}_h$ satisfies the equations 3.4, that is, if it is weakly incompressible, then $\mathbf{w}=\mathbb{P} \mathbf{u}_h$ is exactly incompressible. These results are gathered in the next result and are given in terms of the Piola transformation, which maps any vector field $\widehat{\mathbf{v}}$ on the reference triangle $\widehat{K}$ into

$$\begin{equation*} \mathsf{P}_{K} \widehat{\mathbf{v}} = \det (DF_K)^{-1}\, DF_K\,\widehat{\mathbf{v}}\circ F^{-1}_K,\qquad K\in {\mathcal{T}}_h, \end{equation*}$$

and of the BDM projection on $\widehat{K}$, $\mathbb{P}^{\mathrm{BDM}}_{\widehat{K}}$; see 4.

Proposition 3.2.

We have the following results.

(1)

$\mathbb{P}\mathbf{u}$ is well defined and $\mathbb{P}\mathbf{u}$ is in the space $\widetilde{\mathbf{V}}_h = \{ \mathbf{v}\in \mathbf{V}_h : \nabla \cdot \mathbf{v}\in L^2(\Omega )\}.$

(2)

If $\mathbf{u}\in H^1_0(\Omega )^2$ and $K\in {\mathcal{T}}_h$, then ${\mathbb{P}}_K \mathbf{u} = \mathsf{P}_{K}\, \mathbb{P}^{\mathrm{BDM}}_{\widehat{K}}\, \mathsf{P}^{-1}_{K} \mathbf{u}$.

(3)

If $\mathbf{u}\in \mathbf{V}_h$ satisfies 3.4, then $\nabla \cdot \mathbb{P}\mathbf{u}=0$ in $\Omega$ and $\mathbb{P}\mathbf{u}\in \mathbf{J}({\mathcal{T}}_h)$.

Proof.

The proof of the first assertion is straightforward, and that of the second easily follows from the definitions of the projections and from the fact that if $\mathbf{u}\in H^1_0(\Omega )^2$, then $\widehat{\mathbf{u}}^p=\mathbf{u}$.

To prove the third assertion, we first observe that $\nabla \cdot {\mathbb{P}}\mathbf{u}\in Q_h$. This is due to the fact that $\nabla \cdot {\mathbb{P}}\mathbf{u}|_K\in {P}_{k-1}(K)$ for all $K\in {\mathcal{T}}_h$ and

$$\int _{\Omega }\, \nabla \cdot {\mathbb{P}}\mathbf{u}\,d\mathbf{x} =\int _{\Gamma }\, {\mathbb{P}}\mathbf{u}\cdot \mathbf{n}\, ds =\int _{\Gamma }\, \widehat{\mathbf{u}}^p\cdot \mathbf{n}\, ds =0,$$

in view of the definitions of ${\mathbb{P}}$ and $\widehat{\mathbf{u}}^p$ in 3.9.

Now, let $\mathbf{u}\in \mathbf{V}_h$ satisfy 3.4. For $q\in Q_h$, we obtain

$$\begin{eqnarray*} \int _\Omega \, \nabla \cdot {\mathbb{P}}\mathbf{u}\,q \, d\mathbf{x} &=&\sum _{K\in {\mathcal{T}}_h}\left( - \int _K\, {\mathbb{P}}\mathbf{u}\cdot \nabla q\, d\mathbf{x} +\int _{\partial K}\, {\mathbb{P}}\mathbf{u}\cdot \mathbf{n}_K q \, ds \right)\\ &=&\sum _{K\in {\mathcal{T}}_h}\left( -\int _K\, \mathbf{u}\cdot \nabla q \, d\mathbf{x} +\int _{\partial K}\, \widehat{\mathbf{u}}^p \cdot \mathbf{n}_K q \, ds \right)\\ &=&0. \end{eqnarray*}$$

Here, we have used integration by parts, the properties of ${\mathbb{P}}$ and 3.4. Thus, we have $\nabla \cdot {\mathbb{P}}\mathbf{u}\equiv 0$ in $\Omega$. It follows that ${\mathbb{P}}\mathbf{u}\in \mathbf{J}({\mathcal{T}}_h)$.

■

Remark 3.3.

For the LDG method using the div-conforming space $\widetilde{\mathbf{V}}_h$ in 3.10, it can be readily seen that a field $\mathbf{u}\in \widetilde{\mathbf{V}}_h$ satisfying 3.4 already is exactly incompressible and belongs to $\mathbf{J}({\mathcal{T}}_h)$. Hence, for this particular LDG method, we can take ${\mathbb{P}}$ as the identity.

3.4. The mixed setting of the LDG method

Next we recast the LDG method under consideration in a classical mixed setting, not only to facilitate its analysis but to be able to state our main results in a more precise way. Thus, we eliminate the auxiliary variable $\underline{\sigma }_h$ and show that the approximation $(\mathbf{u}_h,p_h)\in \mathbf{V}_h\times Q_h$ given by the LDG method satisfies

$$\begin{alignat}{3} & A_h(\mathbf{u}_h,\mathbf{v})+ {O}_h(\mathbf{w};\mathbf{u}_h,\mathbf{v}) &&+ B_h(\mathbf{v},p_h) &&= \int _\Omega \, \mathbf{f}\cdot \mathbf{v}\,d\mathbf{x},{\tag{3.11}}\\ -&B_h(\mathbf{u}_h,q) & & & &= 0, {\tag{3.12}} \end{alignat}$$

for all $(\mathbf{v},q)\in \mathbf{V}_h\times Q_h$ where

$$\begin{equation} \mathbf{w}={\mathbb{P}}\mathbf{u}_h. {\tag{3.13}} \end{equation}$$

Here the forms $A_h$, ${O}_h$, and $B_h$ are associated to the discretization of the Laplacian, the convective term, and the incompressibility constraint, respectively. We proceed in several steps.

Step 1: Solving for $\underline{\sigma }_h$ in terms of $\mathbf{u}_h$

To be able to eliminate the auxiliary variable $\underline{\sigma }_h$, we introduce the lifting operator $\underline{{\mathcal{L}}}:\mathbf{V}_h\to \underline{\Sigma }_h$ by

$$\int _{\Omega }\, \underline{{\mathcal{L}}}(\mathbf{v}):\underline{\tau }\, d\mathbf{x} = \sum _{e\in {\mathcal{E}}_h}\int _e\,[\![\mathbf{v}\otimes \mathbf{n}]\!]: \{\!\!\{\underline{\tau }\}\!\!\}\, ds\qquad \forall \underline{\tau }\in \underline{\Sigma }_h.$$

It is now easy to see that the equation defining $\underline{\sigma }_h$ in terms of $\mathbf{u}_h$, equation 3.2, can be rewritten as

$$\begin{equation} \underline{\sigma }_h =\nu \big [ \nabla _h\mathbf{u}_h - \underline{{\mathcal{L}}}(\mathbf{u}_h)\big ], {\tag{3.14}} \end{equation}$$

with $\nabla _h$ denoting the elementwise gradient. Note that, $\underline{\sigma }_h$ can be computed from $\mathbf{u}_h$ in an element-by-element fashion. Using this identity, it is possible to eliminate $\underline{\sigma }_h$ from the equations, as we show next.

Step 2: Eliminating $\underline{\sigma }_h$

To eliminate $\underline{\sigma }_h$ from equation 3.3, we make use of a second lifting operator ${\mathcal{M}}:\mathbf{V}_h\to Q_h$ given by

$$\int _{\Omega }\,{\mathcal{M}}(\mathbf{v})\; q\,d\mathbf{x} =\sum _{e\in {\mathcal{E}}_h}\int _e \, \{\!\!\{q\}\!\!\}[\![\mathbf{v}\cdot \mathbf{n}]\!]\, ds\qquad \forall q\in Q_h.$$

If we insert the expression of $\underline{\sigma }_h$ into equation 3.3 and use the definitions of the numerical fluxes and the lifting operators, we readily get (see 2, 8, and 19)

$$\begin{equation*} A_h(\mathbf{u}_h,\mathbf{v})+ {O}_h(\mathbf{w};\mathbf{u}_h,\mathbf{v}) + B_h(\mathbf{v},p_h) = \int _\Omega \,\mathbf{f}\cdot \mathbf{v}\,d\mathbf{x}, \end{equation*}$$

where

$$\begin{alignat*}{2} & A_h(\mathbf{u},\mathbf{v}) &&:= \int _{\Omega }\, \nu \,\big [\nabla _h \mathbf{u} -\underline{{\mathcal{L}}}(\mathbf{u})\big ]: \big [\nabla _h \mathbf{v} -\underline{{\mathcal{L}}}(\mathbf{v})\big ]\, d\mathbf{x}\\ &&& \qquad \qquad \qquad \qquad \qquad +\sum _{e\in {\mathcal{E}}_h}\int _e \,\kappa \, [\![\mathbf{u}\otimes \mathbf{n}]\!]: [\![\mathbf{v}\otimes \mathbf{n}]\!]\, ds, \\ & {O}_h(\mathbf{w};\mathbf{u},\mathbf{v}) &&:= -\sum _{K\in {\mathcal{T}}_h}\int _K\, \mathbf{u}\cdot \nabla \cdot (\mathbf{v}\otimes \mathbf{w}) \, d\mathbf{x} + \sum _{K}\int _{\partial K\setminus \Gamma _-}\, \mathbf{w}\cdot \mathbf{n}_K \,\widehat{\mathbf{u}}^\mathbf{w}\cdot \mathbf{v}\, ds, \\ & B_h(\mathbf{v},q) &&:= -\int _\Omega \, q\, \nabla _h\cdot \mathbf{v} \, d\mathbf{x} +\int _{\Omega }\, q {\mathcal{M}}(\mathbf{v}) \,d\mathbf{x}. \end{alignat*}$$

This completes the elimination of the auxiliary variable $\underline{\sigma }_h$ from the equations defining the LDG method. Note that exactly the same form $B_h$ is also used in the mixed DG approaches of 11, 23, 19, 10.

Step 3: Rewriting the incompressibility constraint

Finally, it is a simple exercise to see that equation 3.4 can be rewritten as

$$\begin{equation*} -B_h(\mathbf{u}_h,q) = 0\qquad \forall q\in Q_h. \end{equation*}$$

This shows that the LDG method in 3.2–3.4 can be cast in the form given in 3.11–3.12.

4. The main results

In this section, we state and discuss the main results of this paper.

4.1. Preliminaries

We consider LDG methods with a very specific stabilization function $\kappa$. To define it, we introduce on the edges the local mesh-size function ${\mathtt{h}}$ by ${\mathtt{h}}|_e := h_e$ for all $e\in {\mathcal{E}}_h$, with $h_e$ denoting the length of the edge $e$. We then set

$$\begin{equation} \kappa :=\nu \kappa _0{\mathtt{h}}^{-1}, {\tag{4.1}} \end{equation}$$

with $\kappa _0>0$ independent of the mesh-size and the viscosity $\nu$.

The results will be stated in terms of norms we introduce next. We consider the space

$$\begin{equation} \mathbf{V}(h):=H^1_0(\Omega )^2+\mathbf{V}_h, {\tag{4.2}} \end{equation}$$

endowed with the norm

$$\begin{equation*} \|\mathbf{v}\|^2_{1,h}:=\sum _{K\in {\mathcal{T}}_h} \|\nabla \mathbf{v}\|_{0,K}^2+ \sum _{e\in {\mathcal{E}}_h}\int _e \, \kappa _0 {\mathtt{h}}^{-1} \,|[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2\, ds. \end{equation*}$$

For this norm, we have the Poincaré inequality

$$\begin{equation} \|\mathbf{v}\|_0 \leq C_p \|\mathbf{v}\|_{1,h}\qquad \forall \mathbf{v}\in \mathbf{V}(h), {\tag{4.3}} \end{equation}$$

for a constant $C_p>0$ independent of the mesh-size; see, e.g., 13. Finally, the space $Q_h$ is equipped with the $L^2$-norm $\|\cdot \|_{0}$.

4.2. Stability properties of the bilinear forms

Here we collect crucial stability properties of the forms that are used to define the LDG method. Our main results will be stated in terms of the corresponding stability constants.

Continuity

First we study the continuity of the forms involved in the LDG method. We observe that the lifting operators $\underline{{\mathcal{L}}}$ and ${\mathcal{M}}$ can be extended to operators $\underline{{\mathcal{L}}}:\mathbf{V}(h)\to \Sigma _h$ and ${\mathcal{M}}:\mathbf{V}(h)\to Q_h$, respectively, using the same definitions. It is then well known (from, e.g., 16, Section 3 and 19, Lemma 7.2) that

$$\begin{alignat}{2} \|\underline{{\mathcal{L}}}(\mathbf{v})\|^2_0 & \leq C_{\operatorname {lift}}^2 \sum _{e\in {\mathcal{E}}_h} \int _e\,\kappa _0{\mathtt{h}}^{-1}|[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2\, ds, & \qquad & \mathbf{v}\in \mathbf{V}(h),{\tag{4.4}}\\ \|{\mathcal{M}}(\mathbf{v})\|^2_0 & \leq C_{\operatorname {lift}}^2 \sum _{e\in {\mathcal{E}}_h}\int _e\, \kappa _0{\mathtt{h}}^{-1}|[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2\, ds, & \qquad & \mathbf{v}\in \mathbf{V}(h), {\tag{4.5}} \end{alignat}$$

for a constant $C_{\operatorname {lift}}$ only depending on $\kappa _0$, the shape-regularity of the mesh, and the polynomial degree $k$. As a consequence, the forms $A_h$ and $B_h$ are well defined and continuous on $\mathbf{V}(h)\times \mathbf{V}(h)$ and $\mathbf{V}(h)\times L^2(\Omega )$.

Proposition 4.1.

We have that

$$\begin{alignat*}{2} |A_h(\mathbf{u},\mathbf{v})| & \leq \nu \,C_a \|\mathbf{u}\|_{1,h} \|\mathbf{v}\|_{1,h} & \qquad &\forall \mathbf{u},\mathbf{v}\in \mathbf{V}(h),\\ |B_h(\mathbf{v},q)| & \leq C_b \|\mathbf{v}\|_{1,h} \|q\|_0 & \qquad &\forall (\mathbf{v},q)\in \mathbf{V}(h)\times L^2(\Omega ), \end{alignat*}$$

for continuity constants $C_a$ and $C_b$ that are independent of the mesh-size.

Proof.

The continuity properties of $A_h$ and $B_h$ follow from 4.4, 4.5, and the Cauchy-Schwarz inequality; see 16, Proposition 3.1 and 19, Lemma 7.5.

■

Proposition 4.2.

Let $\mathbf{w}_1,\mathbf{w}_2\in \mathbf{J}({\mathcal{T}}_h)$, $\mathbf{u}\in \mathbf{V}(h)$, and $\mathbf{v}\in \mathbf{V}_h$. Then we have

$$|O_h(\mathbf{w}_1;\mathbf{u},\mathbf{v})-O_h(\mathbf{w}_2;\mathbf{u},\mathbf{v})|\, \leq \, C_o \, \|\mathbf{w}_1-\mathbf{w}_2\|_{1,h} \,\|\mathbf{u}\|_{1,h}\, \|\mathbf{v}\|_{1,h},$$

for a continuity constant $C_o$ that is independent of the mesh-size.

The proof is given in subsection 5.1. It is obtained by using the embedding and trace theorems of 13 and 10.

Coercivity and inf-sup condition

Next, we discuss the coercivity properties of the forms $A_h$ and $O_h$. We have the following result.

Proposition 4.3.

Let $\kappa$ be given by 4.1. Then, for any $\kappa _0 > 0$, there exists a constant $\alpha >0$ independent of the mesh-size such that

$$\begin{equation*} A_h(\mathbf{v},\mathbf{v})\geq \nu \alpha \|\mathbf{v}\|^2_{1,h} \qquad \forall \, \mathbf{v}\in \mathbf{V}_h. \end{equation*}$$

Furthermore, for $\mathbf{w}\in \mathbf{J}({\mathcal{T}}_h)$ and $\mathbf{v}\in \mathbf{V}_h$, it holds that

$${O}_h(\mathbf{w};\mathbf{v},\mathbf{v}) = \frac{1}{2}\sum _{e\in {\mathcal{E}}_h^{\mathcal{I}}}\int _e\, |\mathbf{w}\cdot \mathbf{n}|\,|[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2\,ds +\frac{1}{2}\int _{\Gamma }\, |\mathbf{w}\cdot \mathbf{n}| |\mathbf{v}|^2\,ds.$$

In the integrals over edges in $\mathcal{E}_h^{\mathcal{I}}$, we denote by $\mathbf{n}$ any unit normal to the edge $e$ under consideration.

For a proof of the coercivity of the form $A_h$, we refer to 2 or 16. A similar coercivity result involving also the discrete velocity gradient was used in 9 and 8. We further note that, for the similar symmetric interior penalty forms $A_h$ used in the DG approach of 11, the parameter $\kappa _0$ has to be chosen large enough. The proof of the second assertion is standard; see, e.g., 8.

Finally, we have the following inf-sup condition for the form $B_h$.

Proposition 4.4 (11).

There exists a constant $\beta >0$ independent of the mesh-size such that

$$\sup _{\mathbf{0}\neq \mathbf{v}\in \mathbf{V}_h}\frac{B_h(\mathbf{v},q)}{\|\mathbf{v}\|_{1,h}}\geq \beta \,\|q\|_{0}\qquad \forall q\in Q_h.$$

Extensions of this result to the $hp$-version of the finite element method and to quadrilateral meshes can be found in 23, 19.

Remark 4.5.

A careful inspection of the proof in 11 reveals that the discrete inf-sup condition in Proposition 4.4 also holds for the smaller space $\widetilde{\mathbf{V}}_h$ in 3.10. Consequently, all the stability results of this section hold for the particular LDG method in Remark 3.1.

Stability of the post-processing operator

The next result states that the operator ${\mathbb{P}}$ is a bounded linear operator from $\mathbf{V}_h$ to $\widetilde{\mathbf{V}}_h$ with respect to the norm $\|\cdot \|_{1,h}$. It is one of the main technical results needed to analyze the locally conservative LDG method 3.11–3.13.

Proposition 4.6.

Let $\mathbf{v}\in \mathbf{V}_h$. Then we have

$$\|\,{\mathbb{P}}\mathbf{v}\,\|_{1,h}\leq C_{\mathrm{stab}} \|\,\mathbf{v}\,\|_{1,h},$$

with a stability constant $C_{\mathrm{stab}}>0$ that is independent of the mesh-size.

The proof of this proposition is given in subsection 5.2 below. We remark that, for the LDG method in Remark 3.1, we have $C_{\mathrm{stab}}=1$ since ${\mathbb{P}}$ is chosen as the identity.

We are now ready to state our main results.

4.3. The results

Our first result states that, under a smallness condition similar to the one for the exact solution 1.7, the LDG method 3.11–3.13 defines a unique discrete approximation. Moreover, it actually gives an efficient way to compute it.

Theorem 4.7 (Existence and uniqueness of discrete solutions).

Assume that

$$\begin{equation} \mu :=\frac{C_o\,C_{\mathrm{stab}}\,C_p \|\mathbf{f}\|_0}{\nu ^2\,\alpha ^2} < 1. {\tag{4.6}} \end{equation}$$

Then the LDG method 3.11–3.13 defines a unique solution $(\mathbf{u}_h,p_h)\in \mathbf{V}_h\times Q_h$. It satisfies the bounds

$$\begin{alignat}{2} &\|\mathbf{u}_h\|_{1,h}&&\leq \frac{C_p\|\mathbf{f}\|_0}{\nu \,\alpha },{\tag{4.7}}\\ &\|p_h\|_0 &&\leq \beta ^{-1} \left(\frac{C_a+2\alpha }{\alpha }\right) \,C_p\,\|\mathbf{f}\|_0, {\tag{4.8}} \end{alignat}$$

as well as

$$\begin{equation} \sum _{e\in {\mathcal{E}}_h^{\mathcal{I}}}\int _e\, |{\mathbb{P}}\mathbf{u}_h\cdot \mathbf{n}|\,|[\![\mathbf{u}_h\otimes \mathbf{n}]\!]|^2\,ds +\int _{\Gamma }\, |{\mathbb{P}}\mathbf{u}_h\cdot \mathbf{n}| |\mathbf{u}_h|^2\,ds \leq \frac{C_p^2 \|\mathbf{f}\|_0^2}{\nu \,\alpha }. {\tag{4.9}} \end{equation}$$

Moreover, if $(\mathbf{u}^{\ell +1}_{h},p_h^{\ell +1})$ is the approximate solution given by the LDG method in 3.11–3.12 for the Oseen equations with $\mathbf{w}=\mathbb{P} \mathbf{u}^\ell _h$, $\ell \ge 0$, then

$$\begin{alignat*}{2} &\|\mathbf{u}_h-\mathbf{u}^{\ell +1}_{h}\|_{1,h}&&\le \;2\, \left( \frac{C_p\|\mathbf{f}\|_0}{\nu \alpha }\right)\, \frac{\mu ^\ell }{(1-\mu )},\\ &\|p_h-p_h^{\ell +1}\|_0 &&\le \;2\,\beta ^{-1}\;\left(\frac{C_a+2\alpha }{\alpha }\right) \,C_p\|\mathbf{f}\|_0\, \frac{\mu ^{\ell }}{(1-\mu )}, \end{alignat*}$$

for any initial guess $(\mathbf{u}^{0}_{h},p_h^{0})\in \mathbf{V}_h\times Q_h$.

This result, whose proof is given in subsection 5.3, states that the LDG method in 3.11–3.13 is well defined and that we can compute its approximate solution by solving a sequence of Oseen problems. Since the parameter $\mu$ is independent of the mesh-size, the convergence of that sequence is always exponential and so computationally efficient.

Note that if we set

$$\begin{equation} \alpha _\star =\min \{1,\alpha \},\qquad C_{\mathrm{poinc}} = \max \{C_P,C_p\},\qquad C_O = \max \{ C_\Omega , C_o C_{\mathrm{stab}}\}, {\tag{4.10}} \end{equation}$$

then both the smallness assumptions in 1.7 and 4.6 are satisfied if we have that

$$\frac{C_O C_{\mathrm{poinc}} \|\mathbf{f}\|_0}{\nu ^2\, \alpha _\star ^2} < 1.$$

Hence, both the Navier-Stokes equations and their LDG approximations are uniquely solvable. Under a smallness condition that is slightly more restrictive, we obtain the following estimates.

Theorem 4.8 (Error estimates).

Assume that

$$\begin{equation} \frac{C_O C_{\mathrm{poinc}} \|\mathbf{f}\|_0}{\nu ^2\,\alpha _\star ^2} \leq \frac{1}{2}, {\tag{4.11}} \end{equation}$$

and that the exact solution $(\mathbf{u},p)$ of the Navier-Stokes equations 1.1 satisfies

$$\begin{equation} \mathbf{u}\in H^{s+1}(\Omega )^2,\qquad p\in H^s(\Omega ),\qquad s\geq 1. {\tag{4.12}} \end{equation}$$

Then

$$\begin{alignat*}{3} &\|\mathbf{u}-\mathbf{u}_h\|_{1,h} &\le &\,C_{\mathbf{u}}\,{C}_{\mathrm{app}}\, \big [\|\mathbf{u}\|_{s+1}+\nu ^{-1}\,\|p\|_s \big ]\,&&h^{\min \{k,s\}}, \\ &\|\mathbf{u}-\mathbb{P}\mathbf{u}_h\|_{1,h} &\le &\,C_{\mathbf{w}}\,{C}_{\mathrm{app}}\, \big [\|\mathbf{u}\|_{s+1}+\nu ^{-1}\,\|p\|_s \big ]\,&&h^{\min \{k,s\}}, \\ & \|p-p_h\|_0 &\le &\,C_{p}\, {C}_{\mathrm{app}}\, \big [ \nu \,\|\mathbf{u}\|_{s+1}+\|p\|_s \big ]\,&&h^{\min \{k,s\}}, \\ & \|\underline{\sigma }-\underline{\sigma }_h\|_0 &\le &\,C_{\sigma }\,{C}_{\mathrm{app}}\, \big [\nu \,\|\mathbf{u}\|_{s+1}+\|p\|_s \big ]\,&&h^{\min \{k,s\}}, \end{alignat*}$$

where ${C}_{\mathrm{app}}$ only depends on the regularity of the mesh and the polynomial degree $k$, and

$$\begin{alignat}{1} {C}_{\mathbf{u}} &=\max \left\{ \big (\frac{\beta +C_b}{\beta }\big ) \big (\frac{2\,C_a+3\,\alpha _\star }{\alpha _\star }\big ), \frac{1+C_{\mathrm{stab}}}{C_{\mathrm{stab}}}, \frac{2\,C_b}{\alpha _\star }, \frac{2}{\alpha _\star } \right\}, {\tag{4.13}}\\ {C}_{\mathbf{w}} &=(1+C_{\mathrm{stab}}) + C_{\mathrm{stab}} \,{C}_{\mathbf{u}}, {\tag{4.14}}\\ {C}_p &=\max \left\{ \big ( \frac{C_a+\alpha _\star }{\beta }\big )\, C_\mathbf{u}, \frac{\alpha _\star \,(1+C_{\mathrm{stab}})}{2\beta \,C_{\mathrm{stab}}}, \frac{\beta +C_b}{\beta }, \frac{1}{\beta } \right\}, {\tag{4.15}}\\ C_\sigma &=(1+C_{\mathrm{lift}}) {C}_{\mathbf{u}}. {\tag{4.16}} \end{alignat}$$

This result, whose proof is given in subsection 5.4, states that the LDG method under consideration converges with optimal order. Note also that, since the function $\mathbf{w}=\mathbb{P}\mathbf{u}_h$ is exactly divergence free, it provides an optimally convergent globally solenoidal approximation to the velocity!

Let us briefly discuss some extensions of these results.

• First, we point out that all the results of this section are valid verbatim for the LDG method in Remark 3.1 where $\mathbf{V}_h$ is replaced by the div-conforming space $\widetilde{\mathbf{V}}_h$ in 3.10 and ${\mathbb{P}}$ is chosen to be the identity.

• Although here we only considered the case of triangular meshes, the results of this paper can be straightforwardly extended to simplicial meshes in three dimensions. Furthermore, the LDG approach we propose here can be easily extended to $Q_k-Q_k-P_{k-1}$ elements on quadrilateral or hexahedral affine meshes by using a post-processing operator ${\mathbb{P}}$ that is given by a slight modification of the BDM projection on quadrilaterals or hexahedra; see 5. The results in this section then hold true for this LDG method as well. This fact is actually verified in our numerical experiments for which we have used square meshes and $Q_1-Q_1-P_0$ elements. However, the extension of our results to $Q_k-Q_k-Q_{k-1}$ elements on quadrilateral meshes is not straightforward. Although, by using a post-processing operator ${\mathbb{P}}$ that is a slight modification of the standard Raviart-Thomas projection, it is easy to define a solenoidal velocity field $\mathbf{w}$ that belongs to the anisotropic polynomial space $Q_{k-1,k}\times Q_{k,k-1}$. The approximation properties in this space give rise to only suboptimal convergence rates. If, on the other hand, the polynomial degree of the post-processed velocity is increased, the field $\mathbf{w}$ cannot be shown to be solenoidal, as $\nabla \cdot \widetilde{\mathbf{V}}_h\not \subset Q_h$.

• Finally, let us remark that here we have used the LDG approach to discretize the viscous terms. However, our results remain valid for any other DG discretization of these terms whose primal bilinear form $A_h$ is both coercive and continuous, such as, e.g., the interior penalty form. For detail, we refer the reader to the discussions in 2 and 19.

5. Proofs

In this section, we provide the proofs of our main results.

5.1. Proof of Proposition 4.2

We begin by proving Proposition 4.2. To do this, note that, if we insert the definition of the upwinding numerical flux into the form $O_h$, we have

$$\begin{eqnarray*} O_h(\mathbf{w};\mathbf{u},\mathbf{v})&=&-\sum _{K\in {\mathcal{T}}_h}\int _K\,\mathbf{u}\cdot \nabla \cdot (\mathbf{v}\otimes \mathbf{w})\,d\mathbf{x}\\ &&+\sum _{K\in {\mathcal{T}}_h}\, \int _{\partial K}\, \bigg [\mathbf{w}\cdot \mathbf{n}_K\{\!\!\{\mathbf{u}\}\!\!\}-\frac{1}{2}|\mathbf{w}\cdot \mathbf{n}_K|(\mathbf{u}^{\mathrm{ext}}-\mathbf{u})\bigg ] \cdot \mathbf{v}\, ds. \end{eqnarray*}$$

Here, $\mathbf{u}^{\mathrm{ext}}$ denotes the exterior trace of $\mathbf{u}$ taken over the edge under consideration and set to zero on the boundary. If we perform now a simple integration by parts, we get

$$\begin{eqnarray*} O_h(\mathbf{w};\mathbf{u},\mathbf{v})&=&\sum _{K\in {\mathcal{T}}_h}\int _K\,\nabla \mathbf{u}:(\mathbf{v}\otimes \mathbf{w})\,d\mathbf{x}\\ &&+\sum _{K\in {\mathcal{T}}_h}\, \int _{\partial K}\, \bigg [\frac{1}{2}\mathbf{w}\cdot \mathbf{n}_K(\mathbf{u}^{\mathrm{ext}}-\mathbf{u})-\frac{1}{2}|\mathbf{w}\cdot \mathbf{n}_K|(\mathbf{u}^{\mathrm{ext}}-\mathbf{u})\bigg ] \cdot \mathbf{v}\, ds. \end{eqnarray*}$$

This implies that

$$\begin{equation*} O_h(\mathbf{w}_1;\mathbf{u},\mathbf{v})-O_h(\mathbf{w}_2;\mathbf{u},\mathbf{v})=T_1 +T_2, \end{equation*}$$

where

$$\begin{alignat*}{2} &T_1 =&& \sum _{K\in {\mathcal{T}}_h} \int _K\,\nabla \mathbf{u}:\big (\mathbf{v}\otimes (\mathbf{w}_1-\mathbf{w}_2)\big )\,d\mathbf{x}, \\ & T_2=&& \sum _{K\in {\mathcal{T}}_h}\int _{\partial K}\, \frac{1}{2}(\mathbf{w}_1-\mathbf{w}_2)\cdot \mathbf{n}_K (\mathbf{u}^{\mathrm{ext}}-\mathbf{u})\cdot \mathbf{v}\,ds\\ &&& - \sum _{K\in {\mathcal{T}}_h}\int _{\partial K}\, \frac{1}{2} \bigg ( |\mathbf{w}_1\cdot \mathbf{n}_K|-|\mathbf{w}_2\cdot \mathbf{n}_K| \bigg ) (\mathbf{u}^{\mathrm{ext}}-\mathbf{u})\cdot \mathbf{v}\, ds. \end{alignat*}$$

To bound the term $T_1$, we recall the following embedding result proved in 13, Proposition 4.5 for smooth and convex domains and in 10, Lemma 6.2 for general polygons:

$$\|\mathbf{v}\|_{L^4(\Omega )}\leq C \|\mathbf{v}\|_{1,h},\qquad \mathbf{v}\in \{\mathbf{w}\in L^2(\Omega )^2: \mathbf{w}|_K\in H^1(K)^2,\ K\in \mathcal{T}_h\},$$

with a constant independent of the mesh-size. (We point out that the broken $H^1$-norm used in 13 is slightly different than the one we use here. However, a careful inspection of the proof of Proposition 4.5 therein shows that, in fact the result holds for our definition of $\|\cdot \|_{1,h}$.) It is then clear that we can use Hölder’s inequality to obtain

$$|T_1| \leq \|\mathbf{w}_1-\mathbf{w}_2\|_{L^4(\Omega )} \|\mathbf{u}\|_{1,h} \|\mathbf{v}\|_{L^4(\Omega )} \leq C \|\mathbf{w}_1-\mathbf{w}_2\|_{1,h} \|\mathbf{u}\|_{1,h} \|\mathbf{v}\|_{1,h}.$$

It remains to estimate the term $T_2$. Using the Lipschitz continuity of the function $x\mapsto |x|$, we get

$$|T_2|\leq \sum _{K\in {\mathcal{T}}_h}\int _{\partial K}\, |\mathbf{w}_1\cdot \mathbf{n}_K-\mathbf{w}_2\cdot \mathbf{n}_K| \, |[\![\mathbf{u}\otimes \mathbf{n}]\!]|\, |\mathbf{v}|\, ds,$$

and, proceeding as in the proof of 13, Proposition 4.5, we obtain

$$|T_2|\leq C \|\mathbf{w}_1-\mathbf{w}_2\|_{1,h} \|\mathbf{u}\|_{1,h} \|\mathbf{v}\|_{1,h},$$

with a constant independent of the mesh-size. This completes the proof of Proposition 4.2.

5.2. Proof of Proposition 4.6

We prove Proposition 4.6 by first establishing local stability bounds over patches of elements and then by summing up these local results.

Let $\mathbf{v}\in \mathbf{V}_h$ be fixed. We proceed in several steps.

Step 1: Local bounds in the interior. Let $e=\partial K_1 \cap \partial K_2$ be an interior edge shared by two elements $K_1$ and $K_2$. We wish to establish a local stability bound over the patch formed by $K_1$ and $K_2$. Namely, by defining the local seminorm

$$|\mathbf{v}|_e^2 = \|\nabla \mathbf{v}\|_{0,K_1}^2+ \|\nabla \mathbf{v}\|_{0,K_2}^2 +\int _e \, {\mathtt{h}}^{-1} |[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2\, ds,$$

we claim that

$$\begin{equation} |{\mathbb{P}}\mathbf{v}|_e^2 \, \leq \, C \, \big [\, |\mathbf{v}|_e^2 + \int _{\partial K_1}\, {\mathtt{h}}^{-1} |(\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_{K_1}|^2\, ds +\int _{\partial K_2}\, {\mathtt{h}}^{-1}|(\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_{K_2}|^2\, ds \big ], {\tag{5.1}} \end{equation}$$

with a constant $C$ independent of the mesh-size.

To prove 5.1, it is enough to consider the case where $K_1$ and $K_2$ form a reference patch of unit size. The general case then follows from a scaling argument and the shape-regularity assumption on the mesh, by mapping $K_1\cup K_2$ onto the reference patch using elementwise Piola transforms.

By the triangle inequality, we have

$$\begin{equation} |{\mathbb{P}}\mathbf{v}|_e \leq |\mathbf{v}-{\mathbb{P}}\mathbf{v}|_e + |\mathbf{v}|_e. {\tag{5.2}} \end{equation}$$

It remains to bound $|\mathbf{v}-{\mathbb{P}}\mathbf{v}|_e$. To do so, we note that, for any element $K\in {\mathcal{T}}_h$, the restriction $(\mathbf{v}-{\mathbb{P}}\mathbf{v})|_K$ belongs to $P_k(K)^2$ and is uniquely defined by the moments

$$\begin{alignat*}{2} \int _e \, (\mathbf{v}-{\mathbb{P}}\mathbf{v})\cdot \mathbf{n}_K\, \varphi \, ds & = \int _e \, (\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_K\, \varphi \, ds & \qquad & \forall \varphi \in {P}_{k}(e), \ e\subset \partial K,\\ \int _K \, (\mathbf{v}-{\mathbb{P}}\mathbf{v}) \cdot \nabla \varphi \, d\mathbf{x} & = 0 & \qquad & \forall \varphi \in P_{k-1}(K),\\ \int _K \, (\mathbf{v}-{\mathbb{P}}\mathbf{v}) \cdot \mathbf{\Psi }\, d\mathbf{x} & = 0 & \qquad & \forall \mathbf{\Psi } \in \mathbf{\Psi }_{k}(K). \end{alignat*}$$

Hence, the restriction of $\mathbf{v}-{\mathbb{P}}\mathbf{v}$ to the patch $K_1\cup K_2$ belongs to the space

$$\begin{equation} \begin{split} \mathbf{\mathsf{V}}_e = \bigg \{ \mathbf{\mathsf{v}}\in L^2(K_1\cup K_2)^2:\ & \mathbf{\mathsf{v}}|_{K_i}\in P_k(K_i)^2,\\ & \int _{K_i}\,\mathbf{\mathsf{v}}|_{K_i}\cdot \nabla \varphi \,d\mathbf{x}=0, \ \varphi \in P_{k-1}(K_i),\\ & \int _{K_i}\,\mathbf{\mathsf{v}}|_{K_i}\cdot \mathbf{\Psi }\,d\mathbf{x}=0, \ \mathbf{\Psi }\in \mathbf{\Psi }_k(K_i),\ i=1,2\, \bigg \}. \end{split} {\tag{5.3}} \end{equation}$$

Furthermore, it can be easily seen that the mappings $\mathbf{\mathsf{v}}\mapsto \|\mathbf{\mathsf{v}}\|_e$ and $\mathbf{\mathsf{v}}\mapsto |\!|\!| {\mathbf{\mathsf{v}}} |\!|\!|_e$, given by

$$\begin{eqnarray*} \|\mathbf{\mathsf{v}}\|_e^2&=& |\mathbf{\mathsf{v}}|_e^2+ \int _{\partial K_1\setminus e}\, {\mathtt{h}}^{-1} |\mathbf{\mathsf{v}}\cdot \mathbf{n}_{K_1}|^2\, ds +\int _{\partial K_2\setminus e}\, {\mathtt{h}}^{-1}|\mathbf{\mathsf{v}}\cdot \mathbf{n}_{K_2}|^2\, ds,\\ |\!|\!| {\mathbf{\mathsf{v}}} |\!|\!|_e^2 &=& \int _{\partial K_1} \, {\mathtt{h}}^{-1}|\mathbf{\mathsf{v}} \cdot \mathbf{n}_{K_1}|^2\, ds + \int _{\partial K_2} \, {\mathtt{h}}^{-1}|\mathbf{\mathsf{v}}\cdot \mathbf{n}_{K_2}|^2\, ds, \end{eqnarray*}$$

define norms on $\mathbf{\mathsf{V}}_e$. By the equivalence of all norms on a finite dimensional space, there holds

$$\|\mathbf{\mathsf{v}}\|_e \leq C |\!|\!| {\mathbf{\mathsf{v}}} |\!|\!|_e \qquad \forall \mathbf{\mathsf{v}}\in \mathbf{\mathsf{V}}_e,$$

with a constant only depending on the the polynomial degree. Thus, we obtain that

$$|\mathbf{v}-{\mathbb{P}}\mathbf{v}|_e \leq \|\mathbf{v}-{\mathbb{P}}\mathbf{v}\|_e \leq C |\!|\!| {\mathbf{v}-{\mathbb{P}}\mathbf{v}} |\!|\!|_e.$$

On the other hand, since, for $i=1,2$,

$$(\mathbf{v}-{\mathbb{P}}\mathbf{v})\cdot \mathbf{n}_{K_i} = (\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_{K_i} \qquad \text{on $\partial K_i$},$$

we conclude that

$$\begin{eqnarray*} |\mathbf{v}-{\mathbb{P}}\mathbf{v}|_e^2 &\leq & C |\!|\!| {\mathbf{v}-{\mathbb{P}}\mathbf{v}} |\!|\!|_e^2\\ &=& C \int _{\partial K_1}\, {\mathtt{h}}^{-1} |(\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_{K_1}|^2\, ds + C \int _{\partial K_2}\, {\mathtt{h}}^{-1}|(\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_{K_2}|^2\, ds. \end{eqnarray*}$$

This estimate and the inequality in 5.2 prove the local stability bound in 5.1.

Step 2: Local bounds on the boundary. The analogous stability result holds on the boundary. Let $K$ be an element on the boundary and $e\subset \partial K$ a boundary edge. By setting

$$|\mathbf{v}|_e^2 = \|\nabla \mathbf{v}\|_{0,K}^2 +\int _e \, {\mathtt{h}}^{-1} |\mathbf{v}\otimes \mathbf{n}|^2\, ds,$$

there exists a constant $C$ independent of the mesh-size such that

$$\begin{equation} |{\mathbb{P}}\mathbf{v}|_e^2 \leq C \big [ |\mathbf{v}|_e^2 + \int _{\partial K}\, {\mathtt{h}}^{-1} |(\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_{K}|^2\, ds \big ]. {\tag{5.4}} \end{equation}$$

Step 3: Summing up the local bounds. We complete the proof of Proposition 4.6 by summing up the local stability estimates established in 5.1 and 5.4. To this end, we first note that $\mathbf{v}-\widehat{\mathbf{v}}^p=\frac{1}{2}(\mathbf{v}-\mathbf{v}^{\mathrm{ext}})$ on interior edges and $\mathbf{v}-\widehat{\mathbf{v}}^p=\mathbf{v}$ on boundary edges. Here, we write $\mathbf{v}^{\mathrm{ext}}$ to denote the exterior trace of $\mathbf{v}$ over the edge under consideration. Therefore, we have for any edge $e\in {\mathcal{E}}_h$

$$\int _e \, {\mathtt{h}}^{-1} |(\widehat{\mathbf{v}}^p-\mathbf{v})\cdot \mathbf{n}_K|^2\, ds \leq \int _e \, {\mathtt{h}}^{-1} |[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2 \, ds.$$

Using the local bounds in 5.1 and 5.4 and the above estimate for $\mathbf{v}-\widehat{\mathbf{v}}^p$, we obtain

$$\begin{eqnarray*} \|{\mathbb{P}}\mathbf{v}\|_{1,h}^2 & \leq & C \sum _{e\in {\mathcal{E}}_h} |{\mathbb{P}}\mathbf{v}|_{e}^2\\ &\leq & C \sum _{e\in {\mathcal{E}}_h} |\mathbf{v}|_{e}^2+ C \sum _{K\in {\mathcal{T}}_h} \int _{\partial K}\, {\mathtt{h}}^{-1}|(\mathbf{v}-\widehat{\mathbf{v}}^p)\cdot \mathbf{n}_K|^2\, ds\\ &\leq & C \sum _{e\in {\mathcal{E}}_h} |\mathbf{v}|_{e}^2+ C \sum _{e \in {\mathcal{E}}_h}\int _e \, {\mathtt{h}}^{-1} |[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2 \, ds\\ &\leq & C \|\mathbf{v}\|_{1,h}^2, \end{eqnarray*}$$

with constants $C$ independent of the mesh-size. This completes the proof of Proposition 4.6.

5.3. Proof of Theorem 4.7

To prove Theorem 4.7, we proceed as follows. First, we eliminate the pressure from the problem by restricting ourselves to the weakly divergence-free subspace of $\mathbf{V}_h$,

$$\begin{equation} \mathbf{Z}_h = \left\{\mathbf{v}\in \mathbf{V}_h : B_h(\mathbf{v},q) = 0\ \quad \forall q\in Q_h\right\}. {\tag{5.5}} \end{equation}$$

The approximate velocity is thus characterized as the only function $\mathbf{u}_h\in \mathbf{Z}_h$ such that

$$\begin{equation} A_h(\mathbf{u}_h,\mathbf{v})+ {O}_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v}) = \int _\Omega \, \mathbf{f}\cdot \mathbf{v}\,d\mathbf{x} \qquad \forall \mathbf{v}\in \mathbf{Z}_h. {\tag{5.6}} \end{equation}$$

Then, we construct a contractive mapping $S$ defined on a ball of $\mathbf{Z}_h$ whose only fixed point is precisely the above velocity. The properties for the corresponding pressure $p_h$ then follow from its characterization

$$\begin{equation*} B_h(\mathbf{v},p_h)=\int _\Omega \, \mathbf{f}\cdot \mathbf{v}\,d\mathbf{x} -A_h(\mathbf{u}_h,\mathbf{v}) - {O}_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v}) \qquad \forall \mathbf{v}\in \mathbf{V}_h/\mathbf{Z}_h, \end{equation*}$$

and from the $\inf$-$\sup$ condition for the incompressibility form $B_h$.

We proceed in several steps.

Step 1: The operator $S$. We begin by introducing the operator $S$. For $\overline{\mathbf{u}}\in \mathbf{Z}_h$, $\mathbf{u}=S(\overline{\mathbf{u}})$ denotes the solution of the following problem. Find $\mathbf{u}\in \mathbf{Z}_h$ such that

$$A_h(\mathbf{u},\mathbf{v})+ {O}_h({\mathbb{P}}\overline{\mathbf{u}};\mathbf{u},\mathbf{v}) = \int _\Omega \, \mathbf{f}\cdot \mathbf{v}\,d\mathbf{x} \qquad \forall \mathbf{v}\in \mathbf{Z}_h.$$

Note that since $\overline{\mathbf{u}}\in \mathbf{Z}_h$ we have, by Proposition 3.2, that $\mathbb{P}\overline{\mathbf{u}}\in \mathbf{J}_h({\mathcal{T}}_h)$. As a consequence, this problem is uniquely solvable.

Furthermore, by the coercivity of the form $A_h$ and $O_h$ in Proposition 4.3,

$$\begin{alignat*}{2} & \nu \alpha \|\mathbf{u}\|_{1,h}^2 && \le A_h(\mathbf{u},\mathbf{u})+{O}_h({\mathbb{P}}\overline{\mathbf{u}};\mathbf{u},\mathbf{u})\\ &&& = \int _\Omega \, \mathbf{f}\cdot \mathbf{u}\,d\mathbf{x}\\ &&& \le \|\mathbf{f}\|_0 \|\mathbf{u}\|_0. \end{alignat*}$$

By the Poincaré inequality in 4.3, we obtain

$$\nu \alpha \|\mathbf{u}\|_{1,h}^2 \leq C_p \|\mathbf{f}\|_0 \|\mathbf{u}\|_{1,h}.$$

Hence, the solution $\mathbf{u}$ to the above problem satisfies

$$\begin{equation} \|\mathbf{u}\|_{1,h} \leq \frac{C_p \|\mathbf{f}\|_0}{\nu \alpha }. {\tag{5.7}} \end{equation}$$

This implies that $S$ maps ${\mathcal{K}}_h$ into ${\mathcal{K}}_h$, where

$${\mathcal{K}}_h = \{ \mathbf{v}\in \mathbf{Z}_h: \ \|\mathbf{v}\|_{1,h}\leq \frac{C_p \|\mathbf{f}\|_0}{\nu \alpha } \}.$$

Step 2: The operator $S$ is a contraction. Next, we show that $S$ is a contraction on ${\mathcal{K}}_h$ under the smallness condition 4.6. To do so, let $\overline{\mathbf{u}}_1,\overline{\mathbf{u}}_2$ be in ${\mathcal{K}}_h$, and set $\mathbf{u}_1=S(\overline{\mathbf{u}}_1)$, $\mathbf{u}_2=S(\overline{\mathbf{u}}_2)$. Then

$$\begin{alignat*}{2} &\nu \alpha \|\mathbf{u}_1-\mathbf{u}_2\|_{1,h}^2 \le && A_h(\mathbf{u}_1-\mathbf{u}_2,\mathbf{u}_1-\mathbf{u}_2). \end{alignat*}$$

Since

$$A_h(\mathbf{u}_1-\mathbf{u}_2,\mathbf{v}) + {O}_h({\mathbb{P}}\overline{\mathbf{u}}_1;\mathbf{u}_1,\mathbf{v}) -{O}_h({\mathbb{P}}\overline{\mathbf{u}}_2;\mathbf{u}_2,\mathbf{v})=0,$$

for any $\mathbf{v}\in \mathbf{Z}_h$, taking $\mathbf{v}=\mathbf{u}_1-\mathbf{u}_2$, we get

$$\begin{alignat*}{2} &\nu \alpha \|\mathbf{u}_1-\mathbf{u}_2\|_{1,h}^2 &\le &\, - {O}_h({\mathbb{P}}\overline{\mathbf{u}}_2;\mathbf{u}_1-\mathbf{u}_2,\mathbf{u}_1-\mathbf{u}_2)\\ &&& + O_h({\mathbb{P}}\overline{\mathbf{u}}_2;\mathbf{u}_1,\mathbf{u}_1-\mathbf{u}_2) - O_h({\mathbb{P}}\overline{\mathbf{u}}_1;\mathbf{u}_1,\mathbf{u}_1-\mathbf{u}_2) \\ &&=:&\; T_1 + T_2. \end{alignat*}$$

By the coercivity property of the form ${O}_h$ in Proposition 4.3,

$$\begin{equation*} T_1 \le 0. \end{equation*}$$

Moreover, by the continuity property of $O_h$ in Proposition 4.2, the bound 5.7 and the continuity of the post-processing operator $\mathbb{P}$ in Proposition 4.6,

$$\begin{alignat*}{2} & T_2&\le &\; C_o \|{\mathbb{P}}\overline{\mathbf{u}}_1-{\mathbb{P}}\overline{\mathbf{u}}_2\|_{1,h} \|\mathbf{u}_1\|_{1,h}\|\mathbf{u}_1-\mathbf{u}_2\|_{1,h} \\ &&\leq &\; \frac{C_o C_p \|\mathbf{f}\|_0}{\nu \alpha } \|{\mathbb{P}}\overline{\mathbf{u}}_1-{\mathbb{P}}\overline{\mathbf{u}}_2\|_{1,h} \|\mathbf{u}_1-\mathbf{u}_2\|_{1,h} \\ &&\leq &\;\frac{C_o C_{\mathrm{stab}} C_p \|\mathbf{f}\|_0}{\nu \alpha } \|\overline{\mathbf{u}}_1-\overline{\mathbf{u}}_2\|_{1,h} \|\mathbf{u}_1-\mathbf{u}_2\|_{1,h} \\ &&= &\; \nu \alpha \mu \, \|\overline{\mathbf{u}}_1-\overline{\mathbf{u}}_2\|_{1,h} \|\mathbf{u}_1-\mathbf{u}_2\|_{1,h}. \end{alignat*}$$

This implies that

$$\begin{equation*} \|\mathbf{u}_1-\mathbf{u}_2\|_{1,h} \leq \;\mu \; \|\overline{\mathbf{u}}_1-\overline{\mathbf{u}}_2\|_{1,h}, \end{equation*}$$

and so, if $\mu <1$, that is, if the smallness condition 4.6 is satisfied, the mapping $S$ is a contraction. Hence, $S$ has a unique fixed point $\mathbf{u}_h\in {\mathcal{K}}_h$, which is the solution to the problem 5.6.

Step 3: Recovering the pressure. Now that the velocity $\mathbf{u}_h$ has been computed, the pressure is the solution $p_h\in Q_h$ of

$$\begin{equation} B_h(\mathbf{v},p_h)=\int _\Omega \, \mathbf{f}\cdot \mathbf{v}\,d\mathbf{x} -A_h(\mathbf{u}_h,\mathbf{v}) - {O}_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v})\qquad \forall \mathbf{v}\in \mathbf{V}_h/\mathbf{Z}_h. {\tag{5.8}} \end{equation}$$

Due to Proposition 4.1, Proposition 4.2, and the Poincaré inequality in 4.3, the right-hand side defines a continuous linear functional on $\mathbf{V}_h/\mathbf{Z}_h$. The inf-sup condition in Proposition 4.4 then guarantees the existence of a unique solution $p_h$ to the above problem. It can then be easily seen that the tuple $(\mathbf{u}_h,p_h)$ is the unique solution to the LDG method in 3.11 and 3.12 with $\mathbf{w}={\mathbb{P}}\mathbf{u}_h$.

Step 4: The stability bounds. Next, let us show the stability bounds for $(\mathbf{u}_h,p_h)$ in Theorem 4.7. The bound for $\|\mathbf{u}_h\|_{1,h}$ in 4.7 follows in a straightforward way since $\mathbf{u}_h\in {\mathcal{K}}_h$. To obtain the bound for the upwind term in 4.9, note that

$$\nu \alpha \|\mathbf{u}_h\|_{1,h}^2 + O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{u}_h) \leq C_p \|\mathbf{f}\|_0 \|\mathbf{u}_h\|_{1,h}^2 \leq \frac{1}{2} \frac{C_p^2 \|\mathbf{f}\|_0^2}{\nu \alpha } + \frac{1}{2} \nu \alpha \|\mathbf{u}_h\|_{1,h}^2.$$

Similarly to the previous arguments, here we have used the coercivity of $A_h$, equation 5.6 with test function $\mathbf{v}=\mathbf{u}_h$, and the Poincaré inequality 4.3. Bringing the term $\frac{1}{2} \nu \alpha \|\mathbf{u}_h\|_{1,h}^2$ to the left-hand side and observing the coercivity of $O_h$ give the stability bound in 4.9.

Moreover, using the inf-sup condition in Proposition 4.4, the Poincaré inequality in 4.3, the continuity properties in Proposition 4.1 and Proposition 4.2, and the stability of ${\mathbb{P}}$ in Proposition 4.6, we have from 5.8

$$\begin{equation*} \beta \|p_h\|_0 \leq \sup _{\mathbf{0}\neq \mathbf{v}\in \mathbf{V}_h} \frac{B_h(\mathbf{v},p_h)}{\|\mathbf{v}\|_{1,h}} \leq C_p\|\mathbf{f}\|_0 + \nu \,C_a \|\mathbf{u}_h\|_{1,h}+C_o C_{\mathrm{stab}} \|\mathbf{u}_h\|_{1,h}^2. \end{equation*}$$

Taking into account the stability bound for $\mathbf{u}_h$ and assumption 4.6 gives

$$\begin{eqnarray*} \|p_h\|_0 & \leq & \beta ^{-1} \bigg ( C_p\|\mathbf{f}\|_0 + \frac{C_a C_p \|\mathbf{f}\|_0}{\alpha } +\frac{C_o C_{\mathrm{stab}} C_p \|\mathbf{f}\|_0}{\nu ^2\alpha ^2} C_p \|\mathbf{f}\|_0 \bigg )\\ & \leq & \beta ^{-1} C_p\|\mathbf{f}\|_0 \left( 2+\frac{C_a}{\alpha }\right). \end{eqnarray*}$$

This gives the desired bound 4.8 for $p_h$.

Step 5: The convergence estimates. It remains to prove the error estimates for the sequence $\{(\mathbf{u}^\ell _h,p_h^\ell )\}_{\ell \ge 0}$. Let us begin with that of the velocity. From Step 3, we have that $\mathbf{u}_h^{\ell +1}= S(\mathbf{u}_h^\ell )$. As a consequence, since $S$ is a contraction with Lipschitz constant $\mu$, we immediately get

$$\begin{equation*} \|\mathbf{u}_h-\mathbf{u}^{\ell +1}_{h}\|_{1,h} \le \left( \frac{\mu ^\ell }{(1-\mu )}\right)\, \|\mathbf{u}^2_h-\mathbf{u}^1_{h}\|_{1,h}. \end{equation*}$$

The result now follows from the fact that, by the stability bound 5.7,

$$\begin{equation*} \|\mathbf{u}^m_h\|_{1,h} \leq \frac{C_p \|\mathbf{f}\|_0}{\nu \alpha }, \end{equation*}$$

for $m \ge 1$.

To obtain the estimate for the pressure, we proceed as follows. First, we note that, from Step 4, we have

$$\begin{equation*} B_h(\mathbf{v},p^{\ell +1}_h) =\int _\Omega \, \mathbf{f}\cdot \mathbf{v}\,d\mathbf{x} -A_h(\mathbf{u}_h^{\ell +1},\mathbf{v}) - {O}_h({\mathbb{P}}\mathbf{u}^\ell _h;\mathbf{u}^{\ell +1}_h,\mathbf{v}) \qquad \forall \mathbf{v}\in \mathbf{V}_h. \end{equation*}$$

This implies that

$$\begin{eqnarray*} B_h(\mathbf{v},p_h-p^{\ell +1}_h) &=& -A_h(\mathbf{u}_h-\mathbf{u}_h^{\ell +1},\mathbf{v}) - {O}_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v}) + O_h({\mathbb{P}}\mathbf{u}_h^\ell ;\mathbf{u}_h,\mathbf{v})\\ && - O_h({\mathbb{P}}\mathbf{u}_h^\ell ;\mathbf{u}_h,\mathbf{v}) + {O}_h({\mathbb{P}}\mathbf{u}^\ell _h;\mathbf{u}^{\ell +1}_h,\mathbf{v})\\ &=& -A_h(\mathbf{u}_h-\mathbf{u}_h^{\ell +1},\mathbf{v}) - {O}_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v}) + O_h({\mathbb{P}}\mathbf{u}_h^\ell ;\mathbf{u}_h,\mathbf{v})\\ && \phantom {+} {O}_h({\mathbb{P}}\mathbf{u}^\ell _h;\mathbf{u}^{\ell +1}_h-\mathbf{u}_h,\mathbf{v}). \end{eqnarray*}$$

We insert this expression in the inf-sup condition for $B_h$, use the stability properties of $A_h$, ${O}_h$, and $\mathbb{P}$, take into account the bounds for $\|\mathbf{u}_h\|_{1,h}$, $\|\mathbf{u}_h^\ell \|_{1,h}$, and the contraction property of $S$ to obtain

$$\begin{eqnarray*} \|p_h-p^{\ell +1}_h\|_0 &\leq & \beta ^{-1} \left(\nu C_a \mu + \frac{C_o C_{\mathrm{stab}} C_p \|\mathbf{f}\|_0}{\nu \alpha } (1+\mu )\right) \|\mathbf{u}_h-\mathbf{u}_h^{\ell }\|_{1,h}\\ &\leq & \beta ^{-1}\mu \left(\nu C_a + 2 \nu \alpha \right) \|\mathbf{u}_h-\mathbf{u}_h^{\ell }\|_{1,h}. \end{eqnarray*}$$

The desired bound for $\|p_h-p_h^{\ell +1}\|_0$ then follows from the bound for $\|\mathbf{u}_h-\mathbf{u}_h^{\ell }\|_{1,h}$.

This completes the proof of Theorem 4.7.

5.4. Proof of Theorem 4.8

Here, we derive the error estimates in Theorem 4.8. To do that, we modify the approach used in 8 to get error estimates for the LDG method for the Oseen problem in two ways. First, we use the nonconforming approach introduced in 16 and later used in 19, and consider the expression

$$\begin{equation*} {\mathcal{R}}_h(\mathbf{u},p):=\sup _{\mathbf{0}\neq \mathbf{v}\in \mathbf{V}_h} \frac{|{R}_h(\mathbf{u},p;\mathbf{v})|}{\|\mathbf{v}\|_{1,h}}, \end{equation*}$$

where

$$\begin{equation*} {R}_h(\mathbf{u},p;\mathbf{v}):= A_{h}(\mathbf{u},\mathbf{v})+O_h(\mathbf{u};\mathbf{u},\mathbf{v})+B_h(\mathbf{v},p)- \int _{\Omega }\,\mathbf{f}\cdot \mathbf{v}\,d\mathbf{x}, \qquad \mathbf{v}\in \mathbf{V}_h. \end{equation*}$$

The second modification is, of course, due to the presence of the convective nonlinearity.

Proof.

We proceed in several steps.

Step 1: The abstract estimate for $\|\mathbf{u}-\mathbf{u}\|_{1,h}$. Let us begin with the estimate for the error $\|\mathbf{u}-\mathbf{u}\|_{1,h}$. We claim that we have

$$\begin{alignat}{1} \|\mathbf{u}-\mathbf{u}_h\|_{1,h} \leq {C}_{\mathbf{u}} \; \bigg [ &\inf _{\mathbf{v}\in \mathbf{V}_h} \|\mathbf{u}-\mathbf{v}\|_{1,h}+ \inf _{\widetilde{\mathbf{v}}\in \widetilde{\mathbf{V}}_h}\|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} {\tag{5.9}}\\ & +\inf _{q\in Q_h} \frac{1}{\nu }\|p-q\|_0+\frac{1}{\nu }{\mathcal{R}}_h(\mathbf{u},p)\bigg ], \end{alignat}$$

where $C_{\mathbf{u}}$ is given by 4.13.

To prove this result, we proceed as in the error analysis of standard mixed methods (see, e.g., 5) and consider first an element $\mathbf{v}\in \mathbf{Z}_h$ where $\mathbf{Z}_h$ is the kernel in 5.5. We then have the trivial inequality

$$\begin{equation} \|\mathbf{u} -\mathbf{u}_h\|_{1,h}\le \|\mathbf{u} -\mathbf{v}\|_{1,h} +\|\mathbf{v} -\mathbf{u}_h\|_{1,h}. {\tag{5.10}} \end{equation}$$

Next, we obtain an estimate of $\mathbf{v} -\mathbf{u}_h$. By the coercivity property of the form $A_h$ in Proposition 4.3, we have

$$\begin{equation} \nu \, \alpha _{\star }\;\|\mathbf{v} -\mathbf{u}_h\|^2_{1,h} \le \nu \, \alpha \; \|\mathbf{v} -\mathbf{u}_h\|^2_{1,h} \le A_h(\mathbf{v} -\mathbf{u}_h,\mathbf{v} -\mathbf{u}_h), {\tag{5.11}} \end{equation}$$

and since

$$\begin{equation*} A_h(\mathbf{u}-\mathbf{u}_h,\mathbf{\mathsf{v}})+O_h(\mathbf{u};\mathbf{u},\mathbf{\mathsf{v}}) -O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{\mathsf{v}})+B_h(\mathbf{\mathsf{v}},p-p_h) =R_h(\mathbf{u},p;\mathbf{\mathsf{v}}) \end{equation*}$$

for any $\mathbf{\mathsf{v}}\in \mathbf{V}_h$, for $\mathbf{\mathsf{v}}=\mathbf{v} -\mathbf{u}_h$, we get that

$$\begin{alignat*}{2} &A_h(\mathbf{v} -\mathbf{u}_h,\mathbf{v} -\mathbf{u}_h) &= &\; A_h(\mathbf{v} -\mathbf{u},\mathbf{v} -\mathbf{u}_h) \\ &&&\;-O_h(\mathbf{u};\mathbf{u},\mathbf{v} -\mathbf{u}_h) +O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v} -\mathbf{u}_h) \\ &&&\; -B_h(\mathbf{v} -\mathbf{u}_h,p-p_h) \\ &&&\;+R_h(\mathbf{u},p;\mathbf{v} -\mathbf{u}_h) \\ &&=:&\; T_1+T_2+T_3+T_4. \end{alignat*}$$

Hence,

$$\begin{equation*} \alpha _{\star }\,\nu \;\|\mathbf{v} -\mathbf{u}_h\|^2_{1,h}\le T_1+T_2+T_3+T_4. \end{equation*}$$

The terms $T_1$, $T_3$, and $T_4$ can be easily estimated as follows. First, by the continuity of the form $A_h$,

$$\begin{equation*} T_1 \leq \nu \,C_a \|\mathbf{u}-\mathbf{v}\|_{1,h}\|\mathbf{v}-\mathbf{u}_h\|_{1,h}. \end{equation*}$$

Then, by the definition of $R_h$ and ${\mathcal{R}}_h$,

$$\begin{equation*} T_4 \leq {\mathcal{R}}_h(\mathbf{u},p)\;\|\mathbf{v}-\mathbf{u}_h\|_{1,h}. \end{equation*}$$

Finally, since $\mathbf{v}-\mathbf{u}_h\in \mathbf{Z}_h$, we have

$$\begin{equation*} B_h(\mathbf{v}-\mathbf{u}_h,p-p_h)=B_h(\mathbf{v}-\mathbf{u}_h,p)=B_h(\mathbf{v}-\mathbf{u}_h,p-q), \end{equation*}$$

for any $q\in Q_h$. From Proposition 4.1, it follows that

$$\begin{equation*} T_3\leq C_b\|\mathbf{v}-\mathbf{u}_h\|_{1,h} \|p-q\|_0\qquad \forall \, q\in Q_h. \end{equation*}$$

It remains to estimate the term $T_2$. To do that, consider the identity

$$\begin{alignat*}{2} &T_2&=&\;-O_h({\mathbb{P}}\mathbf{v};\mathbf{u},\mathbf{v}-\mathbf{u}_h) +O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u},\mathbf{v}-\mathbf{u}_h) \\ &&&\; -O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{v}-\mathbf{u}_h,\mathbf{v}-\mathbf{u}_h) \\ &&&\; +O_h({\mathbb{P}}\mathbf{v};\mathbf{u},\mathbf{v}-\mathbf{u}_h) -O_h(\mathbf{u};\mathbf{u},\mathbf{v}-\mathbf{u}_h) \\ &&& +O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{v}-\mathbf{u},\mathbf{v}-\mathbf{u}_h) \\ \ &&=:&\; T_{21}+ T_{22}+ T_{23}+ T_{24}. \end{alignat*}$$

Note that due to Propositions 4.2 and 4.6, the stability bound for $\mathbf{u}$ in 1.5, and the definitions of the parameters in 4.10, we have

$$\begin{eqnarray*} T_{21} &\leq & C_o\, C_{\mathrm{stab}}\|\mathbf{u}\|_1 \|\mathbf{v}-\mathbf{u}_h\|_{1,h}^2\\ &\leq & \frac{C_o\, C_{\mathrm{stab}}\, C_P\, \|\mathbf{f}\|_0}{\nu }\|\mathbf{v}-\mathbf{u}_h\|_{1,h}^2\\ &\leq & \frac{C_O\, C_{\mathrm{poinc}}\|\mathbf{f}\|_0}{\nu \alpha _\star } \|\mathbf{v}-\mathbf{u}_h\|_{1,h}^2\\ &\leq &\frac{1}{2}\nu \alpha _\star \; \|\mathbf{v}-\mathbf{u}_h\|_{1,h}^2, \end{eqnarray*}$$

by the smallness condition 4.11.

Next, by the coercivity property of ${O}_h$ in Proposition 4.3,

$$\begin{equation*} T_{22}\le 0, \end{equation*}$$

and, by Propositions 4.2 and 4.6, and the bound for $\mathbf{u}_h$ in Theorem 4.7,

$$\begin{alignat*}{2} &T_{24} & \leq &\; \frac{C_o C_{\mathrm{stab}} C_p \|\mathbf{f}\|_0}{\nu \alpha } \|\mathbf{u}-\mathbf{v}\|_{1,h}\|\mathbf{v}-\mathbf{u}_h\|_{1,h} \\ &&\le &\;\frac{C_O C_{\mathrm{poinc}}\|\mathbf{f}\|_0}{\nu \alpha _\star } \|\mathbf{u}-\mathbf{v}\|_{1,h}\|\mathbf{v}-\mathbf{u}_h\|_{1,h} \\ &&\le &\frac{1}{2}\nu \alpha _\star \; \|\mathbf{u}-\mathbf{v}\|_{1,h}\|\mathbf{v}-\mathbf{u}_h\|_{1,h}, \end{alignat*}$$

by the smallness condition 4.11.

Finally, by the Lipschitz property of the form $O_h$ in Proposition 4.2, we have

$$\begin{alignat*}{2} &T_{23} &\leq & C_o \|\mathbf{u}-{\mathbb{P}}\mathbf{v}\|_{1,h} \|\mathbf{u}\|_1 \|\mathbf{v}-\mathbf{u}_h\|_{1,h}\\ && \leq & \frac{C_o C_P \|\mathbf{f}\|_0}{\nu } \|\mathbf{u}-{\mathbb{P}}\mathbf{v}\|_{1,h}\|\mathbf{v}-\mathbf{u}_h\|_{1,h} \\ && \leq & \frac{\nu \alpha _\star }{2\,C_{\mathrm{stab}}}\, \|\mathbf{u}-{\mathbb{P}}\mathbf{v}\|_{1,h}\|\mathbf{v}-\mathbf{u}_h\|_{1,h}, \end{alignat*}$$

by the bound for $\mathbf{u}$ in 1.5 and the smallness condition 4.11. Now, take an arbitrary function $\widetilde{\mathbf{v}}$ in $\widetilde{\mathbf{V}}_h$. Since ${\mathbb{P}}$ reproduces functions in $\widetilde{\mathbf{V}}_h$, we have ${\mathbb{P}}\widetilde{\mathbf{v}}=\widetilde{\mathbf{v}}$, and so

$$\begin{eqnarray*} \|\mathbf{u}-{\mathbb{P}}\mathbf{v}\|_{1,h} &\leq & \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} +\|{\mathbb{P}}\widetilde{\mathbf{v}}-{\mathbb{P}}\mathbf{v}\|_{1,h}\\ &\leq & \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} + C_{\mathrm{stab}} \|\widetilde{\mathbf{v}}-\mathbf{v}\|_{1,h}\\ &\leq & (1+C_{\mathrm{stab}}) \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h}+C_{\mathrm{stab}}\|\mathbf{u}-\mathbf{v}\|_{1,h}, \end{eqnarray*}$$

by Proposition 4.6. This implies that

$$\begin{alignat*}{2} &T_{23} &\leq & \frac{\nu \alpha _\star }{2} \left(\left(\frac{1+C_{\mathrm{stab}}}{C_{\mathrm{stab}}}\right)\, \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h}+\|\mathbf{u}-\mathbf{v}\|_{1,h} \right)\,\|\mathbf{v}-\mathbf{u}_h\|_{1,h}. \end{alignat*}$$

Thus, gathering all the estimates above and inserting them in the right-hand side of inequality 5.11, and bringing $T_{21}$ to the left-hand side, we obtain

$$\begin{alignat*}{1} \|\mathbf{v}-\mathbf{u}_h\|_{1,h} &\leq 2\,\left(\frac{C_a+\alpha _\star }{\alpha _\star }\right) \|\mathbf{u}-\mathbf{v}\|_{1,h}+ \left(\frac{1+C_{\mathrm{stab}}}{C_{\mathrm{stab}}}\right) \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} \\ &\quad + \left(\frac{2\,C_b}{\nu \alpha _\star }\right) \|p-q\|_0 +\left(\frac{2}{\nu \alpha _\star }\right) {\mathcal{R}}_h(\mathbf{u},p). \end{alignat*}$$

Inserting this estimate in the right-hand side of inequality 5.10, we get

$$\begin{equation} \begin{split} \|\mathbf{u}-\mathbf{u}_h\|_{1,h} \leq & \left(\frac{2 C_a+3 \alpha _\star }{\alpha _\star }\right) \|\mathbf{u}-\mathbf{v}\|_{1,h}+ \left(\frac{1+C_{\mathrm{stab}}}{C_{\mathrm{stab}}}\right) \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} \\ & + \left(\frac{2\,C_b}{\nu \alpha _\star }\right) \|p-q\|_0 +\left(\frac{2}{\nu \alpha _\star }\right) {\mathcal{R}}_h(\mathbf{u},p), \end{split} {\tag{5.12}} \end{equation}$$

for any $\mathbf{v}\in \mathbf{Z}_h$, $\widetilde{\mathbf{v}}\in \widetilde{\mathbf{V}}_h$, and $q\in Q_h$.

It remains to replace $\mathbf{v}\in \mathbf{Z}_h$ by an arbitrary function in $\mathbf{v}\in \mathbf{V}_h$. To this end, fix $\mathbf{v}\in \mathbf{V}_h$ and consider the problem: Find $\mathbf{r}\in \mathbf{V}_h$ such that

$$B_h(\mathbf{r},\mathsf{q})=B_h(\mathbf{u}-\mathbf{v},\mathsf{q}) \qquad \forall \mathsf{q}\in Q_h.$$

The inf-sup condition in Proposition 4.4 guarantees that such a solution $\mathbf{r}$ exists. Furthermore, it can be easily seen that we have

$$\|\mathbf{r}\|_{1,h}\leq \beta ^{-1} C_b \|\mathbf{u}-\mathbf{v}\|_{1,h},$$

in view of the inf-sup condition for $B_h$ and the continuity of the form $B_h$. By construction and since $B_h(\mathbf{u},\mathsf{q})=0$ for any $\mathsf{q}\in Q_h$, we further have that $\mathbf{r}+\mathbf{v}\in \mathbf{Z}_h$. Inserting $\mathbf{r}+\mathbf{v}$ in 5.12, employing the triangle inequality, and taking into account the above bound for $\mathbf{r}$ yield the abstract error estimate 5.9 for the velocity.

Step 2: The abstract estimate for $\|\mathbf{u}-\mathbb{P}\mathbf{u}_h\|_{1,h}$. As a consequence of the approximation result 5.9, we obtain the estimate of the error between $\mathbf{u}$ and its globally solenoidal approximation $\mathbb{P}\mathbf{u}_h$

$$\begin{equation} \|\mathbf{u}-\mathbb{P}\mathbf{u}_h\|_{1,h} \le (1+C_{\mathrm{stab}}) \inf _{\widetilde{\mathbf{v}}\in \widetilde{\mathbf{V}}_h} \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h}+C_{\mathrm{stab}} \|\mathbf{u}-\mathbf{u}_h\|_{1,h}. {\tag{5.13}} \end{equation}$$

To see this, note that

$$\begin{eqnarray*} \|\mathbf{u}-\mathbb{P}\mathbf{u}_h\|_{1,h} &\leq & \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h}+ \|\widetilde{\mathbf{v}}-\mathbb{P}\mathbf{u}_h\|_{1,h}\\ &\leq & \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h}+ C_{\mathrm{stab}} \|\widetilde{\mathbf{v}}-\mathbf{u}_h\|_{1,h}\\ &\leq & (1+C_{\mathrm{stab}}) \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} +C_{\mathrm{stab}}\|\mathbf{u}-\mathbf{u}_h\|_{1,h}, \end{eqnarray*}$$

where $\widetilde{\mathbf{v}}$ is any element of $\widetilde{\mathbf{V}}_h$. Here we have used the stability bound in Proposition 4.6 and the fact that $\mathbb{P}$ reproduces polynomials in $\widetilde{\mathbf{V}}_h$. This shows the inequality 5.13.

Step 3: The abstract estimate for the pressure. Now, let us obtain the estimate for the pressure. We claim that the error in the pressure satisfies

$$\begin{alignat}{1} \|p-p_h\|_{0} \leq {C}_p\; \bigg [ &\inf _{\mathbf{v}\in \mathbf{V}_h} \nu \;\|\mathbf{u}-\mathbf{v}\|_{1,h}+ \inf _{\widetilde{\mathbf{v}} \in \widetilde{\mathbf{V}}_h} \nu \;\|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} {\tag{5.14}}\\ &+\inf _{q\in Q_h} \|p-q\|_0+{\mathcal{R}}_h(\mathbf{u},p)\bigg ], \end{alignat}$$

where $C_p$ is given by 4.15.

To see this, we proceed in a way similar to the one used to deal with the velocity. Thus, we begin by noting that for $q\in Q_h$

$$\begin{alignat*}{2} &\|p-p_h\|_{0} &\; \leq &\; \|q-p_h\|_0 + \|p-q\|_0 \\ && \leq &\; \beta ^{-1} \sup _{\mathbf{v}\in \mathbf{V}_h} \frac{B_h(\mathbf{v},q-p_h)}{\|\mathbf{v}\|_{1,h}} +\|p-q\|_0 \\ &&\leq &\; \beta ^{-1} \sup _{\mathbf{v}\in \mathbf{V}_h} \frac{B_h(\mathbf{v},p-p_h)}{\|\mathbf{v}\|_{1,h}} + \beta ^{-1}\sup _{\mathbf{v}\in \mathbf{V}_h} \frac{B_h(\mathbf{v},q-p)}{\|\mathbf{v}\|_{1,h}}+\|p-q\|_0, \end{alignat*}$$

where we have used the inf-sup condition in Proposition 4.4. Therefore,

$$\begin{equation} \|p-p_h\|_0\leq \beta ^{-1} \sup _{\mathbf{v}\in \mathbf{V}_h} \frac{B_h(\mathbf{v},p-p_h)}{\|\mathbf{v}\|_{1,h}} + (1+\beta ^{-1} C_b)\, \|p-q\|_0, {\tag{5.15}} \end{equation}$$

by the continuity of the form $B_h$.

To bound the first term on the right-hand side of 5.15, we note that

$$\begin{equation*} B_h(\mathbf{v},p-p_h) = -A_h(\mathbf{u}-\mathbf{u}_h,\mathbf{v})-O_h(\mathbf{u};\mathbf{u},\mathbf{v}) +O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h,\mathbf{v}) +R_h(\mathbf{u},p;\mathbf{v}), \end{equation*}$$

for any $\mathbf{v}\in \mathbf{V}_h$, and proceed as in the previous step to obtain

$$\begin{alignat*}{2} & B_h(\mathbf{v},p-p_h) \;&= &\; A_h(\mathbf{u}-\mathbf{u}_h,\mathbf{v}) + O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u}_h-\mathbf{u},\mathbf{v}) \\ & &&+O_h({\mathbb{P}}\mathbf{u}_h;\mathbf{u},\mathbf{v})-O_h(\mathbf{u};\mathbf{u},\mathbf{v}) +R_h(\mathbf{u},p;\mathbf{v})\\ &&\leq &\bigg [ \, (\nu \,C_a+ \frac{1}{2}\nu \,\alpha _\star )\,\|\mathbf{u}-\mathbf{u}_h\|_{1,h}\\ &&& + \frac{\nu \,\alpha _\star }{2\,C_{\mathrm{stab}}}\, \|\mathbf{u}-{\mathbb{P}}\mathbf{u}_h\|_{1,h} +{\mathcal{R}}_h(\mathbf{u},p)\bigg ]\; \|\mathbf{v}\|_{1,h}, \end{alignat*}$$

and since

$$\begin{equation*} \|\mathbf{u}-{\mathbb{P}}\mathbf{u}_h\|_{1,h} \leq (1+C_{\mathrm{stab}}) \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} +C_{\mathrm{stab}}\|\mathbf{u}-\mathbf{u}_h\|_{1,h}, \end{equation*}$$

from Step 2, we get

$$\begin{alignat*}{1} B_h(\mathbf{v},p-p_h) \leq & (C_a +\alpha _\star )\,\nu \, \|\mathbf{u}-\mathbf{u}_h\|_{1,h}\\ & \qquad \qquad +\frac{\alpha _\star }{2} \, \left(\frac{1+C_{\mathrm{stab}}}{C_{\mathrm{stab}}}\right) \,\nu \, \|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} +{\mathcal{R}}_h(\mathbf{u},p). \end{alignat*}$$

Inserting this inequality in 5.15 and using the bound 5.9 for $\|\mathbf{u}-\mathbf{u}_h\|_{1,h}$ from Step 1, we immediately obtain the abstract estimate 5.14 for the pressure.

Step 4: The abstract estimate for the velocity gradient. To derive the error estimate for $\underline{\sigma }-\underline{\sigma }_h$, we note that from 3.14 we have

$$\underline{\sigma }-\underline{\sigma }_h=\nu \big [ \nabla \mathbf{u}-\nabla _h\mathbf{u}_h + \underline{{\mathcal{L}}}(\mathbf{u}_h)\big ].$$

Hence,

$$\|\underline{\sigma }-\underline{\sigma }_h\|_0 \leq \nu \|\mathbf{u}-\mathbf{u}_h\|_{1,h}+\nu \|\underline{{\mathcal{L}}}(\mathbf{u}_h)\|_0.$$

Using the stability bound 4.4 for the lifting operator $\underline{\mathcal{L}}$ yields

$$\|\underline{{\mathcal{L}}}(\mathbf{u}_h)\|^2_0\leq C_{\mathrm{lift}}^2 \sum _{e\in {\mathcal{E}}_h}\int _e\, \kappa _0{\mathtt{h}}^{-1}|[\![\mathbf{u}_h\otimes \mathbf{n}]\!]|^2\, ds \leq C_{\mathrm{lift}}^2 \|\mathbf{u}-\mathbf{u}_h\|^2_{1,h}.$$

The last inequality follows as the jumps of the exact solution vanish. This shows that

$$\begin{equation} \|\underline{\sigma }-\underline{\sigma }_h\|_0 \leq \nu \;(1+C_{\mathrm{lift}})\; \|\mathbf{u}-\mathbf{u}_h\|_{1,h}. {\tag{5.16}} \end{equation}$$

Step 5: Approximation estimates. Under the regularity assumption 4.12, the standard approximation property

$$\inf _{\mathbf{v}\in \mathbf{V}_h} \nu \,\|\mathbf{u}-\mathbf{v}\|_{1,h} +\inf _{q\in Q_h} \|p-q\|_{0}\leq C_{\mathrm{app}} h^{\min \{k,s\}} \big [\nu \,\|\mathbf{u}\|_{s+1}+\|p\|_s\big ]$$

holds. Moreover, from the results in 5 (see also 11, Section 3) we have

$$\inf _{\widetilde{\mathbf{v}}\in \widetilde{\mathbf{V}}_h}\|\mathbf{u}-\widetilde{\mathbf{v}}\|_{1,h} \leq C_{\mathrm{app}} h^{\min \{k,s\}} \|\mathbf{u}\|_{s+1}.$$

Finally, we have that

$${\mathcal{R}}_h(\mathbf{u},p)\leq C_{\mathrm{app}} h^{\min \{k,s\}} \big [ \nu \,\|\mathbf{u}\|_{s+1} + \|p\|_s \big ],$$

with a constant $C_{\mathrm{app}}$ independent of the mesh-size.

To see the estimate of the residual, we proceed as follows. For $\mathbf{v}\in \mathbf{V}_h$, it is easy to see that ${R}_h(\mathbf{u},p;\mathbf{v})$ is given by

$${R}_h(\mathbf{u},p;\mathbf{v})=\sum _{e\in {\mathcal{E}}_h} \int _e\, \{\!\!\{\nu \nabla \mathbf{u}- \underline{P}_h(\nu \nabla \mathbf{u})\}\!\!\}:[\![\mathbf{v}\otimes \mathbf{n}]\!]\, ds -\sum _{e\in {\mathcal{E}}_h} \int _e\, \{\!\!\{p- P_h p\}\!\!\}[\![\mathbf{v}\cdot \mathbf{n}]\!]\, ds,$$

with $\underline{P}_h:L^2(\Omega )^{2\times 2}\to \underline{\Sigma }_h$ and $P_h:L^2(\Omega )/\mathbb{R}\to Q_h$ denoting the $L^2$-projections onto $\underline{\Sigma }_h$ and $Q_h$, respectively. The desired estimate then follows by proceeding as the proof of 19, Proposition 8.1 and using standard approximation results for $\underline{P}_h$ and $P_h$.

Step 6: Conclusion. It is now a simple matter to see that the error estimates of Theorem 4.8 follow by inserting the approximation estimates obtained in the previous step into the abstract bounds for the velocity 5.9, for its globally solenoidal post-processing 5.13, the pressure 5.14, and the velocity gradient 5.16. This completes the proof of Theorem 4.8.

■

6. Numerical results

In this section, we present numerical experiments that show that the theoretical rates of convergence are sharp. We also display the behavior of the iterative method as a function of the Reynolds number.

As a reference solution in our tests, we take the analytical solution $(\mathbf{u},p)$ of the incompressible Navier-Stokes equations that was obtained by Kovasznay in 14. For a given viscosity $\nu$, this solution is given by

$$\begin{align*} u_1(x,y) &= 1 - e^{\lambda x} \cos (2\pi y),\\ u_2(x,y) &= \frac{\lambda }{2\pi } e^{\lambda x} \sin (2\pi y), \\ p(x,y) &= -\frac{1}{2} e^{2\lambda x} + \bar{p}, \end{align*}$$

where

$$\begin{gather*} \lambda = \frac{-8\pi ^2}{\nu ^{-1}+\sqrt {\nu ^{-2}+64\pi ^2}}. \end{gather*}$$

It solves 1.1 with a suitably chosen right-hand side $\mathbf{f}$. Here, the constant $\bar{p}$ is such that $\int _{\Omega }p\, d\mathbf{x}=0$. Further we use the values of the exact solution $\mathbf{u}=(u_1,u_2)$ to prescribe inhomogeneous Dirichlet boundary data $\mathbf{g}$ for the velocity on the whole boundary of the computational domain, which we take to be $\Omega = (-\frac{1}{2},\frac{3}{2})\times (0,2)$. Note that in this case, the numerical fluxes (on edges $e$ lying on the boundary) must be modified as

$$\begin{alignat*}{1} & \widehat{\mathbf{u}}_h^{\mathbf{w}}(\mathbf{x}) = \begin{cases} \lim _{\epsilon \downarrow 0} \mathbf{u}_h \left(\mathbf{x}-\epsilon \mathbf{w}(\mathbf{x})\right), & \qquad \mathbf{x}\in e\setminus \Gamma _-,\\ \mathbf{g}(\mathbf{x}), & \qquad \mathbf{x}\in e\cap \Gamma _-, \end{cases} \\ & \widehat{\underline{\sigma }}_h =\underline{\sigma }_h-\kappa (\mathbf{u}_h-\mathbf{g})\otimes \mathbf{n},\qquad \widehat{\mathbf{u}}_h^\sigma =\mathbf{g}, \\ &\mathbf{\widehat{u}}^p_h=\mathbf{g},\qquad \widehat{p}_h=p_h. \end{alignat*}$$

We consider square meshes that are generated by refining the single grid cell $(-\frac{1}{2},\frac{3}{2})\times (0,2)$ uniformly. Therefore, a mesh on “level” $L$ consists of $2^L$ by $2^L$ cells. All computations are performed with bilinear shape functions for $\underline{\sigma }_h$ and $\mathbf{u}_h$ and piecewise constants for $p_h$, according to the remarks in subsection 4.3. The stabilization function $\kappa$ is chosen as in 4.1 with $\kappa _0=4$. Our analysis then predicts first-order convergence for $\mathbf{u}_h$ in the norm $\|.\|_{1,h}$ and for the pressure in the $L^2$-norm.

In Table 1 we show the errors and convergence rates in $p$, $\mathbf{u}$ and $\underline{\sigma }$ obtained for $\nu =0.1$. The errors in $p$ and $\underline{\sigma }$ are measured in the $L^2$-norm while $\mathbf{u}-\mathbf{u}_h$ and $\mathbf{u}-{\mathbb{P}}\mathbf{u}_h$ are evaluated in the norm $\|.\|_{1,h}$. We observe the predicted first-order convergence for all the error components, in full agreement with the results of Theorem 4.8. Note that we have scaled the $L^2$-error in $\underline{\sigma }$ by $\nu ^{-1}$ so that this error can be directly compared to the $H^1$-errors in $\mathbf{u}-\mathbf{u}_h$ and $\mathbf{u}-\mathbb{P}\mathbf{u}_h$. These three errors are all of the same magnitude, with a slight advantage for the post-processed solution.

In Table 2 we show the seminorm of the errors $\mathbf{u}-\mathbf{u}_h$ and $\mathbf{u}-{\mathbb{P}}\mathbf{u}_h$ which measures their jumps. Note that it superconverges with order $3/2$. This means that the relative contribution of this seminorm to the $\|\cdot \|_{1,h}$-norm diminishes as $h$ decreases. An analysis of this phenomenon remains to be carried out.

In Table 3 we show the $L^2$-errors in the velocities and their corresponding convergence orders. In the first column we observe that the velocities converge with second-order accuracy. In the second column, we notice that by post-processing the error is reduced by a factor of roughly 3/2. Therefore, the post-processed solution should be used as the best approximation obtained by our scheme. Furthermore, we show the $L^\infty$-norms of the divergence of $\mathbb{P} \mathbf{u}_h$ (evaluated at the points of a 4-by-4 Gauss formula on each cell). These are of the order of the residual of the nonlinear iteration, confirming that the post-processed solution is indeed divergence free.

The convergence of the nonlinear iteration under consideration is illustrated in Table 4. The number of steps required to reduce the start residual by a factor of $10^{7}$ is displayed. The initial guess for the iterations is the vector $\mathbf{u}_h^0 = \mathbf{0}$.

The linear system in each step is solved by a preconditioned GMRES method up to a relative accuracy of $10^{-4}$, using the preconditioners described in 12. Therefore, the error of the linear iterations is small enough to be neglected. Table 4 shows that the number of iteration steps is not only bounded independently of the mesh-size, but in fact decreasing. This is in perfect agreement with our theoretical results in Theorem 4.7. If we decrease the viscosity, the increase of the number of iteration steps for convergence is quite moderate on fine grids. Of course, this only holds as long as there is convergence. With $\nu =10^{-3}$, the nonlinear iteration does not converge anymore, probably because the stationary solution is not stable in this case.

7. Concluding remarks

In this paper, we described and analyzed a new LDG method for the approximation of the two-dimensional incompressible Navier-Stokes equations on triangular meshes. The approximation is based on discontinuous $P_k-P_k-P_{k-1}$ elements for the approximation of the velocity, the velocity gradient, and the pressure, respectively, and on a post-processing procedure. Alternatively, the approximation of the velocity field can be based on div-conforming BDM elements for which this procedure can be omitted. These are the only DG methods that are locally conservative and stable for the incompressible Navier-Stokes equations. They are also the only methods that provide a systematic and simple way to obtain an approximate velocity that is exactly divergence free for polynomials of degree greater than or equal to 1; see 20.

As discussed in subsection 4.3, the results of this paper can be extended in a straightforward way to simplicial meshes in three dimensions and to $Q_k-Q_k-P_{k-1}$ elements on quadrilateral or hexahedral affine meshes. Future work will be devoted to extending the methods to the case in which the numerical flux $\widehat{\mathbf{u}}_h^p$ depends also on the pressure. This happens, for example, when all the unknowns are approximated with the same polynomial space; see 9, 8. Indeed, in this case, the operator $\mathbb{P}$ depends not only on the velocity but also on the pressure. As a consequence, the analysis of the convergence of the fixed point iteration is much more delicate. It will be carried out in a forthcoming paper.

Table 1.

Errors and orders of convergence for $\nu =0.1$.

$L$	$\|p-p_h\|_0$		$\|\mathbf{u}-\mathbf{u}_h\|_{1,h}$		$\|\mathbf{u}-\mathbb{P} \mathbf{u}_h\|_{1,h}$		$\nu ^{-1}\|\underline{\sigma }-\underline{\sigma }_h\|_0$
3	2.2e+0	—	1.2e+1	—	8.1e+0	—	7.0e-0	—
4	1.0e+0	1.12	5.4e+0	1.11	3.2e+0	1.33	3.4e-0	1.05
5	4.8e-1	1.10	2.4e+0	1.16	1.4e+0	1.18	1.6e-0	1.07
6	2.3e-1	1.04	1.1e+0	1.18	6.8e-1	1.06	7.8e-1	1.04
7	1.2e-1	1.01	4.7e-1	1.17	3.4e-1	1.02	3.9e-1	1.02
8	5.8e-2	1.00	2.2e-1	1.13	1.7e-1	1.01	1.9e-1	1.02

Table 2.

Errors and orders of convergence for $\nu =0.1$ in the seminorm $|\mathbf{v}|_h:=\left\{\sum _{e\in {\mathcal{E}}_h}\int _e \, \kappa _0 {\mathtt{h}}^{-1} \,|[\![\mathbf{v}\otimes \mathbf{n}]\!]|^2\, ds\right\}^{1/2}$.

$L$	$|\,\mathbf{u}-\mathbf{u}_h\,|_{h}$		$|\,\mathbf{u}-\mathbb{P} \mathbf{u}_h\,|_{h}$
3	9.1e+0	—	4.8e+0	—
4	4.2e+0	1.11	1.5e+0	1.72
5	1.8e+0	1.23	4.7e-1	1.64
6	7.2e-1	1.32	1.6e-1	1.54
7	2.8e-1	1.39	5.6e-2	1.52
8	1.0e-1	1.44	2.0e-2	1.51

Table 3.

$L^2$-errors and orders of convergence in the velocity and $L^\infty$-norm of the divergence of the post-processed solution $\mathbb{P} \mathbf{u}_h$ for $\nu =0.1$.

$L$	$\|\mathbf{u}-\mathbf{u}_h\|_0$		$\|\mathbf{u}-\mathbb{P} \mathbf{u}_h\|_0$		$\| \nabla \cdot \mathbb{P} \mathbf{u}_h\|_\infty$
3	6.4e-1	—	4.9e-1	—	1.4e-12
4	1.6e-1	2.03	1.1e-1	2.22	1.4e-12
5	3.3e-2	2.22	2.0e-2	2.37	3.2e-12
6	7.1e-3	2.24	4.2e-3	2.27	1.5e-11
7	1.6e-3	2.19	9.8e-4	2.12	1.8e-12
8	3.5e-4	2.13	2.4e-4	2.04	2.9e-11

Table 4.

Number of iterations for convergence of the nonlinear iteration.

$L$	$\nu =1$	$\nu =0.1$	$\nu = 0.01$
3	14	33	865
4	10	21	106
5	8	16	54
6	6	12	29
7	5	10	18
8	6	10	13