2 Deformation and breakup statistics

The model of Maffettone & Minale (Reference Maffettone and Minale1998) assumes that both the fluid of which the drop is composed and the fluid in which it is immersed are Newtonian. The drop is neutrally buoyant and is transported passively (i.e. it does not affect the surrounding flow), it is ellipsoidal at all times, and its volume is preserved. In addition, the flow about the drop is incompressible and linear. This latter assumption is appropriate for turbulent flows if the size of the drop is smaller than $\ell _{K}$ . The volume fraction is very low, so that hydrodynamic interactions between drops are negligible and attention can be directed to the dynamics of a single drop.

The shape and the orientation of the drop are described by a second-rank symmetric positive definite tensor $\unicode[STIX]{x1D648}$ , whose eigenvectors are the semi-axes of the drop and whose eigenvalues $m_{1}^{2}\geqslant m_{2}^{2}\geqslant m_{3}^{2}$ yield the squared lengths of the same semi-axes. The centre of mass of the drop evolves as a tracer, while the Lagrangian evolution of $\unicode[STIX]{x1D648}$ is given by the following equation:

(2.1) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D648}}=\unicode[STIX]{x1D642}\unicode[STIX]{x1D648}+\unicode[STIX]{x1D648}\unicode[STIX]{x1D642}^{\top }-{\displaystyle \frac{f_{1}(\unicode[STIX]{x1D707})}{\unicode[STIX]{x1D70F}}}[\unicode[STIX]{x1D648}-g(\unicode[STIX]{x1D648})\unicode[STIX]{x1D644}], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D642}=f_{2}(\unicode[STIX]{x1D707})\unicode[STIX]{x1D64E}+\unicode[STIX]{x1D734}$ is an effective velocity gradient; $\unicode[STIX]{x1D734}=[\unicode[STIX]{x1D735}\boldsymbol{u}-(\unicode[STIX]{x1D735}\boldsymbol{u})^{\top }]/2$ and $\unicode[STIX]{x1D64E}=[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\top }]/2$ are the vorticity and rate-of-strain tensors evaluated at the centre of mass of the drop, respectively. Note that here $(\unicode[STIX]{x1D735}\boldsymbol{u})_{ij}=\unicode[STIX]{x2202}_{j}u_{i}$ . The coefficients $f_{1}(\unicode[STIX]{x1D707})$ and $f_{2}(\unicode[STIX]{x1D707})$ depend on the ratio $\unicode[STIX]{x1D707}$ of the viscosity of the drop and that of the external fluid and were chosen in such a way as to match theoretical predictions for small capillary numbers (Maffettone & Minale Reference Maffettone and Minale1998):

(2.2a,b ) $$\begin{eqnarray}\displaystyle f_{1}(\unicode[STIX]{x1D707})={\displaystyle \frac{40(\unicode[STIX]{x1D707}+1)}{(2\unicode[STIX]{x1D707}+3)(19\unicode[STIX]{x1D707}+16)}},\quad f_{2}(\unicode[STIX]{x1D707})={\displaystyle \frac{5}{2\unicode[STIX]{x1D707}+3}}. & & \displaystyle\end{eqnarray}$$

Note that $f_{2}(1)=1$ and hence, for $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D642}=\unicode[STIX]{x1D735}\boldsymbol{u}$ . The last term in (2.1) describes the capillary relaxation to the spherical shape with a time scale $\unicode[STIX]{x1D70F}$ . Thanks to an appropriate choice of the function $g(\unicode[STIX]{x1D648})$ , the same term enforces that $\det \unicode[STIX]{x1D648}$ is constant in time and hence the volume of the drop is preserved. The function $g(\unicode[STIX]{x1D648})$ has the form $g(\unicode[STIX]{x1D648})=3\text{III}_{M}/\text{II}_{M}$ , where $\text{II}_{M}$ and $\text{III}_{M}$ are the second and third invariants of $\unicode[STIX]{x1D648}$ , i.e. $\text{II}_{M}=[(\operatorname{tr}\unicode[STIX]{x1D648})^{2}-\operatorname{tr}\unicode[STIX]{x1D648}^{2}]/2$ and $\text{III}_{M}=\det \unicode[STIX]{x1D648}$ . Maffettone & Minale (Reference Maffettone and Minale1998) also proposed an improved expression of $f_{2}$ that depends on the capillary number and more accurately describes the deformations observed in experiments for large strains and high viscosity ratios. For the sake of simplicity, here we use the coefficients given in (2.2); the improved version of the model of Maffettone & Minale (Reference Maffettone and Minale1998) is discussed in § 4.

This section provides insight into the statistics of drop deformation and breakup in three-dimensional homogeneous isotropic turbulence. We obtain such a turbulent flow by performing direct numerical simulations of the three-dimensional Navier–Stokes equation

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\boldsymbol{u}+\boldsymbol{F} & & \displaystyle\end{eqnarray}$$

for the velocity field $\boldsymbol{u}$ (and pressure $p$ ) augmented with the incompressibility condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ . We use the standard, fully de-aliased pseudo-spectral method on a cubic domain of size $2\unicode[STIX]{x03C0}$ with $512^{3}$ collocation points and periodic boundary conditions. By using these boundary conditions, we do not take into account the interaction of the drops with the walls that confine the fluid. The flow is driven to a non-equilibrium steady state by an external force $\boldsymbol{F}$ with a fixed energy input $\unicode[STIX]{x1D716}$ . Our choice of $\unicode[STIX]{x1D716}$ and kinematic viscosity $\unicode[STIX]{x1D708}$ ensures a Taylor-scale Reynolds number $Re_{\unicode[STIX]{x1D706}}\approx 111$ .

