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

ContentslistsavailableatScienceDirect

International Journal of Heat and Mass Transfer

journalhomepage:www.elsevier.com/locate/hmt

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

Giandomenico Lupo

^a,∗

, Andrea Gruber

^b,c

, Luca Brandt

^a,c

, Christophe Duwig

^a

aLinné FLOW Centre and SeRC (Swedish e-Science Research Centre), Department of Mechanics, Royal Insitute of Technology (KTH), Stockholm, Sweden

bSINTEF Energy Research, Thermal Energy Department, Trondheim, Norway

cDepartment of Energy and Process Engineering, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

a rt i c l e i n f o

Article history:

Received 9 March 2020 Revised 18 June 2020 Accepted 8 July 2020 Available online 19 July 2020 Keywords:

Spray Fuel Droplet

Turbulent multiphase flow Evaporation

Phase change

Direct numerical simulation

a b s t r a c t

Weperformadirectnumericalsimulation(DNS)of14081“cold” sphericaldropletsevaporatingina“hot”

fully-developedturbulentchannelflow.Thiseffortisthefirstextensivecomputationthatemploysfour- waycouplingofthedropletmotionwiththeturbulentcarrierphaseandinterface-resolvedevaporation dynamics,foraflowconfigurationthatapproachesconditionsencounteredinspraycombustionapplica- tions.Thecomplexinteractionofmomentum,heat,speciestransferandphasechangethermodynamics isexplored.Large-scaledropletmotion,modulationofthecarrierphaseturbulence,andinfluenceofthe meanandturbulentmasstransportontheevaporationdynamicsareobservedandquantified.Basedon thedataset,phenomenologicalexplanationsoftheshear-inducedmigrationofthedispersedphaseand oftheeffectofturbulentmasstransportontheevaporationareprovided.ThetransientnatureoftheDNS isexploitedtogenerateanoveldatabasethatsamplesarangeofturbulenceandevaporationtimescales, fromwhichamodelfortheenhancementoftheevaporationratebytheambientturbulenceisextracted.

© 2020TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

In recent years, the growth of computational power is being exploitedinscientificcomputing, andexascalesimulations areon track to be the next major breakthrough in the field. As a con- sequence, it is now affordable in computational fluid dynamics (CFD) to employ computationally intensive direct numericalsim- ulations (DNS)oncertaincomplexproblems thatwere previously onlytreatable withcoarse-grainedsimplifiedmodels.Sprayevap- orationinturbulentenvironmentfitsthiscategory, beingamulti- componentandmultiphaseflowwithcomplexinteractionofmass, momentum,energytransportandphasechangethermodynamics.

Evaporation of liquid fuel spray is common in energy and propulsionsystems,withadecisiveimpactonflame stabilization, combustionefficiencyandemissionsofpollutants.Withtheriseof alternative fuels for powergeneration andtransport applications, newchallengesareencountered inthedesign andengineeringof stable,reliable,andlow-emissionfuelinjectorsandburners[1,2].

The verylimitedavailability andthe highcost,per batch pro- duced, of novel alternative fuels pose serious challengesin con- ventionalapproachestothedevelopmentofindustrialcombustion

∗ Corresponding author.

E-mail address: [email protected] (G. Lupo).

devicesthatoftenrequirerelativelylargeamounts offuelforfull- scale testing. In this context, reliable modelling tools would be greatlybeneficialinthe developmentofsuch devices,butacom- prehensive understanding of the physics of the fundamentalun- derlying processesthat governfuelevaporation(andcombustion) in turbulent flows is still lacking. To date, numerical studies in thefield havemostlybeenlimitedto theLargeEddySimulation- Lagrangian Particle Tracking (LES-LPT) approach, which relies on parametrizedclosuresfor theturbulence,theinter-phase interac- tions,andthephasechange.Morefundamentalinvestigations,like DNS,havebeenrareandalsolimitedinoneormoreaspects.

AreviewofrecenteffortsinthefieldofDNSofturbulentflows withdropletshasbeencarriedout byElghobashi[3].Therein,the worksof Miller & Bellan [4], Le Clercq& Bellan [5], Russo et al.

[6], Kuerten&Vreman [7]are concerned withsprayevaporation.

OtherstudiesinthesameveinarethosebyReveillon&Demoulin [8], Weiss et al. [9], Dalla Barba & Picano [10]. Several impor- tantinsightsontheinteractionofturbulence,dispersedphasemo- tionandmasstransferaregainedfromtheseworks.However,the workscited all rely ona material point description (LPT)for the evaporating droplets, usingclosure models forthe phase change.

Therefore, the full implications of the interactions between the phase change andthe other flow phenomena are not taken into account. Earlier efforts aimed at direct solution ofdroplet phase https://doi.org/10.1016/j.ijheatmasstransfer.2020.120184

0017-9310/© 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/ )

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.

changeinaturbulentflowarelimitedtoasingleorafewdroplets [11].

In this work, we venturefor the first time into DNSof spray evaporationinturbulentflowwithfour-waycouplingofthecarrier anddispersed phases, andinterface resolved phasechangeatthe droplet level, under conditions that approach those encountered inspray combustionapplications, using themethod proposed by Lupoetal.[12].Thisallowsforthedirectsolutionof14081“cold”

dropletsevaporatingina“hot” turbulentchannel flow,by acom- binationofphysicalassumptions(droplet sphericity)andnumeri- caltechniques(immersedboundarytreatmentofthegas-liquidin- terface).Theresultspresentedhereprovideunprecedentedinsight aboutthe influence oflarge-scale motion ofthe dispersed phase andturbulentmasstransportontheevaporationprocess.

An advantageous outcome of direct numerical simulations is thegenerationofdatabasesthatcanbe usedtoparametrizephe- nomenological models,which are employed in coarser modelling approachesandlessdemandingcomputations.InthewordsofEl- ghobashi:“Since large-eddy simulationwill be usedforthe fore- seeable future to predict turbulent multiphase flows at practical Reynoldsnumbers,accuratesubgridscale(SGS)modelsneedtobe developedandvalidatedbyDNSresults[...].AccurateSGSmodels donotpresentlyexist”[3].Wecontributetothisaimbyextracting fromour DNSdata a model forthe enhancing effect of ambient turbulenceonthedropletevaporationrate.

The paper is structured as follows: Section 2 is devoted to the statement of the physical problem andgoverning 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 dropletmotion,turbulencemodulationandevaporationdynamics.

Section6presentsasummaryofthemainfindingsandanoutlook tofuturework.

