Influence of density and viscosity on deformation, breakage, and coalescence of bubbles in turbulence

We investigate the effect of density and viscosity differences on a swarm of large and deformable bubbles dispersed in a turbulent channel flow. For a given shear Reynolds number, Re=300, and a constant bubble volume fraction, Phi=5.4%, we perform a campaign of direct numerical simulations of turbulence coupled with a phase-field method accounting for interfacial phenomena. For each simulation, we vary the Weber number (We, ratio of inertial to surface tension forces), the density ratio (r, ratio of bubble density to carrier flow density) and the viscosity ratio (e, ratio of bubble viscosity to carrier flow viscosity). Specifically, we consider two Weber numbers, We=1.50 and We=3.00, four density ratios, from r=1 down to r=0.001, and five viscosity ratios, from e=0.01 up to e=100. Our results show that density differences have a negligible effect on breakage and coalescence phenomena, while a much stronger effect is observed when changing the viscosity of the two phases. Increasing the bubble viscosity with respect to the carrier fluid viscosity damps turbulence fluctuations, makes the bubble more rigid, and strongly prevents large deformations, thus reducing the number of breakage events. Local deformations of the interface, on the contrary, depend on both density and viscosity ratios. The opposite effect is observed for increasing bubble viscosities. We report that these effects are mostly visible for larger Weber numbers, where surface forces are weaker. Finally, we characterize the flow inside the bubbles; as the bubble density is increased, we observe, as expected, an increase in the turbulent kinetic energy (TKE) inside the bubble, while as the bubble viscosity is increased, we observe a mild reduction of the TKE inside the bubble and a strong suppression of turbulence.

We numerically investigate the effect of density and viscosity differences on a swarm of large and deformable bubbles dispersed in a turbulent channel flow.For a given shear Reynolds number, Reτ = 300, and a constant bubble volume fraction, Φ ≃ 5.4%, we perform a campaign of direct numerical simulations (DNS) of turbulence coupled with a phase-field method (PFM) accounting for interfacial phenomena.For each simulation, we vary the Weber number (W e, ratio of inertial to surface tension forces), the density ratio (ρr, ratio of bubble density to carrier flow density) and the viscosity ratio (ηr, ratio of bubble viscosity to carrier flow viscosity).Specifically, we consider two Weber numbers, W e = 1.50 and W e = 3.00, four density ratios, from ρr = 1 down to ρr = 0.001 and five viscosity ratios from ηr = 0.01 up to ηr = 100.Our results show that density differences have a negligible effect on breakage and coalescence phenomena, while a much stronger effect is observed when changing the viscosity of the two phases.Increasing the bubble viscosity with respect to the carrier fluid viscosity damps turbulence fluctuations, makes the bubble more rigid and strongly prevents large deformations, thus reducing the number of breakage events.Local deformations of the interface, on the contrary, depend on both density and viscosity ratios: as the bubble density is increased, a larger number of small-scale deformations, small dimples and bumps, appear on the interface of the bubble.The opposite effect is observed for increasing bubble viscosities: the interface of the bubbles become smoother.We report that these effects are mostly visible for larger Weber numbers, where surface forces are weaker.Finally, we characterize the flow inside the bubbles; as the bubble density is increased, we observe, as expected, an increase in the turbulent kinetic energy (TKE) inside the bubble, while as the bubble viscosity is increased, we observe a mild reduction of the TKE inside the bubble and a strong suppression of turbulence.

I. INTRODUCTION
Interactions among turbulence and deformable interfaces are common in many physical instances, from ocean waves formation [1,2] to atomization processes [3], as well as drops and bubbles entrained in a turbulent flow [4][5][6].The outcome of these interactions is of fundamental importance as it controls the exchanges of heat, mass, and momentum across the interface and thus between the two phases.The study of turbulence-interface interactions, however, is a non-trivial task as these interactions are governed by a physics acting at very different spatio-temporal scales: from the largest problem scale, down to the Kolmogorov scale of turbulence and further down to the molecular scale of the interface.This multi-scale nature makes the investigation of multiphase turbulence very challenging.In particular, experimental investigations using optical techniques are usually limited to small volume fractions due to the difficulty of accessing phases with heterogeneous optical properties [7][8][9] and limited range of length scales that can be possibly measured.In this scenario, despite some limitations, numerical simulations represent an essential tool to investigate multiphase flows as they allow to access detailed space-and time-resolved information on the flow field and dispersed phase.Specifically, direct numerical simulation (DNS), in which all relevant scales of turbulence are resolved, proved to be a tool of paramount importance for a deeper understanding of single-phase [10,11] and multiphase turbulence [4][5][6].
In this work we focus on the interactions of a swarm of large and deformable bubbles or drops (bubbles hereinafter without any loss of generality) with wall-bounded turbulence (turbulent channel flow).This setup has been widely used in the past to investigate different aspects of bubbly flows, from bubbles shape, deformation and clustering to the flow modifications produced by the bubbles themselves.In the pioneering works of Lu & Tryggvason [12,13], the effects of the bubble size and deformability were investigated: they observed that as bubbles become more deformable, they move towards the middle of the channel and have a relatively small effect on the flow-rate.Scarbolo et al. [14,15], considering a matched density and viscosity system, investigated the effect of the surface tension, observing that surface tension forces play a key role in determining the dispersed phase topology.Roccon et al. [16] studied the effect of the bubble viscosity, finding that for small surface tension values, larger internal viscosities reduce the drop deformability.Recently, Soligo et al. [17,18], considering also the presence of a soluble surfactant, investigated the surfactant effects on drop morphology [17] and flow behavior [18].Finally, Hasslberger et al. [19] analyzed the coherent structures obtained in a bubble-laden turbulent channel flow while Cannon et al. [20] investigated the role played by droplets coalescence on drag in turbulent channel flows.
The foremost goal of this paper is to improve the fundamental understanding of bubble-bubble and bubbleturbulence interactions.Indeed, bubbles transported by a turbulent flow are characterized by complex dynamics, as they collide, coalesce and break apart.This behavior is governed by the forces generated by the surrounding continuous phase, acting on the surface of the bubbles with shear and normal stresses, and by the response of bubbles, which depends on their surface tension and their density and viscosity.The ultimate competition among these forces determines the number, shape, and deformation of the bubbles.In this work, we want to extend our previous works [14,16] and provide a comprehensive analysis on the effects of density ratio (ratio between the density of the bubble phase over the dispersed phase), viscosity ratio (ratio between the viscosity of the bubble phase over the dispersed phase) and surface tension (controlled by the Weber number, ratio of inertial over surface tension contributions) on the multiphase system.Specifically, the first objective is to investigate the effects of these parameters on the dispersed phase topology and its topological modifications (coalescence and breakage events), and to characterize the shape and deformation of the bubbles.The second objective of this work is to characterize the global and local flow modifications produced by bubbles on the turbulent channel flow behavior.To this aim, we build and analyze a database of direct numerical simulations of turbulent channel flows laden with deformable bubbles, considering different values of density ratios, viscosity ratios, and surface tension.The numerical framework of the simulations relies on a direct solution of the Navier-Stokes equations coupled with a phase-field method.Direct solutions of the Navier-Stokes equations are used to accurately resolve all the relevant turbulence scales, while the phase-field method [21,22] -an interface capturing approach that relies on an order parameter to define the local concentration of each phase -is used to describe in a thermodynamically consistent manner the motion of the deformable interface and its topological modifications (i.e.coalescence and breakage events).
The paper is organized as follows: in section II, we introduce the numerical method, the simulation setup, and the parameters of the simulations.Then, in section III, we present the results obtained from the analysis of the simulations database.First, we focus on the effects of density and viscosity ratios and surface tension values on the topology of the dispersed phase and its topological changes (breakage and coalescence).Secondly, we evaluate the effects of these parameters on the overall interfacial area and curvature of the bubble interface.Thirdly, we study the effects of density and viscosity ratios and Weber number on the mean velocity profiles and on the turbulent kinetic energy (TKE) of the bubbles.Finally, we summarize the results and draw our conclusions in section IV.

II. METHODOLOGY
We consider the case of a swarm of bubbles injected in a turbulent channel flow with a rectangular cross-section.The dispersed and carrier phases are characterized by density ρ d and ρ c , and viscosity η d and η c , where the subscripts d and c identify the dispersed and carrier phase, respectively.We define the density ratio and viscosity ratio as ρ r = ρ d /ρ c and η r = η d /η c respectively.The interface that separates the two phases is characterized by a constant and uniform value of the surface tension, σ.To describe the dynamics of the system, direct numerical simulation (DNS) of the Navier-Stokes equations, used to describe the flow field, are coupled with a phase-field method (PFM), used to describe interfacial phenomena [21,22].

