Error Analysis of Non Inf-sup Stable Discretizations of the time-dependent Navier–Stokes Equations with Local Projection Stabilization

# Error Analysis of Non Inf-sup Stable Discretizations of the time-dependent Navier–Stokes Equations with Local Projection Stabilization

Javier de Frutos Instituto de Investigación en Matemáticas (IMUVA), Universidad de Valladolid, Spain. Research supported by Spanish MINECO under grants MTM2013-42538-P (MINECO, ES) and MTM2016-78995-P (AEI/FEDER, UE). (frutos@mac.uva.es)    Bosco García-Archilla Departamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research supported by Spanish MINECO under grant MTM2015-65608-P (bosco@esi.us.es)    Volker John Weierstrass Institute for Applied Analysis and Stochastics, Leibniz Institute in Forschungsverbund Berlin e. V. (WIAS), Mohrenstr. 39, 10117 Berlin, Germany. Freie Universität Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany.    Julia Novo Departamento de Matemáticas, Universidad Autónoma de Madrid, Spain. Research supported by Spanish MINECO under grants MTM2013-42538-P (MINECO, ES) and MTM2016-78995-P (AEI/FEDER, UE) (julia.novo@uam.es)
July 26, 2019
###### Abstract

This paper studies non inf-sup stable finite element approximations to the evolutionary Navier–Stokes equations. Several local projection stabilization (LPS) methods corresponding to different stabilization terms are analyzed, thereby separately studying the effects of the different stabilization terms. Error estimates are derived in which the constants in the error bounds are independent of inverse powers of the viscosity. For one of the methods, using velocity and pressure finite elements of degree , it will be proved that the velocity error in decays with rate in the case that , with being the dimensionless viscosity and the mesh width. In the analysis of another method, it was observed that the convective term can be bounded in an optimal way with the LPS stabilization of the pressure gradient. Numerical studies confirm the analytical results.

## 1 Introduction

Let , , be a bounded domain with polyhedral and Lipschitz boundary . The incompressible Navier–Stokes equations model the conservation of linear momentum and the conservation of mass (continuity equation) by

 ∂tu−νΔu+(u⋅∇)u+∇p =f in  (0,T]×Ω, ∇⋅u =0 in  (0,T]×Ω, (1) u(0,⋅) =u0(⋅) in  Ω,

where is the velocity field, the kinematic pressure, the kinematic viscosity coefficient, a given initial velocity, and represents the external body accelerations acting on the fluid. The Navier–Stokes equations (1) are equipped with homogeneous Dirichlet boundary conditions on .

This paper studies approximations to the Navier–Stokes equations (1) with non inf-sup stable mixed finite elements in space and the implicit Euler method in time. We use the so-called local projection stabilization (LPS) method to stabilize the pressure (since non inf-sup stable elements are used) plus other stabilization terms which aim at allowing to derive error estimates where the constants do not depend explicitly on inverse powers of the viscosity but only implicitly through norms of the solution of (1). This kind of bounds are called semi-robust or quasi-robust in the literature, see for example [4].

In the literature, one can find already investigations of LPS methods for approximating the solution of (1). LPS methods for inf-sup stable elements are analyzed in [3]. The derived error bounds depend explicitly on inverse powers of the viscosity parameter , unless the grids are becoming sufficiently fine (, where is the mesh width), see also [13] where error bounds for the Oberbeck-Boussinesq model model are obtained with an assumption on the regularity of the finite element solution. In [2], the authors consider non inf-sup stable mixed finite elements with LPS stabilization. The so called term-by-term stabilization is applied, see [11]. This method is a particular type of a LPS method that is based on continuous functions, it does not need enriched finite element spaces, and an interpolation operator replaces the standard projection operator of the classical LPS methods. As in the present paper, a fully discrete scheme with the implicit Euler method as time integrator is considered. A fully discrete LPS method for inf-sup stable pairs of finite element spaces and a pressure-projection scheme is analyzed in [4].

Our analysis starts as in [2], but there are several major differences in the formulation of the discrete problem as well as in the obtained results. First of all, as an important result which was not achieved in [2], we are able to derive error bounds in which the constants do not depend on inverse powers of the diffusion parameter. Also, contrary to [2], where only one method is analyzed (with LPS stabilizations of the pressure, the divergence, and the convective term), we consider several methods, because our aim is to study separately the effects of the different stabilization terms. For all of them, error bounds with constants independent on inverse powers of the diffusion parameter are achieved with the smallest possible number of stabilization terms. Also, in contrast to [2], only moderate assumptions on the smallness of the time step are needed, like in the error analysis of the pressure, while in [2] the smallness assumption on the mesh width is required.

