A robust interface method for drop formation and breakup simulation at high density ratio using an extrapolated liquid velocity

A two-phase ﬂow formulation for atomisation modelling is presented, with a Coupled Level Set/Volume Of Fluid (CLSVOF) technique adopted for interface-tracking. In order to achieve stable numerical solution at high density ratios, an extrapolated liquid velocity ﬁeld is constructed and used in discretisation of the momentum equations. Solution accuracy is also improved when this ﬁeld is also used in the scalar (VOF and Level Set) advection equations. A divergence-free algorithm is proposed to ensure satisfaction of the continuity condition for the extrapolated liquid velocity. The density and viscosity across the interface are treated sharply as a function of the Level Set to maintain the physical discontinuity. The developed method is shown to accurately predict drop formation in low Re liquid jets and the deformation and breakup morphology of a single droplet in uniform air ﬂow at different Weber numbers (from 3.4 to 96). The mechanism for droplet breakup is determined based on an analysis of the simulation results. The computed Rayleigh–Taylor instability wavelength extracted from the acceleration of the simulated liquid droplet agrees well with experimental measurements and theoretical analysis, conﬁrming that Rayleigh– Taylor instability dominates single drop breakup in the Weber number range studied. Finally, the inﬂu- ence of liquid viscosity on droplet breakup is numerically investigated; the critical Weber number sepa-rating deformation and breakup regimes is well predicted at different Ohnesor ge numbers in comparison with the experimental data.


