Droplet-turbulence interactions and quasi-equilibrium dynamics in turbulent emulsions

TL;DR

This paper uses advanced simulations to study how droplets interact with turbulence in emulsions, revealing a quasi-equilibrium cycle of coalescence and breakup, and highlighting the unique flow topologies caused by droplet interfaces.

Contribution

It demonstrates the capability of the pseudopotential lattice-Boltzmann method for stable, long-term emulsion simulations and uncovers the dynamic limit-cycle behavior of emulsions in turbulence.

Findings

Droplet size distribution follows a d^{-10/3} scaling.

Emulsions evolve into a quasi-equilibrium cycle of coalescence and breakup.

Flow topology in emulsions shows dominant vortex compression and axial straining.

Abstract

We perform direct numerical simulations (DNSs) of emulsions in homogeneous, isotropic turbulence using a pseudopotential lattice-Boltzmann (PP-LB) method. Improving on previous literature by minimizing droplet dissolution and spurious currents, we show that the PP-LB technique is capable of long, stable simulations in certain parameter regions. Varying the dispersed phase volume fraction $\phi$, we demonstrate that droplet breakup extracts kinetic energy from the larger scales while injecting energy into the smaller scales, increasingly with higher $\phi$, with the Hinze scale dividing the two effects. Droplet size ($d$) distribution was found to follow the $d^{-10/3}$ scaling (Deane & Stokes 2002). We show the need to maintain a separation of the turbulence forcing scale and domain size to prevent the formation of large connected regions of the dispersed phase. For the first time, we show that turbulent emulsions evolve into a quasi-equilibrium cycle of alternating coalescence and breakup dominated processes. Studying the system in its state-space comprising kinetic energy $E_k$, enstrophy $\omega^2$ and the droplet number density $N_d$, we find that their dynamics resemble limit-cycles with a time delay. Extreme values in the evolution of $E_k$ manifest in the evolution of $\omega^2$ and $N_d$ with a delay of $\sim0.3\mathcal{T}$ and $\sim0.9\mathcal{T}$ respectively (with $\mathcal{T}$ the large eddy timescale). Lastly, we also show that flow topology of turbulence in an emulsion is significantly more different than single-phase turbulence than previously thought. In particular, vortex compression and axial straining mechanisms become dominant in the droplet phase, a consequence of the elastic behaviour of droplet interfaces. Revised and extended version now published in the Journal of Fluid Mechanics: https://doi.org/10.1017/jfm.2019.654

Tables

Table 1. Table 1 : Simulations parameters for all cases. Here viscosity ν 𝜈 \nu and interfacial tension γ 𝛾 \gamma are in lattice units [lu], along with length and time measured as multiples of Δ ​ x Δ 𝑥 \Delta x and Δ ​ t Δ 𝑡 \Delta t . The density and viscosity ratio between the components is kept at unity. The turbulence forcing is distributed over the range of wavenumbers from k a subscript 𝑘 𝑎 k_{a} to k b subscript 𝑘 𝑏 k_{b} centered at k f subscript 𝑘 𝑓 k_{f} . The fluid densities are initialized to ρ α , β in = 4.0 superscript subscript 𝜌 𝛼 𝛽 in 4.0 \rho_{\alpha,\beta}^{\mathrm{in}}=4.0 ρ α , β out = 0.77 superscript subscript 𝜌 𝛼 𝛽 out 0.77 \rho_{\alpha,\beta}^{\mathrm{out}}=0.77 . The average kinetic energy ⟨ E k ⟩ = ( ∑ k E ​ ( k ) ) / N delimited-⟨⟩ subscript 𝐸 𝑘 subscript 𝑘 𝐸 𝑘 𝑁 \left\langle{E_{k}}\right\rangle=(\sum_{k}E(k))/N , and the average rate of energy dissipation ⟨ ϵ ⟩ = ( ∑ k 2 ​ ν ​ k 2 ​ E ​ ( k ) ) / N delimited-⟨⟩ italic-ϵ subscript 𝑘 2 𝜈 superscript 𝑘 2 𝐸 𝑘 𝑁 \left\langle{\epsilon}\right\rangle=(\sum_{k}2\nu k^{2}E(k))/N .

Sim	$N$	$\nu$	$G_{\alpha \beta }$	$\varphi$	$A$	$k_a,k_f,k_b$	$\gamma$	$⟨E_k⟩$	$⟨ϵ⟩$	${\tau }_k$	$R{e}_{\lambda }$	$\eta$
SP	${256}^3$	0.0047	-	-	0.0005	2.0	-	$1.8\times {10}^{-3}$	$5.0\times {10}^{-7}$	97	95	0.7
P1	${256}^3$	0.0047	0.015	0.01	0.0005	$1,2,8$	0.017	$2.0\times {10}^{-4}$	$2.21\times {10}^{-8}$	$461$	$51$	$1.47$
P2	${256}^3$	0.0047	0.015	0.06	0.0005	$1,2,8$	0.017	$2.0\times {10}^{-4}$	$2.10\times {10}^{-8}$	$474$	$53$	$1.49$
P3	${256}^3$	0.0047	0.015	0.15	0.0005	$1,2,8$	0.017	$1.7\times {10}^{-4}$	$1.93\times {10}^{-8}$	$493$	$47$	$1.52$
P4	${256}^3$	0.0047	0.015	0.2	0.0005	$1,2,8$	0.017	$1.5\times {10}^{-4}$	$1.75\times {10}^{-8}$	$518$	$45$	$1.56$
P5	${256}^3$	0.0047	0.015	0.45	0.0005	$1,2,8$	0.017	$1.3\times {10}^{-4}$	$1.65\times {10}^{-8}$	$534$	$39$	$1.58$
T1	${256}^3$	0.0047	0.015	0.10	0.00025	$1,1.5,8$	0.017	$8.4\times {10}^{-5}$	$4.90\times {10}^{-9}$	$980$	$44$	$2.14$
T2	${256}^3$	0.0047	0.015	0.10	0.0005	$1,1.5,8$	0.017	$2.4\times {10}^{-4}$	$1.87\times {10}^{-8}$	$502$	$64$	$1.54$
T3	${256}^3$	0.0047	0.015	0.10	0.00075	$1,1.5,8$	0.017	$4.6\times {10}^{-4}$	$4.78\times {10}^{-8}$	$313$	$78$	$1.22$
T4	${256}^3$	0.0047	0.015	0.10	0.001	$1,1.5,8$	0.017	$6.5\times {10}^{-4}$	$8.36\times {10}^{-8}$	$237$	$84$	$1.06$
T5	${256}^3$	0.0047	0.016	0.10	0.001	$1,1.5,8$	0.04	$6.6\times {10}^{-4}$	$7.70\times {10}^{-8}$	$247$	$90$	$1.08$
D1	${128}^3$	0.0047	0.015	0.15	0.0005	$1,1.5,6$	0.017	$1.2\times {10}^{-4}$	$1.80\times {10}^{-8}$	$511$	$34$	$1.55$
D2	${256}^3$	0.0047	0.015	0.15	0.0005	$1,3,6$	0.017	$1.1\times {10}^{-4}$	$1.77\times {10}^{-8}$	$514$	$30$	$1.56$
D3	${384}^3$	0.0047	0.015	0.15	0.0005	$2,4.5,8$	0.017	$1.1\times {10}^{-4}$	$1.81\times {10}^{-8}$	$509$	$30$	$1.55$
D4	${512}^3$	0.0047	0.015	0.15	0.0005	$3,6,9$	0.017	$1.1\times {10}^{-4}$	$1.87\times {10}^{-8}$	$501$	$30$	$1.54$
D5	${512}^3$	0.0047	0.015	0.2	0.0005	$1,1.5,6$	0.017	$4.2\times {10}^{-4}$	$1.80\times {10}^{-8}$	$511$	$118$	$1.55$

Equations

$$d_{\max }=0.725({\rho }^c/\gamma {)}^{-3/5}{ϵ}^{-2/5}$$

$$N_x\gg \mathcal{L}\gg d\gg \eta \gg \zeta$$

$$N_x>\mathcal{L}\gg d\gg \eta \approx \zeta$$

$$f_i^{\sigma }\left(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t\right)$$

$$f_i^{\mathrm{eq},\sigma }=w_i\rho \left(1+\frac{{\mathbf{u}}^{\mathrm{eq}}\cdot {\mathbf{c}}_i}{c_s^2}+\frac{{\left({\mathbf{u}}^{\mathrm{eq}}\cdot {\mathbf{c}}_i\right)}^2}{2c_s^4}-\frac{{\mathbf{u}}^{\mathrm{eq}}\cdot {\mathbf{u}}^{\mathrm{eq}}}{2c_s^2}\right)$$

$${\mathbf{u}}^{\mathrm{eq}}={\mathbf{u}}^{\prime }+\frac{{\tau }^{\sigma }{\mathbf{F}}^{\sigma }}{{\rho }^{\sigma }}$$

$${\mathbf{u}}^{\prime }=\frac{{\sum }_{\sigma }\frac{{\rho }^{\sigma }{\mathbf{u}}^{\sigma }}{{\tau }^{\sigma }}}{{\sum }_{\sigma }\frac{{\rho }^{\sigma }}{{\tau }^{\sigma }}}$$

$${\mathbf{u}}^{\sigma }=\frac{1}{{\rho }^{\sigma }}\sum _i{f}_i^{\sigma }{\mathbf{c}}_i$$

$${\mathbf{F}}^{\mathrm{SC},\sigma }\left(\mathbf{x}\right)=-G_{\sigma \overline{\sigma }}{\psi }^{\sigma }\left(\mathbf{x}\right)\sum _{\sigma \ne \overline{\sigma }}{\psi }^{\overline{\sigma }}\left(\mathbf{x}+{\mathbf{c}}_i\Delta t\right){\mathbf{c}}_i{w}_i\Delta t$$

$$p=c_s^2\sum _{\sigma }{\rho }^{\sigma }+\frac{c_s^2\Delta t^2}{2}\sum _{\sigma ,\overline{\sigma }}{G}_{\sigma \overline{\sigma }}{\psi }^{\sigma }{\psi }^{\overline{\sigma }}$$

$$F_x^{\sigma }$$

$$F_y^{\sigma }$$

$$F_z^{\sigma }$$

$$A(k)=A\exp \left(-\frac{{\left(k-k_f\right)}^2}{c}\right)$$

$$P=P_1+P_2=\frac{1}{2}\overline{f_k{f}_k}\Delta t+\overline{u_k{f}_k}$$

$$\mathcal{L}\sim \frac{2\pi }{k_f{k}_{\min }}=\frac{N_x}{k_f}$$

$$\lambda ={\left(\frac{15\nu {u^{\prime }}^2}{ϵ}\right)}^{1/2}$$

$$R{e}_{\lambda }=\frac{u^{\prime }\lambda }{\nu }$$

$${\tau }_k={\left(\frac{ϵ}{\nu }\right)}^{-1/2}$$

$$S(k,t)=\frac{{\sum }_k|{\hat{\varphi }}_{\mathbf{k}}{|}^2}{{\sum }_{k}1}$$

$$L(t)=2\pi \frac{{\sum }_{k}S(k,t)}{{\sum }_{k}kS(k,t)}$$

$$\mathrm{Corr}=\frac{\overline{⟨E_k(t)⟩⟨{\omega }^2(t+\delta t)⟩}}{\overline{⟨E_k⟩}\overline{⟨{\omega }^2⟩}}$$

$$\cos (\theta )=\frac{\nabla \psi }{|\nabla \psi |}\cdot \frac{\mathbf{\omega }}{|\mathbf{\omega }|}$$