In order to study the deformation of droplets in a turbulent flow, we seed our turbulent, statistically steady, flow with Lagrangian tracers and follow their trajectories, by using a trilinear-interpolation scheme to obtain the tracer velocity from the Eulerian velocity evaluated from (2.3); such trajectories define the motion of the centre of mass of the droplets. We refer the reader to James & Ray (Reference James and Ray2017) for a more detailed description of our numerical procedure.

The capillary number is defined as $Ca=\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}$ , where $\unicode[STIX]{x1D706}$ is the Lyapunov exponent of the flow. This latter represents the average stretching rate in a turbulent flow and provides a measure of the ability of the flow to deform a drop. We calculate $\unicode[STIX]{x1D706}$ by using the fluid velocity gradients along the trajectories and obtain $\unicode[STIX]{x1D706}\approx 4.22\approx 0.15\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}^{-1}$ (where $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ is the Kolmogorov time scale associated with the flow), consistent with earlier results (Bec et al. Reference Bec, Biferale, Boffetta, Cencini, Lanotte, Musacchio and Toschi2006). (Note that Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) defined the capillary number in terms of the root mean square of $\unicode[STIX]{x2202}_{x}u_{x}$ instead of $\unicode[STIX]{x1D706}$ . However, this fact only leads to a different proportionality factor in the definition of $\mathit{Ca}$ , since in isotropic turbulence $\sqrt{\langle (\unicode[STIX]{x2202}_{x}u_{x})^{2}\rangle }=1/\sqrt{15}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ and hence in our case $\sqrt{\langle (\unicode[STIX]{x2202}_{x}u_{x})^{2}\rangle }=1.72\unicode[STIX]{x1D706}$ .)

Equation (2.1) is integrated by using the second-order Adam–Bashforth method with the same time step as for the Navier–Stokes equation. The integration of (2.1) must preserve the positive definite character of $\unicode[STIX]{x1D648}$ . This is achieved by adapting to (2.1) the Cholesky-decomposition method proposed by Vaithianathan & Collins (Reference Vaithianathan and Collins2003) (see appendix A for details). Unless otherwise stated, the initial condition is $\unicode[STIX]{x1D648}(0)=\unicode[STIX]{x1D644}$ . As in Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014), it is assumed that drops break when their aspect ratio $|m_{1}/m_{3}|$ exceeds a threshold value $\unicode[STIX]{x1D6FC}$ . In view of the fact that we are only interested in the dynamics up to the first breakup and do not consider secondary breakup events, drops are removed from the flow as soon as they break. In the simulations presented below, the initial number of drops $N(0)=10^{6}$ .

The deformation of a drop is described in terms of the statistics of $m_{1}^{2}$ , i.e. the squared length of the semi-major axis. Let $p(m_{1}^{2},t)$ be the probability density function (p.d.f.) of $m_{1}^{2}$ and ${\mathcal{P}}(m_{1}^{2})\equiv \int _{0}^{\infty }p(m_{1}^{2},t)\,\text{d}t$ its time integral. Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) showed that ${\mathcal{P}}(m_{1}^{2})$ behaves as a power law for values of $m_{1}^{2}$ smaller than its maximum value (it is easy to check that the conditions $m_{1}^{2}\geqslant m_{2}^{2}\geqslant m_{3}^{2}$ , $m_{1}^{2}m_{2}^{2}m_{3}^{2}=1$ and $|m_{1}/m_{3}|\leqslant \unicode[STIX]{x1D6FC}$ imply that $m_{1}^{2}$ is bounded and, more precisely, $m_{1}^{2}\leqslant \unicode[STIX]{x1D6FC}^{4/3}$ ). The slope increases as a function of $\mathit{Ca}$ for small $\mathit{Ca}$ and saturates to $-1$ when $\mathit{Ca}$ exceeds a critical value $\mathit{Ca}_{c}$ , which for $\unicode[STIX]{x1D707}=1$ was found to be $\mathit{Ca}_{c}=f_{1}(\unicode[STIX]{x1D707})/2$ . Figure 1(a) shows that the behaviour of ${\mathcal{P}}(m_{1}^{2})$ is accurately reproduced in our simulations. The power law is even clearer when the p.d.f of $m_{1}^{2}/m_{3}^{2}$ is considered (figure 1 b).

Figure 1. Time-integrated p.d.f. of (a) the largest eigenvalue of $\unicode[STIX]{x1D648}$ and (b) the ratio of the largest and the smallest eigenvalue of $\unicode[STIX]{x1D648}$ for $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D6FC}=10^{3}$ and different values of $\mathit{Ca}$ . For this value of $\unicode[STIX]{x1D707}$ , $Ca_{c}=0.23$ (see § 3). The p.d.f.s are artificially translated in order to render their power-law behaviours more easily visible.

It should be noted, however, that because of the breakups the total number of drops decays in time. The fraction of drops surviving at time $t$ , $N(t)/N(0)\equiv \int p(m_{1}^{2},t)\,\text{d}m_{1}^{2}$ , indeed decreases exponentially, the decay rate growing rapidly when $\mathit{Ca}$ exceeds $\mathit{Ca}_{c}$ (figure 2 a). Accordingly, the statistics of drop sizes is not stationary and the p.d.f. of $m_{1}^{2}$ varies in time (see figure 2 b, where the p.d.f.s are translated vertically in order to facilitate the comparison at different times). At long times $p(m_{1}^{2},t)$ reaches an asymptotic shape, but this does not show any definite power-law behaviour. The power law observed by Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) is thus recovered only when the time-integrated p.d.f is considered; indeed the distributions shown in Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) were obtained by averaging over both the Lagrangian trajectories and time.

