Modulation of homogeneous and isotropic turbulence in emulsions

Abstract

We present a numerical study of emulsions in homogeneous and isotropic turbulence (HIT) at $Re_\lambda =137$. The problem is addressed via direct numerical simulations, where the volume of fluid is used to represent the complex features of the liquid–liquid interface. We consider a mixture of two iso-density fluids, where fluid properties are varied with the goal of understanding their role in turbulence modulation, in particular the volume fraction ($0.03<\alpha <0.5$), viscosity ratio ($0.01<\mu _d/\mu _c<100$) and large-scale Weber number ($10.6< We_\mathcal {L}<106.5$). The analysis, performed by studying integral quantities and spectral scale-by-scale analysis, reveals that energy is transported consistently from large to small scales by the interface, and no inverse cascade is observed. Furthermore, the total surface is found to be directly proportional to the amount of energy transported, while viscosity and surface tension alter the dynamic that regulates energy transport. We also observe the $-10/3$ and $-3/2$ scaling on droplet size distributions, suggesting that the dimensional arguments that led to their derivation are verified in HIT conditions.

1. Introduction

Emulsions are multiphase flows of two immiscible (totally or partially) liquid phases with similar densities. Such flows are extremely common in industrial applications such as pharmaceuticals (Nielloud 2000; Spernath & Aserin 2006), food processing (McClements 2015), oil production (Kokal 2005; Mandal et al. 2010; Kilpatrick 2012) and waste treatment. Emulsions are also relevant for environmental flows such as oil spilling in oceans, when the oil droplets distribution becomes fundamental for quantifying environmental damage (Li & Garrett 1998; French-McCay 2004; Gopalan & Katz 2010). Many studies have been performed on the rheological behaviour of emulsions in the past (Einstein 1906, 1911; Jansen, Agterof & Mellema 2001; Pal 2000, 2001; De Vita et al. 2019), while the current knowledge on their behaviour in turbulent flows is limited (Yi, Toschi & Sun 2021).

The two fluids are usually referred to as continuous phase (or carrier phase in case of strong advection) and dispersed phase (or droplet phase) depending on whether the volume fraction $\alpha$ is respectively greater or lower than $0.5$; the system is denoted as binary flow when $\alpha = 0.5$. As the density ratio is usually considered to be close to 1, gravity effects are negligible with respect to the stirring and advection needed to sustain turbulence in the flow. For this reason, four dimensionless numbers can be used to describe these flows, namely the Reynolds number $Re$, the Weber number $We$, the volume fraction of the dispersed phase, and the viscosity contrast. Depending on the specific configuration under investigation, the definitions of these numbers can change, yet they completely define the case studied provided that the two fluid have the same density.

Several aspects of fundamental importance in emulsions, such as turbulence modulation, droplet size distributions and interphase energy fluxes, are not fully understood. We therefore aim to partially fill this gap by means of numerical simulations. In the following, we provide an overview of the main results available in the literature. Results for bubble/droplet-laden flows are also discussed when relevant to the present work.

1.1. Observations on droplet size distribution

The droplet size distribution (DSD) is a key aspect of emulsions, as its prediction becomes fundamental in most applications. In his early seminal work, Kolmogorov (1949) discussed the criteria under which a droplet undergoes breakup when subject to surrounding turbulence. Kolmogorov first proposed a dimensional argument according to which surface tension forces need to be balanced locally by turbulent energy fluctuations. This idea was later addressed in Hinze (1955) and translated into a critical Weber number $We_c$ of order 1 at which breakup occurs, leading to the definition of the Hinze scale $d_H$ as the minimum droplet diameter at which breakup may occur due to pressure fluctuations. A general definition for this scale is

(1.1)\begin{equation} d_H = \left(\frac{We_c}{2}\right)^{3/5} \left(\frac{\sigma}{\rho_c}\right)^{3/5}\varepsilon^{{-}2/5}, \end{equation}

where $\sigma$ is the surface tension coefficient, $\rho _c$ is the carrier phase density and $\varepsilon$ is the energy dissipation rate. This estimate proved valid for bubbles (Chan et al. 2021; Masuk, Salibindla & Ni 2021) and emulsions (Perlekar et al. 2012; Mukherjee et al. 2019; Rosti et al. 2020; Yi et al. 2021). Different ${O}(1)$ values have been reported for $We_c$ in numerical (Rivière et al. 2021) and experimental works (Deane & Stokes 2002; Lemenand et al. 2017), from $0.5$ up to $5$; for dilute emulsions in turbulence $We_c\approx 1.17$, according to the values from both numerical (Perlekar et al. 2012) and experimental (Yi et al. 2021) data.

For bubbles larger than the Hinze scale, Garrett, Li & Farmer (2000) found that in isotropic turbulent conditions, droplets break with a cascade process, and the diameter distribution follows a $d^{-10/3}$ power law. This deterministic process can describe accurately bubble size distributions in breaking waves obtained in experiments (Garrett et al. 2000; Deane & Stokes 2002; Qi, Mohammad Masuk & Ni 2020) and numerical simulations (Deike, Melville & Popinet 2016; Chan et al. 2021). The same power law has also been proposed for emulsions, based on diffuse-interface numerical simulations (Skartlien, Sollum & Schumann 2013; Mukherjee et al. 2019; Soligo, Roccon & Soldati 2019). For bubbles smaller than the Hinze scale, Deane & Stokes (2002) suggested the existence of a fragmentation process; in this case, a $d^{-3/2}$ power law is used to fit experimental data accurately. Agreement with this empirical power law has been observed in homogeneous and isotropic turbulence (HIT) for both bubbles (Rivière et al. 2021) and emulsions (Mukherjee et al. 2019). The transition between the two power laws is defined by the Hinze scale. A consequence of this transition is that droplets with $d\gg d_H$ generate both local and non-local bubble/droplet production, as they can fragment in droplets both larger and smaller than the Hinze scale (Rivière et al. 2021). Although both power laws have been derived under the hypothesis of dilute conditions ($\alpha \lesssim 0.05$), they have been observed recently in HIT studies of dense emulsions (Mukherjee et al. 2019), raising the question on the effective role of coalescence in the process.

The connection between bubbles and emulsions is non-trivial and deserves special attention. Hinze (1955) discussed how $We_c$ depends on the fluid properties of the dispersed phase. He assumed that $We_c=C[1-f(N_{Vi})]$, with $f$ a generic function of the viscosity group $N_{Vi}=\mu _d/\sqrt {\rho _d\sigma d}$, where $\mu _d$ is the dispersed phase viscosity. On the other hand, $d_H$ was derived under the assumption of a dilute emulsion, hence the density in (1.1) refers to the carrier phase, as the phase where the energy dissipation rate $\varepsilon$ could be measured in experiments. This allows the direct application of the Hinze criteria in flows where density/viscosity ratios are significant as in air–water flows. However, significant uncertainties are discussed in the literature about the properties of the function $f$, and the role of the dispersed phase properties remains mostly unknown (Masuk et al. 2021). Also unknown is the role of turbulence inhomogeneity and anisotropy, which, according to Hinze (1955), may be a further source of nonlinear effects in the determination of $We_c$. In fact, in flows where the energy dissipation rate shows strong spatial variations, $We_c$ varies for each bubble/droplet, and it assumes meaning only in an average sense, making it difficult to disentangle the effects of turbulence anisotropy and property contrast. Despite all these uncertainties, correlations from Hinze (1955), Garrett et al. (2000) and Deane & Stokes (2002), derived for isotropic turbulent conditions, apply in most studies with strong property contrasts and large-scale anisotropy. This is likely due to the underlying assumption that the breakup process is purely inertial, as it depends only on $\varepsilon$ (Garrett et al. 2000). Thus bubble breakup studies become relevant also for the present study.

Finally, it is worth noticing that the flow configuration appears to have a significant impact on DSD, and experimental observations in shear flows can depart quite substantially from the discussed power-law behaviours. The recent work of Yi et al. (2021) presents strong experimental evidences of gamma/log-normal DSD in Taylor–Couette flow, confirming the previous findings of Pacek, Man & Nienow (1998). These configurations are characterized by strong anisotropy, making the comparison with data obtained for emulsions and bubbles in HIT difficult. On the other hand, Soligo et al. (2019) studied breakup and coalescence of emulsions dynamic in a turbulent channel flow. These authors observed the appearance of the $-10/3$ power law for the DSD in the presence of surfactants. It is interesting to observe that in this numerical study, the scaling from Garrett et al. (2000) seems to apply in anisotropic configurations. Fortunately, there has been a significant effort in recreating local HIT conditions in experiments in recent years (Debue et al. 2018; Dubrulle 2019; Knutsen et al. 2020), and new studies are expected to provide new insights on these aspects.

1.2. Studies of two-fluid turbulence

With the advent of more powerful computational resources, a significant number of studies have considered droplets in turbulent flows, yet almost all only through diffuse-interface methods that may display significant mass loss. In their study of emulsions in HIT turbulence, Perlekar et al. (2012) show that a statistical stationary state can be reached for the DSD. In the study, the authors used the pseudo-potential lattice Boltzmann method (Biferale et al. 2011), which compensates mass losses (due to droplets dissolution) by artificially re-inflating existing ones. Simulations of the Cahn–Hilliard–Navier–Stokes formulation are presented in Perlekar et al. (2014) for binary fluids. These authors found that enforcing large-scale HIT arrests coarsening. This result is particularly significant for emulsions (of which binary fluids represent a special case) as it shows that turbulence is the main factor to determine the droplet size. Furthermore, these authors report modified energy spectra for the mixtures, with a crossover in correspondence to the Hinze scale.

