The effect of droplet coalescence on drag in turbulent channel flows

We study the effect of droplet coalescence on turbulent wall-bounded flows, by means of direct numerical simulations. In particular, the volume-of-fluid and front-tracking methods are used to simulate turbulent channel flows containing coalescing and non-coalescing droplets, respectively. We find that coalescing droplets have a negligible effect on the drag, whereas the non-coalescing ones steadily increase drag as the volume fraction of the dispersed phase increases: indeed, at 10\% volume fraction, the non-coalescing droplets show a 30\% increase in drag, whereas the coalescing droplets show less than 4\% increase. We explain this by looking at the wall-normal location of droplets in the channel and show that non-coalescing droplets enter the viscous sublayer, generating an interfacial shear stress which reduces the budget for viscous stress in the channel. On the other hand, coalescing droplets migrate towards the bulk of the channel forming large aggregates, which hardly affect the viscous shear stress while damping the Reynolds shear stress. We prove this by relating the mean viscous shear stress integrated in the wall-normal direction to the centreline velocity.


I. INTRODUCTION
Two-fluid turbulent flows are found in many cases in industry and nature (Balachandar and Eaton, 2010), such as human arteries, industrial pipelines, and the injection of bubbles to enable drag reduction of ships (Ceccio, 2010). In all of these cases, surfactants are known to have dramatic effects on the flow, often by preventing coalescence (Takagi and Matsumoto, 2011). However, due to the multi-scale nature of the problems, the mechanisms by which coalescence affects drag are not fully known and understood yet. Thus, the objective of this work is to explain how coalescence affects drag in wall-bounded flows.
Many experimental studies of surfactants in multiphase flow have been made. Frumkin and Levich (1947) were the first to describe the mechanism by which the rising speed of bubbles in water is reduced by surfactants (see (Levich, 1962) for English version). Descamps et al. (2008) measured the wall shear stress in pipe flows of air bubbles in water, and found that larger bubbles produced less drag. Duineveld (1997) studied pairs of bubbles rising in a vertical channel; he showed that coalescence is prevented when the surfactant concentration is above a critical value. As well as preventing coalescence, surfactants produce other effects on bubbles, such as clustering (Takagi, Ogasawara, and Matsumoto, 2008), reduction of rising velocity (Frumkin and Levich, 1947;Levich, 1962), and reduction of shear-induced lift forces (Takagi and Matsumoto, 2011). Since all these effects can happen at the same time, the effect of different coalescence rate is difficult to highlight; on the other hand, simulations allow us to eliminate these effects, and focus solely on the impact of coalescence.
The majority of numerical multiphase flow studies have been made using interface-tracking methods, such as the front-tracking (FT) method (Unverdi and Tryggvason, 1992). Front-tracking simulations of homogeneous-isotropic flows a) Corresponding author: marco.rosti@oist.jp (Druzhinin and Elghobashi, 1998) are well suited for investigating the effect of droplet size on the turbulent length scales, such as bubble arrays Tryggvason, 1998, 1999) or channel flows (Lu, Biswas, and Tryggvason, 2006;Dabiri, Lu, and Tryggvason, 2013;Tryggvason and Lu, 2015;Tryggvason, Ma, and Lu, 2016;Lu, Muradoglu, and Tryggvason, 2017;Ahmed et al., 2020). An advantage of shear flow and channel-flow simulations is the ability to measure the effective viscosity and flow rate, which can then be compared with experiments. In the case of interface-tracking simulations of channel flows, Lu, Biswas, and Tryggvason (2006) simulated laminar bubbly upflows and downflows, Dabiri, Lu, and Tryggvason (2013) showed that more deformable bubbles produced lower drag, and Lu, Muradoglu, and Tryggvason (2017) modelled bubbles with insoluble surfactant, and Ahmed et al. (2020) with soluble surfactant, showing their main effects. However, none of the interface-tracking studies cited here includes a model for the breakup or coalescence of droplets, with only a few recent works tackling these phenomena Tryggvason, 2018, 2019).
Interface-capturing methods, such as the volume-of-fluid (VOF) method (Noh and Woodward, 1976), naturally allow coalescence and breakup of droplets (Elghobashi, 2019). Interface-capturing simulations of homogenous isotropic turbulence (Dodd and Ferrante, 2016;Perlekar et al., 2012;Komrakova, Eskin, and Derksen, 2015;Bolotnov, 2013) and shear flows  have shed some light on the effect of coalescence on turbulence. Notably, Dodd and Ferrante (2016) and Maxey (2017) showed that coalescence is a source of turbulent kinetic energy, while breakup is a sink. Scarbolo, Bianco, and Soldati (2015) investigated the effect of Weber number on breakup and coalescence, Soligo, Roccon, and Soldati (2019) modelled surfactant laden drops in turbulent channel flows, while Bolotnov et al. (2011) used the level-set method to simulate bubbly channel flows. Roccon et al. (2017) investigated the coalescence and breakup of large droplets in channel flow using the phase field method. Interface capturing methods are known to over-predict coalescence rates, because numerical coalescence occurs whenever the film thickness is less than the numerical grid spacing. In contrast, in real fluids film rupture occurs at molecular length-scales, which are in the tens of nanometres, orders of magnitude smaller than the Kolmogorov length Soligo, Roccon, and Soldati, 2019). A number of methods have been used to reduce the coalescence rate of interface capturing methods, such as adaptive grid refinement (Innocenti et al., 2021), film drainage models (Thomas, Esmaeeli, and Tryggvason, 2010), coupling to molecular dynamics simulations (Chen et al., 2004), and artificial forces .
In this paper, we use the front-tracking method to make simulations of droplets which cannot break up or coalesce, and we use the volume-of-fluid method to make simulations of droplets that easily break up and coalesce. As we are interested in the effects of coalescence, we do not use any methods to reduce the volume-of-fluid coalescence rate. The two methods give idealized models of a mixture saturated with surfactants (FT), and completely clean mixture (VOF). Aside from coalescence and breakup, the physical properties (surface tension, viscosity, density, etc.) of the fluids in the two methods are identical. To the authors' knowledge, this is the first direct comparison of coalescing and non-coalescing droplets in a turbulent channel flow.
The manuscript is organised as follows. First, in section II, we describe the mathematical model governing the problem at hand and the numerical techniques used to solve them numerically. In particular, we describe our chosen interface-tracking and interface-capturing methods in more detail. Section III reports the values of the parameters explored in our simulations. In section IV, we present statistics of the flow to elucidate how coalescence is affecting drag in the channel. Finally, section V gives conclusions and places them in the context of the current literature.

