Large Eddy Simulation of a Turbulent Spray Jet Flame Using Filtered Tabulated Chemistry

+is work presents Large Eddy Simulations of the unconfined CORIA Rouen Spray Burner, fed with liquid n-heptane and air. Turbulent combustionmodeling is based on the Filtered TAbulated Chemistrymodel for LES (F-TACLES) formalism, designed to capture the propagation speed of turbulent stratified flames. Initially dedicated to gaseous combustion, the filtered flamelet model is challenged for the first time in a turbulent spray flame configuration. Two meshes are employed. +e finest grid, where both flame thickness and wrinkling are resolved, aims to challenge the chemistry tabulation procedure. At the opposite the coarse mesh does not allow full resolution of the flame thickness and exhibits significant unresolved contributions of subgrid scale flame wrinkling. Both LES solutions are extensively compared against experimental data. For both nonreacting and reacting conditions, the flow and spray aerodynamical properties are well captured by the two simulations. More interesting, the LES predicts accurately the flame lift-off height for both fine and coarse grid conditions. It confirms that the modeling methodology is able to capture the filtered turbulent flame propagation speed in a two-phase flow environment and within grid conditions representative of practical applications. Differences, observed for the droplet temperature, seem related to the evaporation model assumptions.


Introduction
Aeronautical engines are operated with liquid fuel directly injected in the combustor. Two-phase combustion is extremely difficult to understand as it requires a simultaneous access to a large number of highly-correlated thermophysical properties [1]. e Large Eddy Simulation (LES) approach, which represents, nowadays, the best compromise between cost and accuracy to simulate complex reactive flows, is especially attractive for computing realistic gas turbine combustors [2,3]. Despite recent impressive progress, many efforts are still being put by the combustion modeling community to develop and validate LES for turbulent spray flame computational strategies [4][5][6][7][8]. Model comparison against accurate experimental data is crucial to properly assess the ability of numerical strategies to recover the turbulent spray flame properties. It includes the flow velocity, the droplets characteristics and the flame structure.
Flame stabilization and pollutant formation require a fine description of the interactions between combustion kinetics and turbulence [9].
is is especially true in two phase combustors, where fuel-air mixing and finiterate kinetics phenomena must be carefully modeled at the subgrid scale to capture the stabilization physical process [10]. Tabulated chemistry methodologies have been developed during the last decades to account for detailed chemistry effects at a reduced CPU cost [11,12]. Among them, the Filtered Tabulated Chemistry for LES (F-TACLES), has been especially developed to incorporate complex chemistry effects in an LES formalism [13]. It consists in tabulating the chemical ingredients needed by the LES in a filtered lookup table. F-TACLES has been applied to complex gaseous turbulent flames such as stratified [14] and nonadiabatic [15,16] configurations. However, constrained by severe assumption of a low-dimensional manifold reduction, the suitability of such LES-flamelet approach for two-phase reactive flows remains to be demonstrated [17]. e suitability of F-TACLES to turbulent spray flames simulations, which has never been addressed, is the main objective of this article. e present work presents the first application of the filtered tabulated chemistry model F-TACLES in a turbulent spray combustion configuration. High-fidelity databases have to be specifically designed to provide technical performance metrics for model LES validation [18]. e configuration retained here is a new well-instrumented experimental turbulent spray flame that has been designed and operated at the CORIA laboratory [19]. Simulations are conducted on two different grids: a coarse one, representative of meshing constrains encountered in industrial applications, and a fine one for which the size of the cells within the reaction zones has been chosen so that both flame thickness and subgrid flame wrinkling are fully resolved. e fine grid simulation will challenge the ability of the chemistry tabulation to retrieve the spray flame structure [20], whereas the coarse LES will also test the suitability of F-TACLES to capture unresolved interactions between the spray flame and turbulence. Experimental and numerical data are compared and analyzed in terms of gas velocity, spray diameter distribution and velocity, flame structure, and spray temperature.