Komrakova, Eskin & Derksen (2015) used a free-energy lattice Boltzmann method to simulate numerically emulsion breakup in HIT, induced by an external large-scale linear forcing. Their findings show that energy spectra present deviations with respect to the single-phase configuration and that the numerical method employed may alter the small-scale dynamics of the flow. Finally, increased coalescence is found for volume fractions $\alpha >0.05$ also owing to the nature of the diffuse-interface method.

Droplet interactions with turbulence have been studied by Dodd & Ferrante (2016) in decaying isotropic turbulence. Amongst several observations, these authors discuss the effects of droplet breakup and coalescence on the turbulent kinetic energy budget. Droplet coalescence lowers the total amount of area, hence decreases the surface energy and consequently increases the kinetic energy locally, while the opposite occurs in the case of breakup. More recently, Mukherjee et al. (2019) have studied emulsions in HIT conditions using a pseudo-potential lattice Boltzmann method, discussing droplet statistics and their correlation with the surrounding turbulence. They confirm the findings of Perlekar et al. (2014) for energy spectra pivoting at the Hinze scale, demonstrating that energy is subtracted from large scales and injected at small scales, while no direct observation of the underlying mechanism is presented. These authors also show that the droplet generation can be described through the Weber number spectra. In the same work, Mukherjee and co-workers discuss and demonstrate the need for using a forcing scale smaller than the turbulent-box size in order to achieve a polydisperse droplet distribution. It is important to note that Mukherjee et al. (2019) used a pseudo-potential lattice Boltzmann method, which leads to a significant loss of the dispersed mass during the simulation, as discussed fairly by the authors.

As concerns binary fluids, Perlekar (2019) shows how the presence of interfaces leads to a different energy transfer mechanism, confirming the conclusions in Dodd & Ferrante (2016). Perlekar (2019) uses the scale-by-scale (SBS) energy balance to show that the energy absorption at larger scales is given mainly by the interface source term in the Cahn–Hilliard equation used by the author to describe the multiphase nature of the flow. Furthermore, Perlekar (2019) shows that small-scale statistics are almost unchanged when changing $We$, while they are affected by the Reynolds number. This study complements the previous findings in binary fluids (Perlekar et al. 2014; Perlekar, Pal & Pandit 2017), where coarsening was analysed in three- and two-dimensional turbulence by means of a spinoidal decomposition. Rosti et al. (2020) study droplets in homogeneous shear turbulence, focusing on the effect of the droplet initial diameter and the shear rate magnitude; the results show that a statistically stationary regime (i.e. balance of coalescence and breakup events, and energy balance convergence) can be reached, while the Taylor-scale Reynolds number $Re_\lambda$ decreases with increasing surface tension.

Despite the growing literature on the subject, many issues are yet to be fully understood. In particular, most of the studies have been carried out using diffuse-interface approaches, which cannot exactly represent the surface terms effects yet are key in many situations. In this sense, our work complements the very recent one by Rivière et al. (2021), focused on the bubble breakup dynamics.

1.3. Objectives of the present study

In the present work, we use direct numerical simulations (DNS) to study the effects of viscosity ratio, volume fraction and surface tension on the emulsion turbulent behaviour. The chosen set-up is tri-periodic HIT, with turbulence sustained throughout the simulation time. The analysis is performed at $Re_\lambda \approx 137$, which is larger than that attained in previous interface-resolved numerical studies of multiphase turbulent flows. Also, volume fraction, viscosity ratio and surface tension are varied to cover most relevant applications (Jansen et al. 2001). We analyse the turbulence through global and phase-averaged energy balance, energy spectra, SBS energy budget, and probability density functions (p.d.f.s) for the intermittency analysis. Furthermore, we discuss DSD for all cases. In summary, we will show that: (i) the energy balance is significantly altered by the properties of the dispersed phase; (ii) surface tension forces induce an additional mechanism for energy transfer from larger scale towards the energy dissipation range; (iii) the modified energy transport mechanism alters the energy spectra; (iv) the presence of the interface increases intermittency and alters the small-scale statistics; (v) the DSD displays both the $d^{-3/2}$ and $d^{-10/3}$ power laws, with remarkable accuracy also for $d< d_H$.

2. Methodology

2.1. Governing equations and numerical method

We consider an incompressible flow obeying the continuity and Navier–Stokes equations

(2.1a)\begin{gather} \partial_i u_i = 0, \end{gather}

(2.1b)\begin{gather}\rho(\partial_t u_i+u_j\,\partial_j u_i) ={-}\partial_i p + \partial_i \left[ \mu (\partial_i u_j + \partial_j u_i) \right] +f^\sigma_i +f^T_i , \end{gather}

where $u_i$ is the velocity in the $i$th direction, $p$ is the pressure, and $\rho$ and $\mu$ are the local density and viscosity. The forcing term $f^\sigma _i = \sigma \xi \delta _s n_i$ represents the surface tension force, where $\sigma$ is the surface tension, $\xi$ is the local interface curvature, $n_i$ is the $i$th component of the surface normal vector, and $\delta _s$ the Dirac delta function that ensures that the surface force is applied only at the interface (Tryggvason, Scardovelli & Zaleski 2011). The last term in (2.1b) is the forcing needed to sustain turbulence by injecting energy at the large scales; among the several algorithms available to force sustained homogeneous and isotropic turbulence (e.g. Eswaran & Pope 1988; Rosales & Meneveau 2005; Mallouppas, George & van Wachem 2013; Bassenne et al. 2016), we use here the Arnold–Beltrami–Childress (ABC) forcing (Mininni, Alexakis & Pouquet 2006)

(2.2)\begin{equation} \left.\begin{gathered} f^T_x = A \sin \kappa_0 z + C\cos \kappa_0 y, \\ f^T_y = B\sin \kappa_0 x + A\cos \kappa_0 z, \\ f^T_z = C\sin \kappa_0 y + B\cos \kappa_0 x, \end{gathered}\right\} \end{equation}

with $x,y,z\in [0,2{\rm \pi} ]$. As reported by Podvigina & Pouquet (1994), the ABC forcing creates an unstable single-phase flow for $1/\nu >20$, with $\nu$ the kinematic viscosity and $\kappa _0$ the forcing wavelength.

The description of the code and the algorithm used can be found in Rosti, De Vita & Brandt (2019) and Rosti et al. (2020), together with several validations. The method is therefore described only briefly here; see also Costa (2018) for references to the code structure. The equations are discretized on a staggered uniform Cartesian mesh; the spatial derivatives are computed using second-order centred finite differences, and a second-order Adam–Bashford scheme is used for the time integration. The pressure splitting method presented in Dodd & Ferrante (2014) is used to obtain a constant-coefficient Poisson equation, which we then solve with the direct fast Fourier transform based pressure solver presented in Costa (2018).

The interface between the two fluids is described with the volume of fluid (VOF) method, in particular the Multi-dimensional Tangent Hyperbola Interface Capturing (MTHINC) algorithm developed by Ii et al. (2012). The advection equation for the VOF can be written in divergence form as

(2.3)\begin{equation} \partial_t \phi + \partial_i u_i\,\mathcal{H} = \phi\,\partial_i u_i, \end{equation}

where $\mathcal {H}$ is the colour function assuming the values 0 and 1 in each of the fluids, and $\phi$ is the cell-averaged value of $\mathcal {H}$. In the MTHINC method, the function $\mathcal {H}$ is approximated locally using the hyperbolic tangent:

(2.4)\begin{equation} \mathcal{H}(x',y',z') \approx \tfrac{1}{2}(1+\tanh(\beta(P(x',y',z')+d))), \end{equation}

where $(x',y',z')\in [0,1]$ is the cell-centred local coordinate system, $\beta$ is a sharpness parameter (equal to 1 in the current work), $d$ is a normalization factor, and $P$ is the three-dimensional surface function, assumed here to be quadratic (Ii et al. 2012). The advantage of the method is that (2.4) allows us to solve the fluxes in (2.3) by semi-analytical integration. Once the VOF function $\phi$ is known, we evaluate the local fluid properties as

(2.5)\begin{equation} \left.\begin{gathered} \rho = \rho_d \phi + \rho_c(1-\phi), \\ \mu = \mu_d \phi + \mu_c(1-\phi), \end{gathered}\right\} \end{equation}

where the subscripts $c$ and $d$ indicate carrier and dispersed phase, respectively. Finally, the continuum surface force (CSF) model is used to compute the surface tension force (Brackbill, Kothe & Zemach 1992), with the normal evaluated with Young's method and the curvature as in Ii et al. (2012).

2.2. Flow configuration

All the simulations are performed using the same ABC forcing, injecting energy at wavenumber $\kappa _0=2{\rm \pi} /\mathcal {L}=2$, with $A=B=C=1$, corresponding to $Re_\lambda \approx 137$ for the single-phase flow (see Appendix B for the characteristics of the reference single-phase flow). As reported in the literature (Komrakova et al. 2015; Mukherjee et al. 2019), forcing the second wavelength is recommended in order to avoid coalescence induced by large turbulent structures in periodic domains.

