Progress in the Smoothed Particle Hydrodynamics Method to Simulate and Post-process Numerical Simulations of Annular Airblast Atomizers

This paper illustrates recent progresses in the development of the smoothed particle hydrodynamics (SPH) method to simulate and post-process liquid spray generation. The simulation of a generic annular airblast atomizer is presented, in which a liquid sheet is fragmented by two concentric counter swirling air streams. The accent is put on how the SPH method can bridge the gap between the CAD geometry of a nozzle and its characterization, in terms of spray characteristics and dynamics. In addition, the Lagrangian nature of the SPH method allows to extract additional data to give further insight in the spraying process. First, the sequential breakup events can be tracked from one large liquid blob to very fine stable droplets. This is herein called the tree of fragmentation. From this tree of fragmentation, abstract quantities can be drawn such as the breakup activity and the fragmentation spectrum. Second, the Lagrangian coherent structures in the turbulent flow can be determined easily with the finite-time Lyapunov exponent (FTLE). The extraction of the FTLE is particularly feasible in the SPH framework. Finally, it is pointed out that there is no universal and ultimate non-dimensional number that can characterize airblast primary breakup. Depending on the field of interest, a non-dimensional number (e.g. Weber number) might be more appropriate than another one (e.g. momentum flux ratio) to characterize the regime, and vice versa.


Introduction
Liquid atomization has a wide range of applications, such as in agriculture with the dispersion of fertilizer, in automotive industry with coating deposition or in turbomachinery with fuel injection. Even though this phenomenon is widely used, the core mechanism is 1 3 not always well understood. The development of new nozzles always requires many iterations in the chain of design/manufacture/characterize processes. In addition, some applications involving liquid atomization, especially in the field of aeroengines and energy, require extreme conditions of pressure and temperature that strongly limit the instrumentation and the characterization of the nozzles. The numerical tool is a promising candidate to accelerate the design of new nozzles, in two manners. First, it can simulate configurations in extreme conditions, with access to all flow quantities in the whole numerical domain. Second, the increase of computational capabilities allows finer resolutions, thus allowing to resolve smaller and smaller length scales of atomization , providing insight in the phenomenon and increasing the scientific knowledge.
In this work, the focus is given on annular airblast atomizers. This type of atomizer is mainly used to inject the liquid fuel into the combustion chamber of gas turbines. Most of the literature concerning the numerical simulation of annular airblast atomizers based on first principles is dedicated to coaxial flows where the liquid is injected as an axial jet enclosed in a high-speed airstream. Moureau and Desjardins (2008) investigated this configuration by coupling accurate conservative level-set (ACLS) method ) to the flow solver by a second-order ghost fluid method (GFM). Xiao et al. (2014) performed a LES of the same type of configuration with a coupled level set (LS) and volume of fluid (VOF) method (CLSVOF), with a density ratio of 830. Ma et al. (2016) performed a LES of a hybrid airblast atomizer with a VOF/Lagrangian spray framework. Müller et al. (2017) used OpenFOAM with VOF to investigate the influence of the nozzle design on the breakup of a cylindrical viscous jet. Odier et al. (2018) investigated the flapping mechanism of a coaxial liquid jet at a moderate density ratio (10) with the ACLS method from Desjardins et al. (2008). To the authors knowledge, there is no publications on the simulation based on first principles, of airblast atomizers where the liquid is injected as an annular sheet. This is the configuration of the present study, which will be simulated with the use of the smoothed particle hydrodynamics (SPH) method.
The SPH method was originally developed for astrophysics by Gingold and Monaghan (1977) and Lucy (1977), and its applicability to simulate fluid dynamics cases was already stated in Monaghan (1982). A weakly-compressible approach was proposed by Monaghan (1994) to approximate incompressible free surface flows and finally, multiphase flows (Monaghan and Kocharyan 1995). Nowadays the SPH method is successfully applied in many fields of fluid mechanics such as free surface flows (Peng et al. 2017), turbulent flows (Monaghan 2011), fluid-structure interactions (Tofighi et al. 2015) and compressible flows (Federrath et al. 2010). An exhaustive review of SPH industrial applications is provided in Shadloo et al. (2016). The successful application of SPH to liquid atomization in various configuration was recently shown by Koch et al. (2017), Chaussonnet et al. (2019) and Braun et al. (2019). The SPH method, as a Lagrangian method, naturally convects fluid particles and their physical properties inside the numerical domain. Therefore the convective acceleration is not discretized but is inherently taken into account by the motion of the particles. In the simulation of multiphase flows with the SPH method, each phase is represented by particles of different types (e.g. liquid and gas particles), carrying the corresponding physical properties. All phases are computed by the same formalism. The phase interface is naturally taken into account by the motion of liquid particles compared to gaseous particles. Therefore no computation expense is necessary to track or reconstruct the interface.
This paper is organized as follows. The numerical method and its implementation aspects are presented in Sects. 2 and 3. These include the description of the numerical method, the conversion tools to generate a particle initial solution from a CAD model 1 3 and to post-process the SPH simulations. Then the geometry and the numerical setup are presented in Sects. 4 and 5. The single simulation is commented in Sect. 6. Sections 7 to 11 are dedicated to the analysis of the liquid phase.