A. Modeling of interfacial phenomena
The phase-field method uses an order parameter, the phase field ϕ, to identify the two phases: the order parameter is uniform in the bulk of each phase (ϕ = ±1) and undergoes a smooth transition across the interface.Indeed, the sharp interface is replaced by a thin transition layer.The transport of the phase field ϕ is described by the Cahn-Hilliard equation, which in dimensionless form reads as: where u = (u, v, w) is the velocity vector, P e is the Péclet number, µ is the chemical potential and f p is the penalty flux introduced with the profile-corrected formulation of the phase-field method [23][24][25].The Péclet number is defined as follows: where u τ = τ w /ρ c is the friction velocity (being τ w the shear stress at the wall and ρ c the carrier phase density), h is the channel half-height, M is the mobility parameter 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 variational derivative of a Ginzburg-Landau free-energy functional, the expression of which is selected to represent an immiscible binary mixture of isothermal fluids [17,25,26].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 , is a mixing 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 thin interfacial layer between the two fluids (ξ is the physical thickness of the interface).
From equation (3), the expression of the chemical potential can be derived as the functional derivative with respect to the order parameter: At the equilibrium, the chemical potential will be constant throughout all the domain.The equilibrium profile for a flat interface can thus be obtained solving ∇µ = 0, hence obtaining: where s is a coordinate normal to the interface.Finally, f p is the penalty-flux employed in the profile-corrected formulation of the phase-field method.This formulation is an improvement to the standard phase-field formulation: it allows to better maintain the equilibrium interfacial profile and it overcomes the drawbacks of the method (e.g.mass leakages among the phases and misrepresentation of the interfacial profile [23,27]).This penalty flux is defined as: where the numerical parameter λ can be set via the scaling λ = 0.0625/Ch [25].Before proceeding, it is worth to briefly discuss the main capabilities and limitations of interface-resolved simulations in describing topological modifications of the interface [6,17,28].The numerical description of breakages and coalescences is indeed one of the most challenging aspects of interface-resolved simulation methods.A fully-resolved simulation of topological changes would require resolving all the scales, from the molecular scale of the interface [29] up to the largest scales of the flow.This type of simulation, however, is way beyond the capabilities of any existing supercomputing facility.The common choice is to avoid resolving the small interfacial scales and to find a way to approximate their dynamics on a much larger scale.Here, following a similar approach, we limit the resolved range to the scales of turbulence: from the Kolmogorov length scale up to the problem size.Thus, phenomena occurring at scales smaller than Kolmogorov are smeared out on the smallest resolved scale.This choice however influences the description of coalescence and breakage events.For coalescences, a part of the physics involved in the coalescence process [30] (i.e.film drainage and rupture) cannot be directly resolved.As a result, regardless of the approach employed to describe coalescence (models for interface tracking methods [31,32] or implicit description for interface capturing methods [6,33]), numerical simulations struggle in predicting physical coalescence, with this inaccuracy referred to as numerical coalescence.For breakages, the picture is different and their numerical description is less troublesome.Indeed, being breakage a very quick phenomenon, it can be well approximated without resolving the dynamics at the molecular scale and there is evidence that the Navier-Stokes equations alone provide an adequate description of a breakage event [34].Besides, the small time scale of a breakage limits the impact of the approximation on the overall flow dynamics [32,35].Therefore, the description of breakages on turbulence-resolved grids is considered to be rather accurate, although in the pinch-off region the smallest interfacial features, characterized by high curvature, may not be perfectly resolved.

B. Hydrodynamics
To describe the hydrodynamics of the multiphase system, the Cahn-Hilliard equation is coupled with the Navier-Stokes equations.The presence of a deformable 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 having different densities and viscosities, we use here the formulation of continuity and Navier-Stokes equations proposed by Dong & Shen [36].The resulting governing equations for the hydrodynamics read as follow: where u = (u, v, w) is the velocity vector, p is the pressure field, T c is the Korteweg tensor and ρ(ϕ) and η(ϕ) are the density and viscosity fields, respectively.The density and viscosity fields are dimensionless scalar functions that account for the local value of density and viscosity respectively [37][38][39]; the carrier phase properties are used to make these fields dimensionless.The local density and viscosity are assumed to be linear functions of the phase field: where ρ r and η r are the density and viscosity ratios, respectively.
The Korteweg tensor [40], used to account for the surface tension forces, is defined as follows: The dimensionless groups appearing in the Navier-Stokes equations are the shear Reynolds number, Re τ , and the Weber number, W e, which are defined as: The Reynolds number represents the ratio between inertial and viscous forces, while the Weber number is the ratio between inertial and surface tension forces.Both Reynolds and Weber numbers are defined using the carrier phase properties (ρ c and η c ).

C. Numerical method
The governing equations ( 1), ( 8) and ( 9) are solved using a pseudo-spectral method, which uses Fourier series along the periodic directions (streamwise and spanwise) and Chebyshev polynomials along the wall-normal direction.The Navier-Stokes and continuity equations are solved using the velocity-vorticity formulation: equation ( 9) is rewritten as a 4 th order equation for the wall-normal component of the velocity u z and a 2 nd order equation for the wallnormal component of the vorticity ω z [11,41].Equation ( 1) is also split into two 2 nd order equations [22]; this way the governing equations are recasted as a coupled system of Helmholtz equations, which can be readily solved.The governing equations are time advanced using an implicit-explicit scheme.For the Navier-Stokes equations, the non-linear term is first rewritten as the sum of a linear and a non-linear contribution [42].Then, the linear part is integrated using a Crank-Nicolson implicit scheme, while the non-linear part is integrated using an Adams-Bashforth explicit scheme.Likewise, for the Cahn-Hilliard equation, the linear term is integrated using an implicit Euler scheme, while the non-linear term is integrated in time using an Adams-Bashforth scheme.The adoption of the implicit Euler scheme helps damping unphysical high-frequency oscillations that could arise from the steep gradients of ϕ.

D. 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 wall (z/h = ±1): For the Cahn-Hilliard equation, no-flux boundary conditions are applied at the two walls, yelding the following boundary conditions: 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 enforces mass conservation of the entire system but does not guarantee the conservation of the mass of each phase [25,27], as some leakages between the phases may occur.This drawback is rooted in the phase-field method and is here mitigated with the adoption of the profile-corrected formulation.In the present cases, mass leakages are limited to at most 8% of the dispersed phase mass and occur only in the initial transient phase; once the statistically-stationary condition is reached, the mass of each phase keeps constant.