2. Governingequations

The evaporation of a dilute spray is a multiphase and mul- ticomponentflow characterized by a dispersed phase (the liquid droplets)and a carrier phase (the gas mixture). Momentum and energyare exchanged betweenthe two phases, aswell as trans- portedineachphaseseparately.Themassofthevaporizingchem- ical species is exchanged between the two phases owing to the phasechange,andtransportedthroughthecarrierphase,whereit mixeswiththeinertgas.Nospeciestransportoccursintheliquid phasewhenthe dropletcomposition ispure. Inorderto simplify thetreatmentoftheproblem,werelyonsomeassumptions,which arelisted inTable 1.The readerisreferred to [12]fora detailed discussionoftherangeofvalidityofeachassumption.

Giventheseassumptions,theconservationofglobalmass,mo- mentum, energyand mass of the vaporizing chemical species in thecarrierphasecanbeexpressedbythecontinuity,Navier-Stokes, temperature, andvapour mass fraction equations,which in their

non-dimensionalformread:

∇

^{·u=0; (1)}

∂

^u

∂

^{t =−u·}

∇

^u−

∇

^p+_Re¹

∇

^{2u; (2)}

∂

^T

∂

^{t =−u·}

∇

^T+_ReP¹_r

∇

^2T+

φ

cp

ReSc

∇

^T·

∇

^{Y; (3)}

∂

^Y

∂

^{t =−u·}

∇

^Y+ 1

ReSc

∇

^2Y; (4)

where u is the carrierphase velocity, p its mechanical pressure, T its temperature, andY the massfraction ofthe vapour species in the carrier phase. The parameters Re, Pr, Sc are the Reynolds, PrandtlandSchmidtnumberrespectively.Thecross-transportterm inEq. (3)representsthenetenthalpytransport byspeciesdiffu- sion,whichoccursduetothedifferentheatcapacityofthevapour andinertgasspecies,specifiedbytheparameter

φ

cp.

The a priori knowledge of the droplet shape (assumption A4) allowstowritethemassandenergyconservationofeachindivid- ual liquiddropletin theform ofglobalbalance equations,which intheirnon-dimensionalformread:

dr_d

dt =− m˙_d

4

π

^r2d

φ

ρ; (5)

dT_s

dt =g₁(^rd,r˙_d,Re,Pr,

φ

α,Ts,t)−g⁻₂¹(^rd,Re,Pr,

φ

α,t)³

_q˙

φcdp +^m_Ste^˙d

4

π

^r3d

φ

ρ . (6) Here,r_d isthedropletradiusandT_sisthedropletsurfacetem- perature, considered uniform(assumption A6). The rates of heat and vapour mass exchange between the droplet and the carrier phaseareq˙_dandm˙_d,respectively.Theparameters

φ

ρ and

φ

cp rep- resentthedensityandspecificheatcapacityratiooftheliquidto thegasmixture.TheStefannumberSteisanon-dimensionalfunc- tionofthelatentheatofvaporization

λ

l.

Thefunctionsg₁andg₂arecorrectionsthataccountforthefact thattemperatureisnotuniforminsidethedroplet,anddependon theheatdiffusivityratiooftheliquidtothegasmixture

φ

α:

g1=

6

RePr

φ

α

(

^Ts,0−Ts

)

1 r²_d−_r²3

d drd

dtt

^∞

n=1

e−_RePr¹

_nπ

rd

2

φ

αt

g2 ; (7a)

g₂=1−^∞

n=1

6

(

ⁿ

π )

^2e−

1 RePr

_nπ

rd

2

φ

αt

. (7b)

Theireffectisonlysignificantforbigdroplets,andvanishes(i.e.

g₁→0andg₂→1)asthedropletBiotnumbertendstozero.More detailscanbefoundin[12].

Table 2

Non-dimensional flow parameters.

Parameter Name Definition Value

Re Reynolds number ^UL_ν^z 5600

Pr Prandtl number ^ν_α 0.739

Sc Schmidt number _D^ν_vap 2.07

φρ density ratio ^ρ_ρ^l 32

φcp heat capacity ratio ^c_c^pl_p 2.14 φα heat diffusivity ratio ^α_α^l 0.031 φcp,vap vapour to gas mixture ^cp,_c^v_p^ap 2.47

heat capacity ratio

φcp normalized vapour and ^cp,vap−c_cp^p,inert 1.4857 inert gas heat capacity

difference

Ste Stefan number ^cpl_λ^T∞

l 5.3135

Asummary ofthenon-dimensional parametersthat appearin thegoverningequationsisprovidedinTable2.

Sincethefluidmotioninsidethedropletisneglected(assump- tionA5),thedropletmotionconsistsofrigidbodytranslationand rotation. The momentum exchange between the droplet and the carrier phase is thus described by the Newton-Euler equations, whichintheirnon-dimensionalformread:

dxd

dt =u_d; (8)

4

3

πφ

ρr³_ddu_d dt =

S

−pI+ 1 Re

∇

^u+

∇

^uT

·ndS; (9)

8

15

πφ

ρr_d⁵d

ω

d

dt =

S

(

^rdn

)

^×

−pI+ 1 Re

∇

^u+

∇

^uT

·ndS; (10)

wherex_d,u_d,and

ω

daretheposition,velocity,andangularveloc- ityof the dropletcentroid, respectively.The force and torquein- tegrals onthe righthand side are takenoverthe 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,andis basedon the asymptoticexpansion oftheStoke- siananalyticalsolution,asgivenbyJeffrey[13].Weuseanactiva- tion thresholdof d/r_d≤0.025 fordroplet-dropletinteractions and d/r_d≤0.05 for droplet-wall interactions, where d is the gap dis- tance betweenthe two droplets orbetween the droplet andthe wall, andr_d is the larger droplet radius. When d=0, the lubri- cation force isswitched off, and a soft-spherecollision model is activated.Inthismodel,thenormalandtangentialcollisionforces arecalculatedindependentlyusingsimplespring-dashpotsystems, with the addition of a Coulomb friction for the tangential force.

More details on the lubrication and soft-sphere collision models andtheirparameterscanbefoundin[14,15].

The boundary conditions for the carrier phase at the droplet surface(

|

^x−x_d

|

=r_d)are:

u=ud+

( ω

d×n

)

^rd−

φ

ρ−1

drd

dt n; (11)

T=Ts; (12)

Y=Ys= P^sat

(

^Ts

)

Ptot . (13)

