Energy balance in lubricated drag-reduced turbulent channel flow

Abstract We use direct numerical simulation (DNS) to study drag reduction in a lubricated channel, a flow instance in which a thin layer of lubricating fluid is injected in the near-wall region so as to favour the transportation of a primary fluid. In the present configuration, the two fluids have equal density but different viscosity, so that a viscosity ratio $\lambda =\eta _1/\eta _2$ can be defined. To cover a meaningful range of possible situations, we consider five different $\lambda$ in the range $0.25 \le \lambda \le 4$. All DNS are run using the constant power input (CPI) approach, which prescribes that the flow rate is adjusted according to the actual pressure gradient so as to keep constant the power injected into the flow. The CPI approach has been purposely extended here for the first time to the case of multiphase flows. A phase-field method is used to describe the dynamics of the liquid–liquid interface. We unambiguously show that a significant drag reduction (DR) can be achieved for $\lambda \le 2.00$. Reportedly, the observed DR is a non-monotonic function of $\lambda$ and, in the present case, is maximum for $\lambda = 1.00$ (${\simeq }13\,\%$ flow-rate increase). Upon a detailed analysis of the energy budgets, we are able to show the existence of two different DR mechanisms. For $\lambda =1.00$ and $\lambda =2.00$, DR is purely due to the effect of the surface tension – a localized elasticity element that separates the two fluids – which, decoupling the wall-normal momentum transfer mechanisms between the primary and the lubricating layer, suppresses turbulence in the lubricating layer (laminarization) and reduces the overall drag. For $\lambda < 1.00$, turbulence can be sustained in the lubricating layer, because of the increased local Reynolds number. In this case, DR is simply due to the smaller viscosity of the lubricating layer that acts to decrease directly the corresponding wall friction. Finally, we show evidence that an upper bound for $\lambda$ exists, for which DR cannot be observed: for $\lambda =4.00$, we report a slight drag enhancement, thereby indicating that the turbulence suppression observed in the lubricating layer cannot completely balance the increased friction due to the larger viscosity.


Introduction
There is firm evidence (Joseph, Renardy & Renardy 1984;Hu, Lundgren & Joseph 1990;Joseph et al. 1997;Kim & Choi 2018;Roccon, Zonta & Soldati 2019) that the injection of a small amount of a low viscosity fluid into a pipeline used to transport a high viscosity fluid produces drag reduction (DR). The key physical principle at the heart of the observed DR mechanism is the natural tendency for the low viscosity fluid to migrate to the pipe wall and to lubricate the flow (Joseph et al. 1997). Known since the beginning of the last century (Isaac & Speed 1904;Looman 1916;Clark & Shapiro 1949), this DR mechanism has received a lot of attention, precisely because of its potential applicability to the strategically and industrially relevant case of water-lubricated oil pipelines (Russel & Charles 1959;Charles, Govier & Hodgson 1961;Hasson, Mann & Nir 1970).
Literature in the field is vast (see Joseph et al. (1997), Ghosh, Mandal & Das (2009) for reviews on the topic), but it is essentially limited to theoretical studies focused on the stability and persistency of the proposed flow configuration (Yih 1967;Hickox 1971;Hooper & Boyd 1983;Joseph et al. 1984), or to experimental studies measuring the overall performance and effectiveness of the targeted DR strategy (Oliemans et al. 1987;Bai, Chen & Joseph 1992;Arney et al. 1993;Bannwart 2001). Detailed experimental measurements of the flow field near the walls and near the liquid/liquid interface remain extremely challenging, in particular when optical techniques cannot easily be applied, as for example when opaque fluids (i.e. oils) are employed (Reinecke et al. 1998;Lindken, Gui & Merzkirch 1999;Heindel 2011). In this context, accurate simulations can be considered as a valuable alternative tool that, giving access to the entire velocity/stress field and to the corresponding liquid/liquid interface deformation, can be used -by dissecting the relevant flow phenomena occurring in the near-wall region and in the proximity of the interface -to fully characterize the underlying physics of the DR mechanisms. Not surprisingly, direct numerical simulations (DNS) of lubricated channels/pipelines have been performed more frequently in recent years (Ahmadi et al. 2018a,b;Kim & Choi 2018;Roccon et al. 2019), and have demonstrated the importance of both viscosity and surface tension on selectively modulating turbulence so as to give the observed DR. Reportedly, in these numerical studies, the flow is forced to move inside the channel/pipeline either by imposing a constant pressure gradient (CPG), or by imposing a constant flow rate (CFR). With the CPG approach, the mean pressure gradient that drives the flow is constantand it is usually set via the shear Reynold number, Re τ -while the volume flow rate can exhibit fluctuations that ultimately depend on the turbulent nature of the flow. In contrast, with the CFR approach, the flow rate is constant -usually set via the bulk Reynolds number, Re b -while the mean pressure gradient that pushes the flow is adapted at each time step so as to give the prescribed value of the volume flow rate. The effectiveness of the CPG or CFR approach in evaluating the DR performance is, however, questionable (Hasegawa, Quadrio & Frohnapfel 2014), and precisely because, by keeping the pressure gradient constant and increasing the total flow rate (respectively by keeping the flow rate constant and decreasing the pressure drop), the power injected into the flow -proportional to the product between the pressure gradient and the flow rate -increases (respectively decreases).
Motivated by this idea, we decided to reconsider the DR problem in lubricated channels by running a completely new set of finely resolved DNS using the so-called constant power input (CPI) approach (Hasegawa et al. 2014;Gatti et al. 2018). Within the CPI approach, the pressure gradient used to drive the flow is dynamically adjusted depending on the measured instantaneous flow rate so as to keep constant the overall power injected into the system. Thanks to this feature, the CPI approach paves the way for a more meaningful comparison among the different DR techniques, since the obtained results are no longer biased by the different powers injected into the system (Hasegawa et al. 2014).
In this work we focus on a simplified, yet practically relevant, configuration in which a thick layer of a primary fluid (also referred to as the primary layer hereinafter, and characterized by density ρ 2 , viscosity η 2 , thickness h 2 ) and a thin layer of a lubricating fluid (also referred to as lubricating layer hereinafter, and characterized by density ρ 1 , viscosity η 1 , thickness h 1 ) are forced to flow inside a channel lying one on top of the other, in a way such that the lubricating layer wets the top wall of the channel only. The primary and the lubricating fluids, which are immiscible, are characterized by the same density (ρ 1 = ρ 2 = ρ) but different viscosity, so that a viscosity ratio λ = η 1 /η 2 can be introduced. To cover a meaningful range of all possible situations, we consider five different values of λ in the range 0.25 ≤ λ ≤ 4, indicating that the lubricating fluid can be either less or more viscous than the primary fluid. We unambiguously show that a significant DR can be achieved. In particular, the observed DR is a non-monotonic function of λ, and appears to be maximized for λ = 1.00. However, and counter to intuition, we demonstrate that DR occurs not only for λ ≤ 1.00 -i.e. when the viscosity of the lubricating layer is smaller than that of the primary layer -but also for λ = 2.00 -i.e. when the viscosity of the lubricating layer is twice that of the primary layer. The situation changes on further increasing λ, and for λ = 4.00 we report a slight drag enhancement. These observations seem to suggest the presence of different turbulence modulation mechanisms, which we try to characterize by looking at the mean and turbulent kinetic energy budget, and at the corresponding energy fluxes. To do this, and precisely to evaluate the different energy fluxes -injected, transferred and dissipated by the systemwe adopt the energy-box representation (introduced by Ricco et al. (2012), Gatti et al. (2018) for single-phase flows), and we purposely extend it to the case of a liquid-liquid multiphase flow.
The paper is organized as follows: in § 2 we present the governing equations, the numerical method and the CPI approach; the main results of the simulations are shown in § 3: first, we compute and discuss the main global parameters of the flow, including the overall pressure gradient and the flow rates of the two fluid layers, as well as the corresponding flow structures; then, we study the mean and turbulent kinetic energy budget and we graphically explain the obtained results employing the modified energy-box representation. Finally, conclusions are outlined in § 4.

