Investigations of Evaporative Cooling and Turbulence Flame Interaction Modeling in Ethanol Turbulent Spray Combustion Using Tabulated Chemistry

Evaporative cooling effects and turbulence flame interaction are analyzed in the large eddy simulation (LES) context for an ethanol turbulent spray flame. Investigations are conducted with the artificially thickened flame (ATF) approach coupled with an extension of the mixture adaptive thickening procedure to account for variations of enthalpy. Droplets are tracked in a Euler–Lagrangian framework, in which an evaporation model accounting for the inter-phase non-equilibrium is applied. The chemistry is tabulated following the flamelet generated manifold (FGM) method. Enthalpy variations are incorporated in the resulting FGM database in a universal fashion, which is not limited to the heat losses caused by evaporative cooling effects. The relevance of the evaporative cooling is evaluated with a typically applied setting for a flame surface wrinkling model. Using one of the resulting cases from the evaporative cooling analysis as a reference, the importance of the flame wrinkling modeling is studied. Besides its novelty, the completeness of the proposed modeling strategy allows a significant contribution to the understanding of the most relevant phenomena for the turbulent spray combustion modeling.


Introduction
Turbulent spray combustion features the unsteady interaction of many strongly coupled physical phenomena varying in a broad range of time and length scales. The complexity resulting from such interacting phenomena can be appraised with the studies presented for instance by Jenny et al. [1], Gutheil [2], and Sacomano Filho [3]. Heat and mass transfer between phases and the effects of the turbulence on a flame, as well as their reciprocal, stand out among these coupled phenomena. Although both interactions are strongly coupled, they are not always comprehensively addressed in numerical simulations. The cooling down of the gas mixture caused by the evaporation process, namely the evaporative cooling (EC), is often neglected [4][5][6][7][8] or simplified [9,10] in the context of tabulated chemistry. Herein, the consideration of heat losses is not trivial and typically involves high computational efforts [11][12][13]. Commonly, the low volume fraction of diluted sprays are used to justify simulations without heat losses [14,15]. Nevertheless, as analyzed in [16] for one-dimensional flames propagating in droplet mists, the inclusion of the evaporative cooling is quite relevant in the spray combustion modeling. Yet the importance of the proper modeling of the turbulence-flame interaction is However, further investigations are still necessary to turn this approach more consistent for general turbulent spray flame simulations. On the other hand, analysis of the influence of the model parameter β are needed for spray flames, which have not been comprehensively addressed in cases where heat losses are considered.
In the present work, the effectiveness of the evaporative cooling and the FSW modeling to characterize a diluted turbulent spray flame is analyzed. The chemistry tabulation method proposed by van Oijen and de Goey [11] and comprehensively tested for spray flames in [16] is employed here. Namely, effects of the evaporative cooling interfere with reaction rates, and the consideration of heat losses are not limited to the evaporative cooling in the flamelet generated manifold (FGM) table.
Simulations are conducted following the large eddy simulation (LES) approach coupled with the ATF method. Besides the novelty associated with such a consistent consideration of heat losses in a turbulent spray flame, the mixture adaptive thickening approach [21] is extended to include effects of enthalpy variations. The flame EtF5 of the Sydney diluted spray flame burner [31] is selected as a benchmark to the subsequent analysis. Detailed chemistry effects are included by 57 species and 379 intermediary reactions mechanism of ethanol-air proposed by Marinov [32] by means of the FGM methodology. The unsteadiness arising from the turbulent dispersion of evaporating droplets are captured by a Euler-Lagrangian spray module relying on the LES approach. An evaporation model accounting for the inter-phase non-equilibrium is applied to describe the droplet evaporation process. A parametric study has been performed to evaluate the effectiveness of the evaporative cooling and the turbulence flame interaction modeling consideration.
The remaining structure of this paper is divided into three parts. An overview of the theoretical background is described in Section 2. In Section 3, simulation results are presented. Analyses of effects of the evaporative cooling and its coupling with the flame surface wrinkling modeling are systematically addressed. In the last part, final remarks and the main conclusions are summarized.

Materials and Methods
A Euler-Lagrangian approach was adopted to represent the two-phase flow. Herein, a full inter-phase two-way coupling was accounted for. Carrier gas-phase quantities were interpolated into droplets positions, while influences of the dispersed phase were introduced through source terms in computational cells.

Gas Phase
The turbulent motions of the carrier phase are described in the LES context following a variable-density low Mach number formulation. According to this approach, mass and momentum equations are described by ∂ρ ∂t ∂ρ u i ∂t The dependent filtered variables were obtained from spatial filtering as ψ = ψ + ψ" with ψ = ρψ/ρ. Over-bars and tildes express spatially filtered and density-weighted filtered values with a filter width ∆ mesh , respectively, while double prime represents sub-grid scale (SGS) fluctuations. ρ is the mixture density , t time, u j components of velocity in j (j = 1, 2, 3) direction, p pressure, x j Cartesian coordinate in j direction, µ the dynamic viscosity, g i the component of gravitational acceleration, δ ij the Kronecker's delta, and S ij the strain rate. The term S m corresponds to the introduction of mass from the droplets into the carrier phase, while S u,i is the source term of momentum due to the presence of the dispersed phase. Both follow the implementations presented by Chrigui et al. [33]. The SGS stress tensor τ sgs ij is closed by means of the Smagorinsky model with the dynamic procedure of Germano et al. [34]. More details about the mathematical treatment given here to mass and momentum equations can be found in Sacomano Filho et al. [21].

