Direct Numerical Simulations of the Evaporation of Dilute Sprays in Turbulent Swirling Jets

The effects of swirled inflows on the evaporation of dilute acetone droplets dispersed in turbulent jets are investigated by means of direct numerical simulation. The numerical framework is based on a hybrid Eulerian–Lagrangian approach and the point-droplet approximation. Phenomenological and statistical analyses of both phases are presented. An enhancement of the droplet vaporization rate with increasing swirl velocities is observed and discussed. The key physical drivers of this augmented evaporation, namely dry air entrainment and swirl-induced centrifugal forces acting on the droplets, are isolated with the aid of additional simulations in which the inertial properties of the droplets are neglected. The correlation between swirl and dry air entrainment rate is found to be responsible for the increase of the global evaporation rate and the spray penetration length reduction, while swirl-induced centrifugal forces are found to be effective only in the jet shear layer, close to the injection orifice, for the analyzed cases.


Introduction
Turbulent sprays are complex multiphase flows playing an essential role in several technological devices as well as in natural processes. This class of flows involves unsteady and multiscale processes, such as turbulent flow motions and dispersed phase transition processes. The main challenges in numerically tackling this problem arise from the mutual interaction of two distinct phases, consisting of mass, momentum, and energy exchange (Elghobashi 2019). Although a satisfactory comprehension of the turbulent spray dynamics has not yet been achieved, the research progress in this field is a key enabling technology for several industrial applications towards the increase of efficiency, emissions reduction, and control. As an example, in combustion engines for aeronautical applications, liquid fuel is directly injected into the combustion chamber, 1 3 where the vaporization of fuel droplets occurs together with chemical reactions within the turbulent, gaseous environment. The formation of pollutants in turbulent flows is related to multiscale phenomena that involve fluctuations of temperature and reactants concentrations (Attili et al. 2014).
A physical description of the liquid-gaseous phase interaction, driving the evolution of bulk liquid jets into sprays, can be found in the work of Lefebvre (2017), while a phenomenological description of the spray dynamics and the corresponding numerical description techniques are reported in the review of Jenny et al. (2012) The atomization dynamics are controlled by the competition between the liquid surface tension, and the gaseous stresses acting on the two-phase interface. Firstly, the liquid jet interacting with the gaseous environment experiences the so-called primary atomization. In this phase, interface instabilities, such as Kelvin-Helmholtz and Rayleigh-Taylor, fragmentize the jet into large drops and liquid ligaments (Lefebvre 2017). The aerodynamic forces acting on the liquid phase, arising from the relative inter-phase velocities, promote additional instability mechanisms that result in a further disintegration of drops and ligaments into smaller drops, namely, secondary breakup. The spray is now referred to as in dilute regime (Chen et al. 2006;Jenny et al. 2012). The atomization process terminates when the surface tension prevails on aerodynamic stresses preventing further fragmentation. In the dilute regime, droplet mutual interactions, such as collisions and coalescence, are negligible. Nonetheless, the effect of dispersed phase on the gaseous phase cannot be discarded, as demonstrated via direct numerical simulation (DNS) by Ferrante and Elghobashi (2003) and Gualtieri et al. (2015), among the others. In contrast with the dense regime, in which the surface-to-volume ratio is low enough to make the vaporization rate negligible, most of the liquid phase evaporates in the dilute regime, where the vaporization rate becomes significant. In this regime, the surface tension dominates on the aerodynamic stresses, this resulting in a spherical shape of the small droplets. Moreover, the droplet size is usually comparable, or smaller, than the smallest scales of the turbulent flow, so that the point-droplet approximation can be successfully employed in simulations as demonstrated by Elghobashi (1994). Nonetheless, in dilute conditions, droplets exert a non-negligible effect on the gaseous flow in terms of mass, momentum, and energy balance that has to be accounted for Elghobashi (1994).
For the reasons mentioned above, a hybrid Eulerian-Lagrangian description results to be the appropriate candidate for the mathematical description of droplet-laden flows in dilute conditions. The Eulerian-Lagrangian description (Jemison et al. 2014;Dalla Barba and Picano 2018) ensures the two-way coupling between the two phases by resolving the classical Navier-Stokes equations for the gaseous carrier flow, while using additional sinksource terms to represent the mass, momentum and energy exchanges between the Eulerian phase and Lagrangian point-droplets. The dynamics of the droplets is described through Lagrangian variables. These are evolved simultaneously with the Eulerian fields describing the dynamics of the carrier phase. The aerodynamic forces governing the point-droplet motion and the thermodynamic variables driving the vaporization rate are computed as a function of the continuous Eulerian variables by means of interpolating kernels. In the framework of DNS, the hybrid Eulerian-Lagrangian approach can be achieved through well-established models meant to take into account the mutual effects of the two phases (Mashayek 1998;Miller and Bellan 1999;Mashayek and Pandya 2003;Bukhvostova et al. 2014;Sardina et al. 2015). Nonetheless, when dealing with large eddy simulation (LES) or Reynolds Averaged Navier Stokes (RANS) approaches, in which only the larger scales are resolved, the interaction of the small-scale turbulent motion with the dispersed drops requires proper modeling (Pozorski and Apte 2009;Li et al. 2019). In this context, a-priori DNS studies can help to isolate the main features of the small-scale fluid-drops interaction and guide the formulation of proper closure sub-grid models for both RANS and LES Eulerian-Lagrangian formulations.
Among the small-scale processes that crucially influence the evaporation rate of particle-laden flows, the preferential segregation of the dispersed phase plays a particularly determinant role, i.e., the inhomogeneous localization of droplets throughout the flow. In a cluster of droplets a rapid increase of the local vapor mass fraction can be observed due to the high local concentration of droplets and, in turn, a reduction of the local evaporation rate takes place (Reveillon and Demoulin 2007;Weiss et al. 2018). A fundamental study on preferential segregation in evaporating turbulent jet sprays was carried out by Dalla Barba and Picano (2018). The study analyzes the DNS simulation of a particle-laden turbulent jet, obtained using a two-way coupling approach between the two phases. Two different mechanisms were found to be the key drivers of droplet preferential segregation. The former is the small-scale inertial clustering, which is a result of the dissipative time-scale of the flow being of the same order of the droplet inertial time-scale. The latter is related to the intermittency of the turbulent/non-turbulent interface dynamics in the jet mixing layer. This phenomenon is due to entrainment of the environmental gas surrounding the jet, which alternately engulfs dry-air regions into high-vapor-concentration structures originated in the core and protracting outwards. This latter mechanism was found to be crucial in the outer part of the jet core, where the evaporation rate peaks, strongly impacting the overall vaporization process.
Besides, practical devices which exploit evaporating sprays are frequently equipped with swirled inflows. The effect of swirling on primary atomization was inquired by Fuster et al. (2009), who numerically described the formation of the hollow cone generated as a consequence of jet rotation. Experimental and numerical investigations on non-Newtonian liquids injected by pressure swirl atomizers were performed by Renze et al. (2011). Kant et al. (2019) examined the sensitivity to mesh resolution of the numerical prediction accuracy of a swirled liquid injection. Gaseous swirled jets originated by axially rotating turbulent pipe flows were inquired experimentally and numerically by Facciolo and Alfredsson (2004); Facciolo et al. (2007), showing that the rotation strongly affects both turbulence and the mean flow properties. Aggarwal and Park (1999) inquired numerically about the effects of swirl and vaporization on droplet dispersion in a two-phase, swirling, heated jet in a low-Reynolds axisymmetric configuration. Wicker and Eaton (2001) studied the effect of swirl on the motion of solid particles dispersed in a recirculating coaxial free jet, focusing on the particle concentration field. Noh et al. (2018) employed LES to assess the performance of various evaporation models in a turbulent droplet-laden flame obtained injecting n-heptane liquid fuel by means of a pressure-swirl atomizer.
Although different studies have examined the dynamics of swirling jets laden with solid/ liquid particles, archival literature lacks in works focusing on the detailed investigation of the evaporation process and preferential segregation of liquid droplets dispersed in swirling jet sprays, at least with the richness of details achievable in a 3D, DNS framework.
The present work aims at extending the study on the influence of small-scale processes on evaporation of Dalla Barba and Picano (2018) by investigating the sensitivity of the jet spray evaporation rate to a swirled inflow velocity. The presence of a swirled motion introduces additional processes that have effects on the particles evaporation rate, which is found to be enhanced with respect to a purely axial jet. Specifically, two concurring phenomena are found to be responsible for a shortening of the penetration length: (1) a centrifugal force acting on the droplets and driving them radially towards the low-saturation shear layer of the jet proportionally to their inertia, and (2) increased dry-air entrainment. A raising of the overall jet spread and the enhancement of the mass entrainment rate have 1 3 been reported for swirling jets both experimentally (Örlü and Alfredsson 2008;Panda and McLaughlin 1994;Park and Shin 1993) and numerically (McIlwain and Pollard 2002). Nonetheless, the effects of the enhanced entrainment rate and the radial forces induced by the swirled motion on the evaporation dynamics of dispersed liquid droplets need to be further investigated. This work aims at characterizing and quantifying the relative importance of these two processes on the evaporation rate and thus provide detailed physical insights, which are expected to be relevant for the development of RANS and LES sub-grid models to be employed in numerical simulations of power plants and aeronautical engines combustion chambers, in which a swirl motion is given to the reactants inflow.

