Droplet oscillations in a turbulent flow

The oscillations of an initially unperturbed spherical droplet immersed in a homogeneous and isotropic turbulent background flow are investigated through spherical harmonic decomposition. As suggested in the literature, the shape oscillations under turbulent conditions are related to the frequency of droplets oscillating in a fluid without background flow. A series of direct numerical simulations (DNS) of droplets with single deformation modes in a fluid at rest are first performed. The frequency and damping rate are compared with weakly viscous linear theory. Then, a database of 220 droplets deformed under turbulent conditions for a single Weber and Reynolds number is generated with an identical numerical set-up. Each spherical harmonic coefficient shows an oscillatory motion with comparable frequency to the single deformation mode simulations. The power spectrum of the coefficients provides the amount of surface of each mode. After a transient regime, the surface area reaches a stationary saturation level. The saturation level of each mode is linked to the turbulence and the energy stored at the interface. Droplets after a high deformation are studied with and without background flow. As expected, the physics of relaxation is driven by capillary forces.

1 Introduction

Droplet and bubble oscillations have been a subject of research for over a century due to their relevance in transport processes across the interface. As the droplet departs from its spherical equilibrium state, its surface increases, inducing gradients of velocity, temperature, and species mass concentration, influencing the transport of momentum, heat, and mass [1]. Droplet deformation and breakup are important phenomena in the atomization process due to their impact in the final droplet size distribution. Many authors attribute the breakup of droplets due to turbulence to their free oscillations, in particular for low turbulence levels as a result of a resonant mechanism [2–4]. First studies of droplet oscillations are recorded from Lord Rayleigh [5], where he investigated the vibrations of a liquid mass about a spherical shape finding the well-known solutions of the fundamental modes of motion and calculated the corresponding frequencies. This work was later generalized by Lamb [6], adding the influence of the viscosity of the liquid sphere. Non-linear oscillations of inviscid drops and bubbles have been numerically studied for modes m = 2, 3, 4 [7,8], showing that finite viscosity has a strong effect on mode-coupling phenomena and demonstrating that a drop undergoing mode 2 oscillations spends a longer part of each period in a prolate form rather than in an oblate form. This conclusion was later reinforced through a weak stability analysis and numerical simulations of both the inviscid and viscous cases [9,10], along with quantifying the effects of viscosity in the damping and its role in enhancing the coupling of different oscillation modes. A common idea in various models in the literature is the existence of a critical Weber number (the ratio between turbulent and surface tension forces) for breakup. The problem with this criterion is that it is not universal as it is not suitable for cases where the residence time of the drop is large compared to the period of its shape oscillations [11]. Direct numerical simulation (DNS) has been extensively used to study interface–droplet/bubble turbulence interactions from the seminal works of Qian et al. [12] or Perlekar et al. [13] to the most recent works of Perrard et al. [14] or Vela-Martín and Avila [15]. Different numerical approaches have been used such as volume-of-fluid (VOF) [14,16], level set [17,18], coupled level-set method and volume-of-fluid (CLSVOF) [19], and Lattice–Boltzmann [13,20] approaches. A review of these methods can be found in Tryggvason et al. [21].

The current work is motivated by the previous literature, in which a spherical harmonic decomposition has been used as a mean to describe the deformed surface. The study of a liquid globe oscillating in a resting fluid [22,23] introduced the study of mode oscillations. Later on, the damping isolating the spherical harmonic coefficient for mode 2 was investigated [24], and the study of the signature of the modes of deformation in order to infer bubble dynamics in a turbulent flow was presented [3]. The development of a similar framework to describe bubble deformation [14] has also motivated this work. We wish to generalize the classical Lamb theory in order to account for more complex geometries that better represent realistic droplets with the introduction of spherical harmonics, being able to trace the oscillations not only in the azimuthal but also in the polar direction, as most of the literature assumes an axisymmetric behavior of the oscillations of droplets, which, in our case of study, is simply not realistic.

