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

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 𝒪(Re𝑒^Re), which is an improvement with respect to the continuous finite element method where the dependence is 𝒪(Re𝑒^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 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 𝒪(Re^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 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])𝜙=𝐺𝜙inΩ,𝜙=𝟎on𝜕Ω, (1.1)(1.2) where 𝐺⁻¹=−𝛿²Δ+𝐼, Ω 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 𝑁=0,1,…, it computes an approximate solution 𝐮_𝑁 by N steps of a fixed point iteration for the fixed point problem (cf. Bertero and Boccacci [10]) given𝐮solve𝐮=𝐮+{𝐮−𝐺𝐮},for𝐮. (1.3) The deconvolution approximation is then computed as follows.

Algorithm 1.1 (van Cittert approximate deconvolution algorithm). Consider that 𝐮₀=𝐮, for 𝑛=1,2,…,𝑁−1, perform 𝐮_𝑛+1=𝐮_𝑛+{𝐮−𝐺𝐮_𝑛}. (1.4)

By eliminating the intermediate steps, it is easy to find an explicit formula for the 𝑁th deconvolution operator 𝐺_𝑁 given by𝐺_𝑁𝐮∶=^𝑁_𝑛=0(𝐼−𝐺)^𝑛𝐮. (1.5) For example, the approximate deconvolution operator corresponding to 𝑁=0,1,2 is 𝐺₀𝐮=𝐺𝐮,₁𝐮=2𝐮−𝐺𝐮,₂𝐮=3𝐮−3𝐮+𝐮. (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, 12] and arises by adding a linear, lower order time regularization term, 𝜒𝐮^′=𝜒(𝐮−𝐺_𝑁𝐮), to the Navier-Stokes equations (NSE). This term involves 𝐮^′ 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𝐮_𝑡+𝐮⋅∇𝐮−𝜈Δ𝐮+∇𝑝+𝜒(𝐮−𝐺_𝑁𝐮)=𝐟in(0,𝑇]×Ω,∇⋅𝐮=0in[0,𝑇]×Ω,𝐮=𝟎in[0,𝑇]×𝜕Ω,𝐮(0,⋅)=𝐮₀(⋅)inΩ, (1.7)(1.8)(1.9)(1.10) where Ω⊂ℝ², is a convex bounded regular domain with boundary 𝜕Ω, 𝐮 is the fluid velocity, 𝑝 is the fluid pressure and 𝐟 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 𝐮₀. A pressure normalization condition ∫_Ω𝑝=0 is also needed for uniqueness of the pressure. The time relaxation coefficient 𝜒 has units 1/𝑡𝑖𝑚𝑒. 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 𝐮 and pressure 𝑝 are discontinuous piecewise polynomials of degree 𝑟 and 𝑟−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. Section 2 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 ℎ 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 ℰ_ℎ of Ω, consisting of triangles of maximum diameter ℎ, is introduced. Let ℎ_𝐸 denote the diameter of a triangle 𝐸 and 𝜌_𝐸 the diameter of its inscribed circle. Regulary, it is meant that there exists a parameter 𝜁>0, independent of ℎ, such that ℎ_𝐸𝜌_𝐸=𝜁_𝐸≤𝜁,∀𝐸∈ℰ_ℎ. (2.1) This assumption will be used throughout this work. Γ_ℎ denotes the set of all interior edges of ℰ_ℎ. Let 𝑒 denote a segment of Γ_ℎ shared by two triangles 𝐸^𝑘 and 𝐸^𝑙 (𝑘<𝑙) of ℰ_ℎ; it is associated with 𝑒 a specific unit normal vector 𝐧_𝑒 directed from 𝐸^𝑘 to 𝐸^𝑙 and the jump and average of a function 𝜙 on 𝑒 is formally defined by[𝜙]=(𝜙|_𝐸𝑘)|_𝑒−(𝜙|_𝐸𝑙)|_𝑒1,{𝜙}=2(𝜙|_𝐸𝑘)|_𝑒+12(𝜙|_𝐸𝑙)|_𝑒. (2.2) If 𝑒 belongs to the boundary 𝜕Ω, then 𝐧_𝑒 is the unit normal 𝐧 exterior to Ω and the jump and the average of 𝜙 on 𝑒 coincide with the trace of 𝜙 on 𝑒.

Here, for any domain 𝒪, 𝐿²(𝒪) is the classical space of square-integrable functions with inner-product (𝑓,𝑔)_𝒪=∫_𝒪𝑓𝑔 and norm ‖⋅‖_0,𝒪. The space 𝐿²₀(Ω) is the subspace of functions of 𝐿²(Ω) with zero mean value 𝐿²₀(Ω)={𝐯∈𝐿²(Ω)∶_Ω𝐯=0}. (2.3) The standard Sobolev spaces are denoted by 𝑊^𝑟_𝑝(Ω) where 𝑊^𝑟₂(Ω) is the 𝐻^𝑟(Ω), with norm ‖⋅‖_𝑟,Ω and seminorm |⋅|_𝑟,Ω.

Next, the discrete velocity and pressure spaces are defined to be consisting of discontinuous piecewise polynomials. For any positive integer 𝑟, the corresponding finite-dimensional spaces are 𝐗^ℎ={𝐯∈(𝐿²(Ω))²∶∀𝐸∈ℰ_ℎ,𝐯∈(ℙ_𝑟(𝐸))²𝑄},_ℎ={𝑞∈𝐿²₀(Ω)∶∀𝐸∈ℰ_ℎ,𝑞∈ℙ_𝑟−1(𝐸)}, (2.4) where 𝑃_𝑟(𝐸)=span{𝑥^𝑖𝑦^𝑗∶𝑖+𝑗≤𝑟} is defined as the span of polynomials of order 𝑟 over triangle 𝐸.

Denoting by |𝑒| the measure of 𝑒, the following norms are associated for the spaces 𝐗^ℎ and 𝑄^ℎ‖𝐯‖_𝑋=(|‖∇𝐯‖|²_0,Ω+_{𝑒∈Γℎ∪𝜕Ω}1|𝑒|‖[𝐯]‖²_0,𝑒)^1/2,‖𝑞‖_𝑄=‖𝑞‖_0,Ω, (2.5) where |‖𝐯‖|_0,Ω is the broken norm defined by |‖𝐯‖|_0,Ω=(_{𝐸∈ℰℎ}‖𝐯‖²_0,𝐸)^1/2. (2.6)

