Droplets in homogeneous shear turbulence

Abstract

We simulate the flow of two immiscible and incompressible fluids separated by an interface in a homogeneous turbulent shear flow at a shear Reynolds number equal to 15 200. The viscosity and density of the two fluids are equal, and various surface tensions and initial droplet diameters are considered in the present study. We show that the two-phase flow reaches a statistically stationary turbulent state sustained by a non-zero mean turbulent production rate due to the presence of the mean shear. Compared to single-phase flow, we find that the resulting steady-state conditions exhibit reduced Taylor-microscale Reynolds numbers owing to the presence of the dispersed phase, which acts as a sink of turbulent kinetic energy for the carrier fluid. At steady state, the mean power of surface tension is zero and the turbulent production rate is in balance with the turbulent dissipation rate, with their values being larger than in the reference single-phase case. The interface modifies the energy spectrum by introducing energy at small scales, with the difference from the single-phase case reducing as the Weber number increases. This is caused by both the number of droplets in the domain and the total surface area increasing monotonically with the Weber number. This reflects also in the droplet size distribution, which changes with the Weber number, with the peak of the distribution moving to smaller sizes as the Weber number increases. We show that the Hinze estimate for the maximum droplet size, obtained considering break-up in homogeneous isotropic turbulence, provides an excellent estimate notwithstanding the action of significant coalescence and the presence of a mean shear.

1 Introduction

The understanding of liquid–liquid emulsions is important in many industrial processes e.g. hydrocarbon separation, suspension crystallization and emulsion polymerization. These flows are characterized by density and viscosity ratios of the order of unity (e.g. water and oil mixtures) and a source of agitation (e.g. an impeller) that creates a turbulent two-phase mixture consisting of a dispersed phase of droplets and a continuous phase. The resulting turbulence in the carrier phase is altered directly by the droplet feedback on the surrounding fluid and indirectly by droplet–droplet interactions. Many aspects of the complex interaction of the dispersed phase with the continuous phase are not well understood. In particular, there are questions related to the topological changes and to the role of the surface tension of the dispersed phase, the stationarity of the turbulent statistics and the kinetic energy budget.

Liquid–liquid emulsions have been the subject of numerous experimental (Berkman & Calabrese 1988; Pacek, Man & Nienow 1998; Lovick et al. 2005) and computational studies (Perlekar et al. 2012; Skartlien, Sollum & Schumann 2013; Komrakova, Eskin & Derksen 2015; Scarbolo, Bianco & Soldati 2015; Dodd & Ferrante 2016). The computational studies can be broadly categorized as forced homogeneous isotropic turbulence (Perlekar et al. 2012; Skartlien et al. 2013; Komrakova et al. 2015), decaying homogeneous isotropic turbulence (Dodd & Ferrante 2016) and turbulent wall flows (Scarbolo et al. 2015). Forced homogeneous isotropic turbulence has the advantage of producing a statistically homogeneous and isotropic flow field that, in time, can reach a statistically stationary state. However, in forced homogeneous isotropic turbulence, the turbulent kinetic energy must be induced artificially via a forcing term in the Navier–Stokes equations. This is in contrast to a natural forcing mechanism that produces turbulent kinetic energy from finite Reynolds stresses interacting with a mean velocity gradient. While forcing homogeneous isotropic turbulence may be appropriate for studying the droplet size distributions, it has been argued that artificial forcing is inappropriate for studying two-way coupling effects (Elghobashi 2019). Therefore, for studying the turbulent kinetic energy budget, either decaying isotropic turbulence or turbulent shear flow might be preferable.

In decaying isotropic turbulence, it was shown that the presence of finite-size droplets always enhances the decay rate of the turbulent kinetic energy (Dodd & Ferrante 2016). Also, the deformation, break-up or coalescence of the droplets introduces an additional term to the turbulent kinetic energy equation – the power of the surface tension – termed $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}$ by Dodd & Ferrante (2016), which describes the rate of change of the interfacial energy, balancing the kinetic energy transfer between the external fluid and the flow inside the droplets. Correct identification of these pathways for the turbulent kinetic energy exchange is fundamental to understand the turbulence modulation by the droplets and then to model it.

Building upon previous studies, we consider finite-size bubbles/droplets of Taylor length scale in homogeneous shear turbulence (Tavoularis & Corrsin 1981a , b ; Pumir 1996; Mashayek 1998; Sekimoto, Dong & Jimenez 2016). Homogeneous shear turbulence flow is conceivably the simplest case in which the flow remains statistically homogeneous in all spatial directions. Moreover, compared to forced isotropic turbulence, it has a natural energy production mechanism via a mean velocity gradient. We note that ideal homogeneous shear turbulence is self-similar, implying an unbounded energy growth within infinite domains (Sukheswalla, Vaithianathan & Collins 2013). This condition limits any numerical simulations to relatively short times, concerning only the initial shearing of isotropic turbulence (Rogers & Moin 1987; Lee, Kim & Moin 1990; Sukheswalla et al. 2013). However, as demonstrated by Pumir (1996) and Sekimoto et al. (2016) in single-phase flow, the finite computational box introduces a large-scale confinement effect similar to that enforced by a wall; thus, a meaningful statistically stationary state can be reached over long periods, termed statistically stationary homogeneous shear turbulence (SS–HST). In particular, Sekimoto et al. (2016) showed that long-term simulations of HST are ‘minimal’ in the sense of containing on average only a few large-scale structures: all the one-point statistics agree well with those of the logarithmic layer in turbulent channel flows, particularly when scaled with the friction velocity derived from the measured Reynolds stresses. The same holds for the wall-parallel spectra of the wall-normal velocity. The authors concluded that the similarities between the steady-state homogeneous shear turbulence and other shear flows, particularly with the logarithmic layer of wall turbulence, make it a promising system to study shear turbulence in general. These observations, combined with the insights recently gained in the droplet–turbulence interaction in decaying homogeneous isotropic turbulence, motivate us to further investigate turbulence modulation due to droplets/bubbles in steady-state homogeneous shear turbulence.

In this paper, we present a direct numerical simulation of an emulsion created by droplets dispersed in homogeneous shear turbulence. By changing the initial size of the dispersed phase and the Weber number, we aim to answer the following questions:

1. (i) Can a statistically stationary state be reached when the dispersed phase actively undergoes break-up and coalescence in homogeneous shear turbulence?

2. (ii) If so, what determines the steady-state size distribution of the dispersed phase?

3. (iii) How does the dispersed phase change the turbulent kinetic energy budget?

Homogeneous shear turbulence shares many similarities with other shear flows, including turbulent wall flows (Sekimoto et al. 2016); therefore, by answering these questions, we expect to improve our understanding of the droplet–turbulence interaction and, hopefully, help future modelers gain intuition about more complex conditions.

To capture the complex phenomena accurately in a direct numerical simulation of turbulent two-phase flow, we need a numerical method that is reliable and possess the following properties: (i) discrete mass, momentum and kinetic energy conservation, (ii) ability to handle large jumps in density, (iii) ability to handle complex topologies and separation of scales and (iv) accurate surface tension implementation (Mirjalili, Jain & Dodd 2017). In the present work, we choose to use an algebraic volume of fluid method known as THINC (tangent of hyperbola for interface capturing) method, which is a sharp-interface method. This method is relatively new and has been demonstrated to be as accurate as and also cost effective compared to the well-known geometric volume of fluid methods in canonical test cases (Xie, Ii & Xiao 2014), which makes it a good alternative. However, Mirjalili et al. (2017) indicate that large-scale realistic simulations of turbulent two-phase flows using THINC methods are still lacking in the literature and are crucial to fully evaluate the capabilities of these methods (see Rosti, De Vita & Brandt (2019) for the use of the THINC method for low Reynolds number flows). Hence we choose to use this method in the current study, which will serve as an evaluation of the robustness of THINC methods for complex realistic simulations.

