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

Authors: Adrien Chatelier; Benoît Fiorina; Vincent Moureau; Nicolas Bertier
Journal: Journal of Combustion
 (2020)

---

## Abstract

This work presents Large Eddy Simulations of the unconfined CORIA Rouen Spray Burner, fed with liquidn-heptane and air. Turbulent combustion modeling is based on the Filtered TAbulated Chemistry model 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. The 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.


---

## Body

## 1. 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 thermo-physical properties [1]. The 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–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]. This is especially true in two phase combustors, where fuel-air mixing and finite-rate 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]. The suitability of F-TACLES to turbulent spray flames simulations, which has never been addressed, is the main objective of this article.The 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]. The 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. The 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.


## 2. Turbulent Spray Combustion Modeling

### 2.1.N-Heptane Air Combustion Chemistry

Liquidn-heptane is used in the targeted experimental configuration. Three n-heptane/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. The ability of the three chemical schemes to reproduce experimental laminar flame burning velocity measurements [25] is shown in Figure 1. The 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.Figure 1
Comparison of laminar flame speed between experimental data from [25] (green stars) and numerical simulations (lines).A tabulated chemistry method is retained to drastically reduce the CPU cost of the chemistry model [11]. The 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 Yc and a mixture fraction Z as follows:(1)φ=φTABYc,z,where TAB indicates that the variable φ is stored in a look-up table. The progress variable is defined as Yc=∑k=1NspnkYk, where nk is the weight associated to species mass fraction Yk. φ may include chemical reaction rates, species mass fractions, density but also thermodynamical and transport properties such as the heat capacity cp and thermal conductivity λ. The suitability of tabulated chemistry to two-phase 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. The 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]).


### 2.2. Turbulent Combustion Modeling

The 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]. The modification proposed by Mercier et al. [15] to account for the impact of differential diffusion on the flame consumption speed is retained. The F-TACLES model assumes that the chemical structure of the filtered flame front is captured by an ensemble of 1-D filtered flame elements. The 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(2)φ˜=φ˜TABY˜c,z˜,Δ,where Y˜c and z˜ are the filtered progress variable and mixture fraction, respectively. The 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.The 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 temperatureT˜ and the density ρ¯ are tabulated in the filtered chemical look-up table given by equation (2), as any other thermochemical variables φ˜ [28]. z˜ and Y˜c are solutions of the two following balance equations:(3)∂ρ¯z˜∂t+∂∂xiρ¯u˜iz˜=∂∂xiλ¯cp+μtSct∂z˜∂xi+ω˙˜evap,(4)∂ρ¯Y˜c∂t+∂∂xiρ¯u˜iY˜c=∂∂xiΞΔαYcY˜c,z˜ρ0D0∂Y˜c∂xi+ΞΔΩYcY˜c,z˜+ΞΔρ¯ω˙˜YcY˜c,z˜,where ρ is the density, μt is the turbulent viscosity, Sct 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, αYc is the progress variable diffusion factor, ρ0 is the density in fresh gases, D0 is the diffusion coefficient in fresh gases, ΩYc is the progress variable unresolved convective fluxes due to thermal expansion, and ω˙˜Yc is the progress variable reaction rate.The functionsαYc, ΩYc, and ω˙˜Yc in equation (4) are designed to model the subgrid scale (SGS) laminar contributions to molecular diffusion, convection, and chemical reaction, respectively. These terms are computed from 1-D filtered premixed flamelet solutions and stored in the F-TACLES look-up table as follows:(5)αYcY˜c,z˜,Δ=−∑k=1Nspnkρ∗Yk∗Vk∗¯ρ0D0∂Y˜c∗/∂x∗,ΩYcY˜c,z˜,Δ=ρ0∗Sl∗∂Yc∗∂x∗¯−ρ0∗Sl∗¯∂Y˜c∗∂x∗,ω˙˜YcY˜c,z˜,Δ=ω˙˜Yc∗,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 speedST,Δ:(6)ST,Δ=ΞΔSl0,where Sl0 is the adiabatic consumption speed of a freely propagating laminar premixed flame. The model for ΞΔ is modeled using the Charlette et al. formulation [29]:(7)ΞΔ=1+minmax0,Δδl0−1,ΓΔΔδl0,uΔ′Sl0,ReΔuΔ′Sl0,where ReΔ=uΔ′Δ/ν and uΔ′ are the subgrid scale Reynolds number and velocity fluctuations, respectively, while δl0 is the laminar flame thickness. The efficiency function ΓΔ ([29]) estimates the net straining effect of all turbulent scales smaller than Δ. The exponent β is set constant and equal to β=0.5 as initially prescribed in [29].