Mixture Formation and Combustion Modeling
In order to account for the evaporative cooling and general heat exchanges, three scalar quantities are used to characterize the mixture following the FGM method: the mixture fraction Z, the reaction progress variable Y pv , and the absolute enthalpy of the gas mixture h. The transport equation for them can be written in terms of a general variable ψ within the ATF modeling as where µ t is the turbulent viscosity. For the transport equations of Z and Y pv , σ ψ and σ t,ψ respectively represent the laminar Sc and the turbulent Sc t Schmidt numbers (Sc = Sc t = 0.7), while for h both respectively correspond to the laminar Pr and the turbulent Pr t Prandtl numbers (Pr = Pr t = 0.7). It is important to highlight that, with these values for Sc and Pr the unitary Le = Pr/Sc approach is maintained. The quantity F corresponds to the thickening factor, E * ∆ to the efficiency function, and Ω denotes the flame sensor. Details about these quantities are addressed below.
The termω ψ =ω ψ Z, Y pv , h corresponds to the reaction rate for the Y pv whereas it is set to zero for the mixture fraction and absolute enthalpy equation. Similarly to the reaction progress variable source term, ρ = ρ Z, Y pv , h and µ = µ Z, Y pv , h are obtained from the employed FGM table. Therefore, these are expressed as functions of the transported scalar quantities Z, Y pv , and h. The source term S ψ consists of the source of vapor introduced by the dispersed phase in the transport equation for Z, specifically S Z = S m . Considering that the mass fraction of fuel is not present in the combination used to define the reaction progress variable (see Equation (15)), as well as no isolated droplet burning model is included in the employed approach, S ψ does not contribute (therefore it is set to zero) for the transport of Y pv . For the transport equation of h, S h is given by whereṁ d,p denotes the mass of vapor released by the parcel p into the control volume V, m d,p is the droplet mass, N p the number of real droplets in parcel p, N the total number of tracked parcels, c l the specific sensible heat of liquid, T ref a reference temperature (298 K), T τ d the droplet temperature at time step τ, h f the formations enthalpy, and L v the heat of vaporization. More description about how the phase coupling source terms are computed in the ATF context can be found in [3,21].
Following the ATF method, the flame thickening is performed by means of a dynamic procedure. Accordingly, only the flame region is thickened and no interferences of the ATF with the pre-vaporization zone occur. The flame sensor Ω used by Aschmoneit [35] to simulate partially premixed flames is employed here. The quantityω pv,max Z, h in Equation (5) is associated to the maximum value of the source termω pv Z, Y pv , h at the same mixture composition and enthalpy level.
where Ω c is the flame sensor proposed by Durand and Polifke [36] as implemented in [29]. In this context, a progress variable in its normalized form is needed, which is defined as c = Y pv /Y eq pv . Herein, Y eq pv denotes the equilibrium value of Y pv for a specific mixture composition.
In order to facilitate the access ofω pv,max during the simulation time, an extra look-up table is stored with its values mapped on Z and h. According to the dynamic method, the thickening factor F is defined as F = 1 + (F max − 1)Ω. (6) in which F max is defined according to an extension of the mixture adaptive thickening proposed in [21] to account for enthalpy variations as in which V denotes the cell volume and ∆ x,max is the maximum cell size necessary to capture the laminar flame speed with less than 15 % error in one-dimensional simulations [29]. In agreement with Kuenne et al. [29], ∆ x,max ≈ 0.3δ 0 l ( Z, h). Additionally, two safety factors are also considered to this maximum cell size definition. The first of 1/ √ 3 is included to account for the worst scenario of flame propagating in cell diagonal direction for three-dimensional cases, the second 1/α is used to account for the non-uniformity of cells in F. Applying both safety factors, the maximum cell size is defined here as with α = 1.0, since cubic cells are employed in the jet core region (see Section 2.3).
To recover the modifications introduced by the ATF on the turbulence-flame interaction while including the effects of the unresolved FSW, the efficiency function E ∆ follows the fractal model. It is written as where ∆ = F · δ l is the combustion LES filter width, δ 0 l is the laminar flame thickness defined from the temperature gradient, β a model exponent, and δ l an equivalent laminar flame thickness in terms of Gaussian LES filter as proposed in [37] (see more details in the sequel). This is a reasonable approach in the simulations presented here, since we could notice that already for flames presenting a smaller Reynolds number (Re) in [3]. Considering the minimum operator in this inequality shows that the saturation holds. Equation (9) is connected with the efficiency function E * ∆ by means of Equation (12) is used to enforce that E * ∆ → 1 outside the flame region, since an efficiency function is no longer necessary in this part of the domain. Hence, the term FE * ∆ µ/σ ψ in Equation (3) reduces to µ/σ ψ , whereas (1 − Ω)µ t /σ t,ψ becomes µ t /σ t,ψ . Since σ = σ t = 0.7, both terms are merged in which the molecular viscosity is added with the turbulent one and the effective viscosity µ eff = µ + µ t is obtained. In contrast to this scenario, within the flame region Ω → 1 and the term (1 − Ω)µ t /σ t,ψ is canceled, where molecular and turbulent transport in the flame are modeled by FE ∆ µ/σ ψ .
The determination of ∆ follows a formulation that can be directly applied for future implementations accounting for the dynamic modeling of β. Accordingly, the equivalent laminar flame thickness used to define the LES combustion filter is given by δ l = ζδ 0 l (in agreement with [37]), in which ζ is an adjusting factor to enforce that the efficiency function goes to one for laminar flames. Values of ζ are supposed to vary according to the mixture state, however only the variations associated with the composition are accounted for in this work. As a result, the actual formulation of the mixture adaptive thickening which accounts for enthalpy variations is considered to achieve In regions where the cell sizes deliver F greater than 1, Equation (13) reduces to

