On the influence of wavy riblets on the slip behaviour of viscous fluids

On the influence of wavy riblets on the slip behaviour of viscous fluids Matthieu Bonnivard, Francisco Suarez-Grau, Giordano Tierra To cite this version: Matthieu Bonnivard, Francisco Suarez-Grau, Giordano Tierra. On the influence of wavy riblets on the slip behaviour of viscous fluids. 2015. <hal-01161634> HAL Id: hal-01161634 https://hal.archives-ouvertes.fr/hal-01161634 Submitted on 9 Jun 2015 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destin´ee au d´epˆot et `a la diffusion de documents scientifiques de niveau recherche, publi´es ou non, ´emanant des ´etablissements d’enseignement et de recherche fran¸cais ou ´etrangers, des laboratoires publics ou priv´es. On the influence of wavy riblets on the slip behaviour of viscous fluids Matthieu Bonnivard∗,Francisco J. Su´ arez-Grau†&Giordano Tierra‡ June 7, 2015 Abstract In this work, we use the homogenization theory to investigate the capability of wavy riblet patterns to influence the behaviour of a viscous flow near a ribbed boundary. Starting from perfect slip conditions on the wall, we show that periodic oscillations of wavy riblets in the lateral direction may induce a friction effect in the direction of the flow, contrary to what happens with straight riblets. Finally, we illustrate this effect numerically by simulating riblet profiles that are widely used in experimental studies: the V-shape, U-shape, and blade riblets. Keywords: viscous fluids, slip condition, rough boundary, straight riblets, wavy riblets, homogenization. Mathematics subject classifications (2000): 35B27, 35Q35, 74A55. Contents 1 Introduction 2 2 Notation 5 3 Setting of the problem and main result 6 4 Proofs of the results 10 5 Numerical results 21 5.1 Numerical approximation of auxiliary system (3.14) ............. 21 5.2 Definition of the geometry of the wavy riblets ................. 22 5.3 Results and discussion .............................. 23 ∗Laboratoire Jacques-Louis Lions, Universit´e Paris Diderot - Paris 7, Paris, 75013, France, email: [email protected] †Dpto. Ecuaciones Diferenciales y An´alisis Num´erico, Universidad de Sevilla, Seville, 41012, Spain, e-mail: [email protected] ‡Mathematical Institute, Faculty of Mathematics and Physics, Charles University, Prague 8, 186 75, Czech Republic, e-mail: [email protected] 1 1 Introduction During the last decades, the reduction of the skin friction associated with a viscous flow in the vicinity of a solid boundary has been the subject to an abundant literature. From an experimental point of view, an efficient method to reduce the friction consists in using a certain type of rough boundaries, that are called riblets or straight riblets. These riblets are characterized by fast oscillations in the spanwise direction, with a low amplitude, and by their constancy in the streamwise direction; they are essentially onedimensional perturbations of the boundary of the solid. The spanwise geometry of the riblets is generally described by the periodic repetition of a concave pattern, typically, the shape of a V or a U. The aim of the use of riblets is to prevent vortices to appear in the neighborhood of the solid wall, and thus to reduce the momentum transfer from the vortices to the solid boundary. Recently, new methods have been investigated to enhance the efficiency of the riblets. One of these methods consists in a combination of riblets with an oscillating wall, that is achieved by tilting riblet blades sinusoidally. This approach has improved the drag reduction from 8.6% to 11.1% [30]. However, this result needs to be put in perspective, since the temporal oscillations of the boundary require a certain input of energy into the system. Another possible approach consists in using spatial oscillations of the riblets in the streamwise direction, thus replacing a temporal, active process by a static, passive one. This method was studied experimentally in [29], with a net power reduction increased up to 23%. The corresponding devices are called wavy riblets; they consist in straight riblets that are deflected sinusoidally in the spanwise direction, as a function of the streamwise position. Wavy riblets were studied in [25], both numerically and experimentally, in the case of a turbulent flow; however, the authors were not able to conclude regarding the efficiency of the device. In the present paper, we study the effect of wavy riblets on the slip behaviour of a viscous flow, when the amplitude and the periods of the spanwise and streamwise oscillations are small. Although the main field of application of wavy riblets concerns high Reynolds number flows, the present work deals with the case of low Reynolds number, in which the flow can be modeled via the linear Stokes equations. We believe that the presented results could be instrumental in the mathematical study of certain models for turbulence used in Large Eddy Simulation, that we aim to address in a future work. Definition of straight and wavy riblets. The ribbed walls that we consider are geometric perturbations of a planar boundary, which is identified with the plane R2× {0}. In the cartesian system of coordinates (x1, x2, x3), x3denotes the wall-normal coordinate. We define x2as the direction of the flow, or streamwise direction; x1is thus the lateral, or spanwise direction. We assume that the geometry of the riblets is associated with a small parameter ε, that controls the (spatial) oscillations of the boundary. More precisely, we name wavy riblets a parametric surface defined by a relation of the form x3=−Ψε(x1, x2), where Ψε is a positive function given by Ψε(x1, x2) := a1 εΨ1x1 p1 ε +a2 ε p1 ε Ψ2x2 p2 ε,(1.1) Ψ1, Ψ2being smooth and periodic functions of period 1. Parameters a1 ε, p1 ε>0 (resp. 2 a2 ε, p2 ε>0) are the amplitude and period of oscillation in the x1-x3axes (resp., in the x1-x2axes). Let us stress that the definition of wavy riblets (1.1) is a generalization of straight riblets, which correspond to the choice Ψ2≡0 (see Figure 1). Figure 1: Views from above of wavy riblets (left) and straight riblets (right) In order to study the rugosity effect of straight riblets endowed with perfect slip conditions, several theoretical works related to homogenization theory have been carried out. It was shown in [7] that, starting with perfect slip condition at the ribbed boundary given by Ψε(x1, x2) = εΨx1 ε, in the limit ε→0, no-slip condition appears in the spanwise direction while perfect slip still holds in the streamwise direction. This means that straight riblets tend to prevent the fluid from slipping laterally, whereas the motion in the streamwise direction is allowed with no constraint. The former work was generalized in [15] for riblets with low amplitude, i.e. Ψε(x1, x2) = δεΨx1 ε, where δεε. In that case, the limit velocity field satisfies a perfect slip condition in the streamwise direction, but its slippage behaviour in the spanwise direction is dependant on the scaling of δεwith respect to the critical scaling ε3 2: - if δεε3 2, then there is no roughness effect and perfect slip is achieved in the spanwise direction as well; - if δε≈ε3 2, with δε/ε3 2→λ, 0 < λ < +∞, then partial slip appears in the spanwise direction, which is associated with a constant friction coefficient depending on the 3 limit λand on the riblet geometry; - if δεε3 2, then the same consequences given in [7] hold: no-slip condition appears in the spanwise direction while perfect slip holds in the streamwise one. The common feature of these previous works is that the rugosity effect induced by straight riblets can only appear in the spanwise direction. Thus, our goal in the present paper is to analyse if the periodic displacement of wavy riblets in the lateral direction may produce an effect on the streamwise component of the velocity. For the sake of simplicity, in this work we will consider the same amplitudes and periods in both the streamwise and spanwise directions, by setting a1 ε=a2 ε=δε, p1 ε=p2 ε=ε, where δεεas εtends to 0. The function Ψεwill thus be defined by (see Figure 2) Ψε(x1, x2) = δεΨ1x1 ε+δε εΨ2x2 ε.(1.2) x" 2" x" 1" ε" δ" ε" ε" δ" ε" x" 1" x" 3" Figure 2: Sketch of sinusoidal displacement. View of wavy riblet surface from above (left) and wavy riblet cross-section (right) For such configuration, we prove that the asymptotic effects of wavy riblets also depend on the scaling of δεwith respect to ε; however, two critical regimes need to be considered, namely δε≈ε3 2and δε≈ε5 4, and therefore five different regimes. As it was the case in [15], the critical regime δε≈ε3 2is associated with the effect of the riblets in the spanwise direction. The novelty of the model we consider is that the combination of streamwise and spanwise oscillations of wavy riblets give rise to a new critical regime δε≈ε5 4, that rules the asymptotic effect of the geometry on the streamwise component of the fluid velocity. Thus, two additional asymptotic regimes need to be considered: - if δε≈ε5 4, with δε/ε5 4→µ, 0 <µ<+∞, then no-slip holds in the spanwise direction, while partial slip described by a constant friction coefficient Mµappears in the streamwise direction; - if δεε5 4, then the periodical displacement of wavy riblets in the spanwise direction enforces the fluid to satisfy no-slip condition in both directions. Let us stress that the case δεε5 4, where perfect slip condition transforms asymptotically in no-slip condition, is an example of the strongest rugosity effect: the roughness 4 prevents the fluid from slipping on the wall. This effect was identified in the pioneering paper [13], in the case of a periodic, self-similar boundary with C2-regularity. This regularity was weakened to ”almost Lipschitz” boundaries in [6]. Other results obtained in [9], [10] (see also [23], [26]) give a quite complete understanding of the rugosity effect for arbitrary boundaries. The paper is organized as follows. In Section 2, we introduce the notation. In Section 3, we formulate the problem and state our main result, which is proved in Section 4by means of an adaptation of the unfolding method (see e.g. [3], [11], [14], [16]), that was proposed in [15]. Finally, in Section 5, we illustrate numerically the asymptotic effect of wavy riblets on the streamwise component of the velocity field. To this aim, we simulate the friction coefficient Mµinduced by the critical regime δε≈ε5 4, and the corresponding homogenized system, using different riblet profiles that are widely used in the experimental studies: the V-shape, the U-shape, and the blade riblets. 2 Notation We consider a viscous flow in a domain delimited by two horizontal plates. For the sake of simplicity, we assume that the flow is periodic with respect to the horizontal coordinates (x1, x2). We denote by T2:= R2/Z2(resp. T1:= R/Z) the two-dimensional (resp. one-dimensional) torus. We fix two functions Ψ1,Ψ2∈W2,∞(T1) and define Ψε(x1, x2) := δεΨ1x1 ε+δε εΨ2x2 ε,a.e. (x1, x2)∈ T2,(2.1) where ε > 0 is a small parameter and δεsatisfies lim ε→0δε= 0,lim ε→0 δε ε= 0.(2.2) We consider the spatial domain Ωε:= (x1, x2, x3)∈ T2×R,−Ψε(x1, x2)< x3<1.(2.3) The two components of the boundary ∂Ωεare denoted by Γtop := T2×{1},(2.4) Γε:= (x1, x2, x3)∈ T2×R, x3=−Ψε(x1, x2).(2.5) With these definitions, the upper surface Γtop is assumed to be plane, while the lower one, Γε, is composed of wavy riblets (see Figure 3). We also define the limit domain Ω and its plane bottom surface Γ by Ω := T2×(0,1),Γ := T2×{0}. For k0∈Z2, we denote Ck0 ε=εk0+εT2, Qk0 ε= Ωε∩(Ck0 ε×(−∞,1)) . We define κ:T2→Z2by κ(x0) = k0⇔x0∈Ck0 1. 5 Figure 3: Schematic view of domain Ωε Remark that κis well defined up to a set of zero measure in T2(the set ∪k0∈Z2∂Ck0 1). Moreover, for every ε > 0, we have κx0 ε=k0⇔x0∈Ck0 ε. For a.e. x0∈ T 2we define Cε(x0) as the square Ck0 εsuch that x0belongs to Ck0 ε. For every ε > 0, we define Iεby Iε={k0∈Z2:Ck0 ε∩T26=∅}. We define the set b Q:= T2×(0,+∞), and for every M > 0, we write b QM=T2×(0, M) and we decompose each element x∈b QMin x= (x0, x3), with obvious notation. We denote by Vthe space of functions v:b Q→Rsuch that v∈H1(b QM)∀M > 0,∇v∈L2(b Q)3. Vis a Hilbert space endowed with the norm k·kVdefined by kvk2 V=kvk2 L2(T2×{0})+k∇vk2 L2( b Q). For a bounded measurable set Θ ⊂RN, we denote by L2 0(Θ) the space of functions of L2(Θ) with null integral. We denote by Oεa generic real sequence which tends to zero with εand can change from line to line. We denote by Ca generic positive constant which can change from line to line. 3 Setting of the problem and main result Let Ωεbe the domain defined by (2.3). We consider the Stokes equations   −∆uε+∇pε=fin Ωε, div uε= 0 in Ωε, (3.1) 6 where the unknown functions uεand pεare the velocity and the pressure of the fluid, respectively. Note that we have tacitly assumed that the driving force fis defined on all domains Ωε, say, fis a restriction of a fixed function belonging to the class L∞(T2×R)3. We assume for simplicity that the fluid adheres to the top boundary, that means, uε=0on Γtop .(3.2) The aim is to study the roughness-induced effects of wavy riblets at the bottom of the domain. In view of that, we assume impermeability of the lower surface uε·ν= 0 on Γε,(3.3) together with the perfect slip condition ∂uε ∂ντ = 0 on Γε,(3.4) where νdenotes the outside unitary normal vector to Ωεon Γε, and [ ]τdenotes the tangential component. It is well known (see e.g. [7], [15]) that problem (3.1) has a unique solution (uε, pε)∈ H1(Ωε)3×L2 0(Ωε), and moreover there exists C > 0, which does not depend on ε, such that kuεkH1(Ωε)3+kpεkL2(Ωε)≤C . (3.5) We are now in a position to state the main result of this paper, which describes the asymptotic behaviour of the solution (uε, pε) of problem (3.1)-(3.4). The problem satisfied by the limit (u, p) still possesses the same structure (see formula (3.7)), because the roughness-induced effect appears only on the boundary condition on the limit lower surface Γ. This effect depends on the asymptotic behaviour of δε, with two critical regimes, namely δε≈ε3 2, and δε≈ε5 4, and therefore five different regimes. Theorem 3.1. Let (uε, pε)∈H1(Ωε)3×L2 0(Ωε)be the solution to problem (3.1) such that (3.5) holds. Then, there exists (u, p)∈H1(Ω)3×L2 0(Ω) such that, up to a subsequence, uε*uin H1(Ω)3, pε* p in L2(Ω).(3.6) The pair (u, p)satisfies the Stokes system   −∆u+∇p=fin Ω, div u= 0 in Ω, (3.7) together with the boundary condition u=0on Γtop .(3.8) Moreover, it also satisfies the following boundary condition on Γdepending on the asymptotic behaviour of δε: (i) If δεε3 2, then −∂3u1= 0 ,−∂3u2= 0 , u3= 0 ,on Γ.(3.9) 7 (ii) If δε≈ε3 2, with δε/ε3 2→λ,0<λ<+∞, we define (b φλ,bqλ)∈ V3×L2(b Q), as the solution of                  −∆yb φλ+∇ybqλ=0in R2×R+, divyb φλ= 0 in R2×R+, b φλ 3(y0,0) = λ ∂1Ψ1(y1)on R2×{0}, −∂3b φλ 1= 0,−∂3b φλ 2= 0 on R2×{0}, (3.10) and Mλ∈Rby Mλ=Zb Q|Dyb φλ|2dy . (3.11) Then −∂3u1+Mλu1= 0,−∂3u2= 0, u3= 0 on Γ.(3.12) (iii) If ε3 2δεε5 4, then u1= 0 ,−∂3u2= 0 , u3= 0 ,on Γ.(3.13) (iv) If δε≈ε5 4, with δε/ε5 4→µ,0< µ < +∞, we define (b φµ,bqµ)∈ V3×L2(b Q), as the solution of                −∆yb φµ+∇ybqµ=0in R2×R+, divyb φµ= 0 in R2×R+, b φµ 3(y0,0) = µ2∂1Ψ1(y1)∂2Ψ2(y2)on R2×{0}, −∂3b φµ 1= 0 ,−∂3b φµ 2= 0 ,on R2×{0}, (3.14) and Mµ∈Rby Mµ=Zb Q|Dyb φµ|2dy. (3.15) Then u1= 0 ,−∂3u2+Mµu2= 0 , u3= 0 ,on Γ.(3.16) (v) If δεε5 4, then u1= 0 , u2= 0 , u3= 0 ,on Γ.(3.17) Remark 3.2. Observe that in the first three cases, i.e. when δεε5 4, roughness effects alter only the spanwise component of the velocity, u1, while total slip is satisfied by the streamwise component, u2. However, if δε'ε5 4, the wavy riblets have an influence on both components of the velocity; in the spanwise direction, the no-slip condition u1= 0 is enforced, and - If δε≈ε5 4, partial slip appears in the streamwise direction, associated with a constant friction coefficient Mµ; - If δεε5 4, the no-slip condition u2= 0 is imposed asymptotically. In this case, the limit velocity field is completely adherent to the wall Γ. 8 in T2, we get ZT2×T 2 δε ε3 2 ∂y1Ψ1y1+δε εΨ2(y2)buε,1(x0, y0,0) −δε εZT2 ∂y1Ψ1y1+δε εΨ2(y2)buε,1(x0, z0,0) −¯uε,1(x0) √εdz0 +δ2 ε ε5 2 ∂y1Ψ1y1+δε εΨ2(y2)∂y2Ψ2(y2)buε,2(x0, y0,0) −δ2 ε ε2ZT2 ∂y1Ψ1 y1+δε εΨ2(y2)∂y2Ψ2(y2)buε,2(x0, z0,0) −¯uε,2(x0) √εdz0 +buε,3(x0, y0,0) −¯uε,3(x0) √ε 2 dx0dy0≤Cδε ε→0. (4.25) Depending on the values of λand µ, limits of δε/ε3 2and δε/ε5 4respectively, we deduce: - When δε≈ε3 2, then (4.25) shows that δε ε3 2 ∂y1Ψ1y1+δε εΨ2(y2)buε,1(x0, y0,0) is bounded in L2(T2×T2).(4.26) We can pass to the limit by using (4.23) which implies that ∂y1Ψ1y1+δε εΨ2(y2)buε,1(x0, y0,0) →∂y1Ψ1(y1)u1(x0,0) in L2(T2×T 2). This gives bu3(x0, y0,0) = −λ∂y1Ψ1(y1)u1(x0,0) a.e. (x0, y0)∈ T2×T2. - When δεε3 2, then from (4.26) and (4.23) we deduce ∂y1Ψ1y1+δε εΨ2(y2)buε,1(x0, y0,0) →0 in L2(T2×T 2), which gives u1(x0,0)∂y1Ψ1(y1) = 0 a.e. (x0, y0)∈ T2×T2. - When δε≈ε5 4, then (4.25) shows that δ2 ε ε5 2 ∂y1Ψ1y1+δε εΨ2(y2)∂y2Ψ2(y2)buε,2(x0, y0,0) is bounded in L2(T2×T2). (4.27) We can pass to the limit by using (4.23) which implies that ∂y1Ψ1y1+δε εΨ2(y2)∂y2Ψ2(y2)buε,2(x0, y0,0) →∂y1Ψ1(y1)∂y2Ψ2(y2)u2(x0,0) in L2(T2×T 2). This gives bu3(x0, y0,0) = −µ2∂y1Ψ1(y1)∂y2Ψ2(y2)u2(x0,0) a.e. (x0, y0)∈ T 2×T2. 15 - When δεε5 4, then from (4.27) and (4.23) we deduce ∂y1Ψ1y1+δε εΨ2(y2)∂y2Ψ2(y2)buε,2(x0, y0,0) →0 in L2(T2×T 2), which gives u2(x0,0)∂1Ψ1(y1)∂2Ψ2(y2) = 0 a.e. (x0, y0)∈ T2×T2. Finally, using these properties, we easily deduce the statements of Lemma 4.4. Proof of Theorem 3.1.Thanks to (3.5), there exists a subsequence of ε, still denoted by ε, and (u, p)∈H1(Ω)3×L2(Ω) such that (3.6) holds. On the other hand, we observe that (uε, pε) satisfies the variational equation     ZΩε Duε:Dvε−ZΩε pεdiv vεdx =ZΩε f·vεdx , ∀vε∈H1(Ωε)3,vε·ν= 0 on Γε,vε=0on Γtop . (4.28) The proof of Theorem 3.1 is similar to the one given in [15], and it will be carried out using suitable test functions vεdepending on the values of δε. In particular, the cases (i)−(iii) have already been addressed in [15]. Here, we will develop in more detail the new critical case (iv), that is when δε≈ε5 4, with δε/ε5 4→µ, 0 <µ<+∞. Finally, case (v) will be described briefly. Step 1 Case δε≈ε5 4, with δε/ε5 4→µ, 0 <µ<+∞. We consider v∈C1 c(T2×(−1,1))3and b v∈C1 c(T2;C(b Q)3), with Dyb v(x0, y) = 0 a.e. in {y3> M}, for some M > 0, such that          v(x) = v(x0,0) if x3≤0, v3(x0,0) = 0 , v1(x0,0)∂y1Ψ1(y1) = 0 a.e. (x0, y0)∈ T 2×T2,   b v(x0, y) = b v(x0, y0,0) if y3≤0, bv3(x0, y0,0) = −µ2∂y1Ψ1(y1)∂2Ψ2(y2)v2(x0,0) . (4.29) Besides, we take ζ∈C∞(R) such that ζ(x3) =      1 if x3<1 3, 0 if x3>2 3, (4.30) and Rε>0 such that Rε→+∞, Rε"δ2 ε ε5 2−µ22 +δ2 ε ε2#→0.(4.31) 16 Then, we define vε∈H1(Ωε)3by                                                                    vε,1(x) = v1(x) + √εbv1x0,x εζ(x3) vε,2(x) = v2(x) + √εbv2x0,x εζ(x3) vε,3(x) = v3(x) + √ε"bv3x0,x εζ(x3)+ +µ2ζx3 εRε∂1Ψ1x1 ε∂2Ψ2x2 εv2(x0,0) −δ2 ε ε5 2 ζx3 εRε∂1Ψ1x1 ε+δε εΨ2x2 ε∂2Ψ2x2 εv2(x0,0) −δ2 ε ε2ζx3 εRε∂1Ψ1x1 ε+δε εΨ2x2 ε∂2Ψ2x2 εbv2x0,x0 ε,0# −δε εζx3 ε2Rε∂1Ψ1x1 ε+δε εΨ2x2 εv1(x0,0) −δε √εζx3 εRε∂1Ψ1x1 ε+δε εΨ2x2 εbv1x0,x0 ε,0. (4.32) Taking into account properties given in (4.29), the sequence vεsatisfies vε= 0 on Γtop,vε·ν= 0 on Γε. Thus, we can take such vεas test function in (4.28). The problem is to pass to the limit in every term that appears in (4.28). To simplify the computations, we will first estimate the derivatives of vε. Taking into account that Dyb v= 0 a.e. in {y3> M}and that ζ= 1 a.e. on {x3<1 3}, we have Dvε(x) = Dv(x) + 1 √εDyb vx0,x ε+hε(x),(4.33) where, using that v,b vand ζare bounded with bounded derivatives, the function hε∈C0(Ωε)3×3satisfies |hε| ≤ C√ε+C"1 √εRε +1 √ε δ2 ε ε5 2−µ2+δ2 ε ε2 +δ2 ε ε5 2 δε ε3 2 +δ2 ε ε2+δ2 ε ε2 δε ε3 2 +δε ε3 2#χ{x3<εRε} +C"δε ε 1 ε2Rε +δε ε 1 ε+δε ε δε ε2+δε ε#χ{x3<ε2Rε} +C"δε ε3 2 +δ2 ε ε5 2 +δε ε1 2 +δε ε1 2 1 εRε#χ{x3<εRε}, a.e. in Ωε. Using that Rεtends to infinity and that δε/ε5 4is bounded, we get 17 |hε| ≤ C√ε+C"1 √ε δ2 ε ε5 2−µ2+δ2 ε ε2+δε ε3 2 +δ2 ε ε2+δ2 ε ε2 δε ε3 2#χ{x3<εRε} +C"δε ε 1 ε+δε ε δε ε2+δε ε#χ{x3<ε2Rε} +C"δε ε3 2 +δ2 ε ε5 2 +δε ε1 2#χ{x3<εRε}. Therefore, by (4.31) we have ZΩε|hε|2dx ≤Oε+CRε"δ2 ε ε5 2−µ22 +δ2 ε ε2+δ4 ε ε3+δ6 ε ε6+δ4 ε ε4+δ2 ε# ≤Oε+CRε"δ2 ε ε5 2−µ22 +δ2 ε ε2#=Oε. (4.34) Taking vεas a test function in (4.28) and using that kuεkH1(Ωε)3,kpεkL2(Ωε)are bounded, kvε−vkC0(¯ Ωε)3tends to zero, (4.33) and (4.34), we get ZΩε (Duε:Dv−pεdiv v)dx +1 √εZΩεDuε:Dyb vx0,x ε−pεdivyb vx0,x εdx +Oε=ZΩε f·vdx . (4.35) In last equality, we have used the estimate ZΩε\ΩDv+1 √εDyb vx0,x ε 2 dx ≤C ε|Ωε\Ω| ≤ Cδε ε≤C√ε . Therefore, (4.35) can be rewritten as ZΩ (Duε:Dv−pεdiv v)dx +1 √εZΩDuε:Dyb vx0,x ε−pεdivyb vx0,x εdx +Oε= 0 . Taking into account (3.6), we deduce ZΩ (Du:Dv−pdiv v)dx +1 √εZΩDuε:Dyb vx0,x ε−pεdivyb vx0,x εdx +Oε=ZΩ f·vdx , (4.36) with vand b vsatisfying (4.29). In order to estimate the second term in (4.36), we introduce the sequences b uε,bpεrespectively defined by (4.1) and (4.2). By (3.5) 18 and Lemmas 4.2 and 4.4, we can assume that there exist bp∈L2(T2×b Q) and b u∈L2(T2;V3) which satisfy (4.5), (4.8), (4.9) and (4.11). This yields: 1 √εZΩDuε:Dyb vx0,x ε−pεdivyb vx0,x εdx =X k0∈IεZQk0 ε1 √εDuε(x) : Dyb vx0,x ε−1 √εpε(x) divyb vx0,x εdx =X k0∈Iε ε3ZQk0 ε1 √εDuεε(k0+y0), εy3:Dyb vx0, y−1 √εpεε(k0+y0), εy3divyb vx0, ydy =X k0∈Iε ε2ZQk0 εε √εDuεε(k0+y0), εy3:Dyb vx0, y−√εpεε(k0+y0), εy3divyb vx0, ydy =X k0∈Iε ε2ZQk0 ε1 √εDyb uε(x0, y) : Dyb vx0, y−√εbpε(x0, y)divyb vx0, ydy =ZT2Zb QM1 √εDyb uε(x0, y) : Dyb v(x0, y)−√εbpε(x0, y)divyb v(x0, y)dx0dy =ZT2× b QM (Dyb u:Dyb v−bpdivyb v)dx0dy +Oε. Substituting in (4.36) we set ZΩ (Du:Dv−pdiv v)dx +ZT2× b QM (Dyb u:Dyb v−bpdivyb v)dx0dy =ZΩ f·vdx, (4.37) for every v∈C1 c(T2×(−1,1))3,b v∈C1 c(T2;C1(b Q)3), with Dyb v(x0, y) = 0 a.e. in {y3> M}, for some M > 0, and such that (4.29) is satisfied. By density, this equality holds true for every v∈H1(Ωε)3, and every b v∈L2(T2;V3) such that v= 0 on Γtop , v3(x0,0) = 0,bv3(x0, y0,0) = −µ2∂y1Ψ1(y1)∂y2Ψ2(y2)v2(x0,0) , v1(x0,0)∂y1Ψ1(y1)=0,a.e. (x0, y0)∈ T 2×T2. Let us now obtain an equation for ueliminating b uand bpin (4.37). For this purpose, we take v=0in (4.37). This proves that (b u,bp) is a solution of                        −∆yb u+∇ybp=0in R2×R+, divyb u= 0 in R2×R+, (b u,bp)∈ V3×L2(b Q), bu3(x0, y0,0) = −µ2∂y1Ψ1(y1)∂y2Ψ2(y2)u2(x0,0) on R2×{0}, −∂y3bu1= 0,−∂y3bu2= 0 on R2×{0}, (4.38) a.e. in T2. Defining (b φµ,bqµ) by (3.14), we deduce by linearity and uniqueness Dyb u(x0, y) = −u2(x0,0)Dyb φµ(y) a.e. in R2×(0,+∞),(4.39) 19 bp(x0, y) = u2(x0,0)bqµ(y) a. e. in R2×(0,+∞).(4.40) Now, consider v∈H1(Ω)3, with v= 0 on Γtop,v3= 0 on Γ. We take vand b v(x0, y) = −v2(x0,0)b φµ(y) as test functions in (4.37); taking into account (4.39), we get ZΩ (Du:Dv−pdiv v)dx +ZΓ Mµu2v2dx0=ZΩ f·vdx . (4.41) By the arbitrariness of v, this proves that (u, p) is a solution of problem (3.7)-(3.8)- (3.16). Step 2 The case δεε5 4. As in Step 1, we consider v∈C1 c(T2×(−1,1))3, with v(x0, x3) = v(x0,0) if x3≤0, v3= 0 on Γ. Then, for ζ∈C∞(R), which satisfies (4.30), and Rεsatisfying (4.31), we define vε∈H1(T2×(−1,1))3by                        vε,1(x) = v1(x) vε,2(x) = v2(x) vε,3=v3(x)−δε εζx3 ε2Rε∂1Ψ1x1 ε+δε εΨ2x2 εv1(x) −δ2 ε ε2ζx3 ε2Rε∂1Ψ1x1 ε+δε εΨ2x2 ε∂2Ψ2x2 εv2(x). The sequence vεsatisfies vε·ν= 0 on Γ and vε= 0 on Γtop. We can take vεas test function in (4.28). Proceeding as in Step 1, we can prove that Dvε=Dv+hε, with hεsatisfying ZΩε|hε|2dx ≤Oε. Then, passing to the limit in εand reasoning by density we get            ZΩ Du:Dvdx −ZΩ pdiv vdx = 0 ∀v∈H1(Ω)3, v3= 0 on Γ ,v=0on Γtop , v1(x0,0)∂1Ψ1(y1) = 0, v2(x0,0)∂1Ψ1(y1)∂2Ψ2(y2) = 0 a.e. (x0, y0)∈ T 2×T2. This is equivalent to problem (3.7)-(3.8)-(3.17). Proof of Proposition 3.4.It results from Lemma 4 in [2], applied to (b φλ−Vλ,bpλ), where Vλis a constant vector defined by Vλ i=ZT2b φλ i(y0,0) dy0, i = 1,2, V λ 3= 0 . 20 This translation ensures that the vector field b φλ−Vλsatisfies the following property: ZT2 (b φλ(y0,0) −Vλ)dy0= 0 . Notice that, by periodicity, the boundary condition satisfied by b φλ 3on y3= 0 yields RT2b φλ 3(y0,0) dy0= 0 . 5 Numerical results In this section, we simulate numerically the asymptotic effect of wavy riblets on a viscous fluid, in the streamwise direction. To this aim, we consider a flow driven by a constant external force f= (0,1,0) and modeled by Stokes equations (3.1) associated with boundary conditions (3.2), (3.3), (3.4), and periodic in the lateral coordinates. In order to capture the influence of fast oscillations of the ribbed boundary, we simulate the homogenized system described by Theorem 3.1. Since we focus on the streamwise component of the fluid velocity, we consider the second critical case, namely, δε∼µ ε5/4 with 0 < µ < +∞. In that case, the homogenized velocity field uis the solution to Stokes equations (3.7) in the unit cube Ω, completed with boundary conditions (3.8), (3.16). Boundary condition (3.16), satisfied by the streamwise component of the velocity, introduces an extra friction term Mµdefined by (3.15). As a result, the numerical approximation of uis determined by the resolution of the auxiliary problem (3.14), that allows for the computation of this friction coefficient. This section is organized as follows. In Subsection 5.1, we present and justify the numerical approximation of the solution to the auxiliary problem (3.14). In Subsection 5.2, we define the geometry of the wavy riblets that we consider in the simulations. Finally, Subsection 5.3 is devoted to the statement and discussion of the results. 5.1 Numerical approximation of auxiliary system (3.14) In order to simulate system (3.14), we fix a truncature parameter H > 0 and we approximate the functions (b φµ,bqµ) on T2×(0, H), as follows. We define (b φµ H,bqµ H)∈ H1(R2×(0, H),R3)×L2(R2×(0, H)), as the solution of                                      −∆yb φµ H+∇ybqµ H=0in R2×(0, H), divyb φµ H= 0 in R2×(0, H), b φµ H,3(y0,0) = µ2∂1Ψ1(y1)∂2Ψ2(y2) on R2×{0}, −∂3b φµ H,1= 0,−∂3b φµ H,2= 0 on R2×{0}, ∂3b φµ H,1(y0, H)=0, ∂3b φµ H,2(y0, H)=0, ∂3b φµ H,3−bqµ H= 0 on R2×{H}, (5.1) and the corresponding approximate value of the coefficient Mµ,H by Mµ,H =ZT2×(0,H)Dyb φµ H2dy . (5.2) 21 Justification of the approximation. The validity of the approximate system (5.1) is strongly related to the exponential decay properties (3.22). Assume, for instance, that H > 1. We denote by FH∈C∞(T2)3the function defined by FH(y0) := ∂3b φλ(y0, H)−bqλ(y0, H)e3∀(y1, y2)∈ T 2. Applying property (3.22) with α= 1, t = 1, we obtain the existence of a constant C > 0 such that kFHkL2(T2)3≤C µ2k∂1Ψ1kL2(T1)k∂2Ψ2kL2(T1)e−τH. In order to compare the exact solution (b φµ,bqµ) with the approximation (b φµ H,bqµ H), we notice that (b φµ,bqµ) is a solution to system (5.1) in T2×(0, H), but with a slightly modified boundary condition on y3=H, namely, ∂3b φµ−bqµe3=FHon y3=H . As a result, the difference (b φµ−b φµ H,bqµ−bqµ H) is a solution to a Stokes system in T2×(0, H), with the following boundary conditions: −∂3(b φµ 1−b φµ H,1)=0,−∂3(b φµ 2−b φµ H,2) = 0 on y3= 0 , b φµ 3−b φµ H,3= 0 on y3= 0 , ∂3(b φµ−b φµ H)−(bqµ−bqµ H)e3=FHon y3=H . It follows that there exists a constant C > 0 such that kD(b φµ−b φµ H)kL2(T2×(0,H)) ≤C µ2k∂1Ψ1kL2(T1)k∂2Ψ2kL2(T1)e−τH .(5.3) Numerical method. In order to compute the solution to system (5.1), we use a finite element discretization of the problem, using FreeFem++ software [21]. 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 [20], [5]. The three-dimensional mesh of the fluid domain T2×(0, H) is obtained by a constrained Delaunay tetrahedralization. Finally, the incompressibility condition is approximated by the standard penalty formulation ZΩdivb φµ+ηbqµbr= 0 ∀br∈Qh, where Qhis the discrete approximation space for the pressure, and η= 10−6is the penalty parameter [20], [22]. 5.2 Definition of the geometry of the wavy riblets The geometry of the wavy riblets is determined by the choice of the profile functions Ψ1,Ψ2. To perform comparisons, and in view of estimate (5.3), we normalize functions Ψ1,Ψ2by setting k∂1Ψ1kL2(T1)= 1 ,k∂2Ψ2kL2(T1)= 1 . Choice of Ψ2.We consider a sinusoidal function Ψ2(y2) = K2sin(2πy2), where K2is the normalization constant. 22 Choice of Ψ1.The functions Ψ1that we consider are obtained by regularizing riblet profiles that are widely used in the experimental studies [18]: (i) the V-shape riblets, or sawtooth riblets; (ii) the U-shape riblets, or scallop riblets; (iii) the blade riblets. For each of these profiles, we have applied a parabolic regularization in the vicinity of each singular point, obtaining functions Ψ1of class C2(see Figure 4). Finally, Figure 5 provides an outline of the associated profiles Ψε, restricted to the line R×{0}. Figure 4: Functions Ψ1associated with V-shaped, U-shaped and blade riblets (from left to right) Figure 5: Functions Ψεassociated with V-shaped, U-shaped and blade riblets (from left to right) 5.3 Results and discussion For each type of geometry of riblets described in Subsection 5.2, we compute an approximate value of the friction coefficient Mµ, for different values of µ. To this aim, we start by solving system (5.1) in a domain of height H= 3, for a fixed value µ= 1. The corresponding results are given in Table 1. Then, for any value of µ, we set Mµ,H =µ4M1,H , which is the numercial counterpart of the relation Mµ=µ4M1, satisfied by the exact friction coefficient as a direct consequence of (3.15). Then, for each value of Mµ,H, we solve by a finite element method (see Subsection 5.1 for details) the homogenized Stokes system 23 V-shape U-shape Blade Mµ,H 311.516 1270.89 1865.94 Table 1: Approximate values of the friction coefficient Mµ,H, with µ= 1, and H= 3, for the U-shape, V-shape and blade riblets.                −∆u+∇p=fin Ω , div u= 0 in Ω , u= 0 on Γtop , u1= 0 ,−∂3u2+Mµ,H u2= 0 , u3= 0 on Γ , (5.4) where Ω = T2×(0,1), Γ = T2×{0}, Γtop =T2×{1}, in the particular case f= (0,1,0). We define the average streamwise velocity u2by u2=ZΩ u2=ZΩ f·u. The results are plotted in Figure 6. We observe that for each riblet profile, there exists a transition between the perfect slip regime, associated with the limit µ→0, and the no-slip regime µ→+∞. During this transition, the values of the average streamwise velocities that we observe for the V-shape, the U-shape and the blade riblets, are ordered in a way that is consistent with the values of the approximate friction coefficients given in Table 1. The V-shape riblet, associated with the smallest friction coefficient, corresponds to the highest velocity of the flow; on the opposite, the U-shape and the blade riblets, whose friction coefficients at µ= 1 are 4 to 6 times larger, have a stronger slowing effect on the flow. This significant difference between the V-shape riblet on the one hand, and the U- shape and blade riblets on the other hand, can be put in perspective with the well-known experimental results on straight riblets, that demonstrate the efficiency of thin protrusions in the reduction of the drag for turbulent flows (see [25] and the references therein). In our study, such fast and well-localized variations of the boundary give rise to a high L1 norm of the derivative ∂1Ψ1; combined with sinusoidal oscillations of the boundary along the streamwise direction, they generate high friction coefficients Mµ. Such high values of the friction coefficients could play a decisive role in the stabilization of the flow in the vicinity of the wall, in the case of flows at larger Reynolds numbers. This question will be addressed in a future work, by applying our homogenization results to certain turbulent models used in Large Eddy Simulation. Acknowledgements. The research of Francisco J. Su´arez-Grau has been partially supported by the project MTM2011-24457 of the “Ministerio de Econom´ıa y Competitividad”, by the project FQM309 of the “Junta de Andaluc´ıa”, and by the program “Research in Paris 2014” of Mairie de Paris . The research of Giordano Tierra has been partially supported by ERC-CZ project LL1202 (Ministry of Education, Youth and Sports of the Czech Republic). 24