II. GOVERNING EQUATIONS AND FLOW GEOMETRY
We consider turbulent channel flows such as those shown in figure 1.
The numerical domain has size × × = 6 × 2 × 3 , where is the half-height of the channel. The flow is laden with an ensemble of droplets, initially spherical with radius = ∕8 and randomly arranged. We impose periodic boundary conditions in the streamwise ( ) and spanwise ( ) directions, while the nonslip and non-penetration boundary conditions are enforced at the two walls = 0 and = 2 . An imposed pressure gradient , uniform throughout the domain and constant in time, sustains the flow in the direction. Balancing the forces on the fluid in the direction, we obtain an expression for the shear stress at the wall, ≡ ⟨ | =0 ⟩ = , showing that remains constant in time. Note that, here and in the rest of the manuscript, we use angled brackets to represent an average over the subscripted directions.
The Cartesian components of the fluid velocity field ( 1 , 2 , 3 ) ≡ ( , , ) are found by solving the incompress-ible multiphase Navier-Stokes equations at each location , where , ∈ {1, 2, 3}. Throughout this article, we use Einstein notation (Einstein, 1916) where repeated indices are summed over, and the subscript comma denotes partial differentiation, i.e.,  , ≡  . The scalar is the pressure field used to enforce the incompressibility constraint stated in equation (2). The density and dynamic viscosity are the local weighted averages among the two phases, i.e., = + (1 − ) and = + (1 − ) , where subscripts and denote properties of the dispersed and continuum phases respectively. In the above, represents the volume fraction of the dispersed phase in each computational cell of the domain, with = 1 in the dispersed phase and = 0 in the continuum phase. The Kronecker delta is used to ensure that the pressure gradient is imposed in the direction. The last term on the right hand side of equation (1) is the volumetric formulation of the surface tension (Popinet, 2018); it is the product of the surface tension coefficient , the interface local curvature , and the unit normal to the interface . Note that we used ( ) in equation (1) to represent the surface delta function, which is zero everywhere except for the surface at the interface between the two phases. ( ) has dimensions of inverse length.

