Droplets in homogeneous shear turbulence

TL;DR

This study simulates turbulent shear flow with dispersed droplets, revealing how the dispersed phase influences turbulence characteristics, energy spectra, and droplet size distribution, with findings applicable to multiphase flow modeling.

Contribution

It provides new insights into the behavior of droplets in shear turbulence, including energy transfer, size distribution, and validation of the Hinze estimate under shear conditions.

Findings

Dispersed phase reduces Reynolds number compared to single-phase flow.

Energy at small scales increases due to droplets, diminishing with higher Weber number.

Hinze estimate accurately predicts maximum droplet size despite coalescence and shear.

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 15200. 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 droplets 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 breakup in homogeneous isotropic turbulence, provides an excellent estimate notwithstanding the action of significant coalescence and the presence of a mean shear.

Tables

Table 1. Table 1: Summary of the direct numerical simulations performed with different initial droplet sizes D 0 subscript 𝐷 0 D_{0} , numbers of droplets 𝒩 0 subscript 𝒩 0 \mathcal{N}_{0} and surface tension σ 𝜎 \sigma , all at a fixed Reynolds number R ​ e z = 15200 𝑅 subscript 𝑒 𝑧 15200 Re_{z}=15200 and volume fraction Φ = 5 % Φ percent 5 \Phi=5\% .

Case	Symbol	$D_0/L_z$	${𝒩}_0$	$W{e}_{{𝒮}_0}$	$W{e}_{{\text{rms}}_0}$	$W{e}_0$	$W{e}_{\lambda }$	$R{e}_{\lambda }$
$1$	-	$-$	$-$	$-$	$-$	$-$	$-$	$145$
$2$	$\bullet$	$0.36$	$4$	$0.2$	$0.04$	$0.5330$	$0.0220$	$83$
$3$	$\bullet$	$0.16$	$51$	$0.2$	$0.2$	$0.8000$	$0.0776$	$101$
$4$	$\bullet$	$0.08$	$564$	$0.2$	$1$	$2.0943$	$0.9339$	$111$
$5$	$■$	$0.36$	$4$	$1$	$0.2$	$2.0944$	$0.6773$	$113$
$6$	$■$	$0.16$	$51$	$1$	$1$	$4.0156$	$0.7536$	$117$
$7$	$■$	$0.08$	$564$	$1$	$5$	$10.4717$	$4.9313$	$132$
$8$	$▲$	$0.36$	$4$	$5$	$1$	$4.1890$	$2.0103$	$122$
$9$	$▲$	$0.16$	$51$	$5$	$5$	$7.9999$	$4.0868$	$131$
$10$	$▲$	$0.08$	$564$	$5$	$25$	$20.9432$	$13.3057$	$142$

Equations

$$u_i\left(x_1+L_x,x_2,x_3\right)$$

$$u_i\left(x_1,x_2+L_y,x_3\right)$$

$$u_i\left(x_1,x_2,x_3+L_z\right)$$

$$\rho \left(\frac{\partial u_i^{\prime }}{\partial t}+\frac{\partial u_i^{\prime }{u}_j^{\prime }}{\partial x_j}+\mathcal{S}{x}_2\frac{\partial u_i^{\prime }}{\partial x_1}+\mathcal{S}{u}_2^{\prime }{\delta }_{i1}\right)$$

$$\frac{\partial u_i^{\prime }}{\partial x_i}$$

$$\frac{\partial \varphi }{\partial t}+\frac{\partial u_{i}H}{\partial x_i}=\varphi \frac{\partial u_i}{\partial x_i},$$

$$H\left(X,Y,Z\right)\approx\widehat{H}\left(X,Y,Z\right)=\frac{1}{2}\bigg{(}1+\tanh\big{(}\beta\left(P\left(X,Y,Z\right)+d\right)\big{)}\bigg{)},$$

$${u_i^{\prime }}^{\ast }={u_i^{\prime }}^n+\Delta t\left(\frac{3}{2}{\text{rhs}}_i^n-\frac{1}{2}{\text{rhs}}_i^{n-1}\right),$$

