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

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 isORee^Re, which is an improvement with respect to the continuous finite element method where the dependence isORee^Re3.

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/4in 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 differential filtercf. Germano8

φGφ in Ω, 1.1

φ0 on ∂Ω, 1.2

whereG⁻¹ −δ²Δ 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 eachN0,1, . . . ,it computes an approximate solutionuNby N steps of a fixed point iteration for the fixed point problemcf. Bertero and Boccacci10

givenu solve uu{u−Gu}, foru. 1.3

The deconvolution approximation is then computed as follows.

Algorithm 1.1van Cittert approximate deconvolution algorithm. Consider thatu0 u, for n1,2, . . . , N−1, perform

u_n1un{u−Gun}. 1.4

By eliminating the intermediate steps, it is easy to find an explicit formula for theNth deconvolution operatorG_Ngiven by

G_Nu :^N

n0

I−Gⁿu. 1.5

For example, the approximate deconvolution operator corresponding toN0,1,2 is G0uu,

G₁u2u−u, G₂u3u−3uu.

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_Nu, 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

utu· ∇u−νΔu∇pχu−G_Nu f in0, T×Ω, 1.7

∇ ·u0 in 0, T×Ω, 1.8

u0 in 0, T×∂Ω, 1.9

u0,· u0· inΩ, 1.10

whereΩ⊂R², is a convex bounded regular domain with boundary∂Ω,u is the fluid velocity, pis the fluid pressure andf 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 u0. A pressure normalization condition

Ωp0 is also needed for uniqueness of the pressure.

The time relaxation coefficientχ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 degreer andr −1, respectively. Because of the lack of continuity constraint between elements, the Discontinuous Galerkin DG methods offer 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 ashand δ → 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 triangulationsEhofΩ, consisting of triangles of maximum diameterh, is introduced. Leth_Edenote the diameter of a triangleE andρ_Ethe diameter of its inscribed circle. Regulary, it is meant that there exists a parameter ζ >0, independent ofh, such that

hE

ρ_E ζ_E≤ζ, ∀E∈ Eh. 2.1

This assumption will be used throughout this work.Γhdenotes the set of all interior edges of Eh. Letedenote a segment ofΓhshared by two trianglesE^kandE^lk < lofEh; it is associated withea specific unit normal vectorne directed fromE^ktoE^land the jump and average of a functionφoneis formally defined by

φ

φ

E^k

e− φ

E^l

e, φ

1 2

φ

E^k

e1 2

φ

E^l

e. 2.2

Ifebelongs to the boundary∂Ω, thenneis the unit normaln exterior toΩand the jump and the average ofφonecoincide with the trace ofφone.

Here, for any domainO,L²Ois the classical space of square-integrable functions with inner-productf, g_O

Ofgand norm · _0,O. The spaceL²₀Ωis the subspace of functions ofL²Ωwith zero mean value

L²₀Ω

v∈L²Ω:

Ωv0

. 2.3

The standard Sobolev spaces are denoted byW_p^rΩwhereW₂^rΩis theH^rΩ, with norm · _r,Ωand seminorm| · |_r,Ω.

Next, the discrete velocity and pressure spaces are defined to be consisting of discontinuous piecewise polynomials. For any positive integerr, the corresponding finite- dimensional spaces are

X^h

v∈

L²Ω2

:∀E∈ Eh, v∈PrE²

, Q^h

q∈L²₀Ω:∀E∈ Eh, q∈Pr−1E ,

2.4

wherePrE span{xⁱy^j :ij ≤ r}is defined as the span of polynomials of orderr over triangleE.

Denoting by|e|the measure ofe, the following norms are associated for the spacesX^h andQ^h

v _X

| ∇v |²_0,Ω

e∈Γh∪∂Ω

1

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

, q

Q q

0,Ω,

2.5

where| v |_0,Ωis the broken norm defined by

| v |_0,Ω

E∈Eh

v ²_0,E 1/2

. 2.6

Finally, some trace and inverse inequalities are recalled, that hold true on each element EinEh, with diameterh_E, the constantCis independent ofh_E

v _0,e≤C

