category: literaturenote citekey: gugliettadeformationellipsoidaldropletshomogeneous2025 title: Deformation of ellipsoidal droplets in homogeneous and isotropic turbulence authors: "Guglietta, Fabio; Taglienti, Diego; Sbragaglia, Mauro" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.09692v1" zotero_key: 9ZQF9B2E zotero_storage: 9ZQQBX9I collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Guglietta, Fabio"

status: converted

DEFORMATION OF ELLIPSOIDAL DROPLETS IN HOMOGENEOUS AND ISOTROPIC TURBULENCE

Fabio Guglietta

Department of Physics & INFN, University of Rome "Tor Vergata", Via della Ricerca Scientifica 1, 00133, Rome, Italy. fabio.guglietta@roma2.infn.it

Diego Taglienti

Department of Physics & INFN, University of Rome "Tor Vergata", Via della Ricerca Scientifica 1, 00133, Rome, Italy.

Mauro Sbragaglia

Department of Physics & INFN, University of Rome "Tor Vergata", Via della Ricerca Scientifica 1, 00133, Rome, Italy.

ABSTRACT

We study the statistics of deformation of neutrally buoyant droplets in homogeneous isotropic turbulence (HIT), wherein the characteristic droplet size R is smaller than the characteristic Kolmogorov scale η of the turbulent flow. We systematically focus on the characterization of droplet statistics obtained with various phenomenological ellipsoidal models (EMs) – assuming that the droplet preserves the ellipsoidal shape at all times – with the droplet moving as a passive tracer in the turbulent flow. The predictions of the EMs are compared with ground-truth data obtained with three-dimensional fully resolved simulations (FRSs) without any ad-hoc assumption on the droplet shape. Our work helps in elucidating the applicability of the EMs in describing droplet deformation in HIT at changing the capillary number Ca = τσ/τη, weighting the relative importance of the droplet characteristic time τ^σ with respect to the turbulent flow characteristic time τη.

1 Introduction

Turbulent flows laden with droplets or bubbles are widely encountered in several settings, including atmospheric clouds [1], combustion chambers [2], nuclear reactors [3], oceanic transport [4] and ocean spray aerosols [5, 6] just to cite a few examples. In these situations, droplets and bubbles move in a complex fluid dynamics environment and deform due to the action of the outer surrounding flow. Various reviews have been written on the topic [7–15]. The process of deformation depends on how the hydrodynamic and surface forces balance at the droplet interface, which we consider in this paper solely characterized by a homogeneous surface tension σ. Such balance is quantified by the dimensionless capillary number Ca = τσ/τη, where τ^η is a characteristic time of the flow and τ^σ = µR/σ the characteristic droplet time, with R a representative length scale for the droplet [we choose R as the radius at rest, see Fig. 1(a)] and µ the outer flow viscosity. The time τ^σ represents the time it takes for the droplet to relax after deformation or, equivalently, to adapt to an outer perturbation with characteristic time τη. If τ^σ ≪ τ^η (i.e., small Ca), we have a linear regime where small deviations from sphericity are promoted by the outer flow; when τ^σ ∼ τ^η (i.e., large Ca), we have a non-linear response resulting in larger deformations [see Fig. 1(b)]. Furthermore, in the paradigmatic set-up of homogeneous and isotropic turbulence (HIT), the flow is typically generated by energy injection at large scales. Energy is then transported via inertial (non-dissipative) contributions towards the small scales, where the energy spectrum is cut by viscous effects at the Kolmogorov scale η [16]. Therefore, droplet dynamics crucially depends on the ratio between the droplet scale R and the Kolmogorov length scale η: for R ≫ η, inertial effects play a role in droplet dynamics whereas, for R ≪ η, they are negligible and one can describe the coupling between the droplet or bubble and the fluid via the Stokes equations [17]. In this paper, we focus on the latter case.

Figure 1: Panel (a): fully resolved simulations (FRSs) of sub-Kolmogorov droplet dynamics along a homogeneous isotropic turbulence (HIT) trajectory in an incompressible flow. The colorbar indicates the local magnitude of the velocity gradient ∇u along the trajectory. FRSs are performed with the numerical methodology described in Ref. [32], capable of reproducing droplets with a structured triangular mesh (see zoom on the interface in the center-top panel) at the interface undergoing their dynamics in HIT, wherein the characteristic droplet size (R) falls below the characteristic Kolmogorov scale (η) of the HIT. Droplet morphology is quantified via the construction of an equivalent ellipsoid with semiaxes L, B, W (L > W > B). Statistics of L, B , W over an ensemble of HIT trajectories is computed and compared with corresponding ellipsoidal models (EMs) predictions [see Eq. (1) and text for more details], assuming ellipsoidal shapes at all times. Panel (b): illustration of the effect of increasing deformation on droplet morphology.

Due to the complexity of the problem under study, numerical simulations are frequently used to unravel many questions [14]. Different numerical techniques have been employed to investigate the dynamics of sub-Kolmogorov droplets and bubbles, like the boundary integral method [17, 18], the diffuse-interface lattice Boltzmann method [19], approaches based on point-like (i.e., non-deformable) particle assumption [20 –27], point-like particle approaches coupled to subgrid models [28], dynamics of neutrally buoyant sub-Kolmogorov droplets [29–31], hybrid immersed boundary-lattice Boltzmann methods [32]. In some of these studies [28–31], instead of directly considering the Stokes equations for the deformable droplets, it is assumed a reduced order model relying on simpler ordinary differential equations regulating the evolution of a morphological tensor S used to describe the droplet interface [33]. These models are based on the assumption of small deformations (Ca ≪ 1), implying that the droplet retains an ellipsoidal shape at all times [33, 34] [see Fig. 1(b)], hence the name of ellipsoidal models (EMs). As shown in Refs. [28–31], EMs can be used to infer relevant features regarding droplet rheology, such as deformation, orientation, morphology, and breakup. Furthermore, this approach allows for efficient statistical analyses, for instance, when considering a collection of Lagrangian trajectories followed and probed with sub-Kolmogorov droplets. In light of experimental observations for stationary flows [35–37], droplets are found to retain an ellipsoidal shape even for finite to high Ca values, so that one may ask whether extending the use of EMs in the high capillary regime to complex turbulent flows could be effective to reproduce rheological data. The aim of this paper is to address this question by considering the statistics of neutrally buoyant sub-Kolmogorov droplet deformation in HIT, and comparing results from EMs with those of fully resolved simulations (FRSs). Specifically, we extract the information on the velocity gradients probed from direct numerical simulations of HIT [38], and we use it to both numerically integrate the EMs and to perform ab-initio FRSs [32], thereby obtaining quantitative information on the droplet's principal axes of deformation [see Fig. 1] and enabling a direct comparison between EMs and FRSs.