N-Heptane Air Combustion Chemistry.
Liquid n-heptane is used in the targeted experimental configuration. ree nheptane/air chemical schemes are considered: the detailed chemical mechanism POLIMI [21] which includes 106 species and 1738 reactions, the two-step global scheme 2S [6] fitted by using the methodology proposed in [22], and an Analytically Reduced Scheme ARC developed in [23] by applying methodology from [24] which includes 24 transported species, 32 species in quasi-steady state and 217 reactions. e ability of the three chemical schemes to reproduce experimental laminar flame burning velocity measurements [25] is shown in Figure 1. e global step chemistry fails to reproduce the flame speed over rich conditions and is, therefore, not retained in this study. Both POLIMI and ARC scheme fairly capture the experimental measurements, but the number of species and the stiffness of the schemes remain too important for a direct coupling with an LES flow solver.
A tabulated chemistry method is retained to drastically reduce the CPU cost of the chemistry model [11]. e chemical subspace accessed by a spray flame is here approximated by an ensemble 1-D premixed flamelet trajectories, following FPI [26] or FGM [12] approaches. Each thermochemical variable φ expresses in terms of a progress variable Y c and a mixture fraction Z as follows: where TAB indicates that the variable φ is stored in a look-up table. e progress variable is defined as Y c � N sp k�1 n k Y k , where n k is the weight associated to species mass fraction Y k . φ may include chemical reaction rates, species mass fractions, density but also thermodynamical and transport properties such as the heat capacity c p and thermal conductivity λ. e suitability of tabulated chemistry to twophase reactive flows has been investigated by Franzelli et al. [7,20]. FPI tabulated chemistry actually reproduces accurately the temperature and heat release profiles over a wide range of spray conditions. e chemical table is built from a library of laminar freely propagating n-heptane/air premixed flamelet computed with the REGATH code [27] and by using the POLIMI detailed mechanism reactions ( [21]).

Turbulent Combustion Modeling.
e premixed flamelet manifold is coupled to LES using the F-TACLES formalism, developed first for premixed combustion [13] and then extended to stratified flames [14]. e modification proposed by Mercier et al. [15] to account for the impact of differential diffusion on the flame consumption speed is retained. e F-TACLES model assumes that the chemical structure of the filtered flame front is captured by an ensemble of 1-D filtered flame elements. e premixed flamelets used to build-up the FPI manifold are here filtered in pysical space at a size Δ. Filtered thermochemical variables φ are, therefore, stored in terms of Y c , z, and Δ in a chemical look-up table such as where Y c and z are the filtered progress variable and mixture fraction, respectively. e filter size Δ is chosen to broaden the flame so that the filtered reactive layer is resolved on the LES grid. As demonstrated in [13], about 4-5 nodes are needed to ensure a proper filtered flame front propagation without introducing numerical artifacts. e flow is given by the solution of the filtered Navier-Stokes equations. As low Mach number flow assumption is made in this work, the filtered temperature T and the density ρ are tabulated in the filtered chemical lookup table given by equation (2), as any other thermochemical variables φ [28]. z and Y c are solutions of the two following balance equations: where ρ is the density, μ t is the turbulent viscosity, Sc t is the turbulent Schmidt number, _ ω evap is the source term of mixture fraction due to the evaporation of the spray, Ξ Δ is the subgrid scale flame wrinkling, α Y c is the progress variable diffusion factor, ρ 0 is the density in fresh gases, D 0 is the diffusion coefficient in fresh gases, Ω Y c is the progress variable unresolved convective fluxes due to thermal expansion, and _ ω Y c is the progress variable reaction rate. e functions α Y c , Ω Y c , and _ ω Y c in equation (4) are designed to model the subgrid scale (SGS) laminar contributions to molecular diffusion, convection, and chemical reaction, respectively. ese terms are computed from 1-D filtered premixed flamelet solutions and stored in the F-TACLES look-up table as follows: where * denotes quantities issued from the computations of 1-D unstrained laminar premixed flames. By construction, this model propagates the resolved flame front at the subgrid scale turbulent flame speed S T,Δ : where S 0 l is the adiabatic consumption speed of a freely propagating laminar premixed flame. e model for Ξ Δ is modeled using the Charlette et al. formulation [29]: where Re Δ � (u Δ ′ Δ)/] and u Δ ′ are the subgrid scale Reynolds number and velocity fluctuations, respectively, while δ 0 l is the laminar flame thickness. e efficiency function Γ Δ ( [29]) estimates the net straining effect of all turbulent scales smaller than Δ. e exponent β is set constant and equal to β � 0.5 as initially prescribed in [29].