The utility of the present framework comes from its ability to observe the effect of the different length scales on the liquid bulk by isolating each mode and allowing us to describe the droplet’s oscillation frequency as an ensemble of the different mode frequencies present since we are able to monitor the oscillations in all directions.

This framework is used in a series of independent droplets immersed in a turbulent background flow generated by the DNS of the Navier–Stokes equations. The in-house code Archer is chosen for its capabilities to reproduce Lamb theory oscillations [25,26] and the level-set implicit description that allows an accurate computation of the spherical harmonic coefficients [27]. In the following work, we consider carrier and disperse phases with similar density and viscosity as a particular case between droplets and bubbles. We consider neither mass exchange nor mixing between both fluids. In addition, DNS of the individual mode of deformations in a fluid initially at rest, initialized with amplitudes comparable to turbulent mode oscillations, has been performed in order to compare their corresponding frequencies with those obtained from the study under turbulence effects. A good agreement between the obtained results and both experimental and numerical literature is observed, particularly exhibiting a mode 2 deformation dominance throughout the droplet lifetime and a strong coupling between the even modes (i.e., 2, 4, and 6). Through a study of the decay following high-amplitude deformations in both turbulent and quiescent flows, it is confirmed that the oscillator is driven by capillary forces.

This paper is organized as follows: first, in section 1.1, a brief review of the existing theory of droplet oscillations is presented. Section 2 focuses on the different aspects of the methodology: the numerical solver, the generation of the turbulent flow, and the spherical harmonic framework. In section 3, we focus on the analysis of the numerical results. First, we present different reference cases where initial droplets are deformed with a specific shape in a fluid at rest. Second, we analyze the evolution of the deformation of a droplet immersed in a turbulent background flow, showing the dominance of mode 2 in the deformation dynamics in agreement with the experimental and numerical literature. Then, attention is set on the initial exponential growth of the droplet deformation and this deformation saturation level. A specific subsection is carried out on the relaxation time after high deformation. Finally, in section 4, we close our analysis with discussions and conclusions.

