On the interaction of Taylor length-scale size droplets and homogeneous shear turbulence

Abstract

The main objective of the present work is to explain the physical mechanisms occurring in droplet-laden homogeneous shear turbulence (HST) with a focus on the modulation of turbulence kinetic energy (TKE) caused by the droplets. To achieve such an objective, first, we performed direct numerical simulations (DNS) of HST laden with droplets of initial diameter approximately equal to twice the Taylor length scale of turbulence, droplet-to-fluid density and viscosity ratios equal to ten and a 5 % droplet volume fraction. We investigated the effects of shear number and Weber number on the modulation of TKE for $Sh \approx 2$ and $4$, and $0.02 \le {{We_{rms}}} \le 0.5$. Then, we derived the TKE equations for the two-fluid, carrier-fluid and droplet-fluid flow in HST and the relationship between the power of surface tension and the rate of change of total droplet surface area, providing the pathways of TKE for two-fluid incompressible HST. Our DNS results show that, for ${{We_{rms}}} = 0.02$, the rate of change of TKE is increased with respect to the single-phase cases, for ${{We_{rms}}} = 0.1$, the rate of change of TKE oscillates near the value for the single-phase cases and, for ${{We_{rms}}} = 0.5$, the rate of change of TKE is decreased with respect to the single-phase cases. Such modulation is explained from the analysis of production, dissipation and power of surface tension in the carrier-fluid and droplet-fluid flows. Finally, we explain the effects of droplets on the production and dissipation rate of TKE through the droplet ‘catching-up’ mechanism, and on the power of the surface tension through the droplet ‘shearing’ mechanism.

1. Introduction

The interaction of dispersed droplets and turbulence is important in many natural and industrial processes, e.g. rain formation (Shaw 2003), liquid–liquid emulsion (Berkman & Calabrese 1988), spray cooling (Qin et al. 2014) and spray atomization in combustors (Sirignano 1983; Faeth, Hsiang & Wu 1995). In these flows the droplet volume fraction is typically of the order of 1–10 % such that the turbulence is altered by droplet feedback on the surrounding fluid and by droplet–droplet interactions, placing the flow in the four-way coupling regime (Elghobashi 1994). A review on the state-of-the-art of direct numerical simulations (DNS) of turbulent flows laden with droplets or bubbles is provided by Elghobashi (2019).

The main objective of the present work is to explain the physical mechanisms occurring in droplet-laden homogeneous shear turbulence (DLHST) with a focus on the modulation of turbulence kinetic energy (TKE) caused by the droplets when compared with single-phase homogeneous shear turbulence (HST). Kida & Tanaka (1992) explained the physical mechanisms of TKE production in single-phase HST via DNS. Mashayek (1998) used DNS to study the modulation of HST at low Mach number by droplets of size smaller than the Kolmogorov length scale and found that the presence of non-evaporating droplets decreases the TKE of the carrier phase. Ahmed & Elghobashi (2000) explained the physical mechanisms responsible for the modulation of TKE budget in HST by sub-Kolmogorov solid particles via DNS and found that the presence of particles can decrease TKE production. Nicolai et al. (2014) conducted both DNS and experiments for the one-way coupling regime of HST laden with solid particles of the size of the Kolmogorov length scale and reported the preferential concentration and orientation of particle clusterings. Tanaka & Teramoto (2015) and Tanaka (2017) performed DNS of HST laden with finite-size particles of diameter ten times the Kolmogorov length scale ($D_0 \sim 10\eta$) and reported enhanced dissipation near the particle surface, in accordance with the findings of Lucci, Ferrante & Elghobashi (2010) for particle-laden decaying isotropic turbulence with particles from 16 to 35 Kolmogorov length scales. Kasbaoui, Koch & Desjardins (2019a) studied clustering of sub-Kolmogorov particles in HST via DNS and found three mechanisms leading to significant particle clustering. Kasbaoui (2019) performed DNS of particle-laden HST in the two-way coupling regime and found that the ratio of TKE production to dissipation increases or decreases with respect to that of the single-phase case depending on the particle mass loading.

