Numerical analysis of a finite element projection-based VMS turbulence model with wall laws

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ Esta es la versión aceptada del artículo publicado en: This is an accepted manuscript of a paper published in: Computer Methods in Applied Mechanics and Engineering (2014): 13/11/2014 DOI: Copyright: El acceso a la versión publicada del artículo puede requerir la suscripción de la revista. Access to the published version may require subscription. “This is an Accepted Manuscript of an article published by Elsevier in Computer Methods in Applied Mechanics and Engineering on 13 November 2014, available at: http://dx.doi.org/10.1016/j.cma.2014.11.023” Accepted Manuscript Numerical analysis of a finite element projection-based VMS turbulence model with wall laws Tom´ as Chac´ on Rebollo, Macarena G´ omez M´ armol, Samuele Rubino PII: S0045-7825(14)00449-6 DOI: http://dx.doi.org/10.1016/j.cma.2014.11.023 Reference: CMA 10473 To appear in: Comput. Methods Appl. Mech. Engrg. Received date: 26 June 2014 Revised date: 9 November 2014 Accepted date: 13 November 2014 Please cite this article as: T. Chac´ on Rebollo, M. G´ omez M´ armol, S. Rubino, Numerical analysis of a finite element projection-based VMS turbulence model with wall laws, Comput. Methods Appl. Mech. Engrg. (2014), http://dx.doi.org/10.1016/j.cma.2014.11.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. • We develop a stabilized projection-based FE-VMS turbulence model including wall laws. • We perform the numerical analysis of the model (stability, convergence, error estimates). • We prove the asymptotic energy balance of the system for slightly smooth flows. • Smooth flows are solved with optimal accuracy. • Good accuracy is obtained with benchmark turbulent flow problems on coarse grids. Highlights (for review) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Numerical analysis of a finite element projection-based VMS turbulence model with wall laws Tom´as Chac´on Rebollo ∗ , Macarena G´omez M´armol † , Samuele Rubino ‡ November 9, 2014 Abstract This paper deals with the numerical analysis of a finite element projection-based VMS turbulence model that includes general non-linear wall laws. Only a single mesh and interpolation operators on a virtual coarser mesh are needed to implement the model. We include a projection-stabilization of pressure to use the same interpolation for velocity and pressure. Good accuracy is obtained with benchmark turbulent flow problems on coarse grids, that justify the interest of this approach. Also, the model solves smooth flows with optimal accuracy. Keywords: Variational Multi-Scale Methods, Interpolation Operators, Wall Laws, Turbulence Modeling, Stabilized Methods, Numerical Analysis. 2010 Mathematics Subject Classification: 76D05, 76F06, 76N20, 65M60, 65M12. 1 Introduction This paper deals with the numerical analysis of the approximation of incompressible flows in turbulent regime by means of the Variational Multi-Scale (VMS) models. The VMS models are increasingly used as a valid alternative to Large Eddy Simulation (LES) models, as they provide a similar accuracy and avoid some drawbacks. Indeed, no commutation error between the variational projection and differential operators occurs, that arises for LES models because the averaged/filtered equations do not satisfy the boundary conditions (Cf. [38,39]). On the contrary, VMS are intrinsically discrete models, and no approximation of an intermediate averaged model is needed. We focus here on the projection-based VMS-LES turbulence models (Cf. [33,34,35,36]). These are three-level methods with large, sub-filter (resolved) scales and small un-resolved scales. The multi-scale setting clarifies the use of sub-grid eddy viscosity to model the viscosity interaction between the sub-filter scales and the small un-resolved scales. In particular, we will address a multi-scale Smagorinsky modeling of the eddy viscosity, which contains the restriction to the sub-filter scales through a projection/interpolation operator. The simulation of wall-bounded flows with VMS models, however, may become very expensive in terms of computational resources due to the computation of boundary layers, as this requires very fine meshes in the normal direction to the wall (see, for instance, John ∗Dpto. EDAN & IMUS, Universidad de Sevilla, C/Tarfia, s/n. 41012 Sevilla, Spain [email protected] †Dpto. EDAN, Universidad de Sevilla, C/Tarfia, s/n. 41012 Sevilla, Spain [email protected] ‡Dpto. EDAN & IMUS, Universidad de Sevilla, C/Tarfia, s/n. 41012 Sevilla, Spain [email protected] 1 *Manuscript Click here to download Manuscript: Steady_VMS-SGM_8.tex Click here to view linked References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 and Kindl [41], Bazilevs et al. [2]). A way to overcome this difficulty, recently applied to VMS models, is to weakly impose no-slip boundary conditions (Cf. [3,4]). An alternative is the use of wall laws, which in their turn replace the usual no-slip boundary conditions by modeled conditions that set the stress of the flow at some distance from the wall. Wall laws are widely used in engineering simulation of turbulence, usually in RANS models. In this paper, we focus on the combined use of VMS-LES models with general non-linear wall-law boundary conditions, in the context of the Finite Element Method (FEM). The most common wall law is the logarithmic law, introduced by Prandtl in 1925 (Cf. [49]), and derived by similarity laws by Von Karm´an in 1930 (Cf. [58]). The mathematical and numerical analysis of wall laws, as boundary conditions for the Smagorkinsky LES turbulence model, was introduced by Par´es (Cf. [47]). A serious difficulty linked to the approximation of wall laws is the discretization of slip boundary conditions. A strong imposition of these laws only applies to polyhedric domains. The more general treatment, introduced by Verf¨urth (Cf. [56,57]), is by duality. This kind of discretization is quite involved, as it requires a specific boundary finite element space for the multipliers. In [48], Par´es introduces a weak formulation of the slip condition. However, the analysis of these techniques has not been extended to mixed boundary conditions, in which other kinds of boundary conditions are imposed, in combination with wall laws, in different parts of the boundary. This is due to a lack of density results by smooth functions for the functional spaces involved. Here, we will analyze the case of mixed boundary conditions including wall laws, replacing the lacking density result by smooth functions by a similar one with finite element functions, for polyhedric domains. In this paper, we perform the numerical analysis of a finite element projection-based VMS model that only needs a (fine) grid and interpolation operators on a virtual coarser grid. The large scales are represented in the coarse grid, while the sub-filter scales are their complement into the fine grid. The eddy diffusion term has a projection structure to filter out the large scales and let the eddy diffusion act only on the sub-filter resolved scales. We use high-order term-by-term stabilization to stabilize each single term that could lead to unstable discretizations (e.g. convection, pressure gradient), with high accuracy (Cf. [16,17,20]). This allows in particular to use polynomials of the same degree to interpolate velocity and pressure. The used stabilization procedure perfectly fits into the VMS framework. The model includes mixed Dirichlet - wall-law boundary conditions to take into account inflow and solid wall boundaries at the same time. We perform a numerical analysis of this approximation in steady regime. We prove stability and weak convergence for solutions with the natural minimal regularity (u, p)∈ [H1(Ω)]3×L2(Ω). Moreover, we perform an error analysis, which strengthen the fact that the proposed model is suitable both for laminar and turbulent flows. In particular, for diffusion-dominated flows, we recover optimal convergence rates. The error analysis permits to prove the strong convergence of the proposed model for slightly smooth flows, and a subsequent asymptotic energy balance of the system. Finally, we perform numerical tests. On one hand, we show the good numerical performances of the proposed model confirming the theoretical convergence expectations for a 2D smooth steady flow (Test 1). On another hand, we present simulations of an equilibrium 3D channel flow at a friction-velocity Reynolds number Reτ= 180 on coarse grids, to analyze the basic numerical performances of the proposed model applied to the computation of turbulent flows, with and without wall-law boundary conditions (Test 2). Similar error levels are obtained for first and second-order statistics with respect to a residual-based VMS model on an equivalent coarse grid. The use of wall laws permits to maintain a similar good accuracy 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 (at least for the leading stream-wise velocity, as expected), with a significant reduction of the computational cost in comparison to the standard application of no-slip boundary conditions at the walls. The structure of the paper is as follows: In Section 2, we present the variational formulation of the continuous and discrete problems we work with, and we state their main properties. Section 3is devoted to the numerical analysis of the proposed discrete model (stability, convergence, error estimates), and to the study of the asymptotic energy balance of the system. The numerical tests are presented in Section 4. Finally, Section 5states the main conclusions of the paper. 2 The continuous and discrete problems We introduce a mixed boundary value problem for the steady incompressible Navier-Stokes equations, that includes a wall-law boundary condition in combination with an inflow boundary condition. Let Ω be a bounded polyhedric connected domain in Rd,d= 2 or 3, with a Lipschitz boundary split into Γ = ΓD∪Γn. We suppose Γ split as the union of the sides Σ1,...,Σr, that we assume to be closed (d−1)−dimensional sets (straight segments when d= 2 or polygons when d= 3), in such a way that ΓD=Sk−1 i=1 Σi, Γn=Sr i=kΣi, for some integer k∈ {2, . . . , r}. We impose a Dirichlet inflow boundary condition on ΓDand a wall-law boundary condition on Γn. The problem reads: Find u: Ω −→ Rdand p: Ω −→ Rsuch that: (2.1)            ∇·(u⊗u)−2ν∇·D(u) + ∇p=fin Ω, ∇·u= 0 in Ω, −[n·2νD(u)]τ=g(u)τon Γn, u·n= 0 on Γn, u=uDon ΓD, where u⊗uis the tensor function of components uiuj,D(u) is the symmetric deformation tensor given by D(u) = (1/2)(∇u+ (∇u)t), nis the outer normal to Γ, the notation τ represents the tangential component with respect to Γ defined as uτ=u−(u·n)n, and g:Rd→Rdis a given function, that determines the wall law (see below). The unknowns are the velocity uand the pressure pof the incompressible fluid. The data are the source term f, that represents a body force per mass unit (typically the gravity), the kinematic viscosity νof the fluid, that is a positive constant, and the Dirichlet data uD. Remark 2.1. The extension of the analysis performed in this paper to free-normal-tension (also called “do-nothing”) boundary conditions on outflow boundaries is similar to the case of Navier-Stokes equations, with the essential difficulty of defining the proper projection operator to compute error estimates. Remark 2.2. The analysis performed in the paper strongly uses the assumption that Ωis a polyhedric domain, to approximate the slip boundary condition u·n= 0 on Γn. There exist well-established techniques to solve this difficulty for domain with curved boundaries, introduced by Verf¨urth. For instance, the slip condition may be considered as a restriction, and implemented through a saddle-point problem approach (Cf. [57]). Another possible remedy is to use isoparametric finite elements to fit the curved parts of the boundary (Cf. [56]). We do not consider here this situation, to avoid nonessential complexities that have been treated elsewhere. 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 2.1 Variational formulation of the continuous problem We consider the Sobolev spaces Hs(Ω), s∈R, and Wm,p(Ω), m∈N, 1 ≤p≤ ∞, equipped with the standard norms. We denote by ||·||m,p,Ωthe standard norm in Wm,p(Ω). In order to give a variational formulation of problem (2.1), let us consider the velocity space: W(Ω) = nw∈[H1(Ω)]d:w=0on ΓD,w·n= 0 on Γno. This is a closed linear subspace of [H1(Ω)]d, and thus a Hilbert space endowed with the [H1(Ω)]d-norm. Indeed, it is well known that the trace is a continuous mapping from H1(Ω) into L2(Γ). Also, the normal component of wexists a.e. on Γ, as n∈[L∞(Γ)]d, and since w∈[H1(Ω)]d, by trace theorem and Sobolev injection (Cf. [9]), its trace on Γ belongs to [L4(Γ)]d. Then, the mapping w∈W(Ω) 7→ w·n∈L4(Γn) is continuous, as: kw·nk0,4,Γn≤ kwk0,4,Γn≤Ckwk1,2,Ω. Thanks to Korn’s inequalities (Cf. [32]), the [H1(Ω)]d-norm is equivalent on W(Ω) to the norm: kwkW(Ω) =kD(w)k0,2,Ω. We assume the Dirichlet boundary data admissible, in the sense that there exists a divergence-free lifting UD∈[H1(Ω)]dsuch that UD|ΓD=uDand UD= 0 on Γn. Such a lifting exists if uD∈[H1/2 00 (ΓD)]d, as we assume hereafter. In this way, we search for a continuous solution u=u0+UD, with u0∈W(Ω) divergence-free. We shall consider the following variational formulation of (2.1): Find (u, p)∈[UD+W(Ω)] ×L2 0(Ω) such that: (2.2) b(u;u,v) + a(u,v)−(p, ∇·v)Ω+hG(u),vi=hf,vi, (∇·u, q)Ω= 0, for any (v, q)∈W(Ω) ×L2 0(Ω), where h·,·i stands for the duality pairing between W(Ω) and its dual [W(Ω)]0. The forms b,aand Gare given by: b(w;u,v) = 1 2[(w·∇u,v)Ω−(w·∇v,u)Ω],(2.3) a(u,v)=2ν(D(u), D(v))Ω,(2.4) hG(u),vi= (g(u),v)Γn,(2.5) for u,v,w∈[H1(Ω)]d. Semicolons (; ) are used for forms that are non-linear with respect to its first argument. The function gis given in implicit form as: g(u) =      u |u|(uτ)2if |u|>0, 0if |u|= 0, where uτ=uτ(|u|) is the wall-friction velocity, computed as unique solution of the algebraic equation: (2.6) u+=L(y+),with u+=|u| uτ and y+=uτy ν. Here, u+is a friction non-dimensional velocity, Lis the wall-law function, obtained from an asymptotic analysis in the boundary layer, y+denotes a friction non-dimensional normal 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 distance to the solid wall, and ydenotes the normal distance to the solid wall. We suppose that the boundary layer is divided into two sub-layers (Cf. [21]): T+ 1= Γn×[0, y+ 0], T+ 2= Γn×[y+ 0, A+], where y+ 0denotes a fixed friction non-dimensional normal distance to the solid wall. The most common wall-law function is the logarithmic law of Prandtl [49] and Von K´arm´an [58]: (2.7) L(y+) =      y+if y+∈[0, y+ 0], 1 C1 log(y+) + C2if y+∈[y+ 0, A+], where C1≃0.41 and C2≃5.5 are constants, calculated from experimental measurements, and y+ 0is chosen by preserving the continuity of L(y+ 0≃11.5). The law (2.7) does not take into account the transition zone between the viscous and logarithmic sub-layer, called the buffer layer. Actually, there exist other several possible settings of L(e.g., the Spalding’s wall law [54]) which model the three boundary sub-layers by a single formula. In all cases, the wall-law function Lis non-negative, strictly increasing and continuous, L0admits a finite number of discontinuities, and there exist two positive constants K1and K2such that: (2.8) lim z+→0+ L(z+) z+=K1,lim z+→∞ L(z+) log z+=K2. This ensures that the associated mapping Gis well defined from W(Ω) into its dual (Cf. Par´es [47]). 2.2 Finite element spaces This section focuses on the construction of finite element (FE) spaces that approximate the slip condition u·n= 0 on Γn. Let {Th}h>0be a family of affine-equivalent and conforming (i.e., without hanging nodes) triangulations of Ω, formed by triangles or quadrilaterals (d= 2), tetrahedra or hexaedra (d= 3). We shall assume that the family of triangulations {Th}h>0is also admissible in the following sense: Definition 2.3. The family of triangulations {Th}h>0is admissible if ΓDand Γnare the union of whole sides of elements of Th. Given an integer l≥0, and an element K∈ Th, denote by Rl(K) either Pl(K) (i.e., the space of Lagrange polynomials of degree ≤l, defined on K), if the grids are formed by triangles (d= 2) or tetrahedra (d= 3), or Ql(K) (i.e., the space of Lagrange polynomials of degree ≤lon each variable, defined on K), if the family of triangulations is formed by quadrilaterals (d= 2) or hexaedra (d= 3). We consider the following FE spaces for the velocity: (2.9)              Yl h=Vl h(Ω) = {vh∈C0(Ω) : vh|K∈Rl(K),∀K∈ Th}, Yl h= [Yl h]d={vh∈[C0(Ω)]d:vh|K∈[Rl(K)]d,∀K∈ Th}, Xh={vh∈Yl h:vh=0on ΓD,vh·ni= 0 on Σi, i =k,...,r} ⊂ Yl h, 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 where niis the outer normal to Σifor i=k,...,r, and we recall that Γn=Sr i=kΣi. Hereafter, Yl h(resp., Yl h) will constitute the discrete foreground vectorial (resp., scalar) spaces in which we will work on. We prove that the family of spaces {Xh}h>0is effectively an internal approximation of W(Ω), i.e. a family of finite-dimensional sub-spaces of W(Ω) such that for any v∈W(Ω): lim h→0d1,2,Ω(v,Xh)=0, where: d1,2,Ω(v,Xh) = inf vh∈Xhkv−vhk1,2,Ω. To do it, let us consider the uniformly stable and convergent Bernardi-Maday-Rapetti (BMR, [6]) interpolation operator IPhfrom [H1(Ω)]don Yl has follows. Let us denote by Ahthe set of Lagrange interpolation nodes for space Yl h. Then: (2.10) IPhv=X α∈Ah vαλα(x) for x∈Ω, where λαare the canonic basis functions of the Lagrange interpolation, given by: λα∈Yl h, λα(β) = δα,β for all α, β ∈ Ah, with δα,β the Kronecker delta and vαan averaged value of vin a neighborhood of node α. Following Chac´on and Lewandowski [22], Sect. A.3, it may be proved that if the family of triangulations is admissible, then the values vαmay be chosen to preserve both the no-slip and slip boundary conditions: If v∈W(Ω), then vα·n|F= 0 for any F∈∂Th(α) if α∈ Ah∩Γn, vα=0if α∈ Ah∩ΓD, where: ∂Th(α) = {F⊂Γ : Fis a side of some element of Thsuch that α∈F}, and n|Fdenotes the outer normal to Ω on F. This permits to prove the following: Lemma 2.4. Assume that the family of triangulations {Th}h>0is admissible. Then, IPhv∈Xhif v∈W(Ω). The proof of this Lemma can be found in [22], so that we omit it for brevity. Lemma 2.4 and the convergence in H1(Ω) of the BMR interpolation operator IPhpermits easily to conclude that the family {Xh}h>0is an internal approximation of W(Ω) for regular triangulations. 2.3 A projection-based VMS turbulence model We approximate the weak formulation (2.2) of the boundary value problem (2.1) for the steady incompressible Navier-Stokes equations by a projection-based eddy viscosity multiscale model. To state it, let us introduce the space: (2.11) Xh={vh∈Yl−1 h:vh=0on ΓD,vh·ni= 0 on Σi, i =k,...,r}, 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 where the function e:R+→R+is defined by: (3.14) e(|u|) =            [uτ(|u|)]2 |u|if |u|>0, 1 (y/ν)K1 if |u|= 0, and K1is given by (2.8). We recall that uτis the unique solution of the implicit equation (2.6). We prove that eis continuous. Indeed, let us re-write equation (2.6) as: |u|=uτL(uτy/ν) = F(uτ). The wall-law function Lis strictly increasing and continuous, then Fis strictly increasing and continuous. Due to (2.8), lim uτ→0+F(uτ) = 0 and lim uτ→∞ F(uτ) = +∞. It implies that uτ=F−1(|u|) is continuous, and so eis continuous at any |u| 6= 0. To prove the continuity at |u|= 0, we set t=uτ(|u|). Then, |u|=tL(λt), with λ=y/ν, and: lim |u|→0e(|u|) = lim t→0 t2 tL(λt)= lim t→0 t L(λt)=1 λK1 , where the last identity follows from (2.8). We deduce that e(|wh|) is bounded on Ω: max x∈Ω e(|wh(x)|)≤M(kwhk0,∞,Ω),where M(r) = max x∈[0,r]e(x). The estimate for the form b Gbecomes: hb G(u0h),vhi ≤ M(kwhk0,∞,Ω)ku0hk0,2,Γnkvhk0,2,Γn ≤CM(kwhk0,∞,Ω)kD(u0h)k0,2,ΩkD(vh)k0,2,Ω.(3.15) Step 2: Existence of solution of problem (3.13). Problem (3.13) is equivalent to a linear system with dim(Xh) + dim(Mh) unknowns and equations. Then, the existence of solutions is equivalent to its uniqueness. To prove uniqueness, let us assume that problem (3.13) admits a solution (u0h, ph)∈Xh×Mh, that we next estimate in terms of the data. •Velocity estimate. Take vh=u0hand qh=phin (3.13) as test functions. This yields: 2νkD(u0h)k2 0,2,Ω+c0(wh,u0h,u0h)+bsconv(u0h,u0h)+hb G(u0h),u0hi+spres(ph, ph) = hef,u0hi, as the form bis antisymmetric. Since the forms c0and b Gare non-negative, we have: (3.16) 2νkD(u0h)k2 0,2,Ω+bsconv(u0h,u0h) + spres(ph, ph)≤ hef,u0hi. From the boundedness of the forms band ain [H1(Ω)]d, and (3.5), it follows: hef,u0hi ≤ (kfk[W(Ω)]0+δ)kD(u0h)k0,2,Ω+C0RkD(UD)k0,2,ΩkD(u0h)k0,2,Ω +C1[diam(Ω)]2−d/2(R+kD(UD)k0,2,Ω)kD(UD)k0,2,Ω+|bsconv(UDh,u0h)|, 13 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 where δ= (C0kD(UD)k0,2,Ω+ 2ν)kD(UD)k0,2,Ω, and we have used kD(UDh)k0,2,Ω≤ CkD(UD)k0,2,Ω. By applying Hypothesis 2.5 and Young’s inequality, we have: |bsconv(UDh,u0h)| ≤ C2 X K∈Th h2 K 2kσ∗ h(wh·∇UDh)k2 0,2,K +1 2bsconv(u0h,u0h) = I+II. The second summand II passes to the left hand side of inequality (3.16). Using the (L2-)local stability property of the interpolation operator σhon the first summand I, we obtain: I≤C3X K∈Th h2 K(kw0h·∇UDhk2 0,2,ωK+kUDh ·∇UDhk2 0,2,ωK) ≤C3X K∈Th h2 K(kw0hk2 0,∞,ωK+kUDhk2 0,∞,ωK)kD(UDh)k2 0,2,ωK, where we have used H¨older’s inequality in the last line. By a local inverse estimate (Cf. [6]), the regularity of the grid, and the Sobolev embedding theorem, we can write: I≤C4X K∈Th h2(1−d/4) K(kw0hk2 0,4,ωK+kUDhk2 0,4,ωK)kD(UDh)k2 0,2,ωK ≤C5h2(1−d/4)(R2+kD(UD)k2 0,2,Ω)kD(UD)k2 0,2,Ω≤M(R), with M(R) = C5[diam(Ω)]2(1−d/4)(R2+kD(UD)k2 0,2,Ω)kD(UD)k2 0,2,Ω. Collecting all these estimates, from (3.16) we obtain: 2νkD(u0h)k2 0,2,Ω+1 2bsconv(u0h,u0h) + spres(ph, ph) ≤hkfk[W(Ω)]0+e δ(R)ikD(u0h)k0,2,Ω+M(R), where we have denoted: e δ(R) = δ+{C0R+C1[diam(Ω)]2−d/2(R+kD(UD)k0,2,Ω)}kD(UD)k0,2,Ω. Using Young’s inequality: νkD(u0h)k2 0,2,Ω+1 2bsconv(u0h,u0h) + spres(ph, ph)≤hkfk[W(Ω)]0+e δ(R)i2 4ν+M(R). Finally, we conclude: kD(u0h)k0,2,Ω≤hkfk[W(Ω)]0+e δ(R)i 2ν+rM(R) ν;(3.17) kσ∗ h(wh·∇u0h)kτν≤hkfk[W(Ω)]0+e δ(R)i √2ν+p2M(R);(3.18) kσ∗ h(∇ph)kτp≤hkfk[W(Ω)]0+e δ(R)i 2√ν+pM(R).(3.19) Observe that the bound (3.17) may be written as: (3.20) kD(u0h)k0,2,Ω≤ K+A 2kD(UD)k0,2,ΩR, 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 with: (3.21) K=kfkW(Ω)0+δ 2ν+C1+ 2√C5ν[diam(Ω)]2−d/2 2νkD(UD)k2 0,2,Ω, and: (3.22) A=C0+C1+ 2√C5ν[diam(Ω)]2−d/2 ν. As we want w0hin the same ball of u0h, we impose R= 2K, so that it must be: (3.23) kD(u0h)k0,2,Ω≤ K+AkD(UD)k0,2,ΩK ≤ 2K, that is satisfied if: kD(UD)k0,2,Ω≤ A−1. •Pressure estimate. Take qh= 0 in (3.13). Using the previous estimates, this yields: |(ph,∇·vh)Ω| ≤ |−hef,vhi+b(wh,u0h,vh) + a(u0h,vh) + c0(wh,u0h,vh) +bsconv(u0h,vh) + hb G(u0h),vhi| ≤Cnkefk[W(Ω)]0+h1 + h2−d/2+ 2ν+M(kwhk0,∞,Ω)ikD(u0h)k0,2,ΩokD(vh)k0,2,Ω ≤Cnkefk[W(Ω)]0+ [1 + ν+M(kwhk0,∞,Ω)] KokD(vh)k0,2,Ω, where the second inequality follows from the boundedness of the forms band ain [H1(Ω)]d, and we have used properties (3.5), (3.7) and (3.15) to estimate respectively the forms c0, bsconv and b G, while the last inequality follows from (3.23). By the discrete inf-sup condition (3.3) and (3.19), we obtain: kphk0,2,Ω≤C sup vh∈Xh |(∇·vh, ph)Ω| kD(vh)k0,2,Ω +kσ∗ h(∇ph)kτp! ≤Cnkefk[W(Ω)]0+1 + ν+√ν+M(kwhk0,∞,Ω)Ko. Step 3: Existence of solution of problem (2.15). We use Brouwer fixed point theorem to prove existence of solution of problem (2.15). Let us define the mapping F:Xh7→ Xh, that transforms w0h∈Xhinto u0h, the unique solution of problem (3.13). The previous estimates and the uniqueness of solution of (3.13) allow to prove the continuity of F. Indeed, let {wn}n∈N⊂Xhbe a sequence convergent to w∈[H1(Ω)]d. Let us consider the sequence of images {F(wn)}n∈N={un}n∈N⊂Xh. The previous estimates allow to extract a sub-sequence, that we denote in the same way, strongly convergent (as Xhis of finite dimension) to t∈Xh. Also, there exists a subsequence {pn}n∈N⊂Mhstrongly convergent to p∈Mh. We take the limit n→+∞ in: BT(wn,(un, pn),(vh, qh)) = hef,vhi,∀(vh, qh)∈Xh×Mh. 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Note that, as we are working in finite dimension, BTis a continuous function, as all terms appearing in the definition of BTare continuous functions of the d.o.f. of their arguments. We can conclude that: BT(w,(t, p),(vh, qh)) = hef,vhi,∀(vh, qh)∈Xh×Mh. Thus, F(w) = t. As the limit problem satisfied by tadmits a unique solution, then the whole sequence {F(wn)}n∈Nconverges to it, by reductio ad absurdum. The previous estimates (3.23) also prove that Ftransforms the closed ball BXh(0,2K) into itself. Then, by Brouwer fixed point theorem, the mapping Fadmits a fixed point. This fixed point is a solution of (2.15), satisfying the estimates (3.10), (3.11). Step 4: Estimate (3.12). It remains just to prove that the estimate for kphk0,2,Ωdoes not depend on hfor the nonlinear problem (2.15). Note that the dependence on hin the linear case is due to estimate (3.15) of the form b G, where Mdepends on kwhk0,∞,Ω, and thus on h. Instead, in the non-linear case, we can directly apply (3.8) to estimate the form G, and thus the pressure estimate does not depend on h.  Remark 3.7. The estimate (3.11) for the convective and pressure stabilizing terms guarantees an extra-control on the high frequencies of the convective derivative and pressure gradient, which is not obtained by standard projection-based VMS methods (Cf. [41]), for which only the sub-grid eddy viscosity term of Smagorinsky type is added to the standard Galerkin discretization. 3.3 Convergence analysis The convergence analysis is based upon the theory developed in [15], that enables to extend to stabilized methods the standard techniques for the numerical analysis of mixed methods. We shall need a technical result that allows to represent formulation (2.2) as an internal approximation of an “augmented” variational formulation. Definition 3.8. A FE space Zh, constructed on a triangulation Th, is called a bubble FE space if, for all bh∈Zh, for all K∈ Th,bh∈H1 0(K). A similar definition applies for vectorial bubble FE spaces. Lemma 3.9. There exists a family {Zh}h>0of bubble FE sub-spaces of [H1 0(Ω)]dand a family {Sh}h>0of bilinear uniformly continuous and uniformly coercive forms on [H1 0(Ω)]d such that: (3.24) spres(ph, qh) = Sh(Rh(σ∗ h(∇ph)),Rh(σ∗ h(∇qh))),∀ph, qh∈Mh, where Rh: [H−1(Ω)]d−→ Zhis the “static condensation” operator on Zhdefined as follows: Given ϕ∈[H−1(Ω)]d,Rh(ϕ)is the only element of Zhthat satisfies: Sh(Rh(ϕ), zh) = hϕ, zhi,∀zh∈Zh. This result is proved in [15]. We shall also need the following property of bubble FE spaces (Cf. [15]): 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Lemma 3.10. If a sequence {Zh}h>0of bubble FE sub-spaces of [H1 0(Ω)]dis uniformly bounded in [H1 0(Ω)]d, then it weakly converges to zero in [H1 0(Ω)]d. We now state the weak convergence of solutions provided by method (2.15) to a weak solution of the Navier-Stokes boundary value problem (2.2). Theorem 3.11. Assume that Hypothesis 2.5 holds. Then, the sequence {(uh, ph)}h>0 of solutions of the VMS-S approximation (2.15) contains a sub-sequence which is weakly convergent in [H1(Ω)]d×L2(Ω) to a solution of the steady Navier-Stokes equations (2.2). If this solution is unique, then the whole sequence converges to it. Proof. The proof is divided into various steps. Step 1: Extracting sub-sequences. Due to estimates (3.10)-(3.12), the sequence {(u0h, ph)}h>0is uniformly bounded in the space W(Ω) ×L2 0(Ω), which is a Hilbert space. Then, it contains a sub-sequence, that we still denote in the same way, weakly convergent in that space to some pair (u0, p). As the injection of H1(Ω) in Lq(Ω) is compact for 1 ≤q < q∗= 2d/(d−2), by the Rellich-Kondrachov compactness theorem (Cf. [9]) we may assume that the sub-sequence {u0h}h>0is strongly convergent in [Lq(Ω)]d, and so, in particular, in [L4(Ω)]d. Also, the operator Gis compact from W(Ω) to its dual [W(Ω)]0, by Lemma 3.5. Then, we may assume that the sub-sequence {G(u0h)}h>0is strongly convergent in [W(Ω)]0. We recall that UDh is the Stokes projection of UDon Yl+1 h, thus it strongly converges to UDin [H1(Ω)]d. Let us prove that (u, p) is a solution of problem (2.2). Step 2: Limit of convection terms. Let us consider a pair of test functions (v, q) such that v∈W(Ω), q∈ D(Ω) ∩L2 0(Ω), where D(Ω) is the space of C∞ 0(Ω)−functions (i.e., smooth functions with compact support in Ω), that is dense in L2(Ω) (Cf. [44]). As Xhis an internal approximation of W(Ω) (see Subsection 2.2), then there exists a sequence {vh}h>0∈Xhstrongly convergent to vin W(Ω). Moreover, as Mhis an internal approximation of H1(Ω) ∩L2 0(Ω), then there exists a sequence {qh}h>0∈Mhstrongly convergent to qin particular in L2 0(Ω). Then, it is standard to prove that: lim h→0b(uh;uh,vh) = b(u;u,v). Moreover, by (3.7), for the convection stabilizing term we have: |sconv(uh;uh,vh)| ≤ Ch2−d/2kD(uh)k3 0,2,ΩkD(vh)k0,2,Ω. Since the sequences {uh}h>0and {vh}h>0are bounded in [H1(Ω)]d, we deduce: lim h→0sconv(uh;uh,vh)=0. Step 3: Limit of diffusion terms. As ais bilinear and continuous: lim h→0a(uh,vh) = a(u,v). Next, since the sequences {uh}h>0and {vh}h>0are bounded in [H1(Ω)]d, we deduce from Lemma 3.3: lim h→0c0(uh;uh,vh)=0. 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Step 4: Limit of pressure terms. Since {∇·uh}h>0is weakly convergent in L2(Ω) to ∇·uand {qh}h>0is strongly convergent in L2(Ω) to q: lim h→0(∇·uh, qh)Ω= (∇·u, q)Ω. Also, we obviously have: lim h→0(ph,∇·vh)Ω= (p, ∇·v)Ω. The pressure stabilizing term also vanishes in the limit. To prove this, we use the representation formula (3.24). This yields: kσ∗ h(∇ph)k2 τp=spres,h(ph, ph) = Sh(Rh(σ∗ h(∇ph)),Rh(σ∗ h(∇ph))) ≥µSkRh(σ∗ h(∇ph))k2 H1 0(Ω), using the uniform coercivity of the forms Sh. Then, using the fact that kσ∗ h(∇ph)kτp≤ 2√νK, we deduce that the sequence {Rh(σ∗ h(∇ph))}h>0is uniformly bounded in [H1 0(Ω)]d and, by Lemma 3.10, as: spres(ph, qh) = Sh(Rh(σ∗ h(∇ph)),Rh(σ∗ h(∇qh))) = hσ∗ h(∇qh),Rh(σ∗ h(∇ph))i, we conclude: lim h→0spres(ph, qh)=0. Consequently, the pair (u, p) is a weak solution of the Navier-Stokes equations (2.2). Step 7: Uniqueness. As this weak convergence follows from a compactness argument, it is standard to prove, by reductio ad absurdum, that if the limit is unique, then the whole sequence converges to it, in the same weak sense.  Remark 3.12. Following a standard strategy to prove the strong convergence of the velocities for the stabilized VMS-S method (2.15), one can only conclude: lim h→02νkD(u0h)k2 0,2,Ω+kσ∗ h(∇ph)k2 τp= 2νkD(u0)k2 0,2,Ω, so that we cannot achieve a strong convergence result for solutions with natural minimal regularity, neither the corresponding asymptotic energy balance of the system (which is the case, instead, for mixed methods). Remark 3.13. Assuming a slightly increased regularity of the solution for the steady Navier-Stokes equations (2.1) leads to a strong convergence result also in the context of the stabilized VMS-S method (2.15) through an error estimates analysis (See Appendix). This will guarantee the corresponding asymptotic energy balance of the system, as we will state in the next subsection. Also, the error analysis performed in the Appendix proves that the convergence order of the VMS-S method (2.15) is optimal with respect to the polynomial interpolation for laminar smooth flows. 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 3.4 Asymptotic energy balance The error estimates result stated in Theorem 5.3 in Appendix, that implies the strong convergence in particular of the velocities (See Remark 5.4), contains as a sub-product the asymptotic energy balance of VMS-S model (2.15). Indeed, let us define the deformation energy ED, the boundary friction energy EF, the sub-grid eddy dissipation energy ES, and the energy ESC and ESP respectively corresponding to the convection and the pressure stabilizing terms by: ED(u) = a(u,u) = 2νkD(u)k2 0,2,Ω, EF(u) = hG(u),ui=ZΓn g(u)·uds, ES(uh) = c(uh;uh,uh) = 2 X K∈Th (CShK)2ZK|D(uh0)|3dx, ESC(uh) = sconv,h(uh;uh,uh) = kσ∗ h(uh·∇uh)k2 τν, ESP (ph) = spres,h(ph, ph) = kσ∗ h(∇ph)k2 τp. Then, it holds: Corollary 3.14. Let {(uh, ph)}h>0be a sequence of solutions of the VMS-S model (2.15) strongly convergent in [H1(Ω)]d×L2(Ω) to a solution (u, p)of the Navier-Stokes equations (2.2) with regularity [H2(Ω)]d×H1(Ω) (at least). Then: lim h→0[ED(uh) + EF(uh) + ES(uh) + ESC(uh) + ESP (ph)] = ED(u) + EF(u), as: lim h→0ED(uh) = ED(u),lim h→0EF(uh) = EF(u), and: lim h→0ES(uh) = lim h→0ESC(uh) = lim h→0ESP (ph)=0. Thus, the total energy balance is asymptotically maintained in such a way that the deformation energy and the energy dissipated at the wall pass to the limit. In addition, the dissipated eddy energy so as the sub-grid energy due to stabilizing terms asymptotically vanish. 4 Numerical examples We have performed numerical experiments by using the projection-based VMS model (2.15), on one hand to test the theoretical convergence order predicted by the numerical analysis, and stated in Appendix (Theorem 5.3). This is performed in Test 1, where we have considered a 2D smooth steady flow with analytical solution to check the expected optimal convergence order with respect to the polynomial interpolation. On another hand, to analyze the basic numerical performances of the proposed model applied to the computation of turbulent flows, with and without wall-law boundary conditions. This is performed in Test 2, where we address simulations of a fully developed 3D turbulent channel flow at Reτ= 180, for which there exists an extensive literature providing reliable numerical results. 19 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 4.1 Test 1: Convergence order This test is aimed to confirm the optimal convergence order of model (2.15) for laminar smooth flows. We have considered a 2D steady flow in the domain Ω = (0, π)2, whose analytical solution is given by: •Horizontal velocity: u1= 2(sin x1)2sin x2cos x2; •Vertical velocity: u2=−2 sin x1(sin x2)2cos x1; •Pressure: p= cos x1cos x2. The velocity field vanishes on ∂Ω. As physical property, we have taken ν= 10−2. The force term fis computed to match the exact solution. We have used regular meshes with grid size hranging from 0.9 to 0.06, and triangular finite elements of degree 2, 3, and 4 for velocity and pressure (l= 2,3,4 = degree of the polynomial interpolation). We have used setting (2.17) for the multi-scale eddy viscosity term, with nodal interpolation operator Πhtaking values in Yhdefined by (2.12). In this way, only a single mesh is needed for implementation. The stabilization coefficients are given by the adapted Codina’s form (2.19), with constants c1= 4, c2=√c1= 2 (Cf. [26]). In Figure 1, we have plotted the slopes of the estimated error curves for the velocity (in H1−norm) and the pressure (in L2−norm), related to the grid size h. We observe that the convergence rate is optimal, so that we obtain an excellent agreement with respect to the theoretical predictions. 4.2 Test 2: Turbulent channel flow The 3D channel flow is one of the most popular test problems for the investigation of wall bounded turbulent flows. It was pioneered as a LES test problem by Moin and Kim (Cf. [45]), and more recently has been extensively used to test several versions of LES models. Let us mention, for instance, the simulation carried out by Iliescu and Fischer [37] using a Rational LES (RLES) model, in the context of a spectral element code. A number of numerical studies on turbulent channel flow has also been performed applying VMS methods. Among others, Bazilevs et al. [2] tested a Residual-based VMS (RB-VMS) turbulence model on a channel flow, employing linear, quadratic and cubic NURBS (Non- Uniform Rational B-Splines), while John and Kindl [41] compared the performances of two types of VMS methods, where one uses bubble functions to model resolved small scales (Residual Free Bubble-based method), whereas the other one contains the definition of the resolved small scales by an explicit projection in its set of equations (Projection-based method). For the setup of our numerical simulations, we choose to follow the guidelines given by Gravemeier in [30], where different scale-separating operators for a VMS-LES of a turbulent channel flow in the context of a Finite Volume Method (FVM) are tested. As a benchmark, we will use the fine Direct Numerical Simulation (DNS) of Moser, Kim and Mansour [46]. •Setup for numerical simulations. We test model (2.15) with the following different settings of the eddy viscosity term: •SMA model: The Smagorinsky setting, given by c0(uh;uh,vh) = 2(νT(uh)D(uh), D(vh))Ω; 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 10−1 100 10−4 10−3 10−2 10−1 100 h Estimated convergence order P2 h2 u H1−norm p L2−norm 10−1 100 10−6 10−5 10−4 10−3 10−2 10−1 100 h Estimated convergence order P3 h3 u H1−norm p L2−norm 10−1 100 10−8 10−6 10−4 10−2 100 h Estimated convergence order P4 h4 u H1−norm p L2−norm Figure 1: Estimated convergence orders using P2, P3, and P4 FE for velocity and pressure. 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 •VMS-S model: The Small-Small VMS-Smagorinsky setting, given by c0(uh;uh,vh) = 2(νT(u0 h)D(u0 h), D(v0 h))Ω; •VMS-B model: The Berselli-Iliescu-Layton setting of Ref. [7], in which: c0(uh;uh,vh) = 2(νT(e Π∗ hD(uh))e Π∗ hD(uh),e Π∗ hD(vh))Ω, where e Π∗ h=Id −e Πh, and we have denoted by e Πhthe [L2(Ω)]d×d-orthogonal projection on the space: Lh=D(Yh) = {dh∈[L2(Ω)]d×d:dh=D(wh),for some wh∈Yh}. In the numerical experiments, we perform simulations with a slight modification of this method, where we replace the [L2(Ω)]d×d-orthogonal projection e Πhon space Lhby an interpolation operator on a coarser (e.g., P0) FE space, much faster to compute. We use a setup similar to the one of Gravemeier [30]. The computational domain is Ω = (0, L1)×(−δ, δ)×(0, L3), with δ= 1 (wall-normal direction), L1= 2π(stream-wise direction), and L3= (4/3)π(span-wise direction). The boundary conditions are periodic in both the stream-wise and span-wise directions, commonly referred to as homogeneous directions. We perform a comparison between the application of no-slip and wall-law boundary conditions at the walls. The viscosity is ν= 1/180 = 5.5×10−3. The turbulent wall-shear velocity uτ=√τ, where τdenotes the wall-shear stress, and the channel halfwidth δdefine the Reynolds number: Reτ=uτδ ν, besides the kinematic viscosity ν. The Reynolds number based on a unit friction velocity reachable at a steady state is Reτ= 180. Our strategy is as follows: to reach a statistically steady state, we use an evolution approach starting by an initial parabolic velocity profile perturbed by a random velocity fluctuation. We first run a simulation with no-slip boundary conditions at the walls, in order to stabilize uτnear a unitary value. Random velocity fluctuations of 10%-amplitude of the bulk mean stream-wise velocity perturbs the initial condition for the velocity field:    u1(y, t = 0) = u1,c(1 −y2)+0.1u1,mψran, u2(y, t = 0) = 0.1u1,mψran, u3(y, t = 0) = 0.1u1,mψran, where u1,c denotes the stream-wise velocity at the centerline of the channel, u1,m the bulk mean stream-wise velocity, and ψran ∈[−1,1] a random number. We choose u1,c = 25, and hence u1,m =Zδ 0 u1,c(1 −y2)dy = 2u1,c/3 = 16.7. So, the corresponding Reynolds number based on the bulk mean stream-wise velocity (Rem=u1,m2δ/ν, see [43]) is Rem= 6 012. The flow is driven by a constant forcing f= (fp,0,0) = (1,0,0), that models an imposed pressure gradient in the stream-wise direction. The specific choice of a unit value for fpaims at obtaining a unit value for uτin the statistically steady state, subject to the relation uτ=pfpδ(Cf. [29]). We choose to work with Van Driest Damping [55], so that the Smagorinsky constant is changed to the expression CS(1 −exp(−y+/A+)), where 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 y Normalized Reynolds shear stress Reτ = 180 DNS (Moser et al.) VMS−S (NO−SLIP BC) VMS−S (WALL−LAW BC) VMS−B (NO−SLIP BC) SMA (NO−SLIP BC) STAB (NO−SLIP BC) STAB (WALL−LAW BC) Figure 4: Normalized Reynolds shear stress in global coordinates y. Methods e√h˜u2 1i 0e√h˜u2 2i 0e√h˜u2 3i 0eh˜u1˜u2i 0 VMS-S (NO-SLIP BC) 0.2252 0.1652 0.1108 0.1162 VMS-B (NO-SLIP BC) 0.2881 0.2018 0.1246 0.1706 SMA (NO-SLIP BC) 0.3002 0.2236 0.1597 0.1249 STAB (NO-SLIP BC) 0.3781 0.2536 0.1955 0.1708 RB-VMS (Akkerman) 0.5694 0.1753 0.1331 - Table 3: L2-norm of the deviation from the DNS profiles for the second-order statistics. 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure 5: Snapshot of stream-wise velocity contours at t= 5 (top) and t= 10 (bottom). We highlight also some considerations on the computational cost corresponding to the various methods. The computing times related to the statistical period from t= 5 to t= 10 time units (1 250 time steps) for all methods are given in Table 4. These are referred to the sequential execution of the numerical code on a same machine. It can be seen that effectively the use of wall laws for VMS-S method provides a significant reduction of the computing time, of about 34%, while the CPU time of both VMS approaches is comparable with the one of the standard Smagorinsky approach, that saves at most a 4% of computing time. Also, the CPU time for STAB method is comparable with the computing times of the other methods in which the multi-scale Smagorinsky term is added, and again the use of wall laws provides a reduction of the computing time, of about 35%. 5 Conclusions We have developed stabilized projection-based FE-VMS formulations of the incompressible Navier-Stokes equations incorporating wall-law boundary conditions in combination with inflow boundary conditions. We have performed the numerical analysis of the model introduced, proving stability and convergence for flows with minimal regularity. We have replaced the lacking density results by smooth functions by an ad-hoc internal approximation result by finite elements. We also have proved the asymptotic energy balance that holds for flows slightly smoother. Further, we have proved that the model is useful to solve both smooth and turbulent flows. Smooth flows are solved with optimal accuracy, as confirmed by numerical simulations too. The numerical experiments on a turbulent channel flow confirm on one hand that the 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Methods CPU time VMS-S (NO-SLIP BC) 398 207 (134) VMS-S (WALL-LAW BC) 297 524 (100) VMS-B (NO-SLIP BC) 401 128 (135) SMA (NO-SLIP BC) 389 260 (131) STAB (NO-SLIP BC) 369 872 (135) STAB (WALL-LAW BC) 273 714 (100) Table 4: Computing times (in seconds); in parentheses: percentage to computing time with VMS-S (WALL-LAW BC). application of wall-law boundary conditions could provide (at least for the leading streamwise component of the velocity) similar results to those obtained by the standard approach based on the use of no-slip boundary conditions, a refined mesh towards the walls and the Van Driest damping improvement, with a noticeable reduced computational cost. On another hand, they show that the VMS-S method gives quite good results for both first and second-order statistics, in the worst condition of a very coarse basic discretization. Note that when the multi-scale Smagorinsky term is neglected (STAB method), we obtain worse results, even if they are quite close to the ones obtained by adding the sub-grid eddy viscosity. Indeed, STAB method is not strongly consistent, and the addition of the multiscale Smagorinsky term helps to counter-balance the consistency error. The numerical results clarify that the differences are due to the lower numerical diffusion of the STAB method, and that VMS models are more accurate due to the turbulent diffusion term. Summarizing, the proposed VMS-S method provides a good compromise between accuracy and computational complexity, which is an important feature in the context of its practical performances. Appendix Error estimates We derive error estimates for the VMS-S discretization (2.15) for diffusion-dominated flows. The interest of this analysis is to highlight the fact that the VMS-S method proposed may be used to approximate both laminar and turbulent flows. About the former, if they are regular enough, we obtain convergence of optimal order, and the order decreases with the regularity. As a consequence, we also obtain an asymptotic energy balance of the system for slightly smooth flows, as stated in Subsection 3.4 in the body of the paper. Similar error estimates may be obtained in a more general framework, when the solution of Navier-Stokes equations is located in a branch of non-singular Reynolds numbers, in the sense that at these Reynolds numbers there are no bifurcations to more complex flows (Cf. [10,11,12]). We start by setting a condition that ensures the uniqueness of solutions of the steady Navier-Stokes equations (2.15). Theorem 5.1. Assume that: (5.1) 2ν > β(2K+kD(UD)k0,2,Ω), 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 where Kis the quantity defined by (3.21), and: β:= sup z,v,w∈W(Ω) b(z;v,w) kD(z)k0,2,ΩkD(v)k0,2,ΩkD(w)k0,2,Ω . Then, the solution of the steady Navier-Stokes equations (2.2) is unique. Remark 5.2. The condition (5.1) means that the flow is diffusion-dominated. The viscosity νis large enough to balance the convection effects relatively to the data fand UD. A rather standard but quite lengthy and technical error analysis (that we skip for brevity) proves that (Cf. [50]): Theorem 5.3. Assume that Hypothesis 2.5 and estimate (5.1) hold, and that, for smooth enough data, the (unique) solution (u, p)of the steady Navier-Stokes equations (2.15) has augmented regularity, i.e. (u, p)∈[Hs+1(Ω)]d×Hs(Ω),2≤s≤l. Then, the following error estimates for a solution (uh, ph)of the VMS-S method (2.15) hold: (5.2) kD(u−uh)k0,2,Ω+kp−phk0,2,Ω≤Chs+h2−d/2+2η, for some constant Cinversely proportional to νand independent of h, and η≤(l−1) if Πhtakes values in Yhdefined by (2.12), while η≤lif Πhtakes values in Yhdefined by (2.13), ldenoting the degree of the polynomial interpolation. Remark 5.4. We note that, in any case, the convergence order of the VMS-S method (2.15) is optimal with respect to the polynomial interpolation for laminar smooth flows. Decreasing the regularity of the solution (u, p)of the steady Navier-Stokes equations (at most) to [H2(Ω)]d×H1(Ω) implies that the convergence order of the VMS-S method is limited to one, due to the estimate of the pressure stabilizing term. However, this slightly augmented regularity already guarantees the strong convergence of the VMS-S method. Acknowledgements: Research partially supported by the Spanish Government project MTM2012-36124-C02-01. References [1] I. Akkerman. Adaptative variational multiscale formulations using the discrete Germano approach. PhD thesis, Delft University of Technology, 2009. [2] Y. Bazilevs, V. M. Calo, J. A. Cottrell, T. J. R. Hughes, A. Reali, and G. Scovazzi. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Engrg., 197(1-4):173–201, 2007. [3] Y. Bazilevs, C. Michler, V. M. Calo, and T. J. R. Hughes. Weak Dirichlet boundary conditions for wall-bounded turbulent flows. Comput. Methods Appl. Mech. Engrg., 196(49-52):4853–4862, 2007. [4] Y. Bazilevs, C. Michler, V. M. Calo, and T. J. R. Hughes. Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes. Comput. Methods Appl. Mech. Engrg., 199(13- 16):780–790, 2010. 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [5] C. Bernardi, T. Chac´on Rebollo, and D. Yakoubi. Finite element discretization of the Stokes and Navier-Stokes equations with boundary conditions on the pressure. Submitted (http://hal.upmc.fr/hal-00961653), 2014. [6] C. Bernardi, Y. Maday, and F. Rapetti. Discr´etisations variationnelles de probl`emes aux limites elliptiques, volume 45 of Math´ematiques & Applications. Springer-Verlag, 2004. [7] L. C. Berselli, T. Iliescu, and W. J. Layton. Mathematics of large eddy simulation of turbulent flows. Scientific Computation. Springer-Verlag, Berlin, 2006. [8] M. Braack and E. Burman. Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method. SIAM J. Numer. Anal., 43(6):2544–2566, 2006. [9] H. Brezis. Analyse fonctionnelle: Th´eorie et applications. Masson, 1983. [10] F. Brezzi, J. Rappaz, and P.-A. Raviart. Finite-dimensional approximation of nonlinear problems. I. Branches of nonsingular solutions. Numer. Math., 36(1):1–25, 1980/81. [11] F. Brezzi, J. Rappaz, and P.-A. Raviart. Finite-dimensional approximation of nonlinear problems. II. Limit points. Numer. Math., 37(1):1–28, 1981. [12] F. Brezzi, J. Rappaz, and P.-A. Raviart. Finite-dimensional approximation of nonlinear problems. III. Simple bifurcation points. Numer. Math., 38(1):1–30, 1981/82. [13] E. Burman. Interior penalty variational multiscale method for the incompressible Navier-Stokes equation: monitoring artificial dissipation. Comput. Methods Appl. Mech. Engrg., 196(41-44):4045–4058, 2007. [14] T. Chac´on Rebollo. A term by term stabilization algorithm for finite element solution of incompressible flow problems. Numer. Math., 79(2):283–319, 1998. [15] T. Chac´on Rebollo. An analysis technique for stabilized finite element solution of incompressible flows. M2AN Math. Model. Numer. Anal., 35(1):57–89, 2001. [16] T. Chac´on Rebollo, M. G´omez M´armol, V. Girault, and I. S´anchez Mu˜noz. A high order term-by-term stabilization solver for incompressible flow problems. IMA J. Numer. Anal., 33(3):974–1007, 2013. [17] T. Chac´on Rebollo, M. G´omez M´armol, and M. Restelli. Numerical analysis of penalty stabilized finite element discretizations of evolution Navier-Stokes equation. J. Sci. Comput., 61(1):1–28, 2014. [18] T. Chac´on Rebollo, M. G´omez M´armol, and S. Rubino. Derivation of the Smagorinsky model from a Galerkin discretization. In Mascot11 Proc., volume 17 of IMACS Series in Comp. and Appl. Math., pages 61–70, 2013. [19] T. Chac´on Rebollo, M. G´omez M´armol, and S. Rubino. TurboMathS: a FreeFem++ numerical code for VMS turbulence models. Technical report, Research Group “Modelado Matem´atico y Simulaci´on de Sistemas Medioambientales”, University of Seville, 2014. 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [20] T. Chac´on Rebollo, F. Hecht, M. G´omez M´armol, G. Orzetti, and S. Rubino. Numerical approximation of the Smagorinsky turbulence model applied to the primitive equations of the ocean. Math. Comput. Simulation, 99:54–70, 2014. [21] T. Chac´on Rebollo and R. Lewandowski. A variational finite element model for largeeddy simulations of turbulent flows. Chin. Ann. Math. Ser. B, 34(5):667–682, 2013. [22] T. Chac´on Rebollo and R. Lewandowski. Mathematical and numerical foundations of turbulence models and applications. Birkh¨auser, 2014. [23] H. Choi and P. Moin. Effects of the computational time step on numerical solutions of turbulent flow. J. Comput. Phys., 113(1):1–4, 1994. [24] R. Codina. Comparison of some finite element methods for solving the diffusionconvection-reaction equation. Comput. Methods Appl. Mech. Engrg., 156(1-4):185– 210, 1998. [25] R. Codina. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Methods Appl. Mech. Engrg., 190(13- 14):1579–1599, 2000. [26] R. Codina. A stabilized finite element method for generalized stationary incompressible flows. Comput. Methods Appl. Mech. Engrg., 190(20-21):2681–2706, 2001. [27] R. Codina, J. Principe, and S. Badia. Dissipative structure and long term behavior of a finite element approximation of incompressible flows with numerical subgrid scale modeling. In Multiscale methods in computational mechanics, volume 55 of Lect. Notes Appl. Comput. Mech., pages 75–93. Springer, 2011. [28] J. W. Deardorff. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J. Fluid Mech., 41:453–465, 1970. [29] T. Dubois, F. Jauberteau, and R. Temam. Dynamic multilevel methods and the numerical simulation of turbulence. Cambridge University Press, 1999. [30] V. Gravemeier. Scale-separating operators for variational multiscale large eddy simulation of turbulent flows. J. Comput. Phys., 212(2):400–435, 2006. [31] F. Hecht. New development in freefem++. J. Numer. Math., 20(3-4):251–265, 2012. [32] C. O. Horgan. Korn’s inequalities and their applications in continuum mechanics. SIAM Rev., 37(4):491–511, 1995. [33] T. J. R. Hughes, G. R. Feij´oo, L. Mazzei, and J.-B. Quincy. The variational multiscale method—a paradigm for computational mechanics. Comput. Methods Appl. Mech. Engrg., 166(1-2):3–24, 1998. [34] T. J. R. Hughes, L. Mazzei, and K. E. Jansen. Large eddy simulation and the variational multiscale method. Comput. Vis. Sci., 3(1-2):47–59, 2000. [35] T. J. R. Hughes, L. Mazzei, A. A. Oberai, and A. Wray. The multiscale formulation of large eddy simulation: Decay of homogeneous isotropic turbulence. Phys. Fluids, 13(2):505–512, 2001. 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [36] T. J. R. Hughes, A. A. Oberai, and L. Mazzei. Large eddy simulation of turbulent channel flows by the variational multiscale method. Phys. Fluids, 13(6):1784–1799, 2001. [37] T. Iliescu and P. F. Fischer. Large eddy simulation of turbulent channel flows by the rational large eddy simulation model. Phys. Fluids, 15(10):3036–3047, 2003. [38] V. John. On large eddy simulation and variational multiscale methods in the numerical simulation of turbulent incompressible flows. Appl. Math., 51(4):321–353, 2006. [39] V. John and S. Kaya. A finite element variational multiscale method for the Navier- Stokes equations. SIAM J. Sci. Comput., 26(5):1485–1503, 2005. [40] V. John, S. Kaya, and A. Kindl. Finite element error analysis for a projection-based variational multiscale method with nonlinear eddy viscosity. J. Math. Anal. Appl., 344(2):627–641, 2008. [41] V. John and A. Kindl. Numerical studies of finite element variational multiscale methods for turbulent flow simulations. Comput. Methods Appl. Mech. Engrg., 199(13- 16):841–852, 2010. [42] W. P. Jones and M. Wille. Large eddy simulation of a jet in a cross-flow. In 10th Symposium on Turbulent Shear Flows, volume 4, pages 1–6, 1995. [43] J. Kim, P. Moin, and R. Moser. Turbulence statistics in fully developed channel flow at low Reynolds number. J. Fluid Mech., 177:133–166, 1987. [44] J. L. Lions. Quelques m´ethodes de r´esolution des probl`emes aux limites non lin´eaires. Dunod, 2002. [45] P. Moin and J. Kim. Numerical investigation of turbulent channel flow. J. Fluid Mech., 118:341–377, 1982. [46] R. Moser, J. Kim, and N. N. Mansour. Direct numerical simulation of turbulent channel flow up to Reτ= 590. Phys. Fluids, 11(4):943–945, 1999. [47] C. Par´es. Existence, uniqueness and regularity of solution of the equations of a turbulence model for incompressible fluids. Appl. Anal., 43(3-4):245–296, 1992. [48] C. Par´es. Approximation de la solution des ´equations d’un mod`ele de turbulence par une m´ethode de Lagrange-Galerkin. Rev. Mat. Apl., 15(2):63–124, 1994. [49] L. Prandtl. ¨ Uber die ausgebildeten Turbulenz. Zeitschrift f¨ur angewandte Mathematik und Mechanik, 5:136–139, 1925. [50] S. Rubino. Numerical modeling of turbulence by Richardson number-based and VMS models. PhD thesis, Univeristy of Seville, 2014. [51] Y. Saad and M. H. Schultz. GMRES - A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comp., 7(3):856–869, 1986. [52] R. L. Scott and S. Zhang. Finite element interpolation of non-smooth functions satisfying boundary conditions. Math. Comput., 54(190):483–493, 1990. 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 [53] J. Smagorinsky. General circulation experiment with the primitive equations: I. The basic experiment. Mon. Weather Rev., 91(3):99–164, 1963. [54] D. B. Spalding. A single formula for the “law of the wall”. J. Appl. Mech., 28(3):455– 458, 1961. [55] E. R. Van Driest. On turbulent flow near a wall. J. Aerosp. Sci., 23(11):1007–1011, 1956. [56] R. Verf¨urth. Finite element approximation of steady Navier-Stokes equations with mixed boundary conditions. RAIRO Mod´el. Math. Anal. Num´er., 19(3):461–475, 1985. [57] R. Verf¨urth. Finite element approximation of incompressible Navier-Stokes equations with slip boundary condition. Numer. Math., 50(6):697–721, 1987. [58] T. Von K´arm´an. Mechanische ¨ Ahnlichkeit und Turbulenz. Nachr. Ges. Wiss. Gottingen, Math. Phys. Klasse, 58, 1930. 36