Prof Stéphane Zaleski

Samta Kothari

Gaurav Prabhudesai

Anastasia Loukrezi

Christina Agapopoulou

François Chantriaux

Georg Lopez

The completion of the present body of

List of Figures

acceleration. The flow is along the positive X direction, with gravity along the opposite direction. The red contour indicates the isosurface of the volume fraction field corresponding to a value of 0.5, whereas the black contours surrounding the drop represent isosurfaces of the magnitude of vorticity.

The droplet has a relatively poor numerical resolution, corresponding to D/h 16, where D is the drop diameter and h is the grid size. The top row illustrates the phenomenon of "artificial atomization" that occurs while using the standard version of our numerical method, which does not maintain consistency between the discrete transport of mass and momentum. In contrast, the implementation of consistent mass and momentum transport enables us to suppress the massive interfacial deformations of the drop, as evidenced by the bottom row of figures depicting the same flow configuration. The development of methods concerning the consistent transport of mass and momentum is one of the central themes of this study. . . A D schemetic of the Eulerian (geometric) flux calculation using the Weymouth-Yue [ ] scheme for the advection substep along the horizontal direction, with the interface reconstructed using the volume fraction field at the start of the substep. The colour fraction of 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 -1, j and the red one corresponding to that exiting i, j into i + 1, j. The geometric flux calculations are made on the basis of the interfacial positions at the start of the substep, and the face centered velocities of the cell in question. . . . . . . . . . . A D schematic of the staggered spatial configuration of the pressure and velocity variables. The pressure p i, j is based on the center of its control volume (light colour area); the horizontal velocity component u 1;i+1/2, j is defined in the middle of the right edge of the pressure control volume and centered on its control volume (dash-dotted area); the vertical velocity component u 2;i, j+1/2 is defined in the middle of the top edge and centered on its corresponding control volume (dashed area). The volume fraction is defined at the same location as that of the pressure, resulting in pressure and density being centered on identical control volumes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The reference control volume Ω for the advected velocity component φ u q is shown. A horizontal advection is here considered and both the advecting velocity u f and the advected velocity require an interpolation for their value on the left face f 1-: (a) the value ūq φ 0 (full circle) is interpolated from the values φ u q on the nodes (open circles); (b) a more sophisticated interpolation predicts the value φ( x) where x is at the center of the "donating" region Ω D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A D schematic of the reconstruction procedure used in the shifted fractions method. In this case, advection of the horizontal component of momentum is considered. The half-fractions from the centered control volumes Ω i, j and Ω i+1, j (dash-dotted area) are added together to obtain the volume fraction in the shifted (staggered) volume Ω i+1/2, j . Similarly, the half-fractions from Ω i-1, j and Ω i, j (dotted area) are combined to obtain the volume fraction Ω i-1/2, j . Once this horizontally shifted volume fraction field is reconstructed, the corresponding density, momentum and "advected" velocity fields can be computed in order to carry out consistent transport of mass and momentum on this staggered configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A bird's eye view of the shifted fractions method, highlighting the operations performed at each time step. This method achieves consistency in the discrete transport of mass and momentum on the staggered control volumes. In the interest of brevity, we present a D case for the density on the grid i, j and horizontal velocity u 1 on the staggered grid i + 1/2, j. The evolution of the velocity component u 2 on the staggered grid i, j + 1/2 is similar. The initial variables ρ n , u n 1 , u n 2 are inside the ellipses. The interpolated "advecting" components u 1 and u 2 have superscript n. The shifted density ρ n,0 is constructed with the shifted fractions of C to initialize the momentum component (ρu 1 ) n,0 . The first split advection is along the x direction to variables with superscript n, 1, the second one is along the y direction to variables with superscript n, 2.

The updated density is ρ n,2 ρ n+1 while the horizontal velocity u n,2 . A D schematic of the arrangement of primary variables on the coarse and sub-grid levels. The control volume is decomposed into the two figures, the one on the left corresponding to the coarse grid variables, while the one on the right corresponding to the sub-grid. Superposition of these two figures, one on top of each other, represent the spatial relationships between the two sets of variables. The velocities on the sub-grid level are simple first or zeroth order interpolations of the coarse grid velocities. . . . . . . . . . . . . . . . . . . . . . A D illustration of the restriction operations involved between the different grid resolutions, where the sub-grid grid boundaries are denoted by the dashed lines, and that of the overlapping coarse grid by solid lines. Restriction operations in the form of simple volume integrals are performed on the subgrid volume fraction in order to compute the staggered and centered volume fraction fields at the coarse grid level. Restriction operators in the form of surface integrals are also performed in order to compute the volume fluxes at the boundaries of the coarse staggered control volumes. The red shaded areas correspond to the volume fractions, whereas the green shaded area corresponds to that of volume fluxes. The interface segments as depicted at the coarse grid level are just for purpose of illustration, and are not 'reconstructed'.

Unlike the shifted fractions method, the interfacial reconstructions on the staggered volumes are not necessary, due to the fact that sub-grid volume fluxes are simply combined to compute the staggered fluxes. . . . . . . . . . . . A bird's eye view of the sub-grid method, highlighting the operations performed at each time step. and additionally ensures conservation of momentum. For sake of clarity, present a D case for the density on the grid i, j and horizontal velocity ũ on the staggered grid i + 1/2, j. The variables in red are all defined on the coarse level, and those in green are at the sub-grid level. 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 convery that the variables are centered on the staggered control volumes at the coarse grid level, whereas the 'bar' on top conveys that the variables are staggered with respect to the sub-grid. The starting variables are enclosed in elliptical borders, 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. . . . . . . . . . . time-scales for different density-ratios and Laplace numbers. The currents seem to initially decay quickly for all higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For combinations of large ρ l /ρ g and large La, the spurious currents seem to grow back to an order of magnitude (10 -8 ) which is quite far from that of machine precision (10 -14 ). . MSHIFT Decay of normalized spurious currents as a function of viscous dissipation time-scales for different density-ratios and Laplace numbers. The currents seem to initially decay quickly for all higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For combinations of large ρ l /ρ g and large La, the spurious currents seem to grow back to an order of magnitude (10 -8 ) which is quite far from that of machine precision (10 -14 ). No considerable improvement is observed with respect to STD. . . .

. MSUB Decay of normalized spurious currents as a function of viscous dissipation time-scales for different density-ratios and Laplace numbers. The currents seem to decay very quickly in the case of higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For all combinations of ρ l /ρ g and La numbers, the decayed spurious currents are not observed to grow back as in the cases of STD and MSHIFT, and hover around values close to machine precision (10 -14 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Second-order spatial convergence for the spurious current error norms corresponding to the most stringent parameter combination (ρ l /ρ g 1000 , La 12000) . Both of the norms (L ∞ and L 2 ) seem to demonstrate a roughly second order rate of spatial convergence with each of the methods tested. However, MSUB has a marginally lower L 2 error compared to both STD and MSHIFT for all resolutions tested. There is negligible difference observed in the shape errors between STD and MSHIFT in both of the norm definitions. . . Schematic of the droplet of dense fluid advected in a surrounding medium of lighter fluid. A 50 × 50 grid is employed to spatially discretize the domain, which is spatially periodic in the direction of droplet advection. . . . . . . . . . STD Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers. The currents seem to hover around 10 -3 , with a larger Laplace number corresponding to a higher error for all density-ratios. We 0.4 for all the cases presented. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MSHIFT Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers. There seems to be no appreciable difference from the evolution seen in the case of STD (fig. . ). The currents seem to hover around 10 -3 , with a larger Laplace number corresponding to a higher error for all density-ratios. We 0.4 for all the cases presented. . . . . . . . . . . . . . . . . . MSUB Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers. In terms of the errors observed in STD and MSHIFT, 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 -4 , with larger Laplace numbers corresponding to larger errors for all density-ratios. We 0.4 for all the cases presented. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . First-order (approximately) spatial convergence of the maximum of the spurious current error norms in the frame of reference of the moving droplet, for the most stringent parameter combination (ρ l /ρ g 1000 , La 12000, We 0.4). Methods STD and MSHIFT display similar convergence properties, whereas MSUB leads to significantly lower errors even though it doesn't quite follow the first-order convergence rate. . , with the width of the box corresponding to the size of the perturbed wavelength. . . . . STD Time evolution of the amplitude of the planar interface undergoing damped capillary oscillations, comparing the solution obtained by our numerical method with the closed-from Prosperetti solution. More or less good agreement with theory is observed for all the density-ratios tested. . . . . . .

. MSHIFT Time evolution of the amplitude of the planar interface undergoing damped capillary oscillations, comparing the solution obtained by our numerical method with the closed-from Prosperetti solution. Behavior is quite similar to STD, with good agreement with theory for all the density-ratios tested. . . . MSUB Time evolution of the amplitude of the planar interface undergoing damped capillary oscillations, comparing the solution obtained by our numerical method with the closed-from Prosperetti solution. Slightly better agreement with theory when comparing to STD and MSHIFT, for all density-ratios tested. . Comparison of spatial convergence for the case of ρ l /ρ g 1, La 3000, for our class of methods. There is no viscosity jump across the interface. All methods seem to demonstrate approximately second-order convergence. There seems to be no appreciable difference in the behavior of STD and MSHIFT, with MSUB displaying marginally lower errors compared to the others. . . . . . . . Comparison of spatial convergence for the case of ρ l /ρ g 10, La 3000, for our class of methods. Again, there is no viscosity jump across the interface. All methods seem to demonstrate approximately second-order convergence upto L/∆x 32, beyond which there is a slight saturation in the rate of convergence. Qualitatively, STD and MSHIFT demonstrate similar behavior, with MSUB delivering the lowest errors. In case of MSUB, the errors are lower by a factor of 2 compared to STD and MSHIFT for resolutions above L/∆x 16. . . . . . Comparison of spatial convergence for the Air-Water case corresponding to ρ l /ρ g 1000.0/1.2, µ l /µ g 1.003 • 10 -3 /1.8 • 10 -5 , La 3000, for our class of methods. All methods seem to demonstrate approximately secondorder convergence. No appreciable difference is observed between STD and MSHIFT, with MSUB delivering slightly lower errors although there is some saturation in the convergence rate at higher resolutions. . . . . . . . . . . . . . . A D schematic of the numerical setup for the falling raindrop. 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 of air. We apply a uniform inflow velocity condition with U 0 (t) and an outflow velocity condition at the top which corresponds to zero normal gradient. Boundary conditions on the side walls correspond to those of impenetrable free slip (no shear stress). . . . . . . . . . . . . . . . . . . The origin of the pressure peak in the front of the droplet. (a) The profile of the pressure on the axis a few timesteps after initialisation with the standard, nonconsistent method (red curve) and the consistent method based on the shifted fractions (MSHIFT) strategy (green curve). Much larger pressure gradients are present across the interface using the non-consistent method. (b) The pressure distribution immediately after the start of the simulation using the standard, non-consistent method. The pressure peak has not yet resulted in the formation of a dimple. In all figures the droplet resolution corresponds to D/h 16. The simulations are carried out with the CIAM advection method, in conjunction with the Superbee limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Comparison of the temporal evolution of droplet kinetic energy. The standard non-consistent method displays spikes in the kinetic energy that are approximately orders of magnitude larger than that of the consistent method, leading to rapid destabilization. (b) The pressure distribution immediately after the start of the simulation using the consistent method based on the shifted fractions (MSHIFT) strategy. In all figures the droplet resolution corresponds to D/h 16. The simulations are carried out with the CIAM advection method, in conjunction with the Superbee limiter. . . . . . . . . . . . . . . . . . . . . .

. Flow field around the 3mm droplet with D/h 16 immediately after the start of the simulation with the consistent method (MSHIFT), demonstrating the contours of the norm of the vorticity field (black lines) . The D cross-section in these figures corresponds to the mid-plane slice along the Z axis, where the inflow is along the X axis and gravity opposite to it. (a) The velocity magnitude.

As one can observe, the boundary layer is resolved by only -cells. (b) The velocity component in the Y direction, perpendicular to the flow. As the flow develops further, a marked separation of the boundary layers is observed with a more complex vortical region in the wake. The figures correspond to simulations carried out with the CIAM advection method, in conjunction with the Superbee limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A comparison of the temporal evolution between the standard non-consistent method (top row) and the consistent method based on the shifted fractions (MSHIFT) strategy (bottom row), using the combination of CIAM advection scheme coupled with the Superbee flux limiter. The flow is along the positive X direction, with gravity along the opposite direction. The red contour indicates the isosurface of the volume fraction field corresponding to a value of 0.5, whereas the black contours surrounding the drop represent isosurfaces of the magnitude of vorticity. The raindrop with the non-consistent method displays massive deformations leading to artificial breakup as a result of rapidly growing numerical instabilities. The droplet resolution for both methods is D/h 16 . The temporal evolution in case of the sub-grid based consistent method (MSUB) produces almost indistinguishable figures (bottom row), hence are not shown here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A juxtaposition of the different manifestations of 'artificial' atomization one encounters while using the standard method (STD) that does not involve consistency between the discrete transport of mass and momentum. The red contours indicate the isosurface of the volume fraction field corresponding to a value of 0.5, acting as a good proxy for the exact interfacial position.

One can clearly observe that the un-physical fragmentation of the raindrop is symptomatic of the non-consistent method, systematically across all combinations of flux limiters and advection schemes. The symbol "WY" in the legend corresponds to those run using the Weymouth-Yue advection scheme, and "CIAM" corresponds to the CIAM scheme. Brief descriptions of these advection schemes can be found in chapter . The implementation of the non-linear flux limiters i.e WENO, ENO, QUICK, Superbee, Verstappen correspond to that of well established methods developed to deal with hyperbolic conservation laws, for more details refer to the studies of Leveque [ ] and Sweby [ ] . Comparison of droplet velocity as a function of time, for different droplet resolutions. The simulations were carried out with a consistent scheme based on the shifted fractions (MSHIFT) method, using the WY advection scheme with the QUICK limiter. The droplet velocities correspond to that of their respective center of masses. Inset : Convergence of the droplet acceleration as a function of its resolution, computed using the best linear fit over the temporal variation of their respective velocities. The error bars signify the asymptotic standard error (least-squares) corresponding to the linear fits. . . . . . . . . . . Comparison of droplet velocity as a function of time, for different droplet resolutions. The simulations were carried out with a consistent scheme based on the sub-grid (MSUB) method, using the WY advection scheme with the QUICK limiter. The droplet velocities correspond to that of their respective center of masses. Inset : Convergence of the droplet acceleration as a function of its resolution, computed using the best linear fit over the temporal variation of their respective velocities. The error bars signify the asymptotic standard error (least-squares) corresponding to the linear fits. . . . . . . . . . . . . . . . . Relative change in the mass of the droplet as a function of time in plots (a), (b) and (c). The simulations are carried out with a resolution of D/h 16, for a total time of 0.1 milliseconds using the consistent scheme based on the shifted fractions (MSHIFT) method. The symbol "WY" in the legend corresponds to those run using the Weymouth-Yue advection scheme in combination with the QUICK limiter, and "CIAM" corresponds to the CIAM scheme in combination with the Superbee limiter. (a) Mass conservation properties of the scheme as a function of the Poisson solver tolerance is tested. For the WY-QUICK combination, using a stricter tolerance results in better mass conservation, whereas the conservation properties of the CIAM-Superbee combination seems to be independent of the tolerance. (b) Dependence of the mass conservation on the clipping parameter ( ) employed. The CIAM-Superbee combination again seems to be impervious to changes in clipping, whereas WY-QUICK seems to perform slightly better by lowering the parameter. . A schematic of the computational domain. An infinitely long and axisymmetric corrugated ligament of mean radius R is placed in a square box of side L. The bottom side of the box acts as the axis of symmetry for the ligament, with spatial periodicity along the horizontal direction. Inset : A close up view of the corrugated profile of the ligament, where the local radius is defined as the sum of the unperturbed (mean) radius R and the local perturbation (x).

The material properties of the liquid and gas phases are denoted with the subscripts l and g respectively, which in our case corresponds to an air-water system with the surface tension coefficient σ. . . . . . . . . . . . . . . . . . . . . Radial profile of the noisy and the subsequently low-pass filtered surface of the ligament, corresponding to one spatial period along its axis. Inset : Spectral densities of the initial noisy signal and the corresponding filtered signal. Only a handful of small wavenumbers are chosen to create the surface. . . . . . . . . The dispersion relation of the Rayleigh-Plateau instability [ , ], expressing the growth rate of the small amplitude perturbations as a function of the reduced wavenumber ( . ). The vertical dashed line A represents the upper limit to the perturbation wavelength, corresponding to the domain size L. At the opposite end of the spectrum, E represents smallest wavelength limited by the grid size 1/2 M , where M is the maximum level of resolution for our octree grid. B corresponds to the wavenumber with the maximum growth based on linear theory, whereas C represents the limit beyond which all perturbations are stable. The number of discrete modes that fall within the unstable part of the spectrum (∆K) is a linear function of the ligament aspect ratio. . . . . . . . Dynamically adapted octree meshes in the periphery of the interface location, with the refinement criteria based on limiting gradients of volume fraction and velocity. The interface is represented by the white contours, the colormap on the left half is based on the axial velocity component, whereas the one on the right corresponds to that of vorticity. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. . . . . Illustration of the effects of our "droplet removal" feature, with the figure on the left corresponding to simulations that filter out structures that are resolved by cells or less, with the one on the right showing results with no droplet removal in place. The box (dashed line) in the right hand side figure displays the extremely small structures that consume a disproportionate amount of computational resources, without having any meaningful impact on the dynamics that we are interested in. The colormap on the left side of the axis of symmetry corresponds to the axial component of velocity, whereas the one on the right refers to the level of numerical resolution. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. . The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes. . . The tail of the distribution fitted with the function in the exponential form ( . ), using the optimal parameter set that corresponds to the lowest standard deviation in the parameter estimate. The uncertainty in the bin heights of the underlying distribution is taken into account while conducting the least-squares procedure. The bottom figure is a closeup of the region in the vicinity of the tail. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes. . . The variation in the optimal parameter m as a function of the number of points in the window (see . ). The least-squares procedure is performed on the candidate function in the logarithmic ( . ) form. The error bars represent the standard deviation of the parameter estimate obtained from the leastsquares routine, where the value enclosed in the red box corresponds to the minimum error. The figure on top (a) corresponds to the least-squares procedure where the uncertainty in the bin heights is taken into account, whereas the bottom figure (b) corresponds to those carried out with zero uncertainty in the underlying distribution. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes. . . . . . . . . . . . . . . . . . . . The tail of the distribution fitted with the function in the logarithmic form ( . ), using the optimal parameter set that corresponds to the lowest standard deviation in the parameter estimate. The uncertainty in the bin heights of the underlying distribution is not taken into account while conducting the least-squares procedure. The bottom figure is a closeup of the region in the vicinity of the tail. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes. . . The tail of the distribution fitted with the function in the logarithmic form ( . ), using the optimal parameter set that corresponds to the lowest standard deviation in the parameter estimate. The uncertainty in the bin heights of the underlying distribution is taken into account while conducting the least-squares procedure. The bottom figure is a closeup of the region in the vicinity of the tail. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes. .

Introduction

The dynamics of liquid-gas interfacial flows play a critical role in several processes in nature, as well as in myriad industrial applications. The key elements of surface tension dominated flows such as droplets and bubbles constitute the fundamental mechanisms governing the exchange of heat and mass at the ocean-atmosphere interface [ , ], mixing/separation in the context of metallurgical processes [ , ], conventional modes of heat transfer [ , ] and ever so importantly, the transmission of pathogens [ , ]. One of the most fascinating features of multiphase flows is the process of atomization, in which a liquid volume breaks into smaller fragments via a series of topological changes of varying complexity, ultimately resulting in the emergence of drops of various sizes driven by the action of capillary forces at the interface separating the fluids. Such processes are ubiquitous in a diverse range of applications spanning from combustion related processes ([ , ]) to agricultural irrigation ([ , ]). In view of the broad spectrum of liquid fragmentation phenomena that pique our scientific interest, the present body of work is divided into two parts :

Development of numerical methods that can reproduce the dynamics of liquid-gas interfacial flows at low to moderate computational cost (spatial resolution), aimed towards flow configurations involving significant contrasts in material properties across the interface. Application of numerical methods in an effort to quantify the influence of certain topological characteristics of liquid structures on the resulting drop size distributions.

In what follows, we take a closer look at the challenges associated to the above mentioned themes.

Challenges in Numerical Modeling

A substantially large subset of all surface tension dominated flows involve significant disparities in the material properties across the interface, the most common example being flow configurations corresponding to air-water systems, where the densities and viscosities of the fluids are separated by (approximately) and orders or magnitude, respectively. The development of numerical methods that attempt to model such interfacial flows involving marked contrasts in density face several challenges, key amongst them being the transport of mathematical discontinuities that arise out of the aforementioned contrasts. Extremely small numerical errors are ubiquitous as a consequence of the numerous approximations involved at each and every step of the algorithm (e.g. propagation of the interface, curvature computation, surface tension modeling etc.) . In the context of such flows, such "numerical errors" may result in physically inconsistent mass and momentum transfer across the interface, often from the denser phase towards the lighter phase as a : By marked density contrasts, we mean large | log r |, where r is the ratio of the densities of the two fluids. consequence of inadequate numerical resolution. The presence of large density contrasts tend to amplify the growth of these cascading numerical errors, eventually leading to significant (often catastrophic) interfacial deformations followed rapidly by a loss of numerical stability.

A prime example of such a situation is the numerical simulation of a raindrop falling in air under the action of gravity, at speeds near terminal velocity. For a regime corresponding to low Weber numbers (≈ 3 in our case), there is abundant experimental evidence (e.g. see [ ]) demonstrating that the drops do not undergo any significant deformation, apart from some minor shape oscillations driven by the interplay between surface tension and inertia. Fig. . shows the results from a moderately resolved numerical simulation, where the rapid and uncontrolled growth of numerical errors lead to severe deformations of the interface, finally resulting in the un-physical (read "artificial") atomization of the raindrop into numerous smaller fragments.

In typical liquid atomization processes, the scale separation between the sizes corresponding to the largest and smallest liquid structures are often several orders of magnitude, which inevitably lead to the smaller structures often being poorly resolved due to the inherent limitations in computational power. Thus, one can consider that the case of the (less than well resolved) falling raindrop to be representative of many of the small structures that emerge out of liquid fragmentation phenomena. In that regard, the occurrence of "artificial" atomization may significantly distort our understanding of the distribution of drop sizes, as there is no meaningful way to distinguish between drops generated via physically consistent processes and those produced as a result of cascading numerical errors.

Ever since such numerical instabilities were first observed with the "SURFER" [ ] code in the context of planar jets and rising bubbles, considerable efforts have been made towards the design of numerical methods to specifically deal with flows involving such marked density contrasts . The underlying principle behind these endeavours is that the governing equations for the transport of mass and momentum are modeled using a conservative formulation (divergence of fluxes), instead of standard non-conservative forms which themselves were adapted directly from techniques developed originally for single-phase flows. This formulation enables one to render the discrete transport of momentum consistent with respect to the discrete transport of mass. Such a tight coupling of the propagation of errors between the discrete mass and momentum fields enables alleviation of many of the issues that plague such numerical methods, especially in the context of low to moderate spatial resolutions (see Fig.

. ).

A difference in the drop size distributions obtained via consistent and non-consistent methods has been observed by Yue Ling and coworkers [ , ] in the context of planar atomizing mixing layers. Use of the non-consistent method yielded in an exponential distribution for the large sizes, whereas with the consistent method they recovered a Log-Normal distribution.

Exacerbating the already complicated nature of the discrete transport of material discontinuities is the role of capillary forces on the evolution of the interface. They are commonly modeled as singular source terms in the momentum balance equation that governs the evolution of the velocity field, with the capillary force itself being proportional to the third derivative of the interfacial position. Thus, a secondary but important mitigating factor is the advancements made in the modeling of capillary forces, resulting in the adoption of consistent and well-balanced surface tension formulations. Consistency in the context of surface tension models : Mass advection was carried out using using algebraic flux reconstructions. refers to the ability of methods to progressively achieve more accurate estimations of interfacial curvature as a result of increasing spatial resolution, whereas well-balanced refers to the ability to recover certain static equilibrium solutions pertaining to surface tension dominated flows without the perpetual presence of parasitic or spurious currents in the velocity fields. We refer the reader to influential works of Popinet [ , ] to get a better understanding of the issues surrounding different surface tension implementations.

Figure . :

Numerical simulations of a 3mm raindrop falling in air under gravitational acceleration. The flow is along the positive X direction, with gravity along the opposite direction. The red contour indicates the isosurface of the volume fraction field corresponding to a value of 0.5, whereas the black contours surrounding the drop represent isosurfaces of the magnitude of vorticity. The droplet has a relatively poor numerical resolution, corresponding to D/h 16, where D is the drop diameter and h is the grid size. The top row illustrates the phenomenon of "artificial atomization" that occurs while using the standard version of our numerical method, which does not maintain consistency between the discrete transport of mass and momentum. In contrast, the implementation of consistent mass and momentum transport enables us to suppress the massive interfacial deformations of the drop, as evidenced by the bottom row of figures depicting the same flow configuration. The development of methods concerning the consistent transport of mass and momentum is one of the central themes of this study.

The first study to address the issue of consistency between mass and momentum transport was the seminal work of Rudman [ ] [ ]: Rudman ( ), 'A volume-tracking method for incompressible multifluid flows with large density variations' . The fundamental hurdle in the implementation of mass-momentum consistent transport for staggered configurations of primary variables (pressure and velocity) is the inherent difficulty in reconstructing mass (defined on centered control volumes) and its corresponding fluxes onto the staggered control volumes on which momentum is defined. Rudman introduced the strategy of carrying out mass advection on a grid twice as fine as that of momentum, thereby enabling a 'natural' and intuitive way to reconstruct mass and its fluxes onto staggered momentum control volumes. However, their method uses a VOF based convolution technique for curvature computation, which is neither consistent nor well-balanced. ), 'Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method' , followed by Ghods and Hermann [ ]. In the former, the consistency problem is tackled by means of a semi-Lagrangian approach, computing geometric level set derived fluxes at two different time intervals, whereas in the latter, a collocated arrangement is used. Nonetheless, both methods face certain drawbacks, notably the applicability only to D in the case of Raessi and Pitsch, as well as a lack of well-balanced surface tension models for both these methods. , they use mesh-decoupled height functions to compute curvature while coupling it with a well-balanced surface tension model. However, their semi-Lagrangian flux computation procedure involving streak tubes and flux polyhedra are extremely convoluted in D.

Bussmann et al. [ ]

Recent advances concerning

