Song Zhao

Zakaria Bouali

Arnaud Mura
email: arnaud.mura@isae-ensma.fr

Computational investigation of weakly turbulent flame kernel growths in iso-octane droplet clouds in CVC conditions

Keywords: Turbulent combustion, Two-phase flow, Flame kernel growth, Constant volume combustion, Turbulence-scalar interaction

Numerical simulations of turbulent flame kernel growths in monodisperse clouds of iso-octane liquid droplets are conducted in conditions relevant to constant volume combustors. The simulations make use of a low-Mach number Navier-Stokes solver and a thermodynamic pressure evolution model has been implemented to reproduce the pressure variation that may be issued from either experiments or from a standard (i.e., analytical) compression law. Chemistry is described with a representative skeletal mechanism featuring 29 species and 48 elementary reaction steps. The computational results clearly confirm the enhancement of flame propagation in constant volume combustion conditions. The impact of the droplet diameter on the turbulent flame development is scrutinized for two distinct values of the Stokes number St equal to 0.1 and 1.0. Significant influence on the flame dynamics is put into evidence. This is a direct outcome of the equivalence ratio and temperature heterogeneities, which are themselves very sensitive to the choice of the Stokes number value. Then, small-scale turbulence-scalar interactions (TSI) are studied by analyzing the fields of the scalar gradients and strain-rate. Their dynamics is investigated for both non-reactive and reactive two-phase flows conditions. The TSI analysis is performed on the basis of time evolution equations written for quantities that characterize the couplings between the velocity gradient tensor and scalar gradients vectors. Special emphasis is placed on the possible influence of mass exchange terms between the liquid and gaseous phases.

Introduction

Constant-volume combustion (CVC) of fuel sprays appears as a standard combustion mode that is relevant for a wide range of applications including not only conventional piston engines but also more recently developed aeronautical engine concepts based on isochoric combustion and associated Humphrey cycle [1]. Each stage involved in such two-phase flow combustion process deserves detailed investigations: dynamics of liquid droplets, mixing between fuel and oxidizer, ignition, and turbulent flame propagation in inhomogeneous stratified media. To this purpose, conducting direct numerical simulations (DNS) of turbulent flame kernel growth in fuel sprays appears as an interesting way to analyze and understand these phenomena as well as their possible interactions. Such numerical investigations may constructively complete experimental investigations and provide additional insights onto the dynamics of the turbulent reactive flow at the smallest characteristic scales.

In this respect, the occurence of significant pressure variations is undoubtedly one of the most relevant specificities of CVC. The unburnt mixture is compressed and heated through the corresponding pressure evolution, which may influence the turbulent flame development and structure as well as auto-ignition probabilities [2]. The evaporation processes and subsequent mixing between fuel vapor and oxidizer may also be markedly altered. This is however not easy to represent all these processes in DNS of flame kernel growth. The restricted size of DNS computational domains is indeed not compatible with the resolution of the large-scale thermal and compositional stratification observed in actual combustion chambers. Thus, as regards the heat released by chemical reactions, such computations may lead to an unrealistic pressure rise.

As far as a possible volumetric compression is concerned, a classical strategy to mimic its influence consists in adding an artificial mass source term so as to regulate the pressure evolution by injecting (resp. removing) mass to (resp. from) the computational domain [2,3,4]. In the present study, this strategy is extended to low-Mach number flow solvers. In fact, it seems indeed better suited to discuss this specific issue within the low-Mach number framework where the so-called thermodynamic component p 0 (t) of the pressure field, which is associated to volumetric compression, can be easily discriminated from its dynamic counterpart p 1 (x, t), the gradient of which is involved in the Navier-Stokes equations. By considering the evolution equation of the corresponding thermodynamic component p 0 , the mechanism of pressure variation is scrutinized, leading to a computational procedure that will not alter the local flow structures and asssociated dynamics.

In turbulent two-phase flows such as those considered herein, the Stokes number (St), which is presently defined as the ratio between the particle relaxation time and the turbulence integral time scale, measures the ability of droplets to adapt to the carrier-phase velocity variations [5,6]. As previously pointed out in reference [7], the value of this dimensionless quantity does not only characterize the possible particle clustering, it also influences the topology and segregation-rate of the composition field. In the present study, some emphasis is thus placed on the impact of the Stokes number on the topology of the turbulent reactive two-phase flows. Especially, the influence of St on the combustion regime, flame propagation, and turbulent reactive flow structure is analyzed. A representative skeletal mechanism is considered to represent chemical kinetics. Since a wide range of characteristic length and time scales are involved in turbulent spray combustion, the consideration of multi-step chemistry remains quite challenging to perform such DNS of turbulent reactive two-phase flows. Among the relevant lenght scales are the integral length scale of the turbulence (l t ), the Kolmogorov length scale (η k ), the flame thickness (δ f ), and the droplet diameter (a). At this level, it is noteworthy that the description of the liquid-gas interface may also require the resolution of scales much smaller than a. As depicted in Fig. 1, these various length scales generally differ significantly from each others and simulations resolving directly all these scales with representative molecular transport and chemistry models require tremendous computational resources and therefore remain quite challenging. In the present study, the typical size L of the computational domain is chosen to be approximately equal to 10 times the integral length scale l t , while the typical size of the computational cell (i.e., the mesh resolution) is constrained according to the flame thickness δ f and Kolmogorov length scale η k . The liquid phase is represented within the discrete particle simulation (DPS) framework [7,[START_REF] Reveillon | Proceedings of the Combustion Institute[END_REF][START_REF] Kah | [END_REF]10,11]. In comparaison with computational strategies that resolve the boundary layers at the liquid-gas interfaces (e.g., [12]), such a framework allows to scrutinize the influence of large-scale dispersion processes and associated Stokes number effects while maintaining reasonable computational costs [7,10,13,14]. It is however obvious that, due to the lack of resolution at the droplet characteristic size, this computational strategy is not suited to recover certain combustion patterns, such as those associated to single droplet combustion.