LES Equations for Two-Phase Flow.
e two phase flow is modeled by a Euler-Lagrange approach. Filtered governing equations for continuity, momentum and energy are solved in addition to balance equations for filtered progress variable and mixture fraction given by equations (3) and (4), respectively. e diluted spray is described with a Lagrangian point-force approach, which is two-way coupled to the gaseous phase. e following transport equations are solved for each droplet: where x p is the particle position, u p is the particle velocity, u is the gas velocity, m p is the particle mass, C D is the drag coefficient, ] is the kinematic viscosity, ρ p is the particle density, and Re p is the particle Reynolds number. e evaporation of the spray is modelled with the classical approach derived by Spalding [30]. e droplet mass transfer equation reads where d p is the particle diameter, D is the diffusion coefficient, Sh is the Sherwood number, and B M is the Spalding mass number. e term _ ω evap in the mixture fraction equation is obtained by adding the mass transfer contribution of all the droplets around each node of the mesh: where V node is the volume around the node. e other droplet parameters are derived by integrating either the droplet mass or energy equations. Droplet temperature T p and diameter d p are obtained by solving the following set of equations: where τ p is the thermal characteristic time of the Spalding model, T ∞ is the gas temperature in the far field, L v is the latent heat of vaporization of the fuel, B T is the Spalding thermal number, C p,1/3 is the heat capacity at a classical reference state assuming a one-third/two-third equilibrium Journal of Combustion between the far field and the droplet surface, μ 1/3 is the dynamic viscosity at the same reference state, and S c is the Schmidt number.

Experimental Configuration
e experimental configuration is an n-heptane spray/air jet burner experimented at CORIA by [19]. It is operated at atmospheric pressure and 298 K. e air injection is performed from a plenum to a non-swirling injector in order to generate the coflow where the liquid fuel is atomized. e air mass flow rate of is 6 g·s −1 . e injection of liquid n-heptane comes from a simplex injector that generates a hollow cone with a mass flow rate of 0.28 g·s −1 . A general view of the configuration geometry is shown in Figure 2.
Several experimental measurements have been performed. e Phase Doppler Anemometry (PDA) gives access to the gas and spray velocity and the spray diameter distribution. e flame structure is determined thanks to OH Planar Laser-Induced Fluorescence (PLIF). Finally, the Global Rainbow Technique (GRT) ( [31]) provides the spray temperature, which is rarely available in experimental diagnostics. Further details about these measurements can be found in [19]. e flame structure shown in Figure 3 by the OH-PLIF measurement exhibits a double branch. e inner flame front corresponds to a premixed flame where the small droplets are vaporized rapidly and the high levels of turbulence favor the air/fuel mixing, forming a highly wrinkled flame front. e outer flame front is closer to a diffusion flame, where air located outside reacts with rich hot gases still containing a large amount of unburnt gaseous n-heptane. OH-PLIF also shows that the flame is lifted from the injection plane.