Figure 2. (a) Fraction of surviving drops as a function of time for $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D6FC}=10^{3}$ and different values of $\mathit{Ca}$ . Time is rescaled by the Kolmogorov time ( $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ) of the flow. The inset shows the exponential decay rate of the fraction of surviving drops versus the capillary number rescaled by its critical value. (b) Time-dependent p.d.f. of $m_{1}^{2}$ for $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D6FC}=10^{3}$ , $\mathit{Ca}=1$ and increasing time instants. In the legend, the fraction of drops surviving at time $t$ is also indicated. The red curve is the time-integrated p.d.f. ${\mathcal{P}}(m_{1}^{2})$ corresponding to the same parameters. For the sake of comparison, the p.d.f.s are translated vertically.

Since the dynamics of drops is not statistically stationary, ${\mathcal{P}}(m_{1}^{2})$ may depend on the initial shape of drops, namely on the value of the aspect ratio at time $t=0$ . We thus performed simulations in which the initial shape tensor is $\unicode[STIX]{x1D648}(0)=\operatorname{diag}(\unicode[STIX]{x1D70C}_{0},1,\unicode[STIX]{x1D70C}_{0}^{-1})$ , where $\unicode[STIX]{x1D70C}_{0}>1$ is both the aspect ratio and the largest eigenvalue of $\unicode[STIX]{x1D648}$ at $t=0$ . Two different behaviours are observed depending on the value of $\mathit{Ca}$ . For small $\mathit{Ca}$ , the shape of ${\mathcal{P}}(m_{1}^{2})$ is not affected significantly by the value of $\unicode[STIX]{x1D70C}_{0}$ (not shown). By contrast, for large $\mathit{Ca}$ , the interval over which ${\mathcal{P}}(m_{1}^{2})\sim (m_{1}^{2})^{-1}$ shrinks as $\unicode[STIX]{x1D70C}_{0}$ is increased and the drop volume is kept constant. In this case, indeed, ${\mathcal{P}}(m_{1}^{2})\sim (m_{1}^{2})^{-1}$ only for $m_{1}^{2}\gg \unicode[STIX]{x1D70C}_{0}$ (figure 3 a). The $(m_{1}^{2})^{-1}$ behaviour may therefore be difficult to detect when $\unicode[STIX]{x1D70C}_{0}$ approaches the critical aspect ratio for breakup. In fact, when $\unicode[STIX]{x1D70C}_{0}$ is sufficiently large a second power law emerges for $m_{1}^{2}\ll \unicode[STIX]{x1D70C}_{0}$ whose slope increases as a function of $\mathit{Ca}$ and can turn from negative to positive at large $\mathit{Ca}$ (figure 3 b).

Figure 3. (a) Time-integrated p.d.f. of the largest eigenvalue of $\unicode[STIX]{x1D648}$ for $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D6FC}=10^{3}$ , $\mathit{Ca}=1$ and different values of $\unicode[STIX]{x1D70C}_{0}$ . The p.d.f.s are normalized to 1 to facilitate the comparison. The dashed line is proportional to $(m_{1}^{2})^{-1}$ . (b) Time-integrated p.d.f. of the largest eigenvalue of $\unicode[STIX]{x1D648}$ for $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D70C}_{0}=10^{3}$ and different values of $\mathit{Ca}$ .

The dependence of the deformation and breakup statistics on $\unicode[STIX]{x1D707}$ is shown in figure 4. For small values of $\mathit{Ca}$ , the slope of ${\mathcal{P}}(m_{1}^{2})$ varies with $\unicode[STIX]{x1D707}$ and is steeper for larger viscosity ratios (figure 4 a). It saturates to $-1$ beyond the critical capillary number, but the transition to the supercritical regime is slower for larger values of $\unicode[STIX]{x1D707}$ . These results differ somewhat from those of Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014). The discrepancy may be explained by considering the time scales associated with the breakup process. Whereas the time-integrated statistics displays a weak dependence on the viscosity ratio, the time scale over which breakup occurs depends strongly on $\unicode[STIX]{x1D707}$ , and the breakup process considerably slows down as $\unicode[STIX]{x1D707}$ increases (figure 4 b). For large values of $\unicode[STIX]{x1D707}$ , very long Lagrangian trajectories therefore need to be considered in order to compute ${\mathcal{P}}(m_{1}^{2})$ ; otherwise small deformations are privileged and the slope of ${\mathcal{P}}(m_{1}^{2})$ may be steeper than it actually should be. This point is elucidated further in § 3.

Figure 4. (a) Time-integrated p.d.f. of the largest eigenvalue of $\unicode[STIX]{x1D648}$ for (from bottom to top) $\mathit{Ca}=0.21,0.32,0.51$ and different values of $\unicode[STIX]{x1D707}$ . The p.d.f.s corresponding to different values of $\mathit{Ca}$ are translated vertically. (b) Exponential decay rate of the number of surviving drops as a function of $\unicode[STIX]{x1D707}$ for $\mathit{Ca}=0.6$ .

Figure 5. Mean lifetime as a function of the initial aspect ratio for (a) $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D6FC}=10^{2}$ , $\mathit{Ca}<\mathit{Ca}_{c}$ and (b) $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D6FC}=10^{3}$ , $\mathit{Ca}>\mathit{Ca}_{c}$ . In (a) the $\mathit{Ca}=0.206$ and $\mathit{Ca}=0.137$ curves are multiplied by a factor 3 and $1/3$ , respectively, in order to make the three curves distinguishable.

Finally, we consider the mean lifetime of a drop, $\overline{T}(\unicode[STIX]{x1D70C}_{0})$ , i.e. the mean time it takes for a drop of initial aspect ratio $\unicode[STIX]{x1D70C}_{0}$ to break. Two different behaviours are found depending on whether $\mathit{Ca}$ exceeds or not its critical value. The mean lifetime $\overline{T}(\unicode[STIX]{x1D70C}_{0})$ decreases as a power law of $\unicode[STIX]{x1D70C}_{0}$ if $\mathit{Ca}<\mathit{Ca}_{c}$ and logarithmically if $\mathit{Ca}>\mathit{Ca}_{c}$ (figure 5).