$${\lambda }^3+P{\lambda }^2+Q\lambda +R=0$$

Full text

1 Introduction

An emulsion consists of a dense suspension of droplets of one fluid (the dispersed phase) suspended in another fluid (the continuous phase), and is formed due to turbulent mixing of these two immiscible fluids. Emulsions are found (both desirably and undesirably) in a wide range of industries. For instance, in food processing, diverse products depend on the stability and texture of emulsions (McClements, 2015). In biotechnology, emulsions can serve as miniature laboratories where living cells can be compartmentalized into individual droplets (Griffiths & Tawfik, 2006). They are also known to cause various losses in crude oil production (Kokal, 2005), or to the contrary, enable enhanced oil recovery (Banat, 1995). The formation of an emulsion requires shearing of droplets which can occur both in laminar and turbulent flows, although the latter may be a more common occurrence. Turbulent emulsions can be said to form a particular class of droplet laden turbulent flows where there is close interplay between turbulence and the dynamics of the dispersed fluid. An accurate description of these systems hence involves the dynamics of deforming interfaces, while allowing for coalescence and breakup of droplets, resolution of a range of length and time scales of turbulent flow and the possible presence of surface active agents (surfactants) that can alter the interfacial dynamics.

The primary effect of turbulence on droplets in these flows is to cause fragmentation, where an initially large connected volume of the dispersed phase is broken into smaller droplets. Under sustained turbulence, there is a supposed equilibrium between coalescence and breakup which leads to a droplet spectrum around a theoretical maximum stable diameter, known as the Hinze scale (Hinze, 1955). This droplet distribution is also known to follow a $d^{-10/3}$ slope (Deane & Stokes, 2002), where $d$ is the droplet diameter. The dispersed phase influences the continuous phase turbulence by drawing turbulent kinetic energy (TKE) from the flow, which partially goes into the difference between the surface energy of parent and daughter droplets, while the rest is stored in the deformation of interfaces. This reduces the effective turbulent kinetic energy (TKE), which has consequences on the turbulence cascade and spectrum, noticeably at scales comparable to droplet sizes. Coalescing droplets in turn set finer flow structures into motion, where interfacial tension releases the energy stored in droplet deformations back as TKE into the flow (Dodd & Ferrante, 2016) at scales smaller than the droplet sizes.

Literature review

In this paper, we are interested in simulating an emulsion under sustained turbulence. Simulations are key here, as experiments as yet are incapable of revealing the underlying dynamical processes - limited by emulsions being optically opaque, while interfacial dynamics is inherently three-dimensional, whereby experiments give mostly statistical or phenomenological results. There have been only a few numerical studies devoted to turbulent emulsion dynamics, some of which have been detailed in the recent review by Elghobashi (2019) on DNS simulations of turbulent flows laden with droplets or bubbles. We refer interested readers to it for a general overview, while we shall discuss the current state of simulating turbulent emulsions, highlighting those aspects that we intend to address in our work.

In one of the first studies, Derksen & Van den Akker (2007) simulated a turbulent liquid-liquid dispersion using a free-energy based lattice-Boltzmann (LB) method. They modeled a fluid packet as it passes by the impeller in a stirred vessel, hence experiencing a burst of turbulence, before entering a quiescent zone. They show evolution of the droplet distribution in the dispersion under first constant and then decaying turbulence, also reporting the modification to the kinetic energy spectra at a crossover scale.

Perlekar et al. (2012) simulated droplet breakup in homogeneous, isotropic turbulence using a pseudopotential (PP) LB method, showing that the distribution of droplet diameters has a finite width around the Hinze scale. Since Hinze’s criterion does not account for droplet coalescence or coagulation, deviation from it was found at higher volume fractions. Further, droplet breakup was attributed to peaks in the local energy dissipation rate. The study reported on the method being originally incapable of attaining steady state simulations due to droplet dissolution, which was remedied by a mass correction scheme to artificially re-inflate droplets which helped maintain a steady volume fraction (Biferale et al., 2011). Later, Perlekar et al. (2014) simulated turbulent spinodal decomposition to show coarsening arrest in a symmetric binary fluid mixture (which is similar to an emulsion, although the morphology is distinctly different). The presence of turbulence was shown to inhibit the coarsening dynamics at droplet sizes larger than the Hinze scale.

Skartlien et al. (2013) simulated a surfactant laden emulsion under weak turbulence ($Re_{\lambda}\leq 20$) using a free-energy LB method, and reproduced a $d^{-10/3}$ droplet distribution as predicted by Deane & Stokes (2002) for air bubbles above the Hinze scale. They did not find any influence of the surfactant in altering the coalescence rates in the considered range of surfactant activity and turbulence intensity. Also using a free-energy LB method Komrakova et al. (2015a) simulated turbulent liquid-liquid dispersions at varying volume fractions, focusing on the resolution of droplets with respect to the Kolmogorov scale. They found that droplet dissolution was a significant issue, which made it impossible to obtain a steady state droplet distribution at low phase fractions, while at higher phase fractions ($\phi>0.2$), despite breakup, most droplets coalesce to form a single connected region with multiple smaller satellite droplets. Increasing the resolution of the Kolmogorov scale remedied droplet dissolution to some extent, and a log-normal droplet distribution was shown from transient simulations, as has been experimentally found for turbulent liquid-liquid dispersions (Pacek et al., 1998; Lovick et al., 2005). The multiphase energy spectra could not be reproduced due to spurious currents which caused unphysical energy gain at high wavenumbers, whose magnitude was found to be close to the turbulent velocity scale $u^{\prime}$.

In their detailed study on droplet-turbulence interaction, Dodd & Ferrante (2016) simulated a large number of initially spherical droplets ($\phi=0.05$) in decaying homogeneous, isotropic turbulence using a mass conserving volume-of-fluid method. They considered a wide range of density and viscosity ratios between the droplet and carrier fluid, and showed an enhanced rate of energy dissipation for increasing droplet Weber number ($W\!e$). Introducing the TKE equations, they show that breakup and coalescence act as source and sink terms of TKE. Roccon et al. (2017) studied the influence of viscosity on breakup and coalescence in a swarm of droplets ($\phi=0.18$) in wall bounded turbulent flow using a coupled Cahn-Hillard Navier-Stokes solver. They report a slight drag reduction in the flow due to the presence of droplets, and show that a higher interfacial tension or droplet viscosity favours coalescence, and the number of droplets rapidly decreases to $1-10\%$ of its initial value. At low viscosity, where breakup dominates, around $50\%$ of the droplets remain separated and their sizes follow Hinze’s $\left\langle D\right\rangle\propto W\!e^{-3/5}$ criterion.

Recently, using a mass conserving level-set method, Shao et al. (2018) studied interface-turbulence interactions in droplet breakup simulations. They showed that vortical structures tend to align with large scale interfaces before breakup. They also show that there is a slight increase in axial straining and vortex compression in the flow topology in the presence of droplets, in comparison to single-phase turbulence.

Our study

In this study, we resolve several of the issues faced in previous work, and report new findings from direct numerical simulations of turbulent emulsions. We use the PP-LB method for a multicomponent fluid system without phase change to simulate the formation of a dispersion starting from a single large drop in the center of a periodic box. PP-LB is well suited for the simulation of multiphase flows comprising deformable droplets due to the spontaneous formation of interfaces emerging from simplified interparticle repulsion forces (Shan & Chen, 1993, 1994; Shan & Doolen, 1995). The method has been used before for simulating droplets in turbulence (Perlekar et al., 2012, 2014; Albernaz et al., 2017), an in general has been applied many more times for non-turbulent flows with droplets. It allows for coalescence and breakup to occur naturally, all without the need for interface tracking or models for film drainage. However, it comes with a caveat that due to interfaces being diffuse, coalescence is favourable when interfaces overlap. This makes the resolution of interface width relative to droplet sizes, i.e. the Cahn number, an important criterion (Shardt et al., 2013). The diffuse interface also leads to dissolution of small droplets as has been noted before (Perlekar et al., 2012; Komrakova et al., 2015a). We show that droplet dissolution can be limited to a minor effect in certain parameter ranges, and that a mass correction scheme as used in Biferale et al. (2011); Perlekar et al. (2012) is not requisite for simulating droplets in turbulence.

Additionally, multiphase LB simulations suffer from spurious currents ($u^{\mathrm{sp}}$) which are velocities arising from anisotropy in the discretization of inter-particle forces. Although $u^{\mathrm{sp}}$ in pseudopotential LB have been shown to be much lower than in conventional finite volume techniques like the volume-of-fluid method (Kamali & Van den Akker, 2013; Mukherjee et al., 2018), in the free-energy LB method they were found strong enough to dominate the multiphase kinetic energy spectra at high wavenumbers (Komrakova et al., 2015a). Further, in LB, the characteristic fluid velocity (here the large scale velocity $\mathcal{U}$) should be kept smaller than the lattice speed of sound $c_{s}$, such that the flow Mach number $M\!a=\mathcal{U}/c_{s}$ is low (where traditionally $M\!a<0.3$ is considered incompressible) and hence the flow being simulated obeys the incompressible Navier-Stokes equations. Hence, the velocities should scale as $c_{s}>\mathcal{U}\gg u^{\mathrm{sp}}$, which we maintain in our work.

We simulate a dispersion in a periodic box, employing a forcing scheme to generate homogeneous, isotropic turbulence. Here we particularly study the influence of varying the dispersed phase volume fraction ($\phi$) and turbulence intensity ($Re_{\lambda}$) on the properties of the dispersion formed. We show the influence of the dispersed phase on the multiphase kinetic energy spectra which has not been systematically presented before, or was not possible due to the limitations of the numerical method (Komrakova et al., 2015a). We show that $\phi$, $Re_{\lambda}$ and the interfacial tension $\gamma$ together determine the dispersion morphology, and that droplets of a particular characteristic length can be generated by varying these parameters. Investigating local flow topology, we show that the effect of the dispersed phase is significant and more pronounced than previously stated (Shao et al., 2018), with a sharp increase in vortex compression and axial straining in the droplet regions. We also present, for the first time, an analysis of the equilibrium dynamics of a droplet laden isotropic turbulent flow, showing that the system evolution in its state-space is akin to a limit-cycle with alternating dominance of coalescence and breakup as the system oscillates between different dispersion morphologies.

Length scales

Through this study we highlight a few considerations that have not been discussed in previous work and are crucial to simulating droplets in turbulence. First is numerically resolving to a sufficient degree the several length scales that govern different aspects of these simulations. The main length scale is the typical droplet diameter, which can be taken to be the Hinze scale and is given as (Hinze, 1955)

[TABLE]

where $\rho^{c}$, $\gamma$ and $\epsilon$ are the carrier fluid density, interfacial tension and rate of energy dissipation, respectively, while it is now accepted that the local variations in $\epsilon$ (intermittency) set a local Hinze scale, and an entire spectrum of droplets centered around $d_{\mathrm{max}}$ tends to arise. A closely associated length scale is the interface width $\zeta$, which in physical systems can be of the order of nanometers for micron to millimeter size droplets. However, as a limitation of our simulation technique (and every other diffuse interface method), the interface width extends over a few computational grid cells. The ratio between $\zeta$ and the droplet diameter $d$ is termed the Cahn number $Ch=\zeta/d$ (Komrakova et al., 2015b), and extreme values of $Ch$ are undesirable. Hence the relative separation between $d$ and $\zeta$ needs to be considered.