Here, Ptot is the total thermodynamic pressure of the system, and P^sat is the saturation pressure of the vaporizing chemical species,evaluatedatthedropletsurfacetemperatureTs.

Theheatandmasstransferratesarespecifiedbyintegratingthe heatandspeciesfluxesoverthedropletsurfaceS(t) inthecarrier phase:

˙ q_d=

S

− 1

RePr

∇

^T

·ndS; (14)

˙ m_d=

S

− 1

ReSc

∇

^Y+uY

·ndS. (15)

3. Numericalmethod

Thegoverning equationsofthe carrierphase 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 carrierphase PDEs, andfor the dis- persedphaseODEsEqs.5to(10).WithintheRunge-Kuttascheme, a pressurecorrection scheme isused forthe Navier-Stokesequa- tion.

ThecalculationoftherighthandsideintegralsoftheNewton- Euler equations Eqs. 9 and (10), the enforcement of the carrier phaseboundaryconditionsatthedropletinterfaceEqs.11to(13), andthecalculationoftheheatandmasstransferratesEqs.14and (15)areall performedontheLagrangian meshpointsthat repre- sentthedropletsurface,whichistreatedasan immersedbound- ary[16,17].The Lagrangian pointsmove rigidlytogether withthe dropletcentroid,anddonotconformtotheEulerianmeshcells,al- thoughtheir spatialresolutionmatchesapproximatelythat ofthe Eulerianmesh. Interpolation fromthe Eulerian tothe Lagrangian meshandspreadingfromtheLagrangiantotheEulerianmeshare performedbymeansoftheregularizedDiracdeltafunction

δ

din- troducedbyRomaetal.[18].

Thepresenceofthephasechangeimpliesthat thedropletsur- faceinjectsmassandenergyintothecarrierphase.Thereforethe followingsourceterms havebeenaddedinorderto mimicmass, energy,andvapourspeciesinjections,consistentlywiththeinter- faceboundaryconditions:

s_{i jk,U}=−dr_d dt

φ

ρ r_d

1+cos

π

^rr^{i jk}d

1−_π⁶2

; (16)

s_{i jk,T}=s_{i jk,U}

1−

φ

cp,v^ap

T; (17)

si jk,Y=si jk,U; (18)

wherer_{i jk}=

|

^xi jk−xc

|

^{foreachEuleriancell.Theparameter}

φ

c_p,vap

(see Table 2) is the specific heat capacity ratio of the vapour speciestothegasmixture.Thespecificformofthesourceterms, 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 zeroforr→r_d,andthattheirintegralisconsistentwiththephase changemassandenergyfluxes.

Adetaileddescriptionoftheimmersedboundaryimplementa- tionandofthesolutionalgorithmforthegoverningequationscan befoundin[12].

4. Flowconfiguration

Thenumericalsetupreproducesaturbulentchannelflowofn- heptanesprayinair.Thetemperatureandpressureconditionsap- proach those commonly encountered in internal combustion en- gines,andthephysicalandtransportpropertieshavebeencalcu- latedfromsaidtemperatureandpressureusingtheaveragingrule knownas“1/3 rule” [19],commonlyemployed inconstant prop- erties evaporation models. The channel walls are adiabatic, with

Table 3

Computational domain and operating conditions.

Descritpion Value Units

d 0 Droplet initial diameter 80 μm

N d Droplet number 14081 -

L x/ d 0 Streamwise channel length 96 - L y/ d 0 Wall normal channel length 32 - L z/ d 0 Spanwise channel length 48 - N x Streamwise resolution 2304 cells N y Wall normal resolution 768 cells N z Spanwise resolution 1152 cells N l Lagrangian points per droplet 1721 points T ∞ Initial gas temperature 741 K T s,0 Initial droplet temperature 342.85 K P tot Thermodynamic gas pressure 43.57 atm

homogeneousNeumannboundaryconditionsforTandY,andthe domainisperiodicinthestreamwiseandspanwisedirections.The streamwiseflowismaintained atconstantbulk Reynoldsnumber byabody force inthemomentumequation that linearlycorrects thestreamwisebulk velocity atevery time step.Inorder tovent out of the periodic domain the additional mass injected by the evaporating droplets (through Eq. 16), a non-zero normal veloc- ityisprescribed atthe walls. Its value is equal tothe integrated volumesource fromall thedroplets, divided by thewall surface, reachingapeak of2%ofthe meanstreamwisevelocitywhen the evaporationisstrongest.Thetangentialvelocitiesaresettozeroat thewalls.

Table 3showsthedetails ofthecomputationalsetup.The val- uesofthenon-dimensional parametersthat characterizethe mo- mentum,species,andenergybalancesareshowninTable2.

The bulk Reynolds number of 5600 corresponds to a friction Reynoldsnumber,basedonhalf thechannel width,ofRe_τ=180, forastatisticallystationarysinglephasechannelflow.

The Eulerian gridresolution andnumber ofLagrangian points on the surfaceof each droplet are chosen in orderto guarantee that thedroplet initial diameter isresolved by 24computational cells.This resolution ensures that boththe turbulencein the gas phase, whose smallest scale is the Kolmogorov scale

η

^{, andthe}

masstransportatthedropletinterface,whosesmallestscaleisthe Batchelorscale

λ

B=

η

/√

Sc,arefullyresolved.Theviscousscaleat thewallisalsoresolved(y⁺≈0.6).

The number of droplets is chosen in order to have an initial liquidvolumeloadingof

φ

0=0.05.

Thepressureandtemperatureconditionsaretypicalofinternal combustionengines.Undertheseconditions,adropletdiameterof 80μmensuresthat theaverageWeber numberofthedroplets is belowunity,sothatdropletdeformationandbreakupisnegligible, andassumptionA4isreasonable.

The flow of the carrier phase is first initialized without the droplets. Following Henningson & Kim [20], a vortex pair is su- perimposedon aplane laminarPoiseuille profile,andallowed to breakdownintoa turbulentspot, asitis carrieddownstreamby themean flow. The flow eventually becomesfully turbulent and developsintoastatisticallystationarystate,withRe_τ=180.

Atthispoint,thedispersedphaseisintroducedintheflow.The dropletsareinitializedinquasi-randompositionswithaLatinHy- percubealgorithm,suchthattheliquidvolume loading

φ

0 isspa- tiallyuniform.It is found that the initial droplet velocity hasno significanteffectonthesimulationresults,owingtotherelatively lowdensityratio

φ

ρ,andcorrespondingdropletStokesnumber¹

