Droplet-laden, turbulent sprays are complex multiphase flows where a dispersed liquid phase is transported by a turbulent gaseous one. The distinguished phases mutually interact exchanging mass, momentum and energy in the turbulent regime. These complex flows play an important role in many industrial applications as well as in a large variety of natural and environmental processes. Despite the technical and scientific relevance of the problem, the understanding and modeling of the fundamental mechanisms regulating the evaporation and dispersion of droplets in turbulent sprays are still a challenging matter due to the unsteady, multi-scale and multiphase processes involved in these complex flows [23] . In industry, for instance, an improved understanding of such processes would contribute significantly to the development of high-efficiency and low-emission internal combustion engines. In these applications, the liquid fuel is directly injected into a hot-temperature and high-pressure combustion chamber where the liquid fuel undergoes a phase transition, via primary and secondary atomization, followed by a fierce chemical reaction process within the turbulent gaseous environment. Such a complex process completes in a very short time inside the combustion chamber that multi-scale, multiphase and unsteady phenomena are involved. A typical example is the formation of pollutants which is related to the fluctuations of temperature and reactant concentrations. In particular, soot forms through a pyrolysis process in fuelrich regions where high temperature subsists without enough oxidizer to react [26, 2, 41] . This phenomenology can be observed within droplet clusters, where the concentration of the fuel droplets can be even a thousand times higher than its bulk value giving rise to a peak in the fuel vapor concentration. Another recent instance relies on the transmission of respiratory infections, such as COVID-19, through virusladen droplets that are formed in the respiratory tract of an infected person and expelled from the mouth and nose during spasmodic events such as coughing and sneezing [36, 4] . In this context, to prepare ourselves to tackle the outbreak of such respiratory diseases, it is mandatory to improve our knowledge and modeling capability concerning the fundamental mechanisms which govern the dispersion range of virus-laden droplets ejected into the environment within the exhaled turbulent vapor-air mixture.
A phenomenological description of the overall evolution of a turbulent spray can be found in the review of Jenny et al. [23] . The whole process starts with the jet breakup, or atomization, as a high-velocity liquid flow is injected from a duct into a gaseous atmosphere. During the so-called primary atomization, Kelvin-Helmholtz and Rayleigh-Taylor instabilities develop at the gas-liquid interface promoting the jet fragmentation and giving rise to a system of large droplets and ligaments enclosed within the environmental gas [32] . Further downstream, secondary atomization follows: the aerodynamic forces, arising as a consequence of the relative inter-phase velocities, causes interface instabilities that lead to a further disintegration of droplets and ligaments into even smaller fragments [32, 23] . The breakup process occurs in a so-called dense regime and terminates when the surface tension prevails on aerodynamic stresses preventing a further fragmentation [23] . At this stage, a dilute regime establishes where the dispersed phase volume fraction is small enough that the droplet mutual interactions, such as collisions and coalescence, can be neglected [17] . Whereas in the dense regime the vaporization rate is negligible, the vaporization process becomes significant in dilute conditions and most of the liquid evaporates at this stage of the turbulent spray evolution [16] . In the dilute regime, small droplets evolve preserving a spherical shape due to the dominance of the surface tension on the aerodynamic stresses [15, 17, 23] . Even though the presence of the liquid phase still exerts a considerable effect on the flow in terms of mass, momentum, and energy balance, in dilute conditions, the droplet size is usually below, or comparable to, the dissipative scales of the turbulent flow. In these conditions the point-droplet approximation is widely accepted and the hybrid Eulerian-Lagrangian approach is considered suitable and reliable for the numerical description of problem [15] . In this framework, the mass, momentum, and energy exchanges between the Eulerian carrier phase and Lagrangian point-droplets can be accounted for via distributed sink-source terms in the right-hand side of the balance equations governing the Eulerian phase while the droplet kinematic and thermodynamic variables are evolved in the Lagrangian frame [15, 33, 8, 14, 12, 49] .
In order to gain an insight into the fundamental physics of turbulent sprays, various theoretical and numerical studies, which complement traditional experimental techniques, have been conducted with a rapid growth of their applications. One of the former attempts to describe the vaporization process of spherical droplets dragged by a gaseous phase flow relies on the well-known d-square law, that has been formerly proposed by Spalding [47] and Godsave [18] . Fixing the thermodynamical properties (temperature, vapor concentration, etc.) in a quiescent environment, they found that the droplet surface decreases linearly with time. Essentially based on the d-square law, Law [30] and Law and Sirignano [29] proposed the rapid-mixing model: the temperature is assumed to be uniform inside the droplet but temporary varying. The model applies well under the assumption that the thermal conductivity of the liquid phase is high enough to permit a rapid stabilization of the droplet internal temperature in relation to local changes of the environmental temperature. Later, Abramzon and Sirignano [1] proposed an improved droplet vaporization model applicable to non-uniform and time-dependent environmental conditions, which takes into account forced convection, molecular diffusion, and the Stefan flow contribution to the vapor transport from the droplet surface to the neighboring environment. The author pointed out that the Ranz-Marshall correlations [42] for the Nusselt and Sherwood numbers, and ℎ, respectively, adopted in the rapid-mixing model and d-square law may lead to a physically incorrect super-sensitivity of the transfer rates to the small turbulent velocity fluctuations in the vicinity of the zero Reynolds number.Considering the motion of a finite size spherical particle, or droplet, in turbulent flows, Maxey [34] derived a remarkable equation which accounts for the Stokes drag, added-mass effect and gravitational force under low-Reynolds number circumstance. All these approaches, among many other reported in archival literature, have lead to the development of well-established and reliable pointdroplet equations in the hybrid Eulerian-Lagrangian frame that have proved to be effective in reproducing droplet evaporation dynamics in both the direct numerical simulation and large eddy simulation frames [15, 35, 6, 8, 49, 9] .
Along with the previously reported models for droplet dynamics in turbulent gaseous flows, various numerical techniques have been proposed and many are currently in use to tackle the simulation of dilute sprays laden with dispersed droplets. The hybrid Eulerian-Lagrangian numerical approaches based on the Reynolds averaged Navier-Stokes are widely used to perform prediction for industrial applications due to their simplicity and low computational requirements. The main limitation of these methods, which is partially removed in the Large Eddy Simulation frame, is incapable to capture the scales of time-dependent and unsteady motion occurring in turbulent sprays [37] . In contrast, the Direct Numerical Simulation (DNS), that resolves all the length and time scales present in a turbulent flow, has proved its capability to capture the whole physics of the spray vaporization process providing an accurate insight into the complex phenomena involved. Mashayek [33] conducted one of the first DNS of evaporating droplets in turbulent flows. The author neglects the effect of the dispersed phase on the carrier one (one-way coupling) under the assumption of incompressible flow. Bukhvostova et al. [8] compared the performances of the incompressible and low-Mach number asymptotic formulation of the Navier-Stokes equation in reproducing the dynamics of a turbulent evaporating spray, considering the effects of changes in the total mass density caused by evaporation of droplets and condensation of water vapor. Although the results don't show significant differences in mild initial conditions, the usage of the low-Mach number formulation is found to be crucial in order to obtain a reliable quantitative prediction of heat and mass transfer. Reveillon and Demoulin [43] studied the impact of preferential segregation of droplets on the mixture fraction field with one-way coupling method. They found that the evaporation process evolves in three different stages in time: single droplet mode in the early stage, cluster mode in the intermediate stage and a gaseous mode in the late stage. Recently, Dalla Barba and Picano [14] investigated the dynamics of droplet vaporization within a three-dimensional turbulent spatial developing jet in dilute, non-reacting conditions, using DNS framework based on a hybrid Eulerian-Lagrangian approach and the point-droplet approximation. The preferential segregation of droplets is found to be driven by two distinct mechanisms: the inertial small-scale clustering in the jet core and the intermittent dynamics of the jet across the turbulent-non-turbulent interface in the mixing layer, where dry air entrainment occurs.
Although different numerical studies, in conjunction with theoretical exploitation, have been already conducted, it is still premature to draw a conclusion that a satisfying comprehension of turbulent spray dynamics has been achieved. The capabilities of existing models of reproducing the involved phenomena are not thoroughly confirmed. Moreover, most part of the numerical dataset and studies available in archival literature in the frame of direct numerical simulation is limited to relatively low Reynolds numbers due to the high computational requirements of such simulations. To cover these lacks and to improve model capabilities for applications, the availability of data from three-dimensional direct numerical simulation of evaporating jet sprays at high Reynolds number is necessary. In this context, the present paper aims to address the numerical study of the basic mechanisms governing the evaporation and preferential segregation of droplets in a three dimensional turbulent jet-spray at high Reynolds number. The problem is addressed via a three-dimensional direct numerical simulation. In particular, the focus is on the effect of the Reynolds number of the jet on different Lagrangian and Eulerian observables, such as evaporation length and droplet size spectrum. A comparison between the data produced by the present simulation and the previous DNS data by Dalla Barba and Picano [14] at = 6, 000 is provided. In both cases a similar and strongly inhomogeneous distribution of droplets is observed. Nonetheless, as the Reynolds number is increased from 6, 000 to 10, 000, the vaporization length grows and the mean evaporation rate is lower. Different mechanisms providing an explanation for these differences are discussed in detail. The authors believe that the data reported in this paper could contribute to improve our knowledge about the dynamics of turbulent evaporating sprays, providing, at the same time, a benchmarking test case for the future low-order modeling of the phenomena involved.
The numerical results reported in this paper have been computed employing a MPI-parallel code that has been previously used in other numerical studies and has undergone an extensive validation and testing campaign [44, 38, 14, 12] . The numerical algorithm is based on a hybrid Eulerian-Lagrangian approach and the point-droplet approximation. The exchanges of mass, momentum and energy between the dispersed phase and the carrier one are accounted for, whereas collisions and coalescence within the dispersed phase are neglected in a two-way coupling frame [17] . A brief description of the numerical methodology is provided in the following for the self-consistency of the paper, the reader being referred to the references for additional details [44, 38, 14, 12] . The conservation equations for the gas phase relies in a low Mach number asymptotic formulation of the Navier-Stokes equations in an open environment:
The Eulerian fields , , and are the density, velocity, temperature and hydrodynamic pressure of the carrier phase, respectively. In the frame of the low-Mach number asymptotic expansion considered in this paper, the thermodynamic pressure, referred to as 0 , results to be spatially uniform and constant in time, due to the open environment conditions Majda and Sethian [31] . The Eulerian field = ∕ is the vapor mass fraction, with the partial density of the vapor. The viscous stress tensor is , , with the dynamic viscosity of the carrier phase. The parameter = ∕ is the specific heat ratio of the carrier mixture, where and are its specific heat capacities at constant pressure and volume, respectively. The flux vectors and represent the thermal and mass diffusion and are computed according to the Fourier's law and Fick's laws, respectively, the parameters and being the thermal conductivity of the carrier mixture and the binary diffusivity of the vapor into the gas. A reference temperature 0 = 0 is fixed and the assumption of calorically perfect chemical species is considered to estimate the sensible enthalpy. Under this hypothesis, 0 is the latent heat of vaporization of the liquid phase evaluated at the reference temperature, 0 . The carrier vapor-gas mixture is assumed to be governed by the state equation for perfect gases (5) , where = ∕ is the specific gas constant of the carrier mixture, its molar mass and  the universal gas constant. Consistently with previous studies in this field [33, 8] , the effect of the dispersed phase on the carrier one is accounted for by employing three sink-source coupling terms in the right-hand side of the mass, momentum, and energy equations, , , ,and , respectively:
The Lagrangian variables , , , , and are the position, velocity, mass and temperature of the droplets, respectively, whereas the parameter is the specific heat of the liquid phase. The summations are taken over the whole droplet population located within the domain (index ). The delta function, ( − ), accounts for the fact that the sinksource terms act only at the Eulerian positions occupied, at a given time, by the point-droplets. In the numerical algorithm, the Eulerian terms (6)- (8) are calculated, in correspondence of each node of the computational grid, by volume averaging the mass, momentum, and energy sources arising from all the droplets located within the grid cell pertaining to the considered grid node.
As mentioned above, the point-droplet approximation is adopted to describe the dispersed phase: the droplets are treated as rigid evaporating spheres and the liquid properties (e.g. temperature) are assumed to be uniform inside each droplet. Since the focus of this paper is on the spray dynamics restricted to dilute conditions, mutual interactions among droplets (e.g. collisions, and coalescence) are neglected. For a discussion about the validity of the preceding assumptions the reader is referred to [14] . Under these hypotheses, the dynamics of the droplets is completely described by the following Lagrangian equations:
where the parameter is the density of the liquid phase, , is the heat capacity at constant pressure of the gaseous component of the carrier mixture and the latent heat of vaporization of the liquid phase evaluated at the droplet temperature. The variable is the droplet relaxation time, = 2 2 ∕(9 ), with = 2 ‖ − , ‖ ∕ is the droplet Reynolds number. The Schiller-Naumann correlation is adopted in equation (10) to account for the effect of the finite Reynolds number of the droplets on the drag. In equations (11) and (12) , the Schmidt number, = ∕( ), and Prandtl number, = ∕( ) are employed to compute the mass diffusivity and thermal conductivity, respectively. Besides, the Nusselt number, 0 , and Sherwood number, ℎ 0 , are estimated as a function of the droplet Reynolds number through the Frössling correlation:
ℎ 0 = 2 + 0.552
A correction is then applied to 0 and ℎ 0 to account for the effect of Stefan flow Abramzon and Sirignano [1] :
where , is the specific heat capacity of the pure vapor at constant pressure, is the vapor mass fraction evaluated at the droplet position, whereas , is the vapor mass fraction evaluated at the droplet surface. The latter corresponds to the mass fraction of the vapor, in a fully saturated vapor-gas mixture, evaluated at the droplet temperature. In order to estimate , we employed the Clausius-Clapeyron relation:
where , is the vapor molar fraction evaluated at the droplet temperature and 0 the thermodynamic pressure. The parameter is the saturated vapor pressure evaluated at the reference temperature , whereas = ∕ is the specific constant of the vapor. The saturated vapor mass fraction is, then,
The Eulerian computational domain is a cylinder. The droplet-laden jet-spray is injected through a circular orifice of radius located at the centre of the lower base of the domain and streams out towards the opposite one. The domain extends for 2 × 20 × 70 in the azimuthal, , radial, and axial, , directions and is discretized by means of a staggered, cylindrical grid of × × = 192 × 211 × 1152 nodes. An Eulerian algorithm advances in time the Eulerian fields by solving the Low-Mach number formulation of the Navier-Stokes equations (1)-(5). Second-order, central finite difference schemes are employed for the spatial discretization, whereas the temporal evolution is performed by a low-storage, third-order Runge-Kutta algorithm. The computational grid is uniform along the azimuthal direction whereas it is stretched along the radial and axial ones. The grid spacing is maintained of the same order of the Kolmogorov length scale for the whole downstream evolution of the jet-spray. Details about grid resolution and adequacy of grid spacing can be found in [14] for bulk Reynolds number, = 0 ∕ = 6, 000. In the present case, the grid spacing has been obtained by a scale-up of the computational grid used in [14] in order to maintain an adequate resolution down to the dissipative length scales, but at a higher Reynolds number, = 10, 000. The mean grid spacing ratios between the present and reference case, Δ 10,000 ∕Δ 6,000 , are about 0.67 and 0.55 along the azimuthal and axial directions, respectively. In these directions, the grid spacing ratio is of the same order of the dissipative length-scale ratio 10,000 ∕ 6,000 ≃ (6, 000∕10, 000) 4∕3 ≃ 0.68. Meanwhile, to get benefits of high-resolution in the shear layer, a stretched and dense mesh instead of a uniform one is implemented along the radial direction near the pipe in which the present mesh spacing is about 0.44 times of the reference one, though the total grid points don't change too much. A convective boundary condition is adopted at the outlet section located on the upper base of the domain. An adiabatic, traction-free boundary condition is prescribed at the side boundary of the domain making the entrainment of external fluid possible, which, in the present case, consists of dry air. Timedependent and fully turbulent boundary conditions are prescribed on the inflow section by means of a companion DNS reproducing a fully-developed, periodic pipe flow. A fully turbulent velocity field is assigned on the jet inflow by a Dirichlet condition. This two-dimensional field is computed on a cross-sectional slice of the turbulent pipe. Excluding the circular inflow, the remaining part of the domain base is impermeable and adiabatic. The turbulent pipe extends for 2 × 1 × 6 in the azimuthal, , radial, and axial, directions. The pipe domain is discretized by a staggered grid containing × × = 192 × 91 × 128 nodes in order to match the corresponding jet-grid nodes at the pipe discharge. A sketch of the cylindrical domain together with the grid structure and the turbulent, periodic pipe is provided in figure 1. Concurrently with the Eulerian solver, a synchronous Lagrangian solver evolves the droplet mass, momentum, and temperature laws (9)-(12) by using the same Runge-Kutta scheme adopted for the Eulerian solver. Numerical stability issues, that can arise when the droplet size becomes too small, are avoided by setting a threshold radius , ℎ corresponding to a mass-loss due to evaporation of 99.95% of the initial droplet mass. When the radius of a droplet decreases under the fixed stability threshold, the droplet is considered as completely evaporated and is removed from the simulation.
The present paper addresses the numerical simulation of liquid acetone droplets dispersed within a gaseous jet consisting of mixture of air and acetone vapor. The jet streams out from the inflow section into an open environment filled by dry air. The absolute pressure of the environment is 0 = 101300 while the temperature is set to 0 = 275.15 . The radius of the inlet section is set to = 4.9 ⋅ 10 −3 while the bulk inflow velocity of the jet is 0 = 13.9 ∕ . A monodisperse population of liquid acetone droplets of initial radius ,0 = 6 is randomly distributed over the inflow section at each time-step. The injection temperature is fixed to 0 = 275. 15 for both the droplets and the carrier mixture. The injection flow rate of the gaseous phase is kept constant by fixing a bulk Reynolds number = 2 0 ∕ = 10, 000, with = 1.35 ⋅ 10 −5 2 ∕ the kinematic viscosity. A nearly saturated condition is prescribed at the inflow for the airacetone vapor mixture, = ∕ , = 0.99, with the actual vapor mass fraction and , ( 0 , 0 ) the vapor mass fraction at saturation, evaluated at the actual inflow temperature and thermodynamic pressure. The acetone mass flow rate is set by means of the ratio Φ =̇ ∕̇ = 0.28, wherė =̇ , +̇ , is the sum of the liquid and vapor acetone mass flow rates, whilė is the dry-air flow rate. The correspondent bulk volume fraction of the liquid phase is Ψ = 8.0 ⋅ 10 −5 . All the thermodynamic and physical properties of the vapor, gas, and liquid phases are reported in table 1. The thermodynamic conditions at the inlet section are comparable to that adopted in the well-controlled experiments on dilute coaxial sprays published by the group of Masri and coworkers Chen et al. [10] . The time step is set to Δ ∕ 0 = 0.001 where the reference time scale is time scales in order to achieve a statistical steady condition for the two-phase evaporating flow before collecting data. All the statistics presented in the following are computed considering around 100 samples separated in time by ∕ 0 = 1.
Concerning the reliability of the simulation, additional information and validation benchmarks about the turbulent periodic pipe, single-phase jet simulations as well as evaporation model can be found in [14, 12] . Table 1 Thermodynamic and physical properties of acetone and dry air. 
In this section the outcomes from the present DNS case, = 10, 000, is shown and compared with results at = 6, 000 [14] . The essential difference between the two cases taken into consideration is the bulk Reynolds number of the injected carrier phase, all of the other parameters being the same. From physical point of view, we are comparing two jet sprays with different bulk velocities. It is worth noting that, the initial Stokes number of droplets injected into the domain is larger in the = 10, 000 case, since
In particular, the value of the initial Stokes numbers are 0 ≃ 1.04 and 0 ≃ 0.62 for the = 10, 000 and = 6, 000 cases, respectively.
A macroscopic overview of the turbulent jets can be given visualizing coherent structures and enstrophy. Several methods have been proposed in order to identify, quantify and visualize the three-dimensional coherent structures in incompressible turbulent flows [22, 11, 24, 51] . We address the visualization of vortical structure by means of the -criterion [22] :
with
where Ω and are the vorticity tensor and the rate-ofstrain tensor, respectively. By this approach, a vortex is defined as a spatial region where > 0, . . the Euclidean (c) Radial-axial slices of the turbulent sprays at = 10, 000 and = 6, 000. The black points represent a subset of the whole droplet population constituted by droplets located within a distance Δ∕ = 0.1 to the slice plane. Each point size is proportional to the corresponding droplet radius with a scale factor 100. The carrier phase is contoured according to the instantaneous vapor mass fraction field, , which is bounded between 0 and 0.18, the former value corresponding to the dry condition while the latter to the 99% saturation level prescribed at inflow section.
norm of the vorticity tensor dominates that of the rate of strain [20] . In general, larger values are associated with more intensive vortexes. Figure 2a shows a comparison of the instantaneous vortical structure of the flow for the two Reynolds number cases visualized by the iso-surfaces = 1 contoured by the magnitude of the Eulerian velocity field. Moving downstream from the inflow section, the vortical structures spread out into the far-field while subjected to a decay process. As the Reynolds number is increased, vortical tubes becomes denser and smaller, but also the vortices spread is prolonged downstream [3] . This aspect suggests a longer turbulent spray jet evolution for the higher Reynolds number. Figure 2b shows a two-dimensional radial-axial snapshot of the instantaneous distribution of droplets and the instantaneous enstrophy field 2 = ‖∇ × ⃖ ⃗ ‖. As expected, in the near-field, high-intensity vorticity regions are mainly located in the shear layer whereas vorticity is lower in the jet core. It should be noted that present jets are generated by fully turbulent pipe flows and so fluctuations are present also in the core. Hence the so-called potential core region should be named as unperturbed core region where the turbulence statistics remain closer to the inlet pipe flow conditions and not to those of jet far-field. Moving downstream, the jet decay leads to a different behavior. This dynamics, that is intrinsically related to the physics of turbulent jets, is observed in both the cases, the high-Reynolds jet presenting vorticity-related features that are shifted downstream compared to the lower Reynolds one. The decay of turbulent structure of the jet observed in this paper in consistent with the experimental observation by [27, 28] . Droplets develop clustering and appear less frequently in high vorticity regions especially in the shear layer. We attribute this effect to two different causes. First, inertial droplets tend to escape from high vorticity regions which is an established explanation of small-scale clustering [48] . Second, near high vorticity regions in the shear layer strong mixing events occur where dry environmental air depleted of droplets is mixed with inner saturated gas laden of droplets. This mixed (non-saturated) regions lead to a fast droplet evaporation that disappear. The combination of these two phenomena originate a strong clustering. This anti-correlation between droplets and vapor phase is better highlighted in figure 2c where the instantaneous vapor mass fraction field together with the instantaneous droplet distribution is provided. It is worth noticing that most part of the droplets are located in the core regions of the turbulent jets in which a strongly inhomogeneous preferential distribution is observable over the whole downstream evolution of the two turbulent sprays. Nonetheless, in the near-field region, this non-uniform characteristic of droplet dispersion emerges mainly in the mixing layer. This is a source of clustering which propagates also in the downstream evolution. We also emphasize that, in the near-field, a self-preserving core exists with a nearly saturated vapor environment which prevents the droplet to evaporate. In both cases, moving downstream the inhomogeneous distribution of the droplets becomes more apparent and can be observed both in the near-axis region and in the jet shear layers. Nonetheless, with the bulk Reynolds number increasing from = 6, 000 to = 10, 000, the onset of this intensive preferential segregation of the droplets dispersion appears to be slightly shifted in the downstream direction. In addition, we also underline how long is the evaporation length for the higher Reynolds number case where several droplet clusters can be appreciated at a higher distance from the inlet. and are the mass of the liquid acetone and the mass of gaseous phase evaluated inside each mesh cell. A spray vaporization length is defined as the axial distance from the inflow section to where the droplets lost, in average, 99% of their initial mass (Dalla Barba and Picano [14] ). According to this definition, the vaporization process ends at about ∕ ≃ 44 and ∕ ≃ 50 for = 6, 000 and = 10, 000, respectively. The vaporization length of the high-Reynolds jet is about 13% longer than the low-Reynolds one. An additional comparison of the spatial distribution of the average liquid mass fraction between the two cases is provided in figure 3b that shows the centerline mean liquid mass fraction versus the axial coordinate, . The centerline behavior exhibits a near-field hump: the liquid mass fraction increases along the jet axis from the inlet section up to a peak located at ∕ ≃ 9 for = 6, 000 and ∕ ≃ 11 for = 10, 000. Further downstream, Φ reduces along the jet axis until droplets inside the turbulent spray totally evaporate. The centerline behavior of Φ keeps higher values in the = 10, 000 case than in its counterpart. This behavior is consistent with the observations by [39, 27, 28] in which the centerline concentration of inertial particles is found to increase above the exit value. The phenomenon can be explained considering the interplay of the droplet/particle inertia and the decay of the mean fluid velocity. The inertia induces a delay of the particle velocity to adapt to the slower (decayed) flow velocity creating a local concentration peak, see Picano et al. [39] for more details.
Fluid and particle mean velocities are displayed in fig-ure 4. In the panels 4a and 4b for each case is shown the mean gas velocity in the half-left image and the mean droplet velocity in the half-right image. As apparent, the dispersed phase velocity is slightly higher than the corresponding fluid velocity, but no significant differences emerge comparing the jets at the two Reynolds numbers. The centerline values of the same quantities, shown in panel 4c, confirm this overall impression. Given that, the macroscopic flow behaviors of present cases is not affected by the Reynolds number, so the difference found in the evaporation lengths is not induced by a deviation of the gas and particle velocity fields.
In a turbulent spray, as the turbulent core spreads and slowly decays, the dry and irrotational environmental air surrounding the jet is entrained continuously diluting the vapor concentration and permitting the overall vaporization process to advance. Since the inner core fluid does not reach the outer jet dry region, but is just diluted by the entrainment, the centerline region presents a higher saturation level over the whole downstream evolution of the flow. Mean saturation field, , is provided in figure 5a. Both flows are almost saturated near the inlet section, as prescribed by the inlet conditions. The saturation level gradually decreases in the downstream direction, maintaining a sharp gradient towards the outer jet region. No significant differences are found between the two cases. Some small discrepancies are present in the intermediate and far-field region, ∕ > 20. The case at = 6, 000 (left panel) shows a relatively higher saturation level in the region 20 < ∕ < 30 while the opposite occur beyond ∕ ≃ 30. As will be further addressed in the following, at lower Reynolds number we observe a faster evaporation rate in the first part of the spray, which explains the higher saturation level in this region. Conversely, in the far-field there are more dispersed evaporating droplets at the higher Reynolds number, so we expect a higher saturation level. Figure 6 provides a comparison between the two cases for the mean droplet radius, ⟨ ⟩. In both cases, in correspondence of each axial position, ∕ , larger droplets are always located in the jet core region whereas smaller droplets can be found towards the mixing layer. This behavior is expected since the saturated jet core reduces the vaporization rate of the droplets near the centerline. On the other hand, the entrainment of environmental dry air dilutes the acetone vapor concentration in the jet mixing layer enhancing the vaporization process.At the higher Reynolds number, we find larger droplets denoting a slower average evaporation rate. A more quantitative view assessing this aspect is provided in panel 6b where the centerline value of the mean droplet radius is shown. The mean evaporation rate is shown in figure 7a where the droplet evaporation rate is normalized by the flow time scale 0 and the initial droplet mass ,0 . Evaporation peaks in the mixing layer and is more intense in the near field. The observed behavior is consistent with the mean droplet radius previously discussed. From quantitative point of view, the case at = 10, 000 shows less intense values which explain the longer evaporation length. Dealing with dimension-less quantities, droplet mass transfer is proportional to 1∕ ∝ 1∕ , implying a relatively slower mass (a) (b) Figure 6 : (a) Average droplet radius, ⟨ ⟩, scaled by the droplet initial radius ,0 = 6 . (b) Mean droplet radius distribution near the center axis with the range 0 < ∕ < 0.2 transport with respect to advection time scale. Nonetheless, the dissipative length-scale, as well as the finer length-scales of turbulent mixing-layer structures, become smaller in the high-Reynolds case, leading to a faster time-scale of the entrainment process arising in the mixing layer. This peculiar feature makes the (non-dimensional) entrainment, spreading and decay rates of turbulent jets independent of the Reynolds number [40] . Being the dispersed droplets inertial, they are not able to get fully benefit of these fast features and globally show a slower evaporation rate. To better characterize this aspect, in figure 7b we show the mean droplet evaporation rate normalized by the initial droplet relaxation time 0 = 0 and the initial droplet mass. This normalization takes into account the droplet inertia and would remove the Reynolds number direct dependence. Hence, we observe that at higher Reynolds number the mean ( 0 -dimensionless) evaporation rate appears a bit higher with respect to the case at = 6, 000. The more intense fluctuations of the higher Re-number case fasten the evaporation rate normalized with the droplet-relaxation time, i.e. accounting for the proper droplet time-scale. However, from physical point of view, this small growth is not sufficient for compensating the faster advection time-scale 0 induced by the higher inflow velocity of = 10, 000. Hence, this explains the longer jet evaporation rate which, nevertheless, does not scale proportionally to the Reynolds number. As previously stated, a strong preferential segregation of the dispersed phase rises, both in the low and high Reynolds number cases. Several mechanisms governing the preferential concentration of droplets in a turbulent flow have been proposed, e.g. the small-scale clustering [48] , the sweepstick mechanism [19] , the accumulation of droplets along jet axis [39, 27, 28] and the intermittent dynamics of the jet mixing layer [14] . Independently from the phenomena giving rise to the observed inhomogeneous and preferentially-segregated spatial distribution of droplets, within clusters the local vapor concentration is significantly higher than its bulk counterpart. Consequently, the evaporation process of droplets located into a cluster is significantly slowed down. The vaporization may be even locally stopped, if the vapor concentration reaches the saturation level, , ( 0 , ), producing a nonevaporating core around the cluster [43] . In literature, several approaches have been proposed to measure the intensity of the preferential segregation of the dispersed phase in droplet or particle laden multiphase flows [46] . Here, in order to provide a spatial map we employ the clustering index [5, 14] :
To compute in equation above, (21) , the Eulerian domain is discretized employing an uniform and equispaced Cartesian grid of cubic sampling cells. The edge size of the cells is set to ∕ = 0.2. The variables and ( ) 2 are the mean and the variance of the number of droplets located into each sampling cell, respectively. The clustering index, , is vanishing for any cell where the distribution of the droplets is random (Poisson process), whereas becomes positive when the variance exceeds the mean value due to the existence of a preferential concentration of droplets within the considered cell. In this sense, large positive values of correspond to a strong segregation of the dispersed phase whereas low values to a uniform spatial distribution of droplets. The results are presented in figure 8a , showing the contour plot of the mean clustering index computed by averaging the instantaneous and local values of along the azimuthal direction and over time. The distribution and the intensity of the preferential segregation of the dispersed phase present only a weak dependence from the bulk Reynolds number of the jet, since no strong differences exist in the spatial distribution of , at least for the cases considered in this paper. In the immediate proximity of the inflow section, assumes positive values mainly in the mixing layer, whereas it is only weakly positive in the core. On the other hand, in the intermediate and farfield, assumes quite large positive values over the whole turbulent jet, both in the core and mixing layer. The peak values for the clustering index occur in the jet core, between ∕ ≃ 10 and ∕ ≃ 20 in both low and high Reynolds cases, the latter presenting a small shift in the downstream direction of the -peaking-region.
To determine the leading mechanism driving the preferential segregation of droplets in different spray regions, we consider the droplet Stokes number evaluated at the dissipative time scale. Figure 8b provides the trend of the mean Stokes number , defined as the ratio of the droplet response time, , and the characteristic time of the dissipative scale, = ( ∕ ) 1∕2 , computed along the jet axis and plotted versus ∕ . The local dissipative time scale is evaluated in correspondence of each cell of the computational Eulerian grid, located within the radial range 0 < ∕ < 0.2. The computed time scale is then adopted to estimate the Stokes number of the droplets located in the correspondent cell. Finally, the mean Stokes number is computed by averaging its instantaneous and local values along the azimuthal direction, within the considered radial range as well as over time, concurrently. The axial trend of shows a distinct hump around ∕ ≃ 15, which is followed by a gradually decrease along the streamwise direction until a unity value of is achieved around ∕ ≃ 30 and ∕ ≃ 40, at = 6, 000 and = 1000, respectively. The peak values for the Stokes number are achieved in correspondence of the peak of the clustering index. Small-scale inertial clustering rises from the competition between the droplet inertia and the Stokes drag. The drag tends to force the droplets to follow the highly convoluted paths of the turbulent motion, whereas the finite inertia of the droplets prevents them to exactly move along these material-paths. By this mechanism, largest droplets, for which ≫ 1, act as ballistic particles with respect to the smallest scales of turbulence. These droplets move across the smallest turbulent structures being only weakly perturbed and do not tend to accumulate in clusters. On the contrary, the smallest droplets, for which ≪ 1, act as passive tracers which follow exactly the path of the local turbulent motion. Also these droplets do not contribute significantly to inertial clustering. Finally, droplets for which ≃ 1 manifest an intermediate dynamics, accumulating into the small-scale interstitial vortical regions and giving rise to strong inhomogeneities of the spatial distribution of the dispersed phase. Based on the trend of , since the Stokes number keeps values significantly higher than unity in the intermediate field, we assess that inertial clustering must play a minor role in the preferential segregation of the liquid phase in this region, whereas it becomes dominant in the far field of the jet being not-negligible in the near-field, as well. Different mechanism, rather than small-scale inertial clustering, should be considered to explain the clustering of droplets in the near and intermediate field of the jet. We attribute the rising of the preferential concentration of droplets in these regions to the effect of the intermittent dynamics of the mixing layer of the jet. In the outer-spray region, high-intensity vortical structures separate the turbulent core, populated by droplets, by the dry environment, depleted of droplets. By an Eulerian point of view, this region is characterized by large temporal and spatial fluctuations of the local droplet concentration because of the intermittent dynamics of the turbulent-nonturbulent interface. The entrainment causes the envelopment of dry-air bubbles, depleted of droplets, within the turbulent jet. These dry and droplet-depleted bubbles mix with high-saturated gaseous regions full of droplets, giving rise, instant by instant, to large inhomogeneities of the droplet distributions which are, then, advected downstream. To highlight this source of clustering in the mixing layer we show in figure 8c the mean vapour mass fraction variance ⟨ ′ 2 ⟩ for the two cases. It is possible to note that the location of the peak of and ⟨ ′ 2 ⟩ in the near field region is similar. Hence we believe that the initial source of droplet clustering is induced by the intermittent dynamics of the mixing layer which create clusters that are advected downstream. Then in the intermediate and far-field regions inertial clustering becomes more important.
To quantify the importance of the droplet preferential segregation on the evaporation process, we provide a comparison between the conditional mean vapor concentration sampled by the droplets (droplet-conditional), , , and the unconditional Eulerian one, , . The latter is computed by averaging the vapor mass fraction field along the azimuthal direction and over time, concurrently. The former is defined as the average of the Eulerian vapor mass fraction field, con-ditioned to the presence of a droplet at a given point. The statistics is computed on the Eulerian grid, associating to each grid cell, the mean vapor mass fraction sampled by the droplets instantaneously located within the considered cell. Figure 9 shows the radial profiles of , and , at two different axial distances from the jet inflow, ∕ = 20 and ∕ = 30, respectively. Globally, the vapor mass fraction sampled by droplets results to be higher than the correspondent unconditional value, regardless of the Reynolds number. As already reported in [14] , this preferential sampling of the vapor mass fraction field actuated by the droplets can be related to different mechanisms. The primary contribution originates from the vapor self-produced by droplets with a tracer behavior ( ≪ 1). These droplets exactly follow the path of the turbulent motion of the surrounding vapor atmosphere, increasing the local vapor concentration in their surrounding material volume, and hence sampling higher vapor concentration. Secondarily, within a droplet cluster, independently from the mechanisms driving its formations, a highly-saturated cloud establishes in its proximity. Hence, in the presence of preferential segregation, an oversampling of the vapor concentration field is expected also by the droplets evaporating inside clusters. The radial profiles of the mean vapor concentration are similar inside the jet core, Fig.9a , whereas significant differences appear in the mixing layer.
In the core region the jet presents a nearly uniform vapor mass fraction field that nullifies the oversampling actuated by the droplets. On the other hand, in the mixing layer, the entrainment of bubble of dry air depleted of droplets causes the Eulerian, unconditional statistics to assume lower values than the conditional one. This aspect is strictly connected to the clustering mechanism previously discussed. Further downstream, the difference between the droplet-conditional vapor concentration and the unconditional one gradually increases, Fig.9b , even in the jet core. In the far-field, the inertial small-scale clustering plays a dominant role in the preferential segregation of droplets as shown in figure 8 . This lead to a considerable increment of the vapor mass fraction field sampled by the droplets grouped in clusters with respect to the unconditional statistics. Besides, no evident disagreement is observable for the unconditional mean vapor concentration profiles in = 10, 000 and = 6, 000 cases. ,and the unconditional Eulerian one, , . The enstrophythreshold conditional average is calculated by sampling the vapor mass fraction only over turbulent core events ( = 1), that is, when local enstrophy exceeds a fixed threshold. , is the vapor concentration field obtained by a conditional average on the droplet presence in a given point.
To better analyze the mixing layer of jet spray, we remind that the inner core and the irrotational outer region are separated by a nearly-sharp fluctuating interface, which is highly convoluted over a wide range of vortical scales [13] . Different approaches can be employed to identify this interface [7, 50, 21, 25] . Here, we consider the local enstrophy, 2 = ‖∇ × ⃖ ⃗ ‖. The inner turbulent core is characterized by large fluctuation of enstrophy, while in the outer region enstrophy is null. Thus, fixing a threshold, 2 ℎ , it is possible to distinguish if a given Eulerian position at fixed time is located in the turbulent jet core or in the outer irrotational environment employing the following index:
with the Heaviside function. = 1 denotes a turbulent event at Eulerian position , while = 0 an irrotational one. The conditionally average profile of vorticity magnitude has been proven to weakly depend on the threshold value in the range ℎ ≈ 0.7 0 ∕ , with 0 the mean velocity scale of the wake and the shear layer thickness [7, 13] . Hence, the index, employed here to distinguish between the rotational jet core and irrotational environment, is weakly dependent from the threshold value ℎ . In Dalla Barba and Picano [14] , it has been observed no relevant differences on the by changing the value of ℎ = 0.6 0 ∕ by a factor 2. Consistently, the same threshold value is adopted in the present work. An enstrophy conditional average for the vapor mass fraction field, , , can be defined as the mean vapor mass fraction field obtained by only sampling the instantaneous values of ( , ) linked to turbulent events, that is only for time and positions where = 1. Figure 9 provides, together , and , , the radial profile of the enstrophy conditioned , . The unconditional and both the conditional statistics upstream ∕ ≃ 20 are similar in the core, whereas in the mixing layer the radial profile of , is similar to that of , , both the statistics assuming higher values than that of , . This difference, confined to the outer spray regions, confirms the contribution of the intermittent dynamics of the mixing layer to the non-homogeneous spatial distribution of the dispersed phase discussed above. Droplets moving towards the outer region are enclosed by a high concentration vapor cloud which is expelled from the turbulent jet core. Simultaneously, dry air without droplets is engulfed in the core enhancing the fluctuation of droplets distribution and vapor concentration. Further downstream, ∕ ≃ 30, the dropletconditional vapor concentration profile shows the same shape as the enstrophy-conditional one, except for an almost constant offset, both in the core and mixing layer. This confirms how, in the far field, both the small-scale inertial clustering and the intermittent dynamics of the mixing layer contribute to the oversampling of the vapor mass fraction field actuated by the dispersed phase, the former mechanism becoming the dominant effect in the jet core. The combination of these two mechanisms, leading to the observed oversampling, strongly reduces the vaporization rate with respect to the one that would be estimated employing unconditioned statistics. The phenomenologies described above are similar for both the cases considered in this study. In particular, no significant difference is observable for unconditional vapor concentration profiles between = 6, 000 and = 10, 000.
Turbulent fluctuations are responsible for extremely different Lagrangian histories of the evaporating droplets. The joint probability density function (JPDF) of the droplet vaporization lengths and time are reported for both cases in figure 10 . The vaporization length for a single droplet is defined as the axial distance from the inlet, ∕ , where the droplet radius is decreased from ,0 to a threshold radius, , ℎ = 0.1 ,0 (99.9% of the mass evaporated). The correspondent time is defined as the droplet vaporization time, ∕ 0 . As it can be appreciated from the figure, for both = 6, 000 and = 10, 000 cases, there is a strong linear correlation between the droplet vaporization length and time, whereas different slope rates are observable. It is also remarkable how different the droplets histories are: half of ejected droplets are still present around ∕ ≃ 32 for = 6, 000 and ∕ ≃ 36 for = 10, 000, where 90% of the droplet mass is evaporated (see figure 3) . Correspondingly, the median vaporization times are about ∕ 0 = 65 and ∕ 0 = 90 for = 6, 000 and = 10, 000, respectively. Consistently with previous discussions, we note a slower average vaporization rate for the higher Reynolds number. In addition we note that the actual droplet evaporation time varies of more than 50% of its mean value indicating how heterogeneous are the dynamics. In addition, as the droplets in the = 10, 000 case have a relative large initial Stokes number, they tend to stay longer within the nearly saturated turbulent jet core, which can result in an additional delay for the vaporization process in the near field. Figure 11 (a) compares the probabil-ity density function (PDF) of the droplet radius at different ∕ positions for the two Reynolds numbers considered. We note that all PDF profiles end sharply. This is due to the fact that droplets with a radius < 0.5 are removed from the simulation, because of the numerical stability of the evolution of tiny droplets (more than 99.9% of the liquid mass is evaporated before). It is worth remarking how fast the turbulent spray is in promoting polydispersity. Though starting from an identical monodisperse status, the droplet radius distribution spreads over a wide range at = 6, 000 even at five jet radii from the inlet, whereas a slightly narrower PDF is observed at = 10, 000. Further downstream, i.e. ∕ = 20 and ∕ = 40, a similar qualitative behavior can be observed: the flat PDFs are related to a nearly constant probability to find droplets of arbitrary sizes. We attribute this fast widening of the droplet size spectrum to the strong inhomogeneous conditions of local vapour concentration. In the mixing layer, the entrainment process originates spatial regions with very low levels of saturation, which promote a fast droplet vaporization. Concurrently, because of clustering, clouds of droplets originating in the jet core maintain nearly zero vaporization rates for long time. Focusing on panel (b) of figure 11 , it is possible to appreciate how flat is the PDF of the saturation level of the droplet atmospheres. This reflects in a strongly varying droplet vaporization rate, whose PDF is shown in panel (c). The droplet vaporization rate is slightly higher for the low Reynolds number case, consistently with the slower jet evaporation length and previous discussion on the number effect. Finally, we aim to highlight that in the near and intermediate fields some rare condensation events occur. When the droplet evaporation proceeds too fast, which is more common in near-field region close to mixing layer, droplet temperature decreases resulting in the appearing of saturation conditions which can block the evaporation even in a region without high values of vapor concentration. Then, if the droplets move into a higher vapor concentration zone due to its inertia, supersaturation conditions could arise. These inertial effects are more worthy of attention for droplets with a higher Stokes number, which is the reason of the relatively broader distribution range of the saturation shown in the = 10, 000 case. The strongly polydisperse behavior shown by the evaporating droplets affects the vaporization length and time of the droplets within the spray.
To predict the temporal evolution of the droplet dynamics, the d-square law is often adopted. In this context, a linear behavior for the square diameter of each droplet is expected,
where, in many practical and experimental application, the constant can only be estimated based on the reference environmental state, e.g. bulk vapor concentration. We provide in Fig. 12 the joint probability density function of the droplet square diameter and flight time, 2 ( , 2 ). As expected, droplets do not follow a linear time behavior for the square diameter as predicted by the d-square law, otherwise all data should collapse on a single straight line. In the same figure, we also report the average droplet square diameter conditioned to the flight time, ⟨ 2 ∕ 2 ,0 ⟩( ∕ 0 ), and the droplets mean flight time conditioned to the square droplet diameter,
,0 ), for both Reynolds number cases. These two quantities have been directly extracted from the JPDF, 2 , as:
Although the two quantities are strictly related, their significance is different. The mean square diameter, eq. (25), pro-(a) (b) Figure 12 : Joint probability density function of droplet square diameter and flight time, 2 ( , 2 ), for for different Reynolds numbers. The mean square droplet diameter conditioned to the flight time (triangle symbol), ⟨ 2 ∕ 2 ,0 ⟩( ∕ 0 ), and the mean droplet flight time conditioned to the square droplet diameter (square symbol), ⟨ ∕ 0 ⟩( 2 ∕ 2 ,0 ). In the upper abscissa the time is normalized by the initial particle relaxation time, whereas the lower axis provides the non-dimensional time, ∕ 0 . vides the mean square diameter at a fixed flight time, . . after a fixed time from droplet injection. The latter, the mean flight time at fixed mean square diameter, highlights the amount of time needed, in average, for a droplet to reach a fixed diameter. In both the considered cases, the two observables assume similar values for relatively large droplets ( 2 ∕ 2 ,0 > 0.25), while strongly different for small droplets. The mean flight time conditioned to a droplet diameter shows a finite average time for a full evaporation, i.e. 2 ∕ 2 ,0 → 0: We will name it mean evaporation flight time, . Consistently with previous findings, the mean evaporation flight time increases with the Reynolds number when scaled with the advection reference time, being ∕ 0 ≃ 66 and ∕ 0 ≃ 89, for = 6, 000 and = 10, 000 respectively. We also note that, by scaling the mean flight time for the full droplet evaporation with the droplet initial relaxation time, the two cases show closer values, ∕ ,0 ≃ 106 and ∕ ,0 ≃ 86, with a slightly faster evaporation rate for the higher Reynolds number case. This behavior confirms that, by increasing the Reynolds number, two different competing effects are present: a stronger turbulence intensity which tends to fasten the evaporation rate and a higher droplet relative inertia which tends to slow down the process. We highlight how different are the behavior of the two curves, Eq. (25) versus Eq. (26), for small droplet sizes. We note that for time longer than the mean evaporation flight time, , we still have a finite average droplet diameter. This peculiar behavior can be explained considering that droplets in clusters require very long time to evaporate, hence some droplets are able to survive for very long time displaying a finite average diameter. This last aspect could be crucial in the dispersion modeling of infectious droplets [4] . To better characterize how the d-square law is able to describe present cases we provide in figure 13 the evolution of droplet square diameter predicted by the d-square law, eq. (23), estimating the constant via the bulk environmental conditions, , together with ⟨ 2 ∕ 2 ,0 ⟩( ∕ 0 ) and ⟨ ∕ 0 ⟩( 2 ∕ 2 ,0 ). The estimates for the two cases in dimensionless units are ≃ 0.115 and ≃ 0.069 for = 6, 000 and = 10, 000, respectively. In both cases the predicted droplet evaporation time is much lower than the actual one, . . the estimate of the evaporation rate is much faster than the actual one. We attribute this disagreement to the effects of clustering and preferential sampling of vapor concentration previously discussed, that is not accounted by the d-square law. Since the mean flight time, eq. (26), shows a linear behavior for small droplets and long time, we provide also the linear fitting in the form ⟨ ∕ 0 ⟩( 2 ∕ 2 ,0 ) ≃ − ( ∕ 0 ) restricted to ∕ 0 > 30. In this case the fitting constants are ≃ 1.25 and ≃ 0.0190 for = 6000 whereas ≃ 1.164 and ≃ 0.0131 for = 10000. Using these parameters we can reformulate an effective d-square law as we assume that droplets do not evaporate for some time after the injection when they travel in the saturated core and then follow an effective d-square law. The initial time when droplets are assumed not evaporating is similar for the two cases, i.e. ∼ 13 0 , which corresponds to a traveled distance around 13 jet radii.
Then the two estimated values of appear six times lower than than those estimated at environemental conditions, , for both cases. The large discrepancies between the values of and further confirm the crucial importance of improving bulk spray dispersion models accounting for the complex dynamics arising from the strong inhomogeneity of droplet and vapor mass fraction distribution in turbulent jet-sprays.
A direct numerical simulation has been conducted to investigate the motion and evaporation of inertial droplets in (a) (b) Figure 13 : Mean square droplet diameter conditioned to the flight time (triangle symbol), ⟨ 2 ∕ 2 ,0 ⟩( ∕ 0 ), and the mean droplet flight time conditioned to the square droplet diameter (square symbol), ⟨ ∕ 0 ⟩( 2 ∕ 2 ,0 ). In the upper axis the time is normalized by the initial particle relaxation time, ,0 , whereas in the lower one by the advecting time, 0 . The blue line shows the evolution of droplet square diameter predicted by employing the d-square law provided in Eq. (23) and by estimating the constant via the bulk environmental conditions. The black, dashed line provides the linear fitting of ⟨ ∕ 0 ⟩( 2 ∕ 2 ,0 ) restricted to ∕ 0 > 30. a high-Reynolds-number turbulent jet-spray ( = 10, 000). The Eulerian-Lagrangian point-droplet framework combined with a two way coupling concept, including the mass, momentum and energy exchanges between two phases, is considered for the present acetone-air spray. Monodisperse acetone droplets are continuously injected within the turbulent gaseous phase at a bulk Reynolds number = 2 0 ∕ = 10, 000. Then, a systematic and comprehensive dataset of both the instantaneous and mean fields from Eulerian and Lagrangian observables are collected, analyzed and compared to a dataset at lower Reynolds number, = 6, 000, which from physical point of view corresponds to a lower jet velocity, [14] . This comparison allowed to study thenumber effect on the turbulent evaporation of droplets in dilute sprays. Both cases show a strong evaporation rate in the mixing layer where the entrainment of dry air dilutes the saturated jet core. However, a longer evaporation length is observed for the higher Reynolds number case. This feature is related to the slower average vaporization rate observed at = 10, 000 when scaled by the jet advection reference time scale, 0 = ∕ 0 . This slower rate is attributed to the higher droplet relative inertia (Stokes number) which slows-down the vapour-liquid mass exchange.
An intense droplet clustering is apparent in both cases. This originates in the mixing layer of the near field and propagates downstream. In the mixing layer, we identify two states separated by the turbulent/non-turbulent interface: the entrained dry air depleted of droplets and saturated gas clouds full of droplets coming from the saturated core. Since droplets entering in dry entrained regions fast evaporate, while aggregates in vapor-saturated clouds cannot evaporate, these dynamics tends to intensify the clustering in the downstream evolution. In addition, in the far-field evolution small-scale clustering is also promoted by a Kolmogorov Stokes number order one for present cases. These dynamics strongly impact the Lagrangian evolution of the droplets creating very different histories. The droplet size spectrum becomes extremely wide even starting from a monodisperse distribution. This widening is associated to also wide PDF of the evaporation rate and saturation level felt by the droplets. These findings could play a role in explaining the fast widening of the droplet spectrum in warm clouds [45] . Finally, we aim to comment on the accuracy of the d-square law in approximating the droplet evolution for present cases. First, we note that the actual droplet evaporation time varies of more than 50% of its mean value indicating that actual droplet history differs from the d-square law. Second, while the mean evaporation (flight) time, , could be approximated by a d-square law, the slope appears much lower than the value extracted from environmental conditions. In addition, we note that after time much longer than a full evaporation of droplets does not occur, highlighted by a finite mean droplet size (square diameter). This peculiar aspect has been attributed to the clustering dynamics. Aggregate of droplets with a nearly saturated atmosphere may survive for longer time before a full evaporation. Hence, some droplets may show relatively large sizes for long time and distances. This behavior should be considered in the modeling of the dispersion of infectious droplets since some droplets persist much longer than what expected using the d-square law especially if based on environmental conditions. Perspective of this study is to assess the ability of Large-Eddy-Simulation models to capture this complex dynamics at relatively high Reynolds number.
This study has been funded by the research project promoted by China Scholarship Council (CSC) Grant #20180625 0023. We also thank CINECA for providing computational resources via the ISCRA C project: HiReS (HP10C3YNLC).