Certain methods attempt to combine the qualities of both VOF and level set methodologies (CLSVOF), as proposed in the works of Vaudor et al. [ ] [ ]: Vaudor et al. ( ), 'A consistent mass and momentum flux computation method for two phase flows. Application to atomization process' , and more recently by Zuzio et al. [ ]. They both tackle the consistency issue by means of projecting the direction-split geometric fluxes onto a twice finer grid, which are subsequently recombined to reconstruct consistent fluxes for mass and momentum for the staggered control volumes. This approach allows them to bypass the requirement of conducting mass advection on a twice finer grid (as in the original Rudman method), thereby deriving the benefits of a sub grid without the added computational cost of doubly refined mass transport. In addition, both Vaudor et al. [ ] and Zuzio et al. [ ] [ ]: Zuzio et al. ( ), 'A new efficient momentum preserving Level-Set/VOF method for high density and momentum ratio incompressible two-phase flows' adopt well-balanced surface tension models with consistent level set based curvature estimation. However, the purported advantages of both these methods with regard to reduced computational costs is not quite evident, as additional complexities are introduced due to the projection (reconstruction) of fluxes onto a the twice finer mesh, which would not be necessary in the first place if mass transport had been carried out on the twice finer mesh itself. the additional mass 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 constructions being an open question .

Patel and Natarajan [ ]

To get a bird's eye view of the numerous features employed by the methods in existing literature, we refer the reader to tables . and . , which respectively provide systematic overviews of the VOF and level set based approaches. In the present body of work, we start by precisely defining the essential (desired) attributes of a numerical scheme that ensures discrete consistency between mass and momentum transport :

The discontinuity of the mass (volume-fraction derived density field) should propagate at the exact 'numerical' speed as that of the discontinuity of the momentum field. The numerical transport of momentum should be performed in a manner consistent with the transport of mass for each direction, implying that the momentum fluxes must be obtained directly from geometrically computed fluxes of mass (volume).

In order to tackle the challenge of consistent transport on staggered control volume configurations, we have developed two different strategies, namely, the shifted fractions method and the sub-grid method. The former uses geometrical reconstructions to derive a shifted volume fraction field which is centered on the staggered control volumes, whereas the latter adopts the Rudman [ ] strategy of volume advection on a twice finer grid in order to enable consistent reconstruction of mass and momentum on the staggered control volumes. Another key contribution of this body of work is to extend the conservative direction-split mass transport algorithm of Weymouth and Yue [ ]

Ligament Mediated Paradigm

Liquid atomization is basically the transformation of a compact volume into drops. However, this simplistic view masks the intricate interplay between inertia, viscosity and capillarity across different length and time scales spanning several orders of magnitude. Such non-trivial interactions are responsible for the abundant variety in the sequence of topological progressions that eventually lead to the formation of stable drops. The most basic transition involves the breakup of a cylindrical filament structure (see Fig. . ) at approximately regular intervals, driven primarily by the growth of long wavelength perturbations due to capillary forces. A slightly more complicated transition involves liquid sheets (refer to Fig.

. ), where the inertial expansion of the sheet is opposed by the capillary deceleration of the edges, resulting in the formation of liquid rims due to volume accumulation at the edges, the subsequent destabilization of which leads to the generation of drops. Sticking to the breakup of liquid sheets, another important transition involves the appearance of perforations or holes (see Fig. . ) when the thickness of the sheet reaches a certain limit. The ensuing rapid expansion of such perforations due to capillary retraction results in the generation of drops. Arguably, the most convoluted route to drop formation in encountered when : Refer to Wang and Bourouiba [ ] for a discussion. macroscopic liquid structures are subjected to shear-driven instabilities (see Fig. . ) that arise at the interface due to the differential gas and liquid velocities. These primary instabilities of the Kelvin-Helmholtz [ ] variety lead to the creation of a plethora of secondary structures, many of them corresponding to the aforementioned topologies like filaments, expanding sheets with rims, expanding holes in thin sheets etc. amongst numerous others.

(i) (ii)
The development of quantitative models geared towards predicting the statistical properties of certain attributes of the drops (sizes, shapes, velocities, etc.) generated from such fragmentation phenomena has drawn considerable scientific interest throughout the past few decades. Numerous investigations (both experimental and numerical) into the statistical description of drop sizes have led to the popularization of three distinct classes of probability density functions, namely the Log-normal, Gamma and Poisson distributions [ ]. In addition to the aforementioned families, distributions such as the Weibull [ ], Exponential [ ] and Beta [ ] have also received significant attention, particularly (but not limited to) in the context of raindrop sizes. In terms of physical interpretations, the Log-normal model [ ] is based on a sequential cascade of breakups, the Gamma family [ ] based on the competing effects of aggregation and coalescence of liquid volumes, and the Poisson model [ ] based on the instantaneous and random splitting of a volume into smaller constituents. Even though the aforementioned models have been applied to a wide variety of flow configurations with varying degrees of success, there is a general lack of consensus regarding their validity, primarily due to the markedly different trajectories followed by the initial liquid structures towards eventual drop formation. Nevertheless, the one common feature that unites these seemingly disparate fragmentation processes we have encountered thus far is the transformation of the initial liquid topology into transient columnar structures called ligaments. It is the evolution and concomitant disintegration of these ligaments that finally lead to the production of drops, therefore we entertain the possibility that the dynamics of such structures may provide us with clues towards a "universal" mechanism that can explain the polidispersity in drop sizes.

The topological change from the thread-like structure of the ligaments to the (approximately) spherical geometry of drops can proceed along different paths, depending on whether the drop forms at the free ends or through the rupture of the main liquid column at different locations. The former mode is referred to as end-pinching [ , ], which generally occurs for regimes (Oh < O(10 -2 )) where the capillary driven retraction of the ligament tips cannot be rapidly damped by the viscous forces. Additionally, if the ligament is free at both ends and not slender enough (small elongation ratios), there is a possibility that the capillary retraction contracts the entire volume into a single drop [ ]. The dynamics of the end-pinching mode is further complicated by the interaction of additional parameters such as the properties of the surrounding fluid and rate of intertial stretching. Despite the richness in the dynamics of end-pinching, the mechanism determining the size of the resulting drops is quite well understood, with robust scaling laws establishing the drop diameters to be a linear function of the width of the original ligament [ ] (although there is some weak (W e -1/7 ) dependence on the inertial stretching rate [ ] ). This turns our attention solely towards the destabilization of the ligament bulk.

The mechanism of drop formation via the rupture of a liquid column is fundamentally non-linear in nature [ , ], corresponding to finite time singularities of the Navier-Stokes equations. The fluid flow in the vicinity of the point of breakup along the liquid column is essentially governed by self-similar solutions [ ], whose local nature (independent of initial conditions) is attributed to the intrinsic non-linearities of the underlying equations. The exact form of such solutions depends on the relative strengths of inertia, capillary forces and viscosity at that particular length scale, which can be succinctly characterized by a point in the phase space of Weber and Ohnesorge numbers. The regime in which inertial forces balance capillarity is commonly referred to as the inertial regime [ , ] and the one in which viscosity balances the capillary forces is termed as the viscous regime [ , ], where both regimes correspond to the evolution of the minimum thickness as self-similar solutions of the second kind ([ ]) due to their dependence on an external length scale. Finally, at moments very close to the pinch-off, neither of the above regimes are valid as the thickness of the thread approaches the intrinsic length scale (corresponding to Oh 1) of the fluid where inertia, capillary and viscous forces are in perfect balance [ ]. Here one observes the "universal" inertio-viscous regime [ ], although the exact trajectory towards this final shape before pinch-off involves complicated to and fro transitions between the inertial and viscous regimes on extremely fast time scales [ ]. Despite the fact that the non-linear dynamics near pinchoff provides us with insights into the mechanism of drops formation, including the production of "satellite" drops [ -], it still does not explain the dispersion in the sizes of drops formed via rupture of the ligament along several points. Prior to the onset of breakup, the flow inside the ligament is driven by the several factors such as the inertial stretching, propagation of capillary waves along the surface, shear-driven interactions at the interface etc. This brings us to the second theme of the present study, where we focus on quantifying the impact of the rearrangement of liquid volumes within the ligament on the resulting polydispersity in the size of the drops generated. Towards this objective, our approach involves conducting well resolved direct numerical simulations of spatially periodic ligaments, where the periodicity enables us to suppress competition from the endpinching mode. The specifications of the initial conditions involve length scales, material properties of the fluids, and most importantly, generation of random, but unique geometrical shapes of the ligament surface using sets of randomly overlapping waves. We incorporate stochasticity into our numerical investigation by carrying out ensembles of such numerical simulations of slender ligaments with randomly perturbed shapes, in which each realization corresponds to a random but unique initial configuration.

An outline of the ensuing chapters is as follows:

Chapter provides details about the governing equations, basic numerical techniques regarding the transport of interfaces (VOF), and a thorough exposition of our class of numerical methods that maintain consistency between the discrete transport of mass and momentum in the context of staggered Cartesian grids. Chapter evaluates the performance and accuracy of the different numerical methods we have developed thus far, presenting several quantitative comparisons against certain benchmark flow configurations. Chapter demonstrates the robustness and stability of our numerical methods in the context of a challenging flow configuration in which a raindrop falls in the air under gravitational acceleration. Chapter provides a thorough description of our numerical setup for the spatially periodic ligaments, the methodology behind the generation of randomly perturbed (corrugated) surfaces and the impact of initial conditions on the drops generated due to the ligament breakup. Chapter presents a statistical description of drop sizes, including theoretical predictions based on random Gaussian processes and detailed explorations at the large size limits of the distributions. In this chapter, we describe our numerical methodology behind modeling the dynamics of immiscible incompressible liquid-gas interfacial flows under isothermal conditions. The numerical implementation is based on finite volume discretizations on uniform Cartesian grids, utilizing state of the art methods in interfacial reconstruction coupled with geometric transport of the corresponding fluxes, curvature computation and surface tension modeling. A detailed exposition of our class of mass-momentum consistent numerical methods is provided, which are specifically designed to circumvent or suppress the uncontrolled and rapid growth of numerical instabilities that arise when dealing with flows entailing marked density contrasts . The implementations of the algorithms are developed on the free and open-source numerical platform called 'PARIS Simulator' [ ], with the detailed descriptions of the general capabilities of the solver to be found in the previously cited reference.

. Governing Equations

We use the one-fluid formulation for our system of governing equations, thus solving the incompressible Navier-Stokes equations throughout the whole domain including regions of variable density and viscosity, which itself depend on the explicit location of the interface separating the two fluids. In the absence of mass transfer, the velocity field is continuous across the interface at the incompressible limit, with the interface evolving according to the local velocity vector. Thus, the equations are as follows :

∂ρ ∂t + ∇ • ρu 0 ( . ) ∂ ∂t ρu + ∇ • ρu ⊗ u -∇p + ∇ • 2µD + σκδ s n + ρ g ( . )
with ρ and µ being the density and dynamical viscosity respectively. The volumetric sources are modeled by the acceleration g, and the deformation rate tensor D used to model the viscous stresses defined as:

D 1 2 ∇u + (∇u) T ( . )
The term σκδ s n models the surface tension forces in the framework of the continuum surface-force (CSF) method [ ]. The normal vector to the interface is n, σ the coefficient of surface tension and κ the local interfacial curvature. The operator δ s is the Dirac delta function, the : Generally, χ is assigned values of 0 in one phase and 1 in the other.

: VOF based methods implicitly track the evolution of the discontinuous density field, which is not the case in level set based methods.

: The differentiable nature of the level set function lends itself to straightforward curvature computation routines. numerical approximation of which allows us to map the singular surface force distribution along the interface onto their volumetric equivalents for our Cartesian control volumes. At the incompressible limit, the advection of mass given by equation . can be treated as equivalent to that of the advection of volume.

Material Properties

Within the framework of interface capturing schemes, the temporal evolution of the interface separating the two fluids can be tracked by the following advection equation :

∂χ ∂t + u∇χ 0 ( . )
where χ is the phase-characteristic function, that has different values in each phase . Mathematically, the function χ is equivalent to a Heaviside function in space and time. At the macroscopic length scales under consideration, the interface evolution as described by equation

. is modeled as having infinitesimal thickness under the continuum hypothesis. The coupling of the interfacial evolution with the equations of fluid motion as described in . and . is provided by :

ρ ρ 1 χ + (1 -χ) ρ 2 ( . ) µ µ 1 χ + (1 -χ) µ 2 ( . )
where ρ 1 , ρ 2 are the densities of fluids and respectively, likewise for viscosities µ 1 and µ 2 . For certain flow configurations, it might be beneficial to opt for a weighted harmonic mean description of the variable dynamic viscosity [ ] .

Interfacial Transport : VOF

Our numerical studies are based on the Volume-of-Fluid methodology.

We refer to the discontinuous approximation to the Heaviside function χ as the volume fraction field or colour function interchageably, which is defined below in the context of finite volume discretization :

C i jk (t) 1 ∆V ∫ ∆V χ(x, t) dx ( . )
where C is the colour function with its values lying between 0 and 1, with i,j and k being the indices to the corresponding discretized control volume of volume ∆V. There are two steps involved in the VOF method, the reconstruction of the interface and its subsequent propagation (advection). We present a brief overview of the two steps in the following sections, as going into detailed descriptions of the reconstruction and propagation procedures are not the focus of the present body of work .

PLIC representation

We use geometric reconstructions to explicitly define the interface location using the discrete colour function information.

Flux Computation

Once the geometric PLIC reconstructions have been carried out, the interface segments are advected using the the velocity field. This entails computation of fluxes of the colour function, which can be computed via : The reader can refer to [ , ] for further details : A detailed exposition of this problem along with a noteworthy solution can be found in the work by Weymouth and Yue [ ] : For details regarding geometric flux reconstruction, refer to [ , ] algebraic transport schemes (generally less accurate), or by using geometric reconstructions in either Eulerian, Lagrangian or hybrid frameworks. In the context of our numerical platforms ('PARIS' and 'Basilisk'), stateof-the-art geometrical flux reconstruction procedues are utilised. The temporal integration of the fluxes could be carried out either as a series of one dimensional 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 of the one dimensional algorithm to D), but suffer from lack of conservation (to the order of machine precision) when it comes to D . Multidimensional (unsplit) methods have an advantage in that respect due to the fact that they are conservative by nature of their design, but are inherently more complicated to develop and implement, with no straightforward extension from D to D. For a more detailed and nuanced evaluation of the comparative advantages of interfacial transport methods, we refer the reader to the recent review by Mirjalili et al. [ ] on the given subject. The propagation of the interface can be described by the evolution of the colour function (volume fraction field) as

∂C ∂t + u • ∇C 0 ( . )
We can express the left hand side of ( . ) in the conservative form as

∂C ∂t + ∇ • (Cu) C (∇ • u) ( . )
As one can observe, the "compression" term on the right hand side of equation . 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 the direction-split framework. Discretization of the above equation results in :

C n,d+1 i, j,k C n,d i, j,k -F n,d + (C) -F n,d -(C) + C n,d i, j,k ∆u q ∆x n,d i, j,k ( . )
The above equation represents an advection substep, the superscripts n and d refer to the timestep and direction of integration respectively.

The notation d 0 refers to the field at the n th timestep, before any integration is performed along any direction. The fluxes F n,d ± (C) in equation . are derived through geometrical reconstructions . The + andsubcripts refer to the orientation with respect to the central cell (i, j, k). The subscript q refers to the direction corresponding to that given advection substep i.e either X, Y or Z. In our case, our numerical platforms are based on cubic (regular Cartesian) grids, consequently there is no requirement for a subscript with ∆x. After each substep, the interface is reconstructed once again with the updated volume fraction : The classification of Lagrangian explicit and Eulerian implicit are in accordance with the paper by Aulisa et al. [ ] field in order to compute the fluxes for the next advection substep. Finally, the volume fraction field for the next timestep is given by

C n+1 i, j,k C n,3 i, j,k ( . )
The interpretation and numerical approximation of the prefactor C n,d i, j,k to the directional divergence, as well as the fluxes F n,d ± (C), depend on the exact nature of the geometrical advection scheme in question, which in our context is either CIAM (Lagrangian explicit) or Weymouth-Yue (Eulerian implicit) . A brief outline of these two methods is presented in the subsequent sections. As one can infer from the geometrical representation in figure . , the fluxes F n,d ± (C) correspond to the volumes Φ -and Φ + . Thus, the updated field C n,d+1 i, j,k is the sum of the three contributions Φ -, Φ 0 and Φ + . We can rewrite equation . specifically for the CIAM scheme as

Lagrangian Explicit

C n,d+1 i, j,k C n,d i, j,k 1 + ∆u q ∆x n,d i, j,k -F n,d + (C) -F n,d -(C) ( . )
where the compression coefficient C n,d i, j,k is simply equal to the value of the colour function C n,d i, j,k at the start of the corresponding advection substep. Although the flux terms cancel upon integration throughout the whole domain, one can clearly see that the compression terms do not sum up to zero due to the changing prefactor in front of the directional divergences. This precise issue brings us to the next advection scheme.

: By 'discrete conservation' we mean that the sum of the directional divergences sum upto zero, to the accuracy of machine precision.

: Compression coefficient is used as a short-hand version of 'prefactor to the compression term'.

: In numerical terms, we can only ensure that they sum up to the accuracy of the Poisson solver ∼ 10 -3 -10 -6 , with the limiting factor being the level of machine accuracy (∼ 10 -14 -10 -17 ).

: For a proof of discrete volume conservation subject to certain CFL criteria, refer to the appendix of Weymouth and Yue [ ].

Eulerian Implicit This advection scheme was developed by Weymouth and Yue [ ] in order to specifically tackle the problem of discrete conservation when it comes to direction-split geometrical advection schemes . The scheme fundamentally employs a forward Eulerian method in order to carry out temporal integration of the fluxes, with the fluxes themselves computed as the quantity of the substance entering or exiting a given control volume through its fixed surfaces, as shown in figure . . This is in contrast with the flux computation method in the case of CIAM, where the interface segments are propagated forward in time in a Lagrangian fashion.

Figure . :

A D schemetic of the Eulerian (geometric) flux calculation using the Weymouth-Yue [ ] scheme for the advection substep along the horizontal direction, with the interface reconstructed using the volume fraction field at the start of the substep. The colour fraction of 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 -1, j and the red one corresponding to that exiting i, j into i + 1, j. The geometric flux calculations are made on the basis of the interfacial positions at the start of the substep, and the face centered velocities of the cell in question.

The subtle but important tactic used in this scheme lies in the manner in which the prefactor to the compression term is treated, with its definition being :

C n,d i, j,k H C n,0 i, j,k -1/2 ( . )
where H is a one-dimensional Heaviside function. This renders the compression coefficient independent of the direction of the advection substep, consequently enabling the three discrete directional divergences to sum up to zero Therefore, the scheme is able to demonstrate volume conservation, subject to local CFL restrictions . To summarise, we can rewrite equation . for the Weymouth-Yue scheme as

C n,d+1 i, j,k C 1 + ∆u q ∆x n,d i, j,k -F n,d + (C) -F n,d -(C) ( . )
where C is a constant with a value of either 0 or 1, determined by the value of C n,0 i, j,k according to equation . .

: also referred to as the strong form, necessitating certain orders of smoothness of the primitive variable : These methods are descendants of the class of numerical schemes used to solve hyperbolic partial differential equations.

: Refer to [ ] for a more detailed discussion.

. Consistent Transport of Mass and Momentum

Generally, we have a choice regarding how to discretize the convective operator of the incompressible Navier-Stokes equations. There is a well established corpus of numerical methods tailored specifically to deal with the non-conservative form of the convective operator that appears in the transport equations of mass and momentum , which perform quite well in the context of single phase flows. However, in interfacial flows we often deal with sharp discontinuities that arise as a consequence of the contrast in material properties between the two fluids. Therefore, even though the velocity field remains continuous throughout the domain, the otherwise smooth density and momentum fields contain sharp jumps localized at the interfacial position. Therefore, in regions of constant density and and viscosity, the standard non-conservative form of the Navier-Stokes equations are discretized, whereas the discretization is reformulated in a conservative form at locations in the vicinity of the interface.

Spatiotemporal Discretization

We start by describing the spatial arrangement of our primary variables i.e. pressure (p) and velocity (u q ), where q 1, 2, 3 represents the x,y and z componenets respectively. The control volume is in the form of a cube ( D) or a square ( D). The pressure and velocity variables are defined in a staggered arrangement, which means that the control volumes corresponding to the velocity components u 1 and u 2 are shifted with respect to the control volume for pressure. The staggered grid with the primary variables are represented in Figure . . The use of staggered control volumes has the advantage of suppressing neutral modes (pressure oscillations) often observed in collocated methods but leads to more complex discretizations .

In order to describe the overall numerical algorithm for the one-fluid Navier-Stokes equations with variable density and viscosity, we choose to the reframe our equations in a more convenient operator form, as presented below :

∂ ∂t ρu L ρ, u -∇p ( . )
The operator L in the above expression can be decomposed as :

L L adv + L µ + L σ + L g ( . )
where the L adv represents the conservative advection, L µ represents the diffusive forces generated by viscous stresses, L σ represents the capillary forces arising from the surface tension model and finally L g represents the volumetric (body forces) source term.

u 2;i,j-1/2 u 1;i-1/2,j u 1;i+1/2,j u 2;i,j+1/2 p i,j-1 p i+1,j p i-1,j p i,j p i,j+1 u 1;i-1/2,j+1
Figure . : A D schematic of the staggered spatial configuration of the pressure and velocity variables. The pressure p i, j is based on the center of its control volume (light colour area); the horizontal velocity component u 1;i+1/2, j is defined in the middle of the right edge of the pressure control volume and centered on its control volume (dash-dotted area); the vertical velocity component u 2;i, j+1/2 is defined in the middle of the top edge and centered on its corresponding control volume (dashed area). The volume fraction is defined at the same location as that of the pressure, resulting in pressure and density being centered on identical control volumes.

: These high-order spatial schemes are based on well established methods developed to deal with hyperbolic conservation laws, for more details refer to the studies of Leveque [ ] and Sweby [ ].

We apply the spatially discretized versions of these operators (denoted by the superscript h) onto the primary variables (C, u), and march forward in time using a small, possibly variable time-step τ such that t n+1 t n + τ.

In the first part of the algorithm, the volume fraction field C n is updated to the next timestep , with the superscript n signifying discretization in time. The operation can be written as follows

C n+1 C n + τL h vof (C n , u n ) ( . )
The subscripts i, j, k are dropped in equations ( . ) and ( . ), with the understanding that the operators in equation . apply uniformly to all control volumes.

The temporal evolution of the volume fraction field represented above by the operator L h

vof is in accordance with the Lagrangian explicit or Eulerian implicit advection schemes, as described in the previous sections. Once we have obtained the updated field C n+1 , we can move on to the temporal update of our momentum field given by

ρ n+1 • u * ρ n • u n + τL h adv (C n , u n ) + τ L h µ C n+1 , u n + L h σ C n+1 + L h g C n+1 ( . )
In regions of constant density. the advection operator L h adv is implemented using higher order spatial schemes coupled with a choice of non-linear flux limiters such as QUICK, ENO, WENO, Superbee, Verstappen and BCG. . For control volumes in the vicinity of the interface location, we revert to lower order schemes due to the sharp jumps in the material properties across the interface. The functionality of the operator L h adv near the interface is tighly coupled to that of L h vof from equation . , so as to ensure consistency in the discrete advection of mass and momentum. The details regarding this coupling shall be the focus of this section.

: These fluxes of (uH) are computed via geometrical reconstructions, and not via high-order non-linear interpolation schemes which are the standard (in the absence of discontinuities) .

: Component of the velocity normal to the control surface f .

Advection of Conserved Quantities

We start with the advection of the interface position, which in the VOF framework is given by -In the VOF approach, the approximation to the interface Heaviside function contains a sharp discontinuity, in contrast to the level set approach where the approximation is continuous in the discrete sense.

∂H ∂t + ∇ • (uH) H (∇ • u) ( . )
Integrating this equation in time after carrying out spatial discretization, one obtains -

C n i, j,k
is the colour function or volume fraction field at time step n obtained through the finite volume discretization of the interface Heaviside function.

C n+1 i, j,k -C n i, j,k - faces f F (c) f + ∫ t n+1 t n dt ∫ Ω H (∇ • u) dx , ( . )
where the first term on the right-hand side is the summation over the cell faces f of the fluxes

F (c) f
. As one can clearly observe, the "compression" term on the right-hand side of equation ( . ) disappears for incompressible flow, however it is essential in the context of directionsplit geometric advection schemes. The definition of the geometric fluxes in mathematical form is -

F (c) f ∫ t n+1 t n dt ∫ f u f (x, t)H(x, t) dx , ( . )
where u f u • n f . Once an approximation for the evolution of (uH) during the time step is chosen, a four-dimensional integral remains to be computed in equation ( . ). The advection schemes used are CIAM and Weymouth-Yue, which have been described in the previous section. Directional splitting results in the breakdown of equation ( . ) into three equations, one for each advection substep -

C n,l+1 i, j,k -C n,l i, j,k -F (c) m--F (c) m+ + c m ∂ h m u m , ( . )
The superscript l 0, 1, 2 is the substep

index, i.e. C n,0 i, j,k C n i, j,k and C n,3 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-≥ 0 if the flow is locally from right to left. A similar reasoning applies to the "right" face m+.

After each advection substep ( . ), the interface is reconstructed with the updated volumes C n,l+1 i, j,k , then the fluxes

F (c)
f are computed for the next substep. Importantly, we have approximated the compression term in ( . ) by -

∫ t n+1 t n dt ∫ Ω H∂ m u m dx c m ∂ h m u m ( . )
There is no implicit summation carried out over m, and the superscript h denotes the spatial discretization of the operator.

On the right hand sides of ( . ) and ( . ) the flux terms

F (c)
f and the partial derivative ∂ m u m must be evaluated using identical discretized velocities. The expression ∂ h m u m is a finite volume approximation of the : The preservation of 0 ≤ C i, j,k ≤ 1 is referred to as C-bracketing.

spatial derivative corresponding to the mth component of the velocity vector along direction m, and the "compression coefficient" c m approximates the color fraction. Its exact expression is dependent on the advection method and it also entails the desirable property of C-bracketing . Due to the possible dependence of the compression coefficient c m on the Cartesian direction corresponding to the advection substep, the sum m c m ∂ h m u m may not necessarily vanish even if the flow is incompress- ible. As mentioned previously, the main appeal of the Weymouth-Yue method is its ability to keep this sum to zero, within the limits of machine precision. In order to achieve consistent transport of the discontinuous fields of mass and momentum, we start by trying to understand the advection of a generic conserved scalar quantity φ by a continuous velocity field -

∂φ ∂t + ∇ • φ u 0 ( . )
where the field φ is smoothly varying except at the interface position, where it may be discontinuous. The essence of this study lies in the search of a scheme that propagates this discontinuity (φ) at the same speed as that of the advection of volume fraction (C).

The smoothness of the advected quantity away from the interface is verified for fields such as density ρ, momentum ρu and internal energy ρe.

Temporal integration of the spatially discretized version of ( . ) gives us -