Theoretical and Numerical Formulation
The results presented hereby are obtained through the DNS solution of evaporating sprays in a hybrid Eulerian-Lagrangian framework. The Eulerian gaseous phase is numerically described by means of a low-Mach number asymptotic formulation of the Navier-Stokes equations in an open environment. This approach allows us to account for significant density variations in the gaseous phase discarding the effects of acoustics (Majda and Sethian 1985). Consistently with previous studies in this field (Miller and Bellan 1999;Bukhvostova et al. 2018;Mashayek 1998), the effects of the dispersed phase on the gaseous carrier phase are accounted for through three sink-source coupling terms in the right-hand side of the mass, momentum and energy equations. The point-droplets are treated as rigid evaporating spheres, and the liquid phase properties, such as temperature, are assumed to be uniform inside each droplet. Consistently with the dilute spray conditions, droplet mutual interactions, such as collisions and coalescence, are neglected (Dalla Barba and Picano 2018). Moreover, the gravitational forces are intentionally neglected to inquire about the effects of swirl on the turbulent spray dynamics independently of the physical spatial orientation of the jet. It should be noted that in applications the droplet settling is usually negligible because of the high jet speed. The details of the computational framework, including the resolution of the computational mesh, are provided in the work of Dalla Barba and Picano (2018) and are recalled here for the sake of self-consistency of the document.
The Navier-Stokes equations employed for the present study are written as: where , , and T are the velocity, density, and temperature of the carrier mixture, respectively, while Y v = v ∕ is the vapor mass fraction field, v being the vapor partial density. The viscous stress tensor is = (∇ + ∇ T ) − (2∕3) (∇ ⋅ ) , with the dynamic viscosity of the carrier mixture. Assuming calorically perfect chemical species and a reference temperature T 0 = 0 K, we denote L 0 v as the latent heat of vaporization of the liquid phase evaluated at the reference temperature T 0 . The thermodynamic and hydrodynamic pressure fields are denoted as p 0 and P, respectively. It is worth to remind that, in a low Mach number frame, the zero-order pressure, p 0 , is uniform over the computational domain and constant over time due to the free convection boundary conditions adopted in the present study. The thermal conductivity of the carrier mixture and the binary mass diffusion coefficient of the vapor are k and D . The parameter = c p ∕c v is the specific heat ratio of the carrier mixture, c p and c v being its constant-pressure and constant-volume specific heat capacities, respectively. The carrier phase is assumed to be governed by the ideal gas law 5, where R m =R∕W m is the specific gas constant of the mixture, W m is its molar mass and R the universal gas constant.
Equation (4) is derived by combining the continuity Eq. (1), the equation of state (5) and the energy equation. Starting from the continuity equation, the velocity divergence reads since the equation of state (5) is a functional dependence = (T) , is the thermal expansion coefficient, and for the ideal gas case = 1∕T = ∕p 0 . Then, substituting the energy equation in thermal form: in the Eq. (6) of velocity divergence, we obtain Eq. (4).
The forcing of the dispersed phase on the carrier one is accounted for through sinksource coupling terms in the right-hand sides of the mass, momentum and energy equations, S m , p and S e , respectively. These latter are provided below in discrete notation: where the Lagrangian variables d,i , d,i , m d,i , and T d,i , represent the position, velocity, mass and temperature of the i th droplet, respectively, while c l is the specific heat of the liquid, dispersed phase. The summations are taken over the whole population of droplets located in the computational domain while the discrete delta function, ( − d,i ) , accounts for the fact that each sink-source term acts only at the domain locations occupied by each point-droplet. In this frame, the kinematics, the dynamics, and the thermodynamics of the dispersed phase are completely described by the following Lagrangian equations: where c p,g is the constant-pressure specific heat capacity of the gas, and L v is the latent heat of vaporization of the liquid acetone evaluated at the droplet temperature. The coefficient d = 2 l r 2 d ∕(9 ) is the droplet relaxation time, l being the liquid density. The Schiller-Naumann correlation is used in Eq. 12 to account for the effect of the finite Reynolds number of the droplet, Re d = 2 || − d ||r d ∕ , on the drag force. The Schmidt number, Sc = ∕( D) , and the Prandtl number Pr = ∕(c p k) appearing in Eqs. (13) and (14) account for the mass diffusivity and the thermal conductivity at the droplet surfaces. The Nusselt number, Nu 0 , and the Sherwood number, Sh 0 , are computed as a function of the Reynolds number of the droplets according to the Frössling correlation, The effects of the Stefan flow (Abramzon and Sirignano 1989) are also accounted for by the following corrections applied to Nu 0 and Sh 0 , where H m = ln(1 + B m ) and H t = ln(1 + B t ) with B m and B t the Spalding mass transfer number and heat transfer number, respectively. These latter read: where c p,v is the constant-pressure specific heat capacity of the vapor, Y v is the Eulerian vapor mass fraction field evaluated at the droplet centroid, while Y v,s is the vapor mass fraction in a fully saturated vapor-gas mixture evaluated at the droplet temperature and ambient pressure. This latter is computed by using the Clausius-Clapeyron relation, with v,s the vapor molar fraction evaluated at the droplet surface.The parameters p ref and T ref are arbitrary reference pressure and temperature while R v =R∕W l is the specific gas constant of the vapor. Then, the saturated vapor mass fraction reads: where W g and W l are the molar mass of the gas and liquid phases respectively. The Eulerian Eqs.
(1)-(5) are discretized on a cylindrical staggererd mesh by means of second-order, central finite differences schemes. The equations are integrated in time adopting a low-storage, third-order, Runge-Kutta scheme, see e.g. (Battista et al. 2014;Rocco et al. 2015; Dalla Barba and Picano 2018) for more details. To avoid unphysical oscillations for the mass fraction, Y, in Eq. (2), the convective term of the scalar quantities is discretized using a bounded central difference scheme designed to avoid spurious oscillation, as detailed in Waterson and Deconinck (2007). Moreover, in the limit of dilute systems, the variations of volume fraction in the dynamical Eulerian equations are neglected as in Mashayek (1998) and Miller and Bellan (1999).
It is worth recalling that, in a low-Mach number framework, the divergence constraint provided by Eq. (4) can be prescribed similarly to the way Chorin projection methods impose the divergence-free condition to the velocity field in incompressible flows. The overall procedure for the solution of Eqs. (1)-(5) requires first the integration of the mass conservation equation to compute the fluid density at the time level n + 1 , n+1 . Equation (3) is then solved without the pressure gradient term, ∇P , to compute the unprojected velocity field, n+1 . The latter does not satisfy Eq. (4), but can be corrected via the gradient of a potential field, , defined such that, in semi-discrete notation: By applying the divergence operator, the preceding equation can be recast as: the latter being an elliptical Poisson equation where all the terms are known with the exception of n+1 and n+1 , while the term ∇ ⋅ n+1 are provided by Eq. (4).
The Algorithm 1 provides the iterative procedure used to solve for n+1 and n+1 starting from the estimates ∇ 2 s=0 and s=0 . In the algorithm, the superscript n refers to the temporal levels whilst the s one to the iteration-levels. First, the velocity field, n+1 is projected enforcing for s the local value of the divergence, ∇ ⋅ n+1 , computed according to Eq. (4). Then, ∇ 2 s is re-computed using the updated velocity field, s . Finally, the Poisson equation is solved for the updated value of s by means of a Fast Fourier Transform (FFT) algorithm in the azimuthal direction together with a LU decomposition method for each axial-radial plane. The parameter is used as a convergence indicator; the summation is taken over the entire computational domain. The number of iterations required for the convergence of the algorithm is proportional to the intensity of the density gradients in the flow. In the present case, only a couple of iterations are sufficient to solve for n+1 and n+1 with a residual smaller than t ≃ 10 −8 . Additional information is provided in the references Battista et al. (2014) and Rocco et al. (2015).
The droplet mass, momentum, and temperature laws are evolved by means of a Lagrangian approach. The temporal integration uses the same Runge-Kutta scheme adopted by the Eulerian algorithm. The Eulerian quantities at the droplet positions are computed resorting to a second-order accurate polynomial interpolation.

