The classical droplet oscillation theory considers the shape oscillations of an incompressible liquid drop of radius R ₀, viscosity μ _l , and density ρ _l . The drop is considered to be close to the spherical stable condition due to the surface tension σ. The surrounding fluid is considered to be initially at rest, with density ρ _g and viscosity μ _g . Both fluids are considered isothermal.

Any shape that is star-shaped (radially convex set) can be described in a spherical coordinate reference (see Figure 1). Since the considered droplet is close to a sphere, this framework is used. Any function based on the spherical coordinates (θ, φ), where 0 < θ < π is the polar angle (colatitude) and 0 < φ < 2π is the azimuthal angle (longitude), can be expressed as an expansion of spherical harmonic functions. The distance from the center of the reference frame to the droplet interface can, thus, be described by the following equation: R θ , φ , t = R 0 ∑ m = 0 ∞ ∑ l = − m m a m , l t Y m l θ , φ , (1) where a _m,l (t) corresponds to the amplitude of the deformation expressed as a spherical harmonic coefficient and Y m l is the spherical harmonic function. Finally, m and l correspond to the spherical harmonic degree and order, respectively, with m ≥ 0 and −m ≤ l ≤ m. Degrees (which are referred to in the literature as modes) describe the number of lobes in the polar coordinate. Order refers to the orientation of these lobes in the reference frame system. The spherical harmonic function is defined as follows: H m l θ , φ = 2 m + 1 m − l ! m + l ! P m l c o s ⁡ θ e i l φ , (2) with P m l ( cos ⁡ θ ) corresponding to the Legendre polynomials, m ∈ N and l ∈ [−m, m]. Then, the real form of the spherical harmonics is given by the following equation: Y m l θ , φ = i 2 H m l − − 1 m H m − l if  m < 0 H m 0 if  m = 0 i 2 H m − l + − 1 m H m l if  m > 0 . (3)

The equation of motion of the free surface can be obtained from the linearization of the Navier–Stokes equations [28]. For the particular case, in which the surrounding fluid is of negligible dynamical effects, the evolution in time of the coefficients a _m,l are solutions of the independent second-order differential equation [24]: a ̈ m , l + 2 β m 0 a ̇ m , l + ω m 0 2 a m , l = 0 , (4) where ^⋅ corresponds to the time derivative, β _m0 is the damping coefficient of mode m, and ω _m0 is the angular frequency associated with mode m. It was shown that for this particular case, the angular frequency and damping rate are given by the following equations: ω m 0 = m m + 1 m − 1 m + 2 σ R 0 3 m ρ g + m + 1 ρ l , (5) β m 0 = 2 m + 1 m − 1 μ l ρ l R 0 2 . (6)

For the studied case of liquid–liquid droplets, the solution given by Prosperetti [24] is considered. In the limit of weak viscous effects, this equation can then be simplified by estimations obtained from an asymptotic development [29], which draws the following solution for the angular frequency: ω m = ω m 0 1 − 2 m + 1 2 ρ ̂ μ ̂ 2 2 ω m 0 ρ l μ l R 0 m ρ ̂ + m + 1 1 + ρ ̂ μ ̂ ⏟ A . (7) Then, the expression for the damping rate is defined as follows: β m = ω m 0 A − 2 A 2 + 2 + m 2 m 2 − 1 + m + 2 μ ̂ − m − 1 ρ ̂ μ ̂ + 2 m m + 2 ρ ̂ μ ̂ 2 2 m ρ ̂ + m + 1 1 + ρ ̂ μ ̂ 2 μ l R 0 2 ρ l , (8) where ρ ̂ = ρ g / ρ l is the density ratio and μ ̂ = μ g / μ l . It should be noted that the subscripts m are used to describe the solution of these equations for the associated mode m. In the present paper, we consider the case where the density ratio between two phases is one.

The zeroth coefficient a _0,0 evolves in time in order to ensure mass conservation [30]. According to these authors, the first mode should not report any oscillation. Since Eq. 4 is independent, the present theory does not consider any coupling between the coefficients.