h^−1/2_E v _0,Eh^1/2_E ∇v _0,E

, ∀e∈∂E, ∀v∈

H¹E₂

, 2.7

∇v _0,e≤C

h^−1/2_E ∇v _0,Eh^1/2_E ∇²v

0,E

, ∀e∈∂E, ∀v∈

H²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_E ∇v _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 defineda:X^h×X^h → R, andJ0:X^h×X^h → Rby

az,v

E∈Eh

E

∇z :∇v−

e∈Γh∪∂Ω

e

{∇z}ne·v a

e∈Γh∪∂Ω

e

{∇v}ne·z,

J0z,v

e∈Γh∪∂Ω

σ

|e|

e

z·v.

3.1

The parameter _a takes the value−1,0 or 1: this will yield different 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 · _X. In the case where a −1, the bilinear formais symmetric; otherwise it is nonsymmetric. We remark that the formau,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 sufficiently large to obtain coercivity of asee Lemma 3.1. If a1, then the jump parameterσis taken equal to 1.

The incompressibility condition1.8 is enforced by means of the bilinear formb : X^h×Q^h → Rdefined by

b v, q

−

E∈Eh

E

q∇ ·v

e∈Γh∪∂Ω

e

q

v·ne. 3.2

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

c^zu; v,t

E∈Eh

E

u· ∇v·t1 2

E

∇ ·uv·t

−1 2

e∈Γh∪∂Ω

e

u·ne{v·t}

E∈Eh

∂E−

|{u} ·nE|

v^int−v^ext

·t^int,

3.3

where

∂E₋{x∈∂E:{z} ·nE<0}, 3.4

the superscriptz denotes the dependence of∂E₋onz and the superscript intresp., extrefers to the trace of the function on a side ofEcoming from the interior ofEresp., coming from the exterior ofEon that side. When the side ofEbelongs 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 toz, but linear with respect to u,v andt.

Some important properties satisfied by the formsa, b, ccf. Wheeler20, and Girault et al.16,17are recalled.

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

av,v J0v,v≥κ v ²_X, ∀v∈X^h. 3.5

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

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

q∈Qinf^hsup

v∈X^h

b v, q v _Xq

0,Ω

≥β. 3.6

Lemma 3.3Positivity. One has

c^vv,z,z≥0, ∀v,z∈

t∈

L²Ω₂

: t|_E∈

H²E₂

∀E∈ Eh

. 3.7

The discrete form of the differential filter1.1is defined following the work of Manica and Merdan22.

Definition 3.4Discrete differential filter. Givenv ∈ L²Ω, for a given filtering radiusδ >

0, G^h:L²Ω → X^h,v^hG^hv where v^his the unique solution inX^hof

δ² a

v^h, φ J₀

v^h, φ

v^h, φ

v, φ

∀φ∈X^h. 3.8

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

Definition 3.6. The discrete van Cittert deconvolution operatorsG^h_Nare

G^h_Nv :^N

n0

Πh−Ghⁿv, 3.9

whereΠh:L²Ω → X^his theL²projection.

Forv∈X^h, the discrete deconvolution operator forN0,1,2 is G^h₀vv,

G^h₁v2v−v^h, G^h₂v3v−3v^hv^hh.

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. GNis a bounded, self-adjoint positive operator.GNis anOδ^2N2asymptotic inverse to the filterG. Specifically, for smoothφand asδ → 0,

φGNφ −1^N1δ^2N2Δ^N1G^N1φ. 3.11 Some basic facts about discrete differential filters and deconvolution operators are presented next.

Lemma 3.8. For v∈L²Ω, one has the following bounds for the discretely filtered and approximately deconvolvedv:

v^h

0,Ω≤ v _0,Ω, 3.12

G^h_Nv^h

0,Ω≤CN v _0,Ω. 3.13

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

δ²κv^h2

Xv^h2

0,Ω≤ v,v^h

Ω≤ 1

2 v ²_0,Ω1 2

v^h2

0,Ω. 3.14

The termδ²κ v^{h 2}_Xis positive, so it will be dropped, which yields 1

2 v^h2