Finally, some trace and inverse inequalities are recalled, that hold true on each element 𝐸 in ℰ_ℎ, with diameter ℎ_𝐸, the constant 𝐶 is independent of ℎ_𝐸‖𝐯‖_0,𝑒≤𝐶(ℎ_𝐸^−1/2‖𝐯‖_0,𝐸+ℎ_𝐸^1/2‖∇𝐯‖_0,𝐸),∀𝑒∈𝜕𝐸,∀𝐯∈(𝐻¹(𝐸))²,‖∇𝐯‖_0,𝑒≤𝐶(ℎ_𝐸^−1/2‖∇𝐯‖_0,𝐸+ℎ_𝐸^1/2‖∇²𝐯‖_0,𝐸),∀𝑒∈𝜕𝐸,∀𝐯∈(𝐻²(𝐸))²,‖𝐯‖_0,𝑒≤𝐶ℎ_𝐸^−1/2‖𝐯‖_0,𝐸,∀𝑒∈𝜕𝐸,∀𝐯∈𝐗^ℎ,‖∇𝐯‖_0,𝑒≤𝐶ℎ_𝐸^−1/2‖∇𝐯‖_0,𝐸,∀𝑒∈𝜕𝐸,∀𝐯∈𝐗^ℎ. (2.7)(2.8)(2.9)(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 𝑎∶𝐗^ℎ×𝐗^ℎ→ℝ, and 𝐽₀∶𝐗^ℎ×𝐗^ℎ→ℝ by 𝑎(𝐳,𝐯)=_{𝐸∈ℰℎ𝐸}∇𝐳∶∇𝐯−_{𝑒∈Γℎ∪𝜕Ω𝑒}{∇𝐳}𝐧_𝑒⋅[𝐯]+𝜖_{𝑎𝑒∈Γℎ∪𝜕Ω𝑒}{∇𝐯}𝐧_𝑒𝐽⋅[𝐳],₀(𝐳,𝐯)=_{𝑒∈Γℎ∪𝜕Ω}𝜎|𝑒|_𝑒[𝐳]⋅[𝐯]. (3.1)

The parameter 𝜖_𝑎 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 ‖⋅‖_𝑋. In the case where 𝜖_𝑎=−1, the bilinear form 𝑎 is symmetric; otherwise it is nonsymmetric. We remark that the form 𝑎(𝐮,𝐯) is the standard primal DG discretization of the operator −Δ𝐮. Finally, if 𝜖_𝑎 is either −1 or 0, the jump parameter 𝜎 should be chosen sufficiently large to obtain coercivity of 𝑎 (see Lemma 3.1). If 𝜖_𝑎=1, then the jump parameter 𝜎 is taken equal to 1.

The incompressibility condition (1.8) is enforced by means of the bilinear form 𝑏∶𝐗^ℎ×𝑄^ℎ→ℝ defined by 𝑏(𝐯,𝑞)=−_{𝐸∈ℰℎ𝐸}𝑞∇⋅𝐯+_{𝑒∈Γℎ∪𝜕Ω𝑒}{𝑞}[𝐯]⋅𝐧_𝑒. (3.2) Finally, the DG discretization of the nonlinear convection term 𝐰⋅∇𝐰, which was introduced in [16] by Girault et al. and studied extensively in [16, 17] by the same authors, is recalled as follows: 𝑐^𝐳(𝐮;𝐯,𝐭)=_{𝐸∈ℰℎ}(_𝐸1(𝐮⋅∇𝐯)⋅𝐭+2_𝐸1(∇⋅𝐮)𝐯⋅𝐭)−2_{𝑒∈Γℎ∪𝜕Ω𝑒}[𝐮]⋅𝐧_𝑒{𝐯⋅𝐭}+_{𝐸∈ℰℎ𝜕𝐸−}|{𝐮}⋅𝐧_𝐸|(𝐯^int−𝐯^ext)⋅𝐭^int, (3.3) where 𝜕𝐸₋={𝐱∈𝜕𝐸∶{𝐳}⋅𝐧_𝐸<0}, (3.4) the superscript 𝐳 denotes the dependence of 𝜕𝐸₋ on 𝐳 and the superscript int (resp., ext) refers to the trace of the function on a side of 𝐸 coming from the interior of 𝐸 (resp., coming from the exterior of 𝐸 on that side). When the side of 𝐸 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 𝐳, but linear with respect to 𝐮,𝐯 and 𝐭.

Some important properties satisfied by the forms 𝑎,𝑏,𝑐 (cf. Wheeler [20], and Girault et al. [16, 17]) are recalled.

Lemma 3.1 (Coercivity). If 𝜖_𝑎=1, assume that 𝜎=1. If 𝜖_𝑎∈{−1,0}, assume that 𝜎 is sufficiently large. Then, there is a constant 𝜅>0, independent of ℎ, such that 𝑎(𝐯,𝐯)+𝐽₀(𝐯,𝐯)≥𝜅‖𝐯‖²_𝑋,∀𝐯∈𝐗^ℎ. (3.5)

It is clear that 𝜅=1 if 𝜖_𝑎=1. Otherwise, 𝜅 is a constant that depends on the polynomial degree of 𝐯 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 ℎ such that inf_{𝑞∈𝑄ℎ}sup_{𝐯∈𝐗ℎ}𝑏(𝐯,𝑞)‖𝐯‖_𝑋‖𝑞‖_0,Ω≥𝛽. (3.6)

Lemma 3.3 (Positivity). One has 𝑐^𝐯(𝐯,𝐳,𝐳)≥0,∀𝐯,𝐳∈{𝐭∈(𝐿²(Ω))²∶𝐭|_𝐸∈(𝐻²(𝐸))²∀𝐸∈ℰ_ℎ}. (3.7)

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). Given 𝑣∈𝐿²(Ω), for a given filtering radius 𝛿>0,𝐺^ℎ∶𝐿²(Ω)→𝐗^ℎ, 𝐯^ℎ=𝐺^ℎ𝐯 where 𝐯^ℎ is the unique solution in 𝐗^ℎ of 𝛿²(𝑎(𝐯^ℎ,𝜙)+𝐽₀(𝐯^ℎ,𝜙))+(𝐯^ℎ,𝜙)=(𝐯,𝜙)∀𝜙∈𝐗^ℎ. (3.8)

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 operators 𝐺^ℎ_𝑁 are 𝐺^ℎ_𝑁𝐯∶=^𝑁_𝑛=0(Π_ℎ−𝐺_ℎ)^𝑛𝐯, (3.9) where Π_ℎ∶𝐿²(Ω)→𝐗^ℎ is the 𝐿² projection.

For 𝐯∈𝐗^ℎ, the discrete deconvolution operator for 𝑁=0,1,2 is𝐺^ℎ₀𝐺𝐯=𝐯,_ℎ1𝐯=2𝐯−𝐯^ℎ,𝐺^ℎ₂𝐯=3𝐯−3𝐯^ℎ+𝐯^ℎℎ. (3.10)

𝐺_𝑁 was shown to be an 𝑂(𝛿^2𝑁+2) approximate inverse to the filter operator 𝐺 in Lemma 2.1 of Dunca and Epshteyn [26], recalled next.

Lemma 3.7. 𝐺_𝑁 is a bounded, self-adjoint positive operator. 𝐺_𝑁 is an 𝑂(𝛿^2𝑁+2) asymptotic inverse to the filter 𝐺. Specifically, for smooth 𝜙 and as 𝛿→0, 𝜙=𝐺_𝑁𝜙+(−1)^(𝑁+1)𝛿^2𝑁+2Δ^𝑁+1𝐺^(𝑁+1)𝜙. (3.11)