The deformation and breakup statistics presented above is derived analytically in the next section.

3 Analytical predictions

For large deformations, the model of Maffettone & Minale (Reference Maffettone and Minale1998) is statistically equivalent to a vector model proposed by Olbricht, Rallison & Leal (Reference Olbricht, Rallison and Leal1982). The assumptions on the drop and on the external fluid are essentially the same, and the semi-major axis $\boldsymbol{r}$ of the drop satisfies the equation

(3.1) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{r}}=\unicode[STIX]{x1D642}\,\boldsymbol{r}-{\displaystyle \frac{f_{1}(\unicode[STIX]{x1D707})}{2\unicode[STIX]{x1D70F}}}\,\boldsymbol{r}+\sqrt{\frac{r_{eq}^{2}f_{1}(\unicode[STIX]{x1D707})}{\unicode[STIX]{x1D70F}}}\,\unicode[STIX]{x1D743}(t), & & \displaystyle\end{eqnarray}$$

where $r_{eq}$ is the drop equilibrium size and $\unicode[STIX]{x1D743}(t)$ is white noise describing thermal fluctuations. Although thermal noise does not appear in the original model of Olbricht, Rallison & Leal (Reference Olbricht, Rallison and Leal1982), it is included in (3.1) in order to regularize the p.d.f. of $\boldsymbol{r}$ at $r=0$ . It is in any case expected to play a minor role when the flow is turbulent or when deformations larger than $r_{eq}$ are considered. In this linear model, the condition for drop breakup is expressed in terms of $r=|\boldsymbol{r}|$ , i.e. it is assumed that a drop breaks if $r$ exceeds a threshold size $\ell$ .

Equation (3.1) is closely related to (2.1). Indeed, from the vector $\boldsymbol{r}$ one can form the second-rank tensor $\unicode[STIX]{x1D648}=\langle \boldsymbol{r}\otimes \boldsymbol{r}\rangle _{\unicode[STIX]{x1D709}}$ ( $\langle \cdot \rangle _{\unicode[STIX]{x1D709}}$ denotes the average over the realizations of  $\unicode[STIX]{x1D743}(t)$ ), which evolves according to the equation (Olbricht, Rallison & Leal Reference Olbricht, Rallison and Leal1982)

(3.2) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D648}}=\unicode[STIX]{x1D642}\unicode[STIX]{x1D648}+\unicode[STIX]{x1D648}\unicode[STIX]{x1D642}^{\top }-{\displaystyle \frac{f_{1}(\unicode[STIX]{x1D707})}{\unicode[STIX]{x1D70F}}}[\unicode[STIX]{x1D648}-r_{eq}^{2}\unicode[STIX]{x1D644}]. & & \displaystyle\end{eqnarray}$$

The only difference between (2.1) and (3.2) is in the coefficient of the identity, which in (2.1) preserves the volume of the drop whereas it does not enjoy this property in (3.2). Notwithstanding, this term is negligible in both models when the drops are sufficiently deformed. Moreover, when $r$ is large $\unicode[STIX]{x1D648}\approx \boldsymbol{r}\otimes \boldsymbol{r}$ , so $r^{2}$ is the largest eigenvalue of $\unicode[STIX]{x1D648}$ and $\boldsymbol{r}$ the associate eigenvector. The statistics of large drop deformations can therefore be deduced from (3.1), and potential discrepancies between the two approaches are only expected for small deformations (see Vincenzi et al. (Reference Vincenzi, Perlekar, Biferale and Toschi2015) for a more detailed discussion of this point in the $\unicode[STIX]{x1D707}=1$ case).

Let us introduce the Kubo number $\mathit{Ku}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}_{c}$ , where $\unicode[STIX]{x1D70F}_{c}$ is the correlation time of $\unicode[STIX]{x1D735}\boldsymbol{u}$ . In three-dimensional homogeneous and isotropic turbulence $\mathit{Ku}\approx 0.6$ (Girimaji & Pope Reference Girimaji and Pope1990; Bec et al. Reference Bec, Biferale, Boffetta, Cencini, Lanotte, Musacchio and Toschi2006; Watanabe & Gotoh Reference Watanabe and Gotoh2010). However, for $\unicode[STIX]{x1D707}=1$ , it was shown in Musacchio & Vincenzi (Reference Musacchio and Vincenzi2011) that as long as $\mathit{Ku}\lesssim 1$ , the p.d.f. of $r$ does not depend upon $\mathit{Ku}$ appreciably. Furthermore, in Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) the qualitative features of the statistics of drop deformation were found not to depend significantly on the intermittency of the turbulent flow. To make analytical progress, we therefore study (3.1) under the assumption that $\unicode[STIX]{x1D735}\boldsymbol{u}$ has a Gaussian statistics and $\unicode[STIX]{x1D70F}_{c}$ vanishes. More specifically, we assume that $\unicode[STIX]{x1D734}$ and $\unicode[STIX]{x1D64E}$ are zero-mean Gaussian processes with correlations $\langle \unicode[STIX]{x1D6FA}_{ij}(t)\unicode[STIX]{x1D6FA}_{pq}(t^{\prime })\rangle =(d+2)C(\unicode[STIX]{x1D6FF}_{ip}\unicode[STIX]{x1D6FF}_{jq}-\unicode[STIX]{x1D6FF}_{jp}\unicode[STIX]{x1D6FF}_{iq})\unicode[STIX]{x1D6FF}(t-t^{\prime })$ and $\langle S_{ij}(t)S_{pq}(t^{\prime })\rangle =dC(\unicode[STIX]{x1D6FF}_{ip}\unicode[STIX]{x1D6FF}_{jq}+\unicode[STIX]{x1D6FF}_{iq}\unicode[STIX]{x1D6FF}_{jp}-2\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}_{pq}/d)\unicode[STIX]{x1D6FF}(t-t^{\prime })$ , where $d$ is the spatial dimension of the flow and $C>0$ determines the amplitude of the fluctuations of $\unicode[STIX]{x1D735}\boldsymbol{u}$ . In this setting, $\unicode[STIX]{x1D642}(t)$ is a multiplicative noise and is interpreted in the Stratonovich sense (Falkovich, Gawȩdzki & Vergassola Reference Falkovich, Gawȩdzki and Vergassola2001). The form of the correlations ensures that the flow is incompressible and statistically isotropic (e.g. Brunk & Koch Reference Brunk and Koch1998). In addition, the Lyapunov exponent of this flow is $\unicode[STIX]{x1D706}=Cd(d-1)$ (Le Jan Reference Le Jan1984, Reference Le Jan1985).