The paper is organized as follows: in Sec. 2, FRSs and EMs are reviewed and described; then in Sec. 3 we present the numerical framework employed to perform the FRSs; results are presented in Sec. 4; finally, conclusions are provided in Sec. 5.

2 Fully Resolved Simulations and Ellipsoidal Models

The dynamics of droplets in turbulence can be investigated using two complementary approaches. On the one hand, fully resolved simulations offer high-fidelity, ab-initio solutions by explicitly resolving the fluid–structure interactions

at the droplet interface. On the other hand, reduced-order models provide an efficient way to approximate droplet deformation based on simplified physics. In the following, we briefly describe both approaches used in this study.

Fully Resolved Simulation (FRSs)

FRSs performed with the immersed boundary-lattice Boltzmann method [39] provide a powerful tool to generate ground-truth data for droplet dynamics in HIT. Such methodology has been recently employed to simulate the dynamics of droplets and soft particles in a variety of setups [32, 40–48]. In particular, for the simulations of droplets in turbulence, we employ the numerical method described in Ref. [32]. In a nutshell, the lattice Boltzmann method is used to simulate the fluid both inside and outside the droplet, and the immersed boundary method couples the fluid with the droplet interface, the latter being represented by a triangular mesh [see Fig. 1(a)] and characterized by a tunable surface tension σ. Numerical simulations are performed with a mesh-smoothing technique to regularize the triangular mesh [32]. We refer the reader interested in technical details to Ref. [32].

Ellipsoidal models (EMs)

EMs are reduced-order models designed to describe the evolution of droplet shape, which is represented by a secondorder morphological tensor S. These models offer two main advantages. First, the evolution equation for S is governed solely by the velocity gradient tensor ∇u and its interactions with S, without requiring full knowledge of the flow field. Second, the resulting ordinary differential equation can be integrated efficiently using standard numerical methods, such as Runge–Kutta schemes, instead of performing ab-initio numerical simulations. It is important to stress that EMs assume the droplet to be always ellipsoidal: while this is a good approximation for small degrees of deformation, it does not hold anymore when the deformation is large, e.g., close to the breakup [see Fig. 1(b)]. Information about droplet deformation can be obtained from the eigenvalues of the tensor S, which correspond to the squares of the ellipsoid's semiaxes. Taking the square roots of these eigenvalues yields the lengths of the droplet's principal semiaxes, which we denote as the length (L), breadth (B), and width (W), conventionally ordered such that L > W > B (see Fig. 1).

2.1 Maffettone and Minale (MM) model.

One of the most celebrated EMs is that proposed by Maffettone and Minale (MM) [34]. The model equation describing the time evolution of S is:

$$\frac{dS}{dt} - (\Omega \cdot S - S \cdot \Omega) = -\frac{f1}{\tau{\sigma}} [S - g(S)I] + f_2 (E \cdot S + S \cdot E) , \qquad (1)$$

where

$$\Omega = \frac{1}{2} (\nabla u - \nabla u^T), \qquad E = \frac{1}{2} (\nabla u + \nabla u^T)$$ (2)

are the vorticity and the rate-of-strain tensors, respectively. The non-linear function g (S) = ³III^S II^S is used to enforce volume conservation, with II^S = 2 tr(S) ² − tr S 2 i and III^S = det(S) being the second and third invariant of the tensor S, respectively [34]. The two parameters f¹ = f1(λ) and f² = f2(λ) are set by the viscosity ratio between the inner and outer fluid, λ, and are retrieved in the framework of the perturbation theory for small deformations (Ca ≪ 1) applied to Stokes equations [34, 49]. Indeed, expanding the morphological tensor S up to second order (excluded) in capillary number: S = IR² + S (1)Ca + O(Ca² ), and plugging this result in Eq. (1), we obtain [34]:

$$\frac{d\mathbf{S}}{dt} = -\frac{f1}{\tau{\sigma}} \left[ \mathbf{S} - \mathbf{I} R^2 \right] + 2R^2 f_2 \mathbf{E} , \qquad (3)$$

where terms of O(Ca² ) have been neglected. We will refer to the model descibed by Eq. (3) as linear theory (LT). By comparing Eq. (3) with the analytical results in the perturbation theory aforementioned [34, 49], Maffettone and Minale obtained [34]:

$$f_1^{\text{MM}}(\lambda) = \frac{40(\lambda + 1)}{(2\lambda + 3)(19\lambda + 16)}, \qquad f_2^{\text{MM}}(\lambda) = \frac{5}{2\lambda + 3}.$$ (4)

We remark that the model in Eq. (1) exploits the linear nature of the Stokes equations to superimpose the deformation effects given by the flow with the counteracting surface tension effects; the function g(S) cannot be derived from linear theory, and is thus manually inserted into the model to enforce volume conservation by requiring the determinant of S to be constant [34]. It is worth noting that the LT model has a trivial volume conservation that reads g(S) = R² . This conservation is essential according to the primitive equations, where the divergence of the velocity field is zero.

2.2 Modified Maffettone and Minale (MM)

