domain: nlin spi

text:
Cristian R Constante-Amores

Lyes Kahouadji

Assen Batchvarov

Seungwon Shin

Jalel Chergui

Damir Juric

Omar K Matar

Direct numerical simulations of transient turbulent jets: vortex-interface interactions

The breakup of an interface into a cascade of droplets and their subsequent coalescence is a generic problem of central importance to a large number of industrial settings such as mixing, separations, and combustion. We study the breakup of a liquid jet introduced through a cylindrical nozzle into a stagnant viscous phase via a hybrid interface-tracking/level-set method to account for the surface tension forces in a threedimensional Cartesian domain. Numerical solutions are obtained for a range of Reynolds (Re) and Weber (We) numbers. We find that the interplay between the azimuthal and streamwise vorticity components leads to different interfacial features and flow regimes in Re-We space. We show that the streamwise vorticity plays a critical role in the development of the three-dimensional instabilities on the jet surface. In the inertiacontrolled regime at high Re and We, we expose the details of the spatio-temporal development of the vortical structures affecting the interfacial dynamics. A mushroomlike structure is formed at the leading edge of the jet inducing the generation of a liquid sheet in its interior that undergoes rupture to form droplets. These droplets rotate inside the mushroom structure due to their interaction with the prevailing vortical structures. Additionally, Kelvin-Helmholtz vortices that form near the injection point deform in the streamwise direction to form hairpin vortices, which, in turn, trigger the formation of interfacial lobes in the jet core. The thinning of the lobes induces the creation of holes which expand to form liquid threads that undergo capillary breakup to form droplets.

Introduction

The breakup of a dispersed fluid in a stagnant phase is a classical problem in multiphase flows, as it encompasses a multitude of interfacial singularities, which exemplify situations wherein the interface undergoes topological transitions, e.g., liquid-threads breakup into drops, and merging of drops/bubbles due to coalescence. These transitions involve the development of singularities where interfacial distances vanish and velocity fields diverge, and the system is controlled by a combination of capillary, inertial, and viscous forces. The complex topological nature of the jet phenomenon has fascinated the scientific community for decades which has led to numerous, comprehensive reviews, see for example [START_REF] Lin | Drop and spray formation from a liquid jet[END_REF]; [START_REF] Eggers | Physics of liquid jets[END_REF]; [START_REF] Lasheras | Liquid jet instability and atomization in a coaxial gas stream[END_REF].

Jet breakup is influenced by a range of multi-scale physics; these include the interaction of turbulence with interfaces (creating cascades of motion featuring a large separation of scales), capillarity, potentially complex rheology, the presence of fields (e.g. gravitational, electromagnetic), as well as heat transfer and phase change. [START_REF] Reitz | Mechanisms of breakup of round liquid jets[END_REF] proposed four jet breakup regimes depending on the appearance of the jet far downstream from the injection point. In the 'Rayleigh regime', the onset of breakup is the Rayleigh-Plateau instability, in which the growth of the linear modes on the jet-surface leads to the formation of large droplets with respect to the jet-nozzle. In the 'first/second windinduced' regime the resulting droplets have roughly the same scale or smaller as the jet-nozzle. Finally, in the 'atomisation regime', the generated turbulence helps to create spatio-temporal chaos resulting in droplets which are up to four orders of magnitude smaller than the size of the injection nozzle. The physics of droplet generation in turbulent jets remains partially-understood despite the significant scientific attention it has received over the years [START_REF] Dombrowski | A photographic investigation into the disintegration of liquid sheets[END_REF][START_REF] Hoyt | Waves on water jets[END_REF][START_REF] Lasheras | Liquid jet instability and atomization in a coaxial gas stream[END_REF]Marmottant & Villermaux 2004;[START_REF] Villermaux | Ligament-mediated spray formation[END_REF][START_REF] Eggers | Physics of liquid jets[END_REF]. The rest of this introduction aims to provide an up-to-date summary of the experimental and numerical efforts to study the jet breakup phenomenon.

The ground-breaking experiments of [START_REF] Hoyt | Waves on water jets[END_REF] and [START_REF] Taylor | Water jet photography techniques and methods[END_REF] showed the complex topological features of the surface waves formed on the jet. They revealed that these waves are responsible for the transition from laminar to turbulent flow, and found underlying similarities of the occurring instabilities to inviscid linear theory. Marmottant & Villermaux (2004) with their pioneering experiments scrutinised the various stages of the jet dynamics: from the growth of linear modes (through a Kelvin-Helmholtz instability) that characterises the early-time dynamics, to the development of nonlinearities leading to 'primary' and subsequent 'secondary' breakup events (through long filament pinchoff modulated by a Rayleigh-Taylor instability), and the formation of a cascade of droplet sizes. These authors found that mean droplet-size is proportional to the wavelength selected during the RayleighTaylor instability (also observed by [START_REF] Varga | Initial breakup of a small-diameter liquid jet by a high-speed gas stream[END_REF]). In the same premise, [START_REF] Kooij | What determines the drop size in sprays?[END_REF] showed that the droplet size distribution is also a function of the nozzle geometry and the surrounding pressure. More recently, [START_REF] Ibarra | On the near-field interfaces of homogeneous and immiscible round turbulent jets[END_REF] presented the results of an experimental study of the spatial evolution of turbulent immiscible liquid-liquid jets; however, a detailed account of the spatio-temporal development, and critical mechanisms leading to droplet generation remains outstanding.

The multi-scale nature of the flow, and the complex interfacial topology complicate the experimental scrutiny of the different physical mechanisms occurring across the scales. Thus, elucidating the fundamental physics of this problem has also relied on highfidelity simulations exemplified by the work of [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF]; [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF]; [START_REF] Gorokhovski | Modeling primary atomization[END_REF]; [START_REF] Desjardins | Detailed numerical investigation of turbulent atomization of liquid jets[END_REF]; [START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF]; [START_REF] Herrmann | On simulating primary atomization using the refined level set grid method[END_REF]; [START_REF] Chenadec | A stable and conservative framework for detailed numerical simulation of primary atomisation[END_REF]; [START_REF] Desjardins | Direct numerical and largeeddy simulation of primary atomization in complex geometries[END_REF]; [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF]; [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF]; [START_REF] Agbaglah | Numerical simulation of the initial destabilization of an air-blasted liquid layer[END_REF]; [START_REF] Zandian | Understanding liquid-jet atomization cascades via vortex dynamics[END_REF]; [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF]; [START_REF] Zandian | Vorticity dynamics in a spatially developing liquid jet inside a co-flowing gas[END_REF]. The first step towards the use of numerical simulations to understand the physical mechanisms at play was conducted by [START_REF] Desjardins | Detailed numerical investigation of turbulent atomization of liquid jets[END_REF], who were able to identify a sequence of essential steps during the spatio/temporal interfacial development of a planar liquid jet-segment: the formation of initial corrugations on the surface, followed by the development of ligaments whose capillary instability leads to droplet formation. They also showed that the early interfacial corrugations are a consequence of the turbulent eddies which carry enough kinetic energy to overcome capillary forces.