### 2.3. LES Equations for Two-Phase Flow

The 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. The diluted spray is described with a Lagrangian point-force approach, which is two-way coupled to the gaseous phase. The following transport equations are solved for each droplet:(8)dxpdt=up,mpdupdt=mpup−u3CDRepρν4ρpdp2withRep=dpup−uν,where xp is the particle position, up is the particle velocity, u is the gas velocity, mp is the particle mass, CD is the drag coefficient, ν is the kinematic viscosity, ρp is the particle density, and Rep is the particle Reynolds number.The evaporation of the spray is modelled with the classical approach derived by Spalding [30]. The droplet mass transfer equation reads(9)m˙p=−πdpρDShlog1+BM,where dp is the particle diameter, D is the diffusion coefficient, Sh is the Sherwood number, and BM is the Spalding mass number. The 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:(10)ω˙˜evap=−1Vnode∑droplet∈nodem˙p,where Vnode is the volume around the node. The other droplet parameters are derived by integrating either the droplet mass or energy equations. Droplet temperature Tp and diameter dp are obtained by solving the following set of equations:(11)dTpdt=−1τpTp−T∞−LvBTCp,1/3,ddp2dt=−2Shμ1/3log1+BMdpρpSc,τp=ρpdp26ScSh⋅μ1/3Cp,kCp,1/3BTlog1+BM,where τp is the thermal characteristic time of the Spalding model, T∞ is the gas temperature in the far field, Lv is the latent heat of vaporization of the fuel, BT is the Spalding thermal number, Cp,1/3 is the heat capacity at a classical reference state assuming a one-third/two-third equilibrium between the far field and the droplet surface, μ1/3 is the dynamic viscosity at the same reference state, and Sc is the Schmidt number.


## 2.1.N-Heptane Air Combustion Chemistry

Liquidn-heptane is used in the targeted experimental configuration. Three n-heptane/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. The ability of the three chemical schemes to reproduce experimental laminar flame burning velocity measurements [25] is shown in Figure 1. The 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.Figure 1
Comparison of laminar flame speed between experimental data from [25] (green stars) and numerical simulations (lines).A tabulated chemistry method is retained to drastically reduce the CPU cost of the chemistry model [11]. The 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 Yc and a mixture fraction Z as follows:(1)φ=φTABYc,z,where TAB indicates that the variable φ is stored in a look-up table. The progress variable is defined as Yc=∑k=1NspnkYk, where nk is the weight associated to species mass fraction Yk. φ may include chemical reaction rates, species mass fractions, density but also thermodynamical and transport properties such as the heat capacity cp and thermal conductivity λ. The suitability of tabulated chemistry to two-phase 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. The 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]).


## 2.2. Turbulent Combustion Modeling