Chemistry
To construct our non-adiabatic FGM table, the strategy proposed by van Oijen and de Goey [11] is applied. Particularly, the reference table described in [16] is employed. It is based on the combination of freely propagating adiabatic flamelets for a mixture composition span of Z ∈ [0.050, 0.194] with fresh gas temperature varying from 250 K to 900 K with burner stabilized and extrapolated flamelets. The maximum value of this temperature determines the upper enthalpy value, while its minimum defines the lowest temperature (T low ) allowed in the computational domain. Interpolations are done [38] between the limiting values of Z up to pure air and ethanol at 300 K, which corresponds to the ambient temperature of the investigations performed in this work. To compute both kinds of flamelets, the chemical mechanism proposed by Marinov [32] is employed. It represents the oxidation of ethanol in the air by means of 57 species and 379 intermediate reactions.
The inclusion of burner-stabilized and extrapolated flamelets in our table allows a more universal description of the heat losses. Burner-stabilized flamelets mimic premixed flames burning attached to a porous medium with controlled temperature. By keeping the temperature of this medium constant and adjusting the inflow velocity, heat losses can be introduced in the domain resulting in different enthalpy levels for the flame. In this particular step, the temperature is set to the same value as for the unburned gas of the coldest adiabatic flamelet. When the lowest enthalpy levels with the burner stabilized flamelets is reached, the burnt gas is still far away from T low . To cover the remaining range of physical states, step-wise extrapolations of the thermo-chemical data are done to the lowest enthalpy level. This is computed by assuming the equilibrium mixture composition at T low . Finally, once a manifold is generated, the look-up table is build up over the controlling variables Z, Y pv , and h. Following the strategy of Ketelheun et al. [12], an additional table is created to facilitate the establishment of the boundary conditions in terms of temperature.
The definition of Y pv follows the combination of mass fractions presented in [16]. This was empirically defined in view of a valid progress variable for the full range of enthalpy levels encountered in the employed table. Namely, a monotonically evolving variable from fresh to burnt gases along all flamelets used to construct the employed FGM table. The Y pv can be written as where M k is the molar mass of species k. For more details about the definition of a progress variable in FGM context, the reader is referred to [11]).

Liquid Phase
The computation of the motion of non-rotating droplets considers only drag and buoyancy forces. Once that the density ratio of liquid ethanol and the gaseous mixture has an order of 10 3 , complementary forces are neglected.
Heat and mass exchanges are computed following the model proposed by Miller et al. [39] with the correction procedure derived in [21] to include effects of the flame thickening onto the dispersed phase. Both phenomena are described by and with T the temperature, Nu the Nusselt number, f 2 a correction factor due to evaporation computed as in [39], Pr the Prandtl number, θ 1 the ratio of the gaseous and liquid specific sensible heat (c p /c l ), τ p = ρ p d 2 p /18µ expresses the particle relaxation time,ṁ p = dm p /dt, Sh the Sherwood number, and H M represents the specific driving potential for mass transfer according to [39]. It is important to highlight that vapor pressures, necessary to determine H M are computed here with the Wagner equation as in [40].
For the vast majority of droplets tracked in our simulations, the relaxation time exceeded the scales within the subgrid. As a consequence, SGS turbulent dispersion and micro-mixing modeling for the evaporation process had little effect.

