Volume 2010, Article ID 419021,21pages doi:10.1155/2010/419021

*Research Article*

**Discontinuous Time Relaxation Method for** **the Time-Dependent Navier-Stokes Equations**

**Monika Neda**

*Department of Mathematical Sciences, University of Nevada Las Vegas, Las Vegas, NV 89154-4020, USA*

Correspondence should be addressed to Monika Neda,monika.neda@unlv.edu Received 17 July 2010; Accepted 16 September 2010

Academic Editor: William John Layton

Copyrightq2010 Monika Neda. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

A high-order family of time relaxation models based on approximate deconvolution is considered.

A fully discrete scheme using discontinuous finite elements is proposed and analyzed. Optimal
velocity error estimates are derived. The dependence of these estimates with respect to the
Reynolds number Re isORe*e*^{Re}, which is an improvement with respect to the continuous finite
element method where the dependence isORe*e*^{Re}^{3}.

**1. Introduction**

Turbulence is a phenomenon that appears in many processes in the nature, and it is connected
with many industrial applications. Based on the Kolmogorov theory1, Direct Numerical
Simulation DNS of turbulent flow, where all the scales/structures are captured, requires
the number of mesh points in space per each time step to beORe^{9/4}in three-dimensional
problems, where Re is the Reynolds number. This is not computationally economical and
sometimes not even feasible. One approach is to regularize the flow and one such type of
regularization is the time relaxation method, where an additional term, the so-called time
relaxation term, is added to the Navier-Stokes equationscf. Adams and Stolz2and Layton
and Neda3. The contribution to the Navier-Stokes equations from the time relaxation term
induces an action on the small scales of the flow in which these scales are driven to zero. This
time relaxation term is based on filtering and deconvolution methodology. In general, many
spacial filtering operators associated with a length-scale*δ*are possiblecf. Berselli et al.4,
John5, Geurts6, Sagaut7and Germano8. First, consider the equations of diﬀerential
filtercf. Germano8

*φGφ* in Ω, 1.1

*φ***0 on** *∂Ω,* 1.2

where*G*^{−1} −δ^{2}Δ *I,*Ωis the domain and *∂Ω*is its boundary. Here*δ >* 0 represents the
averaging radius, in general, chosen to be of the order of the mesh size.

The deconvolution algorithm that it is considered herein was studied by van Cittert in
1931, and its use in Large Eddy SimulationLESpioneered by Stolz and Adamscf. Stolz
and Adams2,9. For each*N*0,1, . . . ,it computes an approximate solution**u***N**by N steps*
of a fixed point iteration for the fixed point problemcf. Bertero and Boccacci10

given**u solve uu**{u−*Gu},* for**u.** 1.3

The deconvolution approximation is then computed as follows.

*Algorithm 1.1*van Cittert approximate deconvolution algorithm. Consider that**u**0 **u, for***n*1,2, . . . , N−1, perform

**u**_{n1}**u***n*{u−*Gu**n*}. 1.4

By eliminating the intermediate steps, it is easy to find an explicit formula for the*Nth*
deconvolution operator*G** _{N}*given by

*G*_{N}**u :**^{N}

*n0*

I−*G*^{n}**u.** 1.5

For example, the approximate deconvolution operator corresponding to*N*0,1,2 is
*G*0**uu,**

*G*_{1}**u**2u−**u,**
*G*_{2}**u**3u−3u**u.**

1.6

Many fluid models that are based on numerical regularization and computational
stabilizations have been explored in computational fluid dynamics. One such regularization
and the most recent has been proposed by Stolz et al.11,12and arises by adding a linear,
lower order time regularization term, *χu*^{} *χu*−*G*_{N}**u, to the Navier-Stokes equations**
NSE. This term involves **u**^{} which represents the part of the velocity that fluctuates on
scales less than order*δ* and it is added to the NSE with the aim of driving the unresolved
fluctuations of the velocity field to zero. The time relaxation family of models, under the no-
slip boundary condition, is then defined by

**u***t***u**· ∇u−*νΔu*∇p*χu*−*G*_{N}**u f in**0, T×Ω, 1.7

∇ ·**u**0 in 0, T×Ω, 1.8

**u0 in** 0, T×*∂Ω,* 1.9

**u0,**· **u**0· inΩ, 1.10

whereΩ⊂R^{2}, is a convex bounded regular domain with boundary*∂Ω,***u is the fluid velocity,**
*p*is the fluid pressure and**f is the body force driving the flow. The kinematic viscosity***ν >*0
is inversely proportional to the Reynolds number of the flow. The initial velocity is given by
**u**0. A pressure normalization condition

Ω*p*0 is also needed for uniqueness of the pressure.

The time relaxation coeﬃcient*χ*has units 1/time. The domain is two-dimensional, but the
numerical methods and the analysis can be generalized to three-dimensional domains, as
stated in13for the case of Stokes and Navier-Stokes problems.

Existence, uniqueness and regularity of strong solutions of these models are discussed in3. Even though there are papers on the simulation of the models for incompressible and compressible flows, there is little published work in the literature on the numerical analysis of the models. In14, a fully discrete scheme using continuous finite elements and Crank- Nicolson for time discretization is analyzed and the energy cascade and joint helicity-energy cascades are studied in3,15, respectively.