In addition to the Reynolds number, the emulsion flows are characterized by four non-dimensional parameters: the volume fraction $\alpha =\mathcal {V}_d/\mathcal {V}$, defined as the ratio between the volume occupied by the dispersed phase $\mathcal {V}_d$ and the total volume $\mathcal {V}=(2{\rm \pi} )^3$; the viscosity ratio $\mu _d/\mu _c$, where the subscripts $d$ and $c$ indicate the dispersed and carrier phases; the Weber number $We_\mathcal {L}=\rho _c\mathcal {L} u_{rms} ^2/\sigma$, where $u_{rms}$ is the space–time average of the root-mean-square velocity of the single-phase case (which can be related to the forcing amplitude $A=B=C$); and the scale of the ABC forcing $\mathcal {L}$. Finally, the density ratio $\rho = \rho _c/\rho _d$ is kept constant, equal to 1 in this study.

Here, we will vary the dispersed phase volume fraction, the viscosity ratio and the Weber number; the parameters pertaining to the different simulations discussed below are presented in table 1. A convergence study motivating the choice of the resolution in the table is reported in Appendix B. Finally, the resolution $N=512$ has been chosen as it allows us to resolve adequately all the different cases. Note, finally, that the table also indicates the integration time $N_\mathcal {T}$ required to reach statistical convergence of the turbulent quantities and DSD in units of large eddy turnover times, $\mathcal {T}=\mathcal {L} u_{rms}$ (Mininni et al. 2006). The simulations are considered at convergence when global energy production balances dissipation (see § 2.3 and (2.7) for their definitions) with an error of less than $4\,\%$, also implying that the average of the time derivative of the interfacial area is negligible (see Appendix B for further details). Interestingly, $N_\mathcal {T}$ varies significantly with the physical configuration. In particular, starting with the reference cases BE1 and BE2, we observe that increasing $\mu _d/\mu _c$ means longer times are needed to reach a statistically stationary state, which we will attribute to a decrease of the breakup rate. A similar behaviour is observed when decreasing $We_\mathcal {L}$, when higher surface tension forces decrease the probability of breakup events. Finally, large structures become unavoidable when increasing the volume fraction $\alpha$ (Komrakova et al. 2015; Mukherjee et al. 2019), which implies longer simulation times.

Table 1. Parameter settings for the simulations considered in this study: number of grid points in each direction $N$, viscosity ratio $\mu _d/\mu _c$, Weber number $We_\mathcal {L}$, with surface tension $\sigma$, volume fraction $\alpha$ and integration time to reach statistical convergence $N_\mathcal {T}$. All simulations are performed with $\mu _c=0.006$ and the same ABC forcing. Each case is denoted by a letter indicating the parameter that is varied: V for viscosity ratio, C for volume fraction, and W for Weber number. SP are the single-phase flows, and BE are configurations that recur in different parametrizations (base emulsions).

Visualizations of the transient phase to reach the final steady state are reported in figure 1 for the reference case BE1 with $\alpha =0.03$. The simulation starts at $t_0$ using the fully developed single-phase HIT field from case SP2. The dispersed phase is initialized as an ensemble of spheres, which soon deform in the flow as shown in figure 1(b), pertaining to time $t_1=\mathcal {T}/4$. At statistical convergence, $t \approx 10\mathcal {T}$, when statistics are collected, we observe a polydispersed distribution of asymmetric droplets. Note finally that for $\alpha \leqslant 10\,\%$, the simulations are initialized using spherical droplets of size $d_0\approx 0.12L$, while a single spherical droplet of initial size $d_0=(6\alpha L^3/{\rm \pi} )^{1/3}$ was used for larger values of $\alpha$. We have checked that the initial distribution has no effect on the final DSD, as also reported in Mukherjee et al. (2019) for a similar configuration.

Figure 1. Initial evolution of the emulsion flow (example reported for case BE1). The droplets are initialized at $t_0$ in a developed turbulent field. As turbulence is maintained, breakup and coalescence start occurring ($t_1$), and statistical convergence in the DSD is achieved after a few turnover times ($t_2$): (a) $t_0$, (b) $t_1=\mathcal {T}/4$, (c) $t_{2}=10\mathcal {T}$.

2.3. Observables, phase-averaged energy balance and scale-by-scale budget

In this subsection, we introduce the theoretical tools and the physical observable that will be discussed throughout the study. The full mathematical derivation, as well as details on the computation of the energy spectrum and DSD, can be found in Appendix A. The objective of this study is to understand the turbulence modulations induced by a second phase, focusing on comparisons of the energy spectra and the SBS analysis. In particular, we will consider the Taylor-scale Reynolds number, the energy spectra, the phase-averaged energy balance, p.d.f.s of velocity fluctuations, dissipation and vorticity, and the SBS budgets. The Taylor-scale Reynolds number is defined as

(2.6)\begin{equation} Re_\lambda = \left(\frac{u_i u_i}{3}\right)^{1/2} \frac{\lambda}{\nu}, \end{equation}

where $\lambda =(5\nu u_i u_i/\varepsilon )^{1/2}$ is the Taylor scale, with the energy dissipation rate computed as $\varepsilon =\nu \,\partial _i u_j\,\partial _j u_i$ for the reference single-phase flow $Re_\lambda =137$. Here, we compute $\varepsilon$ and all the relevant observables at each computational grid point and then average in space and time. This procedure is required due to material properties discontinuities when $\mu _d/\mu _c$ is varied. Note that from now on, $\varepsilon$ will denote the space–time averaged value.

Further insight on the global behaviour of multiphase flows can be gained through the phase-averaged energy balance:

(2.7a)\begin{gather} \rho\,\partial_t k_m = \mathcal{P}_m - \varepsilon_m + \mathcal{T}^\nu_m + \mathcal{T}^p_m, \end{gather}

(2.7b)\begin{gather}k_m=\langle u_i u_i/2 \rangle_m, \quad \mathcal{P}_m = \rho\langle u_i f^T_i\rangle_m, \quad \varepsilon_m = \langle \nu\,\partial_j u_i\,\partial_i u_j\rangle_m, \end{gather}

(2.7c)\begin{gather}\mathcal{T}^\nu_m = \langle \partial_j \mu u_i \left(\partial_j u_i + \partial_i u_j\right)\rangle_m, \quad \mathcal{T}^p_m ={-}\langle\partial_i u_i p\rangle_m. \end{gather}

Here, $\mathcal {P}_m$ and $\varepsilon _m$ indicate, respectively, production rate and viscous dissipation rate per unit volume in each phase. The terms $\mathcal {T}^\nu _m$ and $\mathcal {T}^p_m$ are the viscous and pressure transport densities, respectively, and represent the coupling between the two phases; when the sum of these two, $\mathcal {T}_m = \mathcal {T}^\nu _m +\mathcal {T}^p_m$, is positive, energy is absorbed from phase $m$, when negative energy is transferred to the other phase. Again, in statistical stationary conditions, $\partial _t k_m\approx 0$. Further details on the derivation of (2.7) can be found in Appendix A.

We now move to spectral space and present the SBS balance. This is derived for the two-fluid flows following the formulation in Olivieri et al. (2020a,b); for more details, the reader is referred to Frisch (1995) and Alexakis & Biferale (2018). Taking the Fourier transform of the momentum equations (2.1b), we obtain

(2.8)\begin{equation} \partial_t \tilde{u_i} + \tilde{G}_i = -\mathrm{i}\kappa\widetilde{p/\rho} - \tilde{V}_i + \tilde{f^\sigma_i} +\tilde{f_i^T}. \end{equation}

Denoting the Fourier transform of a quantity $J(x_i,t)$ as $\tilde {J}(\kappa _i,t)=\mathscr {F}\{J(x_i,t)\}$, with $\kappa _i$ the $i$th component of the wavelength vector, in the expression above we have $\tilde {G}_i = \mathscr {F}\{u_j\,\partial _j u_i\}$ and $\tilde {V}_i = \mathscr {F}\{\partial _i(\nu [\partial _i u_j + \partial _j u_i])\}$. Note that as the viscosity $\mu$ is a function of space and time, we actually compute the dissipation term in physical space to avoid a convolution in the spectral space. We next multiply (2.8) with the complex conjugate of the velocity $\tilde {u}_i^*$, and drop the pressure term by imposing the incompressibility condition $\kappa _i\tilde {u}_i=0$, as in this work $\rho _c=\rho _d=1$. Multiplying the complex conjugate of (2.8) by $\tilde {u}_i$, summing the equations obtained for $\tilde {u}$ and $\tilde {u}^*$, and averaging in time, we finally obtain

(2.9)\begin{equation} \partial_t E(\kappa_i) = T(\kappa_i) + \mathcal{D}(\kappa_i) +\mathcal{S}_\sigma(\kappa_i) + \mathcal{F}(\kappa_i), \end{equation}

where

1. (i) $E = \langle \tilde {u}_i\tilde {u}_i^*\rangle _t$ is the time-averaged kinetic energy in the spectral domain, whose time derivative is zero at statistical steady state;

2. (ii) $T = -\langle \tilde {G}_i\tilde {u}_i^* + \tilde {G}_i^*\tilde {u}_i\rangle _t$ is the time-averaged energy transfer due to the nonlinear term;

3. (iii) $\mathcal {D} = -\langle \tilde {V}_i\tilde {u}_i^* + \tilde {V}_i^*\tilde {u}_i\rangle _t$ is the time-averaged viscous dissipation;

4. (iv) $\mathcal{S}_\sigma = -\langle\tilde{f_i^\sigma}\tilde{u}_i^* + \tilde{f_i^\sigma}^*\tilde{u}_i\rangle_t$ is the time-averaged work of the surface tension force at the different scales;