$${\text{rhs}}_i=-\mathcal{S}{u}_2^{\prime }{\delta }_{i1}-\frac{\partial u_i^{\prime }{u}_j^{\prime }}{\partial x_j}+\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial x_j}.$$

$${u_i^{\prime }}^{\ast \ast }\left(x_1,x_2,x_3\right)={u_i^{\prime }}^{\ast }\left(x_1-\Delta t\mathcal{S}{x}_2,x_2,x_3\right).$$

$${u_i^{\prime }}^{\ast \ast \ast }={u_i^{\prime }}^{\ast \ast }+\Delta t\frac{f_i^{n+1}}{\rho }.$$

$$\frac{{\partial }^2{p}^{n+1}}{\partial x_j\partial x_j}=\frac{\rho }{\Delta t}\frac{\partial {u_i^{\prime }}^{\ast \ast \ast }}{\partial x_i},$$

$${u_i^{\prime }}^{n+1}={u_i^{\prime }}^{\ast \ast \ast }-\Delta t\frac{1}{\rho }\frac{\partial p^{n+1}}{\partial x_i}.$$

$$R{e}_z=\frac{\mathcal{S}{L}_z^2}{\nu },$$

$$W{e}_{{\mathcal{S}}_0}=\frac{\rho {\mathcal{S}}^2{D}_0^3}{\sigma },$$

$$R{e}_{\lambda }={\left(\frac{2\mathcal{K}}{3}\right)}^{1/2}\frac{\lambda }{\nu }={\left(\frac{5}{3\nu \epsilon }\right)}^{1/2}2\mathcal{K},$$

$$S^{\ast }=\frac{2\mathcal{S}\mathcal{K}}{\epsilon },$$

$$W{e}_{{\text{rms}}_0}=\frac{2\rho \mathcal{K}{D}_0}{\sigma }.$$

$$W{e}_0=\frac{\rho {\left(\sqrt{2\mathcal{K}}+\mathcal{S}{D}_0\right)}^2{D}_0}{\sigma }.$$

$$W{e}_{\lambda }=\frac{\rho {\left(\sqrt{2\mathcal{K}}+\mathcal{S}\lambda \right)}^2\lambda }{\sigma }.$$

$$\rho \left(\frac{\partial u_i^{\prime }{u}_i^{\prime }/2}{\partial t}+\frac{\partial u_i^{\prime }{u}_i^{\prime }{u}_j^{\prime }/2}{\partial x_j}+\mathcal{S}{x}_2\frac{\partial u_i^{\prime }{u}_i^{\prime }/2}{\partial x_1}+\mathcal{S}{u}_1^{\prime }{u}_2^{\prime }\right)=-\frac{\partial u_i^{\prime }p}{\partial x_i}+u_i^{\prime }\frac{\partial {\tau }_{ij}}{\partial x_j}+u_i^{\prime }{f}_i.$$

$$u_i^{\prime }\frac{\partial {\tau }_{ij}}{\partial x_j}=\frac{\partial u_i^{\prime }{\tau }_{ij}}{\partial x_j}-{\tau }_{ij}\frac{\partial u_i^{\prime }}{\partial x_j}=\frac{\partial u_i^{\prime }{\tau }_{ij}}{\partial x_j}-{\tau }_{ij}{\mathcal{D}}_{ij},$$

$$\rho \left(\frac{\partial u_i^{\prime }{u}_i^{\prime }/2}{\partial t}+\frac{\partial u_i^{\prime }{u}_i^{\prime }{u}_j^{\prime }/2}{\partial x_j}+\mathcal{S}{x}_2\frac{\partial u_i^{\prime }{u}_i^{\prime }/2}{\partial x_1}+\mathcal{S}{u}_1^{\prime }{u}_2^{\prime }\right)=-\frac{\partial u_i^{\prime }p}{\partial x_i}+\frac{\partial u_i^{\prime }{\tau }_{ij}}{\partial x_j}-{\tau }_{ij}{\mathcal{D}}_{ij}+u_i^{\prime }{f}_i.$$

$$⟨\cdot ⟩=\frac{1}{\mathcal{V}}{\int }_{\mathcal{V}}\cdot \mathrm{d}\mathcal{V},$$