In this work, a class of discontinuous finite element methods for solving high-
order time relaxation family of fluid models1.7–1.10is formulated and analyzed. The
approximations of the averaged velocity **u and pressure** *p* are discontinuous piecewise
polynomials of degree*r* and*r* −1, respectively. Because of the lack of continuity constraint
between elements, the Discontinuous Galerkin DG methods oﬀer several advantages
over the classical continuous finite element methods: i local mesh refinement and
derefinement are easily implemented several hanging nodes per edge are allowed; ii
the incompressibility condition is satisfied locally on each mesh element;iiiunstructured
meshes and domains with complicated geometries are easily handled. In the case of DNS,
DG methods have been applied to the steady-state NSEcf. Girault et al. 16and to the
time-dependent NSE cf. Girault et al. 17 where they are combined with an operator-
splitting technique. Another discontinuous Galerkin method for the NSE based on a mixed
formulation are considered in 18 by Cockburn et al.. For high Reynolds numbers, the
numerical analysis of a DG scheme combined with a large eddy simulation turbulence model
subgrid eddy viscosity modelis derived in19by Kaya and Rivi`ere.

This paper is organized in the following way. Section 2 introduces some notation
and mathematical properties. InSection 3, the fully discrete schemes are introduced and it
is proved that the schemes solutions are computable. A priori velocity error estimates are
derived inSection 4. The family of models1.7–1.10is regularization of the NSE. Thus, the
correct question is to study convergence of discretizations of1.7–1.10to solutions of the
NSE as*h*and *δ* → 0rather than to solution of 1.7–1.10. This is the problem studied
herein. Conclusions are given in the last section.

**2. Notation and Mathematical Preliminaries**

To obtain a discretization of the model a regular family of triangulationsE*h*ofΩ, consisting
of triangles of maximum diameter*h, is introduced. Leth** _{E}*denote the diameter of a triangle

*E*and

*ρ*

*the diameter of its inscribed circle. Regulary, it is meant that there exists a parameter*

_{E}*ζ >*0, independent of

*h, such that*

*h**E*

*ρ*_{E}*ζ** _{E}*≤

*ζ,*∀E∈ E

*h*

*.*2.1

This assumption will be used throughout this work.Γ*h*denotes the set of all interior edges of
E*h*. Let*e*denote a segment ofΓ*h*shared by two triangles*E** ^{k}*and

*E*

*k < lofE*

^{l}*h*; it is associated with

*e*a specific unit normal vector

**n**

*e*directed from

*E*

*to*

^{k}*E*

*and the jump and average of a function*

^{l}*φ*on

*e*is formally defined by

*φ*

*φ*

*E*^{k}

*e*−
*φ*

*E*^{l}

*e**,* *φ*

1 2

*φ*

*E*^{k}

*e*1
2

*φ*

*E*^{l}

*e**.* 2.2

If*e*belongs to the boundary*∂Ω, then***n***e*is the unit normal**n exterior to**Ωand the jump and
the average of*φ*on*e*coincide with the trace of*φ*on*e.*

Here, for any domainO,*L*^{2}Ois the classical space of square-integrable functions with
inner-productf, g_{O}

O*fg*and norm · _{0,O}. The space*L*^{2}_{0}Ωis the subspace of functions
of*L*^{2}Ωwith zero mean value

*L*^{2}_{0}Ω

**v**∈*L*^{2}Ω:

Ω**v**0

*.* 2.3

The standard Sobolev spaces are denoted by*W*_{p}* ^{r}*Ωwhere

*W*

_{2}

*Ωis the*

^{r}*H*

*Ω, with norm ·*

^{r}*and seminorm| · |*

_{r,Ω}*.*

_{r,Ω}Next, the discrete velocity and pressure spaces are defined to be consisting of
discontinuous piecewise polynomials. For any positive integer*r, the corresponding finite-*
dimensional spaces are

**X**^{h}

**v**∈

*L*^{2}Ω2

:∀E∈ E*h**,* **v**∈P*r*E^{2}

*,*
*Q*^{h}

*q*∈*L*^{2}_{0}Ω:∀E∈ E*h**, q*∈P*r−1*E
*,*

2.4

where*P**r*E span{x^{i}*y** ^{j}* :

*ij*≤

*r*}is defined as the span of polynomials of order

*r*over triangle

*E.*

Denoting by|e|the measure of*e, the following norms are associated for the spaces***X*** ^{h}*
and

*Q*

^{h} v _{X}

| ∇v |^{2}_{0,Ω}

*e∈Γ**h*∪∂Ω

1

|e| v ^{2}_{0,e}
_{1/2}

*,*
*q*

*Q* *q*

0,Ω*,*

2.5

where| v |_{0,Ω}is the broken norm defined by

| v |_{0,Ω}

*E∈E**h*

v ^{2}_{0,E}
1/2

*.* 2.6

Finally, some trace and inverse inequalities are recalled, that hold true on each element
*E*inE*h*, with diameter*h** _{E}*, the constant

*C*is independent of

*h*

_{E} v _{0,e}≤*C*

*h*^{−1/2}* _{E}* v

_{0,E}

*h*

^{1/2}

*∇v*

_{E}_{0,E}

*,* ∀e∈*∂E,* ∀v∈

*H*^{1}E_{2}

*,* 2.7

∇v _{0,e}≤*C*

*h*^{−1/2}* _{E}* ∇v

_{0,E}

*h*

^{1/2}

*∇*

_{E}^{2}

**v**

0,E

*,* ∀e∈*∂E,* ∀v∈

*H*^{2}E2

*,* 2.8

v _{0,e}≤*Ch*^{−1/2}* _{E}* v

_{0,E}

*,*∀e∈

*∂E,*∀v∈

**X**

^{h}*,*2.9 ∇v

_{0,e}≤

*Ch*

^{−1/2}

*∇v*

_{E}_{0,E}

*,*∀e∈

*∂E,*∀v∈

**X**

^{h}*.*2.10

**3. Numerical Methods**

In this section, the DG scheme is introduced and the existence of the numerical solution is
shown. First, the bilinear forms are defined*a*:**X*** ^{h}*×

**X**

*→ R, and*

^{h}*J*0:

**X**

*×*

^{h}**X**

*→ Rby*

^{h}*az,***v **

*E∈E**h*

*E*

∇z :∇v−

*e∈Γ**h*∪∂Ω

*e*

{∇z}n*e*·v *a*

*e∈Γ**h*∪∂Ω

*e*

{∇v}n*e*·z,

*J*0z,**v **

*e∈Γ**h*∪∂Ω

*σ*

|e|

*e*

z·v.

3.1

The parameter * _{a}* takes the value−1,0 or 1: this will yield diﬀerent schemes that are
slight variations of each other. It will be showed that all the resulting schemes are convergent
with optimal convergence rate in the energy norm ·

*. In the case where*

_{X}*a*−1, the bilinear form

*a*is symmetric; otherwise it is nonsymmetric. We remark that the form

*au,*

**v**is the standard primal DG discretization of the operator −Δu. Finally, if

*a*is either −1 or 0, the jump parameter

*σ*should be chosen suﬃciently large to obtain coercivity of

*a*see Lemma 3.1. If

*a*1, then the jump parameter

*σ*is taken equal to 1.

The incompressibility condition1.8 is enforced by means of the bilinear form*b* :
**X*** ^{h}*×

*Q*

*→ Rdefined by*

^{h}*b*
**v, q**

−

*E∈E**h*

*E*

*q∇ ·***v**

*e∈Γ**h*∪∂Ω

*e*

*q*

v·**n***e**.* 3.2

Finally, the DG discretization of the nonlinear convection term**w**· ∇w, which was introduced
in16by Girault et al. and studied extensively in16,17by the same authors, is recalled as
follows:

*c*** ^{z}**u; v,

**t**

*E∈E**h*

*E*

u· ∇v·**t**1
2

*E*

∇ ·**uv**·**t**

−1 2

*e∈Γ**h*∪∂Ω

*e*

u·**n***e*{v·**t}**

*E∈E**h*

*∂E*−

|{u} ·**n***E*|

**v**^{int}−**v**^{ext}

·**t**^{int}*,*

3.3

where

*∂E*_{−}{x∈*∂E*:{z} ·**n***E**<*0}, 3.4

the superscript**z denotes the dependence of***∂E*_{−}on**z and the superscript int**resp., extrefers
to the trace of the function on a side of*E*coming from the interior of*E*resp., coming from
the exterior of*E*on that side. When the side of*E*belongs to*∂Ω, the convention is the same*
as for defining jumps and average, that is, the jump and average coincide with the trace of
the function. Note that the form c is not linear with respect to**z, but linear with respect to u,v**
and**t.**

Some important properties satisfied by the forms*a, b, c*cf. Wheeler20, and Girault
et al.16,17are recalled.

**Lemma 3.1**Coercivity. If *a**1, assume thatσ1. If* *a*∈ {−1,0}, assume that*σis suﬃciently*
*large. Then, there is a constantκ >0, independent ofh, such that*

*av,***v ***J*0v,**v**≥*κ v *^{2}_{X}*,* ∀v∈**X**^{h}*.* 3.5

It is clear that*κ*1 if *a*1. Otherwise,*κ*is a constant that depends on the polynomial
degree of**v and of the smallest angle in the mesh. A precise lower bound for***σ*is given in21
by Epshteyn and Rivi`ere.