Numerical Setup
is experiment has been previously studied numerically by Shum-Kivan et al. [6] by using global two-step chemistry [32] combined with the TFLES approach [33,34]. e flow velocity, as well as the droplet size distribution and velocity have been well predicted. However, an underestimation of the flame lift-off has been observed, which is probably due to the limitation of the reduced two-step chemistry model. Other approaches were tested on this configuration, for example, with the stochastic field method [35]. e computational domain defined in [6] is also used in the present study. Two cases (A and B) are considered. Case A features an unstructured mesh composed of 53 million elements and 10.5 million nodes, identical to [6]. Case B is performed on a coarser mesh of 17 million elements and 3.5 million nodes. Case A is sufficiently resolved so that artificial broadening of the flame front is not required. Indeed, the mesh size in the reaction zone is less than 0.1 mm, whereas the minimum possible flame thickness, given by a laminar stoichiometric premixed freely propagating flames, is about 0.5 mm. With 5 nodes across the flame front, the resolution of the chemical layer is therefore sufficient to ensure the proper propagation of the flame without introducing numerical artifact in both premixed [13] and stratified [14] mixtures. e flame front resolution in Case B is more representative of LES conditions encountered in industrial configurations. e mesh size in the reaction zone, around 0.5 mm, is not sufficient to resolve the flame front. e filter size Δ associated to the flame is therefore chosen to artificially enlarge the filtered reactive layer front is therefore required. In addition, the subgrid scale flame wrinkling cannot be neglected and requires modeling. e modeling challenges are to recover the flame dynamic on case B, where the subgrid scale turbulent combustion model is of importance. e chemical table is built from a library of laminar freely propagating n-heptane/air premixed flamelet computed with the REGATH code [27] and by using the POLIMI 106 detailed mechanism made of 106 species and 1738 reactions [21]. For case A simulation, as the flame is fully resolved on the LES mesh, this look-up table is directly used to close equation (4) without being filtered (Δ � 0). Consequently, by assuming the flamelet regime, the flame wrinkling is also fully resolved on the LES grid and Ξ Δ � 1. At the opposite, the flamelet library is filtered in Case B by using a filter width Δ � 3.5 mm so that the resolved filtered flame thickness is sufficient to capture the flame consumption speed on the coarse mesh. Subgrid scale flame wrinkling is modelled as in Charlette et al. [29] given by equation (7). e combustion model properties used for case A and B are summarized in Table 1.
e YALES2 flow solver is used [36]. e time integration relies on a low-Mach number projection method for variable density flows. e temporal integration and spatial discretization are of fourth order. e subgrid scale Reynolds stresses are closed with the SIGMA model [37]. e injected spray is polydispersed in size, following a two-parameter Rosin-Rammler distribution [38]   obtained with the Liquid Injection for Swirl Atomizers (LISA) formalism [39] to obtain the desired swirled hollow cone spray. Parameters of droplet distribution in size are empirically adjusted to fit measurements at 10 mm above the burner exit as shown in

Results and Analysis
e two cases A and B are computed in both nonreactive and reactive configurations. erefore, four simulations are presented in the following sections. e nonreacting cases are appended with the suffix -NR and the reacting ones with -R. Figure 5 shows the positions of the profiles that are used for comparing experimental and numerical results. e temperature field of case A-R is shown in transparency to indicate the position of the flame in reacting cases. Figures 6 and 7 show for cases A-NR and B-NR, the instantaneous and mean axial velocity fields in the central vertical plane. e mean flow topology is very similar for both meshes. Several zones are identified. First, the flow is accelerated up to 30 m/s downstream the injector before exiting into the atmosphere. A recirculation zone also appears at the exit of the injector, where the liquid injection is made. e effect of the injection of the droplets is visible in this zone, with a local increase of the axial velocity. Finally, a mixing layer appears between the      fast air that exits the injector and the air at rest in the atmosphere. Figures 8 and 9 show for cases A-R and B-R, the instantaneous and mean axial velocity fields in the central vertical plane. e general flow topology is similar to nonreacting cases. e main difference is linked to the presence of the flame, which enlarges the width of the jet through thermal expansion. Grid effects are visible in the instantaneous axial velocity field of radial velocity shown in Figures 6 and 8, where smaller vortices appear for the fine mesh.