Experimental Configuration and Numerical Setup
The flame EtF5 of the Sydney diluted spray flame burner [31] is used as a benchmark for the performed LES. This flame belongs to a series of flames (EtF1, EtF2, EtF5, and EtF7), where the premixed flame mode is dominant. This is an important feature in view of the chosen modeling approach, which has been originally developed for premixed flames. EtF5 was produced by a round jet burner with a diameter D of 10.5 mm fueled with ethanol. The flame stabilization was achieved with a series of pilot flames disposed of in a concentric arrange to the main jet nozzle with an annular diameter of 25 mm. A co-flow surrounded both streams with an annular diameter of 104 mm and a velocity of 4.5 m/s, which corresponds to the air velocity of the wind tunnel where the burner is mounted. The bulk velocities of main jet and pilot were 48 m/s and 11.6 m/s, respectively. At the nozzle exit, a mixture of droplets, air, and fuel vapor was formed with a Re of 38,200 and an equivalence ratio φ = 0.15 for the gaseous phase, which is found below the lower flammability limit (φ low ).
The computational domain is constructed with a Cartesian grid in order to keep the cell growth rate smooth and to avoid large values of the factors F and α (see Equations (3) and (14), respectively). Specifically, mean values of F are about 7.0 but no bigger than approximately 9.0, while α = 1.0 due to the cubic cell arrangement in the jet core region. When compared to the "O-grid"arrangement (as employed in [3]), which would be an alternative to run in the used block-structured code, the Cartesian arrangement avoids the high cell diagonal growth rate. Figure 1 presents a sketch of the employed grid, which is composed of approximately 5.4 million cells distributed over 468 blocks. In this figure black lines represent block limits, whereas gray lines are cell sides. The dashed orange box indicates the jet core region, in which cells have uniform size with a side of 0.75 mm. A contour plot of the instantaneous value ofω pv is included to show that the main reaction zone is entirely found within the jet core. The computational domain starts at the exit of the experimental burner and extends up to 70 D (735 mm). The most upstream part of this domain has a diameter of 10 D that linearly grows up its maximum with 20 D. The employed grid can be seen as an improvement of the coarse grid used in [3], since it preserves similar characteristics of the previous one (e.g., longitudinal and radial cell sizes are the same, however with higher homogeneity here). In [3] a comprehensive grid study has been conducted for reactive and non-reactive flows.
Regarding the boundary conditions, attention is given to the description of the pilot flame. This flame is fueled by a mixture of acetylene and hydrogen at the hydrogen to carbon proportion of the stoichiometric reaction of ethanol. Under adiabatic assumptions, such a boundary condition can be easily represented by adjusting the mixture composition to the stoichiometry. This procedure is previously adopted for simulations presented in [3] and it is commonly employed by diverse research groups (e.g., [7,33,41]) which also used data from the Sydney diluted spray flame burner for modeling validation. Concerning a non-adiabatic formulation, the introduction of the enthalpy transport allows a more accurate representation of the pilot flame. It is important to notice that, despite the same proportion of hydrogen and carbon elements, the presence of the hydroxyl (OH) into the ethanol molecules does not allow the mixture of acetylene and hydrogen to produce the same composition of products. This characteristic becomes clear comparing the stoichiometric reaction of ethanol and the corresponding mixture of acetylene and hydrogen with air Since the main function of a pilot flame consists in the production of hot products to stabilize the main flame, the mixture fraction is adjusted to match the products ratio generated by Equation (19). Namely, the ratio formed by the sum Y CO 2 + Y H 2 O considering only products. As a result, the following reaction is used to define the mixture composition of the pilot flame for the non-adiabatic simulations performed here Once that the absolute enthalpy of the reaction given by Equation (20) does not match with that obtained from Equation (19), the reactants are lightly pre-heated in order to pair up with the desired value. It is important to highlight that with a tabulation method where heat losses are exclusively linked to the evaporation process as in [13,22], such a modeling feature would not be realizable. The relevance of the specification for this boundary condition has been analyzed in [3] for the flame EtF2.
The turbulent boundary condition follows the synthetic turbulence generator proposed by Klein et al. [42]. To account for the pre-vaporized mass of fuel exiting the injection's tube, the mixture fraction receives the explicit value of 0.0165 (corresponding to φ over = 0.15). Effects of the heat transfer throughout the injection tube are considered. For that, the mixture composed by dispersed fuel droplets and air entering at 300 K into the burner exchanges heat with the injections tube wall at 300 K, which results in a mixture exiting the nozzle with 298 K. This arises from our assumption that the mixture of vapor and air with φ = 0.15 and Re = 38,200 results from the evaporation of liquid in an air stream where both were initially at 300 K. This mixture was then heated up by tube walls at 300 K assuming Re = 38,200. Considering the remaining boundary conditions for the enthalpy, walls are treated as adiabatic and incoming streams of air are kept at 300 K.
With respect to the boundary conditions of the diluted phase (i.e., droplet size distributions, conditioned velocities on size, injected mass, and injection position), these have been extracted from experimental measurements presented in [31]. An exception occured for the radial component of droplets velocities, for which the procedure suggested by De and Kim [6] was employed here. More details about this procedure are presented in Section 3.1. The initial droplet temperature has been set to 298 K, namely droplets and the pre-vaporized mass of fuel exiting the injection's tube had the same temperature.
Calculations for this research were conducted with the coupled version of the academic software FASTEST (see [3]). Different from previous works, some improvements in numerical methods and algorithms are included. For instance, the ordinary differential equations of the evaporation modeling are integrated in time here with the fifth-order Runge-Kutta scheme of [43].
Simulations are initialized from reactive single-phase computations obtained after 0.055 s, which coincides with one flow-through-time (FTT). The FTT of this configuration corresponds to the time necessary for the gas flow to fill a volume defined from the droplets injection position to 30 D downstream the burner exit, where the last measurement plane is set. Namely, a cylinder with a diameter of 10 D (to account for the burner coflow) and a length of 30 D. Thereafter, parcels start to be injected and a stable number of them is achieved after approximately 1 FTT. Once a stable number of parcels is reached, simulations run for an additional FTT until the averaging begins. Statistics (mean and rms) are collected for 2 FTT for all cases, which amounts circa 35,300 core-hours on 80 CPUs and corresponds to a reduced computational time (RCT) of ≈ 7.11 × 10 −8 s [RCT = elapsed time × number of CPU/(number of cells × number of iterations)].
To perform the time integration, a second-order, three steps Runge-Kutta scheme (see details in [44]) has been employed. The time-step is maintained constant at 2.4 × 10 −6 s. Particles are coupled at each three-time steps. As a result, particles were tracked in a coupling time step of 7.2 × 10 −6 s. Despite this discrepancy, the corresponding Courant-Friedrichs-Lewy (CFL) number is kept below 1 respectively for droplets and the gas phase. This coupling strategy alleviates the computational costs of the particle tracking, which amounts in our calculations approximately 50% of the entire computational time.