φ n+1 i, j,k -φ n i, j,k - faces f F (φ) f ( . )
The sum on the right-hand side is the sum over faces f of cell i, j, k of the fluxes

F (φ)
f of φ, which are defined as the color function fluxes

F (c) f in ( . ).
where the fluxes are defined as -

F (φ) f ∫ t n+1 t n dt ∫ f u f (x, t) φ(x, t) dx ( . )
In order to "extract" the discontinuity we introduce the interface Heaviside function H(x, t)

F (φ) f ∫ t n+1 t n dt ∫ f u f H φ + u f (1 -H) φ dx ( . )
Therefore, the flux can be decomposed into two components as -

F (φ) f φ1 ∫ t n+1 t n dt ∫ f u f H dx + φ2 ∫ t n+1 t n dt ∫ f u f (1 -H) dx , ( . )
where the face averages φs , s 1, 2, are defined as -

F (φ) f φ1 F (c) f + φ2 F (1-c) f . ( . )

Advection of Density

The density field ρ(x, t) follows the temporal evolution of the generic conserved quantity ( . ) by setting φ ρ. At the incompressible limit the velocity field is solenoidal (divergence-free), with constant densities in each phase. We can extract the density trivially from the integrals ( . ) to exactly obtain ρs ρ s . The flux definitions corresponding to ρ become

F (ρ) f ρ 1 F (c) f + ρ 2 F (1-c) f ( . )
We use the term 'conservative' to express the fact that eq. ( .) computes the temporal evolution of ρ as a difference of fluxes, instead of velocity times its gradient.

Using the above definitions for density fluxes, one can employ any VOF method to construct fluxes of the color function in order to obtain conservative transport for ρ In principle, this should result in the conservation of total mass. However, as we have already pointed out the fact that the CIAM method does not conserve volume exactly , consequently, the advection of the density field (mass) is not consistent with that of the advection of C.

This apparent paradox is resolved by including the compression term on the right hand side of ( . ), in order to ensure consistency . The advection equation for the conserved quantity now becomes -

∂φ ∂t + ∇ • (φ u) φ (∇ • u) ( . )
This equation can be decomposed into advection substeps corresponding to the direction-split integration framework as follows -

The direction-split integration steps operations for φ mirror that of the volume fraction C.

φ n,l+1 i, j,k -φ n,l i, j,k -F (φ) m--F (φ) m+ + φm 1 c (1) m + φm 2 c (2) m ∂ h m u m ( . ) c (1) m
c m is the compression coefficient corresponding to the particular VOF advection method for volume fraction C , while c

(2) m

1c m corresponds to that of the symmetric color fraction 1 -C.

The fluxes F (φ)

m± follow the definition given in ( . ), , while the cell averages φm s are defined as -: This usually corresponds to the tolerance of the Poisson solver, with a typical value being 10 -6 . φm

s ∫ t n+1 t n dt ∫ Ω φH s ∂ h m u m dx ∫ t n+1 t n dt ∫ Ω H s ∂ h m u m dx ( . )
Rewriting the direction-split integration operation for ρ we get -

In case of CIAM, the compression coefficient is defined as c m C n,l i, j,k , whereas in Weymouth-Yue the coefficient is independent of the advection substep l, defined as c m H C n i, j,k -1/2 , where H is the Heaviside function.

ρ n,l+1 i, j,k -ρ n,l i, j,k -F (ρ) m--F (ρ) m+ + C (ρ) m , ( . )
where the fluxes are given by ( . ) and the compression term is Again, there is no implicit summation rule on m.

C (ρ) m ρ 1 c (1) m + ρ 2 c (2) m ∂ h m u m ( . )
For the WY method, the compression terms eventually cancel upon summation over the substeps, therefore resulting in the conservation of mass at the same accuracy as the discrete incompressibility condition

3 m 1 ∂ h m u m 0 is verified .

Advection of Momentum

Within the framework of advection of a generic conserved scalar, we consider momentum advection as transport of the scalar quantities φ ρu q , where q 1, 2, 3 is the component index. Using the definition given in ( . ), we obtain the expression -

The expressions with the 'bar' refer to face weighted averages of the varible in question. Here, φs ρu q s ρu q s ρ s ūq,s ( . )

where ūq,s is termed as the "advected interpolated velocity", whose precise definition is given asūq,s

∫ t n+1 t n dt ∫ f u q u f H s dx ∫ t n+1 t n dt ∫ f u f H s dx ( . )
The subscript s denotes the phase or fluid, and f represents the normal components defined on the face centers.

Thus, the direction-split integration of the momentum can be represented as -

(ρu q ) n,l+1 i, j,k -(ρu q ) n,l i, j,k -F (ρu) m--F (ρu) m+ + ρ 1 ũm q,1 c (1) m + ρ 2 ũm q,2 c (2) m ∂ h m u m ( . )
where the momentum fluxes are constructed in the following manner -

F (ρu) f ρ 1 ūq,1 F (c) f + ρ 2 ūq,2 F (1-c) f , ( . )
The expression for the "central interpolated velocity" corresponding to the averages φm s of ( . ) is -Notice that "cloning" the advected velocities ūq,1 and ūq,2 would make it easier to advect a velocity field with a jump on the interface, but would render the overall algorithm quite complicated. ũm q,s

∫ t n+1 t n dt ∫ Ω u q H s ∂ h m u m dx ∫ t n+1 t n dt ∫ Ω H s ∂ h m u m dx ( . )
The superscript m is intentionally omitted for the velocities ũm q in order to avoid cumbersome and complicated notations.

Due to the viscous effects and the absence of phase change in our fluid dynamics model, the velocity filed maintains continuity across the interface.

As a reasonable approximation we choose to put ūq ūq,1 ūq,2 for the "advected interpolated velocity" and ũq ũq,1 ũq,2 for the "central interpolated velocity". An important simplification which serves as the central model in our development is given by -To add clarity to the notion of "central interpolated velocity", one can understand this as the face-centered interpolations of the velocity field (component normal to control surface), which is required to compute the fluxes of volume.

F (ρu) f ūq F (ρ) f ( . )
Therefore, the advection substep for the momentum can finally we written as -

(ρu q ) n,l+1 i, j,k -(ρu q ) n,l i, j,k -ūq F (ρ) m--ūq F (ρ) m+ + ũq C (ρ) m ( . )
where the density fluxes are defined in ( . ) and the compression term C (ρ) in ( . ).

In the above expression the face-weighted average velocities ūq are defined using ( . ) on the corresponding left face mor right face m+.

Although the weighted averages ūq and ũq have been defined, the manner in which they are estimated shall be covered in the following section.

In the context of the CIAM advection scheme, the compression coefficient for the volume fraction field is C n,l , similarly, for the central interpolated velocity we take ũq u n,l q . Due to the summation of the directional divergences not cancelling out, the resulting transport of the momentum field is not exactly conservative. On the other hand, the compression coefficient is independent of the advection substep in the context of the Weymouth-Yue advection scheme. The final expression for the compression coefficient becomes -

The Weymouth-Yue coefficient is a con- stant value c, that is independent of direc- tion m. C (ρ) m ρ 1 c + ρ 2 (1 -c) ∂ h m u m . ( . )
Since there is no bracketing on any velocity component, we take ũq u n q , which is independent of the substep l.

The compression terms sum up to the accuracy to which the solenoidal nature of the velocity field is discretely verified .

The final expression after cancellation of the compression terms, having undergone three advection substeps ( . ) is : The "advecting" velocity computes the fluxes of volume fraction.

: Figure . has the same variables arrangement that is found in D on a plane perpendicular to the z-axis and through the pressure point p i, j,k .

: In particular , f 1-.

(ρu q ) n,3 i, j,k -(ρu q ) n i, j,k - faces f ūq F (ρ) f . ( . )
Therefore, the extension of the Weymouth-Yue mass (volume) advection scheme to the consistent transport of momentum allows the discrete transport to be exactly conservative.

Velocity Interpolations

The direction-split momentum transport given by ( . ) can be carried out in either of the bulk of the phases, or at close proximity to the interface location. The treatment away from the interface is expressed as -

In the bulk, the expression simplifies considerably due to the density and volume fraction being constant, hence resulting in the cancellation of the spurious compression terms. .

u n,3 q -u n q - faces f ūq u f ( . )
Here, we distinguish between an "advecting" velocity u f u • n f , and an "advected velocity" component ūq, f , which is basically an average over the face f . Due to the momentum transport being carried out on staggered control volumes, both these velocity components require some interpolations from their original positions onto the necessary positions on the staggered grid. Using this new nomenclature, ( . ) becomes -

u n,3 q -u n q - faces f ū(advected) q u (advecting) f ( . )
Now, moving onto the approximation in the neighbourhood of the interface, we get -

In the previous section we had derived a new expression for the momentum fluxes and the compression term, i.e. the RHS of ( . ), that is consistent with the discrete transport of the volume fraction. .

(ρu q ) n,3 -(ρu q ) n - faces f ū(advected) q F (ρ) f + 3 m 1 ũq C (ρ) m ( . )
The crux of our model for momentum advection is described by these two equations (( . ) & ( . )). We start by pointing out the staggered spatial arrangement of our primary variables (presssure, velocity, volume fraction), as illustrated by Fig.

. in a D representation.

To estimate the advecting velocities u

(advecting) f we use a centered scheme.

Consider a face perpendicular to the horizontal direction 1 There are two cases, the first corresponding to when the advected component is not aligned with the face normal, which is the case with q 2. The u 2 control volume in Fig.

. is centered on i, j + 1/2, k, and face f 1-is then centered on i -1/2, j + 1/2, k.

The advecting velocity u (advecting) 1-is not given on this point, therefore it has to be interpolated as -

u (advecting) 1;i-1/2, j+1/2,k 1 2 u 1;i-1/2, j,k + u 1;i-1/2, j+1,k ( . ) face f Ω φ φ Ω Ω face f φ (b) D (a) x q -3/2 -1/2 advected 1/2 φ u = f u > 0 advecting 0

Figure . :

The reference control volume Ω for the advected velocity component φ u q is shown. A horizontal advection is here considered and both the advecting velocity u f and the advected velocity require an interpolation for their value on the left face f 1-: (a) the value ūq φ 0 (full circle) is interpolated from the values φ u q on the nodes (open circles); (b) a more sophisticated interpolation predicts the value φ( x) where x is at the center of the "donating" region Ω D .

: φ 0 serves as an approximation of ūq given in ( . ).

In the second case the advected component is aligned with the face normal, corresponding to q 1. The u 1 control volume in Fig.

. is centered on i + 1/2, j, k and face f 1-is then centered on i, j, k.

In this case, the interpolation becomes -

u (advecting) 1;i, j,k 1 2 u 1;i-1/2, j,k + u 1;i+1/2, j,k ( . )
Thus far, we have covered only interpolations concerning the "advecting" velocity, which is basically used to compute fluxes of the volume fraction. Now we turn our attention to interpolations of the advected velocity ūq which appears in ( . ).

The interpolants we use in this case are one-dimensional and operate on the velocities u q , on the center of their control volumes, that are regularly spaced on a segment aligned with the direction of the advection, that is also perpendicular to face f . In Fig. . , we demonstrate the principle behind our interpolation scheme along the horizontal direction 1, with the lighter notation φ u q . The objective is to interpolate the "advected" velocity on the left face of the reference control volume Ω. In case of advected velocity u 1 and the advecting velocity ( . ) on face f 1-, the relationship with the φ values in Fig. . is -

φ -3/2 u 1;i-3/2, j,k , φ -1/2 u 1;i-1/2, j,k , φ 1/2 u 1;i+1/2, j,k , • • • ( . )
whereas in the case of the "advected" velocity u 2 and the advecting velocity ( . ), the relationship turns out to be -

φ -3/2 u 2;i-2, j+1/2,k , φ -1/2 u 2;i-1, j+1/2,k , φ 1/2 u 2;i, j+1/2,k , • • • ( . )
In order to get a prediction of φ 0 on the face f 1-in Fig. .  we use an interpolation function f that computes this value as a function of the four nearest points, and in an upwind manner based on the sign of the advecting velocity u f . Thus, the interpolation can be expressed in functional form as -

The extension of these approximations to advection along the other two directions q 2, 3 is straightforward.

φ 0 f φ -3/2 , φ -1/2 , φ 1/2 , φ 3/2 , sign(u f ) ( . )
Throughout this body of work, we have extensively tested and used two variants of such interpolation functions :

The reader can refer to influential works of LeVeque [ ] and Sweby [ ] regaring the role of non-linear flux limiters such as QUICK, Superbee, WENO/ENO etc in the context of numerical methods for hyperbolic conservation laws.

: The reconstructions involve geometrical operations, thus are not easily translated into equations.

: These estimations utilize the set of same set of geometric operations that are used in the reconstruction of the interface and fluxed volumes. . : These provisional velocities are required to compute the "advected" velocities using the non-linear flux limiters, as described in ( . ).

. A scheme based on the QUICK third order interpolant in the bulk, away from the interface and a simple first order upwind flux near the interface. . A scheme based on the Superbee slope limiter [ ] [ ]: Roe ( ), 'Some contributions to the modelling of discontinuous flows' for the flux in the bulk and a modified Superbee limiter tuned to a shifted interpolation point near the interface.

. Advection on Staggered Grids

The development of the consistent transport schemes should be integrated into the broader context of our numerical algorithm which deals with the coupling between the conservative formulation of the one-fluid Navier-Stokes equations and the geometric transport of the interface. In this section, we describe two distinct strategies that we have implemented in order to achieve discrete consistency between mass (VOF) and momentum transport, specifically in the context of our staggered Cartesian grid (see Fig. . ).

Shifted Fractions Method

This strategy tackles the problem of staggered control volumes by reconstructing the centered volume fraction field onto the staggered control volumes , thereby resulting in a field of shifted fractions [ ].

A D schematic of this reconstruction procedure is given in Fig. . , considering the advection of the horizontal component of momentum. At the beginning of the operations summarized by L adv , each velocity (momentum) control volume overlaps two pressure/VOF (centered) control volumes i.e. Ω i+1/2, j overlaps Ω i, j and Ω i+1, j Therefore, the two half-fractions from Ω i, j and Ω i+1, j are added together to obtain an estimate of the shifted fraction C i+1/2, j in Ω i+1/2, j . Once this key reconstruction step is carried out, the mass and momentum fields, as well as their corresponding fluxes can all be derived using the information from the shifted volume fraction field. These operations performed at each time step are summarized in the following steps :

. Interface reconstruction at time t n using data C n .

q 1, 2, 3 is the component index, e.g ρ 1 in Ω i+1/2, j,k for the horizontal momentum component ρ 1 u 1 .
. Computation of half-fractions as shown in Fig. . , in order to obtain the "shifted" fraction field C n q and ρ n q in the staggered control volumes. . Computation of all three momentum components (ρ q u q ) n at t n .

The fluxes of momentum are derived from those of mass using the relation ( . ) , which in turn are obtained in a consistent fashion from the fluxes corresponding to 'shifted' volume fraction field.

. Advection of all three momentum components along one coordinate direction, say x direction, using ( . ) to obtain the updated momentum components (ρ q u q ) n,1 after the first substep. . Advection of the "shifted" density field ρ n q on the staggered volumes using VOF consistent advection, along the x direction in order to obtain the updated volume fractions after the first substep. . Extraction of the provisional velocity components u n,1 q after the first substep, i.e. u n,1 q (ρ q u q ) n,1 /ρ n,1 q .

Staggered Centered

Figure . : A D schematic of the reconstruction procedure used in the shifted fractions method. In this case, advection of the horizontal component of momentum is considered. The half-fractions from the centered control volumes Ω i, j and Ω i+1, j (dash-dotted area) are added together to obtain the volume fraction in the shifted (staggered) volume Ω i+1/2, j . Similarly, the half-fractions from Ω i-1, j and Ω i, j (dotted area) are combined to obtain the volume fraction Ω i-1/2, j . Once this horizontally shifted volume fraction field is reconstructed, the corresponding density, momentum and "advected" velocity fields can be computed in order to carry out consistent transport of mass and momentum on this staggered configuration.

: At each time step, the sequence x, y, z is permuted, in order to avoid any systematic biases in error propagation.

. The above operations are repeated for the different momentum components, shifted volume fractions and densities, and velocity components for the next two substeps with split advections along the y and z directions. Finally, we obtain (ρ q u q ) n+1 (ρ q u q ) n,3 and ρn+1 q ρ n,3 q . . In parallel, the centered volume fraction field C n is advected in order to obtain C n+1 C n,3 using the VOF advection method.

An important point to note from the above algorithm is that the advected velocity components u q are updated at each substep, while the advecting velocities u f are interpolated from the initial velocity field u n at time t n . The shifted fractions of Fig.

. are computed using routines that are already part of our geometrical toolbox, used primarily in flux computation.

In essence, we have a consistent transport of density (mass) and momentum on the staggered grids, while simultaneously carrying out advection of the volume fraction field on the centered grid.

A naive attempt to achieve momentum conservation could entail always using the shifted fields C n q and evolve them by the VOF method on the staggered cells. Such an approach would render the scheme conservative, but would result in the independent evolution of the three staggered grids, eventually leading to divergence.

As a consequence, the momentum components (ρ q u q ) n+1 (at time step t n+1 ) that are derived from the corresponding "shifted" volume fractions C n+1 q and densities ρ n+1 q (obtained from C n+1 ) are different from the densities ρn+1 q ρ n,3 q , which were computed during the previous time step by directly advecting the "shifted" volume fractions C n q . The origin of this discrepancy is the approximate nature of the linear reconstruction of the interfaces, which is not even continuous on the boundary of its control volume. Therefore, this implies that the momentum is not exactly conserved between two time steps, even though the Weymouth-Yue method ensures conservative transport of the momentum on the staggered control volumes.

shifted fractions

Figure . : A bird's eye view of the shifted fractions method, highlighting the operations performed at each time step. This method achieves consistency in the discrete transport of mass and momentum on the staggered control volumes. In the interest of brevity, we present a D case for the density on the grid i, j and horizontal velocity u 1 on the staggered grid i + 1/2, j. The evolution of the velocity component u 2 on the staggered grid i, j+1/2 is similar. The initial variables ρ n , u n 1 , u n 2 are inside the ellipses. The interpolated "advecting" components u 1 and u 2 have superscript n. The shifted density ρ n,0 is constructed with the shifted fractions of C to initialize the momentum component (ρu 1 ) n,0 . The first split advection is along the x direction to variables with superscript n, 1, the second one is along the y direction to variables with superscript n, 2. The updated density is ρ n,2 ρ n+1 while the horizontal velocity u n,2 1 u * 1 enters the RHS of the Poisson-like equation . . : The grid is twice as fine when compared to the control volume for momentum.

: In the shifted fractions method, flux computations on the staggered cells necessitate another round of interfacial reconstructions based on the shifted volume fraction field, whereas the sub-grid volume fraction information enables one to circumvent the reconstructions and obtain the fluxes directly using simple surface integrals.

Sub-Grid Method

In this strategy, the difficulty associated with consistent transport on staggered control volumes is resolved by advecting the volume fraction (mass) on a twice finer grid , very much in the spirit of Rudman's [ ] original work.

The primary motivation regarding the development of this method is the issue of conservation (lack of) concerning the discrete transport of momentum in case of the shifted fractions method. Basically, the issue arises quite simply due to the fact that the volume fraction C n+1 used to construct the staggered momentum field at the start of a new time step is not the identical to that which is obtained via advection on the staggered control volumes itself (C n+1 q ). The sub-grid method is able to evade these complications that arise due to the parallel evolution of two different volume fraction fields by conducting volume fraction on solely on the twice refined grid. An additional advantage of the having the sub-grid volume fraction field is that it facilitates a more 'natural' computation of not only the staggered volume fraction field C n+1 q , but more importantly its fluxes. This method is currently only compatible with the implicit fluxes computed in the Weymouth-Yue VOF advection method.

Considering the overall approach, the key differences between our present algorithm and the implementations in the original works of Rudman and Weymouth-Yue are :

Rudman (

)[ ] The use of a conservative geometrical mass (volume fraction) transport framework instead of an algebraic one in the original method, which is applied to sub-grid volume fraction field.

Weymouth & Yue (

) [ ] The extension of the original method for conservative direction-split mass transport to the momentum field, culminating in the discrete consistency between mass momentum transport in the context of staggered Cartesian grids in D.

Spatial Configuration : Coarse & Sub-Grid Variables

Figure . : A D schematic of the arrangement of primary variables on the coarse and sub-grid levels. The control volume is decomposed into the two figures, the one on the left corresponding to the coarse grid variables, while the one on the right corresponding to the sub-grid. Superposition of these two figures, one on top of each other, represent the spatial relationships between the two sets of variables. The velocities on the sub-grid level are simple first or zeroth order interpolations of the coarse grid velocities.

To start, we describe the spatial orientation of the different variables on the two refinement levels of the grid.

The notations for the variables are slightly different that those used in the previous sections (Fig. . .) The velocity component in the horizontal and vertical directions are now denoted as U and V respectively, and the pressure by p.

At the coarse grid level, the pressure and velocity fields are defined in a staggered configuration, with pressures at the cell centers (centroids in D) and velocities on the cell face centers. In Fig. . , we demonstrate the arrangement of variables, which one can easily extrapolate to D during implementation, but for the sake of clarity we illustrate its D equivalent. The volume fraction C is only defined on a grid that is twice as fine than that of the pressure/velocity grid, which we refer to as the sub-grid. Therefore in D, each cubic control volume is divided into eight constituent smaller cubic volumes at the centroids of which, the volume fraction field is centered.

Multiscale Coupling : Restriction & Prolongation Operations

As is paradigmatic for dual-grid methods in the context of Navier-Stokes solvers, we carry out mass transport at the sub-grid level, and momentum transport at the coarse level. In pressure-projection based methods such as ours, the predominant bottleneck in terms of computational speed is the iterative solution to the discrete Poisson equation.

Therefore, solving the mass advection equation on the fine grid allows us not only to obtain more accurate solutions to the flow physics involved, but also enables us to avoid the significantly higher computational costs associated with solving a Poisson problem at a twice finer resolution.

Prolongation In order to advect the volume fraction at the sub-grid level, we need to reconstruct a velocity field at that level of resolution using information from the coarse grid velocity field.

The prolongation operator has to be used at the start of each new time step, in order to compute the sub-grid velocity field.

In order to achieve this, we have implemented a prolongation operator which asserts the following relations between the velocities on the two refinement levels : This is a direct consequence of the solenoidal nature of the coarse grid velocity field.

: The use of higher order interpolation schemes would necessitate finding the solutions of local Poisson problems, in order to render the sub-grid velocity field discretely divergence-free.

u 2i-1 2 ,2 j u 2i-1 2 ,2j+1 U i-1 2 , j ( . ) u 2i+ 3 2 ,2 j u 2i+ 3 2 ,2j+1 U i+ 1 2 , j ( . ) u 2i+ 1 2 ,2 j u 2i+ 1 2 ,2j+1 U i+ 1 2 , j + U i-1 2 , j /2 ( . ) v 2i,2 j-1 2 v 2i+1,2 j-1 2 V i, j+ 1 2 ( . ) v 2i,2 j+ 3 2 v 2i+1,2 j+ 3 2 V i, j-1 2 ( . ) v 2i,2 j+ 1 2 v 2i+1,2 j+ 1 2 V i, j+ 1 2 + V i, j-1 2 /2 ( . )
The notations in the above equations refer to the description of the coarse and sub-grid variables in Fig. . . The first-order interpolation applied to the coarse grid velocity field ensures discrete incompressibility at the sub-grid level. The choice of interpolation order used in the prolongation operator is identical to that of the originial method of Rudman ([ ]). In the current version of the sub-grid method, we have chosen to eschew the added complexity of solving local Poisson problems through our prolongation operators, therefore sticking with the much simpler Rudman inspired interpolations. Restriction In order to compute the required density and momentum fields on the staggered grid (coarse level), we implement restriction operators that use information from the sub-grid volume fraction field. For sake of clarity, present a D case for the density on the grid i, j and horizontal velocity ũ on the staggered grid i + 1/2, j. The variables in red are all defined on the coarse level, and those in green are at the sub-grid level. 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 convery that the variables are centered on the staggered control volumes at the coarse grid level, whereas the 'bar' on top conveys that the variables are staggered with respect to the sub-grid. The starting variables are enclosed in elliptical borders, and the variables at the end of the directionsplit integrations are u * and ρ n+1 , where u * is subsequently fed into the RHS of the Poisson problem.

Fine Grid -→ Coarse Centered Grid

Density & Viscosity The sub-grid volume fraction field is restricted to compute a centered volume fraction field at the coarse grid level, which is subsequently used to derive the centered density and viscosity fields at the coarse grid level.

Algorithm

To summarize the operations performed at each time step pertaining to the advection operator in our one-fluid Navier-Stokes framework, we present Fig. . , which illustrates a D version of the sub-grid method that ensures consistent and conservative mass-momentum transport, highlighting the interactions between the variables defined on the different grids. The algorithm is also summarized in the following steps :

. Interface reconstruction at time t n , at the sub-grid level, using data c n .

: These provisional velocities are required to compute the "advected" velocities using the non-linear flux limiters, as described in ( . ).

: At each time step, the sequence x, y is permuted, in order to avoid any systematic biases in error propagation. q 1, 2 is the component index, e.g ρ 1 in Ω i+1/2, j for the horizontal momentum component ρ 1 u 1 .

. Restriction of the sub-grid volume fraction in order to compute the "shifted" fraction fields C n q and ρ n q in the staggered control volumes.

. Computation of all momentum components (ρ q u q ) n at t n .

The fluxes of momentum are derived from those of mass using the relation ( . ) , which in turn are obtained via restriction operations applied to the sub-grid volume fraction fluxes.

. Advection of the sub-grid volume fraction c n along one coordinate direction, say x direction, to obtain c n+ 1 2 , using the Weymouth-Yue explicit advection method.

. Advection of all momentum components along the x direction, in sync with the sub-grid VOF advection, using ( . ) to obtain the updated momentum components (ρ q u q ) n+ 1 2 after the first substep. . Restriction of the sub-grid field c n+ (ρ q u q ) n+ 1 2 /ρ n+ 1 2 q . . The above operations are repeated for the different momentum components, in sync with the direction-split advection of the subgrid volume fraction along the remaining direction. Finally, we obtain u * q (ρ q u q ) n+1 /ρ n+1 q on the staggered grid. . In parallel, the updated centered volume fraction field C n+1 is obtained at the coase grid level, by applying the restriction operator on the updated sub-grid volume fraction c n+1 .

The resulting method not only maintains discrete consistency between mass and momentum transport on the staggered grid, but also ensures that the transport of momentum is consistent in D. The other time-split terms in equation ( . ) that arise from the source term discretizations, as well as the terms in the projection step ( . ) are solved in a standard non-conservative way, and are briefly covered in what follows. The density on the faces of the central cells are estimated using the restriction operations

Generally, the face-centered density field is estimated using a simple average i.e. ρ i+1/2, j,k (ρ i, j,k + ρ i+1, j,k )/2.