The 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]. The modification proposed by Mercier et al. [15] to account for the impact of differential diffusion on the flame consumption speed is retained. The F-TACLES model assumes that the chemical structure of the filtered flame front is captured by an ensemble of 1-D filtered flame elements. The 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(2)φ˜=φ˜TABY˜c,z˜,Δ,where Y˜c and z˜ are the filtered progress variable and mixture fraction, respectively. The 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.The 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 temperatureT˜ and the density ρ¯ are tabulated in the filtered chemical look-up table given by equation (2), as any other thermochemical variables φ˜ [28]. z˜ and Y˜c are solutions of the two following balance equations:(3)∂ρ¯z˜∂t+∂∂xiρ¯u˜iz˜=∂∂xiλ¯cp+μtSct∂z˜∂xi+ω˙˜evap,(4)∂ρ¯Y˜c∂t+∂∂xiρ¯u˜iY˜c=∂∂xiΞΔαYcY˜c,z˜ρ0D0∂Y˜c∂xi+ΞΔΩYcY˜c,z˜+ΞΔρ¯ω˙˜YcY˜c,z˜,where ρ is the density, μt is the turbulent viscosity, Sct 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, αYc is the progress variable diffusion factor, ρ0 is the density in fresh gases, D0 is the diffusion coefficient in fresh gases, ΩYc is the progress variable unresolved convective fluxes due to thermal expansion, and ω˙˜Yc is the progress variable reaction rate.The functionsαYc, ΩYc, and ω˙˜Yc in equation (4) are designed to model the subgrid scale (SGS) laminar contributions to molecular diffusion, convection, and chemical reaction, respectively. These terms are computed from 1-D filtered premixed flamelet solutions and stored in the F-TACLES look-up table as follows:(5)αYcY˜c,z˜,Δ=−∑k=1Nspnkρ∗Yk∗Vk∗¯ρ0D0∂Y˜c∗/∂x∗,ΩYcY˜c,z˜,Δ=ρ0∗Sl∗∂Yc∗∂x∗¯−ρ0∗Sl∗¯∂Y˜c∗∂x∗,ω˙˜YcY˜c,z˜,Δ=ω˙˜Yc∗,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 speedST,Δ:(6)ST,Δ=ΞΔSl0,where Sl0 is the adiabatic consumption speed of a freely propagating laminar premixed flame. The model for ΞΔ is modeled using the Charlette et al. formulation [29]:(7)ΞΔ=1+minmax0,Δδl0−1,ΓΔΔδl0,uΔ′Sl0,ReΔuΔ′Sl0,where ReΔ=uΔ′Δ/ν and uΔ′ are the subgrid scale Reynolds number and velocity fluctuations, respectively, while δl0 is the laminar flame thickness. The efficiency function ΓΔ ([29]) estimates the net straining effect of all turbulent scales smaller than Δ. The exponent β is set constant and equal to β=0.5 as initially prescribed in [29].


## 2.3. LES Equations for Two-Phase Flow

The 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. The diluted spray is described with a Lagrangian point-force approach, which is two-way coupled to the gaseous phase. The following transport equations are solved for each droplet:(8)dxpdt=up,mpdupdt=mpup−u3CDRepρν4ρpdp2withRep=dpup−uν,where xp is the particle position, up is the particle velocity, u is the gas velocity, mp is the particle mass, CD is the drag coefficient, ν is the kinematic viscosity, ρp is the particle density, and Rep is the particle Reynolds number.The evaporation of the spray is modelled with the classical approach derived by Spalding [30]. The droplet mass transfer equation reads(9)m˙p=−πdpρDShlog1+BM,where dp is the particle diameter, D is the diffusion coefficient, Sh is the Sherwood number, and BM is the Spalding mass number. The 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:(10)ω˙˜evap=−1Vnode∑droplet∈nodem˙p,where Vnode is the volume around the node. The other droplet parameters are derived by integrating either the droplet mass or energy equations. Droplet temperature Tp and diameter dp are obtained by solving the following set of equations:(11)dTpdt=−1τpTp−T∞−LvBTCp,1/3,ddp2dt=−2Shμ1/3log1+BMdpρpSc,τp=ρpdp26ScSh⋅μ1/3Cp,kCp,1/3BTlog1+BM,where τp is the thermal characteristic time of the Spalding model, T∞ is the gas temperature in the far field, Lv is the latent heat of vaporization of the fuel, BT is the Spalding thermal number, Cp,1/3 is the heat capacity at a classical reference state assuming a one-third/two-third equilibrium between the far field and the droplet surface, μ1/3 is the dynamic viscosity at the same reference state, and Sc is the Schmidt number.


## 3. Experimental Configuration

The experimental configuration is ann-heptane spray/air jet burner experimented at CORIA by [19]. It is operated at atmospheric pressure and 298 K. The air injection is performed from a plenum to a non-swirling injector in order to generate the coflow where the liquid fuel is atomized. The air mass flow rate of is 6 g·s−1. The 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.Figure 2
Experimental setup: air path in red and spray injection in blue, from [6].Several experimental measurements have been performed. The Phase Doppler Anemometry (PDA) gives access to the gas and spray velocity and the spray diameter distribution. The 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].The flame structure shown in Figure3 by the OH-PLIF measurement exhibits a double branch. The 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. The 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.Figure 3
Instantaneous (a) and mean (b) OH-PLIF shots, from [19].
(a)
(b)


