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

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 TikhonovLavrentiev 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.


Approximate Deconvolution for Turbulence Modeling
Numerical simulations of complex flows present many challenges. The Navier-Stokes NS equations NSE , given by the following Advances in Numerical Analysis 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: time relaxation w t w · ∇w − Re −1 Δw ∇p χ w − w f, 1.5 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 replacing w 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: find w : Ω × 0, T → R 3 and π : Ω × 0, T → R satisfying the following w t D w · ∇w − Re −1 Δw ∇π χ w − D w f, ∇ · w 0, 1.6 Advances in Numerical Analysis 3 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 .

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 ,

Advances in Numerical Analysis
In fact, the updates {ω j } j≥0 satisfying 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 w h of the general deconvolution turbulence model Problem 1 satisfies Theorem 4.6

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.

Function Spaces and Approximations
Let the flow domain Ω ∈ R 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 L 2 -norm and inner product, respectively.

2.1
Denote the pressure and velocity spaces by Q : 2 Ω and equipped with the norm Fix f ∈ X and Re > 0. In this setting, we consider strong NS solutions: find u ∈ Let V : {v ∈ X : ∇ · v 0}. Restricting test functions v ∈ V reduces 2.3 -2.5 to find u : 0, T → V satisfying u t , v u · ∇u, v Re −1 ∇u, ∇v f, v , a.e. t ∈ 0, T , ∀v ∈ V 2.6 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:

Discrete Function Setting
Fix h > 0. Let T h be a family of subdivisions e.g., triangulation of Ω ⊂ R d satisfying Ω E∈T h E so that diameter E ≤ h and any two closed elements E 1 , E 2 ∈ T 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, T h consists of triangles for d 2 or tetrahedra for d 3 that are nondegenerate as h → 0.

Advances in Numerical Analysis
Let X h ⊂ X and Q h ⊂ Q be a conforming velocity-pressure mixed FE space. For example, let X h and Q h be continuous, piecewise on each E ∈ T h polynomial spaces. The discretely divergence-free space is given by Note that in general V h / ⊂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 V h / ⊂V , we introduce the explicitly skew-symmetric convective term b h u, v, w : Note that b h u, v, w u · ∇v, w when u ∈ V . Moreover, the trilinear from b h ·, ·, · is continuous and skew-symmetric on X × X × X.
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, let 0 t 0 < t 1 < · · · < t M−1 T < ∞ be a discretization of the time interval 0, T for a constant time step Δt t n 1 −t n . Write t n 1/2 : t n 1 t n /2, z n z t n and, if z ∈ C 0 t n , t n 1 , z n 1/2 z n 1 z n /2. Define u l q m 1 ,m 2 ;W k p Ω : for any 0 ≤ n m 1 , m 1 1, . . . , m 2 ≤ M. Write ||u|| l q W k p Ω ||u|| l q 0,M;W k p Ω . We say that u ∈ l q m 1 , m 2 ; W k p Ω if the associated norm defined above stays finite as Δt → 0. The discrete Gronwall inequality is essential to the convergence analysis in Section 4.2. Suppose that for all n, Δtκ n < 1 and set g n 1 − Δtκ n −1 . Then, Proof. The proof follows from 26 .

Approximation Theory
Let C > 0 be a generic constant independent of h → 0 . Preserving an abstract framework for the FE spaces, we assume that X h ×Q h inherit several fundamental approximation properties.
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 timedependent problems: for any n 0, 1, . . . , M − 1, where θ ∈ H 1 H k Ω , θ ∈ H 2 H k Ω , and θ ∈ H 3 H k Ω is required, respectively, for some k ≥ −1. Each estimate 2.18 -2.20 is a result of a Taylor expansion with integral remainder.

Advances in Numerical Analysis
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 p 0 is the solution of the well-posed Neumann problem 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 : Replacing 2.21 with 2.21 a , 2.22 defines an overdetermined Neumann-type problem. Condition 2.22 is a nonlocal condition relating u 0 and f 0 . Condition 2.22 is satisfied for several practical applications including start up from rest with zero force, u 0 0, f 0 0. In general, however, condition 2.22 cannot be verified. In this case, it is shown that, for We finish with an approximation property of the L 2 -projection. Indeed, Assumption 2.4 holds for smooth enough Ω.