Some basic facts about discrete differential filters and deconvolution operators are presented next.

Lemma 3.8. For 𝐯∈𝐿²(Ω), one has the following bounds for the discretely filtered and approximately deconvolved 𝐯:‖𝐯^ℎ‖_0,Ω≤‖𝐯‖_0,Ω,‖𝐺^ℎ_𝑁𝐯^ℎ‖_0,Ω≤𝐶(𝑁)‖𝐯‖_0,Ω. (3.12)(3.13)

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 𝜙=𝐯^ℎ. Then, using coercivity result 𝛿²𝜅‖𝐯^ℎ‖²_𝑋+‖𝐯^ℎ‖²_0,Ω≤(𝐯,𝐯^ℎ)_Ω≤12‖𝐯‖²_0,Ω+12‖𝐯^ℎ‖²_0,Ω. (3.14) The term 𝛿²𝜅‖𝐯^ℎ‖²_𝑋 is positive, so it will be dropped, which yields 12‖𝐯^ℎ‖²_0,Ω≤12‖𝐯‖²_0,Ω. (3.15) 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 𝐺^ℎ_𝑁.

Lemma 3.9. For smooth 𝜙 the discrete approximate deconvolution operator satisfies ‖𝜙−𝐺^ℎ_𝑁𝜙^ℎ‖_0,Ω≤𝐶₁𝛿^2𝑁+2‖𝐺𝜙‖_2𝑁+2,Ω+𝐶₂(𝛿ℎ^𝑟+ℎ^𝑟+1)(^𝑁+1_𝑛=2|𝐺^𝑛𝜙|_𝑟+1,Ω). (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: ‖𝜙−𝐺^ℎ_𝑁𝜙^ℎ‖_0,Ω≤‖𝜙−𝐺_𝑁𝜙‖_0,Ω+‖𝐺_𝑁𝜙−𝐺^ℎ_𝑁𝜙‖_0,Ω+‖𝐺^ℎ_𝑁𝜙−𝐺^ℎ_𝑁𝜙^ℎ‖_0,Ω. (3.17) Lemma 3.7 gives ‖𝜙−𝐺_𝑁𝜙‖_0,Ω≤𝐶𝛿^2𝑁+2‖𝜙‖_2𝑁+2,Ω. (3.18) The standard discontinuous finite element bound for (3.8) is given by (cf. Rivière [13]) ‖𝜙−𝜙^ℎ‖_0,Ω≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)|𝜙|_𝑟+1,Ω. (3.19)
Lemma 3.8 gives for the third term in (3.17) that ‖𝐺^ℎ_𝑁𝜙−𝐺^ℎ_𝑁𝜙^ℎ‖≤𝐶‖𝜙−𝜙^ℎ‖. Then, (3.19) is applied.
Now, it is left to bound the second term from (3.17). First, note that for 𝑁=0, ‖𝐺₀𝜙^ℎ−𝐺^ℎ₀𝜙^ℎ‖_0,Ω=0. Based on Definition 3.6 of continuous and discrete deconvolution operators and their expansion (see (3.10)), 𝐺_𝑁 is a polynomial of degree 𝑁 in 𝐺 (and 𝐺^ℎ_𝑁 in 𝐺^ℎ as well). Thus, the second term in (3.17) can be written as‖𝐺_𝑁𝜙−𝐺^ℎ_𝑁𝜙‖_0,Ω=‖^𝑁_𝑛=0𝛼_𝑛(𝐺^𝑛𝜙−(𝐺^ℎ)^𝑛𝜙)‖_0,Ω≤^𝑁_𝑛=0𝛼_𝑛‖𝐺^𝑛𝜙−(𝐺^ℎ)^𝑛𝜙‖_0,Ω. (3.20)
For 𝑂(1) coefficients 𝛼_𝑛 and for 𝑁=1, the result (3.19) gives‖𝐺𝜙−(𝐺^ℎ)𝜙‖_0,Ω=‖𝜙−𝜙^ℎ‖_0,Ω≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)(|𝜙|_𝑟+1,Ω). (3.21)
For 𝑁=2, the results (3.19) and (3.12) give ‖𝐺₂𝜙−(𝐺^ℎ₂)𝜙‖_0,Ω=‖𝜙−𝜙^ℎℎℎ‖_0,Ω≤‖𝜙−𝜙^ℎ‖_0,Ω+‖𝜙^ℎ−𝜙^ℎℎℎ‖_0,Ω≤‖𝜙−𝜙^ℎ‖_0,Ω+‖𝜙−𝜙^ℎ‖_0,Ω≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)(|𝜙|_𝑟+1,Ω+|𝜙|_𝑟+1,Ω). (3.22)
Inductively, ‖𝐺_𝑁𝜙−𝐺^ℎ_𝑁𝜙‖_0,Ω≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)(^𝑁+1_𝑛=2|𝐺^𝑛𝜙|_𝑟+1,Ω). (3.23)
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, |(𝐺)^𝑛𝜙|_𝑟+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 𝑘|𝐺𝜙|_𝑟+1,Ω≤𝐶|𝜙|_𝑟+1,Ωandthus|𝐺^𝑛𝜙|_𝑟+1,Ω≤𝐶|𝜙|_𝑟+1,Ω. (3.24) The nonperiodic case can be more delicate. Suppose 𝜕Ω∈𝐶^𝑟+3 and 𝜙=0 on 𝜕Ω (i.e., 𝜙∈𝐻¹₀(Ω)∩𝐻^𝑟+1(Ω)). Call 𝐺𝜙=𝜙 so 𝜙 satisfies −𝛿²Δ𝜙+𝜙=𝜙inΩ,𝜙=0on𝜕Ω. (3.25) Then it is known that 𝜙∈𝐻^𝑟+3(Ω)∩𝐻¹₀(Ω), and Δ𝜙=0 on 𝜕Ω. Further, ‖𝜙‖_𝑗,Ω≤𝐶‖𝜙‖_𝑗,Ω,𝑗=0,1,2. (3.26) So, (3.24) holds for 𝑟=−1,0,+1. It also holds for higher values of 𝑟 provided additionally Δ^𝑗𝜙=0 on 𝜕Ω for 0≤𝑗≤[(𝑟+1)/2]−1.
Now consider the second term 𝑛=2, that is, 𝐺²𝜙=𝜙. We know from elliptic theory for 𝜙∈𝐻^𝑟+1(Ω)⋂𝐻¹₀(Ω), that 𝜙∈𝐻^𝑟+3(Ω)⋂𝐻¹₀(Ω), (as noted above) Δ𝜙=0 on 𝜕Ω and−𝛿²Δ𝜙+𝜙=𝜙inΩ,𝜙=Δ𝜙=0on𝜕Ω. (3.27) Theorem 1.1 in [28] then implies, uniformly in 𝛿, ‖𝜙‖_𝑗,Ω≤𝐶‖𝜙‖_𝑗,Ω,𝑗=0,1,2,3,4. (3.28) This extends directly to 𝐺^𝑛𝜙.

Extending Lemma 3.9, the following assumption will be made.