A general review of the different theories can be found in [23]. More recent developments can be found in [10,31].

For the sake of clarity, in the current paper, we use the frequency f ₂ = 2πω ₂ computed from Eq. 7 to non-dimensionalize the time, considering that mode 2 is expected to dominate the deformation.

The in-house high-performance computing code Archer is used to perform all DNS presented in this work. Archer was one of the first codes undertaking the simulation of liquid-jet atomization under a realistic diesel injection configuration [32]. It solves on a staggered Cartesian mesh the one-fluid formulation of the incompressible Navier–Stokes equation. ∇ ⋅ u ⃗ = 0 , (9) ∂ t ρ u ⃗ + ∇ ⃗ ⋅ ρ u ⃗ ⊗ u ⃗ = − ∇ ⃗ p + ∇ ⃗ ⋅ 2 μ D + f ⃗ + σ κ δ s n ⃗ , (10) where ρ s the density, p is the pressure field, μ is the dynamic viscosity, D is the strain rate tensor, f ⃗ is a source term, σ is the surface tension, n ⃗ is the unit normal vector to the interface, κ is its mean curvature, and δ _s is the Dirac function characterizing the location of the interface. For solving Eq. 10, the convective term is written in the conservative form and solved using the improved Rudman’s technique [33] presented in [34]. The latter allows mass and momentum to be transported consistently, thereby enabling flows with large liquid/gas density ratios to be simulated accurately. The viscosity term is computed following the method presented by [35]. To ensure the incompressibility of the velocity field, a Poisson equation is solved. The latter accounts for the surface tension force and is solved using a multi-grid preconditioned conjugate gradient algorithm (MGCG) [36] coupled with the ghost fluid method [37].

For transporting the interface, a CLSVOF method is used, in which the level-set function accurately describes the geometric features of the interface (its normal and curvature) and the VOF function ensures mass conservation. The density is calculated from the VOF, ψ, as ρ = ρ _l ψ + ρ _g (1 − ψ), and the dynamic viscosity (μ _l or μ _g ) depends on the sign of the level-set function. In cells containing both a liquid and gas phase, a specific treatment is performed to evaluate the dynamic viscosity, following the procedure of [35]. The code has been validated and tested in a variety of atomization conditions, such as high-density ratios and varying Weber and Reynolds numbers [32,34]. The case of an oscillating droplet has been previously validated for the Archer code [25,26], finding good agreement with the linear theory.

To produce a turbulent flow in a physical space-based code, linear forcing has been applied. This forcing was first developed for a single-phase flow by [38] and has been adapted to two-phase flows with particles [39] and interfacial flows [40,41]. This subsequent implementation is used in the present research. The linear forcing consists in introducing a volume force in the Navier–Stokes equation proportional to the velocity fluctuations, f ⃗ = A u ⃗ ′ . Here, u ⃗ ′ = u ⃗ − u ⃗ is computed at the end of the prediction step, and the parameter A is updated each time step in order to obtain a target total energy. Limited by the number of nodes, the present simulations are performed at a fixed turbulent Reynolds number, R e λ = ρ l u ′ λ μ l ≃ 45 , where λ is the Taylor length scale.

The impact of the forcing scheme is currently under study [42] and will be presented in a future communication.

A database of independent droplets immersed in a homogeneous and isotropic turbulent background flow was generated as a part of a previous study ([43,44]—forthcoming). We consider the liquid–liquid interaction to experience no mixing or mass exchange. The Weber number, defined as W e = ρ l u ′ 2 2 R 0 σ , remained constant over all cases (We = 0.9). This Weber number was set to ensure the droplet would oscillate without experiencing an immediate breakup. The viscosity ratio (μ _l /μ _g ) and density ratio (ρ _l /ρ _g ) are set to 1. In order to focus on the interactions between the interface and the turbulence, it was decided to work with the same density and the same dynamic viscosity for both fluids [13]. Turbulence forcing can introduce additional problems, in particular when the density ratio is large [18]; therefore, by setting this parameter to unity, the influence of the density difference on the carrier turbulent flow is avoided. The numerical simulation was performed over a 64³ grid. A ratio between the droplet radius and domain size of R 0 L = 0.233 and an Ohnesorge number ( O h = μ l σ ρ l R 0 ) of 0.00465 were considered. Each droplet is generated in the following way: first, a single-phase isotropic and homogeneous turbulent flow is generated in the domain and run until it converges. After this, a spherical solid body is located inside the domain with the surrounding fluid to capture the non-slip condition. The background flow is then left to evolve around the sphere, and once convergence has been reached respecting the no-slip condition, the spherical solid body is changed for a liquid bulk. The numerical simulation was set to continue until the droplet experienced breakup. For the results presented in this article, over 220 droplets were analyzed. The database generation is fully presented in Deberne et al. (2023—forthcoming).