5. (v) $\mathcal{F} = \langle\tilde{f^T_i}\tilde{u}_i^* + \tilde{f^T_i}^*\tilde{u}_i\rangle_t$ is the time-averaged energy input due to the large-scale forcing.

All of the above are three-dimensional fields in spectral space. Note that at steady state when the total interfacial area is constant, $\mathcal {S}_\sigma$ integrates to zero (Dodd & Ferrante 2016), so this term can be seen effectively as an energy transport due to the surface tension.

Finally, we perform a spherical-shell integral in spectral space and express each term in the budget as a function of the magnitude of the wavevector $\kappa$. This operation results in e.g.

(2.10)\begin{equation} T(\kappa) = \sum_{\kappa<|\kappa_i|<\kappa+1} T(\kappa_i). \end{equation}

This term represents the shell-to-shell energy transfer function for the nonlinear term of the momentum equation, and similarly for the other terms above.

3. Results

3.1. Emulsions at different volume fractions

We first examine the influence of the dispersed-phase volume fraction on the turbulent flow, cases BEx and Cxx in table 1, corresponding to increasing values of $\alpha$ from 3 % to 50 %. A render of the cases discussed here is shown in figure 2, where the iso-contours of VOF fields are shown for volume fractions $0.06$, $0.2$ and $0.5$.

Figure 2. Render of the two-fluid interface (corresponding to the value of the VOF function $\phi =0.5$) for different values of the volume fraction $\alpha$: (a) $\alpha =0.06$, (b) $\alpha =0.2$, (c) $\alpha =0.5$. The vorticity fields are shown on the box faces on a planar view.

The modulation of the turbulence is first quantified in terms of integral quantities. Figure 3 shows $Re_\lambda$, computed according to (2.6), versus the volume fraction $\alpha$. Here, $Re_\lambda$ increases almost linearly with $\alpha$, by approximately 15 % for $\alpha =0.5$. A similar trend is found for $\lambda$, as shown in the inset of figure 3. Considering that the average of $\varepsilon$ and $k$ is approximately constant in all cases (variations of ${\pm }3\,\%$), the increase of $Re_\lambda$ and $\lambda$ is therefore due to the local variations of the ratio $k/\varepsilon$.

Figure 3. Taylor-scale Reynolds number $Re_\lambda$ versus the dispersed-phase volume fraction $\alpha$, for viscosity ratio $\mu =1$ and density ratio $\rho =1$. The inset shows the Taylor scale, $\lambda$, versus the different values of $\alpha$ under investigation.

In particular, the increased values of $k/\varepsilon$ for similar averaged values of the two quantities is attributed to the increased correlation between regions of strong turbulent kinetic energy and low dissipation. Graphical evidence is presented in figure 4, where we show the instantaneous ratio $k/\varepsilon$ for the single-phase (case SP2 in panel a) and multiphase flows (case BE2 with $\alpha =0.1$ in panel b) in logarithmic scale. The figure shows that when the dispersed phase is present, large regions of fluid with higher $k/\varepsilon$ are observed far from the droplet interface (denoted with a white line). This can be explained as follows: as the total dissipation is constant, the local increase of $\varepsilon$ near the interface, as also observed in Dodd & Ferrante (2016), corresponds to a decrease of the dissipation rate in large portions of the fluid, those far from an interface. Considering that the turbulent kinetic energy is less affected by the presence of the interface, the ratio $(k/\varepsilon )$ increases in average. To support this statement, figures 4(c,d) depict the instantaneous energy dissipation rates for the same planes. In the emulsion (panel d), higher values of $\varepsilon$ are found close to the droplet interface and to the clustering regions, while for the single-phase flow (panel c), no specific pattern is observed.

Figure 4. (a,b) Contours of the ratio $k/\varepsilon$ with logarithmic scale in two planes. (c,d) Energy dissipation rate $\varepsilon$. (a,c) Results for the single-phase case SP2. (b,d) Results for case BE2 ($\alpha =0.1$). The white lines represent the VOF iso-contours for $\phi =0.5$.

The one-dimensional energy spectra $E(\kappa )$ multiplied by $\kappa ^{5/3}$, i.e. the so-called compensated spectra, are displayed in figure 5 for the different $\alpha$ considered. The Taylor scale of the single-phase flow is indicated by the dot-dashed line, while the vertical dotted black line is used for the wavenumber $\kappa _H$ corresponding to the Hinze scale, defined as

(3.1)\begin{equation} d_H = 0.725\varepsilon^{{-}2/5}(\rho_c/\sigma)^{{-}3/5}. \end{equation}

Figure 5. Compensated energy spectra for simulations at different volume fractions $\alpha$; the dot-dashed line represents Taylor scale $\lambda$, while the dotted line represents the Hinze scale $d_H$.

Note that the prefactor $0.725$ in (3.1) is set according to the original work of Hinze (1955) for emulsions in HIT conditions, corresponding to $We_c=1.17$.

The data in figure 5 reveal that the presence of the dispersed phase reduces the energy with respect to the single-phase case (SP1) for $\kappa <\kappa _H$. At the same time, the energy content increases at the smaller scales, $\kappa >\kappa _H$, in the dissipative range of the spectra. As noted in previous studies (Mukherjee et al. 2019), the amount of energy subtracted to the large scales is proportional to the volume fraction $\alpha$. Interestingly, the wavenumber at which the curves cross over from reduced to increased energy content corresponds to the Hinze scale. For brevity, we will denote as the pivoting point the wavelength where the spectra of the multiphase cases intersect the one from the single-phase reference case.

Pivoting points were not observed clearly in some previous studies on emulsions (Mukherjee et al. 2019; Perlekar 2019), while they are clearly visible in others (Perlekar et al. 2014; Dodd & Ferrante 2016; Rosti et al. 2020). This is possibly due to the different methods used to simulate the dispersed phase: the ability of the VOF to accurately resolve the interface reduces the energy dissipation by the surface tension term in the dissipative range. Such energy dissipation is indeed observed clearly by Perlekar (2019), who present results obtained by solving the Cahn–Hilliard equation in a diffuse-interface formulation. As mentioned in Appendix B, these numerical artefacts do not have significant effects on the dynamics at the inertial range, while they affect the dissipative range. This aspect will also be discussed later in this section.

Insight into the energy transfer among the different scales is gained by using the SBS analysis. The full SBS energy budget, i.e. the contributions from the different terms in (2.9), is displayed in figure 6(a) for case BE2, chosen as an illustrative example with an intermediate value $\alpha =0.1$. The external forcing is injecting energy at $\kappa =2$, which is absorbed by the nonlinear transfer term $T$ for a large majority, and by the surface tension term $\mathcal {S}_\sigma$ for a small part. The nonlinear term transfers energy towards smaller scales, larger values of $\kappa$. The surface tension term, $\mathcal {S}_\sigma$, acts as a dissipative process at large scales, where it absorbs approximately the same energy as the dissipative term $\mathcal {D}$. However, for $10<\kappa <20$, we observe a significant change in the energy transport mechanism: $\mathcal {S}_\sigma$ becomes positive, hence contributing to transferring energy towards the small scales, similarly to $T$, a process active until $\kappa _{max}$. It is important to note that the surface tension transport remains active also at small scales in the dissipative range, consequently extending the range of wavelengths where the dissipation term remains active. These observations confirm the previous findings obtained in Perlekar (2019) for binary mixtures, while showing that dissipation at small scales is due only to $\mathcal {D}$ in the case of sharp interface methods.

Figure 6. Scale-by-scale energy budget for different volume fractions $\alpha$. (a) Full energy balance for the case BE2 with $\alpha =0.1$; (b) the energy transfer $T$ due to the nonlinear terms; (c) the energy transfer $\mathcal {S}_\sigma$ associated with the surface tension term; (d) energy dissipation rate $\mathcal {D}$.

The details of the effect of the volume fraction $\alpha$ on each term of the SBS balance are displayed in figures 6(b–d). We first analyse the nonlinear transfer term $T$ in panel (b). As $\alpha$ increases, $T$ absorbs progressively less energy at the injection frequency $\kappa =2$. Consequently, less energy is transferred towards smaller scales by nonlinear advection. The energy flux $\varPi$ (not shown) does not display an inverse cascade for any $\alpha$. Furthermore, we notice that no energy is transferred to the far end of the dissipative range, which is resolved over a large range of scales in all cases.

The contribution from the surface tension $\mathcal {S}_\sigma$ (see figure 6c) confirms that interfacial stresses absorb part of the energy injected into the domain at $\kappa =2$. The energy absorbed by the surface tension term at large scales is approximately proportional to $\alpha$. The surface tension term becomes positive at smaller scales, where energy is released. The positive peak is reached at approximately the Hinze scale for all cases. As for the energy absorption, the magnitude of the peak scales proportionally to $\alpha$. We also notice that for any $\alpha$, the surface tension terms act also in the dissipative range, where the nonlinear term $T$ is zero.

The behaviour observed so far for $T$ and $\mathcal {S}_\sigma$ provides a clear explanation for the previous observations on the energy spectra. At small wavenumber, the energy cascade produced by the nonlinear energy transfer is partially inhibited by the presence of the interfacial forces. For high wavenumbers, $T$ reaches zero progressively, but the energy previously subtracted by the interfacial stresses at large scale is redistributed at small scales, which can be seen in figure 5 as an energy increase at high wavenumbers.

