Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier-Stokes equations

S  eMA Journal n o 60(2012), 51–74 ERROR ANALYSIS OF A SUBGRID EDDY VISCOSITY MULTI-SCALE DISCRETIZATION OF THE NAVIER–STOKES EQUATIONS CHRISTINE BERNARDI∗, TOM ´ AS CHAC ´ ON REBOLLO†, MACARENA G ´ OMEZ M ´ ARMOL‡ ∗Laboratoire Jacques-Louis Lions, C.N.R.S & Universit´ e Pierre et Marie Curie – Paris 6, 4 Place Jussieu, 75252 Paris Cedex 05, France †Departamento de Ecuaciones Diferenciales y An´ alisis Num´ erico & IMUS, Apdo. de correos 1160, Universidad de Sevilla, 41080 Sevilla, Spain ‡Departamento de Ecuaciones Diferenciales y An´ alisis Num´ erico, Apdo. de correos 1160, Universidad de Sevilla, 41080 Sevilla, Spain [email protected], [email protected], [email protected] Abstract We propose a finite element discretization of the Navier–Stokes equations that relies on the variational multi-scale approach together with the addition of a Smagorinsky type viscosity, in order to take into account possible subgrid turbulence. We recall that the discrete problem admits a solution and prove a priori error estimates. Next we perform the a posteriori analysis of the discretization. Some numerical experiments justify the interest of this approach. 1 Introduction This paper deals with the numerical approximation of incompressible flows in turbulent regime by means of grid adaptation techniques. Grid adaptation techniques are currently used to solve fluid flow problems, providing large savings of computational complexity. The numerical analysis of this technique is largely based upon a well-sound mathematical analysis of the problem considered. However, the mathematical analysis of many standard turbulence models is not well developed up-to-date. Even more, some commonly used models, in particular the k-εone, do not seem to be well posed from the mathematical point of view. An emerging class of turbulence models that is increasingly used due to its accuracy and simplicity is provided by the Variational Multi-Scale (VMS) setting (see [16] for a general description of VMS models). This is a fully discrete model, that does not require a continuous modeling step. The VMS procedure yields a discrete equation for a finite element approximation of the flow, where the eddy viscosity only acts on the Received: April 10, 2012. Accepted: August 14, 2012. This research was partially supported by Spanish Government MTM 2009-07719 Grants. 51 52 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol small resolved scales of the flow. A simple modeling of the eddy viscosity acting on these scales (the Smagorinsky model) is used. This kind of combination of variational multi-scale and Smagorinsky models for small scales yields numerical results similar to those provided by up-to-date Large Eddy Simulation (LES) models (see [17], [18], [19]). We consider in this paper a VMS-Smagorinsky turbulence model for which we perform a numerical analysis that extends the usual one for standard discretizations of Navier–Stokes equations. Let Ωbe a bounded connected domain in Rd,d= 2 or 3, with a Lipschitzcontinuous boundary ∂Ω. We are interested in the finite element discretization of the Navier–Stokes equations in this domain relying on the variational multi-scale method, in cases where turbulence phenomena may occur. More precisely, (i) For both the velocity and the pressure, we introduce two spaces of discretization, one called coarse and the other one called fine: Indeed, the fine space either is defined from a mesh which is refined from the mesh used for the coarse space or involves higher degree polynomials. (ii) A further nonlinear viscosity term is added. In standard turbulence models, this viscosity can depend on other unknowns which are the solution of convectiondiffusion equations, e.g. the turbulent kinetic energy (see [2] for instance), or the temperature (see [3] for instance). Here, we have chosen to work with the well-known nonlinear Smagorinsky viscosity, introduced in [22]. We refer to [8] (Chap. 2) for a complete description of this discretization which brings to light its interest for the approximation of turbulent flows. We thus consider a variational multi-scale approximation of Navier–Stokes equations in turbulent regime for which the Smagorinsky eddy viscosity acts only on the small resolved scales of the flow. In a rather general finite element framework, relying on standard arguments for nonlinear problems, we perform the a priori and a posteriori analysis of the discrete problem. This leads to optimal error estimates. A few numerical experiments are in good coherence with the theoretical results. An outline of the paper is as follows: •In Section 2, we present the continuous and discrete problems we work with and recall their main properties. •Section 3 and 4 are devoted to the a priori and a posteriori analysis of the discretization, respectively. •Numerical experiments are presented in Section 5. •The final section states the main conclusions of the paper. Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 53 2 The continuous and discrete problems We have decided for simplicity to work with homogeneous no-slip boundary conditions. In this case, the Navier–Stokes equations in Ωread    −ν∆u+ (u·∇)u+∇p=fin Ω, ∇·u= 0 in Ω, u=0on ∂Ω. (1) The unknowns are the velocity uand the pressure pof the fluid. The data are only the distribution fwhich represents a density of body forces, while the viscosity νof the fluid is a positive constant. 2.1 Variational formulation of the continuous problem We consider the full scale of Sobolev spaces Hs(Ω),s∈R, and Wm,p(Ω),m∈N, 1≤p≤ ∞, equipped with the standard norms and seminorms. In order to write a variational formulation of problem (1), we also introduce the space H1 0(Ω) of functions in H1(Ω) vanishing on ∂Ωand its dual space H−1(Ω). We finally need the space L2 0(Ω) = {q∈L2(Ω); ∫Ω q(x)dx= 0}. Standard density results yield that system (1) (where the first two lines are satisfied in the distribution sense) is fully equivalent to the following variational problem: Find (u, p) in H1 0(Ω)d×L2 0(Ω) such that {a(u,v) + c(u,u,v) + b(v, p) = ⟨f,v⟩,∀v∈H1 0(Ω)d, b(u, q) = 0,∀q∈L2 0(Ω),(2) where ⟨·,·⟩stands for the duality pairing between H−1(Ω)dand H1 0(Ω)d. The bilinear forms a(·,·)and b(·,·)are defined by a(u,v) = ν∫Ω (∇u)(x) : (∇v)(x)dx, b(v, q) = −∫Ω (∇·v)(x)q(x)dx, while the trilinear form c(·,·,·)is given by c(w,u,v) = ∫Ω ((w· ∇)u) (x)·v(x)dx. We now recall the main properties of this problem from [13](Chap. IV, Sect. 2.1) and [23](Chap. 2, Sect. 1), and also the regularity results on its solution from [14](Section 7.3.3) and [9]. The form a(·,·)is continuous on H1(Ω)d×H1(Ω)dand elliptic on H1 0(Ω)d. The form b(·,·)is continuous on H1(Ω)d×L2(Ω) and satisfies the following inf-sup condition, for a constant β0>0, ∀q∈L2 0(Ω),sup v∈H1 0(Ω)d b(v, q) ∥v∥H1(Ω)d≥β0∥q∥L2(Ω).(3) 54 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol The form c(·,·,·)is continuous on H1(Ω)d×H1(Ω)d×H1(Ω)dowing to Sobolev imbeddings and satisfies the following anti-symmetry property, valid for any divergence-free function win H1(Ω)d, ∀v∈H1 0(Ω)d, c(w,v,v) = 0.(4) By combining all these properties with Brouwer’s fixed-point theorem, the following results can be derived. Proposition 1 (i) For any datum fin H−1(Ω)d, problem (2) admits at least a solution. Moreover this solution satisfies ∥u∥H1(Ω)d+∥p∥L2(Ω) ≤C∥f∥H−1(Ω)d where the constant Conly depends on Ω,νand β0. (ii) If the datum fand the viscosity νsatisfy N ν2∥f∥H−1(Ω)d<1, where Nis the norm of c(·,·,·), this solution is unique. (iii) If the datum fbelongs to L2(Ω)d, the solution (u, p)belongs to Hs+1(Ω)d× Hs(Ω) for all s≤s0, with 1 2≤s0≤1in the general case, 1 2< s0≤1when Ω is a polygon or a polyhedron, and s0= 1 when Ωis convex. 2.2 The discrete problem From now on, we assume that Ωis a polygon (d= 2) or a polyhedron (d= 3). We introduce a regular family (Th)hof triangulations of Ωby triangles or tetrahedra, in the sense that, for each h: •Ωis the union of all elements of Th; •The intersection of two different elements of Th, if not empty, is a vertex or a whole edge or a whole face of both of them; •The ratio of the diameter hKof any element Kof Thto the diameter of its inscribed circle or sphere is smaller than a constant σindependent of h. As usual, hstands for the maximum of the diameters hK,K∈ Th. For each h, we introduce a pair of finite element spaces Xhand Mhassociated with the triangulation Thon which we make the following assumption, in order to work with a conforming discretization Xh⊂H1 0(Ω)dand Mh⊂L2 0(Ω). In view of a variational multi-scale discretization, we also consider another pair of finite element spaces Xh′and Mh′which will be finer than the previous ones, in a sense which is made precise later on. There also, we assume that Xh′⊂H1 0(Ω)dand Mh′⊂L2 0(Ω). Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 55 We make the further assumption that both intersections Xh∩Xh′and Mh∩Mh′are reduced to {0}and set Xh=Xh⊕Xh′and Mh=Mh⊕Mh′.(5) Remark 1 Most often, the spaces Xhand Mhare constructed from Xhand Mhin one of the following ways: (i) They are built with polynomials of the same degree as Xhand Mhbut associated with a triangulation Th′constructed from Thby a refinement; (ii) They are still associated with the triangulation Thbut higher degrees of polynomials are used on each Kin Th, In both cases, there is not a unique way to build the spaces Xh′and Mh′to have (5). Space Xh′, for instance, may be constructed by means of a surjective linear operator of restriction or interpolation Πh:Xh→Xh, by Xh′= (Id −Πh)Xh, where Id is the identity operator. Space Mh′is built from Mhsimilarly. In this way Xh′and Mh′do not need to be explicitly constructed. We model the eddy viscosity by means of the Smagorinsky model: We associate with each function vin H1 0(Ω)dits eddy viscosity νS(v)defined by ∀K∈ Th, νS(v)|K= (CShK)2|∇v|, where |·|here denotes the Euclidean norm on Rd×d. The quantity CSis called Smagorinsky constant, and several values of it have been proposed. It is typically equal to 0.18 (see Germano [11], [12]) although it can be dynamically adapted in a time-dependent computation (see Lilly [20]). It can also be adjusted close to solid walls, in order to avoid over-diffusion (see Van Driest [24]). In the sequel we shall assume that the interpolation operator Πhis defined from H1 0(Ω)donto Xhand satisfies the following stability property: ∀v∈H1 0(Ω)d,∥Πhv∥H1(Ω)d≤C∥v∥H1(Ω)d.(6) We are thus in a position to write the discrete problem. It reads: Find (uh, ph) in Xh×Mh such that {a(uh,vh) + c(uh,uh,vh) + b(vh, ph) + aS(uh,vh) = ⟨f,vh⟩,∀vh∈Xh, b(uh, qh) = 0,∀qh∈Mh, (7) where the “Smagorinsky” eddy viscosity form aS(·,·)is now defined by aS(u,v) = ∫Ω νS(uh′)(x) (∇uh′)(x) : (∇vh′)(x)dx 56 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol with uh′= (Id −Πh)u,vh′= (Id −Πh)v. It can be noted that, up to the eddy viscosity term aS(uh,vh), this problem is constructed from (2) by the Galerkin method. This term models the sub-grid eddy viscosity effects, that are taken into account by means of the Smagorinsky term with a projection term (Id −Πh) that filters out the action of eddy viscosity on the large resolved scales. Remark 2 The standard Smagorinsky model models the eddy viscosity effects by the form faSdefined by faS(u,v) = ∫Ω νS(u)(x) (∇u)(x) : (∇v)(x)dx. This form faSdoes not include the projection term, so the eddy viscosity acts on both large and small resolved scales. This produces an over-diffusive effect, that the projection term intends to correct. Remark 3 Let us define the residual of the Navier–Stokes equations by duality: For any triplet (v;v, q)in H1(Ω)d×H1(Ω)d×L2(Ω), the quantity R(v;v, q)belongs to H−1(Ω)dand satisfies ⟨R(v;v, q),w⟩=a(v,w) + c(v,v,w) + b(w, q)−⟨f,w⟩,∀w∈H1 0(Ω)d. Then, problem (7) is equivalent to the following variational multi-scale method Find (uh, ph) satisfying uh=uh+uh′, ph=ph+ph′,(8) where the pair (uh, ph)is a solution of the problem Find (uh, ph) in Xh×Mh such that {a(uh,vh) + c(uh,uh,vh) + b(vh, ph) = −⟨R(uh;uh′, ph′),vh⟩,∀vh∈Xh, b(uh, qh) = −b(uh′, qh),∀qh∈Mh, (9) and the pair (uh′, ph′)is a solution of the problem Find (uh′, ph′) in Xh′×Mh′ such that    a(uh′,vh′) + c(uh,uh′,vh′) + b(vh′, ph′) + aS(uh′,vh′) = −⟨R(uh;uh, ph),vh′⟩,∀vh′∈Xh′, b(uh′, qh′) = −b(uh, qh′),∀qh′∈Mh′. (10) Indeed, observe that, when taking vhequal to vh(in Xh) in (7), the term aS(uh,vh) vanishes. It follows that if (uh, ph)satisfies (7), then (uh, ph)satisfies (9). Also, (10) follows from (7) by taking vh=vh′and ph=ph′as test functions. Finally (7) follows from (9)-(10) by summing up these equations, and using that from (8), uh′= (Id −Πh)uh′. Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 57 This new formulation brings to light the fact that the Smagorinsky eddy viscosity term only acts on the small scales of the discretization. Moreover, the method (7) is a method with two grids, but it only needs a grid and an interpolation operator on a virtual coarser grid to be programmed. In fact, method (7) includes three grid levels: Large resolved scales (those of Xh), small resolved scales (those of Xh′) and un-resolved scales (the remaining scales, that are taken into account by means of the eddy diffusion term). There also, the existence of a solution to problem (7) can be derived by using Brouwer’s fixed point theorem. However we prefer to postpone this proof to the next section where a more precise result is established. 3 A priori analysis As now standard for nonlinear problems, the a priori analysis of the discrete problem (7) is performed thanks to the discrete implicit function theorem due to Brezzi, Rappaz and Raviart [7]. This requires some further notation that we now introduce. We are led to make two further hypotheses: (i) There exists a constant βindependent of hsuch that ∀qh∈Mh,sup vh∈Xh b(vh, qh) ∥vh∥H1(Ω)d≥β∥qh∥L2(Ω).(11) (ii) For each nonnegative integer k, let Pk(K)denote the space of restrictions to Kof polynomials with dvariables and total degree ≤k. Then, Xhcontains the space X1 h={vh∈H1 0(Ω)d;∀K∈ Th,vh|K∈P1(K)d},(12) and Mhcontains either the space M1 h={qh∈H1(Ω) ∩L2 0(Ω); ∀K∈ Th, qh|K∈P1(K)},(13) or the space M0 h={qh∈L2 0(Ω); ∀K∈ Th, qh|K∈P0(K)}. This last assumption is satisfied by all the finite element spaces we work with. In what follows, C,C′,... stand for generic constants that may vary from line to line, but are always independent of h. 3.1 Some notation Let Sdenote the Stokes operator, more precisely the operator which associates with any datum fin H−1(Ω)dthe part uof the solution (u, p)of the problem Find (u, p) in H1 0(Ω)d×L2 0(Ω) such that {a(u,v) + b(v, p) = ⟨f,v⟩,∀v∈H1 0(Ω)d, b(u, q)=0,∀q∈L2 0(Ω).(14) 58 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol It follows from the properties of the forms a(·,·)and b(·,·)stated in Section 2.1 that this problem is well-posed. So, the operator Sis well defined and continuous from H−1(Ω)dinto H1 0(Ω)d. Owing to the inf-sup condition (3), it is readily checked that problem (2) can equivalently be written as a fixed-point equation F(u) = u−SG(u) = 0,with G(u) = f−(u· ∇)u.(15) Similarly, let Shbe the discrete Stokes operator, i.e., the operator which associates with any datum fin H−1(Ω)dthe part uhof the solution (uh, ph)of the problem Find (uh, ph) in Xh×Mh such that {a(uh,vh) + b(vh, ph) = ⟨f,vh⟩,∀vh∈Xh, b(uh, qh) = 0,∀qh∈Mh.(16) There also, owing to (11) this operator is well-defined. Let us now introduce the mapping associated with the Smagorinsky form: It is defined from H1 0(Ω)dinto H−1(Ω)dby ⟨AS(u),v⟩=aS(u,v),∀u∈H1 0(Ω)d,∀v∈H1 0(Ω)d Thus, problem (7) can equivalently be written Fh(uh) = uh−ShGh(uh) = 0,with Gh(u) = G(u)−AS(u).(17) To go further, we need some properties of the operator Sh. 3.2 Basic properties of the discrete Stokes operator The stability property of the operator Shis immediately derived by taking vhequal to uhin problem (16). Lemma 2 The operator Shsatisfies the following stability property: For all fin H−1(Ω)d, ∥Shf∥H1(Ω)d≤Csup vh∈Xh ⟨f,vh⟩ ∥vh∥H1(Ω)d .(18) We refer to [6](§IV.2) and [13](Chap. II) among others for the following convergence properties which require the assumptions on Xhand Mhmade in the beginning of the section. Lemma 3 The operator Shsatisfies the following convergence property: For all fin L2(Ω)d, ∥(S −Sh)f∥H1(Ω)d≤C hs0∥f∥L2(Ω)d, where the real number s0is introduced in Proposition 1. From Lemmas 2 and 3, we easily derive that, for all fin H−1(Ω)d, lim h→0∥(S −Sh)f∥H1(Ω)d= 0,(19) which will be of great use in what follows. Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 59 3.3 Preliminary lemmas From now on, we assume that fbelongs to L2(Ω)dand that (u, p)is a nonsingular solution of problem (2), in the sense made precise in [7], see also [13](Chap. IV, Sect. 3.1): The operator DF(u) = Id − SDG(u)is an isomorphism of H1 0(Ω)d (where Ddenotes the differential operator with respect to u). It can be noted that this assumption is much less restrictive than the global uniqueness of the solution (see part (ii) of Proposition 1) since it only ensures its local uniqueness. We denote by L(E, F)the space of linear continuous mappings from a normed space Eon another normed space Fand by Ethe space of endomorphisms of H1 0(Ω)d: E=L(H1 0(Ω)d, H1 0(Ω)d). Lemma 4 If property (6) holds, there exists an h0>0such that, for all h≤h0, the operator DFh(u)is an isomorphism of H1 0(Ω)dand the norm µof its inverse is bounded independently of h. Proof: We have DFh(u) = DF(u) + (S −Sh)DG(u) + Sh(DG(u)−DGh(u)). Since DF(u)is an isomorphism of H1 0(Ω)d, the desired property will be established if the last two terms in the previous expansion tend to zero. We study successively these two terms. 1) We have DG(u)·w=−(u· ∇)w−(w· ∇)u.(20) Since s0>1 2, see Proposition 1, the mapping: w7→ DG(u)·wsends the unit sphere of H1 0(Ω)dinto a compact subset of H−1(Ω)d. Thus, since for all ε > 0, a compact subset admits an overlap by a finite number of balls with centre fiand radius ε, applying (19) to all these fiyields that lim h→0∥(S −Sh)DG(u)∥E= 0. 2) Since DG(u)−DGh(u)is equal to DAS(u)and thanks to (18), we have to evaluate the quantity, for wrunning through the unit ball of H1 0(Ω)dand zh running through the unit ball of Xh, BS(u;w,zh) = C2 S∑ K∈Th h2 K(∫K|∇u∗|(x)(∇w∗)(x):(∇z∗ h)(x)dx +∫K(∇u∗:∇w∗ |∇u∗|)(x) (∇u∗)(x) : (∇z∗ h)(x)dx), where for brevity, for any function v, we denote by v∗the function (Id −Πh)v. By using Cauchy–Schwarz inequalities and also property (6), we obtain |BS(u;w,zh)| ≤ C2 S∑ K∈Th h2 K∥∇u∥L2(K)d×d∥∇w∥L2(K)d×d∥∇zh∥L∞(K)d×d. 66 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol Proof: We bound successively the three terms in ηK. 1) If ψKdenotes the bubble function on K(equal to the product of the Lagrange coordinates associated with the vertices of K), we set: vK={(fh+ν∆uh−(uh· ∇)uh−∇ph)ψKon K, 0elsewhere. Taking vequal to vK(which belongs to H1 0(Ω)d) in (32) gives ∥(fh+ν∆uh−(uh· ∇)uh−∇ph)ψ 1 2 K∥2 L2(K) ≤(∥u−uh∥H1(K)d+∥p−ph∥L2(K))∥vK∥H1(K)d+∥f−fh∥L2(K)d∥v∥L2(K)d. (35) Using standard inverse inequalities, see [25] (Lemma 3.3) for instance, and multiplying by hKthus yield hK∥fh+ν∆uh−(uh· ∇)uh−∇ph∥L2(K)d ≤C(∥u−uh∥H1(K)d+∥p−ph∥L2(K)+hK∥f−fh∥L2(K)d).(36) 2) To bound the terms on the edges or faces and, for each ein EK, we introduce a fixed lifting operator Le,K •that maps functions on evanishing on ∂e into functions on Kvanishing on ∂K \e, •and is constructed from a fixed lifting operator on the reference triangle or tetrahedron. Next, for each edge ein EKshared by the two elements Kand K′of Th, we set: ve={Le,κ([ν ∂nuh−phn]eψe)on κ∈ {K, K′} 0elsewhere, where ψeis now the bubble function on e. By taking vequal to vein (32) and using the same inverse inequalities as previously and other ones, see once more [25](Lemma 3.3), together with estimate (36), we derive h 1 2 e∥[ν ∂nuh−phn]e∥L2(e)d ≤C∑ κ∈{K,K′}(∥u−uh∥H1(κ)d+∥p−ph∥L2(κ)+hκ∥f−fh∥L2(κ)d). (37) 3) Finally, we set: qK= (∇·uh)χK, where χKis now the characteristic function of K. Taking qequal to qKin (33) immediately yields ∥∇·uh∥L2(K)≤ ∥u−uh∥H1(K)d.(38) The proposition is now a direct consequence of (36), (37) and (38). Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 67 Estimate (34) is fully local. But a direct consequence of it is that (∑ K∈T h η2 K)1 2≤C(∥u−uh∥H1(Ω)d+∥p−ph∥L2(Ω)+(∑ K∈T h h2 K∥f−fh∥2 L2(K)d)1 2). When compared with (31), this last estimate proves the optimality of the family of indicators (ηK)K∈T h. Moreover, since estimate (34) is local, it can be hoped that these indicators are an efficient tool for adapting the mesh. Remark 5 If the space Xhcontains the space X1 hdefined in (12), the following upper bound can be proved for the indicators ηSK ,K∈ T h, ∑ K∈T h η 3 2 SK ≤C(h hmin )1−d 4∑ K∈T h (ηK+hK∥f−fh∥L2(K))3 2,(39) where hmin stands for the minimum of the hK,K∈ Th. Even if the family of triangulations is uniformy regular (which means that hmin ≥c h), estimate (39) is not fully optimal with respect to (29) since it involves a bad power of the indicators. Proving local estimates seems rather difficult for general finite elements. 4.4 A look at the standard subgrid eddy viscosity discretization For a while, we consider the simple mono-scale discrete problem Find (uh, ph) in Xh×Mh such that {a(uh,vh) + c(uh,uh,vh) + b(vh, ph) + faS(uh,vh) = ⟨f,vh⟩,∀vh∈Xh, b(uh, qh) = 0,∀qh∈Mh, where the new Smagorinsky eddy viscosity form faS(·,·)is now defined by faS(u,v) = ∫Ω νS(u)(x) (∇u)(x):(∇v)(x)dx. By the same arguments as in Section 3, it is readily checked that this problem has a unique solution (uh, ph)in a neighbourhood of a nonsingular solution of problem (2) which still satisfies the a priori error estimates (22) and (23). On the other hand, for each Kin Th, we introduce the modified error indicator linked to the new Smagorinsky term gηSK = (CShK)2∥∇uh∥2 L4(K)d×d. The same arguments as for Proposition 10 lead to the estimate ∥u−uh∥H1(Ω)d+∥p−ph∥L2(Ω) ≤C(∑ K∈T h (η2 K+gηSK 2+h2 K∥f−fh∥2 L2(K)d))1 2, and the upper bound (34) still holds for the indicators ηK. However, concerning the indicators gηSK , we can only prove the analogue of (39), which is not optimal. This lack of optimality brings to light the interest of using a multi-scale discretization. 68 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol 5 Numerical experiments The computations are performed with the code FreeFem++ due to Hecht and Pironneau [15]. We have decided here to work with the final spaces associated with Taylor– Hood finite elements, see [6](Chap. VI, Example 3.5) or [13](Chap. II, Sect. 4.2) for instance. More precisely, these spaces are defined by Xh={vh∈H1 0(Ω)d;∀K∈ Th,vh|K∈P2(K)d}, Mh={qh∈H1(Ω) ∩L2 0(Ω); ∀K∈ Th, qh|K∈P1(K)}. Note that the space Mhcoincides with the space M1 hdefined in (13). Moreover it is well-known that the inf-sup condition (11) holds for these spaces. We take the space Xhequal to the space X1 hintroduced in (12), namely Xh={vh∈H1 0(Ω)d;∀K∈ Th,vh|K∈P1(K)d}, and finally the operator Πhequal to the Lagrange interpolation operator at all vertices of elements of Thwhich are not on the boundary ∂Ωwith values in Xh. Nothing more is needed to implement the discrete problem (7). As Xhcontains (in fact, is equal to) the space X1 h, owing to the error estimate (4.10) we only use the estimators ηKto perform the grid adaptations. We have decided to present two numerical experiments. The first one deals with a flow with known smooth solution in order to test the efficiency of the error estimators ηKand the ability of the grid adaptation process to obtain accurate solutions with reduced computational time. The second one involves a more realistic case (the step flow) to check the ability of the grid adaptation process to accurately solve a flow with low smoothness at large Reynolds number. In both cases we consider laminar flows, the application to fully turbulent flows is in progress. 5.1 Case of a given solution This test analyzes the efficiency of the error indicators ηKdefined by (24). We check whether the error ∥u−uh∥H1(Ω)dis proportional to the Hilbertian sum of the indicators (∑ K∈Th η2 K)1/2 asymptotically as h→0for small values of the laminar diffusion ν. We consider the square domain Ω =]0, π[×]0, π[, and the pair (velocity, pressure) given by u(x, y) = (eysin(x),−eycos(x)), p(x, y) = ex+y+x2+c, where the constant cis set to have a zero-mean pressure. The pair (u, p)is solution of the Navier–Stokes equations (1) with an appropriate r.h.s. fand non-homogeneous Dirichlet boundary conditions. The solution of the discrete problem (7) has been computed through a single fixed point iteration based upon linearization of this problem. This simple procedure Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 69 converges without difficulty for all values of νtested, ν= 10−1,ν= 10−2,ν= 10−3 and ν= 2 ×10−4. The initialization is the solution of the Stokes problem with the same ν. The procedure is assumed to have converged when the relative error between the L2(Ω)-norm of two consecutive iterates is smaller than 10−8. We have used nonstructured meshes to avoid possible super-convergence effects. We first check the efficiency of the error indicator by computing the efficiency index, i.e., the ratio of the relative error to the Hilbertian sum of the indicators. Tables 1 to 4 display these results corresponding to the four values of the Reynolds numbers here equal to 1 ν, respectively, for increasing numbers of degrees of freedom. Degrees of freedom CPU Time Relative Error Sum Indicators Efficiency 1192 18.44 0.076481 19.4276 0.003937 4460 68.05 0.031235 4.04622 0.007720 9897 125.05 0.016780 1.931050 0.008690 17728 195.66 0.010158 1.030970 0.009853 27305 325.99 0.006380 0.641225 0.009950 39159 440.47 0.005159 0.479480 0.010761 52534 573.74 0.004103 0.368128 0.011147 68708 685.35 0.003551 0.293501 0.012100 85665 768.25 0.002753 0.241493 0.011398 108184 967.91 0.002078 0.183102 0.011347 Table 1: Efficiency for Re = 10 Degrees of freedom CPU Time Relative Error Sum Indicators Efficiency 1192 25.62 0.137147 23.2329 0.005903 4460 105.19 0.076050 5.269140 0.014433 9897 185.15 0.047807 2.853680 0.016753 17728 350.90 0.035531 1.781520 0.019944 27305 517.89 0.025477 1.169900 0.021777 39159 696.21 0.020265 0.893653 0.022676 52534 879.45 0.015772 0.700153 0.022526 68708 1123.80 0.014649 0.620200 0.023620 85665 1500.64 0.011526 0.491341 0.023458 108184 1790.11 0.009919 0.423223 0.023437 Table 2: Efficiency for Re = 100 In all cases the efficiency index tends to a constant as the number of degrees of freedom increases. This is made apparent in Figure 1, that displays the efficiency index for all cases considered. We also have compared the CPU time required to achieve an error below a certain tolerance for all considered Reynolds numbers, between the direct calculation in a fixed grid and the adaptive calculation. The results are displayed in Table 5. The CPU 70 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol Degrees of freedom CPU Time Relative Error Sum Indicators Efficiency 1192 26.65 0.151984 23.9895 0.006335 4460 113.71 0.100590 5.990380 0.016792 9897 210.73 0.070909 3.692480 0.019204 17728 396.32 0.060863 2.759200 0.022058 27305 588.75 0.052434 2.190070 0.023942 39159 853.88 0.044483 1.791170 0.024835 52534 1086.01 0.036521 1.454960 0.025101 68708 1389.95 0.034326 1.351890 0.025391 85665 1864.85 0.028890 1.184070 0.024399 108184 2168.61 0.026960 1.063480 0.025350 Table 3: Efficiency for Re = 1000 Degrees of freedom CPU Time Relative Error Sum Indicators Efficiency 1192 27.23 0.153598 24.0634 0.006383 4460 121.41 0.104726 6.089710 0.017197 9897 225.29 0.076238 3.858290 0.019759 17728 430.01 0.066984 2.968460 0.022565 27305 684.72 0.061366 2.504370 0.024504 39159 965.48 0.053840 2.197200 0.024504 52534 1366.51 0.046629 1.805160 0.025831 68708 1611.31 0.044679 1.716120 0.026035 85665 2248.91 0.040041 1.520500 0.026334 108184 2962.29 0.039491 1.580815 0.024981 Table 4: Efficiency for Re = 5000 Re ν Relative Error CPU Time Relative Error CPU Time Gain direct direct with with adaptation adaptation 10 1.e-1 0.002078 967.91 0.001785 248.68 3.9 100 1.e-2 0.020265 696.21 0.020538 311.16 2.2 100 1.e-2 0.009919 1790.11 0.008803 469.35 3.8 1000 1.e-3 0.036521 1086.01 0.037762 370.63 2.9 1000 1.e-3 0.026960 2168.61 0.017130 1104 2.0 5000 2.e-4 0.046629 1366.51 0.04766 431.62 3.2 5000 2.e-4 0.039491 2962.29 0.028350 1548.02 1.9 Table 5: Comparison of CPU times, direct versus adaptive solution Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 71 Figure 1: Efficiency of the indicators time gain ranges from 2 to 4, essentially depending on the number of grid adaptations required to achieve the prescribed tolerance. 5.2 A more realistic case In this test we analyze the performance of our adaptive strategy to solve a well known, but rather hard to solve, problem. This is the backward-facing step flow. This flow takes place in a non-convex domain, and consequently the velocity and pressure have a low accuracy. In this flow the challenge is to accurately compute the reattachment length, here denoted by x1, which is the length of the main recirculating region behind the step. We refer to [1] and [10] for a more detailed description of the difficulties linked to this problem. The characteristics of the backward-facing step flow considered in this study are shown in Figure 2, where hl=5.2, H=10.1, h=4.9. The Reynolds number is defined as Re =UD ν, where Uis the inlet mean velocity or in other words two thirds of the maximum inlet velocity and Dis the hydraulic diameter of the inlet channel D= 2hl. To change the Reynolds number, we set the value of the laminar viscosity to ν= 15, and we re-scale the inlet velocity profile. The single fixed point iteration used in Section 5.1 fails here, even for moderate Reynolds numbers. Instead, we have used a Newton method with a continuation strategy with respect to the Reynolds number: The initialization to compute Re = 100 is the solution for Re = 10, and so on for all values considered: Re = 100,Re = 300, 72 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol Figure 2: Geometric characteristics of backward-facing step problem Re Our x1Our x1Ref. [1] (without adapt.) (with adapt.) (experimental) 100 3.01114 3.00439 2.98 300 6.38367 6.66069 6.66 500 9.68886 9.91935 10.10 1000 17.55904 17.70518 17.86 Table 6: Comparison of reattachment points Re = 500 and Re = 1000. For larger values the flow becomes unsteady, up to Re ≃8.000, where the flow becomes fully turbulent. The computation for non-adapted grids also uses a continuation strategy with respect to the Reynolds number, where the initialization for the current Reynolds number is the solution obtained with the adapted grid with the preceding Reynolds number. If instead of this initialization we use the solution with the non-adapted grid for the preceding Reynolds number, the Newton’s method does not converge. Then, the comparison of CPU times is not meaningful. Let us say, anyhow, that these times are similar, but the precision obtained with the adaptive procedure is better than the one obtained with the direct solution (see Table 6). 6 Conclusions In this paper we have developed the a posteriori error analysis of a sub-grid eddy viscosity - VMS method. This is a method formally with two grids that only needs a Figure 3: Streamlines of the adaptive solution for Re = 1000 Error analysis of a subgrid eddy viscosity multi-scale discretization of the Navier–Stokes equations 73 grid and an interpolation operator on a virtual coarser grid to be programmed. The subgrid eddy viscosity effects are taken into account by means of the Smagorinsky model with a projection term that filters the action of eddy viscosity on the large resolved scales. We have proved that the standard a posteriori error analysis for Navier–Stokes equations may be extended to the VMS model considered. This mainly arises because the eddy viscosity term is locally Lipschitz in H1norm. We obtain the same family of error indicators as for the Navier–Stokes equations, due to the projection structure of the eddy viscosity term. This error estimator is optimal. A more general error estimator is needed for the standard Smagorinsky turbulence model, that does not seem to be optimal. Our numerical tests yield quite satisfying results for analytic and step flow problems in laminar regime. We obtain a remarkable computing time saving with a very good accuracy. The application to fully turbulent flows is in progress. References [1] B.F. Armaly, F. Durst, J.C.F. Pereira and D. Sch¨ oung, Experimental and theoretical investigation of backward-facing step flow, J. Fluid. Mech 127 (1983), 473–496. [2] C. Bernardi, T. Chac´ on Rebollo, F. Hecht and R. Lewandowski, Automatic insertion of a turbulence model in the finite element discretization of the Navier– Stokes equations, Math. Models and Methods in Applied Sciences 19 (2009), 1139–1183. [3] C. Bernardi, F. Coquel and P.-A. Raviart, Automatic coupling and finite element discretization of the Navier–Stokes and heat equations, submitted. [4] C. Bernardi, Y. Maday and F. Rapetti, Discr´ etisations variationnelles de probl` emes aux limites elliptiques, Collection “Math´ ematiques et Applications” 45, Springer-Verlag (2004). [5] L.C. Berselli, T. Iliescu and W.J. Layton, Mathematics of Large Eddy Simulation of Turbulent Flows, Springer (2006). [6] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer- Verlag (1991). [7] F. Brezzi, J. Rappaz and P.-A. Raviart, Finite dimensional approximation of nonlinear problems, Part I: Branches of nonsingular solutions, Numer. Math. 36 (1980), 1–25. [8] T. Chac´ on Rebollo and R. Lewandowski, Mathematical and Numerical Foundations of Turbulence Models and Applications, in preparation. [9] M. Dauge, Stationary Stokes and Navier–Stokes systems on two- or threedimensional domains with corners I, SIAM Journal on Mathematical Analysis 20 (1989), 74–97. 74 C. Bernardi, T. Chac´ on Rebollo, M. G´ omez M´ armol [10] E. Erturk, Numerical solutions of 2-D steady incompressible flow over a backward-facing step, Part I: High Reynolds number solutions, Computer & Fluids 37 (2008), 633–655. [11] M. Germano, Differential filters for the large eddy numerical simulation of turbulent flows, Phys. Fluids 29(6) (1986), 1755–1757. [12] M. Germano, Differential filters of elliptic type, Phys.Fluids 29(6) (1986), 1757– 1758. [13] V. Girault and P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations, Theory and Algorithms, Springer–Verlag (1986). [14] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman (1985). [15] F. Hecht and O. Pironneau, FreeFem++, second edition, v 3.15, Universit´ e Pierre et Marie Curie, Paris (2011), on the web at http://www.freefem.org/ff++/ftp/freefem++doc.pdf. [16] T. J. R. Hughes, G. Feij´ oo, L. Mazzei and J-B Quincy, The variational multiscale method: a paradigm for computational mechanics, Comp. Meth. Appl. Mech. Engrg. 166 (1-2) (1998), 3–24. [17] T. J. R. Hughes, L. Mazzei and K. E. Jansen, Large eddy simulation and the variational multiscale method, Comput. Vis. Sci. 3(2000), 47–59. [18] T. J. R. Hughes, L. Mazzei, A. Oberai and A. Wray, The multiscale formulation of large eddy simulation, Phys. Fluids 13(2) (2001), 505–512. [19] T. J. R. Hughes, A. Oberai and L. Mazzei, Large eddy simulation of turbulent channel flows by the variational multiscale method, Phys. Fluids 13(6) (2001), 1784–1799. [20] D. K. Lilly, A proposed modification of the Germano subgrid closure method, Physics of Fluids A,4(3) (1992), 633–635. [21] J. Pousin and J. Rappaz, Consistency, stability, a priori and a posteriori errors for Petrov-Galerkin methods applied to nonlinear problems, Numer. Math. 69 (1994), 213–231. [22] J. Smagorinsky, General circulation experiment with the primitive equations. I. The basic experiment, Monthly Weather Review 91 (1963), 99–164. [23] R. Temam, Navier–Stokes Equations, AMS Chelsea (2001). [24] E. R. VanDriest, On turbulent flow near a wall, J. Aerosp. Sci. 23 (1956),1007. [25] R. Verf¨ urth, A Review of A Posteriori Error Estimation and Adaptive Mesh- Refinement Techniques, Wiley & Teubner (1996).