Evaporation of dilute droplets in a turbulent jet: clustering and entrainment effects
Federico Dalla Barba
Francesco Picano
Droplet evaporation in turbulent sprays involves unsteady, multiscale and multiphase processes which make its comprehension and model capabilities still limited. The present work aims to investigate droplet vaporization dynamics within a turbulent spatial developing jet in dilute, non-reacting conditions. We address the problem using a Direct Numerical Simulation of jet laden with acetone droplets using an hybrid Eulerian/Lagrangian approach based on the point droplet approximation. A detailed statistical analysis of both phases is presented. In particular, we show how crucial is the preferential sampling of the vapour phase induced by the inhomogeneous localization of the droplets through the flow. The preferential segregation of droplets develops suddenly downstream the inlet both within the turbulent core and in the mixing layer. Two distinct mechanisms have been found to drive these phenomena, the inertial small-scale clustering in the jet core and the intermittent dynamics of droplets across the turbulent/non-turbulent interface in the mixing layer where dry air entrainment occurs. These phenomenologies strongly affect the overall vaporization process and lead to a spectacular widening of droplets size and vaporization rate distributions in the downstream evolution of the turbulent spray.
paper.tex
%\documentclass[aip,pof,reprint,twocolumn]{revtex4-1}
\documentclass[aip,reprint,twocolumn]{revtex4-1}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{natbib}
\usepackage{subfigure}
\usepackage{color}
\def\fp#1{{\textcolor{red}{\bf ? #1 ?}}}  
\begin{document}

\title{Evaporation of dilute droplets in a turbulent jet: clustering and entrainment effects}

\author{Federico Dalla Barba}
%\affiliation{University of Padova, Department of Industrial Engineering, Via Venezia 1, 35131 Padova, Italy}
\author{Francesco Picano}
\affiliation{University of Padova, Department of Industrial Engineering, Via Venezia 1, 35131 Padova, Italy}
%\corresp[cor1]{federico.dallabarba@unipd.it}
%\corresp{francesco.picano@unipd.it}

\date{\today}

\begin{abstract}

Droplet evaporation in turbulent sprays involves unsteady, multiscale and multiphase processes which make its comprehension and model capabilities still limited. The present work aims to investigate droplet vaporization dynamics within a turbulent spatial developing jet in dilute, non-reacting conditions. We address the problem using a Direct Numerical Simulation of jet laden with acetone droplets using an hybrid Eulerian/Lagrangian approach based on the point droplet approximation. A detailed statistical analysis of both phases is presented. In particular, 
we show how crucial is the preferential sampling of the vapour phase induced by the inhomogeneous localization of the droplets through the flow. 
The preferential segregation of droplets develops suddenly downstream the inlet both within the turbulent core and in the mixing layer.
Two distinct mechanisms have been found to drive these phenomena, the inertial small-scale clustering in the jet core and the intermittent dynamics of droplets across the turbulent/non-turbulent interface in the mixing layer where dry air entrainment occurs. 
These phenomenologies strongly affect the overall vaporization process and lead to a spectacular widening of droplets size and vaporization rate distributions
 in the downstream evolution of the turbulent spray. 
\end{abstract}

\pacs{}
\maketitle  

\section{Introduction}

%

Turbulent sprays are complex multiphase flows involving unsteady and multi-scale phenomena such as turbulence coupled with phase transition processes. The presence of two distinguished phases which mutually interact exchanging mass, momentum and energy makes the description of the problem extremely challenging. If combustion is considered, chemical reactions and heat release add some clear complexities. In this scenario, a satisfactory comprehension of turbulent spray dynamics has not yet been achieved and existing models capabilities for applications are still limited~\cite{jenny2012,kruger2016effect}. Nevertheless, the progress of the research in this field is crucial for several industrial applications as well by an environmental point of view. A typical example can be found in the development of high efficiency and low emission internal combustion engines. In these applications the liquid fuel can be directly injected  into the combustion chamber where the vaporization of fuel droplets occurs together with chemical reactions within the turbulent gaseous environment.  The pollutants formation in turbulent spray combustion is related to complex multi-scale phenomenons that involve fluctuations of temperature and reactants concentrations. In particular, the soot formation occurs through a pyrolysis process in fuel-rich regions that experience high temperature without enough oxidizer to react ~\cite{raman2016modeling,attili2014formation,kennedy1997models}. This can be observed within droplets clusters, where the concentration of fuel droplets can be even thousand times higher then its bulk value leading to a peak in the fuel vapor concentration. Hence, in order to predict and model soot formation an improved understanding of the mechanisms that govern the distribution of droplets and fuel/oxidizer mixture within a turbulent jet spray is mandatory.

A phenomenological description of the overall evolution of the spray dynamics can be found in the review of Jenny et al. ~\cite{jenny2012}. The process starts with the \emph{primary atomization} of a high velocity liquid jet. As the liquid flow is ejected from a duct into a gaseous environment, interface instabilities such as Kelvin-Helmholtz and Rayleigh-Taylor, fragmentize the jet into large drops and ligaments ~\cite{marmottant2004spray}. In the downstream evolution the liquid ligaments and the drops are subjected to aerodynamic forces induced by the velocity difference at liquid/gas interface. The stresses induced by the aerodynamic forces produce a further brake-up of the liquid phase (\emph {secondary atomization}) giving origin to a system of small droplets dispersed within the gaseous turbulent phase. The atomization process occurs in a so called \emph{dense regime} and terminate when surface tension prevails on aerodynamic stresses preventing further fragmentation. At this step a \emph{dilute regime} establishes: droplets mutual interactions (e.g. collisions and coalescence) are negligible but the effect of droplets on the carrier flow is still significant~\cite{ferrante2003physical,gualtieri2015exact}. While in the dense regime the vaporization rate is negligible, in dilute conditions the vaporization process becomes significant. In this phase of the turbulent spray evolution, the most part of the liquid evaporates. Moreover, the small droplets evolve preserving a spherical shape due to the dominance of surface tension on aerodynamic stresses. Even if the presence of droplets exerts a significant effect on the flow in terms of mass, momentum and energy balance, at this step the dimension of droplets is below or comparable to the smallest scales of the turbulent flow so point-droplet approximation applies~\cite{elghobashi1994predicting}. Hence, in dilute conditions the mathematical description of droplet-laden flows lends itself particularly well to an hybrid Eulerian/Lagrangian description. The Navier-Stokes equations govern the continuous phase dynamics if distributed sink-source terms are considered in order to represent the mass, momentum and energy exchange between the Eulerian carrier phase and Lagrangian point-droplets. 
%In addition, droplets dynamic is completely described by three Lagrangian equations for position, radius and temperature.
%One of the most peculiar phenomenology that characterize droplets-laden flows in dilute conditions lay in the preferential segregation of droplets as a result of the interaction of droplets inertia with the carrier phase turbulent dynamics (see e.g. ~\cite{toschi2009lagrangian}). The mechanisms at the origin of droplets preferential segregation in free flows have been widely investigated in literature both for solid particles and liquid droplets. Bec et al. ~\cite{bec2007heavy} report the direct numerical simulations of an incompressible, particle-laden flow over a triply periodic domain. The results confirm that the intensity of droplet segregation at smaller scale of turbulent motion is maximum when the particle relaxation time, $\tau_p$, equals the Kolmogorov time scale, $\tau_eta$, such that the Stokes number, $St_\eta$ is unitary. Reveillon \& Vervisch \cite{reveillon2007effects} observe the same mechanisms for evaporating droplet segregation in statistically stationary, isotropic homogeneous turbulence. Concerning turbulent sprays, a mean accumulation of the dispersed phase has been observed at specific distances from the inflow both experimentally (~\cite{au2016effect}) and numerically (~\cite{casciola2010dynamics}). Even if droplets preferential concentration has been well characterized, the effect of this phenomenology on the overall vaporization process within turbulent sprays is still object of research and constitute one of the main focuses of this paper.  

One of the most peculiar phenomenology that characterize dispersed multiphase turbulent flows in dilute conditions is the preferential segregation of particles/droplets as a result of the interaction of their inertia with the carrier phase turbulent dynamics (see e.g. ~\cite{toschi2009lagrangian}). The mechanisms at the origin of preferential segregation in free flows have been widely investigated in literature both for solid particles and liquid droplets. The intensity of small scale segregation of solid particles in homogeneous turbulence is found to be driven by the Stokes number $St_{\eta}=\tau_p/\tau_{\eta}$, with $\tau_p$ the particle relaxation time and $\tau_{\eta}$ the Kolmogorov time scale. More specifically, the intensity of small scale clustering is maximum when the particle relaxation time is of the order of the Kolmogorov time scale such that the Stokes number, $St_\eta\simeq1$ \cite{calzavarini2008quantifying, bec2007heavy, saw2008inertial, gualtieri2009anisotropic}. The same behavior is observed for evaporating droplets which behave as inertial particles\cite{reveillon2007effects}. Concerning turbulent jets, a mean accumulation of the dispersed phase has been observed at specific distances from the inflow both experimentally ~\cite{lau2016effect} and numerically ~\cite{picano2010anomalous}. Even if  preferential concentration of a dispersed phase has been well characterized in homogeneous turbulence, the effect of this phenomenology on the overall vaporization process within turbulent jet sprays is still object of research and constitute one of the main focuses of this paper.  