Once the droplets are introduced, the evaporation is imme- diately activated by the vapour pressure build-up at the droplet

1The droplet Stokes number is defined as St = ^τ_τ^d

f, where τd= ρld ²/ 18 μ^{is the} droplet relaxation time, and τfis a characteristic time scale of the flow.

surface, and is sustained by the heat provided by the hot gas flow.Thesimulationiscarriedonuntilthecarrierphasegetssat- urated with vapour, whereupon the average droplet evaporation rateissignificantlyreduced, andcondensationeventsstartto be- comeimportant.Thiscorresponds toan averagedropletdiameter ofd/d₀=0.97,whichensures thattheresolutionsoftheEulerian andLagrangianmeshesarestillapproximatelymatchingattheend ofthecalculation,arequirementfortheaccuracyoftheimmersed boundarymethod[21].

5. Results

Fig. 1 shows three instantaneous snapshots of the flow. The backgroundcarrierphaseiscolouredbythevapourmassfraction, 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 thevapourspeciesduringthecourseofthesimulationcanbeob- servedinthevisualizations,andwillbequantifiedinthefollowing.

Basedontheresults,thetransientevolutionofeachofthesephe- nomenacan be roughlydivided intothree stages.We choosethe followingsequence,basedonthedropletevaporationdynamics,as areference:

• StageI:0≤t/t_end࣠0.11;

• StageII:0.11࣠t/t_end࣠0.33;

• StageIII:0.33࣠t/t_end≤1;

wheret_endisthetimeattheendofthesimulation.

Giventhestrongcouplingbetweenthephysicalphenomenain- volved inourcase, it willbe shownthat thisdivision isrelevant forthedescriptionofthedropletmigrationandturbulencemodu- lationaswell.

5.1. Large-scaledropletmotion

The dispersed phase is initially distributed uniformly in the channel,andissettomotionbythebackgroundcarrierflow.Fig.2 showsthelocaldroplet numberdistribution,averagedinthe two statistically homogeneous directions (streamwise x and spanwise z), as a function ofwall distance y andtime t. A large-scale mi- grationofthedropletstowardsthechannelcentrelineisobserved, andbytheendofthesimulationaround50%ofthedropletshave gatheredina0.23lywide bandacrossthechannelcentreline.The average droplet Stokes number, based on the local Kolmogorov time scaleofthecarrierphase flow,isfound tobeSt≈30,estab- lishingthecasewithintheregimedescribedasfour-waycoupling byElghobashi[22].Four-waycouplingofthedispersedphasemo- tion with the carrier phase implies that the carrier phase mean flow, its turbulent fluctuations, the droplet-droplet and droplet- wall collisions, and the localmodification of flow streamlines by the droplet excluded volume, all contribute to determine each droplet’sindividual trajectory,aswellasthelarge-scalemigration ofthedispersedphase.

Thelocalliquidvolumefraction

φ

^{isshowninFig.3,asafunc-}

tion ofthe distance fromthe wall. The quantity is ensemble av- eraged inthe two statisticallyhomogeneous directions (xandz), andtimeaveraged overStagesI,IIandII.DuringStage I,theini- tiallyuniformdropletdistributionisperturbedbythebackground flow;StageIIrepresentsaperiodofongoingbulkmigrationofthe droplets; duringStage III a quasi-steadydistribution isfinally at- tained.Thefinaldistributionischaracterizedby dropletclustering atthecentreline,wherethelocalvolumefractionreachesthemax- imumvalue

φ

=0.11(attheendofStageIII),i.e.morethandouble its initial uniformvalue. The clusteringtapers off gradually away fromthe centreline. A layer ofdroplets is also observedclose to thewall. From thethreeprofiles ofFig.3,it isobservedthat the

Fig. 1. Instantaneous snapshots of the flow for three successive times. Three orthogonal planes xy, xz, yz are shown, with colour contours of the vapour mass fraction Y . For clarity, only the droplets lying within the xy and yz planes are shown, and they are coloured with their surface temperature T s, normalized with the initial gas temperature T ∞.

formationofthewalllayerisafasterprocessthanthebulkmigra- tiontowardsthecentreline.

Thedispersed phasebehaviourobservedinthepresentcaseis similar to theone reported by Fornarietal. [23],fora turbulent channel flowladenwithsolid particles,withconditionssimilarto the presentcase(Re_τ=180,

φ

⁰=0.05,

φ

ρ=10,Ly/d0=18).The final droplet distribution,with thecoexistence ofdroplet cluster- ing atthecentreline andadroplet walllayer, isexplainedasthe

resultoftwoconcurrent mechanisms: shear-inducedmigrationto- wardsthechannelcentreline,andturbophoresistowardsthechan- nelwalls.

Shear-induced migration is the net displacement of the dis- persed phase from regions of high shear-ratetowards regions of low shear-rate, i.e. away from the wallsin a wall boundedflow.

Thisis theresultof irreversibleshort-range interactions between the droplets in the presence of shear. While classical literature

Fig. 2. Droplet number distribution along the wall normal coordinate, as a function of time.

Fig. 3. Local liquid volume fraction φ^{as a function of the wall distance.} Comparison between Stages I, II and III.

onshear-induced migration focuseson Stokes flow [24–26], it is reasonable to expect that the phenomenon can occur at higher Reynoldsnumbers andalsointurbulentconditions,providedthat theStokes numberishighenough forthe dropletsto escape the carrierphasemeanflowstreamlinesafteracollisionevent.Infact,

it isreportedby Fornariet al. [23] that increasingcentreline mi- gration is observed for increasing values of the dispersed phase density(i.e.forincreasingStokesnumber),whenthecarrierphase turbulenceanddispersedphasevolumefractionarekeptfixed.The work ofFornarietal. [23] dealswithstatisticallystationary flow, whereas in thefollowing we take advantage of thetransient na- ture of our flow configuration, and present dynamical evidence in favour of the hypothesis that the Stokesian theory of shear- induced migration isapplicable in turbulent conditions,forsuffi- cientlylargedispersedphaseStokesnumber.

Theotherobservedmechanism, turbophoresis,isthemigration ofthedispersedphasetowardsregionsoflowturbulenceintensity [27].Inawallboundedflowtheseregionscorrespondtothevicin- ityofthewall.ItisarguedbySardinaetal.[28]thatturbophoresis isalarge-scalemanifestation, inthepresence ofturbulenceinho- mogeneity,ofsmall-scaleclusteringofthedispersedphaseintore- gionsofhighstrainandlowvorticity,amoregeneralphenomenon thatisalsoobservedinhomogeneousturbulence[29–31].