Introduction
An accurate method for atomisation prediction is of utmost significance in many industrial applications. Substantial research has therefore been carried out to understand the atomisation process. State of the art reviews relevant to this important phenomenon have appeared recently, for example [1] and [2] . These reviews have outlined the complex numerical challenges introduced by the need to deal with accurate liquid/gas interface tracking, discontinuous fluid properties across the interface, and surface tension effects.
Three popular interface capturing methods have been proposed for flows involving a liquid/gas interface. The first was the Volume of Fluid (VOF) method, where the VOF function F is defined as the liquid volume fraction within any element of space. It was first proposed by Hirt and Nichols [3] , and further developed/applied in [4][5][6][7][8][9][10][11][12][13] . VOF can inherently conserve the liquid volume, but it needs significant numerical effort to extract the interface geometrical properties (normal and curvature) to 2nd order accuracy from F due to its discontinuous nature, especially in 3D [8,11,13,14] . An alternative choice for interface tracking, proposed by Sussman et al. [15] , is the Level Set (LS) method, with recent developments reported in [16][17][18] . Unlike F the LS function φ is a continuous variable; the interface is defined by the contour φ = 0 , and φ represents the signed distance from the interface, with φ > 0 liquid and φ < 0 gas. The LS method thus provides a numerically more convenient representation of the interface, making it straightforward to locate interface position and calculate interface normal and curvature. However, the LS method can induce considerable error of liquid mass as the calculation advances in time [19] , requiring special procedures to alleviate this problem [20][21][22] . To combine the advantages of VOF and LS methods into a single algorithm, a hybrid or coupled LS and VOF method (CLSVOF) has been proposed by Sussman and Puckett [23] , with applications to two-phase flow simulations in [24][25][26] .
Simulation of atomisation at high liquid/gas density ratios is very challenging as it is prone to numerical instabilities [2,27] . Since the non-conservative form of the governing equations are commonly solved for incompressible two-phase flows, large momentum discretisation errors can occur in computational cells in the interface vicinity when cell face flux interpolations are carried out. Momentum flux interpolation is not straightforward near an interface since both density and velocity gradient (due to the viscosity discontinuity) change sharply; conventional linear interpolation practices can lead to significant errors that cause numerical instability as the density ratio increases. Rudmann [28] , Raessi [29] , and Desjardins and Moureau [30] describe various techniques to solve this problem, in particular special interpolation practices for cell face density are suggested. It would be preferable if the sharp jump in density (and viscosity) across the interface were retained as far as possible, so a more promising method has been proposed by Sussman et al. [24] ; this introduces the concept of an extra liquid velocity field 'extrapolated' onto the gas phase side of the interface, and a modified version of this is incorporated in the present methodology.
Since it is computationally extremely expensive to carry out Direct Numerical Simulation of atomisation, a lot of effort has been put into the development of two-phase Large Eddy Simulation (LES). Additional complexities arise when applied to two-phase flows in connection with extra terms requiring sub-grid-scale (SGS) modelling. The surface tension term in the momentum equations is non-linear and hence filtering will lead to an additional SGS term; similarly, extra SGS fluxes will also appear due to the non-linear convection terms in the scalar equations which determine the interface behaviour (VOF and Level Set). In the vast majority of LES studies for two-phase flows so far, all specific SGS terms associated with the presence of an interface have been ignored, and this has been described by Gorokhovski and Herrmann [1] and Bianchi et al. [10] as quasi-DNS/LES as explained in the next section, this corresponds to an under-resolved DNS of interface tracking combined with an LES of the two-phase flow equations without any explicit inclusion of SGS surface tension modelling.
Different modelling approaches have been put forward in the last few years to address the extra SGS terms. The work of Toutant et al. [31,32] for example proposes a fundamentally different filtering approach for two-phase flows. This aims to account for the inevitable under-resolution of the interface in an LES calculation and leads to extra SGS terms modelled using a scale similarity principle. So far, however, this approach has only been analysed via explicit filtering of an a-priori DNS calculation to extract the extra terms introduced after LES filtering and compare these with proposed models; so far no actual LES calculations have been reported following this approach.
An alternative approach specifically relevant to surface tension SGS modelling has been proposed by Herrmann and Gorokhovski [33,34] and was recently investigated by Aniszewski et al. [35] . Calculation of the modelled surface tension SGS term requires evaluation of an expression which involves not only the filtered (resolved) surface normal n i , but also the instantaneous surface normal n i . To calculate n i requires the transport equations for the instantaneous interface-determining variables ( F and φ) to be solved using the instantaneous velocity field. To reconstruct this from the LES-filtered velocity field the Approximate Deconvolution Model (ADM) of Stolz et al. [36] was used. The test cases used in [35] to assess the importance of the surface tension SGS term and the accuracy of the proposed model had low density ratios ( ≤ 10), far removed from the typical (water/air) level of 800 of interest in the current work. Thus, it seems that SGS modelling fully extended to two-phase flow is still in its infancy, and, for this reason, the approach adopted in the present study was the quasi-DNS/LES approach.
In order to elucidate the physical atomisation mechanism, many experimental and numerical studies on the deformation/breakup of a single drop (often referred to as secondary atomisation) have been carried out [37] . The Weber number ( W e = ρ G U 2 G D 0 /σ ) which represents the ratio of the disintegrating aerodynamic force to the stabilising surface tension force is the most important characteristic parameter in single drop breakup. As We increases, several modes of behaviour are observed experimentally: a pure deformation mode, bag breakup, bag-stamen breakup, multimode breakup, sheet-thinning (or shear) breakup and shear-induced entrainment (or catastrophic breakup) [26,[37][38][39] . Theofanous et al. [40] and Zhao et al. [41] demonstrated that the Rayleigh-Taylor instability determines the drop breakup morphology at low Weber number ( We < 80) by comparing their theory with their own experimental results. The mechanism for droplet breakup has been more disputed for We higher than 80. The shear-stripping mechanism was proposed by Ranger and Nicholls [42] and was widely adopted in the last century [38] . In the shear-stripping mechanism, it is postulated that a liquid boundary layer is developed adjacent to the interface inside the drop under the action of shear stress from the gas flow. As the liquid boundary layer becomes unstable, liquid mass is stripped at the drop periphery. A lot of doubt has been cast on this shear-stripping mechanism by Liu and Reitz [43] , Lee and Reitz [44] , Guildenbecher et al. [37] and Theofanous et al. [39] [45] . The shear-stripping model suggests that the breakup mode should be a function of Re , which is contradictory to Liu and Reitz's experimental findings [43] . Thus, Liu and Reitz [43] proposed a sheetthinning breakup mechanism (for 80 < We < 350) which is consistent with their experimental results. In this sheet-thinning breakup mechanism, the droplet first deforms into a disc-like shape with the thickness growing thinner from the center to the edge; then the periphery of the flattened drop is bent in the direction of the flow due to its low inertia, forming a liquid sheet which disintegrates into ligaments and droplets. Since numerical methods can provide more flow details which can help understand the atomisation mechanism, many simulations of droplet breakup have been carried out [27,[46][47][48][49][50] . However, most of these published numerical studies to date have been limited to liquid/gas density ratios only of order 1-100. Since the majority of atomisation experiments are carried out in air at atmospheric pressure with high density liquids, experimental data are mainly available for density ratios a factor of 10 greater, and thus quantitative comparison between numerical modelling and experiment is quite rare. In the present simulations of droplet breakup, water is used as the liquid and air at atmospheric pressure is used as the gas, resulting in a high density ratio of 830. The drop breakup mechanism will then be analysed and the breakup mechanism proposed based on the simulation results.
The effect of liquid viscosity is to retard the drop deformation process and thus hinder breakup, and when viscous effects are significant, this can affect the critical Weber number We cr that separates breakup and deformation modes. In order to characterise this effect, the Ohnesor ge number ( Oh = μ L / ρ L D 0 σ ) may be introduced as the ratio of the liquid viscous force to the surface tension force. Empirical correlations between We cr and Oh have been proposed by Pilch and Erdman [51] ( W e cr = W e cr0 (1 + 1 . 077 Oh 1 . 6 ) ) and Gelfand [52] ( W e cr = W e cr0 (1 + 1 . 5 Oh 0 . 74 ) ) based on experimental data. Cohen [53] has also proposed a semi-empirical correlation based on analysis of energy transfer in secondary breakup ( W e cr = W e cr0 (1 + C Oh ) where C has a value between 1.0 and 1.8). The proposed empirical correlations differ significantly from each other due to inaccuracies in the experimental measurements. An initiation time -defined as the time required for a drop to deform beyond oblate ellipsoid shape -has been identified for a range of Oh values [37,51] , with several correlation functions proposed by Pilch and Erdman [51] , Hsiang and Faeth [54] and Gelfand et al. [55] . However, significant discrepancies may be observed between these correlations. Therefore, the opportunity is taken here to use the developed numerical technique to investigate the influence of liquid viscosity on drop deformation and breakup.
The present two-phase flow solver is developed using an existing multi-block structured mesh code for LES of single phase constant density turbulent flow [56][57][58] . In the following sections, the two-phase flow governing equations and numerical methods are first described. Validation test cases covering interface instability and drop formation in laminar liquid jets are presented next. Finally, the mechanism of droplet breakup and influence of liquid viscosity are examined using the developed numerical technique.

Formulation of two-phase flow governing equations
The philosophy for the current approach to two-phase flow simulation is to adopt the usual spatially filtered LES formulation [59] , with an overbar used to represent spatial filtering. Thus for any variable , its spatially filtered value is given by: where G is the filter kernel. In the present formulation a box filter is used [59] .
Since both liquid and gas are assumed to be incompressible and immiscible, the continuity equation is then everywhere the same as in single phase flows, and its filtered version reads: where U i is the instantaneous velocity. After filtering, the nonlinearity of the convection terms leads to the appearance of a residual or sub-grid-scale (SGS) stress tensor ( τ SGS i j as defined below) in the filtered momentum equations: where P is pressure, g i is gravitational acceleration, τ mol i j represents the molecular viscous stress. ˜ F ST i is the singular surface tension force located on the filtered interface which is constructed by the surface-averaging filter kernel proposed by Pitsch [60] for a flame front. A simple Smagorinsky eddy viscosity approach is used here, where is used to represent the filter width (taken here as the cube root of the local computational cell volume, with the value of the Smagorinsky constant C S set in all calculations to be 0.1); the full expressions for the diffusion terms become: S i j is the resolved strain rate tensor, with magnitude S . Since the interface is tracked explicitly in the current formulation, density and viscosity are treated sharply to maintain the physical discontinuity across the interface, and thus are set to be the properties of liquid/gas depending on the local value of the resolved LS variable ˜ φ ( ˜ φ is the LS representation of the filtered interface as detailed below): H( ˜ φ) is the Heaviside function; subscripts G and L indicate gas and liquid properties respectively. For the filtered momentum equations, it remains only to express the filtered surface tension in terms of resolved variables. As noted above, SGS surface tension modelling is so far rather immature and still under development; furthermore, it has not been well validated against experimental data particularly for the high liquid/gas density ratios of interest here. For these reasons, no SGS component of the surface tension force has been included, and the resolved surface tension force is computed directly from the morphology of the filtered interface: Here, σ is the surface tension coefficient. ˜ κ and ˜ n i are respectively the curvature and normal vector (pointing from the liquid phase into the gas) of the filtered interface.
Since both liquid and gas phases are viscous, the jump condition on the interface is (readers are referred to [61,62] for more details): Finally, the LES version of the scalar advection equations which determine the interface movement is derived. Oberlack et al. [63] showed that the classical Reynolds ensemble averaging and some LES SGS models violated the generalized scaling symmetry of the G-equation (which reduces to LS equation when the laminar burning velocity is set to be 0). A consistent formulation of the G-equation for the filtered flame front based on a new filtering technique was proposed by Pitsch [60,64] . Following Pitsch's procedure [60,64] , the LS equation which governs the evolution of the filtered interface is given as follows: where ˜ φ is the LS representation of the filtered interface rather than the filtered LS field, and is the signed distance to the filtered interface.
In order to reproduce a sharp interface by the VOF method in the LES formulation, the spatial filtering operation defined by Eq.
(1) can not be applied to the VOF field and VOF advection equation. To be consistent with the LS representation and keep a sharp interface, the resolved VOF field ˜ F in the LES is the liquid volume fraction determined by the filtered interface rather than the spatially filtered VOF field. ˜ F is evolved by the filtered velocity field, and the contribution of the SGS term is neglected, resulting in a scalar advection equation with the same form as the LS equation: Some further justification for neglect of SGS surface tension and SGS interface dynamics here was provided by the observation that the liquid flow is laminar in the jet/drop breakup simulations reported below and thus no SGS interface deformation is induced by any SGS velocity in the liquid field. Marmottant and Villermaux demonstrate in their experiments ( Fig. 42 in [65] ) that a smooth interface and axisymmetric deformation are observed for an injected laminar liquid jet while irregular interface distortions are observed right after a turbulent liquid jet exits the nozzle. This is also confirmed in the LES of a liquid jet in coaxial air flow by Xiao et al. [66] . In the liquid jet/drop deformation period, both experiments ( Fig. 1 in [67] , Fig. 7 in [41] ) and the current simulations give smooth axisymmetric deformation without any irregular interface distortions, implying that the velocity field in the liquid phase for the simulations presented here is laminar.
Although the gas flow is turbulent, the SGS gaseous eddies do not have enough energy to directly distort the interface as the liquid has a high inertia arising from its high density in comparison with the gas. LES of a liquid jet in air crossflow by Xiao et al. [68] shows that the large eddies in the incoming air flow do not cause any deformation of the liquid column due to the high liquid/gas density ratio. Therefore, the less energetic SGS eddies in the gas phase are unable to contribute to the interface deformation.
Since solution of the Level Set Eq. (10) does not guarantee satisfaction of the signed distance property ( ∇ ˜ φ = 1 ), a reinitialisation equation is solved: Re-initialisation is carried out for a pseudo-level set variable ϕ in pseudo-time τ . ϕ is initialised ( τ = 0 ) to equal the solution of the physical (resolved) level set Eq. (10) at time t, i.e., ϕ 0 = ϕ( In what follows, for simplicity the overbar indicating spatial filtering has been omitted, but all variables represent LES resolved quantities.

Numerical details
The present two-phase flow solver was developed using a Cartesian staggered mesh. For convenience the methodology and discretisation scheme is illustrated here in 2D for simplicity though it is actually implemented in 3D in the code, with details provided in [26] . The grid and variable arrangement are shown in Fig. 1 ; pressure, φ, F, ρ, and μ are located at cell centres; the velocity components are located at corresponding faces. u and v are the velocity components obtained after solution of the governing equations outlined in the previous section, while u L and v L are the components of a liquid velocity field constructed by an extrapolation and divergence-free approach as detailed below.
The implementation of the present CLSVOF method is described in [26,66,68,69] . A detailed description of the overall algorithm, including the way in which F and φ solutions are combined to deliver the transient dynamics of the interface geometry (interface reconstruction, and F / φ field evolution), has been outlined in full in [26] . In brief, 2nd order accurate operator split methods [24] are used in both F and φ equations. The interface normal vector is calculated from the level set gradient ( d φ) by a proper choice of forward difference ( d + ), central difference ( d c ), and backward difference ( d − ): when the level set is determined by one in- Then the interface is translated along the normal direction in order to ensure compatibility with both F and φ solutions within each cell. Numerical tests (see [26] ) have shown that, when combined with the CLSVOF approach, this allows good accuracy of interface tracking simultaneously with ensuring low liquid mass errors. The extrapolated liquid velocity is used also in the LS and VOF advection equations following [24] . Emphasis is placed on providing full details of the creation of the extrapolated velocity field, and its use in the discretisation of the governing equations in the following description.

Temporal discretisation
The LES code for single phase flow used as the starting point for the present work [56] had adopted the classical techniques of a centred 2nd order method for all spatial discretisation and a 2nd order Adams-Bashforth scheme for temporal advancement. In single phase flows, convection and diffusion terms are continuous in both time and space. However, in the present application these terms are discontinuous across the liquid/gas interface. Convection and diffusion terms will be discontinuous in time in any cell where the phase changes from gas to liquid (or vice-versa) during a time step. Since the Adams-Bashforth scheme is based on Taylor series expansion and requires continuity of convection and diffusion terms, this may cause local errors for two-phase flow. For caution's sake therefore, a simple 1st order forward Euler projection method was used for temporal discretisation of the two-phase flow equations. Since explicit time-marching was involved, and the associated CFL max constraint in LES calculations means the time step is usually very small (typically for the simulations presented below CFL max was O (0.1)), any associated error should be negligible, as was confirmed by a time-step sensitivity study for single droplet breakup [69] .
sion of the solution of momentum equations which includes convection, molecular and SGS diffusion, and gravitational terms (Note that surface tension is treated as part of the pressure term using the approach to be described in sub-Section 3.4 ): The intermediate velocity field U * i is updated using a pressure gradient to obtain the velocity field n + 1 : Since the velocity field at time step n + 1 must satisfy the continuity equation, a pressure Poisson equation may be derived by taking the divergence of the above equation, whose solution allows P n +1 to be calculated: P n +1 is then used in Eq. (14) to update the intermediate velocity field to establish U n +1 i .

Extrapolated liquid velocity and divergence-free approach
The philosophy of using an extrapolated liquid velocity field is based on the strong discontinuity of the velocity gradient across the interface observed in two-phase flow with high density and viscosity ratio. For the two-phase shear flow in equilibrium demonstrated in Fig. 1 , the high liquid/gas viscosity ratio ( O (100)) indicates that the velocity gradient in the gas is much larger than that in the liquid as the shear stress is the same across the interface. Therefore, the spatial discretisation of the governing equations in the vicinity of the interface must be carefully designed to tackle this discontinuity.
In the following, the discretisation of the momentum equation for u in x-momentum CVs is examined. When the level set value at a CV node is positive, this CV is referred to as a liquid CV, otherwise, it is referred to as a gas CV. For example, since φ i +1 / 2 , j > 0 , the x-momentum CV i +1 / 2 , j is a liquid CV, and it is treated as if it contains only liquid: the resolved velocity u i +1 / 2 , j thus represents a liquid velocity, and the density in this x-momentum CV is set to be liquid density (i.e. ρ i +1 / 2 , j = ρ L ). Similarly, due to φ i −1 / 2 , j > 0 , the x-momentum CV i −1 / 2 , j is considered as a gas CV: the resolved velocity u i −1 / 2 , j represents the gas velocity, and the density in this x-momentum CV is set to be gas density (i.e. When solving for the gas velocity u i −1 / 2 , j from the momentum equation in the gas CV i −1 / 2 , j , the gas momentum flux at the right face ρ G u i, j u i, j needs to be computed. When solving for the liquid velocity u i +1 / 2 , j from the momentum equation in the liquid CV i +1 / 2 , j , the liquid momentum flux at the left face ρ L u i, j u i, j needs to be computed. In this sense, the borderline between these two CVs is numerically treated as the two-phase interface, and u i, j should therefore represent the interface velocity in the numerical approach. As a consequence of the fact that the interface velocity is much closer to that in neighbouring liquid CVs than in neighbouring gas CVs, a good approximation to u i, j is to extrapolate the velocity in the neighbouring liquid CV to this point. For convenience, the extrapolated liquid velocity at the gas CV node u L i −1 / 2 , j (indicated by the red arrow labelled u L in Fig. 1 ) is first calculated, and then u i, j is computed from arith- This issue can also be explained in another way. Since the liquid/gas density ratio is large ( O (10 0 0)), any error in u i, j can result in a much larger error in the momentum flux term ρuu in the liquid phase than in the gas phase. Therefore, it is more important to find a proper value of u i, j so that the calculated convection term ( in the conventional discretisation approach would considerably overpredict the velocity gradient in the liquid CV because of the use of the gas velocity u i −1 / 2 , j , resulting in significant numerical error in momentum. Use of an extrapolated liquid velocity approach is a better option in this case as demonstrated below. To implement the extrapolated liquid velocity technique, a separate array U L i was created (with components u L and v L in 2D) as follows: i. U L i at liquid phase nodes ( φ > 0) was set equal to the momentum-equation-deduced velocity U i : ii. U L i at gas phase nodes ( φ < 0) was computed by outwards extrapolation along the interface normal from liquid into gas. For example, the liquid velocity component u L at gas phase nodes close to the interface was calculated by solving the following extrapolation equation to steady state: A forward Euler scheme was used for discretisation of this extrapolation equation in pseudo-time τ : where the pseudo time step

This equation is solved for 8 time steps
to create an extrapolated liquid velocity in a two-cell thick layer on the gas phase side of the interface. And a first order upwind scheme was used for spatial discretisation, e.g.: Since the extrapolated liquid velocity field will not automatically be divergence-free, the liquid velocity field calculated at the gas nodes must be corrected to satisfy the continuity equation. First, a liquid velocity source term in a typical cell ( i, j ) was computed: The extrapolated liquid velocity was then corrected by an upwind scheme: where A is the projected cell face area and the coefficients are given by: Eqs. (20) , (21) , and (22) were solved for sufficient time steps (typically 4) to create a continuity-satisfying extrapolated liquid velocity in a two-cell thick layer on the gas phase side of the interface.

Spatial discretisation for momentum equation
In general, a centered 2nd order approximation is followed using a classical finite-volume approach, leading to the need to evaluate cell face convective and diffusive fluxes. For the momentum equations the approach adopted deviates from the classical form in terms of how the liquid extrapolated velocity is used in the convective flux (this then also applies to the F and φ equations), and how the cell face density and effective viscosity are chosen in the diffusive flux. The practices adopted are illustrated here by considering a single convective and diffusive contribution as examples. Each momentum CV is treated as either liquid or gas depending on the sign of the Level Set at the node for that CV. For example, if the φ value at the u i −1 / 2 , j node of the x-momentum CV i −1 / 2 , j ( Fig. 1 ) is greater than zero, it is a treated as a liquid CV, the resolved velocity at the node represents a liquid velocity and the density for this CV is set to the liquid density; otherwise gas density is used, thus: Considering for illustration purposes a typical x-momentum convective term in a cell intersected by the interface, to reduce momentum errors, the extrapolated velocity is introduced into the discretisation following the rules indicated below: Similarly, considering a typical diffusion term in the xmomentum equation, special care must be taken to respect the discontinuous nature of density and effective viscosity across the gas/liquid interface: where, for example: μ e f f E is the effective viscosity at the east face of the x-momentum control volume i −1 / 2 , j . This face is co-incident with the pressure node ( i, j ). Thus, for most cells the value of μ eff calculated at ( i, j ) can be used for μ e f f E ; however, care must be taken for momentum nodes which lie close to the interface. The practice followed is: an effective viscosity μ e f f i, j is calculated and stored at the nodes of pressure CVs: Where S i, j and S L i, j are the magnitude of the resolved strain rate tensor at node ( i, j ) evaluated using momentum-equation-deduced velocity field and liquid velocity field respectively. Next, the eddy viscosity at face E of the momentum cell is determined. In Eq. (34) , velocities u i −1 / 2 , j and u i +1 / 2 , j are used to calculate the effective stress τ e f f xx i, j . When both x-momentum CVs i −1 / 2 , j and i +1 / 2 , j are liquid Readers are referred to [26] for calculation of other diffusion terms.

Surface tension term -Ghost Fluid method
In the current simulations, the interfacial pressure jump arising from surface tension is incorporated into the discretisation of the pressure gradient following the Ghost Fluid Method [70][71][72] . If the two-phase interface is located between nodes (i − 1 , j) and ( i, j ), the pressure gradient at face (i − 1 / 2 , j) is discretised as: Here, [ P ] σ is the pressure jump across the interface due to the surface tension: The curvature at the interface κ is calculated using linear interpolation: where κ i, j is calculated from derivatives of the Level Set function [16] .

Pressure Poisson equation
The Poisson Eq. (15) is discretised by integration over i, j (see Fig. 1 ), and has the form in 2D: The pressure Laplace operator is discretised via: Pressure gradients are discretised as described in sub-Section 3.4 , incorporating surface tension effects via the ghost fluid method. The density is treated sharply as in Eq. (28) . The source term is discretised as: ρ ∇P is continuous across the interface although there is a jump in ∇P due to the density discontinuity. The standard multigrid method for solving elliptic PDEs, in particular the bilinear interpolation operator, implicitly relies on continuity of ∇P . A more appropriate interpolation operator is one that can exploit the continuity of 1 ρ ∇P . Operator-induced interpolation is the optimum route to achieve this and was implemented in the Box multigrid method (BoxMG) by Dendy [73,74] . The BoxMG code is available (see [75] ), and its implementation into the code used here is given in [26] . For two-phase flows containing complex interfaces and strong discontinuity, combining the multigrid method with a preconditioned conjugate method improves both robustness and scalability [76] , and this is the solution route adopted here.

Preliminary tests of extrapolated liquid velocity approach
Four preliminary tests are carried out in this subsection to demonstrate the superiority of the extrapolated liquid velocity approach, with water used as the liquid and atmospheric air as the gas. The air density ρ G and dynamic viscosity μ G were 1.272 kg / m 3 and 1 . 86 × 10 −5 Pa · s ; the water density ρ L and dynamic viscosity μ L were 1002 kg / m 3  An illustration of the improved performance when the liquid velocity extrapolation technique was used is provided in Fig. 2 , for an initially static water droplet in a uniform airflow. The drop diameter was 3.1 mm with drop centre located initially at (8,0,0) mm, and the gas velocity was 15.7 m/s. To exclude the pressure jump across the interface, in and only in this test the surface tension was set to zero .
The snapshot was taken 200 time steps ( t = 2 × 10 −7 second) after the simulation was initialised, when the drop is still nearly spherical. Fig. 2 a shows the predicted pressure on a plane through the drop mid-section with conventional centred 2nd order accurate spatial discretisation schemes for the momentum equations [26] ; the pressure field is discontinuous across the interface and contains isolated 'spikes', caused by numerical error at interface cells. It is observed in the simulations that the velocity gradient in the gas phase is much larger than that in the liquid phase due to the high density and viscosity ratio. Therefore, considerable momentum error can be introduced if the neighboring gas velocity is used for the convection term discretisation in the liquid momentum CVs. In complete contrast, a smooth pressure distribution ( Fig. 2 b) was correctly predicted across the interface when using the extrapolated liquid velocity approach.
The momentum error observed above with conventional discretisation can lead to numerical breakup as shown in Fig. 3 (left).
The surface tension coefficient σ was set to 0.072 N / m in this and the following two preliminary tests. This figure thus illustrates the case of an initially static water droplet exposed to a uniform air flow at a Weber number of 3.4; experiments show for this condition that no breakup but just oscillatory deformation of the drop should be observed (the detailed morphology of droplet behaviour at various We will be presented below in Section 4 ). The correct oscillatory deformation mode can be numerically reproduced only after inclusion of a liquid velocity extrapolation technique as shown in Fig. 3 (right).
The importance of the extrapolated liquid velocity satisfying the continuity condition is demonstrated by simulating a moving water drop. The velocities of both the water drop and air flow are 9 m/s, resulting in a Weber number of 0, and thus the drop should therefore remain spherical. The drop diameter was 3.1 mm with drop centre located initially at (4,0,0) mm. Fig. 4 shows the morphology of the simulated drop after 1.2 ms. Without the divergence-free step for the extrapolated liquid velocity, the drop undergoes unphysical breakup due to numerical errors. With the   divergence-free step in the simulation, the drop almost retains its spherical shape with its center correctly moving to (14.8,0,0) mm, confirming that the implemented divergence-free algorithm can achieve satisfaction of the continuity condition for the extrapolated liquid velocity.
Finally, whilst it was momentum equation errors which motivated the introduction and design of the liquid velocity extrapolation approach adopted, it was observed that for consistency the ex-trapolated liquid velocity should also be used in the F and φ equations which determine the interface convection. Evidence to support this is given in Fig. 5 showing the predicted interface topology for the same case as in Fig. 3 of an initially static drop at the low Weber number of W e = 3 . 4 , where only oscillatory deformation should be seen. Fig. 5 indicates that when the momentumequation-deduced velocity U i is used for F and φ advection a considerably distorted interface is obtained, whereas with the  reconstructed liquid velocity U L i for interface advection a smooth interface was predicted as in the experiments. The underlying reason for this difference is examined in the following. Fig. 6 shows the momentum-equation-deduced velocity and constructed liquid velocity fields (after extrapolation). Both velocity fields satisfy the continuity equation while only the momentum equation deduced velocity field is physical as it also satisfies the momentum equation. However, the physical momentum-equation-deduced velocity represents the gas velocity rather than liquid velocity in cell ( i, j ) in the numerical calculation. Since the shear stress is the same across the interface, the velocity gradient in the liquid is much smaller than in the gas phase because of the large liquid/gas viscosity ratio. Therefore, the velocity of the liquid phase in cell ( i, j ) is much closer to the velocity in the neighbouring liquid cell (i, j − 1) . If the momentum-equation-deduced velocity is used for VOF function advection in cell ( i, j ), the liquid in cell ( i, j ) moves at the physical gas velocity, resulting in a large error. If the extrapolated liquid velocity u L is used, much better accuracy can be obtained in the calculation of liquid volume flux.
Based on the above evidence, the liquid velocity extrapolation procedure was adopted for all two-phase flow calculations reported below.

Plateau-Rayleigh instability
Plateau-Rayleigh instability is a surface tension influenced instability occurring when a stationary liquid cylinder with initial ra-dius R 0 is perturbed by a wave-like axially varying disturbance of initial amplitude η 0 ( η 0 sin ( kz ) where z is the cylinder axial direction and k is the perturbation wavenumber, related to wavelength λ by k = 2 π /λ). The perturbation will grow and the cylinder deform and eventually break up due to the Plateau-Rayleigh instability if kR 0 < 1. The perturbation grows according to η( t )sin ( kz ), with η(t) = η 0 e ωt ( ω is the growth rate). This problem was calculated using the fluid properties of air and water, an initial cylinder radius of R 0 = 0 . 14 m and a perturbation wavelength λ = 9 R 0 ( kR 0 = 0 . 698 ) and η 0 = R 0 / 28 = 0 . 5 , where is the uniform mesh spacing. The simulations were run with a domain size of 5 R 0 × 5 R 0 × 9 R 0 (uniform Cartesian mesh of 70 × 70 × 126). Periodic boundary conditions were used in the z -direction, and zerogradient conditions in the x and y directions. The predicted deformation and breakup of the liquid cylinder for this wavenumber is illustrated in Fig. 7 (i), showing an initial linear phase followed by a non-linear phase leading to breakup and ligament/droplet formation. The numerically predicted growth of perturbation magnitude in the initial phase could be accurately fitted with an exponential. Four further wavenumbers were calculated and Fig. 7 (ii) shows that the predicted growth rate agrees very well with the linear theory dispersion equation [77] .

Laminar liquid jet breakup in stagnant air
Transition from a dripping to a jetting mode for a laminar liquid jet has been studied experimentally by Clanet and Lasheras [67] . Water was injected downward into stagnant air at a velocity V J under gravity g, through a round tube. As the injected liquid velocity was increased, periodic dripping (PD), chaotic dripping (CD), and jetting (J) modes were observed: PD at very low velocity, liquid drops detach from the tube at a constant frequency, resulting in drops with constant mass. The detachment point is ∼1 diameter downstream of the tube exit. CD as liquid velocity increases over a first threshold, drops with different masses detach in a chaotic manner, instability waves are observed from the tube exit, and the detachment point moves downstream to a few diameters from tube exit. J as liquid velocity increases further, the detachment point moves suddenly to a downstream distance of greater than 10 diameters, and a smooth jet is formed upstream of the break point.
In the current present CLSVOF calculation of this problem, the tube diameter was D = 1 .  Fig. 8 shows that periodic dripping, chaotic dripping and jetting modes were indeed predicted respectively at these three velocities. These results indicate that the transition from dripping to jetting was correctly predicted by the present two-phase LES formulation.
The liquid jet breakup length was computed for the jetting mode in two further simulations for liquid velocities equal to 0.817 m/s and 0.999 m/s. Since the three cases at the higher exit velocities are all in the jetting regime, the breakup length should be proportional to the liquid velocity (see [78] ), and Fig. 9 indicates this was correctly predicted.

Resolution study
A water droplet in a uniform air flow is simulated for the temporal and spatial resolution study. The diameter of the initially spherical drop D 0 was 3.1 mm; the air density ρ G and dynamic viscosity μ G were 1.272 kg / m 3 and 1 . 86 × 10 −5 Pa · s ; the liquid density ρ L and dynamic viscosity μ L were 1002 kg / m 3    to reduce the computational cost. In study of droplet deformation and breakup, a characteristic time scale is commonly defined as t * = √ ρ L / ρ G D 0 / U G following [37,42] ; the dimensionless time is then defined as T = t /t * . First, temporal resolution was investigated on mesh M 1 . Three simulations were run with a time step of 1 × 10 −7 s, 2 × 10 −7 s, and 4 × 10 −7 s, corresponding to a CFL number of ∼ 0.025, ∼ 0.05, and ∼ 0.1 respectively. Fig. 10 shows the predicted shape of the deformed drop at T = 0 . 43 ( t = 2 . 4 ms ). The predicted drop shapes in the three simulations effectively collapse onto each other. The predicted cross-stream diameter at this moment is measured in the zoomed-in view as shown in Fig. 10 , resulting in 3.5156 mm, 3.5132 mm, and 3.508 mm. The cross-stream diameter difference between simulations with time steps of 2 × 10 −7 s and 4 × 10 −7 s is nearly twice that between simulations with time steps of 1 × 10 −7 s  and 2 × 10 −7 s, indicating that simulations with the current temporal discretisation scheme converge in first order as the time step decreases. This is reasonable since the first order forward Euler projection method was used. It is also observed that the crossstream diameter difference between simulations with time steps of 1 × 10 −7 s and 4 × 10 −7 s is very small ( ∼ 0.2%). Therefore, the time step corresponding to a CLF number of 0.1 was used in all following simulations.
In order to study the spatial resolution, simulations were carried out on meshes M 1 , M 2 , and M 3 respectively. Fig. 11 shows that the calculated drops deformed into the maximal liquid disc shape at T = 1 . 36 on meshes M 1 , M 2 , and M 3 . This implies that the initiation time (period of the deformation stage) is predicted well on all three grids, although the cross-stream dimension of the liquid disc shows some difference. Fig. 12 shows the temporal growth of drop cross-stream dimension ( D c ) predicted by LES on meshes M 1 , M 2 , and M 3 . The difference of calculated D c between meshes M 1 and M 2 is around twice that between meshes M 2 and M 3 , implying that D c converges in second order on mesh refinement.
Since the Rayleigh-Taylor instability (as the liquid drop is accelerated by the lighter gas) determines the drop breakup mode as detailed below, the drag which is correlated to the drop acceleration must be correctly reproduced in LES. Therefore, a convergence study is carried out for the drag coefficient C D which is related to the drag and the drop acceleration by: where F D is the drag, m is the drop mass, A and a is the cross section area and acceleration of the drop at the simulation initiation when the drop is still spherical. We first run the code by setting the velocity in the liquid to be 0 as if the drop is frozen until the gas flow around the drop fully developed. Then the two-phase simulation started running normally, and the drop acceleration at this moment is extracted from LES and used to calculate C D . The experimental value for the drag of a solid sphere at the corresponding [79] .
Four simulations were respectively run on grid M 0 , M 1 , M 3 , and M 4 . Fig. 13 shows the deviation of the predicted C D from C D, exp for increasing grid resolution. It is shown that C D converges by around first order (from M 0 to M 1 ) at low grid grid resolution and by nearly second order (from M 3 to M 4 ) at high grid resolution.

Effect of Weber number
The effect of the Weber number on the breakup morphology is investigated here for a single water droplet in a uniform air flow, with the underlying physical mechanism examined. 3D illustrations of simulated oscillatory deformation ( We = 3.4), bag breakup ( We = 13.5) and sheet-thinning breakup ( We = 100) have been presented in [69] . Four further simulations at Weber numbers of 12. The simulated droplet undergoes bag breakup at We = 12.5 as for the case with We = 13.5. A bag-stamen breakup mode is predicted for the droplet at We = 22 and 25, with a 3D view of the breakup process illustrated in Fig. 14 for the case We = 22. As in bag breakup, the drop first deforms into a disc, with a noticeably larger dimension due to increased aerodynamic forces, see Figs. 14 a-c. Fig. 14 d indicates that the liquid film center is thick and a bag now grows between the rim and center. Eventually the bag bursts, leaving a liquid rim and a central liquid stamen ( Fig. 14 e-f).   [80] . Then the drop further deforms into a liquid disk ( Fig. 15 e), and a thin liquid sheet forms at the liquid disk periphery ( Fig. 15 f). Under the action of aerodynamic forces, thin liquid sheet is blown downstream ( Fig. 15 g), and the thin liquid film disintegrates ( Fig. 15 h), forming a large number of ligaments aligning to the air-stream direction ( Fig. 15 i). As the liquid ligaments disintegrates from the liquid disk, the centre part of the liquid disk evolves into a liquid core ( Fig. 15 j). The whole simulated multimode breakup process agrees well with the experimental observations for the case of We = 49 and We = 53 from Zhao et al. [41] [80] .
According to experiments reported in [51] and [41] , oscillatory deformation happens for 2.5 < We < 12, bag breakup when 12 < We < 16, bag-stamen breakup when 16 < We < 28, multimode breakup when 41 < We < 80, and sheet-thinning breakup when We > 80. For the Weber numbers chosen for the above test cases, the deformation/breakup morphology of the simulated droplet agrees well with experimental observations.
In order to explore further the mechanism of single drop deformation and breakup, a detailed analysis of velocity and pressure fields has been carried out. Fig. 16 shows predicted velocity vectors and pressure contours for We = 13.5 at the early time of T = 0.036, when the drop is still nearly spherical with only small induced liquid velocities. Similar to the flow around a solid sphere, the gas velocity reduces to zero at the front stagnation point, resulting in the highest gas phase pressure ( Fig. 16 a). The gaseous flow around the drop periphery accelerates, indicated by the low pressure zones. Vortices develop in the low pressure wake of the drop. The liquid velocity vectors and liquid in-plane streamtraces are visualised in Fig. 16 b. Since the droplet is still nearly spherical, the liquid pressure on the interface is the sum of gas pressure and the almost constant pressure due to the jump arising from surface tension. Thus, the spatial pressure distribution within the liquid phase is directly determined by the gas pressure field. The pressure gradients inside the drop induce the liquid velocity, which moves liquid radially outward from front and rear stagnation regions to the drop periphery as indicated by the streamtraces. Fig. 17 shows a similar picture as in Fig. 16 b but now for the highest Weber number of 96, and at a later time of T = 0.29 so that considerable drop shape distortion has occurred. On close examination, no evidence of any boundary layer can be seen in the liquid phase. This contradicts arguments previously put forward for the physical mechanism behind the shear-stripping mechanism,    which postulated that a liquid boundary layer developed adjacent to the interface under the action of shear from the gas flow. The predicted liquid velocity vectors are effectively mainly controlled by the gaseous pressure distribution around the drop, with shear forces playing an insignificant role even at this high Weber number case. Fig. 18 shows the velocity and pressure fields for the case We = 13.5 at the later time of T = 1.3 when the drop has deformed into a disc. Due to the large aerodynamically induced front/back pressure difference, the high-density liquid disc is accelerated considerably by the low-density air flow, resulting in a Rayleigh-Taylor (RT) instability [81,82] . At this Weber number, the RT wavelength is larger than the maximum cross stream diameter, producing a bag in the disc centre ( Figs. 19 d-e). As We increases to 22, the RT wavelength becomes shorter than the maximum cross stream diameter, resulting in a liquid disc with a greater thickness ventrally; the bag then forms between rim and thicker liquid centre  ( Figs. 20 e-f). When We grows to 50, the periphery of the formed liquid disk is thinner than the disk center as shown in Fig. 21 e, which is significantly different from the bag breakup mode (disk center is thinner than periphery as shown in fig. 19 c) and bagstamen breakup (disk center is nearly as thick as periphery as shown in Fig. 20 c). The pressure contour lines in Fig. 22 show that the pressure gradient in the liquid disk periphery is higher than in the disk center. Therefore, the disk periphery is subject to stronger acceleration than the disk center, resulting in that the disk periphery bends downstream and the RT wave develops at the disk periphery as shown in Figs. 21 f-g. As We increases further to 96, the RT wavelength decreases even more, and two waves of liquid films form and disintegrate sequentially from the main drop before leaving a central liquid core behind ( Fig. 23 ). The rim of the RT wave becomes thinner and holds less mass due to its decreased wavelength as We increases; thus, when We = 96, liquid wave rims are more prone to downstream deflection due to the smaller inertia of the rim liquid.
The droplet breakup process is divided into a deformation stage and a breakup stage. In the first stage, deformation into a liquid disc occurs under the pressure imbalance created by the gas flow. In the second, because the liquid disc is accelerated by the lighter gas phase, an RT instability is induced, forming a thin bag (bag or bag-stamen breakup) or a wave/ridge of liquid sheet (sheetthinning breakup) on the drop periphery; it is these liquid structures which subsequently disintegrate into ligaments and droplets. The elapsed time of the deformation period (which ends when the    disc reaches its maximum dimension ) is defined as the breakup initiation time T ini , and the disc diameter at this instant is defined as the maximum cross-stream dimension D max as in [41] . Figs. 19 ac demonstrate that the liquid drop deforms into a liquid disc at the first three time snapshots whilst a bag can just be observed forming at the third. At T = 1.36, the liquid disc reached its maximum cross-stream dimension D max /D 0 = 1 . 73 , and this is used to define the simulated initiation time. T ini for a bag-stamen breakup mode is defined as the time instant shown in Fig. 20 c; for multimode breakup, the initiation time is the moment when the drop evolves into a liquid disk as shown in Fig. 21 e; for sheet-thinning breakup, the end of the deformation period is indicated in Fig. 23 d.
Predicted initiation time T ini and maximum cross-stream dimension D max at different Weber numbers are compared with experiments in Figs. 24 and 25 . Fig. 24 indicates that the deformation period calculated from the simulations reduces as We increases, which is consistent with the tendency observed by [51] in their review of available experiments. The predicted T ini lies between Pilch and Erdman's data and Hsiang and Faeth's data [54] for We larger than ∼ 15. As the drop accelerates to the freestream velocity, the aerodynamic force exerted on the drop is decreasing. The drop will either experience breakup in a finite time or only undergo deformation. Pilch and Erdman's data for Weber numbers approaching the critical value We cr are doubtful since they imply infinite initiation time. [55] studied bag breakup at Weber numbers close to critical value in detail, and reported a non-dimensional initiation time of 1.42 at Oh = 1 . 9 × 10 −3 . T ini predicted by the current methodology agrees well with this experimental result when We approaches the critical Weber number when breakup first begins We cr . The predicted maximum cross-stream dimension at different Weber numbers obtained is presented in Fig. 25 . Zhao et al.s experiments [41] found that D max grows as We increases in the bag breakup regime, and this is reproduced well by the current algorithm. For bag-stamen and sheet-thinning breakup, the predicted D max is approximately constant, located between the experiments of [41] and [83] .
In order to calculate the RT instability wavelength, the acceleration of the drop is first computed as follows: • The velocity of the centre of mass is computed from: The acceleration of the drop over time is shown in Fig. 26 , and the drag coefficient derived from the drop acceleration agrees well with the experimental measurements as shown in [26,69] . The acceleration is 230 m / s 2 at the initiation time T ini = 1 . 36 when the drop reaches its maximum cross stream dimension D max = 5 . 363 mm . The wavelength of the most unstable RT wave is calculated from λ max = 2 π 3 σ ρ L a L = 6 . 08 mm Following the definition in [41] , the nondimensional RT wave number in the maximum cross stream dimension is computed as:  [41] that bag breakup occurs when 1 / √ 3 < N RT < 1 and bag-stamen breakup when 1 < N RT < 2. This also confirms the descriptions provided above that it is a Rayleigh-Taylor instability that dominates the drop deformation and breakup process for the Weber number range studied.

Effect of Ohnesorge number ( Oh
The present numerical methodology was used to investigate the effect of Oh on We cr . Simulations were run at different Ohnesorge numbers by adjusting liquid viscosity. When water viscosity ( 0 . 892 × 10 −3 Pa · s ) was used as liquid viscosity as in the above subsections, Oh = 1 . 9 × 10 −3 . Simulations were run here for three Ohnesorge numbers 0.1, 0.7, and 2, by setting the liquid dynamic viscosity to 0.0473 Pa · s , 0.331 Pa · s , and 0.946 Pa · s respectively. Since the drop breakup time grows as liquid viscosity increases, the breakup occurs further downstream. Thus, a simulation domain with a larger streamwise dimension than used above was needed. The domain size was altered to [0 115 . 2] × [ −9 9] × [ −9 9] mm in the x, y , and z directions respectively. In order to reduce computational cost, a uniform fine mesh was used in the region [0 115 . 2] × [ −4 . 8 4 . 8] × [ −4 . 8 4 . 8] mm with cell size 0.12 mm, and an expanding mesh was used elsewhere. Since the liquid flow inside the drop remains laminar in the drop deformation and breakup process, the computed SGS eddy viscosity must be small enough if it is not to interfere with the physical processes. Fig. 27 illustrates the computed Smagorinsky SGS eddy viscosity nondimensionalised by local gas/liquid molecular viscosity for the case with Oh = 1 . 9 × 10 −3 . In the gas phase, the SGS eddy viscosity took effect mainly in the wake region behind the drop and in the region adjacent to the interface to dissipate the smallest resolved eddy energy, ensuring the robustness of the proposed methodology. In the liquid phase, the SGS eddy viscosity was one order of magnitude smaller than the liquid molecular viscosity. At the low Oh of 1 . 9 × 10 −3 , the liquid molecular viscosity will not influence the drop deformation/breakup, and this is also the case for the smaller SGS eddy viscosity in the liquid. As the liquid molecular viscosity grows to satisfy Oh = 0 . 1 , the velocity gradient inside the drop decreases as shown in Fig. 28 , and thus the computed SGS eddy viscosity declines in the liquid phase. Fig. 28 shows that the computed SGS eddy viscosity in the liquid was three orders of magni-  tude smaller than the high liquid molecular viscosity for the case with Oh = 0 . 1 . For the two cases with Oh = 0 . 7 and Oh = 2 , the non-dimensionalised SGS eddy viscosity in the liquid phase will be even smaller. Thus the effects of the liquid molecular viscosity on the drop deformation and breakup were correctly captured in the simulations without the pollution of SGS eddy viscosity. Fig. 29 shows that two simulation cases were run for each Ohnesorge number: the first always displayed oscillatory deformation, and the second, at a slightly higher Weber number, displayed bag breakup. Based on these results, a line fitted between the simulations ( W e cr = 12 . 3 ( 1 + 1 . 1 Oh ) ) agreed very well with the measurements of Lane [84] , Hinze [85] , Hanson et al. [86] , Loparev [87] , and Hsiang and Faeth [38] , with the fitted correlation line matching Cohen's energy transfer analysis [53] . Finally, Fig. 30 shows the initiation time predicted in the four current LES  breakup cases of Fig. 29 . Pilch and Erdman [51] , Hsiang and Faeth [54] and Gelfand et al. [55] have all proposed empirical correlations as shown in Fig. 30 ; the present results support Gelfand's correlation very closely since the data used by Gelfand were measured for bag breakup at approximately critical Weber numbers as is the case in the four simulations. These results confirm the accuracy of the present methodology for predicting We cr and the effects of Oh .

Summary and conclusions
A quasi-DNS/LES algorithm for two-phase flows based on a CLSVOF formulation has been developed and applied in the present paper. Most importantly, the two-phase flow governing equations are discretised using an extrapolated liquid velocity approach, with sharp treatment of the density and molecular and SGS eddy viscosity. Several test cases were presented to demonstrate that the proposed method is accurate and robust, even at high liquid/gas density ratio. For the low-speed liquid jet, the transition from dripping to jetting mode is correctly predicted by the current twophase flow solver. Simulations of breakup of a single drop showed that the breakup morphology can be correctly reproduced at corresponding We by the developed two-phase solver. The calculated velocity and pressure fields indicates that it is the Rayleigh-Taylor instability that determines the breakup mode in the studied We range, with the calculated Rayleigh-Taylor wavelength in good agreement with the experimental measurements. Finally, the effect of liquid viscosity on drop deformation and breakup is properly captured by the current simulations, showing that the critical Weber number required for breakup mode increases linearly with Ohnesorge number at high Oh . The initiation time characterising the start of droplet breakup is accurately predicted at different We and Oh numbers.
The present algorithm is developed for two-phase flow with high viscosity and density ratio where the velocity gradient in the liquid is much smaller than in the gas. For the two-phase flow with low viscosity and density ratio ( O (1)) where the velocity gradient across the interface is in the same order of magnitude, considerable error will be introduced as the interface velocity is approximated by the neighbouring liquid velocity in the current method. Further work is required to develop a proper approximation to the interface velocity for all cases.