described previously, with these densities appearing via the 1/ρ q pre-factor in front of all the terms in the momentum transport (e.g surface tension, viscous diffusion) that are discretized in a nonconservative manner. Additionally, the sub-grid volume fraction can be used to get better estimates of the curvature field, which itself is defined on the staggered control volumes. At the time of writing, two different approaches are being developed and tested in order to get more accurate curvature fields. The detailed descriptions of the methods used in our numerical platforms to deal with surface tension, viscosity and body forces have already been carried out in [ , , ], therefore we briefly touch upon certain aspects of the operators in question, in particular, their interaction with the volume fraction field.

Surface Tension

We use the Continuum Surface Force method (CSF) as our model for surface tension, coupled with height functions for curvature computation. The height functions used in our implementation were first introduced in [ ] , subsequently tested, revised and improved in [ , ]. In general, the height functions are used to compute the curvature field based on secondorder finite differences applied to the heights. Although, in regions of : These are regions where the local radius of curvature is comparable to the grid size, in fact for certain cases the height functions perform poorly even for regions resolved by -grid points.

: The exact implementation differs slightly between 'Basilisk' and 'PARIS Simulator'.

poor interfacial resolution , the method reverts to certain fallbacks, one of which is curve fitting instead of height functions. The resulting curvature field is coupled with a well-balanced discretization with respect to the discrete pressure gradient, with the same discretization stencil applied to the volumetric (body) force term as well.

Viscous Diffusion

We use second-order spatial discretizations of the viscous stresses, using centered differences . The (variable) dynamic viscosity is computed based on the volume fraction field via weighted arithmetic or harmonic averaging (equations . ). The temporal treatment of the viscous term can be either in explicit or semi-implicit fashion, , but in the context of the present study we will be sticking exclusively with the explicit version.

Pressure-Poisson Projection

The velocity field is evolved using a classical time-splitting projection method as described in the seminal work of Chorin [ ], which involves predicting an 'intermediate' velocity field u * as given by equation . , followed by a correction step as follows

u n+1 u * - τ ρ n+1 ∇ h p n+1 ( . )
The discrete pressure field required to correct the intermediate velocity is determined by imposing the conservation of mass, which in our incompressible framework reduces to necessitating the resulting velocity field to be divergence-free (solenoidal)

∇ h • u n+1 0 ( . )
Thus, combining equations . and . , we are left with a variable coefficient Poisson equation for the pressure :

∇ h • τ ρ n+1 ∇ h p n+1 ∇ h • u * ( . )
The Poisson solver used in 'PARIS Simulator' to invert the elliptic operator appearing in eqn. . is a red-black Gauss-Seidel (GS) solver with overrelaxation [ ] [ ]: Briggs ( ), 'A Multigrid Tutorial, SIAM'

. There is also has an in-house implementation of a multigrid solver for structured grids with 2 n number of points per direction, utilizing a fully parallelized V-Cycle scheme [ ]. Relaxation operations are applied starting from the finest to the coarsest first, and then from the coarsest to the finest, with the number of relaxation operations being a user-adjustable parameter. Having a native multigrid solver allows for an efficient solution of the Poisson equation without the necessity of having external libraries/pre-conditioners (e.g. HYPRE) installed on the system. When it comes to 'Basilisk', an atypical multigrid solver is implemented using a "half" V-cycle in order to deal with the spatial inhomogeity of the grid size arising due to adaptive mesh refinement. For more details regarding the differences between the multigrid solver of 'Basilisk' and the classical implementation of a multigrid, one can refer to [ ] [ ]: Popinet ( ), 'Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries' .

The whole set of operations described up to this point, constitutes a temporal integration scheme of the first-order, which can be expressed as

C n+1 , u n+1 L 1 (C n , u n ) ( . )
where L 1 is the operator consisting of all the steps described so far, applied to the primary fields C and u. Therefore, a second-order time integration can easily be computed by using L 1 to get a first prediction

(C * * , u * * ) L 1 (C n , u n ) ( . )
The superscript * * refers to our first order prediction of the primary variables. Therefore, the second-order estimate can be obtained via averaging

C n+1 , u n+1 1 2 [(C * * , u * * ) + L 1 (C * * , u * * )] ( . )
In the following chapter, we shall take a look at some performance aspects of these mass-momentum consistent methods, using both quantitative and qualitative comparisons. In the upcoming sections, we demonstrate the robustness and accuracy of our class of mass-momentum consistent numerical methods when applied to challenging flow configurations involving marked density contrasts, primarily in comparison with the version of our method which does not maintain consistency between the mass and momentum advection. Most of the standard tests that exist in the current literature concerning numerical methods to tackle liquid-gas flows such as the decay of spurious currents in static and moving droplets, viscous damping of capillary waves etc., are carried out in the absence of any density jump (or viscosity jump) across the interface separating the fluids. In this chapter, we shall take a closer look in detail at the behavior of our methods when dealing with difficulties that arise due to the non-linear coupling between interfacial deformation/propagation, capillary and viscous forces, especially in the regime where the material properties across the interface are separated by orders of magnitude, particularly in which the flow features in question are poorly resolved.

Numerical

In order to assess the performance of the different methods, we shall use an easier nomenclature to describe the different methods, which are as follows :

STD : The standard method which does not maintain consistency between mass and momentum transport. MSHIFT : The method that ensures consistency between mass and momentum transport, but is not exactly momentum conservative.

It is based on the shifted fractions strategy, as detailed in the previous chapter. MSUB : The method that ensures both consistency between mass and momentum transport, as well as ensuring exactly conservative transport. It uses the sub-grid strategy as discussed in the preceding chapter.

. Static Droplet

A popular numerical benchmark in the existing literature relevant to surface tension dominated flows is the case of a spherical droplet of the denser fluid immersed in a quiescent surrounding medium of the lighter fluid. In the hydrostatic limit of the Navier-Stokes equations, the droplet should stay in equilibrium, with a curvature induced pressure jump across the interface corresponding to Laplace's equilibrium. In practice however, numerically reproducing such a trivial equilibrium condition is not as straightforward, as there exists a slight difference between the initial numerical interface and the exact analytical shape of the sphere, thereby resulting in the generation of the well documented 'spurious' or 'parasitic' currents of varying intensity in the velocity field . A lot or progress has been made since in the context of wellbalanced surface tension formulations, that ensure consistency between the numerical stencils used for the discretization of the pressure gradient and the Heaviside approximation (nδ s ) that projects the the surface force distribution onto the control volumes [ , ]. A significant contribution to the interpretation of these parasitic currents within the well-balanced framework was made by Popinet [ ] which demonstrated that given sufficient time (of the order of viscous dissipation time-scales), a wellbalanced method will relax to the 'numerical' equilibrium shape through the damping of the 'physically consistent' numerical capillary waves, therefore allowing us to recover the exact (to machine precision) Laplace equilibrium condition. The key difference in our implementation of this classic test case from that of Popinet [ ] [ ]: Popinet ( ), 'An accurate adaptive solver for surface-tension-driven interfacial flows' is that we consider the effect of density contrast across the interface separating the fluids. As we have previously discussed, a sharp density jump across the interface may have an amplification effect on the numerical errors incurred as a result of interfacial reconstructions, curvature estimation and various other truncations, thereby rendering the method unstable. We demonstrate that in our framework of mass consistent momentum transport coupled with a well-balanced surface tension discretization, density-ratios as large as 1000 : 1 can be simulated : The viscous time-scale corresponding to the droplet length-scale is the most commonly used in literature.

Setup

: Similar criteria are defined on the basis of the viscous and advection operators as well, with the smallest amongst the three selecting the numerical time-step without loss of numerical stability, in conjuction with the ability to recover the exact numerical equilibrium through the dissipation of spurious currents within relevant time-scales .

We consider a circular droplet of size D placed at the centre of a square domain 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 surface tension coefficient (fig. . ). The ratio of the droplet size to the box is chosen as D/L 0.4, coupled with a numerical resolution of D/∆x 16 (where ∆x is the grid size). As for boundary conditions, we use symmetry conditions on all sides of the square domain.

The problem incorporates two natural time-scales, the capillary oscillation scale and the viscous dissipation scale, which are defined below :

T σ ρ l D 3 σ 1/2 , T µ ρ l D 2 µ l ( . )
The ratio of these time-scales give us -

T µ T σ ρ l σD/µ l √ La ( . )
where La is the Laplace number based upon the heavier fluid. In the present study, we introduce the density-ratio ρ l /ρ g as another important parameter. In order to rescale our 'parasitic' velocity field, we define a velocity scale based on capillary oscillations as -

U σ σ/ρ l D ( . )
Additionally, the time-step in our numerical simulation must be smaller than the oscillation period corresponding to the grid wavenumber (fastest capillary wave with a time period ∼ ρ l ∆x 3 /σ

1/2
) as a stability criterion , as our surface tension model is explicit in time. For the scope of the present study, we shall not consider any viscosity contrast between the two fluids while varying the density-ratio, therefore µ l /µ g 1 for all the cases under study.

Decay of Spurious Currents

In figures . to . , we illustrate the decay of the root-mean-square of the spurious currents as a function of time, in the case of four different densityratios, with three different Laplace numbers for each ratio. The first figure ( . ) refers to simulations carried out without consistency between the momentum-mass transport (STD), the second ( . ) corresponds to that of the consistent but not conservative method (MSHIFT), and final one ( . ) refers to that of the consistent and conservative method (MSUB). The time is rescaled by the viscous dissipation scale, and the spurious currents by the capillary velocity scale. We have two main observations, the rapid decay of the rescaled spurious currents for all combinations of density-ratios and Laplace numbers within approximately 0.2T µ , and the slower re-growth of the currents in question for combinations of nonunity density-ratios and large Laplace numbers, in all simulations except those carried out with MSUB. With method MSUB, the decayed currents keep hovering around levels of machine precision for remainder of time. Although there is a re-growth of the currents using the consistent method (MSHIFT) after 0.2T µ , the behavior is not quite alarming as the rate of this re-growth is quite low. Therefore, out of all the methods tested, the consistent and conservative method (MSUB) does seem to demonstrate the desired performance, especially when it comes to combinations of large density contrasts coupled with large Laplace numbers. The currents seem to initially decay quickly for all higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For combinations of large ρ l /ρ g and large La, the spurious currents seem to grow back to an order of magnitude (10 -8 ) which is quite far from that of machine precision (10 -14 ).

Spatial Convergence

Once the solution relaxes to a numerical equilibrium curvature (spurious currents are approximately at the order of machine precision), there still exists a difference between the numerical curvature and the exact analytical curvature corresponding to the spherical (circular) shape. We use the definitions of the shape errors as introduced in the seminal work of Popinet [ ] to assess the convergence of our class of methods to the exact (analytical) curvature as we increase spatial resolution. The norms are defined as follows : The currents seem to initially decay quickly for all higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For combinations of large ρ l /ρ g and large La, the spurious currents seem to grow back to an order of magnitude (10 -8 ) which is quite far from that of machine precision (10 -14 ). No considerable improvement is observed with respect to STD.

Figure . : MSUB Decay of normalized spurious currents as a function of viscous dissipation time-scales for different

density-ratios and Laplace numbers. The currents seem to decay very quickly in the case of higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For all combinations of ρ l /ρ g and La numbers, the decayed spurious currents are not observed to grow back as in the cases of STD and MSHIFT, and hover around values close to machine precision (10 -14 ).

L 2 i C i -C exact i 2 i , L ∞ max i |C i -C exact i | ( . )
where C i is the volume fraction of a cell after the solution has relaxed to the numerical equilibrium curvature, and C exact i is the volume fraction corresponding to the exact circular shape which was initialized at the start of the simulation. Fig. . demonstrates the behavior of the shape errors defined in eqn. . for the case of the most stringent parameter combination ( ρ l /ρ g 1000 , La 12000 ) as a function of the droplet resolution. As one can clearly observe, all the methods tested display a roughly second-order convergence in space for both the error norms. In terms of the L 2 norm, the consistent and conservative method (MSUB) does indeed achieve smaller errors as compared to both STD and MSHIFT for all spatial resolutions. As a minor remark, there is not much to discern in terms of shape error when it comes to comparing the performances of the consistent (MSHIFT) method with the non-consistent one (STD). Shape Error

L ∞ : MSUB L ∞ : MSHIFT L ∞ : STD L 2 : MSUB L 2 : MSHIFT L 2 : STD Second Order
Figure . : Second-order spatial convergence for the spurious current error norms corresponding to the most stringent parameter combination (ρ l /ρ g 1000 , La 12000) . Both of the norms (L ∞ and L 2 ) seem to demonstrate a roughly second order rate of spatial convergence with each of the methods tested. However, MSUB has a marginally lower L 2 error compared to both STD and MSHIFT for all resolutions tested. There is negligible difference observed in the shape errors between STD and MSHIFT in both of the norm definitions.

. Moving Droplet

An incisive numerical setup that enables us to evaluate the accuracy of the coupling between interfacial propagation and surface tension discretization was first proposed by Popinet [ ] . The manner in which this test differs from that of the static droplet is the addition of a uniform background velocity field, therefore serving as a better representation of droplets in complex surface tension dominated flows where they might be advected by the mean flow. In terms of the Laplace equilibrium, the hydrostatic solution is still valid in the frame of reference of the moving droplet. The point at which the solution in the moving reference frame diverges from that of the static droplet ( . ) is through the continuous injection of noise at the scale of the grid size. This 'numerical' : In Popinet [ ], a resolution of D/∆x 25.6 corresponding to a grid of 64 × 64 is used noise emanates from the perturbations to the curvature estimates, which are in turn induced by the interfacial reconstructions carried out to propagate the interface (temporal integration) . These fluctuating errors act as source terms for the momentum, thereby transforming the problem into that of viscous dissipation in the presence of continuous forcing (in the reference frame of the moving drop). In the present study, we evaluate our class of methods using the advection of a droplet in a spatially periodic domain using an identical setup as [ ], but with the important difference of including sharp density jumps across the interface as well as using lower spatial resolutions. As previously discussed ( . ), large density contrasts tend to amplify the fluctuations induced by the myriad numerical approximations (interface reconstruction, curvature estimation etc) involved in the algorithm.

Setup

We consider a circular droplet of diameter D placed at the centre of a square domain 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 surface tension coefficient (fig. . ). A uniform velocity field U is initialized on the entire domain (only a horizontal component). The ratio of the droplet size to the box is D/L 0.4, with D/∆x 20 (∆x being the grid size. ). As for boundary conditions, we use symmetry conditions on the top and bottom sides, and periodic boundary conditions on the horizontal direction (along which advection by U takes place). We characterize by problem by introducing the following adimensional parameters (based on the heavier fluid) :

La ρ l σD µ 2 l , We ρ l U 2 D σ ( . )
In addition to the capillary and viscous time-scales for the static case (eqns. . ), we have an additional scale defined as :

T u D/U ( . )
which is the time-scale of advection. In our subsequent analysis, we shall use T u and U as the time and velocity scales, repectively. The currents seem to hover around 10 -3 , with a larger Laplace number corresponding to a higher error for all density-ratios. We 0.4 for all the cases presented.

Evolution of Spurious Currents

Figures . to . : MSHIFT Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers. There seems to be no appreciable difference from the evolution seen in the case of STD (fig. . ). The currents seem to hover around 10 -3 , with a larger Laplace number corresponding to a higher error for all density-ratios. We 0.4 for all the cases presented.

carried out without consistency between the momentum-mass transport (STD), the second ( . ) corresponds to that of the consistent but not conservative method (MSHIFT), and final one ( . ) refers to that of the consistent and conservative method (MSUB). We again have a couple of important observations, the first being that spurious currents do not decay to machine precision as in the static droplet case for all of the combinations and methods tested, instead they oscillate around a mean value of the order of (0.1 -0.01)% of the constant field U. The second observation is regarding the significantly smaller error (almost by one order of magnitude) in the case of the consistent and conservative method (MSUB) when compared to that of STD and MSHIFT. As a minor remark, in case of large Laplace numbers, the MSUB method displays a slight upward trend in the error evolution, which is not the case in either STD or MSHIFT. This is not too worrisome as the growth is over a time-scale much larger than T u , with the oscillations corresponding to a time-scale of the order U/∆x. All of the plots in figures . to . correspond to We 0.4, alongside an additional simplification of equal viscosities across the interface i.e µ l /µ g 1 .

As evindenced by the persistence of these spurious currents due to the addition of grid-level noise emanating from interfacial reconstructions, further advancements should be made with respect to the combined performace of the interfacial transport, curvature computation and the surface tension model. Nonetheless, all the methods tested do seem to be quite numerically stable when dealing with the high density-ratios, and are not subject to rapid uncontrollable amplifications of the interfacial Although an upward trend is observed for large Laplace numbers, the growth rate is quite low. The currents seem to hover slightly above 10 -4 , with larger Laplace numbers corresponding to larger errors for all density-ratios. We 0.4 for all the cases presented.

: In existing literature, convergence rates have only been studied in case of equal density fluids across the interface perturbations even for high Laplace numbers.

Spatial Convergence

In order to evaluate the performance of our class of methods at different resolutions, we define the errors as the maximum values of the norms L ∞ and L 2 of the rescaled field U rms /U over time (5 times T u ). In fig.

. , we show the scaling of the error as a function of spatial resolution for the most stringent case of ρ l /ρ g 1000 , La 12000, for each of our different methods. As similarly observed in section . , in terms of both L ∞ and L 2 norms, there is no appreciable difference in the behaviors of STD and MSHIFT. For MSUB, we do observe significantly lower maximum errors compared to other two methods, but at a cost of slightly less than first-order convergence. The overall convergence behavior of the class of methods we have tested seems to be consistent with earlier studies of Popinet [ ] and others .

Error Dependence : Laplace & Weber numbers

As the final point of inquiry into the performance of our class of methods, figures . and . demonstrate the influence of the Laplace and Weber numbers on the behavior of the maximum error norm, carried out for the largest density-ratio (ρ l /ρ g 1000). We present results obtained using the consistent and conservative method (MSUB), for a resolution corresponding to D/∆x 25.6. As we can observe, the error (both L ∞ and L 2 ) scales as We -1/3 over orders of magnitude, which is different from the We -1/2 scaling observed by Popinet [ ] . In terms of Laplace numbers, the errors scale as La 1/6 over two orders of magnitude, which is the same as that observed in [ ] (for equal densities). .

Capillary Wave