**Lemma 3.2**Inf-sup condition. There exists a positive constant*β, independent ofhsuch that*

*q∈Q*inf* ^{h}*sup

**v∈X**^{h}

*b*
**v, q**
v _{X}*q*

0,Ω

≥*β.* 3.6

**Lemma 3.3**Positivity. One has

*c*** ^{v}**v,

**z,z**≥0, ∀v,

**z**∈

**t**∈

*L*^{2}Ω_{2}

: **t|*** _{E}*∈

*H*^{2}E_{2}

∀E∈ E*h*

*.* 3.7

The discrete form of the diﬀerential filter1.1is defined following the work of Manica and Merdan22.

*Definition 3.4*Discrete diﬀerential filter. Given*v* ∈ *L*^{2}Ω, for a given filtering radius*δ >*

0, G* ^{h}*:

*L*

^{2}Ω →

**X**

*,*

^{h}**v**

^{h}*G*

^{h}**v where v**

*is the unique solution in*

^{h}**X**

*of*

^{h}*δ*^{2}
*a*

**v**^{h}*, φ*
*J*_{0}

**v**^{h}*, φ*

**v**^{h}*, φ*

**v, φ**

∀φ∈**X**^{h}*.* 3.8

*Remark 3.5. An attractive alternative is to define the diﬀerential filter by a discrete Stokes*
problem so as to preserve incompressibility approximately23–25.

*Definition 3.6. The discrete van Cittert deconvolution operatorsG*^{h}* _{N}*are

*G*^{h}_{N}**v :**^{N}

*n0*

Π*h*−*G**h*^{n}**v,** 3.9

whereΠ*h*:*L*^{2}Ω → **X*** ^{h}*is the

*L*

^{2}projection.

For**v**∈**X*** ^{h}*, the discrete deconvolution operator for

*N*0,1,2 is

*G*

^{h}_{0}

**vv,**

*G*^{h}_{1}**v**2v−**v**^{h}*,*
*G*^{h}_{2}**v**3v−3v^{h}**v**^{h}^{h}*.*

3.10

*G** _{N}* was shown to be an

*Oδ*

^{2N2}approximate inverse to the filter operator

*G*in Lemma 2.1 of Dunca and Epshteyn26, recalled next.

**Lemma 3.7.** *G**N**is a bounded, self-adjoint positive operator.G**N**is anOδ*^{2N2}*asymptotic inverse*
*to the filterG. Specifically, for smoothφand asδ* → *0,*

*φG**N**φ* −1^{N1}*δ*^{2N2}Δ^{N1}*G*^{N1}*φ.* 3.11
Some basic facts about discrete diﬀerential filters and deconvolution operators are
presented next.

* Lemma 3.8. For v*∈

*L*

^{2}Ω, one has the following bounds for the discretely filtered and approximately

*deconvolved*

**v:**v^{h}

0,Ω≤ v _{0,Ω}*,* 3.12

G^{h}_{N}**v**^{h}

0,Ω≤*CN v *_{0,Ω}*.* 3.13

*Proof. The proof of*3.12follows from the standard finite element techniques applied on the
discretized equation3.8of the filter problem1.1. Pick*φ***v*** ^{h}*. Then, using coercivity result

*δ*^{2}*κ*v^{h}^{2}

*X*v^{h}^{2}

0,Ω≤
**v,v**^{h}