Next, the two length scales characterizing turbulence are the energy injection scale $\mathcal{L}$ which is determined by the forcing scheme, and the smallest (or Kolmogorov) scale $\eta$ which is determined by the viscosity $\nu$ and the dissipation rate $\epsilon$. A wide separation between $\mathcal{L}$ and $\eta$ means a higher Reynolds number $Re$, which can be expressed as $Re\approx\left(\mathcal{L}/\eta\right)^{4/3}$. A final length scale of importance in simulations is the size of the simulation domain, which along one spatial direction can be considered to be $N_{x}$, and this is generally chosen to be close to $\mathcal{L}$. As droplets will break up due to extension under turbulent stresses, the domain size $N_{x}$ should be sufficiently larger than the maximum droplet elongation before breakup to yield meaningful results (particularly for simulations on periodic domains, where large droplets would begin to interact with images of themselves). Here a particular caveat is also the simplistic description of highly deformed droplets, where an equivalent droplet diameter $d=(6V/\pi)^{1/3}$ gives the impression of $N_{x}\gg d$, whereas in the form of long, slender filaments, droplets can linearly extend across the entire domain. This can give rise to elongated droplets that remain connected due to periodicity, and this is more prone to occur at high volume fractions under weak turbulence, as for instance can be seen in Skartlien et al. (2013).

Comparing these length scales, the required spatial separation between them for simulating droplets in the inertial range, at least from a stance of reasoning, would follow as

[TABLE]

while $N_{x}>\mathcal{L}$ may also be sufficient, and most studies currently are limited to $N_{x}\approx\mathcal{L}$. Also, $d$ can vary over a range of values, extending upto $d\sim\eta$ if the Kolmogorov scale is over-resolved. Upon conceding to limitations of modeling, current simulations can at best reproduce

[TABLE]

We try to maintain such a separation of scales in the study, except that we have $\zeta>\eta$. Lastly, having $\eta>d$ would mean sub-Kolmogorov droplets. These droplets can also deform and breakup due to the action of viscous stresses instead of inertial stresses (Elghobashi, 2019).

We begin with a description of the numerical method in section 2, followed by a brief validation of the turbulence forcing scheme. We then present results from turbulent emulsions in section 4, for varying volume fraction in section 4.2 and varying turbulence intensity in section 4.3. Section 4.4 discusses the importance of sufficient resolution of the largest scales and section 4.5 shows the influence of the turbulence forcing wavenumber on the dispersion morphology. Finally, in section 5 we discuss some general results regarding emulsion dynamics, with the quasi-equilibrium limit-cycle presented in section 5.1, droplet-vorticity alignment in section 5.2 and influence of droplets on local flow topology in section 5.3, after which we end with the conclusions.

2 Numerical Method

2.1 Lattice Boltzmann Method

Each component $\sigma\in\left\{\alpha,\beta\right\}$ obeys the standard LBGK equation with a single relaxation time which can be written as (Krüger et al., 2017)

[TABLE]

where $f_{i}^{\sigma}$ is the distribution function of component $\sigma$ along the discrete velocity direction $\mathbf{c}_{i}$. Here $\tau^{\sigma}$ is the lattice relaxation time towards local equilibrium which relates to the macroscopic component viscosity $\nu^{\sigma}=c_{s}^{2}(\tau^{\sigma}-1/2)$ where $c_{s}=1/\sqrt{3}$ is the lattice speed of sound (the mixture viscosity is a more complex expression when the components have different $\tau$). The equilibrium distribution $f_{i}^{\mathrm{eq},\sigma}$ is given by the local Maxwellian as

[TABLE]

where $w_{i}$ are the LB weights in each direction $i$, and $\mathbf{u}^{\mathrm{eq}}$ is the equilibrium velocity which is given as

[TABLE]

Here $\mathbf{F}^{\sigma}$ incorporates all the forces (here the inter-component interactions and the turbulence forcing), into the common fluid velocity $\mathbf{u}^{\prime}$ between the two components which is given as

[TABLE]

where $\mathbf{u}^{\sigma}$ is the bare component velocity. This is calculated in its usual form

[TABLE]

For details see Succi (2001); Krüger et al. (2017). The inter-component interaction force, $\mathbf{F}^{\mathrm{SC}}$, is modeled using the method of Shan & Doolen (1995), which can be written as

[TABLE]

where $\psi^{\sigma}$ is the pseudopotential function for component $\sigma$ and in this study we have chosen $\psi^{\sigma}=\rho^{\sigma}$ (while other definitions are possible). This force between the components is kept to be repulsive, hence the interaction strength parameter $G_{\sigma\overline{\sigma}}$ should have a positive value. It should be noted that the fluids remain partially miscible, and essentially the final composition consists of $\alpha-$rich and $\beta-$rich regions, while a small amount of one component remains dissolved in the other. A higher magnitude of $G_{\sigma\overline{\sigma}}$ results in lower solubility and gives rise to a higher interfacial tension. The equation of state for this multicomponent system is (Krüger et al., 2017)

[TABLE]

Lastly, the interfacial tension $\gamma$ can be calculated using the Laplace law $\Delta p=2\gamma/r$, where $\Delta p$ is the pressure difference across the interface of a spherical droplet.

The simulations here have been performed on a $D3Q19$ lattice, i.e. a three-dimensional lattice with a set of $19$ discrete velocity directions. Further, the lattice spacing $\Delta x$ and time step $\Delta t$ are both set equal to $1$, and consequently all quantities are expressed in dimensionless lattice units [lu].

2.2 Turbulence Forcing

To generate and sustain turbulence in the fluid, a constant source of energy is required, which is constantly being dissipated by viscosity at the smallest scales (i.e. the Kolmogorov scales). This is done by setting the largest scales of flow into motion, and if the fluid viscosity is low enough, these large structures become unstable and give rise to successively smaller scales. One of the ways to achieve this numerically is by employing a low wavenumber spectral forcing, as given by Alvelius (1999), while alternative techniques could also be used (Eswaran & Pope, 1988; Rosales & Meneveau, 2005). This forcing was also implemented by Ten Cate et al. (2004) in LB to simulate the response of clouds of spherical solid particles to homogeneous isotropic turbulence. A very similar form of the forcing is used by Perlekar et al. (2012), which is constructed directly in real space but could be made to have a similar effective spectral form as (Ten Cate et al., 2004, 2006), albeit with less control over output parameters, as we do in this study. The forcing is divergence free by construction and can be written as

[TABLE]

Here $\rho^{\mathrm{tot}}=\sum_{\sigma}\rho^{\sigma}$ is the total density considering both components and each $\phi_{i}(k)$ is a unique random phase. Alternatively, $\phi_{i}(k)$ can be evolved as a stochastic process, as done in Perlekar et al. (2012), but in our approach $\phi_{i}(k)$ (and hence the forcing) varies as white noise in time. This ensures that the force is not related to any timescale of turbulent motion, and is a choice also made in Ten Cate et al. (2006). The force is distributed over a small range of wavenumbers $k_{a}\leq k\leq k_{b}$, while the contribution of each of these wavenumbers is determined by $A(k)$ which centers the Gaussian around $k_{f}$ in Fourier space, given as

[TABLE]

where $k_{f}$ is the central forcing wavenumber, $c$ is a width over which to distribute the force amplitude and is set to $c=1.25$, and $A$ is a forcing magnitude. This method ensures that there is a dominant central wavenumber $k_{f}$ (which can also be a fraction) in the forcing scheme, while neighbouring wavenumbers also contain some energy, which makes the scheme more stable (Ten Cate et al., 2006). Lastly, the total power input to the fluid can be written as the sum of two terms as follows

[TABLE]

where the two terms are the force-force and force-velocity correlations respectively, and $u_{k},f_{k}$ refer to the volumetric velocity and force fields. The force-velocity correlation, $P_{2}$, should be [math] to avoid an uncontrolled growth of energy in the fluid (Alvelius, 1999), and it is achieved by varying the force term at each time step. This is computationally expensive, hence some studies (Ten Cate et al., 2004, 2006) vary the force by choosing randomly from a pre-computed set of force fields at each time step. This was found to introduce a non-zero contribution from the $P_{2}$ term, where the steady state kinetic energy was roughly $10$ times larger than with a unique random force at each time step - hence in this study we adhere to the latter approach.

The largest scale in the system is given by the domain size $N_{x}$, which sets the minimum wavenumber $k_{\mathrm{min}}=2\pi/N_{x}$. All other wavenumbers are integer multiples of $k_{\mathrm{min}}$, with the maximum wavenumber being $k_{\mathrm{max}}=k_{\mathrm{min}}N_{x}/2=\pi$. The smallest scale of turbulence (Kolmogorov scale) is calculated as $\eta\sim\left(\nu^{3}/\epsilon\right)^{1/4}$ where $\nu$ and $\epsilon$ are the kinematic viscosity and energy dissipation rate respectively. The criterion for a resolved DNS simulation is that $k_{\mathrm{max}}\eta>1$ (Moin & Mahesh, 1998), and the Kolmogorov scale should obey $\eta>0.318$ [lu] (Ten Cate et al., 2006). We shall mention the forcing wavenumber $k_{f}$ and the wavenumber bounds as multiples of $k_{\mathrm{min}}$ in this study. For a central forcing wavenumber $k_{f}$, the associated large scale length then becomes

[TABLE]

Further, the Taylor microscale is calculated as

[TABLE]

where $u^{\prime}$ is the root mean square velocity along one direction, and $u^{\prime}_{x}=u^{\prime}_{y}=u^{\prime}_{z}$ in isotropic turbulence. The rate of energy dissipation $\left\langle{\epsilon}\right\rangle$ can be found in two ways, as $\epsilon\approx\nu\left\langle{\omega^{2}}\right\rangle\approx\sum_{k}2\nu k^{2}E(k)/N_{x}^{3}$ where $\left\langle{\omega^{2}}\right\rangle$ is the average enstrophy and $E(k)$ is the kinetic energy spectrum. Using $\lambda$, the Taylor Reynolds number is calculated as

[TABLE]

Lastly, the Kolmogorov timescale is given as

[TABLE]

For eddies in the inertial range with a size $l$, the velocity $u(l)$ and timescale $\tau(l)$ are determined uniquely by $\epsilon$ and $l$ alone as $u(l)=(\epsilon l)^{1/3}\sim\mathcal{U}(l/\mathcal{L})^{1/3}$ and $\tau(l)=(l^{2}/\epsilon)^{1/3}\sim\mathcal{T}(l/\mathcal{L})^{2/3}$, where $\mathcal{L}$, $\mathcal{T}$ and $\mathcal{U}$ are the characteristic length, time and velocity of the largest eddies (with $\mathcal{T}=\mathcal{L}/\mathcal{U}$). We consider $\mathcal{U}\approx{\left\langle{E_{k}}\right\rangle}^{1/2}$ as the largest eddies contain most of the kinetic energy, and generally $u^{\prime}<\mathcal{U}$. The characteristic velocity at a particular length scale can also be found from the kinetic energy spectrum as $u(l)\approx\sqrt{E(k_{l})}$ where $k_{l}=2\pi/l$.

3 Single-phase turbulence