One of the fundamental features of immiscible multiphase flows involving interfaces are the presense and propagation of capillary waves. Therefore, a robust and accurate numerical method should not only be able to adequately resolve, but also accurately emulate the spatio-temporal evolution of such surface tension induced oscillations. A brief outline on the state-of-the-art numerical implementations of capillary waves In the present study, we evaluate the accuracy of our class of methods by comparing with an analytical solution of damped capillary oscillations. Generally, analytical solutions exist only for cases corresponding to extremely small initial perturbations, that too either in the inviscid limit (Lamb [ ] We consider a rectangular domain of dimensions L × 3L, where L corresponds to the wavelength of our initial perturbation. The densities of the heavier and lighter phases are ρ l and ρ g respectively, likewise for the viscosities µ l and µ g , and σ being the surface tension coefficient (fig. . ). An intial perturbation amplitude of L/100 is used, coupled with a numerical resolution given by L/∆x 32 (∆x being the grid size). Symmetry conditions are applied on the top and bottom sides, with periodic conditions along the horizontal direction. We use the following adimensional parameters to characterize our problem :

T 0 Tω 0 , La ρ l σL µ 2 l ( . )
where La is the Laplace number based on the heavier fluid, and ω 0 is defined using the dispersion relation used in Popinet [ ] given as :

ω 2 0 σk 3 2ρ l , where k 2π L ( . )
The dispersion relation is obtained via linear stability analysis at the inviscid limit [ ] [ ]: Lamb ( ), Hydrodynamics

. In order to evaluate the influence of density-ratio on the performance of our class of methods, we use three different numerical setups keeping the same Laplace number (La 3000) as follows :

ρ l /ρ g 1 , µ l /µ g 1 (Popinet [ ])
ρ l /ρ g 10 , µ l /µ g 1 ρ l /ρ g 1000.0/1.2 , µ l /µ g 1.003 • 10 -3 /1.8 • 10 -5 (Air-Water)

The final setup corresponds to that of an air-water interface (physical properties corresponding to 20 o Celsuis), which is the most stringent due to the significant density and viscosity jumps.

Comparison with Prosperetti Solution

The theoretical solution to this configuration corresponds to the closedform expressions of the planar interface shape evolution established by Prosperetti [ , ], which takes into account the finite time-scales at which the vorticity (generated due to interface oscillations) diffuses into the bulk medium. These closed-form expressions are subsequently integrated using a fourth-order Runge-Kutta time integrator (details of which not described here), and used to assess the accuracy of the results obtained by our class of numerical methods.

As we can in figures . to . , solutions from our class of numerical methods (circles) are compared to that of the theoretical (Prosperetti) solution (black curves), where the amplitude is normalized by the initial value (A 0 ) and the time rescaled by T 0 . The first figure ( . ) refers to simulations carried out by the non-consistent method, the second ( . ) corresponds to that of the consistent method, and the final one ( . ) refers to that of the consistent and conservative method. We observe that there is hardly any appreciable qualitative difference between the results obtained via the different methods STD, MSHIFT and MSUB, although MSUB does seem to perforn marginally better when it comes to the most stringent case (air-water configuration). Surprisingly, even the non-consistent method (STD) does not seem to show any un-physical interfacial deformations for all the density-ratios tested, and that it is difficult to distinguish between the different methods for the lower density-ratios.

Spatial Convergence

The next step in our evaluation would be to quantify the accuracy of our numerical results to the Prosperetti solution using an integral (in time) error norm, the same as defined in [ ] :

L 2 1 L ω 0 25 ∫ T t 0 (h -h exact ) 2 ( . )
where h is the maximum inteface height obtained using our numerical simulations, and h e xact being the maximum height obtained via time integration of the Prosperetti solution. In figures . to . we demonstrate the rate of spatial convergence of the L 2 error norms for different density-ratios, simultaneously comparing the behavior of the different methods STD, MSHIFT and MSUB at each density-ratio. In all the results presented, we maintain La 3000 for all density-ratios, spatial resolutions and methods tested.

In figure . we observe roughly second-order spatial convergence when it comes to equal densities across the interface, with STD and MSHIFT displaying nearly identical behavior, whereas MSUB does marginally better with lower errors for all resolutions. When it comes to ρ l /ρ g 10 in figure . , we observe a saturation in the initial second-order convergence rate irrespective of whichever method is used, however MSUB performs much better in terms of error when compared STD and MSHIFT. Finally, figure . demonstrates the roughly second-order convergence of all three methods when it comes to the air-water configuration, again, with MSUB performing marginally better with lower errors. Not surprisingly, the largest errors arise for the air-water configuration errors across all methods.

Second Order

Figure . : Comparison of spatial convergence for the case of ρ l /ρ g 1, La 3000, for our class of methods. There is no viscosity jump across the interface. All methods seem to demonstrate approximately second-order convergence. There seems to be no appreciable difference in the behavior of STD and MSHIFT, with MSUB displaying marginally lower errors compared to the others.

2 3 2 4 2 5 2 6
L/∆x The focus of the present chapter shall be to take a closer look at the issues plaguing numerical methods that try to deal with flows entailing significant density contrasts between the fluids. Arguably, the most common instances of such flows are those involving the interaction between air and water, primarily in the form of droplets and bubbles that play important roles in both natural and industrial processes. As for the numerical platform, we use 'PARIS Simulator' to carry out our simulations on this topic. A flow configuration that combines the complexities of large density-ratios with the interaction between capillary, viscous and inertial stresses is that of a water droplet falling in air under the influence of gravitational acceleration.

. Problem Setup

The problem is characterized by a combination of Reynolds, Weber and Bond numbers, the definitions of which are as follows :

We

ρ g U 2 D σ , Re ρ l UD µ g , Bo ρ l -ρ g gD 2 σ ( . )
As the density and viscosity ratios are corresponding to that of air-water systems at degree Celcius, an equivalent characterization could be done using the Morton number Mo

gµ 4 g (ρ l -ρ g ) ρ 2 g σ 3
in place of the Bond number.

The subscripts l and g represent liquid and gas phases respectively. In our particular numerical setup, We 3.2, Re 1455 and Bo 1.2, thus corresponding to that of a 3mm diameter raindrop (a relatively large one) falling in the air at an approximate terminal velocity of 8 m/s (interpolated from empirical data, refer to Gunn and Kinzer [ ]). This choice of length scale of the drop is motivated by the paradigmatic value of a near-spherical raindrop simulation, and by the fact that the corresponding Weber number (∼ 3) is the same as in a similar air-water setup corresponding to suddenly-accelerated-droplet (or "secondary atomization" ) simulations in the studies [ , ]. For such a low Weber number the capillary forces dominate and the droplet should remain intact, and definitely not undergo any subsequent atomization. The parameters in the problem setup are given in Table . , and the schematic diagram given in Fig. . . The droplet is initially placed at the center of a cubic domain ( D), where the length of the side is times the diameter of the drop.

ρ g ρ l µ g µ l σ D g k g/m 3 k g/m 3 (Pa s) (Pa s) (N/m) (m) (m/s 2 )
. 0.9982 × 10 3 1.98 × 10 -5 8.9 × 10 -4 0.0728 3 × 10 -3 9.81

: Even for simulations with points across the droplet diameter, the boundary layer has only by -cells across it. : In the Dodd and Ferrante [ ] study , a droplet was allowed to fall from rest along a domain with a length corresponding to times the droplet diameter, therefore necessitating a problem size of approximately million cells.

: Droplet resolutions between D/h 16 to D/h 64 are considered to be moderate, where D is the droplet diameter and h is the grid size.

In order to properly reproduce and analyse the dynamics of a relatively large drop (high Reynolds flow) such as in our case, the numerical method has to accurately resolve the thin boundary layers , the interaction of such layers with the capillary forces and finally the non-linear feedback of the complex D vortical structures present in the wake behind the droplet. Arguably, the most natural type of computational setup would involve using a large domain filled with air at rest, with zero inflow velocity and to let the droplet fall from the top of the domain. Such an undertaking was attempted by Dodd and Ferrante [ ], in which they managed to delineate the different regimes concerning the behavior of the wake behind the droplet, although at relatively lower Reynolds numbers corresponding to smaller drops (the maximum Reynolds tested was 500, whereas in our case it is 1500). We on the other hand use a significantly smaller domain (L/D 4 where L is the domain size), with a constant value of inflow velocity (close to 8 m/s), thus resulting in the drop exiting the domain after a certain amount of time. This setup proves to be quite convenient for relatively short-time investigations. Therefore, our objective behind the demonstration of this particular test case is not to develop a high fidelity model of a raindrop , but rather carry out a stringent evaluation of the robustness of our class of numerical methods (shifted fractions and sub-grid strategies) compared to the standard version of the method that is not mass-momentum consistent. 4. The liquid properties (ρ l , µ l ) correspond to that of water, and the gas properties (ρ g ,µ g ) correspond to that of air. We apply a uniform inflow velocity condition with U 0 (t) and an outflow velocity condition at the top which corresponds to zero normal gradient. Boundary conditions on the side walls correspond to those of impenetrable free slip (no shear stress).

. Numerical Instabilities

Numerical simulations using a fixed inflow velocity with U 0 (t) 8 m/s were carried out for very short times (of the order of ms) at moderate resolution . The simulations carried out with the standard version of our numerical method result in the catastrophic deformations of the droplet as illustrated in Fig. . , which we describe as "fictitious" or "artificial" atomization. They display marked peaks or spikes in kinetic energy as a function of time, associated with massively deformed interface shapes (see figure . ). Additionally, our studies suggest that certain combinations of the advection scheme and the flux limiter are numerically more robust than others ,in particular, the most stable combinations are that of CIAM advection with Superbee limiter, and the WY advection with QUICK limiter.

: As the simulation starts with uniform flow in the gas and a zero velocity in the drop, there is a sudden large acceleration in the gas, resulting in the development of a dipolar velocity field akin to that of potential flow around a cylinder (sphere). : In the thesis of Xiao, the focus was on the analysis of droplet breakup corresponding to different Weber number regimes.

Instability Mechanism

We propose the following explanation in order to account for such numerical artifacts. To start with, we neglect gravity and viscous effects at this relatively large Reynolds number. Also, we are interested in steadystate flow . On the axis and near the hyperbolic stagnation point at the front of the droplet one has u 2 0 for the transverse (radial) velocity and for the axial momentum balance

u 1 ∂ 1 u 1 - 1 ρ ∂ 1 p. ( . )
Due to the large viscosity and density ratios, it is not possible for the air flow to immediately entrain the water, so the fluid velocity is significantly smaller inside the droplet. In the air the acceleration near the stagnation point is of the order U 2 /D, whereas the pressure gradient is

∂ 1 p ∼ ρ g U 2 /D. ( . )
The pressure gradient in the liquid is much smaller, however, in the case of a mixed cell the water density may multiply (due to numerical errors) with the gas acceleration U 2 /D, so that

∂ 1 p ∼ ρ l U 2 /D, ( . )
then a large pressure gradient results in the mixed cell or cells. This large pressure gradient results in pressure spikes inside the droplet near the front stagnation point, as shown in Figure . . Such pressure gradients are balanced by surface tension only for a sufficiently large curvatures near the droplet front. This explains the presence of a "dimple" often observed in low resolution simulations of the falling drop. This artifact had been observed by Xiao [ ] in a similar case involving the sudden interaction of a droplet at rest with a uniform gas flow. The resulting large un-physical pressure gradients across the interface eventually lead to its rapid destabilization and concomitant breakup.

Visualization of the flow around the droplets (Figure . ) illustrates the challenging nature of the flow configuration, even for such a seemingly simple physical problem. As one can observe, the thin boundary layers are poorly resolved. Therefore, even though the velocity field is continuous across the interface (in the discrete sense) in the absence of mass transfer, there is the appearance of strong velocity variations at the scale of the grid size for such coarse levels of grid refinement. : A prerequisite of methods that use consistent mass-momentum transport is the formulation of the convective operators in the conservative manner i.e divergence of fluxes instead of gradients of the primary (discontinuous) variables.

. Stabilization Strategies

The cascading nature of the numerical instabilities that lead to the eventual (un-physical) fragmentation of a droplet is not just a cause of concern in the context of modeling a raindrop, but it has important implications within the more broader scope of atomization and other complex fragmentation phenomena. Arguably, the most important objective of numerical investigations of physical phenomena involving fragmentation is the quantification of the size of the features that result from a series of topological changes A typical example is the statistical distributions of drop sizes. As demonstrated in the present study, using the standard non-consistent method one observes artificially atomized drops thus leading to the generation of a large number of smaller fragments, which if counted, will skew the resulting droplet size distributions towards the smaller sizes. Therefore, we assert that the falling raindrop case serves as a faithful representation of the numerous under-resolved features that are typically omnipresent in numerical studies of complex fragmentation phenomena. The inability to discern between the drops resulting from physically consistent mechanisms and those that result from numerical fragmentation is exactly why there is a need to develop numerical schemes that circumvent the occurrence of such instabilities, especially as there will always be some constraints in computational power.

Consistent Mass-Momentum Transport

The most common and brute force approach that one can apply in order to suppress or circumvent such numerical instabilities is by using a combination of extremely refined meshes coupled with large domains [ ]. A more computationally efficient approach might be to use a conservative formulation of the Navier-Stokes equations, in order to achieve consistency in the discrete transport of mass and momentum. A thorough and detailed exposition of the principles behind the consistent methods and the different strategies towards achieving consistent transport has already been covered in the chapter . We observe that the use of our class of consistent methods enables us to stabilize a majority of the numerical instabilities and bring a systematic improvement over wide range of flux limiters (WENO, ENO, Superbee, QUICK, Verstappen) and CFL numbers, as evidenced by comparing the figures . and . .

Convergence Study

For the next set of simulations, we use a fixed inflow velocity setup but with smaller initial velocity. We systematically vary the resolution from D/h 8, 16, 32 and 64. Despite using our consistent methods, simulations at D/h 8 are sometimes unstable, so we use a workaround and use a lower fixed inflow velocity of U 0 5 m/s, which differs from the expected long term terminal velocity U t m/s used in the previous set of simulations. But, it does offer a milder initial condition and allows us to observe the first phase of the (physical) acceleration towards the final statistical steady state. The simulations are carried out for ms (for reasons of computational cost and also in order to avoid the droplet getting too close to the domain boundaries) and the convergence : A juxtaposition of the different manifestations of 'artificial' atomization one encounters while using the standard method (STD) that does not involve consistency between the discrete transport of mass and momentum. The red contours indicate the isosurface of the volume fraction field corresponding to a value of 0.5, acting as a good proxy for the exact interfacial position. One can clearly observe that the un-physical fragmentation of the raindrop is symptomatic of the nonconsistent method, systematically across all combinations of flux limiters and advection schemes. The symbol "WY" in the legend corresponds to those run using the Weymouth-Yue advection scheme, and "CIAM" corresponds to the CIAM scheme. Brief descriptions of these advection schemes can be found in chapter . The implementation of the non-linear flux limiters i.e WENO, ENO, QUICK, Superbee, Verstappen correspond to that of well established methods developed to deal with hyperbolic conservation laws, for more details refer to the studies of Leveque [ ] and Sweby [ ]. The time stamps are indicative of the moments at which the interface is the most deformed, and do not necessarily correspond to the moment at which the code crashes. properties of the numerical system are examined within this time frame. A brief review of the relevant time scales are as follows :

The time scale t a D/U 0 of the air flow around the droplet, around . milliseconds, much shorter than the simulation time. The time t w L/[2(U t -U 0 )] 3ms that the droplet would take to travel by half the domain the domain once it had reached the terminal velocity. This time is not relevant here since one needs to wait first for the next two times before terminal velocity is reached.

The time scale t c 15.1ms [ ] of capillary oscillations of the droplet shape. The time scale t i of relaxation to terminal velocity. Using the dynamics ( . ) below, this time is (ρ l /ρ g ) • (D/U t ) 215ms, which is much longer than the simulation time. The time scale t µ needed to entrain the internal vortical motion of the liquid under the action of the gas. An estimate this time is D 2 /ν l ReD/U t 400ms.

The time of relaxation may be estimated using a simple square-velocity drag law for the droplet. We model the droplet motion as a onedimensional dynamics under the effect of gravity and drag as The simulations were carried out with a consistent scheme based on the shifted fractions (MSHIFT) method, using the WY advection scheme with the QUICK limiter. The droplet velocities correspond to that of their respective center of masses. Inset : Convergence of the droplet acceleration as a function of its resolution, computed using the best linear fit over the temporal variation of their respective velocities. The error bars signify the asymptotic standard error (least-squares) corresponding to the linear fits.

Figure . :

Comparison of droplet velocity as a function of time, for different droplet resolutions. The simulations were carried out with a consistent scheme based on the sub-grid (MSUB) method, using the WY advection scheme with the QUICK limiter. The droplet velocities correspond to that of their respective center of masses. Inset : Convergence of the droplet acceleration as a function of its resolution, computed using the best linear fit over the temporal variation of their respective velocities. The error bars signify the asymptotic standard error (least-squares) corresponding to the linear fits.

ρ l πD 3 6 dU dt ρ l πD 3 6 g -C D ρ g πD 2 8 U 2 , ( . ) hence dU dt - 3 4 r D (U 2 -U 2 t ) - U -U t t i , ( .
)
where for U U t , thus giving us

t i 2 3 D rU t . ( . )
We illustrate the performance of the consistent methods through the results of our simulations in figures . and . for the shifted fractions (MSHIFT) and sub-grid (MSUB) methods respectively. In case of both the methods, we use the WY advection scheme in combination with the QUICK flux limiter, while keeping the same value for the inflow velocity boundary condition. The quantities of interest while examining the robustness of the method are the temporal evolution of the droplet kinetic energy (figures . . where D is the computational domain and x m is the distance of the interface relative to the center of mass of the droplet. The droplet kinetic energy is defined relative to the droplet center-of-mass, given by

E k ρ(x, t)||u(x, t) -u CM || 2 ( . )
where . is the spatial averaging operator over the entire domain and u CM is the droplet center of mass.

The kinetic energy of the droplet evolves in a relatively smooth manner, without the presence of sudden spikes and falls which are emblematic of the standard non-consistent method (refer to Fig. . ). Such abrupt changes in kinetic energy of the droplet have been found to be associated with instants when the droplet undergoes 'artificial' atomization or breakup, henceforth resulting in catastrophic loss of stability for the numerical method. We observe a systematic drop in the amount of the droplet kinetic energy as we increase resolution, with the most probable explanation being that of the suppression of spurious interfacial jitter which is rampant at lower resolutions. There is also a component of the kinetic energy of the droplet associated with the internal coherent vortical structures generated due to the interaction of aerodynamic shear at the interface, evidenced by the non-zero value of the kinetic energy even for the most highly resolved droplets. Finally, the moments of inertia of the droplet appear to evolve in a smooth manner for both our methods, for all droplet resolutions.

In figures . and . , we show the velocity of the center of mass of the droplet (in the frame of reference of the box enclosing the droplet) as a function of time, and its behavior as we increase the droplet resolution. We have to keep in mind that the velocity inflow condition of 5 m/s does not correspond to the terminal velocity field of an actual falling raindrop might experience, hence the droplet in our numerical setup has some near constant acceleration for very small times once the simulation are started. The temporal variation in the droplet velocity is fitted to a linear polynomial in order to evaluate the droplet acceleration by means of a standard least-squares approach, for each droplet resolution. The observed accelerations for the MSHIFT and MSUB methods are dU dt 5.8 ± 0.1 m/s 2 and dU dt 4.5 ± 0.4 m/s 2 respectively, the error estimates being obtained from the difference between the D/∆x 2 5 and the D/∆x 2 6 resolutions. The disparity between the two methods concerning the acceleration of the droplet center of mass is approximately 20%, the exact cause of this discrepancy is yet to be ascertained. The decrease of the error, obtained as above from the difference between succesive resolutions using the MSHIFT method is plotted on Fig. . (d), which shows the lack of any apparent order, although the second-order convergence line is plotted as a reference.

In Fig. . we show the mass conservation properties of the shifted fractions method (MSHIFT) corresponding to the two most stable combinations, that is the WY advection with the QUICK flux limiter and the CIAM advection with the Superbee limiter. The mass conservation of the WY advection is strongly dependent on how accurately the divergencefree condition is enforced, which in turns is determined by the Poisson's solver tolerance. A clipping procedure, with no redistribution, affects mass conservation as well. The WY combination is thus rather sensitive to these two parameters but overall performs much better than the CIAM combination which is inherently not mass-conserving and not very sensitive to the Poisson's solver tolerance and the clipping parameter . Mass conservation is also plotted with various resolutions D/h in Fig. . c. It is seen that going from D/h 32 to D/h 64 the mass variation becomes larger for both methods, and grows over time at high resolution for WY, pointing to an accumulation of many machine precision errors.

To sum it up, the results obtained in this chapter strongly suggest that our class of consistent numerical methods can be used to get relatively good estimates of the underlying flow features of the droplet without observing any un-physical evolution due to the discretization errors at low to moderate resolutions. Ligaments constitute the penultimate stage in the complex sequence of capillary-driven topological changes that are typical of liquid fragmentation processes, finally resulting in the generation of polydisperse collections of drops. In the context of the present body of work, we shall only concern ourselves with the dynamics of Newtonian ligaments, which are well described by the Navier-Stokes equations at the incompressible and isothermal limits. The rearrangement of the liquid volumes that constitute the ligaments play a key role in determining the size of the droplets that emerge immediately after the disintegration of the thread-like structure (Fig. . ).

P F

The dynamics of these rearrangements prior to breakup are governed by non-linear interactions between several physical mechanisms such as the growth and propagation of capillary waves along the ligament surface, remnants of the internal flow, stretching induced either by the surrounding gas flow or by acceleration of the liquid into the surrounding medium itself, not to mention the dissipative effects of viscosity. Although we present a convoluted picture of ligament dynamics prior to breakup, it is ultimately the capillary force that drives the eventual topology change from the thread-like structure to that of the drop. Often the role played by both inertia and viscous forces is to oppose or dampen the destabilizing effects of the capillarity induced deformations. An additional mechanism to consider is the merger or coalescence of the newly created droplets (post-breakup) along the axis of the former ligament structure, or with droplets originating from other sources in the near vicinity. The generation of drops whose diameters are significantly larger than the width of original ligament are predominantly driven by the aforementioned coalescences.

. Shape Characterization

Linear theory based on the Rayleigh breakup of infinitely long liquid cylinders in a quiescent medium predicts the disintegration of the thread along regular intervals, the length of which is governed by the fastest growing spatial frequency [ , ]. As an extrapolation of this model, one can consider the size of the drops formed to be determined by the liquid volume contained within one wavelength corresponding to the fastest growing mode, thereby resulting in a uniform drop size distribution. This is in sharp contrast to the significant variances in the size distributions observed in reality . Even in the simple case of decaying liquid jets, non-linearities near the breakup zone kick into effect beyond the initial exponential growth phase predicted by linear theory, eventually resulting in the jet breaking up into "main" drops and significantly smaller "satellite" droplets [ ]. Thus, simply taking into account the non-linear effects driving the breakup would lead us to expect a bimodal distribution of droplet size, but this explanation still fails to account for the broad distribution of sizes observed in the outcomes of the majority of natural [ ] and industrial [ ] processes. asserts that the variance in the droplet sizes is strongly correlated to the initially corrugated shape of the ligament from which the drops originate. In this model, the corrugated ligament is represented as a collection of 'blobs' (see Fig. . ), with the continuous interaction between such blobs throughout the destabilization phase accounting for the rearrangement and concomitant aggregation of the liquid volumes prior to breakup.

There are several hypotheses in existing literature [ ]

The corrugation-coalescence mechanisn has been popularized in several studies such as [ , , -], and more recently in [ , ].

The key parameter concerning the predictions of this model is n, which characterizes the corrugations in the initial geometrical shape of the ligament, defined as

n ≡ 1 ( d 2 -d 2 ) / d 2 , ( .
)
where the quantities in the brackets represent the mean and the diameters d i correspond to that of the blobs as shown in Fig. . . These studies essentially demonstrate that the degree of polydispersity in the drop sizes can be simply explained by examining the degree of "smoothness" in the initial geometry of the ligaments, with the standard deviation or distribution width of the size of the emerging drops scaling as 1/

√

n. As a corollary to this prediction, a broader distribution of drop sizes arise out of the disintegration of strongly corrugated (rough) ligaments, whereas smoother ligaments with uniform thickness lead to a narrower (monodisperse) distributions.

Over the last two decades, a handful of studies have been conducted with the aim of carrying out detailed observations surrounding the aforementioned "coalescence cascade", whereby the drops constitutive of the ligament tend to coalesce and form larger drops as and when they detach from the initial thread-like structure. Most notable amongst those studies is the one by Marmottant and Villermaux [ ], wherein they illustrate the successive stages of the disintegration of liquid ligaments formed when capillary tubes are rapidly retracted from the surface of a liquid pool (refer to The fundamental limitation that is common to experiments cited above is the inability to precisely control the initial conditions. In other words, generating ligaments conforming 'exactly' to a specified geometrical shape is practically impossible due to the inherent difficulties associated with the control of free-surface flows. In the absence of such precise control and reproducibility of the initial ligament corrugations, one has to resort to a posteriori correlations between the width of the final drop size distributions and qualitative descriptions of the initial surface profiles.

: The characteristic time-scale of the flow near the vicinity of the breakup is given by t σ ∼ ρr 3 /σ 1/2 , thus leading to faster flows and contraction of the time-scales as the local thickness (r) decreases to zero. In our numerical methods, we can adaptively select the time interval over which the integration of the governing equations are carried out, to accommodate such nonlinearly accelerating flows. A central theme of the current chapter is concerning the design and conception of 'numerical' experiments that lend themselves to accurate and repeatable specifications of the initial conditions of the ligaments in question. An additional advantage offered by our high-fidelity numerical approach is the access to a much wider range (compared to physical experiments) within the space of all control parameters that influence the dispersion in the drop size distributions, also including higher levels of temporal resolution on offer that enable better tracking of the non-linear dynamics near pinch-off . The combination of these two features enables us to quantitatively map the influence that the degree of smoothness of the initial ligament has on the ensuing liquid rearrangements, as well as the coalescence dynamics following the ligament.

. Numerical Setup

Our numerical experiments entail direct numerical simulations of the two-phase Navier-Stokes equations, employing the one-fluid model which solves only a single set of governing equations throughout the whole domain, where the domains are composed of regions of variable density and viscosity. In order to reproduce the free-surface dynamics, a Volume-of-Fluid (VOF) approach is used to capture the interface, the explicit location of which marks the separation between the two fluids. The details regarding our VOF approach, governing equations etc are covered in Chapter .

A more detailed evaluation of the comparative advantages of different interfacial tracking approaches (e.g. VOF, level set, phase field etc) can be found in the comprehensive review by The bottom side of the box acts as the axis of symmetry for the ligament, with spatial periodicity along the horizontal direction. Inset : A close up view of the corrugated profile of the ligament, where the local radius is defined as the sum of the unperturbed (mean) radius R and the local perturbation (x). The material properties of the liquid and gas phases are denoted with the subscripts l and g respectively, which in our case corresponds to an air-water system with the surface tension coefficient σ.

: We use inertial forces as an umbrella term to describe all the forces that might contribute to the "stretching" of the ligament, thus opposing the capillary deformations. These primarily include the remnant flow within the ligament bulk, the aerodynamic shear at the interface and the acceleration of the ligament itself into the surrounding gas medium. created due to bursting bubbles [ , ] modeling of viscoelastic droplet splashes [ ] and many more.

Computational Setup

As discussed in the previous sections, the dynamics of the liquid volume rearrangements within the ligaments are quite involved due to the nontrivial interactions between the inertial, viscous and capillary effects at different length and time scales that span over several orders of magnitude. In order to isolate the direct influence of initial geometrical shape of the ligament on the subsequent dynamics, we have decided to exclude the inertial forces with respect to the initial conditions. Therefore, we opt for a particular setup, with infinitely long ligaments having zero initial kinetic energy and spatial periodicity in the geometrical shape along the axis. We carry out the simulations in an axisymmetric framework, therefore excluding any azimuthal variations in the initial shape as a consequence of this simplification. Fig. . illustrates the a schematic of the computational setup, where the domain is a square of side L. The bottom side of the box acts as the axis of symmetry for the corrugated ligament (detailed view in the inset of Fig. . ), which has an unperturbed (mean) radius R. Due to the periodicity of the ligaments in the axial direction, we impose periodic boundary conditions for the primary variables on the left and right faces of the domain. Symmetry conditions are imposed on the bottom side, with impenetrable free-slip conditions applied on the top side. : Radial profile of the noisy and the subsequently low-pass filtered surface of the ligament, corresponding to one spatial period along its axis. Inset : Spectral densities of the initial noisy signal and the corresponding filtered signal. Only a handful of small wavenumbers are chosen to create the surface.

: As one can infer from the inset of Fig.

. , the creation of our corrugated profiles are based on the popular model in information theory termed as white Gaussian noise, where white refers to the flat spectrum in frequency space, and Gaussian due to the normal distribution in physical space.

Random Surface Generation

A key ingredient of this investigation on the impact of corrugations is the ability to precisely control the quantitative aspects of the initial corrugated surface profile of the ligament. Towards that objective, we construct a surface using the superposition of random overlapping waves, corresponding to an overall spatial period L of the resulting profile. The aspect of "randomness" is achieved by using a state of the art white noise generation algorithm [ ] . A robust and open source implementation of the aforementioned algorithm is available in the GNU Scientific Library toolkit [ ], and thus readily incorporated into our numerical platform Basilisk.

The discretization of the radial profile R i can be expressed as the sum of the mean radius (unperturbed) and some random perturbation such that

R i R + λ i , λ i ∼ N(0, 2 0 ) , ( . )
where i represents the discretization of the continuous profile along the axis of symmetry, and λ i corresponds to random numbers generated by our white noise algorithm , with the numbers normally distributed with zero mean and 2 0 as variance . The exact sequence of numbers λ i is uniquely determined by the value of the seed given as input to the white noise generator. This ensures that there are no unforeseen correlations between the shapes of any two ligament surfaces, even though some long range correlations are introduced within the individual ligaments themselves due to the filtering process. The initial noisy surface generated by the white noise algorithm contains extremely sharp edges corresponding to infinite values of interfacial curvature (see Fig. . ), which are present due to the inclusion of perturbations corresponding to large wavenumbers.

It is interesting to note that we can independently control the dispersion (standard deviation) in the final ligament profile via two means, first is the variance of the initial noise ( 0 ) and second is the threshold used for the low-pass filter.

In order to avoid the "smoothening" of these edges by our numerical method, we remove these high frequency modes a priori in a controllable fashion. As demonstrated in Fig. . , the noisy surface is treated by a low-pass filter in order to obtain the final description of the ligament's radial profile.

: The viscous length is defined as l µ µ 2 /ρσ. It is intrinsic to the system of the liquid and its free-surface, therefore is not influenced by initial conditions or any global parameters specific to the problem. : In our particular case, only the modes corresponding to the the longest n c wavelengths (L, L/2, L/3 ... L/n c ) are left with any excitation once the initial noisy profile is passed through the low-pass filter.

Adimensional Groups

The initial conditions of the computational setup can be characterized by the Ohnesorge number, as defined below Oh µ/ ρσR.

( . )

The Ohnesorge number is simply the square-root of the ratio of the viscous length scale (l µ ) to the characteristic length scale of the problem (R). The ligaments have zero initial kinetic energy, but a part of the surface potential from the static initial conditions is immediately converted into liquid inertia, hence we can interpret that our setup is characterized by We 1, where We is the Weber number. The geometrical shape of the ligament can be characterized by the combination of a measure of the initial corrugation amplitude (ε), cut-off wavenumber (K) and the aspect ratio (Λ), defined as

ε 0 /R , K n c • (2πW/L) , Λ L/W, ( . )
where W 2R denotes the mean width of the ligament, with only the n c smallest (discrete) wavenumbers passing through the low-pass filter .

The volume of the corrugated ligament per unit spatial period is controlled by the aspect ratio, thus indirectly setting the upper limit to the size of the largest drop that could thoeretically form if all the volume of the ligament were to congregate into a single shape.

It is also important to note that the normalized corrugation amplitude is defined based on the dispersion of the initial white noise signal ( 0 ), and not on the dispersion of the "filtered" shape of the final profile. As in the case with most surface tension dominated flows, the capillary time scale (t σ (ρR 3 /σ) 1/2 ) is the most natural choice of time interval by which we can rescale our physical time. Thus, our non-dimensional time is given by

T t • ρR 3 /σ -1/2 . ( . )

Active Wavenumbers

As a consequence of the discrete nature of space in the context of numerical simulations, we can inject energy only into certain spatial frequencies . The general and important result that a liquid cylinder is unstable only to perturbations of wavelength larger than π times the diameter is was first demonstrated in the original work of Plateau [ ].

The subsequent evolution of these low frequency modes is determined by the dynamical response of the system in question, subject to the underlying governing equations. In the context of inviscid flows, the dispersion relation resulting from linear stability analysis of the potential flow inside a liquid cylinder subject to sinusoidal perturbations is given by ], expressing the growth rate of the small amplitude perturbations as a function of the reduced wavenumber ( . ). The vertical dashed line A represents the upper limit to the perturbation wavelength, corresponding to the domain size L. At the opposite end of the spectrum, E represents smallest wavelength limited by the grid size 1/2 M , where M is the maximum level of resolution for our octree grid. B corresponds to the wavenumber with the maximum growth based on linear theory, whereas C represents the limit beyond which all perturbations are stable. The number of discrete modes that fall within the unstable part of the spectrum (∆K) is a linear function of the ligament aspect ratio.

ω ω 0 I 1 (x) I 0 (x) 1 -x 2 , ( . )
: Beyond the exponential and independent growth phase of the modes in the linear regime, non-linear effects eventually kick into place, hence the interaction between the different modes. I 0 and I 1 are modified Bessel's functions of the first kind, corresponding to order 0 and 1 respectively. Detailed derivations of ( . ) and ( . ) can be found in the comprehensive review by Eggers [ ]. where ω 0 σ/ρR 3 and x kR is the reduced wavenumber of the initial perturbations. The relation ( . ) describing the classical Rayleigh-Plateau instability [ , ] demonstrates that only the perturbations corresponding to x ≤ 1 undergo an exponential growth phase, while the remaining perturbations are immediately damped out.

A similar analysis for the full Navier-Stokes equations was carried out by Chandrasekhar [ ]. The corresponding dispersion relation for viscous flows incorporating the long wavelength approximation of Weber [ ], at the limit of low Reynolds is given by

ω ω µ 6 1 -x 2 , ( .
)
where ω µ σ/µR. A similar conclusion can be drawn from ( . ), that all the perturbations corresponding to x ≥ 1 are stable with respect to linear theory. In Fig.

. , we illustrate this principle using the dispersion relation of the inviscid Rayleigh-Plateau instability as a toy model. As one can observe, the number of discrete modes in the unstable part of the spectrum amongst the n c modes initially excited is given by ∆k Λ/π -1.

( . )

Thus, the aspect ratio (Λ) of our ligaments is the parameter responsible for selecting the number of modes actually playing a role in the destabilization of the ligament into the eventual drops. To this effect, in the present study we opt to study the dynamics of slender ligaments i.e. large aspect ratios, so that we excite a large enough number of discrete frequencies that are unstable with respect to linear theory. The expectation is that the competition and interactions between these unstable modes will play a key role in the complex liquid rearrangements within the ligament bulk.

Numerical Fine Tuning

As we have discussed in the previous part of the thesis, one has to pay particular attention to several aspects pertaining to the numerical stability of solvers when dealing with surface tension dominated flows involving significant contrasts in material properties across the interface (such as our system of water ligaments in air). In our case with 'Basilisk', we opt to use the standard non-conservative formulation of the Navier-Stokes equations, which does not maintain consistency between mass and momentum transport. But, in this case we are able to circumvent effects such as 'artificial atomization' mainly due to two reasons :

The use of adaptive mesh refinement based on the octree grids in Basilisk allows us to attain much higher levels of numerical resolution for the smaller/thinner liquid structures, while using significantly coarser meshes away from the interface. For example, in our current simulations (see Fig. . ), the number of cells in the domain are of the order 10 4 , whereas an equivalent uniform mesh would result in a problem size of order 10 6 , thereby differing by a factor of . The absence of the interaction of aerodynamic shear forces with the interface, in our particular problem configuration. This prohibits the 'numerical mixing' of the gas velocity with the liquid density in the vicinity of poorly resolved interfacial cells, the concomitant unphysical momentum transfer resulting in the rapid destabilization of the interface. Furthermore, we have to tune the 'knobs' of the solver concerning the criteria and corresponding thresholds regarding the dynamic mesh refinement. In all of the simulations/results demonstrated from this point onwards, the mesh adaptation is made according to the value of the local gradients of the volume fraction field, as well as local gradients of the velocity field components. The corresponding values of the thresholds are 10 -3 and 10 -1 for the volume fraction and velocity components respectively. The time step for the numerical integration is in itself dependent on the size of the smallest cells, as well as the local velocities. Therefore, we opt for a fixed CFL criteria in order to choose the time step, with the upper bound of the local CFL number throughout the domain being set to 0.05.

Another numerical parameter that has a big influence on the performance of the numerical method is what we refer to as "droplet removal", which essentially removes extremely small and poorly resolved structures that float around in the domain. In our experience, these small 'debris' like structures (see Fig. . ) hamper the simulations in several ways e.g. demanding extra computational resources in the form of additionally refined cells, or by imposing extremely small time steps (due to restrictions of local CFL number based on cell size) used for the numerical integration. The removal of these structures do not have any effect on the underlying dynamics of the breakups and coalescences, and their sizes are verified to be orders of magnitude below those corresponding to satellite drops. Our threshold for removal corresponds to structures that are resolved by cells or less. With regards to the iterative Poisson solver, we set a tolerance of 10 -3 , which in our operator-split framework (see ( . )) corresponds to the accuracy by which the incompressibility condition ( . ) is verified.

Axisymmetric D versus Full D

An important question to ask ourselves before we conduct any large scale study is the feasibility of carrying out full D simulations. In Basilisk, we have the option of carrying out axisymmetric simulations that resolve the underlying governing equations on cylindrical coordinates. Due to the axial symmetry in the construction of our ligaments, carrying out axisymmetric D simulations instead of full D ones seem to offer a more computationally efficient approach, especially as we are interested in studying the dynamics of thousands of such ligaments. We start with ligaments that correspond to the following parameters Oh 10 -2 , ε 0.8 , ( . ) Λ 50 , K 2π .

( . )

The numerical resolution in our feasibility studies (in both axisymmetric D and full D) correspond to W/∆x min 10.24, where ∆x min is the size of the smallest grid cell in the octree mesh. The computations were carried out on processors in parallel, within the framework of parallelization through domain decomposition. The computational times required to resolve the dynamics up to T 30 are 54 seconds and 976 seconds for the axisymmetric D and full D cases, respectively. In figures . and . , we illustrate the differences in the temporal evolution between a full D simulation and its axisymmetric D (axisymmetric) equivalent, on fast and slow time scales respectively. As one can observe, the major qualitative difference lies in the excitation of certain azimuthal modes in the full D case, which cannot be captured in the axisymmetric D case due to the imposition of symmetry about the central axis. Such a subtle difference in the dynamics may explain the slightly higher tendency towards coalescence in the axisymmetric D case, due to the fact that some of the energy from the axial modes is transferred into the azimuthal modes in the full D case. Apart from that, the full D simulations also resolve some finer thread like structures that appear in the vicinity of the breakup region, the details of which are missing in the axisymmetric D simulations. Overlooking some finer aspects captured by the full D simulations, there seem to be no major discernible differences between the axisymmetric D and full D simulations, with the axisymmetric D simulations acting as an accurate reduced order model concerning the breakup and coalescence dynamics over both short and long time scales. Looking at the feasibility in terms of computational resources, the full D simulations are substantially more expensive than their axisymmetric D counterparts, by almost a factor of 20 in most cases. Considering the fact that we plan to carry out thousands of such simulations, the slightly more accurate dynamics modeled by the full D simulations do not justify the approximately fold increase in computational cost. Hence, from this point onwards, unless mentioned otherwise, all of the simulations or results we present correspond to the axisymmetric D axisymmetric framework.

: These simulations are carried out in parallel on processors.

Spatial Resolution

Arguably the most important parameter on which the fidelity of our numerical simulations hinges on is the spatial resolution or grid size. The octree structure in our solver entails successive refinements in the size of the cubic cell (square cell in axisymmetric D) by factors of 2 in each direction.

In our simulations, the square domain of size L has a minimum level of refinement corresponding to a grid size of L/∆x 2 4 , where ∆x is the size of the square cell.

As we have discussed before, we use the adaptive mesh refinement capabilities of Basilisk to add greater resolution to regions in the vicinity of the interface, therefore we have to impose an upper limit to the level of refinement as well. Henceforth, we refer to this upper limit as "maxlevel", where a "maxlevel" of 10 corresponds to a grid size of L/∆x 2 -10 , In Fig. . we illustrate the differences between simulations carried out using three different "maxlevels", corresponding to identical ligaments that are characterized by the non-dimensional parameters given in ( . ). The computational cost of the most resolved simulation characterized by a "maxlevel" of 11 (W/∆x min 40.96) is 516 seconds , for "maxlevel" of (W/∆x min 20.48) is 156 seconds, and finally for "maxlevel" of (W/∆x min 10.24) is 54 seconds. Comparing the simulations at the different resolutions, we do not observe any substantial differences in the dynamics at intermediate as well as long time scales.

As a minor remark, we point out that coalescences take place marginally faster for the lower resolutions compared to the one as the "maxlevel" of 11.

Judging by the qualitative and quantitative similarities between the simulations corresponding to "maxlevels" 10 and 11, we conclude that our simulations are converged (more or less), and that a "maxlevel" of 10 serves as an optimal level of refinement for our future simulations.

. Impact of Initial Conditions

The characteristic length scale of the problem is the radius R of our ligaments, which are in turn determined by the choice of Ohnesorge number, given that we are dealing with air-water systems. On the other hand, the global dynamics of the ligament breakup take place over a multitude of different time scales, in particular the penultimate stages of the non-linear breakup which occur on successively shrinking (faster) time scales. We direct our attention towards understanding the dynamics over much larger time scales rather than those associated with the self-similar collapse of the liquid bridges that lead to the formation of drops. In this study, we loosely define short time scales as intervals ∆t corresponding to ∆t ≤ 2T, intermediate ones as 2 < ∆t ≤ 5T and long time scales as intervals ∆t > 5T. In what follows, we briefly describe the influence of the adimensional groups which precisely define the initial conditions for our ligament setup.

Corrugation Amplitude

The parameter ε acts as a quantitative measure of the "smoothness" of the initial surface profile, being defined as the standard deviation of the un-filtered white noise signal, rescaled by the mean radius of the ligament. Figures . and . illustrates the differences brought about in the temporal evolution of ligaments (over both short and long time scales) with identical surface profiles due to the differences in the average strength of the perturbations corresponding to that particular set of frequencies. As one can clearly observe, the stronger the corrugations, the lesser is the time available for rearrangement of the liquid volumes within the ligament, which completely changes the trajectory of the system dynamics involving the subsequent coalescences.

Ohnesorge Number

The Ohnesorge number selects the characteristic size R of our ligaments, where Oh O(10 -2 ) and Oh O(10 -3 ) corresponds to millimeter and centimeter scale water ligaments in air, respectively. Fig. . outlines the temporal evolution of ligaments with the exact same initial shapes, but corresponding to different length scales.

The interfacial shapes become more compact as one goes towards larger Ohnesorge numbers, also leading to slower flows and lower amounts of vorticity concentration at the interface. These trends are quite expected as the effects of visous dissipation increase in influence as one goes towards smaller length scales (R ∼ Oh -2 ).

Aspect Ratio

The parameter Λ gives us a measure of the "slenderness" of the ligaments, and more importantly governs the total volume of liquid contained within the structure. Fig.

. compares the temporal evolution of ligaments described by identical shapes and mean radii, but differing in their lengths by a factor of 2. As expected, the longer ligament has a greater number of drops formed (possibly with a larger dispersion in sizes as well), compared to the shorter one over similar time scales. In our particular setup, large aspect ratios are necessary to obtain drops at the large size limit of the probability distributions, as the size of the largest possible drop is set by the ligament volume.

Cut-Off Frequency

The cut-off parameter K controls the threshold for the low-pass filter which is used to treat the initial white noise signal. Fig. . shows the time evolution of ligaments differing in their initial geometries due to the addition of energy into successively higher frequency modes. Beyong a certain point, adding more high frequency perturbations into the initial geometric shape do not seem to have any significant impact on the ensuing dynamics, neither prior or post breakup. This can be partly explained by the fact that the number of waves that actively play a role in the destabilization of the ligament are determined by the aspect ratio, as described in the previous section.

Thus far, we have described our methodology concerning the creation of ligaments with randomly corrugated surfaces, where the corrugations are subject to precise quantitative descriptions. We have also briefly illustrated the differences in the dynamics of the ligament destabilization as well as the ensuing coalescences, brought about by the differences in initial conditions of the ligaments in question. Nonetheless, it remains a difficult task to quantify the influence of the above parameters on the resulting drop sizes, especially as we have a very small sample of drops obtained via breakup of the individual ligaments. Thus, in the following chapter we utilise the current framework in order to generate large ensembles of such randomly corrugated ligaments, the subsequent destabilization of which will provide us with drops corresponding to adequately large sample sizes, hence enabling us to carry out proper statistical descriptions pertaining to the influence of the initial conditions. The destabilization and eventual breakup of ligaments are triggered by the growth of undulations that define the initial geometrical shape. In the preceding chapter, we have presented our numerical setup pertaining to ligaments that are initially at rest, the axisymmetric shapes of which are generated in a precisely controllable (and reproducible) manner, entailing the superposition of randomly overlapping waves. The resulting individual ligaments are characterised primarily by the characteristic length scale, aspect ratio and the variance of the initial white noise signal used to create the randomly corrugated shape. Direct numerical simulations (DNS) performed on these ligament shapes provide us with important insights, primarily regarding the differences in the trajectories of the destabilization and subsequent coalescence dynamics as a function of the set of initial conditions used. In principle we can generate an infinite number of unique initial corrugated shapes using our white noise algorithm, each shape characterized by an identical set of adimensional parameters (Oh, ε, Λ and K) . The properties of the drops generated via the ligament breakup using the DNS is in strict correspondence to the exact initial shape of the ligament due to the deterministic nature of the underlying governing equations. Therefore, even a slight change in initial conditions lead us to a different collection of drops, which may or may not be representative of the total population of drops associated to that particular combination of dimensionless parameters.

In an effort to address this issue, stochasticity is introduced into the mix by conducting an ensemble of direct numerical simulations of our slender corrugated ligaments, where each individual realization corresponds to a random but unique initial configuration. This approach can be viewed as the application of Monte Carlo methods to DNS, where uncertainties in the initial conditions are leveraged to provide a statistical description of problems which are deterministic in principle.

. Millimeter Scale Ensembles

We now turn our attention towards slender ligaments whose widths are of the order of millimeters, a scale at which a multitude of fragmentation phenomena take place that are more or less visible to the naked eye.

In particular, for air-water systems ( degree Celsius) the millimeter length scale corresponds to Oh ∼ O(10 -2 ), where Oh is the Ohnesorge number based on the characteristic length scale of the ligaments. The characterization of our corrugated ligament setup using the set of adimensional (equation ( . )) numbers as discussed in the previous chapter is extended to describe our ensemble of ligaments, such that the exact set of parameters correspond to a unique point in the phase space Φ of all possible combinations. Coalescence that intermittently occurs after the initial breakup phase of the thread-like structure not only change the number of drops, but also impacts the mean and dispersion of the size distribution due to the increase in the number of drops significantly larger than the width of the ligament of origin. This aspect of temporal dependence of the drop sizes corresponding to a particular ligament ensemble is taken into account by specifying the time at which the statistics are recorded. Thus, our droplet ensembles associated with the millimeter scale ligaments can be uniquely described by a combination of the time T and Φ 0 , where

Φ 0 ≡ Oh 10 -2 , K 2π, 1.0, Λ 50 . ( . )
The definition of the adimensional groups are identical to that used in the previous chapter. Statistical descriptions of the drop sizes are commonly carried out using particular probability density functions (PDF) defined as follows Gaussian : P x; µ, σ

1 σ √ 2π exp - 1 2 x -µ σ 2 , ( . )
Log-Normal : P x; µ, σ

1 xσ √ 2π exp - 1 2 log x -µ σ 2 , ( . ) Gamma : P (x; n, k) k n Γ(n) x n-1 exp (-kx) , ( .
)
where x is the variable (diameter, volume etc) in question. The Gamma distribution can be rearranged so that x represents the variable rescaled by its mean (n/k), which becomes

P (x; n) n n Γ(n) x n-1 exp (-nx) . ( . )
The change of variable renders the function dependent solely on the parameter n, identical to the form first used in [ ] to represent the distribution of the normalized diameter. At the large size limit of the distributions, the leading order terms scale as P(x) ∼ exp(-x 2 ) for the Gaussian, P(x) ∼ exp(-log x) for the Log-Normal and finally P(x) ∼ exp(-x) in case of the Gamma function. In essence, out of the three density functions, the Guassian predicts the lowest frequency of rare events (very large drops), the Log-Normal predicts the highest frequency whereas the Gamma function prediction is somewhere in between the other two.

Another important aspect to consider is the bin-width of the histograms used to construct the probability distributions. There are several choices regarding the criteria used to determine the optimal bin-width, taking into consideration the size, variability and skewness of the underlying data set. In the interest of simplicity, we restrict our focus to histograms with uniform bin-width. In order to ascertain the dependence (if any) of the distribution on the number of bins (uniform size), we use 4 different estimators which can be found in the Numpy library [ ]. The simplest : The error is defined as the L 2 norm of the differences between the bin hights of the sample and the largest sample, where both samples have identical bin-width.

of the four estimators selects the number of bins as the square root of the sample size. Another estimator is also based solely on the size of the dataset, but the number of bins scales as the cube root of the data size. A more robust estimator that takes into account the variability in the dataset is the Freedman Diaconis estimator, in which the bin width is proportional to the inter quartile range, and inversely proportional to the cube root of the data size. Finally, we also use an improved version of the Sturges estimator, that performs well for non-normal datasets due to the fact that it takes into account the skewness.

In the subsequent sections, we present a statistical picture of the drop sizes generated from our ligament ensemble Φ 0 , where the sample size is of the order of 64, 000. The droplet data is primarily recorded on slow time scales, so that the data reflects a considerable amount of coalescences after the initial destabilization of the ligaments. In order to extract the diameter and volumes of the drops generated in our axisymmetric simulations, we for the commonly used area based estimation (see Fig. . ).

Diameter Distributions

Figures . and . illustrate the probability distributions of the normalized droplet diameter corresponding to T 15 and T 30 respectively. The diameters are normalized using the means of the associated samples. The bin widths are always based on the properties of the largest sample. As one can observe, there are fewer drops in all samples at T 30 compared to the corresponding ones at T 15, due to the additional coalescences. The drop diameters are found to be more or less symmetrically distributed about the mean, although there is an small peak at the lower end of the distribution. The choice of binning criteria does not seem to have any noticeable effect on the overall shape of the distribution, regardless of the size of the sample N. A Gaussian probability density function ( . ) appears to be a good approximation to the underlying distribution. It is important to note that the Gaussian curve is not fitted to the heights of the histogram bins, instead it is plotted on top of the distributions using the mean and variance corresponding to the largest sample size, therefore rendering it independent of the choice of bin-width and free from any additional fitting parameters.

In Fig. . , we present a comparison of the different probability density functions applied to the largest ensemble of drops characterized by Φ 0 and T 30. The Gaussian is found to be the best approximation to the normalized diameter distribution, followed by the Gamma and Log-Normal functions. The Gaussian and Log-Normal curves are generated using the mean and variance of the dataset, whereas the free parameter n of the Gamma function is determined according to the best least-squares fit over the histogram. A noteworthy point demonstrated by the inset of Fig. . is that the error scales as N -1/2 , N being the sample size. Additionally, the scaling is not sensitive to the choice of bin-width. This strongly suggests the absence of any unforeseen correlations between the individual realizations (corrugated ligaments) that form the ensemble.

Volume Distributions

We now shift our attention towards the distribution of droplet volumes. Figures . and. show the PDFs of the normalized droplet volume corresponding to T 15 and T 30 respectively. The volumes are normalized by the means of the associated samples. As in the diameter distributions, the bin widths are always based on the properties of the largest sample. In stark contrast to the case of the diameters, the drop volumes are found to be heavily skewed towards the smaller sizes. Similar to the diameter PDFs, the choice of bin-width seems to not have any considerable impact on the overall shape of the distribution, and is also independent of the sample size. In this case though, the single parameter Gamma density function ( . ) appears to be the approximation to the volume distribution, with the best fit corresponding to values of n ≈ 1.68 and n ≈ 1.58 corresponding to T 15 and T 30 respectively. There is a slight dependence of the parameter n on the bin-width, which is expected as n itself is determined according to the best fit to the histogram in question. ) nature of the diameter (D) distributions, a simple variable transformation leads us to volume (D 3 ) distributions of the form e -x 2/3

x 2/3 , which would result in a similar shape as that of the Gamma distribution. Similar to the diameter distribution fit, the Gaussian and Log-Normal curves are generated using the mean and variance of the dataset, whereas the free parameter n of the Gamma function is determined according to the best least-squares fit over the histogram. As seen before, the error (inset of Fig. . ) scales as N -1/2 , regardless of the binning criteria used.

Influence of Corrugation Amplitude

As discussed in the previous chapter, the "smoothness" of the initial geometrical shape of the ligament is expected to have a significant influence on the destabilization of the liquid thread and the concomitant coalescence dynamics. In order to ascertain the impact of the strength of the initial perturbations on the resulting drop size distribution, we conduct simulations of two different millimeter scale ensembles Φ 1 and Φ 2 , defined below as

Φ 1 ≡ Oh 10 -2 , K 2π, ε 0.4, Λ 50 , ( . ) Φ 2 ≡ Oh 10 -2 , K 2π, ε 1.6, Λ 50 . ( . )
For the case of weakly perturbed ligaments, Fig. . demonstrates the temporal evolution of the droplet diameter distributions, where the diameter is rescaled by the (mean) width W of the initial ligaments. At small time scales, the distribution is bimodal in nature, including a large number of drops that are smaller than the initial ligament width. This may be a consequence of the non-linear effects [ , ] that kick in immediately following the initial linear (exponential) growth phase, often resulting in the liquid thread disintegrating into a collection of main and satellite drops (see Fig. . in the preceding chapter). At large times, the peak corresponding to the smaller sizes is considerably diminished, whereas the larger size peak develops into more of a "plateau" type shape. This shift in the distribution in time can be explained by successive coalescence, where the smallest drops coalesce with the largest ones. By largest drops, we mean the ones situated to the right of the large size peak, corresponding to the region of the distribution which eventually gives rise to the plateau in the droplet distribution.

Coming to the strongly perturbed ligament shapes, Fig. . illustrates the unimodal character of the overall distribution shape. There seems to be a lack of any significant qualitative change when comparing the distributions at small and large time scales. The large amplitude perturbations of the ligament surface do not allow for any significant amount of liquid rearrangement within the bulk prior to disintegration, thus the collection of drops formed immediately following the breakup of the thread-like structure are quite disperse when with respect to their sizes. The mean of the distribution displays a slight shift towards the larger sizes, with the passage of time. This shift can be explained by the effect of random coalescences amongst droplets corresponding to all sizes.

To summarize, in Fig. . we present a side-by-side comparison of the aforementioned temporal evolutions of the droplet diameter distributions. We infer that the primary influence of the corrugation strength is that it changes the nature of the resulting distribution from bimodal to unimodal, as one increases the amplitude of the initial perturbations. Although the dynamics of the droplet coalescences (post-breakup) follow different trajectories in accordance with the corrugation strength (see Fig.

. ), the distributions at large times end up looking rather similar.

T : The expected value of n i (where n i ∼ B(N, p i )) is given by N p i , which is equal to the height of bin i. The corresponding variance of n i is given by N p i (1p i ). .

Description of Large Sizes

The approximately symmetrical diameter distributions that we have seen thus far seem to be relatively well described by a Gaussian probability density function. We shift our focus to the shape of the distribution towards the right hand side of the peak, which represents the large drops emerging from the coalescences. Although it is quite tempting to conclude that the PDF of the large sizes scale as e -x 2 (where x D/ D ), a quantitatively rigourous procedure is required to ascertain the exact shape in the vicinity of the PDF tail. In the distributions we have shown thus far, we have not taken into account the uncertainties corresponding to the bin heights, which might play a crucial role in determining which among the candidate scalings (such as e -x , e -x 2 etc) best describe the rare events. In order to compute the error bars associated to the bin heights, we define a Bernoulli trial in which a particular drop from the sample size of N either falls into bin i or not, with an associated probability p i of falling into that particular bin. Therefore, n i drops falling into bin i from a sample of N can be described by a binomial distribution. Now, we use a simple model for the bin height uncertainty by setting the length of the error bars as the standard deviation of the corresponding binomially distributed variable. Our candidate functions for the large size descriptions are of the form

y A • exp [-nx m ] , where x D/ D . ( . )
The goal is to search for the optimal parameter m such that the function ( . ) is the best fit to the bin heights corresponding to the tail region of the distribution. A non-linear least squares fitting method based on the Levenberg-Marquardt algorithm [ ] [ ]: Moré ( ), 'The Levenberg-Marquardt algorithm: implementation and theory' is utilized, through its open-source implementation in the SciPy library [ ]. This fitting routine can find the optimal parameters while taking into account the aforementioned uncertainty in the bin heights computed using the binomial variance. An important point to note is that even though we can simply find the best fit to a given set of points, the quality of the fit is also a function of the number of points used to minimize the least-squares error.

Therefore, we use a moving window of points (see Fig. . ) and quantify the quality of the fit as a function of the number of points used to estimate the optimal parameters. The "quality" of fit is quantified by the standard deviation of the estimate for "m", which is returned by the non-linear leastsquares fitting algorithm. Finally, the value of the optimal parameters that minimize the error (standard deviation of parameter estimate) over the different window sizes is deemed to be the best fit to our set of bin heights near the tail region. Additionally, we reframe the fitting problem on logarithmic scales, so as to make sure the optimal parameters found are not dependent on the rescaling of the axes. Rearranging equation ( . ) by taking successive logarithms, we obtain ) where Y log(-log y) and X log x, hence transforming it into a linear regression problem. The overall fitting procedure can be summarized in the following steps :