Ω≤ 1

2 v ^{2}_{0,Ω}1
2

v^{h}^{2}

0,Ω*.* 3.14

The term*δ*^{2}*κ v*^{h}^{2}* _{X}*is positive, so it will be dropped, which yields
1

2
v^{h}^{2}

0,Ω ≤ 1

2 v ^{2}_{0,Ω}*.* 3.15

Multiplying by 2 and taking the square root yields the estimate3.12. Equation3.13follows
immediately from3.12and the definition of*G*^{h}* _{N}*.

**Lemma 3.9. For smooth**φthe discrete approximate deconvolution operator satisfies*φ*−*G*^{h}_{N}*φ*^{h}

0,Ω ≤*C*_{1}*δ*^{2N2}*Gφ*

2N2,Ω*C*_{2}

*δh*^{r}*h*^{r1}*N1*

*n2*

*G*^{n}*φ*

*r1,Ω*

*.* 3.16

*Proof. The proof follows the same arguments as in* 27 for the case of continuous finite
element discretization of the filter problem. The error is decomposed in the following way:

*φ*−*G*^{h}_{N}*φ*^{h}

_{0,Ω}≤*φ*−*G*_{N}*φ*

0,Ω*G*_{N}*φ*−*G*^{h}_{N}*φ*

0,Ω

*G*^{h}_{N}*φ*−*G*^{h}_{N}*φ*^{h}

_{0,Ω}*.* 3.17

Lemma 3.7gives

*φ*−*G*_{N}*φ*

0,Ω ≤*Cδ*^{2N2}*φ*

2N2,Ω*.* 3.18

