Oguz Emre
email: oguz.emre@ecp.fr

Rodney O Fox
email: rofox@iastate.edu

Marc Massot
email: marc.massot@ecp.fr

Stéphane De Chaisemartin
email: stephane.de-chaisemartin@ifpen.fr

Stéphane Jay
email: stephane.jay@ifpen.fr

Frédérique Laurent
email: frederique.laurent@ecp.fr

Eulerian modeling of a polydisperse evaporating spray under realistic internal-combustion-engine conditions

Keywords: Eulerian models, polydispersity, high-order moment methods

come

high-order (in size) moment model has proved to be well adapted for injection simulations with moving geometries. Moreover, it requires less computational effort as compared to existing methods, with a single section for the size phase space. Academic test cases have demonstrated its great potential for industrial applications using one-way coupling. In order to draw comparisons with experimental data, a two-way coupling framework accounting for the droplet-gas turbulence interactions is developed and validated through homogeneous test-cases under both an academic framework and realistic internal combustion engine conditions.

Introduction

In internal combustion engines (ICE), the direct injection of the liquid fuel jet inside the combustion chamber has a great influence on both fuel consumption and pollutant production. The dynamics of the jet being very fast, one observes rapid temporal and spatial variations of mass, momentum and energy of the flow inside the chamber. Moreover, the phase change prompted by the evaporation of the fuel along with the turbulent character of the flow further complicates the physics. The liquid jet is composed of a dense zone near the injector and a dispersed zone with a cloud of droplets, also called a spray, downstream of the injector. Many previous works have addressed the modeling of the two-phase flow composed of the gas and liquid droplets. For such simulations, Lagrangian methods have been widely adopted since they provide high numerical efficiency [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF]. Furthermore, their implementation in computing software is also quite straightforward. However, the high number of particles required for statistical convergence increases the simulation cost. Lagrangian methods also introduce numerical difficulties related to the coupling with the Eulerian grid used for the gas phase around the droplets and the dense zone of the spray. Moreover, they face difficulties on parallel architectures due to possible heterogeneous load balancing between the processors [START_REF] Garcia | Développement et validation du formalisme Euler-Lagrange dans un solveur parallle non-structur pour la simulation aux grandes chelles[END_REF]. Given these shortcomings of Lagrangian methods, the Eulerian description is considered as a promising alternative. However, the precise description of the polydispersity and of the two-way turbulent interaction of evaporating droplets with the surrounding gas phase at a reasonable cost remain a challenge for Eulerian methods. Moreover, satisfying the moment method stability constraints is complicated due to the rapid variations occurring inside the flow domain. Nevertheless, some pioneering work has been carried out in the context of one-way coupling [START_REF] Kah | Taking into account polydispersity for the modeling of liquid fuel injection in internal combustion engines[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF].

The Eulerian Multi-Size Moment (EMSM) method, derived from the Williams-Boltzmann kinetic model [START_REF] Williams | Spray combustion and atomization[END_REF] and developed in [START_REF] Kah | Taking into account polydispersity for the modeling of liquid fuel injection in internal combustion engines[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], has shown promising potential for fuel-injection applications based on academic configurations [START_REF] Kah | Modeling of polydisperse sprays using a high order size moment method for the numerical simulation of advection and evaporation[END_REF]. Its dedicated numerical scheme, respecting the moment stability constraints, treats precisely the evaporating droplets and requires less computational effort as compared to other methods [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF]. Moreover, two major advances towards industrial applications have been recently achieved [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF]. The first is the adaptation of EMSM method to an unstructured, staggered moving grid under the Arbitrary Lagrangien Eulerian (ALE) numerical formalism. The second is the development of a stable and accurate numerical strategy for treating the polydisperse two-way coupling of the evaporating spray with its surrounding gas, while respecting the conservation of moment space, i.e., the stability requirement of the method. All these developments have been integrated into the industrial computational fluid dynamics (CFD) software IFP-C3D dedicated to compressible reactive flows under ICE conditions. Verification tests between the Eulerian and Lagrangian descriptions of the spray under injection conditions with a d 2 -constant evaporation law and two-way coupled polydisperse droplets have been successfully carried out [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF]. However, the issue of turbulent droplet dispersion and its impact on two-way coupling dynamics still needs to be addressed to move towards more realistic engine conditions. These topics are the focus of this work.

The reminder of this paper is organized by the following way. The second section discusses the derivation of an Eulerian-Eulerian model, within the framework of laminar two-phase flows composed of an evaporating polydisperse spray and a compressible gas. The correct behavior of the energy partition in the spray phase for the turbulence modeling requires taking into account the granular temperature effect (also called uncorrelated motion), as highlighted first in [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF]. This accounts for considering a poly-kinetic velocity distribution at the kinetic level. Based on this idea, first the original monokinetic EMSM model is extended to poly-kinetic in the context of laminar flow, through a transport equation for the granular temperature. The gas phase is modeled with the compressible Navier-Stokes equation. The third section is devoted to a new Reynolds-averaged (RA) turbulence model derived from the two-phase model presented in the second section. This is based on the same philosophy introduced in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF] for two-way coupled monodisperse flows. Here, one must deal with new terms and equations that arise due to size moment equations of the polydisperse evaporating spray and the gas-phase internal energy equation. To overcome this bottleneck, new closure models are provided and discussed. The fourth section is dedicated to homogeneous test cases. First, the new model is qualitatively validated as compared to the test case of [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF] for one-way coupling and then the extension to two-way coupling is studied for both evaporating and non-evaporating sprays. Next, the model is investigated under the conditions typical of high-pressure direct injection in ICE applications. Finally, some relevant conclusions along with several insights on future work are discussed in section five.

Modeling approach

A cloud of droplets undergoing Brownian motion (i.e., a spray) can be described using the statistical formalism originally proposed by [START_REF] Williams | Spray combustion and atomization[END_REF] for combustion and atomization applications. This formalism is appropriate for the disperse-flow regime where inertial forces, leading to a nonzero relative velocity between a droplet and the gas, are not very large compared to droplet surface-tension forces, characterized by a small Weber number We < 12, and is appropriate for dilute sprays, which have a large droplet mean free path with respect to the characteristic length of the flow (i.e., Knudsen number Kn > 0.1). In this context, droplets can be assumed to be roughly spherical, and their number concentration can be described using a number density function (NDF). The NDF f is defined such that f (t, x, u, S, T )dxdudSdT represents the number of droplets residing in the small volume [x, x + dx] having velocities between [u, u + du], sizes (i.e., surface area) between [S, S + dS], and temperatures between [T, T + dT ] at time t. The evolution equation for the NDF, first introduced in [START_REF]Combustion Theory (Combustion Science and Engineering Series[END_REF], is

∂ t f + ∇x • (uf ) -∂ S (R S f ) + ∇u • (F f ) = 0 (1)
where F is the drag force per unit mass, and R S ≥ 0 is the drift velocity due to evaporation. In this work, the spray is assumed to be collisionless at the far downstream of the injector since the spray volume fraction is very small Φv < 10 -3 . Moreover, the assumption of a d 2 -evaporation law makes the heat-exchange term negligible. Therefore Eq. ( 1) contains no additional terms (e.g., for coalescence or heat transfer). For the sake of simplicity, the drag term is assumed to obey Stokes law:

F = 1 τ d (ug -u) with τ d = ρ d S 18πµg (2)
where τ d is the dynamic time scale associated with droplets of size S, ug is the gas velocity seen by droplets, ρ d is the material density of the fuel droplets, which is assumed constant, and µg is the gas dynamic viscosity.

As in [START_REF] Laurent | Multi-fluid modeling of laminar polydispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], an assumption on the form of f for closure at the kinetic level must be made. Here we assume, as done in [START_REF] Massot | Eulerian multi-fluid models for polydisperse evaporating sprays[END_REF][START_REF] Reveillon | Analysis and modeling of the dispersion of vaporizing polydispersed sprays in trubulent flows[END_REF] , that the velocity dispersion around the mean velocity u d is independent of size such that f takes the following form:

f (t, x, S, u) = n (t, x, S) φ (u -u d (t, x)) ( 3
)
where φ is an isotropic Gaussian distribution:

φ(u) = 1 (4πΘ/3) 3/2 exp - 3u 2 4Θ . (4)
The macroscopic variable Θ(t, x), which describes isotropic velocity dispersion due to the nonzero Stokes number of the droplets, is expressed as a function of the droplet velocity as

Θ = 1 2 u 2 φ (u) du. (5)
We refer to Θ as the granular temperature of the droplets. In the framework of the moment approach, semi-kinetic equations are derived, integrating the moments of the NDF over droplet velocity phase space. Taking velocity moments of order 0, 1 and 2 of Eq. ( 1), the following system is obtained:

∂ t n + ∇x • nu d = R S ∂ S n, (6a)
∂ t nu d + ∇x • n (u d ⊗ u d + P d I) = n ug -u d τ d + R S ∂ S nu d , (6b)
∂ t n u 2 d 2 + Θ + ∇x • nu d u 2 d 2 + Θ + P d = n ug • u d -2Θ -u 2 d τ d + R S ∂ S n u 2 d 2 + Θ . ( 6c
)
The pressure tensor P d I is a second-order isotropic tensor, which is related to Θ by

P d = 2 3 Θ ( 7
)
where the value of Θ is found from Eq. (6c).

System [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF] still contains terms that are functions of the size phase space variable S. To retain information on the size distribution, we follow the classic approach of moment methods and establish conservation equations for the macroscopic quantities defined as

m k = S k n(t, x, S)dS for 0 ≤ k ≤ 3 ( 8
)
where k is an integer and m k denotes moments of order k of the size distribution. Through successive integrations in S k of Eq. (6a) and in S of Eq. (6b), we can obtain conservation equations for a polydisperse spray [START_REF] Kah | Taking into account polydispersity for the modeling of liquid fuel injection in internal combustion engines[END_REF]. Note that the term Θ can be closed by integrating in S over Eq. (6c). The final spray governing equations are given by

∂ t m 0 + ∇x • m 0 u d = R S n(t, x, S = 0), ( 9a
)
∂ t m 1 + ∇x • m 1 u d = M, (9b)
∂ t m 2 + ∇x • m 2 u d = -2R S m 1 , (9c)
∂ t m 3 + ∇x • m 3 u d = -3R S m 2 , ( 9d
) ∂ t m 1 u d + ∇x • m 1 u 2 d + m 1 P d I = A + Mu d , (9e)
∂ t m 1 E + ∇x • (m 1 u d E + m 1 P d u d ) = E + ME (9f)
where total energy and the source terms on the right-hand side are defined by

E = u 2 d 2 + Θ, (10a)
M = -R S m 0 , (10b)
A = m 0 ug -u d τ * d , (10c)
E = m 0 ug • u d -2E τ * d . ( 10d
)
where τ * d = τ d /S. Remark that Eq. (9f) represents the total energy including both the kinetic energy and the granular temperature contributions. Also remark that a separate equation for the granular temperature can be obtained after manipulating system (9):

∂ t m 1 Θ + ∇x • m 1 u d Θ + m 1 P d ∇x • u d = U + MΘ (11)
with U = - 2Θm 0 τ * d . (12)
In any case, system (9) represents a closed mesoscale description of the fuel spray that is coupled to governing equations for the gas phase.

Let us recall that assuming Θ = 0 in system (9) yields the conservation equations for the classical monokinetic EMSM model [START_REF] Kah | Taking into account polydispersity for the modeling of liquid fuel injection in internal combustion engines[END_REF] with a single unclosed term n(S = 0) in Eq. (9a). This flux term represents the disappearance of droplets due to evaporation. For an accurate evaluation of this flux, a continuous representation of n(S) must be found from data for the integer moments. Reconstructing this profile requires solving the finite Hausdorff moment problem on the interval [0, 1], i.e., finding a positive NDF n belonging to the moment sequence (m 0 , . . . , m k ) t [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF].

The existence of such an NDF is called the realizability condition. Otherwise the moments of the NDF are said to be corrupted leading to the immediate crash of the simulation. Although the set of all possible moments, called moment space, has a complex geometry in R N , a solution of this problem, preserving the moment space, has been developed in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. This solution is based on reconstructing the NDF ñ(m 0 , . . . , m k , S) using entropy maximisation of the moment sequence (m 0 , . . . , m k ) t and an associated numerical scheme for the evaporation. Using this method, Eq. (9a) is closed and the evaporative flux of a polydisperse spray can be evaluated using a kinetic scheme that preserves the moment space.

The description of the gas phase in the context of compressible, multi-species reactive flows is given by the classical Navier-Stokes conservation equations for mass, species, momentum and total energy [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]:

∂ t ρg + ∇x • ρgug = S ρ (13a)
∂ t ρgug + ∇x • (ρgug ⊗ ug + PgI) = ∇x • ρgτ + S ρu (13b)
∂ t ρgEg + ∇x • (ρgEgug + Pgug) = ∇x • ρgτ • ug -∇x • q + S ρE (13c)
where ρg is the gas density, Eg the total energy of the gas, Pg the gas pressure, τ the gas viscous-stress tensor, and q the energy flux. System (13) is closed using the equation of state for an ideal gas:

Pg = ρgRTg. ( 14
)
The viscous-stess tensor is closed in the context of the Navier-Stokes equation for Newtonian fluids using

τ = νg ∇xug + (∇xug) t - 2 3 νg (∇x • ug) I ( 15
)
where νg is the kinematic viscosity of the gas. The energy flux is given by

q = λg∇xTg + ρg ns k=1 h k D k ∇xY k ( 16
)
where ns is the number of chemical species, λg the thermal conductivity, h k the specific enthalpy associated with species k, and D k the gas mass diffusion coefficient. The internal energy evolution of the gas is governed by

∂ t ρgeg + ∇x • ρgegug = -Pg∇x • ug + ρgτ : ∇xug -∇x • q + S ρe (17)
where eg = Eg -1 2 u 2 g + Pg/ρg is the internal energy of the gas. In the literature, the tensor τ is shown to contain the gas velocity fluctuations due to droplet wakes leading to a pseudo-turbulent kinematic viscosity [START_REF] Tenneti | Direct numerical simulation of gas-solid suspensions at moderate reynolds number: Quantifying the coupling between hydrodynamic forces and particle velocity fluctuations[END_REF]. However, we assume the latter to be null in this work since we deal with a linear Stokes drag law that does not take into account wakes created by droplets. This assumption is appropriate for small droplets and is consistent with our treatment of velocity dispersion.

System [START_REF] Laurent | Multi-fluid modeling of laminar polydispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and Eq. ( 17) are written in the framework of a two-phase flow in the presence of a cloud of droplets. The disperse phase is assumed to be dilute enough such that its influence on the gas phase can be described by source terms in the gas-phase governing equations. Indeed, the kinetic model in Eq. (1) provides these source terms, which represent the variation of mass density due to evaporation, variation of momentum of the spray due to evaporation and the drag force, and the variation of the total/internal energy due to spray evaporation with uniform droplet temperature. These source terms are given, respectively, by

S ρ = 3 2 R S m 1/2 , (18a)
S ρu = -m 1/2 ug -u d τ * d + 3 2 R S m 1/2 u d , (18b)
S ρE = -m 1/2 ug • u d -2E τ * d + 3 2 R S m 1/2 E, (18c)
from which (using Eqs. (13a) and (13b)) the internal energy source term is found to be

S ρe = 2 τ * d + 3 2 R S m 1/2 Θg ( 19
)
where the gas-phase velocity dispersion is defined by

Θg = 1 2 (u d -ug) 2 + Θ. ( 20
)
Note that S ρe ≥ 0 so that internal energy is transferred from the droplet phase to the gas phase due both to a nonzero relative velocity and to velocity dispersion. Once in the gas phase, this internal energy contributes directly to increasing the gas temperature, as opposed to generating gas-phase turbulent kinetic energy (TKE). Systems ( 9) and ( 13), along with the coupling terms in [START_REF] Pope | Turbulent Flows[END_REF], provide a complete mesoscale description of the two-phase flow. In practice, the solution to these systems under ICE conditions will be highly turbulent and thus the computational cost of solving the mesoscale model will be very high. Therefore, to make the model tractable for realistic ICE conditions, it will be necessary to introduce a turbulence model, which is the topic of the next section.

3 Turbulence models for spray and gas phases

The aim of this section is to apply the Reynolds-averaging (RA) philosophy, originally introduced for the two-way coupled monodisperse flow in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF], to the two-way coupled polydisperse flow [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF]. The exact derivation of RA transport equations both for the compressible Navier-Stokes equation and the moment equations will be given below. As done in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF], let us define the following RA and phase-average (PA) quantities appearing in the RA equations.

-RA moments: m 0 , m 1 , m 2 , m 3 .

-RA number of disappearing droplets: n(t, x, 0) .

-PA velocities:

m 1 u d = m 1 u d d and ρgug = ρg ug g . -Spray-phase PA granular temperature: m 1 Θ = m 1 Θ d .
-Gas-phase PA internal energy: ρgeg = ρg eg g .

-Spray-phase PA Reynolds-stress tensor:

m 1 u ′′ d u ′′ d = m 1 u ′′ d u ′′ d d . -Gas-phase PA Reynolds-stress tensor: ρgu ′′′ g u ′′′ g = ρg u ′′′ g u ′′′ g g .
-Spray-phase PA total granular energy:

m 1 E d = m 1 E d = m 1 1 2 u d • u d d + Θ d = m 1 K + k d + Θ d (21) with K = 1 2 u d d • u d d the mean kinetic energy and k d = 1 2 u ′′ d • u ′′ d d the spray-phase TKE.
-Gas-phase PA total energy:

ρgEg = ρg Eg g = ρg 1 2 ug • ug g + eg g = ρg Kg + kg + eg g (22) with Kg = 1 2 ug g • ug g the mean kinetic energy and kg = 1 2 u ′′′ g • u ′′′
g g the gas-phase TKE.

The next step is consider the RA versions of systems ( 9) and ( 13) with appropriate turbulence closures.

Reynolds-average equations for the spray

RA moment equations

The RA moment equations are found by applying the RA to Eqs. (9a), (9b), (9c) and (9d). Let us state at the outset that the modeling of physical phenomena related to turbulent diffusive fluxes coming from size-velocity correlations is beyond the scope of this work. Thus, we assume that the following relation holds for the size moments:

m ′ k A m k = m ′ 1 A m 1 (23)
with A is any random quantity belonging to either the gas or spray phase. In words, this assumption means that whatever the index k, the correlation can expressed in terms of m 1 . This assumption implies that the size-velocity correlations leading to turbulent diffusive fluxes appearing in the RA moment equations (see Eq. ( 105) in the appendix) are null. Therefore, the RA moment equations reduce to

∂ t m 0 + ∇x • m 0 u d d = R S n(t, x, 0) , (24a)
∂ t m 1 + ∇x • m 1 u d d = M , (24b)
∂ t m 2 + ∇x • m 2 u d d = -2R S m 1 , (24c)
∂ t m 3 + ∇x • m 3 u d d = -3R S m 2 , (24d)
and thus the spatial fluxes of the RA moments depend only on the spray-phase PA velocity u d d .

Spray-phase mean momentum equation

Taking the RA of (9e) yields

∂ t m 1 u d d + ∇x • m 1 u d d ⊗ u d d + u ′′ d u ′′ d d + 2 3 Θ d I = A + Mu d . (25)
The spray-phase Reynolds-stress tensor is closed using a turbulent-viscosity model:

u ′′ d u ′′ d d = -2ν d,t S d - 1 3 ∇x • u d d I + 2 3 k d I (26)
where

S d = 1 2 ∇x u d d + (∇x u d d ) t . (27)
For later use, we also define

Sd = S d - 1 3 ∇x • u d d I. (28)
The turbulent viscosity of the spray phase is defined by

ν d,t = C d,µ k 2 d ε d . ( 29
)
Models for the spray-phase TKE k d and TKE dissipation rate ε d will be described below.

The full expression for the source terms on the right-hand side of Eq. ( 25) is given in Eqs. ( 107) and (108). Note that correlations containing the turbulent drag term (1/τ * d ) ′ can be neglected in the framework of Stokes drag with a small relative velocity between the phases. Moreover, the use of the model in ( 23) reduces the RA momentum source terms to the following expressions:

A = m 0 1 τ * d ug g -u d d + m ′ 1 u ′ g m 1 - ρ ′ g u ′ g ρg (30)
and

Mu d = -R S m 0 u d d . (31)
Remark that further simplifications can be introduced to reduce the number of unclosed correlations in (30). The term ρ ′ g u ′ g / ρg , referred to as the turbulent gas density flux, can be rewritten by adopting the low-Mach-number flow approximation, i.e., neglecting correlations between the compressibility and the turbulence in the gas phase:

α ′ d u ′ g 1 -α d = - ρ ′ g u ′ g ρg . (32)
Combining relations (30) and (32), replacing m 1 by the moment of order 3/2 through expression [START_REF]Combustion Theory (Combustion Science and Engineering Series[END_REF] to get the spray volume fraction α d , the following relation results

m ′ 1 u ′ g m 1 - ρ ′ g u ′ g ρg = α ′ d u ′ g 1 -α d α d . (33)
The term α ′ d u ′ g , referred to as the turbulent drift due to preferential segregation of droplets, can be closed through a turbulent-flux model:

α ′ d u ′ g = - ν g,t Sc g,d ∇x α d + Cg α d 1 -α d u d d -ug g (34)
where ν g,t is the turbulent viscosity of the gas, defined by

ν g,t = Cg,µ k 2 g εg , (35)
and Sc g,d the turbulent Schmidt number of the spray in the presence of the gas phase, taken as a constant in the case of small droplets but depending on the Stokes number for inertia droplets. Models for the gas-phase TKE kg and TKE dissipation rate εg will be described below. Equation (30) can then be rewritten in closed form as

A = m 0 1 τ * d ug g -u d d + uc (36)
where the drift velocity is defined by

uc = Cg u d d -ug g - ν g,t ∇x α d Sc g,d α d 1 -α d . ( 37
)
Note that since 0 ≤ Cg < 1, the mean drag is lower due to preferential concentration in a turbulent flow.

The final form for the RA spray momentum equation is

∂ t m 1 u d d + ∇x • m 1 u d d ⊗ u d d -2ν d,t Sd + 2 3 (k d + Θ d ) I = m 0 1 τ * d ug g -u d d + uc -R S m 0 u d d (38)
where we have used the equality

P d d = 2 Θ d /3.
Note that turbulence introduces the diffusive flux with coefficient ν d,t on the left-hand side, and modifies the momentum exchange terms of the right-hand side of Eq. ( 38). Likewise, in a turbulent flow, the total effective pressure in the droplet phase is proportional to

( Θ d + k d ).

Spray-phase granular temperature equation

The RA granular temperature equation is

∂ t m 1 Θ d +∇x• m 1 Θ d u d d + Θ ′′ u ′′ d d = -m 1 2 3 Θ d ∇x• u d d -m 1 P d ∇x • u ′′ d d +S Θ . (39)
where the source term is given by

S Θ = U + MΘ . ( 40
)
The term P d ∇x • u ′′ d d , which corresponds to the spray TKE dissipation, leads to

P d ∇x • u ′′ d d = -ε d (41)
where ε d is the spray-phase TKE dissipation rate. Following the modeling considerations introduced above, the source term for the granular temperature

S Θ becomes U = -2 m 0 1 τ * d Θ d (42) and MΘ = -R S m 0 Θ d . (43)
The term Θ ′′ u ′′ d d in Eq. ( 39) is the turbulent granular-temperature flux. Its closure can be achieved using a gradient-diffusion model:

Θ ′′ u ′′ d d = - ν d,t Pr d,t ∇x Θ d (44)
where Pr d,t is the turbulent Prandtl number.

The first term on the right-hand side of (39) is the production term due to mean velocity gradients whereas the second term is the production term coming from TKE dissipation. The final equation for the granular temperature is

∂ t m 1 Θ d + ∇x • m 1 Θ d u d d - ν d,t Pr d,t ∇x Θ d = -m 1 2 3 Θ d ∇x • u d d + m 1 ε d -2 m 0 1 τ * d Θ d -R S m 0 Θ d . (45)
Note that in a turbulent flow the principal production term for PA granular temperature is the one involving ε d , which increases with increasing Reynolds number.

Spray-phase total granular energy equation

The RA total granular equation found from (9f) is

∂ t m 1 E d +∇x• m 1 E d u d d + Eu ′′ d d + 2 3 Θ d u d d = -∇x• m 1 2 3 Θu ′′ d d + S E (46)
where S E is the RA source term. There are several unclosed term in this expression, such as the turbulent total granular energy flux Eu ′′ d d . Note that using the properties of PA, this term can be rewritten as

Eu ′′ d d = u ′′ d u ′′ d d • u d d + 1 2 u ′′ d u ′′ d • u ′′ d d + Θ ′′ u ′′ d d .
(47) 1 2

Terms u ′′ d u ′′ d • u ′′
u ′′ d u ′′ d • u ′′ d d + 2 3 Θu ′′ d d = - ν d,t σ d,t ∇xk d - 2ν d,t 3Pr d,t ∇x Θ d (48)
where σ d,t = 5/3 is a model constant (Rumsey 2009).

The source term contributions due to drag and evaporation on the right-hand side of Eq. (46) lead to

S E = m 0 1 τ * d u ′′′ g • u ′′ d d -2k d + U + A • u d d -R S m 0 (K + k d )+ MΘ . (49) The covariance term u ′′′ g • u ′′ d d
appearing in the drag term is important for capturing the TKE exchange between phases. In the case of point particles, Fox [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF] proposes the following model (in agreement with the DNS of Fevrier et al. [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF]):

u ′′′ g • u ′′ d d = 2β(St) (kgk d ) 1/2 (50)
where β ≈ 1 in the case of low to moderate Stokes number droplets. The final RA total granular energy equation is

∂ t m 1 E d + ∇x • m 1 E d u d d + 2 3 ( Θ d + k d ) u d d -u d d • 2ν d,t Sd = ∇x • m 1 5ν d,t 3Pr d,t ∇x Θ d + ν d,t σ d,t ∇xk d -R S m 0 (K + k d + Θ d ) + m 0 2 τ * d β (kgk d ) 1/2 -k d -Θ d + m 0 1 τ * d ug g -u d d + uc • u d d .
(51)

Note that the TKE dissipation rate does not appear on the right-hand side of this expression because the total granular energy includes the sum of the TKE and the granular energy and the role of TKE dissipation rate is to transform spray-phase TKE into Θ d .

Spray-phase mean kinetic energy

Multiplying Eq. ( 38) by the mean spray velocity u d d gives rise to the following mean kinetic energy equation after some manupulations:

∂ t m 1 K + ∇x • m 1 K u d d + 2 3 ( Θ d + k d ) u d d -u d d • 2ν d,t Sd = -Π K -R S m 0 K + m 0 1 τ * d ug g -u d d + uc • u d d (52)
where the spray-phase fluctuating energy production due to the mean spray velocity is defined by

Π K = m 1 2ν d,t Sd : Sd - 2 3 ( Θ d + k d ) ∇x • u d d . (53)
Equation ( 52) can be used to find the transport equation for the spray-phase TKE.

Spray-phase fluctuating energy

As pointed out above, both the velocity dispersion and the uncorrelated motion of droplets contribute to the fluctuations around the mean velocity. The fluctuating energy is thus the sum of the spray-phase TKE and granular energy:

κ = Θ d + k d . (54)
The expression for the fluctuating energy transport is straightforward. Subtracting (52) from (51) gives rise to the fluctuating energy equation:

∂ t m 1 κ + ∇x • m 1 κ u d d - 5ν d,t 3Pr d,t ∇x Θ d - ν d,t σ d,t ∇xk d = Π K -R S m 0 κ + m 0 2 τ * d β (kgk d ) 1/2 -κ . ( 55
)
The right-hand side of this expression has production due to the mean flow and exchange terms with the gas phase due to evaporation and drag.

Spray-phase turbulent kinetic energy

As noted above, the spray-phase fluctuating energy has two contributions. We therefore obtain the spray-phase TKE equation by subtracting (45) from (55):

∂ t m 1 k d + ∇x • m 1 k d u d d - ν d,t σ d,t ∇x ln k d - ν d,t Pr d,t M d ∇x ln Θ d = Π k d -m 1 ε d -R S m 0 k d + m 0 2 τ * d β (kgk d ) 1/2 -k d ( 56
)
where

M d = 2 Θ d 3k d , (57)
and the spray-phase TKE production due to the mean spray velocity is

Π k d = m 1 2ν d,t Sd : Sd - 2 3 k d ∇x • u d d .
(58)

Note that the spray-phase TKE dissipation rate now appears with the correct sign on the right-hand side of (56). The remaining terms represent exchanges with the gas phase due to evaporation and drag. As discussed in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF], there are no TKE production terms for the spray phase due to drag or evaporation. The ratio M d will be small at high-Reynolds numbers for small-Stokes-number droplets, and can often be neglected.