Section 3 considers a method with LPS stabilization for the pressure and a global grad-div stabilization term. The global grad-div stabilization term was proposed to reduce the violation of mass conservation of finite element methods, but there are already investigations which show that this term also stabilizes dominant convection. In [14], semi-robust error estimates are proved for the standard Galerkin method plus grad-div stabilization in the case of inf-sup stable elements, both for the continuous-in-time case and for the fully discrete case. Paper [14] considers both, the regular case and the situation in which nonlocal compatibility conditions for the solution are not assumed. The results of Section 3 can be seen as an extension of some of the results from [14] to the case of non inf-sup stable elements and also as an improvement of the results from [2]. Error bounds of order are obtained for a sufficiently smooth solution, where , being the regularity index of the solution and being the degree of the polynomials used. The error is bounded in a norm that includes the norm of the velocity at the final time step and the norm of the divergence. This rate of convergence is the same as obtained in [14] for a similar norm and also the same rate as proved in [2]. However, as we pointed out above, in [2] more terms are included in the method, the bound depends explicitly on , and the restriction is assumed. For the error bound of the pressure, we get the optimal order . However, following the ideas of [2], we are able to bound the error of the norm of a discrete in time primitive of the pressure instead of the stronger discrete in time norm of the pressure. Although Section 3 studies the term-by-term stabilization, the analysis also holds for the standard one-level LPS method, see [15, 21], with slight modifications.

In Section 4, we analyze a method with LPS stabilization for the pressure and LPS stabilization with control of the fluctuations of the gradient. For this section, the use of term-by-term stabilization is necessary since in the error analysis we need to have the same polynomial spaces for the velocity and the pressure. A key ingredient in the error analysis is the application of [8, Theorem 2.2]. This result was already applied in the error analysis in [10], where the authors proved semi-robust error bounds for the evolutionary Navier–Stokes equations and a continuous interior penalty (CIP) method in space assuming enough regularity of the solution. For the method studied in Section 4, the convective term is estimated in an optimal way (with constants independent on inverse powers of the diffusion parameter) with the help of the LPS stabilization of the gradient of the pressure. This LPS term was introduced in [6] to account for the violation of the discrete inf-sup condition by the used pair of finite elements.

Following the analysis of the previous section, Section 5 presents analogous error bounds for a method with both LPS stabilization for the pressure and the divergence.

For the methods analyzed in Sections 35, error estimates with constants independent on inverse powers of the diffusion parameter are derived with the help of stabilization terms that were not proposed for stabilizing dominant convection but to account for the non-satisfaction of the discrete inf-sup condition or the violation of the mass conservation (note that the LPS term of the velocity gradient of the method from Section 4 was not utilized for estimating the convective term). The deeper reasons for this behavior are not yet understood and their explanation is formulated as an open problem in [19].

In Section 6, it is shown that the rate of decay of the velocity error in the situation can be improved for the method from Section 4 by choosing different values of the stabilization parameters and increasing the regularity assumption for the pressure. Concretely, a bound of order is proved for an error which contains the error of the velocity. This is the same order that was obtained for the CIP method in [10] under the same regularity assumptions. We are not aware of any other paper where this order is proved and it is still an open question whether the optimal expected order for the error of the velocity can be achieved or not, see [19].

Finally, Section 7 presents numerical studies that confirm the analytical results.

## 2 Preliminaries and notation

Throughout the paper, will denote the Sobolev space of real-valued functions defined on the domain with distributional derivatives of order up to in . These spaces are endowed with the usual norm denoted by . If is not a positive integer, is defined by interpolation [1]. In the case it is . As it is standard, will be endowed with the product norm and, since no confusion can arise, it will be denoted again by . The case will be distinguished by using to denote the space . The space is the closure in of the set of infinitely differentiable functions with compact support in . For simplicity, (resp. ) is used to denote the norm (resp. semi norm) both in or . The exact meaning will be clear by the context. The inner product of or will be denoted by and the corresponding norm by . For vector-valued functions, the same conventions will be used as before. The norm of the dual space of is denoted by . As usual, is always identified with its dual, so one has with compact injection. The following Sobolev’s embedding [1] will be used in the analysis: For let be such that . There exists a positive constant , independent of , such that

 ∥v∥Lq′(Ω)≤C∥v∥Ws,p(Ω),1q′≥1q,v∈Ws,p(Ω). (2)