To close the SBS balance, we examine the viscous dissipation term $\mathcal {D}$ (see figure 6d). First, we note that only a small amount of the injected energy (less than $5\,\%$ for all cases) is absorbed by the dissipation term at the scale of the forcing, $\kappa = 2$. The overall effect of the dispersed phase is to shift the energy dissipation towards smaller scales. This constitutes the natural reaction of the system to the increased activity in the dissipative range caused by the surface tension term. This behaviour becomes more evident as $\alpha$ increases and progressively enhances dissipation at those small scales where the single-phase dissipation is negligible.

Summarizing, the surface tension introduces an alternative path for energy transmission from large towards small scales, as discussed for binary flows in Perlekar (2019). The amount of energy transferred by the surface tension is directly proportional to the total droplets surface area $\mathcal {A}$, as shown in figure 7(a), where we display the maximum energy transferred via surface tension, $\max (\sum _{|\kappa _i|<\kappa } \mathcal {S}_\sigma (\kappa ))$, and the total area of the dispersed phase $\mathcal {A}$ for the different volume fractions under consideration with a linear fit to the data. This observation reinforces our previous conclusion that the interface transfers energy among different scales by disrupting larger turbulent structures and creating smaller ones, hence affecting the canonical $-5/3$ slope of the turbulence spectra. Note also that while in monodispersed flows this results in a deviation at a specific spectral frequency (Dodd & Ferrante 2016), in polydispersed flows this behaviour is seen at all scales.

Figure 7. (a) Correlation between the maximum surface tension term, $\max (\sum _{|\kappa _i|<\kappa } \mathcal {S}_\sigma (\kappa ))$, and the total surface area $\mathcal {A}$ for the different volume fractions $\alpha$. The dashed black line is the linear fit to the data. (b) P.d.f. of the DSD for different values of $\alpha$. The dashed black line indicates the $d^{-3/2}$ law from Deane & Stokes (2002), the continuous black line the $d^{-10/3}$ law from Garrett et al. (2000), and the dotted black line the Hinze scale $d_H$. The droplet size is normalized by the Kolmogorov scale of simulation SP2, $\eta _{sp}$.

We next consider the dynamics of the dispersed phase. We first examine the DSD for all the values of volume fraction studied – see figure 7(b), where we display the droplet diameters normalized by the single-phase (SP1) Kolmogorov scale $\eta _{sp}$ . The dashed black line depicts the $d^{-3/2}$ law by Deane & Stokes (2002), and the solid line depicts the $d^{-10/3}$ law by Garrett et al. (2000), valid for larger droplets. For small droplets, the $-3/2$ law is well captured also for marginally resolved droplets (with $d/\eta _{sp}<6$). For droplets larger than the Hinze scale, the $-10/3$ law is also a very good fit, with increasing accuracy for increasing values of $\alpha$. Our data are in agreement with the findings by Mukherjee et al. (2019) and explained by higher coalescence probability at higher volume fractions, leading to a bigger population of larger droplets. Interestingly, the Hinze scale turns out to define approximately the transition between the $-3/2$ and $-10/3$ scalings as proposed in Deane & Stokes (2002), although for higher values of $\alpha$, the onset of the $d^{-10/3}$ power law occurs at larger diameters. As the droplet distributions can be, to a good approximation, represented by these two laws, it follows that $\mathcal {A}\propto \alpha$, explaining why $\mathcal {A}$, $\mathcal {S}_\sigma$ and $\alpha$ are linearly correlated (see figure 7a).

We now consider the phase-averaged energy budget, introduced in § 2.3. The different terms of (2.7) – production, dissipation and transport by pressure and viscous forces – are shown in figure 8, normalized by the single-phase dissipation. We first observe that the total production and dissipation is

(3.2)\begin{equation} \mathcal{P}= \alpha \mathcal{P}_d + (1- \alpha) \mathcal{P}_c \approx\varepsilon\approx 0.95\varepsilon_{sp} \end{equation}

for $\alpha <0.5$. The energy production density $\mathcal {P}_m$ (green symbols in figure 8a) is higher in the dispersed phase for low volume fractions, while it is comparable to that of the carrier phase for $\alpha >0.1$. The energy dissipation rate per unit volume in the dispersed phase $\varepsilon _d$ (red symbols in figure 8a) is also larger at low volume fractions and decreases monotonically with increasing $\alpha$. The dissipation in the carrier phase, $\varepsilon _c$, also decreases, as it compensates for the energy transport $\mathcal {T}_m$ from the carrier flow towards the dispersed phase. The viscous transport (see blue symbols in figure 8b) is significantly lower than its pressure-induced counterpart, $\mathcal {T}^p_m$, although they exhibit similar behaviours: they first increase until $\alpha =0.1$, and then decrease to reach zero for a binary mixture. Note, again, that the total transport term is zero, i.e. the sum of $\mathcal {T}^p$ and $\mathcal {T}^{\nu }$ from both phases. The case $\alpha =0.5$ deserves a specific mention. In this case, production and dissipation in the two phases are equal, hence the transport term satisfies $\mathcal {T}_m=0$. Intuitively, it is not possible to define unambiguously a carrier and dispersed phase in binary mixtures; while the energy is locally transported from one phase to the other, the global average is zero for both pressure and viscous transfer.

Figure 8. Phase-averaged energy balance versus the dispersed phase volume fraction $\alpha$; see (2.7). In each plot, coloured triangles represent the dispersed phase ($m=d$), while circles represent the carrier phase ($m=c$). Each term is normalized by the single-phase energy dissipation $\varepsilon _{sp}$ computed for case SP2. (a) Energy production $\mathcal {P}_m$ and energy dissipation $\varepsilon _m$; (b) viscous energy transport $\mathcal {T}^\nu _m$ and pressure energy transport $\mathcal {T}^p_m$.

Finally, we analyse the p.d.f.s of velocity fluctuations $u_n=u/\sigma _u$, vorticity fluctuations $\omega _n=(|\omega _i|-\langle |\omega _i| \rangle )/\sigma _\omega$, and energy dissipation $\varepsilon _n=(\varepsilon -\langle \varepsilon \rangle )/\sigma _\varepsilon$, normalized by their standard deviations. In figure 9(a), we observe that while the p.d.f. of the velocity fluctuation remains symmetric, the tails deviate strongly from the typical pseudo-Gaussian behaviour of single-phase turbulence (Sreenivasan & Antonia 1997; Jimenez 2000; Ishihara, Gotoh & Kaneda 2009). Considering the vorticity in figure 9(b), no deviation is observed in the Gaussian core (as defined in Sreenivasan & Antonia 1997). However, the distributions of the multiphase flows depart strongly from the single-phase case in the tails. In particular, the exponentially decaying tails have a higher exponent in the case of emulsions, indicating more events with strong vorticity. Interestingly, while increasing the volume fraction does not influence the value of the exponent, increasing $\alpha$ induces deviations in the distributions already at lower values of $\omega _n$. We observe a similar behaviour for $\varepsilon _n$ in figure 9(c): the intermittency of the single-phase flow is amplified by the presence of the interface. As for the vorticity, departures from the single-phase distributions are observed at lower values of $\varepsilon$ when increasing the volume fraction $\alpha$. As a final general remark, the analysis of the p.d.f.s reveals that strong deviations are induced by the presence of the interface, already at low volume fractions, overall increasing the intermittent behaviour of the flow. As no collapse is observed for the normalized variables, it can be inferred that the small-scale statistics are affected by the presence of the interface.

Figure 9. P.d.f.s of (a) velocity fluctuations $u$, (b) vorticity $\omega$, and (c) dissipation $\varepsilon$. All quantities are normalized as standard score.

3.2. Influence of viscosity ratio

We consider now the influence of the viscosity ratio on the flow turbulence, i.e. cases BEx, V1x and V2x in table 1. The viscosity ratios analysed span the range $0.01<\mu _d/\mu _c<100$, while $We_\mathcal {L}=42.7$ for all cases. Two values of the volume fractions are considered, $\alpha =0.03$ (series V1x) and $\alpha =0.1$ (series V2x). A render of the two-fluid interface (corresponding to the value of the VOF function $\phi =0.5$) is shown in figures 10(a–c) for cases V11, BE1 and V14, respectively. As $\mu _d/\mu _c$ increases, larger droplets appear; at low viscosity ratios, we find a significantly higher number of droplets.

Figure 10. Render of the two-fluid interface (corresponding to the value of the VOF function $\phi =0.5$) for different values of the viscosity ratio $\mu _d / \mu _c$: (a) $0.01$, (b) $1$, (c) $100$. The vorticity fields are shown on the box faces on a planar view. All simulations are performed at $\alpha =0.03$ and $We_\mathcal {L}=42.6$.

We start by examining the Taylor-scale Reynolds number of the emulsion flows for the different viscosity ratios under investigation. Figures 11(a) and 11(b) show the variation of $Re_\lambda$ versus the viscosity ratio for the two volume fractions considered, $\alpha =0.03$ and $\alpha =0.1$. As expected, $Re_\lambda$ decreases with the viscosity ratio. Significant variations in $Re_\lambda$ are observed already for small volume fractions, the effects being amplified for $\alpha =0.1$. In the insets of figure 11, we can observe that $\lambda$ (i.e. the local variations of $k/\varepsilon$) does not increase linearly with $\mu _d/\mu _c$, indicating increased velocity fluctuations for the dispersed phase at lower viscosity.

Figure 11. Taylor-scale Reynolds number of the emulsion flows for the different viscosity ratios examined. $Re_\lambda$ is shown versus $\mu _d/\mu _c$. (a) Cases BE1 and V1x, with volume fraction $\alpha =0.03$; (b) cases BE2 and V2x, with $\alpha =0.1$. The insets show the evolution of $\lambda$ with the viscosity ratio.