In turbulent reactive two-phase flows of practical interest, small-scale scalar mixing (i.e., mixing at the molecular scales) is undoubtedly one of the key phenomena that requires to be modelled. Especially, in either non-premixed or two-phase flow conditions, the associated modelling issues are crucial since mixing between fuel and oxidizer at the molecular scales is a prerequisite before combustion occurs. Whatever the scalar is chemically-reactive (e.g., progress variable) or not (e.g., mixture fraction), its variance as well as the shape of its probability density function (PDF) are indeed settled, to a large extent, by the scalar dissipation rate (SDR), see for instance [15,16,[START_REF] Dopazo | Recent developments in PDF methods[END_REF][START_REF] Mura | [END_REF]. The SDR is the product of the scalar diffusivity with the square of the scalar gradient norm [19,20,21] and, the studies reported in references [21,[START_REF] Beishuizen | Pdf modelling and particle-turbulence interaction of turbulent spray flames[END_REF][START_REF] Xia | [END_REF]24,25,26] excepted, it is noteworthy that there were quite few efforts directed towards its detailed analysis in turbulent reactive two-phase flows. Moreover, it should be emphasized that most of the seldom studies that have been conducted in this direction were focused on the direct effect of the additional source terms associated with the evaporation process, see for instance [START_REF] Xia | [END_REF]24,25]. However, considering previous orders of magnitude analyses as well as DNS investigations of the SDR transport equation [27,[START_REF] Chakraborty | Turbulent Premixed Flames[END_REF], one may expect that the SDR is also controlled to a large extent by the turbulencescalar interaction (TSI) contribution [27,[START_REF] Chakraborty | [END_REF]. The corresponding TSI term is defined as the third-order correlation between the velocity gradient and the small-scale scalar anisotropy tensor. A classical approach to analyzing this contribution is to study the orientation of the scalar gradient in the eigenframe of the strain-rate tensor [30,31,[START_REF] Tsinober | An informal conceptual introduction to turbulence[END_REF]. In the present work, the kinematics of the orientation is studied under the influence of both evaporation and combustion processes. As far as the evaporation processes are concerned, a previous analysis has concluded that the direct influence of evaporation terms present in the mixture fraction SDR transport equation may be negligible [START_REF] Xia | [END_REF]. This does not, however, preclude the possible influence of evaporation processes on scalar orientation and associated TSI term in the SDR transport equation. The present study aims at scrutinizing this specific issue. The corresponding analysis is conducted on the basis of a direct inspection of the orientation vector transport equation. Such an equation has been recently considered to scrutinize the combustion of fully premixed reactants in reference [START_REF] Zhao | [END_REF] but, to the best of the authors' knowledge, this is the first time that such a numerical investigation is conducted within the context of turbulent reactive two-phase flows.

The present manuscript is organized as follows, in section 2, some numerical and computational details are provided, especially the strategy retained to handle the pressure evolution is presented. Section 3 is devoted to the analysis of turbulent spray dispersion and evaporation in homogeneous isotropic turbulence (HIT). The structure of the turbulent flames, as obtained from different spray and pressure evolution conditions, is discussed in section 4. It is followed by section 5, which is concerned with turbulence-scalar interactions (TSI) in reactive two-phase flows. Finally, conclusions and perspectives for future works are gathered in section 6.

Problem description

Computational model

The databases that have been generated within the course of the present study have been constructed with Asphodele [34]. This reactive flow solver, the main features of which have been already presented elsewhere -see for instance references [START_REF] Kah | [END_REF]10,[START_REF] Zhao | [END_REF]34,35,36] for further details -is based on a dilatable low-Mach number framework. The pressure field is decomposed into the sum of a thermodynamic and a dynamic contribution. The thermodynamic pressure p 0 (t) is spatially uniform. The dynamic pressure p 1 (x, t) is much smaller than p 0 (t): it is of the order of the square of the Mach number compared to the thermodynamic contribution. The former contribution p 0 is relevant to the equation of state (EoS), while the gradient of p 1 explicitly appears in the RHS of the momentum equation. Within such a low-Mach number framework, it is possible to perform the computations at a constant thermodynamic pressure by holding p 0 constant during the whole simulation time: p 0 (t) = cst, as previously done in references [START_REF] Zhao | [END_REF]36]. It is also possible to impose a given evolution for the thermodynamic pressure component, i.e., d p 0 /dt, as further discussed below. Spatial derivatives are evaluated on regular meshes by making use of high-order precision finite difference schemes while time integration is carried out with a minimal storage third-order explicit Runge-Kutta scheme. The pressure-velocity coupling is handled within a standard projection-correction framework but it is noteworthy that the numerical efficiency of the solver has been very recently improved with the implementation of the HYPRE library [START_REF] Falgout | HYPRE: a library of high performance preconditioners[END_REF]. Species mass fractions transport equations are solved for all chemical species involved in the retained chemical scheme with the molecular transport represented thanks to mixture-averaged transport properties. Finally, chemical source terms ωk are evaluated using the Sundials CVODE solver [START_REF] Hindmarsh | [END_REF].

The two phase-flows are handled within the particle-source in cell (PSI-cell) framework [39]. The carrier gas phase is described from an Eulerian viewpoint based on the standard set of Navier-Stokes equations (NSE) within the low Mach number framework described above. Conservation equations are considered for the mass (i.e., density of the mixture ρ), momentum components (ρu i ), mass fraction of each chemical species (Y α ), and sensible enthalpy

h s = α=Nsp α=1 Y α h s,α = α=Nsp α=1 Y α T T0 C p,α (θ)dθ (1)
with h s,α and C p,α the sensible enthalpy and constant pressure heat capacity of the αth species, respectively, and T 0 = 298.0 K a reference temperature. Thus, in the gaseous phase, the corresponding set of conservation equations may be written as follows

∂ρ ∂t + ∂ ∂x j ρu j = ḋρ ( 2a
)
∂ρu i ∂t + ∂ ∂x j ρu i u j = - ∂p 1 ∂x i + ∂τ i j ∂x j + ḋρui , for i = 1, . . . , 3 (2b)
∂ρh s ∂t + ∂ ∂x j ρh s u j = d p 0 dt + ωT - ∂q j ∂x j + τ i j ∂u i ∂x j + ḋρhs (2c) ∂ρY α ∂t + ∂ ∂x j ρY α u j = - ∂J α j ∂x j + ωα + ḋρYα for α = 1, . . . , N sp (2d)
where the substantial derivative of the pressure is involved in the RHS of the sensible enthalpy budget (i.e., Eq. 2c), while the spatial derivative of its dynamic component p 1 is present in the momentum transport equation (Eq. 2b).

In the above equation, the expression of the viscous stress tensor is

τ i j = µ ∂u i ∂x j + ∂u j ∂x i - 2 3
∂u k ∂x k δ i j , (3)
where µ denotes the viscosity of the mixture. The quantity J α j = ρY α V α j in the species mass fraction transport equation (i.e., Eq. 2d) denotes the component in direction j of the mass flux of chemical species α, the vector α being the molecular diffusion velocity. The components of this diffusion velocity are determined from a Fickian approximation

ρY α V α j = -ρD α ∂Y α ∂x j + ρY α β=Nsp β=1 D β ∂Y β ∂x j (4)
with D α the molecular diffusion coefficient of species α. The second term in the RHS of Eq. 4 is a correction that ensures the compatibility between the species mass fraction transport equation and the continuity equation (i.e., total mass conservation). The heat flux term q j in Eq. 2c is determined from

q j = -λ ∂T ∂x j + α=Nsp α=1 J α j ∂h s ∂x j (5)
with λ the conductivity of the reactive mixture. Finally, the above system is complemented with the perfect gas equation of state (EoS)

p 0 = ρRT M , 1 M = α=Nsp α=1 Y α M α (6)
where standard notations have been retained: R is the universal gas constant, M is the molecular weight of the multi-component mixture, and M α , the one of the chemical species α.

There are also some source terms present in the above set of conservation equations. For instance, chemical sources terms ωα and ωT do appear. The former corresponds to the production/destruction rate of chemical species α and the latter denotes the resulting heat release rate (HRR). In the present study, the combustion of a primary reference fuel (PRF), iso-octane (C 8 H 18 ), is considered. Its chemistry is described with a skeletal mechanism featuring 29 chemical species and 48 elementary reactions [40].

The other terms ḋρ , ḋρui , ḋρhs , and ḋρYα represent the mass, momentum, energy, and species mass fraction sources that are associated to droplet effects (momentum couplings and evaporation). The last contribution (i.e., ḋρYα ) is non-zero only if the chemical species α corresponds to the chemical species present in the liquid (i.e., C 8 H 18 ). These various source terms are evaluated from the contributions of each liquid droplet through a volumetric weighted average performed on each computational cell in a way similar to the one previously described in reference [7].

Liquid droplets are modeled as spherical particles within a Lagrangian framework. Their position (X d ), velocity (V d ), temperature (T d ), and diameter (a) are computed from the following set of evolution equations:

dX d dt = V d (7a) dV d dt = 1 τ Vd (u(X d , t) -V d ) ( 7b
)
dT d dt = 1 τ Td T (X d , t) -T d - B T L v C p (X d , t) (7c)
da 2 dt = - a 2 τ a (7d)
where u, T , and C p denote the velocity, temperature, and constant pressure heat capacity of the gas, which can be determined from the characteristics of the gaseous mixture resolution according to the set of Navier-Stokes transport equations (i.e., Eq. 2). The quantities B T and L v denote the thermal transfer number and the latent heat of evaporation, τ Vd and τ Td are the instantaneous kinetic and temperature relaxation time of each droplet, respectively, while τ a denotes the characteristic time of evaporation. In practice, the kinetic relaxation time τ Vd may be evaluated from

τ Vd = ρ d a 2 18C u µ(X d , t) (8)
where ρ d is the liquid phase density, µ(X d , t) the carrier phase viscosity, and C u = 1 + Re d 2/3 /6 is a corrective coefficient that is introduced to vary the drag factor according to local droplet Reynolds number value Re d = ρa|u(X d , t) -V d |/µ. The heating and evaporation of droplets are modeled by Eq. 7c and Eq. 7d with the temperature and evaporation relaxation times computed as

τ Td = Pr 6Nu c ρ d a 2 µ R Cp B m ln (1 + B m ) (9a) τ a = Sc 4Sh c ρ d a 2 µ 1 ln (1 + B m ) (9b)
where Sc and Pr are the Schmidt and Prandtl numbers, Sh c and Nu c denote the convective Sherwood and Nusselt numbers respectively, R Cp denotes the ratio of normalized heat capacities of liquid and gaseous phase. The quantity B m represents the normalized mass flux of liquid through the droplet surface, which is evaluated from

B m = (Y s F -Y F (X d , t))/(1 -Y s F ), with Y F (X d , t
) the fuel mass fraction in the vicinity of the droplets and Y s F its value at the liquid-gas interface, which is deduced from the Clausius-Clapeyron saturation law.

Because of their small variations with respect to temperature, density (ρ l ) and heat capacity of the liquid (C p,l ) are considered as constant values. As regards the latent heat of vaporization, it is estimated from the following expression:

L v,l = L v,b
T cr -T l T cr -T b 0.38 (10) where L v,b is a function of the temperature (T cr ) and pressure (P cr ) at the critical point, and boiling temperature (T b ). The values of these properties are gathered in Table 1. The thermodynamic pressure variation models applied in the current study will be presented below. In practice, the corresponding pressure variations are processed through a supplementary source term ṁp , which is added to the right hand side (RHS) of the mass conservation equation (i.e., Eq. 2a). The corresponding quantity is an additional space-independent artificial mass source term, which aims at regulating the thermodynamic pressure variation. As emphasized above, within the low-Mach number framework, the pressure field p(x, t) can indeed be split into a space-independent thermodynamic pressure field p 0 (t) plus a dynamic pressure field p 1 (x, t):

p(x, t) = p 0 (t) + γ ∞ Ma 2 ∞ p 1 (x, t) + O(Ma ∞ 4 ) (11)
where γ ∞ and Ma ∞ are the reference heat capacity ratio and Mach number. It is worth recalling that the temporal variation of the pressure, which does appear in the RHS of the enthalpy conservation Eq. 2c, is zero in open domains. From a more general viewpoint, see Appendix A, it can be explicitly derived and evaluated in any volume V from the multicomponent mixture EoS (Eq. 6):

d p 0 dt = 1 V V                       (γ -1) ρ d s dt (I) + α=Nsp α=1 γr α T -(γ -1) h s,α ρ dY α dt (II) + γrT ṁp + ḋρ - ∂ρu j ∂x j (III)                       dV , (12)
In the above expression, γ is the local heat capacity ratio, r = R/M (note that r α = R/M α ) and the quantity ρd s /dt represents the sum of chemical heat release, heat transfer, and energy transfer induced by evaporation, i.e., ρd s /dt ≡ ωT -∂q j /∂x j + ḋρhs . Unsurprisingly, three possible thermodynamic pressure variation mechanisms are found in the RHS of Eq. 12. The first (I) is related to pressure variation induced by the heat release, heat transfer, and evaporation terms. The second (II) is associated with the variation of thermodynamic properties of the mixture as a consequence of composition changes. Finally, the third and last contribution (III) reflects the effects of density variation, associated to mass exchange between the two phases, volumetric expansion, etc.

Based on this equation, different thermodynamic pressure evolution models can be proposed. On the one hand, a constant pressure (CP) model may be recovered by setting

d p 0 dt CP ≡ 0 . ( 13
)
and the corresponding artificial mass source can be deduced from Eq. 12,

ṁp CP = - 1 (γrT )V V         (γ -1) ρ d s dt + α=Nsp α=1 γr α T -(γ -1) h s,α ρ dY α dt + γrT ḋρ - ∂ρu j ∂x j         dV (14)
and applied in the mass conservation Eq. 2a. It is worth noting that, for a combustion process, this artificial mass source term should be negative so as to balance the pressure increment that is induced by the heat released by chemical reactions. On the other hand, a constant volume (CV) evolution model can be achieved by setting

ṁp CV ≡ 0 . (15)
In this respect, it is worth noting that, although the artificial mass injection is set to zero, the mass sources issued from liquid evaporation ( ḋρ ) are still considered for two-phase flow applications. Thus, the thermodynamic pressure varies according to

d p 0 dt CV = 1 V V         (γ -1) ρ d s dt + α=Nsp α=1 γr α T -(γ -1) h s,α ρ dY α dt + γrT ḋρ - ∂ρu j ∂x j         dV . (16)
Finally, an imposed pressure (PT) model may be designed to regulate the thermodynamic pressure in the domain according to a prescribed pressure history. Under such conditions, the pressure variation may be denoted d p 0 dt PT [START_REF] Dopazo | Recent developments in PDF methods[END_REF] and the artificial source term is

ṁp (t) PT = 1 γrT         d p 0 dt PT - 1 V V         (γ -1) ρ d s dt + α=Nsp α=1 γr α T -(γ -1) h s,α ρ dY α dt + γrT ḋρ - ∂ρu j ∂x j         dV         (18)
Thanks to this model, simulations following a given pressure variation, e.g., an experimental pressure trace, can be achieved.

Computational setups

The database gathers several configurations, the main characteristics of which have been summarized in Table 2. Simulations 3D_NR_CP_ST10 and 3D_NR_CP_ST01 correspond to three-dimensional non-reactive cases associated to two distinct initial values of the Stokes number St = τ Vd /τ T (i.e., two values of the initial droplet diameter1 ). It is worth noting that the present definition of the Stokes number differs from the one considered by Wang and Maxey [5], which is based on the Kolmogorov time scale. Simulations 3D_NR_CP_ST10 and 3D_NR_CP_ST01 are used to analyze the effects of droplet dispersion and evaporation. Computations labelled 2D_R_CP_ST10, 2D_R_CV_ST10, 2D_R_CP_ST01, and 2D_R_CV_ST01 correspond to two-dimensional simulations of two-phase flow flame kernel growths in HIT performed in either constant volume or constant pressure conditions. They are performed to evaluate the impact of both the Stokes number St and pressure model on the flame. Finally, the database also includes three-dimensional constant pressure computations of flame kernel growths in HIT (3D_R_CP_ST10 and 3D_R_CP_ST01). The main objectives with these three-dimensional reactive cases are (i) to assess the possible impact of three-dimensional effects and (ii) to proceed with a refined analysis of turbulence-scalar interactions (TSI). The computational domain dimension is the same in the three directions, i.e.,

L x = L x1 = L x2 = L x3 .
It is approximately 10 times the integral length scale of the turbulence and it is discretized with N x = N x1 = N x2 = N x3 grid points along each direction, in such a manner the computational mesh resolution ∆x = L x /N x is kept constant and equal to ∆x = 15.625 µm for all simulations. The initial fluctuating velocity field of the gas carrier phase is generated on the basis of the classical Rogallo's procedure [START_REF] Rogallo | Numerical experiments in homogeneous turbulence[END_REF] according to a Passot-Pouquet spectrum [42]. A constant-energetics physical-space forced homogeneous isotropic turbulence (HIT) is maintained by resorting to the procedure described in reference [43]. The present computations have been conducted with a Courant and Fourier number values set to 0.1. The integration time step ∆t is restricted by molecular diffusion and it is of the order of 10 -4 ms. This leeds to high computational costs of around 490 CPU hours for a two-dimensional and 60,000 CPU hours for a three-dimensional reactive simulation performed on Intel Sandy Bridge E5-2680 cores. The liquid fuel is initialized as a mono-disperse spray of iso-octane droplets that is homogeneously distributed in the gas. The initial droplet temperature is uniform and set equal to T l =300 K. The total number of droplets in the computational domain is determined in such a manner that the global, i.e., overall, fuel/air equivalence ratio is around 0.9. In this respect, it is noteworthy that, consistently with the use of the DPS framework, the average liquid volume fraction remains rather small: its maximum level corresponds to its initial value, which does not exceed 0.01%. The iso-octane vapor concentration in the gas phase at the initial time is Y C8H18 (t = 0) = 0.026, which corresponds to an equivalence ratio Φ = 0.4. At this level, it should be emphasized that, from a practical viewpoint, it is possible to proceed with the ignition of a mixture featuring a zero concentration of iso-octane in the gas phase and, provided that a sufficient deposit of energy is considered, such an ignition procedure may indeed be successful. However, it has been experienced that the required amount of energy is so large that the homogeneity of the droplet distribution in the computational domain and the subsequent propagation of the flame were significantly altered. Therefore, the initial fuel vapor concentration has been set to 0.026, a value that corresponds to the lean limit of C 8 H 16 -air flame propagation (i.e., Φ = 0.4). In these conditions, successful ignition can be reached even for a limited amount of energy deposit. The computations are conducted according to the following procedure: in the first part of the simulation, only a one-way momentum coupling is considered between the two phases (i.e., no heat and mass transfer): the liquid particles disperse within a constant-energetics physical-space forced HIT. After about 10 integral time scales, the averaged turbulent flowfield reaches its steady-state in terms of velocity statistics, droplet clustering, etc. The main features of the corresponding flowfield are gathered in Table 3. In this respect, it is important to emphasize that, for all computations of the present database, the characteristics of the turbulent flowfield are the same. The computational grid size is about one half of the Kolmogorov length scale at this time, thus ensuring a satisfactory level of resolution of the velocity field. The two-way coupling of mass, heat, and momentum are then activated and the simulation is pursued without any forcing so as to simulate dispersion and evaporation of the droplet cloud in a self-decaying HIT. Moreover, chemical reactions are activated in the two-dimensional computations (2D_R_CP_ST10, 2D_R_CV_ST10, 2D_R_CP_ST01, and 2D_R_CV_ST01) and in the three-dimensional simulations 3D_R_CP_ST10 and 3D_R_CP_ST01, see Table 2. These different cases indeed correspond to reactive conditions for which a small flame kernel is initialized in the center of the computational domain. The generation of this flame kernel follows the same procedure as the one previously described in reference [START_REF] Zhao | [END_REF].

Spray dispersion and evaporation

Based on simulations cases 3D_NR_CP_ST10 and 3D_NR_CP_ST01, a preliminary study of the dispersion and evaporation processes can be conducted. As emphasized in the previous section, the only difference between these two cases is related to the characteristic droplet size. In configuration 3D_NR_CP_ST10, the initial droplet diameter is 6.1 µm while in case 3D_NR_CP_ST01 it is 1.9 µm. As shown in Fig. 2, after being dispersed by the same constant-energetics forced turbulence, very different droplet cloud topologies are obtained. This is a direct consequence of the difference in the particle relaxation time. Droplet clouds featuring unity Stokes number produce smaller clusters, i.e., higher local particle densities and larger droplet segregation levels, compared to the other cases. First row: evolutions plotted versus time t (physical units). Second row: evolutions plotted versus t/τ a with τ a the evaporation time scale evaluated from the initial characteristics of the two-phase flow.

Y F, RMS /Y F (t = 0) (b) 0.0 0.5 1.0 1.5 2.0 2.5 t/τ a 1.0 1.2 1.4 1.6 1.8 2.0 Y F /Y F (t = 0) (c) 0.0 0.5 1.0 1.5 2.0 2.5 t/τ a 0.0 0.2 0.4 0.6 0.8 1.0 Y F, RMS /Y F (t = 0) (d)
After droplet clustering, the evaporation process is started. Evaporation and segregation processes impact the evolutions of the mean value of fuel vapor mass fraction and the root-mean-square (RMS) of associated fluctuations. Indeed, the averaged fuel concentration increases much faster in case 3D_NR_CP_ST01 compared to case 3D_NR_CP_ST10 (see Fig. 3(a)), which confirms that much higher evaporation rate levels are reached within the spray that features the smaller characteristic droplet diameter (3D_NR_CP_ST01). Moreover, a larger RMS of the fuel vapor mass fraction is obtained from simulation 3D_NR_CP_ST10 as shown in Fig. 3(b). This high segregation level also lasts over longer time in case 3D_NR_CP_ST10.

In this respect, it seems worth emphasizing again that the single-droplet characteristic evaporation time is approximately ten times larger in case 3D_NR_CP_ST10 compared to case 3D_NR_CP_ST01, which unambiguously explains why evaporation is much faster in case 3D_NR_CP_ST01. However, it is noteworthy that, even once normalized by the single-droplet evaporation characteristic time τ a , which is about 0.2 ms for 3D_NR_CP_ST10 and 0.02 ms for 3D_NR_CP_ST01, these evolutions show that the droplet spray featuring a unity St number (3D_NR_CP_ST10) still has a longer evaporation time and leads to a higher fuel vapor segregation level. This is a direct outcome of droplet clustering, which leads to higher local fuel concentrations in the gas, thus slowing down the evaporation process in case 3D_NR_CP_ST10 compared to case 3D_NR_CP_ST01. This set of observations agrees well with conclusions previously drawn in the literature [7,[START_REF] Reveillon | Proceedings of the Combustion Institute[END_REF].

In addition to the mean value of fuel vapor concentration and RMS of associated fluctuations, the integral characteristic length scale of the composition field may also significantly impact the subsequent turbulent combustion process. Normalized energy spectra of the fluctuations of the composition field are reported in Fig. 4. For the purpose of the present analysis, the composition field is followed with a scalar field ζ, which is defined as a normalized fuel mass fraction, i.e.,

ζ ≡ Y C8H18 -Y min C8H18 / Y max C8H18 -Y min C8H18 .
The associated characteristic length scale l ζ can be deduced from the most energetic wave number (i.e., K max ) of the fluctuating scalar energy spectrum: l ζ ≡ 2π/K max . In fact, the determination of this characteristic length scale may raise some difficulties since two local maximum can be identified during the early stages of the evaporation processes, see Fig. 4. The largest one may be associated to the fluctuating velocity field (i.e., the integral length scale of turbulence), while the second one corresponds to a characteristic dispersion length scale of the droplet cloud. Unsurprisingly, as evaporation proceeds, the impact of droplet dispersion is filtered out by the turbulence, and scalar spectra similar to velocity spectra are finally achieved, as shown in Fig. 4(d).

The behaviour of the characteristic scalar length scale l ζ (not reported for the sake of conciseness) is such that, except at the very early stages of the simulation, it remains larger in the droplet spray featuring the smallest characteristic size (i.e., case 3D_NR_CP_ST01). However, this is in contrast with the behaviour observed when the temporal range is restricted (for both conditions 3D_NR_CP_ST10 and 3D_NR_CP_ST01) to one single characteristic evaporation time. In this case, the condition 3D_NR_CP_ST10, which is associated to larger droplet size, leads to a larger length scale value l ζ . This behaviour can be explained as follows: at the beginning of the evaporation process, l ζ is mainly controlled by the initial droplet distribution. For case 3D_NR_CP_ST10, as the droplets are much clustered, see Fig. 2(a), the value of l ζ is mainly settled by the size of the corresponding clusters. This is in contrast with case 3D_NR_CP_ST01 where droplets are almost homogeneously distributed. In this case, l ζ is mainly controlled by the droplet characteristic size itself (i.e., the droplet diameter). This conclusion is supported by the presence of a peak in the ζ-energy spectrum, which lies at a value that is almost the Nyquist cut-off length scale K/K Lx = 128 = N x /2 and corresponds to the mesh resolution.

Before embarking on the description of the turbulent reactive two-phase flows cases, the possible turbulent combustion regimes are discussed on the basis of a preliminary inspection of the present non-reactive evaporating two-phase flows (i.e., cases 3D_NR_CP_ST01 and case 3D_NR_CP_ST10). This is a worthwhile discussion since, as emphasized in the introduction section, single droplet combustion regime should be avoided as much as possible because they cannot be handled properly within the standard DPS numerical framework. To this purpose, the spray separation number S p is considered so as to characterize the possible two-phase flow flame patterns. This dimensionless number may be defined as the ratio between the mean droplet spacing (δ p ) and the diffusion flame radius surrounding a single droplet in a quiescent oxidizer (δ r f ). Provided that the dilution level is sufficiently high, the mean droplet spacing can be approximated from δ p = n -1/3 , with n the droplet density (i.e., the number of droplets per unit volume). In this respect, it should be emphasized that the product n 1/3 • δ r f (which is 1/S p ) has been early considered in reference [44]. It is sometimes referred to as the Kerstein-Law parameter [START_REF] Borghi | Symposium (International) on Combustion[END_REF]. Based on the percolation theory, there exist two critical values of 1/S p that identify three possible combustion regimes [START_REF] Borghi | Modélisation et théorie des flammes[END_REF]. On the one hand, as long as 1/S p does not exceed a first critical value approximately equal to 0.43, the spray remains rather dispersed and droplets are expected to burn either individually or within small groups. Such conditions correspond to isolated droplet or group combustion regimes. On the other hand, when 1/S p exceeds another critical value (approximately equal to 0.77), group droplet combustion (and individual droplet burning as its limiting case) is no longer possible, and the corresponding combustion mode is termed pocket combustion regime. Finally, the combustion regimes between these two extreme conditions are referred to as percolation combustion regimes. Considering the large variation of local droplet density that results from the clustering induced by the turbulent flowfield, the value 1/S p = δ r f /δ p is presently evaluated for each droplet. The local droplet spacing δ p is calculated as the distance to the nearest neighboring droplet, and the individual flame radius δ r f is determined from its standard expression [START_REF] Borghi | Modélisation et théorie des flammes[END_REF] with the actual fuel mass fraction on the droplet surface. The PDFs of 1/S p obtained for both cases 3D_NR_CP_ST10 and 3D_NR_CP_ST01 at the beginning of the evaporation phase are reported in Fig. 5. From these results, it can be concluded that either the group combustion or the isolated droplet combustion regimes are quite unlikely to take place for the present set of conditions, which brings further support to the use of the retained numerical strategy. The analysis of the corresponding turbulent two-phase flow flames is the subject of the next section.

Turbulent two-phase flow flames

In this section, turbulent flame developments in sprays characterized by different Stokes number values are studied for two pressure variation models. Instantaneous snapshots of turbulent velocity, equivalence ratio, and flame front issued from simulations 2D_R_CP_ST10, 2D_R_CV_ST10, 2D_R_CP_ST01, and 2D_R_CV_ST01 are presented in Fig. 6. This figure provides an overview of the flame kernel growth development in the droplet clouds. At the initial time, see the first row of Fig. 6, droplets have already been clustered by the forced HIT, so that the sprays featuring a unity Stokes number form smaller clusters, which are located outside the instantaneous vortical structures. At time t = τ T , as shown in the second row of Fig. 6, for unity Stokes number value simulations (2D_R_CP_ST10, 2D_R_CV_ST10), droplets start to interact with the flame front, thus altering its propagation. Some local extinctions do appear because of large equivalence ratio variations. This is in contrast with the simulations characterized by a smaller Stokes number value (2D_R_CP_ST01, 2D_R_CV_ST01), in these cases, sprays have already been fully evaporated and flames propagate in nearly homogeneous gaseous mixtures.

Although there are still some flame droplet direct interactions at this time in simulations featuring unity Stokes number (2D_R_CP_ST10, 2D_R_CV_ST10), the resulting turbulent flames are mostly premixed flames. This is confirmed from the inspection of the PDFs of normalized premixedness index2 reported in Fig. 7.

From Fig. 7, it is clear that, even in simulations featuring the biggest droplets, the amount of diffusion combustion remains relatively small at this time. In practice, the computation of the above index is conditioned to a non-negligible value of the heat release rate, which is presently set to one percent of the corresponding maximum. At time t = 3τ T and t = 4τ T , see the third and fourth rows of Fig. 6, all flame kernels go on expanding in the fresh gas mixture. However, smaller propagation velocities and larger flame curvature levels are observed for situations featuring a unity Stokes number (St=1.0) compared to cases with St=0.1. Finally, the flame kernel development is faster with the CV model in comparison with the growth obtained with the CP model. The influence of (i) spray characteristic Stokes number and (ii) thermodynamic pressure variation on the flame propagation are quantified by considering the temporal evolution of its characteristic radius as evaluated from the area occupied by burnt gases. This area allows to estimate the flame radius in laminar cases or an equivalent flame radius corresponding to a spherical flame (i.e., circular in two dimensions) featuring the same amount of burnt gases in turbulent cases. This evolution is plotted together with the PDF of the flame front local displacement speed in Fig. 8. The time evolution of this equivalent radius reflects the average consumption flame speed S c . The burnt gas condition has been determined from c > 0.9 with c the progress variable. In the present study, this variable is evaluated from c = Y c /Y eq c , with Y c the sum of the mass fractions of carbon monoxide and dioxide, i.e., Y c = Y CO2 + Y CO and Y eq c the corresponding value at equilibrium. It seems worth emphasizing that the flame front displacement speed reported in Fig. 8(b) has been evaluated on the isoline c = 0.05 from the following expression

S d = 1 |∇c| dc dt = 1 |∇c|           1 Y eq c dY c dt + Y c Y eq c 2 dY eq c dt           (19)
Considering the dependency of the composition at chemical equilibrium on the pressure p, temperature T , and chemical species mass fraction Y k , the substantial derivative of Y eq c can be expanded as follows The internal structure of local flame front and the way it may be altered by the droplet cloud and turbulence in either CP or CV conditions is now inspected. In this respect, it seems worth recalling that the computational mesh has been designed in such a manner that ∆x is of the order of 10 µm, which allows for a satisfactory level of resolution of the flame thickness. Scatterplots of the magnitude of the corresponding progress variable gradient |∇c| are plotted versus the progress variable c in Fig. 9 together with the corresponding conditional averages. These scatterplots show that the turbulent flames associated to sprays featuring unity Stokes number are much more dispersed, thus revealing that their inner structures is significantly altered compared to those associated to smaller droplets (St=0.1). In addition to this, consistently with the results obtained for the flame displacement speed, CV simulations are found to lead to thinner flames, i.e., larger progress variable gradients inside the flame.

The total derivatives of the primitive variables that remain in the expression of S d are directly saved from the computations, while the partial derivatives of Y eq c with respect to p, T , and Y k can be deduced from chemical equilibrium computations. Thus, the flame displacement speed is fully determined. The results reported in Fig. 8(a) undoubtedly confirms that, in situations featuring a smaller droplet size (or a constant volume model), the flame kernels development is faster than in cases associated to larger droplets (or constant pressure model). The same conclusion can be drawn from the inspection of the local flame displacement speed PDFs gathered in Fig. 8(b).

The influence of the Stokes number can be readily explained: sprays with St=1 are associated with higher segregation levels. Moreover, as shown in section 3, they also display smaller characteristic length scale of fluctuations l ξ compared to computations performed with St=0.1. Thus, during the early stages of the flame kernel development, the corresponding variations of local equivalence ratio can lead to flame extinction, as it has been shown above in the second row of Fig. 6.

In CV conditions, the enhancement of flame propagation results from the increase of the unburnt gases temperature, which is a direct consequence of the compression induced by the production of (light) burnt gas within the constant volume. The time evolution of the thermodynamic pressure component and fresh reactants temperature has been reported in Fig. 10. It can be seen that, as an outcome of the heat losses induced by liquid evaporation, both CP and CV computations first exhibit a pressure (and a temperature) drop. In this respect, it is noteworthy that the simulations conducted with smaller droplets (2D_R_CP_ST01 and 2D_R_CV_ST01) display a faster decrease. Then, as the flame kernels develop, the heat released by chemical reactions leads to an increase of p 0 and to an associated increase of enthalpy (i.e., temperature) of the fresh reactants. This is in contrast with the CP model, for which p 0 is set constant: heat release by chemical reactions cannot affect the fresh reactant temperature, thus leading to a smaller flame propagation velocity compared to CV conditions. Configurations 3D_R_CP_ST10 and 3D_R_CP_ST01 correspond to three-dimensional computations of flame kernel growth in droplet clouds featuring St=0.1 and St=1.0. Figure 11 displays a snapshot of the turbulent reactive two-phase flows obtained at time t = τ T . At this time, the evaporation processes are finalized in case 3D_R_CP_ST01 and they are in their final stage in case 3D_R_CP_ST10. The corresponding premixedness index, local flame displacement speed, and flame internal structure statistics are reported in Fig. 12. In comparison with two-dimensional computations discussed above, there is now almost no diffusive combustion as indicated by the PDFs of premixedness index in Fig. 12(a). The PDFs of the flame displacement speed reported in Fig. 12(b) is quite similar to those obtained in the corresponding two-dimensional cases (Fig. 8). Finally, the flame inner structures issued from the three-dimensional simulations seems to be less influenced by the Stokes number compared to their two-dimensional couterparts, see Fig. 9(a) and 9(c). Considering that more sampling points are available from the three-dimensional computations, three-dimensional turbulence leads to enhanced flame curvature and scalar mixing, the observed differences remain quite reasonable. In the next section, the corresponding set of three-dimensional data will be used to analyze in more details the processes of turbulence-scalar interactions in reactive two-phase flows.

Small-scale turbulence-scalar interactions

The database is now used to scrutinize the turbulence-scalar interaction (TSI) in reactive two-phase flows. The TSI term denotes one of the leading-order source terms that appears in the scalar dissipation rate (SDR) transport equation. It is defined as the third-order correlation between the velocity gradient and the small-scale scalar anisotropy tensor. It should be emphasized that, with only of few exceptions, e.g. [12,[START_REF] Xia | [END_REF][START_REF] Ahmed | Proceedings of the Combustion Institute[END_REF], the impact of either the evaporation or the combined action of evaporation and chemical reactions on scalar mixing and associated TSI term has been quite seldom studied in the literature. This is in sharp contrast with the amount of turbulencescalar interaction studies conducted with fully premixed reactants. The most significant conclusions drawn from these studies have been gathered in a book chapter by Chakraborty et al. [START_REF] Chakraborty | Turbulent Premixed Flames[END_REF]. They are also summarized in review papers by Lipatnikov [START_REF] Lipatnikov | [END_REF] and by Sabelnikov and Lipatnikov [49]. Considering the difficulties associated to the simultaneous characterization of the velocity and scalar fields at the smallest scales of turbulent reactive flows, most of these studies were based on DNS data (e.g., [50]) but some of them also rely on experimental investigations. For instance, turbulence-scalar interactions have been studied in premixed flames stabilized on a bluff body by Hartung et al. [51]. It has been followed by the investigation of turbulent premixed Bunsen burner flames conducted by Steinberg et al. [52]. In the latter, by following individual turbulent structures in a Lagrangian manner, it has been shown that the flame affects the scalar alignment through the rotation of the structures as they approach the flame front 3 . The influence of combustion on principal strain-rate magnitudes has been also investigated in reference [55]. In this respect, the conclusions that were drawn from these studies [52,55] have been recently confirmed through the inspection of data [START_REF] Zhao | [END_REF].

Let us consider for instance a scalar quantity ϕ that is bounded between zero and unity, e.g., a mixture fraction or progress variable. Denoting by g ϕ = ∇ϕ its local gradient and D its molecular diffusivity, the TSI term writes

(TSI) ≡ -2ρD g T ϕ • S • g ϕ = 2ρD| g ϕ | 2 -n T ϕ • S • n ϕ = 2ρD|g ϕ | 2         - i=3 i=1 λ i n 2 ϕ,i         (21)
where n ϕ is the unit vector associated to the scalar gradient direction, i.e., n ϕ ≡ g ϕ /|g ϕ |. The quantity λ i denotes the i th eigenvalue of the strain-rate tensor S = ∇u + ∇u T /2 and n ϕ is the projection of n ϕ on the eigenframe of S, which is evaluated from

n ϕ,i ≡ n ϕ • e i = |n ϕ | |e i | cos(θ i ) = cos(θ i ) ( 22
)
with e i the i th unit eigenvector of S and θ i the angle between the scalar gradient g ϕ and eigenvector e i . It is noteworthy that, by convention, the principal strain-rates values (or eigenvalues) are sorted in descending order such that λ 1 > λ 2 > λ 3 . The corresponding eigenvectors are associated to the most extensive (e 1 ), intermediate (e 2 ), and compressive (e 3 ) principal directions. Clearly, the sign of the TSI is given by the normalized TSI term, i.e., -n T ϕ • S • n ϕ , and it is thus settled by the orientation of scalar gradient in the eigenframe, n ϕ and eigenvalues λ i . According to Eq. 21, its magnitude does also involve the SDR level

N ϕ = D|g ϕ | 2 .
The geometrical analysis of the scalar orientation has concentrated a large amount of work over the last recent years. In standard (i.e., incompressible and homogeneous) turbulence, the scalar gradient is known to be mostly aligned with the compression principal direction e 3 , see for instance [31,56,57], and this leads the TSI to be positive, it increases the SDR. However, as far as turbulent reactive two-phase flows are concerned, the question of the possible influence of evaporation combined with combustion remains mostly unanswered. Therefore, the evolution of the scalar gradient orientation is presently scrutinized during the evaporation stage and some attention is also subsequently paid to the possible influence of combustion-induced thermal expansion.

In a first step of the analysis, attention is focused on the sole influence of the evaporation process. As discussed in the previous section, the distributions of equivalence ratio and temperature in the fresh reactants drive to a large extent the evolution of the turbulent reative flow. The following analysis is therefore aimed at analyzing the turbulent scalar mixing during the evaporation process on the basis of the three-dimensional cases 3D_NR_CP_ST10 and 3D_NR_CP_ST01. The scalar under consideration is a normalized value of the ratio

y CO = Y C /Y O with Y C
and Y O the mass fractions of C and O atoms, respectively 4 . It is defined as follows ξ = y CO /y CO,sty CO /y CO,st min y CO /y CO,st maxy CO /y CO,st min [START_REF] Xia | [END_REF] where y CO,st = 8M C /(25M O ). It is should be emphasized that ξ is not a passive scalar in the presence of the evaporation process since there is a source of carbon atoms in the gaseous phase. PDFs of the normalized TSI term, scalar gradient orientation, and eigenvalues of the strain-rate tensor are reported for different simulation times in Fig. 13. At the very beginning of the evaporation process (see Fig. 13(a)), there is no striking difference in the normalized TSI terms issued from both configurations featuring different Stokes numbers. The obtained results do not exhibit any special or preferential alignment of the scalar gradient in the eigenframe. The PDFs of strain-rate eigenvalues issued from both simulations are also almost identical. Since both databases are based on the same initial turbulent flowfield, these results are quite reasonable. However, at t = 0.5τ T , although the eigenvalues remains almost identical, as it can be seen in the third column of Fig. 13(b), the normalized TSI associated to case 3D_NR_CP_ST01 becomes larger than the one of case 3D_NR_CP_ST10. The preferential alignment between the scalar gradient and the compression direction e 3 is more significant in case 3D_NR_CP_ST01 than in case 3D_NR_CP_ST10. Thus, the normalized TSI, which is the sum of the square of n i weighted by the eigenvalue λ i , shows a higher probability to be positive in case 3D_NR_CP_ST01. Finally, at time t = τ T (Fig. 13(c)), the differences in the PDFs of scalar gradient alignments are significantly lessened (see the second column of Fig. 13(c)), the normalized TSI terms issued from both simulation cases tend to recover some similarity. From the evolution reported in Fig. 13, one may ascribe the observed differences to the evolution of the alignment that takes place the during the evaporation process. A similar conclusion can be drawn from Fig. 14, which displays the temporal evolution of the volume-averaged values of the normalized TSI and alignment components. It is clear that the averaged TSI increases faster in case 3D_NR_CP_ST01, which is consistent with the evolution observed for the alignment vector n. The evolution discussed above motivates a detailed inspection of the kinematics of the scalar gradient alignment in the strain-rate eigenframe. Therefore, a transport equation of n is considered. The details of its derivation are provided in Appendix C, and the resulting form can be written as

∂ n ∂t = -TR + TS + TW + GR + WN (24)
In this equation, TR ≡ u•∇ n denotes the transport term, TS ≡ -R T Inn T Sn = n T Λ n I -Λ n corresponds to the effects of the strain-rate tensor, TW ≡ -R T (n × ω) /2 refers to the contribution associated to vorticity, and

GR ≡ |g| -1 I -n n T R T ∇ (Dξ/Dt) ( 25
)
is the contribution relevant to the evolution of the scalar field itself, i.e., molecular diffusion and productiondestruction term induced by evaporation (or combustion processes if a reactive scalar is considered). Finally, the rotation of the eigenframe itself contributes to the variation of the alignments, and the corresponding term is denoted WN. In the above expressions, I is the identity matrix, Λ denotes the diagonal matrix made up of the three eigenvalues of S, i.e Λ ≡ diag(λ 1 , λ 2 , λ 3 ), and R is the eigenmatrix of S with R ≡ [e 1 |e 2 |e 3 ]. Fig. 15: Temporal evolution of the transport equations for the three components of the orientation vector i.e., the projection of the scalar gradient in the eigenframe of the strain-rate tensor, as given by Eq. 24. Solid line: case 3D_NR_CP_ST10, dashed line: case 3D_NR_CP_ST01.

Equation 24 has been scrutinized in its averaged form for both 3D_NR_CP_ST10 and 3D_NR_CP_ST01 cases. The corresponding budgets are displayed in Fig. 15. The substantial derivatives of the three components of the orientation vector are reported in Fig. 15(a). The corresponding evolution is consistent with Fig. 14(b) and confirms that the 3D_NR_CP_ST01 case displays a much faster evolution towards e 3 than 3D_NR_CP_ST10, especially at the early stages of the evaporation process. Then, the evolutions of the various terms present in the averaged budgets of each component of the alignment vector are plotted in Fig. 15(b)-15(d) 5 . From the corresponding results, the following set of conclusions can be drawn:

1. In the two cases, it is found that the contribution TS 1 is negative, while TS 3 is positive. The corresponding evolutions are also similar in both cases and the strain-rate term TS seems to be the most important term in the scalar gradient orientation evolution towards the compression direction e 3 . 2. The contribution of the vorticity (TW) is less important than the one of the strain-rate and it seems to oppose to the alignment with the direction of compression (i.e., e 3 ). Instead, it drives the scalar gradient towards an alignment with the extension direction (e 1 ). Moreover, the term TW appears to be larger in simulation 3D_NR_CP_ST01 compared to 3D_NR_CP_ST10. This may be an outcome of the clustering effects, which are observed when varying the Stokes number value: liquid particles better follow the high-vorticity regions in case 3D_NR_CP_ST01.

3. The term GR is far from negligible, especially during the evaporation process. Its averaged contribution tends to drive the scalar gradient towards an alignment with the extension direction (e 1 ) whatever the evaporation process is proceeding or not. The resulting effects is thus opposite to those of strain. It is noteworthy that in simulation 3D_NR_CP_ST01, the droplets evaporate until time t = 0.1τ T is reached, and there is a strong decay of GR during this period. This confirms that droplets evaporation pilots the alignments towards e 1 or, at least, prevents the scalar to be aligned with the compression direction. This behaviour is fully consistent with the conclusions that were previously drawn in reference [12]. 4. It is more difficult to conclude as regards the impact of term WN, which represents the effects of the rotation of the eigenframe of S. This term is rather large only during the early stage of the evaporation process. Moreover, in case 3D_NR_CP_ST01, the sign of this contribution even changes at the end of the evaporation process. As a consequence, the resulting impact of this term is quite difficult to assess. Fig. 17: Temporal evolution of the transport equations for the three components of the orientation vector, i.e., the projection of the scalar gradient ξ in the eigenframe of strain-rate tensor, as given by Eq. 24. Solid line: case 3D_R_CP_ST10, dashed line: case 3D_R_CP_ST01.

In a second step of the analysis, the possible impact of the flame is now studied by considering reactive simulations 3D_R_CP_ST10 and 3D_R_CP_ST01. Attention is focused on the results obtained at time t = τ T = 0.2 ms, as previously shown in Fig. 3, at this time, the end of the evaporation process has been reached for case 3D_R_CP_ST01, and evaporation is almost finalized in case 3D_R_CP_ST10. The scalar ξ can thus be considered as a passive scalar. The PDF of the TSI term, orientation components, and eigenvalues of S obtained in the fresh reactants are provided in Fig. 16(a). The obtained results are almost identical to those associated to the non-reactive simulations, see Fig. 13(c): the flame kernel does not significantly alter the scalar and velocity fields in the fresh reactants. However, this conclusion does not hold for the PDF of the same quantities obtained in the burnt gases, which are reported in Fig. 16(b). Especially, the PDF of the orientations obtained in the burnt gases significantly differ from one case to another case, which may be ascribed to differences in the flame structure. Moreover, in comparison to those obtained in fresh reactants, the PDF of TSI and eigenvalues are much narrower as a consequence of the temperature (and resulting viscosity) increase, which significantly alters the turbulence kinetic energy. The influence of chemical reactions on the TSI term are further analyzed by computing conditional averages for different progress variable levels across the flame, as shown in Fig. 16(c).

The profiles of the conditional TSI has a parabolic shape with a minimum value obtained for c approximately equal to 0.5 in case 3D_R_CP_ST10, and equal to 0.4 in case 3D_R_CP_ST01. It is noteworthy that the values of the normalized TSI obtained in case 3D_R_CP_ST01 are mostly negative and significantly smaller than those obtained in case 3D_R_CP_ST10. This suggests that, inside the flame, the reactive flowfield topology favors the decay of the scalar gradient norm instead of enhancing it. The conditional averages of the orientations obtained in case 3D_R_CP_ST10, which are depicted with solid lines in the second column of Fig. 16(c), are relatively stable, while their counterparts obtained for case 3D_R_CP_ST01 display stronger variations. Within the range 0.2 < c < 0.6, an alignment of the scalar gradient with e 1 is favored in case 3D_R_CP_ST01, whereas the scalar gradient preferentially aligns with direction e 3 in case 3D_R_CP_ST10. These differences in the behaviour of the orientation vector explain those observed at the level of the normalized TSI term.

The budget of the orientation vector transport equation, Eq. 24, is analyzed and the corresponding results are reported in Fig. 17. The substantial derivative displayed in Fig. 17(a) confirms the tendency that was observed in Fig. 16(c). The budgets for each component of the orientation vector are plotted in the next three sub-figures. It seems that the effects of the turbulent flowfield, i.e., terms TS and TW, remain relatively stable across the flame: the strain (TS) tends to drive the alignment towards e 3 , while the vorticity (TW) acts to favor an alignment with e 1 . The behaviour of the term associated to scalar transport, i.e., GR, is more oscillatory. It is generally larger in case 3D_R_CP_ST01 compared to the case 3D_R_CP_ST10, especially in the vicinity of the burnt gases. The impact of the rotation of the eigenframe WN displays a similar behaviour across the flame. It is more significant in case 3D_R_CP_ST01 in the ranges 0.1 < c < 0.3 and 0.5 < c < 0.8, where WN is the dominant term of the budget and explains the change of alignment from e 3 to e 1 and then e 1 from to e 3 . Such a behaviour is similar to the one recently reported for fully premixed conditions by Zhao et al. [START_REF] Zhao | [END_REF]. The above analysis has been focused on the dynamics of the gradient of a passive scalar ξ that follows the fuel-air equivalence ratio variations. It would be also interesting to take a closer look at the interaction between a reactive scalar and the turbulent flowfield. To this end, we consider the progress variable c. Conditional averages of the normalized TSI term and alignments between the progress variable gradient and the eigenvectors of S are reported in Fig. 18. Compared to the behaviour of the passive scalar, see Fig. 16(c), the variation of both the TSI and alignments are more pronounced for the reactive scalar, especially in the case 3D_R_CP_ST10. The reactive scalar gradient is less preferentially aligned with the compression direction e 3 and more preferentially aligned with the extensive direction e 1 in the flame, which leads to a negative TSI contribution in both cases 3D_R_CP_ST01 and 3D_R_CP_ST10. Such a rotation of the orientation vector in the vicinity of the flame has been previously reported in experimental and numerical studies conducted under different conditions [START_REF] Chakraborty | [END_REF][START_REF] Chakraborty | Proceedings of the Combustion Institute[END_REF]52,26]. To assess the physical processes that are responsible for this change of orientation, the budget of the vector n is considered in Fig. 19. In contrast to the passive scalar behaviour that has been discussed above, WN is found to be the leading-order term while the other contributions (TS, TW, and GR) remains relatively stable and small in comparison to this contribution.

The cross-scalar dissipation rate, which is defined as

N ξc ≡ D∇ξ • ∇c = D|∇ξ| |∇c|(n ξ • n c
), is also an essential ingredient in the modelling of stratified or partially premixed turbulent flames [START_REF] Mura | [END_REF][START_REF] Bray | [END_REF]60,61,62,63]. The influence of the mixture fraction gradient on laminar flames has been early documented, e.g., [START_REF] Da Cruz | Proceedings of the Combustion Institute[END_REF][START_REF] Renou | [END_REF] and, from a general viewpoint, this may lead to a strong correlation between the two scalar quantities and associated gradients 6 . Thus, the cross-scalar dissipation rate cannot be simply deduced from the dissipation rates of both scalar and a modelled transport equation should be considered. In this respect, the DNS investigation conducted by Malkeson and Chakraborty in reference [63] confirmed that the cross scalar turbulence interaction (CSTI) term 6 One may however expect that the corresponding correlation is decreased for sufficiently large values of the velocity fluctuations RMS normalized by the laminar flame velocity S 0 L . is a leading order term in the cross scalar dissipation rate transport equation. The corresponding term is inspected below for the cases of turbulent spray flames 3D_R_CP_ST10 and 3D_R_CP_ST01.

Before proceeding with this analysis, the possible correlations between the reactive and passive scalars deserve further discussions. In this respect, it should be emphasized that strong correlations between the reactive scalar c and the passive scalar ξ as well as in the magnitude of their gradients (|∇ξ| and |∇c|) have been reported in a very recent study of the cross-dissipation rate in turbulent reactive two-phase flow [START_REF] Ahmed | Proceedings of the Combustion Institute[END_REF]. One may therefore expect that the analysis of the alignment between their gradients, n ξ • n c may also bring some interesting and useful insights into the behaviour of the cross scalar dissipation rate. The alignment characteristics indeed settles the sign of the cross-dissipation rate.

Conditional averages of n ξ • n c are presented in Fig. 20(b) for different values of the progress variable. It is worth noting that, within the flame (0.2 < c < 0.6), the two vectors are almost perfectly aligned (n ξ • n c = -1), especially in case 3D_R_CP_ST01. Since the orientations of ∇ξ and ∇c in the strain-rate eigenframe are quite similar in this progress variable range, as indicated by Fig. 16(c) and 18(b), such an alignment between the two vectors could have been anticipated. Moreover, such a conclusion is fully consistent with the PDFs of n ξ • n c that were previously reported by Malkeson and Chakraborty in reference [63].

Finally, the cross scalar turbulence interaction (CSTI) term is analyzed in its normalized form, -n ξ T •S• n c , see Fig. 20(a). One can see that this term remains relatively small and approximately constant across the flame in case 3D_R_CP_ST10, whereas it reaches a rather large (and positive) value within the flame in case 3D_R_CP_ST01. By expanding the CSTI term as follows: -n ξ T • S • n c = -i=3 i=1 n ξ,i λ i n c,i , one can see that this larger variation of the CSTI term in case 3D_R_CP_ST01 (compared to case 3D_R_CP_ST10) can be explained by the larger variations of the orientations observed across the flame in case 3D_R_CP_ST01 for both the passive and the reactive scalar gradients.

Discussion and conclusions

In this manuscript, turbulent combustion of a liquid spray is studied in constant volume conditions on the basis of a new set of computational data. Concerning the droplet dispersion and evaporation processes, the influence of the spray Stokes number is emphasized by considering the differences that are observed on the topology of the droplet clusters and the resulting scalar fields. The temporal evolutions of the energy spectra of the composition fluctuations are compared and analyzed for two distinct values of the Stokes number St. The impact of this spray Stokes number St on the turbulent flame development is also investigated by considering flame kernel growth simulations. The turbulent flames propagate faster in the spray characterized by a smaller Stoke number because of its faster evaporation and less segregated composition field. The flame kernel growth enhancement that may be induced by the consideration of a confined setup is highlighted thanks to a specific pressure variation model. This behaviour may be readily explained by the coupling existing between the heat released by the flame and the resulting fresh reactant temperature increase, which is induced by the uniform thermodynamic pressure increase. Further studies can be envisioned from this point, for instance by using the imposed pressure model (i.e., the PT model) so as to reproduce experimental pressure trace for future studies of flame propagation and/or auto-ignition.

In a second step of the analysis, turbulence-scalar interactions are studied for various droplet cloud combustion scenarios. Especially, the impact of the St number on the variation of the non-reactive scalar gradient orientation is inspected for different cases featuring both evaporation and combustion. To the best of the authors' knowledge, this is the first time that the alignment of the scalar gradient in the strain-rate tensor eigenframe is analysed through a direct inspection of the alignment transport equation in such conditions (i.e., with evaporation and combustion). The contributions of the strain, vorticity, scalar transport, and rotation of the strain-rate tensor eigenframe are quantified. In a first step of the analysis, it is found that the evaporation process tends to delay the alignment of the scalar gradient with the principal direction of compression e 3 . In reactive cases, provided that the fresh reactants are considered, the standard and preferred alignment of the non-reactive scalar gradient ∇ξ with e 3 as well as the behaviour of the TSI term do not seem to be altered by the flame kernel development. This is in contrast with the results obtained inside the flame and on the burnt gas side, which display significant departures from the behaviour observed in standard non-reactive and incompressible turbulence, especially in the cases featuring the smallest value of the Stokes number.

Such effects are more pronounced when a reactive scalar is under consideration. In fact, for sufficiently small values of the Karlovitz number, the evolution of the alignment that is observed across the turbulent premixed flame can be explained by the influence of the local premixed flame structure [27]. Indeed, as long as the turbulent kinetic energy remains moderate enough, the local strain-tensor is almost dominated by the thermal expansion effect and the local laminar flame structure imposes a strong correlation between the velocity and reactive scalar gradients. This leads to a preferential alignment of the scalar gradient with the principal direction of extension e 1 . Moreover, as far as the alignment vector is concerned, the rotation of the eigenframe has been recently shown to be essential in piloting the evolution of the orientation of the reactive scalar across a turbulent premixed flame brush [START_REF] Zhao | [END_REF].

As shown in Fig. 12, the alignment with e 1 is more pronounced in 3D_R_CP_ST01 than in 3D_R_CP_ST10 because the local flame is less perturbed by the two-phase flows processes in the conditions St=0.1. Indeed, in contrast to the flamelet regime of turbulent premixed combustion, it is quite less obvious to define a relevant local one-dimensional picture (featuring well-defined boundary conditions), from which the correlation between the scalar and velocity gradient can be estimated. Moreover, it should be acknowledged that the particle-source in cell (PSI) representation that is retained in the present study cannot provide all the details of the scalar-velocity interactions that takes place in the direct vicinity of the liquid-gas interface. Complementary numerical studies including the resolution of the liquid/gas boundary layer effects should therefore be performed to assess our understanding of the behaviour observed in reactive two-phase flows. This will be the subject of our future works. A first step in this direction has been reported in reference [12] but it is now necessary to take into account the effect of chemical reaction. Finally, it should be acknowledged that the present analysis has been conducted for moderate values of the turbulent Reynolds number. The corresponding values are relevant to Laboratory flames and practical devices operating at moderate levels of turbulence. It is however clear that some conclusions of the present study might be modified for higher Reynolds numbers and this offers an interesting direction for future complementary analyses.

∂ρ + ∂t + + ∂(ρ + u + j ) ∂x + j = ṁ+ p + ḋ+ ρ ρ ∞ u re f x -1 re f (26a) ∂ρ + u + i ∂t + + ∂ ρ + u + i u + j ∂x + j = - ∂p + 1 ∂x + i + ∂τ + i j ∂x + j + ṁ+ p u + i + ḋ+ ρui ρ ∞ u 2 re f x -1 re f (26b) ∂ρ + h + s ∂t + + ∂ ρ + h + s u + j ∂x + j = γ ∞ -1 γ ∞ d p + 0 dt + + ω+ T - ∂q + j ∂x + j + Ma 2 ∞ (γ ∞ -1) τ + i j ∂u + i ∂x + j + ḋ+ ρhs + ṁ+ p h + s ρ ∞ C p∞ T ∞ u re f x -1 re f (26c) ∂ρ + Y α ∂t + + ∂ ρ + Y k u + j ∂x + j = - ∂J α+ j ∂x + j + ḋ+ ρYα + ω+ α + ṁ+ p Y α ρ ∞ u re f x -1 re f (26d)
In the above equations, variables with superscript ( + ) are normalized quantities, and the subscript ( re f ) and ( ∞ ) denote either normalization (reference) parameters or values of physical quantities taken at infinity. The quantities that are reported at the end of each line (in brackets) correspond to values that are used to normalize the corresponding equations. A standard normalization procedure for low-Mach conditions is retained:

x i = x re f x + i ; u i = u re f u + i ; t = (x re f /u re f ) t + ; (27a) ρ = ρ ∞ ρ + ; ṁp = ρ ∞ (u re f /x re f ) ṁ+ p ; T = T ∞ T + ; (27b)

p = p ∞ p + = ρ ∞ (RT ∞ /M ∞ ) p + ; M = M ∞ M + ; h s = C p∞ T ∞ h + s ; (27c) γ ∞ = C p∞ C p∞ -R/M ∞ ; Ma ∞ = u re f γ ∞ T ∞ R/M ∞ ; (27d)
ωT = ρ ∞ C p∞ T ∞ (u re f /x re f ) ω+ T ; ωk = ρ ∞ (u re f /x re f ) ω+ k ; (27e) τ i j = ρ ∞ C p∞ T ∞ τ + i j ; q j = ρ ∞ C p∞ T ∞ u re f q + j ; J α j = ρ ∞ u re f J α + j ; (27f) ḋρ = ρ ∞ (u re f /x re f ) ḋ+ ρ ; ḋρui = ρ ∞ u re f (u re f /x re f ) ḋ+ ρui ; (27g) ḋρhs = ρ ∞ C p∞ T ∞ (u re f /x re f ) ḋ+ ρhs ; ḋρYα = ρ ∞ (u re f /x re f ) ḋ+ ρYα .

(27h) and the normalized EoS can be written as follows

p + 0 = ρ + T + M + ; (28a) h + s = 1 C p∞ T ∞ Nsp α=1 Y α T T0 C p,α (θ)dθ . ( 28b
)
On the one hand, the derivative analysis of Eq. 28a gives

d p + 0 = ρ + M + dT + + ρ + T + d 1 M + + T + M + dρ + = γ ∞ γ ∞ -1         ρ + r + dT + + ρ + T + Nsp α=1 r + α dY α + r + T + dρ +         , (29)
where

r + α ≡ 1 C p∞ R M α , r + ≡ 1 C p∞ Nsp α=1 RY α M α .
On the other hand, the derivative of the sensible enthalpy definition (i.e., Eq. 28b) reads

dh + s = C + p dT + + Nsp α=1 h + s,α dY α , (30)
Using Eq. 30 and considering r + /C + p = (γ -1)/γ, Eq. 29 becomes

γ ∞ -1 γ ∞ γd p + 0 = (γ -
where dh + s /dt + must be further expanded using the energy conservation Eq. 26c, because it includes a pressure derivative in its RHS. Then, we define ρ + d + s /dt + ≡ ω+ T -∂q + j /∂x + j + Ma 2 ∞ (γ ∞ -1) τ + i j ∂u + i /∂x + j + ḋ+ ρhs -ḋ+ ρ h + s , where it is noteworthy that the second term, which is associated to the dissipation function, can be neglected within the low Mach number framework.

The governing equation of pressure can be expressed as

γ ∞ -1 γ ∞ d p 0 dt = (γ -1) ρ d s dt + Nsp α=1 γr α T -(γ -1) h s,α ρ dY α dt + γrT dρ dt . ( 33
)
For the sake of simplicity, the super/subscript for normalized ( + ) and leading order ( 0 ) are no longer considered in the above expression. It is worth noting that, this equation is valid everywhere in the calculation domain. However, due to the spatial uniformity of the thermodynamic pressure, an averaged form of this equation over a volume V is preferred so as to avoid numerical difficulties. Finally, the thermodynamic pressure evolution will be determined from where α denotes the molecular diffusion velocity of species α. The above expression can also be written as follows: ξ p = (1 + n F • n O ) /2 and, since Fickian diffusion has been assumed, the unit vector n k , which gives the diffusion fluxes direction, may be expressed from the species mass fraction isolines: n α = ∇Y α /|∇Y α |. This allows to relate the premixedness index ξ p to the Takeno index G FO = ∇Y F • ∇Y O [67]:

d p 0 dt = γ ∞ γ ∞ -1 1 V V                       (γ -1)
ξ p = 1 2 1 + G FO |∇Y F | • |∇Y O | (36)
In one-dimensional flame structures described with single-step chemistry, the above index is zero in non-premixed conditions while it is unity in premixed conditions [67]. In contrast to the standard definition of reference [67], which was based on single-step chemistry, Y F and Y O do not presently correspond to the mass fractions of fuel and oxidizer. The definition has been indeed generalized herein by considering the mass fractions of unburnt carbon and oxygen atoms

Y F ≡ α CO2,H2O α M C M K Y α , Y O ≡ α CO2,H2O α M O M K Y α (37)
Appendix C: Derivation of transport equations of scalar alignments

For a scalar field ξ(x, t) the transport equation of which can be written as

∂ξ ∂t + u • ∇ξ = 1 ρ ∇ • (ρD∇ξ) + ωξ (38)
with D the diffusion coefficient of the scalar and ωξ the source term, the transport equation for its gradient g ≡ ∇ξ can be obtained

D g Dt = -S g - 1 2 g × ω + ∇(RHS ξ ) ( 39
)
where S is the strain tensor, ω is the vorticity and RHS ξ ≡ ∇ • (ρD∇ξ)/ρ + ωξ denotes the right hand side (RHS) of the transport equation of the scalar. Then, the transport equation for its norm direction vector n ≡ g/| g| can be deduced from Eq. ( 39)

Dn Dt = 1 | g| IN D g Dt = -INSn - 1 2 IN (n × ω) + IN |g| ∇(RHS ξ ) (40)
where the matrix IN is defined as follows

IN ≡ I -nn T (41)
with I the identity matrix. The projection of n in the eigenframe of S can be characterized with the vector

n ≡           n1 n2 n3           =           n • e 1 n • e 2 n • e 3           = R T n ( 42
)
with R is the eigenmatrix of S R ≡ [e 1 |e 2 |e 3 ] ( 43)

where e i are the eigenvectors of S corresponding to eigenvalues λ i with a descending order, i.e., λ 1 > λ 2 > λ 3 . Finally, the transport equation for n can be achieved

D n Dt = D(R T n) Dt = R T Dn Dt + DR T Dt n = -R T INSn - 1 2 R T IN (n × ω) + R T IN |g| ∇(RHS ξ ) + DR T Dt n = n T Λ n I -Λ n - R T 2 (n × ω) + I -n n T |g| R T ∇(RHS ξ ) + DR T Dt n ( 44
)
where Λ is the diagonal matrix based on the eigenvalues of S,

Λ ≡           λ 1 λ 2 λ 3           . ( 45
)
Since we have

n = 3 j=1 n j e j , (46)
the last term in Eq. ( 44) can be expanded as follows

DR T Dt n = W n ( 47
)
where the anti-symmetric tensor W denotes the rate of rotation of the principal axes

W i j = De i Dt • e j . (48)
In fact, this tensor can be evaluated from the transport equation of the strain-rate tensor in eigenframe (D(R T SR)/Dt), it reads

W i j = De i Dt • e j = 1 -δ i j λ i -λ j       - 1 4 ωi ω j + 1 2 R T ∇ Du Dt + ∇ Du Dt T R ji       (49)
where δ i j is the Kronecker delta tensor and the substantial derivative of the velocity is given by

Du Dt = 1 ρ (-∇p + ∇ • τ + f ) . (50)
In summary, the transport equation for the projection of n in the eigenframe of the strain-tensor reads: where the vector TW and GR are defined as

Dn 1 Dt =         3 j=1 λ j n2 j -λ 1         n1 + TW 1 + ñ2 W 12 + n3 W 13 + (GR) 1 (51a)
TW ≡ - R T 2 (n × ω) , GR ≡ I -n n T |g| R T ∇(RHS ξ ) . (52)

Compliance with ethical standards

Funding: As stated in the acknowledgements section, this work has been funded by the CAPA Program (CNRS, ENSMA, SAFRAN, and MBDA) of the ANR (Agence Nationale de la Recherche). Conflict of Interest: The authors declare that they have no conflict of interest.

Fig. 1 :

 1 Fig. 1: Illustration of the typical range of length scales involved in a droplet cloud combustion scenario.

Fig. 2 :

 2 Fig. 2: Visualization of droplet clusters in the middle cut-plane after dispersion by the constant-energetics physical-space forced HIT. Left: case 3D_NR_CP_ST10 (St = 1.0). Right: case 3D_NR_CP_ST01 (St = 0.1).

Fig. 3 :

 3 Fig. 3: Time evolutions of the average and RMS values of the iso-octane mass fraction Y C8H18 . Solid line: case 3D_NR_CP_ST10, dashed line: case 3D_NR_CP_ST01. Left: normalized mean value. Right: normalized RMS.First row: evolutions plotted versus time t (physical units). Second row: evolutions plotted versus t/τ a with τ a the evaporation time scale evaluated from the initial characteristics of the two-phase flow.

1 EFig. 4 :

 14 Fig. 4: Spectral density of the energy of scalar fluctuations (normalized by the total energy) plotted versus the wave number K (normalized by L x ). Solid line: 3D_NR_CP_ST10. Dashed line: 3D_NR_CP_ST01.

Fig. 5 :

 5 Fig. 5: PDFs of 1/S p . Solid line: 3D_NR_CP_ST10, dashed line: 3D_NR_CP_ST01, abscissa axis is in log scale.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Turbulent flame kernel growths in droplet sprays. Four rows corresponding (from top to bottom) to four distinct times t = 0.0, τ T , 3.0 • τ T and 4.0 • τ T ; four columns corresponding (from left to right) to simulations 2D_R_CP_ST10, 2D_R_CV_ST10, 2D_R_CP_ST01, and 2D_R_CV_ST01. The flowfield is colored by the equivalence ratio and the flame front is visualized by the heat release rate.

Fig. 8 :

 8 Fig. 8: Left: temporal evolution of the equivalent radius, right: PDF of local flame displacement speed at t = 2τ T =0.4 ms normalized by the reference laminar flame speed S 0 L = 1.19 m/s.

Fig. 9 :Fig. 10 :

 910 Fig. 9: Scatterplots of the magnitude of the progress variable gradient versus the progress variable at time t = 2τ T . Red curves corresponds to conditional averages |∇c| c .

Fig. 11 :

 11 Fig. 11: Snapshots of the three-dimensional simulations of the turbulent flame kernel growth obtained at t = τ T . In the tridimensional visualisations (a,b), the flame front is put into evidence with the iso-surface c = 0.5 colored by the local HRR. The blue surface is an iso-contours of the Q-criterion (Q = 5 × 10 8 s -2 ) colored by the equivalence ratio. Colored contours of the HRR in a cut-plane are provided in (c) and (d).

Fig. 12 :

 12 Fig. 12: PDFs of the premixedness index (a), PDFs of the flame displacement speed (b), and scatterplots of the magnitude of the progress variable gradient plotted versus the progress variable (c,d) at time t = τ T for the three-dimensional reactive simulations.

3 (Fig. 13 :

 313 Fig.13: PDFs of the normalized TSI terms (left), orientation of the scalar gradient in the eigenframe of the strain-rate tensor (middle), and eigenvalues of the strain-rate tensor (right) issued from non-reactive simulations at different times. Solid line: 3D_NR_CP_ST10, dashed line: 3D_NR_CP_ST01.

3 Fig. 14 :

 314 Fig. 14: Temporal evolution of the averaged normalized TSI term (a) and orientation of the scalar gradient in the strain-rate eigenframe (b). Solid line: 3D_NR_CP_ST10, dashed line: 3D_NR_CP_ST01.

D n 1 3 (

 13 /Dt D n 2 /Dt D n 3 /Dt (a) Temporal evolution D n /Dt d) Budget of D n 3 /Dt

3 ( 3 (Fig. 16 :D n 1 /

 33161 Fig.16: PDFs of the normalized TSI terms (left), PDF of the orientations between the scalar gradient and eigenvectors of strain-rate tensor (middle), and PDF of the eigenvalues of strain-rate tensor (right) in the fresh reactant and burnt gases (a,b). Profiles of their conditional averages plotted versus the progress variable (c). Solid line: 3D_R_CP_ST10, dashed line: 3D_R_CP_ST01.

Fig. 18 :

 18 Fig. 18: Conditional averages at time t = τ T of the turbulence-reactive scalar interactions (a) and alignments of the reactive scalar gradient the eigenframe of the strain-rate tensor (b). Solid line: 3D_R_CP_ST10, dashed line: 3D_R_CP_ST01.

Fig. 19 :

 19 Fig.19: Conditional averages at time t = τ T of the substantial derivative of the alignments of the progress variable gradient in the eigenframe of the strain-rate tensor across the flame front (a), and budgets of its transport equation (Eq. 24) (b,c,d). Solid line: 3D_R_CP_ST10, dashed line: 3D_R_CP_ST01.

Fig. 20 :

 20 Fig.20: Conditional averages at time t = τ T of the cross scalar turbulence interaction (a) and alignments between the gradients of the reactive scalar c and the non-reactive scalar ξ (b). Solid line: 3D_R_CP_ST10, dashed line: 3D_R_CP_ST01.

F

  T -(γ -1) h s,α ρ dY α dt (II)

  TW 2 + n1 W 21 + n3 W 23 + (TW 3 + n1 W 31 + n2 W 32 + (GR) 3 (51c)

Table 1 :

 1 Liquid iso-octane properties

	ρ l (kg/m 3 ) C p,l (kJ/kg.K)	P cr (bar) T cr (K)	T b (K)
	690.0	2122.8	25.7	543.9	372.4

Table 2 :

 2 Main characteristics of the numerical simulations. CP and CV denote constant pressure and constant volume model, respectively (see Eq. 13 to Eq. 16).

		L x (mm)	N x	St number Reactive Pressure Model
	3D_NR_CP_ST10	4.000	256	1.0	No	CP
	3D_NR_CP_ST01	4.000	256	0.1	No	CP
	2D_R_CP_ST10	4.000	256	1.0	Yes	CP
	2D_R_CV_ST10	4.000	256	1.0	Yes	CV
	2D_R_CP_ST01	4.000	256	0.1	Yes	CP
	2D_R_CV_ST01	4.000	256	0.1	Yes	CV
	3D_R_CP_ST10	3.125	200	1.0	Yes	CP
	3D_R_CP_ST01	3.125	200	0.1	Yes	CP

Table 3 :

 3 Initial flowfield characteristics (carrier phase)

	T 0	700.0 K	Initial temperature
	p 0	5.0•10 5 Pa	Initial pressure
	ν	1.3•10 -5 m 2 /s Gas dynamic viscosity
	u rms	1.5 m/s	Velocity fluctuation RMS
	l t	0.300 mm	Integral length scale
	η k	0.023 mm	Kolmogorov length scale
		7790 m 2 /s 3	Turbulence dissipation rate
	τ T	0.200 ms	Integral time scale
	τ k	0.041 ms	Kolmogorov time scale
	Re λ	27.0	Taylor Reynolds number
	Re T	34.0	Turbulent Reynolds number

  1) ρ + dh + s + + -(γ -1) h + s,α ρ + dY α + γr + T + dρ +

				Nsp
					γr + α T (31)
				α=1
	from which one may obtain
	γ ∞ -1 γ ∞	γ	d p + 0 dt + = (γ -1) ρ + dh + s dt + +	Nsp α=1	γr + α T

+ -(γ -1) h + s,α ρ + dY α dt + + γr + T + dρ + dt + ,

The droplet diameter is approximately 6.1 µm for cases associated to St = 1.0. It is approximately 1.9 µm for those corresponding to St = 0.1, the values of the kinetic and evaporation relaxation times are thus significantly smaller (i.e., 10 times smaller) than those associated with simulations conducted at St = 1.0

The retained definition of the premixedness index is provided in Appendix B.

A similar behaviour was also put into evidence in the theoretical work of Gonzalez and Paranthoën[53,54], which was devoted to the effects variable mass density on the kinematics of scalar gradient.

Once normalized by its value at stoichiometry, i.e., y CO,st , the carbon atom to oxygen atom mass fractions ratio, i.e., y CO /y CO,st , can be thought as the equivalence ratio.

The transport term TR is sufficiently small to be discarded from the present analysis.

Acknowledgements This work is funded by the CAPA Program (CNRS, ENSMA, SAFRAN, and MBDA) of the ANR (Agence Nationale de la Recherche). A part of this work has been performed using high-performance computing (HPC) resources from GENCI (Grand Equipement National de Calcul Intensif)-[CCRT/CINES/IDRIS] under Grant A0052B07456.

Appendix A: Thermodynamic pressure variation with artificial mass sources and chemical heat release

The thermodynamic pressure evolution is deduced from the analysis of the normalized EoS together with a non-dimensional form of the NSE featuring an artificial mass source term ṁp , see Eq. 2 in section 2. Since Asphodele resolves the non-dimensional form of NSE, the thermodynamic pressure evolution is considered within the corresponding non-dimensional framework:

Appendix B: Evaluation of the premixedness index

The premixedness index is evaluated from the generalized definition introduced in reference [66]: