Enstrophy transport rates determine the Kolmogorov-Hinze scale in turbulent fragmentation of droplets

PHYSICAL REVIEW FLUIDS 10, 044301 (2025) Enstrophy transport rates determine the Kolmogorov-Hinze scale in turbulent fragmentation of droplets Mahdi Saeedipour *and Simon Schneiderbauer Department of Particulate Flow Modelling, Johannes Kepler University, 4040, Linz, Austria (Received 25 November 2024; accepted 8 April 2025; published 28 April 2025) This study investigates the physical connection between turbulence statistics and the fragmentation cascade of droplets in incompressible nondecaying homogeneous isotropic turbulence. Using interface-resolved direct numerical simulation and based on the concept of vorticity transport across turbulent length scales, we analyze the spectral rate of enstrophy generation by two major production mechanisms: vortex stretching and surface tension. We demonstrate that these spectral rates always equate at a characteristic length scale, implying a balance between disruptive and resistant mechanisms in turbulent fragmentation. At this length scale, the droplet-laden energy cascade shows pivoting from its single-phase counterpart, and the size distribution of drops yields a drastic slope change, indicating where turbulent fragmentation stops. This observation holds for different Taylorscale Reynolds numbers and provides a more deterministic interpretation of the largest stable droplet in turbulence, also known as the Kolmogorov-Hinze (K-H) scale. It is also observed that by increasing the surface tension coefficient, and consequently decreasing the Taylor-scale Weber number toward unity, the K-H scale is determined mainly by the spectral rate of surface tension contribution. The present study carries a direct implication for the size prediction of drop and/or bubble-laden turbulent flows. DOI: 10.1103/PhysRevFluids.10.044301 I. INTRODUCTION Fragmentation and emulsification are multiphase processes that generate dispersed droplets and bubbles recurring in countless real-life applications from fuel atomization and bubble plumes to food emulsions and oil transport. In most of these applications, the multiphase flow becomes turbulent, and the turbulence, in turn, predominantly affects the fragmentation process. Therefore, an in-depth understanding of turbulent fragmentation and the ability to estimate and control the characteristic size of drops or bubbles are essential. The seminal works of Kolmogorov [1] and Hinze [2]on characterizing the size of maximum stable droplets in turbulence, known as Kolmogorov-Hinze (KH) theory, provide a meaningful estimation of the dropand/or bubble sizes in turbulent multiphase flows. However, the theory was developed based on certain assumptions, such as homogeneous isotropic turbulence (HIT) with no viscosity and density ratios between the dispersed and carrier phases [3], and rather dilute conditions. Thus, quantification of the underlying mechanisms of turbulence-interface interactions toward a general theory has been the subject of experimental and computational turbulence research. Recently, Crialesi-Esposito et al. [4] carried out an extensive *Contact author: [email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4. 0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. 2469-990X/2025/10(4)/044301(14) 044301-1 Published by the American Physical Society MAHDI SAEEDIPOUR AND SIMON SCHNEIDERBAUER overview of previous attempts, seeking a connection between turbulence statistics and fragmentation cascades. Various direct numerical simulation (DNS) studies have revealed the scale-dependent nature of surface tension in HIT laden with drops, bubbles, or interfaces: it consumes energy at large scales of turbulence because of fragmentation while it reverts the energy to the small scales through coalescence and surface minimization effects. This entails a deviation in the turbulence energy cascade of interfacial HIT compared with the single-phase flows, and the pivoting scale is linked to the K-H scale [5–8]. The drops smaller than the K-H scale are dominated by coalescence and follow a d−3/2power law in their size distribution [9], whereas the drops larger than this scale tend to break up more and form a sharper slope of d−10/3in their size distribution spectrum [10]. These power laws are observed in various DNS studies [11–13] and confirm the balance between the breakup and coalescence in statistically stationary conditions. Even though the outcome of turbulent fragmentation is governed by the turbulence-interface interactions and interfacial physics, relying on power laws to identify the boundary between coalescence-dominated and breakup-dominated scales implies a statistical nature of the K-H scales, while a deterministic interpretation would be more useful in practice. Therefore, the interpretation of size distribution remains an important active topic in the literature. Recently, Crialesi-Esposito et al. [14] have come up with a more general interpretation for the K-H scale beyond the pivoting point in the energy or drop size spectrum. They analyzed the scale-by-scale energy fluxes in a nondecaying drop-laden HIT at a Taylor-scale Reynolds number of Reλ=140. They concluded that the K-H scale is the scale at which the net contribution of surface tension to the energy spectrum is zero. This definition is used as the basis for proposing an extension to Hinze’s correlation. In our previous work [8], we performed spectral analysis of the enstrophy transport across the scales in a decaying interfacial HIT with Reλ=74. Similar to Crialesi-Esposito et al. [14], we identified a length scale at which the surface tension contribution to the spectral rate of enstrophy changes sign, indicating the change between destruction and generation of enstrophy. Nevertheless, we also identified another characteristic length scale corresponding to the size of eddies at which the rates of vortex stretching and surface tension become equal. This length scale matches the pivoting wave number in energy spectra, and is the same length scale where the drop size distribution changes slope, thus recovering the K-H scale by its definition. This observation was not reported in Crialesi-Esposito et al. [14], which could be due to their higher range of surface tension coefficient, which leads to smaller Weber numbers, reducing the probability of fragmentation events [15]. Besides, only the slope change in the probability density function (PDF) of drops was used as the reference for the interpretation. On the other hand, our previous study focused on a decaying HIT that limits conclusions on the statistically stationary droplet-laden turbulence. In the present study, we further investigate the validity of this observation by extending our enstrophy-based analysis to the nondecaying droplet-laden HIT with higher surface tension coefficients. We explore the connection between enstrophy rates and both energy spectrum and size distribution to identify the K-H scale. II. METHODOLOGY A. Spectral enstrophy transport equation The vorticity transport equation (ωi=ijk∂jUk) is derived by taking the curl of the momentum equation for the one-fluid formulation. The transport equation for vorticity in an incompressible interfacial flow of two immiscible fluids whose interface is tracked by a volume fraction scalar function (α) reads ∂ωi ∂t+Uk ∂ ∂xk ωi=ωk ∂Ui ∂xk+ijk ρ ∂2τkl ∂xj∂xl−ijk ρ2 ∂ρ ∂xj ∂τkl ∂xl+ijk ρ2 ∂ρ ∂xj ∂p ∂xk+ijk∂jσ ρkniδs,(1) 044301-2 ENSTROPHY TRANSPORT RATES DETERMINE THE … where Uiis the mixture velocity vector, pis the pressure, and τkl =μ(∂Uk ∂xl+∂Ul ∂xk) is the viscous stress tensor. The one-fluid density and viscosity are determined by ρ=αρ1+(1 −α)ρ2and μ=αμ1+(1 −α)μ2with ρand μbeing the density and viscosity of each fluid indicated by subscripts. In this volume of fluid-based formulation, the surface tension force is computed by the Continuous Surface Force method [16]. Accordingly, σis the surface tension coefficient, k is the interface curvature, niis the interface normal vector, and δs≡√(∂α/∂xi)(∂α/∂xi)isthe mathematical delta function; thus, niδs=∂α/∂xi. The first two terms on the right-hand side of the Eq. (1) are the vortex stretching (VS) and viscous dissipation (D) terms, similar to single-phase flows, and the latter reduces to ν∂2ωj ∂xj∂xj, for viscosity-matched setups with νbeing kinematic viscosity [8]. The third and fourth terms on the right-hand side are the vorticity production stemming from the misalignment of the density gradient vector with viscous stresses (V) and pressure gradient (B), respectively, and both vanish with a density ratio of unity. Finally, the last term is the curl of surface tension force (σ), indicating the vorticity production by the existence of the fluid-fluid interface. Due to the alignment of density and volume fraction gradients, and given the fact that the curl of a gradient of a scalar is zero, this term reduces to ijkσ ρ ∂k ∂xj ∂α ∂xk, which is solely dependent on interface geometrical changes, that is, derivatives of α. The present analysis is based on the transport of enstrophy (E=ωiωi) across the scales and its connection to the turbulent energy and/or enstrophy spectrum. For this purpose, the spectral enstrophy transport equation is derived by the complex dot product of the Fourier-transformed vorticity vector with Eq. (1) in Fourier space, that reads DE(κ) Dt =β(κ,t),(2) where the contribution of each vorticity generation and/or destruction mechanism across the wave numbers is defined as β(κ,t) with β∈{VS,D,V,B,σ}, and is computed by β(κ,t)= κ<κ<κ+1 R{ˆω i(κ,t)ˆ fβi(κ,t)},(3) where fβirepresents each term on the right-hand side of Eq. (1), “ˆ” denotes the Fourier transform, “” is the complex conjugate, and R{} is the real part. This scale-by-scale analysis is carried out based on the shell-integration process for the spherical shells with a radius of κ[17]. From a physical point of view, represents the spectral contribution of enstrophy generation and/or destruction by each vorticity mechanism and has the physical unit of [ T−3] implying the rate of enstrophy generation associated with each eddy size in turbulence. For a surface tension-only interfacial HIT, the relevant transport rates are VS,D, and σ, which could be reconstructed and analyzed in an a posteriori manner based on the interface-resolved simulation data [18]. In Saeedipour [8], we plotted VS versus σ, herein referred to as -plot, to analyze the turbulence-interface interactions. As mentioned earlier, such plots could disclose certain characteristic wave numbers in interfacial turbulence: κ+ swhere σ∼ =0 and κ+ cwhere VS ∼ =σ.The former is equivalent to the K-H wave number reported by Crialesi-Esposito et al. [14], while the latter matches the definition of the K-H scale in Saeedipour [8]. A similar strategy is followed in the present study to analyze the statistically stationary droplet-laden HIT and pursue a general conclusion about the K-H scale. B. Fully resolved simulation of nondecaying interfacial HIT To investigate the connection between the K-H scale and the -plot, a nondecaying HIT flow laden with interfacial structures should be established. Such a setup can be achieved by implanting a sheet or several large drops in a nondecaying HIT box [4]. Although this entails a transient period until both forced turbulence and the fragmentation process reach a statistically stationary 044301-3 MAHDI SAEEDIPOUR AND SIMON SCHNEIDERBAUER state, it can be assumed that, after a sufficiently long time, the fragmentation outcome will be independent of the initial configuration. Because what matters in a sustained interfacial HIT is the fulfillment of the homogeneity assumption, that is, the drops or bubbles should be dispersed uniformly all over the domain. With this assumption, similar to previous works [8,18,19], we consider a triply periodic cubic box with a length of L=2πm, which is initially filled with a thin sheet of fluid-1 in the center of the box surrounded by fluid-2. In this surface tension-only setup, no density and viscosity contrast is assumed between the fluids (ρ1=ρ2=1kg/m3,μ1= μ2=0.003 Pa ·s). First, we define a reference condition in which the surface tension coefficient and the sheet thickness are set to σ=0.00855 N/m and δ=L/20, respectively. The latter ensures a total volume fraction of =5% for the dispersed phase after the forced turbulence fully fragments the sheet. It must be noted that the domain-averaged turbulent kinetic energy (K) and dissipation rate (ε)are used to compute the Taylor microscale (λ=√10 νK/ε), and consequently, Reλ=ρ√2/3Kλ/μ and Weλ=ρ(2/3K)λ/σ . The interfacial HIT box gradually reaches the statistically stationary state at two different pairs of (Reλ,Weλ): (51, 6) and (121, 45). Hereafter, these two sustained fragmenting systems are referred to as REF-1 and REF-2 cases, respectively. We carried out interface-resolved simulations using a geometric volume of fluid (VOF) method [20] which solves the Navier-Stokes equations together with the VOF equation on a grid number of 5123 that is sufficient for the fully resolved simulation at this range of Reλ. The simulations were run for about 30 eddy turnover time (τe=CeK/ε with Ce=0.15). The VOF simulation procedure and numerical schemes are similar to our previous DNS studies [8,21], and are not repeated here. The HIT is enforced using the method of Eswaran and Pope [22], allowing efficient forcing in parallel finite-volume simulations. More recently, this method has been successfully applied to force HIT of gas-particle flows [23]. First, the forcing method is expressed in Fourier space, where the forcing term is denoted by ˆ fi(κ,t) with the wave number vector κ=(κ1,κ 2,κ 3). Forcing is only applied to the smallest wave numbers, that is, κ<2.8, while the origin κi=0 is omitted. The complex Fourier coefficients ˆ fi(κ,t) are computed from independent Uhlenbeck-Ornstein processes for each of the forced wave numbers. This random process is given by ˆ bi(κ,t+t)=ˆ bi(κ,t)+ei(κ,t)2S2t TL1/2 ,(4) where tdenotes the numerical time step and ei(κ,t) is a complex random number sampled from a standardized Gaussian distribution with zero mean and variance unity at each time step for each wave number. Thus, in total, six scalar random numbers are required (for real and imaginary parts) for each time step and wave number. TLrepresents the characteristic timescale of the random process and S2its variance. Interestingly, in the case of TL=tthe random process gets uncorrelated in time [22]. Finally, the volume forcing term ˆ fi(κ,t) is determined by the orthogonal projection of ˆ bi, which reads ˆ fi(κ,t+t)=ˆ bi(κ,t+t)−κi κjˆ bji(κ,t+t) κjκj ,(5) which ensures zero divergence of the forcing. The forcing parameters TLand ∗=S2TLare set to achieve certain combinations Reλand Weλ. For more details, the reader is referred to Eswaran and Pope [22] and Chouippe and Uhlmann [23]. It remains to discuss the implementation of the complex forcing term, ˆ fi, which is expressed in Fourier space, into a finite-volume solver. Since a complete Fourier transform in massively parallel computations is costly, we evaluated the Fourier series directly on each processor. This procedure does not include any processor communication except distributing ˆ fiacross the pro044301-4 ENSTROPHY TRANSPORT RATES DETERMINE THE … FIG. 1. Instantaneous snapshots of interfacial structures in nondecaying HIT for different cases at t/τe> 20 visualized by iso-contours of α=0.5. (a) REF-1: Reλ=51, Weλ=6 (b) REF-2: Reγ=121, Weλ=45. cessors and evaluating the Fourier series in the local memory of each parallel process is rather cost-efficient [23]. III. ANALYSIS AND DISCUSSION A. Macroscopic characteristics of the nondecaying HIT Fully resolved simulation of the nondecaying HIT exhibits different macroscopic flow characteristics between the investigated cases. Figure 1displays the instantaneous three-dimensional snapshots of the interfacial flow in the cases REF-1 and REF-2, after the statistically stationery state is reached, demonstrating the size, shape, and fragmentation pattern of the drops. Additionally, Fig. 2presents two consecutive pairs of zoomed-in snapshots of arbitrary coalescing and fragmenting drops from REF-1, which demonstrate a drastic change in vorticity magnitudes during these interfacial phenomena. This underlines the importance of vorticity transport in the context of fragmentation and coalescence in turbulence and explains the motivation for conducting a spectral analysis of vorticity and enstrophy mechanisms. From a physical viewpoint, the fragmentation processes such as atomization and emulsification increase the total interfacial area (A=Vδsdv) of the dispersed phase while the total volume remains the same. On the other hand, coalescence or the act of surface minimization decreases the total interfacial area. Thus, the temporal variation of the interfacial area, that is, dA/dt is a macroscopic indicator of equilibrium conditions in interfacial systems. Increasing both Reλand Weλreduces drop size and increases fragmentation due to higher inertia, generating a greater number of smaller drops and thus a higher total interfacial area. Analyzing the cases REF-1 and REF-2 shows a transition period until domain-averaged kinetic energy reaches a constant level, as depicted in Fig. 3(a). During this period, the total interfacial area increases sharply due to initial sheet fragmentation. This indicates a positive slope in the temporal evolution of the interfacial area, that is, dA/dt >0. Nevertheless, the positive slope decreases gradually with time and dA/dt →0 for both cases [Fig. 3(b)]. This implies a balance between fragmentation and surface minimization through coalescence. Such a statistically stationary condition allows further spectral analysis of energy, enstrophy, and fragmentation cascades, which were carried out for t/τe>20 in both cases. B. Analysis of the K-H scale The primary objective of this study is to establish a physical connection between the statistics of nondecaying HIT and the fragmentation process through the concept of enstrophy. In the absence 044301-5 MAHDI SAEEDIPOUR AND SIMON SCHNEIDERBAUER FIG. 2. Consecutive pairs of zoomed-in snapshots of exemplary fragmenting and coalescing drops from REF-1. The interfacial structures are colored by ω∗=|ω| |ω| (the magnitude of the local vorticity vector normalized by the domain-averaged vorticity magnitude). The time interval between each pair of interfacial events is 0.5τe. of material property gradients, vortex stretching and surface tension misalignment are the main mechanisms for enstrophy generation and/or destruction. Their spectral rates, VS and σ[Eq. (3)], are compared across the scales to identify characteristic wave numbers associated with the size of the largest stable droplets in sustained turbulence [8]. The results are presented as the -plots, energy spectra, and droplet size distributions for each case separately across the wave number spectrum, which is normalized by the Taylor microscale, that is, κ+=κ κλ. Figures 4(a) and 5(a) depict the -plots for both REF-1 and REF-2 at two different instants of time, confirming a sustained distribution across the scales. The scale-dependent nature of σis evident, changing sign at κ+ sfrom negative (fragmentation-dominated large scales) to positive (coalescence-dominated small scales). Nevertheless, in both REF-1 and REF-2, there are certain wave numbers with σ>0 for which, still, the rate of enstrophy generation by vortex stretching is higher than surface tension. This implies that the eddies in this range are still sufficiently energetic to continue the cascade and promote fragmentation in the form of drop formation. Thus, another characteristic wave number (κ+ c)exists at which both enstrophy rates are equal. Both, κ+ sand κ+ care smaller for REF-1, indicating sustained dispersed characteristics with larger droplets, which is consistent with the overall picture in Fig. 1. Figures 4(b) and 5(b) show the normalized kinetic energy spectra for the same instants and confirm that the energy cascade does not decay with time. They further show the deviation from a single-phase (droplet-free) HIT at large wave numbers with the same forcing parameters. This deviation starts at a wave number identical to κ+ cfor both cases, whereas no specific conclusion can be drawn for κ+ s. This important observation aligns with the energy-based description of turbulent fragmentation: energetic eddies stretch and counteract surface tension to induce fragmentation. Beyond a certain wave number (κ+>κ + c), the eddies are no longer able to form drops, and energy begins to rise as surface tension predominantly reverts energy due to the coalescence and surface minimization. While the -plot and energy spectra correspond at these major characteristic scales, they are both computed from turbulent properties. Thus, further analysis of fragmentation statistics in terms of size distribution is necessary. Figures 4(c) and 5(c) show the number density as the PDF of droplets, computed for the same instants. The distributions converge to a unique PDF over time, 044301-6 ENSTROPHY TRANSPORT RATES DETERMINE THE … (a) ( b ) FIG. 3. Temporal variation of the domain-averaged quantities of cases REF-1 and REF-2: (a) kinetic energy and (b) normalized total interfacial area. indicating sustained fragmentation statistics. The 3/2 and 10/3 power law slopes are also examined, and the wave numbers associated with the slope change are identified. As an essential observation, this wave number is in agreement with κ+ cfrom -plot. It must be noted that determining the slope change may introduce uncertainty, depending on how the PDFs are generated. For certain cases, like REF-2, the 10/3 power law is not maintained for larger droplets, which is consistent with findings in other studies [24]. Nevertheless, the change in power law serves as a common criterion to identify the K-H scale, distinguishing between fragmentation-dominated and coalescence-dominated scales in size distribution [25], and is therefore employed here as well. We also analyzed the average of statistical quantities for REF-2 to ensure that a sustained statistical steady state is reached. Figure 6presents the time averages of both -plots and droplet PDFs, computed from five different instances for t/τe>20, including their corresponding standard deviations. These results reveal that the droplet size distribution as well as the -terms converge to unique values and remain constant over time, which is particularly observed around the wave numbers, where the slope change of the droplet size distribution occurs, as well as where VS ∼ =σ. This demonstrates the reliability of the present cases for drawing general conclusions, despite acknowledging the potential for under-resolved very small droplets due to the universal limitation of the VOF method. The present spectral analysis of investigated cases demonstrates a systematic connection between turbulence and fragmentation via features of -plots and determines the K-H scale as the κ+ c. Nevertheless, this observation is different from the recent work of Crialesi-Esposito et al. [14], who proposed κ+ sas the basis for refining the K-H correlation. Despite similarities in most physical 044301-7 MAHDI SAEEDIPOUR AND SIMON SCHNEIDERBAUER (a) (b) (c) FIG. 4. (a) The -plots for REF-1 at two instants of time: VS (solid curves), σ(dashed curves), and D(thin curves) divided by the domain-averaged enstrophy. The green curves show the results of single-phase simulations. (b) The normalized kinetic energy spectra compared with their single-phase similitudes (dashed curve). (c) Droplet size distribution as a function of normalized equivalent wave number (κ+ eq =2π d·1 κλ) with the power law slopes. The green vertical line displays the PDF slope change as determined by the derivative of the data. In all plots, the normalized characteristic wave numbers κ+ sand κ+ care displayed by vertical gray and red lines, respectively. and simulation parameters, the surface tension coefficients used in Crialesi-Esposito et al. [14]are up to two orders of magnitude larger than in the present study. This affects the Taylor-scale Weber number and has a direct impact on the characteristic length scales in turbulent fragmentation. C. The effect of surface tension coefficient The influence of the surface tension coefficient is examined through Weλ. Only REF-2 is considered for the new simulations with surface tension coefficients 2, 3, 5, 10, and 30 times larger than the initial value, while keeping other physical and simulation parameters constant. Figure 7 displays the -plot energy spectra, and the PDF of the drops for some of the new cases, together with the original simulation. The notable effect of increased surface tension is observed in the drop sizes, which, in turn, affects the -plot by shifting the characteristic wave numbers κ+ sand κ+ c toward larger scales. This indicates a reduction in the range of fragmenting scales (where σ<0) and a broader range of coalescing structures across the spatial spectrum, which is also evident in the energy spectra, showing larger deviation from the single-phase spectrum at the highest σ.In addition, the change in power law slopes becomes more evident in the PDFs, as depicted in Fig. 7(c). This observation is aligned with the findings of Vela-Martín and Avila [15], where they described the turbulent fragmentation in HIT as a memoryless process dependent on the Weber number. As the Weber number decreases, the probability of droplet survival increases, and the fragmentation rate 044301-8 ENSTROPHY TRANSPORT RATES DETERMINE THE … (a) (b) (c) FIG. 5. (a) The -plots for REF-2 at two instants of time: VS (solid curves), σ(dashed curves), and D(thin curves) divided by the domain-averaged enstrophy. The green curves show the results of single-phase simulations. (b) The normalized kinetic energy spectra compared with their single-phase similitudes (dashed curve). (c) Droplet size distribution as a function of normalized equivalent wave number with power law slopes. The green vertical line displays the PDF slope change as determined by the derivative of the data. In all plots, the normalized characteristic wave numbers κ+ sand κ+ care displayed by vertical gray and red lines, respectively. decreases. Figure 8displays the slopes of the size distribution with respect to the equivalent wave numbers, which confirms this observation. In this figure, the slopes are computed by the derivatives of the data shown in Fig. 7(c), that is, log(pdf(d))/ log(κ+ eq ) using a finite difference approach, and the corresponding κ+ sand κ+ cextracted from the -plots are plotted for better comparison. It is evident that for all the cases, the slope change occurs in a range between κ+ sand κ+ c. In the case of lower surface tension, the slope shifts from 3/2to10/3 immediately at a wave number similar to κ+ c, confirming the observation for REF-2. However, increasing the surface tension alters this pattern. For the high surface tension coefficient, the scaling law with the slope of 3/2 still holds for κ+ cand changes to 10/3 in a wave number closer to κ+ s. This observation confirms the significance of the surface tension coefficient in determining the K-H scale from the enstrophy plots. For similarly forced HITs at the same Reλ,alargeWeλresults in a very high fragmentation rate, meaning that VS > σcan occur at almost all small scales until it stops at VS ∼ =σand determines the K-H scale. This can even occur for wave numbers larger than κ+ sdue to the higher frequency of intense eddying motions [15]. But increasing surface tension reduces Weλ, implying that the interface provides a larger resistance to fragmentation and reduces the range of eddies from the inertial subrange of turbulence capable of producing droplets of comparable size. Consequently, fragmentation induced by eddying motions may stop at a size comparable to λ. It is established that the survival probability of these droplets is extremely high, and they could remain stable for a long time [15]. In other words, finding droplet breakup events at scales below where σ∼ =0 becomes 044301-9