The time scale ofshear-induced migrationis governed by the rate of non-hydrodynamic droplet-droplet interactions, which is indirectly affected by the hydrodynamics of the carrier phase throughthe meanshearrategradient[26].Turbophoresis, onthe other hand,is a purely hydrodynamiceffect,whosetime scale is governedbythecorrelationofthegradientofthevelocityfluctua- tion[27],i.e.itisrelatedtotheTaylormicro-scale[32].Separation of time scales between the two phenomena, as observed in the presentcase,isthereforelinkedtothevaluesoffourflowparam- eters,namelythecarrierphaseReynoldsnumber,determiningthe ratio of droplet diameter to turbulent length scale, the liquid to gasdensityratio

φ

ρ,determiningthe dropletStokes number,the liquid volume fraction

φ

0,determining the droplet numberden- sityandthusthefrequencyofdroplet-dropletinteractions,andthe ratioofchannel widthto dropletdiameterL_y/d₀,determiningthe frequencyofdropletinteractionswiththewalls.

Thedynamicsofthelarge-scalemigrationare showninFig.4, in terms of the droplet wall normal mean velocity, averaged in thetwohomogeneous directions,andtherootmeansquare ofits (a)

(b)

Fig. 4. Distribution along the wall normal coordinate of the droplet average wall normal velocity component (a) and its root mean square fluctuation (b), as functions of time.

Fig. 5. Instantaneous snapshots of the local liquid volume fraction φ, as a function of the wall normal coordinate, averaged in the two homogeneous directions, and filtered with a spatial window equal to (5 L y/ N y) along the wall normal direction. The shear-induced macroscopic collision event that triggers the centreline migration, in the form of a collision wave, is indicated.

fluctuation, asfunctionsof wall distanceandtime. A large num- ber ofdroplets undergoes a velocityreversal, marked by theyel- low dashed line in Fig. 4(a): they drift towards the wall in the beginning,asaresultofturbophoresis,buttheyturntowardsthe centreline aftera majorcollision eventtakesplace inthevicinity ofthe dropletwall layerandtriggers shear-inducedmigration.In Fig. 5,displaying eightsuccessive instantaneous snapshots ofthe localliquid volumefractionprofile alongthe wallnormal coordi- nate,we marktheabovementioned collisionevent, andthesub- sequent“collision wave” that travelstowards thechannel centre- line.TheprofilesofFig.5havebeenfilteredwithaspatialwindow equalto(5Ly/Ny)forbetterclarity.

Thespeedofthecollisionwaveisestimatedastheslopeofthe yellowdashedlineinFig.4(a),as

v

wav^e/(

ν

^Re/d₀)=2.2×10⁻³.Itis interestingtonote that,accordingtotheclassicaltheory ofshear- induced migration in Stokes flow [26], the shear-induced migra- tion velocitycan be estimatedas

v

migr∼ −d₀^2∂_∂²_y^U₂,whereU isthe streamwisevelocityofthecarrierphase.Forthepresentcase,after

calculatingthecarrierphasemeanvelocityastheensembleaver- ageinthetwohomogeneousdirections,wespace-averageitssec- ondorder gradientin thewall normaldirectionacross thechan- nel 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 startingtime of the collision wave, as seen inFig. 5. This givesan estimation of

v

migr/(

ν

^Re/d₀)=2.3× 10⁻³, falling remarkably close to the calculated value of vwave. This finding substantiates our hypothesis that Stokesian shear- induced migration dynamics are still active in turbulent flow, provided that the dispersed phase Stokes number is sufficiently large.

Thereducedwallnormalvelocityfluctuationsatthecentreline, visiblein Fig.4(b), arean effectofthe increasedconfinement of thedispersedphaseinthe centrelineregion,wherethelocalvol- umefractionreachesits peak value,andis stronglycorrelated to theturbulence attenuation inthe sameregion, aswill be shown inSection5.2.The quenchedwallnormalvelocityfluctuationsef-

(a)

(b)

Fig. 6. Comparison of turbulence statistics of the carrier phase, between the present case (solid lines) and a single-phase channel flow with the same Re τ= 180 (dashed lines, from [34] ). Root mean square velocity fluctuation profiles along the wall normal coordinate on the left, profiles of the first non-diagonal component of the Reynolds stress tensor along the wall normal coordinate on the right. The quantities are normalized with the friction velocity. Time averages over Stage I (a), Stage II (b), first half of Stage III (c), second half of Stage III (d) are shown. (Continues on the next page).

fectivelytrapthedropletsthathavereachedtheregionaroundthe centreline,stabilizingthecentraldropletcluster.

We shallseeinSection 5.3that thelarge-scaledropletmotion inthechannel,inparticular theshear-inducedmigrationtowards thecentreline,hasamacroscopiceffectontheevaporationdynam- ics.

5.2.Modulationofturbulence

Inthepresentflow conditions(Re=5600,

φ

0=0.05,

φ

ρ=32, Ly/d₀=32), the ratio of droplet size to Kolmogorov length scale isapproximatelyd₀/

η

≈20,andthe ratioofdropletsize toTaylor lengthscaleisapproximatelyd₀/

λ

t≈4,whenthedispersedphase isreleasedintothefullydevelopedcarrierphaseturbulence.It is thereforeexpectedforthedispersedphase tohaveastronginflu- enceonthecarrierphaseturbulence,asthedropletsizefallsinto theinertialsubrangeofturbulentlengthscales[31,33].

In the following, turbulencestatistics of the carrierphase are analysed.Thestatistics,conditionedtothecarrierphase,havebeen calculatedneglectingtheregionsoccupiedbytheliquidphase,and ensembleaveraginginthetwohomogeneousdirections.

Profiles of thesquare rootof thediagonal components ofthe turbulentReynoldsstresstensor(rootmeansquarevelocityfluctu- ations)andofthefirstoff-diagonalcomponent,alongthewallnor-

malcoordinate,areshowninFig.6,normalizedwiththeinstanta- neous friction velocity u_τ. They are time averaged over Stages I, II andIII, withStage III dividedinto two sub-stages,to better il- lustratetheevolution oftheReynoldsshearstress.The quantities are comparedwiththe corresponding valuesfroma singlephase channelflowwithRe_τ=180(datafrom[34]).