By a theoretical and numerical point of view one of the first description of the vaporization of spherical droplets dragged by a gaseous phase flow was addressed in the seminal works of Spalding and Godsave ~\cite{spalding1953, godsave1953}. Fixing the environmental vapor concentration, they found that droplets surface decreases linearly with time ($D^2$ law). Maxey \& Riley, ~\cite{maxey1983equation} report an equation for the motion of a finite size spherical particle/droplets at low droplet Reynolds number. The equation accounts for the Stokes drag,  added-mass effect and buoyancy force. Dealing with a dispersed phase with a density much higher than that of the fluid the Stokes drag is sufficient to describe its dynamics ~\cite{armenio2001importance, olivieri2014effect}.  Abramzon \& Sirignano ~\cite{abramzon1989droplet} proposed an improved model for droplets vaporization, considering non-uniform and time-dependent environmental conditions, taking into account forced convection, molecular diffusion and the Stefan flow contribution to the vapor transport from droplet surface to neighboring environment. Even if several models ~\cite{jenny2012, marchioli2008some} can be found in literature for the simulation of turbulent evaporating sprays in Reynolds Averaged Navier-Stokes (RANS) or Large-Eddy Simulation (LES) frameworks, these models lack in capabilities to accurately reproduce complex phenomena such as droplets small scale clustering~\cite{kruger2016effect}. Despite the highly demanding computational resources, the use of Direct Numerical Simulation (DNS) allows to capture the whole physics of the spray vaporization process in order to understand the complex phenomenons involved. In this context,  Mashayek ~\cite{mashayek1998direct} adopted an Eulerian/Lagrangian approach in order to perform a DNS of low Mach number, homogeneous shear turbulent flows laden with droplets . Miller \& Bellan ~\cite{miller1999direct} report a DNS of a confined three-dimensional, temporally developing gas mixing layer laden with evaporating hydrocarbon droplets at subsonic Mach number. Reveillon and co-workers ~\cite{reveillon2007effects} studied the effect of preferential droplet accumulation on the evaporation in isotropic turbulence showing that different regime takes place depending on droplet concentration. %Droplets preferential segregation is observed at interface between the free-stream and the droplet laden flow. 
Recently, Bukhvostova et al. ~\cite{bukhvostova2014comparison} consider the DNS of a turbulent channel flow of a mixture of air and water vapor laden with water droplets. The work focuses on the comparison between the performances of an incompressible and a low Mach number asymptotic formulation in reproducing the flow dynamics. Even if the two formulations show good qualitative agreement, the low Mach formulation is found to be crucial in order to obtain a reliable quantitative predictions of heat and mass transfer.

The prototypal flow for an evaporating spray is constituted by a turbulent free jet which is characterized by the effect of environmental gas entrainment. In more details, a turbulent jet is constituted by a rotational turbulent core which is continuously entrained by the surrounding irrotational fluid ~\cite{da2014interfacial}. In sprays, the entrained dry flow dilutes the vapor concentration and controls the vaporization process. This phenomenology was found to be of critical importance also in natural phenomena. One example consists in the effect of entrainment in stratocumulus clouds where it constitutes a driving parameters in the determination of cloud lifetimes and  in turns even regulating planetary-scale properties ~\cite{mellado2017cloud}. The fast grow rate of the droplet size spectrum in warm cloud is a challenging, still not understood, problem in cloud physics ~\cite{falkovich2002acceleration, sardina2015continuous}, despite its importance in determining raining conditions. 

To the best of authors knowledge, a fundamental study on the effects of the entrainment in an evaporating turbulent jet spray together with an analysis of  the preferential segregation effect is still missing. This work aims to cover this lack considering DNS data of an evaporating turbulent spray using a 2-way coupling approach between the two phases and accounting for the entrainment effect. The  numerical algorithm adopts a hybrid Eulerian/Lagrangian approach and point droplets approximation. In addition, the effect of density variation is accounted by a low Mach number formulation of the Navier-Stokes equations.  A strong preferential segregation of droplets is observed over the whole downstream evolution of the spray which induces a preferential sampling of vapor concentrated regions. Two different mechanisms are found to drive this process, the former is due to  inertial clustering, the second is related to the dynamics of the jet 
entrainment. This last mechanism is found to be crucial in the outer part of the jet core where the evaporation peaks and strongly impacts the vaporization dynamics which is characterized by a strong widening of the droplet size spectrum. 

\section{Numerical Method}
In this paper we report a direct numerical simulation of a turbulent evaporating spray in an hybrid Eulerian/Lagrangian framework adopting the point droplet approximation in the 2-way coupling conditions. The governing equations for the Eulerian gaseous phase consist in a low Mach number formulation of the Navier-Stokes equations in an open environment where arbitrary density variations can be accounted neglecting acoustics~\cite{majda1985derivation}. Consistently with previous studies ~\cite{miller1999direct, bukhvostova2014comparison}, the effect of the dispersed phase on the gaseous phase is accounted by sink/source coupling terms appearing in the right hand side of the mass, momentum and energy equations, $S_\rho$, $\bf{S_m}$ and $S_e$ (see ~\cite{mashayek1998direct}),

\begin{align}
\label{eq:1}
& \frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \mathbf{u})=S_\rho
\\
& \frac{\partial}{\partial t}(\rho Y_V)+\nabla \cdot (\rho Y_V \mathbf{u}) = \nabla \cdot (\rho {\cal D} \nabla Y_V)+S_\rho
\\
& \frac{\partial}{\partial t}(\rho \mathbf{u})+\nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = \nabla \cdot \pmb{\tau} - \nabla {P} + \mathbf{S_m} 
\\
& \nabla \cdot \mathbf{u} = \frac{1}{p_0} \left[ \frac{\gamma-1}{\gamma} \nabla \cdot (k \nabla T) + \ S_e \right]
\\
\label{eq:2}
& p= \rho \  R_G (1+M \ Y_V) T
\end{align}

where $\rho$, $\bf u$, $p_0$, $P$ and $T$ are respectively the density, the velocity, the thermodynamic pressure, the hydrodynamic pressure and the temperature of the carrier vapor/gas mixture that will be denoted in the following as the carrier phase. The vapor mass fraction, $Y_V$, is defined as the ratio of the vapor partial density and the total density of the carrier phase, $Y_V=\rho_V/\rho$. The vapor/gas binary diffusion coefficient, the thermal conductivity of the mixture and the specific constant of gas are denoted respectively as $\cal D$, $k$ and $R_G$. The parameter $M$ is defined as $M=W_G/W_L-1$ with $W_G$ and $W_L$ the molar weight of the gas and liquid phases. The ratio of the constant pressure coefficient, $C_P$, and the constant volume coefficient, $C_V$, of the mixture is denoted by $\gamma$. The viscous stress tensor is $\pmb{\tau} = 2 \mu (\nabla \mathbf{u} + {\nabla \mathbf{u}}^T)-2\mu_b/3 \nabla \cdot \mathbf{u} \mathbf{I}$, with $\mu$ and $\mu_b$ the dynamic and bulk viscosities.  It should be remarked that the thermodynamic pressure $p_0$ is constant in space, due to the low-Mach number asymptotic expansion~\cite{majda1985derivation} and in time, due to the open space conditions.

Droplets are treated as rigid evaporating spheres and the liquid phase properties (e.g. temperature) are assumed to be uniform inside the droplets. Droplet rotation, distortion and mutual interactions (e.g. collision, coalescence) are neglected considering dilute volume fractions. 

\begin{align}
\label{eq:3}
& \frac{d \mathbf{u_d}}{dt}=\frac{\mathbf{u}-\mathbf{u_d}}{\tau_d}
\\
& \frac{dm_d}{dt}=-\frac{1}{3}\frac{m_d}{\tau_d} \frac{Sh}{Sc} \ ln(1+B_m)   
\\
\label{eq:4}
& \frac{dT_d}{dt} = \frac{1}{3 \ \tau_d} \left(\frac{Nu}{Pr} \frac{C_{P,G}}{C_L}(T-T_d)- \frac{Sh}{Sc} \frac{L_V}{C_L} \ ln(1+B_m)\right)
\end{align}

where $\mathbf{u_d}$, $m_d$ and $T_d$ are droplet velocity, mass and temperature, $C_{P,G}$ and $C_L$ are the constant pressure coefficient of the gas and the liquid specific heat and $\tau_d=2 \rho_L r_d^2 / (9 \mu)$ is the droplet relaxation time, with $\rho_L$ the liquid phase density. The mass diffusivity and the thermal conductivity are accounted through the Schmidt and Prandtl numbers respectively, $Sc=\mu / (\rho \cal D)$ and $Pr=\mu / (C_p k)$. The Nusselt number, $Nu_0$, and the Sherwood number, $Sh_0$, are estimated as a function of the droplets Reynolds number, $Re_d=\rho || \pmb{u} - \pmb{u}_d || / \mu $, according to the Fr\"{o}ssling correlation:

\begin{align}
& Nu_0 = 2+0.552 \ Re_d^{\frac{1}{2}} \ Pr^{\frac{1}{3}} 
\\
& Sh_0 = 2+0.552 \ Re_d^{\frac{1}{2}} \ Pr^{\frac{1}{3}}
\end{align}

A correction is then applied to $Nu_0$ and $Sh_0$ in order to account for the Stefan flow \cite{abramzon1989droplet}:

\begin{align}
& Nu = 2+\frac{(Nu_0-2)}{F_T},  & F_T= \frac{(1+B_t)^{0.7}}{B_t} \ ln(1+B_t)
\\
& Sh = 2+\frac{(Sh_0-2)}{F_M},  & F_M= \frac{(1+B_m)^{0.7}}{B_m} \ ln(1+B_m)
\end{align}

