A diffuse interface approach for disperse two-phase flows involving dual-scale kinematics of droplet deformation based on geometrical variables

The purpose of this contribution is to derive a reduced-order two-phase flow model in- cluding interface subscale modeling through geometrical variables based on Stationary Action Principle (SAP) and Second Principle of Thermodynamics in the spirit of [6, 14]. The derivation is conducted in the disperse phase regime for the sake of clarity but the resulting paradigm can be used in a more general framework. One key issue is the definition of the proper potential and kinetic energies in the Lagrangian of the system based on geometrical variables (Interface area density, mean and Gauss curvatures...), which will drive the subscale kinematics and dissipation, and their coupling with large scales of the flow. While [14] relied on bubble pulsation, that is normal deformation of the interface with shape preservation related to pressure changes, we aim here at tackling inclusion deformation at constant volume, thus describing self-sustained oscillations. In order to identify the proper energies, we use Direct Numerical Simulations (DNS) of oscillating droplets using ARCHER code and recently devel- oped library, Mercur(v)e, for mean geometrical variable evaluation and analysis preserving topological invariants. This study is combined with historical analytical studies conducted in the small perturba- tion regime and shows that the proper potential energy is related to the surface difference compared to the spherical minimal surface. A geometrical quasi-invariant is also identified and a natural definition of subscale momentum is proposed. The set of Partial Differential Equations (PDEs) including the conservation equations as well as dissipation source terms are eventually derived leading to an original two-scale diffuse interface model involving geometrical variables.


Introduction
The simulation of two-phase flows in many industrial processes requires to account for physical phenomena that operate over a large range of different scales. For example, in the case of jet atomization in sub-critical conditions, the interface between the liquid and the gas experiences a dramatic change of topology from a separate-phase flow regime to a disperse flow regime. Although Direct Numerical Simulations (DNS) have provided impressive results in this field [3,4,5,6], they remain too costly in terms of computing resources to be applied in an industrial context or conduct parametrical studies. Therefore, the derivation of averaged models is still an important question and it has to cope with the predictive modeling of small scales of interfacial flows, which can not be resolved. Several authors have proposed various means to enrich these reduced-order models by introducing additional flow parameters that are reminiscent of small scale features that cannot be described by the bulk variables such as the density of interfacial area [7], or local curvature of the interface to account for small scale dynamics connected to capillary effects [8]. Works in this direction are also found in [9,10] where a unified model accounting for micro-inertia and micro-viscosity associated to bubble pulsation is proposed.
In this paper, we propose to focus on the coupling in a disperse flow between small and large scales through the kinematics description by revisiting the approach proposed in [2,10] in a different physical setting. More specifically, we propose to study means to supplement the bulk velocity that allows to represent large scale motions with additional parameters that describe small scale properties of the interface. Such approach relies on several hypotheses related to the evolution of the small scale features of the interface. The first assumption concerns the choice of kinetic and potential energies associated with the small scale motion. The second assumption deals with the connections between the small scale parameters. A set of constraints expressed by Partial Differential Equations (PDEs) has been proposed in [2] by considering that small scale perturbations behave in the same way as small amplitude motions along the normal to the interface. In this contribution, we propose to examine a complementary type of small scale motions by considering ellipsoidal deformed droplets oscillating around a spherical shape. An analytic study of this problem has been achieved in [11] thanks to a harmonic analysis based on the work of Taylor and was also reproduced in [12] about the well-known TAB model for secondary break-up. More recently, [8] considered a harmonic oscillator involving another set of variables, including mean curvature. Revisiting the theory allows to show that the various approaches are equivalent in the limit of small perturbations but we aim at deriving a model valid with large deformations. We thus need to propose the proper set of kinetic and potential energies related to the droplet deformations and the proper variables in order to describe such physics. We propose to rely on a DNS of a single droplet performed with the code ARCHER and a post-processing library Mercur(v)e, which is designed for mean geometrical variable evaluation and analysis preserving topological invariants. The results of this simulation suggest a specific choice for both the constraints to be applied on the small scale dynamics and the small scale kinetic and potential energies defined as function of the proper variable, which is the surface variation compared to the spherical minimal surface of the inclusions. The system of PDEs is eventually derived and involves both large scale and small scales dynamics as well as their coupling through convection and dissipation.
The paper is organized as follows: we first briefly recall the main results obtained by an harmonic analysis of the motion of an ellipsoidal droplet performed in [11] in the limit of an infinitely small perturbation and extend the derivations to obtain analytical expressions for several geometric variables involved in the dynamics beyond the kinetic energy and the potential energy. We then introduce the DNS results obtained with the code ARCHER and post-processed by the library Mercur(v)e in order to estimate the evolution of geometric variables and associated potential energies. We identify the proper variable for the potential energy to remain harmonic in the large deformation regime, thus identifying the best candidate to use the SAP. Finally, we propose a new set of PDEs that relates parameters accounting for small scale motion of the interface and the new form of potential and kinetic energies. The proposed form of the equations and identification of the various contribution in the Lagrangian naturally allow to separate the properties of the interface between large scales and small scales in the two-phase system. This is a specific case of a more general approach where a systematic framework is proposed in order to conduct such a study in the general case of separated and disperse phases [13]. Thanks to the Stationary Action Principle, we obtain a set of evolution equations for the medium without dissipation that enables us to assess the hyperbolic behaviour of the system. Eventually, we show that this system of equations can be complemented with a mathematical entropy evolution equation that is compatible with the Second Principle of Thermodynamics. This helps bringing up additional terms to the system that account for internal dissipation mechanisms.