Spray-phase turbulent kinetic energy dissipation

In the model for ε d , empirical constants originated from DNS and experimental studies on classical single-phase flows need to be carefully modified in the case of two-phase flows. The reason is that the turbulence associated with each phase is produced at different integral scales. For example, the spray-phase TKE is not only produced by the gradient of the mean velocity but also through the coupling term u ′′′ g • u ′′ d d . This implies that the spray-phase TKE dissipation rate can be either higher or lower compared to the gas-phase turbulent scales. It is therefore important to choose reasonable values for the constants in the model for ε d . By analogy with Eq. (56), we model the TKE dissipation using

∂ t m 1 ε d + ∇x • m 1 ε d u d d - ν d,t σ d,ǫ ∇x ln ε d - ν d,t Pr d,t M d ∇x Θ d = C 1 d,ǫ ε d k d Π k d -C 2 d,ǫ m 1 ε 2 d k d -C 5 d,ǫ R S m 0 ε d + C 3 d,ǫ m 0 2 τ * d βǫ (εgε d ) 1/2 -ε d ( 59
)
where

C 1 d,ǫ , C 2 d,ǫ , C 3 d,ǫ , C 5 d,ǫ
, βǫ and σ d,ǫ are constants. Note that in this model, the evaporation and drag terms are written in terms of the TKE dissipation rates, and not as source terms. This modification from the classical two-phase turbulence model leads to a more robust formulation [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF].

Reynolds-average equations for gas phase

The RA source terms seen by the gas phase contain numerous correlations that can lead to confusion when deriving RA gas equations. We thus prefer to work with some new expressions for these source terms given in [START_REF] Pope | Turbulent Flows[END_REF]. Let us first define

Mg = 3 2 R S m 1/2 , Ag = -m 1/2 ug -u d τ * d , (60)
and then we rewrite the gas-phase governing equations using these definitions.

