A variational finite element model for Large-Eddy simulations of turbulent flows

arXiv:1302.2753v1 [math-ph] 12 Feb 2013 A Variational Finite Element Model for Large-Eddy Simulations of Turbulent Flows Tom´as Chac´on Rebollo ∗ , Roger Lewandowski† February 13, 2013 Abstract We introduce a new Large Eddy Simulation model in a channel, based on the projection on finite element spaces as filtering operation in its variational form, for a given triangulation {Th}h>0. The eddy viscosity is expressed in terms of the friction velocity in the boundary layer due to the wall, and is of a standard sub grid-model form outside the boundary layer. The mixing length scale is locally equal to the grid size. The computational domain is the channel without the linear sub-layer of the boundary layer. The no slip boundary condition (BC) is replaced by a Navier (BC) at the computational wall. Considering the steady state case, we show that the variational finite element model we have introduced, has a solution (vh, ph)h>0that converges to a solution of the steady state Navier-Stokes Equation with Navier BC. MCS Classification : 76D05, 76F65, 65M60, 1 Introduction Numerical simulations of incompressible turbulent flows cannot be performed from the evolutionary Navier-Stokes Equations (NSE), (1.1) ∂tv+ (v· ∇)v−ν∆v+∇p=f, ∇ · v= 0, because of a great computational complexity due to the structure of the turbulence [30]. This is why various mathematical models derived from the NSE are used to simulate some features of turbulent flows, such as their statistical means or their large scales motions, this last way being known as ”Large-Eddy Simulation” (LES), which is our concern in the present paper. LES has attracted much attention these last two decades, especially because of the increasing of computational ressources, enabling to enlarge the range of scales that LES models might simulate. Basically, LES aims at computing filtered fields such as v=G ⋆ v,G being a smooth transfer function [25,30,31,32]. The filtering operation also might be carried out by solving PDE’s [6,7,14,16,23,24]. ∗Dpto. EDAN & IMUS, Universidad de Sevilla, C/Tarfia, s/n. 41012 Sevilla, Spain & BCAM-Basque Center for Applied Mathematics, Avda. Mazarredo, 14. 48009 Bilbao, Spain; [email protected] †Mathematical Research Institute of Rennes, UMR CNRS 6625, European University of Brittany & Fluminance team, INRIA, Bat 22, Campus Beaulieu, 35042 Rennes cedex, France; [email protected] 1 Stresses that appear by filtering the non linear term (v· ∇)vin the NSE, are considered to be diffusive, therefore often modeled by a turbulent diffusion term such as −∇ ·(νt∇v), where νtis an eddy viscosity. One challenge of the LES is the determination of νt. In this paper, we study the case of a channel flow, periodic in the x1−x2axis for simplicity. The first idea is that the projection on finite element spaces, based on a given triangulation {Th}h>0, is a natural filtering operation, so that we seek for vhinstead of v, where vhis the projection of von a suitable finite element space Wh. The second idea is that one can specifically model the eddy viscosity on the boundary layer by means of wall laws. Indeed, following Kolmogorov theory [22], we consider the turbulence to be isotropic at scales small enough inside the flow domain. This assumption yields to take the eddy viscosity of a Kolmogorov-Prandtl-Smagorinsky form, νt=h2|Dvh|1,hbeing the mixing length, which is the standard sub-grid model (SGM) [11]. However, near the wall, turbulence is not isotropic and complexity is higher than far from the wall [29], so that standard SGM cannot be used there. Usual methods proceed as follows. First one uses the known structure of the turbulent boundary layer, as initially described by von K´arm´an [21] and fully developed by Schlichting [34]. Basically, the boundary layer may be split into two sub-layers, the linear sub-layer where the mean velocity profile is linear, and next, the log sub-layer where the mean velocity profile is specified by a log function. Notice that one can consider more sophisticated models to model the boundary layer [35], nevertheless always involving a log law. In all cases, those models involve an essential quantity which is the friction velocity u⋆(see (2.4) in 2.1.ii below). Next, one splits the domain into two subdomains, the boundary layer, and the computational domain which is the domain’s part not containing the boundary layer. One then uses non linear boundary conditions at boundaries of the computational domain such as wall laws [26,27,28]. Based on the fact that today more computational resources are available to increase accuracy for simulating the mean flow inside the log layer, we take as computational domain the domain’s part without the linear sub -layer, using an eddy viscosity of the form νt=hu⋆ inside de log layer, deduced from standard dimensional analysis [11,26]. To conclude the modeling process, it remains to: i) specify how u⋆is calculated, ii) specify boundary conditions (BC) at computational domain boundaries, iii) fix the choice of the mixing length scale. i) We assume that log law holds inside the boundary layer. Thanks to invertibility of the non linear profil, we can define u⋆as u⋆(v,x), that satisfies suitable estimates (see 2.2.iii and estimate (2.17)). ii) As the thickness of the linear sub-layer is very small compared to other scales involved in the problem, a Taylor expansion allows to deduce from the no slip condition at the flow domain boundary a Navier BC at the computational walls ( (3.6) in 3.1.i). This is as if the linear sub-layer would exert a friction over the log sub-layer. iii) The mesh yields natural numerical length scales hK, where hKis a diameter of any K∈ Th. Therefore, one takes νtof the form νt=h2 K|Dvh|on K∈ Thinside the computational domain, and νt=hKu⋆(v,x) on K∈ Thin the log layer (see 3.2.ii). 1Dvh= (1/2)(∇vh+∇vt h) 2 Once this modeling process is completed, we get a model expressed in its variational form over finite element space Wh×Mh, as described in 3.2.iv. So far as we know, this model is totally new, and can be generalized to more complex and realistic geometries thanks to a careful differential geometry analysis, which is a work under progress. We consider all over the paper the steady-state case, which is in coherence with the fact that in a permanent regime and for a developed turbulence, mean fields are steady, which is not in contradiction with the fact that fluctuations might be time dependent. We prove that this variational problem has a solution (vh, ph)∈Wh×Mh(Theorem 4.3) which converges to a solution (v, p) of the steady-state Navier Stokes Equation (NSE) with Navier BC (Theorem 4.4). The paper is organized as follows. We start with general setting. Then we derive from the NSE a description of the boundary layer, introducing the friction velocity. We specify the computational domain and Navier BC, and next we perform the finite element setting and get the model. Finally we state and prove Theorem 4.3 and Theorem 4.4. Acknowledgements: R. Lewandowski thanks Professor Li Tatsien and the ISFMA at Fudan University, Shangha¨ı, China, for their hospitality during spring 2012, where part of this work was performed. He also thanks E. Memin and the Fluminance team at INRIA Rennes for their support during the second term of the academic year 2012-2013. The research work of T. Chac´on was partially supported by the Spanish Government Grant MTM2012-36124-C02-01. 2 General framework 2.1 Channel flow 2.1.i. Geometry, equations and boundary conditions. Let Ωfbe a channel periodic in the x1axis and x2axis, of height 1 + 2din the x3-axis, for a small parameter d << 1, (2.1) Ωf=x= (x1, x2, x3)∈T2×IR3s. t. −d < x3<1 + d, where T2is the two dimensional torus defined by T2=IR2 T2 where T2=2πZZ2 L, and L > 0 is a given length scale. Let Γfdenote (2.2) Γf={x∈T2×IR3s. t. x3=−dor x3= 1 + d}. The steady-state Navier-Stokes equations with the no-slip boundary condition are as follows, (2.3)    (v· ∇)v−ν∆v+∇p=fin Ωf, ∇ · v= 0 in Ωf, v=0on Γf. The source term fis a body force per mass unit, typically the gravity. Assuming f∈ L2(Ωf)3=L2(Ωf), we know that this equation has a solution (v, p)∈W2,3/2(Ω)3× W1,3/2(Ω) (see in [36]), whose norms are bounded by constants that only depend on ν, 3 ||f||0,2,Ωfand d, and pis defined up to a constant. Uniqueness is known when ||f||0,2,Ωf/ν2 is small enough. 2.1.ii. Friction velocity. Let (v, p) be any solution of 2.3. We still denote by vthe trace of von Γf. We deduce from trace Theorems and Sobolev Theorem that v∈W1,3(Γf)3= W1,3(Γf). Therefore, it makes sense to consider Dv·non Γf, where ndenotes the outward-pointing unit normal vector at Γf,Dv= (1/2)(∇v+∇vt). We split the vector Dv·ninto its tangential part and its normal part, (2.4) Dv·n= (Dv·n)τ+ ((Dv·n)·n)n. Let v⋆∈L6(Γf) be defined on Γfby (2.5) v⋆=v⋆(v)(x) = (ν|(Dv·n)τ(x)|)1 2, called the friction velocity associated to vat x∈Γf. 2.2 Boundary layer description 2.2.i. Length scale. Condition of uniqueness to system (2.3) is not satisfied in a steadystate turbulent regime. Let Sbe the set of solutions, which is a closed subset in L2(Ωf). According to [11], one can construct a probability measure µon S. We consider the following velocity friction w⋆∈L6(Γf) defined by, (2.6) w⋆=ZS v⋆(v)dµ(v). We finally define the meanfriction velocity by (2.7) u⋆=1 L||w⋆||0,2,Γf∈IR to which is associated the typical length scale λthat characterises the boundary layer, (2.8) λ=ν u⋆ , assuming u⋆6= 0. We conjecture that u⋆→ ∞ when ||f||0,2,Ωf→ ∞. 2.2.ii. Main assumption about the boundary layer structure. We focus on the bottom of Ωf,{x3=−d}, assuming that the boundary layer at the top {x3= 1 + d}has a similar structure. According to experiments (see in [34]), we assume that in the boundary layer, the mean fluid velocity has a constant direction and only depends on the variable x3, which means v(x) = v(x3)e, for some fixed unit vector e. Without loss of generality, we can assume that e=e1is the unit vector pointing along the x1-axis. Notice that any plane Pof the form P={x3=h}included in the boundary layer, and any vector Northogonal at Pbeing given, our assumption yields in particular v·N= 0 at P. 2.2.iii. Log law. Experiments and suitable assumptions about turbulence [11,34] indicate that the boundary layer can be decomposed into two sub layers: •near the boundary where the velocity profile vis linear (linear sub layer), 4 •the next sub-layer specified by a log profile (log layer). To be more specific, we introduce the dimensionless variable (2.9) z+=x3 λ, and we consider the following continuous function defined on [0, z+ max] by (2.10) L(z+) =    z+if 0 ≤z+≤z+ 0 1 κlog z+ z+ 0+z+ 0if z+ 0≤z+≤z+ max, where κ≈0,41 is the Von K´arm´an constant. In practical calculations, one takes z+ 0≈20, and z+ max ≈100, that measures the thickness of the logarithmic boundary layer, taken to be equal to 100λ. According to experiments [34], boundary layer thickness goes to zero as the Reynolds number goes to infinity. The profile vin the boundary layer at the bottom of Ωfis given by the formula (2.11) v(x3) = u⋆Lx3 λ. A similar description applies to the boundary layer at the top of Ωf,{z= 1 + d}. 2.2.iv. Friction velocity expressed as a function of the velocity. We still focus on the bottom. Any x3>0 being given, Let (2.12) F(β) = βL(αβ), α =x3 ν. With this notation, equation (2.11) may be written as (2.13) v=F(u⋆), thanks to definition (2.8). Lemma 2.1. Let F: [0,+∞)→[0,+∞)be defined by by (2.12). The function Fis invertible, so that equation (2.13) can be written as u⋆=F−1(v)at each given x3. Proof. We observe that the function Lsatisfies (2.14) lim x→0+ L(x) x=C1, (2.15) lim x→∞ L(x) log x=C2, where C1and C2are non-zero constants. As Lis strictly increasing and continuous in (0,+∞), then Fis strictly increasing and continuous in (0,+∞). Also, by (2.14)Fis continuous at β= 0 with F(0) = 0. Moreover, by (2.15) lim x→∞ F(x) = +∞. Then Fis bijective from [0,+∞) onto [0,+∞), which yields the invertibility of Fas claimed.  Lemma 2.2. Denote h=F−1. Then there exist a constant C=C(x3)>0, bounded, such that (2.16) ∀γ > 0, h(γ)≤C(x3) (1 + γ)≤C(1 + γ), C = sup C(x3). 5 Proof. Then h: [0,+∞)7→ [0,+∞) is bijective and continuous. Also, lim γ→∞ h(γ) γ= lim t→∞ t F(t)= lim t→∞ 1 L(αt)= lim t→∞ 1 log(αt) log(αt) L(αt)= 0. The conclusion is a consequence of the continuity of h. We deduce from Lemma 2.1, inequality (2.16) and because top and bottom layers have the same structure, that the friction velocity can be calculated at each x∈BL from the velocity v, and satisfies the estimate (2.17) 0 < u⋆=u⋆(v,x)≤C(1 + |v|). 3 Turbulence model 3.1 Geometry and meshing 3.1.i. Calculation domain. From now we assume that the boundary layer is included in the union of two strips, (3.1) BL ={−d≤x3≤D/2−d} ∪ {1 + d−D/2≤x3≤1 + d}, where d < D << 1, dbeing the order of the linear sub layer, Dthe thikness of the global boundary layer. Standard numerical simulations are carried out in a sub domain of the flow domain that does not include the boundary layer at all, using a wall law [11,27,28] at artificial boundaries (walls). Our model includes the log layer, using a Navier BC based on a Taylor expansion as shown below. The computational domain is (3.2) Ω = x= (x1, x2, x3)∈T2×IR3s. t. 0 < x3<1, the artificial wall being defined by (3.3) Γw={x∈T2×IR3s. t. x3= 0 or x3= 1}. 3.1.ii. Boundary conditions. As in above, we focus on the bottom layer. By a Taylor expension we get (3.4) 0 = v|x3=−d≈v|x3=0 −d∂v ∂x3 |x3=0. From the view point of the domain Ω, v=vτ|Γw, and ∂/∂x3=−∂/∂nat Γw, where vτis the tangential part of v, defined by (3.5) v=vτ+ (v·n)n, by still denoting vthe trace of vat Γw, so far no risk of confusion occurs. Therefore, by remarks in 2.2.ii together with (3.4), we get (3.6) v·n|Γw= 0,∂vτ ∂n|Γw=−1 dvτ, 6 which is a Navier boundary condition at the artificial wall, that expresses in some sense that the linear sub-layer exerts a friction on the log layer. Hence, system (2.3) becomes in Ω, (3.7)          (v· ∇)v−ν∆v+∇p=fin Ω, ∇ · v= 0 in Ω, v·n=0on Γw, −∂vτ ∂n=1 dvτon Γw. Navier-Stokes equations with Navier boundary conditions was studied before [1,2,3,10, 37], and existence of a solution to system (3.7) is already ensured. 3.1.iii. Variational formulation. Let us define the spaces W(Ω) = {w∈H1(Ω),w·n|Γw= 0}, M(Ω) = {q∈L2(Ω),ZΩc q dx= 0}, by reminding that H1(Ω) = H1(Ω)3. Strictly speaking, The space M(Ω) is isomorphic to the quotient space L2(Ω)/IR, endowed with the usual quotient norm (3.8) || . p||M= inf p∈. p||p||0,2,Ω. It also may be viewed as a closed subspace of L2(Ω) endowed with the L2(Ω) norm. The space W(Ω) is endowed with the H1norm, denoted || · ||1,2,Ω. As a consequence of Korn’s inequality, the following usefull estimate holds, (3.9) ∀v∈W(Ω),||v||1,2,Ω≤C(||Dv||0,2,Ω+||v||0,2,Γw), of proof of which being carried out in [10]. Let a,band Gthe forms defined by a(v,w) = ν(Dv, Dw)Ω,(3.10) b(z;v,w) = 1 2[((z· ∇)v,w)Ω−((z· ∇)w,v)Ω](3.11) G(v,w) = ν d(vτ,wτ)Γw,(3.12) for z,v,w∈H1(Ω). Recall that when z,v,w∈W(Ω) and ∇ · z= 0, then b(z;v,w) = ((z·∇)v,w)Ω, and (∇z,∇w)Ω= (Dz, Dw)Ω. Also remark that when v∈b, then v=vτ at Γw. We say that a pair (v, p)∈W(Ω) ×M(Ω) is a weak solution of the boundary value (3.7) if it satisfies (3.13) b(v;v,w) + a(v,w)−(p, ∇ · w)Ω+G(v,w) = hf,wi, (∇ · v, q)Ω= 0, for any (w, q)∈W(Ω) ×M(Ω). 3.1.iv. A priori estimate and existence result. Assume f∈W(Ω)′. Let (v, p) be any solution of Problem (3.13), and take v=win (3.13). From the standard formula b(v;v,v) = 0 that holds since ∇ · v= 0 and v·n= 0 at Γw, we get (3.14) ν||Dv||0,2,Ω+ν d||v||0,2,Γw=hf,vi, 7 from where we deduce (3.15) ||v||1,2,Ω≤Cκ−1||f||W(Ω)′, κ = min ν, ν d, by using (3.9). 3.2 Finite element setting 3.2.i. Triangulation. Let D⊂IR3denotes the sample box D= [0, L]2×[0,1]. The computational domain Ω may be viewed as the periodic reproduction of Din the x1− x2axes. Let {Th}(h>0) be a regular familly of triangulation of D, compatible with the periodicity of the domain: The restriction of the grid to the planes x1= 0 and x1=D is the same, and the restriction of the grid to the planes x2= 0 and x2=Dis the same. Reproducing this triangulation by periodicity, we get a regular triangulation of Ω, still denoted by {Th}(h>0). In the following, for each K∈ Th,hK=diam(K) denotes the diameter of K, and (3.16) h= max K∈Th hK. 3.2.ii. Eddy viscosities. We assume isotropy of the turbulence inside the domain defined by Ωin = Ω\BL, the boundary layer BL being defined by (3.1). This yields to consider in Ωin the eddy viscosity νt,in to be of the same form as in usual Sub-Grid-Models of Prandtl- Kolmogorov-Smagorinsky type, where following [11], we take in each Kthe length scale equal to hK, leading to consider νt,in to be of the form (3.17) νt,in(v) = C2 sX K∈Th h2 K1K|Dv|, Cs>0 being an empirical constant, 1Adenotes the characteristic function for any set A. In the boundary layer part, Ωw=BL ∩Ω, turbulence is no longer isotropic and depends on the friction velocity. Taking again hKas typical length scale and by a dimensional analysis argument [11], we define the eddy viscosity νt,w in Ωwby (3.18) νt,w(v) = CwX K∈Th hK1Ku⋆(v,x), where Cw>0 is an empirical constant and u⋆is expressed in 2.2.iv. Finally, the eddy viscosity we consider is of the form (3.19) νt=νt(v) = 1Ωin νt,in(v) + 1Ωwνt,w(v), 3.2.iii. Finite element spaces. The model is a mixed formulations, based upon pairs of finite element spaces (Wh, Mh)⊂W(Ω) ×M(Ω), associated to the family of regular triangulations {Th}h>0of Ω in the sense of Ciarlet [12]. We assume that the family of pairs of spaces {(Wh, Mh)}h>0satisfies the following hypothesis: Hypothesis 1: The family of spaces {Wh×Mh}h>0is an internal approximation of W(Ω) ×M(Ω): For all (w, p)∈W(Ω) ×M(Ω) there exists a sequence {(vh, ph)}h>0such that (vh, ph)∈Wh×Mh, and lim h→0(kv−vhk1,2,Ω+kp−phk0,2,Ω) = 0. 8 Hypothesis 2: The family of pairs of spaces {(Wh, Mh)}h>0satisfies the uniform discrete inf-sup condition : There exists a constant α > 0 such that (3.20) αkqhk0,2,Ω≤sup wh∈Wh (∇ · wh, qh)Ω kwhk1,2,Ω ,for all qh∈Mh There is a wide literature about finite element spaces satisfying those properties (Cf. [5], [9], [15]). 3.2.iv. The model. Our LES model is expressed by the following variational problem: Find (vh, ph)∈Wh×Mhsuch that for all (wh, qh)∈Wh×Mh, (3.21)    b(vh;vh,wh) + a(vh,wh) + c(vh;wh) + G(vh,wh)−(ph,∇ · wh)Ω=hf,whi, (∇ · vh, qh)Ω= 0; the form cbeing defined by (3.22) c(v;w) = (νt(v)Dv, Dw)Ωor c(v;w) = (νt,in(v)Dv, Dw)Ωin + (νt,w(v)∂3v, ∂3w)Ωw. The second expression neglects the tangential eddy viscosity in the boundary layer, which is very small compared to the normal one. 4 Analysis of the model 4.1 Technical results We state in this sub-section some technical results concerning the eddy viscosities and the associated turbulent diffusion form c, that are needed by our analysis. 4.1.i. L∞eddy viscosties estimates. Lemma 4.1. There exists a constant C > 0depending only on the aspect ratio of the family of triangulations {Th}h>0such that (4.1) kνt(vh)k0,∞,Ω≤C h1/2kvhk1,2,Ω,for all vh∈Wh. Proof. We start with the internal part of the eddy viscosity νt,in. Consider vh∈Wh. As ∇vhis piecewise continuous, there exists K∈ Thsuch that kνt,in(vh)k0,∞,Ω=kνt,in(vh)k0,∞,K ≤C2 Sh2 Kk∇vhk0,∞,K. By a standard finite element inverse estimate (Cf. [4]), k∇vhk0,∞,K ≤C h−3/2 Kk∇vhk0,2,K for some constant C > 0 depending only on the aspect ratio of the family of triangulations. Then, (4.2) kνt,in(vh)k0,∞,Ω≤CC2 Sh2−3/2 Kk∇vhk0,2,K ≤CC2 Sh1/2k∇vhk0,2,Ω, 9 [24] Layton, W., Lewandowski, R.: A simple and stable scale similarity model for large eddy simulation: energy balance and existence of weak solutions. Applied Maths. Letters, 16, 1205-1209, (2003) [25] Lesieur, M., M´etais, O., Comte, P. Large-eddy simulations of turbulence. Cambridge University Press, (2005) [26] Lewandowski, R.: Analyse Math´ematique et Oc´eanographie. Masson (1997) [27] Mohammadi, M., Pironneau, O.: Analysis of the k-epsilon turbulence model. Masson, Paris, (1994) [28] Par´es, C.: Existence, uniqueness and regularity of solution of the equations of a turbulence model for incompressible fluids. Applicable Analysis 43, 245-296 (1992) [29] Piomelli, U., Balaras, E.: Wall-layer models for large-eddy simulations. Annu. Fluid Mech. 34, 349-374, (2002) [30] Pope, S.B.: Turbulent flows. Cambridge University Press, (2000) [31] Rodi, W., Ferziger, J.H., Breuer, M., Pourqui´e, M.: Status of large-eddy simulation. Workshop on LES of flows past bluff bodies (Rottach-Egern, Tegernsee, 1995). ASME J. Fluid Engng., 119, 248–262, (1997) [32] Sagaut, P.: Large Eddy Simulation for Incompressible Flows. Springer, Belin (2002) [33] Scott, R., Zhang, S.: Finite element interpolation of non-smooth functions satisfying boundary conditions. Math. Comput. 54, 483–493 (1990) [34] Schlichting, H. Boundary Layer Theory, 8th. Edition. Springer (2000) [35] Spalding, D.B.: A single formula for the law of the wall. Transactions of the ASME, Series E: Journal of Applied Mechanics, 28, 455–458 (1961) [36] Temam, R.: Navier-Stokes equations. AMS Chelsea (2001) [37] Verf¨urth, R.: Finite Element approximation of incompressible Navier-Stokes equations with slip boundary conditions. Numer. Math. 50, 697–721 (1987) 16