This paper is organized as follows. In § 2, we first discuss the flow configuration and the governing equations and then present the numerical methodology used. The results on the fully developed two-phase homogeneous shear turbulent flow are presented in § 3, where we answer the questions discussed above based on our observations. In particular, we first show how the turbulent flow is modified by the droplets and how the droplets evolve in the turbulent flow, and then explain how these modifications occur by studying the turbulent kinetic energy balance in the two-phase flow. Finally, all the main findings and conclusions are summarized in § 4.

Figure 1. Sketch of the computational domain and of the Cartesian coordinate system. The visualization pertains to flow at $Re_{z}\approx 15\,000$ with $5\,\%$ volume fraction of the dispersed phase at $We_{\unicode[STIX]{x1D706}}\approx 0.75$ . The blue colour is used to depict the surfaces of the droplets.

2 Methodology

2.1 Governing equations and numerical methods

We consider the flow of two immiscible incompressible fluids in a periodic box subject to a uniform mean shear ${\mathcal{S}}$ . Figure 1 shows a sketch of the geometry and the Cartesian coordinate system, where $x$ , $y$ and $z$ ( $x_{1}$ , $x_{2}$ and $x_{3}$ ) denote the streamwise, shear and spanwise coordinates, and $u$ , $v$ and $w$ ( $u_{1}$ , $u_{2}$ and $u_{3}$ ) denote the respective components of the velocity field. Standard periodic conditions are applied in $x$ and $z$ , and a shear-periodic boundary condition is enforced in $y$ , i.e. 

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}(x_{1}+L_{x},x_{2},x_{3})=u_{i}(x_{1},x_{2},x_{3}), & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}(x_{1},x_{2}+L_{y},x_{3})=u_{i}(x_{1}-{\mathcal{S}}tL_{2},x_{2},x_{3}), & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}(x_{1},x_{2},x_{3}+L_{z})=u_{i}(x_{1},x_{2},x_{3}). & \displaystyle\end{eqnarray}$$

The total velocity field $u_{i}$ can be decomposed for convenience into the sum of a mean component $\langle u_{i}\rangle _{xz}$ generated by the imposed shear ${\mathcal{S}}$ , i.e.  $\langle u_{i}\rangle _{xz}={\mathcal{S}}x_{2}\unicode[STIX]{x1D6FF}_{1i}$ where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta, and a fluctuating part $u_{i}^{\prime }$ ( $u_{i}^{\prime }=u_{i}-\langle u_{i}\rangle _{xz}$ ). In this article we indicate the spatial average in the $x$ and $z$ directions with $\langle \cdot \rangle _{xz}$ , fluctuations with the prime symbol ( $^{\prime }$ ), and the average in the full volume with $\langle \cdot \rangle$ . The time evolution of the fluctuating velocity $u_{i}^{\prime }$ is described by

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{j}^{\prime }}{\unicode[STIX]{x2202}x_{j}}}+{\mathcal{S}}x_{2}{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{1}}}+{\mathcal{S}}u_{2}^{\prime }\unicode[STIX]{x1D6FF}_{i1}\right)=-{\displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}}+{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}+f_{i}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{i}}}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the fluid density, $p$ is the pressure, $\unicode[STIX]{x1D70F}_{ij}=2\unicode[STIX]{x1D707}\unicode[STIX]{x1D60B}_{ij}$ with $\unicode[STIX]{x1D707}$ the dynamic viscosity and $\unicode[STIX]{x1D60B}_{ij}$ the strain rate tensor ( $\unicode[STIX]{x1D60B}_{ij}=(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})/2$ ) and $f_{i}$ is the surface tension force defined as $f_{i}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}n_{i}\unicode[STIX]{x1D6FF}$ , where $\unicode[STIX]{x1D6FF}$ is the Dirac delta function at the interface, $\unicode[STIX]{x1D70E}$ the interfacial surface tension, $\unicode[STIX]{x1D705}$ the interface curvature and $n_{i}$ the normal to the interface. This equation is written in the so-called one-fluid formulation (Tryggvason, Sussman & Hussaini 2007) so that only one set of equations is solved in both phases. The problem is solved by introducing an indicator function $H$ to identify each fluid phase so that $H=1$ in the region occupied by the suspended dispersed fluid (fluid $1$ ) and $H=0$ in the carrier phase (fluid $2$ ). Considering that both fluids are transported by the flow velocity, we update $H$ in the Eulerian framework by the following advection equation written in divergence form,

(2.6) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}H}{\unicode[STIX]{x2202}x_{i}}}=\unicode[STIX]{x1D719}{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ is the cell-averaged value of the indicator function.

The above governing equations are solved numerically. First, the transport equation for $\unicode[STIX]{x1D719}$ is updated following the methodology described by Ii et al. (2012) and Rosti et al. (2019) in order to obtain $\unicode[STIX]{x1D719}^{n+1}$ , which is used to update the density and viscosity of the fluids. In particular, the indicator function $H$ is approximated as

(2.7) $$\begin{eqnarray}\displaystyle H(X,Y,Z)\approx \widehat{H}(X,Y,Z)={\textstyle \frac{1}{2}}(1+\tanh (\unicode[STIX]{x1D6FD}(P(X,Y,Z)+d))), & & \displaystyle\end{eqnarray}$$

where $X,Y,Z\in [0,1]$ is a centred local coordinate system defined in each cell, $P$ is a three-dimensional quadratic curved surface function determined algebraically by imposing the correct value of the three normal components and the six components of the Cartesian curvature tensor in each cell, $d$ is a normalization parameter used to enforce that the integral of the indicator function in each cell equals $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6FD}$ is a sharpness parameter; $\unicode[STIX]{x1D6FD}$ is set equal to $1$ in the present work, the smallest value allowed by the method which ensures the sharpest possible interface for a given mesh size. Second, the momentum equation and the incompressibility constraint are solved following the method proposed by Gerz, Schumann & Elghobashi (1989) and recently adopted by Tanaka (2017), in which the third term on the left-hand side of the momentum equation (2.4), i.e. the advection due to the mean shear flow, is solved separately using a Fourier approximation. In particular, the second-order Adams–Bashforth method is applied for the convection and viscous terms in (2.4) to obtain an intermediate velocity

(2.8) $$\begin{eqnarray}\displaystyle {u_{i}^{\prime }}^{\ast }={u_{i}^{\prime }}^{n}+\unicode[STIX]{x0394}t\left({\textstyle \frac{3}{2}}\text{rhs}_{i}^{n}-{\textstyle \frac{1}{2}}\text{rhs}_{i}^{n-1}\right), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x0394}t$ is the time step from time $t^{n}$ to $t^{n+1}$ and

(2.9) $$\begin{eqnarray}\displaystyle \text{rhs}_{i}=-{\mathcal{S}}u_{2}^{\prime }\unicode[STIX]{x1D6FF}_{i1}-{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{j}^{\prime }}{\unicode[STIX]{x2202}x_{j}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}. & & \displaystyle\end{eqnarray}$$

The time step $\unicode[STIX]{x0394}t$ is chosen such that the Courant–Friedrichs–Lewy number $U_{max}\unicode[STIX]{x0394}t/\unicode[STIX]{x0394}x$ is smaller than unity, where $U_{max}={\mathcal{S}}L_{y}$ is the maximum of the mean shear flow velocity inside the computational domain. The advection due to the mean shear flow is then solved separately using a Fourier approximation as

(2.10) $$\begin{eqnarray}\displaystyle {u_{i}^{\prime }}^{\ast \ast }(x_{1},x_{2},x_{3})={u_{i}^{\prime }}^{\ast }(x_{1}-\unicode[STIX]{x0394}t{\mathcal{S}}x_{2},x_{2},x_{3}). & & \displaystyle\end{eqnarray}$$

Note that Tanaka (2017) modified the approach of Gerz et al. (1989) by performing a similar additional step for the pressure. Our tests suggest that the original form by Gerz et al. (1989) is numerically more stable and physically consistent with the incompressibility constraint because the pressure is not a transported quantity. The surface tension term $f_{i}$ is then taken into account by updating the velocity field: we use the continuum surface force model by Brackbill, Kothe & Zemach (1992) to compute the surface tension force where the normals are obtained with the well-known Youngs approach (Youngs 1982), i.e.  $f_{i}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}x_{i}$ , thus giving