Results
The evaluation of the proposed modeling strategies starts with the establishment of a reference scenario. This corresponds to simulation performed with a fixed exponent β (i.e., β = 0.5) of Equation (9) and the inclusion of evaporative cooling effects. The value 0.5 resembles the one used in the first investigations of Charlette et al. [28] and is a typical value assumed by different authors (e.g., [9,20,29,44]), therefore it is chosen as a reference here. Altogether three cases are studied, besides the mentioned reference case, one neglecting effects of the evaporative cooling and another considering a parametric variation of the model parameter β while including EC effects. In this last case, the value of β is defined as 0.2.
The influence of the evaporative cooling onto the flow is expressed in terms of the heat necessary to vaporize an amount of liquid (i.e., that associated with the latent heat of vaporization) and the heat-up of droplets. In order to investigate such an influence, EC effects are neglected by removing the terms associated with the latent heat of vaporization and the heat-up of droplets from Equation (4). As a result, only the term including the vapor formation's enthalpy is preserved The source term presented in Equation (21) is necessary for our approach since the energy equation is expressed in terms of the absolute enthalpy. In cases where sensible enthalpy is transported S h = 0.
The remainder of this section is divided into two parts. In the first one, simulations results are compared with available experimental data following a validation process. In the second part, analysis of the simulated flames characteristics are conducted in a qualitative fashion.