The standard discontinuous finite element bound for3.8is given bycf. Rivi`ere13
*φ*−*φ*^{h}

0,Ω≤*C*

*δh*^{r}*h*^{r1}*φ*

*r1,Ω**.* 3.19

Lemma 3.8gives for the third term in3.17that G^{h}_{N}*φ*−*G*^{h}_{N}*φ** ^{h}* ≤

*C φ*−

*φ*

*. Then, 3.19is applied.*

^{h}Now, it is left to bound the second term from 3.17. First, note that for *N* 0,
G0*φ** ^{h}*−

*G*

^{h}_{0}

*φ*

^{h}_{0,Ω}0. Based on Definition 3.6 of continuous and discrete deconvolution operators and their expansionsee3.10,

*G*

*is a polynomial of degree*

_{N}*N*in

*G*and

*G*

^{h}*in*

_{N}*G*

*as well. Thus, the second term in3.17can be written as*

^{h}G*N**φ*−*G*^{h}_{N}*φ*

0,Ω

*N*
*n0*

*α**n*

*G*^{n}*φ*−
*G*^{h}_{n}

*φ*
0,Ω

≤^{N}

*n0*

*α**n*G^{n}*φ*−
*G*^{h}_{n}

*φ*

0,Ω*.* 3.20

For*O1*coeﬃcients*α** _{n}*and for

*N*1, the result3.19gives Gφ−

*G*^{h}*φ*

0,Ω

*φ*−*φ*

*h*

0,Ω

≤*C*

*δh*^{r}*h*^{r1}*φ*

*r1,Ω*

*.*

3.21

For*N*2, the results3.19and3.12give
G2*φ*−

*G*^{h}_{2}
*φ*

0,Ω
*φ*−*φ*^{h}

*h**h*
0,Ω

≤
*φ*−*φ*

*h*

0,Ω

*φ*

*h*

−*φ*^{h}

*h**h*
0,Ω

≤
*φ*−*φ*

*h*

0,Ω

*φ*−*φ*^{h}

0,Ω

≤*C*

*δh*^{r}*h*^{r1}*φ*

*r1,Ω*

*φ*

_{r1,Ω}

*.*

3.22

Inductively,

G*N**φ*−*G*^{h}_{N}*φ*

0,Ω≤*C*

*δh*^{r}*h*^{r1}*N1*

*n2*

*G*^{n}*φ*

*r1,Ω*

*.* 3.23

The proof is completed by combining the derived bounds for the terms in3.17.

*Remark 3.10. There remains the question of uniform in* *δ* bound of the last term,

|G^{n}*φ|** _{r1,Ω}*, in3.16. This is a question about uniformregularity of an elliptic-elliptic singular
perturbation problem and some results are proven in28by Layton. To summarize, in the
periodic case it is very easy to show by Fourier series that for all

*k*

*Gφ** _{r1,Ω}*≤

*Cφ*

*and thus*

_{r1,Ω}*G*

^{n}*φ*

*≤*

_{r1,Ω}*Cφ*

_{r1,Ω}*.*3.24 The nonperiodic case can be more delicate. Suppose

*∂Ω*∈

*C*

*and*

^{r3}*φ*0 on

*∂Ω i.e.,φ*∈

*H*

_{0}

^{1}Ω∩

*H*

*Ω. Call*

^{r1}*Gφφ*so

*φ*satisfies

−δ^{2}Δφ*φφ* inΩ, *φ*0 on*∂Ω.* 3.25

Then it is known that*φ*∈*H** ^{r3}*Ω∩

*H*

_{0}

^{1}Ω, andΔφ0 on

*∂Ω. Further,*φ

*j,Ω*≤*Cφ*_{j,Ω}*,* *j*0,1,2. 3.26

So,3.24holds for*r* −1,0,1. It also holds for higher values of *r* provided additionally
Δ^{j}*φ*0 on*∂Ω*for 0≤*j*≤r1/2−1.

Now consider the second term*n* 2, that is,*G*^{2}*φ* *φ. We know from elliptic theory*
for*φ*∈*H** ^{r1}*Ω

*H*_{0}^{1}Ω, that*φ*∈*H** ^{r3}*Ω

*H*_{0}^{1}Ω,as noted above Δφ0 on*∂Ω*and

−δ^{2}Δφ*φφ* inΩ, *φ* Δφ0 on *∂Ω.* 3.27

Theorem 1.1 in28then implies, uniformly in*δ,*
*φ*

*j,Ω* ≤*Cφ*

*j,Ω**,* *j*0,1,2,3,4. 3.28

This extends directly to*G*^{n}*φ.*

ExtendingLemma 3.9, the following assumption will be made.

*Assumption DG1. The*|G* ^{n}*φ|

*terms in3.16are independent of*

_{r1,Ω}*δ*and

*φ*−

*G*

^{h}

_{N}*φ*

^{h}≤*C*1*δ*^{2N2}*φ*

2N2,Ω*C*2

*δh*^{r}*h*^{r1}*φ*

*r1,Ω**.* 3.29

The minimal conditions that are assumed throughout are that thediscretefilter and discrete deconvolution used satisfy the following consistency conditions of Stanculescu 29.

*Assumption DG2.* *G** ^{h}*andI−

*G*

^{h}

_{N}*G*

*are symmetric, positive definiteSP Doperators.*

^{h}These have been proven to hold for van Cittert deconvolutioncf. Stanculescu29,
Manica and Merdan22and Layton et al.27. For the DG method, the second assumption
restricts our parameter *a* −1 for the discretization of the filter problem3.8, so that the
bilinear form*a·,*·is symmetric.

The numerical scheme that uses discontinuous finite elements in space and backward
Euler in time is derived next. For this, letΔtdenote the time step such that*M* *T/Δt*is a
positive integer. Let 0 *t*_{0} *< t*_{1} *<* · · · *< t*_{M}*T* be a subdivision of the interval 0, T. The
function*φ*evaluated at the time*t** _{m}*is denoted by

*φ*

*. With the above forms, the fully discrete scheme is: findu*

_{m}

^{h}

_{n}*, p*

^{h}

_{n}*∈*

_{n≥0}**X**

*×*

^{h}*Q*

*such that:*

^{h}1 Δt

**u**^{h}* _{n1}*−

**u**

^{h}

_{n}*,*

**v**

Ω*ν*
*a*

**u**^{h}_{n1}*,***v**
*J*_{0}

**u**^{h}_{n1}*,***v**
*χ*

**u**^{h}* _{n1}*−

*G*

^{h}

_{N}**u**

^{h}

_{n1}

^{h}*,*

**v**

Ω

*c*^{u}^{h}^{n}

**u**^{h}* _{n}*;

**u**

^{h}

_{n1}*,*

**v**

*b*

**v, p**_{n1}^{h}

f*n1**,***v**_{Ω} ∀v∈**X**^{h}*,*

3.30

*b*

**u**^{h}_{n1}*, q*

0 ∀q∈*Q*^{h}*,* 3.31
**u**^{h}_{0}*,***v**

Ω u0*,***v**_{Ω} ∀v∈**X**^{h}*.* 3.32
*Remark 3.11. The time relaxation term can be treated explicitly such that the optimal accuracy*
and stability are obtained and this would make the scheme much easier to compute30.

The consistency result of the semidiscrete scheme is showed next.

**Lemma 3.12**Consistency. Letu, p*be the solution to*1.7–1.10, thenu, p*satisfies*

u*t**,***v**_{Ω}*νau,***v ***J*0u,**v ***χ*

**u**−*G*^{h}_{N}**u**^{h}*,***v**

Ω*c*** ^{u}**u; u,

**v**

*b*

**v, p**f,

**v**

_{Ω}

*E*

**u, p,f; v**

∀v∈**X**^{h}*,*

3.33

*b*
**u, q**

0 ∀q∈*Q*^{h}*,* 3.34
u0,**v**_{Ω} u0*,***v**_{Ω} ∀v∈**X**^{h}*,* 3.35

*where the consistency errorEu,***v ***χu*−*G*^{h}_{N}**u**^{h}*,***v**_{Ω}−*χu*−*G**N***u,v**_{Ω}*. Furthermore,*
*Eu,***v**≤*C*

*δh*^{r}*h*^{r1}

|u|* _{r1,Ω}* v

_{0,Ω}

*.*3.36

*Proof. Equations*3.34and3.35are clearly satisfied because of1.8,1.9, and1.10and
the regularity of**u. Next, we multiply**1.7by**v and integrate over one mesh element***E*

u*t**,***v*** _{E}* ∇ ·uu,

**v**

*−*

_{E}*νΔu,*

**v**

_{E}*χu*−

*G*

_{N}**u,v**

_{E}∇p,**v**

*E* f,**v**_{E}*.* 3.37

Summing over all elements*E*

*E∈E**h*

u*t**,***v*** _{E}* ∇ ·uu,

**v**

_{Ω}−

*ν*

*E∈E**h*

Δu,**v**_{E}*χu*−*G*_{N}**u,v**_{Ω}

∇p,**v**

Ω f,**v**_{Ω}*.* 3.38

By Green’s formula

−ν

*E∈E**h*

Δu,**v*** _{E}*−νΔu,

**v**

_{Ω}

*ν*

*E∈E**h*

∇u,∇v* _{E}*−

*ν*

*e∈Γ**h*∪∂Ω

*e*

∇un*e*·v. 3.39

The regularity of**u then yields**

−ν

*E∈E**h*

Δu,**v**_{E}*νau,***v ***J*_{0}u,**v.** 3.40

Note that Green’s formula yields

∇p,**v**

Ω*b*
**v, p**

*,* 3.41

and that the incompressibility condition with the regularity of**u yields**

u· ∇u,**v**_{Ω} *c*** ^{u}**u; u,

**v.**3.42

The final result is obtained by bounding the consistency error*Eu,***v. For***N*0, we have
*Eu,***v**≤u−**u**^{h}

0,Ω v _{0,Ω} ≤*C*

*δh*^{r}*h*^{r1}

|u|* _{r1,Ω}* v

_{0,Ω}

*.*3.43

For*N*1,

*Eu,***v**≤Gu−G* ^{h}*u

^{h}0,Ω v _{0,Ω}

≤

**u**−**u**^{h}^{h}

0,Ω v _{0,Ω}

≤

**u**−**u*** ^{h}*
0,Ω

**u*** ^{h}*−

**u**

^{h}*0,Ω*

^{h}
v _{0,Ω}

≤

**u**−**u**^{h}

0,Ωu−**u**^{h}

0,Ω

v _{0,Ω}

≤*C*

*δh*^{r}*h** ^{r1}*u

*r1,Ω*|u|_{r1,Ω}

v _{0,Ω}*.*

3.44

The bound for*Eu,***v**is obtained by applying an induction argument andRemark 3.10.

The existence and uniqueness of the discrete solution is stated next.

**Proposition 3.13. Assume that**Lemma 3.1*holds. Then, there exists a unique solution to* 3.30–

3.32.

*Proof. The existence of***u**^{h}_{0} is trivial. Given**u**^{h}* _{n}*, the problem of finding a unique

**u**

^{h}*satisfying 3.30-3.31is linear and finite-dimensional. Therefore, it suﬃces to show uniqueness of the solution. Consider the problem restricted to the subspace*

_{n1}**V**

*defined by*

^{h}**V**^{h}

**v**∈**X*** ^{h}*:

*b*

**v, q**

0∀q∈*Q*^{h}

*.* 3.45

Let**u**^{h}* _{n1}*and

**u**

^{h}*be two solutions and let*

_{n1}

**θ***n1*

**u**

^{h}*−*

_{n1}**u**

^{h}*. Then,*

_{n1}

**θ***n1*satisfies:

1

Δtθ_{n1}*,***v**_{Ω}*νaθ*_{n1}*,***v ***J*_{0}θ_{n1}*,***v ***c*^{u}^{h}^{n}

**u**^{h}* _{n}*;

**θ**

_{n1}*,*

**v**

*χ*

**θ***n1*−*G*^{h}_{N}**θ***n1**h*

*,***v**

0 ∀v∈**V**^{h}*.*

3.46

Choosing**v** **θ*** _{n1}* and using the coercivity result3.5, positivity result3.7and positivity
of the operator

*I*−

*G*

^{h}

_{N}*G*

*given in Assumption DG2, we obtain:*

^{h}1

Δt θ*n1* 2

0,Ω*νκ θ**n1* 2

*X*≤0, 3.47

which yields that**θ***n1***0. The existence and uniqueness of the pressure***p*^{h}* _{n1}*is then obtained
from the inf-sup condition3.6.

Some approximation properties of the spaces **X*** ^{h}* and

*Q*

*are recalled next. From Crouzeix and Raviart 31, and Girault et al. 16, for each integer*

^{h}*r*≥ 1, and for any

**v**∈H

_{0}

^{1}Ω

^{2}, there is a unique discrete velocity

**v**∈

**X**

*such that*

^{h}*b*

**v**−**v, q**

0 ∀q∈*Q*^{h}*.* 3.48

Furthermore, if**v**∈H_{0}^{1}Ω^{2}∩H* ^{r1}*Ω

^{2}, there is a constant

*C*independent of

*h*such that

v−**v ** * _{X}* ≤

*Ch*

*|v|*

^{r}

_{r1,Ω}*,*3.49

|v−**v|*** _{m,Ω}* ≤

*Ch*

*|v|*

^{r1−m}

_{r1,Ω}*,*

*m*0,1. 3.50

For the pressure space, we use the approximation given by the*L*^{2}projection. For any
*q*∈*L*^{2}_{0}Ω, there exists a unique discrete pressure*q*∈*Q** ^{h}*such that

*q*−*q, z*

Ω 0 ∀z∈*Q*^{h}*.* 3.51
In addition, if*q*∈*H** ^{r}*Ω, then

*q*−*q*

*m,E*≤*Ch*^{r−m}*q*

*r,E**,* ∀E∈ E*h**, m*0,1,2. 3.52

The discrete Gronwall’s lemma plays an important role in the following analysis.

**Lemma 3.14**Discrete Gronwall’s Lemmacf. Heywood and Rannacher32. LetΔt,*H,and*
*a**n**,b**n**,c**n**,γ**n**(for integersn*≥0*be nonnegative numbers such that*

*a**l* Δt^{l}

*n0*

*b**n*≤Δt^{l}

*n0*

*γ**n**a**n* Δt^{l}

*n0*

*c**n**H* *forl*≥0. 3.53

*Suppose that*Δtγ*n**<1, for alln, and setσ**n* 1−Δtγ*n*^{−1}*. Then,*

*a**l* Δt^{l}

*n0*

*b**n*≤exp

Δt^{l}

*n0*

*σ**n**γ**n*

Δt^{l}

*n0*

*c**n**H*

*forl*≥0. 3.54

The family of time relaxation models is a regularization of NSE, and therefore it is natural to investigate the finite element error between the discretized model and the NSE.

To that end, we will assume that the solution to the Navier-Stokes equationsw, Pthat is approximated is a strong solution and in particular satisfiescf. Rivi`ere13