## 4. Numerical Setup

This 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]. The 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].The 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. The flame front resolution in Case B is more representative of LES conditions encountered in industrial configurations. The mesh size in the reaction zone, around 0.5 mm, is not sufficient to resolve the flame front. The 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. The modeling challenges are to recover the flame dynamic on case B, where the subgrid scale turbulent combustion model is of importance.The chemical table is built from a library of laminar freely propagatingn-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.5mm 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). The combustion model properties used for case A and B are summarized in Table 1.Table 1
Studied cases.
Case
A
B
Elements (million)
53
17
F-TACLES filter sizeΔ (mm)
0.0
3.5
Subgrid flame wrinklingΞΔ
1.0
[29]The YALES2 flow solver is used [36]. The time integration relies on a low-Mach number projection method for variable density flows. The temporal integration and spatial discretization are of fourth order. The subgrid scale Reynolds stresses are closed with the SIGMA model [37].The injected spray is polydispersed in size, following a two-parameter Rosin–Rammler distribution [38] with a Sauter Mean Diameter (SMD) d32 of 31 microns and a spread parameter q of 2.3. The form of the injected spray is 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 Figure 4.Figure 4
Particle size distribution. Experiments: grey bars, Rosin–Rammler distribution: black diamonds.


## 5. Results and Analysis

The two cases A and B are computed in both nonreactive and reactive configurations. Therefore, four simulations are presented in the following sections. The 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. The temperature field of case A-R is shown in transparency to indicate the position of the flame in reacting cases.Figure 5
Positions of the profiles of the experimental database. Red lines: radial profiles atZ = 10, 20, and 40 mm. Blue lines: axial profiles at r = 0, 5, 10, and 15 mm.


### 5.1. Flow Topology and Gas Velocity

Figures6 and 7 show for cases A-NR and B-NR, the instantaneous and mean axial velocity fields in the central vertical plane. The 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. The 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.Figure 6
Instantaneous axial velocity fields for cases A-NR (a) and B-NR (b) in the central vertical plane.
(a)
(b)Figure 7
Mean axial velocity fields for cases A-NR (a) and B-NR (b) in the central vertical plane.
(a)
(b)Figures8 and 9 show for cases A-R and B-R, the instantaneous and mean axial velocity fields in the central vertical plane. The general flow topology is similar to nonreacting cases. The 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.Figure 8
Instantaneous axial velocity fields for cases A-R (a) and B-R (b) in the central vertical plane.
(a)
(b)Figure 9
Mean axial velocity fields for cases A-R (a) and B-R (b) in the central vertical plane.
(a)
(b)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 Figures10 and 11, respectively. The 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. The main difference is the underestimation of the maximal axial velocity around 10 mm above the injection plane. The 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. The effect of the thermal expansion from the flame is visible on the profiles at Z=40mm. The axial velocity in nonreacting conditions drops rapidly to 0 m/s between r=5mm and r=20mm while in reacting conditions, the axial velocity decreases slowly between r=5mm and r=20mm. The RMS is correctly captured for both meshes. The 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=10mm.Figure 10
Radial profiles of mean axial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: Case A dashed line: Case B.
(a)
(b)Figure 11
Radial profiles of axial velocity RMS for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)The results for the radial velocity component are plotted in Figure12 (mean) and 13 (RMS). The simulations capture the general features of the mean radial velocity profiles. The RMS are also rather well predicted by the simulations. At Z=10mm, the peaks at r=0 and 10 are correctly located in the simulations. The general flow topology prediction by the simulations is satisfactory, for both meshes and for both nonreacting and reacting conditions.Figure 12
Radial profiles of mean radial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)Figure 13
Radial profiles of radial velocity RMS for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)


### 5.2. Flame Topology