Maffettone and Minale noted that the parameter f² = f MM ² defined in Eq. (4) underestimates the deformation and overestimates the values of the critical capillary number Cacr (i.e., the capillary number at which the droplet shows breakup), which can be computed as a function of f¹ and f² [34]. For such reasons, they propose a different choice of f2, allowing the latter to depend also on Ca (we will refer to this model extension as MM model):

$$f_2^{\overline{\text{MM}}}(\lambda, \text{Ca}) = \frac{5}{2\lambda + 3} + \frac{3\text{Ca}^2}{2 + 6\text{Ca}^{2+\delta}} \frac{1}{1 + \epsilon\lambda^2},$$ (5)

with δ and ϵ being small positive numbers. However, Maffettone and Minale stated that if λ and Ca are not large, then δ = ϵ = 0 (we adopt the same choice in this work). Note that f ² → f MM in the limit Ca → 0.

2.3 Extended Maffettone and Minale (EMM).

We already showed that Maffettone and Minale [34] matched the parameters f¹ = f MM and f² = f MM ² with the analytical results obtained in the perturbation theory [8, 34, 49], thus making them intrinsically built on the assumption of small deformation. Indeed, the MM model underestimates the deformation, in particular for high values of capillary number Ca, even for steady flows [34, 41]. In our previous work [41], we proposed an extension of the MM model (EMM) by providing a polynomial expression in Ca for f¹ and f2:

$$f{1,2}^{\text{EMM}}(\lambda, \text{Ca}) = \sum{i} a_i^{(1,2)}(\lambda) \text{Ca}^i,$$ (6)

where the coefficients a (1,2) ⁱ were fitted with data from numerical simulations. Values of a (1,2) i for λ = 1 are reported here 1 .

2.4 Quasi-static (QS) Approximation.

In the limit of small deformations (i.e., Ca → 0), the droplet dynamics is expected to be quasi-static (QS), dS/dt ≈ 0, due to the fact that the characteristic droplet time is much smaller than the characteristic time of the outer flow [18]. From Eq. (1), in the limit Ca → 0, we have g (S) → R² , (Ω · S − S · Ω) → 0 and (E · S + S · E) → 2R2E, so that QS approximation implies:

$$S_{\rm QS} = IR^2 + \frac{2f2\tau\sigma}{f_1}R^2E$$ (7)

In other words, in the limit Ca → 0, there is no need to integrate Eq. (1), since the eigenvalues of the morphological tensor S can be reconstructed from the eigenvalues Λ1,2,³ of the rate-of-strain tensor E (see Fig. 2).

3 Numerical Setup

To assess the accuracy of the EMs introduced in Sec. 2 within the context of HIT, we compare their predictions against FRSs under identical flow conditions. In both cases, droplet dynamics are driven by the same set of time-dependent velocity gradients ∇u(t), extracted from direct numerical simulations of HIT [38]. The focus is on quantifying the statistical features of droplet deformation across a range of capillary numbers and evaluating the impact of different closure strategies for the EM model coefficients f¹ and f2. The HIT flow field is characterized by a Taylor-scale Reynolds number Re ∼ 130 and is drawn from the TURB-Lagr open-source database [38], which provides over 3×10⁵ particle trajectories. Each trajectory lasts for roughly an eddy turnover time τeddy ∼ 200τ^η [38]. For completeness, in Fig. 2, we provide details on the PDF of the eigenvalues Λ1,2,³ of the rate-of-strain tensor E extracted from the HIT trajectories that we utilized [50–57]. Consistently with previous literature results, we have positively/negatively skewed PDFs for the maximum/minimum eigenvalues Λ1,3, and a positively skewed PDF for the intermediate eigenvalue Λ2, this latter fact being crucially related to a non-zero strain production of the turbulent flow [58]. Following earlier investigations, we treat neutrally buoyant droplets with a unitary viscosity ratio (λ = 1) and a unitary density ratio between the droplet fluid and the outer fluid. Droplets are treated as passive tracers, following the Lagrangian trajectories of an outer turbulent flow. The outer flow of the FRS domain is adapted to accommodate the turbulent gradient matrix ∂iu^j (t) at all times under the assumption that we are dealing with a sub-Kolmogorov droplet, hence the outer flow

¹The coefficients a (1,2) (λ = 1) for i ∈ [0, 5] are reported: a ⁰ = 0.457, a ¹ = 0.235, a ² = −4.546, a ³ = 20.536, a ⁴ = −34.798; and a (2) ⁰ = 1, and a ⁱ = 0 for i ∈ [1, 5].

Figure 2: PDFs of the eigenvalues Λ1,2,³ of the rate of strain tensor E of the HIT turbulent trajectories used for the analysis of droplet deformation. Eigenvalues are normalized with the norm ||E|| of the rate-of-strain tensor E.

matches a linear flow constructed from the turbulent gradients [17, 29]. The analysis for droplet deformation is conducted for three different capillary numbers: Ca = 0.025, Ca = 0.15, and Ca = 0.4. For each Ca studied, simulations for one thousand independent trajectories are performed. Concerning the FRSs, the fluid domain consists of a cubic box with size L^x = L^y = L^z = 10R; the droplet is simulated with a structured mesh composed of 20480 triangles (see Fig. 1). The radius of the droplet at rest is chosen as R/η = 0.472, which guarantees a fair convergence to the sub-Kolmogorov dynamical equations [32]. Since FRSs are performed within a finite fluid domain, a threshold is imposed on the droplet elongation: when L/R exceeds 5, the corresponding trajectory is excluded from the statistics and replaced by a new one (this effect is particularly prominent at high capillary numbers, such as Ca = 0.4). Such a condition mimics the onset of breakup in finite domains, consistent with observations from Ref. [17]. All simulations are run using an in-house GPU-accelerated code on Nvidia A100 GPUs, each requiring roughly three hours of computation time. To generate EM data, we use the same set of turbulent trajectories from the HIT database. Instead of resolving the fluid–structure interaction, we integrate the tensorial evolution Eq. (1), where the time-dependent tensors E(t) and Ω(t) are dynamically updated along each trajectory. Integration is performed with a standard explicit Runge–Kutta scheme [59]. We consider four EMs: the linear theory (LT) model [Eq. (3)], the original MM model [Eqs. (1),(4)], the MM model with modified f² [Eq. (5)], and the extended model EMM [Eq. (6)]. These comparisons aim to evaluate how each model captures the deformation statistics, particularly across different regimes of Ca, and clarify the role of nonlinearities and volume conservation in improving model fidelity.

