Subhandu Rawat

Agathe Chouippe

Rémi Zamansky

Dominique Legendre
email: dominique.legendre@imft.fr

Eric Climent

Drag modulation in turbulent boundary layers subject to different bubble injection strategies

Keywords:

The aim of this study is to investigate numerically the interaction between a dispersed phase composed of micro-bubbles and a turbulent boundary layer flow. We use the Euler-Lagrange approach based on Direct Numerical Simulation of the continuous phase flow equations and a Lagrangian tracking for the dispersed phase. The Synthetic Eddy Method (SEM) is used to generate the inlet boundary condition for the simulation of the turbulent boundary layer. Each bubble trajectory is calculated by integrating the force balance equation accounting for buoyancy, drag, added-mass, pressure gradient, and the lift forces. The numerical method accounts for the feedback effect of the dispersed bubbles on the carrying flow. Our approach is based on local volume average of the two-phase Navier-Stokes equations. Local and temporal variations of the bubble concentration and momentum source terms are accounted for in mass and momentum balance equations. To study the mechanisms implied in the modulation of the turbulent wall structures by the dispersed phase, we first consider simulations of the minimal flow unit laden with bubbles. We observe that the bubble effect in both mass and momentum equations plays a leading role in the modification of the flow structures in the near wall layer, which in return generates a significant increase of bubble volume fraction near the wall. Based on these findings, we discussed the influence of bubble injection methods on the modulation of the wall shear stress of a turbulent boundary layer on a flat plate. Even for a relatively small bubble volume fraction injected in the near wall region, we observed a modulation in the flow dynamics as well as a reduction of the skin friction.

Introduction

It is a well-established fact that any object having nonzero relative velocity with respect to the upstream fluids experiences viscous drag at large Reynolds number due to the thin boundary layer developing on the object's surface, where the velocity of the fluid gradually reduces to that of the object surface. The reduction of viscous drag is one of the most active areas of research in fluid mechanics by virtue of its importance in numerous technological applications, such as ship locomotion, aircraft flight, fuel transportation through pipelines and vehicle aerodynamics. Plethora of experimental and numerical studies have been conducted in order to manipulate shear induced drag. In the present study we have used direct numerical simulations to study the drag modification caused by the interaction of micro-bubbles with the near-wall structures. Ever since the first DNS was conducted by [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF] , it has become a fundamental tool to a better understanding of turbulence. With the rapid advancement in data storage capability and processor speed researchers can do DNS with a large number of bubbles. Whether the researches use experiments or DNS, most of the drag reduction techniques focus on the manipulation of near-wall cycle of turbulence regeneration using either active or passive mechanisms. The near-wall structures enhance the mixing and thus the transfer of momentum towards the wall, hence are directly responsible for the production of drag.

Near-wall regeneration cycle

The flow visualization experiments performed by [START_REF] Kline | The structure of turbulent boundary layers[END_REF] revealed structures of near-wall turbulence called streaks. Streaks are quasistreamwise regions where the velocity is less than the local mean velocity at the same height from the wall. The average spanwise streak spacing in the viscous sublayer and the buffer layer scales in wall units and corresponds approximately to λ + z = 100 , in wall units. These coherent structures are ubiquitous in all turbulent shear flows. Numerous experimental and numerical investigations observed that these structures are self-sustained through a cycle of breaking up and regeneration. The flow structures are randomly located in space and time, making very difficult to extract a simple mechanism to explain the breakdown and regeneration of them. With the minimal-flow strategy of [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF] this randomness is reduced and it becomes possible to identify a simple regeneration cycle. Following this approach [START_REF] Hamilton | Regeneration mechanisms of near-wall turbulence structures[END_REF] investigated the minimum Reynolds number and domain size needed to sustain a turbulent regime. The flow visualization and careful analysis of the energy associated with various Fourier modes evidenced a three-step process. First, the interaction of the spanwise modulation of the streamwise velocity component (i.e., near wall streaks) with the mean shear, leads to the lift-up effect. This lift-up effect is responsible for the redistribution of streamwise momentum through counter-rotating quasistreamwise vortices. This redistribution of momentum leads to the transient growth of the high speed and low speed streamwise streaks. Second, when these streaks reach high amplitude they become unstable to secondary perturbation via inflectional type instability. Finally, the unstable modes of the streaks selfinteract to re-energize the streamwise vortices. This phenomenon of lift up and breakdown is called bursting. [START_REF] Kline | The structure of turbulent boundary layers[END_REF] reported that a strong favorable pressure gradient suppresses the bursting process and flow return to laminar state. This gave striking evidence that bursting plays a crucial role in turbulence generation near the wall.

Mechanism of drag reduction

Bursting of near-wall coherent structures intermittently produce and eject turbulent fluctuations into the core region of the domain (see [START_REF] Itano | The dynamics of bursting process in wall turbulence[END_REF] ). In between the two bursts there is a quiescent state where the streaks get produced, become unstable and breakdown. Numerous attempts have been made to reduce the viscous drag under the premise that the bursting cycle can be controlled by the modification and alteration of near-wall coherent structures. It is observed that streamwise riblets reduce the drag by stabilizing the coherent streaks (see [START_REF] Chu | A direct numerical simulation of laminar and turbulent flow over riblet-mounted surfaces[END_REF] ). [START_REF] Sirovich | Turbulent drag reduction by passive mechanism[END_REF] used specific patterns of protrusions dispersed on the wall, which interfere with the coherent structures such that it reduces the rate at which energy is dissipated. They found that the streak spacing increased by 10% and fluctuations of wall-normal velocity component is diminished as compared to the smooth wall which was consistent with drag reduction using polymer additives close to the wall [START_REF] Yiqing | Suppressing wall turbulence by means of a transverse traveling wave[END_REF] . used transverse traveling waves in their direct numerical simulations to suppress wall turbulence. They used finely tuned transverse force confined to the viscous sublayer to induce a traveling wave which consequently eliminates the near-wall streaks [START_REF] Fransson | Experimental study of the stabilization of Tollmien-Schlichting waves by finite amplitude streaks[END_REF] . induced streaks by generating optimal vortices using equally spaced cylindrical roughness elements just before the leading edge of the boundary layer. It was shown that the streaks suppress the instability of Tollmien-Schlichting (TS) waves thereby delaying the transition to turbulence and reducing viscous drag.