Filters and Deconvolution
We prescribe the essential properties our filter G and deconvolution operator D in this section.
Let G G δ > 0 be a linear, bounded, compact operator on X representing a generic smoothing filter with filter radius δ > 0: Advances in Numerical Analysis 9 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 A h −1 Definition 3.3 , for implementation of our numerical scheme Section 5.2 .
It is well known that A −1 and A h −1 are each linear and bounded, A −1 is compact, and the spectrum of A and A h on X and X h , resp. is contained in 1, ∞ and spectrum of A −1 and A h −1 on X and X h , resp. is contained in 0, 1 so that For more detailed exposition on these operators, see 13 .

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 that D : X → X is linear, bounded, spd, and commutes with G so that for some constant d 1 > 0. Moreover, suppose that D is parametrized by α > 0, δ > 0 so that 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 H 1 -seminorm of DG is required for the convergence analysis in Section 4.2.
Assumption 3.5 prescribes properties of the discrete analogue D h : X h → X h corresponding to the continuous deconvolution operator D : X → X Assumptions 3.4 .
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.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 D h G h . We next derive several important consequences of D and D h under Assumptions 3.4, 3.5 required in the forthcoming analysis.
Proof. For the continuous operator, Advances in Numerical Analysis 11 Then, 3.10 a follows from Assumption 3.4, and 3.10 b is derived similarly applying Assumption 3.5 instead.
As a consequence, Proof. Assumptions 3.4, 3.5 guarantee that the spectral radius ρ DG ≤ 1 and ρ D h G h ≤ 1. Also, D > 0 and G > 0 and commute so that have spectrum contained in 0, 1 which ensures the non-negativity of both || · || and || · || h .

Iterated Deconvolution
One can show, by eliminating intermediate steps in the definition of the iterated regularization operator D j in 1.9 with base operator D satisfying Assumption 3.4, that Similarly, the discrete iterated regularization operator D h j with discrete base operator D h satisfying Assumption 3.5, is given by the following We next show that D j and D h j for j > 0 inherit several important properties from D and D h , respectively, via Assumption 3.5. Proposition 3.10. Fix j ∈ N. Then D j : X → X defined by 3.14 satisfies Assumption 3.4. In particular, D j > 0 is linear, bounded, commutes with G and satisfies 3.6 (a). Estimate 3.6 (b) is replaced by the following Advances in Numerical Analysis for some constant d 1,j > 0. Estimate 3.7 is replaced by the following Proof. First notice that D j 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 D j . Next, D j is a sum of spd and snn operators D > 0, D I − DG i ≥ 0. Hence, D j > 0. Next, notice that Letting λ k · denote the kth ordered eigenvalue of a given operator, we can characterize the spectrum of D j by summing the resulting finite geometric series 3.18 to get Then under Assumption 3.4, Lemma 3.9 with 3.19 implies that 0 ≤ λ k D j G ≤ ||D j G|| ≤ 1. Hence, D j 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 a i in DG, so that since |Dv| 1 ≤ d 1 |v| 1 can be applied successfully. Therefore 3.16 follows with d 1,j j i 0 d i 1 . Next, start with 3.18 to get for any v ∈ X ∩ H k 1 Ω for some k ≥ 0. Moreover, c 2,j ≤ β j c 2 for some constant β β j > 0.
Proof. The first two assertions follow similarly as in the previous proof of Proposition 3.10. To prove 3.23 , we start by writing and then subtract 3.24 from 3.18 to get Then taking norms across 3.25 , we get Notice that ||I − DG|| ≤ 1 so that ||Λ j || ≤ j 1. Moreover, || DG − D h G h v|| ≤ c 2 ||v|| k 1 via Assumption 3.5. Next, using the binomial theorem and factoring, we get

Tikhonov-Lavrentiev Regularization
We provide two examples of discrete deconvolution operators D h 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
Definition 3.12 weak modified Tikhonov-Lavrentiev deconvolution . Fix α > 0. Let G A −1 . For any w ∈ X, let ω 0 : D α,0 w ∈ X be the unique solution of The iterated modified Tikhonov-Lavrentiev operator for linear, compact G > 0 is obtained from the Tikhonov-Lavrentiev operator with updates via 1.9 :

Well Posedness of the Fully Discrete Model
We now state the proposed algorithm.