The parameter $B_m$ and $B_t$ are the Spalding mass and heat transfer number respectively, the former being the driven parameter for the vaporization rate,
%Each variable is non-dimensional according to the the following reference scales evaluated at the jet inlet section: $\rho_{ref}$, $p_{ref}$, $U_{ref}$, $T_{ref}$, $\mu_{ref}$, $k_{ref}$, $C_{P,ref}$, $R_{ref}$ and $D_{ref}$ which are respectively the reference density, pressure, velocity (inlet bulk velocity), temperature, dynamic viscosity, thermal conductivity, heat capacity, length (inlet radius) and mass diffusivity. 
%All fields appearing in the equations consists of the zero order terms of the following expansion:
%\begin{equation}
%\bf f(\bf{x},t,Ma)=\bf f_0(\bf{x},t)+ \bf f_1(\bf{x},t) \ Ma + \bf f_2(\bf{x},t) \ Ma^2 +\textit{O}(Ma^3)
%\end{equation}
%where $\bf f$ is a generic scalar or vectorial field. The only exception is the term $p_2$ which is the second order term of the pressure expansion, being $\nabla p_0=0$ and $\nabla p_1=0$. The subscript zero has been omitted in equations. 
%The main non-dimensional groups are the bulk Reynolds number, $Re=\rho_{ref} U_{ref} R_{ref} / \mu_{ref}$, the Schmidt number, $Sc=\mu_{ref} / (\rho_{ref} D_{ref})$, the Prandtl number, $Pr=\mu_{ref} C_{P,ref} / k_{ref}$ and the Stokes number, $St=\tau_d/t_{ref}$, with $\tau_d=2 \rho_L r_d^2 / \mu$ the droplet relaxation time and $t_ref=R_ref/U_ref$ the reference time scale based on jet inlet bulk velocity and radius. The carrier gaseous phase is constituted by a mixture of vapor and gas can be characterized by the vapor mass fraction field, $Y_v=\rho_v/\rho$, with $\rho_v$ the vapor partial density and $\rho$ the total mixture density. $Y_v$ is bounded between zero and $Y_{v,s}$ corresponding to the dry and saturation level, respectively.
%The remaining terms are the vapor mass fraction, $Y_V$, the specific heat of the gaseous and liquid phases, $C_{P,G}$ and $C_L$, the latent heat of vaporization of the liquid phase, $L_V$, $M=W_G/W_L-1$ with $W_G$ and $W_L$ the molar weight of the gas and liquid phase, and the reference heat capacity ratio, $\gamma=C_{P,ref}/C_{V,ref}$. The Sherwood and the Nusselt numbers are evaluated as that for a rigid sphere by the Fr\"{o}ssling correlation using droplets Reynolds number, $Re_d=\rho || \mathbf{u} - \mathbf{u_d} || 2 r_d / \mu$ \cite{}. 

\begin{align}
& B_m=\frac{(Y_{V,s}-Y_V)}{(1-Y_{V,s})}
\\
& B_t=\frac{C_{P,V}}{L_V}(T-T_d)
\end{align}

where $C_{P,V}$ is the constant pressure coefficient of the vapor, $Y_V$ is the vapor mass fraction in the carrier phase evaluated at droplet center and $Y_{V,s}$ is the vapor mass fraction evaluated at droplet surface. This latter corresponds to mass fraction of vapor  in a saturated vapor/gas mixture at droplet temperature. In order to estimate $Y_{V,s}$ we assume the equilibrium hypothesis such that the Clausius-Clapeyron relation applies:
 
\begin{equation}
\chi_{V,s}=\frac{p_{ref}}{p_0} exp \left[ { \frac{L_V}{R_V} \big( \frac{1}{T_{ref}}-\frac{1}{T_d} \big)} \right]
\end{equation}

with $\chi_{V,s}$ the vapor molar fraction at droplet surface, p the thermodynamic pressure, $p_{ref}$ and $T_{ref}$ arbitrary reference pressure and temperature and $R_V$ the vapor gas constant. The saturated vapor mass fraction is then:

\begin{equation}
Y_{V,s}=\frac{\chi_{V,s}}{\chi_{V,s}+(1-\chi_{V,s}) \frac{W_g}{W_L}}
\end{equation} 

%The source/sink coupling terms can be expresses as in the following ~\cite{}:

%\begin{align}
%\label{eq:5}
%& S_{\rho}=-\sum_i \frac{d}{dt}(m_{d,i}) \delta(\pmb{x}-\pmb{x}_{d,i})
%\\
%& \pmb{S_m}=-\sum_i\frac{d}{dt}(m_{d,i} \pmb{u}_{d,i}) \delta(\pmb{x}-\pmb{x}_{d,i})
%\\
%\label{eq:6}
%& S_e=-\sum_i \left[ \frac{d}{dt} (m_{d,i}  \ C_L \ T_{d,i}) - L \ \frac{d}{dt}(m_{d,i}) \right]  \delta(\pmb{x}-\pmb{x}_{d,i})
%\end{align} 

%where the summation is taken over the whole droplets in the domain.

We have performed the Direct Numerical Simulation of the evaporation of liquid acetone droplets dispersed within a turbulent air/acetone vapor mixture. 
 The numerical code is constitute by two different modules. An Eulerian algorithm directly evolves the gaseous phase dynamics solving the Low-Mach number formulation of the Navier-Stokes equations \eqref{eq:1}-\eqref{eq:2} (see e.g. \cite{picano2011dynamics,rocco2015curvature} and references therein for validation and tests). A second order central finite differences scheme is adopted on the staggered grid for space discretization, while temporal evolution is performed by a low-storage third order Runge-Kutta scheme. A Lagrangian solver evolves droplets mass, momentum and temperature laws \eqref{eq:3}-\eqref{eq:4}. The temporal integration uses the same Runge-Kutta scheme of the Eulerian phase and second-order accurate polynomial interpolations are used
to calculate Eulerian quantities at droplet positions. We have preliminarily tested the evaporation dynamics in the numerical code in two different test cases. The former case concerns a liquid water droplet carried by a laminar dry jet. In this extremely dilute conditions the Spalding $D^2$ law is a valid analytical solution for the droplet radius evolution over time (Fig. \ref{fig:1:a}). In the latter case, a water droplets freely falling in wet air is considered and the numerical solution for the droplet temperature evolution is compared to an experimental dataset (Fig. \ref{fig:1:b}). 

\begin{figure}[h]
\subfigure[]{\includegraphics[width=0.47\columnwidth]{art/drop_radius_vs_time2.pdf} \label{fig:1:a}} \
\subfigure[]{\includegraphics[width=0.47\columnwidth]{art/drop_temp_vs_time2.pdf} \label{fig:1:b}} \
 \caption{(a) Evolution over time of the radius of a single water droplet in a dry air laminar jet. The ambient pressure and temperature are set to $p= 101300 Pa$ and $T=273.15 K$. The inlet radius and the bulk velocity of the jet are respectively $R=5 \ 10^{-3} m$ and $U=1.7 m/s$. The droplet is injected at the local carrier phase velocity. The initial droplet radius is set to $r_{d,0}=5 \mu m$ and the temperature to $T_{d,0}=273.15 K$.  In the figure, the droplet radius and the time are scaled respectively by the droplet initial radius and the reference time scale $t_0=R/U$. The continuous line represents the analytical solution computed by the Spalding $D^2$ law, $r_d=\sqrt{r_{d,0}^2-k t}$, where $k=2 \rho {\cal D} \ ln(1+B_m) / \rho_L$. Dots represents the numerical results. (b) Temperature evolution over time of a water droplet freely falling in air at pressure $p=101300 Pa$, temperature $T=301.45 K$ and relative humidity $\chi=0.22$. Droplet is initially at rest with an initial temperature, $T_{d,0}$, equal to the environmental air one. In the figure the temperature is scaled by the initial droplet temperature, $T_{d,0}$, while time is scaled by droplet thermal relaxation time, $\tau_t$. This latter is defined as the time required by droplet temperature to change by the 63 \% of its total change between initial temperature and regime temperature. The continuous line represents the result of the simulation while dots report an experimental dataset ~\cite{kinzer1951}. The regime temperature of droplet corresponds to the ventilated wet bulb temperature at prescribed environmental pressure and actual temperature, which is $T_{wb}=288.15 K$.}
\label{fig:1}
\end{figure}

The computational domain consists of a cylinder extending for $2\pi\times22\,R\times70R$ in the azimuthal, $\theta$, radial, $r$ and axial, $z$, directions. The domain has been discretized by $N_\theta\times N_r\times N_z=128\times225\times640$ points using a non-equispaced, staggered mesh in the radial and axial direction. The mesh has been stretched in order to be of the order of the Kolmogorov length in the downstream evolution. The flow is injected at the center of one base of the cylindrical domain and  streams out at the other end. Time-dependent inflow boundary conditions are prescribed. A fully turbulent velocity is assigned at the jet inflow section (Dirichlet condition) by means of a cross-sectional slice of a fully developed companion pipe flow DNS. The flow is injected through a center orifice while the remaining part of the domain base is impermeable and adiabatic.  A convective  condition is adopted at the outlet and  an adiabatic traction-free condition is prescribed at the side boundary. This side boundary condition makes the entrainment of external fluid possible which
in the present case is dry air.
The gas/vapor mixture is injected at a bulk velocity $U_0=9 m/s$ through a nozzle of radius $R=5 \ 10^{-3} m$. The ambient pressure is set to $p=101300 Pa$ while the initial temperature is fixed to $T=275.15 K$. The injection flow rate of the gas is kept constant fixing a bulk Reynolds number $Re=2 U_0 R/\nu=6000$, with $\nu=1.5 10^{-5} m^2 / s$ the  kinematic viscosity. At the inflow section a near saturation condition is prescribed for the air/acetone vapor mixture, $S=Y_V/Y_{V,s}=0.99$, with $Y_{V,s}$ the vapor mass fraction saturation level at the actual temperature. The acetone mass flow rate is set by the mass flow rate ratio $\Phi_m=\dot m_{act}/\dot m_{air}=0.23$, with $\dot m_{act}=\dot m_{act,L}+\dot m_{act,V}$ the sum of liquid and gaseous acetone mass flow rates and $\dot m_{air}$ the gaseous one. Liquid  monodisperse droplets  with radius $r_{d,0}=6 \mu m$ are injected within the saturated vapor carrier phase. All the droplet characteristics have been chosen to reproduce acetone liquid. The injected droplets are distributed randomly over the inflow section with initial velocity equals to the local turbulent gas phase velocity. Before the injection of droplets, the simulation is started considering only the single-phase flow until statistical steady conditions have been attained. Then the simulation is run for about $200R/U_0$ time scales in order to reach a statistical steady condition for the two-phase evaporating flow before to collect the dataset. The statistics considers around one hundred samples separated in time by $R/U_0=1$. 


\section{Results and discussions}

\begin{figure}[t]
\subfigure[]{\includegraphics[width=0.248\textwidth]{art/spray_z_208000_cmp.png} \label{fig:2:a}} \
\subfigure[]{\includegraphics[width=0.198\textwidth]{art/zoom_2.png} \label{fig:2:b}} \
 \caption{(a) A radial-axial slice of the turbulent spray. The black points represent a subset of droplets formed by 1/5 of the whole population. Only droplets located inside a slice of width $w/R=0.01$ are visible. Each point size is proportional to the corresponding droplet radius (scale factor 100). The carrier phase is contoured according to the instantaneous vapor mass fraction field, $Y_V$, which  is bounded between $0$ and $0.18$, the former corresponding to the dry condition and the latter to the 99\% saturation level prescribed at inlet. (b) Enlargements of two different jet regions centered at $z/R=10$ (lower panel) and $z/R=20$ (upper panel).}