Test Cases Description
All the presented DNS computations reproduce a turbulent air-acetone vapor jet laden with liquid acetone droplets in the dilute regime. The flow conditions at the inlet section of the non-swirled jet flow, from now on referred to as non-swirled baseline case, are comparable to those adopted in the well-controlled experiments on dilute coaxial sprays published by the group of Chen et al. (2006) and Villermaux et al. (2017), with the only exception of a lower Reynolds number. The experiments use acetone droplets dispersed in air at the temperature of 275.15 K in both non-reactive and reactive conditions.
In the present numerical study, the gas-vapor mixture is injected into an open environment through an orifice of radius R = 5 ⋅ 10 −3 m at a bulk axial velocity U z,0 = 8.1 m/s. The distribution of droplets on the inflow section consists of a random position and velocity distribution of monodisperse, liquid-acetone droplets with an initial radius r d,0 = 6 m. The ambient pressure is set to p 0 = 101,300 Pa, the injection temperature is fixed to T 0 = 275.15 K for both the carrier and the dispersed phases. The injection flow rate of the gaseous phase is kept constant fixing a bulk Reynolds number Re = 2U z,0 R∕ = 6000 , = 1.35 ⋅ 10 −5 m 2 /s being the kinematic viscosity. A nearly-saturated condition is prescribed for the air-acetone vapor mixture at the inflow section, S = Y v ∕Y v,s = 0.99 , where S, is the saturation, Y v is the actual vapor mass fraction on the inflow section and Y v,s (p 0 , T 0 ) is the vapor mass fraction in fully-saturated condition evaluated at the inflow temperature and pressure. The acetone to air mass-flow-rate ratio is set to =ṁ act ∕ṁ air = 0.28 , ṁ act =ṁ act,l +ṁ act,v being the sum of the liquid, ṁ act,l and gaseous, ṁ act,l , acetone mass-flow-rates. This configuration corresponds to a bulk volume fraction of the liquid phase = 8 ⋅ 10 −5 .
The non-swirled baseline case (Dalla Barba and Picano 2018) is characterized by a purely axial inflow velocity. The thermodynamic and physical properties of the vapor gas and liquid phases are summarized in Table 1. The inflow velocity conditions are obtained by assigning a fully turbulent velocity at the jet inflow section via a Dirichlet condition. The two-dimensional inflow velocity field is generated by means of a companion DNS reproducing a fully developed, turbulent pipe flow. The prescribed inflow condition is extracted from a cross-sectional slice of the pipe. The flow is injected through a center orifice in the lower domain base, the remaining part of it being impermeable and adiabatic. In addition to the non-swirled baseline case, two additional configurations are investigated. The simulation parameters for these latter test cases are identical to the ones used in the baseline configuration, with the only exception of the swirl number, defined as S w = U t,0 ∕U z,0 . The velocity U z,0 is the mean, bulk axial velocity computed from the pipe axial flow rate, whilst U t,0 is the mean azimuthal velocity related to the mean azimuthal flow rate of the rotating pipe (Facciolo and Alfredsson 2004), The integrands, ⟨u t ⟩(r) , and ⟨u z ⟩(r) , are the mean axial and azimuthal velocities in the pipe averaged concurrently over time and along the z and directions. The rotating pipe technique was chosen in accordance with the experimental studies reported in Facciolo et al. (2007). As reported in Table 2, the swirl numbers of the Table 1 Thermodynamic and physical properties of acetone and dry air employed in the numerical simulations The pressure p 0 is constant over both time and space due to the low Mach number formulation and open environmental conditions. The temperature, T 0 , the bulk velocity, U z,0 , and dynamic and kinematic viscosities, and , are evaluated at the inflow section. The specific heat capacity at constant pressure of the gas and vapor are c p,g and c p,v , while the liquid specific heat capacity is c l . The molar masses and the thermal conductivities of the gas and the liquid are W g , W l , k g , and k l , respectively. D is the mass binary diffusion coefficient of the vapor into the gas. l and L v are the liquid density and latent heat of vaporization evaluated at the inflow temperature, respectively. The Reynolds number, the orifice radius, and the droplets initial radius are Re, R, and r d,0 , while =ṁ act ∕ṁ air and are the acetone to air mass-flow-rate ratio and the bulk volume fraction of the liquid phase respectively 6 ⋅ 10 −6 Re = U z,0 R∕ 3000 three simulations are set to 0, 0.4, and 0.95 respectively. The three cases will be referred to as Zero Swirl (ZS), Medium Swirl (MS) and High Swirl (HS), respectively. The rotating turbulent periodic pipe employed to generate the inflow boundary condition extends for 2 × 1R × 8R in the azimuthal, , radial, r and axial, z, directions. The domain is discretized with a staggered mesh containing N × N r × N z = 128 × 80 × 128 nodes in order to match the corresponding jet computational grid at the pipe discharge. In Fig. 1 is shown a snapshot of the cylindrical domain for the simulation Sw = 0 and Sw = 0.95 (Fig. 1, left and center panels), along with the turbulent periodic pipe for the case at Sw=0.95 (Fig. 1, right panel). The computational domain consists of a cylinder extending for 2 × 22R × 70R in the azimuthal, , radial, r, and axial, z, directions. The domain is discretized using N × N r × N z = 128 × 225 × 640 nodes distributed on a staggered mesh. The flow is injected at the center of one base of the cylindrical domain and streams out towards the other base. A fixed velocity condition is imposed at the pipe walls. A convective condition is adopted at the outlet, while an adiabatic traction-free condition is prescribed on the side surface of the cylindrical domain to mimic an open environment and make the entrainment of external fluid possible.
The computed ratios ∕ , with = 3 √ r r z , are lower than 4.7 close to the pipe inflow for all the simulations. The values decrease below 3 for the main part of the downstream evolution of the flow, and below 2 in the spray far field. In order to validate the dataset, we report in Fig. 2 the axial distributions of the non-dimensional, mean axial velocity reciprocal, ⟨U z,0 ⟩∕⟨U z ⟩ , computed on the centerline of the domain. The decay of the latter, in the far-field of a round jet in zero pressure gradient conditions, follows the well-known trend ⟨U z,0 ⟩∕⟨U z ⟩ ≃ 0.5(z∕R − z 0 ∕R)∕B , where B is a decay constant and z 0 is the virtual origin of the jet (Hussein et al. 1994). The mean axial velocity reciprocal shows the expected linear behaviour in the far-field, confirming that the boundary   (2014) and Xu and Antonia (2002). The non-dimensional mean axial velocity ⟨U z,0 ⟩∕⟨U z ⟩ multiplied by non-dimensional axial coordinate z/R is reported in the inset of the right-side panel of Fig. 2. At a sufficiently large distance from the inflow, the quantity (⟨U z,0 ⟩∕⟨U z ⟩)(z∕R) exhibits the expected almost constant behavior for all the test cases, confirming that the solution is not affected by mesh resolution and boundary conditions.