w*t**,***v**_{Ω}*νaw,***v ***J*_{0}w,**v ***c*** ^{w}**w; w,

**v**

*bv, P*f,

**v**

_{Ω}∀v∈

**X**

^{h}*,*3.55

*b*

**w, q**

0 ∀q∈*Q*^{h}*,* 3.56
**w**^{h}_{0}*,***v**

Ω u0*,***v**_{Ω} ∀v∈**X**^{h}*.* 3.57

**4. A Priori Error Estimates**

In this section, convergence of the scheme3.30–3.32is proved. Optimal error estimates in the energy norm are obtained.

* Theorem 4.1. Assume that w* ∈

*l*

^{2}0, T;H

*Ω*

^{r1}^{2}∩

*l*

^{2}0, T;H

^{2N2}Ω

^{2}

*,*

**w**

*t*∈

*l*

^{2}0, T; H

*Ω*

^{r1}^{2}∩L

^{∞}0, T×Ω,

**w**

*tt*∈

*L*

^{2}0, T;H

^{1}Ω

^{2},

*p*∈

*l*

^{2}0, T;

*H*

*Ω*

^{r}*and*

**u**0∈H

*Ω*

^{r1}^{2}

*.*

*Assume also that the coercivityLemma 3.1holds and that*3.24. If

