Role of viscosity in turbulent drop break-up

Abstract

We investigate drop break-up morphology, occurrence, time and size distribution, through large ensembles of high-fidelity direct-numerical simulations of drops in homogeneous isotropic turbulence, spanning a wide range of parameters in terms of the Weber number $We$, viscosity ratio between the drop and the carrier flow $\mu _r=\mu _d/\mu _l$, where d is the drop diameter, and Reynolds ($Re$) number. For $\mu _r \leq 20$, we find a nearly constant critical $We$, while it increases with $\mu _r$ (and $Re$) when $\mu _r > 20$, and the transition can be described in terms of a drop Reynolds number. The break-up time is delayed when $\mu _r$ increases and is a function of distance to criticality. The first break-up child-size distributions for $\mu _r \leq 20$ transition from M to U shape when the distance to criticality is increased. At high $\mu _r$, the shape of the distribution is modified. The first break-up child-size distribution gives only limited information on the fragmentation dynamics, as the subsequent break-up sequence is controlled by the drop geometry and viscosity. At high $We$, a $d^{-3/2}$ size distribution is observed for $\mu _r \leq 20$, which can be explained by capillary-driven processes, while for $\mu _r > 20$, almost all drops formed by the fragmentation process are at the smallest scale, controlled by the diameter of the very extended filament, which exhibits a snake-like shape prior to break-up.

1. Introduction

1.1. Context and motivation

Breakage of dispersed liquid drops in a carrier flow is of importance in engineering, from industrial manufacturing of petrochemicals, polymers or pharmaceuticals, to transport of mixtures of oil, water and gas in pipelines (Galinat et al. 2005), in the environment with oil-spill mitigation strategies (Li et al. 2017); while gas bubbles breaking in turbulence controls mass exchange at the air–water interface, from lakes to the ocean (Deike 2022). Quantifying mass exchange in such flows requires a model for drop break-up in turbulence over a wide range of parameters. Drops that disperse in a turbulent carrier flow typically lead to a distribution of drop diameters from $\mathrm {\mu }$m to cm for oil-in-water problems.

In homogeneous turbulence, the transient evolution of the drop-size distribution is described by a population balance equation, whose source terms account for a set of physical phenomena, in particular, the rate of consumption of drops due to break-up, and the rate of production of drops due to the break-up of larger drops (Martinez-Bazan et al. 2010; Qi, Masuk & Ni 2020). The break-up kernel then determines the rate of production of drops, from one size to another, and can be decomposed into a parent break-up frequency and child-size distribution (Martinez-Bazan et al. 2010; Aiyer et al. 2019; Gaylo, Hendrickson & Yue 2021; Ruth et al. 2022). A large variety of break-up mechanisms have been introduced (Balachandar & Eaton 2010), with most models for droplet break-up in a turbulent flow considering a single mechanism with drop interaction with a single turbulent eddy, or a continuum of turbulence scales, with fluctuations at the size of the drop being the most effective part of the process (Hinze 1955; Luo & Svendsen 1996; Liao & Lucas 2009; Lalanne, Masbernat & Risso 2019). Binary break-up is often assumed, while highly viscous or energetic break-up presents a rapid sequence with multiple time scales and processes (Roccon et al. 2017; Rivière et al. 2022; Ruth et al. 2022). Several models consider dimensional units (Martinez-Bazan et al. 2010; Qi et al. 2020), complicating the comparisons between different configurations and missing some hidden controlling parameters, while empirical formula, with important assumptions in the models not rigorously verified, have been applied outside their initial fitting range (Tsouris & Tavlarides 1994; Andersson & Andersson 2006; Martinez-Bazan et al. 2010). Recent computational work has investigated various aspects of the break-up of drops in turbulent flow, including long bubble lifetime near critical conditions and limits of the Hinze-scale concept (Vela-Martín & Avila 2022), the interaction of turbulence cascade and the break-up dynamics in steady-state turbulent emulsion configuration (Crialesi-Esposito et al. 2022; Crialesi-Esposito, Chibbaro & Brandt 2023), while experimental work has discussed the role of small-scale eddies in bubble break-up (Qi et al. 2022) or drop break-up in a von Kármán flow (Ravichandar et al. 2022). Separately, atomization of jets has received considerable attention, with analysis of the subsequent instabilities controlling the scale selection, break-up and drop-size distribution (e.g. Eggers & Villermaux 2008; Villermaux 2020) or analysis of the flow vortical structure (e.g. Zandian, Sirignano & Hussain 2018).

At this stage, no universal formulation exists over a wide range of controlling non-dimensional numbers, comparing inertial, capillary and viscous effects, and spanning various viscosity ratios and turbulence intensities. When the viscosity ratio is significantly larger than one, surface tension is less important and the break-up process can take much longer than the turbulence time scales, with viscous breakage only occurring when the drop has been substantially stretched, similar to low-$Re$ break-up (Stone 1994; Andersson & Andersson 2006; Roccon et al. 2017). We aim to describe the occurrence of break-up, the break-up time and the number and child sizes, as a function of a set of well-chosen non-dimensional numbers defined as a function of the liquids properties.

