Direct numerical simulation of spray droplet evaporation in hot turbulent channel flow

We perform a direct numerical simulation (DNS) of 14081 “cold” spherical droplets evaporating in a “hot” fully-developed turbulent channel flow. This effort is the first extensive computation that employs fourway coupling of the droplet motion with the turbulent carrier phase and interface-resolved evaporation dynamics, for a flow configuration that approaches conditions encountered in spray combustion applications. The complex interaction of momentum, heat, species transfer and phase change thermodynamics is explored. Large-scale droplet motion, modulation of the carrier phase turbulence, and influence of the mean and turbulent mass transport on the evaporation dynamics are observed and quantified. Based on the data set, phenomenological explanations of the shear-induced migration of the dispersed phase and of the effect of turbulent mass transport on the evaporation are provided. The transient nature of the DNS is exploited to generate a novel database that samples a range of turbulence and evaporation timescales, from which a model for the enhancement of the evaporation rate by the ambient turbulence is extracted. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
In recent years, the growth of computational power is being exploited in scientific computing, and exascale simulations are on track to be the next major breakthrough in the field. As a consequence, it is now affordable in computational fluid dynamics (CFD) to employ computationally intensive direct numerical simulations (DNS) on certain complex problems that were previously only treatable with coarse-grained simplified models. Spray evaporation in turbulent environment fits this category, being a multicomponent and multiphase flow with complex interaction of mass, momentum, energy transport and phase change thermodynamics.
Evaporation of liquid fuel spray is common in energy and propulsion systems, with a decisive impact on flame stabilization, combustion efficiency and emissions of pollutants. With the rise of alternative fuels for power generation and transport applications, new challenges are encountered in the design and engineering of stable, reliable, and low-emission fuel injectors and burners [1,2] .
The very limited availability and the high cost, per batch produced, of novel alternative fuels pose serious challenges in conventional approaches to the development of industrial combustion * Corresponding author.
E-mail address: gianlupo@mech.kth.se (G. Lupo). devices that often require relatively large amounts of fuel for fullscale testing. In this context, reliable modelling tools would be greatly beneficial in the development of such devices, but a comprehensive understanding of the physics of the fundamental underlying processes that govern fuel evaporation (and combustion) in turbulent flows is still lacking. To date, numerical studies in the field have mostly been limited to the Large Eddy Simulation-Lagrangian Particle Tracking (LES-LPT) approach, which relies on parametrized closures for the turbulence, the inter-phase interactions, and the phase change. More fundamental investigations, like DNS, have been rare and also limited in one or more aspects.
A review of recent efforts in the field of DNS of turbulent flows with droplets has been carried out by Elghobashi [3] . Therein, the works of Miller & Bellan [4] , Le Clercq & Bellan [5] , Russo et al. [6] , Kuerten & Vreman [7] are concerned with spray evaporation. Other studies in the same vein are those by Reveillon & Demoulin [8] , Weiss et al. [9] , Dalla Barba & Picano [10] . Several important insights on the interaction of turbulence, dispersed phase motion and mass transfer are gained from these works. However, the works cited all rely on a material point description (LPT) for the evaporating droplets, using closure models for the phase change.  Table 1 Underlying assumptions of the physical model.

A1
All physical and transport properties are constant. A2 The flow is incompressible. A3 Gravity is neglected. A4 Droplets remain spherical. A5 The fluid motion inside the droplets is neglected. A6 Temperature is uniform on the droplet surface. A7 The gas phase is ideal. A8 The inert gas is insoluble in the liquid phase. A9 Thermodynamic equilibrium prevails at the droplet surface. A10 The surface tension effect on vapor pressure (Kelvin effect) is neglected. A11 Viscous dissipation of energy is neglected. A12 Soret and Dufour effects are neglected.
change in a turbulent flow are limited to a single or a few droplets [11] .
In this work, we venture for the first time into DNS of spray evaporation in turbulent flow with four-way coupling of the carrier and dispersed phases, and interface resolved phase change at the droplet level, under conditions that approach those encountered in spray combustion applications, using the method proposed by Lupo et al. [12] . This allows for the direct solution of 14081 "cold" droplets evaporating in a "hot" turbulent channel flow, by a combination of physical assumptions (droplet sphericity) and numerical techniques (immersed boundary treatment of the gas-liquid interface). The results presented here provide unprecedented insight about the influence of large-scale motion of the dispersed phase and turbulent mass transport on the evaporation process.
An advantageous outcome of direct numerical simulations is the generation of databases that can be used to parametrize phenomenological models, which are employed in coarser modelling approaches and less demanding computations. In the words of Elghobashi: "Since large-eddy simulation will be used for the foreseeable future to predict turbulent multiphase flows at practical Reynolds numbers, accurate subgrid scale (SGS) models need to be developed and validated by DNS results [... ]. Accurate SGS models do not presently exist" [3] . We contribute to this aim by extracting from our DNS data a model for the enhancing effect of ambient turbulence on the droplet evaporation rate.
The paper is structured as follows: Section 2 is devoted to the statement of the physical problem and governing equations, Section 3 to a brief presentation of the numerical algorithm, Section 4 to the description of the DNS flow configuration and parameters. We analyse the results in Section 5 , focusing on droplet motion, turbulence modulation and evaporation dynamics. Section 6 presents a summary of the main findings and an outlook to future work.