Y m • X + C, ( .
. The candidate function is chosen either in the exponential ( . ) or the logarithmic form ( . ). . A starting window size (greater than points) is chosen. . The least-squares routine is run using the points falling within the window, either by taking into account the uncertainty in the bin heights or not. . An optimal set of parameters (A, n and m) is found, along with the variance of the parameter estimate. . The above operations are repeated by successively adding more points into the window, till the edge of the window reaches the peak of the distribution. . The parameter m is selected corresponding to the number of points in the window for which the variance in the parameter estimate is minimum.

In figures . and . we illustrate the determination of the optimal parameter set using the moving window of points as descibed previously.

The optimal parameter set corresponding to the lowest standard deviation of the parameter estimate gives us m ≈ 6, with the robustness of the result established by the fact that it is insensitive to not only the form of the candidate function used, but also does not depend on whether the uncertainty in the bin heights are taken into account or not. This leads us to the rather unusual result that even though the global shape of the distribution varies as e -x 2 , at the limit of large x (see figures . , . , . and . ) the distribution tail is most accurately described by a function of the form e -x 6

. In what follows, we attempt to explain the origin of this scaling using a random process theory for near-monochromatic waves.

Figure . :

The tail of the distribution fitted with the function in the logarithmic form ( . ), using the optimal parameter set that corresponds to the lowest standard deviation in the parameter estimate. The uncertainty in the bin heights of the underlying distribution is not taken into account while conducting the least-squares procedure. The bottom figure is a closeup of the region in the vicinity of the tail. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes.

Figure . :

The tail of the distribution fitted with the function in the logarithmic form ( . ), using the optimal parameter set that corresponds to the lowest standard deviation in the parameter estimate. The uncertainty in the bin heights of the underlying distribution is taken into account while conducting the least-squares procedure. The bottom figure is a closeup of the region in the vicinity of the tail. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes.

. Theoretical Development

In the preceding section, we have obtained an exponential scaling for the probability distributions of drop diameter of the from exp[-d 6 ], at the limit of large sizes. Such a scaling predicts a significantly lower frequency of rare events than all of our candidate probability density functions such as the Gaussian, Log-Normal and Gamma. In an effort to explain this result, we develop a theoretical model based on the linearization of the Navier-Stokes equations for weakly corrugated ligaments, where the ligament destabilizes via the growth of perturbations corresponding to a random set of waves, with the assumption that the wavelengths are more or less normally distributed around the optimal wavelength of the Rayleigh-Plateau-Savart instability [ , ]. Note that the mathematical symbols, notations and variable names used in this section are unique to this section only.

In the surface tension dominated regime, it is convenient to rescale all lengths by the initial unperturbed radius a of the ligaments and times by the capillary time τ (ρa 3 /γ) 1/2 where γ is the surface tension. Our numerical setup is associated to a finite domain such that 0 < x < L. The system is initialized at zero velocity and with a local radius r(x, 0) 1 + h(x, 0)

h(x, 0) N n -N
A n e ik n x , ( . )

where k n 2πn/L and A n A -n . The amplitudes A n are Gaussian independent variables generated by a pseudo random number generator with identical variance v 0 A n A n N, ( . ) for all n. The cutoff frequency defined as k N k c 2π/L c allows us to obtain a smooth initial surface profile. The resulting perturbation h(x, 0) is also Gaussian with variance v 0 , which means that there is a finite probability

p i exp[-1/(2v 2 0 )] ( . )
that the radius corresponding to the initial condition reaches negative values, which is of course unphysical. As a result we take care to exclude initializations where at least in one location, the local perturbation h(x, t) is close to -1 or smaller.

Gaussian random process theory for near-monochromatic waves

In the regime where h 1 the Navier-Stokes equations may be linearized. After a phase of linear growth, we expect non linear effects to take over : Refer to [ ] Chapter for a good introduction.

: Even in the absence of this cutoff the strong damping of the linear theory would still annihilate the large k modes thereby giving us identical results. and the thread to pinch wherever a minimum of radius has developed. A simple model accounting for the linear growth is

h t Lh, ( . )
where L is a linear operator on the function h(x, t).

In Fourier space the equation becomes ĥk,t s(k) ĥk , ( . )

where s(k) has a maximum s m s(k m ) at the optimal wavelength k m 2π/λ m of the Rayleigh-Plateau-Savart instability, with λ m ∼ 9.02a.

As a result one can expand

s(k) s m -C(k -k m ) 2 + O((k -k m ) 3 ), ( .
)
where s m , C and k m are O(1) constants in dimensionless variables. The initial condition is a low-pass filtered noisy shape. We use the theory of Gaussian processes . We represent the initial solution as

h(x, 0) ∫ ∞ -∞ A(k)e ikx dk, ( . )
where k M k m is the cutoff defined above in the numerical section. The complex amplitudes A(k) form also a Gaussian process [ ] with

A(k)A(k ) 1 2 E(k)δ(k -k ) ( . )
where E(k) is the power spectrum of the initial conditions and • is the averaging operator. The initial condition is a random function characterised by its variance v 0 , which can be shown using the Wiener-Khinchin theorem to be

v 0 h(x, 0) 2 ∫ ∞ 0 E(k)dk. ( . )
In the numerics we use a low pass filter with E(k) v 0 /(2k M )

for k < k M and E(k) 0 otherwise. From ( . ) one has A(k, t) A(k) exp[s m t -C(kk m ) 2 t] and at finite time t one has the shape h(x, t) ∫ ∞ -∞ A(k)e s m t-C(k-k m ) 2 t e ikx dk ( . )

which can be rewritten as h(x, t) e s m t+ik m x ∫ ∞ -∞ B (κ)e -Cκ 2 t e iκx dκ, ( . )

where B(κ) A(k m + κ) A(-k m -κ). The integral on the RHS is a slowly varying modulation of the monochromatic wave e ik m x . The long lengthscale for this slow variation is -1 (Ct) 1/2 and the long space variable is X -1 x, with the corresponding small wavenumber K κ so that the slow modulation can be rewritten in slow variables as h(x, t) e s m t+ik m x H(X) + c.c., ( . ) where c.c. denotes the complex conjugate and the complex amplitude is

H(X) ∫ ∞ -∞
B(K)e -K 2 e iKX dK, ( . )

where B(K) B (κ). The slow modulation can be written as a product of a slowly varying amplitude and the exponential of a slow phase H(X) ρ(X) exp[iΦ(X)]. The throats of the ligament are located at the minima of h(x, t) that is at places where the total phase φ(x) k m x+Φ(X) obeys the discrete condition φ(x n ) 2nπ. The problem is translationally invariant so without loss of generality the origin can be picked at a minimum where Φ(X) 0 so that the next minimum is at x 1 such that φ(x 1 ) 2π. Thus we obtain k m x 1 + Φ( x 1 ) 2π, ( . ) which expresses that x 1 is a small perturbation of the optimal Rayleigh-Plateau-Savart wavelength x R 2π/k m and

x 1 x R + Φ( x 1 )/k m .

( . )

Since the variable X is slow at we are looking at small distances X 1 ∼ x 1 in the slow variables one has ρ(X 1 ) ρ 0 and Φ(X) αX where the amplitude ρ 0 and the slope α are random variables. Hence we use the approximation

x 1 λ m -2πα , ( . ) ( .
) B(K)B(K )e -K 2 -K 2 K KdK, ( . )

ρ[1 + iαX]
By the central limit theorem, since B(K)) is a Gaussian process, its integrals ρ and α are Gaussian variables.