Numerical Method
The fundamental principle of the SPH method is to express any field (mass, momentum and energy) and its derivatives at a given particle location (subscript a) with an interpolation from its neighbor particles (subscript b) (Monaghan 2005): where V b is the volume of the neighbor particles. The function W is referred to as the kernel and depends on the inter particle distance r b − r a and a characteristic length scale h called the smoothing length. It promotes the influence of closer neighbors (Fig. 1 top). The smoothing length defines a sphere of influence Ω a (Fig. 1 bottom) so that the neighbors located outside of it are not taken into account. In this study, the kernel is a quintic spline, given by: where r = |r a − r b |∕h is the normalized distance between particle a and b. The smoothing length h is equal to the initial mean particle spacing Δx, and the radius of Ω is R = 3h. In the following, f (r a ), f (r b ) and W(r b − r a , h) are abbreviated by f a , f b and W ab , respectively.
In the present study, the motion of the fluid phases (gas and liquid) is described by the weakly-compressible Navier-Stokes equations. The continuity is expressed algebraically by calculating the particle volume and density (Español and Revenga 2003): where m a is the constant mass of particle a. Equation 3 exactly conserve mass. This expression of the density shows the advantages (i) to conserve mass exactly, and (ii) to avoid any numerical diffusion of density at phases interface (Hu and Adams 2006). The conservation of linear momentum on the particle a is expressed by: where u is the particle velocity and the terms f p,a , f v,a , f a,g and f st,a are the forces due to pressure, viscosity, gravity and surface tension, respectively. They write: The viscous term (Eq. 5b, proposed by Cleary 1998) conserves locally both linear and angular momentum. The inter-particle viscosity ab is proposed by Szewc et al. (2012) and reads: where = ∕ is the kinematic viscosity. This formulation is in practice more stable for airblast atomization. Surface tension (Eq. 5d) is adapted from the continuum surface force (CSF) model (Brackbill et al. 1992), as proposed by Hu and Adams (2006) and improved by Adami et al. (2010). According to the weakly compressible assumption, the pressure p a at particle a is expressed by the equation of state originally proposed by Cole (1965): where 0 , and p back are the nominal particle density, the polytropic ratio and the background pressure, respectively. The artificial speed of sound (c) must satisfy c ⩾ 10 u max to respect the weakly compressible condition, ensuring a density variation lower than 1% (Monaghan 1994). In the present case the maximum Mach number of the gas phase is 0.38. Therefore, as the weakly compressible assumption is not perfectly fulfilled, the compressibility effects must be accounted for. This is possible in the present framework because Eq. 7 can describe an isothermal compressible flow (Colagrossi and Landrini 2003). Considering the gas as an ideal gas, the increase of temperature due to compressibility can be estimated  to 4.2%, which leads to variations of density and 1 3 viscosity of 4.2 and 1.7%, respectively. As this is below 5%, the flow is considered as isothermal, and the pressure is proportional to the density. Moreover, by defining the effective speed of sound as c eff = √ p∕ and using the equation of state defined in Eq. 7, it can be shown that: Therefore, is set to one and the artificial speed of sound is set to 340 m/s for the gas phase in order to retrieve the real speed of sound. For the liquid phase, the polytropic ratio l is classically set to 7 to mimic the incompressibility of water, and the speed of sound in the liquid is also set to 340 m/s to increase the time step and keep the Mach number low. A blending function is applied for the pressure gradient at the interface as in Chaussonnet et al. (2019) to damp numerical instabilities.
Equation 4 is integrated in time by a predictor-corrector scheme. The time step of the numerical simulation is controlled by convection, viscosity, surface tension and the overall acceleration. It is computed at each iteration step: where Γ a is the magnitude of the particle acceleration and C 1 to C 4 are constants equal to 0.5, 0.25, 0.25 and 0.5, respectively.
At the boundaries a special treatment is required to complete the sphere of influence of near-wall particles. Boundaries are classified in two categories, permeable boundaries such as inlet or outlet, and impermeable boundaries such as solid walls. Permeable boundaries are modeled with two additional types of particles (Fig. 2). First, markers are stationary virtual particles located at the physical prescribed boundary, and act as sensors and prescriptors of the flow conditions. They contain the physical boundary conditions (e.g. inlet velocity for inlet) as well as the normal vector. They record also the other flow conditions that are not imposed (e.g. pressure for inlet). Second, buffer particles are located in layers upstream the physical inlet, whose thickness dr is equal to the radius of the sphere of influence. Their purpose is to complete the sphere of influence of particles that are injected at the boundary. Source terms are not computed for these particles. In this buffer layer, Fig. 2 Buffer (gray) and regular (white) particles separated by markers (squares), from Braun et al. (2015) the particles physical properties are imposed by the markers, i.e. prescribed and recorded properties from the markers. They are convected in the direction of the flow. At the inlet, when they cross the virtual line of the markers, they become regular particles. At the outlet, when regular particles cross the markers line, they become buffer particles. They are deleted from memory when they exceed the distance of the sphere influence radius. An ad hoc method was developed to damp incoming pressure waves at permeable boundaries. It is thoroughly explained in Braun et al. (2015). For solid walls, several layers of stationary particles are disposed at the boundary in order to complete the sphere of influence, as depicted in Fig. 3 (left). Based on the approach of Takeda et al. (1994) and Morris et al. (1997), the no-slip condition is enforced during the calculation of the viscous terms: a virtual velocity is applied at wall particles, pointing in the opposite direction of the flow, thus leading to a shear-stress that decelerates the particle in the vicinity of the wall (Fig. 3 right). This methods relies on the assumption of the linear velocity distribution in the vicinity of the wall. It was found to be second order convergent (Maciá et al. 2011).

