1. Introduction
The interaction of dispersed droplets and turbulence is important in many natural and industrial processes, e.g. rain formation (Shaw Reference Shaw2003), liquid–liquid emulsion (Berkman & Calabrese Reference Berkman and Calabrese1988), spray cooling (Qin et al. Reference Qin, Loth, Li, Simon and Van de Ven2014) and spray atomization in combustors (Sirignano Reference Sirignano1983; Faeth, Hsiang & Wu Reference Faeth, Hsiang and Wu1995). 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 Reference Elghobashi1994). 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 (Reference Elghobashi2019).
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 (Reference Kida and Tanaka1992) explained the physical mechanisms of TKE production in single-phase HST via DNS. Mashayek (Reference Mashayek1998) 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 (Reference Ahmed and Elghobashi2000) 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. (Reference Nicolai, Jacob, Gualtieri and Piva2014) 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 (Reference Tanaka and Teramoto2015) and Tanaka (Reference Tanaka2017) 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 (Reference Lucci, Ferrante and Elghobashi2010) for particle-laden decaying isotropic turbulence with particles from 16 to 35 Kolmogorov length scales. Kasbaoui, Koch & Desjardins (Reference Kasbaoui, Koch and Desjardins2019a) studied clustering of sub-Kolmogorov particles in HST via DNS and found three mechanisms leading to significant particle clustering. Kasbaoui (Reference Kasbaoui2019) 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 (Reference Dodd and Ferrante2016) 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. (Reference Rosti, Ge, Jain, Dodd and Brandt2019). 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 Reference Tavoularis and Karnik1989). 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 (Reference Sekimoto, Dong and Jiménez2016) 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 Reference Elghobashi and Truesdell1993; Elghobashi Reference Elghobashi2019; Ferrante & Elghobashi Reference Ferrante and Elghobashi2022, p. 93). Secondly, Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019) 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 Reference Schumann, Elghobashi and Gerz1986; Kasbaoui et al. Reference Kasbaoui, Patel, Koch and Desjardins2017). Although no instability was reported by Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019), 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. (Reference Rosti, Ge, Jain, Dodd and Brandt2019) 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 (Reference Dodd and Ferrante2016) 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
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
$$\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
\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,
\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 (Reference Dodd and Ferrante2016) illustrates the direction of the interface normal 
$\boldsymbol {n}$ and the sign of the interface curvature 
$\kappa$.
2.2. Numerical method
In Dodd & Ferrante (Reference Dodd and Ferrante2016) we employed a new pressure-correction method for simulating incompressible two-fluid flows called FastP
$^*$ (Dodd & Ferrante Reference Dodd and Ferrante2014). 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. Reference Schumann, Elghobashi and Gerz1986; Kasbaoui et al. Reference Kasbaoui, Patel, Koch and Desjardins2017). Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017) 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 Reference Dodd and Ferrante2014) with FastRK3 (Aithal & Ferrante Reference Aithal and Ferrante2020). 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 Reference Aithal and Ferrante2020). Also, Aithal, Tipirneni & Ferrante (Reference Aithal, Tipirneni and Ferrante2022) 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 (Reference Dodd and Ferrante2014).
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
\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
\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 Reference Brackbill, Kothe and Zemach1992),
\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 Reference Cummins, Francois and Kothe2005) with improvements developed by López et al. (Reference López, Zanzi, Gómez, Zamora, Faura and Hernández2009).
The solution algorithm proceeds by calculating three intermediate velocities for the three stages of the RK3 algorithm using the FastRK3 method of Aithal & Ferrante (Reference Aithal and Ferrante2020) as
$$\begin{gather} \boldsymbol{u}_1^* = \boldsymbol{u}^n + \tfrac{1}{3}\Delta t \boldsymbol{\mathcal{M}}(\boldsymbol{u}^n), \end{gather}$$
$$\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}$$
$$\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
\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
\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
\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 Reference Dodd and Ferrante2014)
\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 Reference Dong and Shen2012) as
\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 Reference Schmidt, Schumann and Volkert1984). Finally, we update the velocity field by applying the pressure correction to 
$\boldsymbol {\check {u}}_3^*$ as
\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 Reference Youngs1982). 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 Reference Miller and Colella2002) and Youngs’ method (Youngs Reference Youngs1982) known as the mixed-Youngs-centred method (Aulisa et al. Reference Aulisa, Manservisi, Scardovelli and Zaleski2007).
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:
\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
\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 (Reference Weymouth and Yue2010), which is mass conserving, and wisps are redistributed and suppressed using the method of Baraldi, Dodd & Ferrante (Reference Baraldi, Dodd and Ferrante2014).
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
\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 Reference Schumann1977).
The initial energy spectrum at time 
$t=0$ is prescribed as (Pope Reference Pope2000, § 6.5.3)
\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
\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
\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 Reference Pearson1959; Moffatt Reference Moffatt1965; Kasbaoui, Koch & Desjardins Reference Kasbaoui, Koch and Desjardins2019b). 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 (Reference Ahmed and Elghobashi2000). 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
\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 Reference Freund and Ferrante2019). 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
\begin{equation} \frac{\mathrm{d} k}{\mathrm{d} t} = {{\mathcal{P}}} -\varepsilon + {{\varPsi'_{\sigma}}}, \end{equation}where
$$\begin{gather} k(t)\equiv \tfrac{1}{2} \langle \rho u_j u_j \rangle, \end{gather}$$
$$\begin{gather}{{\mathcal{P}}}(t) \equiv{-}S \langle \rho u w \rangle, \end{gather}$$
$$\begin{gather}\varepsilon (t) \equiv \frac{1}{{\textit{Re}}} \langle {\mathsf{T}}'_{ij}{\mathsf{S}}'_{ij} \rangle, \end{gather}$$
$$\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
\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
\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
\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
\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
\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
\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).
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
\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 (Reference Kida and Tanaka1992) 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 D2 and D4. 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 (Reference Dodd and Ferrante2016). 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
\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,
\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
\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$:
\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
\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 (Reference Dodd and Ferrante2016) 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 Reference Ferrante2022), as, for example, has been done by Freund & Ferrante (Reference Freund and Ferrante2021) with their mixed artificial neural network approach based on the work by Dodd & Ferrante (Reference Dodd and Ferrante2016) and Freund & Ferrante (Reference Freund and Ferrante2019).
In order to perform the DNS study of DLHST, first we developed FastRK3P
$^*$ (§ 2.2.1) by combining FastP
$^*$ (Dodd & Ferrante Reference Dodd and Ferrante2014) with FastRK3 (Aithal & Ferrante Reference Aithal and Ferrante2020). 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 Reference Weymouth and Yue2010; Baraldi et al. Reference Baraldi, Dodd and Ferrante2014), 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 (Reference Dodd and Ferrante2016), 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.
(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$.(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.
(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).
(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).(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).
(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 }}}$.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.647.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Jump conditions at the droplet interface and integral form of the kinetic energy equation for a two-fluid flow
We consider the incompressible flow of two immiscible fluids separated by a common interface in the absence of any body forces and without phase change and gravity, starting from the approach used in Appendix A by Dodd & Ferrante (Reference Dodd and Ferrante2016) for isotropic turbulence to extend it to HST. The geometry of the control volume we consider, which is adapted from Joseph (Reference Joseph1976, p. 242), is illustrated in figure 23. The control volume, 
$\mathcal {V} (t)$, is a material volume, i.e. fluid elements can not cross its boundaries; 
$\mathcal {V}(t)$ consists of two volumes of fluid, 
$\mathcal {V}_c(t)$ and 
$\mathcal {V}_d(t)$ (e.g. the carrier and droplet fluid), separated by an interface 
$\varSigma (t)$, such that 
$\mathcal {V}=\mathcal {V}_c \cup \mathcal {V}_d$. The volumes 
$\mathcal {V}_c$ and 
$\mathcal {V}_d$ are bounded by 
$\partial \mathcal {V}_c(t)$ and 
$\partial \mathcal {V}_d(t)$, respectively, the boundary of 
$\mathcal {V}(t)$ is 
$\partial \mathcal {V} = \partial \mathcal {V}_c \cup \partial \mathcal {V}_d - \varSigma$ and the interface satisfies 
$\varSigma = \partial \mathcal {V}_c \cap \partial \mathcal {V}_d$. The unit normals to 
$\partial \mathcal {V}_c$ and 
$\partial \mathcal {V}_d$ are 
$\boldsymbol {n}_c$ and 
$\boldsymbol {n}_d$, respectively. Here, 
$\boldsymbol {n}$ is a unit vector normal to 
$\varSigma$ that is directed from the carrier fluid to the droplet fluid, and consequently, 
$\boldsymbol {n}=\boldsymbol {n}_c$ for 
$\boldsymbol {x} \in \varSigma$. Here, 
$\partial \varSigma$ is a contact line satisfying 
$\partial \varSigma = \varSigma \cap \partial \mathcal {V}$; 
$\boldsymbol {t}_S$ is a unit vector tangent to 
$\partial \varSigma$; 
$\boldsymbol {m}$ is a unit vector perpendicular to 
$\partial \varSigma$ and pointing out of 
$\varSigma$; 
$\boldsymbol {n}$, 
$\boldsymbol {t}_S$ and 
$\boldsymbol {m}$ are defined such that they form an orthonormal set (e.g. 
$\boldsymbol {m}= \boldsymbol {n} \times \boldsymbol {t}_S$).
Figure 23. Control volume 
$\mathcal {V}(t)$ containing an interface 
$\varSigma (t)$ separating two immiscible volumes of fluid, 
$\mathcal {V}_c(t)$ and 
$\mathcal {V}_d(t)$.
Note: in the following subsections we refer to quantities with ‘
$d$’ subscript as droplet quantities; however, we have made no assumptions about the density ratio and viscosity ratio, and therefore the following equations not only hold for droplet-laden flows but also for bubble-laden flows and, in general, for the incompressible flow of two immiscible fluids separated by an interface.
Following Reynolds decomposition, the fluid velocity, 
$\boldsymbol {U}$, can be decomposed as
\begin{equation} \boldsymbol{U} = \bar{\boldsymbol{U}} + \boldsymbol{u}, \end{equation}
where 
$\bar {\boldsymbol {U}}$ is the mean fluid velocity and 
$\boldsymbol {u}$ is the fluctuating velocity about the mean. For HST,
\begin{equation} \left.\begin{gathered}
\bar{U}_x = S z, \\ \bar{U}_y = 0, \\
\bar{U}_z = 0,
\end{gathered}\right\}\end{equation}
and the fluctuating velocity 
$\boldsymbol {u}$ satisfies the two-fluid incompressible Navier–Stokes equations (2.1a) and (2.1b).
A.1. Conservation of mass
The principle of conservation of mass states that the mass of a material volume does not change, i.e.
\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} \int_\mathcal{V} \rho \mathrm{d} \mathcal{V} = 0. \end{equation}
By taking the limit 
$\partial \mathcal {V} \to \varSigma$, and assuming that there is no mass flux across the interface, one obtains that the normal components of velocity are equal (Aris Reference Aris1989, p. 236), i.e.
\begin{equation} \left.\begin{gathered} \boldsymbol{U}_d \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{U}_c \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{U}_\varSigma \boldsymbol{\cdot} \boldsymbol{n}, \quad \boldsymbol{x} \in \varSigma, \\ ( \bar{\boldsymbol{U}} + \boldsymbol{u}_d ) \boldsymbol{\cdot} \boldsymbol{n} = ( \bar{\boldsymbol{U}} + \boldsymbol{u}_c ) \boldsymbol{\cdot} \boldsymbol{n} = ( \bar{\boldsymbol{U}} + \boldsymbol{u}_\varSigma ) \boldsymbol{\cdot} \boldsymbol{n}, \quad\boldsymbol{x} \in \varSigma, \\ \boldsymbol{u}_d \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{u}_c \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{u}_\varSigma \boldsymbol{\cdot} \boldsymbol{n}, \quad\boldsymbol{x} \in \varSigma. \end{gathered}\right\} \end{equation}Also for viscous fluids under standard operating conditions, it is an experimentally observed fact that the two fluids do not slip and, therefore, the velocity is continuous across the interface,
\begin{equation} \left.\begin{gathered} \boldsymbol{U}_d = \boldsymbol{U}_c = \boldsymbol{U}_\varSigma , \quad \boldsymbol{x} \in \varSigma, \\ \bar{\boldsymbol{U}} + \boldsymbol{u}_d = \bar{\boldsymbol{U}} + \boldsymbol{u}_c = \bar{\boldsymbol{U}} + \boldsymbol{u}_\varSigma , \quad \boldsymbol{x} \in \varSigma, \\ \boldsymbol{u}_d = \boldsymbol{u}_c = \boldsymbol{u}_\varSigma , \quad \boldsymbol{x} \in \varSigma, \end{gathered}\right\}\end{equation}
which we rewrite using jump notation (i.e. 
$[\kern-1pt[ \phi ]\kern-1pt] = \phi _c - \phi _d$),
\begin{equation} \left.\begin{gathered} {[\kern-1pt[ \boldsymbol{U} ]\kern-1pt]} =0, \\ [\kern-1pt[ \bar{\boldsymbol{U}} + \boldsymbol{u} ]\kern-1pt] =0, \\ [\kern-1pt[ \boldsymbol{u} ]\kern-1pt] =0. \end{gathered}\right\}\end{equation}A.2. Conservation of momentum
The conservation equation for the linear momentum of 
$\mathcal {V}$ is
\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \int_\mathcal{V} \rho \boldsymbol{U} \,\mathrm{d} \mathcal{V} = \int_{\partial \mathcal{V}} \boldsymbol{\tau} \boldsymbol{n} \, \mathrm{d} \mathcal{A} + \int_{\partial \varSigma} \sigma \boldsymbol{m} \, \mathrm{d} \ell, \end{equation}
where 
$\boldsymbol {\tau }$ is the fluid stress tensor, which for an incompressible Newtonian fluid is
\begin{align} \boldsymbol{\tau} &={-}p \boldsymbol{\mathsf{I}} + \mu ( \boldsymbol{\nabla} \boldsymbol{U} + (\boldsymbol{\nabla} \boldsymbol{U})^{\rm T} ) \nonumber\\ &={-}p \boldsymbol{\mathsf{I}} + 2 \mu \boldsymbol{\mathsf{S}}, \end{align}
where 
$\boldsymbol{\mathsf{I}}$ is the identity tensor, 
$\boldsymbol{\mathsf{S}} \equiv \tfrac {1}{2} ( \boldsymbol {\nabla } \boldsymbol {U} + (\boldsymbol {\nabla } \boldsymbol {U})^\textrm {T})$ is the strain-rate tensor and 
$\mathrm {d} \ell$ is an infinitesimal arc length (not to be confused with the integral length scale of turbulence, 
$\ell$). On the left-hand side of (A7) is the rate of change of momentum, and on the right-hand side the two terms represent, respectively, the force acting on the boundary due to fluid stress and the force due to surface tension. We note that the last term in (A7) is a line integral along the contact line 
$\partial \varSigma$. This term can be converted from a line integral to a surface integral by using the theorem for curved surfaces, i.e.
\begin{equation} \int_{\partial \varSigma} \sigma \boldsymbol{m} \, \mathrm{d} \ell = \int_{\varSigma} \boldsymbol{\nabla}_s \sigma \, \mathrm{d} \mathcal{A} + \int_{\varSigma} \kappa \sigma \boldsymbol{n} \, \mathrm{d} \mathcal{A} , \end{equation}
where 
$\boldsymbol {\nabla }_s$ is the surface gradient defined as
\begin{equation} \boldsymbol{\nabla}_s = \boldsymbol{\nabla} - \boldsymbol{n} ( \boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{\nabla}). \end{equation}
Combining (A7) and (A9), using the divergence theorem, and accounting for a discontinuity in 
$\boldsymbol {\tau }$ yields
\begin{align} 0 &= \int_\mathcal{V} \left ( \rho \frac{\mathrm{d} \boldsymbol{U}}{\mathrm{d} t} -\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau} \right ) \mathrm{d} \mathcal{V} \nonumber\\ &\quad + \int_\varSigma \left ( [\kern-1pt[ \boldsymbol{\tau} ]\kern-1pt] \boldsymbol{n} - \boldsymbol{\nabla}_s \sigma - \kappa \sigma \boldsymbol{n} \right ) \mathrm{d} \mathcal{A} . \end{align}Using (A1) and (A2), (A11) can be rewritten as
\begin{align} 0 &= \int_\mathcal{V} \left ( \rho \left ( \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} + S z \frac{\partial \boldsymbol{u}}{\partial x} + Su_3 \hat{\boldsymbol{i}} \right ) - \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau} \right ) \mathrm{d} \mathcal{V} \nonumber\\ &\quad + \int_\varSigma \left ( [\kern-1pt[ \boldsymbol{\tau} ]\kern-1pt] \boldsymbol{n} - \boldsymbol{\nabla}_s \sigma - \kappa \sigma \boldsymbol{n} \right ) \mathrm{d} \mathcal{A} . \end{align}
By taking the limit 
$\partial \mathcal {V} \to \varSigma$ and noting that 
$\varSigma$ was chosen arbitrarily, we obtain the jump condition for momentum at 
$\varSigma$,
\begin{equation} {[\kern-1pt[ \boldsymbol{\tau} ]\kern-1pt]} \boldsymbol{n} = \boldsymbol{\nabla}_s \sigma + \kappa \sigma \boldsymbol{n}, \quad \boldsymbol{x} \in \varSigma. \end{equation}Using (A8), the normal stress boundary condition is
\begin{equation} [\kern-1pt[ -p + \boldsymbol{n} \boldsymbol{\cdot} 2 \mu \boldsymbol{\mathsf{S}} \boldsymbol{n} ]\kern-1pt] = \sigma \kappa \end{equation}and the tangential stress boundary condition is
\begin{equation} [\kern-1pt[ \boldsymbol{t}_k \boldsymbol{\cdot} 2 \mu \boldsymbol{\mathsf{S}} \boldsymbol{n} ]\kern-1pt] = \boldsymbol{t}_k \boldsymbol{\cdot} \boldsymbol{\nabla}_s \sigma, \ \ \ k=1,2, \end{equation}
where 
$\boldsymbol {t}_1$ and 
$\boldsymbol {t}_2$ are two unit vectors that are tangent to 
$\varSigma$ and orthogonal to 
$\boldsymbol {n}$. When the surface tension coefficient is constant, (A15) simplifies to
\begin{equation} [\kern-1pt[ \boldsymbol{t}_k \boldsymbol{\cdot} 2 \mu \boldsymbol{\mathsf{S}} \boldsymbol{n} ]\kern-1pt] = 0, \quad k=1,2, \end{equation}
which shows that the tangential stress is continuous across the interface and, if 
$\mu _c \ne \mu _d$, then the rate of strain 
$\boldsymbol{\mathsf{S}}$ is discontinuous at the interface.
A.3. Balance equation of kinetic energy
To derive the kinetic energy balance equation for 
$\mathcal {V}$, we begin with the momentum conservation equation (A12) for 
$\mathcal {V}_c$ and 
$\mathcal {V}_d$,
\begin{equation} \rho \left ( \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} + S z \frac{\partial \boldsymbol{u}}{\partial x} + Su_3 \hat{\boldsymbol{i}} \right ) = \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau} , \quad \boldsymbol{x} \in \mathcal{V}_c \ \text{or}\ \mathcal{V}_d. \end{equation}Using (A1) and (A2), the divergence of the fluid stress tensor can be rewritten in terms of only the fluctuating velocity field to obtain
\begin{align}
\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau} & =
\boldsymbol{\nabla}\boldsymbol{\cdot} ({-}p
\boldsymbol{\mathsf{I}} + \mu (
\boldsymbol{\nabla} \boldsymbol{U} + (\boldsymbol{\nabla}
\boldsymbol{U})^{\rm T} ) )\nonumber\\ & =
\boldsymbol{\nabla}\boldsymbol{\cdot} ({-}p
\boldsymbol{\mathsf{I}} + \mu (
\boldsymbol{\nabla} ( \bar{\boldsymbol{U}} +
\boldsymbol{u} ) + (\boldsymbol{\nabla} (
\bar{\boldsymbol{U}} + \boldsymbol{u} ))^{\rm T}
) )\nonumber\\ & =
\boldsymbol{\nabla}\boldsymbol{\cdot} ({-}p
\boldsymbol{\mathsf{I}} + \mu (
\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla}
\boldsymbol{u})^{\rm T} ) ) +
\boldsymbol{\nabla}\boldsymbol{\cdot} ( \mu (
\boldsymbol{\nabla} \bar{\boldsymbol{U}} + (
\boldsymbol{\nabla} \bar{\boldsymbol{U}} ) ^{\rm T}
) ) \nonumber\\ & =
\boldsymbol{\nabla}\boldsymbol{\cdot} ({-}p
\boldsymbol{\mathsf{I}} + \mu (
\boldsymbol{\nabla} \boldsymbol{u} + ( \boldsymbol{\nabla}
\boldsymbol{u})^{\rm T} ) ) +
\boldsymbol{\nabla} \boldsymbol{\cdot} ( \mu S ( \boldsymbol{e}_1 \boldsymbol{e}_3 + \boldsymbol{e}_3
\boldsymbol{e}_1 ) )\nonumber\\ & =
\boldsymbol{\nabla} \boldsymbol{\cdot} ({-}p
\boldsymbol{\mathsf{I}} + \mu (
\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla}
\boldsymbol{u})^{\rm T} ) ) \nonumber\\ & =
\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau'}.
\end{align}
Here, 
$\boldsymbol {\tau '}$ is the fluid stress tensor of the fluctuating velocity field defined as
\begin{align}
\boldsymbol{\tau'} &={-}p \boldsymbol{\mathsf{I}} + \mu
( \boldsymbol{\nabla} \boldsymbol{u} +
(\boldsymbol{\nabla} \boldsymbol{u})^{\rm T} )\nonumber\\
&={-}p \boldsymbol{\mathsf{I}} + 2 \mu
\boldsymbol{\mathsf{S}}',
\end{align}
where 
$\boldsymbol{\mathsf{S}}' \equiv \frac {1}{2} (\boldsymbol {\nabla } \boldsymbol {u} + (\boldsymbol {\nabla } \boldsymbol {u} ) ^\textrm {T} )$ is the strain-rate tensor of the fluctuating velocity field. Using (A1), the relation between 
$\boldsymbol {\tau }$ and 
$\boldsymbol {\tau '}$ is
\begin{equation} \boldsymbol{\tau} = \langle \boldsymbol{\tau} \rangle + \boldsymbol{\tau'}, \end{equation}
where 
$\langle \boldsymbol {\tau } \rangle = \mu \left ( \boldsymbol {\nabla } \bar {\boldsymbol {U}} + \left ( \boldsymbol {\nabla } \bar {\boldsymbol {U}} \right ) ^\textrm {T} \right )$. Combining (A17) and (A18) and taking the dot product with 
$\boldsymbol {u}$ yields
\begin{align} &\rho \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} ((\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u}))/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right )\nonumber\\ &\quad = (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau'}) \boldsymbol{\cdot} \boldsymbol{u} , \quad \boldsymbol{x} \in \mathcal{V}_c\ \text{or} \ \mathcal{V}_d. \end{align}
We integrate (A21) over 
$\mathcal {V}_c$ and 
$\mathcal {V}_d$, use the identity 
$\boldsymbol {u} \boldsymbol {\cdot } (\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {\tau '})= \boldsymbol {\nabla }\boldsymbol {\cdot } (\boldsymbol {\tau '} \boldsymbol {u}) - \boldsymbol {\tau '} : \boldsymbol {\nabla } \boldsymbol {u}$ and use the divergence theorem to obtain
\begin{equation} \left.\begin{gathered} \int_{\mathcal{V}_c} \rho_c \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) \mathrm{d} \mathcal{V} \\ = \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \, \mathrm{d} \mathcal{A} - \int_{\mathcal{V}_c} \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \, \mathrm{d} \mathcal{V}, \\ \int_{\mathcal{V}_d} \rho_d \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) \mathrm{d} \mathcal{V} \\ = \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \,\mathrm{d} \mathcal{A} - \int_{\mathcal{V}_d} \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \,\mathrm{d} \mathcal{V}. \end{gathered}\right\} \end{equation}
Adding together the kinetic energy equations for the carrier and droplet fluid in (A22) yields the evolution equation for the kinetic energy of 
$\mathcal {V}$,
\begin{align} &\int_\mathcal{V} \rho \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) \mathrm{d} \mathcal{V}\nonumber\\ &\quad = \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \,\mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \,\mathrm{d} \mathcal{A} - \int_{\mathcal{V}} \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \,\mathrm{d} \mathcal{V} . \end{align}
We then use the following transformation to account for the jump in 
$\boldsymbol {\tau '}$ at the interface:
\begin{align} \int_{\partial \mathcal{V}} \boldsymbol{\tau'} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \,\mathrm{d} \mathcal{A} &= \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \,\mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \, \mathrm{d} \mathcal{A} \nonumber\\ &\quad - \int_{\varSigma} ( \boldsymbol{\tau'}\boldsymbol{u} )_c \boldsymbol{\cdot} \boldsymbol{n}_c \, \mathrm{d} \mathcal{A} - \int_{\varSigma} ( \boldsymbol{\tau'}\boldsymbol{u} )_d \boldsymbol{\cdot} \boldsymbol{n}_d \,\mathrm{d} \mathcal{A}\nonumber\\ &= \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \,\mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \,\mathrm{d} \mathcal{A}\nonumber\\ &\quad - \int_{\varSigma} ( \boldsymbol{\tau'}\boldsymbol{u} )_c \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A} + \int_{\varSigma} ( \boldsymbol{\tau'}\boldsymbol{u} )_d \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A}\nonumber\\ &= \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \, \mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \, \mathrm{d} \mathcal{A} - \int_{\varSigma} [\kern-1pt[ \boldsymbol{\tau'} \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A}. \end{align}
In the second line we have used the fact that 
$\boldsymbol {n}_c=\boldsymbol {n}$ and 
$\boldsymbol {n}_d=-\boldsymbol {n}$ for 
$\boldsymbol {x} \in \varSigma$. Combining (A23) and (A24) yields
\begin{align} &\int_\mathcal{V} \rho \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) \mathrm{d} \mathcal{V}\nonumber\\ &\quad = \int_{\partial \mathcal{V}} \boldsymbol{\tau'} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \,\mathrm{d} \mathcal{A} - \int_{\mathcal{V}} \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \, \mathrm{d} \mathcal{V} + \int_{\varSigma} [\kern-1pt[ \boldsymbol{\tau'} \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} \,\mathrm{d} \mathcal{A}. \end{align}
The work due to surface tension contributes to the last term on the right-hand side of (A25). This is made clear by dotting (A13) with 
$\boldsymbol {u}$,
\begin{equation} \left.\begin{gathered} {[\kern-1pt[ \boldsymbol{\tau} \boldsymbol{u} ]\kern-1pt]} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}_s \sigma + \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}, \ \ \ \boldsymbol{x} \in \varSigma, \\ [\kern-1pt[ \boldsymbol{\tau'} \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}_s \sigma + \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} - [\kern-1pt[ \langle \boldsymbol{\tau} \rangle \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n}, \quad \boldsymbol{x} \in \varSigma, \\ [\kern-1pt[ \boldsymbol{\tau'} \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}_s \sigma + \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}, \end{gathered}\right\} \end{equation}
where we have used (A5), (A20) and the fact that 
$[\kern-1pt[ \langle \boldsymbol {\tau } \rangle \boldsymbol {u} ]\kern-1pt]$ is zero because 
$\langle \boldsymbol {\tau } \rangle$ is constant and 
$\boldsymbol {u}$ is continuous across the interface. Combining (A25) and (A26) and using the divergence theorem yields the integral form of the kinetic energy equation for a two-fluid flow,
\begin{align} & \int_\mathcal{V} \left( \rho \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) - \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\tau'} \boldsymbol{u}) + \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \right) \mathrm{d} \mathcal{V} \nonumber\\ &\quad = \int_{\varSigma} \left ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}_s \sigma + \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \right ) \mathrm{d} \mathcal{A}. \end{align}Appendix B. Turbulence kinetic energy equations in DLHST
We now derive the balance equations for the TKE of the carrier fluid 
$\mathcal {V}_c$, droplet fluid 
$\mathcal {V}_d$ and combined fluid 
$\mathcal {V}$, starting from the approach used in Appendix B by Dodd & Ferrante (Reference Dodd and Ferrante2016) for isotropic turbulence to extend it to HST. We consider decaying homogeneous isotropic turbulence laden with droplets with a constant surface tension coefficient, 
$\sigma$, in a periodic domain. Up until now, we have considered only two volumes of fluid, 
$\mathcal {V}_c$ and 
$\mathcal {V}_d$, separated by an interface as depicted in figure 23. When deriving TKE equations for droplet-laden turbulence, we take the control volume 
$\mathcal {V}$ to be the two-fluid flow, which includes 
$N_d$ droplets with volumes 
$\mathcal {V}_{d}^{(1)}, \mathcal {V}_{d}^{(2)}, \ldots, \mathcal {V}_{d}^{(N_d)}$ immersed in the carrier-fluid volume 
$\mathcal {V}_c$. An example of this configuration with 
$N_d=4$ is shown in figure 41 of Dodd & Ferrante (Reference Dodd and Ferrante2016).
The two-fluid TKE, 
$k$, is defined as the spatial average of the kinetic energy (per unit volume) of the fluctuating velocity field, 
$\boldsymbol {u}(\boldsymbol {x},t)$, (i.e. 
$\boldsymbol {u}=\boldsymbol {U}-\bar {\boldsymbol {U}}$, where 
$\boldsymbol {U}$ is the total velocity and 
$\bar {\boldsymbol {U}}$ is its mean),
\begin{equation} k \equiv \frac{1}{2} \frac{1}{\mathcal{V}} \int_\mathcal{V} \rho \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \, \mathrm{d} \mathcal{V} = \frac{1}{2}\langle \rho \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \rangle, \end{equation}
where 
$\langle \cdots \rangle$ denotes spatial averaging of the enclosed quantity. Likewise, we define the carrier-fluid TKE as
\begin{equation} k_c \equiv \frac{1}{2} \frac{1}{\mathcal{V}_c} \int_{\mathcal{V}_c} \rho_c \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \,\mathrm{d} \mathcal{V} = \frac{1}{2} \langle \rho_c \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \rangle_c \end{equation}and the droplet-fluid TKE as
\begin{equation} k_d \equiv \frac{1}{2} \frac{1}{\mathcal{V}_{d}} \sum_{n=1}^{N_d} \int_{\mathcal{V}_{d}^{(n)}} \rho_d \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \,\mathrm{d} \mathcal{V} = \frac{1}{2} \langle \rho_d \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \rangle_d , \end{equation}
where 
$\mathcal {V}_{d}^{(n)}$ is the volume of the 
$n$th droplet and 
$\mathcal {V}_{d}$ is the total volume of the droplet fluid (i.e. 
$\mathcal {V}_{d} = \sum _{n=1}^{N_d} \mathcal {V}_{d}^{(n)}$). The summation is performed over the total number of droplets 
$N_d$.
We first derive the TKE evolution equation for the carrier fluid. From (A22) and invoking incompressibility (
$\textrm {d}\rho /\textrm {d}t=0$) yields
\begin{align} &\int_{\mathcal{V}_c} \frac{\partial (\rho_c u_i u_i)/2}{\partial t} + \frac{\partial (\rho_c u_i u_i u_j)/2}{\partial x_j} + S z \frac{\partial (\rho_c u_i u_i)/2}{\partial x} + S \rho_c u_1 u_3 \,\mathrm{d} \mathcal{V}\nonumber\\ &\quad = \int_{\partial \mathcal{V}_c} \tau'_{ij} u_i n_{j,c} \, \mathrm{d} \mathcal{A} - \int_{\mathcal{V}_c} \tau'_{ij} \frac{\partial u_i}{\partial x_j} \, \mathrm{d} \mathcal{V} . \end{align}Then, by applying the divergence theorem to the first term on the right-hand side of (B4) and using (A19) we obtain
\begin{align} & \int_{\mathcal{V}_c} \frac{\partial (\rho_c u_i u_i)/2}{\partial t} + \frac{\partial (\rho_c u_i u_i u_j)/2}{\partial x_j} + S z \frac{\partial (\rho_c u_i u_i)/2}{\partial x} + S \rho_c u_1 u_3\,\mathrm{d} \mathcal{V}\nonumber\\ &\quad = \int_{\mathcal{V}_c} \left [ \frac{\partial ( \boldsymbol{\mathsf{T}}'_{ij} u_i)}{\partial x_j} - \frac{\partial(u_i p)}{\partial x_i} \right ] \mathrm{d} \mathcal{V} - \int_{\mathcal{V}_c} \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \,\mathrm{d} \mathcal{V} , \end{align}
where 
$\boldsymbol{\mathsf{T}}'_{ij}=2 \mu \boldsymbol{\mathsf{S}}'_{ij}$. For immiscible fluids, there is no convective transport of TKE between the carrier fluid and droplet fluid, and, therefore, the second and third terms on the left-hand side of (B5) are zero. By dividing (B5) by 
$\mathcal {V}_c$ and rewriting the resulting equation in non-dimensional form, we obtain
\begin{equation}
\underbrace{\frac{\mathrm{d} k_c}{\mathrm{d} t}}_{
\begin{smallmatrix} \textit{rate of change of} \\
\textit{carrier-fluid TKE} \end{smallmatrix}} =
\underbrace{{{\mathcal{P}}}_c}_{ \begin{smallmatrix}
\textit{carrier-fluid} \\ \textit{production}
\end{smallmatrix}} -\underbrace{\varepsilon_c}_{
\begin{smallmatrix} \textit{carrier-fluid} \\
\textit{dissipation rate} \end{smallmatrix}} +\underbrace{
T_{\nu,c}}_{ \begin{smallmatrix} \textit{carrier-fluid} \\
\textit{viscous power} \end{smallmatrix}} +\underbrace{
T_{p,c}}_{ \begin{smallmatrix} \textit{carrier-fluid} \\
\textit{pressure power} \end{smallmatrix}},
\end{equation}
where
\begin{align} {{\mathcal{P}}}_c \equiv{-} S \langle \rho_c u_1 u_3 \rangle_c, \quad \varepsilon_c \equiv \frac{1}{{\textit{Re}}} \langle \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \rangle_c, \quad T_{\nu,c} \equiv \frac{1}{{\textit{Re}}} \frac{\partial \langle \boldsymbol{\mathsf{T}}'_{ij} u_i \rangle_c }{\partial x_j}, \quad T_{p,c} \equiv{-} \frac{\partial \langle u_i p \rangle_c }{\partial x_i} . \end{align}The terms in (B6) are, from left to right, the rate of change of carrier-fluid TKE, the production of TKE in the carrier fluid, the dissipation rate of TKE in the carrier fluid, the pressure power of the carrier fluid (transport of TKE due to pressure) and the viscous power of the carrier fluid (transport of TKE due to viscous stresses). The TKE equation for the droplet fluid is found by writing the TKE equation for each droplet and then summing over all droplets. The final result is
\begin{equation}
\underbrace{\frac{\mathrm{d} k_d}{\mathrm{d} t}}_{
\begin{smallmatrix} \textit{rate of change of} \\
\textit{droplet-fluid TKE} \end{smallmatrix}} =
\underbrace{{{\mathcal{P}}}_d}_{ \begin{smallmatrix}
\textit{droplet-fluid} \\ \textit{production}
\end{smallmatrix}} -\underbrace{\varepsilon_d}_{
\begin{smallmatrix} \textit{droplet-fluid} \\
\textit{dissipation rate} \end{smallmatrix}} +\underbrace{
T_{\nu,d}}_{ \begin{smallmatrix} \textit{droplet-fluid} \\
\textit{viscous power} \end{smallmatrix}} +\underbrace{
T_{p,d}}_{ \begin{smallmatrix} \textit{droplet-fluid} \\
\textit{pressure power} \end{smallmatrix}},
\end{equation}
where
\begin{align} {{\mathcal{P}}}_d \equiv{-} S \langle \rho_d u_1 u_3 \rangle_d, \quad \varepsilon_d \equiv \frac{1}{{\textit{Re}}} \langle \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \rangle_d, \quad T_{\nu,d} \equiv \frac{1}{{\textit{Re}}} \frac{\partial \langle \boldsymbol{\mathsf{T}}'_{ij} u_i \rangle_d }{\partial x_j}, \quad T_{p,c} \equiv{-} \frac{\partial \langle u_i p \rangle_d }{\partial x_i} . \end{align}The TKE equation for the two-fluid flow, which includes the interface, is found by adding the two equations of (A22), resulting in
\begin{align} &\int_\mathcal{V} \rho \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) \mathrm{d} \mathcal{V}\nonumber\\ &\quad = \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \, \mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \, \mathrm{d} \mathcal{A} - \int_{\mathcal{V}} \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \, \mathrm{d} \mathcal{V} . \end{align}
Because 
$\mathcal {V}$ is periodic, the left-hand side of (A24) is zero such that
\begin{equation} \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \, \mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \, \mathrm{d} \mathcal{A} = \int_{\varSigma} [\kern-1pt[ \boldsymbol{\tau'} \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A}. \end{equation}
For constant 
$\sigma$, integration of (A26) over 
$\varSigma$ gives
\begin{equation} \int_{\varSigma} [\kern-1pt[ \boldsymbol{\tau'} \boldsymbol{u} ]\kern-1pt] \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A} = \int_{\varSigma} \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A}. \end{equation}Substituting (B12) in (B11) yields
\begin{equation} \int_{\partial \mathcal{V}_c} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_c \, \mathrm{d} \mathcal{A} + \int_{\partial \mathcal{V}_d} \boldsymbol{\tau'}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n}_d \, \mathrm{d} \mathcal{A} = \int_{\varSigma} \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A}, \end{equation}and substituting (B13) into (B10) gives
\begin{align} &\int_\mathcal{V} \rho \left ( \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 + S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} + S u_1 u_3 \right ) \mathrm{d} \mathcal{V}\nonumber\\ &\quad = \int_{\varSigma} \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A} - \int_{\mathcal{V}} \boldsymbol{\tau'} : \boldsymbol{\nabla}\boldsymbol{u} \, \mathrm{d} \mathcal{V}. \end{align}We transform the surface integral term in (B14) to a volume integral as
\begin{equation} \int_\varSigma \kappa \sigma \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} \mathcal{A} = \frac{1}{{\textit{We}}} \int_\mathcal{V} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{f}_\sigma \, \mathrm{d} \mathcal{V}, \end{equation}
where we recall that 
$\boldsymbol {f}_\sigma \equiv \kappa \delta (s) \boldsymbol {n}$ is the surface tension force, and that 
$\sigma = 1/{\textit {We}}$ is constant. By using (B15) in (B14), invoking incompressibility and noting that because 
$\mathcal {V}$ is periodic, the convective terms from the fluctuating velocity field in (B14) are equal to zero,
\begin{equation} \left.\begin{gathered} \int_{\mathcal{V}} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2 \, \mathrm{d} \mathcal{V} =0, \\ \int_{\mathcal{V}} S z \frac{\partial (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u})/2}{\partial x} \, \mathrm{d} \mathcal{V} =0, \end{gathered}\right\} \end{equation}we obtain, after reorganizing,
\begin{equation} \int_{\mathcal{V}} \frac{\partial (\rho u_i u_i)/2}{\partial t} \, \mathrm{d} \mathcal{V} ={-} \int_{\mathcal{V}} S \rho u_1 u_3 \ \mathrm{d} \mathcal{V} - \int_{\mathcal{V}} \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \, \mathrm{d} \mathcal{V} + \frac{1}{{\textit{We}}} \int_{\mathcal{V}} u_j f_{\sigma,j} \, \mathrm{d} \mathcal{V}. \end{equation}
Because 
$\mathcal {V}$ was not chosen arbitrarily in this derivation, the integrals in (B17) must be retained. By dividing (B17) through by 
$\mathcal {V}$, rewriting the equation in non-dimensional form and noting that the time differentiation and integration commute, we obtain the following TKE evolution equation for the two-fluid flow:
\begin{equation}
\underbrace{\frac{\mathrm{d} k}{\mathrm{d} t}}_{
\begin{smallmatrix} \textit{rate of change} \\ \textit{of TKE}
\end{smallmatrix}} = \underbrace{{{\mathcal{P}}}}_{
\textit{production}} -\underbrace{\varepsilon}_{
\textit{dissipation rate}} +\underbrace{
{{\varPsi'_{\sigma}}}}_{ \begin{smallmatrix} \textit{power of
the} \\ \textit{surface tension} \end{smallmatrix}} .\end{equation}
Here,
\begin{equation} {{\mathcal{P}}} \equiv{-}S \langle \rho u_1 u_3 \rangle, \quad \varepsilon \equiv \frac{1}{{\textit{Re}}} \langle \boldsymbol{\mathsf{T}}'_{ij}\boldsymbol{\mathsf{S}}'_{ij} \rangle, \quad {{\varPsi'_{\sigma}}} \equiv \frac{1}{{\textit{We}}} \langle u_j\, f_{\sigma,j} \rangle. \end{equation}The terms in (B18) are, from left to right, the rate of change of TKE for the two-fluid flow, the production of TKE for the two-fluid flow, the dissipation rate of TKE for the two-fluid flow and the power of the surface tension due to the fluctuating velocity. To summarize, we have derived the TKE evolution equations for the
(i) carrier fluid
(B20)\begin{equation} \frac{\mathrm{d} k_c}{\mathrm{d} t} = P_c -\varepsilon_c + T_{\nu,c} + T_{p,c}, \end{equation}(ii) droplet fluid
(B21)\begin{equation} \frac{\mathrm{d} k_d}{\mathrm{d} t} = P_d -\varepsilon_d + T_{\nu,d} + T_{p,d}, \end{equation}(iii) two-fluid flow
(B22)\begin{equation} \frac{\mathrm{d} k}{\mathrm{d} t} = P -\varepsilon + {{\varPsi'_{\sigma}}}. \end{equation}
A useful identity, which comes from using (B15) in (B13) and invoking the divergence theorem, is
\begin{equation} \mathcal{V} {{\varPsi'_{\sigma}}} = \mathcal{V}_c ( T_{\nu,c} + T_{p,c} ) + \mathcal{V}_d ( T_{\nu,d} + T_{p,d} ) , \end{equation}
and after dividing (B23) throughout by 
$\mathcal {V}$, we obtain
\begin{equation} {{\varPsi'_{\sigma}}} = (1- \phi_v) ( T_{\nu,c} + T_{p,c} ) + \phi_v ( T_{\nu,d} + T_{p,d} ), \end{equation}
where 
$\phi _v \equiv \mathcal {V}_d / \mathcal {V}$ is the droplet volume fraction. Equation (B24) shows that the power of the surface tension due to the fluctuating velocity is equal to the sum of the volume weighted carrier and droplet-fluid viscous and pressure powers.
Appendix C. Relationship between the rate of change of total interfacial surface area and the power of the surface tension due to the total velocity for closed surfaces
In this section we derive the relationship between the rate of change of interfacial energy (
$({\mathrm {d}}/{\mathrm {d} t})\int _\varSigma \sigma \,\mathrm {d} \mathcal {A}$) and the power of the surface tension due to the total velocity (
${{\varPsi _{\sigma }}}$), starting from the approach used in Dodd & Ferrante (Reference Dodd and Ferrante2016, Appendix C) for isotropic turbulence to extend it to HST. We begin be deriving the relationship for a single droplet and then generalize it for an arbitrary number of droplets 
$N_d$. Starting from the analog of the Reynolds transport theorem for a surface (Aris Reference Aris1989, p. 230), under the assumption that fluid parcels do not cross the interface, Joseph (Reference Joseph1976, p. 243) derived the following transport equation for the surface tension:
\begin{align} \frac{\mathrm{d}}{\mathrm{d} t} \int_{\varSigma^{(n)}} \sigma\, \mathrm{d} \mathcal{A} = \int_{\varSigma^{(n)}} \left ( \frac{\partial \sigma}{\partial t} - \boldsymbol{U}_\varSigma \boldsymbol{\cdot}\boldsymbol{\nabla}_s \sigma - \kappa \sigma \boldsymbol{U}_\varSigma \boldsymbol{\cdot} \boldsymbol{n} \right ) \mathrm{d} \mathcal{A} + \int_{\partial \varSigma^{(n)}} \sigma \boldsymbol{U}_\varSigma \boldsymbol{\cdot} \boldsymbol{m} \, \mathrm{d} \ell. \end{align}
In our case, 
$\varSigma ^{(n)}$ is the interface of the 
$n$th droplet and 
$\partial \varSigma ^{(n)}$ is the contact line as defined in Appendix A and figure 23. Because 
$\varSigma ^{(n)}$ forms a closed surface, 
$\partial \varSigma ^{(n)}$ is inexistent, and therefore the last term in (C1) is null. Also, if the surface tension coefficient is constant in space and time (i.e. 
$\sigma (\boldsymbol {x},t)=\sigma$), then (C1) reduces to
\begin{equation} \sigma \frac{\mathrm{d}
A^{(n)}}{\mathrm{d} t} = \sigma \int_{\varSigma^{(n)}} -
\kappa \boldsymbol{U}_\varSigma \boldsymbol{\cdot}
\boldsymbol{n} \, \mathrm{d} \mathcal{A}, \end{equation}
where 
$A^{(n)}=\int _{\varSigma ^{(n)}} \,\mathrm {d} \mathcal {A}$ is the surface area of 
$\varSigma ^{(n)}$. We can also transform the right-hand side of (C2) from a surface integral to a volume integral to obtain
\begin{equation}
\sigma \frac{\mathrm{d} A^{(n)}}{\mathrm{d} t} ={-}\sigma\int_{\mathcal{V}^{(n)}}
\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{f}_\sigma \, \mathrm{d} \mathcal{V},
\end{equation}
which relates the rate of change of surface area of droplet 
$n$ to the power of the surface tension due to the total velocity. To derive a relationship for the more general case of multiple droplets, we sum over 
$N_d$ droplets,
\begin{equation} \sum_{n=1}^{N_d}
\biggl(\sigma \frac{\mathrm{d} A^{(n)}}{\mathrm{d} t} \biggr)
={-} \sigma \sum_{n=1}^{N_d} \left(\int_{\mathcal{V}^{(n)}}
\boldsymbol{U} \boldsymbol{\cdot}\boldsymbol{f}_\sigma \, \mathrm{d}\mathcal{V}\right).
\end{equation}
Interchanging the summation and differentiation and defining 
$A \equiv \sum _{n=1}^{N_d} A^{(n)}$ and 
$\mathcal {V} \equiv \sum _{n=1}^{N_d} \mathcal {V}^{(n)}$ in (C4) yields
\begin{equation} \sigma \frac{\mathrm{d} A}{\mathrm{d} t} ={-}\sigma \int_{\mathcal{V}} \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{f}_\sigma \, \mathrm{d} \mathcal{V}. \end{equation}
Using 
$\langle \cdots \rangle$ to denote spatial averaging and dividing (C5) by 
$\mathcal {V}$ yields
\begin{equation} \sigma \frac{1}{\mathcal{V}} \frac{\mathrm{d} A}{\mathrm{d} t} ={-}\sigma \langle \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{f}_\sigma \rangle. \end{equation}
Recalling that 
$\sigma =1/{\textit {We}}$ (§ 2.1), and by defining the power of the surface tension due to the total velocity as
\begin{equation} {{\varPsi_{\sigma}}} \equiv \frac{1}{{\textit{We}}} \langle \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{f}_\sigma \rangle, \end{equation}then, (C6) gives
\begin{equation} {{\varPsi_{\sigma}}}
={-}\frac{1}{{\textit{We}}}\frac{1}{\mathcal{V}}
\frac{\mathrm{d} A}{\mathrm{d} t}.
\end{equation}
Open access