**w satisfies***δis chosen of the order*

*ofh, and*Δt <

*1, there exists a constantC, independent ofhand*Δt

*but dependent onν*

^{−1}

*such that*

*the following error bound holds, for any 1*≤

*m*≤

*M:*

w*m*−**u**^{h}_{m}

0,Ω

*νκΔt*^{m}

*n1*

w*n*−**u**^{h}_{n}^{2}

*X*

_{1/2}

≤*Ch*^{r}

*ν*^{−1}*νχ*1

*CΔtCχh*^{2N2}*.*
4.1

*Remark 4.2. The dependence of these error estimates with respect to the Reynolds number Re*

∼1/νisORe*e*^{Re}, which is an improvement with respect to the continuous finite element
method where the dependence isORe*e*^{Re}^{3}.

*Proof. Defining***e***n***u*** ^{h}*t

*−*

^{n}**wt**

*and subtracting3.55from3.30, we have*

^{n}1

Δte*n1*−**e***n**,***v**_{Ω}−w*t*t*n1*,**v**_{Ω}*νaJ*0e*n1**,***v ***χ*

**e***n1*−*G*^{h}_{N}**e***n1**h**,***v**

Ω

*c*^{u}^{h}^{n}

**u**^{h}* _{n}*;

**u**

^{h}

_{n1}*,*

**v**

−*c*^{w}* ^{n1}*w

*n1*;

**w**

*n1*

*,*

**v**

*b*

**v, p**^{h}* _{n1}*−

*P*

*− 1*

_{n1}Δtw*n1*−**w***n**,***v**_{Ω}−*χ*

**w***n1*−*G*^{h}_{N}**w***n1**h**,***v**

Ω ∀v∈**X**^{h}*,*

4.2

*b*
**e***n1**, q*

0 ∀q∈*Q*^{h}*.* 4.3

We now decompose the error**e***n* *φ** _{n}*−

**η***, where*

_{n}*φ*

_{n}**u**

^{h}*−*

_{n}**w**

*n*and

**η***is the interpolation*

_{n}error**η**_{n}**w***n*−**w***n*. Choosing**v***φ** _{n1}*in the equation above, using the coercivity result3.5
and positivity of the operator

*I*−

*G*

^{h}

_{N}*G*

*, we obtain for4.2*

^{h}1 2Δt

*φ*_{n1}^{2}

0,Ω−*φ*_{n}^{2}

0,Ω

*νκφ*_{n1}^{2}

*X*

*c*^{u}^{h}^{n}

**u**^{h}* _{n}*;

**u**

^{h}

_{n1}*, φ*

_{n1}−*c*^{w}^{n1}

**w*** _{n1}*;

**w**

_{n1}*, φ*

_{n1}≤

**η*** _{t}*t

*n1*, φ

*n1*

Ω

*νaJ*_{0}

**η**_{n1}*, φ*_{n1}

**w***t*t*n1*− 1

Δt**w***n1*−**w***n*, φ*n1*

Ω

*χ*

**η*** _{n1}*−

*G*

^{h}

_{N}

**η**

_{n1}

^{h}*, φ*

_{n1}Ω−*χ*

**w***n1*−*G*^{h}_{N}**w***n1**h**, φ*_{n1}

Ω*b*

*φ*_{n1}*, P** _{n1}*−

*p*

^{h}

_{n1}*.*4.4

Consider now the nonlinear terms from the above equation. First note that since **w is**
continuous, we can rewrite

*c*^{w}^{n1}

**w*** _{n1}*;

**w**

_{n1}*, φ*

_{n1}*c*

^{u}

^{h}

^{n}**w*** _{n1}*;

**w**

_{n1}*, φ*

_{n1}*,* 4.5

so, for readability, the superscript **w**^{h}* _{n}* in the

*c*form is dropped. Therefore, adding and subtracting the interpolant

**w**

*n1*yields

*c*^{u}^{h}^{n}

**u**^{h}_{n}*,***u**^{h}_{n1}*, φ*_{n1}

−*c*^{u}^{h}^{n}

**w***n1**,***w***n1**, φ*_{n1}*c*

**u**^{h}_{n}*, φ*_{n1}*, φ*_{n1}

−*c*

*φ**n**, η*

_{n1}*, φ*

_{n1}*c*

*φ**n**,***w***n1**, φ*_{n1}

−*c*

**η**_{n}*,***w**_{n1}*, φ*_{n1}

−*c*

**w***n**, η*

_{n1}*, φ*

_{n1}−*c*

**w*** _{n1}*−

**w**

*n*

*,*

**w**

_{n1}*, φ*

_{n1}*.*

4.6

Thus, the error equation4.4is rewritten as 1

2Δt

*φ*_{n1}^{2}

0,Ω−*φ*_{n}^{2}

0,Ω

*νκφ*_{n1}^{2}

*X**c*

**u**^{h}_{n}*, φ*_{n1}*, φ*_{n1}

≤*c*

*φ*_{n}*, η*

_{n1}*, φ*

_{n1}*c*

*φ*_{n}*,***w**_{n1}*, φ*_{n1}*c*

**η**_{n}*,***w**_{n1}*, φ*_{n1}*c*

**w***n**, η*

_{n1}*, φ*

_{n1}*c*

**w***n1*−**w***n**,***w***n1**, φ*_{n1}**η*** _{t}*t

*n1*, φ

*n1*

Ω

**w***t*t*n1*− 1

Δt**w***n1*−**w***n*, φ*n1*

Ω

*νaJ*0

**η**_{n1}*, φ** _{n1}*
χ

**η*** _{n1}*−

*G*

^{h}

_{N}

**η**

_{n1}

^{h}*, φ*

_{n1}Ω

χ

**w*** _{n1}*−

*G*

^{h}

_{N}**w**

_{n1}

^{h}*, φ*

_{n1}Ω

b

*φ*_{n1}*, P** _{n1}*−