E. Simulation set-up
We consider a turbulent channel flow at a shear Reynolds number Re τ = 300 for all the cases.The computational domain has dimensions L x × L y × L z = 4πh × 2πh × 2h, which corresponds to L + x × L + y × L + z = 3770 × 1885 × 600 wall units.The domain is discretized with N x × N y × N z = 512 × 256 × 513 grid points; the computational grid has uniform spacing in the homogenous directions, while Chebyshev-Gauss-Lobatto points are used in the wall-normal direction.The flow is driven by an imposed constant pressure gradient in the streamwise direction.We consider two surface tension values, which are set via the Weber number: W e = 1.50 (higher surface tension) and W e = 3.00 (lower surface tension).The selected values are characteristics of air/water mixtures [43].For each surface tension value (i.e. for each Weber number), we first keep a unitary density ratio and we analyze the effect of different viscosity ratios: from η r = 0.01 (less viscous bubbles) up to η r = 100 (more viscous bubbles).Then, we keep a unitary viscosity ratio and we consider different density ratios: from ρ r = 1 (matched density bubbles) down to ρ r = 0.001 (lighter bubbles).Finally, to evaluate the combined effect of density and viscosity differences, we consider a case in which both bubble density and viscosity are smaller than those of the carrier fluid: ρ r = 0.1 and η r = 0.1.In addition, we perform a single-phase flow simulation as a reference case and to provide initial velocity fields for the multiphase simulations.It is worthwhile noting that when different properties (i.e.density and viscosity) are considered, the local value of the Reynolds number changes as well as the range of spatiotemporal scales that needs to be resolved to fulfill the DNS requirements.These modifications can be appreciated from table II E in which we show an estimate of the turbulence length scale inside the dispersed phase (computed from the definition of the Kolmogorov length scale), η + k,d , the grid resolution, the final average bubble-size, ⟨d + eq ⟩, and its root mean square value, RMS(d + eq ), for all the different combination of density and viscosity ratios considered as well as for the reference single-phase case.The bubble size has been characterized using the equivalent diameter, d + eq , i.e. the diameter of an equivalent spherical bubble with the same volume as the bubble considered [17]: where V + is the volume of the bubble.All dimensions are reported in wall units (based on the carrier flow shear Reynolds number) and refer to the channel centre, where most bubbles are located.The Kolmogorov scale, which is used here to provide an estimate of the smallest length scale inside the bubbles, has been computed as follows: where ϵ + is the dissipation at the channel center evaluated in the region characterized by ϕ ≥ 0 (i.e.inside the bubbles), η r is the viscosity ratio and Re τ is the shear Reynolds number.We can observe that for almost all the cases presented here, the estimated Kolmogorov scale is of the order of the grid spacing thus ensuring a correct resolution of all the relevant flow scales.Only for the cases with η r ≤ 0.1 (most critical cases due to the largest local Reynolds number increase), the smallest flow scales (which are found inside the bubbles) cannot be completely resolved.From table II E, we can also observe that the average bubble size is always at least one order of magnitude larger than the grid spacing.
For the phase field, the Cahn number is set to Ch = 0.02.This value is selected based on the grid resolution: at least three grid points are required across the interface to accurately describe the steep gradients present [26].The phase field Péclet number has been set according to the scaling P e = 1/Ch = 50, to achieve convergence to the sharp interface limit [27,44].More refined grids allow to reduce the thickness of the interface and to adopt smaller Cahn numbers.However, the resulting computational cost is much larger: grid resolution needs to be refined along all three directions, as the orientation of the interfacial layer is arbitrary, and the time step has to be reduced as well to satisfy the Courant-Friedrichs-Lewy condition.Overall, the computational cost of a simulation with an halved Cahn number is roughly 16 times larger: grid refinement makes the simulation eight times more expensive and the time step limitation makes the simulation twice as expensive.
At the beginning of each simulation, a regular array of 256 spherical droplets with diameter d = 0.4h (corresponding to d + = 120 wall units) is initialized in a fully-developed single-phase turbulent channel flow.The total volume fraction of the dispersed phase is Φ = V d /(V c + V d ) ≃ 5.4%, being V d and V c the volume of the dispersed and carrier phase, respectively.As the array of spherical bubbles is suddenly released in a single-phase turbulent flow, turbulent velocity fluctuations strongly perturb the interfacial profile; during this initial coupling phase, mass leakages among the phases may occur [25,27] After this initial transient, the mass of each phase keeps constant over time.While the initial condition chosen for the dispersed phase may seem unphysical, after a short transient, memory of the initial condition is completely lost and the results are not affected by the initial condition selected [17] ), for all the different simulations performed.All dimensions are reported in wall units; Kolmogorov scale is measured at the channel centre.Single-phase flow values at the channel centre have been also reported as a reference.
been tested (e.g., the injection of a thin liquid sheet at the channel center) and the same statistically statisticallystationary results were obtained.We selected the current initial configuration as it reduces the time required to reach statistically-stationary conditions.

III. RESULTS
We present here the results obtained from the analysis of the simulation database, starting from the effects of the density ratio, viscosity ratio, and Weber number on the topology of the dispersed phase (number of bubbles) and on its topological changes (coalescence and breakage rates).Then we evaluate the effects of these parameters on the shape and deformation of the bubbles studying the local curvature of the interface and the time evolution of the interfacial area.Finally, we investigate the flow modifications produced by the bubbles by analyzing the mean velocity profiles and the turbulent kinetic energy inside the bubbles.All the results will be presented according to the following color code: a red-colors scale is used to show the density ratio variations and a blue-colors scale to show the viscosity ratio variations.The case with both non-matched density and viscosity is represented in green, while the reference case (matched density and matched viscosity) is shown in black.
A. Bubbles: number and topological modifications

Number of bubbles
The topology of the dispersed phase is the direct consequence of the ultimate competition between breakage and coalescence events.To obtain a first qualitative insight of the effects of density ratio, viscosity ratio and Weber number on the statistically-stationary number of bubbles (i.e.once the effect of the initial condition is completely dissipated), we can consider figure 1. Panel (a) refers to W e = 1.5, while panel (b) to W e = 3.0.In each panel of figure 1, four snapshots of the multiphase system at statistically-stationary are arranged in a plot according to the values of density (horizontal axis) and viscosity (vertical axis) ratio of each case.The surface of the bubbles, identified as the

System
Reτ  50 and W e = 3.00.For each Weber number, we consider four density ratios: from ρr = 0.001 up to ρr = 1.000; five viscosity ratios: from ηr = 0.01 up to ηr = 100 and a combined case ρr = 0.1 and ηr = 0.1.In addition, a single-phase flow simulation has been also conducted.
iso-contour ϕ = 0, is reported at the time instant t + = 4000 (statistically-stationary conditions); in the background the contour map of the turbulent kinetic energy, TKE= (ρ/ρ c )(u ′2 + v ′2 + w ′2 )/2 (where ρ identifies the local density value, ρ d in the bubbles and ρ c in the carrier phase), on a x + − y + plane located at the channel centre is shown.Among all cases, we select those with the extreme values of the density (ρ r = 0.001 -η r = 1) and viscosity ratio (ρ r = 1 -η r = 100 and ρ r = 1 -η r = 0.01).As a reference, also the matched density and viscosity case (ρ r = 1η r = 1) is shown.We can observe that for W e = 1.5 (figure 1a), the number of bubbles remains almost unchanged when both density and viscosity contrasts are introduced in the system.For W e = 3.0 (figure 1b), the number of bubbles is higher in all the cases, compared to W e = 1.5.If we then look along the density axis (namely to the pictures in the central row) of figure 1b, we see that the number of bubbles is quite similar in the two cases, suggesting a negligible effect of density for the range of values considered here.By opposite, looking along the viscosity axis (thus to the pictures on the right column), we notice that viscosity does play an important role, as the number of bubbles significantly reduces from η r = 0.01 to η r = 100, with a more marked difference between η r = 1 and η r = 100, than between η r = 0.01 and η r = 1, thus hinting that the viscosity difference among the phases may actually be the relevant factor, rather than the viscosity ratio.
To evaluate these results more quantitatively, we compute at each time the number of bubbles, N (t + ), normalized by the initial bubbles number, N 0 .Figure 2 shows the results obtained for all the combination of density and viscosity ratios considered, and for the two Weber numbers as well.Left column refers to W e = 1.5 (figure 2a,c,e), while the right column to W e = 3.0 (figure 2b,d,f ).The top, middle and bottom rows show, in order, the effects of the density ratio, viscosity ratio and of their combination.
We start by analyzing the effect of Weber number solely and we consider the matched density and viscosity case (black lines in figure 2a-d ).For W e = 1.5, the number of bubbles decreases monotonically: coalescence events dominate the initial transient phase (up to t + = 2000).Then a balance between breakage and coalescence events is attained and the number of bubbles settles on a stationary value, N (t + )/N 0 ≃ 0.1.Likewise, for W e = 3.0, an initial transient mainly characterized by coalescence events can be also observed.However, this phase ends at an earlier time (about t + = 500) and is followed by a statistically-stationary condition where breakups and coalescences alternately prevail on each other.
Comparing simultaneously the plots at W e = 1.5 (figure 2a,c,e), we can observe that the effects of both density and viscosity ratios (and of their combination) are very small.This behavior can be traced back to the dominant role played by surface tension forces.The Weber number quantifies the relative importance of surface tension forces with respect to inertial forces: the lower is the Weber number, the stronger is the action of surface tension in controlling bubbles dynamics.Thus, for W e = 1.5, the surface tension forces are dominant and are those determining the topology of the dispersed phase (i.e.number of bubbles).For the higher Weber, surface tension forces are weaker in comparison, and density and viscosity ratios effects become more significant.In particular, for W e = 3.0 (figure 2b,d,f ), the statistically-stationary value obtained for the number of bubbles shows a marked dependence on the viscosity ratio.
As the dispersed phase dynamics for the cases at W e = 1.5 are dominated by surface tension forces, we focus on the cases at W e = 3.0 to investigate the effects of density and/or viscosity ratios.First, we consider the effects of the density ratio solely.Figure 2b shows the time evolution of the number of bubbles for different density ratios (from ρ r = 1.0 down to ρ r = 0.001) and a fixed unitary viscosity ratio.We notice that the influence of the density ratio on the number of bubbles is small: the red-colors lines do not depart in average from the black reference line, nor from each other.Hence, no significant modifications are introduced in the topology of the dispersed phase when density contrasts are present between the phases (with respect to a two-phase system with uniform density).This behavior suggests that, for the range of density ratios considered, the external inertial forcing is the main factor that determines the bubble size and thus the dispersed phase topology.In contrast, the density (and thus the inertia) of the bubble plays a negligible role in determining the dispersed phase topology.
On the other hand a marked effect of the viscosity ratio alone can be observed, figure 2d.We observe in this case a much clearer trend: after the initial transient the curves depart from each other and set on different equilibrium values once statistically-stationary conditions are reached.In particular, as the viscosity ratio is increased, the statisticallystationary number of bubbles is reduced.For high viscosity ratio (η r > 1) fragmentation is prevented, coalescence dominates and only a few bubbles are present in the channel.By opposite, for low viscosity ratio (η r < 1) breakups are favored, the average bubble size decreases, and the resulting number of bubbles is slightly larger when smaller viscosity ratios are considered.Hence, it is evident that viscosity acts as a stabilizing factor, in a similar way as surface tension does.Indeed, it is interesting to observe that the behavior of the number of bubbles for η r = 100 at W e = 3.0 (high viscosity) resembles those of the cases at W e = 1.5 (high surface tension, figure 2c).This suggests that a very high viscosity ratio can compensate a low surface tension and produce similar results in terms of topology.A physical argument that can explain the action of viscosity is related to the deformations that the external turbulent flow is able to induce on the bubble.When the internal viscosity is larger than the external one, the larger internal viscous dissipation damps all the turbulent fluctuations produced by the external flow.This hinders large deformations of the bubble surface and, as a consequence, it reduces the possibility of bubble breakage.
Finally, we analyze the combined effects of density and viscosity ratios.In figure 2f, we report the results obtained from the case ρ r = 0.1 and η r = 0.1 and from two cases with one matched property and one non-matched property, ρ r = 0.1 and η r = 1 (red line) and ρ r = 1 and η r = 0.1 (blue line).We can first note that these two latter cases, where only one property is non-matched, exhibit a very similar behavior for the entire duration of the simulation.This is consistent with our previous observation: the influence of the density ratio is almost negligible (figure 2b) and the effects of the viscosity ratio are relatively small for η r = 0.1 (figure 2d ).Then, we observe that the combined case (green line) does not deviate largely from the other two cases.This indicates that a simultaneous reduction of the density and viscosity ratios does not remarkably modify the general picture for the range of density and viscosity ratios here tested.Nevertheless, it is interesting to observe that the green line lies above the red and blue lines for a longer timespan, indicating that the number of bubbles for the combined case is slightly higher than in the other two cases.