Governing equations
The evaporation of a dilute spray is a multiphase and multicomponent flow characterized by a dispersed phase (the liquid droplets) and a carrier phase (the gas mixture). Momentum and energy are exchanged between the two phases, as well as transported in each phase separately. The mass of the vaporizing chemical species is exchanged between the two phases owing to the phase change, and transported through the carrier phase, where it mixes with the inert gas. No species transport occurs in the liquid phase when the droplet composition is pure. In order to simplify the treatment of the problem, we rely on some assumptions, which are listed in Table 1 . The reader is referred to [12] for a detailed discussion of the range of validity of each assumption.
Given these assumptions, the conservation of global mass, momentum, energy and mass of the vaporizing chemical species in the carrier phase can be expressed by the continuity, Navier-Stokes, temperature, and vapour mass fraction equations, which in their non-dimensional form read: where u is the carrier phase velocity, p its mechanical pressure, T its temperature, and Y the mass fraction of the vapour species in the carrier phase. The parameters Re, Pr, Sc are the Reynolds, Prandtl and Schmidt number respectively. The cross-transport term in Eq. (3) represents the net enthalpy transport by species diffusion, which occurs due to the different heat capacity of the vapour and inert gas species, specified by the parameter φ c p .
The a priori knowledge of the droplet shape (assumption A4) allows to write the mass and energy conservation of each individual liquid droplet in the form of global balance equations, which in their non-dimensional form read: Here, r d is the droplet radius and T s is the droplet surface temperature, considered uniform (assumption A6). The rates of heat and vapour mass exchange between the droplet and the carrier phase are ˙ q d and ˙ m d , respectively. The parameters φ ρ and φ c p represent the density and specific heat capacity ratio of the liquid to the gas mixture. The Stefan number Ste is a non-dimensional function of the latent heat of vaporization λ l .
The functions g 1 and g 2 are corrections that account for the fact that temperature is not uniform inside the droplet, and depend on the heat diffusivity ratio of the liquid to the gas mixture φ α : Their effect is only significant for big droplets, and vanishes (i.e. g 1 → 0 and g 2 → 1) as the droplet Biot number tends to zero. More details can be found in [12] . A summary of the non-dimensional parameters that appear in the governing equations is provided in Table 2 .
Since the fluid motion inside the droplet is neglected (assumption A5), the droplet motion consists of rigid body translation and rotation. The momentum exchange between the droplet and the carrier phase is thus described by the Newton-Euler equations, which in their non-dimensional form read: 8 15 (10) where x d , u d , and ω d are the position, velocity, and angular velocity of the droplet centroid, respectively. The force and torque integrals on the right hand side are taken over the droplet surface S ( t ).
When the droplets come close to each other, or to the wall, a lubrication force, normal to the surface, is activated. The force accounts for the gas film drainage in the gap between the two surfaces, and is based on the asymptotic expansion of the Stokesian analytical solution, as given by Jeffrey [13] . We use an activation threshold of d / r d ≤ 0.025 for droplet-droplet interactions and d / r d ≤ 0.05 for droplet-wall interactions, where d is the gap distance between the two droplets or between the droplet and the wall, and r d is the larger droplet radius. When d = 0 , the lubrication force is switched off, and a soft-sphere collision model is activated. In this model, the normal and tangential collision forces are calculated independently using simple spring-dashpot systems, with the addition of a Coulomb friction for the tangential force. More details on the lubrication and soft-sphere collision models and their parameters can be found in [14,15] .
The boundary conditions for the carrier phase at the droplet Here, P tot is the total thermodynamic pressure of the system, and P sat is the saturation pressure of the vaporizing chemical species, evaluated at the droplet surface temperature T s .
The heat and mass transfer rates are specified by integrating the heat and species fluxes over the droplet surface S ( t ) in the carrier phase:

Numerical method
The governing equations of the carrier phase Eqs. 1 to (4) are discretized in space with second order central finite differences, on a uniform ( x = y = z) staggered Cartesian grid (Eulerian mesh). Time integration is performed with a three-step Runge-Kutta scheme, both for the carrier phase PDEs, and for the dispersed phase ODEs Eqs. 5 to (10) . Within the Runge-Kutta scheme, a pressure correction scheme is used for the Navier-Stokes equation.
The calculation of the right hand side integrals of the Newton-Euler equations Eqs. 9 and (10) , the enforcement of the carrier phase boundary conditions at the droplet interface Eqs. 11 to (13) , and the calculation of the heat and mass transfer rates Eqs. 14 and (15) are all performed on the Lagrangian mesh points that represent the droplet surface, which is treated as an immersed boundary [16,17] . The Lagrangian points move rigidly together with the droplet centroid, and do not conform to the Eulerian mesh cells, although their spatial resolution matches approximately that of the Eulerian mesh. Interpolation from the Eulerian to the Lagrangian mesh and spreading from the Lagrangian to the Eulerian mesh are performed by means of the regularized Dirac delta function δ d introduced by Roma et al. [18] .
The presence of the phase change implies that the droplet surface injects mass and energy into the carrier phase. Therefore the following source terms have been added in order to mimic mass, energy, and vapour species injections, consistently with the interface boundary conditions: where Table 2 ) is the specific heat capacity ratio of the vapour species to the gas mixture. The specific form of the source terms, shown in the right-hand-sides of Eqs. 16 -18 , ensures that the mass, energy and species injections are distributed inside the droplet volume and vanish in the gas phase, going smoothly to zero for r → r d , and that their integral is consistent with the phase change mass and energy fluxes. A detailed description of the immersed boundary implementation and of the solution algorithm for the governing equations can be found in [12] .