1.2. Non-dimensional numbers

Considering an initially spherical drop of diameter $d$, viscosity $\mu _d$, and density $\rho _d$ in a turbulent carrier flow, characterized by velocity $U$, density $\rho _l$, viscosity $\mu _l$ and surface tension $\sigma$, we have seven controlling parameters, with three units, leading to four non-dimensional groups. Early theoretical discussion on drop break-up in turbulent flow (Hinze 1955) has mostly been based on the competition between inertial and surface forces, and as such the most important relevant dimensionless quantity is the Weber number

(1.1)\begin{equation} We = {\rho_d U^2 d}/{\sigma}={2\rho_d \varepsilon^{2/3}d^{5/3}}/{\sigma}, \end{equation}

where $\varepsilon$ is the turbulence dissipation rate, used to characterize the inertial range of the turbulent flow. When considering a turbulent flow, the characteristic velocity is given by the root mean square fluctuations, $u_{rms}$, the energy-injection scale is the integral length scale, $L_{int}$, and the turbulence intensity is characterized by the turbulent dissipation rate, $\varepsilon$, with $\varepsilon = C_{\varepsilon }u_{rms}^3/L_{int}$ where $C_{\varepsilon } \approx 2$ is a non-dimensional number universal at high Reynolds number (Dimotakis 2005; Vassilicos 2015). The turbulence is characterized by the Reynolds number, which can be defined at the integral length scale: $Re = {u_{{rms}}L_{int}}/{\nu _l} \propto {k^2}/{\nu _l \varepsilon }$, with $\nu _l=\mu _l/\rho _l$ the kinematic viscosity and $k \propto u_{{rms}}^2$ the turbulent kinetic energy. Another commonly used characterization of the turbulent flow is the Taylor-scale Reynolds number, which compares the strength of velocity fluctuations within the inertial range

(1.2)\begin{equation} Re_{\lambda}= {u_{rms} \lambda}/{\nu_l}, \end{equation}

where $\lambda$ is the Taylor microscale: $\lambda = \sqrt {{15\nu _l}/{\varepsilon }}u_{rms}=\sqrt {15\nu _l}({L_{int}}/{C_{\varepsilon }})^{1/3}\varepsilon ^{-1/6}$, such that $Re \propto Re_{\lambda }^2$. For $Re_{\lambda }>100$, the mixing properties become mostly independent of $Re$ (Dimotakis 2005; Mostert, Popinet & Deike 2022).

For drops of low viscosity (or bubbles), break-up occurs above a certain critical Weber number, $We_c^0$, which is of order one (Hinze 1955; Risso & Fabre 1998; Rivière et al. 2022), and the associated Hinze scale can be defined: $d_h^0=({W e_c^0}/{2})^{3 / 5}({\sigma }/{\rho })^{3 / 5} \varepsilon ^{-2 / 5}$. As discussed in Rivière et al. (2021, 2022), the critical Weber number, or Hinze scale, should be considered as a soft limit as the turbulent flow presents large fluctuations, leading to a broad range of break-up times for the same turbulent conditions, especially close to stable conditions (Vela-Martín & Avila 2022), which makes estimations of the critical Weber number challenging. For bubbles and low viscosity drops, variations in the experimentally reported values of the critical Weber typically go from 0.5 to 5 (Risso & Fabre 1998; Andersson & Andersson 2006; Liao & Lucas 2009; Martinez-Bazan et al. 2010; Vejražka, Zedníková & Stanovskỳ 2018). The critical Weber number, or Hinze scale, is defined in a statistical sense, for a given time (and window) of observation, typically corresponding to the conditions where half of the droplets will break. Such a definition inherently depends on an observation time constrained by experimental or computational constraints. The range of critical Weber numbers observed can be related to variability in the experimental configurations with different large-scale shear or spatial variations of the dissipation rate.

When the viscosity of the liquid is increased, viscous effects will come into play and we introduce a dimensionless ‘viscosity group’ (Hinze 1955), which can be defined as the Ohnesorge number: $Oh = {\mu _d}/{\sqrt {\rho _d \sigma d}}$. At low $Oh$ (typically $Oh <0.1$), viscous effects are not needed to describe break-up properties, while they become essential at higher $Oh$. This dimensionless viscosity group can be replaced by another choice of non-dimensional number, such as the viscosity ratio

(1.3)\begin{equation} \mu_r={\mu_d}/{\mu_l}. \end{equation}

Finally, the fourth non-dimensional parameter can be chosen as the density ratio, which is taken to be unity in the present study, $\rho _d/\rho _l \approx 1$, relevant for several oil-in-water problems.

The paper is organized as follows: § 2 presents the computational approach and the large ensemble of simulations performed. Section 3 discusses the critical break-up conditions as a function of the appropriate viscosity group, and characterizes the associated time scale of break-up and size distribution of drops. Section 4 presents conclusions and perspectives for future work.

2. Approach

We perform direct numerical simulations (DNS) of the two-phase, incompressible Navier–Stokes equations with surface tension using the open-source software Basilisk (Popinet 2015, 2018). We solve the full two-phase Navier–Stokes equations based on a quad/octree adaptive spatial discretization and a multilevel Poisson solver. The interface between the two liquids is reconstructed by a sharp geometric volume-of-fluid (VOF) method, momentum-conserving scheme and an accurate, well-balanced surface tension model. The methodology has been demonstrated to be accurate for capturing the multi-scale nature of turbulent multiphase flows involving break-up, leveraging adaptive mesh refinement, where the spatial discretization is adjusted to follow the spatial and temporal evolution of flow structures. This has been amply demonstrated in recent work on breaking waves (Deike, Melville & Popinet 2016; Mostert & Deike 2020; Mostert et al. 2022), bubble bursting (Berny et al. 2020), bubble deformation and break-up in turbulence (Perrard et al. 2021; Rivière et al. 2021, 2022) and bubble gas transfer (Farsoiya, Popinet & Deike 2021; Farsoiya et al. 2023).

The simulations are performed in two steps. First, we perform a precursor simulation to generate a homogeneous, isotropic turbulent flow in a single phase, with well-defined properties. Second, we insert the drop with initially no velocity and analyse its evolution in the forced turbulent flow.

The turbulent flow is generated by a volumetric force that is locally proportional to the velocity at every point of the physical space. This approach has been introduced by Rosales & Meneveau (2005) and produces a well-characterized, homogeneous, isotropic turbulent flow with properties that closely resemble those obtained with a spectral code. This approach has been used to study bubble rising (Loisy & Naso 2017), bubble deformation (Perrard et al. 2021), gas exchange (Farsoiya et al. 2021, 2023) and break-up (Rivière et al. 2021) in turbulence. The turbulence statistics, including kinetic energy, dissipation rate and turbulent Reynolds number $Re_{\lambda }$ as a function of time, have been verified to reach a statistically stationary regime, while the second-order structure functions display a well-defined inertial range, in good agreement with the homogeneous, isotropic turbulence scaling described in the literature; as demonstrated in Perrard et al. (2021), Rivière et al. (2021) and Farsoiya et al. (2021, 2023). Once the statistically stationary regime is reached, we extract the velocity fields at different times, termed as precursors, and utilize them as initial conditions for subsequent simulations where we introduce a drop in the central region. A central drop of diameter $d$ and density $\rho _d$ is placed in the domain, and the flow in the inner phase is initially set to be zero. The solver relaxes within a few time steps after the initialization and the velocity and pressure become unaffected by the specific initialization details well before deformations become significant. Each precursor is spaced by at least one eddy-turnover time at the scale of the drop to ensure statistical independence between the precursor flow. A detailed description of the numerical configuration can be found in Perrard et al. (2021), Rivière et al. (2021) and Farsoiya et al. (2021, 2023).

The maximum level of refinement, $L$, sets the smallest grid size in the domain and allows us to compare with a fixed grid, equivalent to $(2^{L})^3$ number of grid points. Adaptive mesh refinement allows us to reach a wide separation of scales and obtain grid convergence for child-size data. In the simulations in this study, the flow variables of interest are the VOF and velocity fields, as in Rivière et al. (2021) and Mostert et al. (2022), and a full discussion on the adaptive mesh refinement algorithm can be found in van Hooft et al. (2018) and Mostert et al. (2022). Grid convergence is verified in multiple ways. First, the present numerical set-up is inherited from our work on bubble break-up, where grid convergence on break-up time, number of children and size distribution was verified for a wide range of Weber numbers (Rivière et al. 2021). Similarly, all metrics discussed in the present paper (break-up occurrence, time and size distribution of child droplets) are verified to be grid independent in an ensemble average sense, while individual realizations also demonstrate grid convergence (see Appendix A). Adaptive mesh refinement allows us to reduce the number of grid points, reducing memory requirements and run time and allowing us to perform a wide sweep of parameters at moderate cost. The typical number of grid points for L10 is between $3\times 10^6$ and $6\times 10^6$ depending on the values of $We$ and viscosity ratios (to be compared with $(2^{10})^3\approx 10^9$ grid points with a fixed grid, i.e. an ${\approx }200$-fold reduction in the number of grid points). Reaching a higher Reynolds number is feasible but significantly increases the number of grid points to resolve the finer scales of the turbulent flow, hence the computational cost increases. Here, we focus on the wide sweep in viscosity ratio and Weber number.

The ratio of the initial drop radius to the box size is 0.067 so that the volume ratio is 0.0012; and we can consider that we are in dilute regime. The drop size compared with the box size is not varied in the present work, and the drop is always in the inertial range (with slight variations of the ratio of the initial drop size with the turbulence length scales as the Reynolds number is changed; see Farsoiya et al. (2021, 2023)). Because we are in the dilute regime, coalescence will be negligible statistically.