Using a similar approach to Desjardins & Pitsch ( 2010), but for a cylindrical liquid jet-segment surrounded by an outer gas phase, [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF] provided a comprehensive study of the flow structures in terms of the vortex-surface interaction. They showed that 'hole formation' of a liquid-sheet is an essential requirement to trigger the formation of droplets, and the thinning of the liquid-sheet is driven by the superposition of hairpin-vortices near the interface rather than by capillarity action. Similar findings have been reported for a planar liquid jet-segment surrounded by an outer gas phase [START_REF] Zandian | Understanding liquid-jet atomization cascades via vortex dynamics[END_REF], and for the transient dynamics of a cylindrical liquid jet surrounded by a coaxial air phase [START_REF] Zandian | Vorticity dynamics in a spatially developing liquid jet inside a co-flowing gas[END_REF]. [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF] performed simulations of a two-phase mixing layer between parallel gas and liquid streams to investigate the interfacial dynamics, and the statistics for the multiphase turbulence.

They also observed that the formation of ligaments, and subsequently-formed droplets, are triggered by the hole-induced perforation of the liquid sheets.

In the previous numerical studies, interface-capturing capabilities were used to account for the surface tension forces in the absence of intermolecular forces (i.e., disjoining pressure). For static sheets, the disjoining pressure is neglected as the minimum computational cell is larger than the film sheet thickness in which the intermolecular forces will drive its perforation. Therefore, the hole formation is an outcome of the numerical cut-off interfacial length scale, i.e., minimum mesh size (O(10 -6 )m). Nevertheless, recent experiments [START_REF] Kooij | What determines the drop size in sprays?[END_REF][START_REF] Marston | Crown sealing and buckling instability during water entry of spheres[END_REF][START_REF] Néel | The spontaneous puncture of thick liquid films[END_REF]) demonstrate the existence of hole formation in dynamic sheets with a characteristic film thickness on the order of microns.

In this study, we aim to provide a comprehensive explanation of the physical mechanisms governing the interfacial dynamics of turbulent jets focusing on the less well-studied liquid-liquid systems. We will perform high-resolution three-dimensional direct numerical simulations using a hybrid interface-tracking/level-set approach to resolve the interfacial dynamics. We will demonstrate how the interaction between the vortical structures, which accompany the development of the flow, and the interface influence the mechanisms underlying droplet generation over a wide range of Reynolds and Weber numbers.

The rest of this paper is organised as follows. In Section 2, we present the governing equations along with the numerical technique used to carry out the simulations. In Section 3, we begin with the presentation of a regime map in Re -W e space, classifying jet spatio-temporal development, followed by an in-depth discussion of the vortex-surface interactions linked to the topological changes; in addition, we elucidate the role of hole formation as a precursor to droplet generation. Finally, concluding remarks are provided in Section 4.

Problem formulation, and numerical techniques

Figure 1 shows a three-dimensional representation of the interface highlighting several specific features that are discussed in detail in the present work: Kelvin Helmholtz (KH) waves close to the injection nozzle, outer lobes formed on the main body of the jet, a leading edge 'mushroom'-like structure, and droplets resulting from the atomisation process. The numerical framework is based on solving the two-phase incompressible Navier-Stokes equations in a three-dimensional Cartesian domain x = (x, y, z). The surface tension force in the momentum equation is treated by using a hybrid fronttracking/level-set method presented by [START_REF] Shin | Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity[END_REF]. The governing flow equations

∇ • u = 0 , ρ ∂u ∂t + u • ∇u = -∇p + ∇ • µ(∇u + ∇u T ) + F s , (2.1)
where t, u, p, and F s stand for time, velocity, pressure, and the surface tension force, respectively. The density and viscosity are expressed using the following formulation ρ = ρ so + (ρ wρ so ) H (x, t), µ = µ so + (µ wµ so ) H (x, t).

(2.2) wherein H (x, t) represents a smoothed Heaviside function, which is zero in the dispersed phase (water) and unity in the stagnant phase (silicone-oil), while the subscripts 'w' and 'so' refer to the individual phases. The surface tension force F s is defined by using the hybrid formulation as in [START_REF] Shin | A hybrid interface method for three-dimensional multiphase flows based on front-tracking and level set techniques[END_REF] and [START_REF] Shin | A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows[END_REF])

F s = σκ H ∇H, (2.3)
where σ refers to surface tension, and κ H is twice the mean interface curvature calculated from the the Eulerian grid by using

κ H = F L • G σG • G, (2.4)
where

F L = Γ (t) σκnδ(x -x f )ds, G = Γ (t) nδ(x -x f )ds. (2.5)
Here, x f is the parametrisation of the interface, Γ (t), and δ(xx f ) is a Dirac delta function which vanishes everywhere except at the interface; n is the outward-pointing unit normal vector to the interface, and ds is the length of the surface element.

