Stochastic simulation of droplet breakup in turbulence

Abstract This study investigates single droplet breakup from a theoretical perspective and addresses whether breakup in turbulent flows can be studied using highly-resolved simulations. Transient and three-dimensional turbulent flow simulations are performed to investigate if the apparent stochastic outcome from the droplet breakup can be predicted. For a given turbulent dissipation rate the breakup events were simulated for various detailed turbulence realizations. For this purpose, a well-characterized system widely used for kernel development is utilized to validate the simulations with respect to the key characteristics of stochastic breakup, including droplet deformation time, the number of fragments, and the specific breakup rate. The statistical validations show very good agreement with all the quantitative properties relevant to the breakup dynamics. Necklace breakup is also observed in line with patterns found in experiments. Evidence is found that the rate of energy transfer is positively correlated with higher order fragmentation. This can allow development of more accurate breakup kernels compared to the ones that only relies on the maximum amount of energy transfer. It is concluded that the simulation method provides new data on the stochastic characteristics of breakup. The method also provides a means to extract more details than experimentally possible since the analysis allows better spatial and temporal resolutions, and 3D analysis of energy transfer which provides better accuracy compared to experimental 2D data.

• Analysis of stochastic breakup and validation with experimental data.
• Qualitative and quantitative agreements between simulations and experiments.
• Higher order fragmentation linked to intensity of vortex interaction. This study investigates single droplet breakup from a theoretical perspective and addresses whether breakup in turbulent flows can be studied using highly-resolved simulations. Transient and three-dimensional turbulent flow simulations are performed to investigate if the apparent stochastic outcome from the droplet breakup can be predicted. For a given turbulent dissipation rate the breakup events were simulated for various detailed turbulence realizations. For this purpose, a well-characterized system widely used for kernel development is utilized to validate the simulations with respect to the key characteristics of stochastic breakup, including droplet deformation time, the number of fragments, and the specific breakup rate. The statistical validations show very good agreement with all the quantitative properties relevant to the breakup dynamics. Necklace breakup is also observed in line with patterns found in experiments. Evidence is found that the rate of energy transfer is positively correlated with higher order fragmentation. This can allow development of more accurate breakup kernels compared to the ones that only relies on the maximum amount of energy transfer. It is concluded that the simulation method provides new data on the stochastic characteristics of breakup. The method also provides a means to extract more details than experimentally possible since the analysis allows better spatial and temporal resolutions, and 3D analysis of energy transfer which provides better accuracy compared to experimental 2D data.