(2.11) $$\begin{eqnarray}\displaystyle {u_{i}^{\prime }}^{\ast \ast \ast }={u_{i}^{\prime }}^{\ast \ast }+\unicode[STIX]{x0394}t{\displaystyle \frac{f_{i}^{n+1}}{\unicode[STIX]{x1D70C}}}. & & \displaystyle\end{eqnarray}$$

Then, we enforce the zero divergence of the velocity field by solving the following Poisson equation

(2.12) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}^{2}p^{n+1}}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}}={\displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x0394}t}}{\displaystyle \frac{\unicode[STIX]{x2202}{u_{i}^{\prime }}^{\ast \ast \ast }}{\unicode[STIX]{x2202}x_{i}}}, & & \displaystyle\end{eqnarray}$$

which is solved with a standard fast-Fourier-transform-based solver by exploiting the periodic and shear-periodic boundary conditions as detailed in Tanaka (2017). Finally, we correct the velocity with $p^{n+1}$ to enforce the incompressibility constraint

(2.13) $$\begin{eqnarray}\displaystyle {u_{i}^{\prime }}^{n+1}={u_{i}^{\prime }}^{\ast \ast \ast }-\unicode[STIX]{x0394}t{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}p^{n+1}}{\unicode[STIX]{x2202}x_{i}}}. & & \displaystyle\end{eqnarray}$$

Note that our numerical scheme discretely conserves both momentum and kinetic energy (in the absence of viscosity and surface tension) since we use a second-order centred finite difference on a staggered mesh and the divergence form of the convective terms (Morinishi et al. 1998).

2.2 Set-up

The problem is governed by several dimensionless parameters, which define the problem under consideration. First, the computational box is defined by two aspect ratios $A_{xz}=L_{x}/L_{z}$ and $A_{yz}=L_{y}/L_{z}$ which are fixed equal to $2.05$ and $1.025$ , respectively. These values have been chosen according to what was proposed by Sekimoto et al. (2016) as ‘acceptable’ in the sense that they fall within the range of parameters in which the flow is as free as possible from box effects and can thus be used as a model of shear-driven turbulence in general. Indeed, homogeneous shear turbulence in an infinite domain evolves towards larger and larger length scales while simulations in a finite box are necessarily constrained to some degree by the box geometry. These authors noticed that the effect of the geometry can be reduced by ensuring that $L_{z}$ is the main constraint, thus resulting in the flow being ‘minimal’ in the spanwise direction. Next, once the size of the numerical box is fixed, to fully characterize the problem we define the shear Reynolds number based on the box width

(2.14) $$\begin{eqnarray}\displaystyle Re_{z}={\displaystyle \frac{{\mathcal{S}}L_{z}^{2}}{\unicode[STIX]{x1D708}}}, & & \displaystyle\end{eqnarray}$$

the Weber number based on the initial droplet diameter $D_{0}$

(2.15) $$\begin{eqnarray}\displaystyle We_{{\mathcal{S}}_{0}}={\displaystyle \frac{\unicode[STIX]{x1D70C}{\mathcal{S}}^{2}D_{0}^{3}}{\unicode[STIX]{x1D70E}}}, & & \displaystyle\end{eqnarray}$$

and the ratio of the initial droplet diameter to the box size $A_{Dz}=D_{0}/L_{z}$ . In the following, we consider one case of single-phase flow as reference and nine cases of two-phase flows, all at the same Reynolds number equal to $15\,200$ ; in the multiphase cases, we vary the ratio $A_{Dz}$ and $We_{{\mathcal{S}}_{0}}$ , as summarized in table 1. Note that the Weber number here is mainly determined by the interfacial surface tension $\unicode[STIX]{x1D70E}$ . Two other non-dimensional parameters are the density and viscosity ratios, which are fixed to unity to study the individual effect of the Weber number (interfacial surface tension).

Table 1. Summary of the direct numerical simulations performed with different initial droplet sizes $D_{0}$ , numbers of droplets ${\mathcal{N}}_{0}$ and surface tension $\unicode[STIX]{x1D70E}$ , all at a fixed Reynolds number $Re_{z}=15\,200$ and volume fraction $\unicode[STIX]{x1D6F7}=5\,\%$ . r.m.s., root mean square.

Besides the parameters just defined and based on the geometrical dimensions and initial and boundary conditions alone, in the following discussion we will use other non-dimensional numbers because they turned out to be more relevant to understand the problem at hand; in particular, the two non-dimensional parameters which characterize the single-phase homogeneous shear turbulent flows, the Taylor-microscale Reynolds number $Re_{\unicode[STIX]{x1D706}}$ and the shear-rate parameter $S^{\ast }$ , defined as

(2.16) $$\begin{eqnarray}\displaystyle Re_{\unicode[STIX]{x1D706}}=\left({\displaystyle \frac{2{\mathcal{K}}}{3}}\right)^{1/2}{\displaystyle \frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D708}}}=\left({\displaystyle \frac{5}{3\unicode[STIX]{x1D708}\unicode[STIX]{x1D700}}}\right)^{1/2}2{\mathcal{K}}, & & \displaystyle\end{eqnarray}$$

and

(2.17) $$\begin{eqnarray}\displaystyle S^{\ast }={\displaystyle \frac{2{\mathcal{S}}{\mathcal{K}}}{\unicode[STIX]{x1D700}}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}=\sqrt{10\unicode[STIX]{x1D708}{\mathcal{K}}/\unicode[STIX]{x1D700}}$ is the Taylor microscale (Sekimoto et al. 2016), ${\mathcal{K}}=\langle \unicode[STIX]{x1D70C}u_{i}^{\prime }u_{i}^{\prime }\rangle /2$ is the turbulent kinetic energy per unit volume and $\unicode[STIX]{x1D700}=\unicode[STIX]{x1D707}\langle \unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}\rangle$ is the dissipation rate of the fluctuating energy. These two non-dimensional numbers can be interpreted as the ratio of the eddy-turnover time $\unicode[STIX]{x1D70F}_{0}=(2{\mathcal{K}})^{1/2}/\unicode[STIX]{x1D700}$ and the Kolmogorov time scale $\unicode[STIX]{x1D70F}_{K}=(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D700})^{1/2}$ and the mean shear time scale $\unicode[STIX]{x1D70F}_{{\mathcal{S}}}=1/{\mathcal{S}}$ , respectively.

Figure 2. The ratio of the two Weber numbers introduced here, one based on the mean shear $We_{{\mathcal{S}}_{0}}$ and one on the velocity fluctuations $We_{rms_{0}}$ , as a function of the Weber number based on the initial droplet size, $We_{0}$ . The circle, square and triangle symbols are used to distinguish cases with different surface tension but same ratio $We_{{\mathcal{S}}_{0}}/We_{rms_{0}}$ , while the brown, green and blue colours represent cases with the ratio $We_{{\mathcal{S}}_{0}}/We_{rms_{0}}$ equal to $1/5$ , $1$ and $5$ , respectively. These symbols and colour scheme will be used throughout the rest of the paper.

Weber numbers can be defined in several ways. In (2.15) we defined the Weber number based on the mean shear, but it can also be defined based on the velocity fluctuations, thus giving

(2.18) $$\begin{eqnarray}\displaystyle We_{rms_{0}}={\displaystyle \frac{2\unicode[STIX]{x1D70C}{\mathcal{K}}D_{0}}{\unicode[STIX]{x1D70E}}}. & & \displaystyle\end{eqnarray}$$

Note that the latter definition is the one usually used in homogeneous isotropic turbulent flows in the absence of a mean flow (Dodd & Ferrante 2016). Both the Weber numbers $We_{{\mathcal{S}}_{0}}$ and $We_{rms_{0}}$ are of interest since they are based on two different mechanisms that may affect the droplets dynamics: on large scales (large droplets) the effect of the mean shear is dominant, while on small scales (small droplets) the flow is mainly dominated by the isotropic turbulent fluctuations. Our set of parameters is chosen such that the ratio of these two Weber numbers $We_{{\mathcal{S}}_{0}}/We_{rms_{0}}$ equals $1/5$ , $1$ and $5$ , as reported in figure 2. In general both the mechanisms are present together and hence we can define a Weber number which incorporates both effects as