Validation of Modeling Strategies
Radial profiles of mean and RMS values of the resolved longitudinal and radial components of droplet velocities are compared with experimental data in Figure 2. As a Cartesian coordinate system is employed the radial direction is represented by the coordinate Y, while the longitudinal by X. Focusing on the first column, one can observe that the proposed strategy is able to recover the mean values of droplets longitudinal velocities reasonably well. Deviations from experimental data are more pronounced downstream of the positions 25 X/D for all cases. Comparing the different simulated cases, further deviations are noticed between positions 5 X/D to 15 X/D. Similar to the observations in [3], it can be clearly noticed that the FSW modeling affects the estimation of mean droplet velocities. Specifically, when β = 0.2, results better agree with experimental data. Considering this aspect and the fact that such deviations are located in a radial region far from the jet centerline (i.e., Y/D > 0.5), it is possible to state that deviations occur in regions where droplets interact with reaction zones. Nevertheless, due to the strong coupling among the diverse modeling approaches it is not straightforward to draw an exact explanation for it. Influences of velocity fluctuations, thickening factor, and mixture stratification are also supposed to contribute to it, as pointed out throughout the discussions presented in this section. Concerning both cases computed with the same value for the model parameter β, marginal deviations are observed for U d between them. Close to the jet core (i.e., Y/D < 0.5) both models deliver the same results. Influences are supposed to be caused by the variations of the flame speed caused by heat losses in the flame front as marginal deviations are noticed when Y/D > 0.5 between positions 5 X/D to 15 X/D.
Considering the second column of Figure 2, it is important to highlight that despite the absence of an SGS turbulent dispersion model, droplet velocity fluctuations agree reasonably well with the experimental data. The discrepancies noticed for the positions 10-15 X/D are in great parts assigned to the flame-dispersed phase interaction, which directly influences the droplet tracking. This aspect is supposed to justify the inflection observed when 0.5 < Y/D < 1.0 in the positions 10 X/D and 15 X/D. By comparing the different cases, it is again noticed that the flame wrinkling model also contributes to deviations.
According to the radial profiles of V d , the experimental values are different from zero at Y/D ≈ 0.0, i.e., the jet centerline. As already discussed in [3], the measurements of this velocity component expose a problem on the experimental configuration. For vertical jets, the time averaged radial velocity along the centerline is always zero. However, this is not observed for all comparison planes, including the injection position. This result possibly indicates that the experimental spray is not strictly vertical. Due to this characteristic, the procedure suggested by De and Kim [6] is used here to define the radial velocity of droplets. Namely, the radial velocities of the carrier phase are used to define the starting conditions for droplets. This option justifies the disparity observed in the boundary condition's plane ( Figure 2) for V d . Altogether, simulations results follow a similar trend as the measured data. The existence of non zero values of V d in the proximity of the jet centerline motivates the addition of a constant offset for all points (for which a zero velocity is attained at the centerline) of the experimental data. Considering them as a reference, a better agreement is noticed for the computed sprays and even better for the case with β = 0.2. However, because of the anomalies perceived in the experimental measurements (see e.g., spray jet anomaly in [45,46]), further discussions about the radial components of mean droplets' velocity are difficult to be conducted. The same is valid for the last column of Figure 2, where the RMS values of the fluctuations of such velocity components are exhibited. These data are rather used here to give support to the discussions about the different simulated cases. An important quantity to observe the mass transfer between phases is the liquid volumetric flux (LVF), whose radial profiles at different longitudinal positions are presented in the first column of Figure 3. Both cases computed with β = 0.5 show lower values for LVF downstream of the position 10 X/D when compared with the case with β = 0.2. This is a straightforward outcome caused by the increased flame surface wrinkling obtained with β = 0.5. As the total FSW increases, the consumption flame speed also increases and the overall flame becomes shorter in the longitudinal direction while broader in the radial direction. Consequently, the spray inner core faces a gas flow with higher temperatures than the case computed with β = 0.2 (see Figure 4). Such higher temperatures increase the evaporation rate and reduce the corresponding LVF. The hot mixture that intensifies the evaporation process is an outcome of the combustion process, which is found in reaction and oxidation zones. As droplets cross a reaction zone, interference with the reaction progress is expected to occur. The amount of droplets interacting with reaction zones has not been quantified in the present work. Nevertheless, studies conducted in [9,16,21] show that FGM tables constructed with premixed flamelets are able to recover the effects of variations of the mixture fraction throughout a flame with good accuracy.The effects of the temperature reduction of the jet inner core and the consequent increase of the simulated LVF are also noticed between both cases computed with β = 0.5. The inclusion of EC effects indeed decreases the mass transfer between phases as observed by the higher values of LVF downstream of the position 15 X/D. Sacomano Filho et al. [16] indicated by means of flames propagating in droplet mists that the evaporative cooling effects also contribute to the reduction of the laminar flame speed. Bearing this aspect in mind, the increase of LVF and decrease of gas temperature (see Figure 4) when including EC are combined features from the reduction of the inner core temperature and the flame speed. However, the contribution given by the inclusion of EC over the LVF is not so pronounced as that given by the reduced value of the model parameter β. It is important to highlight that the case computed with β = 0.2 agrees much better with the experimental data than the other cases. Again, this gives an indication of the relevance of the turbulence-flame interaction modeling.
Considering the droplet size distribution, reasonably good agreement with the experimental data can be perceived for the Sauter mean diameter (SMD-D 32 ) in the inner part of the jet (i.e., Y/D < 0.5) and up to 25 X/D. In the farther downstream regions, the SMD becomes slightly smaller than the measured values in all radial positions. A probable cause for it is the penetration of droplets in higher temperature zones, which also agrees with the observed evolution of LVF radial profiles. The interaction of the dispersed phase with reaction zones may justify the deviations noticed for the SMD for radial positions far from the jet centerline (i.e., Y/D > 0.5) from 5 X/D to 20 X/D.  In the last two columns in Figure 3, radial profiles of the mean and the RMS values of the resolved longitudinal component of the carrier phase velocity are shown. Experimental data acquired for droplets smaller than 10 µm are used as a reference. In fact, both quantities are not strictly comparable, but the experimental data is used here to guide the subsequent discussions. Many aspects observed for the dispersed phase in Figure 2 are noticed in the gas flow. As for the spray, deviations occur close to reaction zones and in a similar fashion among the different cases. The correct prediction of the flame topology is supposed to contribute to these observed deviations. Concerning the velocity fluctuations, it is important to bear in mind that comparisons performed here do not take into account the contribution of SGS quantities since these are not available. Veynante and Knikker [47] showed that this lack of information does not modify mean but RMS values. As a consequence, the inclusion of these SGS contributions should modify the scenario regarding the rms values.  To conclude the validation process, radial profiles of the mean gas temperature are compared with the available experimental data in Figure 4. Results obtained with β = 0.2 better agree with experiments when compared with those from the remaining cases. Effects of the evaporative cooling are noticed in all positions, but these are more pronounced as the distance from the burner nozzle increases. This is expected since the heat removed from the gas phase increases due to the cumulative impact of the evaporation. However, the effects of FSW modeling are more pronounced than the sole consideration of EC while keeping β = 0.5. Specifically, deviations between cases computed with β = 0.2 and β = 0.5 are more pronounced at 10 X/D and 20 X/D. At 30 X/D, similar results can be seen between cases including EC. This overall behavior regarding FSW modeling is explained considering the flame structures presented in the next section. The FSW modeling is strictly connected with the source termω pv (see Equation (3)). According to the qualitative contour plots presented in Figure 8, source terms of the progress variable are predominantly found upstream 25 X/D. Hence, the effects of the FSW modeling are more pronounced upstream of this region. Downstream of it, unburned combustion products are mixed together with fuel vapor released from droplets that are able to cross this main reaction zone. The resulting mixture burns further in a diffusion flame mode, which is less affected by the FSW model. Therefore, different approaches for the computation of β deliver similar results at 30 X/D.
Analyzing the temperature results with more details, simulations computed with β = 0.5 show broader regions with high-temperature values than the results obtained with β = 0.2 at 10 X/D. In fact, simulations computed with β = 0.2 agree quite good with the experimental data at this position. Contributions from the inclusion of EC are seen at the region where the temperature peaks, but these are not high. As previously mentioned, effects of the EC are better seen at 20 X/D and 30 X/D throughout the coordinate Y. At 20 X/D, the ascending part of the temperature profile in Y direction for the case in which β = 0.2 agrees better with experiments than the other cases. The performance of this case at these two positions explain the more accurate prediction of LVF in Figure 3.