4 Results

In Fig. 3, we analyze a representative turbulent trajectory and report the time evolution of the normalized droplet semiaxes, L/R, B/R, and W/R [see also Fig. 1(a)], for the smallest and largest investigated capillary numbers Ca. Time is normalized with the eddy turnover time τeddy [38]. Data from FRSs are reported with markers, while data from EMs are reported with solid lines. Different colors represent different semiaxes. At Ca = 0.025, all EMs exhibit excellent agreement with the FRSs. In this regime, deviations from sphericity are minimal (i.e., S → IR² ), and even the LT prediction [see Eq. (3)] provides accurate results. The small Ca limit implies that the droplet relaxation time is much shorter than the characteristic timescale of the flow, making the coupling between Ω and S, which appears on the l.h.s. of Eq. (1) and is known to scale as O(Ca² ) [34, 49], negligible. Notably, the LT model can be derived analytically from the Stokes equations in the asymptotic limit Ca → 0 [49]. Thus, the analysis at small Ca also serves as a benchmark for the numerical framework, validating both the accuracy of the FRSs and extending the validation of Ref. [32] to include all droplet semiaxes. For Ca = 0.4, data from FRSs reveal (see t/τeddy ≈ 0.62, in correspondence of the gray-shaded region) a correlation between enhanced peaks in L/R and sudden compressions on the orthogonal semiaxes (W and B). Such behavior was also observed in Ref. [29], where the authors identified a transition first from an oblate shape (i.e., L/R ∼ W/R and B/R < 1) to a prolate shape (i.e., L/R ≫ 1 and W/R ∼ B/R < 1), and then the droplet recovers the spherical shape (i.e., L/R ∼ W/R ∼ B/R ∼ 1). It is worth noting that, in the specific example of Fig. 3, among the EMs considered, the MM and EMM models best capture the high values of L/R. In particular, MM tends to slightly underestimate the FRS data, while EMM slightly overestimates it. This trend will also

Figure 3: Comparison between fully resolved simulations (FRSs) and different ellipsoidal models (EMs) for two values of capillary number, Ca = 0.025 (top panels) and Ca = 0.4 (bottom panels), along a representative turbulent trajectory. Panels (a–d) correspond to Ca = 0.025, showing the normalized semiaxes L/R, B/R, and W/R as functions of normalized time t/τeddy, as predicted by: (a) the linear theory (LT), (b) the Maffettone–Minale (MM) model, (c) the modified MM model with f² given in Eq. (5) (MM), and (d) the extended MM model (EMM). Panels (e–h) show the same analysis for Ca = 0.4, where droplet deformation is significantly more pronounced. In each panel, symbols represent FRSs data (blue circles for L, red triangles for B, and yellow squares for W), while solid lines show the predictions from the corresponding EM. The vertical gray band highlights a region of peak deformation (see Sec. 4).

be observed in the statistical analysis reported below (see later Fig. 7). We also stress that when the droplets undergo strong deformations, the ellipsoidal assumption does not hold anymore, as can be seen in Fig. 1, panel (b): the reported droplet shapes come from the FRSs, while the orange dashed lines serve as a reference of the ellipsoidal shape. For Ca = 0.4, variations in droplet length L are not strongly pronounced in the LT model nor in the MM and MM models, if seen in comparison to FRS data; among all EMs analyzed, data from EMM model display more pronounced variations in L, a fact that we can attribute on the way the coefficients f EMM and f EMM 2 are computed (see Sec. 2). Overall, apart from rare events during which the droplet undergoes quick, strong deformations, there are time lapses over the eddy turnover time analyzed where the MM, MM, and EMM predictions show a good agreement with FRSs. We remark, however, that data reported in Fig. 3 just refer to a single turbulent trajectory; hence, it is dutiful to assess the robustness of these findings via some statistical analysis on the various observables on an ensemble of droplet trajectories.

Figure 4: Probability distribution functions (PDFs) of normalized droplet semiaxes L/R, B/R, and W/R for a statistical ensemble of turbulent trajectories at Ca = 0.025. The statistics is collected from an ensemble of 1000 independent droplet trajectories in HIT. Panels (a–d) compare the fully resolved simulation (FRS) data (symbols) with the predictions of various ellipsoidal models (EMs) (solid lines): (a) the linear theory (LT), (b) the Maffettone–Minale (MM) model, (c) the modified MM model with f² given in Eq. (5) (MM), and (d) the extended MM model (EMM). Colors and markers denote the three principal semiaxes: blue circles and lines for L, red triangles and lines for B, and yellow squares and lines for W. All PDFs are shown on a logarithmic scale to emphasize tail behavior.

Figure 5: Probability distribution functions (PDFs) of normalized droplet semiaxes L/R, B/R, and W/R for a statistical ensemble of turbulent trajectories at Ca = 0.025. The statistics is collected from an ensemble of 1000 independent droplet trajectories in HIT. Comparison between the fully resolved simulation (FRS) data (symbols) with the predictions of the Maffettone and Minale (MM) model (solid lines) and the quasi-static (QS, dashed lines) approximation (see Eq. 7). Colors and markers denote the three principal semiaxes: blue circles and lines for L, red triangles and lines for B, and yellow squares and lines for W. All PDFs are shown on a logarithmic scale to emphasize tail behavior.