Methodology
We consider the case of two immiscible fluid layers flowing one on top of the other inside a rectangular flat channel. The channel has dimensions L x × L y × L z = 4πh × 2πh × 2h along the streamwise (x), spanwise (y) and wall-normal (z) directions. The bottom part of the channel is occupied by the primary fluid layer, having thickness h 2 , density ρ 2 and viscosity η 2 , while the top part of the channel is occupied by the thin lubricating layer, having thickness, h 1 , density ρ 1 and viscosity η 1 . The two layers have the same density (ρ 1 = ρ 2 = ρ) but different viscosity, so that a viscosity ratio λ = η 1 /η 2 can be defined. The interface separating the two phases is characterized by a constant and uniform value of the surface tension, σ . To capture the dynamics of the system, we couple direct numerical simulation (DNS) of the Navier-Stokes equations, used to describe the flow field, with a phase-field method (PFM), used to describe the interfacial motion (Jacqmin 1999;Badalassi, Ceniceros & Banerjee 2003;Roccon et al. 2017;Soligo, Roccon & Soldati 2019b).

Phase-field method
In the framework of the PFM, the sharp interface between the two fluids is replaced by a thin transition layer where the interfacial forces are applied. The basic idea of the PFM is to introduce an order parameter, the phase field φ, that is uniform in the bulk of the two phases (φ = 1 in the lubricating layer and φ = −1 in the primary layer) and that changes rapidly, yet smoothly, across the thin interfacial layer that separates the two phases. The transport of the phase field φ is described by a Cahn-Hilliard equation, which in dimensionless form reads as where u i is the ith component of the velocity vector, Pe Π is the Péclet number and μ is the chemical potential. The Péclet number is defined as follows: where u Π is the characteristic velocity (that will be properly introduced and discussed later, see (2.11) for details), h is the channel half-height, M is the mobility and β is a positive constant introduced to make the chemical potential dimensionless. The Péclet number identifies the ratio between the diffusive time scale, h 2 /Mβ, and the convective time scale, h/u Π of the interface. The chemical potential μ is defined as the functional derivative of a Ginzburg-Landau free-energy functional, the expression of which is chosen to represent an immiscible binary mixture of isothermal fluids (Soligo, Roccon & Soldati 2019a, 2020aSoligo et al. 2019b). The functional is composed by the sum of two different contributions: the first contribution, f 0 , accounts for the tendency of the system to separate into the two pure stable phases, while the second contribution, f mix (mixing energy), is a non-local term accounting for the energy stored at the interface. The mathematical expression of the functional is where Ω is the domain considered and Ch is the Cahn number, which represents the dimensionless thickness of the interface between the two fluids and is defined as where ξ is the physical thickness of the interface. From (2.3), we can obtain the expression of the chemical potential as (2.5) Note that, upon substitution of this expression into (2.1), a fourth-order equation for the phase field variable is obtained.

Hydrodynamics
To describe the hydrodynamics of the multiphase system, the Cahn-Hilliard equation is coupled with the Navier-Stokes equations. The presence of an interface (and of the corresponding surface tension forces) is accounted for by introducing an interfacial term in the Navier-Stokes equations. Recalling that in the present study we consider two fluids with the same density (ρ = ρ 1 = ρ 2 ) but different viscosity (η 1 / = η 2 ), continuity and Navier-Stokes equations become where u i is the ith component of the velocity vector, p is the pressure field, p x is the mean pressure gradient driving the flow, δ ij the Kronecker delta, η(φ) is the viscosity map accounting for the viscosity contrast between the two phases and defined as (Ahmadi et al. 2018a,b) Finally, τ c ij is the Korteweg tensor (Korteweg 1901), which is used to account for the surface tension forces, and is defined as follows: (2.9) The dimensionless groups appearing in the Navier-Stokes equations are the power Reynolds number, Re Π , and the Weber number, We Π , which are defined as The Reynolds number represents the ratio between inertial and viscous forces and is defined based on the viscosity of the primary layer η 2 , while the Weber number is the ratio between inertial and surface tension forces.
In contrast with our previous works (Ahmadi et al. 2018a,b;Roccon et al. 2019), in which we applied a constant mean pressure gradient to drive the flow through the channel (Lyons, Hanratty & McLaughlin 1991;Soldati & Banerjee 1998), in the present work we employ a CPI approach Hasegawa et al. 2014;Gatti et al. 2018 based on driving the flow by an imposed constant pumping power, P p . Naturally, to keep the pumping power constant over time, the mean pressure gradient is dynamically adjusted according with the overall flow rate, Q t (see appendices A and B for details).
Within the CPI approach, instead of using the commonly adopted friction velocity u τ (or the bulk velocity u b ), the following velocity is introduced as reference: where, as stated above, η 2 is the viscosity of the primary layer, while B and D are two coefficients used to account for the presence of a thin lubricating layer with different viscosity in the upper part of the channel (see details in appendix A). From a physical point of view, the characteristic velocity u Π represents the bulk velocity (average velocity across the channel section) of the actual two-phase flow configuration (i.e. two immiscible fluid layers having different viscosity and flowing inside a channel under the action of a pumping power P p ), but in laminar conditions. When the viscosity of the two layers is the same (single-phase case and λ = 1.00 case), the coefficients B and D are unitary and the characteristic velocity reduces to clearly matching the standard definition of the reference velocity under CPI conditions (Hasegawa et al. 2014).