(2.19) $$\begin{eqnarray}\displaystyle We_{0}={\displaystyle \frac{\unicode[STIX]{x1D70C}(\sqrt{2{\mathcal{K}}}+{\mathcal{S}}D_{0})^{2}D_{0}}{\unicode[STIX]{x1D70E}}}. & & \displaystyle\end{eqnarray}$$

Finally, we can define a Weber number based on $\unicode[STIX]{x1D706}$ as

(2.20) $$\begin{eqnarray}\displaystyle We_{\unicode[STIX]{x1D706}}={\displaystyle \frac{\unicode[STIX]{x1D70C}(\sqrt{2{\mathcal{K}}}+{\mathcal{S}}\unicode[STIX]{x1D706})^{2}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D70E}}}. & & \displaystyle\end{eqnarray}$$

The choice of using $\unicode[STIX]{x1D706}$ in the definition of the Weber number instead of a dimension associated with the suspended phase is due to the fact that the interface is not only deforming, thus losing its original spherical shape, but also actively undergoing merging and break-up processes, which makes the definition of a unique dimension difficult. Therefore, we propose to rely on a fluid length scale, which, as shown below in the results, yields a good collapse of our data. In the following discussion, we use $We_{\unicode[STIX]{x1D706}}$ to discuss the results; the value of $We_{0}$ is reported in order to fully characterize the initial conditions of the present simulations.

2.3 Code validation

The numerical code used in this work has been extensively validated in the past for multiphase turbulent flows simulations (Rosti & Brandt 2017; Rosti et al. 2018a , b ). Here, we provide one more comparison with literature results for the specific case of HST. The single-phase homogeneous shear turbulence has been validated by reproducing one of the cases investigated by Pumir (1996); in particular, we simulated the Run No.  $2$ in that paper. The initial condition at $t=0$ is a homogeneous isotropic turbulent field at $Re_{\unicode[STIX]{x1D706}}=50.8$ , obtained in a square computational box of size $2\unicode[STIX]{x03C0}$ discretized with $256$ grid points in each direction. From the time history of the turbulent kinetic energy ${\mathcal{K}}$ and of the enstrophy $\unicode[STIX]{x1D6FA}=\langle \unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}\rangle$ , shown in figure 3(a), we observe a first transient phase for $0\leqslant t{\mathcal{S}}\leqslant 30$ , where the kinetic energy and enstrophy grow rapidly, followed by a statistically stationary state characterized by a cyclic succession of turbulent kinetic energy peaks rapidly followed by a peak in enstrophy with a time lag of approximately $5{\mathcal{S}}$ . This behaviour is well captured in our simulation. A quantitative validation is performed first by comparing the mean components of the velocity anisotropy tensor, $\unicode[STIX]{x1D623}_{ij}=\langle u_{i}u_{j}/u_{k}^{\prime }u_{k}^{\prime }-\unicode[STIX]{x1D6FF}_{ij}/3\rangle$ computed in our simulations ( $b_{11}=0.231$ , $b_{22}=0.129$ , $b_{33}=0.101$ , $b_{12}=0.147$ ) with the data reported by Pumir (1996), and we found that the differences are below $5\,\%$ . A further comparison is shown in figure 3(b) where the normalized histograms of the streamwise and shear components of the velocity obtained with the present simulations are compared with the results reported in the literature (Pumir 1996); again we observe a very good agreement.

Figure 3. (a) Time history of the turbulent kinetic energy ${\mathcal{K}}=\langle u_{i}^{\prime }u_{i}^{\prime }\rangle /2$ (black line) and enstrophy $\unicode[STIX]{x1D6FA}$ (grey line), normalized by their mean values. (b) Normalized histogram of the streamwise (red) and shear (orange) components of the velocity fluctuations, $u^{\prime }$ and $v^{\prime }$ . The lines and symbols are used to distinguish our results (lines) from those by Pumir (1996).

3 Results

3.1 Statistically stationary state

We start our analysis by considering the single-phase flow at $Re_{z}=15\,200$ . The problem is solved numerically on a computational mesh of $1312\times 640\times 624$ grid points and the simulation is run for approximately $250{\mathcal{S}}$ time units. Note that the grid spacing is chosen sufficiently small for good resolution of the smallest turbulent scales as indicated by $\unicode[STIX]{x0394}x/\unicode[STIX]{x1D702}\approx 0.7$ , where $\unicode[STIX]{x1D702}$ is the Kolmogorov scale defined as $\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C})^{3/4}/\unicode[STIX]{x1D700}^{1/4}$ . The initial flow field is fully developed single-phase homogeneous isotropic turbulence, and the mean shear ${\mathcal{S}}$ is applied from $t=0$ . As shown in figure 4(a), once the shear is applied, the flow undergoes an initial transient characterized by a strong increase in the production of turbulent kinetic energy, which is not in balance with the dissipation rate. After some time, however, the turbulent kinetic energy ${\mathcal{K}}$ decreases owing to an increase in the dissipation, reaching a new statistically steady state where, on average, the production balances the dissipation ( ${\mathcal{P}}\approx \unicode[STIX]{x1D700}$ ). This state, called steady-state shear turbulence, was first found and characterized by Pumir (1996) and later investigated by others (e.g. Sekimoto et al. 2016). The resulting Taylor-microscale Reynolds number at the steady state is equal to $Re_{\unicode[STIX]{x1D706}}\approx 145$ with the averaged spectrum of the turbulent kinetic energy reported in figure 4(b). Owing to the high Reynolds number, a clear $k^{-5/3}$ regime develops at intermediate scales. We also observe that the spectra of each individual component of the velocity are different at small wavenumbers because of the large-scale anisotropy, while all spectra coincide at higher wavenumbers, consistently with what was observed by Pumir (1996).

Figure 4. (a) Time history of the ratio between the turbulent production ${\mathcal{P}}=-\langle u^{\prime }v^{\prime }\rangle \text{d}\langle u\rangle /\text{d}y$ and the turbulent dissipation rate $\unicode[STIX]{x1D700}=\unicode[STIX]{x1D707}\langle \unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}\rangle$ . The black and green lines represent the single and multiphase flows ( $D_{0}=0.16L_{z}$ and $We_{\unicode[STIX]{x1D706}}\approx 0.75$ ), respectively. (b) Spectrum of the mean turbulent kinetic energy (black solid line) and its three spatial components (black dashed, dotted, and dashed-dotted lines) for the single-phase flow. The other three coloured solid lines (blue, green and brown) are used for the spectra of the two-phase flows with $We_{\unicode[STIX]{x1D706}}=0.02$ , $0.75$ and $5$ . The grey line is $\propto k^{-5/3}$ , and the three vertical dashed lines represent the initial size of the droplets. The spectra are normalized by multiplying by $\unicode[STIX]{x1D700}^{-2/3}$ .

We now consider the multiphase problem. After around $100{\mathcal{S}}$ , when the single-phase flow has already reached a statistically steady state, we inject spherical droplets into the domain at random locations, globally enclosing a volume fraction of the carrier phase of $5\,\%$ . The initial droplet diameter $D_{0}$ is in the inertial range, as shown in figure 4(b) with the vertical dashed lines. In particular, three different initial diameters are chosen, $D_{0}/L_{z}\approx 0.08$ (brown), $0.16$ (green) and $0.32$ (blue), corresponding to approximately $1.1$ , $2.5$ and $5.6$ times the single-phase Taylor microscale $\unicode[STIX]{x1D706}$ . After the introduction of the dispersed phase, a new short transient arises lasting approximately $50{\mathcal{S}}$ , eventually leading to a new statistically steady state, as depicted in figure 4(a). Also, in the multiphase case, we observe that, at steady state, the turbulent production balances on average the dissipation rate ( ${\mathcal{P}}\approx \unicode[STIX]{x1D700}$ ).