Flow configuration
The numerical setup reproduces a turbulent channel flow of nheptane spray in air. The temperature and pressure conditions approach those commonly encountered in internal combustion engines, and the physical and transport properties have been calculated from said temperature and pressure using the averaging rule known as "1/3 rule" [19] , commonly employed in constant properties evaporation models. The channel walls are adiabatic, with homogeneous Neumann boundary conditions for T and Y , and the domain is periodic in the streamwise and spanwise directions. The streamwise flow is maintained at constant bulk Reynolds number by a body force in the momentum equation that linearly corrects the streamwise bulk velocity at every time step. In order to vent out of the periodic domain the additional mass injected by the evaporating droplets (through Eq. 16 ), a non-zero normal velocity is prescribed at the walls. Its value is equal to the integrated volume source from all the droplets, divided by the wall surface, reaching a peak of 2% of the mean streamwise velocity when the evaporation is strongest. The tangential velocities are set to zero at the walls. Table 3 shows the details of the computational setup. The values of the non-dimensional parameters that characterize the momentum, species, and energy balances are shown in Table 2 .
The bulk Reynolds number of 5600 corresponds to a friction Reynolds number, based on half the channel width, of Re τ = 180 , for a statistically stationary single phase channel flow.
The Eulerian grid resolution and number of Lagrangian points on the surface of each droplet are chosen in order to guarantee that the droplet initial diameter is resolved by 24 computational cells. This resolution ensures that both the turbulence in the gas phase, whose smallest scale is the Kolmogorov scale η, and the mass transport at the droplet interface, whose smallest scale is the Sc , are fully resolved. The viscous scale at the wall is also resolved ( y + ≈ 0 . 6 ).
The number of droplets is chosen in order to have an initial liquid volume loading of φ 0 = 0 . 05 .
The pressure and temperature conditions are typical of internal combustion engines. Under these conditions, a droplet diameter of 80 μm ensures that the average Weber number of the droplets is below unity, so that droplet deformation and breakup is negligible, and assumption A4 is reasonable.
The flow of the carrier phase is first initialized without the droplets. Following Henningson & Kim [20] , a vortex pair is superimposed on a plane laminar Poiseuille profile, and allowed to break down into a turbulent spot, as it is carried downstream by the mean flow. The flow eventually becomes fully turbulent and develops into a statistically stationary state, with Re τ = 180 .
At this point, the dispersed phase is introduced in the flow. The droplets are initialized in quasi-random positions with a Latin Hypercube algorithm, such that the liquid volume loading φ 0 is spatially uniform. It is found that the initial droplet velocity has no significant effect on the simulation results, owing to the relatively low density ratio φ ρ , and corresponding droplet Stokes number 1 Once the droplets are introduced, the evaporation is immediately activated by the vapour pressure build-up at the droplet 1 The droplet Stokes number is defined as surface, and is sustained by the heat provided by the hot gas flow. The simulation is carried on until the carrier phase gets saturated with vapour, whereupon the average droplet evaporation rate is significantly reduced, and condensation events start to become important. This corresponds to an average droplet diameter of d/d 0 = 0 . 97 , which ensures that the resolutions of the Eulerian and Lagrangian meshes are still approximately matching at the end of the calculation, a requirement for the accuracy of the immersed boundary method [21] . Fig. 1 shows three instantaneous snapshots of the flow. The background carrier phase is coloured by the vapour mass fraction, the droplets are coloured by their temperature. Migration of the majority of the droplets towards the channel centre, modulation of the flow turbulence, and saturation of the carrier phase with the vapour species during the course of the simulation can be observed in the visualizations, and will be quantified in the following. Based on the results, the transient evolution of each of these phenomena can be roughly divided into three stages. We choose the following sequence, based on the droplet evaporation dynamics, as a reference:

Results
where t end is the time at the end of the simulation.
Given the strong coupling between the physical phenomena involved in our case, it will be shown that this division is relevant for the description of the droplet migration and turbulence modulation as well.

