# Convergence Analysis of a Fully Discrete Family of Iterated Deconvolution Methods for Turbulence Modeling with Time Relaxation

Authors: R. Ingram; C. C. Manica; N. Mays; I. Stanculescu
Journal: Advances in Numerical Analysis
 (2012)

---

## Abstract

We present a general theory for regularization models of the Navier-Stokes equations based on the Leray deconvolution model with a general deconvolution operator designed to fit a few important key properties. We provide examples of this type of operator, such as the (modified) Tikhonov-Lavrentiev and (modified) Iterated Tikhonov-Lavrentiev operators, and study their mathematical properties. An existence theory is derived for the family of models and a rigorous convergence theory is derived for the resulting algorithms. Our theoretical results are supported by numerical testing with the Taylor-Green vortex problem, presented for the special operator cases mentioned above.


---

## Body

## 1. Approximate Deconvolution for Turbulence Modeling

Numerical simulations of complex flows present many challenges. The Navier-Stokes (NS) equations (NSE), given by the following(1.1)ut+u·∇u-νΔu+∇p=f,∇·u=0,inΩ×(0,T)
is an exact model for the flow of a viscous, incompressible fluid, [1]. For turbulent flows (characterized by Reynold’s number Re≫1), it is infeasible to properly resolve all significant scales above the Kolmogorov length scale O(Re-3/4) by direct numerical simulation. Thus, numerical simulations are often based on various regularizations of NSE, rather than NSE themselves. Accordingly, regularization methods provide a (computationally) efficient and (algorithmically) simple family of turbulence models. Several of the most commonly applied regularization methods include:
(1.2)(Leray)wt+w¯·∇w-Re-1Δw+∇p=f,∇·w=0,(1.3)(NS-α)wt+w¯×(∇×w)-Re-1Δw+∇P=f,(1.4)(NS-ω¯)wt+w×(∇×w¯)-Re-1Δw+∇P=f,(1.5)(timerelaxation)wt+w·∇w-Re-1Δw+∇p+χ(w-w¯)=f,
where ∇·w=0 in each case and w¯ is an averaged velocity field w, p is pressure, and P is the Bernoulli pressure. More details about these models can be found for instance in [2–5] and references therein. Although these regularization methods achieve high theoretical accuracy and perform well in select practical tests, those models do not provide a fully-developed numerical solution for decoupling the scales in a turbulent flow. In fact, results show that only time relaxation regularization truncates scales sufficiently for practical computations. Indeed, it is shown that time relaxation term χ(w-w¯) for χ>0 damps unresolved fluctuations over time [5, 6]. Note that the choice of χ is an active area of research and that solutions are very sensitive to variations in χ.Deconvolution-based regularization is also an active area of research obtained, for example, by replacingw¯ by D(w¯) in each (1.2)–(1.5) for some deconvolution operator D. In [7], Dunca proposed the general Leray-deconvolution problem (D(w¯) instead of w¯) as a more accurate extension to Leray’s model [8]. Leray used the Gaussian filter as the smoothing (averaging) filter G, denoted above by overbar. In [9], Germano proposed the differential filter (approximate-Gaussian) G=(-δ2Δ+I)-1 where δ>0 is the filter length. The differential filter is easily modeled in the variational framework of the finite element (FE method (FEM)). We provide a brief overview of continuous and discrete operators (Section 3).The deconvolution-based models have proven themselves to be very promising. However, among the very many known approximate deconvolution operators from image processing, for example, [10], so far only few have been studied for turbulence modeling, for example, the van Cittert and the modified Tikhonov-Lavrentiev deconvolution operators. Their success suggests that it is time to develop a general theory for regularization models of the NSE as a guide to development of models based on other, possibly better, deconvolution operators and refinement of existing ones.Herein, we present a general theory for regularization models of the NSE based on the Leray deconvolution model with a general deconvolution operator. We prove energetic stability (and hence existence) and convergence of an FE (in space) and Crank-Nicolson (CN) (in time) discretization of the following family of Leray deconvolution regularization models with time relaxation: findw:Ω×(0,T]→ℝ3 and π:Ω×(0,T]→ℝ satisfying the following
(1.6)wt+D(w¯)·∇w-Re-1Δw+∇π+χ(w-D(w¯))=f,∇·w=0,
for some appropriate boundary and initial conditions (the fully discrete model is presented in Problem 1 with energetic stability and well posedness proved in Theorem 4.5).


### 1.1. Improving Accuracy of Approximate Deconvolution Methods

The fundamental difficulties corresponding to regularization methods applied as a viable turbulence model include ensuring that:(i)
scales are appropriately truncated (model microscale = filter radius = mesh width)(ii)
smooth parts of the solution are accurately approximated (D(w¯)≈w for smooth w)(iii)
physical fidelity of flow is preserved.Due to the nonlinearity in (1.6), different choices of the filter and deconvolution operator yield significant changes in the solution of the corresponding model. Implementation concerns for deconvolution methods, for example, Tikhonov-Lavrentiev regularization given by D=(G+αI)-1, include selection of deconvolution parameter α>0. Iterated deconvolution methods reduce approximation sensitivity relative to α-selection and, hence, allow a conservatively large α-selection for stability with updates (fixed number of iterations) used to recover higher accuracy (Section 3.3). For example we prove, under usual conditions, (Proposition 5.3),
(1.7)Modified  Tikhonov(j=0)error(w-D(w¯))≤𝒪(αδ2),Iterated  Modified  Tikhonov(j>0)error(w-D(w¯))≤𝒪((αδ2)j+1),
so that iterated modified Tikhonov regularization gives geometric convergence with respect to the update number j. In either case, Tikhonov-Lavrentiev or iterated Tikhonov-Lavrentiev regularization, we prove, under usual conditions, (Proposition 5.4)
(1.8)error(Dh(w¯h)-D(w¯))≤𝒪(hk),
where Dh and w¯h represent the discrete deconvolution operator and discrete filter, respectively, and k is order of FE polynomial space. We propose minimal properties for a general family of deconvolution operators D and filters G (Section 3.1) satisfying (Assumptions 3.4, 3.5), for example,(i)
||DG||≤1, forces spectrum of DG in [0,1](ii)
||∇DGv||≤d1||∇v||, controls size of ∇DG(iii)
||(I-DG)v||≤c1(α,δ)||Δv||→0 as α, δ→0 (for smooth v), ensures convergence of method(iv)
||(DG-DhGh)v||≤c2(h,α,δ)→0 as h→0 (for smooth v), ensures convergence of method.In fact, the updates{ωj}j≥0 satisfying
(1.9)ω0:=D(w¯),ωj-ωj-1:=D(w-ωj-1)¯,
inherit the properties assumed for the base operator D in Assumptions 3.4, 3.5 (Propositions 3.10, 3.11). These iterates represent defect correction generalization of iterated Tikhonov regularization operator [11]. We prove that the FE-CN approximation wh of the general deconvolution turbulence model (Problem 1) satisfies (Theorem 4.6)
(1.10)error(wh-uNSE)≤C(hk+Δt2+c1(α,δ,j)+c2(h,α,δ,j))→0,ash,δ,α→0
for smooth enough solutions uNSE of the NSE (see variational formulation (2.3)–(2.5)). We show that c1=C(αδ2)β (for any 0≤β≤j+1 and ΔβuNSE∈L2(Ω)) (Proposition 5.3) and c2=Chk (for D(w¯)∈Hk+1(Ω)) (Proposition 5.4) for modified iterated Tikhonov-Lavrentiev regularization (see Corollary 5.5 for corresponding error estimate). We conclude with a numerical test that verifies the theoretical convergence rate predicted in Theorem 4.6 (Section 5.2).


### 1.2. Background and Overview

One of the most interesting approaches to generate turbulence models is via approximate deconvolution or approximate/asymptotic inverse of the filtering operator. Examples of such models include: Approximate Deconvolution Models (ADM) and Leray-Tikhonov Deconvolution Models. Layton and Rebholz compiled a comprehensive overview and detailed analysis of ADM [12] (see also references therein). Previous analysis of the ADM with and without the time-relaxation term used van Cittert deconvolution operators [5, 6]; although easily programmed, van Cittert schemes can be computationally expensive [5]. Tikhonov-Lavrentiev regularization is another popular regularization scheme [13]. Determining the appropriate value of α to ensure stability while preserving accuracy is challenging, see for example, [14–19]. Alternatively, iterated Tikhonov regularization is well known to decouple stability and accuracy from the selection of regularization parameter α, see for example, [11, 20–22]. Iterated Tikhonov regularization is one special case of the general deconvolution operator we propose herein.


## 1.1. Improving Accuracy of Approximate Deconvolution Methods

The fundamental difficulties corresponding to regularization methods applied as a viable turbulence model include ensuring that:(i)
scales are appropriately truncated (model microscale = filter radius = mesh width)(ii)
smooth parts of the solution are accurately approximated (D(w¯)≈w for smooth w)(iii)
physical fidelity of flow is preserved.Due to the nonlinearity in (1.6), different choices of the filter and deconvolution operator yield significant changes in the solution of the corresponding model. Implementation concerns for deconvolution methods, for example, Tikhonov-Lavrentiev regularization given by D=(G+αI)-1, include selection of deconvolution parameter α>0. Iterated deconvolution methods reduce approximation sensitivity relative to α-selection and, hence, allow a conservatively large α-selection for stability with updates (fixed number of iterations) used to recover higher accuracy (Section 3.3). For example we prove, under usual conditions, (Proposition 5.3),
(1.7)Modified  Tikhonov(j=0)error(w-D(w¯))≤𝒪(αδ2),Iterated  Modified  Tikhonov(j>0)error(w-D(w¯))≤𝒪((αδ2)j+1),
so that iterated modified Tikhonov regularization gives geometric convergence with respect to the update number j. In either case, Tikhonov-Lavrentiev or iterated Tikhonov-Lavrentiev regularization, we prove, under usual conditions, (Proposition 5.4)
(1.8)error(Dh(w¯h)-D(w¯))≤𝒪(hk),
where Dh and w¯h represent the discrete deconvolution operator and discrete filter, respectively, and k is order of FE polynomial space. We propose minimal properties for a general family of deconvolution operators D and filters G (Section 3.1) satisfying (Assumptions 3.4, 3.5), for example,(i)
||DG||≤1, forces spectrum of DG in [0,1](ii)
||∇DGv||≤d1||∇v||, controls size of ∇DG(iii)
||(I-DG)v||≤c1(α,δ)||Δv||→0 as α, δ→0 (for smooth v), ensures convergence of method(iv)
||(DG-DhGh)v||≤c2(h,α,δ)→0 as h→0 (for smooth v), ensures convergence of method.In fact, the updates{ωj}j≥0 satisfying
(1.9)ω0:=D(w¯),ωj-ωj-1:=D(w-ωj-1)¯,
inherit the properties assumed for the base operator D in Assumptions 3.4, 3.5 (Propositions 3.10, 3.11). These iterates represent defect correction generalization of iterated Tikhonov regularization operator [11]. We prove that the FE-CN approximation wh of the general deconvolution turbulence model (Problem 1) satisfies (Theorem 4.6)
(1.10)error(wh-uNSE)≤C(hk+Δt2+c1(α,δ,j)+c2(h,α,δ,j))→0,ash,δ,α→0
for smooth enough solutions uNSE of the NSE (see variational formulation (2.3)–(2.5)). We show that c1=C(αδ2)β (for any 0≤β≤j+1 and ΔβuNSE∈L2(Ω)) (Proposition 5.3) and c2=Chk (for D(w¯)∈Hk+1(Ω)) (Proposition 5.4) for modified iterated Tikhonov-Lavrentiev regularization (see Corollary 5.5 for corresponding error estimate). We conclude with a numerical test that verifies the theoretical convergence rate predicted in Theorem 4.6 (Section 5.2).


## 1.2. Background and Overview

One of the most interesting approaches to generate turbulence models is via approximate deconvolution or approximate/asymptotic inverse of the filtering operator. Examples of such models include: Approximate Deconvolution Models (ADM) and Leray-Tikhonov Deconvolution Models. Layton and Rebholz compiled a comprehensive overview and detailed analysis of ADM [12] (see also references therein). Previous analysis of the ADM with and without the time-relaxation term used van Cittert deconvolution operators [5, 6]; although easily programmed, van Cittert schemes can be computationally expensive [5]. Tikhonov-Lavrentiev regularization is another popular regularization scheme [13]. Determining the appropriate value of α to ensure stability while preserving accuracy is challenging, see for example, [14–19]. Alternatively, iterated Tikhonov regularization is well known to decouple stability and accuracy from the selection of regularization parameter α, see for example, [11, 20–22]. Iterated Tikhonov regularization is one special case of the general deconvolution operator we propose herein.


## 2. Function Spaces and Approximations