Qualitative Analysis of Simulated Cases
To facilitate the interpretation of the mixture composition, Figure 4 presents time-averaged and RMS values of the equivalence ratio (obtained from the resolved fields of mixture fraction as in [48]). The time averaged values of φ point out how strong the FSW model and the inclusion of EC affect the flame structure. Particularly, this variable is connected with the LVF seen in Figure 3, as the vapor released from droplets directly affects the mixture composition. Overall deviations among cases are more pronounced between different definitions of β than the two treatments reserved for the EC. Deviations between both cases computed with β = 0.5 occur at the position 20 X/D. Observing the time-averaged Lagrangian source terms of mass S m and enthalpy S h respectively in Figures 5 and 6, we see that such effect is a direct outcome of the inclusion of EC. Both cases show similar S m profiles at 20 X/D and upstream this position, while clearly different profiles of S h . Therefore, the inclusion of EC effects clearly interferes with the mixture composition, as it is a product of higher transported temperatures of the gas phase (see Figure 4) through the higher droplet evaporation rates observed in terms of LVF profiles in Figure 3. With respect to the FSW modeling, the case computed with β = 0.2 already interferes with the mixture composition at 10 X/D. As seen in Figure 4 and reinforced in Figures 5 and 6, the main contribution to this aspect stems from the region where the reaction zone is found (Figure 8 points out the position of reaction zones by means of the time-averagedω pv ), i.e., 0.5 < Y/D < 1.0. Figures 5 and 6 indicate that small impact of this modeling approach is found in the jet core region, namely Y/D < 0.5. As the flame evolves, more interactions with the spray at the jet inner core occur at a downstream position as seen through the Lagrangian source terms at 20 X/D and 30 X/D. The topology of the flame is directly related to this observation. As seen in Figures 7 and 8 upstream of 20 X/D both flames computed with β = 0.5 depict similar topology. Downstream of this position small flame interactions occur with the spray, whereas droplets that are able to cross the flame tip interact with combustion products coming out of the main flame. The Lagrangian source term profiles indicate that differences between models are also happening at 30 X/D. Nevertheless, these are less intense as no significant impact is noticed in time-averaged profiles of φ in Figure 4 at this position. Considering the RMS of φ, higher values are found in the region that coincides with the descending part of the temperature profiles. Such an outcome can be seen as a combined effect of the interaction of droplets with high-temperature combustion products and radial flow coming out of the main jet core as seen through the radial profiles of V d in Figure 2. Deviations among cases are analogous to those found in the time-averaged profiles of the gas temperature. These are significant between FSW modeling approaches at 20 X/D, whereas between EC treatments at 30 X/D. Regarding the RMS of gas temperature, similar trends are noticed among cases, except at 20 X/D. At this position, higher values are found for the case with β = 0.2 in the proximities of the jet centerline. This characteristic can be explained considering the instantaneous plots of the source terṁ ω pv presented in Figure 7. Different from the cases computed with β = 0.5, the case with β = 0.2 deliver a longer flame in which pockets are frequently formed and released close to the flame tip. This occurs close to 20 X/D, where high fluctuations appear for β = 0.2 in Figure 4. An additional aspect that can be extracted from the contour plots ofω pv presented in Figures 7 and 8 refers to the mixture stratification. Such a feature can be observed in terms of the different intensities oḟ ω pv throughout the flame. The stratification is also noticed in φ radial profiles of Figure 4. The variation of the mixture composition within the computational domain is a direct outcome of the evaporation of liquid droplets. Due to the free jet configuration, the amount of released vapor increases as the spray penetrates into the domain. The opposite behavior is noticed in the liquid phase by means of the LVF profiles in Figure 3.
Contour plots of time-averaged values of the gas temperature give a more complete overview of this quantity for the different cases in Figure 9 than the radial profiles shown in Figure 4. The higher penetration length of the flame computed with β = 0.2 is also observed here considering the penetration of the cold jet core. Additionally, the higher temperatures observed for the case where EC effects are neglected (using the other case computed with β=0.5 as a reference) are also noticed. The extension of the slice plane used in Figure 9 allows the visualization of another import aspect. Considering the jet centerline (i.e., Y/D = 0.0), as the distance from the nozzle increases, the temperature initially grows from the cold-core level up to a first pick which coincides with the main flame tip position. Downstream of it, the temperature slightly reduces and grows up again for positions above 30 X/D. Particularly, at the region located between 20 X/D and 35 X/D for cases computed with β = 0.5, whereas between 24 X/D and 38 X/D for the one computed with β = 0.2, a local cold-core is identified. A similar structure is also perceived, but with much less contrast, in Figure 8 in terms of time-averaged values of the source termω pv . The appearance of it is a result of the combustion of unburnt products mixed with the vapor released from droplets that are able to cross the main reaction zone. In Figure 7 isolated flame structures are found downstream of the main reaction zone. Besides being quite sparse these structures are preferentially located on the jet borders, where the mixture of products and vapor diffuses with the co-flowing air stream. Accordingly, the occurrence of these secondary reactions resemble processes controlled by the flow mixing. This characteristic aligns with the definition of diffusion flames, namely mixing controlled combustion reactions [49]. It is important to highlight that such a multi-regime flame differs from those observed in [10], where the main premixed flame is enveloped by a diffusion-reaction layer. Such an overall behavior is supposed to do not significantly affect the conducted analysis in the present work since the focus has been given to positions upstream of the second cold-core region. Nevertheless, the inclusion of the method recently proposed in [40] to improve the evaporation modeling of droplets interacting with combustion products, as well as the adoption of a multi-regime combustion model (e.g., [50]), may improve the characterization of this post main flame phenomena which may not be reserved to this configuration.