Oscillating inclusion literature review
In [14], Sir G.I. Taylor studied secondary breakup through droplet deformation under a non-homogeneous pressure field and surface tension. He noticed an oscillating phenomena among the droplets, which he saw as an important source of dissipation for the system. This phenomenon leads to a higher pressure difference than expected before breakup of the droplets, as the energy from the pressure wave gets converted into mechanical oscillating energy and thus is rapidly dissipated. He concluded that in order to get a satisfactory model of secondary breakup, the frequency of the incoming pressure waves should be compared to the oscillating frequency of the droplets since some frequencies may trigger droplet oscillation rather than breakup.
Later on, in [11], S. Chandrasekhar considered an object with a deformation around its original spherical shape governed by a spherical harmonic function. He solved the velocity field inside a spherical droplet that is only subjected to its own weight, and then drew a formal analogy with a droplet experiencing surface tension. The exact resolution of the velocity field is conducted through solving Laplace's equation using the potential nature of gravity, and using Kelvin's formulas in the limit of small displacements. In this approach, the key variable leading to a harmonic oscillation is the amplitude of the spherical harmonic function, the frequency of oscillation being related to the characteristics of the spherical harmonic perturbation and droplet physical parameters. The derivations were extended for bubbles in [15] and the reader is referred to [16] for a more modern view on the literature on which we will come back later.

Analytical energy balance of an isolated oscillating inclusion
The analytic derivation of the flow field arising in a droplet of radius R, oscillating due to capillarity, the surface tension coefficient of which is denoted σ, may be found in Chandrasekhar [11,17] and we briefly recall here some of the main results.
In the limit of small initial perturbations and considering one deformed incompressible spherical inclusion, the governing equations can be linearized as leading to the harmonic pressure field ∆p = 0.
Under the assumption that the initial surface deformation of the droplet is described by spherical harmonic functions Y m l , the flow field may be explicitly calculated as a harmonic function of time for both inviscid and viscous cases and thus yields the evolution of kinetic K and potential U energies per unit of volume during the oscillation of the droplet.
We choose to focus on a small axisymmetric ellipsoidal perturbation with of non-dimensionalized magnitude η << 1 (η = p/q − 1, where p and q are the large and small axis of the ellipse -see Appendix A), which is applied to a spherical liquid droplet. Thus, the initial perturbation surface function of r/R may be approximated by the Y 0 2 spherical harmonic, as detailed in Appendix A. It leads to the following expression for the kinetic and potential energies per unit of volume where the reference interfacial density area Σ 0 is defined as the reference surface S 0 = 4πR 2 divided by the volume V of interest (which will be the volume of the simulation in the DNS configuration). It is then clear that the potential energy, in this limit, is a harmonic function, reaching zero every time a spherical shape is recovered, the frequency of which is clearly identified and coherent with various contributions in the literature. Let us underline that the system naturally admits an invariant quantity in the limit of small perturbations Σ + 3 H α = 0 (see Appendix A), where α = V 0 /V is a constant, V 0 = 4πR 3 /3, Σ is the surface area density and H the surface averaged mean curvature of the droplet. However, at this level, the key issue is the functional dependency of the potential energy since several quantities are harmonic oscillators in the limit of small perturbations.
Indeed, in the asymptotic case of small deformations, two expressions are eligible for the potential energy of the oscillating droplet. The first is directly related to the excess surface of the perturbed droplet and reads: the second one is, following Herrmann [8], related to the surface averaged mean curvature H of the droplet and reads where Σ and H expressions are given in Appendix A and where H 0 = −1/R. Furthermore, in [12], a harmonic oscillator model was derived providing a new modeling strategy called TAB (for Taylor Analogy Breakup). This model allows for a better handling of secondary breakup through modeling high Weber number spherical droplets as ellipsoids with axisymmetric forms, a framework closely related to our problem. The key parameter in this approach leading to harmonic oscillations is the displacement of the equator of the droplet from its equilibrium position noted x.
In the following, a DNS-based study will be conducted in order to discriminate what should be the proper variable underlying the definition of the potential energy, keeping in mind that several variables are harmonic functions in the limit of small displacement, but we aim at defining a subscale model relying on an oscillator that is valid beyond the regime of infinitely small perturbations.