We begin with a single-phase turbulence simulation to show that the forcing scheme is able to maintain a statistically stationary turbulent flow (simulation “SP” in table 1) and to compare it with results available in literature. A domain of $256^{3}$ lattice nodes representing a length $(2\pi)^{3}$ is initialized with a uniform initial density of $\rho^{\alpha}=4.0$ [lu]. The relaxation time of the two fluids is set to $\tau=0.5141$ which gives a viscosity of $\nu=0.0047$ [lu] (Perlekar et al. (2012) use a similar value with $\tau=0.515$), which is a low enough viscosity to sustain turbulence while still being numerically stable. The forcing is concentrated around $k_{f}=2k_{\mathrm{min}}$ and is distributed in the range of $k=k_{\mathrm{min}}$ to $8k_{\mathrm{min}}$. Further, $A=0.0005$, which generates a turbulent flow with a Taylor microscale of $\lambda=13$ [lu], $Re_{\lambda}=95$, $\tau_{k}=97$ [lu], $\eta=0.7$ [lu] ($k_{\mathrm{max}}\eta=2.2$) and $\left\langle{\epsilon}\right\rangle\approx 5\times 10^{-7}$ [lu], which are calculated a posteriori. The simulation is performed for $10^{5}$ $\Delta t$, which corresponds to $1000\tau_{k}$.

Figure 1 shows the evolution of $\left\langle{E_{k}}\right\rangle$ and $\left\langle{\omega^{2}}\right\rangle$ which attain their steady state values around $100\tau_{k}$ and continue to oscillate around this value. The crests and troughs of the $\left\langle{E_{k}}\right\rangle$ evolution show up in the $\left\langle{\omega^{2}}\right\rangle$ evolution with a slight delay (as seen in the inset of figure 1, where the quantities have been normalized with their time averaged values over the latter $3/4$th of the simulation duration). This has been observed before, and ascribed to the energy cascading mechanism (Pearson et al., 2004; Biferale et al., 2011) while Tsinober (2009) acknowledges this feature without invoking a cascade.

Figure 2 shows typical velocity and enstrophy field snapshots from a planar cross-section in the center of the domain at $500\tau_{k}$. The velocity field shows motions across various scales, while the enstrophy field (which is the square of the vorticity and relates directly to the rate of energy dissipation) shows typical small scale localized structures. Also note that $\omega^{2}$ assumes values as much as $10$ times the average $\left\langle{\omega^{2}}\right\rangle$ (while at higher $Re_{\lambda}$, more extreme values are found), showing that intermittency is well reproduced in the simulations. This patchy structure of enstrophy (and hence dissipation) is an important factor to consider in simulations of turbulent dispersions, as the local rate of energy dissipation sets the local maximum stable droplet diameter.

The kinetic energy spectrum is shown in figure 3, along with a benchmark spectrum from the Johns Hopkins Turbulence Database (Li et al., 2008) for a homogeneous isotropic turbulence simulation with $Re_{\lambda}=433$ (on a grid of $1024^{3}$, generated with a spectral solver). The energy $E(k)$ has been normalized by the total energy $\sum_{k}E(k)$, and the wavenumber is normalized to show multiples of $k_{\mathrm{min}}$, which is done to compare the two spectra. A well developed inertial range is seen to exist, following the $k^{-5/3}$ spectral slope, which falls off around $k=30k_{\mathrm{min}}$ in our simulation. Lastly, in this simulation $u^{\prime}=0.034$ [lu], and since the speed of sound is $c_{s}=1/\sqrt{3}$ [lu], the flow Mach number is $M\!a=0.06$ which is well within the incompressibility limit.

4 Turbulent emulsions

4.1 Simulation setup

The turbulent emulsion simulations are initialized with two fluids, which we denote by $\alpha$ (the carrier fluid) and $\beta$ (the droplet fluid). For a chosen volume fraction $\phi$ of fluid $\beta$, a single spherical droplet (a $\beta$-rich region) is initialized in the center of the domain which is otherwise $\alpha$-rich. The droplet density is denoted by $\rho_{\beta}^{\mathrm{in}}$, i.e the density of $\beta$ in the $\beta$-rich region, while $\rho_{\beta}^{\mathrm{out}}$ denotes the dissolved amount of component $\beta$ in the $\alpha$-rich region (i.e. the continuous phase), and likewise for component $\alpha$. Further, $\rho_{\beta}^{\mathrm{avg}}$ is used to refer to the average density of component $\beta$ in the entire domain. During the simulation, these density values can change to some extent depending on the $G_{\alpha\beta}$ parameter, though due to the symmetry of the model we have $\rho_{\beta}^{\mathrm{in}}/\rho_{\alpha}^{\mathrm{in}}\approx 1$ and $\rho_{\beta}^{\mathrm{out}}/\rho_{\alpha}^{\mathrm{out}}\approx 1$, which well represents many oil in water emulsions. We also keep $\nu_{\beta}/\nu_{\alpha}=1$ (with $\nu=0.0047$ [lu]). Spurious velocities ($u^{\mathrm{sp}}$) in these simulations have been limited to values sufficiently smaller than the physical velocity, so that their influence on the results is negligible. Typically, the large scale velocity $\mathcal{U}\sim O(10^{-2})$ while $u^{\mathrm{sp}}_{\mathrm{mean}}\sim O(10^{-4})$ and $u^{\mathrm{sp}}_{\mathrm{max}}\sim O(10^{-3})$ (for the range of $G_{\alpha\beta}$ values used in this study). Given that the speed of sound in these simulations $c_{s}=1/\sqrt{3}$, we maintain that $u^{\mathrm{sp}}\ll\mathcal{U}\ll c_{s}$, which is in line with our recent findings for emulsion droplets simulated with PP-LB (Mukherjee et al., 2018; Berghout & Van den Akker, 2019).

We carried out three sets of simulations, the details of which are mentioned in table 1. In all these simulations, the turbulence force is applied starting at $t=0$. The turbulence energy density $\left\langle{E_{k}}\right\rangle$ in an emulsion, for the same forcing amplitude, can be an order of magnitude lower than in single-phase turbulence. The Kolmogorov scale values have been calculated using the scaling $\eta\approx(\nu^{3}/\overline{\left\langle{\epsilon}\right\rangle})^{1/4}$ where $\overline{\left\langle{\epsilon}\right\rangle}$ is the spatio-temporally averaged dissipation rate (with $\overline{\left\langle{.}\right\rangle}$ denoting time averaging after the first quarter of the simulation time, during which the flow is well developed). We report $\eta$ upto two decimal places that follow from this scaling. The three sets are divided as follows

• •

Set 1 (P1-P5): In these simulations, only the dispersed phase volume fraction has been changed (from $\phi=0.01$ to $\phi=0.45$). Here $\eta$ is found to increase in simulations P1-P5, which is because the turbulence forcing scale $\mathcal{L}$ remains the same while $Re_{\lambda}$ decreases, hence reducing the separation between the largest and smallest scales.

• •

Set 2 (T1-T5): In these simulations, the turbulence force amplitude is varied to change $Re_{\lambda}$ (at a fixed volume fraction $\phi=0.10$). For case T5, the interfacial tension has also been increased. Due to increasing $Re_{\lambda}$ in these simulations, since $\mathcal{L}$ is kept constant, $\eta$ is found (as expected) to decrease.

• •

Set 3 (D1-D5): In these simulations, the domain size is increased while keeping the forcing lengthscale $\mathcal{L}$, amplitude and volume fraction ($\phi=0.15$) fixed, which keeps the turbulence energy density (or $Re_{\lambda}$) fixed. An additional simulation, D5, has been performed where the turbulence intensity and volume fraction have been increased for comparison with case D4. For all cases, $\eta$ remains almost constant as $Re_{\lambda}$ is kept constant by varying $\mathcal{L}$ (so that the ratio $\mathcal{L}/\eta$ is constant). In simulation D5, $Re_{\lambda}$ is increased fourfold in comparison to D1-D4, yet $\eta$ is the same as the increase in $Re_{\lambda}$ is achieved by the added scale separation due to a fourfold decrease in the forcing wavenumber in D5 ($k_{f}=1.5$) as opposed to D4 ($k_{f}=6.0$).

To study the droplet characteristics in these simulations, we segment the droplets in space (also known as clustering) by thresholding the droplet density field at a cutoff value $\rho^{c}/\rho_{\beta}^{\mathrm{in}}=0.57$ (which is effectively the density along the interface where $\rho^{c}\approx\rho_{\alpha}=\rho_{\beta}$) based on the algorithm used in Siebesma & Jonker (2000). This allows us to identify and mark all lattice points within individual droplets, which gives the droplet volume $V$, which in turn is used to calculate an effective diameter $d=(6V/\pi)^{1/3}$. Estimating the surface area of these droplets, which are in voxel form, requires more care. Often, the ‘GNU triangulation surface’ (GTS) library (Popinet & Jones, 2004) is used in studies due to its efficient surface splitting operations (without the need for volumetric droplet segmentation). However, it was not used in this study as it did not provide a straightforward way of identifying droplets cut-off at domain edges due to periodicity (an issue implicitly resolved by our segmentation algorithm). Also, the GTS library was found to underpredict the surface area of a sphere by around $10\%$. Instead, we use the method proposed by Windreich et al. (2003) (originally developed for medical MRI data) to calculate surface area directly from voxels using a look-up table which divides surface voxels into 9 classes, and each class has a weighted contribution to the surface area. Using only the first 4 of these 9 classes, the area estimation error for a sphere was found to decrease to $1\%$, which was sufficiently accurate for our study.

4.2 Effect of volume fraction

We now show results from simulations with varying dispersed phase volume fractions $\phi\in\left\{0.01,0.06,0.15,0.2,0.45\right\}$ under identical turbulence forcing conditions (corresponding to P1-P5 in table 1). These simulations are performed for $10^{5}$ time steps. Figure 4 shows the dispersion formation process at various time instances starting from the initial spherical droplet of component $\beta$ shown as the iso-surfaces representing $\rho_{\beta}=\rho_{\alpha}$. The droplet begins to deform under the turbulent stresses, eventually breaking up to form a dispersion with a characteristic distribution.

Of the various volume fractions considered, $\phi=0.06$ and $0.15$ are most emulsion-like, i.e. they have a profusion of small droplets with a few large connected filaments. At $\phi=0.01$, the dispersed phase is too dilute to be considered an emulsion, although the droplet dynamics is interesting as the number of droplets $N_{d}$ and their characteristic diameters $d$ is small, and hence most of the droplets remain dispersed with relatively few coalescence events, and when droplets do coalesce, they break up soon after. At $\phi\geq 0.2$, most of the fluid volume remains connected, which is aggravated by the enhanced coalescence inherent to diffuse interface methods (Komrakova et al., 2015a; Roccon et al., 2017). This in turn is due to insufficient resolution of the interface with respect to the droplet sizes (Shardt et al., 2013), an effect we discuss more in depth in section 4.4. At higher turbulence intensity, the large connected regions can be expected to break into smaller droplets, and any coalescence will generate droplets of sizes larger than the maximum stable diameter, which will again breakup.

Phase fraction evolution

Figure 5 shows the evolution of the dispersed phase volume fraction $\phi$ normalized by the initial volume fraction $\phi_{0}$. There is a clear decrease over time (upto around $100\tau_{k}$) in the relative volume fraction, beyond which the value plateaus to a level around which it continues to oscillate (this will be confirmed subsequently from simulations T1-T5 in section 4.3 which were performed for a five times longer duration). This relative reduction in $\phi$ is more pronounced at lower $\phi$ values (up to around $30\%$) than at higher $\phi$ (around $2-5\%$). Note that this is not a mass conservation issue, as the total component mass is perfectly conserved in the system, and only the amount of component $\beta$ present as the dispersed phase reduces, which gets dissolved in the $\alpha$-rich (continuous phase) region. This is also why the relative decrease in $\phi$ is strongest for $\phi=0.01$, as the dissolution of $\beta$ into the continuous phase is provided by a very low number of droplets.