To compute their variances we use

B(K)B(K ) E(k m ) 2 δ(K -K ). ( . )
Besides, from ( . )

ρ v 0 /2. ( . )
From all the aforementioned relations

ρ 2 ρ 2 + ρ 2 √ π 2 √ 2 E(k m ) , ( . )
and

ρ 2 α 2 √ π 8 √ 2 E(k m ) . ( . ) Since ρ 2 α 2 ρ 2 α 2 , ( . )
one obtains Var(α) 1 4 and thus from ( . ) Var(x 1 ) π 2 2 . The liquid mass in the interval of lenght x 1 will be captured in a breaking droplet once the nonlinear effects set in, at time t NL . This will happen once the amplitude has reached an order level that is when t t NL with t NL ∼ 1 2s m ln v 0 .

( . )

Once the breakup happens, the volume of the droplet captured between the throats is V e πx 1 , where the equivalent diameter of the resulting drops is πd 3 /6 πx 1 and d 3 6x 1 so that the variable d 3 is Gaussian with standard deviation The mechanism described thus far is surface-tension dominated. Another mode for breakup arises when the initial distribution reaches at some points the nonlinear regime. This occurs randomly along the ligament length and is thus a Poisson process with probability of occurence per unit length λ exp[-1/(2v 2 0 )]/l c , ( . )

where l c ∼ 1/k N is the correlation length of the initial condition. There is an exponential distribution of breaking lengths x 1 and droplet diameters d with rate λ , thus resulting in a distribution of diameters . . The d 3 -Gaussian does not seem to be the best fit, although the mode of the distribution seems to be well approximated by it. Revisiting the description of the tail of the drop diameter distributions from the previous section, the value of n in the function y A • exp[-nx m ] is related to the standard deviation of the d 3 -Gaussian (( . )). We find that the theoretical predictions of n differ significantly from what we obtain through our numerical results corresponding to the moderately perturbed ligament ensemble. In addition, we do not observe the transition from a e -d 6 to a e -d 3 scaling as predicted in the theory, at least in the few ensembles we have analysed thus far.

We are currently in the process of applying successive refinements to the theoretical development, as well as comparing the dispersion in the droplet sizes resulting from the disintegration of large enough ensembles of weakly corrugated ligaments (the d 3 -Gaussian prediction is valid only for very small amplitudes).

Conclusions & Perspectives

In this thesis, numerical investigations pertaining to two different aspects of liquid-gas interfacial flows with marked density contrasts were conducted using advanced geometrical reconstruction based Volume-of-Fluid methods. The first aspect involved understanding and eventually suppressing certain types of instabilities that plague standard numerical models due to the cascading growth of discretization errors that are particularly rampant at low to moderate numerical resolution. On the other hand, the second aspect of the investigation deals with a statistical description of drop sizes, where the drops were generated by carrying out large ensembles of well resolved direct numerical simulations of the capillary-induced breakup of randomly corrugated ligaments.

In the first part of the thesis, a detailed exposition into the design and development of our numerical methods is provided, where the principle of maintaining consistency between the discrete transport of the discontinuous mass and momentum fields was implemented using a conservative discretization of the two phase Navier-Stokes equations on staggered Cartesian grids. The shifted fractions and sub-grid methods were developed specifically in order to address the challenge of consistent advection of volume fraction, mass and momentum on staggered grids, coupled with the extension of the conservative direction-split mass transport algorithm of Weymouth and Yue [ ] to the direction-split transport of momentum. Several performance aspects of the aforementioned methods were presented using standard benchmark cases such as the behavior of spurious currents in static and moving droplets, and the propagation of damped capillary waves. Each of these tests were conducted for density contrasts spanning several orders of magnitude, including a realistic air-water configuration in the case of capillary waves. Both of our consistent methods demonstrate moderate to significant levels of improvement in accuracy when compared to the standard version of our method that is not consistent with respect to mass-momentum transport. A second-order error convergence rate is observed in the case of the static droplet and the capillary wave, with the rate reverting to first-order in the case of droplet advection. Finally, the robustness and stability of the present method were demonstrated using the case of a raindrop falling in air, which involves complex interactions of capillary forces, viscosity and inertia combined with the difficulties associated with flows involving marked density contrasts. A key feature of our class of consistent methods is the ability to circumvent or suppress the un-physical interfacial deformations that are rampant in the case of the standard non-consistent method, which typically results in the "artificial" atomization of the raindrop . Our methods are able to deliver relatively good estimates of the flow features especially at low to moderate resolutions, as evidenced by the smooth evolution of the droplet kinetic energy and the moments of inertia describing the average shape. Furthermore, both methods are able to reproduce the physically consistent early phase (almost constant) acceleration of the droplet towards the final statistical steady state.

Many of the complications in our methods are due to the staggered arrangement of the primary variables, thus opening a possible avenue for exploration of different grid configurations when it comes to consistent transport. Another perspective is an assessment of the computational speed-up of our methods compared to standard non-consistent methods which necessitate significantly higher spatial resolutions in order to ensure numerical stability. An interesting possibility exists in the context of the sub-grid method, where volume fraction information from the doubly refined grid could be used to obtain more accurate estimations of the interface curvature defined on the coarse grid level. An interesting test would be to verify if both the consistent and non-consistent methods converge to identical solutions at sufficiently large mesh resolution, owing to the fact that they are based on identical governing equations and differ mainly due to the choice of discretization scheme. In view of the pervasiveness of shear-driven primary instabilities in a wide range of fragmentation phenomena, a long term perspective would be to thoroughly evaluate and subsequently improve the accuracy of our consistent methods in comparison to the theoretical growth rate predictions (linear theory) of the Kelvin-Helmholtz instability for flows involving marked density contrasts, in a manner akin to our recent study [ ].

Moving onto the second part of the thesis, a thorough description of the methodology behind the generation of slender periodic ligaments is provided, where the ligament surfaces are created via superposition of random overlapping waves. A numerical setup was developed that enables precise quantitative control over the statistical properties of the randomly corrugated initial surface profiles, therefore setting up the stage for obtaining controllable, reliable and reproducible ensembles of droplets that are generated via disintegration of the original ligaments. High fidelity simulations of ligaments corresponding to air-water systems were carried out, delineating the influence of important factors such as the length scale, time scale, aspect ratio and strength of the initial corrugations on the subsequent set of drops formed. In the case of moderately perturbed millimeter scale water ligaments in air, the diameters of the resulting drops were found to be (approximately) symmetrically distributed around the mean size, where the overall shape is relatively well described by a Gaussian probability density function devoid of any free parameters. The shape of the diameter distributions is fairly robust, displaying negligible sensitivity to the choice of bin-width used to generate the histograms. In addition, the distributions converge towards a fixed shape as we increase the sample size N, where the error (defined with respect to the largest sample) scales as N -1/2 , therefore strongly suggesting the absence of any hidden correlations between the individual realizations constituting the sample. Amongst the initial conditions that influence the shape of the diameter distribution, the strength of the corrugations acts as a bifurcation parameter, changing the bimodal nature of the distribution shape for weakly perturbed ligaments to a unimodal one as we go towards stronger initial perturbations. Interestingly, in the case of moderately perturbed ligaments, the tail of the diameter distribution corresponding to the large drops is best described by an exponential scaling of the form e -d 6 even though the overall distribution seems to scale as e -d 2 . This scaling result was found to be independent of whether uncertainty in the bin heights was included or not, and further verified by reconstructing the histograms on "log(log) vs. log" scales. An interesting perspective is to compare the quality of the fit associated with the exponential scaling to alternate forms such as log-normal, power law etc. Finally, a Gaussian random process theory for near-monochromatic waves was developed in an effort to explain the origins of the newfound e -d 6 scaling. In the surface tension dominated regime associated to weakly perturbed initial ligament shapes, the theory predicts a resulting drop diameter distribution in the form of a d 3 -Gaussian, as well as a functional relationship between the dispersion in the size of the drops and the variance of the initial perturbations used to create the corrugated ligament surface. The theoretical predictions significantly differ from the observations derived from our numerical experiments corresponding to moderately perturbed ligament shapes, therefore opening the perspective of verifying the predictions within its scope of validity i.e. adequately large ensembles of weakly perturbed ligaments.

As we have devised a numerical setup that allows us to efficiently generate large numbers of drops using high-fidelity simulations, an avenue for further investigation lies in the exploration of ligament ensembles corresponding to larger aspect ratios, centimeter length scales and a broader range of initial corrugation strength. Such an undertaking may allow us to quantify the variation in the mean drop sizes as a function of the point in the parameter space Φ, opening up the possibility that the combination of size distributions characterized by different points in Φ could lead to more complicated and asymmetric probability distributions, thereby providing an alternate interpretation to the skewness observed in the drop size distributions typical of more complex fragmentation processes. On a final note, towards the objective of rendering our ligament model more representative of realistic fragmentation scenarios, a long term perspective would be to include the added complexity of inertial stretching as part of our initial conditions (We > 1).