0,Ω ≤ 1

2 v ²_0,Ω. 3.15

Multiplying by 2 and taking the square root yields the estimate3.12. Equation3.13follows immediately from3.12and the definition ofG^h_N.

Lemma 3.9. For smoothφthe discrete approximate deconvolution operator satisfies φ−G^h_Nφ^h

0,Ω ≤C₁δ^2N2Gφ

2N2,ΩC₂

δh^rh^r1N1

n2

Gⁿφ

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^rh^r1φ

r1,Ω. 3.19

Lemma 3.8gives for the third term in3.17that G^h_Nφ−G^h_Nφ^h ≤ C φ−φ^h . Then, 3.19is applied.

Now, it is left to bound the second term from 3.17. First, note that for N 0, G0φ^h−G^h₀φ^h _0,Ω 0. Based on Definition 3.6 of continuous and discrete deconvolution operators and their expansionsee3.10,G_N is a polynomial of degree Nin GandG^h_N inG^has well. Thus, the second term in3.17can be written as

GNφ−G^h_Nφ

0,Ω

N n0

αn

Gⁿφ− G^h_n

φ 0,Ω

≤^N

n0

αnGⁿφ− G^h_n

φ

0,Ω. 3.20

ForO1coefficientsα_nand forN1, the result3.19gives Gφ−

G^h φ

0,Ω

φ−φ

h

0,Ω

≤C

δh^rh^r1 φ

r1,Ω

.

3.21

ForN2, the results3.19and3.12give G2φ−

G^h₂ φ

0,Ω φ−φ^h

hh 0,Ω

≤ φ−φ

h

0,Ω

φ

h

−φ^h

hh 0,Ω

≤ φ−φ

h

0,Ω

φ−φ^h

0,Ω

≤C

δh^rh^r1 φ

r1,Ω

φ

_r1,Ω

.

3.22

Inductively,

GNφ−G^h_Nφ

0,Ω≤C

δh^r h^r1N1

n2

Gⁿφ

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ⁿφ|_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 allk

Gφ_r1,Ω≤Cφ_r1,Ω and thusGⁿφ_r1,Ω≤Cφ_r1,Ω. 3.24 The nonperiodic case can be more delicate. Suppose∂Ω ∈ C^r3 andφ 0 on∂Ω i.e.,φ ∈ H₀¹Ω∩H^r1Ω. CallGφφsoφsatisfies

−δ²Δφφφ inΩ, φ0 on∂Ω. 3.25

Then it is known thatφ∈H^r3Ω∩H₀¹Ω, andΔφ0 on∂Ω. Further, φ

j,Ω≤Cφ_j,Ω, j0,1,2. 3.26

So,3.24holds forr −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 termn 2, that is,G²φ φ. We know from elliptic theory forφ∈H^r1Ω

H₀¹Ω, thatφ∈H^r3Ω

H₀¹Ω,as noted above Δφ0 on∂Ωand

−δ²Δφφφ inΩ, φ Δφ0 on ∂Ω. 3.27

Theorem 1.1 in28then implies, uniformly inδ, φ

j,Ω ≤Cφ

j,Ω, j0,1,2,3,4. 3.28

This extends directly toGⁿφ.

ExtendingLemma 3.9, the following assumption will be made.

Assumption DG1. The|Gⁿφ|_r1,Ωterms in3.16are independent ofδand φ−G^h_Nφ^h

≤C1δ^2N2φ

2N2,ΩC2

δh^rh^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^handI−G^h_NG^hare symmetric, positive definiteSP Doperators.

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 forma·,·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 thatM T/Δtis a positive integer. Let 0 t₀ < t₁ < · · · < t_M T be a subdivision of the interval 0, T. The functionφevaluated at the timet_mis denoted byφ_m. With the above forms, the fully discrete scheme is: findu^h_n, p^h_nn≥0∈X^h×Q^hsuch that:

1 Δt

u^h_n1−u^h_n,v

Ων a

u^h_n1,v J₀

u^h_n1,v χ

u^h_n1−G^h_Nu^h_n1^h,v

Ω

c^uhn