Well Posedness
We establish existence of w at each time step of 4.3 by Leray-Schauder's fixed-point theorem.

4.4
for any y ∈ X and θ h , v h ∈ V h . Suppose that D h satisfies Assumption 3.5. Then a ·, · : V h × V h → R is a continuous and coercive bilinear form and l y · : V h → R is a linear, continuous functional.
Proof. Linearity for l y · 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.
Let T : X → V h be such that, for any y ∈ X , θ h : T y solves Then T is a well-defined, linear, bounded operator.

Advances in Numerical Analysis
Proof. Linearity is clear. The results of Lemma 4.1, and the Lax-Milgram theorem prove the rest.

4.6
Then N θ h is well-defined, bounded, and continuous.

linear functional apply
Hölder's inequality and 2.11 . Since V h is a Hilbert space, we conclude that N θ h is well defined, by the Riesz-Representation theorem. Moreover, N θ h is bounded on V h and since the underlying function space is finite dimensional, continuity follows.
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.

Convergence Analysis
Under usual regularity assumptions, we summarize the main convergence estimate in Theorem 4.6. Suppose that D represents deconvolution with J-updates.  Suppose that u, p are strong solutions to 2.3 , 2.4 , 2.5 and that G, G h , D, D then,

4.17
where E > 0 is given in 4.52 .

4.19
Advances in Numerical Analysis

19
The consistency error for the time-discretization τ 1 n u, p; v h and regularization/timerelaxation error τ 2 n u, p; v h are given by, for n 0, 1, . . . , M − 1,

4.20
where v h ∈ X h . Write τ n : τ 1 n τ 2 n . Using 4.20 , rewrite 4.19 in a form conducive to analyzing the error between the continuous and discrete models:

4.21
Let v h U h n be the L 2 -projection of u ·, t n so that η n 1 − η n , φ h

20
Advances in Numerical Analysis Apply 3.12 and duality estimate on X × X to get

4.25
We bound the convective terms next. First, u ∈ H 2 Ω and estimate 2.

4.27
Next, rewrite the remaining nonlinear term

4.31
Time discretization error : First, apply duality estimate on X × X to get

Advances in Numerical Analysis
Majorize either directly or with 2.7 a to get

4.35
or with 2.7 b , 2.7 c and Hölder's inequality in time applied to 4.34 to get

4.48
Lastly, the triangle inequality and approximation theory estimates 2.23 , 2.16 along with 3.23 applied to 4.47 give

4.49
It remains to bound Δt n E n . Theorem 4.6 : Suppose that ∂ t u ∈ L 2 X ∩ l 2 X . The triangle inequality and 2.18 gives

Applications
We show that the iterated modified Tikhonov regularization operator satisfied Assumption 3.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.

Iterated (Modified) Tikhonov-Lavrentiev Regularization
We will prove that D α,J , D h α,J 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 d 1,j , c 1,j , c 2,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 of D α,0 , D h α,0 .
Proof. The functions f, 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 that D α,0 and D h α,0 satisfy part of Assumptions 3.4, 3.5.

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 30 Advances in Numerical Analysis  rates are tested against the Taylor-Green vortex problem 13, 33-35 . We use a domain of Ω 0, 1 × 0, 1 and take u u 1 , u 2 , where u 1 x, y, t − cos nπx sin nπy e −2n 2 π 2 t/τ , u 2 x, y, t sin nπx cos nπy e −2n 2 π 2 t/τ , p x, y, t − 1 4 cos nπx cos nπy e −2n 2 π 2 t/τ .

5.20
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 Problem 1. 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

5.21
We summarize the results in Tables 1 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 − w h ∞,0 and ∇ u − w h 2,0 tend to the predicted convergence rate O 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 − w h ∞,0 and ∇ u − w h 2,0 tend to the predicted convergence rate O h 2 .
Advances in Numerical Analysis 31

Conclusion
It is infeasible to resolve all persistent and energetically significant scales down to the Kolmogorov microscale of O 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 α,0 u O αδ 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 O h k h √ αδ 2 h s 1 Δt 2 αδ 2 J 1 . We also examined the Taylor-Green vortex problem using Problem 1 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 O h J 1 for J 0 and J 1. The convergence rates calculated correspond to those predicted, that is O h 1 for J 0 and O h 2 for J 1.