A novel momentum-conserving, mass-momentum consistent method for interfacial flows involving large density contrasts
HHighlights
A novel momentum-conserving, mass-momentum consistent methodfor interfacial flows involving large density contrasts
Sagar Pal, Daniel Fuster, St´ephane Zaleski • Conservative formulation of Navier Stokes with interfaces using the Volume-of-Fluid method. • Geometrical interface and flux reconstructions on a twice finer grid en-abling discrete consistency between mass and momentum on staggereduniform Cartesian grids. • Conservative direction-split time integration of geometric fluxes in 3D,enabling discrete conservation of mass and momentum. • Quantitative comparisons with standard benchmarks for flow configura-tions involving large density contrasts. • High degree of robustness and stability for complex turbulent interfacialflows, demonstrated using the case of a falling raindrop. a r X i v : . [ phy s i c s . c o m p - ph ] J a n novel momentum-conserving, mass-momentumconsistent method for interfacial flows involving largedensity contrasts Sagar Pal a, ∗ , Daniel Fuster a , St´ephane Zaleski a a Institut Jean le Rond ∂ ’Alembert, Sorbonne Universit´e and CNRS, Paris, France Abstract
We propose a novel method for the direct numerical simulation of interfacialflows involving large density contrasts, using a Volume-of-Fluid method. Weemploy the conservative formulation of the incompressible Navier-Stokes equa-tions for immiscible fluids in order to ensure consistency between the discretetransport of mass and momentum in both fluids. This strategy is implementedon a uniform 3D Cartesian grid with a staggered configuration of primitive vari-ables, wherein a geometrical reconstruction based mass advection is carried outon a grid twice as fine as that for the momentum. The implementation is in thespirit of Rudman (1998) [41], coupled with the extension of the direction-splittime integration scheme of Weymouth & Yue (2010) [46] to that of conservativemomentum transport. The resulting numerical method ensures discrete con-sistency between the mass and momentum propagation, while simultaneouslyenforcing conservative numerical transport to arbitrary levels of precision in3D. We present several quantitative comparisons with benchmarks from exist-ing literature in order to establish the accuracy of the method, and henceforthdemonstrate its stability and robustness in the context of a complex turbulentinterfacial flow configuration involving a falling raindrop in air.
Keywords: large density contrasts, consistent discrete mass-momentumtransport, sub-grid advection, geometric Volume-of-Fluid, staggered Cartesiangrids, 3D conservative direction-split advection
1. Introduction
The dynamics of liquid-gas interfacial flows play a critical role in severalprocesses in nature, as well as in myriad industrial applications. The key ele-ments of such flows are droplets and bubbles, that constitute the fundamental ∗ Corresponding author: Sagar Pal
Email addresses: [email protected] (Sagar Pal), [email protected] (Daniel Fuster), [email protected] (St´ephane Zaleski)
Preprint submitted to Journal of Computational Physics January 13, 2021 echanisms governing the exchange of heat and mass at the ocean-atmosphereinterface [12, 42], mixing/separation in the context of metallurgical processes[21, 49], conventional modes of heat transfer [11, 22] and ever so importantly,the transmission of pathogens [6, 16].A substantially large subset of all surface tension dominated flows involvesignificant disparities in the material properties across the interface, the mostcommon example being flow configurations corresponding to air-water systems,where the densities and viscosities of the fluids are separated by (approximately)3 and 2 orders or magnitude, respectively. The development of numerical meth-ods that attempt to model such interfacial flows involving marked contrasts indensity face several challenges, key amongst them being the transport of math-ematical discontinuities that arise out of the aforementioned contrasts. Ex-tremely small numerical errors are ubiquitous as a consequence of the numerousapproximations involved at each and every step of the algorithm (e.g. inter-face reconstruction, flux computation etc.) . In the context of such flows, such“numerical errors” may result in physically inconsistent mass and momentumtransfer across the interface, often from the denser phase towards the lighterphase as a consequence of inadequate numerical resolution. The presence oflarge density contrasts tend to amplify the growth of these cascading numericalerrors, eventually leading to significant (often catastrophic) interfacial deforma-tions, followed rapidly by a loss of numerical stability.Since such numerical instabilities were first observed with the “SURFER”[23] code in the context of planar jets and rising bubbles, considerable effortshave been made towards the design of numerical methods to specifically dealwith flows involving such marked density contrasts (i.e. large | log r | , where r is the ratio of the densities of the two fluids). The underlying principle behindthese endeavours is that the governing equations for the transport of mass andmomentum are solved using a conservative formulation (divergence of fluxes),instead of standard non-conservative forms, which themselves were adapted di-rectly from techniques developed originally for single-phase flows. This for-mulation enables one to render the discrete transport of momentum consistent with respect to the discrete transport of mass. Such a tight coupling of thepropagation of errors between the discrete mass and momentum fields enablesalleviation of many of the issues that plague such numerical methods, especiallyin the context of low to moderate spatial resolutions.Exacerbating the already complicated nature of the discrete transport ofmaterial discontinuities is the role of surface tension on the evolution of theinterface. They are commonly modeled as singular source terms in the mo-mentum balance equation that governs the evolution of the velocity field, withthe capillary force itself being proportional to the third derivative of the in-terfacial position. Thus, a secondary but important mitigating factor is theadvancements made in the modeling of capillary forces, resulting in the adop-tion of consistent and well-balanced surface tension formulations. Consistencyin the context of surface tension models refers to the ability of methods to pro-2ressively achieve more accurate estimations of interfacial curvature as a resultof increasing spatial resolution, whereas well-balanced refers to the ability torecover certain static equilibrium solutions pertaining to surface tension domi-nated flows without the perpetual presence of parasitic or spurious currents inthe velocity fields. We refer the reader to influential works of Popinet [34, 35]to get a better understanding of the issues surrounding different surface tensionimplementations.The two most popular approaches that deal with interfacial transport arethe Volume-of-Fluid (VOF) method first developed by Hirt and Nichols [20],and the level set class of methods pioneered by Osher and Sethian [29]. Eachclass of methods has its own set of merits (and demerits) relative to each other.Generally speaking, Volume-of-Fluid based methods display superior mass con-servation whereas in terms of interface curvature computation, level set basedmethods hold an advantage. For a more detailed and nuanced evaluation of thecomparitive advantages of interfacial transport methods, we refer the reader tothe recent review by Mirjalili et al. [27] on the given subject. An additionalfeature one has to consider in the context of Volume-of-Fluid based methods isthat of conservative mass transport in 3D. Fluxes can be computed via alge-braic transport schemes (generally less accurate), or by using geometric recon-structions in either Eulerian, Lagrangian or hybrid frameworks. The temporalintegration of the fluxes could be carried out either as a series of one dimen-sional propagations along each of the spatial directions, termed as direction-split, or carried out in one single sweep, termed as multidimensional or unsplit.Direction-split methods are more intuitive and easier to develop (extension to3D in particular), but generally suffer from lack of conservation (to the orderof machine precision) when it comes to 3D. A notable exception to the abovepoint is the conservative algorithm originally developed by Weymouth & Yue[46], which maintains conservative direction-split VOF transport in 3D (subjectto local CFL restrictions). Multidimensional (unsplit) methods have an advan-tage in that respect due to the fact that they are conservative by nature of theirdesign, but are inherently more complicated to develop and implement, with nostraightforward extension from 2D to 3D.The first study to address the issue of consistency between mass and mo-mentum transport was the seminal work of Rudman [41]. The fundamentalhurdle in the implementation of mass-momentum consistent transport for stag-gered configurations of primary variables (pressure and velocity) is the inherentdifficulty in reconstructing mass (defined on centered control volumes) and itscorresponding fluxes onto the staggered control volumes on which momentumis defined. Rudman introduced the strategy of carrying out mass advection(although using algebraic flux reconstructions) on a grid twice as fine as thatof momentum, thereby enabling a ‘natural’ and intutive way to reconstructmass and its fluxes onto staggered momentum control volumes. However, themethod uses a VOF based convolution technique for curvature computation,which is neither consistent nor well-balanced. Bussmann et al. [9] were ableto circumvent the issue surrounding staggered grids altogether by using a collo-3ated arrangement in the context of hexahedral unstructured meshes, coupledwith an unsplit Eulerian flux computation method. The study though makesno mention of any surface tension model.Level set based methods in the context of mass-momentum consistent trans-port were implemented first by Raessi and Pitsch [39], followed by Ghods andHermann [15]. In the former, the consistency problem is tackled by means ofa semi-Lagrangian approach, computing geometric level set derived fluxes attwo different time intervals, whereas in the latter, a collocated arrangement isused. Nonetheless, both methods face certain drawbacks, notably the applica-bility only to 2D in case of Raessi and Pitsch, as well as a lack of well-balancedsurface tension models for both these methods. Recent advances concerningVolume-of-Fluid based methods that employ unsplit (conservative) geometricflux reconstructions were made by LeChenadec and Pitsch [25], and later byOwkes and Desjardins [31]. LeChenadec and Pitsch utilize a Lagrangian remapmethod in order to construct consistent mass-momentum fluxes for the stag-gered control volumes, while Owkes and Desjardins use mass advection on adoubly refined grid (same principle as Rudman) to achieve consistency. Al-though [25] implements a well-balanced surface tension model, the VOF convo-lution based curvature computation is not consistent. In case of [31], they usemesh-decoupled height functions to compute curvature while coupling it witha well-balanced surface tension model. However, their semi-Lagrangian fluxcomputation procedure involving streak tubes and flux polyhedra are extremelyconvoluted in 3D.Certain methods attempt to combine the qualities of both Volume-of-Fluidand level set methodologies (CLSVOF), as proposed in the works of Vaudor etal. [45], and more recently by Zuzio et al. [50]. They both tackle the consis-tency issue by means of projecting the direction-split geometric fluxes onto atwice finer grid, which are subsequently recombined to reconstruct consistentfluxes for mass and momentum for the staggered control volumes. This ap-proach allows them to bypass the requirement of conducting mass advectionon a twice finer grid (as in the original Rudman method), thereby deriving thebenefits of a sub grid without the added computational cost of doubly refinedmass transport. In addition, both methods adopt well-balanced surface tensionmodels with consistent level set based curvature estimation. However, the pur-ported advantages of both these methods with regards to reduced computationalcosts is not quite evident, as additional complexities are introduced due to theprojection (reconstruction) of fluxes onto a the twice finer mesh, which wouldnot be necessary in the first place if mass transport had been carried out on thetwice finer mesh itself. Patel and Natarajan [32] developed a hybrid staggered-collocated approach to solve the consistency issue on polygonal unstructuredmeshes, complemented with a well-balanced surface tension model. Neverthe-less, the VOF advection is based on algebraic transport, not to mention theuse of a VOF convolution based curvature computation, which is inherently notconsistent. More recently, Nangia et al. [28] developed a CSLVOF method fordynamically refined staggered Cartesian grids. They utilize Cubic Upwind Inter-4olation (CUI) schemes to reconstruct consistent mass and momentum fluxeson the staggered control volumes, using the information from the additionalmass advection equation they solve alongside the level set function. However,the reconstruction of mass fluxes using CUI schemes are inherently algebraic,with their comparative advantage against fluxes computed via geometric con-structions being an open question (refer to Mirjalili et al. [27]) . In our previousstudy [3], the issue of consistent mass-momentum advection on staggered gridswas tackles by using shifted fractions . Although in this approach the shifted (mass and momentum) fields are transported in a conservative direction-splitfashion, the method is not exactly conservative when it comes to momentum.A comparative summary of all the methods in existing literature can be foundin Appendix A 6 .In the present study, we adopt the strategy of mass advection on a twicefiner grid in order to enable consistent reconstruction of mass and momentumon the staggered control volumes, in the spirit of Rudman [41]. A geometricrepresentation of the interface within the Volume-of-Fluid framework is utilized, along with geometric reconstruction of fluxes whose temporal integration is acarried out in a direction-split manner. The conservative direction-split mass(volume) transport algorithm of Weymouth and Yue [46] is extended to thedirection-split transport of momentum. The implementations of the algorithmsare developed on the free and open-source numerical platform called “PARISSimulator” [2], with the detailed descriptions of the general capabilities of thesolver to be found in the reference provided. The organization of the paper isquite simple, in section 2 we describe our set of governing equations, followedby their numerical implementation in section 3. In the numerical benchmarkssection 4 we assess the accuracy of our consistent and conservative method bymeans of canonical test cases which are well-established in existing literature.Finally we present the numerically challenging case of a raindrop falling inair. In the process, we evaluate the robustness, stability and accuracy of thepresent method compared to the standard version of our method which doesnot ensure consistent mass-momentum transport. We wrap up this study withsome concluding remarks and future perspectives.
2. Governing Equations
In this section, we describe our methodology behind modeling the dynamicsof immiscible incompressible liquid-gas interfacial flows under isothermal con-ditions.
We use the one-fluid formulation for our system of governing equations,thus solving the incompressible Navier-Stokes equations throughout the wholedomain including regions of variable density and viscosity, which itself depend onthe explicit location of the interface separating the two fluids. In the absenceof mass transfer, the velocity field is continuous across the interface at the5ncompressible limit, with the interface evolving according to the local velocityvector. Thus, the equations are ∂ρ∂t + ∇ · ( ρ u ) = 0 , (1) ∂∂t ( ρ u ) + ∇ · ( ρ u ⊗ u ) = −∇ p + ∇ · (2 µ D ) + σκδ s n + ρ g , (2)with ρ and µ being the density and dynamical viscosity respectively. Thevolumetric sources are modeled by the acceleration g , and the deformation ratetensor D used to model the viscous stresses defined as D = 12 (cid:104) ∇ u + ( ∇ u ) T (cid:105) . (3)The term σκδ s n models the surface tension forces in the framework of thecontinuum surface-force (CSF) method [7]. The normal vector to the interface is n , σ the coefficient of surface tension and κ the local interfacial curvature. Theoperator δ s is the Dirac delta function, the numerical approximation of whichallows us to map the singular surface force distribution along the interface ontoits volumetric equivalents. At the incompressible limit, the advection of massgiven by equation 1 can be treated as equivalent to that of the advection ofvolume. Within the framework of interface capturing schemes, the temporal evolu-tion of the interface separating the two fluids can be tracked by the advectionequation ∂χ∂t + u · ∇ χ = 0 , (4)where χ is the phase-characteristic function, that has different values in eachphase. Generally, χ is assigned a value of 0 in one phase and 1 in the other.Mathematically, the function χ is equivalent to a Heaviside function in spaceand time. At the macroscopic length scales under consideration, the interfaceevolution as described by equation 4 is modeled as having infinitesimal thicknessunder the continuum hypothesis. The coupling of the interfacial evolution withthe equations of fluid motion as described in 1 and 2 is provided by ρ = ρ χ + (1 − χ ) ρ , (5) µ = µ χ + (1 − χ ) µ , (6)where ρ , ρ are the densities of fluids 1 and 2 respectively, likewise forviscosities µ and µ . For certain flow configurations, it might be beneficial toopt for a weighted harmonic mean description of the variable dynamic viscosity[4], instead of the weighted arithmetic mean as in equation 6.6 . Numerical Methodology The governing equations are numerically solved using finite volume dis-cretizations on uniform Cartesian grids, utilizing state of the art methods ininterfacial reconstruction coupled with geometric transport of the correspond-ing fluxes, curvature computation and surface tension modeling. A detailedexposition of our class of mass-momentum consistent numerical methods is pro-vided, which are specifically designed to circumvent or suppress the uncontrolledand rapid growth of numerical instabilities that arise when dealing with flowsentailing marked density contrasts .
Our numerical studies are based on the Volume-of-Fluid methodology. Werefer to the discontinuous approximation to the Heaviside function χ as thevolume fraction field or colour function interchageably, which is defined belowin the context of finite volume discretization as C i,j,k ( t ) = 1∆ V (cid:90) ∆ V χ ( x , t ) d x , (7) χ = 1 χ = 0 A B volume fraction fi eld PLIC representation Smooth Discrete
Figure 1: Exlicit definition of the interface location using the volume-of-fluid approach. Thesmooth interfacial representation is discretely represented on the Cartesian grid using the vol-ume fraction (VOF) field, which is defined in (7). The interface is geometrically reconstructedusing the volume fraction field as piecewice continuous line segments (PLIC), details of whichcan be found in [17, 44]. where C is the colour function with its values lying between 0 and 1, i , j and k being the indices to the corresponding discretized control volume ∆ V . Thereare two steps involved in the VOF method, the reconstruction of the interface,and its subsequent propagation (advection). The interface is represented bydisjointed line segments under the PLIC (piecewise linear interface construc-tion) framework. Such reconstructions involve the determination of interfacenormals using the Mixed Youngs Centered method (refer to [17, 44]). Once thegeometric PLIC reconstructions have been carried out, the interface segments7re advected using the the velocity field. This entails computation of the fluxesin a geometric manner, and subsequent integration of the fluxes in a direction-split manner. Towards that objective, we employ the Eulerian implicit schemeoriginally developed by Weymouth and Yue [46] in order to specifically tacklethe problem of discrete conservation when it comes to direction-split schemes.The scheme is basically a forward Eulerian method with respect to the tem-poral integration of fluxes i.e. the fluxes are computed as the quantity of fluidentering or exiting a given control volume through its fixed surfaces as shownin figure We start with the advection of the interface position by the velocityfield u , using the conservative form of (4) as - ∂χ∂t + ∇ · ( u χ ) = χ ∇ · u . (8)As one can observe, the “compression” term on the right hand side of (8)equals to zero in the context of incompressible flows without mass transfer,but it is important to keep this term in our numerical formulation within thedirection-split framework. Integrating this equation in time after carrying outspatial discretization, one obtains C n +1 i,j,k − C ni,j,k = − (cid:88) faces f F ( C ) f + (cid:90) t n +1 t n d t (cid:90) Ω χ ∇ · u d x , (9)where the first term on the right-hand side is the summation over the cellfaces f of the fluxes F ( c ) f . These fluxes of ( u χ ) are computed via geometrical re-constructions, and not via high-order non-linear interpolation schemes which arethe standard in the absence of discontinuities. The definition of the geometricfluxes in mathematical form is F ( C ) f = (cid:90) t n +1 t n d t (cid:90) f u f ( x , t ) χ ( x , t ) d x , (10)where u f = u · n f is the component of the velocity normal to the controlsurface f . Directional splitting results in the decomposition of equation (9) intothree equations, one for each advection substep as C n,l +1 i,j,k − C n,li,j,k = − F ( C ) m − − F ( C ) m + + c m ∂ hm u m . (11)The superscript l = 0 , , C n, i,j,k = C ni,j,k and C n, i,j,k = C n +1 i,j,k . The face with subscript m − is the “left” face in direction m with F ( C ) m − ≥ m +. After each advection substep (11), the interface is recon-structed with the updated volumes C n,l +1 i,j,k , then the fluxes F ( C ) f are computed8 -1, j i , j i +1, j Figure 2: A 2D schemetic of the Eulerian (geometric) flux calculation using the Weymouth-Yue [46] algorithm for the advection substep along the horizontal direction, with the interfacereconstructed using the volume fraction field at the start of the substep. The colour fractionof the central cell ( i, j ) is updated during this substep through the addition of the fluxes(coloured regions), with the green polygon corresponding to the volume entering the cell i, j from the i − , j and the red one corresponding to that exiting i, j into i + 1 , j . The geometricflux calculations are made on the basis face centered velocities of the cell i, j at the start ofthe advection substep. for the next substep. Importantly, we have approximated the compression termin (9) by (cid:90) t n +1 t n d t (cid:90) Ω χ∂ m u m d x (cid:39) c m ∂ hm u m . (12)One must note that in the above equation there is no implicit summationcarried out over m , and the superscript h denotes the spatial discretization of theoperator. On the right hand sides of (11) and (12) the flux terms F ( c ) f and thepartial derivative ∂ m u m must be evaluated using identical discretized velocities.The expression ∂ hm u m is a finite volume approximation of the spatial derivativecorresponding to the m th component of the velocity vector along direction m ,and the “compression coefficient” c m approximates the color fraction. The ex-act expression of this coefficient depends on the advection method and it alsoentails the desirable property of C -bracketing (preservation of 0 ≤ C i,j,k ≤ c m on the Cartesian direction corresponding tothe advection substep might render the discrete sum (cid:80) m c m ∂ hm u m to be non-zero, even if the flow is solenoidal. The primary appeal of the Weymouth-Yuemethod lies in the subtle but important tactic used to estimate the prefactor tothe compression term, with its definition being c m = H (cid:16) C n, i,j,k − / (cid:17) , (13)where H is a one-dimensional Heaviside function. This renders the com-pression coefficient independent of the direction of the advection substep, con-sequently enabling the three discrete directional divergences to sum up to zeroi.e. (cid:80) m c m ∂ hm u m = 0. For a formal proof, we refer the reader to the appendixof the original study [46], which demonstrates that volume conservation is guar-anteed, given certain local CFL restrictions. From a practical point of view, we9an only ensure that they sum up to the tolerance of the Poisson solver (con-vergence criteria usually set between 10 − − − ), and furthermore, the sum isultimately limited by the level of machine precision ( ∼ − − − ) availablefor representing floating point numbers. In the subsections that follow, we detailthe extension of this exactly conservative direction-split transport framework tothe advection of mass and momentum. In order to achieve consistent transport of the discontinuous fields of massand momentum, we start by trying to understand the advection of a generic conserved scalar quantity φ by a continuous velocity field ∂φ∂t + ∇ · ( φ u ) = 0 , (14)where the field φ is smoothly varying except at the interface position, whereit may be discontinuous. The smoothness of the advected quantity away fromthe interface is verified for fields such as density ( ρ ) and momentum ( ρ u ). Amajor theme of this study lies in the search for a scheme that propagates thisdiscontinuity ( φ ) at the same speed as that corresponding to the advection ofvolume fraction ( C ). Temporal integration of the spatially discretized versionof (14) gives us φ n +1 i,j,k − φ ni,j,k = − (cid:88) faces f F ( φ ) f . (15)The sum on the right-hand side is the sum over faces f of cell i, j, k of thefluxes F ( φ ) f of φ . The flux definitions are given by F ( φ ) f = (cid:90) t n +1 t n d t (cid:90) f u f ( x , t ) φ ( x , t ) d x . (16)In order to “extract” the discontinuity we introduce the interface Heavisidefunction χ ( x , t ), thus giving us F ( φ ) f = (cid:90) t n +1 t n d t (cid:90) f [ u f χ φ + u f (1 − χ ) φ ] d x . (17)Therefore, the flux can be decomposed into two components as F ( φ ) f = ¯ φ (cid:90) t n +1 t n d t (cid:90) f u f χ d x + ¯ φ (cid:90) t n +1 t n d t (cid:90) f u f (1 − χ ) d x , (18)where the face averages ¯ φ s , s = 1 ,
2, are defined as10 φ s = (cid:82) t n +1 t n d t (cid:82) f φ u f χ s d x (cid:82) t n +1 t n d t (cid:82) f u f χ s d x , (19)and χ = χ , χ = 1 − χ . The total flux can be rearranged as a sum of theconstituents corresponding to the different “fluids” F ( φ ) f = ¯ φ F ( C ) f + ¯ φ F (1 − C ) f . (20) Moving forward, the density field ρ ( x , t ) follows the temporal evolution of thegeneric conserved quantity (14) by simply setting φ = ρ . At the incompressiblelimit the velocity field is solenoidal (divergence-free), with constant densitiesin each phase. We can extract the density trivially from the integrals (19) to exactly obtain ¯ ρ s = ρ s . The flux definitions corresponding to ρ thus become F ( ρ ) f = ρ F ( C ) f + ρ F (1 − C ) f . (21)Using the above definitions for density fluxes, we use our VOF method toconstruct fluxes of the color function in order to obtain conservative transportfor ρ . In principle, the discretization of (14) while setting φ = ρ should resultin conservation of total mass during the advection step. However, as we canobserve, the discrete transport of the density field (15) does not have the “com-pression term” which is present in the discrete transport of volume fraction (11).Thus we have to make several adjustments in order to render the transport ofmass consistent with that of volume fraction at the discrete level. We start byadding the divergence term on the right hand side of the continuous form ofmass transport (14), which gives us ∂φ∂t + ∇ · ( φ u ) = φ ( ∇ · u ) . (22)This equation can be decomposed into advection substeps corresponding tothe direction-split integration framework as φ n,l +1 i,j,k − φ n,li,j,k = − F ( φ ) m − − F ( φ ) m + + (cid:16) ˜ φ m c (1) m + ˜ φ m c (2) m (cid:17) ∂ hm u m . (23)The term c (1) m = c m is the compression coefficient corresponding to the VOFadvection method, while c (2) m = 1 − c m corresponds to that of the symmetriccolor fraction 1 − C . The fluxes F ( φ ) m ± follow the definition given in (20) , whilethe cell averages ˜ φ ms are defined as 11 φ ms = (cid:82) t n +1 t n d t (cid:82) Ω φχ s ∂ hm u m d x (cid:82) t n +1 t n d t (cid:82) Ω χ s ∂ hm u m d x . (24)The direction-split integration operation for ρ is rewritten as ρ n,l +1 i,j,k − ρ n,li,j,k = − F ( ρ ) m − − F ( ρ ) m + + C ( ρ ) m , (25)where the fluxes are given by (21), and the corresponding compression termis C ( ρ ) m = (cid:16) ρ c (1) m + ρ c (2) m (cid:17) ∂ hm u m . (26)The above equations (25,26) ensures that the discrete directional divergenceterms of the mass (15) and volume fraction (11) transport equations are con-sistent at the discrete level. As we have discussed before, the Weymouth-Yuemethod ensures that the “directional” compression terms eventually cancel uponsummation over the substeps, therefore resulting in the conservation of mass atthe same accuracy as the discrete incompressibility condition (cid:80) m =1 ∂ hm u m = 0is verified. Thus far, we have made the direction-split (conservative) advection of volumeconsistent with that of density (mass). We now consider momentum advectionwithin the framework of the conserved scalar transport, by setting φ = ρu q ,where q = 1 , , ρu qs = ρ s ¯ u q,s , (27)where ¯ u q,s is termed as the “advected interpolated velocity”, whose precisedefinition is given as ¯ u q,s = (cid:82) t n +1 t n d t (cid:82) f u q u f χ s d x (cid:82) t n +1 t n d t (cid:82) f u f χ s d x . (28)As a reminder, the subscript s denotes the phase or fluid, and f represents thenormal components defined on the face centers. Thus, we obtain the expressionfor the direction-split integration of the momentum as( ρu q ) n,l +1 i,j,k − ( ρu q ) n,li,j,k = − F ( ρu ) m − − F ( ρu ) m + + (cid:16) ρ ˜ u mq, c (1) m + ρ ˜ u mq, c (2) m (cid:17) ∂ hm u m , (29)12here the momentum fluxes are constructed using F ( ρu ) f = ρ ¯ u q, F ( C ) f + ρ ¯ u q, F (1 − C ) f . (30)The expression for the “central interpolated velocity” corresponding to theaverages ˜ φ ms of (24) is˜ u mq,s = (cid:82) t n +1 t n d t (cid:82) Ω u q H s ∂ hm u m d x (cid:82) t n +1 t n d t (cid:82) Ω H s ∂ hm u m d x . (31)The superscript m is intentionally omitted for the velocities ˜ u mq in order toavoid cumbersome and complicated notations. As a reasonable approximation,we choose to put ¯ u q = ¯ u q, = ¯ u q, for the “advected interpolated velocity”and ˜ u q = ˜ u q, = ˜ u q, for the “central interpolated velocity”. To add clarityto the notion of “central interpolated velocity”, one can interpret this as theface-centered interpolations of the velocity field (component normal to controlsurface), which is required to compute the volume fluxes. This leads us to animportant simplification given by F ( ρu ) f = ¯ u q F ( ρ ) f . (32)This model is central in our effort towards consistent mass-momentum trans-port. Therefore, the advection substep for the momentum can finally we writtenas ( ρu q ) n,l +1 i,j,k − ( ρu q ) n,li,j,k = − ¯ u q F ( ρ ) m − − ¯ u q F ( ρ ) m + + ˜ u q C ( ρ ) m , (33)where the density fluxes are defined in (21) and the compression term C ( ρ ) in (26). In the above expression, the face-weighted average velocities ¯ u q aredefined using (28) on the corresponding left face m − or right face m +.As we have referred to in section 3.1, the compression coefficient is inde-pendent of the directional substep due to the adoption of the Weymouth-Yuealgorithm. The final expression for the compression coefficient becomes C ( ρ ) m = (cid:16) ρ C n,l + ρ (1 − C n,l ) (cid:17) ∂ hm u m . (34)Since there is no bracketing on any velocity component, we take ˜ u q = u nq ,which is independent of the substep l . The final expression after cancellation ofthe compression terms, having undergone three advection substeps (33) is( ρu q ) n, i,j,k − ( ρu q ) ni,j,k = − (cid:88) faces f ¯ u q F ( ρ ) f . (35)13herefore, the extension of the Weymouth-Yue algorithm for exact mass(volume) conservation has been extended to incorporate the transport of mo-mentum, thus ensuring consistency at the discrete level between mass and mo-mentum transport. The manner in which the velocity interpolations concerningthe weighted averages ¯ u q and ˜ u q are carried out (near the interface and in thebulk), are identical to that covered in our previous study [3], hence in the in-terest of brevity the reader can refer to the appendix of the aforementionedreference. Figure 3: A 2D schematic of the staggered spatial configuration of the pressure and velocityvariables. The pressure p i,j is based on the center of its control volume (light red); thehorizontal velocity component U i +1 / ,j is defined in the middle of the right edge of the pressurecontrol volume and centered on its control volume (blue dashed box); the vertical velocitycomponent V i,j +1 / is defined in the middle of the top edge and centered on its correspondingcontrol volume (green dashed box). We start by describing the spatial arrangement of our primary variablesi.e. pressure ( p ) and velocity u , where U , V and W represent the x , y and z componenets of u respectively. The control volume is in the form of a cube (3D)14r a square (2D). The pressure and velocity variables are defined in a staggeredarrangement, which is illustrated in Figure 3. The use of staggered controlvolumes has the advantage of suppressing neutral modes (pressure oscillations)often observed in collocated methods but leads to more complex discretizations(refer to [44] for a detailed discussion) . In order to describe the overall numericalalgorithm for the one-fluid Navier-Stokes equations with variable density andviscosity, we choose to the reframe our equations in a more convenient operatorform, given as ∂∂t ( ρ u ) = L ( ρ, u ) − ∇ p . (36)The operator L in the above expression can be decomposed as L = L adv + L µ + L σ + L g , (37)where the L adv represents the conservative advection, L µ represents the dif-fusive forces generated by viscous stresses, L σ represents the capillary forcesarising from the surface tension model and finally L g represents the volumetric(body forces) source term. We apply the spatially discretized versions of theseoperators (denoted by the superscript h ) onto the primary variables ( c, u ), andmarch forward in time using a small, possibly variable time-step τ such that t n +1 = t n + τ . The volume fraction ( c ) is defined on the sub-grid (covered inthe next subsection), whereas pressure and momentum are defined at the coursegrid level. In the first part of the algorithm, the volume fraction field c n is up-dated to the next timestep , with the superscript n signifying discretization intime. The operation can be written as c n +1 = c n + τ L h vof ( c n , u n ) . (38)The subscripts i, j, k are dropped in equations (38) and (39), with the un-derstanding that the operators in equation 37 apply identically to all controlvolumes. The temporal evolution of the volume fraction field represented aboveby the operator L h vof is in accordance with the Eulerian implicit Weymouth-Yueadvection scheme, as described in section 3.1. Once we have obtained the up-dated field c n +1 , we can move on to the temporal update of our momentum fieldgiven by ρ n +1 · u ∗ = ρ n · u n + τ L h adv ( C n , u n ) + τ (cid:2) L hµ (cid:0) C n +1 , u n (cid:1) + L hσ (cid:0) C n +1 (cid:1) + L h g (cid:0) C n +1 (cid:1)(cid:3) . (39)The C field is defined on the coarse grid, reconstructed using informationfrom the sub-grid volume fraction ( c ) field. In regions of constant density.the advection operator L h adv is implemented using higher order spatial schemes15oupled with a choice of non-linear flux limiters such as QUICK, ENO, WENO,Superbee, Verstappen and BCG. These high-order spatial schemes are basedon well established methods developed to deal with hyperbolic conservationlaws, for more details refer to the studies of Leveque [26] and Sweby [43]. Forcontrol volumes in the vicinity of the interface location, we revert to lower ordervariants of these schemes due to the sharp jumps in the material propertiesacross the interface (refer to [3]). The functionality of the operator L h adv nearthe interface is tighly coupled to that of L h vof from equation 38, thus representingthe consistent volume-mass-momentum transport scheme which was covered inthe preceding section. The development of the consistent transport schemes should be integratedinto the broader context of our numerical algorithm which deals with the cou-pling between the conservative formulation of the one-fluid Navier-Stokes equa-tions and the geometric transport of the interface. In our particular approach,the difficulty associated with consistent transport on staggered control volumesis resolved by advecting the volume fraction (mass) on a twice finer grid, ,verymuch in the spirit of Rudman’s [41] original work. In the method developedin our previous study [3] which we refer to as the shifted fractions method,flux computations on the staggered cells necessitate another round of interfacialreconstructions based on the shifted volume fraction field after each advectionsubstep, whereas the sub-grid volume fraction information enables us to cir-cumvent the reconstructions and obtain the fluxes directly using simple surfaceintegrals. The sub-grid method allows us to evade complications that arise dueto the parallel evolution of two different volume fraction fields by conductingvolume fraction advection solely on the twice refined grid. From this point on-wards, we use the terms sub-grid and fine grid interchangeably. An additionaladvantage of the having the sub-grid volume fraction field is that it facilitatesa more “natural” computation of not only the staggered volume fraction field,but more importantly its fluxes. In contrast to the shifted fractions methodin which one can use both Lagrangian and Eulerian fluxes, the implementationof the sub-grid lends itself compatible only with an Eulerian flux computationmethod (e.g. Weymouth-Yue scheme). Considering the overall approach, the key differences between our present algorithm and the implementations in theoriginal works of Rudman and Weymouth & Yue are :
Rudman (1998)[41]
The use of geometrical mass (volume fraction) transportinstead of the algebraic transport in the original method.
Weymouth & Yue (2010) [46]
The extension of the original method for con-servative direction-split mass transport to the momentum field, culminat-ing in the discrete consistency between mass momentum transport.
To start, we describe the spatial orientation of the different variables on thetwo refinement levels of the grid. At the coarse grid level, the pressure and16 igure 4: A 2D schematic of the arrangement of primary variables on the coarse and sub-gridlevels. The control volume on the left hand side is identical to the red colored square in Fig. ?? , corresponding to the grid level for the momentum and pressure. Thus, the left hand sidecontrol volume corresponds to the coarse grid, while the one on the right corresponds to thesub-grid . Superposition of these two figures, one on top of each other, represent the spatialrelationships between the two sets of variables. The velocities on the sub-grid level are simplefirst or zeroth order interpolations ((40)) of the coarse grid velocities. velocity fields are defined in a staggered configuration, with pressures at thecell centers (centroids in 3D) and velocities on the cell face centers. In Fig. ,we demonstrate the arrangement of variables, which one can easily extrapolateto 3D during implementation, but for the sake of clarity we illustrate its 2Dequivalent. The volume fraction C is only defined on a grid that is twice asfine than that of the pressure/velocity grid, which we refer to as the sub-grid.Therefore in 3D, each cubic control volume is divided into eight constituentsmaller cubic volumes at the centroids of which, the volume fraction field iscentered. As is paradigmatic for dual-grid methods in the context of Navier-Stokessolvers, we carry out mass transport at the sub-grid level, and momentum trans-port at the coarse level. In pressure-projection based methods such as ours, thepredominant bottleneck in terms of computational speed is the iterative solu-tion to the discrete Poisson equation. The problem size of such elliptical partialdifferential equations is governed by the number of cells/points discretizing thevelocity (or pressure) field. If we choose to solve the problem on a mesh withtwice the resolution in 3D, that would lead to an 8 fold increase in the prob-lem size for the discrete pressure-Poisson problem. Therefore, solving the massadvection equation on the fine grid allows us not only to obtain more accuratesolutions to the flow physics involved, but also enables us to avoid the signifi-cantly higher computational costs associated with solving a Poisson problem ata twice finer resolution. 17 rolongation . In order to advect the volume fraction at the sub-grid level, weneed to reconstruct a velocity field at that level of resolution using informationfrom the coarse grid velocity field. The prolongation operator has to be used atthe start of each new time step, in order to compute the sub-grid velocity field.The relationships between the coarse and fine grid discrete velocity fields as aresult of the prolongation operator are given below u i − , j = u i − , j +1 = U i − ,j ,u i + , j = u i + , j +1 = U i + ,j ,u i + , j = u i + , j +1 = (cid:16) U i + ,j + U i − ,j (cid:17) / ,v i, j − = v i +1 , j − = V i,j + ,v i, j + = v i +1 , j + = V i,j − ,v i, j + = v i +1 , j + = (cid:16) V i,j + + V i,j − (cid:17) / . (40)The notations in the above equations refer to the description of the coarseand sub-grid variables in Fig. The first-order interpolation applied to the coarsegrid velocity field ensures discrete incompressibility at the sub-grid level, whichis a direct consequence of the solenoidal nature of the coarse grid velocity field.The choice of interpolation order used in the prolongation operator is identicalto that of the originial method of Rudman ([41]). The use of higher orderinterpolation schemes would necessitate finding the solutions to local Poissonproblems, in order to render the sub-grid velocity field discretely divergence-free. In the present version of our method, we have chosen to eschew the addedcomplexity of solving local Poisson problems , therefore sticking with the muchsimpler interpolations. Restriction . In order to compute the required density and momentum fieldson the staggered grid (coarse level), we implement restriction operators that useinformation from the sub-grid volume fraction field. The restriction operatorswork in an identical manner as in Rudman ([41]), which are nothing but simplevolume integrals of the sub-grid field mapped onto the domains corresponding tothe coarse grid control volumes. In Fig. , we demonstrate the different volumeintegrals for the centered and staggered fields, as well as the surface integrals fortheir corresponding fluxes. A description of the different functions performedby the restriction operator is as follows : • Fine Grid (cid:55)−→
Coarse Staggered Grid
Density & Momentum
The sub-grid volume fraction field is restrictedto obtain a staggered density field at the coarse grid level (see Fig. 5),which on combining with the appropriate component of the velocityfield produces the staggered momentum field.18 oarse Staggered GridCoarse Centered GridSub-grid staggeredrestrictioncenteredrestriction
Figure 5: A 2D illustration of the restriction operations involved between the different gridresolutions, where the sub-grid grid boundaries are denoted by the dashed lines, and that ofthe overlapping coarse grid by solid lines. Restriction operations in the form of simple volumeintegrals are performed on the sub-grid volume fraction in order to compute the staggered andcentered volume fraction fields at the coarse grid level. Restriction operators in the form ofsurface integrals are also performed in order to compute the volume fluxes at the boundariesof the coarse staggered control volumes. The red shaded areas correspond to the volumefractions, whereas the green shaded area corresponds to that of volume fluxes. The interfacesegments as depicted at the coarse grid level are just for purpose of illustration, and are not‘reconstructed’.
Fluxes
The geometric fluxes from the sub-grid volume fraction field arerestricted in order to obtain coarse grid fluxes for the corresponding19taggered volumes, which on combining with relation (32) gives usthe corresponding momentum fluxes. • Fine Grid (cid:55)−→
Coarse Centered Grid
Density & Viscosity
The sub-grid volume fraction field is restricted tocompute a centered volume fraction field at the coarse grid level,which is subsequently used to derive the centered density and viscos-ity fields at the coarse grid level.
To summarize the operations performed at each time step pertaining to theadvection operator in our one-fluid Navier-Stokes framework, we present Fig. 6,which illustrates a 2D version of the sub-grid method that ensures consistent andconservative mass-momentum transport, highlighting the interactions betweenthe variables defined on the different grids. The algorithm is also summarizedin the following steps :1. Interface reconstruction at time t n , at the sub-grid level, using data c n .2. Restriction of the sub-grid volume fraction in order to compute the stag-gered fraction fields C nq and ρ nq in the coarse grid control volumes, where q = 1 , ρ q u q ) n at t n .4. Advection of the sub-grid volume fraction c n along one coordinate direc-tion, say x direction, to obtain c n + , using the Weymouth-Yue explicitadvection method.5. Advection of all momentum components along the x direction, in sync withthe sub-grid VOF advection, using (33) to obtain the updated momentumcomponents ( ρ q u q ) n + after the first substep.6. Restriction of the sub-grid field c n + as shown in Fig. 5 , in order tocompute the “shifted” density field ρ n + q on the staggered volumes.7. Extraction of the provisional velocity components. u n + q after the firstsubstep, i.e. u n + q = ( ρ q u q ) n + /ρ n + q . These provisional velocities arerequired to compute the “advected” velocities using the non-linear fluxlimiters.8. The above operations are repeated for the different momentum compo-nents, in sync with the direction-split advection of the sub-grid volumefraction along the remaining direction. At each time step, the sequence x, y is permuted, in order to avoid any systematic biases in error propa-gation. Finally, we obtain u ∗ q = ( ρ q u q ) n +1 /ρ n +1 q on the staggered grid.9. In parallel, the updated centered volume fraction field C n +1 is obtained atthe coase grid level, by applying the restriction operator on the updatedsub-grid volume fraction c n +1 . 20 P R R P R R R Figure 6: A bird’s eye view of the present method, highlighting the operations performedat each time step. For sake of clarity, present a 2D case for the density and the horizontalvelocity ˜ U n which is defined on the staggered grid i + 1 / , j . The variables in red are alldefined on the coarse level, and those in green are at the sub-grid level (below the blue dashedline). The operators R and P denote the restriction and prolongation operations respectively.The “tilde” on top of the variables (e.g. ˜ ρ n ) are used to convey that the variables are definedon the staggered control volumes at the coarse grid level, whereas those without it are definedon the centered control volumes at the coarse grid level. The operator Ψ applied to thevelocity field (e.g. Ψ( ˜ U n )) represents the non-linear slope (flux) limiters that operate onlocal stencils , in order to estimate the momentum per unit mass exiting the control volumeduring the advection substep. The central model ((32)) of our mass-momentum consistentmethod is represented in the operation boxes (e.g. Ψ( ˜ U n ) · F ( ρ ) x ), thus allowing us to contructmomentum fluxes from the corresponding mass fluxes. The starting variables are enclosed inthe gray elliptical boxes, and the variables at the end of the direction-split integrations are ˜ U ∗ and ρ n +1 , where ˜ U ∗ is subsequently fed into the RHS of the Poisson problem. The centereddensity field (coarse grid) at the end of the time step ( ρ n +1 ) is used by the surface tensionand viscous diffusion operators. The resulting method not only maintains discrete consistency between massand momentum transport on the staggered grid, but also ensures that the trans-port of momentum is conservative in 3D. The other time-split terms in equation(39) that arise from the source term discretizations, as well as the terms inthe projection step (41) are solved in a standard non-conservative way, and arebriefly covered in what follows. The density on the faces of the central cellsare estimated using the restriction operations described previously, with thesedensities appearing via the 1 /ρ q pre-factor in front of all the the other termsinvolved in momentum transport (e.g surface tension, viscous diffusion).21 .5. Surface Tension and Viscous Diffusion The detailed descriptions of the methods used in our numerical platformPARIS Simulator [2] to deal with surface tension, viscosity and body forceshave already been carried out in [2, 33, 34], therefore we only briefly touch uponcertain aspects of the operators in question, in particular, their interaction withthe volume fraction field.
Surface Tension.
We use the Continuum Surface Force method (CSF) as ourmodel for surface tension, coupled with height functions for curvature compu-tation. The height functions used in our implementation were first introducedin [34] , subsequently tested, revised and improved in [5, 30]. In general, theheight functions are used to compute the curvature field based on second-orderfinite differences applied to the heights. The height functions are derived fromthe centered volume fraction field at the coarse grid level, which itself is com-puted via the previously discussed restriction operations. Although, in regionsof poor interfacial resolution where the local radius of curvature is comparableto the grid size, the method reverts to certain fallbacks such as curve fitting.The resulting curvature field is coupled with a well-balanced discretization withrespect to the discrete pressure gradient, with the same discretization stencilapplied to the volumetric (body) force term as well. Interestingly, the sub-gridvolume fraction can be used to get better estimates of the curvature field, whoseestimates are required on staggered control volumes. At the time of writing, twodifferent approaches are being developed and tested in order to get more accu-rate curvature fields, but in the context of the present study we refrain fromany further complexity.
Viscous Diffusion.
We use second-order spatial discretizations of the viscousstresses, using centered differences. The (variable) dynamic viscosity computa-tions are based on the centered volume fraction field at the coarse grid level, asexpressed in (6). The temporal treatment of the viscous term can be either inexplicit or semi-implicit fashion, , but in the context of the present study wewill be sticking exclusively with the explicit version.
The velocity field is evolved using a classical time-splitting projection methodas described in the seminal work of Chorin [10], which involves predicting an “in-termediate” velocity field u ∗ as given by equation (39), followed by a correctionstep given by u n +1 = u ∗ − τρ n +1 ∇ h p n +1 . (41)The discrete pressure field required to correct the intermediate velocity isdetermined by imposing the conservation of mass, which in our incompressibleframework reduces to necessitating the resulting velocity field to be divergence-free (solenoidal) 22 h · u n +1 = 0 . (42)Thus, combining equations (41) and (42), we are left with a variable coeffi-cient Poisson equation for the pressure, given as ∇ h · (cid:18) τρ n +1 ∇ h p n +1 (cid:19) = ∇ h · u ∗ . (43)The Poisson solver used in PARIS Simulator to invert the elliptic operatorappearing in eqn. (43) is a red-black Gauss-Seidel (GS) solver with overrelax-ation [8]. There is also has an in-house implementation of a multigrid solver forstructured grids with 2 n number of points per direction, utilizing a fully paral-lelized V-Cycle scheme [8]. Relaxation operations are applied starting from thefinest to the coarsest first, and then from the coarsest to the finest, with the num-ber of relaxation operations being a user-adjustable parameter. Having a nativemultigrid solver allows for an efficient solution of the Poisson equation with-out the necessity of having external libraries/pre-conditioners (e.g. HYPRE)installed on the system.
4. Numerical Benchmarks
In this section we carry out quantitative comparisons between the solutionsobtained by our present mass-momentum consistent and conservative method,and certain analytical, semi-analytical, equilibrium or asympotic solutions whichcorrespond to well-established benchmark configurations. Generally, these es-tablished test cases in existing literature (e.g. decay of spurious currents in staticand moving droplets, viscous damping of capillary waves etc.), are carried out inthe absence of any density jump (or viscosity jump) across the interface. Thus,we are distinctly interested in the behavior of our present method as it dealswith the non-linear coupling between interfacial deformation/propagation, cap-illary and viscous forces, especially in regimes where fluid densities are separatedby orders of magnitude across the interface. In the interest of clarity, we shalladopt the following nomenclature from this point onwards • STD : Standard version of our method based on the non-conservativeformulation of Navier-Stokes equations, does not maintain discrete consis-tency between mass and momentum transport. • MSUB : Present method which ensures consistent and conservative trans-port of mass-momentum.
A popular numerical benchmark in the existing literature relevant to surfacetension dominated flows is the case of a spherical droplet of a denser fluid23mmersed in a quiescent surrounding medium of lighter fluid. In the hydrostaticlimit of the Navier-Stokes equations, the droplet should stay in equilibrium,with a curvature induced pressure jump across the interface corresponding toLaplace’s equilibrium condition. In practice however, numerically reproducingsuch a trivial equilibrium condition is not as straighforward, as there exists aslight difference between the initial numerical interface and the exact analyticalshape of the sphere, thereby resulting in the generation of the well documented’ spurious ’ or ’ parasitic ’ currents of varying intensity in the velocity field [19,23, 36]. Considerable progress has been made by the adoption of well-balanced surface tension formulations, that ensure consistency between the numericalstencils used for the discretization of the pressure gradient and nδ s ([14, 34]).In this context, Popinet [34] demonstrated that given sufficient time (orderof viscous dissipation time scales), a well-balanced method will relax to the’ numerical ’ equilibrium shape through damping of the ’physically consistent’numerical capillary waves, therefore allowing us to recover the exact (to machineprecision) Laplace equilibrium condition. Figure 7: Schematic of the static droplet of dense fluid surrounded by a quiescent medium oflighter fluid. A 40 ×
40 grid is employed to spatially discretize the domain.
The key difference in our implementation of this classic test case from that ofPopinet [34] is that we consider the effect of density contrast across the interfaceseparating the fluids. A sharp density jump across the interface may exacerbatethe effect of numerical errors incurred as a result of interfacial reconstructions,curvature estimation and various other truncations. We consider a circulardroplet of size D placed at the centre of a square domain of side L . The densitiesof the heavier and lighter phases are ρ l and ρ g respectively, likewise for theviscosities µ l and µ g , and σ being the surface tension coefficient (fig. 7). Theratio of the droplet size to the box is chosen as D/L = 0 .
4, and the numericalresolution is D/ ∆ x = 16 (∆ x is the grid size). We use symmetry conditionson all sides of our square domain. The problem incorporates two natural time-scales, the capillary oscillation scale and the viscous dissipation scale, which aredefined as 24 σ = (cid:18) ρ l D σ (cid:19) / , T µ = ρ l D µ l . (44)The ratio of these time-scales give us T µ T σ = (cid:112) ρ l σD/µ l = √ La , (45)where La is the Laplace number based upon the heavier fluid. In the presentstudy, we introduce the density-ratio ρ l /ρ g as another important parameter.We rescale our ”parasitic” velocity field using a velocity scale as defined below U σ = (cid:112) σ/ρ l D . (46)Additionally, due to the explicit nature of our surface tension model, thetime-step in our numerical simulation must be smaller than the oscillation periodcorresponding to the grid wavenumber (fastest capillary wave ∼ (cid:0) ρ l ∆ x /σ (cid:1) / ). For the scope of the present study, we shall not consider any viscosity con-trast between the two fluids while varying the density contrast, therefore setting µ l /µ g = 1 for all the cases.In figures 8 and 9, we illustrate the decay of the root-mean-square of thespurious currents as a function of time, in the case of four different density con-trasts, and with three different Laplace numbers for each ratio. Figure 8 refersto simulations carried out with the non-consistent ( STD ) method, and Fig. 9corresponds to that using the present method (
MSUB ). Time is rescaled bythe viscous dissipation scale, and the spurious currents by the capillary velocityscale.We have two main observations, first being the rapid decay of the rescaledspurious currents for all combinations of density contrasts and Laplace numberswithin approximately 0 . T µ . Secondly, there is a slower re-growth of the currentsfor combinations of non-unity density contrasts and large Laplace numbers,for all simulations except those carried out with the present method ( MSUB )With the present method, the once decayed currents keep hovering around levelsof machine precision for remainder of time. Therefore, the present methoddemonstrates the desired performance, systematically for all combinations oflarge density contrasts coupled with large Laplace numbers.Once the solution relaxes to a numerical equilibrium curvature, there stillexists a difference between the numerical curvature and the exact analyticalcurvature corresponding to the spherical (circular) shape. We use the definitionsof the shape errors as introduced in the seminal work of Popinet [34] to assess theconvergence to the exact (analytical) curvature as we increase spatial resolution.The norms are defined as 25 . . . . . . T /T µ − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 Figure 8:
STD : Decay of normalized spurious currents as a function of viscous dissipationtime-scales for different density-ratios and Laplace numbers. The currents seem to initiallydecay quickly for larger density-ratios, and relax to the numerical equilibrium curvature within0 . · T µ . For combinations of large ρ l /ρ g and large La, the spurious currents seem to grow backto an order of magnitude (10 − ) which is quite far from that of machine precision (10 − ). L = (cid:115) (cid:80) i ( C i − C exact i ) (cid:80) i , L ∞ = max i (cid:0) | C i − C exact i | (cid:1) , (47)26 . . . . . . T /T µ − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − − U r m s U σ (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 Figure 9:
MSUB : Decay of normalized spurious currents as a function of viscous dissipationtime-scales for different density-ratios and Laplace numbers. The currents seem to decay veryquickly in the case of large density-ratios, and relax to the numerical equilibrium curvaturewithin 0 . · T µ . For all combinations of ρ l /ρ g and La numbers, the decayed spurious currentsare not observed to grow back as in the cases of STD , and hover around values close tomachine precision (10 − ). where C i is the volume fraction of a cell after the solution has relaxed tothe numerical equilibrium curvature, and C exact i is the volume fraction corre-27 D/ ∆ x − − − − S h a p e E rr o r L ∞ : MSUB L ∞ : STD L : MSUB L : STDSecond Order Figure 10: Second-order spatial convergence for the spurious current error norms correspond-ing to a stringent parameter combination ( ρ l /ρ g = 1000 , La = 12000) . Both of the norms( L ∞ and L ) seem to demonstrate a roughly second order rate of spatial convergence witheach of the methods tested. However, the L errors of MSUB are smaller by approximatelya factor of 2 compared to
STD for all resolutions tested. sponding to the exact circular shape. The fields C i and C exact i correspond tothe coarse grid level, and are obtained via restriction of the sub-grid volumefraction field. Fig. 10 demonstrates the behavior of the shape errors defined in(47) for the case of the most stringent parameter combination ( ρ l /ρ g = 1000 ,La = 12000 ), as a function of the droplet resolution. One can clearly observethat both methods display a roughly second-order convergence in space for theerror norms. In terms of the L norm, the present method ( MSUB ) achievessmaller errors (roughly by a factor of 2) as compared to the non-consistentmethod (
STD ) for all resolutions tested.
An incisive numerical setup which evaluates the accuracy of the couplingbetween interfacial propagation and surface tension discretization was first pro-posed by Popinet [34], subsequently employed in the comparative study ofAbadie et al. [1]. This test differs from that of the static droplet due to thepresence of a uniform background velocity field, therefore serving as a betterrepresentation of droplets in realistic surface tension dominated flows wherethey might be advected by the mean flow.In terms of the Laplace equilibrium, the hydrostatic solution is still valid inthe frame of reference of the moving droplet. The solution in the moving refer-ence frame diverges from that of the static droplet as a result of the continuous28 igure 11: Schematic of the droplet of dense fluid advected in a surrounding medium of lighterfluid. A 50 ×
50 grid is employed to spatially discretize the domain, which is spatially periodicin the direction of droplet advection. injection of noise at the grid scale. This noise originates from the perturbationsto the curvature estimates, which are themselves born out of the interfacialreconstructions needed to advect the interface. This transforms the probleminto that of viscous dissipation with continuous forcing (in the moving referenceframe).In the present study, we evaluate our present method using the advection ofa droplet in a spatially periodic domain in the same setup as [34], but with theimportant difference of including sharp density jumps across the interface.We consider a circular droplet of diameter D placed at the centre of a squaredomain of side L . The densities of the heavier and lighter phases are ρ l and ρ g respectively, likewise for the viscosities µ l and µ g , and σ being the surfacetension coefficient (fig. 11). A uniform (horizontal) velocity field U is initializedon the entire domain. The ratio of the droplet size to the box is D/L = 0 . D/ ∆ x = 20. For reference, Popinet [34] used a resolution of D/ ∆ x = 25 . ×
64. We use symmetry conditions on the top andbottom sides, and periodic boundary conditions along the horizontal direction.The problem is characterized by adimensional parameters based on the heavierfluid , given by La = ρ l σDµ l , We = ρ l U Dσ . (48)In addition to the capillary and viscous time-scales for the static case (eqns.44), we have an additional scale defined as T U = D/U , (49)which is the time-scale of advection. In our subsequent analysis, we shalluse T u and U as the time and velocity scales, repectively. Figures 12 ( STD . . . . . . T /T u − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 Figure 12:
STD : Time evolution of normalized spurious currents as a function of advectiontime-scales ( T u ) for different combinations of density-ratio and Laplace numbers. The currentsseem to hover around 10 − , with a larger Laplace number corresponding to a higher error forall density-ratios. We = 0 . method ) and 13 ( MSUB ) depict the evolution of the root-mean-square (RMS)error of the velocity field in the moving frame of reference, as a function ofdifferent Laplace numbers, spanning over several density contrasts. We againhave a couple of important observations, the first being that spurious currents30 . . . . . . T /T u − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 . . . . . . − − − U r m s U (cid:18) ρ l ρ g = (cid:19) La = 120 La = 1200 La = 12000 Figure 13:
MSUB : Time evolution of normalized spurious currents as a function of advectiontime-scales ( T u ) for different combinations of density-ratio and Laplace numbers. In termsof the errors observed in STD , we observe a decrease of roughly one order of magnitude.Although an upward trend is observed for large Laplace numbers, the growth rate is quite low.The currents seem to hover slightly above 10 − , with larger Laplace numbers correspondingto larger errors for all density-ratios. We = 0 . do not decay to machine precision as in static droplet case for all combinationsand methods tested, instead they oscillate around a mean value of the order of31 − Weber − − − − − M a x i m u m E rr o r L ∞ L We − / Laplace − − − − M a x i m u m E rr o r L ∞ L La / Figure 14: Maximum of the error norms as a function of different
W e and La numbers for ourpresent method (
MSUB ), ρ l /ρ g = 1000 for all the cases presented. (a) Maximum error normas a function of Weber (La = 12000, ρ l /ρ g = 1000). (b) Maximum error norm as a functionof Laplace (We = 0 . ρ l /ρ g = 1000). . − .
01% of the constant field U . The second observation is regarding thesignificantly smaller error (roughly by one order of magnitude) in the case of thepresent method ( MSUB ) when compared to the non-consistent one (
STD ). Asa minor remark, in case of large Laplace numbers, the
MSUB method displaysa slight upward trend in the error evolution, which is not the case in
STD . Thisis not too worrisome as the growth is over a time-scale much larger than T u ,with the faster oscillations corresponding to a time-scale of the order U/ ∆ x . Allof the plots correspond to We = 0 .
4, coupled with an additional simplificationof having equal viscosities across the interface i.e µ l /µ g = 1 . As evindencedby the persistence of these spurious currents due to the addition of grid-levelnoise emanating from interfacial reconstructions, further advancements shouldbe made with respect to the combined performace of the interfacial transport,curvature computation and the surface tension model.In order to evaluate the performance of our methods at different resolutions,we define the errors as the maximum values of the norms L ∞ and L of therescaled field U rms /U over time (5 times T u ). In Fig. 15, we show the scaling ofthe error as a function of spatial resolution for the most stringent case of ρ l /ρ g =1000 , La = 12000. For our present method ( MSUB ), we observe significantlylower maximum errors compared to the non-consistent method (
STD ), but ata cost of slightly less than first-order convergence rate. The overall convergencebehavior of the methods we have tested seem to be consistent with earlier studiesof Popinet [34], although in that study only equal densites across the interfacewere considered.As the final point of inquiry, Fig. 14 shows the influence of the Laplaceand Weber numbers on the behavior of the maximum error norm, carried outfor the largest density contrast ( ρ l /ρ g = 1000). We only present the resultsobtained using the present method ( MSUB ), for a resolution corresponding to32 D/ ∆ x − − − − M a x i m u m E rr o r First Order L ∞ : STD L ∞ : MSUB L : STD L : MSUB Figure 15: First-order (approximately) spatial convergence of the maximum of the spuriouscurrent error norms in the frame of reference of the moving droplet, for the most stringent pa-rameter combination ( ρ l /ρ g = 1000 , La = 12000, We = 0 . MSUB ) leads to significantly lower errors. D/ ∆ x = 25 .
6. We observe that the error (both L ∞ and L ) scales as W e − / over 4 orders of magnitude, which is different from the We − / scaling observedby Popinet [34] (for equal densities). In case of Laplace numbers, the errorsscale as La / over two orders of magnitude, which is identical to the that in[34] (again, for equal densities). One of fundamental features of surface tension dominated flows are the prop-agation of capillary waves. A numerical method should not only be able toadequately resolve, but also accurately emulate the spatio-temporal evolutionof such capillary oscillations. A brief outline on the state-of-the-art numericalimplementations of capillary waves (and surface tension models in general) isprovided by Popinet [35].In the present study, we evaluate the accuracy of our method by comparingwith an analytical solution of damped capillary oscillations. Generally, suchsolutions exist only for extremely small initial perturbations, that too eitherin the inviscid limit (Lamb [24]) or the asymptotic limit of vanishing viscos-ity (Prosperetti [37, 38]). For our purposes, we use the configuration of theviscosity-damped capillary oscillations of a planar interface, as implemented byPopinet & Zaleski [36]. 33 igure 16: Schematic of the initially perturbed planar interface separating two immisciblefluids of different densities and viscosities. A spatial resolution of 32 ×
96 is used for spatialdiscretization (compared to 64 ×
192 in Popinet [34]), with the width of the box correspondingto the size of the perturbed wavelength.
We consider a rectangular domain of dimensions L × L , where L correspondsto the wavelength of our initial perturbation. The densities of the heavier andlighter phases are ρ l and ρ g respectively, likewise for the viscosities µ l and µ g ,and σ being the surface tension coefficient (fig. 16). An intial perturbationamplitude of L/
100 is used, coupled with a numerical resolution of L/ ∆ x = 32. Symmetry conditions are applied on the top and bottom sides, with periodicconditions along the horizontal direction. The problem is characterized by theadimensional parameters given as T = T ω , La = ρ l σLµ l , (50)where La is the Laplace number based on the heavier fluid, and ω is definedusing the dispersion relation [34] given as ω = σk ρ l , where k = 2 πL . (51)The above relation is obtained using linear stability analysis at the inviscidlimit [24]. In order to evaluate the influence of density contrast on the per-34 T ω . . . . . . A / A ( ρ l / ρ g = ) SimulationTheory . . . . . . A / A ( ρ l / ρ g = ) SimulationTheory . . . . . . A / A ( A i r - W a t e r ) SimulationTheory
Figure 17:
STD : Time evolution of the amplitude of the planar interface undergoing dampedcapillary oscillations, comparing the solution obtained by the non-consistent method with theclosed-from Prosperetti solution. A relatively good agreement with theory is observed for allexcept the largest density-ratio tested. formance of our method, we use 3 different numerical setups keeping the sameLaplace number (La = 3000) • ρ l /ρ g = 1 , µ l /µ g = 1 (Popinet [34]). • ρ l /ρ g = 10 , µ l /µ g = 1 . • ρ l /ρ g = 1000 . / . µ l /µ g = 1 . · − / . · − (Air-Water) .The final setup corresponds to that of an air-water interface, which is alsothe most stringent due to the large density and viscosity contrasts.The analytical solution for this configuration corresponds to closed-form ex-pressions of the planar interface shape evolution given by Prosperetti [37, 38],35 T ω . . . . . . A / A ( ρ l / ρ g = ) SimulationTheory . . . . . . A / A ( ρ l / ρ g = ) SimulationTheory . . . . . . A / A ( A i r - W a t e r ) SimulationTheory
Figure 18:
MSUB : Time evolution of the amplitude of the planar interface undergoingdamped capillary oscillations, comparing the solution obtained by the consistent and conser-vative method with the closed-from Prosperetti solution. A marginally better agreement withtheory is observed for all the density-ratios tested, compared to
STD especially for the largestdensity-ratio. which takes into account the finite time-scales at which the vorticity (generateddue to interface oscillations) diffuses into the bulk medium. These expressionsare numerically solved using using a fourth-order Runge-Kutta time integrator.The amplitude is normalized by the initial value ( A ) and the time rescaled by T . As we can in figures 17 ( STD ) and 18 (
MSUB ), both numerical methodsprovide relatively good agreement with the Prosperetti solution (black curves).There is hardly any appreciable qualitative difference between the results ob-tained via the different methods
STD and
MSUB , although the present method(
MSUB ) seems to perforn marginally better for the most stringent case (air-water configuration). Next, we evaluate the accuracy of our numerical results36 L/ ∆ x − − − − L E rr o r STD : ρ l /ρ g = 1 STD : ρ l /ρ g = 10 STD : Air-WaterMSUB : ρ l /ρ g = 1 MSUB : ρ l /ρ g = 10 MSUB : Air-WaterSecond Order
Figure 19: Comparison of spatial convergence for all three setups, between the non-consistentand the present method. Both methods seem to demonstrate approximately second-orderconvergence. The consistent and conservative method (
MSUB ) delivers slightly lower errorsin all the setups tested, although there is some saturation in the convergence rate at higherresolutions for the stringent air-water setup. to the Prosperetti solution using an integral (in time) error norm, the same asdefined in Popinet [34], given as L = 1 L (cid:115) ω (cid:90) Tt =0 ( h − h exact ) , (52)where h is the maximum inteface height obtained using our numerical simula-tions, and h exact being the maximum height obtained via time integration of theProsperetti solution. In Fig. 19 we show the rate of spatial convergence of the L error norms for different density contrasts, simultaneously comparing the behav-ior of the different methods STD and
MSUB . We maintain La = 3000 for alldensity contrasts, spatial resolutions and methods tested. We observe roughlysecond-order spatial convergence when it comes to equal densities across theinterface, with
STD and
MSUB displaying nearly identical behavior, although
MSUB performs slightly better with lower ( L ) errors for all resolutions. Whenit comes to ρ l /ρ g = 10 , we observe a saturation in the initial second-order con-37ergence rate irrespective of whichever method is used, however MSUB againhas lower errors. Finally, both methods demonstrate roughly second-order con-vergence when it comes to the stringent air-water configuration, with
MSUB performing better with lower errors at the moderate resolutions. Not surpris-ingly, the largest errors arise for the air-water configuration errors irrespectiveof the method.
5. Falling Raindrop
We turn our attention towards the issues plaguing standard non-consistentnumerical methods when dealing with realistic flows involving large density con-trasts (e.g air and water systems), primarily in the form of droplets and bubbles.A flow configuration that combines the difficulties of large density contrasts inconjunction with complexities involved due to the coupling of inertia, capillary,viscous stresses is that of a water droplet falling in air subject to gravity.
The problem is characterised by a combination of Reynolds, Weber and Bondnumbers, the definitions of which are as follows :We = ρ g U Dσ ,
Re = ρ g U Dµ g , Bo = ( ρ l − ρ g ) gD σ . (53)The density and viscosity ratios are corresponding to that of air-water sys-tems at 20 degree Celcius. An equivalent characterization could be done usingthe Morton number (cid:16) Mo = gµ g ( ρ l − ρ g ) ρ g σ (cid:17) instead of the Bond number. The sub-scripts l and g represent liquid and gas phases respectively. In our particularnumerical setup, We (cid:39) .
2, Re (cid:39) (cid:39) .
2, thus corresponding tothat of a 3 mm diameter raindrop (a relatively large one) falling in the air at anapproximate terminal velocity of 8 m/s (interpolated from empirical data, referto Gunn and Kinzer [18]). This choice of length scale of the drop is motivatedby the paradigmatic value of a near-spherical raindrop simulation, and by thefact that the corresponding Weber number ( ∼
3) is the same as in a similarair-water setup corresponding to suddenly-accelerated-droplet (or “secondaryatomization” ) simulations in the studies [47, 48]. For such a low Weber num-ber, the capillary forces should keep the interface shape more or less intact. Theparameters in the problem setup are given in Table 1, and a schematic diagramgiven in Fig. 20. The droplet is initially placed at the center of a cubic domain(3D), where the length of the side is 4 times the diameter of the drop.In order to properly reproduce and analyse the dynamics of a relativelylarge drop (high Reynolds flow) such as in our case, the numerical method hasto accurately resolve the thin boundary layers, the interaction of such layerswith the capillary forces, and finally the non-linear feedback of the complex 3Dvortical structures in the wake. In our experience, even for simulations with 64points across the droplet diameter, the resulting boundary layer has only by 3-438 g ρ l µ g µ l σ D g (cid:0) kg/m (cid:1) (cid:0) kg/m (cid:1) ( P a s ) (
P a s ) (
N/m ) ( m ) ( m/s )1.2 0 . × . × − . × − . × − . Table 1: Parameter values used in the simulation of a falling water droplet in air. cells across it. Arguably, the most natural type of computational setup wouldinvolve using a large domain filled with air at rest, with zero inflow velocityand to let the droplet fall from the top of the domain. Such an undertakingwas attempted by Dodd and Ferrante [13], in which they managed to delineatethe different regimes concerning the behavior of the wake behind the droplet,although at relatively low Reynolds numbers corresponding to smaller drops(the maximum Reynolds tested was (cid:39) (cid:39)
Figure 20: A 2D schematic of the falling raindrop computational setup. A droplet of diameter D is placed at the center of a cubic domain of side L and L/D = 4. The liquid properties( ρ l , µ l ) correspond to that of water, and the gas properties ( ρ g , µ g ) correspond to that ofair. We apply a uniform inflow velocity condition with U and an outflow velocity conditionat the top which corresponds to zero normal gradient. Boundary conditions on the side wallscorrespond to those of impenetrable free slip (zero shear stress). We use a significantly smaller domain (
L/D = 4 where L is the domainsize), with a constant value of inflow velocity (close to 8 m/s), implying thatthe drop will exit the domain after some time. This setup proves to be quiteconvenient and computationally efficient for relatively short-time investigations.Therefore, our objective behind the demonstration of this particular case is not
39o develop a high fidelity model of a raindrop, but instead carry out a stringentevaluation of the robustness of our present method compared to the standardnon-consistent version.
Numerical simulations using a fixed inflow velocity with U = 8 m/s werecarried out for very short times (a few milliseconds) at moderate resolutioncorresponding to D/ ∆ x = { , , , } , where D is the diameter and ∆ x the grid size. Simulations carried out using the standard non-consistent version( STD ) result in catastrophic deformations of the droplet as catalogued in Fig.21, which we describe as “fictitious” or “artificial” atomization. They displaymarked peaks or spikes in kinetic energy as a function of time, associated withmassively deformed interface shapes (see figure 22).
CIAMWY ENO WENO QUICK Superbee Verstappen
Figure 21: A catalogue of the different manifestations of “artificial” atomization encounteredwhile using the standard non-consistent method (
STD ). The red contours indicate the iso-surface of the volume fraction field corresponding to a value of 0 .
5, which is a good proxy forthe exact interfacial position. The resolution of the droplet in all the cases is D/ ∆ x = 16.The symbol “WY” in the legend corresponds to those run using the Weymouth-Yue advec-tion scheme, and “CIAM” corresponds to the a Lagrangian explicit flux advection scheme, asdetailed in [2, 3, 17]. The implementation of the non-linear flux limiters i.e WENO, ENO,QUICK, Superbee, Verstappen are identical to that of well established methods in the contextof hyperbolic conservation laws, the reader can refer to [3, 26, 43] for more details. The timestamps are indicative of the moments at which the interface appears most deformed, and donot necessarily refer to moments at which the solver crashes. We can clearly observe thatthe un-physical fragmentation of the raindrop is symptomatic of the non-consistent method( STD ), systematically across all combinations of flux limiters and advection schemes.
40 simplified model based on inviscid steady-state flow presented in our pre-vious study [3] explains the origins of this numerical instability. In essence,simulations carried out with the standard method (for moderate resolutions)lead to mixing of the gas velocity with the liquid density at mixed cells nearthe hyperpolic stagnation point at the droplet front. This results in extremelylarge (unphysical) pressure gradients which can only be balanced by the surfacetension force by invoking sufficiently large curvatures. These pressure “spikes”(see [3]) across the interface eventually lead to its rapid destabilization andconcomitant breakup.
Figure 22: Comparison of the temporal evolution of droplet kinetic energy between the non-consistent (
STD ) method and the present (
MSUB ) method. The kinetic energy is computedaccording to (55), and the droplet resolutions are D/ ∆ x = 8 for both sets of simulations. Inthe case of the non-consistent method, the kinetic energy of the drop climbs climbs rapidly overseveral orders of magnitude before the drop finally breaks up into many smaller fragments, Incontrast, kinetic energy of the drop evolves in a smooth manner using our present method. The most common and brute force approach that one can apply in order tosuppress or circumvent such numerical instabilities is by using a combination ofextremely refined meshes coupled with large domains [13]. In our study, we haveobserved that even a resolution of D/ ∆ x = 64 is not sufficient to suppress the in-stabilities in case of the standard method, although some combinations of advec-tion schemes (CIAM, Weymouth-Yue) and flux limiters (WENO,ENO,Superbeeetc) seem to be more stable than others. We observe that the use of our presentmethod enables us to stabilize the aforementioned numerical instabilities (com-pare Fig. 23 and Fig. 21), that too in a systematic manner that spans a widerange of flux limiters (WENO, ENO, Superbee, QUICK, Verstappen) and CFLnumbers. The kinetic energy of the droplet (using MSUB ) evolves in a rela-tively smooth manner, without the presence of sudden spikes and falls which41re emblematic of the non-consistent method (refer to Fig. 22). Such abruptchanges in kinetic energy of the droplet have been found to be associated withinstants when the droplet undergoes “artificial” atomization or breakup, hence-forth resulting in catastrophic loss of stability for the non-consistent (
STD )method. (i) (ii)(iii) (iv)
T = 3ms T = 5msT = 7ms T = 9ms
Figure 23: Numerical simulations using our present method (
MSUB ) of a 3 mm raindropfalling in air under gravitational acceleration. The flow is along the positive X direction witha constant inflow velocity of 8m/s, gravity is along the opposite direction, and with “QUICK”as the flux limiter. The blue contour indicates the isosurface of the volume fraction fieldcorresponding to a value of 0 .
5, and the instantaneous streamlines are colored according tothe velocity magnitude. The droplet has a relatively poor numerical resolution of D/ ∆ x = 16. Within the short time scales under consideration, the droplet shape undergoes oscillations(corresponding to the 2 nd Lamb mode [24]) along the flow direction, with the structure of thewake becoming more chaotic with time. This demonstrates that our present method (
MSUB )can easily suppress the numerical instabilties and massive interfacial deformations that arerampant in the case of the standard (
STD ) method. .2.1. Convergence Study For the next set of simulations, we use a fixed inflow velocity setup but withsmaller initial velocity. We systematically vary the resolution from D/ ∆ x =8 , ,
32 and 64. Despite using our present method, simulations at D/ ∆ x = 8are sometimes unstable, so we use a workaround and use a lower fixed inflowvelocity of U = 5 m/s. This simplification acts as a milder initial condition andenables us to observe the first phase of the (physical) acceleration towards thefinal statistical steady state. The simulations are carried out for 5 ms, in orderto minimize computational cost as well as to avert situations in which dropletapproaches too close to the domain boundaries. There are several time scalesin this problem, the important ones are • t a = D/U ≈ • t w = L/ [2( U t − U )] = 3 ms : The time taken by the drop to approachthe domain boundary, assuming terminal velocity. • t c (cid:39) . • t i = ( ρ l /ρ g ) · ( D/U t ) = 215 ms : The time scale to reach statisticallystationary state (terminal velocity), This estimate is obtained using asimple 1D model of the droplet dynamics (refer to [3]) using a square-velocity drag law. • t µ = Re · ( D/U t ) = 400ms : The time scale for the gas inertia to entrainvortical motion inside the drop (liquid side) due to the action of shear atthe interface.We show the results of our convergence study in figures 24 and 25. Thequantities of interest while examining the robustness of the present methodare the temporal evolution of the droplet kinetic energy (Fig. 24. (a) ) , andthe moments of inertia of the droplet along the three directions (Figures 24.(b),(c),(d)) . The inflow is along the X direction with gravity opposite to it.The moment of inertia is used as a descriptor of the “average” shape of thedrop, with their definitions along the different axes I m given by I m = (cid:90) D χx m d x , where , ≤ m ≤ , (54)where D is the computational domain and x m is the distance of the interfacerelative to the center of mass of the droplet. The droplet kinetic energy is definedrelative to the droplet center-of-mass, given by E k = (cid:104) ρ l · χ ( x, t ) || u ( x, t ) − u CM || (cid:105) , (55)43 Time ( × − s) D r o p l e t K i n e t i c E n e r g y ( K g . m / s ) × − (a) D/ ∆ x = 8 D/ ∆ x = 16 D/ ∆ x = 32 D/ ∆ x = 64 Time ( × − s) . . . . D r o p l e t M o m e n t o f I n e r t i a ( I xx ) × − (b) D/ ∆ x = 8 D/ ∆ x = 16 D/ ∆ x = 32 D/ ∆ x = 64 Time ( × − s) . . . . . . D r o p l e t M o m e n t o f I n e r t i a ( I yy ) × − (c) D/ ∆ x = 8 D/ ∆ x = 16 D/ ∆ x = 32 D/ ∆ x = 64 Time ( × − s) . . . . . . D r o p l e t M o m e n t o f I n e r t i a ( I zz ) × − (d) D/ ∆ x = 8 D/ ∆ x = 16 D/ ∆ x = 32 D/ ∆ x = 64 Figure 24: Temporal evolution of quantities of interest to evaluate the performance of ourconsistent scheme based on the sub-grid (
MSUB ) method, for different spatial resolutions.(a) Kinetic energy relative to the droplet center-of-mass as defined in (55). (b) Moment ofinertia of the droplet along the flow (X) direction. (c) and (d) Moment of inertia of the dropletalong the directions perpendicular to flow (Y,Z), evolution of I yy seems to be more or lessidentical to I zz . where (cid:104) . (cid:105) is the spatial averaging operator over the entire domain and u CM is the droplet center of mass. We observe a systematic drop in the amount ofthe droplet kinetic energy as we increase resolution, with the most probableexplanation being that of the suppression of spurious interfacial jitter which isrampant at lower resolutions.There is also a component of the kinetic energy of the droplet associatedwith the internal coherent vortical structures generated due to the interactionwith aerodynamic shear at the interface, evidenced by the non-zero value of thekinetic energy even for the most highly resolved droplets. Finally, the momentsof inertia of the droplet evolve in a smooth manner, and seem to converge as wego from D/ ∆ x = 32 to D/ ∆ x = 64, for each of the 3 directions.In fig. 25, we plot the the velocity of the center of mass of the dropletas a function of time, and its behavior for increasing droplet resolutions. Thetemporal variation in the droplet velocity is fitted to a linear polynomial inorder to evaluate the droplet acceleration by means of a standard least-squaresapproach. The acceleration corresponding to the finest droplet resolution is44 . . . . . . . . Time ( × − s) − − D r o p l e t V e l o c i t y ( × − m / s ) MSUB : D/ ∆ x = 8 MSUB : D/ ∆ x = 16 MSUB : D/ ∆ x = 32 MSUB : D/ ∆ x = 64 MSHIFT : D/ ∆ x = 642 Resolution ( D/ ∆ x ) − − − − − A cce l e r a t i o n ( m / s ) MSUBMSHIFT
Figure 25: Comparison of droplet (center of mass) velocity as a function of time, for differentdroplet resolutions. The simulations were carried out with a consistent scheme based onthe present (
MSUB ) method, using the WY advection scheme with the QUICK limiter.For comparison with our previous study [3], we have added the results obtained with the shifted fractions method (
MSHIFT , WY and QUICK), corresponding to D/ ∆ x = 64. Inset: Droplet acceleration estimates as a function of its resolution, computed using the best linearfit over the interval [1ms , d U d t (cid:39) . ± . m/s , where the uncertainty is estimated from the differencebetween D/ ∆ x = 2 and the D/ ∆ x = 2 . The acceleration estimate obtainedfrom our previous study [3] which uses the shifted fractions ( MSHIFT ) methodwith an identical setup, is d U d t (cid:39) . ± . m/s . The discrepancy between theacceleration estimates are clearly visible in the inset of fig. 25, where the shiftedfractions method displays some oscillatory behaviour with increasing resolution.The present method ( MSUB ) on the other hand displays a monotonic trendtowards a converged estimate, and also has lower fitting errors compared to
MSHIFT (inset fig. 25), indicating that the
MSHIFT method is “noisier”than our present method.Fig. 26 shows the convergence rate of the methods. A relative error estimate45 Resolution ( D/ ∆ x ) − − R e l a t i v e E rr o r E s t i m a t e MSUBMSHIFTSecond Order
Figure 26: Variation in the error estimates of droplet acceleration as a function of dropletresolution. The
MSHIFT method displays a better convergence rate, but that is due tothe wildly oscillating behaviour when it comes to the acceleration estimate as a functionof resolution. The
MSUB method on the other hand, displays a monotonic decrease inthe acceleration estimate with resolution, although with a saturating convergence rate. Thesecond order convergence rate is plotted as a reference. at resolution level n (i.e. D/ ∆ x = 2 n ) is defined as the difference between theacceleration estimates (inset fig. 25) at the n and n + 1 levels, normalized bythe value at level n . The underlying cause of the difference between MSUB and
MSHIFT methods is most likely due to the differing levels of dissipationbetween the methods. The shifted fractions method is more dissipative, whichmight lead to faster (and greater) momentum diffusion. This might explainthe differences we have observed in the (not shown) temporal development ofthe boundary layers around the raindrop, between the
MSUB and
MSHIFT methods. A more detailed and quantitative anylysis is warranted to verify theabove assertions, but it is beyong the scope of our present study. To sum it up,the results obtained in the case of the falling droplet strongly suggest that ourpresent numerical method (
MSUB ) can be used to get relatively good estimatesof the underlying flow features of the drop, without observing any un-physicalevolution due to the discretization errors at low to moderate resolutions.
6. Conclusions & Perspectives
We have presented a numerical method designed to simulate complex liquid-gas interfacial flows involving significant contrasts in density. The conservativeformulation of governing equations are solved with momenutum and density asthe primary variables, on staggered Cartesian grids. Rudman’s [41] strategy46f mass advection on a doubly refined grid is adopted using geometric flux re-constructions so as to maintain consistent mass-momentum advection on thestaggered grid. The Weymouth and Yue algorithm [46] is extended to momen-tum transport, in order to ensure discrete consistency between the conservativedirection-split (3D) transport of volume, mass and momentum.The performance of the method is assessed using standard benchmark casessuch as the behavior of spurious currents in static and moving droplets, andthe propagation of damped capillary waves. All of the test cases tests wereconducted for density contrasts spanning 3 orders of magnitude, including arealistic air-water configuration for capillary waves. Our method demonstratesincreased accuracy compared to the standard non-consistent version. A second-order error convergence rate is observed for the static droplet and the capillarywave, with the rate reverting to first-order in the case of droplet advection.Finally, the robustness and stability of our present method was demonstratedusing the case of a raindrop falling in air under gravitational acceleration. A keyfeature of the method is its ability to avoid un-physical interfacial deformations(artificial atomization of the raindrop), which are rampant in the case of thenon-consistent method. The method delivers good estimates of the flow fea-tures especially at low resolution, as evidenced by the smooth (and convergent)evolution of the droplet kinetic energy and the moment of inertia. A short timeinvestigation of the falling drop was able to reproduce the physically consistentinitial acceleration phase (towards the final statistically stationary state) of theraindrop.In terms of future perspectives, the present method can be further adaptedto derive more accurate curvature estimates by using the volume fraction infor-mation available from the twice finer grid, and we are currently working towardsfinding the optimal strategy for the such curvature restrictions. The methodis also going to be used to investigate several flow configurations of scientificinterest such as atomization of shear-layers and Kelvin-Helmholtz instabilties,secondary atomization of drops etc.
Acknowledgements
This work has been supported by the grant SU-17-R-PER-26-MULTIBRANCHfrom Sorbonne Universit´e and the ERC Advanced Grant TRUFLOW. This workwas granted access to the HPC resources of TGCC- CURIE, TGCC-IRENE andCINES-Occigen under the allocations A0032B07760 and 2020225418, made byGENCI and PRACE respectively. We would like to thank G. Tryggvason, R.Scardovelli, Dr. Y. Ling, Dr. W. Aniszewski, Dr. S. Dabiri, Dr. J. Lu andDr. P. Yecko for their contribution to the development of the code “PARISSimulator”. Finally, the simulation data are visualized by the software VisIt,developed by the Lawrence Livermore National Laboratory.47 ppendix A :
Brief review of mass-momentum consistent methods
A bird’s eye view of the numerous features employed by the methods inexisting literature, we refer the reader to tables 2 and 3, which respectivelyprovide systematic overviews of the VOF based and level set based approaches.48 a b l e : Su mm a r y o f V o l u m e - o f - F l u i d b a s e d m a ss - m o m e n t u m c o n s i s t e n t a d v ec t i o n m e t h o d s f o r i n c o m p r e ss i b l e fl o w s o f i mm i s c i b l e flu i d s . R ud m a n ( I J N M F ) B u ss m a nn e t a l. ( A S M E ) L e C h e n a d ec & P i t s c h ( J C P ) O w k e s & D e s j a r d i n s ( J C P ) P a t e l & N a t a r a j a n ( J C P ) A rr u f a t e t a l. ( C A F ) B a s i c C o nfi g u r a t i o n D C a r t e s i a n , s t agg e r e d D h e x a h e d r a l un s t r u c t u r e d , c o ll o c a t e d D C a r t e s i a n , s t agg e r e d D C a r t e s i a n , s t agg e r e d D p o l y go n a l un s t r u c t u r e d , s t agg e r e d D C a r t e s i a n s t agg e r e d I n t e r f a c e R e p r e s e n t a t i o n V O F , P i ece w i s e L i n e a r V O F , P i ece w i s e L i n e a r V O F , P i ece w i s e L i n e a r V O F , P i ece w i s e L i n e a r V O F , n o n e V O F , P i ece w i s e L i n e a r F l u x C o m pu t a t i o n s p li t , a l g e b r a i c F l u x C o rr ec t e d T r a n s p o r t , E u l e r i a nun s p li t , g e o m e t r i c , E u l e r i a nun s p li t , g e o m e t r i c , s e m i - L ag r a n g i a nun s p li t , g e o m e t r i c , s e m i - L ag r a n g i a n a l g e b r a i c C ub i c U p w i nd , E u l e r i a n s p li t , g e o m e t r i c , E u l e r i a n Su r f a c e T e n s i o n C o n t i nuu m Su r f a ce - F o r ce n o t s p ec i fi e d G h o s t - F l u i d M e t h o d , w e ll - b a l a n ce d C o n t i nuu m Su r f a ce - F o r ce , w e ll - b a l a n ce d C o n t i nuu m Su r f a ce - F o r ce , w e ll - b a l a n ce d C o n t i nuu m Su r f a ce - F o r ce , w e ll - b a l a n ce d C u r v a t u r e E s t i m a t i o n V O F c o n v o l u t i o n ( s m oo t h e d ) n o t s p ec i fi e d V O F c o n v o l u t i o n ( s m oo t h e d ) h y b r i d m e s h - d ec o up l e d h e i g h t f un c t i o n s V O F c o n v o l u t i o n ( s m oo t h e d ) h y b r i dh e i g h t f un c t i o n s a nd c u r v e fi tt i n g V i s c o u s S t r e ss e s e x p li c i t , h a r m o n i c a v e r ag i n g n o t s p ec i fi e d e x p li c i t , h a r m o n i c a v e r ag i n g n o t s p ec i fi e d e x p li c i t , h a r m o n i c a v e r ag i n g e x p li c i t , a r i t h m e t i c a v e r ag i n g P r e ss u r e - P o i ss o n S o l v e r p r ec o nd i t i o n e d m u l t i g r i d , G a u ss - S e i d e l n o t s p ec i fi e dn o t s p ec i fi e dn o t s p ec i fi e dp r ec o nd i t i o n e d G M R E S , K r y l o v s ub s p a ce p r ec o nd i t i o n e d m u l t i g r i d , G a u ss - S e i d e l a b l e : Su mm a r y o f L e v e l S e t a nd C L S V O F b a s e d m a ss - m o m e n t u m c o n s i s t e n t a d v ec t i o n m e t h o d s f o r i n c o m p r e ss i b l e fl o w s o f i mm i s c i b l e flu i d s R a e ss i & P i t s c h ( C A F ) G h o d s & H e rr m a nn ( P h y s i c a S c r i p t a 2013 ) V a ud o r e t a l. ( C A F ) N a n g i a e t a l. ( J C P ) Z u z i o e t a l. ( J C P ) B a s i c C o nfi g u r a t i o n D C a r t e s i a n , s t agg e r e d D h e x a h e d r a l un s t r u c t u r e d , c o ll o c a t e d D C a r t e s i a n , s t agg e r e d D C a r t e s i a n s t agg e r e d D C a r t e s i a n s t agg e r e d I n t e r f a c e R e p r e s e n t a t i o n L e v e l S e t L e v e l S e t C o up l e d L e v e l S e t - V O F L e v e l S e t C o up l e d L e v e l S e t - V O F F l u x C o m pu t a t i o n L e v e l S e t d e r i v e d , s e m i - L ag r a n g i a n L e v e l S e t d e r i v e d , E u l e r i a n s p li t , g e o m e t r i c , E u l e r i a n a l g e b r a i c C ub i c U p w i nd , E u l e r i a n s p li t , g e o m e t r i c , E u l e r i a n Su r f a c e T e n s i o n G h o s t - F l u i d M e t h o d C o n t i nuu m Su r f a ce - F o r ce G h o s t - F l u i d M e t h o d , w e ll - b a l a n ce d C o n t i nuu m Su r f a ce - F o r ce , w e ll - b a l a n ce d G h o s t - F l u i d M e t h o d C u r v a t u r e E s t i m a t i o n L e v e l S e t b a s e d L e v e l S e t b a s e d L e v e l S e t b a s e d L e v e l S e t b a s e d L e v e l S e t b a s e d V i s c o u s S t r e ss e s i m p li c i t , h a r m o n i c a v e r ag i n g e x p li c i t , a r i t h m e t i c a v e r ag i n g s e m i - i m p li c i t , h a r m o n i c a v e r ag i n g e x p li c i t , a r i t h m e t i c a v e r ag i n g e x p li c i t , a r i t h m e t i c a v e r ag i n g P r e ss u r e - P o i ss o n S o l v e r p r ec o nd i t i o n e d m u l t i g r i d , K r y l o v s ub s p a ce n o t s p ec i fi e dp r ec o nd i t i o n e d m u l t i g r i d , C o n j u ga t e G r a d i e n t p r ec o nd i t i o n e d G M R E S , K r y l o v s ub s p a ce p r ec o nd i t i o n e d m u l t i g r i d , C o n j u ga t e G r a d i e n t eferences [1] Abadie, T., Aubin, J., Legendre, D., 2015. On the combined effects of sur-face tension force calculation and interface advection on spurious currentswithin volume of fluid and level set frameworks. Journal of ComputationalPhysics 297, 611–636.[2] Aniszewski, W., Arrufat, T., Crialesi-Esposito, M., Dabiri, S., Fuster, D.,Ling, Y., Lu, J., Malan, L., Pal, S., Scardovelli, R., et al., 2019. Parallel,robust, interface simulator (paris) .[3] Arrufat, T., Crialesi-Esposito, M., Fuster, D., Ling, Y., Malan, L., Pal, S.,Scardovelli, R., Tryggvason, G., Zaleski, S., 2020. A mass-momentum con-sistent, volume-of-fluid method for incompressible flow on staggered grids.Computers & Fluids , 104785.[4] Boeck, T., Li, J., L´opez-Pag´es, E., Yecko, P., Zaleski, S., 2007. Ligamentformation in sheared liquid–gas layers. Theoretical and ComputationalFluid Dynamics 21, 59–76.[5] Bornia, G., Cervone, A., Manservisi, S., Scardovelli, R., Zaleski, S., 2011.On the properties and limitations of the height function method in two-dimensional cartesian geometry. Journal of Computational Physics 230,851–862.[6] Bourouiba, L., Dehandschoewercker, E., Bush, J.W., 2014. Violent expi-ratory events: on coughing and sneezing. Journal of Fluid Mechanics 745,537–563.[7] Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method formodeling surface tension. Journal of computational physics 100, 335–354.[8] Briggs, W.L., 1987. A multigrid tutorial, siam. Carver, MB (1984)“Numer-ical Computation of Phase Separation in Two Fluid Flow,” ASME Journalof Fluids Engineering 106, 147–153.[9] Bussmann, M., Kothe, D.B., Sicilian, J.M., 2002. Modeling high densityratio incompressible interfacial flows, in: ASME 2002 Joint US-EuropeanFluids Engineering Division Conference, American Society of MechanicalEngineers. pp. 707–713.[10] Chorin, A.J., 1969. On the convergence of discrete approximations to thenavier-stokes equations. Mathematics of computation 23, 341–353.[11] Deckwer, W.D., 1980. On the mechanism of heat transfer in bubble columnreactors. Chemical Engineering Science 35, 1341–1346.[12] Deike, L., Lenain, L., Melville, W.K., 2017. Air entrainment by breakingwaves. Geophysical Research Letters 44, 3779–3787.5113] Dodd, M.S., Ferrante, A., 2014. A fast pressure-correction method forincompressible two-fluid flows. Journal of Computational Physics 273, 416–434.[14] Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M.,Williams, M.W., 2006. A balanced-force algorithm for continuous andsharp interfacial surface tension models within a volume tracking frame-work. Journal of Computational Physics 213, 141–173.[15] Ghods, S., Herrmann, M., 2013. A consistent rescaled momentum transportmethod for simulating large density ratio incompressible multiphase flowsusing level set methods. Physica Scripta 2013, 014050.[16] Gilet, T., Bourouiba, L., 2015. Fluid fragmentation shapes rain-inducedfoliar disease transmission. Journal of the Royal Society Interface 12,20141092.[17] Gueyffier, D., Li, J., Nadim, A., Scardovelli, R., Zaleski, S., 1999. Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows. Journal of Computational physics 152, 423–456.[18] Gunn, R., Kinzer, G.D., 1949. The terminal velocity of fall for waterdroplets in stagnant air. Journal of Meteorology 6, 243–248.[19] Harvie, D.J., Davidson, M., Rudman, M., 2006. An analysis of parasiticcurrent generation in volume of fluid simulations. Applied mathematicalmodelling 30, 1056–1066.[20] Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (vof) method for thedynamics of free boundaries. Journal of computational physics 39, 201–225.[21] Johansen, S., Boysan, F., 1988. Fluid dynamics in bubble stirred ladles:Part ii. mathematical modeling. Metallurgical Transactions B 19, 755–764.[22] Kantarci, N., Borak, F., Ulgen, K.O., 2005. Bubble column reactors. Pro-cess biochemistry 40, 2263–2283.[23] Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S., Zanetti, G., 1994.Modelling merging and fragmentation in multiphase flows with surfer. Jour-nal of Computational Physics 113, 134–147.[24] Lamb, H., 1993. Hydrodynamics. Cambridge university press.[25] Le Chenadec, V., Pitsch, H., 2013. A monotonicity preserving conservativesharp interface flow solver for high density ratio two-phase flows. Journalof Computational Physics 249, 185–203.[26] LeVeque, R.J., 1996. High-resolution conservative algorithms for advectionin incompressible flow. SIAM Journal on Numerical Analysis 33, 627–665.5227] Mirjalili, S., Jain, S.S., Dodd, M., 2017. Interface-capturing methods fortwo-phase flows: An overview and recent developments. Center for Turbu-lence Research Annual Research Briefs 2017, 117–135.[28] Nangia, N., Griffith, B.E., Patankar, N.A., Bhalla, A.P.S., 2019. A robustincompressible navier-stokes solver for high density ratio multiphase flows.Journal of Computational Physics 390, 548–594.[29] Osher, S., Sethian, J.A., 1988. Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations. Jour-nal of computational physics 79, 12–49.[30] Owkes, M., Desjardins, O., 2015. A mesh-decoupled height function methodfor computing interface curvature. Journal of Computational Physics 281,285–300.[31] Owkes, M., Desjardins, O., 2017. A mass and momentum conserving un-split semi-lagrangian framework for simulating multiphase flows. Journalof Computational Physics 332, 21–46.[32] Patel, J.K., Natarajan, G., 2017. A novel consistent and well-balanced al-gorithm for simulations of multiphase flows on unstructured grids. Journalof computational physics 350, 207–236.[33] Popinet, S., . Basilisk, a free-software program for the solution of par-tial differential equations on adaptive cartesian meshes (2018). URLhttp://basilisk. fr .[34] Popinet, S., 2009. An accurate adaptive solver for surface-tension-driveninterfacial flows. Journal of Computational Physics 228, 5838–5866.[35] Popinet, S., 2018. Numerical models of surface tension. Annual Review ofFluid Mechanics 50, 49–75.[36] Popinet, S., Zaleski, S., 1999. A front-tracking algorithm for accurate repre-sentation of surface tension. International Journal for Numerical Methodsin Fluids 30, 775–793.[37] Prosperetti, A., 1980. Free oscillations of drops and bubbles: the initial-value problem. Journal of Fluid Mechanics 100, 333–347.[38] Prosperetti, A., 1981. Motion of two superposed viscous fluids. The Physicsof Fluids 24, 1217–1223.[39] Raessi, M., Pitsch, H., 2012. Consistent mass and momentum transport forsimulating incompressible interfacial flows with large density ratios usingthe level set method. Computers & Fluids 63, 70–81.[40] Rayleigh, L., 1879. On the capillary phenomena of jets. Proc. R. Soc.London 29, 71–97. 5341] Rudman, M., 1998. A volume-tracking method for incompressible multi-fluid flows with large density variations. International Journal for numericalmethods in fluids 28, 357–378.[42] Seinfeld, J.H., Pandis, S.N., 1998. From air pollution to climate change.Atmospheric Chemistry and Physics , 1326.[43] Sweby, P.K., 1984. High resolution schemes using flux limiters for hyper-bolic conservation laws. SIAM journal on numerical analysis 21, 995–1011.[44] Tryggvason, G., Scardovelli, R., Zaleski, S., 2011. Direct numerical simu-lations of gas–liquid multiphase flows. Cambridge University Press.[45] Vaudor, G., M´enard, T., Aniszewski, W., Doring, M., Berlemont, A., 2017.A consistent mass and momentum flux computation method for two phaseflows. application to atomization process. Computers & Fluids 152, 204–216.[46] Weymouth, G.D., Yue, D.K.P., 2010. Conservative volume-of-fluid methodfor free-surface simulations on cartesian-grids. Journal of ComputationalPhysics 229, 2853–2865.[47] Xiao, F., 2012. Large eddy simulation of liquid jet primary breakup. Ph.D.thesis. ©©