Breakage and coalescence rates
The evolution of the number of bubbles provides useful insights on the time behavior of the dispersed phase topology, although it only shows the net outcome of the competition between breakage and coalescence events.To evaluate whether density and viscosity differences among the phases affect breakage and coalescence dynamics, we compute the instantaneous number of breakage and coalescence events.Evaluating these effects is not only crucial to better understand the involved physics, but is also extremely important for the development of accurate coalescence and breakage kernels [45].The time behavior of the breakage and coalescence is directly linked to the number of bubbles present in the channel, as hinted by the balance population equation [46]: where N (t + ) is the number of bubbles and Ṅb (t + ) and Ṅc (t + ) are respectively the breakage and coalescence rates.We compute the breakage and coalescence rates counting the number of breakage or coalescence events that occur within a set temporal window ∆t + (see Appendix A for details): where the temporal window has been chosen equal to ∆t + = 300.As the number of breakage and coalescence events that occur in a certain temporal window is also influenced by the number of bubbles present in the channel [17], we normalize the breakage and coalescence rates, Ṅb (t + ) and Ṅc (t + ), by the instantaneous number of bubbles N (t + ).Being the description of coalescence and breakage events in numerical simulations influenced by grid resolution [3,6,17,28], a convergence study has been also performed to ensure that the grid employed is sufficient to obtain convergent results, please refer to Appendix B for details.
Figure 3 shows the results obtained for all cases examined: breakage rate is plotted over time as a positive quantity, while coalescence rate as a negative quantity, being them related to an increase and decrease of the number of bubbles, respectively.We will first discuss the effect of the Weber number comparing the left column (figure 3a,c,e) with the right column (figure 3b,d,f ).For W e = 1.5 (left column), the breakage and coalescence rates behave nearly in the same way for all the combinations of density and viscosity ratios.After the initial transient where the behavior of the rates is influenced by the selected initial condition for the phase-field, an equilibrium is reached at about t + = 1000 where both rates set on a constant value.At this stage, bubbles keep on breaking and coalescing, but with the same rate, thus maintaining their number in statistical equilibrium.This value of the Weber number does not allow density and viscosity contrasts to substantially modify the evolution of bubbles topology, as a good correspondence among the curves can be noticed in all the plots.Indeed, when a low Weber number is considered the deformability, which is a crucial factor for coalescence and breakage events, is mainly determined by surface tension forces that dominate over density and viscosity contributions.For W e = 3.0 (right column), the results are qualitatively and quantitatively different: breakage and coalescence rates reach in general larger values, and some significant deviations among the curves are visible.This is a direct consequence of the larger Weber number: surface tension forces, which are smaller in magnitude, weakly counteract turbulent velocity gradients, that can more easily deform and break the bubbles.Thus, we observe a larger number of breakage and coalescence events due to the larger deformability of the bubbles, as can be appreciated from figure 2b,d,f.In addition, for this larger Weber number, we can clearly observe how the density and viscosity ratios play a much more important role in the dynamics of breakage and coalescence events (with respect to W e = 1.5).
For this reason, we move now to discuss the effect of non-matched density or viscosity on the cases at W e = 3.0 in more detail.Figure 3b shows the breakage and coalescence rate for different values of the density ratios.In the first transient phase, all cases manifest a very high frequency of both breakage and coalescence events, slightly larger for coalescences at the very beginning (coherently with the evolution of the number of bubbles shown in figure 2b).Later on, both rates stabilize and set on two equal (in magnitude) stationary values.Although a clear trend among the different density ratios cannot be observed, it is worth noticing that all the rates seem slightly larger when sub-unitary density ratios are considered (especially in the early stage of simulations).Overall, these observations suggest that density differences between the phases do not introduce remarkable changes in the dispersed phase topology and on its modifications: the number of bubbles and breakage and coalescence rates are weakly influenced by changing the density ratio.
Moving now to the effect of the viscosity ratio, figure 3d depicts the time evolution of the breakage and coalescence rates obtained for different viscosity ratios (and a fixed unitary density ratio).Again, once the initial transient is finished, a statistically-stationary phase can be distinguished for all cases.From a qualitative viewpoint, coalescence is predominant at the beginning of the transient (consistently with the behavior reported in figure 2d ); then relatively high values for both rates are maintained during the rest of the transient, until they stabilize on steady values.The cases, however, deeply differ from a quantitative point of view.We see in this case that the rates significantly change when the viscosity ratio is changed: both breakage and coalescence rates decrease in magnitude as the viscosity ratio is increased (i.e. when bubble viscosity is increased).This modification of the breakage and coalescence rates is clear when the case η r = 100 is considered: the statistically-steady value of both rates is smaller than the one attained by the other cases.A similar trend was experimentally measured by Eastwood et al. [47] for the breakup of immiscible fluid particles in a turbulent jet: it was observed that the breakage rate of the droplets scales inversely with the inner bubble capillary number (ratio between bubble viscous forces and surface tension forces).Present results seem to confirm this finding: bubble viscosity and the corresponding viscous forces, acting as a damper of external velocity fluctuations [16], make bubbles less deformable and the probability of breakage and coalescence decreases.
Finally, we discuss the combination of density and viscosity contrasts (figure 3f ).The three curves do not deviate considerably from each other and a clearcut trend cannot be appreciated.As the density effect is generally unimportant and the viscosity one shall be small for η r = 0.1, the case ρ r = 0.1 -η r = 0.1 does not give us clear information on how density and viscosity effects combine together.

