Improving Accuracy in α-Models of Turbulence through Approximate Deconvolution

In this report, we present several results in the theory of α -models of turbulence with improved accuracy that have been developed in recent years. The α -models considered herein are the Leray- α model, the zeroth Approximate Deconvolution Model (ADM) turbulence model, the modified Leray- α and the Navier–Stokes- α model. For all of the models from above, the accuracy is limited to α 2 in smooth flow regions. Better accuracy requires decreasing the filter radius α , which, in turn, requires a smaller mesh width that will lead in the end to a higher computational cost. Instead, one can use approximate deconvolution (without decreasing the mesh size) to attain better accuracy. Such deconvolution methods have been considered recently in many studies that show the efficiency of this approach. For smooth flows, periodic boundary conditions and van Cittert deconvolution operator of order N, the expected accuracy is α 2 N + 2 . In a bounded domain, such results are valid only in case special conditions are satisfied. In more general conditions, the author has recently proved that, in the case of the ADM, the expected accuracy of the finite element method with Taylor–Hood elements and Crank–Nicolson time stepping method is Δ t 2 + h 2 + K N α 2 , where the constant K < 1 depends on the ratio α / h , which is assumed constant. In this study, we present the extension of the result to the rest of the models.


Introduction
At a high Reynolds number, turbulence is not efficient to simulate using the Navier-Stokes equations because they require a very fine mesh and a high computational cost [1].Therefore, in practical applications, reduced order turbulence models, such as the α-models of turbulence discussed herein, have to be used.
The α-models of turbulence are regularizations of the Navier-Stokes equations (NSE) whose goal is to allow stable computations on coarser grids than NSE.They are given by in Ω, w(x, 0) = w 0 (x) in Ω. (1) Here, the nonlinearity N is N(w) = w • ∇w in the Leray-α model [2], N(w) = w • ∇w in the case of the modified Leray-α model [3], N(w) = w • ∇w in the case of the zeroth Approximate Deconvolution Model (ADM) model [4] (here, the ADM is formulated as in [5]), and N(w) = (∇ × w) × w in the case of the Navier-Stokes-α model.
One property that all four models have in common is that their accuracy ||u NSE − w α−model || is limited to O(α 2 ) for smooth flows.In practice, this may lead to over-diffusivity and a lack of accuracy that requires smaller filter radius α and, as a consequence, a smaller mesh size that in turn causes higher computational costs.One way to correct this deficiency is to use an approximate deconvolution operator D having the approximation property u ≈ Du.One example of such deconvolution operator is the van Cittert deconvolution operator that approximates u by D N u defined as For N = 0, 1, 2, the velocity field u will be approximated by The finite element analysis of the van Cittert deconvolution method has been studied in [21].The α-models enhanced using approximate deconvolution become where the nonlinearity DN is DN(w) = Dw • ∇w in the Leray-deconvolution model, DN(w) = w • ∇Dw in the case of the modified Leray-deconvolution model, DN(w) = Dw • ∇Dw in the case of the ADM model (formulated as in [21]) and DN(w) = (∇ × w) × Dw in the case of the enhanced Navier-Stokes-α model.
Most of the analysis for the deconvolution-enhanced α-models has been carried out in the context of van Cittert deconvolution.For periodic boundary conditions and assuming a smooth NSE solution, the modeling error is where w is the solution of the deconvolution enhanced Leray-α or the ADM or the deconvolution enhanced Navier-Stokes-α model [5,11,19,24,34].Here, N is the order of the van Cittert deconvolution operator D N .In a bounded domain, such estimates are valid only if special boundary conditions are satisfied by the exact NSE solution u [35].
The finite element analysis has been carried out for the Leray-deconvolution and the deconvolution enhanced Navier-Stokes-α in [19,22] with Crank-Nicolson time discretization and P2/P1 elements and the main error estimate is in the case of periodic boundary conditions.Here, the filter radius α satisfies α = O(h).In the above formula, the term α 2N+2 is induced by the deconvolution modeling.The case N = 0 corresponds to the classical α turbulence models.In the case of a bounded domain, the formula above is no longer valid [35], and one can only prove that, for all N, the error satisfies Recently, in [29], the following estimate has been proved for the ADM: In the above estimate, N is the order of the van Cittert deconvolution operator and K < 1 is a constant that only depends on the ratio α/h (which is assumed constant) and the sequence of quasiuniform meshes used in the computation.The term K N α 2 is due to the deconvolution and it can be made small by increasing N and decreasing the ratio α/h [21].This estimate supports the claim that high order deconvolution operators indeed improve accuracy, a behavior that has been observed in the numerical tests in [22,29].
We believe that this is the general behavior of deconvolution enhanced α-models of turbulence in case the van Cittert deconvolution procedure is used.Using similar techniques such as those in [29], similar estimates can be obtained for the Leray-deconvolution model, the Navier-Stokes α deconvolution model and the modified Leray-α deconvolution model.