Effects of Swirl on the Jet Topology
The effect of the swirl intensity on the turbulent evaporating jet is here discussed. Different statistical quantities are provided. Each mean quantity should be intended as a mean taken over both the tangential direction and time on the Eulerian grid. The Lagrangian mean quantities are computed by assigning the Lagrangian observable to be averaged to the related grid cell on the Eulerian grid. The averages are computed after statistically steady conditions have been attained. Figure 3 shows the contour lines of the mean liquid mass fraction, M , defined on the computational grid as M = m l ∕m g , where m l and m g are the mean mass of liquid acetone and air inside each mesh cell, respectively. An accurate inspection of Fig. 3 reveals that the shape of the region where evaporation occurs becomes both shorter and broader for higher swirl numbers. This effect is particularly pronounced in the first axial region, for z∕R < 10 , where the iso-contours of M are bent downwards for higher swirl intensity because of the stronger centrifugal forces acting on the carrier phase. An overall jet vaporization length is defined as the axial distance from the inflow section, where the 99% of the injected liquid mass has transitioned to the vapor phase. The overall vaporization length appears to move Fig. 2 Centerline values of the reciprocal of the mean axial velocity non-dimensionalized by the bulk velocity, ⟨U z,0 ⟩∕⟨U z ⟩ . The inset provides the nondimensional mean axial velocity ⟨U z ⟩∕⟨U z,0 ⟩ multiplied by nondimensional axial coordinate z/R 1 3 upstream for increasing swirl number, as further reported by Fig. 4, where the non-dimensional vaporization length z v ∕R is plotted as a function of the swirl number.
A contour plot of the mean droplet radius, normalised by the droplet injection radius r d,0 = 6 m , is reported in Fig. 5. The plot suggests that, in the proximity of the inflow, because of the swirl, larger drops are radially advected. As expected, a progressive decrease of the mean droplet radius is observed moving along the jet axis due to the evaporation process. Consistently with the findings of previous studies (Dalla Barba and Picano 2018), a turbulent jet spray can be described as a turbulent core decaying along the axial direction and being surrounded and entrained by dry air. In the region where, due to the entrainment, dry air meets the spray mixture, the vapor concentration diminishes, enhancing vaporization, and reducing the radius of the droplets. Since the effect of dry air entrainment plays a crucial role in the overall vaporization process, the inner jet core shows higher saturation The swirled inflow velocity acts on the spatial location of the mixing layer, increasing the jet spread angle. This line of reasoning is confirmed by the mean droplet vaporization rate distribution, displayed in Fig. 6 for the three inflow conditions. The plot shows how the vaporization is enhanced in the shear layer, attaining maximum values close to the inflow orifice. Moreover, the spread angle of the evaporation region increases with the swirl intensity. In all the test cases, the peak value is found to be located in the shear layer, immediately downstream the inflow section, where large droplets enter in direct contact with the dry environmental air. This region is less elongated in the axial direction, and more stretched in the radial direction for higher swirl numbers. Although for all test cases, for a given axial distance, the radius of the droplets diminishes moving away from the jet axis, the spray core extension in the axial direction diminishes with the swirl intensity. This means that the spray core, which is the region where the vaporization process is slowed down by the high vapor concentration, is shortened in the axial direction by the swirl.
The axial distribution of the non-dimensional mean axial velocity is provided in Fig. 7a. The plot shows the strong influence of the swirl on the centerline velocity decay. Although the pipe exit velocity increases with the swirl number (Facciolo et al. 2007), before z = 10R , the trend is completely overturned. The droplet evaporation is enhanced by the increase of the swirl number, as shown by the centerline distribution of the droplet radius normalized with the inlet radius, see Fig. 7b. On the other hand, the swirl increases the mixing of the vapor phase with the external gas, as shown by the vapor mass fraction axial distribution reported in Fig. 7c, where faster decay of the vapor mass fraction is observed. Figure 8 shows the non-dimensional mean axial, radial, and tangential velocity distributions as functions of the radial distance from the jet axis. The quantities are plotted at three different axial locations z∕R = 10 , z∕r = 20 , and z∕R = 30 , for the three different swirl cases. The data support the evidence of a wider and shorter jet with increasing Swirl number: the mean axial velocity profile flattens and widens, the radial mean velocity profile widens, suggesting higher entrainment of dry air, and the tangential velocity appears intensified at every radial and axial location, as expected.
The sensitivity of the droplet dimension to the swirl number is highlighted by the Probability Density Function (PDF) of the droplet radius evaluated at different axial distances from the origin, as reported in Fig. 9 for the cases at Sw=0 and 0.95. Since the inflow condition is a monodisperse suspension, in all the test cases, the PDF at the inflow section is a Dirac delta function centered at r d ∕r d,0 =1. Nonetheless, consistently with experimental observations in turbulent sprays (Chen et al. 2006), even in the case of zero swirl, shown in Fig. 9a, it is observed a radius distribution spanning around one decade after 2 jet radii away from the inlet section. This intense spread in the droplet radius statistical distribution is enhanced by the increase of the swirl number, as shown in Fig. 9b, where the Sw =0.95 case is shown. As an example, the PDF of the droplet distribution at z∕R = 20 of the highswirl test cases (green curve in Fig. 9b), is very similar to the one observed at z∕R = 40 of the zero-swirl test cases (yellow curve in Fig. 9b).
The previous observations confirm how the swirl, modifying the jet shape and the mixing layer spread angle, promotes the evaporation process leading to the higher spread of droplet dimensions in the regions closer to the inflow. Nonetheless, the identification (c) Fig. 7 Axial distributions of the mean axial velocity, mean droplet radius, and mean vapor mass fraction computed on the centerline of the domain. Each quantity is non-dimensional, the reference scales being set to the bulk inflow velocity U z,0 and the initial droplet radius r d,0 of the physical processes underlying this enhancement of the overall evaporation rate requires further analysis. In the conditions of the test cases taken into consideration, the ambient temperature and pressure are kept constant, and the droplets are injected into a nearly saturated flow. Hence, the only local condition that can lead to a non-zero evaporation rate for a single droplet is a decrease in the vapor mass fraction in the carrier mixture it is surrounded by. Since in all the test cases the droplets are injected into nearly saturated air-acetone mixture, a droplet can evaporate because of one, or a combination, of the following two causes: (1) it escapes the original saturated environment, due to inertial effects, ending up in a less saturated region; (2) it experiences a decrease of the surrounding saturation, because of either laminar or turbulent mixing between the carrier mixture and the surrounding dry air.