Numerical Infrastructure and Pre/Post-processing
In this section, the infrastructure of the code, the pre-processing, and the post processing strategy are presented.
The SPH code is written in C++ and uses the Message Passing Interface (MPI) library for the parallelization. The interprocessor communication is based on the exchange of halo particles located on the boundary of the subdomains via non-blocking communications (Braun et al. 2016). The performances were optimized at the core-level by the use of a data structure of arrays. The particles are tracked using a list-search algorithm based on a background grid with a complexity of O(n) where n is the number of particles. The particles are sorted in a manner to enhance the data locality, which promotes a cache-friendly data access ). The code can handle up to 10 billion particles and shows acceptable strong-scaling characteristics up to 10,000 cores (Braun et al. 2016).
Even if SPH does not require any meshing step, the numerical domain must be imported from a CAD file and converted into wall particles. The fluid volume is defined by subtracting the walls from the overall cavity volume with the CAD software. Then, it is exported as a tessellated surface mesh into MeshLab (Cignoni et al. 2008) where it is checked and Fig. 3 Sketch of no-slip boundary condition at a wall moving at u w possibly repaired. In order to create the layer of wall particles from the liquid volume, an additional non-intrusive surface mesh is created by imposing a shift normal to the walls from the initial mesh. The volume contained between the two surface meshes defines the wall volume. Thereafter, the two volume meshes are converted into particle lattices by a in-house code (CAD2SPH), where particles are placed on the mesh nodes. At this point, wall surfaces may be aliased due to the initial mesh. To smooth out the wall surface, a particle packing algorithm is applied, which reorders the wall particles in the constrained wall volume. The whole process is thoroughly explained in Dauch et al. (2017). The order of magnitude of the time to conduct the preprocessing workflow is one day, which is to be compared with the usual time to design and generate a mesh in a complex geometry. Concerning the post-processing, Eulerian fields can be extracted by projecting the particle variables on a stationary grid. This enables to compute time averages, which are appropriate for quantitative validation purposes. These can be the quantities directly carried by the particle (density, velocity, pressure) but also fields related to multiphase flows such as the liquid volume fraction. With respect to liquid atomization, clusters of particles that constitute liquid droplets are identified by a connected component labeling (CCL) technique (Rosenfeld and Pfaltz 1966), which gives access to diameter, velocity, position and sphericity of each Lagrangian droplets. Such fields can be also converted to Eulerian fields. For a more qualitative post-processing of the liquid phase, especially to be used in a rendering engine, the α-shape algorithm (Edelsbrunner et al. 1983) is applied to the particles of liquid type, which compute a tessellated surface of the liquid bulk. Note that this method leads to aliased surface and hence it is required to apply a spatial smoothing filter to obtain a more realistic rendering.
The Lagrangian nature of SPH, and the continuous tracking of each particles enables to extract information in a straightforward manner compared to traditional grid-based methods. The only requirement is to identify each particle by a unique ID number. Thanks to that, it is possible to tag each cluster of liquid particles (i.e. each liquid droplet) during the simulation. This enables us to record the historical connection between all droplets, especially during cascade processes such as breakup. This connection will be herein referred to as the tree of fragmentation. It is illustrated in Fig. 4.
From this tree of fragmentation, it is possible to determine the time and the location of the breakup event. We define N the breakup activity as the number of breakup events occurring in a volume dV and during a time dt. It is computed as: where B p is a breakup event.
Finally, the Lagrangian aspect of the SPH and the ID tracking allow also to extract Lagrangian coherent structures (LCS) to investigate the vortical structures of the flow. LCS (10) are distinct closed surfaces containing fluid particles which are coherent in time, i.e. whose deformation is lower than a given level. As shown by Haller (2005), the finite-time Lyaponov exponent (FTLE) is an objective quantity to identify LCS from a flow field. It is based on a flow map t t 0 that links any fluid particle located at the position x 0 at t 0 , to its position at t according to the motion of the flow. The gradient of t t 0 at x 0 gives an indication on the deformation of the continuum in the vicinity of x 0 between t 0 and t. The singular values of ∇ t t 0 , labeled i , are therefore an objective measure of the local expansion or contraction at x 0 between t 0 and t. The FTLE between t 0 and t, noted t t 0 (x), is based on the largest singular value n : For t > t 0 the loci of local maximum of t t 0 shows the region of large repulsion whereas for t 0 > t the loci of large values of t t 0 reveal the region of attraction of the flow. In the following, only the regions of attraction are shown. The core of FTLE computation is to determine the new location of a fluid particle after a period Δt = t 0 − t. With Eulerian methods, this is a challenging task that can be fulfilled, for example, by resolving the transport of a passive scalar. With Lagrangian methods, it is instantaneously achieved, provided that particles own a unique identifier. Extraction of FTLE with the SPH method was pioneered by Sun et al. (2016) and efficiently implemented using GPUs by Dauch et al. (2018).