Flow Topology and Gas Velocity.
Axial velocity component from LES results are compared against the measurements at 10, 20, and 40 mm high above the burner exit. Nonreacting (left) and reacting (right) mean and RMS quantities are plotted in Figures 10 and 11, respectively. e solutions of both cases A and B agree with the experimental data, meaning that the flow statistics are well captured, even on the coarse grid. e main difference is the underestimation of the maximal axial velocity around 10 mm above the injection plane. e origin of this discrepancy can be attributed to the resolution of the boundary layer in the injector. Indeed, a wall-law approach is chosen for this configuration, and the boundary layer velocity profile is not fully resolved. A finer mesh close to the injector walls would improve the prediction of the peak of velocity. e effect of the thermal expansion from the flame is visible on the profiles at Z � 40 mm. e axial velocity in nonreacting conditions drops rapidly to 0 m/s between r � 5 mm and r � 20 mm while in reacting conditions, the axial velocity decreases slowly between r � 5 mm and r � 20 mm. e RMS is correctly captured for both meshes. e effect of the mixing layers (between the recirculation zone and the main flow and between the main flow and the air at rest) is visible as the two peaks of axial velocity RMS at Z � 10 mm.
e results for the radial velocity component are plotted in Figure 12 (mean) and 13 (RMS). e simulations capture the general features of the mean radial velocity profiles. e RMS are also rather well predicted by the simulations. At Z � 10 mm, the peaks at r � 0 and 10 are correctly located in the simulations. e general flow topology prediction by the simulations is satisfactory, for both meshes and for both nonreacting and reacting conditions. Figure 14(a) shows an instantaneous normalized OH mass fraction field for each simulated case and an instantaneous snapshot of OH-PLIF measurements. It gives a qualitative analysis of the instantaneous flame structure, which is challenging to compute as the stabilization processes are very sensitive to finite-rate chemistry effects. e inner flame front, by a highly wrinkled by the turbulence, is qualitatively reproduced by the LES. e outer diffusion flame observed in the experiments, featuring a large and unwrinkled reaction zone, is also present. e mean normalized OH mass fraction field for cases A and B are compared against the mean shot of OH-PLIF measurements in Figure 14(a) . e mean inner flame front position is well captured by the simulations and is located at |π| ≈ 15 mm up to z � 80 mm. e mean OH-PLIF measurements show that the outer flame front extends up to |x| � 40 mm.