The presence of the droplets modifies the flow profoundly. The averaged spectrum of the turbulent kinetic energy in both phases in the two-phase case is reported in figure 4(b), where we observe that the interface mostly affects the large wavenumbers (small scales) for which higher levels of energy are evident, while slightly lower energy is present at the large scales. Note that the result is analogous to what was observed in decaying homogeneous isotropic turbulence for solid particles (Lucci, Ferrante & Elghobashi 2010) and bubbles (Dodd & Ferrante 2016); the increased energy at high wavenumbers has been explained by the break-up of large eddies due to the presence of the suspended phase and the consequent creation of new eddies of smaller scale. In the same figure we can also observe that the effect of the droplets decreases as the Weber number increases; in other words, the spectra of the multiphase cases approach the single phase one as $We$ increases, while for low $We$ the spectra depart from the single-phase case at smaller and smaller wavenumbers.

Figure 5. (a) Weber numbers based on the Taylor microscale, $We_{\unicode[STIX]{x1D706}}$ , as a function of the initial Weber number $We_{0}$ and (b) Reynolds numbers based on the Taylor microscale, $Re_{\unicode[STIX]{x1D706}}$ , as a function of $We_{\unicode[STIX]{x1D706}}$ . The grey solid line in (a) is a fit to our data in the form of $We_{\unicode[STIX]{x1D706}}\propto We_{0}^{2}$ , while the grey solid line in (b) represents the Taylor-microscale Reynolds number $Re_{\unicode[STIX]{x1D706}}$ of the single-phase flow.

As already discussed above, $We_{0}$ is the Weber number based on the initial droplet size, but since the droplets break-up or coalesce, this measure is not fully representative of the final state of the multiphase problem; because of this, in the following sections we prefer to use the Weber number based on a flow length scale, $We_{\unicode[STIX]{x1D706}}$ , reported in figure 5(a) as a function of $We_{0}$ . We can observe that the two Weber numbers are well correlated, with $We_{\unicode[STIX]{x1D706}}$ scaling approximately as the square of $We_{0}$ , i.e.  $We_{\unicode[STIX]{x1D706}}\propto We_{0}^{2}$ . The good level of correlation between the two definitions is a further demonstration that for the parameter range considered here the Weber number variations are mainly due to the changes of the interfacial surface tension rather than the chosen length scale.

We quantify the turbulence modulation by examining the resulting $Re_{\unicode[STIX]{x1D706}}$ , shown for all our simulations in figure 5(b) as a function of the Weber number based on the Taylor microscale $We_{\unicode[STIX]{x1D706}}$ , and also reported in table 1. We can observe that the Reynolds number grows with $We_{\unicode[STIX]{x1D706}}$ and that all the two-phase flow cases exhibit lower Taylor-microscale Reynolds numbers than the single-phase case. Moreover, we observe that the difference decreases as the Weber number increases, with the two-phase flow cases approaching the single-phase one as $We_{\unicode[STIX]{x1D706}}$ increases, consistently with what was already observed in figure 4(b). Indeed, the Reynolds number for the case with the most rigid droplets ( $We_{\unicode[STIX]{x1D706}}\approx 0.02$ ) is approximately half the single phase value ( $-41\,\%$ ), while the difference with the single-phase flow is only $2\,\%$ in the most deformable case ( $We_{\unicode[STIX]{x1D706}}\approx 13$ ). Note that, in the context of unbounded forced turbulent flows, such as homogeneous isotropic turbulence and homogeneous shear turbulence, a reduction of the Reynolds number can be interpreted as a drag increase, contrary to what is usually found in wall-bounded flows with constant flow rates where a reduction in the friction Reynolds number leads to drag decrease.

As a first noteworthy result, the above data demonstrate that a statistical stationarity is not unique to single-phase homogeneous shear turbulent flows, but it is also realizable in the presence of a second, dispersed phase. Here, we have defined the stationary state in terms of the statistical properties of the flow averaged over both phases, but since the droplets can also break-up or coalesce, it is natural to ask what the steady-state size distributions are and how that relates to the turbulence features. These questions are answered in the following sections.

3.2 Size distribution

We now study the transient and steady-state property of the interface separating the two fluids. Figure 6 shows instantaneous snapshots of the two-phase flow at the statistically steady state, which is characterized by droplets with different sizes and shapes: in general we can observe that small droplets are approximately spherical, while the largest ones have very anisotropic shapes and show a preferential alignment with the direction of the mean shear. Also, as the Weber number decreases, the droplets size increases and larger droplets can sustain the spherical shape.

Figure 6. Visualization in the $x{-}y$ plane of the interface in the homogeneous shear turbulent flow for different $We_{\unicode[STIX]{x1D706}}$ : (a) $We_{\unicode[STIX]{x1D706}}\approx 0.02$ , (b) $0.08$ , (c) $0.8$ , (d) $4$ , (e) $5$ and (f) $13$ . In the panels the flow is from left to right.

Figure 7. (a) Time history of the number of droplets ${\mathcal{N}}$ in the domain for different Weber numbers. The rhombus symbols at $t=0$ represent the initial number of droplets. (b) The mean number of droplets ${\mathcal{N}}_{s}$ at the statistically steady state as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ . The grey solid line in (b) is a fit to our data in the form of ${\mathcal{N}}_{s}\propto We_{\unicode[STIX]{x1D706}}$ .