u^h_n;u^h_n1,v b

v, p_n1^h

fn1,v_Ω ∀v∈X^h,

3.30

b

u^h_n1, q

0 ∀q∈Q^h, 3.31 u^h₀,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.12Consistency. Letu, pbe the solution to1.7–1.10, thenu, psatisfies

ut,v_Ωνau,v J0u,v χ

u−G^h_Nu^h,v

Ωc^uu; 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_Nu^h,v_Ω−χu−GNu,v_Ω. Furthermore, Eu,v≤C

δh^r h^r1

|u|_r1,Ω v _0,Ω. 3.36

Proof. Equations3.34and3.35are clearly satisfied because of1.8,1.9, and1.10and the regularity ofu. Next, we multiply1.7byv and integrate over one mesh elementE

ut,v_E ∇ ·uu,v_E−νΔu,v_Eχu−G_Nu,v_E

∇p,v

E f,v_E. 3.37

Summing over all elementsE

E∈Eh

ut,v_E ∇ ·uu,v_Ω−ν

E∈Eh

Δu,v_Eχu−G_Nu,v_Ω

∇p,v

Ω f,v_Ω. 3.38

By Green’s formula

−ν

E∈Eh

Δu,v_E−νΔu,v_Ων

E∈Eh

∇u,∇v_E−ν

e∈Γh∪∂Ω

e

∇une·v. 3.39

The regularity ofu then yields

−ν

E∈Eh

Δu,v_Eνau,v J₀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 ofu yields

u· ∇u,v_Ω c^uu; u,v. 3.42

The final result is obtained by bounding the consistency errorEu,v. ForN0, we have Eu,v≤u−u^h

0,Ω v _0,Ω ≤C

δh^rh^r1

|u|_r1,Ω v _0,Ω. 3.43

ForN1,

Eu,v≤Gu−G^hu^h

0,Ω v _0,Ω

≤

u−u^hh

0,Ω v _0,Ω

≤

u−u^h 0,Ω

u^h−u^hh 0,Ω

v _0,Ω

≤

u−u^h

0,Ωu−u^h

0,Ω

v _0,Ω

≤C

δh^rh^r1u

r1,Ω|u|_r1,Ω

v _0,Ω.

3.44

The bound forEu,vis obtained by applying an induction argument andRemark 3.10.

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

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

3.32.

Proof. The existence ofu^h₀ is trivial. Givenu^h_n, the problem of finding a uniqueu^h_n1satisfying 3.30-3.31is linear and finite-dimensional. Therefore, it suffices to show uniqueness of the solution. Consider the problem restricted to the subspaceV^hdefined by

V^h

v∈X^h:b v, q

0∀q∈Q^h

. 3.45

Letu^h_n1andu^h_n1be two solutions and letθn1u^h_n1−u^h_n1. Then,θn1satisfies:

1

Δtθ_n1,v_Ωνaθ_n1,v J₀θ_n1,v c^uhn

u^h_n;θ_n1,v χ

θn1−G^h_Nθn1h

,v

0 ∀v∈V^h.

3.46

Choosingv θ_n1 and using the coercivity result3.5, positivity result3.7and positivity of the operatorI−G^h_NG^hgiven in Assumption DG2, we obtain:

1

Δt θn1 2

0,Ωνκ θn1 2

X≤0, 3.47

which yields thatθn10. The existence and uniqueness of the pressurep^h_n1is then obtained from the inf-sup condition3.6.

Some approximation properties of the spaces X^h and Q^h are recalled next. From Crouzeix and Raviart 31, and Girault et al. 16, for each integer r ≥ 1, and for any v∈H₀¹Ω², there is a unique discrete velocityv∈X^hsuch that

b

v−v, q

0 ∀q∈Q^h. 3.48

Furthermore, ifv∈H₀¹Ω²∩H^r1Ω², there is a constantCindependent ofhsuch that

v−v _X ≤Ch^r|v|_r1,Ω, 3.49

|v−v|_m,Ω ≤Ch^r1−m|v|_r1,Ω, m0,1. 3.50