The reason for the reduction in $\phi$ is twofold. First is the dissolution of small droplets due to a finite interface width, which is an issue inherent to most diffuse interface methods. This can be characterized by the Cahn number $Ch$. Hence, if $Ch\sim O(1)$ (or greater), the droplet becomes unstable and is prone to dissolution. On the other hand, Shardt et al. (2013) showed for droplet collision in shear flow that coalescence is inhibited with decreasing $Ch$ number. In the limit of $Ch\to 0$, coalescence would cease to occur, while at larger $Ch$ numbers, coalescence is favoured. These considerations mandate having a finite $Ch$ number in the range $0<\ll Ch\ll 1$(for all droplet sizes in the system) for achieving steady state simulations while allowing for both coalescence and breakup. The second reason for the reduction in $\phi$ is its sensitivity to the segmentation threshold. In appendix A we demonstrate that only this result, i.e. the evolution of the volume fraction, depends on the choice of the segmentation threshold. Part of the droplet phase fraction goes into constituting the increased interfacial region (i.e. roughly the total surface area of all droplets $S_{A}$ multiplied by the interface width $\zeta$). Slightly varying the segmentation threshold to lower values (so that it is closer to $\rho_{\beta}^{\mathrm{out}}$), the volume fraction loss is reduced (which may indicate that $\rho^{\mathrm{c}}\neq(\rho^{\mathrm{out}}+\rho^{{in}})/2$), although the exact choice of $\rho^{\mathrm{c}}$ does not change our results. Further, the reduction in $\phi$ is also not monotonic, as mass of component $\beta$ dissolved in the $\alpha$-rich region can eventually accumulate inside other droplets.

Droplet dissolution can be a debilitating numerical issue, where for instance Biferale et al. (2011); Perlekar et al. (2012) had to resort to artificially inflating droplets to maintain a constant phase fraction and Komrakova et al. (2015a) reported that they could not attain steady state simulations with the free-energy LB method at low volume fractions as all droplets dissolved away into the continuous phase. In our PP-LB simulations, this issue is due to an interplay of three main factors - (i) the liquid-liquid repulsion $G_{\alpha\beta}$ which keeps the two components demixed, (ii) the turbulence intensity which breaks large droplets into smaller ones and (iii) the phase fraction which at low values makes $\rho_{\beta}^{\mathrm{out}}\approx\rho_{\beta}^{\mathrm{avg}}$ (i.e. at low $\phi$, phase segregation can become weaker). Despite being present, droplet dissolution is limited to a minor effect in our simulations. More precisely, the PP-LB method employed in this study can be used to reasonably simulate certain regions of the turbulent emulsions parameter space where droplet dissolution is not significant. Namely, for a given turbulence intensity ($Re_{\lambda}$), there will be a critical lower bound on the interfacial tension $\gamma_{c}$ such that droplets with $\gamma>\gamma_{c}$ can be simulated. For increasing $Re_{\lambda}$, $\gamma_{c}$ would increase as well, and its exact dependence on $Re_{\lambda}$ could be investigated by numerically mapping the phase space which is out of the scope of the current study. Similarly, there will be a lower bound on the value of $\phi$, below which all droplets will dissolve due to weak phase segregation when $\rho_{\beta}^{\mathrm{out}}\approx\rho_{\beta}^{\mathrm{avg}}$. Considering these related effects, we restrict ourselves to a parameter range where we can attain long, stable simulations to collect meaningful statistics pertaining to the droplet coalescence and breakup equilibrium.

Droplet number density evolution

Figure 6 shows the evolution of the number of droplets ($N_{d}$) in the system for varying $\phi$. $N_{d}$ begins to increase following the first breakup events around $25\tau_{k}$ and rises steadily to its characteristic value around $75\tau_{k}$, around which it continues to oscillate. These oscillations in $N_{d}$ are indicative of competing coalescence and breakup dynamics. The falls in the $N_{d}$ evolution profiles are due to coalescence events, which generate droplets of large sizes that are unstable. These droplets then break up under turbulent stresses and $N_{d}$ increases again. Breakup is delayed for $\phi=0.01$ as compared to the other cases and $N_{d}$ only begins to increase around $50\tau_{k}$. This is because the size of initial droplet is much smaller ($\sim 64$ [lu]) than the forcing wavelength ($\sim 128$ [lu]), and the droplet starts to advect initially, as seen from figure 4. When smaller scales are generated (around $50\tau_{k}$, as can be seen from the enstrophy evolution in figure 1), the droplet begins to shear and break. The evolution of $N_{d}$ does not show large fluctuations for $\phi=0.01$ due to relatively fewer coalescence and breakup events in this case, which is because the droplets are smaller and more distant from each other than in higher $\phi$ cases.

Although $\phi=0.15$ and $0.2$ simulations have a larger volume of fluid $\beta$, the number of droplets generated is lower than $\phi=0.06$. This is because of a higher propensity for coalescence in these systems which generates large connected regions and smaller satellite droplets. This is most prominently seen for $\phi=0.45$, where $N_{d}$ is even lower than $\phi=0.01$, as most of the fluid forms extended filaments that remain connected across the periodic boundaries. Increasing the turbulence intensity can be expected to generate more droplets at higher $\phi$, and hence for a given $Re_{\lambda}$, there will be a specific $\phi$ that maximizes the number of droplets formed and hence produce a more emulsion-like droplet size distribution.

Droplet size distribution

Figure 7 shows the distribution of the equivalent droplet diameter $d=(6V/\pi)^{1/3}$ (where $V$ is the droplet volume) for varying $\phi$ (calculated with $25000-35000$ droplets identified between times $75\tau_{k}$ to $200\tau_{k}$, sampled at each $\tau_{k}$). Case (a) $\phi=0.01$ shows a peak around $d/\eta\approx 10$, beyond which the distribution rapidly falls off due to the dispersion being dilute (see 4th panel in the top row of figure 4). Due to infrequent coalescence, large droplets are not formed very often. This was also reflected in the $N_{d}$ evolution (figure 6) which does not fluctuate as much as higher $\phi$ simulations. Cases (b), (c) and (d) with $\phi\geq 0.06$ follow a $d^{-10/3}$ power law in an intermediate droplet size range $15<d/\eta<30$, showing the formation of larger droplets. This is in accordance to the prediction of Garrett et al. (2000) and Deane & Stokes (2002) for droplets in the inertial range of turbulence, which was also found by Skartlien et al. (2013) in their simulations.

Also, for $\phi\geq 0.15$, a secondary peak appears at high $d/\eta$, which is due to a few large connected regions forming due to coalescence, which remain connected despite occasional satellite droplets breaking off. Hence in these simulations, droplets in an intermediate range are less frequent, as upon formation they would soon coalesce with the larger connected region. This is first a consequence of having a high volume fraction at a lower turbulence intensity. At higher $Re_{\lambda}$, the large region would be unstable and hence break apart forming droplets with a range of diameters. Secondly, the formation of this larger connected region also depends on $Ch$. If a simulation is performed on a much larger domain for the same volume fraction $\phi=0.20$ and turbulence intensity $Re_{\lambda}=45$, due to an increased separation between $d$ and $\zeta$ (lower $Ch$), coalescence would be inhibited. It should be noted that the uncertainty in determination of $d$ is around $10\%$, as shown in appendix A.

Further, $\eta\approx 1.5$ [lu] here and given that the interface width $\zeta\approx 5-6$ [lu], the $Ch$ for these droplets is approximately in the range $0.03<Ch<1.5$. Droplets towards the higher side where $Ch\sim O(1)$ will be unstable and prone to dissolution, which is reflected in the distributions falling off to the left of $d/\eta\approx 5$. In physical systems, small droplets are stable and can only be destroyed by coalescence. Resolving droplets in this range of diameters (where $d/\eta\sim O(1)$) will require over-resolving the Kolmogorov scale (to decrease the relative $Ch$), as was done by Komrakova et al. (2015a). Lastly, the length scales are ordered as $N_{x}>\mathcal{L}\gg d\gg\zeta>\eta$ for cases P1-P3 while $N_{x}>\mathcal{L}>d\gg\zeta>\eta$ for cases P4 and P5 (where due to higher $\phi$, the long droplet filaments can be of length $\sim\mathcal{L}$).

Multiphase kinetic energy spectra

Figure 8 shows the kinetic energy spectra for the droplet laden simulations, in comparison to the single-phase turbulence simulation with identical forcing. The first effect to note is the suppression of the inertial range (i.e. deviation from the $k^{-5/3}$ law) which is seen more clearly in the inset figure, and has been found previously (Perlekar et al., 2014). For increasing $\phi$, the spectra between $1<k/k_{f}<10$ shift away from the inertial range scaling and the single-phase spectrum, which shows that the cascading mechanism becomes weaker. Interestingly, the spectra pass through a single point, which is marked by the vertical line in the inset figure. This point is very close to the inverse of the Hinze length scale given by $d_{\mathrm{max}}=0.725(\rho/\gamma)^{-3/5}\epsilon^{-2/5}$. Beyond this point, the higher $\phi$ simulations contain higher energy at the smaller scales (large wavenumbers). This is due to coalescence which generates small scale eddies and is more frequent at higher $\phi$. Two or more droplets coalescing add kinetic energy to the flow by loss of surface energy due to a reduction in overall surface area. The crossover of the multiphase spectra (for $\phi\geq 0.15$ cases) with the single-phase spectrum shows that the dissipation range has higher energy in the presence of droplets, as was also reported by Perlekar et al. (2014). Interestingly, Ten Cate et al. (2004) also found such a spectral crossover at increasing volume fractions for solid spherical particles in turbulence.

The $\phi=0.01$ simulation has the lowest energy at high wavenumbers, as coalescence events are rare, and the droplet sizes are smaller which derive energy from eddies corresponding to higher wavenumbers. Lastly, a small jump in the spectra at $k/k_{f}\approx 50$ is consistently seen for all cases, which corresponds precisely with the interface width in our simulations (i.e $5-6$ [lu]). The extra energy there is due to the spurious currents present in the system, which are found to be much weaker that the physical velocity scales. Komrakova et al. (2015a) reported that spurious currents completely dominated the higher wavenumbers of the kinetic energy spectra in their turbulent dispersion simulations, due to which the spectra could not be well represented. Our work does not suffer from this problem, and although spurious currents are present, they do not adversely influence our results.

4.3 Effect of turbulence intensity

As mentioned earlier, the idea behind applying turbulence is to cause fragmentation of the dispersed phase, and the number of droplets thus formed depends upon the intensity of turbulence. We now keep the volume fraction fixed at $\phi=0.1$ and increase the turbulence intensity by increasing the forcing amplitude. These are simulations T1-T5 in table 1, and are run for $t=0.5$ million time steps each, though the simulations will have different $\tau_{k}$. Figure 9 shows the evolution of the normalized phase fraction over time, and as expected, at higher turbulence intensities (which leads to a higher $Re_{\lambda}$), $\phi/\phi_{0}$ reduces over time to an individual stable value. For the case T4, all the droplets dissolve within $10^{5}\ \Delta t$ , which shows that for this combination of parameters (refer to table 1), turbulence forcing undesirably outclasses the PP-LB phase segregation. The small droplets formed in this system are subsequently unstable (due to $Ch\sim O(1)$), which causes complete dissolution of the dispersed phase. Upon increasing the liquid-liquid repulsion parameter $G_{\alpha\beta}$ (hence also changing the fluid composition and dimensionless numbers that include interfacial tension, like the Weber or Ohnesorge number) in case T5, we see that for the same turbulence intensity as case T4, $\phi/\phi_{0}$ remains stable. This reaffirms that, with the original PP-LB method, certain regions of the turbulent emulsions parameter space can be simulated properly, while in other cases (case T4 and to some degree also case T3) simulations may require additional numerical remedies like the mass correction scheme of Biferale et al. (2011); Perlekar et al. (2012) or an enhanced Kolmogorov scale resolution (to achieve higher Cahn numbers) as done by Komrakova et al. (2015a).