Assumption DG1. The |𝐺^𝑛(𝜙)|_𝑟+1,Ω terms in (3.16) are independent of 𝛿 and ‖𝜙−𝐺^ℎ_𝑁𝜙^ℎ‖≤𝐶₁𝛿^2𝑁+2‖𝜙‖_2𝑁+2,Ω+𝐶₂(𝛿ℎ^𝑟+ℎ^𝑟+1)|𝜙|_𝑟+1,Ω. (3.29)

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. 𝐺^ℎ and (𝐼−𝐺^ℎ_𝑁𝐺^ℎ) are symmetric, positive definite (𝑆𝑃𝐷) 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 𝜖_𝑎=−1 for the discretization of the filter problem (3.8), so that the bilinear form 𝑎(⋅,⋅) is symmetric.

The numerical scheme that uses discontinuous finite elements in space and backward Euler in time is derived next. For this, let Δ𝑡 denote the time step such that 𝑀=𝑇/Δ𝑡 is a positive integer. Let 0=𝑡₀<𝑡₁<⋯<𝑡_𝑀=𝑇 be a subdivision of the interval (0,𝑇). The function 𝜙 evaluated at the time 𝑡_𝑚 is denoted by 𝜙_𝑚. With the above forms, the fully discrete scheme is: find (𝐮^ℎ_𝑛,𝑝^ℎ_𝑛)_𝑛≥0∈𝐗^ℎ×𝑄^ℎ such that:1Δ𝑡(𝐮^ℎ_𝑛+1−𝐮^ℎ_𝑛,𝐯)_Ω+𝜈(𝑎(𝐮^ℎ_𝑛+1,𝐯)+𝐽₀(𝐮^ℎ_𝑛+1,𝐯))+𝜒(𝐮^ℎ_𝑛+1−𝐺^ℎ_𝑁𝐮^ℎℎ_𝑛+1,𝐯)_Ω+𝑐^𝐮ℎ𝑛(𝐮^ℎ_𝑛;𝐮^ℎ_𝑛+1,𝐯)+𝑏(𝐯,𝑝^ℎ_𝑛+1)=(𝐟_𝑛+1,𝐯)_Ω∀𝐯∈𝐗^ℎ,𝑏(𝐮^ℎ_𝑛+1,𝑞)=0∀𝑞∈𝑄^ℎ,(𝐮^ℎ₀,𝐯)_Ω=(𝐮₀,𝐯)_Ω∀𝐯∈𝐗^ℎ. (3.30)(3.31)(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 compute [30].

The consistency result of the semidiscrete scheme is showed next.

Lemma 3.12 (Consistency). Let (𝐮,𝑝) be the solution to (1.7)–(1.10), then (𝐮,𝑝) satisfies (𝐮_𝑡,𝐯)_Ω+𝜈(𝑎(𝐮,𝐯)+𝐽₀(𝐮,𝐯))+𝜒(𝐮−𝐺^ℎ_𝑁𝐮^ℎ,𝐯)_Ω+𝑐^𝐮(𝐮;𝐮,𝐯)+𝑏(𝐯,𝑝)=(𝐟,𝐯)_Ω+𝐸(𝐮,𝑝,𝐟;𝐯)∀𝐯∈𝐗^ℎ,𝑏(𝐮,𝑞)=0∀𝑞∈𝑄^ℎ,(𝐮(0),𝐯)_Ω=(𝐮₀,𝐯)_Ω∀𝐯∈𝐗^ℎ, (3.33)(3.34)(3.35) where the consistency error 𝐸(𝐮,𝐯)=𝜒(𝐮−𝐺^ℎ_𝑁𝐮^ℎ,𝐯)_Ω−𝜒(𝐮−𝐺_𝑁𝐮,𝐯)_Ω. Furthermore, 𝐸(𝐮,𝐯)≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)|𝐮|_𝑟+1,Ω‖𝐯‖_0,Ω. (3.36)

Proof. Equations (3.34) and (3.35) are clearly satisfied because of (1.8), (1.9), and (1.10) and the regularity of 𝐮. Next, we multiply (1.7) by 𝐯 and integrate over one mesh element 𝐸(𝐮_𝑡,𝐯)_𝐸+(∇⋅(𝐮𝐮),𝐯)_𝐸−𝜈(Δ𝐮,𝐯)_𝐸+𝜒(𝐮−𝐺_𝑁𝐮,𝐯)_𝐸+(∇𝑝,𝐯)_𝐸=(𝐟,𝐯)_𝐸. (3.37) Summing over all elements 𝐸_{𝐸∈ℰℎ}(𝐮_𝑡,𝐯)_𝐸+(∇⋅(𝐮𝐮),𝐯)_Ω−𝜈_{𝐸∈ℰℎ}(Δ𝐮,𝐯)_𝐸+𝜒(𝐮−𝐺_𝑁𝐮,𝐯)_Ω+(∇𝑝,𝐯)_Ω=(𝐟,𝐯)_Ω. (3.38)
By Green's formula−𝜈_{𝐸∈ℰℎ}(Δ𝐮,𝐯)_𝐸=−𝜈(Δ𝐮,𝐯)_Ω=𝜈_{𝐸∈ℰℎ}(∇𝐮,∇𝐯)_𝐸−𝜈_{𝑒∈Γℎ∪𝜕Ω𝑒}(∇𝐮)𝐧_𝑒⋅[𝐯]. (3.39) The regularity of 𝐮 then yields −𝜈_{𝐸∈ℰℎ}(Δ𝐮,𝐯)_𝐸=𝜈(𝑎(𝐮,𝐯)+𝐽₀(𝐮,𝐯)). (3.40) Note that Green's formula yields (∇𝑝,𝐯)_Ω=𝑏(𝐯,𝑝), (3.41) and that the incompressibility condition with the regularity of 𝐮 yields (𝐮⋅∇𝐮,𝐯)_Ω=𝑐^𝐮(𝐮;𝐮,𝐯). (3.42) The final result is obtained by bounding the consistency error 𝐸(𝐮,𝐯). For 𝑁=0, we have 𝐸(𝐮,𝐯)≤‖𝐮−𝐮^ℎ‖_0,Ω‖𝐯‖_0,Ω≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)|𝐮|_𝑟+1,Ω‖𝐯‖_0,Ω. (3.43) For 𝑁=1, 𝐸(𝐮,𝐯)≤‖𝐺𝐮−(𝐺^ℎ)𝐮^ℎ‖_0,Ω‖𝐯‖_0,Ω≤‖𝐮−𝐮^ℎℎ‖_0,Ω‖𝐯‖_0,Ω≤(‖𝐮−𝐮^ℎ‖_0,Ω+‖𝐮^ℎ−𝐮^ℎℎ‖_0,Ω)‖𝐯‖_0,Ω≤(‖𝐮−𝐮^ℎ‖_0,Ω+‖𝐮−𝐮^ℎ‖_0,Ω)‖𝐯‖_0,Ω≤𝐶(𝛿ℎ^𝑟+ℎ^𝑟+1)(|𝐮|_𝑟+1,Ω+|𝐮|_𝑟+1,Ω)‖𝐯‖_0,Ω. (3.44) The bound for 𝐸(𝐮,𝐯) 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 Lemma 3.1 holds. Then, there exists a unique solution to (3.30)–(3.32).