To better quantify the variations of the flow gradients, we show the phase-averaged (see (A2)) enstrophy $\omega ^2_m$ in figure 12, normalized by the single-phase values from SP1. The viscosity ratio affects enstrophy strongly in the dispersed phase, while the magnitude in the carrier phase is almost constant. Further, smaller variations can be observed when changing the volume fraction from 0.03 to 0.1. For $\mu _d/\mu _c\leqslant 1$, the enstrophy in the dispersed phase goes approximately as $\omega ^2_c\propto -\log (\mu _d/\mu _c)$. As the viscosity of the dispersed phase becomes larger ($\mu _d>\mu _c$), $\omega _d$ decreases below the average value of the single-phase flow and tends towards zero, as high viscosity dampens velocity fluctuations in the dispersed phase. It is worth noting that for incompressible flows, the energy dissipation rate can be defined as $\varepsilon \equiv \nu |\omega _i|^2$; however, when phase averaging, the two formulations differ for by a term proportional to $\partial _{ii}p$.

Figure 12. Phase-averaged enstrophy $\omega ^2_m$ (normalized by its value in the single-phase case SP1) for different viscosity ratios $\mu _d/\mu _d$. Triangles indicate the dispersed phase ($m=d$), while circles indicate the carrier phase ($m=c$). (a) Results for $\alpha =0.03$; (b) results for $\alpha =0.1$.

We now discuss the influence of viscosity ratio on the compensated energy spectra, shown in figures 13(a,b) for $\alpha =0.03$ and $\alpha =0.1$, respectively. Similarly to previous observations for figure 5, the Hinze scale shows, to a good approximation, the pivoting point, below which energy increases with respect to the single-phase spectra. Differences in the inertial subrange are hardly observable for $\alpha =0.03$, while they become more prominent when the volume fraction is increased (see figure 13b). Analysis of the data for $\kappa <\kappa _H$ reveals that the simulations with a dispersed phase present less energy than the single-phase case. In the dissipative range, the trend emerges more clearly. As $\kappa >\kappa _H$, the less viscous the dispersed phase, the more energy is injected in the smaller scales. As discussed in the previous section, energy reduces at large scales and increases at small scales when increasing the volume fraction $\alpha$.

Figure 13. One-dimensional compensated energy spectra for (a) $\alpha =0.03$, and (b) $\alpha =0.1$, and different values of the viscosity ratio $\mu _b/\mu _c$. The vertical dotted line indicates the Hinze scale wavelength, $\kappa _H$.

Figure 14 shows the DSD for all configurations with different viscosity ratios. As for the data in § 3.1, we also display the $-3/2$ power law, which well describes the distribution of small droplets, and the $d^{-10/3}$ law from Garrett et al. (2000) for larger droplets. In this range, $d>d_H$, the $-10/3$ law is observed only in a limited region of the spectrum. As noted previously, this is most likely due to the low volume fraction considered. The variation of $\mu _d$ has an influence on large droplets, as higher viscosity in the dispersed phase increases the probability of formation of these large droplets. This was also observed qualitatively in figure 10 and confirms previous findings (Roccon et al. 2017).

Figure 14. P.d.f.s of the DSD for different values of $\mu$, at (a) $\alpha =0.03$, and (b) $\alpha =0.1$. The dashed line represents the $d^{-3/2}$ scaling from Deane & Stokes (2002), while the continuous black line shows the $d^{-10/3}$ law from Garrett et al. (2000).

We present the SBS energy budget for $\alpha =0.1$ in figure 15. The results for $\alpha =0.03$ show similar trends; see Appendix C for the details. Following the same scheme as in the previous section, we depict in figure 15(a) the energy balance for case V22, when $\mu _b/\mu _c=0.1$. Similarly to previous observations, the dispersed phase absorbs energy at large scales and redistributes it to small scales; that is, the presence of the interface provides an alternative path for energy transfer from small to large wavenumbers and no inverse cascade is observed. The nonlinear energy transfer $T$ (see figure 15b) displays a weak sensitivity to the viscosity ratio (almost negligible for $\alpha =0.03$, as shown in figure 27 in Appendix C). Thus the differences in the $Re_\lambda$ and energy spectra discussed above are not associated with an extension of the inertial range. For wavenumbers larger than that of the forcing, the nonlinear energy transfer is higher at large scales and lower at small scales than for the single-phase case.

Figure 15. Scale-by-scale energy budget for different viscosity ratios $\mu _d/\mu _c$ at $\alpha =0.1$. (a) Complete energy balance for case V22 with $\mu _d/\mu _c=0.1$; (b) nonlinear energy transfer $T$; (c) the term $\mathcal {S}_\sigma$ associated with the surface tension; and (d) the energy dissipation transfer function $\mathcal {D}$.

Figure 15(c) shows the energy transport due to the surface tension term, $\mathcal {S}_\sigma$. As $\mu _d/\mu _c$ increases, the wavelength where the positive energy transport is maximum shifts to larger scales. This behaviour is possibly due to increased coalescence for high $\mu _d/\mu _c$ (as discussed later in this section). We also observe that with decreasing viscosity ratio, $\mu _d/\mu _c<1$, the curves tend to collapse, as the data for $\mu _d/\mu _c=0.1$ and $\mu _d/\mu _c=0.01$ are approximately overlapping. At the injection scale, $\kappa =2$, almost all cases behave similarly. At intermediate wavelengths, the lower the viscosity ratio, the higher the energy absorbed by the surface tension forces. As previously observed, the Hinze scale represents, to a good approximation, the point where the energy transfer towards small scales by the surface tension term $\mathcal {S}_\sigma$ is maximum. All these observations apply to the two values of $\alpha$ considered; see also Appendix C.

A note should be made on the flow with the highest viscosity ratio: in this case, the energy is not transferred down to the dissipative range. A qualitative explanation is given by the following scenario. When the interface interacts with a sufficiently large vortex in the carrier phase, it tends to deform, and in doing so, it absorbs energy through the work of the interfacial stresses. The deformation of the interface induces shear in the dispersed phase that is opposed by viscous forces. A higher viscosity in the dispersed phase will therefore dump larger and more energetic structures, reducing the energy available at small scales.

Finally, figure 15(d) shows the transfer function of the energy dissipation term $\mathcal {D}$. We observe that simulations with higher viscosity of the dispersed phase dissipate more energy at large scales, hence dumping turbulence in the inertial range, as expected for more viscous flows. This trend is maintained until the dissipative range, where, instead, a lower viscosity ratio produces higher dissipation. This causes the apparently paradoxical situation that despite there being limited energy transport by the nonlinear terms, dissipation is still active at small scales because of the energy brought by the interfacial stresses; this may suggest the need for a specific definition of dissipative range for multiphase flows.

Next, we discuss the influence of the viscosity ratio on the phase-averaged energy budget, shown in figure 16 for $\alpha =0.1$. As for the SBS balance, the same analysis for $\alpha =0.03$ can be found in Appendix C, as the variation of volume fraction does not change the underlying physical process significantly. The production density (green symbols in figure 16a) shows only slight variations with viscosity ratios in the carrier phase, whereas it increases in the dispersed phase when its viscosity increases; in particular, $\mathcal {P}_d<\mathcal {P}_c$ when $\mu _d<\mu _c$. A similar trend is observed for the dissipation rate (red symbols in figure 16a), when the differences between dispersed and carrier phases become more evident. In this case, the dissipation in the dispersed case increases with its viscosity until $\mu _d/\mu _c=100$, when it decreases because of the lower energy transferred to smaller scales inside the droplets. The transport terms, $\mathcal {T}^\mu$ and $\mathcal {T}^p$ in figure 16(b), indicate that energy is always transferred from the carrier to the dispersed phase. Both terms increase in magnitude when decreasing the viscosity ratio, indicating that energy needs to be supplied to the dispersed phase to sustain turbulence when viscous forces are increasing. The pressure transport is the preferential path for energy transfer from the carrier to the dispersed phase for low and moderate values of $\mu$. For the case with largest viscosity of the dispersed phase, the transfer due to pressure forces becomes lower than that associated with viscous forces.

Figure 16. Phase-averaged energy balance versus the emulsion viscosity ratio; see definitions of the terms in (2.7). Coloured triangles represent the dispersed phase ($m=d$), while circles are used for the carrier phase ($m=c$). Each term is normalized by the single-phase energy dissipation $\varepsilon _{sp}$ computed for case SP2. The energy production $\mathcal {P}_m$ and energy dissipation $\varepsilon _m$ are reported in (a), while viscous energy transport $\mathcal {T}^\nu _m$ and the pressure energy transport $\mathcal {T}^p_m$ are shown in (b).

Finally, we consider the effect of the viscosity ratio on the p.d.f.s of velocity, vorticity and dissipation rate; see figure 17 for $\alpha =0.1$, while data for $\alpha =0.03$ can be found in Appendix C. A low viscosity in the dispersed phase generates larger velocity fluctuations (see also figure 12), hence the tails of the p.d.f.s are more evident for small values of $\mu _d/\mu _c$ in figure 17(a). Interestingly, when the viscosity ratio increases above unity, velocity fluctuations in the dispersed phase are quenched, the standard deviation decreases, and the statistics are closer to those of the single-phase reference case. The distributions of the normalized vorticity are shown in figure 17(b). As for the velocity fluctuations, a higher viscosity in the dispersed phase decreases the intermittency, and the distributions approach the single-phase values. For $\mu _d/\mu _c<1$, intermittency increases and the tails of the distribution are more evident; nonetheless, they can still be fitted with decaying exponentials. As observed previously for varying $\alpha$, the pseudo-Gaussian part of the vorticity p.d.f. collapses for all cases.