Large-scale droplet motion
The dispersed phase is initially distributed uniformly in the channel, and is set to motion by the background carrier flow. Fig. 2 shows the local droplet number distribution, averaged in the two statistically homogeneous directions (streamwise x and spanwise z ), as a function of wall distance y and time t . A large-scale migration of the droplets towards the channel centreline is observed, and by the end of the simulation around 50% of the droplets have gathered in a 0.23 l y wide band across the channel centreline. The average droplet Stokes number, based on the local Kolmogorov time scale of the carrier phase flow, is found to be St ≈ 30, establishing the case within the regime described as four-way coupling by Elghobashi [22] . Four-way coupling of the dispersed phase motion with the carrier phase implies that the carrier phase mean flow, its turbulent fluctuations, the droplet-droplet and dropletwall collisions, and the local modification of flow streamlines by the droplet excluded volume, all contribute to determine each droplet's individual trajectory, as well as the large-scale migration of the dispersed phase.
The local liquid volume fraction φ is shown in Fig. 3 , as a function of the distance from the wall. The quantity is ensemble averaged in the two statistically homogeneous directions ( x and z ), and time averaged over Stages I, II and II. During Stage I, the initially uniform droplet distribution is perturbed by the background flow; Stage II represents a period of ongoing bulk migration of the droplets; during Stage III a quasi-steady distribution is finally attained. The final distribution is characterized by droplet clustering at the centreline, where the local volume fraction reaches the maximum value φ = 0 . 11 (at the end of Stage III), i.e. more than double its initial uniform value. The clustering tapers off gradually away from the centreline. A layer of droplets is also observed close to the wall. From the three profiles of Fig. 3 , it is observed that the formation of the wall layer is a faster process than the bulk migration towards the centreline.
The dispersed phase behaviour observed in the present case is similar to the one reported by Fornari et al. [23] , for a turbulent channel flow laden with solid particles, with conditions similar to the present case ( Re τ = 180 , φ 0 = 0 . 05 , φ ρ = 10 , L y /d 0 = 18 ). The final droplet distribution, with the coexistence of droplet clustering at the centreline and a droplet wall layer, is explained as the result of two concurrent mechanisms: shear-induced migration towards the channel centreline, and turbophoresis towards the channel walls.
Shear-induced migration is the net displacement of the dispersed phase from regions of high shear-rate towards regions of low shear-rate, i.e. away from the walls in a wall bounded flow. This is the result of irreversible short-range interactions between the droplets in the presence of shear. While classical literature  on shear-induced migration focuses on Stokes flow [24][25][26] , it is reasonable to expect that the phenomenon can occur at higher Reynolds numbers and also in turbulent conditions, provided that the Stokes number is high enough for the droplets to escape the carrier phase mean flow streamlines after a collision event. In fact, it is reported by Fornari et al. [23] that increasing centreline migration is observed for increasing values of the dispersed phase density (i.e. for increasing Stokes number), when the carrier phase turbulence and dispersed phase volume fraction are kept fixed. The work of Fornari et al. [23] deals with statistically stationary flow, whereas in the following we take advantage of the transient nature of our flow configuration, and present dynamical evidence in favour of the hypothesis that the Stokesian theory of shearinduced migration is applicable in turbulent conditions, for sufficiently large dispersed phase Stokes number.
The other observed mechanism, turbophoresis, is the migration of the dispersed phase towards regions of low turbulence intensity [27] . In a wall bounded flow these regions correspond to the vicinity of the wall. It is argued by Sardina et al. [28] that turbophoresis is a large-scale manifestation, in the presence of turbulence inhomogeneity, of small-scale clustering of the dispersed phase into regions of high strain and low vorticity, a more general phenomenon that is also observed in homogeneous turbulence [29][30][31] .
The time scale of shear-induced migration is governed by the rate of non-hydrodynamic droplet-droplet interactions, which is indirectly affected by the hydrodynamics of the carrier phase through the mean shear rate gradient [26] . Turbophoresis, on the other hand, is a purely hydrodynamic effect, whose time scale is governed by the correlation of the gradient of the velocity fluctuation [27] , i.e. it is related to the Taylor micro-scale [32] . Separation of time scales between the two phenomena, as observed in the present case, is therefore linked to the values of four flow parameters, namely the carrier phase Reynolds number, determining the ratio of droplet diameter to turbulent length scale, the liquid to gas density ratio φ ρ , determining the droplet Stokes number, the liquid volume fraction φ 0 , determining the droplet number density and thus the frequency of droplet-droplet interactions, and the ratio of channel width to droplet diameter L y / d 0 , determining the frequency of droplet interactions with the walls.
The dynamics of the large-scale migration are shown in Fig. 4 , in terms of the droplet wall normal mean velocity, averaged in the two homogeneous directions, and the root mean square of its  fluctuation, as functions of wall distance and time. A large number of droplets undergoes a velocity reversal, marked by the yellow dashed line in Fig. 4 (a): they drift towards the wall in the beginning, as a result of turbophoresis, but they turn towards the centreline after a major collision event takes place in the vicinity of the droplet wall layer and triggers shear-induced migration. In Fig. 5 , displaying eight successive instantaneous snapshots of the local liquid volume fraction profile along the wall normal coordinate, we mark the above mentioned collision event, and the subsequent "collision wave" that travels towards the channel centreline. The profiles of Fig. 5 have been filtered with a spatial window equal to (5 L y / N y ) for better clarity.
The speed of the collision wave is estimated as the slope of the yellow dashed line in Fig. 4 (a), as v wa v e / ( νRe/d 0 ) = 2 . 2 × 10 −3 . It is interesting to note that, according to the classical theory of shearinduced migration in Stokes flow [26] , the shear-induced migration velocity can be estimated as v migr ∼ −d 2 where U is the streamwise velocity of the carrier phase. For the present case, after calculating the carrier phase mean velocity as the ensemble average in the two homogeneous directions, we space-average its second order gradient in the wall normal direction across the channel width to find its mean value in the region affected by the migration, and perform a time average starting at t/t end = 0 . 11 , which is the approximate starting time of the collision wave, as seen in Fig. 5 . This gives an estimation of v migr / ( νRe/d 0 ) = 2 . 3 × 10 −3 , falling remarkably close to the calculated value of v wave . This finding substantiates our hypothesis that Stokesian shearinduced migration dynamics are still active in turbulent flow, provided that the dispersed phase Stokes number is sufficiently large.
The reduced wall normal velocity fluctuations at the centreline, visible in Fig. 4 (b), are an effect of the increased confinement of the dispersed phase in the centreline region, where the local volume fraction reaches its peak value, and is strongly correlated to the turbulence attenuation in the same region, as will be shown in Section 5.2 . The quenched wall normal velocity fluctuations ef- fectively trap the droplets that have reached the region around the centreline, stabilizing the central droplet cluster.
We shall see in Section 5.3 that the large-scale droplet motion in the channel, in particular the shear-induced migration towards the centreline, has a macroscopic effect on the evaporation dynamics.