Owing to statistical isotropy, at long times the p.d.f. of $r$ , $p(r,t)$ , satisfies the Fokker–Planck equation

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}p=-\unicode[STIX]{x2202}_{r}(D_{1}p)+\unicode[STIX]{x2202}_{r}^{2}(D_{2}p), & & \displaystyle\end{eqnarray}$$

where time has been rescaled by $2\unicode[STIX]{x1D70F}/f_{1}(\unicode[STIX]{x1D707})$ (with a slight abuse of notation we continue to denote the rescaled time by $t$ ) and

(3.4a,b ) $$\begin{eqnarray}\displaystyle D_{1}(r)=\left[{\displaystyle \frac{2(d+1)\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}}{d}}-1\right]r+(d-1){\displaystyle \frac{r_{eq}^{2}}{r}},\quad D_{2}(r)={\displaystyle \frac{2\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}}{d}}\,r^{2}+r_{eq}^{2} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})=f_{2}(\unicode[STIX]{x1D707})/f_{1}(\unicode[STIX]{x1D707})$ . This equation can be obtained from the $\unicode[STIX]{x1D707}=1$ case (see Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2005) by noting that the vorticity tensor does not contribute to the time evolution of $p(r,t)$ . The assumptions that $r$ is a positive quantity and drops break at $r=\ell$ are implemented by imposing a reflecting boundary condition at $r=0$ ( $D_{1}p-\unicode[STIX]{x2202}_{r}(D_{2}p)=0$ at $r=0$ ) and an absorbing boundary condition at $r=\ell$ ( $p(\ell ,t)=0$ ), respectively.

The form of the coefficients $D_{1}(r)$ and $D_{2}(r)$ shows that changing the viscosity ratio merely rescales $\mathit{Ca}$ by a factor of $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ . Also note that $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ depends weakly upon $\unicode[STIX]{x1D707}$ , since it varies from $\unicode[STIX]{x1D6FE}(0)=2$ to $\unicode[STIX]{x1D6FE}(\infty )=19/8=2.375$ .

3.1 Time-integrated distribution of drop sizes

Equation (3.3) can be used to derive the power-law behaviour of the time-integrated p.d.f ${\mathcal{P}}(r)$ as well as its dependence on the initial drop-size distribution. From (3.3), ${\mathcal{P}}(r)$ satisfies

(3.5) $$\begin{eqnarray}\displaystyle -{\displaystyle \frac{\text{d}}{\text{d}r}}(D_{1}{\mathcal{P}})+{\displaystyle \frac{\text{d}^{2}}{\text{d}r^{2}}}(D_{2}{\mathcal{P}})=-p(r,0), & & \displaystyle\end{eqnarray}$$

with boundary conditions $D_{1}{\mathcal{P}}-\unicode[STIX]{x2202}_{r}(D_{2}{\mathcal{P}})=0$ at $r=0$ and ${\mathcal{P}}(\ell )=0$ . To obtain (3.5), we have used the fact that in the presence of an absorbing boundary $\lim _{t\rightarrow \infty }p(r,t)=0$ for all $r$ . It is now assumed that $p(r,0)=\unicode[STIX]{x1D6FF}(r-r_{0})$ with $r_{0}>r_{eq}$ , i.e. a monodisperse initial distribution. Integrating (3.5) from 0 to $r$ and using the reflecting boundary condition at $r=0$ yields

(3.6a,b ) $$\begin{eqnarray}\displaystyle -D_{1}(r){\mathcal{P}}(r)+{\displaystyle \frac{\text{d}}{\text{d}r}}[D_{2}(r){\mathcal{P}}(r)]=\left\{\begin{array}{@{}ll@{}}0 & \quad \text{if }0\leqslant r<r_{0},\\ -1 & \quad \text{if }r_{0}<r\leqslant \ell .\end{array}\right. & & \displaystyle\end{eqnarray}$$

The solution of (3.6) takes the form (Risken Reference Risken1989)

(3.7) $$\begin{eqnarray}\displaystyle {\mathcal{P}}(r)\propto \left\{\begin{array}{@{}ll@{}}\text{e}^{-\unicode[STIX]{x1D6F7}(r)}[\unicode[STIX]{x1D711}(\ell )-\unicode[STIX]{x1D711}(r_{0})]\quad & \text{if}~0\leqslant r\leqslant r_{0},\\ \text{e}^{-\unicode[STIX]{x1D6F7}(r)}[\unicode[STIX]{x1D711}(\ell )-\unicode[STIX]{x1D711}(r)]\quad & \text{if}~r_{0}<r\leqslant \ell ,\end{array}\right. & & \displaystyle\end{eqnarray}$$

with

(3.8a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}(r)=\ln D_{2}(r)-\int _{r_{1}}^{r}\frac{D_{1}(\unicode[STIX]{x1D701})}{D_{2}(\unicode[STIX]{x1D701})}\,\text{d}\unicode[STIX]{x1D701},\quad \unicode[STIX]{x1D711}(r)=\int _{r_{1}}^{r}{\displaystyle \frac{\text{e}^{\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D701})}}{D_{2}(\unicode[STIX]{x1D701})}}\,\text{d}\unicode[STIX]{x1D701}, & & \displaystyle\end{eqnarray}$$