Figure 17. P.d.f.s of (a) velocity fluctuations $u$, (b) vorticity $\omega$, and (c) energy dissipation. All quantities are normalized by their standard deviations. The data pertain to cases V2x in table 1, with $\alpha =0.1$.

The p.d.f. of the energy dissipation shows almost no alteration between cases with different viscosity ratios, while intermittency is increased strongly with respect to the single-phase case. Due to the normalization with $\sigma _\varepsilon$, the curve collapse indicates that variations induced by $\mu _d/\mu _c$ of the small-scale dynamics are negligible.

To conclude, the turbulence is significantly affected by variations of the viscosity ratio already at small volume fractions. Higher viscosity of the dispersed phase dampens the small-scale structures because of higher viscous dissipation at all wavelengths. For emulsions with viscosity of the dispersed phase lower than that of the carrier phase, the activity at small scales increases and so does intermittency. The surface tension term $\mathcal {S}_\sigma$ contributes significantly to the transfer of energy to the smallest scales in this case.

3.3. Influence of Weber number

The influence of the surface tension coefficient, expressed through the large-scale $We_{\mathcal {L}}$ number, is examined in this subsection. As discussed in the literature (Komrakova et al. 2015; Roccon et al. 2017; Mukherjee et al. 2019), the combination of volume fraction, surface tension coefficient and energy injection scale, $\mathcal {L}$, has to be chosen accurately because the HIT configuration is very sensitive to coalescence. Furthermore, high $We_{\mathcal {L}}$ may generate an excess of unresolved droplets, affecting the results significantly. Therefore, all the simulations discussed in this subsection are performed at $\alpha =0.03$, while the forcing is maintained at $\kappa _0=2$. The cases discussed in this subsection are BE1 and the series W1x with reference to table 1, covering a significant range of $We_{\mathcal {L}}$, from 10.6 to 106.5. In figure 18, we show a render of the flow for different values of $We_\mathcal {L}$. As expected, at low $We_\mathcal {L}$, we observe the appearance of large liquid structures due to higher surface tension forces. At high $We_\mathcal {L}$, on the other hand, the dispersed phase undergoes severe fragmentation (see Appendix B for the DSD convergence study of case W13). The presence of small droplets resulting from fragmentation can be observed in all cases.

Figure 18. Render of the two-fluid interface (corresponding to the value of the VOF function $\phi =0.5$) for different values of the Weber number $We_\mathcal {L}$: (a) $10.6$, (b) $42.6$, and (c) $106.5$. The vorticity fields are shown on the box faces on a planar view. All simulations are performed at $\alpha =0.03$ and $\mu _d/\mu _c=1$.

We start by studying the global behaviour of the flow through the integral quantities in figure 19. In figure 19(a), $Re_\lambda$ shows an almost linear increment with both the Taylor scale $\lambda$ and $We_{\mathcal {L}}$ (represented with colours). Unlike previous observations in § 3.1 where the increase of $Re_{\lambda }$ was due mostly to local variations of the ratio $k/\varepsilon$, decreasing surface tension also lowers the volume-averaged energy dissipation, as shown in the inset. These findings are in agreement with the results on turbulent emulsions in Rosti et al. (2020). As the viscosity ratio $\mu _d=\mu _c$ is constant, the decrease of the dissipation is caused by lower enstrophy levels, as can be appreciated from the data for the carrier phase in figure 19(b). The behaviour of the enstrophy of the dispersed phase is less intuitive, exhibiting non-monotonicity; this will be addressed later when discussing the phase-averaged energy balance.

Figure 19. (a) $Re_\lambda$ versus the Taylor scale $\lambda$; the inset shows the global energy dissipation $\varepsilon$ (normalized by its single-phase value) as a function of $We_{\mathcal {L}}$. (b) Phase-averaged enstrophy versus the Weber number $We_{\mathcal {L}}$; coloured triangles represent the dispersed phase ($m=d$), while circles indicate data pertaining to the carrier phase ($m=c$).

Figure 20(a) shows the compensated energy spectra at different $We_\mathcal {L}$. As we are varying the surface tension, the Hinze scale varies in each case (see vertical dotted lines of corresponding colours). As mentioned before, the Hinze scale defines with good approximation the spectra pivoting point.

Figure 20. (a) One-dimensional compensated energy spectra for different large-scale Weber numbers $We_\mathcal {L}$; the wavelengths corresponding to the Hinze scale of each spectra are plotted with vertical dotted lines of corresponding colours. (b) DSD for different $We_\mathcal {L}$; the Hinze scale $d_H$ is reported with dotted lines of corresponding colours. The continuous black line represents the region where the $-10/3$ power law applies.

As discussed previously, energy is reduced at larger scales in the inertial range, and increases at smallest scales. With increasing $We_\mathcal {L}$, higher energy is observed at high wavelengths.

Figure 20(b) shows the DSD for all the $We_\mathcal {L}$ under investigation. As for figure 20(a), we show the Hinze scale for each case with vertical dotted lines. Again we observe that the $-10/3$ power law from Garrett et al. (2000) provides a reasonable description for the largest droplets, $d>d_H$; as we increase $\sigma$, i.e. low $We_{\mathcal {L}}$, larger droplets may appear, as expected because of the increased cohesion forces. In this case, the energy required to break up large droplets is available only in large eddies. As their turnover time is of the order of $\mathcal {T}$, large droplet breakup becomes a rare event and the distributions are more noisy, so it is more difficult to identify a clear trend. In addition, by reducing $We_\mathcal {L}$, the distribution becomes more irregular for $d< d_H$ as most of the dispersed phase is in large droplets.

The effects of $We_{\mathcal {L}}$ can be described better by the scale-by-scale analysis, shown in figure 21. The complete energy balance is shown for case W11 ($We_\mathcal {L}=10.6$) in figure 21(a). Unlike the cases shown previously for $We_\mathcal {L}=42.6$, the surface tension energy transfer $\mathcal {S}_\sigma$ is more uniform through the different scales, and its effects are less evident globally. To deepen the analysis, we display the nonlinear energy transfer function $T$ for each case at different $We_\mathcal {L}$ in figure 21(b). At the injection wavelength $\kappa =2$, no major differences are observed when varying the surface tension. At small wavelengths, $\kappa >2$, we observe that the energy transfer by the nonlinear term increases with $We_\mathcal {L}$, compensating for the effect of the energy absorption from the surface tension. The energy transfer at smaller scales, after the peak, increases with $\sigma$, approaching the values of the single-phase flow.

Figure 21. Scale-by-scale energy budget for different large-scale Weber numbers $We_\mathcal {L}$. (a) Full energy balance for case W11, with $We_\mathcal {L}=10.6$; (b) energy transfer $T$ due to the nonlinear term; (c) energy flux $\mathcal {S}_\sigma$ associated with the surface tension term; (d) energy dissipation rate $\mathcal {D}$.

The energy transfer via the interfacial stresses, $\mathcal {S}_{\sigma }$, is shown in figure 21(c). The energy is again absorbed at large scales and distributed at small scales. Flows with small $We_\mathcal {L}$ absorb more energy at small wavenumbers, and the transmission of energy (i.e. positive $\mathcal {S}_\sigma$) is smeared over a higher range of scales, hence the peak ($\max (\mathcal {S}_\sigma )$) is also less evident. For all $We_{\mathcal {L}}$ investigated, the surface tension term $\mathcal {S}_{\sigma }$ transfers energy also within the dissipative range at small scales, where the transport from nonlinear terms has become negligible.

The energy dissipation $\mathcal {D}$, in figure 21(d), decreases with $We_{\mathcal {L}}$ at large and intermediate scales, as energy is absorbed partially by $\mathcal {S}_{\sigma }$. However, the amplitude of the dissipation rates becomes almost independent of the Weber number at the smallest scales. Further, as previously observed, the presence of the dispersed phase delays the onset of the dissipative range.

The phase-averaged energy balance from simulations with different Weber numbers is shown in figure 22. Both production and dissipation (figure 22a) are found to decrease in the carrier phase when increasing $We_\mathcal {L}$, while the former increases and then decreases in the carrier phase. Possibly, this can be related to the DSD: decreasing the droplet size increases the internal dissipation, which may explain the behaviour at the lower Weber number examined. On the other hand, high deformability decreases the dissipation close to the interface, which may explain the decrease at the largest $We_\mathcal {L}$. For all values of $We_\mathcal {L}$ considered, the dispersed phase extracts kinetic energy from the carrier phase, as $\mathcal {T}_c>0$ (figure 22b). The decrease of surface tension forces results in a monotonic decrease of the viscous transfer and an increase of the pressure transport for the dispersed phase. Consistently, dissipation is always higher in the dispersed phase, while it decreases in the carrier phase with increasing $We_\mathcal {L}.$

Figure 22. Phase-averaged energy balance versus the emulsion Weber number; see definitions of the terms in (2.7). Coloured triangles represent the dispersed phase ($m=d$), while circles indicate data pertaining to the carrier phase ($m=c$). Each term is normalized by the single-phase energy dissipation $\varepsilon _{sp}$ computed for case SP2. The energy production $\mathcal {P}_m$ and energy dissipation $\varepsilon _m$ are reported in (a), while viscous energy transport $\mathcal {T}^\nu _m$ and the pressure energy transport $\mathcal {T}^p_m$ are shown in (b).

