A posteriori error estimates for the Large Eddy Simulation applied to incompressible fluids

. We study the two dimensional time dependent Large Eddy Simulation method applied to the incompressible Navier–Stokes system with Smagorinsky’s eddy viscosity model and a filter width that depends on the local mesh size. The discrete model is based on the implicit Euler scheme and a conforming finite element method for the time and space discretizations, respectively. We establish a reliable and efficient a posteriori error estimation between the numerical LES solution and the exact solution of the original Navier–Stokes system, which involves three types of error indicators respectively related to the filter and to the discretizations in time and space. Numerical results show the effectiveness of adaptive simulations


Introduction.
Let Ω be a bounded connected open domain in IR 2 , with a Lipschitz-continuous connected boundary Γ = ∂Ω.Let [0, T ] denote an interval in IR where T is a positive constant.For a positive constant viscosity ν, we consider the following time-dependent Navier-Stokes system, in which we choose a vanishing initial condition for simplification only: (t, x) − ν∆u(t, x) + u(t, x) • ∇u(t, x) + ∇p(t, x) = f (t, x) in]0, T [×Ω, div u(t, x) = 0 in[0, T ] × Ω, where f represents a density of body forces and supposed to be in [L 2 ((0, T ) × Ω)] 2 (although a more general setting would be to consider f ∈ L 2 (0, T ; X ′ ) = L 2 (0, T ; H −1 (Ω) 2 ), where X ′ is the dual of the Sobolev space X = H 1 0 (Ω) 2 ).The unknowns are the velocity u and the pressure p of the fluid.Note that the hypothesis of vanishing initial condition is only here for the sake of simplicity and may easily be removed.
Since an exact solution of System (P ) is in general out of reach, one may resort to a numerical approximation of it.A first method of simulation in fluid mechanics is the direct numerical simulation (DNS) of the flow at all significant length scales.Since it is necessary to capture all fluctuations of the velocity and pressure fields having an impact on the flow, DNS is very expensive and, for high Reynolds numbers, is even not achievable.Under these circumstances, one often prefers the large eddy simulation (LES) method, which consists in solving the large scales and in modelling the influence of small scales by adding a supplementary non-linear diffusion term in the momentum equation.In LES, large scales are defined by a spatial average of velocity, pressure and external forcing terms.A common method is to define this spatial average by convolution of these quantities with an appropriate filter function of width denoted by δ.The velocity field will be decomposed as: u(t, x) = ū(t, x) + u ′ (t, x), with ū(t, x) the filtered part of the velocity field and u ′ (t, x) = u(t, x) − ū(t, x), the residue whose effect on the large scale motions is the main issue in LES.Indeed, the convolution of the non-linear term u • ∇u is not equal to the non-linear term applied to the convolution ū • ∇ū and the difference between these two terms needs to be modeled.There is vast literature on this issue and we refer to monographs [27,10,21] for a mathematically oriented review.The Smagorinksy model is the simplest model of LES that uses the assumption of local balance between production and dissipation of turbulent kinetic energy to express the turbulent viscosity according to the large scales (see [21]).The resulting filtered problem is expressed as: a 2 ij ) 1/2 for any matrix A. The value of the constant c s is an input of the model and a debated question among specialists (see, e.g., [10, p. 76]).Moreover, finding relevant boundary conditions for ū in bounded domains is a central issue in LES (see, e.g., [10,Part IV]); in ( P ), we have denoted by B a general boundary operator that we do not discuss here.The boundary condition issue is also closely related to that of filtering with non constant width δ, which is commonly used in practice (the value of δ is typically adjusted to the local size of the mesh cells in numerical simulations) but introduces further problems since there is now a commutation error between the operations of filtering and differentiating.
In view of all these questions, our approach here is not to further refine LES models, but rather to consider a discretization of the Smagorinksy model with turbulent viscosity adjusted to the local size of the mesh cells and to derive a posteriori error estimates between its solution and the exact solution of (P ), which will allow us to locally refine or coarsen the mesh and/or the time steps according to local error indicators.Adaptive LES is an active field of research in the quest for accurate simulations with an affordable computational cost.Among works using mesh adaptation in the context of the Smagorinksy model for incompressible flows, we may cite for example [1,3,19,26].However, these works are either based on local error metrics driven by the Hessian of the numerical solution, a process that is independent of the model and numerical method, or only mention a posteriori estimates, without indicating their actual expression, nor giving a proof of their derivation, reliability and efficiency.
The present work provides a rigorous theory based on a posteriori error estimates for adaptivity in LES simulations of incompressible flows.For a general introduction to a posteriori error estimation, we refer for example to the books of Verfürth [30] or Ainsworth and Oden [2].As far as time-dependent models are concerned, a large number of contributions may be found.To cite only a few of them, we can refer, for example, to Ladevèze [22] for constitutive relation error estimators for time-dependent non-linear FE analysis, Verfürth [31] for the heat equation, Bernardi and Verfürth [9] for the time dependent Stokes equations, Bernardi and Süli [8] for the time and space adaptivity for the second-order wave equation, Bergam, Bernardi and Mghazli [4] for some parabolic equations, Ern and Vohralïk [16] for estimation based on potential and flux reconstruction for the heat equation, and Bernardi and Sayah [7] for the time dependent Stokes equations with mixed boundary conditions.In [25], Nassreddine and Sayah treated the time dependent Navier-Stokes equations in two dimensions.They use [6,7] for the linear (Stokes) part of the problem and apply the continuous and discrete Gronwall lemma to treat the non-linear term.In this contribution, we extend the method of [25] to the LES problem, for which we shall obtain three types of computable error indicators; the first one is linked to the time discretization, the second one to the filter and the last one to the space discretization.This leads to an adaptive strategy in which we sometimes change the time step and sometimes adapt the mesh.
The outline of the paper is as follows: Section 2 is devoted to some functional analysis tools and to the variational formulation of the continuous Navier-Stokes system.In section 3, we introduce the discrete filtered problem and we recall its main properties.In sections 4 and 5 we study the a posteriori error indicators and derive quasi-optimal estimates, in the sense that we need to require higher regularity of the exact solution.In section 6, we present numerical results in which uniform and adaptive refinement strategies are applied and comparisons with a reference DNS simulation show the superiority of the LES adaptive strategy.