Geometry
The numerical domain (Fig. 5) is made of an airblast atomizer discharging into a cylindrical cavity of length 15 mm and diameter 12.4 mm. The airblast atomizer (Fig. 6) injects an annular liquid sheet enclosed between two swirling flows. The flows are counter-rotating to increase the shearing in the liquid sheet and to avoid a blasted breakdown configuration of the central recirculation zone (Falese et al. 2014) which would impose to simulate a larger domain. The swirl splitter is added to simplify the boundary conditions for the swirl motion in two manners. First, it artificially creates a core for the vortex. Second, it provides a solid separation on the central axis, which avoids that the streams of opposite azimuthal directions meet. This would lead to numerical instabilities otherwise. At the exit of the second stage, the outer lip is curved towards the center to counteract the centrifugal effects and to force the gas to impact the liquid sheet. Dimensions of the nozzle are given in Table 1 and in Fig. 6.

Numerical Setup
The boundary conditions are annotated in Figs. 5 and 6. A turbulent and laminar axial mean velocity u z profile is imposed for the gaseous and liquid ducts, respectively. No time fluctuations are added to the mean profile in case of turbulent inlet. For the gaseous stages of the nozzle, a circumferential component u is imposed to generate a swirling flow. It follows a Rankine vortex with a constant angular momentum C = r u (r) in the free flow and a linearly decreasing u (r) in the boundary layers. The component u is determined with the swirl number, defined as: where Q m is the mass flow rate and C = R out u (R out ) is the angular momentum taken at the outer radius of the stage. Capital letters (e.g. U z ) represent bulk quantities. An entrainment flow with a large velocity is imposed around the nozzle to enhance the convection of the liquid and minimize the liquid impact on the cavity wall, without triggering additional breakup events. At the outlet a free flow condition with a constant pressure of 11 bar is imposed. The elevated pressure has an impact on the air density only. Inlets and outlet velocities are corrected to damp pressure waves following the method presented in Braun et al. (2015). A no-slip boundary condition is imposed on the injector walls. In order to reduce the size of the computational domain, the walls of the cavity are located closer to the nozzle than in typical industrial configurations. Therefore, to limit the influence of the cavity walls, a slip boundary condition is imposed. The initial interparticle distance Δx is 33.33 μm leading to 68 millions particles. The simulations are run on 2000 cores for a physical time of 13 ms with a mean time step of 30 ns, which leads to 440,000 time steps and 391,000 CPU hours per simulation.
The properties of the gas are a density g and a dynamic viscosity g of 13.2 kg∕m 3 and 18.61 μPa s, respectively. The properties of the liquid mimic a Newtonian viscous slurry which highlights a large viscosity. Its characteristics are a density l , a dynamic viscosity l and a surface tension of 1233 kg∕m 3 , 0.2 μPa s and 63.6 g∕s 2 , respectively. One single-flow and two two-phase flows will be investigated in the present study, their boundary conditions are recalled in Table 2 where G stands for gas and L stands for liquid. The order of magnitude of the velocities (1 and 100 m/s for the liquid and gas, respectively) corresponds to the one found in aeronautical burners. Several non-dimensional groups are used to describe the breakup phenomenon. They can be sorted into two categories, local and global parameters. The local parameters are the Reynolds number, the Weber number, the Ohnesorge number and the momentum flux ratio M . The Reynolds number represents the competition between inertia and viscosity, and the Weber number represents the competition between inertia and surface tension. Note that the Weber number can be expressed based on the liquid inertia ( We l ) or the gas inertia ( We g ), as explained below. These numbers are expressed as: where D h,l is the hydraulic diameter of the liquid inlet equal to 2(R out − R in ) and Both of these numbers can be regarded as the ratio of destabilizing by stabilizing effects, even though they characterize different regimes. The Reynolds number delimits the transition between laminar and turbulent flows and can be equally applied to gas and liquid whereas the Weber number describes the stability of a liquid structure where the surface tension coefficient can be defined. As the liquid structure can be destabilized by its own inertia (droplet impact on walls or Rayleigh-Plateau instability) or by gas momentum transfer (airblast or air-assisted atomization), the Weber number can be expressed based on the liquid or gas inertia, as shown in Eq. 13. Note that the total velocity U g,tot and not only the axial velocity U z is considered for We g because the azimuthal component of the gas carries a significant amount of momentum and hence it plays a role in the atomization.
The momentum flux ratio is expressed as M = g U z U g,tot ∕ l U 2 l and represents the local amount of gas momentum available to atomize a given mass flux of liquid of a given momentum. This parameter is recognized as one of the most important parameter to characterize airblast or air-assisted atomization. The Ohnesorge number is expressed as: and represents the combined effect of viscosity, surface tension and inertia. It is used to distinguish between conditions that lead to different modes of liquid jet breakup. Since no velocity appears in its expression, the Ohnesorge number only depends on the mechanical properties of the liquid structure and on its size. In the present study, the Ohnesorge number is equal to 1.01, based on the height of the liquid duct.
The global parameters are constructed by integrating local quantities over their flow section. The gas-to-liquid-ratio (GLR) is defined as the ratio of mass flow rate of the gas to the one of the liquid. It is expressed here as: The two terms in the numerator correspond to the two gas ducts of the present configuration. The GLR was identified as an important parameter in pioneering works on airblast atomization (Rizkalla and Lefebvre 1975;El-Shanawany and Lefebvre 1980). Moreover, when chemical reactions occur downstream of the nozzle, the GLR also represents the global chemical equilibrium. The momentum flow rate ratio M F is similar to the (13) Re = U g D h and We l = l U 2 l D h,l and We g = g U 2 g,tot D h,g , momentum flux ratio M , except that the fluxes are integrated on the flow sections. It is thus expected to bear a global characterization of the atomization. It is expressed as: Note that due to the swirling flow, a significant amount of the azimuthal momentum is used to atomize the liquid. Therefore, we consider the flow rate through the exit section ( U z ) of the total momentum ( U g,tot ), and not only the axial momentum. The values of the nondimensional parameters depending on the case is given in Table 3. It is interesting to note that the Reynolds and Weber numbers are different between L 1 and L 2 whereas the momentum flux/flow rate ratios and the gas-to-liquid ratio are identical. In the following when we will compare different output fields of the simulation between L 1 and L 2 , similarities (resp. differences) in the fields will promote the importance of M , GLR, M F (resp. We l and We g ) to characterize the state of the system relative to the field.