Finally, the p.d.f.s for velocity, vorticity and dissipation are shown in figure 23(a–c). Strong variations are induced in all p.d.f.s, showing that indeed a more rigid interface favours the appearance of extreme events. Since a more deformable interface offers lower resistance to the propagation of velocity disturbances from one phase to the other, less modification of the p.d.f.s with respect to the single-phase case can be expected at higher $We_{\mathcal {L}}$ (see also Rosti et al. 2020). This is indeed observed in all p.d.f.s, where the distributions are seen to approach the single-phase distribution when increasing $We_\mathcal {L}$. Nevertheless, rare events are still evident also at the largest Weber considered, especially for the energy dissipation. Vorticity shows, again, that the pseudo-Gaussian part of the distribution is identical for all $We_\mathcal {L}$, while the exponentially decaying tails display strong variations.

Figure 23. P.d.f.s of (a) velocity fluctuations $u$, (b) vorticity $\omega$, and (c) energy dissipation. All quantities are normalized by their standard deviations. The data pertain to cases W1x, BE1 and SP1 in table 1.

4. Conclusions

In this work we discuss how volume fraction, viscosity ratio and Weber number influence HIT in emulsions. The analyses are performed at different levels of detail, spanning from phase-averaged balances to SBS energy transfer in spectral space. Some observations are common to all configurations and highlight some fundamental physical effects introduced by the dispersed phase. Here, we first consider these different aspects and then discuss the modulation introduced by the variation of material properties.

4.1. Spectra and SBS energy balance

In all simulations with a dispersed phase, the energy decreases at large scales and increases at small scales, corroborating previous findings (Ten Cate et al. 2004; Perlekar et al. 2014; Dodd & Ferrante 2016; Mukherjee et al. 2019; Rosti et al. 2020; Olivieri et al. 2020a). Interestingly, this behaviour applies to both solid and liquid dispersed phases in HIT. Furthermore, the pivoting point of the energy spectra is found to be described, with a good approximation, by the Hinze scale. This has also been observed in binary mixtures (Perlekar et al. 2014) and emulsions (Mukherjee et al. 2019), and is extended here to several operating conditions.

In general, the mechanisms of energy transport are modified as follows. The transfer by the nonlinear advection terms decreases, as the surface tension forces absorb energy at large scales. In an emulsion, energy is transferred to small scales also by the surface tension force, well within the dissipative range of the corresponding single-phase flow, forcing the viscous dissipation to be active at even smaller scales. No inverse cascade has been observed in the present simulations.

The general idea, according to which coalescence and breakup are responsible for modifications of the energy spectra, seems to only partially explain our observations. In fact, according to this hypothesis, significant deviations should be observed when comparing spectra for different volume fractions. Here, instead, we observe the largest deviations in the energy spectra, in particular at small scales, when varying the viscosity ratio.

From the discussion above, it appears that the amount of kinetic energy at small wavenumbers may not be used to predict modifications of the various terms of the SBS budget. In fact, it may be hypothesized that the SBS budget shows a more direct connection with the DSD. In particular, the SBS term $\mathcal {S}_\sigma$ changes with the total surface area and interface deformation. In other words, strong topological changes, breakup and coalescence, change the energy transport by surface tension. Future studies may address the relation between the energy at small scales and the energy cascade by numerically inhibiting, reducing or controlling coalescence and breakup (e.g. using front-tracking methods).

4.2. Effects of the dispersed phase on the dissipative range

The classical ‘far’ dissipative range ($\kappa \sim \kappa _{max}$), where both nonlinear energy transport and energy dissipation of the SBS budget are zero, is lost when a dispersed phase is introduced. In multiphase flows, despite the nonlinear energy transfer vanishing at certain small scales, the energy dissipation does not because energy is brought to these smaller scales by the action of the surface tension. As discussed above, energy dissipation is thus forced to extend towards smaller scales, overall increasing the range of wavelengths where there is activity. In other words, this increased activity at small scale translates also into an extension of the dissipative range, with the nonlinear transport substituted by the surface tension transport. This is an important finding in the present work, which was not clearly observed in previous studies with diffuse-interface methods (Perlekar 2019).

It is important to understand how the scaling in the inertial range might be affected by these modifications of the dissipation range. From a practical viewpoint, the results in Appendix B show that increasing the mesh resolution does not result in significant alterations of the inertial range, indicating that a relevant analysis of the inertial range dynamics is still possible even in simulations where the surface tension terms are slightly under-resolved at small scale. Nevertheless, resolving the dissipative range is important for a complete discussion of the SBS budget and e.g. the DSD; understanding turbulence at small scales in multiphase flows remains, therefore, a relevant question also from a fundamental point of view.

4.3. Flow intermittency

We have observed that the presence of a dispersed phase increases intermittency, unless the dispersed phase is highly viscous. The probability of detecting rare events increases, mainly for energy dissipation and vorticity, as shown here by the p.d.f. analysis. This finding represents a significant observation for multiphase turbulence and needs to be addressed further in future studies.

In particular, at higher volume fractions and constant $\mu _d/\mu _c$ and $We_\mathcal {L}$, the exponent describing the distribution tail exponential decay is independent of the volume fraction $\alpha$. The onset of the exponential tail (hence the probability of an extreme event) is, on the other hand, affected by $\alpha$, proving that these events are occurring mostly at the interface. This is a confirmation of the observations in Dodd & Ferrante (2016) on the increased energy dissipation at the interface. The variation of the exponential tail for both energy dissipation rate and vorticity at different $\mu _d/\mu _c$ and $We_\mathcal {L}$ reveals that intermittency is affected significantly by the fluid properties. In cases with high $\mu _d$ and low surface tension, the vorticity intermittency is attenuated and similar to the single-phase cases. On the other hand, the dissipation seems to be always affected by the multiphase nature of the flow.

4.4. Droplet statistics

In all the conditions analysed, the DSD shows both the $-3/2$ exponential scaling from Deane & Stokes (2002) for the small droplets and the $-10/3$ exponential scaling from Garrett et al. (2000) for the larger ones, confirming and extending the previous findings of Mukherjee et al. (2019) to a significant number of different configurations. Moreover, employing a VOF approach, and its known mass conserving properties, allows us to extend the $-3/2$ scaling to significantly small droplets, which was not observed clearly in previous studies.

The power law $d^{-10/3}$ well describes the distributions of larger droplets when the volume fraction is below $10\,\%$, with only a small loss in accuracy for higher values of $\alpha$, in agreement with the assumption of negligible coalescence in Garrett et al. (2000). Although this power law was obtained under the assumption of a dilute dispersed phase, recent works based on a diffuse-interface approach report the same scaling in the presence of coalescence (Mukherjee et al. 2019; Soligo et al. 2019). However, Deike et al. (2016) estimate through accurate sharp-interface simulations a similar exponent, $-3$, so that it might be difficult to have a clear distinction on the different effects. Finally, we show that the estimate of the Hinze scale as a transition point between the two power laws is less accurate for $\alpha >0.1$, suggesting that different model coefficients may be needed when coalescence is relevant.

4.5. Role of the fluid properties

Our analyses demonstrate that the volume fraction $\alpha$ is the parameter that modifies mostly the energy fluxes in the flow; yet, increasing the volume of the dispersed phase does not change the underlying physics. This is documented notably in § 3.1, where we show that the amount of total interface area determines the energy transport across scales. Moreover, the simulation data reveal that the energy transfer via surface tension forces is enhanced at low viscosity ratios, while high viscosity in the carrier phase inhibits the propagation of vortices through the interface, hence reducing the overall energy transport. Changing the Weber number amounts to modulating the pivoting frequency below which energy transfer through surface tension is directed towards smaller scales. In particular, as the dispersed phase is less deformable, the energy absorption from the dispersed phase occurs at larger scales, and turbulence is reduced progressively. In fact, as the surface tension increases, more energy is required to deform the droplets, an energy that can be found only in large-scale eddies.

To study the role of the viscosity ratio (see § 3.2), we consider values ranging from $10^{-2}$ (a value typical of bubbles) to $10^{2}$ (typical of droplets). The analysis reveals that for $\mu _d/\mu _c\leqslant 1$, $Re_\lambda$ increase significantly, due to the lower viscosity in the dispersed phase. The scale-by-scale energy budget shows that the interfacial and nonlinear transport terms are not strongly affected at these low viscosity ratios. For $\mu _d/\mu _c> 1$, on the other hand, the turbulence in the dispersed phase is reduced, which implies a significantly smaller $Re_\lambda$, below the value of the single-phase case. In these cases, the energy transfer induced by the interfacial stresses is reduced significantly, suggesting that large differences may be found in liquid–gas and gas–liquid emulsions. The DSD does not show strong differences, although larger droplets are more likely to be generated by a more viscous dispersed phase. Note, as discussed above, that the viscosity ratio has a significant impact on the flow intermittency.

Finally, we have examined the role of the large-scale Weber number $We_{\mathcal {L}}$. At low $We_{\mathcal {L}}$, coalescence is more likely to occur, hence there is a higher probability to find large droplets. Nevertheless, the Hinze scale proves to be an accurate estimation of the transition between $-3/2$ and $-10/3$ for all the cases analysed. Changing $We_\mathcal {L}$ and thus the DSD also affects the energy transport across scales by the surface tension forces. Specifically, when decreasing $We_\mathcal {L}$, the energy injection from interfacial tension moves to larger scales.