Numerical method
The governing equations (2.1)-(2.6) and (2.7) are solved using a pseudo-spectral method based on transforming the field variables into wavenumber space via a combination of Fourier series (along the periodic streamwise and spanwise directions) and Chebyshev polynomials (along the inhomogeneous wall-normal direction). In particular, (2.7) is rewritten as a fourth-order equation for the wall-normal component of the velocity u z and a second-order equation for the wall-normal component of the vorticity ω z (Kim, Moin & Moser 1987;Speziale 1987) while (2.1) is split into two second-order equations (Badalassi et al. 2003).
The governing equations are advanced in time using a mixed implicit-explicit scheme. For the Navier-Stokes equations, the nonlinear diffusive term is first rewritten as the sum of a linear and a nonlinear contribution (Zonta, Marchioli & Soldati 2012a;Ahmadi et al. 2018a,b). The linear part is then integrated using a Crank-Nicolson scheme (second-order accurate) while the nonlinear part, together with the nonlinear convective terms, is integrated using an Adams-Bashforth scheme (second-order accurate). Similarly, for the Cahn-Hilliard equation, the linear term is integrated using an implicit Euler scheme (first-order accurate), while the nonlinear term is integrated in time using an Adams-Bashforth scheme (second-order accurate). The adoption of the implicit Euler scheme helps damping unphysical high-frequency oscillations that could arise from the steep gradients of φ. Note that the applied mean pressure gradient is also updated at each time step, so as to keep constant the power injected into the system (see appendix B). Further details on the code implementation, parallelization and validation can be found in Soligo et al. (2019bSoligo et al. ( , 2020a.

Boundary conditions
The resulting set of governing equations is complemented by suitable boundary conditions. For the Navier-Stokes equations, no-slip boundary conditions are enforced at the top and bottom walls (z/h = ±1) u i (z/h = ±1) = 0. (2.13) For the Cahn-Hilliard equation, no-flux boundary conditions are applied at the two walls, and finally give ∂φ ∂z Along the streamwise and spanwise directions (x and y), periodic boundary conditions are imposed for all variables (Fourier discretization). The adoption of these boundary conditions leads to the conservation of the phase field over time This gives the mass conservation of the entire system but does not guarantee the mass conservation of each of the two phases φ = +1 and φ = −1 (Yue, Zhou & Feng 2007;Soligo, Roccon & Soldati 2019c) and some small leakages between the phases may occur.
In the present case, mass leakage is always below 1 %.

Simulations set-up
We run six different simulations: a single-phase reference simulation and five simulations of lubricated channels, each characterized by a different value of the viscosity ratio λ. We consider both the case of λ < 1, i.e. the lubricating fluid is less viscous than the primary fluid, and of λ > 1, i.e. the lubricating fluid is more viscous than the primary fluid. In particular, we set: λ = 0.25, λ = 0.50, λ = 1.00, λ = 2.00 and λ = 4.00. All simulations are performed injecting into the system the same physical power P p (CPI approach). For the single-phase case and λ = 1.00 (uniform viscosity), this leads to a power Reynolds number equal to Re Π = 12220 (approximately corresponding to a shear Reynolds number Re τ 300). However, when the viscosity is not uniform (λ / = 1), the characteristic velocity u Π changes, leading to slightly different power Reynolds numbers: from Re Π = 14830 (λ = 0.25) down to Re Π = 11240 (λ = 4.00).
The surface tension value of the liquid-liquid interface is constant for all cases and is set via the Weber number so as to be representative of an oil/water interface (Than et al. 1988). The resulting Weber number is We Π = 830 for λ = 1.00, while it changes -because of the different characteristic velocity -for the other cases.
The grid resolution is chosen to fulfil requirements imposed by DNS and at the same time to guarantee a proper resolution of the thin interface between the two fluid layers. For the single-phase reference case, we use N x × N y × N z = 512 × 256 × 257 grid points, while for the lubricated channel cases we use N x × N y × N z = 1024 × 512 × 513. The Cahn number is set to Ch = 0.01 to allow for the accurate description of the steep gradients present at the interface (Soligo et al. 2019b). The Péclet number (or more specifically the mobility) is chosen following previous investigations (see for example Yue, Zhou & Feng 2010), which prescribe an optimal value Pe = 3/Ch to obtain an asymptotic convergence to the sharp-interface limit. The resulting Péclet number is Pe Π = 12 220 for λ = 1.00 and changes slightly for the other cases (different characteristic velocity). Please refer to table 1 for an overview of the simulation parameters.
For all simulations, the initial condition is taken from a preliminary DNS of a single-phase fully developed turbulent channel flow at Re τ = 300 (performed using a CPG approach), and complemented by a proper definition of the initial distribution of the phase field φ, so that the liquid-liquid interface is at the beginning flat and located at distance h 1 = 0.15h from the top wall. Specifically, the initial condition used for the phase field is (2.16) 2.6. CPI scaling and flow field decomposition In the following, results are presented using the CPI scaling system, which employs u Π as reference velocity, h as reference length and h/u Π as reference time. Unless otherwise mentioned, and for the sake of comparison, all results are normalized by the single-phase reference velocity, u sp Π . Angular brackets, · , indicate average in space along the two homogeneous directions (x and y), while square brackets, [·], indicate average in space (along x and y) and in time. The flow field, pressure field and Korteweg tensor are decomposed into a mean and a fluctuating component using a standard Reynolds decomposition (Reynolds 1895): the mean component is a function of the wall-normal coordinate (z) and time (t), while the fluctuating component depends on all the three spatial coordinates (x, y, z) and time (t). Therefore, a generic field, a (x, y, z, t), is decomposed as follows: (2.17) Note that the mean pressure gradient, used to drive the flow, has been already separated from the pressure field (see (2.7)).

Results
First, results will be discussed on a qualitative basis, by looking at instantaneous flow visualizations. Then, suitable flow statistics such as the flow rate of the two liquid layers, the mean pressure gradient and the corresponding mean velocity, will be used to quantify more closely the overall turbulence modulation in the lubricated channel.
To detail the different mechanisms that contribute more significantly to the observed flow behaviour, we will focus on the mean and turbulent kinetic energy budget, which we will also describe with the help of the so-called energy-box representation (Ricco

Qualitative behaviour of the multiphase flow
We start our discussion by looking at the qualitative structure of the flow in the statistically steady state that develops when the initial transient -required to absorb the initial conditions -is finished. Figure 1 shows a map of the TKE, TKE = u i u i /2, on a cross-section of the channel (y-z) located at x = L x /2. Each panel refers to a different case: single phase (a), λ = 0.25 (b), λ = 0.50 (c), λ = 1.00 (d), λ = 2.00 (e) and λ = 4.00 ( f ). The position of the interface (iso-level φ = 0) is indicated by a white line.
The qualitative results presented here confirm and extend our previous observations (Roccon et al. 2019). To discuss the flow modifications induced by the thin lubricating layer on the overall flow behaviour, we conveniently keep the single-phase case (figure 1a), for which classical near-wall turbulence structures are observed near the top and bottom boundaries, as reference. We immediately observe that, for all the viscosity stratified cases, the symmetry about the channel centre is broken, and the top and bottom boundaries behave differently based on the specific value of λ. We will focus on the bottom wall first, on the top wall later.
At the bottom wall (z/h = −1) the turbulence structure for all viscosity stratified cases (figure 1b-f ) is similar to that of the single-phase case (figure 1a). The thin lubricating layer is too far to influence the near-wall turbulence regeneration cycle at the bottom wall (Lam & Banerjee 1992;Jiménez 2013). Nevertheless, a slight enhancement of the turbulence activity is observed for λ = 0.25, λ = 0.50 and λ = 1.00. We anticipate here, but it will become clear later by looking at the fluid statistics, that this increase can be traced back to the larger gradients of the mean velocity profile and to the associated enhancement of the TKE production (Mansour, Kim & Moin 1988;Zonta et al. 2012a).
At the top wall the situation is completely different and depends on the specific value of λ. For ease of discussion, we focus first on λ = 1.00. Results for this case, which are given in figure 1(d), clearly show that turbulence is almost completely suppressed in the lubricating layer. This result, which confirms the observations made in our previous work (Roccon et al. 2019), can be directly attributed -as will be shown below -to the lubricating layer being too thin (or, alternatively to the local Reynolds number being too small) to sustain the near-wall turbulence cycle (Jiménez & Moin 1991;Jiménez & Pinelli 1999;Jiménez 2013;Roccon et al. 2019). From a physical viewpoint, turbulence suppression is in this case due to the presence of the liquid-liquid interface -an elasticity element with the corresponding surface tension -that limits the wall-normal vertical transport of momentum and decouples the dynamics of the primary layer from that of the lubricating layers, which is then too thin to sustain turbulence itself. Mild velocity fluctuations might appear occasionally, and are induced by the motion of the liquid-liquid interface. The turbulence suppression process observed in the lubricating layer for λ = 1.00 appears magnified for λ ≥ 1 (figure 1e, f ), since the effect of the larger viscosity of the lubricating layer sums up with the action of the surface tension forces (Roccon et al. 2017). By contrast, for λ = 0.25 and λ = 0.50 (figure 1b,c), the viscosity in the lubricating layer is small enough -and the corresponding local Reynolds number large enough (Pecnik & Patel 2017) -to sustain turbulence in the lubricating layer.
To appreciate better the different flow structures in the lubricating layer, we turn now our attention to the distribution of TKE on a x − y horizontal plane located in the proximity of the top wall, at distance d/h = 0.03. Results are shown in figure 2 and each panel refers to a different case: single phase (a), λ = 0.25 (b), λ = 0.50 (c) and λ = 1.00 (d). Note that results for λ = 2.00 and λ = 4.00 are not shown since velocity fluctuations are largely suppressed and are therefore not visible. For the single-phase case (figure 2a), we recover the classical near-wall turbulence structure in which elongated high and low speed streaks appear side by side and line up along the streamwise direction (Lam & Banerjee 1992;Schoppa & Hussain 2002;Chernyshenko & Baig 2005). The picture changes quite significantly for the lubricated channel. In particular, for λ = 1.00 (figure 2d), we see no evidence of turbulence activity, as turbulence is almost completely suppressed by the presence of the interface (Verschoof et al. 2016;Roccon et al. 2019). For λ = 0.25 and λ = 0.50 (figure 2b,c) turbulence appears reactivated in the lubricating layer. Yet, compared to the single-phase case, turbulence is less homogeneous, with the presence of banded structures of turbulent-laminar patches closely recalling those observed in other important flow instances (for example in thermally stratified turbulence, see García-Villalba & Del Álamo 2011;Zonta, Onorato & Soldati 2012b;Zonta & Soldati 2018). The coexistence of turbulent and laminar patches seems to be due to the deformation of the underneath interface separating the primary and the lubricating layer: crests and troughs (see appendix C for details) that naturally develop at the interface sometimes make the lubricating layer thin enough to locally suppress the turbulence activity. As a side observation -but we will not further comment on it -we report the change of size of turbulence structures at the wall, which appear smaller for smaller viscosity (i.e. larger local Reynolds number).
From a quantitative viewpoint, the presence of two distinct flow regimes (laminar/turbulent) in the lubricating layer can be justified using a semi-local scaling (Pecnik & Patel 2017;Roccon et al. 2019) and computing the thickness of the lubricating 911 A37-10 where h 1 is the lubricating layer thickness in outer units, Re τ is the equivalent shear Reynolds number (approximately 300 here) and τ w,1 and τ w,2 are the wall-shear stresses at the top and bottom wall, respectively. In the present case, h + 1 30 w.u. for λ = 1.00, while it increases up to h + 1 = 138 w.u. for λ = 0.25. These results suggest that the lubricating layer is too small to sustain turbulence for λ = 1.00, whereas it becomes large enough to sustain turbulence for λ < 1. Note that Jiménez & Pinelli (1999) established that the minimum dimensionless channel height to maintain a self-sustained near-wall turbulence cycle is h + 60 w.u.

Flow rates and pressure gradient
The qualitative flow changes observed above naturally reflect into corresponding changes of the macroscopic flow parameters. Here, we focus on the time-averaged flow rate through the entire cross-section, [Q t ], on the time-averaged flow rate of the primary layer, [Q 2 ], and on the time-averaged mean pressure gradient, [p x ]. These results, which are represented in figure 3, are normalized by the corresponding values for the single-phase case, [Q sp ] and [p sp x ], respectively. The region coloured in light grey in the diagram refers to flow configurations for which the viscosity of the lubricating layer is smaller than that of the primary layer (i.e. λ < 1), while the region coloured in dark grey refers to flow configurations for which the viscosity of the lubricating layer is larger than that of the primary layer (i.e. λ > 1).
Low-viscosity lubricating layer High-viscosity lubricating layer We immediately observe that, for λ ≤ 2.00, [Q 2 ]/[Q sp ] > 1 (empty circles) and . In other words, by keeping the applied power constant, the transferred flow rate is increased. This is a clear indication of DR, which is maximum for λ = 1.00 ( 13 % of flow-rate increase), and is only slightly reduced for λ = 0.50 and λ = 0.25 ( 11 % of flow-rate increase). The reduction appears more intense for λ = 2.00 ( 5 % of flow-rate increase). When the viscosity of the lubricating layer is increased beyond λ > 2, we observe a drop of the flow rate that, giving [Q 2 ]/[Q sp ] < 1 for λ = 4.00 ( 5 % of flow-rate decrease), marks the occurrence of a drag increase.
Consistently with the previous observations, for λ ≤ 2.00 the normalized pressure gradient is [p x ]/[p sp x ] < 1, so as to balance the increased flow rate and to keep the total power -product between the flow rate and the applied pressure gradient, in physical units -constant. The pressure gradient, which attains a minimum for λ = 1.00 x ] 0.9 for λ = 0.5 and λ = 0.25); on the other hand, it increases more vigorously for λ > 1, up to the point at which, for λ = 4.00, it becomes slightly larger than unity.
Note that, since simulations are performed at slightly different power Reynolds numbers (see table 1), the product between the total flow rate and the mean pressure gradient, which is by construction constant in physical units, is slightly different when considered in dimensionless units.

Mean velocity profiles
Linked to the previous analysis on the flow-rate and pressure gradient modulations in the lubricated channel, we focus now on the behaviour of the mean streamwise velocity profile, [u x ], as a function of the wall-normal coordinate. Results are shown in figure 4 according to the following colour code: red is used for λ = 0.25, yellow for λ = 0.50, green for λ = 1.00, light blue for λ = 2.00 and dark blue λ = 4.00. The single-phase case is represented by the thin black line, while the reference position of the interface (z/h = 0.85) is identified by the vertical dashed black line. As expected, the shape of the mean velocity profile appears skewed by the presence of the liquid-liquid interface.
In the primary layer, and consistently with the previous analysis on the overall flow rate, the behaviour of [u x ] is non-monotonic and is maximized for λ = 1.00. A reduction of the viscosity ratio from λ = 1.00 to λ = 0.50 and λ = 0.25 -which almost overlap -induces only a small decrease of [u x ]. On the other hand, an increase of λ beyond unity, λ > 1.00, produces stronger modification. For λ = 2.00, the mean velocity profile appears largely reduced -although in any case well above the single-phase limit -compared to λ = 1.00. Moving to λ = 4.00, we reach a condition at which the velocity profile is very close to that of the single-phase case in a large proportion of the primary layer, but suddenly drops in the proximity of the interface.
In the lubricating layer, 0.85 < z/h < 1.00, the trend of the mean velocity turns out to be monotonic with the viscosity ratio λ: [u x ] increases consistently with decreasing λ. A similar behaviour characterizes the mean velocity gradient that, being minimum for λ = 4.00 and maximum for λ = 0.25, reflects the different flow structure in the lubricating layer -with turbulence progressively suppressed for increasing λ (see also figures 1 and 2). To have a classical view of the velocity profiles, they must be rescaled with the local friction velocity (i.e. evaluated at the corresponding wall). This comparison is offered in appendix D.

Energy box: preliminary concepts
To provide a sound characterization of the energy fluxes in the lubricated channel, we decide to use the so-called energy-box technique (Ricco et al. 2012;Gatti et al. 2018). This technique, starting from the mean and TKE balance equations, provides a thorough and meaningful representation of the energy fluxes in the system. Note that, and precisely because of the adoption of a CPI approach, all simulations are performed driving the system with the same physical power input P p . Therefore, a direct comparison between the different cases is possible (Hasegawa et al. 2014).
To compute the energy budget, the total kinetic energy of the flow is decomposed into mean kinetic energy (MKE), MKE = 1 2 u i u i , and TKE, TKE = 1 2 u i u i . In the following, 911 A37-13 and building on top of the MKE and TKE definitions, we will present and discuss the energy-box technique for single and multiphase flows.

Energy box for a single-phase flow
The starting points of the energy-box technique are the MKE and TKE transport equations, obtained multiplying the Navier-Stokes equations by the mean velocity field -MKE balance equation -and by the fluctuating velocity field -TKE balance equation (Mansour et al. 1988;Kasagi, Tomita & Kuroda 1992;Iwamoto, Suzuki & Kasagi 2002). For the reference single-phase case, the ensemble average balance equation for the MKE becomes ( 3.2) where u i u j are the Reynolds stresses. While the left-hand side represents the material rate of change of MKE, the different terms on the right-hand side represent the power injected in the system via the mean pressure gradient, Π m , the production of TKE, P k , the work done by the Reynolds stresses, T m , the viscous diffusion of MKE, D m , and the mean flow viscous dissipation, m . The ensemble-averaged balance equation for the TKE reads as ( 3.3) The left-hand side represents the material rate of change of TKE. The terms on the right-hand side represent the pressure diffusion, Π k , the production of TKE, P k , the turbulent diffusion, T k , the viscous diffusion of TKE, D k , and the turbulent viscous dissipation, k . We can now proceed to discuss the different terms of each equation from a physical point of view. For (3.2) -MKE balance -the left-hand side is zero, since the flow is statistically steady. On the right-hand side, the power input, Π m , is a source term and represents the power injected in the system via the mean pressure gradient. This power is then partially dissipated by the mean flow viscous dissipation, m , and partially used to generate turbulent fluctuations via the production term, P k . Therefore, m and P k represent a sink of energy in the MKE balance equation. The two remaining terms, the energy transport by the Reynolds stress, T m , and the viscous diffusion, D m , redistribute the energy across the channel and thus act only as internal transport mechanisms that do not bear a net contribution. In other words, they are neither sinks nor sources.
Considering (3.3) -TKE balance -the left-hand side is also zero, always in light of the statistically steady condition. On the right-hand side, the TKE production term, P k , accounts for the generation of velocity fluctuations via mean shear. This term, which was a sink in the MKE equation, represents a source in the TKE equation. The power injected via P k is entirely dissipated by the turbulent dissipation, k , which is the only sink term in (3.3). The pressure diffusion, Π k , the turbulent diffusion, T k , and the viscous diffusion, D k , are redistribution terms with no net contribution.
To evaluate the overall importance of the different terms, we compute the integral of the two balance equations along the wall-normal direction. With reference to a generic term A m/k , and using an overbar to indicate the integral over the wall-normal direction z, we obtain (3.4) Clearly, this integral is zero for the internal transport (redistribution) terms, while it is positive (respectively negative) for the source (respectively sink) terms. Therefore, the integral counterparts of the MKE and TKE balance equations, (3.2) and (3.3), read as respectively. Combining (3.5) and (3.6), we obtain a balance equation for the total kinetic energy which states that the entire power injected into the system, Π m , is ultimately dissipated by viscosity, both via the viscous dissipation term of the mean field, m , and via the viscous dissipation term of the fluctuating field, k . We can now visualize the energy fluxes in the present configuration using the energy-box representation, as shown in figure 5. The two adjacent boxes represent the MKE (left) and TKE (right) content of the flow, respectively. Arrows are used to represent energy source (green), energy exchange (blue) and energy sink (red). The magnitude of each contribution, normalized by the power input, is given near the corresponding arrow. By looking at the left box, it is apparent that the applied power, Π m , which is injected into the mean flow (MKE), is in part ( 54 %) dissipated by the mean flow itself, m , and in part ( 46 %) transferred to TKE, P k . The power transferred to TKE is then dissipated by turbulent dissipation, k . When lumped together, the MKE and TKE boxes represent the energy balance of the total kinetic energy, (3.7).

Energy box for the lubricated channel
We now extend the energy-box representation to the case of a lubricated channel. As done before, to derive the modified energy balance equations, we multiply the Navier-Stokes equations by the mean and fluctuating velocity field and we average the resulting equations in space and time. For the MKE, we get where we can observe the appearance of an additional term, ψ m , which represents the work of the surface tension forces on the mean flow.
For the TKE, we get (3.9) Even in this case, there is an additional term, ψ k , which represents the work exchanged, via the surface tension forces, between the interface and the fluctuating field (Li & Jaberi 2009). To derive the macroscopic balance equations for the lubricate channel and hence to set the corresponding energy-box representation, we integrate over the wall-normal direction the obtained MKE/TKE balance equations. For the MKE, we have while for the TKE, we have At this point, it is worth further elaborating on the energetics associated with the presence of an interface -an elasticity element -inside the flow. Following Joseph (1976, p. 242), and under the assumption of constant and uniform surface tension, we can obtain an energy balance equation for the liquid-liquid interface (see appendix E for details) which states that an increase of the interfacial area (i.e. an increase of the interfacial energy) must correspond to a negative work of the surface tension forces (i.e. velocity and surface tension forces act in opposite directions). In statistically steady state conditions, the average value of the interfacial area remains constant and thus indicating that the total power -sum of mean and fluctuating contributionsabsorbed/released by the surface tension forces is null (Joseph 1976;Aris 1989;Dodd & Ferrante 2016). Equation (3.13) highlights the pure elastic behaviour of the interface and thus the absence of net energy dissipation associated with its deformation. Combining (3.13) with (3.10) and (3.11), we get the balance equation for the total kinetic energy Π m + m + k = 0. (3.14) Therefore, very much like the TKE production term, surface tension forces act as a transfer of energy between MKE and TKE. We anticipate here that the surface tension contribution is a sink term in the MKE equation, and a source term in the TKE equation. However, it must be also mentioned that as customary when a one-fluid approach is employed to analyse a multiphase system (Prosperetti & Tryggvason 2009;Popinet 2018), and precisely because the surface tension forces are smeared out on a thin interfacial layer, the two interfacial contributions do not exactly match and |ψ m | ≥ |ψ k |. Nevertheless, for the cases presented here, this difference is always smaller (in the worst condition) than 0.4 % of the total power injected in the system. For ease of notation, in the energy-box representation reported below, the surface tension contribution is conventionally identified by the value of ψ m . A mechanistic view of the role of the surface tension forces on the energy transfer between the main and the lubricating layer is given in appendix F. The energy-box representation for the lubricated channel is shown in figure 6. Each panel refers to a different viscosity ratio: λ = 0.25 (a), λ = 1.00 (b) and λ = 4.00 (c). The results for λ = 0.50 and λ = 2.00 are not shown since they represent only intermediate cases that do not add to the discussion. As set above, source terms are represented by green arrows, transport terms by blue arrows and sink terms by red arrows. The additional contribution of the interfacial terms is rendered by a yellow arrow connecting the MKE box to the TKE box via an additional box that represents the interface and its dynamicsan elastic component absorbing and releasing energy. We will discuss the case λ = 1.00 (figure 6b) first, the case λ = 0.25 (figure 6a) and λ = 4.00 (figure 6c) later. This strategy helps explicating the impact of surface tension and viscosity ratio on the energy fluxes in a more systematic way, focusing on the contribution of each single property one independently to each other. For λ = 1.00, and compared to the single-phase case, we observe an increase of the mean flow viscous dissipation, m , from 54 % up to 58 %. This result indicates that the reduced velocity gradients in the lubricating layer, which are due to the documented turbulence suppression therein (see § 3.1), are completely compensated by the increased velocity gradients at the bottom wall.
In addition, and precisely because of the observed turbulence suppression in the lubricating layer, we notice a significant decrease of TKE production, P k , which goes from 46 % (single phase) down to 39 % (λ = 1.00). The contribution of the surface tension forces, ψ m , is in this case rather small ( 2 %). The sum of the TKE production and of the interfacial contribution gives the amount of power used to generate TKE, which is ultimately dissipated via turbulent dissipation, k . Examined under this perspective, the DR turns out to be induced by the surface tension forces that, blocking the wall-normal transfer of momentum, promote the laminarization of the lubricating layer (which is characterized by a small thickness, hence a small local Reynolds number) and induce a sharp decrease of the overall turbulence production. As a consequence, the power injected in the system is more efficiently used to drive the motion of the two liquid layers.
Moving to λ = 0.25 (figure 6a), we now observe remarkable differences. The mean flow viscous dissipation, m , reduces to 48 %, while the TKE production increases up to 49 %. Accordingly, and considering that the contribution of the interfacial term, ψ m , remains rather small ( 3 %), we notice a corresponding increase of the turbulent dissipation up to 51 %. Present results, reporting a decrease of m and an increase of k , seem to be in contrast to the observed DR. Note indeed that, in the framework of the CPI approach and using passive DR techniques, DR is usually associated with an increase of m and a decrease of k (although, strictly speaking, this latter conclusion is based on the specific procedure used in Gatti et al. (2018) to model the wall-normal behaviour of the wall-shear stress). To reconcile present results with previous predictions, an extended version of the energy-box technique -as will be discussed in § 3.8 -must be introduced.
For λ = 4.00 (figure 6c), the energy box is similar to that observed for λ = 1.00 (figure 6a), with an increase of the mean flow dissipation up to 57 %, and a decrease of the TKE production down to 38 % -again due to the turbulence suppression in the lubricating layer. Interestingly, we also report a slight increase of the contribution of the interfacial term, ψ m , which becomes 5 %. This latter observation seems to be related to the dynamics of interfacial waves: when the viscosity of the lubricating layer is increased, the larger shear rate induces a larger deformation of the interface and this leads to a larger ψ m (see also appendix F). The increase of ψ m (in magnitude) is also reflected by the release of TKE near the interface: increasing the viscosity of the lubricating layer, the liquid-liquid interface becomes less and less compliant and behaves almost like a solid wavy interface that promotes the generation of velocity fluctuations (Breugem, Boersma & Uittenbogaard 2006;Rosti & Brandt 2017) (see also the higher TKE levels close to the interface -on the primary layer side -shown in figure 1f ). Overall, for λ = 4.00, the analysis of the energy budgets seems to indicate DR (increasing m and decreasing k ), which it is in fact not observed in practice (see § 3.2).  lubricating and primary layers, respectively. Each arrow refers to a different energy flux, whose magnitude is reported (near the arrow) normalized by the power input value. Compared to the previous energy boxes (ensemble averaged), an additional subscript is used to distinguish between the lubricating and primary layer contributions. The power injected in the system, Π m , is represented by a green arrow; the mean flow viscous dissipations, m,1 and m,2 , by red arrows (left), and the TKE production terms, P k,1 and P k,2 , by blue arrows. Finally, the turbulent dissipations, k,1 and k,2 , are represented with red arrows (right). Note the appearance of dark blue arrows, identifying energy fluxes exchanged between the primary and the lubricating layer, F m and F k .
Altogether, these results for λ / = 1 show that the original energy-box technique, when applied straightforwardly to a multiphase flow, leads to conclusions which seem to be in contrast with previous findings (Gatti et al. 2018). To overcome these limitations, the energy-box technique must be extended -and a phase discrimination procedure must be included -to account for the multiphase nature of the system. This is the objective of the next sections.

Virtually lubricated channel
Before proceeding with the derivation of the phase-averaged approach applied to the lubricated channel, and precisely to evaluate a representative reference case against which current results can be conveniently compared, we introduce the concept of a virtually lubricated channel. A virtually lubricated channel is a single-phase flow that is a posteriori split -using a virtual interface located at a distance 0.15h far from the top wall. This virtual interface has no surface tension and therefore does not have any effects on the velocity field, since the flow is obtained by the solution of the Navier-Stokes equations in a single-phase system with uniform density and viscosity. The energy-box representation for the virtually lubricated channel is given in figure 7. Since the different terms refer now to a specific subregion of the domain (lubricating or primary layer), an additional subscript is used to distinguish between the virtual lubricating layer -subscript 1 -and the virtual primary layer -subscript 2 -contributions. Note also the presence of dark blue arrows, which identify the energy fluxes (transport/redistribution terms, referred to as F m and F k ) exchanged between the two virtual fluid layers -redistribution terms are indeed null only when integrated over the entire volume.
By looking at figure 7 we notice that the viscous dissipation of MKE has a similar value in both layers, m,1 m,2 , while there is a sharp difference in the viscous dissipation of TKE, k,1 / = k,2 . This representation gives a clearcut idea of the energy that is transferred and dissipated in the region virtually corresponding to the primary layer of the channel, and favours the introduction of an energy transfer efficiency, which we conveniently define as This parameter represents the ratio between the mean energy dissipation in the primary layer, and the total energy dissipation, i.e. after deduction of the energy transferred to the lubricating layer, in the primary layer.
In the following, we will make specific reference to this configuration.
3.8. Phase-averaged energy box for viscosity stratified flow Based on the above discussion (see § 3.6), it is apparent that the original energy-box technique must be properly extended to analyse the case of a lubricated channel. Because of the multiphase nature of the flow, an averaging procedure applied to the entire volumesuch as the one employed above -does not provide a separation between the primary and the lubricating layer contributions and, more importantly, does not grant access to the evaluation of the energy fluxes between these two layers.
To investigate these aspects, we compute the energy budget for each layer in isolation (Dodd & Ferrante 2016;Rosti et al. 2019). The obtained results are then represented in an opportunely modified energy box. To derive the phase-averaged MKE/TKE balance equations, we use the phase-field variable φ and we define a local concentration, 0 ≤ c l ≤ 1, for each liquid layer where, as anticipated before, the subscript l = 1, 2 is used to identify the lubricating and primary layer local concentrations, respectively. Multiplying MKE and TKE balance equations by the two local concentrations, and averaging in space (x, y) and in time as explained in § 3.6, we obtain the energy balance equations for each liquid layer. In particular, the phase-averaged MKE for each layer l, [MKE] l -counterpart of (3.8)becomes ( 3.17) Similarly, the phase-averaged TKE for each liquid layer l, [TKE] l -counterpart of (3.9)becomes From a physical point of view, the meaning of the different MKE/TKE terms remains unchanged. However, the contributions refer now to a specific layer; for instance, P k,1 is the TKE production in the lubricating layer, while P k,2 is the TKE production in the primary layer. Upon integration of the MKE/TKE balance equations along the wall-normal direction, a set of macroscopic balance equations is obtained for each liquid layer. For the MKE, we get (3.20) in the lubricating and in the primary layer, respectively. We recall that, compared to (3.10), the integral of the work done by the Reynolds stress, T m,l , and the integral of the viscous diffusion, D m,l , is not anymore zero since these contributions are now evaluated only in a portion of the domain (Dodd & Ferrante 2016;Rosti et al. 2019). Combining (3.19) and (3.20), and taking into account that the contribution of the internal transport terms -which represent the energy fluxes between the two phases -cancels out, we have P k,1 + P k,2 P k from which, employing mass conservation c 1 + c 2 = 1, we recover the total MKE balance, (3.8).
In a similar fashion, for the TKE, we get (3.23) in the lubricating and in the primary layer, respectively. Even in this case, the integral of the pressure diffusion, Π k,l , turbulent diffusion, T k,l , and viscous diffusion, D m,l , is not anymore zero since it is evaluated only in a portion of the domain. Naturally, combining (3.22) and (3.23), these contributions cancel out and we have −P k,1 − P k,2 −P k which represents the phase-averaged counterpart of the global TKE balance, (3.9).
The resulting energy-box representation is given in figure 8. Each panel refers to a different case: λ = 0.25 (a), λ = 1.00 (b) and λ = 4.00 (c). As done for figure 6, the cases λ = 0.50 and λ = 2.00 are not shown since they are intermediate configurations that do not add to the discussion. For completeness, we briefly recall the convention used in this representation. The left dashed rectangle refers to the MKE flow content, while the right dashed rectangle refers to the TKE flow content. The surface tension contribution is represented by a yellow arrow connecting the MKE box to the TKE box via an additional box (interface). The power input, which is in part injected into the primary layer, and in part injected into the lubricating layer, is represented by green arrows. The mean flow and turbulent dissipations of each layer are represented by red arrows, while the TKE production terms are represented by blue arrows. Vertical arrows, coloured in dark blue, connect the lubricating and primary layer boxes and identify the mean and TKE fluxes between the two phases (Dodd & Ferrante 2016;Rosti et al. 2019). Specifically, for the MKE box (left), these arrows account for the work done by the Reynolds stress and the MKE viscous diffusion -lumped together into the quantity F m -while, for the TKE box (right), they account for the energy fluxes due to pressure, viscous and turbulent diffusion -lumped together into the quantity F k .
As was explicitly discussed at the end of § 3.6, the value of the energy fluxes in the primary and in the lubricating layer cannot be directly compared with that observed in the canonical single-phase flow, since the latter is characterized by a different volume. Yet, they can be conveniently compared to those of the virtually lubricated channel, the conceptual framework in which the virtual (primary and lubricating) layers are characterized by the same nominal volume as the actual layers in the lubricated channel.
We consider first the case λ = 1.00 (figure 8b). Not surprisingly, most of the injected energy goes into the primary layer ( 96 %) while only a small amount of it goes into the lubricating layer ( 4 %). Focusing on the primary layer, we observe that a larger proportion of the injected power is dissipated by the mean flow ( 41 %). As expected in case of sustained turbulence, production of TKE is also large ( 38 %) and is entirely dissipated by viscous dissipation ( 39 %). In the lubricating layer the situation changes completely. While the mean flow dissipation remains very important ( 17 %), the production of TKE drops down drastically ( 1 %), as does TKE viscous dissipation ( 2 %). These results are consistent with the turbulence suppression mechanism observed in the lubricating layer (see also figures 1d and 2d). The contribution of the interfacial terms is rather small ( 2 %), whereas the contribution of the energy redistribution terms (dark blue arrows), in particular for the MKE, is not negligible ( 8 %). These observations confirm our previous intuition that, for λ = 1.00 and thanks to the action of the surface tension forces, the lubricating layer becomes laminar and reduces the irreversible energy dissipation. As a consequence, a larger amount of energy is available for the transportation of the primary layer (surface tension-driven DR).
For the case λ = 0.25 (figure 8a), and compared to the case λ = 1.00, in the primary layer we do not observe significant differences in the energy transfer rates, which appear only slightly reduced. What does change significantly is the energy transfer in the lubricating layer, which is characterized by slightly smaller mean flow dissipation ( 13 %) but a definitely larger TKE production ( 14 %), and TKE viscous dissipation ( 14 %). This clearly indicates a sustained turbulence activity in the lubricating layer. As for the case λ = 1.00, the interfacial term is rather small ( 3 %), whereas the energy redistributionin particular for the MKE -becomes progressively significant ( 15 %). This latter observation, which is particularly important, suggests that a considerable fraction of the total energy is removed from the primary layer, and spent into the lubricating layer. In this case, DR decreases compared to λ = 1.00.   and λ = 4.00 (c). The left dashed box represents the MKE of the flow while the right one represents the TKE of the flow. Inside each dashed box, the top and bottom rectangles identify the lubricating and the primary layer, respectively. The power injected in the system, Π m , is represented by a green arrow. The mean flow viscous dissipations, m,1 and m,2 , and turbulent dissipations, k,1 and k,2 , are represented by red arrows, and the TKE production terms, P k,2 and P k,1 , by blue arrows. Finally, surface tension contribution, ψ m , is represented by a yellow arrow. The dark blue arrows linking the lubricating and the primary layer boxes (inside each dashed box) represent the energy fluxes of MKE and TKE exchanged between the two layers, F m and F k , respectively.
We turn now to the case λ = 4.00 (figure 8c). While in the primary layer the dynamics does not change much, and the amount of mean flow dissipation, TKE production and TKE dissipation are only slightly different compared to λ = 0.25, in the lubricating layer the situation changes drastically, becoming somehow similar to that observed for λ = 1.00. Indeed, and consistently with a situation in which turbulence is almost entirely suppressed, in the lubricating layer the mean flow dissipation is very large ( 25 %) while the TKE production and dissipation are very small ( 1 %) and ( 4 %), respectively. The interfacial terms are only slightly larger ( 5 %), whereas the energy redistribution -in particular for the MKE -becomes important ( 16 %). Despite the fact that the distribution of energy fluxes for λ = 4.00 is seemingly close to that for λ = 1.00, the final outcome of these two cases is completely different, with the case λ = 4.00 being characterized by drag increase rather than DR.
A physically sound interpretation of these results can only be given upon comparison with the virtually lubricated channel (see figure 7 and corresponding discussion). Recalling the transfer efficiency parameter defined before, H = m,2 /( m,2 + k,2 ), and applying it for λ = 1.00, λ = 0.25 and λ = 4.00, we obtain: H λ=1.00 /H sp = 1.11, H λ=0.25 /H sp = 1.06 and H λ=4.00 /H sp = 0.99. Written in these terms, previous results can be interpreted as follows. For λ = 1.00, and precisely because of the combined effect of the turbulence suppression in the lubricating layer, and of the small energy fraction transferred from the primary to the lubricating layer, the amount of energy that remains available in the primary layer is larger compared to the virtually lubricated channel. DR is therefore observed. For λ = 0.25, the lubricating layer becomes turbulent, and the associated TKE production and dissipation increase. However, the viscosity of the lubricating layer is smaller, and the corresponding energy transfer from the primary to the lubricating layer remains small enough so that DR can still be observed. For λ = 4.00, turbulence is suppressed in the lubricating layer, with corresponding increase of the mean dissipation and decrease of TKE production and dissipation. However, in this case, the energy extracted from the primary layer to help driving the motion of the more viscous lubricating layer becomes important to such an extent that a slight drag increase is observed, and the overall energy loss is marginally larger than that observed for the single-phase case.
3.9. Theoretical prediction of the DR performance Before moving to the conclusions, we would like to present a simplified theoretical approach that can be used to predict the effects of the control parameters -lubricating layer thickness, Reynolds number and viscosity ratio -on the DR performance. To this aim, we can consider a simplified version of (3.1), which gives an estimate of the lubricating layer thickness in wall units where h 1 is the lubricating layer thickness in outer units, Re τ the equivalent shear Reynolds number and λ the viscosity ratio. Comparing the value obtained from (3.25) with the minimum value require to generate a self-sustained near-wall turbulence cycle (h + 1 60 w.u., Jiménez & Moin 1991;Jiménez & Pinelli 1999), we can predict the flow regime in the lubricating layer and thus the DR performance.
For λ < 1 (viscosity-driven DR), increasing h 1 or the Reynolds number (or both), the lubricating layer remains large enough (h + 1 > 60 w.u.) to sustain turbulence. As a result, the main driving mechanism for the DR is the viscosity difference between the two fluids, and large modifications of the DR are not expected. It must be also noted that the larger inertial forces that characterize the flow at larger Reynolds numbers can lead to the breakage of the thin lubricating layer and to a consequent loss of the DR performance. By contrast, decreasing h 1 or the Reynolds number (Ahmadi et al. 2018a,b), we expect a partial laminarization of the lubricating layer and thus an increase of the lubricating layer velocity. This can potentially lead to a further improvement of the DR performance.
For λ ≥ 1 (surface tension-driven DR), increasing h 1 or the Reynolds number (or both), the lubricating layer could at some point become turbulent. The transition to turbulenceor the presence of turbulence patches -will be more pronounced for λ = 1.00 and less pronounced for the larger viscosity ratios. Overall, this can lead to a weakening of the surface tension-driven DR mechanism and thus to a reduction of the DR performance. Decreasing h 1 or the Reynolds number, we do not expect large modifications of the DR performance since the flow in the lubricating layer will remain laminar.
Finally, note that decreasing the thickness of the lubricating layer -which is in principle a good strategy to improve the DR performances -might have some drawbacks, since the waves generated at the liquid/liquid interface could reach the top wall, thereby destroying the lubricating layer configuration and the possibility of having DR.

Conclusions
In this work, we have used DNS to analyse the problem of DR in a lubricated channel, a flow configuration in which a thin layer of a lubricating fluid is injected in the near-wall region of a plane channel so as to favour the transportation of a primary fluid.
To examine a wide range of possible situations, we have considered not only the case of a lubricating layer that is less viscous than the primary fluid (λ < 1), but also the case of a more viscous lubricating layer (λ > 1). In particular, we have tested five different viscosity ratios, from λ = 0.25 up to λ = 4.00. A PFM has been used to describe the dynamics of the interface, which is characterized by a uniform value of the surface tension. An important aspect of the present study is the use of the CPI approach to perform the simulations. This implies that the channel flow rate is modified depending on the actual value of the pressure gradient, in a way that the total power injected into the flow is kept constant. The CPI approach offers a well-defined theoretical framework for the analysis of the different DR techniques since the obtained results are not anymore influenced by the amount of power injected into the flow. Originally developed for single-phase flows, the CPI approach has been extended here for the first time to the case of multiphase flows.
Our results have clearly shown that a significant DR can be obtained in a lubricated channel flow. The observed DR is a non-monotonic function of λ, and reaches a maximum for λ = 1 -a case for which a flow-rate increase of 13 % is observed. This result directly points to the crucial role played by the surface tension in the DR process. The surface tension -an elasticity element that hinders the exchange of momentum between the primary and the lubricating layer -decouples the wall-normal momentum transfer mechanisms between the primary and the lubricating layer, suppresses turbulence in the lubricating layer (laminarization process) and naturally reduces the overall drag (surface-tension-driven DR). It is well known in the literature that the introduction of an elasticity factor inside a flow -for example, polymers or wall elasticity -can induce DR (White & Mungal 2008). However, and differently from previous investigations, in the present work, the elasticity is concentrated in a deformable surface inside the channel. When the viscosity of the lubricating layer is reduced (λ ≤ 1.00), turbulence is sustained in the lubricating layer, because of the increased local Reynolds number. In this case, DR is attributed to the smaller viscosity of the lubricating layer that directly decreases the corresponding wall friction (viscosity-driven DR). Interestingly, we have shown that DR can be obtained also for λ = 2.00 -i.e. when the viscosity of the lubricating fluid is twice the size of that of the primary fluid. In this case, DR is the outcome of the combined surface tension/large viscosity effect, which, reducing the vertical momentum and leading to a full relaminarization of the lubricating layer, completely balances the enhanced friction induced by the larger viscosity. We have also demonstrated that an upper limit for λ exists and, for λ = 4.00 we report a slight drag enhancement.
A clearcut explanation of all present results have been offered by looking at the mean and TKE budgets, which has been graphically rendered with the help of the energy-box representation. The energy-box representation has been presented and applied to the single-phase case first, and extended to the multiphase case later (phase-average technique). Thanks to this approach, we have been able to characterize the energy fluxes exchanged between the two phases, and to explain from a quantitative viewpoint the observed DR. In particular, introducing the energy transfer efficiency parameter, Hwhich quantifies the ratio between the mean and the total energy dissipation in the primary layer -we have been able to show that DR is observed when H/H sp > 1 (with H sp the energy transfer efficiency for the single phase flow), whereas drag enhancement is expected when H/H sp < 1.
Funding. Vienna Scientific Cluster (Vienna, Austria) and CINECA supercomputing center (Bologna, Italy) are gratefully acknowledged for generous allowance of computer resources under grants 71026, HP10BCFP82 and HP10BAX0FY.
The authors acknowledge the TU Wien University Library for financial support through its Open Access Funding Program.
Declaration of interest. The authors report no conflict of interest.
Appendix A. Extension of the CPI approach to viscosity stratified multiphase flows We present here the extension of the CPI approach to viscosity stratified flows. We start by briefly recapping its derivation for a single-phase flow. In particular, we consider the flow inside a plane channel bounded by two walls located at z = ±h, figure 9(a). The fluid is characterized by uniform viscosity η. A constant pumping power per unit area, P p , is used to drive the flow along the streamwise direction. Assuming the flow laminar, the following velocity profile is obtained: where (p x ) is the pressure gradient along the streamwise direction. The bulk velocity (i.e. the average velocity across the channel section) can be computed integrating the velocity profile along the wall-normal direction Viscosity η 1 Viscosity η 2 Figure 9. Sketch of the channel configurations used to derive the laminar flow solution: single phase (a) and viscosity stratified (b). For the single phase, the viscosity, η, is uniform and thus a symmetric velocity profile, u x , is obtained. For the viscosity stratified case, the thin lubricating layer (mh < z < h) has viscosity η 1 while the primary layer (−h < z < mh) has viscosity η 2 and an asymmetric velocity profile (u x,1 in the lubricating layer and u x,2 in the primary layer) is obtained. For both panels, the maximum velocity, u m , and the bulk velocity, u b , have been highlighted.
The power dissipated by the viscous forces can be computed as where the expression of the bulk velocity has been highlighted. At equilibrium (i.e. fully developed flow), the power dissipated by the viscous forces is equal to the power injected in the system. Matching the expression of P p with (A3), we can identify a velocity scale, u Π , for the problem To extend the CPI approach to viscosity stratified flows, we can proceed in a similar way. We consider now a fully developed channel flow in which two liquid layers flow one on top of the other, as sketched in figure 9(b). The top layer is characterized by viscosity η 1 while the bottom (primary) layer is characterized by viscosity η 2 . The flow is laminar and the interface between the two layers is located at z = mh (note that m = 0.85 in the present work). No-slip boundary conditions are applied at the two walls, while continuity of velocity and stress is enforced at the interface (Leal 1992, p. 198). The resulting velocity profile is where u x,1 is the velocity in the thin lubricating layer and u x,2 is the velocity in the primary layer. These two velocities are defined as follows: As previously done for the single-phase flow case, we can highlight the expression of the bulk velocity Consistently with the previous derivation, the bulk velocity is used as reference velocity scale and its expression can be obtained by matching the expression of the viscous dissipation (equation A18) with that of the pumping power. This gives Under CPI conditions, the coefficient B 2 /D is used to account for the presence of a lubricating layer of different viscosity near the top wall (equation A4).  Figure 11. Mean velocity profiles at the bottom wall (a) and at the top wall (b) rescaled in wall units using the local friction scale at the corresponding wall. The classical law of the wall: u + = z + and u + = (1/k) log(z + ) + 5 (where k = 0.41 is the von Kármán constant) is also reported as a reference. At the bottom wall, all profiles collapse on the classical law of the wall while at the top wall, for λ > 1, the velocity profiles depart significantly from the law of the wall and gradually approach the laminar behavior.
going from λ = 1 to λ = 4, the PDF widens. This trend seems to be due to the larger shear rate observed for λ ≥ 2.00, which acts to deform more the interface. Indeed, for λ ≥ 2.00, the interfacial contribution in the energy box (ψ m ), whose magnitude depends on the combined effect of the interface shape and mean velocity profile, is larger.

Appendix D. Rescaling of the mean velocity profiles in wall units
The mean velocity profiles reported in figure 4 using the CPI scaling system, can be also represented in wall units using a semi-local scaling (Pecnik & Patel 2017;Roccon et al. 2019) that relies on the local value of the wall-shear stress and viscosity (local friction velocity). The resulting friction velocity is then used to rescale the different mean velocity profiles, see figure 11. Panel (a) refers to the bottom wall (primary layer) while panel (b) refers to the top wall (lubricating layer). The classical law of the wall, u + = z + and u + = (1/k) log(z + ) + 5, with k = 0.41 the von Kármán constant (Von Kármán 1931), is also shown by a thin black line for comparison purposes.
At the bottom wall (figure 11a), all the velocity profiles (for both the single-phase and the multiphase case) collapse on the classical law of the wall. This indicates that the presence of the interface induces only negligible effects on the near-wall turbulence cycle in the lower part of the channel.
At the top wall (figure 11b), the mean velocity profiles are remarkably different. For λ < 1.00, the velocity profiles match the laminar behaviour u + = z + only very close the wall, but they seem to approach soon a logarithmic behaviour (although with a different steepness). For λ ≥ 1.00, the velocity profiles depart significantly from the law of the wall and gradually collapse onto the laminar behaviour u + = z + . Please note that, due to the rescaling (which is based on the local friction scale), the nominal interface position -which also corresponds to the lubricating layer thickness expressed in wall units h + 1 -is different among the various cases. In particular, the nominal interface position for the different cases is: z + 138 w.u. (λ = 0.25), z + 70 w.u. (λ = 0.50), z + 30 w.u. (λ = 1.00), z + 19 w.u. (λ = 2.00) and z + 12 w.u. (λ = 4.00).

Appendix E. Energy balance equation for surface tension
To derive the energy balance equation for the surface tension, we follow Joseph (1976, p. 242). Under the assumption of constant and uniform surface tension -namely, dσ/dt = 0 and ∇ s σ = 0, with ∇ s the surface gradient operator (see Slattery, Sagis & Oh 2007, p. 624) -the balance equation at the liquid-liquid interface can be written as where A is the interfacial area, u i is the ith component of the velocity vector and f σ i is the ith component of the surface tension forces. Since the interface has a finite, yet very small, thickness, the right-hand side of (E1) can be rewritten as a volume integral Recalling the mean value theorem for integrals, the right-hand side of (E2) can be rewritten as  This expression can be further simplified decomposing also the Korteweg tensor into a mean and a fluctuating component where, thanks to the properties of the average operator, the cross-terms vanish. Using the definition introduced in (3.4), we obtain

Appendix F. Physical interpretation of the interfacial contribution
To obtain a physical interpretation of the interfacial contribution, ψ m , we can consider the wall-normal behaviour of the mean velocity profile and of the mean component of the surface tension forces. Indeed, once integrated along the wall-normal direction, the product between these two quantities represents the interfacial contribution in the MKE balance equation. In figure 12, we report the qualitative behaviour of the mean velocity profile (left side), the mean component of the surface tension forces (centre) and a sketch of the interfacial waves (right side) and of the corresponding surface tension forces (red arrows). In addition, the nominal interface position (z/h = 0.85) is reported with a dashed horizontal line. As can be appreciated from the sketch, the waves are not symmetric (with respect to the nominal interface position axis) and the resulting contribution of the surface tension forces (i.e. the mean component of the surface tension forces) is positive above the nominal interface position and negative below it. These two contributions (positive and negative), which have the same magnitude, are then multiplied by the mean velocity and integrated along the wall-normal direction to compute the interfacial contribution ψ m . Since the velocity is larger below the nominal interface position than above it, the negative part dominates and we obtain ψ m < 0. From an energetic point of view, a negative value for the interfacial term identifies a situation in which surface tension forces absorb power from the mean flow. A corresponding amount of power is necessarily released in the flow in the form of turbulence fluctuations (TKE) by the term ψ k so that the overall power absorbed/released by the interface is zero (Joseph 1976;Aris 1989;Dodd & Ferrante 2016).
Although the present discussion is only qualitative, the very same conclusions can be obtained using a more quantitative approach.