In comparison to solid particles, droplets can deform, develop internal circulation, break up and coalesce with other droplets. Thus, with respect to the modulation of HST with solid particles, the interaction of finite-size droplets and HST is expected to reveal new physical mechanisms. For decaying isotropic turbulence laden with droplets of initial diameter of Taylor length-scale size, via DNS, Dodd & Ferrante (2016) explained the physical mechanisms of droplet/turbulence interaction and the pathways of TKE between droplets, carrier fluid and the interface between the two. Their results showed that the droplet-carrier-fluid interface represents a sink or source of TKE through the power of the surface tension due to the fluctuating velocity, ${{\varPsi '_{\sigma }}}$, which acts as a sink (source) of TKE when the total surface area of the interface increases (decreases). In decaying isotropic turbulence, the absence of mean shear translates to the absence of production of TKE. Thus, the next step of complexity in our understanding of droplet/turbulence interaction, including the effects of shear on droplets and the effects of droplets on the production of TKE, is studying DLHST.

The DNS of droplet-laden statistically stationary homogeneous shear turbulence (SS-HST) has been studied by Rosti et al. (2019). In our opinion, this work has three weaknesses, which we discuss herein. Firstly, in § 1 of their study the following question was posed as one of their three objectives: ‘How does the dispersed phase change the turbulent kinetic energy budget?’ In experiments, HST exhibits unbounded growth of length scales and of TKE (Tavoularis & Karnik 1989). Statistically stationary HST artificially constrains the growth of the large scales of the turbulent flow to the domain size, which produces ‘bursting’ events, i.e. sudden reductions of TKE. These sudden modulations of TKE are not due to the droplets, and affect the droplet dynamics as well. Thus, the bursting events have an effect on the rate of change of TKE, which may mask the effects of droplets on the TKE budget. Thus, while for single-phase flows, Sekimoto, Dong & Jiménez (2016) found similarities between SS-HST and the logarithmic layer in wall turbulence, we discourage its use for studying two-way or four-way coupling effects in particle-, droplet- or bubble-laden turbulent flows. This is analogous to the critique of studying two-way coupling effects for particle-laden forced isotropic turbulence, which forces the turbulence to a statistically stationary state, instead of decaying isotropic turbulence (Elghobashi & Truesdell 1993; Elghobashi 2019; Ferrante & Elghobashi 2022, p. 93). Secondly, Rosti et al. (2019) used the standard second-order Adams–Bashforth (AB$_2$) scheme to integrate the governing equations in time. This scheme is weakly unstable for simulations of HST performed with higher resolutions and longer simulation times (Schumann, Elghobashi & Gerz 1986; Kasbaoui et al. 2017). Although no instability was reported by Rosti et al. (2019), the AB$_2$ scheme can cause a spurious increase of the TKE energy spectrum at high wavenumbers, as shown herein in § 2.2.1. Finally, in their § 3.3, Rosti et al. (2019) included the relationship ${{\varPsi '_{\sigma }}} = (-\sigma /\mathcal {V}_m)\,\mathrm {d} A/\mathrm {d} t$ between the power of the surface tension due to the fluctuating velocity and the rate of change of the total droplet surface area. While such a relationship was derived by Dodd & Ferrante (2016) for isotropic turbulence, such an equation is not applicable to HST due to the presence of a mean velocity with shear. The equations for the power of surface tension for HST are derived from basic principles in Appendix C and reported and analysed in § 3.3.4.

In the present work we consider finite-size droplets larger than the Taylor length-scale size at the time of release ($D_0 \sim 2\lambda _0$, where $\lambda$ is the Taylor length scale and the subscript 0 means at droplet release time) in HST without gravity. We ensure that the simulation is physically meaningful by monitoring the expansion of the length scales, and we ensure that the two-point velocity autocorrelation in the $x$ direction diminishes to zero in less than half the length of the computational domain. We perform a parametric study of DNS of DLHST in which we vary the Weber number based on the root-mean-square (r.m.s.) velocity of turbulence and the shear number.

The paper is organized as follows. The mathematical description is presented in § 2, which includes the governing equations (§ 2.1) and the numerical method FastRK3P$^*$ (§ 2.2) that solves the issue of weak instability for simulating HST while being computationally efficient. Next, the results are presented and discussed in § 3, starting with a description of the initial conditions and the droplet parameters (§ 3.1). We introduce the TKE budget for droplet-laden flows in § 3.2, which is derived in Appendix B. We compare the time evolution of TKE in single-phase HST to that of DLHST, and explain the physical mechanisms of the droplet/turbulence interaction and the modulation of TKE due to the droplets in § 3.3. Finally, we summarize the findings of this work in § 4.

2. Mathematical description

2.1. Governing equations

The non-dimensional governing equations for an incompressible flow of two immiscible fluids with mean shear in the absence of gravity are

(2.1a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\dfrac{\partial \boldsymbol{u}}{\partial t} ={-} \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{u}\boldsymbol{u} \right ) - Sz \dfrac{\partial \boldsymbol{u}}{\partial x} - Sw \hat{\boldsymbol{i}} + \dfrac{1}{\rho} \left[ - \boldsymbol{\nabla} p + \dfrac{1}{{\textit{Re}}} \boldsymbol{\nabla}\boldsymbol{\cdot} \left ( 2 \mu \boldsymbol{\mathsf{S}}' \right ) + \dfrac{1}{{\textit{We}}}\,\boldsymbol{f}_{\boldsymbol{\sigma}} \right] , \end{gather}$$

where $\boldsymbol {u}=\boldsymbol {u}(\boldsymbol {x},t)$ is the fluid fluctuating velocity, $S = {\partial \bar {U}}/{\partial z}$ is the mean shear rate where $\bar {U}$ is the mean velocity, $p=p(\boldsymbol {x},t)$ is the pressure, $\rho =\rho (\boldsymbol {x},t)$ is the density, $\mu =\mu (\boldsymbol {x},t)$ is the dynamic viscosity, $\boldsymbol{\mathsf{S}}'=\boldsymbol{\mathsf{S}}'(\boldsymbol {x},t)$ is the strain-rate tensor of the fluctuating velocity ($\boldsymbol{\mathsf{S}}'=\tfrac {1}{2}[\boldsymbol {\nabla } \boldsymbol {u} + (\boldsymbol {\nabla } \boldsymbol {u})^{{\rm T}}]$). Here, ${\textit {Re}}$ and ${\textit {We}}$ are the Reynolds and Weber numbers, respectively, which are defined as

(2.2a,b)\begin{equation} {\textit{Re}}=\frac{\tilde{U} \tilde{L} \tilde{\rho}_c}{\tilde{\mu}_c}, \quad {\textit{We}}=\frac{\tilde{\rho}_c \tilde{U}^2 \tilde{L}}{\tilde{\sigma}}, \end{equation}

where $\tilde {U}$, $\tilde {L}$, $\tilde {\rho }_c$, $\tilde {\mu }_c$ and $\tilde {\sigma }$ denote, in order, the reference dimensional velocity, length, carrier-fluid density, carrier-fluid dynamic viscosity and surface tension coefficient used to non-dimensionalize the governing equations (2.1a) and (2.1b). The subscripts $c$ and $d$ indicate the carrier fluid and droplet fluid, respectively. Throughout the paper, all quantities are dimensionless unless they are accented with $\sim$. Also, note that ${\textit {Re}}=1/\nu _c$, where $\nu _c=\mu _c/\rho _c$ and ${\textit {We}}=1/\sigma$; thus, we may use ${\textit {Re}}^{-1}$ or ${\textit {We}}^{-1}$ instead of $\nu _c$ or $\sigma$ throughout the paper. We have chosen to non-dimensionalize the density and dynamic viscosity in (2.1b) by choosing the carrier fluid as the reference phase, such that $\rho _c = 1$ and $\mu _c = 1$. Here, $\boldsymbol {f}_\sigma =\boldsymbol {f}_\sigma (\boldsymbol {x},t)$ is the force per unit volume due to surface tension,

(2.3)\begin{equation} \boldsymbol{f}_\sigma = \kappa \delta(s)\boldsymbol{n}, \end{equation}

where $\kappa =\kappa (\boldsymbol {x},t)$ is the curvature of the droplet interface, $\boldsymbol {n}=\boldsymbol {n}(\boldsymbol {x},t)$ is the unit vector that is normal to the interface and directed towards the interior of the droplet, $\delta$ is the Dirac $\delta$ function that is needed to impose $\boldsymbol {f}_\sigma$ only at the interface position and $s$ is a normal coordinate centred at the interface, such that $s=0$ at the interface. Figure 1 of Dodd & Ferrante (2016) illustrates the direction of the interface normal $\boldsymbol {n}$ and the sign of the interface curvature $\kappa$.

2.2. Numerical method

In Dodd & Ferrante (2016) we employed a new pressure-correction method for simulating incompressible two-fluid flows called FastP$^*$ (Dodd & Ferrante 2014). This method reduces the variable coefficient Poisson equation that arises in solving the incompressible Navier–Stokes equations for two-fluid flows to a constant coefficient equation, which, depending on the boundary conditions, e.g. for periodic boundary conditions, can be solved with a fast Fourier transform (FFT)-based, fast Poisson solver rather than multigrid. FastP$^*$ uses the AB$_2$ scheme to integrate the governing equations in time. This scheme is known to be weakly unstable for simulating HST, particularly for higher resolutions and longer simulation times (Schumann et al. 1986; Kasbaoui et al. 2017). Kasbaoui et al. (2017) showed that this instability arises from using solutions from previous time steps in flux calculations. In order to solve this issue, we have developed a new numerical method for simulating DLHST called FastRK3P$^*$ that combines FastP$^*$ (Dodd & Ferrante 2014) with FastRK3 (Aithal & Ferrante 2020). FastRK3 is a third-order Runge–Kutta (RK3) pressure-correction method for solving the incompressible Navier–Stokes equations, which requires solving the Poisson equation of pressure only once per time step versus three times for a standard RK3 methodology (Aithal & Ferrante 2020). Also, Aithal, Tipirneni & Ferrante (2022) have shown that FastRK3 preserves the temporal accuracy of the underlying standard RK3 methodology even if the Poisson equation for pressure is solved only once per time step versus three for standard RK3. Thus, by combining these two methodologies, FastRK3P$^*$ has two main qualities: first, it does not use the solution from the previous time step to advance the solution in time, which is required by AB$_2$, and, second, it only requires one solution of the Poisson equation for pressure per time step. The first quality ensures that the issue of weak instability for simulating HST is solved, and the second makes the solver faster than the standard RK3 or Crank–Nicholson methods that require solving the Poisson equation multiple times per time step. FastRK3P$^*$ can be seen as the FastRK3 methodology extended to two-fluid immiscible flows, or as FastP$^*$ methodology using FastRK3 time integration instead of AB$_2$.

In § 2.2.1 we describe the FastRK3P$^*$ method that is used to solve numerically the two-fluid governing equations (2.1a) and (2.1b). This method is coupled to the volume-of-fluid (VoF) method presented in § 2.2.2, which is used to capture the motion of the droplet interface analogously to Dodd & Ferrante (2014).

2.2.1. FastRK3P$^*$

We solve the governing equations (2.1a) and (2.1b) throughout the whole computational domain, including the interior of the droplets. The domain is a rectangular prism with side lengths $(L_x,L_y,L_z)=(2\mathcal {L},\mathcal {L},\mathcal {L})$, where $\mathcal {L} = 1$. The governing equations are discretized in space in an Eulerian framework using the second-order central difference scheme on a uniform staggered mesh.

The solution algorithm begins by advecting the volume fraction of the droplet fluid, $C(\boldsymbol {x},t)$, based on the known velocity field $\boldsymbol {u}^n$. The volume fraction has value $C=0$ in the carrier fluid, $C=1$ in the droplet fluid and $0< C<1$ in cells containing the droplet interface. After computing $C^{n+1}$ (§ 2.2.2), the density and viscosity can be computed at time level $n+1$ as

(2.4) \begin{equation} \left.\begin{gathered} \rho^{n+1}(\boldsymbol{x}) = \rho_d C^{n+1}(\boldsymbol{x}) + \rho_c [ 1 - C^{n+1}(\boldsymbol{x})], \\ \mu^{n+1}(\boldsymbol{x}) = \mu_d C^{n+1}(\boldsymbol{x}) + \mu_c [ 1 - C^{n+1}(\boldsymbol{x})]. \end{gathered}\right\}\end{equation}

Runge–Kutta methods are a family of multi-step iterative methods that construct approximate velocities at intermediate time steps, starting with the velocity at time level $n$, to obtain the velocities at time level $n+1$. First, the computation of the approximate velocity omits the pressure term in (2.1b) and the second term on the right-hand side in (2.1b) is omitted. This term represents the advection of momentum by the mean velocity and is accounted for later in the solution algorithm by a ‘shear-remapping’ operation. The momentum operator for the right-hand side of (2.1b) with the omitted terms is defined as

(2.5)\begin{align} \boldsymbol{\mathcal{M}}(\boldsymbol{u}) ={-} \boldsymbol{\nabla}\boldsymbol{\cdot} \left ( \boldsymbol{u}\boldsymbol{u} \right ) - Sw \hat{\boldsymbol{i}} + \dfrac{1}{{\textit{Re}}} \left[ \dfrac{1}{\rho^{n+1}} \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( 2 \mu^{n+1} \boldsymbol{\mathsf{S}}' \right ) \right] + \dfrac{1}{{\textit{We}}} \left[ \frac{\kappa^{n+1} \boldsymbol{\nabla} C^{n+1} }{\bar{\rho}} \right], \end{align}

where the surface tension force, $\boldsymbol {f}_\sigma$, of (2.1b) has been substituted by using Brackbill's continuum surface force approach (Brackbill, Kothe & Zemach 1992),

(2.6)\begin{equation} \boldsymbol{f}_\sigma = \frac{\rho}{\bar{\rho}} \kappa \boldsymbol{\nabla} C, \end{equation}

where $\bar {\rho } \equiv \tfrac {1}{2}(\rho _1 + \rho _2)$. The interface curvature $\kappa ^{n+1}$ is computed using the height-function method (Cummins, Francois & Kothe 2005) with improvements developed by López et al. (2009).

The solution algorithm proceeds by calculating three intermediate velocities for the three stages of the RK3 algorithm using the FastRK3 method of Aithal & Ferrante (2020) as

(2.7)$$\begin{gather} \boldsymbol{u}_1^* = \boldsymbol{u}^n + \tfrac{1}{3}\Delta t \boldsymbol{\mathcal{M}}(\boldsymbol{u}^n), \end{gather}$$

(2.8)$$\begin{gather}\boldsymbol{u}_2^* = \boldsymbol{u}^n + \Delta t \left[ -\boldsymbol{\mathcal{M}}(\boldsymbol{u}^n) + 2 \boldsymbol{\mathcal{M}}\left (\boldsymbol{u}_1^* - \dfrac{1}{3}\Delta t \dfrac{\boldsymbol{\nabla} \phi}{\rho^{n+1}}\right ) \right], \end{gather}$$

(2.9)$$\begin{gather}\boldsymbol{u}_3^* = \boldsymbol{u}^n + \Delta t \left[\dfrac{3}{4}\boldsymbol{\mathcal{M}}\left (\boldsymbol{u}_1^* - \dfrac{1}{3}\Delta t \dfrac{\boldsymbol{\nabla} \phi}{\rho^{n+1}}\right ) + \dfrac{1}{4}\boldsymbol{\mathcal{M}}\left (\boldsymbol{u}_2^* - \dfrac{1}{3}\Delta t \dfrac{\boldsymbol{\nabla}\phi}{\rho^{n+1}}\right ) \right], \end{gather}$$

where the $\boldsymbol {\nabla } \phi$ terms represent a pressure-like field that correct $\boldsymbol {u}_1^*$ and $\boldsymbol {u}_2^*$ to be approximately divergence-free. For FastRK3P$^*$, these terms are defined as

(2.10)\begin{equation} \dfrac{\boldsymbol{\nabla}\phi}{\rho^{n+1}} = \left[ \dfrac{\boldsymbol{\nabla} p}{\rho} \right]^n . \end{equation}

The right-hand side of (2.10) is computed and stored at each time step according to the FastP$^*$ pressure splitting

(2.11)\begin{equation} \left[ \dfrac{\boldsymbol{\nabla} p}{\rho} \right]^n = \dfrac{1}{\rho_0} \boldsymbol{\nabla} p^n + \left ( \dfrac{1}{\rho^n} - \dfrac{1}{\rho_0} \right ) \boldsymbol{\nabla} p^*. \end{equation}

Next, the advection by the mean velocity is accounted for by the ‘shear-remapping’ operator that maps local values of velocity to values computed upstream according to the magnitude of the local mean velocity by using Fourier interpolation. The advection of mean velocity is, thus, applied to $\boldsymbol {u}_3^*$ with the ‘shear-remapping’ operator as

(2.12)\begin{equation} \boldsymbol{\check{u}}_3^* = \boldsymbol{u}_3^*(\boldsymbol{x} - \Delta t Sz \hat{\boldsymbol{i}}). \end{equation}

The pressure is computed by solving the Poisson equation (Dodd & Ferrante 2014)

(2.13)\begin{equation} \nabla^2 p^{n+1}=\boldsymbol{\nabla}\boldsymbol{\cdot} \left [ \left ( 1 - \frac{\rho_0}{\rho^{n+1}} \right ) \boldsymbol{\nabla} p^* \right ] + \frac{\rho_0}{\Delta t}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\check{u}}_3^*, \end{equation}

where we have split the pressure gradient term (Dong & Shen 2012) as

(2.14)\begin{equation} \frac{1}{\rho^{n+1}}\boldsymbol{\nabla} p^{n+1} \rightarrow \frac{1}{\rho_0}\boldsymbol{\nabla} p^{n+1} + \left ( \frac{1}{\rho^{n+1}} -\frac{1}{\rho_0} \right ) \boldsymbol{\nabla} p^*, \end{equation}

where $\rho _0=\min (\rho _1, \rho _2)$ and $p^*=2 p^n - p^{n-1}$. The advantage of using (2.14) is that it yields a constant coefficient Poisson equation (2.13) that can be solved efficiently using direct methods. Equation (2.13) is solved directly using a combination of a two-dimensional FFT in the $x$–$y$ plane and Gauss elimination in the $z$ direction (Schmidt, Schumann & Volkert 1984). Finally, we update the velocity field by applying the pressure correction to $\boldsymbol {\check {u}}_3^*$ as

(2.15)\begin{equation} \boldsymbol{u}^{n+1} = \boldsymbol{\check{u}}_3^* - \Delta t \left[ \dfrac{1}{\rho_0} \boldsymbol{\nabla} p^{n+1} + \left ( \dfrac{1}{\rho^{n+1}} - \dfrac{1}{\rho_0} \right ) \boldsymbol{\nabla} p^* \right]. \end{equation}

Figure 1 shows the difference in the TKE spectra when using AB$_2$ versus the FastRK3 method to simulate HST with shear number ($Sh = S/(u_{rms}/l)$) $Sh_0 \approx 2$, i.e. case A$_2$ (see table 1). The TKE spectrum from the AB$_2$ method shows unphysical fluctuations at higher wavenumbers, while the spectrum from the FastRK3 method decays as expected at high wavenumbers.

Figure 1. Spectrum of TKE, $E(\kappa )$, for single-phase HST at $tS = 2$ for ${\textit {Re}}_{\lambda 0} = 52$ and $Sh_0 \approx 2$, using the AB$_2$ (dotted line) and FastRK3 (solid line) time-integration schemes. The wavenumber, $\kappa$, is normalized by $\kappa _0 = 2{\rm \pi}/\mathcal {L}$.

Table 1. Flow parameters (dimensionless) at initial time ($t=0$), shear activation time ($t = 0.1$), droplet release time ($t_r = 0.5$ for case A$_2$, and $t_r = 0.3$ for case A$_4$) and at the final non-dimensional time ($t = 1.7$ for case A$_2$, and $t = 0.9$ for case A$_4$). Here, $t^*$ is defined in (3.4). Cases A$_2$ and A$_4$ are the single-phase HST flow with $Sh_0 \approx 2$ and $Sh_0 \approx 4$, respectively (see table 2).

2.2.2. Volume-of-fluid method

In the VoF method the sharp interface between the two immiscible fluids is determined using the VoF colour function, $C$, which represents the volume fraction of the droplet fluid in each computational cell. In our VoF method the interface between the two fluids is reconstructed using a piecewise linear interface calculation (Youngs 1982). The interface reconstruction in each computational cell consists of two steps: the computation of the interface normal $\boldsymbol {n} = (n_x,n_y,n_z)$ and the computation of the interface location. The algorithm that we use to evaluate the interface normal is a combination of the centred-columns method (Miller & Colella 2002) and Youngs’ method (Youngs 1982) known as the mixed-Youngs-centred method (Aulisa et al. 2007).

If we consider a characteristic function $\chi$ that has value $\chi =1$ in the droplet fluid and $\chi =0$ in the carrier fluid, $\chi$ is governed by the following advection equation:

(2.16)\begin{equation} \frac{\partial \chi}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \chi = 0. \end{equation}

The volume fraction $C_{i,j,k}$ of grid cell $i,j,k$ is related to the characteristic function $\chi$ by the integral relation

(2.17)\begin{equation} C_{i,j,k}(t) = \frac{1}{V_0} \iiint_{V_0} \chi(\boldsymbol{x},t) \,{\rm d}\kern0.7pt\boldsymbol{x}, \end{equation}

where $V_0$ is the volume of the $i,j,k$ cell. The volume fraction $C$ is advanced in time using the advection algorithm of Weymouth & Yue (2010), which is mass conserving, and wisps are redistributed and suppressed using the method of Baraldi, Dodd & Ferrante (2014).

2.2.3. Shear-periodic boundary conditions

In HST, periodic boundary conditions are applied in the streamwise $x$ direction and spanwise $y$ direction. In the $z$ direction, in which the mean carrier flow velocity varies ($\bar {U}(z)$, figure 2), the shear ($S = {\partial \bar {U}}/{\partial z}$) requires shear-periodic boundary conditions that, for a generic dependent variable $f$, are expressed as

(2.18)\begin{equation} f(x, y, L_z, t) = f(x - tS L_z, y, 0, t). \end{equation}

Figure 2. Schematic showing the shear-periodic boundary conditions in the $z$ direction.

Depending on the choice of $S$ and time step $\Delta t$, the $x$ position ($x - tSL_z$) on the right-hand side of (2.18) may fall in between grid points. The boundary values in the $z$ direction of velocity and pressure are computed using Fourier interpolation. The VoF variables, such as the interface normal, plane constant and curvature, are discontinuous and, thus, computing their boundary values via Fourier interpolation would be inaccurate. The way that we impose shear-periodic boundary conditions for the VoF variables is explained next. All VoF variables are located at cell centres along with the pressure field, while velocities are located at the staggered cell faces. FastRK3P$^*$ computes the momentum operator at staggered grid locations. In order to solve (2.1b) numerically, the surface tension term, $\boldsymbol {f}_\sigma$, must be computed on the staggered cell faces by averaging the values at the two nearest cell centres. In order to compute $\boldsymbol {f}_\sigma$ at the $z$ boundaries, the shear-periodic boundary conditions need to fill the values of the VoF variables in a number of ‘ghost cells’ in the $z$ direction next to the bottom and top boundaries in a two-step process. First, the VoF variables from a slab of four cells in the $z$ direction are copied from the interior, next to the bottom (and top) boundary, to the ghost cells next to the top (and bottom) boundary at the same $x,y$ locations. Next, the VoF advection algorithm is employed to shift the values of the ghost cells in the $x$ direction by the corresponding distance $\Delta x = S t L_z$ in accordance with (2.18). Next, for both the top and bottom $z$ boundaries, from the VoF variables in the four ghost cells, the interfaces are reconstructed and the curvature is computed, such that $\boldsymbol {f}_\sigma$ can be computed at the cell centres in the first ghost cells according to (2.3). Finally, $\boldsymbol {f}_\sigma$ is interpolated from the cell centres to the staggered cell faces at the $z$ boundaries.

3. Results and discussion

3.1. Initial conditions and droplet properties

3.1.1. Carrier flow parameters and initial conditions

The initial turbulent velocity field is generated by prescribing the TKE spectrum, $E(\kappa )$, and ensuring that the velocity field is isotropic, divergence free with respect to the discretized form of the continuity equation and that the velocity cross-correlation spectra, $R_{ij}(\kappa )$, satisfy the realizability constraint (Schumann 1977).

The initial energy spectrum at time $t=0$ is prescribed as (Pope 2000, § 6.5.3)

(3.1)\begin{equation} E(\kappa) = 1.5 [ \varepsilon_0^{2/3} \kappa^{{-}5/3} f_L(\kappa L) f_\eta(\kappa \eta) ], \end{equation}

where $\kappa$ is the wavenumber, $\varepsilon _0$ is the initial dissipation rate of TKE, $L \equiv k_0^{3/2}/\varepsilon _0$, where $k_0$ is the initial TKE, $f_L$ is given by

(3.2)\begin{equation} f_L(\kappa L) = \left ( \frac{\kappa L}{[(\kappa L)^2 + c_L]^{1/2}} \right ) ^{11/3}, \end{equation}

and $f_\eta$ is given by

(3.3) \begin{equation} f_\eta(\kappa \eta)=\exp \left\{{-}5.2 \left\{\left[(\kappa \eta)^4 + c_\eta^{4}\right]^{1/4} -c_\eta \right\} \right\}, \end{equation}

where $c_L=3.579$ and $c_\eta =0.440$. The constants $c_L$ and $c_\eta$ are calculated such that $E(\kappa )$ and $2 {\textit {Re}}^{-1} \kappa ^2 E(\kappa )$ integrate to $k_0$ and $\varepsilon _0$, respectively. The values of the dimensionless parameters at $t=0$ were $k_0=4.867 \times 10^{-2}$, $\varepsilon _0=1.243 \times 10^{-1}$ and ${\textit {Re}}=1.27 \times 10^{4}$. These parameters yield an initial Reynolds number based on the Taylor length scale of ${\textit {Re}}_{\lambda 0}=40$ (${\textit {Re}}_\lambda = \lambda (k^2/3)^{1/2}/\nu$). The non-dimensional time step used is $\Delta t = 0.1\Delta x/(SL_z)$.

The initial velocity field is allowed to develop with periodic boundary conditions and without shear (i.e. as decaying isotropic turbulence), until the skewness of the velocity derivative $S_u$ has reached $\approx -0.50$. At that time, a constant mean velocity gradient $S=5$ or $S=10$, which corresponds to an initial shear number $Sh_0 \approx 2$ or $Sh_0 \approx 4$, respectively, is imposed to the flow field. These values of $Sh$ are below the strong shearing regime ($Sh > 20$) that can be described using rapid distortion theory (Pearson 1959; Moffatt 1965; Kasbaoui, Koch & Desjardins 2019b). In order to ensure that our simulations are physically meaningful, we check that $\eta \kappa _{max} \ge 1$ at all times, where $\kappa _{max} = {\rm \pi}N$ is the maximum resolved wavenumber and $N = 600$ is the number of grid points in the $y$ and $z$ directions, while $N_x=2N$. Additionally, we check that the two-point Eulerian velocity autocorrelation in the $x$ direction diminishes to zero in less than half the length of $L_x = 2 \mathcal {L}$ at all times. To satisfy this condition, the domain length in the $x$ direction is double its length in the $y$ and $z$ directions.

Table 1 shows the dimensionless flow parameters at different times for the droplet-free flows (cases A$_2$ and A$_4$): $\ell$ and $\tau _\ell$ are the integral length and time scales, respectively; ${\textit {Re}}_\ell$ is the Reynolds number based on $\ell$; $\lambda$ is the Taylor length scale; $\eta$ and $\tau _\eta$ are the Kolmogorov length and time scales, respectively.

3.1.2. Droplet properties

We perform two simulations of single-phase flow, A$_2$ and A$_4$, and eight simulations of DLHST (table 2). Cases A$^{*}_{2}$ and A$^{*}_{4}$ are limiting cases in which the viscosity and density ratios are unity and the Weber number is infinity. We analyse the effects of varying the shear number $Sh = S/(u_{rms}/l)$ and the initial droplet Weber number based on the r.m.s. of velocity fluctuations ${{We_{rms}}} = D_0 u_{rms}^2 \rho _c/\sigma$, where $l$ is the integral length scale of turbulence and $D_0$ is the initial droplet diameter. In cases A$_2$–D$_2$, $Sh_0 \approx 2$ and in cases B$_2$–D$_2$, ${{We_{rms}}}$ increases from 0.02 to 0.5. In cases A$_4$–D$_4$, $Sh_0 \approx 4$ and in cases B$_4$–D$_4$, ${{We_{rms}}}$ increases from 0.02 to 0.5. These Weber numbers were selected, from a larger set of Weber numbers investigated, because they produced different effects on the evolution of TKE with respect to single-phase HST. The values of shear number were selected based on the simulations of Ahmed & Elghobashi (2000). The density and viscosity ratios for all droplet-laden cases are set to be $\varphi = 10$ and $\gamma = 10$, respectively. These properties were selected for their engineering relevance to spray combustion devices. For all cases, the initial number of droplets is $N_d = 1258$ and the initial droplet diameter is $D_0 = 0.0533$, for which the resulting droplet volume fraction and droplet mass fraction are, respectively, $\phi _v = 0.05$ and $\phi _m = 0.5$.

Table 2. Simulation properties (dimensionless) at droplet release.

The flow field evolves free of droplets until $tS = 2$, which corresponds to one flow through of the mean shear. To compare $Sh_0 \approx 2$ and $Sh_0 \approx 4$ cases, we introduce a new time quantity, defined as

(3.4)\begin{equation} t^* = t - t_{r}, \end{equation}

where $t_{r} = 0.5$ and $t_{r} = 0.3$ are the droplet release time for $Sh\approx 2$ and $Sh\approx 4$ cases, respectively. After droplets are released, all cases advance in time for three flow throughs of the mean shear, i.e. $0 \leq {t^*S} \leq 6$. Equal values of ${t^*S}$ between $Sh\approx 2$ and $Sh\approx 4$ cases correspond to equal shifts in the boundary conditions due to the mean shear, allowing for better comparison between different values of $Sh$. At ${t^*S} = 0$, droplets are randomly seeded in the domain under the constraint that the distance between droplet centres must be at least $2.1D_0$ and by setting the fluctuating velocity in the interior of the droplets to zero. Figure 3 shows that at ${t^*S} = 6$ the spectra of cases A$_2$ and A$_4$ are nearly identical to the spectra of cases A$^*_2$ and A$^*_4$, respectively, which indicates that setting the fluctuating velocity to zero in the droplet interior has a negligible effect on the spectra of HST. Wavelet-spectral analysis would be needed in order to accurately interpret the spectra of droplet-laden cases (Freund & Ferrante 2019). We also tested different initial droplet positions and found that for all droplet-laden cases, the values of $\mathrm {d} k /\mathrm {d} t$ match within $3\,\%$ for $5.25 < {t^*S} < 6$, and the values of $k$ match within $1.5\,\%$ at ${t^*S}=6$. Thus, we conclude that the results are nearly independent of the initial positions of the droplets.

Figure 3. Spectra of the TKE at ${t^*S} = 6$ in (a) $Sh_0 \approx 2$ cases and (b) $Sh_0 \approx 4$ cases. The wavenumber, $\kappa$, is normalized by $\kappa _0 = 2{\rm \pi}/ \mathcal {L}$.

3.2. Turbulence kinetic energy equations

In order to explain the fundamental physical mechanisms of the interactions of droplets with HST, we start by analysing the evolution equation of TKE, $k(t)$, for the two-fluid flow, $k_c(t)$ for the carrier-fluid flow and $k_d(t)$ for the droplet-fluid flow.

The evolution equation of $k(t)$ is derived in Appendix B as

(3.5)\begin{equation} \frac{\mathrm{d} k}{\mathrm{d} t} = {{\mathcal{P}}} -\varepsilon + {{\varPsi'_{\sigma}}}, \end{equation}

where

(3.6a)$$\begin{gather} k(t)\equiv \tfrac{1}{2} \langle \rho u_j u_j \rangle, \end{gather}$$

(3.6b)$$\begin{gather}{{\mathcal{P}}}(t) \equiv{-}S \langle \rho u w \rangle, \end{gather}$$

(3.6c)$$\begin{gather}\varepsilon (t) \equiv \frac{1}{{\textit{Re}}} \langle {\mathsf{T}}'_{ij}{\mathsf{S}}'_{ij} \rangle, \end{gather}$$

(3.6d)$$\begin{gather}{{\varPsi'_{\sigma}}} (t) \equiv \frac{1}{{\textit{We}}} \langle u_j f_{\sigma,j} \rangle, \end{gather}$$

where $\langle \cdots \rangle$ denotes instantaneous volume averaging over the entire computational domain. Here, $\boldsymbol{\mathsf{T}}'_{ij} = 2\mu {\mathsf {S}}'$ is the viscous stress tensor and ${\mathsf {S}}'_{ij}$ is the strain-rate tensor of the fluctuating velocity defined in § 2.1. In (3.5) and (3.6), ${{\mathcal {P}}}(t)$ is the production of $k(t)$, $\varepsilon (t)$ is the dissipation rate of $k(t)$ and ${{\varPsi '_{\sigma }}} (t)$ is the power of the surface tension due to the fluctuating velocity.

The evolution equation for the TKE of the carrier-fluid flow, $k_c(t)$, is

(3.7)\begin{equation} \frac{\mathrm{d} k_c}{\mathrm{d} t} = {{\mathcal{P}}}_c -\varepsilon_c + T_{\nu,c} + T_{p,c}, \end{equation}

and the evolution equation for the TKE of droplet-fluid flow, $k_d(t)$, is

(3.8)\begin{equation} \frac{\mathrm{d} k_d}{\mathrm{d} t} = {{\mathcal{P}}}_d -\varepsilon_d + T_{\nu,d} + T_{p,d}. \end{equation}

The terms in (3.7) and (3.8) are defined as

(3.9)\begin{equation} \left.\begin{gathered} k_{c}(t) \equiv \frac{1}{2} \langle \rho u_j u_j \rangle_c, \quad {{\mathcal{P}}}_c(t) \equiv{-}S \langle \rho uw \rangle_c, \quad \varepsilon_c(t) \equiv \frac{1}{{\textit{Re}}} \langle \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \rangle_c, \\ T_{\nu,c}(t) \equiv \frac{1}{{\textit{Re}}} \frac{\partial \langle \boldsymbol{\mathsf{T}}'_{ij} u_j \rangle_c }{\partial x_i}, \quad T_{p,c}(t) \equiv{-} \frac{\partial \langle u_j p \rangle_c }{\partial x_j} , \end{gathered}\right\} \end{equation}

and

(3.10)\begin{equation} \left.\begin{gathered} k_d(t) \equiv \frac{1}{2} \langle \rho u_j u_j \rangle_d, \quad {{\mathcal{P}}}_d(t) \equiv{-}S \langle \rho uw \rangle_d, \quad \varepsilon_d(t) \equiv \frac{1}{{\textit{Re}}} \langle \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \rangle_d, \\ T_{\nu,d}(t) \equiv \frac{1}{{\textit{Re}}} \frac{\partial \langle \boldsymbol{\mathsf{T}}'_{ij} u_j \rangle_d }{\partial x_i}, \quad T_{p,d}(t) \equiv{-} \frac{\partial \langle u_j p \rangle_d }{\partial x_j} , \end{gathered}\right\} \end{equation}

where $\langle \ldots \rangle _{c}$ and $\langle \ldots \rangle _d$ denote instantaneous volume averaging over the carrier fluid and droplet fluid, respectively. In (3.7)–(3.10), ${{\mathcal {P}}}_c$ and ${{\mathcal {P}}}_d$ are the productions of $k_c$ and $k_d$, $\varepsilon _c$ and $\varepsilon _d$ are the dissipation rates of $k_c$ and $k_d$, $T_{\nu,c}$ and $T_{\nu,d}$ are the viscous powers and $T_{p,c}$ and $T_{p,d}$ are the pressure powers. The power terms are related through the identity

(3.11)\begin{equation} {{\varPsi'_{\sigma}}} = (1- \phi_v) [ T_{\nu,c} + T_{p,c} ] + \phi_v [ T_{\nu,d} + T_{p,d} ], \end{equation}

which is also derived in Appendix B. We also analyse the modulation of the interfacial surface energy by the mean flow via the power of the surface tension due to the mean velocity, defined as

(3.12)\begin{equation} {{\bar{\varPsi}_\sigma}} (t) \equiv \frac{1}{{\textit{We}}} \langle \bar{U}_j f_{\sigma,j} \rangle, \end{equation}

which is discussed in more detail in § 3.3.4 and Appendix C.

The derived equations, (3.5), (3.7), (3.8) and (3.12), are summarized schematically in figure 4, which depicts the pathways for TKE exchange in DLHST, and, more generally, in two-fluid (liquid–liquid or gas–liquid) incompressible HST. All terms responsible for the evolution of $k$ (3.5), $k_c$ (3.7) and $k_d$ (3.8) are represented. The rectangles from left to right encompass the mean flow kinetic energy, the interfacial surface energy, the TKE of the two-fluid flow $k$ and the internal energy. In the current work the mean shear is prescribed and kept constant in time, which means that the mean flow kinetic energy is constant in time and that the modulation of the mean flow by the droplets is not allowed. This is indicated by the solid line boundary of the leftmost rectangle, as opposed to the dashed line boundaries of the other rectangles that represent energies that change in time. The light purple arrows represent the production, ${{\mathcal {P}}}$, of TKE in the carrier and droplet fluids due to the mean shear. The red arrows represent TKE of the carrier fluid and droplet fluid being transformed into internal energy by viscous dissipation, $\varepsilon$. The dark purple arrow represents mean flow kinetic energy being converted to interfacial surface energy by the power of the surface tension due to the mean velocity, ${{\bar {\varPsi }_\sigma }}$. The blue arrow represents TKE being exchanged for interfacial surface energy and vice versa by the power of the surface tension due to the fluctuating velocity, ${{\varPsi '_{\sigma }}}$. The power (or transport) terms $T_{\nu,c}$, $T_{p,c}$, $T_{\nu,d}$, $T_{p,d}$ (green arrows) act to redistribute TKE between the carrier fluid and droplet fluid or into interfacial surface energy via three bidirectional pathways: (i) carrier fluid $\leftrightarrow$ droplet fluid, (ii) carrier fluid $\leftrightarrow$ interface and (iii) droplet fluid $\leftrightarrow$ interface. This relationship is expressed mathematically by (3.11).

Figure 4. Schematic showing the pathways for TKE exchange in DLHST, or, in general, for two-fluid incompressible HST, summarizing the results of (3.5)–(3.12).

3.3. Comparison of TKE budget for single-phase and droplet-laden turbulence

In this section we present the effects of droplets on HST relative to the single-phase cases by analysing the terms of the TKE budget equation (3.5) and, then, we explain the underlying physical mechanisms.

3.3.1. Two-fluid TKE budget

Figure 5 shows the temporal evolution of $k(t)$ normalized by its initial value at droplet release time, $k/k_0$, for all cases. The average rates of change of TKE after ${t^*S} > 5$ are calculated and shown. For cases B$_2$ and B$_4$, the rate of change of TKE is increased with respect to the single-phase cases (A$_2$ and A$_4$). For cases C$_2$ and C$_4$, the rate of change of TKE oscillates near the value for the single-phase cases. For cases D$_2$ and D$_4$, the rate of change of TKE is decreased with respect to the single-phase cases. For all droplet-laden cases, $\textrm {d}(k/k_0)/\textrm {d}t$ is smaller for cases with larger values of ${{We_{rms}}}$.

Figure 5. Temporal evolution of the TKE, $k$, normalized by its initial value (a) $k_{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $k_{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

To explain why droplets modify the rate of change of $k$, we analyse the temporal evolution of the terms on the right-hand side of (3.5), which are ${{\mathcal {P}}}$, $\varepsilon$ and ${{\varPsi '_{\sigma }}}$. Figure 6 shows the temporal evolution of the production of TKE normalized by the initial dissipation rate of TKE, ${{\mathcal {P}}}/\varepsilon _0$. For cases B$_2$ and B$_4$, the production is increased with respect to the single-phase cases. For cases C$_2$ and C$_4$, the production closely matches that of the single-phase cases. For cases D$_2$ and D$_4$, the production is reduced with respect to the single-phase cases. For all droplet-laden cases, ${{\mathcal {P}}}$ is smaller for cases with larger values of ${{We_{rms}}}$.

Figure 6. Temporal evolution of the production of TKE, ${{\mathcal {P}}}$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

Figure 7 shows the temporal evolution of the normalized dissipation rate of TKE, $\varepsilon /\varepsilon _0$. For all droplet-laden cases, the dissipation rate is enhanced compared with the single-phase cases, with a larger increase in dissipation for cases with smaller values of ${{We_{rms}}}$.

Figure 7. Temporal evolution of the dissipation rate of TKE, $\varepsilon$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

Figure 8 shows the temporal evolution of the power of the surface tension due to the fluctuating velocity normalized by the initial dissipation rate of TKE, ${{\varPsi '_{\sigma }}} /\varepsilon _0$. For cases B$_2$ and B$_4$, ${{\varPsi '_{\sigma }}}$ oscillates around roughly $200\,\%$ of the initial dissipation rate, $\varepsilon _0$, which corresponds to $30\,\%$ of the instantaneous values of the dissipation rate, $\varepsilon$, at ${t^*S} = 6$. Therefore, in cases B$_2$ and B$_4$, ${{\varPsi '_{\sigma }}}$ represents a significant source of TKE for ${t^*S} > 3$. For cases C$_2$ and C$_4$, ${{\varPsi '_{\sigma }}}$ initially exhibits oscillations around zero up to $80\,\%$ of $\varepsilon _0$ (case C$_2$) and $200\,\%$ (case C$_4$), which decay to less than $30\,\%$ of $\varepsilon _0$ for ${t^*S} > 3$. For cases C$_2$ and C$_4$, ${{\varPsi '_{\sigma }}}$ represents a moderate source or sink of TKE for $0 < {t^*S} < 3$, and has a less significant role in the time evolution of the TKE for ${t^*S} > 3$. For cases D$_2$ and D$_4$, ${{\varPsi '_{\sigma }}}$ is limited to $\pm 20\,\%$ of $\varepsilon _0$, thus playing a less significant role in the time evolution of the TKE.

Figure 8. Temporal evolution of the power of the surface tension due to the fluctuating velocity, ${{\varPsi '_{\sigma }}}$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

3.3.2. Production of TKE

To explain why for cases B$_2$ and B$_4$, ${{\mathcal {P}}}$ is increased with respect to the single-phase cases, but for cases D$_2$ and D$_4$, ${{\mathcal {P}}}$ is reduced with respect to the single-phase cases, we analyse the contributions to ${{\mathcal {P}}}$ from the carrier-fluid production, ${{\mathcal {P}}}_c$, and the droplet-fluid production, ${{\mathcal {P}}}_d$, represented as

(3.13)\begin{equation} {{\mathcal{P}}} = (1 - \phi_v){{\mathcal{P}}}_c + \phi_v{{\mathcal{P}}}_d. \end{equation}

Figure 9 shows that, for droplet-laden cases, production is decreased in the carrier fluid compared with the single-phase cases, and figure 10 shows that, for droplet-laden cases, production is increased in the droplet fluid compared with the single-phase cases. The relative importance of these effects for the different cases is explained next.

Figure 9. Temporal evolution of the carrier-fluid contribution to the production of TKE, $(1-\phi _v){{\mathcal {P}}}_c$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

Figure 10. Temporal evolution of the droplet-fluid contribution to the production of TKE, $\phi _v{{\mathcal {P}}}_d$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

Figure 9 shows that ${{\mathcal {P}}}_c$ is smaller for all droplet-laden cases when compared with A$^*_2$ and A$^*_4$, and is smaller for the cases with larger ${{We_{rms}}}$. For single-phase HST, Kida & Tanaka (1992) explain how, on average, vortical structures are first elongated and then inclined by about $20^\circ$ to the streamwise direction by the mean shear. Pairs of inclined counter-rotating vortical structures cause a negative correlation of $uw$ in the region between them, and therefore, positive local production, ${{\mathcal {P}}}'= -S\rho uw$. The presence of the droplets interrupts this mechanism due to the droplets’ higher inertia with respect to the surrounding fluid, thereby reducing the regions of positive ${{\mathcal {P}}}'$ in the carrier fluid compared with cases A$^{*}_{2}$ and A$^{*}_{4}$. Figure 11 shows that the total droplet surface area, $A(t)$, decreases with decreasing ${{We_{rms}}}$, and that $A(t)$ is largest for cases D$_2$ and D$_4$. The droplets in cases D$_2$ and D$_4$ interrupt the carrier-fluid flow in the regions between pairs of counter-rotating vortical structures more than the droplets in cases B$_2$ and B$_4$ due to their larger total surface area. This explains why ${{\mathcal {P}}}_c$ is lowest for cases D$_2$ and D$_4$ among the cases studied.

Figure 11. Temporal evolution of the total surface area of the droplets, $A$, normalized by its initial value, $A_0$.

Figure 10 shows that ${{\mathcal {P}}}_d$ is larger for all droplet-laden cases when compared with A$^*_2$ and A$^*_4$, and is smaller for the cases with larger ${{We_{rms}}}$. Droplets with smaller ${{We_{rms}}}$, such as for cases B$_2$ and B$_4$, tend to deform less than droplets with larger ${{We_{rms}}}$, such as for cases D₂ and D₄. Due to the fact that the mean shear is positive, the mean velocity of a droplet whose centre of mass is higher in the $z$ direction tends on average to be larger than the mean velocity of a droplet whose centre of mass is at lower $z$. Because of this, droplets with lower ${{We_{rms}}}$ are more likely to keep their original spherical shape, catch up with droplets at lower $z$ on their path and collide such that their centres of mass are aligned along the northwest-southeast direction as depicted in figure 12. As pairs of droplets coalesce and return toward a spherical shape, the surface tension force squeezes the droplet fluid in the northwest-southeast direction, corresponding to a negative correlation of $uw$, and therefore, positive ${{\mathcal {P}}}' \equiv -S \rho u w$ in the droplet fluid as shown in figures 12 and 13. Figure 12 depicts a schematic of this droplet ‘catching-up’ mechanism. Figure 13 and supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.647 show the instantaneous results of two catching-up droplets obtained from a simulation placing two droplets in a shear flow with zero fluctuations and all droplet properties, numerical viscosity and mean shear matching those of case B$_4$. This effect only occurs when the Weber number is small enough to keep the shape of the droplets closer to their initial spherical shape. Figure 14 and supplementary movie 2 show that, for larger ${{We_{rms}}}$, droplets equivalent to those of case D$_4$ are deformed by the shear instead of returning toward a spherical shape. This explains why for cases B$_2$ and B$_4$ with smaller ${{We_{rms}}}$, there is a larger increase in ${{\mathcal {P}}}_d$ when compared with those of cases D$_2$ and D$_4$, respectively. Figure 15 and supplementary movie 3 show a contour plot of ${{\mathcal {P}}}'$ in the $x$–$z$ plane of case B$_4$. These figures demonstrate several instances of two droplets colliding in a similar fashion to the laminar two-droplet simulations. Figure 16 shows that more droplet collisions occur in cases B$_2$ and B$_4$ when compared with other cases, further showing that the droplet ‘catching-up’ mechanism increases ${{\mathcal {P}}}_d(t)$.

Figure 12. Schematic showing the droplet ‘catching-up’ mechanism.

Figure 13. Two droplets demonstrating the droplet ‘catching-up’ mechanism in laminar shear flow. All droplet properties, the numerical viscosity and the mean shear are equal to those in case B$_4$. Droplet interfaces are black lines, velocity vectors deviation from the mean velocity field are black arrows, colour contours of ${{\mathcal {P}}}'=-S\rho uw$ and temporal evolution of ${{\mathcal {P}}}_d = \langle {{\mathcal {P}}}' \rangle _d$ in insert. Results are shown for (a) ${t^*S}=0.5$; (b) ${t^*S}=0.9$; (c) ${t^*S}=1.4$; (d) ${t^*S}=1.9$.

Figure 14. Two droplets demonstrating the droplet ‘catching-up’ mechanism in laminar shear flow. All droplet properties, the numerical viscosity and the mean shear are equal to those in case D$_4$. Droplet interfaces are black lines, velocity vectors deviation from the mean velocity field are black arrows, colour contours of ${{\mathcal {P}}}'=-S\rho uw$ and temporal evolution of ${{\mathcal {P}}}_d = \langle {{\mathcal {P}}}' \rangle _d$ in insert. (a) ${t^*S}=0.0$; (b) ${t^*S}=2.5$; (c) ${t^*S}=5.0$; (d) ${t^*S}=7.5$.

Figure 15. Two instantaneous colour contours in the $x$–$z$ plane of ${{\mathcal {P}}}' = -S \rho uw$ and black lines for droplet interfaces highlighted within green circles where the droplet ‘catching-up’ mechanism is occurring for case B$_4$. (a) Case B$_4$, ${t^*S}=1.5$. (b) Case B$_4$, ${t^*S}=1.8$.

Figure 16. Temporal evolution of the total number of droplets.

For cases B$_2$ and B$_4$, the increase in ${{\mathcal {P}}}_d$ due to the droplet ‘catching-up’ mechanism is greater than the decrease in ${{\mathcal {P}}}_c$ when compared with the single-phase cases, which explains the overall increase in ${{\mathcal {P}}}$. For cases C$_2$ and C$_4$, both effects are relatively balanced, which explains why ${{\mathcal {P}}}$ closely matches the single-phase results. For cases D$_2$ and D$_4$, the decrease in ${{\mathcal {P}}}_c$ is the most significant of all droplet-laden cases and, additionally, the ‘catching-up’ mechanism does not cause a larger ${{\mathcal {P}}}_d$, which explains the overall decrease in ${{\mathcal {P}}}$. It should be noted that the VoF method used in the present work will always produce coalescence when the interfaces of two droplets come to occupy the same computational cell. Thus, the droplet ‘bouncing’ regime which may occur in droplet–droplet collisions is not captured by our VoF method, resulting in more coalescence events and no ‘bouncing’ regime.

3.3.3. Dissipation rate of TKE

To explain why $\varepsilon (t)$ is greater in all droplet-laden cases compared with single-phase cases, figure 17 shows the instantaneous two-dimensional contours of $\varepsilon '\equiv {\textit {Re}}^{-1} (\boldsymbol{\mathsf{T}}'_{ij} \boldsymbol{\mathsf{S}}'_{ij})$ of the computational domain at ${t^*S} = 3$. Figure 17 shows that $\varepsilon '$ is enhanced near the droplet interface for droplet-laden cases, and in the droplet interior in case B$_4$. This is explained by two separate mechanisms.

Figure 17. Instantaneous contours in the $x$–$z$ plane of $\varepsilon ' = Re^{-1} (\boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij})$ at ${t^*S} = 3$ for cases A$_4$, B$_4$ and D$_4$. (a) Case $\textrm {A}_4$, ${t^*S}=3.0$. (b) Case $\textrm {B}_4$, ${t^*S}=3.0$. (c) Case $\textrm {D}_4$, ${t^*S}=3.0$.

Firstly, the increased $\varepsilon '$ in the carrier phase near the droplet interface is due to the local increase of $\boldsymbol{\mathsf{S}}'_{ij}$ that is due to the local increase of the velocity gradient ($\partial u_i / \partial x_j$). Such increase in $\partial u_i / \partial x_j$ is caused by the droplet trajectories deviating from the motion of the carrier fluid because of the larger density of the droplets that, due to their higher inertia with respect to the carrier fluid, force the surrounding flow to move around them. Note that the droplet trajectories deviate from the trajectories of the turbulent eddies (both large and small scales of motion) because the droplet Stokes numbers, based on either the integral time scale or the Kolmogorov time scale, are both much larger than unity. This had been observed also in droplet-laden decaying isotropic turbulence by Dodd & Ferrante (2016). Figure 18 shows the contribution from the carrier fluid to $\varepsilon (t)$. For all droplet-laden cases, $\varepsilon _c$ is larger for cases with smaller ${{We_{rms}}}$, and significantly larger for case B$_2$. Figure 17 shows regions of large $\varepsilon '$ in the carrier fluid surrounding droplets that have coalesced via the ‘catching-up’ mechanism. The significant increase in $\varepsilon _c$ for case B$_2$ is explained by the greater amount of coalescence events, since case B$_2$ has the least amount of droplets at the end of the simulation (figure 16).

Figure 18. Temporal evolution of the carrier-fluid contribution to the dissipation rate of TKE, $(1-\phi _v)\varepsilon _c$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

Secondly, the increased $\varepsilon '$ in the droplet interior in cases B$_2$ and B$_4$ is due to the droplet ‘catching-up’ mechanism. Figure 7 shows that, in cases B$_2$ and B$_4$, $\varepsilon (t)$ is greatly enhanced compared with all other cases. To explain this increase in magnitude, figure 19 shows the contribution from the droplet fluid to $\varepsilon (t)$. For all droplet-laden cases, $\varepsilon _d$ is larger for cases with smaller ${{We_{rms}}}$, and significantly larger in cases B$_2$ and B$_4$. As explained in § 3.3.2, in cases B$_2$ and B$_4$, the droplet ‘catching-up’ mechanism causes droplets to coalesce and return towards a spherical shape. Figure 20 and supplementary movie 4 show a contour plot of $\varepsilon '$ in the $x$–$z$ plane of two droplets in a shear flow with initial zero velocity fluctuations in the flow field for which all the droplet properties and the shear match those of case B$_4$, with an inset of the temporal evolution of the quantity $\varepsilon _d$. These figures show that after coalescence, the elastic return towards a spherical shape creates large velocity gradients ($\partial u_i / \partial x_j)$ in the flow inside the coalesced droplets near the region of coalescence, and, therefore, enhanced $\varepsilon '$ in the droplet fluid.

Figure 19. Temporal evolution of the droplet-fluid contribution to the dissipation rate of TKE, $\phi _v\varepsilon _d$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

Figure 20. Two droplets demonstrating the droplet ‘catching-up’ mechanism in laminar shear flow. All droplet properties, the numerical viscosity and the mean shear are equal to those in case B$_4$. Results are shown for (a) ${t^*S}=0.4$; (b) ${t^*S}=0.8$; (c) ${t^*S}=1.2$.

Finally, the contribution to $\varepsilon (t)$ from the carrier-fluid dissipation, $\varepsilon _c$, and the droplet-fluid dissipation, $\varepsilon _d$, is represented as

(3.14)\begin{equation} \varepsilon = (1 - \phi_v)\varepsilon_c + \phi_v\varepsilon_d. \end{equation}

The local increase of $\varepsilon '$ inside the droplets and near the droplet interfaces increases $\varepsilon (t)$ because $\varepsilon (t) = \langle \varepsilon ' \rangle$.

3.3.4. Power of the surface tension due to the fluctuating velocity

In § 3.3.1 the results have shown that the power of the surface tension due to the fluctuating velocity, ${{\varPsi '_{\sigma }}} (t) = ({1}/{{\textit {We}}})\langle u_j f_{\sigma,j} \rangle$, acts as a sink or source of TKE initially for $0 \leq {t^*S} \leq 3$, and acts as a source of TKE at later times for $3 \leq {t^*S} \leq 6$. We now explain in more detail the behaviour of ${{\varPsi '_{\sigma }}}(t)$, and why it changes for varying ${{We_{rms}}}$. In order to do so, we introduce the power of the surface tension due to the total velocity,

(3.15)\begin{equation} {{\varPsi_{\sigma}}} (t) \equiv \frac{1}{{\textit{We}}} \langle U_j\, f_{\sigma,j} \rangle. \end{equation}

Since $U_j = \bar {U}_j + u_j$, ${{\varPsi _{\sigma }}} (t)$ can be decomposed into the power of the surface tension due to the mean velocity, ${{\bar {\varPsi }_\sigma }} (t)$, and the power of the surface tension due to the fluctuating velocity, ${{\varPsi '_{\sigma }}}(t)$, as

(3.16)\begin{equation} {{\varPsi_{\sigma}}}(t) = {{\bar{\varPsi}_\sigma}}(t) + {{\varPsi'_{\sigma}}}(t). \end{equation}

In Appendix C we derive the following relationship between ${{\varPsi _{\sigma }}}(t)$ and the rate of change of the total droplet surface area $\mathrm {d} A(t) / \mathrm {d} t$:

(3.17)\begin{equation} {{\varPsi_{\sigma}}}(t) ={-} \frac{1}{{\textit{We}}} \frac{1}{\mathcal{V}} \frac{\mathrm{d} A(t)}{\mathrm{d} t}. \end{equation}

Here, $\mathcal {V}$ is the volume of the domain ($\mathcal {V}=2$) and $A(t)$ is the total droplet surface area defined as

(3.18)\begin{equation} A(t) \equiv \sum_{n=1}^{N_d(t)} \iint_{\partial \mathcal{V}_{c}^{(n)}(t)} \,\mathrm{d} \mathcal{A} = \sum_{n=1}^{N_d(t)} A^{(n)}(t), \end{equation}

where $N_d(t)$ is the instantaneous number of droplets, $\partial \mathcal {V}_{c}^{(n)}(t)$ is the instantaneous control surface that bounds the $n$th droplet interface from the carrier-fluid side and $A^{(n)}(t)$ is the instantaneous surface area of the $n$th droplet. Here, ${{\varPsi _{\sigma }}}(t)$ is directly proportional to the rate of change of droplet surface area (with opposite sign) and the constant of proportionality is the non-dimensional surface tension coefficient ${\textit {We}}$. Physically, ${{\varPsi _{\sigma }}}(t)$ represents the transfer of interfacial surface energy to the total kinetic energy of the flow, i.e. mean kinetic energy plus TKE. To analyse the behaviour of ${{\varPsi '_{\sigma }}}$, we first use (3.17) to determine how the interfacial surface energy evolves, and then consider the contributions to the evolution of the interfacial surface energy from the mean and fluctuating velocities.

Figure 11 shows that in cases B$_2$, C$_2$, B$_4$ and C$_4$, $A(t)$ increases slightly after the droplets are released, and then decreases below the initial total droplet surface area, $A(t=0)$, with small oscillations. These changes in total droplet surface area, $A(t)$, correspond to an overall transfer of energy from the interfacial surface energy to the total kinetic energy of the flow, with some transfer in the opposite direction, from the kinetic energy of the flow to the interfacial surface energy, due to the oscillations for $0 \leq {t^*S} \leq 3$. For cases B$_2$, C$_2$, B$_4$ and C$_4$, the total droplet surface area decreases in time. The only mechanism that can cause $A(t) < A(0)$ in the present flow is the prevalence of droplet coalescence over break ups and large deformations. Figure 16 shows that the total number of droplets decreases in time. For cases D$_2$ and D$_4$, instead $A(t)$ increases to $5\,\%$ of $A(0)$ and then remains larger than $A(0)$ for all times. These changes in total droplet surface area correspond to an overall transfer of energy from the total kinetic energy to the interfacial surface energy.

Here, ${{\bar {\varPsi }_\sigma }}(t)$ represents the transfer of energy from the interfacial surface energy to mean flow kinetic energy. Figure 21 shows that ${{\bar {\varPsi }_\sigma }}(t)$ is negative for all cases, with increasing magnitude as ${{We_{rms}}}$ decreases. Negative ${{\bar {\varPsi }_\sigma }}(t)$ indicates that the mean flow kinetic energy is being transferred to the interfacial surface energy. Physically, the negative contribution of ${{\bar {\varPsi }_\sigma }}(t)$ to ${{\varPsi _{\sigma }}}(t)$ means that the mean flow is acting to deform the droplets; thus, it is a source of interfacial surface energy.

Figure 21. Temporal evolution of the power of the surface tension due to the mean velocity, ${{\bar {\varPsi }_\sigma }}$, normalized by the initial value of the dissipation rate (a) $\varepsilon _{0_{Sh\approx 2}}$ for $Sh\approx 2$ cases, and (b) $\varepsilon _{0_{Sh\approx 4}}$ for $Sh\approx 4$ cases.

In all cases except D$_2$ and D$_4$, after energy is transferred from the mean flow kinetic energy to the interfacial surface energy by deforming the droplets, the surface tension force is large enough to bring the droplets back to a more spherical shape. When the droplets return to a more spherical shape, $A(t)$ decreases and the interfacial surface energy is transferred to TKE via ${{\varPsi '_{\sigma }}}(t)$. This explains why ${{\varPsi '_{\sigma }}}(t)$ increases for decreasing ${{We_{rms}}}$ (increasing surface tension force). When droplets with a larger surface tension force coalesce, they tend to come back to a more spherical shape compared with droplets with a smaller surface tension force. Figure 22 depicts what we have named the droplet ‘shearing’ mechanism, which summarizes the above described roles of the three powers of surface tension of (3.16).

Figure 22. Schematic showing the droplet ‘shearing’ mechanism.

4. Conclusions

We have performed DNS of HST laden with deformable droplets, whose diameter is approximately equal to twice the Taylor length scale of turbulence at the time the droplets are released in the flow field. The goal of this study was to extend the work by Dodd & Ferrante (2016) on droplet-laden decaying isotropic turbulence, to include and explain the role of the mean shear on the physical mechanisms of droplet/turbulence interaction for moderate initial shear numbers of $Sh_0 \approx 2$ and $Sh_0 \approx 4$. Understanding these mechanisms is a prerequisite for developing predictive, physics-based turbulence models (Ferrante 2022), as, for example, has been done by Freund & Ferrante (2021) with their mixed artificial neural network approach based on the work by Dodd & Ferrante (2016) and Freund & Ferrante (2019).

In order to perform the DNS study of DLHST, first we developed FastRK3P$^*$ (§ 2.2.1) by combining FastP$^*$ (Dodd & Ferrante 2014) with FastRK3 (Aithal & Ferrante 2020). FastRK3P$^*$ has two main qualities: first, it does not use the solution from the previous time step to advance the solution in time, which is required by AB$_2$, and, second, it only requires one solution of the Poisson equation for pressure per time step. These qualities make FastRK3P$^*$ computationally more efficient than projection methods using standard RK3 or Crank–Nicholson schemes for time integration, or using multigrid for solving the variable-density Poisson equation for pressure, while solving the issue of weak instability inherent to AB$_2$ for time integration.

Then, we performed DNS of DLHST using FastRK3P$^*$ (§ 3) for five cases at $Sh_0 \approx 2$, and for five cases at $Sh_0 \approx 4$. For the droplet-laden cases, we released 1258 droplets in HST at an initial Reynolds number based on the Taylor length scale of $Re_\lambda =52$ for $Sh_0 \approx 2$ cases, and $Re_\lambda =53$ for $Sh_0 \approx 4$ cases (table 1). For the droplet-laden cases, we varied the Weber number ($0.02 \le {{We_{rms}}} \le 0.5$), while the volume fraction was set to 5 %, and the droplet-to-fluid density and viscosity ratios were set to $\rho _d/\rho _c = 10$ and $\mu _d/\mu _c = 10$, respectively (table 2). The governing equations were discretized and solved in time, using FastRK3P$^*$ coupled with the VoF method to capture the droplet interface (Weymouth & Yue 2010; Baraldi et al. 2014), in a domain using $1200 \times 600 \times 600$ grid points, and each droplet was initially resolved by 32 grid points across its diameter. The new findings of this study are summarized below for the modulation of the TKE budget (§ 4.1) and the physical mechanisms that explain such modulation (§ 4.2).

4.1. Modulation of the TKE budget

In order to explain the modulation of TKE in DLHST, we derived the TKE budget equations, (3.5), (3.7) and (3.8), of DLHST for the two fluids, the carrier fluid and the droplet fluid (Appendix B). Compared with the TKE equations for decaying isotropic turbulence derived by Dodd & Ferrante (2016), the TKE equations for DLHST each have an additional term of production for $k$, $k_c$ and $k_d$, called ${{\mathcal {P}}}$, ${{\mathcal {P}}}_c$ and ${{\mathcal {P}}}_d$, respectively, which are due to the presence of the mean shear $S$. Additionally, we derived the equations relating the rate of change of total droplet interfacial area with the power of surface tension, (3.16) and (3.17). For decaying isotropic turbulence, due to the absence of a mean velocity the power of the surface tension is due only to the fluctuating velocity, ${{\varPsi '_{\sigma }}}$, whereas for DLHST, the power of the surface tension has the additional contribution due to the mean velocity, ${{\bar {\varPsi }_\sigma }}$, which also modulates the interfacial surface energy. These equations allowed us to summarize the pathways of TKE exchange in DLHST and, in general, for two-fluid incompressible turbulent shear flows in figure 4. Our main findings on the modulation of TKE budget terms in DLHST are summarized next.

1. (a) For ${{We_{rms}}} = 0.02$, the rate of change of TKE is increased with respect to the single-phase cases. For ${{We_{rms}}} = 0.1$, the rate of change of TKE oscillates near the value for the single-phase cases. For ${{We_{rms}}} = 0.5$, the rate of change of TKE is decreased with respect to the single-phase cases. For the droplet-laden cases, $\textrm {d}(k/k_0)/\textrm {d}t$ is smaller for cases with larger values of ${{We_{rms}}}$ (figure 5). Table 3 summarizes the described results, i.e. the effects of droplets of different ${{We_{rms}}}$ on $\mathrm {d} k/\mathrm {d} t$, ${{\mathcal {P}}}$, ${{\mathcal {P}}}_c$, ${{\mathcal {P}}}_d$, $\varepsilon$, $\varepsilon _c$, $\varepsilon _d$ and ${{\varPsi '_{\sigma }}}$ when compared with the single-phase cases. Red up arrows indicate an increase of the listed quantity when compared with the single-phase case, black tildes indicate a similar value to the single-phase case, and blue down arrows indicate a decrease in the quantity when compared with the single-phase case. Two arrows indicates an increase in magnitude of more than $100\,\%$, when compared with the case with the next largest magnitude at ${t^*S} = 6$.

2. (b) For TKE-increasing droplets (${{We_{rms}}} = 0.02$), the budget terms that contribute the most to the increase of $\mathrm {d} k/\mathrm {d} t$ are ${{\mathcal {P}}}_d$ and ${{\varPsi '_{\sigma }}}$. The significant increase of ${{\mathcal {P}}}_d$ results in an increase of ${{\mathcal {P}}}$. For TKE-neutral droplets (${{We_{rms}}} = 0.1$), the increase of ${{\mathcal {P}}}_d$ and the moderate contribution from ${{\varPsi '_{\sigma }}}$ are balanced by the increase of $\varepsilon$. For TKE-reducing droplets (${{We_{rms}}} = 0.5$), the increase of ${{\mathcal {P}}}_d$ is less than for other droplet-laden cases, leading to an overall decrease of ${{\mathcal {P}}}$ that results in a decrease of $\mathrm {d} k/\mathrm {d} t$. These results have shown that the flow inside the droplets in response to their dynamics has a dominant role in the modulation of the budget terms of TKE and, thus, of TKE, e.g. figures 9 and 10 for ${{\mathcal {P}}}_c$ vs ${{\mathcal {P}}}_d$ and figures 18 and 19 for $\varepsilon _c$ vs $\varepsilon _d$.

Table 3. Summary of ${{We_{rms}}}$ effects on $\mathrm {d} k/\mathrm {d} t$, ${{\mathcal {P}}}$, ${{\mathcal {P}}}_c$, ${{\mathcal {P}}}_d$, $\varepsilon$, $\varepsilon _c$, $\varepsilon _d$ and ${{\varPsi '_{\sigma }}}$ compared with the single-phase cases.

4.2. Droplet ‘catching-up’ and ‘shearing’ mechanisms

The main findings of this work on the physical mechanisms that explain the modulation of the TKE budget in DLHST are summarized next.

1. (a) The droplet ‘catching-up’ mechanism explains how droplets with lower ${{We_{rms}}}$ tend to coalesce and return toward a spherical shape. Droplets with lower ${{We_{rms}}}$ tend to deform less and, thus, stay more spherical. Due to the mean shear, spherically shaped droplets are more likely to catch up and collide with droplets in their path, so that their centres are aligned along the northwest-southeast direction (figure 12).

   1. (i) For ${{We_{rms}}} = 0.02$, as pairs of droplets coalesce and return toward a spherical shape, the surface tension force squeezes the droplet fluid in the northwest-southeast direction, generating a region of flow with negative $uw$ inside the newly formed droplet, and, therefore, positive ${{\mathcal {P}}}' \equiv -S \rho u w$ within the droplet fluid (figure 10).

   2. (ii) For ${{We_{rms}}} = 0.02$, the droplet ‘catching-up’ mechanism creates large velocity gradients ($\partial u_i / \partial x_j)$ in the flow inside the coalesced droplets, and, therefore, enhanced $\varepsilon '$ in the droplet fluid (figure 19).

2. (b) The droplet ‘shearing’ mechanism explains how the mean shear transfers kinetic energy from the mean flow to interfacial surface energy of the droplets, and, then, how the interfacial surface energy is transferred to TKE via ${{\varPsi '_{\sigma }}}$ (figure 22). The mean shear deforms the droplets, and so the interfacial surface energy is increased (${{\bar {\varPsi }_\sigma }} < 0$). For ${{We_{rms}}} = 0.02$, droplets coalesce, e.g. through the ‘catching-up’ mechanism, and then tend to return towards a spherical shape, such that the total droplet surface area is reduced. In that process, interfacial surface energy is transferred to TKE via ${{\varPsi '_{\sigma }}}$.