Effect of Droplet Inertia on Evaporation
The jet axis distribution of the mean droplet Stokes number based on the Kolmogorov dissipative scale, St = d ∕ for the three swirl numbers is shown in Fig. 10. The turbulent dissipation is computed in correspondence of each mesh node in the range 0.5 < r∕R < 1 . The local dissipative time scale is then adopted to estimate the Stokes number of droplets located within the correspondent cell. The average of these values is considered as the mean droplet Stokes number on the jet axis at a given distance from the inflow section, z/R. The droplets Stokes number is higher than one in the region close to the inlet with the maximum decreasing in value and receding toward the inflow for higher swirl numbers, meaning that the swirl motion increases the inertial effects of the droplets in the first region out the pipe. Further downstream for higher Swirl numbers the droplets evaporate earlier, and St behaves accordingly reaching values lower than those closer to the inflow. The effectiveness of the inertial effects on evaporation of the droplets has been inquired through an additional set of three simulations, whose setup is identical to the previous ones, except for the lack of inertial effects on the droplets. The point-droplets are treated as Lagrangian tracers, with a Stokes number St = 0 ( d = 0 ). Under this assumption, Eq. (12) takes the form d = | = d , while the right-hand side term in Eq. (9) is imposed to zero. In this new set of simulations, hereinafter referred to as tracer simulations, the droplet motion is slaved to the local motion of the carrier mixture.
A comparison between the mean distribution of the droplet non-dimensional radius in the baseline simulations with those obtained for St = 0 is reported in Fig. 11 for the three swirl numbers. The role of the droplet inertial effects on the evaporation is almost negligible in the case of the zero swirl but becomes prominent as the Sw number increases. In the swirled test cases, where the inertial effects are included, the droplet radii decrease faster with respect to the tracer solutions. Nonetheless, this behavior is observed only in the region enclosed between the jet inflow and z∕R = 30 . Further downstream, the droplet radius distributions in the tracers and in the baseline simulations are similar. The regions in which the droplets are smaller than half of the initial radius are almost identical in the simulation sets, also for the highest swirl number. The confinement close to the inflow orifice of the inertial effects is expected, and due to the evaporation process itself. Indeed, the progressive decrease of the droplets' mass, along with an increase of the turbulent time-scales in the downstream evolution of the flow, leads to a substantial diminution of the droplet Stokes number in the downstream regions (Dalla Barba and Picano 2018).
To rigorously quantify the effects of the swirled inflow on the droplet dynamics, a Swirl Stokes number is introduced as: is the droplet relaxation time, while sw = x d,r ∕U t is a rotational velocity time scale, being U t the tangential velocity of the carrier phase at the droplet location and x d,r the radial position of the droplet with respect to the jet axis. The mean values of the swirl-based Stokes number along with the non-dimensional droplet radius distribution are reported in Fig. 12.
The maximum values of the swirl-based Stokes number St sw,max are located where both the mean angular velocity and the mean droplet radius peak, i.e., close to the injector pipe walls, where x d,r = R , and the droplet diameter is maximum, i.e., d d = d d,max . Additionally, St sw,max can be written as a function of the jet swirl number Sw = U t ∕U z,0 , and jet Reynolds number Re = U z R∕ , as well as the Swirl number definition Sw = U t ∕U z,0 , as follows: In the above equation, the direct proportionality of St sw,max to the droplet dimension, the Reynolds number, and the swirl number Sw is outlined. Moreover, the order of magnitude of the swirl number can be estimated as St sw,max ≈ 0.768Sw , assuming that l ∕ ≈ 8 ⋅ 10 2 and the maximum nondimensional droplet diameter d d,max = 2.4 ⋅ 10 −3 . The previous estimates are confirmed, since the ensemble average of St sw is found to be uniformly zero in the non-swirled simulation, Sw = 0 and having peak values of 0.2 and 0.5 in the medium, Sw = 0.4 and high swirl, Sw = 0.95 simulations, respectively.
In the regions where the St sw peaks, the effect of centrifugal forces on the droplets is expected to result in an acceleration in the radial direction, with the effect of observing droplets with a radial velocity being higher than the one of the carrier phase. It is possible to define the ensemble average of the non-dimensional, radial velocity increment u = ⟨u d,r − u r ⟩ as the average of the radial velocity difference between the droplets and the carrier fluid. This latter is provided in Fig. 13. The correspondence between the St sw and the u distribution is evident, with stronger velocity increments in the higher swirl cases.  An inertial effect is present and quantified in terms of a swirl-based Stokes number and a radial velocity increment of the droplets. Nonetheless, it is also evident that these effects are present only in the proximity of the velocity shear layer, where the rotational velocity is the highest, and only in the regions close to the injection pipe, where the droplets are larger and heavier. In the presented test cases, the consequences of the inertial effects are significant until the mean droplet radius is higher than half of the initial one. The differences in the distribution of the mean droplet radius become negligible when the swirl motion is absent, this implying that the inertial effects are related to rotational velocity. Thus, the presence of regions where the inertial effects are relevant while St sw and u are negligible, must be a consequence of the inertial effects occurring in the region close to the injection pipe, where St sw and u are significant.