A wide range of parameters relevant to the physical problem is swept (see table 1), with $\mu _r$ between $10^{-2}$ and 350, equivalent to $Oh$ from $10^{-5}$ to 28, in order to assess the role of viscosity in break-up frequency, child-size distribution and break-up morphology. The value of $We$ is varied from 1 to 30, spanning regimes including no break-up, break-up threshold and violent break-up. For each condition, we perform an ensemble of simulations (10 to 100 realizations) in order to obtain statistical convergence to the metrics of interest, following the approach successfully employed for bubble break-up (Perrard et al. 2021; Rivière et al. 2021, 2022). Thanks to the adaptive mesh refinement approach, each individual run has a low computational cost making the very large ensemble possible.

Table 1. List of simulations for various $We$, $Re$ and $\mu _r$ indicating the level of resolution on the interface, the number of runs and the initial drop-to-grid-size ratio. The minimal length of the simulation is $t/t_c = 20$ for cases that do not break, where $t_c$ is the eddy-turnover time at the scale of drop diameter and is given by $t_c=d^{2/3}\varepsilon ^{-1/3}$. The simulation time is large compared with the Kolmogorov time scale $\tau _{\eta }$ with $t/\tau _{\eta } \approx 200$, where $\tau _{\eta }=\sqrt {\nu _l/\epsilon }$. Small ensembles are used to obtain the break-up phase diagram (figures 2 and 3), while large ensembles are used to obtain the child numbers and size distributions (figure 4).

We will define the critical Weber as in Perrard et al. (2021) and Rivière et al. (2021) in terms of the size of our ensembles (ranging from 10 to 150 runs per conditions) and the duration chosen for the simulations (at least 20 eddy-turnover times at the scale of the drop, equivalent to approximately 200 times the Kolmogorov time scale, which is long enough to capture significant statistical variability in the break-up time but is still manageable in terms of computational cost for the large range of parameters we consider) and explore the dependence of the obtained critical Weber number on the viscosity ratio and turbulence Reynolds number. We expect the proposed scalings for the mean lifetime to be independent of the definition while an analysis of the lifetime distribution would be sensitive to the statistical definitions (as in any other studies) used for the critical conditions. We note that all results presented here are statistically converged in terms of the size of the sample (i.e. results do not change if we only consider half of the ensemble per condition).

3. Drop break-up: occurrence, morphology and time

3.1. Drop morphology

Examples of drop break-up in turbulent flow for increasing Weber number and viscosity ratios are shown in figure 1 and illustrate the complexity of the processes at play. The top row shows a low viscosity drop ($\mu _r=0.01$) close to critical break-up conditions ($We=4\gtrsim We_c^0$), with a mild break-up occurring on a time scale controlled by the turbulent flow, and leading to two children of similar sizes. The second row shows the effect of increasing the droplet viscosity, for the same turbulent flow and $We=4\gtrsim We_c^0$. As the viscosity of the drop increases, break-up is delayed and the formation of a thin thread prior to break-up is observed. This leads to an increased number of children with increasing $\mu _r$. Note that when further increasing the viscosity ratio at $We=4$, no break-up is observed within the observation time ($t/t_c=20$).

Figure 1. Break-up event and drop morphology for $Re_{\lambda }=77$, $We = 4$, $\mu _r = 0.01,10$ (a–f), $Re_{\lambda }=77$, $We=8$, $\mu _r = 0.01,150$ (g–l) and $Re_{\lambda }=77$, $We=20$, $\mu _r = 150$ (m–o). Hinze scale $d_h$ is shown, considering a constant $We_c^0$, and indicates the scale of the images. The time $t/t_c$ corresponds to the time since injection of the drop in the turbulent flow. The colour variation is due to the ray tracing of the oil droplet.

The next two rows show the $We=8$ dynamics for increasing viscosity ratios. For low to intermediate viscosity ratios, a scenario somewhat similar to bag break-up is observed, reminiscent of the role of turbulent fluctuations at the drop scale at high $We$ far from critical conditions ($We\gg We_c$). At very high viscosity ($\mu _r=150$) and $We=8$, which is close to break-up threshold, we observe a very thin thread that characterizes the break-up, with a break-up time more than an order of magnitude longer due to the viscous effects. Finally, at $\mu _r=150$ and further from break-up threshold ($We=20$, last row), the break-up happens faster but also leads to a very thin filament leading to many tiny drops. It is clear from these examples that the surrounding flow leading to break-up is different, with high stretching being required to create the elongated structure necessary for break-up with a possible importance of the internal flow of the drop.

3.2. Critical Weber number