Conclusions
Investigations of the effects caused by the consistent consideration of the evaporative cooling and the turbulence flame interaction demonstrated that both phenomena do interfere not only with the estimation of the flame structure, but also with the prediction of spray properties. As previously pointed out in [3] for a similar flame but with lower Reynolds number, such an outcome reveals the strong interaction among turbulence, flame, and dispersed phase. The necessity of an appropriate computation of each of these phenomena is again reinforced in the present study.
The combination of evaporative cooling effects through the FGM method and the flame surface wrinkling model considering β = 0.2 demonstrated to have the best performance in the validation process. The corresponding simulated case showed the best overall agreement with available experimental data and could recover the main phenomena observed experimentally. For the validation process, the experimentally measured flame EtF5 of the Sydney diluted spray flame burner [31] was used as a benchmark. The other two cases have been used to allow a parametric study of the influence of the EC effects and FSW modeling. The raised hypothesis in that EC effects could be quite important [16] was accessed by comparing two cases, with and without consideration of EC effects at a rather standard of the FSW model parameter β. Although EC effects were shown to be important for the correct flame prediction, the influence of the turbulence-flame interaction modeling demonstrated to be more pronounced.
Thanks to the phenomena explicitness that LES and ATF approaches allow, the main features observed in the simulation results could be recalled by analyzing further computed quantities. Such a clear interpretation of different phenomena could not be easily done if the proposed modeling extensions to account for heat losses were not consistently done.
Places for improvements are found as well as pointed out throughout the result presentation and discussion. It is important to highlight that the proposed modeling strategy is still limited since part of the underlying phenomena to turbulent spray combustion is not accounted for. Indeed, some effects were not comprehensively studied or even understood to allow to infer on their relevance in the simulations. Two of them are the inclusion of actual mixture composition on droplet evaporation (see e.g., [40]) and effects of differential diffusion in spray flames (see e.g., [16]). Additional effects that may improve the predicting ability of the proposed modeling strategy is, on the one hand, the inclusion of models capable to represent multi-regime combustion, as studied by Knudsen et al. [13] and Hu and Kurose [50], and the incorporation of single droplet burning effects (see [51]) on the other. Nevertheless, it worths to be noticed that the proposed methods following the inclusion of the evaporative cooling effects in the chemistry and the selected FSW modeling converge to an optimum combination of strategies which better approaches the experimental observations. Models able to represent the above-listed phenomena are desired and should improve the capability of the presented module.