$$\frac{d\mathcal{K}}{dt}=\mathcal{P}-\epsilon +{\Psi }_{\sigma },$$

$$\mathcal{K}=⟨\rho u_i^{\prime }{u}_i^{\prime }⟩/2,\mathcal{P}=-\mathcal{S}⟨\rho u_1^{\prime }{u}_2^{\prime }⟩,\epsilon =⟨{\tau }_{ij}{\mathcal{D}}_{ij}⟩,{\Psi }_{\sigma }=⟨u_i^{\prime }{f}_i⟩.$$

$$⟨\cdot {⟩}_m=\frac{1}{{\mathcal{V}}_m}{\int }_{{\mathcal{V}}_m}\cdot \mathrm{d}\mathcal{V},$$

$$\frac{d{\mathcal{K}}_m}{dt}={\mathcal{P}}_m-{\epsilon }_m+{\mathcal{T}}_m^{\nu }+{\mathcal{T}}_m^p,$$

$$\mathcal{K}=⟨\rho u_j^{\prime }{u}_j^{\prime }{⟩}_m/2,{\mathcal{P}}_m=-\mathcal{S}⟨\rho u_1^{\prime }{u}_2^{\prime }{⟩}_m,\epsilon =⟨{\tau }_{ij}{\mathcal{D}}_{ij}{⟩}_m,{\mathcal{T}}_m^{\nu }={⟨\frac{\partial u_i^{\prime }{\tau }_{ij}}{\partial x_j}⟩}_m,{\mathcal{T}}_m^p=-{⟨\frac{\partial u_i^{\prime }p}{\partial x_i}⟩}_m.$$

$$\mathcal{Q}=\frac{{\mathcal{D}}^2-{\Omega }^2}{{\mathcal{D}}^2+{\Omega }^2},$$