\label{fig:2}
\end{figure}

%The main mechanism driving this phenomenon relies in the fact that droplets finite inertia prevents them to follow highly convoluted local turbulent motion. More specifically, in a turbulent droplet-laden flow vortexes act as centrifuges which eject droplets heavier then the fluid away from vortexes core such that lighter droplets are entrapped in coherent vortical structures while heavier ones tend to accumulate inside convergence regions. The small-scale droplets distribution is governed by the Stokes number, $St_\eta$ , which is defines as the ratio of droplets response time,$\tau_d$, with the characteristic time of the small turbulent eddies $\tau_eta$. 
%In order to highlight by a phenomenological point of view some of the most relevant phenomenons that affect the evolution of the spray dynamics, the instantaneous vapor mass fraction field and droplets population distribution are considered. 
%In the following the analysis of the DNS dataset is reported evaluating both the average Eulerian fields of the gaseous carrier phase and the distribution of average droplets Lagrangian variables. To this purpose an ensemble average over both time and azimuthal direction is considered. The averaging of the Eulerian fields is computed as the arithmetic average of fields evaluated on the mesh nodes over the whole range of azimuthal extension of the domain and time steps after the establishment of steady conditions ,  $[0,2\pi]$ and $[t_{steady}, t_{max}]$. Concerning the averaging of Lagrangian variables the same ranges are considered. Droplets are located within mesh cells, then the same process is applied evaluating Lagrangian variables as a propriety of each mesh cell.   
%The main mechanism driving droplets preferential segregation relies in the competition between inertial forces and aerodynamic forces induced on droplets by the turbulent motion. The aerodynamic forces tend to trail droplets according to the highly convoluted local turbulent structures while droplets finite inertia prevents them to follow exactly the turbulent flow motion. By this mechanism droplets heavier then the fluid tend to be ejected from vortexes core while lighter droplets are entrapped in coherent vortical structures~\cite{}. The small-scale droplets distribution is governed by the Stokes number, $St_\eta$ , which is defined as the ratio of droplets response time, $\tau_d$, and the characteristic time of the small turbulent eddies, $\tau_eta$. Droplets with $St_\eta \gg 1$ act as ballistic particles that can move across turbulent structures being only weakly perturbed. These droplets are subjected to a weak preferential segregation in low vorticity regions. On the opposite, droplets with  $St_\eta \ll1$ act as passive tracer which move according to the local turbulent motion. These droplets tend to segregate in high vorticity region, nevertheless turbulent micro-mixing diffuses them constituting dispersed clouds. The intensity of preferential segregation is maximum when $St_\eta = 1$. In this condition droplets entrapped into low-vorticity regions are not heavy enough to move away, while the effect of micro-mixing is not strong enough to disperse droplets entrapped into high vorticity areas. 

A general overview of the instantaneous vapor mass fraction field and droplet population distribution is provided in Fig. \ref{fig:2}. Droplets populate only the turbulent jet core, while are not present in the outer dry region. The distribution of droplets is strongly inhomogeneous and clustering is also apparent. In particular, droplets seem to preferentially segregate in regions characterized by high vapor concentration while only few isolated droplets can be found in poorly saturated areas (see the enlargements shown in Fig. \ref{fig:2:b}). There are several mechanisms driving droplets preferential concentration in a turbulent flow, the most relevant of them being the small-scale clustering ~\cite{eaton1994preferential,toschi2009lagrangian}, droplets accumulation along jet axis ~\cite{picano2010anomalous,lau2016effect} and as we will show the droplet dynamics in the mixing layer. The small scale clustering concerns the interaction of droplets distribution with the smallest structures of turbulence that promotes the segregation if the dissipative flow time scale is of the same order of the particle/droplet (inertial) time scale. In turbulent jets or sprays, a mean accumulation of the dispersed phase has been observed at specific distances from the inflow. This location is determined by the matching of particle time scale and the local large scale jet time which quadratically decreases with the downstream axial distance ~\cite{picano2010anomalous,lau2016effect}. 

%In the downstream evolution of the spray droplets evaporate reducing their inertia, hence droplets dynamics tends to progressively change from that of ballistic particles to that of passive tracers. This evolution leads to an accumulation of droplets in the spray core at some axial distance downstream the jet inflow \cite{picano2011particle}. 
%Finally, turbophoresis consists in the tendency of inertial particles to move away from regions characterized by high turbulent kinetic energy thus accumulating in lower energy regions.
Independently by the mechanisms driving droplets segregation, the vapor mass fraction increases rapidly inside clusters of evaporating droplets due to their high concentration. As the vapor concentration increases, the local evaporation rate is reduced. The vaporization process may be even completely blocked if the vapor concentration reach the saturation level, $Y_{V,s}$. In this case a non evaporating, fully saturated core appears inside the cluster ~\cite{reveillon2007evaporating}. It is then clear how the clustering phenomenon can strongly affect the overall vaporization process (e.g. evaporation length) by reducing locally the rate of vaporization. 

The Fig. \ref{fig:3:a} provides the average distribution of liquid mass fraction over the spray, $\Psi=m_L/m_{G}$, where $m_L$ is the mean mass of liquid acetone and $m_{G}$ is the mean mass of the gaseous carrier phase inside an arbitrary small control volume, $\Delta V$. The overall vaporization length can be defined as the axial distance from inlet where the 99$\%$ of the injected liquid mass has transit to vapor phase. According to this definition, the vaporization process is terminated at about $z/R\simeq43$. It should be noted that the liquid phase mass fraction is significantly higher in the spray core then in the outer layer, such that the turbulent gaseous phase dynamics is especially affected by the presence of the dispersed phase in this region. This is consistent with the distribution of the average saturation field, $S$, reported in Fig. \ref{fig:3:b}. We note an almost-saturated flow near the inlet caused by the prescribed inflow conditions, then the saturation level gently decreases in the downstream evolution while sharply towards the outer jet region. The turbulent spray is constituted by a spreading and slowly decaying turbulent core which is surrounded by the dry and irrotational environmental air. The turbulent core is continuously entrained by the environmental air which mixes up with the turbulent air/vapor mixture thus reducing the vapor concentration. Since the inner core fluid cannot reach the outer region, the spray core shows higher saturation level over the whole downstream evolution of the flow. The effect of dry air entrainment is  crucial on the overall vaporization process. The dilution of vapor concentration is indeed fundamental in order to allow the vaporization process to advance.

\begin{figure}[t]
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/drop_mass_frac_2.png} \label{fig:3:a}} \
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/Sv_2.png} \label{fig:3:b}} \
\caption{(a) Average liquid phase mass fraction, $\Psi=m_L/m_{gas}$ where $m_L$ is the mean mass of liquid acetone and $m_{G}$ is the mean mass of the gaseous carrier phase.% inside an arbitrary small control volume, $\Delta V$. 
The labels show different distances from the jet inlet, $z/R$, in correspondence of which the 50\%, 90\%, 99\% and 99.9\% of the injected liquid mass is evaporated. (b) Average saturation field, $S=Y_{V}/Y_{V,s}$, where $Y_V$ is the actual vapor mass fraction field and $Y_{V,s}=Y_{V,s}(p,T)$ is the value of vapor mass fraction corresponding to the local saturation condition.}
\label{fig:3}
\end{figure}

The average droplets radius and evaporation rate distributions are reported respectively in Fig. \ref{fig:4:a} and Fig. \ref{fig:4:b}. According to the discussed entrainment effect, the vaporization rate is maximum in the mixing layer separating the jet outer and core regions. The peak value is reached in the shear layer immediately downstream the inflow section, where large droplets enter in direct contact with the dry environmental air. Consequently, at each axial distance form inlet, larger droplets can be found in the spray core where the vaporization process is slowed down by the high vapor concentration, while smaller droplets can be found towards the outer region where the vaporization proceeds faster.

\begin{figure}[t]
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/drop_r_2.png} \label{fig:4:a}} \
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/drop_evp_rate_2.png} \label{fig:4:b}} \
\caption{(a) Average droplet radius rescaled by the initial droplets radius $r_{d,0}$. (b) Average droplet vaporization rate divided by the reference scale defined as $\dot{m_{d,0}}=m_{d,0}/\tau_{d,0}$  with $m_{d,0}$ the initial droplet mass and $\tau_{d,0}$ the initial droplet relaxation time.}
\label{fig:4}
\end{figure}

The transition of liquid phase to vapor phase requires an amount of energy per unit mass equal to the latent heat of vaporization of acetone so that the power required by the vaporization is proportional to the evaporation rate. The overall energy required by the vaporization process is provided by the internal energy of both the gaseous carrier phase and the liquid dispersed phase, thus resulting in an overall cooling of the spray in the downstream evolution. The average distribution of the carrier phase and the droplets temperature is reported in Fig. \ref{fig:5}. In the outer spray region the smaller droplets surrounded by low-saturated gas are colder than core droplets due to the higher evaporation rate. Nevertheless, the carrier phase shows an opposite behavior: the spray core is sensibly colder then the shear layer and the minimum gas temperature can be observed between $z/R=25$ and $z/R=45$. This phenomenon is due to the distribution of the liquid phase mass fraction. In the spray core the liquid mass represent a significant part of the overall spray mass. Hence, the cooling effect due to vaporization is much more intense in this region where a large amount of droplets slowly evaporate.

\begin{figure}[t]
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/T_2.png} \label{fig:5:a}} \
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/drop_temp_2.png} \label{fig:5:b}} \
\caption{(a) Average gas phase temperature, $T$, rescaled by the injection temperature, $T_0$. (b) Average droplet temperature, $T_d$ rescaled by droplet initial temperature, $T_{d,0}$. The reference temperature scales are equal $T_d=T_{d,0}$.}
\label{fig:5}
\end{figure}