where the specific value of $r_{1}$ is irrelevant. To examine the behaviour of ${\mathcal{P}}(r)$ for $r_{eq}\ll r\ll \ell$ , we now insert the limiting forms of $D_{1}(r)$ and $D_{2}(r)$ for $r_{eq}\rightarrow 0$ into (3.8) and obtain $\text{e}^{\unicode[STIX]{x1D6F7}(r)}\sim r^{\unicode[STIX]{x1D6FD}}$ and $\unicode[STIX]{x1D711}(r)\sim r^{\unicode[STIX]{x1D6FD}-1}$ with $\unicode[STIX]{x1D6FD}=1-d+d/2\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}$ . Therefore, there exists a critical value of the capillary number, $\mathit{Ca}_{c}=1/2\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ , such that for $\mathit{Ca}<\mathit{Ca}_{c}$

(3.9a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{P}}(r)\sim \left\{\begin{array}{@{}ll@{}}(\ell ^{\unicode[STIX]{x1D6FD}-1}-r_{0}^{\unicode[STIX]{x1D6FD}-1})r^{-\unicode[STIX]{x1D6FD}} & \quad \text{if }r_{eq}\ll r\ll r_{0},\\ \ell ^{\unicode[STIX]{x1D6FD}-1}r^{-\unicode[STIX]{x1D6FD}} & \quad \text{if }r_{0}\ll r\ll \ell ,\end{array}\right. & & \displaystyle\end{eqnarray}$$

whereas for $\mathit{Ca}>\mathit{Ca}_{c}$

(3.10a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{P}}(r)\sim \left\{\begin{array}{@{}ll@{}}(\ell ^{\unicode[STIX]{x1D6FD}-1}-r_{0}^{\unicode[STIX]{x1D6FD}-1})r^{-\unicode[STIX]{x1D6FD}} & \quad \text{if}~r_{eq}\ll r\ll r_{0},\\ r^{-1} & \quad \text{if}~r_{0}\ll r\ll \ell .\end{array}\right. & & \displaystyle\end{eqnarray}$$

(The exact form of ${\mathcal{P}}(r)$ over the entire interval $0\leqslant r\leqslant \ell$ may be calculated by using the full expressions of $D_{1}(r)$ and $D_{2}(r)$ in (3.4) and involves a hypergeometric function; the details, however, are omitted.) Since $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ depends weakly on $\unicode[STIX]{x1D707}$ , the same holds true for $\mathit{Ca}_{c}$ . The above value of $\mathit{Ca}_{c}$ was also found by Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) for more general flows; they applied a criterion based on the statistics of the finite-time Lyapunov exponents of the flow that was previously used to study the deformation of flexible polymers (Balkovsky, Fouxon & Lebedev Reference Balkovsky, Fouxon and Lebedev2001). Likewise, the prediction of Biferale, Meneveau & Verzicco (Reference Biferale, Meneveau and Verzicco2014) for the exponent $\unicode[STIX]{x1D6FD}$ in the $\mathit{Ca}<\mathit{Ca}_{c}$ case reduces to the expression above when $\unicode[STIX]{x1D735}\boldsymbol{u}$ has the properties considered here.

The scaling of ${\mathcal{P}}(r^{2})$ can be deduced from that of ${\mathcal{P}}(r)$ by using ${\mathcal{P}}(r^{2})=(1/2)r^{-1}{\mathcal{P}}(r)$ . The above power-law behaviours thus reproduce the numerical results shown in figures 1 and 3 for the time-integrated p.d.f.s of the squared length of the semi-major axis. It should be noted that whereas the power-law behaviour of ${\mathcal{P}}(r)$ for small $r$ is specific to a monodisperse initial distribution, the large- $r$ power law holds for any $p(r,0)$ that vanishes beyond a given $r_{\star }<\ell$ . Integrating (3.5) from 0 to $r>r_{\star }$ indeed yields (3.6b ) and hence (3.9b ) or (3.10b ) depending on the value of $\mathit{Ca}$ . If, by contrast, the initial size of the drops can approach $\ell$ , in general ${\mathcal{P}}(r)$ does not display a power-law behaviour.

3.2 Time-dependent distribution of drop sizes and breakup frequency

The eigenfunctions of the Fokker–Planck operator that satisfy the reflecting boundary condition at $r=0$ are of the form $f_{\unicode[STIX]{x1D708}}(r)=r^{d-1}\text{}_{2}F_{1}(c_{\unicode[STIX]{x1D708}}^{-},c_{\unicode[STIX]{x1D708}}^{+},d/2,-\unicode[STIX]{x1D716}r^{2})$ (Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2005), where $\text{}_{2}F_{1}$ is the Gauss hypergeometric function with $\unicode[STIX]{x1D716}=2\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}/dr_{eq}^{2}$ and

(3.11) $$\begin{eqnarray}\displaystyle c_{\unicode[STIX]{x1D708}}^{\pm }={\displaystyle \frac{d}{4}}\left[{\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}}}+1\right]\pm {\displaystyle \frac{1}{4}}\sqrt{d^{2}\left[{\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}}}-1\right]^{2}-{\displaystyle \frac{2d\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})\mathit{Ca}}}}. & & \displaystyle\end{eqnarray}$$