Given the assumptions made concerning the velocity-size correlations in [START_REF]Combustion Theory (Combustion Science and Engineering Series[END_REF], working with (60) is entirely equivalent to working with the original moments.

Gas-phase continuity equation

Applying the RA to Eq. (13a) yields

∂ t ρg + ∇x • ρg ug g = Mg . ( 61
)
This equation is closed.

Gas-phase momentum equation

The gas-phase PA momentum equation is found from (13b):

∂ t ρg ug g +∇x• ρg ug g ⊗ ug g + u ′′′ g u ′′′ g g + Pg I = ∇x• ρg τ g + S ρu (62) with S ρu = Mgu d + Ag . (63)
The RA viscous-stress tensor is

τ g = νg g ∇x ug g + ∇x ug g t - 2 3 ∇x • ug g I + ν ′′′ g ∇xu ′′′ g + ∇xu ′′′ g t - 2 3 ν ′′′ g ∇x • u ′′′ g I g ( 64
)
where the second contribution is neglected in compressible turbulence models. For the sake of simplicity, we define Sg = Sg-1 3 (∇x• ug g )I and Sg = 1 2 ∇x ug g + (∇x ug g ) t . In the literature, it is common to combine the RA viscous-stress tensor with the gas-phase Reynolds-stress tensor u ′′′ g u ′′′ g g , for which closure is achieved through a turbulent-viscosity model. Therefore, it is convenient to express the total stress tensor as

u ′′′ g u ′′′ g g -τ g = -2 νg g + ν g,t Sg + 2 3 kgI. ( 65
)
The source terms accounting for drag and evaporation are, respectively,

Ag = -m 1/2 1 τ * d ug g -u d d + uc (66)
and

Mgu d = R S m 1/2 u d d . (67)
The final form for the gas-phase RA momentum equation is

∂ t ρg ug g +∇x• ρg ug g ⊗ ug g + 2 3 kg + R Tg g I -2 νg g + ν g,t Sg = 3 2 R S m 1/2 u d d -m 1/2 1 τ * d ug g -u d d + uc . ( 68
)
Note that the exchange terms on the right-hand side have the opposite sign as the corresponding terms in the spray-phase momentum equation. By construction, the total mean momentum of the two-phase system is conserved.

Gas-phase total energy equation

Reynolds averaging the total energy equation for the gas phase (13c) yields the following expression:

∂ t ρg Eg g + ∇x • ρg Eg g ug g + Egu ′′′ g g + Pg ug g + Pgu ′′′ g = + ∇x • ρg τ g ug g + τ ′′′ u ′′′ g g -∇x • q + S ρE . ( 69
)
The turbulent total energy flux can be decomposed as follows:

Egu ′′′ g g = ug g • u ′′′ g u ′′′ g g + 1 2 u ′′′ g u ′′′ g • u ′′′ g g + e ′′′ g u ′′′ g g - Pgu ′′′ g ρg . (70)
Consistent with the velocity flux, the turbulent fluxes can be grouped together and modeled as 1 2

u ′′′ g u ′′′ g • u ′′′ g g -τ ′′′ u ′′′ g g = -νg + ν g,t σ g,t ∇xkg (71)
where σ g,t is a constant whose value depends on the turbulence model. The final expression for the gas-phase PA total energy is

∂ t ρg Eg g + ∇x • ρg Eg g ug g - ν g,t
Pr g,t ∇x eg g -νg + ν g,t σ g,t ∇xkg

+ ∇x • ρg R Tg g + 2 3 kg ug g -2 ν g,t + νg g Sg • ug g + ∇x • q = 3 2 R S m 1/2 (K + k d + Θ d ) -m 1/2 2 τ * d β (kgk d ) 1/2 -k d -Θ d -m 1/2 1 τ * d ug g -u d d + uc • u d d . ( 72
)
As with the mean momentum, the total energy of the two-phase system is conserved by the exchange terms.

Gas-phase internal energy equation

Taking the RA of Eq. ( 17) leads to

∂ t ρg eg g +∇x• ρg eg g ug g + e ′′′ g u ′′′ g g = -Pg ∇x• ug g - ρ ′ g u ′ g ρg -P ′ g ∇x • u ′ g + ρg τ g : ∇x ug g + τ ′′′ : ∇xu ′′′ g g -∇x • q + S ρe . ( 73
)
The internal energy turbulent flux is closed through a gradient-diffusion model:

e ′′′ g u ′′′ g g = - ν g,t Pr g,t ∇x eg g (74)
where Pr g,t is the gas-phase turbulent Prandtl number. The gas-phase TKE dissipation rate produces a source term in the internal energy equation. This is due to P ′ g ∇x • u ′ g -ρg τ ′′′ : ∇xu ′′′ g g whose closure is achieved by the following equality

P ′ g ∇x • u ′ g -ρg τ ′′′ : ∇xu ′′′ g g = -ρg εg ( 75
)
where εg is the TKE dissipation rate in the gas phase. The RA source term is

S ρe = 2 τ * d + 3 2 R S m 1/2 K -ug g • u d d + Kg + k d -u ′′′ g • u ′′ d d + kg + Θ d - 2 τ * d + 3 2 R S m 1/2 α ′ d u ′′′ g α d 1 -α d • u d d -ug g - 1 2 α ′ d u ′′′ g • u ′′′ g α d 1 -α d . (76)
The triple correlation α ′ d u ′′′ g • u ′′′ g in (76) represents correlations between the spray volume fraction and the gas-phase TKE. This term is likely to be negative (i.e., the more inertial the droplets, the faster they segregate away from regions of high gas vorticity), but we will assume it is small and neglect it.

Introducing the turbulence models gives rise to final gas-phase RA internal energy equation:

∂ t ρg eg g + ∇x • ρg eg g ug g - ν g,t
Pr g,t ∇x eg g + q ρg

= Π e g + ρg εg

+ 2 τ * d + 3 2 R S m 1/2 K -ug g • u d d + Kg + k d -2β (kgk d ) 1/2 + kg + Θ d - 2 τ * d + 3 2 R S m 1/2 uc • u d d -ug g (77)
where

Π e g = -Πe g -ρg R Tg g ∇x • ug g + 2 ρg νg g Sg : Sg (78)
is the production of gas-phase internal energy through the mean velocity, and

Πe g = -ρg R Tg g ∇x • α d uc (79)
is the contribution due to the drift velocity generated by preferential concentration. Note that the contribution of the mean velocities in the exchange term can be rewritten as

K -ug g • u d d + Kg = 1 2 ug g -u d d 2 , (80)
which is strictly positive. Likewise, the TKE contributions can be rewritten as (k

1/2 d -k 1/2 d ) 2 + 2(1 -β)(k d k d ) 1/2 with β ≤ 1.
Thus, the exchange terms act to increase the gas-phase PA internal energy.

Gas-phase mean kinetic energy equation

The gas-phase mean kinetic energy is known from the gas-phase PA velocity. It can be found by multiplying the gas-phase PA momentum equation (68) by ug g :

∂ t ρg Kg+∇x• ρg Kg + R Tg g + 2 3 kg ug g -2 ν g,t + νg g Sg • ug g = ρg R Tg g + 2 3 kg I -2 ν g,t + νg g Sg : ∇x ug g -m 1/2 1 τ * d ug g -u d d + uc • ug g + 3 2 R S m 1/2 u d d • ug g -Kg .
(81)

This expression can now be used to find the equation for the gas-phase TKE.

Gas-phase turbulent kinetic energy equation

Using the equality Eg g = eg g + Kg + kg, the gas-phase TKE can be easily found by subtracting the gas-phase internal energy equation ( 77) and the gas-phase mean kinetic energy equation (81) from the total energy equation (72):

∂ t ρg kg + ∇x • ρg kg ug g -νg + ν g,t σ g,t ∇x ln kg = Π k g + Πe g +Π k c -ρg εg+ m 1/2 2 τ * d + 3R S β (kgk d ) 1/2 -kg + 3 2 R S m 1/2 kg (82)
where

Π k g = ρg ν g,t Sg : Sg - 2 3 ρg kg∇x • ug g (83)
is the gas-phase TKE production due to mean velocity gradients, and

Π k c = m 1/2 1 τ * d + 3 2 R S max 0, uc • u d d -ug g (84)
is a TKE production term due to two-way coupling and the preferential concentration of droplets in a turbulent flow. Note that the coupling terms in (82) are asymmetric with respect to those in (56) for the spray phase. This is caused by the preferential concentration of droplets: under such circumstances the gas seen by droplets has different statistics than the gas itself [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF]. Nevertheless, in shear-driven flows such as fuel jets in ICE injectors, the main TKE production term will be Π k g .

Gas-phase turbulent kinetic energy dissipation equation

By analogy with Eq. ( 82), the gas-phase TKE dissipation rate is found from

∂ t ρg εg + ∇x • ρg εg ug g -νg + ν g,t σg,ǫ ∇x ln εg = εg kg C 1 g,ǫ Π k g + C 6 g,ǫ Πe g + ε d k d C 4 g,ǫ Π k c -C 2 g,ǫ ρg ε 2 g kg + C 3 g,ǫ m 1/2 2 τ * d + 3R S βǫ (εgε d ) 1/2 -εg + C 5 g,ǫ 3 2 R S m 1/2 εg (85)
where σg,ǫ, C 1 g,ǫ , C 2 g,ǫ , C 3 g,ǫ , C 4 g,ǫ , C 5 g,ǫ and C 6 g,ǫ are model constants.

Final remarks

The variables and the corresponding transport equations solved for each phase in the CFD code are as follows.

-Spray phase: As in the case of laminar flow, the RA moments are used to reconstruct the RA droplet-size NDF n(t, x, S) , which is needed to close the phase-space flux term in (24a) and to find other moments such as m 1/2 . In the turbulence models, the model constants are set at the standard values used in free-shear flows such as turbulent jets.

4 Homogeneous turbulence of two-phase polydisperse flows

The this section we apply the multiphase turbulence model derived in Sec. 2 to a simple test case, namely, homogeneous two-phase turbulence. The most important point is to capture the correct behavior of the total energy partition between the evaporating spray and the compressible gas in the context of two-way coupling.

In [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF] the effectiveness of the closures proposed for the turbulent velocity correlations between phases u ′′′ g • u ′′ d d = 2β(kgk d ) 1/2 has been shown in the context of one-way and two-way coupling of collisionless monodisperse droplets with an incompressible gas. The model showed good agreement with the total fluctuating energy expression first given in [START_REF] Tchen | Mean value and correlation problems connected with the motion of small particles suspended in a turbulent fluid[END_REF], and provided results comparable with [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF] for the partition between the spray-phase TKE k d and the granular temperature Θ d . Here we focus on the extension of these results to a more complicated framework with a polydisperse spray and two-way coupling. In order to simplify the analysis of the different flow regimes, we first present a dimensionless homogeneous system of equations. We then focus on selected test cases, which highlight the turbulence model predictions for different flow conditions relevant to ICE sprays. In the first case, we are interested in the energy partition for a polydisperse, nonevaporating spray with frozen gas turbulence, similar to [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF][START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF]. The second case considers two-way coupling effects on the energy partition for both evaporating and non-evaporating sprays. After these preliminary studies, the third case considers conditions relevant to direct-injection ICE applications.

Dimensionless equations for homogeneous turbulent flow

A summary of the equations needed for the two-phase turbulence model with polydisperse evaporating sprays is given in Sec. 3.3. By neglecting all terms representing spatial transport, a simplified set of equation suitable for time-evolving homogeneous flow is found. We can define the following reference quantities: turbulent kinetic energy k∞, material density of the liquid inside the droplets ρ l,∞ , material density of the gas ρ l,∞ , kinematic viscosity of the gas νg,∞, maximum droplet size S 0 , number per volume of droplets n 0 , and the integral time scale of the gas-phase turbulence τg = k∞/εg(0) where εg(0) is the value the gas-phase TKE dissipation rate at time zero. The dimensionless quantities appearing in the homogeneous model are as follows:

kg = kg k∞ , k d = k d k∞ , εg = εgτg k∞ , ε d = ε d τg k∞ , ρ g = ρg ρg,∞ , x = x √ k∞τg , t = t τg , ug = ug g √ k∞ , u d = u d d √ k∞ , S = S S 0 , R S = R S τg S 0 , νg = νg νg,∞ , Sd =τg Sd , Sg =τg Sg, eg = eg g k∞ , ν g,t = ν g,t νg,∞Re , ν d,t = ν d,t νg,∞Re , ρ d = ρ d ρ d,∞ , m k = m k S k 0 n 0 , n = n S 0 n 0
(86) where Re = k∞τg/νg,∞ is the turbulence Reynolds number [START_REF] Pope | Turbulent Flows[END_REF].

In the following, we work with dimensionless RA equations, omitting the bars on variables of the mean two-phase flow as well as the RA and PA brackets (e.g., m k = m k ) for clarity. The dimensionless RA moment equations are given by

d t m 0 = -n(t, x, 0)R S , (87a)
d t m 1 = -m 0 R S , (87b)
d t m 2 = -m 1 R S , (87c)
d t m 3 = -m 2 R S . ( 87d
)
The dimensionless continuity equation for the gas phase is

d t ρg = 3 2 Φmm 1/2 R S ( 88
)
with the dimensionless mass loading Φm = ρ d,∞ S 3/2 0 n 0 /ρg,∞. The dimensionless RA mean momentum equations are

d t Φmm 3/2 u d = Φmm 1/2 St (1 -Cg) (ug -u d ) - 3 2 Φmm 1/2 R S u d , (89a)
d t ρgug = - Φmm 1/2 St (1 -Cg) (ug -u d ) + 3 2 Φmm 1/2 R S u d ( 89b
)
where St is the Stokes number relative to the droplet size S 0 . The reader can note that the total mean momentum of the two-phase system is constant. The dimensionless RA internal-energy equation for the gas phase is

d t ρgeg = Π e g Re +ρgεg+ 2 St + 3 2 R S Φmm 1/2 1 2 -Cg (ug -u d ) 2 + kg -2β (kgk d ) 1/2 + κ ( 90
)
where the production term due to molecular dissipation of mean shear Π e g will be non-zero for free-shear flow but otherwise is null. Note that the right-hand side of this equation should always non-negative since the gas-phase internal energy never decreases, which implies that Cg ≤ 1/2.

The dimensionless system of conserved quantities representing fluctuations in the spray phase is

d t m 3/2 κ =Π K + 2 St m 1/2 β (kgk d ) 1/2 -κ - 3 2 R S m 1/2 κ (91a) d t m 3/2 k d =Π k d -m 3/2 ε d + 2 St m 1/2 β (kgk d ) 1/2 -k d - 3 2 R S m 1/2 k d (91b) d t m 3/2 ε d =C 1 d,ǫ ε d k d Π k d -C 2 d,ǫ m 3/2 ε 2 d k d + C 3 d,ǫ 2 St m 1/2 βǫ (εgε d ) 1/2 -ε d -C 5 d,ǫ 3 2 R S m 1/2 ε d (91c)
where Π K = Π k d are the dimensionless production terms due to turbulent dissipation of mean shear in the spray phase. The dimensionless system of conserved quantities representing the turbulent fluctuations in the gas phase is

d t ρgkg =Π k g -ρgεg + 2 St + 3R S Φmm 1/2 β (kgk d ) 1/2 -kg + 1 St + 3 2 R S Φmm 1/2 Cg (ug -u d ) 2 + 3 2 R S Φmm 1/2 kg (92a) d t ρgεg =C 1 g,ǫ εg kg Π k g -C 2 g,ǫ ρg ε 2 g kg + C 3 g,ǫ 2 St + 3R S Φmm 1/2 βǫ (εgε d ) 1/2 -εg + C 4 g,ǫ ε d k d 1 St + 3 2 R S Φmm 1/2 Cg (ug -u d ) 2 + C 5 g,ǫ 3 2 R S Φmm 1/2 εg (92b)
where Π k g is the dimensionless production term due to turbulent dissipation of mean shear in the gas phase. In the following, we consider

C 1 g,ǫ = C 1 d,ǫ = 1.44, C 2 g,ǫ = C 2 d,ǫ = 1.92, C 3 g,ǫ = C 3 d,ǫ = 1.55, C 4 g,ǫ = 2.11, C 5 g,ǫ = C 5 d,ǫ = 1.
55, and Cg = 0.3. The values of the remaining production terms are set equal to zero unless noted otherwise.

Example results for fluctuating energy partition with decaying turbulence

In order to validate the spray-phase energy partition for polydisperse droplets, we use the three one-way-coupling cases in [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF], originally conducted for monodisperse droplets. These cases use three sets of initial conditions: (i) κ(0) = 1, k d (0) = 1; (ii) κ(0) = 0, k d (0) = 0; and (iii) κ(0) = 0.83, k d (0) = 0. The gas-phase turbulence is frozen and set as kg = 1. The mean velocities are equal ug = u d , and remain so for all time due to the abscence of gravity. For the case with initial value k d (0) = 1, we use the initial value ε d (0) = 2. Otherwise we use ε d (0) = 0 when k d (0) = 0 for the sake of consistency. For monodisperse droplets in [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF], a Stokes number equal to Stm = 0.81 is employed. This value corresponds to the Stokes number to the time scale relative to the mean droplet size of the polydisperse distribution over the (constant) gas integral time scale τg. For convenience, we set τg = 1, leading to a gas-phase turbulence energy dissipation of εg = 1.

For the following simulations, the values for β and βǫ are fixed to 0.8, as suggested in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF], since this provides the correct steady-state results for κ as compared to [START_REF] Tchen | Mean value and correlation problems connected with the motion of small particles suspended in a turbulent fluid[END_REF] and the correct energy partition as found in [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF]. Moreover, the following Rosin-Rammler function is used to characterize the droplet size distribution in the spray:

n (S) = 1 2 qrr16 qrr/2 S qrr 2 -1 exp -(16S) qrr/2 (93)
with the constant qrr = 3.5 determining the sharpness of the distribution and the dimensionless size phase variable S ∈ [0, 1]. According to this distribution, the largest droplet size is equal to S 0 = 16Smean where Smean is the mean droplet size. This yields a characteristic Stokes number of St = 16Stm.

Recall that for one-way coupling of polydisperse droplets with a frozen gas phase, only the system (91) with Π K = Π k d = 0 is solved. Results are plotted in Fig. 1. The spray energy partition is a good agreement with Fig. 3 in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF], and the small differences in the steady-state values are due to the polydispersity of the droplets in Fig. 1. In summary, the turbulence model proposed in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF] has been successfully extended to the Eulerian high-order moment method for polydisperse droplets.

Two-way coupled interactions require the complete turbulence model given by ( 91) and (92) for a non-evaporating spray and (91), and (92) and (87) for an evaporating spray. As with one-way coupling, production terms involving mean gradients are set to zero, i.e.,

Π K = Π k d = Π k g = 0.
The mass loading number is set to Φm = 0.4 for both non-evaporating and evaporating cases, and the evaporation is taken as R S = 10St for the evaporating case. Results for a non-evaporating spray are plotted in Figs. 2 and3. In contrast to the spray energy partition observed for one-way coupling, the normalized values in Fig. 2 show that the granular energy decreases towards zero while the sprayphase TKE moves towards the total fluctuating energy. Moreover, the gas-phase turbulence dissipates into gas-phase internal energy since no production term due to the mean flow are considered. Furthermore, the spray mass is not high enough to generate significant turbulence inside the gas phase. However, let us underline the small effect of initial spray fluctuating energy on the gas phase as observed in the internal energy profiles, which are not the same for different spray initial conditions and increase with κ. 1 Dynamics of the dimensionless non-evaporating spray-phase energy components with frozen gas-phase turbulence (one-way coupling). Curves correspond to three initial conditions: (i) solid lines, (ii) dashed lines, and (iii) stars. Curves corresponding to the fluctuating energy κ, the spray-phase kinetic energy k d and the granular temperature Θ are respectively denoted through circle, square and triangle symbols. 2 Dynamics of the normalized non-evaporating spray-phase energy components with twoway coupling. Curves correspond to three initial conditions: (i) solid black lines, (ii) dashed blue lines, (iii) red dots. Curves corresponding to the fluctuating energy κ, the spray-phase kinetic energy k d and the granular temperature Θ are respectively denoted through circle, square and triangle symbols.

0 2 4 6 8 t / τ g 0 0,5 1 κ / k g , k d / k g , Θ d / k g Fig.
Results for an evaporating spray are plotted in Figs. 4, 5 and 6. In this case, time t = 8 corresponds to the characteristic evaporation time for the mean droplet radius of the distribution. The spray-phase energy partition displayed in Fig. 5 is very similar to the case with a non-evaporating spray despite some small differences. When it comes to the gas phase, although the internal energy profiles do not change much, the final values decrease for each initial condition, since the spray evaporation heats the gas phase. Moreover, it is observed that the gas-phase TKE remains the same as compared to the non-evaporating case. This proves that the 3 Dynamics of the gas-phase energy components for non-evaporating droplets with twoway coupling. Curves correspond to gas-phase TKE (top) and gas-phase internal energy (bottom) for three initial conditions: (i) solid black lines, (ii) dashed blue lines, (iii) red dots. turbulence generation in the gas due to the droplet presence is going down since droplets are disappearing but the evaporation creates some turbulence in the gas, counterbalancing the former.