Introduction
Physical understanding of droplet breakup benefits the design and optimization of multiphase processes. It is important to predict how the size distribution of the dispersed phase evolves, how fast the breakup process occurs, and how many daughter droplets are produced in each breakup (i.e. binary, tertiary, or higher order fragmentation). improved control of the interfacial area for mass, heat and momentum transfers [1,2]. The mathematical framework widely used to study multiphase systems involving the breakup of the dispersed phase is computational fluid dynamics (CFD) coupled with the population balance model (PBM) [3][4][5][6]. Attention has been given to the development of mathematical models for the breakup process in the past few decades [7][8][9] since the breakup kernel is a determining factor for the accuracy of the CFD-PBM approach [10]. A common strategy is to formulate the breakup as the encounter between a fluid particle and a turbulent structure with regard to the frequency and probability of their interactions. This methodology, breakup based on the encounter between fluid particle and vortices, has been widely used to develop various breakup rate models [11][12][13][14][15][16][17][18]. Furthermore, breakup kernels are typically evaluated by validation against experimental data for daughter size distributions or breakup rate measurements. There is a twofold hurdle to overcome for comprehensive validation of breakup rate models. The number of data sets for validation that includes breakup rates, the number of fragments, and daughter size distributions is scarce. Moreover, characterizations of hydrodynamic properties, such as turbulent dissipation rate and turbulent kinetic energy, often introduce uncertainty into model validations. For instance, Castellano et al. [19] introduced an alternative zero-dimensional method to account for the flow inhomogeneities without performing 3D CFD simulations. Few breakup measurements in the literature attempt to determine the breakup properties inside different devices, such as stirred tanks, pipe flows, and nozzles [20][21][22][23][24][25][26][27][28][29]. Eastwood et al. [21] for instance, observed that fluid particles can deform to a length scale comparable to the integral length scale of the turbulent field, whereas Andersson and Andersson [30] reported the deformation ratio (i.e., the ratio between the long axis and the initial droplet diameter) of up to 20.
Other good examples are works of Maaβ et al. [22,23]. They proposed an extended beta function for the daughter size distribution that depends on Weber number and droplet diameter for liquid-liquid dispersion in stirred tank. Later, they studied breakup of toluene and petroleum droplets by measuring the breakage time and breakage probability with a constant turbulent dissipation rate. Solsvik et al. [31] studied the breakup of single droplet in a stirred tank with an average value for the turbulent dissipation rate. They differentiated between the initial and the entire cascade of breakup that could result in 47 daughter droplets for a 1 mm toluene droplet. For the same droplet diameter, however, the initial breakup generated on average 2.39 daughters. It is inferred from the examples above that the stochastic nature of breakup significantly varies the research outcomes. In addition to that the hydrodynamics of the primary phase (e.g. the turbulent dissipation rate) is difficult to measure locally. For instance, Vankova et al. [32] applied CFD simulations to obtain the turbulent dissipation rate values when studying the effects of different parameters (e.g., dissipation rate, interfacial tension, and viscosity) on the breakup rate and daughter size distribution of emulsions. Furthermore, the number of stable droplets or bubbles (i.e., where the cohesive stresses dominate the disruptive stresses) has not been measured, and the breakup rate has been indirectly calculated as a solution for the population balance equation. Thus, the uncertainty associated with the approximations adopted in previous works has led to inaccurate models for breakup rate predictions. We have chosen two sets of experiments that report both flow hydrodynamics and breakup dynamics, including the number of droplets that are stable and do not break. The reason is to avoid introducing undesirable errors due to measurements and establish a concrete validation study. More details about these measurements are presented in the original articles [15,33].
An attractive alternative to validation with experiments is validation using highly-resolved simulations, since these simulations provide detailed information on the physics governing the characteristics of the fluid particle breakup, limited only by the computational resources available. There are essential difficulties in conducting experiments for validating a breakup kernel. For instance, one can only investigate the breakup phenomenon by monitoring the behavior of the fluid particle. The associated vortex (or vortices), on the other hand, that transfers energy, deforms the fluid particle, and causes the fragmentation is not captured. There is an explicit coupling between data resolution and the interrogation window size when high video imaging is utilized. In other words, the width of the interrogation window is determined by the deformation time and the breakup time of the fluid particle. The window should be large enough that the entire deformation and the droplet breakup event take place within the observational zone. Otherwise, breakup occurrences might be unnoticed outside this region, which leads to erroneous breakup rate measurements. This effect becomes even more significant when the multiphase system functions under different local Re numbers conditions. In other words, a careful selection of interrogation window size should be conducted with respect to the local Re number covering the entire breakup phenomenon (i.e., deformation and breakage). Secondly, the light sources must be intense enough to allow short exposure time to capture sharp images of the drop dynamics. Thirdly, image resolution should be enough to identify the smallest fragments formed during fragmentations. Thus, the interrogation window should be selected according to the different Re numbers and this requirement is troublesome during experimentation. By contrast, CFD simulations decouple the resolution of the data obtained from the domain size, thereby circumventing the limitation that is present in experimental data. With increasing computational power, the additional cost due to the variation of domain size is feasible. Table 1 summarizes major differences between the two available possibilities for the validation of breakup kernels. The data in the table indicate that once the highly-resolved CFD simulation method has been validated, it can replace the experiments, and builds confidence in the validation of breakup rate models. Highly-resolved CFD simulations, built on Direct Numerical Simulation (DNS) or Large Eddy Simulation (LES) in conjunction with a method to capture the deforming interface, have been performed to investigate the dynamics of fluid particles in multiphase systems. The interface topology is quite complex and typically two approaches are employed to obtain the interface deformations. The first method tracks an indicator function to capture different phases on Eulerian grid (e.g., Volume of Fluid, VOF, and level-set methods), while the second approach connects marker points to track the evolution of the interface. One of the earliest works is the study of buoyant bubbles with deformed interfaces in a periodic domain using DNS [34,35]. Low and moderate Reynolds number values have been investigated to explain the breakup behaviors of two-and three-dimensional representations of bubbles. To avoid numerical diffusion of the density and viscosity at the interface, they used a hybrid method between front capturing and front tracking that reconstructed the grid at each step. Unverdi and Tryggvason [36] explicitly tracked the interface for multi-fluid flows with sharp interfaces using a front-tracking method and included the influence of surface tension. The drawback of this method is that when two interfaces are close enough, they can cancel each other. Li et al. [37] have numerically confirmed the application of VOF in addressing the deformation of high viscosity droplets in a shear flow, while the DNS-VOF method has been validated with experiments to study the collision of liquid droplets [38]. The LES-VOF framework has also been utilized to understand the interactions between turbulence and deforming interfaces in a downward injection of air into an infinite water pool [39]. The VOF method conserves the liquid volume but extracting geometrical properties requires numerical efforts. However, the level-set method by Sussman et al. [40] facilitated convenient computation of interface normal and curvature. Recent advances and applications of level-set methods have been reviewed by Gibou et al. [41]. To combine the conservative benefits of VOF and the easier estimation of curvature in level-set method CLSVOF technique was proposed by Sussman and Puckett [42]. An example of this approach is the work of Xiao et al. [43] on the droplet formation and breakup for atomization modeling. Baraldi et al., presented a VOF/projection method for droplet-laden turbulent flows. They applied Eulerian-implicit-Eulerian-Algebraic-Lagrangian-explicit algorithm as the VOF advection scheme which only required three advections and reconstruction per time step. An Initially spherical object was simulated by this method and then it was coupled with a projection method using the continuous surface model to compute the surface tension force. They argued that the benefit of their method was fewer steps for advection and reconstruction. Komrakova et al. [44] performed 3D simulations of a periodic cube using a free energy lattice Boltzmann equation to study drop breakup and coalescence. They used an interface-tracking method and stated that the drawback of this method was the dissolution of small droplets. Recently, Shao et al. [45] have applied DNS with a modified version of the level-set method to investigate the effect of the We number on the breakup process. Their simulations showed a sudden bursting of the droplet without deformation that could be related to the initialization of the simulations. Nevertheless, for typical chemical engineering applications and kernel developments, the droplet may experience stretching and deformation through interaction with one or a pair of vortices. In the case two vortices simultaneously act to deform and break a drop it provides larger probability for breakup. Considering that the probability of simultaneous interaction with two vortices of sizes, both in the range influencing breakup, is lower than interaction with individual vortices, this can possibly explain experiments where size distribution changes over long time (hours). On the other hand, we found no evidence that energy is accumulated in between which is logical since there is a simultaneous dissipation of energy. The mentioned works show examples of DNS or LES methodologies for studying the dynamics of fluid particle deformations and breakups. For more examples, interested readers can refer to [46][47][48][49][50]. Another industrial context in which the interaction between the continuous and the dispersed phases is crucial for process performance is atomizing two-phase flows (e.g. fuel injection and aerosol spray). DNS is a more recent option in this field to interpret the unsteady mechanisms of two-phase phenomena that occur during primary breakup. For instance, efforts have been made to study the various features of interface dynamics, such as the formation of ligaments throughout the primary breakup in gas-liquid jets. The efforts have focused on the direct solution of Navier-Stokes equations incorporated into an interface capturing/tracking method [43,[51][52][53][54][55].
As seen, highly-resolved simulations when combined with an efficient method for capturing the interface topology have become a feasible option for studying multiphase systems and the apparent stochastic processes involved. Thus, in this work, we continue the previous work of Andersson and Helmi [56] by performing highly-resolved large eddy simulations combined with VOF. The objective is to construct a validated numerical methodology to study the breakup phenomenon. The detailed data obtained using highly-resolved simulations can be employed for developing closure terms for PBM. Therefore, this work evaluates the alternative approach of highly-resolved simulations for the development of breakup kernels. This article also elaborates on the dynamics of breakup using the information extracted from the simulations to understand the apparent stochastic process. This method ultimately provides a better understanding of the different breakup dynamics, including deformation time, rate of energy transfer, specific breakup rate, and the characteristics of daughter particles. The novelty lies in the statistical validation of a stochastic phenomenon such as single droplet breakup in turbulence. The statistical validation of the single droplet breakup has not been demonstrated before. The flow field close to the homogeneous turbulence in the mixing element of a flow reactor was simulated with space resolution much finer than the dispersed phase length scale, and with a time resolution finer than the breakup timescales. Quantitative and qualitative comparisons were made with a well characterized water-dodecane system [15] to validate the specific breakup rate (breakup frequency), the number of fragments, and the deformation time. The main contribution of this work is to build, analyze, and validate a high-fidelity tool for studying the single droplet breakup phenomenon.

Methodology
LES models are based on the idea that intermediate to large turbulent structures are explicitly resolved, while smaller eddies, that are more isotropic in nature, are modeled using sub-grid-scale models. This necessitates a spatial filtering operation to decompose the velocity field U x t ( , ) into a filtered one, U x t ( , ), and an unresolved or instantaneous component, u x t ( , ) ' [57]. Applying the filtering operation to the governing equations of two-phase incompressible flow results in: (1) In Eq. (2), g denotes the gravitational acceleration, and µ is the continuous phase viscosity. To close the momentum equation, the two unknown terms including the sub-grid-scale turbulent viscosity, µ t , and interfacial tension force, F i , must be formulated. The former is modeled in this work using the dynamic Smagorinsky-Lilly model [58]. Liovic and Lakehal [39,59] have formulated the filtered momentum equation with four additional sub-grid terms to resolve the interface properties. An indicator function for volume fraction is solved to account for the effect of surface tension. The incompressibility assumption for the flow assures that the density of the fluid particle is invariable throughout its course. Thus, a standard partial differential equation (with accumulation and convection) can be constructed for the indicator function. The indicator function ( ) approximates the position of the interface between the droplet and the continuous phase. The function is, in fact, a discontinuous step that is zero for the droplet, one for the continuous phase, and < < 0 1 for the transition region. The function is further coupled to each fluid, and the function propagates with them. The transport equation for reads as: The local density ( ) and the viscosity (µ) are averaged, accordingly: The subscript c in Eqs. (4) and (5) is for the continuous phase and d is for the droplet. The surface tension force can be applied using the Continuum Surface Force (CSF) approach, which adds a momentum source term to the balance equation [60]. Alternative methods are also available to determine the effect of surface tension. A recent review paper by Popinet [61] has summarized the different alternatives and discussed their accuracy and robustness.
The CSF method defines a normal vector to the interface, which is zero everywhere except in the transitional area between the fluids. The gradient of defines the normal vector to the interface: And the curvature of interface ( ) can be expressed as: Surface tension force has a normal component for curved interfaces such as a deformed droplet. This component for fluids in equilibrium is mechanically balanced by the pressure jump at the interface. Therefore, the interfacial tension force can be written as: The governing equations are now closed, and the equations can be solved using the finite volume method. The governing equations were solved using Ansys Fluent 18 on a high-performance computing cluster. In this work the temporal discretization relies on the second order accurate implicit transient formulation, while the spatial discretization of the momentum transport equations is based on central differencing. In this study the CICSAM discretization is used for the spatial discretization of the indicator function using a low value, × 1 10 6 , of the volume fraction cut-off and high resolution of the computational mesh, , to ensure volume conservation of the dispersed phase. The high grid resolution requires using short time step, = t µs 1 which means the simulations are computationally expensive. Considering also the fact that the system is simulated at conditions where drops do not always break during observation, some simulations provide no information on the breakup statistics. On the other hand, this data is very important as it provides the opportunity to judge the balance between disruptive and cohesive stresses, and ensure the breakup is not due to numerical effects. Additionally, the cumulative data on droplet breakup and no breakup events allows quantification of the breakup frequency. This is a key property of the breakup kernels and this has never been studied using highly resolved simulations.
The computational domain includes an individual mixing plate of a multipurpose reactor that produces a turbulent field close to homogeneous. The assembly of this flow reactor permits the continuous production and dissipation of turbulence due to repeated small mixing components [62]. PIV measurements for single phase was performed to obtain the turbulent properties of the system. This reactor delivers an ideal environment to study the dynamics of single droplet breakup both numerically and experimentally. The reactor has transparent walls allowing excellent opportunity for using intense illumination along with observing the droplet without optical distortion. An 8-bit CCD camera with a high frequency of data sampling rates (between 1 and 4 kHz depending on the system) was used to capture the breakup for dodecane droplets with density and viscosity of 750 kg m 3 and × 1.5 10 Pa s 3 , respectively Measurements were made using deionized water for the continuous phase and dodecane for the dispersed phase. In all visualization experiments, the temperature was kept constant at 20°C. The fluid properties relevant to the study are given in Table 2.
Additional details about the experimental setup are provided in Table 3.
The experimental data including the flow hydrodynamics [63] and the characteristics of the droplet breakup have also been reported earlier [15,56]. The influence of coalescence is negligible due to the low dispersed phase hold-up during the experiments (i.e. single droplet), and the continuous phase turbulence is not modified by the presence of droplets.
Short exposure time, µs 10 , was used to obtain sharp images of the breakup event. A linear velocity of 1.1 m s means the translational distortion is limited to µ 11 m, which compared to the mother droplet size of 1mm is only 1.1%. However, the main uncertainty comes from the image resolution which was 309, 000 pixels cm 2 . For 1mm mother droplets the error with respect to the resolution corresponds to an uncertainty of approximately 4%. Drop fragments smaller than µ 50 m cannot be identified. However, breakup measurements revealed only 3% fragments smaller than µ 100 m thereby the influence on mean number of fragments is very low. The characteristics of the hydrodynamics are directly controlled by the accuracy of the flowrate and fluid temperature control. Calibration of the flow rate ensured an error less than 1%. Temperature was kept constant at°20 C and thermocouple measurements are within ±°0. 1 C, which means that the uncertainty of fluid density and kinematic viscosity have uncertainties of 0.002% and 0.2%, respectively.
The single droplet was patched randomly within the domain. To avoid the dependency of results on the initialization method, the droplet was analyzed after the simulation was run a period equal to the ). The simulation and data analysis methodology are summarized by the following 8 steps.
(1) Simulation of the time filtered turbulent flow field using a turbulence model (e.g. k-ω) that allows near wall resolution of = + y 1.
(2) Superimpose synthetic turbulence on the mean flow. Note that on the short time scale < t l , the flow field is deterministic. This is because the energy is attributed to turbulent structures on different scales. This allows mesh independence to be studied. On the other hand, for longer time, turbulence is stochastic which means that a drop introduced at a different time (step 7 above) will interact with other turbulent structures. Sometimes these interactions are sufficiently strong to break the drop, and sometimes the disruptive stress does not exceed the cohesive stress and breakup does not occur.
The surface tension coefficient between water and droplet is = 0.053 N m 1 . Moreover, the Taylor-scale Reynolds number (Re ) for this system is 77, and the Taylor length scale ( T ) equals µ 320 m. The interface cells, where there are high gradients of , were dynamically refined, and the mesh for the bulk region of the flow was coarsened. All fine scale changes in the interface down to a length of 2% of the mother diameter are resolved. Hence the flow also inside the drop is resolved in detail which means even contributions from internal viscosity on drop stabilization is accounted for. However, for this system this has no influence, since the Oh number is low, but the resolution is still needed to resolve the smallest fragments correctly. The dynamic mesh refinement around the interface provides a feasible model that balances computational time, memory requirements, and solution accuracy. The quality of LES data was further verified using the ratio of the resolved turbulent kinetic energy to the total turbulent kinetic energy. The ratio was greater than 85%, and it has been suggested that the simulations are well-resolved if the ratio is greater than 80% [1,64]. Furthermore, vortices are identified as coherent fluid regions that has a positive Q-criterion (the second invariant of u), and the pressure at the core of the turbulent structure is lower than ambient pressure [65]. According to the expression of the Q-criterion (i.e. Eq. (9)), this quantity shows a balance between shear strain rate (S) and vorticity magnitude ( ).