Proof. The existence of 𝐮^ℎ₀ is trivial. Given 𝐮^ℎ_𝑛, the problem of finding a unique 𝐮^ℎ_𝑛+1 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 𝐕^ℎ defined by 𝐕^ℎ={𝐯∈𝐗^ℎ∶𝑏(𝐯,𝑞)=0∀𝑞∈𝑄^ℎ}. (3.45) Let 𝐮^ℎ_𝑛+1 and 𝐮^ℎ_𝑛+1 be two solutions and let 𝜽_𝑛+1=𝐮^ℎ_𝑛+1−𝐮^ℎ_𝑛+1. Then, 𝜽_𝑛+1 satisfies: 1Δ𝑡(𝜽_𝑛+1,𝐯)_Ω+𝜈(𝑎(𝜽_𝑛+1,𝐯)+𝐽₀(𝜽_𝑛+1,𝐯))+𝑐^𝐮ℎ𝑛(𝐮^ℎ_𝑛;𝜽_𝑛+1,𝐯)+𝜒(𝜽_𝑛+1−𝐺^ℎ_𝑁𝜽^ℎ_𝑛+1,𝐯)=0∀𝐯∈𝐕^ℎ. (3.46) Choosing 𝐯=𝜽_𝑛+1 and using the coercivity result (3.5), positivity result (3.7) and positivity of the operator 𝐼−𝐺^ℎ_𝑁𝐺^ℎ given in Assumption DG2, we obtain: 1Δ𝑡‖𝜽_𝑛+1‖²_0,Ω+𝜈𝜅‖𝜽_𝑛+1‖²_𝑋≤0, (3.47) which yields that 𝜽_𝑛+1=𝟎. The existence and uniqueness of the pressure 𝑝^ℎ_𝑛+1 is then obtained from the inf-sup condition (3.6).

Some approximation properties of the spaces 𝐗^ℎ and 𝑄^ℎ are recalled next. From Crouzeix and Raviart [31], and Girault et al. [16], for each integer 𝑟≥1, and for any 𝐯∈(𝐻¹₀(Ω))², there is a unique discrete velocity 𝐯∈𝐗^ℎ such that 𝑏(𝐯−𝐯,𝑞)=0∀𝑞∈𝑄^ℎ. (3.48) Furthermore, if 𝐯∈(𝐻¹₀(Ω))²∩(𝐻^𝑟+1(Ω))², there is a constant 𝐶 independent of ℎ such that‖𝐯−𝐯‖_𝑋≤𝐶ℎ^𝑟|𝐯|_𝑟+1,Ω,|𝐯−𝐯|_𝑚,Ω≤𝐶ℎ^{𝑟+1−𝑚}|𝐯|_𝑟+1,Ω,𝑚=0,1. (3.49)(3.50)