Gaseous Phase
The single phase simulation is presented in this part. The streamlines issuing from a vertical line at the inlet are shown in Fig. 7. They are colored by u to highlight the two counter rotating flows, whose angular momentum decreases towards the domain outlet as they mix together. As a consequence, the swirl jet is not in a blasted configuration and a weak precessing vortex core (PVC) is observed. This will reduce the liquid impaction on the cavity wall.
(16) M F = g (U z,1 U g,tot,1 A 1 + U z,2 U g,tot,2 A 2 ) l U 2 l A l . The time statistics of the axial velocity on the center plane are shown in Fig. 8. The opening of the jet due to centrifugal effect is observed, and leads to the typical central recirculation zone (CRZ) which is somewhat smaller here than compared to the experimental results from Vashahi et al. (2017). The most probable reason is that in our case the swirl numbers of each vanes are exactly opposite (1 and −1 ) whereas they are 0.767 and −0.911 in Vashahi et al. (2017). Hence in our configuration the relative flux of angular momentum (given by the swirl number) is balanced between the two vanes, and ideally the CRZ should disappear. The RMS of the axial velocity is large in the shear layer between the two counter rotating streams, and on the boundary of the CRZ. The mean and RMS fields of the circumferential velocity on the center plane are depicted in Fig. 9. The mean of u (Fig. 9 top) shows that the circumferential velocity of the two stages sharply decreases as they mix together. The RMS (Fig. 9 bottom) of u also show large fluctuations in the shear layer between the two streams, at the outer part of the outer stream and also at the front of the CRZ. Finally, the axial velocity was measured 3 mm downstream the nozzle at two locations: on the central axis (i.e. in the CRZ) and at a radius of 1.3 mm. Their cross-spectral density (CSD) was derived, using zero padding to increase the frequency resolution to 22 Hz. They are reported in Fig. 10. The largest peak is located at f 1 ≈ 6037 Hz and corresponds to the frequency of the PVC.

Breakup Sequence of the Liquid Sheet
Time series of the liquid breakup occurring at the atomizer lip are shown in Fig. 11. Since the gas velocity is larger in case L 2 than in case L 1 , the breakup frequency is also found larger in L 2 . Hence, the two time series in Fig. 11 are shown with a different time increment. It is to be highlighted that the volume of the liquid in Fig. 11 is significantly decreased due to the conversion from particles to an Eulerian smoothed surface, as explained in Chaussonnet et al. (2019). Nevertheless, the breakup mechanism is well observable in both cases. It occurs as follows, with annotation in Fig. 11 (left). Due to the curved lip of the outer gas stage, the flow is deviated towards the central axis (a), which in turns, deviates the liquid issuing from the duct towards the axis (b). This reduces the flow section of the inner stage (c) and tends to locally accelerate the inner flow, leading to a larger shear of the liquid. The liquid sheet then starts to fragment itself (d) and is pushed in the radial outer direction (e). This leads to ligaments bent towards larger radii, and to a liquid sheet in a shape of a trumpet (f). This liquid sheet is then pinched off by the stream of the outer stage (g) that locally increases the dynamic pressure on the liquid surface. The torn liquid blobs are then further fragmented by the gas in streamwise stretched ligaments. This mechanism leads to a flapping motion whose frequency depends on, among all, the gas velocity.
To strengthen this scenario, the time signals of the axial velocity was measured at the exit of the injectors in the two different stages, and their CSD was derived. They are plotted in Figs. 12 and 13 for L 1 and L 2 , respectively. Windowing (Blackman) and zero padding were applied, the frequency resolution is 25 Hz. A peak in the CSD is observed at 4003 and 7780 Hz, which leads to a characteristic time of 250 and 158 ns, which corresponds to the time series of breakup event depicted in Fig. 11. Note that in the case of L 2 , the frequency of the PVC is not visible in the spectrum. The time cross-correlation of the axial velocity (not depicted here) is zero for L 1 , meaning that the velocities in the two stages are not correlated, whereas for L 2 , the cross-correlation is ≈ 50 ns, which roughly corresponds to the half time of a breakup sequence. This is consistent with the aforementioned scenario in the conditions of L 2 .
Finally, the turbulent intensity u � ∕U z of the gas is reported in Fig. 14, superimposed with the iso contour of the zero axial velocity. L 1 and L 2 show a similar qualitative map, with a maximum of fluctuation between the liquid sheet and the outer stage. In the case of L 2 , the influence of the liquid phase on the gaseous phase can be assessed by a qualitative comparison with Figs. 8 and 9. The most striking feature is the disappearance of the CRZ, which is due to the reduction of the swirl of stage 1 due to momentum transfer. This is shown later. Also, the field of turbulent intensity is more homogenous with the presence of liquid. Indeed, the flapping motion of the liquid sheet leads to strong flow deflections that broaden the distribution of fluctuations.