A. Discretisation of the Navier-Stokes equations
For simulations of coalescing and non-coalescing droplets, we use near-identical numerical methods to solve the momentum and continuity equations (eqs. 1 and 2). This ensures that any difference in our results is due to the droplets, not the integration scheme.
Equations 1 and 2 are solved numerically using a finite difference method on a fixed Eulerian grid with a staggered arrangement, i.e., fluid velocities are located on the cell faces and all other variables (pressure, density, viscosity, volume-offluid, etc.) are located at the cell centres. All the spatial derivatives appearing in the equations are discretised with secondorder central differences, except for the convective terms in the FT simulations where the QUICK scheme (Leonard, 1979) is used instead. In the single-phase (SP) and VOF simulations, time integration is performed with the Adams-Bashforth method. In the FT simulations, time integration is performed with a predictor-corrector method, in which the first-order solution (Euler method) serves as a predictor which is then corrected by the trapezoidal rule (Tryggvason et al., 2001;Farooqi et al., 2019). Both schemes are second-order in time. Finally, in regards to the pressure solver, the fractional step technique (Kim and Moin, 1985) presented by Dong and Shen (2012) and Dodd and Ferrante (2014) is adopted, allowing for direct solution of a constant-coefficient Poisson equation using an FFT-based solver, even in the presence of density differences among the carrier and dispersed phases.

B. Volume-of-fluid method
We use the volume-of-fluid (VOF) method to simulate droplets undergoing topological changes, i.e., coalescence and break-up. This is an Eulerian-Eulerian technique in which the fluid phases are tracked using the local volume fraction scalar field . Since Noh and Woodward (1976), a number of variants of the VOF method have been developed (Youngs, 1982(Youngs, , 1984Puckett et al., 1997;Rider and Kothe, 1998;Xiao, Honma, and Kono, 2005;Yokoi, 2007). Here we use the multi-dimensional tangent of hyperbola for interface capturing (MTHINC) method, developed by Ii et al. (2012). In this method, we use a smooth hyperbolic tangent function to approximate the interface, where is a parameter controlling the sharpness of the interface, and a normalisation parameter to enforce ∭ = in each cell. is a three-dimensional function in the cell, with the same normal and curvature as the interface. Normals are evaluated using the Youngs approach (Youngs, 1982), while the surface tension force appearing in the momentum equation (1) is computed using the continuum surface force (CSF) approach (Brackbill, Kothe, and Zemach, 1992). See Rosti, De Vita, and Brandt (2019) for a detailed description of the volume-of-fluid code employed in this work, and in several other works De Vita et al., 2019). See Ii et al. (2012) andDe Vita et al. (2020) for validations against numerical benchmarks and experiments.

C. Front-tracking method
We use the front-tracking (FT) method to simulate droplets that can deform, but cannot break up or coalesce. This is an Eulerian-Lagrangian scheme in which the interface between the phases is tracked by a "front", composed of triangular elements. The method was introduced by Unverdi and Tryggvason (1992), with many refinements over the past 30 years (Tryggvason et al., 2001;Tryggvason, Scardovelli, and Zaleski, 2011), including techniques to correct for errors in volume conservation of the phases (Takeuchi and Tryggvason, 2020). The surface tension force acting on the th element is a volume integral of the surface tension force from equation (1), where and are the area and perimeter of the th element and is the tangent to the perimeter. The force is then interpolated onto the Eulerian grid by means of a conservative weighting function and used to update the fluid velocity, which in turn is used to update the position of the interface. As the interface evolves, the unstructured grid can greatly deform, resulting in a non-uniform grid. Thus, periodical restructuring of the Lagrangian grid is performed to maintain a nearly uniform size, comparable to the Eulerian grid size. See Muradoglu and Tryggvason (2014) for a detailed description and validation of the front-tracking code employed in this work, and used in several other works (Izbassarov and Muradoglu, 2015;Lu, Muradoglu, and Tryggvason, 2017;Ahmed et al., 2020). Extensive tests of the front tracking method are shown in Tryggvason et al. (2001).