For the pressure space, we use the approximation given by theL²projection. For any q∈L²₀Ω, there exists a unique discrete pressureq∈Q^hsuch that

q−q, z

Ω 0 ∀z∈Q^h. 3.51 In addition, ifq∈H^rΩ, then

q−q

m,E≤Ch^r−mq

r,E, ∀E∈ Eh, m0,1,2. 3.52

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

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

al Δt^l

n0

bn≤Δt^l

n0

γnan Δt^l

n0

cnH forl≥0. 3.53

Suppose thatΔtγn<1, for alln, and setσn 1−Δtγn⁻¹. Then,

al Δt^l

n0

bn≤exp

Δt^l

n0

σnγn

Δt^l

n0

cnH

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

wt,v_Ωνaw,v J₀w,v c^ww; w,v bv, P f,v_Ω ∀v∈X^h, 3.55 b

w, q

0 ∀q∈Q^h, 3.56 w^h₀,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²0, T;H^r1Ω² ∩ l²0, T;H^2N2Ω², wt ∈ l²0, T; H^r1Ω²∩L^∞0, T×Ω,wtt∈L²0, T;H¹Ω²,p∈l²0, T;H^rΩandu0∈H^r1Ω². Assume also that the coercivityLemma 3.1holds and thatw satisfies3.24. Ifδis chosen of the order ofh, andΔt <1, there exists a constantC, independent ofhandΔtbut dependent onν⁻¹such that the following error bound holds, for any 1≤m≤M:

wm−u^h_m

0,Ω

νκΔt^m

n1

wn−u^h_n²

X

_1/2

≤Ch^r

ν⁻¹νχ1

CΔtCχh^2N2. 4.1

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

∼1/νisORee^Re, which is an improvement with respect to the continuous finite element method where the dependence isORee^Re3.

Proof. Definingenu^htⁿ−wtⁿand subtracting3.55from3.30, we have

1

Δten1−en,v_Ω−wttn1,v_ΩνaJ0en1,v χ

en1−G^h_Nen1h,v

Ω

c^uhn

u^h_n;u^h_n1,v

−c^wn1wn1;wn1,v b

v, p^h_n1−P_n1 − 1

Δtwn1−wn,v_Ω−χ

wn1−G^h_Nwn1h,v

Ω ∀v∈X^h,

4.2

b en1, q

0 ∀q∈Q^h. 4.3

We now decompose the erroren φ_n−η_n, whereφ_n u^h_n−wnandη_nis the interpolation

errorη_nwn−wn. Choosingvφ_n1in the equation above, using the coercivity result3.5 and positivity of the operatorI−G^h_NG^h, we obtain for4.2

1 2Δt

φ_n1²

0,Ω−φ_n²

0,Ω

νκφ_n1²

X

c^uhn

u^h_n;u^h_n1, φ_n1

−c^wn1

w_n1;w_n1, φ_n1

≤

η_ttn1, φn1

Ω

νaJ₀

η_n1, φ_n1

wttn1− 1

Δtwn1−wn, φn1

Ω

χ

η_n1−G^h_Nη_n1^h, φ_n1

Ω−χ

wn1−G^h_Nwn1h, φ_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^wn1

w_n1;w_n1, φ_n1 c^uhn

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 interpolantwn1yields

c^uhn

u^h_n,u^h_n1, φ_n1

−c^uhn

wn1,wn1, φ_n1 c

u^h_n, φ_n1, φ_n1

−c

φn,η_n1, φ_n1 c

φn,wn1, φ_n1

−c

η_n,w_n1, φ_n1

−c

wn,η_n1, φ_n1

−c

w_n1−wn,w_n1, φ_n1 .

4.6

Thus, the error equation4.4is rewritten as 1

2Δt

φ_n1²

0,Ω−φ_n²

0,Ω

νκφ_n1²

Xc

u^h_n, φ_n1, φ_n1

≤c

φ_n,η_n1, φ_n1c

φ_n,w_n1, φ_n1c

η_n,w_n1, φ_n1 c

wn,η_n1, φ_n1c

wn1−wn,wn1, φ_n1η_ttn1, φn1

Ω