We have highlighted the existence of a strong preferential segregation of droplets, focusing in particular on the effect of this inhomogeneous distribution on the overall vaporization process dynamics. The intensity of droplets segregation can be measured in different ways, e.g.\cite{shaw2002towards}. 
%We will start the analysis using the clustering index $K$ which is useful to determine if clustering is present in each point of a inhomogeneous field\cite{}.  
We will measure the intensity of clustering in each point of the inhomogeneous turbulent jet spray\cite{battista2011intermittent} using the clustering index $K$, which is defined as

\begin{equation}
K=\frac{\overline{(\delta n)^2}}{\overline{n}}-1
\end{equation}

where $\overline n$ and $\overline{(\delta n)^2}$ are the mean and variance of the number of droplets in given small volume $\Delta V$. If droplets are
completely randomly located, their distribution is determined by a Poisson process in which mean and variance coincide. Hence if clustering is not present
and particles are random distributed $K=0$. %if segregation is present $K>0$. 
%
%by the comparison of the actual distribution of droplets population with a reference distribution of random independent points, this latter being described by a Poisson distribution:
%%
%%\begin{equation}
%%p(n,\overline{n})=\overline{n}^n e^{-n} / n! 
%%\end{equation}
%
%For a single event of a Poisson process, p is the probability to find n random points located inside a domain $\Omega$ of volume $\Delta V$, $\overline n$ being the average number of points over the domain $\Omega$.  
%For a Poisson process the variance $\overline(\delta)^2$ is equal to the mean $\overline{n}$. 
%Where $K=0$ droplets are distributed according to a Poisson distribution and no preferential segregation can be observed. 
On the opposite, if $K > 0$  the variance exceeds the mean value indicating that droplets preferentially segregate in clusters. Fig. \ref{fig:6:a} shows the clustering index $K$ computed over the whole spray domain. The large positive value of $K$ corresponds  to strong preferential segregation of droplets. 
We note that droplets are initially random distributed and then develop clustering. In particular, near the inflow $K$ assumes positive value only in the mixing 
layer where the local droplet concentration is intermittent because of the fluctuation of the turbulent/non-turbulent interface which separates the turbulent core populated by droplets and the outer region without droplets. %As it will be shown soon, this dynamics is crucial to determine the vaporization rate. 
It should be noted that the air regions entrained from the environment in the core are almost droplet-free and enhance the fluctuation level of 
the droplet concentration even in the jet core, see the snapshots reported in Fig~\ref{fig:2}. 
Downstream the clustering appears in the whole turbulent jet core. We attribute this phenomenon also to  the developing of small-scale turbulent clustering.  

The main mechanism driving the small scale clustering relies in the competition between inertia and Stokes drag. The drag tends to trail droplets according to the highly convoluted local turbulent structures while droplets finite inertia prevents them to follow exactly the turbulent flow motion. By this mechanism droplets heavier then the fluid tend to be ejected from vortex cores ~\cite{bec2007heavy}. The small-scale droplet distribution is governed by the Stokes number, $St_\eta$ , which is defined as the ratio of droplets response time, $\tau_d=2 \rho_L r_d^2/(9 \rho \nu)$, and the characteristic time of the dissipative scales, $\tau_\eta=(\nu/\varepsilon)^{1/2}$. Droplets with $St_\eta \gg 1$ act as ballistic particles that move across turbulent structures being only weakly perturbed and showing a negligible preferential segregation. On the opposite, droplets with  $St_\eta \ll1$ act as passive tracers which move according to the local turbulent motion without exhibiting clustering. Preferential segregation is maximum when $St_\eta \sim 1$. 
%In this condition droplets entrapped into low-vorticity regions are not heavy enough to move away, while the effect of micro-mixing is not strong enough to disperse droplets entrapped into high vorticity areas ~\cite{reveillon2007effects}. 
The Fig. \ref{fig:6:b} provides the Stokes number, $St_\eta$, of droplets located within a radial distance $r/R=0.2$ from the jet centerline. In this region the small-scale clustering should be the most significant with respect to the preferential segregation effect induced by the intermittency of the external shear layer.  The Stokes number decreases in the downstream evolution of the spray assuming unity value around $z/R\simeq25$.
Hence it appears that droplets show an intense clustering promoted both by the small-scale turbulent clustering and by the entrainment process associated to the fluctuation of the interface which separates the inner core region rich of droplet and the outer region which is depleted of droplets.

\begin{figure}[t]
\subfigure[]{\includegraphics[width=0.45\columnwidth]{art/cli.png} \label{fig:6:a}} \
\subfigure[]{\includegraphics[width=0.50\columnwidth]{art/stokes_true.pdf} \label{fig:6:b}}
\caption{(a) Droplet clustering index, $K$. %The value of the clustering index is estimated averaging the number of droplets, $N$,  inside fixed control volume, $\Delta V$, defined on a Cartesian grid. 
(b) Evolution of the mean droplets Stokes number, $St_\eta=\tau_d/\tau_\eta$, based on the Kolmogorov dissipative scales on the jet axis. 
%Only droplets located within $r/R<0.2$ have been considered.
}
\label{fig:6}
\end{figure}

To quantify the importance of the droplet clustering in the evaporation process we compare the mean vapor concentration field \emph{felt} by the droplets, 
$Y_{V,dc}$, and the unconditioned Eulerian one, $Y_V$. $Y_{V,dc}$ is the vapor concentration field obtained by a conditional average on the droplet presence in a given point. Hence if droplets preferentially accumulate in locations where the vapor concentration is relatively high it results $Y_{V,dc}>Y_V$. 
Figure~\ref{fig:9} reports the radial profiles of  $Y_{V,dc}$ and $Y_V$ at different axial distances $z/R$. 

\begin{figure}[ht]
\subfigure[]{\includegraphics[width=0.48\columnwidth]{art/slice_z10.pdf} \label{fig:9:a}} 
\subfigure[]{\includegraphics[width=0.48\columnwidth]{art/slice_z20.pdf} \label{fig:9:b}}\\
\subfigure[]{\includegraphics[width=0.48\columnwidth]{art/slice_z30.pdf} \label{fig:9:c}} 
\subfigure[]{\includegraphics[width=0.48\columnwidth]{art/slice_z40.pdf} \label{fig:9:d}}
\caption{The figure report the radial profiles of the average vapor mass fraction field at four different axial distances from the origin: (a) $z/R=10$, (b) $z/R=20$, (c) $z/R=30$ and (d) $z/R=40$. Each plot shows the enstrophy-threshold conditional average, $Y_{v,ec}$, the droplet-presence conditional average, $Y_{v,dc}$ and the unconditional Eulerian one, $Y_v$. The enstrophy-threshold conditional average is calculated by sampling the vapor mass fraction only over turbulent core events ($I=1$), that is when local enstrophy exceeds a fixed threshold. $Y_{v,dc}$ is the vapor concentration field obtained by a conditional average on the droplet presence in a given point.}
\label{fig:9}
\end{figure}

The vapor mass fraction felt by droplets is usually higher then the correspondent unconditional value. At $z/R=10$ the droplet conditioned and unconditioned
vapor concentration are similar with the exception of the outer part. This behaviour is expected since we have observed the clustering to be small 
near the inflow with the exception of the mixing layer. The clustering associated to the mixing layer, separating the outer and inner jet regions, 
will be discussed in details in the following. At higher $z/R$ the preferential sampling of the vapor phase operated by segregating droplet is significant  
with an oversampling of about $10\div40\%$ more the unconditioned value even in the inner jet core. 
To characterize droplet dynamics in the mixing layer, we need to discern between the inner turbulent jet core and the irrotational outer region.
The two regions are separated by an almost sharp fluctuating layer, so-called turbulent/non-turbulent interface~\cite{da2014interfacial}, that is highly convoluted over a wide range of turbulent scales. The most used observable to characterize the two regions is 
 the local enstrophy, $\zeta^2=  || \nabla \times \vec u || ^2$.
The inner turbulent core is characterized by large fluctuations of enstrophy while in the outer region the enstrophy is null.  Thus, fixing an entrophy threshold, $\zeta^2_{th}$, it is possible to distinguish if a point is located into the turbulent region or not:

\begin{equation}
I(\vec x, t)=H[\zeta^2(\vec x, t) -\zeta^2_{th}]
\end{equation}
 
\noindent with $H$ the heaviside function. $I=1$ denotes a turbulent event, while $I=0$ an irrotational one. 
The value of $\zeta_{th}$ has been shown to play a 
weak influence~\cite{da2014interfacial}. Using $I$ we can define an enstrophy-threshold conditional average for the vapor concentration, 
$Y_{V,ec}$, by sampling the vapor mass fraction only over turbulent core events. Besides the unconditioned and the droplet conditioned statistics, 
Fig. \ref{fig:9} provides also $Y_{V,ec}$ which 
can be seen as the mean concentration field of the turbulent core region. In the mixing layer, this turbulent conditional average shows an excellent agreement
with the vapor concentration felt by the droplets which is significantly higher than the unconditioned value. Since the unconditioned value $Y_V$ is determined 
both by irrotational dry outer and turbulent vapor-concentrated events, the present analysis indicates that the droplet dynamics in the mixing layer is mainly 
determined by turbulent events. 
In other words, in the mixing layer, 
droplets moving towards the outer region from  the inner turbulent core are surrounded by highly concentrated vapor gas ejected with the droplets. On the contrary, in a point of the mixing
 layer, a low vapor concentration event is associated to an engulfment of entrained dry air which is depleted of droplets. 
 Because of this dynamics, in average, droplets
 evaporating in the mixing layer do not feel the unconditioned mean vapor concentration, but a higher level.     
At further downstream distance the effect
of small-scale clustering previously discussed adds its contribution to this dynamics. 