Figure 7(a) shows the temporal evolution of the number of droplets ( ${\mathcal{N}}$ ) under various $We_{\unicode[STIX]{x1D706}}$ and initial sizes $D_{0}$ . The counting of the droplets is conducted automatically by checking the connectivity of the local volume of fluid field ( $\unicode[STIX]{x1D719}$ ) using a $n$ -dimensional image processing library (scipy.ndimage, https://docs.scipy.org/doc/scipy/reference/ndimage.html). We observe that ${\mathcal{N}}$ has an initial transient phase of same duration as the fluid transient phase observed previously in figure 4(a) ( $t{\mathcal{S}}\lesssim 50$ ), before the droplets count approaches a statistically steady value for all the cases considered, consistently with the statistically stationarity of the averaged global flow quantities. Note that the final state is a statistically steady state since the number of droplets ${\mathcal{N}}$ is not constant but continuously varies and oscillates around a mean value, denoted later on as ${\mathcal{N}}_{s}$ . From the figure we observe also that the initial transient phase differs among the cases, with three distinct behaviours evident: (i) in most of the cases, ${\mathcal{N}}$ increases rapidly after the injection (within $t{\mathcal{S}}\approx 10$ ); however, the growth slows down and ${\mathcal{N}}$ reaches its final steady-state value almost monotonically; (ii) cases $4$ and $5$ exhibit a significant overshoot of the number of droplets ${\mathcal{N}}$ for short times before ${\mathcal{N}}$ reduces to the final regime values due to the coalescence; and (iii) case $3$ shows an initial decrease of the number of droplets followed by an increase. Notwithstanding the different behaviours, in all the cases the final number of droplets is always larger than the initial one.

The steady-state value of the number of droplets ${\mathcal{N}}_{s}$ as a function of $We_{\unicode[STIX]{x1D706}}$ is reported in figure 7(b); we observe that ${\mathcal{N}}_{s}$ grows monotonically with $We_{\unicode[STIX]{x1D706}}$ (see also the visualizations in figure 6) and that the growth is nearly linear over the three decades spanned in the present study, i.e. a fit to our data produces ${\mathcal{N}}_{s}\propto We_{\unicode[STIX]{x1D706}}$ with an exponent of $1$ . Since a high Weber number corresponds to a low surface energy, we conjecture that ${\mathcal{N}}_{s}$ grows indefinitely with $We_{\unicode[STIX]{x1D706}}$ . Note also that cases $5$ and $6$ which have different initial droplet diameters, have almost the same final count of droplets ${\mathcal{N}}_{s}$ as well as $We_{\unicode[STIX]{x1D706}}$ . This provides additional evidence that the droplet statistics are better defined by the Weber number $We_{\unicode[STIX]{x1D706}}$ based on the flow quantities rather than by that based on the initial droplet size $We_{0}$ . These results suggest that the relative strength between the break-up and coalescence reflects the history of the flow features, and at equilibrium measurable quantities depend only on the global physical parameters.

Figure 8. (a) Normalized cumulative volume distributions ${\mathcal{V}}/{\mathcal{V}}_{tot}$ of the dispersed phase at the steady state as a function of the equivalent spherical droplet diameters $D$ . The horizontal grey line corresponds to the level ${\mathcal{V}}=0.95{\mathcal{V}}_{tot}$ . (b,c) Contour of the temporal evolution of the normalized cumulative volume distributions of the dispersed phase as a function of the equivalent spherical droplet diameter for cases $4$ (b) and $5$ (c).

Next, we aim to characterize the steady-state size distribution of the emulsion. Thus, we first examine the cumulative volume, ${\mathcal{V}}$ , as a function of the equivalent spherical diameter $D$ defined as the diameter of the sphere occupying the same volume, see figure 8(a). Specifically, figure 8(a) shows the steady-state distributions of all cases, where each point on the curves represents the total volume of the droplets with equivalent diameter lower than $D$ . In the figure, both ${\mathcal{V}}$ and $D$ are normalized by the global maximal values so that the curves are bounded uniformly from above by  $1$ . The figure shows that the cumulative volume distribution only has one inflection point ( $\text{d}^{2}{\mathcal{V}}/\text{d}D^{2}=0$ ), thus indicating that the probability density plot ( $\text{d}{\mathcal{V}}/\text{d}D$ ) is single peaked. In figure 8(a) the Weber number $We_{\unicode[STIX]{x1D706}}$ grows from right to left, as indicated by the list of symbols, suggesting that small droplets tend to be more common at high Weber numbers. Additionally, the range of the droplet diameters also narrows with increasing $We_{\unicode[STIX]{x1D706}}$ , since the cumulative volume grows faster to unity, as visually confirmed in figure 6. Case $2$ , blue line with circle, is the only simulation exhibiting a double peak (i.e.  $\text{d}{\mathcal{V}}/\text{d}D$ has two local maxima): this is due to the presence of very small droplets together with few large ones, as can be seen in figure 6(a). Nevertheless, the overall trend of decreasing size for increasing Weber number is still consistent with the linear scaling between ${\mathcal{N}}_{s}$ and $We_{\unicode[STIX]{x1D706}}$ , as already observed in figure 7(b). Figures 8(a) and 8(c), are contours of ${\mathcal{V}}/{\mathcal{V}}_{tot}$ as a function of the equivalent diameter $D$ and time, and can thus be interpreted as a cumulative spectrogram with most of the droplets centred in the region where the gradient of the colour is the largest. In particular, we selected two specific cases, with same initial Weber number $We_{0}\approx 2$ , but different initial droplet size and surface tension, thus leading to different $We_{\unicode[STIX]{x1D706}}$ . The two figures show the transient behaviour for cases $4$ and $5$ , respectively: in figure 8(b) the mean size distribution remains relatively unchanged over time but it is subject to strong fluctuations, while figure 8(c) shows a clear shift of the population from large droplets to small ones, with a statistically steady state characterized by small fluctuations.

Another important parameter related to the size distribution is the largest droplet size, $D_{max}$ . Assuming break-up of droplets due to the dynamic pressure ( ${\sim}\unicode[STIX]{x1D70C}U^{2}$ ), Hinze (1955) proposed that the largest possible droplet in a turbulent emulsifier is determined by the velocity fluctuation across $D_{max}$ , i.e. one can define a critical Weber number $We_{crit}=\unicode[STIX]{x1D70C}\overline{u^{\prime 2}}D_{max}/\unicode[STIX]{x1D70E}$ , above which the droplet breaks up. Hinze (1955) showed that simple dimensional analysis leads to $D_{max}\propto \unicode[STIX]{x1D700}^{-2/5}$ , if isotropy prevails and the scaling by Kolmogorov (1941) is assumed valid. $D_{max}$ can be in general approximated by the diameter of the equivalent droplet occupying $95\,\%$ of the total dispersed volume, i.e.  $D_{max}\approx D_{95}$ , which is represented in figure 8(a) with the dashed grey line. The symbols in the same figure provide the values of $D_{95}$ for our data. Figure 9(a) shows the normalized $D_{95}$ as a function of the scaled energy input, and indeed we can observe that our data scale with an approximately $-2/5$ slope. We remark that, although Hinze developed his theory considering only isotropic turbulent flows dominated by the break-up process and neglected the coalescence, he hypothesized that the same scaling law might still hold for non-isotropic flows provided that the droplet sizes fall within the inertial range, such as in all our cases. More importantly, the success of the Hinze theory relies on the central assumption that break-up results from the dynamic pressure force, corresponding to a fixed critical Weber number. This is clearly shown in figure 9(b), which shows the Weber number based on $D_{95}$ as a function of $We_{\unicode[STIX]{x1D706}}$ . For all our cases, we obtain that $We_{crit}\approx 1$ . Our results thus confirm that the $-2/5$ scaling between the maximum droplet diameter and the turbulence dissipation applies not only to isotropic turbulence, but also to the homogeneous shear turbulence that we have analysed.

Figure 9. (a) Normalized maximum droplet size $D_{95}$ as a function of the energy input $\unicode[STIX]{x1D700}$ . The grey solid line is the relation $\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}D_{95}/\unicode[STIX]{x1D707}^{2}=0.725(\unicode[STIX]{x1D707}^{5}\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}^{4})^{-2/5}$ proposed by Hinze (1955). (b) Critical Weber number $We_{D_{95}}$ based on the maximum droplet size $D_{95}$ for all the cases considered.

Figure 10. Total interfacial area ${\mathcal{A}}$ as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ . The grey solid line is a fit to our data in the form of ${\mathcal{A}}\propto We_{\unicode[STIX]{x1D706}}^{1/3}$ .

Finally, we can further characterize the size distribution of the emulsion by inspecting the total surface area ${\mathcal{A}}$ of the dispersed phase. This quantity is very important when studying multiphase flows with interfaces, since the rate of work due to the surface tension is equal to the product of the surface tension coefficient and the rate of change in interfacial surface area (Dodd & Ferrante 2016); also, for many industrial applications, the total surface area is often the most important parameter as surfactants tend to reside on the interface or it determines the chemical reaction rate. Figure 10 reports the steady state surface area ${\mathcal{A}}$ as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ and clearly shows that the surface area increases monotonically with the Weber number. As we have shown above, ${\mathcal{N}}\propto We_{\unicode[STIX]{x1D706}}$ , combined with mass conservation, i.e.  ${\mathcal{N}}D^{3}\propto 1$ , leads to the following relation for the total area: ${\mathcal{A}}\propto {\mathcal{N}}D^{2}\propto We_{\unicode[STIX]{x1D706}}^{1/3}$ . In other words, the surface area of the droplets shall also increase with the Weber number defined by the Taylor length of the flow, with a slope of $1/3$ . Figure 10 verifies this scaling. We remark that in the derivation above, we have assumed that the droplets are spherical, which is not always true in our cases. However, provided the linear scaling between ${\mathcal{N}}$ and $We_{\unicode[STIX]{x1D706}}$ remains valid, we expect the $1/3$ scaling law to hold for a wide range of emulsions.

3.3 Turbulent kinetic energy budget

We now study how the multiphase nature of the problem affects the turbulent kinetic energy. To do so, we derive the turbulent kinetic energy evolution equation by first multiplying the momentum conservation equation (2.4) by the velocity fluctuation $u_{i}^{\prime }$ ,

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{i}^{\prime }/2}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{i}^{\prime }u_{j}^{\prime }/2}{\unicode[STIX]{x2202}x_{j}}}+{\mathcal{S}}x_{2}{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{i}^{\prime }/2}{\unicode[STIX]{x2202}x_{1}}}+{\mathcal{S}}u_{1}^{\prime }u_{2}^{\prime }\right)=-{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }p}{\unicode[STIX]{x2202}x_{i}}}+u_{i}^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}+u_{i}^{\prime }f_{i}. & & \displaystyle\end{eqnarray}$$