In summary, the test cases without turbulence source terms due to mean velocity gradients but with different initial energy distributions illustrate that the two-phase turbulence model with evaporating droplets behaves as might be expected from the behavior observed with frozen gas-phase turbulence. In general, the turbulence decay observed for the gas phase is only slightly altered by the presence of the spray. Nevertheless, the numerical method used to reconstruct the NDF from the RA moments works equally well as it did for the laminar cases in our earlier work. 6 Dynamics of the gas-phase energy components with evaporating droplets. Curves correspond to gas-phase TKE (top) and gas-phase internal energy (bottom) for three initial conditions: (i) solid black lines, (ii) dashed blue lines, (iii) red dots.

Example results for typical ICE flow conditions

The complete turbulence model will eventually be implemented in an industrial CFD code [START_REF] Bohbot | IFP-C3D: an unstructured parallel solver for reactive compressible gas flow with spray[END_REF] and validated for high-pressure direct injection simulations of ICE applications. The scope of this section is to verify the feasibility of using the turbulence model under typical ICE operating conditions [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF]. In the combustion chamber of a typical diesel engine, the gas-phase turbulence is mainly produced by the high-speed direct injection of the fuel spray. Due to the significant velocity difference between the gas and the injected liquid fuel, the spray phase is nearly laminar near the nozzle and the associated spray-phase TKE is very small. However, the droplet flow regime, having strong interactions with the turbulent gas phase, becomes turbulent farther from the injector nozzle. Therefore, we will focus on a two-phase flow regime composed of a fully-developed gas-phase turbulent jet interacting with a nearly laminar spray phase composed of polydisperse droplets. Injection velocities for each phase and other relevant quantities needed to initialize the time-evolving simulations are representative of ICE conditions.

To obtain the characteristic initial data for ICE conditions, an injection test case, representing Spray-H conditions in [START_REF]Engine combustion network data archive[END_REF], has been computed using the IFP-C3D software [START_REF] Bohbot | IFP-C3D: an unstructured parallel solver for reactive compressible gas flow with spray[END_REF], adopting a Lagrangian description of polydisperse droplets [START_REF] O'rourke | Collective drop effects on vaporizing liquid sprays[END_REF]. The relevant characteristic values are the gas-phase TKE kg, gas-phase TKE dissipation εg, the spray-phase volume fraction α d , the mean droplet radius rsmr, characteristic evaporation time scale τ S,smr and the drag time scale τ d,smr associated with the spray. These are obtained from the IFP-C3D simulation results and are used as initial conditions for our time-evolving homogeneous simulation. Table 1 provides the relevant data. These values are employed to compute the dimensionless numbers given in Table 2 characterizing the disperse-flow regime. Note that the drag and integral time scales are comparable, while the evaporation time scale is two orders of magnitude larger. Also note that the Stokes number St is nearly unity, indicating that particle trajectory crossing (represented by Θ in our model) will be important at the turbulence integral scale. In a self-similar turbulent jet, the turbulence Reynolds number (Re = k 2 g /(νgεg)) is constant and τg increases along the centerline (i.e., with time in our model). Thus, with one-way coupling, the Stokes numbers would decrease in the non-evaporating case due to the change in τg. Obviously, the decrease in Stokes number will be more rapid with evaporating droplets due to the decrease in τ d,smr .

The conceptual difference between the Lagrangian description of the droplets used to characterize the turbulent disperse-phase flow and the Eulerian turbulent model developed in this work do not allow us to extract directly the values for spray-phase TKE k d and its dissipation rate ε d from the IFP-C3D simulation results. Nevertheless, the integral time scale associated with the spray-phase TKE should be comparable to the gas-phase integral time scale τg in the injection zone where the velocity difference between phases is significant. Therefore, we set k d = κ = 1 and ε d = 10 as the initial state of the spray-phase turbulence. Note that these values are much smaller than in the gas phase, and correspond to nearly laminar flow. The initial polydisperse droplet size distribution follows the function defined by Eq. ( 93).

Although we are considering a homogeneous case, energy production terms due to mean velocity gradients in the gas phase need to be taken into consideration since the direct injection process is modeled as a free-shear flow for the gas-phase turbulence [START_REF] Pope | Turbulent Flows[END_REF]. It is therefore appropriate to provide closures for these production terms, considering a constant turbulence Reynolds number as is observed for singlephase turbulent jets [START_REF] Pope | Turbulent Flows[END_REF]. For this purpose, we deduce an analytical expression for Π k g in Eq. (92) and Eq. ( 91), respectively, corresponding to a self-similar turbulent jet:

Π k g = 2 -C 2 g,ǫ 2 -C 1 g,ǫ ρgεg, Π k d = 2 -C 2 d,ǫ 2 -C 1 d,ǫ m 3/2 ε d . (94)
In the absence of two-way coupling, this expression for Π k g will yield a constant turbulence Reynolds number as expected for a turbulent jet [START_REF] Pope | Turbulent Flows[END_REF]. However, due two-way coupling, the turbulence Reynolds number under ICE conditions will vary with time (or distance from the nozzle). The model for Π k d is taken in analogy to the gas phase.

In the following, three time-dependent homogeneous cases are considered, based on the following initial conditions: (i) single-phase gas flow, (ii) two-phase flow of a compressible gas and a non-evaporating polydisperse spray and (iii) two-phase flow composed of a compressible gas and an evaporating polydisperse spray. Case (i) is a reference case for a self-similar turbulent jet where the data corresponds to the centerline turbulence statistics moving with the mean velocity. Case (ii) uses the same initial conditions as case (i) but with R S = 0, and thus illustrates the effect of two-way coupling due to drag on the turbulence statistics. Finally, case (iii) shows the additional effect of droplet evaporation on two-way coupling. For reference, the mean velocities and the spray density concentration are shown in Fig. 7 as a function of dimensionless time. (In all figures, the time is made dimensionless by the initial gas-phase integral time scale.) Note that the mean momentum equations (see, e.g., Eq. (89)) do not contain terms corresponding to the spread of the turbulent jet, which would lower the mean velocities. Thus, the mean velocities plotted in Fig. 7 can be thought of as normalized values relative to the centerline jet velocity. In any case, from Fig. 7 we can observe that for the two spray cases the mean velocities are nearly equal for t ≥ 4τg. The fact that the two spray cases have nearly the same velocity profiles is due to the relatively slow evaporation rate for case (iii) as seen from the slow decrease in the spray density concentration.

Figures 8 and9 present the evolution of the turbulence Reynolds number and gas-phase energies, respectively, for each case. As expected due to the modeling assumptions for Π k g and Π k d , in case (i) the turbulence Reynolds number remains constant although the gas-phase TKE decreases with time. For cases (ii) and (iii), the turbulence Reynolds number first increases and decreases sharply at the beginning of the simulation and then slowly decreases. This behavior comes from the fact that the presence of droplets increases significantly the gas-phase TKE kg due to two-way coupling (i.e., the source terms in Eq. (92) involving Cg). It can also be remarked that the gas-phase internal energy increases up to ∼ 8 from 0.5 for cases (ii) and (iii), whereas it remains quite low for case (i). This difference comes from the two-way coupling source terms, especially due to the mean velocity difference between the two phases shown in Fig. 7, which has a significant effect on the gas-phase turbulence.

In Fig. 10, the spray-phase energy partition and the TKE dissipation rates are given for cases (ii) and (iii). It can be observed that the complete relaxation between the gas and spray does not take place since the the normalized spray-phase TKE k d remains below one. Nonetheless, the energy partition seems reasonable since the granular temperature Θ tends to diminish while the spray-phase TKE k d and the total fluctuating energy κ become closer and approach kg with time. Note the differences between cases (ii) and (iii) are minimal due to the relatively slow evaporation of droplets in case (iii). Overall, the spray-phase turbulence is driven by the energy exchange terms from the gas phase due to drag, although as seen in Fig. 9 the presence of droplets significantly increases the gas-phase turbulence and, hence, leads to turbulence in the spray phase.

Conclusions

Within the context of two-phase turbulence modeling for ICE applications, the main purpose of the present contribution is to apply the Reynolds-averaging (RA) philosophy, originally introduced for two-way coupled monodisperse flows in [START_REF] Fox | On multiphase turbulence models for collisional fluid-particle flows[END_REF], to two-way coupled polydisperse flows with evaporating droplets and a compressible gas [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF]. For clarity, we have considered only the simplest forms for drag model through Stokes law and d 2 -constant evaporation law for the droplets. However, the generalization to more complex models should be straightforward.

In the context of homogeneous two-phase turbulence, the new model has been investigated for polydispere droplets with both one-way and two-way coupling. The realizability condition of EMSM method has been successfully respected, validating the adequacy of the new turbulence modeling approach with high-order moment methods in the presence of spray evaporation. Moreover, it has been qualitatively validated as compared to work done in [START_REF] Fevrier | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF] within the framework of one-way coupling. The correct behavior of the energy partition in the two-phase flow has been also observed for two-way coupled evaporating and non-evaporating sprays. As far as industrial applications are concerned, the model has been studied under realistic characteristic scales and values provided from the data of 3D injection simulations (spray-h conditions [START_REF]Engine combustion network data archive[END_REF]) in the industrial CFD code IFP-C3D [START_REF] Bohbot | IFP-C3D: an unstructured parallel solver for reactive compressible gas flow with spray[END_REF]. The latter work gave some significant insights on the underlying physics although it lacks of quantitative validations.