Hence, since droplet vaporization rate is driven by the vapor concentration sampled by droplets, the oversampling of the vapor concentration field slows down the overall vaporization process, thus increasing the overall vaporization length and time. We find that droplets preferential sampling is the results of two different contributions, originated by distinct mechanisms. The first mechanism is induced by the fluctuation of the turbulent/non-turbulent interface in the mixing layer. In this area droplets are entrapped in turbulent structures with high vapor concentration that protract into the irrotational, droplet-free ambient gas. 
%This is responsible for a part of the entire droplets preferential sampling which is much more evident in the shear layer: in this region the enstrophy and droplets conditional average are comparable but dominates over the unconditional vapor mass fraction average.
The second contribution is given 
by the inertial small scale clustering. Clusters of evaporating droplets move together with their own highly saturated atmosphere induce an oversapling of
the vapor concentration with respect to the neighboring environment. 
 The contribution of this latter mechanism is more evident in the spray core and tend to increase in intensity in the downstream evolution of the spray, while
 the former is dominant in the mixing layer.


%We measure the effect of preferential segregation on droplets vaporization dynamics, by the comparison between three different conditional and unconditional statistics for the vapor mass fraction field. The turbulent jet is constituted by a rotational vortical core which is bounded by the irrotational ambient fluid. The time-evolving interface between these two region is sharp and highly convoluted over a wide range of turbulent scale. The local enstrophy, or enstrophy density, can be defiend as:
%
%\begin{equation}
%\zeta= || \nabla \times \vec u || ^2
%\end{equation}
%We show that the values of the vapor mass fraction field samples by droplets population are significantly higher then their correspondent unconditional average values over the whole axial and radial extension of the turbulent spray. Since droplet vaporization rate is driven by the vapor concentration sampled by droplets, this mechanism slows down the overall vaporization process, thus increasing the overall vaporization length and time. We find that droplets preferential sampling is the results of three different contributions, originated by different preferential segregation mechanisms. The first mechanisms relies in droplets preferential segregation induced by the fluctuations of the turbulent/non-turbulent interface region. In this area droplets are entrapped in turbulent structure that protract into the irrotational, droplet-free ambient gas. This is responsible for a part of the entire droplets preferential sampling which is much more evident in the shear layer: in this region the enstrophy and droplets conditional average are comparable but dominates over the unconditional vapor mass fraction average.
%The second contribution is given by small isolated droplets moving together with their own highly saturated atmosphere, and thus being constantly surrounded by higher vapor concentration then that in the neighboring environment.  Finally, the last contribution is given by the small scale clustering. Clusters can be trailed by the gaseous carrier phase such that droplets inside them are maintained in a highly saturated environment. The contribution of these two latter mechanisms on the overall preferential sampling is more evident in the spray core, and tend to increase in intensity in the downstream evolution of the spray.

\subsection {Probability density function of droplet observables}
 
In order to further characterize droplet vaporization dynamics, we consider the probability density function of the vaporization length and time computed over the whole droplets population. In analogy with the overall vaporization length definition, one single droplet vaporization length can be defined as the axial distance from inlet, $z_e$, necessary for the the droplets radius to decrease from $r_{d,0}$ to a threshold radius $r_{d,th}=.01 r_{d,0}$. The vaporization time $t_e$ is the corresponding amount of time. The PDFs for droplets evaporation length and time are reported in Fig. \ref{fig:7:a} and \ref{fig:7:b} respectively. The mean, standard deviation, skewness and kurtosis are reported in table \ref{tab:1} and show a nearly Gaussian behavior with significant standard deviations. The Gaussian behavior in turbulent flows is usually associated to fluctuations induced by the large-scale motions.  It is remarkable how different are the histories of the droplets: half of the injected droplets is still present at about $z/R\simeq32$ where about $90\%$ of the injected liquid mass fraction is evaporated. This aspect is connected to the high polydispersity developed by the droplets.

\begin{table}[h]
\centering
\begin{tabular}{|l|c|c|c|c|}
\hline
                            & $\mu$  & $\sigma$ & $K$  & $S$  \\ \hline
\ \ $z_e/R$   \ \    & 30.64 & 10.22 & 2.39 & -0.010      \\ \hline
\ \ $t_e/t_0$  \ \   & 60.22 & 14.06 & 2.71 & -0.378      \\ \hline
\end{tabular}
\label{tab:1}
\caption{The table provides the mean, $\mu$, standard deviation, $\sigma$ , kurtosis, $K$ and skewness, $S$, of the PDFs of droplets evaporation length and time with $\mu=E[X]$,  $\sigma=\sqrt{E[(X-\mu)^2]}$, $K=E[(X-\mu)^4]/E[(X-\mu)^2]^2$ and $S=E[(X-\mu)^3]/E[(X-\mu)^2]^{3/2}$. All variables are non-dimensional.}
\end{table}

\begin{figure}[ht]
\centering
\subfigure[]{\includegraphics[width=0.47\columnwidth]{art/pdf_el2.pdf} \label{fig:7:a}} \
\subfigure[]{\includegraphics[width=0.47\columnwidth]{art/pdf_et2.pdf} \label{fig:7:b}} \
\caption{ (a) Probability density function of non-dimensional droplet vaporization length, $z_e/R$, with R the jet inlet radius. (b) Probability density function of non-dimensional droplet vaporization time, $t_e/t_0$, with $t_0$ the reference time scale, $t_0=R/U_b$. The PDFs are computed over the entire droplets population injected into the computational domain.}
\label{fig:7}
\end{figure}

\begin{figure*}[ht]
\subfigure[]{\includegraphics[width=0.45\textwidth]{art/pdf_r.pdf} \label{fig:8:a}} \
\subfigure[]{\includegraphics[width=0.45\textwidth]{art/pdf_s.pdf} \label{fig:8:b}} \\ 
\subfigure[]{\includegraphics[width=0.45\textwidth]{art/pdf_temp.pdf} \label{fig:8:c}} \ 
\subfigure[]{\includegraphics[width=0.45\textwidth]{art/pdf_evp.pdf} \label{fig:8:d}}
\caption{Probability Density Function (PDF) of Lagrangian variables. (a):  PDF of non-dimensional droplets radius, $r_d/r_{d,0}$, where $r_{d,0}$ is the initial radius of injected droplets. (b): PDF of the saturation field at droplets surface, $S_d=Y_{v,s}/Y_{v,d}$, where $Y_{v,s}=Y_{v,s}(T_d,p)$ is the vapor mass fraction at saturation computed as a function of droplet actual temperature and the carrier phase thermodynamics pressure, p. $Y_{v,ds}$ is the actual vapor mass fraction in the carrier gaseous mixture evaluated at droplet center position. (c): PDF of non-dimensional droplets temperature, $T_d/T_{d,0}$, where $T_{d,0}$ is the initial temperature of injected droplets. (d): PDF of non-dimensional droplets vaporization rate, $-\dot{m_d} \tau_{d,0} / m_{d,0}$, where $\tau_{d,0}$ and $m_{d,0}$ are the initial relaxation time and mass of injected droplets. The PDF plots (a), (b) and (c) are log-linear, while plot (d) is log-log.}
\label{fig:8}
\end{figure*}


Fig. \ref{fig:8:a} shows the probability density function of the droplet radius at different axial distances from the origin. Even starting from a monodisperse suspension, we suddenly observe a radius distribution which spans for around one decade after 10 jet radii from the inlet. It should be remarked that this quantity amounts in differences of droplet volumes of about $10^3$. This intense spread of the droplets size spectrum may be attributed to the complex preferential segregation dynamics which has been previously discussed. Indeed, the evolution of droplets, which is made up of aggregates of different size and dynamics, induces extremely different \emph{surrounding} conditions on droplets themselves as can be observed by the vapor saturation level \emph{felt} by the droplets (see Fig. \ref{fig:8:b}). Near the inlet at $z/R=0.25$, droplets sample the almost saturated vapor phase. From $z/R=10$ the PDF shows a wide range of sampled value caused by the complex droplet dynamics previously discussed. We also observe that few droplets show condensation as denoted by saturation values above 1.  Further downstream $z/R=20$ the spreading trend inverts and droplets are subjected to progressively more uniform saturation levels. The different evaporation dynamics caused by the saturation level felt by the particles induces a similar statistical behavior for the temperature PDF, see Fig. \ref{fig:8:c}.  
The high polydispersity combined to the wide spectrum of saturation levels sampled by the droplets induces a non-trivial behavior of the vaporization rate PDF. Even though we cannot provide arguments, we find that the PDFs of the vaporization rate, shown in Fig.~\ref{fig:8:d}, appears to follow a power-law with exponent about $-3$, independently by the axial distance from inlet. 

\section{Final remarks}
%The dynamics of a turbulent evaporating spray in an open environment is reproduced in a Direct Numerical Simulation framework in order to take advantage of the richness of details that this approach has the capabilities to capture. The 3D simulation reproduces a spray dynamics in dilute, non-reacting conditions accounting for the full coupling between the two phases due to mass, momentum and energy exchange. During this phase of the spray downstream evolution the most part of the injected liquid transit to vapor. Hence, the phenomenology that characterized the droplets dynamics in these conditions strongly affect the overall vaporization process and constitute the main objectives of study of this paper. Due to the prescribed dilute conditions, droplets mutual interactions are neglected and the point-droplet approximation is assumed. An hybrid Eulerian/Lagrangian approach is considered in order to describe the physics of the system. The evolution of the gaseous carrier phase is reproduced accounting for the effect of the density variations by a low Mach number formulation of the Navier-Stokes equations, while the dynamics of droplets is described by three Lagrangian equations for position, radius and temperature. Three sink/source coupling terms fully describe the mass, momentum and energy exchange between the two phases. The simulation is performed over a cylindrical domain discretized with a staggered grid of N \fp{numero!!!!} million nodes. The entrainment of external fluid is allowed by the traction-free boundary conditions prescribed at the side of the boundary. 