We make use of

(3.2) $$\begin{eqnarray}\displaystyle u_{i}^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}={\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}-\unicode[STIX]{x1D70F}_{ij}{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}}}={\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}-\unicode[STIX]{x1D70F}_{ij}\unicode[STIX]{x1D60B}_{ij}, & & \displaystyle\end{eqnarray}$$

to obtain

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{i}^{\prime }/2}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{i}^{\prime }u_{j}^{\prime }/2}{\unicode[STIX]{x2202}x_{j}}}+{\mathcal{S}}x_{2}{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }u_{i}^{\prime }/2}{\unicode[STIX]{x2202}x_{1}}}+{\mathcal{S}}u_{1}^{\prime }u_{2}^{\prime }\right)=-{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }p}{\unicode[STIX]{x2202}x_{i}}}+{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}-\unicode[STIX]{x1D70F}_{ij}\unicode[STIX]{x1D60B}_{ij}+u_{i}^{\prime }f_{i}.\quad \,\, & & \displaystyle\end{eqnarray}$$

Equation (3.3) can then be either volume averaged over both phases to obtain the total kinetic energy equation, or phase averaged over the phase $m$ (e.g. carrier or dispersed phase) to obtain the turbulent kinetic energy evolution equation for one phase only.

The equation for the two-fluid mixture is obtained by applying the volume averaging operator

(3.4) $$\begin{eqnarray}\displaystyle \langle \cdot \rangle ={\displaystyle \frac{1}{{\mathcal{V}}}}\int _{{\mathcal{V}}}\cdot \,\text{d}{\mathcal{V}}, & & \displaystyle\end{eqnarray}$$

leading to

(3.5) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\text{d}{\mathcal{K}}}{\text{d}t}}={\mathcal{P}}-\unicode[STIX]{x1D700}+\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}, & & \displaystyle\end{eqnarray}$$

where the different terms indicate the rate of change of turbulent kinetic energy ${\mathcal{K}}$ , the turbulent production rate ${\mathcal{P}}$ , the dissipation rate $\unicode[STIX]{x1D700}$ and the power of the surface tension $\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D70E}}$ , defined as

(3.6a-d ) $$\begin{eqnarray}\displaystyle {\mathcal{K}}=\langle \unicode[STIX]{x1D70C}u_{i}^{\prime }u_{i}^{\prime }\rangle /2,\quad {\mathcal{P}}=-{\mathcal{S}}\langle \unicode[STIX]{x1D70C}u_{1}^{\prime }u_{2}^{\prime }\rangle ,\quad \unicode[STIX]{x1D700}=\langle \unicode[STIX]{x1D70F}_{ij}\unicode[STIX]{x1D60B}_{ij}\rangle ,\quad \unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}=\langle u_{i}^{\prime }f_{i}\rangle . & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}$ is the rate of work performed by the surface tension force on the surrounding fluid. It represents exchange of turbulent kinetic energy and interfacial surface energy and can be either positive or negative and thus a source or sink of turbulent kinetic energy. In particular, $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}$ is proportional to the rate at which surface area is decreasing, i.e.  $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}\propto -\text{d}A/\text{d}t$ (Dodd & Ferrante 2016), and therefore decreasing (increasing) interfacial area through droplet restoration (deformation) or coalescence (break-up) is associated with $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}$ being a source (sink) of turbulent kinetic energy. Note that all the transport terms in (3.3) vanish due to the homogeneity of the domain. On the other hand, if we apply the phase average operator

(3.7) $$\begin{eqnarray}\displaystyle \langle \cdot \rangle _{m}={\displaystyle \frac{1}{{\mathcal{V}}_{m}}}\int _{{\mathcal{V}}_{m}}\cdot \,\text{d}{\mathcal{V}}, & & \displaystyle\end{eqnarray}$$

we obtain

(3.8) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\text{d}{\mathcal{K}}_{m}}{\text{d}t}}={\mathcal{P}}_{m}-\unicode[STIX]{x1D700}_{m}+{\mathcal{T}}_{m}^{\unicode[STIX]{x1D708}}+{\mathcal{T}}_{m}^{p}, & & \displaystyle\end{eqnarray}$$

where the different terms now indicate the rate of change of turbulent kinetic energy ${\mathcal{K}}_{m}$ , the turbulent production rate ${\mathcal{P}}_{m}$ , the dissipation rate $\unicode[STIX]{x1D700}_{m}$ and the viscous ${\mathcal{T}}_{m}^{\unicode[STIX]{x1D708}}$ and pressure ${\mathcal{T}}_{m}^{p}$ powers of the phase $m$ , defined as

(3.9a-c ) $$\begin{eqnarray}\displaystyle {\mathcal{K}}=\langle \unicode[STIX]{x1D70C}u_{j}^{\prime }u_{j}^{\prime }\rangle _{m}/2,\quad {\mathcal{P}}_{m}=-{\mathcal{S}}\langle \unicode[STIX]{x1D70C}u_{1}^{\prime }u_{2}^{\prime }\rangle _{m},\quad \unicode[STIX]{x1D700}=\langle \unicode[STIX]{x1D70F}_{ij}\unicode[STIX]{x1D60B}_{ij}\rangle _{m}, & & \displaystyle\end{eqnarray}$$

(3.9d,e ) $$\begin{eqnarray}\displaystyle {\mathcal{T}}_{m}^{\unicode[STIX]{x1D708}}=\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}}\right\rangle _{m},\quad {\mathcal{T}}_{m}^{p}=-\left\langle {\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }p}{\unicode[STIX]{x2202}x_{i}}}\right\rangle _{m}. & & \displaystyle\end{eqnarray}$$

In this case, the viscous and pressure transport terms are retained to account for a net flux of turbulent kinetic energy from one phase to the other caused by the coupling between the droplets and the carrier fluid (this physical interpretation can be seen more clearly by applying Gauss’s theorem to rewrite the terms as surface integrals, thus resulting in surface integration over the droplet surface). Note finally that the convective transport terms are zero because the two fluids are immiscible and therefore turbulent eddies cannot transport turbulent kinetic energy across the interface.

Figure 11. (a) Turbulent kinetic energy production ${\mathcal{P}}$ and (b) dissipation $\unicode[STIX]{x1D700}$ rates averaged over both phases as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ , normalized by their value in the single-phase flow ( ${\mathcal{P}}_{0}$ and $\unicode[STIX]{x1D700}_{0}$ ).

First, we focus on the equation for ${\mathcal{K}}$ obtained by averaging over the whole volume and over both phases (3.5). At steady state, the rate of change of ${\mathcal{K}}$ is obviously zero and the remaining terms are the production and dissipation rates and the power of surface tension. Figure 11 shows the production ${\mathcal{P}}$ and dissipation $\unicode[STIX]{x1D700}$ rates, normalized by their single-phase values ${\mathcal{P}}_{0}$ and $\unicode[STIX]{x1D700}_{0}$ , for all the simulations performed in the present study as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ . We observe that both the normalized production and dissipation rates are greater than unity and decrease monotonically as the $We_{\unicode[STIX]{x1D706}}$ increases, indicating that the presence of the droplets leads to turbulence augmentation. As $We_{\unicode[STIX]{x1D706}}$ decreases, the droplets become increasingly rigid, and therefore they exert a blocking effect on the surrounding turbulent flow. This effect abruptly re-orients the turbulent eddies leading to an increase in the magnitude of the Reynolds stress, $\langle u_{1}^{\prime }u_{2}^{\prime }\rangle$ , causing an increase in ${\mathcal{P}}$ , which also leads to an increase in the magnitude of the velocity gradients $\unicode[STIX]{x1D60B}_{ij}$ , associated with an increase in $\unicode[STIX]{x1D700}$ relative to the single-phase flow, as shown in figure 11. Moreover, the two quantities have approximately the same value (the difference is less than $3\,\%$ ), thus indicating that at steady state the production balances the dissipation and that the power of surface tension is on average zero (i.e.  ${\mathcal{P}}\approx \unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}\approx 0$ ). These results are consistent with what was previously observed in figure 4(a) and indirectly confirm the relation $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}=-\unicode[STIX]{x1D70E}/{\mathcal{V}}_{m}\,\text{d}{\mathcal{A}}/\text{d}t$ derived by Dodd & Ferrante (2016). Indeed, this relation implies that at steady state $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70E}}$ is zero since the rate of change of ${\mathcal{A}}$ is null.