To study the shape evolution of a drop immersed in a turbulent flow, we introduce the spherical harmonic functions Y m l ( θ , φ ) (Eq. 2). The representation of the spherical harmonic function as a double series is a generalized Fourier series. Therefore, if normalized, it represents a complete orthonormal system of functions for all degrees m and orders l. Due to this property, we can then compute the spherical harmonic coefficients of a function with the integral: a m , l = ∫ Ω R θ , φ Y m l θ , φ d Ω , (11) where dΩ corresponds to the differential surface area on the unit sphere sin θdθ dϕ.

The values of R(θ, φ) are obtained through a loop that computes the values of the radius across the interface (Figure 2A), following the spherical coordinates reference frame over an array of angles, whose size is given by the maximum number of modes, M, we wish to evaluate. More specifically, the array (as shown in Figure 2B) will be of shape 2(M + 1) × 2(M + 1). The coefficients are then computed from this array with the help of pyshtools [45]. The numerical error for computing the radius R(θ, φ) is less than 0.1%. Similarly, the numerical error for computing the coefficients a _m,l is less than 0.1% for m ≤ 6. A detailed overview of the spherical harmonic expansion framework applied here has been determined as part of this research [27].

For the studied case, we consider the computation of the spherical harmonic modes for m ≤ 6 as we observe a rapid decrease of the amplitude of the coefficient a _m,l with m. Since we aim to use this framework to characterize the deformation on the droplet surface, we are interested in reducing as much as possible the effect of mode 1 (a _1,l ≈ 0) since it describes the translatory motion of the droplet. The center of mass was opted as the center of space after considering different approaches [27]. For the results shown here, 90% of the computed mode 1 coefficients are below 0.05. An exception to this rule occurs when the droplet is subject to high-amplitude deformation.

Our goal is to relate the spherical harmonic expansion with the surface deformation. We recall that the total energy stored at the interface is proportional to the total surface area η Ω 2 . Parseval’s theorem, also known as Rayleigh’s energy theorem, relates the integral of the square of a function to the sum of the square of its transform. Due to the orthogonality properties of the spherical harmonic functions, this can be written as follows: 1 4 π ∫ Ω R 2 θ , ϕ d Ω = η Ω 2 , (12) where η Ω 2 = ∑ m = 0 ∞ ∑ l = − m m a m , l 2 . (13)

The power spectrum a m 2 is computed by the following equation: a m 2 = ∑ l = − m m a m , l 2 . (14)

Since the coefficients a _m,l of the same degree m are associated with the same oscillatory frequency, the power spectrum provides the energy contained in each mode and relates it to the total deformation. This parameter will always be positive by definition and independent of the orientation of the reference frame.

The spherical harmonic coefficients are dependent on the reference frame of the system. This means that the results from the expansion are directly affected by the center of the coordinate system and the reference frame in which we define the angles θ and φ. Due to the nature of the surrounding fluid in this study, the droplet will not remain static, but it will experience translation and rotation. In the present, an origin at the center of mass of the body is considered, with a reference frame aligned to the Cartesian mesh, except for the analysis conducted, as shown in section 3.4.