The absorbing boundary condition $f_{\unicode[STIX]{x1D708}}(\ell )=0$ selects a discrete set of acceptable eigenfunctions. The p.d.f. of $r$ can thus be expanded as $p(r,t)=\sum _{n=1}^{\infty }a_{n}\text{e}^{-\unicode[STIX]{x1D708}_{n}t}f_{\unicode[STIX]{x1D708}_{n}}(r)$ , and hence

(3.12) $$\begin{eqnarray}\displaystyle p(r,t)\sim \text{e}^{-\unicode[STIX]{x1D708}_{1}t}f_{\unicode[STIX]{x1D708}_{1}}(r)\quad \text{as}~t\rightarrow \infty . & & \displaystyle\end{eqnarray}$$

This result confirms that at long times $p(r,t)$ approaches an asymptotic shape, but this does not show a power-law behaviour (figure 2 a). From (3.12), the fraction of drops surviving at time $t$ decays as

(3.13) $$\begin{eqnarray}\displaystyle N(t)/N(0)\equiv \int _{0}^{\ell }p(r,t)\,\text{d}r\sim \text{e}^{-\unicode[STIX]{x1D708}_{1}t}\quad \text{as}~t\rightarrow \infty , & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{1}$ is the smallest solution of the equation $f_{\unicode[STIX]{x1D708}_{n}}(\ell )=0$ . Figure 6(a) shows that the decay rate of the drop number increases rapidly as a function of $\mathit{Ca}$ when $\mathit{Ca}$ exceeds its critical value, whereas it decreases as a function of $\unicode[STIX]{x1D707}$ . In addition, although $\mathit{Ca}_{c}$ depends weakly on $\unicode[STIX]{x1D707}$ , the transition to the supercritical regime is much steeper at small $\unicode[STIX]{x1D707}$ (see also figure 6(b)). These results reproduce the behaviours observed in the numerical simulations (figures 2(b) and 4).

Figure 6. Exponential decay rate of the number of drops for $d=3$ and $\ell =10^{3}$ (a) as a function of $\unicode[STIX]{x1D707}$ and $\mathit{Ca}$ and (b) as a function of the capillary number rescaled by its critical value for fixed $\unicode[STIX]{x1D707}=0.1,1,10$ (from top to bottom).

3.3 Mean lifetime of a drop

The average time it takes for a drop of initial size $r_{0}$ to break can be calculated from ${\mathcal{P}}(r)$ as follows. Consider the transition probability $p(r,t|r_{0},0)$ , which is the solution of (3.3) that satisfies the initial condition $p(r,0|r_{0},0)=\unicode[STIX]{x1D6FF}(r-r_{0})$ . Let $T(r_{0})$ be the time it takes for the drop to break in a given realization of the flow and of thermal noise, and let $\mathbb{P}(r_{0},t)$ be the probability of $T(r_{0})$ taking the value $t$ . Note that $\mathbb{P}(r_{0},t)=-\unicode[STIX]{x2202}_{t}F$ , where $F(r_{0},t)=\int _{t}^{\infty }\mathbb{P}(r_{0},s)\,\text{d}s$ is the probability that $T(r_{0})\geqslant t$ and can be written as $F(r_{0},t)=\int _{0}^{\ell }p(r,t|r_{0},0)\,\text{d}r$ . Therefore, the average of $T(r_{0})$ is (Gardiner Reference Gardiner1983)

(3.14) $$\begin{eqnarray}\displaystyle \overline{T}(r_{0})=\int _{0}^{\infty }t\,\mathbb{P}(r_{0},t)\,\text{d}t=-\int _{0}^{\infty }t[\unicode[STIX]{x2202}_{t}F(r_{0},t)]\,\text{d}t=\int _{0}^{\infty }F(r_{0},t)\,\text{d}t, & & \displaystyle\end{eqnarray}$$

where we used $\lim _{t\rightarrow \infty }F(r_{0},t)=0$ , a consequence of the absorbing boundary condition for $p(r,t|r_{0},0)$ . By changing the order of integration, we finally obtain

(3.15) $$\begin{eqnarray}\displaystyle \overline{T}(r_{0})=\int _{0}^{\ell }{\mathcal{P}}(r)\,\text{d}r, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{P}}(r)$ is the solution of (3.5) corresponding to the initial condition $p(r,0)=\unicode[STIX]{x1D6FF}(r-r_{0})$ . Inserting now the asymptotic behaviours (3.9b ) and (3.10b ) into (3.15) yields