III. SETUP
Due to the different nature of the numerical schemes used to describe the presence of the interface, the numerical domain is discretised on two different sets of grids, both verified to  Table I: Details of each turbulent channel flow simulation performed in the present study. The first column gives a unique name to each run for ease of reference, and the second describes the colours and markers that are used in the following figures. Input variables are shown in the subsequent columns in the middle, and output statistics are shown in the three rightmost columns.
provide grid-independent results. The non-coalescing-droplet simulations use a uniform grid in the homogenous directions and a non-uniform grid in the wall-normal direction, with minimum spacing Δ = 3 × 10 −3 at the channel wall. The minimum spacing in wall units is Δ + ≡ Δ ∕ = 0.5, where and are defined later in this section. Overall, the grid size for the non-coalescing droplet simulations (FT) is × × = 576 × 240 × 288, which is comparable to that used in Dabiri and Tryggvason (2015), and gives around 24 Eulerian grid points per droplet diameter. Due to periodic restructuring, we also have around 24 Lagrangian grid points per droplet diameter. The single-phase and coalescing-droplet simulations (VOF) use a cubic uniform grid with spacing Δ + = 0.8, and total size × × = 1296×432×648. This grid has 108 points per initial droplet diameter. We use more grid points in the VOF simulations in order to accurately resolve breakup and coalescence events throughout the domain.
The values of the non-dimensional parameters used in the simulations are shown in table I. We consider a total volume fraction of the dispersed phase in the range 0% ≤ Φ ≤ 10%, with the continuum phase being denser and more viscous than the droplet phase, as the density ratio is fixed equal to ∕ = 50 and the dynamic viscosity to ∕ = 50 for all runs. Thus, the kinematic viscosity ≡ ∕ has ratio ∕ = 1 for all runs. The problem approaches the density and viscosity ratios of air in water ( ∕ ≈ 830, ∕ ≈ 55) while still being numerically tractable. The friction Reynolds number ≡ ∕ is set to 180 for all runs, where ≡ √ ∕ is the friction velocity. We define the capillary number as 0 ≡ 0 ∕ (where 0 is the bulk velocity of the single-phase channel flow) for which two values are considered, 0 = 0.05 and 0.10. Based on these capillary numbers, the friction Weber number ≡ 2 ∕ assumes values smaller or larger than unity. Finally, is the number of droplets at the start of the simulation, which are initially identical spheres in a random arrangement.
The three rightmost columns in table I report three output statistics: the bulk Reynolds number, ≡ ∕ , where ≡ ⟨ ⟩ is the bulk velocity; the bulk Weber number, ≡ 2 ∕ ; the centreline velocity in plus units + ≡ ⟨ | = ⟩ ∕ . In the next section, we present these and other statistics of the channel flows, and discuss their implications.