Flame Topology.
is comparison shows that even if the instantaneous flame structure seems qualitatively well retrieved by the simulations, the mean outer flame front position is not perfectly captured by the simulations. Indeed, both simulations on coarse and fine grids predict that the outer flame front extends up to |x| ≈ 30 mm mm and quickly merges with the inner flame front for z > 50 mm. e inner flame front is located in a region of high velocity while the outer one is located in a low-velocity region, as shown in Figure 10. erefore, the amount of flow-through times simulated differs between the two flame fronts. e statistics are well converged for the inner flame front because the velocity is much higher. On the contrary, as the velocity in the outer flame front is low, the simulated physical time (tens of milliseconds) may not be sufficient to capture the dynamics of the outer flame front that was found in the experiments, where the OH-PLIF shots are averaged over a much longer period of time (several seconds). e lift-off of the flame is a critical aspect of this flame. In order to assess the lift-off height in the simulations, Figures 15 and 16 show a contour of temperature in transparency for both meshes. ese views demonstrate that the lift-off height is fairly constant for both meshes. Figures 17 and 18 show a clip in the central vertical plane of the contour of temperature presented above. e influence of the mesh is visible in Figure 17 where the flame wrinkling is more resolved in case A (fine mesh) than in case B (coarse mesh). e lift-off height is defined experimentally as the closer position of the flame front from the burner exit. e flame front position is defined from the maximum value isoline given by the mean OH-PLIF signal shown in Figure 14. e lift-off of the flame is estimated similarly from the simulations. is height depends on the angular position since the flame is not perfectly axisymmetric. e circumferential mean and RMS of the lift-off position are therefore computed. e experimental value is 25 ± 3 mm while case A recovers a lift-off of 22 ± 1 mm and case B a lift-off of 24 ± 1 mm. Comparison between case A and B shows that the F-TACLES approach is able to model fairly well unresolved flame turbulence interaction on a coarse mesh representative of practical industrial conditions.
Previously published computations with a global twostep mechanism ( [6]) underpredict the flame lift-off h lo by approximately 20%. Surprisingly, simulations conducted with a reduced analytical scheme involving 24 transported species, 32 quasisteady state species, and 217 reactions also  Journal of Combustion did not succeed to retrieve the flame lift-off, with a CPU cost 10 times higher ( [23]) than F-TACLES. Note that the flame is rather controlled by front propagation than autoignition for two reasons. First, there is no hot stream which could increase sufficiently the fresh gas temperature to reach self-ignition conditions. Second, results obtained in [6,23] with an analytically reduced scheme including 56 species do not evidenced the presence of radical species characteristics of autoignition downstream the flame base. Such a configuration is favourable for an F-TACLES model which has been designed to capture flame propagations with or without subgrid scale wrinkling.
Indeed, with the F-TACLES tabulated chemistry method, the flame lift-off height is recovered for both  Journal of Combustion 13 statistics, normalized by the global scheme computation, is also indicated. Figures 19 and 20 show the particles in the central vertical plane colored by their diameter for the cold and reacting cases, respectively. For both conditions, the distribution of diameter is similar. e smaller droplets are located in the central part of the flow while the larger droplets are located on the outer part of the spray. e influence of the flame in Figure 20 is the low density of particles above z � 20 mm, especially on the outer region. Figure 21 compares at 10, 20, and 40 mm high above the burner exit, the mean spray diameter as a function of the radial coordinates for the cold and reacting cases, respectively. e LES results show a correct evolution of the radial stratification in droplet diameter for both cases A and B. e small droplets follow the streamlines because of their small Stokes number and are therefore located at the center of the flow. e larger droplets, characterized by a higher Stokes number, follow a ballistic trajectory and are located on the outer rim of the spray, as a result of the hollow cone injection. e profiles are similar in both reacting and nonreacting cases between 0 and 20 mm, as flame is located  14 Journal of Combustion   Journal of Combustion Despite a significant computed flow-through time (equal to 3 and 5 for cases A and B, respectively), a lack of statistics is observed in Figure 21 at high radial values for both reactive and nonreactive cases. It causes discrepancies between numerical and experimental solutions, which are attributed to the number of large droplets in the outer part of the jet being too small to ensure the statistical convergence of the Lagrangian phase. Figures 22 and 23 show the particles in the central vertical plane colored by their axial velocity for the cold and reacting cases, respectively. In Figure 22, the small droplets reach high axial velocity (up to 30 m/s), carried by the surrounding gas while the large droplets velocity decreases because of drag. In Figure 23, the droplets have the same behavior. Some large droplets are not entering the flame and are not consumed at the extremity of the spray.