Figure 6: Probability distribution functions (PDFs) of normalized droplet semiaxes L/R, B/R, and W/R for a statistical ensemble of turbulent trajectories at Ca = 0.15. The statistics is collected from an ensemble of 1000 independent droplet trajectories in HIT. Panels (a–d) compare the fully resolved simulation (FRS) data (symbols) with the predictions of various ellipsoidal models (EMs) (solid lines): (a) the linear theory (LT), (b) the Maffettone–Minale (MM) model, (c) the modified MM model with f² given in Eq. (5) (MM), and (d) the extended MM model (EMM). Colors and markers denote the three principal semiaxes: blue circles and lines for L, red triangles and lines for B, and yellow squares and lines for W. All PDFs are shown on a logarithmic scale to emphasize tail behavior.

Figure 7: Same as Fig. 6, for Ca = 0.4.

We then study the probability distribution functions (PDFs) of the ellipsoid normalized semiaxes L/R, B/R, and W/R, collecting data for 1000 independent trajectories at changing Ca. In Fig. 4, we report data for Ca = 0.025. Apart from rare events in the tails, we observe a good matching between the PDF extracted from the FRSs and those extracted from the EMs; thus, the agreement observed in Fig. 3 for Ca = 0.025 and for a single trajectory is rather robust throughout

the ensemble of turbulent trajectories. The tails of the distributions for L/R and B/R predicted by the LT model exhibit larger discrepancies compared to those from the other EMs. As discussed in Sec. 2.1, the LT model can be derived from the MM model [see Eq. (1)] when the capillary number approaches zero (Ca → 0). Moreover, we remark that the LT model has a trivial volume conservation in the context of the EMs, i.e., g(S) = IR² [see Eq. (3)]. Although data in Fig. 4 corresponds to a relatively small value of Ca (i.e., Ca = 0.025), the observed deviations suggest that the effects of a finite capillary number start to emerge. We then expect that the LT model will be the first among the considered EMs to lose predictive accuracy as Ca increases. As already stated in Sec. 2, in the limit of small deformations (Ca → 0), the droplet dynamics is expected to be quasi-static [see Eq. (7)]. This behavior is quantitatively confirmed in Fig. 5, where we report the PDFs of the normalized semiaxes for the FRSs and the MM model (previously shown in Fig. 4), along with the predictions obtained from the morphological tensor in Eq. (7), namely LQS/R, WQS/R, and BQS/R. The excellent agreement among these datasets supports the validity of the QS approximation in this low-Ca regime. Since in the QS limit, the morphological tensor describing the deviation from sphericity is the strain matrix E multiplied by a factor [i.e., SQS − IR² ∼ E, see Eq. (7)], then the droplet deformation statistics inherits the structural features of the turbulent strain field. In other words, the morphological tensor SQS becomes "enslaved" to the strain matrix E. In particular, the PDFs of the normalized semiaxes reflect the positive skewness associated with the intermediate eigenvalue Λ² of the rate-of-strain tensor E, as highlighted in Fig. 2. The fact that this statistical feature reappears in the droplet deformation analysis underscores a meaningful connection between turbulent strain dynamics and the behavior of sub-Kolmogorov deformable droplets. It is worth noting that the tails of the QS model for B/R and L/R show similar discrepancies to those of the LT model at Ca = 0.025 [see Fig. 4, panel (a)]. Indeed, both the QS and LT models are valid in the Ca → 0 limit.

In Figs. 6 and 7 we extend our investigation on the PDFs at higher Ca regime showing results for Ca = 0.15 and Ca = 0.4, respectively. Regarding the statistics of the normalized length L/R, Fig. 6 (i.e., Ca = 0.15) shows that the LT model [panel (a)] significantly underestimates the probability of large deformations compared to FRS. This limitation becomes even more pronounced at Ca = 0.4 [Fig. 7, panel (a)], where the MM model [panel (b)] also begins to deviate substantially in the tail, failing to capture extreme elongations. In contrast, both the MM [panel (c)] and EMM [panel (d)] models follow the FRS data more closely, especially in reproducing the long-tailed distribution of L/R, demonstrating improved robustness at higher deformation regimes. For the normalized breadth B/R, the discrepancies between EMs and FRS data are less severe than those observed for L/R: while the LT model still shows visible deviations at both Ca = 0.15 and Ca = 0.4, the other EMs maintain overall a better alignment with the FRS statistics. Concerning the normalized width W/R, the evolution of the PDF shape as Ca increases provides additional insight. As remarked earlier, at lower Ca, the distribution shows a slight positive skewness, consistent with the QS regime and the inherent asymmetry of the turbulent strain field (see Fig. 2). As Ca increases, however, the skewness of W/R reverses and becomes negative. This transition indicates the growing influence of large deformation events associated with configurations approaching breakup, where droplets acquire cigar-like topologies elongated in the direction of maximum extension. Such behavior has been reported in previous works [29], where negative skewness in the distribution of minor semiaxis has been linked to precursors of breakup events. At Ca = 0.4, the negative skewness of W/R becomes prominent in the FRS data, yet it is poorly captured by the LT model. The MM model offers a partial improvement, while the MM and EMM models provide a closer approximation of the observed asymmetry. Importantly, at Ca = 0.4, FRS data reveal that droplet elongation can reach values large enough to exceed the domain boundaries. For this reason, in the case of the FRS, a threshold on droplet deformation must be introduced: without it, excessively deformed droplets would exceed the boundaries of the finite-size numerical domain. This issue does not arise for EMs, which are based on the integration of ordinary differential equations and do not involve spatial confinement, allowing, in principle, for unbounded deformation. To ensure a fair comparison, we applied the same deformation threshold to both FRS and EM data. Once this threshold is enforced, we observe that FRS trajectories exceed it more frequently than those predicted by the MM model. However, this is not generally true for all EMs, as the EMM model exhibits threshold-crossing events at a frequency comparable to that of FRSs. This enhanced stability shown by the MM model is consistent with classical findings on droplet breakup thresholds. For example, in simple shear flow with viscosity ratio λ = 1, breakup is observed at Cacr ≈ 0.41, whereas the MM model predicts no breakup, i.e., Cacr = ∞ [34]. Similarly, in axisymmetric extensional flow, numerical simulations report Cacr ≈ 0.06 [60], while the MM model yields a significantly higher threshold, Cacr ≈ 0.228 [34]. These results confirm that the MM model systematically underestimates the onset of breakup, yielding artificially stable droplet dynamics under turbulent forcing.