Modulation of turbulence
In the present flow conditions ( Re = 5600 , φ 0 = 0 . 05 , φ ρ = 32 , L y /d 0 = 32 ), the ratio of droplet size to Kolmogorov length scale is approximately d 0 / η ≈ 20, and the ratio of droplet size to Taylor length scale is approximately d 0 / λ t ≈ 4, when the dispersed phase is released into the fully developed carrier phase turbulence. It is therefore expected for the dispersed phase to have a strong influence on the carrier phase turbulence, as the droplet size falls into the inertial subrange of turbulent length scales [31,33] .
In the following, turbulence statistics of the carrier phase are analysed. The statistics, conditioned to the carrier phase, have been calculated neglecting the regions occupied by the liquid phase, and ensemble averaging in the two homogeneous directions.
Profiles of the square root of the diagonal components of the turbulent Reynolds stress tensor (root mean square velocity fluctuations) and of the first off-diagonal component, along the wall nor-mal coordinate, are shown in Fig. 6 , normalized with the instantaneous friction velocity u τ . They are time averaged over Stages I, II and III, with Stage III divided into two sub-stages, to better illustrate the evolution of the Reynolds shear stress. The quantities are compared with the corresponding values from a single phase channel flow with Re τ = 180 (data from [34] ).
It is observed in Fig. 6 (a) that, when the droplets are released into the flow, and for an initial transient period that lasts approximately until the beginning of the droplet migration towards the centreline (Stage I), the carrier phase turbulence is greatly enhanced compared to the single phase case. In particular, the root mean square velocity fluctuations exhibit a flat profile along the channel width, aside from the near-wall region where they drop to zero. This reflects the roughly uniform droplet distribution in the channel, during this period. During Stage II, i.e. the main period of droplet migration, the carrier phase turbulence starts to decrease, especially in the region across the centreline, as a comparison between Fig. 6 (b) and Fig. 6 (a) shows. During Stage III, the bulk droplet migration slows down until the quasi-steady distribution of the dispersed phase shown in Fig. 3 is attained. At this stage, the flow is clearly divided into two regions with distinct turbulent features, as evident from Figs. 6 (c,d). When t = t end , the region around the centreline (0.15 ࣠ y / L y ࣠ 0.85), where around 93% of the droplets have gathered, has a much lower turbulent  Fig. 4 (b) shows that the diminished wall normal velocity fluctuation at the centreline is also strongly correlated to the small value of the wall normal droplet velocity fluctuation, hence the quasi-steady distribution of the dispersed phase is coupled to the turbulence attenuation. The value of the streamwise root mean square velocity fluctuation u rms in Fig. 6 (d) is still much larger than it is in the single phase case; however we interpret this result as the sampling by the carrier phase velocity of the many laminar-like droplet wakes present in this region, rather than as actual turbulent fluctuations of the carrier phase flow. This interpretation was first suggested by Parthasarathy & Faeth [35] , who pointed out that the contribution of mean streamwise velocity cannot be entirely separated from the contribution of turbulence in the particle wakes, since particle arrivals are random. Outside of the central region ( y / L y ࣠ 0.15 and y / L y ࣡ 0.85), only around 7% of the droplets remain at t = t end , which makes the local liquid volume fraction φ ≈ 0.001. Thus it is not surprising that the turbulence characteristics in this region are similar to those of the single phase channel, but squeezed in the wall normal coordinate, as the region occupies only ∼ 30% of the channel volume. In particular, as evident from Fig. 6 (d), v rms , w rms reach at y / L y ≈ 0.15 the value that they attain at the centreline in the single phase case, while u v has an inflection point that marks the boundary between the two flow regions at approximately the same location. The slightly augmented velocity fluctuations and Reynolds stress in the immediate vicinity of the walls ( y / L y ࣠ 0.05 and y / L y ࣡ 0.95), noticeable in Fig. 6 (d), can be explained by the excluded volume effect of the droplet centreline cluster, which squeezes part of the carrier phase flow in the near wall region, locally increasing the Reynolds number. These results are in line with the findings of Fornari et al. [23] for the case of turbulent channel flow laden with solid particles at similar condi- Asymptotically ( t > 0.67 t end ), it is found that the overall turbulence level is damped by the droplets. In fact, for a dispersed phase which is denser than the carrier phase, and with size larger than the Kolmogorov scale and smaller than the integral scale of turbulence, the higher specific inertia ( φ ρ > 1) and the enhanced dissipation in the inertial subrange arising from the dispersed phase drag act as sinks of the carrier phase turbulent kinetic energy [36] . The diminished turbulent transport has a considerable effect on the evaporation dynamics: the wall normal profiles of the turbulent vapour species wall normal flux, during Stages I, II, and III, are displayed in Fig. 8 , showing that the flux is reduced by one order of magnitude over the course of the simulation.