Full text

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 on 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 et al., 1998; Lovick et al., 2005) and computational studies (Perlekar et al., 2012; Skartlien et al., 2013; Komrakova et al., 2015; Scarbolo et al., 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, breakup, or coalescence of the droplets introduces an additional term to the turbulent kinetic energy equation - the power of the surface tension - termed $\Psi_{\sigma}$ 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 et al., 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 et al., 2013). This condition limits any numerical simulations to relatively short times, concerning only the initial shearing of isotropic turbulence (Rogers & Moin, 1987; Lee et al., 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 DNS 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. (a)

Can a statistically stationary state be reached when the dispersed phase actively undergoes breakup and coalescence in homogeneous shear turbulence? 2. (b)

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

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 et al., 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 and also cost effective compared to the well known geometric volume of fluid methods in canonical test cases (Xie et al., 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 et al., 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 section 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 section 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 section 4.

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.,

[TABLE]

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}\delta_{1i}$ where $\delta_{ij}$ is the Kronecker delta, and a fluctuating part $u^{\prime}_{i}$ ($u^{\prime}_{i}=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 (′), and the average in the full volume with $\langle\cdot\rangle$. The time evolution of the fluctuating velocity $u^{\prime}_{i}$ is described by

[TABLE]

where $\rho$ is the fluid density, $p$ is the pressure, $\tau_{ij}=2\mu\mathcal{D}_{ij}$ with $\mu$ the dynamic viscosity and $\mathcal{D}_{ij}$ the strain rate tensor ($\mathcal{D}_{ij}=\left(\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i}\right)/2$), and $f_{i}$ is the surface tension force defined as $f_{i}=\sigma\kappa n_{i}\delta$, where $\delta$ is the Dirac delta function at the interface, $\sigma$ the interfacial surface tension, $\kappa$ the interface curvature and $n_{i}$ the normal to the interface. This equation is written in the so-called one-fluid formulation (Tryggvason et al., 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,

[TABLE]

where $\phi$ is the cell-averaged value of the indicator function.

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

[TABLE]

where $X,Y,Z\in\left[0,1\right]$ is a centered 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 $\phi$ and $\beta$ is a sharpness parameter. $\beta$ 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 et al. (1989) and recently adopted by Tanaka (2017), in which the third term on the left-hand side of the momentum equation (equation (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 equation (4) to obtain an intermediate velocity

[TABLE]

where $\Delta t$ is the time step from time $t^{n}$ to $t^{n+1}$ and

[TABLE]

The time step $\Delta t$ is chosen such that the Courant–Friedrichs–Lewy (CFL) number $U_{max}\Delta t/\Delta 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

[TABLE]

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 et al. (1992) to compute the surface tension force where the normals are obtained with the well known Youngs approach (Youngs, 1982), i.e., $f_{i}=\sigma\kappa\partial\phi/\partial x_{i}$, thus obtaining

[TABLE]

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

[TABLE]

which is solved with a standard FFT-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

[TABLE]

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

2.2 Setup

The problem is governed by several dimensionless parameters, which define the problem under consideration. First, the computational box is defined by two aspect ratios $\AR_{xz}=L_{x}/L_{z}$ and $\AR_{yz}=L_{y}/L_{z}$ which are fixed equal to $2.05$ and $1.025$ respectively. These values have been chosen accordingly to what 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

[TABLE]

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

[TABLE]

and the ratio of the initial droplet diameter to the box size $\AR_{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 $15200$; in the multiphase cases, we vary the ratio $\AR_{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 $\sigma$. Two other nondimensional parameters are the density and viscosity ratios, which are fixed to unity to study the individual effect of the Weber number (interfacial surface tension).

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 nondimensional 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_{\lambda}$ and the shear-rate parameter $S^{*}$, defined as

[TABLE]

and

[TABLE]

where $\lambda=\sqrt{10\nu\mathcal{K}/\varepsilon}$ is the Taylor microscale (Sekimoto et al., 2016), $\mathcal{K}=\langle\rho u^{\prime}_{i}u^{\prime}_{i}\rangle/2$ is the turbulent kinetic energy per unit volume, and $\varepsilon=\mu\langle\partial u^{\prime}_{i}/\partial x_{j}\partial u^{\prime}_{i}/\partial 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 $\tau_{0}=\left(2\mathcal{K}\right)^{1/2}/\varepsilon$ and the Kolmogorov time scale $\tau_{K}=\left(\nu/\varepsilon\right)^{1/2}$ and the mean shear time scale $\tau_{\mathcal{S}}=1/\mathcal{S}$, respectively.

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

[TABLE]

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_{\textrm{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 the effects as

[TABLE]

Finally, we can define a Weber number based on $\lambda$ as

[TABLE]

The choice of using $\lambda$ in the definition of the Weber number instead of a dimension associated to 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_{\lambda}$ 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_{\lambda}=50.8$, obtained in a square computational box of size $2\pi$ discretised with $256$ grid points in each direction. From the time history of the turbulent kinetic energy $\mathcal{K}$ and of the enstrophy $\Omega=\langle\omega_{i}\omega_{i}\rangle$, shown in figure 3(a), we observe a first transient phase for $0\leq t\mathcal{S}\leq 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 behavior is well captured in our simulation. A quantitative validation is performed first by comparing the mean components of the velocity anisotropy tensor, $b_{ij}=\langle u_{i}u_{j}/u^{\prime}_{k}u^{\prime}_{k}-\delta_{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.

3 Results

3.1 Statistically stationary state

We start our analysis by considering the single-phase flow at $Re_{z}=15200$. 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 $\Delta x/\eta\approx 0.7$, where $\eta$ is the Kolmogorov scale defined as $\eta=\left(\mu/\rho\right)^{3/4}/\varepsilon^{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\varepsilon$). 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_{\lambda}\approx 145$ with the averaged spectrum of the TKE 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 wave numbers because of the large-scale anisotropy, while all spectra coincide at higher wave numbers, consistently with what observed by Pumir (1996).

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 $\lambda$. 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 regime, the turbulent production balances on average the dissipation rate ($\mathcal{P}\approx\varepsilon$).

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 wave numbers (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 et al., 2010) and bubbles (Dodd & Ferrante, 2016); the increased energy at high wave numbers has been explained by the breakup 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.

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 that, in the following sections we prefer to use the Weber number based on a flow length scale, $We_{\lambda}$, reported in figure 5(a) as a function of $We_{0}$. We can observe that the two Weber numbers are well correlated, with $We_{\lambda}$ scaling approximately as the square of $We_{0}$, i.e., $We_{\lambda}\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_{\lambda}$, shown for all our simulations in figure 5(b) as a function of the Weber number based on the Taylor microscale $We_{\lambda}$, and also reported in table 1. We can observe that the Reynolds number grows with $We_{\lambda}$ 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_{\lambda}$ increases, consistently with what was already observed in figure 4(b). Indeed, the Reynolds number for the case with the most rigid droplets ($We_{\lambda}\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_{\lambda}\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 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 7(a) shows the temporal evolution of the number of droplets ($\mathcal{N}$) under various $We_{\lambda}$ and initial sizes $D_{0}$. The counting of the droplets is conducted automatically by checking the connectivity of the local VOF field ($\phi$) using a $n$-dimensional image processing library111scipy.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 behaviors 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; iii) case $3$ shows an initial decrease of the number of droplets followed by an increase. Notwithstanding the different behaviors, 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_{\lambda}$ is reported in figure 7(b); we observe that $\mathcal{N}_{s}$ grows monotonically with $We_{\lambda}$ (see also the visualisations 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_{\lambda}$ 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_{\lambda}$. 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_{\lambda}$. This provides additional evidence that the droplet statistics are better defined by the Weber number $We_{\lambda}$ 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 breakup and coalescence reflects the history of the flow features, and at equilibrium measurable quantities depend only on the global physical parameters.

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 ($d^{2}\mathcal{V}/dD^{2}=0$), thus indicating that the probability density plot ($d\mathcal{V}/dD$) is single peaked. In figure 8(a) the Weber number $We_{\lambda}$ 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_{\lambda}$, 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., $d\mathcal{V}/dD$ 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(left). Nevertheless, the overall trend of decreasing size for increasing Weber number is still consistent with the linear scaling between $\mathcal{N}_{s}$ and $We_{\lambda}$, as already observed in figure 7(b). The two bottom panels in figure 8, are contours of $\mathcal{V}/\mathcal{V}_{\textrm{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 centered in the region where the gradient of the color 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_{\lambda}$. The two figures show the transient behavior 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_{\textrm{max}}$. Assuming breakup of droplets due to the dynamic pressure ($\sim\rho U^{2}$), Hinze (1955) proposed that the largest possible droplet in a turbulent emulsifier is determined by the velocity fluctuation across $D_{\textrm{max}}$, i.e., one can define a critical Weber number $We_{\textrm{crit}}=\rho\overline{u^{\prime 2}}D_{\textrm{max}}/\sigma$, above which the droplet breaks up. Hinze (1955) showed that simple dimensional analysis leads to $D_{\textrm{max}}\propto\varepsilon^{-2/5}$, if isotropy prevails and the scaling by Kolmogorov (1941) is assumed valid. $D_{\textrm{max}}$ can be in general approximated by the diameter of the equivalent droplet occupying $95\%$ of the total dispersed volume, i.e., $D_{\textrm{max}}\approx D_{95}$, which is represented in figure 8(top) 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 normalised $D_{95}$ as a function of the scaled energy input, and indeed we can observe that our data scales with an approximately $-2/5$ slope. We remark that, although Hinze developed his theory considering only isotropic turbulent flows dominated by the breakup process and neglecting 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 breakup 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_{\lambda}$. 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 analyzed.

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_{\lambda}$ and clearly shows that the surface area increases monotonically with the Weber number. As we have shown above that $\mathcal{N}\propto We_{\lambda}$, combining 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_{\lambda}^{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_{\lambda}$ 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 equation (4) by the velocity fluctuation $u^{\prime}_{i}$,

[TABLE]

We make use of

[TABLE]

to obtain

[TABLE]

Equation (23) 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

[TABLE]

leading to

[TABLE]

where the different terms indicate the rate of change of turbulent kinetic energy $\mathcal{K}$, the turbulent production rate $\mathcal{P}$, the dissipation rate $\varepsilon$ and the power of the surface tension $\psi_{\sigma}$, defined as

[TABLE]

$\Psi_{\sigma}$ 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, $\Psi_{\sigma}$ is proportional to the rate at which surface area is decreasing, i.e., $\Psi_{\sigma}\propto-dA/dt$ (Dodd & Ferrante, 2016), and therefore decreasing (increasing) interfacial area through droplet restoration (deformation) or coalescence (breakup) is associated with $\Psi_{\sigma}$ being a source (sink) of turbulent kinetic energy. Note that all the transport terms in equation (23) vanish due to the homogeneity of the domain. On the other hand, if we apply the phase average operator

[TABLE]

we obtain

[TABLE]

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 $\varepsilon_{m}$ and the viscous $\mathcal{T}^{\nu}_{m}$ and pressure $\mathcal{T}^{p}_{m}$ powers of the phase $m$, defined as

[TABLE]

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 the 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 can not transport turbulent kinetic energy across the interface.

First, we focus on the equation for $\mathcal{K}$ obtained by averaging over the whole volume and over both phases (equation (25)). 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 $\varepsilon$ rates, normalized by their single-phase values $\mathcal{P}_{0}$ and $\varepsilon_{0}$, for all the simulations performed in the present study as a function of the Weber number $We_{\lambda}$. We observe that both the normalized production and dissipation rates are greater than unity and decrease monotonically as the $We_{\lambda}$ increases, indicating that the presence of the droplets leads to turbulence augmentation. As $We_{\lambda}$ 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 $\mathcal{D}_{ij}$, associated with an increase in $\varepsilon$ 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\varepsilon$ and $\Psi_{\sigma}\approx 0$). These results are consistent with what was previously observed in figure 4(a) and indirectly confirm the relation $\Psi_{\sigma}=-\sigma/\mathcal{V}_{m}~{}d\mathcal{A}/dt$ derived by Dodd & Ferrante (2016). Indeed, this relation implies that at steady state $\Psi_{\sigma}$ is zero since the rate of change of $\mathcal{A}$ is null.

Next, we focus on the equation obtained by phase averaging in one of the two fluids (equation (28)). 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}^{\nu}_{m}$ and $\mathcal{T}^{p}_{m}$. Figure 12 shows histograms of the production $\mathcal{P}_{m}$ and dissipation $\varepsilon_{m}$ rates in the two phases for three selected Weber numbers $We_{\lambda}$ (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}^{\nu}_{m}+\mathcal{T}^{p}_{m}$ 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 $\varepsilon_{m}$ decreases with $We_{\lambda}$.

Finally, figure 13 shows the decomposition of the total transport term $\mathcal{T}_{m}$ into its pressure and viscous contributions, $\mathcal{T}^{p}_{m}$ and $\mathcal{T}^{\nu}_{m}$. In the dispersed phase shown in the left panel, 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}^{\nu}_{d}$. On the other hand, an opposite behavior is evident in the carrier phase shown in the right panel: the pressure transport term $\mathcal{T}^{p}_{c}$ 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.

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

[TABLE]

where $\mathcal{D}^{2}=\mathcal{D}_{ij}\mathcal{D}_{ji}$ and $\Omega^{2}=\Omega_{ij}\Omega_{ji}$, being $\Omega_{ij}$ the rate of rotation tensor, $\Omega_{ij}=(\partial u_{i}/\partial x_{j}-\partial u_{j}/\partial 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 (pdf) 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}=15200$, where the two-phase nature of 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 coalescence mechanism and the presence of a mean shear which were not considered in the original formulation by Hinze (1955).

Acknowledgments

M.R. and L.B. acknowledge financial support by the European Research Council Grant no. ERC-2013-CoG-616186, TRITOS and by the Swedish Research Council Grant no. VR 2014-5001. S.S.J. was supported by the Franklin P. and Caroline M. Johnson Fellowship. The authors acknowledge computer time provided by the Swedish National Infrastructure for Computing (SNIC), the National Infrastructure for High Performance Computing and Data Storage in Norway (project no. NN9561K), and the Certainty cluster awarded by the National Science Foundation to CTR. Finally, special thanks are given to the CTR for hosting M.R.  Z.G. and L.B. in the 2018 CTR Summer Program and to the other participants of the Multiphase Flows group for useful discussions.