Effect of Dry Air Entrainment on Evaporation
Observing the distribution of the non-dimensional, mean droplet radius provided in Fig. 11, it is clearly visible how moving away from the inflow section, the differences between the baseline and the tracer simulations vanish. Drawing from the conclusions above and with reference to Fig. 11, it can be inferred that the swirl motions affect the global characteristics of the jet mainly because of a decreasing of the mean vapor mass fraction due to an increasing effect of the dry air entrainment. We introduce a measure of the entrainment as the axially averaged mass-flow-rate occurring in the radial direction, which reads where ⟨ u r ⟩� r=R ∞ is the mean radial momentum of the carrier phase at fixed large r = R ∞ . The quantity ⟨Q r ⟩ formally depends on both z and R ∞ , however it becomes essentially constant for sufficiently high values of the two quantities. This measure, non-dimensionalized by the radial mass-flow-rate in the zero swirl case, is reported in Fig. 14. The entrained flow rate increases with the swirl number, with a ∼50% increase in the High Swirl case. The increase of the entrainment with the swirl number is also consistent with the observation of Fig. 8 panels b, e, and h, in which it is shown that, at a sufficient radial distance from the jet, an increase of the swirl number corresponds to an increase in the radial velocity magnitude. The enhancement of the entrainment intensity is also linked to a faster decrease of the axial velocity, as shown in Fig. 7. The simultaneous observation of Figs. 11 and 14 reveals that the effect of swirl-induced small-scale inertia on the evaporation rate is negligible in the whole field, exception made for the first part of the jet and that the process responsible for the globally enhanced evaporation is indeed the increased dry air entrainment. This physical insight is expected to be a valuable piece of information when dealing with the dispersion and evaporation closure models in RANS and LES context. The role of the entrainment and droplet inertia on evaporation gives also a clear indication for the design of burners, suggesting that swirl motion has an impact on droplets evaporation only if the swirl-induced inertial effects are significant, i.e., large droplets at high swirl number, while the entrainment has a beneficial effect on the evaporation efficiency in any case even for small tracer-like droplets.