Itis observedinFig.6(a)that, whenthedroplets arereleased into the flow, and for an initial transient period that lasts ap- proximately untilthe beginning ofthedroplet migration towards thecentreline(StageI),thecarrierphaseturbulenceisgreatlyen- hancedcompared tothe singlephasecase. In particular,theroot mean square velocity fluctuationsexhibit a flat profile along the channel width, asidefromthe near-wall region where they drop to zero. This reflects the roughly uniform droplet distribution in thechannel, duringthisperiod.DuringStage II, i.e.themain pe- riod of droplet migration, the carrier phase turbulence starts to decrease,especiallyintheregion acrossthecentreline,asacom- parison between Fig. 6(b) and Fig. 6(a) shows. During Stage III, thebulkdropletmigrationslowsdownuntilthequasi-steadydis- tribution of the dispersed phase shown in Fig. 3 is attained. At this stage, the flow is clearly divided into two regions with dis- tinctturbulentfeatures,asevidentfromFigs.6(c,d).Whent=t_end, theregionaroundthecentreline(0.15࣠y/Ly࣠0.85),wherearound 93% of the droplets have gathered, has a much lower turbulent

(c)

(d)

Fig. 7. Comparison of turbulence statistics of the carrier phase, between the present case (solid lines) and a single-phase channel flow with the same Re τ= 180 (dashed lines, from [34] ). Root mean square velocity fluctuation profiles along the wall normal coordinate on the left, profiles of the first non-diagonal component of the Reynolds stress tensor along the wall normal coordinate on the right. Quantities are normalized with the friction velocity. Time averages over Stage I (a), Stage II (b), first half of Stage III (c), second half of Stage III (d) are shown.

intensity than the corresponding single phase case: in particu- lar the values of

v

_rms, w_rms and u

v

are very low and suggest flow re-laminarization. A comparison with Fig. 4(b) shows that the diminished wall normal velocity fluctuation atthe centreline is also strongly correlated to the small value of the wall nor- mal droplet velocity fluctuation, hence the quasi-steadydistribu- tion ofthedispersedphase iscoupledtotheturbulenceattenua- tion.The valueofthe streamwiserootmeansquarevelocity fluc- tuationu_rmsinFig.6(d)isstillmuchlargerthanitisinthesingle phase case; howeverwe interpretthis resultas thesampling by thecarrierphasevelocity ofthemanylaminar-likedropletwakes presentinthisregion,ratherthanasactualturbulent fluctuations ofthecarrierphaseflow.Thisinterpretationwasfirstsuggestedby Parthasarathy&Faeth[35],whopointedoutthat thecontribution of mean streamwise velocity cannot be entirely separated from the contribution ofturbulence in theparticle wakes, since parti- cle arrivalsare random. Outside ofthe central region (y/Ly࣠0.15 andy/Ly࣡0.85),onlyaround7%ofthedropletsremainatt=t_end, which makes the local liquid volume fraction

φ

≈0.001. Thus it is notsurprisingthat theturbulencecharacteristics inthisregion are similar tothose ofthesinglephase channel, butsqueezedin thewall normalcoordinate,astheregionoccupiesonly ∼30%of the channelvolume. Inparticular, asevident fromFig.6(d),

v

rms,

w_rms reachaty/Ly≈0.15the value thatthey attain atthecentre- line in the single phase case, while u

v

has an inflection point that marks the boundary between the two flow regions at ap- proximately the same location. The slightly augmented velocity fluctuationsandReynolds stress inthe immediatevicinity of the walls(y/Ly࣠0.05 andy/Ly࣡0.95), noticeable in Fig. 6(d), can be explainedbytheexcluded volumeeffectofthedropletcentreline cluster,whichsqueezes partofthecarrierphase flowinthenear wallregion,locallyincreasing theReynoldsnumber.Theseresults arein linewiththe findings ofFornarietal. [23] for thecaseof turbulentchannelflowladenwithsolidparticlesatsimilarcondi- tions(Re=5600,

φ

0=0.05,

φ

ρ=10,Ly/d₀=18).

Asymptotically (t>0.67t_end), it is found that the overall tur- bulence level is damped by the droplets. In fact, fora dispersed phasewhichisdenserthanthecarrierphase,andwithsizelarger than the Kolmogorov scale and smaller than the integral scale of turbulence, the higher specific inertia (

φ

ρ>1) and the en- hanced dissipation in the inertial subrange arising from the dis- persed phasedrag act assinks ofthe carrierphase turbulent ki- netic energy[36]. The diminished turbulent transport hasa con- siderable effect on the evaporation dynamics: the wall normal profiles of the turbulent vapour species wall normal flux, dur- ing Stages I,II, and III, are displayedin Fig. 8, showing that the

Fig. 8. Profiles of the turbulent wall normal flux of the vapour species, along the wall normal coordinate. Time averages over Stages I, II and III are shown.

fluxisreducedbyone orderofmagnitudeoverthecourse ofthe simulation.

5.3.Evaporationdynamics

Inthepresentconfiguration,theevaporationisactivesincethe beginningofthedroplets’ lifein theflow. Thereforetheevapora- tiondynamics evolve simultaneously with the other phenomena, namely the droplet migration and development of quasi-steady droplet distribution across the channel, the increase and succes- siveattenuationofthecarrierphaseturbulenceintensity,andthe progressivesaturationofthecarrierphasewithvapour.Theevolu- tionoftheindividualdropletdiameterandsurfacetemperatureis showninFig.9(a,b),together withtheaverageoverall the14081 droplets(black line).The threedistinct phasesoftheevaporation dynamics,whichwereidentifiedasStagesI,IIandII,canbechar- acterizedasfollows:

1.Stage I, 0≤t/t_end࣠0.11. The droplet temperature rises, first sharplyandthenmoregently,asthedropletsareheatedbythe surrounding gas, and the droplet average diameter decreases asaresultoffastevaporation. Werecallthat thedropletcen- treline migration is triggered at t/t_end≈0.11, asshown in the fourthpanelofFig.5.

2. StageII,0.11࣠t/t_end࣠0.33.Thedroplettemperatureslowlyde- creases,astheimbalancebetweenthesensibleheattransferred to the droplets by the gas (q˙_d/

φ

^cp<0 in Eq. 6) and the la- tent heat of evaporation (m˙_d/Ste>0 in Eq. 6) shifts towards thelatter.As a consequence,theevaporationslows downand thedropletdiameterdecreasesatalowerrate.

3. Stage III, 0.33࣠t/t_end≤1. The mean droplet temperature and meandroplet diameterevolve towardsan asymptotic value:a dynamic equilibrium is eventually established between evap- oration and condensation events. During this period, the fi- nal droplet distribution of Fig. 3 is established, the carrier phase becomes saturated with vapour, and its turbulent fea- turesevolvetothefinalconfigurationofFig.6(d).