Momentum Transfer Between Liquid and Gas
The time average of the axial momentum mean flux along the axis was computed for the gas and liquid phases, to highlight the global exchange of momentum between phases. It was obtained by integrating u 2 over time and over slices perpendicular to the main axis every millimeter on the z-axis. It thus intended to give a global analysis of the momentum transfer at each axial position without focus on the local exchanges. They are shown in Fig. 15. For the case L 2 , the profile from Case G is also given. Note that it is more computationally efficient to compute the momentum flux over a surface than integrating the momentum over a volume. First, the total axial momentum decrease along the axis, which correspond to dissipation of the momentum by friction. For both operating points, the increase/decrease of the liquid/gas axial momentum clearly depict the momentum transfer from the gas to the liquid. However, it is observed that the momentum flow rate of the liquid reaches a maximum and eventually decreases. This is a consequence of the radial deflection of the liquid blobs/droplets due the flapping effect. Indeed, when these blobs/droplets of moderate velocity reach larger radii where the axial velocity is lower, they are decelerated due to drag, thus decreasing the liquid momentum flux. For L 1 , the liquid momentum flux increases smoothly until z = 7 mm whereas it increases more steeply until z = 3 mm for L 2 . This suggests that the characteristic length L MT over which the most of the momentum is transferred depends on the case. However, as seen later, this length of Fig. 14 Turbulent intensity u � ∕U z for L 1 (top) and L 2 (bottom) superimposed with the contour of zero axial velocity. Walls are indicated by black dashed lines Fig. 15 Axial profile of the axial momentum for L 1 (left) and L 2 (right) momentum transfer is not correlated to the position at which the spray is fully formed. This is observed on the videos of the simulation and confirmed later. Also, the residual difference between the gas and liquid momentum flux at the domain exit depicts the efficiency of the momentum transfer. This difference is rather independent of the case, equal to ≈ 0.1 N. Finally, the axial momentum flux in the single phase simulation (Fig. 15 right) also significantly decreases whereas there is no momentum transfer to the liquid phase. First, this is due the high shear stress between the two counter rotating channels that leads to intense friction effects. Second, the CRZ reduce the flow section in the cavity, thus leading locally to high velocity zones where dissipation plays a larger role.
The angular momentum flux along the axis is given in Fig. 16. The transfer of momentum is also visible, somewhat smaller than the axial momentum transfer. This is due to the cylindrical geometry of the case and to the ballistic trajectory of the liquid blobs. Indeed, when a liquid blob moves in a straight line above the central axis, a part of its azimuthal momentum becomes accounted as a radial momentum. This leads to decrease the angular momentum for the single blob. Also, the same difference of the length of momentum transfer L MT is observed in the case of angular momentum with a sharper increase of liquid angular momentum flux for L 2 . The angular momentum flux of the single phase simulation is equal to the total angular momentum flux of case L 2 because the liquid is injected with no angular momentum. A significant amount of gas angular momentum is transferred to the liquid. In turn, the swirl motion becomes too weak to trigger the recirculation zone.