Direct Numerical Simulation Comparison and Validation
The PhD of R. Di Battista [18] at CMAP laboratory proposes a methodology to perform post-processing of DNS simulations exploiting geometrical properties of the interface between two phases. In particular, a library called Mercur(v)e has been developed and allows the characterization of two-phase interface geometrical features thanks to a triangulation of the interface and the calculation of geometrical properties from the discretized surface, while preserving topological and geometrical invariants of the considered objects. We rely on this tools to study droplet pulsation or collisions. We propose in this work to apply the same methodology in order to investigate small axisymmetric ellipsoidal deformations of a liquid droplet in order to gain insight on the evolution of the droplet geometrical features. Therefore, we will first perform DNS simulations of the axisymmetric ellipsoidal deformation of a spherical droplet with ARCHER and then use Mercur(v)e to compute its geometrical features.

Two-phase Interface Model and Discretization Strategy of the code ARCHER
The code ARCHER is developed at the CORIA laboratory. It has been successfully used for simulating the atomization of a liquid-jet under a realistic diesel injection configuration [19]. The two-phase medium considered in ARCHER is composed by two incompressible fluids separated by a sharp interface. The position of the interface and the description of its geometrical characteristics is ensured by a Coupled Level Set/Volume of Fluid (CLSVOF) method: a Level Set function is used for evaluating the normal and the curvature of the interface while the Volume of Fluid approach ensures the mass conservation of each component. A key element of the numerical strategy implemented in ARCHER relies on solving a Poisson equation for the dynamic pressure that is performed using a MultiGrid preconditioned Conjugate Gradient algorithm (MGCG) [20] coupled with a Ghost-Fluid method [21] to take into account the pressure jump due to the presence of surface tension. The time-integration is performed with a second-order Runge-Kutta scheme. For more information about the ARCHER solver, we refer the reader to [19,22,23,24].

Interface triangulation and geometrical properties estimation in Mercur(v)e
Mercur(v)e library takes as input the volumetric Level-Set field provided by ARCHER and performs an isocontouring procedure using a Marching Cube approach [25,26], generating a triangulated and discretized representation of the interface. For each vertex of the triangles composing the discretized interface, geometrical properties are estimated as 1-ring surface-averaged quantities as presented in [27], where the 1-ring is a portion of surface around a vertex that is equivalent to a Voronoi region if the triangles composing it are non-pathological. For pathological cases an error-bounding fix is possible. This computational approach guarantees the conservation of the Gauss-Bonnet theorem in discrete form. For more information the reader can refer to [18,28].