Figure14(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. The inner flame front, by a highly wrinkled by the turbulence, is qualitatively reproduced by the LES. The outer diffusion flame observed in the experiments, featuring a large and unwrinkled reaction zone, is also present.Figure 14
Normalized OH mass fraction (Case A top, case B middle). Experiments (bottom): OH-PLIF, from [19]. (a) Instantaneous signal, (b) mean signal.
(a)
(b)The mean normalized OH mass fraction field for cases A and B are compared against the mean shot of OH-PLIF measurements in Figure14(a) . The mean inner flame front position is well captured by the simulations and is located at π≈15mm up to z=80mm. The mean OH-PLIF measurements show that the outer flame front extends up to x=40mm. This 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≈30mm mm and quickly merges with the inner flame front for z>50mm.The 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 Figure10. Therefore, the amount of flow-through times simulated differs between the two flame fronts. The 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).The lift-off of the flame is a critical aspect of this flame. In order to assess the lift-off height in the simulations, Figures15 and 16 show a contour of temperature in transparency for both meshes. These views demonstrate that the lift-off height is fairly constant for both meshes.Figure 15
Contour of instantaneous temperatureTinst=1300K colored by the instantaneous velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)Figure 16
Contour of mean temperatureTmean=1300K colored by the mean velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)Figures17 and 18 show a clip in the central vertical plane of the contour of temperature presented above. The 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).Figure 17
Contour of instantaneous temperatureTinst=1300K colored by the instantaneous velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)Figure 18
Contour of mean temperatureTmean=1300K colored by the mean velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)The lift-off height is defined experimentally as the closer position of the flame front from the burner exit. The flame front position is defined from the maximum value isoline given by the mean OH-PLIF signal shown in Figure14. The lift-off of the flame is estimated similarly from the simulations. This height depends on the angular position since the flame is not perfectly axisymmetric. The circumferential mean and RMS of the lift-off position are therefore computed.The experimental value is25±3mm while case A recovers a lift-off of 22±1mm and case B a lift-off of 24±1mm. 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 two-step mechanism ([6]) underpredict the flame lift-off hlo by approximately 20%. Surprisingly, simulations conducted with a reduced analytical scheme involving 24 transported species, 32 quasisteady state species, and 217 reactions also 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 meshes and for a CPU cost even lower than the global mechanism since there are only two transport equations for the chemistry (the progress variable and the mixture fraction) compared to six transported species. The good performances of F-TACLES are attributed to its ability to retrieve the flame propagation speed in turbulent stratified mixture ([14]), even on coarse grid where the flame front is not fully resolved. The edge flame propagation is, however, not influenced by the diffusion branch. The errors expected in the diffusion flame regions by the F-TACLES model which is more adapted to turbulent weakly stratified flame fronts ([14]), do not affect the lift-off height prediction in this configuration. Table 2 compares against experiments the flame lift-off height predicted by global, analytical, and tabulated chemistry on the investigated spray flame configuration. The CPU cost required to obtain reactive flow statistics, normalized by the global scheme computation, is also indicated.Table 2
Comparison between chemistry modeling strategies.
Experiment
Two-step scheme [3]
Analytical scheme [3]
F-TACLES (case A)
F-TACLES (case B)
Grid
—
Fine
Fine
Fine
Coarse
Lift-off (mm)
25
±
3
20
±
1
20
±
1
22
±
1
24
±
1
Estimated relative CPU cost
—
1
10
0.5
0.1


### 5.3. Droplets Diameter

Figures19 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. The smaller droplets are located in the central part of the flow while the larger droplets are located on the outer part of the spray. The influence of the flame in Figure 20 is the low density of particles above z=20mm, especially on the outer region.Figure 19
Droplet diameter for cases A-NR (a) and B-NR (b).
(a)
(b)Figure 20
Droplet diameter for cases A-R (a) and B-R (b).
(a)
(b)Figure21 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. The LES results show a correct evolution of the radial stratification in droplet diameter for both cases A and B. The small droplets follow the streamlines because of their small Stokes number and are therefore located at the center of the flow. The 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. The profiles are similar in both reacting and nonreacting cases between 0 and 20 mm, as flame is located further downstream. The smaller diameters encountered at 40 mm in the reacting case are the result of the stronger evaporation process due to the presence of the flame. This phenomenon is well captured by the F-TACLES model, even in case B where subgrid scale contributions are significant.Figure 21
Radial profiles of droplet diameter for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)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 Figure21 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.