We analyse the break-up threshold as a function of the various non-dimensional numbers. Figure 2(a) shows the break-up probability obtained from ensembles of simulations at various $We$ and $\mu _r$, for $Re_{\lambda }=77$. Each point corresponds to an ensemble of ten runs with various precursors, and the condition is considered in the break-up region if more than 50 % of the cases broke before 20$t_c$ after insertion. The critical $We$ ($We_c$) is obtained from the contours of the break-up probability and is shown as a solid line, showing that $We_c$ is a function of $\mu _r$ for large viscosity ratio, at a given $Re$. Grid resolution is verified by running the same ensemble with higher refinement level. The simulation length of 20$t_c$ was chosen as a trade-off between a long enough sequence to properly define the break-up occurrence and a reasonable computational cost to perform large ensembles for a wide range of conditions. This time corresponds to 200$\tau _\eta$, the Kolmogorov time scale, so that the time we systematically explore is much longer than the times tested in Håkansson & Brandt (2022), giving strong confidence in the resolution of the break-up boundaries. Note that the time for pinching is smaller than the eddy-turnover time, so when no break-up is observed at high viscosity ratio, it is not due to an incomplete pinching but simply because the conditions leading to break-up have not been met. As summarized in table 1, we performed close to 10 000 individual simulations at relatively high resolutions (interface resolution equivalent to $512^3$ to $1024^3$ grid points), systematically exploring high viscosity ratios, which is the focus on the study, and not the study of rare events, which would require a much larger number of eddy-turnover times.

Figure 2. Drop break-up phase diagram, as a function of $We$ and viscosity group: $Oh$ (b); $\mu _r$ (c); and inverse of a drop Reynolds number (d). (a) Example of data analysis leading to extraction of the critical-$We_c$ line for $Re_{\lambda }=77$ and ${\rm \Delta} x/d_0=64$ (L9). Each symbol is for an ensemble of simulations obtained using different precursors, and the percentage of cases broken up within $20 t_c$ is colour-coded. The critical-$We_c$ line is extracted as the 50 % contour. (b) Critical-$We_c$ contours for various resolutions (circles $d_0/{\rm \Delta} x=64$, L9; squares $d_0/{\rm \Delta} x=128$, L10; and crosses $d_0/{\rm \Delta} x=32$, L8) and increasing $Re_{\lambda }$. Grid convergence is observed between the various resolutions at a given $Re_{\lambda }$, as shown by the overlapping symbols. At low $Oh$ (typically $Oh < 0.1$), a constant $We_c$ is observed, close to 2.5. At high $Oh$ ($Oh > 1$), $We_c$ becomes a function of $Oh$ and $Re_{\lambda }$. (c) (3.1) shown as dashed lines introduces a characteristic viscosity ratio $\mu _r^0$ above which viscosity controls the break-up, with $\mu _r^0\approx 1.5 Re_{\lambda }$ (inset). (d) Break-up boundaries $We_c$ can be rescaled by the inverse of a drop Reynolds number $(({\rho _l}/{\rho _d})({\rho _d u_{rms} \lambda }/{\mu _d}))$.

Figure 2(b) shows $We_c$ contours as a function of $Oh$, with $Oh$ ranging from 0.01 to 10, and $We$ from 0.5 to 30. We observe that, at low $Oh$ (low viscosity), break-up occurs at $We$ close to $We_c\approx 2.5$, in agreement with values reported in the literature for low viscosity and bubble conditions (Rivière et al. 2021). At high $Oh$ (higher viscosity), we observe that $We_c$ is now a function of $\mu _r$. Additionally, higher $We$ needs to be reached to obtain break-up, and there is a weaker influence of the Reynolds number $Re_{\lambda }$. Here, $We_c$ is grid converged as shown in figure 2(b) by the different overlapping markers corresponding to increasing resolution.

Figure 2(c) shows $We_c$ contours as a function $\mu _r$, for different $Re_{\lambda }$. We observe that $We_c$ is a function of $\mu _r$ alone, with the critical-Weber contour described by a constant value at low $\mu _r$, $We_c=We_c^0$, with $We_c^0\approx 2.5$, which can be seen as an absolute value at low viscosity. Note that the particular value might be slightly sensitive to the choice of averaging, i.e. the size of the ensemble, the 50 % threshold and duration of the simulations; here, 20$t_c$ as for any statistical definition of critical break-up conditions in a turbulent flow. However, the critical value is within the bounds of the literature, and the important physical result is the existence of the critical value at low viscosity ratio; while when the viscosity ratio increases, $We_c$ becomes a steep function of $\mu _r$ for $\mu _r > 10$. Overall the data can be described by an empirical function

(3.1)\begin{equation} We_c = {We_c^0} \exp{\left({\mu_r}/{\mu_r^0}\right)}, \end{equation}