The incompressible Navier-Stokes equations (2.1) are solved by using a second-order finite-difference method on a staggered grid [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de navier-stokes[END_REF]). The computational domain is then discretised by a fixed, regular, Eulerian grid, and the spatial derivatives are approximated by standard centred difference discretisation, except for the non-linear term, which makes use of a second-order essentially non-oscillatory (ENO) scheme [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. As for the viscous term, we use a second-order centred difference scheme. We have used the projection method to handle the incompressibility condition combined with a multigrid iterative method for solving the elliptic pressure Poisson equation [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Kwak | Multigrid algorithm for the cell-centered finite difference method ii: Discontinuous coefficient case[END_REF]. The numerical method uses an additional adaptive Lagrangian grid based on a hybrid front-tracking/levelset method [START_REF] Shin | High order level contour reconstruction method[END_REF][START_REF] Shin | A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows[END_REF] to track the interface location. The geometrical information, the pressure and velocity variables are exchanged between the adaptive Lagrangian mesh and the fixed Eulerian grid following the Immersed Boundary Method of [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF]. The physical elements that form the Lagrangian interface are advected according to dx f /dt = V, where V is the interface velocity interpolated numerically from the fixed Eulerian grid, using a second-order Runge-Kutta method. Finally, the numerical method uses a domain-decomposition technique for its parallelisation and Message Passing Interface (MPI) for exchanging information between adjacent subdomains. More information regarding the full implementation of the numerical method can be found in [START_REF] Shin | A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows[END_REF][START_REF] Shin | A hybrid interface tracking -level set technique for multiphase flow with soluble surfactant[END_REF] 2.1. Numerical configuration and physical parameters

Figure 1 shows the three-dimensional computational domain, which is a rectangular box of size 20D × 4D × 4D, where D stands for the inner diameter of the nozzle (e.g., 4mm). The jet is produced when water leaves the cylindrical nozzle to enter progressively into the stagnant silicone oil. The physical properties of the fluids are given in table 1. An inflow boundary condition is applied to the nozzle on the left of the domain, i.e., at (x = 0), which follows a simplified power-law turbulent velocity profile:

u (r, t) = 15 14 U 1 - r D/2 28 (1 + A sin (2πf t)) .
(2.6)

Here, U , A and f , stand for the average injection velocity, amplitude and frequency, respectively, of the external pulsatile perturbation. The radial distance, r, within the jet measured from its centreline (y 0 , z 0 ) is r = (yy 0 ) 2 + (zz 0 )

2 . The values for A and f are informed by the previous work of [START_REF] Ling | Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model[END_REF] (e.g., A = 0.05 m/s and f = 20 Hz). The numerical setup closely follows other computational studies; for example, we impose a free boundary condition on the walls of the computational domain to let the fluid freely enter or leave the boundaries [START_REF] Taub | A direct numerical simulation study of higher order statistics in a turbulent round jet[END_REF][START_REF] Ling | Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model[END_REF][START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF]. A pressure outflow boundary condition is applied on the right surface of the domain to allow the fluid to exit the domain. The solid nozzle is treated as a no-slip surface [START_REF] Asadi | On the scaling of propagation of periodically generated vortex rings[END_REF]. The distance x, velocity u, time t, and pressure p in equation (2.1) are rendered dimensionless using the following characteristic scales, D, U , D/U , and ρ w U 2 , respectively. Hence, the dimensionless control parameters governing the phenomena we will study are given by:

Re = ρ w U D µ w , W e = ρ w U 2 D σ , (2.7)
where Re stands for the Reynolds number (i.e. the ratio of inertial to viscous forces), and W e represents the Weber number (i.e. the ratio of inertial to capillary forces). All the variables appearing in the equations and boundary conditions are rendered dimensionless using the aforementioned scalings, unless stated otherwise.

The first part of the results section corresponds to a discussion of a regime map of jet dynamics in the Re -W e space. In this instance, we have used the M2-mesh (see table 4 in the Appendix) to perform the fully three-dimensional simulations. The selection of this mesh is dictated by the need to map out parameter space relatively rapidly before focusing on elucidating the details of the dynamics using the M3-mesh, which provides higher resolution, for four cases; the Re -W e combinations for these cases are listed in table 2. Information regarding the mesh-refinement study, resolution considerations, and the validation for this work are detailed in the Appendix.

Results

Interfacial dynamics: phase diagram in Re -W e space

We start the discussion of the results by presenting a phenomenological picture of the interfacial dynamics in a phase diagram in Re -W e space. The Reynolds and Weber numbers range between 10 3 -10 4 and 30 -900, respectively. Figure 2 shows the regime map in terms of the interfacial dynamics predicted from our numerical simulations. Several features emerge from this figure, which we have divided into four phenomenological regions based on the appearance of the interfacial structures. We aim to quantify the different jet behaviours by close inspection of the interplay between the vorticity field, ω = ∇×u, and the interface. Such in-depth analysis is carried out following the introduction of each region, utilising the higher resolution M3-mesh examples, as shown in table 2.

Region 'A' in figure 2 is characterised by low Reynolds and Weber numbers and defined by the dominance of the capillary over inertial forces, yielding an axisymmetric behaviour of the jet. The formation of neither a leading mushroom-like structure nor interfacial lobes is observed. Figure 3 shows the interplay between the azimuthal and streamwise vorticity components, ω θ and ω x , respectively for Case 1 (see table 2). We observe that ω θ exceeds ω x by two orders of magnitude for the entirety of the jet, which explains the lack of deformation of the jet-core, and its axisymmetric shape. Region 'B' is defined by low Reynolds and high Weber numbers. As shown in figure 2, the snapshots of the interface in this region of parameter space reveal the development of interfacial waves as well as the formation of a 'mushroom'-like structure at the jet leading edge. This is due to the large capillary pressure generated as a result of the interfacial curvature at this leading edge, which induces radial flow that is retarded by the viscous resistance from the stagnant phase. Figure 4 shows that although the magnitude of ω θ exceeds that of ω x by almost two orders of magnitude for Case 2, the relative significance of ω x has increased in comparison with Case 1. This is correlated with the corrugations which are observed on the main body of the jet in Case 2 that appear to be largely absent in Case 1. These corrugations are related to the development of Kelvin-Helmholtz (KH) instabilities that arise from the velocity difference between the injected jet and the surrounding, initially stagnant phase.As shown in figure 4, the mushroom structure is also accompanied by the formation of droplets within it. More information regarding the mechanisms which induce droplet formation is provided in Section 3.3.

Region 'C' is characterised by low Weber and high Reynolds numbers. To elucidate the mechanisms at play in this region, we refer to Case 3 (see table 2), as shown in figure 5. It is seen clearly that the magnitude of ω θ is an order of magnitude larger than that of ω x . The higher levels of inertia in Case 3, in comparison to Case 2, leads to formation of pronounced KH instabilities with vortical rings close to the interface, yielding spanwise interfacial deformations. Close inspection of figure 5 reveals that the mushroom structure, which is also present in Case 3, undergoes significant deformation leading to the formation of a toroidal sheet, which envelopes that main body of the jet.

In region 'D', the flow is characterised by both high Reynolds and Weber numbers. Here, it is seen from the snapshots of the interface in figure 2 that the interfacial dynamics in this region of Re -W e space are the most complex. The mushroom structures suffer severe deformation as does the main body of the jet which also features the formation of lobes arising from the KH-induced corrugations. It is also evident that flow in region D is also accompanied by droplet generation. In Section 3.2 below, we focus on Case 4 in figure 2 and provide an extensive explanation of the mechanisms linking droplet formation to vortex-interface interaction.

Interfacial dynamics explained through vortex-surface interaction

This section focuses on the inertia-controlled Region D and Case 4 as a principal example of the flow in this region. In figure 6, we examine the spatio-temporal interfacial dynamics overlaid with the magnitude of ω θ and ω x in the xz plane. During the early stages of the injection, the mushroom-shape structure develops due to the capillary-driven radial flow arising from the large interfacial curvature at the jet leading edge, which is resisted by the stagnant phase. As time evolves, the entrainment and formation of a toroidal liquid sheet of stagnant phase are observed inside the mushroom structure (see for example figure 6a,b). Upstream of the mushroom, it is seen that the KH instability develops (see figure 6c) amplified by the pulsatile-injection into waves which cause a local adverse pressure gradient by virtue of the local interfacial curvature (see figure 6d). With increasing time, we observe the formation of outer lobes as a result of the entrainment of the stagnant phase in the jet-core (see figure 6e).

Figure 6 also highlights the spatio-temporal development of the azimuthal and streamwise components of the vorticity. At early times, vorticity generation coincides with the velocity boundary layer attached to the interface, which corresponds to strong tangential flow near the interface. The roll-up of the shear layers, and subsequently, the vortex rollup gives rise to the formation of KH vortex rings close to the interface (see figure 6c). As time increases, the vorticity boundary layer is convected towards the jet-core (see figure 6e). Additionally, vorticity dissipates strongly in the stagnant phase due to the damping effect of the viscosity.

Attention is now turned towards the competition between ω θ and ω x . In the region adjacent to the injection point, the azimuthal vorticity component dominates over its streamwise counterpart by two orders of magnitude; this dominance is reflected by the fact that the KH vortex rings shown in figure 6c are essentially axisymmetric. As the flow develops downstream (see figure 6d,e), ω x becomes comparable in magnitude to ω θ leading to streamwise stretching of the KH vortex rings to form hairpin-shaped vortices. Hairpin vortical or 'horseshoe' structures were proposed by [START_REF] Theodorsen | Mechanism of turbulence[END_REF] to describe Further details about the specific alignment of the vortical structure is provided below. Figure 7a presents a three-dimensional representation of the jet at t = 28.97 alongside five spatial locations which will be used to show the evolution of vorticity and velocity across the jet core. The specific locations were chosen as follows: panel 'A' is representative of the dynamics near the injection point, where there is little disturbance to the interface; panel 'B' displays the cross-section of the jet in a jet-ring (i.e., the part of the interface where we observe the capillary wave crest), whereas panel 'C' depicts a jet-braid (i.e., the capillary wave trough); panels 'D' and 'E' portray the jet core dynamics where lobe formation is present. Figure 7b shows the velocity and vorticity line profiles at the y = 2 plane alongside the streamwise vorticity in the y-z plane for each of those locations. Near the injection point (location 'A'), ω θ >> ω x , which is in agreement with the lack of deformation of the vortical structures in the streamwise direction (see also figure 6f). The tangential motion of the fluid close to the surface results in the emergence of the velocity boundary layer. In turn, this causes the formation of two large peaks of opposite signs for ω θ . As the flow evolves downstream, ω θ still dominates the physics, although, its value has been reduced in favour of an increase in ω x , and consequently the local formation of corrugations on the jet-core (location 'B'). Additionally, the ω θ peaks have widened in their base due to the growth of the velocity boundary layer as the jet moves downstream [START_REF] Liepmann | The role of streamwise vorticity in the near-field entrainment of round jets[END_REF]. Further downstream (location 'D' and 'E'), the interface has lost its cylindrical shape due to the entrainment of the stagnant phase as ω x becomes comparable in magnitude to ω θ . Several peaks can be seen in the vorticity profiles as the probe line passes through several lobes. By inspecting the instantaneous velocity-fields in locations 'D' and 'E', we observe that the main component of velocity is associated with the direction of injection, where there is a decay of its mean centreline value as flow evolves downstream, which is in agreement with [START_REF] Pope | Turbulent Flows[END_REF] for single-phase jets. The velocity decay leads to a reduction in the velocity difference between the injected and stagnant phases, which attenuates the interfacial shearing, supporting the development of the vortical rings close to the surface of the jet. Notice that the y and z velocitycomponents are positive or negative close to the surface, showing the entrainment of injected fluid towards the stagnant viscous phase or vice versa.

Next our attention is turned towards the physical mechanisms which lead to the alignment of the hairpin-like vortical structures along the free-surface. Figure 7c shows cross-section panels with ω x of the jet at the streamwise locations shown in figure 7a. As depicted in panel 'A', ω x is mainly confined inside of the injected stream, where hairpin vortex legs are observed (see first panel of figure 7c). In panel 'B', the outer layer of vorticity comes from the braid located immediately downstream, which starts to roll-up around the ring. Further analysis of panel 'B' shows the existence of additional pairs of vortex legs which are 180 • out-of-phase with respect to the inner layer, and subsequently they correspond to a second hairpin vortex layer with opposite direction (e.g., upstream direction). Similarly, in panel 'C' we observe that the inner layer of the vorticity comes from the ring located upstream, whereas the outer layer comes from the ring located immediately downstream. Therefore, we can conclude that the alignment of the hairpin vortices is a result of the existence of a vortex-induction mechanism which causes the reorientation of vortical structure within the ring and braid skeleton. This phenomenon is in agreement with [START_REF] Brancher | Direct numerical simulations of round jets: Vortex induction and side jets[END_REF] for homogeneous jets, [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF] for the coaxial atomisation of a round liquid jet, [START_REF] Zandian | Understanding liquid-jet atomization cascades via vortex dynamics[END_REF] for coaxial atomisation of planar jets, and [START_REF] Bernal | Streamwise vortex structure in plane mixing layers[END_REF] for plane mixing layers. The analysis between the panels 'B' and 'C', can be extended to other locations of the jet where the core has undergone further interfacial development. For example, in panels 'D' and 'E' we observe the same distribution of inner and outer layer of hairpin vortex legs. The alignment of vortical structures has a detrimental effect on the interfacial dynamics, which will be shown in section 3.3.

Additionally, we have used the Q-criterion to present a three-dimensional visualisation of the vortices in the present study. The Q-criterion was described by [START_REF] Hunt | Eddies, streams, and convergence zones in turbulent flows[END_REF] as a quantity which measures the dominance of vorticity ω over strain s, Q = 1/2( ω 2s 2 ). Figure 8 shows the spatial development of the coherent structures for Q = 0.1. Near the injection point, axisymmetric KH vortex rings are located close to the surface (a magnified view is shown). As explained previously, when ω x and ω θ are of comparable magnitude, the KH vortices are stretched downstream or upstream. The topological shape of these vortical structures resembles the instantaneous hairpin-like vortical structures reported in experiments and numerical simulations by [START_REF] Head | New aspects of turbulent boundary-layer structure[END_REF]; To conclude, we have identified three different zones in the transient flow field which are associated with different vortical dynamics. We present the different regions at t = 28.97 (see figure 8). Zone 1 starts from the injection point and extends up to the deformation of the free-surface (x ∼ 1.8). This zone is characterised by the dominance of ω θ . Zone 2, which extends from x ∼ 1.8 up to x ∼ 11.5 (e.g., behind the mushroom-like structure). This region is characterised by the interfacial deformation of the jet-core by KH vortex rings and their posterior deformation by the competition between ω x and ω θ . As the flow moves downstream, ω x becomes responsible for the entrainment of the stagnant phase to form interfacial lobes (this agrees with [START_REF] Liepmann | The role of streamwise vorticity in the near-field entrainment of round jets[END_REF] for homogeneous jets). In this region, vortex rings pair up and merge together. Additionally, the ending of the visible-potential core of the jet and the beginning of the mixing region is also observed. Zone 3, expands from the end of zone 2 up to the leading edge of the mushroom-like structure. A large vortex-ring dominates the dynamics in this region (socalled 'cap-vortex', [START_REF] Zandian | Vorticity dynamics in a spatially developing liquid jet inside a co-flowing gas[END_REF]), wherein its interaction with upstream vortex structures via means of velocity induction leads to the formation of the neck, connecting the mushroom-like structure with the cylindrical body of the jet. The narrowing of the neck as a result of the streamwise convection of vortical structures is in agreement with the work of [START_REF] Asadi | On the scaling of propagation of periodically generated vortex rings[END_REF].

Cascade mechanism for droplet formation: hole-formation genesis

In the previous section, we discussed the vortex dynamics of the transient jet. This section focuses on the effect of vortices on the formation of inner/outer lobes, the thinning of which leads to genesis of droplets. Figure 9 shows a close view of the interfacial jet dynamics where outer lobes are present. The formation of lobes is observed along the spanwise direction of the jet. Close inspection of the interface shows it has assumed the shape of the surrounding hairpin vortical structures, where a change of the vorticity direction is observed for subsequent lobes (similar observations have been made by [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF]; [START_REF] Zandian | Understanding liquid-jet atomization cascades via vortex dynamics[END_REF]). Consequently, the alignment of the hairpin vortices induces the thinning of the outer lobes, leading to formation of holes, which expand and eventually give rise to droplets.

Figure 10 shows the temporal stretching of an outer lobe to form a ligament, and eventually droplets. The ligament orientation is linked to the vortical structures, as suggested above. A bulbous-tip is formed at the edge of the ligament driven by capillarity. The ligament posterior breaks at its base due to thinning then retracts under the action of surface tension forces to form droplets from both its up-and downstream ends an end-pinching mechanism [START_REF] Notz | Dynamics and breakup of a contracting liquid filament[END_REF]. The characteristic radial length scale (8.84, 9.96, 10.10, 10.92, 11.36, 11.41, 12.06, 13.61, 16.38), respectively. (l) Illustration of the entrapped-droplets inside of the mushroom-like structure at t = 7.33 from the back (probe 'A') and front (probe 'B') of the structure, respectively. A magnified region of the structure is also shown to illustrate the hole-formation behind the rim-edge. of the ligament is influenced by the external velocity field, the ratio of densities, and the surface tension. The ligament is characterised by having a small thickness prior to its detachment from the base. However, the droplets formed after the breakup process have a larger length scale (in agreement with [START_REF] Villermaux | Fragmentation[END_REF]).

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) A B (l)
Figure 11 shows the formation of droplets inside the jet core, and the mushroom structure. Specifically, figures 11b-k depict the formation and elongation of the inner lobe in the streamwise direction. As before, the stretching of the interface gives rise to a thin sheet as a result of the mutual induction of neighbouring hairpin-vortices. The hairpin vortex located above the liquid sheet induces flow downwards, whereas the vortex located under the sheet induces flow upwards. Under the joint action of the vortical structures, the liquid sheet thins until it is perforated (see figure 11e). Following the formation of the hole, retraction of the liquid sheet is radially driven by capillarity (see figures 11f and 11g), and as shown in figure 11i, it subsequently gives rise to the formation of a curved liquid thread or ligament. Once more, ligament retraction is driven by capillarity where bulbous ends that form initially undergo 'end-pinching' to form smaller droplets, as shown in figures 11j-k. The inertia-induced mechanism for the sheet-thinning presented here has been previously reported by [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF].

The analysis will now focus on the dynamics underneath the mushroom-like structure. Figure 11 shows entrapped droplets inside of the leading-edge structure. The formation of these droplets is also linked to the interfacial rupture of the toroidal liquid sheet. The liquid sheet retracts driven by surface tension to accumulate liquid in a rim at its edge. The rim experiences a spanwise destabilisation (i.e., Rayleigh-Plateau type), which leads to a non-uniform rim-radius. The film retraction is also accompanied by the formation of interfacial capillary waves that precede the rim. The capillary waves vary the film thickness, and consequently induce the perforation of the film adjacent to the rim in regions where the film thickness is sufficiently small (in agreement with [START_REF] Mirjalili | High fidelity simulations of microbubble shedding from retracting thin gas films in the context of liquid-liquid impact[END_REF]).The radial expansion of multiple holes yields the formation of liquid threads which experience a capillary instability to produce droplets. Following their formation, the entrapped droplets rotate inside the leading structure due to their interaction with the vorticity field (not shown here). The rotation leads to complex interfacial phenomena such as coalescence or collision (similar to that reported by [START_REF] Desjardins | Detailed numerical investigation of turbulent atomization of liquid jets[END_REF]). The coalescence has been observed not only between droplets, but also between droplets and ligaments or droplets and the jet core (see figure 12). Next, we aim to explain the validity of the hole formation mechanism as a result of liquid film thinning. Our numerical method does not consider the destabilising effect brought about by the van der Waals attractive forces as the film thickness tends to zero. [START_REF] Churaev | Surface Forces[END_REF] reported that while the thinning takes place, the dynamics of the film enters into an asymptotic behaviour leading to its rupture, depending only on a balance between viscous, van der Waals, and surface tension forces. After the film puncture, the circular hole expansion is driven by capillarity, and fluid is accumulated in a circular rim which grows in size as it moves away from the puncture point. The retraction speed, V T C can be estimated by the classical 'Taylor-Culick' theory, V T C = (2σ/ρh) 1/2 , where h is the film thickness. Prior to the puncture h ∼ 40 µm, giving an estimate of the retraction velocity as V T C ∼ 1.2 m/s. The measured retraction speed of the holes in our simulations ranges from 0.75 m/s to 1.2 m/s. This ensures that although our simulations cannot predict the exact location of the puncture, and their formation is mesh-dependent, its expansion is well predicted, ensuring that the physics is fully-resolved following the formation of the hole. Additionally, the film thickness of the sheet prior to its rupture is on the order of h ∼ 10 -6 m, which agrees with the findings of [START_REF] Marston | Crown sealing and buckling instability during water entry of spheres[END_REF], [START_REF] Lhuissier | Soap films burst like flapping flags[END_REF][START_REF] Lhuissier | Soap films burst like flapping flags[END_REF][START_REF] Néel | The spontaneous puncture of thick liquid films[END_REF]Villermaux (2018) who observed interface rupture for film sheets of the same order of magnitude.

Droplet size distribution

This section draws attention to the size distribution of the droplets, formed as a result of the ruptures of the liquid threads. Attention is focused on Case 4, which is in the inertia-dominated regime characterised by high Re and W e; as discussed above, the flow in this case is accompanied by significant drop creation. During the early stages, the formation of droplets is only observed inside the leading-edge structure owing to film rupture via hole formation. At later times, both entrapped and dispersed droplets coexist in the computational domain. We have identified all the droplets inside the entire domain, and each droplet diameter is calculated through knowledge of its volume, and the assumption of a spherical shape.

Figure 13 shows the probability density function (pdf) of the entrapped and dispersed droplets normalised by the mean droplet diameter at different temporal stages: t = 9.37, 20.37 and 28.52. The shape of the distribution for the entrapped droplets (see figure 13a) remains largely unchanged with increasing time except for the development of a tail; this represents the creation of larger numbers of smaller droplets due to ligament breakup as explained previously. It is possible to obtain information regarding the characteristics of the droplets in the domain; for instance, at t = 9.37, the mean and median diameters of the entrapped droplets are approximately 342 µm and 246 µm, respectively, with a standard deviation of 296 µm. Figure 13b shows that the pdf for the injected droplets has a similar shape to that of the entrapped ones depicted in 13a, and the distribution remains essentially unaltered for t = 20.37 and 28.52.

Conclusions

Three-dimensional numerical simulations using a hybrid interface-tracking/level-set approach were carried out for turbulent water jets entering a stagnant and more viscous silicon oil phase. Particular attention has been paid to the temporal interfacial dynamics that arise as a result of the vortex-surface interaction. Using less computationally expensive simulations, we explored the Re -W e space, identifying four distinct regions based on the appearance of the interfacial structures. On this basis, we carried out high-fidelity simulations representative of each region which showed that the streamwise vorticity plays a major role in the development of the three-dimensional instabilities on the jet surface.

At low Re and W e numbers (i.e., in the capillary-controlled regime), the jet interface remains axisymmetric and no surface corrugations are observed due to strong capillarity. Moreover, this region is also characterised by the fact that the streamwise vorticity, ω x , never becomes comparable to the azimuthal vorticity, ω θ : the azimuthal vorticity always remains two orders of magnitude larger than the streamwise vorticity. At low Re and high W e numbers, the reduction of capillarity enhances the formation of a mushroom-like structure at the leading edge of the jet, and the formation of corrugations in the jet core; however, ω x is still almost two orders of magnitude smaller than ω θ . At high Re and low W e numbers, streamwise vorticity gains in importance resulting in further deformation of the jet core.

In the inertia-controlled regime (i.e., for high Re and W e numbers), the streamwise vorticity becomes comparable to the azimuthal vorticity triggering the streamwise deformation of KH vortices. We show the formation of hairpin vortices near the interface, which are aligned in the streamwise direction forming layers of inner-and outer-hairpinvortices. The formation of inner and outer lobes in the jet core are closely linked to their neighbouring hairpin-vortices. The alignment of vortices enhances the perforation of the lobes to form holes which expand radially by capillarity, and ultimately give rise to the formation of droplets. Another main feature of the inertia-controlled regime is the creation of a thin toroidal sheet inside of the leading-edge structure. The thinning of this sheet leads to its perforation to form droplets. The droplets rotate as a result of their interaction with the vorticity field, and further topological transitions occur (e.g., coalescence).

While in the present study we have focused on the spatial and temporal development of the interfacial dynamics as a result of the vortex-surface interaction, further investigations should be carried out on the statistical response of the turbulent multiphase flow, similar to the work presented by [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF] for a mixing layer between parallel gas and liquid streams. Additionally, our study assumes a constant value of surface tension, but it is known that streams are usually contaminated with deliberately-placed or naturally-occurring surfactants, which reduce the surface tension and give rise to surface tension gradients and Marangoni stresses. As shown recently by Constante-Amores et al. et al. (2017, 2019) performed simulations of a two-phase mixing layer between parallel gas and liquid streams to investigate the interfacial dynamics and the statistics of the multiphase turbulence, estimating the Kolmogorov length scale to be η ∼ 0.945µm, which leads to ∆x/η ∼ 3. Through their simulations, they estimated that η = 3 -4.5µm and also showed that their lower resolution mesh, i.e., ∆x/η ∼ 6 (see their figure 19d), was capable of predicting similar results in terms of turbulence dissipation. For our highest Reynolds number (case-5), our simulation does not meet Pope's criterion. But as shown by [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF], the actual η could be larger. Additionally, the atomisation of the injected phase in a stagnant viscous phase alleviates the range of relevant physical scales.

The second biggest challenge of computational atomisation is 'numerical breakup' of liquid threads. A coarse grid would trigger the formation of thicker numerical threads, and consequently the formation of larger droplets (this problem has been previously reported by [START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF][START_REF] Herrmann | On simulating primary atomization using the refined level set grid method[END_REF], [START_REF] Gorokhovski | Modeling primary atomization[END_REF][START_REF] Gorokhovski | Modeling primary atomization[END_REF]Jarrahbashi &[START_REF] Jarrahbashi | Vorticity dynamics for transient high-pressure liquid injection[END_REF]). [START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF] stated that the mesh should be refined up to the point where the dynamics of the thread are solely governed by surface tension forces. These forces would trigger the formation of capillary waves during the retraction of these ligaments, giving rise to the onset of the Rayleigh-Plateau instability (i.e., the 'end-pinching' mechanism). After those mechanisms are initialised, the threaddynamics enter in an asymptotic behaviour towards the interfacial singularity (i.e., a refined mesh would not affect the size of the resulting droplets). On this basis, the numerical resolution regarding the interfacial length scales will be assessed following the methodology proposed by [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF] and [START_REF] Desjardins | Detailed numerical investigation of turbulent atomization of liquid jets[END_REF], who used a 'grid-based Weber number', defined by W e ∆xmin = ρ w U 2 ∆x min /σ. This equation provides us with the smallest interfacial length scale that the simulation is capable of resolving by assuming that the smallest interfacial structure is equal to the minimum mesh size of the computational domain. [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF] suggested that no further breakup is observed for values under 10. For a Reynolds number corresponding to Re ∼ 10 4 and M3-type mesh, the grid-based Weber number is W e ∆xmin ∼ 4.12, which meets the above criterion, suggesting that all capillary singularities would be resolved.

Therefore, we have proved that the M3 mesh is capable of resolving turbulent-scales and interfacial singularities, and consequently detailed analysis of interfacial and vortical structures is performed using a M3 mesh-type (unless stated otherwise). Figure 14 shows the temporal evolution of the kinetic energy, E k = V (ρu 2 )/2dv, the interfacial area and the maximum axial location of the jet tip (i.e., leading edge) for different meshes. Additionally, the vorticity profiles were checked between the M2-and M3-mesh, and no significant differences were found.

Scalings laws for the capillary breakup of liquid threads

As we have presented in the introduction, direct numerical simulations must be capable of predicting the developing two-phase fluid interfacial dynamics featuring interface breakup, and droplet coalescence. In light of this, we aim to show the capabilities of our numerical framework in the prediction of the scaling laws for the capillary singularity of liquid threads. As the point of singularity approaches, the system is driven locally by the large interfacial curvature, and its interfacial dynamics depends solely upon the physical properties of the liquid. [START_REF] Lister | Capillary breakup of a viscous thread surrounded by another viscous fluid[END_REF] have suggested that the pinchoff of a viscous thread (of radius r(z, t), density ρ, viscosity µ and surface tension σ) surrounded by another viscous fluid transits between different dynamical regimes as r(t) → 0 towards the breakup time τ (see figure 1 of [START_REF] Lister | Capillary breakup of a viscous thread surrounded by another viscous fluid[END_REF]). The thinning of a water thread, surrounded by air, transitions from an inertial-capillary regime (r ∼ τ 2/3 and u ∼ τ -1/3 ) to an inertial-viscous regime (r ∼ τ and u ∼ τ -1/2 ), more details can be found in [START_REF] Day | Self-similar capillary pinchoff of an inviscid fluid[END_REF], [START_REF] Eggers | Universal pinching of 3d axisymmetric free-surface flow[END_REF] and [START_REF] Lister | Capillary breakup of a viscous thread surrounded by another viscous fluid[END_REF].

Therefore, we consider the capillary singularity of a water thread (of initial radius R) surrounded by air in absence of gravity as the interfacial dynamics are locally driven by a balance of viscous-capillary and inertial forces, which can be expressed by the Ohnesorge number Oh = µ/(ρσR). As time evolves, the thinning of the thread is driven by capillarity, and the predicted thread radius towards the singularity agrees with the theoretical scalings (see figure 15) as the thread transits between different regimes (i.e., r(t) → 0). After the singularity point, the formation of a satellite ligament is observed which has an initially cylindrical shape (with Oh = 2.62×10 -3 and L o = 11.55); however, it undergoes retraction driven by capillarity to form more spherical droplets. This process known as 'end-pinching' has been previously well described by [START_REF] Notz | Dynamics and breakup of a contracting liquid filament[END_REF].

By performing this analysis we have proved the capability of our numerical technique to predict the dynamics of the capillary singularity of a liquid thread and its post-breakup events (table 5 shows a summary of the different meshes evaluated for this study). Additionally, the accuracy and validation of the numerical method has been previously addressed to other complex interfacial phenomena. These phenomena include breakup and recoiling of liquid threads (Constante-Amores et al. 2020b), falling film flows (Batchvarov et al. 2020a), propagation of elongated bubbles in channels (Batchvarov et al. 2020b), bubbles undergoing bursting (Constante-Amores et al. 2020a), and drops coalescing partially or completely with deformable interfaces.

Time step

Finally, the temporal integration scheme is based on a second-order Gear method, with implicit solution of the viscous terms of the velocity components. The time step ∆t is adaptive to ensure stability, and it is defined through the criterion: ∆t = min {∆t cap , ∆t vis , ∆t CF L , ∆t int } (4.1)

where ∆t cap , ∆t vis , ∆t CF L , ∆t int stand for the capillary time step, the viscous time step, the Courant-Friedrichs-Lewy (CFL) time step, and interfacial CFL time step, respectively. Those terms are defined by ∆t vis = min where ∆x min = min j (∆x j ). In our simulations the adaptive time-step is restricted by ∆t int (e.g., O(10 -5 )s).

Figure 1 :

 1 Figure 1: Injection of a turbulent water jet into a stagnant silicone-oil phase. Threedimensional representation of the interfacial shape for Case 4 of table 2 (i.e., Re = 6530 and W e = 303) at t = 28.97 with definitions of particular regions and specific features discussed in this work.

Figure 2 :Figure 3 :

 23 Figure 2: Regime map of the phenomenological interfacial dynamics in the Re-W e space using the M 2-mesh (see table4in Appendix for mesh details). Four different regimes and their boundaries are identified. A snapshot of the flow corresponding to the threedimensional representation of the interface for each simulated point (i.e., black marker) is shown. The red squares refer to the cases presented in table 2.

Figure 4 :Figure 5 :

 45 Figure 4: Two-dimensional representation of the interfacial location together with ω θ (top panel) and ω x (bottom panel) in the xz plane (y = 2) for Case 2 in table 2 (Re = 1000 and W e = 100) at t = 26.87. The colour represents the respective vorticity fields, where appropriate scales are shown in each panel.

Figure 6 :

 6 Figure 6: Two-dimensional representation of the spatio-temporal evolution of the interface for Case 4 in table 2 (Re = 6530 and W e = 303) together with ω θ (top panel of each sub-figure) and ω x (bottom panel of each sub-figure) in the xz plane (y = 2) for the dimensionless times shown in each panel. The colour represents the respective vorticity fields, where appropriate scales are shown in each panel.

Figure 7 :Figure 8 :

 78 Figure 7: Three-dimensional representation of the interface for Case 4 (Re = 6530 and W e = 303) at t = 28.97, showing the spatial locations 'A'-'E', (a). The streamwise location of the 'A'-'E' probe lines correspond to x = (1.50, 2.50, 2.75, 4.25, 5.25), respectively. (b) Velocity (top panel) and vorticity profiles (middle panel) in the y = 2 plane for each probe location. (c) Streamwise vorticity in the yz plane for each sampling location. The arrows show examples of identified hairpin-vortex legs. Solid and dashed arrows correspond to inner and outer hairpin vortex legs, respectively. The colour represents the streamwise vorticity field, ω x .

Figure 9 :

 9 Figure 9: Three-dimensional representation of the surface of the jet for Case 4 (Re = 6530 and W e = 303) at t = 28.97 showing the spatial formation of outer lobes. The colour represents the streamwise vorticity field, ω x .

  Figure 10: Three-dimensional representation of the interface showing the spatio-temporal evolution of the formation of droplets via the stretching of the outer lobe for Case 4 (Re = 6530 and W e = 303). Dimensionless times are shown in each panel.

Figure 11 :

 11 Figure 11: (a) Three-dimensional representation of the interface for Case 4 (Re = 6530 and W e = 303) at t = 7.33 showing a selected region in the jet core where inner lobes will form, together with two probe locations in the leadingedge structure, where the arrows indicate the direction of the view. (b-k) Cascade mechanism for the formation of entrapped droplets within the jet core at t =(8.84, 9.96, 10.10, 10.92, 11.36, 11.41, 12.06, 13.61, 16.38), respectively. (l) Illustration of the entrapped-droplets inside of the mushroom-like structure at t = 7.33 from the back (probe 'A') and front (probe 'B') of the structure, respectively. A magnified region of the structure is also shown to illustrate the hole-formation behind the rim-edge.

  Figure 12: Three-dimensional representation of the interface for Case 4 (Re = 6530 and W e = 303) at dimensionless times shown in each panel. A magnified region of the structure is also shown to illustrate droplet-droplet coalescence and the droplet-rotation inside of the leading-edge structure.

Figure 13 :

 13 Figure 13: Probability density function (p.d.f.) for the entrapped and injected droplets for Case 4 (Re = 6530 and W e = 303) at t = 9.37, 20.37, and 28.52 shown in panels (a) and (b), respectively. The distributions are normalised by the calculated mean diameter for each time, d m .

  (2020b), surfactants are capable of inhibiting capillary singularities and rigidifying the interface, subsequently changing the fate of the atomisation. Hence, the presence of surfactants may change the present results, providing us with an exciting avenue of future research. Declaration of Interests. The authors report no conflict of interest.

Figure 14 :

 14 Figure 14: Mesh study for Case 4 (Re = 6530 and W e = 303). The panels highlight the temporal evolution of the kinetic energy E k (left panel), the interfacial area (middle panel) and the maximum axial location of the jet tip (right panel).

Figure 15 :

 15 Figure15: Scaling laws for the capillary singularity of a water-thread surrounded by air for different numerical resolutions (see table5). Here, the Ohnesorge number is 1.178 × 10 -3 . The minimum thread radius, (a), and the maximum streamwise velocity, (b), versus the remaining time to breakup τ , agree with the inertial-capillary and inertial-viscous regimes presented by[START_REF] Eggers | Universal pinching of 3d axisymmetric free-surface flow[END_REF];[START_REF] Day | Self-similar capillary pinchoff of an inviscid fluid[END_REF]. Additionally, a three-dimensional representation of the interface is shown to depict the retraction of the satellite ligament to produce daughter droplets via the 'end-pinching' mechanism (in agreement with[START_REF] Notz | Dynamics and breakup of a contracting liquid filament[END_REF]).

Table 1 :

 1 Density and viscosity of the fluids used throughout this work[START_REF] Ibarra | Horizontal and low-inclination oil-water flow investigations using laser-based diagnostic techniques[END_REF]).

	ρ w kg/m 3 kg/m 3 mP a.s mP a.s ρ so µ w µ so
	998	824	1.0	5.4