(3.16) $$\begin{eqnarray}\displaystyle \overline{T}(r_{0})\sim \left\{\begin{array}{@{}ll@{}}\left({\displaystyle \frac{\ell }{r_{eq}}}\right)^{\unicode[STIX]{x1D6FD}-1}-\left({\displaystyle \frac{r_{0}}{r_{eq}}}\right)^{\unicode[STIX]{x1D6FD}-1} & \text{if}~\mathit{Ca}<\mathit{Ca}_{c},\\[12.0pt] \ln \bigg({\displaystyle \frac{\ell }{r_{0}}}\bigg) & \text{if}~\mathit{Ca}>\mathit{Ca}_{c},\end{array}\right. & & \displaystyle\end{eqnarray}$$

as seen in figure 5.

4 Improved version of the model of Maffettone & Minale (Reference Maffettone and Minale1998)

Figure 7. (a) Dependence of the critical capillary number on the viscosity ratio in the original model (red solid curve) and in the improved model (blue dashed curve). (b) Exponential decay rate of the number of drops for $d=3$ and $\ell =10^{3}$ as a function of $\mathit{Ca}$ and $\unicode[STIX]{x1D707}$ for the improved model.

Maffettone & Minale (Reference Maffettone and Minale1998) proposed a modification of their model that improves the agreement with experimental data for high viscosity ratios and large capillary numbers. In the modified model, the coefficient in front of the strain tensor also depends on $\mathit{Ca}$ , i.e. $f_{2}(\unicode[STIX]{x1D707})$ is replaced with

(4.1) $$\begin{eqnarray}\displaystyle \tilde{f}_{2}(\unicode[STIX]{x1D707},\mathit{Ca})=f_{2}(\unicode[STIX]{x1D707})+{\displaystyle \frac{3(\unicode[STIX]{x1D70E}\mathit{Ca})^{2}}{2+6(\unicode[STIX]{x1D70E}\mathit{Ca})^{2}}}, & & \displaystyle\end{eqnarray}$$

where the coefficient $\unicode[STIX]{x1D70E}$ accounts for the fact that here $\mathit{Ca}$ is defined in terms of $\unicode[STIX]{x1D706}$ (hence in our case $\unicode[STIX]{x1D70E}=1.72$ ). The above expression still reproduces the theoretical limits, for both small $\mathit{Ca}$ and large $\unicode[STIX]{x1D707}$ , as well as the affine deformation of the drop when $\unicode[STIX]{x1D707}=1$ and $\mathit{Ca}\rightarrow \infty$ . Using $\tilde{f}_{2}$ instead of $f_{2}$ yields significantly more accurate predictions for a pure shear; for an elongational flow, the effect is much weaker (Maffettone & Minale Reference Maffettone and Minale1998).

It should be noted that the original model and the improved one can be mapped into each other by suitably modifying the viscosity ratio and the capillary relaxation time. The original model with parameters $\unicode[STIX]{x1D707}$ , $\unicode[STIX]{x1D70F}$ is indeed the same as the improved one with parameters $\unicode[STIX]{x1D707}^{\prime }$ , $\unicode[STIX]{x1D70F}^{\prime }$ , where $\unicode[STIX]{x1D707}^{\prime }$ and $\unicode[STIX]{x1D70F}^{\prime }$ are the solutions of the system

(4.2a,b ) $$\begin{eqnarray}\displaystyle f_{1}(\unicode[STIX]{x1D707}^{\prime })/\unicode[STIX]{x1D70F}^{\prime }=f_{1}(\unicode[STIX]{x1D707})/\unicode[STIX]{x1D70F},\quad \tilde{f}_{2}(\unicode[STIX]{x1D707}^{\prime },\unicode[STIX]{x1D70E}\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}^{\prime })=f_{2}(\unicode[STIX]{x1D707}). & & \displaystyle\end{eqnarray}$$

Therefore, for fixed $\unicode[STIX]{x1D707}$ and $\mathit{Ca}$ , the results described in the previous sections also hold for the improved model, provided the parameters are suitably adjusted. The reader should note that such a nonlinear transformation of the parameters leads to a slight variation, quantitatively, of the results without changing the overall picture. It is nonetheless important to examine the effect of the modified coefficient $\tilde{f}_{2}$ on quantities such as the critical capillary number, the rate of decay of the drop fraction and the exponent $\unicode[STIX]{x1D6FD}$ that defines the power-law behaviour of both the p.d.f. of the size and the mean lifetime. This is achieved by replacing $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ in § 3 with $\tilde{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D707},\mathit{Ca})=\tilde{f}_{2}(\unicode[STIX]{x1D707},\mathit{Ca})/f_{1}(\unicode[STIX]{x1D707})$ . Thus, the differences between the two versions of the model are mainly due to the fact that $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D707})$ varies weakly with $\unicode[STIX]{x1D707}$ and is bounded for $\unicode[STIX]{x1D707}\rightarrow \infty$ , whereas $\tilde{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D707},\mathit{Ca})\rightarrow \infty$ in the same limit. It is shown below that, for a turbulent flow, these differences impact our predictions only marginally.

When the improved model is considered, the critical capillary number is the solution of the cubic equation $\tilde{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D707},\mathit{Ca}_{c})\mathit{Ca}_{c}=1/2$ . It can be checked that the discriminant of this equation is negative for all values of $\unicode[STIX]{x1D707}$ and hence there is only one real root. Figure 7(a) compares the critical capillary number in the original model and in the improved one. In both cases, $\mathit{Ca}_{c}$ is a decreasing function of $\unicode[STIX]{x1D707}$ . The main difference is that, for $\unicode[STIX]{x1D707}\rightarrow \infty$ , $\mathit{Ca}_{c}$ tends to the asymptotic value $4/19\approx 0.21$ in the original model, whereas it tends to zero in the improved one. Nevertheless, for the values of $\unicode[STIX]{x1D707}$ typically found in experiments, $\mathit{Ca}_{c}$ does not differ considerably in the two models.

In the improved model, the rate of decay of the number of drops is slightly greater and decreases less rapidly as a function of $\unicode[STIX]{x1D707}$ (compare figures 6 a and 7 b). The exponent $\unicode[STIX]{x1D6FD}$ takes the form $\unicode[STIX]{x1D6FD}=1-d+d/2\tilde{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D707},\mathit{Ca})\mathit{Ca}$ ; it is smaller than in the original model and varies more with $\unicode[STIX]{x1D707}$ , as a consequence of the different behaviour of $\tilde{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D707},\mathit{Ca})$ (figure 8). Both for $\unicode[STIX]{x1D6FD}$ and the decay rate, the discrepancies between the two models are, however, small.

In conclusion, despite some quantitative differences, for realistic values of $\unicode[STIX]{x1D707}$ and Ca the qualitative behaviour of the model of Maffettone & Minale (Reference Maffettone and Minale1998) in a turbulent flow is largely insensitive to the use of either $f_{2}(\unicode[STIX]{x1D707})$ or $\tilde{f}_{2}(\unicode[STIX]{x1D707},\mathit{Ca})$ .

Figure 8. Exponent $\unicode[STIX]{x1D6FD}$ as a function of the capillary number for $d=3$ (a) in the original model and (b) in the improved model.