wttn1− 1

Δtwn1−wn, φn1

Ω

νaJ0

η_n1, φ_n1 χ

η_n1−G^h_Nη_n1^h, φ_n1

Ω

χ

w_n1−G^h_Nw_n1^h, φ_n1

Ω

b

φ_n1, P_n1−p^h_n1

≤ |T0||T1|· · ·|T10|.

4.7

From property3.7, the termcw^h_n;φ_n1, φ_n1in the left-hand side of4.7is positive and therefore it will be dropped. For the other terms of the formc·,·,·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 constantCis a generic constant that is independent ofh, ν andΔt, and that takes different values at different places

|T0|c

φn,η_n1, φ_n1≤ νκ

20φ_n1²

XC νφn²

0,Ω,

|T1|c

φ_n,w_n1, φ_n1≤ νκ

20φ_n1²

XC νφ_n²

0,Ω,

|T2|c

η_n,wn1, φ_n1≤ νκ

20φ_n1²

XC

νh^2r|wn|²_2r,Ω,

|T3|c

wn,η_n1, φ_n1≤ νκ

20φ_n1²

X C

νh^2r|wn|²_2r,Ω,

|T4|c

wn1−wn,wn1, φ_n1≤ νκ

20φ_n1²

XC

νΔt² wt 2

L^∞tn,tn1×Ω.

4.8

Therefore, we have

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

20 φ_n1²

XCν⁻¹φn²_0,ΩCν⁻¹h^2r|wn|²_r1,ΩCν⁻¹Δt² wt 2

L^∞tn,t_n1×Ω. 4.9

To boundT5, Cauchy-Schwarz’s inequality, Young’s inequality and the approximation result3.50forwtare applied

|T5| ≤φ_n1

0,Ωη_tt_n1

0,Ω

≤ 1

8φ_n1²

0,ΩCh^2r2wt

tⁿ¹²

k1,Ω.

4.10

To bound the termT₆, a Taylor expansion with integral remainder is used

wnwn1−Δtwttn1 1 2

_tn1

tn

s−tnwttsds. 4.11

This implies that

wtt_n1−w_n1−wn

Δt

²_0,Ω ≤ Δt 6

_tn1

tn

wtts ²_0,Ωds. 4.12

Thus, with3.50, we have

|T6| ≤ 1

8φ_n1²

0,ΩCΔt _tn1

tn

wtts ²_0,Ωds

≤ 1

8φ_n1²_0,ΩCΔt _tn1

tn

wtts ²_0,Ωds.

4.13

Next, the termT7is expanded as

|T7| ≤ ν

E∈Eh

E

∇η_n1:∇φ_n1

ν

e∈Γh∪∂Ω

e

∇η_n1 ne·

φ_n1

ν _a

e∈Γh∪∂Ω

e

∇φ_n1 ne·

η_n1 νJ₀

η_n1, φ_n1 |T71||T72||T73||T74|.

4.14

The term T₇₁ is bounded using Cauchy-Schwarz inequality, Young’s inequality and the approximation result3.49

|T71| ≤νφ_n1

Xη_n1

X

≤ νκ

20φ_n1²

XCνη_n1²

X

≤ νκ

20φ_n1²

XCνh^2r|wn1|²_r1,Ω.

4.15

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

|T72| ≤ν

e∈Γh∪∂Ω

{∇η_n1}ne

0,e

e∈Γh∪∂Ω

φ_n1

0,e

≤Cν

e∈Γh∪∂Ω

1

|e|φ_n1²

0,e

_1/2∇η_n1

0,Ωh∇²η_n1

0,Ω

≤ νκ

20φ_n1²

XCνh^2r|wn1|²_r1,Ω.

4.16

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

|T73| ≤ν

e∈Γh∪∂Ω

∇φ_n1 ne²

0,e

_1/2

e∈Γh∪∂Ω

η_n1²

0,e

_1/2

≤Cνφ_n1

X

e∈Γh∪∂Ω

1

|e|η_n1²

0,e

1/2

≤ νκ

20φ_n1²

XCνh^2r|wn1|²_r1,Ω.

4.17