.

  Liquid fragmentation triggered by the growth of perforations in thin liquid sheets. (i) Secondary atomization (bag mode) of a drop in a crossflow, driven by the rapid capillary expansion of the hole. Images reproduced from Opfer et al. [ ]. (ii) Effervescent atomization of expanding thin liquid sheets driven by the expansion of multiple perforations, image reproduced from Dombrowski and Fraser [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical simulations of a 3mm raindrop falling in air under gravitational

  . . . . . . . . . . . . . . . . . . . . . . . . . Decay of a liquid jet into droplets driven by the growth of capillary instabilities corresponding to different excitation frequencies. Image reproduced from Rutland and Jameson [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Examples of steady and unsteady liquid fragmentation, demonstrating the disintegration of the radially expanding liquid sheets driven by the capillary deceleration of the rims forming the edges of the sheet. (i) Steady fragmentation of radially expanding sheets driven by jet impact, image reproduced from Bremond et al. [ ]. (ii) Unsteady fragmentation of liquid sheets following drop impact, image reproduced from Wang and Bourouiba [ ]. . . . . . . . . . Experimental and numerical investigations of liquid atomization in which the primary stages of topological change are driven by shear instabilities. (i) Liquid jet disintegration by a high speed coaxial gas flow, image reproduced from Lasheras and Hopfinger [ ]. (ii) Atomization of a two-phase mixing layer between parallel liquid and gas streams, image reproduced from Ling et al. [ ]. . Exlicit definition of the interface location using the volume-of-fluid approach. These images are reproduced with permission from Scardovelli and Zaleski [ ]. (a) The exact discrete representation of a circular arc on a regular Cartesian grid using the colour function field (volume fraction). (b) The piecewise linear (PLIC) approximation to the smooth circular arc shown in (a), which entails second-order spatial accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Lagrangian transport of the interface segments using the CIAM scheme, the image is reproduced with permission from Gueyffier et al. [ ]. A D schemetic of the geometric calculation of the fluxes of the volume fraction field is shown, for the advection substep along the horizontal direction. The central cell (i, j, k) undergoes a net compression during this substep. The fluxes Φ -, Φ 0 and Φ + are the volumes under the advected interface segments advected by the interpolated velocity field, intersected by the (i, j, k) cell boundaries. . . . . .

  of the Poisson-like equation . . . . . . . . . . . . . . . . . . .

  Schematic of the static droplet of dense fluid surrounded by a quiescent medium of lighter fluid. A 40 × 40 grid is employed to spatially discretize the domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . STD Decay of normalized spurious currents as a function of viscous dissipation

  . . . . . . . . . . . . . . . . . . . . . . . . . . Scaling of the maximum error norm as a function of Weber (La 12000, ρ l /ρ g 1000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Scaling of the maximum error norm as a function of Laplace (We 0.4, ρ l /ρ g 1000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Schematic of the initially perturbed planar interface separating two immiscible fluids of different densities and viscosities. A spatial resolution of 32 × 96 is used for spatial discretization (compared to 64 × 192 in Popinet [ ])

  . The time stamps are indicative of the moments at which the interface is the most deformed, and do not necessarily correspond to the moment at which the code crashes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temporal evolution of quantities of interest to evaluate the performance of our consistent scheme based on the shifted fractions (MSHIFT) method, for different spatial resolutions. (a) Kinetic energy relative to the droplet center-ofmass as defined in ( . ). (b) Moment of inertia of the droplet along the flow (X) direction. (c) and (d) Moment of inertia of the droplet along the directions perpendicular to flow (Y,Z), evolution of I y y seems to be more or less identical to I zz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temporal evolution of quantities of interest to evaluate the performance of our consistent 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 ( . ). (b) Moment of inertia of the droplet along the flow (X) direction. (c) and (d) Moment of inertia of the droplet along the directions perpendicular to flow (Y,Z), evolution of I y y seems to be more or less identical to I zz . . . . . . .

  (c) Dependence of the mass conservation properties on the droplet resolution. The WY-QUICK combination displays better results for all except the lowest resolution. (d) Error estimation of the droplet acceleration in the frame of reference of the static box The corresponding accelerations are plotted in the inset of fig. . . . . Fragmentation of stretched liquid (Newtonian) ligaments formed immediately after retraction of a capillary tube. Image reproduced from Marmottant & Villermaux [ ]. The complex rearrangement of the liquid volumes inside the ligament plays a key role in its subsequent disintegration into droplet of various sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Representation of the liquid volumes in an isolated ligament as 'blobs' of sizes d matching the corresponding the local thicknesses, just before the destabilization into droplets. Image reproduced from Villermaux [ ]. . . . . . . . . . . . . . . Successive stages in the breakup of a ligament into drops. The system shown here is that of a perforated liquid sheet, where the growth of such perforations due to capillary retraction of the liquid rims lead to the formation of networks of connected ligaments. Image reproduced from Lhuissier and Villermaux [ ]. . A schematic of the octree grid used in Basilisk, and its corresponding logical structure (quadtree in D). Image reproduced from Agbaglah et al. [ ]. . . . . A demonstration of multiscale capabilties of Basilisk, wherein the mesh dynamically adapts in order to resolve the smaller scale features of the interface. Image reproduced from Agbaglah et al. [ ]. . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over intermediate time scales. The top row representssimulations carried out at a "maxlevel" of 11, the middle and bottom rows represents those carried out at "maxlevels" of 10 and 9 respectively. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the level of numerical resolution (from levels 4 to 11). The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. As one goes from the middle to the top row, there are no significant distinctions between the interface shapes (black curves), especially at long time scales such as T 20. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dynamics of the ligament destabilization, subsequent breakup and eventual coalescences over short time scales. The physical parameters correspond to those defined in ( . ). The top row (light blue) represents the geometrical reconstructions of the interface shape for the axisymmetric D simulations, and the bottom row (deep blue) represents those corresponding to the full D simulations. The discrepancy between the full D and axisymmetric D cases in terms of distances between the neighbouring structures is merely an effect of the distortion of lengths in the full D case, due to the use of a single point perspective (non-parallel) projection in order to visualize the interface shapes. . Dynamics of the ligament destabilization, subsequent breakup and eventual coalescences over long time scales. The physical parameters correspond to those defined in ( . ). The top row (light blue) represents the geometrical reconstructions of the interface shape for the axisymmetric D simulations, and the bottom row (deep blue) represents those corresponding to the full D simulations. The discrepancy between the full D and axisymmetric D cases in terms of distances between the neighbouring structures is merely an effect of the distortion of lengths in the full D case, due to the use of a single point perspective (non-parallel) projection in order to visualize the interface shapes. . Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over short time scales. Simulations of identical ligament surface profiles subject to different initial corrugation amplitudes, with the smallest amplitude represented by the top row and the largest amplitude by the bottom row. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the black contours. . . Dynamics of the ligament destabilization, subsequent breakup and eventual coalescences over long time scales. Simulations of identical ligament surface profiles subject to different initial corrugation amplitudes, with the smallest amplitude represented by the top row and the largest amplitude by the bottom row. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the black contours. . . . . . . . . . . . . . . Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over intermediate time scales. Simulations of identical ligament surface profiles at different length scales in accordance with the Ohnesorge numbers that characterize them, with the smallest length scale represented by the top row (Oh 0.1) and the largest length scale by the bottom row (Oh 0.001). The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the white contours. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dynamics of the ligament destabilization, subsequent breakup and eventual coalescences over intermediate time scales. Simulations for identical ligament surface profiles corresponding to different aspect ratios (having identical radii), with the shorter ligament represented by the top row and the longer one by the bottom row. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the black contours. . . Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over intermediate time scales. Simulations of ligament surface profiles that differ due to the different set of random waves used to generate them, with the profile corresponding to the lowest cut-off frequency represented by the top row and the highest cut-off by the bottom row. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the black contours. . . . . . . . . . . . . . . The shaded area A represented by the dotted and dashed lines is used to estimate the diameter of the droplet in question. The resulting diameter is computed as D 4A/π, and the corresponding volume is given by V πD 3 /6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Probability distribution functions of the droplet diameter at time T 15. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The diameters are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gaussian functions are characterized by the mean and variance of the original dataset, hence they are plotted alongside the histograms instead of being fitted to the bin heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Probability distribution functions of the droplet diameter at time T 30. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The diameters are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gaussian functions are characterized by the mean and variance of the original dataset, hence they are plotted alongside the histograms instead of being fitted to the bin heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Probability distribution functions of the droplet diameter at time T 30. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The diameters are normalized by the mean of the corresponding sample. The distribution is generated using bins of uniform size, represented in green. The Gaussian and Log-Normal functions are characterized by the mean and variance of the original dataset, therefore are plotted alongside the histogram, whereas, the Gamma function is fitted to the bin heights, with n 11.63 being the value corresponding to the best fit. . . . . . . . . . . . . . . . . . . . . . . . Probability distribution functions of the droplet volume at time T 30. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The diameters are normalized by the mean of the corresponding sample. The distribution is generated using bins of uniform size, represented in green. The Gaussian and Log-Normal functions are characterized by the mean and variance of the original dataset, therefore are plotted alongside the histogram, whereas, the Gamma function is fitted to the bin heights, with n 1.57 being the value corresponding to the best fit. . . . . . . . . . . . . . . . . . . . . . . . Probability distribution functions of the droplet volume at time T 15. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The volumes are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gamma functions are characterized by the parameter n, the value of which is determined by that which provides the best (least-squares) fit to the corresponding bin heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Probability distribution functions of the droplet volume at time T 30. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The volumes are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gamma functions are characterized by the parameter n, the value of which is determined by that which provides the best (least-squares) fit to the corresponding bin heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The fate of two ligaments whose initial surfaces are created using the same set of random overlapping waves, but differ by the strength of the perturbations. The ligament on the left is part of ensemble Φ 1 and the one on the right is part of Φ 2 . The droplets enclosed in the boxes (dashed lines) are characteristic of the distintegration of weakly perturbed ligaments, wherein their sizes are smaller than the typical size at least by a factor of 2. . . . . . . . . . . . . . . . . . . . . . Temporal evolution of the probability distribution functions of the droplet diameter. The diameter is rescaled by the initial width W of the ligaments. The ensemble is characterized by Φ 1 ≡ Oh 10 -2 , K 2π, ε 0.4, Λ 50 , thus consisting of weakly perturbed initial ligaments shapes. . . . . . . . . . . . . . . Temporal evolution of the probability distribution functions of the droplet diameter. The diameter is rescaled by the initial width W of the ligaments. The ensemble is characterized by Φ 2 ≡ Oh 10 -2 , K 2π, ε 1.6, Λ 50 , thus consisting of strongly perturbed initial ligaments shapes. . . . . . . . . . . . . . Comparison between the temporal evolution of the probability distribution functions of the droplet diameter, corresponding to the strongly (ε 0.4) and weakly (ε 1.6) perturbed ligament ensembles. The diameter is rescaled by the width of the initial ligament. The ensembles are characterized by Φ {1,2} ≡ Oh 10 -2 , K 2π, ε {0.4, 1.6}, Λ 50 . The thick dashed lines correspond to the mean of the respective distribution, and the thinner dashed lines located on both sides of the thicker line corresponding to the respective interquartile ranges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A representation of the moving window of points which is used in our search for the optimal parameters. The window starts from the extreme tail end of the distribution, extending towards the peak of the curve by including more bin heights on the way. In this particular figure, the window contains 6 points. The minimum size of the window should be greater than the number of free parameters in the fitting function. . . . . . . . . . . . . . . . . . . . . . . . . . . The variation in the optimal parameter m as a function of the number of points in the window (see . ). The least-squares procedure is performed on the candidate function in the exponential ( . ) form. The error bars represent the standard deviation of the parameter estimate obtained from the leastsquares routine, where the value enclosed in the red box corresponds to the minimum error. The figure (a) on top corresponds to the least-squares procedure where the uncertainty in the bin heights is taken into account, whereas the bottom figure (b) corresponds to those carried out with zero uncertainty in the underlying distribution. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes. . . . . . . . . . . . . . . . . . . . The tail of the distribution fitted with the function in the exponential form ( . ), using the optimal parameter set that corresponds to the lowest standard deviation in the parameter estimate. The uncertainty in the bin heights of the underlying distribution is not taken into account while conducting the least-squares procedure. The bottom figure is a closeup of the region in the vicinity of the tail.

Figure . :

 . Figure . : Liquid fragmentation triggered by the growth of perforations in thin liquid sheets. (i) Secondary atomization (bag mode) of a drop in a crossflow, driven by the rapid capillary expansion of the hole. Images reproduced from Opfer et al. [ ]. (ii) Effervescent atomization of expanding thin liquid sheets driven by the expansion of multiple perforations, image reproduced from Dombrowski and Fraser [ ].

  [ ]:Bussmann et al. ( ), 'Modeling high density ratio incompressible interfacial flows' were able to circumvent the issue surrounding staggered grids altogether by using a collocated arrangement in the context of hexahedral unstructured meshes, coupled with an unsplit Eulerian flux computation method. The study though makes no mention of any surface tension model.Level set based methods in the context of mass-momentum consistent transport were implemented first byRaessi and Pitsch [ ] [ ]:Raessi et al. (

  staggered-collocated approach to solve the consistency issue on polygonal unstructured meshes, complemented with a well-balanced surface tension model. Nevertheless, the VOF advection is based on algebraic transport, not to mention the use of a convolution based curvature computation, which is inherently not consistent. More recently, Nangia et al. [ ] : refer to Mirjalili et al. [ ].

Figure . :

 . Figure . : Decay of a liquid jet into droplets driven by the growth of capillary instabilities corresponding to different excitation frequencies. Image reproduced from Rutland and Jameson [ ].

  [ ]: Weymouth et al. ( ), 'Conservative volume-of-fluid method for free-surface simulations on cartesian-grids' to the direction-split transport of momentum.

Figure . :

 . Figure . : Examples of steady and unsteady liquid fragmentation, demonstrating the disintegration of the radially expanding liquid sheets driven by the capillary deceleration of the rims forming the edges of the sheet. (i) Steady fragmentation of radially expanding sheets driven by jet impact, image reproduced from Bremond et al. [ ]. (ii) Unsteady fragmentation of liquid sheets following drop impact, image reproduced from Wang and Bourouiba [ ].

Figure . :

 . Figure . : Experimental and numerical investigations of liquid atomization in which the primary stages of topological change are driven by shear instabilities. (i) Liquid jet disintegration by a high speed coaxial gas flow, image reproduced from Lasheras and Hopfinger [ ]. (ii) Atomization of a two-phase mixing layer between parallel liquid and gas streams, image reproduced from Ling et al. [ ].

  [ ]: Boeck et al. ( ), 'Ligament formation in sheared liquid-gas layers' , instead of the weighted arithmetic mean as in equation . . The two main (and most popular) approaches in the context of interface capturing schemes are the volume-of-fluid (VOF) method first developed by Hirt and Nichols [ ], and the level set class of methods pioneered by Osher and Sethian [ ]. The principal difference between the two approaches lies in the manner in which the Heaviside function χ is modeled in the discrete sense, either as a smooth differentiable field in the case of level sets, or as a sharp discontinuous field in the volumeof-fluid (VOF) context. Each class of methods has its own set of merits (and demerits) relative to each other. Generally speaking, volume-offluid based methods display superior mass conservation whereas in terms of interface curvature computation, level set based methods hold an advantage . A detailed exposition into the different classes of interfacial transport methods can be found in the seminal monograph by Tryggvason, Scardovelli and Zaleski [ ] : A comprehensive discussion about the different types of approximations to the interface Heaviside function can be found in Popinet [ ] : In-depth explanations into these numerical techniques can be found in [ -].

Figure . :

 . Figure . : Exlicit definition of the interface location using the volume-of-fluid approach. These images are reproduced with permission from Scardovelli and Zaleski [ ]. (a) The exact discrete representation of a circular arc on a regular Cartesian grid using the colour function field (volume fraction). (b) The piecewise linear (PLIC) approximation to the smooth circular arc shown in (a), which entails second-order spatial accuracy.

  This scheme was originally described in the work of Li [ ], 'CIAM' being an abbreviation for the French title 'Calcul d'interface affine par morceaux', which can be thought of as a straightforward Lagrangian transport of the interface Heaviside function. After the interface segments are reconstructed from the discrete colour function at the start of the time-step, the interfacial points are transported by the component of the velocity field corresponding to the direction of transport. A geometrical interpretation of the scheme is illustrated in figure . , reproduced from the seminal work of Gueyffier et al. [ ].

Figure . :

 . Figure . : Lagrangian transport of the interface segments using the CIAM scheme, the image is reproduced with permission from Gueyffier et al. [ ]. A D schemetic of the geometric calculation of the fluxes of the volume fraction field is shown, for the advection substep along the horizontal direction. The central cell (i, j, k) undergoes a net compression during this substep. The fluxes Φ -, Φ 0 and Φ + are the volumes under the advected interface segments advected by the interpolated velocity field, intersected by the (i, j, k) cell boundaries.

:

  The volume fraction field C is denoted as c in Fig... : The problem size of such elliptical partial differential equations is governed by the number of cells/points discretizing the velocity (or pressure) field. If we choose to solve the problem on a mesh with twice the resolution in D, that would lead to an fold increase in the problem size for the discrete pressure-Poisson problem.

  Figure . :A D illustration of the restriction operations involved between the different grid resolutions, where the subgrid grid boundaries are denoted by the dashed lines, and that of the overlapping coarse grid by solid lines. Restriction operations in the form of simple volume integrals are performed on the sub-grid volume fraction in order to compute the staggered and centered volume fraction fields at the coarse grid level. Restriction operators in the form of surface integrals are also performed in order to compute the volume fluxes at the boundaries of the coarse staggered control volumes. The red shaded areas correspond to the volume fractions, whereas the green shaded area corresponds to that of volume fluxes. The interface segments as depicted at the coarse grid level are just for purpose of illustration, and are not 'reconstructed'. Unlike the shifted fractions method, the interfacial reconstructions on the staggered volumes are not necessary, due to the fact that sub-grid volume fluxes are simply combined to compute the staggered fluxes.

Figure . :

 . Figure .: A bird's eye view of the subgrid method, highlighting the operations performed at each time step. and additionally ensures conservation of momentum. For sake of clarity, present a D case for the density on the grid i, j and horizontal velocity ũ on the staggered grid i + 1/2, j. The variables in red are all defined on the coarse level, and those in green are at the sub-grid level. 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 convery that the variables are centered on the staggered control volumes at the coarse grid level, whereas the 'bar' on top conveys that the variables are staggered with respect to the sub-grid. The starting variables are enclosed in elliptical borders, and the variables at the end of the directionsplit integrations are u * and ρ n+1 , where u * is subsequently fed into the RHS of the Poisson problem.

  The restriction operators work in an identical manner as in Rudman ([ ]), which are nothing but simple volume integrals of the sub-grid field mapped onto the domains corresponding to the coarse grid control volumes.In Fig. . , we demonstrate the different mappings between the two grid levels, which essentially boils down to volume integrals for the centered fields, and surface integrals for their corresponding fluxes. A description of the different functions performed by the restriction operator is as follows : Fine Grid -→ Coarse Staggered Grid Density & Momentum The sub-grid volume fraction field is restricted to obtain a staggered density field at the coarse grid level (see Fig. . ), which on combining with the appropriate component of the velocity field produces the staggered momentum field. Fluxes The geometric fluxes from the sub-grid volume fraction field are restricted in order to obtain coarse grid fluxes for the corresponding staggered volumes, which on combining with relation ( . ) gives us the corresponding momentum fluxes.

[

Figure . :

 . Figure . : Schematic of the static droplet of dense fluid surrounded by a quiescent medium of lighter fluid. A 40 × 40 grid is employed to spatially discretize the domain.

Figure . :

 . Figure .: STD Decay of normalized spurious currents as a function of viscous dissipation time-scales for different density-ratios and Laplace numbers. The currents seem to initially decay quickly for all higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For combinations of large ρ l /ρ g and large La, the spurious currents

Figure . :

 . Figure . : MSHIFT Decay of normalized spurious currents as a function of viscous dissipation time-scales for different density-ratios and Laplace numbers.The currents seem to initially decay quickly for all higher density-ratios, and relax to the numerical equilibrium curvature even within 0.2 • T µ . For combinations of large ρ l /ρ g and large La, the spurious currents

  , and subsequently employed in the comparative study of Abadie et al. [ ] [ ]: Abadie et al. ( ), 'On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks'

Figure . :

 . Figure . : Schematic of the droplet of dense fluid advected in a surrounding medium of lighter fluid. A 50 × 50 grid is employed to spatially discretize the domain, which is spatially periodic in the direction of droplet advection.

Figure . :

 . Figure . : STD Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers.The currents seem to hover around 10 -3 , with a larger Laplace number corresponding to a higher error for all density-ratios. We 0.4 for all the cases presented.

  Figures . to .depict the evolution of the root-mean-square (RMS) error of the velocity field in the moving frame of reference, as a function of different Laplace numbers, spanning over density ratios separated by orders of magnitude. The first figure ( . ) refers to simulations

Figure .

  Figure .: MSHIFT Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers. There seems to be no appreciable difference from the evolution seen in the case of STD (fig..). The currents seem to hover around 10 -3 , with a larger Laplace number corresponding to a higher error for all density-ratios. We 0.4 for all the cases presented.

Figure . :

 . Figure .: MSUB Time evolution of normalized spurious currents as a function of advection time-scales (T u ) for different combinations of density-ratio and Laplace numbers. In terms of the errors observed in STD and MSHIFT, 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 -4 , with larger Laplace numbers corresponding to larger errors for all density-ratios. We 0.4 for all the cases presented.

Figure . :Figure . :

 .. Figure . : Scaling of the maximum error norm as a function of Weber (La 12000, ρ l /ρ g 1000).

  Figure . : Schematic of the initially perturbed planar interface separating two immiscible fluids of different densities and viscosities. A spatial resolution of 32×96 is used for spatial discretization (compared to 64 × 192 in Popinet [ ]), with the width of the box corresponding to the size of the perturbed wavelength.

  [ ]:Lamb ( ), Hydrodynamics) or the asymptotic limit of vanishing viscosity(Prosperetti [ , ] [ ]: Prosperetti ( ), 'Free oscillations of drops and bubbles: the initial-value problem' [ ]: Prosperetti ( ), 'Motion of two superposed viscous fluids'). For our purposes, we use the configuration of the viscosity-damped capillary oscillations of a planar interface, as was first implemented and popularized byPopinet & Zaleski [ ].

Figure . :

 . Figure . : STD Time evolution of the amplitude of the planar interface undergoing damped capillary oscillations, comparing the solution obtained by our numerical method with the closed-from Prosperetti solution. More or less good agreement with theory is observed for all the densityratios tested.

Figure . :

 . Figure . : MSHIFT Time evolution of the amplitude of the planar interface undergoing damped capillary oscillations, comparing the solution obtained by our numerical method with the closed-from Prosperetti solution. Behavior is quite similar to STD, with good agreement with theory for all the density-ratios tested.

Figure . :

 . Figure . : MSUB Time evolution of the amplitude of the planar interface undergoing damped capillary oscillations, comparing the solution obtained by our numerical method with the closed-from Prosperetti solution. Slightly better agreement with theory when comparing to STD and MSHIFT, for all density-ratios tested.

Figure . :

 . Figure . : A D schematic of the numerical setup for the falling raindrop. A droplet of diameter D is placed at the center of a cubic domain of side L and L/D4. The liquid properties (ρ l , µ l ) correspond to that of water, and the gas properties (ρ g ,µ g ) correspond to that of air. We apply a uniform inflow velocity condition with U 0 (t) and an outflow velocity condition at the top which corresponds to zero normal gradient. Boundary conditions on the side walls correspond to those of impenetrable free slip (no shear stress).

Figure . :

 . Figure . :The origin of the pressure peak in the front of the droplet. (a) The profile of the pressure on the axis a few timesteps after initialisation with the standard, non-consistent method (red curve) and the consistent method based on the shifted fractions (MSHIFT) strategy (green curve). Much larger pressure gradients are present across the interface using the non-consistent method. (b) The pressure distribution immediately after the start of the simulation using the standard, non-consistent method. The pressure peak has not yet resulted in the formation of a dimple. In all figures the droplet resolution corresponds to D/h 16. The simulations are carried out with the CIAM advection method, in conjunction with the Superbee limiter.

  Figure . : (a) Comparison of the temporal evolution of droplet kinetic energy. The standard non-consistent method displays spikes in the kinetic energy that are approximately orders of magnitude larger than that of the consistent method, leading to rapid destabilization. (b)The pressure distribution immediately after the start of the simulation using the consistent method based on the shifted fractions (MSHIFT) strategy. In all figures the droplet resolution corresponds to D/h 16. The simulations are carried out with the CIAM advection method, in conjunction with the Superbee limiter.

Figure . :

 . Figure . : Flow field around the 3mm droplet with D/h 16 immediately after the start of the simulation with the consistent method (MSHIFT), demonstrating the contours of the norm of the vorticity field (black lines) . The D cross-section in these figures corresponds to the mid-plane slice along the Z axis, where the inflow is along the X axis and gravity opposite to it. (a) The velocity magnitude. As one can observe, the boundary layer is resolved by only -cells. (b) The velocity component in the Y direction, perpendicular to the flow. As the flow develops further, a marked separation of the boundary layers is observed with a more complex vortical region in the wake. The figures correspond to simulations carried out with the CIAM advection method, in conjunction with the Superbee limiter.

Figure . :

 . Figure . :A comparison of the temporal evolution between the standard non-consistent method (top row) and the consistent method based on the shifted fractions (MSHIFT) strategy (bottom row), using the combination of CIAM advection scheme coupled with the Superbee flux limiter. The flow is along the positive X direction, with gravity along the opposite direction. The red contour indicates the isosurface of the volume fraction field corresponding to a value of 0.5, whereas the black contours surrounding the drop represent isosurfaces of the magnitude of vorticity. The raindrop with the non-consistent method displays massive deformations leading to artificial breakup as a result of rapidly growing numerical instabilities. The droplet resolution for both methods is D/h 16 . The temporal evolution in case of the sub-grid based consistent method (MSUB) produces almost indistinguishable figures (bottom row), hence are not shown here.

Figure .

  Figure . :A juxtaposition of the different manifestations of 'artificial' atomization one encounters while using the standard method (STD) that does not involve consistency between the discrete transport of mass and momentum. The red contours indicate the isosurface of the volume fraction field corresponding to a value of 0.5, acting as a good proxy for the exact interfacial position. One can clearly observe that the un-physical fragmentation of the raindrop is symptomatic of the nonconsistent method, systematically across all combinations of flux limiters and advection schemes. The symbol "WY" in the legend corresponds to those run using the Weymouth-Yue advection scheme, and "CIAM" corresponds to the CIAM scheme. Brief descriptions of these advection schemes can be found in chapter . The implementation of the non-linear flux limiters i.e WENO, ENO, QUICK, Superbee, Verstappen correspond to that of well established methods developed to deal with hyperbolic conservation laws, for more details refer to the studies of Leveque [ ] and Sweby [ ]. The time stamps are indicative of the moments at which the interface is the most deformed, and do not necessarily correspond to the moment at which the code crashes.

Figure . :

 . Figure . : Temporal evolution of quantities of interest to evaluate the performance of our consistent scheme based on the shifted fractions (MSHIFT) method, for different spatial resolutions. (a) Kinetic energy relative to the droplet center-of-mass as defined in ( . ). (b) Moment of inertia of the droplet along the flow (X) direction. (c) and (d) Moment of inertia of the droplet along the directions perpendicular to flow (Y,Z), evolution of I y y seems to be more or less identical to I zz .

Figure . :

 . Figure . : Temporal evolution of quantities of interest to evaluate the performance of our consistent 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 ( . ). (b) Moment of inertia of the droplet along the flow (X) direction. (c) and (d) Moment of inertia of the droplet along the directions perpendicular to flow (Y,Z), evolution of I y y seems to be more or less identical to I zz .

Figure . :

 . Figure . :Comparison of droplet velocity as a function of time, for different droplet resolutions. The simulations were carried out with a consistent scheme based on the shifted fractions (MSHIFT) method, using the WY advection scheme with the QUICK limiter. The droplet velocities correspond to that of their respective center of masses. Inset : Convergence of the droplet acceleration as a function of its resolution, computed using the best linear fit over the temporal variation of their respective velocities. The error bars signify the asymptotic standard error (least-squares) corresponding to the linear fits.

Figure . :

 . Figure . : Relative change in the mass of the droplet as a function of time in plots (a), (b) and (c).The simulations are carried out with a resolution of D/h 16, for a total time of 0.1 milliseconds using the consistent scheme based on the shifted fractions (MSHIFT) method. The symbol "WY" in the legend corresponds to those run using the Weymouth-Yue advection scheme in combination with the QUICK limiter, and "CIAM" corresponds to the CIAM scheme in combination with the Superbee limiter. (a) Mass conservation properties of the scheme as a function of the Poisson solver tolerance is tested. For the WY-QUICK combination, using a stricter tolerance results in better mass conservation, whereas the conservation properties of the CIAM-Superbee combination seems to be independent of the tolerance. (b) Dependence of the mass conservation on the clipping parameter ( ) employed. The CIAM-Superbee combination again seems to be impervious to changes in clipping, whereas WY-QUICK seems to perform slightly better by lowering the parameter. (c) Dependence of the mass conservation properties on the droplet resolution. The WY-QUICK combination displays better results for all except the lowest resolution. (d) Error estimation of the droplet acceleration in the frame of reference of the static box The corresponding accelerations are plotted in the inset of fig. . .

:

  The accelerations correspond to the finest droplet resolution (D/∆x 64) in figures . and . for the MSHIFT and MSUB methods respectively.

Figure . :

 . Figure . : Fragmentation of stretched liquid (Newtonian) ligaments formed immediately after retraction of a capillary tube. Image reproduced from Marmottant & Villermaux [ ]. The complex rearrangement of the liquid volumes inside the ligament plays a key role in its subsequent disintegration into droplet of various sizes.

Figure . :

 . Figure . : Representation of the liquid volumes in an isolated ligament as 'blobs' of sizes d matching the corresponding the local thicknesses, just before the destabilization into droplets. Image reproduced from Villermaux [ ].

Figure . :

 . Figure . : Successive stages in the breakup of a ligament into drops. The system shown here is that of a perforated liquid sheet, where the growth of such perforations due to capillary retraction of the liquid rims lead to the formation of networks of connected ligaments. Image reproduced from Lhuissier and Villermaux [ ].

  model the underlying physical mechanisms that are responsible for the selection of droplet size. A model first proposed by Villermaux and coworkers

  Fig. . ). Another set of observations was presented by Lhuissier and Villermaux [ ], where the disintegration phase sets in starting from a transient configuration of the ligaments arranged in a spider-web like network (Fig. . ). A more recent experiment is that of Keshavarz et al. [ ], which is inspired by the capillary tube retractions of [ ] but directed towards the behaviour of non-Newtonian fluids.

Figure . :

 . Figure . : A schematic of the octree grid used in Basilisk, and its corresponding logical structure (quadtree in D). Image reproduced from Agbaglah et al. [ ].

Figure . :

 . Figure . : A demonstration of multiscale capabilties of Basilisk, wherein the mesh dynamically adapts in order to resolve the smaller scale features of the interface. Image reproduced from Agbaglah et al. [ ].

  Scardovelli and Zaleski [ ], or in a more recent exposition byMirjalili et al. [ ] on the same topic. The immiscible nature of the flows under consideration translates to the absence of mass transfer across the interface, thereby rendering the velocity field continuous across the interface at the incompressible limit. The study is primarily focused on isothermal air-water systems, with the material properties corresponding to that at degrees Celcius.NumericalPlatform : Basilisk We solve our set of governing equations using the open-source toolbox Basilisk [ ], which couples finite-volume discretization on octree meshes (Fig. . ) to solve the two-phase Navier-Stokes equations with surface tension. The use of adaptive mesh refinement (AMR) in conjunction with the octree structure allows an extremely efficient (compared to uniform grids) utilization of computational resources especially when dealing with the multiscale nature (see Fig. . ) of surface tension dominated flows, where the bulk of the computational power is diverted towards resolving regions in close proximity to the interface location. The solver utilizes a robust and accurate implementation of the CSF [ ] (Continuum Surface Force) surface tension model that governs the volumetric distribution of the capillary forces on the dynamically refined grid (refer to [ ] for more details). The numerical methods used in Basilisk have been extensively validated through studies conducted on a variety of free-surface flows, recent examples being on the topic of air entrainement through wave breaking [ , ], jets

Figure . :

 . Figure . :A schematic of the computational domain. An infinitely long and axisymmetric corrugated ligament of mean radius R is placed in a square box of side L. The bottom side of the box acts as the axis of symmetry for the ligament, with spatial periodicity along the horizontal direction. Inset : A close up view of the corrugated profile of the ligament, where the local radius is defined as the sum of the unperturbed (mean) radius R and the local perturbation (x). The material properties of the liquid and gas phases are denoted with the subscripts l and g respectively, which in our case corresponds to an air-water system with the surface tension coefficient σ.

  Figure . :Radial profile of the noisy and the subsequently low-pass filtered surface of the ligament, corresponding to one spatial period along its axis. Inset : Spectral densities of the initial noisy signal and the corresponding filtered signal. Only a handful of small wavenumbers are chosen to create the surface.

Figure . :

 . Figure . :The dispersion relation of theRayleigh-Plateau instability [ , ], expressing the growth rate of the small amplitude perturbations as a function of the reduced wavenumber ( . ). The vertical dashed line A represents the upper limit to the perturbation wavelength, corresponding to the domain size L. At the opposite end of the spectrum, E represents smallest wavelength limited by the grid size 1/2 M , where M is the maximum level of resolution for our octree grid. B corresponds to the wavenumber with the maximum growth based on linear theory, whereas C represents the limit beyond which all perturbations are stable. The number of discrete modes that fall within the unstable part of the spectrum (∆K) is a linear function of the ligament aspect ratio.

Figure . :

 . Figure . : Dynamically adapted octree meshes in the periphery of the interface location, with the refinement criteria based on limiting gradients of volume fraction and velocity. The interface is represented by the white contours, the colormap on the left half is based on the axial velocity component, whereas the one on the right corresponds to that of vorticity. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps.

Figure . :

 . Figure . :Illustration of the effects of our "droplet removal" feature, with the figure on the left corresponding to simulations that filter out structures that are resolved by cells or less, with the one on the right showing results with no droplet removal in place. The box (dashed line) in the right hand side figure displays the extremely small structures that consume a disproportionate amount of computational resources, without having any meaningful impact on the dynamics that we are interested in. The colormap on the left side of the axis of symmetry corresponds to the axial component of velocity, whereas the one on the right refers to the level of numerical resolution. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps.

Figure . :Figure . :Figure . :Figure . :Figure . :

 ..... Figure . : Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over intermediate time scales.The top row represents simulations carried out at a "maxlevel" of 11, the middle and bottom rows represents those carried out at "maxlevels" of 10 and 9 respectively. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the level of numerical resolution (from levels 4 to 11). The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. As one goes from the middle to the top row, there are no significant distinctions between the interface shapes (black curves), especially at long time scales such as T 20.

Figure . :Figure . :

 .. Figure . : Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over intermediate time scales. Simulations of identical ligament surface profiles at different length scales in accordance with the Ohnesorge numbers that characterize them, with the smallest length scale represented by the top row (Oh 0.1) and the largest length scale by the bottom row (Oh 0.001). The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the white contours.

Figure . :

 . Figure . : Dynamics of the ligament destabilization (closeup), subsequent breakup and eventual coalescences over intermediate time scales. Simulations of ligament surface profiles that differ due to the different set of random waves used to generate them, with the profile corresponding to the lowest cut-off frequency represented by the top row and the highest cut-off by the bottom row. The colormap is kept the same across all figures, where the left side of the axis of symmetry in each figure corresponds to the axial component of velocity (values from -1 to 1), whereas the one on the right refers to the vorticity (values from -1 to 1) perpendicular to the plane. The colors red and blue correspond to the higher and lower end values respectively, in case of both colormaps. The interface shapes are represented by the black contours.

AFigure . :

 . Figure . :The shaded area A represented by the dotted and dashed lines is used to estimate the diameter of the droplet in question. The resulting diameter is computed as D 4A/π, and the corresponding volume is given by V πD 3 /6.

Fig

  Fig. . shows a comparison of the different probability density functions applied to the largest ensemble of drops characterized by Φ 0 and T 30. The Gamma density function seems to be the best representation of the volume distributions. Considering the Gaussian (∼ e -x 2) nature of the diameter (D) distributions, a simple variable transformation leads us

Figure . :Figure . :

 .. Figure . : Probability distribution functions of the droplet diameter at time T 15.The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The diameters are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gaussian functions are characterized by the mean and variance of the original dataset, hence they are plotted alongside the histograms instead of being fitted to the bin heights.

Figure . :

 . Figure . : Probability distribution functions of the droplet volume at time T 30.The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The diameters are normalized by the mean of the corresponding sample. The distribution is generated using bins of uniform size, represented in green. The Gaussian and Log-Normal functions are characterized by the mean and variance of the original dataset, therefore are plotted alongside the histogram, whereas, the Gamma function is fitted to the bin heights, with n 1.57 being the value corresponding to the best fit.

Figure . :

 . Figure . :The fate of two ligaments whose initial surfaces are created using the same set of random overlapping waves, but differ by the strength of the perturbations. The ligament on the left is part of ensemble Φ 1 and the one on the right is part of Φ 2 . The droplets enclosed in the boxes (dashed lines) are characteristic of the distintegration of weakly perturbed ligaments, wherein their sizes are smaller than the typical size at least by a factor of 2.

Figure . :Figure . :Figure . :

 ... Figure . : Temporal evolution of the probability distribution functions of the droplet diameter. The diameter is rescaled by the initial width W of the ligaments. The ensemble is characterized by Φ 1 ≡ Oh 10 -2 , K 2π, ε 0.4, Λ 50 , thus consisting of weakly perturbed initial ligaments shapes.

  Figure . :A representation of the moving window of points which is used in our search for the optimal parameters. The window starts from the extreme tail end of the distribution, extending towards the peak of the curve by including more bin heights on the way. In this particular figure, the window contains 6 points. The minimum size of the window should be greater than the number of free parameters in the fitting function.

Figure . :Figure . :Figure . :

 ... Figure . : The variation in the optimal parameter m as a function of the number of points in the window (see . ). The least-squares procedure is performed on the candidate function in the exponential ( . ) form. The error bars represent the standard deviation of the parameter estimate obtained from the least-squares routine, where the value enclosed in the red box corresponds to the minimum error. The figure (a) on top corresponds to the least-squares procedure where the uncertainty in the bin heights is taken into account, whereas the bottom figure (b) corresponds to those carried out with zero uncertainty in the underlying distribution. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes.

Figure . :

 . Figure . : The variation in the optimal parameter m as a function of the number of points in the window (see . ). The least-squares procedure is performed on the candidate function in the logarithmic ( . ) form. The error bars represent the standard deviation of the parameter estimate obtained from the least-squares routine, where the value enclosed in the red box corresponds to the minimum error. The figure on top (a) corresponds to the least-squares procedure where the uncertainty in the bin heights is taken into account, whereas the bottom figure (b) corresponds to those carried out with zero uncertainty in the underlying distribution. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 and T 30, thus consisting of moderately perturbed initial ligaments shapes.

Figure . :

 . Figure . : Probability distribution functions of the droplet diameter at time T 30.The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The distribution is generated using bins of uniform size, represented in green. The Gaussian and Log-Normal functions are characterized by the mean and variance of the original dataset, therefore are plotted alongside the histogram, whereas, the Gamma and d 3 -Gaussian functions are fitted to the bin heights.

  P d,P (d) 3d 2 λe -λd 3 . ( . ) For large lengths and small rates λ, both types of breakup have small probability and the distributions combine. A crossover occurs for d c of small v 0 one has d m and thus d c ∼ . At d d c the Poisson distribution dominates and P d (d) ∼ P d,P . For d d c the d 3 -Gaussian distribution dominates. Therefore, for ligaments corresponding to weakly perturbed initial conditions (small corrugation amplitude), a exp[-d 6 ] scaling is expected by virtue of the d 3 -Gaussian ( . ) predicted by our theory. At the time of writing, we do not yet have a large enough ensemble of drops corresponding to such weakly corrugated ligaments, thus for the time being we compare the theoretical prediction to our ensemble of moderately perturbed ligaments (Φ 0 at T 30) as shown in Fig.

Table . :

 . Parameter values used in the simulation of a falling water droplet in air.

The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The volumes are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gamma functions are characterized by the parameter n, the value of which is determined by that which provides the best (least-squares) fit to the corresponding bin heights. The ensemble is characterized by Φ 0 ≡ Oh 10 -2 , K 2π, ε 1.0, Λ 50 . The volumes are normalized by the mean of the corresponding sample. The distributions are generated using datasets corresponding to four different sample sizes, including four different choices of (uniform) bin width. The Gamma functions are characterized by the parameter n, the value of which is determined by that which provides the best (least-squares) fit to the corresponding bin heights.