Spray Quantities
When the liquid blobs are torn apart from the sheet, they are identified by our post-processing tool ) as a detached liquid structure, and therefore their position, velocity, equivalent diameter and sphericity are monitored. Hence, it is possible to draw a map of the spray quantities inside the cavity. Figures 17, 18 and 19 show timeaveraged spray quantities. On these figure the (x, y) axes correspond to the (z, r) axis of the geometry, i.e. radial versus axial coordinates. The fields are obtained by averaging the quantity of interest over an elementary ring of radius r at position z and of cross section drdz over all the time steps (see Fig. 18 in Chaussonnet et al. 2019). Because the sampling time to build these fields was much smaller than the typical flow through time, the velocity of the liquid was taken into account in the average construction in order to increase the weight of fast droplets (which appear on fewer frames) versus slow droplets. In these figures, the origin of the map corresponds to the center of the injector exit. Fig. 16 Axial profile of the angular momentum for L 1 (left) and L 2 (right) Figure 17 shows the mean liquid fraction α L of the liquid phase downstream the nozzle exit for L 1 (left) and L 2 (right). The global spatial distribution is almost identical in both cases. This could be surprising at first sight because the liquid mass flow rate is twice larger in L 2 than in L 1 , and hence a larger concentration could be expected. But since the gas velocity is also twice larger, the liquid structures are almost accelerated by a factor of two, thus leading to cloud of droplets more expanded in space, which leads to approximately the same volume fraction. The opening of the spray is ≈ 54 • until z = 3 mm, then ≈ 40 • . The dense regime ( α L > 10 -3 ) extends up to z = 6 mm in both cases. Comparing these differences with the non-dimensional numbers defined in Sect. 5, the Weber number has no influence on the spatial distribution of the liquid volume fraction, and the momentum flux or flow rate seems a better indicator.
Then, the droplet concentration is shown in Fig. 18. In this case as well the global spatial distribution is the same between L 1 and L 2 . In the near-nozzle region ( z < 6 mm), case L 2 shows a larger concentration of droplet which suggests a better atomization. At the nozzle exit where the gas streams enter in contact with the liquid phase (z = 0 mm, r ≈ 1.8 mm) some droplets are already created. This comes from the high shear stress imposed on the liquid surface that leads to the striping of the liquid sheet and generates very small droplets. The droplet concentration at this location is already higher for L 2 than for L 1 . This is because the shear stress on the liquid surface depends on the relative velocity (much larger for L 2 ) and not on the velocity ratio. Hence, at this location, the Weber number is more representative than the momentum flux ratio to characterize the breakup regime.
The global SMD is presented in Fig. 19. At each sampling location (r, z), the total droplet volume V d is summed over all the time step, and divided by the total droplet surface S d summed over all time steps: First, the oversampling in-time is well visible with large droplets locally increasing the SMD and being convected towards the exit, highlighting their trajectory in dotted line. Second, the case L 1 produces larger droplets in the outer part of the spray. These large droplets are created by the ligament breakup, not directly at the nozzle exit, but later downstream at location ≈ (2, 2) in the (r, z) map. Then they follow a ballistic trajectory. The same type of droplets are visible with L 2 , but in a much smaller concentration, thus leading to a rather smaller SMD in the outer part of the spray. This difference between L 1 and L 2 suggests that the Weber number is more appropriate than momentum flux (or flow rate) ratio to characterize the primary breakup.
The spray properties were space-averaged on slices of 200 μm along the axial coordinate. A Fourier-transform (FT) was applied to each slice, so that the spectral properties of the spray quantities along the axial coordinate can be monitored. Figure 20 shows the evolution of the CSD between the SMD and the droplet concentration along z. The zone up to 2 mm shows a white noise because at this stage, most of the liquid is contained in the sheet and there are few droplets. For L 1 , a peak is observed at 4 kHz for z between 2 and 4 mm, and above 7 mm. This means that the spray quantities fluctuate at the frequency of the flapping of the liquid sheet. A peak is also observable at 8 kHz for z between 4 and 6 mm and could be a harmonic. For L 2 , a peak is found at the . Fig. 20 Maps of CSD between the SMD and droplet concentration for case L 1 (left) and L 2 (right) breakup frequency of 7.8 kHz for z in the range of 2-4 mm only. This suggests that for z > 4 mm, the turbulent dispersion smooths the spray fluctuations out. In both cases, there is a low frequency peak at ≈ 1 and 2 kHz for L 1 and L 2 , respectively. Even though this frequency is proportional to the gas bulk velocity in the two cases, it was not detected in the single phase simulation. In addition, in the case of L 2 , the frequency of the PVC (6037 Hz) from the single phase simulation is not visible in the spray quantities. Figure 21 shows the axial profile of the time averaged global SMD. It corresponds to the field shown in Fig. 19 averaged over the radial axis. Directly downstream of the injector ( z ≈ 2 mm), the liquid phase consist of a mix of ligaments (not detected as droplets), and very fine droplets generated by the peeling-off of the liquid sheet, which leads to a small SMD. Then, some ligaments are detached from the liquid sheet, and are accounted as droplets, increasing the global SMD up to z ≈ 4 mm. Downstream this point, all larges structures are detached from the liquid sheet and start to be fragmented in a cascade process, thus decreasing the SMD. Note that in our study, we dismiss droplets smaller than 2 dx (= 67 μ m) because of a too low resolution. Therefore the absolute value of the SMD is overestimated, but the trend is still valid. It is observed in Fig. 21 that even though the global parameters are constant between L 1 and L 2 , the SMD is lower for L 2 , confirming the superior representativeness of the Weber number to describe a spray in airblast atomization.

