ORIGINAL RESEARCH article
Sec. Interdisciplinary Physics
Volume 11 - 2023 | https://doi.org/10.3389/fphy.2023.1173521
Droplet oscillations in a turbulent flow
- CORIA-UMR 6614 - Normandie Université, CNRS-Université and INSA de Rouen, Saint Etienne du Rouvray, France
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.
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 . 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 , 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 , 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 . Direct numerical simulation (DNS) has been extensively used to study interface–droplet/bubble turbulence interactions from the seminal works of Qian et al.  or Perlekar et al.  to the most recent works of Perrard et al.  or Vela-Martín and Avila . 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) , and Lattice–Boltzmann [13,20] approaches. A review of these methods can be found in Tryggvason et al. .
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 , and the study of the signature of the modes of deformation in order to infer bubble dynamics in a turbulent flow was presented . The development of a similar framework to describe bubble deformation  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 . 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 R0, 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:
where am,l(t) corresponds to the amplitude of the deformation expressed as a spherical harmonic coefficient and
The equation of motion of the free surface can be obtained from the linearization of the Navier–Stokes equations . For the particular case, in which the surrounding fluid is of negligible dynamical effects, the evolution in time of the coefficients am,l are solutions of the independent second-order differential equation :
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:
For the studied case of liquid–liquid droplets, the solution given by Prosperetti  is considered. In the limit of weak viscous effects, this equation can then be simplified by estimations obtained from an asymptotic development , which draws the following solution for the angular frequency:
Then, the expression for the damping rate is defined as follows:
The zeroth coefficient a0,0 evolves in time in order to ensure mass conservation . 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.
For the sake of clarity, in the current paper, we use the frequency f2 = 2πω2 computed from Eq. 7 to non-dimensionalize the time, considering that mode 2 is expected to dominate the deformation.
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 . It solves on a staggered Cartesian mesh the one-fluid formulation of the incompressible Navier–Stokes equation.
where ρ s the density, p is the pressure field, μ is the dynamic viscosity, D is the strain rate tensor,
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 . 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  and has been adapted to two-phase flows with particles  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,
The impact of the forcing scheme is currently under study  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
2.4 Spherical harmonic decomposition
To study the shape evolution of a drop immersed in a turbulent flow, we introduce the spherical harmonic functions
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 . The numerical error for computing the radius R(θ, φ) is less than 0.1%. Similarly, the numerical error for computing the coefficients am,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 .
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 R0.
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 am,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 (a1,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 . 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
The power spectrum
Since the coefficients am,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.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 am,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. Oscillations of a droplet initialized with a spherical harmonic mode coefficient m = 2–6. (A) a2,0=0.36R0, (B) a3,0=0.035R0, (C) a4,0=0.1R0, (D) a5,0=0.04R0, and (E) a6,0=0.016R0. The coefficients am,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 . This can be seen comparing the time spent on positive (i.e., prolate) and negative (i.e., oblate) am,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 ; Basaran ; Meradji et al. ; Zrnić et al. . 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  and Alonso . 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., am,l ≪ R0), assumption that is not verified in the current simulations. Nonetheless, these values will serve as a reference for the rest of the article.
3.2 Oscillations of droplets immersed in a turbulent flow
We investigate the oscillations by computing the coefficients am,l for each droplet of the database in order to correlate their oscillations with deformation. Since the frequency of each coefficient am,l for the same m is homologous, we make use of the power spectrum coefficient
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 tf2 of one droplet from the database. In the bottom, coefficients a2,l as a function of the dimensionless time tf2 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 am,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 tf2 > 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 , 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  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 tf2. The ensemble averages defined by
FIGURE 5. (A) Evolution of the deformation for each individual spherical harmonic mode coefficient am as a function of the dimensionless time tf2 (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.  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 (—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. . 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 tf2 > 6. The obtained average, defined by
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 a2,0. A minimization Broyden–Fletcher–Goldfarb–Shanno algorithm is used to find the two Euler angles that maximize a2,0. The corresponding reference frame is often close to the inertial tensor of the droplet (see ) 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 a2(tpeak) > 0.3R0 and do not experience breakup in the 4tf2 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.  and confirmed in the current database (—forthcoming). In other words, these large deformations are often followed by breakup. In Figure 6, the evolution in time of the deformation coefficient am,0 following a high-amplitude deformation for individual realizations is plotted, along with the ensemble average of the chosen droplets.
FIGURE 6. Evolution of the spherical harmonic coefficients am,0 following a high-amplitude deformation a2>0.3 (individual realizations from the database in solid lines) with its respective ensemble average
First, the evolution of coefficient a2,0 provides an idea of the total deformation evolution. A rapid increase of the am,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 f2. Similarly to the observations of , 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 a4,0 component. In other words, the orientation of modes 2 and 4 deformations is aligned. The increase of coefficient a6,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.  in the research for a breakup criterion. The selected droplet is present in Figure 6 (blue line), and its shape at t − tpeak is shown in Figure 7A. The reference frame is set as before using the maximization of a2,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. (A) Snapshot of the droplet in its maximum deformation time. (B) t − tpeak spectrum coefficient am as a function of the dimensionless time tf2 of a realization transferred into a quiescent flow. (C) Evolution of the spherical harmonic coefficients a2,l after rotating the reference frame to maximize the coefficient a2,0. (D) Evolution of the spherical harmonic coefficient a4,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 a2,0 and a4,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 a4,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 (a4,2 and a4,−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 a4,0.
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 a2,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.
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.
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.
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.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2023.1173521/full#supplementary-material
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
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
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
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
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).
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
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
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
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
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
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
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
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
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
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
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.
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
Keywords: drops and bubbles, turbulence, two-phase flows, linear theory, CFD
Citation: Roa I, Renoult M-C, Dumouchel C and Brändle de Motta JC (2023) Droplet oscillations in a turbulent flow. Front. Phys. 11:1173521. doi: 10.3389/fphy.2023.1173521
Received: 24 February 2023; Accepted: 24 April 2023;
Published: 30 May 2023.
Edited by:Carole Planchette, Graz University of Technology, Austria
Reviewed by:Marc Pradas, The Open University, United Kingdom
Michele La Rocca, Roma Tre University, Italy
Copyright © 2023 Roa, Renoult, Dumouchel and Brändle de Motta. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jorge César Brändle de Motta, email@example.com