Homogenization of a large eddy simulation model for turbulent fluid motion near a rough wall

Homogenization of a Large Eddy Simulation Model for Turbulent Fluid Motion near a Rough Wall Matthieu Bonnivard∗Francisco J. Su´arez-Grau† September 6, 2017 Abstract In this paper, we are interested in the numerical simulation of a turbulent viscous flow in the vicinity of a periodic rough wall. To this aim, we consider a Large Eddy Simulation (LES) model, and perform its asymptotic analysis, letting the period εof the roughness tend to zero. Relying on an homogenization method closely related to the two-scale convergence method [2, 38], we exhibit a critical scaling associated with parameter λ:= limε→0aε/ε5/3, where aεstands for the amplitude of the roughness. In the critical regime 0 < λ < ∞, the roughness effect is described by an additional, nonlinear friction term, that depends locally on the solution to a nonlinear auxiliary problem of boundary layer type. In the particular case of riblets, which is of major practical interest, we show that the directional invariance of the geometry results in a simplified expression of this extra term, and we propose a numerical method to simulate the homogenized LES model, based on that observation. Finally, we study the influence of parameter λon the flow approximated by the homogenized model, using different commonly used riblet shapes. Our simulations demonstrate that the extra friction term can be handled numerically, and may have an important influence on the results of LES in presence of roughness. AMS Classification: 35B27, 76F65. Keywords: turbulence model; large-eddy simulation; roughness effect; homogenization. 1 Introduction The study of roughness effects in fluid-structure interaction has seen considerable development in the last decades. Its main objective is to explain how the fine structure of solid walls may influence viscous flows in their vicinity, and how the effects of surface irregularities should be taken into account in order to enhance the precision of the existing models. One of the first mathematical results in this direction concerns the justification of the no-slip boundary condition, which is commonly used in continuum mechanics models and prescribes a null velocity for the fluid particles along the solid boundaries. In their pioneering paper [16], Casado-D´ıaz et al. have explained that boundary condition using an asymptotic analysis based on the presence of small-scale irregularities on the wall, and on a minimal boundary condition – the non ∗Univ. Paris Diderot, Sorbonne Paris Cit´e, Laboratoire Jacques-Louis Lions, UMR 7598, UPMC, CNRS, F-75205 Paris, France †Dpto. de Ecuaciones Diferenciales y An´alisis Num´erico, Fac. de Matem´aticas, Univ. de Sevilla 1 penetration condition – expressing the impossibility for the fluid to get through the solid. In a simplified manner, that phenomenon could be interpreted as follows: the rough geometrical structure of the wall causes the normal nto the boundary to oscillate with high spatial frequency, turning non penetration condition u·n= 0 into no-slip condition u= 0, where ustands for the fluid velocity. The study of such rugosity effects is of particular importance in the field of microfluidics. Indeed, when fluid-structure interaction occurs at micro or nanoscale, surface effects become prominent due to miniaturization, and may have a strong impact on the behaviour of the fluid. At those scales, positive slip lengths have been measured experimentally, which justify the introduction of Navier boundary conditions in the Stokes models used to describe that type of configuration. We refer to [32] for a review of such experimental results. In the framework of fluid-solid contact with slip, numerous mathematical results have been obtained, that draw a rather complete picture of the effects that rough walls may produce on viscous flows at low Reynolds number. The stronger effect consists in forcing the fluid velocity to vanish, in the limit where the characteristic dimensions of the roughness go to zero, as was previously described. The rugosity effect can also be directional, when the surface texture forces the fluid to flow in a specific direction, tangential to the wall. Last, rough walls may produce an extra friction coefficient (or in other words, reduce the slip length), potentially anisotropic. We refer to [14] for a general mathematical description of those effects, that can also be combined. By specifying particular types of rough patches, it is possible to quantify precisely such effects, especially by computing friction coefficients. One typical example is the case of periodic distributions of roughness, that can be treated by computating a friction matrix, whose coefficients are determined by solving boundary layer problems [17, 28]. Among such periodical roughness profiles, the riblets, which are assemblies of the same elongated structure, are of considerable interest. They are widely used in engineering, extensively studied in experiments and mechanical models (see [7,8,25,30,33]) and have been the subject to many mathematical studies (see, for instance, [3,13, 29]). However, in many applications such as aeronautics [45] or shark-skin drag reduction [20], rugosity effects (and especially, riblet effects) occur in turbulent regime. At present, a rigourous mathematical analysis of fluid-structure interaction at high Reynolds number remains out of reach, due to the well-known, tremendous difficulties posed by the three-dimensional Navier-Stokes equations in that context. From a numerical point of view, the direct numerical simulation of those equations, at high Reynolds number, can be very expensive, since it requires solving a very wide range of time and length scales – whose presence is characteristic of turbulence. As a result, lots of tools and models have been developed to describe and simulate turbulent flows, that are commonly used in mechanics [39]. In this paper, we focus on the class of Large Eddy Simulation (LES) models, that were introduced by Smagorinsky in 1963 [42] and are widely used in computational fluid dynamics. The principle of LES is to apply a spatial convolution operator to the Navier-Stokes equations, in order to decompose the velocity, pressure and body forces into two parts: an averaged part describing the global, large scale behaviour of the fluid, and a remainder corresponding to the small scales. The goal of LES models is to compute only the large scale structure of the flow. However, due to the presence of the nonlinear convective term in the moment equation, there is an interaction between small and large scales that one needs to model in order to close the average Navier Stokes system of equations. This well-known closure problem has led to a large variety of LES models, whose derivation can be motivated by physical or by mathematical considerations. The mathematical approach is based on the use of the Fourier transform, as will be detailed in Section 2. We refer to [46] and the references therein for a description of the main LES models based on this approach. In this work, we address the problem of simulating numerically the turbulent flow of a viscous fluid, in the vicinity of a rough wall, using one of the simplest and most popular LES model: the Smagorinsky 2 model. As was proved by Coletti in [19], this model enjoys nice mathematical properties, in particular, global existence and uniqueness of a weak solution. In order to develop an asymptotic analysis of the model, with respect to the characteristic length of the rough pattern, we rely on a remark by Galdi and Layton [24], which is also developed in [46] (Chap. 5). Starting from the observation that the noslip boundary condition, which is the commonly accepted boundary condition for viscous flows at high Reynolds number, cannot be preserved by spatially filtering the velocity using convolution, the previous authors have proposed to equip LES models with Navier boundary conditions. As stated in the beginning of this Introduction, such boundary conditions are well adapted to the study of rugosity effects. As a matter of fact, we are able to generalize the multi-scale analysis developed in [17], to a time-dependant system of equations, involving a nonlinear stress tensor, in the case of a periodic rough wall. We identify a new critical scaling between the amplitude and the period of the oscillations, and obtain an effective system set in the homogenized domain, that expresses the effect of the rough patches on the LES model, by adding a nonlinear interaction term to the initial boundary condition of Navier type. In this effective system, the homogenized boundary condition is determined locally by solving a boundary layer problem that involves only the nonlinear part of the stress tensor – that is, the part modelling the large scales/small scales interactions. Since our work is motivated by the numerical simulation of turbulent flow near a rough wall, we conclude by addressing a particular geometric configuration, which is widely studied in the experimental fluid mechanics literature, and turns out to be well-adapted for simulating the limit system: the case of riblets. Indeed, using their directional invariance, we are able to obtain a simplified expression of the extra friction term, that requires the solving of only two different auxiliary problems. Then, we propose a numerical method to solve the limit LES model, based on a fixed point algorithm for treating the nonlinear part of the stress tensor, and the nonlinear friction term accounting for the rugosity effect. We implement that method to simulate the steady-state of the LES model, using three different riblet shapes that are commonly used in the literature, and discuss the influence of parameter λfor these geometries. Our simulations demonstrate the extra friction term emerging from the asymptotic analysis of the model, can be taken into account in the applications, and have a non neglictable impact on the physical properties of the flow approximated by Large Eddy Simulation. The paper is organized as follows. The Smagorinsky model is presented in Sec. 2. The main results are stated in Sec. 3, and proved in Sec. 4. Finally, Sec. 5 is dedicated to numerics. 2 The Smagorinsky model for turbulent flow We consider the flow of a viscous and homogeneous fluid, driven by an external force, in a threedimensional domain delimited by two fixed horizontal plates. We assume that the upper plate is perfectly smooth, but that the lower one presents a periodic rough surface. We work in nondimensionalized coordinates and denote by Re the Reynolds number associated with the flow. Let us start with the geometric notation that will be used throughout the paper. 2.1 Geometric notation We denote by (e1, e2, e3) the canonical basis of R3and by x= (x1, x2, x3) the spatial coordinates. The top smooth boundary occupies the plane x3= 1, with normal e3. For the sake of simplicity, we assume that the flow is periodic with respect to the horizontal coordinates (x1, x2), which motivates the introduction of the tori T2:= R2/Z2and T1:= R/Z. To describe the geometry of the lower (rough) plate, we fix a positive function Ψ ∈W2,∞(T2) and define Ψε(x1, x2) = aεΨx1 ε,x2 ε,a.e. (x1, x2)∈ T 2, 3 where ε > 0 is a small parameter and aεsatisfies lim ε→0aε= 0,lim ε→0 aε ε= 0.(2.1) The rough boundary Γεis then defined by Γε:= {(x1, x2, x3)∈ T 2×R:x3=−Ψε(x1, x2)}.(2.2) By the periodicity of Ψ, Ψεis periodic of period εin the x1and x2directions, and by construction, the maximum value of Ψεis proportional to aε. Thus, hypotheses (2.1) mean that the amplitude of the roughness is negligible with respect to the period, as the period tends to zero. Let νthe be the outside unitary normal vector to Ωεon Γε. The orthogonal projection on the tangent space of ∂Ωεwill be denoted by T, i.e. Tξ =ξ−(ξ·ν)ν, ∀ξ∈R3,a.e. on Γε. We denote by Γ0the upper (smooth) boundary, defined by Γ0:= T2× {1}, and introduce the fluid domain Ωεdefined by Ωε:= {(x1, x2, x3)∈ T 2×R:−Ψε(x1, x2)< x3<1}. Since the sequence (Ψε)ε>0converges uniformly to zero, the sequence (Γε)ε>0converges in the sense of Hausdorff to the flat wall Γ defined by Γ = T2× {0}.(2.3) Accordingly, we introduce the limit domain Ω = T2×(0,1),(2.4) and define Λ by Λ = T2×(−1,1), so that Λ contains all the rough domains Ωε. 2.2 Derivation of the LES model Let uεbe the nondimensionalized velocity of the fluid and pεthe nondimensionalized pressure. We assume that the external force is given by the restriction to Ωεof a vector field f∈L2(Λ,R3). The fluid motion is described by the Navier-Stokes equations ∂tuε+ div (uε⊗uε)−1 Re∆uε+∇pε=fin ]0, T[×Ωε,(2.5) div uε= 0 in ]0, T[×Ωε,(2.6) completed with the initial condition uε(x, 0) = u0 ε(x), x ∈Ωε(2.7) 4 and no-slip boundary conditions on the plates: uε= 0 on Γ0,(2.8) uε= 0 on Γε.(2.9) As it is well known, in the three-dimensional case, existence and uniqueness of the solution of (2.5)– (2.7) has only been proved for small Reynolds number [31]. For large values of Re, analytical solutions of Navier-Stokes equations become unstable; this is called turbulence. Our aim is to investigate numerically the influence of the rough boundary on turbulent states of the system. Although several numerical procedures have been developed to solve equations (2.5)-(2.6), they suffer the same instability issues as the analytical solutions, when the Reynolds number becomes large. To enhance the stability of the numerical solutions, one needs to consider a very fine mesh of the space domain, which increases the storage and computational costs dramatically, especially when the boundary of the domain presents fast oscillations. As an alternative to direct numerical simulation (DNS), several models of turbulent flows have been developed. We refer to [37] for an introduction to the main models of turbulence. In order to simulate flows at high Reynolds number, we will apply the space filtered approach to large eddy simulation (LES), introduced in 1970 by Deardorff [21]. Let us recall the main steps of the derivation of this model; for more details on the derivation of classical LES models, we refer to [9,19,40]. The spatial averaging operator is defined as a convolution with the Gaussian filter gδ(x) := γ π3/21 δ3exp[−γ|x|2/δ2], where δ > 0 represents an averaging radius and γ > 0 is a constant. All dependent variables being extended by zero to T2×R, the average of any variable r(t, x) is defined by r(t, x)=(gδ∗r)(t, x) = ZR3 gδ(x−y)r(t, y)dy. Applying this averaging operator to system (2.5)-(2.7) and commuting spatial partial derivatives with convolutions, we formally obtain the system ∂tuε+ div (uε⊗uε)−1 Re∆uε+∇pε=fin ]0, T[×Ωε,(2.10) div uε= 0 in ]0, T[×Ωε,(2.11) uε(x, 0) = u0 ε(x), x ∈Ωε.(2.12) However, the nonlinear term div (uε⊗uε) cannot be expressed as a function of the new unknown uε. This notorious closure problem prevents us from solving the averaged problem (2.10)-(2.12) directly. To overcome this difficulty, it is required to derive an approximate expression of uε⊗uε, involving only the average velocity field uε. Introducing the small-scale velocity u0 εdefined by uε=uε+u0 ε, one starts with the following development: uε⊗uε=uε⊗uε+uε⊗u0 ε+u0 ε⊗uε+u0 ε⊗u0 ε.(2.13) The principle of the approximation is to apply the Fourier transform to every term on the right hand side of (2.13), to exploit its commutation properties with respect to convolution, and to use a formal 5 asymptotic development of the the result with respect to the small-scale δ, as δtends to zero, before returning to space variables by applying the inverse Fourier transform. Denoting by Fthe Fourier transform on T2×Rand by ξthe dual variable, the Fourier transform of the Gaussian filter is defined by bgδ(ξ) := [F(gδ)] (ξ) = exp −δ2 4γ|ξ|2, and using the properties of the convolution, there holds F(uε⊗uε) = bgδ(ξ)F(uε⊗uε) = e−δ2 4γ|ξ|2(b uε∗b uε). Using a Taylor series approximation of the exponential, we can write the asymptotic development F(uε⊗uε) = 1−δ2 4γ|ξ|2(b uε∗b uε) + O(δ4),as δ→0. Applying the same treatment to the three other terms in the right-hand side of (2.13), taking the inverse Fourier transform and neglecting every term of order O(δ4) or higher, one obtains the approximation uε⊗uε≈uε⊗uε+δ2 2γDuε(Duε)T, where Duεstands for the Jacobian matrix of uε. This leads to a classical space-filtered LES model that is used in many studies, see for instance [1,5,6]. However, following the previous approach, the turbulent fluctuations u0 ε⊗u0 εin (2.13) turn out to be formally of order O(δ4), and thus neglictable. In the present work, we choose to model u0 ε⊗u0 εby a Smagorinsky term of the form −Cs|D(uε)|D(uε), where Cs>0 is of order O(δ2). This term was introduced in [42]. As a result, the model that we consider relies on the following approximation: uε⊗uε≈uε⊗uε+δ2 2γDuε(Duε)T−Cs|Duε|Duε. Defining (wε, qε) as an approximation of (uε, pε), we obtain the system ∂twε+ (wε· ∇)wε−div S(Dwε) + ∇qε= 0 in ]0, T [×Ωε,(2.14) div wε= 0 in ]0, T[×Ωε,(2.15) wε(x, 0) = w0 ε(x), x ∈Ωε,(2.16) where the stress tensor S(Dwε) is defined by S(Dwε) = 1 Re +Cs|Dwε|Dwε−δ2 2γDwε(Dwε)T.(2.17) In (2.16), the initial datum w0 εis obtained by setting w0 ε(x) := u0 ε(x) a.e. in Ωε. Treatment of the boundary conditions. To complete system (2.14)–(2.16), we need to impose boundary conditions. Formally, the approximation wεshould satisfy wε=uεon the boundary of the domain, so we need to replace such boundary conditions by suitable approximations in order to close the system. Since we want to focus on the interaction between the fluid and the rough boundary Γε, for simplicity, we impose no-slip boundary condition on the top boundary wε= 0 on Γ0.(2.18) 6 To derive a boundary condition for wεon the rough plate Γε, we follow a strategy proposed by Galdi and Layton [24]. In their paper, these authors suggest that the averaged velocity field, which is a model for the motion of the so-called large eddies, should satisfy a Navier slip boundary condition on solid boundaries, which can be seen as the most general boundary condition for a fluid. Some numerical evidence of this fact are cited by Serrin in [41], for high Reynolds flows of gases. In the context of the nonlinear stress tensor (2.17), such boundary conditions read wε·ν= 0 on Γε, [S(Dwε)ν]·τi+β wε·τi= 0 on Γεi= 1,2. (2.19) In these expressions, νdenotes the outward unit vector on Γε,{τ1, τ2}is an orthonormal system of tangent vectors on Γεand β=β(Re, δ) is a nonnegative friction coefficient, that depends typically on the Reynolds number Re and on the averaging radius δ. In particular, for consistency, the model should impose that βgoes to infinity as δgoes to 0, so that formally, boundary conditions (2.19) turn into no slip conditions as the averaging radius vanishes, conformly to the physical boundary conditions (2.9) satisfied by the fluid. In the present study, we make the assumption that, δand Re being fixed, βremains of order 1 with respect to the roughness parameter ε, as εtends to zero. Thus, we assume for simplicity that β(Re, δ) is independant on ε. As it will be clear in the proofs, as long as the friction coefficient remains bounded when εgoes to zero, it does not play any significant role in the asymptotic analysis of system (2.14)–(2.19). Remark 2.1. In the computations, we will often divide the tensor S(Dwε)in two parts: the Newtonian and non-Newtonian parts. Namely, we set S(Dwε) = 1 ReDwε | {z } Newtonian part +b S(Dwε) | {z } non-Newtonian part , with b S(Dwε) = Cs|Dwε|Dwε−δ2 2γDwε(Dwε)T. Remark 2.2. Tensor S(ξ)satisfies the homogeneity condition S(t ξ) = t Reξ+t2b S(ξ),∀ξ∈R3×3,∀t∈R+.(2.20) Weak formulation of system (2.14)–(2.19).We introduce the functional spaces Hε=L2(Ωε,R3), Wε=w∈W1,3(Ωε,R3), w ·ν= 0 on Γε, w = 0 on Γ0,(2.21) Vε={w∈Wε,div w= 0 in Ωε}. Wεand Vεare closed subspaces of W1,3(Ωε,R3), and thus Banach spaces for the same norm kDwkL3(Ωε,R9), which is equivalent to the usual W1,3norm by Poincar´e inequality. We define Jε=H1(]0, T[, Hε)∩L3(]0, T[, Vε).(2.22) 7 Assume that w0 ε∈Vεis given. We say that wεis a weak solution to problem (2.14)–(2.19) if wε∈Jεand satisfies initial condition (2.16), and if the integral identity ZT 0ZΩε ∂twε·ϕ dx dt +ZT 0ZΩε [(wε· ∇)wε]·ϕ dx dt +ZT 0ZΩε S(Dwε) : Dϕ dx dt +βZT 0ZΓε wε·ϕ dH2dt =ZT 0ZΩε f·ϕ dx dt (2.23) holds for every test function ϕsuch that ϕ∈C1 c(]0, T[×Ωε,R3),div ϕ= 0 in Ωε, ϕ ·ν= 0 on Γε, ϕ = 0 on Γ0. In the surfacic integral on Γε,H2stands for the bidimensional Hausdorff measure in R3. Notice that, by Sobolev embedding theorem, the space Jεis continuously embedded in C0([0, T], Hε), so that initial condition (2.16) is well defined in Hε. Our main objective is to study the asymptotic behaviour of weak solutions to system (2.14)–(2.19) when εtends to zero, and to identify the homogenized model that accounts for the roughness induced effects. 3 Main results The first result ensures the existence and uniqueness of a weak solution wε∈Jεto problem (2.23), and states the expected energy estimates. It is adapted from Coletti [19]. Theorem 3.1. Let T > 0and f∈H1(]0, T[, L2(Λ,R3)). Assume that Cs>δ2 γ. Then for every ε > 0and every initial datum w0 ε∈Vε, there exists a unique weak solution wε∈Jεto problem (2.14)– (2.19). Moreover, there exists a constant C=C(Re, Cs, δ, γ)such that this solution satisfies the energy inequalities: sup 0≤t≤TZΩε |wε|2dx≤ZΩε |w0 ε|2dx +CZT 0ZΩε |f|2dx dt, (3.1) 1 Re ZT 0ZΩε |Dwε|2dx dt + 2 ZT 0ZΓε β|wε|2dH2dt + 2 Cs−δ2 2γZT 0ZΩε |Dwε|3dx dt (3.2) ≤ZΩε |w0 ε|2dx +CZT 0ZΩε |f|2dx dt, ZT 0ZΓε |∂twε|2dσ dt ≤CZΩε |w0 ε|2dx +ZΩε |Dw0 ε|2dx +ZΩε |Dw0 ε|3dx +ZT 0ZΩε |f|2dx dt +ZT 0ZΩε |∂tf|2dx dt.(3.3) In the next Proposition, we construct a pressure field qε∈W−1,∞(]0, T[, L3/2 0(Ωε)), associated with the divergence free constraint on wε. 8 Proposition 3.2. Let T > 0,f∈H1(]0, T[, L2(Λ,R3)) and assume that Cs>δ2 γ. Let ε > 0,w0 ε∈Vε and let wε∈Jεbe the weak solution to problem (2.14)–(2.19). Then, there exists a distribution qε∈ W−1,∞(]0, T[, L3/2 0(Ωε)) such that for every ϕ∈C1 c(]0, T[, Wε), ZT 0ZΩε ∂twε·ϕ dx dt +ZT 0ZΩε [(wε· ∇)wε]·ϕ dx dt +ZT 0ZΩε S(Dwε) : Dϕ dx dt +βZT 0ZΓε wε·ϕ dH2dt −ZT 0ZΩε qεdiv ϕ dx dt =ZT 0ZΩε f·ϕ dx dt . (3.4) Our main contribution is the obtention of the asymptotic behaviour of the velocity/pressure couple (wε, qε) satisfying (3.4), as εtends to zero. The statement of the effective system requires some extra notation to define the auxiliary (cell) problem, that we introduce in the next paragraph. Notation related to the cell problem. By Y0, we denote the unitary cube of R2,Y0= (−1 2,1 2)2, and by b Qthe set b Q=Y0×(0,+∞). For every M > 0 we write b QM=Y0×(0, M). We use the index # to mean periodicity with respect Y0, for example Lp #(Y0) denotes the space of functions u∈Lp loc(R2) which are Y0-periodic, while Lp #(b Q) denotes the space of functions bu∈Lp loc(R2×(0,+∞)) such that Zb Q |bu|pdy < +∞,bu(y0+k0, y3) = bu(y),∀k0∈Z2,a.e. y∈R2×(0,+∞). We denote by Vthe space of functions bv:R2×(0,+∞)→Rsuch that ∀M > 0bv∈W1,3 #(b QM,R),and Dbv∈L3 #(b Q, R3). It is a Banach space endowed with the norm k·kVdefined by kbvk3 V=kbvk3 L3 #(Y0×{0},R)+kDbvk3 L3 #( b Q,R3). Asymptotic behaviour of (wε, qε).Depending on the scaling of the amplitude aεwith respect to ε, the asymptotic behaviour of system (2.14)–(2.19) is described by three different regimes, that we call critical, sub-critical and super-critical. The relevant scaling is controlled by parameter λ, which is a nonnegative number, possibly infinite, defined by λ:= lim ε→0 aε ε5 3 .(3.5) Here, we state the result in the critical case; the other cases will be discussed in a Remark 3.7. Theorem 3.3. Let T > 0,f∈H1(]0, T[, L2(Λ,R3)) and assume that Cs>δ2 γ. We consider a family of initial data w0 ε∈Vε, satisfying the uniform bound ∀ε > 0kw0 εkWε≤C, where C > 0is a constant. Without loss of generality, we may assume the existence of w0∈W1,3(Ω,R3) such that w0 ε|Ωconverges weakly in W1,3(Ω,R3), to w0. For every ε > 0, let wε∈Jεbe the weak solution to problem (2.14)–(2.19) and qε∈W−1,∞(]0, T[, L3/2 0(Ωε)) be the associated pressure. Assume aε≈ε5 3with aε/ε5 3→λ,0<λ<+∞. Then, there exists (w, q)∈L3(]0, T [, W 1,3(Ω,R3)) ×W−1,∞(]0, T [, L 3 2 0(Ω)) such that wε|Ω* w in L3(]0, T[, W1,3(Ω,R3)), qε|Ω* q in W−1,∞(]0, T[, L3 2(Ω)),(3.6) 9 Remark that κis well defined up to a set of measure zero in R2(the set ∪k0∈Z2∂Ck0 1). Moreover, for every ρ > 0 and k0∈Z2, we have κx0 ρ=k0⇔x0∈Ck0 ρ. Finally, for a.e. x0∈R2we define Cε(x0) as the square Ck0 εsuch that x0belongs to Ck0 ε. The first step in this approach is to introduce a microscopic variable y, as follows. Given (wε, qε)∈ L3(]0, T[, W1,3(Ωε,R3)) ×L3 2(]0, T[, L 3 2 0(Ωε,R)), with wε= 0 on Γ0, a.e. in ]0, T[, we define ( bwε,bqε) by bwε(t, x0, y) = wεt, εκ x0 ε+εy0, εy3(4.14) bqε(t, x0, y) = qεt, εκ x0 ε+εy0, εy3,(4.15) for a.e. (t, x0, y)∈]0, T[×T 2×b Yε, with b Yε=y∈Y0×R:−aε εΨ(y0)< y3<1 ε. Remark 4.3. For k0∈Z2the restriction of (bwε,bqε)to ]0, T [×Ck0 ε×b Yεdoes not depend on x0, while as a function of yit is obtained from (wε, qε)by using the change of variables y0=x0−εk0 ε, y3=x3 ε,(4.16) which transforms Qk0 εinto b Yε. In order to capture the asymptotic behaviour of (wε, qε) near ]0, T[×Γε, let us first study, in Lemmas 4.4 and 4.5, the asymptotic behaviour of ( bwε,bqε) in ]0, T [×T 2×b QM, for every M > 0. These results represent a generalization to the evolutive case of Lemmas 4.5 and 5.6 in [43], originally stated for the stationary case in the particular case p= 3. Lemma 4.4. Let πεbe a bounded sequence in L∞(]0, T [, L3/2(Ωε,R)). For a.e. (t, x0, y)∈]0, T[×T 2×b Yε we define bπε(t, x0, y) = πεt, εκ x0 ε+εy0, εy3. Then there exists bπ∈L∞(]0, T [, L3/2(T2×b Q, R)) such that ∀M > 0ε2/3bπε*bπweakly- ∗in L∞(]0, T [, L3/2(T2×b QM,R)).(4.17) Moreover, defining qε=−∂tπεand bq=−∂tbπ, there holds for every bϕ∈C1 c(]0, T[×T 2, C1 #(b Q, R3)): lim ε→0"ε−1/3Z(0,T)×Ω qε(t, x) (divybϕ)(t, x0,x ε)dx dt#=Z(0,T)×T 2× b QMbqdivybϕ dx0dy dt. (4.18) 16 Proof. Let M > 0 and ε > 0. We write for a.e. t∈]0, T[, εZT2× b QM |bπε(t)|3/2dx0dy =εZT2× b QMπεt, εκ x0 ε+εy0, εy3 3/2 dx0dy =ε3X k0∈KεZb QMπεt, εk0+εy0, εy33/2dy =X k0∈KεZCk0 ε×(0,εM)πεt, x0, x33/2dx ≤ZΩε |πε(t, x)|3/2dx. Taking the essential supremum in the former inequality, we deduce that kε2/3bπεkL∞(]0,T [,L3/2(T2× b QM,R)) ≤C. Since L3(T2×b QM,R) is separable, the space L1(]0, T[, L3(T2×b QM,R)) is separable. Hence, we obtain the convergence (4.17) by applying Banach Alaoglu theorem and a diagonal argument. Finally, for any bϕ∈C1 c(]0, T[×T 2, C1 #(b Q, R3)), the limit (4.18) comes from the following computation: lim ε→0"ε−1/3Z(0,T)×Ωε qε(t, x) (divybϕ)(t, x0,x ε)dx dt# = lim ε→0"ε−1/3Z(0,T)×Ω πε(t, x) (divy(∂tbϕ))(t, x0,x ε)dx dt# = lim ε→0"Z(0,T)×T 2× b QM ε2/3bπε(divy(∂tbϕ)) dx0dy dt# =Z(0,T)×T 2× b QMbπdivy(∂tbϕ)dx0dy dt =Z(0,T)×T 2× b QMbqdivybϕ dx0dy dt. Lemma 4.5. We consider a sequence of functions wε∈Jε, where Jεis defined by (2.22). Assume that there exists a constant C > 0, such that kwεkH1(]0,T[,L2(Ωε,R3)) ≤C, kwεkL3(]0,T [,W1,3(Ωε,R3)) ≤C. Assume that wε|Ωconverges weakly to a function win L3(]0, T[, W1,3(Ω,R3)) (this always holds for a subsequence). Then, the third component w3of wvanishes on ]0, T[×Γ. Moreover, the following result holds. (i) If aε≈ε5 3, with aε/ε5 3→λ,0<λ<+∞, then there exists bw∈L3(]0, T[, L3(T2,V3)) with bw3(t, x0, y0,0) = −λ∇Ψ(y0)·w0(t, x0,0),a.e. (t, x0, y0)∈]0, T [×T 2×Y0,(4.19) 17 such that for every M > 0, the sequence bwεdefined by (4.14) satisfies ε−2 3Dybwε* Dybwin L3(]0, T[×T 2×b QM,R9).(4.20) Besides, if div wε= 0 in ]0, T [×Ωε, then divybw= 0 in ]0, T[×T 2×b Q. (4.21) (ii) If aεε5 3, then w0(t, x0,0) · ∇Ψ(y0) = 0 a.e. (t, x0, y0)∈]0, T[×T 2×Y0.(4.22) Proof. We divide the proof in five steps. Step 1. Let us first prove that w3vanishes on Γ. Since wε·ν= 0 on ]0, T[×Γε, for every ϕ∈C1 c(Λ,R) we have ZΩε wε∇ϕ dx =−ZΩε (div wε)ϕ dx (4.23) in D0(0, T). Using ZΩε\Ω wε∇ϕ dx≤ZΩε |wε|3dx1 3 ZΩε\Ω |∇ϕ|3 2dx!2 3 →0, ZΩε\Ω (div wε)ϕ dx≤ZΩε |div wε|3dx1 3 ZΩε\Ω |ϕ|3 2dx!2 3 →0, and the weak convergence of wεto win L3(]0, T[, W1,3(Ω,R3)), we can pass to the limit in (4.23) to deduce ZΩ w· ∇ϕ dx =−ZΩ (div w)ϕ dx, in D0(0, T), and then ZΓ w3ϕ dx0= 0,∀ϕ∈C1 c(Λ,R), which proves w3= 0 on ]0, T[×Γ. Step 2. Let us obtain some estimates for the sequence bwεgiven by (4.14). For M > 0, the definition (4.14) of bwεproves for every ε > 0 small enough ZT2× b QM |Dybwε(t, x0, y)|3dx0dy =ε5X k0∈KεZb QM |Dwε(t, ε(k0+y0), εy3)|3dy =X k0∈Kε ε2ZCk0 ε×(0,εM) |Dwε|3dx ≤ε2ZΩε |Dwε|3dx. (4.24) 18 Integrating between 0 and T, and taking into account (3.1)-(3.2), we have ZT 0ZT2× b QM |Dybwε(x0, y)|3dx0dydt ≤Cε2.(4.25) On the other hand, defining ¯wε(t, x0) = 1 ε2ZCε(x0) wε(t, τ0,0) dτ0=ZY0bwε(t, x0, y0,0) dy0,(4.26) and using the inequality Zb QM |bwε(t, x0, y)−¯wε(t, x0)|3dy ≤CMZb QM |Dybwε|3dy, a.e. (t, x0)∈]0, T[×T 2, where CMdoes not depend on εand integrating between 0 and Twe have ZT 0Zb QM |bwε(x0, y)−¯wε(x0)|3dy ≤CMZT 0Zb QM |Dybwε|3dy, a.e. x0∈ T 2,(4.27) Taking into account (4.25), we deduce that c Wε=bwε(t, x0, y)−¯wε(t, x0) ε2 3 is bounded in L3(]0, T[×T 2, W1,3(b QM,R3))),∀M > 0.(4.28) Thus, there exists bw: ]0, T[×T 2×b Q→R3, such that, up to a subsequence, c Wε*bwin L3(]0, T[×T 2, W1,3(b QM,R3)),∀M > 0,(4.29) and then ε−2 3Dybwε* Dybwin L3(]0, T[×T 2×QM,R9),∀M > 0.(4.30) Passing to the limit by semicontinuity in inequalities (4.25) and (4.27) (this latest one after integration in T2), we get ZT 0ZT2× b QM |Dybw|3dx0dydt ≤C, ZT 0ZT2× b QM |bw|3dx0dydy ≤CM, and by the arbitrariness of M, once we prove the Y0-periodicity of bwin y0(Step 3), then bw∈L3(]0, T[×T 2,V3).(4.31) Moreover, if we also assume that div wε= 0 in ]0, T[×Ωε, then by definition (4.14) of bwε, we have divybwε= 0 in ]0, T[×T 2×b QM, which together with (4.30) proves divybw= 0 in ]0, T[×T 2×b Q. (4.32) Step 3. Let us prove that bwis Y0-periodic in y0. We observe that by definition (4.14) of bwε, for every M > 0, we have bwε(t, x1+ε, x2,−1 2, y2, y3) = bwε(t, x0,1 2, y2, y3),a.e. (t, x0, y2, y3)∈]0, T[×T 2×(−1 2,1 2)×(0, M). 19 Therefore the sequence c Wεsatisfies c Wε(t, x1+ε, x2,−1 2, y2, y3)−c Wε(t, x0,1 2, y2, y3) = −¯wε(t, x1+ε, x2) + ¯wε(t, x0) ε2 3 ,(4.33) a.e. (t, x0, y2, y3)∈]0, T[×T 2×(−1 2,1 2)×(0, M). where ¯wεis defined by (4.26). Using the fact that wε(t, ·) is bounded in L3(Γ)3, we can proceed as the proof of Lemma 4.6 in [43] to deduce that the right-hand side of this equality tends to zero in the sense of distributions in T2. Therefore, passing to the limit in (4.33) by (4.29), and taking into account the arbitrariness of Mwe get bw(t, x0,−1 2, y2, y3)−bw(t, x0,1 2, y2, y3) = 0 a.e. (t, x0, y2, y3)∈]0, T[×T 2×(−1 2,1 2)×R. Analogously, we can prove bw(t, x0, y1,−1 2, y3)−bw(t, x0, y1,1 2, y3) = 0 a.e. (t, x0, y1, y3)∈]0, T[×T 2×(−1 2,1 2)×R. These equalities prove that bwis periodic with respect to Y0. Step 4. Since wεand ∂twεare uniformly bounded in L3(]0, T[, W1,3(Ω,R3)) and L2(]0, T[×Γ,R3) respectively, by compactness of the embedding of W1,3(Ω) into L2(Γ), one can apply Aubin-Lions theorem to deduce that wεis precompact in L3(]0, T[, L2(Γ,R3)). Up to a subsequence, this implies the strong convergence of ¯wεto w(t, x0,0) in L3(]0, T[, L2(Γ,R3)). Thus, by (4.28), we deduce bwε(t, x0, y)→w(t, x0,0) in L3(]0, T[×T 2, W1,3(b QM,R3)),∀M > 0.(4.34) Step 5. Using the change of variables (4.16), which defines bwε, in the equality wεν= 0 on Γε, we get aε ε∇Ψ(y0)·bw0 ε(t, x0, y0,−aε εΨ(y0)) + bwε,3(t, x0, y0,−aε εΨ(y0)) = 0,(4.35) a.e. in ]0, T [×T 2×Y0. Thanks to (4.35) and (4.25), we have then aε ε∇Ψ(y0)·bw0 ε(t, x0, y0,0) + bwε,3(t, x0, y0,0) ≤Z0 −aε εΨ(y0)aε ε∇Ψ(y0)·∂y3bw0 ε(t, x0, y0, s) + ∂y3bwε,3(t, x0, y0, s)ds ≤Caε ε2 3 Z0 −aε εΨ(y0) |∂y3bwε(t, x0, y0, s)|3ds!1 3 a.e. (t, x0, y0)∈]0, T[×T 2×Y0. Taking the power 3, integrating in T2×Y0and using (4.25) we then deduce ZT2×Y0aε ε∇Ψ(y0)·bw0 ε(t, x0, y0,0) + bwε,3(t, x0, y0,0)3dx0dy0≤Ca2 ε, 20 a.e. in ]0, T [, which implies ZT2×Y0 aε ε∇Ψ(y0)·bw0 ε(t, x0, y0,0) + bwε,3(t, x0, y0,0) −ZY0aε ε∇Ψ(z0)·bw0 ε(t, x0, z0,0) + bwε,3(t, x0, z0,0)dz0 3 dx0dy0≤Ca2 ε, a.e. in ]0, T[. Integrating between 0 and Tand dividing by ε2, and taking into account that ∇Ψ has mean value zero in Y0, we get ZT 0ZT2×Y0 aε ε5 3 ∇Ψ(y0)·bw0 ε(t, x0, y0,0) −aε εZY0 ∇Ψ(z0)·bw0 ε(t, x0, z0,0) −¯w0 ε(t, x0) ε2 3dz0 +bwε,3(t, x0, y0,0) −¯wε,3(t, x0) ε2 3 3 dx0dy0dt ≤Caε ε2→0. (4.36) Thanks to (2.1) and taking into account the definition (3.5) of λ, we deduce: •If λ= +∞, statement (4.36) shows that aε ε5 3 ∇Ψ(y0)·bw0 ε(t, x0, y0,0) is bounded in L3(]0, T[×T 2×Y0,R), and then that ∇Ψ(y0)·bw0 ε(t, x0, y0,0) tends to zero in L3(]0, T [×T 2×Y0,R). By (4.34), this proves (4.22). •If λ∈(0,+∞), proceeding as in Step 4 we can pass to the limit in (4.36) to deduce (4.19). Proof of Theorem 3.3. First of all, thanks to (3.2), convergences (3.6) hold. The proof will be carried out by using suitable test functions ϕεin the variational formulation (3.4). We will divide the proof in several steps. Step 1. Obtaining the limit system. We consider φ∈C1 c([0, T],R), ϕ∈C1 c(Λ,R)3,ϕ3(x0,0) = 0, bϕ∈C1 c(T2, C1 #(b Q, R3)), with Dybϕ(x0, y) = 0 a.e. in {y3> M}, for some M > 0, such that    ϕ(x) = ϕ(x0,0) if x3≤0, ϕ3(x0,0) = 0,  bϕ(x0, y) = bϕ(x0, y0,0) if y3≤0, bϕ3(x0, y0,0) = −λ∇Ψ(y0)·ϕ0(x0,0). (4.37) Besides, we take ζ∈C∞(R) such that ζ(x3) =      1 if x3<1 3 0 if x3>2 3, (4.38) and Rε>0 such that Rε→ ∞, Rε  aε ε5 3 −λ 3 +aε ε3!→0.(4.39) 21 Then, we define ϕε∈W1,3(Ωε)3by                      ϕ0 ε(x) = ϕ0(x) + ε2 3bϕ0(x0,x ε)ζ(x3) ϕε,3(x) = ϕ3(x) + ε2 3"bϕ3(x0,x ε)ζ(x3)+ +ζ(x3 εRε )∇Ψ(x0 ε)·λ−aε ε5 3ϕ0(x0,0) −aε εbϕ0(x0,x0 ε,0)#. Thanks to ϕ(x) and bϕ(x0, y) equal zero for x0outside a compact subset of T2and (4.37), the sequence vε(t, x) = φ(t)ϕε(x) satisfies vε= 0 on ]0, T[×T0, vε·ν= 0 on ]0, T[×Γε. Thus, we can take such vεas test function in (3.4). The problem is to pass to the limit in the different terms which appears in (3.4). Before, we remark that since Dybϕ= 0 a.e. in {y3> M}and that ζ= 1 a.e. on {x3<1/3}, we have Dϕε(x) = Dϕ(x) + ε−1 3Dybϕ(x0,x ε) + hε(x),(4.40) where, using that ϕ,bϕand ζare bounded and have bounded derivatives, the function hε∈C0([0, T]× ¯ Ωε,R9) satisfies |hε| ≤ Cε2 3+C" 1 ε1 3Rε +1 ε1 3! aε ε5 3 −λ+aε ε+ε2 3 aε ε5 3 −λ#χ{x3<2 3εRε}, a.e. in Ωε. Using that Rεtends to infinity, we get |hε| ≤ Cε2 3+C1 ε1 3 aε ε5 3 −λ+aε ε1+1 3χ{x3<2 3εRε},a.e. in Ωε. Therefore, by (4.39), we have ZΩε |hε|3dx ≤Oε+CRε  aε ε5 3 −λ 3 +aε ε3!=Oε.(4.41) Taking vεas test function in (3.4) and using that kwεkL3(]0,T [,W1,3(Ωε,R3)) and kqεkL∞(]0,T [,L 3 2(Ωε,R)) are bounded, kϕε−ϕkC0(¯ Ωε,R)3tends to zero, (4.40) and (4.41), we get ZT 0 φ(t)ZΩε ∂twε(t)·ϕ dxdt +ZT 0ZΩε (wε(t)· ∇)wε(t)ϕ dxdt +ZT 0 φ(t)ZΩε S(Dwε(t)) : Dϕ +ε−1 3Dybϕ(x0,x ε)dxdt −ZT 0 φ(t)ZΩε qεdiv ϕ+ε−1 3divybϕ(x0,x ε)dxdt +βZT 0 φ(t)ZΓε wε·ϕ dH2dt =ZT 0 φ(t)ZΩε f·ϕ dxdt +Oε. (4.42) 22 Let us now pass to the limit in every term of the above equality. •First term of (4.42). Since |Ωε\Ω| → 0 and ∂twε∈L∞(]0, T[, L2(Ωε,R3)), then we have ZΩε ∂twε·ϕ dx =ZΩ ∂twε·ϕ dx +Oε. Moreover, we have ZΩε ∂twε(t)·ϕ dx =d dt ZΩε wε(t)·v dx, and so d dt ZΩε wε(t)·ϕ dxφ(t) = d dt φ(t)ZΩε wε(t)·v dx−∂tφ(t)ZΩε wε(t)·v dx. Thus, we get ZT 0 φ(t)ZΩε ∂twε(t)·v dx dt =−ZT 0 ∂tφ(t)ZΩε wε(t)·ϕ dx dt. Since wεconverges weakly-∗to win L∞(]0, T[, L2(Ω,R3)) and taking ϕ ∂tφ(t)∈L1(]0, T[, L2(Λ,R3)), we have ZT 0 φ(t)ZΩε ∂twε(t)·v dx dt =−φ(0) ZΩε w0·ϕ dx −ZT 0 ∂tφ(t)ZΩε w(t)·ϕ dx dt +Oε =ZT 0 φ(t)ZΩε ∂tw(t)·v dx dt +Oε. (4.43) •Second term of (4.42). From div wε= 0 and wεν= 0 on Γε, we have ZT 0ZΩε (wε(t)· ∇)wε(t)φ(t)ϕ dx =−ZT 0 φ(t)ZΩε (wε(t)· ∇ϕ)wεdxdt. Taking into account that |Ωε\Ω| → 0, and wε∈Jε, we have ZT 0ZΩε (wε(t)· ∇)wε(t)φ(t)ϕ dx =−ZT 0 φ(t)ZΩ (wε(t)· ∇ϕ)wε(t)dxdt +Oε. Since wε∈L3(]0, T[, W1,3(Ω,R3)) and ∂twε∈L2(]0, T[×Ω,R3), by applying Aubin-Lions theorem we deduce that wεis compact in L3(]0, T[, L3(Ω,R3)) wich implies the strong convergence of wεto win L3(]0, T[, L3(Ω,R3)). This together with weak convergence of wεto win L3(]0, T[, W1,3(Ω,R3)) implies ZT 0ZΩε (wε(t)· ∇)wε(t)φ(t)ϕ dx =−ZT 0 φ(t)ZΩ (w(t)· ∇ϕ)w(t)dxdt +Oε.(4.44) •Third term of (4.42). Taking into account that ZΩε\ΩDϕ +ε−1 3Dybϕ(x0,x ε)3dx ≤C ε|Ωε\Ω| ≤ Caε ε=Oε,(4.45) 23 we divide the term in two integrals ZT 0 φ(t)ZΩε S(Dwε(t)) : Dϕ +ε−1 3Dybϕ(x0,x ε)dxdt =ZT 0 φ(t)ZΩ S(Dwε(t)) : Dϕ dxdt +ε−1 3ZT 0 φ(t)ZΩ S(Dwε(t)) : Dybϕ(x0,x ε)dxdt +Oε. (4.46) Since wε∈L3(]0, T[, W1,3(Ωε)3), then S(Dwε) is bounded in L3 2(]0, T[, L3 2(Ωε,R9)), and so there exists ξ∈L3 2(]0, T[, L3 2(Ω,R9)) such that S(Dwε)* ξ in L3 2(]0, T[, L3 2(Ω,R9)).(4.47) Then the first term in the right-hand side of (4.46) reads ZT 0 φ(t)ZΩ S(Dwε) : Dϕ dxdt =ZT 0 φ(t)ZΩ ξ:Dϕ dxdt +Oε.(4.48) In the second term of the right-hand side of (4.46), we introduce the sequence bwεdefined by (4.14). By wε∈L3(]0, T[, W1,3(Ωε,R3)) and Lemma 4.5, we have that b S(ε−2 3Dybwε) is bounded in L3 2(]0, T[, L3 2(T2× b QM,R9)) ∀M > 0, and so there exists b ξ∈L3 2(]0, T[, L3 2(T2×b Q, R9)) such that b S(ε−2 3Dybwε)*b ξin L3 2(]0, T[, L3 2(T2×b Q, R9)). Moreover, we have ε−1 3ZT 0 φ(t)ZΩ S(Dwε) : Dybϕ(x0,x ε)dxdt =ε−1 3 Re ZT 0 φ(t)ZΩ Dwε:Dybϕ(x0,x ε)dxdt +ε−1 3ZT 0 φ(t)ZΩb S(Dwε) : Dybϕ(x0,x ε)dxdt =ε−1 3 Re ZT2ZT 0 φ(t)Zb QM Dybwε:Dybϕ dydx0dt +ε2 3ZT 0 φ(t)ZT2Zb QMb S(ε−1Dybwε) : Dybϕ dydx0dt +Oε =ε1 3 Re ZT 0 φ(t)ZT2Zb QM ε−2 3Dybwε:Dybϕ dydx0dt +ZT 0 φ(t)ZT2Zb QMb S(ε−2 3Dybwε) : Dybϕ dydx0dt +Oε =ZT 0 φ(t)ZT2Zb QMb S(ε−2 3Dybwε) : Dybϕ dydx0dt +Oε =ZT 0 φ(t)ZT2× b Qb ξ:Dybϕ dx0dydt +Oε. Thus, (4.48) and the above equality imply ZT 0 φ(t)ZΩε S(Dwε) : Dϕ +ε−1 3Dybϕ(x0,x ε)dxdt =ZT 0 φ(t)ZΩ ξ:Dϕ dxdt +ZT 0 φ(t)ZT2× b Qb ξ:Dybϕ dx0dydt +Oε. (4.49) 24 •Fourth term of (4.42). Similarly to the previous term, taking into account (4.45), we split this term as follows ZT 0ZΩε qε(t)div ϕ+ε−1 3divybϕ(x0,x ε) dxdt =ZT 0ZΩ qε(t) div ϕ dxdt +ε−1 3ZT 0ZΩ qε(t) divybϕ(x0,x ε)dxdt +Oε. (4.50) From (3.6) in the first term in the right-hand side of (4.50), and (4.18) for the second term, then we have ZT 0ZΩ qε(t)div ϕ+ε−1 3divybϕ(x0,x ε)dx =ZT 0ZΩ q(t) div ϕ dx +ZT 0ZT2× b Qbq(t) divybϕ dx0dy +Oε. (4.51) •Fifth term of (4.42). In order to identify the limit of the boundary term, we write βZT 0 φ(t)ZΓε wε(t)·ϕ dH2 =βZT 0 φ(t)ZT2 wε(t)t, x0,−aεΨ(x0 ε)·ϕx0,−aεΨ(x0 ε)s1 + aε ε2∇Ψ(x0 ε) 2 dx0. (4.52) The inequality ZT2wεt, x0,−aεΨ(x0 ε)−wε(t, x0,0) 3 dx0=ZT2Z0 −aεΨ( x0 ε) ∂3wε(t, x0, s)ds 3 dx0 ≤Ca2 εZΩε |∂3wε(t)|3dx =Oε, (4.53) the strong convergence of wε(t, x0,0) to w(t, x0,0) in L3(]0, T [, L3(T2,R3)) and the fact that aε/ε tends to zero, imply βZT 0 φ(t)ZΓε wε(t)·ϕ dH2dt =βZT 0 φ(t)ZΓ w(t)·ϕ dH2dt +Oε.(4.54) •Last term of (4.42). Since |Ωε\Ω| → 0, we have ZT 0 φ(t)ZΩε f·ϕ dxdt =ZT 0 φ(t)ZΩ f·ϕ dxdt +Oε.(4.55) Therefore, substituting (4.43), (4.44), (4.49), (4.51), (4.54) and (4.55) in (4.42), we get that the limit problem is equivalent to the following one ZT 0ZΩ ∂tw(t)·ϕ dxdt +ZT 0ZΩ (w(t)· ∇)w(t)ϕ(t)dxdt +ZT 0ZΩ ξ(t) : Dϕ(t)dxdt +ZT 0ZΩ q(t) div ϕ(t)dxdt +ZT 0ZT2× b Qb ξ(t) : Dybϕ(t)dx0dydt +ZT 0ZT2× b Qbq(t) divybϕ(t)dx0dydt +βZT 0ZΓ w(t)·ϕ(t)dH2dt =ZT 0ZΩ f·ϕ(t)dxdt , (4.56) 25 Proof. Let ξ0 1, ξ0 2∈R2. Developing the scalar product and using definition (3.9), we obtain (R(ξ0 1)−R(ξ0 2)) ·(ξ0 1−ξ0 2) = (R1(ξ0 1)−R1(ξ0 2))(ξ0 1−ξ0 2)1+ (R2(ξ0 1)−R2(ξ0 2))(ξ0 1−ξ0 2)2 =Zb Qb S(Dyb φξ0 1)−b S(Dyb φξ0 2):(ξ0 1−ξ0 2)1Dyb φe1+ (ξ0 1−ξ0 2)2Dyb φe2dy =Zb Qb S(Dyb φξ0 1)−b S(Dyb φξ0 2):Dy(ξ0 1−ξ0 2)1b φe1+ (ξ0 1−ξ0 2)2b φe2dy. Now, defining G∈(V/R)2×V ×L 3 2 #(b Q, R) as the combination G= (ξ0 1−ξ0 2)1b φe1+(ξ0 1−ξ0 2)2b φe2, we notice that Gsatisfies the same boundary condition as the difference b φξ0 1−b φξ0 1, namely, G3(y0,0) = −∇Ψ(y0)·(ξ0 1−ξ0 2) = b φξ0 1 3(y0,0) −b φξ0 2 3(y0,0) for a.e. y0∈ T 2. Consequently, G−(b φξ0 1−b φξ0 1) is an admissible test function for problem (3.10), with ξ0∈ {ξ0 1, ξ0 2}. Arguing as in the proof of Proposition 3.8, we further deduce that (R(ξ0 1)−R(ξ0 2)) ·(ξ0 1−ξ0 2) = Zb Qb S(Dyb φξ0 1)−b S(Dyb φξ0 2):Dy(b φξ0 1−b φξ0 2)dy =Zb Qb S(Dyb φξ0 1)−b S(Dyb φξ0 2):Dyb φξ0 1−Dyb φξ0 2dy, which is nonnegative by monotonicity of b S. 5 Numerical results In this section, we investigate numerically the behaviour of the solution of the LES model (2.14)–(2.16), in the case where the lower plate is composed of riblets, which are very commonly used in experimental and theoretical studies of the roughness effects on viscous flows, in particular for turbulent regimes (see, for instance, [8, 22, 25]). In that case, in our framework, the sequence of functions Ψεdescribing the geometry of the rough wall, is defined by Ψε(x0) = aεΨ(x1/ε), x0∈ T 2,(5.1) where Ψ is a given, positive 1-periodic and regular function of variable x1, and where the sequence (aε)ε>0satisfies assumptions (2.1). By Corollary 3.6, the effective velocity and pressure w, q satisfy the limit system (3.7), completed with the boundary conditions (3.11) on the (idealized) flat lower plate Γ, that express the effect of the asperities on the LES model, at main order in ε. In order to stress the influence of the parameter λ, defined by (3.5), we simulate stationary solutions of the limit system (3.7), which allows us to aggregate numerous numerical results associated with a wide spectrum of values of λ. In the case of the riblet geometry described by (5.1), we seek the velocity wand the pressure qsatisfying the following effective system:                      (w· ∇)w−div (S(Dw)) + ∇q=fin Ω, div w= 0 in Ω, w= 0 on Γ0, w3= 0 on Γ, −S(Dw)1,3+β w1+λ3w2 1R1(sgn(w1)e1) = 0 on Γ, −S(Dw)2,3+β w2= 0 on Γ. (5.2) 32 To investigate the influence of the shape of the riblets on the effective flow, we will focus on three different profiles for the function Ψ, that correspond to the V-shaped, the U-shaped and the blade-shaped riblets. These shapes are widely used in the engineering literature, and are drawn in Table 1. Notice that, for each profile, we have applied a parabolic regularization in the vicinity of each singular point, obtaining functions Ψ of class C2. Table 2 gives an example of transversal sections of ribbed walls Γε, associated with those different profiles Ψ by the definition Γε=(x1, x2,−aεΨ(x1/ε)) ,(x1, x2)∈ T 2. We stress that, in the case of riblets, according to boundary conditions (3.11), one needs only to compute the coefficient R1(ξ0) for ξ0=e1and ξ0=−e1, and to have access to the sign of w1(x0,0), to define the boundary condition at point x0∈ T 2. Let us first address the numerical approximation of coefficients R1(e1), R1(−e1). Figure 1: Functions Ψ associated with V-shaped, U-shaped and blade riblets (from left to right). Figure 2: Sections of ribbed walls Γε, associated with a function Ψ modelling the V-shaped, U-shaped and blade-shaped riblet profile (from left to right). 5.1 Numerical methods and algorithms Resolution of the auxiliary system (3.10).In order to compute coefficients R1(e1), R1(−e1), one needs to solve system (3.10), with ξ0∈ {e1,−e1}. We fix a truncature parameter H > 0 and we consider the 33 following approximated problem: find (b φξ0 H,bπξ0 H)∈W1,3(T2×(0, H),R3)×L3/2(T2×(0, H)) satisfying                        −divy(b S(Dyb φξ0 H)) + ∇ybπξ0 H= 0 in T2×(0, H), divyb φξ0 H= 0 in T2×(0, H), b φξ0 H,3(y0,0) = −∇Ψ(y0)·ξ0on T2× {0}, b S(Dyb φξ0 H)i,3= 0, i = 1,2,on T2× {0}, b S(Dyb φξ0 H)i,3= 0, i = 1,2,3,on T2× {H}. (5.3) Then, the approximate values Rapp 1(ξ0) of R1(ξ0), for ξ0∈ {e1,−e1}, will be given by Rapp 1(ξ0) = ZT2×(0,H)b S(Dyb φξ0 H) : Dyb φe1 Hdy. Such an approximation can be justified by standard arguments, relying on Proposition 3.8 that states the exponential decay of the gradient of b φξ0, along the y3-direction. Since problem (5.3) is nonlinear (due to the tensor b S), we introduce a fixed point algorithm, that is defined as follows. Here, we drop the ξ0exponents and the Hsubscripts to lighten the notation. Assume that b φk∈W1,3(T2×(0, H),R3), and satisfies divyb φk= 0 in T2×(0, H),b φ3,k(y0,0) = −∇Ψ(y0)·ξ0on T2× {0}.(5.4) Then, b φk+1 is determined as the only solution in W1,3(T2×(0, H),R3), satisfying the same incompressibility and boundary condition (5.4), to the following variational formulation: ZT2×(0,H)Cs|Db φk|b Dφk+1 −δ2 2γDb φk(Db φk+1)T:Dζ dy = 0,(5.5) for any test function ζ∈W1,3(T2×(0, H),R3) satisfying divyζ= 0 in T2×(0, H), ζ3(y0,0) = 0 on T2× {0}. We stop the iterations once the stopping criterion kDbϕk+1 −DbϕkkL3(T2×(0,H),R3) kDbϕkkL3(T2×(0,H),R3) <  is fulfilled, where is a small number (set at = 10−2in the numerical simulations); then, we use the approximation b φξ0≈b φk+1. Resolution of the effective problem (5.2).Denoting by Vthe functional space V={w∈W1,3(Ω,R3),div w= 0 in Ω, w3= 0 on Γ, w = 0 on Γ0}, endowed with the norm kwkV:= kDwk1/3 L3(Ω,R3×3) the variational formulation of problem (5.2) reads: for every ψ∈V, ZΩ(w· ∇)w·ψ+S(Dw) : Dψdx +ZΓβ w ·ψ+λ3w2 1R1(sgn(w1)e1)ψ1dH2=ZΩ f·ψ dx. (5.6) 34 To deal with the nonlinear terms appearing in the previous formulation, we proceed in an analogous way as in the treatment of system (5.3), by introducing a fixed point algorithm. Starting from an initial guess w0∈V, we construct a sequence of functions wk∈Vdefined recursively, as follows. For k∈N,wkbeing known, we define wk+1 ∈Vas the solution to the following variational formulation: ZΩ(wk· ∇)wk+1 ·ψ+S(Dwk+1) : Dψdx +ZΓβ wk+1 ·ψ+λ3wk,1wk+1,1R1(sgn(wk,1)e1)ψ1dH2 =ZΩ f·ψ dx, ∀ψ∈V. (5.7) Numerically, we stop the iterations if the relative error kwk+1 −wkkV/kwkkVis smaller than the threshold defined in the previous paragraph. Spatial discretization by finite element methods. In order to compute the solutions to the variational formulations (5.5) and (5.7), we rely on mixed formulations of the problems, that are discretized by finite element methods using FreeFem++ software [27]. For both systems of equations, we consider the Taylor- Hood approximation for the velocity-pressure pair, i.e.,P2elements for the velocity field and P1elements for the pressure. It is well known that this choice is compatible with the Babuˇska-Brezzi condition [12,26]. The three-dimensional meshes of the cell T2×(0, H) and the fluid domain Ω are obtained by constrained Delaunay tetrahedralizations. 5.2 Numerical results Definition of the parameters associated with the turbulence model. In order to illustrate the effect of coefficient λ, which acts directly on the component of the fluid velocity that is transverse to the riblet direction, we consider a flow driven by an external force f=e1. We set the Reynolds number to Re = 10000, a value for which the use of Large Eddy Simulation may be justified in practice. For instance, in Reynolds’ pipe experiment, the flow is turbulent if the Reynolds number exceeds about 4000 (see [39], p. 6). We choose δas the grid size of the tetrahedralization of domain Ω, and define Cs=δ2and γ= 6, which appears to be standard for the numerical simulation of Smagorinsky’s turbulent model (see for instance [46]). In particular, such definition of δis coherent with the modelling of the turbulence developed in Sec. 2, since δis a small averaging radius and the grid size is the smallest length scale accessible in the mesh. To set the value of parameter β, we follow an argument developed by Berselli et al. in [10], and define β=L Re δ ,(5.8) where Lis a characteristic length of the domain (in our case, we take L= 1). To obtain formula (5.8), the previous authors have adapted Maxwell’s theory of gases [36], considering the micro length scale δas the analog of the mean free path of gas particles. Computation of the coefficients R1(e1), R1(−e1)associated with the different shapes of riblets. We have collected in Table 1 the numerical values obtained for the coefficient R1(e1), R1(−e1), for the different riblet geometries. In our simulations, we have defined H= 3 – a relatively small value that is justified by the exponential decay of the derivatives of b φin the vertical direction. In order to be consistent with estimate (3.12) stated in Proposition 3.8, we have normalized every function Ψ so that ZT1 |Ψ0(y1)|3dy1= 1. We stress that the absolute values of those coefficients depend strongly on the geometry of the riblets, through the choice of the function Ψ, even though the L3-norm of Ψ0is normalized. For instance, the 35 V-shape coefficient is about 40 times smaller as the U-shape coefficient, and 80 times smaller as the blade geometry one. Another interesting feature is that, even if the value R1(−e1) is not numerically equal to the opposite of R1(e1), the difference appears only at the fourth digit (or at the sixth digit in the V-shape case). Such result may be explained by our choice of parameters, for which parameter CSis 12 times bigger than the factor δ2/(2γ). If one neglects that term, one can easily verify that the solution to system (5.3) satisfies b φ−e1=−b φe1, hence R1(−e1) = −R1(e1). So, we hypothesize that, with our choice of parameters, the nonlinear tensor b S(b φξ) is dominated by the term Cs|φξ|φξ, which explains such apparent symmetry. V-shape U-shape Blade R1(e1) 0.20911 4.8220 8.6947 R1(−e1) -0.20911 -4.8294 -8.6903 Table 1: Approximate values of the coefficients R1(e1), R1(−e1), for the U-shape, V-shape and blade riblets. Numerical results and comments. In order to establish some qualitative comparisons between the solutions to the effective system (5.2), let us introduce certain quantities that are classically associated with the stationary flow of a viscous fluid, in the vicinity of a solid wall where friction occurs: •the kinetic energy 1 2RΩ|w|2dx ; •the dissipated energy 1 2ZΩ1 Re|Dw|2+Cs|Dw|3−δ2 2γDw(Dw)T:Dwdx +ZΓβ|w|2+λ3w3 1R1(sgn(w1)e1)dH2, that we decompose into two contributions: a) the viscous dissipation 1 2ZΩ1 Re|Dw|2+Cs|Dw|3−δ2 2γDw(Dw)T:Dwdx, b) the dissipation by friction 1 2ZΓβ|w|2+λ3w3 1R1(sgn(w1)e1)dH2; •the work of the body force f, which in our case is simply defined by ZΩ w1dx, and will be refered to as the average streamwise velocity. We have plotted the kinetic and dissipated energies, the viscous and frictional dissipation, and the average streamwise velocity, with respect to the value of λ, and considering V-shape, U-shape and blade riblets, in Figures 3 to 7. 36 0 0.511.5 2 0 10 20 30 40 50 λ Kinetic energy V-shape U-shape blade Figure 3: Kinetic energy with respect to λ, for the different riblet geometries. We observe that for every riblet geometry, increasing λleads to decreasing the kinetic energy of the fluid, as well as the dissipated energy and the streamwise velocity. For λ= 0, the geometry of the asperity has no effect in the model, and everything happens as if the lower wall was flat, and produced a friction measured by the coefficient β. In that case, the kinetic energy is above 50. Then, as λincreases, the kinetic energy diminishes drastically, before stabilizing around the value 10, which is fifth time smaller as its initial value. The profiles of the dissipated energy, its viscous part, and the streamwise velocity, are similar, but their values are only divided by a factor 2 between λ= 0 and λ= 2.3. We notice that for every mechanical quantity that we have computed, the numerical values obtained for the U-shape and the blade riblets are very close, whereas the V-shape riblet behaves differently. In particular, the decreasing of the former quantities that appears for small values of λ, is much slower for the V-shape riblets than for the other ones. That can be explained by the proximity of the R1coefficients associated with the U-shape and the blade riblets: they only differ by a factor 2 (see Table 1), whereas the V-shape riblet produces R1coefficients 40 or 80 smaller in magnitude. Consequently, the effect of an increasing of λis less pronounced for the V-shape riblet, especially since the corresponding term in the variational formulation is weighted by λ3, which is increasing very slowly with λ. Among all the mechanical quantities introduced above, the frictional dissipation is the only one that presents a different behaviour. Considering the case of the V-shape riblet, we observe on Fig. 6 that the frictional dissipation start from about 0.3, and then increases to reach its maximal value (around 1) for λ≈0.5; from that point, the value decreases, mimicing the behaviour that the kinetic or dissipated energies presents from λ= 0. The origin of that initial increase of the frictional dissipation is difficult to analyze; however, Fig. 6 is an indicator of the complexity of the interaction between the fluid and the ribbed wall, which is described by a monotone term (see Remark 4.6), depending locally on the sign of the first component of the velocity. Still, it is worth noticing that in the present simulations, the viscous dissipation dominates the frictional one, so that the dissipated energy behaves similarly to the viscous one with respect to λ. Concluding remarks. In the case of straight riblets, where the boundary conditions associated with the effective system (3.7) are given by (3.11), we have proposed a methodology that allows to compute 37 0 0.5 1 1.5 2 2 3 4 λ Dissipated Energy V-shape U-shape blade Figure 4: Dissipated Energy with respect to λ, for the different riblet geometries. 0 0.511.5 2 2 3 4 λ Viscous Dissipation V-shape U-shape blade Figure 5: Viscous dissipation with respect to λ, for the different riblet geometries. 38 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 λ Frictional Dissipation V-shape U-shape blade Figure 6: Frictional dissipation with respect to λ, for the different riblet geometries. 0 0.5 1 1.5 2 4 6 8 10 λ Streamwise Velocity V-shape U-shape blade Figure 7: Average streamwise velocity w1, with respect to λ, for the different riblet geometries. 39 the required coefficients R1(e1), R1(−e1), and to solve the variational formulation (5.6) by a fixed point method. Such procedure could be easily extended to treat the case of a time-dependant body force f, adapting standard time discretization schemes for Navier-Stokes equations. Hence, the new boundary conditions that the asymptotic analysis detailed in Sec. 4 has pointed, are usable in practice, at least in the case of a ribbed geometry. Besides, the analysis of the characteristics of the flow with respect to the value of coefficient λ, for different types of riblets, reveals that the effective system can be strongly affected by the extra friction term associated with R1and λin boundary conditions (3.11). By construction, that term is a non physical one, that is, a term of interaction between the fluid and the ribbed wall resulting from the averaging procedure leading to the Large Eddy Simulation model, and not from the mechanical interaction ifself. As a matter of fact, that new friction term has the same mathematical properties as the tensor b Sthat results from the spatial averaging of the Navier-Stokes equations; that is clear from the computations following Remark 4.6, and could explain the large influence of that extra friction term on the macroscopic features of the simulations. Hence, such effective boundary conditions as (3.11) need certainly to be taken into account in the simulation of homogenized LES models, in the vicinity of a ribbed wall. Acknowledgements F.J. Su´arez-Grau has been supported by Ministerio de Econom´ıa y Competitividad (Spain), Proyecto Excelencia MTM2014-53309-P. References [1] N. A. Adams and S. Stolz, A Subgrid-Scale Deconvolution Approach for Shock Capturing, J. Comput. Phys. 178 (2002), no. 2, 391–426. [2] G. Allaire, Homogenization and Two-Scale Convergence, SIAM Journal on Mathematical Analysis 23 (1992), no. 6, 1482–1518. [3] Y. Amirat and J. Simon, Riblets and Drag Minimization, Contempory Mathematics 209 (1997), 9–17. [4] T. Arbogast, J. Douglas Jr., and U. Hornung, Derivation of the Double Porosity Model of Single Phase Flow via Homogenization Theory, SIAM Journal on Mathematical Analysis 21 (1990), no. 4, 823–836. [5] T. G. Bagwell, Stochastic Estimation of Near Wall Closure in Turbulence Models., Ph.D. thesis, University of Illinois at Urbana-Champaign, 1994. [6] G. I. Barenblatt and A. J. Chorin, New Perspectives in Turbulence: Scaling Laws, Asymptotics, and Intermittency, SIAM Review 40 (1998), 265–291. [7] D. Bechert and M. Bartenwerfer, The Viscous Flow on Surfaces with Longitudinal Ribs, J. Fluid Mech. 206 (1989), no. 1, 105–129. [8] D. W. Bechert, G. Hoppe, J. G. Th. van der Hoeven, and R. Makris, The Berlin Oil Channel for Drag Reduction Research, Experiments in Fluids 12 (1992), no. 4, 251 – 260. [9] L. C. Berselli, G. P. Galdi, T. Iliescu, and W. J. Layton, Mathematical Analysis for the Rational Large Eddy Simulation Model, Mathematical Models and Methods in Applied Sciences 12 (2002), no. 08, 1131–1152. 40 [10] L. C. Berselli, T. Iliescu, and W. J. Layton, Mathematics of Large Eddy Simulation of Turbulent Flows, Springer-Verlag Berlin Heidelberg, 2006. [11] H. Br´ezis, Analyse fonctionnelle : th´eorie et applications, Collection Math´ematiques appliqu´ees pour la maˆıtrise, Masson, Paris, New York, Barcelone, 1983. [12] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag New York, Inc., New York, NY, USA, 1991. [13] Dorin Bucur, Eduard Feireisl, and ˇ S´arka Neˇcasov´a, On the Asymptotic Limit of Flows Past a Ribbed Boundary, Journal of Mathematical Fluid Mechanics 10 (2008), no. 4, 554–568. [14] , Boundary Behavior of Viscous Fluids: Influence of Wall Roughness and Friction-driven Boundary Conditions, Archive for Rational Mechanics and Analysis 197 (2010), no. 1, 117–138. [15] J. Casado-D´ıaz, Exponential Decay for the Solutions of Nonlinear Elliptic Systems Posed in Unbounded Cylinders, Journal of Mathematical Analysis and Applications 328 (2007), no. 1, 151 – 169. [16] J. Casado-D´ıaz, E. Fern´andez-Cara, and J. Simon, Why Viscous Fluids Adhere to Rugose Walls: a Mathematical Explanation, Journal of Differential Equations 189 (2003), no. 2, 526 – 537. [17] J. Casado-D´ıaz, M. Luna-Laynez, and F. J. Su´arez-Grau, Asymptotic Behavior of a Viscous Fluid with Slip Boundary Conditions on a Slightly Rough Wall, Mathematical Models and Methods in Applied Sciences 20 (2010), no. 01, 121–156. [18] D. Cioranescu, A. Damlamian, and G. Griso, Periodic Unfolding and Homogenization, C. R. Acad. Sci. Paris 335 (2002), no. 1, 99 – 104. [19] P. Coletti, Analytical and Numerical Analysis of k-e and Large Eddy Simulation Turbulence Models, Ph.D. thesis, Dipartimento di Matematica, Universit`a di Trento, 1998. [20] B. Dean and B. Bhushan, Shark-Skin Surfaces for Fluid-Drag Reduction in Turbulent Flow: a Review, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 368 (2010), no. 1933, 5737–5737. [21] J. W. Deardorff, A Numerical Study of Three-Dimensional Turbulent Channel Flow at Large Reynolds Numbers, Journal of Fluid Mechanics 41 (1970), 453–480. [22] C.D. Dritselis, Large Eddy Simulation of Turbulent Channel Flow with Transverse Roughness Elements on One Wall, International Journal of Heat and Fluid Flow 50 (2014), 225 – 239. [23] G. P. Galdi, An Introduction to the Mathematical Theory of the Navier-Stokes Equations : Steady- State Problems, Springer Monographs in Mathematics, Springer, New York, NY, 2011. [24] G. P. Galdi and W. J. Layton, Approximation of the Larger Eddies in Fluid Motions II: a Model for Space-Filtered Flow, Mathematical Models and Methods in Applied Sciences 10 (2000), no. 03, 343–350. [25] R. Garc´ıa-Mayoral and J. Jim´enez, Drag Reduction by Riblets, Phil. Trans. R. Soc. A 369 (2011), 1412 – 1427. 41