Preliminaries
In this section, we begin by some notations and definitions which will be used later on and we recall the continuous and discrete Gronwall Lemmas.We denote by [L p (Ω)] 2 the space of measurable functions v such that |v| p is integrable.For v ∈ [L p (Ω)] 2 , the norm is defined by Let α = (α 1 , α 2 ) be a couple of non negative integers and |α| = α 1 + α 2 .We define the partial derivative ∂ α by Then, for any positive integer m and number p ≥ 1, we recall the classical Sobolev space equipped with the following semi-norm and norm : When p = 2, this space is the Hilbert space [H m (Ω)] 2 .In view of the boundary conditions in system (P ), we thus consider the space We denote by M the space of functions in L 2 (Ω) with a zero mean-value on Ω.
Lemma 2.1.For any 1 ≤ p < +∞, there exists a positive constant S p only depending on Ω such that Lemma 2.2.(See [29, page 291]) We have the following inequality for every v ∈ X Remark 2.
, which then prevents from pursuing like in dimension 2.
We introduce the kernel of the divergence operator which is a closed subspace of X and coincides with Remark 2.4.The spaces M and X satisfy a uniform inf-sup condition (see [18]): There exists a constant As usual, for handling time-dependent problems, it is convenient to consider functions defined on a time interval ]a, b[ with values in a separable functional space, say Y .More precisely, let || • || Y denote the norm of Y ; then for any r, 1 ≤ r ≤ ∞, we define Definition 2.5.We introduce the trilinear form c defined by : Lemma 2.6.For every u, v, w ∈ X we have As far as Problem (P ) is concerned, we recall classical existence and uniqueness results that may be found, for example in [11,Th. V.1.4]and [29,Th. 3.1 and Th. 3.2] which require the definition , where γ n v is the normal trace operator on Γ.
Lemma 2.9.(Generalized Gronwall Lemma) [32, p. 292] and [15, p. 252].Let: (1) f, g and k be an integrable functions defined from 4) gk is an integrable function on IR + . then: Lemma 2.10.(Discrete Gronwall Lemma) [32, p. 294] let (y n ) n , (f n ) n and (g n ) n be three positive sequences that verify: Then, we have: Finally, we shall make an extensive use of the following inequalities, the first of which is Young's inequality with ε and is used three times to obtain the second: In particular, we shall use these inequalities as follows: For any (a, b, c, d) ∈ R 4 + , and any

The discrete problem
From now on, we assume that Ω is a polyhedron.In order to describe the time discretization with an adaptive choice of local time steps, we introduce a partition of the interval [0, T ] into subintervals , by τ the N-tuple (τ 1 , . . ., τ N ), by |τ | the maximum of the τ n , 1 ≤ n ≤ N , and finally by σ τ the regularity parameter In what follows, we work with a regular family of partitions, i.e. we assume that σ τ is bounded independently of τ .This has the practical implication that the time step should not be modified too abruptly in the adaptive algorithm described in Section 6.We now describe the space discretization.For each n, 0 ≤ n ≤ N , we consider (T nh ) h , a partition of Ω by triangles that belongs to a regular family of triangulations in the usual sense: • Ω is the union of all elements of T nh ; • the intersection of two different elements of T nh , if not empty, is a vertex or a whole edge of both of them; this, in particular, implies that meshes remain conforming (no hanging nodes) during the adaptive refinement process.• the ratio of the diameter h κ of an element κ in T nh to the diameter of its inscribed circle is bounded by a constant independent of n and h.
As usual, h denotes the maximal diameter of the elements of all T nh , 0 ≤ n ≤ N , while for each n, h n denotes the maximal diameter of the elements of T nh .For each κ in T nh and each non-negative integer k, we denote by P k (κ) the space of restrictions to κ of polynomials of 2 variables and total degree at most k.
In what follows, c, c ′ , C . . .stand for generic constants which may vary from line to line but are always independent of h and n.From now on, we call finite element space associated to T nh a space of functions such that their restrictions to any element κ of T nh belong to a space of polynomials of fixed degree.
For each n and h, we associate with T nh two finite element spaces X nh and M nh which are contained in X and M , respectively, and such that the following inf-sup condition holds for a constant β > 0, which is independent of n and h, There exist many examples of finite element spaces satisfying this condition.We give an example involving continuous discrete pressures.Velocities are discretized with the "Mini-Element" where the space P b (κ) is spanned by functions in P 1 (κ) and the bubble function on κ (for each element κ, the bubble function is equal to the product of the barycentric coordinates associated with the vertices of κ).The pressure is discretized with classical continuous finite elements of order one As usual, we denote by V nh the kernel of the divergence As discussed in the Introduction, the fully discrete problem which we consider is a variational version of ( P ) with a local choice of the filter size, in which we use homogeneous Dirichlet boundary conditions for the velocity: where the turbulent viscosity is locally defined on each mesh cell κ by and where the trilinear form d on X 3 is defined by We initialize the process by choosing ū0 h = 0. We also choose x ∈ Ω, where the source term f results from a local average of the original source term f ; moreover we suppose that Note that this problem is linear, since we discretize all non-linear terms semi-explicitly.
We begin by showing a bound for the solution ūn with values in X nh × M nh .This solution satisfies, for m = 1, . . ., N , Proof.For ūn−1 h ∈ X (n−1)h , it is clear that Problem (F V n,h ) has a unique solution (ū n h , pn h ) as a consequence of the coercivity of the corresponding bilinear form on X nh × X nh and the inf-sup condition (8).The proof of (9) follows by testing (F V n,h ) with ūn h and summing the resulting inequalities.□