Evaporation dynamics
In the present configuration, the evaporation is active since the beginning of the droplets' life in the flow. Therefore the evaporation dynamics evolve simultaneously with the other phenomena, namely the droplet migration and development of quasi-steady droplet distribution across the channel, the increase and successive attenuation of the carrier phase turbulence intensity, and the progressive saturation of the carrier phase with vapour. The evolution of the individual droplet diameter and surface temperature is shown in Fig. 9 (a,b), together with the average over all the 14081 droplets (black line). The three distinct phases of the evaporation dynamics, which were identified as Stages I, II and II, can be characterized as follows: 1. Stage I, 0 ≤ t / t end ࣠ 0.11. The droplet temperature rises, first sharply and then more gently, as the droplets are heated by the surrounding gas, and the droplet average diameter decreases as a result of fast evaporation. We recall that the droplet centreline migration is triggered at t / t end ≈ 0.11, as shown in the fourth panel of Fig. 5 . 2. Stage II, 0.11 ࣠ t / t end ࣠ 0.33. The droplet temperature slowly decreases, as the imbalance between the sensible heat transferred to the droplets by the gas ( ˙ q d /φ c p < 0 in Eq. 6 ) and the latent heat of evaporation ( ˙ m d /Ste > 0 in Eq. 6 ) shifts towards the latter. As a consequence, the evaporation slows down and the droplet diameter decreases at a lower rate. 3. Stage III, 0.33 ࣠ t / t end ≤ 1. The mean droplet temperature and mean droplet diameter evolve towards an asymptotic value: a dynamic equilibrium is eventually established between evaporation and condensation events. During this period, the final droplet distribution of Fig. 3 is established, the carrier phase becomes saturated with vapour, and its turbulent features evolve to the final configuration of Fig. 6 (d).
The final distributions of the droplet diameter and droplet surface temperature are shown in Fig. 9 (c,d). The droplet surface temperature distribution ( Fig. 9 (d)) is bimodal, having two very close but distinct peaks. The negative skewness of both distributions reflects the asymptotic droplet distribution across the channel ( Fig. 3 ), suggesting that the droplet cluster around the centreline includes a group of droplets that have experienced stronger evaporation than the majority, and end up being smaller and colder. This is confirmed by the distribution of these tail droplets (i.e. those with d / d 0 < 0.96) along the wall normal coordinate, shown in Fig. 10 .
The strong coupling between the large-scale migration of the dispersed phase and the evaporation dynamics is evident from Fig. 11 , which shows the distribution of the droplet evaporation rate along the wall normal coordinate, time averaged over Stages I, II and III. The mean vapour mass fraction in the carrier phase, as a function of the wall normal coordinate, is shown on the side. The evaporation rate K = d d 2 d t is here defined as the dimensionless rate of change of the droplet surface area; it is negative for evaporation ( K < 0) and positive for condensation ( K > 0).
Before the migration starts, the evaporation rate is strongest, as already noted, and distributed equally across the channel ( Fig. 11 (a)). The vapour mass fraction in the carrier phase is also homogeneous and low, with a small dip close to the wall, corresponding to the narrow region depleted of droplets visible in the first profile of Fig. 3 ( y / L y ࣠ 0.05 and y / L y ࣡ 0.95).  As the evaporation advances ( Fig. 11 (b)), the vapour mass fraction builds up in the carrier phase; however it does so at a slower pace outside of the centreline region ( y / L y ࣠ 0.1 and y / L y ࣡ 0.9), owing to the large-scale migration towards the centreline. As a consequence, the evaporation rate of the droplets that linger outside of the central region stays higher the closer they are to the wall. During this stage, the droplet wall layer is found to have little to no influence on the evaporation rate distribution. The central region (0.1 ࣠ y / L y ࣠ 0.9) is characterized by a uniform distribution of the evaporation rate across its width, and a homogeneous vapour mass fraction in the carrier phase, despite the fact that the droplet distribution itself is not uniform ( Fig. 3 ). Condensa- tion events also start to appear where the liquid volume fraction is large. Similar considerations apply to Stage III ( Fig. 11 (c)), when a quasi-steady state is reached, with quasi-steady droplet distribution across the channel and dynamic equilibrium of evaporation and condensation events. The centreline cluster is more compact, thus its region of influence is now smaller (0.3 ࣠ y / L y ࣠ 0.7). A significant amount of condensation events takes place in this region: this is reflected in the noticeable local drop of the vapour mass fraction in the carrier phase. The region outside of the influence of the centreline cluster is now bigger ( y / L y ࣠ 0.3 and y / L y ࣡ 0.7): as before, the tail of the evaporation rate distribution bends towards stronger evaporation as it approaches the wall, but the influence of the droplet wall layer has become appreciable, and the evaporation rate distribution at the wall has now a large variance owing to the appearance of condensation events.
By showing side by side the time development of various mean quantities of the carrier and dispersed phases, Fig. 12 illustrates the influence of large-scale droplet motion, carrier phase turbulence and carrier phase mass transport on the evaporation dynamics. Times t/t end = 0 . 11 and t/t end = 0 . 33 , which delimit Stages I, II and III, are marked in the figure.
Stage I, characterized by droplet heating and fast evaporation, ends at t / t end ≈ 0.11, coinciding with the beginning of the large-scale centreline migration, as evidenced by the mean droplet velocity reversal of Fig. 12 (f). With regard to the interaction of turbulence and evaporation, this stage can be divided into two sub-stages, displaying some complex trends. In the begin-ning ( t / t end ࣠ 0.02), turbulence is increasing (growth of the mean Reynolds stress of Fig. 12 (a)), augmenting the heat transfer to the droplets and hence the droplet heating ( Fig. 9 (b)). Consequently, even though the vapour mass fraction Y starts to build up in the bulk carrier phase, the droplet temperature rise causes a faster increase of the vapour mass fraction at the droplet surface Y s ( Fig. 12 (c)). Since the difference ( Y s − Y ) is the driving force of the evaporation, a growth of the evaporation rate is observed during this initial sub-stage, and a peak value is reached at t / t end ≈ 0.02. Turbulence attenuation begins at t / t end ≈ 0.02 ( Fig. 12 (a)): the heat transfer from the gas to the liquid phase is reduced and the droplet heating stops ( Fig. 9 (b)), which halts the increase of vapour mass fraction at the droplet surface Y s ( Fig. 12 (c)), so that the evaporation rate starts to decrease ( Fig. 12 (d)). Interestingly, there is a lag between the average turbulent vapour species flux in the carrier phase and the evolution of turbulence: Fig. 12 (b) shows that, for 0.02 ࣠ t / t end ࣠ 0.11, v Y increases and then saturates, while u v is decreasing. This may be due to the fact that velocity fluctuations decay faster than mass fraction fluctuations ( Sc > 1).
During Stage II, 0.11 ࣠ t / t end ࣠ 0.33, the droplets accelerate towards the centreline ( Fig. 12 (f)). Simultaneously, the turbulent fluxes in the carrier phase decrease ( Fig. 12 (a,b)), the droplets cool down ( Fig. 9 (b)), and the vapour mass fraction at the droplet surface decreases while it keeps accumulating in the bulk ( Fig. 12 (c)). Thus, we observe a further decrease of the evaporation rate ( Fig. 12 (d)).
At t / t end ≈ 0.33, the large-scale migration reaches its peak velocity ( Fig. 12 (f)), and the droplets start to decelerate and attain a quasi-steady distribution ( Fig. 3 ). This marks the beginning of Stage III, when the carrier phase turbulent momentum flux decays to its asymptotic value ( Fig. 12 (a)), the turbulent vapour species flux is progressively extinguished ( Fig. 12 (b)), the carrier phase becomes saturated with vapour ( Fig. 12 (c)), and the evaporation rate further decreases to almost zero ( Fig. 12 (d)). Fig. 12 (e) shows the time evolution of the average vaporization Damköhler number of the dispersed phase. The concept of a vaporization Damköhler number was first introduced by Gökalp et al. [37] , in analogy with turbulent combustion, as a way to characterize the influence of turbulence on the vaporization dynamics. It is defined as the ratio of the characteristic time scale of turbulent eddies to the characteristic time scale of vaporization:

Turbulent scaling of the evaporation rate
The relevant turbulent time scale is that of the eddies having a length scale of the order of the instantaneous droplet diameter d . When the droplet size falls within the inertial subrange of turbulence, τ eddy can be estimated as: where ε is the local dissipation rate of turbulent kinetic energy.
The vaporization time scale is estimated as the advection time of the vapour across the mass transfer boundary layer on the droplet surface: The boundary layer thickness δ M is estimated from the film theory of Abramzon & Sirignano [38] : where  , spanning about two decades. Thus, a sizeable sample is available to suggest a meaningful correlation between turbulence and evaporation. Fig. 13 relates the mean evap- oration rate enhancement by turbulence, K/K l , defined by normalizing the evaporation rate with the evaporation rate under laminar conditions, to the mean vaporization Damköhler number Da v . The laminar evaporation rate is estimated as: Sh ReSc Data points for Stage I ( t / t end < 0.11) have been deemed unreliable with respect to the calculation of K l , owing to inaccurate evaluation of the Spalding mass transfer number B M , and are not reported in Fig. 13 . For t / t end ≥ 0.11, it is found by a least squares fit that the turbulent enhancement of the evaporation rate decreases with the Damköhler number as: This trend is remarkably similar to the correlation K/K l = 0 . 771 Da −0 . 111 v , found by Wu et al. [40] , and obtained after a series of evaporation experiments of a single liquid fuel droplet for various free-stream turbulence intensities and scales [40,41] . It is thus shown for the first time that a Damköhler power law holds even when the droplets have a strong feedback on the turbulence field, and in the presence of a denser spray for which the influence of neighbouring droplets on the evaporation dynamics cannot be neglected. Thus, Eq. 27 reveals a more fundamental mechanism of mass transfer than what was previously reported by Wu et al. [40,41] , whose experiment was limited to isolated droplets in an externally controlled turbulent stream.
It is the hope of the authors that the present result attracts new attention to the distinction between turbulence effects (Damköhler number) and mean flow effects (Sherwood number) on the evaporation dynamics, a fundamental point that is still neglected by most evaporation models.