If the above relation is valid for . A similar embedding inequality holds for vector-valued functions.

Using the function spaces

 V=H10(Ω)d,Q=L20(Ω)={q∈L2(Ω):(q,1)=0},

the weak formulation of problem (1) is as follows: Find such that for all ,

 (∂tu,v)+ν(∇u,∇v)+((u⋅∇)u,v)−(∇⋅v,p)+(∇⋅u,q)=(f,v), (3)

and .

The Hilbert space

 Hdiv={u∈L2(Ω)d ∣ L2(Ω)∋∇⋅u=0,u⋅n|∂Ω=0}

will be endowed with the inner product of and the space

 Vdiv={u∈V ∣ ∇⋅u=0}

with the inner product of .

In the error analysis, the Poincaré–Friedrichs inequality

 ∥v∥0≤CPF∥∇v∥0∀v∈V (4)

will be used.

## 3 Local projection stabilization with global grad-div stabilization.

Let be a family of triangulations of . Given an integer and a mesh cell we denote by the space of polynomials of degree less or equal to . We consider the following finite element spaces

 Ylh = {vh∈C0(¯¯¯¯Ω)∣vh∣K∈Pl(K),∀K∈Th}, l≥1, Ylh = (Ylh)d,Xh=Ylh∩(H10)d, Qh = Ylh∩L20.

It will be assumed that the family of meshes is quasi-uniform and that the following inverse inequality holds for each , e.g., see [12, Theorem 3.2.6],

 ∥vh∥Wm,p(K)≤Cinvhn−m−d(1q−1p)K∥vh∥Wn,q(K), (5)

where , , and is the size (diameter) of the mesh cell .

We consider the approximation of (1) with the implicit Euler method in time and a LPS method with grad-div stabilization in space. Given , find such that

 (un+1h−unhΔt,vh)+ν(∇un+1h,∇vh)+b(un+1h,un+1h,vh)−(pn+1h,∇⋅vh) +Sh(un+1h,vh) = (fn+1,vh)∀vh∈Xh, (6) (∇⋅un+1h,qh)+spres(pn+1h,qh) = 0∀qh∈Qh,

where

 Sh(u,v) = μ(∇⋅u,∇⋅v), b(u,v,w) = (B(u,v),w)∀u,v,w∈H10(Ω)d, B(u,v) = (u⋅∇)v+12(∇⋅u)v∀u,v∈H10(Ω)d, spres(pn+1h,qh) = ∑K∈Thτp,K(σ∗h(∇pn+1h),σ∗h(∇qh))K,

and and are the grad-div and pressure stabilization parameters, respectively. In addition, , where is a locally stable projection or interpolation operator from on , that is, there exists a constant such that for any

 ∥σjh(v)∥L2(K)≤C∥v∥L2(ωK),∀v∈L2(Ω)d, (7)

where is the union of all mesh cells whose intersection with is not empty. It will be assumed that the number of mesh cells in each set is bounded independently of the triangulation and of . From (7), also the stability of follows. The operator can be chosen as a Bernardi–Girault [7], Girault–Lions [16], or the Scott–Zhang [22] interpolation operator in the space (for a proof of (7) in the case of the last two operators see [9]). The following bound holds for ,

 ∥v−σjh(v)∥L2(K)≤ChsK∥v∥Hs(ωK),1≤s≤j+1 (8)

from which it can be deduced that

 ∥v−σjh(v)∥0≤Chs|v|s,1≤s≤j+1 (9)

see [22, 7, 9]. Bounds (8) and (9) will be applied for .

In the sequel, we will assume that

 α1h2K≤τp,K≤α2h2K (10)

for some positive constants independent of . In addition, the notations

 (f,g)τp=∑K∈Thτp,K(f,g)K,∥f∥τp=(f,f)1/2τp (11)

are used.

The following inf-sup condition holds (see [2, Lemma 4.2]).

###### Lemma 1

The following inf-sup condition holds

 ∥qh∥0≤β0(supvh∈Xh(∇⋅vh,qh)∥∇vh∥0+∥σ∗h(∇qh)∥τp)∀qh∈Qh.

Along the paper we will use the following discrete Gronwall inequality whose proof can be found in [17].

###### Lemma 2