Interfacial area
A bubble released in a turbulent flow is constantly subjected to deformations due to the action of turbulent fluctuations [48,49].Turbulence fluctuations deform and stretch the bubble and, if strong enough, can lead to breakage of the bubble.The result of turbulence actions in terms of deformation can be evaluated by computing the total interfacial area.This quantity gives a general indication of the average bubble deformation and also provides a quantification of the amount of energy stored at the interface [33,50,51].Indeed, in the hypothesis of constant surface tension (as in the present case), surface tension energy is proportional to the amount of interfacial area available [33,50,51].
With the aim of evaluating the effects of the simulations parameters (density ratio, viscosity ratio, and Weber number) on the interfacial energy, we compute the time behavior of total interfacial area, A(t + ), for all cases considered.The results are presented normalized by the initial value A 0 .In figure 4, the results are shown using the same arrangement adopted in the previous figures.To correctly interpret these results, it is necessary to make a preliminary remark.The area of the interface between the dispersed phase and the carrier fluid evolves in time depending on two factors: the evolution of the number of bubbles and the modifications of the shape of the bubbles.This concept can be explained by considering the following example: to have a minimal interface area, the dispersed phase should consist of a unique spherical bubble, since, for a given volume, the spherical shape is the one that minimizes the surface area.If we split this bubble into several smaller spherical bubbles the total interface area will increase, being the total volume constant.If these smaller spherical bubbles are then deformed and elongated the area will further increase, as for each bubble the same mass will be redistributed in a way that makes it more exposed to the external flow.Thus, when we look at the evolution of the total interface area we are simultaneously observing the effect of the number of bubbles and of their deformation.We start by analyzing the effects of the density ratio for the cases at W e = 1.5, figure 4a.We notice an initial transient that is characterized by a nearly monotonic decrease of A(t + )/A 0 , for all the considered cases.In particular, during this transient, the curves corresponding to sub-unitary density ratios are superposed, while a remarkable discrepancy is visible between them and ρ r = 1.As soon as the flow reaches a steady behavior, all the curves differentiate and a trend becomes visible, where the higher is the density ratio the larger is the total interface area.Considering that for W e = 1.5 the number of bubbles is almost unaffected by the density ratio (figure 2a), this indicates that the trends observed in figure 4a are mainly caused by the bubble deformation: when smaller density ratios are considered, bubbles tend to be less deformed with respect to the case ρ r = 1.The origin of this behavior can be traced back to the local Reynolds number (i.e.evaluated using the bubble proprieties): as the density ratio is decreased, the inertial forces become smaller, the local Reynolds number decreases and less deformed bubbles are obtained.
For W e = 3.0 (figure 4b), we notice a similar but more irregular behavior.For all density ratios, the normalized interfacial area decreases and sets on stationary values that are higher than the final stationary values obtained for W e = 1.5 (figure 4a).This is coherent with the fact that increasing the Weber number, the number of bubbles increases, and so does the interfacial area.For this larger Weber number, the trend among the different density ratios is now less clear and the differences between the curves are slightly smaller.Nevertheless, consistently with the results obtained for W e = 1.5 (figure 4a), the matched density case (ρ r = 1) is clearly above all the other curves (ρ r < 1) for almost the entire time range of the simulations.Being the number of bubbles similar for all the cases shown in figure 4b, this seems to confirm that for smaller density ratios the overall interfacial area is reduced.
The viscosity effect can be appreciated in figure 4c,d.For W e = 1.5 (panel c), the total interface area is practically independent on the viscosity ratio and no significant changes can be observed.As the number of bubbles is similar for all cases (figure 2c), this indicates that no significant effects on the average bubble deformation are observed.Even though bubble viscosity does not play an important role in the average bubble deformation, we can anticipate that it still plays a role when more local quantities are analyzed (e.g.local curvature), see section III B 2. For W e = 3.0, a remarkable difference is present between η r = 100 (larger bubble viscosity) and all the other cases.This is consistent with the time evolution of the number of bubbles (figure 2d ).Indeed, when the statistically-stationary configuration is reached, the number of bubbles for η r = 100 is much lower than that obtained for the other ratios.As a result, the interfacial area is much lower than the other cases.For the other cases (from η r = 10 down to η r = 0.01), a clear trend cannot be observed thus suggesting that no large modifications of the average bubble deformation are obtained for η r < 10.However, as already anticipated for W e = 1.5, larger modifications are observed when local quantities are analyzed, see next section for details.
Finally, we discuss the combined effect of density and viscosity ratios (figure 4e,f ).For W e = 1.5, the case with both non-matched density and viscosity (green line) overlaps the case with non-matched density (red line) during the transient and in the final steady configuration, while in the first steady part it is intermediate between the two other cases, ρ r = 0.1 -η r = 1 and ρ r = 1 -η r = 0.1.On average the combined case is therefore closer to the non-matched density case, suggesting that the density ratio has a larger influence on the total interfacial area (and thus on the stretching of the bubbles) with respect to the viscosity ratio.This is confirmed by the plot for W e = 3.0, where the green line shows again values that on average are much closer to the non-matched density case (i.e.ρ r = 0.1).