Conclusions and outlook
This work presents the results of a direct numerical simulation of spray evaporation in turbulent channel flow. The 14081 spray droplets are interface resolved, and their motion is four-way coupled to the carrier phase. Thus, this effort represents one of the largest numerical studies of phase change flow to date, both in scope and in the detail of the physical phenomena described.
The conditions explored in the simulation lead to the emergence of several interrelated phenomena. The droplet size lies in the inertial subrange of the carrier phase turbulence, which leads to significant turbulence modulation. The droplet Stokes number is intermediate ( St ≈ 30), which allows the coexistence and interaction of hydrodynamic and non-hydrodynamic effects, such as turbophoresis and shear-induced migration. The choice of the standard turbulent channel flow represents a convenient and widely adopted idealized configuration. While of limited direct relevance for applications, it enables the simultaneous exploration of droplet migration, turbulence modulation and change in evaporation dynamics, thereby providing the fundamental physical insights required for the design and optimization of more complex configurations.
The dynamics of the shear-induced migration are investigated: the triggering event is identified in a large-scale collision close to the wall; the migration velocity is measured and found to be remarkably close to the value predicted by the classical theory for Stokes flows. We then suggest that Stokes flow migration dynamics can be applicable in turbulent conditions depending on the Stokes number of the dispersed phase. Further investigation is necessary in order to validate this claim; in particular a parametric study on the dispersed phase Stokes number, volume loading φ 0 , and size ratio L y / d 0 is desirable, despite the large computational cost.
While the evaporation rates are not strong enough to have a feedback on the carrier phase mean flow and turbulence, and on the dispersed phase motion, the reverse interactions are strong: turbulence modulation and large-scale droplet migration are found to have a crucial influence on the evaporation dynamics. The details of these interactions are analysed, and reveal the important role of correlated velocity and vapour mass fraction fluctuations, and of the local liquid loading and gas saturation.
The influence of the carrier phase turbulence on the droplet evaporation rate is quantified in the form of a correlation between the evaporation enhancement and the vaporization Damköhler number, which is found to follow a power law with exponent −0 . 1 . This result reveals a fundamental difference between turbulent and mean flow effects on mass transfer enhancement, which is shown to be significant even when the droplets have a strong feedback on the turbulence field, and in the presence of a dense spray for which the influence of neighbouring droplets on the evaporation dynamics cannot be neglected.
While the present study is the first one of its kind where a large number of spray droplets undergoing phase change in a turbulent flow is resolved at all scales, several simplifying assumptions have been made in order to keep the computational cost feasible. It is therefore desirable to determine which assumptions are the most critical and in which regimes they are acceptable. To this end, we are developing a Volume of Fluid numerical framework that incorporates low-Mach treatment of the continuous phase, with temperature and composition dependent density and transport coefficients, based on the method illustrated in [42] . The effects of density and transport coefficients variations, non-instantaneous relaxation of the vapour velocity to the surrounding gas mixture value, liquid phase circulation, non-rigid rotation, non-uniform temperature, deformation, break-up and coalescence, will be captured by the Volume of Fluid framework, in contrast to the present Immersed Boundary implementation, which forgoes these features in order to enable the computation of a statistically large sample of droplets.
Another limitation of the present flow configuration lies in the boundary conditions, which lead to saturation of the gas phase with the vaporizing species over the course of the simulation. While this grants that a wide range of evaporation rates, and therefore Damköhler numbers, is sampled, on the other hand it forces the evaporation to reach a dynamical equilibrium in a rather short time, and therefore keeps the droplet size distribution narrow. A configuration in which saturation does not occur will allow to explore the dynamical effect of size distribution on the physics of the problem.
Finally, an important future development will be the inclusion of multi-species thermodynamics, to simulate the evaporation of liquid mixtures.

Declaration of Competing Interest
We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.