Tree of Fragmentation
In the present study, the breakup activity was averaged on a crown around the central axis, so that it can be represented on a map depending on r and z. It is shown in Fig. 22 for L 1 (left) and L 2 (right). The key features of a breakup activity map annular airblast atomization are found for both cases. First, most of the breakup events occurs significantly downstream the nozzle, and not directly at the nozzle exit. There might be some peeling-off the liquid sheet directly where the high-speed air stream meets the liquid surface which would produce very small droplets, but since we neglect droplets smaller than 67 μ m, it is not visible in the breakup map where almost no breakup events are recorded for z < 1 mm. Second, the hollow zone delimited by black dashed line is free of breakup event because the flow prevents the liquid structures to reach this zone. Third, the opening of the jet and the turbulent dispersion spread the breakup events in the radial direction. Also, it is noticeable that even though the flow is swirling and tends to eject liquid structures of large Stokes number towards larger radii, there is still a significant number of breakup events close the z-axis. Finally the flapping of the liquid sheet is visible by trumpet-opening-like envelop of the breakup map (contour of N for z ∈ [1, 3] mm), depicted as grey plain line in Fig. 22.
For z > 3 mm, the opening envelope of the breakup events locus for L 1 is less smooth than for L 2 . This suggests for L 1 a more scarce atomization at larger radii, most presumably due to a lower turbulent dispersion and a lower aerodynamic stress.
Another interesting quantity that can be derived from the tree of fragmentation is the fragmentation spectrum, defined by probability density function (PDF) of the child-tomother droplet diameter ratio q = D child ∕D mother . Breakup is identified by q < 1 while coalescence corresponds to q > 1. In Fig. 23 (left), the whole fragmentation spectrum is shown, i.e. for all values of q. Coalescence is observed, but in negligible proportions compared to atomization. The fragmentation spectrum restricted to atomization is shown in Fig. 23 (right). The spectra for L 1 and L 2 have a similar linear trend for q ∈ [0.002, 0.4] (zone delimited by dashed lines), which corresponds to a fragmentation spectrum in power law, as also suggested by Brown (1989). For q ∈ [0.5, 1], the spectrum of L 1 slightly deviates from the linear trend, and the one of L 2 strongly decreases. This means that when a droplet breaks up, the diameter of the children droplets is usually more than twice smaller as its initial diameter, which corresponds to typical observations at high Weber number. Since the spectrum of L 2 is above of the one of L 1 in the range 0.002-0.5, it is expected that the proportion of shearing breakup events is larger in L 2 than in L 1 . This trend is corroborated by the fact that the Weber number for L 2 is four times larger than for L 1 .
As highlighted earlier, the parameters M , GLR and M F are the same for L 1 and L 2 whereas there are noticeable differences in their fragmentation spectrum. Once again, this suggests that the Weber number is more appropriate to describe the cascade process that occurs in secondary atomization.

Lagrangian Coherent Structures
The field of t t+Δt (Δt = 28 μ s) for L 1 is depicted in Fig. 24, on the mid-plane and in the nozzle area. The integration time of Δt = 28 μ s was found by an try and error sequence. If Δt was too small it would lead to a noisy field, and if Δt was too large, it leads to decrease the contrast between coherent and converging zone. In the channel (a), the liquid moves much more slowly than the gas, leading to low values of t t+Δt . The gas of the second stage meets the stream of the first stage at the injector exit at (b). In the recirculation outside the injector (c) the flow has a motion similar to a solid rotation. Concerning the liquid sheet, it is deflected towards the outer radius at (d), which leads to a contraction of the gaseous channel and a deflection of the flow. This deflection leads to a local increase of pressure and to the breakup of the liquid sheet (e). The flow is further deflected in the cavity (f). The remaining tip of the liquid sheet is pushed towards the central axis (g) where it partially blocks the first stage air stream. The flow section of the gas channel increases again (h) whereas it is the contrary on the other side (i). The liquid tip is now fragmented by the first stage air stream (j) in a mode similar to a bag breakup. Note that the central swirling zone is deflected toward the bottom of the picture (k). In addition, it is observed here that the breakup mechanism and the flapping of the liquid sheet on the two side of the injector are not purely in antiphase.
The same FTLE ( Δt = 28 μ s) are shown for L 2 case in Fig. 25. The same features inside the nozzle are observed, whereas a peak of FTLE is visible in the recirculation zone (a), highlighting a zone of large convergence, probably due to the vortex shedding. The same flapping as in L 1 is observed, with the difference that in this case, the sheet is not large enough to firmly obstruct the gas channel. Hence, the sheet is not fragmented in a bag breakup mode, but closer to a shearing breakup where the liquid is peeled-off. Consequently, when the tip is deflected towards the central axis, it does not block the air stream of the first stage, which leads to another peeling-off of the tip (c), and the tip is once again deflected towards the second stage (d). Note that Fig. 25 depicts a complete cycle of the breakup for the top part of the lip while no such event occurs at the bottom lip. This illustrates the random location of the breakup event. To conclude this part, it was shown that the FTLE are of interests to observe coherent flow features inside the nozzle, and to understand the coupling between the air flow field and the dynamics of the liquid sheet.

Conclusion
In this work a generic airblast atomizer operating high pressure with a highly viscous fluid was simulated by means of the SPH method. The geometry was directly converted from a CAD model to a particle file with the in-house tool CAD2SPH. The simulation is conducted in a HPC framework and the tree of fragmentation as well as the FTLE were efficiently extracted thanks to Lagrangian aspect of SPH. From an engineering point of view, it was shown that SPH can extract the same post-processed fields as traditional mesh-based methods, and can even inherently highlight additional flow features which are difficult to obtain for mesh-based methods. The momentum transfer from the gas to the liquid was quantified, the coupling between the airflow and the liquid sheet in annular airblast atomization was revealed by the FTLE fields, and the fragmentation spectrum was extracted for the first time in such a configuration. From an academic point of view, it was observed 1 3 that there is no best parameter to characterize airblast atomization, and the superiority of a non-dimensional parameter group over the other, namely the Weber numbers versus the momentum flux/flow rate ratios depends on the considered field, as summarized in Table 4.