### 5.4. Spray Velocity

Figures22 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.Figure 22
Droplet axial velocity for cases A-NR (a) and B-NR (b).
(a)
(b)Figure 23
Droplet axial velocity for cases A-R (a) and B-R (b).
(a)
(b)Droplet axial velocity is reported in Figure24 for the cold and reacting cases, respectively. The 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. The 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. This discrepancy is attributed to the method of injection (from [39]) that may overestimate the large droplets velocity.Figure 24
Radial profiles of droplet axial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)Droplet radial velocity is reported in Figure25 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.Figure 25
Radial profiles of droplet radial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)


### 5.5. Spray Temperature

Figures26 and 27 show the particles in the central vertical plane colored by their temperature for the cold and reacting cases, respectively. The scale is 280K<Tp<300K for the cold case and 280K<Tp<370K for the reacting cases. In Figure 26, the small droplets temperature decreases rapidly to ≈280K as they are convected downstream. This evolution is due to the evaporation. The 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 ≈370K because of the heat released by the flame. The small droplets located in the center of the flow are progressively heated by the hot gases until they are fully evaporated.Figure 26
Droplet temperature for cases A-NR (a) and B-NR (b).
(a)
(b)Figure 27
Droplet temperature for cases A-R (a) and B-R (b).
(a)
(b)The droplet temperature predicted by the LES is now compared with the Global Rainbow Technique (GRT) measurements, whose uncertainty is±3K [19]. Figure 28 presents radial profiles of temperature for the cold (left) and reacting (right) configurations.Figure 28
Radial profiles of droplet temperature for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)The experimental data highlight two zones. Forr>5mm, 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. This trend is not captured by the simulation, which predicts the wet bulb temperature for all droplet positions. The wet bulb temperature is defined as the equilibrium temperature reached by evaporating a liquid to saturation in a gas. This difference between simulations and experiments could be explained by limitations of the evaporation model ([40]).The thermal characteristic time of the Spalding model, notedτth, is expressed as(12)τth=ρpdp26ScShμ1/3Cp,kCp,1/3BTlog1+BM,where ρp is the droplet density; dp, its diameter; Sc, the Schmidt number; Sh, the Sherwood number; Cp,k, the heat capacity at constant pressure of the n-heptane; Cp,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); BT, the thermal Spalding number; and BM, 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.Figure29 presents axial profiles of temperature for the cold (left) and reacting (right) configurations. For r=0mm, the droplets (which are small at this radial position) temperature drops quickly to ≈282K. As the radial distance r increase, the mean droplet diameter growth as discussed previously, and the droplets temperature decreases. This tendency is consistent with the Spalding model assumptions.Figure 29
Axial profiles of droplet temperature for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)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 atr>10mm and z>20mm, the droplet temperature rises quickly due to the high gas temperature. This 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.


## 5.1. Flow Topology and Gas Velocity

Figures6 and 7 show for cases A-NR and B-NR, the instantaneous and mean axial velocity fields in the central vertical plane. The 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. The 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.Figure 6
Instantaneous axial velocity fields for cases A-NR (a) and B-NR (b) in the central vertical plane.
(a)
(b)Figure 7
Mean axial velocity fields for cases A-NR (a) and B-NR (b) in the central vertical plane.
(a)
(b)Figures8 and 9 show for cases A-R and B-R, the instantaneous and mean axial velocity fields in the central vertical plane. The general flow topology is similar to nonreacting cases. The 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.Figure 8
Instantaneous axial velocity fields for cases A-R (a) and B-R (b) in the central vertical plane.
(a)
(b)Figure 9
Mean axial velocity fields for cases A-R (a) and B-R (b) in the central vertical plane.
(a)
(b)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 Figures10 and 11, respectively. The 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. The main difference is the underestimation of the maximal axial velocity around 10 mm above the injection plane. The 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. The effect of the thermal expansion from the flame is visible on the profiles at Z=40mm. The axial velocity in nonreacting conditions drops rapidly to 0 m/s between r=5mm and r=20mm while in reacting conditions, the axial velocity decreases slowly between r=5mm and r=20mm. The RMS is correctly captured for both meshes. The 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=10mm.Figure 10
Radial profiles of mean axial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: Case A dashed line: Case B.
(a)
(b)Figure 11
Radial profiles of axial velocity RMS for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)The results for the radial velocity component are plotted in Figure12 (mean) and 13 (RMS). The simulations capture the general features of the mean radial velocity profiles. The RMS are also rather well predicted by the simulations. At Z=10mm, the peaks at r=0 and 10 are correctly located in the simulations. The general flow topology prediction by the simulations is satisfactory, for both meshes and for both nonreacting and reacting conditions.Figure 12
Radial profiles of mean radial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)Figure 13
Radial profiles of radial velocity RMS for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)