A posteriori error analysis
We now intend to prove a posteriori error estimates between the exact solution (u, p) of Problem (F V ) and the numerical solution of Problem (F V n,h ).The main result in this respect is given by Corollary 4.16.
Several steps are needed to obtain this result.We begin by constructing the a posteriori indicators in Section 4.1 and then bound the error by the estimates in Section 4.2.
4.1.Construction of the error indicators.We first introduce the space and, for 1 ≤ n ≤ N , we fix an approximation f n h of the data f n in Z nh .This is a technicality which is needed in Section 5 in which we shall make use of Properties 4.1 and 4.2 below, which are valid only for polynomial functions.
Next, for every element κ in T nh , we denote by • ε κ the set of edges of κ that are not contained in ∂Ω, • ∆ κ the union of elements of T nh that share at least one common vertex with κ, • h e the diameter of the edge e ∈ ε κ , • [•] e the jump through e for each edge e in ε κ (specifying its sign is not necessary).
• n κ the unit outward normal vector to κ on ∂κ.
For the proofs of the next theorems, we introduce for an element κ of T nh , the bubble function ψ κ (resp.ψ e for the edge e) which is equal to the product of the 3 barycentric coordinates associated with the vertices of κ (resp. of the 2 barycentric coordinates associated with the vertices of e).We also consider a lifting operator L e defined on polynomials on e that vanish on ∂e into polynomials on the at most two elements (κ, κ ′ ) containing e and vanishing on ∂κ \ e and ∂κ ′ \ e.This lifting operator is constructed by affine transformation from a fixed operator on the reference element.We recall the next results from [30,Lemma 3.3].Property 4.1.Denoting by P r (κ) the space of polynomials of degree smaller than r on κ, we have Denoting by P r (e) the space of polynomials of degree smaller than r on e, we have , and, for all polynomials v in P r (e) vanishing on ∂e, if κ is an element which contains e, e ∥v∥ L 2 (e) .We also introduce a Clément type regularization operator C nh [13] which has the following properties, see [5, section IX.3]:For any function w in H 1 (Ω) 2 , C nh w belongs to the continuous affine finite element space and satisfies for any κ in T nh and e in ε κ , Note that we use the variant of C nh which ensures that C nh w belongs to H 1 0 (Ω) d when w belongs to H 1 0 (Ω) d (see [13]).Furthermore, we introduce the Scott-Zhang operator F nh [12] which has the following properties : For any function v ∈ H s (Ω) 2 , we have where s ∈] 1 2 , 1[ and t ∈ [0, s] (see [12]).For the a posteriori error studies, we consider the piecewise affine function ūh which takes in the interval and the piecewise constant function ph equal to pn h on the interval ]t n−1 , t n ].
Property 4.3.For any dimension d and for any non negative integer r, there exists a constant c such that for any polynomial function v h of degree r on κ We prove optimal a posteriori error estimates by using the norm: To prove the upper bound, we follow the idea used by Bernardi and Verfurth [9] or Bernardi and Sayah [7] for the Stokes problem in order to uncouple time and space errors.But in this work, the non linear term coming from the Navier-Stokes system requires more sophisticated calculations.
We introduce an auxiliary problem corresponding to the time discretization and calculate upper bounds for the errors between the corresponding solution and the exact solution firstly and the discrete filtered solution of (F V n,h ) secondly.Finally, we combine the obtained errors to derive the desired upper bound for the a posteriori error estimation.
Setting u 0 = 0 and setting dt for almost all x ∈ Ω, we introduce the following time semi-discrete problem: Lemma 4.4.Problem (P aux ) has a unique solution because the bilinear form is coercive (owing to the fact that div u n−1 = 0) and because of the inf-sup condition.Furthermore, we have: We define the piecewise affine function u τ by its value on the interval [t n−1 , t n ]: and we define p τ as the piecewise constant function equal to p n on the interval ]t n−1 , t n ].An easy calculation leads to the following lemma.
Lemma 4.5.By combining Problems (F V ) and (P aux ), we observe that the pair (17) where Lemma 4.6.By combining Problems (P aux ) and (F V n,h ), we observe for v ∈ X and v h ∈ X nh : and where the residuals R h,1 (ū h ) and R h,2 (ū h ) are defined by: where ν th is a piecewise constant approximation of ν t (a technicality needed in Section 5 in which we shall use Properties 4.1 and 4.2 valid for polynomial functions) defined on each κ and for every v ∈ X by: and where In order to derive the upper bounds, we use the integration by parts formula to rewrite the residual operators R τ (u h )(t), R h,1 (u h )(t) and R h,2 (u h )(t) in the following forms: All this leads to the following definition of the error indicators: Remark 4.8.Even if these indicators are a little complex, each term in them is easy to compute since it only depends on the discrete solution and involves (usually low degree) polynomials.
Lemma 4.9.The following estimates hold for 1 ≤ n ≤ N , (1) For all v in X and v h = C nh v: (2) For all v ∈ X: Proof.We proceed in two steps, one for each estimate.1) We derive the result from formula (25) with v h = C nh v, by using the continuous Cauchy-Schwarz inequality, the properties of C nh given by ( 10) and the discrete Cauchy-Schwarz inequality.
2) From (26), using the continuous and discrete Cauchy-Schwarz inequalities, we derive (31).□ 4.2.Upper bounds of the error.In this section, we establish in Corollary 4.16 the upper bound corresponding to the difference in the solutions of Problems (FV) and (F V n,h ) with a condition on τ n and h n .The main idea is to decompose the error through the introduction of the solution u τ of Problem (P aux ), which leads to the intermediary Theorems 4.10 and 4.12, in which we apply the continuous and discrete Gronwall lemmas.
Theorem 4.10.Let u ∈ L ∞ (0, T, L 3 (Ω) 2 ).The following a posteriori error estimate holds between the velocity u of Problem (F V ) and the velocity u τ associated with the solutions Proof.By taking v = u − u τ and q = p − p τ in (17), we obtain: Let us bound the right-hand side of equation (33).The first and last terms can be bounded by using ( 1) and ( 2), respectively, as well as Lemma 2.11, as follows: Furthermore, the residual term in the right hand side of Equation (33) can be decomposed as: where Using the continuous and then discrete Cauchy-Schwarz inequalities, the fact that t − t n τ n ≤ 1 for all t ∈ [t n−1 , t n ] and using Lemma 2.11, we bound T 1 in the following way: In order to bound T 2 , we insert u(t, x) and u τ : Using that u τ = t − t n−1 τ n (u n − u n−1 ) + u n−1 , the fact that t − t n−1 τ n ≤ 1 and Lemma 2.11, we obtain: Using (2) as well as the fact that (15) and the triangular inequality imply that On the other hand, using Lemma 2.11, it holds that Using the discrete Cauchy-Schwarz inequality, then (2), (3) and Lemma 4.4, we obtain: Since v = u − u τ and using Lemma 2.11 once again, we obtain Now, we bound T 2,3 by taking into consideration that u is supposed to be in L ∞ (0, T, L 3 (Ω)); then we use the discrete Cauchy-Schwarz inequality, Lemmas 2.11 and 2.1 and we obtain In order to bound |T 3 |, we use the fact that and we insert u: By using that v = u τ − u and that c(u, v, v) = 0 the first term in the right-hand side vanishes.The second term is treated exactly like T 2,3 and we obtain Finally, gathering (33), (34), ( 35), ( 36), (37), ( 38), ( 39), ( 40), ( 41) and (42), we obtain on each Then, we simplify (43), take into account the definitions of π τ u and π l,τ u, and integrate equation between 0 and t; we obtain We apply the Gronwall Lemma 2.9 with functions given by the following form : We obtain the following bound: Since f is an increasing function of time, we bound it by f (t m ); we use (3) to bound the integrals of ||u|| 2 X .Applying the resulting inequality at t = t m , we get: To obtain (32), we use the triangle inequality below in the last term in the right-hand side of (44), then use hypothesis (7) and definition ( 27)