Droplet number density evolution

Figure 10 shows the evolution of the number of droplets for cases T1-T5 for varying turbulence forcing amplitudes (excluding case T4 where all the droplets eventually dissolve due to a relatively weaker inter-component repulsion). Increasing $Re_{\lambda}$ increases the average number of droplets in the system (obtained by averaging $N_{d}$ between times $t=1\times 10^{5}$ and $t=5\times 10^{5}$) from around $\overline{N_{d}}=50$ for $Re_{\lambda}=44$ to $\overline{N_{d}}=600$ for $Re_{\lambda}=90$. Further, two interesting features in the evolution of $N_{d}$ can be noted. First is that the variation in $N_{d}$ increases with $Re_{\lambda}$, which results in a larger standard deviation of $\overline{N_{d}}$. This also makes it possible to generate a wider distribution of droplet diameters in the system. The second striking feature is the quasi-periodic rise and fall in the droplet number concentration (with a period of around $8-10\mathcal{T}$), most distinctly seen for the $Re_{\lambda}=90$ simulation (case T5). There seems to be an upper limit to the number of droplets that can be formed, which apart from constraints of resolution and maximum sphere-packing of the domain while keeping the diffuse interfaces apart, indicates also at the underlying physical mechanisms. At its peak, $N_{d}\approx 900$ here, a state corresponding to most droplets being rather small that cannot undergo additional breakup as they would all be well below the Hinze scale. These droplets are advected around by the flow, and they begin to coalesce when they collide, causing $N_{d}$ to drop to its lower limit, where a significant number of droplets will again be larger than the Hinze scale, and they begin to break and this cycle continues. We shall revisit this feature in detail in section 5.1.

Dispersion morphology

The dispersion morphology can be quantified with the concentration spectrum $k^{2}S(k,t)$, a quantity commonly used to describe coarsening dynamics (or spinodal decomposition) (Chen et al., 2000; Perlekar et al., 2014). Here $S(k,t)$ is the shell-averaged structure factor which is obtained using the Fourier transform $\hat{\phi}_{\mathbf{k}}$ of the density-density correlation function $\phi-\overline{\phi}$, where $\phi=(\rho_{\alpha}-\rho_{\beta})$ and $\overline{\phi}$ is the mean value of $\phi$. The quantity $\hat{\phi}_{\mathbf{k}}$ is shell-averaged in wavenumber space to obtain $S(k,t)$ as follows

[TABLE]

Here $\sum$ denotes summation over wavenumber shells $k\in[k-1/2,k+1/2]$ where $k=\sqrt{\mathbf{k}\cdot\mathbf{k}}$. Further, a characteristic length $L(t)$ can be calculated using the first moment of $S(k,t)$ as follows

[TABLE]

Figure 11 shows the concentration spectrum for cases T1-T5, which reveals the influence of the turbulence intensity on the dispersion morphology. As $Re_{\lambda}$ is increased, smaller droplets begin to dominate the system which is seen from the shift towards higher wavenumbers in $k^{2}S(k,t)$. This is also reflected in the time averaged characteristic length $L$ which decreases from $100$ to around $40$ [lu]. Note that cases T3 and T5 almost overlap. This shows that the turbulence intensity and the repulsion parameter $G_{\alpha\beta}$ (or interfacial tension) compete to create a particular morphology, and a similar droplet distribution may be attained by varying the two factors in tandem.

4.4 Effect of domain size

In simulations corresponding to D1-D4 in table 1, we successively increase the domain size $N_{x}$ while keeping the turbulent energy density the same. This essentially creates a separation between the domain size $N_{x}$ and the forcing scale $\mathcal{L}$, and allows for a better resolution of the largest droplet extension before breakup. So far, studies on turbulent dispersions have focused on maximizing the turbulence intensity which is reflected in the general proclivity for achieving higher $Re_{\lambda}$ in DNS simulations with Lagrangian objects like particles or droplets (Toschi & Bodenschatz, 2009). This finds implicit justification in that $Re_{\lambda}$ in real systems where droplets and turbulence interact is typically very high (for instance droplet-turbulence interaction in clouds occurs at $Re>10^{6}$ (Falkovich et al., 2002; Shaw, 2003), where $Re_{\lambda}=\sqrt{15Re}$ for homogeneous, isotropic turbulence). In periodic domain DNSs, a high $Re_{\lambda}$ is achieved by minimally resolving the Kolmogorov scale (the $k_{\mathrm{max}}\eta>1$ condition (Moin & Mahesh, 1998)), while forcing turbulence at the largest possible scales i.e. $\mathcal{L}\approx N_{x}$ or $k_{f}\approx 1-2k_{\mathrm{min}}$. This wide separation of scales manifests a high $Re_{\lambda}$. There are a few connected issues regarding the relative resolution of the various length scales, which is the focus of this section.

The first issue, emphasized by Komrakova et al. (2015a), is the utility of over-resolving the Kolmogorov scale ($\eta\approx 10$ as opposed to $1$ [lu]), which helped remedy the rapid dissolution of droplets in their simulations. The increased resolution of $\eta$ and $d$ can also be seen as a reduction in the size of the interface $\zeta$, i.e. an decrease in the Cahn number $Ch$, since the interface thickness (in terms of the number of lattice spacings) remains constant while smaller droplets and turbulent length scales become better resolved (i.e. they become larger relative to $\zeta$). Droplet dissolution also depends on the relative strengths of turbulence and phase segregation (effectively the interfacial tension), as was demonstrated in section 4.3.

The other issue is that weak large scale forcing introduces a caveat that droplets tend to deform into long, slender filaments that stay connected across the periodic domain. The length scale of the largest droplet extension before breakup $d^{\mathrm{ext}}$ can become comparable to $N_{x}$, which means that breakup cannot be resolved. The dispersion then forms a complex tangled structure, which does not morphologically resemble an emulsion. This issue is aggravated by high volume fractions of the dispersed phase.

In simulations D1-D4, we increase the forcing wavenumber $k_{f}$ by the same factor as the domain size $N_{x}$ (while keeping the forcing amplitude $A$ the same). The upper and lower wavenumber bounds ($k_{a},k_{b}$) are also suitably adjusted to distribute the forcing over a reasonable wavenumber range (and all integer values in the range $k\in\left[k_{a},k_{b}\right]$ are considered). This ensures that the energy density remains the same in these simulations, while larger droplet deformations ($d^{\mathrm{ext}}$) can be resolved accurately. Successively increasing the domain size in this way allows separating $N_{x}$ from $\mathcal{L}$. Note that doing this does not decrease $Ch$ for droplets, as that would entail scaling $\mathcal{L}$ proportionally with $N_{x}$ while weakening the forcing amplitude such that $Re_{\lambda}$ remains constant and $\eta$ is over-resolved (the approach of Komrakova et al. (2015a)). We do not additionally pursue this as droplet dissolution is not significant in most of the parameter range considered in this study.

Figure 12 shows the droplets in the system (volume rendered) at $400\tau_{k}$ for increasing domain sizes. It can be seen that the largest structures in the $128^{3}$ domain span a significant fraction of the domain, whereas for increasing domain sizes the typical large scale structure becomes better resolved in relation to the domain size. The volume averaged droplet number density for these simulations was found to be almost identical.

The domain size limitation becomes apparent when considering the droplet distribution, as shown in figure 13. For the case of $N_{x}=128^{3}$ (D1), the distribution is limited to a small region around the peak, clearly being cut off at a secondary peak emerging at higher $d/\eta$ due to a lack of resolution of larger structures. This case is under-resolved, the issue made acute with the small domain size, significant $\phi$ and moderate $Re_{\lambda}\approx 30$. We include this case to emphasize that the same issue might arise in simulations with higher $Re_{\lambda}$ and $N_{x}$ of high volume fraction dispersions. Upon increasing $N_{x}$, the distribution successively assumes a longer tail which closely follows the $d^{-10/3}$ scaling for droplets larger than the Hinze scale.

Figure 14(a) shows the concentration spectrum for cases D1-D4, which first reflects the proper scaling as the spectra coincide for $k/k_{f}\geq 1$. The importance of resolving the dominant length scales characterizing the dispersion morphology vis-à-vis the domain size $N_{x}$ becomes apparent. The smallest wavenumber (largest length scale) that can be represented depends on $N_{x}$ as $k_{\mathrm{min}}=2\pi/N_{x}$. For case D1, $k_{\mathrm{min}}$ is very close to the wavenumber corresponding to the peak in the concentration spectrum, i.e. the dominant wavenumber $k_{d}$ (or length scale $N_{x}/k_{d}$). If $k_{\mathrm{min}}\approx k_{d}$, two issues would tend to arise. First is that the dominant length scale of the emulsion morphology is comparable to the domain size making its dynamics under-resolved. Secondly, this structure will strongly interact with an image of itself due to periodicity of the domain, which is undesirable. For successively larger domains, the dominant length scale does not change (due to the same energy density across simulations). Further, the separation of $k_{\mathrm{min}}$ and $k_{d}$ is increased, which confirms that the largest structures ($\sim N_{x}/k_{d}$) are well resolved, while even larger structures (in the range of $k<k_{d}$) are formed but not sustained as the peak of $S(k)$ resides at $k_{d}$. The characteristic length evolution in figure 14(b) also shows that the morphology obtained for D1-D4 is similar, and that the typical length scale $L(t)\approx 80$ becomes better resolved in relation to the grid size upon increasing $N_{x}$.

4.5 Effect of forcing wavenumber

To highlight the consequences of forcing turbulence at the largest possible scale i.e. having $\mathcal{L}$ comparable to $N_{x}$ (hence maximizing $Re_{\lambda}$), we performed an additional simulation D5 with $k_{f}=1.5k_{\mathrm{min}}$ and $\phi=0.2$ to compare with D4 ($k_{f}=6k_{\mathrm{min}}$, $\phi=0.15$), while keeping the forcing amplitude the same, which results in $Re_{\lambda}=118$ for case D5 (while $Re_{\lambda}=30$ for D4). Figure 15 shows the typical morphology of the droplets (at a random time instance), where visibly the D4 case seems to have smaller, more spherical droplets, while D5 shows more elongated filaments. Despite the higher $Re_{\lambda}$, the dispersion does not comprise smaller droplets as droplet sizes depend on $\left\langle{\epsilon}\right\rangle$ which remains mostly unchanged. The presence of elongated filaments in D5 reflects the nature of the turbulence forcing. For a long cylindrical filament, a higher wavenumber forcing will generate more curvature variations. This would increase the possibility of filament breakup driven by Rayleigh-Plateau instabilities. A lower wavenumber forcing would generate weaker curvature differences in a long filament, and the timescale of breakup of these filaments might be comparable to the timescale of the large eddies, in which case the filaments will only break when the direction of the large scale shear changes.

We further quantify the differences by calculating the droplet distribution for D4 and D5 (which have slightly different $\phi$), while also comparing simulations D2 (with $k_{f}=3.0$ and $Re_{\lambda}=30$) and P3 ($k_{f}=2.0$ and $Re_{\lambda}=47$) which have the same $\phi$, shown in figure 16. Indeed, the D5 case deviates from the $d^{-10/3}$ distribution above the Hinze scale reflecting the infrequent breakup of the long filaments that would lead to droplets in this range of sizes. This deficit of droplets shows up in a secondary peak at high $d/\eta$, which corresponds to the fewer, larger structures being sustained instead. A similar difference is seen between cases D2 and P3, where the P3 case shows a small peak at high $d/\eta$, again attributed to a lower wavenumber forcing. The same behaviour is reflected in the concentration spectrum as well between the cases (not shown here), where there is a relative increase in concentration at low wavenumbers for cases D5 and P3, although the characteristic length remains similar.