The dynamics of a turbulent evaporating spray is investigated by means of a Direct Numerical Simulation. The  simulation reproduces an acetone/air spray evolving in an open environment considering dilute, non-reacting conditions and accounting for the full coupling between the two phases due to mass, momentum and energy exchange. The entrainment of external dry air is also accounted. Liquid acetone monodisperse droplets are continuously injected within the turbulent gaseous phase at a bulk Reynolds number $Re_R=U R /\nu=6000$. A complete description of both instantaneous and average fields of Eulerian and Lagrangian observables is provided. The distribution of droplets is strongly inhomogeneous with clustering apparent. In particular, droplets seem to preferentially persist in high vapor concentration regions thus being affected by a reduction of the vaporization rate. The intensity of the preferential segregation is estimated by the evaluation of the clustering index. Preferential segregation develops downstream the jet inlet first in the mixing layer and then in the turbulent core. In particular, two different mechanisms driving the inhomogeneous droplet distribution are identified: inertial small scale clustering and droplet segregation induced by the turbulent/non-turbulent interface. The former one is the results of the competition between inertia and Stokes drag and is found to be responsible for droplets preferential accumulation mainly in the spray core and in the far field evolution of the flow. The latter one mainly affects the mixing layer and consists in the entrapment of droplets in  turbulent structures with high vapor concentration which are originated in the core  and protract towards the droplet-free dry environment. Simultaneously, droplet-free dry air regions are engulfed in the jet core enhancing the fluctuation of the droplet concentration. Both these mechanisms affect droplets dynamics and result in an oversampling of the vapor concentration experienced by each droplet, hence affecting the overall vaporization length. Probability Density Function of droplet observable have been reported at different axial distances from the inlet. A spectacular increase of the droplet polydispersity is found to arise in the downstream evolution of the spray resulting in an extreme widening of the droplets size spectrum. This intense spread is attributed to the heavy-tail PDF of the droplet vaporization rate which is the result of the complex dynamics coupling droplet and vapor concentration fields  
This mechanisms is expected to be important in all turbulent flows characterized by a mixing layer with entrainment of dry air, e.g. clouds.\\
The proper modeling of this phenomenon is critical in order to improve LES and RANS model capabilities to accurately reproduce the turbulent vaporization dynamics both for reacting and non-reacting sprays.
\section*{Acknowledgement}

The authors acknowledge financial support through the University of Padova Grant PRAT2015 (CPDA154914), as well as the computer resources provided 
by CINECA ISCRA C project: TaStE (HP10CCB69W).


%\bibliographystyle{aipnum-cp}
\bibliographystyle{plain}
\bibliography{biblio}

\end{document}

reftest4-1.tex
% ****** Start of file reftest.tex ******
%
%   This file is part of the APS files in the REVTeX 4.1 distribution.
%   Version 4.1 of REVTeX, July 1, 2008.
%
%   Copyright (c) 2009 The American Physical Society.
%
%   See the REVTeX 4.1 README file for restrictions and more information.
%
\documentclass{revtex4}%
\begin{document}%

\typeout{Filename: reftest4-1.tex for revtex 4.1i 2009/10/19 (AO)}
% This document, when typeset, prompts for the name of a LaTeX auxiliary (.aux) file 
% and then processes that file, paying attention to 
% keys that have been cited (via \citation)
% or defined (via \bibcite).
% It then prints out diagnostics indicating 
% * any keys that were cited but never defined,
% * any keys that were defined but never cited, and
% * keys that were cited and defined, but not in the same order.
%
% BibTeX takes care of such things automatically,
% so there is little use for such a test when using BibTeX.
%
% A separate question: will this file run under revtex4.cls?
% I can see two lines in my .aux with the potential for problems: 
% \global \mathchardef \firstnote@num\z@ \relax 
% \global\NAT@numberstrue

\makeatletter

\def\g@bbletw@#1#2{}%
\def\g@bble@ne#1{}%

\let\bibdata\g@bble@ne
\let\newlabel\g@bbletw@
\let\@writefile\g@bbletw@
\let\bibstyle\g@bble@ne

\@tempcnta\z@
\@tempcntb\z@

\def\bibcite#1#2{%
 \advance\@tempcnta\@ne
 \expandafter\bib@one\csname bib#1\expandafter\endcsname\expandafter{\the\@tempcnta}{#1}%
}

\def\bib@one#1#2#3{%
 \@ifnotrelax#1{%
  \typeout{*** Tag "#3" used in more than one \string\bibitem.}%
  \typeout{\space\space\space\space You must correct this problem before continuing (Bailing out).}%
  \@@end
 }{%
  \def#1{#2}%
  \expandafter\def\csname bib@num#2\endcsname{#3}%
 }%
}%

\def\citation#1{%
 \@for\@tempa:={#1}\do{%
  \expandafter\cite@one\csname cite\@tempa\expandafter\endcsname\expandafter{\@tempa}%
 }%
}

\def\cite@one#1#2{%
 \expandafter\@ifx\expandafter{\csname cite*\endcsname#1}{%
  \@booleantrue\@cite@star@sw
 }{%
  \@ifnotrelax#1{%
  }{%
   \advance\@tempcntb\@ne
   \expandafter\cite@one@\expandafter#1\expandafter{\the\@tempcntb}{#2}%
  }%
 }%
}%
\expandafter\def\csname cite*\endcsname{*}%
\@booleanfalse\@cite@star@sw

\def\cite@one@#1#2#3{%
 \def#1{#2}%
 \expandafter\def\csname cite@num#2\endcsname{#3}%
}%

\def\strut{\vrule height12pt depth4pt width0pt}%

\typein[\filename]{Type in file name (no extension)}%
\expandafter\input\expandafter{\filename.aux}%

\expandafter\def\expandafter\bibtotal\expandafter{\the\@tempcnta}%
\expandafter\def\expandafter\citetotal\expandafter{\the\@tempcntb}%

\def\app@undefined#1#2#3{%
 \expandafter\@ifnotrelax\csname#2#1\endcsname{}{%
   \appdef#3{#1, }%
  }
}%

\let\uncitedbibs\@empty
\@tempcnta\z@
\@ifnum{\bibtotal>\@tempcnta}{%
 \loopwhile{%
  \advance\@tempcnta\@ne
  \expandafter\app@undefined\csname bib@num\the\@tempcnta\endcsname{cite}\uncitedbibs
  \@ifnum{\bibtotal>\@tempcnta}%
 }%
}%

\let\undefinedcites\@empty
\@tempcnta\z@
\@ifnum{\citetotal>\@tempcnta}{%
 \loopwhile{%
  \advance\@tempcnta\@ne
  \expandafter\app@undefined\csname cite@num\the\@tempcnta\endcsname{bib}\undefinedcites
  \@ifnum{\citetotal>\@tempcnta}%
 }%
}%