Table 2 :

 2 The Reynolds-Weber number combinations for the four cases studied in detail with the M3-mesh (see table 4 in the Appendix).

	Case Re W e
	1 1000 30
	2 1000 100
	3 3260 75
	4 6530 303

Table 3 :

 3 List of computational studies of atomisation showing compliance with the Pope criterion[START_REF] Pope | Turbulent Flows[END_REF].

	Study	Configuration Type of study	Re	∆x/l	∆x/η
	Desjardins & Pitsch (2010)	Planar	Temporal	2000 -3000	0.015	∼ 5 -6
	Herrmann (2011)	Round	Spatial	5000	0.031/0.0078 ∼ 5 -19
	Jarrahbashi et al. (2016)	Round	Temporal	1600 -16000	0.0125	∼ 3 -18
	Zandian et al. (2016, 2018)	Planar	Temporal	2500 -5000	0.025	∼ 9 -15
	Zandian et al. (2019)	Round	Spatial	2000 -3200	0.01	∼ 3 -4
	Ling et al. (2019)	Planar	Temporal	8000	0.0039	∼ 3
	Current study	Round	Spatial	1000 -6530	0.0065	∼ 1 -4.5

Table 4 :

 4 Characteristics of different mesh sizes used to study the jet dynamics in this study.

	Run	Global mesh size (number of cells)	Number of parallel process threads	Minimum mesh size (µm)	Total Comput. hours
	M1	768 × 192 × 192	3 × 3 × 12 = 72	104.1	180
	M2	1536 × 192 × 192	6 × 6 × 24 = 864	58.5	100
	M3	3072 × 384 × 384	48 × 6 × 6 = 1728	26.1	500