Injections of micro-bubbles in a turbulent flow is another robust method for drag reduction. [START_REF] Mccormick | Drag reduction of a submersible hull by electrolysis[END_REF] were the first to demonstrate this phenomenon for fully submerged asymmetric body, where hydrogen bubbles generated via electrolysis in the boundary layer resulted in significant drag reduction. Since then numerous approaches both experimental [START_REF] Gabillet | Experimental study of bubble injection in a turbulent boundary layer[END_REF][START_REF] Guin | Reduction of skin friction by microbubbles and its relation with near-wall bubble concentration in a channel[END_REF][START_REF] Jacob | Drag reduction by microbubbles in a turbulent boundary layer[END_REF][START_REF] Madavan | Reduction of turbulent skin friction by microbubbles[END_REF][START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF][START_REF] Park | Drag reduction promoted by repetitive bubble injection in turbulent channel flows[END_REF] and numerical [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF][START_REF] Ferrante | Reynolds number effect on drag reduction in a microbubble-laden spatially developing turbulent boundary layer[END_REF][START_REF] Kanai | Direct numerical simulation of wall turbulent flows with microbubbles[END_REF][START_REF] Mattson | Simulation of bubble migration in a turbulent boundary layer[END_REF][START_REF] Pang | Numerical study on modulation of microbubbles on turbulence frictional drag in a horizontal channel[END_REF][START_REF] Spandan | Drag reduction in numerical two-phase Taylor-Couette turbulence using an Euler-Lagrange approach[END_REF][START_REF] Xu | Numerical simulation of turbulent drag reduction using micro-bubbles[END_REF][START_REF] Xu | Numerical simulation of turbulent drag reduction using micro-bubbles[END_REF] have been carried out in order to understand this phenomenon and the relative importance of other physical parameters such as diameter and concentration of the bubble, bubble deformation, bubble coalescence and breakdown and the method of bubble injection into the flow, and promising efficiency in actual ships can now be obtained [START_REF] Kumagai | Power-saving device for air bubble generation using a hydrofoil to reduce ship drag: theory, experiments, and application to ships[END_REF] . [START_REF] Xu | Numerical simulation of turbulent drag reduction using micro-bubbles[END_REF] numerically simulated a fully developed turbulent channel flow at Re τ = 135 and injected bubbles of different sizes to achieve drag reduction. They proposed three possible mechanisms: first associated with initial seeding of the bubbles, second is related to the density variation which reduces the turbulence momentum transfer and third is regulated by the correlation between the bubbles and turbulence. They emphasized that unlike the drag reduction techniques of riblets, spanwise traveling waves, "V" shaped protrusions or polymers where the spanwise spacing was increased and streaks were stabilized, here the streak spacing seems to be affected very little and they become more distorted due to the unsteady bubble forcing. [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] numerically simulated turbulence boundary layer seeded with micro-bubbles, they took the concentration of bubbles into account while solving the continuity and momentum equations for the carrier fluid flow. It was reported that micro-bubbles create positive divergence for the fluid and shifts the quasi-streamwise vortical structures away from the wall thereby displacing Reynolds stress u v peak in positive wall normal direction which results in the reduction of production of turbulent kinetic energy. The spanwise gap between the streaks close to the wall was found to be increased, this effect was similar to the drag reduction method of riblets, traveling waves or polymer but different from findings of [START_REF] Xu | Numerical simulation of turbulent drag reduction using micro-bubbles[END_REF] . Furthermore, [START_REF] Jacob | Drag reduction by microbubbles in a turbulent boundary layer[END_REF] conducted experimental investigation of the turbulent boundary layer over a flat plate with "plate on top" configuration injecting bubbles at the leading edge. They found that bubbles get preferentially accumulated in the narrow spatial region close to the wall and the bubble concentration peaks at y + = 25 due to buoyancy force that is acting towards the wall for this configuration. It was asserted this high bubble concentration close to the wall decreases the coherence of the near-wall structures thereby reducing their length scales. This eventually leads to the reduction of momentum flux towards the wall which resulted in the decrease in wall shear stress by 25%.

The case of the vertical channel flow has been investigated experimentally [START_REF] Serizawa | Pseudo-laminar characteristic of bubbly flow with micro bubbles and practical applications of micro bubble[END_REF] and numerically [START_REF] Molin | Turbulence modulation and microbubble dynamics in vertical channel flow[END_REF] , and an increase (respectively decrease) of the wall shear stress has been observed in upflow (respectively downflow) configurations. As described by [START_REF] Molin | Turbulence modulation and microbubble dynamics in vertical channel flow[END_REF] the presence of the dispersed phase will induce, through the forcing term in the momentum equation, an equivalent additional pressure gradient equal to ±φ v (ρ f -ρ b ) g (where φ v is the global void fraction, ρ f and ρ b the density of the fluid and the bubbles and g the gravitational acceleration) that contributes to a significant amount of wall shear stress modification.

Goal of the present study

As described in the section §1.2 significant amount of effort has been put in order to advance the drag reduction techniques in the last three decades. However, it is hard to ignore that there is still no consensus on the mechanism behind drag reduction observed due to bubble injection. For instance, some studies suggest the decrease in the coherence of the near-wall structures due to the interaction of bubbles with near-wall turbulence, whereas others reported near-wall coherent structures becoming more organized as a consequence of streaks stabilization. The mechanism responsible for the drag reduction is dependent on the void fraction, bubble diameter, physical forces like buoyancy and lift forces and initial seeding of bubbles. The goal of our study is to investigate numerically the dispersion of bubbles in turbulent boundary layer flows and the effects of different bubble injection methods on the drag modulation. For this purpose we first consider the minimal-flow unit configuration in order to isolate the influence of the different contributions of the two-way coupling, namely the forcing due to momentum exchange between the bubbles and the flow, the influence of the local concentration in the continuity equation (i.e. divergence effect) and in the momentum equation. We propose four configurations in which those contributions are activated separately. Then the Synthetic Eddy Method (SEM) is used to gener-ate the inlet boundary condition for the simulation of the turbulent boundary layers. A numerical strategy based on two different domains is proposed to save computational time. Different bubble injection methods for a fixed bulk void fraction were tested to know whether they have different im pact on the flow modification and their influence on the near-wall regeneration cycle. Nondeformable spherical bubbles which are typically smaller than the smallest structure of the background flow are considered in this study. Lagrangian tracking of the dispersed phase fits this type of physical configuration. Our approach is in this way similar to the Lagrangian tracking employed by [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF] in Taylor-Couette flow and can reproduce effective compressibility effects due to void fraction variations as described in the work of [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] .

The paper is organized as follows: Section 2 is devoted to a description of the numerical method. We introduce the numerical strategy based on the Synthetic Eddy Method (SEM) to simulate turbulent boundary layer and we detail the model accounting for the feedback effect of the bubbles on the carrying phase. In Section 3 some direct numerical simulations of minimal flow unit configuration were conducted to study the preferential accumulation of bubbles and investigate their influence on the turbulent structures and regeneration cycle. Finally, the main results on the drag modulation related to the injection of bubbles in a boundary layer are reported in Section 4 . More specific details on the SEM is provided in Appendix B followed by a validation of the simulation of the turbulent boundary layer as well as a set of verifications of the implementation of the different contributions of the two-way coupling of the bubble phase on the liquid.

Numerical method

The numerical approach is based on Euler-Lagrange modeling: the continuous phase flow is predicted through direct solution of the Navier-Stokes equations while Lagrangian bubble trajectories are computed by numerical integration of momentum balance equation. Modeling the presence of bubbles is based on volumeaveraging of the momentum and continuity equations including momentum source terms together with spatial and temporal void fraction variations.

Equation for the carrying single phase flow

Assuming constant physical properties (viscosity and density) for a Newtonian incompressible fluid, we consider the Navier-Stokes equations for the continuous (liquid) phase

∇ • u = 0 (1) ∂u ∂t + ∇ • (uu ) = - 1 ρ ∇ p + ∇ • ν ∇u + ( ∇u ) T + g (2)
The system of equations is discretized on a staggered nonuniform grid with a finite volume approach. Spatial derivatives are calculated with second order accuracy and we use semi-implicit Crank-Nicolson scheme for the viscous term and three steps Runge-Kutta scheme for time integration. The JADIM code has already been widely used and validated for laminar and turbulent configurations in single-phase configurations (see for example during the last two decades [START_REF] Calmet | Statistical structure of high-reynolds-number turbulence close to the free surface of an open-channel flow[END_REF][START_REF] Hallez | Interaction between two spherical bubbles rising in a viscous liquid[END_REF][START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF][START_REF] Magnaudet | Accelerated flows past a rigid sphere or a spherical bubble. part 1. steady straining flow[END_REF][START_REF] Merle | Forces on a high-Reynolds-number spherical bubble in a turbulent flow[END_REF] ). Spatially developing flows as a turbulent boundary layer pose some challenges: if one wants to avoid an unnecessary long domain, time-dependent inlet boundary conditions need to be considered. The flow downstream is highly dependent on the inlet condition which should be very close to the real turbulence. Turbulent inlet generation techniques for spatially developing flows have been an active area of research from the last three decades.