□
To derive an a posteriori estimate between the solution u of problem (FV) and the solution ūh corresponding to the solutions ūn h of (F V n,h ), it suffices to get an a posteriori estimate between the solution u τ of Problem (P aux ) and ūh and to apply the triangle inequality using the previous theorem.In order to get an a posteriori error estimate between the solutions u τ and ūh , we introduce the operator Π (see [9]) defined from X into itself as follows: For each v in X, Πv denotes the velocity w of the unique weak solution (w, r) in X × M of the Stokes problem: The next lemma states some properties of the operator Π.We refer to [ The following estimates hold for all v in X, (3) ∀v h ∈ V nh and 1 ≤ n ≤ N : . We are now in a position to prove a posteriori estimates between the solution u τ of Problem (P aux ) and the solution ūh of (F V n,h ).Theorem 4.12.Suppose there exists a positive constant C s such that for all 1 ≤ n ≤ N we have h n ≤ C s τ n .The following a posteriori error estimate holds between the solutions u m and ūm h of Problems (P aux ) and (F V n,h ).
Proof.For abbreviation we set e n = u n − ūn h and ε Adding and subtracting Πe n in both terms in the right-hand side of (47) and since div(e n − Πe n ) = 0, we obtain: By taking v = e n − Πe n in (19), we have for every Next, we evaluate all terms in the right-hand side of (49) separately by repeatedly using both inequalities in Lemma 2.11 for various values of (β, γ, δ).Taking into account that Πe n = −Πū n h and using Lemma 4.11, the first and second terms can be bounded as: and To estimate the third and fourth terms of (49), we take v h = C nh v, use estimation (10) and Lemma 4.11 to derive : The fifth and sixth terms of (49) can be bounded as (cf Lemma 4.9) Finally, we bound the last two terms of (49).We have the relation: We set ).We first bound A 1 by using the L 4 − L 2 − L 4 inequality, the fact that v = e n − Πe n and Πe n = −Πū n h , then (2) and Lemma 4.11; we obtain We bound separately the two terms that result from the expansion of the product in the right-hand side of (55), that we denote by A 1,1 and A 1,2 .First, using Lemma 2.11 we have X Inserting e n−1 in the last term and using that τ n ||u n || 2 X is bounded by (15), we can choose δ such that Secondly, using Lemma 2.11, we have: Using the relation h n ≤ C s τ n and the fact that (15) allows to bound τ 1/2 n ||u n || X by a constant, we have, using Lemma 2.11 twice Finally, summing (56) and (57), we get (58) As far as A 2 is concerned, we first use that div e n−1 = − div ūn−1 h , then use (2), Lemma 2.11 and a triangular inequality.We get, for any β > 0 Using that ).Since β in (59) can be freely chosen, this yields the following bound for A 2 Let us now bound the term B. It can be written as follows: Using the fact that Πe n = −Πū n h , Lemma 4.11, then (2) and Lemma 2.11, we get, for any β > 0 A term very similar to the second one in the right-hand side of (61) was already bounded in (59), replacing u n by ūn−1 h ; since (9) yields similar bounds as (15), choosing β small enough leads to Moreover, using once more the fact that Πe n = −Πū n h , Lemma 2.11, then (2) and Lemma 4.11 we get Using the fact that h n ≤ C s τ n , the fact that τ 1/2 n ||ū n h || X is bounded thanks to (9), Lemma 2.11 and inserting e n−1 , we obtain Using once more that τ n ||ū n h || 2 X is bounded, we may choose β such that Thus, by summing (50), (51), ( 52), ( 53), ( 54), ( 58), ( 60), ( 62) and (63), using Equations ( 47), (49), the fact that e 0 = 0, and summing over n from 1 to m, we obtain: Moreover, we perform a change of indices, we use (7) to bound τ n+1 by σ τ τ n ; we also use the fact that e 0 = 0 and ū0 h = 0 and the fact that || div ūn We apply the Gronwall Lemma 2.10 with the following functions: By remarking that for every n ≤ m it holds that f n ≤ f m , and using the bounds provided by ( 9) and ( 15), the last line in (65) can be bounded by Cf m , and we finally obtain (46).This concludes the proof of Theorem 4.12 where we proved a posteriori error estimate between the velocity u τ of Problem (P aux ) and the velocity ūh of Problem (F V n,h ).□ The following lemma relates the integral of ||(u τ − ūh )|| 2 X with its values at the different time steps Lemma 4.13.It holds that Proof.For the proof of this lemma, we refer to [6] page 15. □ Corollary 4.14.We suppose The following a posteriori error estimate holds between the velocity u solution of problem (F V ) and the solution ūh associated to the solutions ūn h of problem (F V n,h ): Proof.We use the triangle inequality in the left-hand side of (67) by inserting u τ (t m ) = u m in its first term and the function u τ in its second term.Then, the proof is a direct consequence of Theorems 4.10 on the one hand and of the second inequality in (66) and Theorem 4.12 on the other hand.□ Next, the following theorem can be found in Theorem 4.10 in [7] as far as the first inequality is concerned and follows the same steps as Theorem 4.10 in [6] as far as the second inequality is concerned.
Theorem 4.15.The following a posteriori error estimate holds between the solution (u,p) of Problem (FV) and (ū h , π τ p τ ) associated with the solutions of Problem (F V n,h ): and Finally, Definition ( 14) of the [[•]] norm, inequalities (68), ( 67) and (69) lead to the following result: The pressure and the velocity verify the following a posteriori error bound: Remark 4.17.Due to the use of the Gronwall Lemma and of Lemma 2.11, the dependence of the constant C in (70) with respect to ν involves an exponential of 1 ν .This may look very unfavorable as ν is small in turbulent flows.However, from a practical point of view, numerical experiments presented in Section 6.2 show that this dependence is not that dramatic.