Conclusions
The effects of swirled inflows on turbulent jets laden with acetone droplets are investigated by means of direct numerical simulations. The numerical method is based on the low-Mach number expansion of the Navier-Stokes equations for the Eulerian carrier phase, coupled with a Lagrangian description of the dispersed phase based on the point-droplet model. This approach allows taking into account for density variations in the carrier phase as well as for a two-way-coupling between the two phases. The present numerical study focuses on the dilute regime of the overall downstream evolution of jet sprays, where the major part of the liquid phase evaporates. The study provides the outcomes of different DNSs reproducing a turbulent, air-acetone mixture laden with liquid acetone droplets injected into an open environment in nearly-saturated conditions. Two cases with different swirl intensity are considered together with a zero-swirl reference case. It is found that a swirled inflow leads to an enhancement of the overall evaporation rate of the dispersed phase leading to a stronger jet spread and shorter penetration length compared to the reference non-swirled inflow case. This effect is found to be related to two main mechanisms. The former is the droplet inertia: in presence of a swirled inflow, droplets with larger inertia are advected towards the low-saturation mixing layer of the jet, thus increasing their evaporation rate. The latter consists in the enhancement of the mass-flow-rate of the entrained dry air caused by the swirled inflow. The influence of the swirled inflow and the related centrifugal forces on the motion and evaporation of the droplets is inquired employing an additional set of DNSs in which the droplet dynamics are reduced to that of passive Lagrangian tracers. The inertial effects have been quantified by means of a purposely defined Swirl Stokes number. The study highlights as inertial effects are significant only in regions where the Swirl Stokes number is order one or higher. In the present cases, it happens near the injection orifice while they becomes negligible further downstream, where the droplet radius is smaller than half their initial radius at the inflow section. These findings clearly depend on the swirl intensity and droplet sizes. Hence in present conditions, the correlation between the dry air entrainment rate and the swirl intensity is found to be responsible for the global evaporation enhancement and the reduction spray penetration length. An increment of the swirl intensity causes an increase of the jet spreading angle with a consequent decrease of the overall vaporization length. The statistical analysis confirms an enhancement in the overall evaporation rate due to the swirled inflow. The mechanisms discussed in the present study are expected to be of crucial importance for the modeling of turbulent flows characterized by a mixing layer with entrainment of dry air, such as fuel injectors in combustion chambers. The proper modeling of these mechanisms is pivotal for the improvement of LES and RANS model capabilities and to accurately reproduce the turbulent vaporization dynamics for both reacting and non-reacting sprays.