In order to have a reference for the case studied in this article, a series of simulations without turbulence are performed with similar conditions and numerical resolution as for the turbulent case. The selected configuration is to investigate the oscillations of each isolated mode in a fluid initially at rest. Both the drop and the surrounding fluid are of equivalent characteristic values as those used for the generation of the database. Each simulation corresponds to the small-amplitude perturbation of degrees m = 2–6. The spherical harmonic coefficients a _m,0 were chosen in accordance to the highest amplitude of the coefficients found in turbulent flow deformation (presented later in section 3.2). The reference frame was selected so that we would follow the axisymmetric oscillations of the liquid drop (i.e., θ is the angle respecting to the axisymmetric axis.). The temporal evolution of the coefficients and their initial shape are given in Figure 3.

In the simulation, common results are observed. First, the oscillatory behavior is clearly seen. The amplitude reduces in time as expected by the theory. Second, it is noticeable that the droplet tends to spend more time in the prolate state for mode 2 as a known non-linear effect [46]. This can be seen comparing the time spent on positive (i.e., prolate) and negative (i.e., oblate) a _m,0 when considering even modes (i.e. 2, 4, 6). For mode 2 oscillations, it is computed that the droplet spends 56.75% of the total time in the prolate state. This phenomenon has been analyzed in detail in the works of Foote [47]; Basaran [8]; Meradji et al. [48]; Zrnić et al. [10]. For modes 4 and 6, the oscillations appear to be more periodic, not having such a difference when comparing prolate to the oblate state. Finally, even if only a pure mode is initially perturbed, coupling between modes is observed. This is particularly seen for even modes (see Figures 3A, C, E), as suggested by Foote [47] and Alonso [49]. The coupling between orders seems negligible in all the cases, except order l = 4 that appears for modes m = 4–6. Indeed, this order is aligned with the Cartesian grid, and thus, numerical approximations affect this order.

The frequency and damping of each simulation are given in Table 1. The values of ω _m and β _m were computed using a Levenberg–Marquardt fitting algorithm. The covariance of the fitted results was below 0.5% for all cases. As predicted by the theory [6,28], these values increase with the degree. In other words, the deformations with smaller amplitude have larger frequency and are dissipated faster. The computed values differ from the theoretical predictions. This is due to different considerations. First, the selected mesh resolution has been chosen as small as possible in order to realize an extensive database. Second, the droplet is not in an infinite domain since the domain is considered as tri-periodic as in the database. Finally, the theory is applied for small perturbations (i.e., a _m,l ≪ R ₀), assumption that is not verified in the current simulations. Nonetheless, these values will serve as a reference for the rest of the article.

We investigate the oscillations by computing the coefficients a _m,l for each droplet of the database in order to correlate their oscillations with deformation. Since the frequency of each coefficient a _m,l for the same m is homologous, we make use of the power spectrum coefficient a m 2 (Eq. 14), which corresponds to the surface energy contained in each mode m Perrard et al. [14]. The evolution of a m = a m 2 as a function of the dimensionless time tf ₂ for a single droplet is shown in Figure 4. The choice of plotting a _m instead of a m 2 is carried out to highlight the coupling between modes. The amplitude of the coefficient a ₁ is close to zero for most of the duration of the droplet’s run due to the choice of the center of the reference frame, except when we observe high-amplitude deformations. In these extreme cases, the droplet becomes highly elongated with a complex geometry. We can observe a dominance of mode m = 2 in the surface energy deformation, in agreement with the previous literature in bubble dynamics [3,14,50], which have reported this phenomenon. To a lesser extent, the mode m = 4 also contributes to the global deformation and seems to be coupled with the mode m = 2 deformations. Considering the first peak at tf ₂ = 5, the total amount of energy stored at mode m = 2 is equivalent to 10% of the surface energy of a perfect sphere, while the same quantity for mode m = 4 does not reach 0.5%. The coefficients a _2,l are also shown in Figure 4. The coefficients a _m,l are directly dependent on the reference frame. If a different initial reference frame had been chosen, different values for the coefficients a _m,l would have been observed. However, a clear oscillation pattern is observed for each mode, and it is worth noting that the amplitudes of all the exposed coefficients are akin. The same behavior can be observed in the other modes, but with smaller amplitudes (not shown in the figure since larger degrees have too many orders to be plotted).