Direct Numerical Simulation configuration
The computational domain is a box of dimension 40 × 40 × 40 µm that is discretized over a regular Cartesian grid with 128 points per dimension leading to a total of 2 097 152 computational nodes. Symmetric boundary conditions are imposed on each face of the domain. The kinematic viscosity of the fluid is set to ν = 1.7e−8 m s −2 . We consider as initial condition a deformed droplet with an axisymmetric ellipsoidal shape that will evolve due to surface tension effects: the initial shape is defined thanks to the zero-level of a Level Set such that the semi-major axis p over semi-minor axis q ratio is p/q = (1.08, 1.15) and such that the volume of the droplet is 4 3 πR 3 with R = 10 µm. The computational load for a single run is 1280 h CPU which represents 20 h of computation over 64 CPU. We have chosen two deformation amplitudes, one will still be in the small perturbation limit leading to harmonic oscillations of most of the considered quantities, while the second will allow to go beyond this regime and discriminate the proper variables.

Results
One semi-period T of evolution of the droplet deformation is shown in Figure 2. At t/T = 0, Figure 2a, the droplet is initially at a prolate state. Due to capillarity, the droplet potential energy starts to be converted into kinetic energy such that the droplet reaches a spherical shape in Figure 2b, and then continues to deform such that at t/T = T /2 , Figure 2c, the droplet is in an oblate configuration. We now investigate the evolution of geometric quantities over the whole simulation time and integrated over the surface of the droplet for two different initialization of the deformed droplet, p/q = 1.08 and p/q = 1.15 which correspond to η = 0.08 and η = 0.15.
Figures 3b and 5b display the relative mean Gauss curvature ratio, with S 0 G 0 the integral of the Gauss curvature over the spherical droplet, S 0 G 0 = 4π. It is found to be constant and close to zero with three orders of magnitude less than the amplitudes of the other quantities. This result was expected since from the Gauss-Bonnet theorem, G is a topological invariant, but it attests the reliability of the post-processing library Mercur(v)e. Let us underline from the various evaluations, that it is essential to use the Mercur(v)e library for two reasons : 1-we can guarantee the proper evaluation of the topological invariant through a discrete verification of the Gauss-Bonnet theorem, 2-the level of error obtained from the library compared to the evaluation conducted in ARCHER directly from the Level-Set function is much lower and it is essential in this configuration of small perturbation in order to draw some firm conclusions from the investigation. Depending on the cases, since the references quantities used for the relative quantity ratio are exact characteristics of the sphere, mainly functions of its radius, we can have a slight shift in the quantities represented as observed in Figures 3a and 3c. It is due to the level set formulation and to its discretization on the proposed grid, which is not guaranteed to preserve exactly the volume V 0 or the surface S 0 of the sphere.
Figures 3a and 5a show the relative surface ratio evolution over the simulation time, with S 0 = 4πR 2 the theoretical surface of the droplet when reaching a spherical shape. As we can see, once the droplets is left unconstrained in its perturbed configuration, its surface evolves as a harmonic oscillator, as well as the mean curvature, even if we can detect a small departure from harmonic oscillations. In the limit of small amplitude oscillations, where this is valid, we can conduct a direct comparison of the potential energy U with the analytical expression in eq. (3); this is represented in Figure 4, where, up to the initial shift already mentioned and explained, as well as up to the small viscous dissipation still present in the numerical simulations, we have a very good agreement between the theory and the DNS. However in this small amplitude limit, all the variables are quasi-harmonic. This is the reason why we have provided a second case, with larger amplitude, where this departure is more pronounced.
Actually, as plotted in Figures 5c, the relative mean Mean curvature does not show the same behavior for a larger deformation. The integral of the Mean curvature of the droplet at its spherical state divided by the surface reads H 0 = −1/R. Two non-identical periodic peaks are clearly visible, suggesting the surface averaged mean curvature is not harmonic, thus disqualifying the potential energy expressed in terms of H as the proper harmonic oscillator. It confirms that the energy term associated to the "loss of sphericity" expressed in terms of the interfacial density area σ(Σ − Σ 0 ) is reliably representing the underlying harmonic oscillation phenomenon as observed in Figures 5a. This is also coherent with what has been observed numerically and experimentally for bubbles [16], where the harmonic oscillation behavior is very robust, even within the framework of large deformations, an essentiel feature of the model we aim at building.
Eventually, Figure 6 emphasizes the volume conservation over time of the deformed droplet by plotting the relative volume ratio where V 0 is the volume of the droplet when reaching it spherical shape, V 0 = 4/3πR 3 . Again, the order of magnitude is 3 times less than the amplitude of the surface, and similar to the relative ratio of the mean Gauss curvature. These residuals quantify the error made by the post-processing library.  As mentioned before and relying on an analytical derivation presented in Appendix A, we have identified a quasi-invariant quantity, holding exactly in the limit of small perturbations: I = Σ + 3 H α. In Figure 7, we plot the corresponding extensive quantity over the simulation time obtained with the DNS.