IV. RESULTS
We consider turbulent channel flows in which droplets can coalesce, and compare the results with a configuration where coalescence is not allowed. The flow is driven by a constant pressure drop, thus an increase in the flow rate or bulk velocity indicates drag reduction, while its reduction is evidence for drag increase. We start by considering the profile of the streamwise velocity + in the channel, reported in figure 2. The single-phase run SP0 shows the typical velocity profile of a turbulent channel flow, with regions of high shear at the walls and a flattened profile in the channel centre. The runs with coalescing droplets (VOF) mostly collapse onto the single-phase profile, showing only a slight reduction in + toward the centre. Whereas the runs with non-coalescing droplets (FT) show a significant reduction in + , which becomes more pronounced as Φ increases. Also, in the coalescing droplets runs, variation of the capillary number produces little change in + , while in the non-coalescing runs, the change in + with 0 is much more substantial. This is quantified in the inset of figure 2, which shows the bulk velocity in wall units + ≡ ⟨ ⟩ ∕ on the left axis, and the skin-friction coefficient ≡ 2 ∕ 2 on the right axis. We see that, relative to the single-phase run, the coalescing droplets produce a maximum increase of 4% in , whereas the non-coalescing droplets produce a maximum increase of 30%. In the case of non-coalescing droplets, the drag is highly dependent on 0 . The high 0 (i.e., more deformable droplets) runs show little change in whereas the low 0 (i.e., less deformable droplets) runs show a 30% increase in . Notably similar drag increases have been measured for rigid particles in channel flows by Picano, Breugem, and Brandt (2015) and Rosti and Brandt (2020  coalescence of droplets in the channel has a profound effect on the flow. Throughout this section, we present additional statistics of the flows in order to shed light on the mechanisms of this effect. Figure 3 shows the velocity profile again, this time on a semi-log scale in wall units + ≡ ∕ , and + ≡ ∕ , where ≡ ∕ is the viscous lengthscale (Pope, 2000). Away from the wall and the channel centre << << , i.e., the lengthscales affecting the flow are separated, and the single phase flow profile is approximately parallel to a line with constant slope (the dashed line). This is a manifestation of the loglaw for turbulent channel flows (von Kármán, 1930), which can be derived by assuming the quantity + + + has no dependence on + or ∕ (complete similarity). The flow profiles with coalescing droplets in figure 3 are in excellent agreement with the log-law, suggesting that coalescing droplets do not break the scale separation. However, the flow profiles with non-coalescing droplets are not in such good agreement, because these droplets have constant size , and ∼ , so scale separation is prevented, hence + + + shows a dependence on ∕ . To further quantify the effect of coalescence on the flow, we fit a log-law function to each flow profile in the region 30 < + < 100. Our log law function has the form: where 5.89 is the + intercept for run SP0, and Δ + is the shift relative to SP0. The inset of figure 3 shows how the vertical shift Δ + in the log-law region of the channel is affected by the volume fraction Φ and capillary number 0 for the different cases. Again, we see relatively small shifts for simulations with coalescing droplets, and large shifts for simulations with non-coalescing droplets. In particular, Δ + grows in magnitude with Φ, especially for the case with 0 = 0.05. This reinforces our observations of the bulk streamwise ve-locity shown in the inset of figure 2, that the less-deformable, non-coalescing droplets produce a significant drag increase.
To understand what generates the differences observed for configurations of coalescing and non-coalescing droplets, we focus our attention on the total surface area of the droplets. The total interface area is responsible for the overall surface tension stress, and impacts how droplets disperse across the channel. Figure 4 shows how the total interface area at steady state ⟨ ⟩ depends on the total volume fraction Φ of the dispersed phase. The figure shows that the non-coalescing droplets of the FT runs exhibit only 1% increase in area, due to deformation from their initial spherical shape. On the other hand, the coalescing droplets of the VOF runs show more than 80% reduction in interface area, as droplets coalesce and grow in size. In particular, when the volume fraction is large, droplets have a higher likelihood of colliding, and hence more coalescence, leading to a smaller value of ⟨ ⟩ ∕ 0 .
For the coalescing droplets, the interface area ⟨ ⟩ ∕ 0 shows no dependence on capillary number, differently from what was observed by Lu and Tryggvason (2018) and , who found that that as 0 decreases, surface tension increases, the droplets become more stable to perturbations, hence larger, thus leading to a smaller interface area ⟨ ⟩ ∕ 0 . However, in this case, 0 << 1, and the coalescing droplets are limited in size by the channel height, not surface tension. Figure 1b supports this hypothesis, as the coalescing droplets are comparable in size to the channel height.
The inset of figure 4 reports the time history of the interface area: the cases with non-coalescing droplets (FT) rapidly converge to a statistically steady-state, whereas for the coalescing droplets, convergence is reached long after, at about + ≈ 8000. Interestingly, we observe that the coalescing droplet runs with larger capillary number (VOFa) converge to steady-state more rapidly than the smaller capillary number runs (VOFb), i.e., the larger 0 runs show a higher rate of coalescence, although the steady-state areas are roughly the same. This is in contrast with simulations of droplet coales-  table I. For ease of comparison we have moved the Φ = 5% and Φ = 10% volume fraction profiles upwards by + = 5 and + = 10, respectively. In the region 30 < + < 100 shaded in grey, we fit a log-law equation + = ln + 0.41 + 5.89 + Δ + (grey dashed line). Inset: The vertical shift Δ + for each run. Runs with coalescing droplets (VOF) are shown in blue, while runs with non-coalescing droplets (FT) are shown in orange. Runs with coalescing droplets show only small shifts, whereas the runs with non-coalescing, less deformable droplets show significant drag increase. Dependence of the total interface area of the droplets ⟨ ⟩ on the total volume fraction Φ. We have normalised each area by the total initial surface area 0 of the droplets. The VOF runs (blue) show a major reduction in surface area due to coalescence, whereas the FT runs (orange) show a slight increase, due to droplet deformation. Inset: Time history of the total interface area. Each run is plotted according to the colours and markers listed in table I. Note how the coalescing droplets (VOF) reach statistical equilibrium after + ≈ 8000, while the non-coalescing droplets (FT) converge very rapidly because of the absence of topological changes.
cence in simple shear flow in laminar condition by Shardt, Derksen, and Mitra (2013), which show droplet coalescence occurring only below a critical 0 . However, as we shall discuss in the next paragraph, the 0 = 0.1 droplets are more tightly confined in the channel centre than the 0 = 0.05 droplets, thus leading to a higher rate of coalescence. Figure 5 shows how the volume fraction of the dispersed phase depends on the distance from the channel walls. The coalescing droplet profiles (VOF) clearly show a single peak at the channel centre, = : this peak arises as the droplets are driven toward the region of lowest shear ( = ) by a "deformation-induced lift force" (Raffiee, Dabiri, and Ardekani, 2017;Hadikhani et al., 2018;Alghalibi, Rosti, and Brandt, 2019). Confinement in the channel centre leads to coalescence and the formation of large droplets, as seen in figure 1b.
The FT droplets cannot coalesce, and the droplet-droplet interaction produces a volume effect which forces them to spread across the channel: this manifests as an almost flat volume fraction in the region 0.5 < < in figure 5. Also, we see that the volume fraction tends to zero for < = ∕8, as surface tension preserves the droplet radius , and prevents the droplets from fully conforming with the flat channel walls. For all but one of the non-coalescing droplet runs plotted in figure 5, ⟨ ⟩ has a local maximum near the wall, in the region 0.15 < < 0.3 . This phenomenon is due to the "shear-gradient lift force", which is known to act on particles in curved velocity profiles (Ho and Leal, 1974;Martel and Toner, 2014;Hadikhani et al., 2018;Alghalibi, Rosti, and Brandt, 2019). Due to the curvature of the velocity profiles shown in figure 2, the droplets experience different flow velocities on each side, resulting in a lift force toward the wall. From figure 5, we also notice that the more deformable droplets (FT3a, FT5a, and FT10a) produce a maximum which is further from the wall: this is mainly due to (i) an increase in the deformation-induced lift force, and to (ii) a greater elongation of the droplets in the shear direction, producing a wider wall layer.
We are now ready to investigate how droplets affect the turbulent flow, and we start by analysing the second-order statistics of the flow, which tell us how momentum is transferred across different parts of the channel. Figure 6 shows four of the six unique components of the Reynolds stress tensor in wall units ⟨ ′ ′ ⟩ + ≡ ⟨ ′ ′ ⟩ ∕ 2 , with the single-phase (SP0) Reynolds stresses shown in black as reference. The coalescing droplets simulations (VOF) show little change in stresses relative to single-phase flow. Going from single phase to the noncoalescing droplets however, we see a reduction in the stream- FT VOF Φ Φ Figure 5: Dependence of the mean volume fraction of droplets ⟨ ⟩ on the distance from the channel wall. Each run is plotted using the colour and marker listed in table I. The profiles are symmetric about the centreline ( = ), so we have plotted runs with non-coalescing (FT), and coalescing (VOF) droplets on the left and right, respectively. Note that for the runs with coalescing droplets, ⟨ ⟩ peaks in the channel centre, whereas for the non-coalescing droplet runs, ⟨ ⟩ shows a peak near the wall. wise velocity fluctuations ⟨ ′2 ⟩ + , and an increase in the wallnormal ⟨ ′2 ⟩ + and spanwise ⟨ ′2 ⟩ + velocity fluctuations. This shows that the isotropy of the turbulent flow has increased due to the presence of non-coalescing droplets. A similar effect has been observed for particle-laden turbulent channel flows, see e.g. Picano, Breugem, and Brandt (2015), in which particles redistribute energy to a "more isotropic state", inducing an overall drag increase growing with the volume fraction of the dispersed phase. We infer that non-coalescing droplets have a back-reaction on the flow comparable to that of rigid particles, producing an increase in isotropy which correlates with an increase in drag. On the other hand, coalescing droplets produce a weaker back reaction on the flow, which shows little change in isotropy or drag.
When compared to the other components of the Reynolds stresses, the shear stress ⟨ ′ ′ ⟩ + shows only a small change due to the presence of droplets. However, as we shall see next, this shear stress opposes the pressure gradient in the channel, producing a profound impact on the drag. The full shear stress balance for the multiphase problem under investigation can be obtained as follows. We start by taking average of the stream-wise ( = 1) component of equation 1.
In fully developed turbulent channel flows, most of these terms average to zero, and the equation simplifies to where we have moved from the index notation ( 1 , 2 , 3 ) to ( , , ) for the sake of clarity. Hereafter, for brevity we omit the subscripts on angled brackets. Integrating from the wall = 0, to = we obtain The non-penetration boundary conditions at the walls enforce ′ = 0 and with 1 = 0 at the wall, the lower limit of the right hand side is ⟨ , ⟩| =0 = = by the definition of the wall shear stress. We relabel = and obtain  Figure 7: (a) The balance of shear stresses as a function of the distance from the channel wall. The dashed line is the total stress budget. Stresses for run SP0 are shown by solid black lines. The differences between VOF10a and VOF10b stresses are shown in shades of blue, whereas the differences between FT10a and FT10b stresses are shown in shades of orange. We see that + peaks near the wall for the runs with non-coalescing droplets (FT), but is spread across the channel for the coalescing runs (VOF).
The different stress distributions across the channel ultimately lead to different values of drag for coalescing and non-coalescing droplets. (b) Mean shear stresses for all runs. The stacked bars are ⟨ + ⟩ , ⟨ + ⟩ , and ⟨ + ⟩ from bottom to top.
By dividing the equation by , we obtain the following dimensionless expression for the shear stress budget in the channel, where are the viscous, Reynolds, and interfacial shear stresses, respectively. Here, we calculate the viscous stress and Reynolds stress using equations 11 and 12, while the interfacial stress is calculated as the remaining part of the total budget in equation 10 1 . Figure 7a shows the balance of shear stresses from the channel wall ( = 0) to the centre ( = ). In agreement with previous works (Pope, 2000), the single-phase run (SP0) produces a viscous stress + which is highest near the wall where the shear rate is maximum, and a Reynolds stress + which dominates for > 0.1 , where turbulent fluctuations abound.
The coalescing droplet runs (VOF) in figure 7a have an interfacial stress + which peaks around = 0.5 . This stress ⟨ ′ ′ ⟩ −⟨ ⟩⟨ ′ ′ ⟩ for each of the FT runs, and found that the error in shear stress was always less than 3.5% of . occurs due to the droplet interfaces, which resist the deforming effects of turbulent fluctuations, at the detriment of the Reynolds stress. Note that + is larger for the smaller capillary number case (VOF10b compared to VOF10a), because the surface tension coefficient is larger, so surface tension forces are larger.
The non-coalescing droplet runs (FT) in figure 7a, on the other hand, have very little interfacial stress + above > 0.5: instead, the peak of + occurs at roughly the same wall-normal location as the peak in the volume fraction ⟨ ⟩ seen in figure 5. In both figure 5 and figure 7, the peak moves away from the wall when capillary number increases. A similar trend is also observed for the location of the maximum turbulent kinetic energy production (not shown here). The correlation of locations for these three statistics suggests that the "wall layering" and "shear-gradient lift forces discussed above, which produce a peak in ⟨ ⟩ near the channel wall, are also responsible for + generation and kinetic energy generation. The enhanced + close to the wall is compensated in the budget by a reduction in + for the cases of non-coalescing droplets.
The averaged stresses are shown for all runs in figure 7b. The mean stresses are calculated by integrating + , + , and + in the wall-normal direction from 0 to , for example, The averaged form of equation 10 is 0.5 = ⟨ + ⟩ + ⟨ + ⟩ + ⟨ + ⟩ , showing the averaged stresses are also in balance with the wall stress budget. We observe that for coalescing droplets, the dispersed fluid produces an interfacial stress ⟨ + ⟩ which is compensated entirely by a reduction in Reynolds stress ⟨ + ⟩ , with very little change in the viscous stress ⟨ + ⟩ . However, in the case of non-coalescing droplets the increase in interfacial stress ⟨ + ⟩ is compensated by a reduction in both the Reynolds stress ⟨ + ⟩ , and the viscous stress ⟨ + ⟩ . For the single-phase case, the dynamic viscosity is constant throughout the channel, so the mean viscous stress is proportional to the centreline velocity, and hence the variation of ⟨ + ⟩ can be used to quantify drag in the channel, with a larger/smaller ⟨ + ⟩ corresponding to drag reduction/increase. For the multiphase problem, dynamic viscosity is different for the carrier phase and dispersed phases, and we should integrate ⟨ ⟩ ∕ to the centreline, so the relationship between centreline velocity and ⟨ + ⟩ is not exactly linear. However, due to the low volume fraction and low changes of viscosity, we found that considering variation of the material properties ( , ) and variation of the fluid velocity as independent produces only small changes in the averaged statistics. Hence we can still relate the viscous stress to the centreline velocity and thus to the drag changes in the multiphase simulations. Indeed, the three runs with the smallest bulk velocity + in the inset of figure 2 are FT10b, FT5b, and FT3b, and the three runs with the smallest mean viscous stress ⟨ + ⟩ are also FT10b, FT5b, and FT3b (figure 7b). Based on the above discussion, we can now relate the increased drag for non-coalescing droplets to the wall normal location of the droplets: the non-coalescing droplets in runs FT10b, FT5b, and FT3b encroach into the viscous wall region and oppose the shearing flow, reducing the viscous shear stress and thereby increasing drag.

V. CONCLUSIONS
We perform direct numerical simulations of coalescing and non-coalescing droplets in turbulent channel flows to single out the effect of coalescence. Coalescing droplets are simulated using the volume-of-fluid method, and non-coalescing droplets with the front-tracking method. We find that the droplets which are non-coalescing and less deformable produce an increase in drag, whereas the other droplets do not. We explained this by looking at the wall-normal location of droplets in the channel: the coalescing droplets experience a deformation-induced lift force, which drives them away from the shearing flow near the wall, out of the viscous sublayer; this is possible due to the coalescence which allows droplets to accumulate at the centreline. On the other hand, the noncoalescing droplets do not; indeed, non-coalescing droplets roughly behave as particles, uniformly distributing across the channel, forming a wall layer and increasing the isotropy of the flow. In this case, droplets remain in the viscous sublayer, generating an interfacial shear stress, which reduces the budget for viscous shear stress in the channel. From equation 15, we relate a reduction in the viscous shear stress to a reduction in the centreline velocity, and ultimately to an increase in drag.
Our results agree well with the experiments carried out by Descamps et al. (2008), who found that larger bubbles produce less drag; in our study, large droplets are obtained through coalescence, and indeed produce less drag. Our proposed mechanism for drag increase is also similar to that proposed by Dabiri, Lu, and Tryggvason (2013), who showed that less deformable bubbles enter the viscous sublayer, leading to an increase in viscous dissipation and an increase in drag. We offer two main developments. Firstly, we extend the study to coalescing droplets. Secondly, we believe that viscous shear stress is a better predictor of drag than viscous dissipation, as the proportionality between the mean viscous shear stress and centreline velocity (equation 15) is exact for single-phase channel flows, and only slightly affected by the change in material properties. Although we made simulations at a density ratio of ∕ = 50, which is greater than that of oil in water ( ∕ ≈ 1.5), but less than that of air in water ( ∕ ≈ 830), comparison with experimental literature suggests that our current qualitative conclusions still hold for these flows.
Our findings can help to better understand and control multiphase flows in a variety of applications, such as arteries, pipelines or ships. Through numerical experiments, we have been able to fully characterize the effect of coalescence alone, without the interference of other mechanisms which often arise in experiments with surfactants. How these results are affected by surfactant concentrations, will be the topic of future research.