Let the flow domainΩ∈ℝd for d=2,3 be a regular and bounded polyhedral. We use standard notation for Lebesgue and Sobolev spaces and their norms. Let ||·|| and (·,·) be the L2-norm and inner product, respectively. Let ||·||p,k:=||·||Wpk(Ω) represent the Wpk(Ω)-norm. We write Hk(Ω):=W2k(Ω) and ||·||k for the corresponding norm. Let the context determine whether Wpk(Ω) denotes a scalar, vector, or tensor function space. For example let v:Ω→ℝd. Then, v∈H1(Ω) implies that v∈H1(Ω)d and ∇v∈H1(Ω) implies that ∇v∈H1(Ω)d×d. Write Wqm(Wpk(Ω)):=Wqm(0,T;Wpk(Ω)) equipped with the standard norm. For example,
(2.1)∥v∥Lq(Wpk):={(∫0T∥v(·,t)∥p,kqdt)1/q,if1≤q<∞,esssup0<t<T∥v(·,t)∥p,k,ifq=∞.Denote the pressure and velocity spaces byQ:=L02(Ω)={q∈L2(Ω):∫Ω‍q=0} and X:=H01(Ω)={v∈H1(Ω):v|∂Ω=0}, respectively. Moreover, the dual space of X is denoted X′:=W2-1(Ω) and equipped with the norm
(2.2)∥f∥-1:=sup0≠v∈X〈f,v〉X′×X|v|1.Fixf∈X′ and Re>0. In this setting, we consider strong NS solutions: find u∈L2(X)∩L∞(L2(Ω)) and p∈W-1,∞(Q) satisfying
(2.3)(ut,v)+(u·∇u,v)+Re-1(∇u,∇v)-(p,∇·v)=(f,v),a.e.t∈(0,T],∀v∈X,(2.4)(q,∇·u)=0,a.e.t∈(0,T],∀q∈Q,(2.5)u(·,0)=u0inL2(Ω).
Let V:={v∈X:∇·v=0}. Restricting test functions v∈V reduces (2.3)–(2.5) to find u:(0,T]→V satisfying
(2.6)(ut,v)+(u·∇u,v)+Re-1(∇u,∇v)=(f,v),a.e.t∈(0,T],∀v∈V
and (2.5). For smooth enough solutions, solving the problem associated with (2.6), (2.5) is equivalent to (2.3)–(2.5).Control of the nonlinear term is essential for establishing a priori estimates and convergence estimates. We state a selection of inequalities here that will be utilized later:(2.7)|u·∇v,w|≤C(Ω){∥u∥1|v|1∥w∥1∥u∥∥v∥2∥w∥1∀v∈H2(Ω)∥u∥2∥v∥∥w∥1∀u∈H2(Ω)∩V.


### 2.1. Discrete Function Setting

Fixh>0. Let 𝒯h be a family of subdivisions (e.g., triangulation) of Ω¯⊂ℝd satisfying Ω¯=⋃E∈𝒯h‍E so that diameter(E)≤h and any two closed elements E1, E2∈𝒯h are either disjoint or share exactly one face, side, or vertex. See Chapter II, Appendix A in [23] for more on this subject in context of Stokes problem and [24] for a more general treatment. For example, 𝒯h consists of triangles for d=2 or tetrahedra for d=3 that are nondegenerate as h→0.LetXh⊂X and Qh⊂Q be a conforming velocity-pressure mixed FE space. For example, let Xh and Qh be continuous, piecewise (on each E∈𝒯h) polynomial spaces. The discretely divergence-free space is given by
(2.8)Vh={vh∈Xh:(qh,∇·vh)=0∀qh∈Qh}.
Note that in general Vh⊄V (e.g., for Taylor-Hood elements). In order to avoid stability issues arising when FE solutions are not exactly divergence free (i.e., when Vh⊄V), we introduce the explicitly skew-symmetric convective term
(2.9)bh(u,v,w):=12((u·∇v,w)-(u·∇w,v)),
so that
(2.10)bh(u,v,v)=0.
Note that bh(u,v,w)=(u·∇v,w) when u∈V. Moreover, the trilinear from bh(·,·,·) is continuous and skew-symmetric on X×X×X.Lemma 2.1.
Ifu,v,w∈X,
(2.11)|bh(u,v,w)|≤C(Ω){(∥u∥∥u∥1)1/2∥v∥1∥w∥1∥u∥∥v∥2∥w∥1∀v∈H2(Ω).Proof.
The proof of the first inequality can be found in [25]. The second follows from Hölder’s and Poincaré’s inequalities.For the time discretization, let0=t0<t1<⋯<tM-1=T<∞ be a discretization of the time interval [0,T] for a constant time step Δt=tn+1-tn. Write tn+1/2:=(tn+1+tn)/2,zn=z(tn) and, if z∈C0([tn,tn+1]), zn+1/2=(zn+1+zn)/2. Define
(2.12)∥u∥lq(m1,m2;Wpk(Ω)):={(Δt∑n=m1m2∥un∥k,pq)1/q,q∈[1,∞),maxm1≤n≤m2∥un∥k,p,q=∞,
for any 0≤n=m1,m1+1,…,m2≤M. Write ||u||lq(Wpk(Ω))=||u||lq(0,M;Wpk(Ω)). We say that u∈lq(m1,m2;Wpk(Ω)) if the associated norm defined above stays finite as Δt→0.The discrete Gronwall inequality is essential to the convergence analysis in Section4.2.Lemma 2.2.
LetD≥0 and κn, An, Bn, Cn≥0 for any integer n≥0 and satisfy
(2.13)AM+Δt∑n=0M‍Bn≤Δt∑n=0M‍κnAn+Δt∑n=0M‍Cn+D,∀M≥0.
Suppose that for all n, Δtκn<1 and set gn=(1-Δtκn)-1. Then,
(2.14)AM+Δt∑n=0M‍Bn≤exp(Δt∑n=0M‍gnκn)[Δt∑n=0M‍Cn+D],∀M≥0.Proof.
The proof follows from [26].


### 2.2. Approximation Theory

LetC>0 be a generic constant independent of h→0+. Preserving an abstract framework for the FE spaces, we assume that Xh×Qh inherit several fundamental approximation properties.Assumption 2.3.
The FE spacesXh×Qh satisfy:Uniform inf-sup (LBB) condition(2.15)infqh∈Qhsupvh∈Xh(qh,∇·vh)|vh|1∥q∥≥C>0.
FE-approximation
(2.16)infvh∈Xh|u-vh|1≤Chk∥u∥k+1,infqh∈Qh∥p-qh∥≤Chs+1∥p∥s+1fork≥0,s≥-1whenu∈Hk+1(Ω)∩X,p∈Hs+1(Ω)∩Q.
Inverse-estimate
(2.17)|vh|1≤Ch-1∥vh∥,∀vh∈Xh.
The well-known Taylor-Hood mixed FE is one such example satisfying Assumption 2.3.Estimates in (2.18)–(2.20) stated below are used in proving error estimates for time-dependent problems: for any n=0,1,…,M-1,
(2.18)∥θn+1-θnΔt∥k2≤CΔt-1∫tntn+1‍∥θt(t)∥k2dt,(2.19)∥θn+1/2-θ(tn+1/2)∥k2≤CΔt3∫tntn+1‍∥θtt(t)∥k2dt,(2.20)∥1Δt(θn+1-θn)-(θt)n+1/2∥k2≤CΔt3∫tntn+1‍∥θttt(t)∥k2dt,
where θ∈H1(Hk(Ω)), θ∈H2(Hk(Ω)), and θ∈H3(Hk(Ω)) is required, respectively, for some k≥-1. Each estimate (2.18)–(2.20) is a result of a Taylor expansion with integral remainder.These higher-order spatial (k≥2 or s≥1) and temporal estimates (2.18)–(2.20) require that the nonlocal compatibility condition addressed by Heywood and Rannacher in [26, 27] (and more recently, for example, by He in [28, 29] and He and Li in [30, 31]) is satisfied. Suppose, for example, that p0 is the solution of the (well-posed) Neumann problem
(2.21)Δp0=∇·(f0-u0·∇u0),inΩ,∇p0·n^|∂Ω=(Δu0+f0-u0·∇u0)·n^|∂Ω.
In order to avoid the accompanying factor min{t-1,1} in the error estimates contained herein, the following compatibility condition is necessarily required (e.g., see [27, Corollary 2.1]):
(2.22)∇p0|∂Ω=(Δu0+f0-u0·∇u0)|∂Ω.
Replacing (2.21) with (2.21)(a), (2.22) defines an overdetermined Neumann-type problem. Condition (2.22) is a nonlocal condition relating u0 and f0. Condition (2.22) is satisfied for several practical applications including start up from rest with zero force, u0=0, f0=0. In general, however, condition (2.22) cannot be verified. In this case, it is shown that, for example, |ut(·,t)|1,||u(·,t)||3→∞ast→0+.We finish with an approximation property of theL2-projection. Indeed, Assumption 2.4 holds for smooth enough Ω.Assumption 2.4.
Fixw∈V and let wh∈Vh be the unique solution satisfying (w-wh,vh) for all vh∈Vh. Then
(2.23)∥w-wh∥-m≤Chm+1infvh∈Xh|w-vh|1
for m=-1,0,1.Note that the infimum in (2.23) is over all Xh (see intermediate estimate (1.16) of Theorem II.1.1 in [23] for the corresponding estimate relating the spaces Vh and Xh).


## 2.1. Discrete Function Setting

Fixh>0. Let 𝒯h be a family of subdivisions (e.g., triangulation) of Ω¯⊂ℝd satisfying Ω¯=⋃E∈𝒯h‍E so that diameter(E)≤h and any two closed elements E1, E2∈𝒯h are either disjoint or share exactly one face, side, or vertex. See Chapter II, Appendix A in [23] for more on this subject in context of Stokes problem and [24] for a more general treatment. For example, 𝒯h consists of triangles for d=2 or tetrahedra for d=3 that are nondegenerate as h→0.LetXh⊂X and Qh⊂Q be a conforming velocity-pressure mixed FE space. For example, let Xh and Qh be continuous, piecewise (on each E∈𝒯h) polynomial spaces. The discretely divergence-free space is given by
(2.8)Vh={vh∈Xh:(qh,∇·vh)=0∀qh∈Qh}.
Note that in general Vh⊄V (e.g., for Taylor-Hood elements). In order to avoid stability issues arising when FE solutions are not exactly divergence free (i.e., when Vh⊄V), we introduce the explicitly skew-symmetric convective term
(2.9)bh(u,v,w):=12((u·∇v,w)-(u·∇w,v)),
so that
(2.10)bh(u,v,v)=0.
Note that bh(u,v,w)=(u·∇v,w) when u∈V. Moreover, the trilinear from bh(·,·,·) is continuous and skew-symmetric on X×X×X.Lemma 2.1.
Ifu,v,w∈X,
(2.11)|bh(u,v,w)|≤C(Ω){(∥u∥∥u∥1)1/2∥v∥1∥w∥1∥u∥∥v∥2∥w∥1∀v∈H2(Ω).Proof.
The proof of the first inequality can be found in [25]. The second follows from Hölder’s and Poincaré’s inequalities.For the time discretization, let0=t0<t1<⋯<tM-1=T<∞ be a discretization of the time interval [0,T] for a constant time step Δt=tn+1-tn. Write tn+1/2:=(tn+1+tn)/2,zn=z(tn) and, if z∈C0([tn,tn+1]), zn+1/2=(zn+1+zn)/2. Define
(2.12)∥u∥lq(m1,m2;Wpk(Ω)):={(Δt∑n=m1m2∥un∥k,pq)1/q,q∈[1,∞),maxm1≤n≤m2∥un∥k,p,q=∞,
for any 0≤n=m1,m1+1,…,m2≤M. Write ||u||lq(Wpk(Ω))=||u||lq(0,M;Wpk(Ω)). We say that u∈lq(m1,m2;Wpk(Ω)) if the associated norm defined above stays finite as Δt→0.The discrete Gronwall inequality is essential to the convergence analysis in Section4.2.Lemma 2.2.
LetD≥0 and κn, An, Bn, Cn≥0 for any integer n≥0 and satisfy
(2.13)AM+Δt∑n=0M‍Bn≤Δt∑n=0M‍κnAn+Δt∑n=0M‍Cn+D,∀M≥0.
Suppose that for all n, Δtκn<1 and set gn=(1-Δtκn)-1. Then,
(2.14)AM+Δt∑n=0M‍Bn≤exp(Δt∑n=0M‍gnκn)[Δt∑n=0M‍Cn+D],∀M≥0.Proof.
The proof follows from [26].


## 2.2. Approximation Theory

LetC>0 be a generic constant independent of h→0+. Preserving an abstract framework for the FE spaces, we assume that Xh×Qh inherit several fundamental approximation properties.Assumption 2.3.
The FE spacesXh×Qh satisfy:Uniform inf-sup (LBB) condition(2.15)infqh∈Qhsupvh∈Xh(qh,∇·vh)|vh|1∥q∥≥C>0.
FE-approximation
(2.16)infvh∈Xh|u-vh|1≤Chk∥u∥k+1,infqh∈Qh∥p-qh∥≤Chs+1∥p∥s+1fork≥0,s≥-1whenu∈Hk+1(Ω)∩X,p∈Hs+1(Ω)∩Q.
Inverse-estimate
(2.17)|vh|1≤Ch-1∥vh∥,∀vh∈Xh.
The well-known Taylor-Hood mixed FE is one such example satisfying Assumption 2.3.Estimates in (2.18)–(2.20) stated below are used in proving error estimates for time-dependent problems: for any n=0,1,…,M-1,
(2.18)∥θn+1-θnΔt∥k2≤CΔt-1∫tntn+1‍∥θt(t)∥k2dt,(2.19)∥θn+1/2-θ(tn+1/2)∥k2≤CΔt3∫tntn+1‍∥θtt(t)∥k2dt,(2.20)∥1Δt(θn+1-θn)-(θt)n+1/2∥k2≤CΔt3∫tntn+1‍∥θttt(t)∥k2dt,
where θ∈H1(Hk(Ω)), θ∈H2(Hk(Ω)), and θ∈H3(Hk(Ω)) is required, respectively, for some k≥-1. Each estimate (2.18)–(2.20) is a result of a Taylor expansion with integral remainder.These higher-order spatial (k≥2 or s≥1) and temporal estimates (2.18)–(2.20) require that the nonlocal compatibility condition addressed by Heywood and Rannacher in [26, 27] (and more recently, for example, by He in [28, 29] and He and Li in [30, 31]) is satisfied. Suppose, for example, that p0 is the solution of the (well-posed) Neumann problem
(2.21)Δp0=∇·(f0-u0·∇u0),inΩ,∇p0·n^|∂Ω=(Δu0+f0-u0·∇u0)·n^|∂Ω.
In order to avoid the accompanying factor min{t-1,1} in the error estimates contained herein, the following compatibility condition is necessarily required (e.g., see [27, Corollary 2.1]):
(2.22)∇p0|∂Ω=(Δu0+f0-u0·∇u0)|∂Ω.
Replacing (2.21) with (2.21)(a), (2.22) defines an overdetermined Neumann-type problem. Condition (2.22) is a nonlocal condition relating u0 and f0. Condition (2.22) is satisfied for several practical applications including start up from rest with zero force, u0=0, f0=0. In general, however, condition (2.22) cannot be verified. In this case, it is shown that, for example, |ut(·,t)|1,||u(·,t)||3→∞ast→0+.We finish with an approximation property of theL2-projection. Indeed, Assumption 2.4 holds for smooth enough Ω.Assumption 2.4.
Fixw∈V and let wh∈Vh be the unique solution satisfying (w-wh,vh) for all vh∈Vh. Then
(2.23)∥w-wh∥-m≤Chm+1infvh∈Xh|w-vh|1
for m=-1,0,1.Note that the infimum in (2.23) is over all Xh (see intermediate estimate (1.16) of Theorem II.1.1 in [23] for the corresponding estimate relating the spaces Vh and Xh).


## 3. Filters and Deconvolution

We prescribe the essential properties our filterG and deconvolution operator D in this section.Definition 3.1.
LetY be a Hilbert space and T:Y→Y. Write T≥0 if T is self-adjoint T=T′ and (Tv,v)Y≥0 for all v∈Y and call T symmetric nonnegative (snn). Write T>0 if T is self-adjoint T=T′ and (Tv,v)Y>0 for all 0≠v∈Y and call T symmetric positive definite (spd).
LetG=G(δ)>0 be a linear, bounded, compact operator on X representing a generic smoothing filter with filter radius δ>0:
(3.1)G:L2(Ω)→L2(Ω),ϕ¯:=Gϕ.
One example of this operator is the continuous differential filter G=A-1=(-δ2Δ+I)-1 (Definition 3.2), which is used, together with its discrete counterpart (Ah)-1 (Definition 3.3), for implementation of our numerical scheme (Section 5.2).Definition 3.2 (continuous differential filter).
Fixϕ∈L2(Ω). Then ϕ¯∈X is the unique solution of -δ2Δϕ¯+ϕ¯=ϕ with corresponding weak formulation
(3.2)δ2(∇ϕ¯,∇v)+(ϕ¯,v)=(ϕ,v),∀v∈X.
Set A=-δ2Δ+I so that A-1:L2(Ω)→X, defined by ϕ¯:=A-1ϕ, is well defined.Definition 3.3 (discrete differential filter).
Fixϕ∈L2(Ω). Then ϕ¯h∈Xh is the unique solution of the following
(3.3)δ2(∇ϕ¯h,∇vh)+(ϕ¯h,vh)=(ϕ,vh),∀vh∈Xh.
Set Ah=-δ2Δh+Πh so that (Ah)-1:L2(Ω)→Xh, defined by ϕ¯h:=(Ah)-1ϕ, is well defined. Here, Πh:L2(Ω)→Xh is the L2 projection and Δh:X→Xh the discrete Laplace operator satisfying the following
(3.4)(Πhϕ-ϕ,vh)=0,(Δhϕ,vh)=-(∇ϕ,∇vh)∀vh∈Xh.
It is well known thatA-1 and (Ah)-1 are each linear and bounded, A-1 is compact, and the spectrum of A and Ah (on X and Xh, resp.) is contained in [1,∞) and spectrum of A-1 and (Ah)-1 (on X and Xh, resp.) is contained in (0,1] so that
(3.5)A-1>0onX,(Ah)-1>0onXh.
For more detailed exposition on these operators, see [13].


### 3.1. A Family of Deconvolution Operators

We analyze (1.6) for stable, accurate deconvolution D of the smoothing filter G introduced in Section 3 so that DGu accurately approximates the smooth parts of u.Assumption 3.4 (continuous deconvolution operator).
Suppose thatD:X→X is linear, bounded, spd, and commutes with G so that
(3.6)∥DG∥≤1,|Dv¯|1≤d1|v|1∀v∈X,
for some constant d1>0. Moreover, suppose that D is parametrized by α>0, δ>0 so that
(3.7)∥(I-DG)v∥≤c1(δ,α)∥Δv∥→0,asα,δ→0
for smooth enough v∈X.
Note that the first estimate in (3.6) is required so that the spectral radius satisfies ρ(DG)≤1. The second estimate in (3.6) (which controls the H1-seminorm of DG) is required for the convergence analysis in Section 4.2.
Assumption3.5 prescribes properties of the discrete analogue Dh:Xh→Xh corresponding to the continuous deconvolution operator D:X→X (Assumptions 3.4).Assumption 3.5 (discrete deconvolution operator).
LetD satisfy Assumption 3.4. Let Gh:Xh→Xh be a discrete analogue of G that is linear, bounded, spd. Suppose that Dh:Xh→Xh is linear, bounded, spd, and commutes with Gh such that
(3.8)∥DhGh∥≤1,|Dhv¯h|1≤d1|v|1∀v∈X,
for some constant d1>0. Moreover, suppose that Dh is parametrized by α>0,δ>0,h such that D=D(h,δ,α) and
(3.9)∥(DG-DhGh)v∥≤c2(h,δ,α)∥v∥k+1→0,ash,δ,α→0,
for all v∈X∩Hk+1(Ω) for some k≥0.The estimates in (3.8) are motivated by the continuous case of (3.6). The approximation (3.9) is required for the convergence analysis in Section 4.2 (see Theorem 4.6, Corollary 4.7).Remark 3.6.
IfD=f(G) for some continuous map f:ℝ→ℝ, then commutativity is satisfied DG=GD. Tikhonov-Lavrentiev (modified) regularization with G=A-1, Gh=(Ah)-1 given by D=((1-α)A-1+αI)-1, Dh=((1-α)(Ah)-1+αΠh)-1 is one such example with f(x)=((1-α)x+α)-1 and d1=1, c1=αδ2, c2=αδ2hk+hk+1, see [13].Remark 3.7.
Lettingλk(·) denote the kth (ordered) eigenvalue of a given operator, commutativity of D and G provides λk(DG)=λk(D)λk(G) and similarly for the discrete operator DhGh.We next derive several important consequences ofD and Dh under Assumptions 3.4, 3.5 required in the forthcoming analysis.Lemma 3.8.
Suppose thatG, Gh, D, Dh satisfy Assumptions 3.4, 3.5. Then,
(3.10)∥Dv¯∥≤∥v∥,∥Dhv¯h∥≤∥v∥∀v∈L2(Ω).Proof.
For the continuous operator,(3.11)∥Dv¯∥≤∥DG∥∥v∥.
Then, (3.10)(a) follows from Assumption 3.4, and (3.10)(b) is derived similarly applying Assumption 3.5 instead.Lemma 3.9.
Suppose thatG, Gh, D, Dh satisfy Assumptions 3.4, 3.5. Then, the spectrum of both DG and DhGh are contained in [0,1] so that
(3.12)∥I-DG∥≤1,∥I-DhGh∥≤1.
As a consequence,
(3.13)∥v∥⋆2:=(v-Dv-,v)∥v∥⋆h2:=(v-Dhv-h,v)}≥0,∀v∈L2(Ω).Proof.
Assumptions3.4, 3.5 guarantee that the spectral radius ρ(DG)≤1 and ρ(DhGh)≤1. Also, D>0 and G>0 and commute so that DG≥0. Similarly, DhGh≥0. Therefore, the spectrum of DG, DhGh≥0 are each contained in [0,1]. So, I-DG, I-DhGh≥0 have spectrum contained in [0,1] which ensures the non-negativity of both ||·||⋆ and ||·||⋆h.


### 3.2. Iterated Deconvolution

One can show, by eliminating intermediate steps in the definition of the iterated regularization operatorDj in (1.9) with base operator D satisfying Assumption 3.4, that
(3.14)Dj=D∑i=0j‍(FD)i,F:=D-1-G.
Similarly, the discrete iterated regularization operator Djh with discrete base operator Dh satisfying Assumption 3.5, is given by the following
(3.15)Djh=Dh∑i=0j‍(FhDh)i,Fh:=(Dh)-1-Gh.
We next show that Dj and Djh for j>0 inherit several important properties from D and Dh, respectively, via Assumption 3.5.Proposition 3.10.
Fixj∈ℕ. Then Dj:X→X defined by (3.14) satisfies Assumption 3.4. In particular, Dj>0 is linear, bounded, commutes with G and satisfies (3.6)(a). Estimate (3.6)(b) is replaced by the following
(3.16)|Djv¯|1≤d1,j|v|1∀v∈X
for some constant d1,j>0. Estimate (3.7) is replaced by the following
(3.17)∥(I-DjG)v∥≤c1,j(δ,α)∥Δv∥2→0,asα,δ→0.
Moreover, d1,j≤∑i=0j‍d1i and c1,j≤∑i=0j‍c1i.Proof.
First notice thatDj is linear and bounded since it is a linear combination of linear and bounded operators D(FD)i=D(I-DG)i, for i=0,1,…,j. Moreover, since G commutes with D, it follows that G commutes with D(I-DG)i and hence with Dj. Next, Dj is a sum of spd and snn operators D>0, D(I-DG)i≥0. Hence, Dj>0. Next, notice that
(3.18)DjG=(∑i=0j‍(I-DG)i)DG=(I+(I-DG)+⋯+(I-DG)j)DG.
Letting λk(·) denote the kth (ordered) eigenvalue of a given operator, we can characterize the spectrum of Dj by summing the resulting finite geometric series (3.18) to get
(3.19)λk(DjG)=λk(D)λk(G)∑i=0j‍(1-λk(DG))i=(1-(1-λk(DG))j+1).
Then under Assumption 3.4, Lemma 3.9 with (3.19) implies that 0≤λk(DjG)≤||DjG||≤1. Hence, Dj satisfies (3.6)(a). Expanding the terms in (3.18) as powers of DG, we see that (3.18) can be written as a polynomial (with coefficients ai) in DG, so that
(3.20)∇Djv¯=∑i=0j‍ai∇(DG)iv,|Djv¯|1≤∑i=0j‍d1i|v|1,
since |Dv¯|1≤d1|v|1 can be applied successfully. Therefore (3.16) follows with d1,j=∑i=0j‍d1i. Next, start with (3.18) to get
(3.21)∥(I-DjG)v∥=∥((I-DG)v+DG(I-DG)v+⋯+DG(I-DG)jv)∥≤∥(I-DG)v∥+∥DG∥∥(I-DG)v∥+⋯+∥DG∥∥I-DG∥j-1∥(I-DG)v∥.
Estimate (3.17) follows by noting ||DG||≤1, ||I-DG||≤1, and by Assumption 3.5, ||(I-DG)v||≤c1||Δv||.Proposition 3.11.
Fixj∈ℕ. Then Djh:Xh→Xh defined by (3.15) satisfies Assumption 3.5. In particular, Djh>0 is linear, bounded, commutes with Gh and satisfies (3.8)(a). Estimate (3.8)(b) is replaced by the following
(3.22)|Djhv¯h|1≤d1,j|v|1∀v∈X
for some constant d1,j>0. Estimate (3.9) is replaced by the following
(3.23)∥(DjG-DjhGh)v∥≤c2,j(h,δ,α)∥v∥k+1→0,ash,δ,α→0
for any v∈X∩Hk+1(Ω) for some k≥0. Moreover, c2,j≤β(j)c2 for some constant β=β(j)>0.Proof.
The first two assertions follow similarly as in the previous proof of Proposition3.10. To prove (3.23), we start by writing
(3.24)DjhGh=(∑i=0j‍(I-DhGh)i)DhGh,
and then subtract (3.24) from (3.18) to get
(3.25)DjG-DjhGh=Λj(DG-DhGh)+(Λj-Λjh)DhGh,
where
(3.26)Λj=∑i=0j‍(I-DG)i,Λjh=∑i=0j‍(I-DhGh)i.
Then taking norms across (3.25), we get
(3.27)∥(DjG-DjhGh)v∥=∥Λj∥∥(DG-DhGh)v∥+∥DhGh∥∥(Λj-Λjh)v∥.
Notice that ||I-DG||≤1 so that ||Λj||≤j+1. Moreover, ||(DG-DhGh)v||≤c2||v||k+1 via Assumption 3.5. Next, using the binomial theorem and factoring, we get
(3.28)∥(Λj-Λjh)v∥=∥∑i=0j‍j!i!(j-i)!(-1)i[(DG)i-(DhGh)iv]∥=∥∑i=0j‍j!i!(j-i)!(-1)i[∑n=0i‍(DG)n(DhGh)n-i](DG-DhGh)v∥.
Then, applying ||DG||≤1, ||DhGh||≤1 to (3.28) provides
(3.29)∥(Λj-Λjh)v∥=(∑i=0j‍j!i(i)!(j-i)!)∥(DG-DhGh)v∥.
Again, ||(DG-DhGh)||≤c2||v||k+1 via Assumption 3.5. So, we combine these above results to conclude (3.23) with β(j)=∑i=0j‍j!i/(i)!(j-i)!.


### 3.3. Tikhonov-Lavrentiev Regularization

We provide two examples of discrete deconvolution operatorsDh to make the abstract formulation in the previous section more concrete. The Tikhonov-Lavrentiev and modified Tikhonov-Lavrentiev operator (for linear, compact G>0) is given by the following
(3.30)Tikhonov-Lavrentiev⇒Dα,0=(G+αI)-1modified  Tikhonov-Lavrentiev⇒Dα,0=((1-α)G+αI)-1.Definition 3.12 ((weak) modified Tikhonov-Lavrentiev deconvolution).
Fixα>0. Let G=A-1. For any w∈X, let ω0:=Dα,0w¯∈X be the unique solution of
(3.31)αδ2(∇ω0,∇v)+(ω0,v)=(w,v),∀v∈X.Definition 3.13 ((discrete) modified Tikhonov-Lavrentiev deconvolution).
Fixα>0. Let Gh=(Ah)-1 and Dα,0h=((1-α)(Ah)-1+αΠh)-1. For any w∈X, let ω0h:=Dα,0hw¯h∈Xh be the unique solution of
(3.32)αδ2(∇ω0h,∇vh)+(ω0h,vh)=(w,vh),∀vh∈Xh.
The iterated modified Tikhonov-Lavrentiev operator (for linear, compactG>0) is obtained from the Tikhonov-Lavrentiev operator with updates via (1.9):
(3.33)Iterated  Tikhonov-Lavrentiev⇒Dα,j=Dα,0∑i=0j‍(αDα,0)i,Iterated  modified  Tikhonov-Lavrentiev⇒Dα,j=Dα,0∑i=0j‍(α(I-G)Dα,0)i.Definition 3.14 (iterated modified Tikhonov-Lavrentiev deconvolution (weak)).
Fixα>0 and J∈ℕ. Let G=A-1. Define ω-1=0, then for any w∈X and j=0,1,…,J, let ωj:=Dα,jw¯∈X be the unique solution of
(3.34)αδ2(∇ωj,∇v)+(ωj,v)=(w,v)+αδ2(∇ωj-1,∇v),∀v∈X.Definition 3.15 (iterated modified Tikhonov-Lavrentiev deconvolution (discrete)).
Fixα>0 and J∈ℕ. Let Gh=(Ah)-1, and Dα,jh=Dα,0h∑i=0j‍(α(Πh-(Ah)-1)Dα,0h)i. Define ω-1h=0, then for any w∈X and j=0,1,…,J, let ωjh:=Dα,jhw¯h∈Xh be the unique solution of
(3.35)αδ2(∇ωjh,∇vh)+(ωjh,vh)=(w,vh)+αδ2(∇ωj-1h,∇vh),∀vh∈Xh.


## 3.1. A Family of Deconvolution Operators

We analyze (1.6) for stable, accurate deconvolution D of the smoothing filter G introduced in Section 3 so that DGu accurately approximates the smooth parts of u.Assumption 3.4 (continuous deconvolution operator).
Suppose thatD:X→X is linear, bounded, spd, and commutes with G so that
(3.6)∥DG∥≤1,|Dv¯|1≤d1|v|1∀v∈X,
for some constant d1>0. Moreover, suppose that D is parametrized by α>0, δ>0 so that
(3.7)∥(I-DG)v∥≤c1(δ,α)∥Δv∥→0,asα,δ→0
for smooth enough v∈X.
Note that the first estimate in (3.6) is required so that the spectral radius satisfies ρ(DG)≤1. The second estimate in (3.6) (which controls the H1-seminorm of DG) is required for the convergence analysis in Section 4.2.
Assumption3.5 prescribes properties of the discrete analogue Dh:Xh→Xh corresponding to the continuous deconvolution operator D:X→X (Assumptions 3.4).Assumption 3.5 (discrete deconvolution operator).
LetD satisfy Assumption 3.4. Let Gh:Xh→Xh be a discrete analogue of G that is linear, bounded, spd. Suppose that Dh:Xh→Xh is linear, bounded, spd, and commutes with Gh such that
(3.8)∥DhGh∥≤1,|Dhv¯h|1≤d1|v|1∀v∈X,
for some constant d1>0. Moreover, suppose that Dh is parametrized by α>0,δ>0,h such that D=D(h,δ,α) and
(3.9)∥(DG-DhGh)v∥≤c2(h,δ,α)∥v∥k+1→0,ash,δ,α→0,
for all v∈X∩Hk+1(Ω) for some k≥0.The estimates in (3.8) are motivated by the continuous case of (3.6). The approximation (3.9) is required for the convergence analysis in Section 4.2 (see Theorem 4.6, Corollary 4.7).Remark 3.6.
IfD=f(G) for some continuous map f:ℝ→ℝ, then commutativity is satisfied DG=GD. Tikhonov-Lavrentiev (modified) regularization with G=A-1, Gh=(Ah)-1 given by D=((1-α)A-1+αI)-1, Dh=((1-α)(Ah)-1+αΠh)-1 is one such example with f(x)=((1-α)x+α)-1 and d1=1, c1=αδ2, c2=αδ2hk+hk+1, see [13].Remark 3.7.
Lettingλk(·) denote the kth (ordered) eigenvalue of a given operator, commutativity of D and G provides λk(DG)=λk(D)λk(G) and similarly for the discrete operator DhGh.We next derive several important consequences ofD and Dh under Assumptions 3.4, 3.5 required in the forthcoming analysis.Lemma 3.8.
Suppose thatG, Gh, D, Dh satisfy Assumptions 3.4, 3.5. Then,
(3.10)∥Dv¯∥≤∥v∥,∥Dhv¯h∥≤∥v∥∀v∈L2(Ω).Proof.
For the continuous operator,(3.11)∥Dv¯∥≤∥DG∥∥v∥.
Then, (3.10)(a) follows from Assumption 3.4, and (3.10)(b) is derived similarly applying Assumption 3.5 instead.Lemma 3.9.
Suppose thatG, Gh, D, Dh satisfy Assumptions 3.4, 3.5. Then, the spectrum of both DG and DhGh are contained in [0,1] so that
(3.12)∥I-DG∥≤1,∥I-DhGh∥≤1.
As a consequence,
(3.13)∥v∥⋆2:=(v-Dv-,v)∥v∥⋆h2:=(v-Dhv-h,v)}≥0,∀v∈L2(Ω).Proof.
Assumptions3.4, 3.5 guarantee that the spectral radius ρ(DG)≤1 and ρ(DhGh)≤1. Also, D>0 and G>0 and commute so that DG≥0. Similarly, DhGh≥0. Therefore, the spectrum of DG, DhGh≥0 are each contained in [0,1]. So, I-DG, I-DhGh≥0 have spectrum contained in [0,1] which ensures the non-negativity of both ||·||⋆ and ||·||⋆h.


## 3.2. Iterated Deconvolution

One can show, by eliminating intermediate steps in the definition of the iterated regularization operatorDj in (1.9) with base operator D satisfying Assumption 3.4, that
(3.14)Dj=D∑i=0j‍(FD)i,F:=D-1-G.
Similarly, the discrete iterated regularization operator Djh with discrete base operator Dh satisfying Assumption 3.5, is given by the following
(3.15)Djh=Dh∑i=0j‍(FhDh)i,Fh:=(Dh)-1-Gh.
We next show that Dj and Djh for j>0 inherit several important properties from D and Dh, respectively, via Assumption 3.5.Proposition 3.10.
Fixj∈ℕ. Then Dj:X→X defined by (3.14) satisfies Assumption 3.4. In particular, Dj>0 is linear, bounded, commutes with G and satisfies (3.6)(a). Estimate (3.6)(b) is replaced by the following
(3.16)|Djv¯|1≤d1,j|v|1∀v∈X
for some constant d1,j>0. Estimate (3.7) is replaced by the following
(3.17)∥(I-DjG)v∥≤c1,j(δ,α)∥Δv∥2→0,asα,δ→0.
Moreover, d1,j≤∑i=0j‍d1i and c1,j≤∑i=0j‍c1i.Proof.
First notice thatDj is linear and bounded since it is a linear combination of linear and bounded operators D(FD)i=D(I-DG)i, for i=0,1,…,j. Moreover, since G commutes with D, it follows that G commutes with D(I-DG)i and hence with Dj. Next, Dj is a sum of spd and snn operators D>0, D(I-DG)i≥0. Hence, Dj>0. Next, notice that
(3.18)DjG=(∑i=0j‍(I-DG)i)DG=(I+(I-DG)+⋯+(I-DG)j)DG.
Letting λk(·) denote the kth (ordered) eigenvalue of a given operator, we can characterize the spectrum of Dj by summing the resulting finite geometric series (3.18) to get
(3.19)λk(DjG)=λk(D)λk(G)∑i=0j‍(1-λk(DG))i=(1-(1-λk(DG))j+1).
Then under Assumption 3.4, Lemma 3.9 with (3.19) implies that 0≤λk(DjG)≤||DjG||≤1. Hence, Dj satisfies (3.6)(a). Expanding the terms in (3.18) as powers of DG, we see that (3.18) can be written as a polynomial (with coefficients ai) in DG, so that
(3.20)∇Djv¯=∑i=0j‍ai∇(DG)iv,|Djv¯|1≤∑i=0j‍d1i|v|1,
since |Dv¯|1≤d1|v|1 can be applied successfully. Therefore (3.16) follows with d1,j=∑i=0j‍d1i. Next, start with (3.18) to get
(3.21)∥(I-DjG)v∥=∥((I-DG)v+DG(I-DG)v+⋯+DG(I-DG)jv)∥≤∥(I-DG)v∥+∥DG∥∥(I-DG)v∥+⋯+∥DG∥∥I-DG∥j-1∥(I-DG)v∥.
Estimate (3.17) follows by noting ||DG||≤1, ||I-DG||≤1, and by Assumption 3.5, ||(I-DG)v||≤c1||Δv||.Proposition 3.11.
Fixj∈ℕ. Then Djh:Xh→Xh defined by (3.15) satisfies Assumption 3.5. In particular, Djh>0 is linear, bounded, commutes with Gh and satisfies (3.8)(a). Estimate (3.8)(b) is replaced by the following
(3.22)|Djhv¯h|1≤d1,j|v|1∀v∈X
for some constant d1,j>0. Estimate (3.9) is replaced by the following
(3.23)∥(DjG-DjhGh)v∥≤c2,j(h,δ,α)∥v∥k+1→0,ash,δ,α→0
for any v∈X∩Hk+1(Ω) for some k≥0. Moreover, c2,j≤β(j)c2 for some constant β=β(j)>0.Proof.
The first two assertions follow similarly as in the previous proof of Proposition3.10. To prove (3.23), we start by writing
(3.24)DjhGh=(∑i=0j‍(I-DhGh)i)DhGh,
and then subtract (3.24) from (3.18) to get
(3.25)DjG-DjhGh=Λj(DG-DhGh)+(Λj-Λjh)DhGh,
where
(3.26)Λj=∑i=0j‍(I-DG)i,Λjh=∑i=0j‍(I-DhGh)i.
Then taking norms across (3.25), we get
(3.27)∥(DjG-DjhGh)v∥=∥Λj∥∥(DG-DhGh)v∥+∥DhGh∥∥(Λj-Λjh)v∥.
Notice that ||I-DG||≤1 so that ||Λj||≤j+1. Moreover, ||(DG-DhGh)v||≤c2||v||k+1 via Assumption 3.5. Next, using the binomial theorem and factoring, we get
(3.28)∥(Λj-Λjh)v∥=∥∑i=0j‍j!i!(j-i)!(-1)i[(DG)i-(DhGh)iv]∥=∥∑i=0j‍j!i!(j-i)!(-1)i[∑n=0i‍(DG)n(DhGh)n-i](DG-DhGh)v∥.
Then, applying ||DG||≤1, ||DhGh||≤1 to (3.28) provides
(3.29)∥(Λj-Λjh)v∥=(∑i=0j‍j!i(i)!(j-i)!)∥(DG-DhGh)v∥.
Again, ||(DG-DhGh)||≤c2||v||k+1 via Assumption 3.5. So, we combine these above results to conclude (3.23) with β(j)=∑i=0j‍j!i/(i)!(j-i)!.


## 3.3. Tikhonov-Lavrentiev Regularization

We provide two examples of discrete deconvolution operatorsDh to make the abstract formulation in the previous section more concrete. The Tikhonov-Lavrentiev and modified Tikhonov-Lavrentiev operator (for linear, compact G>0) is given by the following
(3.30)Tikhonov-Lavrentiev⇒Dα,0=(G+αI)-1modified  Tikhonov-Lavrentiev⇒Dα,0=((1-α)G+αI)-1.Definition 3.12 ((weak) modified Tikhonov-Lavrentiev deconvolution).
Fixα>0. Let G=A-1. For any w∈X, let ω0:=Dα,0w¯∈X be the unique solution of
(3.31)αδ2(∇ω0,∇v)+(ω0,v)=(w,v),∀v∈X.Definition 3.13 ((discrete) modified Tikhonov-Lavrentiev deconvolution).
Fixα>0. Let Gh=(Ah)-1 and Dα,0h=((1-α)(Ah)-1+αΠh)-1. For any w∈X, let ω0h:=Dα,0hw¯h∈Xh be the unique solution of
(3.32)αδ2(∇ω0h,∇vh)+(ω0h,vh)=(w,vh),∀vh∈Xh.
The iterated modified Tikhonov-Lavrentiev operator (for linear, compactG>0) is obtained from the Tikhonov-Lavrentiev operator with updates via (1.9):
(3.33)Iterated  Tikhonov-Lavrentiev⇒Dα,j=Dα,0∑i=0j‍(αDα,0)i,Iterated  modified  Tikhonov-Lavrentiev⇒Dα,j=Dα,0∑i=0j‍(α(I-G)Dα,0)i.Definition 3.14 (iterated modified Tikhonov-Lavrentiev deconvolution (weak)).
Fixα>0 and J∈ℕ. Let G=A-1. Define ω-1=0, then for any w∈X and j=0,1,…,J, let ωj:=Dα,jw¯∈X be the unique solution of
(3.34)αδ2(∇ωj,∇v)+(ωj,v)=(w,v)+αδ2(∇ωj-1,∇v),∀v∈X.Definition 3.15 (iterated modified Tikhonov-Lavrentiev deconvolution (discrete)).
Fixα>0 and J∈ℕ. Let Gh=(Ah)-1, and Dα,jh=Dα,0h∑i=0j‍(α(Πh-(Ah)-1)Dα,0h)i. Define ω-1h=0, then for any w∈X and j=0,1,…,J, let ωjh:=Dα,jhw¯h∈Xh be the unique solution of
(3.35)αδ2(∇ωjh,∇vh)+(ωjh,vh)=(w,vh)+αδ2(∇ωj-1h,∇vh),∀vh∈Xh.


## 4. Well Posedness of the Fully Discrete Model

We now state the proposed algorithm.Problem 1 (CNFE for Leray-deconvolution).
Let(w0,π0)∈(Xh,Qh). Then, for each n=0,1,…,M-1, find (wn+1h,πn+1h)∈(Xh,Qh) satisfying
(4.1)1Δt(wn+1h-wnh,vh)+bh(φh(wn+1/2h),wn+1/2h,vh)-(πn+1/2h,∇·vh)+Re-1(∇wn+1/2h,∇vh)+χ(wn+1/2h-φh(wn+1/2h),vh)=(fn+1/2,vh),∀vh∈Xh,(4.2)(∇·wn+1h,qh)=0,∀qh∈Qh,
where φh(wn+1/2h)=Dhwn+1/2h¯h.Notice that(qn+1/2h,∇·vh)=0 when vh∈Vh so that the problem of finding wn+1h∈Vh satisfying
(4.3)1Δt(wn+1h-wnh,vh)+bh(φh(wn+1/2h),wn+1/2h,vh)+Re-1(∇wn+1/2h,∇vh)+χ(wn+1/2h-φh(wn+1/2h),vh)=(fn+1/2,vh),∀vh∈Vh.


### 4.1. Well Posedness

We establish existence ofw at each time step of (4.3) by Leray-Schauder’s fixed-point theorem.Lemma 4.1.
Let(4.4)a(θh,vh)=Δt2Re(∇θh,∇vh)+χΔt2(θh-φh(θh),vh),ly(vh)=(y,vh).
for any y∈X′ and θh, vh∈Vh. Suppose that Dh satisfies Assumption 3.5. Then a(·,·):Vh×Vh→ℝ is a continuous and coercive bilinear form and ly(·):Vh→ℝ is a linear, continuous functional.Proof.
Linearity forly(·) is obvious, and continuity follows from an application of Hölder’s inequality. Continuity for a(·,·) also follows from Hölder’s inequality and Assumption 3.5. Coercivity is proven by application of (3.13).Lemma 4.2.
LetT:X′→Vh be such that, for any y∈X′, θh:=T(y) solves
(4.5)a(θh,vh)=ly(vh),∀vh∈Vh.
Then T is a well-defined, linear, bounded operator.Proof.
Linearity is clear. The results of Lemma4.1, and the Lax-Milgram theorem prove the rest.Lemma 4.3.
Fixn=0,1,…,M-1. Let wnh be a solution of Problem 1 and let N:Vh→X′ satisfy, for any θh∈Vh,
(4.6)(N(θh),vh)=-(θh-2wnh,vh)-Δt4bh(φh(θh),θh,vh)+Δt(fn+1/2,vh)=:c(θh,vh),∀vh∈Vh.
Then N(θh) is well-defined, bounded, and continuous.Proof.
For eachθh∈Vh, the map vh∈Vh↦c(θh,vh) is a bounded, linear functional (apply Hölder’s inequality and (2.11)). Since Vh is a Hilbert space, we conclude that N(θh) is well defined, by the Riesz-Representation theorem. Moreover, N(θh) is bounded on Vh and since the underlying function space is finite dimensional, continuity follows.Lemma 4.4.
Fixn∈ℕ. Let F:Vh→Vh be defined such that F(θh)=(T∘N)(θh). Then, F is a compact operator.Proof.
N
(
·
) is a compact operator (continuous on a finite dimensional function space). Thus, F is a continuous composition of a compact operator and hence compact itself.Theorem 4.5 (well posedness).
Fixn=0,1,2,…,M-1<∞. There exists (wnh,πnh)∈Xh×Qh satisfying Problem 1. Moreover,
(4.7)∥wmh∥2+12ReΔt∑n=0m-1‍|wn+1/2h|12+χΔt∑n=0m-1‍∥wn+1/2h∥⋆h2≤∥w0h∥2+ΔtRe2∑n=0m-1‍∥fn+1/2∥-12,
for all integers 1≤m≤M, independent of Δt>0.Proof.
First, assume that(wn+1h,qn+1h) is a solution to (4.1), (4.2). Set vh=wn+1/2h in (4.1) so that skew-symmetry of the nonlinear term provides
(4.8)12Δt(∥wn+1h∥2-∥wmh∥2)+Re-1|wn+1/2h|12+χ(wn+1/2h-φh(wn+1/2h),wn+1/2h)=(fn+1/2,wn+1/2h).
Duality of X′×X with Young’s inequality implies
(4.9)(fn+1/2,wn+1/2h)≤Re2∥fn+1/2∥-12+Re2|wn+1/2h|12.
From (3.13), we have
(4.10)∥wn+1/2h∥⋆h2=(wn+1/2h-φh(wn+1/2h),wn+1/2h)≥0.
Then applying (4.9), (4.10) to (4.8), combining-like terms and simplifying provides
(4.11)12Δt(∥wn+1h∥2-∥wnh∥2)+Re2|wn+1/2h|12+χ∥wn+1/2h∥⋆h2≤Re2∥fn+1/2∥-12.
Summing from n=0 to m-1, we get the desired bound.
Next, letWnh=wn+1h+wnh. Showing that Wnh=F(Wnh) has a fixed point will ensure existence of solutions to (4.3). Indeed, if we can show that W0h=F(W0h), then since w0h is given initial data, existence of w1h is immediate. Induction can be applied to prove existence of (wnh)1≤n≤M. To this end, since F is compact, it is enough to show (via Leray Schauder) that any solution Wn,λh of the fixed-point problem Wn,λh=λF(Wn,λh) is uniformly bounded with respect to 0≤λ≤1. Hence, we consider
(4.12)a(Wn,λh,vh)=λ(N(Wn,λh),vh).
Test withvh=Wn,λh, use skew-symmetry of the trilinear form and properties of Dh given in Assumption 3.5 and (3.13) to get
(4.13)λ∥Wn,λh∥2+Δt2Re|Wn,λh|12+χΔt2∥Wn,λh∥⋆h2≤2λ(wnh,Wn,λh)+λΔt(fn+1/2,Wn,λh).
Duality of X′×X followed by Young’s inequality implies
(4.14)λΔt(fn+1/2,Wn,λh)≤ΔtRe∥fn+1/2∥-12+Δt4Re|Wn,λh|12.
Since wnh∈L2(Ω) from the a priori estimate (4.7), we apply Hölder’s and Young’s inequalities to get
(4.15)2λ(wnh,Wn,λh)≤2∥wnh∥2+λ2∥Wn,λh∥2.
Applying estimates (4.14), (4.15) to (4.13) we get that |Wn,λh|1≤C<∞ independent of λ. By the Leray-Schauder fixed point theorem, given wnh, there exists a solution to the fixed-point theorem Wnh=F(Wnh). By the induction argument noted above, there exists a solution wnh for each n=0,1,2,…,M-1 to (4.3). Existence of an associated discrete pressure follows by a classical argument, since the pair (Xh,Qh) satisfies the discrete inf-sup condition (2.15).


### 4.2. Convergence Analysis

Under usual regularity assumptions, we summarize the main convergence estimate in Theorem4.6. Suppose that D represents deconvolution with J-updates.Theorem 4.6.
Suppose that(u,p) are strong solutions to (2.3), (2.4), (2.5) and that G, Gh, D, Dh satisfy Assumptions 3.4, 3.5. Suppose further that u∈l2(H2(Ω)∩V)∩l∞(X), ut∈l2(X′)∩L2(X′), Δβ+1u∈L2(Ω) for some 0≤β≤J+1, p∈l2(Q). If
(4.16)CReΔt∥un∥22<1,∀n=0,1,…,M
then,
(4.17)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤C(∥u0-w0h∥2+E+∥p∥l2(L2(Ω))2+(∥∇u∥l∞(L2(Ω))2+χ2)∥u∥l2(H1(Ω))2+∥u∥l∞(L2(Ω))2+∥u∥l2(L2(Ω))2+(χ2+∥u∥l∞(H2(Ω))2)c1,J2∥Δβu∥2+c2,j2∥∇u∥l2(L2(Ω))2),
where E>0 is given in (4.52).Corollary 4.7 (convergence estimate).
Under the assumptions of Theorem4.6, suppose further that (u,p) satisfy the assumptions for (2.16) for some k≥1 and s≥0, u∈l∞(Hk(Ω))∩l2(Hk+1(Ω)), ut∈L2(Hk-1(Ω)∩H1(Ω))∩l∞(H1(Ω)), utt∈L2(L2(Ω)), uttt∈L2(X′), and p∈l2(Hs+1(Ω)). If ||u0-w0h||≤C(hk+(αδ2)β) then
(4.18)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤C(h2k+h2s+2+Δt4+c1,J(δ,α)2+c2,J(h,δ,α)2).Proof of Theorem4.6, Corollary 4.7.
Suppose thatu, p satisfying (2.3), (2.4) also satisfy u∈C0(V), ut∈C0(X′), and p∈C0(Q) so that, for each n=n0,n0+1,…,M-1,
(4.19)((ut)n+1/2,v)+12(un+1·∇un+1,v)+12(un·∇un,v)+Re-1(∇un+1/2,∇v)-(pn+1/2,∇·v)=(fn+1/2,v),∀v∈X,∇·un+1=0.
The consistency error for the time-discretization τn(1)(u,p;vh) and regularization/time-relaxation error τn(2)(u,p;vh) are given by, for n=0,1,…,M-1,
(4.20)τn(1)(u,p;vh):=(un+1-unΔt-(ut)n+1/2,vh)+(un+1/2·∇un+1/2,vh)-12(un+1·∇un+1,vh)-12(un·∇un,v),τn(2)(u,p;vh):=-bh(un+1/2-φh(un+1/2),un+1/2,vh)+χ(un+1/2-φh(un+1/2),vh),
where vh∈Xh. Write τn:=τn(1)+τn(2). Using (4.20), rewrite (4.19) in a form conducive to analyzing the error between the continuous and discrete models:
(4.21)(un+1-unΔt,vh)+bh(φh(un+1/2),un+1/2,vh)+Re-1(∇un+1/2,∇vh)-(pn+1/2,∇·vh)+χ(un+1/2-φh(un+1/2),vh)=(fn+1/2,vh)+τn(u,p;vh).
Let v~h=Unh be the L2-projection of u(·,tn) so that (ηn+1-ηn,ϕn+1/2h)=0. Decompose the velocity error
(4.22)en=wnh-un=ϕnh-ηn,ϕnh=wnh-Unh,ηn=un-Unh,
where Unh∈Vh. Fix q~n+1/2h∈Qh. Note that (q~n+1/2h,∇·vh)=0 for any vh∈Vh. Subtract (4.21) from (4.3), apply (3.13)(b), and test with vh=ϕn+1/2h to get
(4.23)12Δt(∥ϕn+1h∥2-∥ϕnh∥2)+Re-1|ϕn+1/2h|12+χ∥ϕn+1/2h∥⋆h2=-(q~n+1/2h-pn+1/2,∇·ϕn+1/2h)-bh(φh(ϕn+1/2h),un+1/2,ϕn+1/2h)=+bh(φh(ηn+1/2),un+1/2,ϕn+1/2h)+bh(φh(wn+1/2h),ηn+1/2,ϕn+1/2h)=+Re-1(∇ηn+1/2,∇ϕn+1/2h)+χ(ηn+1/2-φh(ηn+1/2),ϕn+1/2h)-τn(u,p;ϕn+1/2h).
(Spatial discretization error): Fix ε>0. First, apply Hölder’s and Young’s inequalities to get
(4.24)|Re-1(∇ηn+1/2,∇ϕn+1/2h)+(pn+1/2-q~n+1/2h,∇·ϕn+1/2h)|≤CRe-1|ηn+1/2|12+CRe∥pn+1/2-q~n+1/2h∥2+1εRe|ϕn+1/2h|12.
Apply (3.12) and duality estimate on X×X′ to get
(4.25)|χ(ηn+1/2-φh(ηn+1/2),ϕn+1/2h)|≤Cχ2Re∥ηn+1/2∥-12+1εRe|ϕn+1/2h|1.
We bound the convective terms next. First, u∈H2(Ω) and estimate (2.11)(b) give
(4.26)|bh(φh(ϕn+1/2h),un+1/2,ϕn+1/2h)|≤CRe∥φh(ϕn+1/2h)∥2∥un+1/2∥22+1εRe|ϕn+1/2h|12
and u∈L∞(X) with (2.11)(a) give
(4.27)|bh(φh(ηn+1/2),un+1/2,ϕn+1/2h)|≤CRe∥∇u∥L∞(L2(Ω))2|φh(ηn+1/2)|12+1εRe|ϕn+1/2h|12.
Next, rewrite the remaining nonlinear term
(4.28)bh(φh(wn+1/2h),ηn+1/2,ϕn+1/2h)=bh(φh(un+1/2),ηn+1/2,ϕn+1/2h)-bh(φh(ηn+1/2),ηn+1/2,ϕn+1/2h)+bh(φh(ϕn+1/2h),ηn+1/2,ϕn+1/2h).
Once again, u∈L∞(X) and (2.11)(a) give
(4.29)|bh(φh(un+1/2),ηn+1/2,ϕn+1/2h)+bh(φh(ηn+1/2),ηn+1/2,ϕn+1/2h)|≤CRe∥∇φh(u)∥L∞(L2(Ω))∥φh(u)∥L∞(L2(Ω))|ηn+1/2|12≤+CRe|φh(ηn+1/2)|1|φh(ηn+1/2)|1|ηn+1/2|12+1εRe|ϕn+1/2h|12.
Lastly, estimate (2.11)(a) and inverse inequality (2.17) give (4.30)|bh(φh(ϕn+1/2h),ηn+1/2,ϕn+1/2h)|≤CReh-1∥φh(ϕn+1/2h)∥2|ηn+1/2|12+1εRe|ϕn+1/2h|12.
Apply estimates ||φh(v)||≤||v|| and |φh(v)|1≤d1|v|1 from Assumption 3.5 along with estimates (4.26)–(4.30) to get
(4.31)|bh(φh(ϕn+1/2h),un+1/2,ϕn+1/2h)+bh(φh(ηn+1/2),un+1/2,ϕn+1/2h)+bh(φh(wn+1/2h),ηn+1/2,ϕn+1/2h)|≤4εRe|ϕn+1/2h|12+CRe(∥un+1/2∥22+h-1|ηn+1/2|12)∥ϕn+1/2h∥2≤+Cd1Re(∥∇u∥L∞(L2(Ω))2+|ηn+1/2|12)|ηn+1/2|12.
(Time discretization error): First, apply duality estimate on X×X′ to get
(4.32)|(un+1-unΔt-(ut)n+1/2,ϕn+1/2h)|≤CRe∥un+1-unΔt-(ut)n+1/2∥-12+1εRe|ϕn+1/2h|12.
Taylor-expansion about tn+1/2 with integral remainder gives
(4.33)12(un+1·∇un+1,v)+12(un·∇un,v)=(u(·,tn+1/2)·∇u(·,tn+1/2),v)+12∫tn+1/2tn+1‍(tn+1-t)d2dt2(u(·,t)·∇u(·,t),v)dt+12∫tntn+1/2‍(t-tn)d2dt2(u(·,t)·∇u(·,t),v)dt.
Add/subtract (un+1/2·∇u(·,tn+1/2),v) and apply (4.33) to get
(4.34)(un+1/2·∇un+1/2,v)-12(un+1·∇un+1,v)-12(un·∇un,v)=(un+1/2·∇(un+1/2-u(·,tn+1/2)),v)+((un+1/2-u(·,tn+1/2))·∇u(·,tn+1/2),v)=-12∫tn+1/2tn+1‍(tn+1-t)∫‍(utt·∇u+u·∇utt+2ut·∇ut)·vdt=-12∫tntn+1/2‍(t-tn)∫‍(utt·∇u+u·∇utt+2ut·∇ut)·vdt.
Majorize either directly or with (2.7)(a) to get
(4.35)|(un+1/2·∇un+1/2,v)-12(un+1·∇un+1,v)-12(un·∇un,v)|≤C∥∇u∥l∞(L2(Ω))(|un+1|1+|un|1)|vh|1,
or with (2.7)(b), (2.7)(c) and Hölder’s inequality (in time) applied to (4.34) to get
(4.36)|(un+1/2·∇un+1/2,v)-12(un+1·∇un+1,v)-12(un·∇un,v)|≤C∥u∥l∞(H2(Ω))∥un+1/2-u(·,tn+1/2)∥|vh|1≤+CΔt3/2tn+1/2∥u∥l∞(n,n+1;H2(Ω))(∫tntn+1‍t∥utt(·,t)∥2dt)1/2|vh|1≤+CΔt3/2tn+1/2∥ut∥l∞(n,n+1;L2)(∫tntn+1‍t∥ut(·,t)∥22dt)1/2|vh|1.
Then, to prove Corollary 4.7, apply (4.32), (4.36) with Young’s inequality give
(4.37)|τn(1)(u,p;ϕn+1/2h)|≤2εRe|ϕn+1/2h|12+CRe∥un+1-unΔt-(ut)n+1/2∥-12+CReΔt3∥u∥l∞(H2)2∥utt∥L2(tn,tn+1;L2(Ω))2+CReΔt3∥ut∥l∞(L2(Ω))2∥ut∥L2(tn,tn+1;H2(Ω))2.
We apply (4.35) instead of (4.36) to prove Theorem 4.6:
(4.38)|τn(1)(u,p;ϕn+1/2h)|≤2εRe|ϕn+1/2h|12+CRe∥un+1-unΔt-(ut)n+1/2∥-12+CRe∥∇u∥l∞(L2(Ω))2(|un+1|12+|un|12).
(Deconvolution error): Next, add/subtract φ(un+1/2) we write
(4.39)τn(2)(u,p;vh)=-bh(un+1/2-φ(un+1/2),un+1/2,ϕn+1/2h)-bh((φ(un+1/2)-φh(un+1/2)),un+1/2,ϕn+1/2h)+χ(un+1/2-φ(un+1/2),ϕn+1/2h)+χ(φ(un+1/2)-φh(un+1/2),ϕn+1/2h).
Then duality on X×X′ and Young’s inequalities give
(4.40)|χ(un+1/2-φ(un+1/2),ϕn+1/2h)+χ(φ(un+1/2)-φh(un+1/2),ϕn+1/2h)|≤Cχ2Re∥un+1/2-φ(un+1/2)∥-12≤+Cχ2Re∥φ(un+1/2)-φh(un+1/2)∥-12+1εRe|ϕn+1/2h|12,
and u∈H2(Ω) along with (2.11)(b) give
(4.41)|bh(un+1/2-φ(un+1/2),un+1/2,ϕn+1/2h)+bh(φ(un+1/2)-φh(un+1/2),un+1/2,ϕn+1/2h)|≤CRe∥un+1/2∥22(∥un+1/2-φ(un+1/2)∥2+∥φ(un+1/2)-φh(un+1/2)∥2)+1εRe|ϕn+1/2h|12.
The estimates (4.37), (4.40), (4.41) with identity (4.39) give
(4.42)|τn(1)(u,p;ϕn+1/2h)+τn(2)(u,p;ϕn+1/2h)|≤4εRe|ϕn+1/2h|12+CReEn.
Then estimates (3.17) and (3.23) give
(4.43)En:=(χ2+∥u∥l∞(H2(Ω))2)(c1,j(α,δ)2∥Δj+1un+1/2∥2+c2,j(h,α,δ)2∥un+1/2∥k+12)+∥un+1-unΔt-(ut)n+1/2∥-12+∥u∥l∞(H2)2∥un+1/2-u(·,tn+1/2)∥2+Δt3∥u∥l∞(H2)2∥utt∥L2(tn,tn+1;L2(Ω))2+Δt3∥ut∥l∞(L2(Ω))2∥ut∥L2(tn,tn+1;H2(Ω))2.
Apply estimates from (4.24), (4.25), (4.31), (4.37), (4.42) to (4.23). Set ε=20 and absorb all terms including |ϕn+1/2h|1 from the right into left-hand-side of (4.23). Sum the resulting inequality on both sides from n=0 to n=M-1 to get
(4.44)∥ϕMh∥2+Re-1Δt∑n=0M-1‍|ϕn+1/2h|12+χΔt∑n=0M-1‍∥ϕn+1/2h∥⋆h2≤∥ϕ0h∥2+CReΔt∑n=0M-1‍En+CReΔt∑n=0M-1‍∥pn+1/2-q~n+1/2h∥2≤+CRed12Δt∑n=0M-1‍(∥∇u∥L∞(L2(Ω))2+|ηn+1/2|12)|ηn+1/2|12≤+CReΔt∑n=0M-1‍(∥un+1/2∥22+h-1|ηn+1/2|12)∥ϕn+1/2h∥2≤+Cχ2ReΔt∑n=0M-1‍∥ηn+1/2∥-12+CRe-1Δt∑n=0M-1‍|ηn+1/2|12.
Estimates (2.23), (2.16)(a) imply
(4.45)supn|ηn|1≤C∥∇u∥l∞(L2(Ω)),|ηn|1≤Ch∥un∥2.
These estimates applied to (4.44) give
(4.46)∥ϕMh∥2+Re-1Δt∑n=0M-1‍|ϕn+1/2h|12+χΔt∑n=0M-1‍∥ϕn+1/2h∥⋆h2≤∥ϕ0h∥2+CReΔt∑n=0M-1‍En+CReΔt∑n=0M-1‍∥pn+1/2-q~n+1/2h∥2≤+CReΔt∑n=0M-1‍((d12∥∇u∥L∞(L2(Ω))2+Re-2)|ηn+1/2|12+∥un+1/2∥22∥ϕn+1/2h∥2)≤+Cχ2ReΔt∑n=0M-1‍∥ηn+1/2∥-12.
Suppose that the Δt-restriction (4.16) is satisfied. Then the discrete Gronwall Lemma 2.2 applies to (4.46) and gives
(4.47)∥ϕMh∥2+Re-1Δt∑n=0M-1‍|ϕn+1/2h|12+χΔt∑n=0M-1‍∥ϕn+1/2h∥⋆h2≤GM∥ϕ0h∥2+GMReΔt∑n=0M-1‍En≤+GMReΔt∑n=0M-1‍(∥pn+1/2-q~n+1/2h∥2+(d12∥∇u∥L∞(L2(Ω))2+Re-2)|ηn+1/2|12)≤+GMχ2ReΔt∑n=0M-1‍∥ηn+1/2∥-12,
where
(4.48)GM=Cexp(ReΔt∑n=0M‍gn∥un∥22),gn=(1-CReΔt∥un∥22)-1.
Lastly, the triangle inequality and approximation theory estimates (2.23), (2.16) along with (3.23) applied to (4.47) give
(4.49)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤GM∥e0h∥2+Ch2k(GM∥u0∥k+∥u∥l∞(Hk(Ω))2+Re-1∥u∥l2(Hk+1(Ω))2)≤+GMReΔt∑n=0M-1‍En+GMReh2s+2∥p∥l2(Hs+1(Ω))2≤+GMRe(d12∥∇u∥L∞(L2(Ω))2+Re-2+χ2h4)h2k∥u∥l2(Hk+1(Ω))2.
It remains to boundΔt∑n‍En.(Theorem4.6): Suppose that ∂tu∈L2(X′)∩l2(X′). The triangle inequality and (2.18) gives
(4.50)Δt∑n=n0M-1‍∥un+1-unΔt-(ut)n+1/2∥-12≤C(∥ut∥L2(X′)2+∥ut∥l2(X′)2).
Apply (4.38), instead of (4.37), to derive En in (4.43). Then
(4.51)ReΔt∑n=0M-1‍En≤E,
where
(4.52)E:=CRe(χ2+∥u∥l∞(H2(Ω))2)(c1,j(α,δ)2∥Δj+1un+1/2∥2+c2,j(h,α,δ)2∥un+1/2∥k+12)+CRe(∥ut∥L2(X′)2+∥ut∥l2(X′)2+∥∇u∥l∞(L2(Ω))2∥∇u∥l2(L2(Ω))2).
(Corollary 4.7): Suppose that u∈l∞(H2(Ω)), ut∈l∞(H1(Ω))∩L2(H1(Ω)), utt∈L2(L2(Ω)), and uttt∈L2(X′). Write
(4.53)E:=CRe(χ2+∥u∥l∞(H2(Ω))2)(c1,j(α,δ)2∥Δj+1un+1/2∥2+c2,j(h,α,δ)2∥un+1/2∥k+12):=CRe(∥uttt∥L2(X′)2+∥u∥l∞(H2(Ω))2∥utt∥L2(L2(Ω))2+⋯+∥ut∥l∞(L2(Ω))2∥ut∥L2(H2(Ω))2).
Then apply (2.19), (2.20) to bound En given in (4.43):
(4.54)ReΔt∑n=0M-1‍En≤EΔt4.


## 4.1. Well Posedness

We establish existence ofw at each time step of (4.3) by Leray-Schauder’s fixed-point theorem.Lemma 4.1.
Let(4.4)a(θh,vh)=Δt2Re(∇θh,∇vh)+χΔt2(θh-φh(θh),vh),ly(vh)=(y,vh).
for any y∈X′ and θh, vh∈Vh. Suppose that Dh satisfies Assumption 3.5. Then a(·,·):Vh×Vh→ℝ is a continuous and coercive bilinear form and ly(·):Vh→ℝ is a linear, continuous functional.Proof.
Linearity forly(·) is obvious, and continuity follows from an application of Hölder’s inequality. Continuity for a(·,·) also follows from Hölder’s inequality and Assumption 3.5. Coercivity is proven by application of (3.13).Lemma 4.2.
LetT:X′→Vh be such that, for any y∈X′, θh:=T(y) solves
(4.5)a(θh,vh)=ly(vh),∀vh∈Vh.
Then T is a well-defined, linear, bounded operator.Proof.
Linearity is clear. The results of Lemma4.1, and the Lax-Milgram theorem prove the rest.Lemma 4.3.
Fixn=0,1,…,M-1. Let wnh be a solution of Problem 1 and let N:Vh→X′ satisfy, for any θh∈Vh,
(4.6)(N(θh),vh)=-(θh-2wnh,vh)-Δt4bh(φh(θh),θh,vh)+Δt(fn+1/2,vh)=:c(θh,vh),∀vh∈Vh.
Then N(θh) is well-defined, bounded, and continuous.Proof.
For eachθh∈Vh, the map vh∈Vh↦c(θh,vh) is a bounded, linear functional (apply Hölder’s inequality and (2.11)). Since Vh is a Hilbert space, we conclude that N(θh) is well defined, by the Riesz-Representation theorem. Moreover, N(θh) is bounded on Vh and since the underlying function space is finite dimensional, continuity follows.Lemma 4.4.
Fixn∈ℕ. Let F:Vh→Vh be defined such that F(θh)=(T∘N)(θh). Then, F is a compact operator.Proof.
N
(
·
) is a compact operator (continuous on a finite dimensional function space). Thus, F is a continuous composition of a compact operator and hence compact itself.Theorem 4.5 (well posedness).
Fixn=0,1,2,…,M-1<∞. There exists (wnh,πnh)∈Xh×Qh satisfying Problem 1. Moreover,
(4.7)∥wmh∥2+12ReΔt∑n=0m-1‍|wn+1/2h|12+χΔt∑n=0m-1‍∥wn+1/2h∥⋆h2≤∥w0h∥2+ΔtRe2∑n=0m-1‍∥fn+1/2∥-12,
for all integers 1≤m≤M, independent of Δt>0.Proof.
First, assume that(wn+1h,qn+1h) is a solution to (4.1), (4.2). Set vh=wn+1/2h in (4.1) so that skew-symmetry of the nonlinear term provides
(4.8)12Δt(∥wn+1h∥2-∥wmh∥2)+Re-1|wn+1/2h|12+χ(wn+1/2h-φh(wn+1/2h),wn+1/2h)=(fn+1/2,wn+1/2h).
Duality of X′×X with Young’s inequality implies
(4.9)(fn+1/2,wn+1/2h)≤Re2∥fn+1/2∥-12+Re2|wn+1/2h|12.
From (3.13), we have
(4.10)∥wn+1/2h∥⋆h2=(wn+1/2h-φh(wn+1/2h),wn+1/2h)≥0.
Then applying (4.9), (4.10) to (4.8), combining-like terms and simplifying provides
(4.11)12Δt(∥wn+1h∥2-∥wnh∥2)+Re2|wn+1/2h|12+χ∥wn+1/2h∥⋆h2≤Re2∥fn+1/2∥-12.
Summing from n=0 to m-1, we get the desired bound.
Next, letWnh=wn+1h+wnh. Showing that Wnh=F(Wnh) has a fixed point will ensure existence of solutions to (4.3). Indeed, if we can show that W0h=F(W0h), then since w0h is given initial data, existence of w1h is immediate. Induction can be applied to prove existence of (wnh)1≤n≤M. To this end, since F is compact, it is enough to show (via Leray Schauder) that any solution Wn,λh of the fixed-point problem Wn,λh=λF(Wn,λh) is uniformly bounded with respect to 0≤λ≤1. Hence, we consider
(4.12)a(Wn,λh,vh)=λ(N(Wn,λh),vh).
Test withvh=Wn,λh, use skew-symmetry of the trilinear form and properties of Dh given in Assumption 3.5 and (3.13) to get
(4.13)λ∥Wn,λh∥2+Δt2Re|Wn,λh|12+χΔt2∥Wn,λh∥⋆h2≤2λ(wnh,Wn,λh)+λΔt(fn+1/2,Wn,λh).
Duality of X′×X followed by Young’s inequality implies
(4.14)λΔt(fn+1/2,Wn,λh)≤ΔtRe∥fn+1/2∥-12+Δt4Re|Wn,λh|12.
Since wnh∈L2(Ω) from the a priori estimate (4.7), we apply Hölder’s and Young’s inequalities to get
(4.15)2λ(wnh,Wn,λh)≤2∥wnh∥2+λ2∥Wn,λh∥2.
Applying estimates (4.14), (4.15) to (4.13) we get that |Wn,λh|1≤C<∞ independent of λ. By the Leray-Schauder fixed point theorem, given wnh, there exists a solution to the fixed-point theorem Wnh=F(Wnh). By the induction argument noted above, there exists a solution wnh for each n=0,1,2,…,M-1 to (4.3). Existence of an associated discrete pressure follows by a classical argument, since the pair (Xh,Qh) satisfies the discrete inf-sup condition (2.15).


## 4.2. Convergence Analysis

Under usual regularity assumptions, we summarize the main convergence estimate in Theorem4.6. Suppose that D represents deconvolution with J-updates.Theorem 4.6.
Suppose that(u,p) are strong solutions to (2.3), (2.4), (2.5) and that G, Gh, D, Dh satisfy Assumptions 3.4, 3.5. Suppose further that u∈l2(H2(Ω)∩V)∩l∞(X), ut∈l2(X′)∩L2(X′), Δβ+1u∈L2(Ω) for some 0≤β≤J+1, p∈l2(Q). If
(4.16)CReΔt∥un∥22<1,∀n=0,1,…,M
then,
(4.17)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤C(∥u0-w0h∥2+E+∥p∥l2(L2(Ω))2+(∥∇u∥l∞(L2(Ω))2+χ2)∥u∥l2(H1(Ω))2+∥u∥l∞(L2(Ω))2+∥u∥l2(L2(Ω))2+(χ2+∥u∥l∞(H2(Ω))2)c1,J2∥Δβu∥2+c2,j2∥∇u∥l2(L2(Ω))2),
where E>0 is given in (4.52).Corollary 4.7 (convergence estimate).
Under the assumptions of Theorem4.6, suppose further that (u,p) satisfy the assumptions for (2.16) for some k≥1 and s≥0, u∈l∞(Hk(Ω))∩l2(Hk+1(Ω)), ut∈L2(Hk-1(Ω)∩H1(Ω))∩l∞(H1(Ω)), utt∈L2(L2(Ω)), uttt∈L2(X′), and p∈l2(Hs+1(Ω)). If ||u0-w0h||≤C(hk+(αδ2)β) then
(4.18)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤C(h2k+h2s+2+Δt4+c1,J(δ,α)2+c2,J(h,δ,α)2).Proof of Theorem4.6, Corollary 4.7.
Suppose thatu, p satisfying (2.3), (2.4) also satisfy u∈C0(V), ut∈C0(X′), and p∈C0(Q) so that, for each n=n0,n0+1,…,M-1,
(4.19)((ut)n+1/2,v)+12(un+1·∇un+1,v)+12(un·∇un,v)+Re-1(∇un+1/2,∇v)-(pn+1/2,∇·v)=(fn+1/2,v),∀v∈X,∇·un+1=0.
The consistency error for the time-discretization τn(1)(u,p;vh) and regularization/time-relaxation error τn(2)(u,p;vh) are given by, for n=0,1,…,M-1,
(4.20)τn(1)(u,p;vh):=(un+1-unΔt-(ut)n+1/2,vh)+(un+1/2·∇un+1/2,vh)-12(un+1·∇un+1,vh)-12(un·∇un,v),τn(2)(u,p;vh):=-bh(un+1/2-φh(un+1/2),un+1/2,vh)+χ(un+1/2-φh(un+1/2),vh),
where vh∈Xh. Write τn:=τn(1)+τn(2). Using (4.20), rewrite (4.19) in a form conducive to analyzing the error between the continuous and discrete models:
(4.21)(un+1-unΔt,vh)+bh(φh(un+1/2),un+1/2,vh)+Re-1(∇un+1/2,∇vh)-(pn+1/2,∇·vh)+χ(un+1/2-φh(un+1/2),vh)=(fn+1/2,vh)+τn(u,p;vh).
Let v~h=Unh be the L2-projection of u(·,tn) so that (ηn+1-ηn,ϕn+1/2h)=0. Decompose the velocity error
(4.22)en=wnh-un=ϕnh-ηn,ϕnh=wnh-Unh,ηn=un-Unh,
where Unh∈Vh. Fix q~n+1/2h∈Qh. Note that (q~n+1/2h,∇·vh)=0 for any vh∈Vh. Subtract (4.21) from (4.3), apply (3.13)(b), and test with vh=ϕn+1/2h to get
(4.23)12Δt(∥ϕn+1h∥2-∥ϕnh∥2)+Re-1|ϕn+1/2h|12+χ∥ϕn+1/2h∥⋆h2=-(q~n+1/2h-pn+1/2,∇·ϕn+1/2h)-bh(φh(ϕn+1/2h),un+1/2,ϕn+1/2h)=+bh(φh(ηn+1/2),un+1/2,ϕn+1/2h)+bh(φh(wn+1/2h),ηn+1/2,ϕn+1/2h)=+Re-1(∇ηn+1/2,∇ϕn+1/2h)+χ(ηn+1/2-φh(ηn+1/2),ϕn+1/2h)-τn(u,p;ϕn+1/2h).
(Spatial discretization error): Fix ε>0. First, apply Hölder’s and Young’s inequalities to get
(4.24)|Re-1(∇ηn+1/2,∇ϕn+1/2h)+(pn+1/2-q~n+1/2h,∇·ϕn+1/2h)|≤CRe-1|ηn+1/2|12+CRe∥pn+1/2-q~n+1/2h∥2+1εRe|ϕn+1/2h|12.
Apply (3.12) and duality estimate on X×X′ to get
(4.25)|χ(ηn+1/2-φh(ηn+1/2),ϕn+1/2h)|≤Cχ2Re∥ηn+1/2∥-12+1εRe|ϕn+1/2h|1.
We bound the convective terms next. First, u∈H2(Ω) and estimate (2.11)(b) give
(4.26)|bh(φh(ϕn+1/2h),un+1/2,ϕn+1/2h)|≤CRe∥φh(ϕn+1/2h)∥2∥un+1/2∥22+1εRe|ϕn+1/2h|12
and u∈L∞(X) with (2.11)(a) give
(4.27)|bh(φh(ηn+1/2),un+1/2,ϕn+1/2h)|≤CRe∥∇u∥L∞(L2(Ω))2|φh(ηn+1/2)|12+1εRe|ϕn+1/2h|12.
Next, rewrite the remaining nonlinear term
(4.28)bh(φh(wn+1/2h),ηn+1/2,ϕn+1/2h)=bh(φh(un+1/2),ηn+1/2,ϕn+1/2h)-bh(φh(ηn+1/2),ηn+1/2,ϕn+1/2h)+bh(φh(ϕn+1/2h),ηn+1/2,ϕn+1/2h).
Once again, u∈L∞(X) and (2.11)(a) give
(4.29)|bh(φh(un+1/2),ηn+1/2,ϕn+1/2h)+bh(φh(ηn+1/2),ηn+1/2,ϕn+1/2h)|≤CRe∥∇φh(u)∥L∞(L2(Ω))∥φh(u)∥L∞(L2(Ω))|ηn+1/2|12≤+CRe|φh(ηn+1/2)|1|φh(ηn+1/2)|1|ηn+1/2|12+1εRe|ϕn+1/2h|12.
Lastly, estimate (2.11)(a) and inverse inequality (2.17) give (4.30)|bh(φh(ϕn+1/2h),ηn+1/2,ϕn+1/2h)|≤CReh-1∥φh(ϕn+1/2h)∥2|ηn+1/2|12+1εRe|ϕn+1/2h|12.
Apply estimates ||φh(v)||≤||v|| and |φh(v)|1≤d1|v|1 from Assumption 3.5 along with estimates (4.26)–(4.30) to get
(4.31)|bh(φh(ϕn+1/2h),un+1/2,ϕn+1/2h)+bh(φh(ηn+1/2),un+1/2,ϕn+1/2h)+bh(φh(wn+1/2h),ηn+1/2,ϕn+1/2h)|≤4εRe|ϕn+1/2h|12+CRe(∥un+1/2∥22+h-1|ηn+1/2|12)∥ϕn+1/2h∥2≤+Cd1Re(∥∇u∥L∞(L2(Ω))2+|ηn+1/2|12)|ηn+1/2|12.
(Time discretization error): First, apply duality estimate on X×X′ to get
(4.32)|(un+1-unΔt-(ut)n+1/2,ϕn+1/2h)|≤CRe∥un+1-unΔt-(ut)n+1/2∥-12+1εRe|ϕn+1/2h|12.
Taylor-expansion about tn+1/2 with integral remainder gives
(4.33)12(un+1·∇un+1,v)+12(un·∇un,v)=(u(·,tn+1/2)·∇u(·,tn+1/2),v)+12∫tn+1/2tn+1‍(tn+1-t)d2dt2(u(·,t)·∇u(·,t),v)dt+12∫tntn+1/2‍(t-tn)d2dt2(u(·,t)·∇u(·,t),v)dt.
Add/subtract (un+1/2·∇u(·,tn+1/2),v) and apply (4.33) to get
(4.34)(un+1/2·∇un+1/2,v)-12(un+1·∇un+1,v)-12(un·∇un,v)=(un+1/2·∇(un+1/2-u(·,tn+1/2)),v)+((un+1/2-u(·,tn+1/2))·∇u(·,tn+1/2),v)=-12∫tn+1/2tn+1‍(tn+1-t)∫‍(utt·∇u+u·∇utt+2ut·∇ut)·vdt=-12∫tntn+1/2‍(t-tn)∫‍(utt·∇u+u·∇utt+2ut·∇ut)·vdt.
Majorize either directly or with (2.7)(a) to get
(4.35)|(un+1/2·∇un+1/2,v)-12(un+1·∇un+1,v)-12(un·∇un,v)|≤C∥∇u∥l∞(L2(Ω))(|un+1|1+|un|1)|vh|1,
or with (2.7)(b), (2.7)(c) and Hölder’s inequality (in time) applied to (4.34) to get
(4.36)|(un+1/2·∇un+1/2,v)-12(un+1·∇un+1,v)-12(un·∇un,v)|≤C∥u∥l∞(H2(Ω))∥un+1/2-u(·,tn+1/2)∥|vh|1≤+CΔt3/2tn+1/2∥u∥l∞(n,n+1;H2(Ω))(∫tntn+1‍t∥utt(·,t)∥2dt)1/2|vh|1≤+CΔt3/2tn+1/2∥ut∥l∞(n,n+1;L2)(∫tntn+1‍t∥ut(·,t)∥22dt)1/2|vh|1.
Then, to prove Corollary 4.7, apply (4.32), (4.36) with Young’s inequality give
(4.37)|τn(1)(u,p;ϕn+1/2h)|≤2εRe|ϕn+1/2h|12+CRe∥un+1-unΔt-(ut)n+1/2∥-12+CReΔt3∥u∥l∞(H2)2∥utt∥L2(tn,tn+1;L2(Ω))2+CReΔt3∥ut∥l∞(L2(Ω))2∥ut∥L2(tn,tn+1;H2(Ω))2.
We apply (4.35) instead of (4.36) to prove Theorem 4.6:
(4.38)|τn(1)(u,p;ϕn+1/2h)|≤2εRe|ϕn+1/2h|12+CRe∥un+1-unΔt-(ut)n+1/2∥-12+CRe∥∇u∥l∞(L2(Ω))2(|un+1|12+|un|12).
(Deconvolution error): Next, add/subtract φ(un+1/2) we write
(4.39)τn(2)(u,p;vh)=-bh(un+1/2-φ(un+1/2),un+1/2,ϕn+1/2h)-bh((φ(un+1/2)-φh(un+1/2)),un+1/2,ϕn+1/2h)+χ(un+1/2-φ(un+1/2),ϕn+1/2h)+χ(φ(un+1/2)-φh(un+1/2),ϕn+1/2h).
Then duality on X×X′ and Young’s inequalities give
(4.40)|χ(un+1/2-φ(un+1/2),ϕn+1/2h)+χ(φ(un+1/2)-φh(un+1/2),ϕn+1/2h)|≤Cχ2Re∥un+1/2-φ(un+1/2)∥-12≤+Cχ2Re∥φ(un+1/2)-φh(un+1/2)∥-12+1εRe|ϕn+1/2h|12,
and u∈H2(Ω) along with (2.11)(b) give
(4.41)|bh(un+1/2-φ(un+1/2),un+1/2,ϕn+1/2h)+bh(φ(un+1/2)-φh(un+1/2),un+1/2,ϕn+1/2h)|≤CRe∥un+1/2∥22(∥un+1/2-φ(un+1/2)∥2+∥φ(un+1/2)-φh(un+1/2)∥2)+1εRe|ϕn+1/2h|12.
The estimates (4.37), (4.40), (4.41) with identity (4.39) give
(4.42)|τn(1)(u,p;ϕn+1/2h)+τn(2)(u,p;ϕn+1/2h)|≤4εRe|ϕn+1/2h|12+CReEn.
Then estimates (3.17) and (3.23) give
(4.43)En:=(χ2+∥u∥l∞(H2(Ω))2)(c1,j(α,δ)2∥Δj+1un+1/2∥2+c2,j(h,α,δ)2∥un+1/2∥k+12)+∥un+1-unΔt-(ut)n+1/2∥-12+∥u∥l∞(H2)2∥un+1/2-u(·,tn+1/2)∥2+Δt3∥u∥l∞(H2)2∥utt∥L2(tn,tn+1;L2(Ω))2+Δt3∥ut∥l∞(L2(Ω))2∥ut∥L2(tn,tn+1;H2(Ω))2.
Apply estimates from (4.24), (4.25), (4.31), (4.37), (4.42) to (4.23). Set ε=20 and absorb all terms including |ϕn+1/2h|1 from the right into left-hand-side of (4.23). Sum the resulting inequality on both sides from n=0 to n=M-1 to get
(4.44)∥ϕMh∥2+Re-1Δt∑n=0M-1‍|ϕn+1/2h|12+χΔt∑n=0M-1‍∥ϕn+1/2h∥⋆h2≤∥ϕ0h∥2+CReΔt∑n=0M-1‍En+CReΔt∑n=0M-1‍∥pn+1/2-q~n+1/2h∥2≤+CRed12Δt∑n=0M-1‍(∥∇u∥L∞(L2(Ω))2+|ηn+1/2|12)|ηn+1/2|12≤+CReΔt∑n=0M-1‍(∥un+1/2∥22+h-1|ηn+1/2|12)∥ϕn+1/2h∥2≤+Cχ2ReΔt∑n=0M-1‍∥ηn+1/2∥-12+CRe-1Δt∑n=0M-1‍|ηn+1/2|12.
Estimates (2.23), (2.16)(a) imply
(4.45)supn|ηn|1≤C∥∇u∥l∞(L2(Ω)),|ηn|1≤Ch∥un∥2.
These estimates applied to (4.44) give
(4.46)∥ϕMh∥2+Re-1Δt∑n=0M-1‍|ϕn+1/2h|12+χΔt∑n=0M-1‍∥ϕn+1/2h∥⋆h2≤∥ϕ0h∥2+CReΔt∑n=0M-1‍En+CReΔt∑n=0M-1‍∥pn+1/2-q~n+1/2h∥2≤+CReΔt∑n=0M-1‍((d12∥∇u∥L∞(L2(Ω))2+Re-2)|ηn+1/2|12+∥un+1/2∥22∥ϕn+1/2h∥2)≤+Cχ2ReΔt∑n=0M-1‍∥ηn+1/2∥-12.
Suppose that the Δt-restriction (4.16) is satisfied. Then the discrete Gronwall Lemma 2.2 applies to (4.46) and gives
(4.47)∥ϕMh∥2+Re-1Δt∑n=0M-1‍|ϕn+1/2h|12+χΔt∑n=0M-1‍∥ϕn+1/2h∥⋆h2≤GM∥ϕ0h∥2+GMReΔt∑n=0M-1‍En≤+GMReΔt∑n=0M-1‍(∥pn+1/2-q~n+1/2h∥2+(d12∥∇u∥L∞(L2(Ω))2+Re-2)|ηn+1/2|12)≤+GMχ2ReΔt∑n=0M-1‍∥ηn+1/2∥-12,
where
(4.48)GM=Cexp(ReΔt∑n=0M‍gn∥un∥22),gn=(1-CReΔt∥un∥22)-1.
Lastly, the triangle inequality and approximation theory estimates (2.23), (2.16) along with (3.23) applied to (4.47) give
(4.49)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤GM∥e0h∥2+Ch2k(GM∥u0∥k+∥u∥l∞(Hk(Ω))2+Re-1∥u∥l2(Hk+1(Ω))2)≤+GMReΔt∑n=0M-1‍En+GMReh2s+2∥p∥l2(Hs+1(Ω))2≤+GMRe(d12∥∇u∥L∞(L2(Ω))2+Re-2+χ2h4)h2k∥u∥l2(Hk+1(Ω))2.
It remains to boundΔt∑n‍En.(Theorem4.6): Suppose that ∂tu∈L2(X′)∩l2(X′). The triangle inequality and (2.18) gives
(4.50)Δt∑n=n0M-1‍∥un+1-unΔt-(ut)n+1/2∥-12≤C(∥ut∥L2(X′)2+∥ut∥l2(X′)2).
Apply (4.38), instead of (4.37), to derive En in (4.43). Then
(4.51)ReΔt∑n=0M-1‍En≤E,
where
(4.52)E:=CRe(χ2+∥u∥l∞(H2(Ω))2)(c1,j(α,δ)2∥Δj+1un+1/2∥2+c2,j(h,α,δ)2∥un+1/2∥k+12)+CRe(∥ut∥L2(X′)2+∥ut∥l2(X′)2+∥∇u∥l∞(L2(Ω))2∥∇u∥l2(L2(Ω))2).
(Corollary 4.7): Suppose that u∈l∞(H2(Ω)), ut∈l∞(H1(Ω))∩L2(H1(Ω)), utt∈L2(L2(Ω)), and uttt∈L2(X′). Write
(4.53)E:=CRe(χ2+∥u∥l∞(H2(Ω))2)(c1,j(α,δ)2∥Δj+1un+1/2∥2+c2,j(h,α,δ)2∥un+1/2∥k+12):=CRe(∥uttt∥L2(X′)2+∥u∥l∞(H2(Ω))2∥utt∥L2(L2(Ω))2+⋯+∥ut∥l∞(L2(Ω))2∥ut∥L2(H2(Ω))2).
Then apply (2.19), (2.20) to bound En given in (4.43):
(4.54)ReΔt∑n=0M-1‍En≤EΔt4.


## 5. Applications

We show that the iterated (modified) Tikhonov regularization operator satisfied Assumption3.4, 3.5 in Section 5.1 and verify the theoretical convergence rate predicted by Theorem 4.6, Corollary 4.7 in Section 5.2.


### 5.1. Iterated (Modified) Tikhonov-Lavrentiev Regularization

We will prove thatDα,J, Dα,Jh (Definitions 3.14, 3.15) with the differential filter G=A-1 satisfies Assumptions 3.4, 3.5. Proposition 3.10 implies that it is enough to show that Dα,0 satisfies Assumption 3.4. Additionally, we provide sharpened estimates for d1,j, c1,j, c2,j. The key is that A-1>0 is a continuous function of the Laplace operator -Δ≥0 and hence they commute (on X). Moreover, Dα,0>0 is a continuous function of A-1 so that Dα,0 commutes with A-1 and -Δ (on X).We first characterize the spectrum ofDα,0, Dα,0h.Lemma 5.1.
Fix0<α≤1. Define f:(0,1]→ℝ and g:(0,1]→ℝ by
(5.1)f(x):=1(1-α)x+α,g(x):=x(1-α)x+α.
The maps f and g are continuous and f((0,1])=[1,α-1) and g((0,1])=(0,1].Proof.
The functionsf, g are clearly continuous with f  decreasing and g increasing on (0,1]. Hence, the range of f is [1,α-1) and range of g is (0,1].The next result shows thatDα,0 and Dα,0h satisfy part of Assumptions 3.4, 3.5.Proposition 5.2.
D
α
,
0 and Dα,0h (on X and Xh, resp.) are linear, bounded, spd, and commute with A-1, (Ah)-1 (resp.). Moreover,
(5.2)∥Dα,0A-1∥≤1,|Dα,0u¯|1≤|u|1∀u∈X,∥Dα,0h(Ah)-1∥≤1,|Dα,0huh¯h|1≤|uh|1∀uh∈Xh.
Hence d1=1 in Assumptions 3.4, 3.5.Proof.
It is immediately clear thatDα,0, Dα,0h are linear. As a consequence, since A-1>0 with spectrum in (0,1], then Dα,0=f(A-1) with spectrum contained in [1,α-1) so that D>0. Therefore, Dα,0A-1=g(A-1) with spectrum contained in (0,1]. A similar argument shows that (Ah)-1 has spectrum in (0,1], Dα,0h has spectrum in [1,α-1), and Dα,0h(Ah)-1 has spectrum in (0,1]. Thus Dα,0>0, Dα,0h>0 and ||Dα,0A-1||≤1 and ||Dα,0h(Ah)-1||≤1. Therefore, Dα,0 and Dα,0h are bounded and commute with A-1 and (Ah)-1, respectively, as discussed above.
The second set of inequalities on each line can be proved with an appropriate choice ofv and vh in Definitions 3.2 and 3.3. Starting with Definition 3.2, take ϕ=u and choose v=ΔDα,0u¯. Then integration by parts and the Cauchy-Schwartz inequality give the result. The discrete form is proved using Definition 3.3 and choosing ϕ=uh and v=ΔDα,0huh¯h.It remains to provide estimates forc1 and c2, and sharpened estimates for c1,j and c2,j. Indeed, as a direct consequence of Propositions 5.3, 5.4, we have, for each j=0,1,…,J,
(5.3)c1,j=(αδ2)j+1∥Δj+1v∥,∀v∈Hj+1(Ω),c2,j=C(h+(2jαδ2)1/2)hkmax0≤n≤j∥Dα,nv¯∥k+1,∀v∈Hk+1(Ω).Proposition 5.3.
Letj=0,1,…,J. Then, for some 0≤β≤j+1,
(5.4)∥v-Dα,jv¯∥≤(αδ2)β∥Δβv∥,∀v∈Hβ(Ω).Proof.
Using (1.9), we have
(5.5)Dα,0-1(Dα,Jv¯-Dα,J-1v¯)=v¯-A-1Dα,J-1v¯.
Subtracting (5.5) from the identity
(5.6)Dα,0-1(v-v)=v¯-v¯,
gives us
(5.7)Dα,0-1((I-Dα,JA-1)v-(I-Dα,J-1A-1)v)=-A-1(I-Dα,J-1A-1)v.
Multiplying by Dα,0, rearranging, simplifying, and using A-I=-δ2Δ (Definition 3.2) gives
(5.8)(I-Dα,JA-1)v=[-Dα,0A-1(I-Dα,J-1A-1)+(I-Dα,J-1A-1)]v=(I-Dα,J-1A-1)(I-Dα,0A-1)v=(I-Dα,J-1A-1)Dα,0A-1(Dα,0-1A-I)v=(I-Dα,J-1A-1)Dα,0A-1(((1-α)A-1+αI)A-I)v=(I-Dα,J-1A-1)Dα,0A-1α(A-I)v=-αδ2ΔDα,0A-1(I-Dα,J-1A-1)v.
Applying recursion, we obtain, for any 0≤β<J,(5.9)(I-Dα,JA-1)v=(-αδ2)β(Dα,J-βA-1)βΔβv.
Thus, taking norms and applying ||Dα,J-βA-1||≤1, we get (5.4).Proposition 5.4.
Letj=0,1,…,J. Then
(5.10)∥Dα,jw¯-Dα,jhw¯h∥2≤C(h2+2j+1αδ2)h2kmax0≤n≤j∥Dα,nw¯∥k+12,∀w∈Hk+1(Ω).Proof.
Letv~jh∈Xh be the L2-projection of Dα,jw¯. Take v=vh in (3.34). For j=1,…,J, let ej=Dα,jw¯-Dα,jhw¯h:=ηj-ϕjh, where ηj:=Dα,jw¯-v~jh, and ϕjh:=Dα,jhw¯h-v~jh. Subtract (3.34) and (3.35) to get
(5.11)αδ2(∇ϕjh,∇vh)+(ϕjh,vh)=αδ2(∇ηj,∇vh)+αδ2(∇ej-1,∇vh).
Take vh=ϕjh in (5.11) to get
(5.12)αδ2|ϕjh|12+∥ϕjh∥2=αδ2(∇ηj,∇ϕjh)+αδ2(∇ej-1,∇ϕjh).
Fix ε>0. Apply Hölder’s and Young’s inequalities to (5.12) to get
(5.13)αδ2|ϕjh|12+∥ϕjh∥2≤αδ2|ηj|12+εαδ2|ej-1|12+1εαδ2|ϕjh|12.
Taking ε=1 and ε=2 in (5.13) gives
(5.14)∥ϕjh∥2≤αδ2|ηj|12+αδ2|ej-1|12.(5.15)αδ2|ϕjh|12+2∥ϕjh∥2≤2αδ2|ηj|12+4αδ2|ej-1|12.
The triangle inequality and estimate (5.14) give
(5.16)∥ej∥2≤∥ηj∥2+αδ2|ηj|12+αδ2|ej-1|12.
Backward induction, estimate (5.15), and (2.23) give
(5.17)∥ej∥2≤|e0|12+(h2+αδ2(1+∑i=0j‍2i))max0≤n≤jinfvh∈Xh|Dα,nw¯-vh|12.
It has been shown (Estimate (2.36) in the proof of Lemma 2.7 [25]) that
(5.18)|e0|12≤C(h2+αδ2)h2k|D0w¯|k+12.
Note that ∑i=0j‍2i=2j+1-1. Then, along with application of (2.16), we prove (5.10).Corollary 5.5 (convergence estimate).
Under the assumptions of Corollary4.7, suppose further that, for some J=0,1,…, that G=A-1, Gh=(Ah)-1, D=Dα,J, Dh=Dα,Jh. If Δβu∈L2(Ω) for some 0≤β≤J+1, then
(5.19)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤C((1+h2+2Jαδ2)h2k+h2s+2+Δt4+(αδ2)2β).Proof.
Apply estimates forc1,j, c2,j from (5.3), resulting from Propositions 5.3, 5.4.


### 5.2. Numerical Testing

This section presents the calculation of a flow with an exact solution to verify the convergence rates of the algorithm. FreeFEM++ [32] was used to run the simulations. The convergence rates are tested against the Taylor-Green vortex problem [13, 33–35]. We use a domain of Ω=(0,1)×(0,1) and take u=(u1,u2), where
(5.20)u1(x,y,t)=-cos(nπx)sin(nπy)e-2n2π2t/τ,u2(x,y,t)=sin(nπx)cos(nπy)e-2n2π2t/τ,p(x,y,t)=-14(cos(nπx)+cos(nπy))e-2n2π2t/τ.
The pair (u,p) is a solution the two-dimensional NSE when τ=Re and f=0.We used CN discretization in time and P2-P1 elements in space according to Problem1. That is, we used continuous piecewise quadratic elements for the velocity and continuous piecewise linear elements for the pressure. We chose the spatial discretization elements and parameters n=1, T=0.5, χ=0.1 and Re=10,000 as a illustrative example. We chose h=1/m, dt=(1/4)h, δ=h4 and α=h, where m is the number of mesh divisions per side of [0,1]. These were chosen so that the result of Corollary 5.5 reduces to
(5.21)∥uM-wMh∥+[Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12]1/2≤C(h2+hJ+1).We summarize the results in Tables1 and 2. Table 1 displays error estimates corresponding to no iterations; that is, J=0 in Definition 3.14. For the particular choice of α and δ, the computed errors ∥u-wh∥∞,0 and ∥∇(u-wh)∥2,0 tend to the predicted convergence rate 𝒪(h). Table 2 displays error estimates corresponding to one update; that is, when J=1 in Definition 3.14. Again, for the particular choice of α and δ,the computed errors ∥u-wh∥∞,0 and ∥∇(u-wh)∥2,0 tend to the predicted convergence rate 𝒪(h2).Table 1
Error and convergence rates for Leray-deconvolution withJ=0 for the Taylor-Green vortex with Re=10,000, α=h, and δ=h4. Note the convergence rate is approaching 1 as predicted by (5.21).
m (=1/h)
∥
u
-
w
h
∥
∞
,
0
Rate
∥
∇
(
u
-
w
h
)
∥
2,0
Rate
20
0.038975
1.651230
40
0.024334
0.680
1.468510
0.169
60
0.017751
0.778
1.159840
0.582
80
0.013854
0.862
0.935247
0.748
100
0.011255
0.931
0.774285
0.846Table 2
Error and convergence rates for Leray-deconvolution withJ=1 for the Taylor-Green vortex with Re=10,000, α=h, and δ=h4. Note the convergence rate is approaching 2 as predicted by (5.21).
m (=1/h)
∥
u
-
w
h
∥
∞
,
0
Rate
∥
∇
(
u
-
w
h
)
∥
2,0
Rate
20
0.023384
1.070400
40
0.009739
1.264
0.640360
0.741
60
0.004997
1.646
0.357779
1.436
80
0.002899
1.892
0.212560
1.810
100
0.001915
1.858
0.136724
1.977


## 5.1. Iterated (Modified) Tikhonov-Lavrentiev Regularization

We will prove thatDα,J, Dα,Jh (Definitions 3.14, 3.15) with the differential filter G=A-1 satisfies Assumptions 3.4, 3.5. Proposition 3.10 implies that it is enough to show that Dα,0 satisfies Assumption 3.4. Additionally, we provide sharpened estimates for d1,j, c1,j, c2,j. The key is that A-1>0 is a continuous function of the Laplace operator -Δ≥0 and hence they commute (on X). Moreover, Dα,0>0 is a continuous function of A-1 so that Dα,0 commutes with A-1 and -Δ (on X).We first characterize the spectrum ofDα,0, Dα,0h.Lemma 5.1.
Fix0<α≤1. Define f:(0,1]→ℝ and g:(0,1]→ℝ by
(5.1)f(x):=1(1-α)x+α,g(x):=x(1-α)x+α.
The maps f and g are continuous and f((0,1])=[1,α-1) and g((0,1])=(0,1].Proof.
The functionsf, g are clearly continuous with f  decreasing and g increasing on (0,1]. Hence, the range of f is [1,α-1) and range of g is (0,1].The next result shows thatDα,0 and Dα,0h satisfy part of Assumptions 3.4, 3.5.Proposition 5.2.
D
α
,
0 and Dα,0h (on X and Xh, resp.) are linear, bounded, spd, and commute with A-1, (Ah)-1 (resp.). Moreover,
(5.2)∥Dα,0A-1∥≤1,|Dα,0u¯|1≤|u|1∀u∈X,∥Dα,0h(Ah)-1∥≤1,|Dα,0huh¯h|1≤|uh|1∀uh∈Xh.
Hence d1=1 in Assumptions 3.4, 3.5.Proof.
It is immediately clear thatDα,0, Dα,0h are linear. As a consequence, since A-1>0 with spectrum in (0,1], then Dα,0=f(A-1) with spectrum contained in [1,α-1) so that D>0. Therefore, Dα,0A-1=g(A-1) with spectrum contained in (0,1]. A similar argument shows that (Ah)-1 has spectrum in (0,1], Dα,0h has spectrum in [1,α-1), and Dα,0h(Ah)-1 has spectrum in (0,1]. Thus Dα,0>0, Dα,0h>0 and ||Dα,0A-1||≤1 and ||Dα,0h(Ah)-1||≤1. Therefore, Dα,0 and Dα,0h are bounded and commute with A-1 and (Ah)-1, respectively, as discussed above.
The second set of inequalities on each line can be proved with an appropriate choice ofv and vh in Definitions 3.2 and 3.3. Starting with Definition 3.2, take ϕ=u and choose v=ΔDα,0u¯. Then integration by parts and the Cauchy-Schwartz inequality give the result. The discrete form is proved using Definition 3.3 and choosing ϕ=uh and v=ΔDα,0huh¯h.It remains to provide estimates forc1 and c2, and sharpened estimates for c1,j and c2,j. Indeed, as a direct consequence of Propositions 5.3, 5.4, we have, for each j=0,1,…,J,
(5.3)c1,j=(αδ2)j+1∥Δj+1v∥,∀v∈Hj+1(Ω),c2,j=C(h+(2jαδ2)1/2)hkmax0≤n≤j∥Dα,nv¯∥k+1,∀v∈Hk+1(Ω).Proposition 5.3.
Letj=0,1,…,J. Then, for some 0≤β≤j+1,
(5.4)∥v-Dα,jv¯∥≤(αδ2)β∥Δβv∥,∀v∈Hβ(Ω).Proof.
Using (1.9), we have
(5.5)Dα,0-1(Dα,Jv¯-Dα,J-1v¯)=v¯-A-1Dα,J-1v¯.
Subtracting (5.5) from the identity
(5.6)Dα,0-1(v-v)=v¯-v¯,
gives us
(5.7)Dα,0-1((I-Dα,JA-1)v-(I-Dα,J-1A-1)v)=-A-1(I-Dα,J-1A-1)v.
Multiplying by Dα,0, rearranging, simplifying, and using A-I=-δ2Δ (Definition 3.2) gives
(5.8)(I-Dα,JA-1)v=[-Dα,0A-1(I-Dα,J-1A-1)+(I-Dα,J-1A-1)]v=(I-Dα,J-1A-1)(I-Dα,0A-1)v=(I-Dα,J-1A-1)Dα,0A-1(Dα,0-1A-I)v=(I-Dα,J-1A-1)Dα,0A-1(((1-α)A-1+αI)A-I)v=(I-Dα,J-1A-1)Dα,0A-1α(A-I)v=-αδ2ΔDα,0A-1(I-Dα,J-1A-1)v.
Applying recursion, we obtain, for any 0≤β<J,(5.9)(I-Dα,JA-1)v=(-αδ2)β(Dα,J-βA-1)βΔβv.
Thus, taking norms and applying ||Dα,J-βA-1||≤1, we get (5.4).Proposition 5.4.
Letj=0,1,…,J. Then
(5.10)∥Dα,jw¯-Dα,jhw¯h∥2≤C(h2+2j+1αδ2)h2kmax0≤n≤j∥Dα,nw¯∥k+12,∀w∈Hk+1(Ω).Proof.
Letv~jh∈Xh be the L2-projection of Dα,jw¯. Take v=vh in (3.34). For j=1,…,J, let ej=Dα,jw¯-Dα,jhw¯h:=ηj-ϕjh, where ηj:=Dα,jw¯-v~jh, and ϕjh:=Dα,jhw¯h-v~jh. Subtract (3.34) and (3.35) to get
(5.11)αδ2(∇ϕjh,∇vh)+(ϕjh,vh)=αδ2(∇ηj,∇vh)+αδ2(∇ej-1,∇vh).
Take vh=ϕjh in (5.11) to get
(5.12)αδ2|ϕjh|12+∥ϕjh∥2=αδ2(∇ηj,∇ϕjh)+αδ2(∇ej-1,∇ϕjh).
Fix ε>0. Apply Hölder’s and Young’s inequalities to (5.12) to get
(5.13)αδ2|ϕjh|12+∥ϕjh∥2≤αδ2|ηj|12+εαδ2|ej-1|12+1εαδ2|ϕjh|12.
Taking ε=1 and ε=2 in (5.13) gives
(5.14)∥ϕjh∥2≤αδ2|ηj|12+αδ2|ej-1|12.(5.15)αδ2|ϕjh|12+2∥ϕjh∥2≤2αδ2|ηj|12+4αδ2|ej-1|12.
The triangle inequality and estimate (5.14) give
(5.16)∥ej∥2≤∥ηj∥2+αδ2|ηj|12+αδ2|ej-1|12.
Backward induction, estimate (5.15), and (2.23) give
(5.17)∥ej∥2≤|e0|12+(h2+αδ2(1+∑i=0j‍2i))max0≤n≤jinfvh∈Xh|Dα,nw¯-vh|12.
It has been shown (Estimate (2.36) in the proof of Lemma 2.7 [25]) that
(5.18)|e0|12≤C(h2+αδ2)h2k|D0w¯|k+12.
Note that ∑i=0j‍2i=2j+1-1. Then, along with application of (2.16), we prove (5.10).Corollary 5.5 (convergence estimate).
Under the assumptions of Corollary4.7, suppose further that, for some J=0,1,…, that G=A-1, Gh=(Ah)-1, D=Dα,J, Dh=Dα,Jh. If Δβu∈L2(Ω) for some 0≤β≤J+1, then
(5.19)∥uM-wMh∥2+Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12≤C((1+h2+2Jαδ2)h2k+h2s+2+Δt4+(αδ2)2β).Proof.
Apply estimates forc1,j, c2,j from (5.3), resulting from Propositions 5.3, 5.4.


## 5.2. Numerical Testing

This section presents the calculation of a flow with an exact solution to verify the convergence rates of the algorithm. FreeFEM++ [32] was used to run the simulations. The convergence rates are tested against the Taylor-Green vortex problem [13, 33–35]. We use a domain of Ω=(0,1)×(0,1) and take u=(u1,u2), where
(5.20)u1(x,y,t)=-cos(nπx)sin(nπy)e-2n2π2t/τ,u2(x,y,t)=sin(nπx)cos(nπy)e-2n2π2t/τ,p(x,y,t)=-14(cos(nπx)+cos(nπy))e-2n2π2t/τ.
The pair (u,p) is a solution the two-dimensional NSE when τ=Re and f=0.We used CN discretization in time and P2-P1 elements in space according to Problem1. That is, we used continuous piecewise quadratic elements for the velocity and continuous piecewise linear elements for the pressure. We chose the spatial discretization elements and parameters n=1, T=0.5, χ=0.1 and Re=10,000 as a illustrative example. We chose h=1/m, dt=(1/4)h, δ=h4 and α=h, where m is the number of mesh divisions per side of [0,1]. These were chosen so that the result of Corollary 5.5 reduces to
(5.21)∥uM-wMh∥+[Re-1Δt∑n=0M-1‍|un+1/2-wn+1/2h|12]1/2≤C(h2+hJ+1).We summarize the results in Tables1 and 2. Table 1 displays error estimates corresponding to no iterations; that is, J=0 in Definition 3.14. For the particular choice of α and δ, the computed errors ∥u-wh∥∞,0 and ∥∇(u-wh)∥2,0 tend to the predicted convergence rate 𝒪(h). Table 2 displays error estimates corresponding to one update; that is, when J=1 in Definition 3.14. Again, for the particular choice of α and δ,the computed errors ∥u-wh∥∞,0 and ∥∇(u-wh)∥2,0 tend to the predicted convergence rate 𝒪(h2).Table 1
Error and convergence rates for Leray-deconvolution withJ=0 for the Taylor-Green vortex with Re=10,000, α=h, and δ=h4. Note the convergence rate is approaching 1 as predicted by (5.21).
m (=1/h)
∥
u
-
w
h
∥
∞
,
0
Rate
∥
∇
(
u
-
w
h
)
∥
2,0
Rate
20
0.038975
1.651230
40
0.024334
0.680
1.468510
0.169
60
0.017751
0.778
1.159840
0.582
80
0.013854
0.862
0.935247
0.748
100
0.011255
0.931
0.774285
0.846Table 2
Error and convergence rates for Leray-deconvolution withJ=1 for the Taylor-Green vortex with Re=10,000, α=h, and δ=h4. Note the convergence rate is approaching 2 as predicted by (5.21).
m (=1/h)
∥
u
-
w
h
∥
∞
,
0
Rate
∥
∇
(
u
-
w
h
)
∥
2,0
Rate
20
0.023384
1.070400
40
0.009739
1.264
0.640360
0.741
60
0.004997
1.646
0.357779
1.436
80
0.002899
1.892
0.212560
1.810
100
0.001915
1.858
0.136724
1.977


## 6. Conclusion

It is infeasible to resolve all persistent and energetically significant scales down to the Kolmogorov microscale of𝒪(Re-3/4) for turbulent flows in complex domains using direct numerical simulations in a given time constraint. Regularization methods are used to find approximations to the solution. The modification of iterated Tikhonov-Lavrentiev to the modified iterated Tikhonov-Lavrentiev deconvolution in Definition 3.14 is a highly accurate method of solving the deconvolution problem in the Leray-deconvolution model, with errors u-Dα,0u¯=𝒪((αδ2)J+1) when applied to the differential filter. We use this result to show that under a regularity assumption, the error between the solutions to the NSE and to the Leray deconvolution model with time relaxation using the modified iterated Tikhonov-Lavrentiev deconvolution and discretized with CN in time and FE’s in space are 𝒪(hk(h+αδ2)+hs+1+Δt2+(αδ2)J+1).We also examined the Taylor-Green vortex problem using Problem1 with the deconvolution in Definition 3.14. We use this problem because it has an exact analytic solution to the NSE. The regularization parameters α and δ were chosen so that the convergence of the approximate solution to the error would be 𝒪(hJ+1) for J=0 and J=1. The convergence rates calculated correspond to those predicted, that is 𝒪(h1) for J=0 and 𝒪(h2) for J=1.



---

*Source: 162539-2012-09-27.xml*