Two-scale kinematics Two-phase flow model
We will now follow the lines proposed in [10,29] for instance in order to elaborate a two-phase model that accounts for both large scale kinematics and small scale kinematics using the results of section 1.
We consider two compressible fluids k = 1, 2. Each of them is equipped with a barotropic equation of state (EOS) defined by ρ k → f k (ρ k ), where f k and ρ k are respectively the specific barotropic energy and the density of the fluid k. The partial pressure p k and the sound velocity a k of the component k are defined by p k = ρ 2 k ( df k /dρ k ), a 2 k = dp k /dρ k . We suppose that there is a large scale kinematic equilibrium in the sense that both fluid k have the same bulk velocity v. We note α k the volume fraction of the fluid k and we suppose that two-component medium to be immiscible so that we have α 1 + α 2 = 1. The density ρ of the overall medium is defined by ρ = α 1 ρ 1 + α 2 ρ 2 and the mass fraction y k of the fluid k verifies ρy k = α k ρ k . We set α = α 1 = 1 − α 2 and Y = y 1 = 1 − y 2 . We also suppose that two additional parameters pertaining to the geometry of the interface are defined on the whole computational domain: (x, t) → Σ and (x, t) → H that respectively denote the density of interfacial area of a spray of identical (same area, volume, curvatures) droplets and the mean curvature of the interface. In the following, we will abuse notation by using H to denote H for sake of simplicity.

Hypotheses related to the small scale kinematics
As the droplets are deforming, the variation of their surface, the mean curvature and the volume are connected. In our case, the deformation of the inclusion is expressed by the variations of the parameters Σ, H and α, so that the variations of the latter shall not be independent. In [2], the authors considered a single surface modified by a small displacement nδh oriented along the normal n of interface. In that case, the infinitesimal variations δΣ, δH and δα verify δΣ = −2HΣδh, δα = Σδh.
If the variations of H along the streamlines can be neglected, that is to say D t H = 0, then one obtains D t (Σ + 2Hα) = 0, which means that Σ + 2Hα remains constant along the streamlines.
In the present configuration, the deformations are no more normal to the interface, preventing the use of above geometric constraints. Furthermore in agreement with the volume conservation showed in Figure 6, we suppose that the small scale motion of the droplets is incompressible. As a consequence, if no breakup occurs, the volume of each droplet remains constant. It first implies that D t α = 0, but also that the volumetric density number of the droplet is constant, allowing us to apply the results of Section 1 valid for a single droplet to a spray of droplets. We thus propose to consider that the small scale kinematics of the droplets is now constrained by the relations Thanks to the results of the previous section 1. 3, we can single out parameters that show the oscillating behavior of the system and parameters that are invariant during the small scale motion. According to Figures 3 and 5, the oscillating profile obtained for the density of interfacial area Σ suggests to drive the motion of the small scale using Σ. We express this choice in the definition of our system by considering a small scale kinetic energy K and a potential energy pertaining to small scale capillary effects U as follows where ν Σ = ν Σ (Σ) > 0 and the capillarity coefficient σ > 0 is supposed to be constant. Let us underline that ν Σ can be interpreted as a micro-inertia related to the small-scale momentum along the lines of [10].

Lagrangian energy of the two-phase system
We proceed following [2,10] by postulating a Lagrangian energy L associated with our system. This energy is the difference between the kinetic energy of the system and the potential energy. Concerning the kinetic energy, we simply add the classic kinetic energy associated with the bulk kinematics 1 2 ρ v 2 and the small scale kinetic energy K . We assume that the medium is equipped with a specific bulk potential energy of the form f (ρ, Y, α). This suggests to postulate the following form for L Following the lines of [2,10], we consider the motion of a portion of medium that occupies the volume B(t) for t ∈ [t 0 , t 1 ]. We note (t, X) → ϕ L (t, X) the position of at instant t of a fluid particle that was initially located at X and set Ω = A transformation of the medium will then be defined by the fields (t, x) → (ρ, v, Y, α, H, Σ) and the associated mapping (t, X) → ϕ L (t, X) so that they verify (11) supplemented by the small scale kinematics (8) hypotheses and the following mass conservation equations We now consider a family of transformations medium (t, , λ ∈ [0, 1]} and we suppose that these fields also verify (11), the small scale kinematics hypotheses (8) and the mass conservation equations, that is to say supplemented with the classic boundary constraints Following standard lines, this family of transformation yields a family of infinitesimal transformations defined by for b ∈ {ρ, Y, α, Σ, H}. Applying (15) with the constraints (13) allows to express following relations between the infinitesimal variations