It is worthwhile to summarize the results from the domain size comparison and to draw conclusions. At modest $Re_{\lambda}$ ($<120$ in this study), the turbulence forcing wavelength and domain size influence the morphology. Having $N_{x}>\mathcal{L}\gg d$ (as in case D4) ensures sufficient resolution of the droplet breakup dynamics. While having $N_{x}\approx\mathcal{L}\gg d$ (case D5) causes the formation of longer filaments of the droplet fluid. Spatially, this causes the formation of larger droplets $d/\eta>100$ at the cost of some intermediate droplets $20<d/\eta<100$, for $d/\eta$ above the Hinze scale.

5 Turbulent emulsion dynamics

5.1 A quasi-equilibrium (limit) cycle

Droplet number density plots such as figure 10 show oscillations of $N_{d}$ around a typical mean value which characterizes the dispersion morphology. So far, studies on droplets in turbulence refer to this state as a “steady state” where coalescence and breakup equilibrate. Since these oscillations can be significant (with its extreme values remaining bounded, similar to kinetic energy and dissipation), the dynamics should more accurately be called as a quasi-equilibrium (limit) cycle in the system state space comprising $\left\langle{E_{k}}\right\rangle$, $\left\langle{\omega^{2}}\right\rangle$, $\left\langle{E_{\gamma}}\right\rangle$ (i.e. the specific interfacial area multiplied with the interfacial tension $\gamma$) and $N_{d}$. Coalescence and breakup equilibrate in a statistical sense only, while the instantaneous dynamics is governed by temporal branches of alternating dominance of coalescence and breakup. Note that the term “limit cycle” is used loosely to illustrate the dynamics, since truly closed trajectories in phase space were not found, perhaps primarily due to intermittency and non-periodicity of the numerical solutions.

A dominant mediator of droplet breakup is intense enstrophy (or dissipation $\epsilon$). Since dissipation destroys turbulent kinetic energy, it is interesting to note that its interaction with the dispersed phase is associated with interfacial wrinkling, deformation and breakup - all mechanisms that increase the amount of surface energy in the system at the cost of kinetic energy. This excess energy, however, is still available in the flow field, and true destruction of it (i.e. into heat) must be mediated via kinetic energy dissipation, which occurs by the generation of smaller scales in the flow due to coalescence or damped oscillations of deformed droplet interfaces. A higher globally averaged $\left\langle{\omega^{2}}\right\rangle$ can be expected to increase the chance of droplet breakup (as it also reduces the effective Hinze scale), and vice-versa. Hence the trends seen in the $N_{d}$ evolution should reflect those in the evolution of $\left\langle{\omega^{2}}\right\rangle$, which in turn should follow the peaks and valleys of the kinetic energy $\left\langle{E_{k}}\right\rangle$ evolution.

This hypothesis is found to be true, and is shown in figure 17 as the evolution of $\left\langle{E_{k}}\right\rangle$, $\left\langle{\omega^{2}}\right\rangle$, $N_{d}$ and $\left\langle{E_{\gamma}}\right\rangle$ for case T5, where in the four panels each quantity has been separately highlighted in colour. Each variable has been further normalized by its time averaged value (between $50\tau_{k}$ and $1000\tau_{k}$) so that the evolution profiles become comparable. The vertical lines between the top three panels show two peaks of $\left\langle{E_{k}}\right\rangle$, which are reproduced in the $\left\langle{\omega^{2}}\right\rangle$ evolution and eventually in the $N_{d}$ evolution. A correlation can be calculated between the two signals as

[TABLE]

where $\delta t$ is a time lag and the overbar is a temporal average. This has been done for the different signal pairs and is shown in figure 18. Here $\left\langle{\omega^{2}}\right\rangle$ is found to correlate strongly with $\left\langle{E_{k}}\right\rangle$ with a time delay of $\sim 0.3\mathcal{T}$. $N_{d}$ shows a very strong correlation with $\left\langle{\omega^{2}}\right\rangle$ at a time delay of $\sim 0.6\mathcal{T}$. Consequently, a significant correlation between $N_{d}$ and $\left\langle{E_{k}}\right\rangle$ is found at $\sim 0.9\mathcal{T}$. The converse effect of droplets on turbulence can also be hinted at with this figure, where the valleys of the $N_{d}$ evolution invariably coincide with peaks in the $\left\langle{E_{k}}\right\rangle$ evolution. This shows that when the droplet number density reduces due to coalescence, the excess surface energy is released into the flow as kinetic energy, which has been expounded by Dodd & Ferrante (2016). Since turbulence in our simulations is constantly forced (as opposed to Dodd & Ferrante (2016) who simulate droplets in decaying turbulence)- the variation in $\left\langle{E_{k}}\right\rangle$ in our simulations comes from a more complex confluence of the power input as well as the droplet dynamics.

We also observed this time delayed dynamics of $\left\langle{E_{k}}\right\rangle$ and $N_{d}$ for cases with different parameters like turbulence forcing amplitude and interfacial tension, although for some cases the effect was less explicit. Particularly, for weaker $\gamma$ or lower $Re_{\lambda}$, the $N_{d}$ oscillations were not as extreme as for case T5 (where turbulence intensity and interfacial tension are both relatively stronger forces), although the $\left\langle{E_{k}}\right\rangle$ and $N_{d}$ correlation was found to be strong. Generally, the dynamics can be described as follows. First the large scale structures generate higher velocity gradients at the dissipation scale (which may be due to the energy cascade if such exists) with an initial time lag. This larger dissipation rate is felt by the droplets, which respond by breaking up with a further time delay, increasing the number of droplets in the system. This process (from peaks in $\left\langle{E_{k}}\right\rangle$ to peaks in $N_{d}$) was consistently found to take place with a delay of around $\sim 0.9\mathcal{T}$ across different cases, which is roughly the lifetime of the large eddies. This finding can be important for droplet dynamics models like population balance equations, where breakup kernels rely upon the instantaneous local value of $\epsilon$. If the temporal aspect to droplet populations is important, a relaxation time should separate cause and effect which is not done currently as seen in the various models reviewed by Sajjadi et al. (2013).

We also found that the surface energy $\left\langle{E_{\gamma}}\right\rangle$ peaks prior to $N_{d}$ (last panel in figure 18), which hints at the underlying breakup mechanism. Since generally daughter droplets together have a higher surface area than the parent droplet, $\left\langle{E_{\gamma}}\right\rangle$ attaining its maximum values before $N_{d}$ suggests that droplets before breakup must form a rather elongated fluid filament, which has larger area than the daughter droplets formed after breakup.

In summary, the turbulent emulsion dynamics can also be interpreted as a quasi-periodic evolution in a state space comprising $\left\langle{E_{k}}\right\rangle$, $\left\langle{\omega^{2}}\right\rangle$, $N_{d}$ and $\left\langle{E_{\gamma}}\right\rangle$. Essentially, there are two bounded extrema in the droplet number density at a given turbulent intensity for a certain set of fluid properties. These correspond to a state of low $N_{d}$ which is marked by fewer, relatively large droplets. When dissipation attains a subsequent peak, several of these droplets must be larger than the instantaneous Hinze scale - which leads to accelerated droplet breakup with takes the system to its other extremum - a state marked with high $N_{d}$. Most of the droplets in this state are stable and cannot undergo further breakup. As dissipation reduces, these droplets are advected around, and due to a higher chance of droplet-droplet collisions, coalescence dominates the next part of the state-space evolution. These two states also exhibit slightly different dispersion morphologies, as illustrated in figure 19. The fluctuations in $N_{d}$ are caused by these two phases, where breakup and coalescence alternate in their dominance. In the $E_{k}-E_{\gamma}$ phase space, this can be viewed as (a somewhat erratic) evolution within a bounded region of finite $E_{k}$ and $E_{\gamma}$. We do find signatures of this behaviour, although to more accurately describe the $E_{k}-E_{\gamma}$ phase space requires further work where the contribution from breakup and coalescence are separately accounted for and the surface area is better resolved by simulating larger droplets in weaker turbulence. It should be noted, though, that the dynamics we report would correspond to local dynamics in larger droplet laden systems like stirred vessels or in clouds. When considering these systems as a whole, the equilibrium properties may not fluctuate as much as reported here, as the local fluctuations in different regions of the system would cancel out.

5.2 Vorticity and interface alignment

Figure 20 shows snapshots of enstrophy from a vertical cross-section of the varying $\phi$ simulations (P1-P4), with the droplet contours shown in black. Strong vortical regions are often found in the vicinity of the droplet interface and in the droplet wakes. There is strong interplay between the interfacial dynamics and dissipation, as strong vortical regions align with the interface (Shao et al., 2018) and cause wrinkling, and high local dissipative events can lead to droplet breakup (Perlekar et al., 2012).

The angle between the vorticity vector and the interface normal can be quantified by using the distribution of the cosine of the angle between them. First, the density field $\rho_{\beta}$ is converted to a phase indicator field $\psi=(\rho_{\beta}-\rho_{\beta}^{\mathrm{out}})/(\rho_{\beta}^{\mathrm{in}}-\rho_{\beta}^{\mathrm{out}})$, such that $\psi=1$ in the droplet region, $\psi=0$ in the carrier fluid region, and $0<\psi<1$ at the interface. The typical phase indicator gradient then becomes $\boldsymbol{\nabla}\psi=1/\zeta$, and the cosine of the orientation angle is calculated where $\boldsymbol{\nabla}\psi>0.01\zeta$ (where $0.01$ ensures all the interfacial region is considered while ignoring the bulk regions where $\boldsymbol{\nabla}\psi=0$ by construction) as follows

[TABLE]

where $\boldsymbol{\nabla}\psi/|\boldsymbol{\nabla}\psi|$ gives the interface normal. Recently, Shao et al. (2018) showed using this measure that vorticity tends to align tangentially to droplet interfaces in turbulent flow. Here we extend their result in figure 21 which shows the joint probability distribution of the cosine of the orientation angle $\theta$ and the normalized vorticity vector $\boldsymbol{\omega}/\langle\omega^{2}\rangle^{1/2}$. Stronger vorticity is found to be more prone to align tangentially to the interface, which can be associated to a highly swirling motion in the orthogonal plane ($\cos(\theta)=0$), which causes droplet accretion and subsequent tangential alignment of vorticity with interfaces. Weaker vorticity is incapable of exerting this influence on droplets, and hence would exhibit a uniform random distribution with respect to the interfaces. This result holds for droplets in the inertial range, while sub-Kolmogorov droplets might spin in local shear of the deep dissipation range.

5.3 Effect of droplets on flow topology

Local flow topology is described in terms of the three invariants ($P$, $Q$ and $R$) of the velocity gradient tensor $A_{ij}=\partial u_{i}/\partial x_{j}$, which form the coefficients of its characteristic equation

[TABLE]