## 5.2. Flame Topology

Figure14(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. The inner flame front, by a highly wrinkled by the turbulence, is qualitatively reproduced by the LES. The outer diffusion flame observed in the experiments, featuring a large and unwrinkled reaction zone, is also present.Figure 14
Normalized OH mass fraction (Case A top, case B middle). Experiments (bottom): OH-PLIF, from [19]. (a) Instantaneous signal, (b) mean signal.
(a)
(b)The mean normalized OH mass fraction field for cases A and B are compared against the mean shot of OH-PLIF measurements in Figure14(a) . The mean inner flame front position is well captured by the simulations and is located at π≈15mm up to z=80mm. The mean OH-PLIF measurements show that the outer flame front extends up to x=40mm. This 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≈30mm mm and quickly merges with the inner flame front for z>50mm.The 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 Figure10. Therefore, the amount of flow-through times simulated differs between the two flame fronts. The 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).The lift-off of the flame is a critical aspect of this flame. In order to assess the lift-off height in the simulations, Figures15 and 16 show a contour of temperature in transparency for both meshes. These views demonstrate that the lift-off height is fairly constant for both meshes.Figure 15
Contour of instantaneous temperatureTinst=1300K colored by the instantaneous velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)Figure 16
Contour of mean temperatureTmean=1300K colored by the mean velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)Figures17 and 18 show a clip in the central vertical plane of the contour of temperature presented above. The 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).Figure 17
Contour of instantaneous temperatureTinst=1300K colored by the instantaneous velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)Figure 18
Contour of mean temperatureTmean=1300K colored by the mean velocity magnitude for cases A-R (a) and B-R (b).
(a)
(b)The lift-off height is defined experimentally as the closer position of the flame front from the burner exit. The flame front position is defined from the maximum value isoline given by the mean OH-PLIF signal shown in Figure14. The lift-off of the flame is estimated similarly from the simulations. This height depends on the angular position since the flame is not perfectly axisymmetric. The circumferential mean and RMS of the lift-off position are therefore computed.The experimental value is25±3mm while case A recovers a lift-off of 22±1mm and case B a lift-off of 24±1mm. 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 two-step mechanism ([6]) underpredict the flame lift-off hlo by approximately 20%. Surprisingly, simulations conducted with a reduced analytical scheme involving 24 transported species, 32 quasisteady state species, and 217 reactions also 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 meshes and for a CPU cost even lower than the global mechanism since there are only two transport equations for the chemistry (the progress variable and the mixture fraction) compared to six transported species. The good performances of F-TACLES are attributed to its ability to retrieve the flame propagation speed in turbulent stratified mixture ([14]), even on coarse grid where the flame front is not fully resolved. The edge flame propagation is, however, not influenced by the diffusion branch. The errors expected in the diffusion flame regions by the F-TACLES model which is more adapted to turbulent weakly stratified flame fronts ([14]), do not affect the lift-off height prediction in this configuration. Table 2 compares against experiments the flame lift-off height predicted by global, analytical, and tabulated chemistry on the investigated spray flame configuration. The CPU cost required to obtain reactive flow statistics, normalized by the global scheme computation, is also indicated.Table 2
Comparison between chemistry modeling strategies.
Experiment
Two-step scheme [3]
Analytical scheme [3]
F-TACLES (case A)
F-TACLES (case B)
Grid
—
Fine
Fine
Fine
Coarse
Lift-off (mm)
25
±
3
20
±
1
20
±
1
22
±
1
24
±
1
Estimated relative CPU cost
—
1
10
0.5
0.1


## 5.3. Droplets Diameter