1.1 Droplet oscillations in a resting fluid

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\left(\theta ,\phi ,t\right)=R_0{\displaystyle \sum _{m=0}^{\infty }\sum _{l=-m}^m}{a}_{m,l}\left(t\right)Y_m^l\left(\theta ,\phi \right),(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\left(\theta ,\phi \right)=\sqrt{\left(2m+1\right)\frac{\left(m-l\right)!}{\left(m+l\right)!}}{P}_m^l\left(cos\theta \right)e^{il\phi },(2)$

with $P_m^l(\cos \theta )$ corresponding to the Legendre polynomials, $m\in \mathbb{N}$ and l ∈ [−m, m]. Then, the real form of the spherical harmonics is given by the following equation:

$Y_m^l\left(\theta ,\phi \right)=\left\{\begin{array}{cc}\frac{i}{\sqrt{2}}\left(H_m^l-{\left(-1\right)}^m{H}_m^{-l}\right)\hfill & \text{if\hspace{0.17em}}m<0\hfill \\ H_m^0\hfill & \text{if\hspace{0.17em}}m=0\hfill \\ \frac{i}{\sqrt{2}}\left(H_m^{-l}+{\left(-1\right)}^m{H}_m^l\right)\hfill & \text{if\hspace{0.17em}}m>0.\hfill \end{array}\right.(3)$

FIGURE 1

FIGURE 1. Drop deformation: sketch and definitions.

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]:

${\ddot{a}}_{m,l}+2{\beta }_{m0}{\dot{a}}_{m,l}+{\omega }_{m0}^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:

${\omega }_{m0}=\sqrt{\frac{m\left(m+1\right)\left(m-1\right)\left(m+2\right)\sigma }{R_0^3\left(m{\rho }_g+\left(m+1\right){\rho }_l\right)}},(5)$

${\beta }_{m0}=\left(2m+1\right)\left(m-1\right)\frac{{\mu }_l}{{\rho }_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:

${\omega }_m={\omega }_{m0}\left[1-\underset{A}{\underbrace{\frac{{\left(2m+1\right)}^2\sqrt{\widehat{\rho }\widehat{\mu }}}{2\sqrt{2{\omega }_{m0}\frac{{\rho }_l}{{\mu }_l}}{R}_0\left[m\widehat{\rho }+\left(m+1\right)\right]\left(1+\sqrt{\widehat{\rho }\widehat{\mu }}\right)}}}\right].(7)$

Then, the expression for the damping rate is defined as follows:

${\beta }_m={\omega }_{m0}\left(A-2A^2\right)+\frac{\left(2+m\right)\left[2\left(m^2-1\right)+\left(m+2\right)\widehat{\mu }-\left(m-1\right)\widehat{\rho }\widehat{\mu }+2m\left(m+2\right)\widehat{\rho }{\widehat{\mu }}^2\right]}{2\left(m\widehat{\rho }+m+1\right){\left[1+\sqrt{\widehat{\rho }\widehat{\mu }}\right]}^2}\frac{{\mu }_l}{R_0^2{\rho }_l},(8)$

where $\widehat{\rho }={\rho }_g/{\rho }_l$ is the density ratio and $\widehat{\mu }={\mu }_g/{\mu }_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.

2 Methodology

2.1 Archer solver

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.

$\nabla \cdot \vec{u}=0,(9)$

${\partial }_t\rho \vec{u}+\vec{\nabla }\cdot \left(\rho \vec{u}\otimes \vec{u}\right)=-\vec{\nabla }p+\vec{\nabla }\cdot \left(2\mu \mathbf{D}\right)+\vec{f}+\sigma \kappa {\delta }_s\vec{n},(10)$

where ρ s the density, p is the pressure field, μ is the dynamic viscosity, D is the strain rate tensor, $\vec{f}$ is a source term, σ is the surface tension, $\vec{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.

2.2 Turbulent flow generation

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, $\vec{f}=A{\vec{u}}^{\prime }$. Here, ${\vec{u}}^{\prime }=\vec{u}-⟨\vec{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}_{\lambda }=\frac{{\rho }_{l}u\prime \lambda }{{\mu }_l}\simeq 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.

2.3 Droplets in the turbulent flow database

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 $We=\frac{{\rho }_l{u}^{\prime 2}2R_0}{\sigma }$, 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 $\frac{R_0}{L}=0.233$ and an Ohnesorge number $(Oh=\frac{{\mu }_l}{\sqrt{\sigma {\rho }_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).

2.4 Spherical harmonic decomposition

To study the shape evolution of a drop immersed in a turbulent flow, we introduce the spherical harmonic functions $Y_m^l(\theta ,\phi )$ (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}={\int }_{\mathrm{\Omega }}R\left(\theta ,\phi \right)Y_m^l\left(\theta ,\phi \right)d\mathrm{\Omega },(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].

FIGURE 2

FIGURE 2. Illustration of the interface mapping of a droplet deformed by a turbulent background flow for a spherical harmonic decomposition with M =6. (A) Points are distributed around a revolution axis given by the reference frame. (B) Normalized distances between the center and the interface are stored in a grid. The colorbar corresponds in both figures to the values of the function R(θ, φ) normalized by the reference radius R₀.

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 ${\eta }_{\mathrm{\Omega }}^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:

$\frac{1}{4\pi }{\int }_{\mathrm{\Omega }}{R}^2\left(\theta ,\varphi \right)d\mathrm{\Omega }={\eta }_{\mathrm{\Omega }}^2,(12)$

where

${\eta }_{\mathrm{\Omega }}^2={\displaystyle \sum _{m=0}^{\infty }\sum _{l=-m}^m}{a}_{m,l}^2.(13)$

The power spectrum $a_m^2$ is computed by the following equation:

$a_m^2={\displaystyle \sum _{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.

3 Results

3.1 Single-mode deformation droplets

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.

FIGURE 3

FIGURE 3. Oscillations of a droplet initialized with a spherical harmonic mode coefficient m = 2–6. (A) a_2,0=0.36R₀, (B) a_3,0=0.035R₀, (C) a_4,0=0.1R₀, (D) a_5,0=0.04R₀, and (E) a_6,0=0.016R₀. The coefficients a_m,0 were selected in accordance to the highest amplitude of the coefficients found in turbulent flow deformation (presented in section 3.2) (Supplementary Videos S1–S5 and Supplementary Data S1).

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.

TABLE 1

TABLE 1. Summary of frequencies ω and the damping ratio β normalized by ω₂₀.

3.2 Oscillations of droplets immersed in a turbulent flow

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=\sqrt{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).

FIGURE 4

FIGURE 4. Top image shows the power spectrum of the spherical harmonic decomposition for modes m =1 to m =6 as a function of the dimensionless time tf₂ of one droplet from the database. In the bottom, coefficients a_2,l as a function of the dimensionless time tf₂ of one droplet from the database for l =−2 to l =2 are shown (Supplementary Video S6 and Supplementary Data S2).

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.

3.3 Initial growth and saturation regime

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 ${\alpha }_m(t)=\sqrt{{⟨a_m^2⟩}_{\mathcal{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⟩}_{\mathcal{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⟩}_{\mathcal{N}}$ are not presented in this article.

FIGURE 5

FIGURE 5. (A) Evolution of the deformation for each individual spherical harmonic mode coefficient a_m as a function of the dimensionless time tf₂ (individual realizations from the database in solid lines) with its respective ensemble average α_m (dashed line) for We =0.9. (B) Close up of the initial growth regime (Supplementary Videos S7–S12 and Supplementary Data S3).

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 ${\alpha }_{m,sat}=\sqrt{{⟨a_m^2⟩}_{\mathcal{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.

3.4 Large deformations and damping factor

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.

FIGURE 6

FIGURE 6. Evolution of the spherical harmonic coefficients a_m,0 following a high-amplitude deformation a₂>0.3 (individual realizations from the database in solid lines) with its respective ensemble average ${⟨a_{m,0}⟩}_{\mathcal{N}}$ (dashed line) (Supplementary Video S13 and Supplementary Data S4).

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.

FIGURE 7

FIGURE 7. (A) Snapshot of the droplet in its maximum deformation time. (B) t − t_peak spectrum coefficient a_m as a function of the dimensionless time tf₂ of a realization transferred into a quiescent flow. (C) Evolution of the spherical harmonic coefficients a_2,l after rotating the reference frame to maximize the coefficient a_2,0. (D) Evolution of the spherical harmonic coefficient a_4,l in the same reference frame (Supplementary Video S14 and Supplementary Data S5).

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.

4 Conclusion

A droplet database generated by DNS was used to characterize the oscillations of droplets. A framework for droplet deformation using a spherical harmonic decomposition of the droplet surface has been presented.

As a reference case, droplets with single deformation modes are first considered. The oscillatory motion is consistent with the theoretical framework. The angular frequency and damping rate of the oscillator show a slight discrepancy compared to the linear theory. As expected, the droplets spend more time in the prolate state than in the oblate state, and the coupling between even modes seems to be important.

Droplets immersed in a turbulent background flow are studied. The present work is based on a database that is still being extended for different purposes, in particular to investigate the droplet breakup and the turbulent/interface interactions. This database considers many realizations with similar dimensionless numbers.

It is shown that mode 2 is dominant in the deformation of a droplet undergoing turbulence-driven deformation. Oscillations of each spherical harmonic coefficient are observed, and its frequency is close to the theoretical frequency obtained from the linear theory. These observations extend the existing literature which focuses only on the turbulence coupling with mode 2.

The total deformation is analyzed using the power spectrum. The ensemble average of the power spectrum, which can be seen as the energy stored at the interface for each mode, reaches a statistically converged level. This saturation level is related to the droplets natural damping rate and the turbulent cascade. The scaling of this saturation level has not yet been established.

Finally, a study of high-amplitude deformations is performed. The spherical harmonic reference frame is aligned with a large deformation axis, maximizing the a_2,0 coefficient. For these deformations, the coupling between even modes is observed. The droplet relaxation is related to capillary forces, as shown by the simulation of the same droplet but without any initial velocity field. The damping rate is similar to the one obtained in the single-mode simulations.

Three main perspectives can be drawn and are being explored by different groups. First, it has been seen that the weakly viscous linear theory provides accurate information on how the droplets evolve. Two main questions can be answered by extending the existing theories: i. how modes are coupled in the case of complex shape droplets; and ii. how singular eddies excite each mode modifying the angular frequency oscillations. Second, new simulations with different Weber and Reynolds numbers should be performed in order to provide a more detailed picture of how the turbulence and interface interact at each length scale. From the authors’ point of view, it is the straight way to understand how the saturation level for each mode power spectrum scales with the turbulence. Finally, the extension of the present work to droplets in a gaseous environment and bubbles in a liquid, numerically and experimentally, will improve the understanding of turbulence–interface interactions.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Materials; further inquiries can be directed to the corresponding author.

Author contributions

All authors contributed to the conception and design of the study. IR wrote the first draft of the manuscript. IR and JB performed the simulations presented in this work. All authors contributed to the manuscript revision, and read and approved the submitted version. JB obtained the funding for the project. All authors contributed to the article and approved the submitted version.

Funding

This work is performed under the auspices of the French National Research Agency (Agence Nationale de la Recherche ANR) under project ANR-20-CE46-0002.

Acknowledgments

The authors acknowledge F. Risso and F. Thiesset for the discussions and ideas they shared with them. This work was granted access to the HPC resources of IDRIS, TGCC, and CINES under the allocation A0132B1010 made by GENCI (Grand Equipement National de Calcul Intensif) and the CRIANN (Centre Régional Informatique et d’Applications Numériques de Normandie) under the scientific project N. 2017 004.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2023.1173521/full#supplementary-material

References

1. Germes Martinez L, Duret B, Reveillon J, Demoulin FX. A new dns formalism dedicated to turbulent two-phase flows with phase change. Int J Multiphase Flow (2021) 143:103762. doi:10.1016/j.ijmultiphaseflow.2021.103762

CrossRef Full Text | Google Scholar

2. Sevik M, Park SH. The splitting of drops and bubbles by turbulent fluid flow. ASME. J Fluids Eng (1973) 95 (1):53–60. doi:10.1115/1.3446958

CrossRef Full Text | Google Scholar

3. Risso F, Fabre J. Oscillations and breakup of a bubble immersed in a turbulent field. J Fluid Mech (1998) 372:323–55. doi:10.1017/s0022112098002705

CrossRef Full Text | Google Scholar

4. Lalanne B, Masbernat O, Risso F. A model for drop and bubble breakup frequency based on turbulence spectra. AIChE J (2019) 65(1):347–59. doi:10.1002/aic.16374

CrossRef Full Text | Google Scholar

5. Rayleigh L. On the capillary phenomena of jets. Proc R Soc Lond (1879) 29(30):71–97.

Google Scholar

6. Lamb H. On the vibrations of an elastic sphere. Proc Lond Math Soc (1881) 1(1):189–212. doi:10.1112/plms/s1-13.1.189

CrossRef Full Text | Google Scholar

7. Tsamopoulos JA, Brown RA. Nonlinear oscillations of inviscid drops and bubbles. J Fluid Mech (1983) 127:519–37. doi:10.1017/s0022112083002864

CrossRef Full Text | Google Scholar

8. Basaran OA. Nonlinear oscillations of viscous liquid drops. J Fluid Mech (1992) 241:169–98. doi:10.1017/s002211209200199x

CrossRef Full Text | Google Scholar

9. Zrnić D, Brenn G. Weakly nonlinear shape oscillations of inviscid drops. J Fluid Mech (2021) 923:A9. doi:10.1017/jfm.2021.568

CrossRef Full Text | Google Scholar

10. Zrnić D, Berglez P, Brenn G. Weakly nonlinear shape oscillations of a Newtonian drop. Phys Fluids (2022) 34(4):043103. doi:10.1063/5.0085070

CrossRef Full Text | Google Scholar

11. Risso F. The mechanisms of deformation and breakup of drops and bubbles. Multiphase Sci Technol (2000) 12(1):50. doi:10.1615/multscientechn.v12.i1.10

CrossRef Full Text | Google Scholar

12. Qian T, Wang X-P, Sheng P. A variational approach to moving contact line hydrodynamics. J Fluid Mech (2006) 564:333–60. doi:10.1017/s0022112006001935

CrossRef Full Text | Google Scholar

13. Perlekar P, Biferale L, Sbragaglia M, Srivastava S, Toschi F. Droplet size distribution in homogeneous isotropic turbulence. Phys Fluids (2012) 24(6):065101. doi:10.1063/1.4719144

CrossRef Full Text | Google Scholar

14. Perrard S, Rivière A, Mostert W, Deike L. Bubble deformation by a turbulent flow. J Fluid Mech (2021) 920:A15. doi:10.1017/jfm.2021.379

CrossRef Full Text | Google Scholar

15. Vela-Martín A, Avila M. Memoryless drop breakup in turbulence. Sci Adv (2022) 8(50):eabp9561. doi:10.1126/sciadv.abp9561

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kékesi T, Amberg G, Wittberg LP. Drop deformation and breakup. Int J Multiphase Flow (2014) 66:1–10. doi:10.1016/j.ijmultiphaseflow.2014.06.006

CrossRef Full Text | Google Scholar

17. Trontin P, Vincent S, Estivalezes J, Caltagirone J. Direct numerical simulation of a freely decaying turbulent interfacial flow. Int J Multiphase Flow (2010) 36(11-12):891–907. doi:10.1016/j.ijmultiphaseflow.2010.08.003

CrossRef Full Text | Google Scholar

18. Shao C, Luo K, Yang Y, Fan J. Direct numerical simulation of droplet breakup in homogeneous isotropic turbulence: The effect of the weber number. Int J Multiphase Flow (2018) 107:263–74. doi:10.1016/j.ijmultiphaseflow.2018.06.009

CrossRef Full Text | Google Scholar

19. Duret B, Canu R, Reveillon J, Demoulin F. A pressure based method for vaporizing compressible two-phase flows with interface capturing approach. Int J Multiphase Flow (2018) 108:42–50. doi:10.1016/j.ijmultiphaseflow.2018.06.022

CrossRef Full Text | Google Scholar

20. Montessori A, Rocca ML, Prestininzi P, Tiribocchi A, Succi S. Deformation and breakup dynamics of droplets within a tapered channel. Phys Fluids (2021) 33(8):082008. doi:10.1063/5.0057501

CrossRef Full Text | Google Scholar

21. Tryggvason G, Scardovelli R, Zaleski S. Direct numerical simulations of gas–liquid multiphase flows. Cambridge: Cambridge University Press (2011).

Google Scholar

22. Chandrasekhar S. The oscillations of a viscous liquid globe. Proc Lond Math Soc (1959) 3(1):141–9. doi:10.1112/plms/s3-9.1.141

CrossRef Full Text | Google Scholar

23. Miller C, Scriven L. The oscillations of a fluid droplet immersed in another fluid. J Fluid Mech (1968) 32(3):417–35. doi:10.1017/s0022112068000832

CrossRef Full Text | Google Scholar

24. Prosperetti A. Free oscillations of drops and bubbles: The initial-value problem. J Fluid Mech (1980) 100(2):333–47. doi:10.1017/s0022112080001188

CrossRef Full Text | Google Scholar

25. Aniszewski W, Ménard T, Marek M. Volume of fluid (vof) type advection methods in two-phase flow: A comparative study. Comput Fluids (2014) 97:52–73. doi:10.1016/j.compfluid.2014.03.027

CrossRef Full Text | Google Scholar

26. Cordesse P, Di Battista R, Chevalier Q, Matuszewski L, Ménard T, Kokh S, et al. A diffuse interface approach for disperse two-phase flows involving dual-scale kinematics of droplet deformation based on geometrical variables. In: P Helluy, J-M Hérard, and N Seguin, editors. Esaim: Proceedings and surveys. Cambridge: Cambridge University Press (2020).

CrossRef Full Text | Google Scholar

27. Roa I, Brändle de Motta JC, Poux A. Spherical harmonics decomposition based on Cartesian level-set field. France: Université de Rouen (2023). Technical report. hal-04004128.

Google Scholar

28. Prosperetti A. Viscous effects on perturbed spherical flows. Q Appl Math (1977) 34(4):339–52. doi:10.1090/qam/99652

CrossRef Full Text | Google Scholar

29. Lu H-L, Apfel RE. Shape oscillations of drops in the presence of surfactants. J Fluid Mech (1991) 222:351–68. doi:10.1017/s0022112091001131

CrossRef Full Text | Google Scholar

30. Kornek U, Müller F, Harth K, Hahn A, Ganesan S, Tobiska L, et al. Oscillations of soap bubbles. New J Phys (2010) 12(7):073031. doi:10.1088/1367-2630/12/7/073031

CrossRef Full Text | Google Scholar

31. Brenn G. Analytical solutions for transport processes, mathematical engineering. Salmon: Springer Berlin Heidelberg (2017).

Google Scholar

32. Ménard T, Tanguy S, Berlemont A. Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet. Int J Multiphase Flow (2007) 33(5):510–24. doi:10.1016/j.ijmultiphaseflow.2006.11.001

CrossRef Full Text | Google Scholar

33. Rudman M. A volume-tracking method for incompressible multifluid flows with large density variations. Int J Numer Methods Fluids (1998) 28(2):357–78. doi:10.1002/(sici)1097-0363(19980815)28:2<357::aid-fld750>3.0.co;2-d

CrossRef Full Text | Google Scholar

34. Vaudor G, Ménard T, Aniszewski W, Doring M, Berlemont A. A consistent mass and momentum flux computation method for two phase flows. application to atomization process. Comput Fluids (2017) 152:204–16. doi:10.1016/j.compfluid.2017.04.023

CrossRef Full Text | Google Scholar

35. Sussman M, Smith K, Hussaini M, Ohta M, Zhi-Wei R. A sharp interface method for incompressible two-phase flows. J Comput Phys (2007) 221(2):469–505. doi:10.1016/j.jcp.2006.06.020

CrossRef Full Text | Google Scholar

36. Zhang J. Acceleration of five-point red-black gauss-seidel in multigrid for Poisson equation. Appl Math Comput (1996) 80(1):73–93. doi:10.1016/0096-3003(95)00276-6

CrossRef Full Text | Google Scholar

37. Fedkiw RP, Aslam T, Merriman B, Osher S. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J Comput Phys (1999) 152(2):457–92. doi:10.1006/jcph.1999.6236

CrossRef Full Text | Google Scholar

38. Rosales C, Meneveau C. Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties. Phys Fluids (2005) 17(9):095106. doi:10.1063/1.2047568

CrossRef Full Text | Google Scholar

39. Brändle de Motta JC, Estivalezes J-L, Climent E, Vincent S. Local dissipation properties and collision dynamics in a sustained homogeneous turbulent suspension composed of finite size particles. Int J Multiphase Flow (2016) 85:369–79. doi:10.1016/j.ijmultiphaseflow.2016.07.003

CrossRef Full Text | Google Scholar

40. Duret B, Luret G, Reveillon J, Ménard T, Berlemont A, Demoulin F-X. Dns analysis of turbulent mixing in two-phase flows. Int J Multiphase Flow (2012) 40:93–105. doi:10.1016/j.ijmultiphaseflow.2011.11.014

CrossRef Full Text | Google Scholar

41. Germes Martinez L, Duret B, Reveillon J, Demoulin FX. Vapor mixing in turbulent vaporizing flows. Int J Multiphase Flow (2023) 161:104388. doi:10.1016/j.ijmultiphaseflow.2023.104388

CrossRef Full Text | Google Scholar

42. Aniszewski W, Brändle de Motta JC. A wide-range parameter study for the turbulence-interface interactions. In: 14th European Fluid Mechanics Conference (2022).

Google Scholar

43. Chéron V, Brändle de Motta JC, Blaisot J-B, Ménard T. Analysis of the effect of the 2D projection on droplet shape parameters. Atomization and Sprays (2022) 32(8):59–98. doi:10.1615/atomizspr.2022040525

CrossRef Full Text | Google Scholar

44. Deberne C, Chéron V, Poux A, Brändle de Motta JC. Breakup prediction of oscillating droplets under turbulent flow (2023) Forthcoming.

Google Scholar

45. Wieczorek MA, Meschede M. Shtools: Tools for working with spherical harmonics. Geochem Geophys Geosystems (2018) 19(8):2574–92. doi:10.1029/2018gc007529

CrossRef Full Text | Google Scholar

46. Trinh E, Wang T. Large amplitude drop shape oscillations. In: Proc. Of the 2d intern. Colloq. On drops and bubbles (1982).

Google Scholar

47. Foote GB. A numerical method for studying liquid drop behavior: Simple oscillation. J Comput Phys (1973) 11(4):507–30. doi:10.1016/0021-9991(73)90135-6

CrossRef Full Text | Google Scholar

48. Meradji S, Lyubimova T, Lyubimov D, Roux B. Numerical simulation of a liquid drop freely oscillating. Cryst Res Technol J Exp Ind Crystallogr (2001) 36(7):729–44. doi:10.1002/1521-4079(200108)36:7<729::aid-crat729>3.0.co;2-3

CrossRef Full Text | Google Scholar

49. Alonso CT. The dynamics of colliding and oscillating drops. JPL Proc Intern Colloq Drops Bubbles (1974) 1:139–157.

Google Scholar

50. Ravelet F, Colin C, Risso F. On the dynamics and breakup of a bubble rising in a turbulent flow. Phys Fluids (2011) 23(10):103301. doi:10.1063/1.3648035

CrossRef Full Text | Google Scholar

51. Kang IS, Leal LG. Small-amplitude perturbations of shape for a nearly spherical bubble in an inviscid straining flow (steady shapes and oscillatory motion). J Fluid Mech (1988) 187:231–66. doi:10.1017/s0022112088000412

CrossRef Full Text | Google Scholar

52. Lalanne B. Simulation numérique directe de la déformation, des oscillations et de la rupture d’une bulle en ascension dans un écoulement instationnaire. France: Université de Toulouse (2012). PhD thesis.

Google Scholar

53. Håkansson A, Crialesi-Esposito M, Nilsson L, Brandt L. A criterion for when an emulsion drop undergoing turbulent deformation has reached a critically deformed state. Colloids Surf A: Physicochemical Eng Aspects (2022) 648:129213. doi:10.1016/j.colsurfa.2022.129213

CrossRef Full Text | Google Scholar