\def\addit@bad#1#2{%
 \expandafter\@ifnotrelax\csname bib#1\endcsname{%
  \@ifx{#1#2}{%
   \addit{#1}{OK}%
  }{%
   \advance\@tempcntb\@ne
   \addit{#1}{\bf Out of order}%
  }%
 }{%
  \addit{#1}{\bf Undefined}%
 }%
}%
\def\ordertable{}%
\def\table@line#1#2{%
  \hbox to\hsize{%
   \vrule
   \hbox to.5\hsize{\ #1\hss}%
   \hss\vrule\hss\strut
   \hbox to.5\hsize{\ #2\hss}%
   \vrule
  }%
  \hrule
  \vskip\z@skip
}%
\def\addit#1#2{%
 \appdef\ordertable{\table@line{#1}{#2}}%
}%

\@tempcnta\z@
\@tempcntb\z@
\@ifnum{\citetotal>\@tempcnta}{%
 \loopwhile{%
  \advance\@tempcnta\@ne
  \expandafter\addit@bad\csname cite@num\the\@tempcnta\expandafter\endcsname\csname bib@num\the\@tempcnta\endcsname
  \@ifnum{\citetotal>\@tempcnta}%
 }%
}{}%
\expandafter\def\expandafter\citesdisordered\expandafter{\the\@tempcntb}%

Analysis of \filename.aux, according to rules of \emph{Physical Review}:
All bibitems must occur in the bibliography in order of their first citation in the
file, and all bibitems must be cited in the paper.

\def\@tempa{\textbf{Bibitems that were not explicitly cited:} }%
\@ifx@empty\uncitedbibs{%
 \appdef\@tempa{none.}%
}{%
 \appdef\@tempa{\uncitedbibs\ \@cite@star@sw{were implicitly cited.}{please cite or remove bibitem(s).}}%
}%
\@tempa\expandafter\typeout\expandafter{\@tempa}\par

\def\@tempa{\textbf{Citations not defined in bibliography:} }%
\@ifx@empty\undefinedcites{%
 \appdef\@tempa{none.}%
}{%
 \appdef\@tempa{\undefinedcites\space please add bibitem(s) or remove cite(s).}%
}%
\@tempa\expandafter\typeout\expandafter{\@tempa}\par

\@ifnum{\citesdisordered=\z@}{%
 \def\@tempa{\textbf{The order of the bibitems is OK.}}%
 \@tempa\expandafter\typeout\expandafter{\@tempa}\par
}{%
 \def\@tempa{%
  \textbf{There 
  \@ifnum{\citesdisordered<\tw@}{was an error}{were \citesdisordered\ errors}
  in the order of your bibitems.}%
 }%
 \@tempa\expandafter\typeout\expandafter{\@tempa}
 \@ifx@empty\uncitedbibs{\@ifx@empty\undefinedcites{\true@sw}{\false@sw}}{\false@sw}%
 {(This is in addition to the problems noted above.) }{}%
 Please check the following table to see what the problems were.
 Status in boldface signifies that a correction is needed. \par
}%
 \emph{NOTE:
  This table is a complete listing of all your citations, in order of their first occurrence: 
  your bibitems should be in the order shown in the first column of this table.
 }\par
 \small
 \smallskip
 \hrule
 \table@line{\textbf{Citation key}}{\textbf{Status}}%
 \smallskip
 \hrule
 \ordertable
\end{document}
% ****** End of file reftest.tex ******

paper.bbl
\begin{thebibliography}{10}

\bibitem{abramzon1989droplet}
B~Abramzon and WA~Sirignano.
\newblock Droplet vaporization model for spray combustion calculations.
\newblock {\em International journal of heat and mass transfer},
  32(9):1605--1618, 1989.

\bibitem{armenio2001importance}
Vincenzo Armenio and Virgilio Fiorotto.
\newblock The importance of the forces acting on particles in turbulent flows.
\newblock {\em Physics of Fluids}, 13(8):2437--2440, 2001.

\bibitem{attili2014formation}
Antonio Attili, Fabrizio Bisetti, Michael~E Mueller, and Heinz Pitsch.
\newblock Formation, growth, and transport of soot in a three-dimensional
  turbulent non-premixed jet flame.
\newblock {\em Combustion and Flame}, 161(7):1849--1865, 2014.

\bibitem{battista2011intermittent}
F~Battista, F~Picano, G~Troiani, and Carlo~Massimo Casciola.
\newblock Intermittent features of inertial particle distributions in turbulent
  premixed flames.
\newblock {\em Physics of Fluids}, 23(12):123304, 2011.

\bibitem{bec2007heavy}
Jeremie Bec, Luca Biferale, Massimo Cencini, Alessandra Lanotte, Stefano
  Musacchio, and Federico Toschi.
\newblock Heavy particle concentration in turbulence at dissipative and
  inertial scales.
\newblock {\em Physical review letters}, 98(8):084502, 2007.

\bibitem{bukhvostova2014comparison}
A~Bukhvostova, E~Russo, JGM Kuerten, and BJ~Geurts.
\newblock Comparison of dns of compressible and incompressible turbulent
  droplet-laden heated channel flow with phase transition.
\newblock {\em International journal of multiphase flow}, 63:68--81, 2014.

\bibitem{calzavarini2008quantifying}
Enrico Calzavarini, Massimo Cencini, Detlef Lohse, and Federico Toschi.
\newblock Quantifying turbulence-induced segregation of inertial particles.
\newblock {\em Physical review letters}, 101(8):084504, 2008.

\bibitem{da2014interfacial}
Carlos~B da~Silva, Julian~CR Hunt, Ian Eames, and Jerry Westerweel.
\newblock Interfacial layers between regions of different turbulence intensity.
\newblock {\em Annual review of fluid mechanics}, 46:567--590, 2014.

\bibitem{eaton1994preferential}
John~K Eaton and JR~Fessler.
\newblock Preferential concentration of particles by turbulence.
\newblock {\em International Journal of Multiphase Flow}, 20:169--209, 1994.

\bibitem{elghobashi1994predicting}
S~Elghobashi.
\newblock On predicting particle-laden turbulent flows.
\newblock {\em Applied scientific research}, 52(4):309--329, 1994.

\bibitem{falkovich2002acceleration}
G~Falkovich, A~Fouxon, and MG~Stepanov.
\newblock Acceleration of rain initiation by cloud turbulence.
\newblock {\em Nature}, 419(6903):151--154, 2002.

\bibitem{ferrante2003physical}
Antonino Ferrante and Said Elghobashi.
\newblock On the physical mechanisms of two-way coupling in particle-laden
  isotropic turbulence.
\newblock {\em Physics of fluids}, 15(2):315--329, 2003.

\bibitem{godsave1953}
GAE Godsave.
\newblock Studies of the combustion of drops in a fuel spray?the burning of
  single drops of fuel.
\newblock In {\em Symposium (International) on Combustion}, volume~4, pages
  818--830. Elsevier, 1953.

\bibitem{gualtieri2009anisotropic}
Paolo Gualtieri, F~Picano, and Carlo~Massimo Casciola.
\newblock Anisotropic clustering of inertial particles in homogeneous shear
  flow.
\newblock {\em Journal of Fluid Mechanics}, 629:25--39, 2009.

\bibitem{gualtieri2015exact}
Paolo Gualtieri, F~Picano, Gaetano Sardina, and Carlo~Massimo Casciola.
\newblock Exact regularized point particle method for multiphase flows in the
  two-way coupling regime.
\newblock {\em Journal of Fluid Mechanics}, 773:520--561, 2015.

\bibitem{jenny2012}
Patrick Jenny, Dirk Roekaerts, and Nijso Beishuizen.
\newblock Modeling of turbulent dilute spray combustion.
\newblock {\em Progress in Energy and Combustion Science}, 38(6):846--887,
  2012.

\bibitem{kennedy1997models}
Ian~M Kennedy.
\newblock Models of soot formation and oxidation.
\newblock {\em Progress in Energy and Combustion Science}, 23(2):95--132, 1997.

\bibitem{kinzer1951}
Gilbert~D Kinzer and Ross Gunn.
\newblock The evaporation, temperature and thermal relaxation-time of freely
  falling waterdrops.
\newblock {\em Journal of Meteorology}, 8(2):71--83, 1951.

\bibitem{kruger2016effect}
Jonas Kr{\"u}ger, Nils~EL Haugen, Dhrubaditya Mitra, and Terese L{\o}v{\aa}s.
\newblock The effect of turbulent clustering on particle reactivity.
\newblock {\em Proceedings of the Combustion Institute}, 2016.

\bibitem{lau2016effect}
Timothy~CW Lau and Graham~J Nathan.
\newblock The effect of stokes number on particle velocity and concentration
  distributions in a well-characterised, turbulent, co-flowing two-phase jet.
\newblock {\em Journal of Fluid Mechanics}, 809:72--110, 2016.

\bibitem{majda1985derivation}
ANDFREW MAJDA and James Sethian.
\newblock The derivation and numerical solution of the equations for zero mach
  number combustion.
\newblock {\em Combustion science and technology}, 42(3-4):185--205, 1985.

\bibitem{marchioli2008some}
C~Marchioli, MV~Salvetti, and A~Soldati.
\newblock Some issues concerning large-eddy simulation of inertial particle
  dispersion in turbulent bounded flows.
\newblock {\em Physics of Fluids}, 20(4):040603, 2008.

\bibitem{marmottant2004spray}
Philippe Marmottant and Emmanuel Villermaux.
\newblock On spray formation.
\newblock {\em Journal of fluid mechanics}, 498:73--111, 2004.

\bibitem{mashayek1998direct}
F~Mashayek.
\newblock Direct numerical simulations of evaporating droplet dispersion in
  forced low mach number turbulence.
\newblock {\em International journal of heat and mass transfer},
  41(17):2601--2617, 1998.

\bibitem{maxey1983equation}
Martin~R Maxey and James~J Riley.
\newblock Equation of motion for a small rigid sphere in a nonuniform flow.
\newblock {\em The Physics of Fluids}, 26(4):883--889, 1983.

\bibitem{mellado2017cloud}
Juan~Pedro Mellado.
\newblock Cloud-top entrainment in stratocumulus clouds.
\newblock {\em Annual Review of Fluid Mechanics}, 49:145--169, 2017.

\bibitem{miller1999direct}
Richard~S Miller and J~Bellan.
\newblock Direct numerical simulation of a confined three-dimensional gas
  mixing layer with one evaporating hydrocarbon-droplet-laden stream.
\newblock {\em Journal of Fluid Mechanics}, 384:293--338, 1999.

\bibitem{olivieri2014effect}
Stefano Olivieri, Francesco Picano, Gaetano Sardina, Daniele Iudicone, and Luca
  Brandt.
\newblock The effect of the basset history force on particle clustering in
  homogeneous and isotropic turbulence.
\newblock {\em Physics of fluids}, 26(4):041704, 2014.

\bibitem{picano2011dynamics}
F~Picano, F~Battista, G~Troiani, and Carlo~Massimo Casciola.
\newblock Dynamics of piv seeding particles in turbulent premixed flames.
\newblock {\em Experiments in Fluids}, 50(1):75--88, 2011.

\bibitem{picano2010anomalous}
F~Picano, G~Sardina, Paolo Gualtieri, and Carlo~Massimo Casciola.
\newblock Anomalous memory effects on transport of inertial particles in
  turbulent jets.
\newblock {\em Physics of Fluids}, 22(5):051705, 2010.

\bibitem{raman2016modeling}
Venkat Raman and Rodney~O Fox.
\newblock Modeling of fine-particle formation in turbulent flames.
\newblock {\em Annual Review of Fluid Mechanics}, 48:159--190, 2016.

\bibitem{reveillon2007evaporating}
J~Reveillon and FX~Demoulin.
\newblock Evaporating droplets in turbulent reacting flows.
\newblock {\em Proceedings of the Combustion Institute}, 31(2):2319--2326,
  2007.

\bibitem{reveillon2007effects}
Julien Reveillon and Fran{\c{c}}ois-Xavier Demoulin.
\newblock Effects of the preferential segregation of droplets on evaporation
  and turbulent mixing.
\newblock {\em Journal of Fluid Mechanics}, 583:273--302, 2007.

\bibitem{rocco2015curvature}
G~Rocco, F~Battista, F~Picano, G~Troiani, and Carlo~Massimo Casciola.
\newblock Curvature effects in turbulent premixed flames of h2/air: a dns study
  with reduced chemistry.
\newblock {\em Flow, Turbulence and Combustion}, 94(2):359--379, 2015.

\bibitem{sardina2015continuous}
Gaetano Sardina, Francesco Picano, Luca Brandt, and Rodrigo Caballero.
\newblock Continuous growth of droplet size variance due to condensation in
  turbulent clouds.
\newblock {\em Physical review letters}, 115(18):184501, 2015.

\bibitem{saw2008inertial}
Ewe~Wei Saw, Raymond~A Shaw, Sathyanarayana Ayyalasomayajula, Patrick~Y Chuang,
  and {\'A}rmann Gylfason.
\newblock Inertial clustering of particles in high-reynolds-number turbulence.
\newblock {\em Physical review letters}, 100(21):214501, 2008.

\bibitem{shaw2002towards}
RA~Shaw, AB~Kostinski, and ML~Larsen.
\newblock Towards quantifying droplet clustering in clouds.
\newblock {\em Quarterly journal of the royal meteorological society},
  128(582):1043--1057, 2002.

\bibitem{spalding1953}
Dudley~Brian Spalding.
\newblock The combustion of liquid fuels.
\newblock In {\em Symposium (international) on combustion}, volume~4, pages
  847--864. Elsevier, 1953.

\bibitem{toschi2009lagrangian}
Federico Toschi and Eberhard Bodenschatz.
\newblock Lagrangian properties of particles in turbulence.
\newblock {\em Annual Review of Fluid Mechanics}, 41:375--404, 2009.

\end{thebibliography}