Mathematical Context
The mathematical notations and concepts are similar to the ones used in [19,22,29].Ω ⊂ R d , d = 2, 3, is regular, bounded, polyhedral domain and • and (•, •) is the L 2 (Ω) norm and inner product.Lebesgue spaces and their norms are denoted by L p (Ω), • L p .We use standard notations for the Lebesque and Sobolev spaces and their norms.
In the variational problem, we use the functional spaces [36]: The finite element analysis is carried out on a family of triangulations (T h ) h on Ω that are quasiuniform [37].We also consider an inf-sup stable pair of finite elements (X h , Q h ) where X h ⊂ X, Q h ⊂ Q, which are assumed to satisfy [37]: where Iu h ∈ X h is an interpolant of u and I p h ∈ Q h is an interpolant of p.
The discretely divergence-free space V h is We assume that V h satisfies the first two approximation properties in the inequalities (11) above if u ∈ V.
We let t n = n∆t, t n+1/2 = (n + 1/2)∆t, T := N T ∆t.Here, ∆t is the chosen time step, N T is the number of time steps, and T is the final time.We will also use the notations We will assume that the L 2 projection operator onto V h denoted by The analysis also use the inverse inequality [37]: The constant C T does not depend on h.Definition 1. [30] The discrete Stokes operator A h : One may easily show that A h is a bijective, self-adjoint, positive operator with respect to the inner products (•, The eigenvalues of A h are denoted by One may show using the inverse inequality ( 13) that for v ∈ V h .
Definition 2. (Discrete Stokes average [30]) For given α > 0, we let G : Theorem 1. [29] Due to the spectral mapping theorem, the eigenvalues of the discrete Stokes filter G are .., M, and they satisfy One may also notice that the largest eigenvalue of I − G is bounded by Definition 3. [30] The N-th order discrete van Cittert deconvolution operator D N : V h → V h is given by: Several properties of D N that can be found in [30] are listed below.D N is self-adjoint, positive with respect to the inner products (•, and can be bounded as [29]: Following [19], we define two norms on V h Using inequalities (20), it follows that, see also for more details [29], and We will also use the next estimate of the discrete deconvolution error v − D N v h that has been obtained in [21].
Here, N is the order of the discrete deconvolution operator and K < 1 has been previously defined in (17).
Proof.We recall the main ideas of the proof presented in [21].First, we split the error ||v − D N v h || using the triangle inequality The term ||v − P V h v|| is bounded by Ch k+1 |v| k+1 due to the approximation properties of V h .N+1 , this term is written as

It remains then to estimate ||P
It follows that The discrete Laplacian ||A h v h || can be estimated as The last inequality from above can be found in [38].It follows that This proves the first inequality in inequalities (24).The second inequality in inequalities ( 24) is immediatelly proved using the same arguments combined with inverse inequalities.

Numerical Scheme for the α-Models of Turbulence
In our analysis b * (•, •, •) : X × X × X → IR will denote the standard skew-symmetrized trilinear form [36]: One can show that the trilinear form satisfies [39]: It follows then using inequalities (35) and (36) in [38] that for u, w ∈ X, v ∈ V h .
In the case of the Navier-Stokes in rotational form, the trilinear form is For body force f ∈ L 2 ((0, T], X) and initial velocity u 0 ∈ X, the discrete solution w n+1 h ∈ V h of the accuracy enhanced α model of turbulence at step n + 1 for n = 0, 1, 2, ..., N T − 1 satisfies: where , v), Leray-deconvolution model [22], Thorough the analysis, we assume that α/h is constant.

Lemma 2. (Stability) A finite element solution w n+1
h of the model (34) exists at each time step n = 0, . . .N T − 1 and satisfies the stability estimate where C does not depend on h, α (but it depends on α/h).
Proof.The Leray-deconvolution model case has been studied (for a slightly different filter operator, but the argument is still valid) in [22].The NS-α-deconvolution case in [19].The modified Leray-deconvolution case can be proved using arguments similar to the NS α-deconvolution (since it requires a similar multiplier to cancel the nonlinearity).
The ADM case has been studied in [29] and shares great similarity with the stability proofs for the Leray-deconvolution and the NS-α-deconvolution models mentioned previously.We outline the ideas from [29] here.
To prove the existence of a discrete solution u n+1 h ∈ V h , we consider the operator L : V h → V h , Lϕ = ψ such that ϕ, ψ satisfy the equality We notice that, if ϕ is a fixed point of L, then u n+1 h = 2ϕ − u n h is a solution of the ADM model (34).In order to find a fixed point of ϕ, we apply the Leray-Schauder fixed point theorem.To show that L is a compact operator, we notice that it is the composition of a continuous, linear operator and a nonlinear, continuous and bounded operator N : To show that N is bounded and continuous, we use the equivalence of all norms on the finite element space V h , similar to the proof in [19].
It remains only to show that the set is bounded in the L 2 norm independent of λ.
and set v = ϕ, cancel the nonlinearity and get Next, using the standard Cauchy-Schwartz and Yound inequalities on the right terms, we get that Therefore, by the Leray-Schauder Theorem, it follows that the operator L has a fixed point ϕ = u n+1/2 h and the discrete solution u n+1 h of the ADM system exists.
To obtain the stability estimate from above, we set D N u n+1/2 h h as a test function in the ADM model (34) and the nonlinear term will vanish: Using the standard Cauchy-Schwartz and Young inequalities on the right term, we get Summing up from n = 0 to n = N T − 1, we get the required inequality Here, C = C(Ω).
The norms || • || E , || • || can be replaced by the standard L 2 , respectively H 1 norm via the inequality (22), and this will add a dependence of C on K, (see inequality 17), i.e., on the ratio α/h.
The convergence of the discrete solution w h of the enhanced α-model (34) to the NSE exact solution u is stated in the next theorem.Theorem 3. We let (u,p) be a smooth strong solution of the NSE and w n h , n = 1...N T be the discrete solution of (34).We assume α/h is constant and w 0 = u 0 .Then, if ∆t is picked small enough, there holds sup where C is a constant that does not depend on h, α (though it depends on α/h) and K is given in inequality (17) .
Proof.The convergence proof for the Leray-deconvolution model is done in [22] and for the NS-α deconvolution model in [19], but the two proofs do not use the improved deconvolution error estimate (24) and the error is proved to be in the bounded domain case on the order of In the two proofs in [19,22], one bounds by several terms which can be left unchanged herein and also by the modelling error ( the Intp term in [22], the Interp term in [19]), which itself contains several terms that can be left unchanged but also a term that is bounded by the deconvolution error ||u n+1/2 − D N u n+1/2 ||.Applying the improved deconvolution estimate (24) provides the improved K 2N α 4 instead of α 4 .The proof of the above result for the ADM is done in [29].We outline here the main ideas of the proof presented in [29], which are very similar to the ones in [22] except that the nonlinearity is handled differently.
Similar to [19,22], the NSE is written in the form where Next, one derives the error equation for e n+1 = u n+1 − u n+1 h by subtracting (34) from (50).We let U n+1 ∈ V h be the L 2 projection of u n+1 onto V h , split the error as in the error equation.Upon rearranging terms, we obtain On the right side, all above terms excepting the two trilinear terms appearing in the interpolating term (51) will be estimated exactly as in [22].In all of these terms, wherever needed (to get rid of D N G terms and switch to the norms || • || E , || • || ) to keep the arguments similar to the ones in [22], inequalities (23) will be used.
The two trilinear terms in Interp that are treated differently compared to [22] take the form The first term on the right side above is estimated similar to inequality 3.26 in [22]: For the second term, we first use the inequality (32) to get Using again the argument in inequality (28) to estimate the middle term and Young's inequality, we get further Therefore, both trilinear terms have been bounded by the deconvolution error and some terms that will eventually be hidden on the left side of Equation ( 52).The proof continues as in [22] by applying the discrete Gronwall's inequality.The application of the Gronwall's inequality conditions the time-step to satisfy ∆t ≤ C(ν −3 |||∇u||| ∞,0 + 1) −1 .
For the modified Leray-deconvolution model, the proof follows the steps outlined above down to equality (53) which, in the case of the ML-α, will take the form which is further estimated as in inequality (54) from above.From here on, the proof proceeds as in the ADM case.

A Numerical Experiment
In this section, we verify the theoretical convergence rates on a 2d NSE problem from [29] with exact solution u 1 (x, y, t) = sin(2πy)e −4νπ Herein, we check the Leray-deconvolution model with N = 0, 1 (thus suplementing the numerical results presented in [22]) and the modified Leray-deconvolution model with N = 0, 1.The ADM with N = 0, 1 has been checked in this problem in [29] on a slightly different mesh and is therefore omitted here.The NS-α deconvolution model has been checked numerically in [19].
The scope of the test is to verify the convergence rates and also to check numerically that, as N increases, the corresponding model enters the predicted convergent regime faster.
The test is carried out with the FreeFEM++ package [40].The computational domain is the unit square Ω = (0, 1) × (0, 1).The kinematic viscosity is set equal to ν = 1.0 and the final time T is set equal to T = 0.1.An initial mesh M 0 is generated with five even nodes on the edges x = 0, x = 1, 4 even nodes on y = 0 and 6 even nodes on y = 1.Then, the edges are succesivelly halved to generate the computational meshes M 1 , M 2 , . . ., M 5 , see Figure 1.For the space discretization, we use P2/P1 Taylor-Hood finite elements.In our computation, we set the filter radius α equal to the mesh size h, i.e., α = h.First, we solve the Leray-deconvolution and modified Leray-deconvolution with N = 0, 1 on M 1 with ∆t = 1/80 and eight iterations.When solving the next mesh M 2 , the time step ∆t is halved and we double the number of time iterations.We keep applying this procedure up to the last mesh M 5 .
To stabilize pressure, the model and the filter equation are augmented with an L 2 pressure stabilization term with parameter 10 −9 .The resulting algebraic system corresponding to the discrete model is solved using a fixed point iteration until the L ∞ norm of two successive iterates is less than 10 −12 .The rate presented in the Tables 1-4 is computed as the log 2 of the quotient of two successive errors.The theoretical rate in the L ∞ (L 2 ) and L 2 (H 1 ) norms for the Leray-deconvolution and modified Leray-deconvolution models is 2 and is confirmed by this numerical test.Moreover, the predicted convergent regime is reached faster for N = 1.

Conclusions
This report presents some results on the theory of α models of turbulence that have been obtained in recent years.They show that the accuracy of these models can be improved using higher order approximate deconvolution operators.Some numerical tests are also presented to support the theoretical results.The technique to estimate the deconvolution error can be used in the larger context of models using approximate deconvolution such as the time-relaxation model investigated in [41] or for other models, such as the MHD, that have been recently investigated using approximate deconvolution methods.

Table 3 .
L ∞ (L 2 ) and L 2 (H 1 ) errors and rates for the modified Leray-α model.The final time is T = 0.1.The predicted rate order is 2.