Extremization of the action
Let us now define the Hamiltonian action A associated with the family of transformations (t, x, λ) → (ρ λ , Y λ , α λ , Σ λ ) and (t, X, λ) → ϕ L λ by setting The Least Action Principle states that a physical transformation of the system verifies dA dλ (0) = 0. (18) which will yield the motion equations of the flow. We introduce the following notations Relation (18) reads We follow standards lines: taking into account (16), (14c) and (14d), after thorough calculations relation (20) provides for any δ λ ϕ and δ λ α. It implies that Finally, by gathering the mass conservation equations (12), the small scale kinematics equations (8) and the evolution equations derived by the stationnary action principle, we see that the full system that governs our two-phase model is

Companion Conservation Equation
Manipulating evolution equations (21a) and (21b), we obtain Defining a total energy ρE as

Final system
In order to express the system of equations of our model in a way that is reminiscent of other systems of the literature [10,31,32] we start by making explicit the derivative of L defined by (10). We have We define the pressure p of the two-phase flow medium by and introduce new variable ω Σ that accounts for the pulsation effects by setting Injecting these relations into Equations (22), we obtain the system ∂ t ρ + div (ρv) = 0, (29a) with the extended pressure P = −L = p + 1 2 ν Σ (ρY ω Σ ) 2 − σΣ. The supplementary conservative equation on the total energy E reads with ρE = 1 2 12 The Least Action Principle does not account for irreversible process. Thus, system (29) is purely nondissipative. We shall propose in the next section additional terms that account for inner dissipation processes that are compatible with the Second Principle of Thermodynamics.

Dissipation
Following the lines of [10], we seek an entropy flux G and additional terms in the evolution equation of ω Σ and α so that we can obtain an evolution equation for E of the form or equivalently, noticing that D t Y = 0, Accounting for the mass conservation and the following partial derivatives of E by injecting into (33), we obtain A simple closure consists in choosing and setting where Σ and µ are both positive coefficients.
We finally obtain a new equation for ω Σ and α that takes the form To close the investigated system, we need to specify the thermodynamics of the mixture. We propose to choose the medium barotropic energy as follows where ∆f is a mixing energy that accounts for phase interaction [33,34]. We further assume that the phase internal energies verify the Gibbs equation. In this case, we obtain 13 By noticing that the velocity evolution equation can be expressed as a conservation equation for the momentum ρv, the overall system equipped with its dissipative part reads ∂ t ρ + div (ρv) = 0, (42a) Let us underline that the dissipation coefficients can both be interpreted as a micro-viscosity for the pressure relaxation in this spray configuration along the lines of [10] and, in the same manner, the oscillation microviscosity Σ can be related to the classical dissipation process of a single droplet oscillation described for example in [11,12].

Properties of the system
Considering smooth solutions for one-dimensional problem, the pure convective part of (42) is obtained by canceling the source terms and reads The