state-of-the-art in terms of Pope's criterion. It is evident that Pope's criterion is difficult to meet due to the high computational costs of the simulations. In their recent work, Ling

Table 5 :

 5 Characteristics of different mesh sizes used to study the capillary breakup of a water thread.

	Run	Global mesh size (number of cells)	Number of parallel process threads	Pinchoff time (s)	Total Comput. hours
	V1	64 × 64 × 256	16	0.375	∼ 2
	V2	96 × 96 × 384	54	0.372	∼ 20
	V3	192 × 192 × 768	432	0.369	∼ 48
		(a)		(b)

† Email address for correspondence: o.matar@imperial.ac.uk

Acknowledgements

This work is supported by the Engineering & Physical Sciences Research Council, United Kingdom, through a studentship for C.R C-A. in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (Grant No. EP/L015579/1), and through the EPSRC MEMPHIS (Grant No. EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants. The authors would like to acknowledge the funding and technical support from BP through the BP International Centre for Advanced Materials (BP-ICAM), which made this research possible. O.K.M. also acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics, and the PETRONAS Centre for Engineering of Multiphase Systems. D.J. and J.C. acknowledge support through computing time at the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) Grant No. 2020A0082B06721. We also acknowledge HPC facilities provided by the Research Computing Service (RCS) of Imperial College London for the computing time. The numerical simulations were performed with code BLUE (Shin et al. 2017) and the visualisations have been generated using ParaView.

Appendices

Mesh study and resolution considerations Solving the small scales of the atomisation phenomenon is a challenging process, and in order to ensure the physical validity of the results, we have assessed the grid dependence nature of our results by performing a mesh study for Case 4 (see table 4). The appropriate choice for the minimum grid length-scale will depend on either the vanishing interfacial singularity or the smallest turbulence length-scale (i.e., Kolmogorov-length-scale).

For single-phase jets, [START_REF] Pope | Turbulent Flows[END_REF] suggested that the Kolmogorov-length-scale η is resolved if ∆x/η 2.1, where ∆x represents the minimum computational cell size, and η is estimated by η = Re -3/4 l, where l is the length scale. Previous computational studies of temporal and spatial multiphase jets are presented in table 3, showing the current