For the pressure space, we use the approximation given by the 𝐿² projection. For any 𝑞∈𝐿²₀(Ω), there exists a unique discrete pressure 𝑞∈𝑄^ℎ such that(𝑞−𝑞,𝑧)_Ω=0∀𝑧∈𝑄^ℎ. (3.51) In addition, if 𝑞∈𝐻^𝑟(Ω), then‖𝑞−𝑞‖_𝑚,𝐸≤𝐶ℎ^𝑟−𝑚|𝑞|_𝑟,𝐸,∀𝐸∈ℰ_ℎ,𝑚=0,1,2. (3.52)

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 Δ𝑡, 𝐻, and 𝑎_𝑛, 𝑏_𝑛, 𝑐_𝑛, 𝛾_𝑛 (for integers 𝑛≥0) be nonnegative numbers such that 𝑎_𝑙+Δ𝑡^𝑙_𝑛=0𝑏_𝑛≤Δ𝑡^𝑙_𝑛=0𝛾_𝑛𝑎_𝑛+Δ𝑡^𝑙_𝑛=0𝑐_𝑛+𝐻for𝑙≥0. (3.53) Suppose that Δ𝑡𝛾_𝑛<1, for all 𝑛, and set 𝜎_𝑛=(1−Δ𝑡𝛾_𝑛)⁻¹. Then, 𝑎_𝑙+Δ𝑡^𝑙_𝑛=0𝑏_𝑛≤exp(Δ𝑡^𝑙_𝑛=0𝜎_𝑛𝛾_𝑛){Δ𝑡^𝑙_𝑛=0𝑐_𝑛+𝐻}for𝑙≥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 equations (𝐰,𝑃) that is approximated is a strong solution and in particular satisfies (cf. Rivière [13])(𝐰_𝑡,𝐯)_Ω+𝜈(𝑎(𝐰,𝐯)+𝐽₀(𝐰,𝐯))+𝑐^𝐰(𝐰;𝐰,𝐯)+𝑏(𝐯,𝑃)=(𝐟,𝐯)_Ω∀𝐯∈𝐗^ℎ,𝑏(𝐰,𝑞)=0∀𝑞∈𝑄^ℎ,(𝐰^ℎ₀,𝐯)_Ω=(𝐮₀,𝐯)_Ω∀𝐯∈𝐗^ℎ. (3.55)(3.56)(3.57)

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 that 𝐰∈𝑙²(0,𝑇;(𝐻^𝑟+1(Ω))²)∩𝑙²(0,𝑇;(𝐻^2𝑁+2(Ω))², 𝐰_𝑡∈𝑙²(0,𝑇;(𝐻^𝑟+1(Ω))²)∩𝐿^∞((0,𝑇)×Ω), 𝐰_𝑡𝑡∈𝐿²(0,𝑇;(𝐻¹(Ω))²), 𝑝∈𝑙²(0,𝑇;𝐻^𝑟(Ω)) and 𝐮₀∈(𝐻^𝑟+1(Ω))². Assume also that the coercivity Lemma 3.1 holds and that 𝐰 satisfies (3.24). If 𝛿 is chosen of the order of ℎ, and Δ𝑡<1, there exists a constant 𝐶, independent of ℎ and Δ𝑡 but dependent on 𝜈⁻¹ such that the following error bound holds, for any 1≤𝑚≤𝑀: ‖𝐰_𝑚−𝐮^ℎ_𝑚‖_0,Ω+(𝜈𝜅Δ𝑡^𝑚_𝑛=1‖𝐰_𝑛−𝐮^ℎ_𝑛‖²_𝑋)^1/2≤𝐶ℎ^𝑟(𝜈⁻¹+𝜈+𝜒+1)+𝐶Δ𝑡+𝐶𝜒ℎ^2𝑁+2. (4.1)

Remark 4.2. The dependence of these error estimates with respect to the Reynolds number Re (∼ 1/𝜈) is 𝒪(Re𝑒^Re), which is an improvement with respect to the continuous finite element method where the dependence is 𝒪(Re𝑒^Re3).

Proof. Defining 𝐞_𝑛=𝐮^ℎ(𝑡^𝑛)−𝐰(𝑡^𝑛) and subtracting (3.55) from (3.30), we have 1Δ𝑡(𝐞_𝑛+1−𝐞_𝑛,𝐯)_Ω−(𝐰_𝑡(𝑡_𝑛+1),𝐯)_Ω+𝜈(𝑎+𝐽₀)(𝐞_𝑛+1,𝐯)+𝜒(𝐞_𝑛+1−𝐺^ℎ_𝑁𝐞^ℎ_𝑛+1,𝐯)_Ω+𝑐^𝐮ℎ𝑛(𝐮^ℎ_𝑛;𝐮^ℎ_𝑛+1,𝐯)−𝑐^𝐰𝑛+1(𝐰_𝑛+1;𝐰_𝑛+1,𝐯)+𝑏(𝐯,𝑝^ℎ_𝑛+1−𝑃_𝑛+1)1=−Δ𝑡(𝐰_𝑛+1−𝐰_𝑛,𝐯)_Ω−𝜒(𝐰_𝑛+1−𝐺^ℎ_𝑁𝐰^ℎ_𝑛+1,𝐯)_Ω∀𝐯∈𝐗^ℎ,𝑏(𝐞_𝑛+1,𝑞)=0∀𝑞∈𝑄^ℎ. (4.2)(4.3) We now decompose the error 𝐞_𝑛=𝜙_𝑛−𝜼_𝑛, where 𝜙_𝑛=𝐮^ℎ_𝑛−𝐰_𝑛 and 𝜼_𝑛 is the interpolation error 𝜼_𝑛=𝐰_𝑛−𝐰_𝑛. Choosing 𝐯=𝜙_𝑛+1 in the equation above, using the coercivity result (3.5) and positivity of the operator 𝐼−𝐺^ℎ_𝑁𝐺^ℎ, we obtain for (4.2) 12Δ𝑡(‖𝜙_𝑛+1‖²_0,Ω−‖𝜙_𝑛‖²_0,Ω)+𝜈𝜅‖𝜙_𝑛+1‖²_𝑋+𝑐^𝐮ℎ𝑛(𝐮^ℎ_𝑛;𝐮^ℎ_𝑛+1,𝜙_𝑛+1)−𝑐^𝐰𝑛+1(𝐰_𝑛+1;𝐰_𝑛+1,𝜙_𝑛+1)≤(𝜼_𝑡(𝑡_𝑛+1),𝜙_𝑛+1)_Ω+𝜈(𝑎+𝐽₀)(𝜼_𝑛+1,𝜙_𝑛+1)+(𝐰_𝑡(𝑡_𝑛+11)−(𝐰Δ𝑡^𝑛+1−𝐰_𝑛),𝜙_𝑛+1)_Ω+𝜒(𝜼_𝑛+1−𝐺^ℎ_𝑁𝜼^ℎ_𝑛+1,𝜙_𝑛+1)_Ω−𝜒(𝐰_𝑛+1−𝐺^ℎ_𝑁𝐰^ℎ_𝑛+1,𝜙_𝑛+1)_Ω+𝑏(𝜙_𝑛+1,𝑃_𝑛+1−𝑝^ℎ_𝑛+1). (4.4) Consider now the nonlinear terms from the above equation. First note that since 𝐰 is continuous, we can rewrite 𝑐^𝐰𝑛+1(𝐰_𝑛+1;𝐰_𝑛+1,𝜙_𝑛+1)=𝑐^𝐮ℎ𝑛(𝐰_𝑛+1;𝐰_𝑛+1,𝜙_𝑛+1), (4.5) so, for readability, the superscript 𝐰^ℎ_𝑛 in the 𝑐 form is dropped. Therefore, adding and subtracting the interpolant 𝐰_𝑛+1 yields 𝑐^𝐮ℎ𝑛(𝐮^ℎ_𝑛,𝐮^ℎ_𝑛+1,𝜙_𝑛+1)−𝑐^𝐮ℎ𝑛(𝐰_𝑛+1,𝐰_𝑛+1,𝜙_𝑛+1)=𝑐(𝐮^ℎ_𝑛,𝜙_𝑛+1,𝜙_𝑛+1)−𝑐(𝜙_𝑛,𝜼_𝑛+1,𝜙_𝑛+1)+𝑐(𝜙_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)−𝑐(𝜼_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)−𝑐(𝐰_𝑛,𝜼_𝑛+1,𝜙_𝑛+1)−𝑐(𝐰_𝑛+1−𝐰_𝑛,𝐰_𝑛+1,𝜙_𝑛+1). (4.6) Thus, the error equation (4.4) is rewritten as 12Δ𝑡(‖𝜙_𝑛+1‖²_0,Ω−‖𝜙_𝑛‖²_0,Ω)+𝜈𝜅‖𝜙_𝑛+1‖²_𝑋+𝑐(𝐮^ℎ_𝑛,𝜙_𝑛+1,𝜙_𝑛+1)≤|𝑐(𝜙_𝑛,𝜼_𝑛+1,𝜙_𝑛+1)|+|𝑐(𝜙_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)|+|𝑐(𝜼_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)|+|𝑐(𝐰_𝑛,𝜼_𝑛+1,𝜙_𝑛+1)|+|𝑐(𝐰_𝑛+1−𝐰_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)|+|(𝜼_𝑡(𝑡_𝑛+1),𝜙_𝑛+1)_Ω|+|(𝐰_𝑡(𝑡_𝑛+11)−(𝐰Δ𝑡^𝑛+1−𝐰_𝑛),𝜙_𝑛+1)_Ω|+|𝜈(𝑎+𝐽₀)(𝜼_𝑛+1,𝜙_𝑛+1)|+|𝜒(𝜼_𝑛+1−𝐺^ℎ_𝑁𝜼^ℎ_𝑛+1,𝜙_𝑛+1)_Ω|+|𝜒(𝐰_𝑛+1−𝐺^ℎ_𝑁𝐰^ℎ_𝑛+1,𝜙_𝑛+1)_Ω|+|𝑏(𝜙_𝑛+1,𝑃_𝑛+1−𝑝^ℎ_𝑛+1)|≤|𝑇₀|+|𝑇₁|+⋯+|𝑇₁₀|. (4.7) From property (3.7), the term 𝑐(𝐰^ℎ_𝑛;𝜙_𝑛+1,𝜙_𝑛+1) in the left-hand side of (4.7) is positive and therefore it will be dropped. For the other terms of the form 𝑐(⋅,⋅,⋅) 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 𝐶 is a generic constant that is independent of ℎ,𝜈 and Δ𝑡, and that takes different values at different places |𝑇₀|=|𝑐(𝜙_𝑛,𝜼_𝑛+1,𝜙_𝑛+1)|≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈‖𝜙_𝑛‖²_0,Ω,|𝑇₁|=|𝑐(𝜙_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)|≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈‖𝜙_𝑛‖²_0,Ω,|𝑇₂|=|𝑐(𝜂_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)|≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛|²_2𝑟,Ω,|𝑇₃|=|𝑐(𝐰_𝑛,𝜼_𝑛+1,𝜙_𝑛+1)|≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛|²_2𝑟,Ω,|𝑇₄|=|𝑐(𝐰_𝑛+1−𝐰_𝑛,𝐰_𝑛+1,𝜙_𝑛+1)|≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈Δ𝑡²‖𝐰_𝑡‖²_{𝐿∞([𝑡𝑛,𝑡𝑛+1]×Ω)}. (4.8) Therefore, we have |𝑇₀|+⋯+|𝑇₄|≤5𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈⁻¹‖𝜙_𝑛‖²_0,Ω+𝐶𝜈⁻¹ℎ^2𝑟|𝐰_𝑛|²_𝑟+1,Ω+𝐶𝜈⁻¹Δ𝑡²‖𝐰_𝑡‖²_{𝐿∞([𝑡𝑛,𝑡𝑛+1]×Ω)}. (4.9)
To bound 𝑇₅, Cauchy-Schwarz's inequality, Young's inequality and the approximation result (3.50) for 𝐰_𝑡 are applied|𝑇₅|≤‖𝜙_𝑛+1‖_0,Ω‖𝜼_𝑡(𝑡_𝑛+1)‖_0,Ω≤18‖𝜙_𝑛+1‖²_0,Ω+𝐶ℎ^2𝑟+2|𝐰_𝑡(𝑡^𝑛+1)|²_𝑘+1,Ω. (4.10) To bound the term 𝑇₆, a Taylor expansion with integral remainder is used 𝐰_𝑛=𝐰_𝑛+1−Δ𝑡𝐰_𝑡(𝑡_𝑛+11)+2^𝑡𝑛+1_𝑡𝑛(𝑠−𝑡_𝑛)𝐰_𝑡𝑡(𝑠)𝑑𝑠. (4.11) This implies that ‖𝐰_𝑡(𝑡_𝑛+1𝐰)−^𝑛+1−𝐰_𝑛‖Δ𝑡²_0,Ω≤Δ𝑡6^𝑡𝑛+1_𝑡𝑛‖𝐰_𝑡𝑡(𝑠)‖²_0,Ω𝑑𝑠. (4.12) Thus, with (3.50), we have |𝑇₆1|≤8‖𝜙_𝑛+1‖²_0,Ω+𝐶Δ𝑡^𝑡𝑛+1_𝑡𝑛‖𝐰_𝑡𝑡(𝑠)‖²_0,Ω≤1𝑑𝑠8‖𝜙_𝑛+1‖²_0,Ω+𝐶Δ𝑡^𝑡𝑛+1_𝑡𝑛‖𝐰_𝑡𝑡(𝑠)‖²_0,Ω𝑑𝑠. (4.13) Next, the term 𝑇₇ is expanded as |𝑇₇|≤|𝜈_{𝐸∈ℰℎ𝐸}∇𝜼_𝑛+1∶∇𝜙_𝑛+1|+|𝜈_{𝑒∈Γℎ∪𝜕Ω𝑒}{∇𝜼_𝑛+1}𝐧_𝑒⋅[𝜙_𝑛+1]|+|𝜈𝜖_{𝑎𝑒∈Γℎ∪𝜕Ω𝑒}{∇𝜙_𝑛+1}𝐧_𝑒⋅[𝜼_𝑛+1]|+|𝜈𝐽₀(𝜼_𝑛+1,𝜙_𝑛+1)|=|𝑇₇₁|+|𝑇₇₂|+|𝑇₇₃|+|𝑇₇₄|. (4.14) The term 𝑇₇₁ is bounded using Cauchy-Schwarz inequality, Young's inequality and the approximation result (3.49) |𝑇₇₁|≤𝜈‖𝜙_𝑛+1‖_𝑋‖𝜼_𝑛+1‖_𝑋≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈‖𝜼_𝑛+1‖²_𝑋≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛+1|²_𝑟+1,Ω. (4.15) Using Cauchy-Schwarz's inequality, trace inequality (2.8) and approximation result (3.50), we have |𝑇₇₂|≤𝜈_{𝑒∈Γℎ∪𝜕Ω}‖{∇𝜼_𝑛+1}𝐧_𝑒‖_{0,𝑒𝑒∈Γℎ∪𝜕Ω}‖[𝜙_𝑛+1]‖_0,𝑒≤𝐶𝜈(_{𝑒∈Γℎ∪𝜕Ω}1|𝑒|‖[𝜙_𝑛+1]‖²_0,𝑒)^1/2(|‖∇𝜼_𝑛+1‖|_0,Ω+ℎ|‖∇²𝜼_𝑛+1‖|_0,Ω)≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛+1|²_𝑟+1,Ω. (4.16) Using Cauchy-Schwarz's inequality, trace inequality (2.10), and approximation result (3.49), we have |𝑇₇₃|≤𝜈(_{𝑒∈Γℎ∪𝜕Ω}‖{∇𝜙_𝑛+1}𝐧_𝑒‖²_0,𝑒)^1/2(_{𝑒∈Γℎ∪𝜕Ω}‖[𝜼_𝑛+1]‖²_0,𝑒)^1/2≤𝐶𝜈‖𝜙_𝑛+1‖_𝑋(_{𝑒∈Γℎ∪𝜕Ω}1|𝑒|‖[𝜼_𝑛+1]‖²_0,𝑒)^1/2≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛+1|²_𝑟+1,Ω. (4.17) Using the approximation result (3.49), we have |𝑇₇₄|≤𝜈(_{𝑒∈Γℎ∪𝜕Ω}𝜎|𝑒|‖[𝜙_𝑛+1]‖²_0,𝑒)^1/2(_{𝑒∈Γℎ∪𝜕Ω}𝜎|𝑒|‖[𝜼_𝑛+1]‖²_0,𝑒)^1/2≤𝐶𝜈‖𝜙_𝑛+1‖_𝑋‖𝜼_𝑛+1‖_𝑋≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛+1|²_𝑟+1,Ω. (4.18) Putting together the bounds (4.15), (4.16), (4.17) and (4.18), |𝑇₇|≤4𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈ℎ^2𝑟|𝐰_𝑛+1|²_𝑟+1,Ω. (4.19) For the terms 𝑇₈ and 𝑇₉, Cauchy-Schwarz's inequality, Young's inequality and bounds (3.13) and (3.29) are applied |𝑇₈|=|𝜒(𝜼_𝑛+1−𝐺^ℎ_𝑁𝜼^ℎ_𝑛+1,𝜙_𝑛+1)_Ω|≤𝜒‖𝜼_𝑛+1−𝐺^ℎ_𝑁𝜼^ℎ_𝑛+1‖_0,Ω‖𝜙_𝑛+1‖_0,Ω≤𝐶𝜒²‖𝜼_𝑛+1−𝐺^ℎ_𝑁𝜼^ℎ_𝑛+1‖²_0,Ω+18‖𝜙_𝑛+1‖²_0,Ω≤𝐶(𝑁)𝜒²‖𝜼_𝑛+1‖²_0,Ω+18‖𝜙_𝑛+1‖²_0,Ω≤𝐶(𝑁)𝜒²ℎ^2𝑟+2|𝐰_𝑛+1|²_𝑟+1,Ω+18‖𝜙_𝑛+1‖²_0,Ω,|𝑇₉|=|𝜒(𝐰_𝑛+1−𝐺^ℎ_𝑁𝐰^ℎ_𝑛+1,𝜙_𝑛+1)_Ω|≤𝜒‖𝐰_𝑛+1−𝐺^ℎ_𝑁𝐰^ℎ_𝑛+1‖_0,Ω‖𝜙_𝑛+1‖_0,Ω≤𝐶𝜒²‖𝐰_𝑛+1−𝐺^ℎ_𝑁𝐰^ℎ_𝑛+1‖²_0,Ω+18‖𝜙_𝑛+1‖²_0,Ω≤𝐶𝜒²𝛿^4𝑁+4‖𝐰_𝑛+1‖²_2𝑁+2,Ω+𝐶𝜒²(𝛿²ℎ^2𝑟+ℎ^2𝑟+2)|𝐰_𝑛+1|²_𝑟+1,Ω+18‖𝜙_𝑛+1‖²_0,Ω. (4.20)(1)(2)(3)(4)(5)(6)(7)(8) Using (4.3) with (3.48) and (3.51), the pressure term 𝑇₁₀ is reduced to |𝑇₁₀|=|𝑏(𝜙_𝑛+1,𝑝_𝑛+1−𝑝_𝑛+1)+𝑏(𝜙_𝑛+1,𝑝_𝑛+1−𝑝^ℎ_𝑛+1)|=|𝑏(𝜙_𝑛+1,𝑝_𝑛+1−𝑝_𝑛+1)|=|_{𝑒∈Γℎ𝑒}{𝑝_𝑛+1−𝑝_𝑛+1}[𝜙_𝑛+1]⋅𝐧_𝑒|, (4.21) which is bounded by using Cauchy-Schwarz's inequality, Young's inequality, trace inequality (2.7) and approximation result (3.52) |𝑇₁₀|≤𝐶(_{𝑒∈Γℎ∪𝜕Ω}1|𝑒|‖[𝜙_𝑛+1]‖_0,𝑒)^1/2(‖𝑝_𝑛+1−𝑝_𝑛+1‖_0,Ω+ℎ|‖∇𝑝_𝑛+1−∇𝑝_𝑛+1‖|_0,Ω)≤𝜈𝜅20‖𝜙_𝑛+1‖²_𝑋+𝐶𝜈⁻¹ℎ^2𝑟|𝑝_𝑛+1|²_𝑟,Ω. (4.22) With the bounds (4.9), (4.10), (4.13), (4.19), (4.20), and (4.22), the error equation becomes 12Δ𝑡(‖𝜙_𝑛+1‖²_0,Ω−‖𝜙_𝑛‖²_0,Ω)+𝜈𝜅2‖𝜙_𝑛+1‖²_𝑋≤+𝐶𝜈⁻¹‖𝜙_𝑛‖²_0,Ω+12‖𝜙_𝑛+1‖²_0,Ω+𝐶𝜈⁻¹Δ𝑡²‖𝐰_𝑡‖²_{𝐿∞([𝑡𝑛,𝑡𝑛+1]×Ω)}+𝐶ℎ^2𝑟+2|𝐰_𝑡(𝑡^𝑛+1)|²_𝑘+1,Ω+𝐶Δ𝑡^𝑡𝑛+1_𝑡𝑛‖𝐰_𝑡𝑡(𝑠)‖²_0,Ω𝑑𝑠+𝐶𝜈⁻¹ℎ^2𝑟|𝐰_𝑛|²_𝑟+1,Ω+𝐶𝜈ℎ^2𝑟|𝐰_𝑛+1|²_𝑟+1,Ω+𝐶𝜒²ℎ^2𝑟+2‖𝐰_𝑛+1‖²_𝑟+1,Ω+𝐶(𝜒²(𝛿²ℎ^2𝑟+ℎ^2𝑟+2))|𝐰_𝑛+1|²_𝑟+1,Ω+𝐶𝜒²𝛿^4𝑁+4‖𝐰_𝑛+1‖²_2𝑁+2,Ω+𝐶𝜈⁻¹ℎ^2𝑟|𝑝_𝑛+1|²_𝑟,Ω, (4.23) where 𝐶 are constants independent of ℎ,𝜈 and Δ𝑡. Now, multiply the equation by 2Δ𝑡, sum from 𝑛=0 to 𝑛=𝑚−1, and use the assumption that 𝛿 is chosen of the order of ℎ to obtain (1−Δ𝑡)‖𝜙_𝑚‖²_0,Ω+𝜈𝜅Δ𝑡^𝑚−1_𝑛=0‖𝜙_𝑛+1‖²_𝑋≤‖𝜙₀‖²_0,Ω+𝐶(𝜈⁻¹+1)^𝑚−1_𝑛=0‖𝜙_𝑛‖²_0,Ω+𝐶ℎ^2𝑟Δ𝑡^𝑚−1_𝑛=0(𝜈⁻¹|𝐰_𝑛|²_𝑟+1,Ω+(𝜈+𝜒²)|𝐰_𝑛+1|²_𝑟+1,Ω+|𝐰_𝑡(𝑡_𝑛+1)|²_𝑟+1,Ω+𝜈⁻¹|𝑝_𝑛+1|²_𝑟,Ω)+𝐶Δ𝑡^2𝑇₀‖𝐰_𝑡𝑡(𝑠)‖²_1,Ω𝑑𝑠+𝐶Δ𝑡²𝜈⁻¹‖𝐰_𝑡‖²_{𝐿∞([0,𝑇]×Ω)}+𝐶Δ𝑡𝜒²ℎ^4𝑁+4‖𝐰_𝑛+1‖²_2𝑁+2,Ω. (4.24) Thus, using Gronwall's lemma with Δ𝑡<1, there is a constant 𝐶 independent of ℎ and Δ𝑡, but dependent on 𝜈⁻¹, such that ‖𝜙_𝑚‖²_0,Ω+𝜈𝜅Δ𝑡^𝑚_𝑛=1‖𝜙_𝑛‖²_𝑋≤‖𝜙₀‖²_0,Ω+𝐶ℎ^2𝑟(𝜈⁻¹+𝜈+𝜒²+1)+𝐶Δ𝑡²+𝐶𝜒²ℎ^4𝑁+4. (4.25) The final result is then obtained by noting that the term ‖𝜙₀‖²_0,Ω is of order ℎ^2𝑟 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 𝒪(Re𝑒^Re), whereas the dependence is 𝒪(Re𝑒^Re3) for continuous finite elements (cf. [14]).

The proposed scheme (3.30)–(3.32) contains parameter 𝜖_𝑎∈{−1,0,1} in the bilinear form 𝑎(⋅,⋅) 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 𝜖_𝑎 are preferred for a given mesh size. This is the object of a future paper.