Let be nonnegative numbers such that

 aj+kn∑j=0bj≤kn∑j=0γjaj+kn∑j=0cj+B,forn≥0.

Suppose that , for all j, and set . Then

 aj+kn∑j=0bj≤exp(kn∑j=0σjγj){kn∑j=0cj+B},forn≥0.

### 3.1 Error bound for the velocity

Let us denote by and by . Following [2], we consider an approximation satisfying

 (un−^unh,vh)=0,∀vh∈Yl−1h,n=0,1,…,N. (12)

Let us observe that the above definition for can be applied for any time so that we can consider that is continuous in the variable. The following bound holds, see [2],

 ∥un−^unh∥Wm,p≤Chs+1−m+d/p−d/2|un|s+1,n=0,1,…,N, (13)

for , , .

Let with being the standard interpolation operator. There exists a constant such that

 ∥pn−^pnh∥Wm,p≤Chs−m+d/p−d/2|pn|s,n=0,1,…,N,m=0,1, (14)

see [9].

Let us denote

 ^enh=^unh−un,enh=^unh−unh,^λnh=^pnh−pn,λnh=^pnh−pnh. (15)

Subtracting the discrete problem (3) from the continuous problem (3) yields the error equation

 (en+1h−enhΔt,vh)+ν(∇en+1h,∇vh)+b(^un+1h,^un+1h,vh)−b(un+1h,un+1h,vh) −(λn+1h,∇⋅vh)+(∇⋅en+1h,qh)+spres(λn+1h,qh)+Sh(en+1h,vh) = (ξn+1vh,vh)+(ξn+1qh,qh)+ν(∇^en+1h,∇vh)+spres(^pn+1h,qh) +Sh(^un+1h,vh)−(^λn+1h,∇⋅vh),

for all and . In (3.1), and are defined as follows

 ξn+1vh = ξn+1vh,1+ξn+1vh,2, (17) (ξn+1vh,1,vh) = −(∂tun+1−^un+1h−^unhΔt,vh), (18) (ξn+1vh,2,vh) = −b(un+1,un+1,vh)+b(^un+1h,^un+1h,vh), (19) (ξn+1qh,qh) = (∇⋅^en+1h,qh).

Remark 1  Note that the error equation (3.1) holds even for with . Let and denote by the mean of , then (3.1) gives

 (∇⋅en+1h,qh−m(qh))+spres(λn+1h,qh−m(qh)) = (∇⋅^en+1h,qh−m(qh))+spres(^pn+1h,qh−m(qh)).

Since the terms , , , and vanish, it follows that

 (∇⋅en+1h,qh)+spres(λn+1h,qh)=(∇⋅^en+1h,qh)+spres(^pn+1h,qh)∀ qh∈Ylh.

Setting , rearranging terms, and using the Cauchy–Schwarz inequality and Young’s inequality gives

 ∥en+1h∥202Δt−∥enh∥202Δt+∥en+1h−enh∥202Δt+ν2∥∇en+1h∥20+∥σ∗h(∇λn+1h)∥2τp +Sh(en+1h,en+1h) ≤ ∣∣b(un+1h,un+1h,en+1h)−b(^un+1h,^un+1h,en+1h)∣∣+∥ξn+1vh∥202+∥en+1h∥202 +∣∣Sh(^un+1h,en+1h)∣∣+∣∣(^λn+1h,∇⋅eh)∣∣.

Now, the terms on the right-hand sid of (3.1) will be bounded. We start with the last two terms. Applying the Cauchy–Schwarz inequality, Young’s inequality, and (13) yields

 ∣∣Sh(^un+1h,en+1h)∣∣ =μ∣∣(∇⋅^un+1h,∇⋅en+1h)∣∣≤μ8∥∇⋅en+1h∥20+2μ∥∇⋅^en+1h∥20 ≤18Sh(en+1hen+1h)+Cμh2s∥u∥2L∞(Hs+1). (21)

Similarly, we obtain

 ∣∣(^λn+1h,∇⋅eh)∣∣≤μ8∥∇⋅en+1h∥20+2μ∥^λn+1h∥20≤18Sh(en+1h,en+1h)+Cμh2s∥p∥2L∞(Hs), (22)