Probability density function of mean curvature
The evolution of the total interface area gives us an idea of the overall behavior of the average deformation of the bubbles in presence of density and viscosity contrasts.However, being an average indication, it does not provide a clear indication of the local deformations of the bubbles surface.To obtain a deeper understanding of the deformation, we examine the probability density function (PDF) of the local interface mean curvature in the final statisticallystationary configuration.The mean curvature, K + , can be computed as the divergence of the local normal vector n, which in turn can be defined from the phase variable ϕ [52,53]: We compute the mean curvature, K + , for each point on the surface of the bubbles, corresponding to the points of the iso-level ϕ = 0.The resulting curvature values tell us how much the bubbles deviate from their spherical equilibrium shape, giving rise to small bumps and ripples in the surface when K + is highly positive, or small dimples when K + is highly negative.
From figure 5, we can appreciate the effect of density and viscosity on the mean curvature from a qualitative point of view.The shows for W e = 1.5 (figure 5a) and W e = 3.0 (figure 5b) four top views of the statistically-stationary configurations of the system.Bubbles are colored according to the local value of the mean curvature (blue-low; redhigh).Red areas correspond to bumps and ripples of the interface (positive curvatures), while blue areas to dimples (negative curvatures).
For W e = 1.5 (figure 5a), the effect of the density ratio can be observed by looking at the horizontal sequence of pictures (central row): we notice that moving from ρ r = 1 down to ρ r = 0.001 there is a slight decrease in the extension of both red and blue saturated regions, which correspond to very high and very low curvatures respectively.Therefore a reduction of the density ratio (i.e. a decrease of bubble density), leads to a smoother bubble surface, characterized by fewer ripples and dimples.In the vertical sequence of pictures on the right column, we can appreciate the effect of viscosity.We notice that the shape of the bubbles is qualitatively unchanged increasing the viscosity from η r = 0.01 to η r = 1.However, from η r = 1 to η r = 100 the shape changes remarkably: the irregularities that characterize the bubbles surface at η r = 1 disappear completely at η r = 100, where the surface becomes very smooth and the bubbles shape very closely resembles the spherical shape.Thus, the action of viscosity seems opposite to the one of density in terms of local deformation of the bubble surface: an increase of viscosity prevents the formation of high curvatures values (in magnitude), while an increase of density promotes the formation of large interface deformations.The two opposite trends obtained increasing the density or viscosity ratios can be interpreted in terms of local Reynolds or capillary numbers (i.e.evaluated using the bubble proprieties).An increase of the density ratio leads to an increase of the local Reynolds number and as a consequence, a more irregular surface of the bubbles is obtained.In contrast, an increase of the viscosity ratio, produces a decrease of the local Reynolds number (which also corresponds to an increase of the capillary number) and a smoother surface of the bubbles is attained.Interestingly, the entity of these effects depends on the value of the ratio considered: a slight effect of the density ratio can be observed when it is decreased of three orders of magnitude (from ρ r = 1 down to ρ r = 0.001), as well as for the viscosity ratio when reduced by two orders of magnitude (from η r = 1 down to η r = 0.01), while a more noticeable difference is visible when it is increased of two orders of magnitude (from η r = 1 up to η r = 100).Similar considerations can be obtained from the qualitative results obtained at W e = 3.0 (figure 5b).In this case, we can qualitatively appreciate similar effects for the density and viscosity ratios.These modifications, however, are now reflected on a much larger number of bubbles (larger Weber number).
To confirm these first qualitative observations, we compute the probability density function (PDF) of the mean curvature.Results are reported in figure 6 for different combinations of the density ratio, viscosity ratio, and Weber number.Left column (figure 6a,c,e) refers to W e = 1.5, while right column (figure 6b,d,f ) to W e = 3.0.Before analyzing each curve in detail, we can do some general observations.All curves are centered on a positive value of curvature and present an asymmetry with respect to the null value.Since positive curvatures correspond to convex surfaces and the null curvature corresponds to a flat surface, this is consistent with the fact that bubbles are in average convex, considering an outwards normal vector.Then, comparing the results shown in the left column (cases at W e = 1.5) against those reported in the right column (cases at W e = 3.0), we can appreciate the effect of the Weber number: for W e = 3.0 the curves are extended on a wider range of curvature values with respect to W e = 1.5.In particular, the curves are extended slightly towards negative values and considerably towards positive values, meaning that a higher Weber leads to a higher probability of having irregularities in the surface of the bubbles, especially bump or ripples-like irregularities.The higher probability of having large curvature values is also due to the presence of many smaller bubbles, which are intrinsically more convex (smaller radius) and closer to a spherical shape.
We study now the effects of the density ratio (figure 6a,b).We notice a trend for W e = 1.5 that becomes clearer for W e = 3.0: the cases with ρ r = 0.1, 0.01, 0.001 present a lower probability of having large curvatures (in magnitude) with respect to ρ r = 1.This effect is small for positive curvatures and more pronounced for negative curvatures.We can also observe that while the discrepancy between the reference case (ρ r = 1) and all other cases is clear, there is almost no difference among the cases ρ r = 0.1, 0.01, and 0.001.Interestingly, a similar trend was also reported in a previous work [54] that investigated the rise of bubbles in quiescent liquid.In particular, Cano-Lozano et al. [54] reported that for density ratios smaller than 0.128, a further decrease of the density ratio does not produce significant changes in the shape of the bubbles.This seems to suggest that the modifications produced by the density with respect to the case with ρ r = 1 (matched density case), are likely to be proportional to the density difference between the two phases (i.e.ρ c − ρ d ) rather than their ratio (i.e.ρ d /ρ c ).Further simulations, which consider superunitary density ratios, are however required to confirm this indication.Overall, present results (figure 4) indicate that when sub-unitary density ratios are considered, the probability of having large curvatures values, especially negative, and very stretched bubbles decreases.In other words, when the density of the bubbles is decreased with respect to the carrier density, it becomes more difficult for turbulence fluctuations to locally deform and stretch the bubbles, and in particular, it is difficult to create dimples and concave areas.A possible physical mechanism that supports present is the following: when an external perturbation reaches the deformable interface of a bubble, the bubble surface is modified and the perturbation then propagates to the internal bubble fluid.As bubble density is reduced, however, the propagation of this perturbation to the bubble fluid and thus to the rest of the bubble interface becomes less effective.Indeed, the inertia of the perturbation is modulated by the smaller bubble density and thus the magnitude of the inertial forces is reduced.As a result, viscous and surface tension forces increase their relative importance with respect to inertial forces, and the resulting bubble deformation is reduced.This behavior can be also justified considering the dispersed phase Reynolds number, i.e. the Reynolds number evaluated considering the dispersed phase density.As bubble density is reduced, so does the dispersed phase Reynolds number and the bubbles become less deformable and distorted, as can be also graphically appreciated from figure 1 comparing the case ρ r = 0.001 (orange bubbles) against the case ρ r = 1.000 (white bubbles).
To evaluate the influence of the viscosity, we consider figure 6c,d.A trend can be distinguished for both the Weber numbers: the PDFs become narrower as the viscosity increases.More specifically, the largest effect can be seen for η r = 100, where the range of possible curvatures is significantly reduced.The shrinkage of the pdf is less but still evident for η r = 10, and it becomes almost negligible for η r = 0.1 and η r = 0.01.Unlike density, the impact of viscosity is important for η r = 100 and η r = 10, while it becomes less important for η r = 0.1 and η r = 0.01.Indeed, for these two latter cases, no significant modifications can be appreciated from both Weber numbers.
Finally, the combined effects of the density and viscosity ratio can be evaluated from figure 6e,f.Interestingly, we observe that when both density ratio and viscosity ratios are decreased, the resulting PDF of the mean curvature lies in between the case ρ r = 0.1 (and matched viscosity) and η r = 0.1 (and matched density).This intermediate behavior can be traced back to the two opposite actions of density and viscosity on the mean curvature of the surface of the bubbles: while a decrease of the bubble density (i.e. of the density ratio) makes the bubbles surface more rigid and thus smoother, when bubble viscosity is decreased the bubbles become more deformable and ripples or dimples can be more easily formed on the interface.Thus, when we combine these two effects, these actions balance out and we obtain an intermediate trend.This result is already visible for W e = 1.5 and becomes clearer for W e = 3.0 where, thanks to the higher number of bubbles, a smoother statistic is obtained.

Mean velocity profiles
Once detailed the evolution of the dispersed phase topology, its modifications and the deformation and curvature of the bubbles, we move to analyze the flow modifications produced by the bubbles.We start by analyzing the macroscopic behavior of the multiphase mixture, in terms of flow-rate and mean flow statistics.In particular, we investigate the wall-normal behavior of the mean velocity profiles of the multiphase flow, and we compare them with the single-phase flow statistics at the same Re τ = 300.Indeed, we aim at understanding whether the injection of bubbles in a turbulent flow induces modifications to the mean velocity profile, especially when density or viscosity contrasts are present between the two phases.This aspect is widely studied and a common question that persists in the field concerns the capability of bubbles in generating drag reduction [15,20,[55][56][57][58].
Single-phase 0.1 1 10 100 W e = 3.0 W e = 1.5 FIG. 7. Wall-normal behavior of the streamwise mean velocity profiles.Left column refers to W e = 1.5, while the right column to W e = 3.0.Density ratios effects are shown on the top row for ρr = 0.001, 0.01, 0.1, 1. Viscosity ratio effects are shown on the middle row for ηr = 0.01, 0.1, 1, 10, 100.Finally, the combined effect of the density and viscosity ratios is shown on the bottom row for the case ρr = 0.1 and ηr = 0.1, with respect to the cases where only one effect is considered.As a reference, 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) is also reported with a dashed line.For all cases, with respect to single-phase, we observe a minor increase of the mean velocity.
Figure 7 shows the wall-normal behavior of the mean velocity profiles, computed by averaging the streamwise velocity along the streamwise and spanwise directions in the entire domain (both dispersed and carrier phase).The results are illustrated for all combinations of density and viscosity ratios considered, following the same arrangement of the previously presented statistics.In addition, the velocity profile relative to the single-phase case is shown with a black dashed line, and the classical law of the wall, u + = z + and u + = (1/κ) ln z + + 5.2 [59], is reported as a reference (with κ = 0.41 the Von Kármán constant [60]).We observe that in all the plots the velocity profiles perfectly collapse on each other in the vicinity of the wall, while tiny deviations can be observed in the central part of the channel, where most bubbles are located.In particular, in the core region of the channel, no differences can be appreciated varying the density and viscosity ratios.However, all multiphase cases are characterized by a slightly greater velocity with respect to the single-phase case.As in our simulations a constant mean pressure gradient is used to drive the flow, the observed flow-rate enhancement corresponds to a slight drag reduction.The drag reduction we observe is rather low in all the simulated cases (roughly 1 to 2%), and current results suggest that the presence of density and viscosity contrasts among the phases does not visibly impact it.These results are in agreement with previous works [20,61], which found that drag significantly depends on the bubble size.Specifically, they observe that large and deformable bubbles (obtained allowing bubbles to coalesce) migrate towards the central part of the channel and do not influence the drag significantly [12,13,55,62].By opposite, smaller bubbles (obtained not allowing bubbles to coalesce) move towards the near-wall region and lead to an increase of the drag [12,13,55,62].To support this argument, we can consider figure 8, which shows the scatter plot of the wall-normal location of each bubble over its equivalent diameter.Panel a refers to W e = 1.5 while panel b to W e = 3.0.The bottom and top walls are located at z + = 0 w.u. and z + = 600 w.u.. Two black dashed lines identify the critical condition for which the upper (or lower) part of the bubble interface intercepts the top (or bottom) wall.From a mathematical point of view, this condition can be identified imposing: where z + b is the distance of the center of the mass of the bubble from the closer wall, which can be computed as follows: where z + i is the wall-normal location of the i-th bubble and h + = 300 w.u. is the channel half-height in wall units.Hence, the equations that identify these conditions are: Analyzing the dispersion of the bubbles along the wall-normal direction, we can confirm previous intuitions: smaller bubbles tend to disperse along the entire height of the channel and can get rather close to the two walls while, by opposite, larger bubbles tend to accumulate at the center of the channel and stay farther away from the two walls.
It is worth pointing that despite a few points are located above (or below) the two black dashed lines (i.e. in the gray region), no collisions with the walls are detected.Instead, these points represent bubbles elongated along the streamwise or spanwise directions and thus with a larger d + eq with respect to the actual wall-normal size.Overall, the results presented in figure 7 corroborated by those reported in figure 8 seem to confirm the idea that bubble deformability is a crucial parameter for obtaining drag reduction [12,20,56,58,63].Indeed, bubble deformability plays a central role in determining the preferential distribution of the bubbles [12,32], which is directly linked to drag reduction [15,58].