with the critical viscosity ratio $\mu _r^0$ obtained for each $Re$ by a least-squares fit. In the present simulations, $\mu _r^0$ is found to scale linearly with $Re_{\lambda }$: $\mu _r^0 = 1.5 Re_{\lambda }$; see inset. A similar discussion could be provided using $Oh$. The transition from the inertial to the viscous regime can be interpreted as a constant drop viscosity Reynolds number, since $\mu _r^0 \approx 1.5 Re_{\lambda }$ is equivalent to $1/(({\rho _l}/{\rho _d})({\rho _d u_{rms} \lambda }/{\mu _d}))\approx 1.5$, so that (3.1) becomes

(3.2)\begin{equation} We_c = {We_c^0} \exp{\left(\frac{1}{1.5}\frac{\rho_l}{\rho_d}\frac{\rho_d u_{rms} \lambda}{\mu_d}\right)}, \end{equation}

which rescales all data onto a single curve, as shown in figure 2(d).

Note that the functional form in (3.1) is the simplest one we could find, with the smallest number of parameters, but is not justified by a physical argument at this point and that other functional forms could be applied. We also comment that we did not vary independently the density ratio when demonstrating the rescaling (3.2), and that the range of Reynolds numbers tested remains limited (with $Re_\lambda$ from 40 to 80) due to the higher computational cost for higher Reynolds numbers. While the values of the Taylor-scale Reynolds number are modest, they are high enough to observe a reasonable inertial range, as shown in Perrard et al. (2021), Rivière et al. (2022) and Farsoiya et al. (2021, 2023) and correspond to integral-length-scale Reynolds number from 100 to 400 approximately (with a scaling between the integral-length-scale and Taylor-scale Reynolds numbers being $Re_{int}\propto Re_\lambda ^2$; see Pope 2000).

We can then define a viscosity group–dependent Hinze scale $d_h$ that separates stable drops ($d< d_h$) from those that will break ($d>d_h$) using the viscosity-dependent $We_c$ for the break-up $We_c(\mu _r)$: $d_h=({W e_c(\mu _r)}/{2})^{3 / 5}({\sigma }/{\rho })^{3 / 5} \varepsilon ^{-2 / 5}$.

3.3. Drop break-up frequency

Within the break-up regime shown in figure 2, we extract the (first) break-up time of the drop $T$ for each realization that breaks and compute its ensemble average $\langle T \rangle$. Figure 3 shows the non-dimensional break-up frequency $t_c/\langle T\rangle$ as a function of $We$ for various $\mu _r$ at $Re_{\lambda }=77$. The break-up frequency increases with $We$ while reaching saturation values that increase from high to low $\mu _r$. The overall trend for each $\mu _r$ is similar to what has been observed for bubbles (Martinez-Bazan et al. 2010; Rivière et al. 2021). The data can be described by a formula based on energy arguments (Martinez-Bazan et al. 2010),

(3.3)\begin{equation} \frac{t_c}{\langle T\rangle } = \alpha(\mu_r) \sqrt{\frac{We}{We_c(\mu_r)}-1}, \end{equation}

where $\alpha (\mu _r)$ is a parameter that is obtained by least-squares fitting of the data, while $We_c(\mu _r)$ is given by (3.1). All data can then be rescaled using the parameter $\alpha (\mu _r)$ as shown in figure 3(b) for a single $Re$. The same procedure is applied at each $Re$ (and resolution), and shown in figure 3(c), demonstrating grid convergence (symbols of increasing sizes correspond to increasing resolution and overlap). The values of $\alpha (\mu _r)$ (expressed in terms of the drop Reynolds number) are shown in figure 3(d), with a maximum value of 0.3 at low viscosity and decaying exponentially to reach a value of 0.1 at high viscosity.

Figure 3. (a) Ensemble-averaged drop first break-up frequency $t_c/\langle T \rangle$ made dimensionless using the eddy-turnover time at the scale of the initial drop, for increasing $\mu _r$, at $Re_\lambda =77$ (and L9). Equation (3.3) is fitted to the data for each $\mu _r$ (dashed lines), defining $\alpha (\mu _r)$. (b) Same data showing the rescaled break-up frequency $t_c/\alpha (\mu _r)\langle T \rangle$ as a function of the distance to $We_c(\mu _r)$. (c) Data for L8 (small symbols), L9 (medium) and L10 (large) show grid convergence and various $Re$. All data collapse onto a single curve. (d) The value of $\alpha (\mu _r)$ is a function of $\mu _r$ and $Re_\lambda$, and decreases with the drop Reynolds number.

Note that, as discussed by Håkansson (2020), the use of a break-up frequency within a population model must be made in a consistent way. Here, our definition of the break-up frequency has no ambiguity in the region of the break-up diagram where all cases break before 20$t_c$ (equivalent to 200$\tau _\eta$, significantly larger than any of the break-up times discussed by Håkansson & Brandt 2022), which is almost all data points except the ones on the boundary (see figure 2). For the parameters close to the boundary, the mean break-up time can become close to the window of observation, and less than 100 % of the ensemble leads to break-up.

3.4. Droplet production and size distribution

We now discuss droplet production and size distribution as a function of distance to critical conditions (quantified by $We/We_c(\mu _r)$) and viscosity ratio $\mu _r$.