Upper bounds of the indicators
In this section, we prove three theorems regarding upper bounds for the error indicators proposed in Section 4. We follow exactly the same steps like in [6,Theorems 4.11 and 4.12] by simply changing the form of the non-linear terms.Taken together, these theorems will establish the following general estimate where L wκ is an expression depending on the errors, on the data approximations, and on terms related to the LES viscosity, evaluated on the union of the elements of T nh sharing a face with κ.The discussion is presented in terms of three theorems, each of which providing an upper bound of one term in the left-hand side of (71).In all, suitable regularity assumptions on the exact solution are performed.
Theorem 5.1.Let D(u) ∈ L ∞ (0, T, Ω) and suppose that h n ≤ C s τ n .The following estimate holds: where w κ denotes the union of the elements of T nh that share at least a face with κ and Proof.We proceed in 4 steps: )D(u(t)), then ν t (u(t))D(u(t)) and ν th (u(t))D(u(t)) in the fourth term of (η h,1 n,κ ) 2 defined in (28) and we obtain : We integrate from t n−1 to t n , we bound the first term in the right-hand side of (74) by: By using the inverse inequality ( 13), the fact that h κ ≤ cτ n and equation ( 9), we have: .
By using that ||D(u)|| L ∞ (0,T,Ω) is bounded, we bound the second term in the right-hand side of (74) by: The operator defined in (22) verifies the following property: For every v1 and v2 ∈ Z, we have: Indeed we have By using the Hölder inequality, we obtain (75).This leads to the following bound of the third term in the right-hand side of (74): .
By regrouping the previous inequalities, we obtain: which provides a bound for the fourth term of (η h,1 n,κ ) 2 .Next, the solution (u, p) of problem (F V ) and the solution (ū h , ph ) associated with the solution This equality will help us bound successively the first and second terms of (η h,1 n,κ ) 2 for a given κ.
(2) Regarding the first term in (η h,1 n,κ ) 2 , we set and choose in (77) : where ψ κ is the bubble function which is equal to the product of the barycentric coordinates associated with the vertices of κ.Integrating between t n−1 and t n and using the Cauchy-Schwarz inequality yields The Cauchy-Schwarz inequality allows to write: Multiplying the above inequality with h 2 κ and using Property 4.1, we obtain: Now, we replace v κ by α n h ψ k and use Property 4.1 to get rid of ψ 1/2 k in the norm in the left-hand side and to get rid of ψ k in the norm of v κ in the right-hand side.Moreover, we use Lemma 2.11 and we obtain We insert ν By integrating between t n−1 and t n , and using the inverse inequality (13), the fact h κ ≤ cτ n and the bound (9), we bound the first term in the right-hand side of (79) by: By using that D(u) ∈ L ∞ (0, T, Ω), the second term in the right-hand side of (79) can be bounded by: Integrating (79) between t n−1 and t n , using the above two inequalities and (76), we obtain: We replace (82) into (78) and simplify, we obtain: We replace α n h by its value and we obtain: where we recall that L(•) is defined in (73).
(3) Regarding the second term in (η h,1 n,κ ) 2 , for every e ∈ ε κ , we note by κ ′ the other element containing e.We introduce the function: and, recalling that L e is defined at the beginning of Section 4, we choose in (77) v = v e = L e (R h,1 n,e ψ e ) extended by 0 outside of κ := κ ∪ κ ′ .Integrating by parts the diffusion and pressure terms, we notice that the resulting trace terms on e sum up to e R h, 1  n,e (R h,1 n,e ψ e ); moving all other terms in the right-hand side of the equality, integrating between t n−1 and t n and applying the Cauchy-Schwarz inequality we obtain: Multiplying by h e , and using Property 4.2 for v e = L e (R h,1 n,e ψ e ), we obtain: Using (82), the fact ψ e ≤ 1 implies ψ e ≤ ψ 1 2 e and using Lemma 2.11, we obtain: Then we simplify, use the first property of Proposition 4.2, replace R h,1 n,e by its value, sum over e ⊂ ∂κ using that h e ≤ h κ and then use (83); we obtain: (4) Regarding the third term in (η h,1 n,κ ) 2 , since t → ūh (t, x) is an affine function with value ūn h (x) in t = t n and using div u(t, x) = 0, we have Combining ( 76), ( 83), ( 84) and (85), we obtain the desired result.□ Secondly, we estimate the indicator η h,2 n,κ .
Theorem 5.2.Let D(u) ∈ L ∞ (0, T, Ω) and suppose that h n ≤ C s τ n .We have the following estimate: Proof.We write that τ n (η h,2 n,κ ) 2 = tn tn−1 (η h,2 n,κ ) 2 dt and we insert ν t (u)(t)D(u)(t) and ν t (u)(t)D(ū n h ) into (η h,2 n,κ ) 2 and then use triangular inequalities.Using (80) and ( 81) we obtain the result.□ Remark 5.3.In Theorems 5.1 and 5.2, the term ||ν t (u)D(u)|| 2 L 2 (tn−1,tn,L 2 (κ)) comes from the fact that the LES solves a system that differs from the original Navier-Stokes equations.It is of order h 4 κ because ν t scales like h 2 κ and we remark that it is a higher order term as compared to other error terms.
To conclude, we estimate η τ n,κ .Theorem 5.4.We have the following estimate: Proof.Using the definition of ūh (t), we obtain ∀t ∈ [t n−1 , t n ] : Inserting u(t) in the right-hand side of this equality and taking the gradient of both sides, integrating on κ and using the triangular inequality, we obtain Integrating between t n−1 and t n we obtain the result with C = 6.□ As a final conclusion of Sections 4 and 5, we have obtained the equivalence of the error indicators of Definition 4.7 with the error between the exact solution of the Navier-Stokes equations and the discrete filtered solution, up to terms representing the data oscillation (with respect to time) for f , the change of source term due to the filtering process ( f n instead of f n ), the approximation of f n by f n h , that of ν t by ν th and finally a term related to the additional turbulent diffusion operator.