For bubble laden turbulent boundary layer, [START_REF] Ferrante | A robust method for generating inflow conditions for direct simulations of spatially-developing turbulent boundary layers[END_REF] adapted a spectral method to synthetically generate turbulent initial condition with prescribed energy spectra and Reynolds stress tensor. With this initialization they did a small precursor simulation and combined it with the method of [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF] to generate the inlet velocity field for the main simulation with longer domain. In the present study we implemented Synthetic Eddy Method ( SEM ) proposed by [20] to generate the inlet boundary condition. SEM gives a time dependent two dimensional velocity field at the inlet plane with the prescribed mean and covariance along with two points and two times correlations. This approach along with the numerical implantation details are presented in Appendix A . Note that the recycling technics ( [START_REF] Spalart | Direct numerical simulation of equilibrium turbulent boundary layers[END_REF] , [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Re θ = 1410[END_REF] , [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF] ) are not considered here due to the two-way coupling with the bubbly phase. Indeed, the flow is already modified and therefore cannot be recycled as the flow is not pseudoperiodic anymore.

Lagrangian tracking of bubbles

We utilize two-way coupled Euler-Lagrangian framework where the dispersed phase is composed of N b pointwise spherical bubbles. Bubble trajectories are computed solving Newton's second law, where the rate of change of momentum of the bubble is equal to the total external forces F b acting on it. The external forces consist of drag, fluid acceleration, added-mass and lift forces. The radius, volume and velocity of the bubble are noted as R b , V b = (4 / 3) π R 3 b and v , respectively. The vorticity of the liquid is . The bubble force balance can be written as:

d x d t = v (3)
ρ b V b d v d t = (ρ b -ρ f ) V b g -ρ f V b 3 8 R b C D | v -u | (v -u ) + ρ f V b D u D t + ρ f V b C M D u D t - d v d t -ρ f V b C L (v -u ) × (4)
where C D is the drag coefficient, C M the added mass coefficient, C L is the lift coefficient and the history force is neglected [START_REF] Magnaudet | The motion of high-reynolds-number bubbles in inhomogeneous flows[END_REF] . Bubble shape is assumed to be spherical in this study. The assumption has been a posteriori validated by calculating from the simulations the bubble Weber number W e = ρ d | vu | 2 d/σ based on the slip velocity | vu | . Its maximum value has been found to be of the order of 10 -2 . Hence, the drag coefficient C D is calculated through the correlation of [START_REF] Mei | A note on the history force on a spherical bubble at finite reynolds number[END_REF]

C D = 16 Re b 1 + 8 Re b + 1 2 1 + 3 . 315 Re 1 / 2 b -1 (5)
where the bubble Reynolds number is defined as

Re b = | u -v | 2 R b ν , (6)
the added-mass coefficient is C M = 1 / 2 and the lift coefficient C L is given by the relation of [START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF] :

C L = 1 . 8802 ( 1 + 0 . 2 Re b /Sr ) 3 1 Re b Sr + 1 2 1 + 16 /Re b 1 + 29 /Re b 2 1 / 2 (7)
with Sr the local shear rate nondimensionalized by the relative velocity and the bubble diameter. Direct bubble/bubble interactions are neglected considering that the averaged void fraction in Section 3 and 4 is v = 0 . 02 and 0.001, respectively. However, the effect of concentration peaks on bubble/bubble interaction will be discussed in Section 4.1 . Bubble/wall overlap is prevented by assuming an elastic bouncing with the plane wall. The computation of these forces requires the interpolation of the fluid velocity and its time and space derivatives at each bubble location. We used second-order accuracy linear interpolation scheme. A third-order Runge-Kutta scheme is used for integrating the force balance in time with a typical time step equal to one fifth of the viscous relaxation time of the bubbles. Loop nesting is used when time steps of both phases are widely separated. When the Eulerian time step (based on numerical stability criteria) is smaller than the Lagrangian characteristic time then solutions of both sets of equations are synchronized. For large ratio of these time steps t Euler / t Lagrange , Lagrangian inner loops are integrated in a frozen flow field and the global computing time is dominated by the Lagrangian solver for O (10 5 -10 6 ) bubbles. Typically, we have t Euler / t Lagrange 20.

The Lagrangian solver JADIM has also been validated and used in the past for the simulation of bubbles induced convection [START_REF] Climent | Large-scale simulations of bubble-induced convection in a liquid layer[END_REF] , bubble migration in a pipe flow [START_REF] Legendre | Influence of gravity upon the bubble distribution in a turbulent pipe flow: comparison between numerical simulations and experimental data[END_REF] and bubble preferential accumulation in Couette-Taylor flow [START_REF] Chouippe | Numerical simulation of bubble dispersion in turbulent taylor-couette flow[END_REF][START_REF] Climent | Preferential accumulation of bubbles in couette-taylor flow patterns[END_REF] .

The two-way coupling model

To account for the feedback effect of the dispersed phase on the carrying fluid flow, we average the Navier-Stokes equations in a control volume of fluid populated by many bubbles. The basic principles of this operation can be found in several previous works such as [START_REF] Anderson | Fluid mechanical description of fluidized beds. equations of motion[END_REF][START_REF] Drew | Mathematical modeling of two-phase flow[END_REF][START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] . The averaging procedure is based on the function characterizing the presence of the continuous phase χ c locally which is unity if there is fluid at location x at time t and zero otherwise.

Then the local volume fraction of the continuous phase ε c ( x , t )

ε c (x , t ) = 1 V V χ c (ξ , t ) d 3 ξ = χ c V , (8)
where V is the volume of the cells. We note that the corresponding void fraction is computed as ε d (x , t ) = 1 -ε c (x , t ) . In our study, the void fraction associated with the control volume is estimated by adding the volume of each bubble contained in this control volume, which is estimated through the Lagrangian tracking of the bubbles. The volume averaging process leads to a new system of continuity and momentum equations:

∂ε c ∂t + U c • ∇ε c = -ε c ∇ • U c (9)
ε c ρ c ∂ U c ∂t + ε c ρ c U c • ∇ U c = -∇[ ε c P c ] + ∇[ µ c ε c ∇ U c + (∇ U c ) T ] + F d → c (10)
F d → c = - 1 V N b,l, (i, j,k ) b=1 ρ b V b ( d v b d t -g) . ( 11
)
The influence of the bubbles on the continuous phase is related to the local forcing induced by each bubble via the momentum source term F d → c and to the local evolution of the fluid volume fraction ε c . In those equations, we have neglected momentum transport by small-scale fluctuations (scales smaller than the control volume). The corresponding tensor u ′ i u ′ j c , similar to the Reynolds stress tensor, may have two main physical origins: the classical turbulent fluctuation of the continuous phase, and the bubble induced fluctuations due to finite size effects and wake interactions ( [START_REF] Riboux | Experimental characterization of the agitation generated by bubbles rising at high reynolds number[END_REF] ). The former contribution is not considered since we intend to perform direct numerical simulations of the continuous liquid flow while it is reasonable to neglect the latter when considering dilute bubbly flows. Note that, in some studies (see for example [START_REF] Climent | Large-scale simulations of bubble-induced convection in a liquid layer[END_REF][START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF][START_REF] Lakkaraju | Heat transport in bubbling turbulent convection[END_REF][START_REF] Mazzitelli | The effect of microbubbles on developed turbulence[END_REF] ), the fluid inertia is subtracted from [START_REF] Ferrante | Reynolds number effect on drag reduction in a microbubble-laden spatially developing turbulent boundary layer[END_REF] , because in these studies the volume occupied by the bubbles in the cell is not considered in the momentum balance equation ( 10) through ε c ( x , t ).

ε c (1) = ε c (2) = ε c (3) = ε c (4) .
The system of equation is discretized on a staggered nonuniform grid with a finite volume approach. The pressure and the liquid volume fraction are located at the same node and the velocities are face centered. In the cell ( i , j , k ) of volume V i, j,k , the value of

ε c , ( i , j , k ) is directly calculated from the number of bubbles N b , ( i , j , k ) present in the cell ( i , j , k ) by: ε c, (i, j,k ) = 1 -ε d, (i, j,k ) = 1 - N b, (i, j,k ) b=1 V b V i, j,k (12)
The point particle approach is well suited when the size of the bubbles is smaller than the smallest scale of variation of the fluid flow. Fig. 1 gives a schematic representation of a single bubble close to the wall. Grid size in streamwise and spanwise directions is much larger than the bubble diameter. Eq. 12 is valid for R b ≪ , but close to the wall y > R b . Although beyond the model's assumptions, it is unavoidable to simulate large Reynolds number flows with intermediate bubble size at a reasonable computational cost. The minimum possible distance of the bubble is y + = R + b from the wall. The value of ε c is based on the position of the center of the bubble, as a consequence, the cells close to the wall (where R b > y ) does not feel the presence of the bubbles ( ε c = 0 ). To remedy this problem and prevent the large wall-normal gradient of ε c in the vicinity of the wall we impose a constant value for ε c for the cells located between y = 0 and y = R b . We note that the approximation is rather crude but we tested that without any specific treatment for the value of ε c close to the wall and we obtained qualitatively similar results.

In Appendix B we provide various verification and validation of the numerical development proposed for this paper. Namely, B.1 presents validation of the single phase turbulent boundary layer simulations against the numerical simulations of [START_REF] Jiménez | Turbulent boundary layers and channels at moderate reynolds numbers[END_REF] . Then in B.2 and B.3 , we propose simple original test cases to check the two-way coupling method in order to verify each contribution of the terms induced by the presence of bubbles in the system of equations (9-10) .

Minimal-flow-unit configuration

As seen from the Eqs. 9 and 10 , the presence of bubbles affects the flow through different contributions. There is a void effect in both momentum and continuity equations as well as a direct contribution to the momentum coupling ( Eq. 11 ). The point of the section is to study the effect of each of these contributions (as well as their possible interplays) on the near wall flow structure. We consider in this section the minimal flow unit [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF] in which a dynamical process for the generation of the wall flow structures is extracted. More specifically, to study the interactions of the dispersed bubbly phase with the self-sustained wall structure mechanism, we study the following 5 cases : mfu0: One-way coupled simulation i.e. F d → c = 0 in eq 10 and ε c = 1 in eq 9 and eq 10 mfu1: Only momentum feedback from the dispersed to the fluid phase is considered i.e. F d → c = 0 in eq 10 and ε c = 1 in eq 9 and eq 10 mfu2: Only the void effect is considered in continuity equation ( 9) via ε c with F d → c = 0 in eq 10 and ε c = 1 in eq 10 mfu3: Only the void effect ( ε c = 1) is considered in momentum equation [START_REF] Ferrante | A robust method for generating inflow conditions for direct simulations of spatially-developing turbulent boundary layers[END_REF] with F d → c = 0 in eq 10 and ε c = 1 in eq 9 mfu4: The void effect ( ε c = 1) is considered both in continu- ity and momentum equations and the feedback ( F d → c ) term from the dispersed to continuous is also accounted for.

The description of the domain size and parameters used for the two-phase simulation is given in Table 1 and are similar to the configuration proposed in Jiménez and Moin [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF] . The Stokes number is defined as S t = τ b /τ * , where τ b is the bubble relaxation time defined as τ b = R 2 b / 6 ν and τ * = ν/u 2 τ is the time scale based on friction velocity u τ . We use periodic boundary conditions for the velocity in all three directions, and a mean pressure gradient is imposed in order to obtained the prescribed flow rate. Please note that the case of the mfu4 simulations have also been performed with a twice finer mesh in all directions, and that no sizable differences were observed. Therefore we conclude that the mesh used for the simulations described in the paper is fine enough to capture all the hydrodynamic fluctuations.

In case of one-way coupled minimal-flow-unit ( mfu0 ) the streamwise and the spanwise dimensions of the computational domain are L +

x ∼ 400 and L + z ∼ 180 respectively in wall units. The Reynolds number based on u τ corresponds to Re τ ∼ 180. Spacing between the near-wall structures is found to be around λ + z ∼ 80 -100 therefore mfu0 could have two pairs of high and low speed streaks. This size is chosen so that bubble dynamics could be studied during the interaction of these streaks. We set the bubble volume fraction at v = 0 . 02 along with a bubble radius of R + b ≈ 1 . 3 , corresponding to a Stokes number S t ≈ 0.28.

In Fig. 2 (a) we present the mean bubble concentration profile as a function of the wall normal coordinate. The bubble concentration is computed by averaging over time and directions parallel to the wall. We remark that the profile obtained for mfu0 and mfu1 are very close and decrease towards the wall whereas in the mfu 2, 3, 4 there is a sharp increase in the bubble concentration close to the wall. Since the different bubble distributions are only due to a modification of the carrier flow, we conclude that, at this volume fraction, void effects in both the continuity or momentum equations have much more pronounced effect on the flow than direct interphase momentum coupling. To analyze further the strong modification of the mean bubble concentration we give in figure 2 (b) the mean drag, lift and added mass forces acting on the bubble as a function of wall normal distance for the mfu0 and mfu4 . It appears that in mfu4 , lift and drag forces are modified when the void effect alters the carrier phase flow and there is a negative peak in the drag force close to the wall. As a consequence, there is a net force acting on the bubbles towards the wall which results in high bubble concentration close to the wall.

The second-order turbulent statistics are analyzed to further understand the effect of bubbles on the fluid phase. In Fig. 3 (a) root-mean-square velocity in the streamwise and wall normal direction is shown as a function of wall distance for all cases. In all cases except mfu1 turbulence intensity is higher than single phase case and the peak is slightly shifted towards the wall. Similarly, Reynolds stress profile is also modified due to the presence of the bubble, higher turbulent stress intensity is observed in mfu2,3,4 as compared to mfu0 and mfu1 ( Fig. 3 (b)). The peak of Reynolds stress is moved towards the wall which in turn increases the production of turbulent kinetic energy production near the wall.

Reynolds stress is determined by the contribution of the bursting and sweep events, therefore we investigate the effect of bubble injections on the bursting and sweep cycle. Contribution to the Reynolds stress from the four quadrants ( Fig. 4 ) is computed for mfu0 and mfu4 (see Section 4.4 for details). Negative contribution to the Reynolds stress (u v ) from the second and fourth quadrant accounts for the burst and sweep events. Close to the wall sweep event dominates and bursting is more prominent away from the wall. The crossover takes place at y + ≈ 15 for mfu0 . In mfu4 where drag increase was observed, the negative contribution from both sweep and burst events have increased and the point where they first cross each other is shifted towards the wall ( y + ≈ 14 ).

The positive contribution to -u v from the first and third quadrant is reduced in mfu4 as compared to mfu0 . The above test cases indicated that feedback force F d → c = 0 has minimal effect on the near-wall structure or drag modulation. The most important factor is the void effect in the continuity and momentum transfer equations, which plays a crucial role in the modification of the near-wall coherent structures and significantly increases the bubble concentration near the wall. Drag change is almost same in both mfu2 and mfu4 as reported in Table 2 . This similarity between mfu2 and mfu4 is also clear from rms and Reynolds stress (-u v ) profiles in Fig. 3 . From the above test cases it is clear that bubble concentration close to the wall has a dominant effect on the wall friction.

Modification of a turbulent boundary layer due to bubble injection

We now consider the case of the turbulent boundary layer. The flow is generated with the aid of the Synthetic Eddy Method introduced in Section 2 and detailed in Appendix A . Table 3 summarizes the specificities of the computational domain and discretization used in the current section. For more details on the validation of the simulation of the single-phase turbulent boundary layer the reader is referred to B.1 .

Influence of the bubble injection strategy

Contrary to the minimal flow unit considered in the previous section, in the real wall flows at large Reynolds number, there are multiple near-wall structures along with the large-scale structures. [START_REF] Mattson | Simulation of bubble migration in a turbulent boundary layer[END_REF] in their DNS of turbulent boundary layer with micro-bubbles showed that bubbles migrate away from the wall as they move along the streamwise direction and this decreases their ability to modulate the near-wall coherent structures. In Turbulent Boundary Layer (TBL) simulation of [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] , any time a bubble leaves the domain it is randomly reinjected inside the turbulent boundary layer. With this injection procedure, the bubble concentration profiles peak after y + = 100 and when the average volume fraction of the bubbles is low ( ε d = 0 . 001 ) there are very few bubbles near the wall (see Fig. 5 ). For this reason, the modification of the drag coefficient and the near-wall dynamics is minimal. Therefore, higher bulk bubble volume fraction is required to generate some significant changes in drag. The alternative to enhance the bubble interaction with nearwall coherent structures for a fixed global volume fraction is to inject bubbles at different heights close to the wall.

We consider five cases of bubble injections as summarized in Table 4 . In the first four cases (denoted: 2R b , 4R b , 6R b and 8R b ) bubbles are injected close to the wall at different y + . In the fifth case (named: Random ), the bubbles are spontaneously nucleated at random positions in the whole domain with a velocity equal to the local fluid velocity. In addition those simulations are compared with a single phase flow whose parameters are given in Table 3 . In all cases the overall bubble volume fraction within the boundary layer is kept fixed at v = 0 . 001 . Note that in the case denoted 2R b the bubble radius is one fourth of the radius of the other three cases, this was done to investigate the effect of injecting the bubbles very close to the wall (the injection position has to be larger than the bubble radius). However, note that since the bubble diameter is already smaller than the turbulent wall structures we do not expect a very important effect of the bubble diameter as long as the gas volume fraction is constant. The case Random is similar to the case-B of [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] in which whenever a bubble leaves the computational domain it is reinjected at a random location inside the boundary layer. In other cases when the bubbles leave the domain they are reinjected at random positions on the xz plane located at the height of 2 R b , 4 R b , 6 R b and 8 R b from the wall respectively. Note that for the bubble diameter considered here, it can be checked that for air bubbles in water and for a realistic boundary layer thickness the Weber number is actually vanishingly small, consistently with the assumption of spherical bubbles considered in Section 2.2 .

Figure 5 (a) presents the profile of the average bubble volume fraction of the various cases. The profiles are computed from spatial average in spanwise direction and time averaging and are given for a streamwise location x = 10 δ 0 , for which Re θ ≈ 1100. Bubble concentration profile peaks at the location of the bubble injection and then gradually decreases as we move away from the wall. This explains that the concentration profiles differ from the trend of the literature ( [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF][START_REF] Mattson | Simulation of bubble migration in a turbulent boundary layer[END_REF] ) with much more intense peaks that can be also located in the buffer layer. The peak of bubble concentration in the near wall region suggests the existence of bubble/bubble interactions, which have been neglected in the paper. To discuss the impact on the bubble dynamics of the induced bubble/bubble interactions in the near wall region, we consider the bubble-bubble distance compared to the characteristic length of the turbulent structures. The minimum average inter-bubble distance, in wall unit, is

expressed as δ + b = R + b 4 π 3 ε d 1 / 3
. Taking the value of the local concentration peak, ε d ∼ 0.1 (see fig. 10 ) gives: δ + b ≈ 5 , which remains smaller than the wall structures in the very near-wall region. This indicates that the bubbles interact much more with the turbulent structures than with other bubbles, and that bubble-bubble interactions are not expected to be dominant in this region.

The evolution of the temporal and spanwise average skin friction coefficient C f (x ) = 2 τ w (x ) /ρU 2 in the streamwise direction is shown in figure 5 (b), here τ w = µd U/d y . In the case of 2R b and 4R b the drag is smaller than in the one-way coupled case, for case

Random the drag seems to be unchanged, whereas the drag increases in cases 6R b and 8R b . This increase can be interpreted as the consequence of the wall-normal location of bubble injection.

Note that the small variations between the cases 2R b and 4R b indicate that the main parameter for the alteration of the wall structures is the local volume fraction of the near wall region. The bubble diameter and the bubble number density separately are not dominant. This is consistent with the previous section where it is observed that the bubbles mainly influence the wall through void fraction effect and not direct momentum coupling. This conclusion only holds for small enough bubbles, for large bubbles the momentum forcing caused by each bubble should prevail. The induced effect of the injection position is further detailed in the following considering the Reynolds stress and its analysis using quadrant decomposition.

Turbulent fluctuations and energy-Spectrum

The turbulent statistics ( u rms = u ′ 2 ) computed for the cases discussed above are shown in Figs. 6 (a) and(b). The velocity components are normalized by the free stream velocity U 0 . In cases 2R b , 4R b streamwise and wall-normal turbulent energy fluctuations are slightly lower than the Single-phase case, whereas they are higher in 6R b and 8R b . These findings are consistent with observed modification of the skin-friction coefficient. Similarly, a reduction of the turbulent fluctuation was reported by [START_REF] Ferrante | Reynolds number effect on drag reduction in a microbubble-laden spatially developing turbulent boundary layer[END_REF] and [START_REF] Pang | Numerical study on modulation of microbubbles on turbulence frictional drag in a horizontal channel[END_REF] in their DNS of turbulent boundary layer and channel flow simulations laden with micro-bubbles. However, unlike them we do not observe any noticeable shift in the rms peak along the wallnormal direction. It is possible that because the drag reduction observed here is less than 4 percent, the shift in wall-normal direction is negligible.

To analyze the cause of drag modification due to this injection strategy pre-multiplied spanwise energy spectrum is analyzed for cases 2R b , 4R b , 6R b , 8R b and Single-phase case. Figures 7 (a) and 7 (b) show pre-multiplied spanwise energy spectra for streamwise ( K z E uu ) and spanwise velocity components ( K z E ww ) at three different heights from the wall ( y + = 4 . 0 , 10 , 15 ). In all the cases the spectral energy peaks at the wavelength of λ + z ≈ 80 -100 and this peak shifts gradually towards higher wavelengths as we move away from the wall. This wavelength corresponds to the characteristic spacing between the streaks close to the wall. However, 2R b and 4R b have lower energy content than the Singlephase simulation, whereas 6R b and 8R b has higher energy content. As explained in Section 1.1 near-wall cycle consists of streaks which after gaining sufficient energy become unstable and break down to generate vortices. In cases 6R b and 8R b higher energy streaks are more prone to instability which results in their eventual breakdown into vortices. Increase in instability of the coherent structures leads to more violent burst and sweep events, which results in the increase in turbulent shear stress production (see Fig. 9 ). Opposite could be said about the case 2R b and 4R b where energy in the coher- ent structures is reduced which makes them less unstable and decreases turbulent shear stress production (see figure 9 ).

Shear stress and turbulent energy production

We examine the Reynolds stress ( -u v ) and turbulent energy production term ( P = -u v d U/d y ) in Figs. 8 (

Quadrant analysis: Bursting and sweeping events

Intensity and frequency of bursting and sweep events are the main contributor to the Reynolds stress. In order to understand the modification of Reynolds stress, the influence of micro-bubbles on the bursting cycle is investigated using quadrant analysis. Quadrant analysis computes the fractional contribution of the turbulent fluctuations to the Reynolds stress at each point in the domain. The uv plane is divided into four quadrants; Q1 represents u > 0 , v > 0 , where high-speed fluid rushes away from the wall, Q2 contains u < 0 , v > 0 , where low-speed fluid moves away from the wall, this is generally called bursting, Q3 represents u < 0 , v < 0 events where low speed fluid moves towards the wall and Q4 contains u > 0 , v < 0 events where high-speed fluid rushes towards the wall, usually referred as sweep events. Therefore, Q2 and Q4 contribute to the negative part of Reynolds stress and Q1 and Q3 bring a positive contribution. Contribution from Quadrant 2 and 4 is generally much higher than Quadrant 1 and 3.

In Fig. 9 we show the fractional contribution of these four quadrants as a function of wall-normal coordinate at the streamwise location of x = 6 δ 0 for the first four cases in Table 4 . Very close to the wall, sweeping events (Q4) is a leading contributor to the turbulent Reynolds stress. For the single phase case, at around y + ≈ 14 the contribution of Q2 becomes equal to Q4 and further away bursting events dominate the Reynolds stress. For cases 6R b and 8R b this intersection point for Q4 and Q2 mildly shifted towards the wall to y + ≈ 13 . 6 & 13.4 respectively. On the other hand, for cases 2R b , 4R b this intersection point is moved away from the wall to y + ≈ 14 . 6 & 15.11 respectively. In cases 6R b and 8R b the intensity of sweep and ejection is enhanced as compared to all the other cases. Whereas in cases 2R b , 4R b contribution to the Reynolds stress due to sweeping and bursting events is reduced. The positive contribution to the Reynolds stress from quadrant Q1 and Q3 is enhanced in cases 2R b and 4R b and reduced for cases 6R b and 8R b .

Thus, for cases 6R b and 8R b the bubble injection increases the intensity of bursting and sweeping events thereby increasing the Reynolds shear stress. Whereas, in cases 2R b , 4R b , there is a decrease of the negative contribution from Q2 and Q4 and an increase of the positive from Q1 and Q3, which results in the decrease in total Reynolds stress. As burst and sweeps are located respectively above the low and high speed near-wall streaks i.e Fig. 9. Single-phase , 2R b , 4R b , 6R b and 8R b cases are compared. Panels (I), (II), (III) and (IV) are arranged according to the four quadrants of the u-v plane. Color coding is similar to Fig. 7 . adjacent to the quasi-streamise vortices (see [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] ) change in their intensity must be correlated to the interaction of bubbles with the near-wall coherent structures.

Bubble dispersion: λ 2 criterion

In this section, preferential accumulation of bubbles induced by the coherent structures is examined. Bubbles and vortical structures mutually interact with each other in a turbulent flow. Lift force, local pressure gradient and added mass force encourage the accumulation of bubbles in low pressure regions of the flow [START_REF] Felton | Spherical bubble motion in a turbulent boundary layer[END_REF][START_REF] Mazzitelli | The effect of microbubbles on developed turbulence[END_REF] . Such low pressure regions, which often correspond to the vortex cores, are classically identified by λ 2 criteria ( [START_REF] Jeong | On the identification of a vortex[END_REF] ). It corresponds to negative values of the second-largest eigenvalue of S 2 + 2 where S and are the symmetric and antisymmetric parts of the velocity gradient tensor ∇u . Figure 10 shows the contour of λ 2 < 0 at 10% of its absolute value and the bubble positions at a particular instant. One can observe that bubbles are preferentially accumulated in λ 2 < 0 regions. We further compare the probabil- ity ( P b

Table 4

Bubble radius and distance to the wall of the injection plane for the different run. For all cases the overall bubble volume fraction within the boundary layer is shifts away from the wall, while in case 6R b and 8R b it slightly shifts towards the wall, and for the Random case no such shift was observed. This indicates that the vortical structures are displaced to higher wall normal locations when the drag is reduced and they move towards the wall in case of drag increase. Overall, the probability

v = 0 . 001 . Case 2R b 4R b 6R b 8R b Random R + b 0.3 1.3 1.3 1.3 1.3 y +
P b λ 2 > P f λ 2
demonstrates that bubbles in turbulent flow get preferentially accumulated in the vortex cores, which is consistent with the observations of [START_REF] Climent | Preferential accumulation of bubbles in couette-taylor flow patterns[END_REF] and [START_REF] Chouippe | Numerical simulation of bubble dispersion in turbulent taylor-couette flow[END_REF] .

At the distance to the wall corresponding to the bubble injection there is a sudden variation in P b λ 2 but further away, say y + > 10 , bubbles accumulate in λ 2 < 0 regions. Moreover, the probability distribution P b λ 2 becomes similar away from the wall for all the cases corresponding to bubble injection at the wall ( 2R b , 4R b , 6R b and 8R b ). This is in contrast with the Random case which deviates from the other cases. This is due to spontaneous nucleation of bubble randomly in the whole domain.

Conclusion

We develop further the numerical tools necessary to simulate a spatially developing turbulent boundary layer laden with microbubbles. We have implemented a numerical approach which enables to compute bubble dispersion with two-way coupling based on Euler-Lagrange approach. The formulation of the feedback of the dispersed phase on the continuous one is based on local averaging that modifies the Navier-Stokes equations. To study the physical mechanisms at play in the alteration of the wall flows by the presence of micro-bubbles, we considered two distinct situations. First we considered a minimal flow unit laden with micro-bubbles. This model flow enables to isolate the self-sustain process of generation of the wall structures. Our results demonstrate that the bubbles tend to migrate close to the wall mainly due to the occurrence of sweep events. This migration can be attributed to the relative low Reynolds numbers explored for those configurations and the simplified dynamics of the minimal flow that might emphasize the influence of those sweep events. Moreover, we notice that the presence of the bubbles tends to reinforce the occurrence of such violent events. Also it appears that the bubbles influence the flow mainly through effects caused by the void fraction fluctuations (both in momentum and continuity equations) and not through direct momentum coupling. Since most of the effect on the wall shear stress appears to be due to the bubbles in the close vicinity of the wall, we discussed the influence of bubble injection methods on the modulation of the wall shear stress of a turbulent boundary layer. We observe a change in the flow dynamics as well as a modification of the skin friction, even for a relatively small bubble volume fraction injected in the near wall region. The key parameters seem to be the void fraction and the injection position whereas the bubble size is less important, at least as long as bubbles are smaller than or of the order of the friction length. It seems that the drag is decreased when bubbles are injected in the viscous sub-layer while it is increased when they are injected the buffer layer. We show that this drag behavior is related to a change of the near wall flow structures, and thus of their regeneration cycle, depending on the bubbles injection location. The resulting change in the intensity of the so-called "sweep" and "burst" causes a modification of the mean streamwise momentum flux in the wall normal direction, as confirmed by the Reynolds shear stress profiles resulting in the change of the skin friction. However, it remains an open problem in turbulence how to quantitatively relate the alteration of the regeneration cycle (or at least the evolution of the different quadrants contribution) to the overall wall shear stress.

˜

u i (x ) = 1 N E N E k =1 ǫ k f σ (x ) (x -x k ) (A.1)
where ǫ k accounts for the positive or negative contribution to the synthetic field and has equal probability to take value ± 1. Here the shape function for the eddies is defined as

f σ (x ) (x -x k ) = V B 1 σ x σ y σ z f x -x k σ x f y -y k σ y f z -z k σ z (A.2)
f has a compact support on [ -1 , 1] and satisfies the normalization condition In the present study f is chosen to be a tent function,

f (x ) = 3 2 (1 -| x | ) , x ≤ 1 0 , otherwise (A.4)
Here σ x , σ y , σ z control the size of the eddy generated. It is ad- vised to choose the length scales close to the one found in real turbulent flow as it strongly affects the turbulent fluctuations generated downstream of the inlet (see [START_REF] Pamies | Generation of synthetic turbulent inflow data for large eddy simulation of spatially evolving wall-bounded flows[END_REF] ). However, different choices for length scales have been investigated in [START_REF] Jarrin | Synthetic inflow boundary conditions for wall bounded flows[END_REF] and [START_REF] Jarrin | Reconstruction of turbulent fluctuations for hybrid RANS/LES simulations using a synthetic-eddy method[END_REF] and sensitive dependence of the results on them is discussed in detail. In the present study we have chosen length scales to be isotropic ( σ x = σ y = σ z = σ (y ) ) and to vary with the wall normal direction ( y ) based on Prandtl's mixing-length hypothesis i.e. σ (y ) = k y . However, due to the grid stretching in the wall-normal ( y ), close to the wall the size of the eddy becomes too small to be discretized in spanwise and streamwise direction, therefore the following relation is used to estimate the size of the eddies

x k (t + dt ) = x k (t) + U c dt (A.6)
where dt is the time step of the simulation. When an eddy leaves the box B , a new eddy is generated on the inlet of the box with random spanwise and wall-normal locations.

Once ˜ u is generated on the inlet plane with the appropriate temporal and spatial correlations, one obtains a velocity field u with the prescribed mean and covariance profiles, by linear transformation of ˜ u , (providing that ˜ u i ˜ u j = δ i j ):

u i (y, z, t ) = U i (y ) + A i j (y ) ˜ u j (y, z, t ) (A.7)
where U i ( y ) is the mean velocity profile and A ij ( y ) is the Cholesky decomposition of the Reynolds stress tensor R ij ( y ) profile:

A i j =   R 11 0 0 R 21 /A 11 R 22 -A 2 21 0 R 31 /A 11 (R 22 -A 21 A 31 ) /A 22 R 22 -A 2 21 -A 32   (A.8)
Because with SEM technic only the largest turbulent structures are generated, the flow presents a transient region downstream the inlet plane before reaching the fully realistic turbulent correlations and structures. Eddies of smaller size were found to undergo a long transient and therefore, require a longer computational domain to reach higher Reynolds number. Moreover, away from the wall the smaller structures dissipate rather than to evolve into large-scale structures resulting in unphysical velocity field. However, when Prandtl's mixing-length hypothesis (see Eq. A.5 ) is used to provide inhomogeneity in the wall normal direction, the largescale structures away from the wall break down and energy cascades to smaller scales, thereby proving more physical eddy distributions.

To prevent the interaction of bubbles with unphysical structures due to the transient, we add bubbles after the transient and disregard the transient length in further analysis. To save computational time the simulation is partitioned into two separate simulations ( [START_REF] Ferrante | A robust method for generating inflow conditions for direct simulations of spatially-developing turbulent boundary layers[END_REF] ). They are referred as Domain-A and Domain-B simulations. The advantage of this setup is that we have to perform Domain-A simulation only once to make a database of more realistic inlet boundary conditions. Domain-B , when provided the inlet condition issued from the Domain-A , simulation does not need any transition length to become turbulent. The strategy implemented could be summarized in the following steps:

1. Domain-A utilizes SEM to generate boundary condition at the inlet plane 1 .

2. At every time steps of Domain-A run, the velocity field in the yz plane 2 located at streamwise location 12 δ is saved to the disc. 3. Domain-B reads this velocity field every time step as an inlet boundary condition.

Appendix B. Verification and validation

We present the different test cases used to verify and validate the new implementations in our JADIM code. We first report the validation of the turbulent boundary layer simulation in B.1 . Then in B.2 the momentum forcing term F d → c has been tested without any other void fraction effect by setting ε c = 1 (case 1), and eventually in B.3 the different contributions of the variation ε c in both the continuity and the momentum equations are verified without considering the effect of F d → c (cases 2 to 6).

B1. Validation of the turbulent boundary layer simulation

We report here the validation of the turbulent boundary layer simulation. As explained, in Section Appendix A , the synthetic eddies generated need some streamwise distance to evolve and provide completely physical turbulent structures. Therefore, to avoid

Table B1

Relative error for the pressure gradient as a function of the mean void fraction. the interaction of bubbles with these un-physical structures and save computational time, we have broken the simulation into two parts Domain-A -Domain-B ( Fig. A2 ). In Domain-A we performed a single phase simulation in order to generate an inlet boundary condition database that will be subsequently used for the two-phase flow simulations carried out in the Domain-B .

Table 3 presents the dimension and the mesh resolution of the Domain-B . Note that Domain-A is three times longer in streamwise direction than Domain-B and of the same size in spanwise and wall-normal directions.

Figure B1 shows the statistics of turbulence for the single phase flow compared against the DNS results of [START_REF] Jiménez | Turbulent boundary layers and channels at moderate reynolds numbers[END_REF] at Re θ ≈ 1100 where Re θ = U 0 δ θ /ν. These statistics are computed by temporal and spanwise average at a given streamwise location. The streamwise mean velocity profile is in good agreement with [START_REF] Jiménez | Turbulent boundary layers and channels at moderate reynolds numbers[END_REF] . Similarly, the profiles of the variance of the streamwise and wall-normal velocity components are in good agreement, but, for the spanwise velocity component we observe a slight discrepancy which is likely due to the smaller spanwise width considered in our simulations.

B2. Verification of the implementation of the momentum forcing term F d → c

We first verify the momentum forcing term F d → c by considering the pressure drop generated by a random distribution of fixed identical bubbles. We consider a 3D channel as shown in the figure B2 . The dimensions of the channel are L x , L y and L z in the x , y and z directions, respectively. They are chosen such that L x = 5 L y and L z = 1 . 175 L y . N b bubbles of volume V b are randomly distributed and maintained fixed in the central part of the domain In this case, we do not consider the effect of the void fraction in the system of equations [START_REF] Ferrante | On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles[END_REF][START_REF] Ferrante | A robust method for generating inflow conditions for direct simulations of spatially-developing turbulent boundary layers[END_REF] . By imposing the momentum source term -F b e x for each bubble, the pressure jump between the inlet pressure P in and the outlet pressure P out is given by: where S = L y L z is the channel section. In fact the pressure gradient is constant through the bubble cloud and is proportional to the void fraction α. This is confirmed by Fig. B3 that reports the normalized pressure gradient P/ 0 . 5 ρ f U 2 in as a function of the normalized position x / L x in the channel. The relative error of the pressure jump is reported in Table B1 for the three void fractions considered. As shown agreement is very good.

( P out -P in ) S = -N b F b (B.1) P (x ) -P in 0 . 5 ρ f U 2 in = - L x F b α V b ρ f U 2 in (B.2)
∂ ε c U k ∂x k = 0 2 ρ f ε c ∂ U i ∂t + U k ∂ U i ∂xk = ε c (x ) = ε in / (1 + αx ) U (x ) = U in (1 + αx ) P(x ) = -1 2 ρ f (U 2 (x ) -U 2 (L x )) -∂P ∂xi + µ ∂ ∂x k ∂U i ∂x k + ∂U k ∂xi ∂ ε c ∂t + ∂ ε c U k ∂x k = 0 3 ρ f ε c ∂ U i ∂t + U k ∂ U i ∂x k = ε c (t) = ε in (1 -t/T o ) U (x ) = U in + x t-T o P(x ) = -ρ f 1 2 U 2 -U 2 (x = L x ) -∂P ∂xi + µ ∂ ∂x k ∂U i ∂x k + ∂U k ∂xi -ε in 2 T o (T o -t ) x 2 -L 2 x ∂εc U k ∂x k = 0 4 ρ f ε c ∂ U i ∂t + U k ∂ U i ∂xk = ε c (x ) = ε in / (1 + αx ) U (x ) = U in (1 + αx ) P(x ) = -1 2 ρ f U in ε in [(1 + αx ) 2 -∂ ε c P ∂xi + µ ∂ ∂x k ∂U i ∂x k + ∂U k ∂xi -(1 + αL x ) 2 ](1 + αx ) ∂ εc U k ∂x k = 0 5 ρ f ε c ∂ U i ∂t + ε c U k ∂ U i ∂xk = ε c (x ) = ε in / (1 + αx ) U (x ) = U in (1 + αx ) P(x ) = -ρ f αU 2 in (x -L x )(1 + αx ) -∂ ε c P ∂xi + µ ∂ ∂x k ∂U i ∂x k + ∂U k ∂xi ∂ ε c ∂t + ∂ ε c U k ∂x k = 0 6 ε c U k ∂ U i ∂xk = ε c (x ) = ε in / (1 + αx ) U (x ) = U in (1 + αx ) P(x ) = -ρ f ε in U 2 in α(x -L x ) /ε(x ) -∂ ε c P ∂xi + µ ∂ ∂x k ε c ∂Ui ∂x k + ∂U k ∂xi +2 µαU in ε(x ) -ε L x ε(x )

B3. Verification of the implementation of the void fraction contribution

In order to check the effect of the presence of the dispersed bubbly phase via ε c in the system of equations (9-10) , we have defined several test cases reported in Table B2 . This table gives the system of equations considered, the imposed expressions of ε c and the corresponding exact solutions for both the velocity and the pressure. As shown, the complexity of the tests is progressive. For all the cases, we proceed as following: we consider the configuration of case 1 (same mesh, as well), the inlet velocity is U = U in e x at x = 0 and we imposed simple variations of ε c , in time or in space, so that the exact solution of the velocity and pressure evolutions can be derived and compared to the simulations. With cases 2 and 3, we tested the modification of the continuity equation. The aim of case 4 was to verify the modification of the pressure gradient while the modification of the advective term is tested with case 5. Finally, the modification of the viscous contribution into the momentum equation is verified using case 6. For cases 2 to 5, we neglected the viscous contribution while for case 6 it has been accounted for with Re = 0 . 12 . Fig. B4 gives the evolution of the volume fraction of the fluid, the velocity and the pressure. Table B3 summarizes the errors of the computed fields compared to analytic solutions. The error is always below O (10 -10 ) for the velocity and O (10 -4 ) for the pressure which verifies the numerical implementation.

Fig. 1 .

 1 Fig. 1. Schematic representation of bubble close to the wall. Cell 4 contains the center of the bubble, the contribution due to this bubble to εc is accounted at cell 4 only. Cell 1, 2 and 3 also share some portion of this bubble therefore we set

Fig. 2 .

 2 Fig. 2. (a) mean bubble concentration profiles ( ε d ( y )) as a function of wall normal coordinates for mfu0, 1, 2, 3, 4 . (b) mean lift, drag and added mass forces acting on the bubbles for mfu0 (black) and mfu4 (red) normalized by the friction velocity of the mfu0 case. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .

 3 Fig. 3. (a) Comparison of u rms (y + ) and v rms (y + ) for the mfu0, 1, 2, 3, 4 cases. Each curve is normalized by the friction velocity of the mfu0 case. (b) Reynolds stress profiles -u v normalized by the friction velocity of the mfu0 case.

Fig. 4 .

 4 Fig. 4. Contribution to the Reynolds stress from the four quadrants. black lines: mfu0 and red lines: mfu4 . Each curve is normalized by the friction velocity of the mfu0 case. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .

 5 Fig. 5. Panel (a) on the left compares the mean bubble concentration profiles, panel (b) shows the streamwise evolution of the friction coefficient C f ( x ), the solid black line corresponds to the single phase flow.

Fig. 6 .

 6 Fig. 6. Comparison of the velocity RMS for Single-phase , 2R b , 4R b , 6R b and 8R b cases. (a) variance of the streamwise u rms and (b) wall-normal v rms velocity components normalized by the free stream velocity U 0 as a function of wall-normal coordinate ( y ). Color coding is similar to Fig. 7 .

Fig. 7 .Fig. 8 .

 78 Fig. 7. Pre-multiplied spanwise energy spectrum comparison for Single-phase , 2R b , 4R b , 6R b and 8R b case at y + = 4 . 0 , 10 , 15 . Panel (a) on the left shows pre-multiplied streamwise energy spectra ( K z E uu ) for streamwise velocity component and Panel (b) on the right shows pre-multiplied spanwise energy spectra ( K z E ww ) for streamwise velocity component. To better discern the curves at y + = 10 and 15 for K z E uu we have plotted 5 • 10 3 + K z E uu at y + = 15 . Similarly, 1 • 10 3 + K z E ww is plotted at y + = 15 for K z E ww .

  a), (b) respectively. It is apparent that when compared to Singlephase case the Reynolds stress is reduced for 2R b , 4R b cases and remains almost unchanged for Random case, whereas for 6R b and 8R b it increases. Reynolds stress profile is mildly shifted away from the wall in cases 2R b and 4R b compared to Singlephase case. However, for cases 6R b and 8R b the Reynolds stress profile shifts towards the wall. Reynolds stress profile peaks at y + ≈ 69 in Singlephase case, it shifts to y + ≈ 73 for 2R b and 4R b and to y + ≈ 64 for 6R b and 8R b . The shift in Reynolds stress away from the wall could be attributed to the shift in vortex cores to higher y location. The opposite trend occurs for 6R b and 8R b where the Reynolds stress peak is moved to lower y location inferring the movement of vortical structures towards the wall. Turbulent energy production P term appears to be decreased for cases 2R b and 4R b and increased for 6R b and 8R b when compared to Singlephase case.

Fig. 10 .

 10 Fig. 10. Contour of λ 2 < 0 is shown at 10% of its absolute value in light green color up to streamwise distance of x = 8 δ. Bubbles positions are shown with black spheres. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

λ 2 ) 2 ) 2 and P b λ 2 as

 2222 of bubbles to sample the region of the fluid with λ 2 < 0 against the probability of the fluid ( P f λ having λ 2 < 0. Fig.11 (a)shows the probabilities P f λ a function of wall normal distance for all the cases listed in

Fig. 11 .

 11 Fig. 11. Probability of occurrence of λ 2 < 0 for cases listed in Table 4 and Singlephase case. P f λ2 (y ) and P b λ2 (y ) are shown with lines and points respectively.

20 . 2 peaks at y + ≈ 19 . 5 ,

 202195 More specifically, for the Single-phase case P f λ whereas in cases 2R b , 4R b , 6R b and 8R b the peak respectively shifts to y + ≈ 21 . 2 , 21 . 5 , 16 . 8 and 16.7. For the cases 2R b and 4R b where drag reduction was observed, the peak of the probability distribution P f λ 2

1 - 1 f 2 1 (A. 3 )

 11213 Fig. A1. Pictorial representation of synthetic eddy method. Purple square represent the inlet plane where the initial turbulent condition is supposed to be generated. Black domain is the boundary of the virtual box whose size is dependent on the length scales of the eddies generated.

Fig. A2 .

 A2 Fig. A2. Schematic representation of Domain-A -Domain-B configuration.

  σ (y ) = max { k y, } (A.5) where = max ( x, y, z) and k = 0 . 41 is the Von Kármán constant. The eddies generated inside the box B are convected with constant velocity U c characteristic of the flow using Taylor's frozen turbulence hypothesis. At each time step the new positions of the eddies are given by

Fig. B1 .

 B1 Fig. B1. (a) Mean flow profiles U + (y + ) ; (b) profile of the variance of the three velocity components u + rms , v + rms and w + rms ; (c) Reynolds stress -u ′ v ′ + profile. Comparison with

α 5 . 1 ×

 51 10 -4 5 . 1 × 10 -3 5 . 1 × 10 -2 Relative error for P 5 . 4 × 10 -6 7 . 3 × 10 -6 6 . 3 × 10 -6

Fig. B2 .

 B2 Fig. B2. Sketch of the geometrical configuration used for the verification in the case 1 with N b = 40 , 0 0 0 uniformly distributed bubbles of radius R b = 2 × 10 -4 L y .

Fig. B3 .

 B3 Fig. B3. Evolution in the x direction of (a) the pressure averaged in the channel section, and (b) the same pressure scaled by the global void fraction.

Fig. B4 .

 B4 Fig. B4. Summary of the mean fluid concentration (left), pressure (middle) and velocity (right) for the different cases in the x direction (for case 3 it corresponds to time t = 0 . 05 L y /U in ).

Table 1

 1 Computational domain and discretization parameters used in the minimal-flowunit box at Re τ = 180 .

	Name mfu0	L x 2.2h	L y 2h	L z 1h	N x 32	N y 82	N z 24	S t 0.28	v 0.02	R + b 1.3	Re τ 180

Table 2

 2 Drag modulation in different cases.

	Name	mfu0	mfu1	mfu2	mfu3	mfu4
	C f -C f, 0 C f, 0	0%	0%	8.5%	10.3%	8.5%

Table 3

 3 Computational domains and discretization parameters used in turbulent boundary layer simulation.

	L x	L y	L z	N x	N y	N z	+ x	y + min	+ z
	10 δ 0	2 δ 0	3.56 δ 0	360	66	136	18	0.14	9.0

Table 4 .

 4 Basically, for all the case P

f λ 2 presents a similar profile, with a minimum around

Table B2

 B2 Verification tests and corresponding analytic solutions for non-uniform void fraction distributions.

	Case	Governing Equations	Volume fraction	Velocity	Pressure
			(imposed)	(Exact solution)	(Exact solution)

Table B3

 B3 Errors obtained for the different tests reported in the Table B.6 .

	Case	Max relative error for U	Max relative error for P	Max error for V	Max error for W
	2	4 . 0 × 10 -4	4 . 80 × 10 -4	10 -11	10 -11
	3	3 . 80 × 10 -4	3 . 30 × 10 -3	10 -14	10 -14
	4	4 . 5 × 10 -4	5 . 7 × 10 -4	10 -11	10 -12
	5	1 . 0 × 10 -4	4 . 0 × 10 -4	10 -11	10 -11
	6	4 . 0 × 10 -4	3 . 1 × 10 -4	10 -13	10 -13

Acknowledgments

We acknowledge the support of the french ANR (Project number ANR-12-ASTR-0017). This work was performed using HPC resources from GENCI-CINES/IDRIS and CALMIP center of the University of Toulouse. We also thank Annaïg Pedrono for her technical support.

Appendix A. The Synthetic Eddy Method ( SEM )

The Synthetic Eddy Method SEM technique is based on the hypothesis that turbulence can be modeled as superpositions of coherent structures. [START_REF] Perry | A wall-wake model for the turbulence structure of boundary layers. part 1. extension of the attached eddy hypothesis[END_REF] suggested that superposition of analytical eddies with shapes inspired by real coherent structures provides a decent approximation of lower-order statistics of wall-bounded flows. Therefore, SEM deals with generating synthetic eddies at the inlet plane using specific shape function which contains their spatial and temporal characteristics. The inlet plane is located at x = 0 and has physical dimensions [0, L y ] × [0, L z ] in wall-normal and spanwise directions, while the full domain has the length L x in the streamwise direction. The grid is uniform in the streamwise ( x ) and spanwise ( z ) directions with resolution x and z respectively and stretched in wall-normal direction ( y ) using a tangential stretching function, thus y increases away from the wall.

As a first step a three-dimensional virtual box B containing the eddies is created around the inlet. The dimensions of the box B are [ -l x , l x ] × [ -l y , L y + l y ] × [ -l z , L z + l z ] and eddies are generated at some randomly chosen locations inside this domain. The volume of the box is denoted by V B and the number of eddies produced is N E .

The stochastic velocity signal ˜ u at a point x is assumed to be the superimposition of the contribution from all the eddies N E in the box B ,