Due to the difficulty to extract a pure frequency from the amplitude coefficients a _m,l , the angular frequency of each couple (m, l) was obtained by computing the average oscillation period defined as the peak-to-peak interval time. The statistics were computed for tf ₂ > 6, where the statistically converged regime is reached (see section 3.3). The angular frequencies obtained are given in Table 1. The frequency observed for each mode is of the same order of magnitude as one for pure-deformed oscillations. This result broadens the conclusions of Risso and Fabre [3], which show for bubble in turbulence experiments that the power spectrum of the projected area is dominated by the natural frequencies of modes 2 and 4. These two frequencies appear here for the coefficient corresponding to modes 2 and 4, respectively. This proves that, as expected, the oscillatory frequency is related to the shape deformation given by the same mode. In addition, our results show that all the modes are concerned by these natural oscillations, not only even modes.

The shift between the computed frequency and the natural one is not clear: for mode 2, we observe an increase in the frequency, while for mode 6, we observe a decrease. The origin of these differences is difficult to identify here, but it could be due to the coupling of the interface with the turbulence and the coupling between modes. The seminal work of [51] for bubbles can provide the first answer of this problem. Nevertheless, the extension of this theoretical work to the liquid–liquid case is out of the scope of this paper. Another way to understand the differences consists in considering that mode 2 has a characteristic length scale close to the radius, and in our case, close to the inertial length scale, while the mode 6 characteristic length scale is smaller and close to the Kolmogorov length scale. Thus, the turbulence level seen by each mode is not the same. The turbulence impact on the frequency at each scale should be explored in the future by studying different turbulent Reynolds numbers and different droplet sizes.

We compute the ensemble average deformation by analyzing the entirety of the droplet database. For each droplet, we perform the spherical harmonic decomposition analysis described previously to obtain the spectrum values. In Figure 5, we show the temporal evolution of the spectrum of some realizations as a function of the dimensionless time tf ₂. The ensemble averages defined by α m ( t ) = a m 2 N are described by the black dashed lines. The definition here differs from the definition of Perrard et al. [14] who computes the ensemble average as a m N . Our choice is motivated by the fact that the saturation level should be linked to the total surface. It has been observed that both expressions provide a similar behavior except for mode 6, where the strong short deformations increase the averaged value with our definition. The computed values of a m N are not presented in this article.

If we take a closer look at the earlier time after the drop is left to evolve in the background flow, we observe an exponential growth for each mode. This growth differs from what is observed by Perrard et al. [14] where they noted a linear growth during the initial time of the evolution of the bubble in a similar experiment. The difference can be attributed to a different initialization routine. In their case, the bubble is initialized with a turbulent flow, even at the interface, while in our experiment, the droplet velocity is zero and a no-slip condition is ensured ([44]—forthcoming).

Then, a saturation level is reached after a given time. The time to reach this saturation time reduces with the mode. For all the modes, we observe a first peak of the ensemble average that seems to be well repeatable and certainly related to the initial condition. The time to reach this peak is smaller for larger modes. After this peak, all the modes converge toward a saturation level.

The saturation level has been studied for modes 2–4 by Perrard et al. [14]. In the present paper, we compute this value for modes 2–6. The saturation level here is the ensemble and time average over the 220 droplets and after tf ₂ > 6. The obtained average, defined by α m , s a t = a m 2 N , t f 2 > 6 , is {0.213, 0.0531, 0.0581, 0.0260, 0.0208}R ₀ for m = 2–6. As expected and observed by Perrard et al. [14], the saturation level decreases as the mode increases. That is due to two main phenomena. First, the natural damp rate increases with the mode, and thus, the deformations related to large modes reduce faster. Second, as the mode increases, the characteristic length decreases, and thus, the turbulent energy decreases. The relationship between these two phenomena is not clear, and the scaling of this saturation level has yet to be established. In particular, the difference between odd and even modes should be established. Future work, with different Weber and turbulent Reynolds numbers, will be really helpful in this direction.

After the analysis of regular oscillations, we now consider the study of large deformations.