*p*

^{h}

_{n1}≤ |*T*0||T1|· · ·|T10|.

4.7

From property3.7, the term*cw*^{h}* _{n}*;

*φ*

_{n1}*, φ*

*in the left-hand side of4.7is positive and therefore it will be dropped. For the other terms of the form*

_{n1}*c·,*·,·that appear on the right- hand side of the above error equation we obtain bounds, exactly as in the proof of Theorem 5.2

in19by Kaya and Rivi`ere. The constant*C*is a generic constant that is independent of*h, ν*
andΔt, and that takes diﬀerent values at diﬀerent places

|T0|*c*

*φ**n**, η*

_{n1}*, φ*

*≤*

_{n1}*νκ*

20*φ*_{n1}^{2}

*X**C*
*νφ**n*^{2}

0,Ω*,*

|T1|*c*

*φ*_{n}*,***w**_{n1}*, φ** _{n1}*≤

*νκ*

20*φ*_{n1}^{2}

*X**C*
*νφ*_{n}^{2}

0,Ω*,*

|T2|*c*

*η*_{n}*,***w***n1**, φ** _{n1}*≤

*νκ*

20*φ*_{n1}^{2}

*X**C*

*νh*^{2r}|w*n*|^{2}_{2r,Ω}*,*

|T3|*c*

**w***n**, η*

_{n1}*, φ*

*≤*

_{n1}*νκ*

20*φ*_{n1}^{2}

*X* *C*

*νh*^{2r}|w*n*|^{2}_{2r,Ω}*,*

|T4|*c*

**w***n1*−**w***n**,***w***n1**, φ** _{n1}*≤

*νκ*

20*φ*_{n1}^{2}

*X**C*

*ν*Δt^{2} w*t* 2

*L*^{∞}t*n**,t**n1*×Ω*.*

4.8

Therefore, we have

|T0|· · ·|T4| ≤ 5νκ

20 *φ*_{n1}^{2}

*X**Cν*^{−1}*φ**n*^{2}_{0,Ω}*Cν*^{−1}*h*^{2r}|w*n*|^{2}_{r1,Ω}*Cν*^{−1}Δt^{2} w*t* 2

*L*^{∞}t*n**,t** _{n1}*×Ω

*.*4.9

To bound*T*5, Cauchy-Schwarz’s inequality, Young’s inequality and the approximation
result3.50for**w***t*are applied

|T5| ≤*φ*_{n1}

0,Ω**η*** _{t}*t

_{n1}0,Ω

≤ 1

8*φ*_{n1}^{2}

0,Ω*Ch*^{2r2}w*t*

*t*^{n1}^{2}

*k1,Ω**.*

4.10

To bound the term*T*_{6}, a Taylor expansion with integral remainder is used

**w***n***w***n1*−Δt**w***t*t*n1* 1
2

_{t}_{n1}

*t**n*

s−*t**n***w***tt*sds. 4.11

This implies that

**w***t*t* _{n1}*−

**w**

*−*

_{n1}**w**

*n*

Δt

^{2}_{0,Ω} ≤ Δt
6

_{t}_{n1}

*t**n*

**w***tt*s ^{2}_{0,Ω}*ds.* 4.12

Thus, with3.50, we have

|T6| ≤ 1

8*φ*_{n1}^{2}

0,Ω*CΔt*
_{t}_{n1}

*t**n*

**w***tt*s ^{2}_{0,Ω}*ds*

≤ 1

8*φ*_{n1}^{2}_{0,Ω}*CΔt*
_{t}_{n1}

*t**n*

w*tt*s ^{2}_{0,Ω}*ds.*

4.13

Next, the term*T*7is expanded as

|T7| ≤
*ν*

*E∈E**h*

*E*

∇η* _{n1}*:∇φ

_{n1}*ν*

*e∈Γ**h*∪∂Ω

*e*

∇η_{n1}**n***e*·

*φ*_{n1}

*ν *_{a}

*e∈Γ**h*∪∂Ω

*e*

∇φ_{n1}**n***e*·

**η**_{n1}*νJ*_{0}

**η**_{n1}*, φ** _{n1}*
|T71||T72||T73||T74|.

4.14

The term *T*_{71} is bounded using Cauchy-Schwarz inequality, Young’s inequality and the
approximation result3.49

|T71| ≤*νφ*_{n1}

*X***η**_{n1}

*X*

≤ *νκ*

20*φ*_{n1}^{2}

*X**Cν η*

_{n1}^{2}

*X*

≤ *νκ*

20*φ*_{n1}^{2}

*X**Cνh*^{2r}|w*n1*|^{2}_{r1,Ω}*.*

4.15

Using Cauchy-Schwarz’s inequality, trace inequality2.8 and approximation result3.50, we have

|T72| ≤*ν*

*e∈Γ**h*∪∂Ω

{∇η* _{n1}*}n

*e*

0,e

*e∈Γ**h*∪∂Ω

φ_{n1}

0,e

≤*Cν*

*e∈Γ**h*∪∂Ω

1

|e|*φ*_{n1}^{2}

0,e

_{1/2}∇η_{n1}

0,Ω*h*∇^{2}**η**_{n1}

0,Ω

≤ *νκ*

20*φ*_{n1}^{2}

*X**Cνh*^{2r}|w*n1*|^{2}_{r1,Ω}*.*

4.16

Using Cauchy-Schwarz’s inequality, trace inequality2.10, and approximation result3.49, we have

|T73| ≤*ν*

*e∈Γ**h*∪∂Ω

∇φ_{n1}**n***e*^{2}

0,e

_{1/2}

*e∈Γ**h*∪∂Ω

**η**_{n1}^{2}

0,e

_{1/2}

≤*Cνφ*_{n1}

*X*

*e∈Γ**h*∪∂Ω

1

|e|**η**_{n1}^{2}

0,e

1/2

≤ *νκ*

20*φ*_{n1}^{2}

*X**Cνh*^{2r}|w*n1*|^{2}_{r1,Ω}*.*

4.17