Figure 12. (a) Turbulent kinetic energy production ${\mathcal{P}}_{m}$ and (b) dissipation $\unicode[STIX]{x1D700}_{m}$ rates averaged over the two phases separately as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ for cases $2$ , $6$ and $10$ . The left and right columns are used to distinguish the dispersed and carrier phases, respectively.

Next, we focus on the equation obtained by phase averaging in one of the two fluids (3.8). Again, at steady state the time derivative on the left-hand side is zero and the relation states that the production and dissipation are balanced by the two transport terms ${\mathcal{T}}_{m}^{\unicode[STIX]{x1D708}}$ and ${\mathcal{T}}_{m}^{p}$ . Figure 12 shows histograms of the production ${\mathcal{P}}_{m}$ and dissipation $\unicode[STIX]{x1D700}_{m}$ rates in the two phases for three selected Weber numbers $We_{\unicode[STIX]{x1D706}}$ (cases $2$ , $6$ and $10$ ). We observe that the production rate is lower in the dispersed phase than in the carrier phase, while the dissipation rate is higher in the dispersed fluid than in the carrier fluid. These results indicate that the total transport term ${\mathcal{T}}_{m}={\mathcal{T}}_{m}^{\unicode[STIX]{x1D708}}+{\mathcal{T}}_{m}^{p}$ is positive in the dispersed fluid and negative in the carrier, corresponding to a turbulent kinetic energy transfer from the carrier to the dispersed phase. In other words, the presence of the droplets is overall a sink for the turbulent kinetic energy of the bulk fluid ${\mathcal{K}}_{c}$ . In addition, we observe that the difference in ${\mathcal{P}}_{m}$ and $\unicode[STIX]{x1D700}_{m}$ decreases with $We_{\unicode[STIX]{x1D706}}$ .

Figure 13. (a) Dispersed and (b) carrier transport terms ${\mathcal{T}}_{m}$ , averaged over the two phases separately as a function of the Weber number $We_{\unicode[STIX]{x1D706}}$ for cases $2$ , $6$ and $10$ . The left and right columns are used to distinguish the pressure and viscous contributions, respectively.

Finally, figure 13 shows the decomposition of the total transport term ${\mathcal{T}}_{m}$ into its pressure and viscous contributions, ${\mathcal{T}}_{m}^{p}$ and ${\mathcal{T}}_{m}^{\unicode[STIX]{x1D708}}$ . In the dispersed phase shown in (a), the pressure transport term is very small and almost negligible, with most of the transport of turbulent kinetic energy ( $90{-}95\,\%$ ) due to the viscous contribution ${\mathcal{T}}_{d}^{\unicode[STIX]{x1D708}}$ . On the other hand, an opposite behaviour is evident in the carrier phase shown in (b): the pressure transport term ${\mathcal{T}}_{c}^{p}$ is dominant one and accounting for most of the transport of turbulent kinetic energy ( $65{-}80\,\%$ ), while the pressure contribution is small. Moreover, we can observe that all the transport terms reduce for increasing Weber number, consistently with the discussion concerning figure 12.

Figure 14. Probability density function of the flow topology parameter ${\mathcal{Q}}$ for three different Weber numbers: cases $2$ (blue line), $6$ (green line) and $10$ (brown line), same as figure 13. The solid and dashed lines are used for the dispersed and carrier phase, respectively.

The different mechanism of transport of turbulent kinetic energy between the carrier and dispersed phase is due to the different kind of flow experienced by the two fluids. This is discussed in figure 14 where the so-called flow topology parameter ${\mathcal{Q}}$ (see e.g. De Vita et al. 2018) is presented. The flow topology parameter is defined as

(3.11) $$\begin{eqnarray}\displaystyle {\mathcal{Q}}={\displaystyle \frac{\unicode[STIX]{x1D60B}^{2}-\unicode[STIX]{x1D6FA}^{2}}{\unicode[STIX]{x1D60B}^{2}+\unicode[STIX]{x1D6FA}^{2}}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D60B}^{2}=\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ji}$ and $\unicode[STIX]{x1D6FA}^{2}=\unicode[STIX]{x1D6FA}_{ij}\unicode[STIX]{x1D6FA}_{ji}$ , with $\unicode[STIX]{x1D6FA}_{ij}$ the rate of rotation tensor, $\unicode[STIX]{x1D6FA}_{ij}=(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}-\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})/2$ . When ${\mathcal{Q}}=-1$ the flow is purely rotational, regions with ${\mathcal{Q}}=0$ represent pure shear flow and those with ${\mathcal{Q}}=1$ elongational flow. The distribution of the flow topology parameter for three selected cases is reported in figure 14. Note that in the figure we show the probability density function (p.d.f.) of ${\mathcal{Q}}$ in the two liquid phases separately. We observe that in the carrier fluid (dashed lines) the flow is mostly a shear flow as demonstrated by a single broad peak at ${\mathcal{Q}}=0$ , and that little changes when changing the Weber number. On the other hand, the flow of the dispersed fluid (solid lines) still shows a broad single peak, now shifted towards negative values of ${\mathcal{Q}}$ , meaning that the flow is more rotational. Also, the relevance of the rotational flow is more and more evident as the Weber number increases. This is caused by the increased number of droplets and their consequent reduction in size: indeed, as the droplets size reduces the effect of the shear reduces as well.

4 Conclusions

We perform direct numerical simulations of two-phase homogeneous shear turbulent flows at $Re_{z}=15\,200$ , where the two-phase nature of the problem is tackled numerically using the MTHINC volume of fluid method recently developed. The droplets are initially spheres providing $5\,\%$ volume fraction of the suspended phase and various Weber numbers and droplet initial diameters are investigated.

We show that the two-phase flow is able to reach a statistically steady state as indicated by a balance of turbulent kinetic energy production and dissipation. The results show that the presence of the droplets leads to turbulence augmentation by increasing the dissipation and production rates of the turbulence relative to the droplet-free flow. In particular, we find that as the Weber number decreases (higher droplet surface tension), the dissipation rate increases, causing the Taylor-microscale Reynolds number to decrease. This is explained by the surface tension force exerting a blocking effect on the surrounding turbulent flow. The turbulent production and dissipation rates are on average equal and in balance, with values larger than their single-phase counterparts. Also, the surface tension power is on average zero. The flow modifications are caused by the presence of the dispersed phase, which acts as a sink of turbulent kinetic energy for the carrier phase, with a net flux going from the bulk of the fluid to the dispersed phase where it is dissipated. Moreover, the transport of turbulent kinetic energy in the carrier fluid is mainly due to the pressure transport, while the one inside the dispersed phase is dominated by the viscous contribution. This difference is explained by the different nature of the flow in the two phases: the carrier fluid is mainly a shear flow, while the dispersed fluid is more rotational owing to its smaller length scales where the effect of the mean shear is reduced.

In addition to the flow properties, the droplet distribution eventually reaches a statistically stationary condition. Indeed, we show that the flow reaches a condition where the number of droplets remains almost constant, due to a balance between the break-up and coalescence mechanisms, and that the number of droplets grows approximately linearly with the Weber number. A similar trend is found for the averaged surface area which also grows monotonically with the Weber number, but the growth rate is less than linear (the surface area grows with the Weber number to the power of $1/3$ , at least for moderately large Weber numbers). With the exception of one case, the droplet size distribution is single peaked, with the mean droplet size reducing with the Weber number. Based on the size distribution data, we show that the maximum droplets size scales well with the energy input as proposed by Hinze (1955), although the possibility of a coalescence mechanism and the presence of a mean shear were not considered in the original formulation by Hinze (1955).