We have seen that droplet oscillations in a turbulent flow result in a series of non-axisymmetric oscillations, which are difficult to characterize with the linear stability theory. Large deformation events can be relaxed or can lead to breakup, as suggested by the model of [4,11] which provides a maximum deformation threshold. In the following, we analyze the oscillations following a high-amplitude deformation. In order to bring out axisymmetry as much as possible, we rotate the reference frame, in which the spherical harmonics were computed to the one that maximizes the amplitude of coefficient a _2,0. A minimization Broyden–Fletcher–Goldfarb–Shanno algorithm is used to find the two Euler angles that maximize a _2,0. The corresponding reference frame is often close to the inertial tensor of the droplet (see [52]) that is here used as the initial condition for the minimization algorithm.

As seen in Figure 4, the studied droplets reach high-amplitude deformation several times during their lifetime. A total of 694 of these instances were chosen under the criteria that they have a coefficient a ₂(t _peak ) > 0.3R ₀ and do not experience breakup in the 4tf ₂ after this deformation. The threshold criteria were chosen to be 1.5 α _2,sat and is close to the breakup limit proposed by Lalanne et al. [4] and confirmed in the current database ([44]—forthcoming). In other words, these large deformations are often followed by breakup. In Figure 6, the evolution in time of the deformation coefficient a _m,0 following a high-amplitude deformation for individual realizations is plotted, along with the ensemble average of the chosen droplets.

First, the evolution of coefficient a _2,0 provides an idea of the total deformation evolution. A rapid increase of the a _m,0 coefficients for m = [2, 4] is followed by a contraction back to the previous mean level. In the particular case of mode 2 oscillations, when aligned to the maximum deformation, the fluctuation between growth and decay shows a timescale related to f ₂. Similarly to the observations of [50], this decay can be related directly to the natural damp rate of the droplet.

The coupling between even modes can be also established. First, it is clearly seen that these high-deformation droplets also have an a _4,0 component. In other words, the orientation of modes 2 and 4 deformations is aligned. The increase of coefficient a _6,0 is also to notice and seems to recover the energy relaxed by modes 2 and 4. For modes 3 and 5, there is no clear evolution.

In order to confirm that the physics after a high-amplitude deformation is related to the capillary forces and not to the background flow, we study a given realization and set the velocity to zero, similarly to the work carried out by Håkansson et al. [53] in the research for a breakup criterion. The selected droplet is present in Figure 6 (blue line), and its shape at t − t _peak is shown in Figure 7A. The reference frame is set as before using the maximization of a _2,0. In Figure 7B, the power spectrum evolution in time is given comparing the turbulent realization and the initial quiescent flow case. The initial decaying is clearly similar in both cases, then, for larger times, the turbulent case deviates from the perfect oscillator.

In Figures 7C, D, the decomposition of modes 2 and 4 is given, respectively. The other modes are not given since their amplitude is negligible here. The same conclusion can be drawn: the damping and first oscillations after the large deformation are driven by the capillary.

The present droplet has a complex deformation and is not composed of a pure deformation, as shown in section 3.1. Thus, the coupling between the modes can explain the differences observed between Figures 7C, D and Figures 3A–C. The droplet, which is shown to undergo a high-amplitude mode 2 deformation, spends a slightly larger part of each period in the prolate state (59.15% of the total time). This phenomenon is in accordance with the linear stability analysis theory [10,47]. The frequency and damp factor are given in 1 for mode a _2,0 and a _4,0. Similarly to the individual mode case (3.1), ω _m and β _m were computed by using a Levenberg–Marquardt fitting algorithm. Due to the length of the first oscillation observed for mode 4 (Figure 7D), the curve fit was computed from the a _4,0 amplitude after its first period. The main difference when comparing the pure mode deformation with the complex shape here is on the mode 4 damping rate. This difference can be understood by looking at Figure 7D. It is clearly seen that, in contrast to mode 2, mode 4 deformation is not perfectly aligned with the reference frame (a _4,2 and a _4,−1 are not negligible). Thus, the coupling between orders l plays an important role, and the damping cannot be evaluated only from the time evolution of coefficient a _4,0.