Turbulent Kinetic Energy (TKE) of bubbles
After having analyzed the flow field in terms of mean velocity, we focus on the turbulence behavior inside the bubbles.The characterization of the flow inside the bubbles is of paramount importance in many applications.Indeed, internal circulation controls the transport of heat, mass, momentum and chemical species through the interface [64,65], the motion and deformation of the bubbles [66,67] and particle removal efficiency in scrubbing process [68,69].To characterize the mixing and flow behavior in the dispersed phase, we consider the turbulent kinetic energy (TKE) inside the bubbles.As in the carrier phase no significant modifications of the mean velocity profile (figure 7) and of turbulence statistics are observed, larger modifications are expected in the dispersed phase: the flow inside the bubbles is confined by a deformable interface and continuously forced by the external carrier flow.In addition, fluid properties (density and viscosity) are different.As a results, the magnitude of inertial and viscous forces is changed, as well as the local Reynolds and Weber numbers (i.e.evaluated using the dispersed phase properties).To give a first FIG. 8. Scatter plot of the wall-normal location of each bubble over its size for the different cases considered.The two black dashed lines identify the condition for which the interface of the bubble intercepts the closer wall in the hypothesis of a perfectly spherical bubble.Smaller bubbles tend to disperse along the entire channel height can get rather close to one of the two walls while larger bubbles tend to accumulate at the center of the channel.
qualitative idea of these modifications, we can consider the specific turbulent kinetic energy, TKE, whose definition is here recalled: where ρ is the local density (ρ d in the bubbles and ρ c in the carrier phase).Figure 9 shows the turbulent kinetic energy for two different simulations: panel (a) refers to the case with ρ r = 0.01 and matched viscosity and panel (b) to the case with η r = 0.01 and matched density.Both panels refer to the higher Weber number analyzed (W e = 3.0) and to the time instant t + = 4000, when for both cases a statistically-stationary configuration is attained.The two snapshots illustrate with a white-black scale the contour map of TKE on an x + − y + plane located at the channel center (z + = 0).The interface of the bubbles is marked with a white thin line.We notice that the flow structures in the carrier phase are qualitatively similar in the two pictures, while inside the bubbles the contour maps of TKE look very different and for ρ r = 0.01 and η r = 1 (panel a), low values of TKE inside the bubbles.In evaluating the results presented in panel a, however, it is important to make an important observation: although the energy content of the bubbles is rather low, velocity fluctuations are still present inside the bubbles.Indeed, the low values of TKE in the bubbles obtained for the case ρ r = 0.01 and η r = 1 are due to the low density that characterizes the bubbles: the prefactor ρ/ρ c present in the definition of TKE reduces the values obtained inside the bubbles.Shifting our focus to the case η r = 0.01 and ρ r = 1 (panel b), we can appreciate here the presence of many vortical structures characterized by an energy content similar to that of the carrier phase.Interestingly, the characteristic length scale of these turbulence structures is much smaller than that of the carrier phase.This observation can be traced back to the smaller viscosity of the dispersed phase that results in a larger local Reynolds number, as also observed in other multiphase flow instances [51,70].Turbulence inside the bubbles is the mechanism that can increase or decrease transfer rates across the interface [54,71].To quantify more closely this aspect, we compute the mean value of the specific turbulent kinetic energy inside the bubbles for all simulated cases, except for the combined case, and we collect the results in figure 10.To better evaluate the contribution of density and velocity fluctuations in the resulting TKE values, turbulent kinetic energy is evaluated using the complete definition (equation 25) in panel a while TKE is evaluated considering only the velocity fluctuations contribution in panel b, (i.e.TKE is reported normalized by the local density contribution ρ/ρ c ).The mean values of TKE are reported as a function of the density ratio (scale on the bottom part of the plot), viscosity ratio (scale on the top part of the plot), and Weber number (full circles for W e = 1.5 and empty circles for W e = 3.0).We start by analyzing the effects of the density and viscosity ratios shown in panel a, we can observe two opposite trends: as the viscosity ratio increases, the mean value of TKE inside the bubbles decreases of about one order of magnitude while, by opposite, increasing the density ratio, the mean value of TKE inside the bubbles rapidly increases (of about four orders of magnitude).This behavior reflects the modifications of the inertial and viscous forces inside the bubbles produced by the different dispersed phase density and viscosity.As the viscosity ratio is increased from η r = 0.01 up to η r = 100 (from left to right), viscous forces become dominant over inertial forces and thus local Reynolds number decreases.As a result, for low viscosity ratios, we observe small turbulent structures inside the bubbles characterized by significative TKE levels, while, for viscosity ratios larger than unity, turbulence structures cannot be sustained inside the bubbles (larger viscous dissipation) and bubbles are characterized by a low level of TKE.A similar trend, albeit in a slightly different simulation setup, was reported by Cano-Lozano et al. [54] that investigated the rise of bubbles in still liquid and observed a reduction of the velocity gradients for increasing values of the viscosity ratio.On the other hand, increasing the density ratio from ρ r = 0.001 up to ρ r = 1, inertial forces become dominant over viscous forces, the local Reynolds number increases and the bubbles are characterized by larger TKE values.Interestingly, we observe much stronger action of the density ratio on the mean value of the bubbles TKE.Indeed, if we compute the specific turbulent kinetic energy using equation 25, the resulting TKE values directly depend on the bubble density and, as we can see from panel a, present results roughly follow the ρ r scaling law reported with a dotted line.However, it is worthwhile observing that when the smallest density ratio is considered (ρ r = 0.001), results start to deviate from the ρ r scaling law: as the density ratio is reduced, we observe a reduction int the magnitude of the velocity fluctuations of about one order of magnitude.This deviation can be better appreciated in panel b, where TKE values are reported normalized by the prefactor ρ/ρ c , so that the contribution from velocity fluctuations alone can be better appreciated.The magnitude of velocity fluctuations is roughly constant when considering different density ratios, exception made for the lowest density ration, ρ r = 0.001, thus indicating that the specific TKE scales with the density ratio.Finally, we can consider the effect of the Weber number: increasing the Weber number, thus decreasing the surface tension, the TKE is slightly increased for all the cases.This trend can be attributed to the larger transfer of momentum that occurs when surface tension forces are weaker: as the interface becomes more deformable, the modulation effect of the interface becomes weaker and energy and momentum can be more easily exchanged between the phases.When the surface tension is reduced, in fact, the bubbles become more deformable and reasonably they are more likely to contain a greater amount of TKE.