Numerical results
In this section, we present numerical simulations using the FreeFem++ code (see [20]) and the finite element spaces and time marching scheme discussed in Section 3.They show that LES provides lower error indicators and lower deviations from a reference solution than pure Navier-Stokes simulations; they also show that the efficiency index does not vary much with increasing Reynolds number and, finally, that, compared to uniform refinement, adaptive LES offers a substantial reduction in the number of degrees of freedom needed in the simulations to reach a given error indicator value.
We consider the domain Ω given by Figure 1  Concerning uniform meshes, we divide the edges of ∂Ω into segments of equal lengths, and we define M to be the number of such mesh segments per unit length.For M = 8 this leads to a mesh with around 900 vertices and 1800 triangles, for a total of around 6500 degrees of freedom We consider ν = 1 Re where Re is the Reynolds number and we choose the density of body forces f equal to (−2, 0) in the rectangle LGHK and to (0, 0) elsewhere.This implies, as shown on Fig. 3, that the flow will move from the right to the left at the center of the domain; when the flow hits the left boundary, it splits into an upper flow and a lower flow; the BCDE obstacle will cause recirculations in the upper flow.This will in turn generate turbulent interactions with the main flow at the center of the domain.In all the calculations, we chose the Smagorinsky constant to be c s = 0.1.
From the calculation of (ū n h , pn h ) solution of (F V n,h ) at each time step, we compute the error estimators of Definition 4.7; in the actual calculation we replace ν th by ν t for the sake of simplicity; this has almost no effect on the results.
To study the effect of the addition of the turbulent diffusion term in the LES method, we compare error indicators in time and space for simulations with or without LES for various values of Re ∈ {200, 5000} and M ∈ {8, 16, 32, 64}.In all simulations the final time is set to T = 3.We choose the time step to be ∆t = 1 4M .We define : .
The results are in Table 1 for Re = 200 and in Table 2 for Re = 5000.The total number of degrees of freedom in each simulation is indicated under the column # dof.It is clear that the space and time estimators using the LES method are smaller than those of the Navier-Stokes problem without LES.In addition, improvement is more important for high Reynolds numbers and coarse meshes, which is an indication that the LES method exactly plays the role expected from it.However, these conclusions are based only on observation of indicators; to confirm this trend, we will directly measure the difference between numerical solutions on different meshes and reference solutions obtained on the finest mesh that we could use given the resources at our disposal.This is the subject of the next sub-section.