5 Conclusions

We have studied the statistics of sub-Kolmogorov droplet deformation in homogeneous isotropic turbulence (HIT). With reference to fully resolved simulations (FRSs) performed with the numerical methodology described in Ref. [32], the primary focus of our analysis was on the assessment of the quality of the statistics extracted from some reduced order models that are built on the assumption that the droplet is always ellipsoidal – and for this reason they are referred

as ellipsoidal models (EMs). Among them, we started from the Maffettone and Minale (MM) model, and then we considered other models that can be derived from the MM one, such as: the linear theory (LT) model, which is retrieved in the assumption of Ca → 0; the model proposed again by Maffettone and Minale [34] but with a different choice of f² in order to make the model able to predict better the critical capillary number (which we referred as MM model); the extension of the Maffettone and Minale (EMM) model proposed by Taglienti et al. [41] which can catch the steady-state deformation of a droplet under stationary shear flow for high values of Ca. The comparison of LT with other EMs has been instructive to qualify the importance of the main constitutive elements of the MM, MM, and EMM models (e.g., non-linearity and non-trivial volume conservation). Our analysis highlights that the LT model performs well only in the low capillary number regime, where deformation remains small, and it fails to capture the large deformation statistics observed at higher Ca. The MM model extends the predictive range to moderate capillary numbers but still underestimates rare events and cannot reproduce the negative skewness that develops in the distribution of the normalized width (W/R). Such behavior may be attributed to the MM model's tendency to underestimate droplet deformation, which in turn leads to an overestimation of the critical capillary number. In contrast, both MM and EMM show improved agreement with FRS data across the full range of Ca explored: they reproduce the skewed and heavy-tailed distributions of all three semiaxes, including the transition from positive to negative skewness in W/R, associated with the emergence of cigar-like shapes near breakup [29].

This study opens interesting perspectives for further fruitful research in the future. For example, one could extend our analysis to cases with other Reynolds numbers and/or different viscosity ratios λ, which in the present study have been kept fixed. Moreover, our study can also be considered as a starting point for the exploration of the dynamics of soft particles exhibiting additional complexities at the interface, like for example surface viscosity [61, 62] or elasticity [63, 64]; these interfacial properties can be introduced either in the numerical methodology that we used to get the ground truth data [42, 43, 45] and also in the EMs [61, 62, 65]. All in all, these studies could provide valuable insights to make an assessment of the validity of EMs. Beyond all the physical information that one could gather, we also remark that FRSs may become computationally expensive; hence, having the possibility to generate in shorter times ensembles of data with quantitatively validated EMs is surely welcomed for statistical analysis and/or AI studies [66 –70] on the dynamics of soft particles with complex interfaces in complex flows.

6 Acknowledgments

We wish to acknowledge Fabio Bonaccorso, Luca Biferale and Maurizio Carbone for support and useful discussions. This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 882340). This work was supported by the Italian Ministry of University and Research (MUR) under the FARE program (No. R2045J8XAW), project "Smart-HEART". MS acknowledges the support of the National Center for HPC, Big Data and Quantum Computing, Project CN_00000013 - CUP E83C22003230001, Mission 4 Component 2 Investment 1.4, funded by the European Union - NextGenerationEU. Funding from Tor Vergata University project AI4HEART is also acknowledged.

References