IV. CONCLUSIONS
In this work, we studied the behavior of bubbles in a turbulent channel flow for different values of the density ratio, viscosity ratio, and Weber number.The investigation is based on direct numerical simulation of turbulence coupled with a phase-field method.First, we investigated the topology of the dispersed phase and its modifications.We found that the number of bubbles present in the channel is strongly influenced by the surface tension value (i.e. by the Weber number), in accordance with the results of previous studies [14,16,17].Besides, we observe that an increase of bubble viscosity with respect to the carrier (i.e. an increase of the viscosity ratio) has an important stabilizing role and leads to a remarkable increase of the maximum bubble stable diameter and thus to a decrease of the number of bubbles.By opposite, a reduction of the bubble density (i.e. a reduction of the density ratio), does not remarkably affect the dispersed phase topology.Similar findings are obtained from the analysis of the coalescence and breakage rates: an increase of bubble viscosity or surface tension (i.e. a decrease of the Weber number) leads to a reduction of the breakage and coalescence rates.In contrast, a modification of the density ratio has a marginal effect on the behavior of the breakage and coalescence rates.Secondly, we studied the surface stretching and curvature of the bubbles.We observed that these indicators are influenced by all three parameters investigated.In particular, larger viscosity ratios or lower density ratios or Weber numbers hinder the stretching of the bubbles and as a result the overall amount of interfacial area obtained is lower (with respect to the reference matched density and viscosity cases).
These observations are also reflected in the probability density function of the mean curvature: an increase of bubble viscosity, a decrease of bubble density or a decrease of the Weber number hinder the formation of ripples and dimples on the surface of the bubbles and thus high curvature values are less likely to be found.Thirdly, we evaluated the flow modifications produced by the swarm of bubbles in the background turbulent flow and in the dispersed phase.From a macroscopic point of view, no significant modifications are observed in the wall-normal behavior of the mean velocity profiles and only a minor increase of the flow-rate is observed for all bubbles-laden cases with respect to a single-phase flow, in accordance with previous results [14,16,17].Finally, as bubbles internal circulation play a key role in controlling the transport of heat, mass, momentum through the interface [64,65], we characterized the mixing in the bubbles by studying the turbulent kinetic energy of the bubbles.We observe a clear action of density and viscosity in modulating the turbulent kinetic energy of the bubbles.In particular, a decrease of the bubble density or an increase of the bubble viscosity lead to a remarkable decrease of the turbulent kinetic energy levels in the bubbles.B n+1,1 and B n+1,2 .Once all breakages have been detected, the algorithm looks for coalescence events.A coalescence event is detected whenever two separate droplets at time step n are assigned to the same droplet at time step n + 1.
In particular, referring to figure 12(c) bubbles C n i and C n j are both assigned to bubble C n+1 i , as it is the closest one to their estimated position.So far, only kinematic criteria have been used to determine the trajectory and eventual interactions (coalescences and breakages) of each bubble.Once all the trajectories at the present time step have been determined, the quality index and the balance are computed.In particular, the quality index, Q, is initialized at the beginning of the time step to the number of droplets at the current time step, N n ; every time volume is not conserved (within a certain small threshold) in all the translation, breakages and coalescences, the quality index is reduced by one.At the end of the time step, it is normalized by N n .Recalling the examples of figure 12, three checks on the volume conservation are performed depending on the type of event: To account for numerical errors that could occur in the calculation of the volume of each bubble (that would strongly reduce the quality index of the matching), a small tolerance ε (of the order of few percents of the volume of parent droplet) is used when checking for volume conservation.
The second parameter controlling the quality of the calculated trajectories is the balance, B. The total number of bubbles at each time step is known: N n at the current time step and N n+1 at the following one available.Once the number of breakage and coalescence events is known the balance can be calculated as: where B and N c are respectively the number of breakage and coalescence events that occur between time steps n and FIG.
12. Flow chart of the algorithm used to detect breakage and coalescence events in the post-processing of the simulations.
n + 1.The number of droplets at the current time step, N n , is increased whenever a droplet undergoes breakage into two bubbles and is decreased whenever two bubbles coalesce into one bubble.Here we make the assumption that all breakages are binary breakages and all coalescences involve only two parent droplets at a time.Thus, considering these two parameters, a fair matching of the trajectories is obtained with a quality index Q → 1 and a balance B = 0.This means that the volume is always matched (quality index never or rarely reduced) and no bubble is left out (balance equal to zero).Finally, once known the number of coalescence and breakage rates that occur between each time step n and n + 1, the coalescence and breakage rates, Ṅc and Ṅb , can be computed by counting the overall number of coalescence or breakage occurring in the temporal window ∆t + .Note that the temporal window used to track the trajectories of the bubbles is smaller than the temporal window used to compute the rates.The present algorithm considers only binary breakages and coalescences events.This assumption is not particularly limiting, as binary breakages/coalescences have the highest probability of occurrence [72][73][74].This assumption is also confirmed by the simulations performed: the quality index never reduces below 0.85 (so the volume is matched for at least 85% of all the translation, breakage and coalescence events) and at most few droplets are left unmatched (the balance is almost zero).the trend reported is similar and for all simulations, after an initial transient, both rates set to an equal value (in magnitude).Analyzing the value of the rates obtained, we can notice that some differences are present between the coarser grid resolution (triangles) and the intermediate grid resolution (circles).However, these differences become marginal when the intermediate grid resolution (circles) and the refined grid (squares) results are compared.Overall, present results suggest that the mesh employed is sufficient to investigate breakage/coalescence dynamics.

FIG. 1 .
FIG.1.Top view of four statistically-stationary configurations (t + = 4000) for different combinations of density ratios (ρr = 0.001 and 1) and viscosity ratios (ηr = 0.01, 1 and 100).Panel (a) refers to W e = 1.5, while panel (b) to W e = 3.0.The sub-panels are arranged in a plot using ρr as x-coordinate and ηr as y-coordinate.The effect of density can be appreciated in the sequence of panels on the middle row, while that of viscosity in the right column.The background of the plot shows the turbulent kinetic energy, TKE= (ρ/ρc)(u ′2 + v ′2 + w ′2 )/2 (white-low; black-high), computed on the central x + − y + plane of the channel.

5 NNFIG. 2 .
FIG.2.Time evolution of the number of bubbles, N (t + ), normalized by its initial value N0.Left column refers to W e = 1.5, while the right column to W e = 3.0.Top row: effect of density ratio, for ρr = 0.001, 0.01, 0.1 and 1 (with ηr = 1); Middle row: effect of viscosity ratio, for ηr = 0.01, 0.1, 1, 10 and 100 (with ρr = 1); Bottom row: combined effect of density and viscosity, for the case with ρr = 0.1, the cases ρr = 0.1, ηr = 1 and ρr = 1, ηr = 0.1 are reported for reference.On each line the left plot also includes the color code and a sketch with the definition of the property ratio considered (ρr, ηr or both ratios).

FIG. 5 . 1 −FIG. 6 .
FIG.5.Top view of the mean curvature of the bubble surface, K + , for four different combinations of density ratios (ρr = 0.001 and 1) and viscosity ratios (ηr = 0.01, 1 and 100) once a statistically-stationary configuration is reached (t + = 4000).Panel (a) refers to W e = 1.5 while panel (b) to W e = 3.0.The sub-panels are arranged in a plot using ρr as x-coordinate and ηr as y-coordinate.The effect of density can be appreciated in the sequence of panels on the middle row, while that of viscosity in the right column.Bubble surface (iso-level ϕ = 0) is colored according to the local value of the mean curvature (low-blue; high-red).

FIG. 9 .
FIG. 9. Contour map of the turbulent kinetic energy in a x + − y + plane located at the channel center (z + = 0).Panel (a) refers to the case ρr = 0.01 and ηr = 1 while panel (b) refers to the case ρr = 1 and ηr = 0.01.Both panels refer to the lower surface tension case (W e = 3.0) and to the time instant t + = 4000 (statistically-steady configuration).The interface of the bubbles is highlighted with a white line.For ρr = 0.01, bubbles and characterized by a low and uniform value of the TKE while, for ηr = 0.01, the TKE map is non-uniform and characterized by small scales fluctuations.

FIG. 10 .
FIG.10.Mean value of the turbulent kinetic energy (TKE) inside the bubbles.In panel a, TKE is evaluated using the complete definition of specific TKE (i.e.including the prefactor ρ/ρc) while, in panel b, TKE is evaluated considering only the velocity contribution (i.e.not considering the prefactor ρ/ρc).For both panels, a dashed line (W e = 1.5) and a continuous line (W e = 3.0) are used to show the behavior of TKE as the density or viscosity ratios are changed.Each value of TKE is marked with a circle (empty for W e = 1.5 and filled for W e = 3.0), with a red-color scale for the non-matched density cases and a blue-color scale for the non-matched viscosity cases, while the black color is used for the reference case.

FIG. 11 .
FIG. 11.Possible cases considered for the algorithm: panel (a) corresponds to a translation, panel (b) to a breakage and panel (c) to a coalescence.Red bubbles are at the current time step (n), while blue bubbles are at the next time step available (n + 1).Semi-transparent bubbles show the estimated position, x n+1 i,est .Arrows show the trajectory of the bubbles, ∆T u n i .
• Different initial conditions have TABLE I. Grid resolution, ∆x + , ∆y + and ∆z + c , Kolmogorov scale at the channel centre in the dispersed phase, η + k,d , average equivalent diameter of the bubbles, ⟨d + eq ⟩, and root mean square of the bubble equivalent diameter, RMS(d + eq

TABLE II
. Overview of simulations parameters.Wa analyze two Weber numbers: W e = 1.