where in the last inequality (14) was applied. The nonlinear term in (3.1) can be bounded as in [14] using the skew-symmetric property of

 |b(un+1h,un+1h,en+1h)−b(^un+1h,^un+1h,en+1h)| (23) ≤ |b(en+1h,^un+1h,en+1h)|+|b(un+1h,en+1h,en+1h)| ≤ ∥∇^un+1h∥L∞∥en+1h∥20+12∥∇⋅en+1h∥0∥^un+1h∥L∞∥en+1h∥0 ≤ (∥∇^un+1h∥L∞+∥^un+1h∥2L∞4μ)∥en+1h∥20+μ4∥∇⋅en+1h∥20.

For the fourth term on the right-hand side of (3.1), integrating by parts and using (12), (10), and (13) gives

 |(ξn+1qh,λn+1h)| = |(^en+1h,∇λn+1h)|=|(^en+1h,σ∗h(∇λn+1h))| (24) ≤ ∑K∈Th∥^en+1h∥2L2(K)τp,K+14∥σ∗h(∇λn+1h))∥2τp ≤ Ch2s∥u∥2L∞(Hs+1)+14∥σ∗h(∇λn+1h))∥2τp.

For the fifth term, we use (13) to get

 ν2∥∇^en+1h∥20≤Cνh2s∥u∥2L∞(Hs+1). (25)

To bound the sixth term, the usual inequalities, the definition (11) of , (9), and (14) are utilized

 ∣∣spres(^pn+1h,λn+1h)∣∣ ≤ ∥σ∗h(∇^pn+1h)∥2τp+14∥σ∗h(∇λn+1h)∥2τp (26) ≤ 2∥σ∗h(∇^λn+1h)∥2τp+2∥σ∗h(∇pn+1)∥2τp+14∥σ∗h(∇λn+1h)∥τp ≤ Ch2∥∇^λn+1h∥20+Ch2∥σ∗h(∇pn+1)∥20+14∥σ∗h(∇λn+1h)∥τp ≤ Ch2s∥p∥2L∞(Hs)+14∥σ∗h(∇λn+1h)∥τp.

Inserting now (3.1) – (26) in (3.1) yields

 ∥en+1h∥20−∥enh∥20+Δtν∥∇en+1h∥20+Δt∥σ∗h(∇λn+1h)∥2τp+μΔt∥∇⋅en+1h∥20 ≤ Δt(1+2∥∇^un+1h∥L∞+∥^un+1h∥2L∞2μ)∥en+1h∥20+Δt∥ξn+1vh∥20 +CΔth2s((1+ν+μ)∥u∥2L∞(Hs+1)+(1+μ−1)∥p∥2L∞(Hs)),

such that summing over the discrete times leads to

 ∥enh∥20+Δtνn∑j=1∥∇ejh∥20+Δtn∑j=1∥σ∗h(∇λjh)∥2τp+Δtμn∑j=1∥∇⋅ejh∥20 ≤ ∥e0h∥20+n∑j=1Δt(1+2∥∇^ujh∥L∞+∥^ujh∥2L∞2μ)∥ejh∥20+Δtn∑j=1∥ξjvh∥20 +CTh2s((1+ν+μ)∥u∥2L∞(Hs+1)+(1+μ−1)∥p∥2L∞(Hs)).

Let us bound and , . For the first term, applying (2) and (13) we have

 ∥^ujh∥L∞ ≤ ∥uj∥L∞+∥uj−^ujh∥L∞≤C∥uj∥2+Ch2−d/2∥uj∥2 (27) ≤ C∥u∥L∞(H2).

Using the same argument for the second term, we reach

 ∥∇^ujh∥L∞ ≤ ∥∇uj∥L∞+∥∇uj−∇^ujh∥L∞≤C∥uj∥3+Ch2−d/2∥uj∥3 (28) ≤ C∥u∥L∞(H3).

From (27) and (28) we deduce

 1+2∥∇^ujh∥L∞+∥^ujh∥2L∞2μ≤^Mu,^Mu=1+C⎛⎝2∥u∥L∞(H3)+∥u∥2L∞(H2)2μ⎞⎠. (29)

Let us assume

 Δt^Mu≤12. (30)

Applying the Gronwall lemma, Lemma 2, we get

 ∥enh∥20+Δtνn∑j=1∥∇ejh∥20+Δtn∑j=1∥σ∗h(∇λjh)∥2τp+Δtμn∑j=1∥∇⋅ejh∥20 ≤ e2T^Mu(∥e0h∥20+Δtn∑j=1∥ξjvh∥20) +Ce2T^Mu(Th2s((1+ν+μ)∥u∥2L∞(Hs+1)+(1+μ−1)∥p∥2L∞