6.2.
Comparisons of errors and estimators.In this section, we compare the relative norms of differences between a reference solution u F obtained by a simulation without LES on a fine mesh (M = 256 and ∆t = 0.25/M ) and solutions u C of simulations on coarse meshes without LES on the one hand and with LES on the other hand for different Reynolds numbers.This serves us as a measure of the errors between the numerical solutions obtained on coarse meshes and the exact, unknown solution.We also compare the error estimators and these measures of the errors in order to estimate the unknown constants appearing in the upper bounds of Section 4, see Remark 4.17.The relative difference norm is defined by: .
The relative norms of the differences are in Table 3 and Table 4 for Re ∈ {200, 5000} and T = 3.We note M F = 256 the number of segments per unit length used to cut the edges of ∂Ω in the fine mesh and M C ∈ {8, 16, 32, 64} in the different coarse meshes.In addition, the difference between the solutions with and without LES is also calculated on the fine mesh (this is noted M F = 256 and M C = 256 in the tables below).
The tables show a general trend indicating that the differences between the LES solutions and the reference solution are smaller than those observed without LES.This trend is particularly notable for high Reynolds numbers and coarse meshes, on which solutions without LES are particularly inaccurate.This broadly confirms the observations made from the error estimators in the previous section.Table 3. Relative norm of the difference between the solution of the NS problem on a reference fine mesh and solution of the NS problem with and without LES on coarse meshes for Re = 200 and T = 3.
Figure 2 shows the ratio of e F C in Tables 3 and 4 with estimator η in Tables 1 and 2 (the same vertical scale was chosen on both figures).This figure shows that this ratio, and thus the constant C appearing in (70), seem to increase with the Reynolds number, but not dramatically.On coarse meshes and without LES this constant seems to be larger, but remains moderate with LES.