Figure 4(a,b) shows the child-size distribution of the first break-up obtained from large ensembles ($\approx$100 runs per condition), for $\mu _r=0.01,10,150$ and $We/We_c\approx 1.5$ and 3.2 (corresponding to $We=4$, 8 or 20 depending on $\mu _r$), close to critical conditions (a, $We/We_c\approx 1.5$) and far from critical conditions (b, $We/We_c\approx 3.2$).

Figure 4. Number and size of drops production. (a,b) Child-size distribution of the first break-up event for the six large ensembles (order 100 realizations). Volumes are normalized by the volume of the initial parent droplet. The symbols are DNS data (L10), while the solid (dashed) lines are the probability density function from L10 (L9) data estimated using a Gaussian kernel density estimate, close to (a) and far from (b) critical conditions. (c) Ensemble-averaged number of drops $\langle N\rangle$ produced after the first break-up as a function of time, for the large ensembles, counting drops with more than 4 grid points per diameter at L9 (solid line L10, dashed line L9, showing reasonable grid convergence). (d) Number of drops $\langle N \rangle$ at $(t-{T})/t_c=2$ as a function of $We/We(\mu _r)$ for different viscosity ratios. Circles are from small ensembles (10 realizations) and squares are from L10 large ensembles. (e,f) Drop-size distribution for the large ensembles taken at $(t-T)/t_c=2$. Drops larger than $4 {\rm \Delta} x$ are grid converged, close to (e) and far from (f) critical conditions. The shape of the distribution presents large variations depending on distance to criticality and viscosity ratio.

For low to intermediate $\mu _r$ ($\mu _r=0.01,10$), at $We$, close to critical conditions (a) ($We/We_c(\mu _r)\approx 1.5$), the child-size distribution presents an M shape, compatible with previous descriptions of drop or bubble break-up in turbulence close to critical conditions (e.g. Rivière et al. 2021), corresponding to the production of two children of volumes 0.7$V_0$ and 0.3$V_0$ as the most probable event. At high viscosity, here $\mu _r=150$, break-up at $We=4$ does not happen since $We_c(\mu _r=150)\approx 6.5$. When considering $\mu _r=150$ close to critical conditions ($We/We_c(\mu _r)\approx 1.5$), the child-size distribution clearly changes compared with low viscosity break-up close to critical conditions, with the production of two identical volumes being the most probable outcome, illustrating the change in dynamical regime for the first break-up of viscous drops.

When increasing $We$ moving away from critical conditions (b), $We/We_c(\mu _r)\approx 3.5$, at low viscosity ratio, we observe a transition to a U shape, again compatible with previous descriptions (Rivière et al. 2021; Ruth et al. 2022), with a very small and a very large drop being formed, with similar results for $\mu _r=0.01,10$. However, at a high viscosity ratio, the distribution is close to a bell shape (similar to the shape of the distribution at lower $We$ for the same viscosity ratio) corresponding to two lobes of equal size, just prior to the fragmentation into tiny droplets of the snake like shape seen in figure 1.

Models have been developed to rationalize the observed shape of the child-size distribution of a binary break-up. As summarized by Qi et al. (2020) and Rivière et al. (2021), models for the first break-up child-size distribution lie in three categories: bell shaped (Martínez-Bazán, Montanes & Lasheras 1999), U shaped (Tsouris & Tavlarides 1994; Luo & Svendsen 1996) or M shaped. Bell-shaped models correspond to two equal drops as the most probable outcome, while very unequal break-up that creates a small and a large bubble has a high probability in U-shaped models. Finally, the M-shaped models postulate that there is an unequal break-up that is the most likely. The three models correspond to different assumed physics of the break-up. A mechanism proposed to explain the M-shape distribution is that eddies smaller or equal to the drop size collide with the initial drop, which then excite and break the drop in one eddy-turnover time (see Qi et al. 2020), while the idea of various eddies leading to oscillations at the drop or bubble resonance frequency has also been proposed (Risso & Fabre 1998). Using the M-shape in a population balance model would require physical configurations where the break-ups are independent from one another and always at moderate conditions (close to critical conditions). A mechanism proposed to explain the U-shape distribution can be based on minimizing the surface energy increment, where the difference between the total surface energy before and after break-up controls the probability and a smaller difference means a larger probability for such break-up to happen. The U-shape distributions can be seen as resulting from the tearing of a small part of the initial drop (or bubble) and can correspond to the beginning of a more complicated break-up sequence. The U-shape is therefore often used in population balance models, as it allows us to recreate a rapid sequence as long as the model for the successive break-up times is adequate (see Ruth et al. 2022). A discussion on the bell shape is provided by Martínez-Bazán et al. (1999), and relates to the assumption that the probability of generating a child drop of a specific size is proportional to the product of the excess stress of the two children, where the excess stress is the difference between the turbulent dynamic pressure around the children and the capillary pressure of the parent.