Thefinaldistributionsofthedropletdiameteranddropletsur- facetemperatureareshowninFig.9(c,d).Thedropletsurfacetem- peraturedistribution (Fig.9(d)) isbimodal,havingtwo veryclose but distinct peaks. The negative skewness of both distributions reflects the asymptotic droplet distribution across the channel (Fig. 3), suggestingthat the dropletcluster around the centreline includesagroupofdropletsthathaveexperiencedstrongerevap- oration than the majority, and end up being smaller andcolder.

This is confirmed by the distribution of these tail droplets (i.e.

those withd/d₀<0.96) along the wall normal coordinate,shown inFig.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 ratealong thewall normalcoordinate,time averaged overStages I,IIandIII.Themeanvapourmassfractioninthecarrierphase,as afunctionofthewallnormalcoordinate,isshownontheside.The evaporationrateK=^d_d^d_t² isheredefinedas thedimensionlessrate ofchangeofthedropletsurfacearea;itisnegativeforevaporation (K<0)andpositiveforcondensation(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 fractioninthe carrierphase isalso homogeneous andlow,witha smalldip closetothe wall, corre- spondingto thenarrowregion depletedofdroplets visibleinthe firstprofileofFig.3(y/L_y࣠0.05andy/L_y࣡0.95).

(a) (b)

(c) (d)

Fig. 9. Top panels: (a) Evolution of the droplet diameter and (b) droplet surface temperature in time. The grey lines represent individual droplets (for clarity, a random sample consisting of 1/50 of the total droplets is shown); the thicker black line represents the mean. Bottom panels: final distribution, at t = t end, of the (c) droplet diameter and (d) droplet surface temperature.

Fig. 10. Joint probability density function of the droplet diameter and droplet wall distance, at t = t endfor the droplets with d / d 0< 0.96.

(a)

(b)

(c)

Fig. 11. Joint probability density function of the droplet evaporation rate and droplet wall distance, normalized with the number of droplets in each wall normal bin. Time averages over (a) Stage I, (b) Stage II and (c) Stage III are shown. The corresponding time averages of the mean vapour mass fraction profile in the carrier phase, along the wall normal coordinate, are shown on the left panels.

Astheevaporationadvances(Fig.11(b)),thevapourmassfrac- tionbuildsupinthecarrierphase;howeveritdoessoataslower 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,theevaporationrateofthedropletsthat lingerout- side of thecentral region stays higher thecloser they are tothe

wall.Duringthisstage,thedropletwalllayerisfoundtohavelit- tle tono influenceon theevaporationratedistribution. The cen- tralregion(0.1࣠y/L_y࣠0.9)ischaracterizedbyauniformdistribu- tionoftheevaporationrateacrossits width,andahomogeneous vapour mass fraction in the carrier phase, despite the fact that the droplet distribution itself is not uniform (Fig. 3). Condensa-

(a) (b)

(c) (d)

(e) (f)

Fig. 12. Time development of some mean flow quantities: (a) carrier phase Reynolds stress; (b) carrier phase turbulent vapour species flux; (c) carrier phase and droplet surface vapour mass fraction, (d) droplet evaporation rate; (e) droplet vaporization Damköhler number; (f) droplet wall normal velocity.

tioneventsalsostarttoappearwheretheliquidvolumefractionis large.

Similar considerations apply to Stage III (Fig. 11(c)), when a quasi-steadystate is reached, with quasi-steadydroplet distribu- tion across the channel and dynamic equilibrium of evaporation andcondensation events.Thecentreline cluster ismore compact, thusitsregionofinfluenceisnowsmaller(0.3࣠y/L_y࣠0.7).Asig- nificantamountofcondensationeventstakesplaceinthisregion:

thisis reflected inthe noticeable local drop of the vapour mass fractioninthecarrierphase.Theregionoutsideoftheinfluenceof thecentreline cluster isnow bigger(y/Ly࣠0.3andy/Ly࣡0.7): as before,thetailoftheevaporationratedistributionbendstowards strongerevaporationasitapproachesthewall,buttheinfluenceof thedropletwalllayerhasbecomeappreciable,andtheevaporation ratedistributionatthewallhasnowalargevarianceowingtothe appearanceofcondensationevents.

Byshowingsidebysidethetimedevelopmentofvariousmean quantities of the carrier anddispersed phases, Fig. 12 illustrates the influence of large-scale droplet motion, carrier phase turbu- lenceandcarrierphasemasstransportontheevaporationdynam- ics.Timest/t_end=0.11andt/t_end=0.33,whichdelimitStagesI,II andIII,aremarkedinthefigure.

Stage I, characterized by droplet heating and fast evapora- tion, ends at t/t_end≈0.11, coinciding with the beginning of the large-scalecentrelinemigration,asevidencedbythemeandroplet 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 isincreasing (growth ofthe mean ReynoldsstressofFig.12(a)),augmenting theheattransfertothe droplets and hence the droplet heating (Fig. 9(b)). Consequently, eventhough the vapourmass fractionYstartstobuild up inthe bulk carrier phase, the droplet temperature rise causes a faster increase of the vapour mass fraction at the droplet surface Ys

(Fig.12(c)).Sincethedifference(^Ys−Y)^{isthedrivingforceofthe} evaporation, a growthof the evaporationrateis observedduring thisinitial sub-stage,anda peak value isreachedatt/t_end≈0.02.

Turbulenceattenuationbeginsatt/t_end≈0.02(Fig.12(a)):theheat transferfromthegastotheliquidphaseisreducedandthedroplet heatingstops(Fig.9(b)),which haltstheincreaseofvapour mass fractionatthe dropletsurfaceYs (Fig.12(c)),sothat theevapora- tionratestartstodecrease(Fig.12(d)).Interestingly,thereisalag betweenthe average turbulent vapour species flux in thecarrier phase andthe evolutionof turbulence:Fig. 12(b)shows that,for 0.02࣠t/t_end࣠0.11,

v

Y increasesandthen saturates, whileu

v

is decreasing.Thismay bedueto thefact that velocityfluctuations decayfasterthanmassfractionfluctuations(Sc>1).

During Stage II, 0.11࣠t/t_end࣠0.33, the droplets accelerate to- wards the centreline (Fig. 12(f)). Simultaneously, the turbulent fluxesinthecarrierphasedecrease(Fig.12(a,b)),thedropletscool down(Fig.9(b)),andthevapourmassfractionatthedropletsur- facedecreaseswhileitkeepsaccumulatinginthebulk(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 ve- locity(Fig.12(f)),andthedroplets starttodecelerateandattaina

quasi-steadydistribution(Fig.3).ThismarksthebeginningofStage III,whenthecarrierphaseturbulentmomentumfluxdecaystoits asymptoticvalue (Fig.12(a)),the turbulentvapour species fluxis progressivelyextinguished (Fig.12(b)), thecarrierphase becomes saturatedwithvapour(Fig.12(c)),andtheevaporationratefurther decreasestoalmostzero(Fig.12(d)).

5.3.1. Turbulentscalingoftheevaporationrate

Fig.12(e)showsthetimeevolutionoftheaveragevaporization Damköhler number ofthe dispersed phase. The conceptof a va- porizationDamköhlernumberwasfirstintroducedbyGökalpetal.

[37],inanalogywithturbulentcombustion,asawaytocharacter- izetheinfluenceofturbulenceonthevaporizationdynamics.Itis defined as the ratio of the characteristic time scale of turbulent eddiestothecharacteristictimescaleofvaporization:

Da_v=

τ

eddy

τ

vap. (19)

Therelevantturbulenttimescaleisthatoftheeddieshavinga length scaleoftheorderofthe instantaneousdropletdiameterd. When thedropletsizefallswithintheinertialsubrange ofturbu- lence,

τ

eddycanbeestimatedas:

τ

eddy= d^2/3

ε

^{1/3; (20)}

where

ε

^{is the local} dissipation rate of turbulent kinetic energy.

The vaporizationtime scaleisestimatedastheadvectiontimeof thevapouracrossthemasstransferboundarylayeronthedroplet surface:

τ

v^ap=

δ

M

|

φ

ρ−1

drd

dt

|

^{. (21)}

Theboundarylayerthickness

δ

Misestimatedfromthefilmthe- oryofAbramzon&Sirignano[38]:

δ

M= d

Sh−2; (22)

Sh=2+ BMSh0

(

¹+BM

)

^0.7ln

(

¹+BM

)

^{; (23)} Sh₀=2+0.552Re¹_d^/2Sc^1/3; (24)

B_M=Ys−Y_∞

1−Ys ; (25)

where Shis the vaporization Sherwoodnumber, Sh₀ is theSher- woodnumberformasstransferaroundasphere,accordingtothe correlation by Ranz & Marshall [39], Re_d=

|

^ud

|

^d/

ν

^{is thedroplet}

Reynolds number, B_M isthe Spalding mass transfer number. The free stream vapour mass fractionY_∞ is estimatedas theaverage of Y on a spherical shell centred on the droplet centroid, ofra- dius

r_d+

^x2+^y2+^z2

,where(^{x,y,z)isthesizeof} aEulerianmeshcell.

Fig.12(e)showsthatthetimeevolutionofDa_visverysimilar to that of theturbulent vapour species flux

v

Y, suggestingthat themainmechanismbywhichthecarrierphaseturbulenceinten- sifiestheevaporationistheejectionofvapourfromtheboundary layerbycorrelatedvelocityandmassfractionfluctuations.

Given the transient nature of our simulation, the dispersed phase experiences a wide range of evaporating conditions, and the average vaporization Damköhler number varies in the inter- val [5.4×10⁻³,1.6×10⁻¹],spanning abouttwo decades. Thus, a sizeablesampleisavailabletosuggestameaningfulcorrelationbe- tween turbulenceandevaporation. Fig.13relatesthemeanevap-

Fig. 13. Evaporation rate enhancement by turbulence K/K l, as a function of the va- porization Damköhler number Da v. Data points for t / t end≥0.11 are shown.

orationrateenhancementbyturbulence,K/K_l,definedbynormal- izingtheevaporationratewiththeevaporationrateunderlaminar conditions,tothemeanvaporizationDamköhlernumberDa_v.The laminarevaporationrateisestimatedas:

K_l=4 Sh ReSc

ln

(

^1+BM

)

φ

ρ . (26)

Data pointsforStage I (t/t_end<0.11)have beendeemed unre- liable with respect to the calculation of K_l, owing to inaccurate evaluationoftheSpaldingmass transfernumberB_M,andarenot reportedinFig.13.Fort/t_end≥0.11,itisfoundbyaleastsquaresfit thattheturbulentenhancement oftheevaporationratedecreases withtheDamköhlernumberas:

K/Kl=0.970Da⁻⁰_v ^.1. (27)

This trend is remarkably similar to the correlation K/K_l= 0.771Da⁻⁰_v ^.111, found by Wu et al. [40], and obtained after a se- riesofevaporationexperiments ofasingleliquidfuel dropletfor variousfree-stream turbulenceintensities andscales [40,41].It is thusshown forthefirst time that aDamköhler powerlaw holds evenwhenthedroplets haveastrongfeedbackontheturbulence field, andin the presenceof a denser sprayfor which theinflu- enceofneighbouringdropletsontheevaporationdynamicscannot beneglected.Thus,Eq.27revealsamorefundamentalmechanism ofmasstransfer than whatwaspreviously reportedby Wu etal.

[40,41], whoseexperiment waslimitedto isolated droplets inan externallycontrolledturbulentstream.

Itisthehopeoftheauthorsthatthepresentresultattractsnew attentiontothedistinctionbetweenturbulenceeffects(Damköhler number)andmeanfloweffects(Sherwoodnumber) ontheevap- oration dynamics, a fundamental point that is still neglected by mostevaporationmodels.

6. Conclusionsandoutlook

Thisworkpresentstheresultsofadirectnumericalsimulation of sprayevaporation in turbulent channel flow. The 14081 spray dropletsare interfaceresolved, andtheir motionisfour-waycou- pled to the carrierphase. Thus, this effortrepresents one of the largest numerical studies of phase change flow to date, both in scopeandinthedetailofthephysicalphenomenadescribed.

The conditions explored in the simulation lead to the emer- genceof severalinterrelated phenomena. The dropletsize lies in theinertialsubrangeofthecarrierphaseturbulence,whichleads tosignificantturbulencemodulation.ThedropletStokesnumberis intermediate (St≈30), which allows the coexistence and interac- tionofhydrodynamicandnon-hydrodynamiceffects, suchastur- bophoresis and shear-inducedmigration. The choice of the stan- dard turbulent channel flow represents a convenient andwidely