Conclusion
In the present work, we have recalled and extended the analytical study of [11] within the framework of a spherical liquid droplet subject to an axisymmetric ellipsoidal deformation. We have in particular verified that the system is experiencing a harmonic oscillator behavior with identified frequency and potential damping and have shown that in the small perturbation limit, several candidate in terms of small scale variables can be envisioned since most of the variable experience a quasi-harmonic dynamics in that limit.
In order to discriminate the proper variables, we have conducted a series of DNS configurations with variable deformation amplitudes relying on the code ARCHER for high-fidelity simulations. The evolution of the geometrical characteristics of the droplet surface have been analyzed thanks to a post-processing library Mercur(v)e introduced in [18]. The proper and accurate evaluation of the geometrical properties of the deformed spherical droplet have shown that the excess surface is the eligible variable to relate accurately the harmonic oscillation of the droplet in a large range of amplitudes and thus offer a clear picture on what should be the functional dependency of the potential energy, as well as kinetic energy in order to describe the small scale.
Then following the lines of [2,10], we have proposed a model for disperse flow that can account for small scale oscillations of the liquid inclusions by accounting for small scale kinematics and small scale surface tension effects based on the interfacial density area. One key issued of the derivation is the constraints on the geometric variables. The hypotheses, holding exactly in the limit of the small perturbation, were analyzed through our DNS results and validated. The extension of the present study to large scale and small scale capillarity is the subject of our current research and is conducted in [13].  and reducing the range of parameter θ to θ ∈ [0, π/2] thanks to system's symmetry, the deformation d(θ) may be expressed as DenotingX = X/R the scaling by droplet radius R, the deformation reduces tō where is the radial coordinate in the spherical coordinate system and e is defined as e = e r + tan θe z . For small deformations, the ellipsoid remains close to a sphere and the aspect ratio may be expanded as = 1 +η.
For now we supposeη of small maximal amplitude η with no hypothesis concerning its time dependence. The Taylor expansion of the θ dependance ofd writes cos(θ) 2 + tan 2 θ = cos θ 2η +η 2 + 1/ cos θ 2 = η→0 1 +η cos 2 θ +η Under the hypothesis of incompressible liquid, the volume of the deformed sphere must remains equal to the volume of the sphere, inducing As a consequence, the expression of deformationd (θ) may be expanded with respect to the droplet perturbation η asd (θ) = where Y 2 represents the spherical harmonic Y 0 2 [35]. The initial condition of the weakly deformed ellipsoid droplet is thus similar to the analytical framework developped in [11] such as similar dynamics can be expected asymptotically.
Approximating theη volume-preserving ellipsoidal perturbation by the 2ηY 2 /3 spherical harmonic rises the question of volume conservation. Indeed, performing the exact integration on a 2ηY 2 /3 perturbed droplet lead to the following result sin ϑdϑ, Volume conservation is only valid at first order in the perturbation η -this comes as no surprise given the linearisation -and thus must be, if necessary, enforced a posteriori. This will be done by imposing a spherical radius of the deformed sphere R d slightly smaller than the reference R, This leads to the volume preserving 2ηY 2 /3 expression for the deformed surface r(ϑ, t) = R 1 + 4 15η 2 + 16 945η This approximation is compared to the reference sphere and the exact ellipsoid on Figure 9.

A.2. Kinetic energy
Approximating the initial ellipsoidal droplet by a sphere perturbed by an Y 2 spherical harmonic allows the use of the analytical expressions developped by Chandrasekhar [11] which are briefly recalled hereafter.
In the case of weakly spherical harmonic deformed bubble of non-viscous fluid of density ρ in a vacuum forcefree surroundings only subject to capillarity effects with capillarity coefficient σ, the linearized equation of fluid dynamics can be explicitly solved and yields the following surface tension induced oscillations v ,l = − σ ρ l(l − 1)(l + 2)R 3/2−l−2 A l−1 Y l sin(ω l t), compatible with the given radial deformation r(ϑ, t) = R + AY l (ϑ) cos(ω l t). Each chosen spherical harmonic Y l then selects a precise pulsation given by the relation In the case of the ellipsoid surface approximation, the spherical harmonic Y 2 is of special interest and yields ω 2 = ω = 8σ ρR 3 , and we will supposeη = η cos(ωt) to comply with Chandrasekar's hypothesis set so that the surface deformation will satisfy r(ϑ, t) = η→0 R 1 + 2η 3 Y 2 (ϑ) cos(ωt) + o(η).
Writingη = η sin(ωt) and S 0 = 4πR 2 , we can integrate the velocity field to obtain the kinetic energy

A.3. Potential energy
The dynamics of the oscillating droplet is obtained by Chandrasekhar by considering the Laplace pressure given by the deformation as a boundary condition for the pressure field. Our purpose is here to underline Considering the potential energy derived in terms of the surface, an expression of potential energy in terms of the Mean curvature both complying to the analytical results of Chandrasekhar and the H-based approach of Herrmann [8] would thus write,