6.3.
Comparison between uniform and adapted LES problems.In the previous sections, we proved that solutions obtained with LES are more accurate than those obtained without LES.Therefore, we only use LES in this sub-section.We are now seeking to obtain more accurate simulations by using an algorithm for adapting the mesh and the time step using the associated estimators, seeking to balance the time estimator and the space estimator.We propose an adaptation algorithm described below.We measure its efficiency by comparing the estimators calculated on the one hand on uniform meshes and with  4. Relative norm of the difference between the solution of the NS problem on a reference fine mesh and solution of the NS problem with and without LES on coarse meshes for Re = 5000 and T = 3.   constant time steps and on the other hand on meshes and with time steps adapted using our algorithm, which is driven by a tolerance denoted ε set by the user, chosen all the smaller as the desired precision is high.At each time step, knowing ūn h , the following operations are performed Figure 4 shows the resulting meshes at different values of t, for the particular value ε = 0.08.Since we start from a vanishing velocity at t = 0 and since the momentum source term is piecewise constant in the computational domain, large variations of the velocity and its gradient occur near the discontinuity line of the source term; this is well captured by the mesh which concentrates cells around this line at t = 1.Then, as complex interactions with boundaries occur, the mesh at t = 3 is refined in the vicinity of some parts of the boundary, including, as expected, that of the upper obstacle, since geometric singularities usually affect negatively the regularity of the solution.It is important to note that mesh density is not necessarily driven by velocity modulus, since the estimators rather depend on variations of the velocity gradients.Additionally, many cells are avoided in a wide part of the domain, in particular in the vicinity of a large part of the boundary.This results in huge savings in terms of computational resources: the smallest cell length is around 2.1 × 10 −3 , while the largest is around 5.7 × 10 −1 .For Re = 1000, Figure 5 shows a comparison, in logarithmic scale, of the error indicators for the uniform LES problems and adapted LES problem with respect to the number of degrees of freedom.The uniform curve was obtained by simulations on several uniform meshes with constant time steps.For each of these meshes the number of degrees of freedom is the product of the number of time steps by the number of unknowns discretizing the problem on the considered mesh.The adaptive curve was obtained by choosing several different values of ε (see algorithm above) and for each of these values, we added the number of degrees of freedom at each time step, which may vary during the calculation when the algorithm decides that a mesh adaptation is necessary.This figure shows the effectiveness of the adapted algorithm versus the uniform one.

Conclusions and perspectives
We have obtained reliable and efficient a posteriori error estimators between the numerical solution of LES equations and the exact solution of the Navier-Stokes problem in two space dimensions.The error indicators are of three types: one related to the time discretization, one to the space discretization and finally one to the additional diffusion term introduced by the Smagorinsky model for large eddies.Numerical simulations have shown the usefulness of such indicators to perform adaptive LES.
A first possible extension of this work is to investigate the case of three space dimensions; even under idealized existence and uniqueness assumptions, we have noticed in Remark 2.3 that one of our main tools in two dimensions does not carry over to three dimensions, so that other techniques are to be found.Furthermore, in practice, using homogeneous Dirichlet boundary conditions for the discrete filtered velocity is not the method of choice because this implies resolving the viscous boundary layers, which may be particularly thin for high Reynolds numbers; practitioners often resort to "laws of the wall", which would then need to be included in the analysis in order to account for the additional error that these boundary conditions generate.Another possible additional research path is to consider other types of finite elements, be they conforming like the Taylor-Hood pairs [28] or non-conforming like the Crouzeix-Raviart [14] or the Matthies-Tobiska [23] families.Finally, other, more sophisticated, LES models are routinely used in Computational Fluid Dynamics, like for example that of Germano et al. [17]; it would be worth extending the present work to these alternative LES models.

6. 1 .
Comparison of error indicators for calculations with or without LES.The indicator is defined by:

e
F C without LES e F C with LES M F = 256 and M C = 8 0.743 0.628 M F = 256 and M C = 16 0.359 0.349 M F = 256 and M C = 32 0.164 0.166 M F = 256 and M C = 64 0.078 0.079 M F = 256 and M C = 256 − 5.19 × 10 −4

e
F C without LES e F C with LES M F = 256 and M C = 8 3.052 0.983 M F = 256 and M C = 16 2.617 0.963 M F = 256 and M C = 32 1.612 0.840 M F = 256 and M C = 64 0.975 0.694 M F = 256 and M C = 256 − 0.246 Table of edges per unit length Ratio of error over estimator LES, Re=5000 LES, Re=200 of edges per unit length Ratio of error over estimator NS, Re=5000 NS, Re=200

Figure 2 .
Figure 2. Estimation of ratio of errors over estimators with (left) and without (right) LES.
we calculate η h n and η τ n defined in (88), (3) If η h n + η τ n > ε, we adapt in time or in space : (a) If η h n < η τ n and if the time step τ n is greater than a minimum time step (∆t) min , we adapt in time so that the next estimator in time is slightly less than the estimator in space (b) If η h n < η τ n and τ n = (∆t) min , then (i) if η h n > ε 2 , we adapt in space.(ii) Otherwise, adaptation in time is impossible and adaptation in space useless.(c) If η h n > η τ n (i) If the maximum number of mesh refinements is not reached, we adapt in space (ii) Otherwise we keep the calculation (4) If η h n + η τ n < 0.9ε we can increase the time step; in doing so we try to balance the space and time estimators (5) Otherwise we ensure that the estimators are balanced (a) If η τ n ≪ η h n , the time step is increased to balance the estimators, even if this means having to refine in space at the following iteration.(b) If the estimators are balanced, nothing is done.

Figure 3 .
Figure 3. Velocity modulus for Re = 1000 and t = 1 (upper row) and t = 3 (lower row, with a different color scale).

Figure 5 .
Figure 5.Comparison between uniform and adaptive LES estimators with respect to the total numbers of degrees of freedom in logarithmic scale.

Table 1 .
Comparison between estimators with and without LES for Re = 200 and T = 3.

Table 2 .
Comparison between estimators with and without LES for Re = 5000 and T = 3.