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

Authors: Monika Neda
Journal: Advances in Numerical Analysis
 (2010)

---

## Abstract

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 is𝒪(ReeRe), which is an improvement with respect to the continuous finite element method where the dependence is 𝒪(ReeRe3).


---

## Body

## 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 theory [1], 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 be 𝒪(Re9/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 equations (cf. Adams and Stolz [2] and Layton and Neda [3]). 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 possible (cf. Berselli et al. [4], John [5], Geurts [6], Sagaut [7] and Germano [8]). First, consider the equations of differential filter (cf. Germano [8])(1.1)ϕ¯=GϕinΩ,(1.2)ϕ¯=0on∂Ω,
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 Simulation (LES) pioneered by Stolz and Adams (cf. Stolz and Adams [2, 9]). For each N=0,1,…, it computes an approximate solution uN by N steps of a fixed point iteration for the fixed point problem (cf. Bertero and Boccacci [10]) (1.3)givenu¯solveu=u+{u¯-Gu},foru.
The deconvolution approximation is then computed as follows.Algorithm 1.1 (van Cittert approximate deconvolution algorithm).
Consider thatu0=u¯, for n=1,2,…,N-1, perform
(1.4)un+1=un+{u¯-Gun}.By eliminating the intermediate steps, it is easy to find an explicit formula for theNth deconvolution operator GN given by(1.5)GNu¯:=∑n=0N(I-G)nu¯.
For example, the approximate deconvolution operator corresponding to N=0,1,2 is (1.6)G0u¯=u¯,G1u¯=2u¯-u¯¯,G2u¯=3u¯-3u¯¯+u¯¯¯.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, 12] and arises by adding a linear, lower order time regularization term, χu′=χ(u-GNu¯), 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(1.7)ut+u·∇u-νΔu+∇p+χ(u-GNu¯)=fin(0,T]×Ω,(1.8)∇·u=0in[0,T]×Ω,(1.9)u=0in[0,T]×∂Ω,(1.10)u(0,·)=u0(·)inΩ,
where Ω⊂ℝ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 u0. A pressure normalization condition ∫Ωp=0 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 in [13] for the case of Stokes and Navier-Stokes problems.Existence, uniqueness and regularity of strong solutions of these models are discussed in [3]. 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. In [14], 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 in [3, 15], respectively.In this work, a class of discontinuous finite element methods for solving high-order time relaxation family of fluid models (1.7)–(1.10) is 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 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; (iii) unstructured meshes and domains with complicated geometries are easily handled. In the case of DNS, DG methods have been applied to the steady-state NSE (cf. Girault et al. [16]) and 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 model) is derived in [19] by Kaya and Rivière.This paper is organized in the following way. Section2 introduces some notation and mathematical properties. In Section 3, the fully discrete schemes are introduced and it is proved that the schemes solutions are computable. A priori velocity error estimates are derived in Section 4. The family of models (1.7)–(1.10) is regularization of the NSE. Thus, the correct question is to study convergence of discretizations of (1.7)–(1.10) to solutions of the NSE as h and δ→0 (rather 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 triangulationsℰh of Ω¯, consisting of triangles of maximum diameter h, is introduced. Let hE denote the diameter of a triangle E and ρE the diameter of its inscribed circle. Regulary, it is meant that there exists a parameter ζ>0, independent of h, such that
(2.1)hEρE=ζE≤ζ,∀E∈ℰh.
This assumption will be used throughout this work. Γh denotes the set of all interior edges of ℰh. Let e denote a segment of Γh shared by two triangles Ek and El (k<l) of ℰh; it is associated with e a specific unit normal vector ne directed from Ek to El and the jump and average of a function ϕ on e is formally defined by(2.2)[ϕ]=(ϕ|Ek)|e-(ϕ|El)|e,{ϕ}=12(ϕ|Ek)|e+12(ϕ|El)|e.
If e belongs to the boundary ∂Ω, then ne 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 domain𝒪, L2(𝒪) is the classical space of square-integrable functions with inner-product (f,g)𝒪=∫𝒪fg and norm ∥·∥0,𝒪. The space L02(Ω) is the subspace of functions of L2(Ω) with zero mean value
(2.3)L02(Ω)={v∈L2(Ω):∫Ωv=0}.
The standard Sobolev spaces are denoted by Wpr(Ω) where W2r(Ω) is the Hr(Ω), 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
(2.4)Xh={v∈(L2(Ω))2:∀E∈ℰh,v∈(ℙr(E))2},Qh={q∈L02(Ω):∀E∈ℰh,q∈ℙr-1(E)},
where Pr(E)=span{xiyj:i+j≤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 Xh and Qh(2.5)∥v∥X=(|∥∇v∥|0,Ω2+∑e∈Γh∪∂Ω1|e|∥[v]∥0,e2)1/2,∥q∥Q=∥q∥0,Ω,
where |∥v∥|0,Ω is the broken norm defined by
(2.6)|∥v∥|0,Ω=(∑E∈ℰh∥v∥0,E2)1/2.Finally, some trace and inverse inequalities are recalled, that hold true on each elementE in ℰh, with diameter hE, the constant C is independent of hE(2.7)∥v∥0,e≤C(hE-1/2∥v∥0,E+hE1/2∥∇v∥0,E),∀e∈∂E,∀v∈(H1(E))2,(2.8)∥∇v∥0,e≤C(hE-1/2∥∇v∥0,E+hE1/2∥∇2v∥0,E),∀e∈∂E,∀v∈(H2(E))2,(2.9)∥v∥0,e≤ChE-1/2∥v∥0,E,∀e∈∂E,∀v∈Xh,(2.10)∥∇v∥0,e≤ChE-1/2∥∇v∥0,E,∀e∈∂E,∀v∈Xh.


## 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:Xh×Xh→ℝ, and J0:Xh×Xh→ℝ by
(3.1)a(z,v)=∑E∈ℰh∫E∇z:∇v-∑e∈Γh∪∂Ω∫e{∇z}ne·[v]+ϵa∑e∈Γh∪∂Ω∫e{∇v}ne·[z],J0(z,v)=∑e∈Γh∪∂Ωσ|e|∫e[z]·[v].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 form a is symmetric; otherwise it is nonsymmetric. We remark that the form a(u,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 a (see Lemma 3.1). If ϵa=1, then the jump parameter σ is taken equal to 1.The incompressibility condition (1.8) is enforced by means of the bilinear form b:Xh×Qh→ℝ defined by (3.2)b(v,q)=-∑E∈ℰh∫Eq∇·v+∑e∈Γh∪∂Ω∫e{q}[v]·ne.
Finally, the DG discretization of the nonlinear convection term w·∇w, which was introduced in [16] by Girault et al. and studied extensively in [16, 17] by the same authors, is recalled as follows: (3.3)cz(u;v,t)=∑E∈ℰh(∫E(u·∇v)·t+12∫E(∇·u)v·t)-12∑e∈Γh∪∂Ω∫e[u]·ne{v·t}+∑E∈ℰh∫∂E-|{u}·nE|(vint-vext)·tint,
where (3.4)∂E-={x∈∂E:{z}·nE<0},
the superscript z denotes the dependence of ∂E- on z and the superscript int (resp., ext) refers 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 formsa,b,c (cf. Wheeler [20], and Girault et al. [16, 17]) are recalled.Lemma 3.1 (Coercivity).
Ifϵa=1, assume that σ=1. If ϵa∈{-1,0}, assume that σ is sufficiently large. Then, there is a constant κ>0, independent of h, such that
(3.5)a(v,v)+J0(v,v)≥κ∥v∥X2,∀v∈Xh.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 in [21] by Epshteyn and Rivière.Lemma 3.2 (Inf-sup condition).
There exists a positive constantβ, independent of h such that
(3.6)infq∈Qhsupv∈Xhb(v,q)∥v∥X∥q∥0,Ω≥β.Lemma 3.3 (Positivity).
One has(3.7)cv(v,z,z)≥0,∀v,z∈{t∈(L2(Ω))2:t|E∈(H2(E))2∀E∈ℰh}.The discrete form of the differential filter (1.1) is defined following the work of Manica and Merdan [22].Definition 3.4 (Discrete differential filter).
Givenv∈L2(Ω), for a given filtering radius δ>0,Gh:L2(Ω)→Xh, v¯h=Ghv where v¯h is the unique solution in Xh of
(3.8)δ2(a(v¯h,ϕ)+J0(v¯h,ϕ))+(v¯h,ϕ)=(v,ϕ)∀ϕ∈Xh.Remark 3.5.
An attractive alternative is to define the differential filter by a discrete Stokes problem so as to preserve incompressibility approximately [23–25].Definition 3.6.
The discrete van Cittert deconvolution operatorsGNh are
(3.9)GNhv:=∑n=0N(Πh-Gh)nv,
where Πh:L2(Ω)→Xh is the L2 projection.Forv∈Xh, the discrete deconvolution operator for N=0,1,2 is(3.10)G0hv=v,G1hv=2v-v¯h,G2hv=3v-3v¯h+v¯h¯h.GN was shown to be an O(δ2N+2) approximate inverse to the filter operator G in Lemma 2.1 of Dunca and Epshteyn [26], recalled next.Lemma 3.7.
GN is a bounded, self-adjoint positive operator. GN is an O(δ2N+2) asymptotic inverse to the filter G. Specifically, for smooth ϕ and as δ→0,
(3.11)ϕ=GNϕ¯+(-1)(N+1)δ2N+2ΔN+1G(N+1)ϕ.Some basic facts about discrete differential filters and deconvolution operators are presented next.Lemma 3.8.
Forv∈L2(Ω), one has the following bounds for the discretely filtered and approximately deconvolved  v:(3.12)∥v¯h∥0,Ω≤∥v∥0,Ω,(3.13)∥GNhv¯h∥0,Ω≤C(N)∥v∥0,Ω.Proof.
The proof of (3.12) follows from the standard finite element techniques applied on the discretized equation (3.8) of the filter problem (1.1). Pick ϕ=v¯h. Then, using coercivity result
(3.14)δ2κ∥v¯h∥X2+∥v¯h∥0,Ω2≤(v,v¯h)Ω≤12∥v∥0,Ω2+12∥v¯h∥0,Ω2.
The term δ2κ∥v¯h∥X2 is positive, so it will be dropped, which yields
(3.15)12∥v¯h∥0,Ω2≤12∥v∥0,Ω2.
Multiplying by 2 and taking the square root yields the estimate (3.12). Equation (3.13) follows immediately from (3.12) and the definition of GNh.Lemma 3.9.
For smoothϕ the discrete approximate deconvolution operator satisfies
(3.16)∥ϕ-GNhϕ¯h∥0,Ω≤C1δ2N+2∥Gϕ∥2N+2,Ω+C2(δhr+hr+1)(∑n=2N+1|Gnϕ|r+1,Ω).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:
(3.17)∥ϕ-GNhϕ¯h∥0,Ω≤∥ϕ-GNϕ¯∥0,Ω+∥GNϕ¯-GNhϕ¯∥0,Ω+∥GNhϕ¯-GNhϕ¯h∥0,Ω.
Lemma 3.7 gives
(3.18)∥ϕ-GNϕ¯∥0,Ω≤Cδ2N+2∥ϕ¯∥2N+2,Ω.
The standard discontinuous finite element bound for (3.8) is given by (cf. Rivière [13])
(3.19)∥ϕ¯-ϕ¯h∥0,Ω≤C(δhr+hr+1)|ϕ¯|r+1,Ω.
Lemma3.8 gives for the third term in (3.17) that ∥GNhϕ¯-GNhϕ¯h∥≤C∥ϕ¯-ϕ¯h∥. Then, (3.19) is applied.
Now, it is left to bound the second term from (3.17). First, note that for N=0, ∥G0ϕ¯h-G0hϕ¯h∥0,Ω=0. Based on Definition 3.6 of continuous and discrete deconvolution operators and their expansion (see (3.10)), GN is a polynomial of degree N in G (and GNh in Gh as well). Thus, the second term in (3.17) can be written as(3.20)∥GNϕ¯-GNhϕ¯∥0,Ω=∥∑n=0Nαn(Gnϕ¯-(Gh)nϕ¯)∥0,Ω≤∑n=0Nαn∥Gnϕ¯-(Gh)nϕ¯∥0,Ω.
ForO(1) coefficients αn and for N=1, the result (3.19) gives(3.21)∥Gϕ¯-(Gh)ϕ¯∥0,Ω=∥ϕ¯¯-ϕ¯¯h∥0,Ω≤C(δhr+hr+1)(|ϕ¯¯|r+1,Ω).
ForN=2, the results (3.19) and (3.12) give (3.22)∥G2ϕ¯-(G2h)ϕ¯∥0,Ω=∥ϕ¯¯¯-ϕ¯h¯h¯h∥0,Ω≤∥ϕ¯¯¯-ϕ¯¯¯h∥0,Ω+∥ϕ¯¯¯h-ϕ¯h¯h¯h∥0,Ω≤∥ϕ¯¯¯-ϕ¯¯¯h∥0,Ω+∥ϕ¯¯-ϕ¯¯h∥0,Ω≤C(δhr+hr+1)(|ϕ¯¯¯|r+1,Ω+|ϕ¯¯|r+1,Ω).
Inductively,(3.23)∥GNϕ¯-GNhϕ¯∥0,Ω≤C(δhr+hr+1)(∑n=2N+1|Gnϕ|r+1,Ω).
The proof is completed by combining the derived bounds for the terms in (3.17).Remark 3.10.
There remains the question of uniform inδ bound of the last term, |(G)nϕ|r+1,Ω, in (3.16). This is a question about uniformregularity of an elliptic-elliptic singular perturbation problem and some results are proven in [28] by Layton. To summarize, in the periodic case it is very easy to show by Fourier series that for all k(3.24)|Gϕ|r+1,Ω≤C|ϕ|r+1,Ωandthus|Gnϕ|r+1,Ω≤C|ϕ|r+1,Ω.
The nonperiodic case can be more delicate. Suppose ∂Ω∈Cr+3 and ϕ=0 on ∂Ω (i.e., ϕ∈H01(Ω)∩Hr+1(Ω)). Call Gϕ=ϕ¯ so ϕ¯ satisfies
(3.25)-δ2Δϕ¯+ϕ¯=ϕinΩ,ϕ¯=0on∂Ω.
Then it is known that ϕ¯∈Hr+3(Ω)∩H01(Ω), and Δϕ¯=0 on ∂Ω. Further,
(3.26)∥ϕ¯∥j,Ω≤C∥ϕ∥j,Ω,j=0,1,2.
So, (3.24) holds for r=-1,0,+1. It also holds for higher values of r provided additionally Δjϕ=0 on ∂Ω for 0≤j≤[(r+1)/2]-1.
Now consider the second termn=2, that is,  G2ϕ=ϕ¯¯. We know from elliptic theory for ϕ∈Hr+1(Ω)⋂H01(Ω), that ϕ¯∈Hr+3(Ω)⋂H01(Ω), (as noted above) Δϕ¯=0 on ∂Ω and(3.27)-δ2Δϕ¯¯+ϕ¯¯=ϕ¯inΩ,ϕ¯=Δϕ¯=0on∂Ω.
Theorem 1.1 in [28] then implies, uniformly in δ,
(3.28)∥ϕ¯¯∥j,Ω≤C∥ϕ¯∥j,Ω,j=0,1,2,3,4.
This extends directly to Gnϕ.Extending Lemma3.9, the following assumption will be made.Assumption DG1.
The|Gn(ϕ)|r+1,Ω terms in (3.16) are independent of δ and
(3.29)∥ϕ-GNhϕ¯h∥≤C1δ2N+2∥ϕ∥2N+2,Ω+C2(δhr+hr+1)|ϕ|r+1,Ω.The minimal conditions that are assumed throughout are that the (discrete) filter and (discrete) deconvolution used satisfy the following consistency conditions of Stanculescu [29].Assumption DG2.
Gh and (I-GNhGh) are symmetric, positive definite (SPD) operators.These have been proven to hold for van Cittert deconvolution (cf. Stanculescu [29], Manica and Merdan [22] and Layton et al.  [27]). For the DG method, the second assumption restricts our parameter ϵa=-1 for the discretization of the filter problem (3.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Δt denote the time step such that M=T/Δt is a positive integer. Let 0=t0<t1<⋯<tM=T be a subdivision of the interval (0,T). The function ϕ evaluated at the time tm is denoted by ϕm. With the above forms, the fully discrete scheme is: find (unh,pnh)n≥0∈Xh×Qh  such that:(3.30)1Δt(un+1h-unh,v)Ω+ν(a(un+1h,v)+J0(un+1h,v))+χ(un+1h-GNhun+1h¯h,v)Ω+cunh(unh;un+1h,v)+b(v,pn+1h)=(fn+1,v)Ω∀v∈Xh,(3.31)b(un+1h,q)=0∀q∈Qh,(3.32)(u0h,v)Ω=(u0,v)Ω∀v∈Xh.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 compute [30].The consistency result of the semidiscrete scheme is showed next.Lemma 3.12 (Consistency).
Let(u,p) be the solution to (1.7)–(1.10), then (u,p) satisfies
(3.33)(ut,v)Ω+ν(a(u,v)+J0(u,v))+χ(u-GNhu¯h,v)Ω+cu(u;u,v)+b(v,p)=(f,v)Ω+E(u,p,f;v)∀v∈Xh,(3.34)b(u,q)=0∀q∈Qh,(3.35)(u(0),v)Ω=(u0,v)Ω∀v∈Xh,
where the consistency error E(u,v)=χ(u-GNhu¯h,v)Ω-χ(u-GNu¯,v)Ω. Furthermore,
(3.36)E(u,v)≤C(δhr+hr+1)|u|r+1,Ω∥v∥0,Ω.Proof.
Equations (3.34) and (3.35) are clearly satisfied because of (1.8), (1.9), and (1.10) and the regularity of u. Next, we multiply (1.7) by v and integrate over one mesh element E(3.37)(ut,v)E+(∇·(uu),v)E-ν(Δu,v)E+χ(u-GNu¯,v)E+(∇p,v)E=(f,v)E.
Summing over all elements E(3.38)∑E∈ℰh(ut,v)E+(∇·(uu),v)Ω-ν∑E∈ℰh(Δu,v)E+χ(u-GNu¯,v)Ω+(∇p,v)Ω=(f,v)Ω.
By Green's formula(3.39)-ν∑E∈ℰh(Δu,v)E=-ν(Δu,v)Ω=ν∑E∈ℰh(∇u,∇v)E-ν∑e∈Γh∪∂Ω∫e(∇u)ne·[v].
The regularity of u then yields
(3.40)-ν∑E∈ℰh(Δu,v)E=ν(a(u,v)+J0(u,v)).
Note that Green's formula yields
(3.41)(∇p,v)Ω=b(v,p),
and that the incompressibility condition with the regularity of u yields
(3.42)(u·∇u,v)Ω=cu(u;u,v).
The final result is obtained by bounding the consistency error E(u,v). For N=0, we have
(3.43)E(u,v)≤∥u¯-u¯h∥0,Ω∥v∥0,Ω≤C(δhr+hr+1)|u¯|r+1,Ω∥v∥0,Ω.
For N=1,
(3.44)E(u,v)≤∥Gu¯-(Gh)u¯h∥0,Ω∥v∥0,Ω≤∥u¯¯-u¯h¯h∥0,Ω∥v∥0,Ω≤(∥u¯¯-u¯¯h∥0,Ω+∥u¯¯h-u¯h¯h∥0,Ω)∥v∥0,Ω≤(∥u¯¯-u¯¯h∥0,Ω+∥u¯-u¯h∥0,Ω)∥v∥0,Ω≤C(δhr+hr+1)(|u¯¯|r+1,Ω+|u¯|r+1,Ω)∥v∥0,Ω.
The bound for E(u,v) is obtained by applying an induction argument and Remark 3.10.The existence and uniqueness of the discrete solution is stated next.Proposition 3.13.
Assume that Lemma3.1 holds. Then, there exists a unique solution to (3.30)–(3.32).Proof.
The existence ofu0h is trivial. Given unh, the problem of finding a unique un+1h satisfying (3.30)-(3.31) is linear and finite-dimensional. Therefore, it suffices to show uniqueness of the solution. Consider the problem restricted to the subspace Vh defined by
(3.45)Vh={v∈Xh:b(v,q)=0∀q∈Qh}.
Let un+1h and ûn+1h be two solutions and let θn+1=un+1h-ûn+1h. Then, θn+1 satisfies:
(3.46)1Δt(θn+1,v)Ω+ν(a(θn+1,v)+J0(θn+1,v))+cunh(unh;θn+1,v)+χ(θn+1-GNhθn+1¯h,v)=0∀v∈Vh.
Choosing v=θn+1 and using the coercivity result (3.5), positivity result (3.7) and positivity of the operator I-GNhGh given in Assumption DG2, we obtain:
(3.47)1Δt∥θn+1∥0,Ω2+νκ∥θn+1∥X2≤0,
which yields that θn+1=0. The existence and uniqueness of the pressure pn+1h is then obtained from the inf-sup condition (3.6).Some approximation properties of the spacesXh and Qh are recalled next. From Crouzeix and Raviart [31], and Girault et al. [16], for each integer r≥1, and for any v∈(H01(Ω))2, there is a unique discrete velocity ṽ∈Xh such that
(3.48)b(v-ṽ,q)=0∀q∈Qh.
Furthermore, if v∈(H01(Ω))2∩(Hr+1(Ω))2, there is a constant C independent of h such that(3.49)∥v-ṽ∥X≤Chr|v|r+1,Ω,(3.50)|v-ṽ|m,Ω≤Chr+1-m|v|r+1,Ω,m=0,1.For the pressure space, we use the approximation given by theL2 projection. For any q∈L02(Ω), there exists a unique discrete pressure q̃∈Qh such that(3.51)(q-q̃,z)Ω=0∀z∈Qh.
In addition, if q∈Hr(Ω), then(3.52)∥q-q̃∥m,E≤Chr-m|q|r,E,∀E∈ℰh,m=0,1,2.The discrete Gronwall's lemma plays an important role in the following analysis.Lemma 3.14 (Discrete Gronwall's Lemma (cf. Heywood and Rannacher [32])).
LetΔt, H, and an, bn, cn, γn (for integers n≥0) be nonnegative numbers such that
(3.53)al+Δt∑n=0lbn≤Δt∑n=0lγnan+Δt∑n=0lcn+Hforl≥0.
Suppose that Δtγn<1, for all n, and set σn=(1-Δtγn)-1. Then,
(3.54)al+Δt∑n=0lbn≤exp(Δt∑n=0lσnγn){Δt∑n=0lcn+H}forl≥0.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 equations(w,P) that is approximated is a strong solution and in particular satisfies (cf. Rivière [13])(3.55)(wt,v)Ω+ν(a(w,v)+J0(w,v))+cw(w;w,v)+b(v,P)=(f,v)Ω∀v∈Xh,(3.56)b(w,q)=0∀q∈Qh,(3.57)(w0h,v)Ω=(u0,v)Ω∀v∈Xh.


## 4. A Priori Error Estimates

In this section, convergence of the scheme (3.30)–(3.32) is proved. Optimal error estimates in the energy norm are obtained.Theorem 4.1.
Assume thatw∈l2(0,T;(Hr+1(Ω))2)∩l2(0,T;(H2N+2(Ω))2, wt∈l2(0,T;(Hr+1(Ω))2)∩L∞((0,T)×Ω), wtt∈L2(0,T;(H1(Ω))2), p∈l2(0,T;Hr(Ω)) and u0∈(Hr+1(Ω))2. Assume also that the coercivity Lemma 3.1 holds and that w satisfies (3.24). If δ is chosen of the order of h, and Δt<1, there exists a constant C, independent of h and Δt but dependent on ν-1 such that the following error bound holds, for any 1≤m≤M:
(4.1)∥wm-umh∥0,Ω+(νκΔt∑n=1m∥wn-unh∥X2)1/2≤Chr(ν-1+ν+χ+1)+CΔt+Cχh2N+2.Remark 4.2.
The dependence of these error estimates with respect to the Reynolds number Re (~ 1/ν) is 𝒪(ReeRe), which is an improvement with respect to the continuous finite element method where the dependence is 𝒪(ReeRe3).Proof.
Definingen=uh(tn)-w(tn) and subtracting (3.55) from (3.30), we have
(4.2)1Δt(en+1-en,v)Ω-(wt(tn+1),v)Ω+ν(a+J0)(en+1,v)+χ(en+1-GNhen+1¯h,v)Ω+cunh(unh;un+1h,v)-cwn+1(wn+1;wn+1,v)+b(v,pn+1h-Pn+1)=-1Δt(wn+1-wn,v)Ω-χ(wn+1-GNhwn+1¯h,v)Ω∀v∈Xh,(4.3)b(en+1,q)=0∀q∈Qh.
We now decompose the error en=ϕn-ηn, where ϕn=unh-w̃n and ηn is the interpolation error ηn=wn-w̃n. Choosing v=ϕn+1 in the equation above, using the coercivity result (3.5) and positivity of the operator I-GNhGh, we obtain for (4.2)
(4.4)12Δt(∥ϕn+1∥0,Ω2-∥ϕn∥0,Ω2)+νκ∥ϕn+1∥X2+cunh(unh;un+1h,ϕn+1)-cwn+1(wn+1;wn+1,ϕn+1)≤(ηt(tn+1),ϕn+1)Ω+ν(a+J0)(ηn+1,ϕn+1)+(w̃t(tn+1)-1Δt(w̃n+1-w̃n),ϕn+1)Ω+χ(ηn+1-GNhηn+1¯h,ϕn+1)Ω-χ(wn+1-GNhwn+1¯h,ϕn+1)Ω+b(ϕn+1,Pn+1-pn+1h).
Consider now the nonlinear terms from the above equation. First note that since w is  continuous, we can rewrite
(4.5)cwn+1(wn+1;wn+1,ϕn+1)=cunh(wn+1;wn+1,ϕn+1),
so, for readability, the superscript wnh in the c form is dropped. Therefore, adding and subtracting the interpolant w̃n+1 yields
(4.6)cunh(unh,un+1h,ϕn+1)-cunh(wn+1,wn+1,ϕn+1)=c(unh,ϕn+1,ϕn+1)-c(ϕn,ηn+1,ϕn+1)+c(ϕn,wn+1,ϕn+1)-c(ηn,w̃n+1,ϕn+1)-c(wn,ηn+1,ϕn+1)-c(wn+1-wn,wn+1,ϕn+1).
Thus, the error equation (4.4) is rewritten as
(4.7)12Δt(∥ϕn+1∥0,Ω2-∥ϕn∥0,Ω2)+νκ∥ϕn+1∥X2+c(unh,ϕn+1,ϕn+1)≤|c(ϕn,ηn+1,ϕn+1)|+|c(ϕn,wn+1,ϕn+1)|+|c(ηn,w̃n+1,ϕn+1)|+|c(wn,ηn+1,ϕn+1)|+|c(wn+1-wn,wn+1,ϕn+1)|+|(ηt(tn+1),ϕn+1)Ω|+|(w̃t(tn+1)-1Δt(w̃n+1-w̃n),ϕn+1)Ω|+|ν(a+J0)(ηn+1,ϕn+1)|+|χ(ηn+1-GNhηn+1¯h,ϕn+1)Ω|+|χ(wn+1-GNhwn+1¯h,ϕn+1)Ω|+|b(ϕn+1,Pn+1-pn+1h)|≤|T0|+|T1|+⋯+|T10|.
From property (3.7), the term c(wnh;ϕn+1,ϕn+1) in the left-hand side of (4.7) is positive and therefore it will be dropped. For the other terms of the form 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 in [19] by Kaya and Rivière. The constant C is a generic constant that is independent of h,ν and Δt, and that takes different values at different places
(4.8)|T0|=|c(ϕn,ηn+1,ϕn+1)|≤νκ20∥ϕn+1∥X2+Cν∥ϕn∥0,Ω2,|T1|=|c(ϕn,wn+1,ϕn+1)|≤νκ20∥ϕn+1∥X2+Cν∥ϕn∥0,Ω2,|T2|=|c(ηn,w̃n+1,ϕn+1)|≤νκ20∥ϕn+1∥X2+Cνh2r|wn|2r,Ω2,|T3|=|c(wn,ηn+1,ϕn+1)|≤νκ20∥ϕn+1∥X2+Cνh2r|wn|2r,Ω2,|T4|=|c(wn+1-wn,wn+1,ϕn+1)|≤νκ20∥ϕn+1∥X2+CνΔt2∥wt∥L∞([tn,tn+1]×Ω)2.
Therefore, we have
(4.9)|T0|+⋯+|T4|≤5νκ20∥ϕn+1∥X2+Cν-1∥ϕn∥0,Ω2+Cν-1h2r|wn|r+1,Ω2+Cν-1Δt2∥wt∥L∞([tn,tn+1]×Ω)2.
To boundT5, Cauchy-Schwarz's inequality, Young's inequality and the approximation result (3.50) for wt are applied(4.10)|T5|≤∥ϕn+1∥0,Ω∥ηt(tn+1)∥0,Ω≤18∥ϕn+1∥0,Ω2+Ch2r+2|wt(tn+1)|k+1,Ω2.
To bound the term T6, a Taylor expansion with integral remainder is used
(4.11)w̃n=w̃n+1-Δtw̃t(tn+1)+12∫tntn+1(s-tn)w̃tt(s)ds.
This implies that
(4.12)∥w̃t(tn+1)-w̃n+1-w̃nΔt∥0,Ω2≤Δt6∫tntn+1∥w̃tt(s)∥0,Ω2ds.
Thus, with (3.50), we have
(4.13)|T6|≤18∥ϕn+1∥0,Ω2+CΔt∫tntn+1∥w̃tt(s)∥0,Ω2ds≤18∥ϕn+1∥0,Ω2+CΔt∫tntn+1∥wtt(s)∥0,Ω2ds.
Next, the term T7 is expanded as
(4.14)|T7|≤|ν∑E∈ℰh∫E∇ηn+1:∇ϕn+1|+|ν∑e∈Γh∪∂Ω∫e{∇ηn+1}ne·[ϕn+1]|+|νϵa∑e∈Γh∪∂Ω∫e{∇ϕn+1}ne·[ηn+1]|+|νJ0(ηn+1,ϕn+1)|=|T71|+|T72|+|T73|+|T74|.
The term T71 is bounded using Cauchy-Schwarz inequality, Young's inequality and the approximation result (3.49)
(4.15)|T71|≤ν∥ϕn+1∥X∥ηn+1∥X≤νκ20∥ϕn+1∥X2+Cν∥ηn+1∥X2≤νκ20∥ϕn+1∥X2+Cνh2r|wn+1|r+1,Ω2.
Using Cauchy-Schwarz's inequality, trace inequality (2.8) and approximation result (3.50), we have
(4.16)|T72|≤ν∑e∈Γh∪∂Ω∥{∇ηn+1}ne∥0,e∑e∈Γh∪∂Ω∥[ϕn+1]∥0,e≤Cν(∑e∈Γh∪∂Ω1|e|∥[ϕn+1]∥0,e2)1/2(|∥∇ηn+1∥|0,Ω+h|∥∇2ηn+1∥|0,Ω)≤νκ20∥ϕn+1∥X2+Cνh2r|wn+1|r+1,Ω2.
Using Cauchy-Schwarz's inequality, trace inequality (2.10), and approximation result (3.49), we have
(4.17)|T73|≤ν(∑e∈Γh∪∂Ω∥{∇ϕn+1}ne∥0,e2)1/2(∑e∈Γh∪∂Ω∥[ηn+1]∥0,e2)1/2≤Cν∥ϕn+1∥X(∑e∈Γh∪∂Ω1|e|∥[ηn+1]∥0,e2)1/2≤νκ20∥ϕn+1∥X2+Cνh2r|wn+1|r+1,Ω2.
Using the approximation result (3.49), we have
(4.18)|T74|≤ν(∑e∈Γh∪∂Ωσ|e|∥[ϕn+1]∥0,e2)1/2(∑e∈Γh∪∂Ωσ|e|∥[ηn+1]∥0,e2)1/2≤Cν∥ϕn+1∥X∥ηn+1∥X≤νκ20∥ϕn+1∥X2+Cνh2r|wn+1|r+1,Ω2.
Putting together the bounds (4.15), (4.16), (4.17) and (4.18),
(4.19)|T7|≤4νκ20∥ϕn+1∥X2+Cνh2r|wn+1|r+1,Ω2.
For the terms T8 and T9, Cauchy-Schwarz's inequality, Young's inequality and bounds (3.13) and (3.29) are applied
(4.20)|T8|=|χ(ηn+1-GNhηn+1¯h,ϕn+1)Ω|≤χ∥ηn+1-GNhηn+1¯h∥0,Ω∥ϕn+1∥0,Ω≤Cχ2∥ηn+1-GNhηn+1¯h∥0,Ω2+18∥ϕn+1∥0,Ω2≤C(N)χ2∥ηn+1∥0,Ω2+18∥ϕn+1∥0,Ω2≤C(N)χ2h2r+2|wn+1|r+1,Ω2+18∥ϕn+1∥0,Ω2,|T9|=|χ(wn+1-GNhwn+1¯h,ϕn+1)Ω|≤χ∥wn+1-GNhwn+1¯h∥0,Ω∥ϕn+1∥0,Ω≤Cχ2∥wn+1-GNhwn+1¯h∥0,Ω2+18∥ϕn+1∥0,Ω2≤Cχ2δ4N+4∥wn+1∥2N+2,Ω2+Cχ2(δ2h2r+h2r+2)|wn+1|r+1,Ω2+18∥ϕn+1∥0,Ω2.
Using (4.3) with (3.48) and (3.51), the pressure term T10 is reduced to
(4.21)|T10|=|b(ϕn+1,pn+1-p̃n+1)+b(ϕn+1,p̃n+1-pn+1h)|=|b(ϕn+1,pn+1-p̃n+1)|=|∑e∈Γh∫e{pn+1-p̃n+1}[ϕn+1]·ne|,
which is bounded by using Cauchy-Schwarz's inequality, Young's inequality, trace inequality (2.7) and approximation result (3.52)
(4.22)|T10|≤C(∑e∈Γh∪∂Ω1|e|∥[ϕn+1]∥0,e)1/2(∥pn+1-p̃n+1∥0,Ω+h|∥∇pn+1-∇p̃n+1∥|0,Ω)≤νκ20∥ϕn+1∥X2+Cν-1h2r|pn+1|r,Ω2.
With the bounds (4.9), (4.10), (4.13), (4.19), (4.20), and (4.22), the error equation becomes
(4.23)12Δt(∥ϕn+1∥0,Ω2-∥ϕn∥0,Ω2)+νκ2∥ϕn+1∥X2≤+Cν-1∥ϕn∥0,Ω2+12∥ϕn+1∥0,Ω2+Cν-1Δt2∥wt∥L∞([tn,tn+1]×Ω)2+Ch2r+2|wt(tn+1)|k+1,Ω2+CΔt∫tntn+1∥wtt(s)∥0,Ω2ds+Cν-1h2r|wn|r+1,Ω2+Cνh2r|wn+1|r+1,Ω2+Cχ2h2r+2∥wn+1∥r+1,Ω2+C(χ2(δ2h2r+h2r+2))|wn+1|r+1,Ω2+Cχ2δ4N+4∥wn+1∥2N+2,Ω2+Cν-1h2r|pn+1|r,Ω2,
where C are constants independent of h,ν and Δt. Now, multiply the equation by 2Δt, sum from n=0 to n=m-1, and use the assumption that δ is chosen of the order of h to obtain
(4.24)(1-Δt)∥ϕm∥0,Ω2+νκΔt∑n=0m-1∥ϕn+1∥X2≤∥ϕ0∥0,Ω2+C(ν-1+1)∑n=0m-1∥ϕn∥0,Ω2+Ch2rΔt∑n=0m-1(ν-1|wn|r+1,Ω2+(ν+χ2)|wn+1|r+1,Ω2+|wt(tn+1)|r+1,Ω2+ν-1|pn+1|r,Ω2)+CΔt2∫0T∥wtt(s)∥1,Ω2ds+CΔt2ν-1∥wt∥L∞([0,T]×Ω)2+CΔtχ2h4N+4∥wn+1∥2N+2,Ω2.
Thus, using Gronwall's lemma with Δt<1, there is a constant C independent of h and Δt, but dependent on ν-1, such that
(4.25)∥ϕm∥0,Ω2+νκΔt∑n=1m∥ϕn∥X2≤∥ϕ0∥0,Ω2+Ch2r(ν-1+ν+χ2+1)+CΔt2+Cχ2h4N+4.
The final result is then obtained by noting that the term ∥ϕ0∥0,Ω2 is of order h2r and by using triangle inequality and approximation result.


## 5. Conclusion

In this paper, a numerical scheme for solving the time relaxation family of models based on approximate deconvolution technique for fluid flow problems is formulated and analyzed. The proposed method is convergent with optimal convergence rates with respect to the mesh size. The approximations of the average velocity and pressure are discontinuous piecewise polynomials. One benefit of using discontinuous elements is that the error estimates depend on the Reynolds number as𝒪(ReeRe), whereas the dependence is 𝒪(ReeRe3) for continuous finite elements (cf. [14]).The proposed scheme (3.30)–(3.32) contains parameter ϵa∈{-1,0,1} in the bilinear form a(·,·) for the discretization of the viscous term that yields different numerical approximations. Only numerical simulations of benchmark problems for high Reynolds numbers, will help determine which choices of ϵa are preferred for a given mesh size. This is the object of a future paper.



---

*Source: 419021-2010-10-03.xml*