Given these achievements, the model is now ready to be implemented in an industrial software (i.e. IFP-C3D [START_REF] Bohbot | IFP-C3D: an unstructured parallel solver for reactive compressible gas flow with spray[END_REF]). However, further effort is required for its numerical resolution using the ALE formalism for moving geometries. Let us recall that for laminar two-phase flows, the two-way coupled EMSM method had been successfully implemented in [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF] using a pressure-less gas formalism for its resolution. However, the fact that there are many gradients arising from the turbulence model renders the task more complicated. As soon as the numerical implementation is complete, the turbulence model should be tested under realistic injection conditions and quantitatively validated against experimental data.

A Reynolds and Phase averaging

A.1 Reynolds average

The RA of a product can be decomposed:

AB = (A ′ + A )(B ′ + B ) = A ′ B ′ + A ′ B + B ′ A + A B = A ′ B ′ + A ′ B + B ′ A + A B = A ′ B ′ + A B (95)
and

ABC = (A ′ B ′ + A ′ B + B ′ A + A B )(C ′ + C ) = A ′ B ′ C ′ + A ′ C ′ B + B ′ C ′ A + C ′ A B + A ′ B ′ C + A ′ B C + B ′ A C + A B C = A ′ B ′ C ′ + A ′ C ′ B + B ′ C ′ A + C ′ A B + A ′ B ′ C + A ′ B C + B ′ A C + A B C = A B C + A B ′ C ′ + B A ′ C ′ + C A ′ B ′ + A ′ B ′ C ′ (96)

A.2 Phase average

In the context of dilute disperse flows where the gas phase is modeled through compressible Navier-Stokes equation and polydisperse droplets by a moment method, both the gas density ρg and the spray density concentration denoted by the moment m 3/2 play an important role for the averaging procedure. Therefore, the relation between RA and the PA of an arbitrary variable A for the gas phase is ρgA = ρg A g ,

the polydispersed phase is

m 1 A = m 1 A d . (98)
A useful identity relating the PA of the quantity A to its covariances with respect to the gas density and the moment m 1 can be also derived. Adding Eq. ( 97) to (98), one can obtain the following relation:

A d = A g + ρg + m 1 A -ρg + m 1 A g m 1 (99)
multiplying by ρg both the numerator and the denominator of the second term at the righthand side in Eq. (99) and using the relation, one can easily obtain the following expression:

A d = A g + m ′ 1 A ′ m 1 - ρ ′ g A ′ ρg (100)
Deriving Reynolds averaged moment equations requires a density-weighted statistics of terms m k A where m k is the moment of order k with k is different from 3/2. This accounts for leading the following operation:

m k A = m 1 m k A m 1 = m 1 A m k + A ′ m ′ k m 1 = m 1 A -m ′ 1 A ′ m 1 m k + m ′ k A ′ (101) therefore, m k A = m k A d - m k m 1 m ′ 1 A ′ + m ′ k A ′ (102)
An other difficulty encountered during the derivation phase of RA equations is observed when the RA of a variable A is obtained. The latter need to be defined in function of a densityweighted statistics. This can be done either with respect to the gas phase,

A = A g - ρ ′ g A ′ ρg (103)
or the spray phase,

A = A d - m ′ 1 A ′ m 1 . (104)

B Full RA equations

The aim of this part is to give further details in the derivation of full coupled RA two-phase flow equations. Some source terms contain third-order correlations, which need to be carefully explicated for modeling issues. Let us then focus on each RA equation in the reminder. The source term at the right hand side of Eq. ( 106) still needs to be developed. Through the formulas (102) and (96) given in A, the mean source term due to drag force can be developed as

B.1 RA moment equations

∂t m 0 + ∇x • m 0 u d d - m 0 m 1 m ′ 1 u ′ d + m ′ 0 u ′ d = R S n(t, x, 0) , (105a)
∂t m 1 + ∇x • m 1 u d d = M , ( 105b
)
∂t m 2 + ∇x • m 2 u d d - m 2 m 1 m ′ 1 u ′ d + m ′ 2 u ′ d = -2R S m 1 , (105c)
∂t m 3 + ∇x • m 3 u d d - m 3 m 1 m ′ 1 u ′ d + m ′ 3 u ′ d = -3R S m 2 . (105d)

B.2 RA spray momentum equation

A = m 0 1 τ * d + m ′ 0 ( 1 τ * d ) ′ ug g -u d d + m ′ 1 u ′ d m 1 - ρ ′ g u ′ g ρg + 1 τ * d m ′ 0 u ′ g -m ′ 0 u ′ d + m 0 u ′ g ( 1 τ * d ) ′ -u ′ d ( 1 τ * d ) ′ + m ′ 0 u ′ g ( 1 τ * d ) ′ -m ′ 0 u ′ d ( 1 τ * d ) ′ , (107)
whereas the mean momentum exchange created by the evaporating spray is given through the expression

Mu d = -R S m 0 u d d - m 0 m ′ 1 u ′ d m 1 + m ′ 0 u ′ d . (108)

  d d and Θu ′′ d d together account for the spray-phase energy flux whose closure model needs to be consistent with the spray-phase Reynolds-stress tensor:

  Fig.1Dynamics of the dimensionless non-evaporating spray-phase energy components with frozen gas-phase turbulence (one-way coupling). Curves correspond to three initial conditions: (i) solid lines, (ii) dashed lines, and (iii) stars. Curves corresponding to the fluctuating energy κ, the spray-phase kinetic energy k d and the granular temperature Θ are respectively denoted through circle, square and triangle symbols.

  Fig.3Dynamics of the gas-phase energy components for non-evaporating droplets with twoway coupling. Curves correspond to gas-phase TKE (top) and gas-phase internal energy (bottom) for three initial conditions: (i) solid black lines, (ii) dashed blue lines, (iii) red dots.

Fig. 4

 4 Fig.4Gas-phase density evolution (top) and spray density concentration (bottom) with polydisperse evaporating droplets.

Fig. 5

 5 Fig.5Dynamics of the normalized evaporating spray-phase energy components. Curves correspond to three initial conditions: (i) solid black lines, (ii) dashed blue lines, (iii) red dots. Curves corresponding to the fluctuating energy κ, the spray-phase kinetic energy k d and the granular temperature Θ are respectively denoted through circle, square and triangle symbols.

  Fig.6Dynamics of the gas-phase energy components with evaporating droplets. Curves correspond to gas-phase TKE (top) and gas-phase internal energy (bottom) for three initial conditions: (i) solid black lines, (ii) dashed blue lines, (iii) red dots.

Fig. 7 Fig. 8

 78 Fig.7Mean velocities versus time (top) for cases (ii) (solid black lines) and (iii) (dashed red lines). Spray density concentration m 3/2 versus time for case (iii) (dashed red lines).

Fig. 9 Fig. 10

 910 Fig.9Dimensionless gas-phase TKE kg (top) and gas-phase internal energy eg (bottom) versus time. Cases (i) solid black lines, (ii) dashed red lines, and (iii) blue dots.

u d 2 d

 2 Taking the RA of (9e) yields∂t m 1 u d d + ∇x • m 1 + u ′′ d u ′′ d d + P d d = A + Mu d(106)

Table 1

 1 Parameters characterizing the dense zone of the disperse phase under ICE conditions. Values for inlet velocity u, TKE k, TKE dissipation rate ε, and the volume fraction α are given for each phase whereas rsmr denotes the Sauter mean radius (SMR) of the spray and νg is the gas-phase viscosity.

		u (m/s)	k (m 2 /s 2 )	ε (m 2 /s 3 )	α	rsmr (µm)	νg (m 2 /s)
	gas	71	111	2 × 10 6	1	-	4.15 × 10 -5
	spray	156	1	10	5 × 10 -3	25	-

Table 2

 2 Characteristic inflow values for turbulent disperse-phase flow under ICE conditions. τ S,smr and τ d,smr are the evaporation and drag times, respectively, based on the spray Sauter mean radius (SMR). τη = (νg/εg) 1/2 and Stη = τ d,smr /τη are the Kolmogorov time scale and Stokes number, respectively. St S = τ S,smr /τg is the evaporation number. All values are based on the initial conditions give in Table 1.

	τ S,smr (s)	τ d,smr (s)	τg (s)	τη (s)	Re	St	Stη	St S
	3 × 10 -3	4 × 10 -5	5.55 × 10 -5	4.56 × 10 -6	148	0.721	8.77	54.1