Results and discussion
The main aim of this section is to demonstrate that highly-resolved simulations can be utilized to study the single droplet breakup phenomenon in turbulent flows. The breakup dynamics are divided into three stages namely deformation, initial breakup, and the entire breakup. The statistics of each sub-process were validated with experiments to quantitatively determine whether the simulation data can be used to address the governing physics of the process. This, in turn, led to developing a universal breakup kernel for population balance modeling.

Deformation time
The breakup event consists of an initial deformation phase during which the droplets deform, and the interfacial energy increases when the disruptive stresses overcome the cohesive stresses, followed by a second phase during which droplets separate from the deformed mother drop. Solsvik divided the overall breakup event into initial breakup time and total breakup time [31]. In this study we quantify the deformation time, def , which is defined as the time required for a spherical droplet to deform from spherical shape until it reaches its most deformed state before any fragments are separated from the mother drop i.e.
= t t def d d0 . The dynamics of this stage are specified by the deformation time, def , the time required for the droplet to highly deform before fragmentation occurs. Ashar et al. [33] have introduced Eq. (10) to compute the deformation time of droplets with diameters within the inertial subrange of the turbulent energy spectrum. = l u def (10) In Eq. (10), l is the maximum deformed extent of the droplet. This maximum was experimentally characterized as a function of the aspect ratio, , and the mother droplet diameter, = l d ( 1) 0 . The aspect ratio of 5 has been recommended by Andersson and Andersson [30]. This recommendation is based on sequences of deformed droplets that are captured through high-speed imaging. The denominator in Eq. (10), u , represents the mean turbulent velocity of the eddy with the length scale . For the inertial subrange of turbulence, the Kolmogorov assumption for the mean turbulent fluctuating velocity of turbulent structures, u 2 ( ) 1 3 , relates the velocity to the length scale and the turbulent dissipation rate [57,66,67]. It should be noted the entire energy spectrum of the turbulence can be applied to obtain the turbulent eddy velocity term. The details of this mathematical procedure can be found elsewhere [68,69]. Introducing the concept of deformation time is of particular importance when studying the dynamics of breakup phenomena. This can provide a plausible explanation for the behavior of the breakup kernels when the deformation time is implemented. Most of the commonly applied breakup rate models show a monotonic increase with droplet diameter. However, this increase can be bounded by the application of def , considering that droplets need to deform before they can break. One must pay attention to the range of effective eddy sizes that interacts with the droplet, since the computation of deformation time is conditioned by the minimum and maximum sizes of eddies interacting with the droplet. The minimum vortex size can be associated with the Kolmogorov length scale or, as Andersson and Andersson [15] have suggested, with d 10 0 . Different limits have been employed for the maximum size of effective vortices, including d 3 0 , d 5 0 , and d 10 0 [15,16,[70][71][72]. To show an example of how def determines breakup behaviors, the experimental data reported by Ashar et al. [33] was used. A custom-made rotor-stator mixer was used in their work that enabled access to the mixing zone of the tank. A highspeed camera with 10-bit dynamic resolution, and narrow beam halogen lights guaranteed enough contrast between the droplet and the liquid phase. A plexiglass stator with six-blade rotor of diameter 188 mm were used. The clearance between the rotor and the stator was 0.5 mm with rotor speed ranged between 300 and 500 rpm. Direct breakup rate measurements were conducted in deionized water for single rapeseed oil droplets. The fluid properties are given in Table 4.
To ensure the balance between the energy transfer and heat removal the fluid temperature was kept constant at°20 C. The specific breakup rate reported for the mother droplet with the diameter = d µ 500 m 0 was = 214.3 s s 1 . It is concluded that for the systems studied the fragmentation of a droplet begins once the droplet reaches its maximum deformed state, and thereby deformation time should always be less than or equal to the inverse of the breakup rate (i.e. 1 def s ). For rapeseed oil droplet with = d µ 500 m 0 , the deformation time should be less than 4.7 ms. Fig. 1-a shows the normalized cumulative breakup rate of rapeseed oil droplet with = d µ 500 m 0 as a function of normalized vortex size, which was used to investigate the effective range of vortex sizes for the droplet breakup and deformation. The breakup rates were calculated based on our previous breakup rate model for the entire turbulent energy spectrum [68]. Fig. 1-a shows that the effective size of vortices contributing to the breakup phenomenon falls within the range of d d 0.5 2.6 0 0 for this test case. The associated deformation times for the upper and lower limits of the vortex length scales were 3.6 ms, and 2.2 ms, respectively. This confirms that, for the obtained range of turbulent structures, the droplet had sufficient time to complete the process of surface deformation (i.e. [2.2 3.6] ms 4.7 ms). Following the same procedure, one can apply the deformation time concept to determine the minimum and maximum boundaries for the breakup rate models. An example is shown in Fig. 1-b where the effective vortex sizes were implemented into the Karimi and Andersson [68] model, creating the upper and lower limits for the breakup rate predictions. The significance of this theory can be further emphasized by comparing the numerical predictions of different breakup rate models for a rapeseed oil droplet with = d µ 500 m 0 . The predictions for the breakup rate models of Coulaloglou and Tavlarides [12], Luo and Svendsen [11], Andersson and Andersson [15], and Karimi and Andersson [68] are 1659 s 1 , 1200 s 1 , 312 s 1 , and 188 s 1 , respectively. The corresponding values of 1 s , according to the model predictions, are 0.6 ms, 0.8 ms, 3.2 ms, 5.3 ms. It is evident that the timescale for the first two models was too short to physically accommodate the required time for droplet deformation, whereas the timescales for the latter two models reflect the possibility of a deformation occurrence within the available time window of the breakup. Therefore, one can conclude that obtaining detailed information about the deformation time results in a better understanding of the breakup phenomena which leads to the development of more predictive models for the breakup rate. Hence, it is very important that the high-resolved CFD simulations provide reliable data on the deformation timescale.
Direct comparisons between different observations and simulations are difficult due to the obvious stochastic nature of the phenomenon. This is illustrated in Fig. 2 where two instances of breakup are followed under the same hydrodynamic conditions and with the same fluid properties. The idea for this figure is to illustrate by two examples how the stochastic nature of turbulent fields can affect the dynamics of the breakup for the two cases. The deformation time of the first simulation was 6 ms with a final number of three fragments (Fig. 2-a), whereas the deformation time for the second simulation ( Fig. 2-b) resulted in seven fragments with = 5 ms def . To further explore the characteristics of the two simulations in Fig. 2, one can utilize the energy barrier plots introduced by Helmi and Andersson [56].
The energy barrier plot shows the normalized interfacial area of the deformed droplet as a function of timescale during the breakup process. The uncertainty of the surface integration is less than 2% due to the high resolution in all regions of deformed drops. The increase in the interfacial area is an indication of energy transfer from the turbulent structure to the fluid particle. The most deformed state of the droplet before the breakup is expressed by a maximum in the energy barrier plot. Fig. 3 shows the energy barrier plots for the simulations with three and seven daughter fragments (shown in Fig. 2). There are noticeable differences between the two plots in terms of maximum energy and the maximum derivative of the two curves. The increase in the interfacial energy is significant, and it was used to define breakup criterion in the breakup rate models. Analysis of 2D images, obtained by high-speed imaging showed 34% increase in the energy [15]. The 2D analyses of experimental data relied on simplifying the deformed objects as cylindrical objects with spherical ends. However, numerical simulations provide a more precise estimation of the increase in energy for the deformed 3D objects. The increase can easily be determined by the surface integral of locations in which the value of the phase function equals 0.5. Analysis of the 3D simulations gives a value of approximately 50% with a standard deviation equal to 10%. The maximum required energy obtained by the turbulent structure for the simulation with seven fragments is higher than the second simulation with smaller number of fragments. The maximum slope, or the derivative of the curves, can be interpreted as the rate of energy transfer. It can be seen from the figure that the fast rate of energy transfer from the turbulent vortex to the fluid particle leads to greater number of fragments. Thus, one should incorporate this into kernel development for multi-fragmentations. Currently, there is no rigorous mathematical paradigm to include the energy transfer rate into the breakup rate models. However, the formulation of breakup rate model introduced in our recent publication [68] allows to account for the increase in the surface energy which could be further improved to comprise the dynamics of energy transfer rate.
The droplet breakup occurred with different timescales and generated different numbers of daughter fragments. For this reason, the fragmentation statistics of the CFD simulations that include mean and standard deviations of the deformation times are also given in Table 5.
The uncertainty in determining data is due to the maximum frame rate used, which was Hz 1000 . Therefore, an error of 1 ms was expected, and this could also possibly have an effect on the standard deviations. However, both the camera frame rate and the data export from the transient simulations were done at the same rate; thus, there were no systematic deviations between the data sets. The agreement between the simulations and the experiments is considered good.

Multiple fragmentation
The binary breakage assumption is a broadly applied simplification for breakup rate modeling [11][12][13]. Only limited attempts have been made to formulate higher order breakage kernels [71,73]. However, experimental measurements [22,23,30] have indicated that the probability of multiple fragmentations is significant, and this probability should be integrated into the breakup rate formulation. Owing to the importance of this phenomenon during the fluid particle breakup process, this sub-section reports on how high-resolved simulations enhance our understanding of this intricate phenomenon. Here we report the total number of fragments i.e. all fragments formed as the results of breakup of the mother droplet. This is relevant measure since we seek to confirm the resolution used in simulations is sufficient to capture also the smallest droplets formed during the breakup event.
The number of fragments recognized in experiments is constrained by the resolution and the contrast between the two phases, whereas analyzing the simulation data relies on the indicator function, which makes quantification of the number of fragments accurate. The spatial resolution of the computational domain was made fine enough to detect   all the daughter droplets produced during a breakup. Fig. 4 demonstrates the fragmentation sequences for two and three levels of dynamic mesh refinements on a simulation. This figure aims to demonstrate the required level of mesh refinement for capturing the smallest daughter droplet under the same turbulent field. The results from mesh independence study confirms also that the rate of energy transfer is not changed using the resolutions reported which is also clear from Fig. 4-a and b for time between 0 and 5 ms. The snapshots start to deviate during the last four milliseconds, when the droplet experienced the maximum deformation state that it undergoes through fragmentation. While for the three refinement levels ( Fig. 4-b) at = t 6.0 ms five fragments are distinguishable, Fig. 4-a shows three daughter droplets with two levels of mesh refinements at the same time. In the final stage of the entire breakup, the finer mesh was able to differentiate seven daughter droplets, and the diameter of the smallest was approximately µ 70 m, while two levels of mesh refinement was unable to capture the smallest daughter. In comparison the smallest fragments found in measurement was µ 90 m. A detailed analysis by magnifying the interface confirms it remains intact at each time without formation of excess interface, this leads to well defined interface and conservation of the dispersed phase during the entire breakup event. Results also converge to mesh independent rate of energy transfer, and number of fragments which indicate that the stresses at the interface is not influenced by the numerics. This implies that the spatial resolution of the domain was fine enough to capture such a small droplet. Accurate predictions of satellite fragments are important as they lead to a better prediction of mass and heat transfer rates. Understanding the multiple fragmentation and development of predictive kernels has a potential to improve the quality of PBE-CFD simulations in the future.
The simulations were validated with experimental data for the water-dodecane system [30] by selecting a narrow range of drop sizes < < d 0. 9 1.1 mm. This includes 17 breakup events and the number of fragments formed equals 60. The benefit of using data from a narrow size range and uniform flow is that the drop We number stays in a narrow range which enables the statistics to be calculated also with limited number of observations. In the simulation data drops sizes are < < d 0. 8 1.0 mm and include 8 breakup and the number of fragments formed equals 31. Table 6 compares the mean and standard deviations of the number of fragments obtained using experiments and simulations. The simulation results are close to the measurements, which confirms the accuracy of the methodology. Both experiments and simulations show that the probability of binary breakup is low, and at  each breakup, the most likely outcome is multiple fragmentation with the average of 3.5 for experiment and 3.9 for simulation. The mean number of fragments for the experiments was somewhat smaller than for the simulations which can be explained by limitations on contrast and resolution for experiments that possibly limits the identification of the smallest fragments.
The energy barrier plot verifies that the curve derivative is a more important variable for the number of fragments than the maximum energy increase. An example in Fig. 3 has been discussed for simulations with three and seven daughter fragments. It has been shown that the faster the transfer rate of kinetic energy from turbulent structures to the droplet the more daughter droplets are created. The statistical analysis also establishes this fact by virtue of the correlation coefficient. The correlation coefficient (Pearson correlation coefficient) between the number of fragments and the maximum interfacial energy was +0.62, while this number increased to +0.85 for the number of fragments and the derivative of the energy barrier plot (maximum of dE dt). This means that larger deformations that lead to more fragments occur due to the interaction with more intense or closely located turbulent vortices. While some of the more recent breakup kernels include integration over the energy spectra, this new information could potentially be included in new kernels to go beyond the limitation of binary breakage that the vast majority breakup models are based on.

Specific breakup rate
The final analysis of simulation data comprised validating breakup rates with experiments reported by Andersson and Andersson [30]. The advantage of an explicit comparison with direct measurements is that they are not the solution of an inverse population balance problem. The common method in indirect measurements of breakup rate is that the breakup frequency is obtained by examining the flux of particles in different size bins of a prescribed distribution. Thus, the ad-hoc assumptions of the daughter size distribution and number of fragments introduce uncertainties in the breakup rates and prevents accurate validations. However, the direct measurements of breakup rates do not require an a priori factor and this prevents the introduction of numerical uncertainty. The specific breakup rate (breakup frequency) is given by the number of breakups per time unit. Calculation of the specific breakup rate thereby requires analyses of all the droplets, including the ones that did not break during the observation time. During the observation time for the mother droplet of = d 1.0 mm 0 , the simulations yielded five breakup and seven non-breakup events. The presence of stable droplets confirm that the simulation method did not produce spontaneous burst breakup due to inconsistencies in the initialization.  Table 7. The agreement is considered good, and the small difference can be explained by the fact that the experimental measurements were performed for a range of droplet diameters with an average diameter of 1 mm, while for the simulations this average value was only applied. This could explain the observed differences between the simulations and experiment. To clarify the level of agreement between the stochastic simulations and the experiments, the table also shows the predictions of the commonly used models as well as our recent model for the entire turbulent energy spectrum. The two models (Andersson and Andersson [15] and Karimi and Andersson [68]) that account for the increase in surface energy during the breakup result in closer predictions to the experiment, whereas the other two overestimate the breakup rate. The breakup rate prediction using the stochastic simulation method verified that this approach did not produce a sudden burst of droplets, and stochastic predictions are also able to capture the physical phenomena observed experimentally during the breakup, such as the necklace breakup reported by Maaβ et al. [22,27].

Conclusions
Stochastic droplet breakup in multiphase turbulent systems was successfully simulated using highly-resolved CFD simulations to study the dynamics of breakup. Under diverse turbulence circumstances of a known turbulent dissipation rate different breakup events were investigated. An LES-VOF model with the dynamic mesh refinement approach was used to model the single droplet breakup. The numerical methodology was validated with experiments and showed very good agreement with respect to deformation time, number of fragments, and breakup rates. Highly-resolved simulations, therefore, offer a valuable alternative for validating theoretical models. Furthermore, the importance of the deformation time was demonstrated through comparisons of the results from simulations, experiments, and theoretical calculations. The comparisons suggest that the application of deformation time can provide physical limits for breakup rate predictions. This could explain the commonly seen monotonic increase in breakup rate models.
The results in this work shows that binary breakup is less likely than multiple fragmentation and on average 3-4 fragments is formed at each breakup. Higher order fragmentation should therefore be incorporated into future breakup kernels. The energy barrier plots were assessed to monitor the level of deformation during the breakup process. It was found that the intensity of the interaction between the droplet and the vortex (i.e. the slope of the curve) highly correlated with the number of fragments formed. This suggests that the kinetics of the process should be integrated in conjunction with the energy criterion. Previous experimental measurements have estimated a 34% increase in interfacial energy during the breakup process; however, the simulations showed that the increase should be modified to 50% based on the accurate 3D representation of the deformed droplets. The findings on the dynamic characteristics of interactions between droplets and turbulent structures should be incorporated in the formulation of breakup. This requires more high resolution simulations and validations to derive mathematical formulations for implementation in CFD-PBE approach. The results presented in this article offer a proof-of-concept for studying droplet breakup and developing breakup kernels based on simulation data. It should be noted that the simulations come with large computational costs, which can prevent simulations in larger geometries and more complex systems such as breakup in stirred tank reactors, and rotorstator mixers. There can also be numerical issues in accounting for tangential stress due to Marangoni effects.   [12] 63.0 Model by Luo and Svendsen [11] 29.8 Model by Andersson and Andersson [15] 4.1 Model by Karimi and Andersson [68] 3.8