Figures19 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. The smaller droplets are located in the central part of the flow while the larger droplets are located on the outer part of the spray. The influence of the flame in Figure 20 is the low density of particles above z=20mm, especially on the outer region.Figure 19
Droplet diameter for cases A-NR (a) and B-NR (b).
(a)
(b)Figure 20
Droplet diameter for cases A-R (a) and B-R (b).
(a)
(b)Figure21 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. The LES results show a correct evolution of the radial stratification in droplet diameter for both cases A and B. The small droplets follow the streamlines because of their small Stokes number and are therefore located at the center of the flow. The 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. The profiles are similar in both reacting and nonreacting cases between 0 and 20 mm, as flame is located further downstream. The smaller diameters encountered at 40 mm in the reacting case are the result of the stronger evaporation process due to the presence of the flame. This phenomenon is well captured by the F-TACLES model, even in case B where subgrid scale contributions are significant.Figure 21
Radial profiles of droplet diameter for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)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 Figure21 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.


## 5.4. Spray Velocity

Figures22 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.Figure 22
Droplet axial velocity for cases A-NR (a) and B-NR (b).
(a)
(b)Figure 23
Droplet axial velocity for cases A-R (a) and B-R (b).
(a)
(b)Droplet axial velocity is reported in Figure24 for the cold and reacting cases, respectively. The 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. The 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. This discrepancy is attributed to the method of injection (from [39]) that may overestimate the large droplets velocity.Figure 24
Radial profiles of droplet axial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)Droplet radial velocity is reported in Figure25 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.Figure 25
Radial profiles of droplet radial velocity for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)


## 5.5. Spray Temperature

Figures26 and 27 show the particles in the central vertical plane colored by their temperature for the cold and reacting cases, respectively. The scale is 280K<Tp<300K for the cold case and 280K<Tp<370K for the reacting cases. In Figure 26, the small droplets temperature decreases rapidly to ≈280K as they are convected downstream. This evolution is due to the evaporation. The 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 ≈370K because of the heat released by the flame. The small droplets located in the center of the flow are progressively heated by the hot gases until they are fully evaporated.Figure 26
Droplet temperature for cases A-NR (a) and B-NR (b).
(a)
(b)Figure 27
Droplet temperature for cases A-R (a) and B-R (b).
(a)
(b)The droplet temperature predicted by the LES is now compared with the Global Rainbow Technique (GRT) measurements, whose uncertainty is±3K [19]. Figure 28 presents radial profiles of temperature for the cold (left) and reacting (right) configurations.Figure 28
Radial profiles of droplet temperature for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)The experimental data highlight two zones. Forr>5mm, 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. This trend is not captured by the simulation, which predicts the wet bulb temperature for all droplet positions. The wet bulb temperature is defined as the equilibrium temperature reached by evaporating a liquid to saturation in a gas. This difference between simulations and experiments could be explained by limitations of the evaporation model ([40]).The thermal characteristic time of the Spalding model, notedτth, is expressed as(12)τth=ρpdp26ScShμ1/3Cp,kCp,1/3BTlog1+BM,where ρp is the droplet density; dp, its diameter; Sc, the Schmidt number; Sh, the Sherwood number; Cp,k, the heat capacity at constant pressure of the n-heptane; Cp,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); BT, the thermal Spalding number; and BM, 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.Figure29 presents axial profiles of temperature for the cold (left) and reacting (right) configurations. For r=0mm, the droplets (which are small at this radial position) temperature drops quickly to ≈282K. As the radial distance r increase, the mean droplet diameter growth as discussed previously, and the droplets temperature decreases. This tendency is consistent with the Spalding model assumptions.Figure 29
Axial profiles of droplet temperature for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.
(a)
(b)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 atr>10mm and z>20mm, the droplet temperature rises quickly due to the high gas temperature. This 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.


## 6. Conclusion

The first simulation with the F-TACLES formalism in a spray combustion configuration has been performed. The results show good agreement on the spray diameter and velocity, gas velocity, flame structure and lift-off with respect to experimental data. The 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. The 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. The influence of the droplet evaporation model and of the liquid sheet atomization on the spray temperature should be investigated in the future.



---

*Source: 2764523-2020-03-19.xml*