Spray Velocity.
Droplet axial velocity is reported in Figure 24 for the cold and reacting cases, respectively. e experimental measurements are colored by the diameter of the spray at the considered radial position. Green squares correspond to particle diameters lower than 15 microns, blue squares to diameters between 15 and 35 microns, and red squares to diameters larger than 35 microns. e agreement is good for small-to-medium droplets (below 35 microns), but both LES cases predicts a higher velocity than the experiments for the large droplets. is discrepancy is attributed to the method of injection (from [39]) that may overestimate the large droplets velocity.
Droplet radial velocity is reported in Figure 25 for the cold and reacting cases, respectively. As for the axial velocity, the velocity of the small droplets is well predicted by all the simulations and the velocity of the large droplets is overestimated. Figures 26 and 27 show the particles in the central vertical plane colored by their temperature for the cold and reacting cases, respectively. e scale is 280 K < T p < 300 K for the cold case and 280 K < T p < 370 K for the reacting cases. In Figure 26, the small droplets temperature decreases rapidly to ≈ 280 K as they are convected downstream. is evolution is due to the evaporation. e same process exists for the larger droplets, but much slower. In the reacting case, below the flame, the behavior is the same as in the cold case. When the droplets enter the flame, the ones that are not entirely evaporated are heated rapidly to ≈ 370 K because of the heat released by the flame. e small droplets located in the center of the flow are progressively heated by the hot gases until they are fully evaporated. e droplet temperature predicted by the LES is now compared with the Global Rainbow Technique (GRT) measurements, whose uncertainty is ± 3 K [19]. Figure 28 presents radial profiles of temperature for the cold (left) and reacting (right) configurations. e experimental data highlight two zones. For r > 5 mm, the droplets reach quickly the wet bulb temperature, from the first measured radial profiles, i.e., 20 mm above the burner exit, whereas the liquid spray remains at the injection temperature around the centerline. is trend is not captured by the simulation, which predicts the wet bulb temperature for all droplet positions. e wet bulb temperature is defined as the equilibrium temperature reached by evaporating a liquid to saturation in a gas. is difference between simulations and experiments could be explained by limitations of the evaporation model ( [40]

Spray Temperature.
where ρ p is the droplet density; d p , its diameter; Sc, the Schmidt number; Sh, the Sherwood number; C p,k , the heat capacity at constant pressure of the n-heptane; C p,1/3 and μ 1/3 , the heat capacity at constant pressure and the dynamic viscosity of the mixture according to the 1/3-2/3 rule (see Chapter 1); B T , the thermal Spalding number; and B M , the mass Spalding number. As τ th is proportional to the square of the droplet diameter, temperature will evolve slower for the larger droplets than for the smaller. Figure 29 presents axial profiles of temperature for the cold (left) and reacting (right) configurations. For r � 0 mm, the droplets (which are small at this radial position) temperature drops quickly to ≈ 282 K. As the radial distance r increase, the mean droplet diameter growth as discussed previously, and the droplets temperature decreases.
is tendency is consistent with the Spalding model assumptions.
Another possible explanation would be the choice of the injection model, which, by injecting all droplets from the same point, does not reproduce the spatial distribution of droplets induced by the liquid sheet break-up. Despite a correct prediction of the overall particle size, a local misprediction of the droplet distribution would also impact the mean liquid temperature. A way to overcome this difficulty would be to inject the droplets further downstream and not at the real position of injection. In reacting conditions, in the burnt gases region, located at r > 10 mm and z > 20 mm, the droplet temperature rises quickly due to the high gas temperature. is phenomenon observed in the experiments is fairly tackled by the simulations. However, the droplet temperature measured downstream, between the inner and the outer branch of the flame, reaches a thermal equilibrium around 331 K whereas the numerical simulation predicts 367 K, which is close to the boiling temperature of n-heptane. As discussed in [41], this discrepancy may be also attributed to the Spalding evaporation model, where the limiting value is the boiling temperature. A comparison between the Spalding and Abramzon-Sirignano models, proposed in [42], highlights the differences in droplet temperature predictions.

Conclusion
e first simulation with the F-TACLES formalism in a spray combustion configuration has been performed. e results show good agreement on the spray diameter and velocity, gas velocity, flame structure and lift-off with respect to experimental data. e complex flame structure, which presents an inner premixed flame front and an outer diffusion branch, is well reproduced by the simulation, even on the coarse grid representative of meshing conditions encountered in industrial applications. Fine grid simulations showed that tabulated chemistry based on premixed flamelets is adequate to capture the spray flame chemistry. e good prediction obtained on the coarse grid also demonstrates the ability of F-TACLES to model the unresolved interactions between the spray flame and turbulence. In particular the flame stabilization process is well captured by the turbulent combustion model. As the supplementary CPU cost induced by the combustion model is very low, this method is of interest for the gas turbine engineering community. However, another issue remains to be addressed. Significant discrepancies are indeed found for the droplet temperature. e influence of the droplet evaporation model and of the liquid sheet atomization on the spray temperature should be investigated in the future.

(Nomenclature entries should have the units identified)
A: Reynolds filter of variable Data Availability e numerical data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.