To summarize, the first break-up child-size distribution of drops in turbulence is controlled by two parameters: the distance to critical conditions, with a transition from M-shape to U-shape at low viscosity ratio, and the viscosity ratio itself, with a bell shape both close and far from critical conditions.

However, we caution that the analysis of the first break-up child size only gives part of the story hinted at by the very different geometrical properties of the break-ups (see figure 1) and the significant history effect being present when considering longer time evolution and sequences of break-ups. Indeed, at high viscosity ratios, the first break-up is described by a bell shape, which misses the production of tiny drops related to the snake like geometry.

Figure 4(c) shows the number of drops being produced as a function of time after the first break-up for the same ensembles ($t-T$) (each curve starts with 2 drops at $t-{T}=0$). Close to critical conditions ($We/We_c\approx 1.5$), we observe a mild increase up to 3 to 4 drops, as very few break-up events occur, compatible with the classic Kolmogorov–Hinze cascade of drop break-up, producing children of similar sizes and the cascade ends at the Hinze scale.

Now considering a high viscosity ratio, $\mu _r=150$, still close to critical conditions ($We/We_c(\mu _r)\approx 1.5$, which is now $We=8$ instead of 4) we observe significantly more drops formed over the time period than for the low viscosity cases, with up to 10 drops produced, corresponding to the final break-up of elongated filaments, even close to critical conditions.

As we move away from critical conditions, $We/We_c(\mu _r)=3.2$, for low to intermediate $\mu _r$, we observe a strong increase of the number of drops (similar to that reported for bubbles), with 15 to 20 drops at $(t-T)/t_c=2$ (with the number of drops increasing with increasing viscosity ratios). At large viscosity ratio and similar $We/We_c(\mu _r)=3.2$, we now observe more than 50 droplets related to the extreme stretching of the ligament before break-up.

The number of drops produced is therefore a function of both the distance to critical conditions and the viscosity ratio, with both parameters leading to an independent increase of the number of drops produced. This is confirmed by figure 4(d), which shows $\langle N \rangle$ (taken at $(t-{T})/t_c=2$) as a function of $We/We_c(\mu _r)$. For $\mu _r$ up to 10, we observe a linear increase in the number of drops with $We/We_c(\mu _r)$, independent of the viscosity ratio. The break-up is inertial with no influence of viscosity, similar to the dynamics of bubble break-up described by Rivière et al. (2021, 2022). When the viscosity ratio is further increased, the number of drops increases much more sharply.

The drop-size distributions for the large ensembles are shown in figure 4(e,f), taken at $(t-{T})/t_c=2$. After two eddy-turnover times of the initial drop, the effect of the complete break-up sequence can be discussed on the drop-size distribution, and the role of distance to criticality and viscosity ratio is further confirmed. The time evolution of the drop-size distribution is shown in Appendix B and confirm that analysis after 2$t_c$ is robust and that we have large enough statistical sample to analyse the shape of the drop-size distribution.

At low $\mu _r$ and $We/We_c(\mu _r)\approx 3$, we observe a distribution close to the $d^{-3/2}$ scaling, reminiscent of capillary processes driving fragmentation of large volumes (Rivière et al. 2022; Ruth et al. 2022), with a shape that is not a direct result of the first break-up, as rapid break-ups follow the initial one, leading to a broad sub-Hinze population.

Close to critical conditions, at low viscosity ratios, we observe a size distribution coherent with the child-size distribution of the first break-up with mostly drops close to the Hinze scale. The behaviour of the highest $\mu _r$ is quite different, with a large number of drops being produced, but their sizes are all concentrated at the smallest possible scale being resolved, i.e. the scale of the filament at break-up, which is a consequence of the very strong thinning visible in figure 1. Such a shape could not be inferred from the first break-up, which displays two lobes of similar sizes, as the thin filament has not yet broken into multiple small drops.

4. Conclusions

We performed large ensembles of DNS of drop break-up in isotropic turbulence, spanning orders of magnitude in $\mu _r$ and $We$, while also exploring the sensitivity to $Re$. We report a break-up occurrence diagram, which confirms a constant $We_c$ for low $\mu _r$, while $We_c$ becomes a function of a drop Reynolds number combining $\mu _r$ and $Re$ when $\mu _r$ increases. The first break-up time is found to be a function of the distance to $We_c$, with an increased break-up time as drop viscosity increases. Finally, the drop production and size distribution are analysed. We demonstrate the limits of only looking at the first break-up as the break-up sequences over an extended period of time present very different statistical behaviour, with different dynamical regimes at low and high $\mu _r$, and depending on the distance to critical conditions. This confirms that population balance models for drops should include a multi-time-scale drop-size distribution to account for break-up sequences, especially at high $We$ and high drop Reynolds number. The present work paves the way for a detailed analysis of the break-up mechanisms, including the role of the flow surrounding the drop prior to break-up and the internal drop dynamics for the different regimes identified.