¹R. A. Shaw, "Particle-turbulence interactions in atmospheric clouds", Annual Review of Fluid Mechanics 35, 183–227 (2003).
²M Carl, T. Behrendt, C Fleing, M Frodermann, J Heinze, C Hassa, U Meier, D Wolff-Gaßmann, S Hohmann, and N Zarzalis, "Experimental and numerical investigation of a planar combustor sector at realistic operating conditions", J. Eng. Gas Turbines Power 123, 810–816 (2001).
³A. Poullikkas, "Effects of two-phase liquid-gas flow on the performance of nuclear reactor cooling pumps", Progress in Nuclear Energy 42, 3–10 (2003).
4 I. D. Nissanka and P. D. Yapa, "Calculation of oil droplet size distribution in ocean oil spills: a review", Marine Pollution Bulletin 135, 723–734 (2018).
⁵E. Villermaux, X. Wang, and L. Deike, "Bubbles spray aerosols: certitudes and mysteries", PNAS nexus 1, pgac261 (2022).
⁶L. Deike, "Mass transfer at the ocean–atmosphere interface: the role of wave breaking, droplets, and bubbles", Annual Review of Fluid Mechanics 54, 191–224 (2022).
7 J. Harper, "The motion of bubbles and drops through liquids", Advances in Applied Mechanics 12, 59–129 (1972).
8 J. Rallison, "The deformation of small viscous drops and bubbles in shear flows", Annual Review of Fluid Mechanics 16, 45–66 (1984).
⁹H. A. Stone, "Dynamics of drop deformation and breakup in viscous fluids", Annual Review of Fluid Mechanics 26, 65–102 (1994).
¹⁰P. Fischer and P. Erni, "Emulsion drops in external flow fields—the role of liquid interfaces", Current Opinion in Colloid & Interface Science 12, 196–205 (2007).
¹¹J. Magnaudet and I. Eames, "The motion of high-reynolds-number bubbles in inhomogeneous flows", Annual Review of Fluid Mechanics 32, 659–708 (2000).
¹²V. Cristini and Y.-C. Tan, "Theory and numerical simulation of droplet dynamics in complex flows—a review", Lab on a Chip 4, 257–264 (2004).
¹³D. Lohse, "Bubble puzzles: from fundamentals to applications", Physical review fluids 3, 110504 (2018).
¹⁴S. Elghobashi, "Direct numerical simulation of turbulent flows laden with droplets or bubbles", Annual Review of Fluid Mechanics 51, 217–244 (2019).
¹⁵R. Ni, "Deformation and breakup of bubbles and drops in turbulence", Annual Review of Fluid Mechanics 56, 319–347 (2024).
¹⁶U. Frisch, Turbulence: the legacy of a.n. kolmogorov (Cambridge University Press, Nov. 1995).
¹⁷V. Cristini, J. Bławzdziewicz, M. Loewenberg, and L. R. Collins, "Breakup in stochastic stokes flows: sub-Kolmogorov drops in isotropic turbulence", Journal of Fluid Mechanics 492, 231–250 (2003).
¹⁸V Cristini, S Guido, A Alfani, J Bławzdziewicz, and M Loewenberg, "Drop breakup and fragment size distribution in shear flow", Journal of Rheology 47, 1283–1298 (2003).
¹⁹F. Milan, L. Biferale, M. Sbragaglia, and F. Toschi, "Sub-Kolmogorov droplet dynamics in isotropic turbulence using a multiscale lattice Boltzmann scheme", Journal of Computational Science 45, 101178 (2020).
²⁰S Balachandar and J. K. Eaton, "Turbulent dispersed multiphase flow", Annual review of fluid mechanics 42, 111–133 (2010).
²¹J. G. Kuerten and A. Vreman, "Effect of droplet interaction on droplet-laden turbulent channel flow", Physics of Fluids 27 (2015).
²²E. Russo, J. G. Kuerten, C. W. van der Geld, and B. J. Geurts, "Water droplet condensation and evaporation in turbulent channel flow", Journal of fluid mechanics 749, 666–700 (2014).
²³A. Ferrante and S. Elghobashi, "On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles", Journal of Fluid Mechanics 503, 345–355 (2004).
²⁴I. M. Mazzitelli, D. Lohse, and F. Toschi, "On the relevance of the lift force in bubbly turbulence", Journal of Fluid Mechanics 488, 283–313 (2003).
²⁵P. C. Le Clercq and J. Bellan, "Direct numerical simulation of gaseous mixing layers laden with multicomponent-liquid drops: liquid-specific effects", Journal of Fluid Mechanics 533, 57–94 (2005).
²⁶F. Mashayek, "Droplet–turbulence interactions in low-mach-number homogeneous shear two-phase flows", Journal of Fluid Mechanics 367, 163–203 (1998).
²⁷R. S. Miller and J Bellan, "Direct numerical simulation of a confined three-dimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden stream", Journal of Fluid Mechanics 384, 293–338 (1999).
²⁸V. Spandan, R. Verzicco, and D. Lohse, "Deformable ellipsoidal bubbles in taylor-couette flow with enhanced euler-lagrangian tracking", Physical Review Fluids 2, 104304 (2017).
²⁹L. Biferale, C. Meneveau, and R. Verzicco, "Deformation statistics of sub-Kolmogorov-scale ellipsoidal neutrally buoyant drops in isotropic turbulence", Journal of Fluid Mechanics 754, 184–207 (2014).
³⁰V. Spandan, D. Lohse, and R. Verzicco, "Deformation and orientation statistics of neutrally buoyant sub-Kolmogorov ellipsoidal droplets in turbulent Taylor-Couette flow", Journal of Fluid Mechanics 809, 480–501 (2016).
³¹S. S. Ray and D. Vincenzi, "Droplets in isotropic turbulence: deformation and breakup statistics", Journal of Fluid Mechanics 852, 313–328 (2018).
³²D. Taglienti, F. Guglietta, and M. Sbragaglia, "Droplet dynamics in homogeneous isotropic turbulence with the immersed boundary–lattice boltzmann method", Physical Review E 110 (2024).
³³M. Minale, "Models for the deformation of a single ellipsoidal drop: a review", Rheologica acta 49, 789–806 (2010).
³⁴P. Maffettone and M. Minale, "Equation of change for ellipsoidal drops in viscous flow", Journal of Non-Newtonian Fluid Mechanics 78, 227–241 (1998).
³⁵S. Guido and M. Villone, "Three-dimensional shape of a drop under simple shear flow", Journal of Rheology 42, 395–415 (1998).
³⁶A. Vananroye, P. Van Puyvelde, and P. Moldenaers, "Effect of confinement on droplet breakup in sheared emulsions", Langmuir 22, 3972–3974 (2006).
³⁷K. Feigl, D. Megias Alguacil, P. Fischer, and E. Windhab, "Simulation and experiments of droplet deformation and orientation in simple shear flow with surfactants", Chemical Engineering Science 62, 3242–3258 (2007).
³⁸L. Biferale, F. Bonaccorso, M. Buzzicotti, and C. Calascibetta, "TURB-Lagr. A database of 3D Lagrangian trajectories in homogeneous and isotropic turbulence", arXiv:2303.08662 (2023).
³⁹T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, The lattice Boltzmann method - principles and practice (Springer Cham, Oct. 2016).
⁴⁰P. Li and J. Zhang, "A finite difference method with subsampling for immersed boundary simulations of the capsule dynamics with viscoelastic membranes", International Journal for Numerical Methods in Biomedical Engineering 35, e3200 (2019).
⁴¹D. Taglienti, F. Guglietta, and M. Sbragaglia, "Reduced model for droplet dynamics in shear flows at finite capillary numbers", Physical Review Fluids 8, 013603 (2023).
⁴²F Pelusi, F Guglietta, M Sega, O Aouane, and J Harting, "A sharp interface approach for wetting dynamics of coated droplets and soft particles", Physics of Fluids 35, 082126 (2023).
⁴³F. Guglietta and F. Pelusi, "Analytical prediction for the steady-state behavior of a confined drop with interface viscosity under shear flow", Phys. Rev. Fluids 9, 103603 (2024).
⁴⁴E. Bellantoni, F. Guglietta, F. Pelusi, M. Desbrun, K. Um, M. Nicolaou, N. Savva, and M. Sbragaglia, "Immersed boundary-lattice boltzmann method for wetting problems", arXiv preprint arXiv:2503.20605 (2025).
⁴⁵F. Guglietta, M. Behr, L. Biferale, G. Falcucci, and M. Sbragaglia, "On the effects of membrane viscosity on transient red blood cell dynamics", Soft Matter 16, 6191–6205 (2020).
⁴⁶F. Guglietta, M. Behr, G. Falcucci, and M. Sbragaglia, "Loading and relaxation dynamics of a red blood cell", Soft Matter 17, 5978–5990 (2021).
⁴⁷F. Guglietta, M. Behr, L. Biferale, G. Falcucci, and M. Sbragaglia, "Lattice Boltzmann simulations on the tumbling to tank-treading transition: effects of membrane viscosity", Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 379, 20200395 (2021).
⁴⁸F. Guglietta, F. Pelusi, M. Sega, O. Aouane, and J. Harting, "Suspensions of viscoelastic capsules: effect of membrane viscosity on transient dynamics", Journal of Fluid Mechanics 971, A13 (2023).
⁴⁹J. Rallison, "Note on the time-dependent deformation of a viscous drop which is almost spherical", Journal of Fluid Mechanics 98, 625–633 (1980).
⁵⁰R. M. Kerr, "Histograms of helicity and strain in numerical turbulence", Physical review letters 59, 783 (1987).
⁵¹A Tsinober, E Kit, and T Dracos, "Experimental investigation of the field of velocity gradients in turbulent flows", Journal of Fluid Mechanics 242, 169–192 (1992).
⁵²T. S. Lund and M. M. Rogers, "An improved measure of strain state probability in turbulent flows", Physics of Fluids 6, 1838–1847 (1994).
⁵³L. K. Su and W. J. Dahm, "Scalar imaging velocimetry measurements of the velocity gradient tensor field in turbulent flows. ii. experimental results", Physics of Fluids 8, 1883–1906 (1996).
⁵⁴A. Tsinober, An informal conceptual introduction to turbulence (Springer, 2009).
⁵⁵J. Tom, M. Carbone, and A. D. Bragg, "Exploring the turbulent velocity gradients at different scales from the perspective of the strain-rate eigenframe", Journal of Fluid Mechanics 910, A24 (2021).
⁵⁶D. Buaria, A. Pumir, and E. Bodenschatz, "Generation of intense dissipation in high reynolds number turbulence", Philosophical Transactions of the Royal Society A 380, 20210088 (2022).
⁵⁷M. Carbone and M. Wilczek, "Asymptotic predictions on the velocity gradient statistics in low-reynolds-number random flows: onset of skewness, intermittency and alignments", Journal of Fluid Mechanics 986, A25 (2024).
⁵⁸B. Lüthi, A. Tsinober, and W. Kinzelbach, "Lagrangian measurement of vorticity dynamics in turbulent flow", Journal of Fluid mechanics 528, 87–118 (2005).
⁵⁹J. Butcher, "A history of runge-kutta methods", Applied Numerical Mathematics 20, 247–260 (1996).
⁶⁰N. Singh and V. Narsimhan, "Numerical investigation of the effect of surface viscosity on droplet breakup and relaxation under axisymmetric extensional flow", Journal of Fluid Mechanics 946, A24 (2022).
⁶¹D. Barthes-Biesel and H. Sgaier, "Role of membrane viscosity in the orientation and deformation of a spherical capsule suspended in shear flow", Journal of Fluid Mechanics 160, 119–135 (1985).
⁶²R. W. Flumerfelt, "Effects of dynamic interfacial properties on drop deformation and orientation in shear and extensional flow fields", Journal of Colloid and Interface Science 76, 330–349 (1980).
⁶³D. Barthes-Biesel and J. Rallison, "The time-dependent deformation of a capsule freely suspended in a linear shear flow", Journal of Fluid Mechanics 113, 251–267 (1981).
⁶⁴D. Barthès-Biesel, "Motion and Deformation of Elastic Capsules and Vesicles in Flow", Annual Review of Fluid Mechanics 48, 25–52 (2016).
⁶⁵V. Narsimhan, "Shape and rheology of droplets with viscous surface moduli", Journal of Fluid Mechanics 862, 385–420 (2019).
⁶⁶D. A. de Aguiar, H. L. França, and C. M. Oishi, "Predicting energy budgets in droplet dynamics: a recurrent neural network approach", International Journal for Numerical Methods in Fluids 97, 854–873 (2025).
⁶⁷Y. Jia, X. Cheng, W. Song, Y. Zhou, and H. Zhao, "A droplet state prediction method for inkjet printing under small sample conditions based on the two-stage tradaboost. r2 algorithm", AIP Advances 15, https://doi.org/10. 1063/5.0246942 (2025).
⁶⁸L. Au-Yeung and P. A. Tsai, "Predicting impact outcomes and maximum spreading of drop impact on heated nanostructures using machine learning", Langmuir 39, 18327–18341 (2023).
⁶⁹R. Choudhury, D. M. T. Nasir, D. R. Kaiser, R. Ghimire, D. A. Aliyu, D. B. Sohani, D. J. Atanbori, D. R. Rathaur, and D. R. Mishra, "Uncovering the effect of physical conditions and surface roughness on the maximum spreading factor of impinging droplets using a supervised artificial neural network model", Industrial & Engineering Chemistry Research 62, 21208–21221 (2023).
⁷⁰M. Tembely, D. C. Vadillo, A. Dolatabadi, and A. Soucemarianadin, "A machine learning approach for predicting the maximum spreading factor of droplets upon impact on surfaces with various wettabilities", Processes 10, 1141 (2022).

gugliettadeformationellipsoidaldropletshomogeneous2025.md 58 KB Historia Czysty