where $P=-A_{ii}$, $Q=-A_{ij}A_{ji}/2$ and $R=-A_{ij}A_{jk}A_{ki}/3$. For incompressible flow, $P=0$ (i.e. the sum of the eigenvalues). In the $P=0$ plane (or the $QR$-plane), turbulent flow of diverse kinds produces a teardrop-like profile for the joint probability distribution of $Q$ and $R$ with four distinct flow topologies that have been illustrated in figure 22 (adapted from Ooi et al. (1999)). The curve $D=27R^{2}/4+Q^{3}=0$ divides the region with three real eigenvalues of $A_{ij}$ (below, where $D<0$) from the region with one real and a pair of complex conjugate eigenvalues (above, where $D>0$). The most dominant flow features are stable focus stretching ‘SFS’ (i.e. vortex stretching) and unstable-node/saddle/saddle ‘UN/S/S’ i.e. bi-axial straining (Chacin & Cantwell, 2000). ‘UFC’ corresponds to unstable focus compression (or vortex compression) and ‘SN/S/S’ is stable-node/saddle/saddle (or axial straining).

The presence of droplets or particles which interact with the flow can modify the distribution of flow topologies, which is a modification of turbulence structure at a more local and fundamental level than for instance modifications to the kinetic energy spectrum. This has been well investigated for particle laden turbulence (Rouson & Eaton, 2001; Bijlard et al., 2010) and recently shown for elastic polymers in turbulence by Perlekar et al. (2010). The effect of the latter can be similar to droplets which themselves are elastic objects due to interfacial tension, with the additional complexity of breakup and coalescence. Recently, Shao et al. (2018) showed a mild suppression of bi-axial straining in droplet laden turbulence upon changing the Weber number.

How droplets modify flow topology has not fully been investigated so far. Here, we first show the influence of increasing dispersed phase volume fraction on the $QR$ profiles calculated using simulations P1-P5. Since $Re_{\lambda}$ for these cases varies (and is almost a factor 2 lower than the corresponding single-phase turbulence simulation, see table 1), the normalization factor $\langle Q_{w}\rangle=\langle\omega^{2}\rangle/4$ (Ooi et al., 1999) is calculated for each case separately. This allows us to focus on the modification of flow features alone, without comparing the magnitude of these extreme $QR$ events. Figure 23 shows the $QR$ field sampled over the entire multiphase velocity field. For case (b) $\phi=0.01$, the profile is narrower than for single-phase turbulence, case (a), although the overall shape is similar. This might be due to the $\phi=0.01$ dispersion being dilute, which makes coalescence infrequent. Overall, in this case, the flow field is similar to that in single-phase turbulence, and coalescence generated smaller scale features are rare. This seems likely, as at successively higher volume fractions, cases (c) through (f), the $QR$ profile is influenced more significantly and it tends to become more symmetric across the $R=0$ line. This follows from an increase in the axial straining part of the flow, along with an extension of the profile into the $D>0$ and $R>0$ region which shows a relative increase in vortex compression as opposed to vortex stretching ($D>0$ and $R<0$).

Modification of the $QR$ profile due to an increase in $\phi$ hints that it is a consequence of turbulence being constrained by the dispersed phase. To validate this claim, in figure 24 the $QR$ profiles are shown while being sampled inside and outside the droplet regions (marked as “d” for droplet-phase and “c” for continuous-phase). This has been done for simulations D4 and D5 (which have the highest resolution, and significantly different $Re_{\lambda}=30$ and $118$ respectively). The $QR$ profiles have been sampled at 5 time instances separated by $100\tau_{k}$. The difference between the flow topology in the droplet and continuous phase is striking, where within the droplet region $QR$ profile seems to almost have flipped across the $R=0$ axis. There is a significant increase in axial straining and vortex compression inside the droplets. This may be ascribed to the presence of interfaces surrounding droplets which behave like elastic surfaces. Vortices being stretched inside the droplets will try to elongate the droplet along the stretching axis, and this will be counteracted by interfacial tension which would instead tend to compress vortices. Since vortex compression contributes to energy dissipation (Tsinober, 2009), an enhancement of energy dissipation might be expected inside droplets from these results. Further investigation of this is left for future work. The continuous phase $QR$ profile remains mostly tear-drop like, with minor increase in axial straining and vortex compression.

6 Conclusions

We perform direct numerical simulations of emulsions under homogeneous, isotropic turbulence conditions performed by using the pseudopotential lattice-Boltzmann method. New findings on droplet size distributions, multiphase kinetic energy spectra, coupled kinetic energy and droplet number density dynamics, interface-dissipation interactions and modification of turbulence flow topology in emulsions are reported.

The process of dispersion formation is investigated for varying volume fractions of the dispersed phase and varying turbulence intensities for an emulsion with a density and viscosity ratio of 1. Using an appropriate set of parameters (such that the pseudopotential repulsive force between components dominates the local turbulence force), the effect of droplet dissolution is mitigated, an issue that was found limiting in previous work (Perlekar et al., 2012; Komrakova et al., 2015a). While further maintaining spurious currents to well below the physical velocity scales, the multiphase kinetic energy spectra were shown to exhibit signatures of breakup and coalescence at wavenumbers smaller and larger than the inverse Hinze scale respectively.

At small wavenumbers, energy is primarily extracted from the flow, where a higher dispersed phase volume fraction $\phi$ extracts more energy due to the profusion of interfaces. At large wavenumbers, for successively higher $\phi$, the energy content of the dissipation range increases due to more frequent coalescence which generates smaller scale motions. The droplet distribution is shown to follow the Deane & Stokes (2002) $d^{-10/3}$ scaling above the Hinze (1955) scale.

The importance of the relative resolution between the various length scales that govern turbulence droplet simulations is emphasized. We show that it is important to resolve $N_{x}>\mathcal{L}$ to correctly capture droplet deformation and breakup at relatively weaker turbulence intensities and high volume fractions, where otherwise the droplet fluid can form a complex tangle of elongated filaments as the maximum droplet deformation becomes unresolved. We also maintain that $\mathcal{L}\gg d\gg\eta$, such that the droplets interact mainly with the inertial range of turbulence.

In line with recent results (Shao et al., 2018), vorticity is shown to strongly align tangentially to droplet interfaces. This effect was shown to be stronger for higher vorticity magnitudes. The presence of dispersed phase is also shown to significantly alter the flow topology represented by the joint pdf of $QR$, i.e. the second and third invariants of the velocity gradient tensor, much more acutely than recognized (Shao et al., 2018). The well known tear-drop like profile becomes almost flipped across the $R=0$ axis when sampled inside the droplet in comparison to sampling in the carrier phase. A striking increase in axial straining and vortex compression is found in the droplets, which hints at an interplay of interfacial tension which tries to counteract any extensional vortical motions. This result hints that droplets might cause enhanced dissipation in their interior. The carrier fluid topology retains features of the well known tear-drop profile (Chacin & Cantwell, 2000) with only minor increase in axial straining and vortex compression.

Last but not the least, we show for the first time the dynamics of the quasi-equilibrium between coalescence and breakup under constant energy input to the system which leads to sustained turbulence over very long simulation times (around $100\mathcal{T}$). This state is often called a “steady state”, although the dynamics more closely resembles a limit-cycle in the state-space of kinetic energy $\left\langle{E_{k}}\right\rangle$, enstrophy $\left\langle{\omega^{2}}\right\rangle$, droplet number density $N_{d}$ and surface energy $\left\langle{E_{\gamma}}\right\rangle$. The extreme values of $\left\langle{E_{k}}\right\rangle$ manifest in the $\left\langle{\omega^{2}}\right\rangle$ evolution with a certain time delay, which then again show up in the $N_{d}$ evolution leading to a time-delayed dynamics. The dispersion oscillates between two morphologies, the journey between them being mediated by alternating bouts of dominant breakup and coalescence. Surface energy was found to peak prior to droplet breakup, reflecting the underlying breakup mechanism which involves the stretching of droplet fluid filaments, which have a higher surface area than the subsequently formed daughter droplets.

We believe that this time delayed dynamics will be found in localized regions of much larger droplet laden systems, where the overall system may not exhibit significant fluctuations in state-space variables, as the localized fluctuations would cancel each other. However, in smaller, finite systems (as prevalent in turbulence resolving droplet laden simulations (Elghobashi, 2019)), this can be an important consideration, as the “steady state” can have its own interesting dynamics. These considerations of delayed temporal dynamics may also be relevant to developing more realistic breakup and coalescence kernels which currently correlate state-space variables instantaneously (Sajjadi et al., 2013), which we have not explored given the limits of the current work.

Further investigation of the system evolution in the $\left\langle{E_{k}}\right\rangle-\left\langle{E_{\gamma}}\right\rangle$ phase space would help describe the exact exchange of energy, where the effects of coalescence and breakup would need to be isolated. This may be done by simulating larger droplets in weak turbulence, which would correspond to a detailed view on individual droplets near the dissipation range, and it is something we wish to investigate in the future.

We hope that this paper brings to attention the avenue of considering the details of resolved simulations from different perspectives (as we have attempted, while considering the limitations of our work). This helps reinforce our understanding of the phenomena at different levels. A statistical perspective (looking at spectra, time averaged quantities etc) helps with an overall description, while a dynamical systems perspective on the state-space helps pave the way for deciphering the true mediation of cause and effect like droplet-dissipation interactions and the modification of turbulence due to droplets, which we are only beginning to now understand.

Acknowledgements

SM would like to thank Jason Picardo (International Center for Theoretical Sciences, Bangalore), Luis Portela (TU Delft) and Jos Derksen (University of Aberdeen) for insightful discussions regarding this work. This work is funded by the Institute for Sustainable Process Technology (ISPT), the Netherlands.

Appendix A

In this section we briefly discuss the segmentation of droplets. The simulations output a continuous density field for both components $\alpha$ and $\beta$. As mentioned, the density variation of a component indicates the presence of droplets, where the density of component $\beta$ inside the droplet $\rho_{\beta}^{\mathrm{in}}\approx 4.4$ and that outside the droplet $\rho_{\beta}^{\mathrm{out}}\approx 0.4$ [lu] during the simulations. The droplet identification is done by picking a threshold density value $\rho^{\mathrm{c}}$, and every contiguous region with $\rho_{\beta}>\rho^{\mathrm{c}}$ is identified as a droplet, using a spatial segmentation (or clustering) algorithm previously developed by Siebesma & Jonker (2000). The results, hence, should be independent of $\rho^{\mathrm{c}}$.

Figure 25 shows the relative evolution of the volume fraction over time for the case $\phi=0.06$ for different threshold values shown as a fraction of $\rho^{\mathrm{in}}_{\beta}$, while $\rho^{\mathrm{out}}_{\beta}/\rho^{\mathrm{in}}_{\beta}=0.1$. This makes the useful range of thresholding values $\rho^{\mathrm{c}}/\rho^{\mathrm{in}}_{\beta}\in\left[0.1,1.0\right]$, where $\rho^{\mathrm{c}}/\rho^{\mathrm{in}}_{\beta}=0.55$ is halfway. Lower values of $\rho^{\mathrm{c}}$ create for slightly larger droplet regions, whereby $\phi/\phi_{0}$ increases. The real loss in $\phi$ is due to the dissolution of small droplets, and partially due to a loss in the droplet phase density to compensate for the increase in overall droplet surface area after breakup. This is affirmed by the increase in $\phi$ upon lowering $\rho^{\mathrm{c}}$ so that more of the interfacial region is considered to lie inside the droplets.

Figure 26 shows the evolution of the number of droplets $N_{d}$ in the system for different threshold magnitudes, which is seen to have minimal influence on $N_{d}$. Similarly, the droplet distribution was also found to be virtually unaffected by the choice of $\rho^{\mathrm{c}}$ as long as it lies within the droplet interface. The results reported in this work use $\rho^{\mathrm{c}}/\rho^{\mathrm{in}}_{\beta}=0.57$, which is very close to the halfway value.